Télécharger bbcalc.eso

Retour à la liste

Numérotation des lignes :

  1. C BBCALC SOURCE FANDEUR 10/08/31 21:15:14 6735
  2.  
  3. SUBROUTINE BBCALC(XE,MELE,NBNO,IDIM,NBPGAU,POIGAU,
  4. 1 QSIGAU,ETAGAU,DZEGAU,
  5. 2 NST,LRE,IFOU,A,BB,
  6. 3 NN,SHPTOT,SHP,XDPGE,YDPGE)
  7.  
  8. c=======================================================================
  9. c calcul des composantes de la matrice b-barre-dlatation
  10. c - en 2d: elements ict3,icq4,ict6
  11. c - en 3d: elements ....
  12. c calcul des coefficients de modification de la matrice b-barre-dilatation
  13. c - en 2d: elements icq8
  14. c - en 3d: elements ....
  15. c coefficients {a}={a1,a2,a3}
  16. c=======================================================================
  17. c input
  18. c xe = coordonnees de l'element
  19. c mele = numero de l'element dans nomtp
  20. c nbno = nombre de noeuds
  21. c idim = dimension espace
  22. c nbpgau = nombre de points de gauss pour la rigidite
  23. c poigau = poids des points de gauss
  24. c qsigau,etagau,dzegau = coordonnees des points de gauss
  25. c nst = nombre de composantes de contraintes
  26. c lre = nombre de colonnes de la matrice b
  27. c ifou = ifour dans ccoptio
  28. c nn = numero du mode de fourier
  29. c xdpge,ydpge = coordonnees du point autour duquel se fait
  30. c le mouvement de la section
  31. c
  32. c output
  33. c tableau des coefficients de modification: a(4,*)
  34. c tableau des composantes bbarre-dilatation: bb(3,*)
  35. c=======================================================================
  36.  
  37. IMPLICIT INTEGER(I-N)
  38. IMPLICIT REAL*8 (A-H,O-Z)
  39.  
  40. PARAMETER (XZer=0., XUn=1., X1s2=0.5, X1s4 = 0.25,
  41. & X2s3 = 0.666666666666666666666666666666666666666667,
  42. & X1s3 = 0.333333333333333333333333333333333333333333,
  43. & X1s6 = 0.166666666666666666666666666666666666666667)
  44.  
  45. DIMENSION XE(3,*), POIGAU(*),QSIGAU(*),ETAGAU(*),DZEGAU(*),
  46. & SHP(6,*),SHPTOT(6,NBNO,*),
  47. & A(4,*),BB(3,*)
  48.  
  49. DIMENSION P(4,4),
  50. & XNUM1(8),XNUM2(8),XNUM3(8)
  51.  
  52. *W write(*,*) 'Entree dans le sousprogramme BBCALC'
  53. *W write(*,*) 'Element :',MELE
  54. *W write(*,*) 'Nombre de points Gauss : ', nbpgau
  55.  
  56. IF ( MELE.LT.69 .OR. MELE.GT.78 ) GOTO 999
  57. C=======================================================================
  58. C= Elements incompressibles implementes :
  59. C=======================================================================
  60. C= NOM : ICT3, ICQ4, ICT6, ICQ8, ICC8, ICT4, ICP6, IC20, IC10, IC15
  61. C= MELE : 69 , 70 , 71 , 72 , 73 , 74 , 75 , 76 , 77 , 78
  62. C= Actuellement MELE = 74,75,77,78 elements non traites
  63. C=======================================================================
  64. CALL ZERO(BB,3,60)
  65. GOTO ( 10, 20, 30, 40, 50,666,666, 60,666,666 ) , (MELE-68)
  66. GOTO 999
  67.  
  68. 3 FORMAT (4(A,I1,A,I1,A,F10.4,2X))
  69. 4 FORMAT (3(A,I1,A,I1,A,F10.4,2X))
  70.  
  71. C=======================================================================
  72. C========== Elements MASSIFS INCOMPRESSIBLES BIDIMENSIONNELS ===========
  73. C=======================================================================
  74.  
  75. *-----------------------------------------------------------------------
  76. *----------------------------- Element ICT3 ----------------------------
  77. *-----------------------------------------------------------------------
  78. 10 CONTINUE
  79. * write(*,*) 'Element ICT3',IFOU
  80. * write(*,*) 'Ecriture de la matrice xe[ ]'
  81. * write(*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,3)
  82. IFR = IFOU+4
  83. GOTO (666,101,101,103,666,999),IFR
  84. GOTO 999
  85.  
  86. *--------contraintes planes ou déformations planes (ifour=-2,-1)
  87. 101 CONTINUE
  88. *--------donnee des composantes de bb-dil
  89. r_z = - XE(1,1)*XE(2,3) - XE(1,2)*XE(2,1) + XE(1,2)*XE(2,3)
  90. & + XE(1,1)*XE(2,2) + XE(1,3)*XE(2,1) - XE(1,3)*XE(2,2)
  91. * AIRE = X1s2 * r_z
  92. r_z = XUn / r_z
  93. BB(1,1) = (XE(2,2)-XE(2,3)) * r_z
  94. BB(2,2) = (XE(1,3)-XE(1,2)) * r_z
  95. BB(1,3) = (XE(2,3)-XE(2,1)) * r_z
  96. BB(2,4) = (XE(1,1)-XE(1,3)) * r_z
  97. BB(1,5) = (XE(2,1)-XE(2,2)) * r_z
  98. BB(2,6) = (XE(1,2)-XE(1,1)) * r_z
  99. * write (*,*) 'calcul de aire pour bbdil'
  100. * write (*,*) 'aire = ',AIRE
  101. * write (*,*) 'ecriture de la matrice bb'
  102. * write (*,3) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbno)
  103. GOTO 666
  104. *--------axisymétrique (ifour=0)
  105. 103 CONTINUE
  106. *--------donnee des composantes de bb-dil
  107. r_z = - XE(1,1)*XE(2,3) - XE(1,2)*XE(2,1) + XE(1,2)*XE(2,3)
  108. & + XE(1,1)*XE(2,2) + XE(1,3)*XE(2,1) - XE(1,3)*XE(2,2)
  109. * AIRE = X1s6 * r_z * (XE(1,1)+XE(1,2)+XE(1,3))
  110. r_z = XUn / r_z
  111. BB(1,1) = r_z * (XE(2,2)-XE(2,3))
  112. BB(2,2) = r_z * (XE(1,3)-XE(1,2))
  113. BB(1,3) = r_z * (XE(2,3)-XE(2,1))
  114. BB(2,4) = r_z * (XE(1,1)-XE(1,3))
  115. BB(1,5) = r_z * (XE(2,1)-XE(2,2))
  116. BB(2,6) = r_z * (XE(1,2)-XE(1,1))
  117. *-- Integration analytique des termes BB(3,1),BB(3,3),BB(3,5)
  118. BB(3,1) = XUn / (XE(1,1)+XE(1,2)+XE(1,3))
  119. BB(3,3) = BB(3,1)
  120. BB(3,5) = BB(3,1)
  121. * write(*,*) 'calcul de aire pour bbdil'
  122. * write(*,*) 'aire = ',AIRE
  123. * write(*,*) 'ecriture de la matrice bb'
  124. * write(*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbno)
  125. GOTO 666
  126.  
  127. *-----------------------------------------------------------------------
  128. *----------------------------- Element ICQ4 ----------------------------
  129. *-----------------------------------------------------------------------
  130. 20 CONTINUE
  131. * write(*,*) 'Element ICQ4',IFOU
  132. * write(*,*) 'Ecriture de la matrice xe[ ]'
  133. * write(*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,4)
  134. IFR = IFOU+4
  135. GOTO (666,201,201,203,666,999),IFR
  136. GOTO 999
  137. *--------contraintes planes ou déformations planes (ifour=-2,-1)
  138. 201 CONTINUE
  139. *--------donnée des composantes de bb-dil
  140. r_z = - XE(1,4)*XE(2,3) - XE(1,2)*XE(2,1) + XE(1,1)*XE(2,2)
  141. & - XE(1,1)*XE(2,4) - XE(1,3)*XE(2,2) + XE(1,2)*XE(2,3)
  142. & + XE(1,3)*XE(2,4) + XE(1,4)*XE(2,1)
  143. * AIRE = X1s2 * r_z
  144. r_z = XUn / r_z
  145. BB(1,1) = r_z * (XE(2,2)-XE(2,4))
  146. BB(2,2) = r_z * (XE(1,4)-XE(1,2))
  147. BB(1,3) = r_z * (XE(2,3)-XE(2,1))
  148. BB(2,4) = r_z * (XE(1,1)-XE(1,3))
  149. BB(1,5) = r_z * (XE(2,4)-XE(2,2))
  150. BB(2,6) = r_z * (XE(1,2)-XE(1,4))
  151. BB(1,7) = r_z * (XE(2,1)-XE(2,3))
  152. BB(2,8) = r_z * (XE(1,3)-XE(1,1))
  153. * write (*,*) 'calcul de aire pour bbdil'
  154. * write (*,*) 'aire = ',aire
  155. * write (*,*) 'ecriture de la matrice bb'
  156. * write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbno)
  157. GOTO 666
  158. *
  159. *--------axisymétriques(ifour=0)
  160. 203 CONTINUE
  161. r_z = XE(1,1)*( XE(1,1)*XE(2,2)-XE(1,1)*XE(2,4)-XE(1,4)*XE(2,4)
  162. & +XE(1,4)*XE(2,1)-XE(1,2)*XE(2,1)+XE(1,2)*XE(2,2))
  163. & + XE(1,2)*(-XE(1,2)*XE(2,1)+XE(1,2)*XE(2,3)+XE(1,3)*XE(2,3)
  164. & -XE(1,3)*XE(2,2))
  165. & + XE(1,3)*( XE(1,3)*XE(2,4)+XE(1,4)*XE(2,4)-XE(1,4)*XE(2,3)
  166. & -XE(1,3)*XE(2,2))
  167. & + XE(1,4)*XE(1,4)*(XE(2,1)-XE(2,3))
  168. * AIRE = x1s6 * r_z
  169. r_z = xUn / r_z
  170. BB(1,1) = ( XE(2,2)*(XE(1,2)+XE(1,1))-XE(2,4)*(XE(1,1)+XE(1,4))
  171. & +X1s2*( XE(1,3)*(XE(2,2)-XE(2,4))
  172. & +XE(1,4)*(XE(2,2)+XE(2,3))
  173. & -XE(1,2)*(XE(2,3)+XE(2,4))) ) * r_z
  174. BB(2,2) = ( XE(1,4)*(XE(1,1)+XE(1,4))
  175. & -XE(1,2)*(XE(1,1)+XE(1,2)) ) * r_z
  176. BB(1,3) = ( XE(2,3)*(XE(1,2)+XE(1,3))-XE(2,1)*(XE(1,2)+XE(1,1))
  177. & +X1s2*( XE(1,4)*(XE(2,3)-XE(2,1))
  178. & +XE(1,1)*(XE(2,3)+XE(2,4))
  179. & -XE(1,3)*(XE(2,4)+XE(2,1))) ) * r_z
  180. BB(2,4) = ( XE(1,1)*(XE(1,2)+XE(1,1))
  181. & -XE(1,3)*(XE(1,3)+XE(1,2)) ) * r_z
  182. BB(1,5) = ( XE(2,4)*(XE(1,3)+XE(1,4))-XE(2,2)*(XE(1,3)+XE(1,2))
  183. & +X1s2*( XE(1,1)*(XE(2,4)-XE(2,2))
  184. & +XE(1,2)*(XE(2,4)+XE(2,1))
  185. & -XE(1,4)*(XE(2,1)+XE(2,2))) ) * r_z
  186. BB(2,6) = ( XE(1,2)*(XE(1,2)+XE(1,3))
  187. & -XE(1,4)*(XE(1,4)+XE(1,3)) ) * r_z
  188. BB(1,7) = ( XE(2,1)*(XE(1,4)+XE(1,1))-XE(2,3)*(XE(1,4)+XE(1,3))
  189. & +X1s2*( XE(1,2)*(XE(2,1)-XE(2,3))
  190. & +XE(1,3)*(XE(2,1)+XE(2,2))
  191. & -XE(1,1)*(XE(2,2)+XE(2,3))) ) * r_z
  192. BB(2,8) = ( XE(1,3)*(XE(1,3)+XE(1,4))
  193. & -XE(1,1)*(XE(1,1)+XE(1,4)) ) * r_z
  194. *- Integration analytique des termes BB(3,1),BB(3,3),BB(3,5),BB(3,7)
  195. BB(3,1) = ( XE(1,1)*(XE(2,2)-XE(2,4))+XE(2,1)*(XE(1,4)-XE(1,2))
  196. & +X1s2*( XE(1,2)*(XE(2,3)+XE(2,4))
  197. & +XE(1,3)*(XE(2,4)-XE(2,2))
  198. & -XE(1,4)*(XE(2,3)+XE(2,2))) ) * r_z
  199. BB(3,3) = ( XE(1,2)*(XE(2,3)-XE(2,1))+XE(2,2)*(XE(1,1)-XE(1,3))
  200. & +X1s2*( XE(1,3)*(XE(2,1)+XE(2,4))
  201. & +XE(1,4)*(XE(2,1)-XE(2,3))
  202. & -XE(1,1)*(XE(2,3)+XE(2,4))) ) * r_z
  203. BB(3,5) = ( XE(1,3)*(XE(2,4)-XE(2,2))+XE(2,3)*(XE(1,2)-XE(1,4))
  204. & +X1s2*( XE(1,4)*(XE(2,1)+XE(2,2))
  205. & +XE(1,1)*(XE(2,2)-XE(2,4))
  206. & -XE(1,2)*(XE(2,4)+XE(2,1))) ) * r_z
  207. BB(3,7) = ( XE(1,4)*(XE(2,1)-XE(2,3))+XE(2,4)*(XE(1,3)-XE(1,1))
  208. & +X1s2*( XE(1,1)*(XE(2,2)+XE(2,3))
  209. & +XE(1,2)*(XE(2,3)-XE(2,1))
  210. & -XE(1,3)*(XE(2,1)+XE(2,2))) ) * r_z
  211. * write (*,*) 'calcul de aire pour bbdil'
  212. * write (*,*) 'aire = ',aire
  213. * write (*,*) 'ecriture de la matrice bb'
  214. * write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbno)
  215. GOTO 666
  216.  
  217. *-----------------------------------------------------------------------
  218. *----------------------------- Element ICT6 ----------------------------
  219. *-----------------------------------------------------------------------
  220. 30 CONTINUE
  221. * write(*,*) 'Element ICT6',IFOU
  222. * write (*,*) 'ecriture de la matrice xe[ ]'
  223. * write (*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,6)
  224. IFR = IFOU+4
  225. GOTO (666,301,301,303,666,999),IFR
  226. GOTO 999
  227. *
  228. *--------contraintes planes ou déformations planes(ifour=-2,-1)
  229. 301 CONTINUE
  230. *--------donnée des composantes de bb-dil
  231. r_z = ( XE(1,1)*(XE(2,2)-XE(2,6)) + XE(1,2)*(XE(2,3)-XE(2,1))
  232. & + XE(1,3)*(XE(2,4)-XE(2,2)) + XE(1,4)*(XE(2,5)-XE(2,3))
  233. & + XE(1,5)*(XE(2,6)-XE(2,4)) + XE(1,6)*(XE(2,1)-XE(2,5)) )
  234. & + X1s4 * ( XE(1,1)*(XE(2,5)-XE(2,3))
  235. & + XE(1,3)*(XE(2,1)-XE(2,5))
  236. & + XE(1,5)*(XE(2,3)-XE(2,1)) )
  237. AIRE = X2s3 * r_z
  238. r_z = xUn / r_z
  239. BB(1, 1) = ( (XE(2,2)-XE(2,6)) + X1s4 * (XE(2,5)-XE(2,3)) ) * r_z
  240. BB(2, 2) = ( (XE(1,6)-XE(1,2)) + X1s4 * (XE(1,3)-XE(1,5)) ) * r_z
  241. BB(1, 3) = (XE(2,3)-XE(2,1)) * r_z
  242. BB(2, 4) = (XE(1,1)-XE(1,3)) * r_z
  243. BB(1, 5) = ( (XE(2,4)-XE(2,2)) + X1s4 * (XE(2,1)-XE(2,5)) ) * r_z
  244. BB(2, 6) = ( (XE(1,2)-XE(1,4)) + X1s4 * (XE(1,5)-XE(1,1)) ) * r_z
  245. BB(1, 7) = (XE(2,5)-XE(2,3)) * r_z
  246. BB(2, 8) = (XE(1,3)-XE(1,5)) * r_z
  247. BB(1, 9) = ( (XE(2,6)-XE(2,4)) + X1s4 * (XE(2,3)-XE(2,1)) ) * r_z
  248. BB(2,10) = ( (XE(1,4)-XE(1,6)) + X1s4 * (XE(1,1)-XE(1,3)) ) * r_z
  249. BB(1,11) = (XE(2,1)-XE(2,5)) * r_z
  250. BB(2,12) = (XE(1,5)-XE(1,1)) * r_z
  251.  
  252. * write (*,*) 'calcul de aire pour bbdil'
  253. * write (*,*) 'aire = ',aire
  254. * write (*,*) 'ecriture de la matrice bb'
  255. * write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbno)
  256. GOTO 666
  257. *
  258. *--------axisymétriques (ifour=0)
  259. 303 CONTINUE
  260. *--------donnée des composantes de bb-dil
  261. *-----intégration numérique des composantes bb-dil------
  262. AIRE = XZer
  263. DO i = 1,NBNO
  264. XNUM1(i) = XZer
  265. XNUM2(i) = XZer
  266. XNUM3(i) = XZer
  267. ENDDO
  268. DO 3031 IGAU = 1,NBPGAU
  269. DO i = 1,NBNO
  270. SHP(1,i) = SHPTOT(1,i,IGAU)
  271. SHP(2,i) = SHPTOT(2,i,IGAU)
  272. SHP(3,i) = SHPTOT(3,i,IGAU)
  273. ENDDO
  274. CALL DISTRR(XE,SHP,NBNO,RR)
  275. CALL DISTR2(XE,SHP,NBNO,RR2)
  276. CALL DISTR3(XE,SHP,NBNO,RR3)
  277. CALL DISTZ2(XE,SHP,NBNO,ZZ2)
  278. CALL DISTZ3(XE,SHP,NBNO,ZZ3)
  279. CALL JACOBI(XE,SHP,2,NBNO,DJAC)
  280. r_z = POIGAU(IGAU) * RR
  281. r_zz2 = r_z * ZZ2
  282. r_zz3 = r_z * ZZ3
  283. r_rr2 = r_z * RR2
  284. r_rr3 = r_z * RR3
  285. r_z = POIGAU(IGAU) * DJAC
  286. AIRE = AIRE + (r_z * RR)
  287. DO i = 1,NBNO
  288. XNUM1(i) = XNUM1(i) + ( SHPTOT(2,i,IGAU)*r_zz3
  289. & - SHPTOT(3,i,IGAU)*r_zz2 )
  290. XNUM2(i) = XNUM2(i) + ( SHPTOT(3,i,IGAU)*r_rr2
  291. & - SHPTOT(2,i,IGAU)*r_rr3 )
  292. XNUM3(i) = XNUM3(i) + ( SHPTOT(1,i,IGAU)*r_z )
  293. ENDDO
  294. 3031 CONTINUE
  295. r_z = XUn / AIRE
  296. DO i = 1,NBNO
  297. k = 2*i
  298. BB(1,k-1) = XNUM1(i) * r_z
  299. BB(2,k) = XNUM2(i) * r_z
  300. BB(3,k-1) = XNUM3(i) * r_z
  301. ENDDO
  302. * write (*,*) 'calcul de aire pour bbdil'
  303. * write (*,*) 'aire = ',aire
  304. * write (*,*) 'ecriture de la matrice bb'
  305. * write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbno)
  306. GOTO 666
  307.  
  308. *-----------------------------------------------------------------------
  309. *----------------------------- Element ICQ8 ----------------------------
  310. *-----------------------------------------------------------------------
  311. 40 CONTINUE
  312. * write(*,*) 'Element ICQ8',IFOU
  313. * write (*,*) 'Ecriture de la matrice xe[ ]'
  314. * write (*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,6)
  315.  
  316. IFR = IFOU+4
  317. GOTO (666,401,401,403,666,999),IFR
  318. GOTO 666
  319.  
  320. *--------contraintes planes ou déformations planes (ifour=-2,-1)
  321. 401 CONTINUE
  322. *--------donnée de la matrice symétrique p(3,3) (1 fois par élément)
  323. P(1,1) = X2s3 * ( XE(1,1)*(XE(2,2)-XE(2,8))
  324. 1 + XE(1,2)*(XE(2,3)-XE(2,1))
  325. 1 +XE(1,3)*(XE(2,4)-XE(2,2))+XE(1,4)*(XE(2,5)-XE(2,3))
  326. 2 +XE(1,5)*(XE(2,6)-XE(2,4))+XE(1,6)*(XE(2,7)-XE(2,5))
  327. 3 +XE(1,7)*(XE(2,8)-XE(2,6))+XE(1,8)*(XE(2,1)-XE(2,7)) )
  328. 4 + X1s6 * ( (XE(1,1)-XE(1,5))*(XE(2,7)-XE(2,3))
  329. 5 +(XE(1,3)-XE(1,7))*(XE(2,1)-XE(2,5)) )
  330. P(1,2) = 5.D0/9.D0*(-XE(2,2)*(XE(1,1)+XE(1,3))
  331. 1 +XE(1,2)*(XE(2,1)+XE(2,3))+XE(1,6)*(-XE(2,5)-XE(2,7))
  332. 1 +XE(2,6)*(XE(1,5)+XE(1,7)))
  333. 2 +4.D0/9.D0*(XE(1,8)*(XE(2,7)+XE(2,2)-XE(2,6)-XE(2,1))
  334. 2 +XE(2,8)*(-XE(1,7)-XE(1,2)+XE(1,6)+XE(1,1))
  335. 2 +XE(1,4)*(XE(2,2)-XE(2,6)+XE(2,5)-XE(2,3))
  336. 2 +XE(2,4)*(-XE(1,2)+XE(1,6)-XE(1,5)+XE(1,3)))
  337. 3 +7.D0/45.D0*(XE(1,2)*(XE(2,5)+XE(2,7))
  338. 3 -XE(2,2)*(XE(1,5)+XE(1,7))-XE(1,6)*(XE(2,1)+XE(2,3))
  339. 3 +XE(2,6)*(XE(1,1)+XE(1,3)))
  340. 4 +8.D0/15.D0*(XE(1,6)*XE(2,2)-XE(1,2)*XE(2,6))
  341. 5 +7.D0/90.D0*(XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5)
  342. 5 +XE(1,5)*XE(2,3)-XE(1,1)*XE(2,7))
  343. 6 +1.D0/30.D0*(-XE(1,7)*XE(2,3)+XE(1,3)*XE(2,7)
  344. 6 -XE(1,5)*XE(2,1)+XE(1,1)*XE(2,5))
  345. P(1,3) = 5.D0/9.D0*(-XE(1,8)*(XE(2,7)+XE(2,1))
  346. 1 +XE(2,8)*(XE(1,7)+XE(1,1))
  347. 1 +XE(1,4)*(XE(2,3)+XE(2,5))
  348. 1 -XE(2,4)*(XE(1,3)+XE(1,5)))
  349. 2 +4.D0/9.D0*((XE(1,8)-XE(1,4))*(XE(2,2)+XE(2,6))
  350. 2 +(-XE(2,8)+XE(2,4))*(XE(1,2)+XE(1,6))
  351. 2 +XE(1,2)*(XE(2,1)-XE(2,3))+XE(2,2)*(-XE(1,1)+XE(1,3))
  352. 2 +XE(1,6)*(XE(2,7)-XE(2,5))+XE(2,6)*(-XE(1,7)+XE(1,5)))
  353. 3 +7.D0/45.D0*(-XE(1,8)*(XE(2,3)+XE(2,5))
  354. 3 +XE(2,8)*(XE(1,3)+XE(1,5))
  355. 3 +XE(1,4)*(XE(2,1)+XE(2,7))
  356. 3 -XE(2,4)*(XE(1,1)+XE(1,7)))
  357. 4 +1.D0/30.D0*(-XE(1,7)*XE(2,3)+XE(1,3)*XE(2,7)
  358. 4 +XE(1,5)*XE(2,1)-XE(1,1)*XE(2,5))
  359. 5 +7.D0/90.D0*(XE(1,7)*XE(2,5)-XE(1,3)*XE(2,1)
  360. 5 -XE(1,5)*XE(2,7)+XE(1,1)*XE(2,3))
  361. 6 +8.D0/15.D0*(XE(1,8)*XE(2,4)-XE(1,4)*XE(2,8))
  362. P(1,4) = XZer
  363.  
  364. P(2,1) = P(1,2)
  365. P(2,2) = 16.D0/45.D0*(XE(2,6)*(XE(1,5)-XE(1,7))
  366. 1 +XE(1,6)*(XE(2,7)-XE(2,5))+XE(1,2)*(XE(2,3)-XE(2,1))
  367. 1 +XE(2,2)*(XE(1,1)-XE(1,3)))
  368. 2 +14.D0/45.D0*(XE(2,4)*(XE(1,3)-XE(1,5))
  369. 2 +XE(1,4)*(XE(2,5)-XE(2,3))+XE(1,8)*(XE(2,1)-XE(2,7))
  370. 2 +XE(2,8)*(XE(1,7)-XE(1,1)))
  371. 3 +8.D0/45.D0*(XE(1,4)*(XE(2,2)-XE(2,6))
  372. 3 +XE(1,8)*(XE(2,6)-XE(2,2))+XE(1,2)*(XE(2,8)-XE(2,4))
  373. 3 +XE(1,6)*(XE(2,4)-XE(2,8)))
  374. 4 +2.D0/45.D0*(XE(2,2)*(XE(1,5)-XE(1,7))
  375. 4 +XE(1,2)*(XE(2,7)-XE(2,5)) +XE(1,6)*(XE(2,3)-XE(2,1))
  376. 4 +XE(2,6)*(XE(1,1)-XE(1,3)))
  377. 5 +4.D0/45.D0*(XE(2,4)*(XE(1,1)-XE(1,7))
  378. 5 +XE(1,4)*(XE(2,7)-XE(2,1)) +XE(1,8)*(XE(2,3)-XE(2,5))
  379. 5 +XE(2,8)*(XE(1,5)-XE(1,3)))
  380. 6 +17.D0/90.D0*(XE(1,7)*XE(2,5)+XE(1,3)*XE(2,1)
  381. 6 -XE(1,5)*XE(2,7)-XE(1,1)*XE(2,3))
  382. 7 +1.D0/90.D0*(-XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5)
  383. 7 +XE(1,5)*XE(2,3)+XE(1,1)*XE(2,7))
  384. P(2,3) = 1.D0/3.D0*(XE(1,5)*(XE(2,6)-XE(2,4))
  385. 1 +XE(1,8)*(XE(2,7)+XE(2,1))+XE(1,7)*(XE(2,6)-XE(2,8))
  386. 1 -XE(1,2)*(XE(2,1)+XE(2,3))+XE(1,1)*(XE(2,2)-XE(2,8))
  387. 1 +XE(1,4)*(XE(2,5)+XE(2,3))-XE(1,6)*(XE(2,7)+XE(2,5))
  388. 1 +XE(1,3)*(XE(2,2)-XE(2,4)))
  389. 2 +4.D0/9.D0*(-XE(1,4)*(XE(2,2)+XE(2,6))
  390. 2 -XE(1,8)*(XE(2,6)+XE(2,2))+XE(1,2)*(XE(2,8)+XE(2,4))
  391. 2 +XE(1,6)*(XE(2,4)+XE(2,8)))
  392. 3 +1.D0/9.D0*(XE(1,7)*(XE(2,2)-XE(2,4))
  393. 3 +XE(1,8)*(XE(2,3)+XE(2,5))
  394. 3 +XE(1,5)*(XE(2,2)-XE(2,8))+XE(1,3)*(XE(2,6)-XE(2,8))
  395. 3 -XE(1,2)*(XE(2,5)+XE(2,7))+XE(1,1)*(XE(2,6)-XE(2,4))
  396. 3 +XE(1,4)*(XE(2,1)+XE(2,7))-XE(1,6)*(XE(2,1)+XE(2,3)))
  397. P(2,4) = XZer
  398.  
  399. P(3,1) = P(1,3)
  400. P(3,2) = P(2,3)
  401. P(3,3) = 14.D0/45.D0*(XE(1,6)*(XE(2,7)-XE(2,5))
  402. 1 +XE(2,6)*(XE(1,5)-XE(1,7))+XE(1,2)*(XE(2,3)-XE(2,1))
  403. 1 +XE(2,2)*(XE(1,1)-XE(1,3)))
  404. 2 +16.D0/45.D0*(XE(1,8)*(XE(2,1)-XE(2,7))
  405. 2 +XE(2,8)*(XE(1,7)-XE(1,1))+XE(1,4)*(XE(2,5)-XE(2,3))
  406. 2 +XE(2,4)*(XE(1,3)-XE(1,5)))
  407. 3 +8.D0/45.D0*(XE(1,4)*(XE(2,2)-XE(2,6))
  408. 3 +XE(1,8)*(XE(2,6)-XE(2,2))+XE(1,2)*(XE(2,8)-XE(2,4))
  409. 3 +XE(1,6)*(XE(2,4)-XE(2,8)))
  410. 4 +2.D0/45.D0*(XE(2,4)*(XE(1,7)-XE(1,1))
  411. 4 +XE(1,8)*(XE(2,5)-XE(2,3))+XE(2,8)*(-XE(1,5)+XE(1,3))
  412. 4 +XE(1,4)*(XE(2,1)-XE(2,7)))
  413. 5 +4.D0/45.D0*(XE(2,2)*(XE(1,7)-XE(1,5))
  414. 5 +XE(1,2)*(XE(2,5)-XE(2,7))+XE(1,6)*(XE(2,1)-XE(2,3))
  415. 5 +XE(2,6)*(XE(1,3)-XE(1,1)))
  416. 6 +1.D0/90.D0*(-XE(1,5)*XE(2,7)+XE(1,3)*XE(2,1)
  417. 6 +XE(1,7)*XE(2,5)-XE(1,1)*XE(2,3))
  418. 7 +17.D0/90.D0*(-XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5)
  419. 7 +XE(1,5)*XE(2,3)+XE(1,1)*XE(2,7))
  420. P(3,4) = XZer
  421.  
  422. P(4,1) = XZer
  423. P(4,2) = XZer
  424. P(4,3) = XZer
  425. P(4,4) = XZer
  426.  
  427. * write (*,*) 'ecriture de la matrice p[ ]'
  428. * write (*,4) (('p(',i,',',j,')=',p(i,j),i=1,3),j=1,3)
  429.  
  430. *-----------------donnée des vecteurs {t}={t1,t2,t3}--------------------
  431. A(1,1) = 2.D0/3.D0*(XE(2,2)-XE(2,8))
  432. 1 +1.D0/6.D0*(XE(2,7)-XE(2,3))
  433. A(2,1) = 1.D0/9.D0*(4.D0*XE(2,8)-5.D0*XE(2,2))
  434. 1 +7.D0/90.D0*(-XE(2,7)+2.D0*XE(2,6))
  435. 1 +1.D0/30.D0*XE(2,5)
  436. A(3,1) = 1.D0/9.D0*(-4.D0*XE(2,2)+5.D0*XE(2,8))
  437. 1 +7.D0/90.D0*(XE(2,3)-2.D0*XE(2,4))
  438. 1 -1.D0/30.D0*XE(2,5)
  439.  
  440. A(1,2) = 2.D0/3.D0*(XE(1,8)-XE(1,2))
  441. 1 +1.D0/6.D0*(XE(1,3)-XE(1,7))
  442. A(2,2) = 1.D0/9.D0*(-4.D0*XE(1,8)+5.D0*XE(1,2))
  443. 1 +7.D0/90.D0*(XE(1,7)-2.D0*XE(1,6))
  444. 1 -1.D0/30.D0*XE(1,5)
  445. A(3,2) = 1.D0/9.D0*(4.D0*XE(1,2)-5.D0*XE(1,8))
  446. 1 +7.D0/90.D0*(-XE(1,3)+2.D0*XE(1,4))
  447. 1 +1.D0/30.D0*XE(1,5)
  448.  
  449. A(1,3) = 2.D0/3.D0*(-XE(2,1)+XE(2,3))
  450. A(2,3) = 5.D0/9.D0*(XE(2,1)+XE(2,3))
  451. 1 +7.D0/45.D0*(XE(2,5)+XE(2,7))
  452. 1 -4.D0/9.D0*(XE(2,4)+XE(2,8))
  453. 1 -8.D0/15.D0*XE(2,6)
  454. A(3,3) = 4.D0/9.D0*(XE(2,1)-XE(2,3)+XE(2,4)-XE(2,8))
  455.  
  456. A(1,4) = 2.D0/3.D0*(-XE(1,3)+XE(1,1))
  457. A(2,4) = -5.D0/9.D0*(XE(1,1)+XE(1,3))
  458. 1 -7.D0/45.D0*(XE(1,5)+XE(1,7))
  459. 1 +4.D0/9.D0*(XE(1,4)+XE(1,8))
  460. 1 +8.D0/15.D0*XE(1,6)
  461. A(3,4) = 4.D0/9.D0*(-XE(1,1)+XE(1,3)-XE(1,4)+XE(1,8))
  462.  
  463. A(1,5) = 2.D0/3.D0*(XE(2,4)-XE(2,2))
  464. 1 +1.D0/6.D0*(-XE(2,5)+XE(2,1))
  465. A(2,5) = 1.D0/9.D0*(4.D0*XE(2,4)-5.D0*XE(2,2))
  466. 1 +7.D0/90.D0*(-XE(2,5)+2.D0*XE(2,6))
  467. 1 +1.D0/30.D0*XE(2,7)
  468. A(3,5) = 1.D0/9.D0*(4.D0*XE(2,2)-5.D0*XE(2,4))
  469. 1 +7.D0/90.D0*(-XE(2,1)+2.D0*XE(2,8))
  470. 1 +1.D0/30.D0*XE(2,7)
  471.  
  472. A(1,6) = 2.D0/3.D0*(-XE(1,4)+XE(1,2))
  473. 1 +1.D0/6.D0*(XE(1,5)-XE(1,1))
  474. A(2,6) = 1.D0/9.D0*(-4.D0*XE(1,4)+5.D0*XE(1,2))
  475. 1 +7.D0/90.D0*(XE(1,5)-2.D0*XE(1,6))
  476. 1 -1.D0/30.D0*XE(1,7)
  477. A(3,6) = 1.D0/9.D0*(-4.D0*XE(1,2)+5.D0*XE(1,4))
  478. 1 +7.D0/90.D0*(XE(1,1)-2.D0*XE(1,8))
  479. 1 -1.D0/30.D0*XE(1,7)
  480.  
  481. A(1,7) = 2.D0/3.D0*(-XE(2,3)+XE(2,5))
  482. A(2,7) = 4.D0/9.D0*(XE(2,5)-XE(2,3)+XE(2,2)-XE(2,6))
  483. A(3,7) = 5.D0/9.D0*(XE(2,3)+XE(2,5))
  484. 1 +7.D0/45.D0*(XE(2,1)+XE(2,7))
  485. 1 -4.D0/9.D0*(XE(2,2)+XE(2,6))
  486. 1 -8.D0/15.D0*XE(2,8)
  487.  
  488. A(1,8) = 2.D0/3.D0*(XE(1,3)-XE(1,5))
  489. A(2,8) = 4.D0/9.D0*(-XE(1,5)+XE(1,3)-XE(1,2)+XE(1,6))
  490. A(3,8) = -5.D0/9.D0*(XE(1,3)+XE(1,5))
  491. 1 -7.D0/45.D0*(XE(1,1)+XE(1,7))
  492. 1 +4.D0/9.D0*(XE(1,2)+XE(1,6))
  493. 1 +8.D0/15.D0*XE(1,8)
  494.  
  495. A(1,9) = 2.D0/3.D0*(-XE(2,4)+XE(2,6))
  496. 1 +1.D0/6.D0*(XE(2,3)-XE(2,7))
  497. A(2,9) = 1.D0/9.D0*(-4.D0*XE(2,4)+5.D0*XE(2,6))
  498. 1 +7.D0/90.D0*(XE(2,3)-2.D0*XE(2,2))
  499. 1 -1.D0/30.D0*XE(2,1)
  500. A(3,9) = 1.D0/9.D0*(4.D0*XE(2,6)-5.D0*XE(2,4))
  501. 1 +7.D0/90.D0*(-XE(2,7)+2.D0*XE(2,8))
  502. 1 +1.D0/30.D0*XE(2,1)
  503.  
  504. A(1,10) = 2.D0/3.D0*(XE(1,4)-XE(1,6))
  505. 1 +1.D0/6.D0*(-XE(1,3)+XE(1,7))
  506. A(2,10) = 1.D0/9.D0*(4.D0*XE(1,4)-5.D0*XE(1,6))
  507. 1 +7.D0/90.D0*(-XE(1,3)+2.D0*XE(1,2))
  508. 1 +1.D0/30.D0*XE(1,1)
  509. A(3,10) = 1.D0/9.D0*(-4.D0*XE(1,6)+5.D0*XE(1,4))
  510. 1 +7.D0/90.D0*(XE(1,7)-2.D0*XE(1,8))
  511. 1 -1.D0/30.D0*XE(1,1)
  512.  
  513. A(1,11) = 2.D0/3.D0*(-XE(2,5)+XE(2,7))
  514. A(2,11) = -5.D0/9.D0*(XE(2,7)+XE(2,5))
  515. 1 -7.D0/45.D0*(XE(2,1)+XE(2,3))
  516. 1 +4.D0/9.D0*(XE(2,4)+XE(2,8))
  517. 1 +8.D0/15.D0*XE(2,2)
  518. A(3,11) = 4.D0/9.D0*(XE(2,4)-XE(2,5)-XE(2,8)+XE(2,7))
  519.  
  520. A(1,12) = 2.D0/3.D0*(XE(1,5)-XE(1,7))
  521. A(2,12) = 5.D0/9.D0*(XE(1,7)+XE(1,5))
  522. 1 +7.D0/45.D0*(XE(1,1)+XE(1,3))
  523. 1 -4.D0/9.D0*(XE(1,4)+XE(1,8))
  524. 1 -8.D0/15.D0*XE(1,2)
  525. A(3,12) = 4.D0/9.D0*(-XE(1,4)+XE(1,5)+XE(1,8)-XE(1,7))
  526.  
  527. A(1,13) = 2.D0/3.D0*(-XE(2,6)+XE(2,8))
  528. 1 +1.D0/6.D0*(XE(2,5)-XE(2,1))
  529. A(2,13) = 1.D0/9.D0*(-4.D0*XE(2,8)+5.D0*XE(2,6))
  530. 1 +7.D0/90.D0*(XE(2,1)-2.D0*XE(2,2))
  531. 1 -1.D0/30.D0*XE(2,3)
  532. A(3,13) = 1.D0/9.D0*(-4.D0*XE(2,6)+5.D0*XE(2,8))
  533. 1 +7.D0/90.D0*(XE(2,5)-2.D0*XE(2,4))
  534. 1 -1.D0/30.D0*XE(2,3)
  535.  
  536. A(1,14) = 2.D0/3.D0*(XE(1,6)-XE(1,8))
  537. 1 +1.D0/6.D0*(-XE(1,5)+XE(1,1))
  538. A(2,14) = 1.D0/9.D0*(4.D0*XE(1,8)-5.D0*XE(1,6))
  539. 1 +7.D0/90.D0*(-XE(1,1)+2.D0*XE(1,2))
  540. 1 +1.D0/30.D0*XE(1,3)
  541. A(3,14) = 1.D0/9.D0*(4.D0*XE(1,6)-5.D0*XE(1,8))
  542. 1 +7.D0/90.D0*(-XE(1,5)+2.D0*XE(1,4))
  543. 1 +1.D0/30.D0*XE(1,3)
  544.  
  545. A(1,15) = 2.D0/3.D0*(-XE(2,7)+XE(2,1))
  546. A(2,15) = 4.D0/9.D0*(-XE(2,6)+XE(2,7)+XE(2,2)-XE(2,1))
  547. A(3,15) = -5.D0/9.D0*(XE(2,1)+XE(2,7))
  548. 1 -7.D0/45.D0*(XE(2,3)+XE(2,5))
  549. 1 +4.D0/9.D0*(XE(2,2)+XE(2,6))
  550. 1 +8.D0/15.D0*XE(2,4)
  551.  
  552. A(1,16) = 2.D0/3.D0*(XE(1,7)-XE(1,1))
  553. A(2,16) = 4.D0/9.D0*(XE(1,6)-XE(1,7)-XE(1,2)+XE(1,1))
  554. A(3,16) = 5.D0/9.D0*(XE(1,1)+XE(1,7))
  555. 1 +7.D0/45.D0*(XE(1,3)+XE(1,5))
  556. 1 -4.D0/9.D0*(XE(1,2)+XE(1,6))
  557. 1 -8.D0/15.D0*XE(1,4)
  558.  
  559. * write (*,*) 'ecriture de la matrice t[ ]'
  560. * write (*,4) (('t(',i,',',j,')=',a(i,j),i=1,3),j=1,16)
  561. *--------------résolution matricielle de [p]{a}={t}---------------------
  562. * write (*,*) 'Appel a GAUSSJ'
  563. N = 3
  564. NP = 4
  565. M = 16
  566. MP = 60
  567. CALL GAUSSJ(P,N,NP,A,M,MP)
  568. GOTO 666
  569.  
  570. *--------axisymétriques (ifour=0)
  571. 403 CONTINUE
  572. *--------donnée de la matrice symétrique p(3,3) (1 fois par élément)
  573. CALL ZERO(P,4,4)
  574. DO IGAU = 1,NBPGAU
  575. DO I = 1,NBNO
  576. SHP(1,I) = SHPTOT(1,I,IGAU)
  577. SHP(2,I) = SHPTOT(2,I,IGAU)
  578. SHP(3,I) = SHPTOT(3,I,IGAU)
  579. ENDDO
  580. CALL DISTRR(XE,SHP,NBNO,RR)
  581. CALL JACOBI(XE,SHP,2,NBNO,DJAC)
  582. r_z = POIGAU(IGAU)*RR*DJAC
  583. P(1,1) = P(1,1) + r_z
  584. P(2,1) = P(2,1) + (r_z * QSIGAU(IGAU))
  585. P(3,1) = P(3,1) + (r_z * ETAGAU(IGAU))
  586. P(2,2) = P(2,2) + (r_z * QSIGAU(IGAU) * QSIGAU(IGAU))
  587. P(3,2) = P(3,2) + (r_z * QSIGAU(IGAU) * ETAGAU(IGAU))
  588. P(3,3) = P(3,3) + (r_z * ETAGAU(IGAU) * ETAGAU(IGAU))
  589. ENDDO
  590. P(1,2) = P(2,1)
  591. P(1,3) = P(3,1)
  592. P(2,3) = P(3,2)
  593. * write (*,*) 'ecriture de la matrice p[ ]'
  594. * write (*,4) (('p(',i,',',j,')=',p(i,j),i=1,3),j=1,3)
  595. *-----------calcul des vecteurs {t}={t1,t2,t3}--intégration numérique---
  596. CALL ZERO(A,4,24)
  597. DO 4031 IGAU = 1,NBPGAU
  598. DO I = 1,NBNO
  599. SHP(1,I) = SHPTOT(1,I,IGAU)
  600. SHP(2,I) = SHPTOT(2,I,IGAU)
  601. SHP(3,I) = SHPTOT(3,I,IGAU)
  602. ENDDO
  603. CALL DISTRR(XE,SHP,NBNO,RR)
  604. CALL DISTR2(XE,SHP,NBNO,RR2)
  605. CALL DISTR3(XE,SHP,NBNO,RR3)
  606. CALL DISTZ2(XE,SHP,NBNO,ZZ2)
  607. CALL DISTZ3(XE,SHP,NBNO,ZZ3)
  608. CALL JACOBI(XE,SHP,2,NBNO,DJAC)
  609. K = 0
  610. DO I = 1,NBNO
  611. K = K+1
  612. r_z = POIGAU(IGAU)*RR
  613. & *(SHPTOT(2,I,IGAU)*ZZ3-SHPTOT(3,I,IGAU)*ZZ2)
  614. A(1,K) = A(1,K) + r_z
  615. A(2,K) = A(2,K) + r_z * QSIGAU(IGAU)
  616. A(3,K) = A(3,K) + r_z * ETAGAU(IGAU)
  617. K = K+1
  618. r_z = POIGAU(IGAU)*RR
  619. & *(SHPTOT(2,I,IGAU)*RR3-SHPTOT(3,I,IGAU)*RR2)
  620. A(1,K) = A(1,K) + r_z
  621. A(2,K) = A(2,K) + r_z * QSIGAU(IGAU)
  622. A(3,K) = A(3,K) + r_z * ETAGAU(IGAU)
  623. K = K+2
  624. r_z = POIGAU(IGAU)*SHPTOT(1,I,IGAU)*DJAC
  625. A(1,K) = A(1,K) + r_z
  626. A(2,K) = A(2,K) + r_z * QSIGAU(IGAU)
  627. A(3,K) = A(3,K) + r_z * ETAGAU(IGAU)
  628. ENDDO
  629. 4031 CONTINUE
  630. * write(*,*) 'Ecriture de la matrice t[ ]'
  631. * write(*,4) (('t(',i,',',j,')=',a(i,j),i=1,3),j=1,24)
  632. *------------résolution matricielle de [p]{a}={t}-----------------
  633. N = 3
  634. NP = 4
  635. M = 24
  636. MP = 60
  637. * write (*,*) 'appel a GAUSSJ'
  638. CALL GAUSSJ(P,N,NP,A,M,MP)
  639. * write(*,*) 'Ecriture de la matrice a[ ]'
  640. * write(*,4) (('a(',i,',',j,')=',a(i,j),i=1,3),j=1,24)
  641. GOTO 666
  642.  
  643. C=======================================================================
  644. C========== Elements MASSIFS INCOMPRESSIBLES TRIDIMENSIONNELS ==========
  645. C=======================================================================
  646.  
  647. *-----------------------------------------------------------------------
  648. *----------------------------- Element ICC8 ----------------------------
  649. *-----------------------------------------------------------------------
  650. 50 CONTINUE
  651. * write(*,*) 'Element ICC8',IFOU
  652. * write(*,*) 'Ecriture de la matrice xe[ ]'
  653. * write(*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,3),j=1,8)
  654. IF (IFOU.NE.2) GOTO 999
  655. GOTO 666
  656.  
  657. *-----------------------------------------------------------------------
  658. *----------------------------- Element IC20 ----------------------------
  659. *-----------------------------------------------------------------------
  660. 60 CONTINUE
  661. * write(*,*) 'Element IC20',IFOU
  662. * write(*,*) 'Ecriture de la matrice xe[ ]'
  663. * write(*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,3),j=1,20)
  664. IF (IFOU.NE.2) GOTO 999
  665. GOTO 666
  666.  
  667. C=======================================================================
  668. C======================= ERREUR : Cas non prevus =======================
  669. C=======================================================================
  670. 999 CONTINUE
  671. CALL ERREUR(5)
  672.  
  673. 666 CONTINUE
  674. * write (*,*) 'Sortie de BBCALC'
  675.  
  676. RETURN
  677. END
  678.  
  679.  
  680.  

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