Télécharger pb443.eso

Retour à la liste

Numérotation des lignes :

  1. C PB443 SOURCE CHAT 05/01/13 02:10:40 5004
  2. SUBROUTINE PB443(X,Y,Z,PG,FN,GR,FM,GM,ND,NP,MP,NGG,NPG,NOM2)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. C************************************************************************
  6. C
  7. C CALCULE LES FONCTIONS DE FORME D'UN : Iso-Q2 (iso P1/P0 nc) TE10
  8. C
  9. C
  10. C
  11. C ^ zeta
  12. C |
  13. C |10
  14. C |
  15. C |
  16. C | 9
  17. C | /\ /^ eta
  18. C | / \ /
  19. C |/____\ /
  20. C |7 8 /
  21. C | /
  22. C | 5
  23. C | / \
  24. C | / \
  25. C | 6_____\4
  26. C | / \ /\
  27. C | / \ / \
  28. C |/_____\/____\ ____________________>ksi
  29. C 1 2 3
  30. C
  31. C
  32. C************************************************************************
  33. CHARACTER*4 NOM2
  34. REAL*8 AL,BE
  35. PARAMETER (AL=.5854101966249684D0)
  36. PARAMETER (BE=.1381966011250105D0)
  37. REAL*8 X(NPG),Y(NPG),Z(NPG)
  38. REAL*8 FN(NP,NPG),GR(ND,NP,NPG),PG(NPG)
  39.  
  40. REAL*8 FM(MP,NPG),GM(ND,MP,NPG)
  41. INTEGER CONEK(4,8),L,K,I,N,N1,N2,N3,N4,IGAU,IFN
  42. INTEGER ND,MP,NGG,NP,NPG
  43. REAL*8 KXYZ(3,10),AA(4)
  44. REAL*8 T1X,T1Y,T1Z,C1,P1C
  45. REAL*8 T2X,T2Y,T2Z,C2,P2C
  46. REAL*8 T3X,T3Y,T3Z,C3,P3C
  47. REAL*8 T4X,T4Y,T4Z,C4,P4C
  48. C
  49. C Donnees de la numerotation locale des 4 noeuds de chaque tetraedre
  50. C
  51. DATA CONEK/1,2,6,7,
  52. & 2,8,6,7,
  53. & 2,3,4,8,
  54. & 2,4,6,8,
  55. & 4,5,6,9,
  56. & 8,9,4,6,
  57. & 6,7,8,9,
  58. & 7,8,9,10/
  59. C
  60. C Donnees des coordonnees des points du tetraedre TE10
  61. C X Y Z
  62. DATA KXYZ/0.0, 0.0, 0.0,
  63. & 0.5, 0.0, 0.0,
  64. & 1.0, 0.0, 0.0,
  65. & 0.5, 0.5, 0.0,
  66. & 0.0, 1.0, 0.0,
  67. & 0.0, 0.5, 0.0,
  68. & 0.0, 0.0, 0.5,
  69. & 0.5, 0.0, 0.5,
  70. & 0.0, 0.5, 0.5,
  71. & 0.0, 0.0, 1.0/
  72.  
  73. C
  74. C On initialise les fonctions tests et les gradients a 0
  75. C
  76. DO 10 L=1,NPG
  77. C Poids de Gauss
  78. PG(L) = 0.1875D0
  79. DO 22 K=1,MP
  80. FM(K,L)=0.0D0
  81. DO 22 I=1,ND
  82. GM(I,K,L)=0.0D0
  83. 22 CONTINUE
  84.  
  85. DO 20 K=1,NP
  86. FN(K,L)=0.0D0
  87. DO 20 I=1,ND
  88. GR(I,K,L)=0.0D0
  89. 20 CONTINUE
  90.  
  91. 10 CONTINUE
  92. C
  93. C On traite chaque tetraedre elementaire
  94. C
  95. DO 40 K=1,8
  96. C
  97. C On travaille sur chacun des tetraedres composant le macro
  98. C
  99. N1=CONEK(1,K)
  100. N2=CONEK(2,K)
  101. N3=CONEK(3,K)
  102. N4=CONEK(4,K)
  103. C Calcul des equations de plans oppose a chaque point
  104. C
  105. C Plan 1
  106. T1X=(KXYZ(2,N3)-KXYZ(2,N2))*(KXYZ(3,N4)-KXYZ(3,N2))
  107. & -(KXYZ(3,N3)-KXYZ(3,N2))*(KXYZ(2,N4)-KXYZ(2,N2))
  108. T1Y=(KXYZ(3,N3)-KXYZ(3,N2))*(KXYZ(1,N4)-KXYZ(1,N2))
  109. & -(KXYZ(1,N3)-KXYZ(1,N2))*(KXYZ(3,N4)-KXYZ(3,N2))
  110. T1Z=(KXYZ(1,N3)-KXYZ(1,N2))*(KXYZ(2,N4)-KXYZ(2,N2))
  111. & -(KXYZ(2,N3)-KXYZ(2,N2))*(KXYZ(1,N4)-KXYZ(1,N2))
  112. C1=-(T1X*KXYZ(1,N2)+T1Y*KXYZ(2,N2)+T1Z*KXYZ(3,N2))
  113. P1C=T1X*KXYZ(1,N1)+T1Y*KXYZ(2,N1)+T1Z*KXYZ(3,N1)+C1
  114. T1X=T1X/P1C
  115. T1Y=T1Y/P1C
  116. T1Z=T1Z/P1C
  117. C1=C1/P1C
  118. C Plan 2
  119. T2X=(KXYZ(2,N4)-KXYZ(2,N3))*(KXYZ(3,N1)-KXYZ(3,N3))
  120. & -(KXYZ(3,N4)-KXYZ(3,N3))*(KXYZ(2,N1)-KXYZ(2,N3))
  121. T2Y=(KXYZ(3,N4)-KXYZ(3,N3))*(KXYZ(1,N1)-KXYZ(1,N3))
  122. & -(KXYZ(1,N4)-KXYZ(1,N3))*(KXYZ(3,N1)-KXYZ(3,N3))
  123. T2Z=(KXYZ(1,N4)-KXYZ(1,N3))*(KXYZ(2,N1)-KXYZ(2,N3))
  124. & -(KXYZ(2,N4)-KXYZ(2,N3))*(KXYZ(1,N1)-KXYZ(1,N3))
  125. C2=-(T2X*KXYZ(1,N3)+T2Y*KXYZ(2,N3)+T2Z*KXYZ(3,N3))
  126. P2C=T2X*KXYZ(1,N2)+T2Y*KXYZ(2,N2)+T2Z*KXYZ(3,N2)+C2
  127. T2X=T2X/P2C
  128. T2Y=T2Y/P2C
  129. T2Z=T2Z/P2C
  130. C2=C2/P2C
  131. C Plan 3
  132. T3X=(KXYZ(2,N1)-KXYZ(2,N4))*(KXYZ(3,N2)-KXYZ(3,N4))
  133. & -(KXYZ(3,N1)-KXYZ(3,N4))*(KXYZ(2,N2)-KXYZ(2,N4))
  134. T3Y=(KXYZ(3,N1)-KXYZ(3,N4))*(KXYZ(1,N2)-KXYZ(1,N4))
  135. & -(KXYZ(1,N1)-KXYZ(1,N4))*(KXYZ(3,N2)-KXYZ(3,N4))
  136. T3Z=(KXYZ(1,N1)-KXYZ(1,N4))*(KXYZ(2,N2)-KXYZ(2,N4))
  137. & -(KXYZ(2,N1)-KXYZ(2,N4))*(KXYZ(1,N2)-KXYZ(1,N4))
  138. C3=-(T3X*KXYZ(1,N4)+T3Y*KXYZ(2,N4)+T3Z*KXYZ(3,N4))
  139. P3C=T3X*KXYZ(1,N3)+T3Y*KXYZ(2,N3)+T3Z*KXYZ(3,N3)+C3
  140. T3X=T3X/P3C
  141. T3Y=T3Y/P3C
  142. T3Z=T3Z/P3C
  143. C3=C3/P3C
  144. C Plan 4
  145. T4X=(KXYZ(2,N2)-KXYZ(2,N1))*(KXYZ(3,N3)-KXYZ(3,N1))
  146. & -(KXYZ(3,N2)-KXYZ(3,N1))*(KXYZ(2,N3)-KXYZ(2,N1))
  147. T4Y=(KXYZ(3,N2)-KXYZ(3,N1))*(KXYZ(1,N3)-KXYZ(1,N1))
  148. & -(KXYZ(1,N2)-KXYZ(1,N1))*(KXYZ(3,N3)-KXYZ(3,N1))
  149. T4Z=(KXYZ(1,N2)-KXYZ(1,N1))*(KXYZ(2,N3)-KXYZ(2,N1))
  150. & -(KXYZ(2,N2)-KXYZ(2,N1))*(KXYZ(1,N3)-KXYZ(1,N1))
  151. C4=-(T4X*KXYZ(1,N1)+T4Y*KXYZ(2,N1)+T4Z*KXYZ(3,N1))
  152. P4C=T4X*KXYZ(1,N4)+T4Y*KXYZ(2,N4)+T4Z*KXYZ(3,N4)+C4
  153. T4X=T4X/P4C
  154. T4Y=T4Y/P4C
  155. T4Z=T4Z/P4C
  156. C4=C4/P4C
  157. C
  158. C Boucle sur les points de Gauss
  159. C
  160. DO 50 N=1,(NPG/8)
  161. IF (NPG.EQ.32) THEN
  162. DO 60 L=1,4
  163. AA(L)=BE
  164. 60 CONTINUE
  165. AA(N)=AL
  166. ELSE
  167. DO 70 L=1,4
  168. AA(L)=0.5D1
  169. 70 CONTINUE
  170. ENDIF
  171. C
  172. C Indice globale du point de Gauss
  173. C
  174. IGAU=(K-1)*(NPG/8)+N
  175. C
  176. C Coordonnees du point de Gauss (barycentre des 4 points du tetraedre)
  177. C
  178. X(IGAU)=AA(1)*KXYZ(1,N1)+AA(2)*KXYZ(1,N2)
  179. & +AA(3)*KXYZ(1,N3)+AA(4)*KXYZ(1,N4)
  180. Y(IGAU)=AA(1)*KXYZ(2,N1)+AA(2)*KXYZ(2,N2)
  181. & +AA(3)*KXYZ(2,N3)+AA(4)*KXYZ(2,N4)
  182. Z(IGAU)=AA(1)*KXYZ(3,N1)+AA(2)*KXYZ(3,N2)
  183. & +AA(3)*KXYZ(3,N3)+AA(4)*KXYZ(3,N4)
  184. C
  185. C Boucle sur les fonctions tests
  186. C
  187. C
  188. C Numerotation globale de la fonction test
  189. C
  190. IFN=(K-1)*4
  191. C
  192. C Calcul des fonctions tests et du gradient au point IGAU
  193. C
  194. FN((IFN+1),IGAU)=T1X*X(IGAU)+T1Y*Y(IGAU)
  195. & +T1Z*Z(IGAU)+C1
  196. GR(1,(IFN+1),IGAU)=T1X
  197. GR(2,(IFN+1),IGAU)=T1Y
  198. GR(3,(IFN+1),IGAU)=T1Z
  199. FN((IFN+2),IGAU)=T2X*X(IGAU)+T2Y*Y(IGAU)
  200. & +T2Z*Z(IGAU)+C2
  201. GR(1,(IFN+2),IGAU)=T2X
  202. GR(2,(IFN+2),IGAU)=T2Y
  203. GR(3,(IFN+2),IGAU)=T2Z
  204. FN((IFN+3),IGAU)=T3X*X(IGAU)+T3Y*Y(IGAU)
  205. & +T3Z*Z(IGAU)+C3
  206. GR(1,(IFN+3),IGAU)=T3X
  207. GR(2,(IFN+3),IGAU)=T3Y
  208. GR(3,(IFN+3),IGAU)=T3Z
  209. FN((IFN+4),IGAU)=T4X*X(IGAU)+T4Y*Y(IGAU)
  210. & +T4Z*Z(IGAU)+C4
  211. GR(1,(IFN+4),IGAU)=T4X
  212. GR(2,(IFN+4),IGAU)=T4Y
  213. GR(3,(IFN+4),IGAU)=T4Z
  214. C
  215. C Calcul des fonctions tests pour la pression
  216. C
  217. FM(((K+1)/2),IGAU)=1.0D0
  218. 50 CONTINUE
  219. 40 CONTINUE
  220.  
  221. IF(NOM2.EQ.'MCF1')THEN
  222. CO=6.D0 ** (1.D0/3.D0)
  223. UNSCO=1.D0/CO
  224. XXXX=-1.D0*UNSCO
  225. DO 51 L=1,NPG
  226. FM(1,L)=1.D0-(X(L)+Y(L)+Z(L))*UNSCO
  227. FM(2,L)=X(L)*UNSCO
  228. FM(3,L)=Y(L)*UNSCO
  229. FM(4,L)=Z(L)*UNSCO
  230. C
  231. GM(1,1,L)=XXXX
  232. GM(2,1,L)=XXXX
  233. GM(3,1,L)=XXXX
  234. C
  235. GM(1,2,L)=UNSCO
  236. GM(2,2,L)=0.D0
  237. GM(3,2,L)=0.D0
  238. C
  239. GM(1,3,L)=0.D0
  240. GM(2,3,L)=UNSCO
  241. GM(3,3,L)=0.D0
  242. C
  243. GM(1,4,L)=0.D0
  244. GM(2,4,L)=0.D0
  245. GM(3,4,L)=UNSCO
  246. C
  247. 51 CONTINUE
  248. ENDIF
  249.  
  250. C write(6,*)'X '
  251. C write(6,1008)X
  252. 1008 FORMAT(8(1X,1PE11.4))
  253.  
  254. RETURN
  255. END
  256.  
  257.  
  258.  

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