1 : C ISOVA2 SOURCE BP208322 11/12/20 21:15:02 7231 2 : SUBROUTINE ISOVA2(MELEME,MCHAML,XISO,XTOL,ELEMS) 3 : IMPLICIT REAL*8 (A-H,O-Z) 4 : IMPLICIT INTEGER (I-N) 5 : C*********************************************************************** 6 : C NOM : ISOVA2 7 : C DESCRIPTION : Construit le maillage d'une isovaleur d'un champ par 8 : C éléments 9 : C 10 : C 11 : C 12 : C LANGAGE : ESOPE 13 : C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF) 14 : C mél : gounand@semt2.smts.cea.fr 15 : C*********************************************************************** 16 : C VERSION : v1, 17/12/2008, version initiale 17 : C HISTORIQUE : v1, 17/12/2008, création 18 : C HISTORIQUE : 19/12/2011, bp : ajout cub8 19 : C HISTORIQUE : 20 : C HISTORIQUE : 21 : C*********************************************************************** 22 : -INC CCOPTIO 23 : -INC CCGEOME 24 : -INC SMCOORD 25 : -INC SMELEME 26 : -INC SMCHAML 27 : -INC SMLENTI 28 : * 29 : * Segments ajustables contenant les éléments créés ainsi que leur couleur 30 : * ELEM(1) contient des POI1 31 : * ELEM(2) contient des SEG2 32 : * ELEM(3) contient des TRI3 33 : * ELEM(4) contient des TET4 34 : * ELEM(5) contient des QUA4 35 : * 36 : PARAMETER (NTYEL=5) 37 : SEGMENT ELEMS 38 : POINTEUR ELEM(NTYEL).MLENTI 39 : ENDSEGMENT 40 : * 41 : PARAMETER (NMAXNO=8) 42 : REAL*8 VAL(NMAXNO),VALP(NMAXNO) 43 : LOGICAL LVAL(NMAXNO),LTRI 44 : INTEGER IPERM(NMAXNO),NUMP(NMAXNO) 45 : CHARACTER*8 MCHA 46 : INTEGER IMPR,IRET 47 : 48 : * tableaux descriptifs des cas et des aretes 49 : PARAMETER (NMAXPT=16,NMAXCA=256,npt2=2,NMAXAR=12) 50 : INTEGER MATCAS(NMAXPT,NMAXCA),MATARE(npt2,NMAXAR) 51 : 52 : * 53 : * Fonctions locales 54 : * 55 : LOGICAL LEGISO 56 : LEGISO(X)=(ABS(X-XISO).LT.XTOL) 57 : * 58 : * Executable statements 59 : * 60 : IMPR=0 61 : IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans isova2.eso' 62 : * 63 : C*********************************************************************** 64 : c Boucle sur les zones elementaires 65 : C*********************************************************************** 66 : c 67 : NNO=NUM(/1) 68 : NEL=NUM(/2) 69 : N2=IELVAL(/1) 70 : DO I2=1,N2 71 : 72 : MELVAL=IELVAL(I2) 73 : SEGACT MELVAL 74 : 75 : * Vérification de la cohérence entre la champ et la géométrie 76 : N1PTEL=VELCHE(/1) 77 : N1EL=VELCHE(/2) 78 : IF ((N1PTEL.NE.1.AND.N1PTEL.NE.NNO).OR. 79 : $ (N1EL.NE.1.AND.N1EL.NE.NEL)) THEN 80 : * 148 2 81 : *Le champ par élément et l'objet géométrique correspondant n'ont pas le 82 : * 148 2 83 : *même nombre de points par élément : contacter votre support 83 : CALL ERREUR(148) 84 : RETURN 85 : ENDIF 86 : C 87 : c******* Aiguillage selon type d element ******************************* 88 : c 89 : IF(ITYPEL.eq.1) GOTO 101 90 : IF(ITYPEL.eq.2) GOTO 102 91 : IF(ITYPEL.eq.4) GOTO 104 92 : if(ITYPEL.eq.14) goto 114 93 : if(ITYPEL.eq.23) goto 123 94 : 95 : ******** Cas non traité ********************************************* 96 : MOTERR(1:4)=NOMS(ITYPEL) 97 : * 44 2 98 : *Type d'élément inconnu %m1:4 98 : CALL ERREUR(44) 99 : RETURN 100 : 101 : ******** Cas des POI1 *********************************************** 102 : 101 CONTINUE 103 : DO IEL=1,NEL 104 : VAL(1)=VELCHE(1,MIN(N1EL,IEL)) 105 : IF (LEGISO(VAL(1))) THEN 106 : MLENTI=ELEM(1) 107 : LECT(**)=NUM(1,IEL) 108 : LECT(**)=ICOLOR(IEL) 109 : ENDIF 110 : ENDDO 111 : goto 900 112 : 113 : ******** Cas des SEG2 *********************************************** 114 : 102 CONTINUE 115 : DO IEL=1,NEL 116 : VAL(1)=VELCHE(1,MIN(N1EL,IEL)) 117 : VAL(2)=VELCHE(MIN(N1PTEL,2),MIN(N1EL,IEL)) 118 : VMIN=MIN(VAL(1),VAL(2)) 119 : VMAX=MAX(VAL(1),VAL(2)) 120 : * Y a-t-il une isovaleur ? 121 : IF (((VMIN-XISO).LE.XTOL).AND.((VMAX-XISO).GE.-XTOL)) 122 : $ THEN 123 : NPEQ1=0 124 : DO IVAL=1,2 125 : IF (LEGISO(VAL(IVAL))) THEN 126 : NPEQ1=NPEQ1+1 127 : ENDIF 128 : ENDDO 129 : IF (NPEQ1.EQ.2) THEN 130 : * Le segment entier 131 : MLENTI=ELEM(2) 132 : LECT(**)=NUM(1,IEL) 133 : LECT(**)=NUM(2,IEL) 134 : LECT(**)=ICOLOR(IEL) 135 : ELSEIF (NPEQ1.EQ.1) THEN 136 : * Une des extrémités 137 : MLENTI=ELEM(1) 138 : IF (LEGISO(VAL(1))) THEN 139 : * Le premier point 140 : LECT(**)=NUM(1,IEL) 141 : ELSE 142 : * Le deuxième point 143 : LECT(**)=NUM(2,IEL) 144 : ENDIF 145 : LECT(**)=ICOLOR(IEL) 146 : ELSE 147 : * Un point milieu 148 : MLENTI=ELEM(1) 149 : VAL1=VAL(1) 150 : VAL2=VAL(2) 151 : NUM1=NUM(1,IEL) 152 : NUM2=NUM(2,IEL) 152 : CALL ISOVA3(XISO,VAL1,VAL2,NUM1,NUM2,MLENTI) 153 : LECT(**)=ICOLOR(IEL) 154 : ENDIF 155 : ENDIF 156 : ENDDO 157 : goto 900 158 : 159 : ******** Cas des TRI3 *********************************************** 160 : 104 CONTINUE 161 : DO IEL=1,NEL 162 : VAL(1)=VELCHE(1,MIN(N1EL,IEL)) 163 : VAL(2)=VELCHE(MIN(N1PTEL,2),MIN(N1EL,IEL)) 164 : VAL(3)=VELCHE(MIN(N1PTEL,3),MIN(N1EL,IEL)) 165 : VMIN=MIN(VAL(1),VAL(2),VAL(3)) 166 : VMAX=MAX(VAL(1),VAL(2),VAL(3)) 167 : * Y a-t-il une isovaleur ? 168 : IF (((VMIN-XISO).LE.XTOL).AND.((VMAX-XISO).GE.-XTOL)) 169 : $ THEN 170 : * Permutation permettant d'avoir VAL en ordre croissant 171 : * VAL(IPERM(1)).LE.VAL(IPERM(2)).LE.VAL(IPERM(3)) 172 : DO IVAL=1,3 173 : IPERM(IVAL)=IVAL 174 : ENDDO 174 : CALL ISHELR(3,IPERM,3,VAL,IMPR,IRET) 175 : NPEQ1=0 176 : DO IVAL=1,3 177 : VALP(IVAL)=VAL(IPERM(IVAL)) 178 : NUMP(IVAL)=NUM(IPERM(IVAL),IEL) 179 : LVAL(IVAL)=(LEGISO(VALP(IVAL))) 180 : IF (LVAL(IVAL)) THEN 181 : NPEQ1=NPEQ1+1 182 : ENDIF 183 : ENDDO 184 : C WRITE(IOIMP,*) 'IEL=',IEL 185 : C DO I=1,3 186 : C WRITE(IOIMP,*) I,' ',NUM(I,IEL),' ',VAL(I), 187 : C $ ' ',IPERM(I),' ',LVAL(I) 188 : C ENDDO 189 : IF (NPEQ1.EQ.3) THEN 190 : * Le triangle entier 191 : MLENTI=ELEM(3) 192 : LECT(**)=NUMP(1) 193 : LECT(**)=NUMP(2) 194 : LECT(**)=NUMP(3) 195 : LECT(**)=ICOLOR(IEL) 196 : ELSEIF (NPEQ1.EQ.2) THEN 197 : * Un des côtés 198 : MLENTI=ELEM(2) 199 : IF (LVAL(1)) THEN 200 : LECT(**)=NUMP(1) 201 : ELSE 202 : LECT(**)=NUMP(3) 203 : ENDIF 204 : LECT(**)=NUMP(2) 205 : LECT(**)=ICOLOR(IEL) 206 : ELSEIF (NPEQ1.EQ.1) THEN 207 : * Un des points ou un segment 208 : IF (LVAL(1)) THEN 209 : * Le point P1 210 : MLENTI=ELEM(1) 211 : LECT(**)=NUMP(1) 212 : LECT(**)=ICOLOR(IEL) 213 : ELSEIF (LVAL(2)) THEN 214 : * Un segement reliant P2 à l'intersection sur (P1, P3) 215 : MLENTI=ELEM(2) 216 : LECT(**)=NUMP(2) 216 : CALL ISOVA3(XISO,VALP(1),VALP(3),NUMP(1),NUMP(3) 217 : $ ,MLENTI) 218 : LECT(**)=ICOLOR(IEL) 219 : ELSE 220 : * Le point P3 221 : MLENTI=ELEM(1) 222 : LECT(**)=NUMP(3) 223 : LECT(**)=ICOLOR(IEL) 224 : ENDIF 225 : ELSE 226 : * Un segment reliant (P1,P2) à (P1,P3) 227 : * ou (P2,P3) à (P1,P3) 228 : MLENTI=ELEM(2) 228 : CALL ISOVA3(XISO,VALP(1),VALP(3),NUMP(1),NUMP(3) 229 : $ ,MLENTI) 230 : IF (VALP(2).GT.XISO) THEN 231 : * Un segment reliant (P1,P2) à (P1,P3) 231 : CALL ISOVA3(XISO,VALP(1),VALP(2),NUMP(1),NUMP(2) 232 : $ ,MLENTI) 233 : 234 : ELSE 235 : * Un segment reliant (P2,P3) à (P1,P3) 235 : CALL ISOVA3(XISO,VALP(2),VALP(3),NUMP(2),NUMP(3) 236 : $ ,MLENTI) 237 : ENDIF 238 : LECT(**)=ICOLOR(IEL) 239 : ENDIF 240 : ENDIF 241 : ENDDO 242 : goto 900 243 : 244 : ******** Cas des CUB8 *********************************************** 245 : 114 CONTINUE 246 : * On remplit les tableaux MATCAS et MATARE decrivant les cas et aretes 247 : c write(6,*) 'ISOVA2: cub8: appel a ISOVA4' 247 : CALL ISOVA4(ITYPEL,MATCAS,NMAXCA,NMAXPT,MATARE,2,NMAXAR) 248 : c write(6,*) 'ISOVA2: MATARE=',(MATARE(iou,1),iou=1,2) 249 : c write(6,*) ' ',(MATARE(iou,2),iou=1,2) 250 : c write(6,*) ' ...' 251 : c write(6,*) 'ISOVA2: MATCAS=',(MATCAS(iou,1),iou=1,16) 252 : c write(6,*) ' ',(MATCAS(iou,2),iou=1,16) 253 : c write(6,*) ' ...' 254 : c write(6,*) ' ',(MATCAS(iou,255),iou=1,16) 255 : c write(6,*) ' ',(MATCAS(iou,256),iou=1,16) 256 : 257 : * REM (BP) : 258 : * - on ne fait pas de test avec la tolerance (cf LEGISO autres itypel) et 259 : * on espere que VAL(IVAL).le.XISO au lieu de VAL(IVAL).lt.XISO suffira... 260 : * - on ne construit que des triangles TRI3 en se basant sur l algo : 261 : * "the marching cubes algorithm" (domaine public) 262 : DO 214 IEL=1,NEL 263 : ICAS = 1 264 : DO IVAL=1,8 265 : VAL(IVAL)=VELCHE(MIN(N1PTEL,IVAL),MIN(N1EL,IEL)) 266 : c if(VAL(IVAL).le.XISO) ICAS=ICAS+(2**(IVAL-1)) 267 : ENDDO 268 : c on teste les noeuds (attention: numerotation cub8 cast3m 269 : c differente de celle du marching cubes algorithm) 270 : if(VAL(1).le.XISO) ICAS=ICAS+1 271 : if(VAL(4).le.XISO) ICAS=ICAS+2 272 : if(VAL(3).le.XISO) ICAS=ICAS+4 273 : if(VAL(2).le.XISO) ICAS=ICAS+8 274 : if(VAL(5).le.XISO) ICAS=ICAS+16 275 : if(VAL(8).le.XISO) ICAS=ICAS+32 276 : if(VAL(7).le.XISO) ICAS=ICAS+64 277 : if(VAL(6).le.XISO) ICAS=ICAS+128 278 : c write(6,*) 'ICAS=',ICAS 279 : c write(6,*) 'MATCAS=',(MATCAS(iou,ICAS),iou=1,16) 280 : * Y a-t-il une isovaleur ? 281 : IF (ICAS.gt.1.or.ICAS.lt.256) THEN 282 : MLENTI=ELEM(3) 283 : do iar1=1,15 284 : NUMAR1=MATCAS(iar1,ICAS) 285 : c write(6,*) 'iar1,NUMAR1=',iar1,NUMAR1 286 : if(NUMAR1.eq.0) goto 214 287 : ipoin1 = MATARE(1,NUMAR1) 288 : ipoin2 = MATARE(2,NUMAR1) 289 : VAL1=VAL(ipoin1) 290 : VAL2=VAL(ipoin2) 291 : NUM1=NUM(ipoin1,IEL) 292 : NUM2=NUM(ipoin2,IEL) 293 : c write(6,*) ' ipoin1 et 2 =',ipoin1,ipoin2,NUM1,NUM2 293 : CALL ISOVA3(XISO,VAL1,VAL2,NUM1,NUM2,MLENTI) 294 : c write(6,*) ' VAL1 et 2=',VAL1,VAL2 295 : if(mod(iar1,3).eq.0) LECT(**)=ICOLOR(IEL) 296 : enddo 297 : ENDIF 298 : 214 CONTINUE 299 : c write(6,*) '=> LECT=',(LECT(iou),iou=1,LECT(/1)) 300 : goto 900 301 : 302 : ******** Cas des TET4 *********************************************** 303 : 123 CONTINUE 304 : DO IEL=1,NEL 305 : DO IVAL=1,4 306 : VAL(IVAL)=VELCHE(MIN(N1PTEL,IVAL),MIN(N1EL,IEL)) 307 : ENDDO 308 : VMIN=MIN(VAL(1),VAL(2),VAL(3),VAL(4)) 309 : VMAX=MAX(VAL(1),VAL(2),VAL(3),VAL(4)) 310 : * Y a-t-il une isovaleur ? 311 : IF (((VMIN-XISO).LE.XTOL).AND.((VMAX-XISO).GE.-XTOL)) 312 : $ THEN 313 : * Permutation permettant d'avoir VAL en ordre croissant 314 : * VAL(IPERM(1)).LE.VAL(IPERM(2)).LE.VAL(IPERM(3)).LE.VAL(IPERM(4)) 315 : DO IVAL=1,4 316 : IPERM(IVAL)=IVAL 317 : ENDDO 317 : CALL ISHELR(4,IPERM,4,VAL,IMPR,IRET) 318 : NPEQ1=0 319 : DO IVAL=1,4 320 : VALP(IVAL)=VAL(IPERM(IVAL)) 321 : NUMP(IVAL)=NUM(IPERM(IVAL),IEL) 322 : LVAL(IVAL)=(LEGISO(VALP(IVAL))) 323 : IF (LVAL(IVAL)) THEN 324 : NPEQ1=NPEQ1+1 325 : ENDIF 326 : ENDDO 327 : C WRITE(IOIMP,*) 'IEL=',IEL 328 : C DO I=1,3 329 : C WRITE(IOIMP,*) I,' ',NUM(I,IEL),' ',VAL(I), 330 : C $ ' ',IPERM(I),' ',LVAL(I) 331 : C ENDDO 332 : IF (NPEQ1.EQ.4) THEN 333 : * Isovaleur = tetrahedre 334 : MLENTI=ELEM(4) 335 : DO IVAL=1,4 336 : LECT(**)=NUMP(IVAL) 337 : ENDDO 338 : LECT(**)=ICOLOR(IEL) 339 : ELSEIF (NPEQ1.EQ.3) THEN 340 : * Isovaleur = 1 face du tetrahedre 341 : MLENTI=ELEM(3) 342 : IF (LVAL(1)) THEN 343 : LECT(**)=NUMP(1) 344 : ELSE 345 : LECT(**)=NUMP(4) 346 : ENDIF 347 : LECT(**)=NUMP(2) 348 : LECT(**)=NUMP(3) 349 : LECT(**)=ICOLOR(IEL) 350 : ELSEIF (NPEQ1.EQ.2) THEN 351 : * Isovaleur = 1 segment ou 1 triangle 352 : * appuye sur 2 points du tetrahedre 353 : IF (LVAL(1).OR.LVAL(4)) THEN 354 : * Isovaleur = 1 segment 355 : MLENTI=ELEM(2) 356 : DO IVAL=1,4 357 : IF (LVAL(IVAL)) THEN 358 : LECT(**)=NUMP(IVAL) 359 : ENDIF 360 : ENDDO 361 : LECT(**)=ICOLOR(IEL) 362 : ELSE 363 : * Isovaleur = 1 triangle appuye sur 2 points tetra 364 : MLENTI=ELEM(3) 365 : LECT(**)=NUMP(2) 366 : LECT(**)=NUMP(3) 366 : CALL ISOVA3(XISO,VALP(1),VALP(4),NUMP(1),NUMP(4) 367 : $ ,MLENTI) 368 : LECT(**)=ICOLOR(IEL) 369 : ENDIF 370 : ELSEIF (NPEQ1.EQ.1) THEN 371 : * Isovaleur = 1 point ou 1 triangle 372 : * appuye sur 1 point du tetrahedre 373 : IF (LVAL(1).OR.LVAL(4)) THEN 374 : * Isovaleur = 1 point 375 : MLENTI=ELEM(1) 376 : DO IVAL=1,4 377 : IF (LVAL(IVAL)) THEN 378 : LECT(**)=NUMP(IVAL) 379 : ENDIF 380 : ENDDO 381 : LECT(**)=ICOLOR(IEL) 382 : ELSE 383 : * Isosurface = 1 triangle appuye sur 1 point 384 : * du tetrahedre 385 : MLENTI=ELEM(3) 386 : IF (LVAL(2)) THEN 387 : * Le point appuye est P2 388 : * les autres points intersection (P1,P3) 389 : LECT(**)=NUMP(2) 389 : CALL ISOVA3(XISO,VALP(1),VALP(3),NUMP(1) 390 : $ ,NUMP(3),MLENTI) 391 : ELSE 392 : * Le point appuye est P3 393 : * les autres points intersection (P2,P4) 394 : LECT(**)=NUMP(3) 394 : CALL ISOVA3(XISO,VALP(2),VALP(4),NUMP(2) 395 : $ ,NUMP(4),MLENTI) 396 : ENDIF 397 : * Il y a toujours un point intersection (P1,P4) 397 : CALL ISOVA3(XISO,VALP(1),VALP(4),NUMP(1) 398 : $ ,NUMP(4),MLENTI) 399 : LECT(**)=ICOLOR(IEL) 400 : ENDIF 401 : ELSE 402 : * Isosurface = 1 triangle quelconque section 403 : * du tetrahedre ou un quadrangle (= 2 triangles) 404 : IF (VALP(2).GT.XISO) THEN 405 : * Isosurface entre V1 et V2 406 : * Points intersection (P1,P2) (P1,P3) 407 : * Un seul triangle 408 : MLENTI=ELEM(3) 408 : CALL ISOVA3(XISO,VALP(1),VALP(2),NUMP(1) 409 : $ ,NUMP(2),MLENTI) 409 : CALL ISOVA3(XISO,VALP(1),VALP(3),NUMP(1) 410 : $ ,NUMP(3),MLENTI) 411 : ELSEIF (VALP(3).LT.XISO) THEN 412 : * Isosurface entre V3 et V4 413 : * Points intersection (P3,P4) (P2,P4) 414 : * Un seul triangle 415 : MLENTI=ELEM(3) 415 : CALL ISOVA3(XISO,VALP(3),VALP(4),NUMP(3) 416 : $ ,NUMP(4),MLENTI) 416 : CALL ISOVA3(XISO,VALP(2),VALP(4),NUMP(2) 417 : $ ,NUMP(4),MLENTI) 418 : ELSE 419 : * Isosurface entre V2 et V3 ; 420 : * Points intersection (P1,P3) (P2,P3) (P2,P4) 421 : * Un quadrangle = 2 triangles 422 : MLENTI=ELEM(5) 422 : CALL ISOVA3(XISO,VALP(1),VALP(3),NUMP(1) 423 : $ ,NUMP(3),MLENTI) 423 : CALL ISOVA3(XISO,VALP(2),VALP(3),NUMP(2) 424 : $ ,NUMP(3),MLENTI) 424 : CALL ISOVA3(XISO,VALP(2),VALP(4),NUMP(2) 425 : $ ,NUMP(4),MLENTI) 426 : ENDIF 427 : * Il y toujours un point intersection 428 : * entre P1 et P4 428 : CALL ISOVA3(XISO,VALP(1),VALP(4),NUMP(1) 429 : $ ,NUMP(4),MLENTI) 430 : LECT(**)=ICOLOR(IEL) 431 : ENDIF 432 : ENDIF 433 : ENDDO 434 : goto 900 435 : 436 : 437 : ******** on termine proprement ************************************** 438 : 900 CONTINUE 439 : SEGDES MELVAL 440 : 441 : ENDDO 442 : c fin de Boucle sur les zones elementaires 443 : C*********************************************************************** 444 : 445 : RETURN 446 : * 447 : * End of subroutine ISOVA2 448 : * 449 : END 450 : 451 : 452 : 453 :
© Cast3M 2003 - All rights reserved.
Disclaimer