45
46
47
49
50
51
52#include "implicit_f.inc"
53#include "comlock.inc"
54
55
56
57#include "mvsiz_p.inc"
58
59
60
61 INTEGER NSN, I0,I2SIZE,,IROT, NOINT,NELTST,ITYPTST
62 INTEGER IRECT(4,*),NSV(*),IRTL(*),WEIGHT(*),IADI2(4,*),
63 . NODNX_SMS(*),IVISC
64
66 . visc,dt2t
68 . x(3,*),vr(3,*),v(3,*),a(3,*),ar(3,*),ms(*),crst(2,*),
69 . skew(9,*),dx(3,*),xini(3,*),fini(3,*),fsav(*),fncont(3,*),
70 . stifn(*),stifr(*),stfn(*),stfr(*),fskyi2(i2size,*),
71 . dmint2(4,*),fncontp(3,*) ,ftcontp(3,*)
72 TYPE (H3D_DATABASE) :: H3D_DATA
73
74
75
76#include "com01_c.inc"
77#include "com06_c.inc"
78#include "com08_c.inc"
79#include "scr11_c.inc"
80#include "scr14_c.inc"
81#include "sms_c.inc"
82#include "task_c.inc"
83
84
85
86 INTEGER NIR,I,J,II,JJ,L,W,NN,KK,K,LLT,
87 . IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
88 . NSVG(MVSIZ)
89
91 . s,t,sp,sm,tp,tm,econtt,econvt,e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,
92 . fnorm,flx,fly,flz,fs(3),ddx,ddy,ddz,xsm,ysm,zsm,xm,ym,zm,
93 . x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,x0,y0,z0,xs,ys,zs,stifms,
94 . vx1,vx2,vx3,vx4,vy1,vy2,vy3,vy4,vz1,vz2,vz3,vz4,dlx,dly,dlz,
95 . vx0,vy0,vz0,vsx,vsy,vsz,vmx,vmy,vmz,v1,v2,v3,dtinv,stf,
96 . fxv,fyv,fzv,drx,dry,drz,stbrk,dti,mharm,dkm, dxt
98 . h(4,mvsiz),fn(3),ft(3),fx(4),fy(4),fz(4),fmx(4),fmy(4),fmz(4),
99 . rx(4),ry(4),rz(4),rm(3),rs(3),v0(3),vs(3),vm(3),dxold(3),
100 . stif(mvsiz), vis(mvsiz), stifm(mvsiz),va(3),vb(3),vc(3),vd(3),
101 . hl(4)
103 . viscl
104
105
106 i7kglo = 1
107 econtt = zero
108 econvt = zero
109
110 DO kk=1,nsn,mvsiz
111
112 llt=
min(nsn-kk+1,mvsiz)
113 DO k=1,llt
114
115 ii= kk+k-1
116 i = nsv(ii)
117
118 IF (i > 0) THEN
119 nsvg(k) = i
120 w = weight(i)
121 s = crst(1,ii)
122 t = crst(2,ii)
123 l = irtl(ii)
124
125 ix1(k) = irect(1,l)
126 ix2(k) = irect(2,l)
127 ix3(k) = irect(3,l)
128 ix4(k) = irect(4,l)
129
130 nir= 4
131 sp = one + s
132 sm = one - s
133 tp = fourth*(one + t)
134 tm = fourth*(one - t)
135 h(1,k)=tm*sm
136 h(2,k)=tm*sp
137 h(3,k)=tp*sp
138 h(4,k)=tp*sm
139
140
141
142
143
144
145
146 IF (ix3(k) == ix4(k)) THEN
147 nir = 3
148 h(3,k) = h(3,k) + h(4,k)
149 h(4,k) = zero
150 ENDIF
151
152
153
154 x1 = x(1,ix1(k))
155 y1 = x(2,ix1(k))
156 z1 = x(3,ix1(k))
157 x2 = x(1,ix2(k))
158 y2 = x(2,ix2(k))
159 z2 = x(3,ix2(k))
160 x3 = x(1,ix3(k))
161 y3 = x(2,ix3(k))
162 z3 = x(3,ix3(k))
163 x4 = x(1,ix4(k))
164 y4 = x(2,ix4(k))
165 z4 = x(3,ix4(k))
166 xs = x(1,i)
167 ys = x(2,i)
168 zs = x(3,i)
169 vsx = v(1,i)
170 vsy = v(2,i)
171 vsz = v(3,i)
172 vx1 = v(1,ix1(k))
173 vy1 = v(2,ix1(k))
174 vz1 = v(3,ix1(k))
175 vx2 = v(1,ix2(k))
176 vy2 = v(2,ix2(k))
177 vz2 = v(3,ix2(k))
178 vx3 = v(1,ix3(k))
179 vy3 = v(2,ix3(k))
180 vz3 = v(3,ix3(k))
181 vx4 = v(1,ix4(k))
182 vy4 = v(2,ix4(k))
183 vz4 = v(3,ix4(k))
184
185 CALL i2rep(x1 ,x2 ,x3 ,x4 ,
186 . y1 ,y2 ,y3 ,y4 ,
187 . z1 ,z2 ,z3 ,z4 ,
188 . e1x ,e1y ,e1z ,
189 . e2x ,e2y ,e2z ,
190 . e3x ,e3y ,e3z ,nir )
191
192 IF (nir == 4) THEN
193 xm = x1*h(1,k) + x2*h(2,k) + x3*h(3,k) + x4*h(4,k)
194 ym = y1*h(1,k) + y2*h(2,k) + y3*h(3,k) + y4*h(4,k)
195 zm = z1*h(1,k) + z2*h(2,k) + z3*h(3,k) + z4*h(4,k)
196 x0 = (x1 + x2 + x3 + x4)/nir
197 y0 = (y1 + y2 + y3 + y4)/nir
198 z0 = (z1 + z2 + z3 + z4)/nir
199
200 xm = xm - x0
201 ym = ym - y0
202 zm = zm - z0
203 xs = xs - x0
204 ys = ys - y0
205 zs = zs - z0
206 xsm = xs - xm
207 ysm = ys - ym
208 zsm = zs - zm
209
210 vx0 = (vx1 + vx2 + vx3 + vx4)/nir
211 vy0 = (vy1 + vy2 + vy3 + vy4)/nir
212 vz0 = (vz1 + vz2 + vz3 + vz4)/nir
213 vmx = vx1*h(1,k) + vx2*h(2,k) + vx3*h(3,k) + vx4*h(4,k) - vx0
214 vmy = vy1*h(1,k) + vy2*h(2,k) + vy3*h(3,k) + vy4*h(4,k) - vy0
215 vmz = vz1*h(1,k) + vz2*h(2,k) + vz3*h(3,k) + vz4*h(4,k) - vz0
216 ELSE
217 x0 = (x1 + x2 + x3)/nir
218 y0 = (y1 + y2 + y3)/nir
219 z0 = (z1 + z2 + z3)/nir
220
221 xm = x1*h(1,k) + x2*h(2,k) + x3*h(3,k)
222 ym = y1*h(1,k) + y2*h(2,k) + y3*h(3,k)
223 zm = z1*h(1,k) + z2*h(2,k) + z3*h(3,k)
224
225 xm = xm - x0
226 ym = ym - y0
227 zm = zm - z0
228 xs = xs - x0
229 ys = ys - y0
230 zs = zs - z0
231 xsm = xs - xm
232 ysm = ys - ym
233 zsm = zs - zm
234
235 vx0 = (vx1 + vx2 + vx3)/nir
236 vy0 = (vy1 + vy2 + vy3)/nir
237 vz0 = (vz1 + vz2 + vz3)/nir
238 vmx = vx1*h(1,k) + vx2*h(2,k) + vx3*h(3,k) - vx0
239 vmy = vy1*h(1,k) + vy2*h(2,k) + vy3*h(3,k) - vy0
240 vmz = vz1*h(1,k) + vz2*h(2,k) + vz3*h(3,k) - vz0
241 ENDIF
242 x1 = x1 - x0
243 y1 = y1 - y0
244 z1 = z1 - z0
245 x2 = x2 - x0
246 y2 = y2 - y0
247 z2 = z2 - z0
248 x3 = x3 - x0
249 y3 = y3 - y0
250 z3 = z3 - z0
251 x4 = x4 - x0
252 y4 = y4 - y0
253 z4 = z4 - z0
254 vsx = vsx - vx0
255 vsy = vsy - vy0
256 vsz = vsz - vz0
257
258
259
260 rs(1) = xs*e1x + ys*e1y + zs*e1z
261 rs(2) = xs*e2x + ys*e2y + zs*e2z
262 rs(3) = xs*e3x + ys*e3y + zs*e3z
263 rm(1) = xm*e1x + ym*e1y + zm*e1z
264 rm(2) = xm*e2x + ym*e2y + zm*e2z
265 rm(3) = xm*e3x + ym*e3y + zm*e3z
266
267 rx(1) = e1x*x1 + e1y*y1 + e1z*z1
268 ry(1) = e2x*x1 + e2y*y1 + e2z*z1
269 rz(1) = e3x*x1 + e3y*y1 + e3z*z1
270 rx(2) = e1x*x2 + e1y*y2 + e1z*z2
271 ry(2) = e2x*x2 + e2y*y2 + e2z*z2
272 rz(2) = e3x*x2 + e3y*y2 + e3z*z2
273 rx(3) = e1x*x3 + e1y*y3 + e1z*z3
274 ry(3) = e2x*x3 + e2y*y3 + e2z*z3
275 rz(3) = e3x*x3 + e3y*y3 + e3z*z3
276 rx(4) = e1x*x4 + e1y*y4 + e1z*z4
277 ry(4) = e2x*x4 + e2y*y4 + e2z*z4
278 rz(4) = e3x*x4 + e3y*y4 + e3z*z4
279
280 IF(nir==3)THEN
281 rx(4)=zero
282 ry(4)=zero
283 rz(4)=zero
284 END IF
285
286 vs(1) = vsx*e1x + vsy*e1y + vsz*e1z
287 vs(2) = vsx*e2x + vsy*e2y + vsz*e2z
288 vs(3) = vsx*e3x + vsy*e3y + vsz*e3z
289 vm(1) = vmx*e1x + vmy*e1y + vmz*e1z
290 vm(2) = vmx*e2x + vmy*e2y + vmz*e2z
291 vm(3) = vmx*e3x + vmy*e3y + vmz*e3z
292
293 va(1) = vx1*e1x + vy1*e1y + vz1*e1z
294 va(2) = vx1*e2x + vy1*e2y + vz1*e2z
295 va(3) = vx1*e3x + vy1*e3y + vz1*e3z
296
297 vb(1) = vx2*e1x + vy2*e1y + vz2*e1z
298 vb(2) = vx2*e2x + vy2*e2y + vz2*e2z
299 vb(3) = vx2*e3x + vy2*e3y + vz2*e3z
300
301 vc(1) = vx3*e1x + vy3*e1y + vz3*e1z
302 vc(2) = vx3*e2x + vy3*e2y + vz3*e2z
303 vc(3) = vx3*e3x + vy3*e3y + vz3*e3z
304
305 vd(1) = vx4*e1x + vy4*e1y + vz4*e1z
306 vd(2) = vx4*e2x + vy4*e2y + vz4*e2z
307 vd(3) = vx4*e3x + vy4*e3y + vz4*e3z
308
309 v1 = vs(1) - vm(1)
310 v2 = vs(2) - vm(2)
311 v3 = vs(3) - vm(3)
312
313
314 IF (tt == zero) THEN
315 dx(1,ii) = zero
316 dx(2,ii) = zero
317 dx(3,ii) = zero
318 fini(1,ii) = zero
319 fini(2,ii) = zero
320 fini(3,ii) = zero
321 ENDIF
322
323 CALL i2pen_rot(skew(1,ii) ,tt ,dt1 ,stbrk,
324 . rs ,rm ,v1 ,v2 ,v3 ,
325 . rx ,ry ,rz ,va ,vb ,
326 . vc ,vd)
327
328 dlx = v1*dt1
329 dly = v2*dt1
330 dlz = v3*dt1
331
332 dx(1,ii) = dx(1,ii) + dlx
333 dx(2,ii) = dx(2,ii) + dly
334 dx(3,ii) = dx(3,ii) + dlz
335
336
337
338
339 flx = dx(1,ii) * stfn(ii)
340 fly = dx(2,ii) * stfn(ii)
341 flz = dx(3,ii) * stfn(ii)
342 viscl = visc
343
344 IF (ivisc==1) THEN
345
346 mharm = ms(i)
347 ELSEIF(ms(i)==zero.OR.ms(ix1(k))==zero.OR.
348 . ms(ix2(k))==zero.OR.
349 . ms(ix3(k))==zero.OR.
350 . ms(ix4(k))==zero)THEN
351 mharm = zero
352 viscl = zero
353 ELSEIF(nir==4)THEN
354 mharm = one/ms(i) +
355 . one/ms(ix1(k)) + one/ms(ix2(k)) + one/ms(ix3(k)) + one/ms(ix4(k))
356 mharm = one/mharm
357 ELSE
358 mharm = one/ms(i) +
359 . one/ms(ix1(k)) + one/ms(ix2(k)) + one/ms(ix3(k))
360 mharm = one/mharm
361 END IF
362 dkm = two*stfn(ii)*mharm
363 vis(k) = visc*sqrt(dkm)
364
365 fxv = vis(k) * v1
366 fyv = vis(k) * v2
367 fzv = vis(k) * v3
368
369 dxt = dx(1,ii)**2 + dx(2,ii)**2 + dx(3,ii)**2
370 econtt = econtt + half*stfn(ii)*dxt*w
371
372 econvt = econvt + (fxv*v1 + fyv*v2 + fzv*v3)*dt1*w
373
374 flx = flx + fxv
375 fly = fly + fyv
376 flz = flz + fzv
377
378 DO j=1,4
379 fmx(j) = h(j,k)*flx
380 fmy(j) = h(j,k)*fly
381 fmz(j) = h(j,k)*flz
382 ENDDO
383
384
385 CALL i2loceq( nir ,rs ,rx ,ry ,rz ,
386 . fmx ,fmy ,fmz ,h(1,k) ,stifm(k))
387
388
389
390 DO j=1,4
391 fx(j) = e1x*fmx(j) + e2x*fmy(j) + e3x*fmz(j)
392 fy(j) = e1y*fmx(j) + e2y*fmy(j) + e3y*fmz(j)
393 fz(j) = e1z*fmx(j) + e2z*fmy(j) + e3z*fmz(j)
394 ENDDO
395 fs(1) = zero
396 fs(2) = zero
397 fs(3) = zero
398 DO j=1,nir
399 fs(1) = fs(1) + fx(j)
400 fs(2) = fs(2) + fy(j)
401 fs(3) = fs(3) + fz(j)
402 ENDDO
403 a(1,i) = a(1,i) - fs(1)
404 a(2,i) = a(2,i) - fs(2)
405 a(3,i) = a(3,i) - fs(3)
406
407 IF (ivisc==1) THEN
408
409 dtinv = zero
410 IF (dt1 > zero) dtinv=one/dt1
411 stf = (one+stbrk)*stfn(ii) + two*vis(k)*dtinv
412 ELSE
413 stf = stfn(ii)*(viscl + sqrt(viscl**2 + (one+stbrk)))**2
414 ENDIF
415
416 stifn(i) = stifn(i) + stf
417
418
419 stif(k) = (one+stbrk)*stfn(ii)
420
421
422
423 IF (w == 1) THEN
424 i0 = i0 + 1
425 jj = 1
426 nn = iadi2(jj,i0)
427
428 fskyi2(1,nn) = fx(jj)
429 fskyi2(2,nn) = fy(jj)
430 fskyi2(3,nn) = fz(jj)
431 fskyi2(4,nn) = zero
432 fskyi2(5,nn) = abs(stf*h(jj,k))+stifm(k)*stf
433 IF (iroddl == 1) THEN
434 fskyi2(6,nn) = zero
435 fskyi2(7,nn) = zero
436 fskyi2(8,nn) = zero
437 fskyi2(9,nn) = zero
438 fskyi2(10,nn)= zero
439 ENDIF
440
441 jj = 2
442 nn = iadi2
443 fskyi2(1,nn) = fx(jj)
444 fskyi2(2,nn) = fy(jj)
445 fskyi2(3,nn) = fz(jj)
446 fskyi2(4,nn) = zero
447 fskyi2(5,nn) = abs(stf*h(jj,k))+stifm(k)*stf
448 IF (iroddl == 1) THEN
449 fskyi2(6,nn) = zero
450 fskyi2(7,nn) = zero
451 fskyi2(8,nn) = zero
452 fskyi2(9,nn) = zero
453 fskyi2(10,nn)= zero
454 ENDIF
455
456 jj = 3
457 nn = iadi2(jj,i0)
458 fskyi2(1,nn) = fx(jj)
459 fskyi2(2,nn) = fy(jj)
460 fskyi2(3,nn) = fz(jj)
461 fskyi2(4,nn) = zero
462 fskyi2(5,nn) = abs(stf*h(jj,k))+stifm(k)*stf
463 IF (iroddl == 1) THEN
464 fskyi2(6,nn) = zero
465 fskyi2(7,nn) = zero
466 fskyi2(8,nn) = zero
467 fskyi2(9,nn) = zero
468 fskyi2(10,nn)= zero
469 ENDIF
470
471 jj = 4
472 nn = iadi2(jj,i0)
473 fskyi2(1,nn) = fx(jj)
474 fskyi2(2,nn) = fy(jj)
475 fskyi2(3,nn) = fz(jj)
476 fskyi2(4,nn) = zero
477 fskyi2(5,nn) = abs(stf*h(jj,k))
478 . +stifm(k)*stf*sign(one,abs(h(jj,k)))
479 IF (iroddl == 1) THEN
480 fskyi2(6,nn) = zero
481 fskyi2(7,nn) = zero
482 fskyi2(8,nn) = zero
483 fskyi2(9,nn) = zero
484 fskyi2(10,nn)= zero
485 ENDIF
486 ENDIF
487
488 fini(1,ii) = flx
489 fini(2,ii) = fly
490 fini(3,ii) = flz
491
492
493
494 hl(1:4) = h(1:4,k)
496 . irect(1,l),nir ,fsav ,fncont ,fncontp,
497 . ftcontp ,weight ,h3d_data,i ,hl)
498
499
500 ELSE
501 nsvg(k)= -i
502 l = irtl(ii)
503
504 ix1(k) = irect(1,l)
505 ix2(k) = irect(2,l)
506 ix3(k) = irect(3,l)
507 ix4(k) = irect(4,l)
508 stif(k)= zero
509 vis(k) = zero
510 IF (weight(-i) == 1) THEN
511 i0 = i0 + 1
512 jj = 1
513 nn = iadi2(jj,i0)
514 fskyi2(1,nn) = zero
515 fskyi2(2,nn) = zero
516 fskyi2(3,nn) = zero
517 fskyi2(4,nn) = zero
518 fskyi2(5,nn) = zero
519 IF (iroddl == 1) THEN
520 fskyi2(6,nn) = zero
521 fskyi2(7,nn) = zero
522 fskyi2(8,nn) = zero
523 fskyi2(9,nn) = zero
524 fskyi2(10,nn)= zero
525 ENDIF
526 jj = 2
527 nn = iadi2(jj,i0)
528 fskyi2(1,nn) = zero
529 fskyi2(2,nn) = zero
530 fskyi2(3,nn) = zero
531 fskyi2(4,nn) = zero
532 fskyi2(5,nn) = zero
533 IF (iroddl == 1) THEN
534 fskyi2(6,nn) = zero
535 fskyi2(7,nn) = zero
536 fskyi2(8,nn) = zero
537 fskyi2(9,nn) = zero
538 fskyi2(10,nn)= zero
539 ENDIF
540 jj = 3
541 nn = iadi2(jj,i0)
542 fskyi2(1,nn) = zero
543 fskyi2(2,nn) = zero
544 fskyi2(3,nn) = zero
545 fskyi2(4,nn) = zero
546 fskyi2(5,nn) = zero
547 IF (iroddl == 1) THEN
548 fskyi2(6,nn) = zero
549 fskyi2(7,nn) = zero
550 fskyi2(8,nn) = zero
551 fskyi2(9,nn) = zero
552 fskyi2(10,nn)= zero
553 ENDIF
554 jj = 4
555 nn = iadi2(jj,i0)
556 fskyi2(1,nn) = zero
557 fskyi2(2,nn) = zero
558 fskyi2(3,nn) = zero
559 fskyi2(4,nn) = zero
560 fskyi2(5,nn) = zero
561 IF (iroddl == 1) THEN
562 fskyi2(6,nn) = zero
563 fskyi2(7,nn) = zero
564 fskyi2(8,nn) = zero
565 fskyi2(9,nn) = zero
566 fskyi2(10,nn)= zero
567 ENDIF
568 ENDIF
569 ENDIF
570 ENDDO
571
572 IF(idtmins==2.OR.idtmins_int/=0)THEN
573 dti=dt2t
574 CALL i2sms25(llt ,ix1 ,ix2 ,ix3 ,ix4 ,
575 2 nsvg ,h ,stif ,noint ,
576 3 dmint2(1,kk),nodnx_sms ,vis ,
577 4 stifm ,dti )
578 IF(dti<dt2t)THEN
579 dt2t = dti
580 neltst = noint
581 ityptst = 10
582 ENDIF
583 END IF
584 ENDDO
585
586#include "lockon.inc"
587 econt = econt + econtt
588 econtd = econtd + econvt
589 fsav(26) = fsav(26) + econtt
590 fsav(28) = fsav(28) + econvt
591#include "lockoff.inc"
592
593 RETURN
subroutine i2forces(x, fs, fx, fy, fz, irect, nir, fsav, fncont, fncontp, ftcontp, weight, h3d_data, nsl, h)
subroutine i2loceq(nir, rs, rx, ry, rz, fmx, fmy, fmz, h, stifm)
subroutine i2pen_rot(skew, tt, dt1, stif, rs, rm, v1, v2, v3, rx, ry, rz, va, vb, vc, vd)
subroutine i2rep(x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, nir)
subroutine i2sms25(jlt, ix1, ix2, ix3, ix4, nsvg, h, stif, noint, dmint2, nodnx_sms, vis, stifm, dti)