Télécharger toposurf.procedur

Retour à la liste

Numérotation des lignes :

  1. * TOPOSURF PROCEDUR CB215821 17/12/01 21:15:14 9641
  2. ************************************************************************
  3. ** Extraction procedure of a smoothed surface from a topology
  4. ** Version 1.0
  5. ** Guenhael Le Quilliec (LMR - Polytech Tours)
  6. ** 2016/12/13
  7. ************************************************************************
  8.  
  9. DEBP TOPOSURF Topo0*'TABLE' Filter0/'ENTIER' Iso0/'FLOTTANT'
  10. Elim0/'FLOTTANT' Dec0/'FLOTTANT' Epai0/'FLOTTANT'
  11. Orie0/'LOGIQUE' ;
  12.  
  13. * Default value of the filtering rate to apply
  14. SI (NON (EXIS Filter0)) ;
  15. Filter0 = 1 ;
  16. FINS ;
  17.  
  18. * Default value of the Iso surface to extract
  19. SI (NON (EXIS Iso0)) ;
  20. Iso0 = 0.5 ;
  21. FINS ;
  22.  
  23. * Default minimum size of the output surface elements
  24. SI (NON (EXIS Elim0)) ;
  25. Elim0 = 1.0e-6 ;
  26. FINS ;
  27.  
  28. * Default value of the extra thickness
  29. * This value should be very low compared to the element size
  30. * when the domain is not convex.
  31. * If not, Enve1 may be wrong and GRAD won't work
  32. SI (NON (EXIS Dec0)) ;
  33. Dec0 = 1.0e-3 ;
  34. FINS ;
  35.  
  36. * Default value of the thickness of the output mesh when the input
  37. * modèle is 2D
  38. SI (NON (EXIS Epai0)) ;
  39. Epai0 = 1.0 ;
  40. FINS ;
  41.  
  42. * Orientation of the output surface will be applied by default.
  43. SI (NON (EXIS Orie0)) ;
  44. Orie0 = VRAI ;
  45. FINS ;
  46.  
  47. * Model
  48. Mod0 = Topo0.'MODELE' ;
  49. * Mesh
  50. Mail0 = EXTR Mod0 'MAIL' ;
  51. * Topology
  52. TopoCh0 = Topo0.'TOPOLOGIE_CH' ;
  53. * One of its types of elements
  54. Type0 = EXTR (Mail0 ELEM 'TYPE') 1 ;
  55.  
  56. * Filtering the topology field
  57. SI (Filter0 > 0) ;
  58. REPE loop1 Filter0 ;
  59. TopoCh0 = CHAN 'CHAM' (CHAN 'CHPO' Mod0 TopoCh0 'MOYE')
  60. Mod0 'GRAVITE' ;
  61. FIN loop1 ;
  62. FINS ;
  63. * Node topology field
  64. Ch0 = CHAN 'CHPO' TopoCh0 Mod0 'MOYE' ;
  65. Ch0 = Ch0 * (1.0 / (MAXI Ch0)) ;
  66.  
  67. * Save the options 'DIME', 'MODE', 'ELEM'
  68. Dime0 = VALE 'DIME' ;
  69. Mode0 = VALE 'MODE' ;
  70. ElemOpt0 = VALE 'ELEM' ;
  71.  
  72. * 2D Topology
  73. SI (EXIS (MOTS 'TRI3' 'TRI6' 'QUA4' 'QUA8') Type0) ;
  74. * Switch to 2D
  75. OPTI 'DIME' 2 'MODE' 'PLAN' 'CONT' ;
  76. * If the reasult has to stay in 2D
  77. SI ((ABS Epai0) < 1.0e-9) ;
  78. Type3D = FAUX ;
  79. * Change option 'ELEM' to TRI3
  80. OPTI 'ELEM' TRI3 ;
  81. * Change the element type of the mesh
  82. Mail0 = CHAN 'TRI3' Mail0 ;
  83. * Domain boundary
  84. Cont0 = CONT Mail0 ;
  85. * Shifted boundary
  86. Norm0 = PRES 'MASS' Mod0 (-1.0 * Dec0) Cont0 ;
  87. Norm0 = CHAN 'COMP' (MOTS 'FX' 'FY') (MOTS 'UX' 'UY')
  88. Norm0 ;
  89. Cont1 = Cont0 PLUS Norm0 ;
  90. * External layer of elements
  91. Layer0 = Cont0 REGL 1 Cont1 ;
  92. * Add a field for the shifted boundary TODO 20
  93. Ch0 = Ch0 ET
  94. (MANU 'CHPO' Cont1 1 'SCAL' -1.0e3 'NATURE' 'DIFFUS') ;
  95. * New model
  96. Mod1 = MODE (Mail0 ET Layer0) 'MECANIQUE' 'ELASTIQUE' ;
  97. * Element field
  98. MCh0 = CHAN 'CHAM' Ch0 Mod1 'NOEUD' ;
  99. * Extract the iso surface
  100. Surf0 = (ISOV MCh0 Iso0) ELEM SEG2 ;
  101. * Fuse the nodes
  102. ELIM Surf0 Elim0 ;
  103. * If the surface has to be oriented
  104. SI Orie0 ;
  105. * Fuse any superposed element
  106. Surf1 = TABL ;
  107. REPE loop1 (NBEL Surf0) ;
  108. * Considered element
  109. Elem0 = Surf0 ELEM &loop1 ;
  110. P0 = CHAN POI1 Elem0 ;
  111. SI (EGA (NBEL P0) 2) ;
  112. N1 = NOEU (P0 POIN 1) ;
  113. N2 = NOEU (P0 POIN 2) ;
  114. SI (N1 > N2) ;
  115. N0 = N2 ;
  116. N2 = N1 ;
  117. N1 = N0 ;
  118. FINS ;
  119. Surf1.(CHAI N1 '-' N2) = Elem0 ;
  120. FINS ;
  121. FIN loop1 ;
  122. Ind0 = INDE Surf1 ;
  123. Surf0 = Surf1.(Ind0.(DIME Ind0)) ;
  124. REPE loop1 ((DIME Ind0) - 1) ;
  125. Surf0 = Surf0 ET Surf1.(Ind0.&loop1) ;
  126. FIN loop1 ;
  127. * Create the gradient field giving the orientation
  128. Ch0 = (REDU Ch0 Mail0) ET (MANU 'CHPO' Cont1 1 'SCAL'
  129. -1.0 'NATURE' 'DIFFUS') ;
  130. Mod1 = MODE (Mail0 ET Layer0) 'THERMIQUE' ;
  131. Grad0 = -1.0 * (GRAD Mod1 (CHAN 'COMP' 'T' Ch0)) ;
  132. Grad0 = CHAN 'NOEUD' Mod1 Grad0 ;
  133. * Centers of mass
  134. REPE loop1 (NBEL Surf0) ;
  135. Elem0 = Surf0 ELEM &loop1 ;
  136. SI (&loop1 NEG 1) ;
  137. Pts0 = Pts0 ET (BARY Elem0) ;
  138. SINO ;
  139. Pts0 = BARY Elem0 ;
  140. FINS ;
  141. FIN loop1 ;
  142. * Projection of the gradient on the centers of mass
  143. Grad0 = PROI Pts0 Grad0 ;
  144. * Initialization of the loop on each elements
  145. Inv0 = (0.0 0.0) ET (0.0 0.0) ;
  146. * Loop on each element of the surface
  147. REPE loop1 (NBEL Surf0) ;
  148. * Considered element
  149. Elem0 = Surf0 ELEM &loop1 ;
  150. * Compute the normal direction of the element
  151. P0 = CHAN POI1 Elem0 ;
  152. P1 = P0 POIN 1 ;
  153. P2 = P0 POIN 2 ;
  154. N0 = (P2 MOIN P1) TOUR 90.0 (0.0 0.0) ;
  155. * Center of mass
  156. Bary0 = Pts0 POIN &loop1 ;
  157. * Desired orientation
  158. G1 = EXTR Grad0 'T,X' Bary0 ;
  159. G2 = EXTR Grad0 'T,Y' Bary0 ;
  160. G0 = (G1 G2) ;
  161. * Add the element to inverse
  162. SI ((N0 PSCA G0) < 0.0) ;
  163. Inv0 = Inv0 ET Elem0 ;
  164. FINS ;
  165. FIN loop1 ;
  166. * If there are some elements to invert
  167. SI ((NBEL Inv0) > 2) ;
  168. Inv0 = Inv0 ELEM SEG2 ;
  169. Surf0 = (DIFF Surf0 Inv0) ET (INVE Inv0) ;
  170. FINS ;
  171. FINS ;
  172. * Give information about the output mesh quality
  173. VERM Surf0 ;
  174. SINO ;
  175. Type3D = VRAI ;
  176. * Change option 'ELEM' to QUA4 or QUA8
  177. SI (EXIS (MOTS 'TRI3' 'QUA4') Type0) ;
  178. OPTI 'ELEM' QUA4 ;
  179. SINO ;
  180. OPTI 'ELEM' QUA8 ;
  181. FINS ;
  182. * Force the option 'DIME' to 3D
  183. Dime0 = 3 ;
  184. * Switch to 3D
  185. OPTI 'DIME' 3 ;
  186. * Change option 'ELEM' to CUB8 or CU20
  187. SI (EXIS (MOTS 'TRI3' 'QUA4') Type0) ;
  188. OPTI 'ELEM' CUB8 ;
  189. SINO ;
  190. OPTI 'ELEM' CU20 ;
  191. FINS ;
  192. * Create the second side of the 3D mesh from the 2D mesh
  193. Recto0 ChR0 = Mail0 Ch0 PLUS (0.0 0.0 Epai0) ;
  194. * Create the 3D mesh
  195. Mail0 = Mail0 VOLU 1 Recto0 ;
  196. * Create the 3D topology field
  197. Ch0 = Ch0 ET ChR0 ;
  198. * Create the 3D model
  199. Mod0 = MODE Mail0 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' ;
  200. FINS ;
  201. SINO ;
  202. Type3D = VRAI ;
  203. FINS ;
  204.  
  205. * 3D Topology
  206. SI Type3D ;
  207. * Change option 'ELEM' to CUB8 or CU20
  208. SI (EXIS (MOTS 'CUB8' 'PRI6' 'PYR5' 'TET4') Type0) ;
  209. OPTI 'ELEM' CUB8 ;
  210. SINO ;
  211. OPTI 'ELEM' CU20 ;
  212. FINS ;
  213. * Domain boundary
  214. Enve0 = ENVE Mail0 ;
  215. * Normal direction on the domain boundary
  216. Norm0 = PRES 'MASS' Mod0 (-1.0 * Dec0) Enve0 ;
  217. Norm0 = CHAN 'COMP' (MOTS 'FX' 'FY' 'FZ')
  218. (MOTS 'UX' 'UY' 'UZ') Norm0 ;
  219. * Partitionate the surfaces
  220. NbPart0 = 32 ;
  221. SI (NbPart0 > (NBEL Enve0)) ;
  222. NbPart0 = 16 ;
  223. FINS ;
  224. SI (NbPart0 > (NBEL Enve0)) ;
  225. NbPart0 = 8 ;
  226. FINS ;
  227. PEnve0 = PART 'ARLE 'NbPart0 Enve0 ;
  228. * External layer of elements
  229. Layer0 = TABL 'ESCLAVE' ;
  230. * Loop on each entry of PEnve0
  231. * (since PART does not always respect NbPart0)
  232. REPE loop1 (DIME PEnve0) ;
  233. SI (EXIS PEnve0 &loop1) ;
  234. PEnve1 = PEnve0.&loop1 PLUS
  235. (REDU Norm0 PEnve0.&loop1) ;
  236. Layer0.&loop1 = PEnve0.&loop1 VOLU 1 PEnve1 ;
  237. SI (&Loop1 > 1) ;
  238. Enve1 = Enve1 ET PEnve1 ;
  239. ELIM Enve1 PEnve1 Elim0 ;
  240. SINO ;
  241. Enve1 = PEnve1 ;
  242. FINS ;
  243. FINS ;
  244. FIN loop1 ;
  245. Layer0 = ET Layer0 ;
  246. * Add a field for the shifted boundary
  247. Ch0 = Ch0 ET
  248. (MANU 'CHPO' Enve1 1 'SCAL' -1.0e3 'NATURE' 'DIFFUS') ;
  249. * Create a mesh that will be compatible with ISOV
  250. IsoMail0 = Mail0 ET Layer0 ;
  251. Type0 = IsoMail0 ELEM 'TYPE' ;
  252. IsoMail1 = TABL 'ESCLAVE' ;
  253. REPE loop1 (DIME Type0) ;
  254. Type1 = EXTR Type0 &loop1 ;
  255. SI ((EGA Type1 'TET4') OU (EGA Type1 'CUB8')) ;
  256. IsoMail1.&loop1 = IsoMail0 ELEM Type1 ;
  257. SINO ;
  258. IsoMail1.&loop1 = CHAN 'TET4' (IsoMail0 ELEM Type1) ;
  259. FINS ;
  260. FIN loop1 ;
  261. IsoMail0 = ET IsoMail1 ;
  262. * New model
  263. Mod1 = MODE IsoMail0 'MECANIQUE' 'ELASTIQUE' ;
  264. * Element field
  265. MCh0 = CHAN 'CHAM' Ch0 Mod1 'NOEUD' ;
  266. * Extract the iso surface
  267. Surf0 = ISOV MCh0 Iso0 ;
  268. * Change the surface to TRI3
  269. Type0 = Surf0 ELEM 'TYPE' ;
  270. Surf1 = TABL 'ESCLAVE' ;
  271. REPE loop1 (DIME Type0) ;
  272. Type1 = EXTR Type0 &loop1 ;
  273. SI (EGA Type1 'TRI3') ;
  274. Surf1.&loop1 = Surf0 ELEM Type1 ;
  275. SINO ;
  276. SI (NEG Type1 'POI1') ;
  277. Surf1.&loop1 = CHAN 'TRI3' (Surf0 ELEM Type1) ;
  278. FINS ;
  279. FINS ;
  280. FIN loop1 ;
  281. Surf0 = ET Surf1 ;
  282. * Fuse the nodes
  283. ELIM Surf0 Elim0 ;
  284. * If the surface has to be oriented
  285. SI Orie0 ;
  286. * Fuse any superposed element
  287. Surf1 = TABL ;
  288. REPE loop1 (NBEL Surf0) ;
  289. * Considered element
  290. Elem0 = Surf0 ELEM &loop1 ;
  291. P0 = CHAN POI1 Elem0 ;
  292. SI (EGA (NBEL P0) 3) ;
  293. N1 = NOEU (P0 POIN 1) ;
  294. N2 = NOEU (P0 POIN 2) ;
  295. N3 = NOEU (P0 POIN 3) ;
  296. N0 = ORDO (LECT N1 N2 N3) 'CROI' ;
  297. Surf1.(CHAI (EXTR N0 1) '-' (EXTR N0 2)
  298. '-' (EXTR N0 3)) = Elem0 ;
  299. FINS ;
  300. FIN loop1 ;
  301. Ind0 = INDE Surf1 ;
  302. Surf0 = Surf1.(Ind0.(DIME Ind0)) ;
  303. REPE loop1 ((DIME Ind0) - 1) ;
  304. Surf0 = Surf0 ET Surf1.(Ind0.&loop1) ;
  305. FIN loop1 ;
  306. * Create the gradient field giving the orientation
  307. Ch0 = (REDU Ch0 Mail0) ET
  308. (MANU 'CHPO' Enve1 1 'SCAL' -1.0 'NATURE' 'DIFFUS') ;
  309. Mod1 = MODE (Mail0 ET Layer0) 'THERMIQUE' ;
  310. Grad0 = -1.0 * (GRAD Mod1 (CHAN 'COMP' 'T' Ch0)) ;
  311. Grad0 = CHAN 'NOEUD' Mod1 Grad0 ;
  312. * Centers of mass
  313. REPE loop1 (NBEL Surf0) ;
  314. Elem0 = Surf0 ELEM &loop1 ;
  315. SI (&loop1 NEG 1) ;
  316. Pts0 = Pts0 ET (BARY Elem0) ;
  317. SINO ;
  318. Pts0 = BARY Elem0 ;
  319. FINS ;
  320. FIN loop1 ;
  321. * Projection of the gradient on the centers of mass
  322. Grad0 = PROI Pts0 Grad0 ;
  323. * Initialization of the loop on each elements
  324. Inv0 = (0.0 0.0 0.0) ET (0.0 0.0 0.0) ;
  325. * Loop on each element of the surface
  326. REPE loop1 (NBEL Surf0) ;
  327. * Considered element
  328. Elem0 = Surf0 ELEM &loop1 ;
  329. * Compute the normal direction of the element
  330. P0 = CHAN POI1 Elem0 ;
  331. P1 = P0 POIN 1 ;
  332. P2 = P0 POIN 2 ;
  333. P3 = P0 POIN 3 ;
  334. N0 = (P2 MOIN P1) PVEC (P3 MOIN P2) ;
  335. * Center of mass
  336. Bary0 = Pts0 POIN &loop1 ;
  337. * Desired orientation
  338. G1 = EXTR Grad0 'T,X' Bary0 ;
  339. G2 = EXTR Grad0 'T,Y' Bary0 ;
  340. G3 = EXTR Grad0 'T,Z' Bary0 ;
  341. G0 = (G1 G2 G3) ;
  342. * Add the element to inverse
  343. SI ((N0 PSCA G0) < 0.0) ;
  344. Inv0 = Inv0 ET Elem0 ;
  345. FINS ;
  346. FIN loop1 ;
  347. * If there are some elements to invert
  348. SI ((NBEL Inv0) > 2) ;
  349. Inv0 = Inv0 ELEM TRI3 ;
  350. Surf0 = (DIFF Surf0 Inv0) ET (INVE Inv0) ;
  351. FINS ;
  352. FINS ;
  353. * Give information about the output mesh quality
  354. VERM Surf0 ;
  355. FINS ;
  356.  
  357. * Restore the initial options
  358. OPTI 'DIME' Dime0 'MODE' Mode0 'ELEM' ElemOpt0 ;
  359.  
  360. FINP Surf0 ;
  361.  
  362.  

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