38
39
40
42
43
44
45#include "implicit_f.inc"
46
47
48
49#include "mvsiz_p.inc"
50
51
52
53#include "param_c.inc"
54#include "pt8_c.inc"
55
56
57
58 INTEGER :: IGEO(NPROPGI,*),NGL(*),NGEO(*),LVLOC,,JEUL
60 . vol(*),det(*), veul(lvloc,*),geo(npropg,*),
61 . x1(*), x2(*), x3(*), x4(*), x5(*), x6(*), x7(*), x8(*),
62 . y1(*), y2(*), y3(*), y4(*), y5(*), y6(*), y7(*), y8(*),
63 . z1(*), z2(*), z3(*), z4(*), z5(*), z6(*), z7(*), z8(*),
64 . deltax(*)
65
66
67
68 INTEGER I, J, JPT, IPT
70 . aj11(mvsiz), aj12(mvsiz), aj13(mvsiz), aj21(mvsiz),
71 . aj22(mvsiz), aj23(mvsiz), aj31(mvsiz), aj32(mvsiz),
72 . aj33(mvsiz), ai11(mvsiz), ai12(mvsiz), ai13(mvsiz),
73 . ai21(mvsiz), ai22(mvsiz), ai23(mvsiz), ai31(mvsiz),
74 . ai32(mvsiz), ai33(mvsiz), a1111(mvsiz), a1113(mvsiz),
75 . a1115(mvsiz), a1117(mvsiz), a2111(mvsiz), a2113(mvsiz),
76 . a2115(mvsiz), a2117(mvsiz), a3111(mvsiz), a3113(mvsiz),
77 . a3115(mvsiz), a3117(mvsiz), a1221(mvsiz), a1222(mvsiz),
78 . a1225(mvsiz), a1226(mvsiz), a2221(mvsiz), a2222(mvsiz),
79 . a2225(mvsiz), a2226(mvsiz), a3221(mvsiz), a3222(mvsiz),
80 . a3225(mvsiz), a3226(mvsiz), a1331(mvsiz), a1332(mvsiz),
81 . a1333(mvsiz), a1334(mvsiz), a2331(mvsiz), a2332(mvsiz),
82 . a2333(mvsiz), a2334(mvsiz), a3331(mvsiz), a3332(mvsiz),
83 . a3333(mvsiz), a3334(mvsiz), x12(mvsiz), x34(mvsiz), x56(mvsiz),
84 . x78(mvsiz), y12(mvsiz), y34(mvsiz), y56(mvsiz), y78(mvsiz),
85 . z12(mvsiz), z34(mvsiz), z56(mvsiz), z78(mvsiz), x14(mvsiz),
86 . x23(mvsiz), x58(mvsiz), x67(mvsiz), y14(mvsiz), y23(mvsiz),
87 . y58(mvsiz), y67(mvsiz), z14(mvsiz), z23(mvsiz), z58(mvsiz),
88 . z67(mvsiz), x15(mvsiz), x26(mvsiz), x37(mvsiz), x48(mvsiz),
89 . y15(mvsiz), y26(mvsiz), y37(mvsiz), y48(mvsiz), z15(mvsiz),
90 . z26(mvsiz), z37(mvsiz), z48(mvsiz),
91 . vlinv(mvsiz),vlinc(mvsiz,8), h(8),
92 . px1(mvsiz), px2(mvsiz), px3(mvsiz), px4(mvsiz),
93 . px5(mvsiz), px6(mvsiz), px7(mvsiz), px8(mvsiz),
94 . py1(mvsiz), py2(mvsiz), py3(mvsiz), py4(mvsiz),
95 . py5(mvsiz), py6(mvsiz), py7(mvsiz), py8(mvsiz),
96 . pz1(mvsiz), pz2(mvsiz), pz3(mvsiz), pz4(mvsiz),
97 . pz5(mvsiz), pz6(mvsiz), pz7(mvsiz), pz8(mvsiz)
99 .
poro, smax, vmin, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3
100
101 DO i=1,nel
102 x12(i)=x1(i)-x2(i)
103 x34(i)=x3(i)-x4(i)
104 x56(i)=x5(i)-x6(i)
105 x78(i)=x7(i)-x8(i)
106 x14(i)=x1(i)-x4(i)
107 x23(i)=x2(i)-x3(i)
108 x58(i)=x5(i)-x8(i)
109 x67(i)=x6(i)-x7(i)
110 x15(i)=x1(i)-x5(i)
111 x26(i)=x2(i)-x6(i)
112 x37(i)=x3(i)-x7(i)
113 x48(i)=x4(i)-x8(i)
114
115 y12(i)=y1(i)-y2(i)
116 y34(i)=y3(i)-y4(i)
117 y56(i)=y5(i)-y6(i)
118 y78(i)=y7(i)-y8(i)
119 y14(i)=y1(i)-y4(i)
120 y23(i)=y2(i)-y3(i)
121 y58(i)=y5(i)-y8(i)
122 y67(i)=y6(i)-y7(i)
123 y15(i)=y1(i)-y5(i)
124 y26(i)=y2(i)-y6(i)
125 y37(i)=y3(i)-y7(i)
126 y48(i)=y4(i)-y8(i)
127
128 z12(i)=z1(i)-z2(i)
129 z34(i)=z3(i)-z4(i)
130 z56(i)=z5(i)-z6(i)
131 z78(i)=z7(i)-z8(i)
132 z14(i)=z1(i)-z4(i)
133 z23(i)=z2(i)-z3(i)
134 z58(i)=z5(i)-z8(i)
135 z67(i)=z6(i)-z7(i)
136 z15(i)=z1(i)-z5(i)
137 z26(i)=z2(i)-z6(i)
138 z37(i)=z3(i)-z7(i)
139 z48(i)=z4(i)-z8(i)
140
141 det(i)=zero
142 DO j=1,12
143 veul(j,i)=zero
144 ENDDO
145 DO j=32,51
146 veul(j,i)=zero
147 ENDDO
148
149 ENDDO
150
151 DO jpt=1,8
152 ipt=jpt
153
154
155 CALL basisf (h,p11,p21,p31,ipt)
156
157 DO i=1,nel
158 aj11(i)=p11*x12(i)+p13*x34(i)+p15*x56(i)+p17*x78(i)
159 aj12(i)=p11*y12(i)+p13*y34(i)+p15*y56(i)+p17*y78(i)
160 aj13(i)=p11*z12(i)+p13*z34(i)+p15*z56(i)+p17*z78(i)
161 aj21(i)=p21*x14(i)+p22*x23(i)+p25*x58(i)+p26*x67(i)
162 aj22(i)=p21*y14(i)+p22*y23(i)+p25*y58(i)+p26*y67(i)
163 aj23(i)=p21*z14(i)+p22*z23(i)+p25*z58(i)+p26*z67(i)
164 aj31(i)=p31*x15(i)+p32*x26(i)+p33*x37(i)+p34*x48(i)
165 aj32(i)=p31*y15(i)+p32*y26(i)+p33*y37(i)+p34*y48(i)
166 aj33(i)=p31*z15(i)+p32*z26(i)+p33*z37(i)+p34*z48(i)
167 ENDDO
168
169 DO i=1,nel
170 ai11(i)= aj22(i)*aj33(i)-aj23(i)*aj32(i)
171 ai21(i)=-aj21(i)*aj33(i)+aj23(i)*aj31(i)
172 ai31(i)= aj21(i)*aj32(i)-aj22(i)*aj31(i)
173 vlinc(i,jpt)=aj11(i)*ai11(i)+aj12(i)*ai21(i)+aj13(i)*ai31(i)
174 ENDDO
175
176 DO i=1,nel
177 IF(vlinc(i,jpt)>zero) cycle
179 . msgtype=msgerror,
180 . anmode=aninfo,
181 . i1=ngl(i))
182 ENDDO
183
184 DO i=1,nel
185 vlinv(i)=one/vlinc(i,jpt)
186 ai11(i)=vlinv(i)*ai11(i)
187 ai21(i)=vlinv(i)*ai21(i)
188 ai31(i)=vlinv(i)*ai31(i)
189 ai12(i)=vlinv(i)*(-aj12(i)*aj33(i)+aj13(i)*aj32(i))
190 ai22(i)=vlinv(i)*( aj11(i)*aj33(i)-aj13(i)*aj31(i))
191 ai32(i)=vlinv(i)*(-aj11(i)*aj32(i)+aj12(i)*aj31(i))
192 ai13(i)=vlinv(i)*( aj12(i)*aj23(i)-aj13(i)*aj22(i))
193 ai23(i)=vlinv(i)*(-aj11(i)*aj23(i)+aj13(i)*aj21(i))
194 ai33(i)=vlinv(i)*( aj11(i)*aj22(i)-aj12(i)*aj21(i))
195 ENDDO
196
197 DO i=1,nel
198 a1111(i)=ai11(i)*p11
199 a1113(i)=ai11(i)*p13
200 a1115(i)=ai11(i)*p15
201 a1117(i)=ai11(i)*p17
202 a2111(i)=ai21(i)*p11
203 a2113(i)=ai21(i)*p13
204 a2115(i)=ai21(i)*p15
205 a2117(i)=ai21(i)*p17
206 a3111(i)=ai31(i)*p11
207 a3113(i)=ai31(i)*p13
208 a3115(i)=ai31(i)*p15
209 a3117(i)=ai31(i)*p17
210 a1221(i)=ai12(i)*p21
211 a1222(i)=ai12(i)*p22
212 a1225(i)=ai12(i)*p25
213 a1226(i)=ai12(i)*p26
214 a2221(i)=ai22(i)*p21
215 a2222(i)=ai22(i)*p22
216 ENDDO
217
218 DO i=1,nel
219 a2225(i)=ai22(i)*p25
220 a2226(i)=ai22(i)*p26
221 a3221(i)=ai32(i)*p21
222 a3222(i)=ai32(i)*p22
223 a3225(i)=ai32(i)*p25
224 a3226(i)=ai32(i)*p26
225 a1331(i)=ai13(i)*p31
226 a1332(i)=ai13(i)*p32
227 a1333(i)=ai13(i)*p33
228 a1334(i)=ai13(i)*p34
229 a2331(i)=ai23(i)*p31
230 a2332(i)=ai23(i)*p32
231 a2333(i)=ai23(i)*p33
232 a2334(i)=ai23(i)*p34
233 a3331(i)=ai33(i)*p31
234 a3332(i)=ai33(i)*p32
235 a3333(i)=ai33(i)*p33
236 a3334(i)=ai33(i)*p34
237 ENDDO
238
239 DO i=1,nel
240 px1(i)= a1111(i)+a1221(i)+a1331(i)
241 px2(i)=-a1111(i)+a1222(i)+a1332(i)
242 px3(i)= a1113(i)-a1222(i)+a1333(i)
243 px4(i)=-a1113(i)-a1221(i)+a1334(i)
244 px5(i)= a1115(i)+a1225(i)-a1331(i)
245 px6(i)=-a1115(i)+a1226(i)-a1332(i)
246 px7(i)= a1117(i)-a1226(i)-a1333(i)
247 px8(i)=-a1117(i)-a1225(i)-a1334(i)
248 py1(i)= a2111(i)+a2221(i)+a2331(i)
249 py2(i)=-a2111(i)+a2222(i)+a2332(i)
250 py3(i)= a2113(i)-a2222(i)+a2333(i)
251 py4(i)=-a2113(i)-a2221(i)+a2334(i)
252 py5(i)= a2115(i)+a2225(i)-a2331(i)
253 py6(i)=-a2115(i)+a2226(i)-a2332(i)
254 py7(i)= a2117(i)-a2226(i)-a2333(i)
255 py8(i)=-a2117(i)-a2225(i)-a2334(i)
256 pz1(i)= a3111(i)+a3221(i)+a3331(i)
257 pz2(i)=-a3111(i)+a3222(i)+a3332(i)
258 pz3(i)= a3113(i)-a3222(i)+a3333(i)
259 pz4(i)=-a3113(i)-a3221(i)+a3334(i)
260 pz5(i)= a3115(i)+a3225(i)-a3331(i)
261 pz6(i)=-a3115(i)+a3226(i)-a3332(i)
262 pz7(i)= a3117(i)-a3226(i)-a3333(i)
263 pz8(i)=-a3117(i)-a3225(i)-a3334(i)
264 ENDDO
265
266 DO i=1,nel
267 det(i)=det(i)+vlinc(i,jpt)
268 veul(1,i)=px1(i)*vlinc(i,jpt)+veul(1,i)
269 veul(2,i)=px2(i)*vlinc(i,jpt)+veul(2,i)
270 veul(3,i)=px3(i)*vlinc(i,jpt)+veul(3,i)
271 veul(4,i)=px4(i)*vlinc(i,jpt)+veul(4,i)
272 veul(5,i)=px5(i)*vlinc(i,jpt)+veul(5,i)
273 veul(6,i)=px6(i)*vlinc(i,jpt)+veul(6,i)
274 veul(7,i)=px7(i)*vlinc(i,jpt)+veul(7,i)
275 veul(8,i)=px8(i)*vlinc(i,jpt)+veul(8,i)
276 veul(9,i)=py1(i)*vlinc(i,jpt)+veul(9,i)
277 veul(10,i)=py2(i)*vlinc(i,jpt)+veul(10,i)
278 veul(11,i)=py3(i)*vlinc(i,jpt)+veul(11,i)
279 veul(12,i)=py4(i)*vlinc(i,jpt)+veul(12,i)
280 veul(32,i)=py5(i)*vlinc(i,jpt)+veul(32,i)
281 veul(33,i)=py6(i)*vlinc(i,jpt)+veul(33,i)
282 veul(34,i)=py7(i)*vlinc(i,jpt)+veul(34,i)
283 veul(35,i)=py8(i)*vlinc(i,jpt)+veul(35,i)
284 veul(36,i)=pz1(i)*vlinc(i,jpt)+veul(36,i)
285 veul(37,i)=pz2(i)*vlinc(i,jpt)+veul(37,i)
286 veul(38,i)=pz3(i)*vlinc(i,jpt)+veul(38,i)
287 veul(39,i)=pz4(i)*vlinc(i,jpt)+veul(39,i)
288 veul(40,i)=pz5(i)*vlinc(i,jpt)+veul(40,i)
289 veul(41,i)=pz6(i)*vlinc(i,jpt)+veul(41,i)
290 veul(42,i)=pz7(i)*vlinc(i,jpt)+veul(42,i)
291 veul(43,i)=pz8(i)*vlinc(i,jpt)+veul(43,i)
292 veul(44,i)=veul(44,i)+vlinc(i,jpt)*h(1)
293 veul(45,i)=veul(45,i)+vlinc(i,jpt)*h(2)
294 veul(46,i)=veul(46,i)+vlinc(i,jpt)*h(3)
295 veul(47,i)=veul(47,i)+vlinc(i,jpt)*h(4)
296 veul(48,i)=veul(48,i)+vlinc(i,jpt)*h(5)
297 veul(49,i)=veul(49,i)+vlinc(i,jpt)*h(6)
298 veul(50,i)=veul(50,i)+vlinc(i,jpt)*h(7)
299 veul(51,i)=veul(51,i)+vlinc(i,jpt)*h(8)
300 ENDDO
301
302 ENDDO
303
304 DO i=1,nel
305 xx1 = x1(i) + x2(i) + x3(i) + x4(i)
306 . - x5(i) - x6(i) - x7(i) - x8(i)
307 yy1 = y1(i) + y2(i) + y3(i) + y4(i)
308 . - y5(i) - y6(i) - y7(i) - y8(i)
309 zz1 = z1(i) + z2(i) + z3(i) + z4(i)
310 . - z5(i) - z6(i) - z7(i) - z8(i)
311 xx2 = x1(i) + x2(i) + x5(i) + x6(i)
312 . - x3(i) - x4(i) - x7(i) - x8(i)
313 yy2 = y1(i) + y2(i) + y5(i) + y6(i)
314 . - y3(i) - y4(i) - y7(i) - y8(i)
315 zz2 = z1(i) + z2(i) + z5(i) + z6(i)
316 . - z3(i) - z4(i) - z7(i) - z8(i)
317 xx3 = x1(i) + x4(i) + x5(i) + x8(i)
318 . - x3(i) - x2(i) - x7(i) - x6(i)
319 yy3 = y1(i) + y4(i) + y5(i) + y8(i)
320 . - y3(i) - y2(i) - y7(i) - y6(i)
321 zz3 = z1(i) + z4(i) + z5(i) + z8(i)
322 . - z3(i) - z2(i) - z7(i) - z6(i)
323
324 smax = (yy1 * zz2 - yy2 * zz1)**2
325 . + (zz1 * xx2 - zz2 * xx1)**2
326 . + (xx1 * yy2 - xx2 * yy1)**2
327 smax =
max(smax,(yy1 * zz3 - yy3 * zz1)**2
328 . + (zz1 * xx3 - zz3 * xx1)**2
329 . + (xx1 * yy3 - xx3 * yy1)**2)
330 smax =
max(smax,(yy3 * zz2 - yy2 * zz3)**2
331 . + (zz3 * xx2 - zz2 * xx3)**2
332 . + (xx3 * yy2 - xx2 * yy3)**2)
333 vmin =
min(vlinc(i,1),vlinc(i,2),vlinc(i,3),vlinc(i,4),
334 . vlinc(i,5),vlinc(i,6),vlinc(i,7),vlinc(i,8))
335 deltax(i)=hundred28*vmin/sqrt(smax)
336 ENDDO
337
338 DO j=1,12
339 DO i=1,nel
340 veul(j,i)=veul(j,i)/det(i)
341 ENDDO
342 ENDDO
343
344 DO j=32,43
345 DO i=1,nel
346 veul(j,i)=veul(j,i)/det(i)
347 ENDDO
348 ENDDO
349
350 DO i=1,nel
351 vol(i)=det(i)
352 ENDDO
353
354 IF(jeul /= 0)THEN
355 DO i=1,nel
356 veul(52,i) = vol(i)
357 ENDDO
358 ENDIF
359
360 IF (nint(geo(12,ngeo(1))) == 15) THEN
361 DO i=1,nel
363 veul(44,i)=veul(44,i)*
poro
364 veul(45,i)=veul(45,i)*
poro
365 veul(46,i)=veul(46,i)*
poro
366 veul(47,i)=veul(47,i)*
poro
367 veul(48,i)=veul(48,i)*
poro
368 veul(49,i)=veul(49,i)*
poro
369 veul(50,i)=veul(50,i)*
poro
370 veul(51,i)=veul(51,i)*
poro
372 ENDDO
373 ENDIF
374
375 RETURN
subroutine poro(geo, nodpor, ms, x, v, w, af, am, skew, weight, nporgeo)
subroutine basisf(h, p1, p2, p3, ipt)
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
1.15.0