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