Sources Esope | Cast3M

solpyr
C SOLPYR    SOURCE    PV        17/09/12    21:15:03     9542           
C---------------------------------------------------------------------|
C                                                                     |
       LOGICAL FUNCTION SOLPYR(IF1,IF2,IF3,IF4,IF5)
C                                                                     |
C      CETTE FONCTION LOGIQUE TESTE SI LA PYRAMIDE DECRIT PAR LES     |
C      FACETTES IF1..IF5 EST VIDE DE POINT.                           |
C      LES ENTIERS N1..N5 INDIQUENT L'ORIENTATION DES FACETTES        |
C      SOLPYR EST VRAI SI LA PYRAMIDE EST VIDE (DONC VALIDE)          |
C      SOLPYR EST FAUX SI LA PYRAMIDE CONTIENT UN POINT (ET EST       |
C      DONC INVALIDE)                                                 |
C                                                                     |
C---------------------------------------------------------------------|
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC TDEMAIT
      LOGICAL FACET3,FACET2,vervol
      dimension ift(5)
      ift(1)=if1
      ift(2)=if2
      ift(3)=if3
      ift(4)=if4
      ift(5)=if5
C
C   CHERCHER SI LES FACETTES SONT NOUVELLES OU ANCIENNES
        N1=-1
        N2=-1
        N3=-1
        N4=-1
        N5=-1
         IF (IFAT(IF1).NE.0) N1=1
         IF (IFAT(IF2).NE.0) N2=1
         IF (IFAT(IF3).NE.0) N3=1
         IF (IFAT(IF4).NE.0) N4=1
         IF (IFAT(IF5).NE.0) N5=1
  10    CONTINUE
*  encadrer l'hexaedre
       xmin=1e30
       xmax=-1e30
       ymin=1e30
       ymax=-1e30
       zmin=1e30
       zmax=-1e30
       do 92 j=1,5
       do 90 i=1,3
       xmin=min(xmin,xyz(1,nfc(i,ift(j))))
       xmax=max(xmax,xyz(1,nfc(i,ift(j))))
       ymin=min(ymin,xyz(2,nfc(i,ift(j))))
       ymax=max(ymax,xyz(2,nfc(i,ift(j))))
       zmin=min(zmin,xyz(3,nfc(i,ift(j))))
       zmax=max(zmax,xyz(3,nfc(i,ift(j))))
  90   continue
  92   continue
       iteste=0
*      DO 100 I=1,NPTMAX
       DO 100 I=1,0
          IF (NPF(1,I).EQ.0) GOTO 100
          DO 110 J=1,4
             IF (I.EQ.NFC(J,IF1)) GOTO 100
             IF (I.EQ.NFC(J,IF2)) GOTO 100
             IF (I.EQ.NFC(J,IF3)) GOTO 100
             IF (I.EQ.NFC(J,IF4)) GOTO 100
             IF (I.EQ.NFC(J,IF5)) GOTO 100
110       CONTINUE
       if (iteste.ne.0) then
       if (xyz(1,i).lt.xmin.or.xyz(1,i).gt.xmax) goto 100
       if (xyz(2,i).lt.ymin.or.xyz(2,i).gt.ymax) goto 100
       if (xyz(3,i).lt.zmin.or.xyz(3,i).gt.zmax) goto 100
       endif
       iteste=1
C
       V1=N1*VOL(I,NFC(1,IF1),NFC(2,IF1),NFC(3,IF1))
       V2=N1*VOL(I,NFC(1,IF1),NFC(3,IF1),NFC(4,IF1))
       V3=N2*VOL(I,NFC(1,IF2),NFC(2,IF2),NFC(3,IF2))
       V4=N3*VOL(I,NFC(1,IF3),NFC(2,IF3),NFC(3,IF3))
       V5=N4*VOL(I,NFC(1,IF4),NFC(2,IF4),NFC(3,IF4))
       V6=N5*VOL(I,NFC(1,IF5),NFC(2,IF5),NFC(3,IF5))
       VV=V1+V2+V3+V4+V5+V6
       IF (VV.LE.0.) then
        if(IVERB.EQ.1) write (6,*) ' pyramide volume negatif '
        solpyr=.false.
        return
        GOTO 120
       endif
       IF (V1.LE.-0.01*VV) GOTO 100
       IF (V2.LE.-0.01*VV) GOTO 100
       IF (V3.LE.-0.01*VV) GOTO 100
       IF (V4.LE.-0.01*VV) GOTO 100
       IF (V5.LE.-0.01*VV) GOTO 100
       IF (V6.LE.-0.01*VV) GOTO 100
       GOTO 120
C
C
100    CONTINUE
C
C      IL N'EXISTE PAS DE POINTS INTERNES AU VOLUME
       SOLPYR=.TRUE.
C  RAJOUT PV TEST INTERSECTION DES FACETTES
       IF (N1.EQ.1) SOLPYR=SOLPYR.AND.FACET2(IF1)
          IF (.NOT.SOLPYR) RETURN
       IF (N2.EQ.1) SOLPYR=SOLPYR.AND.FACET2(IF2)
          IF (.NOT.SOLPYR) RETURN
       IF (N3.EQ.1) SOLPYR=SOLPYR.AND.FACET2(IF3)
          IF (.NOT.SOLPYR) RETURN
       IF (N4.EQ.1) SOLPYR=SOLPYR.AND.FACET2(IF4)
          IF (.NOT.SOLPYR) RETURN
       IF (N5.EQ.1) SOLPYR=SOLPYR.AND.FACET2(IF5)
          IF (.NOT.SOLPYR) RETURN
C
C  on teste maintenant qu'on ne recouvre pas un autre volume
       ip1=nfc(1,if1)
       ip2=nfc(2,if1)
       ip3=nfc(3,if1)
       ip4=nfc(4,if1)
       ip5=nfc(1,if2)
       if (ip5.eq.ip1.or.ip5.eq.ip2.or.ip5.eq.ip3.or.ip5.eq.ip4)
     *    ip5=nfc(2,if2)
       if (ip5.eq.ip1.or.ip5.eq.ip2.or.ip5.eq.ip3.or.ip5.eq.ip4)
     *    ip5=nfc(3,if2)
*      write (6,*) ' solpyr sommets ',ip1,ip2,ip3,ip4,ip5
*      write (6,*) ' solpyr coordonnees '
*      write (6,*) xyz(1,ip1),xyz(2,ip1),xyz(3,ip1)
*      write (6,*) xyz(1,ip2),xyz(2,ip2),xyz(3,ip2)
*      write (6,*) xyz(1,ip3),xyz(2,ip3),xyz(3,ip3)
*      write (6,*) xyz(1,ip4),xyz(2,ip4),xyz(3,ip4)
*      write (6,*) xyz(1,ip5),xyz(2,ip5),xyz(3,ip5)
       xbar=(xyz(1,ip1)+xyz(1,ip2)+xyz(1,ip3)+xyz(1,ip4)+xyz(1,ip5))/5.
       ybar=(xyz(2,ip1)+xyz(2,ip2)+xyz(2,ip3)+xyz(2,ip4)+xyz(2,ip5))/5.
       zbar=(xyz(3,ip1)+xyz(3,ip2)+xyz(3,ip3)+xyz(3,ip4)+xyz(3,ip5))/5.
*      write (6,*) ' bary ',xbar,ybar,zbar
       xfac=(xyz(1,ip1)+xyz(1,ip2)+xyz(1,ip3)+xyz(1,ip4))/4.
       yfac=(xyz(2,ip1)+xyz(2,ip2)+xyz(2,ip3)+xyz(2,ip4))/4.
       zfac=(xyz(3,ip1)+xyz(3,ip2)+xyz(3,ip3)+xyz(3,ip4))/4.
       xyz(1,nptmax+1)=xfac+0.1*(xbar-xfac)
       xyz(2,nptmax+1)=yfac+0.1*(ybar-yfac)
       xyz(3,nptmax+1)=zfac+0.1*(zbar-zfac)
       solpyr=solpyr.and.vervol(nptmax+1,ip1,ip2,ip3,ip4,ip5,-1,-1,-1)
       if (.not.solpyr) return
       xfac=(xyz(1,ip1)+xyz(1,ip2)+xyz(1,ip5))/3.
       yfac=(xyz(2,ip1)+xyz(2,ip2)+xyz(2,ip5))/3.
       zfac=(xyz(3,ip1)+xyz(3,ip2)+xyz(3,ip5))/3.
       xyz(1,nptmax+1)=xfac+0.1*(xbar-xfac)
       xyz(2,nptmax+1)=yfac+0.1*(ybar-yfac)
       xyz(3,nptmax+1)=zfac+0.1*(zbar-zfac)
       solpyr=solpyr.and.vervol(nptmax+1,ip1,ip2,ip3,ip4,ip5,-1,-1,-1)
       if (.not.solpyr) return
       xfac=(xyz(1,ip2)+xyz(1,ip3)+xyz(1,ip5))/3.
       yfac=(xyz(2,ip2)+xyz(2,ip3)+xyz(2,ip5))/3.
       zfac=(xyz(3,ip2)+xyz(3,ip3)+xyz(3,ip5))/3.
       xyz(1,nptmax+1)=xfac+0.1*(xbar-xfac)
       xyz(2,nptmax+1)=yfac+0.1*(ybar-yfac)
       xyz(3,nptmax+1)=zfac+0.1*(zbar-zfac)
       solpyr=solpyr.and.vervol(nptmax+1,ip1,ip2,ip3,ip4,ip5,-1,-1,-1)
       if (.not.solpyr) return
       xfac=(xyz(1,ip3)+xyz(1,ip4)+xyz(1,ip5))/3.
       yfac=(xyz(2,ip3)+xyz(2,ip4)+xyz(2,ip5))/3.
       zfac=(xyz(3,ip3)+xyz(3,ip4)+xyz(3,ip5))/3.
       xyz(1,nptmax+1)=xfac+0.1*(xbar-xfac)
       xyz(2,nptmax+1)=yfac+0.1*(ybar-yfac)
       xyz(3,nptmax+1)=zfac+0.1*(zbar-zfac)
       solpyr=solpyr.and.vervol(nptmax+1,ip1,ip2,ip3,ip4,ip5,-1,-1,-1)
       if (.not.solpyr) return
       xfac=(xyz(1,ip4)+xyz(1,ip1)+xyz(1,ip5))/3.
       yfac=(xyz(2,ip4)+xyz(2,ip1)+xyz(2,ip5))/3.
       zfac=(xyz(3,ip4)+xyz(3,ip1)+xyz(3,ip5))/3.
       xyz(1,nptmax+1)=xfac+0.1*(xbar-xfac)
       xyz(2,nptmax+1)=yfac+0.1*(ybar-yfac)
       xyz(3,nptmax+1)=zfac+0.1*(zbar-zfac)
       solpyr=solpyr.and.vervol(nptmax+1,ip1,ip2,ip3,ip4,ip5,-1,-1,-1)
       if (.not.solpyr) return
       RETURN
C
120    CONTINUE
C      LE POINT I EST INTERNE AU VOLUME
       IF (IVERB.EQ.1) WRITE(6,1010)I
 1010   FORMAT(' LE POINT ',I5,' EST INTERNE A la pyramide CREEe')
        IF (IVERB.EQ.1) write (6,*) xyz(1,i),xyz(2,i),xyz(3,i)
       SOLPYR=.FALSE.
       RETURN
       END