1 : C DEDANS SOURCE SH236661 13/11/25 21:15:03 7869 2 : C 3 : SUBROUTINE DEDANS(NP1,IPT1,BOOL1) 4 : C----------------------------------------------------------------------C 5 : C Determine si un point est situe a l'interieur C 6 : C - d'un contour oriente ferme constitue de SEG2 (en 2D) C 7 : C - d'une surface enveloppe orientee fermee constituee de TRI3 (en 3D) C 8 : C C 9 : C NP1 = entier, numero du point C 10 : C IPT1 = entier, numero du maillage contour/enveloppe C 11 : C BOOL1 = booleen egal a .TRUE. si NP1 est dans IPT1 C 12 : C----------------------------------------------------------------------C 13 : IMPLICIT INTEGER(I-N) 14 : IMPLICIT REAL*8 (A-H,O-Z) 15 : -INC CCOPTIO 16 : -INC SMELEME 17 : -INC SMCOORD 18 : -INC CCREEL 19 : PARAMETER (ZERO=1D-10) 20 : LOGICAL BOOL1 21 : C 22 : IDIMP1=IDIM+1 23 : BOOL1=.TRUE. 24 : C 25 : C 26 : IF (IDIM.EQ.2) THEN 27 : X1=XCOOR((NP1-1)*IDIMP1+1) 28 : Y1=XCOOR((NP1-1)*IDIMP1+2) 29 : C On calcule la somme des angles des segments du contour vus par 30 : C le point P1 31 : OMEGA=0D0 32 : DO I=1,IPT1.NUM(/2) 33 : C Coordonnees des points A et B du segment I 34 : NRA=IPT1.NUM(1,I) 35 : NRB=IPT1.NUM(2,I) 36 : XA=XCOOR((NRA-1)*IDIMP1+1) 37 : YA=XCOOR((NRA-1)*IDIMP1+2) 38 : XB=XCOOR((NRB-1)*IDIMP1+1) 39 : YB=XCOOR((NRB-1)*IDIMP1+2) 40 : C Calcul de l'angle oriente entre P1A et P1B 41 : C - le produit scalaire donne l'angle (entre 0 et pi) 42 : C - le produit vectoriel donne le signe 43 : RAX=XA-X1 44 : RAY=YA-Y1 45 : RA=SQRT((RAX**2)+(RAY**2)) 46 : RBX=XB-X1 47 : RBY=YB-Y1 48 : RB=SQRT((RBX**2)+(RBY**2)) 49 : PS=(RAX*RBX)+(RAY*RBY) 50 : PV=(RAX*RBY)-(RAY*RBX) 51 : SIG1=SIGN(1D0,PV) 52 : RD=RA*RB 53 : IF (RD.LE.ZERO) THEN 54 : C Si le denominateur est nul, les points sont confondus 55 : C alors le point P1 coincide avec un sommet du segment 56 : C donc il est sur le contour --> on quitte 57 : OMEGA=2D0*XPI 58 : GOTO 1 59 : ENDIF 60 : IF (ABS(PV).LE.ZERO) THEN 61 : C Si le produit vectoriel est nul, le point P1 est situe 62 : C sur la droite AB, l'angle considere comme nul 63 : OMEGA1=0D0 64 : ELSE 65 : COS1=PS/RD 66 : OMEGA1=SIG1*(ACOS(COS1)) 67 : ENDIF 68 : OMEGA=OMEGA+OMEGA1 69 : ENDDO 70 : 1 CONTINUE 71 : OMEGA=OMEGA/(2D0*XPI) 72 : C Si OMEGA est nul alors P1 est dehors 73 : IF (ABS(OMEGA).LT.ZERO) BOOL1=.FALSE. 74 : C 75 : C 76 : ELSEIF (IDIM.EQ.3) THEN 77 : X1=XCOOR((NP1-1)*IDIMP1+1) 78 : Y1=XCOOR((NP1-1)*IDIMP1+2) 79 : Z1=XCOOR((NP1-1)*IDIMP1+3) 80 : C On calcule la somme des angles solides des triangles de 81 : C l'enveloppe vus par le point P1 82 : OMEGA=0D0 83 : DO I=1,IPT1.NUM(/2) 84 : C Coordonnees des points A,B et C du triangle I 85 : NRA=IPT1.NUM(1,I) 86 : NRB=IPT1.NUM(2,I) 87 : NRC=IPT1.NUM(3,I) 88 : XA=XCOOR((NRA-1)*IDIMP1+1) 89 : YA=XCOOR((NRA-1)*IDIMP1+2) 90 : ZA=XCOOR((NRA-1)*IDIMP1+3) 91 : XB=XCOOR((NRB-1)*IDIMP1+1) 92 : YB=XCOOR((NRB-1)*IDIMP1+2) 93 : ZB=XCOOR((NRB-1)*IDIMP1+3) 94 : XC=XCOOR((NRC-1)*IDIMP1+1) 95 : YC=XCOOR((NRC-1)*IDIMP1+2) 96 : ZC=XCOOR((NRC-1)*IDIMP1+3) 97 : C Calcul de l'angle solide oriente du triangle ABC vu depuis 98 : C le point P1 (formule de Van Oosterom et Strackee) 99 : RAX=XA-X1 100 : RAY=YA-Y1 101 : RAZ=ZA-Z1 102 : RA=SQRT((RAX**2)+(RAY**2)+(RAZ**2)) 103 : RBX=XB-X1 104 : RBY=YB-Y1 105 : RBZ=ZB-Z1 106 : RB=SQRT((RBX**2)+(RBY**2)+(RBZ**2)) 107 : RCX=XC-X1 108 : RCY=YC-Y1 109 : RCZ=ZC-Z1 110 : RC=SQRT((RCX**2)+(RCY**2)+(RCZ**2)) 111 : RD=(RA*RB*RC) + (((RAX*RBX)+(RAY*RBY)+(RAZ*RBZ))*RC) 112 : & + (((RAX*RCX)+(RAY*RCY)+(RAZ*RCZ))*RB) 113 : & + (((RBX*RCX)+(RBY*RCY)+(RBZ*RCZ))*RA) 114 : IF (ABS(RD).LE.ZERO) THEN 115 : C Si le denominateur est nul, alors le point P1 coincide 116 : C avec un sommet du triangle ABC, donc il est sur 117 : C l'enveloppe --> on quitte 118 : OMEGA=4D0*XPI 119 : GOTO 2 120 : ENDIF 121 : RN= (RAX*((RBY*RCZ)-(RBZ*RCY))) 122 : & -(RBX*((RAY*RCZ)-(RAZ*RCY))) 123 : & +(RCX*((RAY*RBZ)-(RAZ*RBY))) 124 : IF (ABS(RN).LE.ZERO) THEN 125 : C Si le numerateur est nul, le point P1 est dans le plan 126 : C du triangle ABC, donc l'angle solide est considere comme 127 : C nul 128 : OMEGA1=0D0 129 : ELSE 130 : AT=ATAN2(RN,RD) 131 : OMEGA1=2D0*AT 132 : ENDIF 133 : OMEGA=OMEGA+OMEGA1 134 : ENDDO 135 : 2 CONTINUE 136 : OMEGA=OMEGA/(4D0*XPI) 137 : C Si OMEGA est nul alors P1 est dehors 138 : IF (ABS(OMEGA).LT.ZERO) BOOL1=.FALSE. 139 : ENDIF 140 : C 141 : C Fin du programme 142 : 999 RETURN 143 : END 144 : 145 : 146 :
© Cast3M 2003 - All rights reserved.
Disclaimer