Télécharger four2tri.procedur

Retour à la liste

Numérotation des lignes :

  1. * FOUR2TRI PROCEDUR BP208322 21/02/02 21:15:00 10876
  2.  
  3. *======================================================================
  4. *
  5. * FOUR2TRI (ou 4 2 3 pour les numericiens purs....) :
  6. * Genere un maillage 3D à partir d'un modele Fourier
  7. * Recombine les CHPO 2D Fourier sur ce maillage 3D
  8. * Creation : D. Combescure, SEMT/DYN, CEA Saclay, Decembre 2006
  9. * Modifs : B. Prabel, 2013/02(#8575), 2017/04(#9396), 2019/03, ...
  10. *
  11. *======================================================================
  12.  
  13. DEBPROC FOUR2TRI TAB1*'TABLE' NUMFOUR/'ENTIER' ;
  14.  
  15.  
  16. * opti debug 1;
  17. ************************************************************************
  18. * LECTURE DES DONNEES D'ENTREE
  19. ************************************************************************
  20. *
  21. * On ne cree le MAILLAGE_3D que s'il n'existe pas dejà
  22. *
  23. SI (NON (EXIS TAB1 'MAILLAGE_3D'));
  24. FLAGMAIL = VRAI;
  25. * il va falloir le creer -> a partir du modele ou d'un maillage 2d
  26. SI (EXIS TAB1 'MAILLAGE');
  27. MESH2D = TAB1 . 'MAILLAGE';
  28. SINON;
  29. SI (EXIS TAB1 'MODELE');
  30. * recup du maillage depuis le modele
  31. MESH2D = EXTR (TAB1 . 'MODELE') 'MAILLAGE';
  32. SINON;
  33. MESS 'Absence de MAILLAGE et de MODELE !';
  34. ERRE 21; QUIT FOUR2TRI;
  35. FINSI;
  36. FINSI;
  37. SINON;
  38. FLAGMAIL = FAUX;
  39. * le maillage existe deja -> on recupere tout
  40. MESH3D = TAB1 . 'MAILLAGE_3D' ;
  41. xy2D = TAB1 . 'COORDONNEES_2D';
  42. xy3D = TAB1 . 'COORDONNEES_3D';
  43. teta3D = TAB1 . 'ANGLE_3D' ;
  44. FINSI;
  45. *
  46. * On cree les champs de DEPLACEMENTS_3D si DEPLACEMENTS est fourni
  47. *
  48. FLAGDEPL = EXIS TAB1 'DEPLACEMENTS';
  49. *
  50. * On cree les champs de EFFORTS_3D si EFFORTS est fourni
  51. *
  52. FLAGFOR = EXIS TAB1 'EFFORTS';
  53. *
  54. * On cree les champs sclaires de EFFORTS_3D si EFFORTS est fourni
  55. *
  56. FLAGCHPS = EXIS TAB1 'CHPO_SYME';
  57. FLAGCHPA = EXIS TAB1 'CHPO_ANTI';
  58. *
  59. *
  60. * NUMFOUR = Harmonique de Fourier
  61. *
  62. *opti debug 1; list (TYPE NUMFOUR); list NUMFOUR;
  63. SI (NEG (TYPE NUMFOUR) 'ENTIER');
  64. * est-on bien en Fourier ?
  65. SI (EGA (VALE 'MODE') (MOT 'FOUR'));
  66. NUMFOUR = VALE 'MODE' 'FOUR';
  67. MESS 'Harmonique de Fourier non fournie :'
  68. ' on prend l harmonique courante' ' ' NUMFOUR ' !';
  69. SINON;
  70. * on n'est pas en Fourier -> on ne peut "que" faire le maillage 3D
  71. SI (FLAGDEPL OU FLAGFOR OU FLAGCHPS OU FLAGCHPA);
  72. ERRE 287; QUIT FOUR2TRI;
  73. FINSI;
  74. FINSI;
  75. FINSI;
  76. *mess 'four2tri: NUMFOUR='NUMFOUR;
  77. *
  78. * ANGLES
  79. *
  80. SI (EXIS TAB1 'ANGLES');
  81. PRANGL = TAB1 . 'ANGLES';
  82. SINON;
  83. MESS 'ANGLES non fourni : on prend tous les 15 degres par defaut !';
  84. PRANGL = PROG 0. PAS 15. 360.;
  85. TAB1 . 'ANGLES' = PRANGL;
  86. FINSI;
  87. * recup
  88. ANGL1 = EXTR PRANGL 1;
  89. NANGL = DIME PRANGL;
  90. ANGLN = EXTR PRANGL NANGL;
  91. * on impose la 1ere valeur a 0°
  92. SI (NEG ANGL1 0.);
  93. MESS 'Le premier angle doit etre 0. !';
  94. ERRE 21; QUIT FOUR2TRI;
  95. FINSI;
  96. * amplitude de l'angle de rotation
  97. ANG_ROTA = ANGLN - ANGL1;
  98. * on impose une progression constante et croissante
  99. DELTA_ANG = (ENLE PRANGL 1) - (ENLE PRANGL NANGL);
  100. MIN_D_ANG = MINI DELTA_ANG;
  101. MAX_D_ANG = MAXI DELTA_ANG;
  102. SI (MIN_D_ANG <EG 0.);
  103. MESS 'Les angles doivent suivre une progression croissante stricte !';
  104. ERRE 21; QUIT FOUR2TRI;
  105. FINSI;
  106. SI ((MAX_D_ANG - MIN_D_ANG) > 0.1);
  107. MESS 'Les angles doivent suivre une progression arithmetique !';
  108. ERRE 21; QUIT FOUR2TRI;
  109. FINSI;
  110. * revolution ?
  111. ZREVOLUTION = EGA ANG_ROTA 360.;
  112. SI (ZREVOLUTION); NROTA = NANGL - 1;
  113. SINON; NROTA = NANGL;
  114. FINSI;
  115.  
  116. *
  117. * Faut-il redresser le maillage et les champs tq (r,\theta,z) = (X,Y,Z)
  118. * ou laisse t'on (r,\theta,z) = (X,-Z,Y) ?
  119. *
  120. SI (EXIS TAB1 'REDRESSE');
  121. REDRESS = TAB1 . 'REDRESSE';
  122. SINON;
  123. REDRESS = FAUX;
  124. TAB1 . 'REDRESSE' = REDRESS;
  125. FINSI;
  126.  
  127. * points central, vecteur de rotation
  128. OPTI DIME 3 ;
  129. pcentre = 0. 0. 0.;
  130. SI REDRESS;
  131. vrota1 = 0. 0. 1.;
  132. vrota2 = 1. 0. 0.;
  133. SINON;
  134. vrota1 = 0. 1. 0.;
  135. FINSI;
  136.  
  137. * tolerance pour ELIM
  138. * rem : si xELIM = 0., alors pas d'ELIM, ce qui peut etre souhaitable
  139. * en cas de noeuds doubles volontaires par ex.)
  140. SI (EXIS TAB1 'ELIM');
  141. xELIM = TAB1 . 'ELIM';
  142. SINON;
  143. xELIM = 1.E-12 ;
  144. MESS 'Tolerance de ELIM non fourni : on prend' ' 'xELIM' par defaut !';
  145. TAB1 . 'ELIM' = xELIM;
  146. FINSI;
  147. FLAGELIM = xELIM '>' 0.;
  148. * xELIM etant utilise par ELIM et par MASQ et IPOL,
  149. * on met une autre valeur si nul
  150. SI (NON FLAGELIM); xELIM = VALE 'PETIT'; FINSI;
  151.  
  152.  
  153. ************************************************************************
  154. *
  155. * CREATION DU MAILLAGE 3D
  156. *
  157. ************************************************************************
  158.  
  159. SI FLAGMAIL;
  160.  
  161. * recup du maillage depuis le modele
  162. * fait+haut! MESH2D = EXTR MOD1 'MAILLAGE';
  163.  
  164. * on supprime les eventuels doublons (pour coque multicouche par ex.)
  165. MESH2D = UNIQ MESH2D;
  166.  
  167. * coordonnees 2D (R,Z) (logiquement, z2D = 0) avant redress !
  168. x2D y2D z2D = COOR MESH2D;
  169. xy2D = (NOMC 'R' x2d 'NATURE' 'DIFFUS')
  170. ET (NOMC 'Z' y2d 'NATURE' 'DIFFUS');
  171.  
  172. * on redresse si besoin
  173. SI REDRESS;
  174. DEPL MESH2D 'TOUR' 90. pcentre vrota2;
  175. FINSI;
  176. *
  177. * elements geometriques surfaciques, lineaires et quadratiques
  178. TYPSURF_LIN = MOTS 'TRI3' 'QUA4';
  179. TYPSURF_QUAD = MOTS 'TRI6' 'QUA8' 'TRI7' 'QUA9';
  180. TYPSURF = TYPSURF_LIN et TYPSURF_QUAD;
  181. TYPSEG = MOTS 'SEG2' 'SEG3';
  182. * split de MESH2D en MESH2D_1 et MESH2D_2 selon le type d'element
  183. TYP2D = ELEM MESH2D 'TYPE';
  184. MESH2D_1 = VIDE 'MAILLAGE';
  185. MESH2D_2 = VIDE 'MAILLAGE';
  186. REPE Bj (DIME TYP2D);
  187. DONE = FAUX;
  188. TYP2Dj = EXTR TYP2D &Bj;
  189. SI (EXIS TYPSURF TYP2Dj);
  190. MESH2D_1 = MESH2D_1 ET (ELEM MESH2D TYP2Dj);
  191. DONE = VRAI;
  192. FINSI;
  193. SI (EXIS TYPSEG TYP2Dj);
  194. MESH2D_2 = MESH2D_2 ET (ELEM MESH2D TYP2Dj);
  195. DONE = VRAI;
  196. FINSI;
  197. SI (NON DONE);
  198. MESS 'Elements de type' ' ' TYP2Dj ' non traite ...';
  199. * par exemple, les elements de type RACCORD ne sont pas taites ici
  200. * (tant mieux car ils ont des noeuds confondus)
  201. FINSI;
  202. FIN Bj;
  203. NEL_1 = NBEL MESH2D_1;
  204. NEL_2 = NBEL MESH2D_2;
  205.  
  206. * TODO: a ameliorer si ef quadratique ? pas sur que ce soit utile ...
  207. OPTI DIME 3 ELEM 'CUB8';
  208.  
  209. * infos utiles en cas d'appels successifs a FOUR2TRI
  210. TAB1 . 'COORDONNEES_2D' = xy2D;
  211.  
  212. * creation effective du maillage 3D
  213. MESH3D = VIDE 'MAILLAGE';
  214.  
  215. * VOLU ROTA pour les elements massifs
  216. SI (NEL_1 > 0);
  217. MESH3D_1 = VOLU 'ROTA' MESH2D_1 NROTA ANG_ROTA pcentre vrota1;
  218. SI FLAGELIM;
  219. ELIM MESH3D_1 xELIM;
  220. MESH3D_1 = REGE MESH3D_1;
  221. FINSI;
  222. MESH3D = MESH3D ET MESH3D_1;
  223. FINSI;
  224.  
  225. * ROTA pour les elements surfaciques
  226. SI (NEL_2 > 0);
  227. MESH3D_2 = ROTA MESH2D_2 NROTA ANG_ROTA pcentre vrota1;
  228. SI FLAGELIM;
  229. ELIM MESH3D_2 xELIM;
  230. MESH3D_2 = REGE MESH3D_2;
  231. FINSI;
  232. MESH3D = MESH3D ET MESH3D_2;
  233. FINSI;
  234.  
  235. * les 1ers ELIM servaient aux noeuds en r=0 et en theta=0=360 degres
  236. * celui-ci aux noeuds communs a MESH3D_1 et MESH3D_2
  237. SI FLAGELIM;
  238. SI ((NEL_1 * NEL_2) > 0);
  239. ELIM MESH3D_1 MESH3D_2 xELIM;
  240. FINSI;
  241. FINSI;
  242.  
  243. * On range le(s) nouveau(x) maillage(s) 3D
  244. TAB1 . 'MAILLAGE_3D_MASSIF' = MESH3D_1;
  245. TAB1 . 'MAILLAGE_3D_COQUE' = MESH3D_2;
  246. TAB1 . 'MAILLAGE_3D' = MESH3D;
  247.  
  248. * Reste a calculer les coordonnees 3D
  249. * pour teta, on evite les pb en r=0 avec un masque
  250. x3D y3D z3D = COOR MESH3D;
  251. SI REDRESS;
  252. r3D = ((x3D**2) + (y3D**2))**0.5;
  253. xy3D = (NOMC 'R' r3D) ET (NOMC 'Z' z3d);
  254. mask0 = MASQ r3D 'EGINFE' xELIM;
  255. teta3D = ATG y3D (x3D + mask0);
  256. SINON;
  257. r3D = ((x3D**2) + (z3D**2))**0.5;
  258. xy3D = (NOMC 'R' r3D) ET (NOMC 'Z' y3d);
  259. mask0 = MASQ r3D 'EGINFE' xELIM;
  260. teta3D = ATG (-1.*z3D) (x3D + mask0);
  261. FINSI;
  262. TAB1 . 'COORDONNEES_3D' = xy3D;
  263. TAB1 . 'ANGLE_3D' = teta3D;
  264.  
  265. FINSI;
  266.  
  267.  
  268. ************************************************************************
  269. *
  270. * CHPO VECTEURS : CHAMPS DE | DEPLACEMENTS
  271. * | EFFORTS
  272. *
  273. ************************************************************************
  274.  
  275. * on evite du copier-coller en introduisant la boucle BITER
  276. NITER =0;
  277. SI (FLAGDEPL ET FLAGFOR);
  278. NITER = 2;
  279. FLAGDEPL = VRAI;
  280. FLAGFOR = FAUX;
  281. SINON;
  282. SI (FLAGDEPL OU FLAGFOR);
  283. NITER =1;
  284. FINSI;
  285. FINSI;
  286.  
  287. REPE BITER NITER;
  288.  
  289. OPTI DIME 3 MODE TRID;
  290. SI FLAGDEPL;
  291. ind2D = mot 'DEPLACEMENTS';
  292. ind3D = mot 'DEPLACEMENTS_3D';
  293. monature= mot 'DIFFUS';
  294. mo2d_1 = mot 'UR';
  295. mo2d_2 = mot 'UT';
  296. mo2d_3 = mot 'UZ';
  297. mo2d_4 = mot 'IUR';
  298. mo2d_5 = mot 'IUT';
  299. mo2d_6 = mot 'IUZ';
  300. mo3d_1 = mot 'UX';
  301. mo3d_2 = mot 'UY';
  302. mo3d_3 = mot 'UZ';
  303. FINSI;
  304. SI FLAGFOR;
  305. ind2D = mot 'EFFORTS';
  306. ind3D = mot 'EFFORTS_3D';
  307. monature= mot 'DISCRET';
  308. mo2d_1 = mot 'FR';
  309. mo2d_2 = mot 'FT';
  310. mo2d_3 = mot 'FZ';
  311. mo2d_4 = mot 'IFR';
  312. mo2d_5 = mot 'IFT';
  313. mo2d_6 = mot 'IFZ';
  314. mo3d_1 = mot 'FX';
  315. mo3d_2 = mot 'FY';
  316. mo3d_3 = mot 'FZ';
  317. FINSI;
  318. SI (EXIS TAB1 ind3D);
  319. MESS 'L indice . ' ind3D 'existe deja : on l ecrase !';
  320. SINON;
  321. TAB1 . ind3D = TABL;
  322. FINSI;
  323.  
  324. * il nous faut les cos() et sin() de \theta et n*\theta
  325. cosA1 = cos teta3D;
  326. sinA1 = sin teta3D;
  327. cosNA = cos (NUMFOUR * teta3D);
  328. sinNA = sin (NUMFOUR * teta3D);
  329.  
  330. *-----------------------------> BOUCLE SUR LES DEPLACEMENTS
  331. NDEP = DIME (TAB1 . ind2D);
  332. REPETER LAB2 NDEP; i2 = &lab2;
  333.  
  334. DEP2D = TAB1 . ind2D . i2;
  335. NUA2D = NUAG ((CHAN 'ATTRIBUT' xy2D 'NATURE' monature) ET DEP2D);
  336. DEP3DFOU = IPOL NUA2D xy3D 'PID' 1. 'ELIM' xELIM;
  337.  
  338. * serie de Fourier
  339. si (ega NUMFOUR 0);
  340. * ici, on considere seulement les ddls symetriques pour n=0
  341. DEP3D0 = (EXCO DEP3DFOU mo2D_1 'NOID' mo2D_1)
  342. + (EXCO DEP3DFOU mo2D_2 'NOID' mo2D_2)
  343. + (EXCO DEP3DFOU mo2d_3 'NOID' mo2d_3);
  344. sinon;
  345. DEP3D0 = ((EXCO DEP3DFOU mo2D_1 'NOID' mo2D_1)*cosNA)
  346. + ((EXCO DEP3DFOU mo2D_2 'NOID' mo2D_2)*sinNA)
  347. + ((EXCO DEP3DFOU mo2d_3 'NOID' mo2d_3)*cosNA)
  348. + ((EXCO DEP3DFOU mo2D_4 'NOID' mo2D_4)*sinNA)
  349. + ((EXCO DEP3DFOU mo2D_5 'NOID' mo2D_5)*cosNA)
  350. + ((EXCO DEP3DFOU mo2D_6 'NOID' mo2D_6)*sinNA);
  351. finsi;
  352. * changement de repere
  353. si REDRESS;
  354. DEP3D1 = ((EXCO DEP3D0 mo2D_1 'NOID' mo3D_1)*cosA1)
  355. + ((EXCO DEP3D0 mo2D_1 'NOID' mo3D_2)*sinA1)
  356. - ((EXCO DEP3D0 mo2D_2 'NOID' mo3D_1)*sinA1)
  357. + ((EXCO DEP3D0 mo2D_2 'NOID' mo3D_2)*cosA1)
  358. + (EXCO DEP3D0 mo2d_3 'NOID' mo3d_3);
  359. DEP3D1 = DEP3D1
  360. + ((EXCO DEP3D0 mo2D_4 'NOID' mo3D_1)*cosA1)
  361. + ((EXCO DEP3D0 mo2D_4 'NOID' mo3D_2)*sinA1)
  362. - ((EXCO DEP3D0 mo2D_5 'NOID' mo3D_1)*sinA1)
  363. + ((EXCO DEP3D0 mo2D_5 'NOID' mo3D_2)*cosA1)
  364. + (EXCO DEP3D0 mo2D_6 'NOID' mo3d_3);
  365. sinon;
  366. DEP3D1 = ((EXCO DEP3D0 mo2D_1 'NOID' mo3D_1)*cosA1)
  367. - ((EXCO DEP3D0 mo2D_1 'NOID' mo3d_3)*sinA1)
  368. - ((EXCO DEP3D0 mo2D_2 'NOID' mo3D_1)*sinA1)
  369. - ((EXCO DEP3D0 mo2D_2 'NOID' mo3d_3)*cosA1)
  370. + (EXCO DEP3D0 mo2d_3 'NOID' mo3D_2);
  371. DEP3D1 = DEP3D1
  372. + ((EXCO DEP3D0 mo2D_4 'NOID' mo3D_1)*cosA1)
  373. - ((EXCO DEP3D0 mo2D_4 'NOID' mo3d_3)*sinA1)
  374. - ((EXCO DEP3D0 mo2D_5 'NOID' mo3D_1)*sinA1)
  375. - ((EXCO DEP3D0 mo2D_5 'NOID' mo3d_3)*cosA1)
  376. + (EXCO DEP3D0 mo2D_6 'NOID' mo3D_2);
  377. finsi;
  378.  
  379. * on stocke le nouveau deplacement 3D
  380. TAB1 . ind3D . i2 = CHAN 'ATTRIBUT' DEP3D1 'NATURE' monature;
  381.  
  382. FIN LAB2;
  383. *-----------------------------> FIN DE BOUCLE SUR LES DEPLACEMENTS
  384. *
  385. SI (EGA NITER 2);
  386. SI (EGA &BITER 1);
  387. * pour la prochaine ITER
  388. FLAGDEPL=FAUX;
  389. FLAGFOR =VRAI;
  390. SINON;
  391. * on sort en remettant les valeurs par acquis de conscience
  392. FLAGDEPL=VRAI;
  393. FLAGFOR =VRAI;
  394. FINSI;
  395. FINSI;
  396.  
  397. FIN BITER;
  398.  
  399.  
  400.  
  401. ************************************************************************
  402. *
  403. * CHPO SCALAIRE (ex : PRESSION, EPAISSEUR ...) : CAS |SYMETRIQUE
  404. * |ANTI-SYMETRIQUE
  405. *
  406. ************************************************************************
  407.  
  408. * on evite du copier-coller en introduisant la boucle BITER
  409. NITER =0;
  410. SI (FLAGCHPS ET FLAGCHPA);
  411. NITER = 2;
  412. FLAGCHPS = VRAI;
  413. FLAGCHPA = FAUX;
  414. SINON;
  415. SI (FLAGCHPS OU FLAGCHPA);
  416. NITER =1;
  417. FINSI;
  418. FINSI;
  419.  
  420. REPE BITER NITER;
  421.  
  422. OPTI DIME 3 MODE TRID;
  423. * il nous faut les cos() ou sin() de n*\theta
  424. SI FLAGCHPS;
  425. ind2D = mot 'CHPO_SYME';
  426. ind3D = mot 'CHPO_SYME_3D';
  427. cosinNA = cos (NUMFOUR * teta3D);
  428. FINSI;
  429. SI FLAGCHPA;
  430. ind2D = mot 'CHPO_ANTI';
  431. ind3D = mot 'CHPO_ANTI_3D';
  432. cosinNA = sin (NUMFOUR * teta3D);
  433. FINSI;
  434. SI (EXIS TAB1 ind3D);
  435. MESS 'L indice . ' ind3D 'existe deja : on l ecrase !';
  436. SINON;
  437. TAB1 . ind3D = TABL;
  438. FINSI;
  439.  
  440. *-----------------------------> BOUCLE SUR LES CHPO
  441. NCHPS = DIME (TAB1 . ind2D);
  442. REPETER LAB3 NCHPS; i3 = &lab3;
  443.  
  444. CHP2D = TAB1 . ind2D . i3 ;
  445. LCOMPU = EXTR CHP2D 'COMP';
  446. CHP2D = EXCO CHP2D LCOMPU NUMFOUR;
  447. NUA2D = NUAG (xy2D ET CHP2D);
  448. CHP3D = IPOL NUA2D xy3D 'PID' 1. 'ELIM' xELIM;
  449.  
  450. * serie de Fourier
  451. CHP3D = CHP3D * cosinNA;
  452. CHP3D = CHAN 'ATTRIBUT' CHP3D 'NATURE' 'DIFFUS';
  453. * pas de changement de repere pour des champs scalaires
  454.  
  455. * on stocke le champs 3D
  456. TAB1 . ind3D . i3 = CHP3D;
  457.  
  458. FIN LAB3;
  459. *-----------------------------> FIN DE BOUCLE SUR LES CHPO
  460.  
  461. SI (EGA NITER 2);
  462. SI (EGA &BITER 1);
  463. * pour la prochaine ITER
  464. FLAGCHPS=FAUX;
  465. FLAGCHPA=VRAI;
  466. SINON;
  467. * on sort en remettant les valeurs par acquis de conscience
  468. FLAGCHPS=VRAI;
  469. FLAGCHPA=VRAI;
  470. FINSI;
  471. FINSI;
  472.  
  473. FIN BITER;
  474.  
  475.  
  476.  
  477. ************************************************************************
  478. * TODO : IL MANQUE LES CHAMPS DE CONTRAINTES ET DE DEFORMATIONS.....
  479. ************************************************************************
  480.  
  481.  
  482. 'FINPROC';
  483. *
  484.  
  485.  
  486.  
  487.  
  488.  
  489.  
  490.  
  491.  

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