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