C HEXA SOURCE JC220346 16/11/29 21:15:16 9221 C---------------------------------------------------------------------| C | SUBROUTINE HEXA(II,JJ,IF1,IF2,IGAGNE) C | C CETTE SUBROUTINE TENTE DE CREER UN HEXAEDRE A PARTIR | C DES QUADRANGLES IF1 ET IF2. | C - IGAGNE=1 EN CAS DE SUCCES | C - IGAGNE=0 EN CAS D'ECHEC | C | C---------------------------------------------------------------------| C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC TDEMAIT -INC PPARAM -INC CCOPTIO LOGICAL REPONS,FACET,SOLHEX,DIAGO,IN2,PLAN,dist2 nfcini=nfcmax nptini=nptmax NPTSAU=NPTMAX C * WRITE(6,1000) 1000 FORMAT(' ----->>> HEXA <<<-----') C * * (1) VERIFICATION DE l'ANGLE du diedre * ANG=TETA(IF1,if2,ii,jj) if (ang.lt.-1.d0) return if (ang.gt.1.d0) goto 1500 IGAGNE=0 L1=0 L2=0 N3=0 N4=0 N5=0 N6=0 IF5=0 IF6=0 ICTF=0 ICTV=0 J1=IPRED(IF1,II) J2=ISUCC(IF1,JJ) K1=ISUCC(IF2,II) K2=IPRED(IF2,JJ) C C RECHERCHE DE LA FACETTE IF3 C --------------------------- C IF3=IFACE3(K1,II,J1) * IF (IF3.NE.0) WRITE(6,1010)IF3 1010 FORMAT(' ** IF3=',I3) C IF (IF3.NE.0) THEN N3=1 IF(NFC(4,IF3).EQ.0) goto 1500 C LA FACETTE DOIT ETRE QUADRANGULAIRE * write (6,*) ' if3 trouvee',nfc(1,if3),nfc(2,if3), * # nfc(3,if3),nfc(4,if3) ANG=TETA(IF1,if3,j1,ii) if (ang.lt.-1.d0) then * write (6,*) ' hexa probleme 1' goto 1500 endif ANG=TETA(IF2,if3,ii,k1) if (ang.lt.-1.d0) then * write (6,*) ' hexa probleme 2' goto 1500 endif ENDIF IF(N3.NE.0) L1=ISUCC(IF3,J1) C C RECHERCHE DE LA FACETTE IF4 C --------------------------- C IF4=IFACE3(J2,JJ,K2) * IF (IF4.NE.0) WRITE(6,1020)IF4 1020 FORMAT(' ** IF4=',I3) C IF (IF4.NE.0) THEN N4=1 IF(NFC(4,IF4).EQ.0) goto 1500 C LA FACETTE DOIT ETRE QUADRANGULAIRE * write (6,*) ' if4 trouvee',nfc(1,if4),nfc(2,if4), * # nfc(3,if4),nfc(4,if4) ANG=TETA(IF1,if4,jj,j2) if (ang.lt.-1.d0) then * write (6,*) ' hexa probleme 3' goto 1500 endif ANG=TETA(IF2,if4,k2,jj) if (ang.lt.-1.d0) then * write (6,*) ' hexa probleme 4' goto 1500 endif ENDIF IF(N4.NE.0) L2=IPRED(IF4,J2) C C C RECHERCHE DE LA FACETTE IF5 C --------------------------- C IF(L1*L2.NE.0) THEN IF5=IFACE4(L1,L2,K1,K2) if (if5.lt.0) then * write (6,*) ' if5 trouve 3pts ' return endif ELSEIF(L1.NE.0) THEN IF5=IFACE3(L1,K1,K2) ELSEIF(L2.NE.0) THEN IF5=IFACE3(L2,K1,K2) ELSE ANG=TETA(IF2,NOISIN(K1,K2,IF2),K1,K2) IF (ANG.GT.-1.d0) IF5=NOISIN(K1,K2,IF2) endif * IF (IF5.NE.0) WRITE(6,1030)IF5 1030 FORMAT(' ** IF5=',I3) C IF (IF5.NE.0) THEN N5=1 IF(NFC(4,IF5).EQ.0) RETURN C LA FACETTE DOIT ETRE QUADRANGULAIRE * write (6,*) ' if5 trouvee',if5,nfc(1,if5),nfc(2,if5), * # nfc(3,if5),nfc(4,if5) ENDIF IF (IF5.NE.0) L1=ISUCC(IF5,K1) IF (IF5.NE.0) L2=IPRED(IF5,K2) C C C C RECHERCHE DE LA FACETTE IF6 C --------------------------- C IF(L1*l2.NE.0) THEN IF6=IFACE4(L1,L2,J1,J2) if (if6.lt.0) then * write (6,*) ' if6 trouve 3pts ' goto 1500 endif ELSEIF(L1.NE.0) THEN IF6=IFACE3(L1,J1,J2) ELSEIF(L2.NE.0) THEN IF6=IFACE3(L2,J1,J2) ELSE ANG=TETA(IF1,NOISIN(J2,J1,IF1),J2,J1) IF (ANG.GT.-1.d0) IF6=NOISIN(J2,J1,if1) endif * IF (IF6.NE.0) WRITE(6,1040)IF6 1040 FORMAT(' ** IF5=',I3) C IF (IF6.NE.0) THEN N6=1 IF(NFC(4,IF6).EQ.0) goto 1500 C LA FACETTE DOIT ETRE QUADRANGULAIRE * write (6,*) ' if6 trouvee',if6,nfc(1,if6),nfc(2,if6), * # nfc(3,if6),nfc(4,if6) ENDIF IF (IF6.NE.0) L1=IPRED(IF6,J1) IF (IF6.NE.0) L2=ISUCC(IF6,J2) C * * (1) VERIFICATION DES ANGLES AVEC FACETTES ADJACENTES (MINI 90) * IF (N6.EQ.0) THEN ANG=TETA(IF1,NOISIN(J2,J1,IF1),J2,J1) * write (6,*) 'facette 1 angle 1',ang IF (ANG.GT.-1.d0) goto 1500 ENDIF IF (N3.EQ.0) THEN ANG=TETA(IF1,NOISIN(J1,II,IF1),J1,II) * write (6,*) 'facette 1 angle 2',ang IF (ANG.GT.-1.d0) goto 1500 ENDIF IF (N4.EQ.0) THEN ANG=TETA(IF1,NOISIN(JJ,J2,IF1),JJ,J2) * write (6,*) 'facette 1 angle 3',ang IF (ANG.GT.-1.d0) goto 1500 ENDIF IF (N5.EQ.0) THEN ANG=TETA(IF2,NOISIN(K1,K2,IF2),K1,K2) * write (6,*) 'facette 2 angle 1',ang IF (ANG.GT.-1.d0) goto 1500 ENDIF IF (N4.EQ.0) THEN ANG=TETA(IF2,NOISIN(K2,JJ,IF2),K2,JJ) * write (6,*) 'facette 2 angle 2',ang IF (ANG.GT.-1.d0) goto 1500 ENDIF IF (N3.EQ.0) THEN ANG=TETA(IF2,NOISIN(II,K1,IF2),II,K1) * write (6,*) 'facette 2 angle 3',ang IF (ANG.GT.-1.d0) goto 1500 ENDIF * a ameliorer plus tard C C C Construction des deux points supplementaires (si necessaire) C IF(L1.EQ.0) THEN L1=NPTMAX+1 NPTMAX=L1 C C DETERMINATION DES COORDONNEES DU NOUVEAU POINT L1 C ------------------------------------------------- xj1=xyz(1,j1)-xyz(1,ii) yj1=xyz(2,j1)-xyz(2,ii) zj1=xyz(3,j1)-xyz(3,ii) vj1=xj1**2+yj1**2+zj1**2 xk1=xyz(1,k1)-xyz(1,ii) yk1=xyz(2,k1)-xyz(2,ii) zk1=xyz(3,k1)-xyz(3,ii) vk1=xk1**2+yk1**2+zk1**2 scal=xj1*xk1+yj1*yk1+zj1*zk1 xyz(1,l1)=(xyz(1,k1)+xj1-scal*xk1/vk1+ * xyz(1,j1)+xk1-scal*xj1/vj1)*0.5d0 xyz(2,l1)=(xyz(2,k1)+yj1-scal*yk1/vk1+ * xyz(2,j1)+yk1-scal*yj1/vj1)*0.5d0 xyz(3,l1)=(xyz(3,k1)+zj1-scal*zk1/vk1+ * xyz(3,j1)+zk1-scal*zj1/vj1)*0.5d0 DO 150 I=1,3 ****** XYZ(I,L1)=(XYZ(I,L1)+XYZ(I,J1)+XYZ(I,K1)-XYZ(I,II))*0.5d0 XYZ(I,L1)= XYZ(I,J1)+XYZ(I,K1)-XYZ(I,II) 150 CONTINUE if (L2.ne.0) then do 151 i=1,3 XYZ(I,L1)=(XYZ(I,L1)+XYZ(I,J1)+XYZ(I,L2)-XYZ(I,J2)+ * XYZ(I,K1)+XYZ(I,L2)-XYZ(I,K2))/3.d0 151 continue endif XYZ(4,L1)=(XYZ(4,J1)+XYZ(4,K1)+XYZ(4,II))/3.d0 CALL DIST(L1,KP,GL,IOK,II,JJ,J1,J2,K1,K2,L2,0,0,0) IF (IOK.EQ.0) THEN * C'est rate nptmax=nptsau goto 1500 C l1=kp C write (6,*) 'hexa point assimile 1 ',kp C if (diago(l1,j1,0.95D0)) then C* write (6,*) ' diago quad',l1,j1 C NPTMAX=NPTsau C goto 1500 C endif C if (diago(l1,k1,0.95D0)) then C* write (6,*) ' diago quad',l1,k1 C NPTMAX=NPTsau C goto 1500 C endif ENDIF * if (dist2(l1)) then * nptmax=nptsau * return * endif XYZ(1,L1+1)=XYZ(1,JJ)+1.35*(XYZ(1,L1)-XYZ(1,JJ)) XYZ(2,L1+1)=XYZ(2,JJ)+1.35*(XYZ(2,L1)-XYZ(2,JJ)) XYZ(3,L1+1)=XYZ(3,JJ)+1.35*(XYZ(3,L1)-XYZ(3,JJ)) XYZ(4,L1+1)=XYZ(4,JJ)+1.35*(XYZ(4,L1)-XYZ(4,JJ)) * IF (.NOT.IN2(jj,L1+1,nfcini)) THEN * write (6,*) ' in incorrect 1' * nptmax=nptsau * goto 1500 * ENDIF * if (dist2(l1+1)) then * nptmax=nptsau * return * endif CALL DIST(L1+1,KP,GL,IOK,II,JJ,J1,J2,K1,K2,L1,L2,0,0) IF (IOK.EQ.0) THEN * C'est encore rate NPTMAX=NPTSAU IF (IVERB.EQ.1) write (6,*) 'point mal place 1.1',kp goto 1500 ENDIF ENDIF * verif que l1 appartient bien a n3 n5 et n6 si ils existent if (n3.ne.0) then if (L1.ne.ISUCC(IF3,J1)) then * write (6,*) ' faces ne correspondent pas ' nptmax=nptsau goto 1500 endif endif if (n5.ne.0) then if (L1.ne.ISUCC(IF5,K1)) then * write (6,*) ' faces ne correspondent pas ' nptmax=nptsau goto 1500 endif endif if (n6.ne.0) then if (L1.ne.IPRED(IF6,J1)) then * write (6,*) ' faces ne correspondent pas ' nptmax=nptsau goto 1500 endif endif IF(L2.EQ.0) THEN L2=NPTMAX+1 NPTMAX=L2 C C DETERMINATION DES COORDONNEES DU NOUVEAU POINT L2 C ------------------------------------------------- xj2=xyz(1,j2)-xyz(1,jj) yj2=xyz(2,j2)-xyz(2,jj) zj2=xyz(3,j2)-xyz(3,jj) vj2=xj2**2+yj2**2+zj2**2 xk2=xyz(1,k2)-xyz(1,jj) yk2=xyz(2,k2)-xyz(2,jj) zk2=xyz(3,k2)-xyz(3,jj) vk2=xk2**2+yk2**2+zk2**2 scal=xj2*xk2+yj2*yk2+zj2*zk2 xyz(1,l2)=(xyz(1,k2)+xj2-scal*xk2/vk2+ * xyz(1,j2)+xk2-scal*xj2/vj2)*0.5d0 xyz(2,l2)=(xyz(2,k2)+yj2-scal*yk2/vk2+ * xyz(2,j2)+yk2-scal*yj2/vj2)*0.5d0 xyz(3,l2)=(xyz(3,k2)+zj2-scal*zk2/vk2+ * xyz(3,j2)+zk2-scal*zj2/vj2)*0.5d0 DO 170 I=1,3 ********** XYZ(I,L2)=(XYZ(I,L2)+XYZ(I,J2)+XYZ(I,K2)-XYZ(I,JJ))*0.5d0 XYZ(I,L2)= XYZ(I,J2)+XYZ(I,K2)-XYZ(I,JJ) 170 CONTINUE if (L1.ne.0) then do 171 i=1,3 XYZ(I,L2)=(XYZ(I,L2)+XYZ(I,J2)+XYZ(I,L1)-XYZ(I,J1)+ * XYZ(I,L1)+XYZ(I,K2)-XYZ(I,K1))/3.d0 171 continue endif XYZ(4,L2)=(XYZ(4,J2)+XYZ(4,K2)+XYZ(4,JJ))/3.d0 CALL DIST(L2,KP,GL,IOK,II,JJ,J1,J2,K1,K2,L1,0,0,0) IF (IOK.EQ.0) THEN * C'est rate nptmax=nptsau goto 1500 C l2=kp C write (6,*) 'hexa point assimile 2 ',kp C if (diago(l2,j2,0.95d0)) then C* write (6,*) ' diago quad',l1,j1 C NPTMAX=NPTsau C goto 1500 C endif C if (diago(l2,k2,0.95d0)) then C* write (6,*) ' diago quad',l2,k2 C NPTMAX=NPTsau C goto 1500 C endif C if (diago(l2,l1,0.95d0)) then C* write (6,*) ' diago quad',l2,l1 C NPTMAX=NPTsau C goto 1500 C endif ENDIF * if (dist2(l2)) then * nptmax=nptsau * return * endif XYZ(1,L2+1)=XYZ(1,II)+1.35*(XYZ(1,L2)-XYZ(1,II)) XYZ(2,L2+1)=XYZ(2,II)+1.35*(XYZ(2,L2)-XYZ(2,II)) XYZ(3,L2+1)=XYZ(3,II)+1.35*(XYZ(3,L2)-XYZ(3,II)) XYZ(4,L2+1)=XYZ(4,II)+1.35*(XYZ(4,L2)-XYZ(4,II)) * IF (.NOT.IN2(ii,L2+1,nfcini)) THEN * write (6,*) ' in incorrect 2' * nptmax=nptsau * goto 1500 * ENDIF * if (dist2(l2+1)) then * nptmax=nptsau * return * endif CALL DIST(L2+1,KP,GL,IOK,II,JJ,J1,J2,K1,K2,L1,L2,0,0) IF (IOK.EQ.0) THEN ** C'est encore rate NPTMAX=NPTSAU IF (IVERB.EQ.1) write (6,*) 'point mal place 2.2',kp goto 1500 ENDIF ENDIF * verif que l2 appartient bien a n4 n5 et n6 si ils existent if (n4.ne.0) then if (L2.ne.IPRED(IF4,J2)) then * write (6,*) ' faces ne correspondent pas ' nptmax=nptsau goto 1500 endif endif if (n5.ne.0) then if (L2.ne.IPRED(IF5,K2)) then * write (6,*) ' faces ne correspondent pas ' nptmax=nptsau goto 1500 endif endif if (n6.ne.0) then if (L2.ne.ISUCC(IF6,J2)) then * write (6,*) ' faces ne correspondent pas ' nptmax=nptsau goto 1500 endif endif IF (DIAGO(L1,L2,0.95D0)) THEN * write (6,*) ' diago quad',L1,L2 NPTMAX=NPTsau goto 1500 ENDIF C C C C CONSTRUCTION DU CUBE C -------------------- C IF (IF3.EQ.0) THEN C C CREATION DE LA FACETTE IF3 C -------------------------- IF3=IFACE4(IPRED(IF1,II),II,ISUCC(IF2,II),L1) if (if3.lt.0) then * write (6,*) ' if3 trouve 3pts ' nfcmax=nfcini nptmax=nptini return elseif (if3.eq.0) then NFCMAX=NFCMAX+1 IF3=NFCMAX ICTF=ICTF+1 C NFC(1,IF3)=IPRED(IF1,II) NFC(2,IF3)=II NFC(3,IF3)=ISUCC(IF2,II) NFC(4,IF3)=L1 C ENDIF ENDIF C C IF (IF4.EQ.0) THEN C C CREATION DE LA FACETTE IF4 C -------------------------- IF4=IFACE4(JJ,ISUCC(IF1,JJ),L2,IPRED(IF2,JJ)) if (if4.lt.0) then * write (6,*) ' if4 trouve 3pts ' nfcmax=nfcini nptmax=nptini return elseif (if4.eq.0) then NFCMAX=NFCMAX+1 IF4=NFCMAX ICTF=ICTF+1 C NFC(1,IF4)=JJ NFC(2,IF4)=ISUCC(IF1,JJ) NFC(3,IF4)=L2 NFC(4,IF4)=IPRED(IF2,JJ) C ENDIF ENDIF C C IF (IF5.EQ.0) THEN C C CREATION DE LA FACETTE IF5 C -------------------------- IF5=IFACE4(l1,k1,k2,l2) if (if5.lt.0) then * write (6,*) ' if5 trouve 3pts ' nfcmax=nfcini nptmax=nptini return elseif (if5.eq.0) then NFCMAX=NFCMAX+1 IF5=NFCMAX ICTF=ICTF+1 C NFC(1,IF5)=L1 NFC(2,IF5)=K1 NFC(3,IF5)=K2 NFC(4,IF5)=L2 C C ENDIF ENDIF IF (IF6.EQ.0) THEN C C CREATION DE LA FACETTE IF6 C -------------------------- IF6=IFACE4(j1,l1,l2,j2) if (if6.lt.0) then * write (6,*) ' if6 trouve 3pts ' nfcmax=nfcini nptmax=nptini return elseif (if6.eq.0) then NFCMAX=NFCMAX+1 IF6=NFCMAX ICTF=ICTF+1 C NFC(1,IF6)=J1 NFC(2,IF6)=L1 NFC(3,IF6)=L2 NFC(4,IF6)=J2 C C ENDIF ENDIF C C ON ENLEVE LES FACETTES IF1, IF2 ET IF3 C -------------------------------------- CALL REPSUB(IF1) CALL REPSUB(IF2) CALL REPSUB(IF3) CALL REPSUB(IF4) CALL REPSUB(IF5) CALL REPSUB(IF6) C C LE VOLUME CREE EST-IL VALIDE ? C ------------------------------ IF (.NOT.PLAN(IF3)) GOTO 160 * write (6,*) ' plan(if3) passe' IF (.NOT.PLAN(IF4)) GOTO 160 * write (6,*) ' plan(if4) passe' IF (.NOT.PLAN(IF5)) GOTO 160 * write (6,*) ' plan(if5) passe' IF (.NOT.PLAN(IF6)) GOTO 160 * write (6,*) ' plan(if6) passe' IF (.NOT.FACET(IF3)) GOTO 160 * write (6,*) ' facet(if3) passe' IF (.NOT.FACET(IF4)) GOTO 160 * write (6,*) ' facet(if4) passe' IF (.NOT.FACET(IF5)) GOTO 160 * write (6,*) ' facet(if5) passe' IF (.NOT.FACET(IF6)) GOTO 160 * write (6,*) ' facet(if6) passe' IF (.NOT.SOLHEX(IF1,IF2,IF3,IF4,IF5,IF6)) GOTO 160 * write (6,*) ' solhex passe' * * VERIFICATION TAILLE IF (N3.EQ.0.AND.N5.EQ.0) THEN KF1=IPRED(IF1,II) KF2=ISUCC(IF2,II) DNORM=(XYZ(1,KF1)-XYZ(1,KF2))**2 # +(XYZ(2,KF1)-XYZ(2,KF2))**2 # +(XYZ(3,KF1)-XYZ(3,KF2))**2 DTEST=tcrit*XYZ(4,KF1)*XYZ(4,KF2) IF (DNORM.GT.DTEST) GOTO 160 ENDIF IF (N4.EQ.0.AND.N5.EQ.0) THEN KF1=IPRED(IF2,JJ) KF2=ISUCC(IF1,JJ) DNORM=(XYZ(1,KF1)-XYZ(1,KF2))**2 # +(XYZ(2,KF1)-XYZ(2,KF2))**2 # +(XYZ(3,KF1)-XYZ(3,KF2))**2 DTEST=tcrit*XYZ(4,KF1)*XYZ(4,KF2) IF (DNORM.GT.DTEST) GOTO 160 ENDIF * write (6,*) 'hexa a cree un cube numero nfacet ',nvol+1,nfacet C C LE VOLUME CREE EST VALIDE | C --------------------------- C MEMORISATION DU VOLUME IF1, IF2, IF3, IF4 ET IF5 C ------------------------------------------------ NVOL=NVOL+1 IF (NFV(1,IF1).EQ.0) NFV(1,IF1)=NVOL IF (NFV(1,IF1).NE.NVOL) NFV(2,IF1)=NVOL IF (NFV(1,IF2).EQ.0) NFV(1,IF2)=NVOL IF (NFV(1,IF2).NE.NVOL) NFV(2,IF2)=NVOL IF (NFV(1,IF3).EQ.0) NFV(1,IF3)=NVOL IF (NFV(1,IF3).NE.NVOL) NFV(2,IF3)=NVOL IF (NFV(1,IF4).EQ.0) NFV(1,IF4)=NVOL IF (NFV(1,IF4).NE.NVOL) NFV(2,IF4)=NVOL IF (NFV(1,IF5).EQ.0) NFV(1,IF5)=NVOL IF (NFV(1,IF5).NE.NVOL) NFV(2,IF5)=NVOL IF (NFV(1,IF6).EQ.0) NFV(1,IF6)=NVOL IF (NFV(1,IF6).NE.NVOL) NFV(2,IF6)=NVOL IVOL(9,NVOL)=20 IVOL(1,NVOL)=II IVOL(2,NVOL)=JJ IVOL(3,NVOL)=J2 IVOL(4,NVOL)=J1 IVOL(5,NVOL)=K1 IVOL(6,NVOL)=K2 IVOL(7,NVOL)=L2 IVOL(8,NVOL)=L1 C * WRITE(6,1100)NVOL,(IVOL(I,NVOL),I=1,9) if (iimpi.eq.1) write (6,1100) nfacet,(ivol(i,nvol),i=1,8) 1100 FORMAT(' HEXA facettes ',i5,' cub8 ',8i5) C * DO 150 J=1,NPTMAX * WRITE(6,1110)J,(NPF(I,J),I=1,40) *1110 FORMAT(I4,4X,40I3) *150 CONTINUE C IGAGNE=1 C RETURN C 160 CONTINUE * write (6,*) ' probleme en validant le cube ' C C LE VOLUME CREE N'EST PAS VALIDE: IL FAUT DONC DETRUIRE LES FACETT C CREEES. --------------------------------------------------------- CALL REPSUB(IF1) CALL REPSUB(IF2) CALL REPSUB(IF3) CALL REPSUB(IF4) CALL REPSUB(IF5) CALL REPSUB(IF6) C NFCMAX=NFCMAX-ICTF C NPTMAX=NPTSAU 1500 continue * maintenant on essaye un prisme NPTMAX=NPTSAU * call prism1(II,JJ,IF1,IF2,IGAGNE) return end