Télécharger pb902.eso

Retour à la liste

Numérotation des lignes :

pb902
  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.  
  45. -INC PPARAM
  46. -INC CCOPTIO
  47.  
  48. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  49. AUX1=1.D0/2.D0+3.D0**(1.D0/2.D0)/6.D0
  50. AUX2=1.D0/2.D0-3.D0**(1.D0/2.D0)/6.D0
  51. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  52. IF(IFOMOD.EQ.0)THEN
  53. alfa=0.50005D0
  54. ELSE
  55. alfa=0.5D0
  56. ENDIF
  57.  
  58. alfai=1.D0/alfa
  59. alfa2=1.D0/alfa/(alfa - 1.D0)
  60. alfa3=1.D0/(alfai - 1.D0)
  61.  
  62. XREF(1,1)=0.D0
  63. XREF(2,1)=0.D0
  64.  
  65. XREF(1,2)=alfa
  66. XREF(2,2)=0.D0
  67.  
  68. XREF(1,3)=1.D0
  69. XREF(2,3)=0.D0
  70.  
  71. XREF(1,4)=1.D0
  72. XREF(2,4)=alfa
  73.  
  74. XREF(1,5)=1.D0
  75. XREF(2,5)=1.D0
  76.  
  77. XREF(1,6)=alfa
  78. XREF(2,6)=1.D0
  79.  
  80. XREF(1,7)=0.D0
  81. XREF(2,7)=1.D0
  82.  
  83. XREF(1,8)=0.D0
  84. XREF(2,8)=alfa
  85.  
  86. XREF(1,9)=alfa
  87. XREF(2,9)=alfa
  88.  
  89. C Verification des coordonnées
  90. c IF(.TRUE.)THEN
  91. IF(.FALSE.)THEN
  92. DO 11 L=1,NP
  93. X(L)=XREF(1,L)
  94. Y(L)=XREF(2,L)
  95. 11 CONTINUE
  96.  
  97. DO 12 L=1,NP
  98. FN(1,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  99. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  100. FN(2,L)=alfa2*X(L)*(X(L)-1.D0)
  101. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  102. FN(3,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  103. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  104. FN(4,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  105. & *alfa2*Y(L)*(Y(L)-1.D0)
  106. FN(5,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  107. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  108. FN(6,L)=alfa2*X(L)*(X(L)-1.D0)
  109. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  110. FN(7,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  111. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  112. FN(8,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  113. & *alfa2*Y(L)*(Y(L)-1.D0)
  114. FN(9,L)=alfa2*X(L)*(X(L)-1.D0)
  115. & *alfa2*Y(L)*(Y(L)-1.D0)
  116.  
  117. write(6,1033)L,FN(1,L),FN(2,L),FN(3,L),FN(4,L),FN(5,L),FN(6,L)
  118. & ,FN(7,L),FN(8,L),FN(9,L)
  119. 12 CONTINUE
  120. 1033 FORMAT(1X,I4,' FN',10(1X,1PD11.4))
  121. ENDIF
  122. C Fin Vérification
  123.  
  124. CALL CALUHG(U,H,NG)
  125. C
  126. A=0.D0
  127. B=1.D0
  128. C=0.D0
  129. D=1.D0
  130. IF(ITYPI.EQ.2)THEN
  131. X(1)=0.D0
  132. Y(1)=0.D0
  133. X(2)=alfa
  134. Y(2)=0.D0
  135. X(3)=1.D0
  136. Y(3)=0.D0
  137. X(4)=1.D0
  138. Y(4)=alfa
  139. X(5)=1.D0
  140. Y(5)=1.D0
  141. X(6)=alfa
  142. Y(6)=1.D0
  143. X(7)=0.D0
  144. Y(7)=1.D0
  145. X(8)=0.D0
  146. Y(8)=alfa
  147. X(9)=alfa
  148. Y(9)=alfa
  149. DO 2 L=1,NPG
  150. PG(L)=1.D0/9.D0
  151. 2 CONTINUE
  152. ELSE
  153. CALL CALG2(A,B,C,D,NG,H,U,X,Y,PG)
  154. ENDIF
  155.  
  156. DO 1 L=1,NPG
  157. C
  158. IF(NOM2.EQ.'PRP0')THEN
  159. FM(1,L)= 1.D0
  160.  
  161. GM(1,1,L)=0.D0
  162. GM(2,1,L)=0.D0
  163.  
  164. ELSEIF(NOM2.EQ.'PRQ1')THEN
  165. FM(1,L)= 3.D0*(X(L)-AUX1)*(Y(L)-AUX1)
  166. FM(2,L)=-3.D0*(X(L)-AUX2)*(Y(L)-AUX1)
  167. FM(3,L)= 3.D0*(X(L)-AUX2)*(Y(L)-AUX2)
  168. FM(4,L)=-3.D0*(X(L)-AUX1)*(Y(L)-AUX2)
  169. C
  170. GM(1,1,L)=3.D0*(Y(L)-AUX1)
  171. GM(2,1,L)=3.D0*(X(L)-AUX1)
  172. GM(1,2,L)=-3.D0*(Y(L)-AUX1)
  173. GM(2,2,L)=-3.D0*(X(L)-AUX2)
  174. GM(1,3,L)=3.D0*(Y(L)-AUX2)
  175. GM(2,3,L)=3.D0*(X(L)-AUX2)
  176. GM(1,4,L)=-3.D0*(Y(L)-AUX2)
  177. GM(2,4,L)=-3.D0*(X(L)-AUX1)
  178.  
  179. ELSEIF(NOM2.EQ.'PRP1')THEN
  180. C? FM(1,L)= (3.D0-4.D0*X(L))/2.D0
  181. C? FM(2,L)= 2.D0*(X(L)-Y(L))
  182. C? FM(3,L)= (4.D0*Y(L)-1.D0)/2.D0
  183. FM(1,L)= (4.D0*X(L)+4.D0*Y(L)-3.D0)/3.D0
  184. FM(2,L)=-1.D0/3.D0*(8.D0*X(L)-4.D0*Y(L)-3.D0)
  185. FM(3,L)= 1.D0/3.D0*(-8.D0*Y(L)+4.D0*X(L)+3.D0)
  186. C
  187. GM(1,1,L)=0.D0
  188. GM(2,1,L)=0.D0
  189. GM(1,2,L)=0.D0
  190. GM(2,2,L)=0.D0
  191. GM(1,3,L)=0.D0
  192. GM(2,3,L)=0.D0
  193. ELSEIF(NOM2.EQ.'PFP1')THEN
  194. FM(1,L)= (X(L)-1.D0)*(Y(L)-1.D0)
  195. FM(2,L)=-X(L)*(Y(L)-1.D0)
  196. FM(3,L)= X(L)*Y(L)
  197. FM(4,L)=-Y(L)*(X(L)-1.D0)
  198.  
  199. GM(1,1,L)=Y(L)-1.D0
  200. GM(2,1,L)=X(L)-1.D0
  201. GM(1,2,L)=-(Y(L)-1.D0)
  202. GM(2,2,L)=-X(L)
  203. GM(1,3,L)=Y(L)
  204. GM(2,3,L)=X(L)
  205. GM(1,4,L)=-Y(L)
  206. GM(2,4,L)=-(X(L)-1.D0)
  207. ENDIF
  208. C
  209. c N1X=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  210. c N1Y=(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  211. c
  212. c N2X=alfa2*X(L)*(X(L)-1.D0)
  213. c N2Y=alfa2*Y(L)*(Y(L)-1.D0)
  214. c
  215. c N3X=alfa3*X(L)*(alfai*X(L)-1.D0)
  216. c N3Y=alfa3*Y(L)*(alfai*Y(L)-1.D0)
  217. C
  218.  
  219. c FN(1,L)=N1X*N1Y
  220. FN(1,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  221. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  222. c FN(2,L)=N2X*N1Y
  223. FN(2,L)=alfa2*X(L)*(X(L)-1.D0)
  224. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  225. c FN(3,L)=N3X*N1Y
  226. FN(3,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  227. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  228. c FN(4,L)=N3X*N2Y
  229. FN(4,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  230. & *alfa2*Y(L)*(Y(L)-1.D0)
  231. c FN(5,L)=N3X*N3Y
  232. FN(5,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  233. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  234. c FN(6,L)=N2X*N3Y
  235. FN(6,L)=alfa2*X(L)*(X(L)-1.D0)
  236. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  237. c FN(7,L)=N1X*N3Y
  238. FN(7,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  239. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  240. c FN(8,L)=N1X*N2Y
  241. FN(8,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  242. & *alfa2*Y(L)*(Y(L)-1.D0)
  243. c FN(9,L)=N2X*N2Y
  244. FN(9,L)=alfa2*X(L)*(X(L)-1.D0)
  245. & *alfa2*Y(L)*(Y(L)-1.D0)
  246.  
  247. C
  248. c FN(1,L)=(X(L)-1.D0)*(Y(L)-1.D0)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)
  249. c FN(2,L)=-4.D0*X(L)*(X(L)-1.D0)*(2.D0*Y(L)-1.D0)*(Y(L)-1.D0)
  250. c FN(3,L)=X(L)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)*(Y(L)-1.D0)
  251. c FN(4,L)=-4.D0*X(L)*Y(L)*(2.D0*X(L)-1.D0)*(Y(L)-1.D0)
  252. c FN(5,L)=X(L)*Y(L)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)
  253. c FN(6,L)=-4.D0*X(L)*Y(L)*(X(L)-1.D0)*(2.D0*Y(L)-1.D0)
  254. c FN(7,L)=Y(L)*(2.D0*Y(L)-1.D0)*(2.D0*X(L)-1.D0)*(X(L)-1.D0)
  255. c FN(8,L)=-4.D0*Y(L)*(Y(L)-1.D0)*(X(L)-1.D0)*(2.D0*X(L)-1.D0)
  256. c FN(9,L)=16.D0*X(L)*Y(L)*(X(L)-1.D0)*(Y(L)-1.D0)
  257. C
  258. c N1X=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  259. c N1Y=(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  260. c
  261. c N2X=alfa2*X(L)*(X(L)-1.D0)
  262. c N2Y=alfa2*Y(L)*(Y(L)-1.D0)
  263. c
  264. c N3X=alfa3*X(L)*(alfai*X(L)-1.D0)
  265. c N3Y=alfa3*Y(L)*(alfai*Y(L)-1.D0)
  266. c
  267. c G1X=2.D0*alfai*X(L)-(1.D0 + alfai)
  268. c G1Y=2.D0*alfai*Y(L)-(1.D0 + alfai)
  269. c
  270. c G2X=alfa2*(2.D0*X(L)-1.D0)
  271. c G2Y=alfa2*(2.D0*Y(L)-1.D0)
  272. c
  273. c G3X=alfa3*(2.D0*alfai*X(L)-1.D0)
  274. c G3Y=alfa3*(2.D0*alfai*Y(L)-1.D0)
  275.  
  276. C
  277. c GR(1,1,L)=G1X*N1Y
  278. GR(1,1,L)=(2.D0*alfai*X(L)-(1.D0 + alfai))
  279. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  280. c GR(2,1,L)=N1X*G1Y
  281. GR(2,1,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  282. & *(2.D0*alfai*Y(L)-(1.D0 + alfai))
  283. c GR(1,2,L)=G2X*N1Y
  284. GR(1,2,L)=(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  285. & *alfa2*(2.D0*X(L)-1.D0)
  286. c GR(2,2,L)=N2X*G1Y
  287. GR(2,2,L)=alfa2*X(L)*(X(L)-1.D0)
  288. & *(2.D0*alfai*Y(L)-(1.D0 + alfai))
  289. c GR(1,3,L)=G3X*N1Y
  290. GR(1,3,L)=alfa3*(2.D0*alfai*X(L)-1.D0)
  291. & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
  292. c GR(2,3,L)=N3X*G1Y
  293. GR(2,3,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  294. & *(2.D0*alfai*Y(L)-(1.D0 + alfai))
  295. c GR(1,4,L)=G3X*N2Y
  296. GR(1,4,L)=alfa3*(2.D0*alfai*X(L)-1.D0)
  297. & *alfa2*Y(L)*(Y(L)-1.D0)
  298. c GR(2,4,L)=N3X*G2Y
  299. GR(2,4,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  300. & *alfa2*(2.D0*Y(L)-1.D0)
  301. c GR(1,5,L)=G3X*N3Y
  302. GR(1,5,L)=alfa3*(2.D0*alfai*X(L)-1.D0)
  303. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  304. c GR(2,5,L)=N3X*G3Y
  305. GR(2,5,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
  306. & *alfa3*(2.D0*alfai*Y(L)-1.D0)
  307. c GR(1,6,L)=G2X*N3Y
  308. GR(1,6,L)=alfa2*(2.D0*X(L)-1.D0)
  309. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  310. c GR(2,6,L)=N2X*G3Y
  311. GR(2,6,L)=alfa2*X(L)*(X(L)-1.D0)
  312. & *alfa3*(2.D0*alfai*Y(L)-1.D0)
  313. c GR(1,7,L)=G1X*N3Y
  314. GR(1,7,L)=(2.D0*alfai*X(L)-(1.D0 + alfai))
  315. & *alfa3*Y(L)*(alfai*Y(L)-1.D0)
  316. c GR(2,7,L)=N1X*G3Y
  317. GR(2,7,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  318. & *alfa3*(2.D0*alfai*Y(L)-1.D0)
  319. c GR(1,8,L)=G1X*N2Y
  320. GR(1,8,L)=(2.D0*alfai*X(L)-(1.D0 + alfai))
  321. & *alfa2*Y(L)*(Y(L)-1.D0)
  322. c GR(2,8,L)=N1X*G2Y
  323. GR(2,8,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
  324. & *alfa2*(2.D0*Y(L)-1.D0)
  325. c GR(1,9,L)=G2X*N2Y
  326. GR(1,9,L)=alfa2*Y(L)*(Y(L)-1.D0)
  327. & *alfa2*(2.D0*X(L)-1.D0)
  328. c GR(2,9,L)=N2X*G2Y
  329. GR(2,9,L)=alfa2*X(L)*(X(L)-1.D0)
  330. & *alfa2*(2.D0*Y(L)-1.D0)
  331. C
  332. c GR(1,1,L)=(Y(L)-1.D0)*(2.D0*Y(L)-1.D0)*(4.D0*X(L)-3.D0)
  333. c GR(2,1,L)=(X(L)-1.D0)*(2.D0*X(L)-1.D0)*(4.D0*Y(L)-3.D0)
  334. c GR(1,2,L)=-4.D0*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)*(Y(L)-1.D0)
  335. c GR(2,2,L)=-4.D0*X(L)*(X(L)-1.D0)*(4.D0*Y(L)-3.D0)
  336. c GR(1,3,L)=(2.D0*Y(L)-1.D0)*(Y(L)-1.D0)*(4.D0*X(L)-1.D0)
  337. c GR(2,3,L)=X(L)*(2.D0*X(L)-1.D0)*(4.D0*Y(L)-3.D0)
  338. c GR(1,4,L)=-4.D0*Y(L)*(Y(L)-1.D0)*(4.D0*X(L)-1.D0)
  339. c GR(2,4,L)=-4.D0*X(L)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)
  340. c GR(1,5,L)=Y(L)*(2.D0*Y(L)-1.D0)*(4.D0*X(L)-1.D0)
  341. c GR(2,5,L)=X(L)*(2.D0*X(L)-1.D0)*(4.D0*Y(L)-1.D0)
  342. c GR(1,6,L)=-4.D0*Y(L)*(2.D0*Y(L)-1.D0)*(2.D0*X(L)-1.D0)
  343. c GR(2,6,L)=-4.D0*X(L)*(X(L)-1.D0)*(4.D0*Y(L)-1.D0)
  344. c GR(1,7,L)=Y(L)*(2.D0*Y(L)-1.D0)*(4.D0*X(L)-3.D0)
  345. c GR(2,7,L)=(2.D0*X(L)-1.D0)*(X(L)-1.D0)*(4.D0*Y(L)-1.D0)
  346. c GR(1,8,L)=-4.D0*Y(L)*(Y(L)-1.D0)*(4.D0*X(L)-3.D0)
  347. c GR(2,8,L)=-4.D0*(X(L)-1.D0)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)
  348. c GR(1,9,L)=16.D0*Y(L)*(Y(L)-1.D0)*(2.D0*X(L)-1.D0)
  349. c GR(2,9,L)=16.D0*X(L)*(X(L)-1.D0)*(2.D0*Y(L)-1.D0)
  350. C
  351. 1 CONTINUE
  352. C
  353. C
  354. C WRITE(6,101)
  355. C WRITE(6,1002)FM
  356. C WRITE(6,1002)GM
  357. C WRITE(6,1002)FN
  358. C WRITE(6,1002)GR
  359. C WRITE(6,101)
  360.  
  361. C write(6,*)' RET PB902 '
  362. RETURN
  363. 1002 FORMAT(10(1X,1PD11.4))
  364. 1001 FORMAT(20(1X,I5))
  365. 101 FORMAT(1X,'... SUBPB902 ... FM,GM,FN,GR ',9(10H..........)/)
  366. C
  367. END
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  
  377.  

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