Sources Esope | Cast3M

soltet
C SOLTET    SOURCE    PV        17/09/12    21:15:04     9542           
C---------------------------------------------------------------------|
C                                                                     |
       LOGICAL FUNCTION SOLTET(IF1,IF2,IF3,IF4,IPIN)
C                                                                     |
C      CETTE FONCTION LOGIQUE TESTE SI LE TETRAEDRE DECRIT PAR LES    |
C      FACETTES IF1..IF4 EST VIDE DE POINT.                           |
C      LES ENTIERS N1..N4 INDIQUENT L'ORIENTATION DES FACETTES        |
C      SOLTET EST VRAI SI LE TETRAEDRE EST VIDE (DONC VALIDE)         |
C      SOLTET EST FAUX SI LE TETRAEDRE 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 VAL,FACET2,FACET3,vervol
C
      dimension ift(4)
      real*8 xo(3),xa(3),xb(3),xc(3),eps
      IPIN=0
      ift(1)=if1
      ift(2)=if2
      ift(3)=if3
      ift(4)=if4
       ip1=nfc(1,if1)
       ip2=nfc(2,if1)
       ip3=nfc(3,if1)
       ip4=nfc(1,if2)
       if (ip4.eq.ip1.or.ip4.eq.ip2.or.ip4.eq.ip3) ip4=nfc(2,if2)
       if (ip4.eq.ip1.or.ip4.eq.ip2.or.ip4.eq.ip3) ip4=nfc(3,if2)
C
C   CHERCHER SI LES FACETTES SONT NOUVELLES OU ANCIENNES
        N1=-1
        N2=-1
        N3=-1
        N4=-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
*  encadrer le tetraedre
       xmin=1e30
       xmax=-1e30
       ymin=1e30
       ymax=-1e30
       zmin=1e30
       zmax=-1e30
       xmin=min(xyz(1,ip1),xyz(1,ip2),xyz(1,ip3),xyz(1,ip4))
       xmax=max(xyz(1,ip1),xyz(1,ip2),xyz(1,ip3),xyz(1,ip4))
       xmin=min(xyz(2,ip1),xyz(2,ip2),xyz(2,ip3),xyz(2,ip4))
       xmax=max(xyz(2,ip1),xyz(2,ip2),xyz(2,ip3),xyz(2,ip4))
       xmin=min(xyz(3,ip1),xyz(3,ip2),xyz(3,ip3),xyz(3,ip4))
       xmax=max(xyz(3,ip1),xyz(3,ip2),xyz(3,ip3),xyz(3,ip4))
       xd=(xmax-xmax)/2
       yd=(ymax-ymax)/2
       zd=(zmax-zmax)/2
       td=max(xd,yd,zd)
       xmin=xmin-td
       xmax=xmax+td
       ymin=ymin-td
       ymax=ymax+td
       zmin=zmin-td
       zmax=zmax+td
       iteste=0
       V1=N1*VOL(1,NFC(1,IF1),NFC(2,IF1),NFC(3,IF1))
       V2=N2*VOL(1,NFC(1,IF2),NFC(2,IF2),NFC(3,IF2))
       V3=N3*VOL(1,NFC(1,IF3),NFC(2,IF3),NFC(3,IF3))
       V4=N4*VOL(1,NFC(1,IF4),NFC(2,IF4),NFC(3,IF4))
       VV=V1+V2+V3+V4
       tm=(xyz(4,ip1)+xyz(4,ip2)+xyz(4,ip3)+xyz(4,ip4))/4
       tv=tm**3*1e-6
       IF (VV.LE.tv) then
        if(IVERB.EQ.1)  write (6,*) ' tetraedre volume negatif '
        GOTO 120
       endif
*      DO 100 I=1,NPTMAX
       DO 100 I=1,-1
          IF (NPF(1,I).EQ.0) GOTO 100
          IF (I.EQ.IP1) GOTO 100
          IF (I.EQ.IP2) GOTO 100
          IF (I.EQ.IP3) GOTO 100
          IF (I.EQ.IP4) GOTO 100
       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
C
       V1=N1*VOL(I,NFC(1,IF1),NFC(2,IF1),NFC(3,IF1))
       V2=N2*VOL(I,NFC(1,IF2),NFC(2,IF2),NFC(3,IF2))
       V3=N3*VOL(I,NFC(1,IF3),NFC(2,IF3),NFC(3,IF3))
       V4=N4*VOL(I,NFC(1,IF4),NFC(2,IF4),NFC(3,IF4))
       IF (V1.LE.-volcri*VV) GOTO 100
       IF (V2.LE.-volcri*VV) GOTO 100
       IF (V3.LE.-volcri*VV) GOTO 100
       IF (V4.LE.-volcri*VV) GOTO 100
       IF (IVERB.EQ.1) WRITE(6,1010)I
 1010   FORMAT(' LE POINT ',I5,' EST INTERNE AU tetraedre CREE |')
       IF (IVERB.EQ.1) write (6,*) ' vv,v1,v2,v3,v4 ',vv,v1,v2,v3,v4
       ipin=i
       GOTO 120
C
C
C
100    CONTINUE
C
C      IL N'EXISTE PAS DE POINTS INTERNES AU VOLUME
       SOLTET=.TRUE.
C  RAJOUT PV TEST ON NE CREE PAS D'ARETE A PLUS DE DEUX FACETTES
*      IF (N1.EQ.1) SOLTET=SOLTET.AND.FACET2(IF1)
*        IF (.NOT.SOLTET) RETURN
*      IF (N2.EQ.1) SOLTET=SOLTET.AND.FACET2(IF2)
*        IF (.NOT.SOLTET) RETURN
*      IF (N3.EQ.1) SOLTET=SOLTET.AND.FACET2(IF3)
*        IF (.NOT.SOLTET) RETURN
*      IF (N4.EQ.1) SOLTET=SOLTET.AND.FACET2(IF4)
*        IF (.NOT.SOLTET) RETURN
C
C  on teste maintenant qu'on ne recouvre pas un autre volume
*      write (6,*) ' soltet sommets ',ip1,ip2,ip3,ip4
*      write (6,*) ' soltet 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)
       xbar=(xyz(1,ip1)+xyz(1,ip2)+xyz(1,ip3)+xyz(1,ip4))/4.
       ybar=(xyz(2,ip1)+xyz(2,ip2)+xyz(2,ip3)+xyz(2,ip4))/4.
       zbar=(xyz(3,ip1)+xyz(3,ip2)+xyz(3,ip3)+xyz(3,ip4))/4.
*      write (6,*) ' bary ',xbar,ybar,zbar
       xfac=(xyz(1,ip1)+xyz(1,ip2)+xyz(1,ip3))/3.
       yfac=(xyz(2,ip1)+xyz(2,ip2)+xyz(2,ip3))/3
       zfac=(xyz(3,ip1)+xyz(3,ip2)+xyz(3,ip3))/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)
       soltet=soltet.and.vervol(nptmax+1,ip1,ip2,ip3,ip4,-1,-1,-1,-1)
       if (.not.soltet) return
       xfac=(xyz(1,ip1)+xyz(1,ip2)+xyz(1,ip4))/3.
       yfac=(xyz(2,ip1)+xyz(2,ip2)+xyz(2,ip4))/3
       zfac=(xyz(3,ip1)+xyz(3,ip2)+xyz(3,ip4))/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)
       soltet=soltet.and.vervol(nptmax+1,ip1,ip2,ip3,ip4,-1,-1,-1,-1)
       if (.not.soltet) return
       xfac=(xyz(1,ip2)+xyz(1,ip3)+xyz(1,ip4))/3.
       yfac=(xyz(2,ip2)+xyz(2,ip3)+xyz(2,ip4))/3
       zfac=(xyz(3,ip2)+xyz(3,ip3)+xyz(3,ip4))/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)
       soltet=soltet.and.vervol(nptmax+1,ip1,ip2,ip3,ip4,-1,-1,-1,-1)
       if (.not.soltet) return
       xfac=(xyz(1,ip3)+xyz(1,ip1)+xyz(1,ip4))/3.
       yfac=(xyz(2,ip3)+xyz(2,ip1)+xyz(2,ip4))/3
       zfac=(xyz(3,ip3)+xyz(3,ip1)+xyz(3,ip4))/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)
       soltet=soltet.and.vervol(nptmax+1,ip1,ip2,ip3,ip4,-1,-1,-1,-1)
       if (.not.soltet) return
*  test sur les facettes touchees
*      do 200 i=1,4
*       if (i.eq.1) ip=ip1
*       if (i.eq.2) ip=ip2
*       if (i.eq.3) ip=ip3
*       if (i.eq.4) ip=ip4
*       do 210 j=1,40
*        if=npf(j,ip)
*        if (if.eq.0) goto 200
*        iprob=0
*        do 220 k=1,3
*         jp=nfc(k,if)
*         if (npf(4,jp).ne.0.and.npf(5,jp).eq.0) iprob=iprob+1
*220     continue
*        if (iprob.eq.3) soltet=.false.
*210    continue
*200   continue
*      if (.not.soltet) write (6,*) ' soltet nouveau test echoue'
*  test sur qualite
*      if (qualt(ip1,ip2,ip3,ip4).lt.5e-2) then
*        write (6,*) ' soltet qualite insuffisante '
*        soltet=.false.
*      endif
       RETURN
C
120    CONTINUE
C      LE POINT I EST INTERNE AU VOLUME
        IF (IVERB.EQ.1) write (6,*) ' facettes ',if1,if2,if3,if4
       SOLTET=.FALSE.
       RETURN
       END