Télécharger continu.procedur

Retour à la liste

Numérotation des lignes :

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

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