43
44
45
46 USE elbufdef_mod
47
48
49
50#include "implicit_f.inc"
51
52
53
54#include "mvsiz_p.inc"
55
56
57
58 INTEGER JFT,JLT,IREP,NLAY,ISH3NFRAM,NEL
60 . dir_a(nel,*), dir_b(nel,*),e1x0(*), e1y0(*), e1z0(*),
61 . e2x0(*), e2y0(*), e2z0(*),e3x0(*), e3y0(*), e3z0(*),offg(*),
62 . e1x(mvsiz), e1y(mvsiz), e1z(mvsiz),
63 . e2x(mvsiz), e2y(mvsiz), e2z(mvsiz),
64 . e3x(mvsiz), e3y(mvsiz), e3z(mvsiz),
65 . x21(mvsiz), y21(mvsiz), z21(mvsiz),
66 . x31(mvsiz), y31(mvsiz), z31(mvsiz),
67 . ecos(*) ,esin(*),
area(mvsiz)
68
69 REAL(kind=8), dimension(mvsiz), INTENT(in) ::x1,x2,x3
70 REAL(kind=8), dimension(mvsiz), INTENT(in) ::y1,y2,y3
71 REAL(kind=8), dimension(mvsiz), INTENT(in) ::z1,z2,z3
72
73
74 TYPE (ELBUF_STRUCT_) :: ELBUF_STR
75
76
77
78 INTEGER I,II,J,N,ILAW, IDRAPE,IGTYP,NPTT,IT,IPT_ALL,IPT
79
80 my_real x32(mvsiz), y32(mvsiz), z32(mvsiz),sum(mvsiz),
81 . e11(mvsiz),e12(mvsiz),e13(mvsiz),
82 . e21(mvsiz),e22(mvsiz),e23(mvsiz)
83 my_real v1,v2,v3,vr,vs,aa,bb,suma
84 my_real,
DIMENSION(:) ,
POINTER :: dir1, dir2
85
86 idrape = elbuf_str%IDRAPE
87 igtyp = elbuf_str%IGTYP
88 DO i=jft,jlt
89 x21(i)=x2(i)-x1(i)
90 y21(i)=y2(i)-y1(i)
91 z21(i)=z2(i)-z1(i)
92 x31(i)=x3(i)-x1(i)
93 y31(i)=y3(i)-y1(i)
94 z31(i)=z3(i)-z1(i)
95 x32(i)=x3(i)-x2(i)
96 y32(i)=y3(i)-y2(i)
97 z32(i)=z3(i)-z2(i)
98 ENDDO
99
100 IF (irep > 0) THEN
101 DO i=jft,jlt
102 e11(i) = x21(i)
103 e12(i) = y21(i)
104 e13(i) = z21(i)
105 e21(i) = x31(i)
106 e22(i) = y31(i)
107 e23(i) = z31(i)
108 ENDDO
109 ENDIF
110
111 DO i=jft,jlt
112 e1x(i)= x21(i)
113 e1y(i)= y21(i)
114 e1z(i)= z21(i)
115 sum(i) = sqrt(e1x(i)*e1x(i)+e1y(i)*e1y(i)+e1z(i)*e1z(i))
116 e1x(i)=e1x(i)/sum(i)
117 e1y(i)=e1y(i)/sum(i)
118 e1z(i)=e1z(i)/sum(i)
119 ENDDO
120
121 DO i=jft,jlt
122 e3x(i)=y31(i)*z32(i)-z31(i)*y32(i)
123 e3y(i)=z31(i)*x32(i)-x31(i)*z32(i)
124 e3z(i)=x31(i)*y32(i)-y31(i)*x32(i)
125 sum(i) = sqrt(e3x(i)*e3x(i)+e3y(i)*e3y(i)+e3z(i)*e3z(i))
126 e3x(i)=e3x(i)/sum(i)
127 e3y(i)=e3y(i)/sum(i)
128 e3z(i)=e3z(i)/sum(i)
129 area(i) = half * sum(i)
130 ENDDO
131
132 DO i=jft,jlt
133 e2x(i)=e3y(i)*e1z(i)-e3z(i)*e1y(i)
134 e2y(i)=e3z(i)*e1x(i)-e3x(i)*e1z(i)
135 e2z(i)=e3x(i)*e1y(i)-e3y(i)*e1x(i)
136 sum(i) = sqrt(e2x(i)*e2x(i)+e2y(i)*e2y(i)+e2z(i)*e2z(i))
137 e2x(i)=e2x(i)/sum(i)
138 e2y(i)=e2y(i)/sum(i)
139 e2z(i)=e2z(i)/sum(i)
140 ENDDO
141
142 DO i=jft,jlt
143 e1x0(i) = e1x(i)
144 e1y0(i) = e1y(i)
145 e1z0(i) = e1z(i)
146 e2x0(i) = e2x(i)
147 e2y0(i) = e2y(i)
148 e2z0(i) = e2z(i)
149 e3x0(i) = e3x(i)
150 e3y0(i) = e3y(i)
151 e3z0(i) = e3z(i)
152 ENDDO
153 IF(ish3nfram ==0 ) THEN
154 ii = 0
156 . x21, y21, z21,
157 . x31, y31, z31,
158 . e1x,e2x,e3x,e1y,e2y,e3y,e1z,e2z,e3z,sum,offg)
159 DO i=jft,jlt
160 ecos(i) = e1x(i)*e1x0(i)+e1y(i)*e1y0(i)+e1z(i)*e1z0(i)
161 aa =
max(zero,one-ecos(i)*ecos(i))
162 esin(i) = sqrt(aa)
163 bb = e2x(i)*e1x0(i)+e2y(i)*e1y0(i)+e2z(i)*e1z0(i)
164 IF (bb >zero) esin(i) = -esin(i)
165 ENDDO
166 END IF
167
168 IF(idrape > 0 .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
169 ipt_all = 0
170 IF (irep == 1) THEN
171 DO n=1,nlay
172 nptt = elbuf_str%BUFLY(n)%NPTT
173
174 DO it = 1,nptt
175 dir1 => elbuf_str%BUFLY(n)%LBUF_DIR(it)%DIRA
176 ipt = ipt_all + it
177 j = 2*(ipt - 1)
178 DO i=jft,jlt
179 aa = dir1(i)
180 bb = dir1(i + nel)
181 v1 = aa*e11(i) + bb*e21(i)
182 v2 = aa*e12(i) + bb*e22(i)
183 v3 = aa*e13(i) + bb*e23(i)
184 vr = v1*e1x(i) + v2*e1y(i) + v3*e1z(i)
185 vs = v1*e2x(i) + v2*e2y(i) + v3*e2z(i)
186 suma=sqrt(vr*vr + vs*vs)
187 dir_a(i,j+1) = vr/suma
188 dir_a(i,j+2) = vs/suma
189 ENDDO
190 ENDDO
191 ipt_all = ipt_all + nptt
192 ENDDO
193 ELSEIF (irep == 2) THEN
194 DO n=1,nlay
195 nptt = elbuf_str%BUFLY(n)%NPTT
196
197 DO it = 1,nptt
198 dir1 => elbuf_str%BUFLY(n)%LBUF_DIR(it)%DIRA
199 dir2 => elbuf_str%BUFLY(n)%LBUF_DIR(it)%DIRB
200 ipt = ipt_all + it
201 j = 2*(ipt - 1)
202 DO i=jft,jlt
203
204 aa = dir1(i)
205 bb = dir1(i + nel)
206 v1 = aa*e11(i) + bb*e21(i)
207 v2 = aa*e12(i) + bb*e22(i)
208 v3 = aa*e13(i) + bb*e23(i)
209 vr = v1*e1x(i) + v2*e1y(i) + v3*e1z(i)
210 vs = v1*e2x(i) + v2*e2y(i) + v3*e2z(i)
211 suma = one /
max( sqrt(vr*vr + vs*vs), em20)
212 dir_a(i,j+1) = vr*suma
213 dir_a(i,j+2) = vs*suma
214
215 aa = dir2(i)
216 bb = dir2(i + nel)
217 v1 = aa*e11(i) + bb*e21(i)
218 v2 = aa*e12(i) + bb*e22(i)
219 v3 = aa*e13(i) + bb*e23(i)
220 vr = v1*e1x(i) + v2*e1y(i) + v3*e1z(i)
221 vs = v1*e2x(i) + v2*e2y(i) + v3*e2z(i)
222 suma = one /
max( sqrt(vr*vr + vs*vs), em20)
223 dir_b(i,j+1) = vr*suma
224 dir_b(i,j+2) = vs*suma
225 ENDDO
226 ENDDO
227 ipt_all = ipt_all + nptt
228 ENDDO
229 ELSEIF (irep == 3) THEN
230
231 DO n=1,nlay
232 ilaw = elbuf_str%BUFLY(n)%ILAW
233 nptt = elbuf_str%BUFLY(n)%NPTT
234 IF (ilaw == 58) THEN
235
236 DO it = 1,nptt
237 dir1 => elbuf_str%BUFLY(n)%LBUF_DIR(it)%DIRA
238 dir2 => elbuf_str%BUFLY(n)%LBUF_DIR(it)%DIRB
239 ipt = ipt_all + it
240 j = 2*(ipt - 1)
241 DO i=jft,jlt
242
243 aa = dir1(i)
244 bb = dir1(i + nel)
245 v1 = aa*e11(i) + bb*e21(i)
246 v2 = aa*e12(i) + bb*e22(i)
247 v3 = aa*e13(i) + bb*e23(i)
248 vr = v1*e1x(i)+ v2*e1y(i) + v3*e1z(i)
249 vs = v1*e2x(i)+ v2*e2y(i) + v3*e2z(i)
250 suma = one /
max( sqrt(vr*vr + vs*vs), em20)
251 dir_a(i,j+1) = vr*suma
252 dir_a(i,j+2) = vs*suma
253
254 aa = dir2(i)
255 bb = dir2(i + nel)
256 v1 = aa*e11(i) + bb*e21(i)
257 v2 = aa*e12(i) + bb*e22(i)
258 v3 = aa*e13(i) + bb*e23(i)
259 vr = v1*e1x(i)+ v2*e1y(i) + v3*e1z(i)
260 vs = v1*e2x(i)+ v2*e2y(i) + v3*e2z(i)
261 suma = one /
max( sqrt(vr*vr + vs*vs), em20)
262 dir_b(i,j+1) = vr*suma
263 dir_b(i,j+2) = vs*suma
264 ENDDO
265 ENDDO
266 ipt_all = ipt_all + nptt
267 ELSE
268 DO it = 1,nptt
269 dir1 => elbuf_str%BUFLY(n)%LBUF_DIR(it)%DIRA
270 ipt = ipt_all + it
271 j = 2*(ipt - 1)
272 DO i=jft,jlt
273 dir_a(i,j+1) = dir1(i)
274 dir_a(i,j+2) = dir1(i+nel)
275 ENDDO
276 ENDDO
277 ipt_all = ipt_all + nptt
278 ENDIF
279 ENDDO
280 ELSEIF (irep == 4) THEN
281
282 DO n=1,nlay
283 ilaw = elbuf_str%BUFLY(n)%ILAW
284 nptt = elbuf_str%BUFLY(n)%NPTT
285
286 IF (ilaw == 58) THEN
287 DO it = 1,nptt
288 dir1 => elbuf_str%BUFLY(n)%LBUF_DIR(it)%DIRA
289 dir2 => elbuf_str%BUFLY(n)%LBUF_DIR(it)%DIRB
290 ipt = ipt_all + it
291 j = 2*(ipt - 1)
292 DO i=jft,jlt
293
294 aa = dir1(i)
295 bb = dir1(i + nel)
296 v1 = aa*e11(i) + bb*e21(i)
297 v2 = aa*e12(i) + bb*e22(i)
298 v3 = aa*e13(i) + bb*e23(i)
299 vr = v1*e1x(i)+ v2*e1y(i) + v3*e1z(i)
300 vs = v1*e2x(i)+ v2*e2y(i) + v3*e2z(i)
301 suma = one /
max( sqrt(vr*vr + vs*vs), em20)
302 dir_a(i,j+1) = vr*suma
303 dir_a(i,j+2) = vs*suma
304
305 aa = dir2(i)
306 bb = dir2(i + nel)
307 v1 = aa*e11(i) + bb*e21(i)
308 v2 = aa*e12(i) + bb*e22(i)
309 v3 = aa*e13(i) + bb*e23(i)
310 vr = v1*e1x(i)+ v2*e1y(i) + v3*e1z(i)
311 vs = v1*e2x(i)+ v2*e2y(i) + v3*e2z(i)
312 suma = one /
max( sqrt(vr*vr + vs*vs), em20)
313 dir_b(i,j+1) = vr*suma
314 dir_b(i,j+2) = vs*suma
315 ENDDO
316 ENDDO
317 ipt_all = ipt_all + nptt
318 ELSE
319
320 DO it = 1,nptt
321 dir1 => elbuf_str%BUFLY(n)%LBUF_DIR(it)%DIRA
322 ipt = ipt_all + it
323 j = 2*(ipt - 1)
324 DO i=jft,jlt
325 aa = dir1(i)
326 bb = dir1(i + nel)
327 v1 = aa*e11(i) + bb*e21(i)
328 v2 = aa*e12(i) + bb*e22(i)
329 v3 = aa*e13(i) + bb*e23(i)
330 vr = v1*e1x(i)+ v2*e1y(i) + v3*e1z(i)
331 vs = v1*e2x(i)+ v2*e2y(i) + v3*e2z(i)
332 suma=sqrt(vr*vr + vs*vs)
333 dir_a(i,j+1) = vr/suma
334 dir_a(i,j+2) = vs/suma
335 ENDDO
336 ENDDO
337 ipt_all = ipt_all + nptt
338 ENDIF
339 ENDDO
340 ENDIF
341 ELSE
342 IF (irep == 1) THEN
343 DO n=1,nlay
344 dir1 => elbuf_str%BUFLY(n)%DIRA
345 j = 2*(n-1)
346 DO i=jft,jlt
347 aa = dir1(i)
348 bb = dir1(i + nel)
349 v1 = aa*e11(i) + bb*e21(i)
350 v2 = aa*e12(i) + bb*e22(i)
351 v3 = aa*e13(i) + bb*e23(i)
352 vr = v1*e1x(i) + v2*e1y(i) + v3*e1z(i)
353 vs = v1*e2x(i) + v2*e2y(i) + v3*e2z(i)
354 suma=sqrt(vr*vr + vs*vs)
355 dir_a(i,j+1) = vr/suma
356 dir_a(i,j+2) = vs/suma
357 ENDDO
358 ENDDO
359 ELSEIF (irep == 2) THEN
360 DO n=1,nlay
361 dir1 => elbuf_str%BUFLY(n)%DIRA
362 dir2 => elbuf_str%BUFLY(n)%DIRB
363 j = (n-1)*nel
364 j = 2*(n-1)
365 DO i=jft,jlt
366
367 aa = dir1(i)
368 bb = dir1(i + nel)
369 v1 = aa*e11(i) + bb*e21(i)
370 v2 = aa*e12(i) + bb*e22(i)
371 v3 = aa*e13(i) + bb*e23(i)
372 vr = v1*e1x(i) + v2*e1y(i) + v3*e1z(i)
373 vs = v1*e2x(i) + v2*e2y(i) + v3*e2z(i)
374 suma = one /
max( sqrt(vr*vr + vs*vs), em20)
375 dir_a(i,j+1) = vr*suma
376 dir_a(i,j+2) = vs*suma
377
378 aa = dir2(i)
379 bb = dir2(i + nel)
380 v1 = aa*e11(i) + bb*e21(i)
381 v2 = aa*e12(i) + bb*e22(i)
382 v3 = aa*e13(i) + bb*e23(i)
383 vr = v1*e1x(i) + v2*e1y(i) + v3*e1z(i)
384 vs = v1*e2x(i) + v2*e2y(i) + v3*e2z(i)
385 suma = one /
max( sqrt(vr*vr + vs*vs), em20)
386 dir_b(i,j+1) = vr*suma
387 dir_b(i,j+2) = vs*suma
388 ENDDO
389 ENDDO
390 ELSEIF (irep == 3) THEN
391
392 DO n=1,nlay
393 ilaw = elbuf_str%BUFLY(n)%ILAW
394 j = (n-1)*nel
395 j = 2*(n-1)
396 IF (ilaw == 58) THEN
397 dir1 => elbuf_str%BUFLY(n)%DIRA
398 dir2 => elbuf_str%BUFLY(n)%DIRB
399 DO i=jft,jlt
400
401 aa = dir1(i)
402 bb = dir1(i + nel)
403 v1 = aa*e11(i) + bb*e21(i)
404 v2 = aa*e12(i) + bb*e22(i)
405 v3 = aa*e13(i) + bb*e23(i)
406 vr = v1*e1x(i)+ v2*e1y(i) + v3*e1z(i)
407 vs = v1*e2x(i)+ v2*e2y(i) + v3*e2z(i)
408 suma = one /
max( sqrt(vr*vr + vs*vs), em20)
409 dir_a(i,j+1) = vr*suma
410 dir_a(i,j+2) = vs*suma
411
412 aa = dir2(i)
413 bb = dir2(i + nel)
414 v1 = aa*e11(i) + bb*e21(i)
415 v2 = aa*e12(i) + bb*e22(i)
416 v3 = aa*e13(i) + bb*e23(i)
417 vr = v1*e1x(i)+ v2*e1y(i) + v3*e1z(i)
418 vs = v1*e2x(i)+ v2*e2y(i) + v3*e2z(i)
419 suma = one /
max( sqrt(vr*vr + vs*vs), em20)
420 dir_b(i,j+1) = vr*suma
421 dir_b(i,j+2) = vs*suma
422 ENDDO
423 ELSE
424 dir1 => elbuf_str%BUFLY(n)%DIRA
425 DO i=jft,jlt
426 dir_a(i,j+1) = dir1(i)
427 dir_a(i,j+2) = dir1(i+nel)
428 ENDDO
429 ENDIF
430 ENDDO
431 ELSEIF (irep == 4) THEN
432
433 DO n=1,nlay
434 ilaw = elbuf_str%BUFLY(n)%ILAW
435 j = (n-1)*nel
436 j = 2*(n-1)
437 IF (ilaw == 58) THEN
438 dir1 => elbuf_str%BUFLY(n)%DIRA
439 dir2 => elbuf_str%BUFLY(n)%DIRB
440 DO i=jft,jlt
441
442 aa = dir1(i)
443 bb = dir1(i + nel)
444 v1 = aa*e11(i) + bb*e21(i)
445 v2 = aa*e12(i) + bb*e22(i)
446 v3 = aa*e13(i) + bb*e23(i)
447 vr = v1*e1x(i)+ v2*e1y(i) + v3*e1z(i)
448 vs = v1*e2x(i)+ v2*e2y(i) + v3*e2z(i)
449 suma = one /
max( sqrt(vr*vr + vs*vs), em20)
450 dir_a(i,j+1) = vr*suma
451 dir_a(i,j+2) = vs*suma
452
453 aa = dir2(i)
454 bb = dir2(i + nel)
455 v1 = aa*e11(i) + bb*e21(i)
456 v2 = aa*e12(i) + bb*e22(i)
457 v3 = aa*e13(i) + bb*e23(i)
458 vr = v1*e1x(i)+ v2*e1y(i) + v3*e1z(i)
459 vs = v1*e2x(i)+ v2*e2y(i) + v3*e2z(i)
460 suma = one /
max( sqrt(vr*vr + vs*vs), em20)
461 dir_b(i,j+1) = vr*suma
462 dir_b(i,j+2) = vs*suma
463 ENDDO
464 ELSE
465 dir1 => elbuf_str%BUFLY(n)%DIRA
466 DO i=jft,jlt
467 aa = dir1(i)
468 bb = dir1(i + nel)
469 v1 = aa*e11(i) + bb*e21(i)
470 v2 = aa*e12(i) + bb*e22(i)
471 v3 = aa*e13(i) + bb*e23(i)
472 vr = v1*e1x(i)+ v2*e1y(i) + v3*e1z(i)
473 vs = v1*e2x(i)+ v2*e2y(i) + v3*e2z(i)
474 suma=sqrt(vr*vr + vs*vs)
475 dir_a(i,j+1) = vr/suma
476 dir_a(i,j+2) = vs/suma
477 ENDDO
478 ENDIF
479 ENDDO
480 ENDIF
481 ENDIF
482
483 IF(ish3nfram ==0 ) THEN
484 DO i=jft,jlt
485 e11(i) = e1x(i)
486 e12(i) = e1y(i)
487 e13(i) = e1z(i)
488 e21(i) = e2x(i)
489 e22(i) = e2y(i)
490 e23(i) = e2z(i)
491 ENDDO
492 DO i=jft,jlt
493 e1x(i) = e1x0(i)
494 e1y(i) = e1y0(i)
495 e1z(i) = e1z0(i)
496 e2x(i) = e2x0(i)
497 e2y(i) = e2y0(i)
498 e2z(i) = e2z0(i)
499 e3x(i) = e3x0(i)
500 e3y(i) = e3y0(i)
501 e3z(i) = e3z0(i)
502 ENDDO
503
504 DO i=jft,jlt
505 e1x0(i) = e11(i)
506 e1y0(i) = e12(i)
507 e1z0(i) = e13(i)
508 e2x0(i) = e21(i)
509 e2y0(i) = e22(i)
510 e2z0(i) = e23(i)
511 ENDDO
512 END IF
513
514 RETURN
subroutine clskew3(jft, jlt, irep, rx, ry, rz, sx, sy, sz, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z, det)
subroutine area(d1, x, x2, y, y2, eint, stif0)