73
74
75
79
80
81
82#include "implicit_f.inc"
83#include "comlock.inc"
84
85
86
87#include "mvsiz_p.inc"
88
89
90
91#include "com01_c.inc"
92#include "com04_c.inc"
93#include "com06_c.inc"
94#include "com08_c.inc"
95#include "scr05_c.inc"
96#include "scr07_c.inc"
97#include "scr11_c.inc"
98#include "scr14_c.inc"
99#include "scr16_c.inc"
100#include "scr18_c.inc"
101#include "units_c.inc"
102#include "parit_c.inc"
103#include "param_c.inc"
104#include "impl1_c.inc"
105#include "sms_c.inc"
106#include "kincod_c.inc"
107
108
109
110 INTEGER NELTST,ITYPTST,JLT,IBCC,IBCM,IBCS,IVIS2,INACTI,IBAG,NIN,
111 . NTY ,NLN,NLG(NLN),NSV(*),
112 . ICODT(*), ITAB(*), ISKY(*), KINET(*),
113 . MFROT, IFQ, NOINT,NEWFRONT,ISECIN, NSTRF(*),
114 . IFPEN(*) ,ICONTACT(*), CAND_N(*),
115 . KINI(*),
116 . ISET, NISKYFI,IADM,INTTH,IFORM, IGAP,JTASK
117 INTEGER IX1L(MVSIZ), IX2L(MVSIZ), IX3L(MVSIZ), IX4L(MVSIZ),
118 . CN_LOC(MVSIZ),CE_LOC(MVSIZ),INDEX(MVSIZ),NSVG(MVSIZ),
119 . NISUB, LISUB(*), ADDSUBS(*), ADDSUBM(*), LISUBS(*),
120 . LISUBM(*),ILAGM,ICURV(3),
121 . ISKYI_SMS(*), NSMS(*), ISENSINT(*),
123 . stiglo,frot_p(*), x(3,*), xa(3,*),dxanc(3,*),
124 . a(3,*), ms(*), va(3,*), fsav(*),fcont(3,*),
125 . cand_fx(*),cand_fy(*),cand_fz(*),alpha0,
126 . gap, fric,visc,viscf,vis,dt2t,stfa(*),stifn(*),
127 . fskyi(lskyi,nfskyi),fsavsub(nthvki,*), fncont(3,*),ftcont(3,*),
128 . mskyi_sms(*)
130 . nx1(mvsiz), nx2(mvsiz), nx3(mvsiz), nx4(mvsiz),
131 . ny1(mvsiz), ny2(mvsiz), ny3(mvsiz), ny4(mvsiz),
132 . nz1(mvsiz), nz2(mvsiz), nz3(mvsiz), nz4(mvsiz),
133 . lb1(mvsiz), lb2(mvsiz), lb3(mvsiz), lb4(mvsiz),
134 . lc1(mvsiz), lc2(mvsiz), lc3(mvsiz), lc4(mvsiz),
135 . p1(mvsiz), p2(mvsiz), p3(mvsiz), p4(mvsiz), stif(mvsiz),
136 . gapv(mvsiz),gapr(mvsiz),secfcum(7,numnod,nsect),
137 . tmp(mvsiz),stifsav(mvsiz), viscn(*),
138 . vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),msi(mvsiz),
139 . x1(mvsiz),y1(mvsiz),z1(mvsiz),
140 . x2(mvsiz),y2(mvsiz),z2(mvsiz),
141 . x3(mvsiz),y3(mvsiz),z3(mvsiz),
142 . x4(mvsiz),y4(mvsiz),z4(mvsiz),
143 . xi(mvsiz),yi(mvsiz),zi(mvsiz),penis(2,*),penim(2,*),
144 . phi(mvsiz), fthe(*),ftheskyi(lskyi),temp(*), tempi(mvsiz),
145 . rstif,fsavparit(nisub+1,11,*)
147 . nod_normal(3,*), rcurvi(*), rcontact(*), acontact(*),
148 . pcontact(*),padm, anglmi(*),gap_s(*),alphak(3,*),cmaj(mvsiz)
149 double precision
150 . daanc6(3,6,*)
151 TYPE(H3D_DATABASE) :: H3D_DATA
152
153
154
155 INTEGER I,J1,IG,J,JG,IM,IS,K0,NBINTER,K1S,K,IL,IE,NN,NI,NA1,NA2,
156 . JSUB,KSUB,,KK,,NSUB,ISIGN,IPROJ,IBID
157 INTEGER IX1G(MVSIZ), IX2G(MVSIZ), IX3G(MVSIZ), IX4G(MVSIZ)
159 . fxr(mvsiz), fyr(mvsiz), fzr(mvsiz),
160 . fxi(mvsiz), fyi(mvsiz), fzi(mvsiz), fni(mvsiz),
161 . fxt(mvsiz),fyt(mvsiz),fzt(mvsiz),
162 . fx1(mvsiz), fx2(mvsiz), fx3(mvsiz), fx4(mvsiz),
163 . fy1(mvsiz), fy2(mvsiz), fy3(mvsiz), fy4(mvsiz),
164 . fz1(mvsiz), fz2(mvsiz), fz3(mvsiz), fz4(mvsiz),
165 . n1(mvsiz), n2(mvsiz), n3(mvsiz
166 . vis2(mvsiz), dtmi(mvsiz), xmu(mvsiz),stif0(mvsiz),
167 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
168 . vx(mvsiz), vy(mvsiz), vz(mvsiz), vn(mvsiz),
169 . st1(mvsiz),st2(mvsiz),st3(mvsiz),st4(mvsiz),stv(mvsiz),
170 . kt(mvsiz),c(mvsiz),cf(mvsiz),
171 . ks(mvsiz),k1(mvsiz),k2(mvsiz),k3(mvsiz),k4(mvsiz),
172 . cs(mvsiz),c1(mvsiz),c2(mvsiz),c3(mvsiz),c4(mvsiz),
173 . p1s(mvsiz),p2s(mvsiz),p3s(mvsiz),p4s(mvsiz),
174 . phi1(mvsiz),phi2(mvsiz),phi3(mvsiz),phi4(mvsiz),
175 . fsavsub1(15,nisub),masm(mvsiz)
177 . vnx, vny, vnz, aa, crit,s2,dist,rdist,
178 . v2, fm2, dt1inv, visca, fac,ff,alphi,
alpha,beta,
179 . fx, fy, fz, f2, mas2, m2sk, dtmi0,dti,ft,fn,fmax,ftn,
180 . facm1, econtt, econvt, h0, la1, la2, la3, la4,
181 . d1,d2,d3,d4,a1,a2,a3,a4,e10, h0d, s2d, sum,
182 . fsav1, fsav2, fsav3, fsav4, fsav5, fsav6, fsav7, fsav8,
183 . fsav9, fsav10, fsav11, fsav12, fsav13, fsav14, fsav15, ffo,
184 . la1d,la2d,la3d,la4d,t1,t1d,t2,t2d,ffd,visd,facd,d1d,
185 . d2d,d3d,d4d,vnxd,vnyd,vnzd,v2d,fm2d,f2d,aad,fxd,fyd,fzd,
186 . a1d,a2d,a3d,a4d,vv,ax1,ax2,ay1,ay2,az1,az2,ax,ay,az,
187 .
area,p,vv1,vv2,v21,dmu, dti2,h00 ,a0x,a0y,a0z,rx,ry,rz,
188 . anx,any,anz,aan,aax,aay,aaz ,rr,rs,aaa,stfr,visr,
189 . prec,ps,xsa,pis,pplus,cx,cy,cfi,aux,tm,ts,impx,impy,impz,bb,
190 . nn1,nn2,nn3,nn4,xn1,yn1,zn1,xn2,yn2,zn2,xn3,yn3,zn3,xn4,yn4,
191 . zn4,dtmini,bid
192
193 DOUBLE PRECISION FX6(6,MVSIZ), FY6(6,MVSIZ), FZ6(6,MVSIZ)
194
195
196 IF (iresp==1) THEN
197 prec = fiveem4
198 ELSE
199 prec = em10
200 ENDIF
201 IF(dt1>zero)THEN
202 dt1inv = one/dt1
203 ELSE
204 dt1inv =zero
205 ENDIF
206 econtt = zero
207 econvt = zero
208 DO i=1,jlt
209 stif0(i) = stif(i)
210 ix1g(i) = nlg(ix1l(i))
211 ix2g(i) = nlg(ix2l(i))
212 ix3g(i) = nlg(ix3l(i))
213 ix4g(i) = nlg(ix4l(i))
214 ENDDO
215
216
217
218
219
220
221 IF(icurv(1) == 3) THEN
222 DO i=1,jlt
223
224 bb = gapv(i)+cmaj
225
226 d1 = sqrt(p1(i))
227 p1(i) =
max(zero, bb - d1)
228
229 d2 = sqrt(p2(i))
230 p2(i) =
max(zero, bb - d2)
231
232 d3 = sqrt(p3(i))
233 p3(i) =
max(zero, bb - d3)
234
235 d4 = sqrt(p4(i))
236 p4(i) =
max(zero, bb - d4)
237
238 a1 = p1(i)/
max(em20,d1)
239 a2 = p2(i)/
max(em20,d2)
240 a3 = p3(i)/
max(em20,d3)
241 a4 = p4(i)/
max(em20,d4)
242 n1(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i)
243 n2(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i)
244 n3(i) = a1*nz1(i) + a2*nz2(i) + a3*nz3(i) + a4*nz4(i)
245 ENDDO
246 ELSE
247 DO i=1,jlt
248
249 d1 = sqrt(p1(i))
250 p1(i) =
max(zero, gapv(i) - d1)
251
252 d2 = sqrt(p2(i))
253 p2(i) =
max(zero, gapv(i) - d2)
254
255 d3 = sqrt(p3(i))
256 p3(i) =
max(zero, gapv(i) - d3)
257
258 d4 = sqrt(p4(i))
259 p4(i) =
max(zero, gapv(i) - d4)
260
261 a1 = p1(i)/
max(em20,d1)
262 a2 = p2(i)/
max(em20,d2)
263 a3 = p3(i)/
max(em20,d3)
264 a4 = p4(i)/
max(em20,d4)
265 n1(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i)
266 n2(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i)
267 n3(i) = a1*nz1(i) + a2*nz2(i) + a3*nz3(i) + a4*nz4(i)
268 ENDDO
269 ENDIF
270
271 DO i=1,jlt
272 IF(ix3g(i)/=ix4g(i))THEN
273 pene(i) =
max(p1(i),p2(i),p3(i),p4(i))
274
275 la1 = one - lb1(i) - lc1(i)
276 la2 = one - lb2(i) - lc2(i)
277 la3 = one - lb3(i) - lc3(i)
278 la4 = one - lb4(i) - lc4(i)
279
280 h0 = fourth *
281 . (p1(i)*la1 + p2(i)*la2 + p3(i)*la3 + p4(i)*la4)
282 h1(i) = h0 + p1(i) * lb1(i) + p4(i) * lc4(i)
283 h2(i) = h0 + p2(i) * lb2(i) + p1(i) * lc1(i)
284 h3(i) = h0 + p3(i) * lb3(i) + p2(i) * lc2(i)
285 h4(i) = h0 + p4(i) * lb4(i) + p3(i) * lc3(i)
286
287 h00 = one/
max(em20,h1(i) + h2(i) + h3(i) + h4(i))
288 h1(i) = h1(i) * h00
289 h2(i) = h2(i) * h00
290 h3(i) = h3(i) * h00
291 h4(i) = h4(i) * h00
292
293 ELSE
294 pene(i) = p1(i)
295 n1(i) = nx1(i)
296 n2(i) = ny1(i)
297 n3(i) = nz1(i)
298 h1(i) = lb1(i)
299 h2(i) = lc1(i)
300 h3(i) = one - lb1(i) - lc1(i)
301 h4(i) = zero
302 ENDIF
303 ENDDO
304
305
306
307 IF(icurv(1)==1)THEN
308
309 na1 = icurv(2)
310 DO i=1,jlt
311 rr = 1.e30
312 a0x = xa(1,na1)
313 a0y = xa(2,na1)
314 a0z = xa(3,na1)
315
316 rx = x1(i)-a0x
317 ry = y1(i)-a0y
318 rz = z1(i)-a0z
319 rr =
min(rr , rx*rx + ry*ry + rz*rz)
320 rx = x2(i)-a0x
321 ry = y2(i)-a0y
322 rz = z2(i)-a0z
323 rr =
min(rr , rx*rx + ry*ry + rz*rz)
324 rx = x3(i)-a0x
325 ry = y3(i)-a0y
326 rz = z3(i)-a0z
327 rr =
min(rr , rx*rx + ry*ry + rz*rz)
328 IF(ix3g(i)/=ix4g(i))THEN
329 rx = x4(i)-a0x
330 ry = y4(i)-a0y
331 rz = z4(i)-a0z
332 rr =
min(rr , rx*rx + ry*ry + rz*rz)
333 ENDIF
334 rx = xi(i)-a0x
335 ry = yi(i)-a0y
336 rz = zi(i)-a0z
337 rs = sqrt(rx*rx + ry*ry + rz*rz)
338 rr = sqrt(rr)
339 IF(rs-rr+gapv(i)<0.)THEN
340 stif(i) = 0.
341 pene(i) = 0.
342 ELSEIF(rs-rr+gapv(i)<pene(i))
THEN
343 pene(i) = rs-rr+gapv(i)
344 ENDIF
345 n1(i) = -rx
346 n2(i) = -ry
347 n3(i) = -rz
348 ENDDO
349 ELSEIF(icurv(1)==2)THEN
350
351 na1 = icurv(2)
352 na2 = icurv(3)
353 DO i=1,jlt
354 rr = 1.e30
355 a0x = xa(1,na1)
356 a0y = xa(2,na1)
357 a0z = xa(3,na1)
358 anx = xa(1,na2)-a0x
359 any = xa(2,na2)-a0y
360 anz = xa(3,na2)-a0z
361 aan = 1. / (anx*anx + any*any + anz*anz)
362
363 aax = x1(i)-a0x
364 aay = y1(i)-a0y
365 aaz = z1(i)-a0z
366 aaa = (aax*anx + aay*any + aaz*anz) * aan
367 rx = aax - aaa * anx
368 ry = aay - aaa * any
369 rz = aaz - aaa * anz
370 rr =
min(rr , rx*rx + ry*ry + rz*rz)
371
372 aax = x2(i)-a0x
373 aay = y2(i)-a0y
374 aaz = z2(i)-a0z
375 aaa = (aax*anx + aay*any + aaz*anz) * aan
376 rx = aax - aaa * anx
377 ry = aay - aaa * any
378 rz = aaz - aaa * anz
379 rr =
min(rr , rx*rx + ry*ry + rz*rz)
380
381 aax = x3(i)-a0x
382 aay = y3(i)-a0y
383 aaz = z3(i)-a0z
384 aaa = (aax*anx + aay*any + aaz*anz) * aan
385 rx = aax - aaa * anx
386 ry = aay - aaa * any
387 rz = aaz - aaa * anz
388 rr =
min(rr , rx*rx + ry*ry + rz*rz)
389 IF(ix3g(i)/=ix4g(i))THEN
390
391 aax = x4(i)-a0x
392 aay = y4(i)-a0y
393 aaz = z4(i)-a0z
394 aaa = (aax*anx + aay*any + aaz*anz) * aan
395 rx = aax - aaa * anx
396 ry = aay - aaa * any
397 rz = aaz - aaa * anz
398 rr =
min(rr , rx*rx + ry*ry + rz*rz)
399 ENDIF
400 aax = xi(i)-a0x
401 aay = yi(i)-a0y
402 aaz = zi(i)-a0z
403
404 aaa = (aax*anx + aay*any + aaz*anz) * aan
405 rx = aax - aaa * anx
406 ry = aay - aaa * any
407 rz = aaz - aaa * anz
408 rs = sqrt(rx*rx + ry*ry + rz*rz)
409 rr = sqrt(rr)
410 IF(rs-rr+gapv(i)<0.)THEN
411 stif(i) = 0.
412 pene(i) = 0.
413 ELSEIF(rs-rr+gapv(i)<pene(i))
THEN
414 pene(i) = rs-rr+gapv(i)
415 n1(i) = -rx
416 n2(i) = -ry
417 n3(i) = -rz
418 ELSEIF(rs-rr-gapv(i)>0.)THEN
419 stif(i) = 0.
420 pene(i) = 0.
421 ELSEIF(rs-rr-gapv(i) < pene(i))THEN
422 xn1 = x1(i) - xi(i)
423 yn1 = y1(i) - yi(i)
424 zn1 = z1(i) - zi(i)
425 xn2 = x2(i) - xi(i)
426 yn2 = y2(i) - yi(i)
427 zn2 = z2(i) - zi(i)
428 xn3 = x3(i) - xi(i)
429 yn3 = y3(i) - yi(i)
430 zn3 = z3(i) - zi(i)
431
432 nn1 = (yn1*zn2-yn2*zn1) * rx +
433 . (zn1*xn2-zn2*xn1) * ry +
434 . (xn1*yn2-xn2*yn1) * rz
435 nn2 = (yn2*zn3-yn3*zn2) * rx +
436 . (zn2*xn3-zn3*xn2) * ry +
437 . (xn2*yn3-xn3*yn2) * rz
438 nn3 = (yn3*zn4-yn4*zn3) * rx +
439 . (zn3*xn4-zn4*xn3) * ry +
440 . (xn3*yn4-xn4*yn3) * rz
441 IF(ix3l(i)/=ix4l(i))THEN
442 xn4 = x4(i) - xi(i)
443 yn4 = y4(i) - yi(i)
444 zn4 = z4(i) - zi(i)
445 nn4 = (yn4*zn1-yn1*zn4) * rx +
446 . (zn4*xn1-zn1*xn4) * ry +
447 . (xn4*yn1-xn1*yn4) * rz
448 ELSE
449 xn4 = zero
450 yn4 = zero
451 zn4 = zero
452 nn4=zero
453 ENDIF
454 IF( nn1>=zero .AND. nn2>=zero
455 . .AND. nn3>=zero .AND. nn4>=zero) THEN
456 iproj = 1
457 ELSEIF( nn1<=zero .AND. nn2<=zero
458 . .AND. nn3<=zero .AND. nn4<=zero) THEN
459 iproj = 1
460 ELSE
461 iproj = 0
462 ENDIF
463
464 IF(iproj == 1)THEN
465 pene(i) = -rs+rr+gapv(i)
466 n1(i) = rx
467 n2(i) = ry
468 n3(i) = rz
469 ENDIF
470 ENDIF
471 ENDDO
472
473 ELSEIF(icurv(1) == 3)THEN
474 CALL i7curv(jlt ,pene ,n1 ,n2 ,
475 1 n3 ,gapv ,xa ,nod_normal,
476 2 ix1l ,ix2l ,ix3l ,ix4l ,
477 3 h1 ,h2 ,h3 ,h4 ,
478 4 x1 ,x2 ,x3 ,x4 ,y1 ,
479 5 y2 ,y3 ,y4 ,z1 ,z2
480 6 z3 ,z4 ,xi ,yi ,zi )
481
482 DO i=1,jlt
483 IF(pene(i)<zero)THEN
484 stif(i) =zero
485 pene(i) =zero
486 END IF
487 END DO
488 ENDIF
489
490 DO i=1,jlt
491 s2 = one/
max(em30,sqrt(n1(i)**2 + n2(i)**2 + n3(i)**2))
492 n1(i) = n1(i)*s2
493 n2(i) = n2(i)*s2
494 n3(i) = n3(i)*s2
495 ENDDO
496
497 DO i=1,jlt
498 vx(i) = vxi(i) - h1(i)*va(1,ix1l(i)) - h2(i)*va(1,ix2l(i))
499 . - h3(i)*va(1,ix3l(i)) - h4(i)*va(1,ix4l(i))
500 vy(i) = vyi(i) - h1(i)*va(2,ix1l(i)) - h2(i)*va(2,ix2l(i))
501 . - h3(i)*va(2,ix3l(i)) - h4(i)*va(2,ix4l(i))
502 vz(i) = vzi(i) - h1(i)*va(3,ix1l(i)) - h2(i)*va(3,ix2l(i))
503 . - h3(i)*va(3,ix3l(i)) - h4(i)*va(3,ix4l(i))
504 vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
505 ENDDO
506
507 DO i=1,jlt
508
509 h0 = -.25*(h1(i) - h2(i) + h3(i) - h4(i))
510 h0 =
min(h0,h2(i),h4(i))
511 h0 =
max(h0,-h1(i),-h3(i))
512 IF(ix3g(i)==ix4g(i))h0 = zero
513 h1(i) = h1(i) + h0
514 h2(i) = h2(i) - h0
515 h3(i) = h3(i) + h0
516 h4(i) = h4(i) - h0
517 ENDDO
518
519
520
521 IF(inacti==5.or.inacti==6)THEN
522
523
524
525
526
527
528
529
530
531
532
533#include "lockon.inc"
534
535 IF(igap > 0)THEN
536 DO i=1,jlt
537 is = cn_loc(i)
538 im = ce_loc(i)
539 nn = nsvg(i)
540 pplus = pene(i) + zep05*(gapv(i)-pene(i))
541 IF(nn > 0) THEN
542 IF (pplus < gap_s(is)) THEN
543 penis(2,is) =
max(penis(2,is),pplus)
544 ELSE
545 penis(2,is) =
max(penis(2,is),gap_s(is))
546 penim(2,im) =
max(penim(2,im),pplus-gap_s(is))
547 END IF
548 ELSE
549 IF (pplus <
gapfi(nin)%P(-nn))
THEN
551 ELSE
554 penim(2,im) =
max(penim(2,im),pplus-
gapfi(nin)%P(-nn))
555 END IF
556 ENDIF
557 ENDDO
558 ELSE
559 DO i=1,jlt
560 im = ce_loc(i)
561 pplus = pene(i) + zep05*(gapv(i)-pene(i))
562 penim(2,im) =
max(penim(2,im),pplus)
563 ENDDO
564 END IF
565
566
567
568
569
570
571
572
573
574
575
576
577
578#include "lockoff.inc"
579 DO i=1,jlt
580 is = cn_loc(i)
581 im = ce_loc(i)
582 nn = nsvg(i)
583 IF(nn > 0) THEN
584 pis = penis(1,is)
585 ELSE
586 pis =
penfi(nin)%P(1,-nn)
587 END IF
588 pene(i) = pene(i) - pis - penim(1,im)
589 pene(i) =
max(pene(i),zero)
590 IF (pene(i) == zero )stif(i)=zero
591 gapv(i) = gapv(i) - pis - penim(1,im)
592 END DO
593 ENDIF
594
595
596 dti = 1.e20
597
598 DO 600 i=1,jlt
599 dist=gapv(i)-pene(i)
600 rdist = half*dist /
max(em30,-vn(i))
602 600 CONTINUE
603
604
605 dtmini=ep20
606
607 IF(dti<=dtmin1(10))THEN
608 dti = 1.e20
609 DO i=1,jlt
610 dist=gapv(i)-pene(i)
611 dti2 = half*dist /
max(em30,-vn(i))
612 IF(dti2<=dtmin1(10))THEN
613#include "lockon.inc"
614 WRITE(iout,'(A,E12.4,A,I10)')
615 . ' **WARNING MINIMUM TIME STEP ',dti2,
616 . ' IN INTERFACE ',noint
617 nn = nsvg(i)
618 IF(nn>0)THEN
619 ni = itab(nn)
620 ELSE
621 ni =
itafi(nin)%P(-nn)
622 ENDIF
623#include "lockoff.inc"
624 IF(idtmin(10)==1)THEN
625#include "lockon.inc"
626 WRITE(iout,'(A,I10)') ' SECONDARY NODE : ',ni
627 WRITE(iout,'(A,4I10)')' MAIN NODES : ',
628 . itab(ix1g(i)),itab(ix2g(i)),itab(ix3g(i)),itab(ix4g(i))
629#include "lockoff.inc"
630 tstop = tt
631 ELSEIF(idtmin(10)==2)THEN
632#include "lockon.inc"
633 WRITE(iout,'(A,I10,A,I10)')' REMOVE SECONDARY NODE ',
634 . ni,' FROM INTERFACE ',noint
635 IF(nn>0) THEN
636 stfa(nsv(cn_loc(i))) = -abs(stfa(nsv(cn_loc(i))))
637 ELSE
639 ENDIF
640#include "lockoff.inc"
641 stif(i) = zero
642 newfront = -1
643 dti = dtmin1(10)
644 ELSEIF(idtmin(10)==5)THEN
645#include "lockon.inc"
646 WRITE(iout,'(A,I10)') ' SECONDARY NODE : ',ni
647 WRITE(iout,'(A,4I10)')' MAIN NODES : ',
648 . itab(ix1g(i)),itab(ix2g(i)),itab(ix3g(i)),itab(ix4g(i))
649#include "lockoff.inc"
650 mstop = 2
651 ELSEIF(idtmin(10)==6.AND.ilagm==2)THEN
652 ig=nsvg(i)
653 IF(kinet(ig)+kinet(ix1g(i))+kinet(ix2g(i))
654 . +kinet(ix3g(i))+kinet(ix4g(i))==0)THEN
655 cand_n(index(i)) = -iabs(cand_n(index(i)))
656 stif(i) = zero
657 dti2 = 1.e20
658#include "lockon.inc"
659 WRITE(iout,'(A,I10)') ' SECONDARY NODE : ',itab(nsvg(i))
660 WRITE(iout,'(A,4I10)')' MAIN NODES : ',
661 . itab(ix1g(i)),itab(ix2g(i)),itab(ix3g(i)),itab(ix4g(i))
662#include "lockoff.inc"
663 ENDIF
665 ENDIF
666 ENDIF
667 ENDDO
668 ENDIF
669 IF(dti<dt2t)THEN
670 dt2t = dti
671 neltst = noint
672 ityptst = 10
673 ENDIF
674
675 IF(impl_s>0)THEN
676 IF(imp_int7==2)THEN
677 DO i=1,jlt
678 IF(stiglo<=zero)THEN
679 stif(i) = half*stif(i)
680 ELSEIF(stif(i)/=zero)THEN
681 stif(i) = stiglo
682 ENDIF
683 fni(i)= -stif(i) * pene(i)
684 ENDDO
685 ELSE
686 DO i=1,jlt
687 fac = gapv(i)/
max( em10,( gapv(i)-pene(i) ) )
688 facm1 = 1./fac
689 IF( (gapv(i)-pene(i))/gapv(i) <prec .AND.
690 . stif(i)>0. ) THEN
691 stif(i) = 0.
692 newfront = -1
693#include "lockon.inc"
694 nn = nsvg(i)
695 IF(nn>0)THEN
696 ni = itab(nn)
697 stfa(nsv(cn_loc(i))) = -abs(stfa(nsv(cn_loc(i))))
698 ELSE
699 ni =
itafi(nin)%P(-nn)
701 ENDIF
702 WRITE(istdo,'(A,I10)')' WARNING INTERFACE ',noint
703 WRITE(istdo,'(A,I10,A)')' NODE ',ni,
704 . ' DE-ACTIVATED FROM INTERFACE'
705 WRITE(iout ,'(A,I10)')' WARNING INTERFACE ',noint
706 WRITE(iout ,'(A,I10,A)')' NODE ',ni,
707 . ' DE-ACTIVATED FROM INTERFACE'
708#include "lockoff.inc"
709 ENDIF
710 IF(stiglo<=zero)THEN
711 econtt = econtt + half*stif(i)*gapv(i)**2 *( facm1 -
712 . one -log(facm1) )
713 stif(i) = half*stif(i) * fac
714 ELSEIF(stif(i)/=zero)THEN
715 econtt = econtt + stiglo*gapv(i)**2 *( facm1 - one -
716 . log(facm1) )
717 stif(i) = stiglo * fac
718 ENDIF
719 fni(i)= -stif(i) * pene(i)
720 ENDDO
721 ENDIF
722 ELSE
723 DO 100 i=1,jlt
724 fac = gapv(i)/
max( em10,( gapv(i)-pene(i) ) )
725 facm1 = 1./fac
726 IF( (gapv(i)-pene(i))/gapv(i) <prec .AND.
727 . stif(i)>0. ) THEN
728 stif(i) = 0.
729 newfront = -1
730#include "lockon.inc"
731 nn = nsvg(i)
732 IF(nn>0)THEN
733 ni = itab(nn)
734 stfa(nsv(cn_loc(i))) = -abs(stfa(nsv(cn_loc(i))))
735 ELSE
736 ni =
itafi(nin)%P(-nn)
738 ENDIF
739 WRITE(istdo,'(A,I10)')' WARNING INTERFACE ',noint
740 WRITE(istdo,'(A,I10,A)')' NODE ',ni,
741 . ' DE-ACTIVATED FROM INTERFACE'
742 WRITE(iout ,'(A,I10)')' WARNING INTERFACE ',noint
743 WRITE(iout ,'(A,I10,A)')' NODE ',ni,
744 . ' DE-ACTIVATED FROM INTERFACE'
745#include "lockoff.inc"
746 ENDIF
747 IF(stiglo<=zero)THEN
748 econtt = econtt + half*stif(i)*gapv(i)**2 *( facm1 - one -
749 . log(facm1) )
750 stif(i) = half*stif(i) * fac
751 ELSEIF(stif(i)/=zero)THEN
752 econtt = econtt + stiglo*gapv(i)**2 *(facm1 - one - log(facm1))
753 stif(i) = stiglo * fac
754 ENDIF
755 fni(i)= -stif(i) * pene(i)
756 100 CONTINUE
757 ENDIF
758
759
760
761 IF(visc/=zero.OR.viscf/=zero)THEN
762 IF(ivis2==0)THEN
763
764
765
766 DO i=1,jlt
767 vis2(i) = two * stif(i) * msi(i)
768 IF(vn(i)<zero) vis2(i) = vis2(i) /
769 . (
max(em10,(gapv(i)-pene(i))/gapv(i)) )
770 ENDDO
771
772 visca = zep4
773 IF(kdtint==0.AND.idtmins/=2)THEN
774 DO i=1,jlt
775 fac = stif(i) /
max(em30,stif(i))
776 vis = sqrt(vis2(i))
777 ff = fac * (
778 . visc * vis +
779 . visca**2 * two* msi(i) *
max(zero,-vn(i)) /
780 .
max((gapv(i) - pene(i)),em10
781 stif(i) = stif(i) * gapv(i) /
max((gapv(i) - pene(i)),em10)
782 stif(i) = stif(i) + ff * dt1inv
783 stif(i) =
max(stif(i) ,fac*sqrt(viscf)*vis*dt1inv)
784 ffo = ff
785 ff = ff * vn(i)
786 fni(i) = fni(i) + ff
787 ENDDO
788 ELSE
789 DO i=1,jlt
790 fac = stif(i) /
max(em30,stif(i))
791 vis = sqrt(vis2(i))
792 c(i)= fac * (
793 . visc * vis +
794 . visca**2 * two * msi(i) *
max(zero,-vn(i)) /
795 .
max((gapv(i) - pene(i)),em10) )
796 stif(i) = stif(i) * gapv(i) /
max((gapv
797 kt(i)= stif(i)
798 stif(i) = stif(i) + c(i) * dt1inv
799 ff = c(i) * vn(i)
800 fni(i) = fni(i) + ff
801 cf(i) = fac*sqrt(viscf)*vis
802 stif(i) =
max(stif(i) ,cf(i
803 ENDDO
804 ENDIF
805 ELSEIF(ivis2==1)THEN
806
807
808
809 DO i=1,jlt
810 masm(i) = ms(ix1g(i))*h1(i)
811 . + ms(ix2g(i))*h2(i)
812 . + ms(ix3g(i))*h3(i)
813 . + ms(ix4g(i))*h4(i)
814 masm(i) = msi(i) * masm(i) /
max(em30,msi(i)+masm(i))
815 vis2(i) = two * stif(i) * masm(i)
816 IF(vn(i)<zero) vis2(i) = vis2(i) /
817 . (
max(em10,(gapv(i)-pene(i))/gapv(i)) )
818 ENDDO
819
820 visca = zep4
821 IF(kdtint==0.AND.idtmins/=2)THEN
822 DO i=1,jlt
823 fac = stif(i) /
max(em30,stif(i))
824 vis = sqrt(vis2(i))
825 ff = fac * (
826 . visc * vis +
827 . visca**2 * two* masm(i) *
max(zero,-vn(i)) /
828 .
max((gapv(i) - pene(i)),em10) )
829 stif(i) = stif(i) * gapv(i) /
max((gapv(i) - pene(i)),em10)
830 stif(i) = stif(i) + ff * dt1inv
831 stif(i) =
max(stif(i) ,fac*sqrt(viscf)*vis*dt1inv)
832 ffo = ff
833 ff = ff * vn(i)
834 fni(i) = fni(i) + ff
835 ENDDO
836 ELSE
837 DO i=1,jlt
838 fac = stif(i) /
max(em30,stif(i))
839 vis = sqrt(vis2(i))
840 c(i)= fac * (
841 . visc * vis +
842 . visca**2 * two * masm(i) *
max(zero,-vn(i)) /
843 .
max((gapv(i) - pene(i
844 stif(i) = stif(i) * gapv(i) /
max((gapv(i) - pene(i)),em10)
845 kt
846 stif(i) = stif(i) + c(i) * dt1inv
847 ff = c(i) * vn(i)
848 fni(i) = fni(i) + ff
849 cf(i) = fac*sqrt(viscf)*vis
850 stif(i) =
max(stif(i) ,cf(i)*dt1inv)
851 ENDDO
852 ENDIF
853 ELSEIF(ivis2==2)THEN
854
855
856
857 DO i=1,jlt
858 vis2(i) = two* stif(i) * msi(i)
859 vis2(i) = vis2(i) /
860 . (
max(em10,(gapv(i)-pene(i))/gapv(i)) )
861 ENDDO
862 visca = half
863 DO i=1,jlt
864 fac = stif(i) /
max(em30,stif(i))
865 vis = sqrt(vis2(i))
866 ff = fac * (
867 . visc * vis +
868 . visca**2 * two * msi(i) * abs(vn(i)) /
869 .
max((gapv(i) - pene(i)),em10) )
870 stif(i) = stif(i) * gapv(i) /
max((gapv(i) - pene(i)),em10)
871 stif(i) = stif(i) + two * ff * dt1inv
873 ff = ff * vn(i)
874 fni(i) = fni(i) + ff
875 ENDDO
876 ELSEIF(ivis2==3)THEN
877
878
879
880 DO i=1,jlt
881 vis2(i) = two * stif(i) * msi(i)
882 ENDDO
883
884 DO i=1,jlt
885 fac = stif(i) /
max(em30,stif(i))
886 vis = sqrt(vis2(i))
887 ff = fac * ( visc * vis
888 .
max((gapv(i) - pene(i)),em10)
889 stif(i) = stif(i) * gapv(i) /
890 .
max((gapv(i) - pene(i)),em10)
891 stif(i) = stif(i) + two* ff * dt1inv
892 stif(i) =
max(stif(i) ,fac*sqrt(viscf)*vis*dt1inv
893 ff = ff * vn(i)
894 fni(i) = fni(i) + ff
895 ENDDO
896 ELSEIF(ivis2==4)THEN
897
898
899
900 DO i=1,jlt
901 vis2(i) = two* stif(i) * msi(i)
902 vis = sqrt(vis2(i))
903 stif(i) = stif(i) * gapv(i) /
904 .
max((gapv(i) - pene(i)),em10)
905 stif(i) =
max(stif(i) ,fac*sqrt(viscf)*vis*dt1inv)
906 ENDDO
907 ELSEIF(ivis2==5)THEN
908
909
910
911
912
913 DO i=1,jlt
914 mas2 = ms(ix1g(i))*h1(i)
915 . + ms(ix2g(i))*h2(i)
916 . + ms(ix3g(i))*h3(i)
917 . + ms(ix4g(i))*h4(i)
918 vis2(i) = two* stif(i) * msi(i)
919 vis = 2. * visc * dt1inv * msi(i) * mas2 /
920 .
max(em30,msi(i)+mas2)
921 stif(i) = stif(i) * gapv(i) /
922 .
max((gapv(i) - pene(i)),em10)
923 stif(i) =
max(stif(i) ,fac*sqrt(viscf*vis2(i))*dt1inv)
924 ff = vis * vn(i)
925 econvt = econvt +
min(zero,ff-fni(i)) * vn(i) * dt1
926
927 ENDDO
928 ELSE
929 ENDIF
930 ELSE
931 DO i=1,jlt
932 vis2(i) = zero
933 stif(i) = stif(i) * gapv(i) /
934 .
max((gapv(i) - pene(i)),em10)
935 ENDDO
936 ENDIF
937
938
939
940#include "lockon.inc"
941 DO i=1,jlt
942 isign=1
943 aaa = one-pene(i)/gapv(i)
944 il = ix1l(i)
945 IF(pene(i)>zero.OR.alphak(2,il)<zero)isign=-1
946 alphak(2,il)=isign*
min(abs(alphak(2,il)),aaa)
947 il = ix2l(i)
948 IF(pene(i)>zero.OR.alphak(2,il)<zero)isign=-1
949 alphak(2,il)=isign
950 il = ix3l(i)
951 IF(pene(i)>zero.OR.alphak(2,il)<zero)isign=-1
952 alphak(2,il)=isign*
min(abs(alphak(2,il)),aaa)
953 il = ix4l(i)
954 IF(pene(i)>zero.OR.alphak(2,il)<zero)isign=-1
955 alphak(2,il)=isign*
min(abs(alphak(2,il)),aaa)
956 IF(nsvg(i)>0) THEN
957 il = nsv(cn_loc(i))
958 IF(pene(i)>zero.OR.alphak(2,il)<zero)isign=-1
959 alphak(2,il)=isign*
min(abs(alphak(2,il)),aaa)
960 ELSE
961
962 il = - nsvg(i)
963 IF(pene(i)>zero.OR.
alphakfi(nin)%P(il)<zero)isign=-1
965 ENDIF
966 ENDDO
967#include "lockoff.inc"
968
969
970
971 fsav1 = zero
972 fsav2 = zero
973 fsav3 = zero
974
975 fsav8 = zero
976 fsav9 = zero
977 fsav10= zero
978 fsav11= zero
979 DO i=1,jlt
980 fxi(i)=n1(i)*fni(i)
981 fyi(i)=n2(i)*fni(i)
982 fzi(i)=n3(i)*fni(i)
983 impx=fxi(i)*dt12
984 impy=fyi(i)*dt12
985 impz=fzi(i)*dt12
986 fsav1 =fsav1 +impx
987 fsav2 =fsav2 +impy
988 fsav3 =fsav3 +impz
989 fsav8 =fsav8 +abs(impx)
990 fsav9 =fsav9 +abs(impy)
991 fsav10=fsav10+abs(impz)
992 fsav11=fsav11+fni(i)*dt12
993 ENDDO
994#include "lockon.inc"
995 fsav(1)=fsav(1)+fsav1
996 fsav(2)=fsav(2)+fsav2
997 fsav(3)=fsav(3)+fsav3
998
999 fsav(8)=fsav(8)+fsav8
1000 fsav(9)=fsav(9)+fsav9
1001 fsav(10)=fsav(10)+fsav10
1002 fsav(11)=fsav(11)+fsav11
1003#include "lockoff.inc"
1004
1005 IF(isensint(1)/=0) THEN
1006 DO i=1,jlt
1007 fsavparit(1,1,i+nft) = fxi(i)
1008 fsavparit(1,2,i+nft) = fyi(i)
1009 fsavparit(1,3,i+nft) = fzi(i)
1010 ENDDO
1011 ENDIF
1012
1013 IF((anim_v(12)+outp_v(12)+h3d_data%N_VECT_PCONT >0.AND.
1014 . ((tt>=tanim .AND. tt<=tanim_stop).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP) .OR.
1015 . (manim>=4.AND.manim<=15).OR.h3d_data%MH3D/=0))
1016 . .OR.h3d_data%N_VECT_PCONT_MAX>0)THEN
1017#include "lockon.inc"
1018 DO i=1,jlt
1019 fncont(1,ix1g(i)) =fncont(1,ix1g(i)) + fxi(i)*h1(i)
1020 fncont(2,ix1g(i)) =fncont(2,ix1g(i)) + fyi(i)*h1(i)
1021 fncont(3,ix1g(i)) =fncont(3,ix1g(i)) + fzi(i)*h1(i)
1022 fncont(1,ix2g(i)) =fncont(1,ix2g(i)) + fxi(i)*h2(i)
1023 fncont(2,ix2g(i)) =fncont(2,ix2g(i)) + fyi(i)*h2(i)
1024 fncont(3,ix2g(i)) =fncont(3,ix2g(i)) + fzi(i)*h2(i)
1025 fncont(1,ix3g(i)) =fncont(1,ix3g(i)) + fxi(i)*h3(i
1026 fncont(2,ix3g(i)) =fncont(2,ix3g(i)) + fyi(i)*h3(i)
1027 fncont(3,ix3g(i)) =fncont(3,ix3g(i)) + fzi(i)*h3(i)
1028 fncont(1,ix4g(i)) =fncont(1,ix4g(i)) + fxi(i)*h4(i)
1029 fncont(2,ix4g(i)) =fncont(2,ix4g(i)) + fyi(i)*h4(i)
1030 fncont(3,ix4g(i)) =fncont(3,ix4g(i)) + fzi(i)*h4(i)
1031 jg = nsvg(i)
1032 IF(jg>0) THEN
1033
1034 fncont(1,jg)=fncont(1,jg)- fxi(i)
1035 fncont(2,jg)=fncont(2,jg)- fyi(i)
1036 fncont(3,jg)=fncont(3,jg)- fzi(i)
1037 ELSE
1038 jg = -jg
1042 ENDIF
1043 ENDDO
1044#include "lockoff.inc"
1045 ENDIF
1046
1047
1048
1049 IF(nisub/=0)THEN
1050 DO jsub=1,nisub
1051 DO j=1,15
1052 fsavsub1(j,jsub)=zero
1053 END DO
1054 ENDDO
1055 DO i=1,jlt
1056 nn = nsvg(i)
1057 IF(nn>0)THEN
1058 in=cn_loc(i)
1059 ie=ce_loc(i)
1060 jj =addsubs(in)
1061 kk =addsubm(ie)
1062 DO WHILE(jj<addsubs(in+1))
1063 jsub=lisubs(jj)
1064 DO WHILE(kk<addsubm(ie+1))
1065 ksub=lisubm(kk)
1066 IF(ksub==jsub)THEN
1067 impx=fxi(i)*dt12
1068 impy=fyi(i)*dt12
1069 impz=fzi(i)*dt12
1070
1071 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1072 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1073 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1074
1075 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1076 fsavsub1(9,jsub) =fsavsub1(9,jsub)
1077 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1078
1079 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1080 kk=kk+1
1081 GO TO 250
1082 ELSE IF(ksub<jsub)THEN
1083 kk=kk+1
1084 ELSE
1085 GO TO 250
1086 END IF
1087 END DO
1088 250 CONTINUE
1089 jj=jj+1
1090 END DO
1091 END IF
1092 END DO
1093
1094 IF(nspmd>1) THEN
1095
1096 DO i=1,jlt
1097 nn = nsvg(i)
1098 IF(nn<0)THEN
1099 nn = -nn
1100 ie=ce_loc(i)
1102 kk =addsubm(ie)
1105 DO WHILE(kk<addsubm(ie+1))
1106 ksub=lisubm(kk)
1107 IF(ksub==jsub)THEN
1108 impx=fxi(i)*dt12
1109 impy=fyi(i)*dt12
1110 impz=fzi(i)*dt12
1111
1112 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1113 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1114 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1115
1116 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1117 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1118 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs
1119
1120 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1121 kk=kk+1
1122 GO TO 150
1123 ELSE IF(ksub<jsub)THEN
1124 kk=kk+1
1125 ELSE
1126 GO TO 150
1127 END IF
1128 END DO
1129 150 CONTINUE
1130 jj=jj+1
1131 END DO
1132 END IF
1133
1134 END DO
1135
1136 END IF
1137 END IF
1138
1139
1140
1141
1142 IF (mfrot==0) THEN
1143
1144 DO i=1,jlt
1145 xmu(i) = fric
1146 ENDDO
1147 ELSEIF (mfrot==1) THEN
1148
1149 DO i=1,jlt
1150 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1151 v2 = (vx(i) - n1(i)*aa)**2
1152 . + (vy(i) - n2(i)*aa)**2
1153 . + (vz(i) - n3(i)*aa)**2
1154 vv = sqrt(
max(em30,v2))
1155 ax1 = x3(i) - x1(i)
1156 ay1 = y3(i) - y1(i)
1157 az1 = z3(i) - z1(i)
1158 ax2 = x4(i) - x2(i)
1159 ay2 = y4(i) - y2(i)
1160 az2 = z4(i) - z2(i)
1161 ax = ay1*az2 - az1*ay2
1162 ay = az1*ax2 - ax1*az2
1163 az = ax1*ay2 - ay1*ax2
1164 area = half*sqrt(ax*ax+ay*ay+az*az)
1166 xmu(i) = fric + (frot_p(1) + frot_p(4)*p ) * p
1167 . +(frot_p(2) + frot_p(3)*p) * vv + frot_p(5)*v2
1168 ENDDO
1169 ELSEIF(mfrot==2)THEN
1170
1171 DO i=1,jlt
1172 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1173 v2 = (vx(i) - n1(i)*aa)**2
1174 . + (vy(i) - n2(i)*aa)**2
1175 . + (vz(i) - n3(i)*aa)**2
1176 vv = sqrt(
max(em30,v2))
1177 ax1 = x3(i) - x1(i)
1178 ay1 = y3(i) - y1(i)
1179 az1 = z3(i) - z1(i)
1180 ax2 = x4(i) - x2(i)
1181 ay2 = y4(i) - y2(i)
1182 az2 = z4(i) - z2(i)
1183 ax = ay1*az2 - az1*ay2
1184 ay = az1*ax2 - ax1*az2
1185 az = ax1*ay2 - ay1*ax2
1186 area = half*sqrt(ax*ax+ay*ay+az*az)
1188 xmu(i) = fric
1189 . + frot_p(1)*exp(frot_p(2)*vv)*p*p
1190 . + frot_p(3)*exp(frot_p(4)*vv)*p
1191 . + frot_p(5)*exp(frot_p(6)*vv)
1192 ENDDO
1193 ELSEIF (mfrot==3) THEN
1194
1195 DO i=1,jlt
1196 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1197 v2 = (vx(i) - n1(i)*aa)**2
1198 . + (vy(i) - n2(i)*aa)**2
1199 . + (vz(i) - n3(i)*aa)**2
1200 vv = sqrt(
max(em30,v2))
1201 IF(vv>=0.AND.vv<=frot_p(5)) THEN
1202 dmu = frot_p(3)-frot_p(1)
1203 vv1
1204 xmu(i) = frot_p(1)+ dmu*vv1*(two-vv1)
1205 ELSEIF(vv>frot_p(5).AND.vv<frot_p(6)) THEN
1206 dmu = frot_p(4)-frot_p(3)
1207 vv1 = (vv - frot_p(5))/(frot_p(6)-frot_p(5))
1208 xmu(i) = frot_p(3)+ dmu * (three-two*vv1)*vv1**2
1209 ELSE
1210 dmu = frot_p(2)-frot_p(4)
1211 vv2 = (vv - frot_p(6))**2
1212 xmu(i) = frot_p(2) - dmu / (one + dmu*vv2)
1213 ENDIF
1214 ENDDO
1215 ELSEIF(mfrot==4)THEN
1216
1217 DO i=1,jlt
1218 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1219 v2 = (vx(i) - n1(i)*aa)**2
1220 . + (vy(i) - n2(i)*aa)**2
1221 . + (vz(i) - n3(i)*aa)**2
1222 vv = sqrt(
max(em30,v2))
1223 xmu(i) = frot_p(1)
1224 . + (fric-frot_p(1))*exp(-frot_p(2)*vv)
1225 xmu(i) =
max(xmu(i),em30)
1226 ENDDO
1227 ENDIF
1228
1229
1230
1231 fsav4 = zero
1232 fsav5 = zero
1233 fsav6 = zero
1234
1235 fsav12= zero
1236 fsav13= zero
1237 fsav14= zero
1238 fsav15= zero
1239
1240 IF (ifq>=10) THEN
1241
1242
1243
1244 IF (ifq==13) THEN
1246 ELSE
1248 ENDIF
1249 DO i=1,jlt
1250 fx = stif0(i)*vx(i)*dt12
1251 fy = stif0(i)*vy(i)*dt12
1252 fz = stif0(i)*vz(i)*dt12
1253
1254 fx = cand_fx(index(i)) +
alpha*fx
1255 fy = cand_fy(index(i)) +
alpha*fy
1256 fz = cand_fz(index(i)) +
alpha*fz
1257
1258 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
1259 fx = fx - ftn*n1(i)
1260 fy = fy - ftn*n2(i)
1261 fz = fz - ftn*n3(i)
1262 ft = fx*fx + fy*fy + fz*fz
1264
1265 fn = fxi(i)**2+fyi(i)**2+fzi(i)
1266
1267 beta =
min(one,xmu(i)*sqrt(fn/ft))
1268
1269 fxt(i) = fx * beta
1270 fyt(i) = fy * beta
1271 fzt(i) = fz * beta
1272
1273 cand_fx(index(i)) = fxt(i)
1274 cand_fy(index(i)) = fyt(i)
1275 cand_fz(index(i)) = fzt(i)
1276 ifpen(index(i)) = 1
1277
1278
1279 fxi(i)=fxi(i) + fxt(i)
1280 fyi(i)=fyi(i) + fyt(i)
1281 fzi(i)=fzi(i) + fzt(i)
1282 econvt = econvt
1283 . + dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
1284 ENDDO
1285
1286
1287
1288 ELSEIF (ifq>0) THEN
1289
1290 IF (ifq==3) THEN
1292 ELSE
1294 ENDIF
1296 DO i=1,jlt
1297 vnx = n1(i)*vn(i)
1298 vny = n2(i)*vn(i)
1299 vnz = n3(i)*vn(i)
1300 vx(i) = vx(i) - vnx
1301 vy(i) = vy(i) - vny
1302 vz(i) = vz
1303 v2 = vx(i)**2 + vy(i)**2 + vz(i)**2
1304 vis2(i) = viscf * vis2(i)
1305 fm2 = (xmu(i)*fni(i))**2
1306 f2 = vis2(i) * v2
1307 a2 =
min(f2,fm2) /
max(em30,f2)
1308 aa = sqrt(a2 * vis2(i))
1309 fx = aa * vx(i)
1310 fy = aa * vy(i)
1311 fz = aa * vz(i)
1312
1313 fxt(i) =
alpha*fx + alphi*cand_fx(index(i))
1314 fyt(i) =
alpha*fy + alphi*cand_fy(index(i))
1315 fzt(i) =
alpha*fz + alphi*cand_fz(index(i))
1316 cand_fx(index(i)) = fxt(i)
1317 cand_fy(index(i)) = fyt(i)
1318 cand_fz(index(i)) = fzt(i)
1319 ifpen(index(i)) = 1
1320
1321 fxi(i) = fxi(i) + fxt(i)
1322 fyi(i) = fyi(i) + fyt(i)
1323 fzi(i) = fzi(i) + fzt(i)
1324 econvt = econvt
1325 . + dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
1326 ENDDO
1327 ELSE
1328
1329
1330
1331 DO i=1,jlt
1332 vnx = n1(i)*vn(i)
1333 vny = n2(i)*vn(i)
1334 vnz = n3(i)*vn(i)
1335 vx(i) = vx(i) - vnx
1336 vy(i) = vy(i) - vny
1337 vz(i) = vz(i) - vnz
1338 v2 = vx(i)**2 + vy(i)**2 + vz(i)**2
1339 vis2(i) = viscf * vis2(i)
1340 fm2 = (xmu(i)*fni(i))**2
1341 f2 = vis2(i) * v2
1342 a2 =
min(f2,fm2) /
max(em30,f2)
1343 aa = sqrt(a2 * vis2(i))
1344 fxt(i) = aa * vx(i)
1345 fyt(i) = aa * vy(i)
1346 fzt(i) = aa * vz(i)
1347
1348 fxi(i)=fxi(i) + fxt(i)
1349 fyi(i)=fyi(i) + fyt(i)
1350 fzi(i)=fzi(i) + fzt(i)
1351 econvt = econvt + aa * v2 * dt1
1352 ENDDO
1353 ENDIF
1354
1355 IF((anim_v(12)+outp_v(12)+h3d_data%N_VECT_PCONT>0.AND.
1356 . ((tt>=tanim .AND. tt<=tanim_stop).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
1357 . (manim>=4.AND.manim<=15).OR.h3d_data%MH3D/=0))
1358 . .OR.h3d_data%N_VECT_PCONT_MAX>0)THEN
1359#include "lockon.inc"
1360 DO i=1,jlt
1361 ftcont(1,ix1g(i)) =ftcont(1,ix1g(i)) + fxt(i)*h1(i)
1362 ftcont(2,ix1g(i)) =ftcont(2,ix1g(i)) + fyt(i)*h1(i)
1363 ftcont(3,ix1g(i)) =ftcont(3,ix1g(i)) + fzt(i)*h1(i)
1364 ftcont(1,ix2g(i)) =ftcont(1,ix2g(i)) + fxt(i)*h2(i)
1365 ftcont(2,ix2g(i)) =ftcont(2,ix2g(i)) + fyt(i)*h2(i)
1366 ftcont(3,ix2g(i)) =ftcont(3,ix2g(i)) + fzt(i)*h2(i)
1367 ftcont(1,ix3g
1368 ftcont(2,ix3g(i)) =ftcont(2,ix3g(i)) + fyt(i)*h3(i)
1369 ftcont(3,ix3g(i)) =ftcont(3,ix3g(i)) + fzt(i)*h3(i)
1370 ftcont(1,ix4g(i)) =ftcont(1,ix4g(i)) + fxt(i)*h4(i)
1371 ftcont(2,ix4g(i)) =ftcont(2,ix4g(i)) + fyt(i)*h4(i)
1372 ftcont(3,ix4g(i)) =ftcont(3,ix4g(i)) + fzt(i)*h4(i)
1373 jg = nsvg(i)
1374 IF(jg>0) THEN
1375
1376 ftcont(1,jg)=ftcont(1,jg)- fxt(i)
1377 ftcont(2,jg)=ftcont(2,jg)- fyt(i)
1378 ftcont(3,jg)=ftcont(3,jg)- fzt(i)
1379 ELSE
1380 jg = -jg
1384 ENDIF
1385 ENDDO
1386#include "lockoff.inc"
1387 ENDIF
1388
1389
1390 DO i=1,jlt
1391 impx=fxt(i)*dt12
1392 impy=fyt(i)*dt12
1393 impz=fzt(i)*dt12
1394 fsav4 =fsav4 +impx
1395 fsav5 =fsav5 +impy
1396 fsav6 =fsav6 +impz
1397 impx=fxi(i)*dt12
1398 impy=fyi(i)*dt12
1399 impz=fzi(i)*dt12
1400 fsav12=fsav12+abs(impx)
1401 fsav13=fsav13+abs(impy)
1402 fsav14=fsav14+abs(impz)
1403 fsav15=fsav15+sqrt(impx*impx
1404 ENDDO
1405#include "lockon.inc"
1406 fsav(4) = fsav(4) + fsav4
1407 fsav(5) = fsav(5) + fsav5
1408 fsav(6) = fsav(6) + fsav6
1409
1410 fsav(12) = fsav(12) + fsav12
1411 fsav(13) = fsav(13) + fsav13
1412 fsav(14) = fsav(14) + fsav14
1413 fsav(15) = fsav(15) + fsav15
1414#include "lockoff.inc"
1415
1416 IF(isensint(1)/=0) THEN
1417 DO i=1,jlt
1418 fsavparit(1,4,i+nft) = fxt(i)
1419 fsavparit(1,5,i+nft) = fyt(i)
1420 fsavparit(1,6,i+nft) = fzt(i)
1421 ENDDO
1422 ENDIF
1423
1424
1425
1426
1427 IF(nisub/=0)THEN
1428 DO i=1,jlt
1429 nn = nsvg(i)
1430 IF(nn>0)THEN
1431 in=cn_loc(i)
1432 ie=ce_loc(i)
1433 jj =addsubs(in)
1434 kk =addsubm(ie)
1435 DO WHILE(jj<addsubs(in+1))
1436 jsub=lisubs(jj)
1437 DO WHILE(kk<addsubm(ie+1))
1438 ksub=lisubm(kk)
1439 IF(ksub==jsub)THEN
1440 impx=fxt(i)*dt12
1441 impy=fyt(i)*dt12
1442 impz=fzt(i)*dt12
1443
1444 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
1445 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
1446 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
1447
1448 impx=fxi(i)*dt12
1449 impy=fyi(i)*dt12
1450 impz=fzi(i)*dt12
1451 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
1452 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
1453 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
1454
1455 fsavsub1(15,jsub)= fsavsub1(15,jsub)
1456 . +sqrt(impx*impx+impy*impy+impz*impz)
1457 kk=kk+1
1458 GO TO 200
1459 ELSE IF(ksub<jsub)THEN
1460 kk=kk+1
1461 ELSE
1462 GO TO 200
1463 END IF
1464 END DO
1465 200 CONTINUE
1466 jj=jj+1
1467 END DO
1468 END IF
1469 END DO
1470
1471 IF(nspmd>1) THEN
1472
1473 DO i=1,jlt
1474 nn = nsvg(i)
1475 IF(nn<0)THEN
1476
1477 nn = -nn
1478 ie=ce_loc(i)
1480 kk =addsubm(ie)
1483 DO WHILE(kk<addsubm(ie+1))
1484 ksub=lisubm(kk)
1485 IF(ksub==jsub)THEN
1486 impx=fxt(i)*dt12
1487 impy=fyt(i)*dt12
1488 impz=fzt(i)*dt12
1489
1490 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
1491 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
1492 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
1493
1494 impx=fxi(i)*dt12
1495 impy=fyi(i)*dt12
1496 impz=fzi(i)*dt12
1497 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
1498 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
1499 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
1500
1501 fsavsub1(15,jsub)= fsavsub1(15,jsub)
1502 . +sqrt(impx*impx+impy*impy+impz*impz)
1503 kk=kk+1
1504 GO TO 300
1505 ELSE IF(ksub<jsub)THEN
1506 kk=kk+1
1507 ELSE
1508 GO TO 300
1509 END IF
1510 END DO
1511 300 CONTINUE
1512 jj=jj+1
1513 END DO
1514 END IF
1515
1516 END DO
1517
1518 END IF
1519#include "lockon.inc"
1520 DO jsub=1,nisub
1521 nsub=lisub(jsub)
1522 DO j=1,15
1523 fsavsub(j,nsub)=fsavsub(j,nsub)+fsavsub1(j,jsub)
1524 END DO
1525 END DO
1526#include "lockoff.inc"
1527 END IF
1528
1529#include "lockon.inc"
1530 econtv = econtv + econvt
1531 econt = econt + econtt
1532#include "lockoff.inc"
1533
1534 IF(kdtint==1)THEN
1535 IF( (visc/=zero.OR.viscf/=zero)
1536 . .AND.(ivis2==0.OR.ivis2==1))THEN
1537 DO i=1,jlt
1538
1539 IF(msi(i)==zero)THEN
1540 ks(i) =zero
1541 cs(i) =zero
1542 stv(i)=zero
1543 ELSE
1544 cx = four*c(i)*c(i)
1545 cy = eight*msi(i)*kt(i)
1546 aux = sqrt(cx+cy)+two*c(i)
1547 stv(i)= kt(i)*aux*aux/
max(cy,em30)
1548 aux = two*cf(i)*cf(i)/
max(msi(i),em20)
1549 IF(aux>stv(i))THEN
1550 ks(i) =zero
1551 cs(i) =cf(i)
1552 stv(i)=aux
1553 ELSE
1554 ks(i)= kt(i)
1555 cs(i) =c(i)
1556 ENDIF
1557 ENDIF
1558
1559 j1=ix1g(i)
1560 IF(ms(j1)==zero)THEN
1561 k1(i) =zero
1562 c1
1563 st1(i)=zero
1564 ELSE
1565 k1(i)=kt(i)*abs(h1(i))
1566 c1(i)=c(i)*abs(h1(i))
1567 cx =four*c1(i)*c1(i)
1568 cy =eight*ms(j1)*k1(i)
1569 aux = sqrt(cx+cy)+two*c1(i)
1570 st1(i)= k1(i)*aux*aux/
max(cy,em30)
1571 cfi = cf(i)*abs(h1(i))
1572 aux = two*cfi*cfi/
max(ms(j1),em20)
1573 IF(aux>st1(i))THEN
1574 k1(i) =zero
1575 c1(i) =cfi
1576 st1(i)=aux
1577 ENDIF
1578 ENDIF
1579
1580 j1=ix2g(i)
1581 IF(ms(j1)==zero)THEN
1582 k2(i) =zero
1583 c2(i) =zero
1584 st2(i)=zero
1585 ELSE
1586 k2(i)=kt(i)*abs(h2(i))
1587 c2(i)=c(i)*abs(h2(i))
1588 cx =four*c2(i)*c2(i)
1589 cy =eight*ms(j1)*k2(i)
1590 aux = sqrt(cx+cy)+two*c2(i)
1591 st2(i)= k2(i)*aux*aux/
max(cy,em30)
1592 cfi = cf(i)*abs(h2(i))
1593 aux = two*cfi*cfi/
max(ms(j1),em20)
1594 IF(aux>st2(i))THEN
1595 k2(i) =zero
1596 c2(i) =cfi
1597 st2(i)=aux
1598 ENDIF
1599 ENDIF
1600
1601 j1=ix3g(i)
1602 IF(ms(j1)==zero)THEN
1603 k3(i) =zero
1604 c3(i) =zero
1605 st3(i)=zero
1606 ELSE
1607 k3(i)=kt(i)*abs(h3(i))
1608 c3(i)=c(i)*abs(h3(i))
1609 cx =four*c3(i)*c3(i)
1610 cy =eight*ms(j1)*k3(i)
1611 aux = sqrt(cx+cy)+two*c3(i)
1612 st3(i)= k3(i)*aux*aux/
max(cy,em30)
1613 cfi = cf(i)*abs(h3(i))
1614 aux = two*cfi*cfi/
max(ms(j1),em20)
1615 IF(aux>st3(i))THEN
1616 k3(i) =zero
1617 c3(i) =cfi
1618 st3(i)=aux
1619 ENDIF
1620 ENDIF
1621
1622 j1=ix4g(i)
1623 IF(ms(j1)==zero)THEN
1624 k4(i) =zero
1625 c4(i) =zero
1626 st4(i)=zero
1627 ELSE
1628 k4(i)=kt(i)*abs
1629 c4(i)=c(i)*abs(h4(i))
1630 cx =four*c4(i)*c4(i)
1631 cy =eight*ms(j1)*k4(i)
1632 aux = sqrt(cx+cy)+two*c4(i)
1633 st4(i)= k4(i)*aux*aux/
max(cy,em30)
1634 cfi = cf(i)*abs(h4(i))
1635 aux = two*cfi*cfi/
max(ms(j1),em20)
1636 IF(aux>st4(i))THEN
1637 k4(i) =zero
1638 c4(i) =cfi
1639 st4(i)=aux
1640 ENDIF
1641 ENDIF
1642 ENDDO
1643
1644 ELSE
1645 DO i=1,jlt
1646 ks(i) =stif(i)
1647 cs(i) =zero
1648 stv(i)=ks(i)
1649 k1(i) =stif(i)*abs(h1(i))
1650 c1(i) =zero
1651 st1(i)=k1(i)
1652 k2(i) =stif(i)*abs(h2(i))
1653 c2(i) =zero
1654 st2(i)=k2(i)
1655 k3(i) =stif(i)*abs(h3(i))
1656 c3(i) =zero
1657 st3(i)=k3(i)
1658 k4(i) =stif(i)*abs(h4(i))
1659 c4(i) =zero
1660 st4(i)=k4(i)
1661 ENDDO
1662 ENDIF
1663 ENDIF
1664
1665 IF(idtmin(10)==1.OR.idtmin(10)==2.OR.
1666 . idtmin(10)==5.OR.idtmin(10)==6)THEN
1667
1668 dtmi0 = ep20
1669 IF(kdtint==0)THEN
1670 DO i=1,jlt
1671 dtmi(i) = ep20
1672 mas2 = two * msi(i)
1673 IF(mas2>zero.AND.stif(i)>zero.AND.
1674 . irb(kini(i))==0.AND.irb2(kini(i))==0)THEN
1675 dtmi(i) =
min(dtmi(i),dtfac1(10)*sqrt(mas2/stif(i)))
1676 ENDIF
1677 mas2 = two
1678 IF(mas2>zero.AND.h1(i)*stif(i)>zero.AND.
1679 . irb(kinet(ix1g(i)))==0.AND.irb2(kinet(ix1g(i)))==0)THEN
1680 dtmi(i) =
min(dtmi(i),dtfac1(10)*sqrt(mas2/(h1(i)*stif(i))))
1681 ENDIF
1682 mas2 = two * ms(ix2g(i))
1683 IF(mas2>zero.AND.h2(i)*stif(i)>zero.AND.
1684 . irb(kinet(ix2g(i)))==0.AND.irb2(kinet(ix2g(i)))==0)THEN
1686 ENDIF
1687 mas2 = two* ms(ix3g(i))
1688 IF(mas2>zero.AND.h3(i)*stif(i)>zero.AND.
1689 . irb(kinet(ix3g(i)))==0.AND.irb2(kinet(ix3g(i)))==0)THEN
1690 dtmi(i) =
min(dtmi(i),dtfac1(10)*sqrt(mas2/(h3(i)*stif(i))))
1691 ENDIF
1692 mas2 = two * ms(ix4g(i))
1693 IF(mas2>zero.AND.h4(i)*stif(i)>zero.AND.
1694 . irb(kinet(ix4g(i)))==0.AND.irb2(kinet(ix4g(i)))==0)THEN
1695 dtmi(i) =
min(dtmi(i),dtfac1(10)*sqrt(mas2/(h4(i)*stif(i))))
1696 ENDIF
1697 dtmi0 =
min(dtmi0,dtmi(i))
1698 ENDDO
1699
1700 ELSE
1701 DO i=1,jlt
1702 dtmi(i) = ep20
1703 mas2 = two * msi(i)
1704 mas2 = two * msi(i)
1705 IF(mas2>zero.AND.stv(i)>zero.AND.
1706 . irb(kini(i))==0.AND.irb2(kini(i))==0)THEN
1708 ENDIF
1709 mas2 = two * ms(ix1g(i))
1710 IF(mas2>zero.AND.st1(i)>zero.AND.
1711 . irb(kinet(ix1g(i)))==0.AND.irb2(kinet(ix1g(i)))==0)THEN
1712 dtmi(i) =
min(dtmi(i),dtfac1(10)*sqrt(mas2/(st1(i))))
1713 ENDIF
1714 mas2 = two * ms(ix2g(i))
1715 IF(mas2>zero.AND.st2(i)>zero.AND.
1716 . irb(kinet(ix2g(i)))==0.AND.irb2(kinet(ix2g(i)))==0)THEN
1717 dtmi(i) =
min(dtmi(i),dtfac1(10)*sqrt(mas2/(st2(i))))
1718 ENDIF
1719 mas2 = two * ms(ix3g(i))
1720 IF(mas2>zero.AND.st3(i)>zero.AND.
1721 . irb(kinet(ix3g(i)))==0.AND.irb2(kinet(ix3g(i)))==0)THEN
1722 dtmi(i) =
min(dtmi(i),dtfac1(10)*sqrt(mas2/(st3(i))))
1723 ENDIF
1724 mas2 = two * ms(ix4g(i))
1725 IF(mas2>zero.AND.st4(i)>zero.AND.
1726 . irb(kinet(ix4g(i)))==0.AND.irb2(kinet(ix4g(i)))==0)THEN
1727 dtmi(i) =
min(dtmi(i),dtfac1(10)*sqrt(mas2/(st4(i))))
1728 ENDIF
1729 dtmi0 =
min(dtmi0,dtmi(i))
1730 ENDDO
1731 ENDIF
1732 IF(dtmi0<=dtmin1(10))THEN
1733 DO i=1,jlt
1734 IF(dtmi(i)<=dtmin1(10))THEN
1735 jg = nsvg(i)
1736 IF(jg>0)THEN
1737 ni = itab(jg)
1738 ELSE
1739 ni =
itafi(nin)%P(-jg)
1740 ENDIF
1741 IF(idtmin(10)==1)THEN
1742#include "lockon.inc"
1743 WRITE(iout,'(A,E12.4,A,I10)')
1744 . ' **WARNING MINIMUM TIME STEP ',dtmi(i),
1745 . ' IN INTERFACE ',noint
1746 WRITE(iout,'(A,I10)') ' SECONDARY NODE : ',ni
1747 WRITE(iout,'(A,4I10)')' MAIN NODES : ',
1748 . itab(ix1g(i)),itab(ix2g(i)),itab(ix3g(i)),itab(ix4g(i))
1749#include "lockoff.inc"
1750 tstop = tt
1751 ELSEIF(idtmin(10)==2)THEN
1752#include "lockon.inc"
1753 WRITE(iout,'(A,E12.4,A,I10)')
1754 . ' **WARNING MINIMUM TIME STEP ',dtmi(i),
1755 . ' IN INTERFACE ',noint
1756 WRITE(iout,'(A,I10,A,I10)')' delete secondary node ',
1757 . NI,' from INTERFACE ',NOINT
1758 WRITE(IOUT,'(a,4i10)
')' main nodes :
',
1759 . ITAB(IX1G(I)),ITAB(IX2G(I)),ITAB(IX3G(I)),ITAB(IX4G(I))
1760 IF(JG>0) THEN
1761 STFA(NSV(CN_LOC(I))) = -ABS(STFA(NSV(CN_LOC(I))))
1762 ELSE
1763 STIFI(NIN)%P(-JG) = -ABS(STIFI(NIN)%P(-JG))
1764 ENDIF
1765#include "lockoff.inc"
1766 NEWFRONT = -1
1767 ELSEIF(IDTMIN(10)==5)THEN
1768#include "lockon.inc"
1769 WRITE(IOUT,'(a,e12.4,a,i10)')
1770 . ' **warning minimum time step ',DTMI(I),
1771 . ' in INTERFACE ',NOINT
1772 WRITE(IOUT,'(a,i10)') ' secondary node : ',NI
1773 WRITE(IOUT,'(a,4i10)
')' main nodes :
',
1774 . ITAB(IX1G(I)),ITAB(IX2G(I)),ITAB(IX3G(I)),ITAB(IX4G(I))
1775#include "lockoff.inc"
1776 MSTOP = 2
1777.AND. ELSEIF(IDTMIN(10)==6ILAGM==2)THEN
1778 IF(KINET(JG)+KINET(IX1G(I))+KINET(IX2G(I))
1779 . +KINET(IX3G(I))+KINET(IX4G(I))==0 )THEN
1780 CAND_N(INDEX(I)) = -IABS(CAND_N(INDEX(I)))
1781 STIF(I) = 0.
1782 FXI(I) = 0.
1783 FYI(I) = 0.
1784 FZI(I) = 0.
1785 ENDIF
1786 ENDIF
1787 ENDIF
1788 ENDDO
1789 ENDIF
1790 ENDIF
1791
1792
1793
1794 DO I=1,JLT
1795 FX1(I)=FXI(I)*H1(I)
1796 FY1(I)=FYI(I)*H1(I)
1797 FZ1(I)=FZI(I)*H1(I)
1798
1799 FX2(I)=FXI(I)*H2(I)
1800 FY2(I)=FYI(I)*H2(I)
1801 FZ2(I)=FZI(I)*H2(I)
1802
1803 FX3(I)=FXI(I)*H3(I)
1804 FY3(I)=FYI(I)*H3(I)
1805 FZ3(I)=FZI(I)*H3(I)
1806
1807 FX4(I)=FXI(I)*H4(I)
1808 FY4(I)=FYI(I)*H4(I)
1809 FZ4(I)=FZI(I)*H4(I)
1810 ENDDO
1811
1812
1813
1814
1815 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FXI, FX6)
1816 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FYI, FY6)
1817 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FZI, FZ6)
1818#include "lockon.inc"
1819
1820 DO I = 1,JLT
1821 IG = NSVG(I)
1822 IF(IG > 0)THEN
1823 IL = NSV(CN_LOC(I))
1824 DO K = 1,6
1825 DAANC6(1,K,IL) = DAANC6(1,K,IL) - FX6(K,I)
1826 DAANC6(2,K,IL) = DAANC6(2,K,IL) - FY6(K,I)
1827 DAANC6(3,K,IL) = DAANC6(3,K,IL) - FZ6(K,I)
1828 ENDDO
1829 ELSE
1830
1831
1832
1833 IL = - IG
1834 DO K = 1,6
1835 DAANC6FI(NIN)%P(1,K,IL) = DAANC6FI(NIN)%P(1,K,IL) - FX6(K,I)
1836 DAANC6FI(NIN)%P(2,K,IL) = DAANC6FI(NIN)%P(2,K,IL) - FY6(K,I)
1837 DAANC6FI(NIN)%P(3,K,IL) = DAANC6FI(NIN)%P(3,K,IL) - FZ6(K,I)
1838 ENDDO
1839 ENDIF
1840 ENDDO
1841
1842
1843
1844 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FX1, FX6)
1845 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FY1, FY6)
1846 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FZ1, FZ6)
1847 DO I = 1,JLT
1848 IL = IX1L(I)
1849 DO K = 1,6
1850 DAANC6(1,K,IL) = DAANC6(1,K,IL) + FX6(K,I)
1851 DAANC6(2,K,IL) = DAANC6(2,K,IL) + FY6(K,I)
1852 DAANC6(3,K,IL) = DAANC6(3,K,IL) + FZ6(K,I)
1853 ENDDO
1854 ENDDO
1855
1856 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FX2, FX6)
1857 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FY2, FY6)
1858 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FZ2, FZ6)
1859 DO I = 1,JLT
1860 IL = IX2L(I)
1861 DO K = 1,6
1862 DAANC6(1,K,IL) = DAANC6(1,K,IL) + FX6(K,I)
1863 DAANC6(2,K,IL) = DAANC6(2,K,IL) + FY6(K,I)
1864 DAANC6(3,K,IL) = DAANC6(3,K,IL) + FZ6(K,I)
1865 ENDDO
1866 ENDDO
1867
1868 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FX3, FX6)
1869 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FY3, FY6)
1870 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FZ3, FZ6)
1871 DO I = 1,JLT
1872 IL = IX3L(I)
1873 DO K = 1,6
1874 DAANC6(1,K,IL) = DAANC6(1,K,IL) + FX6(K,I)
1875 DAANC6(2,K,IL) = DAANC6(2,K,IL) + FY6(K,I)
1876 DAANC6(3,K,IL) = DAANC6(3,K,IL) + FZ6(K,I)
1877 ENDDO
1878 ENDDO
1879
1880 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FX4, FX6)
1881 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FY4, FY6)
1882 CALL FOAT_TO_6_FLOAT(1 ,JLT ,FZ4, FZ6)
1883 DO I = 1,JLT
1884 IL = IX4L(I)
1885 DO K = 1,6
1886 DAANC6(1,K,IL) = DAANC6(1,K,IL) + FX6(K,I)
1887 DAANC6(2,K,IL) = DAANC6(2,K,IL) + FY6(K,I)
1888 DAANC6(3,K,IL) = DAANC6(3,K,IL) + FZ6(K,I)
1889 ENDDO
1890 ENDDO
1891#include "lockoff.inc"
1892
1893
1894
1895
1896
1897 DO I = 1,JLT
1898 IF(GAPV(I) > GAPR(I))THEN
1899 IG = NSVG(I)
1900 IF(IG > 0)THEN
1901 IL = NSV(CN_LOC(I))
1902 XSA = N1(I)*(DXANC(1,IL)-H1(I)*DXANC(1,IX1L(I))
1903 . -H2(I)*DXANC(1,IX2L(I))
1904 . -H3(I)*DXANC(1,IX3L(I))
1905 . -H4(I)*DXANC(1,IX4L(I)))
1906 . + N2(I)*(DXANC(2,IL)-H1(I)*DXANC(2,IX1L(I))
1907 . -H2(I)*DXANC(2,IX2L(I))
1908 . -H3(I)*DXANC(2,IX3L(I))
1909 . -H4(I)*DXANC(2,IX4L(I)))
1910 . + N3(I)*(DXANC(3,IL)-H1(I)*DXANC(3,IX1L(I))
1911 . -H2(I)*DXANC(3,IX2L(I))
1912 . -H3(I)*DXANC(3,IX3L(I))
1913 . -H4(I)*DXANC(3,IX4L(I)))
1914 ELSE
1915
1916
1917
1918
1919
1920 IL = - IG
1921 XSA = N1(I)*(DXANCFI(NIN)%P(1,IL)-H1(I)*DXANC(1,IX1L(I))
1922 . -H2(I)*DXANC(1,IX2L(I))
1923 . -H3(I)*DXANC(1,IX3L(I))
1924 . -H4(I)*DXANC(1,IX4L(I)))
1925 . + N2(I)*(DXANCFI(NIN)%P(2,IL)-H1(I)*DXANC(2,IX1L(I))
1926 . -H2(I)*DXANC(2,IX2L(I))
1927 . -H3(I)*DXANC(2,IX3L(I))
1928 . -H4(I)*DXANC(2,IX4L(I)))
1929 . + N3(I)*(DXANCFI(NIN)%P(3,IL)-H1(I)*DXANC(3,IX1L(I))
1930 . -H2(I)*DXANC(3,IX2L(I))
1931 . -H3(I)*DXANC(3,IX3L(I))
1932 . -H4(I)*DXANC(3,IX4L(I)))
1933 END IF
1934 PS = PENE(I) - XSA - GAPV(I) + GAPR(I)
1935 IF(PS <= ZERO)THEN
1936 STIF(I) = ZERO
1937 FXI(I) = ZERO
1938 FYI(I) = ZERO
1939 FZI(I) = ZERO
1940 FX1(I) = ZERO
1941 FY1(I) = ZERO
1942 FZ1(I) = ZERO
1943 FX2(I) = ZERO
1944 FY2(I) = ZERO
1945 FZ2(I) = ZERO
1946 FX3(I) = ZERO
1947 FY3(I) = ZERO
1948 FZ3(I) = ZERO
1949 FX4(I) = ZERO
1950 FY4(I) = ZERO
1951 FZ4(I) = ZERO
1952 IF (IFQ>0) THEN
1953 CAND_FX(INDEX(I)) = ZERO
1954 CAND_FY(INDEX(I)) = ZERO
1955 CAND_FZ(INDEX(I)) = ZERO
1956
1957 ENDIF
1958 ENDIF
1959 ENDIF
1960 ENDDO
1961
1962
1963
1964
1965.OR. IF(INTTH == 0 IFORM == 0) THEN
1966 DO I=1,JLT
1967 PHI1(I) = ZERO
1968 PHI2(I) = ZERO
1969 PHI3(I) = ZERO
1970 PHI4(I) = ZERO
1971
1972 ENDDO
1973 ELSEIF(IFORM > 0) THEN
1974 DO I=1,JLT
1975 TM = H1(I)*TEMP(IX1G(I)) + H2(I)*TEMP(IX2G(I))
1976 . + H3(I)*TEMP(IX3G(I)) + H4(I)*TEMP(IX4G(I))
1977
1978 TS = TEMPI(I)
1979
1980 AX1 = XA(1,IX3L(I)) - XA(1,IX1L(I))
1981 AY1 = XA(2,IX3L(I)) - XA(2,IX1L(I))
1982 AZ1 = XA(3,IX3L(I)) - XA(3,IX1L(I))
1983 AX2 = XA(1,IX4L(I)) - XA(1,IX2L(I))
1984 AY2 = XA(2,IX4L(I)) - XA(2,IX2L(I))
1985 AZ2 = XA(3,IX4L(I)) - XA(3,IX2L(I))
1986
1987 AX = AY1*AZ2 - AZ1*AY2
1988 AY = AZ1*AX2 - AX1*AZ2
1989 AZ = AX1*AY2 - AY1*AX2
1990
1991 AREA = ONE_OVER_8*SQRT(AX*AX+AY*AY+AZ*AZ)
1992 PHI(I) = AREA* (TM - TS)*DT1 / RSTIF
1993 PHI1(I) = -PHI(I) *H1(I)
1994 PHI2(I) = -PHI(I) *H2(I)
1995 PHI3(I) = -PHI(I) *H3(I)
1996 PHI4(I) = -PHI(I) *H4(I)
1997 ENDDO
1998 ENDIF
1999
2000 IF (NSPMD>1) THEN
2001
2002#include "mic_lockon.inc"
2003 DO I = 1,JLT
2004 NN = NSVG(I)
2005 IF(NN<0)THEN
2006
2007 NSVFI(NIN)%P(-NN) = -ABS(NSVFI(NIN)%P(-NN))
2008 ENDIF
2009 ENDDO
2010
2011#include "mic_lockoff.inc"
2012 ENDIF
2013
2014.OR. IF(IDTMINS==2IDTMINS_INT/=0)THEN
2015 DTI=DT2T
2016 CALL I7SMS2(JLT ,IX1G ,IX2G ,IX3G ,IX4G ,
2017 2 NSVG ,H1 ,H2 ,H3 ,H4 ,STIF ,
2018 3 NIN ,NOINT ,MSKYI_SMS, ISKYI_SMS,NSMS ,
2019 4 KT ,C ,CF ,DTMINI,DTI )
2020 IF(DTI<DT2T)THEN
2021 DT2T = DTI
2022 NELTST = NOINT
2023 ITYPTST = 10
2024 ENDIF
2025 ENDIF
2026
2027 IF(IDTMINS_INT/=0)THEN
2028 STIF(1:JLT)=ZERO
2029 END IF
2030
2031 BID = ZERO
2032 IBID =0
2033 IF(IPARIT==3)THEN
2034 IF(KDTINT==0)THEN
2035 CALL I7ASS3(JLT ,IX1G ,IX2G ,IX3G ,IX4G ,
2036 2 NSVG ,H1 ,H2 ,H3 ,H4 ,STIF ,
2037 3 FX1 ,FY1 ,FZ1 ,FX2 ,FY2 ,FZ2 ,
2038 4 FX3 ,FY3 ,FZ3 ,FX4 ,FY4 ,FZ4 ,
2039 5 FXI ,FYI ,FZI ,A ,STIFN)
2040 ELSE
2041 CALL I7ASS35(JLT ,IX1G ,IX2G ,IX3G ,IX4G ,
2042 2 NSVG ,H1 ,H2 ,H3 ,H4 ,STIF ,
2043 3 FX1 ,FY1 ,FZ1 ,FX2 ,FY2 ,FZ2 ,
2044 4 FX3 ,FY3 ,FZ3 ,FX4 ,FY4 ,FZ4 ,
2045 5 FXI ,FYI ,FZI ,A ,STIFN,VISCN,
2046 6 KS ,K1 ,K2 ,K3 ,K4 ,CS ,
2047 7 C1 ,C2 ,C3 ,C4 )
2048 ENDIF
2049 ELSEIF(IPARIT==0)THEN
2050 IF(KDTINT==0)THEN
2051 CALL I7ASS0(JLT ,IX1G ,IX2G ,IX3G ,IX4G ,
2052 2 NSVG ,H1 ,H2 ,H3 ,H4 ,STIF ,
2053 3 FX1 ,FY1 ,FZ1 ,FX2 ,FY2 ,FZ2 ,
2054 4 FX3 ,FY3 ,FZ3 ,FX4 ,FY4 ,FZ4 ,
2055 5 FXI ,FYI ,FZI ,A ,STIFN ,NIN ,
2056 6 INTTH ,PHI ,FTHE ,PHI1 , PHI2 ,PHI3 ,
2057 7 PHI4 ,BID ,BID ,JTASK,IBID ,IBID )
2058
2059 ELSE
2060
2061 CALL I7ASS05(JLT ,IX1G ,IX2G ,IX3G ,IX4G ,
2062 2 NSVG ,H1 ,H2 ,H3 ,H4 ,
2063 3 FX1 ,FY1 ,FZ1 ,FX2 ,FY2 ,FZ2 ,
2064 4 FX3 ,FY3 ,FZ3 ,FX4 ,FY4 ,FZ4 ,
2065 5 FXI ,FYI ,FZI ,A ,STIFN ,VISCN ,
2066 6 KS ,K1 ,K2 ,K3 ,K4 ,CS ,
2067 7 C1 ,C2 ,C3 ,C4 ,NIN ,INTTH ,
2068 8 PHI ,FTHE ,PHI1 , PHI2 ,PHI3 , PHI4,JTASK,
2069 9 BID ,BID ,IBID ,IBID )
2070 ENDIF
2071
2072 ELSE
2073 IF(KDTINT==0)THEN
2074 CALL I7ASS2(JLT ,IX1G ,IX2G ,IX3G ,IX4G ,
2075 2 NSVG ,H1 ,H2 ,H3 ,H4 ,STIF ,
2076 3 FX1 ,FY1 ,FZ1 ,FX2 ,FY2 ,FZ2 ,
2077 4 FX3 ,FY3 ,FZ3 ,FX4 ,FY4 ,FZ4 ,
2078 5 FXI ,FYI ,FZI ,FSKYI,ISKY ,NISKYFI,
2079 6 NIN ,NOINT ,INTTH,PHI ,FTHESKYI,PHI1,
2080 7 PHI2 ,PHI3 , PHI4,BID ,BID ,
2081 A IBID ,IBID )
2082 ELSE
2083 CALL I7ASS25(JLT ,IX1G ,IX2G ,IX3G ,IX4G ,
2084 2 NSVG ,H1 ,H2 ,H3 ,H4 ,
2085 3 FX1 ,FY1 ,FZ1 ,FX2 ,FY2 ,FZ2 ,
2086 4 FX3 ,FY3 ,FZ3 ,FX4 ,FY4 ,FZ4 ,
2087 5 FXI ,FYI ,FZI ,FSKYI,NISKYFI,NIN ,
2088 6 KS ,K1 ,K2 ,K3 ,K4 ,CS ,
2089 7 C1 ,C2 ,C3 ,C4 ,ISKY ,NOINT ,
2090 8 INTTH ,PHI ,FTHESKYI,PHI1,PHI2 ,PHI3,
2091 9 PHI4 ,BID ,BID ,IBID ,IBID )
2092 ENDIF
2093 ENDIF
2094
2095 IF(ANIM_V(4)+OUTP_V(4)+H3D_DATA%N_VECT_CONT>0)THEN
2096#include "lockon.inc"
2097
2098 DO I=1,JLT
2099 FCONT(1,IX1G(I)) =FCONT(1,IX1G(I)) + FX1(I)
2100 FCONT(2,IX1G(I)) =FCONT(2,IX1G(I)) + FY1(I)
2101 FCONT(3,IX1G(I)) =FCONT(3,IX1G(I)) + FZ1(I)
2102 FCONT(1,IX2G(I)) =FCONT(1,IX2G(I)) + FX2(I)
2103 FCONT(2,IX2G(I)) =FCONT(2,IX2G(I)) + FY2(I)
2104 FCONT(3,IX2G(I)) =FCONT(3,IX2G(I)) + FZ2(I)
2105 FCONT(1,IX3G(I)) =FCONT(1,IX3G(I)) + FX3(I)
2106 FCONT(2,IX3G(I)) =FCONT(2,IX3G(I)) + FY3(I)
2107 FCONT(3,IX3G(I)) =FCONT(3,IX3G(I)) + FZ3(I)
2108 FCONT(1,IX4G(I)) =FCONT(1,IX4G(I)) + FX4(I)
2109 FCONT(2,IX4G(I)) =FCONT(2,IX4G(I)) + FY4(I)
2110 FCONT(3,IX4G(I)) =FCONT(3,IX4G(I)) + FZ4(I)
2111 JG = NSVG(I)
2112 IF(JG>0) THEN
2113
2114 FCONT(1,JG)=FCONT(1,JG)- FXI(I)
2115 FCONT(2,JG)=FCONT(2,JG)- FYI(I)
2116 FCONT(3,JG)=FCONT(3,JG)- FZI(I)
2117 ENDIF
2118 ENDDO
2119
2120#include "lockoff.inc"
2121 ENDIF
2122
2123 IF(ISECIN>0)THEN
2124 K0=NSTRF(25)
2125 IF(NSTRF(1)+NSTRF(2)/=0)THEN
2126 DO I=1,NSECT
2127 NBINTER=NSTRF(K0+14)
2128 K1S=K0+30
2129 DO J=1,NBINTER
2130 IF(NSTRF(K1S)==NOINT)THEN
2131 IF(ISECUT/=0)THEN
2132#include "lockon.inc"
2133 DO K=1,JLT
2134
2135
2136 IF(SECFCUM(4,IX1G(K),I)==1.)THEN
2137 SECFCUM(1,IX1G(K),I)=SECFCUM(1,IX1G(K),I)-FX1(K)
2138 SECFCUM(2,IX1G(K),I)=SECFCUM(2,IX1G(K),I)-FY1(K)
2139 SECFCUM(3,IX1G(K),I)=SECFCUM(3,IX1G(K),I)-FZ1(K)
2140 ENDIF
2141 IF(SECFCUM(4,IX2G(K),I)==1.)THEN
2142 SECFCUM(1,IX2G(K),I)=SECFCUM(1,IX2G(K),I)-FX2(K)
2143 SECFCUM(2,IX2G(K),I)=SECFCUM(2,IX2G(K),I)-FY2(K)
2144 SECFCUM(3,IX2G(K),I)=SECFCUM(3,IX2G(K),I)-FZ2(K)
2145 ENDIF
2146 IF(SECFCUM(4,IX3G(K),I)==1.)THEN
2147 SECFCUM(1,IX3G(K),I)=SECFCUM(1,IX3G(K),I)-FX3(K)
2148 SECFCUM(2,IX3G(K),I)=SECFCUM(2,IX3G(K),I)-FY3(K)
2149 SECFCUM(3,IX3G(K),I)=SECFCUM(3,IX3G(K),I)-FZ3(K)
2150 ENDIF
2151 IF(SECFCUM(4,IX4G(K),I)==1.)THEN
2152 SECFCUM(1,IX4G(K),I)=SECFCUM(1,IX4G(K),I)-FX4(K)
2153 SECFCUM(2,IX4G(K),I)=SECFCUM(2,IX4G(K),I)-FY4(K)
2154 SECFCUM(3,IX4G(K),I)=SECFCUM(3,IX4G(K),I)-FZ4(K)
2155 ENDIF
2156
2157 JG = NSVG(K)
2158 IF(JG>0) THEN
2159
2160 IF(SECFCUM(4,JG,I)==1.)THEN
2161 SECFCUM(1,JG,I)=SECFCUM(1,JG,I)+FXI(K)
2162 SECFCUM(2,JG,I)=SECFCUM(2,JG,I)+FYI(K)
2163 SECFCUM(3,JG,I)=SECFCUM(3,JG,I)+FZI(K)
2164 ENDIF
2165 ENDIF
2166
2167 ENDDO
2168#include "lockoff.inc"
2169 ENDIF
2170
2171 ENDIF
2172 K1S=K1S+1
2173 ENDDO
2174 K0=NSTRF(K0+24)
2175 ENDDO
2176 ENDIF
2177 ENDIF
2178
2179
2180.OR. IF(IBAG/=0IADM/=0)THEN
2181 DO I=1,JLT
2182
2183
2184
2185.OR..OR. IF(FXI(I)/=ZEROFYI(I)/=ZEROFZI(I)/=ZERO)THEN
2186
2187 JG = NSVG(I)
2188 IF(JG>0) THEN
2189
2190 ICONTACT(JG)=1
2191 ENDIF
2192
2193 ICONTACT(IX1G(I))=1
2194 ICONTACT(IX2G(I))=1
2195 ICONTACT(IX3G(I))=1
2196 ICONTACT(IX4G(I))=1
2197 ENDIF
2198 ENDDO
2199 ENDIF
2200
2201 IF(IADM/=0)THEN
2202 DO I=1,JLT
2203 JG = NSVG(I)
2204#include "lockon.inc"
2205 IF(JG>0) THEN
2206
2207 RCONTACT(JG)=MIN(RCONTACT(JG),RCURVI(I))
2208 END IF
2209 RCONTACT(IX1G(I))=MIN(RCONTACT(IX1G(I)),RCURVI(I))
2210 RCONTACT(IX2G(I))=MIN(RCONTACT(IX2G(I)),RCURVI(I))
2211 RCONTACT(IX3G(I))=MIN(RCONTACT(IX3G(I)),RCURVI(I))
2212 RCONTACT(IX4G(I))=MIN(RCONTACT(IX4G(I)),RCURVI(I))
2213#include "lockoff.inc"
2214 END DO
2215 END IF
2216 IF(IADM>=2)THEN
2217 DO I=1,JLT
2218 JG = NSVG(I)
2219#include "lockon.inc"
2220 IF(JG>0) THEN
2221
2222 PCONTACT(JG)=MAX(PCONTACT(JG),PENE(I)/(PADM*GAPV(I)))
2223 ACONTACT(JG)=MIN(ACONTACT(JG),ANGLMI(I))
2224 END IF
2225#include "lockoff.inc"
2226 END DO
2227 END IF
2228
2229 IF(IBCC==0) RETURN
2230
2231 DO 400 I=1,JLT
2232
2233 IF(PENE(I)==ZERO)GOTO 400
2234 IBCM = IBCC / 8
2235 IBCS = IBCC - 8 * IBCM
2236 IF(IBCS>0) THEN
2237 IG=NSVG(I)
2238 IF(IG>0) THEN
2239
2240 CALL IBCOFF(IBCS,ICODT(IG))
2241 ENDIF
2242 ENDIF
2243 IF(IBCM>0) THEN
2244 IG=IX1G(I)
2245 CALL IBCOFF(IBCM,ICODT(IG))
2246 IG=IX2G(I)
2247 CALL IBCOFF(IBCM,ICODT(IG))
2248 IG=IX3G(I)
2249 CALL IBCOFF(IBCM,ICODT(IG))
2250 IG=IX4G(I)
2251 CALL IBCOFF(IBCM,ICODT(IG))
2252 ENDIF
2253 400 CONTINUE
2254
2255 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
subroutine i7curv(jlt, pene, n1, n2, n3, gapv, x, nod_normal, ix1, ix2, ix3, ix4, h1, h2, h3, h4, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi)
type(real_pointer2), dimension(:), allocatable fnconti
type(real_pointer), dimension(:), allocatable stifi
type(real_pointer2), dimension(:), allocatable penfi
type(real_pointer), dimension(:), allocatable alphakfi
type(int_pointer), dimension(:), allocatable lisubsfi
type(real_pointer), dimension(:), allocatable gapfi
type(int_pointer), dimension(:), allocatable addsubsfi
type(real_pointer2), dimension(:), allocatable ftconti
type(int_pointer), dimension(:), allocatable itafi
int main(int argc, char *argv[])
1.15.0