C PYRA3 SOURCE PV 07/10/05 21:16:25 5890 C---------------------------------------------------------------------| C | SUBROUTINE PYRA3(II,JJ,IF1,IF2,IGAGNE) C | C CETTE SUBROUTINE TENTE DE CREER DEUX PYRAMIDES 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 CCOPTIO LOGICAL REPONS,FACET,SOLPYR,SOLTET,DIAGO,IN2,VERDIV,IN C * WRITE(6,1000) 1000 FORMAT(' ----->>> PYRA3 <<<-----') C nptini=nptmax nfcini=nfcmax ICTF=0 ICTV=0 IP1=IPRED(IF1,II) IP2=IPRED(IF1,IP1) JP1=ISUCC(IF2,II) JP2=ISUCC(IF2,JP1) * CREATION POINT MILIEU NPTMAX=NPTMAX+1 XYZ(4,NPTMAX)=(XYZ(4,IP1)+XYZ(4,IP2)+XYZ(4,JP1)+XYZ(4,JP2))/4. * deplacement du point pour l'eloigner de ii jj xn1=(xyz(2,jj)-xyz(2,ii))*((xyz(3,ip1)+xyz(3,ip2))/2-xyz(3,ii))- > (xyz(3,jj)-xyz(3,ii))*((xyz(2,ip1)+xyz(2,ip2))/2-xyz(2,ii)) yn1=(xyz(3,jj)-xyz(3,ii))*((xyz(1,ip1)+xyz(1,ip2))/2-xyz(1,ii))- > (xyz(1,jj)-xyz(1,ii))*((xyz(3,ip1)+xyz(3,ip2))/2-xyz(3,ii)) zn1=(xyz(1,jj)-xyz(1,ii))*((xyz(2,ip1)+xyz(2,ip2))/2-xyz(2,ii))- > (xyz(2,jj)-xyz(2,ii))*((xyz(1,ip1)+xyz(1,ip2))/2-xyz(1,ii)) sn1=sqrt(xn1**2+yn1**2+zn1**2) xn1=xn1/sn1 yn1=yn1/sn1 zn1=zn1/sn1 xn2=((xyz(2,jp1)+xyz(2,jp2))/2-xyz(2,ii))*(xyz(3,jj)-xyz(3,ii))- > ((xyz(3,jp1)+xyz(3,jp2))/2-xyz(3,ii))*(xyz(2,jj)-xyz(2,ii)) yn2=((xyz(3,jp1)+xyz(3,jp2))/2-xyz(3,ii))*(xyz(1,jj)-xyz(1,ii))- > ((xyz(1,jp1)+xyz(1,jp2))/2-xyz(1,ii))*(xyz(3,jj)-xyz(3,ii)) zn2=((xyz(1,jp1)+xyz(1,jp2))/2-xyz(1,ii))*(xyz(2,jj)-xyz(2,ii))- > ((xyz(2,jp1)+xyz(2,jp2))/2-xyz(2,ii))*(xyz(1,jj)-xyz(1,ii)) sn2=sqrt(xn2**2+yn2**2+zn2**2) xn2=xn2/sn2 yn2=yn2/sn2 zn2=zn2/sn2 xn=(xn1+xn2)/2 yn=(yn1+yn2)/2 zn=(zn1+zn2)/2 sn=sqrt(xn**2+yn**2+zn**2) xn=xn/sn yn=yn/sn zn=zn/sn * xmil=(xyz(1,ii)+xyz(1,jj))/2 * ymil=(xyz(2,ii)+xyz(2,jj))/2 * zmil=(xyz(3,ii)+xyz(3,jj))/2 * xv=xyz(1,jj)-xyz(1,ii) yv=xyz(2,jj)-xyz(2,ii) zv=xyz(3,jj)-xyz(3,ii) sv=sqrt(xv**2+yv**2+zv**2) xv=xv/sv yv=yv/sv zv=zv/sv xli1=xv*(xyz(1,ip1)-xyz(1,ii))+yv*(xyz(2,ip1)-xyz(2,ii))+ > zv*(xyz(3,ip1)-xyz(3,ii)) xli2=xv*(xyz(1,ip2)-xyz(1,ii))+yv*(xyz(2,ip2)-xyz(2,ii))+ > zv*(xyz(3,ip2)-xyz(3,ii)) xlj1=xv*(xyz(1,jp1)-xyz(1,ii))+yv*(xyz(2,jp1)-xyz(2,ii))+ > zv*(xyz(3,jp1)-xyz(3,ii)) xlj2=xv*(xyz(1,jp2)-xyz(1,ii))+yv*(xyz(2,jp2)-xyz(2,ii))+ > zv*(xyz(3,jp2)-xyz(3,ii)) xl=(xli1+xli2+xlj1+xlj2+2*sv+2*0)/8 xl=0.5*sv xmil=xyz(1,ii)+xl*xv ymil=xyz(2,ii)+xl*yv zmil=xyz(3,ii)+xl*zv expf = xyz(4,nptmax) xyz(1,nptmax)=xmil+xn*expf*expfac xyz(2,nptmax)=ymil+yn*expf*expfac xyz(3,nptmax)=zmil+zn*expf*expfac * write (6,*) ' pyra3 creation de 2 elements et pt ',nptmax * write (6,*) (xyz(i,ii),i=1,4) * write (6,*) (xyz(i,jj),i=1,4) * write (6,*) (xyz(i,ip1),i=1,4) * write (6,*) (xyz(i,ip2),i=1,4) * write (6,*) (xyz(i,jp1),i=1,4) * write (6,*) (xyz(i,jp2),i=1,4) * write (6,*) (xyz(i,nptmax),i=1,4) * CREATION DES ELEMENTS IPTT=NPTMAX CALL DIST(iptt,nptaux,GL,IOK,ii,jj,ip1,ip2,jp1,jp2,0,0,0,0) IF (IOK.EQ.0) THEN NPTMAX=nptini WRITE (6,*) ' pyra3 DIST ',nptaux return ENDIF IF (gl.lt.xyz(4,iptt)/4) then write (6,*) 'pyra3 GL-1' nptmax=nptini return endif IF (.NOT.IN2(ii,IPTT,nfcini).or..NOT.IN2(jj,IPTT,nfcini)) THEN write (6,*) ' in2 echec ' NPTMAX=nptini return ENDIF * creations des faces NFCMAX=NFCMAX+1 IF3=NFCMAX NFC(1,IF3)=II NFC(2,IF3)=iptt NFC(3,IF3)=IP1 NFC(4,IF3)=0 * NFCMAX=NFCMAX+1 IF4=NFCMAX NFC(1,IF4)=IP1 NFC(2,IF4)=IPTT NFC(3,IF4)=IP2 NFC(4,IF4)=0 * NFCMAX=NFCMAX+1 IF5=NFCMAX NFC(1,IF5)=IP2 NFC(2,IF5)=IPTT NFC(3,IF5)=JJ NFC(4,IF5)=0 * la face commune NFCMAX=NFCMAX+1 IF6=NFCMAX NFC(1,IF6)=JJ NFC(2,IF6)=IPTT NFC(3,IF6)=II NFC(4,IF6)=0 * NFCMAX=NFCMAX+1 IF7=NFCMAX NFC(1,IF7)=JJ NFC(2,IF7)=iptt NFC(3,IF7)=jp2 NFC(4,IF7)=0 * NFCMAX=NFCMAX+1 IF8=NFCMAX NFC(1,IF8)=jp2 NFC(2,IF8)=iptt NFC(3,IF8)=jp1 NFC(4,IF8)=0 * NFCMAX=NFCMAX+1 IF9=NFCMAX NFC(1,IF9)=jp1 NFC(2,IF9)=iptt NFC(3,IF9)=II NFC(4,IF9)=0 * si necessaire verification diago IF (DIAGO(iptt,ii,diacri)) GOTO 275 IF (DIAGO(iptt,jj,diacri)) GOTO 275 IF (DIAGO(iptt,ip1,diacri)) GOTO 275 IF (DIAGO(iptt,ip2,diacri)) GOTO 275 IF (DIAGO(iptt,jp1,diacri)) GOTO 275 IF (DIAGO(iptt,jp2,diacri)) GOTO 275 goto 276 275 continue nptmax=nptini nfcmax=nfcini write (6,*) ' pyra3 echec diago' return 276 continue C CALL REPSUB(IF1) CALL REPSUB(IF3) CALL REPSUB(IF4) CALL REPSUB(IF5) CALL REPSUB(IF6) * verification du premier element element IF (.NOT.SOLPYR(IF1,IF3,IF4,IF5,IF6)) then write (6,*) 'pyra3-1 solpyr invalide' GOTO 160 endif IF (.NOT.FACET(IF3)) then write(6,*) ' pyra3 facet if3 invalide' GOTO 160 ENDIF IF (.NOT.FACET(IF4)) then write(6,*) ' pyra3 facet if4 invalide' GOTO 160 ENDIF IF (.NOT.FACET(IF5)) then write(6,*) ' pyra3 facet if5 invalide' GOTO 160 ENDIF CALL REPSUB(IF2) CALL REPSUB(IF6) CALL REPSUB(IF7) CALL REPSUB(IF8) CALL REPSUB(IF9) * verification du deuxieme element element IF (.NOT.SOLPYR(IF2,IF6,IF7,IF8,IF9)) then write (6,*) 'pyra3-2 solpyr invalide' GOTO 165 endif IF (.NOT.FACET(IF7)) then write(6,*) ' pyra3 facet if7 invalide' GOTO 165 ENDIF IF (.NOT.FACET(IF8)) then write(6,*) ' pyra1 facet if8 invalide' GOTO 165 ENDIF IF (.NOT.FACET(IF9)) then write(6,*) ' pyra1 facet if9 invalide' GOTO 165 ENDIF * OK pour creation elements 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,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)=35 DO 280 I=1,4 IVOL(I,NVOL)=NFC(I,IF1) 280 CONTINUE IVOL(5,NVOL)=iptt if (iimpi.eq.1) write (6,1102) nfacet,(ivol(i,nvol),i=1,5) 1102 FORMAT(' PYRA3-1 facettes ',i5,' pyr5 ',8i5) NVOL=NVOL+1 IVOL(9,NVOL)=35 IF (NFV(1,IF2).EQ.0) NFV(1,IF2)=NVOL IF (NFV(1,IF2).NE.NVOL) NFV(2,IF2)=NVOL IF (NFV(1,IF6).EQ.0) NFV(1,IF6)=NVOL IF (NFV(1,IF6).NE.NVOL) NFV(2,IF6)=NVOL IF (NFV(1,IF7).EQ.0) NFV(1,IF7)=NVOL IF (NFV(1,IF7).NE.NVOL) NFV(2,IF7)=NVOL IF (NFV(1,IF8).EQ.0) NFV(1,IF8)=NVOL IF (NFV(1,IF8).NE.NVOL) NFV(2,IF8)=NVOL IF (NFV(1,IF9).EQ.0) NFV(1,IF9)=NVOL IF (NFV(1,IF9).NE.NVOL) NFV(2,IF9)=NVOL DO 282 I=1,4 IVOL(I,NVOL)=NFC(I,IF2) 282 CONTINUE IVOL(5,NVOL)=iptt if (iimpi.eq.1) write (6,1101) nfacet,(ivol(i,nvol),i=1,5) 1101 FORMAT(' PYRA3-2 facettes ',i5,' pyr5 ',8i5) IGAGNE=1 RETURN 165 continue CALL REPSUB(IF9) CALL REPSUB(IF8) CALL REPSUB(IF7) CALL REPSUB(IF6) CALL REPSUB(IF2) 160 continue CALL REPSUB(IF6) CALL REPSUB(IF5) CALL REPSUB(IF4) CALL REPSUB(IF3) CALL REPSUB(IF1) nptmax=nptini nfcmax=nfcini goto 9000 9000 CONTINUE END