C TIRET3 SOURCE PV 22/04/19 16:18:14 11344 C NOUVELLE VERSION DE TIRET C LA BOUCLE SE FAIT SUR LES ELEMENTS ET NON PLUS SUR LES NOEUDS C CETTE VERSION EST ADAPTEE A GIBI C SUBROUTINE TIRET3(XPROJ,MELEME,ICPR,XMIN,XMAX,YMIN,YMAX,IVU,NELEM, # TMIN,TMAX,MCOUP) C NUMNP NOMBRE DE POINTS A TRACER C NBELEM NOMBRE D'ELEMENTS DANS L'OBJET ELEMENTAIRE C MCOUP <>0 => COUPE LA PLAN DE COUPE EST DANS LA DERNIERE C COMPOSANTE DE MELEME ET SES NOEUDS SONT VUS IMPLICIT INTEGER(I-N) C C TABLEAUX : KTR TRIANGLES PAR ELEMENTS C KTL NOMBRE DE TRIANGLES PAR ELEM C NPZON ZONE A LAQUELLE APPARTIENT UN POINT C INDZON POINTEUR SUR NPOIN C NPOIN TABLEAU DES POINTS CLASSES C IVU TABLEAU VU PAS VU (POINT) C C----------------------------------------------------------------------- C C SEGMENTS C SEGMENT/XPROJ/(X(3,ITE)) SEGMENT IVU(NUMNP) SEGMENT ICPR(0) SEGMENT/TRAV/(NTAUX(NUMNP),INDZON(NUMNP+2),NPZON(NUMNP+2), # IJK(NUMNP+2),NPOIN(NUMNP+2)) SEGMENT MCOUP(0) REAL XMIN,XMAX,YMIN,YMAX C C----------------------------------------------------------------------- -INC PPARAM -INC CCOPTIO -INC CCGEOME -INC CCREEL *- -INC SMELEME NMULT=3 SEGACT ICPR NUMNP=X(/2) NUPLUS=NUMNP+1 * CALL LIRENT(NMULT,0,IRETOU) GENANT EN FULLSCREEN * IF (IRETOU.EQ.1) NMULT=MAX(1,NMULT) SEGINI TRAV C INITIALISATION DE L'HORLOGE XDIST=XMAX-XMIN YDIST=YMAX-YMIN ** write(6,*) 'tiret3 xdist ydist ',xdist,ydist cbp : petites limitations pour les cas XDIST et/ou YDIST tres petits if(XDIST.le.XSPETI.and.YDIST.le.XSPETI) then XDIST=1.D0 YDIST=1.D0 elseif(XDIST.lt.(1.D-5*YDIST)) then XDIST=1.D-5 * YDIST elseif(YDIST.lt.(1.D-5*XDIST)) then YDIST=1.D-5 * XDIST endif SMO=1E-5*XDIST*YDIST/NUMNP SMA=1E-5*(XDIST*YDIST/NUMNP)**1.5 C CALCUL DU DECOUPAGE EN ZONE NBZONE=NUMNP*NMULT RAP=YDIST/XDIST NDEC=MAX(INT(SQRT(REAL(NBZONE)/RAP)),1) MDEC=NBZONE/NDEC NBZONE=NDEC*MDEC IF (IIMPI.NE. 0) WRITE (IOIMP,8933) NDEC,MDEC,NBZONE 8933 FORMAT(' DECOUPAGE EN X ',I6,' EN Y ',I6,' SOIT ',I8,'ZONES ') XDEC=XDIST/NDEC*1.00001 YDEC=YDIST/MDEC*1.00001 IF (IIMPI.NE. 0) WRITE (IOIMP,8935) XMIN,XMAX,YMIN,YMAX,XDEC,YDEC 8935 FORMAT (' XMIN,XMAX,YMIN,YMAX,XDEC,YDEC ',6G12.5) C CLASSEMENT DES POINTS DO 40 I=1,NUPLUS+1 INDZON(I)=0 IJK(I)=0 40 CONTINUE DO 2 I=1,NUMNP IZONE=INT(real((X(1,I)-XMIN)/XDEC))+1+ & INT(real((X(2,I)-YMIN)/YDEC))*NDEC if(izone.lt.1) then ** write(6,*) 'tiret3 izone xmin xdec x',xmin,xdec,x(1,i) endif NTAUX(I)=IZONE IPOINT=IZONE/NMULT+1 INDZON(IPOINT)=INDZON(IPOINT)+1 2 CONTINUE C INDZON= NB DE POINTS PAR SUPER-ZONE DE CLASSEMENT DO 4 I=2,NUPLUS INDZON(I)=INDZON(I-1)+INDZON(I) 4 CONTINUE C INDZON(I)=FIN DE SUPER-ZONE I C RANGEMENT DES POINTS PAR SUPER-ZONE NBPOIN=INDZON(NUPLUS) INDZON(NUPLUS+1)=NBPOIN DO 5 I=1,NUMNP IZONE=NTAUX(I) IPOINT=IZONE/NMULT+1 IAD=INDZON(IPOINT) NPZON(IAD)=I INDZON(IPOINT)=IAD-1 5 CONTINUE C TRI DANS CHAQUE SUPER-ZONE C INDZON(I+1) EST DEVENU LA FIN DE SUPER-ZONE I DO 6 I=1,NUPLUS IDEB=INDZON(I)+1 IFIN=INDZON(I+1) IF (IDEB.GT.IFIN) GOTO 6 C IL FAUT ORDONNER NPZON SUIVANT NTAUX ET METTRE LE RESULTAT DANS NPOIN IJK(1)=IDEB IJK(2)=IFIN IZ=2 8 IZ=IZ-1 IF (IZ.LE.0) GOTO 9 IPB=IJK(IZ*2-1) IPH=IJK(IZ*2) IF (IPB.GE.IPH) GOTO 8 JPB=IPB-1 JPH=IPH+1 C CALCUL DU PIVOT IPV=0 DO 7 J=IPB,IPH IPV=IPV+NTAUX(NPZON(J)) 7 CONTINUE IPV=IPV/(IPH-IPB+1) 42 JPB=JPB+1 IF (JPH.EQ.JPB) GOTO 45 IF (NTAUX(NPZON(JPB)).GE.IPV) GOTO 43 GOTO 42 43 JPH=JPH-1 IF (JPH.EQ.JPB) GOTO 45 IF (NTAUX(NPZON(JPH)).LE.IPV) GOTO 44 GOTO 43 44 IAUX=NPZON(JPB) C CORRECTION JPH AU LIEU DE IPH NPZON(JPB)=NPZON(JPH) NPZON(JPH)=IAUX GOTO 42 C JPH=JPB MAIS ATTENTION PEUT ETRE A L'EXTERIEUR 45 IF (JPB.GE.IPB) GOTO 47 JPB=JPB+1 JPH=JPH+2 GOTO 48 47 IF (JPH.LE.IPH) GOTO 49 JPB=JPB-2 JPH=JPH-1 GOTO 48 49 IF (NTAUX(NPZON(JPB)).GE.IPV) GOTO 50 IF (JPH.EQ.IPH) GOTO 51 52 JPH=JPH+1 GOTO 48 50 IF (JPB.EQ.IPB) GOTO 52 51 JPB=JPB-1 48 IF (JPB.EQ.IPB) GOTO 53 IJK(2*IZ)=JPB IZ=IZ+1 53 IF (JPH.EQ.IPH) GOTO 8 IJK(2*IZ)=IPH IJK(2*IZ-1)=JPH IZ=IZ+1 GOTO 8 9 CONTINUE DO 46 J=IDEB,IFIN NPOIN(J)=NPZON(J) 46 CONTINUE 6 CONTINUE C REMPLISSAGE DE NPZON DO 41 I=1,NUMNP NPZON(I)=NTAUX(NPOIN(I)) 41 CONTINUE C# TMIN=XGRAND TMIN=XSGRAN TMAX=0. DO 30 I=1,NUMNP TMIN=MIN(TMIN,X(3,I)) TMAX=MAX(TMAX,X(3,I)) 30 CONTINUE TDIST=(TMAX-TMIN)/NUMNP*1.00001 C ATTENTION PAS DE CLASSEMENT DES ELEMENTS SELON LA DISTANCE C REMPLISSAGE DE INDZON NBZOR=NBZONE/NMULT+1 DO 23 I=1,NBZOR+1 INDZON(I)=0 23 CONTINUE DO 21 I=1,NUMNP NPOS=NPZON(I)/NMULT+1 IF (INDZON(NPOS).EQ.0) INDZON(NPOS)=I 21 CONTINUE INDZON(NBZOR+1)=NUMNP DO 35 I=1,NBZOR II=NBZOR+1-I IF (INDZON(II).EQ.0) INDZON(II)=INDZON(II+1) 35 CONTINUE C ON COMPTE LES ELEMENTS NELEM=0 IPT1=MELEME SEGACT MELEME DO 3101 IOB=1,MAX(1,LISOUS(/1)) IF (LISOUS(/1).NE.0) THEN IPT1=LISOUS(IOB) SEGACT IPT1 ENDIF NELEM=NELEM+IPT1.NUM(/2) 3101 CONTINUE C MODIF JUIN 1986 IVU=2 POUR LES POINTS DU PLAN DE COUPE (NON CACHABLE IF (MCOUP.NE.0) THEN IPT1=LISOUS(LISOUS(/1)) DO 3200 IEL=1,IPT1.NUM(/2) IF (MOD(MCOUP(IEL)/8,2).EQ.0) GOTO 3200 DO 3201 JNN=1,IPT1.NUM(/1) IVU(ICPR(IPT1.NUM(JNN,IEL)))=2 3201 CONTINUE 3200 CONTINUE ENDIF C C FIN DES PREPARATIFS C MAINTENANT ON PEUT TRAVAILLER C ON BOUCLE SUR LES ELEMENTS EN CACHANT TOUT CE QUI DOIT ETRE CACHE C C ON TRAVAILLE PAR OBJET ELEMENTAIRE IPT1=MELEME DO 3001 IOB=1,MAX(1,LISOUS(/1)) IF (LISOUS(/1).NE.0) IPT1=LISOUS(IOB) SEGACT IPT1 IF (KSURF(IPT1.ITYPEL).EQ.0) GOTO 3004 NBELEM=IPT1.NUM(/2) NBNN=IPT1.NUM(/1) NBFAC=LTEL(1,IPT1.ITYPEL) IAD=LTEL(2,IPT1.ITYPEL)-1 IF (NBFAC.EQ.0) GOTO 3004 DO 101 IFAC=1,NBFAC ITYP=LDEL(1,IAD+IFAC) NPFAC=KDFAC(1,ITYP) JAD=LDEL(2,IAD+IFAC)-1 IDEP=KDFAC(2,ITYP) IFEP=IDEP+3*(KDFAC(3,ITYP)-1) DO 1011 ITRIAN=IDEP,IFEP,3 IAFA=LFAC(JAD+KFAC(ITRIAN)) IBFA=LFAC(JAD+KFAC(ITRIAN+1)) ICFA=LFAC(JAD+KFAC(ITRIAN+2)) DO 110 IEL=1,NBELEM IA=ICPR(IPT1.NUM(IAFA,IEL)) IB=ICPR(IPT1.NUM(IBFA,IEL)) IC=ICPR(IPT1.NUM(ICFA,IEL)) IF (IVU(IA).LT.1.AND.IVU(IB).LT.1.AND.IVU(IC).LT.1) GOTO 110 XA=X(1,IA) XB=X(1,IB) XC=X(1,IC) YA=X(2,IA) YB=X(2,IB) YC=X(2,IC) IZA=INT(real((XA-XMIN)/XDEC))+1+INT(real((YA-YMIN)/YDEC))*NDEC IZB=INT(real((XB-XMIN)/XDEC))+1+INT(real((YB-YMIN)/YDEC))*NDEC IZC=INT(real((XC-XMIN)/XDEC))+1+INT((real(YC-YMIN)/YDEC))*NDEC ZA=X(3,IA) ZB=X(3,IB) ZC=X(3,IC) DISTLI=MIN(ZA,ZB,ZC) DISTLS=MAX(ZA,ZB,ZC) C DECOUPONS LE TRIANGLE IA,IB,IC EN BANDELETTES DE ZONES IRETO = 0 C CLASSEMENT DE A B C IF (IZA.LE.IZB) GOTO 200 IAUX=IZA IZA=IZB IZB=IAUX IAUX=IA IA=IB IB=IAUX 200 IF(IZA.LE.IZC) GOTO 201 IAUX=IZA IZA=IZC IZC=IAUX IAUX=IA IA=IC IC=IAUX C IA EST LE POINT LE PLUS BAS 201 MIA=(IZA-1)/NDEC NIA=IZA-MIA*NDEC MIB=(IZB-1)/NDEC NIB=IZB-MIB*NDEC MIC=(IZC-1)/NDEC NIC=IZC-MIC*NDEC IF (MIA.EQ.MIB) GOTO 300 IF (MIA.NE.MIC) GOTO 202 IAUX=IZB IZB=IZC IZC=IAUX IAUX=IB IB=IC IC=IAUX IAUX=MIC MIC=MIB MIB=IAUX IAUX=NIC NIC=NIB NIB=IAUX GOTO 300 C AUCUN POINT N'EST DANS LA MEME BANDE QUE A 202 PENTB=REAL(NIA-NIB)/REAL(MIA-MIB) PENTC=REAL(NIA-NIC)/REAL(MIA-MIC) IF (PENTC.GT.PENTB) GOTO 203 IAUX=IB IB=IC IC=IAUX IAUX=NIB NIB=NIC NIC=IAUX IAUX=MIB MIB=MIC MIC=IAUX AUX=PENTB PENTB=PENTC PENTC=AUX C IB EST MAITENANT LE POINT DE GAUCHE 203 IBAND=MIA-1 ISTEP=0 IFINB=MAX(MIB,MIC) NP=MIN(NIA,NIB) NG=MAX(NIA,NIC) NMP=MIN(NIB,NIC) NMG=MAX(NIB,NIC) IRETO = 1 C* ASSIGN 205 TO IRETO <- DELETED IF (MIB.EQ.MIC) GOTO 205 PENTA=REAL(NIB-NIC)/REAL(MIB-MIC) 205 IBAND=IBAND+1 IF(IBAND.GT.IFINB) GOTO 400 IF (IBAND.GT.MIB) GOTO 210 PZ=NIA+PENTB*REAL(IBAND-MIA) IZD=INT(PZ-ABS(PENTB)) IZD1=INT(PZ+1+ABS(PENTB)) IZD=MAX(IZD,NP) IF (IBAND.EQ.MIB.AND.ISTEP.EQ.0) GOTO 210 220 IF (IBAND.GT.MIC) GOTO 211 PZ=NIA+PENTC*REAL(IBAND-MIA) IZF=INT(PZ+1+ABS(PENTC)) IZF1=INT(PZ-ABS(PENTC)) IZF=MIN(IZF,NG) IF (IBAND.EQ.MIC.AND.ISTEP.EQ.0) GOTO 211 GOTO 1000 231 CONTINUE PZ=NIA+PENTC*REAL(IBAND-MIA) IZF=INT(PZ+1+ABS(PENTC)) IZF=MIN(IZF,NG) GOTO 230 210 IF (IBAND.EQ.MIC.AND.ISTEP.EQ.0) GOTO 231 ISTEP=1 PZ=NIB+PENTA*REAL(IBAND-MIB) IZD1S=INT(PZ+1+ABS(PENTA)) IZDS=INT(PZ-ABS(PENTA)) IZDS=MAX(IZDS,NMP) IF (IBAND.EQ.MIB) IZDS=MIN(IZD,IZDS) IF (IBAND.EQ.MIB) IZD1S=MAX(IZD1S,IZD1) IZD=IZDS IZD1=IZD1S GOTO 220 211 ISTEP=1 PZ=NIB+PENTA*REAL(IBAND-MIB) IZFS=INT(PZ+1+ABS(PENTA)) IZFS=MIN(IZFS,NMG) IZF1S=INT(PZ-ABS(PENTA)) IF (IBAND.EQ.MIC) IZFS=MAX(IZF,IZFS) IF (IBAND.EQ.MIC) IZF1S=MIN(IZF1,IZF1S) IZF=IZFS IZF1=IZF1S GOTO 1000 230 IZD1=IZF IZF1=IZD IRETO = 2 C* ASSIGN 400 TO IRETO <- DELETED GOTO 1000 300 IBAND=MIA-1 IF (MIA.EQ.MIC) GOTO 311 PENTA=REAL(NIA-NIC)/REAL(MIA-MIC) PENTB=REAL(NIB-NIC)/REAL(MIB-MIC) IRETO = 3 C* ASSIGN 301 TO IRETO <- DELETED NP=MIN(NIA,NIC) NG=MAX(NIB,NIC) 301 IBAND=IBAND+1 IF (IBAND.GT.MIC) GOTO 400 PZ=NIA+PENTA*REAL(IBAND-MIA) IZD=INT(PZ-ABS(PENTA)) IZD1=INT(PZ+1+ABS(PENTA)) IZD=MAX(NP,IZD) PZ=NIB+PENTB*REAL(IBAND-MIB) IZF=INT(PZ+1+ABS(PENTB)) IZF=MIN(NG,IZF) IZF1=INT(PZ-ABS(PENTB)) IF (IBAND.EQ.MIA) GOTO 410 GOTO 1000 410 IZD1=IZF IZF1=IZD GOTO 1000 311 IZD=NIA IZF=MAX(NIB,NIC) IBAND=MIA IRETO = 2 C* ASSIGN 400 TO IRETO <- DELETED GOTO 410 C BOUCLE SUR LES ZONES INTERNES AU TRIANGLE C ON VA TRAVAILLER DANS UNE BANDE ALLANT DE IZD A IZF. C ON SAIT PAR AILLEURS QUE LES ZONES COMPRISES ENTRE IZD1 ET IZF1 SONT C ENTIEREMENT INSCRITE DANS LE TRIANGLE IA IB IC 1000 CONTINUE IOORD=NDEC*IBAND IZD=IZD+IOORD IZD1=IZD1+IOORD IZF=IZF+IOORD IZF1=IZF1+IOORD INDDEB=INDZON(IZD/NMULT+1) ** write(6,*) 'tiret3 inddeb',inddeb C ON SE LIMITE A CHERCHER ENTRE LE DEBUT ET LA FIN DE LA BANDE C POINT DE LA BANDE DO 2000 INDDD=INDDEB,NBPOIN IPZO=NPZON(INDDD) * ON EST AVANT LE DEBUT IF (IPZO.LT.IZD) GOTO 2000 * ON EST APRES LA FIN IF (IPZO.GT.IZF) GOTO 2100 IPOINT=NPOIN(INDDD) IF (IVU(IPOINT).LT.1.OR.IVU(IPOINT).EQ.2) GOTO 2000 IF (IPOINT.EQ.IA) GOTO 2000 IF (IPOINT.EQ.IB) GOTO 2000 IF (IPOINT.EQ.IC) GOTO 2000 ZP=X(3,IPOINT) IF (ZP.LT.DISTLI) GOTO 2000 C DE TOUTE FACON LE POINT EST DEVANT LE TRIANGLE XP=X(1,IPOINT) YP=X(2,IPOINT) VAX=XP-XA VBX=XP-XB VCX=XP-XC VAY=YP-YA VBY=YP-YB VCY=YP-YC DC=VAX*VBY-VAY*VBX DA=VBX*VCY-VBY*VCX IF (DA*DC.LT.0.) GOTO 2000 DB=VCX*VAY-VCY*VAX IF (DA*DB.LT.0.) GOTO 2000 IF (DB*DC.LT.0.) GOTO 2000 IF (IPZO.GT.IZD1.AND.IPZO.LT.IZF1) GOTO 2020 C CAR LE POINT EST STRICTEMENT DANS LE TRIANGLE IF (ABS(DA).GT.SMO) GOTO 2020 IF (ABS(DB).GT.SMO) GOTO 2020 IF (ABS(DC).GT.SMO) GOTO 2020 IF (VAX*VBX.LT.-SMO) GOTO 2020 IF (VAX*VCX.LT.-SMO) GOTO 2020 IF (VAY*VCY.LT.-SMO) GOTO 2020 IF (VAY*VCY.LT.-SMO) GOTO 2020 IF (VBX*VCX.LT.-SMO) GOTO 2020 IF (VBY*VCY.LT.-SMO) GOTO 2020 GOTO 2000 C IPOINT EST INSCRIT A L'INTERIEUR DU TRIANGLE C EST-IL DERRIERE 2020 IF (ZP.GT.DISTLS) GOTO 2030 C ON FAIT UNE INTERPOLATION LINEAIRE ENTRE LES DISTANCES AUX SOMMETS C EN UTILISANT LES COORDONNEES BARYCENTRIQUES QUI NE SONT AUTRES C QUE LE RAPPORT DES SURFACES DES PETITS TRIANGLES AU GRAND S=DA+DB+DC IF (S.EQ.0.) GOTO 2000 DA=DA/S DB=DB/S DC=DC/S S=DA*ZA+DB*ZB+DC*ZC IF (S.GT.ZP) GOTO 2000 2030 CONTINUE C SOMME NOUS VRAIMENT A L'INTERIEUR DE LA FACE C METHODE : NOUS DEVONS APPARTENIR A UN NOMBRE PAIR DE TRIANGLE DE C L'ELEMENT IF (NPFAC.EQ.3) GOTO 2009 IDEDAN=1 DO 2200 ITRIA1=IDEP,IFEP,3 IF (ITRIA1.EQ.ITRIAN) GOTO 2200 IAJ=ICPR(IPT1.NUM(LFAC(JAD+KFAC(ITRIA1)),IEL)) IBJ=ICPR(IPT1.NUM(LFAC(JAD+KFAC(ITRIA1+1)),IEL)) ICJ=ICPR(IPT1.NUM(LFAC(JAD+KFAC(ITRIA1+2)),IEL)) XAJ=X(1,IAJ) XBJ=X(1,IBJ) XCJ=X(1,ICJ) YAJ=X(2,IAJ) YBJ=X(2,IBJ) YCJ=X(2,ICJ) VAX=XP-XAJ VBX=XP-XBJ VCX=XP-XCJ VAY=YP-YAJ VBY=YP-YBJ VCY=YP-YCJ DC=VAX*VBY-VAY*VBX DA=VBX*VCY-VBY*VCX IF (DA*DC.LT.0.) GOTO 2200 DB=VCX*VAY-VCY*VAX IF (DA*DB.LT.0.) GOTO 2200 IF (DB*DC.LT.0.) GOTO 2200 * ON EST DANS LE TRIANGLE IDEDAN=IDEDAN+1 2200 CONTINUE IF (MOD(IDEDAN,2).EQ.0) GOTO 2000 2009 IF (IVU(IPOINT).NE.2) IVU(IPOINT)=-IEL-IOB*NELEM 2000 CONTINUE 2100 CONTINUE GOTO (205,400,301),IRETO C* GOTO IRETO,(205,400,301) <- DELETED 400 CONTINUE C ON A FINI L'ANALYSE D'UN TRIANGLE 110 CONTINUE 1011 CONTINUE 101 CONTINUE 3004 CONTINUE 3001 CONTINUE SEGSUP TRAV C ON MET A 1 LES IVU 2 ( POINTS NON CACHABLES) DO 3400 II=1,IVU(/1) IF (IVU(II).EQ.2) IVU(II)=1 3400 CONTINUE END