1 : C VORO2 SOURCE PASCAL 14/02/17 21:15:10 7965 2 : C 3 : SUBROUTINE VORO2(IPT1,IPT2,DDMAX,MCHPOI,MTAB1) 4 : C----------------------------------------------------------------------C 5 : C Decoupage de la partition de Voronoi selon C 6 : C un contour (2D) ou une enveloppe (3D) C 7 : C C 8 : C IPT1 = maillage de POI1, centres des cellules C 9 : C IPT2 = maillage contour/enveloppe limitant la partition de C 10 : C Voronoi C 11 : C XZERO = Zero relatif a la taille de la partition a creer. C 12 : C Utile pour considerer deux sommets de cellule de Voronoi C 13 : C confondus afin de supprimer les aretes trop petites C 14 : C MTAB1 = table contenant le maillage des cellules de Voronoi non C 15 : C coupees, obtenue avec VORO1. C 16 : C C 17 : C MTAB1 est modifiee pour correspondre a la nouvelle C 18 : C partition de Voronoi restreinte par IPT2 C 19 : C C 20 : C L'algorithme utilise pour la decoupe est issu de : C 21 : C Yan D-M et al. Efficient computation of clipped Voronoi diagram for C 22 : C mesh generation, Computer-Aided Design (2011) C 23 : C C 24 : C----------------------------------------------------------------------C 25 : IMPLICIT INTEGER(I-N) 26 : IMPLICIT REAL*8 (A-H,O-Z) 27 : -INC SMCHPOI 28 : -INC CCOPTIO 29 : -INC SMELEME 30 : -INC SMCOORD 31 : -INC CCGEOME 32 : -INC SMLENTI 33 : -INC SMLREEL 34 : -INC SMTABLE 35 : LOGICAL BOOL1,BOOL2,BOOL3,BOOL4,BOOL5,BOOL6,BOOL7,BOOL8,BOOL9 36 : C Listes d'entiers utilisees un peu partout 37 : POINTEUR LV1.MLENTI,LV2.MLENTI,LV3.MLENTI,LV4.MLENTI,LV5.MLENTI 38 : POINTEUR LV6.MLENTI 39 : POINTEUR LC.MLENTI,LT.MLENTI 40 : C Tables, maillages et listes pour la table de sortie 41 : POINTEUR MTAB11.MTABLE,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 : POINTEUR MTABI2.MTABLE,IPTI2.MELEME,IPTVI2.MELEME 45 : C Tables, maillages et listes intermediaires 46 : POINTEUR MTAB10.MTABLE,IPTG0.MELEME,MTABCC.MTABLE,MTABI0.MTABLE 47 : POINTEUR IPTI0.MELEME,LFI0.MLENTI,IPTVI0.MELEME,MTABFF.MTABLE 48 : POINTEUR MTABF0.MTABLE,IPTF0.MELEME,LA0.MLENTI,MTABAA.MTABLE 49 : POINTEUR IPTA0.MELEME 50 : POINTEUR LA2.MLENTI,LA3.MLENTI 51 : C Tables, maillages et listes pour la table d'entree 52 : POINTEUR MTABCCC.MTABLE,MTABI00.MTABLE,MTABJ00.MTABLE 53 : POINTEUR LFI00.MLENTI,LFJ00.MLENTI,MTABFFF.MTABLE 54 : POINTEUR MTABF100.MTABLE,IPTVI00.MELEME 55 : POINTEUR LA00.MLENTI,LA000.MLENTI,MTABAAA.MTABLE 56 : C Liste de flottants (pour les coordonnees barycentriques des 57 : C points internes) 58 : POINTEUR LZ1.MLREEL 59 : C Tableau destine a decrire les elements adjacents de l'enveloppe 60 : SEGMENT,MADJACEN 61 : INTEGER LEAC(NBL1,IDIM,2) 62 : ENDSEGMENT 63 : POINTEUR ILEA1.MADJACEN 64 : C Tableau d'entiers a deux dimensions 65 : SEGMENT,MVORO 66 : INTEGER LACT(N1,N2) 67 : ENDSEGMENT 68 : POINTEUR ITL1.MVORO,ITARE1.MVORO,ITFAC1.MVORO 69 : C 70 : C Tableau de flottants a deux dimensions 71 : SEGMENT,MVORO2 72 : REAL*8 XACT(N1,N2) 73 : ENDSEGMENT 74 : POINTEUR LTL1.MVORO2,LTL2.MVORO2 75 : C 76 : C 77 : C 78 : IDIMP1=IDIM+1 79 : NBSOUS=0 80 : NBREF=0 81 : C 82 : C ZERO & ZERO2 : tolerances pour eliminer noeuds confondus, 83 : C aretes coplanaires... 84 : XTOL = DDMAX*1.D-10 85 : ZERO=1D-14*DDMAX 86 : ZERO2=1D-10*DDMAX 87 : C 88 : C----------------------------------------------------------------------C 89 : C PARTIE 1 : Pre-traitement C 90 : C - creation d'une table des elements adjacents du maillage C 91 : C contour/enveloppe IPT2 --> ILEA1 C 92 : C----------------------------------------------------------------------C 93 : C Le tableau ILEA1.LEAC comporte 3 dimensions : ILEA1.LEAC(I,J,K) 94 : C dimension 1 : I est le numero de l'element du contour/enveloppe 95 : C dimension 2 : J est le numero du noeud de l'element I 96 : C dimension 3 : K=1 --> renvoie le numero de l'element adjacent a 97 : C l'element I par rapport au noeud J 98 : C K=2 --> renvoie le numero du noeud de l'element 99 : C adjacent situe en vis a vis du noeud J 100 : C On a donc les symetries suivantes : 101 : C si ILEA1.LEAC(I,J,1) = I' et ILEA1.LEAC(I,J,2) = J' 102 : C alors ILEA1.LEAC(I',J',1) = I et ILEA1.LEAC(I',J',2) = J 103 : C 104 : C Initialisation de la table des elements adjacents 105 : NBL1=IPT2.NUM(/2) 106 : NBL2=IPT2.NUM(/1) 107 : SEGINI,ILEA1 108 : JG=IPT2.NUM(/1) 109 : SEGINI,LV1 110 : C Nombre de noeuds en commun a trouver 111 : NNREF=IPT2.NUM(/1)-1 112 : C Somme des numeros des noeuds 113 : IF (IDIM.EQ.2) THEN 114 : NSREF=1+2 115 : ELSE 116 : NSREF=1+2+3 117 : ENDIF 118 : C Boucle sur les elements de IPT2 119 : DO I=1,IPT2.NUM(/2) 120 : C Numeros des noeuds de l'element I 121 : DO J=1,IPT2.NUM(/1) 122 : LV1.LECT(J)=IPT2.NUM(J,I) 123 : ENDDO 124 : C Deuxieme boucle sur les elements de numero superieurs a I 125 : DO II=I+1,IPT2.NUM(/2) 126 : NNC=0 127 : NSOMA=0 128 : NSOMB=0 129 : C Test si l'element II a des noeuds communs a l'element I 130 : DO J=1,IPT2.NUM(/1) 131 : NTEST=IPT2.NUM(J,II) 132 : DO K=1,IPT2.NUM(/1) 133 : IF (NTEST.EQ.LV1.LECT(K)) THEN 134 : NNC=NNC+1 135 : NSOMA=NSOMA+K 136 : NSOMB=NSOMB+J 137 : ENDIF 138 : ENDDO 139 : ENDDO 140 : C Si l'element II est bien adjacent a l'element I 141 : IF (NNC.EQ.NNREF) THEN 142 : NI=NSREF-NSOMA 143 : NII=NSREF-NSOMB 144 : ILEA1.LEAC(I,NI,1)=II 145 : ILEA1.LEAC(I,NI,2)=NII 146 : ILEA1.LEAC(II,NII,1)=I 147 : ILEA1.LEAC(II,NII,2)=NI 148 : ENDIF 149 : ENDDO 150 : C Test si on a bien trouve un voisin pour chaque cote de I 151 : DO J=1,IPT2.NUM(/1) 152 : IF (ILEA1.LEAC(I,J,1).EQ.0) THEN 153 : C Le contour n'est pas reconnu ferme 153 : CALL ERREUR(28) 154 : GOTO 999 155 : ENDIF 156 : ENDDO 157 : ENDDO 158 : C 159 : C 160 : C----------------------------------------------------------------------C 161 : C PARTIE 2 C 162 : C Construction de la surface externe de la partition de Voronoi C 163 : C par decoupage en parcourant la surface enveloppe C 164 : C----------------------------------------------------------------------C 165 : C 166 : IF (IDIM.EQ.2) THEN 167 : C Idee : on parcoure les segments de l'enveloppe et pour chacun 168 : C d'eux on determine son intersection avec la partition de 169 : C Voronoi 170 : C On construit ainsi la liste des couples (Ti;Ci) possedant une 171 : C intersection non vide 172 : C Ti : le numero de segment dans IPT2 173 : C Ci : le numero de la cellule dans IPT1 174 : C La liste des couples (Ti;Ci) a traiter est incrementee au fur 175 : C et a mesure des rencontres 176 : C 177 : C------- Etape 1 : Initialisations 178 : C Initialisation d'une file pour les Ci et Tj 179 : C on utilise deux listes d'entiers LT et LC 180 : JG=(IPT2.NUM(/2))*(IPT1.NUM(/2)) 181 : SEGINI,LT,LC 182 : C On demarre sur le premier segment de l'enveloppe 183 : NLT=1 184 : NT0=1 185 : LT.LECT(NLT)=NT0 186 : C Recherche d'une cellule intersectant le segment NT0, on prend 187 : C celle dont le centre est le plus proche du barycentre du 188 : C segment NT0 189 : NRA=IPT2.NUM(1,NT0) 190 : NRB=IPT2.NUM(2,NT0) 191 : XG=((XCOOR((NRA-1)*IDIMP1+1))+(XCOOR((NRB-1)*IDIMP1+1)))/2. 192 : YG=((XCOOR((NRA-1)*IDIMP1+2))+(XCOOR((NRB-1)*IDIMP1+2)))/2. 193 : DO NCJ=1,IPT1.NUM(/2) 194 : NRCJ=IPT1.NUM(1,NCJ) 195 : XCJ=XCOOR((NRCJ-1)*IDIMP1+1) 196 : YCJ=XCOOR((NRCJ-1)*IDIMP1+2) 197 : DJG=SQRT(((XCJ-XG)*(XCJ-XG))+((YCJ-YG)*(YCJ-YG))) 198 : IF (NCJ.EQ.1) THEN 199 : DJG0=DJG 200 : NCJ0=NCJ 201 : NRCJ0=NRCJ 202 : ENDIF 203 : IF (DJG.LT.DJG0) THEN 204 : DJG0=DJG 205 : NCJ0=NCJ 206 : NRCJ0=NRCJ 207 : ENDIF 208 : ENDDO 209 : LC.LECT(NLT)=NCJ0 210 : C Initialisation d'une table pour les couples (Ti,Ci) 211 : N1=IPT2.NUM(/2) 212 : N2=IPT1.NUM(/2) 213 : SEGINI,ITL1 214 : ITL1.LACT(NT0,NCJ0)=1 215 : C Maillage global IPTG des aretes de la partition coupee 216 : NBNN=2 217 : NBELEM=200 218 : SEGINI,IPTG 219 : IPTG.ITYPEL=2 220 : NBELEG=0 221 : C Liste de tous les points ajoutes a MCOORD suite au decoupage 222 : JG=0 223 : SEGINI,LV1 224 : NLV1=0 225 : C Maillage global IPT6 des centres des cellules dans la 226 : C partition coupee 227 : NBNN=1 228 : NBELEM=IPT1.NUM(/2) 229 : SEGINI,IPT6 230 : IPT6.ITYPEL=1 231 : NBELE6=0 232 : C Table de sortie MTAB11 232 : CALL CRTABL(MTAB11) 232 : CALL CRTABL(MTABC) 232 : CALL CRTABL(MTABA) 232 : CALL ECCTAB(MTAB11,'MOT ',0,0,'VISU',0,0, 233 : & 'MAILLAGE',0,0,0,0,IPTG) 233 : CALL ECCTAB(MTAB11,'MOT ',0,0,'CELL',0,0, 234 : & 'TABLE ',0,0,0,0,MTABC) 234 : CALL ECCTAB(MTAB11,'MOT ',0,0,'ARTS',0,0, 235 : & 'TABLE ',0,0,0,0,MTABA) 236 : C 237 : C------- Etape 2 : Calcul des intersections 238 : C Maintenant que la file (LT,LC) est non vide, on debute 239 : C les iterations. Pour chaque couple, on determine le segment 240 : C d'intersection et on le range dans la table de sortie. 241 : C On incremente egalement la file avec les cellules voisines et 242 : C les segments voisins du segment d'intersection 243 : C Les iterations s'arretent lorsque la file est vide 244 : C 245 : C Tables des cellules et des aretes non coupees 245 : CALL ACCTAB(MTAB1,'MOT ',0,0,'CELL',0,0, 246 : & 'TABLE ',0,0,0,0,MTABCCC) 246 : CALL ACCTAB(MTAB1,'MOT ',0,0,'ARTS',0,0, 247 : & 'TABLE ',0,0,0,0,MTABAAA) 248 : C Debut de la boucle sur la file (LT,LC) 249 : NB1=IPT1.NUM(/2)*IPT2.NUM(/2) 250 : DO I=1,NB1 251 : C Test d'arret de la boucle (la file est vide) 252 : IF (I.GT.NLT) THEN 253 : GOTO 35 254 : ENDIF 255 : C NTI designe le numero de segment dans IPT2 256 : NTI=LT.LECT(I) 257 : C NCI designe le numero la cellule dans IPT1 258 : NCI=LC.LECT(I) 259 : C NRCI designe le numero du point centre de la cellule NCI 260 : NRCI=IPT1.NUM(1,NCI) 261 : C Creation d'un indice NRCI dans la table MTAB11.CELL 262 : C si besoin 262 : CALL ECROBJ('POINT ',NRCI) 262 : CALL ECROBJ('TABLE ',MTABC) 262 : CALL EXIS 262 : CALL LIRLOG(BOOL1,1,IRETOU) 263 : IF (BOOL1) THEN 264 : C Si MTAB11.CELL.NRCI existe, on va le chercher 264 : CALL ACCTAB(MTABC,'POINT ',0,0,' ',0,NRCI, 265 : & 'TABLE ',0,0,0,0,MTABI1) 266 : C Recuperation des indices VISU, ARTS et VOIS 266 : CALL ACCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 267 : & 'MAILLAGE',0,0,0,0,IPTI1) 267 : CALL ACCTAB(MTABI1,'MOT ',0,0,'ARTS',0,0, 268 : & 'LISTENTI',0,0,0,0,LA1) 268 : CALL ACCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, 269 : & 'MAILLAGE',0,0,0,0,IPTVI1) 270 : ELSE 271 : C Sinon, on cree un nouvel indice 271 : CALL CRTABL(MTABI1) 271 : CALL ECCTAB(MTABC,'POINT ',0,0,' ',0,NRCI, 272 : & 'TABLE ',0,0,0,0,MTABI1) 273 : C Initialisation des indices VISU, ARTS et VOIS 274 : NBNN=2 275 : NBELEM=0 276 : SEGINI,IPTI1 277 : IPTI1.ITYPEL=2 277 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 278 : & 'MAILLAGE',0,0,0,0,IPTI1) 279 : JG=0 280 : SEGINI,LA1 280 : CALL ECCTAB(MTABI1,'MOT ',0,0,'ARTS',0,0, 281 : & 'LISTENTI',0,0,0,0,LA1) 282 : NBNN=1 283 : NBELEM=0 284 : SEGINI,IPTVI1 285 : IPTVI1.ITYPEL=1 285 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, 286 : & 'MAILLAGE',0,0,0,0,IPTVI1) 287 : NBELE6=NBELE6+1 288 : IPT6.NUM(1,NBELE6)=NRCI 289 : ENDIF 290 : C Determination de l'intersection entre le segment NTI et 291 : C le polygone NCI 292 : C Clipping selon l'algorithme de Sutherland et Hodgman 293 : C Le segment resultant de l'intersection est decrit par la 294 : C table LTL1 avec : 295 : C LTL1.XACT(i,j) = coordonnee j du i-eme noeud du segment 296 : C i valant 1 ou 2 297 : C --> initialisation de LTL1 aux noeuds du segment NTI 298 : NRTI1=IPT2.NUM(1,NTI) 299 : NRTI2=IPT2.NUM(2,NTI) 300 : N1=2 301 : N2=3 302 : SEGINI,LTL1 303 : XI10=XCOOR((NRTI1-1)*IDIMP1+1) 304 : YI10=XCOOR((NRTI1-1)*IDIMP1+2) 305 : XI20=XCOOR((NRTI2-1)*IDIMP1+1) 306 : YI20=XCOOR((NRTI2-1)*IDIMP1+2) 307 : LTL1.XACT(1,1)=XI10 308 : LTL1.XACT(1,2)=YI10 309 : LTL1.XACT(2,1)=XI20 310 : LTL1.XACT(2,2)=YI20 311 : C --> table du polygone NRCI infini (MTABI0) 311 : CALL ACCTAB(MTABCCC,'POINT ',0,0,' ',0,NRCI, 312 : & 'TABLE ',0,0,0,0,MTABI00) 313 : C --> liste de ses aretes 313 : CALL ACCTAB(MTABI00,'MOT ',0,0,'ARTS',0,0, 314 : & 'LISTENTI',0,0,0,0,LA00) 315 : C --> maillage de ses voisins 315 : CALL ACCTAB(MTABI00,'MOT ',0,0,'VOIS',0,0, 316 : & 'MAILLAGE',0,0,0,0,IPT3) 317 : C Boucle sur toutes les aretes du polygone NCI 318 : DO J=1,LA00.LECT(/1) 319 : C Maillage IPT4 de l'arete J du polygone NCI 320 : NFJ=LA00.LECT(J) 320 : CALL ACCTAB(MTABAAA,'ENTIER ',NFJ,0,' ',0,0, 321 : & 'MAILLAGE',0,0,0,0,IPT4) 322 : NREFA=IPT4.NUM(1,1) 323 : NREFB=IPT4.NUM(2,1) 324 : C Massicotage du segment LTL1 selon le demi plan definit 325 : C la ligne IPT4 et le centre I du polygone de Voronoi 326 : C Coordonnees des points P1 et P2 327 : X1=LTL1.XACT(1,1) 328 : Y1=LTL1.XACT(1,2) 329 : X2=LTL1.XACT(2,1) 330 : Y2=LTL1.XACT(2,2) 331 : C Intersection du segment P1P2 avec le demi plan 331 : CALL INTPSEG2(NREFA,NREFB,NRCI,X1,Y1,X2,Y2, 332 : & IINT,XI1,YI1,XI2,YI2) 333 : C Il y a 2 cas possibles 334 : C -1) l'intersection est vide --> erreur 335 : IF (IINT.EQ.0) THEN 336 : C WRITE(6,*) '> Intersection segment/polygone vide' 337 : C WRITE(6,*) 'Numero du segment concerne:',nti 338 : C WRITE(6,*) 'Point de la cellule concernee',nrci 339 : INTERR(1)=nti 340 : INTERR(2)=nrci 340 : CALL ERREUR(1032) 341 : GOTO 999 342 : ENDIF 343 : C -2) l'intersection n'est pas vide 344 : IF (IINT.EQ.1) THEN 345 : LTL1.XACT(1,1)=XI1 346 : LTL1.XACT(1,2)=YI1 347 : LTL1.XACT(2,1)=XI2 348 : LTL1.XACT(2,2)=YI2 349 : ENDIF 350 : ENDDO 351 : C 352 : C 353 : C Maintenant, on va incrementer la file (LT,LC) 354 : C avec les intersections voisines 355 : C 356 : C Extremite 1 du segment d'intersection LTL1 357 : XI1=LTL1.XACT(1,1) 358 : YI1=LTL1.XACT(1,2) 359 : C 1er cas : P1 correspond l'extremite 1 du segment NTI 360 : C donc le segment adjacent a NTI coupe aussi la 361 : C cellule NCI 362 : IF ((XI1.EQ.XI10).AND.(YI1.EQ.YI10)) THEN 363 : NTVI=ILEA1.LEAC(NTI,2,1) 364 : C On a trouve le segment NTVI coupant la cellule NCI 365 : C On verifie que le couple (NTVI,NCI) n'a pas deja 366 : C ete traite et on l'ajoute a la file 367 : IF (ITL1.LACT(NTVI,NCI).EQ.0) THEN 368 : ITL1.LACT(NTVI,NCI)=1 369 : NLT=NLT+1 370 : LT.LECT(NLT)=NTVI 371 : LC.LECT(NLT)=NCI 372 : ENDIF 373 : C 2eme cas : P1 est un point d'intersection 374 : C donc le polygone voisin coupe le segment NTI 375 : ELSE 376 : C Recherche du polygone voisin en question 377 : C Boucle sur voisins du polygone NCI (on ne s'interesse 378 : C qu'aux aretes internes) 379 : DO J=1,IPT3.NUM(/2) 380 : NRCJ=IPT3.NUM(1,J) 381 : C Liste des aretes du polyedre NRCJ 381 : CALL ACCTAB(MTABCCC,'POINT ',0,0,' ',0,NRCJ, 382 : & 'TABLE ',0,0,0,0,MTABJ00) 382 : CALL ACCTAB(MTABJ00,'MOT ',0,0,'ARTS',0,0, 383 : & 'LISTENTI',0,0,0,0,LA000) 384 : C Recherche du numero de l'arete entre NRCI et NRCJ 385 : C cette arete est a la fois dans LA00 et LA000 386 : NAJ=0 387 : DO K=1,LA00.LECT(/1) 388 : NAJ=LA00.LECT(K) 389 : DO L=1,LA000.LECT(/1) 390 : IF (LA000.LECT(L).EQ.NAJ) GOTO 27 391 : ENDDO 392 : ENDDO 393 : 27 CONTINUE 394 : C Maillage IPT5 de l'arete J 394 : CALL ACCTAB(MTABAAA,'ENTIER ',NAJ,0,' ',0,0, 395 : & 'MAILLAGE',0,0,0,0,IPT5) 396 : C Points AB definissant l'arete IPT5 397 : NRA=IPT5.NUM(1,1) 398 : NRB=IPT5.NUM(2,1) 399 : XA=XCOOR((NRA-1)*IDIMP1+1) 400 : YA=XCOOR((NRA-1)*IDIMP1+2) 401 : XB=XCOOR((NRB-1)*IDIMP1+1) 402 : YB=XCOOR((NRB-1)*IDIMP1+2) 403 : C Test si P1 est situe entre A et B 404 : DA1=SQRT(((XA-XI1)**2)+((YA-YI1)**2)) 405 : DB1=SQRT(((XB-XI1)**2)+((YB-YI1)**2)) 406 : DAB=SQRT(((XA-XB)**2)+((YA-YB)**2)) 407 : DTEST=(DA1+DB1-DAB)/DAB 408 : IF (DTEST.LT.ZERO2) THEN 409 : C NRCJ est bien le voisin de NRCI qui coupe aussi NTI 410 : C Identification de la cellule associee au point NRCJ 411 : NCJ=0 412 : DO K=1,IPT1.NUM(/2) 413 : NRTEST=IPT1.NUM(1,K) 414 : IF (NRTEST.EQ.NRCJ) THEN 415 : NCJ=K 416 : GOTO 28 417 : ENDIF 418 : ENDDO 419 : ENDIF 420 : ENDDO 421 : C Si on passe ici c'est que l'on a pas trouve la cellule 422 : C voisine partageant l'arete 423 : C WRITE(6,*) ' > Cellule voisine non trouvee' 424 : INTERR(1)=NRCI 425 : INTERR(2)=NAJ 425 : CALL ERREUR(1034) 426 : GOTO 999 427 : 28 CONTINUE 428 : C On a trouve la cellule NCJ coupant le segment NTI 429 : C On verifie que le couple (NTI,NCJ) n'a pas deja 430 : C ete traite et on l'ajoute a la file 431 : IF (ITL1.LACT(NTI,NCJ).EQ.0) THEN 432 : ITL1.LACT(NTI,NCJ)=1 433 : NLT=NLT+1 434 : LT.LECT(NLT)=NTI 435 : LC.LECT(NLT)=NCJ 436 : ENDIF 437 : ENDIF 438 : C 439 : C Extremite 2 du segment d'intersection LTL1 440 : XI2=LTL1.XACT(2,1) 441 : YI2=LTL1.XACT(2,2) 442 : C 1er cas : P2 correspond l'extremite 2 du segment NTI 443 : C donc le segment adjacent a NTI coupe aussi la 444 : C cellule NCI 445 : IF ((XI2.EQ.XI20).AND.(YI2.EQ.YI20)) THEN 446 : NTVI=ILEA1.LEAC(NTI,1,1) 447 : C On a trouve le segment NTVI coupant la cellule NCI 448 : C On verifie que le couple (NTVI,NCI) n'a pas deja 449 : C ete traite et on l'ajoute a la file 450 : IF (ITL1.LACT(NTVI,NCI).EQ.0) THEN 451 : ITL1.LACT(NTVI,NCI)=1 452 : NLT=NLT+1 453 : LT.LECT(NLT)=NTVI 454 : LC.LECT(NLT)=NCI 455 : ENDIF 456 : C 2eme cas : P2 est un point d'intersection 457 : C donc le polygone voisin coupe le segment NTI 458 : ELSE 459 : C Recherche du polygone voisin en question 460 : C Boucle sur voisins du polygone NCI (on ne s'interesse 461 : C qu'aux aretes internes) 462 : DO J=1,IPT3.NUM(/2) 463 : NRCJ=IPT3.NUM(1,J) 464 : C Liste des aretes du polyedre NRCJ 464 : CALL ACCTAB(MTABCCC,'POINT ',0,0,' ',0,NRCJ, 465 : & 'TABLE ',0,0,0,0,MTABJ00) 465 : CALL ACCTAB(MTABJ00,'MOT ',0,0,'ARTS',0,0, 466 : & 'LISTENTI',0,0,0,0,LA000) 467 : C Recherche du numero de l'arete entre NRCI et NRCJ 468 : C cette arete est a la fois dans LA00 et LA000 469 : NAJ=0 470 : DO K=1,LA00.LECT(/1) 471 : NAJ=LA00.LECT(K) 472 : DO L=1,LA000.LECT(/1) 473 : IF (LA000.LECT(L).EQ.NAJ) GOTO 29 474 : ENDDO 475 : ENDDO 476 : 29 CONTINUE 477 : C Maillage IPT5 de l'arete J du polygone NCI 477 : CALL ACCTAB(MTABAAA,'ENTIER ',NAJ,0,' ',0,0, 478 : & 'MAILLAGE',0,0,0,0,IPT5) 479 : C Points AB definissant l'arete IPT5 480 : NRA=IPT5.NUM(1,1) 481 : NRB=IPT5.NUM(2,1) 482 : XA=XCOOR((NRA-1)*IDIMP1+1) 483 : YA=XCOOR((NRA-1)*IDIMP1+2) 484 : XB=XCOOR((NRB-1)*IDIMP1+1) 485 : YB=XCOOR((NRB-1)*IDIMP1+2) 486 : C Test si P2 est situe entre A et B 487 : DA2=SQRT(((XA-XI2)**2)+((YA-YI2)**2)) 488 : DB2=SQRT(((XB-XI2)**2)+((YB-YI2)**2)) 489 : DAB=SQRT(((XA-XB)**2)+((YA-YB)**2)) 490 : DTEST=(DA2+DB2-DAB)/DAB 491 : IF (DTEST.LT.ZERO2) THEN 492 : C NRCJ est bien le voisin de NRCI qui coupe aussi NTI 493 : C Identification de la cellule associee au point NRCJ 494 : NCJ=0 495 : DO K=1,IPT1.NUM(/2) 496 : NRTEST=IPT1.NUM(1,K) 497 : IF (NRTEST.EQ.NRCJ) THEN 498 : NCJ=K 499 : GOTO 30 500 : ENDIF 501 : ENDDO 502 : ENDIF 503 : ENDDO 504 : C Si on passe ici c'est que l'on a pas trouve la cellule 505 : C voisine partageant l'arete 506 : C WRITE(6,*) ' > Cellule voisine non trouvee' 507 : INTERR(1)=NRCI 508 : INTERR(2)=NAJ 508 : CALL ERREUR(1034) 509 : GOTO 999 510 : 30 CONTINUE 511 : C On a trouve la cellule NCJ coupant le segment NTI 512 : C On verifie que le couple (NTI,NCJ) n'a pas deja 513 : C ete traite et on l'ajoute a la file 514 : IF (ITL1.LACT(NTI,NCJ).EQ.0) THEN 515 : ITL1.LACT(NTI,NCJ)=1 516 : NLT=NLT+1 517 : LT.LECT(NLT)=NTI 518 : LC.LECT(NLT)=NCJ 519 : ENDIF 520 : ENDIF 521 : C 522 : C Determination des noeuds P1 et P2 de l'arete 523 : C Doit on creer de nouveaux noeuds ? 524 : C On examine la liste LV1 des noeuds crees 525 : C test si P1 est deja present ou non dans LV1 526 : BOOL2=.FALSE. 527 : DO L=1,LV1.LECT(/1) 528 : NRL=LV1.LECT(L) 529 : XL=XCOOR((NRL-1)*IDIMP1+1) 530 : YL=XCOOR((NRL-1)*IDIMP1+2) 531 : IF ((ABS(XL-XI1).LT.XTOL).AND. 532 : & (ABS(YL-YI1).LT.XTOL)) THEN 533 : BOOL2=.TRUE. 534 : NRP1=NRL 535 : GOTO 31 536 : ENDIF 537 : ENDDO 538 : IF (.NOT.BOOL2) THEN 539 : C ajout d'un nouveau noeud dans MCOORD 540 : NRP1=(XCOOR(/1)/IDIMP1)+1 541 : NBPTS=NRP1 542 : SEGADJ,MCOORD 543 : XCOOR((NRP1-1)*IDIMP1+1)=XI1 544 : XCOOR((NRP1-1)*IDIMP1+2)=YI1 545 : NLV1=NLV1+1 546 : IF (NLV1.GT.LV1.LECT(/1)) THEN 547 : JG=NLV1 548 : SEGADJ,LV1 549 : ENDIF 550 : LV1.LECT(NLV1)=NRP1 551 : ENDIF 552 : 31 CONTINUE 553 : C test si P2 est deja present ou non dans LV1 554 : BOOL2=.FALSE. 555 : DO L=1,LV1.LECT(/1) 556 : NRL=LV1.LECT(L) 557 : XL=XCOOR((NRL-1)*IDIMP1+1) 558 : YL=XCOOR((NRL-1)*IDIMP1+2) 559 : IF ((ABS(XL-XI2).LT.XTOL).AND. 560 : & (ABS(YL-YI2).LT.XTOL)) THEN 561 : BOOL2=.TRUE. 562 : NRP2=NRL 563 : GOTO 32 564 : ENDIF 565 : ENDDO 566 : IF (.NOT.BOOL2) THEN 567 : C ajout d'un nouveau noeud dans MCOORD 568 : NRP2=(XCOOR(/1)/IDIMP1)+1 569 : NBPTS=NRP2 570 : SEGADJ,MCOORD 571 : XCOOR((NRP2-1)*IDIMP1+1)=XI2 572 : XCOOR((NRP2-1)*IDIMP1+2)=YI2 573 : NLV1=NLV1+1 574 : IF (NLV1.GT.LV1.LECT(/1)) THEN 575 : JG=NLV1 576 : SEGADJ,LV1 577 : ENDIF 578 : LV1.LECT(NLV1)=NRP2 579 : ENDIF 580 : 32 CONTINUE 581 : C Test sur l'existence de l'arete P1 P2 dans IPTG 582 : IAE=0 583 : DO L=1,NBELEG 584 : NG1=IPTG.NUM(1,L) 585 : NG2=IPTG.NUM(2,L) 586 : IF (((NG1.EQ.NRP1).AND.(NG2.EQ.NRP2)).OR. 587 : & ((NG2.EQ.NRP1).AND.(NG1.EQ.NRP2))) THEN 588 : C si l'arete P1 P2 a deja ete creee, on 589 : C l'indique au moyen de IAE 590 : IAE=1 591 : NAE=L 592 : GOTO 33 593 : ENDIF 594 : ENDDO 595 : 33 CONTINUE 596 : C --> Ajout de l'arete dans MTAB11.VISU 597 : IF (IAE.EQ.0) THEN 598 : NBELEG=NBELEG+1 599 : IPTG.NUM(1,NBELEG)=NRP1 600 : IPTG.NUM(2,NBELEG)=NRP2 601 : IF (NBELEG.EQ.IPTG.NUM(/2)) THEN 602 : NBNN=IPTG.NUM(/1) 603 : NBELEM=IPTG.NUM(/2)+100 604 : SEGADJ,IPTG 605 : ENDIF 606 : NAE=NBELEG 607 : ENDIF 608 : C --> Ajout de l'arete dans MTAB11.ARTS 609 : IF (IAE.EQ.0) THEN 610 : NBNN=2 611 : NBELEM=1 612 : SEGINI,IPTA1 613 : IPTA1.ITYPEL=2 614 : IPTA1.NUM(1,1)=NRP1 615 : IPTA1.NUM(2,1)=NRP2 615 : CALL ECCTAB(MTABA,'ENTIER ',NAE,0,' ',0,0, 616 : & 'MAILLAGE',0,0,0,0,IPTA1) 617 : ENDIF 618 : C --> Ajout de l'arete dans MTAB11.CELL.NRCI si besoin 619 : DO L=1,LA1.LECT(/1) 620 : NAL=LA1.LECT(L) 621 : IF (NAL.EQ.NAE) GOTO 34 622 : ENDDO 623 : C dans le maillage de la cellule 624 : NBNN=IPTI1.NUM(/1) 625 : NBELEM=IPTI1.NUM(/2)+1 626 : SEGADJ,IPTI1 627 : IPTI1.NUM(1,NBELEM)=NRP1 628 : IPTI1.NUM(2,NBELEM)=NRP2 629 : C et la liste des aretes de la cellule 630 : JG=LA1.LECT(/1)+1 631 : SEGADJ,LA1 632 : LA1.LECT(JG)=NAE 633 : 34 CONTINUE 634 : ENDDO 635 : 35 CONTINUE 636 : C Ajustement du maillage global des centres des cellules dans la 637 : C partition coupee IPT6 638 : NBNN=IPT6.NUM(/1) 639 : NBELEM=NBELE6 640 : SEGADJ,IPT6 641 : C Ajustement du maillage global des aretes IPTG 642 : NBNN=IPTG.NUM(/1) 643 : NBELEM=NBELEG 644 : SEGADJ,IPTG 645 : C 646 : C 647 : C Fusion des aretes coplanaires (a faire) 648 : C 649 : C 650 : C 651 : ELSEIF (IDIM.EQ.3) THEN 652 : C Idee : on parcoure les triangles de l'enveloppe et pour chacun 653 : C d'eux on determine son intersection avec la partition de 654 : C Voronoi 655 : C On construit ainsi la liste des couples (Ti;Ci) possedant une 656 : C intersection non vide 657 : C Ti : le numero de triangle dans IPT2 658 : C Ci : le numero de la cellule dans IPT1 659 : C La liste des couples (Ti;Ci) a traiter est incrementee au fur 660 : C et a mesure des rencontres 661 : C 662 : C------- Etape 1 : Initialisations 663 : C Initialisation d'une file pour les Ci et Tj 664 : C on utilise deux listes d'entiers LT et LC 665 : JG=(IPT2.NUM(/2))*(IPT1.NUM(/2)) 666 : SEGINI,LT,LC 667 : C On demarre sur le premier triangle de l'enveloppe 668 : NLT=1 669 : NT0=1 670 : LT.LECT(NLT)=NT0 671 : C Recherche d'une cellule intersectant le triangle NT0, on prend 672 : C celle dont le centre est le plus proche du barycentre du 673 : C triangle NT0 674 : NRA=IPT2.NUM(1,NT0) 675 : NRB=IPT2.NUM(2,NT0) 676 : NRC=IPT2.NUM(3,NT0) 677 : XG=((XCOOR((NRA-1)*IDIMP1+1))+(XCOOR((NRB-1)*IDIMP1+1))+ 678 : & (XCOOR((NRC-1)*IDIMP1+1)))/3. 679 : YG=((XCOOR((NRA-1)*IDIMP1+2))+(XCOOR((NRB-1)*IDIMP1+2))+ 680 : & (XCOOR((NRC-1)*IDIMP1+2)))/3. 681 : ZG=((XCOOR((NRA-1)*IDIMP1+3))+(XCOOR((NRB-1)*IDIMP1+3))+ 682 : & (XCOOR((NRC-1)*IDIMP1+3)))/3. 683 : DO NCJ=1,IPT1.NUM(/2) 684 : NRCJ=IPT1.NUM(1,NCJ) 685 : XCJ=XCOOR((NRCJ-1)*IDIMP1+1) 686 : YCJ=XCOOR((NRCJ-1)*IDIMP1+2) 687 : ZCJ=XCOOR((NRCJ-1)*IDIMP1+3) 688 : DJG=SQRT(((XCJ-XG)*(XCJ-XG))+((YCJ-YG)*(YCJ-YG))+ 689 : & ((ZCJ-ZG)*(ZCJ-ZG))) 690 : IF (NCJ.EQ.1) THEN 691 : DJG0=DJG 692 : NCJ0=NCJ 693 : NRCJ0=NRCJ 694 : ENDIF 695 : IF (DJG.LT.DJG0) THEN 696 : DJG0=DJG 697 : NCJ0=NCJ 698 : NRCJ0=NRCJ 699 : ENDIF 700 : ENDDO 701 : LC.LECT(NLT)=NCJ0 702 : C Initialisation d'une table pour les couples (Ti,Ci) 703 : N1=IPT2.NUM(/2) 704 : N2=IPT1.NUM(/2) 705 : SEGINI,ITL1 706 : ITL1.LACT(NT0,NCJ0)=1 707 : C Maillage global IPTG des aretes de la partition coupee 708 : NBNN=2 709 : NBELEM=200 710 : SEGINI,IPTG 711 : IPTG.ITYPEL=2 712 : NBELEG=0 713 : C Liste de tous les points ajoutes a MCOORD suite au decoupage 714 : JG=0 715 : SEGINI,LV1 716 : NLV1=0 717 : C Maillage global IPT6 des centres des cellules dans la 718 : C partition coupee 719 : NBNN=1 720 : NBELEM=IPT1.NUM(/2) 721 : SEGINI,IPT6 722 : IPT6.ITYPEL=1 723 : NBELE6=0 724 : C Table de sortie MTAB11 724 : CALL CRTABL(MTAB11) 724 : CALL CRTABL(MTABC) 724 : CALL CRTABL(MTABF) 724 : CALL CRTABL(MTABA) 724 : CALL ECCTAB(MTAB11,'MOT ',0,0,'VISU',0,0, 725 : & 'MAILLAGE',0,0,0,0,IPTG) 725 : CALL ECCTAB(MTAB11,'MOT ',0,0,'CELL',0,0, 726 : & 'TABLE ',0,0,0,0,MTABC) 726 : CALL ECCTAB(MTAB11,'MOT ',0,0,'FACS',0,0, 727 : & 'TABLE ',0,0,0,0,MTABF) 727 : CALL ECCTAB(MTAB11,'MOT ',0,0,'ARTS',0,0, 728 : & 'TABLE ',0,0,0,0,MTABA) 729 : C 730 : C------- Etape 2 : Calcul des intersections 731 : C Maintenant que la file (LT,LC) est non vide, on debute 732 : C les iterations. Pour chaque couple, on determine le polygone 733 : C d'intersection et on le range dans la table de sortie. 734 : C On incremente egalement la file avec les cellules voisines et 735 : C les triangles voisins du polygone d'intersection 736 : C Les iterations s'arretent lorsque la file est vide 737 : C 738 : C Tables des cellules et des faces non coupees 738 : CALL ACCTAB(MTAB1,'MOT ',0,0,'CELL',0,0, 739 : & 'TABLE ',0,0,0,0,MTABCCC) 739 : CALL ACCTAB(MTAB1,'MOT ',0,0,'FACS',0,0, 740 : & 'TABLE ',0,0,0,0,MTABFFF) 741 : C 742 : C NF : nombre de faces total dans la partition coupee 743 : NF=0 744 : C ITFAC1 : donne le numero de face pour chaque couple de cellules 745 : N1=IPT1.NUM(/2) 746 : N2=IPT1.NUM(/2) 747 : SEGINI,ITFAC1 748 : C On construit LV5 et LV6 qui donnent le couple de cellule 749 : C associe a chaque numero de face (inverse de ITFAC1) 750 : C LV5(f1)=ci et LV6(f1)=cj --> la face f1 est interne et separe 751 : C les cellules ci et cj 752 : C LV5(f1)=0 et LV6(f1)=0 --> la face f1 est externe 753 : JG=2*IPT2.NUM(/2) 754 : SEGINI,LV5,LV6 755 : C Debut de la boucle sur la file (LT,LC) 756 : NB1=IPT1.NUM(/2)*IPT2.NUM(/2) 757 : DO I=1,NB1 758 : C Test d'arret de la boucle (la file est vide) 759 : IF (I.GT.NLT) THEN 760 : GOTO 1 761 : ENDIF 762 : NF=NF+1 763 : C NTI designe le numero de triangle dans IPT2 764 : NTI=LT.LECT(I) 765 : C NCI designe le numero la cellule dans IPT1 766 : NCI=LC.LECT(I) 767 : C NRCI designe le numero du point centre de la cellule NCI 768 : NRCI=IPT1.NUM(1,NCI) 769 : C Creation d'un indice NRCI dans la table MTAB11.CELL 770 : C si besoin 770 : CALL ECROBJ('POINT ',NRCI) 770 : CALL ECROBJ('TABLE ',MTABC) 770 : CALL EXIS 770 : CALL LIRLOG(BOOL1,1,IRETOU) 771 : IF (BOOL1) THEN 772 : C Si MTAB11.CELL.NRCI existe, on va le chercher 772 : CALL ACCTAB(MTABC,'POINT ',0,0,' ',0,NRCI, 773 : & 'TABLE ',0,0,0,0,MTABI1) 774 : C Recuperation des indices VISU, FACS et VOIS 774 : CALL ACCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 775 : & 'MAILLAGE',0,0,0,0,IPTI1) 775 : CALL ACCTAB(MTABI1,'MOT ',0,0,'FACS',0,0, 776 : & 'LISTENTI',0,0,0,0,LFI1) 776 : CALL ACCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, 777 : & 'MAILLAGE',0,0,0,0,IPTVI1) 778 : ELSE 779 : C Sinon, on cree un nouvel indice 779 : CALL CRTABL(MTABI1) 779 : CALL ECCTAB(MTABC,'POINT ',0,0,' ',0,NRCI, 780 : & 'TABLE ',0,0,0,0,MTABI1) 781 : C Initialisation des indices VISU, FACS et VOIS 782 : NBNN=2 783 : NBELEM=0 784 : SEGINI,IPTI1 785 : IPTI1.ITYPEL=2 785 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 786 : & 'MAILLAGE',0,0,0,0,IPTI1) 787 : JG=0 788 : SEGINI,LFI1 788 : CALL ECCTAB(MTABI1,'MOT ',0,0,'FACS',0,0, 789 : & 'LISTENTI',0,0,0,0,LFI1) 790 : NBNN=1 791 : NBELEM=0 792 : SEGINI,IPTVI1 793 : IPTVI1.ITYPEL=1 793 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, 794 : & 'MAILLAGE',0,0,0,0,IPTVI1) 795 : NBELE6=NBELE6+1 796 : IPT6.NUM(1,NBELE6)=NRCI 797 : ENDIF 798 : C Determination de l'intersection entre le triangle NTI et 799 : C le polyedre NCI 800 : C Clipping selon l'algorithme de Sutherland et Hodgman 801 : C Le polygone resultant de l'intersection est decrit par la 802 : C table LTL1 avec : 803 : C LTL1.XACT(i,j) = coordonnee j du i-eme noeud du polygone 804 : C --> initialisation de LTL1 aux noeuds du triangle NTI 805 : NRTI1=IPT2.NUM(1,NTI) 806 : NRTI2=IPT2.NUM(2,NTI) 807 : NRTI3=IPT2.NUM(3,NTI) 808 : N1=3 809 : N2=3 810 : SEGINI,LTL1 811 : LTL1.XACT(1,1)=XCOOR((NRTI1-1)*IDIMP1+1) 812 : LTL1.XACT(1,2)=XCOOR((NRTI1-1)*IDIMP1+2) 813 : LTL1.XACT(1,3)=XCOOR((NRTI1-1)*IDIMP1+3) 814 : LTL1.XACT(2,1)=XCOOR((NRTI2-1)*IDIMP1+1) 815 : LTL1.XACT(2,2)=XCOOR((NRTI2-1)*IDIMP1+2) 816 : LTL1.XACT(2,3)=XCOOR((NRTI2-1)*IDIMP1+3) 817 : LTL1.XACT(3,1)=XCOOR((NRTI3-1)*IDIMP1+1) 818 : LTL1.XACT(3,2)=XCOOR((NRTI3-1)*IDIMP1+2) 819 : LTL1.XACT(3,3)=XCOOR((NRTI3-1)*IDIMP1+3) 820 : C --> table du polyedre NRCI (MTABI0) 820 : CALL ACCTAB(MTABCCC,'POINT ',0,0,' ',0,NRCI, 821 : & 'TABLE ',0,0,0,0,MTABI00) 822 : C --> liste des faces du polyedre NRCI 822 : CALL ACCTAB(MTABI00,'MOT ',0,0,'FACS',0,0, 823 : & 'LISTENTI',0,0,0,0,LFI00) 824 : C --> maillage IPT3 des voisins du polyedre NRCI 824 : CALL ACCTAB(MTABI00,'MOT ',0,0,'VOIS',0,0, 825 : & 'MAILLAGE',0,0,0,0,IPT3) 826 : C Boucle sur toutes les faces du polyedre NCI 827 : DO J=1,LFI00.LECT(/1) 828 : C Maillage IPT4 de la face J du polyedre NCI 829 : NFJ=LFI00.LECT(J) 829 : CALL ACCTAB(MTABFFF,'ENTIER ',NFJ,0,' ',0,0, 830 : & 'TABLE ',0,0,0,0,MTABF100) 830 : CALL ACCTAB(MTABF100,'MOT ',0,0,'VISU',0,0, 831 : & 'MAILLAGE',0,0,0,0,IPT4) 832 : NREFA=IPT4.NUM(1,1) 833 : NREFB=IPT4.NUM(1,2) 834 : NREFC=IPT4.NUM(1,3) 835 : C Massicotage du polygone LTL1 selon le demi espace definit 836 : C par le plan de IPT4 et le centre I du polyedre de Voronoi 837 : C Boucle sur les aretes du polygone LTL1 838 : C On utilise un polygone temporaire LTL2 839 : N1=LTL1.XACT(/1)+20 840 : N2=LTL1.XACT(/2) 841 : SEGINI,LTL2 842 : NLV2=0 843 : DO K=1,LTL1.XACT(/1) 844 : C Coordonnees des points P1 et P2 845 : X1=LTL1.XACT(K,1) 846 : Y1=LTL1.XACT(K,2) 847 : Z1=LTL1.XACT(K,3) 848 : IF (K.EQ.LTL1.XACT(/1)) THEN 849 : X2=LTL1.XACT(1,1) 850 : Y2=LTL1.XACT(1,2) 851 : Z2=LTL1.XACT(1,3) 852 : ELSE 853 : X2=LTL1.XACT(K+1,1) 854 : Y2=LTL1.XACT(K+1,2) 855 : Z2=LTL1.XACT(K+1,3) 856 : ENDIF 857 : C Intersection du segment P1P2 avec le demi plan ABC,I 857 : CALL INTPSEG(NREFA,NREFB,NREFC,NRCI,X1,Y1,Z1,X2,Y2,Z2, 858 : & IINT,XI1,YI1,ZI1,XI2,YI2,ZI2,ZERO) 859 : C Il y a 4 cas possibles 860 : C -1) le segment P1P2 est entierment hors du demi plan 861 : C on ajoute rien dans LTL2 862 : IF (IINT.EQ.0) THEN 863 : ENDIF 864 : C -2) le segment P1P2 rentre le demi plan 865 : C on ajoute le point d'intersection et P2 dans LTL2 866 : IF (IINT.EQ.1) THEN 867 : NLV2=NLV2+1 868 : LTL2.XACT(NLV2,1)=XI1 869 : LTL2.XACT(NLV2,2)=YI1 870 : LTL2.XACT(NLV2,3)=ZI1 871 : NLV2=NLV2+1 872 : LTL2.XACT(NLV2,1)=XI2 873 : LTL2.XACT(NLV2,2)=YI2 874 : LTL2.XACT(NLV2,3)=ZI2 875 : ENDIF 876 : C -3) le segment P1P2 est entierement dans le demi plan 877 : C on ajoute P2 dans LTL2 878 : IF (IINT.EQ.2) THEN 879 : NLV2=NLV2+1 880 : LTL2.XACT(NLV2,1)=XI2 881 : LTL2.XACT(NLV2,2)=YI2 882 : LTL2.XACT(NLV2,3)=ZI2 883 : ENDIF 884 : C -4) le segment P1P2 quitte le demi plan 885 : C on ajoute le point d'intersection dans LTL2 886 : IF (IINT.EQ.3) THEN 887 : NLV2=NLV2+1 888 : LTL2.XACT(NLV2,1)=XI2 889 : LTL2.XACT(NLV2,2)=YI2 890 : LTL2.XACT(NLV2,3)=ZI2 891 : ENDIF 892 : C Augmentation eventuelle de LTL2 893 : IF (NLV2.GE.(LTL2.XACT(/1)-1)) THEN 894 : N1=NLV2+20 895 : N2=LTL2.XACT(/2) 896 : SEGADJ,LTL2 897 : ENDIF 898 : ENDDO 899 : C Ajustement final de LTL2, ce polygone est l'intersection 900 : C de LTL1 avec le demi espace definit par la face J 901 : IF (NLV2.LT.3) THEN 902 : C WRITE(6,*) ' > Intersection triangle/cellule vide' 903 : C WRITE(6,*) 'NLV2=',nlv2 904 : C WRITE(6,*) 'Numero du triangle concerne:',nti 905 : C WRITE(6,*) 'Points de la face concernee',nrci,nrcj 906 : INTERR(1)=nti 907 : INTERR(2)=nrci 907 : CALL ERREUR(1033) 908 : GOTO 999 909 : ENDIF 910 : N1=NLV2 911 : N2=LTL2.XACT(/2) 912 : SEGADJ,LTL2 913 : C Le polygone LTL1 devient LTL2 puis on passe a la face 914 : C suivante 915 : LTL1=LTL2 916 : ENDDO 917 : C --> Elimination des "petites aretes" selon la tolerance XZERO 918 : C on fusionne les sommets du polygonne LTL1 trop proches 919 : N1=LTL1.XACT(/1) 920 : N2=LTL1.XACT(/2) 921 : SEGINI,LTL2 922 : NLTL2=0 923 : DO J=1,LTL1.XACT(/1) 924 : JJ=J+1 925 : IF (J.EQ.LTL1.XACT(/1)) JJ=1 926 : X1=LTL1.XACT(J,1) 927 : Y1=LTL1.XACT(J,2) 928 : Z1=LTL1.XACT(J,3) 929 : JJ=J+1 930 : IF (J.EQ.LTL1.XACT(/1)) JJ=1 931 : X2=LTL1.XACT(JJ,1) 932 : Y2=LTL1.XACT(JJ,2) 933 : Z2=LTL1.XACT(JJ,3) 934 : D12=SQRT(((X2-X1)**2)+((Y2-Y1)**2)+((Z2-Z1)**2)) 935 : c----------------------------------------------------------------------- 936 : c IF (D12.GE.XZERO) THEN 937 : IF (D12.GE.XTOL) THEN 938 : c----------------------------------------------------------------------- 939 : NLTL2=NLTL2+1 940 : LTL2.XACT(NLTL2,1)=X1 941 : LTL2.XACT(NLTL2,2)=Y1 942 : LTL2.XACT(NLTL2,3)=Z1 943 : ENDIF 944 : ENDDO 945 : N1=NLTL2 946 : N2=LTL2.XACT(/2) 947 : SEGADJ,LTL2 948 : LTL1=LTL2 949 : C 950 : C 951 : C Maintenant, on va incrementer la file (LT,LC) 952 : C avec les intersections voisines 953 : C Au fur et a mesure que l'on identifie les aretes du polygone 954 : C LTL1, on cree un maillage que l'on range dans MTAB11 955 : C 956 : C Nombre d'aretes identifiees (a comparer au nombre de points 957 : C du polygone LTL1 pour verification) 958 : NAI=0 959 : C 960 : C 1er cas : recherche des cellules voisines de NCI coupant le 961 : C triangle NTI 962 : NFF=NF 963 : C Boucle sur voisins du polyedre NCI (on ne s'interesse qu'aux 964 : C faces internes) 965 : DO J=1,IPT3.NUM(/2) 966 : NRCJ=IPT3.NUM(1,J) 967 : C Liste des faces du polyedre NRCJ 967 : CALL ACCTAB(MTABCCC,'POINT ',0,0,' ',0,NRCJ, 968 : & 'TABLE ',0,0,0,0,MTABJ00) 968 : CALL ACCTAB(MTABJ00,'MOT ',0,0,'FACS',0,0, 969 : & 'LISTENTI',0,0,0,0,LFJ00) 970 : C Recherche du numero de la face entre NRCI et NRCJ 971 : C cette face est a la fois dans LFI00 et LFJ00 972 : NFJ=0 973 : DO K=1,LFI00.LECT(/1) 974 : NFJ=LFI00.LECT(K) 975 : DO L=1,LFJ00.LECT(/1) 976 : IF (LFJ00.LECT(L).EQ.NFJ) GOTO 25 977 : ENDDO 978 : ENDDO 979 : 25 CONTINUE 980 : C Maillage IPT5 de la face J du polyedre NCI 980 : CALL ACCTAB(MTABFFF,'ENTIER ',NFJ,0,' ',0,0, 981 : & 'TABLE ',0,0,0,0,MTABF100) 981 : CALL ACCTAB(MTABF100,'MOT ',0,0,'VISU',0,0, 982 : & 'MAILLAGE',0,0,0,0,IPT5) 983 : C Points ABC definissant le plan de la face IPT5 984 : SEGINI,IPT7=IPT5 984 : CALL CHANGE(IPT7,1) 985 : NRA=IPT7.NUM(1,1) 986 : NRB=IPT7.NUM(1,2) 987 : XA=XCOOR((NRA-1)*IDIMP1+1) 988 : YA=XCOOR((NRA-1)*IDIMP1+2) 989 : ZA=XCOOR((NRA-1)*IDIMP1+3) 990 : XB=XCOOR((NRB-1)*IDIMP1+1) 991 : YB=XCOOR((NRB-1)*IDIMP1+2) 992 : ZB=XCOOR((NRB-1)*IDIMP1+3) 993 : C Boucle pour trouver un point C non aligne avec les points 994 : C A et B 995 : DO K=3,IPT7.NUM(/2) 996 : NRC=IPT7.NUM(1,K) 997 : XC=XCOOR((NRC-1)*IDIMP1+1) 998 : YC=XCOOR((NRC-1)*IDIMP1+2) 999 : ZC=XCOOR((NRC-1)*IDIMP1+3) 1000 : C Calcul d'un vecteur normal au plan (ABC) 1001 : XN=((YB-YA)*(ZC-ZA))-((ZB-ZA)*(YC-YA)) 1002 : YN=((ZB-ZA)*(XC-XA))-((XB-XA)*(ZC-ZA)) 1003 : ZN=((XB-XA)*(YC-YA))-((YB-YA)*(XC-XA)) 1004 : DN=SQRT((XN**2)+(YN**2)+(ZN**2)) 1005 : IF (DN.GE.ZERO2) GOTO 24 1006 : ENDDO 1007 : C Si l'on passe ici, c'est que l'on a pas trouve 3 points 1008 : C suffisament espaces pour definir le plan 1009 : C WRITE(6,*) '> Impossible de determiner l''equation du plan' 1010 : C WRITE(6,*) ' points alignes (1)' 1011 : INTERR(1)=NFJ 1012 : INTERR(2)=NRCI 1012 : CALL ERREUR(1035) 1013 : GOTO 999 1014 : 24 CONTINUE 1015 : XN=XN/DN 1016 : YN=YN/DN 1017 : ZN=ZN/DN 1018 : W=-1D0*((XN*XA)+(YN*YA)+(ZN*ZA)) 1019 : C Recherche des aretes du polygone LTL1 dans le plan ABC 1020 : DO K=1,LTL1.XACT(/1) 1021 : IAE=0 1022 : C Coordonnees des points 41 et 42 de l'arete K 1023 : X41=LTL1.XACT(K,1) 1024 : Y41=LTL1.XACT(K,2) 1025 : Z41=LTL1.XACT(K,3) 1026 : IF (K.EQ.LTL1.XACT(/1)) THEN 1027 : X42=LTL1.XACT(1,1) 1028 : Y42=LTL1.XACT(1,2) 1029 : Z42=LTL1.XACT(1,3) 1030 : ELSE 1031 : X42=LTL1.XACT(K+1,1) 1032 : Y42=LTL1.XACT(K+1,2) 1033 : Z42=LTL1.XACT(K+1,3) 1034 : ENDIF 1035 : C Test si P41 et P42 sont dans le plan ABC 1036 : D41=(XN*X41)+(YN*Y41)+(ZN*Z41)+W 1037 : D42=(XN*X42)+(YN*Y42)+(ZN*Z42)+W 1038 : IF ((ABS(D41).LT.ZERO2).AND.(ABS(D42).LT.ZERO2)) THEN 1039 : C L'arete K du polygone LTL1 est dans le plan ABC 1040 : C Identification de la cellule associee au point NRCJ 1041 : NCJ=0 1042 : DO L=1,IPT1.NUM(/2) 1043 : NRTEST=IPT1.NUM(1,L) 1044 : IF (NRTEST.EQ.NRCJ) THEN 1045 : NCJ=L 1046 : GOTO 2 1047 : ENDIF 1048 : ENDDO 1049 : 2 CONTINUE 1050 : C On verifie que le couple (NTI,NCJ) n'a pas deja 1051 : C ete traite et on l'ajoute a la file 1052 : IF (ITL1.LACT(NTI,NCJ).EQ.0) THEN 1053 : ITL1.LACT(NTI,NCJ)=1 1054 : NLT=NLT+1 1055 : LT.LECT(NLT)=NTI 1056 : LC.LECT(NLT)=NCJ 1057 : ELSE 1058 : ENDIF 1059 : C Determination des noeuds P1 et P2 de l'arete 1060 : C Doit on creer de nouveaux noeuds ? 1061 : C LV2 liste des noeuds a examiner = LV1 et les 1062 : C noeuds de la face I,J 1063 : JG=LV1.LECT(/1) 1064 : SEGINI,LV2=LV1 1065 : JG=LV2.LECT(/1)+IPT5.NUM(/2) 1066 : SEGADJ,LV2 1067 : NLV2=LV1.LECT(/1) 1068 : DO L=1,IPT5.NUM(/2) 1069 : NLV2=NLV2+1 1070 : LV2.LECT(NLV2)=IPT5.NUM(1,L) 1071 : ENDDO 1072 : C test si P1 est deja present ou non dans LV2 1073 : BOOL2=.FALSE. 1074 : DO L=1,LV2.LECT(/1) 1075 : NRL=LV2.LECT(L) 1076 : XL=XCOOR((NRL-1)*IDIMP1+1) 1077 : YL=XCOOR((NRL-1)*IDIMP1+2) 1078 : ZL=XCOOR((NRL-1)*IDIMP1+3) 1079 : IF ((ABS(XL-X41).LT.XTOL).AND. 1080 : & (ABS(YL-Y41).LT.XTOL).AND. 1081 : & (ABS(ZL-Z41).LT.XTOL)) THEN 1082 : BOOL2=.TRUE. 1083 : NRP1=NRL 1084 : GOTO 4 1085 : ENDIF 1086 : ENDDO 1087 : IF (.NOT.BOOL2) THEN 1088 : C ajout d'un nouveau noeud dans MCOORD 1089 : NRP1=(XCOOR(/1)/IDIMP1)+1 1090 : NBPTS=NRP1 1091 : SEGADJ,MCOORD 1092 : XCOOR((NRP1-1)*IDIMP1+1)=X41 1093 : XCOOR((NRP1-1)*IDIMP1+2)=Y41 1094 : XCOOR((NRP1-1)*IDIMP1+3)=Z41 1095 : NLV1=NLV1+1 1096 : IF (NLV1.GT.LV1.LECT(/1)) THEN 1097 : JG=NLV1 1098 : SEGADJ,LV1 1099 : ENDIF 1100 : LV1.LECT(NLV1)=NRP1 1101 : ENDIF 1102 : 4 CONTINUE 1103 : C test si P2 est deja present ou non dans LV2 1104 : BOOL2=.FALSE. 1105 : DO L=1,LV2.LECT(/1) 1106 : NRL=LV2.LECT(L) 1107 : XL=XCOOR((NRL-1)*IDIMP1+1) 1108 : YL=XCOOR((NRL-1)*IDIMP1+2) 1109 : ZL=XCOOR((NRL-1)*IDIMP1+3) 1110 : IF ((ABS(XL-X42).LT.XTOL).AND. 1111 : & (ABS(YL-Y42).LT.XTOL).AND. 1112 : & (ABS(ZL-Z42).LT.XTOL)) THEN 1113 : BOOL2=.TRUE. 1114 : NRP2=NRL 1115 : GOTO 5 1116 : ENDIF 1117 : ENDDO 1118 : IF (.NOT.BOOL2) THEN 1119 : C ajout d'un nouveau noeud dans MCOORD 1120 : NRP2=(XCOOR(/1)/IDIMP1)+1 1121 : NBPTS=NRP2 1122 : SEGADJ,MCOORD 1123 : XCOOR((NRP2-1)*IDIMP1+1)=X42 1124 : XCOOR((NRP2-1)*IDIMP1+2)=Y42 1125 : XCOOR((NRP2-1)*IDIMP1+3)=Z42 1126 : NLV1=NLV1+1 1127 : IF (NLV1.GT.LV1.LECT(/1)) THEN 1128 : JG=NLV1 1129 : SEGADJ,LV1 1130 : ENDIF 1131 : LV1.LECT(NLV1)=NRP2 1132 : ENDIF 1133 : 5 CONTINUE 1134 : C Test sur l'existence de l'arete P1 P2 dans IPTG 1135 : DO L=1,NBELEG 1136 : NG1=IPTG.NUM(1,L) 1137 : NG2=IPTG.NUM(2,L) 1138 : IF (((NG1.EQ.NRP1).AND.(NG2.EQ.NRP2)).OR. 1139 : & ((NG2.EQ.NRP1).AND.(NG1.EQ.NRP2))) THEN 1140 : C si l'arete P1 P2 a deja ete creee, on 1141 : C l'indique au moyen de IAE 1142 : IAE=1 1143 : NAE=L 1144 : GOTO 21 1145 : ENDIF 1146 : ENDDO 1147 : 21 CONTINUE 1148 : C --> Ajout de l'arete dans MTAB11.VISU 1149 : IF (IAE.EQ.0) THEN 1150 : NBELEG=NBELEG+1 1151 : IPTG.NUM(1,NBELEG)=NRP1 1152 : IPTG.NUM(2,NBELEG)=NRP2 1153 : IF (NBELEG.EQ.IPTG.NUM(/2)) THEN 1154 : NBNN=IPTG.NUM(/1) 1155 : NBELEM=IPTG.NUM(/2)+100 1156 : SEGADJ,IPTG 1157 : ENDIF 1158 : NAE=NBELEG 1159 : ENDIF 1160 : C --> Ajout de l'arete dans MTAB11.ARTS 1161 : IF (IAE.EQ.0) THEN 1162 : NBNN=2 1163 : NBELEM=1 1164 : SEGINI,IPTA1 1165 : IPTA1.ITYPEL=2 1166 : IPTA1.NUM(1,1)=NRP1 1167 : IPTA1.NUM(2,1)=NRP2 1167 : CALL ECCTAB(MTABA,'ENTIER ',NAE,0,' ',0,0, 1168 : & 'MAILLAGE',0,0,0,0,IPTA1) 1169 : ENDIF 1170 : C --> Ajout de l'arete dans MTAB11.FACS 1170 : CALL ECROBJ('ENTIER ',NF) 1170 : CALL ECROBJ('TABLE ',MTABF) 1170 : CALL EXIS 1170 : CALL LIRLOG(BOOL1,1,IRETOU) 1171 : IF (.NOT.BOOL1) THEN 1172 : C si MTAB11.FACS.NF n'existe pas, on le cree 1172 : CALL CRTABL(MTABF1) 1172 : CALL ECCTAB(MTABF,'ENTIER ',NF,0,' ',0,0, 1173 : & 'TABLE ',0,0,0,0,MTABF1) 1174 : NBNN=2 1175 : NBELEM=0 1176 : SEGINI,IPTF1 1177 : IPTF1.ITYPEL=2 1177 : CALL ECCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 1178 : & 'MAILLAGE',0,0,0,0,IPTF1) 1179 : JG=0 1180 : SEGINI,LA1 1180 : CALL ECCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, 1181 : & 'LISTENTI',0,0,0,0,LA1) 1182 : ENDIF 1183 : C acces a la table de la face NF 1183 : CALL ACCTAB(MTABF,'ENTIER ',NF,0,' ',0,0, 1184 : & 'TABLE ',0,0,0,0,MTABF1) 1185 : C ajout d'une arete dans l'indice VISU 1185 : CALL ACCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 1186 : & 'MAILLAGE',0,0,0,0,IPTF1) 1187 : NBNN=IPTF1.NUM(/1) 1188 : NBELEM=IPTF1.NUM(/2)+1 1189 : SEGADJ,IPTF1 1190 : IPTF1.NUM(1,NBELEM)=NRP1 1191 : IPTF1.NUM(2,NBELEM)=NRP2 1192 : C ajout du numero de l'arete l'indice ARTS 1192 : CALL ACCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, 1193 : & 'LISTENTI',0,0,0,0,LA1) 1194 : JG=LA1.LECT(/1)+1 1195 : SEGADJ,LA1 1196 : LA1.LECT(JG)=NAE 1197 : C --> Ecriture dans MTAB11.CELL.NRCI 1198 : C ajout d'une arete dans l'indice VISU si besoin 1199 : DO L=1,IPTI1.NUM(/2) 1200 : NI1=IPTI1.NUM(1,L) 1201 : NI2=IPTI1.NUM(2,L) 1202 : IF (((NI1.EQ.NRP1).AND.(NI2.EQ.NRP2)).OR. 1203 : & ((NI2.EQ.NRP1).AND.(NI1.EQ.NRP2))) GOTO 15 1204 : ENDDO 1205 : NBNN=IPTI1.NUM(/1) 1206 : NBELEM=IPTI1.NUM(/2)+1 1207 : SEGADJ,IPTI1 1208 : IPTI1.NUM(1,NBELEM)=NRP1 1209 : IPTI1.NUM(2,NBELEM)=NRP2 1210 : 15 CONTINUE 1211 : C ajout du numero de la face dans l'indice FACS 1212 : C si besoin 1212 : CALL ECRENT(NF) 1212 : CALL ECROBJ('LISTENTI',LFI1) 1212 : CALL EXIS 1212 : CALL LIRLOG(BOOL1,1,IRETOU) 1213 : SEGACT,LFI1 1214 : IF (.NOT.BOOL1) THEN 1215 : JG=LFI1.LECT(/1)+1 1216 : SEGADJ,LFI1 1217 : LFI1.LECT(JG)=NF 1218 : ENDIF 1219 : C ajout de NRCJ dans l'indice VOIS si besion 1219 : CALL ECROBJ('MAILLAGE',IPTVI1) 1219 : CALL ECROBJ('POINT ',NRCJ) 1219 : CALL DANS 1220 : SEGACT,IPTVI1 1220 : CALL LIRLOG(BOOL1,1,IRETOU) 1221 : IF (.NOT.BOOL1) THEN 1222 : NBNN=IPTVI1.NUM(/1) 1223 : NBELEM=IPTVI1.NUM(/2)+1 1224 : SEGADJ,IPTVI1 1225 : IPTVI1.NUM(1,NBELEM)=NRCJ 1226 : ENDIF 1227 : C Cette arete appartient aussi a la face interne 1228 : C situee entre NRCI et NRCJ 1229 : C --> Ecriture dans MTAB11.CELL.NRCJ 1230 : MMTABI1=MTABI1 1231 : IIPTI1=IPTI1 1232 : LLFI1=LFI1 1233 : IIPTVI1=IPTVI1 1233 : CALL ECROBJ('POINT ',NRCJ) 1233 : CALL ECROBJ('TABLE ',MTABC) 1233 : CALL EXIS 1233 : CALL LIRLOG(BOOL1,1,IRETOU) 1234 : C si MATB11.CELL.NRCJ n'existe pas on le cree 1235 : IF (.NOT.BOOL1) THEN 1235 : CALL CRTABL(MTABI1) 1235 : CALL ECCTAB(MTABC,'POINT ',0,0,' ',0,NRCJ, 1236 : & 'TABLE ',0,0,0,0,MTABI1) 1237 : NBNN=2 1238 : NBELEM=0 1239 : SEGINI,IPTI1 1240 : IPTI1.ITYPEL=2 1240 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 1241 : & 'MAILLAGE',0,0,0,0,IPTI1) 1242 : JG=0 1243 : SEGINI,LFI1 1243 : CALL ECCTAB(MTABI1,'MOT ',0,0,'FACS',0,0, 1244 : & 'LISTENTI',0,0,0,0,LFI1) 1245 : NBNN=1 1246 : NBELEM=0 1247 : SEGINI,IPTVI1 1248 : IPTVI1.ITYPEL=1 1248 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, 1249 : & 'MAILLAGE',0,0,0,0,IPTVI1) 1250 : NBELE6=NBELE6+1 1251 : IPT6.NUM(1,NBELE6)=NRCJ 1252 : ENDIF 1253 : C acces a la table de la cellule NRCJ 1253 : CALL ACCTAB(MTABC,'POINT ',0,0,' ',0,NRCJ, 1254 : & 'TABLE ',0,0,0,0,MTABI1) 1254 : CALL ACCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 1255 : & 'MAILLAGE',0,0,0,0,IPTI1) 1255 : CALL ACCTAB(MTABI1,'MOT ',0,0,'FACS',0,0, 1256 : & 'LISTENTI',0,0,0,0,LFI1) 1256 : CALL ACCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, 1257 : & 'MAILLAGE',0,0,0,0,IPTVI1) 1258 : C ajout d'une arete dans l'indice VISU si besoin 1259 : DO L=1,IPTI1.NUM(/2) 1260 : NI1=IPTI1.NUM(1,L) 1261 : NI2=IPTI1.NUM(2,L) 1262 : IF (((NI1.EQ.NRP1).AND.(NI2.EQ.NRP2)).OR. 1263 : & ((NI2.EQ.NRP1).AND.(NI1.EQ.NRP2))) GOTO 16 1264 : ENDDO 1265 : NBNN=IPTI1.NUM(/1) 1266 : NBELEM=IPTI1.NUM(/2)+1 1267 : SEGADJ,IPTI1 1268 : IPTI1.NUM(1,NBELEM)=NRP1 1269 : IPTI1.NUM(2,NBELEM)=NRP2 1270 : 16 CONTINUE 1271 : C s'agit il d'une nouvelle face ou pas et si oui 1272 : C laquelle ? 1273 : NFI=ITFAC1.LACT(NCI,NCJ) 1274 : IF (NFI.EQ.0) THEN 1275 : C on renseigne cette nouvelle face dans ITFAC1 1276 : NFF=NFF+1 1277 : ITFAC1.LACT(NCI,NCJ)=NFF 1278 : ITFAC1.LACT(NCJ,NCI)=NFF 1279 : NFI=NFF 1280 : C ainsi que dans LV5 et LV6 1281 : IF (NFF.GT.LV5.LECT(/1)) THEN 1282 : JG=NFF+100 1283 : SEGADJ,LV5,LV6 1284 : ENDIF 1285 : IF (LV5.LECT(NFF).EQ.0) THEN 1286 : LV5.LECT(NFF)=NCI 1287 : LV6.LECT(NFF)=NCJ 1288 : ENDIF 1289 : ENDIF 1290 : C ajout du numero de la face dans l'indice FACS 1291 : C si besoin 1291 : CALL ECRENT(NFI) 1291 : CALL ECROBJ('LISTENTI',LFI1) 1291 : CALL EXIS 1291 : CALL LIRLOG(BOOL1,1,IRETOU) 1292 : SEGACT,LFI1 1293 : IF (.NOT.BOOL1) THEN 1294 : JG=LFI1.LECT(/1)+1 1295 : SEGADJ,LFI1 1296 : LFI1.LECT(JG)=NFI 1297 : ENDIF 1298 : C ajout de NRCI dans l'indice VOIS si besion 1299 : BOOL1=.FALSE. 1300 : DO L=1,IPTVI1.NUM(/2) 1301 : IF (IPTVI1.NUM(1,L).EQ.NRCI) BOOL1=.TRUE. 1302 : ENDDO 1303 : IF (.NOT.BOOL1) THEN 1304 : NBNN=IPTVI1.NUM(/1) 1305 : NBELEM=IPTVI1.NUM(/2)+1 1306 : SEGADJ,IPTVI1 1307 : IPTVI1.NUM(1,NBELEM)=NRCI 1308 : ENDIF 1309 : MTABI1=MMTABI1 1310 : IPTI1=IIPTI1 1311 : LFI1=LLFI1 1312 : IPTVI1=IIPTVI1 1313 : C ajout de l'arete dans MTAB11.FACS.NFI 1314 : MMTABF1=MTABF1 1315 : IIPTF1=IPTF1 1316 : LLA1=LA1 1316 : CALL ECROBJ('ENTIER ',NFI) 1316 : CALL ECROBJ('TABLE ',MTABF) 1316 : CALL EXIS 1316 : CALL LIRLOG(BOOL1,1,IRETOU) 1317 : IF (.NOT.BOOL1) THEN 1318 : C si MTAB11.FACS.NFI n'existe pas, on le cree 1318 : CALL CRTABL(MTABF1) 1318 : CALL ECCTAB(MTABF,'ENTIER ',NFI,0,' ',0,0, 1319 : & 'TABLE ',0,0,0,0,MTABF1) 1320 : NBNN=2 1321 : NBELEM=0 1322 : SEGINI,IPTF1 1323 : IPTF1.ITYPEL=2 1323 : CALL ECCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 1324 : & 'MAILLAGE',0,0,0,0,IPTF1) 1325 : JG=0 1326 : SEGINI,LA1 1326 : CALL ECCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, 1327 : & 'LISTENTI',0,0,0,0,LA1) 1328 : ENDIF 1329 : C acces a la table de la face NFI 1329 : CALL ACCTAB(MTABF,'ENTIER ',NFI,0,' ',0,0, 1330 : & 'TABLE ',0,0,0,0,MTABF1) 1331 : C ajout d'une arete dans la face si besoin 1331 : CALL ACCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 1332 : & 'MAILLAGE',0,0,0,0,IPTF1) 1332 : CALL ACCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, 1333 : & 'LISTENTI',0,0,0,0,LA1) 1333 : CALL ECRENT(NAE) 1333 : CALL ECROBJ('LISTENTI',LA1) 1333 : CALL EXIS 1334 : SEGACT,LA1 1334 : CALL LIRLOG(BOOL1,1,IRETOU) 1335 : IF (.NOT.BOOL1) THEN 1336 : JG=LA1.LECT(/1)+1 1337 : SEGADJ,LA1 1338 : LA1.LECT(JG)=NAE 1339 : NBNN=IPTF1.NUM(/1) 1340 : NBELEM=IPTF1.NUM(/2)+1 1341 : SEGADJ,IPTF1 1342 : IPTF1.NUM(1,NBELEM)=NRP1 1343 : IPTF1.NUM(2,NBELEM)=NRP2 1344 : ENDIF 1345 : MTABF1=MMTABF1 1346 : IPTF1=IIPTF1 1347 : LA1=LLA1 1348 : C On incremente le nombre d'aretes identifiees 1349 : NAI=NAI+1 1350 : ENDIF 1351 : ENDDO 1352 : ENDDO 1353 : C 1354 : C 2eme cas : recherche des triangles voisins de NTI coupant 1355 : C la cellule NCI 1356 : C Boucle sur les cotes du triangle NTI 1357 : DO J=1,3 1358 : IF (J.EQ.1) THEN 1359 : NRA=IPT2.NUM(2,NTI) 1360 : NRB=IPT2.NUM(3,NTI) 1361 : ENDIF 1362 : IF (J.EQ.2) THEN 1363 : NRA=IPT2.NUM(1,NTI) 1364 : NRB=IPT2.NUM(3,NTI) 1365 : ENDIF 1366 : IF (J.EQ.3) THEN 1367 : NRA=IPT2.NUM(1,NTI) 1368 : NRB=IPT2.NUM(2,NTI) 1369 : ENDIF 1370 : XA=XCOOR((NRA-1)*IDIMP1+1) 1371 : YA=XCOOR((NRA-1)*IDIMP1+2) 1372 : ZA=XCOOR((NRA-1)*IDIMP1+3) 1373 : XB=XCOOR((NRB-1)*IDIMP1+1) 1374 : YB=XCOOR((NRB-1)*IDIMP1+2) 1375 : ZB=XCOOR((NRB-1)*IDIMP1+3) 1376 : C La droite (AB) est definie par les deux equations de plan 1377 : C Ax + By + Cy + D = 0 1378 : C Ex + Fy + Gy + H = 0 1379 : C --> pour le plan 1, on definit un point C dans ce plan 1380 : XC=XA+(YB-YA) 1381 : YC=YA+(ZB-ZA) 1382 : ZC=ZA+(XB-XA) 1383 : C --> vecteur normal au plan ABC : n=AB^AC 1384 : A=((YB-YA)*(ZC-ZA))-((ZB-ZA)*(YC-YA)) 1385 : B=((ZB-ZA)*(XC-XA))-((XB-XA)*(ZC-ZA)) 1386 : C=((XB-XA)*(YC-YA))-((YB-YA)*(XC-XA)) 1387 : C --> calcul de D, sachant que A est dans le plan 1388 : D=-1D0*((A*XA)+(B*YA)+(C*ZA)) 1389 : C --> pour le plan 2, on definit un point D dans ce plan 1390 : XD=XA+(ZB-ZA) 1391 : YD=YA+(XB-XA) 1392 : ZD=ZA+(YB-YA) 1393 : C --> vecteur normal au plan ABD : n=AB^AD 1394 : E=((YB-YA)*(ZD-ZA))-((ZB-ZA)*(YD-YA)) 1395 : F=((ZB-ZA)*(XD-XA))-((XB-XA)*(ZD-ZA)) 1396 : G=((XB-XA)*(YD-YA))-((YB-YA)*(XD-XA)) 1397 : C --> calcul de H, sachant que A est dans le plan 1398 : H=-1D0*((E*XA)+(F*YA)+(G*ZA)) 1399 : C Recherche des aretes du polygone LTL1 sur la droite AB 1400 : DO K=1,LTL1.XACT(/1) 1401 : IAE=0 1402 : X41=LTL1.XACT(K,1) 1403 : Y41=LTL1.XACT(K,2) 1404 : Z41=LTL1.XACT(K,3) 1405 : IF (K.EQ.LTL1.XACT(/1)) THEN 1406 : X42=LTL1.XACT(1,1) 1407 : Y42=LTL1.XACT(1,2) 1408 : Z42=LTL1.XACT(1,3) 1409 : ELSE 1410 : X42=LTL1.XACT(K+1,1) 1411 : Y42=LTL1.XACT(K+1,2) 1412 : Z42=LTL1.XACT(K+1,3) 1413 : ENDIF 1414 : C Test si les points 1 et 2 de l'arete K sont sur la 1415 : C droite (AB) 1416 : I41=0 1417 : I42=0 1418 : C Test si P1 est sur la droite (AB) 1419 : V41=(A*X41)+(B*Y41)+(C*Z41)+D 1420 : W41=(E*X41)+(F*Y41)+(G*Z41)+H 1421 : IF ((ABS(V41).LT.ZERO2).AND.(ABS(W41).LT.ZERO2)) THEN 1422 : I41=1 1423 : ELSE 1424 : GOTO 3 1425 : ENDIF 1426 : C Test si P2 est sur la droite (AB) 1427 : V42=(A*X42)+(B*Y42)+(C*Z42)+D 1428 : W42=(E*X42)+(F*Y42)+(G*Z42)+H 1429 : IF ((ABS(V42).LT.ZERO2).AND.(ABS(W42).LT.ZERO2)) THEN 1430 : I42=1 1431 : ELSE 1432 : GOTO 3 1433 : ENDIF 1434 : C Si N41 et N42 sont sur l'arete [AB], on va chercher le 1435 : C numero du triangle adjacent sur cette arete 1436 : IF ((I41.EQ.1).AND.(I42.EQ.1)) THEN 1437 : NTVI=ILEA1.LEAC(NTI,J,1) 1438 : C On verifie que le couple (NTVI,NCI) n'a pas deja 1439 : C ete traite et on l'ajoute a la file 1440 : IF (ITL1.LACT(NTVI,NCI).EQ.0) THEN 1441 : ITL1.LACT(NTVI,NCI)=1 1442 : NLT=NLT+1 1443 : LT.LECT(NLT)=NTVI 1444 : LC.LECT(NLT)=NCI 1445 : ENDIF 1446 : ENDIF 1447 : C Determination des noeuds P1 et P2 de l'arete 1448 : C Doit on creer de nouveaux noeuds ? 1449 : C LV2 liste des noeuds a examiner = LV1 et NRA et NRB 1450 : JG=LV1.LECT(/1) 1451 : SEGINI,LV2=LV1 1452 : JG=LV2.LECT(/1)+2 1453 : SEGADJ,LV2 1454 : NLV2=LV1.LECT(/1) 1455 : NLV2=NLV2+1 1456 : LV2.LECT(NLV2)=NRA 1457 : NLV2=NLV2+1 1458 : LV2.LECT(NLV2)=NRB 1459 : C test si P1 est deja present ou non dans LV2 1460 : BOOL1=.FALSE. 1461 : DO L=1,LV2.LECT(/1) 1462 : NRL=LV2.LECT(L) 1463 : XL=XCOOR((NRL-1)*IDIMP1+1) 1464 : YL=XCOOR((NRL-1)*IDIMP1+2) 1465 : ZL=XCOOR((NRL-1)*IDIMP1+3) 1466 : IF ((ABS(XL-X41).LT.XTOL).AND. 1467 : & (ABS(YL-Y41).LT.XTOL).AND. 1468 : & (ABS(ZL-Z41).LT.XTOL)) THEN 1469 : BOOL1=.TRUE. 1470 : NRP1=NRL 1471 : GOTO 6 1472 : ENDIF 1473 : ENDDO 1474 : IF (.NOT.BOOL1) THEN 1475 : C ajout d'un nouveau noeud dans MCOORD 1476 : NRP1=(XCOOR(/1)/IDIMP1)+1 1477 : NBPTS=NRP1 1478 : SEGADJ,MCOORD 1479 : XCOOR((NRP1-1)*IDIMP1+1)=X41 1480 : XCOOR((NRP1-1)*IDIMP1+2)=Y41 1481 : XCOOR((NRP1-1)*IDIMP1+3)=Z41 1482 : NLV1=NLV1+1 1483 : IF (NLV1.GT.LV1.LECT(/1)) THEN 1484 : JG=NLV1 1485 : SEGADJ,LV1 1486 : ENDIF 1487 : LV1.LECT(NLV1)=NRP1 1488 : ENDIF 1489 : 6 CONTINUE 1490 : C test si P2 est deja present ou non dans LV2 1491 : BOOL1=.FALSE. 1492 : DO L=1,LV2.LECT(/1) 1493 : NRL=LV2.LECT(L) 1494 : XL=XCOOR((NRL-1)*IDIMP1+1) 1495 : YL=XCOOR((NRL-1)*IDIMP1+2) 1496 : ZL=XCOOR((NRL-1)*IDIMP1+3) 1497 : IF ((ABS(XL-X42).LT.XTOL).AND. 1498 : & (ABS(YL-Y42).LT.XTOL).AND. 1499 : & (ABS(ZL-Z42).LT.XTOL)) THEN 1500 : BOOL1=.TRUE. 1501 : NRP2=NRL 1502 : GOTO 7 1503 : ENDIF 1504 : ENDDO 1505 : IF (.NOT.BOOL1) THEN 1506 : NRP2=(XCOOR(/1)/IDIMP1)+1 1507 : NBPTS=NRP2 1508 : SEGADJ,MCOORD 1509 : XCOOR((NRP2-1)*IDIMP1+1)=X42 1510 : XCOOR((NRP2-1)*IDIMP1+2)=Y42 1511 : XCOOR((NRP2-1)*IDIMP1+3)=Z42 1512 : NLV1=NLV1+1 1513 : IF (NLV1.GT.LV1.LECT(/1)) THEN 1514 : JG=NLV1 1515 : SEGADJ,LV1 1516 : ENDIF 1517 : LV1.LECT(NLV1)=NRP2 1518 : ENDIF 1519 : 7 CONTINUE 1520 : C Test sur l'existence de l'arete P1 P2 dans IPTG 1521 : DO L=1,NBELEG 1522 : NG1=IPTG.NUM(1,L) 1523 : NG2=IPTG.NUM(2,L) 1524 : IF (((NG1.EQ.NRP1).AND.(NG2.EQ.NRP2)).OR. 1525 : & ((NG2.EQ.NRP1).AND.(NG1.EQ.NRP2))) THEN 1526 : C si l'arete P1 P2 a deja ete creee, on 1527 : C l'indique au moyen de IAE 1528 : IAE=1 1529 : NAE=L 1530 : GOTO 19 1531 : ENDIF 1532 : ENDDO 1533 : 19 CONTINUE 1534 : C --> Ajout de l'arete dans MTAB11.VISU 1535 : IF (IAE.EQ.0) THEN 1536 : NBELEG=NBELEG+1 1537 : IPTG.NUM(1,NBELEG)=NRP1 1538 : IPTG.NUM(2,NBELEG)=NRP2 1539 : IF (NBELEG.EQ.IPTG.NUM(/2)) THEN 1540 : NBNN=IPTG.NUM(/1) 1541 : NBELEM=IPTG.NUM(/2)+100 1542 : SEGADJ,IPTG 1543 : ENDIF 1544 : NAE=NBELEG 1545 : ENDIF 1546 : C --> Ajout de l'arete dans MTAB11.ARTS 1547 : IF (IAE.EQ.0) THEN 1548 : NBNN=2 1549 : NBELEM=1 1550 : SEGINI,IPTA1 1551 : IPTA1.ITYPEL=2 1552 : IPTA1.NUM(1,1)=NRP1 1553 : IPTA1.NUM(2,1)=NRP2 1553 : CALL ECCTAB(MTABA,'ENTIER ',NAE,0,' ',0,0, 1554 : & 'MAILLAGE',0,0,0,0,IPTA1) 1555 : ENDIF 1556 : C --> Ajout de l'arete dans MTAB11.FACS 1556 : CALL ECROBJ('ENTIER ',NF) 1556 : CALL ECROBJ('TABLE ',MTABF) 1556 : CALL EXIS 1556 : CALL LIRLOG(BOOL1,1,IRETOU) 1557 : IF (.NOT.BOOL1) THEN 1558 : C si MTAB11.FACS.NF n'existe pas, on le cree 1558 : CALL CRTABL(MTABF1) 1558 : CALL ECCTAB(MTABF,'ENTIER ',NF,0,' ',0,0, 1559 : & 'TABLE ',0,0,0,0,MTABF1) 1560 : NBNN=2 1561 : NBELEM=0 1562 : SEGINI,IPTF1 1563 : IPTF1.ITYPEL=2 1563 : CALL ECCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 1564 : & 'MAILLAGE',0,0,0,0,IPTF1) 1565 : JG=0 1566 : SEGINI,LA1 1566 : CALL ECCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, 1567 : & 'LISTENTI',0,0,0,0,LA1) 1568 : ENDIF 1569 : C acces a MTAB11.FACS.NF 1569 : CALL ACCTAB(MTABF,'ENTIER ',NF,0,' ',0,0, 1570 : & 'TABLE ',0,0,0,0,MTABF1) 1571 : C ajout d'une arete dans l'indice VISU 1571 : CALL ACCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 1572 : & 'MAILLAGE',0,0,0,0,IPTF1) 1573 : NBNN=IPTF1.NUM(/1) 1574 : NBELEM=IPTF1.NUM(/2)+1 1575 : SEGADJ,IPTF1 1576 : IPTF1.NUM(1,NBELEM)=NRP1 1577 : IPTF1.NUM(2,NBELEM)=NRP2 1578 : C ajout du numero de l'arete l'indice ARTS 1578 : CALL ACCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, 1579 : & 'LISTENTI',0,0,0,0,LA1) 1580 : JG=LA1.LECT(/1)+1 1581 : SEGADJ,LA1 1582 : LA1.LECT(JG)=NAE 1583 : C --> Ecriture dans MTAB11.CELL.NRCI 1584 : C ajout d'une arete dans l'indice VISU si besoin 1585 : DO L=1,IPTI1.NUM(/2) 1586 : NI1=IPTI1.NUM(1,L) 1587 : NI2=IPTI1.NUM(2,L) 1588 : IF (((NI1.EQ.NRP1).AND.(NI2.EQ.NRP2)).OR. 1589 : & ((NI2.EQ.NRP1).AND.(NI1.EQ.NRP2))) GOTO 17 1590 : ENDDO 1591 : NBNN=IPTI1.NUM(/1) 1592 : NBELEM=IPTI1.NUM(/2)+1 1593 : SEGADJ,IPTI1 1594 : IPTI1.NUM(1,NBELEM)=NRP1 1595 : IPTI1.NUM(2,NBELEM)=NRP2 1596 : 17 CONTINUE 1597 : C ajout du numero de la face dans l'indice FACS 1598 : C si besoin 1598 : CALL ECRENT(NF) 1598 : CALL ECROBJ('LISTENTI',LFI1) 1598 : CALL EXIS 1598 : CALL LIRLOG(BOOL1,1,IRETOU) 1599 : SEGACT,LFI1 1600 : IF (.NOT.BOOL1) THEN 1601 : JG=LFI1.LECT(/1)+1 1602 : SEGADJ,LFI1 1603 : LFI1.LECT(JG)=NF 1604 : ENDIF 1605 : NAI=NAI+1 1606 : 3 CONTINUE 1607 : ENDDO 1608 : ENDDO 1609 : NF=NFF 1610 : C 1611 : C Test pour verifier que le nombre d'aretes identifiees est 1612 : C bien egal au nombre de sommets du polygone d'intersection 1613 : IF (NAI.NE.LTL1.XACT(/1)) THEN 1614 : C WRITE(6,*) '> Le nombre d''aretes identifiees est' 1615 : C WRITE(6,*) ' different du nombre de sommets du polygone' 1616 : C WRITE(6,*) nai,LTL1.XACT(/1) 1617 : C WRITE(6,*) 'Numero du triangle concerne:',nti 1618 : C WRITE(6,*) 'Point de la cellule concernee',nrci 1618 : CALL ERREUR(1036) 1619 : GOTO 999 1620 : ENDIF 1621 : ENDDO 1622 : 1 CONTINUE 1623 : C Ajustement du maillage global des centres des cellules dans la 1624 : C partition coupee IPT6 1625 : NBNN=IPT6.NUM(/1) 1626 : NBELEM=NBELE6 1627 : SEGADJ,IPT6 1628 : C Ajustement du maillage global des aretes IPTG 1629 : NBNN=IPTG.NUM(/1) 1630 : NBELEM=NBELEG 1631 : SEGADJ,IPTG 1632 : C Ajustement de la liste des faces internes LV5 et LV6 1633 : JG=NF 1634 : SEGADJ,LV5,LV6 1635 : C 1636 : C 1637 : C Fusion des faces coplanaires 1638 : C On construit la liste des faces regroupees : LV1 1639 : C si LV1(n1)=n1 la face n1 est seule (interne ou englobante) 1640 : C si LV1(n1)=n2 la face n1 est reliee a la face n2 qui constitue 1641 : C sa face "englobante" et donc LV1(n2)=n2 1642 : JG=NF 1643 : SEGINI,LV1 1644 : C Boucle sur les cellules existantes 1645 : DO I=1,IPT6.NUM(/2) 1646 : NRCI=IPT6.NUM(1,I) 1646 : CALL ACCTAB(MTABC,'POINT ',0,0,' ',0,NRCI, 1647 : & 'TABLE ',0,0,0,0,MTABI1) 1647 : CALL ACCTAB(MTABI1,'MOT ',0,0,'FACS',0,0, 1648 : & 'LISTENTI',0,0,0,0,LFI1) 1649 : C Boucle sur les faces de la cellule NRCI 1650 : DO J=1,LFI1.LECT(/1) 1651 : NFJ=LFI1.LECT(J) 1652 : C Si cette face a deja ete regroupee, on examine la face 1653 : C suivante 1654 : IF (LV1.LECT(NFJ).NE.0) GOTO 9 1655 : C Sinon, on la considere comme englobante 1656 : LV1.LECT(NFJ)=NFJ 1657 : C Mais s'il s'agit d'une face interne, pas besoin de 1658 : C checher des faces coplanaires, on examine la suivante 1659 : IF (LV5.LECT(NFJ).NE.0) GOTO 9 1660 : C A partir d'ici, on va rechercher les autres faces de la 1661 : C cellule NRCI qui sont coplanaires a NFJ 1662 : C Ces faces sont listees dans LV2, utilisee par la suite 1663 : C pour regrouper les maillages des faces 1664 : JG=0 1665 : SEGINI,LV2 1665 : CALL ACCTAB(MTABF,'ENTIER ',NFJ,0,' ',0,0, 1666 : & 'TABLE ',0,0,0,0,MTABF1) 1666 : CALL ACCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 1667 : & 'MAILLAGE',0,0,0,0,IPTF1) 1668 : C Determination de l'equation du plan NFJ avec 3 points 1669 : C A, B et C non alignes 1670 : SEGINI,IPT3=IPTF1 1670 : CALL CHANGE(IPT3,1) 1671 : NRA=IPT3.NUM(1,1) 1672 : NRB=IPT3.NUM(1,2) 1673 : C Boucle pour trouver un point C non aligne avec les points 1674 : C A et B 1675 : DO K=3,IPT3.NUM(/2) 1676 : NRC=IPT3.NUM(1,K) 1677 : XA=XCOOR((NRA-1)*IDIMP1+1) 1678 : YA=XCOOR((NRA-1)*IDIMP1+2) 1679 : ZA=XCOOR((NRA-1)*IDIMP1+3) 1680 : XB=XCOOR((NRB-1)*IDIMP1+1) 1681 : YB=XCOOR((NRB-1)*IDIMP1+2) 1682 : ZB=XCOOR((NRB-1)*IDIMP1+3) 1683 : XC=XCOOR((NRC-1)*IDIMP1+1) 1684 : YC=XCOOR((NRC-1)*IDIMP1+2) 1685 : ZC=XCOOR((NRC-1)*IDIMP1+3) 1686 : C Calcul d'un vecteur normal au plan (ABC) 1687 : XN=((YB-YA)*(ZC-ZA))-((ZB-ZA)*(YC-YA)) 1688 : YN=((ZB-ZA)*(XC-XA))-((XB-XA)*(ZC-ZA)) 1689 : ZN=((XB-XA)*(YC-YA))-((YB-YA)*(XC-XA)) 1690 : DN=SQRT((XN**2)+(YN**2)+(ZN**2)) 1691 : IF (DN.GE.ZERO) GOTO 14 1692 : ENDDO 1693 : C Si l'on passe ici, c'est que l'on a pas trouve 3 points 1694 : C suffisament espaces pour definir le plan 1695 : C WRITE(6,*) ' > Impossible de determiner l''equation du plan' 1696 : C WRITE(6,*) ' points alignes (2)' 1697 : INTERR(1)=NFJ 1698 : INTERR(2)=NRCI 1698 : CALL ERREUR(1035) 1699 : GOTO 999 1700 : 14 CONTINUE 1701 : XN=XN/DN 1702 : YN=YN/DN 1703 : ZN=ZN/DN 1704 : W=-1D0*((XN*XA)+(YN*YA)+(ZN*ZA)) 1705 : C Boucle sur les faces suivantes 1706 : MMTABF1=MTABF1 1707 : IIPTF1=IPTF1 1708 : DO K=J+1,LFI1.LECT(/1) 1709 : NFK=LFI1.LECT(K) 1710 : C S'il s'agit d'une face interne, on ne regarde pas 1711 : C la coplaneite et on examine la face suivante 1712 : IF (LV5.LECT(NFK).NE.0) GOTO 8 1713 : C Si la face NFK a ete traitee (ratachee a une face ou 1714 : C bien considere comme face englobante), on recupere le 1715 : C numero de ratachement 1716 : IF (LV1.LECT(NFK).NE.0) THEN 1717 : NFK=LV1.LECT(NFK) 1718 : ENDIF 1719 : C Si la face NFK est deja reliee a la face NFJ, 1720 : C on ne refait pas le travail 1721 : IF (NFK.EQ.NFJ) GOTO 8 1722 : C Acces au maillage de la face NFK 1722 : CALL ACCTAB(MTABF,'ENTIER ',NFK,0,' ',0,0, 1723 : & 'TABLE ',0,0,0,0,MTABF1) 1723 : CALL ACCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 1724 : & 'MAILLAGE',0,0,0,0,IPTF1) 1725 : C Verification que les points de la face NFK sont dans 1726 : C le plan ABC 1727 : SEGINI,IPT4=IPTF1 1727 : CALL CHANGE(IPT4,1) 1728 : DO L=1,IPT4.NUM(/2) 1729 : NRL=IPT4.NUM(1,L) 1730 : XL=XCOOR((NRL-1)*IDIMP1+1) 1731 : YL=XCOOR((NRL-1)*IDIMP1+2) 1732 : ZL=XCOOR((NRL-1)*IDIMP1+3) 1733 : DTEST=(XN*XL)+(YN*YL)+(ZN*ZL)+W 1734 : IF (ABS(DTEST).GE.ZERO2) GOTO 8 1735 : ENDDO 1736 : C Si on passe ici, c'est que la face NFK est coplanaire 1737 : C a la face NFJ, on le note dans LV2 et on modifie 1738 : C l'indice NFK de LV1 pour rediriger vers la face NFJ 1739 : JG=LV2.LECT(/1)+1 1740 : SEGADJ,LV2 1741 : LV2.LECT(JG)=NFK 1742 : LV1.LECT(NFK)=NFJ 1743 : 8 CONTINUE 1744 : ENDDO 1745 : MTABF1=MMTABF1 1746 : IPTF1=IIPTF1 1747 : C Si on a trouve aucune face coplanaire a NFJ, on itere 1748 : IF (LV2.LECT(/1).EQ.0) GOTO 9 1749 : C On a determine la liste LV2 des faces coplanaires la face 1750 : C NFJ, on incorpore les aretes dans la table MTABF1 dediee 1751 : C a la face NFJ 1752 : MMTABF1=MTABF1 1752 : CALL ACCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, 1753 : & 'LISTENTI',0,0,0,0,LA1) 1754 : C Boucle sur les faces coplanaires a la face NFJ 1755 : DO K=1,LV2.LECT(/1) 1756 : NFK=LV2.LECT(K) 1756 : CALL ACCTAB(MTABF,'ENTIER ',NFK,0,' ',0,0, 1757 : & 'TABLE ',0,0,0,0,MTABF1) 1757 : CALL ACCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, 1758 : & 'LISTENTI',0,0,0,0,LA2) 1759 : C Boucle sur les aretes de cette face NFK 1760 : DO L=1,LA2.LECT(/1) 1761 : NAL=LA2.LECT(L) 1762 : C Si l'arete NAL est presente dans LA1, alors elle 1763 : C est supprimee de la face NFJ car elle est commune a 1764 : C deux faces coplanaires 1765 : DO M=1,LA1.LECT(/1) 1766 : IF ((LA1.LECT(M)).EQ.NAL) THEN 1766 : CALL ENLEV2(LA1,M,LA3) 1767 : SEGACT,LA1,LA3 1768 : LA1=LA3 1769 : GOTO 13 1770 : ENDIF 1771 : ENDDO 1772 : C Si on passe ici, c'est que l'arete NAL est absente 1773 : C de LA1, donc elle est ajoutee et contribue ainsi a 1774 : C la face NFJ 1775 : JG=LA1.LECT(/1)+1 1776 : SEGADJ,LA1 1777 : LA1.LECT(JG)=NAL 1778 : 13 CONTINUE 1779 : ENDDO 1780 : ENDDO 1781 : MTABF1=MMTABF1 1782 : C Ecriture de la nouvelle liste d'aretes pour la face NFJ 1782 : CALL ECCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, 1783 : & 'LISTENTI',0,0,0,0,LA1) 1784 : C Reconstruction du maillage de la face NFJ 1785 : IIPTF1=IPTF1 1785 : CALL ACCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 1786 : & 'MAILLAGE',0,0,0,0,IPTF1) 1787 : NBNN=IPTF1.NUM(/1) 1788 : NBELEM=LA1.LECT(/1) 1789 : SEGINI,IPTF1 1790 : IPTF1.ITYPEL=2 1791 : DO K=1,LA1.LECT(/1) 1792 : NAK=LA1.LECT(K) 1792 : CALL ACCTAB(MTABA,'ENTIER ',NAK,0,' ',0,0, 1793 : & 'MAILLAGE',0,0,0,0,IPTA1) 1794 : IPTF1.NUM(1,K)=IPTA1.NUM(1,1) 1795 : IPTF1.NUM(2,K)=IPTA1.NUM(2,1) 1796 : ENDDO 1796 : CALL ECCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 1797 : & 'MAILLAGE',0,0,0,0,IPTF1) 1798 : IPTF1=IIPTF1 1799 : 9 CONTINUE 1800 : ENDDO 1801 : ENDDO 1802 : C 1803 : C Ainsi, les faces englobees sont identifiees dans LV1, 1804 : C les faces englobantes aussi et leur maillage et liste 1805 : C d'aretes sont mises a jour dans MTABF1 1806 : C Il faut maintenant re-numeroter les faces et ne plus faire 1807 : C apparaitre les faces englobees 1808 : C 1809 : C Re-ecriture de la table de sortie pour ne plus faire figurer 1810 : C les faces/aretes supprimees 1810 : CALL CRTABL(MTAB10) 1810 : CALL CRTABL(MTABCC) 1810 : CALL CRTABL(MTABFF) 1810 : CALL CRTABL(MTABAA) 1811 : NBARETES=0 1812 : NBFACES=0 1813 : C listes faisant la correspondance entre les anciens et les 1814 : C nouveaux elements numerotes : 1815 : C LV3(n)=m l'ancienne arete numero n porte maintenant le numero m 1816 : C LV4(n)=m l'ancienne face numero n porte maintenant le numero m 1817 : C n doit correpondre a une face englobante 1818 : C si n correspond a une face "englobee" alors LV4(n)=0 1819 : JG=IPTG.NUM(/2) 1820 : SEGINI,LV3 1821 : JG=NF 1822 : SEGINI,LV4 1823 : C initialisation du nouveau maillage global des aretes 1824 : NBNN=2 1825 : NBELEM=IPTG.NUM(/2) 1826 : SEGINI,IPTG0 1827 : IPTG0.ITYPEL=2 1827 : CALL ECCTAB(MTAB10,'MOT ',0,0,'VISU',0,0, 1828 : & 'MAILLAGE',0,0,0,0,IPTG0) 1829 : C initialisation de la table des cellules MTABCC 1829 : CALL ECCTAB(MTAB10,'MOT ',0,0,'CELL',0,0, 1830 : & 'TABLE ',0,0,0,0,MTABCC) 1831 : C initialisation de la table des faces MTABF0 1831 : CALL ECCTAB(MTAB10,'MOT ',0,0,'FACS',0,0, 1832 : & 'TABLE ',0,0,0,0,MTABFF) 1833 : C initialisation de la table des aretes MTABA0 1833 : CALL ECCTAB(MTAB10,'MOT ',0,0,'ARTS',0,0, 1834 : & 'TABLE ',0,0,0,0,MTABAA) 1835 : C boucle sur les cellules 1836 : DO I=1,IPT6.NUM(/2) 1837 : NRCI=IPT6.NUM(1,I) 1838 : C initialisation de la table MTABI0 de la cellule NRCI 1838 : CALL CRTABL(MTABI0) 1838 : CALL ECCTAB(MTABCC,'POINT ',0,0,' ',0,NRCI, 1839 : & 'TABLE ',0,0,0,0,MTABI0) 1840 : C initialisation du maillage IPTI0 de la cellule NRCI 1841 : NBNN=2 1842 : NBELEM=0 1843 : SEGINI,IPTI0 1844 : IPTI0.ITYPEL=2 1844 : CALL ECCTAB(MTABI0,'MOT ',0,0,'VISU',0,0, 1845 : & 'MAILLAGE',0,0,0,0,IPTI0) 1846 : C initialisation de la liste LFI0 des faces de NRCI 1847 : JG=0 1848 : SEGINI,LFI0 1848 : CALL ECCTAB(MTABI0,'MOT ',0,0,'FACS',0,0, 1849 : & 'LISTENTI',0,0,0,0,LFI0) 1850 : C ecriture du maillage IPTVI0 des voisins de NRCI (inchange) 1850 : CALL ACCTAB(MTABC,'POINT ',0,0,' ',0,NRCI, 1851 : & 'TABLE ',0,0,0,0,MTABI1) 1851 : CALL ACCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, 1852 : & 'MAILLAGE',0,0,0,0,IPTVI1) 1852 : CALL ECCTAB(MTABI0,'MOT ',0,0,'VOIS',0,0, 1853 : & 'MAILLAGE',0,0,0,0,IPTVI1) 1854 : C boucle sur les faces de la cellule NRCI 1854 : CALL ACCTAB(MTABI1,'MOT ',0,0,'FACS',0,0, 1855 : & 'LISTENTI',0,0,0,0,LFI1) 1856 : DO J=1,LFI1.LECT(/1) 1857 : NFJ=LFI1.LECT(J) 1858 : C s'il s'agit il d'une face englobee, on passe a la face 1859 : C suivante 1860 : IF (LV1.LECT(NFJ).NE.NFJ) GOTO 18 1861 : C la face est elle nouvelle pour cette re-numerotation ? 1862 : NNFJ=LV4.LECT(NFJ) 1863 : BOOL1=NNFJ.EQ.0 1864 : IF (BOOL1) THEN 1865 : C NFJ n'a pas ete re-numerotee, donc on le fait ici 1866 : C NFJ est re-numerotee en NNFJ 1867 : NBFACES=NBFACES+1 1868 : NNFJ=NBFACES 1869 : LV4.LECT(NFJ)=NNFJ 1870 : C creation de la table MTABF0 de la face NNFJ 1870 : CALL CRTABL(MTABF0) 1870 : CALL ECCTAB(MTABFF,'ENTIER ',NNFJ,0,' ',0,0, 1871 : & 'TABLE ',0,0,0,0,MTABF0) 1872 : ELSE 1873 : C NFJ correspond a une face deja re-numerotee, on va 1874 : C chercher sa table 1874 : CALL ACCTAB(MTABFF,'ENTIER ',NNFJ,0,' ',0,0, 1875 : & 'TABLE ',0,0,0,0,MTABF0) 1876 : ENDIF 1877 : C ajout de cette face dans la liste LFI0 de la cellule NRCI 1878 : JG=LFI0.LECT(/1)+1 1879 : SEGADJ,LFI0 1880 : LFI0.LECT(JG)=NNFJ 1881 : C acces a la liste des aretes de cette face 1881 : CALL ACCTAB(MTABF,'ENTIER ',NFJ,0,' ',0,0, 1882 : & 'TABLE ',0,0,0,0,MTABF1) 1882 : CALL ACCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 1883 : & 'MAILLAGE',0,0,0,0,IPTF1) 1883 : CALL ACCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, 1884 : & 'LISTENTI',0,0,0,0,LA1) 1885 : C attention, la liste des aretes LA1 contient les anciens 1886 : C numeros d'arete 1887 : IF (BOOL1) THEN 1888 : C si la face NFJ vient d'etre re-numerotee NNFJ, alors 1889 : C on ecrit le maillage de la face dans MTABF0 1889 : CALL ECCTAB(MTABF0,'MOT ',0,0,'VISU',0,0, 1890 : & 'MAILLAGE',0,0,0,0,IPTF1) 1891 : C initialisation de la liste LA0 des aretes de la face 1892 : C NNFJ que l'on va remplir de suite 1893 : JG=LA1.LECT(/1) 1894 : SEGINI,LA0 1895 : NA0=0 1895 : CALL ECCTAB(MTABF0,'MOT ',0,0,'ARTS',0,0, 1896 : & 'LISTENTI',0,0,0,0,LA0) 1897 : ENDIF 1898 : C boucle sur les aretes de la face 1899 : DO K=1,LA1.LECT(/1) 1900 : C ancien numero d'arete 1901 : NAK=LA1.LECT(K) 1901 : CALL ACCTAB(MTABA,'ENTIER ',NAK,0,' ',0,0, 1902 : & 'MAILLAGE',0,0,0,0,IPTA1) 1903 : NRA=IPTA1.NUM(1,1) 1904 : NRB=IPTA1.NUM(2,1) 1905 : C acces au nouveau numero de cette arete 1906 : NNAK=LV3.LECT(NAK) 1907 : C s'il s'agit d'une nouvelle face 1908 : IF (BOOL1) THEN 1909 : C s'agit il d'une arete non re-numerotee ? 1910 : IF (NNAK.EQ.0) THEN 1911 : NBARETES=NBARETES+1 1912 : NNAK=NBARETES 1913 : LV3.LECT(NAK)=NNAK 1914 : C ecriture du maillage IPTA0 (identique a IPTA1) 1915 : C de l'arete NNAK 1915 : CALL ECCTAB(MTABAA,'ENTIER ',NNAK,0,' ',0,0, 1916 : & 'MAILLAGE',0,0,0,0,IPTA1) 1917 : C ajout cette arete dans le maillage global IPTG0 1918 : IPTG0.NUM(1,NBARETES)=NRA 1919 : IPTG0.NUM(2,NBARETES)=NRB 1920 : ENDIF 1921 : C ajout du numero d'arete dans la liste LA0 1922 : NA0=NA0+1 1923 : LA0.LECT(NA0)=NNAK 1924 : ENDIF 1925 : C ajout de l'arete dans le maillage IPTI0 de la cellule 1926 : C si elle n'y est pas deja 1927 : BOOL2=.FALSE. 1928 : DO L=1,IPTI0.NUM(/2) 1929 : NRL1=IPTI0.NUM(1,L) 1930 : NRL2=IPTI0.NUM(2,L) 1931 : IF (((NRA.EQ.NRL1).AND.(NRB.EQ.NRL2)).OR. 1932 : & ((NRA.EQ.NRL2).AND.(NRB.EQ.NRL1))) THEN 1933 : BOOL2=.TRUE. 1934 : GOTO 23 1935 : ENDIF 1936 : ENDDO 1937 : 23 CONTINUE 1938 : IF (.NOT.BOOL2) THEN 1939 : NBNN=IPTI0.NUM(/1) 1940 : NBELEM=IPTI0.NUM(/2)+1 1941 : SEGADJ,IPTI0 1942 : IPTI0.NUM(1,NBELEM)=NRA 1943 : IPTI0.NUM(2,NBELEM)=NRB 1944 : ENDIF 1945 : ENDDO 1946 : 18 CONTINUE 1947 : ENDDO 1948 : ENDDO 1949 : C Ajustement de IPTG0 1950 : NBNN=IPTG0.NUM(/1) 1951 : NBELEM=NBARETES 1952 : SEGADJ,IPTG0 1953 : C Mise a jour de ITFAC1 avec cette nouvelle numerotation 1954 : DO I=1,LV5.LECT(/1) 1955 : C cellules definissant l'ancienne face I 1956 : NCI=LV5.LECT(I) 1957 : NCJ=LV6.LECT(I) 1958 : IF (NCI.NE.0) THEN 1959 : C ancienne face englobante de I 1960 : NFOLD=LV1.LECT(I) 1961 : C nouveau numero de cette face 1962 : NFNEW=LV4.LECT(NFOLD) 1963 : ITFAC1.LACT(NCI,NCJ)=NFNEW 1964 : ITFAC1.LACT(NCJ,NCI)=NFNEW 1965 : ENDIF 1966 : ENDDO 1967 : C Ecrasement des tables et maillages precedents 1968 : MTAB11=MTAB10 1969 : IPTG=IPTG0 1970 : MTABC=MTABCC 1971 : MTABF=MTABFF 1972 : MTABA=MTABAA 1973 : NF=NBFACES 1974 : ENDIF 1975 : C 1976 : C 1977 : C 1978 : C----------------------------------------------------------------------C 1979 : C PARTIE 3 C 1980 : C Construction des cellules entieres en ajoutant les aretes internes C 1981 : C----------------------------------------------------------------------C 1982 : C 1983 : IF (IDIM.EQ.2) THEN 1984 : C IPT3 : maillage des cellules non coupees de MTAB1 1984 : CALL ACCTAB(MTAB1,'MOT ',0,0,'VISU',0,0, 1985 : & 'MAILLAGE',0,0,0,0,IPT3) 1986 : C ITARE1 : tableau donnant les cellules appuyees sur chaque 1987 : C arete de IPT3 1987 : CALL ACCTAB(MTAB1,'MOT ',0,0,'TAR',0,0, 1988 : & 'MVORO ',0,0,0,0,ITARE1) 1989 : SEGACT,IPT3,ITARE1 1990 : C MTABCC : table des cellules coupees fixee 1991 : SEGINI,MTABCC=MTABC 1992 : C Balayage des aretes non coupees IPT3 et ajout dans MTAB11 si 1993 : C elles sont contenues dans l'enveloppe 1994 : DO I=1,IPT3.NUM(/2) 1995 : C LV1 : liste des cellules appuyees sur l'arete I 1996 : JG=ITARE1.LACT(/1) 1997 : SEGINI,LV1 1998 : DO J=1,ITARE1.LACT(/1) 1999 : NCJ=ITARE1.LACT(J,I) 2000 : IF (NCJ.EQ.0) THEN 2001 : JG=J-1 2002 : GOTO 36 2003 : ENDIF 2004 : LV1.LECT(J)=NCJ 2005 : ENDDO 2006 : 36 CONTINUE 2007 : SEGADJ,LV1 2008 : C Si l'arete s'appuie sur 1 cellule, elle est forcement 2009 : C externe, donc on ne l'examine meme pas ! 2010 : IF (LV1.LECT(/1).LT.1) THEN 2011 : GOTO 39 2012 : ENDIF 2013 : C Extremites de l'arete non coupee 2014 : NPA=IPT3.NUM(1,I) 2015 : NPB=IPT3.NUM(2,I) 2016 : C LV2 : Liste des noeuds situes sur l'arete de A vers B 2017 : JG=2 2018 : SEGINI,LV2 2019 : LV2.LECT(1)=NPA 2020 : LV2.LECT(2)=NPB 2021 : C LV3 : Indique si les noeuds de LV2 sont a inclure dans le 2022 : C maillage de l'arete coupee a sortir 2023 : JG=LV2.LECT(/1) 2024 : SEGINI,LV3 2025 : IA=0 2025 : CALL DEDANS(NPA,IPT2,BOOL1) 2026 : IF (BOOL1) IA=1 2027 : LV3.LECT(1)=IA 2028 : IB=0 2028 : CALL DEDANS(NPB,IPT2,BOOL1) 2029 : IF (BOOL1) IB=1 2030 : LV3.LECT(2)=IB 2031 : C LZ1 : Coordonnees barycentriques des noeuds de LV2 2032 : C 0 pour A, 1 pour B 2033 : JG=LV2.LECT(/1) 2034 : SEGINI,LZ1 2035 : LZ1.PROG(1)=0D0 2036 : LZ1.PROG(2)=1D0 2037 : C La premiere cellule de LV1 est elle une cellule coupee ? 2038 : BOOL1=.FALSE. 2039 : NC1=LV1.LECT(1) 2040 : NRC1=IPT1.NUM(1,NC1) 2040 : CALL ECROBJ('POINT ',NRC1) 2040 : CALL ECROBJ('TABLE ',MTABCC) 2040 : CALL EXIS 2040 : CALL LIRLOG(BOOL1,1,IRETOU) 2041 : IF (BOOL1) THEN 2042 : C La cellule C1 est coupee, donc l'arete I est peut etre 2043 : C coupee. Elle peut etre en plusieurs morceaux, on etablit 2044 : C la liste des noeuds definissant cette arete grace a 2045 : C LV2 et LV3 en ordonnant de PA vers PB 2046 : C Examen des noeuds de MTAB11.CELL.NRC1.VISU alignes 2047 : C sur la droite PA PB 2047 : CALL ACCTAB(MTABC,'POINT ',0,0,' ',0,NRC1, 2048 : & 'TABLE ',0,0,0,0,MTABI1) 2048 : CALL ACCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 2049 : & 'MAILLAGE',0,0,0,0,IPTI1) 2050 : SEGINI,IPT4=IPTI1 2050 : CALL CHANGE(IPT4,1) 2051 : DO J=1,IPT4.NUM(/2) 2052 : NPC=IPT4.NUM(1,J) 2053 : IF ((NPC.EQ.NPA).OR.(NPC.EQ.NPB)) GOTO 37 2054 : C C est il dans le segment [AB] ? 2055 : XA=XCOOR((NPA-1)*IDIMP1+1) 2056 : YA=XCOOR((NPA-1)*IDIMP1+2) 2057 : XB=XCOOR((NPB-1)*IDIMP1+1) 2058 : YB=XCOOR((NPB-1)*IDIMP1+2) 2059 : XC=XCOOR((NPC-1)*IDIMP1+1) 2060 : YC=XCOOR((NPC-1)*IDIMP1+2) 2061 : C Test sur les distances, si AC+CB=AB alors A,C,B 2062 : C sont alignes et se suivent 2063 : AB=SQRT(((XB-XA)**2)+((YB-YA)**2)) 2064 : AC=SQRT(((XC-XA)**2)+((YC-YA)**2)) 2065 : CB=SQRT(((XB-XC)**2)+((YB-YC)**2)) 2066 : c----------------------------------------------------------------------- 2067 : IF ((ABS((AC+CB)-AB))/AB.LT.ZERO) THEN 2068 : c IF (ABS((AC+CB)-AB).LT.ZERO) THEN 2069 : c----------------------------------------------------------------------- 2070 : C Rangement du point C dans LV2 et LV3 selon sa 2071 : C coordonnee barycentrique LC 2072 : GC=AC/AB 2073 : DO K=1,LZ1.PROG(/1)-1 2074 : G1=LZ1.PROG(K) 2075 : G2=LZ1.PROG(K+1) 2076 : IF ((G1.LT.GC).AND.(GC.LT.G2)) THEN 2077 : IEME=K+1 2078 : GOTO 38 2079 : ENDIF 2080 : ENDDO 2081 : 38 CONTINUE 2081 : CALL INSER1(LZ1,IEME,GC) 2082 : SEGACT,LZ1 2082 : CALL INSER2(LV2,IEME,NPC) 2083 : SEGACT,LV2 2084 : IC=0 2084 : CALL DEDANS(NPC,IPT2,BOOL2) 2085 : IF (BOOL2) IC=1 2085 : CALL INSER2(LV3,IEME,IC) 2086 : SEGACT,LV3 2087 : ENDIF 2088 : 37 CONTINUE 2089 : ENDDO 2090 : ELSE 2091 : C La cellule C1 n'est pas coupee, donc l'arete I n'est pas 2092 : C coupee. Est elle a inclure ou bien a exclure ? 2093 : C Test si le centre NRC1 est dans l'enveloppe 2093 : CALL DEDANS(NRC1,IPT2,BOOL2) 2094 : IF (BOOL2) THEN 2095 : C L'arete (LV2,LV3) est a inclure entierement et ne 2096 : C contient que les points A et B 2097 : LV3.LECT(1)=1 2098 : LV3.LECT(2)=1 2099 : ELSE 2100 : C L'arete LV2 est a exclure entierement, on itere 2101 : GOTO 39 2102 : ENDIF 2103 : ENDIF 2104 : C Test si l'arete decrite par LV2 et LV3 est externe 2105 : BOOL1=.FALSE. 2106 : DO J=1,LV3.LECT(/1) 2107 : IF (LV3.LECT(J).NE.0) BOOL1=.TRUE. 2108 : ENDDO 2109 : IF (.NOT.BOOL1) THEN 2110 : GOTO 39 2111 : ENDIF 2112 : C Maillage IPT4 de l'arete suivant les listes LV2 et LV3 2113 : C IPT4 peut contenir plusieurs elements si cette arete est 2114 : C en plusieurs morceaux 2115 : NBNN=2 2116 : NBELEM=0 2117 : SEGINI,IPT4 2118 : IPT4.ITYPEL=2 2119 : BOOL1=.TRUE. 2120 : DO J=1,(LV2.LECT(/1)-1) 2121 : NA=LV2.LECT(J) 2122 : IA=LV3.LECT(J) 2123 : NB=LV2.LECT(J+1) 2124 : IB=LV3.LECT(J+1) 2125 : IF ((IA.EQ.1).AND.(IB.EQ.1).AND.BOOL1) THEN 2126 : C Ajout d'un element de maillage pour l'arete IPT4 2127 : NBNN=IPT4.NUM(/1) 2128 : NBELEM=IPT4.NUM(/2)+1 2129 : SEGADJ,IPT4 2130 : IPT4.NUM(1,NBELEM)=NA 2131 : IPT4.NUM(2,NBELEM)=NB 2132 : BOOL1=.FALSE. 2133 : ELSE 2134 : BOOL1=.TRUE. 2135 : ENDIF 2136 : ENDDO 2137 : C Recuperation du maillage, de la liste des aretes et du 2138 : C maillage des voisins pour les 2 cellules concernees 2139 : NC1=LV1.LECT(1) 2140 : NC2=LV1.LECT(2) 2141 : NRC1=IPT1.NUM(1,NC1) 2142 : NRC2=IPT1.NUM(1,NC2) 2142 : CALL ECROBJ('POINT ',NRC1) 2142 : CALL ECROBJ('TABLE ',MTABC) 2142 : CALL EXIS 2142 : CALL LIRLOG(BOOL1,1,IRETOU) 2143 : IF (.NOT.BOOL1) THEN 2143 : CALL CRTABL(MTABI1) 2143 : CALL ECCTAB(MTABC,'POINT ',0,0,' ',0,NRC1, 2144 : & 'TABLE ',0,0,0,0,MTABI1) 2145 : NBNN=2 2146 : NBELEM=0 2147 : SEGINI,IPTI1 2148 : IPTI1.ITYPEL=2 2148 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 2149 : & 'MAILLAGE',0,0,0,0,IPTI1) 2150 : JG=0 2151 : SEGINI,LA1 2151 : CALL ECCTAB(MTABI1,'MOT ',0,0,'ARTS',0,0, 2152 : & 'LISTENTI',0,0,0,0,LA1) 2153 : NBNN=1 2154 : NBELEM=0 2155 : SEGINI,IPTVI1 2156 : IPTVI1.ITYPEL=1 2156 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, 2157 : & 'MAILLAGE',0,0,0,0,IPTVI1) 2158 : ENDIF 2158 : CALL ACCTAB(MTABC,'POINT ',0,0,' ',0,NRC1, 2159 : & 'TABLE ',0,0,0,0,MTABI1) 2159 : CALL ACCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 2160 : & 'MAILLAGE',0,0,0,0,IPTI1) 2160 : CALL ACCTAB(MTABI1,'MOT ',0,0,'ARTS',0,0, 2161 : & 'LISTENTI',0,0,0,0,LA1) 2161 : CALL ACCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, 2162 : & 'MAILLAGE',0,0,0,0,IPTVI1) 2162 : CALL ECROBJ('POINT ',NRC2) 2162 : CALL ECROBJ('TABLE ',MTABC) 2162 : CALL EXIS 2162 : CALL LIRLOG(BOOL1,1,IRETOU) 2163 : IF (.NOT.BOOL1) THEN 2163 : CALL CRTABL(MTABI2) 2163 : CALL ECCTAB(MTABC,'POINT ',0,0,' ',0,NRC2, 2164 : & 'TABLE ',0,0,0,0,MTABI2) 2165 : NBNN=2 2166 : NBELEM=0 2167 : SEGINI,IPTI2 2168 : IPTI2.ITYPEL=2 2168 : CALL ECCTAB(MTABI2,'MOT ',0,0,'VISU',0,0, 2169 : & 'MAILLAGE',0,0,0,0,IPTI2) 2170 : JG=0 2171 : SEGINI,LA2 2171 : CALL ECCTAB(MTABI2,'MOT ',0,0,'ARTS',0,0, 2172 : & 'LISTENTI',0,0,0,0,LA2) 2173 : NBNN=1 2174 : NBELEM=0 2175 : SEGINI,IPTVI2 2176 : IPTVI2.ITYPEL=1 2176 : CALL ECCTAB(MTABI2,'MOT ',0,0,'VOIS',0,0, 2177 : & 'MAILLAGE',0,0,0,0,IPTVI2) 2178 : ENDIF 2178 : CALL ACCTAB(MTABC,'POINT ',0,0,' ',0,NRC2, 2179 : & 'TABLE ',0,0,0,0,MTABI2) 2179 : CALL ACCTAB(MTABI2,'MOT ',0,0,'VISU',0,0, 2180 : & 'MAILLAGE',0,0,0,0,IPTI2) 2180 : CALL ACCTAB(MTABI2,'MOT ',0,0,'ARTS',0,0, 2181 : & 'LISTENTI',0,0,0,0,LA2) 2181 : CALL ACCTAB(MTABI2,'MOT ',0,0,'VOIS',0,0, 2182 : & 'MAILLAGE',0,0,0,0,IPTVI2) 2183 : C Rangement des elements de IPT4 dans les indices de la table 2184 : C de sortie TAB11 2185 : C --> Ajout de(s) arete(s) dans TAB11.VISU et TAB11.ARTS 2186 : C ainsi que dans les tables des 2 cellules NRC1 et NRC2 2187 : NBELEG0=IPTG.NUM(/2) 2188 : NBNN=IPTG.NUM(/1) 2189 : NBELEM=NBELEG0+IPT4.NUM(/2) 2190 : SEGADJ,IPTG 2191 : DO J=1,IPT4.NUM(/2) 2192 : C pour IPTG 2193 : JJ=J+NBELEG0 2194 : IPTG.NUM(1,JJ)=IPT4.NUM(1,J) 2195 : IPTG.NUM(2,JJ)=IPT4.NUM(2,J) 2196 : C pour la table des aretes MTABA 2197 : NBNN=2 2198 : NBELEM=1 2199 : SEGINI,IPT5 2200 : IPT5.ITYPEL=2 2201 : IPT5.NUM(1,1)=IPT4.NUM(1,J) 2202 : IPT5.NUM(2,1)=IPT4.NUM(2,J) 2202 : CALL ECCTAB(MTABA,'ENTIER ',JJ,0,' ',0,0, 2203 : & 'MAILLAGE',0,0,0,0,IPT5) 2204 : C indice VISU pour la cellule NRC1 2205 : NBNN=IPTI1.NUM(/1) 2206 : NBELEM=IPTI1.NUM(/2)+1 2207 : SEGADJ,IPTI1 2208 : IPTI1.NUM(1,NBELEM)=IPT4.NUM(1,J) 2209 : IPTI1.NUM(2,NBELEM)=IPT4.NUM(2,J) 2210 : C indice ARTS pour la cellule NRC1 2211 : JG=LA1.LECT(/1)+1 2212 : SEGADJ,LA1 2213 : LA1.LECT(JG)=JJ 2214 : C indice VOIS pour la cellule NRC1 2215 : NBNN=IPTVI1.NUM(/1) 2216 : NBELEM=IPTVI1.NUM(/2)+1 2217 : SEGADJ,IPTVI1 2218 : IPTVI1.NUM(1,NBELEM)=NRC2 2219 : C indice VISU pour la cellule NRC2 2220 : NBNN=IPTI2.NUM(/1) 2221 : NBELEM=IPTI2.NUM(/2)+1 2222 : SEGADJ,IPTI2 2223 : IPTI2.NUM(1,NBELEM)=IPT4.NUM(1,J) 2224 : IPTI2.NUM(2,NBELEM)=IPT4.NUM(2,J) 2225 : C indice ARTS pour la cellule NRC2 2226 : JG=LA2.LECT(/1)+1 2227 : SEGADJ,LA2 2228 : LA2.LECT(JG)=JJ 2229 : C indice VOIS pour la cellule NRC2 2230 : NBNN=IPTVI2.NUM(/1) 2231 : NBELEM=IPTVI2.NUM(/2)+1 2232 : SEGADJ,IPTVI2 2233 : IPTVI2.NUM(1,NBELEM)=NRC1 2234 : ENDDO 2235 : 39 CONTINUE 2236 : ENDDO 2237 : C 2238 : C 2239 : ELSEIF (IDIM.EQ.3) THEN 2240 : NF0=NF 2241 : C IPT3 : maillage MAV cellules non coupees de MTAB1 2241 : CALL ACCTAB(MTAB1,'MOT ',0,0,'VISU',0,0, 2242 : & 'MAILLAGE',0,0,0,0,IPT3) 2243 : C ITARE1 : tableau donnant les cellules appuyees sur chaque 2244 : C arete de IPT3 2244 : CALL ACCTAB(MTAB1,'MOT ',0,0,'TAR',0,0, 2245 : & 'MVORO ',0,0,0,0,ITARE1) 2246 : SEGACT,IPT3,ITARE1 2247 : C MTABCC : table des cellules coupees fixee 2248 : SEGINI,MTABCC=MTABC 2249 : C Balayage des aretes non coupees IPT3 et ajout dans MTAB11 si 2250 : C elles sont contenues dans l'enveloppe 2251 : DO I=1,IPT3.NUM(/2) 2252 : C LV1 : liste des cellules appuyees sur l'arete I 2253 : JG=ITARE1.LACT(/1) 2254 : SEGINI,LV1 2255 : DO J=1,ITARE1.LACT(/1) 2256 : NCJ=ITARE1.LACT(J,I) 2257 : IF (NCJ.EQ.0) THEN 2258 : JG=J-1 2259 : GOTO 10 2260 : ENDIF 2261 : LV1.LECT(J)=NCJ 2262 : ENDDO 2263 : 10 CONTINUE 2264 : SEGADJ,LV1 2265 : C Si l'arete s'appuie sur moins de 3 cellules, elle est 2266 : C forcement externe, donc on ne l'examine meme pas ! 2267 : IF (LV1.LECT(/1).LT.3) THEN 2268 : GOTO 11 2269 : ENDIF 2270 : C Extremites de l'arete non coupee 2271 : NPA=IPT3.NUM(1,I) 2272 : NPB=IPT3.NUM(2,I) 2273 : C LV2 : Liste des noeuds situes sur l'arete de A vers B 2274 : JG=2 2275 : SEGINI,LV2 2276 : LV2.LECT(1)=NPA 2277 : LV2.LECT(2)=NPB 2278 : C LV3 : Indique si les noeuds de LV2 sont a inclure dans le 2279 : C maillage de l'arete coupee a sortir 2280 : JG=LV2.LECT(/1) 2281 : SEGINI,LV3 2282 : IA=0 2282 : CALL DEDANS(NPA,IPT2,BOOL1) 2283 : IF (BOOL1) IA=1 2284 : LV3.LECT(1)=IA 2285 : IB=0 2285 : CALL DEDANS(NPB,IPT2,BOOL1) 2286 : IF (BOOL1) IB=1 2287 : LV3.LECT(2)=IB 2288 : C LZ1 : Coordonnees barycentriques des noeuds de LV2 2289 : C 0 pour A, 1 pour B 2290 : JG=LV2.LECT(/1) 2291 : SEGINI,LZ1 2292 : LZ1.PROG(1)=0D0 2293 : LZ1.PROG(2)=1D0 2294 : C La premiere cellule de LV1 est elle une cellule coupee ? 2295 : BOOL1=.FALSE. 2296 : NC1=LV1.LECT(1) 2297 : NRC1=IPT1.NUM(1,NC1) 2297 : CALL ECROBJ('POINT ',NRC1) 2297 : CALL ECROBJ('TABLE ',MTABCC) 2297 : CALL EXIS 2297 : CALL LIRLOG(BOOL1,1,IRETOU) 2298 : IF (BOOL1) THEN 2299 : C La cellule C1 est coupee, donc l'arete I est peut etre 2300 : C coupee. Elle peut etre en plusieurs morceaux, on etablit 2301 : C la liste des noeuds definissant cette arete grace a 2302 : C LV2 et LV3 en ordonnant de PA vers PB 2303 : C Examen des noeuds de MTAB11.CELL.NRC1.VISU alignes 2304 : C sur la droite PA PB 2304 : CALL ACCTAB(MTABC,'POINT ',0,0,' ',0,NRC1, 2305 : & 'TABLE ',0,0,0,0,MTABI1) 2305 : CALL ACCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 2306 : & 'MAILLAGE',0,0,0,0,IPTI1) 2307 : SEGINI,IPT4=IPTI1 2307 : CALL CHANGE(IPT4,1) 2308 : DO J=1,IPT4.NUM(/2) 2309 : NPC=IPT4.NUM(1,J) 2310 : IF ((NPC.EQ.NPA).OR.(NPC.EQ.NPB)) GOTO 22 2311 : C C est il dans le segment [AB] ? 2312 : XA=XCOOR((NPA-1)*IDIMP1+1) 2313 : YA=XCOOR((NPA-1)*IDIMP1+2) 2314 : ZA=XCOOR((NPA-1)*IDIMP1+3) 2315 : XB=XCOOR((NPB-1)*IDIMP1+1) 2316 : YB=XCOOR((NPB-1)*IDIMP1+2) 2317 : ZB=XCOOR((NPB-1)*IDIMP1+3) 2318 : XC=XCOOR((NPC-1)*IDIMP1+1) 2319 : YC=XCOOR((NPC-1)*IDIMP1+2) 2320 : ZC=XCOOR((NPC-1)*IDIMP1+3) 2321 : C Test sur les distances, si AC+CB=AB alors A,C,B 2322 : C sont alignes et se suivent 2323 : AB=SQRT(((XB-XA)**2)+((YB-YA)**2)+((ZB-ZA)**2)) 2324 : AC=SQRT(((XC-XA)**2)+((YC-YA)**2)+((ZC-ZA)**2)) 2325 : CB=SQRT(((XB-XC)**2)+((YB-YC)**2)+((ZB-ZC)**2)) 2326 : c----------------------------------------------------------------------- 2327 : IF ((ABS((AC+CB)-AB))/AB.LT.ZERO) THEN 2328 : c IF (ABS((AC+CB)-AB).LT.ZERO) THEN 2329 : c----------------------------------------------------------------------- 2330 : C Rangement du point C dans LV2 et LV3 selon sa 2331 : C coordonnee barycentrique LC 2332 : GC=AC/AB 2333 : DO K=1,LZ1.PROG(/1)-1 2334 : G1=LZ1.PROG(K) 2335 : G2=LZ1.PROG(K+1) 2336 : IF ((G1.LT.GC).AND.(GC.LT.G2)) THEN 2337 : IEME=K+1 2338 : GOTO 12 2339 : ENDIF 2340 : ENDDO 2341 : 12 CONTINUE 2341 : CALL INSER1(LZ1,IEME,GC) 2342 : SEGACT,LZ1 2342 : CALL INSER2(LV2,IEME,NPC) 2343 : SEGACT,LV2 2344 : IC=0 2344 : CALL DEDANS(NPC,IPT2,BOOL2) 2345 : IF (BOOL2) IC=1 2345 : CALL INSER2(LV3,IEME,IC) 2346 : SEGACT,LV3 2347 : ENDIF 2348 : 22 CONTINUE 2349 : ENDDO 2350 : ELSE 2351 : C La cellule C1 n'est pas coupee, donc l'arete I n'est pas 2352 : C coupee. Est elle a inclure ou bien a exclure ? 2353 : C Test si le centre NRC1 est dans l'enveloppe 2353 : CALL DEDANS(NRC1,IPT2,BOOL2) 2354 : IF (BOOL2) THEN 2355 : C L'arete (LV2,LV3) est a inclure entierement et ne 2356 : C contient que les points A et B 2357 : LV3.LECT(1)=1 2358 : LV3.LECT(2)=1 2359 : ELSE 2360 : C L'arete LV2 est a exclure entierement, on itere 2361 : GOTO 11 2362 : ENDIF 2363 : ENDIF 2364 : C Test si l'arete decrite par LV2 et LV3 est externe 2365 : BOOL1=.FALSE. 2366 : DO J=1,LV3.LECT(/1) 2367 : IF (LV3.LECT(J).NE.0) BOOL1=.TRUE. 2368 : ENDDO 2369 : IF (.NOT.BOOL1) THEN 2370 : GOTO 11 2371 : ENDIF 2372 : C Maillage IPT4 de l'arete suivant les listes LV2 et LV3 2373 : C IPT4 peut contenir plusieurs elements si cette arete est 2374 : C en plusieurs morceaux 2375 : NBNN=2 2376 : NBELEM=0 2377 : SEGINI,IPT4 2378 : IPT4.ITYPEL=2 2379 : BOOL1=.TRUE. 2380 : DO J=1,(LV2.LECT(/1)-1) 2381 : NA=LV2.LECT(J) 2382 : IA=LV3.LECT(J) 2383 : NB=LV2.LECT(J+1) 2384 : IB=LV3.LECT(J+1) 2385 : IF ((IA.EQ.1).AND.(IB.EQ.1).AND.BOOL1) THEN 2386 : C Ajout d'un element de maillage pour l'arete IPT4 2387 : NBNN=IPT4.NUM(/1) 2388 : NBELEM=IPT4.NUM(/2)+1 2389 : SEGADJ,IPT4 2390 : IPT4.NUM(1,NBELEM)=NA 2391 : IPT4.NUM(2,NBELEM)=NB 2392 : BOOL1=.FALSE. 2393 : ELSE 2394 : BOOL1=.TRUE. 2395 : ENDIF 2396 : ENDDO 2397 : C Rangement des elements de IPT4 dans les indices de la table 2398 : C de sortie TAB11 2399 : C --> Ajout de(s) arete(s) dans TAB11.VISU et TAB11.ARTS 2400 : NBELEG0=IPTG.NUM(/2) 2401 : NBNN=IPTG.NUM(/1) 2402 : NBELEM=NBELEG0+IPT4.NUM(/2) 2403 : SEGADJ,IPTG 2404 : DO J=1,IPT4.NUM(/2) 2405 : JJ=J+NBELEG0 2406 : IPTG.NUM(1,JJ)=IPT4.NUM(1,J) 2407 : IPTG.NUM(2,JJ)=IPT4.NUM(2,J) 2408 : NBNN=2 2409 : NBELEM=1 2410 : SEGINI,IPT5 2411 : IPT5.ITYPEL=2 2412 : IPT5.NUM(1,1)=IPT4.NUM(1,J) 2413 : IPT5.NUM(2,1)=IPT4.NUM(2,J) 2413 : CALL ECCTAB(MTABA,'ENTIER ',JJ,0,' ',0,0, 2414 : & 'MAILLAGE',0,0,0,0,IPT5) 2415 : ENDDO 2416 : C --> Ajout de(s) arete(s) pour chaque cellule concernee dans 2417 : C TAB11.CELL.NRC1 2418 : DO J=1,LV1.LECT(/1) 2419 : NC1=LV1.LECT(J) 2420 : NRC1=IPT1.NUM(1,NC1) 2420 : CALL ECROBJ('POINT ',NRC1) 2420 : CALL ECROBJ('TABLE ',MTABC) 2420 : CALL EXIS 2420 : CALL LIRLOG(BOOL1,1,IRETOU) 2421 : IF (.NOT.BOOL1) THEN 2421 : CALL CRTABL(MTABI1) 2421 : CALL ECCTAB(MTABC,'POINT ',0,0,' ',0,NRC1, 2422 : & 'TABLE ',0,0,0,0,MTABI1) 2423 : NBNN=2 2424 : NBELEM=0 2425 : SEGINI,IPTI1 2426 : IPTI1.ITYPEL=2 2426 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 2427 : & 'MAILLAGE',0,0,0,0,IPTI1) 2428 : JG=0 2429 : SEGINI,LFI1 2429 : CALL ECCTAB(MTABI1,'MOT ',0,0,'FACS',0,0, 2430 : & 'LISTENTI',0,0,0,0,LFI1) 2431 : NBNN=1 2432 : NBELEM=0 2433 : SEGINI,IPTVI1 2434 : IPTVI1.ITYPEL=1 2434 : CALL ECCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, 2435 : & 'MAILLAGE',0,0,0,0,IPTVI1) 2436 : ENDIF 2436 : CALL ACCTAB(MTABC,'POINT ',0,0,' ',0,NRC1, 2437 : & 'TABLE ',0,0,0,0,MTABI1) 2438 : C ajout des elements de l'arete IPT4 dans l'indice 2439 : C TAB11.CELL.NRC1.VISU 2439 : CALL ACCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, 2440 : & 'MAILLAGE',0,0,0,0,IPTI1) 2441 : NBELEI10=IPTI1.NUM(/2) 2442 : NBNN=IPTI1.NUM(/1) 2443 : NBELEM=NBELEI10+IPT4.NUM(/2) 2444 : SEGADJ,IPTI1 2445 : DO K=1,IPT4.NUM(/2) 2446 : KK=K+NBELEI10 2447 : IPTI1.NUM(1,KK)=IPT4.NUM(1,K) 2448 : IPTI1.NUM(2,KK)=IPT4.NUM(2,K) 2449 : ENDDO 2449 : CALL ACCTAB(MTABI1,'MOT ',0,0,'FACS',0,0, 2450 : & 'LISTENTI',0,0,0,0,LFI1) 2450 : CALL ACCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, 2451 : & 'MAILLAGE',0,0,0,0,IPTVI1) 2452 : C renseignement des faces/voisins pour le polyedre NRC1 2453 : C boucle sur les autres cellules partageant l'arete IPT4 2454 : DO K=1,LV1.LECT(/1) 2455 : IF (K.EQ.J) GOTO 20 2456 : C la cellule NRC2 partage l'arete IPT4 2457 : NC2=LV1.LECT(K) 2458 : NRC2=IPT1.NUM(1,NC2) 2459 : C mais est elle voisine de NRC1 ? on regarde dans le 2460 : C maillage de ses voisins 2460 : CALL ACCTAB(MTAB1,'MOT ',0,0,'CELL',0,0, 2461 : & 'TABLE ',0,0,0,0,MTABCCC) 2461 : CALL ACCTAB(MTABCCC,'POINT ',0,0,' ',0,NRC1, 2462 : & 'TABLE ',0,0,0,0,MTABI00) 2462 : CALL ACCTAB(MTABI00,'MOT ',0,0,'VOIS',0,0, 2463 : & 'MAILLAGE',0,0,0,0,IPTVI00) 2464 : BOOL1=.FALSE. 2465 : DO L=1,IPTVI00.NUM(/2) 2466 : IF (IPTVI00.NUM(1,L).EQ.NRC2) THEN 2467 : BOOL1=.TRUE. 2468 : GOTO 26 2469 : ENDIF 2470 : ENDDO 2471 : 26 CONTINUE 2472 : IF (BOOL1) THEN 2473 : C si oui, NC1.NC2 forment une face mais laquelle ? 2474 : C on va checher le numero de la face dans ITFAC1 2475 : NF12=ITFAC1.LACT(NC1,NC2) 2476 : IF (NF12.EQ.0) THEN 2477 : C il s'agit d'une nouvelle face a creer 2478 : NF=NF+1 2479 : NF12=NF 2480 : ITFAC1.LACT(NC1,NC2)=NF12 2481 : ITFAC1.LACT(NC2,NC1)=NF12 2482 : ENDIF 2483 : C ajout du numero NF12 a la liste des faces si besoin 2483 : CALL ECRENT(NF12) 2483 : CALL ECROBJ('LISTENTI',LFI1) 2483 : CALL EXIS 2484 : SEGACT,LFI1 2484 : CALL LIRLOG(BOOL2,1,IRETOU) 2485 : IF (.NOT.BOOL2) THEN 2486 : JG=LFI1.LECT(/1)+1 2487 : SEGADJ,LFI1 2488 : LFI1.LECT(JG)=NF12 2489 : ENDIF 2490 : C ajout du point NRC2 au maillage des voisins 2491 : C si besoin 2491 : CALL ECROBJ('MAILLAGE',IPTVI1) 2491 : CALL ECROBJ('POINT ',NRC2) 2491 : CALL DANS 2492 : SEGACT,IPTVI1 2492 : CALL LIRLOG(BOOL2,1,IRETOU) 2493 : IF (.NOT.BOOL2) THEN 2494 : NBNN=IPTVI1.NUM(/1) 2495 : NBELEM=IPTVI1.NUM(/2)+1 2496 : SEGADJ,IPTVI1 2497 : IPTVI1.NUM(1,NBELEM)=NRC2 2498 : ENDIF 2499 : C ajout de(s) arete(s) dans la face NF12 si besoin 2499 : CALL ECRENT(NF12) 2499 : CALL ECROBJ('TABLE ',MTABF) 2499 : CALL EXIS 2499 : CALL LIRLOG(BOOL2,1,IRETOU) 2500 : IF (.NOT.BOOL2) THEN 2500 : CALL CRTABL(MTABF1) 2500 : CALL ECCTAB(MTABF,'ENTIER ',NF12,0,' ',0,0, 2501 : & 'TABLE ',0,0,0,0,MTABF1) 2502 : NBNN=2 2503 : NBELEM=0 2504 : SEGINI,IPTF1 2505 : IPTF1.ITYPEL=2 2505 : CALL ECCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 2506 : & 'MAILLAGE',0,0,0,0,IPTF1) 2507 : JG=0 2508 : SEGINI,LA1 2508 : CALL ECCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, 2509 : & 'LISTENTI',0,0,0,0,LA1) 2510 : 2511 : ENDIF 2511 : CALL ACCTAB(MTABF,'ENTIER ',NF12,0,' ',0,0, 2512 : & 'TABLE ',0,0,0,0,MTABF1) 2512 : CALL ACCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, 2513 : & 'MAILLAGE',0,0,0,0,IPTF1) 2513 : CALL ACCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, 2514 : & 'LISTENTI',0,0,0,0,LA1) 2515 : DO L=1,IPT4.NUM(/2) 2516 : NARETL=NBELEG0+L 2516 : CALL ECRENT(NARETL) 2516 : CALL ECROBJ('LISTENTI',LA1) 2516 : CALL EXIS 2517 : SEGACT,LA1 2517 : CALL LIRLOG(BOOL2,1,IRETOU) 2518 : IF (.NOT.BOOL2) THEN 2519 : JG=LA1.LECT(/1)+1 2520 : SEGADJ,LA1 2521 : LA1.LECT(JG)=NARETL 2522 : NBNN=IPTF1.NUM(/1) 2523 : NBELEM=IPTF1.NUM(/2)+1 2524 : SEGADJ,IPTF1 2525 : IPTF1.NUM(1,NBELEM)=IPT4.NUM(1,L) 2526 : IPTF1.NUM(2,NBELEM)=IPT4.NUM(2,L) 2527 : ENDIF 2528 : ENDDO 2529 : ENDIF 2530 : 20 CONTINUE 2531 : ENDDO 2532 : ENDDO 2533 : 11 CONTINUE 2534 : ENDDO 2535 : ENDIF 2536 : C 2537 : C 2538 : c enregistrer le champ par point 2539 : IF (MCHPOI .ne. 0) then 2539 : CALL ECCTAB(MTAB11,'MOT ',0,0,'POND',0,0, 2540 : & 'CHPOINT',0,0,0,0,MCHPOI) 2541 : ENDIF 2542 : C 2543 : C----------------------------------------------------------------------C 2544 : C PARTIE 4 C 2545 : C On rend la table MTAB11 C 2546 : C----------------------------------------------------------------------C 2547 : C 2548 : C 2549 : MTAB1=MTAB11 2550 : 999 RETURN 2551 : END 2552 : 2553 : 2554 :
© Cast3M 2003 - All rights reserved.
Disclaimer