Sources Esope | Cast3M

dedans
C DEDANS    SOURCE    CB215821  23/01/25    21:15:11     11573          
C
      SUBROUTINE DEDANS(NP1,IPT1,XTOL,BOOL1)
C----------------------------------------------------------------------C
C            Determine si un point est situe a l'interieur             C
C - d'un contour oriente ferme constitue de SEG2 (en 2D)               C
C - d'une surface enveloppe orientee fermee constituee de TRI3 (en 3D) C
C                                                                      C
C     NP1   = entier, numero du point                                  C
C     IPT1  = entier, numero du maillage contour/enveloppe             C
C     XTOL  = flottant, tolerance pour le test de nullite de l'angle   C
C             solide                                                   C
C     BOOL1 = booleen, egal a .TRUE. si NP1 est dans IPT1              C
C----------------------------------------------------------------------C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC SMELEME
-INC SMCOORD
-INC CCREEL
C
C     Zero utilise pour tester la colinearite/coplaneite du point avec
C     les elements du contour/enveloppe
      PARAMETER (ZERO=1D-10)
      LOGICAL BOOL1
C
      IDIMP1=IDIM+1
      BOOL1=.TRUE.
C
C
      IF (IDIM.EQ.2) THEN
         X1=XCOOR((NP1-1)*IDIMP1+1)
         Y1=XCOOR((NP1-1)*IDIMP1+2)
C        On calcule la somme des angles des segments du contour vus par
C        le point P1
         OMEGA=0D0
         DO I=1,IPT1.NUM(/2)
C           Coordonnees des points A et B du segment I
            NRA=IPT1.NUM(1,I)
            NRB=IPT1.NUM(2,I)
            XA=XCOOR((NRA-1)*IDIMP1+1)
            YA=XCOOR((NRA-1)*IDIMP1+2)
            XB=XCOOR((NRB-1)*IDIMP1+1)
            YB=XCOOR((NRB-1)*IDIMP1+2)
C           Calcul de l'angle oriente entre P1A et P1B
C           - le produit scalaire donne l'angle (entre 0 et pi)
C           - le produit vectoriel donne le signe
            RAX=XA-X1
            RAY=YA-Y1
            RA=SQRT((RAX*RAX)+(RAY*RAY))
            RBX=XB-X1
            RBY=YB-Y1
            RB=SQRT((RBX*RBX)+(RBY*RBY))
            PS=(RAX*RBX)+(RAY*RBY)
            PV=(RAX*RBY)-(RAY*RBX)
            SIG1=SIGN(1D0,PV)
            RD=RA*RB
            IF ((RA.LE.ZERO).OR.(RB.LE.ZERO).OR.(RD.LE.ZERO)) THEN
C              Si le denominateur est nul, les points sont confondus
C              alors le point P1 coincide avec un sommet du segment
C              donc il est sur le contour --> on quitte
               OMEGA=2D0*XPI
               GOTO 1
            ENDIF
            IF (ABS(PV).LE.ZERO) THEN
C              Si le produit vectoriel est nul, le point P1 est situe
C              sur la droite AB, l'angle est considere comme nul
               OMEGA1=0D0
            ELSE
               COS1=PS/RD
               IF (COS1.GT. 1.D0) COS1=1D0
               IF (COS1.LT.-1.D0) COS1=-1D0
               OMEGA1=SIG1*(ACOS(COS1))
            ENDIF
            OMEGA=OMEGA+OMEGA1
         ENDDO
1        CONTINUE
         OMEGA=OMEGA/(2D0*XPI)
C        Si OMEGA est nul alors P1 est dehors
         IF (ABS(OMEGA).LT.XTOL) BOOL1=.FALSE.
C
C
      ELSEIF (IDIM.EQ.3) THEN
         X1=XCOOR((NP1-1)*IDIMP1+1)
         Y1=XCOOR((NP1-1)*IDIMP1+2)
         Z1=XCOOR((NP1-1)*IDIMP1+3)
C        On calcule la somme des angles solides des triangles de
C        l'enveloppe vus par le point P1
         OMEGA=0D0
         DO I=1,IPT1.NUM(/2)
C           Coordonnees des points A,B et C du triangle I
            NRA=IPT1.NUM(1,I)
            NRB=IPT1.NUM(2,I)
            NRC=IPT1.NUM(3,I)
            XA=XCOOR((NRA-1)*IDIMP1+1)
            YA=XCOOR((NRA-1)*IDIMP1+2)
            ZA=XCOOR((NRA-1)*IDIMP1+3)
            XB=XCOOR((NRB-1)*IDIMP1+1)
            YB=XCOOR((NRB-1)*IDIMP1+2)
            ZB=XCOOR((NRB-1)*IDIMP1+3)
            XC=XCOOR((NRC-1)*IDIMP1+1)
            YC=XCOOR((NRC-1)*IDIMP1+2)
            ZC=XCOOR((NRC-1)*IDIMP1+3)
C           Calcul de l'angle solide oriente du triangle ABC vu depuis
C           le point P1 (formule de Van Oosterom et Strackee)
            RAX=XA-X1
            RAY=YA-Y1
            RAZ=ZA-Z1
            RA=SQRT((RAX*RAX)+(RAY*RAY)+(RAZ*RAZ))
            RBX=XB-X1
            RBY=YB-Y1
            RBZ=ZB-Z1
            RB=SQRT((RBX*RBX)+(RBY*RBY)+(RBZ*RBZ))
            RCX=XC-X1
            RCY=YC-Y1
            RCZ=ZC-Z1
            RC=SQRT((RCX*RCX)+(RCY*RCY)+(RCZ*RCZ))
            RD=(RA*RB*RC) + (((RAX*RBX)+(RAY*RBY)+(RAZ*RBZ))*RC)
     &                    + (((RAX*RCX)+(RAY*RCY)+(RAZ*RCZ))*RB)
     &                    + (((RBX*RCX)+(RBY*RCY)+(RBZ*RCZ))*RA)
            IF ((RA.LE.ZERO).OR.(RB.LE.ZERO).OR.(RC.LE.ZERO).OR.
     &          (ABS(RD).LE.ZERO)) THEN
C              Si le denominateur est nul, alors le point P1 coincide
C              avec un sommet du triangle ABC, donc il est sur
C              l'enveloppe --> on quitte
               OMEGA=4D0*XPI
               GOTO 2
            ENDIF
            RN= (RAX*((RBY*RCZ)-(RBZ*RCY)))
     &         -(RBX*((RAY*RCZ)-(RAZ*RCY)))
     &         +(RCX*((RAY*RBZ)-(RAZ*RBY)))
            IF (ABS(RN).LE.ZERO) THEN
C              Si le numerateur est nul, le point P1 est dans le plan
C              du triangle ABC, donc l'angle solide est considere comme
C              nul
               OMEGA1=0D0
            ELSE
               AT=ATAN2(RN,RD)
               OMEGA1=2D0*AT
            ENDIF
            OMEGA=OMEGA+OMEGA1
         ENDDO
2        CONTINUE
         OMEGA=OMEGA/(4D0*XPI)
C        Si OMEGA est nul alors P1 est dehors
         IF (ABS(OMEGA).LT.XTOL) BOOL1=.FALSE.
      ENDIF
C
C     Fin du programme
999   RETURN
      END