40
41
42
44
45
46
47#include "implicit_f.inc"
48#include "comlock.inc"
49
50
51
52#include "mvsiz_p.inc"
53
54
55
56#include "com04_c.inc"
57#include "com08_c.inc"
58#include "units_c.inc"
59
60
61
62 INTEGER NDEB, NTC
63 INTEGER KSURF,
64 . KSI(4,*),IACTIV(4,*),NPC(*),KTC(*),
65 . NOINT,NELTST,ITYPTST
66
67 my_real stfac, bufsf(*),x(3,*) , iold(3,*),
68 . hold(3,*) , nold(3,*) ,dold(3,*),ms(*) , v(3,*),
69 . de, pld(*), xtk(*) ,ytk(*) ,ztk(*) ,
70 . ntx(*) ,nty(*) ,ntz(*) ,
71 . penet(*) ,depth(*),
72 . xi(*) ,yi(*) ,zi(*) ,nxi(*) ,nyi(*) ,nzi(*) ,
73 . wnf(3,*) ,wtf(3,*) ,wns(*) ,xp1(3,*) ,xp2(3,*) ,xp3(3,*) ,
74 . fnormx,fnormy,fnormz,ftangx,ftangy,ftangz,
75 . dt2t,vfric
76 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
77
78
79
80 INTEGER ADRBUF, I, IL, IN1, IN2, IN3, IN4, NLS
82 . a, b, c, rot(9), x1, x2, x3, xm, ym, zm,
83 . v1, v2, v3, vxm, vym, vzm, vrx, vry, vrz,
84 . v1x2, v2x1, v1x3, v3x1, v2x3, v3x2,
85 . deltx, delty, deltz, nd, pn, scale, nv,
86 . vxk, vyk, vzk, dvx, dvy, dvz, xh, yh, zh,
87 . dx1, dy1, dz1, dx2, dy2, dz2, dx3, dy3, dz3,
88 . dx, dy, dz,
89 . s1, s2, s3, s,
90 . f1, f11, f12, f13,
91 . ftx, fty, ftz, fnx, fny, fnz, ftn, fmax,fn,
92 . h, dti
94 . fxn(mvsiz),fyn(mvsiz),fzn(mvsiz),stf(mvsiz),
95 . fxt(mvsiz),fyt(mvsiz),fzt(mvsiz),
96 . l1(mvsiz),l2(mvsiz),l3(mvsiz),
97 . vx1(mvsiz),vy1(mvsiz),vz1(mvsiz),
98 . vx2(mvsiz),vy2(mvsiz),vz2(mvsiz),
99 . vx3(mvsiz),vy3(mvsiz),vz3(mvsiz),
100 . vxh(mvsiz),vyh(mvsiz),vzh(mvsiz)
101
102 adrbuf=igrsurf(ksurf)%IAD_BUFR
103 a =bufsf(adrbuf+1)
104 b =bufsf(adrbuf+2)
105 c =bufsf(adrbuf+3)
106 DO i=1,9
107 rot(i)=bufsf(adrbuf+7+i-1)
108 END DO
109
110
111
112 x1=bufsf(adrbuf+16)-bufsf(adrbuf+4)
113 x2=bufsf(adrbuf+17)-bufsf(adrbuf+5)
114 x3=bufsf(adrbuf+18)-bufsf(adrbuf+6)
115 xm=rot(1)*x1+rot(2)*x2+rot(3)*x3
116 ym=rot(4)*x1+rot(5)*x2+rot(6)*x3
117 zm=rot(7)*x1+rot(8)*x2+rot(9)*x3
118
119
120
121 v1=bufsf(adrbuf+19)
122 v2=bufsf(adrbuf+20)
123 v3=bufsf(adrbuf+21)
124 vxm=rot(1)*v1+rot(2)*v2+rot(3)*v3
125 vym=rot(4)*v1+rot(5)*v2+rot(6)*v3
126 vzm=rot(7)*v1+rot(8)*v2+rot(9)*v3
127
128
129
130 v1=bufsf(adrbuf+22)
131 v2=bufsf(adrbuf+23)
132 v3=bufsf(adrbuf+24)
133 vrx=rot(1)*v1+rot(2)*v2+rot(3)*v3
134 vry=rot(4)*v1+rot(5)*v2+rot(6)*v3
135 vrz=rot(7)*v1+rot(8)*v2+rot(9)*v3
136
137
138
139
140 DO 100 nls=ndeb+1,
min(ndeb+mvsiz,ntc)
141 il =ktc(nls)
142 IF (iactiv(1,il)<=0) GOTO 100
143 i =nls-ndeb
144
145 in1=ksi(1,il)
146 v1=v(1,in1)
147 v2=v(2,in1)
148 v3=v(3,in1)
149 vx1(i)=rot(1)*v1+rot(2)*v2+rot(3)*v3
150 vy1(i)=rot(4)*v1+rot(5)*v2+rot(6)*v3
151 vz1(i)=rot(7)*v1+rot(8)*v2+rot(9)*v3
152
153 in2=ksi(2,il)
154 v1=v(1,in2)
155 v2=v(2,in2)
156 v3=v(3,in2)
157 vx2(i)=rot(1)*v1+rot(2)*v2+rot(3)*v3
158 vy2(i)=rot(4)*v1+rot(5)*v2+rot(6)*v3
159 vz2(i)=rot(7)*v1+rot(8)*v2+rot(9)*v3
160
161 in3=ksi(3,il)
162 v1=v(1,in3)
163 v2=v(2,in3)
164 v3=v(3,in3)
165 vx3(i)=rot(1)*v1+rot(2)*v2+rot(3)*v3
166 vy3(i)=rot(4)*v1+rot(5)*v2+rot(6)*v3
167 vz3(i)=rot(7)*v1+rot(8)*v2+rot(9)*v3
168 100 CONTINUE
169
170
171
172 DO 200 nls=ndeb+1,
min(ndeb+mvsiz,ntc)
173 il =ktc(nls)
174 IF (iactiv(1,il)<=0) GOTO 200
175 i =nls-ndeb
176
177 stf(i)=stfac*depth(i)**2/
max((depth(i)-penet(i))**2,em20)
178 fxn(i)=stf(i)*penet(i)*nxi(i)
179 fyn(i)=stf(i)*penet(i)*nyi(i)
180 fzn(i)=stf(i)*penet(i)*nzi(i)
181 200 CONTINUE
182
183
184
185 DO 250 nls=ndeb+1,
min(ndeb+mvsiz,ntc)
186 il =ktc(nls)
187 IF (iactiv(1,il)<=0) GOTO 250
188 i =nls-ndeb
189
190 l1(i)=iold(1,4*(il-1)+1)
191 l2(i)=iold(2,4*(il-1)+1)
192 l3(i)=iold(3,4*(il-1)+1)
193 250 CONTINUE
194
195
196
197#include "vectorize.inc"
198 DO 275 nls=ndeb+1,
min(ndeb+mvsiz,ntc)
199 il =ktc(nls)
200 IF (iactiv(1,il)<=0) GOTO 275
201 i =nls-ndeb
202
203 fxt(i)=zero
204 fyt(i)=zero
205 fzt(i)=zero
206
207
208
209 IF (iactiv(1,il)>2) THEN
210 deltx=dold(1,4*(il-1)+1)
211 delty=dold(2,4*(il-1)+1)
212 deltz=dold(3,4*(il-1)+1)
213 nd =sqrt(deltx*deltx+delty*delty+deltz*deltz)
214 IF (nd/=zero) THEN
215 pn=deltx*nxi(i)+delty*nyi(i)+deltz*nzi(i)
216 deltx=deltx-pn*nxi(i)
217 delty=delty-pn*nyi(i)
218 deltz=deltz-pn*nzi(i)
219 scale=nd/sqrt(deltx*deltx+delty*delty+deltz*deltz)
220 deltx=scale*deltx
221 delty=scale*delty
222 deltz=scale*deltz
223 ENDIF
224 ELSE
225 deltx=zero
226 delty=zero
227 deltz=zero
228 ENDIF
229
230
231
232 xh=hold(1,4*(il-1)+1)
233 yh=hold(2,4*(il-1)+1)
234 zh=hold(3,4*(il-1)+1)
235 v1x2=vrx*(yh-ym)
236 v2x1=vry*(xh-xm)
237 v2x3=vry*(zh-zm)
238 v3x2=vrz*(yh-ym)
239 v3x1=vrz*(xh-xm)
240 v1x3=vrx*(zh-zm)
241 v3 =v1x2-v2x1
242 v1 =v2x3-v3x2
243 v2 =v3x1-v1x3
244 vxh(i)=vxm+v1
245 vyh(i)=vym+v2
246 vzh(i)=vzm+v3
247 IF (iactiv(1,il)>=2) THEN
248 vxk=l1(i)*vx1(i)+l2(i)*vx2(i)+l3(i)*vx3(i)
249 vyk=l1(i)*vy1(i)+l2(i)*vy2(i)+l3(i)*vy3(i)
250 vzk=l1(i)*vz1(i)+l2(i)*vz2(i)+l3(i)*vz3(i)
251 dvx=vxh(i)-vxk
252 dvy=vyh(i)-vyk
253 dvz=vzh(i)-vzk
254 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
255 dvx=dvx-pn*nxi(i)
256 dvy=dvy-pn*nyi(i)
257 dvz=dvz-pn*nzi(i)
258 ELSE
259 dvx=zero
260 dvy=zero
261 dvz=zero
262 ENDIF
263
264
265 dold(1,4*(il-1)+1)=deltx+dvx*dt1
266 dold(2,4*(il-1)+1)=delty+dvy*dt1
267 dold(3,4*(il-1)+1)=deltz+dvz*dt1
268 fxt(i)=stf(i)*dold(1,4*(il-1)+1)
269 fyt(i)=stf(i)*dold(2,4*(il-1)+1)
270 fzt(i)=stf(i)*dold(3,4*(il-1)+1)
271
272 275 CONTINUE
273
274
275
276#include "vectorize.inc"
277 DO 625 nls=ndeb+1,
min(ndeb+mvsiz,ntc)
278 il =ktc(nls)
279 IF (iactiv(1,il)<=0) GOTO 625
280 i =nls-ndeb
281
282 ftx=fxt(i)
283 fty=fyt(i)
284 ftz=fzt(i)
285 fnx=fxn(i)
286 fny=fyn(i)
287 fnz=fzn(i)
288 ftn=sqrt(ftx*ftx+fty*fty+ftz*ftz)
289 fn =sqrt(fnx*fnx+fny*fny+fnz*fnz)
290 fmax=
max(vfric*fn,zero)
291 IF (ftn>fmax) THEN
292 scale =fmax/ftn
293 ftx=scale*ftx
294 fty=scale*fty
295 ftz=scale*ftz
296 ELSE
297 scale=1.
298 ENDIF
299 IF (iactiv(1,il)>1) THEN
300
301 deltx=scale*dold(1,4*(il-1)+1)
302 delty=scale*dold(2,4*(il-1)+1)
303 deltz=scale*dold(3,4*(il-1)+1)
304
305 dold(1,4*(il-1)+1)=deltx
306 dold(2,4*(il-1)+1)=delty
307 dold(3,4*(il-1)+1)=deltz
308 ENDIF
309 fxt(i)=ftx
310 fyt(i)=fty
311 fzt(i)=ftz
312 625 CONTINUE
313
314
315
316
317 DO 650 nls=ndeb+1,
min(ndeb+mvsiz,ntc)
318 il =ktc(nls)
319 IF (iactiv(1,il)<=0) GOTO 650
320 i =nls-ndeb
321
322 dx2=xp2(1,i)-xtk(i)
323 dy2=xp2(2,i)-ytk(i)
324 dz2=xp2(3,i)-ztk(i)
325 dx3=xp3(1,i)-xtk(i)
326 dy3=xp3(2,i)-ytk(i)
327 dz3=xp3(3,i)-ztk(i)
328 dx=dy2*dz3-dz2*dy3
329 dy=dx3*dz2-dz3*dx2
330 dz=dx2*dy3-dy2*dx3
331 s1=half*sqrt(dx*dx+dy*dy+dz*dz)
332
333 dx1=xp1(1,i) -xtk(i)
334 dy1=xp1(2,i) -ytk(i)
335 dz1=xp1(3,i) -ztk(i)
336 dx=dy1*dz3-dz1*dy3
337 dy=dx3*dz1-dz3*dx1
338 dz=dx1*dy3-dy1*dx3
339 s2=half*sqrt(dx*dx+dy*dy+dz*dz)
340
341 dx=dy1*dz2-dz1*dy2
342 dy=dx2*dz1-dz2*dx1
343 dz=dx1*dy2-dy1*dx2
344 s3=half*sqrt(dx*dx+dy*dy+dz*dz)
345
346 s=one/(s1+s2+s3)
347
348 iold(1,4*(il-1)+1)=s1*s
349 iold(2,4*(il-1)+1)=s2*s
350 iold(3,4*(il-1)+1)=s3*s
351 650 CONTINUE
352
353
354
355#include "vectorize.inc"
356 DO 700 nls=ndeb+1,
min(ndeb+mvsiz,ntc)
357 il =ktc(nls)
358 IF (iactiv(1,il)<=0) GOTO 700
359 i =nls-ndeb
360
361 hold(1,4*(il-1)+1)=xi(i)
362 hold(2,4*(il-1)+1)=yi(i)
363 hold(3,4*(il-1)+1)=zi(i)
364
365 nold(1,4*(il-1)+1)=nxi(i)
366 nold(2,4*(il-1)+1)=nyi(i)
367 nold(3,4*(il-1)+1)=nzi(i)
368 700 CONTINUE
369
370
371
372
373 DO 600 nls=ndeb+1,
min(ndeb+mvsiz,ntc)
374 il =ktc(nls)
375 IF (iactiv(1,il)<=0) GOTO 600
376 i =nls-ndeb
377
378 in1=ksi(1,il)
379 in2=ksi(2,il)
380 in3=ksi(3,il)
381
382 wns(in1) =wns(in1)+iold(1,4*(il-1)+1)*stf(i)
383 wns(in2) =wns(in2)+iold(2,4*(il-1)+1)*stf(i)
384 wns(in3) =wns(in3)+iold(3,4*(il-1)+1)*stf(i)
385
386 600 CONTINUE
387
388
389 DO 800 nls=ndeb+1,
min(ndeb+mvsiz,ntc)
390 il =ktc(nls)
391 IF (iactiv(1,il)<=0) GOTO 800
392 i =nls-ndeb
393
394 in1=ksi(1,il)
395 in2=ksi(2,il)
396 in3=ksi(3,il)
397
398 f1=fxn(i)
399 fnormx=fnormx+f1
400 f11=iold(1,4*(il-1)+1)*f1
401 f12=iold(2,4*(il-1)+1)*f1
402 f13=iold(3,4*(il-1)+1)*f1
403 wnf(1,in1)=wnf(1,in1)+f11
404 wnf(1,in2)=wnf(1,in2)+f12
405 wnf(1,in3)=wnf(1,in3)+f13
406
407 f1=fyn(i)
408 fnormy=fnormy+f1
409 f11=iold(1,4*(il-1)+1)*f1
410 f12=iold(2,4*(il-1)+1)*f1
411 f13=iold(3,4*(il-1)+1)*f1
412 wnf(2,in1)=wnf(2,in1)+f11
413 wnf(2,in2)=wnf(2,in2)+f12
414 wnf(2,in3)=wnf(2,in3)+f13
415
416 f1=fzn(i)
417 fnormz=fnormz+f1
418 f11=iold(1,4*(il-1)+1)*f1
419 f12=iold(2,4*(il-1)+1)*f1
420 f13=iold(3,4*(il-1)+1)*f1
421 wnf(3,in1)=wnf(3,in1)+f11
422 wnf(3,in2)=wnf(3,in2)+f12
423 wnf(3,in3)=wnf(3,in3)+f13
424
425 800 CONTINUE
426
427
428 DO 825 nls=ndeb+1,
min(ndeb+mvsiz,ntc)
429 il =ktc(nls)
430 IF (iactiv(1,il)<=0) GOTO 825
431 i =nls-ndeb
432
433 in1=ksi(1,il)
434 in2=ksi(2,il)
435 in3=ksi(3,il)
436
437 f1=fxt(i)
438 ftangx=ftangx+f1
439 f11=iold(1,4*(il-1)+1)*f1
440 f12=iold(2,4*(il-1)+1)*f1
441 f13=iold(3,4*(il-1)+1)*f1
442 wtf(1,in1)=wtf(1,in1)+f11
443 wtf(1,in2)=wtf(1,in2)+f12
444 wtf(1,in3)=wtf(1,in3)+f13
445
446 f1=fyt(i)
447 ftangy=ftangy+f1
448 f11=iold(1,4*(il-1)+1)*f1
449 f12=iold(2,4*(il-1)+1)*f1
450 f13=iold(3,4*(il-1)+1)*f1
451 wtf(2,in1)=wtf(2,in1)+f11
452 wtf(2,in2)=wtf(2,in2)+f12
453 wtf(2,in3)=wtf(2,in3)+f13
454
455 f1=fzt(i)
456 ftangz=ftangz+f1
457 f11=iold(1,4*(il-1)+1)*f1
458 f12=iold(2,4*(il-1)+1)*f1
459 f13=iold(3,4*(il-1)+1)*f1
460 wtf(3,in1)=wtf(3,in1)+f11
461 wtf(3,in2)=wtf(3,in2)+f12
462 wtf(3,in3)=wtf(3,in3)+f13
463
464 825 CONTINUE
465
466
467
468 dti=ep20
469 DO 900 nls=ndeb+1,
min(ndeb+mvsiz,ntc)
470 il =ktc(nls)
471 IF (iactiv(1,il)<=0) GOTO 900
472 i =nls-ndeb
473 dvx=vx1(i)-vxh(i)
474 dvy=vy1(i)-vyh(i)
475 dvz=vz1(i)-vzh(i)
476 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
477 IF (pn<zero) THEN
478 h =xp1(1,i)*nxi(i)+xp1(2,i)*nyi(i)+xp1(3,i)*nzi(i)
479 dti=
min(dti,half*h/
max(em20,-pn))
480 ENDIF
481 dvx=vx2(i)-vxh(i)
482 dvy=vy2(i)-vyh(i)
483 dvz=vz2(i)-vzh(i)
484 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
485 IF (pn<zero) THEN
486 h =xp2(1,i)*nxi(i)+xp2(2,i)*nyi(i)+xp2(3,i)*nzi(i)
487 dti=
min(dti,half*h/
max(em20,-pn))
488 ENDIF
489 dvx=vx3(i)-vxh(i)
490 dvy=vy3(i)-vyh(i)
491 dvz=vz3(i)-vzh(i)
492 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
493 IF (pn<zero) THEN
494 h =xp3(1,i)*nxi(i)+xp3(2,i)*nyi(i)+xp3(3,i)*nzi(i)
495 dti=
min(dti,half*h/
max(em20,-pn))
496 ENDIF
497 IF(dti<dt2t)THEN
498 dt2t = dti
499 neltst = noint
500 ityptst = 10
501 ENDIF
502
503 900 CONTINUE
504
505 IF (dti<=0.) THEN
506 dti=ep20
507 DO 950 nls=ndeb+1,
min(ndeb+mvsiz,ntc)
508 il =ktc(nls)
509 IF (iactiv(1,il)<=0) GOTO 950
510 i =nls-ndeb
511 dvx=vx1(i)-vxh(i)
512 dvy=vy1(i)-vyh(i)
513 dvz=vz1(i)-vzh(i)
514 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
515 IF (pn<zero) THEN
516 h =xp1(1,i)*nxi(i)+xp1(2,i)*nyi(i)+xp1(3,i)*nzi(i)
517 dti=
min(dti,half*h/
max(em20,-pn))
518 ENDIF
519 dvx=vx2(i)-vxh(i)
520 dvy=vy2(i)-vyh(i)
521 dvz=vz2(i)-vzh(i)
522 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
523 IF (pn<zero) THEN
524 h =xp2(1,i)*nxi(i)+xp2(2,i)*nyi(i)+xp2(3,i)*nzi(i)
525 dti=
min(dti,half*h/
max(em20,-pn))
526 ENDIF
527 dvx=vx3(i)-vxh(i)
528 dvy=vy3(i)-vyh(i)
529 dvz=vz3(i)-vzh(i)
530 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
531 IF (pn<zero) THEN
532 h =xp3(1,i)*nxi(i)+xp3(2,i)*nyi(i)+xp3(3,i)*nzi(i)
533 dti=
min(dti,half*h/
max(em20,-pn))
534 ENDIF
535 IF(dti<=zero)THEN
536 in1=ksi(1,il)
537 in2=ksi(2,il)
538 in3=ksi(3,il)
539#include "lockon.inc"
540 WRITE(istdo,'(A,E12.4,A,I8)')
541 . ' **WARNING MINIMUM TIME STEP ',dti,' IN INTERFACE ',noint
542 WRITE(iout ,'(A,E12.4,A,I8)')
543 . ' **WARNING MINIMUM TIME STEP ',dti,' IN INTERFACE ',noint
544 WRITE(iout,'(A,3I8)') ' ELEMENT/SEGMENT NODES :',
545 . in1,in2,in3
546#include "lockoff.inc"
547 tstop = tt
548 ENDIF
549
550 950 CONTINUE
551 ENDIF
552
553 RETURN