C VORO1 SOURCE PV 20/04/01 21:16:44 10569 C SUBROUTINE VORO1(IPT1,IPT2,LNBOIT,ILEA1,ITRC1,DDMAX, & MTAB1) C----------------------------------------------------------------------C C Partition de Voronoi (non coupee) d'un nuage de points (IPT1) a C C partir de la triangulation de Delaunay (IPT2) C C C C IPT1 = pointeur sur le maillage de points, centres des C C polyedres (elements POI1) C C IPT2 = pointeur sur le maillage de la triangulation de C C Delaunay de IPT1 (avec les coins de la boite de C C triangulation) (elements TRI3 ou TET4) C C LNBOIT = tableau contenant les numeros des noeuds formant les C C coins de la boite de triangulation C C ILEA1 = pointeur sur le segment MADJACEN contenant la table des C C connectivite des triangles de IPT1 C C ITRC1 = pointeur sur le segment MCIRCONS contenant la table des C C coordonnees des cercles circonscrits des triangles de C C IPT1 C C XZERO = critere de distance pour considerer deux centres de C C cercles circonscrits confondus C C C C MTAB1 = table de maillage des polyedres non coupes C C----------------------------------------------------------------------C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) c -INC PPARAM -INC CCOPTIO -INC SMELEME -INC SMCOORD -INC CCGEOME -INC SMLENTI -INC SMTABLE DIMENSION LNBOIT(8) DIMENSION XPI1(3) LOGICAL BOOL0,BOOL1,BOOL2,BOOL3 C POINTEUR LCCC1.MLENTI,LAR1.MLENTI,LAR2.MLENTI C Tables, maillages et listes pour la table de sortie POINTEUR IPTG.MELEME,MTABC.MTABLE,MTABI1.MTABLE POINTEUR IPTI1.MELEME,LFI1.MLENTI,IPTVI1.MELEME,MTABF.MTABLE POINTEUR IPTF1.MELEME,LA1.MLENTI,MTABA.MTABLE,IPTA1.MELEME C SEGMENT,MADJACEN INTEGER LEAC(NBL1,IDIM+1,2) ENDSEGMENT POINTEUR ILEA1.MADJACEN C SEGMENT,MCIRCONS REAL*8 TRC(NBE1,4) ENDSEGMENT POINTEUR ITRC1.MCIRCONS C SEGMENT,MVORO INTEGER LACT(N1,N2) ENDSEGMENT POINTEUR ITL1.MVORO,ITCEL1.MVORO,ITCEL2.MVORO,ITARE1.MVORO POINTEUR ITFAC1.MVORO,ITFAC2.MVORO C XZERO = DDMAX*1.D-10 segact mcoord*mod C C 1) Pre-traitement : creation des points representant les centres C cercles circonscrits de chaque triangle (le meme numero est C attribue aux centres de cercles circonscrits superposes) C C Ajustement initial de MCOORD au nombre de cercles circonscrits NBTRI=ITRC1.TRC(/1) NBPTS0=nbpts NBPTS=NBPTS0+NBTRI SEGADJ,MCOORD C LCCC1.LECT(i)=j le triangle i a pour Centre de Cerclce Circonscrit C le point numero j JG=NBTRI SEGINI,LCCC1 I1=0 C Boucle sur les centres des cercles circonscrits des triangles DO I=1,NBTRI DO K=1,IDIM XPI1(K)=ITRC1.TRC(I,K+1) ENDDO DO J=1,NBTRI NCJ=0 IF (J.GE.I) GOTO 101 C Test si les coordonnees de I et J sont egales DIJ1=0. DO K=1,IDIM XPJ=ITRC1.TRC(J,K+1) DIJ1=DIJ1+(XPI1(K)-XPJ)*(XPI1(K)-XPJ) c XDEN=XPI1(K) c IF (XDEN.EQ.0.) XDEN=1. c IF (ABS((XPI1(K)-XPJ)/XDEN).LT.XZERO) NCJ=NCJ+1 ENDDO IF(DIJ1.LT.XZERO) NCJ=IDIM C Si les centres I et J sont superposes IF (NCJ.EQ.IDIM) THEN LCCC1.LECT(I)=LCCC1.LECT(J) GOTO 102 ENDIF 101 CONTINUE ENDDO C Creation d'un nouveau point dans MCOORD I1=I1+1 NREFI=NBPTS0+I1 LCCC1.LECT(I)=NREFI C Ecriture des ses coordonnees DO J=1,IDIM XCOOR(((NREFI-1)*(IDIM+1))+J)=XPI1(J) ENDDO C et de la densite associee, choisie nulle XCOOR(((NREFI-1)*(IDIM+1))+(IDIM+1))=0. 102 CONTINUE ENDDO C Ajustement final de MCOORD (on a ajoute I1 points) IF (I1.NE.NBTRI) THEN NBPTS=NBPTS0+I1 SEGADJ,MCOORD ENDIF C C C 2.0) Initialisations avant creation du maillage des polyedres C IPTG : Maillage des aretes des polyedres (elements SEG2) NBSOUS=0 NBREF=0 NBNN=2 NBELEM=I1*4 SEGINI,IPTG IPTG.ITYPEL=2 NBARETES=0 C Table MTAB1 de sortie contenant les maillages des polyedres CALL CRTABL(MTAB1) CALL ECCTAB(MTAB1,'MOT ',0,0,'VISU',0,0, & 'MAILLAGE',0,0,0,0,IPTG) CALL CRTABL(MTABC) CALL ECCTAB(MTAB1,'MOT ',0,0,'CELL',0,0, & 'TABLE ',0,0,0,0,MTABC) IF (IDIM.EQ.3) THEN CALL CRTABL(MTABF) CALL ECCTAB(MTAB1,'MOT ',0,0,'FACS',0,0, & 'TABLE ',0,0,0,0,MTABF) ENDIF CALL CRTABL(MTABA) CALL ECCTAB(MTAB1,'MOT ',0,0,'ARTS',0,0, & 'TABLE ',0,0,0,0,MTABA) C ITL1 : Table indiquant si la ligne entre deux centres de cercle C circonscrits existe deja ou non N1=I1 N2=I1 SEGINI,ITL1 C ITCEL1 : Table des polyedres voisins d'un polyedre C ITCEL2 : Table des aretes d'un polyedre N1=10 N2=IPT1.NUM(/2) SEGINI,ITCEL1,ITCEL2 C ITARE1 : Table des cellules appuyees sur une arete, on la stocke C dans la table de sortie pour pouvoir la retrouver lors du C decoupage dans VORO2 C ITARE1.LACT(ICE,NLA)=I C NLA : numero (dans IPTG) de l'arete C ICE : index des cellules appuyees sur NLA (de 1 a n, puis 0) C I : numero (dans IPT1) de la ICE-ieme cellule appuyee sur NLA N1=10 N2=IPTG.NUM(/2) SEGINI,ITARE1 CALL ECCTAB(MTAB1,'MOT ',0,0,'TAR',0,0, & 'MVORO ',0,0,0,0,ITARE1) C C C C 2.1) Maillage des cellules de Voronoi en dimension 2 C Boucle sur les centres des polyedres de Voronoi IF (IDIM.EQ.2) THEN DO I=1,IPT1.NUM(/2) NREF_I=IPT1.NUM(1,I) C Initialisation des indices pour la cellule NREF_I CALL CRTABL(MTABI1) CALL ECCTAB(MTABC,'POINT ',0,0,' ',0,NREF_I, & 'TABLE ',0,0,0,0,MTABI1) C --> maillage de la cellule NBNN=2 NBELEM=I1*4 SEGINI,IPTI1 IPTI1.ITYPEL=2 NBELEI1=0 CALL ECCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, & 'MAILLAGE',0,0,0,0,IPTI1) C --> liste des aretes JG=I1*4 SEGINI,LA1 CALL ECCTAB(MTABI1,'MOT ',0,0,'ARTS',0,0, & 'LISTENTI',0,0,0,0,LA1) C --> maillage des voisns NBNN=1 NBELEM=I1*4 SEGINI,IPTVI1 IPTVI1.ITYPEL=1 NVI1=0 CALL ECCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, & 'MAILLAGE',0,0,0,0,IPTVI1) C Recherche du premier triangle appuye sur NREFA NREFA=IPT1.NUM(1,I) DO J=1,IPT2.NUM(/2) DO K=1,IPT2.NUM(/1) IF (IPT2.NUM(K,J).EQ.NREFA) GOTO 211 ENDDO ENDDO 211 NT0=J NA=K NB=NA+1 IF (NB.EQ.4) NB=1 NC=6-NA-NB NREFB=IPT2.NUM(NB,NT0) NREFC=IPT2.NUM(NC,NT0) C Balayage des triangles appuyes sur NREFA en tournant autour C de ce point, a patrir du triangle NT0 C La boucle s'arrete quand on est revenu au triangle NT0 NT1=NT0 NRCC1=LCCC1.LECT(NT1) NCVI=0 NAI=0 DO J=1,I1*4 IJOK=0 C Determination du point NRCC2 a relier a NRCC1 NT2=ILEA1.LEAC(NT1,NB,1) NRCC2=LCCC1.LECT(NT2) C Si la liaison entre CC1 et CC2 existe, on la recupere et C on passe au triangle suivant IF (ITL1.LACT(NRCC1-NBPTS0,NRCC2-NBPTS0).EQ.1) THEN C On determine le numero de l'arete en question DO K=1,ITCEL2.LACT(/1) NVI=0 DO L=1,IPT1.NUM(/2) IF (IPT1.NUM(1,L).EQ.NREFC) NVI=L ENDDO IF (NVI.EQ.0) THEN DO L=1,4 IF (LNBOIT(L).EQ.NREFC) NVI=-1*L ENDDO ENDIF NAR1=ITCEL2.LACT(K,NVI) NPA=IPTG.NUM(1,NAR1) NPB=IPTG.NUM(2,NAR1) IF (((NPA.EQ.NRCC1).AND.(NPB.EQ.NRCC2)).OR. & ((NPA.EQ.NRCC2).AND.(NPB.EQ.NRCC1))) THEN GOTO 212 ENDIF ENDDO ENDIF C Si les centres CC1 et CC2 sont confondus, il n'y a pas de C liaison a creer, on passe au triangle suivant IF (NRCC1.EQ.NRCC2) THEN GOTO 214 ENDIF C Dans les autres cas, il y a bien creation d'une nouvelle C arete dans le maillage IPTG de la cellule I IJOK=1 NBARETES=NBARETES+1 IPTG.NUM(1,NBARETES)=NRCC1 IPTG.NUM(2,NBARETES)=NRCC2 NAR1=NBARETES C Remplissage de la table ILT1 des liaisons entre centres C de cercles circonscrits ITL1.LACT(NRCC1-NBPTS0,NRCC2-NBPTS0)=1 ITL1.LACT(NRCC2-NBPTS0,NRCC1-NBPTS0)=1 C Remplissage de la table ITCEL1 des numeros des cellules C voisines de la cellule I NVI=0 DO K=1,IPT1.NUM(/2) IF (IPT1.NUM(1,K).EQ.NREFC) NVI=K ENDDO IF (NVI.EQ.0) THEN DO K=1,4 IF (LNBOIT(K).EQ.NREFC) NVI=-1*K ENDDO ENDIF 212 IF (NVI.GT.0) THEN NCVI=NCVI+1 IF (NCVI.GT.ITCEL1.LACT(/1)) THEN N1=ITCEL1.LACT(/1)+10 N2=ITCEL1.LACT(/2) SEGADJ,ITCEL1 ENDIF ITCEL1.LACT(NCVI,I)=NVI ENDIF C Remplissage de la table ITCEL2 des numeros d'aretes C appartenant a la cellule I NAI=NAI+1 IF (NAI.GT.ITCEL2.LACT(/1)) THEN N1=ITCEL2.LACT(/1)+10 N2=ITCEL2.LACT(/2) SEGADJ,ITCEL2 ENDIF ITCEL2.LACT(NAI,I)=NAR1 C Remplissage de la table ITARE1 des numeros de cellules C contenant l'arrete creee/recuperee ICE=0 DO K=1,ITARE1.LACT(/1) IF (ITARE1.LACT(K,NAR1).EQ.0) THEN ICE=K GOTO 213 ENDIF ENDDO IF (ICE.EQ.0) THEN ICE=ITARE1.LACT(/1)+1 N1=ITARE1.LACT(/1)+10 N2=ITARE1.LACT(/2) SEGADJ,ITARE1 ENDIF 213 ITARE1.LACT(ICE,NAR1)=I C Ajout d'un element au maillage IPTI1 de la cellule I NBELEI1=NBELEI1+1 IPTI1.NUM(1,NBELEI1)=NRCC1 IPTI1.NUM(2,NBELEI1)=NRCC2 C Ajout du numero d'arete dans la liste de la cellule I LA1.LECT(NBELEI1)=NAR1 C Ajout d'un voisin a la cellule I IF (NVI.GT.0) THEN NVI1=NVI1+1 NREFVI=IPT1.NUM(1,NVI) IPTVI1.NUM(1,NVI1)=NREFVI ENDIF C Ajout de l'arete dans MTABA IF (IJOK.EQ.1) THEN NBNN=2 NBELEM=1 SEGINI,IPTA1 IPTA1.ITYPEL=2 IPTA1.NUM(1,1)=NRCC1 IPTA1.NUM(2,1)=NRCC2 CALL ECCTAB(MTABA,'ENTIER ',NBARETES,0,' ',0,0, & 'MAILLAGE',0,0,0,0,IPTA1) ENDIF C Test d'arret de la boucle (est on revenu sur le triangle C de depart ?) 214 IF (NT2.EQ.NT0) THEN GOTO 215 ENDIF C Nouveau points A, B et C pour le prochain triangle DO K=1,IPT2.NUM(/1) IF (IPT2.NUM(K,NT2).EQ.NREFA) NA=K IF (IPT2.NUM(K,NT2).EQ.NREFC) NB=K ENDDO NC=6-NA-NB NREFB=IPT2.NUM(NB,NT2) NREFC=IPT2.NUM(NC,NT2) NT1=NT2 NRCC1=LCCC1.LECT(NT1) ENDDO 215 CONTINUE C Ajustement des maillages de la cellule I, de ses voisins et C de la liste des aretes NBNN=IPTI1.NUM(/1) NBELEM=NBELEI1 SEGADJ,IPTI1 JG=NBELEI1 SEGADJ,LA1 NBNN=IPTVI1.NUM(/1) NBELEM=NVI1 SEGADJ,IPTVI1 ENDDO C C C C 2.2) Maillage des cellules de Voronoi en dimension 3 ELSEIF (IDIM.EQ.3) THEN C ITFAC1 : tableau indiquant si la face (I,J) existe deja C ou I et J sont des numeros de cellule N1=IPT1.NUM(/2) N2=N1 SEGINI,ITFAC1 C ITFAC2 : tableau indiquant si la face (I,J) existe deja C ou I est un numero de cellule exterieures infinie et J un C numero des coins de la boite de triangulation N1=IPT1.NUM(/2) N2=8 SEGINI,ITFAC2 C Boucle sur les centres des polyedres de Voronoi NBFACES=0 DO I=1,IPT1.NUM(/2) C On va chercher la table de la cellule I NREF_I=IPT1.NUM(1,I) C Initialisation des indices pour la cellule NREF_I CALL CRTABL(MTABI1) CALL ECCTAB(MTABC,'POINT ',0,0,' ',0,NREF_I, & 'TABLE ',0,0,0,0,MTABI1) C --> maillage de la cellule NBNN=2 NBELEM=I1*4 SEGINI,IPTI1 IPTI1.ITYPEL=2 NBELEI1=0 CALL ECCTAB(MTABI1,'MOT ',0,0,'VISU',0,0, & 'MAILLAGE',0,0,0,0,IPTI1) C --> liste des faces JG=I1*4 SEGINI,LFI1 NFI1=0 CALL ECCTAB(MTABI1,'MOT ',0,0,'FACS',0,0, & 'LISTENTI',0,0,0,0,LFI1) C --> maillage des voisns NBNN=1 NBELEM=I1*4 SEGINI,IPTVI1 IPTVI1.ITYPEL=1 IV1=0 CALL ECCTAB(MTABI1,'MOT ',0,0,'VOIS',0,0, & 'MAILLAGE',0,0,0,0,IPTVI1) C Recherche du premier tetraedre appuye sur NREFA NCVI=0 NREFA=IPT1.NUM(1,I) DO J=1,IPT2.NUM(/2) DO K=1,IPT2.NUM(/1) IF (IPT2.NUM(K,J).EQ.NREFA) GOTO 221 ENDDO ENDDO 221 NT0=J NA=K NB=NA+1 NC=NA-1 IF (NB.EQ.5) NB=1 IF (NC.EQ.0) NC=3 ND=10-NA-NB-NC C Liste des aretes appuyees sur NREFA autour desquelles on C va tourner. On incremente ces liste au fur et a mesure que C l'on rencontre de nouveau tetraedres appuyes sur NT0 JG=IPTG.NUM(/2) SEGINI,LAR1,LAR2 C On connait les trois premieres aretes appuyees sur NREFA LAR1.LECT(1)=IPT2.NUM(NB,NT0) LAR1.LECT(2)=IPT2.NUM(NC,NT0) LAR1.LECT(3)=IPT2.NUM(ND,NT0) C et les numeros des premiers tetraedres contenant ces aretes LAR2.LECT(1)=NT0 LAR2.LECT(2)=NT0 LAR2.LECT(3)=NT0 C Balayage des aretes des tetraedres appuyes sur NREFA, C a patrir du tetraedre NT0 (on parcoure LAR1) C Chaque nouvelle arete definit potentielement une nouvelle C face de la partition de Voronoi C La liste LAR1 s'incremente, la boucle s'arrete lorsque l'on C a parcouru toutes les aretes dans LAR1 DO J=1,LAR1.LECT(/1) IJOK=0 C On tourne autour de l'arete J de la liste LAR1 NREFB=LAR1.LECT(J) C Test d'arret de la boucle, si on a fini de parcourir la C liste LAR1, on a fini ! IF (NREFB.EQ.0) THEN GOTO 229 ENDIF C Calcul de NVI, numero de cette cellule voisine NVI=0 DO K=1,IPT1.NUM(/2) IF (IPT1.NUM(1,K).EQ.NREFB) NVI=K ENDDO C ou bien s'il s'agit d'un sommet de la boite de C triangulation IF (NVI.EQ.0) THEN DO K=1,8 IF (LNBOIT(K).EQ.NREFB) NVI=-1*K ENDDO ENDIF C NFIJ est le numero de la face entre I et J C (si elle existe deja) IF (NVI.GT.0) THEN NFIJ=ITFAC1.LACT(NVI,I) ELSE NVII=-1*NVI NFIJ=ITFAC2.LACT(I,NVII) ENDIF C IJOK est mis a 1 quand on identifie une face non deja C existante (mais elle peut etre vide !) IF (NFIJ.EQ.0) IJOK=1 C Dans le cas d'une potentielle nouvelle face, on prepare C son maillage et sa liste des aretes IF (IJOK.EQ.1) THEN NBNN=2 NBELEM=LAR1.LECT(/1) SEGINI,IPTF1 IPTF1.ITYPEL=2 JG=LAR1.LECT(/1) SEGINI,LA1 ENDIF NBELEF1=0 C Le premier tetraedre balaye est NT0 NT0=LAR2.LECT(J) C Calcul de NA,NB,NC et ND dans NT0 DO K=1,IPT2.NUM(/1) IF (IPT2.NUM(K,NT0).EQ.NREFA) NA=K IF (IPT2.NUM(K,NT0).EQ.NREFB) NB=K ENDDO IF ((NA.EQ.1).AND.(NB.EQ.2)) NC=3 IF ((NA.EQ.1).AND.(NB.EQ.3)) NC=2 IF ((NA.EQ.1).AND.(NB.EQ.4)) NC=2 IF ((NA.EQ.2).AND.(NB.EQ.1)) NC=3 IF ((NA.EQ.2).AND.(NB.EQ.3)) NC=1 IF ((NA.EQ.2).AND.(NB.EQ.4)) NC=1 IF ((NA.EQ.3).AND.(NB.EQ.1)) NC=2 IF ((NA.EQ.3).AND.(NB.EQ.2)) NC=1 IF ((NA.EQ.3).AND.(NB.EQ.4)) NC=1 IF ((NA.EQ.4).AND.(NB.EQ.1)) NC=2 IF ((NA.EQ.4).AND.(NB.EQ.2)) NC=1 IF ((NA.EQ.4).AND.(NB.EQ.3)) NC=1 ND=10-NA-NB-NC NREFC=IPT2.NUM(NC,NT0) NREFD=IPT2.NUM(ND,NT0) C Boucle autour de l'arete NREFA-NREFB, en partant du C tetraedre NT0, on construit ainsi la face (I,J) en C tournant autour de ses aretes C La boucle s'arrete quand on est revenu au tetraedre NT0 NT1=NT0 NRCC1=LCCC1.LECT(NT1) DO K=1,I1*4 NT2=ILEA1.LEAC(NT1,NC,1) C On ajoute NREFC et NREFD a la liste des aretes a C balayer plus tard (si non deja presents) DO L=1,LAR1.LECT(/1) IF (LAR1.LECT(L).EQ.NREFC) GOTO 222 IF (LAR1.LECT(L).EQ.0) THEN LAR1.LECT(L)=NREFC LAR2.LECT(L)=NT2 GOTO 222 ENDIF ENDDO 222 CONTINUE DO L=1,LAR1.LECT(/1) IF (LAR1.LECT(L).EQ.NREFD) GOTO 223 IF (LAR1.LECT(L).EQ.0) THEN LAR1.LECT(L)=NREFD LAR2.LECT(L)=NT2 GOTO 223 ENDIF ENDDO 223 CONTINUE NRCC2=LCCC1.LECT(NT2) C C Plusieurs cas sont possibles pour l'arete reliant C CC1 et CC2 : C C CAS 1 : la liaison entre CC1 et CC2 existe C => on la recupere et on passe au tetraedre suivant IF (ITL1.LACT(NRCC1-NBPTS0,NRCC2-NBPTS0).EQ.1) THEN C On cherche le numero de l'arete existante reliant C deja CC1 et CC2 DO L=1,IPTG.NUM(/2) NPA=IPTG.NUM(1,L) NPB=IPTG.NUM(2,L) IF (((NPA.EQ.NRCC1).AND.(NPB.EQ.NRCC2)).OR. & ((NPA.EQ.NRCC2).AND.(NPB.EQ.NRCC1))) THEN NLA=L GOTO 224 ENDIF ENDDO ENDIF C CAS 2 : les centres CC1 et CC2 sont confondus C => il n'y a pas d'arete a faire ! on passe au C tetraedre suivant IF (NRCC1.EQ.NRCC2) THEN GOTO 227 ENDIF C CAS 3 : creation d'une nouvelle arete dans IPTG NBARETES=NBARETES+1 IPTG.NUM(1,NBARETES)=NRCC1 IPTG.NUM(2,NBARETES)=NRCC2 NLA=NBARETES C et ajout de cette arete dans MTABA NBNN=2 NBELEM=1 SEGINI,IPTA1 IPTA1.ITYPEL=2 IPTA1.NUM(1,1)=NRCC1 IPTA1.NUM(2,1)=NRCC2 CALL ECCTAB(MTABA,'ENTIER ',NLA,0,' ',0,0, & 'MAILLAGE',0,0,0,0,IPTA1) C Remplissage de la table ILT1 des liaisons en centres C de cercles circonscrits ITL1.LACT(NRCC1-NBPTS0,NRCC2-NBPTS0)=1 ITL1.LACT(NRCC2-NBPTS0,NRCC1-NBPTS0)=1 C 224 CONTINUE C On incremente le nombre d'arete dans cette face NBELEF1=NBELEF1+1 C Dans le cas d'une potentielle nouvelle face, on C incremente le maillage de la face et la liste de ses C aretes (si l'arete en question n'est pas presente) IF (IJOK.EQ.1) THEN DO L=1,LA1.LECT(/1) IF (LA1.LECT(L).EQ.NLA) GOTO 230 IF (LA1.LECT(L).EQ.0) GOTO 231 ENDDO 231 IF (NBELEF1.GT.LA1.LECT(/1)) THEN JG=LA1.LECT(/1)+10 SEGADJ,LA1 NBNN=IPTF1.NUM(/1) NBELEM=IPTF1.NUM(/2)+10 SEGADJ,IPTF1 ENDIF C ajout dans la liste des aretes pour cette face LA1.LECT(NBELEF1)=NLA C et au maillage de la face IPTF1.NUM(1,NBELEF1)=NRCC1 IPTF1.NUM(2,NBELEF1)=NRCC2 230 CONTINUE ENDIF C Dans le cas d'une potentielle nouvelle face, si le C nombre d'aretes est egal a 3, il s'agit bel et bien C d'une face IF ((NBELEF1.EQ.3).AND.(IJOK.EQ.1)) THEN C on l'ajoute a MTABF NBFACES=NBFACES+1 NFIJ=NBFACES CALL CRTABL(MTABF1) CALL ECCTAB(MTABF,'ENTIER ',NFIJ,0,' ',0,0, & 'TABLE ',0,0,0,0,MTABF1) CALL ECCTAB(MTABF1,'MOT ',0,0,'VISU',0,0, & 'MAILLAGE',0,0,0,0,IPTF1) CALL ECCTAB(MTABF1,'MOT ',0,0,'ARTS',0,0, & 'LISTENTI',0,0,0,0,LA1) ENDIF C Si le nombre d'arete est egal a 3, le point J est bel C et bien un voisin de I IF (NBELEF1.EQ.3) THEN C S'il s'agit d'une face exterieure (J est sommet de C la boite de triangulation) on l'indique juste dans C ITFAC2 IF (NVI.LT.0) THEN NVII=-1*NVI ITFAC2.LACT(I,NVII)=NFIJ C Sinon, on ajoute un voisin et on marque la face C dans ITFAC1 ELSE NREFVI=IPT1.NUM(1,NVI) IV1=IV1+1 IF (IV1.GT.IPTVI1.NUM(/2)) THEN NBNN=IPTVI1.NUM(/1) NBELEM=IPTVI1.NUM(/2)+10 SEGADJ,IPTVI1 ENDIF IPTVI1.NUM(1,IV1)=NREFVI ITFAC1.LACT(I,NVI)=NFIJ ITFAC1.LACT(NVI,I)=NFIJ ENDIF C Ajout du numero global de la face dans la liste des C faces de la cellule NREF_I NFI1=NFI1+1 IF (NFI1.GT.LFI1.LECT(/1)) THEN JG=LFI1.LECT(/1)+10 SEGADJ,LFI1 ENDIF LFI1.LECT(NFI1)=NFIJ ENDIF C Remplissage de la table ITCEL2 des numeros d'aretes C appartenant a la cellule I DO L=1,ITCEL2.LACT(/1) IF (ITCEL2.LACT(L,I).EQ.NLA) THEN GOTO 227 ENDIF IF (ITCEL2.LACT(L,I).EQ.0) GOTO 225 ENDDO 225 NBELEI1=NBELEI1+1 IF (NBELEI1.GT.ITCEL2.LACT(/1)) THEN N1=ITCEL2.LACT(/1)+10 N2=ITCEL2.LACT(/2) SEGADJ,ITCEL2 ENDIF ITCEL2.LACT(NBELEI1,I)=NLA C Ajout de l'arete au maillage de la cellule NREF_I IPTI1.NUM(1,NBELEI1)=NRCC1 IPTI1.NUM(2,NBELEI1)=NRCC2 C Remplissage de la table ITARE1 des numeros de cellules C contenant l'arrete creee/recuperee ICE=0 DO L=1,ITARE1.LACT(/1) IF (ITARE1.LACT(L,NLA).EQ.0) THEN ICE=L GOTO 226 ENDIF ENDDO IF (ICE.EQ.0) THEN ICE=ITARE1.LACT(/1)+1 N1=ITARE1.LACT(/1)+10 N2=ITARE1.LACT(/2) SEGADJ,ITARE1 ENDIF 226 ITARE1.LACT(ICE,NLA)=I C Test d'arret de la boucle (est on revenu sur le C tetraedre de depart ?) 227 IF (NT2.EQ.NT0) THEN C Remplissage de la table ITCEL1 des numeros des C cellules voisines de la cellule I IF (NBELEF1.GE.3) THEN NCVI=NCVI+1 IF (NCVI.GT.ITCEL1.LACT(/1)) THEN N1=ITCEL1.LACT(/1)+10 N2=ITCEL1.LACT(/2) SEGADJ,ITCEL1 ENDIF ITCEL1.LACT(NCVI,I)=NVI IF (IJOK.EQ.1) THEN C ajustement du maillage de la face (I,J) NBNN=IPTF1.NUM(/1) NBELEM=NBELEF1 SEGADJ,IPTF1 C et de la liste de ses aretes JG=NBELEF1 SEGADJ,LA1 ENDIF ENDIF GOTO 228 ENDIF C Nouveau points A, B, C et D pour le prochain tetraedre DO L=1,IPT2.NUM(/1) IF (IPT2.NUM(L,NT2).EQ.NREFA) NA=L IF (IPT2.NUM(L,NT2).EQ.NREFB) NB=L IF (IPT2.NUM(L,NT2).EQ.NREFD) NC=L ENDDO ND=10-NA-NB-NC NREFC=IPT2.NUM(NC,NT2) NREFD=IPT2.NUM(ND,NT2) NT1=NT2 NRCC1=LCCC1.LECT(NT1) ENDDO 228 CONTINUE ENDDO 229 CONTINUE C Ajustement du maillage de la cellule I, de la liste de ses C faces et de ses voisins NBNN=IPTI1.NUM(/1) NBELEM=NBELEI1 SEGADJ,IPTI1 JG=NFI1 SEGADJ,LFI1 NBNN=IPTVI1.NUM(/1) NBELEM=IV1 SEGADJ,IPTVI1 ENDDO SEGSUP,LAR1,LAR2,ITFAC1,ITFAC2 ENDIF C C 2.3) Ajustement final de IPTG et ITARE1 NBNN=IPTG.NUM(/1) NBELEM=NBARETES SEGADJ,IPTG N1=ITARE1.LACT(/1) N2=NBARETES SEGADJ,ITARE1 C C C 3) Menage avant de quitter SEGSUP,LCCC1,ITL1,ITCEL1,ITCEL2,ILEA1,ITRC1 999 RETURN END