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