Télécharger pb702.eso

Retour à la liste

Numérotation des lignes :

  1. C PB702 SOURCE MAGN 10/05/31 21:15:16 6679
  2. SUBROUTINE PB702(XREF,X,Y,PG,FN,GR,FM,GM,ND,NP,MP,NPG,NOM2,ITYPI)
  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 : TRI7
  8. C
  9. C Cf Crouzeix Raviard. Conforming and nonconforming finite element
  10. C methods for solving the stationary Stokes equations I
  11. C R.A.I.R.O. (7e année, décembre 1973,R-3,p.33 à 76)
  12. C
  13. C
  14. C ^ eta
  15. C |
  16. C |n3
  17. C r2 5 1 -> A1
  18. C |\ 2 -> A12 -
  19. C | \ r2=sqrt(2) 3 -> A2
  20. C | \ 4 -> A23
  21. C | \ 5 -> A3
  22. C | \ 6 -> A13
  23. C | \ 7 -> A123
  24. C r2/2 6------4
  25. C | |\ D1=1-x/r2-y/r2
  26. C | | \ D2=x/r2
  27. C | 7 | \ D3=y/r2
  28. C | | \
  29. C | | \
  30. C 1------2-----3--->ksi
  31. C 0 r2/2 r2
  32. C
  33. C Pi = Di(aDi-c)+FD1xD2xD3 1 <= i <= 3
  34. C Pij = bxDi(aDi-c)-GxD1xD2xD3 1 <= i < j <= 3
  35. C P123= QxD1xD2xD3
  36. C Pour l'élément de référence on a a=2 b=4 c=1 F=3 G=12
  37. C Si on veut perturber un peu la position des points il faut néanmoins
  38. C vérifier les propriétés essentielles des fonctions de forme qui outre
  39. C Pi=1 en X=i et Pi=0 ailleur on doit avoir Somme des Pi = 1
  40. C (intégration exacte d'une fonction constante). Ceci conduit à des
  41. C relations entre a,b,c,F,G et Q
  42. C b=2a c=a-1 Q=3(G-F) (Somme Pi = 1) a ne 2 points Aij non au milieux
  43. C des arètes
  44. C
  45. C ITYPI=2 pts d'integration de type Gauss Lobatto i.e. sommets element
  46. C
  47. C REMARQUE :
  48. C En coordonnées cylindriques,pour une distribution régulière des noeuds
  49. C conforme à l'article de Crouzeix Raviard, la matrice masse diagonale
  50. C s'annule sur l'axe, comme pour le QUA9 (voir pb902.eso).
  51. C /1
  52. C | 2 pi x N1 dx = 0.
  53. C /0
  54. C Pour remédier à cela on se sert de la bulle en la décalant un peu par
  55. C rapport à l'élément de référence en jouant sur les constantes F et G.
  56. C On a pas à toucher "a" position des points 2,4 et 6.
  57. C Le choix que l'on a fait conduit à une matrice masse diagonale >0. En
  58. C conséquence ce choix doit être accompagné d'une numérotation dans le
  59. C sens trigo de l'élément courant.
  60. C
  61. C***********************************************************************
  62.  
  63. REAL*8 XREF(ND,NP),X(NPG),Y(NPG)
  64. DIMENSION FN(NP,NPG),GR(ND,NP,NPG),PG(NPG)
  65. DIMENSION FM(MP,NPG),GM(ND,MP,NPG)
  66. REAL*8 A,B,C,D,U(5),H(5)
  67. CHARACTER*4 NOM2
  68. -INC CCOPTIO
  69. C***
  70.  
  71. IF(IFOMOD.EQ.0)THEN
  72. F=3.0001
  73. G=11.9999
  74. F=3.
  75. G=12.
  76. ELSE
  77. F=3.
  78. G=12.
  79. ENDIF
  80. Q=3.*(G-F)
  81. a=2.
  82. b=2.*a
  83. c=a-1.
  84.  
  85.  
  86. R2=SQRT(2.D0)
  87. IPOSP=2
  88.  
  89.  
  90. XREF(1,1)=0.D0
  91. XREF(2,1)=0.D0
  92.  
  93. XREF(1,2)=R2/2.D0
  94. XREF(2,2)=0.D0
  95.  
  96. XREF(1,3)=R2
  97. XREF(2,3)=0.D0
  98.  
  99. XREF(1,4)=R2/2.D0
  100. XREF(2,4)=R2/2.D0
  101.  
  102. XREF(1,5)=0.D0
  103. XREF(2,5)=R2
  104.  
  105. XREF(1,6)=0.
  106. XREF(2,6)=R2/2.D0
  107.  
  108. XREF(1,7)=R2/3.D0
  109. XREF(2,7)=R2/3.D0
  110.  
  111. C Verification des coordonnées
  112. c IF(.TRUE.)THEN
  113. IF(.FALSE.)THEN
  114. DO 11 L=1,NP
  115. X(L)=XREF(1,L)
  116. Y(L)=XREF(2,L)
  117. 11 CONTINUE
  118.  
  119. DO 12 L=1,NP
  120.  
  121. D1=1.D0-X(L)/R2-Y(L)/R2
  122. D2=X(L)/R2
  123. D3=Y(L)/R2
  124.  
  125. FN(1,L)=D1*(a*D1-c)+F*D1*D2*D3
  126. FN(2,L)=b*D1*D2-G*D1*D2*D3
  127. FN(3,L)=D2*(a*D2-c)+F*D1*D2*D3
  128. FN(4,L)=b*D2*D3-G*D1*D2*D3
  129. FN(5,L)=D3*(a*D3-c)+F*D1*D2*D3
  130. FN(6,L)=b*D1*D3-G*D1*D2*D3
  131. FN(7,L)=Q*D1*D2*D3
  132. write(6,1033)L,FN(1,L),FN(2,L),FN(3,L),FN(4,L),FN(5,L),FN(6,L)
  133. & ,FN(7,L)
  134. 12 CONTINUE
  135. 1033 FORMAT(1X,I4,' FN',10(1X,1PD11.4))
  136. ENDIF
  137. C Fin Vérification
  138.  
  139. IF(ITYPI.EQ.2)THEN
  140. X(1)=0.D0
  141. Y(1)=0.D0
  142. X(2)=R2/2.D0
  143. Y(2)=0.D0
  144. X(3)=R2
  145. Y(3)=0.D0
  146. X(4)=R2/2.D0
  147. Y(4)=R2/2.D0
  148. X(5)=0.D0
  149. Y(5)=R2
  150. X(6)=0.D0
  151. Y(6)=R2/2.D0
  152. X(7)=(X(1)+X(3)+X(5))/3.D0
  153. Y(7)=(Y(1)+Y(3)+Y(5))/3.D0
  154. DO 2 L=1,7
  155. PG(L)=1.D0/7.D0
  156. 2 CONTINUE
  157. ELSE
  158. CALL CALUHH(X,Y,PG,NPG)
  159. ENDIF
  160.  
  161. DO 1 L=1,NPG
  162. C
  163. IF(NOM2.EQ.'PRP0')THEN
  164. FM(1,L)=1.D0
  165. GM(1,1,L)=0.D0
  166. GM(2,1,L)=0.D0
  167. ELSEIF(NOM2.EQ.'PRP1')THEN
  168. C? FM(1,L)=R2*(X(L)+Y(L)-R2/2.D0)
  169. C? FM(2,L)=-R2*(X(L)-R2/2.D0)
  170. C? FM(3,L)=-R2*(Y(L)-R2/2.D0)
  171. A1=6.D0/5.D0/R2
  172. FM(1,L)=5.D0/3.D0*(1.D0-A1*X(L)-A1*Y(L))
  173. FM(2,L)=-1.D0/3.D0*(1.D0-6.D0/R2*X(L))
  174. FM(3,L)=-1.D0/3.D0*(1.D0-6.D0/R2*Y(L))
  175. C
  176. GM(1,1,L)=0.D0
  177. GM(2,1,L)=0.D0
  178. GM(1,2,L)=0.D0
  179. GM(2,2,L)=0.D0
  180. GM(1,3,L)=0.D0
  181. GM(2,3,L)=0.D0
  182. ELSEIF(NOM2.EQ.'PFP1')THEN
  183. FM(1,L)=-R2/2.D0*(X(L)+Y(L)-R2)
  184. FM(2,L)=R2/2.D0*X(L)
  185. FM(3,L)=R2/2.D0*Y(L)
  186.  
  187. GM(1,1,L)=-R2/2.D0
  188. GM(2,1,L)=-R2/2.D0
  189. GM(1,2,L)=R2/2.D0
  190. GM(2,2,L)=0.D0
  191. GM(1,3,L)=0.D0
  192. GM(2,3,L)=R2/2.D0
  193. ENDIF
  194.  
  195. C
  196. c IF(.FALSE.)THEN
  197. IF(.TRUE.)THEN
  198. C
  199.  
  200. D1=1.D0-X(L)/R2-Y(L)/R2
  201. D2=X(L)/R2
  202. D3=Y(L)/R2
  203.  
  204. FN(1,L)=D1*(a*D1-c)+F*D1*D2*D3
  205. FN(2,L)=b*D1*D2-G*D1*D2*D3
  206. FN(3,L)=D2*(a*D2-c)+F*D1*D2*D3
  207. FN(4,L)=b*D2*D3-G*D1*D2*D3
  208. FN(5,L)=D3*(a*D3-c)+F*D1*D2*D3
  209. FN(6,L)=b*D1*D3-G*D1*D2*D3
  210. FN(7,L)=Q*D1*D2*D3
  211. c write(6,*)FN(1,L),FN(2,L),FN(3,L),FN(4,L),FN(5,L),FN(6,L),FN(7,L)
  212.  
  213. D1x=-1./R2
  214. D1y=-1./R2
  215.  
  216. D2x=1./R2
  217. D2y=0.
  218.  
  219. D3x=0.
  220. D3y=1/R2
  221.  
  222. GR(1,1,L)=D1x*(a*D1-c)+D1*a*D1x+F*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
  223. GR(2,1,L)=D1y*(a*D1-c)+D1*a*D1y+F*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
  224.  
  225. GR(1,2,L)=b*(D1x*D2+D1*D2x)-G*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
  226. GR(2,2,L)=b*(D1y*D2+D1*D2y)-G*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
  227.  
  228. GR(1,3,L)=D2x*(a*D2-c)+D2*a*D2x+F*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
  229. GR(2,3,L)=D2y*(a*D2-c)+D2*a*D2y+F*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
  230.  
  231. GR(1,4,L)=b*(D2x*D3+D2*D3x)-G*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
  232. GR(2,4,L)=b*(D2y*D3+D2*D3y)-G*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
  233.  
  234. GR(1,5,L)=D3x*(a*D3-c)+D3*a*D3x+F*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
  235. GR(2,5,L)=D3y*(a*D3-c)+D3*a*D3y+F*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
  236.  
  237. GR(1,6,L)=b*(D1x*D3+D1*D3x)-G*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
  238. GR(2,6,L)=b*(D1y*D3+D1*D3y)-G*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
  239.  
  240. GR(1,7,L)=Q*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
  241. GR(2,7,L)=Q*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
  242.  
  243. ELSE
  244.  
  245. c write(6,*)'La base'
  246.  
  247. C
  248. C
  249. FN(1,L)=(1.D0-X(L)/R2-Y(L)/R2)*(1.D0-2.D0*X(L)/R2-2.D0*Y(L)/R2+
  250. & 3.D0/2.D0*X(L)*Y(L))
  251. FN(2,L)=(1.D0-X(L)/R2-Y(L)/R2)*(4.D0/R2-6.D0*Y(L))*X(L)
  252. FN(3,L)=X(L)/R2*(2.D0*X(L)/R2-1.D0+3.D0*(1.D0-X(L)/R2-Y(L)/R2)
  253. & *Y(L)/R2)
  254. FN(4,L)=2.D0*X(L)*Y(L)*(1.D0-3.D0*(1.D0-X(L)/R2-Y(L)/R2))
  255. FN(5,L)=Y(L)*Y(L)-Y(L)*R2/2.D0+3.D0/2.D0*X(L)*Y(L)-
  256. & 3.D0*R2/4.D0*X(L)*X(L)*Y(L)-3.D0*R2/4.D0*X(L)*Y(L)*Y(L)
  257. FN(6,L)=2.D0*R2*(-R2/2.D0*Y(L)*Y(L)+3.D0/2.D0*(Y(L)*Y(L)*X(L)+
  258. & X(L)*X(L)*Y(L))-2.D0*R2*X(L)*Y(L)+Y(L))
  259. FN(7,L)=27.D0/2.D0*(X(L)*Y(L)-R2/2.D0*X(L)*X(L)*Y(L)-
  260. &R2/2.D0*X(L)*Y(L)*Y(L))
  261. C
  262. C
  263. C
  264. GR(1,1,L)=-3.D0*R2/4.D0*Y(L)*Y(L)+2.D0*X(L)+7./2.D0*Y(L)-
  265. * 3.D0*R2/2.D0*X(L)*Y(L)-3.D0*R2/2.D0
  266. GR(2,1,L)=-3.D0*R2/4.D0*X(L)*X(L)+2.D0*Y(L)+7./2.D0*X(L)-
  267. * 3.D0*R2/2.D0*X(L)*Y(L)-3.D0*R2/2.D0
  268. GR(1,2,L)=3.D0*R2*Y(L)*Y(L)+6.D0*R2*X(L)*Y(L)-4.D0*X(L)-8.D0*Y(L)+
  269. * 2.D0*R2
  270. GR(2,2,L)=3.D0*R2*X(L)*X(L)+6.D0*R2*X(L)*Y(L)-8.D0*X(L)
  271. GR(1,3,L)=-3.D0*R2/4.D0*Y(L)*Y(L)-3.D0*R2/2.D0*X(L)*Y(L)+
  272. * 2.D0*X(L)+3.D0/2.D0*Y(L)-R2/2.D0
  273. GR(2,3,L)=R2/2.D0*(-3.D0/2.D0*X(L)*X(L)-3.D0*X(L)*Y(L)+
  274. * 3.D0*R2/2.D0*X(L))
  275. GR(1,4,L)=3.D0*R2*Y(L)*Y(L)+6.D0*R2*X(L)*Y(L)-4.D0*Y(L)
  276. GR(2,4,L)=3.D0*R2*X(L)*X(L)+6.D0*R2*X(L)*Y(L)-4.D0*X(L)
  277. GR(1,5,L)=-3.D0*R2/4.D0*Y(L)*Y(L)-3.D0*R2/2.D0*X(L)*Y(L)
  278. & +3.D0/2.D0*Y(L)
  279. GR(2,5,L)=-3.D0*R2/4.D0*X(L)*X(L)-3.D0*R2/2.D0*X(L)*Y(L)
  280. & +2.D0*Y(L)+
  281. * 3.D0/2.D0*X(L)-R2/2.D0
  282. GR(1,6,L)=3.D0*R2*Y(L)*Y(L)+6.D0*R2*X(L)*Y(L)-8.D0*Y(L)
  283. GR(2,6,L)=3.D0*R2*X(L)*X(L)+6.D0*R2*X(L)*Y(L)-8.D0*X(L)
  284. & -4.D0*Y(L)+2.D0*R2
  285. GR(1,7,L)=27.D0/2.D0*(-R2/2*Y(L)*Y(L)-R2*X(L)*Y(L)+Y(L))
  286. GR(2,7,L)=27.D0/2.D0*(-R2/2*X(L)*X(L)-R2*X(L)*Y(L)+X(L))
  287. C
  288. ENDIF
  289. 1 CONTINUE
  290.  
  291. IF(.FALSE.)THEN
  292. c IF(.TRUE.)THEN
  293. DO 121 L=1,NPG
  294. F1=0.
  295. DO 131 I=1,NP
  296. F1=F1+FN(I,L)
  297. 131 CONTINUE
  298. Write(6,*)' L=',L,F1
  299. 121 CONTINUE
  300.  
  301. DO 122 L=1,NPG
  302. G1=0.
  303. G2=0.
  304. DO 141 I=1,NP
  305. G1=G1+GR(1,I,L)
  306. G2=G2+GR(2,I,L)
  307. 141 CONTINUE
  308. Write(6,*)' L=',L,G1,G2
  309. 122 CONTINUE
  310. ENDIF
  311.  
  312. c stop
  313.  
  314.  
  315. C
  316. C
  317. C WRITE(6,101)
  318. C WRITE(6,1002)FM
  319. C WRITE(6,1002)GM
  320. C WRITE(6,1002)FN
  321. C WRITE(6,1002)GR
  322. C WRITE(6,101)
  323.  
  324. RETURN
  325. 1002 FORMAT(10(1X,1PD11.4))
  326. 1001 FORMAT(20(1X,I5))
  327. 101 FORMAT(1X,'... SUBPB702 ... FM,GM,FN,GR ',9(10H..........)/)
  328. C
  329. END
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  
  338.  

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