1 : C VORO1 SOURCE SH236661 13/11/25 21:15:12 7869 2 : C 3 : SUBROUTINE VORO1(IPT1,IPT2,LNBOIT,ILEA1,ITRC1,DDMAX, 4 : & MTAB1) 5 : C----------------------------------------------------------------------C 6 : C Partition de Voronoi (non coupee) d'un nuage de points (IPT1) a C 7 : C partir de la triangulation de Delaunay (IPT2) C 8 : C C 9 : C IPT1 = pointeur sur le maillage de points, centres des C 10 : C polyedres (elements POI1) C 11 : C IPT2 = pointeur sur le maillage de la triangulation de C 12 : C Delaunay de IPT1 (avec les coins de la boite de C 13 : C triangulation) (elements TRI3 ou TET4) C 14 : C LNBOIT = tableau contenant les numeros des noeuds formant les C 15 : C coins de la boite de triangulation C 16 : C ILEA1 = pointeur sur le segment MADJACEN contenant la table des C 17 : C connectivite des triangles de IPT1 C 18 : C ITRC1 = pointeur sur le segment MCIRCONS contenant la table des C 19 : C coordonnees des cercles circonscrits des triangles de C 20 : C IPT1 C 21 : C XZERO = critere de distance pour considerer deux centres de C 22 : C cercles circonscrits confondus C 23 : C C 24 : C MTAB1 = table de maillage des polyedres non coupes C 25 : C----------------------------------------------------------------------C 26 : IMPLICIT INTEGER(I-N) 27 : IMPLICIT REAL*8 (A-H,O-Z) 28 : c 29 : -INC CCOPTIO 30 : -INC SMELEME 31 : -INC SMCOORD 32 : -INC CCGEOME 33 : -INC SMLENTI 34 : -INC SMTABLE 35 : DIMENSION LNBOIT(8) 36 : DIMENSION XPI1(3) 37 : LOGICAL BOOL0,BOOL1,BOOL2,BOOL3 38 : C 39 : POINTEUR LCCC1.MLENTI,LAR1.MLENTI,LAR2.MLENTI 40 : C Tables, maillages et listes pour la table de sortie 41 : POINTEUR IPTG.MELEME,MTABC.MTABLE,MTABI1.MTABLE 42 : POINTEUR IPTI1.MELEME,LFI1.MLENTI,IPTVI1.MELEME,MTABF.MTABLE 43 : POINTEUR IPTF1.MELEME,LA1.MLENTI,MTABA.MTABLE,IPTA1.MELEME 44 : C 45 : SEGMENT,MADJACEN 46 : INTEGER LEAC(NBL1,IDIM+1,2) 47 : ENDSEGMENT 48 : POINTEUR ILEA1.MADJACEN 49 : C 50 : SEGMENT,MCIRCONS 51 : REAL*8 TRC(NBE1,4) 52 : ENDSEGMENT 53 : POINTEUR ITRC1.MCIRCONS 54 : C 55 : SEGMENT,MVORO 56 : INTEGER LACT(N1,N2) 57 : ENDSEGMENT 58 : POINTEUR ITL1.MVORO,ITCEL1.MVORO,ITCEL2.MVORO,ITARE1.MVORO 59 : POINTEUR ITFAC1.MVORO,ITFAC2.MVORO 60 : 61 : C 62 : XZERO = DDMAX*1.D-10 63 : C 64 : C 1) Pre-traitement : creation des points representant les centres 65 : C cercles circonscrits de chaque triangle (le meme numero est 66 : C attribue aux centres de cercles circonscrits superposes) 67 : C 68 : C Ajustement initial de MCOORD au nombre de cercles circonscrits 69 : NBTRI=ITRC1.TRC(/1) 70 : NBPTS0=XCOOR(/1)/(IDIM+1) 71 : NBPTS=NBPTS0+NBTRI 72 : SEGADJ,MCOORD 73 : C LCCC1.LECT(i)=j le triangle i a pour Centre de Cerclce Circonscrit 74 : C le point numero j 75 : JG=NBTRI 76 : SEGINI,LCCC1 77 : I1=0 78 : C Boucle sur les centres des cercles circonscrits des triangles 79 : DO I=1,NBTRI 80 : DO K=1,IDIM 81 : XPI1(K)=ITRC1.TRC(I,K+1) 82 : ENDDO 83 : DO J=1,NBTRI 84 : NCJ=0 85 : IF (J.GE.I) GOTO 101 86 : C Test si les coordonnees de I et J sont egales 87 : DIJ1=0. 88 : DO K=1,IDIM 89 : XPJ=ITRC1.TRC(J,K+1) 90 : DIJ1=DIJ1+(XPI1(K)-XPJ)*(XPI1(K)-XPJ) 91 : c XDEN=XPI1(K) 92 : c IF (XDEN.EQ.0.) XDEN=1. 93 : c IF (ABS((XPI1(K)-XPJ)/XDEN).LT.XZERO) NCJ=NCJ+1 94 : ENDDO 95 : IF(DIJ1.LT.XZERO) NCJ=IDIM 96 : C Si les centres I et J sont superposes 97 : IF (NCJ.EQ.IDIM) THEN 98 : LCCC1.LECT(I)=LCCC1.LECT(J) 99 : GOTO 102 100 : ENDIF 101 : 101 CONTINUE 102 : ENDDO 103 : C Creation d'un nouveau point dans MCOORD 104 : I1=I1+1 105 : NREFI=NBPTS0+I1 106 : LCCC1.LECT(I)=NREFI 107 : C Ecriture des ses coordonnees 108 : DO J=1,IDIM 109 : XCOOR(((NREFI-1)*(IDIM+1))+J)=XPI1(J) 110 : ENDDO 111 : C et de la densite associee, choisie nulle 112 : XCOOR(((NREFI-1)*(IDIM+1))+(IDIM+1))=0. 113 : 102 CONTINUE 114 : ENDDO 115 : C Ajustement final de MCOORD (on a ajoute I1 points) 116 : IF (I1.NE.NBTRI) THEN 117 : NBPTS=NBPTS0+I1 118 : SEGADJ,MCOORD 119 : ENDIF 120 : C 121 : C 122 : C 2.0) Initialisations avant creation du maillage des polyedres 123 : C IPTG : Maillage des aretes des polyedres (elements SEG2) 124 : NBSOUS=0 125 : NBREF=0 126 : NBNN=2 127 : NBELEM=I1*4 128 : SEGINI,IPTG 129 : IPTG.ITYPEL=2 130 : NBARETES=0 131 : C Table MTAB1 de sortie contenant les maillages des polyedres 131 : CALL CRTABL(MTAB1) 131 : CALL ECCTAB(MTAB1,'MOT ',0,0,'VISU',0,0, 132 : & 'MAILLAGE',0,0,0,0,IPTG) 132 : CALL CRTABL(MTABC) 132 : CALL ECCTAB(MTAB1,'MOT ',0,0,'CELL',0,0, 133 : & 'TABLE ',0,0,0,0,MTABC) 134 : IF (IDIM.EQ.3) THEN 134 : CALL CRTABL(MTABF) 134 : CALL ECCTAB(MTAB1,'MOT ',0,0,'FACS',0,0, 135 : & 'TABLE ',0,0,0,0,MTABF) 136 : ENDIF 136 : CALL CRTABL(MTABA) 136 : CALL ECCTAB(MTAB1,'MOT ',0,0,'ARTS',0,0, 137 : & 'TABLE ',0,0,0,0,MTABA) 138 : C ITL1 : Table indiquant si la ligne entre deux centres de cercle 139 : C circonscrits existe deja ou non 140 : N1=I1 141 : N2=I1 142 : SEGINI,ITL1 143 : C ITCEL1 : Table des polyedres voisins d'un polyedre 144 : C ITCEL2 : Table des aretes d'un polyedre 145 : N1=10 146 : N2=IPT1.NUM(/2) 147 : SEGINI,ITCEL1,ITCEL2 148 : C ITARE1 : Table des cellules appuyees sur une arete, on la stocke 149 : C dans la table de sortie pour pouvoir la retrouver lors du 150 : C decoupage dans VORO2 151 : C ITARE1.LACT(ICE,NLA)=I 152 : C NLA : numero (dans IPTG) de l'arete 153 : C ICE : index des cellules appuyees sur NLA (de 1 a n, puis 0) 154 : C I : numero (dans IPT1) de la ICE-ieme cellule appuyee sur NLA 155 : N1=10 156 : N2=IPTG.NUM(/2) 157 : SEGINI,ITARE1 157 : CALL ECCTAB(MTAB1,'MOT ',0,0,'TAR',0,0, 158 : & 'MVORO ',0,0,0,0,ITARE1) 159 : C 160 : C 161 : C 162 : C 2.1) Maillage des cellules de Voronoi en dimension 2 163 : C Boucle sur les centres des polyedres de Voronoi 164 : IF (IDIM.EQ.2) THEN 165 : DO I=1,IPT1.NUM(/2) 166 : NREF_I=IPT1.NUM(1,I) 167 : C Initialisation des indices pour la cellule NREF_I 167 : CALL CRTABL(MTABI1) 167 : CALL ECCTAB(MTABC,'POINT ',0,0,' ',0,NREF_I, 168 : & 'TABLE ',0,0,0,0,MTABI1) 169 : C --> maillage de la cellule 170 : NBNN=2 171 : NBELEM=I1*4 172 : SEGINI,IPTI1 173 : IPTI1.ITYPEL=2 174 : NBELEI1=0 174 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 175 : & 'MAILLAGE',0,0,0,0,IPTI1) 176 : C --> liste des aretes 177 : JG=I1*4 178 : SEGINI,LA1 178 : CALL ECCTAB(MTABI1,'MOT ',0,0,'ARTS',0,0, 179 : & 'LISTENTI',0,0,0,0,LA1) 180 : C --> maillage des voisns 181 : NBNN=1 182 : NBELEM=I1*4 183 : SEGINI,IPTVI1 184 : IPTVI1.ITYPEL=1 185 : NVI1=0 185 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, 186 : & 'MAILLAGE',0,0,0,0,IPTVI1) 187 : C Recherche du premier triangle appuye sur NREFA 188 : NREFA=IPT1.NUM(1,I) 189 : DO J=1,IPT2.NUM(/2) 190 : DO K=1,IPT2.NUM(/1) 191 : IF (IPT2.NUM(K,J).EQ.NREFA) GOTO 211 192 : ENDDO 193 : ENDDO 194 : 211 NT0=J 195 : NA=K 196 : NB=NA+1 197 : IF (NB.EQ.4) NB=1 198 : NC=6-NA-NB 199 : NREFB=IPT2.NUM(NB,NT0) 200 : NREFC=IPT2.NUM(NC,NT0) 201 : C Balayage des triangles appuyes sur NREFA en tournant autour 202 : C de ce point, a patrir du triangle NT0 203 : C La boucle s'arrete quand on est revenu au triangle NT0 204 : NT1=NT0 205 : NRCC1=LCCC1.LECT(NT1) 206 : NCVI=0 207 : NAI=0 208 : DO J=1,I1*4 209 : IJOK=0 210 : C Determination du point NRCC2 a relier a NRCC1 211 : NT2=ILEA1.LEAC(NT1,NB,1) 212 : NRCC2=LCCC1.LECT(NT2) 213 : C Si la liaison entre CC1 et CC2 existe, on la recupere et 214 : C on passe au triangle suivant 215 : IF (ITL1.LACT(NRCC1-NBPTS0,NRCC2-NBPTS0).EQ.1) THEN 216 : C On determine le numero de l'arete en question 217 : DO K=1,ITCEL2.LACT(/1) 218 : NVI=0 219 : DO L=1,IPT1.NUM(/2) 220 : IF (IPT1.NUM(1,L).EQ.NREFC) NVI=L 221 : ENDDO 222 : IF (NVI.EQ.0) THEN 223 : DO L=1,4 224 : IF (LNBOIT(L).EQ.NREFC) NVI=-1*L 225 : ENDDO 226 : ENDIF 227 : NAR1=ITCEL2.LACT(K,NVI) 228 : NPA=IPTG.NUM(1,NAR1) 229 : NPB=IPTG.NUM(2,NAR1) 230 : IF (((NPA.EQ.NRCC1).AND.(NPB.EQ.NRCC2)).OR. 231 : & ((NPA.EQ.NRCC2).AND.(NPB.EQ.NRCC1))) THEN 232 : GOTO 212 233 : ENDIF 234 : ENDDO 235 : ENDIF 236 : C Si les centres CC1 et CC2 sont confondus, il n'y a pas de 237 : C liaison a creer, on passe au triangle suivant 238 : IF (NRCC1.EQ.NRCC2) THEN 239 : GOTO 214 240 : ENDIF 241 : C Dans les autres cas, il y a bien creation d'une nouvelle 242 : C arete dans le maillage IPTG de la cellule I 243 : IJOK=1 244 : NBARETES=NBARETES+1 245 : IPTG.NUM(1,NBARETES)=NRCC1 246 : IPTG.NUM(2,NBARETES)=NRCC2 247 : NAR1=NBARETES 248 : C Remplissage de la table ILT1 des liaisons entre centres 249 : C de cercles circonscrits 250 : ITL1.LACT(NRCC1-NBPTS0,NRCC2-NBPTS0)=1 251 : ITL1.LACT(NRCC2-NBPTS0,NRCC1-NBPTS0)=1 252 : C Remplissage de la table ITCEL1 des numeros des cellules 253 : C voisines de la cellule I 254 : NVI=0 255 : DO K=1,IPT1.NUM(/2) 256 : IF (IPT1.NUM(1,K).EQ.NREFC) NVI=K 257 : ENDDO 258 : IF (NVI.EQ.0) THEN 259 : DO K=1,4 260 : IF (LNBOIT(K).EQ.NREFC) NVI=-1*K 261 : ENDDO 262 : ENDIF 263 : 212 IF (NVI.GT.0) THEN 264 : NCVI=NCVI+1 265 : IF (NCVI.GT.ITCEL1.LACT(/1)) THEN 266 : N1=ITCEL1.LACT(/1)+10 267 : N2=ITCEL1.LACT(/2) 268 : SEGADJ,ITCEL1 269 : ENDIF 270 : ITCEL1.LACT(NCVI,I)=NVI 271 : ENDIF 272 : C Remplissage de la table ITCEL2 des numeros d'aretes 273 : C appartenant a la cellule I 274 : NAI=NAI+1 275 : IF (NAI.GT.ITCEL2.LACT(/1)) THEN 276 : N1=ITCEL2.LACT(/1)+10 277 : N2=ITCEL2.LACT(/2) 278 : SEGADJ,ITCEL2 279 : ENDIF 280 : ITCEL2.LACT(NAI,I)=NAR1 281 : C Remplissage de la table ITARE1 des numeros de cellules 282 : C contenant l'arrete creee/recuperee 283 : ICE=0 284 : DO K=1,ITARE1.LACT(/1) 285 : IF (ITARE1.LACT(K,NAR1).EQ.0) THEN 286 : ICE=K 287 : GOTO 213 288 : ENDIF 289 : ENDDO 290 : IF (ICE.EQ.0) THEN 291 : ICE=ITARE1.LACT(/1)+1 292 : N1=ITARE1.LACT(/1)+10 293 : N2=ITARE1.LACT(/2) 294 : SEGADJ,ITARE1 295 : ENDIF 296 : 213 ITARE1.LACT(ICE,NAR1)=I 297 : C Ajout d'un element au maillage IPTI1 de la cellule I 298 : NBELEI1=NBELEI1+1 299 : IPTI1.NUM(1,NBELEI1)=NRCC1 300 : IPTI1.NUM(2,NBELEI1)=NRCC2 301 : C Ajout du numero d'arete dans la liste de la cellule I 302 : LA1.LECT(NBELEI1)=NAR1 303 : C Ajout d'un voisin a la cellule I 304 : IF (NVI.GT.0) THEN 305 : NVI1=NVI1+1 306 : NREFVI=IPT1.NUM(1,NVI) 307 : IPTVI1.NUM(1,NVI1)=NREFVI 308 : ENDIF 309 : C Ajout de l'arete dans MTABA 310 : IF (IJOK.EQ.1) THEN 311 : NBNN=2 312 : NBELEM=1 313 : SEGINI,IPTA1 314 : IPTA1.ITYPEL=2 315 : IPTA1.NUM(1,1)=NRCC1 316 : IPTA1.NUM(2,1)=NRCC2 316 : CALL ECCTAB(MTABA,'ENTIER ',NBARETES,0,' ',0,0, 317 : & 'MAILLAGE',0,0,0,0,IPTA1) 318 : ENDIF 319 : C Test d'arret de la boucle (est on revenu sur le triangle 320 : C de depart ?) 321 : 214 IF (NT2.EQ.NT0) THEN 322 : GOTO 215 323 : ENDIF 324 : C Nouveau points A, B et C pour le prochain triangle 325 : DO K=1,IPT2.NUM(/1) 326 : IF (IPT2.NUM(K,NT2).EQ.NREFA) NA=K 327 : IF (IPT2.NUM(K,NT2).EQ.NREFC) NB=K 328 : ENDDO 329 : NC=6-NA-NB 330 : NREFB=IPT2.NUM(NB,NT2) 331 : NREFC=IPT2.NUM(NC,NT2) 332 : NT1=NT2 333 : NRCC1=LCCC1.LECT(NT1) 334 : ENDDO 335 : 215 CONTINUE 336 : C Ajustement des maillages de la cellule I, de ses voisins et 337 : C de la liste des aretes 338 : NBNN=IPTI1.NUM(/1) 339 : NBELEM=NBELEI1 340 : SEGADJ,IPTI1 341 : JG=NBELEI1 342 : SEGADJ,LA1 343 : NBNN=IPTVI1.NUM(/1) 344 : NBELEM=NVI1 345 : SEGADJ,IPTVI1 346 : ENDDO 347 : C 348 : C 349 : C 350 : C 2.2) Maillage des cellules de Voronoi en dimension 3 351 : ELSEIF (IDIM.EQ.3) THEN 352 : C ITFAC1 : tableau indiquant si la face (I,J) existe deja 353 : C ou I et J sont des numeros de cellule 354 : N1=IPT1.NUM(/2) 355 : N2=N1 356 : SEGINI,ITFAC1 357 : C ITFAC2 : tableau indiquant si la face (I,J) existe deja 358 : C ou I est un numero de cellule exterieures infinie et J un 359 : C numero des coins de la boite de triangulation 360 : N1=IPT1.NUM(/2) 361 : N2=8 362 : SEGINI,ITFAC2 363 : C Boucle sur les centres des polyedres de Voronoi 364 : NBFACES=0 365 : DO I=1,IPT1.NUM(/2) 366 : C On va chercher la table de la cellule I 367 : NREF_I=IPT1.NUM(1,I) 368 : C Initialisation des indices pour la cellule NREF_I 368 : CALL CRTABL(MTABI1) 368 : CALL ECCTAB(MTABC,'POINT ',0,0,' ',0,NREF_I, 369 : & 'TABLE ',0,0,0,0,MTABI1) 370 : C --> maillage de la cellule 371 : NBNN=2 372 : NBELEM=I1*4 373 : SEGINI,IPTI1 374 : IPTI1.ITYPEL=2 375 : NBELEI1=0 375 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 376 : & 'MAILLAGE',0,0,0,0,IPTI1) 377 : C --> liste des faces 378 : JG=I1*4 379 : SEGINI,LFI1 380 : NFI1=0 380 : CALL ECCTAB(MTABI1,'MOT ',0,0,'FACS',0,0, 381 : & 'LISTENTI',0,0,0,0,LFI1) 382 : C --> maillage des voisns 383 : NBNN=1 384 : NBELEM=I1*4 385 : SEGINI,IPTVI1 386 : IPTVI1.ITYPEL=1 387 : IV1=0 387 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, 388 : & 'MAILLAGE',0,0,0,0,IPTVI1) 389 : C Recherche du premier tetraedre appuye sur NREFA 390 : NCVI=0 391 : NREFA=IPT1.NUM(1,I) 392 : DO J=1,IPT2.NUM(/2) 393 : DO K=1,IPT2.NUM(/1) 394 : IF (IPT2.NUM(K,J).EQ.NREFA) GOTO 221 395 : ENDDO 396 : ENDDO 397 : 221 NT0=J 398 : NA=K 399 : NB=NA+1 400 : NC=NA-1 401 : IF (NB.EQ.5) NB=1 402 : IF (NC.EQ.0) NC=3 403 : ND=10-NA-NB-NC 404 : C Liste des aretes appuyees sur NREFA autour desquelles on 405 : C va tourner. On incremente ces liste au fur et a mesure que 406 : C l'on rencontre de nouveau tetraedres appuyes sur NT0 407 : JG=IPTG.NUM(/2) 408 : SEGINI,LAR1,LAR2 409 : C On connait les trois premieres aretes appuyees sur NREFA 410 : LAR1.LECT(1)=IPT2.NUM(NB,NT0) 411 : LAR1.LECT(2)=IPT2.NUM(NC,NT0) 412 : LAR1.LECT(3)=IPT2.NUM(ND,NT0) 413 : C et les numeros des premiers tetraedres contenant ces aretes 414 : LAR2.LECT(1)=NT0 415 : LAR2.LECT(2)=NT0 416 : LAR2.LECT(3)=NT0 417 : C Balayage des aretes des tetraedres appuyes sur NREFA, 418 : C a patrir du tetraedre NT0 (on parcoure LAR1) 419 : C Chaque nouvelle arete definit potentielement une nouvelle 420 : C face de la partition de Voronoi 421 : C La liste LAR1 s'incremente, la boucle s'arrete lorsque l'on 422 : C a parcouru toutes les aretes dans LAR1 423 : DO J=1,LAR1.LECT(/1) 424 : IJOK=0 425 : C On tourne autour de l'arete J de la liste LAR1 426 : NREFB=LAR1.LECT(J) 427 : C Test d'arret de la boucle, si on a fini de parcourir la 428 : C liste LAR1, on a fini ! 429 : IF (NREFB.EQ.0) THEN 430 : GOTO 229 431 : ENDIF 432 : C Calcul de NVI, numero de cette cellule voisine 433 : NVI=0 434 : DO K=1,IPT1.NUM(/2) 435 : IF (IPT1.NUM(1,K).EQ.NREFB) NVI=K 436 : ENDDO 437 : C ou bien s'il s'agit d'un sommet de la boite de 438 : C triangulation 439 : IF (NVI.EQ.0) THEN 440 : DO K=1,8 441 : IF (LNBOIT(K).EQ.NREFB) NVI=-1*K 442 : ENDDO 443 : ENDIF 444 : C NFIJ est le numero de la face entre I et J 445 : C (si elle existe deja) 446 : IF (NVI.GT.0) THEN 447 : NFIJ=ITFAC1.LACT(NVI,I) 448 : ELSE 449 : NVII=-1*NVI 450 : NFIJ=ITFAC2.LACT(I,NVII) 451 : ENDIF 452 : C IJOK est mis a 1 quand on identifie une face non deja 453 : C existante (mais elle peut etre vide !) 454 : IF (NFIJ.EQ.0) IJOK=1 455 : C Dans le cas d'une potentielle nouvelle face, on prepare 456 : C son maillage et sa liste des aretes 457 : IF (IJOK.EQ.1) THEN 458 : NBNN=2 459 : NBELEM=LAR1.LECT(/1) 460 : SEGINI,IPTF1 461 : IPTF1.ITYPEL=2 462 : JG=LAR1.LECT(/1) 463 : SEGINI,LA1 464 : ENDIF 465 : NBELEF1=0 466 : C Le premier tetraedre balaye est NT0 467 : NT0=LAR2.LECT(J) 468 : C Calcul de NA,NB,NC et ND dans NT0 469 : DO K=1,IPT2.NUM(/1) 470 : IF (IPT2.NUM(K,NT0).EQ.NREFA) NA=K 471 : IF (IPT2.NUM(K,NT0).EQ.NREFB) NB=K 472 : ENDDO 473 : IF ((NA.EQ.1).AND.(NB.EQ.2)) NC=3 474 : IF ((NA.EQ.1).AND.(NB.EQ.3)) NC=2 475 : IF ((NA.EQ.1).AND.(NB.EQ.4)) NC=2 476 : IF ((NA.EQ.2).AND.(NB.EQ.1)) NC=3 477 : IF ((NA.EQ.2).AND.(NB.EQ.3)) NC=1 478 : IF ((NA.EQ.2).AND.(NB.EQ.4)) NC=1 479 : IF ((NA.EQ.3).AND.(NB.EQ.1)) NC=2 480 : IF ((NA.EQ.3).AND.(NB.EQ.2)) NC=1 481 : IF ((NA.EQ.3).AND.(NB.EQ.4)) NC=1 482 : IF ((NA.EQ.4).AND.(NB.EQ.1)) NC=2 483 : IF ((NA.EQ.4).AND.(NB.EQ.2)) NC=1 484 : IF ((NA.EQ.4).AND.(NB.EQ.3)) NC=1 485 : ND=10-NA-NB-NC 486 : NREFC=IPT2.NUM(NC,NT0) 487 : NREFD=IPT2.NUM(ND,NT0) 488 : C Boucle autour de l'arete NREFA-NREFB, en partant du 489 : C tetraedre NT0, on construit ainsi la face (I,J) en 490 : C tournant autour de ses aretes 491 : C La boucle s'arrete quand on est revenu au tetraedre NT0 492 : NT1=NT0 493 : NRCC1=LCCC1.LECT(NT1) 494 : DO K=1,I1*4 495 : NT2=ILEA1.LEAC(NT1,NC,1) 496 : C On ajoute NREFC et NREFD a la liste des aretes a 497 : C balayer plus tard (si non deja presents) 498 : DO L=1,LAR1.LECT(/1) 499 : IF (LAR1.LECT(L).EQ.NREFC) GOTO 222 500 : IF (LAR1.LECT(L).EQ.0) THEN 501 : LAR1.LECT(L)=NREFC 502 : LAR2.LECT(L)=NT2 503 : GOTO 222 504 : ENDIF 505 : ENDDO 506 : 222 CONTINUE 507 : DO L=1,LAR1.LECT(/1) 508 : IF (LAR1.LECT(L).EQ.NREFD) GOTO 223 509 : IF (LAR1.LECT(L).EQ.0) THEN 510 : LAR1.LECT(L)=NREFD 511 : LAR2.LECT(L)=NT2 512 : GOTO 223 513 : ENDIF 514 : ENDDO 515 : 223 CONTINUE 516 : NRCC2=LCCC1.LECT(NT2) 517 : C 518 : C Plusieurs cas sont possibles pour l'arete reliant 519 : C CC1 et CC2 : 520 : C 521 : C CAS 1 : la liaison entre CC1 et CC2 existe 522 : C => on la recupere et on passe au tetraedre suivant 523 : IF (ITL1.LACT(NRCC1-NBPTS0,NRCC2-NBPTS0).EQ.1) THEN 524 : C On cherche le numero de l'arete existante reliant 525 : C deja CC1 et CC2 526 : DO L=1,IPTG.NUM(/2) 527 : NPA=IPTG.NUM(1,L) 528 : NPB=IPTG.NUM(2,L) 529 : IF (((NPA.EQ.NRCC1).AND.(NPB.EQ.NRCC2)).OR. 530 : & ((NPA.EQ.NRCC2).AND.(NPB.EQ.NRCC1))) THEN 531 : NLA=L 532 : GOTO 224 533 : ENDIF 534 : ENDDO 535 : ENDIF 536 : C CAS 2 : les centres CC1 et CC2 sont confondus 537 : C => il n'y a pas d'arete a faire ! on passe au 538 : C tetraedre suivant 539 : IF (NRCC1.EQ.NRCC2) THEN 540 : GOTO 227 541 : ENDIF 542 : C CAS 3 : creation d'une nouvelle arete dans IPTG 543 : NBARETES=NBARETES+1 544 : IPTG.NUM(1,NBARETES)=NRCC1 545 : IPTG.NUM(2,NBARETES)=NRCC2 546 : NLA=NBARETES 547 : C et ajout de cette arete dans MTABA 548 : NBNN=2 549 : NBELEM=1 550 : SEGINI,IPTA1 551 : IPTA1.ITYPEL=2 552 : IPTA1.NUM(1,1)=NRCC1 553 : IPTA1.NUM(2,1)=NRCC2 553 : CALL ECCTAB(MTABA,'ENTIER ',NLA,0,' ',0,0, 554 : & 'MAILLAGE',0,0,0,0,IPTA1) 555 : C Remplissage de la table ILT1 des liaisons en centres 556 : C de cercles circonscrits 557 : ITL1.LACT(NRCC1-NBPTS0,NRCC2-NBPTS0)=1 558 : ITL1.LACT(NRCC2-NBPTS0,NRCC1-NBPTS0)=1 559 : C 560 : 224 CONTINUE 561 : C On incremente le nombre d'arete dans cette face 562 : NBELEF1=NBELEF1+1 563 : C Dans le cas d'une potentielle nouvelle face, on 564 : C incremente le maillage de la face et la liste de ses 565 : C aretes (si l'arete en question n'est pas presente) 566 : IF (IJOK.EQ.1) THEN 567 : DO L=1,LA1.LECT(/1) 568 : IF (LA1.LECT(L).EQ.NLA) GOTO 230 569 : IF (LA1.LECT(L).EQ.0) GOTO 231 570 : ENDDO 571 : 231 IF (NBELEF1.GT.LA1.LECT(/1)) THEN 572 : JG=LA1.LECT(/1)+10 573 : SEGADJ,LA1 574 : NBNN=IPTF1.NUM(/1) 575 : NBELEM=IPTF1.NUM(/2)+10 576 : SEGADJ,IPTF1 577 : ENDIF 578 : C ajout dans la liste des aretes pour cette face 579 : LA1.LECT(NBELEF1)=NLA 580 : C et au maillage de la face 581 : IPTF1.NUM(1,NBELEF1)=NRCC1 582 : IPTF1.NUM(2,NBELEF1)=NRCC2 583 : 230 CONTINUE 584 : ENDIF 585 : C Dans le cas d'une potentielle nouvelle face, si le 586 : C nombre d'aretes est egal a 3, il s'agit bel et bien 587 : C d'une face 588 : IF ((NBELEF1.EQ.3).AND.(IJOK.EQ.1)) THEN 589 : C on l'ajoute a MTABF 590 : NBFACES=NBFACES+1 591 : NFIJ=NBFACES 591 : CALL CRTABL(MTABF1) 591 : CALL ECCTAB(MTABF,'ENTIER ',NFIJ,0,' ',0,0, 592 : & 'TABLE ',0,0,0,0,MTABF1) 592 : CALL ECCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 593 : & 'MAILLAGE',0,0,0,0,IPTF1) 593 : CALL ECCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, 594 : & 'LISTENTI',0,0,0,0,LA1) 595 : ENDIF 596 : C Si le nombre d'arete est egal a 3, le point J est bel 597 : C et bien un voisin de I 598 : IF (NBELEF1.EQ.3) THEN 599 : C S'il s'agit d'une face exterieure (J est sommet de 600 : C la boite de triangulation) on l'indique juste dans 601 : C ITFAC2 602 : IF (NVI.LT.0) THEN 603 : NVII=-1*NVI 604 : ITFAC2.LACT(I,NVII)=NFIJ 605 : C Sinon, on ajoute un voisin et on marque la face 606 : C dans ITFAC1 607 : ELSE 608 : NREFVI=IPT1.NUM(1,NVI) 609 : IV1=IV1+1 610 : IF (IV1.GT.IPTVI1.NUM(/2)) THEN 611 : NBNN=IPTVI1.NUM(/1) 612 : NBELEM=IPTVI1.NUM(/2)+10 613 : SEGADJ,IPTVI1 614 : ENDIF 615 : IPTVI1.NUM(1,IV1)=NREFVI 616 : ITFAC1.LACT(I,NVI)=NFIJ 617 : ITFAC1.LACT(NVI,I)=NFIJ 618 : ENDIF 619 : C Ajout du numero global de la face dans la liste des 620 : C faces de la cellule NREF_I 621 : NFI1=NFI1+1 622 : IF (NFI1.GT.LFI1.LECT(/1)) THEN 623 : JG=LFI1.LECT(/1)+10 624 : SEGADJ,LFI1 625 : ENDIF 626 : LFI1.LECT(NFI1)=NFIJ 627 : ENDIF 628 : C Remplissage de la table ITCEL2 des numeros d'aretes 629 : C appartenant a la cellule I 630 : DO L=1,ITCEL2.LACT(/1) 631 : IF (ITCEL2.LACT(L,I).EQ.NLA) THEN 632 : GOTO 227 633 : ENDIF 634 : IF (ITCEL2.LACT(L,I).EQ.0) GOTO 225 635 : ENDDO 636 : 225 NBELEI1=NBELEI1+1 637 : IF (NBELEI1.GT.ITCEL2.LACT(/1)) THEN 638 : N1=ITCEL2.LACT(/1)+10 639 : N2=ITCEL2.LACT(/2) 640 : SEGADJ,ITCEL2 641 : ENDIF 642 : ITCEL2.LACT(NBELEI1,I)=NLA 643 : C Ajout de l'arete au maillage de la cellule NREF_I 644 : IPTI1.NUM(1,NBELEI1)=NRCC1 645 : IPTI1.NUM(2,NBELEI1)=NRCC2 646 : C Remplissage de la table ITARE1 des numeros de cellules 647 : C contenant l'arrete creee/recuperee 648 : ICE=0 649 : DO L=1,ITARE1.LACT(/1) 650 : IF (ITARE1.LACT(L,NLA).EQ.0) THEN 651 : ICE=L 652 : GOTO 226 653 : ENDIF 654 : ENDDO 655 : IF (ICE.EQ.0) THEN 656 : ICE=ITARE1.LACT(/1)+1 657 : N1=ITARE1.LACT(/1)+10 658 : N2=ITARE1.LACT(/2) 659 : SEGADJ,ITARE1 660 : ENDIF 661 : 226 ITARE1.LACT(ICE,NLA)=I 662 : C Test d'arret de la boucle (est on revenu sur le 663 : C tetraedre de depart ?) 664 : 227 IF (NT2.EQ.NT0) THEN 665 : C Remplissage de la table ITCEL1 des numeros des 666 : C cellules voisines de la cellule I 667 : IF (NBELEF1.GE.3) THEN 668 : NCVI=NCVI+1 669 : IF (NCVI.GT.ITCEL1.LACT(/1)) THEN 670 : N1=ITCEL1.LACT(/1)+10 671 : N2=ITCEL1.LACT(/2) 672 : SEGADJ,ITCEL1 673 : ENDIF 674 : ITCEL1.LACT(NCVI,I)=NVI 675 : IF (IJOK.EQ.1) THEN 676 : C ajustement du maillage de la face (I,J) 677 : NBNN=IPTF1.NUM(/1) 678 : NBELEM=NBELEF1 679 : SEGADJ,IPTF1 680 : C et de la liste de ses aretes 681 : JG=NBELEF1 682 : SEGADJ,LA1 683 : ENDIF 684 : ENDIF 685 : GOTO 228 686 : ENDIF 687 : C Nouveau points A, B, C et D pour le prochain tetraedre 688 : DO L=1,IPT2.NUM(/1) 689 : IF (IPT2.NUM(L,NT2).EQ.NREFA) NA=L 690 : IF (IPT2.NUM(L,NT2).EQ.NREFB) NB=L 691 : IF (IPT2.NUM(L,NT2).EQ.NREFD) NC=L 692 : ENDDO 693 : ND=10-NA-NB-NC 694 : NREFC=IPT2.NUM(NC,NT2) 695 : NREFD=IPT2.NUM(ND,NT2) 696 : NT1=NT2 697 : NRCC1=LCCC1.LECT(NT1) 698 : ENDDO 699 : 228 CONTINUE 700 : ENDDO 701 : 229 CONTINUE 702 : C Ajustement du maillage de la cellule I, de la liste de ses 703 : C faces et de ses voisins 704 : NBNN=IPTI1.NUM(/1) 705 : NBELEM=NBELEI1 706 : SEGADJ,IPTI1 707 : JG=NFI1 708 : SEGADJ,LFI1 709 : NBNN=IPTVI1.NUM(/1) 710 : NBELEM=IV1 711 : SEGADJ,IPTVI1 712 : ENDDO 713 : SEGSUP,LAR1,LAR2,ITFAC1,ITFAC2 714 : ENDIF 715 : C 716 : C 2.3) Ajustement final de IPTG et ITARE1 717 : NBNN=IPTG.NUM(/1) 718 : NBELEM=NBARETES 719 : SEGADJ,IPTG 720 : N1=ITARE1.LACT(/1) 721 : N2=NBARETES 722 : SEGADJ,ITARE1 723 : C 724 : C 725 : C 3) Menage avant de quitter 726 : SEGSUP,LCCC1,ITL1,ITCEL1,ITCEL2,ILEA1,ITRC1 727 : 999 RETURN 728 : END 729 : 730 :
© Cast3M 2003 - All rights reserved.
Disclaimer