Télécharger bbcal3.eso

Retour à la liste

Numérotation des lignes :

  1. C BBCAL3 SOURCE KICH 18/12/05 21:15:01 10024
  2. SUBROUTINE BBCAL3(I2B,I2GAU,IDIM,NBPGAU,
  3. 1 POIGAU,QSIGAU,ETAGAU,DZEGAU,MELE,MFR,NBNN,LRE,IFOUR,NGRA,NHRM,
  4. 2 A,BB,XE,SHPTOT,SHP,BU,BGR,P)
  5. c SUBROUTINE BBCAL2(I2B,I2GAU,IDIM,NBPGAU,IVACAR,
  6. c 1 POIGAU,QSIGAU,ETAGAU,DZEGAU,MELE,MFR,NBNN,LRE,IFOUR,nstrs,NHRM,
  7. c 2 A,BB,XE,SHPTOT,SHP,BGENE,XDPGE,YDPGE,P)
  8. *******************************************************************
  9. * projection de dN/dx dans l espace de dimension inferieure
  10. c=======================================================================
  11. c calcul des composantes de la matrice b-barre-dlatation
  12. c - en 2d: elements ict3,icq4,ict6
  13. c - en 3d: elements ....
  14. c calcul des coefficients de modification de la matrice b-barre-dilatation
  15. c - en 2d: elements icq8
  16. c - en 3d: elements ....
  17. c coefficients {a}={a1,a2,a3}
  18. c=======================================================================
  19. c input
  20. c xe = coordonnees de l'element
  21. c mele = numero de l'element dans nomtp
  22. c nbnn = nombre de noeuds
  23. c idim = dimension espace
  24. c nbpgau = nombre de points de gauss pour la rigidite
  25. c poigau = poids des points de gauss
  26. c qsigau,etagau,dzegau = coordonnees des points de gauss
  27. c ngra = nombre de composantes du gradient
  28. c lre = nombre de colonnes de la matrice b
  29. c ifour = ifour dans ccoptio
  30. c nhrm = numero du mode de fourier
  31. c bu : tableau utile BGRMAS
  32. c
  33. c output
  34. c tableau des coefficients de modification: a(4,*)
  35. c tableau des composantes bbarre-dilatation: bb(3,*)
  36. c=======================================================================
  37. IMPLICIT INTEGER(I-N)
  38. IMPLICIT REAL*8 (A-H,O-Z)
  39. -INC CCREEL
  40. -INC SMCHAML
  41.  
  42. PARAMETER (XZer=0., XUn=1., X1s2=0.5, X1s4 = 0.25,
  43. & X2s3 = 0.666666666666666666666666666666666666666667,
  44. & X1s3 = 0.333333333333333333333333333333333333333333,
  45. & X1s6 = 0.166666666666666666666666666666666666666667)
  46.  
  47. DIMENSION XE(3,*), POIGAU(*),QSIGAU(*),ETAGAU(*),DZEGAU(*),
  48. & SHP(6,*),SHPTOT(6,NBNN,*),
  49. & A(4,*),BB(3,*),BU(2,*),BGR(ngra,*)
  50.  
  51. DIMENSION P(4,4),XNUM1(8),XNUM2(8),XNUM3(8)
  52.  
  53. SEGMENT MPTVAL
  54. INTEGER IPOS(NS) ,NSOF(NS)
  55. INTEGER IVAL(NCOSOU)
  56. CHARACTER*16 TYVAL(NCOSOU)
  57. ENDSEGMENT
  58.  
  59. if (MELE.EQ.72.or.MELE.GE.76) then
  60. CALL ZERO(A,4,60)
  61. CALL ZERO(P,4,4)
  62. endif
  63.  
  64. IFR = IFOUR+4
  65. if (ifr.le.4.or.ifr.ge.7) then
  66. KINC = 2
  67. else
  68. KINC = 3
  69. endif
  70.  
  71.  
  72. IF ( MELE.LT.69 .OR.
  73. &(MELE.GT.78.AND.MELE.NE.273.AND.MELE.NE.274)) GOTO 999
  74. C=======================================================================
  75. C= Elements incompressibles implementes :
  76. C=======================================================================
  77. C= NOM : ICT3, ICQ4, ICT6, ICQ8, ICC8, ICT4, ICP6, IC20, IC10, IC15,
  78. C= MELE : 69 , 70 , 71 , 72 , 73 , 74 , 75 , 76 , 77 , 78 ,
  79. C= NOM : ICY5, IC13
  80. C= MELE : 273, 274
  81. CC=======================================================================
  82. CALL ZERO(BB,3,60)
  83. IF(MELE.EQ.273.OR.MELE.EQ.274) GOTO 100
  84. GOTO ( 10, 20, 30, 40, 100,100,100,100,100,100) , (MELE-68)
  85. GOTO 999
  86.  
  87. 3 FORMAT (4(A,I1,A,I1,A,F10.4,2X))
  88. 4 FORMAT (3(A,I1,A,I1,A,F10.4,2X))
  89.  
  90. C=======================================================================
  91. C========== Elements MASSIFS INCOMPRESSIBLES BIDIMENSIONNELS ===========
  92. C=======================================================================
  93.  
  94. *-----------------------------------------------------------------------
  95. *----------------------------- Element ICT3 ----------------------------
  96. *-----------------------------------------------------------------------
  97. 10 CONTINUE
  98. * write(*,*) 'Element ICT3',ifour
  99. * write(*,*) 'Ecriture de la matrice xe[ ]'
  100. * write(*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,3)
  101. IFR = ifour+4
  102. GOTO (666,101,101,103,666,999),IFR
  103. GOTO 999
  104.  
  105. *--------contraintes planes ou déformations planes (ifour=-2,-1)
  106. 101 CONTINUE
  107. *--------donnee des composantes de bb-dil
  108. r_z = - XE(1,1)*XE(2,3) - XE(1,2)*XE(2,1) + XE(1,2)*XE(2,3)
  109. & + XE(1,1)*XE(2,2) + XE(1,3)*XE(2,1) - XE(1,3)*XE(2,2)
  110. * AIRE = X1s2 * r_z
  111. r_z = XUn / r_z
  112. BB(1,1) = (XE(2,2)-XE(2,3)) * r_z
  113. BB(2,2) = (XE(1,3)-XE(1,2)) * r_z
  114. BB(1,3) = (XE(2,3)-XE(2,1)) * r_z
  115. BB(2,4) = (XE(1,1)-XE(1,3)) * r_z
  116. BB(1,5) = (XE(2,1)-XE(2,2)) * r_z
  117. BB(2,6) = (XE(1,2)-XE(1,1)) * r_z
  118. * write (*,*) 'calcul de aire pour bbdil'
  119. * write (*,*) 'aire = ',AIRE
  120. * write (*,*) 'ecriture de la matrice bb'
  121. * write (*,3) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn)
  122. GOTO 666
  123. *--------axisymétrique (ifour=0)
  124. 103 CONTINUE
  125. *--------donnee des composantes de bb-dil
  126. r_z = - XE(1,1)*XE(2,3) - XE(1,2)*XE(2,1) + XE(1,2)*XE(2,3)
  127. & + XE(1,1)*XE(2,2) + XE(1,3)*XE(2,1) - XE(1,3)*XE(2,2)
  128. * AIRE = X1s6 * r_z * (XE(1,1)+XE(1,2)+XE(1,3))
  129. r_z = XUn / r_z
  130. BB(1,1) = r_z * (XE(2,2)-XE(2,3))
  131. BB(2,2) = r_z * (XE(1,3)-XE(1,2))
  132. BB(1,3) = r_z * (XE(2,3)-XE(2,1))
  133. BB(2,4) = r_z * (XE(1,1)-XE(1,3))
  134. BB(1,5) = r_z * (XE(2,1)-XE(2,2))
  135. BB(2,6) = r_z * (XE(1,2)-XE(1,1))
  136. *-- Integration analytique des termes BB(3,1),BB(3,3),BB(3,5)
  137. BB(3,1) = XUn / (XE(1,1)+XE(1,2)+XE(1,3))
  138. BB(3,3) = BB(3,1)
  139. BB(3,5) = BB(3,1)
  140. * write(*,*) 'calcul de aire pour bbdil'
  141. * write(*,*) 'aire = ',AIRE
  142. * write(*,*) 'ecriture de la matrice bb'
  143. * write(*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn)
  144. GOTO 666
  145.  
  146. *-----------------------------------------------------------------------
  147. *----------------------------- Element ICQ4 ----------------------------
  148. *-----------------------------------------------------------------------
  149. 20 CONTINUE
  150. * write(*,*) 'Element ICQ4',ifour
  151. * write(*,*) 'Ecriture de la matrice xe[ ]'
  152. * write(*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,4)
  153. IFR = ifour+4
  154. GOTO (666,201,201,203,666,999),IFR
  155. GOTO 999
  156. *--------contraintes planes ou déformations planes (ifour=-2,-1)
  157. 201 CONTINUE
  158. *--------donnée des composantes de bb-dil
  159. r_z = - XE(1,4)*XE(2,3) - XE(1,2)*XE(2,1) + XE(1,1)*XE(2,2)
  160. & - XE(1,1)*XE(2,4) - XE(1,3)*XE(2,2) + XE(1,2)*XE(2,3)
  161. & + XE(1,3)*XE(2,4) + XE(1,4)*XE(2,1)
  162. * AIRE = X1s2 * r_z
  163. r_z = XUn / r_z
  164. BB(1,1) = r_z * (XE(2,2)-XE(2,4))
  165. BB(2,2) = r_z * (XE(1,4)-XE(1,2))
  166. BB(1,3) = r_z * (XE(2,3)-XE(2,1))
  167. BB(2,4) = r_z * (XE(1,1)-XE(1,3))
  168. BB(1,5) = r_z * (XE(2,4)-XE(2,2))
  169. BB(2,6) = r_z * (XE(1,2)-XE(1,4))
  170. BB(1,7) = r_z * (XE(2,1)-XE(2,3))
  171. BB(2,8) = r_z * (XE(1,3)-XE(1,1))
  172. * write (*,*) 'calcul de aire pour bbdil'
  173. * write (*,*) 'aire = ',aire
  174. * write (*,*) 'ecriture de la matrice bb'
  175. * write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn)
  176. GOTO 666
  177. *
  178. *--------axisymétriques(ifour=0)
  179. 203 CONTINUE
  180. r_z = XE(1,1)*( XE(1,1)*XE(2,2)-XE(1,1)*XE(2,4)-XE(1,4)*XE(2,4)
  181. & +XE(1,4)*XE(2,1)-XE(1,2)*XE(2,1)+XE(1,2)*XE(2,2))
  182. & + XE(1,2)*(-XE(1,2)*XE(2,1)+XE(1,2)*XE(2,3)+XE(1,3)*XE(2,3)
  183. & -XE(1,3)*XE(2,2))
  184. & + XE(1,3)*( XE(1,3)*XE(2,4)+XE(1,4)*XE(2,4)-XE(1,4)*XE(2,3)
  185. & -XE(1,3)*XE(2,2))
  186. & + XE(1,4)*XE(1,4)*(XE(2,1)-XE(2,3))
  187. * AIRE = x1s6 * r_z
  188. r_z = xUn / r_z
  189. BB(1,1) = ( XE(2,2)*(XE(1,2)+XE(1,1))-XE(2,4)*(XE(1,1)+XE(1,4))
  190. & +X1s2*( XE(1,3)*(XE(2,2)-XE(2,4))
  191. & +XE(1,4)*(XE(2,2)+XE(2,3))
  192. & -XE(1,2)*(XE(2,3)+XE(2,4))) ) * r_z
  193. BB(2,2) = ( XE(1,4)*(XE(1,1)+XE(1,4))
  194. & -XE(1,2)*(XE(1,1)+XE(1,2)) ) * r_z
  195. BB(1,3) = ( XE(2,3)*(XE(1,2)+XE(1,3))-XE(2,1)*(XE(1,2)+XE(1,1))
  196. & +X1s2*( XE(1,4)*(XE(2,3)-XE(2,1))
  197. & +XE(1,1)*(XE(2,3)+XE(2,4))
  198. & -XE(1,3)*(XE(2,4)+XE(2,1))) ) * r_z
  199. BB(2,4) = ( XE(1,1)*(XE(1,2)+XE(1,1))
  200. & -XE(1,3)*(XE(1,3)+XE(1,2)) ) * r_z
  201. BB(1,5) = ( XE(2,4)*(XE(1,3)+XE(1,4))-XE(2,2)*(XE(1,3)+XE(1,2))
  202. & +X1s2*( XE(1,1)*(XE(2,4)-XE(2,2))
  203. & +XE(1,2)*(XE(2,4)+XE(2,1))
  204. & -XE(1,4)*(XE(2,1)+XE(2,2))) ) * r_z
  205. BB(2,6) = ( XE(1,2)*(XE(1,2)+XE(1,3))
  206. & -XE(1,4)*(XE(1,4)+XE(1,3)) ) * r_z
  207. BB(1,7) = ( XE(2,1)*(XE(1,4)+XE(1,1))-XE(2,3)*(XE(1,4)+XE(1,3))
  208. & +X1s2*( XE(1,2)*(XE(2,1)-XE(2,3))
  209. & +XE(1,3)*(XE(2,1)+XE(2,2))
  210. & -XE(1,1)*(XE(2,2)+XE(2,3))) ) * r_z
  211. BB(2,8) = ( XE(1,3)*(XE(1,3)+XE(1,4))
  212. & -XE(1,1)*(XE(1,1)+XE(1,4)) ) * r_z
  213. *- Integration analytique des termes BB(3,1),BB(3,3),BB(3,5),BB(3,7)
  214. BB(3,1) = ( XE(1,1)*(XE(2,2)-XE(2,4))+XE(2,1)*(XE(1,4)-XE(1,2))
  215. & +X1s2*( XE(1,2)*(XE(2,3)+XE(2,4))
  216. & +XE(1,3)*(XE(2,4)-XE(2,2))
  217. & -XE(1,4)*(XE(2,3)+XE(2,2))) ) * r_z
  218. BB(3,3) = ( XE(1,2)*(XE(2,3)-XE(2,1))+XE(2,2)*(XE(1,1)-XE(1,3))
  219. & +X1s2*( XE(1,3)*(XE(2,1)+XE(2,4))
  220. & +XE(1,4)*(XE(2,1)-XE(2,3))
  221. & -XE(1,1)*(XE(2,3)+XE(2,4))) ) * r_z
  222. BB(3,5) = ( XE(1,3)*(XE(2,4)-XE(2,2))+XE(2,3)*(XE(1,2)-XE(1,4))
  223. & +X1s2*( XE(1,4)*(XE(2,1)+XE(2,2))
  224. & +XE(1,1)*(XE(2,2)-XE(2,4))
  225. & -XE(1,2)*(XE(2,4)+XE(2,1))) ) * r_z
  226. BB(3,7) = ( XE(1,4)*(XE(2,1)-XE(2,3))+XE(2,4)*(XE(1,3)-XE(1,1))
  227. & +X1s2*( XE(1,1)*(XE(2,2)+XE(2,3))
  228. & +XE(1,2)*(XE(2,3)-XE(2,1))
  229. & -XE(1,3)*(XE(2,1)+XE(2,2))) ) * r_z
  230. * write (*,*) 'calcul de aire pour bbdil'
  231. * write (*,*) 'aire = ',aire
  232. * write (*,*) 'ecriture de la matrice bb'
  233. * write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn)
  234. GOTO 666
  235.  
  236. *-----------------------------------------------------------------------
  237. *----------------------------- Element ICT6 ----------------------------
  238. *-----------------------------------------------------------------------
  239. 30 CONTINUE
  240. * write(*,*) 'Element ICT6',ifour
  241. * write (*,*) 'ecriture de la matrice xe[ ]'
  242. * write (*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,6)
  243. IFR = ifour+4
  244. GOTO (666,301,301,303,666,999),IFR
  245. GOTO 999
  246. *
  247. *--------contraintes planes ou déformations planes(ifour=-2,-1)
  248. 301 CONTINUE
  249. *--------donnée des composantes de bb-dil
  250. r_z = ( XE(1,1)*(XE(2,2)-XE(2,6)) + XE(1,2)*(XE(2,3)-XE(2,1))
  251. & + XE(1,3)*(XE(2,4)-XE(2,2)) + XE(1,4)*(XE(2,5)-XE(2,3))
  252. & + XE(1,5)*(XE(2,6)-XE(2,4)) + XE(1,6)*(XE(2,1)-XE(2,5)) )
  253. & + X1s4 * ( XE(1,1)*(XE(2,5)-XE(2,3))
  254. & + XE(1,3)*(XE(2,1)-XE(2,5))
  255. & + XE(1,5)*(XE(2,3)-XE(2,1)) )
  256. AIRE = X2s3 * r_z
  257. r_z = xUn / r_z
  258. BB(1, 1) = ( (XE(2,2)-XE(2,6)) + X1s4 * (XE(2,5)-XE(2,3)) ) * r_z
  259. BB(2, 2) = ( (XE(1,6)-XE(1,2)) + X1s4 * (XE(1,3)-XE(1,5)) ) * r_z
  260. BB(1, 3) = (XE(2,3)-XE(2,1)) * r_z
  261. BB(2, 4) = (XE(1,1)-XE(1,3)) * r_z
  262. BB(1, 5) = ( (XE(2,4)-XE(2,2)) + X1s4 * (XE(2,1)-XE(2,5)) ) * r_z
  263. BB(2, 6) = ( (XE(1,2)-XE(1,4)) + X1s4 * (XE(1,5)-XE(1,1)) ) * r_z
  264. BB(1, 7) = (XE(2,5)-XE(2,3)) * r_z
  265. BB(2, 8) = (XE(1,3)-XE(1,5)) * r_z
  266. BB(1, 9) = ( (XE(2,6)-XE(2,4)) + X1s4 * (XE(2,3)-XE(2,1)) ) * r_z
  267. BB(2,10) = ( (XE(1,4)-XE(1,6)) + X1s4 * (XE(1,1)-XE(1,3)) ) * r_z
  268. BB(1,11) = (XE(2,1)-XE(2,5)) * r_z
  269. BB(2,12) = (XE(1,5)-XE(1,1)) * r_z
  270.  
  271. * write (*,*) 'calcul de aire pour bbdil'
  272. * write (*,*) 'aire = ',aire
  273. * write (*,*) 'ecriture de la matrice bb'
  274. * write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn)
  275. GOTO 666
  276. *
  277. *--------axisymétriques (ifour=0)
  278. 303 CONTINUE
  279. *--------donnée des composantes de bb-dil
  280. *-----intégration numérique des composantes bb-dil------
  281. AIRE = XZer
  282. DO i = 1,nbnn
  283. XNUM1(i) = XZer
  284. XNUM2(i) = XZer
  285. XNUM3(i) = XZer
  286. ENDDO
  287. DO 3031 IGAU = 1,NBPGAU
  288. DO i = 1,nbnn
  289. SHP(1,i) = SHPTOT(1,i,IGAU)
  290. SHP(2,i) = SHPTOT(2,i,IGAU)
  291. SHP(3,i) = SHPTOT(3,i,IGAU)
  292. ENDDO
  293. CALL DISTRR(XE,SHP,nbnn,RR)
  294. CALL DISTR2(XE,SHP,nbnn,RR2)
  295. CALL DISTR3(XE,SHP,nbnn,RR3)
  296. CALL DISTZ2(XE,SHP,nbnn,ZZ2)
  297. CALL DISTZ3(XE,SHP,nbnn,ZZ3)
  298. CALL JACOBI(XE,SHP,2,nbnn,DJAC)
  299. r_z = POIGAU(IGAU) * RR
  300. r_zz2 = r_z * ZZ2
  301. r_zz3 = r_z * ZZ3
  302. r_rr2 = r_z * RR2
  303. r_rr3 = r_z * RR3
  304. r_z = POIGAU(IGAU) * DJAC
  305. AIRE = AIRE + (r_z * RR)
  306. DO i = 1,nbnn
  307. XNUM1(i) = XNUM1(i) + ( SHPTOT(2,i,IGAU)*r_zz3
  308. & - SHPTOT(3,i,IGAU)*r_zz2 )
  309. XNUM2(i) = XNUM2(i) + ( SHPTOT(3,i,IGAU)*r_rr2
  310. & - SHPTOT(2,i,IGAU)*r_rr3 )
  311. XNUM3(i) = XNUM3(i) + ( SHPTOT(1,i,IGAU)*r_z )
  312. ENDDO
  313. 3031 CONTINUE
  314. r_z = XUn / AIRE
  315. DO i = 1,nbnn
  316. k = 2*i
  317. BB(1,k-1) = XNUM1(i) * r_z
  318. BB(2,k) = XNUM2(i) * r_z
  319. BB(3,k-1) = XNUM3(i) * r_z
  320. ENDDO
  321. * write (*,*) 'calcul de aire pour bbdil'
  322. * write (*,*) 'aire = ',aire
  323. * write (*,*) 'ecriture de la matrice bb'
  324. * write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn)
  325. GOTO 666
  326.  
  327. *-----------------------------------------------------------------------
  328. *----------------------------- Element ICQ8 ----------------------------
  329. *-----------------------------------------------------------------------
  330. 40 CONTINUE
  331. * write(*,*) 'Element ICQ8',ifour
  332. * write (*,*) 'Ecriture de la matrice xe[ ]'
  333. * write (*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,6)
  334.  
  335. IFR = ifour+4
  336. GOTO (666,401,401,403,666,999),IFR
  337. GOTO 666
  338.  
  339. *--------contraintes planes ou déformations planes (ifour=-2,-1)
  340. 401 CONTINUE
  341. *--------donnée de la matrice symétrique p(3,3) (1 fois par élément)
  342. P(1,1) = X2s3 * ( XE(1,1)*(XE(2,2)-XE(2,8))
  343. 1 + XE(1,2)*(XE(2,3)-XE(2,1))
  344. 1 +XE(1,3)*(XE(2,4)-XE(2,2))+XE(1,4)*(XE(2,5)-XE(2,3))
  345. 2 +XE(1,5)*(XE(2,6)-XE(2,4))+XE(1,6)*(XE(2,7)-XE(2,5))
  346. 3 +XE(1,7)*(XE(2,8)-XE(2,6))+XE(1,8)*(XE(2,1)-XE(2,7)) )
  347. 4 + X1s6 * ( (XE(1,1)-XE(1,5))*(XE(2,7)-XE(2,3))
  348. 5 +(XE(1,3)-XE(1,7))*(XE(2,1)-XE(2,5)) )
  349. P(1,2) = 5.D0/9.D0*(-XE(2,2)*(XE(1,1)+XE(1,3))
  350. 1 +XE(1,2)*(XE(2,1)+XE(2,3))+XE(1,6)*(-XE(2,5)-XE(2,7))
  351. 1 +XE(2,6)*(XE(1,5)+XE(1,7)))
  352. 2 +4.D0/9.D0*(XE(1,8)*(XE(2,7)+XE(2,2)-XE(2,6)-XE(2,1))
  353. 2 +XE(2,8)*(-XE(1,7)-XE(1,2)+XE(1,6)+XE(1,1))
  354. 2 +XE(1,4)*(XE(2,2)-XE(2,6)+XE(2,5)-XE(2,3))
  355. 2 +XE(2,4)*(-XE(1,2)+XE(1,6)-XE(1,5)+XE(1,3)))
  356. 3 +7.D0/45.D0*(XE(1,2)*(XE(2,5)+XE(2,7))
  357. 3 -XE(2,2)*(XE(1,5)+XE(1,7))-XE(1,6)*(XE(2,1)+XE(2,3))
  358. 3 +XE(2,6)*(XE(1,1)+XE(1,3)))
  359. 4 +8.D0/15.D0*(XE(1,6)*XE(2,2)-XE(1,2)*XE(2,6))
  360. 5 +7.D0/90.D0*(XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5)
  361. 5 +XE(1,5)*XE(2,3)-XE(1,1)*XE(2,7))
  362. 6 +1.D0/30.D0*(-XE(1,7)*XE(2,3)+XE(1,3)*XE(2,7)
  363. 6 -XE(1,5)*XE(2,1)+XE(1,1)*XE(2,5))
  364. P(1,3) = 5.D0/9.D0*(-XE(1,8)*(XE(2,7)+XE(2,1))
  365. 1 +XE(2,8)*(XE(1,7)+XE(1,1))
  366. 1 +XE(1,4)*(XE(2,3)+XE(2,5))
  367. 1 -XE(2,4)*(XE(1,3)+XE(1,5)))
  368. 2 +4.D0/9.D0*((XE(1,8)-XE(1,4))*(XE(2,2)+XE(2,6))
  369. 2 +(-XE(2,8)+XE(2,4))*(XE(1,2)+XE(1,6))
  370. 2 +XE(1,2)*(XE(2,1)-XE(2,3))+XE(2,2)*(-XE(1,1)+XE(1,3))
  371. 2 +XE(1,6)*(XE(2,7)-XE(2,5))+XE(2,6)*(-XE(1,7)+XE(1,5)))
  372. 3 +7.D0/45.D0*(-XE(1,8)*(XE(2,3)+XE(2,5))
  373. 3 +XE(2,8)*(XE(1,3)+XE(1,5))
  374. 3 +XE(1,4)*(XE(2,1)+XE(2,7))
  375. 3 -XE(2,4)*(XE(1,1)+XE(1,7)))
  376. 4 +1.D0/30.D0*(-XE(1,7)*XE(2,3)+XE(1,3)*XE(2,7)
  377. 4 +XE(1,5)*XE(2,1)-XE(1,1)*XE(2,5))
  378. 5 +7.D0/90.D0*(XE(1,7)*XE(2,5)-XE(1,3)*XE(2,1)
  379. 5 -XE(1,5)*XE(2,7)+XE(1,1)*XE(2,3))
  380. 6 +8.D0/15.D0*(XE(1,8)*XE(2,4)-XE(1,4)*XE(2,8))
  381. P(1,4) = XZer
  382.  
  383. P(2,1) = P(1,2)
  384. P(2,2) = 16.D0/45.D0*(XE(2,6)*(XE(1,5)-XE(1,7))
  385. 1 +XE(1,6)*(XE(2,7)-XE(2,5))+XE(1,2)*(XE(2,3)-XE(2,1))
  386. 1 +XE(2,2)*(XE(1,1)-XE(1,3)))
  387. 2 +14.D0/45.D0*(XE(2,4)*(XE(1,3)-XE(1,5))
  388. 2 +XE(1,4)*(XE(2,5)-XE(2,3))+XE(1,8)*(XE(2,1)-XE(2,7))
  389. 2 +XE(2,8)*(XE(1,7)-XE(1,1)))
  390. 3 +8.D0/45.D0*(XE(1,4)*(XE(2,2)-XE(2,6))
  391. 3 +XE(1,8)*(XE(2,6)-XE(2,2))+XE(1,2)*(XE(2,8)-XE(2,4))
  392. 3 +XE(1,6)*(XE(2,4)-XE(2,8)))
  393. 4 +2.D0/45.D0*(XE(2,2)*(XE(1,5)-XE(1,7))
  394. 4 +XE(1,2)*(XE(2,7)-XE(2,5)) +XE(1,6)*(XE(2,3)-XE(2,1))
  395. 4 +XE(2,6)*(XE(1,1)-XE(1,3)))
  396. 5 +4.D0/45.D0*(XE(2,4)*(XE(1,1)-XE(1,7))
  397. 5 +XE(1,4)*(XE(2,7)-XE(2,1)) +XE(1,8)*(XE(2,3)-XE(2,5))
  398. 5 +XE(2,8)*(XE(1,5)-XE(1,3)))
  399. 6 +17.D0/90.D0*(XE(1,7)*XE(2,5)+XE(1,3)*XE(2,1)
  400. 6 -XE(1,5)*XE(2,7)-XE(1,1)*XE(2,3))
  401. 7 +1.D0/90.D0*(-XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5)
  402. 7 +XE(1,5)*XE(2,3)+XE(1,1)*XE(2,7))
  403. P(2,3) = 1.D0/3.D0*(XE(1,5)*(XE(2,6)-XE(2,4))
  404. 1 +XE(1,8)*(XE(2,7)+XE(2,1))+XE(1,7)*(XE(2,6)-XE(2,8))
  405. 1 -XE(1,2)*(XE(2,1)+XE(2,3))+XE(1,1)*(XE(2,2)-XE(2,8))
  406. 1 +XE(1,4)*(XE(2,5)+XE(2,3))-XE(1,6)*(XE(2,7)+XE(2,5))
  407. 1 +XE(1,3)*(XE(2,2)-XE(2,4)))
  408. 2 +4.D0/9.D0*(-XE(1,4)*(XE(2,2)+XE(2,6))
  409. 2 -XE(1,8)*(XE(2,6)+XE(2,2))+XE(1,2)*(XE(2,8)+XE(2,4))
  410. 2 +XE(1,6)*(XE(2,4)+XE(2,8)))
  411. 3 +1.D0/9.D0*(XE(1,7)*(XE(2,2)-XE(2,4))
  412. 3 +XE(1,8)*(XE(2,3)+XE(2,5))
  413. 3 +XE(1,5)*(XE(2,2)-XE(2,8))+XE(1,3)*(XE(2,6)-XE(2,8))
  414. 3 -XE(1,2)*(XE(2,5)+XE(2,7))+XE(1,1)*(XE(2,6)-XE(2,4))
  415. 3 +XE(1,4)*(XE(2,1)+XE(2,7))-XE(1,6)*(XE(2,1)+XE(2,3)))
  416. P(2,4) = XZer
  417.  
  418. P(3,1) = P(1,3)
  419. P(3,2) = P(2,3)
  420. P(3,3) = 14.D0/45.D0*(XE(1,6)*(XE(2,7)-XE(2,5))
  421. 1 +XE(2,6)*(XE(1,5)-XE(1,7))+XE(1,2)*(XE(2,3)-XE(2,1))
  422. 1 +XE(2,2)*(XE(1,1)-XE(1,3)))
  423. 2 +16.D0/45.D0*(XE(1,8)*(XE(2,1)-XE(2,7))
  424. 2 +XE(2,8)*(XE(1,7)-XE(1,1))+XE(1,4)*(XE(2,5)-XE(2,3))
  425. 2 +XE(2,4)*(XE(1,3)-XE(1,5)))
  426. 3 +8.D0/45.D0*(XE(1,4)*(XE(2,2)-XE(2,6))
  427. 3 +XE(1,8)*(XE(2,6)-XE(2,2))+XE(1,2)*(XE(2,8)-XE(2,4))
  428. 3 +XE(1,6)*(XE(2,4)-XE(2,8)))
  429. 4 +2.D0/45.D0*(XE(2,4)*(XE(1,7)-XE(1,1))
  430. 4 +XE(1,8)*(XE(2,5)-XE(2,3))+XE(2,8)*(-XE(1,5)+XE(1,3))
  431. 4 +XE(1,4)*(XE(2,1)-XE(2,7)))
  432. 5 +4.D0/45.D0*(XE(2,2)*(XE(1,7)-XE(1,5))
  433. 5 +XE(1,2)*(XE(2,5)-XE(2,7))+XE(1,6)*(XE(2,1)-XE(2,3))
  434. 5 +XE(2,6)*(XE(1,3)-XE(1,1)))
  435. 6 +1.D0/90.D0*(-XE(1,5)*XE(2,7)+XE(1,3)*XE(2,1)
  436. 6 +XE(1,7)*XE(2,5)-XE(1,1)*XE(2,3))
  437. 7 +17.D0/90.D0*(-XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5)
  438. 7 +XE(1,5)*XE(2,3)+XE(1,1)*XE(2,7))
  439. P(3,4) = XZer
  440.  
  441. P(4,1) = XZer
  442. P(4,2) = XZer
  443. P(4,3) = XZer
  444. P(4,4) = XZer
  445.  
  446. * write (*,*) 'ecriture de la matrice p[ ]'
  447. * write (*,4) (('p(',i,',',j,')=',p(i,j),i=1,3),j=1,3)
  448.  
  449. *-----------------donnée des vecteurs {t}={t1,t2,t3}--------------------
  450. A(1,1) = 2.D0/3.D0*(XE(2,2)-XE(2,8))
  451. 1 +1.D0/6.D0*(XE(2,7)-XE(2,3))
  452. A(2,1) = 1.D0/9.D0*(4.D0*XE(2,8)-5.D0*XE(2,2))
  453. 1 +7.D0/90.D0*(-XE(2,7)+2.D0*XE(2,6))
  454. 1 +1.D0/30.D0*XE(2,5)
  455. A(3,1) = 1.D0/9.D0*(-4.D0*XE(2,2)+5.D0*XE(2,8))
  456. 1 +7.D0/90.D0*(XE(2,3)-2.D0*XE(2,4))
  457. 1 -1.D0/30.D0*XE(2,5)
  458.  
  459. A(1,2) = 2.D0/3.D0*(XE(1,8)-XE(1,2))
  460. 1 +1.D0/6.D0*(XE(1,3)-XE(1,7))
  461. A(2,2) = 1.D0/9.D0*(-4.D0*XE(1,8)+5.D0*XE(1,2))
  462. 1 +7.D0/90.D0*(XE(1,7)-2.D0*XE(1,6))
  463. 1 -1.D0/30.D0*XE(1,5)
  464. A(3,2) = 1.D0/9.D0*(4.D0*XE(1,2)-5.D0*XE(1,8))
  465. 1 +7.D0/90.D0*(-XE(1,3)+2.D0*XE(1,4))
  466. 1 +1.D0/30.D0*XE(1,5)
  467.  
  468. A(1,3) = 2.D0/3.D0*(-XE(2,1)+XE(2,3))
  469. A(2,3) = 5.D0/9.D0*(XE(2,1)+XE(2,3))
  470. 1 +7.D0/45.D0*(XE(2,5)+XE(2,7))
  471. 1 -4.D0/9.D0*(XE(2,4)+XE(2,8))
  472. 1 -8.D0/15.D0*XE(2,6)
  473. A(3,3) = 4.D0/9.D0*(XE(2,1)-XE(2,3)+XE(2,4)-XE(2,8))
  474.  
  475. A(1,4) = 2.D0/3.D0*(-XE(1,3)+XE(1,1))
  476. A(2,4) = -5.D0/9.D0*(XE(1,1)+XE(1,3))
  477. 1 -7.D0/45.D0*(XE(1,5)+XE(1,7))
  478. 1 +4.D0/9.D0*(XE(1,4)+XE(1,8))
  479. 1 +8.D0/15.D0*XE(1,6)
  480. A(3,4) = 4.D0/9.D0*(-XE(1,1)+XE(1,3)-XE(1,4)+XE(1,8))
  481.  
  482. A(1,5) = 2.D0/3.D0*(XE(2,4)-XE(2,2))
  483. 1 +1.D0/6.D0*(-XE(2,5)+XE(2,1))
  484. A(2,5) = 1.D0/9.D0*(4.D0*XE(2,4)-5.D0*XE(2,2))
  485. 1 +7.D0/90.D0*(-XE(2,5)+2.D0*XE(2,6))
  486. 1 +1.D0/30.D0*XE(2,7)
  487. A(3,5) = 1.D0/9.D0*(4.D0*XE(2,2)-5.D0*XE(2,4))
  488. 1 +7.D0/90.D0*(-XE(2,1)+2.D0*XE(2,8))
  489. 1 +1.D0/30.D0*XE(2,7)
  490.  
  491. A(1,6) = 2.D0/3.D0*(-XE(1,4)+XE(1,2))
  492. 1 +1.D0/6.D0*(XE(1,5)-XE(1,1))
  493. A(2,6) = 1.D0/9.D0*(-4.D0*XE(1,4)+5.D0*XE(1,2))
  494. 1 +7.D0/90.D0*(XE(1,5)-2.D0*XE(1,6))
  495. 1 -1.D0/30.D0*XE(1,7)
  496. A(3,6) = 1.D0/9.D0*(-4.D0*XE(1,2)+5.D0*XE(1,4))
  497. 1 +7.D0/90.D0*(XE(1,1)-2.D0*XE(1,8))
  498. 1 -1.D0/30.D0*XE(1,7)
  499.  
  500. A(1,7) = 2.D0/3.D0*(-XE(2,3)+XE(2,5))
  501. A(2,7) = 4.D0/9.D0*(XE(2,5)-XE(2,3)+XE(2,2)-XE(2,6))
  502. A(3,7) = 5.D0/9.D0*(XE(2,3)+XE(2,5))
  503. 1 +7.D0/45.D0*(XE(2,1)+XE(2,7))
  504. 1 -4.D0/9.D0*(XE(2,2)+XE(2,6))
  505. 1 -8.D0/15.D0*XE(2,8)
  506.  
  507. A(1,8) = 2.D0/3.D0*(XE(1,3)-XE(1,5))
  508. A(2,8) = 4.D0/9.D0*(-XE(1,5)+XE(1,3)-XE(1,2)+XE(1,6))
  509. A(3,8) = -5.D0/9.D0*(XE(1,3)+XE(1,5))
  510. 1 -7.D0/45.D0*(XE(1,1)+XE(1,7))
  511. 1 +4.D0/9.D0*(XE(1,2)+XE(1,6))
  512. 1 +8.D0/15.D0*XE(1,8)
  513.  
  514. A(1,9) = 2.D0/3.D0*(-XE(2,4)+XE(2,6))
  515. 1 +1.D0/6.D0*(XE(2,3)-XE(2,7))
  516. A(2,9) = 1.D0/9.D0*(-4.D0*XE(2,4)+5.D0*XE(2,6))
  517. 1 +7.D0/90.D0*(XE(2,3)-2.D0*XE(2,2))
  518. 1 -1.D0/30.D0*XE(2,1)
  519. A(3,9) = 1.D0/9.D0*(4.D0*XE(2,6)-5.D0*XE(2,4))
  520. 1 +7.D0/90.D0*(-XE(2,7)+2.D0*XE(2,8))
  521. 1 +1.D0/30.D0*XE(2,1)
  522.  
  523. A(1,10) = 2.D0/3.D0*(XE(1,4)-XE(1,6))
  524. 1 +1.D0/6.D0*(-XE(1,3)+XE(1,7))
  525. A(2,10) = 1.D0/9.D0*(4.D0*XE(1,4)-5.D0*XE(1,6))
  526. 1 +7.D0/90.D0*(-XE(1,3)+2.D0*XE(1,2))
  527. 1 +1.D0/30.D0*XE(1,1)
  528. A(3,10) = 1.D0/9.D0*(-4.D0*XE(1,6)+5.D0*XE(1,4))
  529. 1 +7.D0/90.D0*(XE(1,7)-2.D0*XE(1,8))
  530. 1 -1.D0/30.D0*XE(1,1)
  531.  
  532. A(1,11) = 2.D0/3.D0*(-XE(2,5)+XE(2,7))
  533. A(2,11) = -5.D0/9.D0*(XE(2,7)+XE(2,5))
  534. 1 -7.D0/45.D0*(XE(2,1)+XE(2,3))
  535. 1 +4.D0/9.D0*(XE(2,4)+XE(2,8))
  536. 1 +8.D0/15.D0*XE(2,2)
  537. A(3,11) = 4.D0/9.D0*(XE(2,4)-XE(2,5)-XE(2,8)+XE(2,7))
  538.  
  539. A(1,12) = 2.D0/3.D0*(XE(1,5)-XE(1,7))
  540. A(2,12) = 5.D0/9.D0*(XE(1,7)+XE(1,5))
  541. 1 +7.D0/45.D0*(XE(1,1)+XE(1,3))
  542. 1 -4.D0/9.D0*(XE(1,4)+XE(1,8))
  543. 1 -8.D0/15.D0*XE(1,2)
  544. A(3,12) = 4.D0/9.D0*(-XE(1,4)+XE(1,5)+XE(1,8)-XE(1,7))
  545.  
  546. A(1,13) = 2.D0/3.D0*(-XE(2,6)+XE(2,8))
  547. 1 +1.D0/6.D0*(XE(2,5)-XE(2,1))
  548. A(2,13) = 1.D0/9.D0*(-4.D0*XE(2,8)+5.D0*XE(2,6))
  549. 1 +7.D0/90.D0*(XE(2,1)-2.D0*XE(2,2))
  550. 1 -1.D0/30.D0*XE(2,3)
  551. A(3,13) = 1.D0/9.D0*(-4.D0*XE(2,6)+5.D0*XE(2,8))
  552. 1 +7.D0/90.D0*(XE(2,5)-2.D0*XE(2,4))
  553. 1 -1.D0/30.D0*XE(2,3)
  554.  
  555. A(1,14) = 2.D0/3.D0*(XE(1,6)-XE(1,8))
  556. 1 +1.D0/6.D0*(-XE(1,5)+XE(1,1))
  557. A(2,14) = 1.D0/9.D0*(4.D0*XE(1,8)-5.D0*XE(1,6))
  558. 1 +7.D0/90.D0*(-XE(1,1)+2.D0*XE(1,2))
  559. 1 +1.D0/30.D0*XE(1,3)
  560. A(3,14) = 1.D0/9.D0*(4.D0*XE(1,6)-5.D0*XE(1,8))
  561. 1 +7.D0/90.D0*(-XE(1,5)+2.D0*XE(1,4))
  562. 1 +1.D0/30.D0*XE(1,3)
  563.  
  564. A(1,15) = 2.D0/3.D0*(-XE(2,7)+XE(2,1))
  565. A(2,15) = 4.D0/9.D0*(-XE(2,6)+XE(2,7)+XE(2,2)-XE(2,1))
  566. A(3,15) = -5.D0/9.D0*(XE(2,1)+XE(2,7))
  567. 1 -7.D0/45.D0*(XE(2,3)+XE(2,5))
  568. 1 +4.D0/9.D0*(XE(2,2)+XE(2,6))
  569. 1 +8.D0/15.D0*XE(2,4)
  570.  
  571. A(1,16) = 2.D0/3.D0*(XE(1,7)-XE(1,1))
  572. A(2,16) = 4.D0/9.D0*(XE(1,6)-XE(1,7)-XE(1,2)+XE(1,1))
  573. A(3,16) = 5.D0/9.D0*(XE(1,1)+XE(1,7))
  574. 1 +7.D0/45.D0*(XE(1,3)+XE(1,5))
  575. 1 -4.D0/9.D0*(XE(1,2)+XE(1,6))
  576. 1 -8.D0/15.D0*XE(1,4)
  577.  
  578. * write (*,*) 'ecriture de la matrice t[ ]'
  579. * write (*,4) (('t(',i,',',j,')=',a(i,j),i=1,3),j=1,16)
  580. * if (IB.EQ.34.and.mele.eq.72) then
  581. * write(*,*)'bbcalc'
  582. * do j = 1,16
  583. * write(*,*) (a(i,j),i=1,3)
  584. * enddo
  585. * do j = 1,4
  586. * write(6,*) (p(i,j),i=1,4)
  587. * enddo
  588. * endif
  589. *--------------résolution matricielle de [p]{a}={t}---------------------
  590. * write (*,*) 'Appel a GAUSSJ'
  591. N = 3
  592. NP = 4
  593. M = 16
  594. MP = 60
  595. CALL GAUSSJ(P,N,NP,A,M,MP)
  596. GOTO 666
  597.  
  598. *--------axisymétriques (ifour=0)
  599. 403 CONTINUE
  600. *--------donnée de la matrice symétrique p(3,3) (1 fois par élément)
  601. CALL ZERO(P,4,4)
  602. DO IGAU = 1,NBPGAU
  603. DO I = 1,nbnn
  604. SHP(1,I) = SHPTOT(1,I,IGAU)
  605. SHP(2,I) = SHPTOT(2,I,IGAU)
  606. SHP(3,I) = SHPTOT(3,I,IGAU)
  607. ENDDO
  608. CALL DISTRR(XE,SHP,nbnn,RR)
  609. CALL JACOBI(XE,SHP,2,nbnn,DJAC)
  610. r_z = POIGAU(IGAU)*RR*DJAC
  611. P(1,1) = P(1,1) + r_z
  612. P(2,1) = P(2,1) + (r_z * QSIGAU(IGAU))
  613. P(3,1) = P(3,1) + (r_z * ETAGAU(IGAU))
  614. P(2,2) = P(2,2) + (r_z * QSIGAU(IGAU) * QSIGAU(IGAU))
  615. P(3,2) = P(3,2) + (r_z * QSIGAU(IGAU) * ETAGAU(IGAU))
  616. P(3,3) = P(3,3) + (r_z * ETAGAU(IGAU) * ETAGAU(IGAU))
  617. ENDDO
  618. P(1,2) = P(2,1)
  619. P(1,3) = P(3,1)
  620. P(2,3) = P(3,2)
  621. * write (*,*) 'ecriture de la matrice p[ ]'
  622. * write (*,4) (('p(',i,',',j,')=',p(i,j),i=1,3),j=1,3)
  623. *-----------calcul des vecteurs {t}={t1,t2,t3}--intégration numérique---
  624. CALL ZERO(A,4,24)
  625. DO 4031 IGAU = 1,NBPGAU
  626. DO I = 1,nbnn
  627. SHP(1,I) = SHPTOT(1,I,IGAU)
  628. SHP(2,I) = SHPTOT(2,I,IGAU)
  629. SHP(3,I) = SHPTOT(3,I,IGAU)
  630. ENDDO
  631. CALL DISTRR(XE,SHP,nbnn,RR)
  632. CALL DISTR2(XE,SHP,nbnn,RR2)
  633. CALL DISTR3(XE,SHP,nbnn,RR3)
  634. CALL DISTZ2(XE,SHP,nbnn,ZZ2)
  635. CALL DISTZ3(XE,SHP,nbnn,ZZ3)
  636. CALL JACOBI(XE,SHP,2,nbnn,DJAC)
  637. K = 0
  638. DO I = 1,nbnn
  639. K = K+1
  640. r_z = POIGAU(IGAU)*RR
  641. & *(SHPTOT(2,I,IGAU)*ZZ3-SHPTOT(3,I,IGAU)*ZZ2)
  642. A(1,K) = A(1,K) + r_z
  643. A(2,K) = A(2,K) + r_z * QSIGAU(IGAU)
  644. A(3,K) = A(3,K) + r_z * ETAGAU(IGAU)
  645. K = K+1
  646. r_z = POIGAU(IGAU)*RR
  647. & *(-SHPTOT(2,I,IGAU)*RR3+SHPTOT(3,I,IGAU)*RR2)
  648. A(1,K) = A(1,K) + r_z
  649. A(2,K) = A(2,K) + r_z * QSIGAU(IGAU)
  650. A(3,K) = A(3,K) + r_z * ETAGAU(IGAU)
  651. K = K+1
  652. r_z = POIGAU(IGAU)*SHPTOT(1,I,IGAU)*DJAC
  653. A(1,K) = A(1,K) + r_z
  654. A(2,K) = A(2,K) + r_z * QSIGAU(IGAU)
  655. A(3,K) = A(3,K) + r_z * ETAGAU(IGAU)
  656. ENDDO
  657. 4031 CONTINUE
  658. * write(*,*) 'Ecriture de la matrice t[ ]'
  659. * write(*,4) (('t(',i,',',j,')=',a(i,j),i=1,3),j=1,24)
  660. *------------résolution matricielle de [p]{a}={t}-----------------
  661. N = 3
  662. NP = 4
  663. M = 24
  664. MP = 60
  665. * write (*,*) 'appel a GAUSSJ'
  666. CALL GAUSSJ(P,N,NP,A,M,MP)
  667. * write(*,*) 'Ecriture de la matrice a[ ]'
  668. * write(*,4) (('a(',i,',',j,')=',a(i,j),i=1,3),j=1,24)
  669. GOTO 666
  670.  
  671. C=======================================================================
  672. C========== Elements MASSIFS INCOMPRESSIBLES TRIDIMENSIONNELS ==========
  673. C=======================================================================
  674.  
  675. 100 CONTINUE
  676. IF (ifour.NE.2) GOTO 999
  677.  
  678. DIM3=1.d0
  679. FACAR=1.D0
  680.  
  681. C= 3.12.3 - Boucle sur les points d'integration
  682. ESTEL=XZero
  683. DO iGau=1,NBPGAU
  684. * IBMN=MIN(IB,MELVA1.VELCHE(/2))
  685. * IGMN=MIN(iGau,MELVA1.VELCHE(/1))
  686. * FACSCA=MELVA1.VELCHE(IGMN,IBMN)
  687. FACSCA = 1.D0
  688.  
  689. * BGR au point IGAU
  690. * if (I2B.EQ.1.and.I2GAU.EQ.1) then
  691. * write(6,*) IGAU,NBPGAU,POIGAU(1),QSIGAU(1),ETAGAU(1),DZEGAU(1)
  692. * endif
  693. CALL BGRMAS(iGau,MELE,NBNN,LRE,IFOUR,NGRA,NHRM,XE,
  694. . XZer,SHPTOT,SHP,BU,BGR,DJAC,0)
  695.  
  696. C= 3.12.3.9 - Calcul de la composante integree en ce point de Gauss
  697. 80 DJAC=ABS(DJAC)*POIGAU(iGau)*FACAR*DIM3
  698. ESTEL=ESTEL+FACSCA*DJAC
  699.  
  700. if (MELE.EQ.72.or.(MELE.ge.76.and.MELE.le.78).or.
  701. &MELE.EQ.274) then
  702. * second membre T
  703. K = 1
  704. jj= 0
  705. DO i = 1,NBNN
  706. jj = jj + 1
  707. A(1,jj) = A(1,jj) + BGR(1,K)*DJAC
  708. A(2,jj) = A(2,jj) + BGR(1,K)*qsigau(igau)*DJAC
  709. A(3,jj) = A(3,jj) + BGR(1,K)*etagau(igau)*DJAC
  710. if(mele.ge.76)
  711. & A(4,jj) = A(4,jj) + BGR(1,K)*dzegau(igau)*DJAC
  712. jj = jj + 1
  713. A(1,jj) = A(1,jj) + BGR(2,K+1)*DJAC
  714. A(2,jj) = A(2,jj) + BGR(2,K+1)*qsigau(igau)*DJAC
  715. A(3,jj) = A(3,jj) + BGR(2,K+1)*etagau(igau)*DJAC
  716. if(mele.ge.76) then
  717. A(4,jj) = A(4,jj) + BGR(2,K+1)*dzegau(igau)*DJAC
  718. jj = jj + 1
  719. A(1,jj) = A(1,jj) + BGR(3,K+2)*DJAC
  720. A(2,jj) = A(2,jj) + BGR(3,K+2)*qsigau(igau)*DJAC
  721. A(3,jj) = A(3,jj) + BGR(3,K+2)*etagau(igau)*DJAC
  722. A(4,jj) = A(4,jj) + BGR(3,K+2)*dzegau(igau)*DJAC
  723. endif
  724. K = K + KINC
  725. ENDDO
  726. * matrice M
  727. P(1,1) = P(1,1) + DJAC
  728. P(1,2) = P(1,2) + qsigau(igau)*DJAC
  729. P(1,3) = P(1,3) + etagau(igau)*DJAC
  730. P(2,2) = P(2,2) + qsigau(igau)*qsigau(igau)*DJAC
  731. P(2,3) = P(2,3) + qsigau(igau)*etagau(igau)*DJAC
  732.  
  733. P(3,3) = P(3,3) + etagau(igau)*etagau(igau)*DJAC
  734.  
  735. if(mele.ge.76) then
  736. P(1,4) = P(1,4) + dzegau(igau)*DJAC
  737.  
  738. P(2,4) = P(2,4) + qsigau(igau)*dzegau(igau)*DJAC
  739.  
  740. P(3,4) = P(3,4) + etagau(igau)*dzegau(igau)*DJAC
  741.  
  742. P(4,4) = P(4,4) + dzegau(igau)*dzegau(igau)*DJAC
  743. endif
  744.  
  745. else if(mele.le.71) then
  746.  
  747. K = 1
  748. DO i = 1,NBNN
  749. bb(1,K) = bb(1,K) + BGR(1,K)*DJAC
  750. bb(2,K+1) = bb(2,K+1) + BGR(2,K+1)*DJAC
  751. bb(3,K) = bb(3,K) + BGR(3,K)*DJAC
  752. K = K + KINC
  753. ENDDO
  754.  
  755. elseif((mele.ge.73.and.mele.le.75)
  756. &.or.mele.eq.273) then
  757.  
  758. K = 1
  759. DO i = 1,NBNN
  760. bb(1,K) = bb(1,K) + BGR(1,K)*DJAC
  761. bb(2,K+1) = bb(2,K+1) + BGR(2,K+1)*DJAC
  762. bb(3,K+2) = bb(3,K+2) + BGR(3,K+2)*DJAC
  763. K = K + KINC
  764. ENDDO
  765.  
  766. endif
  767.  
  768.  
  769. * if (I2B.EQ.34) then
  770. * write(6,*) 'djac',igau,DJAC,POIGAU(IGAU)
  771. * endif
  772.  
  773. ENDDO
  774.  
  775. if (MELE.EQ.72) then
  776. P(1,4) = XZer
  777.  
  778. P(2,1) = P(1,2)
  779. P(2,4) = XZer
  780. P(3,1) = P(1,3)
  781. P(3,2) = P(2,3)
  782. P(3,4) = XZer
  783.  
  784. P(4,1) = XZer
  785. P(4,2) = XZer
  786. P(4,3) = XZer
  787. P(4,4) = XZer
  788.  
  789. elseif(mele.ge.76) then
  790. P(2,1) = P(1,2)
  791. P(3,1) = P(1,3)
  792. P(3,2) = P(2,3)
  793. P(4,1) = P(1,4)
  794. P(4,2) = P(2,4)
  795. P(4,3) = P(3,4)
  796.  
  797. endif
  798.  
  799.  
  800. *
  801. ** synthese
  802. *
  803.  
  804. if (MELE.EQ.72) then
  805. * if (I2B.EQ.34) then
  806. * do j = 1,16
  807. * write(*,*) (a(i,j),i=1,3)
  808. * enddo
  809. * do j = 1,4
  810. * write(6,*) (p(i,j),i=1,4)
  811. * enddo
  812. * endif
  813. N = 3
  814. NP = 4
  815. M = 16
  816. MP = 60
  817. * write (*,*) 'appel a GAUSSJ'
  818. CALL GAUSSJ(P,N,NP,A,M,MP)
  819. * write(*,*) 'Ecriture de la matrice a[ ]'
  820. * write(*,4) (('a(',i,',',j,')=',a(i,j),i=1,3),j=1,24)
  821.  
  822. elseif (mele.ge.76.and.mele.ne.273) then
  823. N = 4
  824. NP = 4
  825. M = NBNN*3
  826. MP = 60
  827. * write (*,*) 'appel a GAUSSJ'
  828. CALL GAUSSJ(P,N,NP,A,M,MP)
  829.  
  830. elseif(mele.le.71) then
  831.  
  832. K = 1
  833. DO i = 1,NBNN
  834. bb(1,K) = bb(1,K)/ESTEL
  835. bb(2,K+1) = bb(2,K+1)/ESTEL
  836. bb(3,K) = bb(3,K)/ESTEL
  837. * if (I2B.EQ.34) then
  838. * write(6,*) 'bb',i,K, bb(1,K),bb(2,K+1),bb(3,K)
  839. * endif
  840. K = K + KINC
  841. ENDDO
  842.  
  843. elseif((mele.ge.73.and.mele.le.75).or.
  844. &mele.eq.273) then
  845.  
  846. K = 1
  847. DO i = 1,NBNN
  848. bb(1,K) = bb(1,K)/ESTEL
  849. bb(2,K+1) = bb(2,K+1)/ESTEL
  850. bb(3,K+2) = bb(3,K+2)/ESTEL
  851. * if (I2B.EQ.34) then
  852. * write(6,*) 'bb',i,K, bb(1,K),bb(2,K+1),bb(3,K+2)
  853. * endif
  854. K = K + KINC
  855. ENDDO
  856.  
  857. endif
  858.  
  859. RETURN
  860.  
  861. C=======================================================================
  862. C======================= ERREUR : Cas non prevus =======================
  863. C=======================================================================
  864. 999 CONTINUE
  865. CALL ERREUR(5)
  866.  
  867. 666 CONTINUE
  868. * write (*,*) 'Sortie de BBCAL2'
  869.  
  870. RETURN
  871.  
  872. END
  873.  
  874.  
  875.  
  876.  

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