Télécharger toposurf.procedur

Retour à la liste

Numérotation des lignes :

  1. * TOPOSURF PROCEDUR FD218221 21/04/14 21:15:09 10965
  2.  
  3. ************************************************************************
  4. ** Extraction procedure of a smoothed surface from a topology
  5. **
  6. ** Author:
  7. ** Guenhael Le Quilliec (LaMe - Polytech Tours)
  8. **
  9. ** Version:
  10. ** 3.1 2020/04/13 Updated to make it compatible with thermic models
  11. ** and fix an issue when OPTI ELEM is not defined
  12. ** 3.0 2018/01/29 Updated to make it independant of TOPOPTIM version,
  13. ** fix some issues with PYR5 and PRI6 elements,
  14. ** compatibility with quadratic elements,
  15. ** the model is no longer an input argument,
  16. ** all input arguments are now provided in a table
  17. ** 2.0 2018/01/26 Updated to make it compatible with TOPOPTIM V2.1
  18. ** 1.0 2016/12/13 Original version compatible with TOPOPTIM V1.0
  19. ************************************************************************
  20.  
  21. DEBP TOPOSURF tab0*'TABLE' ;
  22.  
  23. **********************************************************************
  24. * INPUT *
  25. **********************************************************************
  26.  
  27. * Topology
  28. Topo0 = tab0.'TOPOLOGIE' ;
  29.  
  30. * Model
  31. Mod0 = tab0.'MODELE' ;
  32.  
  33. * Filter rate
  34. SI (NON (EXIS tab0 'TAUX_FILTRAGE')) ;
  35. tab0.'TAUX_FILTRAGE' = 1 ;
  36. FINS ;
  37. Filter0 = tab0.'TAUX_FILTRAGE' ;
  38.  
  39. * Isovalue
  40. SI (NON (EXIS tab0 'ISOVALEUR')) ;
  41. tab0.'ISOVALEUR' = 0.5 ;
  42. FINS ;
  43. Iso0 = tab0.'ISOVALEUR' ;
  44.  
  45. * Default minimum size of the output surface elements
  46. SI (NON (EXIS tab0 'ELIMINATION')) ;
  47. tab0.'ELIMINATION' = 1.0e-6 ;
  48. FINS ;
  49. Elim0 = tab0.'ELIMINATION' ;
  50.  
  51. * Orientation of the output surface will be applied by default
  52. SI (NON (EXIS tab0 'ORIENTATION')) ;
  53. tab0.'ORIENTATION' = VRAI ;
  54. FINS ;
  55. Orie0 = tab0.'ORIENTATION' ;
  56.  
  57. * Default value of the extra thickness
  58. * This value should be very low compared to the element size
  59. * when the domain is not convex.
  60. * If not, Enve1 might be wrong and GRAD won't work
  61. SI (NON (EXIS tab0 'SUREPAISSEUR')) ;
  62. tab0.'SUREPAISSEUR' = 1.0e-3 ;
  63. FINS ;
  64. Shift0 = tab0.'SUREPAISSEUR' ;
  65.  
  66. * Default value of the thickness of the output mesh when the input
  67. * model is 2D
  68. SI (NON (EXIS tab0 'EPAISSEUR')) ;
  69. tab0.'EPAISSEUR' = 1.0 ;
  70. FINS ;
  71. Thick0 = tab0.'EPAISSEUR' ;
  72.  
  73. **********************************************************************
  74. * PREPROCESSING *
  75. **********************************************************************
  76.  
  77. * Mesh
  78. Msh0 = EXTR Mod0 'MAIL' ;
  79.  
  80. * Types of elements
  81. Type0 = Msh0 ELEM 'TYPE' ;
  82.  
  83. * If the model is mechanical
  84. SI (EXIS Mod0 'FORM' 'MECANIQUE') ;
  85. * Create a random material
  86. Mat0 = MATE Mod0 'RHO' 1.0 ;
  87. * Unitary field
  88. un1 = MANU 'CHML' Mod0 'SCAL' 1.0 'TYPE' 'SCALAIRE' 'GRAVITE' ;
  89. * Volume of each element
  90. volEl1 = INTG 'ELEM' un1 Mod0 Mat0 ;
  91. SINO ;
  92. * Random mechanical model
  93. Mod0 = MODE Msh0 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' ;
  94. * Create a random material
  95. Mat0 = MATE Mod0 'RHO' 1.0 ;
  96. * Unitary field
  97. un1 = MANU 'CHML' Mod0 'SCAL' 1.0 'TYPE' 'SCALAIRE' 'GRAVITE' ;
  98. * Volume of each element
  99. volEl1 = INTG 'ELEM' un1 Mod0 Mat0 ;
  100. * Express the input topology element field at the barycenters of the
  101. * elements of the model Mod0
  102. Topo0 = TOPOCHAN Topo0 Mod0 (volEl1 POIN 'SUPE' -1.0) ;
  103. Topo0 = CHAN 'TYPE' Topo0 'SCALAIRE' ;
  104. FINS ;
  105.  
  106. * Make sure that the topology field is not composed of zones where
  107. * the values are constant so that to avoid problems with the function
  108. * POIN 'SUPE'
  109. Topo0 = (INTG 'ELEM' Topo0 Mod0 Mat0) / volEl1 ;
  110.  
  111. * If the mesh is quadratic
  112. SI (EXIS Type0 (MOTS 'TRI6' 'QUA8' 'CU20' 'TE10' 'PR15' 'PY13') 'OU');
  113. * Change it to linear
  114. Msh0 = CHAN 'LINE' Msh0 ;
  115. * Update types of elements
  116. Type0 = Msh0 ELEM 'TYPE' ;
  117. * Create a random model
  118. Mod0 = MODE Msh0 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' ;
  119. * Update the random material
  120. Mat0 = MATE Mod0 'YOUN' 1.0 'NU' 0.3 ;
  121. * Update the unitary field
  122. un1 = MANU 'CHML' Mod0 'SCAL' 1.0 'TYPE' 'SCALAIRE' 'GRAVITE' ;
  123. * Update the volume of each element
  124. volEl1 = INTG 'ELEM' un1 Mod0 Mat0 ;
  125. * Express the input topology element field at the barycenters of the
  126. * elements of the model Mod0
  127. Topo0 = TOPOCHAN Topo0 Mod0 (volEl1 POIN 'SUPE' -1.0) ;
  128. Topo0 = CHAN 'TYPE' Topo0 'SCALAIRE' ;
  129. FINS ;
  130.  
  131. * Filtering the topology field
  132. SI (Filter0 > 0) ;
  133. REPE Loop1 Filter0 ;
  134. Topo0 = CHAN 'CHAM' (CHAN 'CHPO' Mod0 Topo0 'MOYE')
  135. Mod0 'GRAVITE' ;
  136. FIN Loop1 ;
  137. FINS ;
  138.  
  139. * Node topology field
  140. Ch0 = CHAN 'CHPO' Topo0 Mod0 'MOYE' ;
  141. Ch0 = Ch0 * (1.0 / (MAXI Ch0)) ;
  142.  
  143. * Save the options 'DIME', 'MODE', 'ELEM'
  144. Dime0 = VALE 'DIME' ;
  145. Mode0 = VALE 'MODE' ;
  146. ElemOpt0 = VALE 'ELEM' ;
  147.  
  148. * 2D Topology
  149. SI (EXIS Type0 (MOTS 'TRI3' 'QUA4') 'OU') ;
  150. * TODO Make it compatible with 3D in case of shell elements
  151. * Switch to 2D
  152. OPTI 'DIME' 2 'MODE' 'PLAN' 'CONT' ;
  153. * If the reasult has to stay in 2D
  154. SI ((ABS Thick0) < 1.0e-9) ;
  155. Type3D = FAUX ;
  156. * Change option 'ELEM' to TRI3
  157. OPTI 'ELEM' TRI3 ;
  158. * Change the element type of the mesh
  159. Msh0 = CHAN 'TRI3' Msh0 ;
  160. * Domain boundary
  161. Cont0 = CONT Msh0 ;
  162. * Shifted boundary
  163. Norm0 = PRES 'MASS' Mod0 (-1.0 * Shift0) Cont0 ;
  164. Norm0 = CHAN 'COMP' (MOTS 'FX' 'FY') (MOTS 'UX' 'UY')
  165. Norm0 ;
  166. Cont1 = Cont0 PLUS Norm0 ;
  167. * External layer of elements
  168. Layer0 = Cont0 REGL 1 Cont1 ;
  169. * Add a field for the shifted boundary
  170. Ch0 = Ch0 ET
  171. (MANU 'CHPO' Cont1 1 'SCAL' -1.0e3 'NATURE' 'DIFFUS') ;
  172. * New model
  173. Mod1 = MODE (Msh0 ET Layer0) 'MECANIQUE' 'ELASTIQUE' ;
  174. * Element field
  175. MCh0 = CHAN 'CHAM' Ch0 Mod1 'NOEUD' ;
  176. * Extract the iso surface
  177. Surf0 = (ISOV MCh0 Iso0) ELEM SEG2 ;
  178. * Fuse the nodes
  179. ELIM Surf0 Elim0 ;
  180. * If the surface has to be oriented
  181. SI Orie0 ;
  182. * Fuse any superposed element
  183. Surf1 = TABL ;
  184. REPE Loop1 (NBEL Surf0) ;
  185. * Considered element
  186. Elem0 = Surf0 ELEM &Loop1 ;
  187. P0 = CHAN POI1 Elem0 ;
  188. SI (EGA (NBEL P0) 2) ;
  189. N1 = NOEU (P0 POIN 1) ;
  190. N2 = NOEU (P0 POIN 2) ;
  191. SI (N1 > N2) ;
  192. N0 = N2 ;
  193. N2 = N1 ;
  194. N1 = N0 ;
  195. FINS ;
  196. Surf1.(CHAI N1 '-' N2) = Elem0 ;
  197. FINS ;
  198. FIN Loop1 ;
  199. Ind0 = INDE Surf1 ;
  200. Surf0 = Surf1.(Ind0.(DIME Ind0)) ;
  201. REPE Loop1 ((DIME Ind0) - 1) ;
  202. Surf0 = Surf0 ET Surf1.(Ind0.&Loop1) ;
  203. FIN Loop1 ;
  204. * Create the gradient field giving the orientation
  205. Ch0 = (REDU Ch0 Msh0) ET (MANU 'CHPO' Cont1 1 'SCAL'
  206. -1.0 'NATURE' 'DIFFUS') ;
  207. Mod1 = MODE (Msh0 ET Layer0) 'THERMIQUE' ;
  208. Grad0 = -1.0 * (GRAD Mod1 (CHAN 'COMP' 'T' Ch0)) ;
  209. Grad0 = CHAN 'NOEUD' Mod1 Grad0 ;
  210. * Centers of mass
  211. REPE Loop1 (NBEL Surf0) ;
  212. Elem0 = Surf0 ELEM &Loop1 ;
  213. SI (&Loop1 NEG 1) ;
  214. Pts0 = Pts0 ET (BARY Elem0) ;
  215. SINO ;
  216. Pts0 = BARY Elem0 ;
  217. FINS ;
  218. FIN Loop1 ;
  219. * Projection of the gradient on the centers of mass
  220. Grad0 = PROI Pts0 Grad0 ;
  221. * Initialization of the loop on each elements
  222. Inv0 = (0.0 0.0) ET (0.0 0.0) ;
  223. * Loop on each element of the surface
  224. REPE Loop1 (NBEL Surf0) ;
  225. * Considered element
  226. Elem0 = Surf0 ELEM &Loop1 ;
  227. * Compute the normal direction of the element
  228. P0 = CHAN POI1 Elem0 ;
  229. P1 = P0 POIN 1 ;
  230. P2 = P0 POIN 2 ;
  231. N0 = (P2 MOIN P1) TOUR 90.0 (0.0 0.0) ;
  232. * Center of mass
  233. Bary0 = Pts0 POIN &Loop1 ;
  234. * Desired orientation
  235. G1 = EXTR Grad0 'T,X' Bary0 ;
  236. G2 = EXTR Grad0 'T,Y' Bary0 ;
  237. G0 = (G1 G2) ;
  238. * Add the element to inverse
  239. SI ((N0 PSCA G0) < 0.0) ;
  240. Inv0 = Inv0 ET Elem0 ;
  241. FINS ;
  242. FIN Loop1 ;
  243. * If there are some elements to invert
  244. SI ((NBEL Inv0) > 2) ;
  245. Inv0 = Inv0 ELEM SEG2 ;
  246. Surf0 = (DIFF Surf0 Inv0) ET (INVE Inv0) ;
  247. FINS ;
  248. FINS ;
  249. * Give information about the output mesh quality
  250. VERM Surf0 ;
  251. SINO ;
  252. Type3D = VRAI ;
  253. * Force the option 'DIME' to 3D
  254. Dime0 = 3 ;
  255. * Switch to 3D
  256. OPTI 'DIME' 3 ;
  257. * Change option 'ELEM' to CUB8
  258. OPTI 'ELEM' CUB8 ;
  259. * Create the second side of the 3D mesh from the 2D mesh
  260. Recto0 ChR0 = Msh0 Ch0 PLUS (0.0 0.0 Thick0) ;
  261. * Create the 3D mesh
  262. Msh0 = Msh0 VOLU 1 Recto0 ;
  263. * Create the 3D topology field
  264. Ch0 = Ch0 ET ChR0 ;
  265. * Create the 3D model
  266. Mod0 = MODE Msh0 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' ;
  267. FINS ;
  268. SINO ;
  269. Type3D = VRAI ;
  270. FINS ;
  271.  
  272. * 3D Topology
  273. SI Type3D ;
  274. * Change option 'ELEM' to CUB8
  275. OPTI 'ELEM' CUB8 ;
  276. * Domain boundary
  277. Enve0 = ENVE Msh0 ;
  278. * Normal direction on the domain boundary
  279. Norm0 = PRES 'MASS' Mod0 (-1.0 * Shift0) Enve0 ;
  280. Norm0 = CHAN 'COMP' (MOTS 'FX' 'FY' 'FZ')
  281. (MOTS 'UX' 'UY' 'UZ') Norm0 ;
  282. * Partitionate the surfaces
  283. NbPart0 = 32 ;
  284. SI (NbPart0 > (NBEL Enve0)) ;
  285. NbPart0 = 16 ;
  286. FINS ;
  287. SI (NbPart0 > (NBEL Enve0)) ;
  288. NbPart0 = 8 ;
  289. FINS ;
  290. PEnve0 = PART 'OPTI' NbPart0 Enve0 ;
  291. * Since PART does not always respect NbPart0, we redefine it
  292. NbPart0 = DIME PEnve0 ;
  293. REPE Loop1 NbPart0 ;
  294. SI (EXIS PEnve0 NbPart0) ;
  295. QUIT Loop1 ;
  296. FINS ;
  297. NbPart0 = NbPart0 - 1 ;
  298. FIN Loop1 ;
  299. * External layer of elements
  300. Layer0 = TABL 'ESCLAVE' ;
  301. * Loop on each entry of PEnve0
  302. REPE Loop1 NbPart0 ;
  303. PEnve1 = PEnve0.&Loop1 PLUS (REDU Norm0 PEnve0.&Loop1) ;
  304. Layer0.&Loop1 = PEnve0.&Loop1 VOLU 1 PEnve1 ;
  305. SI (&Loop1 > 1) ;
  306. Enve1 = Enve1 ET PEnve1 ;
  307. ELIM Enve1 PEnve1 Elim0 ;
  308. SINO ;
  309. Enve1 = PEnve1 ;
  310. FINS ;
  311. FIN Loop1 ;
  312. Layer0 = ET Layer0 ;
  313. * Add a field for the shifted boundary
  314. Ch0 = Ch0 ET
  315. (MANU 'CHPO' Enve1 1 'SCAL' -1.0e3 'NATURE' 'DIFFUS') ;
  316. * Create a mesh that will be compatible with ISOV
  317. IsoMsh0 = Msh0 ET Layer0 ;
  318. Type0 = IsoMsh0 ELEM 'TYPE' ;
  319. IsoMail1 = TABL 'ESCLAVE' ;
  320. REPE Loop1 (DIME Type0) ;
  321. Type1 = EXTR Type0 &Loop1 ;
  322. SI ((EGA Type1 'TET4') OU (EGA Type1 'CUB8')) ;
  323. IsoMail1.&Loop1 = IsoMsh0 ELEM Type1 ;
  324. SINO ;
  325. * Using CHAN 'TET4' is not a good option as the edges
  326. * of the new elements won't necessary share the edges
  327. * of the neigbor elements and this may create holes
  328. * after using ISOV.
  329. * IsoMail1.&Loop1 = CHAN TET4 (IsoMsh0 ELEM Type1) ;
  330. * We change incompatible elements with ISOV by CUB8
  331. * ****
  332. Volu0 = IsoMsh0 ELEM Type1 ;
  333. Volu1 = VIDE 'MAILLAGE'/'CUB8' ;
  334. SI (EGA Type1 PRI6) ;
  335. * Loop on each element
  336. REPE Loop2 (NBEL Volu0) ;
  337. * Considered element
  338. Elem0 = Volu0 ELEM &Loop2 ;
  339. * Nodes of the element
  340. P0 = CHAN POI1 Elem0 ;
  341. P1 = P0 POIN 1 ;
  342. P2 = P0 POIN 2 ;
  343. P3 = P0 POIN 3 ;
  344. P4 = P0 POIN 4 ;
  345. P5 = P0 POIN 5 ;
  346. P6 = P0 POIN 6 ;
  347. * Manualy change this PRI6 by a CUB8
  348. Volu1 = Volu1 ET
  349. (MANU CUB8 P1 P2 P3 P3 P4 P5 P6 P6) ;
  350. FIN Loop2 ;
  351. FINS ;
  352. SI (EGA Type1 PYR5) ;
  353. * Loop on each element
  354. REPE Loop2 (NBEL Volu0) ;
  355. * Considered element
  356. Elem0 = Volu0 ELEM &Loop2 ;
  357. * Nodes of the element
  358. P0 = CHAN POI1 Elem0 ;
  359. P1 = P0 POIN 1 ;
  360. P2 = P0 POIN 2 ;
  361. P3 = P0 POIN 3 ;
  362. P4 = P0 POIN 4 ;
  363. P5 = P0 POIN 5 ;
  364. * Manualy change this PYR5 by a CUB8
  365. Volu1 = Volu1 ET
  366. (MANU CUB8 P1 P2 P3 P4 P5 P5 P5 P5) ;
  367. FIN Loop2 ;
  368. FINS ;
  369. IsoMail1.&Loop1 = Volu1 ;
  370. * ****
  371. FINS ;
  372. FIN Loop1 ;
  373. IsoMsh0 = ET IsoMail1 ;
  374. * New model
  375. Mod1 = MODE IsoMsh0 'MECANIQUE' 'ELASTIQUE' ;
  376. * Element field
  377. MCh0 = CHAN 'CHAM' Ch0 Mod1 'NOEUD' ;
  378. * Extract the iso surface
  379. Surf0 = ISOV MCh0 Iso0 ;
  380. * Change the surface to TRI3
  381. Type0 = Surf0 ELEM 'TYPE' ;
  382. Surf1 = TABL 'ESCLAVE' ;
  383. REPE Loop1 (DIME Type0) ;
  384. Type1 = EXTR Type0 &Loop1 ;
  385. SI (EGA Type1 'TRI3') ;
  386. Surf1.&Loop1 = Surf0 ELEM Type1 ;
  387. SINO ;
  388. SI (NEG Type1 'POI1') ;
  389. Surf1.&Loop1 = CHAN 'TRI3' (Surf0 ELEM Type1) ;
  390. FINS ;
  391. FINS ;
  392. FIN Loop1 ;
  393. Surf0 = ET Surf1 ;
  394. * Fuse the nodes
  395. ELIM Surf0 Elim0 ;
  396. * Remove any elements with fused nodes
  397. Elem1 = LECT ;
  398. REPE Loop1 (NBEL Surf0) ;
  399. SI (NEG (NBNO (Surf0 ELEM &Loop1)) 3) ;
  400. Elem1 = Elem1 ET (LECT &Loop1) ;
  401. FINS ;
  402. FIN Loop1 ;
  403. Elem0 = LECT 1 'PAS' 1 (NBEL Surf0) ;
  404. SI ((DIME Elem1) > 0) ;
  405. Surf0 = Surf0 ELEM (ENLE Elem0 Elem1) ;
  406. FINS ;
  407. * If the surface has to be oriented
  408. SI Orie0 ;
  409. * Fuse any superposed element
  410. Surf1 = TABL ;
  411. REPE Loop1 (NBEL Surf0) ;
  412. * Considered element
  413. Elem0 = Surf0 ELEM &Loop1 ;
  414. P0 = CHAN POI1 Elem0 ;
  415. SI (EGA (NBEL P0) 3) ;
  416. N1 = NOEU (P0 POIN 1) ;
  417. N2 = NOEU (P0 POIN 2) ;
  418. N3 = NOEU (P0 POIN 3) ;
  419. N0 = ORDO (LECT N1 N2 N3) 'CROI' ;
  420. Surf1.(CHAI (EXTR N0 1) '-' (EXTR N0 2)
  421. '-' (EXTR N0 3)) = Elem0 ;
  422. FINS ;
  423. FIN Loop1 ;
  424. Ind0 = INDE Surf1 ;
  425. Surf0 = Surf1.(Ind0.(DIME Ind0)) ;
  426. REPE Loop1 ((DIME Ind0) - 1) ;
  427. Surf0 = Surf0 ET Surf1.(Ind0.&Loop1) ;
  428. FIN Loop1 ;
  429. * Create the gradient field giving the orientation
  430. Ch0 = (REDU Ch0 Msh0) ET
  431. (MANU 'CHPO' Enve1 1 'SCAL' -1.0 'NATURE' 'DIFFUS') ;
  432. Mod1 = MODE (Msh0 ET Layer0) 'THERMIQUE' ;
  433. Grad0 = -1.0 * (GRAD Mod1 (CHAN 'COMP' 'T' Ch0)) ;
  434. Grad0 = CHAN 'NOEUD' Mod1 Grad0 ;
  435. * Centers of mass
  436. REPE Loop1 (NBEL Surf0) ;
  437. Elem0 = Surf0 ELEM &Loop1 ;
  438. SI (&Loop1 NEG 1) ;
  439. Pts0 = Pts0 ET (BARY Elem0) ;
  440. SINO ;
  441. Pts0 = BARY Elem0 ;
  442. FINS ;
  443. FIN Loop1 ;
  444. * Projection of the gradient on the centers of mass
  445. Grad0 = PROI Pts0 Grad0 ;
  446. * Initialization of the loop on each elements
  447. Inv0 = (0.0 0.0 0.0) ET (0.0 0.0 0.0) ;
  448. * Loop on each element of the surface
  449. REPE Loop1 (NBEL Surf0) ;
  450. * Considered element
  451. Elem0 = Surf0 ELEM &Loop1 ;
  452. * Compute the normal direction of the element
  453. P0 = CHAN POI1 Elem0 ;
  454. P1 = P0 POIN 1 ;
  455. P2 = P0 POIN 2 ;
  456. P3 = P0 POIN 3 ;
  457. N0 = (P2 MOIN P1) PVEC (P3 MOIN P2) ;
  458. * Center of mass
  459. Bary0 = Pts0 POIN &Loop1 ;
  460. * Desired orientation
  461. G1 = EXTR Grad0 'T,X' Bary0 ;
  462. G2 = EXTR Grad0 'T,Y' Bary0 ;
  463. G3 = EXTR Grad0 'T,Z' Bary0 ;
  464. G0 = (G1 G2 G3) ;
  465. * Add the element to inverse
  466. SI ((N0 PSCA G0) < 0.0) ;
  467. Inv0 = Inv0 ET Elem0 ;
  468. FINS ;
  469. FIN Loop1 ;
  470. * If there are some elements to invert
  471. SI ((NBEL Inv0) > 2) ;
  472. Inv0 = Inv0 ELEM TRI3 ;
  473. Surf0 = (DIFF Surf0 Inv0) ET (INVE Inv0) ;
  474. FINS ;
  475. FINS ;
  476. * Give information about the output mesh quality
  477. VERM Surf0 ;
  478. FINS ;
  479.  
  480. * Restore the initial options
  481. OPTI 'DIME' Dime0 'MODE' Mode0 ;
  482. SI (NEG ElemOpt0 ' ') ;
  483. OPTI 'ELEM' ElemOpt0 ;
  484. FINS ;
  485.  
  486. FINP Surf0 ;
  487.  
  488.  
  489.  
  490.  

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