Télécharger aspert.eso

Retour à la liste

Numérotation des lignes :

  1. C ASPERT SOURCE FD218221 20/06/17 21:15:01 10630
  2. SUBROUTINE ASPERT(VAL,NUM,NN,NE,ITYP)
  3. C-----------------------------------------------------------------------
  4. C Calcul de l'aspect ratio des elements d'un maillage elementaire
  5. C Fonctionne pour les elements TRI3,TET4
  6. C
  7. C L'aspect ratio est definit ainsi :
  8. C coef * taille de la plus grande arete / taille de la plus petite hauteur
  9. C
  10. C avec : coef = SQRT(3)/2 pour les elements TRI3
  11. C coef = SQRT(2/3) pour les elements TET4
  12. C
  13. C Ainsi, pour un triangle/tetraedre regulier ASPE = 1
  14. C les bornes : 1 < ASPE < +inf
  15. C
  16. C Entrees :
  17. C NUM(NN,NE) : tableau des numeros des elements du maillage
  18. C NN : nombre de noeuds par element
  19. C NE : nombre d'elements
  20. C ITYP : numero indiquant le type d'element
  21. C Sorties :
  22. C VAL(1,NE) : tableau des valeurs d'aspect ration pour les elements
  23. C du maillage (il s'agit du VELCHE d'un champ par elements)
  24. C-----------------------------------------------------------------------
  25. C
  26. IMPLICIT INTEGER(I-N)
  27. IMPLICIT REAL*8(A-H,O-Z)
  28. -INC PPARAM
  29. -INC CCOPTIO
  30. -INC SMCOORD
  31. C
  32. DIMENSION VAL(1,*),NUM(NN,NE)
  33. C Quelques nombres constants
  34. R23=SQRT(2.D0/3.D0)
  35. R3S2=0.5D0*SQRT(3.D0)
  36. C La dimension de l'espace + 1
  37. IDIMP1=IDIM+1
  38. C Erreur si type d'element non prevu (normalement impossible de passer ici !)
  39. IF ((ITYP.NE.4).AND.(ITYP.NE.23)) THEN
  40. CALL ERREUR(16)
  41. RETURN
  42. ENDIF
  43. C Boucle sur les elements du maillage
  44. DO I=1,NE
  45. C Cas des elements TRI3
  46. IF (ITYP.EQ.4) THEN
  47. C Numeros des 3 sommets
  48. IA=NUM(1,I)
  49. IB=NUM(2,I)
  50. IC=NUM(3,I)
  51. C Coordonnees des 3 sommets
  52. XA=XCOOR((IA-1)*IDIMP1+1)
  53. YA=XCOOR((IA-1)*IDIMP1+2)
  54. XB=XCOOR((IB-1)*IDIMP1+1)
  55. YB=XCOOR((IB-1)*IDIMP1+2)
  56. XC=XCOOR((IC-1)*IDIMP1+1)
  57. YC=XCOOR((IC-1)*IDIMP1+2)
  58. IF (IDIM.EQ.3) THEN
  59. ZA=XCOOR((IA-1)*IDIMP1+3)
  60. ZB=XCOOR((IB-1)*IDIMP1+3)
  61. ZC=XCOOR((IC-1)*IDIMP1+3)
  62. ENDIF
  63. C Tailles des 3 aretes
  64. XBA=XB-XA
  65. YBA=YB-YA
  66. IF (IDIM.EQ.2) THEN
  67. AB=SQRT(XBA*XBA+YBA*YBA)
  68. ELSEIF (IDIM.EQ.3) THEN
  69. ZBA=ZB-ZA
  70. AB=SQRT(XBA*XBA+YBA*YBA+ZBA*ZBA)
  71. ENDIF
  72. XCA=XC-XA
  73. YCA=YC-YA
  74. IF (IDIM.EQ.2) THEN
  75. AC=SQRT(XCA*XCA+YCA*YCA)
  76. ELSEIF (IDIM.EQ.3) THEN
  77. ZCA=ZC-ZA
  78. AC=SQRT(XCA*XCA+YCA*YCA+ZCA*ZCA)
  79. ENDIF
  80. XCB=XC-XB
  81. YCB=YC-YB
  82. IF (IDIM.EQ.2) THEN
  83. BC=SQRT(XCB*XCB+YCB*YCB)
  84. ELSEIF (IDIM.EQ.3) THEN
  85. ZCB=ZC-ZB
  86. BC=SQRT(XCB*XCB+YCB*YCB+ZCB*ZCB)
  87. ENDIF
  88. C AB.AC
  89. SABAC=XBA*XCA+YBA*YCA
  90. IF (IDIM.EQ.3) SABAC=SABAC+ZBA*ZCA
  91. C AB.BC
  92. SABBC=XBA*XCB+YBA*YCB
  93. IF (IDIM.EQ.3) SABBC=SABBC+ZBA*ZCB
  94. C --> hauteur sous A
  95. BH=-SABBC/BC
  96. HA=SQRT(AB*AB-BH*BH)
  97. C --> hauteur sous C
  98. AH=SABAC/AB
  99. HC=SQRT(AC*AC-AH*AH)
  100. C --> hauteur sous B
  101. AH=SABAC/AC
  102. HB=SQRT(AB*AB-AH*AH)
  103. C Valeur de l'aspect ratio
  104. ASPE=R3S2*MAX(AB,AC,BC)/MIN(HA,HB,HC)
  105. C Cas des elements TET4
  106. ELSEIF (ITYP.EQ.23) THEN
  107. C Numeros des 4 sommets
  108. IA=NUM(1,I)
  109. IB=NUM(2,I)
  110. IC=NUM(3,I)
  111. ID=NUM(4,I)
  112. C Coordonnees des 4 sommets
  113. XA=XCOOR((IA-1)*IDIMP1+1)
  114. YA=XCOOR((IA-1)*IDIMP1+2)
  115. ZA=XCOOR((IA-1)*IDIMP1+3)
  116. XB=XCOOR((IB-1)*IDIMP1+1)
  117. YB=XCOOR((IB-1)*IDIMP1+2)
  118. ZB=XCOOR((IB-1)*IDIMP1+3)
  119. XC=XCOOR((IC-1)*IDIMP1+1)
  120. YC=XCOOR((IC-1)*IDIMP1+2)
  121. ZC=XCOOR((IC-1)*IDIMP1+3)
  122. XD=XCOOR((ID-1)*IDIMP1+1)
  123. YD=XCOOR((ID-1)*IDIMP1+2)
  124. ZD=XCOOR((ID-1)*IDIMP1+3)
  125. C Tailles des 6 aretes
  126. XBA=XB-XA
  127. YBA=YB-YA
  128. ZBA=ZB-ZA
  129. AB=SQRT(XBA*XBA+YBA*YBA+ZBA*ZBA)
  130. XCA=XC-XA
  131. YCA=YC-YA
  132. ZCA=ZC-ZA
  133. AC=SQRT(XCA*XCA+YCA*YCA+ZCA*ZCA)
  134. XDA=XD-XA
  135. YDA=YD-YA
  136. ZDA=ZD-ZA
  137. AD=SQRT(XDA*XDA+YDA*YDA+ZDA*ZDA)
  138. XCB=XC-XB
  139. YCB=YC-YB
  140. ZCB=ZC-ZB
  141. BC=SQRT(XCB*XCB+YCB*YCB+ZCB*ZCB)
  142. XDB=XD-XB
  143. YDB=YD-YB
  144. ZDB=ZD-ZB
  145. BD=SQRT(XDB*XDB+YDB*YDB+ZDB*ZDB)
  146. XDC=XD-XC
  147. YDC=YD-YC
  148. ZDC=ZD-ZC
  149. CD=SQRT(XDC*XDC+YDC*YDC+ZDC*ZDC)
  150. C Vecteurs normaux unitaires aux 4 plans
  151. C --> AB^AC
  152. XABAC=YBA*ZCA-ZBA*YCA
  153. YABAC=ZBA*XCA-XBA*ZCA
  154. ZABAC=XBA*YCA-YBA*XCA
  155. C --> normale au plan ABC
  156. SABC=SQRT(XABAC*XABAC+YABAC*YABAC+ZABAC*ZABAC)
  157. XNABC=XABAC/SABC
  158. YNABC=YABAC/SABC
  159. ZNABC=ZABAC/SABC
  160. C --> AB^AD
  161. XABAD=YBA*ZDA-ZBA*YDA
  162. YABAD=ZBA*XDA-XBA*ZDA
  163. ZABAD=XBA*YDA-YBA*XDA
  164. C --> normale au plan ABD
  165. SABD=SQRT(XABAD*XABAD+YABAD*YABAD+ZABAD*ZABAD)
  166. XNABD=XABAD/SABD
  167. YNABD=YABAD/SABD
  168. ZNABD=ZABAD/SABD
  169. C --> AC^AD
  170. XACAD=YCA*ZDA-ZCA*YDA
  171. YACAD=ZCA*XDA-XCA*ZDA
  172. ZACAD=XCA*YDA-YCA*XDA
  173. C --> normale au plan ACD
  174. SACD=SQRT(XACAD*XACAD+YACAD*YACAD+ZACAD*ZACAD)
  175. XNACD=XACAD/SACD
  176. YNACD=YACAD/SACD
  177. ZNACD=ZACAD/SACD
  178. C --> BC^BD
  179. XBCBD=YCB*ZDB-ZCB*YDB
  180. YBCBD=ZCB*XDB-XCB*ZDB
  181. ZBCBD=XCB*YDB-YCB*XDB
  182. C --> normale au plan BCD
  183. SBCD=SQRT(XBCBD*XBCBD+YBCBD*YBCBD+ZBCBD*ZBCBD)
  184. XNBCD=XBCBD/SBCD
  185. YNBCD=YBCBD/SBCD
  186. ZNBCD=ZBCBD/SBCD
  187. C Hauteurs sous les 4 sommets
  188. QBCD=-(XNBCD*XB+YNBCD*YB+ZNBCD*ZB)
  189. HA=ABS(XNBCD*XA+YNBCD*YA+ZNBCD*ZA+QBCD)
  190. QACD=-(XNACD*XA+YNACD*YA+ZNACD*ZA)
  191. HB=ABS(XNACD*XB+YNACD*YB+ZNACD*ZB+QACD)
  192. QABD=-(XNABD*XA+YNABD*YA+ZNABD*ZA)
  193. HC=ABS(XNABD*XC+YNABD*YC+ZNABD*ZC+QABD)
  194. QABC=-(XNABC*XA+YNABC*YA+ZNABC*ZA)
  195. HD=ABS(XNABC*XD+YNABC*YD+ZNABC*ZD+QABC)
  196. C Valeur de l'aspect ratio
  197. ASPE=R23*MAX(AB,AC,AD,BC,BD,CD)/MIN(HA,HB,HC,HD)
  198. ENDIF
  199. C Ecriture dans le tableau
  200. VAL(1,I)=ASPE
  201. ENDDO
  202. C
  203. RETURN
  204. END
  205.  
  206.  

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