Télécharger toposurf.procedur

Retour à la liste

Numérotation des lignes :

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

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