Télécharger @isosurf.procedur

Retour à la liste

Numérotation des lignes :

  1. * @ISOSURF PROCEDUR GOUNAND 15/09/16 21:15:00 8625
  2. ************************************************************************
  3. * Procedure qui extrait les isosurfaces dont les valeurs sont
  4. * listées dans une liste de réels (LIS1) d'un champoint (HANA1)
  5. * appuyé sur un maillage (MASSIF0).
  6. *
  7. * Le résultat final est constitué du maillage surfacique
  8. * regroupant l'ensemble des isosurface MAIF1 et du champoint
  9. * CHPF1 des isovaleurs LIS1 appuyées sur MAIF1.
  10. *
  11. * Postraitement TRAC CACH MAIF1 CHPF1 ;
  12. *
  13. * Syntaxe :
  14. * ---------
  15. * MAIL1 CHPF1 = @ISOSURF MASSIF0 LIS1 HANA1 ;
  16. *
  17. * Entrée :
  18. * ---------
  19. * MASSIF0 : Maillage support du champoint
  20. *
  21. * LIS1 : Liste (LISTREEL) d'isovaleurs à rechercher
  22. *
  23. * HANA1 : Champoint appuye sur MASSIF0
  24. *
  25. * Sortie :
  26. * ---------
  27. * MAIF1 : Maillage de l'ensemble des isosurfaces
  28. *
  29. * CHPF1 : Champoint des isovaleurs LIS1 appuyées sur MAIF1
  30. *
  31. ************************************************************************
  32. * Remarque 1 : Attention la procédure utilise une élimination des points
  33. * doubles des isosurfaces extraites avec un epsilon
  34. * calculé aussi intelligemment que possible.
  35. * (Gounand: l'algo en Gibiane utilise ELIM mais l'algo en
  36. * Esope n'en a plus besoin, et fait l'élimination en
  37. * interne dans l'opérateur ISOV depuis la fiche 8152).
  38. *
  39. * Remarque 2 : pour l'algo Gibiane :
  40. * 6 isovaleurs sur 1016 tetra : cpu = 3.7 s
  41. * 6 isovaleurs sur 6950 tetra : cpu = 39 s
  42. * 6 isovaleurs sur 52000 tetra : cpu = 1315 s
  43. * ---> A reprogrammer en fortran (Gounand : c'est fait !)
  44. *************************************************************************
  45. 'DEBPROC' @ISOSURF MASSIF0*'MAILLAGE' LIS1*'LISTREEL' HANA1*'CHPOINT' ;
  46. *
  47. * Paramètres
  48. * lgibi = FAUX : on utilise l'opérateur ISOV pour construire
  49. * les isosurfaces
  50. * lgibi = VRAI : ancien algorithme en Gibiane pour mémoire
  51. lgibi = FAUX ;
  52. * Garde-fou liste vide
  53. NIS1 = DIME LIS1 ;
  54. SI ( NIS1 '<' 1 ) ;
  55. 'ERREUR' 'Liste vide' ;
  56. 'QUITTER' @ISOSURF ;
  57. FINSI ;
  58. *
  59. * Garde-fou isovaleur non comprise dans le champoint
  60. MAXC1 = 'MAXIMUM' HANA1 ;
  61. MINC1 = 'MINIMUM' HANA1 ;
  62. MAXL1 = 'MAXIMUM' LIS1 ;
  63. MINL1 = 'MINIMUM' LIS1 ;
  64. vmax = 'MAXIMUM' ('PROG' maxc1 maxl1 minc1 minl1) 'ABS' ;
  65. vref = '+' ('*' vmax 1.D-14) 1.D-100 ;
  66. SI ( (MAXL1 '>' ('+' MAXC1 vref)) 'OU' (MINL1 '<' ('-' MINC1 vref)) ) ;
  67. 'ERREUR' 'Isovaleur de la liste hors champs' ;
  68. 'QUITTER' @ISOSURF ;
  69. FINSI ;
  70. *
  71. * Changement du maillage en tétraèdres
  72. *
  73. chp = hana1 ; mail = massif0 ;
  74. mailt = 'CHANGER' mail 'TET4' ;
  75. *
  76. * S'il y a une différence de nombres de noeuds entre mail et mailt
  77. * on va faire une projection
  78. *
  79. 'SI' ('>' ('NBNO' mailt) ('NBNO' mail)) ;
  80. mailp = 'CHANGER' 'POI1' mail ;
  81. mailtp = 'CHANGER' 'POI1' mailt ;
  82. dmp = 'DIFF' mailp mailtp ;
  83. cha = 'CHANGER' 'CHAM' chp mail ;
  84. chp2 = 'PROI' dmp cha ;
  85. chp = chp '+' chp2 ;
  86. 'FINSI' ;
  87. massif0 = mailt ;
  88. hana1 = chp ;
  89. *
  90. * Garde-fou maillage non uniquement constitué de TET4
  91. *gounand Inutile normalement vu les lignes précédentes
  92. *LMOT1 = MASSIF0 ELEM 'TYPE' ;
  93. *NMT1 = DIME LMOT1 ;
  94. *SI (NMT1 > 1) ;
  95. * 'ERREUR' 'Maillage comprenant plusieurs types' ;
  96. * 'QUITTER' @ISOSURF ;
  97. *SINON ;
  98. * MOT1 = EXTR LMOT1 1 ;
  99. * SI ( NEG MOT1 'TET4' ) ;
  100. * 'ERREUR' 'Maillage non constitué de TET4' ;
  101. * 'QUITTER' @ISOSURF ;
  102. * FINSI ;
  103. *FINSI ;
  104. *
  105. LIS1 = ORDO LIS1 ;
  106. *
  107. * Calcul du paramètre d'élimination
  108. *
  109. lmm = 'PROG' 1.D-30 ;
  110. 'REPETER' idim ('VALEUR' 'DIME') ;
  111. cc = 'COORDONNEE' &idim MASSIF0 ;
  112. dmm = '-' ('MAXIMUM' cc) ('MINIMUM' cc) ;
  113. lmm = 'ET' lmm ('PROG' dmm) ;
  114. 'FIN' idim ;
  115. EPSI1 = 0.000001 '*' ('MAXIMUM' lmm) ;
  116. *
  117. 'SI' lgibi ;
  118. *
  119. * Début de l'algorithme codé en gibiane
  120. *
  121. * nombre elements du maillage
  122. NB1 = NBEL MASSIF0 ;
  123. * Boucle sur les isovaleurs
  124. REPETER BOUS1 NIS1 ;
  125. ISV1 = EXTR LIS1 &BOUS1 ;
  126. *
  127. * Nombre elements avec triangle isovaleur
  128. NBIV1 = 0 ;
  129. * Boucle sur les elements constitutifs du maillage
  130. REPETER BOUC1 NB1 ;
  131. * On extrait le nieme element
  132. TET1 = MASSIF0 ELEM &BOUC1 ;
  133. TET2 = CHANGER 'POI1' TET1 ;
  134. * On extrait les valeurs du champoint et les points
  135. P1 = TET2 POIN 1 ;
  136. P2 = TET2 POIN 2 ;
  137. P3 = TET2 POIN 3 ;
  138. P4 = TET2 POIN 4 ;
  139. VAL1 = EXTR HANA1 SCAL P1 ;
  140. VAL2 = EXTR HANA1 SCAL P2 ;
  141. VAL3 = EXTR HANA1 SCAL P3 ;
  142. VAL4 = EXTR HANA1 SCAL P4 ;
  143. LIT1 = PROG VAL1 VAL2 VAL3 VAL4 ;
  144. MAX1 = MAXI LIT1 ;
  145. MIN1 = MINI LIT1 ;
  146. * Rentrer dans element si ISV1 (isovaleur) comprise
  147. SI ((ISV1 > MAX1) OU (ISV1 < MIN1)) ;
  148. * Pas d'isovaleur dans cet element
  149. STRI1 = FAUX ;
  150. SINON ;
  151. * Il y a une isovaleur dans l'element
  152. EXP0 = FAUX ;
  153. STRI1 = FAUX ;
  154. * On ordonne les valeurs aux points (et les points)
  155. REPETER BOUC2 ;
  156. SI (VAL1 &lt;EG VAL2) ;
  157. SI (VAL2 &lt;EG VAL3) ;
  158. SI (VAL3 &lt;EG VAL4) ;
  159. QUITTER BOUC2 ;
  160. SINON ;
  161. VAX1 = VAL4 ;
  162. PAX1 = P4 ;
  163. VAL4 = VAL3 ;
  164. P4 = P3 ;
  165. VAL3 = VAX1 ;
  166. P3 = PAX1 ;
  167. FINSI ;
  168. SINON ;
  169. VAX1 = VAL3 ;
  170. PAX1 = P3 ;
  171. VAL3 = VAL2 ;
  172. P3 = P2 ;
  173. VAL2 = VAX1 ;
  174. P2 = PAX1 ;
  175. FINSI ;
  176. SINON ;
  177. VAX1 = VAL2 ;
  178. PAX1 = P2 ;
  179. VAL2 = VAL1 ;
  180. P2 = P1 ;
  181. VAL1 = VAX1 ;
  182. P1 = PAX1 ;
  183. FINSI ;
  184. FIN BOUC2 ;
  185. * On a fini d'ordonner les valeurs
  186. * On teste si l'isovaleur correspond à un
  187. * noeud du tetrahèdre
  188. NPEQ1 = 0 ;
  189. SI (ISV1 EGA VAL1) ;
  190. NPEQ1 = NPEQ1 + 1 ;
  191. FINSI ;
  192. SI (ISV1 EGA VAL2) ;
  193. NPEQ1 = NPEQ1 + 1 ;
  194. FINSI ;
  195. SI (ISV1 EGA VAL3) ;
  196. NPEQ1 = NPEQ1 + 1 ;
  197. FINSI ;
  198. *gounand SI (ISV1 EGA VAL3) ;
  199. SI (ISV1 EGA VAL4) ;
  200. NPEQ1 = NPEQ1 + 1 ;
  201. FINSI ;
  202. * NPEQ1 indique le nombre de noeuds du tetrahèdre
  203. * qui correspondent à l'isovaleur
  204. SI (NPEQ1 EGA 4) ;
  205. * Isosurface = 4 faces du tetrahedre
  206. NBIV1 = NBIV1 + 1 ;
  207. STRI1 = VRAI ;
  208. MAIL1 = CHANGER 'FACE' TET1 ;
  209. SINON ;
  210. SI (NPEQ1 EGA 3) ;
  211. * Isosurface = 1 face du tetrahedre
  212. STRI1 = VRAI ;
  213. NBIV1 = NBIV1 + 1 ;
  214. SI (ISV1 EGA VAL1) ;
  215. * Point P1
  216. PIN1 = P1 ;
  217. SINON ;
  218. * Point P4 ;
  219. PIN1 = P4 ;
  220. FINSI ;
  221. * Toujours P3 et P2 inclus
  222. PIN2 = P2 ;
  223. PIN3 = P3 ;
  224. SINON ;
  225. SI (NPEQ1 EGA 2) ;
  226. * Isosurface = 1 segment ou 1 triangle
  227. * appuye sur 2 points du tetrahedre
  228. SI ((ISV1 EGA VAL1) OU (ISV1 EGA VAL4)) ;
  229. * Isosurface = 1 segment (on ne fait rien)
  230. SINON ;
  231. * Isosurface = 1 triangle appuye sur 2 points tetra
  232. NBIV1 = NBIV1 + 1 ;
  233. STRI1 = VRAI ;
  234. * point intersection (P1,P4)
  235. PIN1 = P2 ;
  236. PIN2 = P3 ;
  237. VEC1 = P4 MOINS P1 ;
  238. DVE1 = (ISV1 - VAL1) / (VAL1 - VAL4) ;
  239. VEC2 = DVE1 * VEC1 ;
  240. PIN3 = P1 MOINS VEC2 ;
  241. FINSI ;
  242. SINON ;
  243. SI (NPEQ1 EGA 1) ;
  244. * Isosurface = 1 point ou 1 triangle
  245. * appuye sur 1 point du tetrahedre
  246. SI ((ISV1 EGA VAL1) OU (ISV1 EGA VAL4)) ;
  247. * Isosurface = 1 point (on ne fait rien)
  248. SINON ;
  249. * Isosurface = 1 triangle appuye sur 1 point
  250. * du tetrahedre
  251. NBIV1 = NBIV1 + 1 ;
  252. STRI1 = VRAI ;
  253. SI (ISV1 EGA VAL2) ;
  254. * Le point appuye est P2
  255. * les autres points intersection (P1,P3)
  256. PIN1 = P2 ;
  257. VEC1 = P3 MOINS P1 ;
  258. DVE1 = (ISV1 - VAL1) / (VAL1 - VAL3) ;
  259. VEC2 = DVE1 * VEC1 ;
  260. PIN2 = P1 MOINS VEC2 ;
  261. SINON ;
  262. * Le point appuye est P3
  263. * les autres points intersection (P2,P4)
  264. PIN1 = P3 ;
  265. VEC1 = P4 MOINS P2 ;
  266. DVE1 = (ISV1 - VAL2) / (VAL2 - VAL4) ;
  267. VEC2 = DVE1 * VEC1 ;
  268. PIN2 = P2 MOINS VEC2 ;
  269. FINSI ;
  270. * Il y a toujours un point intersection (P1,P4)
  271. VEC1 = P4 MOINS P1 ;
  272. DVE1 = (ISV1 - VAL1) / (VAL1 - VAL4) ;
  273. VEC2 = DVE1 * VEC1 ;
  274. PIN3 = P1 MOINS VEC2 ;
  275. FINSI ;
  276. SINON ;
  277. * Isosurface = 1 triangle quelconque section
  278. * du tetrahedre ou un quadrangle (= 2 triangles)
  279. NBIV1 = NBIV1 + 1 ;
  280. STRI1 = VRAI ;
  281. SI (ISV1 < VAL2) ;
  282. * Isosurface entre V1 et V2
  283. * Points intersection (P1,P2) (P1,P3)
  284. * Un seul triangle
  285. VEC1 = P2 MOINS P1 ;
  286. DVE1 = (ISV1 - VAL1) / (VAL1 - VAL2) ;
  287. VEC2 = DVE1 * VEC1 ;
  288. PIN1 = P1 MOINS VEC2 ;
  289. VEC1 = P3 MOINS P1 ;
  290. DVE1 = (ISV1 - VAL1) / (VAL1 - VAL3) ;
  291. VEC2 = DVE1 * VEC1 ;
  292. PIN2 = P1 MOINS VEC2 ;
  293. SINON ;
  294. SI (ISV1 > VAL3) ;
  295. * Isosurface entre V3 et V4
  296. * Points intersection (P3,P4) (P2,P4)
  297. * Un seul triangle
  298. VEC1 = P4 MOINS P3 ;
  299. DVE1 = (ISV1 - VAL3) / (VAL3 - VAL4) ;
  300. VEC2 = DVE1 * VEC1 ;
  301. PIN1 = P3 MOINS VEC2 ;
  302. VEC1 = P4 MOINS P2 ;
  303. DVE1 = (ISV1 - VAL2) / (VAL2 - VAL4) ;
  304. VEC2 = DVE1 * VEC1 ;
  305. PIN2 = P2 MOINS VEC2 ;
  306. SINON ;
  307. * Isosurface entre V2 et V3 ;
  308. * Points intersection (P1,P3) (P2,P3) (P2,P4)
  309. * Un quadrangle = 2 triangles
  310. EXP0 = VRAI ;
  311. VEC1 = P3 MOINS P1 ;
  312. DVE1 = (ISV1 - VAL1) / (VAL1 - VAL3) ;
  313. VEC2 = DVE1 * VEC1 ;
  314. PIN0 = P1 MOINS VEC2 ;
  315. VEC1 = P3 MOINS P2 ;
  316. DVE1 = (ISV1 - VAL2) / (VAL2 - VAL3) ;
  317. VEC2 = DVE1 * VEC1 ;
  318. PIN1 = P2 MOINS VEC2 ;
  319. VEC1 = P4 MOINS P2 ;
  320. DVE1 = (ISV1 - VAL2) / (VAL2 - VAL4) ;
  321. VEC2 = DVE1 * VEC1 ;
  322. PIN2 = P2 MOINS VEC2 ;
  323. FINSI ;
  324. FINSI ;
  325. * Il y toujours un point intersection
  326. * entre P1 et P4
  327. VEC1 = P4 MOINS P1 ;
  328. DVE1 = (ISV1 - VAL1) / (VAL1 - VAL4) ;
  329. VEC2 = DVE1 * VEC1 ;
  330. PIN3 = P1 MOINS VEC2 ;
  331. FINSI ;
  332. FINSI ;
  333. FINSI ;
  334. FINSI ;
  335. FINSI ;
  336. *
  337. * Fabrication de la surface si elle existe
  338. SI STRI1 ;
  339. * Premier triangle
  340. LS1 = DROIT 1 PIN1 PIN2 ;
  341. LS2 = DROIT 1 PIN2 PIN3 ;
  342. LS3 = DROIT 1 PIN3 PIN1 ;
  343. CONT1 = LS1 ET LS2 ET LS3 ;
  344. MAIL1 = COUL VERT (SURF CONT1 'PLANE') ;
  345. SI EXP0 ;
  346. * Second triangle
  347. LS1 = DROIT 1 PIN3 PIN0 ;
  348. LS2 = DROIT 1 PIN0 PIN1 ;
  349. LS3 = DROIT 1 PIN1 PIN3 ;
  350. CONT1 = LS1 ET LS2 ET LS3 ;
  351. MAIL2 = COUL VERT (SURF CONT1 'PLANE') ;
  352. MAIL1 = MAIL1 ET MAIL2 ;
  353. FINSI ;
  354. SI (NBIV1 EGA 1) ;
  355. MAIT1 = MAIL1 ;
  356. SINON ;
  357. MAIT1 = MAIT1 ET MAIL1 ;
  358. FINSI ;
  359. FINSI ;
  360. FIN BOUC1 ;
  361. * Elimination points multiples du maillage de surface
  362. * obtenu par somme de triangles independants
  363. ELIM EPSI1 MAIT1 ;
  364. CHP1 = MANU 'CHPO' MAIT1 1 'SCAL' ISV1 'NATURE' 'DISCRET' ;
  365. SI (NIS1 > 1) ;
  366. SI (&BOUS1 EGA 1) ;
  367. MAIF1 = MAIT1 ;
  368. CHPF1 = CHP1 ;
  369. SINON ;
  370. MAIF1 = MAIF1 ET MAIT1 ;
  371. CHPF1 = CHPF1 ET CHP1 ;
  372. FINSI ;
  373. SINON ;
  374. MAIF1 = MAIT1 ;
  375. CHPF1 = CHP1 ;
  376. FINSI ;
  377. FIN BOUS1 ;
  378. *
  379. * Fin de l'algorithme en gibiane
  380. * et début de la version fortran
  381. 'SINON' ;
  382. hana2 = 'CHANGER' 'CHAM' hana1 massif0 ;
  383. lmail = vrai ;
  384. 'REPETER' BOUS1 NIS1 ;
  385. ISV1 = EXTR LIS1 &BOUS1 ;
  386. *dbg 'MESSAGE' ('CHAINE' 'isv1=' isv1) ;
  387. MAIT1 = 'ISOV' hana2 ISV1 ;
  388. *dbg ymin = 'MINIMUM' ('COORDONNEE' 2 MAIT1) ;
  389. *dbg 'MESSAGE' ('CHAINE' 'ymin=' ymin) ;
  390. 'SI' ('>' ('NBEL' mait1) 0) ;
  391. *sg: normalement inutile depuis fiche 8152 'ELIM' EPSI1 MAIT1 ;
  392. CHP1 = 'MANU' 'CHPO' MAIT1 1 'SCAL' ISV1 'NATURE' 'DISCRET' ;
  393. 'SI' lmail ;
  394. MAIF1 = MAIT1 ;
  395. CHPF1 = CHP1 ;
  396. lmail = FAUX ;
  397. 'SINON' ;
  398. MAIF1 = MAIF1 ET MAIT1 ;
  399. CHPF1 = CHPF1 ET CHP1 ;
  400. 'FINSI' ;
  401. 'FINSI' ;
  402. 'FIN' BOUS1 ;
  403. 'FINSI' ;
  404. *
  405. 'FINPROC' MAIF1 CHPF1 ;
  406. *
  407. ********************************************************************
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  

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