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