Télécharger pb902.eso

Retour à la liste

Numérotation des lignes :

  1. C PB902 SOURCE MAGN 10/06/02 21:15:00 6681
  2. SUBROUTINE PB902
  3. &(XREF,X,Y,PG,FN,GR,FM,GM,ND,NP,MP,NG,NPG,NOM2,ITYPI)
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. C************************************************************************
  7. C
  8. C CALCULE LES FONCTIONS DE FORME D'UN : QUA9
  9. C eta
  10. C ^
  11. C |
  12. C 7 6 5
  13. C 1. x------x------x
  14. C | |
  15. C | |
  16. C alfa 8x 9x 4x
  17. C | |
  18. C | |
  19. C 0. x------x------x--> ksi
  20. C 1 2 3
  21. C
  22. C 0. alfa 1.
  23. C
  24. C ITYPI=2 pts d'intégration de type Gauss Lobatto i.e. sommets élément
  25. C
  26. C REMARQUE : Pour l'élément quadratique.
  27. C En coordonnées cylindriques si les points milieux (2,4,6,8,9) sont
  28. C pile poil au milieu, la matrice masse diagonale s'annule sur l'axe.
  29. C /1
  30. C | 2 pi x N1 dx = 0.
  31. C /0
  32. C Pour remédier à cela on décale un peu les points mileux sur l'élément
  33. C de référence, alfa=0.50005 au lieu de 0.5 .
  34. C Si on rapprochait les points milieux de l'axe (alfa = 0.45 par exemple)
  35. C la matrice masse diagonale serait négative sur l'axe, ce que l'on ne
  36. C veut pas non plus. En conséquence le choix qu'on a fait doit être
  37. C accompagné d'une numérotation dans le sens trigo de l'élément courant.
  38. C************************************************************************
  39. REAL*8 XREF(ND,NP),X(NPG),Y(NPG)
  40. DIMENSION FN(NP,NPG),GR(ND,NP,NPG),PG(NPG)
  41. DIMENSION FM(MP,NPG),GM(ND,MP,NPG)
  42. CHARACTER*4 NOM2
  43. REAL*8 A,B,C,D,U(5),H(5),AUX1,AUX2
  44. -INC CCOPTIO
  45.  
  46. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  47. AUX1=1.D0/2.D0+3.D0**(1.D0/2.D0)/6.D0
  48. AUX2=1.D0/2.D0-3.D0**(1.D0/2.D0)/6.D0
  49. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  50. IF(IFOMOD.EQ.0)THEN
  51. alfa=0.50005D0
  52. ELSE
  53. alfa=0.5D0
  54. ENDIF
  55.  
  56. alfai=1.D0/alfa
  57. alfa2=1.D0/alfa/(alfa - 1.D0)
  58. alfa3=1.D0/(alfai - 1.D0)
  59.  
  60. XREF(1,1)=0.D0
  61. XREF(2,1)=0.D0
  62.  
  63. XREF(1,2)=alfa
  64. XREF(2,2)=0.D0
  65.  
  66. XREF(1,3)=1.D0
  67. XREF(2,3)=0.D0
  68.  
  69. XREF(1,4)=1.D0
  70. XREF(2,4)=alfa
  71.  
  72. XREF(1,5)=1.D0
  73. XREF(2,5)=1.D0
  74.  
  75. XREF(1,6)=alfa
  76. XREF(2,6)=1.D0
  77.  
  78. XREF(1,7)=0.D0
  79. XREF(2,7)=1.D0
  80.  
  81. XREF(1,8)=0.D0
  82. XREF(2,8)=alfa
  83.  
  84. XREF(1,9)=alfa
  85. XREF(2,9)=alfa
  86.  
  87. C Verification des coordonnées
  88. c IF(.TRUE.)THEN
  89. IF(.FALSE.)THEN
  90. DO 11 L=1,NP
  91. X(L)=XREF(1,L)
  92. Y(L)=XREF(2,L)
  93. 11 CONTINUE
  94.  
  95. DO 12 L=1,NP
  96. FN(1,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  97. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  98. FN(2,L)=alfa2*X(L)*(X(L)-1.D0)
  99. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  100. FN(3,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  101. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  102. FN(4,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  103. & *alfa2*Y(L)*(Y(L)-1.D0)
  104. FN(5,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  105. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  106. FN(6,L)=alfa2*X(L)*(X(L)-1.D0)
  107. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  108. FN(7,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  109. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  110. FN(8,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  111. & *alfa2*Y(L)*(Y(L)-1.D0)
  112. FN(9,L)=alfa2*X(L)*(X(L)-1.D0)
  113. & *alfa2*Y(L)*(Y(L)-1.D0)
  114.  
  115. write(6,1033)L,FN(1,L),FN(2,L),FN(3,L),FN(4,L),FN(5,L),FN(6,L)
  116. & ,FN(7,L),FN(8,L),FN(9,L)
  117. 12 CONTINUE
  118. 1033 FORMAT(1X,I4,' FN',10(1X,1PD11.4))
  119. ENDIF
  120. C Fin Vérification
  121.  
  122. CALL CALUHG(U,H,NG)
  123. C
  124. A=0.D0
  125. B=1.D0
  126. C=0.D0
  127. D=1.D0
  128. IF(ITYPI.EQ.2)THEN
  129. X(1)=0.D0
  130. Y(1)=0.D0
  131. X(2)=alfa
  132. Y(2)=0.D0
  133. X(3)=1.D0
  134. Y(3)=0.D0
  135. X(4)=1.D0
  136. Y(4)=alfa
  137. X(5)=1.D0
  138. Y(5)=1.D0
  139. X(6)=alfa
  140. Y(6)=1.D0
  141. X(7)=0.D0
  142. Y(7)=1.D0
  143. X(8)=0.D0
  144. Y(8)=alfa
  145. X(9)=alfa
  146. Y(9)=alfa
  147. DO 2 L=1,NPG
  148. PG(L)=1.D0/9.D0
  149. 2 CONTINUE
  150. ELSE
  151. CALL CALG2(A,B,C,D,NG,H,U,X,Y,PG)
  152. ENDIF
  153.  
  154. DO 1 L=1,NPG
  155. C
  156. IF(NOM2.EQ.'PRP0')THEN
  157. FM(1,L)= 1.D0
  158.  
  159. GM(1,1,L)=0.D0
  160. GM(2,1,L)=0.D0
  161.  
  162. ELSEIF(NOM2.EQ.'PRQ1')THEN
  163. FM(1,L)= 3.D0*(X(L)-AUX1)*(Y(L)-AUX1)
  164. FM(2,L)=-3.D0*(X(L)-AUX2)*(Y(L)-AUX1)
  165. FM(3,L)= 3.D0*(X(L)-AUX2)*(Y(L)-AUX2)
  166. FM(4,L)=-3.D0*(X(L)-AUX1)*(Y(L)-AUX2)
  167. C
  168. GM(1,1,L)=3.D0*(Y(L)-AUX1)
  169. GM(2,1,L)=3.D0*(X(L)-AUX1)
  170. GM(1,2,L)=-3.D0*(Y(L)-AUX1)
  171. GM(2,2,L)=-3.D0*(X(L)-AUX2)
  172. GM(1,3,L)=3.D0*(Y(L)-AUX2)
  173. GM(2,3,L)=3.D0*(X(L)-AUX2)
  174. GM(1,4,L)=-3.D0*(Y(L)-AUX2)
  175. GM(2,4,L)=-3.D0*(X(L)-AUX1)
  176.  
  177. ELSEIF(NOM2.EQ.'PRP1')THEN
  178. C? FM(1,L)= (3.D0-4.D0*X(L))/2.D0
  179. C? FM(2,L)= 2.D0*(X(L)-Y(L))
  180. C? FM(3,L)= (4.D0*Y(L)-1.D0)/2.D0
  181. FM(1,L)= (4.D0*X(L)+4.D0*Y(L)-3.D0)/3.D0
  182. FM(2,L)=-1.D0/3.D0*(8.D0*X(L)-4.D0*Y(L)-3.D0)
  183. FM(3,L)= 1.D0/3.D0*(-8.D0*Y(L)+4.D0*X(L)+3.D0)
  184. C
  185. GM(1,1,L)=0.D0
  186. GM(2,1,L)=0.D0
  187. GM(1,2,L)=0.D0
  188. GM(2,2,L)=0.D0
  189. GM(1,3,L)=0.D0
  190. GM(2,3,L)=0.D0
  191. ELSEIF(NOM2.EQ.'PFP1')THEN
  192. FM(1,L)= (X(L)-1.D0)*(Y(L)-1.D0)
  193. FM(2,L)=-X(L)*(Y(L)-1.D0)
  194. FM(3,L)= X(L)*Y(L)
  195. FM(4,L)=-Y(L)*(X(L)-1.D0)
  196.  
  197. GM(1,1,L)=Y(L)-1.D0
  198. GM(2,1,L)=X(L)-1.D0
  199. GM(1,2,L)=-(Y(L)-1.D0)
  200. GM(2,2,L)=-X(L)
  201. GM(1,3,L)=Y(L)
  202. GM(2,3,L)=X(L)
  203. GM(1,4,L)=-Y(L)
  204. GM(2,4,L)=-(X(L)-1.D0)
  205. ENDIF
  206. C
  207. c N1X=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  208. c N1Y=(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  209. c
  210. c N2X=alfa2*X(L)*(X(L)-1.D0)
  211. c N2Y=alfa2*Y(L)*(Y(L)-1.D0)
  212. c
  213. c N3X=alfa3*X(L)*(alfai*X(L)-1.D0)
  214. c N3Y=alfa3*Y(L)*(alfai*Y(L)-1.D0)
  215. C
  216.  
  217. c FN(1,L)=N1X*N1Y
  218. FN(1,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  219. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  220. c FN(2,L)=N2X*N1Y
  221. FN(2,L)=alfa2*X(L)*(X(L)-1.D0)
  222. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  223. c FN(3,L)=N3X*N1Y
  224. FN(3,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  225. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  226. c FN(4,L)=N3X*N2Y
  227. FN(4,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  228. & *alfa2*Y(L)*(Y(L)-1.D0)
  229. c FN(5,L)=N3X*N3Y
  230. FN(5,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  231. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  232. c FN(6,L)=N2X*N3Y
  233. FN(6,L)=alfa2*X(L)*(X(L)-1.D0)
  234. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  235. c FN(7,L)=N1X*N3Y
  236. FN(7,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  237. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  238. c FN(8,L)=N1X*N2Y
  239. FN(8,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  240. & *alfa2*Y(L)*(Y(L)-1.D0)
  241. c FN(9,L)=N2X*N2Y
  242. FN(9,L)=alfa2*X(L)*(X(L)-1.D0)
  243. & *alfa2*Y(L)*(Y(L)-1.D0)
  244.  
  245. C
  246. c FN(1,L)=(X(L)-1.D0)*(Y(L)-1.D0)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)
  247. c FN(2,L)=-4.D0*X(L)*(X(L)-1.D0)*(2.D0*Y(L)-1.D0)*(Y(L)-1.D0)
  248. c FN(3,L)=X(L)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)*(Y(L)-1.D0)
  249. c FN(4,L)=-4.D0*X(L)*Y(L)*(2.D0*X(L)-1.D0)*(Y(L)-1.D0)
  250. c FN(5,L)=X(L)*Y(L)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)
  251. c FN(6,L)=-4.D0*X(L)*Y(L)*(X(L)-1.D0)*(2.D0*Y(L)-1.D0)
  252. c FN(7,L)=Y(L)*(2.D0*Y(L)-1.D0)*(2.D0*X(L)-1.D0)*(X(L)-1.D0)
  253. c FN(8,L)=-4.D0*Y(L)*(Y(L)-1.D0)*(X(L)-1.D0)*(2.D0*X(L)-1.D0)
  254. c FN(9,L)=16.D0*X(L)*Y(L)*(X(L)-1.D0)*(Y(L)-1.D0)
  255. C
  256. c N1X=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  257. c N1Y=(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  258. c
  259. c N2X=alfa2*X(L)*(X(L)-1.D0)
  260. c N2Y=alfa2*Y(L)*(Y(L)-1.D0)
  261. c
  262. c N3X=alfa3*X(L)*(alfai*X(L)-1.D0)
  263. c N3Y=alfa3*Y(L)*(alfai*Y(L)-1.D0)
  264. c
  265. c G1X=2.D0*alfai*X(L)-(1.D0 + alfai)
  266. c G1Y=2.D0*alfai*Y(L)-(1.D0 + alfai)
  267. c
  268. c G2X=alfa2*(2.D0*X(L)-1.D0)
  269. c G2Y=alfa2*(2.D0*Y(L)-1.D0)
  270. c
  271. c G3X=alfa3*(2.D0*alfai*X(L)-1.D0)
  272. c G3Y=alfa3*(2.D0*alfai*Y(L)-1.D0)
  273.  
  274. C
  275. c GR(1,1,L)=G1X*N1Y
  276. GR(1,1,L)=(2.D0*alfai*X(L)-(1.D0 + alfai))
  277. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  278. c GR(2,1,L)=N1X*G1Y
  279. GR(2,1,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  280. & *(2.D0*alfai*Y(L)-(1.D0 + alfai))
  281. c GR(1,2,L)=G2X*N1Y
  282. GR(1,2,L)=(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  283. & *alfa2*(2.D0*X(L)-1.D0)
  284. c GR(2,2,L)=N2X*G1Y
  285. GR(2,2,L)=alfa2*X(L)*(X(L)-1.D0)
  286. & *(2.D0*alfai*Y(L)-(1.D0 + alfai))
  287. c GR(1,3,L)=G3X*N1Y
  288. GR(1,3,L)=alfa3*(2.D0*alfai*X(L)-1.D0)
  289. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  290. c GR(2,3,L)=N3X*G1Y
  291. GR(2,3,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  292. & *(2.D0*alfai*Y(L)-(1.D0 + alfai))
  293. c GR(1,4,L)=G3X*N2Y
  294. GR(1,4,L)=alfa3*(2.D0*alfai*X(L)-1.D0)
  295. & *alfa2*Y(L)*(Y(L)-1.D0)
  296. c GR(2,4,L)=N3X*G2Y
  297. GR(2,4,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  298. & *alfa2*(2.D0*Y(L)-1.D0)
  299. c GR(1,5,L)=G3X*N3Y
  300. GR(1,5,L)=alfa3*(2.D0*alfai*X(L)-1.D0)
  301. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  302. c GR(2,5,L)=N3X*G3Y
  303. GR(2,5,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  304. & *alfa3*(2.D0*alfai*Y(L)-1.D0)
  305. c GR(1,6,L)=G2X*N3Y
  306. GR(1,6,L)=alfa2*(2.D0*X(L)-1.D0)
  307. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  308. c GR(2,6,L)=N2X*G3Y
  309. GR(2,6,L)=alfa2*X(L)*(X(L)-1.D0)
  310. & *alfa3*(2.D0*alfai*Y(L)-1.D0)
  311. c GR(1,7,L)=G1X*N3Y
  312. GR(1,7,L)=(2.D0*alfai*X(L)-(1.D0 + alfai))
  313. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  314. c GR(2,7,L)=N1X*G3Y
  315. GR(2,7,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  316. & *alfa3*(2.D0*alfai*Y(L)-1.D0)
  317. c GR(1,8,L)=G1X*N2Y
  318. GR(1,8,L)=(2.D0*alfai*X(L)-(1.D0 + alfai))
  319. & *alfa2*Y(L)*(Y(L)-1.D0)
  320. c GR(2,8,L)=N1X*G2Y
  321. GR(2,8,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  322. & *alfa2*(2.D0*Y(L)-1.D0)
  323. c GR(1,9,L)=G2X*N2Y
  324. GR(1,9,L)=alfa2*Y(L)*(Y(L)-1.D0)
  325. & *alfa2*(2.D0*X(L)-1.D0)
  326. c GR(2,9,L)=N2X*G2Y
  327. GR(2,9,L)=alfa2*X(L)*(X(L)-1.D0)
  328. & *alfa2*(2.D0*Y(L)-1.D0)
  329. C
  330. c GR(1,1,L)=(Y(L)-1.D0)*(2.D0*Y(L)-1.D0)*(4.D0*X(L)-3.D0)
  331. c GR(2,1,L)=(X(L)-1.D0)*(2.D0*X(L)-1.D0)*(4.D0*Y(L)-3.D0)
  332. c GR(1,2,L)=-4.D0*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)*(Y(L)-1.D0)
  333. c GR(2,2,L)=-4.D0*X(L)*(X(L)-1.D0)*(4.D0*Y(L)-3.D0)
  334. c GR(1,3,L)=(2.D0*Y(L)-1.D0)*(Y(L)-1.D0)*(4.D0*X(L)-1.D0)
  335. c GR(2,3,L)=X(L)*(2.D0*X(L)-1.D0)*(4.D0*Y(L)-3.D0)
  336. c GR(1,4,L)=-4.D0*Y(L)*(Y(L)-1.D0)*(4.D0*X(L)-1.D0)
  337. c GR(2,4,L)=-4.D0*X(L)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)
  338. c GR(1,5,L)=Y(L)*(2.D0*Y(L)-1.D0)*(4.D0*X(L)-1.D0)
  339. c GR(2,5,L)=X(L)*(2.D0*X(L)-1.D0)*(4.D0*Y(L)-1.D0)
  340. c GR(1,6,L)=-4.D0*Y(L)*(2.D0*Y(L)-1.D0)*(2.D0*X(L)-1.D0)
  341. c GR(2,6,L)=-4.D0*X(L)*(X(L)-1.D0)*(4.D0*Y(L)-1.D0)
  342. c GR(1,7,L)=Y(L)*(2.D0*Y(L)-1.D0)*(4.D0*X(L)-3.D0)
  343. c GR(2,7,L)=(2.D0*X(L)-1.D0)*(X(L)-1.D0)*(4.D0*Y(L)-1.D0)
  344. c GR(1,8,L)=-4.D0*Y(L)*(Y(L)-1.D0)*(4.D0*X(L)-3.D0)
  345. c GR(2,8,L)=-4.D0*(X(L)-1.D0)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)
  346. c GR(1,9,L)=16.D0*Y(L)*(Y(L)-1.D0)*(2.D0*X(L)-1.D0)
  347. c GR(2,9,L)=16.D0*X(L)*(X(L)-1.D0)*(2.D0*Y(L)-1.D0)
  348. C
  349. 1 CONTINUE
  350. C
  351. C
  352. C WRITE(6,101)
  353. C WRITE(6,1002)FM
  354. C WRITE(6,1002)GM
  355. C WRITE(6,1002)FN
  356. C WRITE(6,1002)GR
  357. C WRITE(6,101)
  358.  
  359. C write(6,*)' RET PB902 '
  360. RETURN
  361. 1002 FORMAT(10(1X,1PD11.4))
  362. 1001 FORMAT(20(1X,I5))
  363. 101 FORMAT(1X,'... SUBPB902 ... FM,GM,FN,GR ',9(10H..........)/)
  364. C
  365. END
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  

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