Télécharger dedans.eso

Retour à la liste

Numérotation des lignes :

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales