Télécharger caljbr.eso

Retour à la liste

Numérotation des lignes :

  1. C CALJBR SOURCE MAGN 08/10/10 21:15:00 6185
  2. SUBROUTINE CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,
  3. *NES,ND,NP,NPG,IAXI,AIRE,AJ,SGN)
  4. C***********************************************************************
  5. C Ce Sp calcule le gradient des fonctions de formes pour un element
  6. C courant a partir des gradients de l'element de reference et des
  7. C coordonnees des points de l'element
  8. C Il calcule le produit PgSq poids d'integration x element d'aire et
  9. C le Jacobien
  10. C Pour les elements coques ou filaires les normales et tangentes
  11. C sont donnees a la place du Jacobien
  12. C les fonctions de forme FN sont donnees en entree pour calculer
  13. C les elements d'aire en axisymetrique
  14. C
  15. C Entree :
  16. C NES Dimension d'espace de l'element de reference
  17. C (different de ND pour les elements coques ou filaires)
  18. C ND Dimension d'espace
  19. C NP Nombre de noeuds de l'element
  20. C NPG Nombre de points d'integration
  21. C IAXI =0 en 3D et en 2D PLAN
  22. C IAXI NE 0 en axi symetrique axe de symetrie x
  23. C FN(NP,NPG) Valeurs des fonctions de forme aux points d'integration
  24. C GR(NES,NP,NPG) Valeurs des gradients dans l'element de reference
  25. C PG(NPG) Poids d'integration
  26. C XYZ(ND,NP) Coordonnees des noeuds de l'element
  27. C
  28. C Sortie
  29. C HR(NES,NP,NPG) Valeurs des gradients dans l'element courant
  30. C PGSQ(NPG) poids d'integration x element d'aire
  31. C ATTENTION en axi : poids d'integration x element d'aire x 2 pi R
  32. C RPG(NPG) Rayon des points d'integration
  33. C AIRE Aire de l'element
  34. C AJ(ND,ND,NPG) Valeurs du jacobien pour chaque point d'integration
  35. C SGN Pour les elements massif =1 si meme orientation element de Ref
  36. C SGN Pour les elements massif =-1 si orientation opposee
  37. C SGN 0 pour les autres types d'elements
  38. C
  39. C
  40. C CE SP TRAITE LES ELEMENTS VOLUMIQUES, SURFACIQUES et FILAIRES
  41. C DANS LES CAS 2D (PLAN et AXI) ET 3D
  42. C
  43. C Dans AJ on stoke l'inverse du Jacobien (AJ=1/J) si l'element est
  44. C volumique
  45. C Sinon pour les elements coques on stoke tangentes et normales
  46. C AJ=|tx ty| |tx ty tz|
  47. C |nx ny| ou |ux uy uz|
  48. C |nx ny nz|
  49. C
  50. C
  51. C CALCUL INTERMEDIAIRE DE L'ELEMENT D'AIRE SQ=DET(J)
  52. C
  53. C***********************************************************************
  54. IMPLICIT INTEGER(I-N)
  55. IMPLICIT REAL*8 (A-H,O-Z)
  56. C
  57. REAL*8 FN(NP,NPG),GR(NES,NP,NPG),HR(ND,NP,NPG)
  58. REAL*8 PG(NPG),XYZ(ND,NP),PGSQ(NPG),RPG(NPG),AJ(ND,ND,NPG)
  59. REAL*8 AG(3,3,25)
  60. -INC CCREEL
  61. C
  62. C***
  63. SGN=0.D0
  64. CALL INITD(RPG,NPG,1.D0)
  65.  
  66. IF(NES.EQ.1.AND.ND.EQ.2)THEN
  67.  
  68. IF(NPG.GT.25)CALL ARRET(0)
  69. AIRE=0.D0
  70. DO 110 L=1,NPG
  71. AJX=0.D0
  72. AJY=0.D0
  73. DO 111 I=1,NP
  74. AJX=AJX+GR(1,I,L)*XYZ(1,I)
  75. AJY=AJY+GR(1,I,L)*XYZ(2,I)
  76. 111 CONTINUE
  77. AJN=(AJX*AJX+AJY*AJY)**0.5D0
  78. C write(6,*)' AJX,AJY=',AJX,AJY,AJN
  79. DET=AJX*AJX+AJY*AJY
  80.  
  81. AJ(1,1,L)=AJX/AJN
  82. AJ(2,1,L)=AJY/AJN
  83. AJ(1,2,L)=-AJY/AJN
  84. AJ(2,2,L)=AJX/AJN
  85.  
  86. AG(1,1,L)=AJX/DET
  87. AG(2,1,L)=-AJY/DET
  88. AG(1,2,L)=AJY/DET
  89. AG(2,2,L)=AJX/DET
  90.  
  91. PGSQ(L)=PG(L)*AJN
  92. AIRE=AIRE+PGSQ(L)
  93. 110 CONTINUE
  94. C write(6,*)' AJ ds CALJBR',AIRE
  95. C write(6,1002)AJ
  96. DO 131 L=1,NPG
  97. DO 131 I=1,NP
  98. DO 131 N=1,ND
  99. U=0.D0
  100. DO 132 M=1,NES
  101. 132 U=U+AG(M,N,L)*GR(M,I,L)
  102. 131 HR(N,I,L)=U
  103. C write(6,*)'GR'
  104. C write(6,1002)gr
  105. C write(6,*)'HR'
  106. C write(6,1002)hr
  107.  
  108. ELSEIF(NES.EQ.1.AND.ND.EQ.3)THEN
  109.  
  110. IF(NPG.GT.25)CALL ARRET(0)
  111. AIRE=0.D0
  112. DO 310 L=1,NPG
  113. AJX=0.D0
  114. AJY=0.D0
  115. AJZ=0.D0
  116. BJX=0.D0
  117. BJY=0.D0
  118. BJZ=0.D0
  119. DO 311 I=1,NP
  120. AJX=AJX+GR(1,I,L)*XYZ(1,I)
  121. AJY=AJY+GR(1,I,L)*XYZ(2,I)
  122. AJZ=AJZ+GR(1,I,L)*XYZ(3,I)
  123. C BJX=BJX+GR(2,I,L)*XYZ(1,I)
  124. C BJY=BJY+GR(2,I,L)*XYZ(2,I)
  125. C BJZ=BJZ+GR(2,I,L)*XYZ(3,I)
  126. 311 CONTINUE
  127. C write(6,*)ajx,ajy,ajz
  128. IF(ABS(AJX).NE.0.D0.AND.ABS(AJY).NE.0.D0)THEN
  129. BJX=AJY
  130. BJY=-AJX
  131. BJZ=AJZ
  132. ELSE
  133. BJX=1.D0
  134. BJY=1.D0
  135. BJZ=0.D0
  136. ENDIF
  137.  
  138. XB=AJY*BJZ-AJZ*BJY
  139. YB=AJZ*BJX-AJX*BJZ
  140. ZB=AJX*BJY-AJY*BJX
  141.  
  142. AJN=(XB*XB+YB*YB*ZB*ZB)**0.5D0+1.D-5
  143. C write(6,*)' AJN=',AJN
  144. PGSQ(L)=PG(L)*AJN
  145. AIRE=AIRE+PGSQ(L)
  146.  
  147. AJ(1,1,L)=AJX/AJN
  148. AJ(2,1,L)=AJY/AJN
  149. AJ(3,1,L)=AJZ/AJN
  150. AJ(1,2,L)=BJX/AJN
  151. AJ(2,2,L)=BJY/AJN
  152. AJ(3,2,L)=BJZ/AJN
  153. AJ(1,3,L)=XB/AJN
  154. AJ(2,3,L)=YB/AJN
  155. AJ(3,3,L)=ZB/AJN
  156.  
  157. AG(1,1,L)=AJX
  158. AG(2,1,L)=AJY
  159. AG(3,1,L)=AJZ
  160. AG(1,2,L)=BJX
  161. AG(2,2,L)=BJY
  162. AG(3,2,L)=BJZ
  163. AG(1,3,L)=XB
  164. AG(2,3,L)=YB
  165. AG(3,3,L)=ZB
  166.  
  167. D11=AG(2,2,L)*AG(3,3,L)-AG(3,2,L)*AG(2,3,L)
  168. D12=AG(1,2,L)*AG(3,3,L)-AG(3,2,L)*AG(1,3,L)
  169. D13=AG(1,2,L)*AG(2,3,L)-AG(2,2,L)*AG(1,3,L)
  170. D21=AG(2,1,L)*AG(3,3,L)-AG(3,1,L)*AG(2,3,L)
  171. D22=AG(1,1,L)*AG(3,3,L)-AG(3,1,L)*AG(1,3,L)
  172. D23=AG(1,1,L)*AG(2,3,L)-AG(2,1,L)*AG(1,3,L)
  173. D31=AG(2,1,L)*AG(3,2,L)-AG(3,1,L)*AG(2,2,L)
  174. D32=AG(1,1,L)*AG(3,2,L)-AG(3,1,L)*AG(1,2,L)
  175. D33=AG(1,1,L)*AG(2,2,L)-AG(2,1,L)*AG(1,2,L)
  176. VINT=AG(1,1,L)*D11-AG(1,2,L)*D21+AG(1,3,L)*D31
  177. RVINT=1.D0/VINT
  178. AG(1,1,L)= RVINT*D11
  179. AG(1,2,L)=-RVINT*D12
  180. AG(1,3,L)= RVINT*D13
  181. AG(2,1,L)=-RVINT*D21
  182. AG(2,2,L)= RVINT*D22
  183. AG(2,3,L)=-RVINT*D23
  184. AG(3,1,L)= RVINT*D31
  185. AG(3,2,L)=-RVINT*D32
  186. AG(3,3,L)= RVINT*D33
  187.  
  188. 310 CONTINUE
  189. C write(6,*)' AJ ds CALJBR'
  190. C write(6,1002)AJ
  191. DO 331 L=1,NPG
  192. DO 331 I=1,NP
  193. DO 331 N=1,ND
  194. U=0.D0
  195. DO 332 M=1,NES
  196. 332 U=U+AG(M,N,L)*GR(M,I,L)
  197. 331 HR(N,I,L)=U
  198. C write(6,*)'GR'
  199. C write(6,1002)gr
  200. C write(6,*)'HR'
  201. C write(6,1002)hr
  202.  
  203. ELSEIF(NES.EQ.2.AND.ND.EQ.3)THEN
  204. C write(6,*)' Je passe ICI'
  205.  
  206. AIRE=0.D0
  207. DO 210 L=1,NPG
  208. AJX=0.D0
  209. AJY=0.D0
  210. AJZ=0.D0
  211. BJX=0.D0
  212. BJY=0.D0
  213. BJZ=0.D0
  214. DO 211 I=1,NP
  215. AJX=AJX+GR(1,I,L)*XYZ(1,I)
  216. AJY=AJY+GR(1,I,L)*XYZ(2,I)
  217. AJZ=AJZ+GR(1,I,L)*XYZ(3,I)
  218. BJX=BJX+GR(2,I,L)*XYZ(1,I)
  219. BJY=BJY+GR(2,I,L)*XYZ(2,I)
  220. BJZ=BJZ+GR(2,I,L)*XYZ(3,I)
  221. 211 CONTINUE
  222.  
  223. XB=AJY*BJZ-AJZ*BJY
  224. YB=AJZ*BJX-AJX*BJZ
  225. ZB=AJX*BJY-AJY*BJX
  226.  
  227. AJN=(XB*XB+YB*YB+ZB*ZB)**0.5D0
  228.  
  229. PGSQ(L)=PG(L)*AJN
  230. AIRE=AIRE+PGSQ(L)
  231.  
  232. AJ(1,1,L)=AJX/AJN
  233. AJ(2,1,L)=AJY/AJN
  234. AJ(3,1,L)=AJZ/AJN
  235. AJ(1,2,L)=BJX/AJN
  236. AJ(2,2,L)=BJY/AJN
  237. AJ(3,2,L)=BJZ/AJN
  238. AJ(1,3,L)=XB/AJN
  239. AJ(2,3,L)=YB/AJN
  240. AJ(3,3,L)=ZB/AJN
  241.  
  242. AG(1,1,L)=AJX
  243. AG(2,1,L)=AJY
  244. AG(3,1,L)=AJZ
  245. AG(1,2,L)=BJX
  246. AG(2,2,L)=BJY
  247. AG(3,2,L)=BJZ
  248. AG(1,3,L)=XB
  249. AG(2,3,L)=YB
  250. AG(3,3,L)=ZB
  251.  
  252. D11=AG(2,2,L)*AG(3,3,L)-AG(3,2,L)*AG(2,3,L)
  253. D12=AG(1,2,L)*AG(3,3,L)-AG(3,2,L)*AG(1,3,L)
  254. D13=AG(1,2,L)*AG(2,3,L)-AG(2,2,L)*AG(1,3,L)
  255. D21=AG(2,1,L)*AG(3,3,L)-AG(3,1,L)*AG(2,3,L)
  256. D22=AG(1,1,L)*AG(3,3,L)-AG(3,1,L)*AG(1,3,L)
  257. D23=AG(1,1,L)*AG(2,3,L)-AG(2,1,L)*AG(1,3,L)
  258. D31=AG(2,1,L)*AG(3,2,L)-AG(3,1,L)*AG(2,2,L)
  259. D32=AG(1,1,L)*AG(3,2,L)-AG(3,1,L)*AG(1,2,L)
  260. D33=AG(1,1,L)*AG(2,2,L)-AG(2,1,L)*AG(1,2,L)
  261. VINT=AG(1,1,L)*D11-AG(1,2,L)*D21+AG(1,3,L)*D31
  262. RVINT=1.D0/VINT
  263. AG(1,1,L)= RVINT*D11
  264. AG(1,2,L)=-RVINT*D12
  265. AG(1,3,L)= RVINT*D13
  266. AG(2,1,L)=-RVINT*D21
  267. AG(2,2,L)= RVINT*D22
  268. AG(2,3,L)=-RVINT*D23
  269. AG(3,1,L)= RVINT*D31
  270. AG(3,2,L)=-RVINT*D32
  271. AG(3,3,L)= RVINT*D33
  272.  
  273. 210 CONTINUE
  274. C
  275. DO 231 L=1,NPG
  276. DO 231 I=1,NP
  277. DO 231 N=1,ND
  278. U=0.D0
  279. DO 232 M=1,NES
  280. 232 U=U+AG(M,N,L)*GR(M,I,L)
  281. 231 HR(N,I,L)=U
  282. C write(6,*)'GR'
  283. C write(6,1002)gr
  284. C write(6,*)'HR'
  285. C write(6,1002)hr
  286.  
  287. ELSE
  288.  
  289. DO 10 L=1,NPG
  290. DO 10 M=1,ND
  291. DO 10 N=1,ND
  292. AJT=0.D0
  293. DO 11 I=1,NP
  294. AJT=AJT+GR(M,I,L)*XYZ(N,I)
  295. 11 CONTINUE
  296. AJ(N,M,L)=AJT
  297. 10 CONTINUE
  298.  
  299. C
  300. DO 20 L=1,NPG
  301. IF(ND.EQ.1)THEN
  302. VINT=AJ(1,1,L)
  303. C VINT=ABS(VINT)
  304. AJ(1,1,L)=1.D0/VINT
  305. ELSEIF(ND.EQ.2)THEN
  306. VINT=AJ(1,1,L)*AJ(2,2,L)-AJ(1,2,L)*AJ(2,1,L)
  307. C VINT=ABS(VINT)
  308. RVINT=1.D0/VINT
  309. D11=AJ(2,2,L)
  310. D12=AJ(1,2,L)
  311. D21=AJ(2,1,L)
  312. D22=AJ(1,1,L)
  313. AJ(1,1,L)= RVINT*D11
  314. AJ(1,2,L)=-RVINT*D12
  315. AJ(2,1,L)=-RVINT*D21
  316. AJ(2,2,L)= RVINT*D22
  317. ELSEIF(ND.EQ.3)THEN
  318. D11=AJ(2,2,L)*AJ(3,3,L)-AJ(3,2,L)*AJ(2,3,L)
  319. D12=AJ(1,2,L)*AJ(3,3,L)-AJ(3,2,L)*AJ(1,3,L)
  320. D13=AJ(1,2,L)*AJ(2,3,L)-AJ(2,2,L)*AJ(1,3,L)
  321. D21=AJ(2,1,L)*AJ(3,3,L)-AJ(3,1,L)*AJ(2,3,L)
  322. D22=AJ(1,1,L)*AJ(3,3,L)-AJ(3,1,L)*AJ(1,3,L)
  323. D23=AJ(1,1,L)*AJ(2,3,L)-AJ(2,1,L)*AJ(1,3,L)
  324. D31=AJ(2,1,L)*AJ(3,2,L)-AJ(3,1,L)*AJ(2,2,L)
  325. D32=AJ(1,1,L)*AJ(3,2,L)-AJ(3,1,L)*AJ(1,2,L)
  326. D33=AJ(1,1,L)*AJ(2,2,L)-AJ(2,1,L)*AJ(1,2,L)
  327. VINT=AJ(1,1,L)*D11-AJ(1,2,L)*D21+AJ(1,3,L)*D31
  328. C VINT=ABS(VINT)
  329. RVINT=1.D0/VINT
  330. AJ(1,1,L)= RVINT*D11
  331. AJ(1,2,L)=-RVINT*D12
  332. AJ(1,3,L)= RVINT*D13
  333. AJ(2,1,L)=-RVINT*D21
  334. AJ(2,2,L)= RVINT*D22
  335. AJ(2,3,L)=-RVINT*D23
  336. AJ(3,1,L)= RVINT*D31
  337. AJ(3,2,L)=-RVINT*D32
  338. AJ(3,3,L)= RVINT*D33
  339. ELSE
  340. VINT=0.D0
  341. ENDIF
  342.  
  343. PGSQ(L)=VINT
  344. C*** CALL DLAIN(ND,ND,AJ(1,1,L),KAUX,B,IER)
  345. C
  346. 20 CONTINUE
  347. C
  348. DO 31 L=1,NPG
  349. DO 31 I=1,NP
  350. DO 31 N=1,ND
  351. U=0.D0
  352. DO 32 M=1,ND
  353. 32 U=U+AJ(M,N,L)*GR(M,I,L)
  354. 31 HR(N,I,L)=U
  355. U=0.D0
  356. V=0.D0
  357. DO 4 L=1,NPG
  358. W=ABS(PGSQ(L))
  359. U=U+PG(L)*W
  360. V=V+PG(L)*PGSQ(L)
  361. 4 PGSQ(L)=PG(L)*W
  362. AIRE=U
  363. SGN=SIGN(1.D0,V)
  364.  
  365. ENDIF
  366.  
  367. IF(IAXI.EQ.0)RETURN
  368. C
  369. DO 45 L=1,NPG
  370. RPGT =0.D0
  371. DO 46 I=1,NP
  372. RPGT=RPGT+XYZ(1,I)*FN(I,L)
  373. 46 CONTINUE
  374. RPG(L)=RPGT
  375. 45 CONTINUE
  376. C
  377. DEUPI=2.D0*XPI
  378. U=0.D0
  379. DO 47 L=1,NPG
  380. C? PGSQ(L)=PGSQ(L)*DEUPI*RPG(L)
  381. U=U+PGSQ(L)*DEUPI*RPG(L)
  382. 47 CONTINUE
  383. AIRE=U
  384. C
  385. RETURN
  386. 1002 FORMAT(10(1X,1PE11.4))
  387. 1001 FORMAT(20(1X,I5))
  388. END
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  

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