40 1 JLT ,A ,V ,IBC ,ICODT ,
41 2 FSAV ,GAP ,FRIC ,MS ,VISC ,
42 3 VISCF ,NOINT ,ITAB ,CS_LOC ,CM_LOC ,
43 4 STIGLO ,STIFN ,STIF ,FSKYI ,ISKY ,
44 5 FCONT ,DT2T ,NRTM ,MSEGTYP ,HS1 ,
45 6 HS2 ,HM1 ,HM2 ,N1 ,N2 ,
46 7 M1 ,M2 ,IVIS2 ,NELTST ,ITYPTST,
47 8 NX ,NY ,NZ ,GAPVE ,INACTI ,
48 9 INDEX ,CAND_P ,NISKYFIE,NEWFRONT,ISECIN ,
49 A NSTRF ,SECFCUM,VISCN ,NEDGE ,MS1 ,
50 B MS2 ,MM1 ,MM2 ,VXS1 ,VYS1 ,
51 C VZS1 ,VXS2 ,VYS2 ,VZS2 ,VXM1 ,
52 D VYM1 ,VZM1 ,VXM2 ,VYM2 ,VZM2 ,
53 E NIN ,NISUB ,LISUB ,ADDSUBE ,ADDSUBM,
54 F LISUBE ,LISUBM ,INFLG_SUBE,INFLG_SUBM,FSAVSUB,
55 G MSKYI_SMS,ISKYI_SMS,NSMS ,JTASK ,ISENSINT,
56 H FSAVPARIT,NFT ,H3D_DATA,INDX1 ,INDX2 ,
57 I ILEV ,MBINFLG, EDGE_ID, NEDGE_REM ,FRICC ,
58 J IFQ ,CAND_FX, CAND_FY, CAND_FZ ,IFPEN ,
59 K TAGNCONT ,KLOADPINTER,LOADPINTER ,LOADP_HYD_INTER,
60 L TYPSUB ,STARTT,NINLOADP ,DGAPLOADINT,S_LOADPINTER)
70#include "implicit_f.inc"
96 INTEGER :: EDGE_ID(2,4*MVSIZ),NEDGE_REM
97 INTEGER NELTST,ITYPTST,JLT,IBC,IVIS2,INACTI,NEDGE,NISKYFIE,NIN,NRTM,ILEV,
99 INTEGER ICODT(*), ITAB(*), ISKY(*),
100 . NOINT,NEWFRONT,ISECIN, NSTRF(*), ISKYI_SMS(*), MSEGTYP(*),
101 . NISUB, LISUB(*), ADDSUBE(*), ADDSUBM(*), LISUBE(*), LISUBM(*),
102 . INFLG_SUBE(*), INFLG_SUBM(*), MBINFLG(*), IFPEN(*),TYPSUB(*)
103 INTEGER N1(*), N2(*), M1(*), M2(*), NSMS(*),
104 . CS_LOC(4*MVSIZ), CM_LOC(4*MVSIZ), JTASK,
105 . ISENSINT(*),NFT,INDEX(*), INDX1(4*MVSIZ), INDX2(4*MVSIZ),
106 . TAGNCONT(NLOADP_HYD_INTER,NUMNOD)
107 INTEGER ,
INTENT(IN) :: NINLOADP,S_LOADPINTER
108 INTEGER ,
INTENT(IN) :: KLOADPINTER(NINTER+1),LOADPINTER(S_LOADPINTER),
109 . LOADP_HYD_INTER(NLOADP_HYD)
112 . A(3,*), MS(*), V(3,*), FSAV(*),FCONT(3,*),
113 . STIFN(*),FSKYI(LSKYI,NFSKYI),FSAVSUB(NTHVKI,*),
114 . MSKYI_SMS(*), GAPVE(*), CAND_P(4,*),
115 . GAP,FRIC,VISC,VISCF,VIS,DT2T,STARTT
117 . hs1(*), hs2(*), hm1(*), hm2(*),
118 . nx(*), ny(*), nz(*), stif(*),
119 . secfcum(7,numnod,nsect), viscn(*),
120 . ms1(*),ms2(*),mm1(*),mm2(*),
121 . vxs1(*),vys1(*),vzs1(*),vxs2(*),vys2(*),
122 . vzs2(*),vxm1(*),vym1(*),vzm1(*),vxm2(*),
123 . vym2(*),vzm2(*),fsavparit(nisub+1,11,*),
124 . fricc(*),cand_fx(4,*),cand_fy(4,*),cand_fz(4,*)
125 my_real ,
INTENT(IN) :: dgaploadint(s_loadpinter)
130 INTEGER I, J1, J , K0,NBINTER,K1S,K, NI, IL, IE, IG, PP, PPL
131 INTEGER NISKYL,NISKYL1,ISIGN
132 INTEGER JSUB,KSUB,NSUB,JJ,KK,ISS1,ISS2,IMS1,IMS2,ITYPSUB,
135 . VX(4*MVSIZ), VY(4*MVSIZ), VZ(4*MVSIZ), VN(4*MVSIZ),
136 . FXI(4*MVSIZ), FYI(4*MVSIZ), FZI(4*MVSIZ), FNI(4*MVSIZ),
137 . fx1(4*mvsiz), fx2(4*mvsiz), fx3(4*mvsiz), fx4(4*mvsiz),
138 . fy1(4*mvsiz), fy2(4*mvsiz), fy3(4*mvsiz), fy4(4*mvsiz),
139 . fz1(4*mvsiz), fz2(4*mvsiz), fz3(4*mvsiz), fz4(4*mvsiz),
140 . fxt(4*mvsiz), fyt(4*mvsiz), fzt(4*mvsiz),
141 . vis2(4*mvsiz),pene(4*mvsiz),dist(4*mvsiz),
142 . vnx, vny, vnz, aa, vmax,s2,
143 . v2, fm2, dt1inv, visca, fac, ff,
144 . fx, fy, fz, f2, mas2, dtmi0,dti,
145 . facm1, econtt, econvt, a2,masm,
146 . fsav1, fsav2, fsav3, fsav4, fsav5, fsav6, fsav7, fsav8,
147 . fsav9, fsav10, fsav11, fsav12, fsav13, fsav14, fsav15,
148 . fsav22, fsav23, fsav24,dgapload,
149 . fsavsub1(24,nisub), impx, impy, impz,ftn ,fn , ft,beta
153 . st1(4*mvsiz),st2(4*mvsiz),st3(4*mvsiz),st4(4*mvsiz),stif0(4*mvsiz),
154 . kt(4*mvsiz),c(4*mvsiz),cf(4*mvsiz),
155 . k1(4*mvsiz),k2(4*mvsiz),k3(4*mvsiz),k4(4*mvsiz),
156 . c1(4*mvsiz),c2(4*mvsiz),c3(4*mvsiz),c4(4*mvsiz),
159 . fx6(6,4*mvsiz), fy6(6,4*mvsiz), fz6(6,4*mvsiz)
182 s2 = sqrt(nx(i)**2 + ny(i)**2 + nz(i)**2)
184 IF(gapve(i)/=zero)
THEN
185 pene(i) =
max(zero,gapve(i) - s2)
190 s2 = one/
max(em30,s2)
198 debug_e2e(edge_id(1,i)==d_em .AND. edge_id(2,i) == d_es,cand_p(indx2(i),indx1(i)))
199 debug_e2e(edge_id(1,i)==d_em .AND. edge_id(2,i) == d_es,pene(i))
200 IF(cand_p(indx2(i),indx1(i))==zero)cand_p(indx2(i),indx1(i))=pene(i)
206 IF(cand_p(indx2(i),indx1(i))<zero)
THEN
208 cand_p(indx2(i),indx1(i))=pene(i)
210 cand_p(indx2(i),indx1(i))=-cand_p(indx2(i),indx1(i))
215 IF(pene(i)/=cand_p(indx2(i),indx1(i)))
216 . cand_p(indx2(i),indx1(i))=
min(cand_p(indx2(i),indx1(i)),
217 . ((one-fiveem2)*cand_p(indx2(i),indx1(i))+fiveem2*pene(i)) )
219 pene(i)=
max(zero,pene(i)-cand_p(indx2(i),indx1(i)))
220 IF( pene(i)==zero ) stif(i) = zero
224 IF(cand_p(indx2(i),indx1(i)) < zero)
THEN
227 cand_p(indx2(i),indx1(i)) = -cand_p(indx2(i),indx1(i))
228 IF(pene(i)/=cand_p(indx2(i),indx1(i)))
229 . cand_p(indx2(i),indx1(i)) =
min(cand_p(indx2(i),indx1(i)),
230 . ((one-fiveem2)*cand_p(indx2(i),indx1(i))+fiveem2*pene(i)) )
231 cand_p(indx2(i),indx1(i)) = -cand_p(indx2(i),indx1(i))
232 IF( pene(i)==zero ) stif(i) = zero
236 IF(pene(i)/=cand_p(indx2(i),indx1(i)))
237 . cand_p(indx2(i),indx1(i))=
min(cand_p(indx2(i),indx1(i)),
238 . ((one-fiveem2)*cand_p(indx2(i),indx1(i))+fiveem2*pene(i)) )
240 pene(i)=
max(zero,pene(i)-cand_p(indx2(i),indx1(i)))
241 IF( pene(i)==zero ) stif(i) = zero
248 vx(i) = hs1(i)*vxs1(i) + hs2(i)*vxs2(i)
249 . - hm1(i)*vxm1(i) - hm2(i)*vxm2(i)
250 vy(i) = hs1(i)*vys1(i) + hs2(i)*vys2(i)
251 . - hm1(i)*vym1(i) - hm2(i)*vym2(i)
252 vz(i) = hs1(i)*vzs1(i) + hs2(i)*vzs2(i)
253 . - hm1(i)*vzm1(i) - hm2(i)*vzm2(i)
254 vn(i) = nx(i)*vx(i) + ny(i)*vy(i) + nz(i)*vz(i)
259 stif(i)= half*stif(i)
260 fni(i) = -stif(i) * pene(i)
261 econvt = econvt+fni(i)*vn(i)*dt1
268 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))
THEN
270 fac = stif(i) /
max(em30,stif(i))
275 vis2(i) = two * stif(i) *
min(mas2,masm)
277 ff = fac * visc * vis
278 stif(i) = stif0(i) + ff * dt1inv
279 stif(i) =
max(stif(i) ,fac*sqrt(viscf)*vis*dt1inv)
281 econvt = econvt + ff * vn(i) * dt1
287 fac = stif(i) /
max(em30,stif(i))
292 vis2(i) = two * stif(i) *
min(mas2,masm)
294 c(i)= fac * visc * vis
296 stif(i) = stif(i) + c(i) * dt1inv
298 econvt = econvt + ff * vn(i) * dt1
300 cf(i) = fac*sqrt(viscf)*vis
301 stif(i) =
max(stif(i) ,cf(i)*dt1inv)
306 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))
THEN
308 fac = stif(i) /
max(em30,stif(i))
318 vis2(i) = two* stif(i) * masm * mas2 /
319 .
max(em30,masm+mas2)
321 ff = fac * visc * vis
322 stif(i) = stif0(i) + ff * dt1inv
323 stif(i) =
max(stif(i) ,fac*sqrt(viscf)*vis*dt1inv)
325 econvt = econvt + ff * vn(i) * dt1
331 fac = stif(i) /
max(em30,stif(i))
336 vis2(i) = two* stif(i) * masm * mas2 /
337 .
max(em30,masm+mas2)
339 c(i)= fac * visc * vis
341 stif(i) = stif(i) + c(i) * dt1inv
343 econvt = econvt + ff * vn(i) * dt1
345 cf(i) = fac*sqrt(viscf)*vis
346 stif(i) =
max(stif(i) ,cf(i)*dt1inv)
355 fac = stif(i) /
max(em30,stif(i))
360 vis2(i) = two * stif(i) *
min(mas2,masm)
362 ff = fac * visc * vis
363 stif(i) = stif0(i) + two * ff * dt1inv
364 stif(i) =
max(stif(i) ,fac*sqrt(viscf)*vis*dt1inv)
366 econvt = econvt + ff * vn(i) * dt1
374 fac = stif(i) /
max(em30,stif(i))
379 vis2(i) = two * stif(i) *
min(mas2,masm)
381 ff = fac * visc * vis
382 stif(i) = stif0(i) + two* ff * dt1inv
383 stif(i) =
max(stif(i) ,fac*sqrt(viscf)*vis*dt1inv)
385 econvt = econvt + ff * vn(i) * dt1
393 fac = stif(i) /
max(em30,stif(i))
398 vis2(i) = two * stif(i) *
min(mas2,masm)
400 stif(i) =
max(stif(i) ,fac*sqrt(viscf)*vis*dt1inv)
409 fac = stif(i) /
max(em30,stif(i))
414 vis2(i) = two* stif(i) * masm * mas2 /
415 .
max(em30,masm+mas2)
417 .
max(em30,masm+mas2)
418 stif(i) =
max(stif0(i) ,fac*sqrt(viscf*vis2(i))*dt1inv)
420 econvt = econvt + ff * vn(i) * dt1
421 fni(i) =
min(fni(i),ff)
439 IF(pene(i) == zero)cycle
441 ims2 = bitget(mbinflg(ie),1)
452 fsav11=fsav11-fni(i)*dt12
459 fsav8 =fsav8 +abs(impx)
460 fsav9 =fsav9 +abs(impy)
461 fsav10=fsav10+abs(impz)
462 IF(isensint(1)/=0)
THEN
464 fsavparit(1,1,i) = -fxi(i)
465 fsavparit(1,2,i) = -fyi(i)
466 fsavparit(1,3,i) = -fzi(i)
468 fsavparit(1,1,i) = fxi(i)
469 fsavparit(1,2,i) = fyi(i)
470 fsavparit(1,3,i) = fzi(i)
476 IF(pene(i) == zero)cycle
486 fsav11=fsav11-fni(i)*dt12
487 fsav8 =fsav8 +abs(impx)
488 fsav9 =fsav9 +abs(impy)
489 fsav10=fsav10+abs(impz)
490 IF(isensint(1)/=0)
THEN
491 fsavparit(1,1,i) = fxi(i)
492 fsavparit(1,2,i) = fyi(i)
493 fsavparit(1,3,i) = fzi(i)
499 fsav(1)=fsav(1)+fsav1
500 fsav(2)=fsav(2)+fsav2
501 fsav(3)=fsav(3)+fsav3
502 fsav(8)=fsav(8)+fsav8
503 fsav(9)=fsav(9)+fsav9
504 fsav(10)=fsav(10)+fsav10
505 fsav(11)=fsav(11)+fsav11
506#include "lockoff.inc"
514 fsavsub1(j,jsub)=zero
519 IF(pene(i) == zero)cycle
524 IF (msegtyp(cm_loc(i)) < 0)
THEN
525 ie= - msegtyp(cm_loc(i))
529 IF(ie > nrtm) ie=ie-nrtm
533 DO WHILE(jj<addsube(il+1))
535 itypsub = typsub(jsub)
537 IF(itypsub == 1 )
THEN
539 iss1 = bitget(inflg_sube(jj),0)
540 iss2 = bitget(inflg_sube(jj),1)
542 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
543 ims1 = bitget(inflg_subm(kk),0)
544 ims2 = bitget(inflg_subm(kk),1)
546 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
547 . (ims2 == 1 .AND. iss1 == 1)))
THEN
557 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
558 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
559 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
560 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
562 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
563 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
564 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
565 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
568 IF(isensint(jsub+1)/=0)
THEN
570 fsavparit(jsub+1,1,i) = -fxi(i)
571 fsavparit(jsub+1,2,i) = -fyi(i)
572 fsavparit(jsub+1,3,i) = -fzi(i)
574 fsavparit(jsub+1,1,i) = fxi(i)
575 fsavparit(jsub+1,2,i) = fyi(i)
576 fsavparit(jsub+1,3,i) = fzi(i)
580 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
581 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
582 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
590 ELSEIF(itypsub == 2 )
THEN
597 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
598 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
599 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
601 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
602 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
603 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
605 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
607 IF(isensint(jsub+1)/=0)
THEN
608 fsavparit(jsub+1,1,i) = fxi(i)
609 fsavparit(jsub+1,2,i) = fyi(i)
610 fsavparit(jsub+1,3,i) = fzi(i)
615 ELSEIF(itypsub == 3 )
THEN
617 iss2 = bitget(inflg_sube(jj),0)
618 iss1 = bitget(inflg_sube(jj),1)
620 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
621 ims2 = bitget(inflg_subm(kk),0)
622 ims1 = bitget(inflg_subm(kk),1)
624 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
625 . (ims2 == 1 .AND. iss1 == 1)))
THEN
632 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
633 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
634 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
635 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
637 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
638 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
639 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
640 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
643 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
644 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
645 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
647 IF(isensint(jsub+1)/=0)
THEN
649 fsavparit(jsub+1,1,i) = -fxi(i)
650 fsavparit(jsub+1,2,i) = -fyi(i)
651 fsavparit(jsub+1,3,i) = -fzi(i)
653 fsavparit(jsub+1,1,i) = fxi(i)
654 fsavparit(jsub+1,2,i) = fyi(i)
655 fsavparit(jsub+1,3,i) = fzi(i)
672 IF (msegtyp(cm_loc(i)) < 0)
THEN
673 ie= - msegtyp(cm_loc(i))
677 IF(ie > nrtm) ie=ie-nrtm
680 DO WHILE(kk<addsube(ie+1))
682 itypsub = typsub(ksub)
683 IF(itypsub == 2 )
THEN
689 fsavsub1(1,ksub)=fsavsub1(1,ksub)+impx
690 fsavsub1(2,ksub)=fsavsub1(2,ksub)+impy
691 fsavsub1(3,ksub)=fsavsub1(3,ksub)+impz
693 fsavsub1(8,ksub) =fsavsub1(8,ksub) +abs(impx)
694 fsavsub1(9,ksub) =fsavsub1(9,ksub) +abs(impy)
695 fsavsub1(10,ksub)=fsavsub1(10,ksub)+abs(impz)
697 fsavsub1(11,ksub)=fsavsub1(11,ksub)-fni(i)*dt12
699 IF(isensint(ksub+1)/=0)
THEN
700 fsavparit(ksub+1,1,i) = -fxi(i)
701 fsavparit(ksub+1,2,i) = -fyi(i)
702 fsavparit(ksub+1,3,i) = -fzi(i)
713 IF(pene(i) == zero)cycle
718 IF (msegtyp(cm_loc(i)) < 0)
THEN
719 ie= - msegtyp(cm_loc(i))
723 IF(ie > nrtm) ie=ie-nrtm
730 itypsub = typsub(jsub)
732 IF(itypsub == 1 )
THEN
736 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
737 ims1 = bitget(inflg_subm(kk),0)
738 ims2 = bitget(inflg_subm(kk),1)
740 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
741 . (ims2 == 1 .AND. iss1 == 1)))
THEN
751 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
752 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
753 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
756 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
757 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
758 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
759 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
762 IF(isensint(jsub+1)/=0)
THEN
764 fsavparit(jsub+1,1,i) = -fxi(i)
765 fsavparit(jsub+1,2,i) = -fyi(i)
766 fsavparit(jsub+1,3,i) = -fzi(i)
768 fsavparit(jsub+1,1,i) = fxi(i)
769 fsavparit(jsub+1,2,i) = fyi(i)
770 fsavparit(jsub+1,3,i) = fzi(i)
774 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
775 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
776 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
784 ELSEIF(itypsub == 2 )
THEN
791 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
792 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
793 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
795 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
796 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
797 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
799 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
801 IF(isensint(jsub+1)/=0)
THEN
802 fsavparit(jsub+1,1,i) = fxi(i)
803 fsavparit(jsub+1,2,i) = fyi(i)
804 fsavparit(jsub+1,3,i) = fzi(i)
809 ELSEIF(itypsub == 3 )
THEN
814 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
815 ims2 = bitget(inflg_subm(kk),0)
816 ims1 = bitget(inflg_subm(kk),1)
818 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
819 . (ims2 == 1 .AND. iss1 == 1)))
THEN
830 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
831 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
832 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
833 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
835 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
836 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
837 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
838 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
841 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
842 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
843 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
845 IF(isensint(jsub+1)/=0)
THEN
847 fsavparit(jsub+1,1,i) = -fxi(i)
848 fsavparit(jsub+1,2,i) = -fyi(i)
849 fsavparit(jsub+1,3,i) = -fzi(i)
851 fsavparit(jsub+1,1,i) = fxi(i)
852 fsavparit(jsub+1,2,i) = fyi(i)
853 fsavparit(jsub+1,3,i) = fzi(i)
890 IF(pene(i) == zero)cycle
892 fx = stif0(i)*vx(i)*dt12
893 fy = stif0(i)*vy(i)*dt12
894 fz = stif0(i)*vz(i)*dt12
896 fx = cand_fx(indx2(i),indx1(i)) + fx
897 fy = cand_fy(indx2(i),indx1(i)) + fy
898 fz = cand_fz(indx2(i),indx1(i)) + fz
900 ftn = fx*nx(i) + fy*ny(i) + fz*nz(i)
904 ft = fx*fx + fy*fy + fz*fz
907 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
908 beta =
min(one,fricc(i)*sqrt(fn/ft))
913 cand_fx(indx2(i),indx1(i)) = fxt(i)
914 cand_fy(indx2(i),indx1(i)) = fyt(i)
915 cand_fz(indx2(i),indx1(i)) = fzt(i)
917 fxi(i)=fxi(i) + fxt(i)
918 fyi(i)=fyi(i) + fyt(i)
919 fzi(i)=fzi(i) + fzt(i)
923 fsav4 = fsav4 + fxt(i)*dt12
924 fsav5 = fsav5 + fyt(i)*dt12
925 fsav6 = fsav6 + fzt(i)*dt12
927 fsav12 = fsav12 + abs(fxi(i)*dt12)
928 fsav13 = fsav13 + abs(fyi(i)*dt12)
929 fsav14 = fsav14 + abs(fzi(i)*dt12)
930 fsav15 = fsav15 + sqrt(fxi(i)*fxi(i)+fyi(i)*fyi(i)+fzi(i)*fzi(i))*dt12
932 . + dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
938 fsav(4) = fsav(4) + fsav4
939 fsav(5) = fsav(5) + fsav5
940 fsav(6) = fsav(6) + fsav6
941 fsav(12) = fsav(12) + fsav12
942 fsav(13) = fsav(13) + fsav13
943 fsav(14) = fsav(14) + fsav14
944 fsav(15) = fsav(15) + fsav15
945#include "lockoff.inc"
954 IF(pene(i) == zero)cycle
959 IF (msegtyp(cm_loc(i)) < 0)
THEN
960 ie= - msegtyp(cm_loc(i))
964 IF(ie > nrtm) ie=ie-nrtm
968 DO WHILE(jj<addsube(il+1))
971 itypsub = typsub(jsub)
973 IF(itypsub == 1 )
THEN
975 iss1 = bitget(inflg_sube(jj),0)
976 iss2 = bitget(inflg_sube(jj),1)
978 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
979 ims1 = bitget(inflg_subm(kk),0)
980 ims2 = bitget(inflg_subm(kk),1)
982 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
983 . (ims2 == 1 .AND. iss1 == 1)))
THEN
992 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
993 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
994 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
999 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
1000 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
1001 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
1003 IF(isensint(jsub+1)/=0)
THEN
1004 fsavparit(jsub+1,4,i) = fxt(i)
1005 fsavparit(jsub+1,5,i) = fyt(i)
1006 fsavparit(jsub+1,6,i) = fzt(i)
1009 fsavsub1(15,jsub)= fsavsub1(15,jsub)
1010 . +sqrt(impx*impx+impy*impy+impz*impz)
1024 ELSEIF(itypsub == 2 )
THEN
1030 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
1031 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
1032 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
1037 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
1038 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
1039 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
1041 IF(isensint(jsub+1)/=0)
THEN
1042 fsavparit(jsub+1,4,i) = fxt(i)
1043 fsavparit(jsub+1,5,i) = fyt(i)
1044 fsavparit(jsub+1,6,i) = fzt(i)
1047 fsavsub1(15,jsub)= fsavsub1(15,jsub)
1048 . +sqrt(impx*impx+impy*impy+impz*impz)
1053 ELSEIF(itypsub == 3 )
THEN
1055 iss2 = bitget(inflg_sube(jj),0)
1056 iss1 = bitget(inflg_sube(jj),1)
1058 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
1059 ims2 = bitget(inflg_subm(kk),0)
1060 ims1 = bitget(inflg_subm(kk),1)
1062 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
1063 . (ims2 == 1 .AND. iss1 == 1)))
THEN
1074 fsavsub1(4,jsub)=fsavsub1(4,jsub)-impx
1075 fsavsub1(5,jsub)=fsavsub1(5,jsub)-impy
1076 fsavsub1(6,jsub)=fsavsub1(6,jsub)-impz
1079 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
1080 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
1088 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
1089 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
1090 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
1092 IF(isensint(jsub+1)/=0)
THEN
1094 fsavparit(jsub+1,4,i) = -fxt(i)
1095 fsavparit(jsub+1,5,i) = -fyt(i)
1096 fsavparit(jsub+1,6,i) = -fzt(i)
1098 fsavparit(jsub+1,4,i) = fxt(i)
1099 fsavparit(jsub+1,5,i) = fyt(i)
1100 fsavparit(jsub+1,6,i) = fzt(i)
1104 fsavsub1(15,jsub)= fsavsub1(15,jsub)
1105 . +sqrt(impx*impx+impy*impy+impz*impz)
1118 IF (msegtyp(cm_loc(i)) < 0)
THEN
1119 ie= - msegtyp(cm_loc(i))
1123 IF(ie > nrtm) ie=ie-nrtm
1126 DO WHILE(kk<addsube(ie+1))
1128 itypsub = typsub(ksub)
1129 IF(itypsub == 2 )
THEN
1135 fsavsub1(4,ksub)=fsavsub1(4,ksub)+impx
1136 fsavsub1(5,ksub)=fsavsub1(5,ksub)+impy
1137 fsavsub1(6,ksub)=fsavsub1(6,ksub)+impz
1142 fsavsub1(12,ksub)=fsavsub1(12,jsub)+abs(impx)
1143 fsavsub1(13,ksub)=fsavsub1(13,jsub)+abs(impy)
1144 fsavsub1(14,ksub)=fsavsub1(14,jsub)+abs(impz)
1146 IF(isensint(ksub+1)/=0)
THEN
1147 fsavparit(ksub+1,4,i) = -fxt(i)
1148 fsavparit(ksub+1,5,i) = -fyt(i)
1149 fsavparit(ksub+1,6,i) = -fzt(i)
1152 fsavsub1(15,ksub)= fsavsub1(15,ksub)
1153 . +sqrt(impx*impx+impy*impy+impz*impz)
1163 IF(pene(i) == zero)cycle
1168 IF (msegtyp(cm_loc(i)) < 0)
THEN
1169 ie= - msegtyp(cm_loc(i))
1173 IF(ie > nrtm) ie=ie-nrtm
1179 itypsub = typsub(jsub)
1181 IF(itypsub == 1 )
THEN
1186 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
1187 ims1 = bitget(inflg_subm(kk),0)
1188 ims2 = bitget(inflg_subm(kk),1)
1190 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
1191 . (ims2 == 1 .AND. iss1 == 1)))
THEN
1200 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
1201 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
1202 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
1207 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
1208 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
1209 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
1211 IF(isensint(jsub+1)/=0)
THEN
1212 fsavparit(jsub+1,4,i) = fxt(i)
1213 fsavparit(jsub+1,5,i) = fyt(i)
1214 fsavparit(jsub+1,6,i) = fzt(i)
1217 fsavsub1(15,jsub)= fsavsub1(15,jsub)
1218 . +sqrt(impx*impx+impy*impy+impz*impz)
1233 ELSEIF(itypsub == 2 )
THEN
1239 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
1240 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
1241 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
1246 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
1247 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
1248 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
1250 IF(isensint(jsub+1)/=0)
THEN
1251 fsavparit(jsub+1,4,i) = fxt(i)
1252 fsavparit(jsub+1,5,i) = fyt(i)
1253 fsavparit(jsub+1,6,i) = fzt(i)
1256 fsavsub1(15,jsub)= fsavsub1(15,jsub)
1257 . +sqrt(impx*impx+impy*impy+impz*impz)
1262 ELSEIF(itypsub == 3 )
THEN
1267 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
1268 ims2 = bitget(inflg_subm(kk),0)
1269 ims1 = bitget(inflg_subm(kk),1)
1271 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
1272 . (ims2 == 1 .AND. iss1 == 1)))
THEN
1283 fsavsub1(4,jsub)=fsavsub1(4,jsub)-impx
1284 fsavsub1(5,jsub)=fsavsub1(5,jsub)-impy
1285 fsavsub1(6,jsub)=fsavsub1(6,jsub)-impz
1288 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
1289 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
1290 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
1296 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
1297 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
1298 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
1300 IF(isensint(jsub+1)/=0)
THEN
1302 fsavparit(jsub+1,4,i) = -fxt(i)
1303 fsavparit(jsub+1,5,i) = -fyt(i)
1304 fsavparit(jsub+1,6,i) = -fzt(i)
1306 fsavparit(jsub+1,4,i) = fxt(i)
1307 fsavparit(jsub+1,5,i) = fyt(i)
1308 fsavparit(jsub+1,6,i) = fzt(i)
1312 fsavsub1(15,jsub)= fsavsub1(15,jsub)
1313 . +sqrt(impx*impx+impy*impy+impz*impz)
1327#include "lockon.inc"
1331 fsavsub(j,nsub)=fsavsub(j,nsub)+fsavsub1(j,jsub)
1333 fsavsub(22,nsub)=fsavsub(22,nsub)+fsavsub1(22,jsub)
1334 fsavsub(23,nsub)=fsavsub(23,nsub)+fsavsub1(23,jsub)
1335 fsavsub(24,nsub)=fsavsub(24,nsub)+fsavsub1(24,jsub)
1337#include "lockoff.inc"
1341#include "lockon.inc"
1342 econtv = econtv + econvt
1343 econt = econt + econtt
1344#include "lockoff.inc"
1349 IF(pene(i) == zero)cycle
1351 fx1(i)=-fxi(i)*hs1(i)
1352 fy1(i)=-fyi(i)*hs1(i)
1353 fz1(i)=-fzi(i)*hs1(i)
1355 fx2(i)=-fxi(i)*hs2(i)
1356 fy2(i)=-fyi(i)*hs2(i)
1357 fz2(i)=-fzi(i)*hs2(i)
1359 fx3(i)=fxi(i)*hm1(i)
1360 fy3(i)=fyi(i)*hm1(i)
1361 fz3(i)=fzi(i)*hm1(i)
1363 fx4(i)=fxi(i)*hm2(i)
1364 fy4(i)=fyi(i)*hm2(i)
1365 fz4(i)=fzi(i)*hm2(i)
1370 stif(i) = two*stif(i)
1374 IF(kdtint==1.OR.idtmins==2)
THEN
1376 . .AND.(ivis2==0.OR.ivis2==1))
THEN
1380 IF(ms1(i)==zero)
THEN
1384 k1(i)=kt(i)*abs(hs1(i))
1385 c1(i)=c(i)*abs(hs1(i))
1386 cx =four*c1(i)*c1(i)
1387 cy =eight*ms1(i)*k1(i)
1388 aux = sqrt(cx+cy)+two*c1(i)
1389 st1(i)= k1(i)*aux*aux/
max(cy,em30)
1390 cfi = cf(i)*abs(hs1(i))
1391 aux = two*cfi*cfi/
max(ms1(i),em20)
1398 IF(ms2(i)==zero)
THEN
1402 k2(i)=kt(i)*abs(hs2(i))
1403 c2(i)=c(i)*abs(hs2(i))
1404 cx =four*c2(i)*c2(i)
1405 cy =eight*ms2(i)*k2(i)
1406 aux = sqrt(cx+cy)+two*c2(i)
1407 st2(i)= k2(i)*aux*aux/
max(cy,em30)
1408 cfi = cf(i)*abs(hs2(i))
1409 aux = two*cfi*cfi/
max(ms2(i),em20)
1416 IF(mm1(i)==zero)
THEN
1420 k3(i)=kt(i)*abs(hm1(i))
1421 c3(i)=c(i)*abs(hm1(i))
1422 cx =four*c3(i)*c3(i)
1423 cy =eight*mm1(i)*k3(i)
1424 aux = sqrt(cx+cy)+two*c3(i)
1425 st3(i)= k3(i)*aux*aux/
max(cy,em30)
1426 cfi = cf(i)*abs(hm1(i))
1427 aux = two*cfi*cfi/
max(mm1(i),em20)
1434 IF(mm2(i)==zero)
THEN
1440 cx =four*c4(i)*c4(i)
1442 aux = sqrt(cx+cy)+two*c4(i)
1443 st4(i)= k4(i)*aux*aux/
max(cy,em30)
1444 cfi = cf(i)*abs(hm2(i))
1445 aux = two*cfi*cfi/
max(mm2(i),em20)
1454 k1(i) =stif(i)*abs(hs1(i))
1456 k2(i) =stif(i)*abs(hs2(i))
1458 k3(i) =stif(i)*abs(hm1(i))
1460 k4(i) =stif(i)*abs(hm2(i))
1468 IF(nintloadp > 0)
THEN
1469 DO k = kloadpinter(nin)+1, kloadpinter(nin+1)
1471 ppl = loadp_hyd_inter(pp)
1472 dgapload = dgaploadint(k)
1474 IF(pene(i) > zero .OR.dist(i) <= dgapload)
THEN
1476 tagncont(ppl,m1(i)) = 1
1477 tagncont(ppl,m2(i)) = 1
1478 IF(cs_loc(i)<=nedge)
THEN
1480 tagncont(ppl,n1(i)) = 1
1481 tagncont(ppl,n2(i)) = 1
1495 2 m2 ,hs1 ,hs2 ,hm1 ,hm2 ,
1496 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,
1497 4 fz2 ,fx3 ,fy3 ,fz3 ,fx4 ,
1498 5 fy4 ,fz4 ,a ,stifn,stif ,
1499 6 nedge,nin ,jtask,pene )
1502 2 m2 ,hs1 ,hs2 ,hm1 ,hm2 ,
1503 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,
1504 4 fz2 ,fx3 ,fy3 ,fz3 ,fx4 ,
1505 5 fy4 ,fz4 ,a ,stifn,nedge,
1506 6 k1 ,k2 ,k3 ,k4 ,c1 ,
1507 7 c2 ,c3 ,c4 ,viscn,nin ,
1513 2 m2 ,hs1 ,hs2 ,hm1 ,hm2 ,
1514 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,
1515 4 fz2 ,fx3 ,fy3 ,fz3 ,fx4 ,
1516 5 fy4 ,fz4 ,fskyi ,isky ,niskyfie,
1517 6 stif ,nedge ,nin ,noint ,pene ,
1521 2 m2 ,hs1 ,hs2 ,hm1 ,hm2 ,
1522 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,
1523 4 fz2 ,fx3 ,fy3 ,fz3 ,fx4 ,
1524 5 fy4 ,fz4 ,isky ,niskyfie,nedge ,
1525 6 k1 ,k2 ,k3 ,k4 ,c1 ,
1526 7 c2 ,c3 ,c4 ,nin , noint,
1534 2 m2 ,hs1 ,hs2 ,hm1 ,hm2 ,
1535 3 stif ,nin ,noint ,mskyi_sms ,iskyi_sms,
1536 4 nsms ,k1 ,k2 ,k3 ,k4 ,
1537 5 c1 ,c2 ,c3 ,c4 ,nedge , edge_id)
1542#include "mic_lockon.inc"
1545 assert(i <= 4*mvsiz)
1546 assert(cs_loc(i) > 0)
1547 printif(cs_loc(i) < 0,i)
1548 printif(cs_loc(i) < 0,cs_loc(i))
1549 IF(cs_loc(i)>nedge)
THEN
1550 ni = cs_loc(i)-nedge
1553 IF(pene(i) /= 0.OR.tagip(i)==1)
THEN
1558#include "mic_lockoff.inc"
1561 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT>0)
THEN
1562#include "lockon.inc"
1566 IF(pene(i) == zero)cycle
1568 IF(cs_loc(i)<=nedge)
THEN
1569 fcont(1,n1(i)) =fcont(1,n1(i)) + fx1(i)
1570 fcont(2,n1(i)) =fcont(2,n1(i)) + fy1(i)
1571 fcont(3,n1(i)) =fcont(3,n1(i)) + fz1(i)
1572 fcont(1,n2(i)) =fcont(1,n2(i)) + fx2(i)
1573 fcont(2,n2(i)) =fcont(2,n2(i)) + fy2(i)
1574 fcont(3,n2(i)) =fcont(3,n2(i)) + fz2(i)
1576 fcont(1,m1(i)) =fcont(1,m1(i)) + fx3(i)
1577 fcont(2,m1(i)) =fcont(2,m1(i)) + fy3(i)
1578 fcont(3,m1(i)) =fcont(3,m1(i)) + fz3(i)
1579 fcont(1,m2(i)) =fcont(1,m2(i)) + fx4(i)
1580 fcont(2,m2(i)) =fcont(2,m2(i)) + fy4(i)
1581 fcont(3,m2(i)) =fcont(3,m2(i)) + fz4(i)
1584#include "lockoff.inc"
1590 IF(nstrf(1)+nstrf(2)/=0)
THEN
1592 nbinter=nstrf(k0+14)
1595 IF(nstrf(k1s)==noint)
THEN
1597#include "lockon.inc"
1600 IF(pene(k) == zero)cycle
1602 IF(cs_loc(k)<=nedge)
THEN
1603 IF(secfcum(4,n1(k),i)==1.)
THEN
1604 secfcum(1,n1(k),i)=secfcum(1,n1(k),i)-fx1(k)
1605 secfcum(2,n1(k),i)=secfcum(2,n1(k),i)-fy1
1606 secfcum(3,n1(k),i)=secfcum(3,n1
1608 IF(secfcum(4,n2(k),i)==1.)
THEN
1609 secfcum(1,n2(k),i)=secfcum(1,n2(k),i)-fx2(k)
1610 secfcum(2,n2(k),i)=secfcum(2,n2(k),i)-fy2(k)
1611 secfcum(3,n2(k),i)=secfcum(3,n2(k),i)-fz2(k)
1614 IF(secfcum(4,m1(k),i)==1.)
THEN
1615 secfcum(1,m1(k),i)=secfcum(1,m1(k),i)-fx3(k)
1616 secfcum(2,m1(k),i)=secfcum(2,m1(k),i)-fy3(k)
1617 secfcum(3,m1(k),i)=secfcum(3,m1(k),i)-fz3(k)
1619 IF(secfcum(4,m2(k),i)==1.)
THEN
1620 secfcum(1,m2(k),i)=secfcum(1,m2(k),i)-fx4(k)
1621 secfcum(2,m2(k),i)=secfcum(2,m2(k),i)-fy4(k)
1622 secfcum(3,m2(k),i)=secfcum(3,m2(k),i)-fz4(k)
1625#include "lockoff.inc"