45
46
47
48 USE elbufdef_mod
49 use element_mod , only : nixtg
50
51
52
53#include "implicit_f.inc"
54
55
56
57#include "mvsiz_p.inc"
58
59
60
61#include "param_c.inc"
62#include "com08_c.inc"
63#include "scr17_c.inc"
64
65
66
67 INTEGER JFT, JLT,NLAY,ISMSTR,IREP,NEL
68 INTEGER IXTG(NIXTG,*),MAT(*),PID(*),NGL(*),IXTG1(4,*),
69 . NVS,IVS(*),IGEO(NPROPGI,*)
70
72 . x(3,*),v(3,*),r(3,*), offg(*), off(*),
73 . r11(*),r12(*),r13(*),r21(*),r22(*),r23(*),
74 . r31(*),r32(*),r33(*),
area(*),area2(*),
75 . n4x(*),n4y(*),n4z(*),n5x(*),n5y(*),n5z(*),n6x(*),n6y(*),
76 . n6z(*),area4(*),area5(*),area6(*),
77 . vlx(mvsiz,2),vly(mvsiz,2),vlz(mvsiz,6),vnz4(mvsiz,3),vnz5(mvsiz,3),vnz6(mvsiz,3),
78 . xl2(*),xl3(*),yl2(*),yl3(*),xl4(*),xl5(*),yl4(*),yl5(*),
79 . xl6(*),yl6(*),zl4(*),zl5(*),zl6(*),
80 . dir_a(nel,*),dir_b(nel,*)
81 TYPE(ELBUF_STRUCT_) :: ELBUF_STR
82 double precision
83 . smstr(*)
84
85
86
87 INTEGER NC1, NC2, NC3, NC4(MVSIZ), NC5(MVSIZ), NC6(MVSIZ),
88 . I,II(13),I1 ,SVS,JVS(MVSIZ),EP, MAT_1
89
91 . vx1(mvsiz), vx2(mvsiz), vx3(mvsiz),
92 . vy1(mvsiz), vy2(mvsiz), vy3(mvsiz),
93 . vz1(mvsiz), vz2(mvsiz), vz3(mvsiz),
94 .
95 .
96 . vx4(mvsiz), vy4(mvsiz), vz4(mvsiz),
97 . vx5(mvsiz), vy5(mvsiz), vz5(mvsiz),
98 . vx6(mvsiz), vy6(mvsiz), vz6(mvsiz),
99 .
100 .
101 . x1(mvsiz),x2(mvsiz),x3(mvsiz),x4(mvsiz),x5(mvsiz),x6(mvsiz),
102 . y1(mvsiz),y2(mvsiz),y3(mvsiz),y4(mvsiz),y5(mvsiz),y6(mvsiz),
103 . z1(mvsiz),z2(mvsiz),z3(mvsiz),z4(mvsiz),z5(mvsiz),z6(mvsiz),
104 . rx(mvsiz), ry(mvsiz), rz(mvsiz),
105 . sx(mvsiz), sy(mvsiz), sz(mvsiz),
106 . vx0, vy0,vz0,off_l,a,b,c ,dt05,exz,eyz,vz21,vz31 ,
107 . ddrx,ddry,v21x,v31x,ddrz1,ddrz2
108
109 DO i=1,13
110 ii(i) = nel*(i-1)
111 ENDDO
112
113 mat_1 = ixtg(1,jft)
114 DO i=jft,jlt
115 mat(i) = mat_1
116 nc1 = ixtg(2,i)
117 nc2 = ixtg(3,i)
118 nc3 = ixtg(4,i)
119 pid(i) = ixtg(5,i)
120 ngl(i) = ixtg(6,i)
121
122
123
124 x1(i)=x(1,nc1)
125 y1(i)=x(2,nc1)
126 z1(i)=x(3,nc1)
127 x2(i)=x(1,nc2)
128 y2(i)=x(2,nc2)
129 z2(i)=x(3,nc2)
130 x3(i)=x(1,nc3)
131 y3(i)=x(2,nc3)
132 z3(i)=x(3,nc3)
133
134
135
136 vx1(i)=v(1,nc1)
137 vy1(i)=v(2,nc1)
138 vz1(i)=v(3,nc1)
139 vx2(i)=v(1,nc2)
140 vy2(i)=v(2,nc2)
141 vz2(i)=v(3,nc2)
142 vx3(i)=v(1,nc3)
143 vy3(i)=v(2,nc3)
144 vz3(i)=v(3,nc3)
145 ENDDO
146
147
148
149 nvs=jft-1
150 svs= nvs
151 DO i=jft,jlt
152 nc4(i) = ixtg1(1,i)
153 nc5(i) = ixtg1(2,i)
154 nc6(i) = ixtg1(3,i)
155 IF (nc4(i) > 0 .AND. nc5(i) > 0 .AND. nc6(i) > 0) THEN
156 nvs=nvs+1
157 ivs(nvs)=i
158 ELSE
159 svs=svs+1
160 jvs(svs)=i
161 ENDIF
162 ENDDO
163 DO i=nvs+1,jlt
164 ivs(i)=jvs(i-nvs)
165 ENDDO
166
167 DO i=jft,jlt
168 rx(i)=x2(i)-x1(i)
169 ry(i)=y2(i)-y1(i)
170 rz(i)=z2(i)-z1(i)
171 sx(i)=x3(i)-x1(i)
172 sy(i)=y3(i)-y1(i)
173 sz(i)=z3(i)-z1(i)
174 ENDDO
175
176
177
178 i1 = 0
180 . rx, ry, rz,
181 . sx, sy, sz,
182 . r11,r12,r13,r21,r22,r23,r31,r32,r33,area2,offg )
183
184 DO i=jft,jlt
185 xl2(i)=r11(i)*rx(i)+r21(i)*ry(i)+r31(i)*rz(i)
186 yl2(i)=r12(i)*rx(i)+r22(i)*ry(i)+r32(i)*rz(i)
187 xl3(i)=r11(i)*sx(i)+r21(i)*sy(i)+r31(i)*sz(i)
188 yl3(i)=r12(i)*sx(i)+r22(i)*sy(i)+r32(i)*sz(i)
189 area(i)=half*area2(i)
190 ENDDO
191
192
193
194 IF (ismstr == 1 .OR. ismstr == 2) THEN
195 DO i=jft,jlt
196 IF(abs(offg(i)) == two)THEN
197 xl2(i)=smstr(ii(1)+i)
198 yl2(i)=smstr(ii(2)+i)
199 xl3(i)=smstr(ii(3)+i)
200 yl3(i)=smstr(ii(4)+i)
201 area2(i)=xl2(i)*yl3(i)-xl3(i)*yl2(i)
202 area(i)=half*area2(i)
203 ELSE
204 smstr(ii(1)+i)=xl2(i)
205 smstr(ii(2)+i)=yl2(i)
206 smstr(ii(3)+i)=xl3(i)
207 smstr(ii(4)+i)=yl3(i)
208 ENDIF
209 ENDDO
210 ENDIF
211 IF (ismstr == 1) THEN
212 DO i=jft,jlt
213 IF (offg(i) == one) offg(i)=two
214 ENDDO
215 ENDIF
216
217
218 DO ep=jft,nvs
219 i =ivs(ep)
220 vx4(i)=v(1,nc4(i))
221 vy4(i)=v(2,nc4(i))
222 vz4(i)=v(3,nc4(i))
223 vx5(i)=v(1,nc5(i))
224 vy5(i)=v(2,nc5(i))
225 vz5(i)=v(3,nc5(i))
226 vx6(i)=v(1,nc6(i))
227 vy6(i)=v(2,nc6(i))
228 vz6(i)=v(3,nc6(i))
229 ENDDO
230 DO ep=jft,nvs
231 i =ivs(ep)
232 x4(i)=x(1,nc4(i))-x1(i)
233 y4(i)=x(2,nc4(i))-y1(i)
234 z4(i)=x(3,nc4(i))-z1(i)
235 x5(i)=x(1,nc5(i))-x2(i)
236 y5(i)=x(2,nc5(i))-y2(i)
237 z5(i)=x(3,nc5(i))-z2(i)
238 x6(i)=x(1,nc6(i))-x3(i)
239 y6(i)=x(2,nc6(i))-y3(i)
240 z6(i)=x(3,nc6(i))-z3(i)
241 ENDDO
242
243 DO ep=nvs+1,jlt
244 i =ivs(ep)
245
246 IF (nc4(i) <= 0) THEN
247 vx4(i)=zero
248 vy4(i)=zero
249 vz4(i)=zero
250 x4(i)=sx(i)
251 y4(i)=sy(i)
252 z4(i)=sz(i)
253 ELSE
254 vx4(i)=v(1,nc4(i))
255 vy4(i)=v(2,nc4(i))
256 vz4(i)=v(3,nc4(i))
257 x4(i)=x(1,nc4(i))-x1(i)
258 y4(i)=x(2,nc4(i))-y1(i)
259 z4(i)=x(3,nc4(i))-z1(i)
260 ENDIF
261 IF (nc5(i) <= 0) THEN
262 vx5(i)=zero
263 vy5(i)=zero
264 vz5(i)=zero
265 x5(i)=x1(i)
266 y5(i)=y1(i)
267 z5(i)=z1(i)
268 ELSE
269 vx5(i)=v(1,nc5(i))
270 vy5(i)=v(2,nc5(i))
271 vz5(i)=v(3,nc5(i))
272 x5(i)=x(1,nc5(i))-x2(i)
273 y5(i)=x(2,nc5(i))-y2(i)
274 z5(i)=x(3,nc5(i))-z2(i)
275 ENDIF
276 IF (nc6(i) <= 0) THEN
277 vx6(i)=zero
278 vy6(i)=zero
279 vz6(i)=zero
280 x6(i)=rx(i)
281 y6(i)=ry(i)
282 z6(i)=rz(i)
283 ELSE
284 vx6(i)=v(1,nc6(i))
285 vy6(i)=v(2,nc6(i))
286 vz6(i)=v(3,nc6(i))
287 x6(i)=x(1,nc6(i))-x3(i)
288 y6(i)=x(2,nc6(i))-y3(i)
289 z6(i)=x(3,nc6(i))-z3(i)
290 ENDIF
291 ENDDO
292
293 DO i=jft,jlt
294 xl4(i)=r11(i)*x4(i)+r21(i)*y4(i)+r31(i)*z4(i)
295 yl4(i)=r12(i)*x4(i)+r22(i)*y4(i)+r32(i)*z4(i)
296 zl4(i)=r13(i)*x4(i)+r23(i)*y4(i)+r33(i)*z4(i)
297 xl5(i)=r11(i)*x5(i)+r21(i)*y5(i)+r31(i)*z5(i)
298 yl5(i)=r12(i)*x5(i)+r22(i)*y5(i)+r32(i)*z5(i)
299 zl5(i)=r13(i)*x5(i)+r23(i)*y5(i)+r33(i)*z5(i)
300 xl6(i)=r11(i)*x6(i)+r21(i)*y6(i)+r31(i)*z6(i)
301 yl6(i)=r12(i)*x6(i)+r22(i)*y6(i)+r32(i)*z6(i)
302 zl6(i)=r13(i)*x6(i)+r23(i)*y6(i)+r33(i)*z6(i)
303 ENDDO
304
305
306
307 IF (ismstr == 1 .OR. ismstr == 2) THEN
308 DO i=jft,jlt
309 IF (abs(offg(i)) == two) THEN
310 xl4(i)=smstr(ii(5)+i)
311 yl4(i)=smstr(ii(6)+i)
312 zl4(i)=smstr(ii(7)+i)
313 xl5(i)=smstr(ii(8)+i)
314 yl5(i)=smstr(ii(9)+i)
315 zl5(i)=smstr(ii(10)+i)
316 xl6(i)=smstr(ii(11)+i)
317 yl6(i)=smstr(ii(12)+i)
318 zl6(i)=smstr(ii(13)+i)
319 ELSE
320 smstr(ii(5)+i) =xl4(i)
321 smstr(ii(6)+i) =yl4(i)
322 smstr(ii(7)+i) =zl4(i)
323 smstr(ii(8)+i) =xl5(i)
324 smstr(ii(9)+i) =yl5(i)
325 smstr(ii(10)+i)=zl5(i)
326 smstr(ii(11)+i)=xl6(i)
327 smstr(ii(12)+i)=yl6(i)
328 smstr(ii(13)+i)=zl6(i)
329 ENDIF
330 ENDDO
331 ENDIF
332
333
334
335 DO i=jft,jlt
336 vx0=r11(i)*vx1(i)+r21(i)*vy1(i)+r31(i)*vz1(i)
337 vy0=r12(i)*vx1(i)+r22(i)*vy1(i)+r32(i)*vz1(i)
338 vz0=r13(i)*vx1(i)+r23(i)*vy1(i)+r33(i)*vz1(i)
339 vx1(i)=vx0
340 vy1(i)=vy0
341 vz1(i)=vz0
342 vx0=r11(i)*vx2(i)+r21(i)*vy2(i)+r31(i)*vz2(i)
343 vy0=r12(i)*vx2(i)+r22(i)*vy2(i)+r32(i)*vz2(i)
344 vz0=r13(i)*vx2(i)+r23(i)*vy2(i)+r33(i)*vz2(i)
345 vx2(i)=vx0
346 vy2(i)=vy0
347 vz2(i)=vz0
348 vx0=r11(i)*vx3(i)+r21(i)*vy3(i)+r31(i)*vz3(i)
349 vy0=r12(i)*vx3(i)+r22(i)*vy3(i)+r32(i)*vz3(i)
350 vz0=r13(i)*vx3(i)+r23(i)*vy3(i)+r33(i)*vz3(i)
351 vx3(i)=vx0
352 vy3(i)=vy0
353 vz3(i)=vz0
354 vx0=r11(i)*vx4(i)+r21(i)*vy4(i)+r31(i)*vz4(i)
355 vy0=r12(i)*vx4(i)+r22(i)*vy4(i)+r32(i)*vz4(i)
356 vz0=r13(i)*vx4(i)+r23(i)*vy4(i)+r33(i)*vz4(i)
357 vx4(i)=vx0
358 vy4(i)=vy0
359 vz4(i)=vz0
360 vx0=r11(i)*vx5(i)+r21(i)*vy5(i)+r31(i)*vz5(i)
361 vy0=r12(i)*vx5(i)+r22(i)*vy5(i)+r32(i)*vz5(i)
362 vz0=r13(i)*vx5(i)+r23(i)*vy5(i)+r33(i)*vz5(i)
363 vx5(i)=vx0
364 vy5(i)=vy0
365 vz5(i)=vz0
366 vx0=r11(i)*vx6(i)+r21(i)*vy6(i)+r31(i)*vz6(i)
367 vy0=r12(i)*vx6(i)+r22(i)*vy6(i)+r32(i)*vz6(i)
368 vz0=r13(i)*vx6(i)+r23(i)*vy6(i)+r33(i)*vz6(i)
369 vx6(i)=vx0
370 vy6(i)=vy0
371 vz6(i)=vz0
372 ENDDO
373
374 DO i=jft,jlt
375
376 a =-yl2(i)*zl4(i)
377 b = xl2(i)*zl4(i)
378 c = xl4(i)*yl2(i)-xl2
379 area4(i)=one/sqrt(a*a+b*b+c*c)
380 n4x(i)=a*area4(i)
381 n4y(i)=b*area4(i)
382 n4z(i)=c*area4(i)
383
384 a =-(yl3(i)-yl2(i))*zl5(i)
385 b = (xl3(i)-xl2(i))*zl5(i)
386 c = xl5(i)*(yl3(i)-yl2(i))-(xl3(i)-xl2(i))*yl5(i)
387 area5(i)=one/sqrt(a*a+b*b+c*c)
388 n5x(i)=a*area5(i)
389 n5y(i)=b*area5(i)
390 n5z(i)=c*area5(i)
391
392 a = yl3(i)*zl6(i)
393 b = -xl3(i)*zl6(i)
394 c = -xl6(i)*yl3(i)+xl3(i)*yl6(i)
395 area6(i)=one/sqrt(a*a+b*b+c*c)
396 n6x(i)=a*area6(i)
397 n6y(i)=b*area6(i)
398 n6z(i)=c*area6(i)
399 ENDDO
400 DO ep=nvs+1,jlt
401 i =ivs(ep)
402 IF (nc4(i) <= 0) THEN
403 n4x(i)=zero
404 n4y(i)=zero
405 n4z(i)=one
406 area4(i)=one/sqrt(area2(i))
407 END IF
408 IF (nc5(i) <= 0) THEN
409 n5x(i)=zero
410 n5y(i)=zero
411 n5z(i)=one
412 area5(i)=one/sqrt(area2(i))
413 END IF
414 IF (nc6(i) <= 0) THEN
415 n6x(i)=zero
416 n6y(i)=zero
417 n6z(i)=one
418 area6(i)=one/sqrt(area2(i))
419 END IF
420 ENDDO
421 DO i=jft,jlt
422 vlz(i,4)=n4x(i)*vx4(i)+n4y(i)*vy4(i)+n4z(i)*vz4(i)
423 vlz(i,5)=n5x(i)*vx5(i)+n5y(i)*vy5(i)+n5z(i)*vz5(i)
424 vlz(i,6)=n6x(i)*vx6(i)+n6y(i)*vy6(i)+n6z(i)*vz6(i)
425
426 vnz4(i,1)=n4x(i)*vx1(i)+n4y(i)*vy1(i)+n4z(i)*vz1(i)
427 vnz4(i,2)=n4x(i)*vx2(i)+n4y(i)*vy2(i)+n4z(i)*vz2(i)
428 vnz4(i,3)=n4x(i)*vx3(i)+n4y(i)*vy3(i)+n4z(i)*vz3(i)
429 vnz5(i,1)=n5x(i)*vx1(i)+n5y(i)*vy1(i)+n5z(i)*vz1(i)
430 vnz5(i,2)=n5x(i)*vx2(i)+n5y(i)*vy2(i)+n5z(i)*vz2(i)
431 vnz5(i,3)=n5x(i)*vx3(i)+n5y(i)*vy3(i)+n5z(i)*vz3(i)
432 vnz6(i,1)=n6x(i)*vx1(i)+n6y(i)*vy1(i)+n6z(i)*vz1(i)
433 vnz6(i,2)=n6x(i)*vx2(i)+n6y(i)*vy2(i)+n6z(i)*vz2(i)
434 vnz6(i,3)=n6x(i)*vx3(i)+n6y(i)*vy3(i)+n6z(i)*vz3(i)
435 ENDDO
436 DO i=jft,jlt
437 vlx(i,1)=vx2(i)-vx1(i)
438 vlx(i,2)=vx3(i)-vx1(i)
439 vly(i,1)=vy2(i)-vy1(i)
440 vly(i,2)=vy3(i)-vy1(i)
441 vlz(i,1)=vz1(i)
442 vlz(i,2)=vz2(i)
443 vlz(i,3)=vz3(i)
444 ENDDO
445
446
447
448 IF (irep > 0) THEN
449 CALL cortdir3(elbuf_str,dir_a,dir_b ,jft ,jlt ,
450 . nlay ,irep ,rx ,ry ,rz ,
451 . sx ,sy ,sz ,r11 ,r21 ,
452 . r31 ,r12 ,r22 ,r32 ,nel )
453 ENDIF
454
455
456
457 dt05 = half*dt1
458 DO i=jft,jlt
459 vz21=vz2(i)-vz1(i)
460 vz31=vz3(i)-vz1(i)
461 exz = yl3(i)*vz21-yl2(i)*vz31
462 eyz = -xl3(i)*vz21+xl2(i)*vz31
463 ddry=dt05*exz/area2(i)
464 ddrx=dt05*eyz/area2(i)
465 v21x = vlx(i,1)
466 v31x = vlx(i,2)
467 ddrz1=dt05*vly(i,1)/xl2(i)
468 ddrz2=dt05*v31x/yl3(i)
469 vlx(i,1) = vlx(i,1)-ddry*vz21-ddrz1*vly(i,1)
470 vlx(i,2) = vlx(i,2)-ddry*vz31-ddrz1*vly(i,2)
471 vly(i,1) = vly(i,1)-ddrx*vz21-ddrz2*v21x
472 vly(i,2) = vly(i,2)-ddrx*vz31-ddrz2*v31x
473 ENDDO
474
475 off_l = zero
476 DO i=jft,jlt
477 off(i) =
min(one,abs(offg(i)))
478 off_l =
min(off_l,offg(i))
479 ENDDO
480 IF(off_l < zero)THEN
481 DO i=jft,jlt
482 IF(offg(i) < zero)THEN
483 vlx(i,1) =zero
484 vlx(i,2) =zero
485 vly(i,1) =zero
486 vly(i,2) =zero
487 vlz(i,1) =zero
488 vlz(i,2) =zero
489 vlz(i,3) =zero
490 vlz(i,4) =zero
491 vlz(i,5) =zero
492 vlz(i,6) =zero
493 vnz4(i,1)=zero
494 vnz4(i,2)=zero
495 vnz4(i,3)=zero
496 vnz5(i,1)=zero
497 vnz5(i,2)=zero
498 vnz5(i,3)=zero
499 vnz6(i,1)=zero
500 vnz6(i,2)=zero
501 vnz6(i,3)=zero
502 ENDIF
503 ENDDO
504 ENDIF
505
506 RETURN
subroutine clskew3(jft, jlt, irep, rx, ry, rz, sx, sy, sz, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z, det)
subroutine cortdir3(elbuf_str, dir_a, dir_b, jft, jlt, nlay, irep, rx, ry, rz, sx, sy, sz, e1x, e1y, e1z, e2x, e2y, e2z, nel)
subroutine area(d1, x, x2, y, y2, eint, stif0)