Télécharger continu.procedur

Retour à la liste

Numérotation des lignes :

  1. * CONTINU PROCEDUR BP208322 20/12/18 21:15:01 10821
  2. *=======================================================
  3. *
  4. * CONTINU
  5. * ANALyse : calcul de la REPONSE_BALOURD NONlineaire par CONTINUATION
  6. * procédure de résolution non-linéaire
  7. * par une méthode de continuation par pseudo-longueur d'arc
  8. * creation : BP 28/01/2014
  9. *
  10. *=======================================================
  11. *
  12. ************************************************************************
  13. *
  14. * CREATION : BP208322 28/01/2014
  15. * MODIF : LX236206 23/06/2015 HBM (n harmoniques), adimensionnement
  16. * analyse de stabilite
  17. * BP208322 08/2016 Mise au propre
  18. *
  19. * OBJET : Procedure de resolution Non-Lineaire
  20. * Methode de continuation par longueur d'arc
  21. *
  22. * ENTREEs et SORTIES : Voir la notice
  23. *
  24. * PROCEDURES APPELEES : CON_CALC, MONOSTA
  25. *
  26. *
  27. *
  28. ************************************************************************
  29.  
  30. DEBPROC CONTINU TAB1*'TABLE';
  31.  
  32.  
  33. mess '__________________________________________________________________________' ;
  34. mess ' .oooooo. ' ;
  35. mess ' d8P `Y8b ' ;
  36. mess '888 o8 o88 ' ;
  37. mess '888 ooooooo oo oooooo o888oo oooo oo oooooo oooo oooo ' ;
  38. mess '888 888 888 888 888 888 888 888 888 888 888 ' ;
  39. mess '`88b ooo 888 888 888 888 888 888 888 888 888 888 ' ;
  40. mess ' `Y8bood8P 88ooo88 o888o o888o 888o o888o o888o o888o 888o88 8o' ;
  41. mess '__________________________________________________________________________' ;
  42.  
  43. ************************************************************************
  44. * *
  45. * VERIFICATION DES DONNEES D'ENTREE + VALEURS PAR DEFAUT *
  46. * *
  47. ************************************************************************
  48. * fldebug = vrai;
  49. fldebug = faux;
  50.  
  51. * definition du probleme a resoudre ************************************
  52.  
  53. * prendre les matrices constantes ou HBM ?
  54. FLHBM = FAUX;
  55. SI (EXIS TAB1 'HBM');
  56. FLHBM = TAB1 . 'HBM';
  57. FINSI;
  58.  
  59. * Alternating Frequency / Time ?
  60. IAFT = 0;
  61. IDIFF = FAUX;
  62. SI (EXIS TAB1 'PROCEDURE_FREQUENCE_TEMPS');
  63. * SI (EGA TAB1 . 'PROCEDURE_FREQUENCE_TEMPS' 'DFT'); IAFT = 1; FINSI;
  64. SI (EGA TAB1 . 'PROCEDURE_FREQUENCE_TEMPS' 'AFT'); IAFT = 2; FINSI;
  65. SI (EGA IAFT 0);
  66. MESS 'PROCEDURE_FREQUENCE_TEMPS non reconnue !!!'; ERRE 21;
  67. FINSI;
  68. SI (NON (EXIS TAB1 'CALC_VITE'));
  69. TAB1 . 'CALC_VITE' = FAUX;
  70. FINSI;
  71. * si AFT, Jacobienne calculee par TFR(Knl*gamma+Cnl*w*gamma°)
  72. * ou par differentiation numerique ?
  73. SI (EGA IAFT 2);
  74. SI (EXIS TAB1 'JACOBIENNE');
  75. IDIFF = VRAI;
  76. SINON;
  77. TAB1 . 'JACOBIENNE' = mot 'TFR';
  78. * TAB1 . 'JACOBIENNE' = 'DIFFERENTIATION';
  79. FINSI;
  80. FINSI;
  81. FINSI;
  82.  
  83. * modele et materiau
  84. FLMOD1 = exis TAB1 'MODELE';
  85. SI (FLMOD1);
  86. MOD1 = TAB1 . 'MODELE';
  87. SI (EXIS TAB1 'CARACTERISTIQUES');
  88. MAT1 = TAB1 . 'CARACTERISTIQUES';
  89. SINON;
  90. MESS 'IL MANQUE LES CARACTERISTIQUES : on poursuit malgré tout...';
  91. FINSI;
  92. SINON;
  93. MESS 'IL MANQUE LE MODELE : on poursuit malgré tout... ';
  94. FINSI;
  95.  
  96. SI (EXIS TAB1 'CHARGEMENT');
  97. CHAR1 = TAB1 . 'CHARGEMENT';
  98. TYPCHAR = extr CHAR1 'COMP';
  99. SINON;
  100. SI (EGA (TYPE CHARMECA) (mot 'PROCEDUR'));
  101. TYPCHAR = mots 'MECA';
  102. SINON;
  103. * MESS 'IL MANQUE LE CHARGEMENT : on poursuit malgré tout...';
  104. MESS 'IL MANQUE LE CHARGEMENT !';
  105. ERRE 21 ;
  106. FINSI;
  107. FINSI;
  108.  
  109. * conditions initiales eventuelles *************************************
  110.  
  111. * s'agit-il d'une reprise ? (reprise si IPAS0 > 0)
  112. IPASM1 = 0;
  113. IPAS0 = 0;
  114. * SI (EXIS TAB1 'TEMPS');
  115. SI (EXIS TAB1 'TEMPS_PROG');
  116. * recherche du dernier pas enregistre
  117. IPAS0 = DIME TAB1 . 'TEMPS_PROG';
  118. repe BPAS0 (dime TAB1 . 'TEMPS_PROG'); IPAS0 = IPAS0 - 1;
  119. si (exis TAB1 . 'TEMPS' IPAS0); quit BPAS0; finsi;
  120. fin BPAS0;
  121. * recherche de l'avant-dernier pas enregistre
  122. IPASM1 = IPAS0;
  123. repe BPASM1; IPASM1 = IPASM1 - 1;
  124. si (exis TAB1 . 'TEMPS' IPASM1); quit BPASM1; finsi;
  125. fin BPASM1;
  126. TIME0 = TAB1 . 'TEMPS' . IPAS0;
  127. MESS 'Reprise du calcul au pas ' IPAS0 ' temps ' TIME0;
  128. SINON;
  129. TIME0 = 0.;
  130. TAB1 . 'TEMPS' = tabl;
  131. TAB1 . 'TEMPS' . 0 = TIME0;
  132. FINSI;
  133. TIME00=TIME0;
  134.  
  135. SI (non (EXIS TAB1 'DEPLACEMENTS'));
  136. TAB1 . 'DEPLACEMENTS' = tabl;
  137. FINSI;
  138. SI (EXIS TAB1 . 'DEPLACEMENTS' IPAS0);
  139. DEPTOT = TAB1 . 'DEPLACEMENTS' . IPAS0;
  140. SINON;
  141. si (neg IPAS0 0); MESS 'pb dans la reprise !!!'; erre 21; finsi;
  142. DEPTOT = VIDE 'CHPOINT'/'DIFFUS';
  143. TAB1 . 'DEPLACEMENTS' . IPAS0 = DEPTOT;
  144. FINSI;
  145.  
  146. SI (FLMOD1);
  147. SI (non (EXIS TAB1 'CONTRAINTES'));
  148. TAB1 . 'CONTRAINTES' = tabl;
  149. FINSI;
  150. SI (EXIS TAB1 . 'CONTRAINTES' IPAS0);
  151. SIGTOT = TAB1 . 'CONTRAINTES' . IPAS0;
  152. SINON;
  153. si (neg IPAS0 0); MESS 'pb dans la reprise !!!'; erre 21; finsi;
  154. SIGTOT = ZERO MOD1 'CONTRAIN';
  155. TAB1 . 'CONTRAINTES' . IPAS0 = SIGTOT;
  156. FINSI;
  157. FINSI;
  158.  
  159.  
  160.  
  161. * plage de pseudo-temps a calculer *************************************
  162.  
  163. SI (EXIS TAB1 'TEMPS_CALCULES');
  164. PRTIME0 = TAB1 . 'TEMPS_CALCULES';
  165. NTIME0 = DIME PRTIME0;
  166. MinTIME0 = mini PRTIME0; MaxTIME0 = maxi PRTIME0;
  167. PRDT0 = (PRTIME0 enle 1) - (PRTIME0 enle NTIME0);
  168. * croissante ou decroissante?
  169. dtimemin = mini PRDT0; dtimemax = maxi PRDT0;
  170. si ((dtimemin * dtimemax) < 0.);
  171. mess 'Les TEMPS_CALCULES doivent etre ordonnées '; erre 21;
  172. finsi;
  173. si (dtimemin < 0.);
  174. mess 'Les TEMPS_CALCULES doivent etre croissantes '; erre 21;
  175. finsi;
  176. * le 1er pas doit etre 0 si non reprise
  177. si (MinTIME0 < 0.);
  178. mess 'Le 1er TEMPS_CALCULES doit etre positif'; erre 21;
  179. finsi;
  180. si ((MinTIME0 > 0.) et (TIME0 ega 0.));
  181. mess 'Le 1er TEMPS_CALCULES est 0. (sauf si reprise)';
  182. PRTIME0 = inse PRTIME0 1 TIME0;
  183. NTIME0 = DIME PRTIME0;
  184. *bp MinTIME0 = mini PRTIME0; MaxTIME0 = maxi PRTIME0;
  185. PRDT0 = (PRTIME0 enle 1) - (PRTIME0 enle NTIME0);
  186. finsi;
  187. * * finesse (=PRDT0) dans le cas croissant
  188. * DTIME0 = (extr PRTIME0 2) - TIME0;
  189. * PRDT0 = (prog DTIME0) et ((PRTIME0 enle 1) - (PRTIME0 enle NTIME0));
  190. *bp* on prolonge pour etre coherent en longueur
  191. *bp PRDT0 = PRDT0 et (extr PRDT0 (NTIME0 - 1));
  192. * pour eviter erreurs d arrondi
  193. * lorsqu'on teste si toute la plage est couverte
  194. dtpetit = maxi (prog ((mini PRDT0 'ABS') * 1.E-3) 1.E-16);
  195. *bp : les CI ne sont pas forcement en equilibre :
  196. PRTIME0 = PRTIME0 enle 1;
  197. MinTIME0 = mini PRTIME0; MaxTIME0 = maxi PRTIME0;
  198. MinTIME0 = MinTIME0 + dtpetit; MaxTIME0 = MaxTIME0 - dtpetit;
  199. mess '* plage de temps calcules : ' MinTIME0 MaxTIME0 ' *';
  200. SINON;
  201. MESS 'IL MANQUE LE LISTREEL DES TEMPS_CALCULES'; ERRE 21;
  202. FINSI;
  203.  
  204.  
  205. * plage de frequence (cas dynamique frequentiel) ***********************
  206.  
  207. SI (EXIS TAB1 'FREQUENCE');
  208. SI (EGA (TYPE TAB1 . 'FREQUENCE') 'MOT');
  209. SI (NEG TAB1 . 'FREQUENCE' (MOT 'INCONNUE'));
  210. MESS 'FREQUENCE doit etre une EVOLUTION w(t) '
  211. 'ou le MOT INCONNUE !'; ERRE 21;
  212. FINSI;
  213. SINON;
  214. SI (NEG (TYPE TAB1 . 'FREQUENCE') 'EVOLUTIO');
  215. MESS 'FREQUENCE doit etre une EVOLUTION w(t) '
  216. 'ou le MOT INCONNUE !'; ERRE 21;
  217. FINSI;
  218. FINSI;
  219. SINON;
  220. MESS '* par defaut, on identifie : w = t';
  221. tp_tmp = ORDO (prog 0. (MinTIME0 - dtpetit) (MaxTIME0 + dtpetit));
  222. evfrqt = EVOL 'MANU' 't' tp_tmp '\w' tp_tmp;
  223. TAB1 . 'FREQUENCE' = evfrqt;
  224. FINSI;
  225.  
  226.  
  227.  
  228. * parametres de la resolution nonlineaire ******************************
  229.  
  230. * nombre d iteration maxi et ideal
  231. SI (EXIS TAB1 'MAXITERATION');
  232. NMAXIT = TAB1 . 'MAXITERATION';
  233. SINON;
  234. * NMAXIT = 50;
  235. NMAXIT = 24;
  236. TAB1 . 'MAXITERATION' = NMAXIT;
  237. FINSI;
  238. SI (EXIS TAB1 'NB_ITERATION');
  239. itideal = TAB1 . 'NB_ITERATION';
  240. si ((itideal > (NMAXIT/2)) ou (itideal <eg 4));
  241. mess 'NB_ITERATION ideal doit etre compris entre 4 et MAXITERATION/2';
  242. erre 21;
  243. finsi;
  244. SINON;
  245. itideal = 6;
  246. TAB1 . 'NB_ITERATION' = itideal;
  247. FINSI;
  248. itmin = (enti (0.5*itideal)) + 1;
  249. itmax = (enti (1.5*itideal)) + 1;
  250.  
  251. * nombre de pas total souhaité
  252. SI (EXIS TAB1 'MAXIPAS');
  253. NPAS = TAB1 . 'MAXIPAS';
  254. SINON;
  255. NPAS = 1000;
  256. TAB1 . 'MAXIPAS' = NPAS;
  257. FINSI;
  258. * nombre de sous-pas total souhaité
  259. SI (EXIS TAB1 'MAXSOUSPAS');
  260. maxsous = TAB1 . 'MAXSOUSPAS';
  261. SINON;
  262. maxsous = 5;
  263. TAB1 . 'MAXSOUSPAS' = maxsous;
  264. FINSI;
  265.  
  266. * precision
  267. SI (EXIS TAB1 'PRECISION');
  268. XPREC = TAB1 . 'PRECISION';
  269. SINON;
  270. XPREC = 1.E-6;
  271. TAB1 . 'PRECISION' = XPREC;
  272. FINSI;
  273. * norme de reference pour le residu
  274. SI (EXIS TAB1 'NORME_RESIDU');
  275. NF20 = TAB1 . 'NORME_RESIDU';
  276. SINON;
  277. NF20 = mot 'NONDEFINI';
  278. FINSI;
  279.  
  280.  
  281. * acceleration toutes les iterations
  282. SI (EXIS TAB1 'ACCELERATION');
  283. nacc = TAB1 . 'ACCELERATION';
  284. SINON;
  285. nacc = 4;
  286. TAB1 . 'ACCELERATION' = nacc;
  287. FINSI;
  288.  
  289. * LINESEARCH ?
  290. SI (EXIS TAB1 'LINESEARCH');
  291. FLAGLS = TAB1 . 'LINESEARCH';
  292. si FLAGLS; mess 'option LINESEARCH debranchee !'; erre 21; finsi;
  293. SINON;
  294. FLAGLS = FAUX;
  295. TAB1 . 'LINESEARCH' = FLAGLS;
  296. FINSI;
  297.  
  298. *Quasi Newton ou Full Newton ?
  299. SI (EXIS TAB1 'NEWTON');
  300. FLFULL0 = ega (TAB1 . 'NEWTON') 'FULL';
  301. SINON;
  302. FLFULL0 = FAUX;
  303. TAB1 . 'NEWTON' = FLFULL0;
  304. FINSI;
  305. FLFULL = FLFULL0;
  306. SI (FLFULL); MESS 'FULL NEWTON'; FINSI;
  307.  
  308. * calcul a omega fixe ?
  309. SI (EXIS TYPCHAR 'MECA');
  310. OMEFIXE = FAUX;
  311. SINON;
  312. SI (EXIS TYPCHAR 'DIMP');
  313. OMEFIXE = VRAI;
  314. SINON;
  315. MESS 'SEULEMENT LEs CHARGEMENTs MECANIQUE et DepIMP sont PREVUs !';
  316. ERRE 21;
  317. FINSI;
  318. FINSI;
  319.  
  320. * verification de la stabilité ?
  321. SI (EXIS TAB1 'STABILITE');
  322. SI (EGA (TYPE TAB1 . 'STABILITE') 'LOGIQUE');
  323. SI (TAB1 . 'STABILITE') ;
  324. SI FLHBM; MOSTAB = MOTS 'DIAG' 'FLOQ';
  325. SINON; MOSTAB = MOTS 'DIAG';
  326. FINSI;
  327. SINON; MOSTAB = MOTS ' ';
  328. FINSI;
  329. SINON;
  330. SI (EGA (TYPE TAB1 . 'STABILITE') 'LISTMOTS');
  331. MOSTAB = TAB1 . 'STABILITE';
  332. SINON;
  333. MESS 'STABILITE doit etre un LOGIQUE ou un LISTMOTS !';
  334. ERRE 21;
  335. FINSI;
  336. FINSI;
  337. si (non (exis TAB1 'RESULTATS_STABILITE'));
  338. TAB1 . 'RESULTATS_STABILITE' = TABL;
  339. finsi;
  340. SINON; MOSTAB = MOTS ' ';
  341. FINSI;
  342. *-via DIAG(K)
  343. FLDIAG = EXIS MOSTAB 'DIAG';
  344. *-via calcul de la matrice de monodromie
  345. FLMONO = EXIS MOSTAB 'MONO';
  346. *-par la theorie de HILL? (exposants de FLOQuet)
  347. FLFLOQ = (EXIS MOSTAB 'HILL') ou (EXIS MOSTAB 'FLOQ');
  348.  
  349.  
  350. * resultats a stocker *************************************************
  351.  
  352. FLRES1 = EXIS TAB1 'RESULTATS';
  353. SI (FLRES1);
  354. SI(FLHBM);
  355. TRES1 = TAB1 . 'RESULTATS_HBM';
  356. SINON;
  357. TRES1 = TAB1 . 'RESULTATS';
  358. FINSI;
  359. FINSI;
  360.  
  361. *pas a sauver
  362. SI (exis TAB1 'PAS_SAUVES');
  363. SI (ega (TYPE TAB1 . 'PAS_SAUVES') 'ENTIER');
  364. ISAUV = TAB1 . 'PAS_SAUVES';
  365. SINON;
  366. si (ega TAB1 . 'PAS_SAUVES' 'TOUS'); ISAUV = 1; finsi;
  367. si (ega TAB1 . 'PAS_SAUVES' 'AUCUN'); ISAUV = NPAS; finsi;
  368. FINSI;
  369. SINON;
  370. ISAUV = 10;
  371. FINSI;
  372. MESS 'ON SAUVE TOUS LES ' ISAUV ' PAS DE CALCULS';
  373.  
  374. *stabilite
  375. *-via DIAG(K)
  376. SI FLDIAG;
  377. SI (NON (EXIS TAB1 . 'RESULTATS_STABILITE' 'DIAG'));
  378. TAB1 . 'RESULTATS_STABILITE' . 'DIAG' = prog;
  379. TAB1 . 'RESULTATS_STABILITE' . 'DIAG_AUGMENTEE' = prog;
  380. FINSI;
  381. FINSI;
  382. *-via Hill (exposants de floquet)
  383. SI FLFLOQ;
  384. SI (NON (EXIS TAB1 . 'RESULTATS_STABILITE' 'FLOQ'));
  385. TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' = tabl;
  386. TFLOQ = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ';
  387. TFLOQ . 'EXPOSANT_REEL' = tabl;
  388. TFLOQ . 'EXPOSANT_IMAG' = tabl;
  389. nfloq=-3;
  390. TFLOQ . 'STABILITE' = lect;
  391. SINON;
  392. TFLOQ = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ';
  393. nfloq = DIME TFLOQ . 'EXPOSANT_REEL' ;
  394. FINSI;
  395. FINSI;
  396.  
  397. * resultats a afficher ************************************************
  398. SI (exis TAB1 'AFFICHAGE_RESULTATS');
  399. laffich = TAB1 . 'AFFICHAGE_RESULTATS';
  400. SINON;
  401. laffich = lect 1 2;
  402. FINSI;
  403. affich1 = extr laffich 1;
  404. affich2 = extr laffich 2;
  405. si(affich1 > affich2);
  406. mess 'les AFFICHAGE_RESULTATS doivent etre par ordre croissant!';
  407. erre 21;
  408. finsi;
  409.  
  410.  
  411.  
  412. * sorties a effectuer *************************************************
  413.  
  414. FLSOR1 = EXIS TAB1 (mot 'SORTIE');
  415. list FLSOR1;
  416. SI (FLSOR1); mess 'FLSOR1 = VRAI';
  417. SINON; mess 'FLSOR1 = faux';
  418. FINSI;
  419.  
  420. SI (FLSOR1);
  421. list TAB1;
  422. TSOR1 = TAB1 . 'SORTIE'; list TSOR1;
  423. SI (EXIS TSOR1 (mot 'FORMAT'));TYPSOR = TSOR1 . 'FORMAT';
  424. SINON; TYPSOR = mot 'VTK';
  425. FINSI;
  426. SI (EXIS TSOR1 (mot 'FICHIER'));FICSOR = TSOR1 . 'FICHIER';
  427. SINON; FICSOR = mot 'AVS/SORTIE_CON';
  428. FINSI;
  429. SI (EXIS TSOR1 (mot 'MAILLAGE'));MAISOR = TSOR1 . 'MAILLAGE' ;
  430. SINON;
  431. SI (FLMOD1); MAISOR = EXTR MOD1 'MAILLAGE';
  432. SINON; MESS 'IL MANQUE LE MAILLAGE de SORTIE'; ERRE 21;
  433. FINSI;
  434. FINSI;
  435. SI (EXIS TSOR1 (mot 'CHAMPS')); LMOSOR = TSOR1 . 'CHAMPS' ;
  436. SINON;
  437. SI (FLMOD1); LMOSOR = mots 'DEPL' 'CONT';
  438. SINON; LMOSOR = mots 'DEPL' ;
  439. FINSI;
  440. FINSI;
  441. * on suppose qu'on sort toujours le deplacement,
  442. * sort-on aussi les contraintes?
  443. FLSOR2 = EXIS LMOSOR (mot 'CONT');
  444. TCHSOR = TABL;
  445. MESS 'SORTIE AU FORMAT ' TYPSOR ' DANS LE FICHIER ' FICSOR;
  446. SINON;
  447. MESS 'PAS DE SORTIE VTK';
  448. FINSI;
  449.  
  450.  
  451. * grands deplacements ? et options associées ***************************
  452.  
  453. SI (EXIS TAB1 'GRANDS_DEPLACEMENTS');
  454. FLGDEP = TAB1 . 'GRANDS_DEPLACEMENTS';
  455. SINON;
  456. MESS 'GRANDS_DEPLACEMENTS = FAUX par defaut'; FLGDEP = FAUX;
  457. TAB1 . 'GRANDS_DEPLACEMENTS' = FLGDEP;
  458. FINSI;
  459.  
  460. SI (EXIS TAB1 'HYPOTHESE_DEFORMATIONS');
  461. HYPDEF = MOT TAB1 . 'HYPOTHESE_DEFORMATIONS';
  462. SINON;
  463. SI FLGDEP; HYPDEF = MOT 'QUADRATIQUE';
  464. SINON; HYPDEF = MOT 'LINEAIRE';
  465. FINSI;
  466. MESS 'HYPOTHESE_DEFORMATIONS = ' HYPDEF' par defaut';
  467. TAB1 . 'HYPOTHESE_DEFORMATIONS' = HYPDEF;
  468. FINSI;
  469.  
  470. SI (EXIS TAB1 'K_SIGMA');
  471. FLKSIG = TAB1 . 'K_SIGMA';
  472. SINON;
  473. MESS (chai 'K_SIGMA = ' FLGDEP ' par defaut'); FLKSIG = FLGDEP;
  474. TAB1 . 'K_SIGMA' = FLKSIG;
  475. FINSI;
  476. * recalcul de K apres la prediction?
  477. SI (EXIS TAB1 'RECALCUL_K');
  478. FLRECA = TAB1 . 'RECALCUL_K';
  479. SINON;
  480. MESS (chai 'RECALCUL_K = FAUX par defaut'); FLRECA = FAUX;
  481. TAB1 . 'RECALCUL_K' = FLRECA;
  482. FINSI;
  483.  
  484.  
  485. * procedures perso1, charmeca ? ****************************************
  486.  
  487. SI (EXIS TAB1 'PROCEDURE_PERSO1');
  488. FLPERSO1 = TAB1 . 'PROCEDURE_PERSO1';
  489. SINON;
  490. FLPERSO1 = FAUX;
  491. TAB1 . 'PROCEDURE_PERSO1' = FLPERSO1;
  492. FINSI;
  493. SI (EXIS TAB1 'PROCEDURE_CHARMECA');
  494. FLSUIV = TAB1 . 'PROCEDURE_CHARMECA';
  495. SINON;
  496. FLSUIV = FAUX;
  497. TAB1 . 'PROCEDURE_CHARMECA' = FLSUIV;
  498. FINSI;
  499.  
  500. * parametres du pilotage par longueur d'arc ****************************
  501.  
  502. *limitation de la variation de l increment relatif en deplacement
  503. *choix d'une norme = 'I' pour Identite
  504. * = 'K' pour la diagonale de la raideur tangente
  505. * = ...
  506. SI (EXIS TAB1 'NORME_DEPLACEMENT');
  507. NONORM = TAB1 . 'NORME_DEPLACEMENT';
  508. SINON;
  509. NONORM = mot 'I';
  510. FINSI;
  511. LISNORM = mots 'I' 'K';
  512. SI (non (exis LISNORM NONORM));
  513. MESS 'NORME_DEPLACEMENT est mal defini (= IDEN ou KDIA)'; erre 21;
  514. FINSI;
  515. *ordre de grandeur du deplacement attendu
  516. SI (EXIS TAB1 'MAXI_DEPLACEMENT');
  517. MAXDEPL = TAB1 . 'MAXI_DEPLACEMENT';
  518. MAXDEP2 = MAXDEPL**2;
  519. MESS 'MAXI_DEPLACEMENT =' MAXDEPL;
  520. SINO;
  521. MESS 'IL MANQUE LE MAXI_DEPLACEMENT'; ERRE 21;
  522. * a estimer + tard
  523. FINSI;
  524.  
  525. * nom des composantes pour le produit scalaire u^T * u
  526. FLMOME = faux;
  527. SI (EXIS TAB1 'COMPOSANTES');
  528. SI(FLHBM);
  529. motCf = TAB1 . 'COMPOSANTES' . 'FORCE_HBM';
  530. motCu = TAB1 . 'COMPOSANTES' . 'DEPLACEMENT_HBM';
  531. SINON;
  532. motCf = TAB1 . 'COMPOSANTES' . 'FORCE';
  533. motCu = TAB1 . 'COMPOSANTES' . 'DEPLACEMENT';
  534. si (exis TAB1 . 'COMPOSANTES' 'ROTATION');
  535. motCm = TAB1 . 'COMPOSANTES' . 'MOMENT';
  536. motCr = TAB1 . 'COMPOSANTES' . 'ROTATION';
  537. FLMOME = vrai;
  538. finsi;
  539. FINSI;
  540. SINO;
  541. MESS 'COMPOSANTES prises par defaut';
  542. si (ega (type MOD1) 'MMODEL');
  543. motCf0 = extr MOD1 'FORC';
  544. motCu0 = extr MOD1 'DEPL';
  545. * on enleve eventuel rotation et moment
  546. motCm0 = mots 'MX' 'MY' 'MZ' 'IMX' 'IMY' 'IMZ' 'MT' 'IMT';
  547. motCr0 = mots 'RX' 'RY' 'RZ' 'IRX' 'IRY' 'IRZ' 'RT' 'IRT';
  548. motCf = mots; motCu = mots; motCm = mots; motCr = mots;
  549. nCu = dime motCu0;
  550. iCu = 0;
  551. repe BCu nCu; iCu = iCu + 1;
  552. UCui = extr motCu0 iCu;
  553. si (exis motCr0 UCui); motCr = motCr et (mots UCui);
  554. sino; motCu = motCu et (mots UCui);
  555. finsi;
  556. FCui = extr motCf0 iCu;
  557. si (exis motCm0 FCui); motCm = motCm et (mots FCui);
  558. sino; motCf = motCf et (mots FCui);
  559. finsi;
  560. fin BCu;
  561. * FLMOME = (dime motCm) > 0;
  562. FLMOME = faux;
  563. sino;
  564. * motCf = mots 'FALF';
  565. * motCu = mots 'ALFA';
  566. si (neg (VALE 'MODE') 'FOUR');
  567. motfx = mots 'FX' 'FY' 'FZ';
  568. motifx = mots 'IFX' 'IFY' 'IFZ';
  569. motux = mots 'UX' 'UY' 'UZ';
  570. motiux = mots 'IUX' 'IUY' 'IUZ';
  571. motCm = mots 'MX' 'MY' 'MZ' 'IMX' 'IMY' 'IMZ';
  572. motCr = mots 'RX' 'RY' 'RZ' 'IRX' 'IRY' 'IRZ';
  573. sino;
  574. motfx = mots 'FR' 'FT' 'FZ';
  575. motifx = mots 'IFR' 'IFT' 'IFZ';
  576. motux = mots 'UR' 'UT' 'UZ';
  577. motiux = mots 'IUR' 'IUT' 'IUZ';
  578. motCm = mots 'MT' 'IMT';
  579. motCr = mots 'RT' 'IRT';
  580. fins;
  581. motCf = motfx et motifx;
  582. motCu = motux et motiux;
  583. fins;
  584. MESS 'FORCE='; list motCf;
  585. MESS 'DEPLACEMENT='; list motCu;
  586. FINSI;
  587. nCu = dime motCu;
  588. * nom de composantes definie dans le bdata
  589. LNOMDD LNOMDU = VALE 'INCO';
  590.  
  591.  
  592.  
  593. ************************************************************************
  594. * *
  595. * INITIALISATIONS *
  596. * *
  597. ************************************************************************
  598.  
  599. *creation (eventuelle) et remplissage de WTAB **************************
  600.  
  601. si (non (exis TAB1 'WTABLE')); TAB1 . 'WTABLE' = tabl; finsi;
  602. WTAB = TAB1 . 'WTABLE';
  603.  
  604. *cas reprise / initialisation de ds = longueur de pas
  605. si (exis WTAB 'DS'); dsP = WTAB . 'DS';
  606. sino; dsP = 1.; WTAB . 'DS' = dsP;
  607. finsi;
  608. si (exis WTAB 'DSMIN'); dsPmin = WTAB . 'DSMIN';
  609. sinon; dsPmin = 1.E-4;
  610. finsi;
  611. si (exis WTAB 'DSMAX'); dsPmax = WTAB . 'DSMAX';
  612. sinon; dsPmax = 1.;
  613. finsi;
  614. si (dsP < dsPmin); dsP=dsPmin; finsi;
  615. si (dsP > dsPmax); dsP=dsPmax; finsi;
  616.  
  617. *parametres lus et ranges dans WTAB (utiles pour CON_CALC , etc ...)
  618. WTAB . 'HBM' = FLHBM;
  619. WTAB . 'AFT' = IAFT;
  620. WTAB . 'DIFFERENTIATION' = IDIFF;
  621. WTAB . 'FREQ_INCONNUE' = (EGA TAB1 . 'FREQUENCE' (MOT 'INCONNUE'));
  622.  
  623. *creation point support pour la frequence si declaree comme une inconnue
  624. si (WTAB . 'FREQ_INCONNUE');
  625. si (non (exis WTAB 'PT_FREQ'));
  626. * ptf = point support de la frequence
  627. SI (EGA (VALE DIME) 2); ptf = poin 0. 0.;
  628. SINO; ptf = poin 0. 0. 0.;
  629. FINSI;
  630. WTAB . 'PT_FREQ' = ptf;
  631. finsi;
  632. * +creation du nom d'inconnue pour la frequence
  633. si (non (exis LNOMDD 'FREQ'));
  634. OPTI 'INCO' 'FREQ' 'FFRE';
  635. finsi;
  636. * +initialisation de la valeur courante de frequence si besoin
  637. si (exis (extr DEPTOT 'COMP') (mot 'FREQ'));
  638. WTAB . 'FREQ' = extr DEPTOT 'FREQ' WTAB . 'PT_FREQ';
  639. sinon;
  640. MESS 'on initialise w a 1 ...';
  641. WTAB . 'FREQ' = 1.;
  642. DEPTOT = DEPTOT
  643. et (manu 'CHPO' 1 WTAB . 'PT_FREQ' 'FREQ' 1. 'NATURE' 'DIFFUS');
  644. finsi;
  645. * +creation d'une condition de phase sur le dual de la frequence
  646. si (exis TAB1 'CHPO_PHASE');
  647. CHPPHA = TAB1 . 'CHPO_PHASE';
  648. sinon;
  649. MESS 'on utilise la condition de phase V2=0';
  650. p2 = TRES1 . 1 . 'POINT_MESURE';
  651. CHPPHA = MANU 'CHPO' 1 p2 'V2' 1.;
  652. finsi;
  653. CLPHASE = MANU 'RIGI' CHPPHA 'LIGN' 'QUEL' ptf 'FFRE';
  654. TAB1 . 'RIGIDITE_HBM' = TAB1 . 'RIGIDITE_HBM' et CLPHASE;
  655. finsi;
  656.  
  657. *preparation des mesures ***********************************************
  658. SI (FLRES1);
  659. ires1 = 0;
  660. TRES1 . ires1 = tabl;
  661. TRES1 . ires1 . 'ITERES' = prog TIME0;
  662. nRES1 = dime TRES1;
  663. repe BRES1 nRES1; ires1 = ires1 + 1;
  664. si (non (exis TRES1 ires1)); iter BRES1; finsi;
  665. moc1 = TRES1 . ires1 . 'COMPOSANTE';
  666. * -mesure de la frequence
  667. si (ega moc1 'FREQ');
  668. x1 = WTAB . 'FREQ';
  669. sinon;
  670. pt1 = TRES1 . ires1 . 'POINT_MESURE';
  671. * -mesure de force
  672. si (exis motCf moc1);
  673. * x1 = resu (redu (exco (REAC RIGTOT0 DEPTOT) moc1) pt1);
  674. * on met les reactions a 0 pour l instant
  675. x1 = 0.;
  676. sinon;
  677. * -mesure de deplacement
  678. si (exis motCu moc1);
  679. si (ega (type pt1) 'POINT');
  680. x1 = extr DEPTOT pt1 moc1 ;
  681. sinon;
  682. x1 = MAXI (REDU (EXCO DEPTOT moc1 'NOID') pt1) 'ABS';
  683. finsi;
  684. * -mesure autre (???)
  685. sinon;
  686. mess 'RESULTATS .' ires1 '. COMPOSANTE incompris !'; erre 21;
  687. finsi;
  688. finsi;
  689. finsi;
  690. * pour la reprise eventuelle
  691. si (non (exis TRES1 . ires1 'RESULTATS'));
  692. * TRES1 . ires1 . 'RESULTATS' = prog x1;
  693. *bp : les CI ne sont pas forcement en equilibre :
  694. TRES1 . ires1 . 'RESULTATS' = prog ;
  695. finsi;
  696. TRES1 . ires1 . 'ITERES' = prog x1;
  697. fin BRES1;
  698. FINSI;
  699.  
  700. *1ere sortie ***********************************************
  701. SI (FLSOR1);
  702. OPTI 'SORT' FICSOR;
  703. TCHSOR . 'DEPL' = DEPTOT;
  704. si FLSOR2;
  705. SIGGRA = CHAN 'GRAVITE' MOD1 SIGTOT;
  706. TCHSOR . 'CONT' = SIGGRA;
  707. SORT TYPSOR MAISOR TCHSOR 'TEMP' TIME0;
  708. finsi;
  709. FINSI;
  710.  
  711.  
  712. *creation (eventuelle) et remplissage de WTAB **************************
  713. * --> déplacé + haut
  714.  
  715.  
  716. *on initialise *********************************************************
  717.  
  718. *le pseudo-temps
  719. *PRTIME = liste des temps reellement calcules
  720. *DTIME0 = t_n - t_n-1
  721. si (exis TAB1 'TEMPS_PROG'); PRTIME = TAB1 . 'TEMPS_PROG';
  722. * sinon; PRTIME = prog TIME0;
  723. *bp : les CI ne sont pas forcement en equilibre :
  724. sinon; PRTIME = prog ;
  725. finsi;
  726. * DTIME0 = +1; sauf si reprise (possibilité <0)
  727. si (IPAS0 > 0);
  728. si (exis WTAB 'DTIME0');
  729. DTIME0 = WTAB . 'DTIME0';
  730. sinon;
  731. DTIME0 = (TAB1 . 'TEMPS' . IPAS0) - (TAB1 . 'TEMPS' . IPASM1);
  732. finsi;
  733. sinon;
  734. DTIME0 = ipol PRTIME0 PRDT0 TIME0;
  735. finsi;
  736. *ptt = point support du pseudo-temps
  737. SI (EGA (VALE DIME) 2); ptt = poin 0. 0.;
  738. SINON; ptt = poin 0. 0. 0.;
  739. FINSI;
  740. WTAB . 'PT_TAU' = ptt;
  741. *+creation du nom d'inconnue pour le pseudo-temps
  742. si (non (exis LNOMDD 'TAU'));
  743. OPTI 'INCO' 'TAU' 'FTAU';
  744. finsi;
  745.  
  746. *le deplacement
  747. *DDEP0 = u_n - u_n-1
  748. si (IPAS0 > 0);
  749. * reprise :
  750. si (exis WTAB 'DDEP0');
  751. DDEP0 = WTAB . 'DDEP0';
  752. sinon;
  753. DDEP0 = DEPTOT - (TAB1 . 'DEPLACEMENTS' . IPASM1);
  754. finsi;
  755. sinon;
  756. * debut du calcul :
  757. DDEP0 = DEPTOT;
  758. finsi;
  759.  
  760. *les parametres numeriques (longueur de pas ...)
  761. itconv = itideal; isconv=1;
  762. *nombre total d'iterations
  763. nittot = 0;
  764.  
  765. *recup de la config initiale (non deformee)
  766. SI (FLGDEP);
  767. GEOM0 = FORM;
  768. WTAB . 'FOR0' = GEOM0;
  769. FINSI;
  770. WTAB . 'DEPLACEMENTS' = DEPTOT;
  771. SI (FLMOD1); WTAB . 'CONTRAINTES' = SIGTOT; FINSI;
  772.  
  773.  
  774. *listmots utiles *******************************************************
  775. *pour appel a CON_CALC
  776. * mopred = mots 'RESI' 'RIGI' 'DRES';
  777. mopred = mots 'RIGI' 'DRES';
  778. mocorr = mots 'RESI' ;
  779.  
  780.  
  781. *petit entete des messages *********************************************
  782. chaline = chai '----'
  783. '--------------------------------------------------------------------';
  784. mess chaline;
  785. MESS (CHAI 'CONTINUATION PAR LONGUEUR D''''ARC ' ' ' );
  786. MESS (chai 'CEA - DEN/DM2S/SEMT/DYN - BP - 26/08/2016');
  787. mess chaline;
  788. mess ' Pas |Sous|iter| Temps |'
  789. ' Resultats | Residu ';
  790. mess (chai ' n |-Pas| (i)| t |'
  791. * ' 1 | 2 | relatif');
  792. affich1*37 ' |' affich2*50 ' | relatif');
  793. mess chaline;
  794.  
  795.  
  796.  
  797. ************************************************************************
  798. * *
  799. *==========================> Boucle BEXTERN sur les pas *
  800. * *
  801. ************************************************************************
  802. *
  803. ipas = IPAS0;
  804. REPETER BEXTERN NPAS; ipas = ipas + 1; WTAB . 'PAS' = IPAS;
  805.  
  806.  
  807. ************************************************************************
  808. * *
  809. * PREDICTION *
  810. * *
  811. ************************************************************************
  812. si(fldebug);
  813. SAUT LIGN ; mess '*** PREDICTION PAS:'ipas 't='TIME0' ***';
  814. finsi;
  815.  
  816. * longueur de pas souhaité : DT^p (a ajouter a t_n = TIME0) ***********
  817. DTP = ipol PRTIME0 PRDT0 TIME0;
  818. * on limite l'augmentation trop brutale de DTP
  819. si (ipas > (IPAS0 + 1));
  820. si (DTP > (1.10 * DTP0)); DTP = 1.10 * DTP0; finsi;
  821. finsi;
  822. DTP0 = DTP;
  823. DTP2 = DTP**2;
  824.  
  825. * Calcul des Force, Matrices et leurs dérivées par rapport a t *********
  826.  
  827.  
  828. * appel a CON_CALC qui calcule tout !
  829. * on se met dans la configuration deformee
  830. si (FLGDEP);
  831. FORM GEOM0;
  832. GEOM1 = FORM DEPTOT;
  833. finsi;
  834. CON_CALC TAB1 mopred TIME0 DTP;
  835.  
  836. * Calcul du Résidu initial
  837. * et verif qu'il est bien < precision --> inutile en general sauf pour le 1er pa
  838. si (&BEXTERN ega 1);
  839. Res2 = WTAB . 'RESIDU' ;
  840. NF2 = maxi ((PSCA Res2 Res2 motCf motCf)**0.5);
  841. mess 'norme du Residu initial =' ((XTX Res2)**0.5);
  842. finsi;
  843.  
  844. * Construction de la matrice a resoudre : RIGTOT0 = RIGTOTF + K^{NL} + ...
  845. * avec RIGTOTF = K + TIME^2*M + TIME*C +/-H
  846. * et raideur NL : K^{NL\ (0)}_n = -\frac{dF^{NL \ (0)}}{dU}
  847. * etc ...
  848. RIGTOT0 = WTAB . 'RAIDEUR' ;
  849. * Derivee de la raideur par rapport à TIME : DRIGTOT = AMOTOT1 + 2*time*M
  850. * Derivee des efforts par rapport à TIME : DRES0
  851. DRES0 = WTAB . 'DERIVEE_RESIDU';
  852.  
  853. * resolution **********************************************************
  854. DDEP2P = RESO RIGTOT0 DRES0;
  855.  
  856. * stabilité en début de pas *******************************************
  857. *-via DIAG
  858. * SI (FLDIAG);
  859. ndiag = DIAG RIGTOT0;
  860. *-via HILL (calcul des exposants de FLOQuet)
  861. SI FLFLOQ;
  862. lrprog liprog = FLOQUET TAB1;
  863. FINSI;
  864.  
  865. * calcul de alpha *****************************************************
  866. * * on peut etre tenté d'introduire le nombre de ddls, mais on ne le fait pas...
  867. * nddl = dime RIGTOT0;
  868. * nddl2 = nddl**2;
  869. * choix de la norme 2 (Identité) pour DDEP2P
  870. si (ega NONORM 'I');
  871. maxdd2 = XTX (DDEP2P enle 'LX');
  872. KDDEP2P = DDEP2P;
  873. finsi;
  874. * choix de la norme diag(K) pour DDEP2P
  875. si (ega NONORM 'K');
  876. diagK = (EXTR RIGTOT0 'DIAG') ENLE 'LX';
  877. * KDDEP2P = * diagK DDEP2P motCu motCu motCu;
  878. * maxdd2 = XTY KDDEP2P DDEP2P motCu motCu;
  879. KDDEP2P = * (diagK**0.5) DDEP2P motCu motCu motCu;
  880. maxdd2 = XTX KDDEP2P ;
  881. finsi;
  882. * calcul de alpha
  883. alpha = (maxdd2/MAXDEP2) + 1.;
  884. alpha = 1./(alpha**0.5);
  885. si OMEFIXE;
  886. salpha= +1.;
  887. sinon;
  888. salpha= ((DDEP0 DDEP2P XTY motCu motCu) / MAXDEP2)
  889. + ((DTIME0 * DTP));
  890. salpha= sign salpha;
  891. finsi;
  892.  
  893. * pas adaptatif dsP
  894. si (itconv > itmax);
  895. dsP = 0.5*dsP;
  896. si (dsP < dsPmin); dsP=dsPmin; finsi;
  897. mess 'pas adaptatif : reduction /2 -> dsP='dsP;
  898. finsi;
  899. si (itconv <eg itmin);
  900. dsP0=dsP;
  901. dsP = 1.2*dsP;
  902. si (dsP > dsPmax);
  903. dsP=dsPmax;
  904. si (dsP > dsP0);
  905. mess 'pas adaptatif : augmentation limitée -> dsP='dsP;
  906. finsi;
  907. sinon;
  908. mess 'pas adaptatif : augmentation *1.2 -> dsP='dsP;
  909. finsi;
  910. finsi;
  911. * finsi;
  912. WTAB . 'DS' = dsP;
  913. alpha0 = salpha * alpha;
  914. * si (ega ipas 3); erre 21; finsi;
  915. alpha = alpha0 * dsP;
  916. si(fldebug);
  917. mess 'dsP alpha DTP=' dsP alpha DTP;
  918. finsi;
  919.  
  920. * Mise A Jour du vecteur inconnu prenant en compte la reduction du predicteur
  921. * deplacement
  922. *dl DDEP2 = ( (1. - alpha) * (DEPTOT EXCO 'LX' 'NOID' 'LX'))
  923. *dl + (alpha * DDEP2P);
  924. DDEP2 = alpha * DDEP2P;
  925. DDEP2P = DDEP2;
  926. *dl DEP2 = (DEPTOT ENLE 'LX') + DDEP2;
  927. DEP2 = DEPTOT + DDEP2;
  928. DEP2P = DEP2;
  929. WTAB . 'DEPLACEMENTS' = DEP2;
  930. * temps
  931. DTIME0 = alpha * DTP;
  932. TIME = TIME0 + DTIME0;
  933. TIMEP = TIME;
  934. * contrainte
  935. SI (FLMOD1);
  936. DEPST = EPSI HYPDEF MOD1 MAT1 (DEP2 - DEPTOT);
  937. DSIGT = ELAS MOD1 DEPST MAT1;
  938. SIG2 = SIGTOT + DSIGT ;
  939. si (FLGDEP);
  940. SIG2 = PICA HYPDEF SIG2 (DEP2 - DEPTOT) MOD1;
  941. finsi;
  942. WTAB . 'CONTRAINTES' = SIG2;
  943. SIG2P = SIG2;
  944. FINSI;
  945.  
  946. * affichage
  947. du0 = alpha * (maxdd2**0.5);
  948. chacha0 = chai ipas*4 ' | ' 1*9 ' | ' 0*14 ' | ' TIME '|' ;
  949.  
  950.  
  951. * recup des mesures
  952. SI (FLRES1);
  953. ires1 = 0;
  954. TRES1 . ires1 . 'ITERES' = prog TIME;
  955. repe BRES1 nRES1; ires1 = ires1 + 1;
  956. si (non (exis TRES1 ires1)); iter BRES1; finsi;
  957. moc1 = TRES1 . ires1 . 'COMPOSANTE';
  958. * -mesure de la frequence
  959. si (ega moc1 'FREQ');
  960. x1 = WTAB . 'FREQ';
  961. * -autres mesures (U, F ...)
  962. sinon;
  963. pt1 = TRES1 . ires1 . 'POINT_MESURE';
  964. si (exis motCf moc1); iter BRES1; finsi;
  965. si (ega (type pt1) 'POINT');
  966. x1 = extr DEP2 pt1 moc1;
  967. sinon;
  968. x1 = MAXI (REDU (EXCO DEP2 moc1) pt1) 'ABS';
  969. finsi;
  970. finsi;
  971. TRES1 . ires1 . 'ITERES' = prog x1;
  972. * pour l'affichage
  973. * si (ires1 <eg 2);
  974. si (exis laffich ires1);
  975. chacha0 = chai chacha0 x1 '|';
  976. finsi;
  977. fin BRES1;
  978. FINSI;
  979.  
  980. si(fldebug); mess '* DTP=' DTP ' alpha =' alpha ; finsi;
  981. * erre 21;
  982.  
  983. ************************************************************************
  984. * *
  985. * CORRECTION *
  986. * *
  987. ************************************************************************
  988.  
  989.  
  990. * initialisation du pas
  991. * NREL2i= 1.E30;
  992. NF2 = NF20;
  993. DUTR1 = 1.E30;
  994. XCONV = faux;
  995. WTAB . 'RECA_K' = FLRECA ou OMEFIXE;
  996. RECA_CO= NON OMEFIXE;
  997. isouspas = 0;
  998. * recup d'info de la prediction (considérée comme l'iteration 0)
  999. DTIME = DTIME0;
  1000. DRES0 = alpha*DRES0;
  1001. Res2i = DRES0;
  1002. DEPTOTi=DEPTOT;
  1003. SI (FLMOD1); SIGTOTi=SIGTOT; FINSI;
  1004.  
  1005. *-----------------------------------------------------> boucle iterative
  1006. REPE BSOUPAS maxsous; isouspas = isouspas + 1;
  1007.  
  1008. si(fldebug); SAUT LIGN ; mess '*** CORRECTION ***' isouspas; finsi;
  1009.  
  1010.  
  1011. * initialisation du pas
  1012. it = 0;
  1013. iacc = 0; itacc = nacc;
  1014. * variables du linesearch : compteur, longueur beta , activation du LS
  1015. nls = 10; ils = 0; beta = 1.; FLLS = FAUX;
  1016.  
  1017. *-----------------------------------------------------> boucle iterative
  1018. REPE BRESONL (nls*NMAXIT); it = it + 1;
  1019.  
  1020. SI (it > NMAXIT) ; QUIT BRESONL; FINSI;
  1021.  
  1022.  
  1023. *---Mise a jour :
  1024. * -de la matrice de raideur :
  1025. SI (WTAB . 'RECA_K' ou FLFULL);
  1026. * appel a CON_CALC qui calcule tout !
  1027. * on se met dans la derniere configuration deformee
  1028. si (FLGDEP);
  1029. FORM GEOM0;
  1030. GEOM2 = FORM DEP2;
  1031. finsi;
  1032. CON_CALC TAB1 mopred TIME DTIME;
  1033. * la matrice a resoudre RIGTOT0 = K + OMEG^2*M + OMEG*C +/-H
  1034. mess 'recalcul de K';
  1035. RIGTOT0 = WTAB . 'RAIDEUR' ;
  1036. DRES0 = WTAB . 'DERIVEE_RESIDU';
  1037. FINSI;
  1038.  
  1039. *---Mise a jour :
  1040. * -des matrices colonne et ligne supplementaires :
  1041. SI (RECA_CO ou FLFULL);
  1042. * on adimensionne DRES0
  1043. dRdT = MANU 'RIGI' ((-1./DTIME)*DRES0) 'COLO' 'QUEL' ptt 'TAU';
  1044. * on fait les corrections dans le plan orthogonal a la prediction
  1045. si (ega NONORM 'I');
  1046. DUPT = MANU 'RIGI' ((DDEP2P) enle 'LX') 'LIGN' 'QUEL' ptt 'FTAU';
  1047. finsi;
  1048. * choix de la norme diag(K) pour DDEP2P
  1049. si (ega NONORM 'K');
  1050. diagK = (EXTR RIGTOT0 'DIAG') ENLE 'LX';
  1051. KDDEP2 = * diagK (DDEP2P) motCu motCu motCu;
  1052. DUPT = MANU 'RIGI' (KDDEP2 enle 'LX') 'LIGN' 'QUEL' ptt 'FTAU';
  1053. finsi;
  1054. DTT = MANU 'RIGI' ptt (mots 'TAU')
  1055. 'DUAL' (mots 'FTAU') (prog DTIME);
  1056. DUPT = (1./MAXDEP2) * DUPT;
  1057. DTT = (1./DTP2) * DTT ;
  1058. * rem: on pourrait sans doute mieux adimensionner ca, voire meme
  1059. * utiliser une expression symetrique de dRdT a la place de DUPT
  1060. FINSI;
  1061. * -de l'operateur final :
  1062. SI ((NON OMEFIXE) et (WTAB . 'RECA_K' ou RECA_CO ou FLFULL));
  1063. RIGTOT2 = RIGTOT0 et dRdT et DUPT et DTT;
  1064. FINSI;
  1065. WTAB . 'RECA_K' = faux;
  1066. RECA_CO= faux;
  1067.  
  1068.  
  1069. *---appel a CON_CALC qui calcule le residu dans la config n+1
  1070. SI (FLGDEP);
  1071. FORM GEOM0;
  1072. GEOM2 = FORM DEP2;
  1073. FINSI;
  1074. CON_CALC TAB1 mocorr TIME DTIME;
  1075. * attention : si on détecte un changement brutal de comportement (mise
  1076. * en contact, etc ..), alors il faut recommencer en recalculant K
  1077. si (WTAB . 'RECA_K');
  1078. it = it - 1;
  1079. ITER BRESONL;
  1080. finsi;
  1081. * on recupere le residu fraichement calculé
  1082. Res2 = WTAB . 'RESIDU' ;
  1083. * puis retour a la config n (en debut de pas)
  1084. SI (FLGDEP);
  1085. FORM GEOM1 ;
  1086. FINSI;
  1087.  
  1088. * FINSI;
  1089.  
  1090. *---calcul des reactions
  1091. *dl Freac2 = REAC KBLO1 DEP2;
  1092. Freac2 = WTAB . 'REACTION';
  1093.  
  1094. *---calcul de la norme NF2 de ce pas
  1095. si (ega NF2 'NONDEFINI');
  1096. NF2 = maxi ((PSCA Res2 Res2 motCf motCf)**0.5);
  1097. NF2r = maxi ((PSCA Freac2 Freac2 motCf motCf)**0.5);
  1098. NF2 = NF2 + NF2r;
  1099. si (NF2 < 1.E-30); mess 'def de NF2 a revoir'; erre 21; finsi;
  1100. si (FLMOME);
  1101. * NF2M = maxi ((PSCA Res2 Res2 motCm motCm)**0.5);
  1102. NF2M = ((XTX Res2)**0.5)
  1103. * ((maxi DDEP2P 'ABS' 'AVEC' motCu)
  1104. /(maxi DDEP2P 'ABS' 'AVEC' motCr));
  1105. si (NF2M < 1.E-30); mess 'def de NF2M a revoir'; erre 21; finsi;
  1106. * mess 'normes pour tester la convergence du Residu =' NF2 NF2M;
  1107. sino;
  1108. * mess 'norme pour tester la convergence du Residu =' NF2;
  1109. finsi;
  1110. finsi;
  1111.  
  1112. * -condition d'equilibre des forces :
  1113. *dl* Equi^{(i)} = Res^{(i)} - C^T*\lambda^{(i)}
  1114. *dl Equi2 = (Res2 + (REAC KBLO1 DEP2)) ENLE 'FLX';
  1115. Equi2 = Res2 ENLE 'FLX';
  1116. Neq2 = maxi ((PSCA Equi2 Equi2 motCf motCf)**0.5);
  1117. NREL2 = Neq2 / NF2;
  1118. chacha = chai chacha0 FORMAT '(E13.6)' NREL2*69;
  1119. si (FLMOME);
  1120. Neq2M = maxi ((PSCA Equi2 Equi2 motCm motCm)**0.5);
  1121. NREL2M = Neq2M / NF2M;
  1122. chacha = chai chacha '|' FORMAT '(E13.6)' NREL2M;
  1123. finsi;
  1124.  
  1125. *---affichage
  1126. mess chacha;
  1127.  
  1128. *---Convergence?
  1129. XCONV = (NREL2 < XPREC);
  1130. si (FLMOME); XCONV = XCONV et (NREL2M < XPREC); finsi;
  1131. SI XCONV;
  1132. isconv = isouspas;
  1133. itconv = it;
  1134. SI (ils > 0);
  1135. si(fldebug); mess 'LineSearch: du^T*R(' beta ')=' DUTR2; finsi;
  1136. FINSI;
  1137. QUIT BRESONL;
  1138. FINSI;
  1139.  
  1140. * SI (NON FLFULL);
  1141. *---Acceleration d'Aitken pour le Residu
  1142. * avec CHR0k = Residu^(i-k) et CHEOk = Equilibre^(i-k)
  1143. si (ega itacc 4); CHR03 = Res2; CHE03 = Equi2; finsi;
  1144. si (ega itacc 3); CHR02 = Res2; CHE02 = Equi2; finsi;
  1145. si (ega itacc 2); CHR01 = Res2; CHE01 = Equi2; finsi;
  1146. si (ega itacc 1); CHR00 = Res2; CHE00 = Equi2; finsi;
  1147. itacc = itacc - 1;
  1148. si (ega itacc 0);
  1149. CORREC = ACT3 CHE02 CHE01 CHE00 CHE03 CHE02 CHE01 CHE00;
  1150. Res2 = Res2 - CORREC;
  1151. itacc = nacc;
  1152. mess 'acceleration Residu';
  1153. finsi;
  1154. *---> ne semble pas tres compatible avec le linesearch...
  1155. * FINSI;
  1156.  
  1157.  
  1158. *---Resolution du pb augmenté :
  1159. si(fldebug); saut lign; mess '*** iteration : ' it; finsi;
  1160. * SI (it ega 2) ; erre 21; FINSI;
  1161. * [K^{tan (i)} C^T dR/dt ] (\Delta U ) (Res^{(i-1)})
  1162. * [C 0 0 ] * (\lambda^{(i)}) = (U^{imp}_n )
  1163. * [Du(i-1)^T 0 Dt(i-1)] (\Delta t ) (0 )
  1164.  
  1165. * si (ega ipas 21); erre 21; finsi;
  1166.  
  1167. si (OMEFIXE);
  1168. DDEP2 = RESO RIGTOT0 Res2;
  1169. DTIME = 0.;
  1170. sinon;
  1171. DDEP2 = RESO RIGTOT2 Res2;
  1172. * -Extraction de DT et de (DU,lambda)
  1173. DTIME = EXTR DDEP2 'TAU' ptt ;
  1174. DDEP2 = DDEP2 ENLE 'TAU';
  1175. * rem : si thermique faudra changer...
  1176. finsi;
  1177.  
  1178.  
  1179. *---Actualisation du vecteur inconnu :
  1180.  
  1181. * -calcul du alpha^(i) de l itéré
  1182. * la longueur de l'itéré doit etre << dsP
  1183. si (ega NONORM 'I');
  1184. maxdd2i = XTX (DDEP2 enle 'LX');
  1185. finsi;
  1186. si (ega NONORM 'K');
  1187. * KDDEP2 = * (diagK) DDEP2 motCu motCu motCu;
  1188. * maxdd2 = XTY KDDEP2 DDEP2 motCu motCu;
  1189. KDDEP2 = * (diagK**0.5) DDEP2 motCu motCu motCu;
  1190. maxdd2i = XTX KDDEP2 ;
  1191. finsi;
  1192. alphai = (maxdd2i/MAXDEP2) + ((DTIME/DTP)**2);
  1193. alphai = dsP / (alphai**0.5);
  1194. si(fldebug); mess 'DTIME= ' DTIME ' alphai= ' alphai ; finsi;
  1195. si (alphai < 1.);
  1196. si (alphai < 1.E-6); alphai = 1.E-6; finsi;
  1197. mess 'reduction de l itéré par ' alphai;
  1198. DTIME = alphai * DTIME;
  1199. *dl DDEP2 = ((1.-alphai) * (exco DEP2 'LX' 'NOID' 'LX'))
  1200. *dl + (alphai * DDEP2);
  1201. DDEP2 = alphai * DDEP2;
  1202. sinon;
  1203. alphai = 1.;
  1204. finsi;
  1205.  
  1206. * -on sauve les vecteurs a literation i + Mise A Jour du vecteur inconnu
  1207. * U^{(i)} = U^{(i-1)} + \Delta U
  1208. * \lambda^{(i)} = 0 + \Delta \lambda
  1209. * deplacement
  1210. DEP2i = DEP2;
  1211. *dl DEP2 = (ENLE DEP2i 'LX') + DDEP2;
  1212. DEP2 = DEP2i + DDEP2;
  1213. WTAB . 'DEPLACEMENTS' = DEP2;
  1214. * temps
  1215. TIMEi = TIME;
  1216. DTIMEi= DTIME;
  1217. TIME = TIMEi + DTIMEi;
  1218. * residu
  1219. * Res2i = Res2;
  1220.  
  1221. * -increment de deformation et de contrainte calculé sur config n
  1222. * rem : pourrait etre calculée sur config intermediaire t_n+0.5
  1223. SI (FLMOD1);
  1224. DEPST = EPSI HYPDEF MOD1 MAT1 (DEP2 - DEPTOTi);
  1225. DSIGT = ELAS MOD1 DEPST MAT1;
  1226. SIG2i = SIG2;
  1227. SIG2 = SIGTOTi + DSIGT ;
  1228. * SIG2 = contrainte PK2 calculée sur config n,
  1229. * car SIGTOT = contrainte de Cauchy au temps n
  1230. * = contrainte de PK au temps n sur config n
  1231. * (cf cours Brunet, p.43)
  1232. si (FLGDEP);
  1233. * on calcul les contraintes de Cauchy (sur config n+1)
  1234. SIG2 = PICA HYPDEF SIG2 (DEP2 - DEPTOTi) MOD1;
  1235. finsi;
  1236. WTAB . 'CONTRAINTES' = SIG2;
  1237. FINSI;
  1238.  
  1239. *---affichage
  1240. * dup = alphai * (maxdd2i**0.5);
  1241. chacha0 = chai ipas*4 ' | ' isouspas*9 ' | ' it*14 ' | ' TIME '|';
  1242.  
  1243. *---recup des mesures
  1244. SI (FLRES1);
  1245. ires1 = 0;
  1246. TRES1 . ires1 . 'ITERES' = TRES1 . ires1 . 'ITERES' et TIME;
  1247. repe BRES1 nRES1; ires1 = ires1 + 1;
  1248. si (non (exis TRES1 ires1)); iter BRES1; finsi;
  1249. moc1 = TRES1 . ires1 . 'COMPOSANTE';
  1250. * -mesure de la frequence
  1251. si (ega moc1 'FREQ');
  1252. x1 = WTAB . 'FREQ';
  1253. * -autres mesures (U, F ...)
  1254. sinon;
  1255. pt1 = TRES1 . ires1 . 'POINT_MESURE';
  1256. si (exis motCf moc1); iter BRES1; finsi;
  1257. si (ega (type pt1) 'POINT');
  1258. x1 = extr DEP2 pt1 moc1;
  1259. sinon;
  1260. x1 = MAXI (REDU (EXCO DEP2 moc1) pt1) 'ABS';
  1261. finsi;
  1262. finsi;
  1263. TRES1 . ires1 . 'ITERES' = TRES1 . ires1 . 'ITERES' et x1;
  1264. * pour l'affichage
  1265. * si (ires1 <eg 2);
  1266. si (exis laffich ires1);
  1267. chacha0 = chai chacha0 x1 '|';
  1268. finsi;
  1269. fin BRES1;
  1270. FINSI;
  1271.  
  1272.  
  1273. FIN BRESONL;
  1274. *<------------------------------------------- fin de la boucle iterative
  1275.  
  1276.  
  1277. ************************************************************************
  1278. * *
  1279. * NON CONVERGENCE *
  1280. * *
  1281. ************************************************************************
  1282. *-on sort si convergence!
  1283. si (XCONV);
  1284. nittot = nittot + itconv;
  1285. quit BSOUPAS;
  1286.  
  1287. *-erreur si non-convergence!
  1288. sino;
  1289. mess 'non convergence en ' NMAXIT 'iterations !';
  1290.  
  1291. si (isouspas >eg maxsous);
  1292. mess 'non convergence en ' isouspas 'sous-pas!';
  1293. * erre 997;
  1294. mess 'On retourne le resultat malgré tout...';
  1295. * si le dernier pas converge n'a pas ete sauve,
  1296. * on le sauve pour eventuelle reprise
  1297. SI (NON (MULT (ipas - 1) ISAUV));
  1298. * sauvegarde du dernier pas converge
  1299. TAB1 . 'TEMPS' . (ipas - 1) = TIME0;
  1300. TAB1 . 'DEPLACEMENTS' . (ipas - 1) = DEPTOTi;
  1301. SI (FLMOD1); TAB1 . 'CONTRAINTES' . ipas = SIGTOTi; FINSI;
  1302. FINSI;
  1303. * doit on faire qqchose sur WTAB ?...?
  1304. quit BEXTERN;
  1305. finsi;
  1306.  
  1307. * si (isouspas < maxsous);
  1308. * si ((isouspas <eg 3) et (dsP > dsPmin));
  1309. si ((isouspas <eg 3));
  1310.  
  1311. si ((dsP > dsPmin));
  1312. * on recommence avec un pas predicitif 1/xred fois + petit
  1313. xred = 0.2;
  1314. dsP = xred * dsP; si (dsP < dsPmin); dsP=dsPmin; finsi;
  1315. alpha = xred * alpha;
  1316. mess ' reduction du pas predicteur de ' xred;
  1317. sinon;
  1318. * on recommence avec un pas predicitif tq ds = 1 !
  1319. alpha = alpha0;
  1320. xred = 1. / dsP;
  1321. dsP = 1;
  1322. mess ' augmentation du pas predicteur de ' xred;
  1323. finsi;
  1324. * deplacement
  1325. *dl DDEP2 = ( (1. - alpha) * (DEPTOT EXCO 'LX' 'NOID' 'LX'))
  1326. *dl + (alpha * DDEP2P);
  1327. DDEP2 = (alpha * DDEP2P);
  1328. DDEP2P = DDEP2;
  1329. *dl DEP2 = (DEPTOT ENLE 'LX') + DDEP2;
  1330. DEP2 = DEPTOT + DDEP2;
  1331. DEP2P = DEP2;
  1332. WTAB . 'DEPLACEMENTS' = DEP2;
  1333. * temps
  1334. DTIME = alpha * DTP;
  1335. DTP = xred * DTP;
  1336. TIME = TIME0 + DTIME;
  1337. TIMEP = TIME;
  1338. * contrainte
  1339. SI (FLMOD1);
  1340. DEPST = EPSI HYPDEF MOD1 MAT1 (DEP2 - DEPTOT);
  1341. DSIGT = ELAS MOD1 DEPST MAT1;
  1342. SIG2 = SIGTOT + DSIGT ;
  1343. si (FLGDEP);
  1344. SIG2 = PICA HYPDEF SIG2 (DEP2 - DEPTOT) MOD1;
  1345. finsi;
  1346. WTAB . 'CONTRAINTES' = SIG2;
  1347. SIG2P = SIG2;
  1348. FINSI;
  1349. * affichage
  1350. * du0 = alpha * (maxdd2**0.5);
  1351. chacha0 = chai ipas*4 ' | ' (isouspas+1)*9 ' | ' it*14 ' | '
  1352. TIME '|';
  1353. *---recup des mesures
  1354. SI (FLRES1);
  1355. ires1 = 0;
  1356. TRES1 . ires1 . 'ITERES' = TRES1 . ires1 . 'ITERES' et TIME;
  1357. repe BRES1 nRES1; ires1 = ires1 + 1;
  1358. si (non (exis TRES1 ires1)); iter BRES1; finsi;
  1359. moc1 = TRES1 . ires1 . 'COMPOSANTE';
  1360. * -mesure de la frequence
  1361. si (ega moc1 'FREQ');
  1362. x1 = WTAB . 'FREQ';
  1363. * -autres mesures (U, F ...)
  1364. sinon;
  1365. pt1 = TRES1 . ires1 . 'POINT_MESURE';
  1366. si (exis motCf moc1); iter BRES1; finsi;
  1367. si (ega (type pt1) 'POINT');
  1368. x1 = extr DEP2 pt1 moc1;
  1369. sinon;
  1370. x1 = MAXI (REDU (EXCO DEP2 moc1) pt1) 'ABS';
  1371. finsi;
  1372. finsi;
  1373. TRES1 . ires1 . 'ITERES' = TRES1 . ires1 . 'ITERES' et x1;
  1374. * pour l'affichage
  1375. * si (ires1 <eg 2);
  1376. si (exis laffich ires1);
  1377. chacha0 = chai chacha0 x1 '|';
  1378. finsi;
  1379. fin BRES1;
  1380. FINSI;
  1381.  
  1382. sinon;
  1383. * autre chance :
  1384. RECA_K = vrai;
  1385. RECA_CO = NON OMEFIXE;
  1386. * Pour le prochain sous-pas, on repart de la position non-convergée
  1387. mess ' on redemarre depuis la derniere position (non-convergée)';
  1388. DTP = TIME - TIME0;
  1389. * TIME0 = TIME;
  1390. DDEP2P = DEP2 - DEPTOT;
  1391. DEPTOTi = DEP2;
  1392. WTAB . 'DEPLACEMENTS' = DEPTOTi;
  1393. SI (FLMOD1);
  1394. SIGTOTi = SIG2;
  1395. WTAB . 'CONTRAINTES' = SIGTOTi;
  1396. FINSI;
  1397. SI (FLGDEP);
  1398. FORM GEOM0;
  1399. GEOM1 = FORM DEPTOTi;
  1400. FINSI;
  1401.  
  1402. finsi;
  1403.  
  1404. fins;
  1405.  
  1406. FIN BSOUPAS;
  1407. *<------------------------------------------- fin de la boucle iterative
  1408.  
  1409.  
  1410. ************************************************************************
  1411. * *
  1412. * CONVERGENCE *
  1413. * *
  1414. ************************************************************************
  1415.  
  1416. *-Ici on a convergé : ******************************
  1417.  
  1418. * Pour le prochain pas
  1419. DTIME0 = TIME - TIME0;
  1420. TIME0 = TIME;
  1421. DDEP0 = DEP2 - DEPTOT;
  1422. DEPTOT = DEP2;
  1423. SI (FLMOD1); SIGTOT = SIG2; FINSI;
  1424.  
  1425. * affichage
  1426. maxdd2 = XTX (DDEP0 enle 'LX');
  1427. ds2 = ((maxdd2 / MAXDEP2) + ((DTIME0 / DTP)**2))**0.5;
  1428. mess 'convergence en ' itconv ' iterations ' ndiag;
  1429. dsPP = dsP;
  1430. si (ds2 < (0.99*dsP)); mods = mot ' < ';
  1431. * si on est trop grand c'est que les corrections sont grandes par
  1432. * rapport a la prediction -> "on tourne" vite -> il faut reduire le pas
  1433. sinon; si (ds2 > (1.2*dsP)); mods = mot ' > ';
  1434. * dsP = 0.5*dsP;
  1435. sinon; mods = mot ' ~ ';
  1436. finsi;
  1437. finsi;
  1438. chaconv = chai 'longueur du pas Ds=' ds2 mods dsPP
  1439. ' (Du=' (maxdd2**0.5) ' Dt=' DTIME0 ')' ;
  1440. mess chaconv;
  1441.  
  1442. * si convergence on enregistre cette frequence
  1443. PRTIME = PRTIME et TIME;
  1444.  
  1445. * sauvegarde pour eventuelle reprise
  1446. WTAB . 'DEPLACEMENTS' = DEPTOT;
  1447. SI (FLMOD1); WTAB . 'CONTRAINTES' = SIGTOT; FINSI;
  1448. WTAB . 'DDEP0' = DDEP0;
  1449. WTAB . 'DTIME0' = DTIME0;
  1450. WTAB . 'TIME0' = TIME0;
  1451.  
  1452. * sauvegarde du pas
  1453. SI (MULT ipas ISAUV);
  1454. TAB1 . 'TEMPS' . ipas = TIME;
  1455. TAB1 . 'DEPLACEMENTS' . ipas = DEPTOT;
  1456. SI (FLMOD1); TAB1 . 'CONTRAINTES' . ipas = SIGTOT; FINSI;
  1457. FINSI;
  1458.  
  1459. * eventuellement recombinaison modale
  1460.  
  1461.  
  1462. * Sauvegarde des Resultats ******************************
  1463.  
  1464. * recup des mesures
  1465. SI (FLRES1);
  1466. nRES1 = dime TRES1;
  1467. ires1 = 0;
  1468. repe BRES1 nRES1; ires1 = ires1 + 1;
  1469. si (non (exis TRES1 ires1)); iter BRES1; finsi;
  1470. moc1 = TRES1 . ires1 . 'COMPOSANTE';
  1471. * -mesure de la frequence
  1472. si (ega moc1 'FREQ');
  1473. x1 = WTAB . 'FREQ';
  1474. sinon;
  1475. * -mesure de U ou F
  1476. pt1 = TRES1 . ires1 . 'POINT_MESURE';
  1477. si (exis motCf moc1);
  1478. x1 = maxi (resu (redu (exco (REAC RIGTOT0 DEPTOT) moc1) pt1));
  1479. sinon;
  1480. si (ega (type pt1) 'POINT');
  1481. x1 = extr DEPTOT pt1 moc1;
  1482. sinon;
  1483. x1 = MAXI (REDU (EXCO DEPTOT moc1) pt1) 'ABS';
  1484. finsi;
  1485. finsi;
  1486. finsi;
  1487. TRES1 . ires1 . 'RESULTATS' = TRES1 . ires1 . 'RESULTATS' et x1;
  1488. fin BRES1;
  1489. FINSI;
  1490.  
  1491. * stockage de la stabilite (calculee avec les matrices en debut de pas)
  1492. *-via DIAG
  1493. SI FLDIAG;
  1494. TAB1 . 'RESULTATS_STABILITE' . 'DIAG' =
  1495. TAB1 . 'RESULTATS_STABILITE' . 'DIAG' et (flot ndiag);
  1496. * En theorie, il faudrait aussi trouver un moyen pour detecter les
  1497. * changements de signe de la jacobienne augmentee
  1498. * et s'assurer que la matrice utilisee soit bien celle en début de pas
  1499. * Ici, on se contente de prendre la valeur ci dessous
  1500. * (dont le sens reste a prouver car matrice non symetrique) :
  1501. ndiag2 = DIAG RIGTOT2;
  1502. TAB1 . 'RESULTATS_STABILITE' . 'DIAG_AUGMENTEE' =
  1503. TAB1 . 'RESULTATS_STABILITE' . 'DIAG_AUGMENTEE' et (flot ndiag2);
  1504. FINSI;
  1505. *-via HILL (calcul des exposants de FLOQuet)
  1506. SI FLFLOQ;
  1507. * lrprog liprog = FLOQUET TAB1; --> fait en debut de pas
  1508. * initialisation des objets en sortie
  1509. si (nfloq ega -3);
  1510. nfloq = dime lrprog;
  1511. ifloq = 0;
  1512. repe Bfloq nfloq; ifloq = ifloq + 1;
  1513. TFLOQ . 'EXPOSANT_REEL' . ifloq = prog ;
  1514. TFLOQ . 'EXPOSANT_IMAG' . ifloq = prog ;
  1515. fin Bfloq ;
  1516. * finsi;
  1517. *bp : les CI ne sont pas forcement en equilibre :
  1518. sinon;
  1519. * remplissage des tables de listreel en sortie
  1520. ifloq = 0;
  1521. repe Bfloq nfloq; ifloq = ifloq + 1;
  1522. TFLOQ . 'EXPOSANT_REEL' . ifloq
  1523. = TFLOQ . 'EXPOSANT_REEL' . ifloq et (extr lrprog ifloq);
  1524. TFLOQ . 'EXPOSANT_IMAG' . ifloq
  1525. = TFLOQ . 'EXPOSANT_IMAG' . ifloq et (extr liprog ifloq);
  1526. fin Bfloq ;
  1527. * on teste la stabilite sur tous les exposants de Floquet
  1528. * OU
  1529. * on ne conserve que les (2*nddls) valeurs centrées autour de 0
  1530. lrprog2 liprog2 = ORDO 'DECROISSANT' lrprog liprog;
  1531. toto liprog lrprog = ORDO 'CROISSANT' (ABS liprog2) liprog2 lrprog2;
  1532. NDDL = ((dime lrprog) / 2) / (2*(TAB1 . 'N_HARMONIQUE') + 1);
  1533. lprem = lect 1 pas 1 (2 * NDDL);
  1534. lrprog = EXTR lrprog lprem;
  1535. minlr = maxi lrprog;
  1536. si (minlr > 1.E-16); ifloq = 1; sinon; ifloq = 0; finsi;
  1537. TFLOQ . 'STABILITE' = TFLOQ . 'STABILITE' et ifloq;
  1538. finsi;
  1539. FINSI;
  1540.  
  1541.  
  1542. * sortie ***********************************************
  1543.  
  1544. SI (FLSOR1);
  1545. * sortie sur la config initiale
  1546. SI (FLGDEP); FORM GEOM0; finsi;
  1547. OPTI 'SORT' FICSOR;
  1548. TCHSOR . 'DEPL' = DEPTOT;
  1549. si FLSOR2;
  1550. SIGGRA = CHAN 'GRAVITE' MOD1 SIGTOT;
  1551. TCHSOR . 'CONT' = SIGGRA;
  1552. * SORT TYPSOR MAISOR DEPTOT 'DEPL' SIGGRA 'SIGM' 'TEMP' TIME;
  1553. * sinon;
  1554. * SORT TYPSOR MAISOR DEPTOT 'DEPL' 'TEMP' TIME;
  1555. finsi;
  1556. SORT TYPSOR MAISOR TCHSOR 'TEMP' TIME0;
  1557. FINSI;
  1558.  
  1559. * trait de fin du pas ******************************************
  1560. mess chaline;
  1561.  
  1562. * On sort si on a fini ! **************************************
  1563. TAB1 . 'TEMPS_PROG' = PRTIME;
  1564. PRTIME1 = PRTIME et TIME00;
  1565. SI ( ((maxi PRTIME1) >eg MaxTIME0) et ((mini PRTIME1) <eg MinTIME0));
  1566. MESS ( CHAI 'Plage temporelle couverte = '
  1567. (mini PRTIME) (maxi PRTIME) );
  1568. mess chaline;
  1569. MESS '----------------------> Succès de CONTINU'
  1570. ' <----------------------';
  1571. QUIT BEXTERN;
  1572. FINSI;
  1573.  
  1574.  
  1575. FIN BEXTERN;
  1576. ************************************************************************
  1577. * *
  1578. *<=== Fin de la Boucle BEXTERN sur les pas de chargement
  1579. * *
  1580. ************************************************************************
  1581.  
  1582.  
  1583. mess (chai 'nombre total d' 'iteration ' ': ' nittot);
  1584. mess chaline;
  1585.  
  1586.  
  1587. ************************************************************************
  1588. * *
  1589. * SAUVEGARDE DES RESULTATS (convergés jusqu'à arret de bextern) *
  1590. * *
  1591. ************************************************************************
  1592.  
  1593. * evolutions
  1594. SI (FLRES1);
  1595. * tables des titres pour l operateur DESSIN
  1596. TRES1 . 'TITRE' = tabl;
  1597. mycoul1 = mots 'VIOL' 'BLEU' 'TURQ' 'VERT' 'OLIV'
  1598. 'ORAN' 'ROUG' 'ROSE' 'AZUR' 'DEFA';
  1599. ncoul1 = dime mycoul1; icoul1 = 0;
  1600. evtot = vide 'EVOLUTIO' ;
  1601. nRES1 = dime TRES1;
  1602. ires1 = 0;
  1603. repe BRES1 nRES1; ires1 = ires1 + 1;
  1604. si (non (exis TRES1 ires1)); iter BRES1; finsi;
  1605. * titre
  1606. si (exis TRES1 . ires1 'TITRE');
  1607. tit1 = TRES1 . ires1 . 'TITRE';
  1608. sinon;
  1609. tit1 = chai 'Resultat' ires1*4;
  1610. tit1 = TRES1 . ires1 . 'TITRE';
  1611. finsi;
  1612. TRES1 . 'TITRE' . ires1 = tit1;
  1613. * couleur
  1614. si (exis TRES1 . ires1 'COULEUR');
  1615. coul1 = TRES1 . ires1 . 'COULEUR';
  1616. sinon;
  1617. icoul1 = icoul1 + 1; si (icoul1 > ncoul1); icoul1 = 1; finsi;
  1618. coul1 = extr mycoul1 icoul1;
  1619. TRES1 . ires1 . 'COULEUR' = coul1;
  1620. finsi;
  1621. ev1 = EVOL coul1 'MANU' 't' PRTIME
  1622. tit1 (TRES1 . ires1 . 'RESULTATS');
  1623. TRES1 . ires1 . 'RESULTATS_EVOL' = ev1;
  1624. evtot = evtot et ev1;
  1625. fin BRES1;
  1626. TRES1 . 'RESULTATS_EVOL' = evtot;
  1627. FINSI;
  1628.  
  1629. * on quitte on revenant sur la config initiale
  1630. SI (FLGDEP);
  1631. FORM GEOM0;
  1632. finsi;
  1633.  
  1634.  
  1635. FINPROC ;
  1636.  
  1637.  
  1638.  
  1639.  
  1640.  
  1641.  
  1642.  
  1643.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales