OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i7for3.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| i7for3 ../engine/source/interfaces/int07/i7for3.f
25!||--- called by ------------------------------------------------------
26!|| i7mainf ../engine/source/interfaces/int07/i7mainf.F
27!||--- calls -----------------------------------------------------
28!|| bitget ../engine/source/interfaces/intsort/i20sto.F
29!|| i7curv ../engine/source/interfaces/int07/i7curv.F
30!|| ibcoff ../engine/source/interfaces/interf/ibcoff.F
31!||--- uses -----------------------------------------------------
32!|| h3d_mod ../engine/share/modules/h3d_mod.F
33!|| output_mod ../common_source/modules/output/output_mod.F90
34!|| tri7box ../engine/share/modules/tri7box.f
35!||====================================================================
36 SUBROUTINE i7for3(OUTPUT,JLT ,A ,V ,IBCC ,ICODT ,
37 2 FSAV ,GAP ,FRIC ,MS ,VISC ,
38 3 VISCF ,NOINT ,STFN ,ITAB ,CN_LOC ,
39 4 STIGLO ,STIFN ,STIF ,FSKYI ,ISKY ,
40 5 NX1 ,NX2 ,NX3 ,NX4 ,NY1 ,
41 6 NY2 ,NY3 ,NY4 ,NZ1 ,NZ2 ,
42 7 NZ3 ,NZ4 ,LB1 ,LB2 ,LB3 ,
43 8 LB4 ,LC1 ,LC2 ,LC3 ,LC4 ,
44 9 P1 ,P2 ,P3 ,P4 ,FCONT ,
45 A IX1 ,IX2 ,IX3 ,IX4 ,NSVG ,
46 B IVIS2 ,NELTST ,ITYPTST ,DT2T ,GAPV ,
47 C INACTI ,CAND_P ,INDEX ,KINET ,NEWFRONT,
48 D ISECIN ,NSTRF ,X ,IRECT ,CE_LOC ,
49 E MFROT ,IFQ ,CAND_FX ,CAND_FY ,
50 F CAND_FZ ,ALPHA0 ,IFPEN ,IBAG ,
51 H ICONTACT ,VISCN ,VXI ,VYI ,VZI ,
52 I MSI ,KINI ,NIN ,NISUB ,LISUB ,
53 J ADDSUBS ,ADDSUBM ,LISUBS ,LISUBM ,FSAVSUB ,
54 K CAND_N ,ILAGM ,ICURV ,NOD_NORMAL ,FNCONT ,
55 L FTCONT ,X1 ,X2 ,X3 ,X4 ,
56 M Y1 ,Y2 ,Y3 ,Y4 ,Z1 ,
57 N Z2 , Z3 ,Z4 ,XI ,YI ,
58 O ZI , IADM ,RCURVI ,RCONTACT ,ACONTACT,
59 P PCONTACT ,ANGLMI ,PADM ,INTTH ,TEMP ,
60 Q TEMPI ,IFORM ,NPC ,TF ,CMAJ ,
61 R DTMINI ,DRAD ,FHEATS ,FHEATM ,EFRICT ,
62 S QFRIC ,FNI ,IFRIC ,
63 T JTASK ,H1 ,
64 U H2 ,H3 ,H4 ,KS ,KT ,
65 V K1 ,K2 ,K3 ,K4 ,C1 ,
66 W C2 ,C3 ,C4 ,CS ,C ,
67 X CF ,TINT ,XFRIC ,FXI ,FYI ,
68 Y FZI ,FX1 ,FY1 ,FZ1 ,FX2 ,
69 Z FY2 ,FZ2 ,FX3 ,FY3 ,FZ3 ,
70 1 FX4 ,FY4 ,FZ4 ,ISENSINT ,FSAVPARIT,
71 5 NFT ,SYM_FLAG_TYPE19,H3D_DATA,FRICC ,VISCFFRIC,
72 6 FRIC_COEFS ,ITIED ,JLT_TIED ,CAND_F ,IORTHFRIC,
73 7 FRIC_COEFS2,FRICC2 ,VISCFFRIC2,NFORTH ,NFISOT ,
74 8 INDEXORTH ,INDEXISOT ,DIR1 ,DIR2 ,TAGNCONT ,
75 9 KLOADPINTER,LOADPINTER,LOADP_HYD_INTER,TYPSUB,INFLG_SUBS,
76 A INFLG_SUBM,NINLOADP ,DGAPLOADINT,S_LOADPINTER,DGAPLOADP,
77 B INTEREFRIC )
78
79C-----------------------------------------------
80C M o d u l e s
81C-----------------------------------------------
82 USE tri7box
83 USE h3d_mod
84 USE output_mod
85C-----------------------------------------------
86C I m p l i c i t T y p e s
87C-----------------------------------------------
88#include "implicit_f.inc"
89#include "comlock.inc"
90C-----------------------------------------------
91C G l o b a l P a r a m e t e r s
92C-----------------------------------------------
93#include "mvsiz_p.inc"
94C-----------------------------------------------
95C C o m m o n B l o c k s
96C-----------------------------------------------
97#include "com01_c.inc"
98#include "com04_c.inc"
99#include "com06_c.inc"
100#include "com08_c.inc"
101#include "scr05_c.inc"
102#include "scr07_c.inc"
103#include "scr11_c.inc"
104#include "scr14_c.inc"
105#include "scr16_c.inc"
106#include "scr18_c.inc"
107#include "sms_c.inc"
108#include "units_c.inc"
109#include "parit_c.inc"
110#include "param_c.inc"
111#include "impl1_c.inc"
112#include "kincod_c.inc"
113C-----------------------------------------------
114C D u m m y A r g u m e n t s
115C-----------------------------------------------
116 TYPE(output_), INTENT(INOUT) :: OUTPUT
117 INTEGER NELTST,ITYPTST,JLT,IBCC,IBCM,IBCS,IVIS2,INACTI,IBAG,NIN,
118 . IFRIC,IFORM,INTFRIC,ITIED,JLT_TIED,IORTHFRIC,NFORTH ,NFISOT ,
119 . ICODT(*), ITAB(*), ISKY(*), KINET(*),
120 . MFROT, IFQ, NOINT,NEWFRONT,ISECIN, NSTRF(*),
121 . IRECT(4,*),IFPEN(*) ,ICONTACT(*), CAND_N(*),
122 . KINI(*),NPC(*),JTASK,TAGNCONT(NLOADP_HYD_INTER,NUMNOD),
123 . typsub(*),
124 . iset, iadm,intth,nft
125 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
126 . CN_LOC(MVSIZ),CE_LOC(MVSIZ),INDEX(MVSIZ),NSVG(MVSIZ),
127 . NISUB, LISUB(*), ADDSUBS(*), ADDSUBM(*), LISUBS(*),
128 . LISUBM(*),ILAGM,ICURV(3), ISENSINT(*),SYM_FLAG_TYPE19,
129 . INDEXORTH(MVSIZ) ,INDEXISOT(MVSIZ),INFLG_SUBS(*), INFLG_SUBM(*)
130 INTEGER , INTENT(IN) :: NINLOADP,S_LOADPINTER
131 INTEGER , INTENT(IN) :: KLOADPINTER(NINTER+1),LOADPINTER(S_LOADPINTER),
132 . LOADP_HYD_INTER(NLOADP_HYD)
133 INTEGER , INTENT(IN) :: INTEREFRIC
134 my_real
135 . STIGLO,CAND_P(*), X(3,*),
136 . A(3,*), MS(*), V(3,*), FSAV(*),FCONT(3,*),
137 . CAND_FX(*),CAND_FY(*),CAND_FZ(*),ALPHA0,
138 . GAP, FRIC,VISC,VISCF,VIS,DT2T,STFN(*),STIFN(*),
139 . FSKYI(LSKYI,NFSKYI),FSAVSUB(NTHVKI,*), FNCONT(3,*),FTCONT(3,*),
140 . DTMINI,FHEATS,FHEATM,
141 . FSAVPARIT(NISUB+1,11,*), CAND_F(8,*)
142 my_real
143 . NX1(MVSIZ), NX2(MVSIZ), NX3(MVSIZ), NX4(MVSIZ),
144 . NY1(MVSIZ), NY2(MVSIZ), NY3(MVSIZ), NY4(MVSIZ),
145 . NZ1(MVSIZ), NZ2(MVSIZ), NZ3(MVSIZ), NZ4(MVSIZ),
146 . LB1(MVSIZ), LB2(MVSIZ), LB3(MVSIZ), LB4(MVSIZ),
147 . LC1(MVSIZ), LC2(MVSIZ), LC3(MVSIZ), LC4(MVSIZ),
148 . P1(MVSIZ), P2(MVSIZ), P3(MVSIZ), P4(MVSIZ), STIF(MVSIZ),
149 . GAPV(MVSIZ),
150 . TMP(MVSIZ),
151 . STIFSAV(MVSIZ), VISCN(*),TF(*),
152 . VXI(MVSIZ),VYI(MVSIZ),VZI(MVSIZ),MSI(MVSIZ),
153 . X1(MVSIZ),Y1(MVSIZ),Z1(MVSIZ),
154 . X2(MVSIZ),Y2(MVSIZ),Z2(MVSIZ),
155 . X3(MVSIZ),Y3(MVSIZ),Z3(MVSIZ),
156 . X4(MVSIZ),Y4(MVSIZ),Z4(MVSIZ),
157 . xi(mvsiz),yi(mvsiz),zi(mvsiz),
158 . temp(*), tempi(mvsiz),
159 . fni(mvsiz), fns(mvsiz)
160 my_real
161 . nod_normal(3,*), rcurvi(*), rcontact(*), acontact(*),
162 . pcontact(*),padm, anglmi(*),cmaj(mvsiz), efrict(mvsiz),
163 . qfric,drad,tint,xfric,fric_coefs(mvsiz,10), fricc(mvsiz),
164 . viscffric(mvsiz),fric_coefs2(mvsiz,10),fricc2(mvsiz),
165 . viscffric2(mvsiz),dir1(mvsiz,3) ,dir2(mvsiz,3)
166 my_real
167 . ks(mvsiz),k1(mvsiz),k2(mvsiz),k3(mvsiz),k4(mvsiz),
168 . cs(mvsiz),c1(mvsiz),c2(mvsiz),c3(mvsiz),c4(mvsiz),
169 . kt(mvsiz),c(mvsiz),cf(mvsiz),
170 . fxi(mvsiz), fyi(mvsiz), fzi(mvsiz),
171 . fx1(mvsiz), fx2(mvsiz), fx3(mvsiz), fx4(mvsiz),
172 . fy1(mvsiz), fy2(mvsiz), fy3(mvsiz), fy4(mvsiz),
173 . fz1(mvsiz), fz2(mvsiz), fz3(mvsiz), fz4(mvsiz)
174 TYPE(h3d_database) :: H3D_DATA
175 my_real , INTENT(IN) :: dgaploadp, dgaploadint(s_loadpinter)
176C-----------------------------------------------
177C L o c a l V a r i a b l e s
178C-----------------------------------------------
179 INTEGER I, J1, IG, J, JG , K0,NBINTER,K1S,K,IL,IE, NN, NI, II,
180 . NA1,NA2,IPROJ,IDTM,PP,PPL,ITYPSUB,ISS1,ISS2,IMS1,IMS2
181 INTEGER JSUB,KSUB,JJ,KK,IN,NSUB
182 my_real
183 . fxr(mvsiz), fyr(mvsiz), fzr(mvsiz)
184 my_real
185 . fxt(mvsiz),fyt(mvsiz),fzt(mvsiz),
186 . n1(mvsiz), n2(mvsiz), n3(mvsiz), pene(mvsiz),
187 . vis2(mvsiz), dtmi(mvsiz), xmu(mvsiz),stif0(mvsiz),
188 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
189 . vx(mvsiz), vy(mvsiz), vz(mvsiz), vn(mvsiz),
190 . xp(mvsiz), yp(mvsiz), zp(mvsiz),efric_l(mvsiz),
191 . vnx, vny, vnz, aa, crit,s2,dist,rdist,
192 . v2, fm2, dt1inv, visca, fac,ff,alphi,alpha,beta,
193 . fx, fy, fz, f2, mas2, m2sk, dtmi0,dti,ft,fn,fmax,ftn,
194 . facm1, econtt, econvt, h0, la1, la2, la3, la4,
195 . d1,d2,d3,d4,a1,a2,a3,a4,econtdt,econttied,
196 . fsav1, fsav2, fsav3, fsav4, fsav5, fsav6, fsav7, fsav8,
197 . fsav9, fsav10, fsav11, fsav12, fsav13, fsav14, fsav15,
198 . fsav22, fsav23, fsav24, ffo,
199 . e10, h0d, s2d, sum,
200 . la1d,la2d,la3d,la4d,t1,t1d,t2,t2d,ffd,visd,facd,d1d,
201 . p1s(mvsiz),p2s(mvsiz),p3s(mvsiz),p4s(mvsiz),
202 . d2d,d3d,d4d,vnxd,vnyd,vnzd,v2d,fm2d,f2d,aad,fxd,fyd,fzd,
203 . a1d,a2d,a3d,a4d,vv,ax1,ax2,ay1,ay2,az1,az2,ax,ay,az,
204 . area,p,vv1,vv2,v21,dmu, dti2,h00 ,a0x,a0y,a0z,rx,ry,rz,
205 . anx,any,anz,aan,aax,aay,aaz ,rr,rs,aaa,stfr,visr ,tm,ts
206 my_real
207 . prec,nn1,nn2,nn3,nn4,xn1,yn1,zn1,xn2,yn2,zn2,xn3,yn3,zn3,
208 . xn4,yn4,zn4,efric_ls,efric_lm
209 my_real
210 . st1(mvsiz),st2(mvsiz),st3(mvsiz),st4(mvsiz),stv(mvsiz),
211 . cx,cy,cfi,aux, fx6(6,mvsiz), fy6(6,mvsiz), fz6(6,mvsiz),
212 . dtm, phim,qfrict,dydx,th,thi(mvsiz),frict(mvsiz),forc_sign,
213 . vt1(mvsiz), vt2(mvsiz),
214 . ft1(mvsiz),ft2(mvsiz),
215 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
216 . t1x(mvsiz),t1y(mvsiz),t1z(mvsiz),
217 . t2x(mvsiz),t2y(mvsiz),t2z(mvsiz), dt05, norminv,
218 . xh, yh, zh, side(mvsiz),xmu2(mvsiz),csa ,sna ,alphaf ,ftt1 ,ftt2,
219 . an(nforth) ,bn(nforth) ,nep1 ,nep2 ,nep ,c11 ,c22 ,ans , ep ,signc
220 my_real
221 . fsavsub1(24,nisub),impx,impy,impz,bb,gapr,dgap,gapp
222 my_real
223 . total
224 my_real distp(mvsiz)
225 INTEGER BITGET
226 EXTERNAL BITGET
227C
228C-----------------------------------------------
229C 1 .. JLT : all impacts
230C 1 .. JLT-JLT_TIED : Standard TYPE7 formulation
231C JLT-JLT_TIED+1 .. JLT : Modified formulation (tied nodes with possible rebound)
232C-----------------------------------------------
233 IF (IRESP==1) then
234 prec = fiveem4
235 ELSE
236 prec = em10
237 ENDIF
238 IF(dt1>zero)THEN
239 dt1inv = one/dt1
240 ELSE
241 dt1inv =zero
242 ENDIF
243 econtt = zero
244 econvt = zero
245 qfrict = zero
246 econtdt = zero
247 econttied = zero
248 DO i=1,jlt
249 stif0(i) = stif(i)
250 ENDDO
251 efrict(1:jlt) = zero
252 efric_l(1:jlt) = zero
253C
254C--------------------------------------------------------
255C For sym type7 of type19 interface normal and tangential forces must have opposite sign
256C--------------------------------------------------------
257 forc_sign = one
258 IF (sym_flag_type19 > 0) forc_sign = -one
259C
260C--------------------------------------------------------
261C case of mixed packets
262C--------------------------------------------------------
263 IF((intth == 0.OR.drad == zero).AND.dgaploadp==zero) THEN ! NO HEAT EXCHANGE/No load pressure
264 IF(icurv(1) == 3) THEN
265 DO i=1,jlt-jlt_tied
266 bb = gapv(i)+cmaj(i)
267C
268 d1 = sqrt(p1(i))
269 p1(i) = max(zero, bb - d1)
270C
271 d2 = sqrt(p2(i))
272 p2(i) = max(zero, bb - d2)
273C
274 d3 = sqrt(p3(i))
275 p3(i) = max(zero, bb - d3)
276C
277 d4 = sqrt(p4(i))
278 p4(i) = max(zero, bb - d4)
279C
280 a1 = p1(i)/max(em20,d1)
281 a2 = p2(i)/max(em20,d2)
282 a3 = p3(i)/max(em20,d3)
283 a4 = p4(i)/max(em20,d4)
284 n1(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i)
285 n2(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i)
286 n3(i) = a1*nz1(i) + a2*nz2(i) + a3*nz3(i) + a4*nz4(i)
287 ENDDO
288 ELSE
289 DO i=1,jlt-jlt_tied
290C
291 d1 = sqrt(p1(i))
292 p1(i) = max(zero, gapv(i) - d1)
293C
294 d2 = sqrt(p2(i))
295 p2(i) = max(zero, gapv(i) - d2)
296C
297 d3 = sqrt(p3(i))
298 p3(i) = max(zero, gapv(i) - d3)
299C
300 d4 = sqrt(p4(i))
301 p4(i) = max(zero, gapv(i) - d4)
302C
303 a1 = p1(i)/max(em20,d1)
304 a2 = p2(i)/max(em20,d2)
305 a3 = p3(i)/max(em20,d3)
306 a4 = p4(i)/max(em20,d4)
307 n1(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i)
308 n2(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i)
309 n3(i) = a1*nz1(i) + a2*nz2(i) + a3*nz3(i) + a4*nz4(i)
310 ENDDO
311 ENDIF
312C
313 DO i=1,jlt-jlt_tied
314 IF(ix3(i)/=ix4(i))THEN
315 pene(i) = max(p1(i),p2(i),p3(i),p4(i))
316C
317 la1 = one - lb1(i) - lc1(i)
318 la2 = one - lb2(i) - lc2(i)
319 la3 = one - lb3(i) - lc3(i)
320 la4 = one - lb4(i) - lc4(i)
321C
322 h0 = fourth *
323 . (p1(i)*la1 + p2(i)*la2 + p3(i)*la3 + p4(i)*la4)
324 h1(i) = h0 + p1(i) * lb1(i) + p4(i) * lc4(i)
325 h2(i) = h0 + p2(i) * lb2(i) + p1(i) * lc1(i)
326 h3(i) = h0 + p3(i) * lb3(i) + p2(i) * lc2(i)
327 h4(i) = h0 + p4(i) * lb4(i) + p3(i) * lc3(i)
328 h00 = one/max(em20,h1(i) + h2(i) + h3(i) + h4(i))
329 h1(i) = h1(i) * h00
330 h2(i) = h2(i) * h00
331 h3(i) = h3(i) * h00
332 h4(i) = h4(i) * h00
333C
334 ELSE
335 pene(i) = p1(i)
336 n1(i) = nx1(i)
337 n2(i) = ny1(i)
338 n3(i) = nz1(i)
339 h1(i) = lb1(i)
340 h2(i) = lc1(i)
341 h3(i) = one - lb1(i) - lc1(i)
342 h4(i) = zero
343 ENDIF
344 ENDDO
345 ELSE ! HEAT EXCHANGE OR LOADP PRESSURE
346 IF(icurv(1) == 3) THEN
347 DO i=1,jlt-jlt_tied
348 bb = gapv(i)+cmaj(i)+dgaploadp
349 gapr = max(bb,drad)
350C
351 d1 = sqrt(p1(i))
352 p1(i) = max(zero, gapr - d1)
353C
354 d2 = sqrt(p2(i))
355 p2(i) = max(zero, gapr - d2)
356C
357 d3 = sqrt(p3(i))
358 p3(i) = max(zero, gapr - d3)
359C
360 d4 = sqrt(p4(i))
361 p4(i) = max(zero, gapr - d4)
362C
363 distp(i) = max(d1,d2,d3,d4)
364 a1 = p1(i)/max(em20,d1)
365 a2 = p2(i)/max(em20,d2)
366 a3 = p3(i)/max(em20,d3)
367 a4 = p4(i)/max(em20,d4)
368 n1(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i)
369 n2(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i)
370 n3(i) = a1*nz1(i) + a2*nz2(i) + a3*nz3(i) + a4*nz4(i)
371 ENDDO
372 ELSE
373 DO i=1,jlt-jlt_tied
374 gapr = max(gapv(i)+dgaploadp,drad)
375C
376 d1 = sqrt(p1(i))
377 p1(i) = max(zero, gapr - d1)
378C
379 d2 = sqrt(p2(i))
380 p2(i) = max(zero, gapr - d2)
381C
382 d3 = sqrt(p3(i))
383 p3(i) = max(zero, gapr - d3)
384C
385 d4 = sqrt(p4(i))
386 p4(i) = max(zero, gapr - d4)
387C
388 distp(i) = max(d1,d2,d3,d4)
389C
390 a1 = p1(i)/max(em20,d1)
391 a2 = p2(i)/max(em20,d2)
392 a3 = p3(i)/max(em20,d3)
393 a4 = p4(i)/max(em20,d4)
394 n1(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i)
395 n2(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i)
396 n3(i) = a1*nz1(i) + a2*nz2(i) + a3*nz3(i) + a4*nz4(i)
397 ENDDO
398 ENDIF
399C
400 DO i=1,jlt-jlt_tied
401 IF(ix3(i)/=ix4(i))THEN
402 pene(i) = max(p1(i),p2(i),p3(i),p4(i))
403 pene(i) = max(zero,pene(i)+gapv(i)-max(gapv(i)+dgaploadp,drad))
404C
405 la1 = one - lb1(i) - lc1(i)
406 la2 = one - lb2(i) - lc2(i)
407 la3 = one - lb3(i) - lc3(i)
408 la4 = one - lb4(i) - lc4(i)
409C
410 h0 = fourth *
411 . (p1(i)*la1 + p2(i)*la2 + p3(i)*la3
412 . + p4(i)*la4)
413 h1(i) = h0 + p1(i) * lb1(i) + p4(i) * lc4(i)
414 h2(i) = h0 + p2(i) * lb2(i) + p1(i) * lc1(i)
415 h3(i) = h0 + p3(i) * lb3(i) + p2(i) * lc2(i)
416 h4(i) = h0 + p4(i) * lb4(i) + p3(i) * lc3(i)
417 h00 = one/max(em20,h1(i) + h2(i) + h3(i) + h4(i))
418C
419 h1(i) = h1(i) * h00
420 h2(i) = h2(i) * h00
421 h3(i) = h3(i) * h00
422 h4(i) = h4(i) * h00
423C
424 ELSE
425 pene(i) = p1(i)
426 pene(i) = max(zero,pene(i)+gapv(i)-max(gapv(i)+dgaploadp,drad))
427 n1(i) = nx1(i)
428 n2(i) = ny1(i)
429 n3(i) = nz1(i)
430 h1(i) = lb1(i)
431 h2(i) = lc1(i)
432 h3(i) = one - lb1(i) - lc1(i)
433 h4(i) = zero
434 ENDIF
435 ENDDO
436 ENDIF ! IF(INTTH == 0)
437C---------------------
438C COURBURE FIXE
439C---------------------
440 IF(icurv(1)==1)THEN
441C spherical (only concave for now)
442 na1 = icurv(2)
443 DO i=1,jlt-jlt_tied
444 rr = 1.e30
445 a0x = x(1,na1)
446 a0y = x(2,na1)
447 a0z = x(3,na1)
448C
449 rx = x1(i)-a0x
450 ry = y1(i)-a0y
451 rz = z1(i)-a0z
452 rr = min(rr , rx*rx + ry*ry + rz*rz)
453 rx = x2(i)-a0x
454 ry = y2(i)-a0y
455 rz = z2(i)-a0z
456 rr = min(rr , rx*rx + ry*ry + rz*rz)
457 rx = x3(i)-a0x
458 ry = y3(i)-a0y
459 rz = z3(i)-a0z
460 rr = min(rr , rx*rx + ry*ry + rz*rz)
461 IF(ix3(i)/=ix4(i))THEN
462 rx = x4(i)-a0x
463 ry = y4(i)-a0y
464 rz = z4(i)-a0z
465 rr = min(rr , rx*rx + ry*ry + rz*rz)
466 ENDIF
467 rx = xi(i)-a0x
468 ry = yi(i)-a0y
469 rz = zi(i)-a0z
470 rs = sqrt(rx*rx + ry*ry + rz*rz)
471 rr = sqrt(rr)
472 IF(rs-rr+gapv(i)<0.)THEN
473 stif(i) = 0.
474 pene(i) = 0.
475 ELSEIF(rs-rr+gapv(i)<pene(i))THEN
476 pene(i) = rs-rr+gapv(i)
477 ENDIF
478 n1(i) = -rx
479 n2(i) = -ry
480 n3(i) = -rz
481 ENDDO
482 ELSEIF(icurv(1)==2)THEN
483 na1 = icurv(2)
484 na2 = icurv(3)
485 DO i=1,jlt-jlt_tied
486 rr = 1.e30
487 a0x = x(1,na1)
488 a0y = x(2,na1)
489 a0z = x(3,na1)
490 anx = x(1,na2)-a0x
491 any = x(2,na2)-a0y
492 anz = x(3,na2)-a0z
493 aan = 1. / (anx*anx + any*any + anz*anz)
494C
495 aax = x1(i)-a0x
496 aay = y1(i)-a0y
497 aaz = z1(i)-a0z
498 aaa = (aax*anx + aay*any + aaz*anz) * aan
499 rx = aax - aaa * anx
500 ry = aay - aaa * any
501 rz = aaz - aaa * anz
502 rr = min(rr , rx*rx + ry*ry + rz*rz)
503 aax = x2(i)-a0x
504 aay = y2(i)-a0y
505 aaz = z2(i)-a0z
506 aaa = (aax*anx + aay*any + aaz*anz) * aan
507 rx = aax - aaa * anx
508 ry = aay - aaa * any
509 rz = aaz - aaa * anz
510 rr = min(rr , rx*rx + ry*ry + rz*rz)
511 aax = x3(i)-a0x
512 aay = y3(i)-a0y
513 aaz = z3(i)-a0z
514 aaa = (aax*anx + aay*any + aaz*anz) * aan
515 rx = aax - aaa * anx
516 ry = aay - aaa * any
517 rz = aaz - aaa * anz
518 rr = min(rr , rx*rx + ry*ry + rz*rz)
519 IF(ix3(i)/=ix4(i))THEN
520 aax = x4(i)-a0x
521 aay = y4(i)-a0y
522 aaz = z4(i)-a0z
523 aaa = (aax*anx + aay*any + aaz*anz) * aan
524 rx = aax - aaa * anx
525 ry = aay - aaa * any
526 rz = aaz - aaa * anz
527 rr = min(rr , rx*rx + ry*ry + rz*rz)
528 ENDIF
529 aax = xi(i)-a0x
530 aay = yi(i)-a0y
531 aaz = zi(i)-a0z
532 aaa = (aax*anx + aay*any + aaz*anz) * aan
533 rx = aax - aaa * anx
534 ry = aay - aaa * any
535 rz = aaz - aaa * anz
536 rs = sqrt(rx*rx + ry*ry + rz*rz)
537 rr = sqrt(rr)
538C --
539 IF(rs-rr+gapv(i)<0.)THEN
540 stif(i) = 0.
541 pene(i) = 0.
542 ELSEIF(rs-rr+gapv(i)<pene(i))THEN
543 pene(i) = rs-rr+gapv(i)
544 n1(i) = -rx
545 n2(i) = -ry
546 n3(i) = -rz
547 ELSEIF(rs-rr-gapv(i)>0.)THEN
548 stif(i) = 0.
549 pene(i) = 0.
550 ELSEIF(rs-rr-gapv(i) < pene(i))THEN
551 xn1 = x1(i) - xi(i)
552 yn1 = y1(i) - yi(i)
553 zn1 = z1(i) - zi(i)
554 xn2 = x2(i) - xi(i)
555 yn2 = y2(i) - yi(i)
556 zn2 = z2(i) - zi(i)
557 xn3 = x3(i) - xi(i)
558 yn3 = y3(i) - yi(i)
559 zn3 = z3(i) - zi(i)
560C --
561 nn1 = (yn1*zn2-yn2*zn1) * rx +
562 . (zn1*xn2-zn2*xn1) * ry +
563 . (xn1*yn2-xn2*yn1) * rz
564 nn2 = (yn2*zn3-yn3*zn2) * rx +
565 . (zn2*xn3-zn3*xn2) * ry +
566 . (xn2*yn3-xn3*yn2) * rz
567 IF(ix3(i)/=ix4(i))THEN
568 xn4 = x4(i) - xi(i)
569 yn4 = y4(i) - yi(i)
570 zn4 = z4(i) - zi(i)
571 nn3 = (yn3*zn4-yn4*zn3) * rx +
572 . (zn3*xn4-zn4*xn3) * ry +
573 . (xn3*yn4-xn4*yn3) * rz
574 nn4 = (yn4*zn1-yn1*zn4) * rx +
575 . (zn4*xn1-zn1*xn4) * ry +
576 . (xn4*yn1-xn1*yn4) * rz
577 ELSE
578 nn3 = (yn3*zn1-yn1*zn3) * rx +
579 . (zn3*xn1-zn1*xn3) * ry +
580 . (xn3*yn1-xn1*yn3) * rz
581 nn4 = nn3
582 ENDIF
583 IF( nn1>=zero .AND. nn2>=zero
584 . .AND. nn3>=zero .AND. nn4>=zero) THEN
585 iproj = 1
586 ELSEIF( nn1<=zero .AND. nn2<=zero
587 . .AND. nn3<=zero .AND. nn4<=zero) THEN
588 iproj = 1
589 ELSE
590 iproj = 0
591 ENDIF
592C --
593 IF(iproj == 1)THEN
594 pene(i) = -rs+rr+gapv(i)
595 n1(i) = rx
596 n2(i) = ry
597 n3(i) = rz
598 ENDIF
599 ENDIF
600 ENDDO
601 ELSEIF(icurv(1) == 3)THEN
602 CALL i7curv(jlt ,pene ,n1 ,n2 ,
603 1 n3 ,gapv ,x ,nod_normal,
604 2 ix1 ,ix2 ,ix3 ,ix4 ,
605 3 h1 ,h2 ,h3 ,h4 ,
606 4 x1 ,x2 ,x3 ,x4 ,y1 ,
607 5 y2 ,y3 ,y4 ,z1 ,z2 ,
608 6 z3 ,z4 ,xi ,yi ,zi )
609 DO i=1,jlt-jlt_tied
610 IF(pene(i)<zero)THEN
611 stif(i) =zero
612 pene(i) =zero
613 END IF
614 END DO
615 ENDIF
616C
617 DO i=1,jlt-jlt_tied
618 s2 = one/max(em30,sqrt(n1(i)**2 + n2(i)**2 + n3(i)**2))
619 n1(i) = n1(i)*s2
620 n2(i) = n2(i)*s2
621 n3(i) = n3(i)*s2
622 ENDDO
623C
624 DO i=1,jlt-jlt_tied
625 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
626 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
627 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
628 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
629 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
630 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
631 vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
632C correction hourglass
633 h0 = -.25*(h1(i) - h2(i) + h3(i) - h4(i))
634 h0 = min(h0,h2(i),h4(i))
635 h0 = max(h0,-h1(i),-h3(i))
636 IF(ix3(i)==ix4(i))h0 = zero
637 h1(i) = h1(i) + h0
638 h2(i) = h2(i) - h0
639 h3(i) = h3(i) + h0
640 h4(i) = h4(i) - h0
641 END DO
642C-------------------------------------------
643C IMPACTS SUIVANTS (Tied)
644C------------------------------------
645 DO i=jlt-jlt_tied+1,jlt
646 ii=index(i)
647 h1(i) = cand_f(4,ii)
648 h2(i) = cand_f(5,ii)
649 h3(i) = cand_f(6,ii)
650 h4(i) = one - h1(i) - h2(i) - h3(i)
651 END DO
652C
653 DO i=jlt-jlt_tied+1,jlt
654C
655 d1 = sqrt(p1(i))
656 p1(i) = max(zero, gapv(i) - d1)
657C
658 d2 = sqrt(p2(i))
659 p2(i) = max(zero, gapv(i) - d2)
660C
661 d3 = sqrt(p3(i))
662 p3(i) = max(zero, gapv(i) - d3)
663C
664 d4 = sqrt(p4(i))
665 p4(i) = max(zero, gapv(i) - d4)
666C
667 IF(ix3(i)/=ix4(i))THEN
668 pene(i) = max(p1(i),p2(i),p3(i),p4(i))
669 ELSE
670 pene(i) = p1(i)
671 ENDIF
672C
673 END DO
674C------------------------------------
675 IF(itied /= 0)THEN
676C
677C Computed for every impact !
678 DO i=1,jlt
679C
680 t1x(i) = x3(i) - x1(i)
681 t1y(i) = y3(i) - y1(i)
682 t1z(i) = z3(i) - z1(i)
683 norminv = one/sqrt(t1x(i)**2+t1y(i)**2+t1z(i)**2)
684 t1x(i) = t1x(i)*norminv
685 t1y(i) = t1y(i)*norminv
686 t1z(i) = t1z(i)*norminv
687C
688 t2x(i) = x4(i) - x2(i)
689 t2y(i) = y4(i) - y2(i)
690 t2z(i) = z4(i) - z2(i)
691C
692 nx(i) = t1y(i)*t2z(i) - t1z(i)*t2y(i)
693 ny(i) = t1z(i)*t2x(i) - t1x(i)*t2z(i)
694 nz(i) = t1x(i)*t2y(i) - t1y(i)*t2x(i)
695 norminv = one/sqrt(nx(i)**2+ny(i)**2+nz(i)**2)
696 nx(i) = nx(i)*norminv
697 ny(i) = ny(i)*norminv
698 nz(i) = nz(i)*norminv
699C
700 t2x(i) = ny(i)*t1z(i) - nz(i)*t1y(i)
701 t2y(i) = nz(i)*t1x(i) - nx(i)*t1z(i)
702 t2z(i) = nx(i)*t1y(i) - ny(i)*t1x(i)
703C
704C SIDE(I)=NX(I)*N1(I)+NY(I)*N2(I)+NZ(I)*N3(I)
705 xh=h1(i)*x1(i)+h2(i)*x2(i)+h3(i)*x3(i)+h4(i)*x4(i)
706 yh=h1(i)*y1(i)+h2(i)*y2(i)+h3(i)*y3(i)+h4(i)*y4(i)
707 zh=h1(i)*z1(i)+h2(i)*z2(i)+h3(i)*z3(i)+h4(i)*z4(i)
708 side(i)=(xi(i)-xh)*nx(i)+(yi(i)-yh)*ny(i)+(zi(i)-zh)*nz(i)
709C
710 END DO
711 END IF
712C
713 DO i=jlt-jlt_tied+1,jlt
714 n1(i) = nx(i)
715 n2(i) = ny(i)
716 n3(i) = nz(i)
717 END DO
718C------------------------------------
719 DO i=jlt-jlt_tied+1,jlt
720 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
721 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
722 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
723 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
724 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
725 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
726 vn(i) = vx(i)*nx(i) + vy(i)*ny(i) + vz(i)*nz(i)
727 vt1(i) = vx(i)*t1x(i) + vy(i)*t1y(i) + vz(i)*t1z(i)
728 vt2(i) = vx(i)*t2x(i) + vy(i)*t2y(i) + vz(i)*t2z(i)
729 END DO
730C---------------------
731C PENE INITIALE
732C---------------------
733 IF(inacti==5)THEN
734 DO i=1,jlt-jlt_tied
735C reduction of the initial penetration
736C CAND_P(INDEX(I))=MIN(CAND_P(INDEX(I)),PENE(I))
737 cand_p(index(i))=min(cand_p(index(i)),
738 . ((one-fiveem2)*cand_p(index(i))+fiveem2*pene(i)) )
739C subtraction of the initial penetration from the penetration and the gap
740 pene(i)=max(zero,pene(i)-cand_p(index(i)))
741 IF( pene(i)==zero ) stif(i) = zero
742 gapv(i)=gapv(i)-cand_p(index(i))
743 ENDDO
744
745 DO i=jlt-jlt_tied+1,jlt
746C
747C CAND_P does not vary during tying (IF ITIED=1, rebound occurs when PENE < CAND_P)
748 pene(i)=max(zero,pene(i)-cand_p(index(i)))
749 gapv(i)=gapv(i)-cand_p(index(i))
750 ENDDO
751
752 ELSE IF(inacti==6)THEN
753 DO i=1,jlt-jlt_tied
754C reduction of the initial penetration
755C CAND_P(INDEX(I))=MIN(CAND_P(INDEX(I)),PENE(I))
756 cand_p(index(i))=min(cand_p(index(i)),
757 . ( (one-fiveem2)*cand_p(index(i))
758 . +fiveem2*(pene(i)+fiveem2*(gapv(i)-pene(i)))) )
759C subtraction of the initial penetration from the penetration and the gap
760 pene(i)=max(zero,pene(i)-cand_p(index(i)))
761 IF( pene(i)==zero ) stif(i) = zero
762 gapv(i)=gapv(i)-cand_p(index(i))
763 ENDDO
764
765 DO i=jlt-jlt_tied+1,jlt
766C
767C CAND_P does not vary during tying (IF ITIED=1, rebound occurs when PENE < CAND_P)
768 pene(i)=max(zero,pene(i)-cand_p(index(i)))
769 gapv(i)=gapv(i)-cand_p(index(i))
770 ENDDO
771
772 ENDIF
773C---------------------
774 IF(itied == 0)THEN
775C
776 dti = 1.e20
777C
778 DO 600 i=1,jlt
779 dist=gapv(i)-pene(i)
780 rdist = half*dist / max(em30,-vn(i))
781 dti = min(rdist,dti)
782 600 CONTINUE
783C
784C Force to DEL
785 IF (dtmini>zero) THEN
786 dtm=dtmini
787 idtm=2
788 ELSE
789 dtm=dtmin1(10)
790 idtm=idtmin(10)
791 END IF
792
793 IF(dti<=dtm)THEN
794 dti = 1.e20
795 DO i=1,jlt
796 dist =gapv(i)-pene(i)
797 dti2 = half*dist / max(em30,-vn(i))
798 IF(dti2<=dtm)THEN
799#include "lockon.inc"
800 WRITE(iout,'(A,E12.4,A,I10,A,E12.4,A)')
801 . ' **WARNING MINIMUM TIME STEP ',dti2,
802 . ' IN INTERFACE ',noint,'(dtmin=',DTM,')'
803 IF ( ISTAMPING == 1) THEN
804 WRITE(ISTDO,'(a)')'the run encountered a problem in an in
805 .terface Type 7.'
806 WRITE(ISTDO,'(a)')'you may need to check if there is enou
807 .gh clearance between the tools,'
808 WRITE(ISTDO,'(a)')'and that they do not penetrate each ot
809 .her during their travel'
810 WRITE(IOUT, '(a)')'the run encountered a problem in an in
811 .terface Type 7.'
812 WRITE(IOUT, '(a)')'you may need to check if there is enou
813 .gh clearance between the tools,'
814 WRITE(IOUT, '(a)')'and that they do not penetrate each ot
815 .her during their travel'
816 ENDIF
817 NN = NSVG(I)
818 IF(NN>0)THEN
819 NI = ITAB(NN)
820 ELSE
821 NI = ITAFI(NIN)%P(-NN)
822 ENDIF
823#include "lockoff.inc"
824 IF(IDTM==1)THEN
825#include "lockon.inc"
826 WRITE(IOUT,'(a,i10)') ' secondary node : ',NI
827 WRITE(IOUT,'(a,4i10)')' main nodes : ',
828 . ITAB(IX1(I)),ITAB(IX2(I)),ITAB(IX3(I)),ITAB(IX4(I))
829#include "lockoff.inc"
830 TSTOP = TT
831 ELSEIF(IDTM==2)THEN
832#include "lockon.inc"
833 WRITE(IOUT,'(a,i10,a,i10)')' remove secondary node ',
834 . NI,' from INTERFACE ',NOINT
835 IF(NN>0) THEN
836 STFN(CN_LOC(I)) = -ABS(STFN(CN_LOC(I)))
837 ELSE
838 STIFI(NIN)%P(-NN) = -ABS(STIFI(NIN)%P(-NN))
839 ENDIF
840#include "lockoff.inc"
841 STIF(I) = ZERO
842 NEWFRONT = -1
843 DTI = DTM
844
845 ELSEIF(IDTM==5)THEN
846#include "lockon.inc"
847 WRITE(IOUT,'(a,i10)') ' secondary node : ',NI
848 WRITE(IOUT,'(a,4i10)')' main nodes : ',
849 . ITAB(IX1(I)),ITAB(IX2(I)),ITAB(IX3(I)),ITAB(IX4(I))
850#include "lockoff.inc"
851 MSTOP = 2
852.AND. ELSEIF(IDTM==6ILAGM==2)THEN
853 IG=NSVG(I)
854 IF(KINET(IG)+KINET(IX1(I))+KINET(IX2(I))
855 . +KINET(IX3(I))+KINET(IX4(I))==0)THEN
856 CAND_N(INDEX(I)) = -IABS(CAND_N(INDEX(I)))
857 STIF(I) = ZERO
858 DTI2 = 1.E20
859#include "lockon.inc"
860 WRITE(IOUT,'(a,i10)') ' secondary node : ',ITAB(NSVG(I))
861 WRITE(IOUT,'(a,4i10)')' main nodes : ',
862 . ITAB(IX1(I)),ITAB(IX2(I)),ITAB(IX3(I)),ITAB(IX4(I))
863#include "lockoff.inc"
864 ENDIF
865 DTI = MIN(DTI2,DTI)
866 ENDIF
867 ENDIF
868 ENDDO
869 ENDIF
870 IF(DTI<DT2T)THEN
871 DT2T = DTI
872 NELTST = NOINT
873 ITYPTST = 10
874 ENDIF
875
876 ENDIF ! IF(ITIED == 0)THEN
877C-------------------------------------------
878C
879C 2:const,1:linear,0:nonlinear
880 IF(IMPL_S>0)THEN
881 IF(IMP_INT7==2)THEN
882 DO I=1,JLT
883 IF(STIGLO<=ZERO)THEN
884 STIF(I) = HALF*STIF(I)
885 ELSEIF(STIF(I)/=ZERO)THEN
886 STIF(I) = STIGLO
887 ENDIF
888 FNI(I)= -STIF(I) * PENE(I)
889 ENDDO
890 ELSE
891 DO I=1,JLT
892 FAC = GAPV(I)/MAX( EM10,( GAPV(I)-PENE(I) ) )
893 FACM1 = 1./FAC
894.AND. IF( (GAPV(I)-PENE(I))/GAPV(I) <PREC
895 . STIF(I)>0. ) THEN
896 STIF(I) = 0.
897 NEWFRONT = -1
898#include "lockon.inc"
899 NN = NSVG(I)
900 IF(NN>0)THEN
901 NI = ITAB(NN)
902 STFN(CN_LOC(I)) = -ABS(STFN(CN_LOC(I)))
903 ELSE
904 NI = ITAFI(NIN)%P(-NN)
905 STIFI(NIN)%P(-NN) = -ABS(STIFI(NIN)%P(-NN))
906 ENDIF
907 WRITE(ISTDO,'(a,i10)')' warning INTERFACE ',NOINT
908 WRITE(ISTDO,'(a,i10,a)')' node ',NI,
909 . ' de-activated from interface'
910 WRITE(IOUT ,'(a,i10)')' warning INTERFACE ',NOINT
911 WRITE(IOUT ,'(a,i10,a)')' node ',NI,
912 . ' de-activated from interface'
913 IF ( ISTAMPING == 1) THEN
914 WRITE(ISTDO,'(a)')'the run encountered a problem in an inter
915 .face Type 7.'
916 WRITE(ISTDO,'(a)')'you may need to check if there is enough
917 .clearance between the tools,'
918 WRITE(ISTDO,'(a)')'and that they do not penetrate each other
919 . during their travel'
920 WRITE(IOUT, '(a)')'the run encountered a problem in an inter
921 .face Type 7.'
922 WRITE(IOUT, '(a)')'you may need to check if there is enough
923 .clearance between the tools,'
924 WRITE(IOUT, '(a)')'and that they do not penetrate each other
925 . during their travel'
926 ENDIF
927#include "lockoff.inc"
928 ENDIF
929 IF(STIGLO<=ZERO)THEN
930 IF(INCONV==1) ECONTT = ECONTT + HALF*STIF(I)*GAPV(I)**2 *
931 . ( FACM1 -ONE -LOG(MAX(TINY(FACM1),FACM1) ))
932 STIF(I) = HALF*STIF(I) * FAC
933 ELSEIF(STIF(I)/=ZERO)THEN
934 IF(INCONV==1)ECONTT = ECONTT + STIGLO*GAPV(I)**2 *
935 . ( FACM1 - ONE -LOG(MAX(TINY(FACM1),FACM1) ))
936 STIF(I) = STIGLO * FAC
937 ENDIF
938 FNI(I)= -STIF(I) * PENE(I)
939 ENDDO
940 ENDIF
941 IF(STIGLO<=ZERO)THEN
942 DO I=1,JLT
943 STIF0(I) = HALF*STIF0(I)
944 ENDDO
945 ENDIF
946 ELSE ! fin impl_s>0
947 DO I=1,JLT-JLT_TIED
948 FAC = GAPV(I)/MAX( EM10,( GAPV(I)-PENE(I) ) )
949 FACM1 = ONE/FAC
950.AND. IF( (GAPV(I)-PENE(I))/GAPV(I) <PREC
951 . STIF(I)>ZERO ) THEN
952 STIF(I) = ZERO
953 NEWFRONT = -1
954#include "lockon.inc"
955 NN = NSVG(I)
956 IF(NN>0)THEN
957 NI = ITAB(NN)
958 STFN(CN_LOC(I)) = -ABS(STFN(CN_LOC(I)))
959 ELSE
960 NI = ITAFI(NIN)%P(-NN)
961 STIFI(NIN)%P(-NN) = -ABS(STIFI(NIN)%P(-NN))
962 ENDIF
963 WRITE(ISTDO,'(a,i10)')' warning INTERFACE ',NOINT
964 WRITE(ISTDO,'(a,i10,a)')' node ',NI,
965 . ' de-activated from interface'
966 WRITE(IOUT ,'(a,i10)')' warning INTERFACE ',NOINT
967 WRITE(IOUT ,'(a,i10,a)')' NODE ',ni,
968 . ' DE-ACTIVATED FROM INTERFACE'
969 IF ( istamping == 1) THEN
970 WRITE(istdo,'(A)')'The run encountered a problem in an inter
971 .face Type 7.'
972 WRITE(istdo,'(A)')'You may need to check if there is enough
973 .clearance between the tools,'
974 WRITE(istdo,'(A)')'and that they do not penetrate each other
975 . during their travel'
976 WRITE(iout, '(A)')'The run encountered a problem in an inter
977 .face Type 7.'
978 WRITE(iout, '(A)')'You may need to check if there is enough
979 .clearance between the tools,'
980 WRITE(iout, '(A)')'and that they do not penetrate each other
981 . during their travel'
982 ENDIF
983#include "lockoff.inc"
984 ENDIF
985 IF(stiglo<=zero)THEN
986 econtt = econtt + half*stif(i)*gapv(i)**2 *( facm1 - one -
987 . log(max(tiny(facm1),facm1)) )
988 stif(i) = half*stif(i) * fac
989 ELSEIF(stif(i)/=zero)THEN
990 econtt = econtt + stiglo*gapv(i)**2 *( facm1 - one -
991 . log(max(tiny(facm1),facm1)) )
992 stif(i) = stiglo * fac
993 ENDIF
994 fni(i)= -stif(i) * pene(i)
995 END DO
996 ENDIF
997C---------------------------------
998 IF(itied/=0)THEN
999 DO i=1,jlt-jlt_tied
1000 fns(i)= fni(i)
1001 END DO
1002 END IF
1003C---------------------------------
1004C Force vs tied impacts (cf TYPE10)
1005C---------------------------------
1006 IF(itied == 1)THEN
1007 DO i=jlt-jlt_tied+1,jlt
1008
1009 ii=index(i)
1010 IF(pene(i)==zero.AND.side(i)*cand_f(8,ii) > zero)THEN
1011C--------------------------------------
1012C REBOUND
1013C--------------------------------------
1014 vn(i) = zero
1015 vt1(i) = zero
1016 vt2(i) = zero
1017 END IF
1018 END DO
1019 END IF
1020
1021 dt05=half*dt1
1022 DO i=jlt-jlt_tied+1,jlt
1023 ii=index(i)
1024
1025 fac = gapv(i)/max( em10,( gapv(i)-min(pene(i),cand_f(7,ii)) ) )
1026 stif(i) = half*stif0(i) * fac**2 ! Tangent Stiffness
1027
1028 econttied = econttied + cand_f(1,ii) * vn(i) * dt05
1029 econttied = econttied + cand_f(2,ii) * vt1(i) * dt05
1030 econttied = econttied + cand_f(3,ii) * vt2(i) * dt05
1031
1032 fni(i) = cand_f(1,ii) + vn(i) * dt1 * stif(i)
1033 ft1(i) = cand_f(2,ii) + vt1(i) * dt1 * stif(i)
1034 ft2(i) = cand_f(3,ii) + vt2(i) * dt1 * stif(i)
1035 END DO
1036
1037 IF(itied == 1)THEN
1038 DO i=jlt-jlt_tied+1,jlt
1039
1040 ii=index(i)
1041 IF(pene(i)==zero.AND.side(i)*cand_f(8,ii) > zero)THEN
1042C--------------------------------------
1043C REBOUND
1044C--------------------------------------
1045 cand_f(1,ii) =zero
1046 fni(i) = zero
1047 ft1(i) = zero
1048 ft2(i) = zero
1049 vn(i) = zero
1050 vt1(i) = zero
1051 vt2(i) = zero
1052 ELSEIF(fni(i)==zero)THEN
1053 cand_f(1,ii) = sign(em30,cand_f(1,ii))
1054 ELSE
1055 IF (inconv==1) cand_f(1,ii) =fni(i) ! CAND_F(1,II) negative vs NX, NY, NZ
1056 ENDIF
1057 IF (inconv==1) THEN
1058 cand_f(2,ii) = ft1(i)
1059 cand_f(3,ii) = ft2(i)
1060 ENDIF
1061 END DO
1062 ELSEIF(itied == 2)THEN
1063 DO i=jlt-jlt_tied+1,jlt
1064
1065 ii=index(i)
1066 IF(fni(i)==zero.AND.pene(i)/=zero)THEN
1067 cand_f(1,ii) = em30
1068 ELSE
1069 IF (inconv==1) cand_f(1,ii) =fni(i) ! CAND_F(1,II) negative vs NX, NY, NZ
1070 ENDIF
1071 IF (inconv==1) THEN
1072 cand_f(2,ii) = ft1(i)
1073 cand_f(3,ii) = ft2(i)
1074 ENDIF
1075 END DO
1076 END IF
1077C
1078 DO i=jlt-jlt_tied+1,jlt
1079 ii = index(i)
1080 econttied = econttied + cand_f(1,ii) * vn(i) * dt05
1081 econttied = econttied + cand_f(2,ii) * vt1(i) * dt05
1082 econttied = econttied + cand_f(3,ii) * vt2(i) * dt05
1083 END DO
1084C---------------------------------
1085C DAMPING + FRIC
1086C---------------------------------
1087C goto 999
1088
1089 IF(visc/=zero)THEN
1090 IF(ivis2==0.OR.ivis2==1)THEN
1091C---------------------------------
1092C VISC QUAD TYPE V227
1093C---------------------------------
1094 DO i=1,jlt-jlt_tied
1095 vis2(i) = two * stif(i) * msi(i)
1096 IF(vn(i)<zero) vis2(i) = vis2(i) /
1097 . ( max(em10,(gapv(i)-pene(i))/gapv(i)) )
1098 ENDDO
1099C---------------------------------
1100 visca = zep4
1101 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
1102 DO i=1,jlt-jlt_tied
1103 fac = stif(i) / max(em30,stif(i))
1104 vis = sqrt(vis2(i))
1105 ff = fac * (
1106 . visc * vis +
1107 . visca**2 * two* msi(i) * max(zero,-vn(i)) /
1108 . max((gapv(i) - pene(i)),em10) )
1109 stif(i) = stif(i) * gapv(i) / max((gapv(i) - pene(i)),em10)
1110 stif(i) = stif(i) + ff * dt1inv
1111 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
1112 ffo = ff
1113 ff = ff * vn(i)
1114 fni(i) = fni(i) + ff
1115 ENDDO
1116 ELSE
1117 DO i=1,jlt-jlt_tied
1118 fac = stif(i) / max(em30,stif(i))
1119 vis = sqrt(vis2(i))
1120 c(i)= fac * (
1121 . visc * vis +
1122 . visca**2 * two * msi(i) * max(zero,-vn(i)) /
1123 . max((gapv(i) - pene(i)),em10) )
1124 stif(i) = stif(i) * gapv(i) / max((gapv(i) - pene(i)),em10)
1125 kt(i)= stif(i)
1126 stif(i) = stif(i) + c(i) * dt1inv
1127 ff = c(i) * vn(i)
1128 fni(i) = fni(i) + ff
1129 cf(i) = fac*sqrt(viscffric(i))*vis
1130 stif(i) = max(stif(i) ,cf(i)*dt1inv)
1131 ENDDO
1132 ENDIF
1133 ELSEIF(ivis2==2)THEN
1134C---------------------------------
1135C VISC QUAD TYPE
1136C---------------------------------
1137 DO i=1,jlt-jlt_tied
1138 vis2(i) = two* stif(i) * msi(i)
1139 vis2(i) = vis2(i) /
1140 . ( max(em10,(gapv(i)-pene(i))/gapv(i)) )
1141 ENDDO
1142C---------------------------------
1143 visca = half
1144 DO i=1,jlt-jlt_tied
1145 fac = stif(i) / max(em30,stif(i))
1146 vis = sqrt(vis2(i))
1147 ff = fac * (
1148 . visc * vis +
1149 . visca**2 * two * msi(i) * abs(vn(i)) /
1150 . max((gapv(i) - pene(i)),em10) )
1151 stif(i) = stif(i) * gapv(i) / max((gapv(i) - pene(i)),em10)
1152 stif(i) = stif(i) + two * ff * dt1inv
1153 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
1154 ff = ff * vn(i)
1155 fni(i) = fni(i) + ff
1156 ENDDO
1157 ELSEIF(ivis2==3)THEN
1158C---------------------------------
1159C VISC QUAD = 0
1160C---------------------------------
1161 DO i=1,jlt-jlt_tied
1162 vis2(i) = two * stif(i) * msi(i)
1163 ENDDO
1164C---------------------------------
1165 DO i=1,jlt-jlt_tied
1166 fac = stif(i) / max(em30,stif(i))
1167 vis = sqrt(vis2(i))
1168 ff = fac * ( visc * vis ) /
1169 . max((gapv(i) - pene(i)),em10)
1170 stif(i) = stif(i) * gapv(i) /
1171 . max((gapv(i) - pene(i)),em10)
1172 stif(i) = stif(i) + two* ff * dt1inv
1173 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
1174 ff = ff * vn(i)
1175 fni(i) = fni(i) + ff
1176 ENDDO
1177 ELSEIF(ivis2==4)THEN
1178C---------------------------------
1179C VISC = 0
1180C---------------------------------
1181 DO i=1,jlt-jlt_tied
1182 vis2(i) = two* stif(i) * msi(i)
1183 vis = sqrt(vis2(i))
1184 stif(i) = stif(i) * gapv(i) /
1185 . max((gapv(i) - pene(i)),em10)
1186 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
1187 ENDDO
1188 ELSEIF(ivis2==5)THEN
1189C---------------------------------
1190C Visc = 2m/dt => For visc <1, stable: Dt <2m/visc = Dt
1191C M = m1*m2/m1+m2 for visc = 1, elastic shock
1192C For visc = 0.5, elastic collision
1193C---------------------------------
1194 DO i=1,jlt-jlt_tied
1195 mas2 = ms(ix1(i))*h1(i)
1196 . + ms(ix2(i))*h2(i)
1197 . + ms(ix3(i))*h3(i)
1198 . + ms(ix4(i))*h4(i)
1199 vis2(i) = two* stif(i) * msi(i)
1200 vis = 2. * visc * dt1inv * msi(i) * mas2 /
1201 . max(em30,msi(i)+mas2)
1202 stif(i) = stif(i) * gapv(i) /
1203 . max((gapv(i) - pene(i)),em10)
1204 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i)*vis2(i))*dt1inv)
1205 ff = vis * vn(i)
1206 econtdt = econtdt + min(zero,ff-fni(i)) * vn(i) * dt1
1207 fni(i) = min(fni(i),ff)
1208 ENDDO
1209 ELSE
1210 ENDIF
1211 ELSE
1212 DO i=1,jlt-jlt_tied
1213 IF(viscffric(i)/=zero) THEN
1214 IF(ivis2==0.OR.ivis2==1)THEN
1215C---------------------------------
1216C VISC QUAD TYPE V227
1217C---------------------------------
1218 vis2(i) = two * stif(i) * msi(i)
1219 IF(vn(i)<zero) vis2(i) = vis2(i) /
1220 . ( max(em10,(gapv(i)-pene(i))/gapv(i)) )
1221C---------------------------------
1222C---------------------------------
1223 visca = zep4
1224 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
1225 fac = stif(i) / max(em30,stif(i))
1226 vis = sqrt(vis2(i))
1227 ff = fac * (
1228 . visc * vis +
1229 . visca**2 * two* msi(i) * max(zero,-vn(i)) /
1230 . max((gapv(i) - pene(i)),em10) )
1231 stif(i) = stif(i) * gapv(i) / max((gapv(i) - pene(i)),em10)
1232 stif(i) = stif(i) + ff * dt1inv
1233!
1234 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
1235 ffo = ff
1236 ff = ff * vn(i)
1237 fni(i) = fni(i) + ff
1238 ELSE
1239 fac = stif(i) / max(em30,stif(i))
1240 vis = sqrt(vis2(i))
1241 c(i)= fac * (
1242 . visc * vis +
1243 . visca**2 * two * msi(i) * max(zero,-vn(i)) /
1244 . max((gapv(i) - pene(i)),em10) )
1245 stif(i) = stif(i) * gapv(i) / max((gapv(i) - pene(i)),em10)
1246 kt(i)= stif(i)
1247 stif(i) = stif(i) + c(i) * dt1inv
1248 ff = c(i) * vn(i)
1249 fni(i) = fni(i) + ff
1250 cf(i) = fac*sqrt(viscffric(i))*vis
1251 stif(i) = max(stif(i) ,cf(i)*dt1inv)
1252 ENDIF
1253 ELSEIF(ivis2==2)THEN
1254C---------------------------------
1255C VISC QUAD TYPE
1256C---------------------------------
1257 vis2(i) = two* stif(i) * msi(i)
1258 vis2(i) = vis2(i) /
1259 . ( max(em10,(gapv(i)-pene(i))/gapv(i)) )
1260C---------------------------------
1261 visca = half
1262 fac = stif(i) / max(em30,stif(i))
1263 vis = sqrt(vis2(i))
1264 ff = fac * (
1265 . visc * vis +
1266 . visca**2 * two * msi(i) * abs(vn(i)) /
1267 . max((gapv(i) - pene(i)),em10) )
1268 stif(i) = stif(i) * gapv(i) / max((gapv(i) - pene(i)),em10)
1269 stif(i) = stif(i) + two * ff * dt1inv
1270 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
1271 ff = ff * vn(i)
1272 fni(i) = fni(i) + ff
1273 ELSEIF(ivis2==3)THEN
1274C---------------------------------
1275C VISC QUAD = 0
1276C---------------------------------
1277 vis2(i) = two * stif(i) * msi(i)
1278C---------------------------------
1279 fac = stif(i) / max(em30,stif(i))
1280 vis = sqrt(vis2(i))
1281 ff = fac * ( visc * vis ) /
1282 . max((gapv(i) - pene(i)),em10)
1283 stif(i) = stif(i) * gapv(i) /
1284 . max((gapv(i) - pene(i)),em10)
1285 stif(i) = stif(i) + two* ff * dt1inv
1286 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
1287 ff = ff * vn(i)
1288 fni(i) = fni(i) + ff
1289 ELSEIF(ivis2==4)THEN
1290C---------------------------------
1291C VISC = 0
1292C---------------------------------
1293 vis2(i) = two* stif(i) * msi(i)
1294 vis = sqrt(vis2(i))
1295 stif(i) = stif(i) * gapv(i) /
1296 . max((gapv(i) - pene(i)),em10)
1297 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
1298 ELSEIF(ivis2==5)THEN
1299C---------------------------------
1300C Visc = 2m/dt => For visc <1, stable: Dt <2m/visc = Dt
1301C M = m1*m2/m1+m2 for visc = 1, elastic shock
1302C For visc = 0.5, elastic collision
1303C---------------------------------
1304 mas2 = ms(ix1(i))*h1(i)
1305 . + ms(ix2(i))*h2(i)
1306 . + ms(ix3(i))*h3(i)
1307 . + ms(ix4(i))*h4(i)
1308 vis2(i) = two* stif(i) * msi(i)
1309 vis = 2. * visc * dt1inv * msi(i) * mas2 /
1310 . max(em30,msi(i)+mas2)
1311 stif(i) = stif(i) * gapv(i) /
1312 . max((gapv(i) - pene(i)),em10)
1313 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i)*vis2(i))*dt1inv)
1314 ff = vis * vn(i)
1315 econtdt = econtdt + min(zero,ff-fni(i)) * vn(i) * dt1
1316 fni(i) = min(fni(i),ff)
1317 ELSE
1318 ENDIF
1319 ELSE
1320cbb initialization of the VIS2 array to avoid problems
1321 vis2(i) = zero
1322 stif(i) = stif(i) * gapv(i) /
1323 . max((gapv(i) - pene(i)),em10)
1324 ENDIF
1325 ENDDO
1326 ENDIF
1327C---------------------------------
1328C Damping vs tied impacts
1329C---------------------------------
1330 IF(itied/=0)THEN
1331 IF(visc/=zero)THEN
1332 IF(ivis2==0.OR.ivis2==1)THEN
1333C---------------------------------
1334C VISC QUAD TYPE V227
1335C---------------------------------
1336 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
1337 DO i=jlt-jlt_tied+1,jlt
1338 fac = gapv(i)/max( em10,( gapv(i)-min(pene(i),cand_f(7,index(i))) ) )
1339 stif(i) = half*stif0(i) * fac
1340
1341 vis = visc * sqrt(two * stif(i) * msi(i))
1342 fni(i) = fni(i) + vn(i) * vis
1343 ft1(i) = ft1(i) + vt1(i) * vis
1344 ft2(i) = ft2(i) + vt2(i) * vis
1345
1346 stif(i) = stif(i)*fac
1347 stif(i) = stif(i) + two * vis * dt1inv
1348
1349 econtdt = econtdt
1350 . + vis * (vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)) * dt1
1351 END DO
1352
1353 ELSE
1354 DO i=jlt-jlt_tied+1,jlt
1355 fac = gapv(i)/max( em10,( gapv(i)-min(pene(i),cand_f(7,index(i))) ) )
1356 stif(i) = half*stif0(i) * fac
1357
1358 vis = visc * sqrt(two * stif(i) * msi(i))
1359 fni(i) = fni(i) + vn(i) * vis
1360 ft1(i) = ft1(i) + vt1(i) * vis
1361 ft2(i) = ft2(i) + vt2(i) * vis
1362 c(i) = vis
1363
1364 stif(i) = stif(i)*fac
1365 kt(i) = stif(i)
1366 stif(i) = stif(i) + two * vis * dt1inv
1367
1368 econtdt = econtdt
1369 . + vis * (vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)) * dt1
1370 END DO
1371 ENDIF
1372 ELSEIF(ivis2==2)THEN
1373C---------------------------------
1374C VISC QUAD TYPE
1375C---------------------------------
1376 DO i=jlt-jlt_tied+1,jlt
1377 fac = gapv(i)/max( em10,( gapv(i)-min(pene(i),cand_f(7,index(i))) ) )
1378 stif(i) = half*stif0(i) * fac
1379
1380 vis = sqrt(two* stif(i) * msi(i))
1381 fni(i) = fni(i) + vn(i) * vis
1382 ft1(i) = ft1(i) + vt1(i) * vis
1383 ft2(i) = ft2(i) + vt2(i) * vis
1384
1385 stif(i) = stif(i)*fac
1386 stif(i) = stif(i) + two * vis * dt1inv
1387
1388 econtdt = econtdt
1389 . + vis * (vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)) * dt1
1390 END DO
1391 ELSEIF(ivis2==3)THEN
1392C---------------------------------
1393C VISC QUAD = 0
1394C---------------------------------
1395 DO i=jlt-jlt_tied+1,jlt
1396 fac = gapv(i)/max( em10,( gapv(i)-min(pene(i),cand_f(7,index(i))) ) )
1397 stif(i) = half*stif0(i) * fac
1398
1399 vis = sqrt(two* stif(i) * msi(i))
1400 fni(i) = fni(i) + vn(i) * vis
1401 ft1(i) = ft1(i) + vt1(i) * vis
1402 ft2(i) = ft2(i) + vt2(i) * vis
1403
1404 stif(i) = stif(i)*fac
1405 stif(i) = stif(i) + two * vis * dt1inv
1406
1407 econtdt = econtdt
1408 . + vis * (vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)) * dt1
1409 END DO
1410 ELSEIF(ivis2==4)THEN
1411C---------------------------------
1412C VISC = 0
1413C---------------------------------
1414 ELSEIF(ivis2==5)THEN
1415C---------------------------------
1416C Visc = 2m/dt => For visc <1, stable: Dt <2m/visc = Dt
1417C M = m1*m2/m1+m2 for visc = 1, elastic shock
1418C For visc = 0.5, elastic collision
1419C---------------------------------
1420 DO i=jlt-jlt_tied+1,jlt
1421 fac = gapv(i)/max( em10,( gapv(i)-min(pene(i),cand_f(7,index(i))) ) )
1422 stif(i) = half*stif0(i) * fac
1423
1424 mas2 = ms(ix1(i))*h1(i)
1425 . + ms(ix2(i))*h2(i)
1426 . + ms(ix3(i))*h3(i)
1427 . + ms(ix4(i))*h4(i)
1428 vis = sqrt(two* stif(i) * msi(i))
1429 fni(i) = fni(i) + vn(i) * vis
1430 ft1(i) = ft1(i) + vt1(i) * vis
1431 ft2(i) = ft2(i) + vt2(i) * vis
1432
1433 stif(i) = stif(i)*fac
1434 stif(i) = stif(i) + two * vis * dt1inv
1435
1436 econtdt = econtdt
1437 . + vis * (vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)) * dt1
1438 END DO
1439 ELSE
1440 DO i=jlt-jlt_tied+1,jlt
1441 fac = gapv(i)/max( em10,( gapv(i)-min(pene(i),cand_f(7,index(i))) ) )
1442 stif(i) = half*stif0(i) * fac
1443 stif(i) = stif(i)*fac
1444 END DO
1445 ENDIF
1446 ELSE
1447 DO i=jlt-jlt_tied+1,jlt
1448 fac = gapv(i)/max( em10,( gapv(i)-min(pene(i),cand_f(7,index(i))) ) )
1449 stif(i) = half*stif0(i) * fac
1450 stif(i) = stif(i)*fac
1451 END DO
1452 END IF
1453 END IF
1454C---------------------------------
1455C SAUVEGARDE DE L'IMPULSION NORMALE
1456C---------------------------------
1457
1458 fsav1 = zero
1459 fsav2 = zero
1460 fsav3 = zero
1461 fsav8 = zero
1462 fsav9 = zero
1463 fsav10= zero
1464 fsav11= zero
1465C
1466 DO i=1,jlt
1467 fxi(i)=n1(i)*fni(i)
1468 fyi(i)=n2(i)*fni(i)
1469 fzi(i)=n3(i)*fni(i)
1470C if type7 sym of type19 opposite sign for normal force -> forc_sign = -1
1471 impx=forc_sign*fxi(i)*dt12
1472 impy=forc_sign*fyi(i)*dt12
1473 impz=forc_sign*fzi(i)*dt12
1474 fsav1 =fsav1 +impx
1475 fsav2 =fsav2 +impy
1476 fsav3 =fsav3 +impz
1477 fsav8 =fsav8 +abs(impx)
1478 fsav9 =fsav9 +abs(impy)
1479 fsav10=fsav10+abs(impz)
1480 fsav11=fsav11+abs(fni(i))*dt12
1481 ENDDO
1482#include "lockon.inc"
1483 fsav(1)=fsav(1)+fsav1
1484 fsav(2)=fsav(2)+fsav2
1485 fsav(3)=fsav(3)+fsav3
1486 fsav(8)=fsav(8)+fsav8
1487 fsav(9)=fsav(9)+fsav9
1488 fsav(10)=fsav(10)+fsav10
1489 fsav(11)=fsav(11)+fsav11
1490C
1491#include "lockoff.inc"
1492C
1493 IF(isensint(1)/=0) THEN
1494 DO i=1,jlt
1495 fsavparit(1,1,i+nft) = forc_sign*fxi(i)
1496 fsavparit(1,2,i+nft) = forc_sign*fyi(i)
1497 fsavparit(1,3,i+nft) = forc_sign*fzi(i)
1498 ENDDO
1499 ENDIF
1500C---------------------------------
1501 IF((anim_v(12)+outp_v(12)+h3d_data%N_VECT_PCONT>0.AND.
1502 . ((tt>=output%TANIM .AND. tt<=output%TANIM_STOP).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
1503 . (manim>=4.AND.manim<=15).OR.h3d_data%MH3D/=0))
1504 . .OR.h3d_data%N_VECT_PCONT_MAX>0)THEN
1505 IF (inconv==1) THEN
1506#include "lockon.inc"
1507 DO i=1,jlt
1508 fncont(1,ix1(i)) =fncont(1,ix1(i)) + fxi(i)*h1(i)
1509 fncont(2,ix1(i)) =fncont(2,ix1(i)) + fyi(i)*h1(i)
1510 fncont(3,ix1(i)) =fncont(3,ix1(i)) + fzi(i)*h1(i)
1511 fncont(1,ix2(i)) =fncont(1,ix2(i)) + fxi(i)*h2(i)
1512 fncont(2,ix2(i)) =fncont(2,ix2(i)) + fyi(i)*h2(i)
1513 fncont(3,ix2(i)) =fncont(3,ix2(i)) + fzi(i)*h2(i)
1514 fncont(1,ix3(i)) =fncont(1,ix3(i)) + fxi(i)*h3(i)
1515 fncont(2,ix3(i)) =fncont(2,ix3(i)) + fyi(i)*h3(i)
1516 fncont(3,ix3(i)) =fncont(3,ix3(i)) + fzi(i)*h3(i)
1517 fncont(1,ix4(i)) =fncont(1,ix4(i)) + fxi(i)*h4(i)
1518 fncont(2,ix4(i)) =fncont(2,ix4(i)) + fyi(i)*h4(i)
1519 fncont(3,ix4(i)) =fncont(3,ix4(i)) + fzi(i)*h4(i)
1520 jg = nsvg(i)
1521 IF(jg>0) THEN
1522C In SPMD: Treatment to be redone after reception node Remote if JG <0
1523 fncont(1,jg)=fncont(1,jg)- fxi(i)
1524 fncont(2,jg)=fncont(2,jg)- fyi(i)
1525 fncont(3,jg)=fncont(3,jg)- fzi(i)
1526 ELSE ! cas noeud remote en SPMD
1527 jg = -jg
1528 fnconti(nin)%P(1,jg)=fnconti(nin)%P(1,jg)-fxi(i)
1529 fnconti(nin)%P(2,jg)=fnconti(nin)%P(2,jg)-fyi(i)
1530 fnconti(nin)%P(3,jg)=fnconti(nin)%P(3,jg)-fzi(i)
1531 ENDIF
1532 ENDDO
1533#include "lockoff.inc"
1534 END IF !(INCONV==1) THEN
1535 ENDIF
1536C
1537C---------------------------------
1538C SORTIES TH PAR SOUS INTERFACE
1539C---------------------------------
1540
1541 IF(nisub/=0)THEN
1542 DO jsub=1,nisub
1543 DO j=1,24
1544 fsavsub1(j,jsub)=zero
1545 END DO
1546 ENDDO
1547 DO i=1,jlt
1548 nn = nsvg(i)
1549 IF(nn>0)THEN
1550 in=cn_loc(i)
1551 ie=ce_loc(i)
1552 jj =addsubs(in)
1553 kk =addsubm(ie)
1554 DO WHILE(jj<addsubs(in+1))
1555 jsub=lisubs(jj)
1556 itypsub = typsub(jsub)
1557
1558 IF(itypsub == 1 ) THEN ! Defining specific inter
1559
1560 ksub=lisubm(kk)
1561
1562
1563 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1564
1565 IF(ksub==jsub)THEN
1566C if type7 sym of type19 opposite sign for normal force -> forc_sign = -1
1567 impx=forc_sign*fxi(i)*dt12
1568 impy=forc_sign*fyi(i)*dt12
1569 impz=forc_sign*fzi(i)*dt12
1570C MAIN side :
1571 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1572 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1573 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1574C
1575 IF(isensint(jsub+1)/=0) THEN
1576 fsavparit(jsub+1,1,i+nft) = forc_sign*fxi(i)
1577 fsavparit(jsub+1,2,i+nft) = forc_sign*fyi(i)
1578 fsavparit(jsub+1,3,i+nft) = forc_sign*fzi(i)
1579 ENDIF
1580C
1581 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1582 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1583 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1584C
1585 fsavsub1(11,jsub)=fsavsub1(11,jsub)+abs(fni(i))*dt12
1586C
1587 END IF
1588
1589 kk=kk+1
1590 ksub=lisubm(kk)
1591
1592 ENDDO
1593 jj=jj+1
1594
1595 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface : secondary side
1596
1597 impx=forc_sign*fxi(i)*dt12
1598 impy=forc_sign*fyi(i)*dt12
1599 impz=forc_sign*fzi(i)*dt12
1600
1601 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1602 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1603 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1604
1605 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1606 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1607 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1608
1609 fsavsub1(11,jsub)=fsavsub1(11,jsub)+abs(fni(i))*dt12
1610
1611 IF(isensint(jsub+1)/=0) THEN
1612 fsavparit(jsub+1,1,i+nft) = forc_sign*fxi(i)
1613 fsavparit(jsub+1,2,i+nft) = forc_sign*fyi(i)
1614 fsavparit(jsub+1,3,i+nft) = forc_sign*fzi(i)
1615 ENDIF
1616
1617 jj=jj+1
1618
1619 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfs
1620
1621 iss2 = bitget(inflg_subs(jj),0)
1622 iss1 = bitget(inflg_subs(jj),1)
1623 ksub=lisubm(kk)
1624 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1625 ims2 = bitget(inflg_subm(kk),0)
1626 ims1 = bitget(inflg_subm(kk),1)
1627 IF(ksub==jsub)THEN
1628 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
1629 . (ims2 == 1 .AND. iss1 == 1))) THEN
1630 kk=kk+1
1631 ksub=lisubm(kk)
1632 cycle
1633 END IF
1634
1635 impx=forc_sign*fxi(i)*dt12
1636 impy=forc_sign*fyi(i)*dt12
1637 impz=forc_sign*fzi(i)*dt12
1638
1639
1640 IF(ims2 > 0)THEN
1641 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
1642 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
1643 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
1644 ELSE
1645 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1646 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1647 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1648 END IF
1649
1650 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1651 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1652 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1653
1654 fsavsub1(11,jsub)=fsavsub1(11,jsub)+abs(fni(i))*dt12
1655
1656
1657 IF(isensint(jsub+1)/=0) THEN
1658 IF(ims2 > 0)THEN
1659 fsavparit(jsub+1,1,i+nft) = -forc_sign*fxi(i)
1660 fsavparit(jsub+1,2,i+nft) = -forc_sign*fyi(i)
1661 fsavparit(jsub+1,3,i+nft) = -forc_sign*fzi(i)
1662 ELSE
1663 fsavparit(jsub+1,1,i+nft) = forc_sign*fxi(i)
1664 fsavparit(jsub+1,2,i+nft) = forc_sign*fyi(i)
1665 fsavparit(jsub+1,3,i+nft) = forc_sign*fzi(i)
1666 END IF
1667 ENDIF
1668 ENDIF
1669 kk=kk+1
1670 ksub=lisubm(kk)
1671 ENDDO
1672 jj=jj+1
1673 ENDIF
1674 END DO
1675
1676 END IF
1677
1678 ie=ce_loc(i)
1679
1680 kk =addsubm(ie)
1681 DO WHILE(kk<addsubm(ie+1))
1682 ksub=lisubm(kk)
1683
1684 itypsub = typsub(ksub)
1685 IF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
1686
1687 impx=-forc_sign*fxi(i)*dt12
1688 impy=-forc_sign*fyi(i)*dt12
1689 impz=-forc_sign*fzi(i)*dt12
1690
1691 fsavsub1(1,ksub)=fsavsub1(1,ksub)+impx
1692 fsavsub1(2,ksub)=fsavsub1(2,ksub)+impy
1693 fsavsub1(3,ksub)=fsavsub1(3,ksub)+impz
1694
1695 fsavsub1(8,ksub) =fsavsub1(8,ksub) +abs(impx)
1696 fsavsub1(9,ksub) =fsavsub1(9,ksub) +abs(impy)
1697 fsavsub1(10,ksub)=fsavsub1(10,ksub)+abs(impz)
1698
1699 fsavsub1(11,ksub)=fsavsub1(11,ksub)+abs(fni(i))*dt12
1700
1701 IF(isensint(ksub+1)/=0) THEN
1702 fsavparit(ksub+1,1,i+nft) = fxi(i)
1703 fsavparit(ksub+1,2,i+nft) = fyi(i)
1704 fsavparit(ksub+1,3,i+nft) = fzi(i)
1705 ENDIF
1706
1707
1708 ENDIF
1709 kk=kk+1
1710 ENDDO
1711 END DO
1712
1713 IF(nspmd>1) THEN
1714C loop split because of a PGI bug
1715 DO i=1,jlt
1716 nn = nsvg(i)
1717 IF(nn<0)THEN
1718 nn = -nn
1719 ie=ce_loc(i)
1720 jj =addsubsfi(nin)%P(nn)
1721 kk =addsubm(ie)
1722 DO WHILE(jj<addsubsfi(nin)%P(nn+1))
1723 jsub=lisubsfi(nin)%P(jj)
1724 itypsub = typsub(jsub)
1725
1726 IF(itypsub == 1 ) THEN ! Defining specific inter
1727
1728 ksub=lisubm(kk)
1729 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1730 IF(ksub==jsub)THEN
1731C if type7 sym of type19 opposite sign for normal force -> forc_sign = -1
1732 impx=forc_sign*fxi(i)*dt12
1733 impy=forc_sign*fyi(i)*dt12
1734 impz=forc_sign*fzi(i)*dt12
1735C MAIN side :
1736 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1737 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1738 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1739C
1740 IF(isensint(jsub+1)/=0) THEN
1741 fsavparit(jsub+1,1,i+nft) = forc_sign*fxi(i)
1742 fsavparit(jsub+1,2,i+nft) = forc_sign*fyi(i)
1743 fsavparit(jsub+1,3,i+nft) = forc_sign*fzi(i)
1744 ENDIF
1745C
1746 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1747 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1748 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1749C
1750 fsavsub1(11,jsub)=fsavsub1(11,jsub)+abs(fni(i))*dt12
1751C
1752 END IF
1753
1754 kk=kk+1
1755 ksub=lisubm(kk)
1756 ENDDO
1757 jj=jj+1
1758
1759 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
1760
1761 impx=forc_sign*fxi(i)*dt12
1762 impy=forc_sign*fyi(i)*dt12
1763 impz=forc_sign*fzi(i)*dt12
1764
1765 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1766 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1767 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1768
1769 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1770 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1771 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1772
1773 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1774
1775 IF(isensint(jsub+1)/=0) THEN
1776 fsavparit(jsub+1,1,i+nft) = forc_sign*fxi(i)
1777 fsavparit(jsub+1,2,i+nft) = forc_sign*fyi(i)
1778 fsavparit(jsub+1,3,i+nft) = forc_sign*fzi(i)
1779 ENDIF
1780
1781 jj=jj+1
1782
1783 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfaces
1784
1785 iss2 = bitget(inflg_subsfi(nin)%P(jj),0)
1786 iss1 = bitget(inflg_subsfi(nin)%P(jj),1)
1787 ksub=lisubm(kk)
1788 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1789 ims2 = bitget(inflg_subm(kk),0)
1790 ims1 = bitget(inflg_subm(kk),1)
1791 IF(ksub==jsub)THEN
1792 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
1793 . (ims2 == 1 .AND. iss1 == 1))) THEN
1794 kk=kk+1
1795 ksub=lisubm(kk)
1796 cycle
1797 END IF
1798
1799 impx=forc_sign*fxi(i)*dt12
1800 impy=forc_sign*fyi(i)*dt12
1801 impz=forc_sign*fzi(i)*dt12
1802
1803 IF(ims2 > 0)THEN
1804 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
1805 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
1806 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
1807 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
1808 ELSE
1809 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1810 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1811 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1812 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1813 END IF
1814C
1815 IF(isensint(jsub+1)/=0) THEN
1816 IF(ims2 > 0)THEN
1817 fsavparit(jsub+1,1,i+nft) = -forc_sign*fxi(i)
1818 fsavparit(jsub+1,2,i+nft) = -forc_sign*fyi(i)
1819 fsavparit(jsub+1,3,i+nft) = -forc_sign*fzi(i)
1820 ELSE
1821 fsavparit(jsub+1,1,i+nft) = forc_sign*fxi(i)
1822 fsavparit(jsub+1,2,i+nft) = forc_sign*fyi(i)
1823 fsavparit(jsub+1,3,i+nft) = forc_sign*fzi(i)
1824 END IF
1825 ENDIF
1826
1827 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1828 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1829 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1830
1831 ENDIF
1832 kk=kk+1
1833 ksub=lisubm(kk)
1834 ENDDO
1835 jj=jj+1
1836
1837 ENDIF
1838 END DO
1839 END IF
1840 END DO
1841 END IF
1842 END IF
1843
1844C------------For /LOAD/PRESSURE tag nodes in contact-------------
1845 IF(ninloadp > 0) THEN
1846 DO k = kloadpinter(nin)+1, kloadpinter(nin+1)
1847 pp = loadpinter(k)
1848 ppl = loadp_hyd_inter(pp)
1849 dgap = dgaploadint(k)
1850 DO i=1,jlt
1851 gapp= gapv(i) + dgap
1852 IF(pene(i) > zero .OR.(distp(i) <= gapp)) THEN
1853 tagncont(ppl,ix1(i)) = 1
1854 tagncont(ppl,ix2(i)) = 1
1855 tagncont(ppl,ix3(i)) = 1
1856 tagncont(ppl,ix4(i)) = 1
1857 jg = nsvg(i)
1858 IF(jg>0) THEN
1859C SPMD : do same after reception of forces for remote nodes
1860 tagncont(ppl,jg) = 1
1861 ENDIF
1862 ENDIF
1863
1864 ENDDO
1865 ENDDO
1866 ENDIF
1867C++++++++++++++++++++++++++++++++++++++++
1868C Friction coefficient computation
1869C++++++++++++++++++++++++++++++++++++++++
1870
1871 IF(iorthfric == 0) THEN
1872C++ Isotropic Friction
1873
1874 IF (mfrot==0) THEN
1875C--- Coulomb friction
1876 DO i=1,jlt-jlt_tied
1877 xmu(i) = fricc(i)
1878 ENDDO
1879 ELSEIF (mfrot==1) THEN
1880C--- Viscous friction
1881 DO i=1,jlt-jlt_tied
1882 ie=ce_loc(i)
1883 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1884 v2 = (vx(i) - n1(i)*aa)**2
1885 . + (vy(i) - n2(i)*aa)**2
1886 . + (vz(i) - n3(i)*aa)**2
1887 vv = sqrt(max(em30,v2))
1888 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1889 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1890 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1891 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1892 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1893 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1894 ax = ay1*az2 - az1*ay2
1895 ay = az1*ax2 - ax1*az2
1896 az = ax1*ay2 - ay1*ax2
1897 area = half*sqrt(ax*ax+ay*ay+az*az)
1898 p = -fni(i)/area
1899 xmu(i) = fricc(i) + (fric_coefs(i,1) + fric_coefs(i,4)*p ) * p
1900 . +(fric_coefs(i,2) + fric_coefs(i,3)*p) * vv + fric_coefs(i,5)*v2
1901 xmu(i) = max(xmu(i),em30)
1902 ENDDO
1903 ELSEIF(mfrot==2)THEN
1904C--- Darmstad Law
1905 DO i=1,jlt-jlt_tied
1906 ie=ce_loc(i)
1907 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1908 v2 = (vx(i) - n1(i)*aa)**2
1909 . + (vy(i) - n2(i)*aa)**2
1910 . + (vz(i) - n3(i)*aa)**2
1911 vv = sqrt(max(em30,v2))
1912 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1913 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1914 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1915 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1916 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1917 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1918 ax = ay1*az2 - az1*ay2
1919 ay = az1*ax2 - ax1*az2
1920 az = ax1*ay2 - ay1*ax2
1921 area = half*sqrt(ax*ax+ay*ay+az*az)
1922 p = -fni(i)/area
1923 xmu(i) = fricc(i)
1924 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
1925 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
1926 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
1927 xmu(i) = max(xmu(i),em30)
1928 ENDDO
1929 ELSEIF (mfrot==3) THEN
1930C--- Renard Law
1931 DO i=1,jlt-jlt_tied
1932 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1933 v2 = (vx(i) - n1(i)*aa)**2
1934 . + (vy(i) - n2(i)*aa)**2
1935 . + (vz(i) - n3(i)*aa)**2
1936 vv = sqrt(max(em30,v2))
1937 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
1938 dmu = fric_coefs(i,3)-fric_coefs(i,1)
1939 vv1 = vv / fric_coefs(i,5)
1940 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
1941 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
1942 dmu = fric_coefs(i,4)-fric_coefs(i,3)
1943 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
1944 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
1945 ELSE
1946 dmu = fric_coefs(i,2)-fric_coefs(i,4)
1947 vv2 = (vv - fric_coefs(i,6))**2
1948 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
1949 ENDIF
1950 xmu(i) = max(xmu(i),em30)
1951 ENDDO
1952
1953 ELSEIF(mfrot==4)THEN
1954C--- Exponential decay model
1955 DO i=1,jlt-jlt_tied
1956 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1957 v2 = (vx(i) - n1(i)*aa)**2
1958 . + (vy(i) - n2(i)*aa)**2
1959 . + (vz(i) - n3(i)*aa)**2
1960 vv = sqrt(max(em30,v2))
1961 xmu(i) = fric_coefs(i,1)
1962 . + (fricc(i)-fric_coefs(i,1))*exp(-fric_coefs(i,2)*vv)
1963 xmu(i) = max(xmu(i),em30)
1964 ENDDO
1965 ENDIF
1966 ELSE
1967C++ Orthotropic Friction
1968
1969 IF (mfrot==0) THEN
1970C--- Coulomb friction
1971#include "vectorize.inc"
1972 DO k=1,nfisot
1973 i = indexisot(k)
1974 xmu(i) = fricc(i)
1975 ENDDO
1976#include "vectorize.inc"
1977 DO k=1,nforth
1978 i = indexorth(k)
1979 xmu(i) = fricc(i)
1980 xmu2(i) = fricc2(i)
1981 IF(xmu(i)<=em30) THEN
1982 xmu(i) = em30
1983 dir1(i,1:3) = zero
1984 an(k) = zero
1985 ELSE
1986 an(k)=xmu(i)**2
1987 an(k)=one/an(k)
1988 ENDIF
1989 IF(xmu2(i)<=em30) THEN
1990 xmu2(i) = em30
1991 dir2(i,1:3) = zero
1992 bn(k) = zero
1993 ELSE
1994 bn(k)=xmu2(i)**2
1995 bn(k)=one/bn(k)
1996 ENDIF
1997 ENDDO
1998
1999
2000 ELSEIF (mfrot==1) THEN
2001C--- Viscous friction
2002#include "vectorize.inc"
2003 DO k=1,nfisot ! isotropic friction couples
2004 i = indexisot(k)
2005 ie=ce_loc(i)
2006 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
2007 v2 = (vx(i) - n1(i)*aa)**2
2008 . + (vy(i) - n2(i)*aa)**2
2009 . + (vz(i) - n3(i)*aa)**2
2010 vv = sqrt(max(em30,v2))
2011
2012 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
2013 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
2014 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
2015 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
2016 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
2017 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
2018 ax = ay1*az2 - az1*ay2
2019 ay = az1*ax2 - ax1*az2
2020 az = ax1*ay2 - ay1*ax2
2021 area = half*sqrt(ax*ax+ay*ay+az*az)
2022 p = -fni(i)/area
2023 xmu(i) = fricc(i) + (fric_coefs(i,1) + fric_coefs(i,4)*p ) * p
2024 . +(fric_coefs(i,2) + fric_coefs(i,3)*p) * vv + fric_coefs(i,5)*v2
2025 xmu(i) = max(xmu(i),em30)
2026 ENDDO
2027
2028#include "vectorize.inc"
2029 DO k=1,nforth ! Orthotropic friction couples
2030 i = indexorth(k)
2031 ie=ce_loc(i)
2032c
2033 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
2034 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
2035 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
2036 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
2037 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
2038 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
2039 ax = ay1*az2 - az1*ay2
2040 ay = az1*ax2 - ax1*az2
2041 az = ax1*ay2 - ay1*ax2
2042 area = half*sqrt(ax*ax+ay*ay+az*az)
2043 p = -fni(i)/area
2044c
2045 v2 = vx(i)*dir1(i,1) +vy(i)*dir1(i,2)+vz(i)*dir1(i,3)
2046 vv = max(em30,v2)
2047 xmu(i) = fricc(i) + (fric_coefs(i,1) + fric_coefs(i,4)*p ) * p
2048 . +(fric_coefs(i,2) + fric_coefs(i,3)*p) * vv + fric_coefs(i,5)*vv*vv
2049
2050 v2 = vx(i)*dir2(i,1) +vy(i)*dir2(i,2)+vz(i)*dir2(i,3)
2051 vv = max(em30,v2)
2052 xmu2(i) = fricc2(i) + (fric_coefs2(i,1) + fric_coefs2(i,4)*p ) * p
2053 . +(fric_coefs2(i,2) + fric_coefs2(i,3)*p) * vv + fric_coefs2(i,5)*vv*vv
2054
2055 xmu(i) = max(xmu(i),em30)
2056 xmu2(i) = max(xmu2(i),em30)
2057 ENDDO
2058
2059#include "vectorize.inc"
2060 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
2061 i = indexorth(k)
2062 IF(xmu(i)<=em30) THEN
2063 xmu(i) = em30
2064 dir1(i,1:3) = zero
2065 an(k) = zero
2066 ELSE
2067 an(k)=xmu(i)**2
2068 an(k)=one/an(k)
2069 ENDIF
2070 IF(xmu2(i)<=em30) THEN
2071 xmu2(i) = em30
2072 dir2(i,1:3) = zero
2073 bn(k) = zero
2074 ELSE
2075 bn(k)=xmu2(i)**2
2076 bn(k)=one/bn(k)
2077 ENDIF
2078 ENDDO
2079
2080 ELSEIF(mfrot==2)THEN
2081C--- Darmstad LAW
2082#include "vectorize.inc"
2083 DO k=1,nfisot ! isotropic friction couples
2084 i = indexisot(k)
2085 ie=ce_loc(i)
2086 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
2087 v2 = (vx(i) - n1(i)*aa)**2
2088 . + (vy(i) - n2(i)*aa)**2
2089 . + (vz(i) - n3(i)*aa)**2
2090 vv = sqrt(max(em30,v2))
2091 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
2092 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
2093 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
2094 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
2095 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
2096 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
2097 ax = ay1*az2 - az1*ay2
2098 ay = az1*ax2 - ax1*az2
2099 az = ax1*ay2 - ay1*ax2
2100 area = half*sqrt(ax*ax+ay*ay+az*az)
2101 p = -fni(i)/area
2102 xmu(i) = fricc(i)
2103 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
2104 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
2105 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
2106 xmu(i) = max(xmu(i),em30)
2107 ENDDO
2108c
2109#include "vectorize.inc"
2110 DO k=1,nforth ! Orthotropic friction couples
2111 i = indexorth(k)
2112 ie=ce_loc(i)
2113c
2114 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
2115 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
2116 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
2117 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
2118 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
2119 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
2120 ax = ay1*az2 - az1*ay2
2121 ay = az1*ax2 - ax1*az2
2122 az = ax1*ay2 - ay1*ax2
2123 area = half*sqrt(ax*ax+ay*ay+az*az)
2124 p = -fni(i)/area
2125c
2126 v2 = vx(i)*dir1(i,1) +vy(i)*dir1(i,2)+vz(i)*dir1(i,3)
2127 vv = max(em30,v2)
2128 xmu(i) = fricc(i)
2129 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
2130 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
2131 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
2132c
2133 v2 = vx(i)*dir2(i,1) +vy(i)*dir2(i,2)+vz(i)*dir2(i,3)
2134 vv = max(em30,v2)
2135 xmu2(i) = fricc2(i)
2136 . + fric_coefs2(i,1)*exp(fric_coefs2(i,2)*vv)*p*p
2137 . + fric_coefs2(i,3)*exp(fric_coefs2(i,4)*vv)*p
2138 . + fric_coefs2(i,5)*exp(fric_coefs2(i,6)*vv)
2139 xmu(i) = max(xmu(i),em30)
2140 xmu2(i) = max(xmu2(i),em30)
2141 ENDDO
2142
2143#include "vectorize.inc"
2144 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
2145 i = indexorth(k)
2146 IF(xmu(i)<=em30) THEN
2147 xmu(i) = em30
2148 dir1(i,1:3) = zero
2149 an(k) = zero
2150 ELSE
2151 an(k)=xmu(i)**2
2152 an(k)=one/an(k)
2153 ENDIF
2154 IF(xmu2(i)<=em30) THEN
2155 xmu2(i) = em30
2156 dir2(i,1:3) = zero
2157 bn(k) = zero
2158 ELSE
2159 bn(k)=xmu2(i)**2
2160 bn(k)=one/bn(k)
2161 ENDIF
2162 ENDDO
2163
2164 ELSEIF (mfrot==3) THEN
2165C--- Renard Law
2166#include "vectorize.inc"
2167 DO k=1,nfisot ! isotropic friction couples
2168 i = indexisot(k)
2169 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
2170 v2 = (vx(i) - n1(i)*aa)**2
2171 . + (vy(i) - n2(i)*aa)**2
2172 . + (vz(i) - n3(i)*aa)**2
2173 vv = sqrt(max(em30,v2))
2174 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
2175 dmu = fric_coefs(i,3)-fric_coefs(i,1)
2176 vv1 = vv / fric_coefs(i,5)
2177 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
2178 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
2179 dmu = fric_coefs(i,4)-fric_coefs(i,3)
2180 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
2181 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
2182 ELSE
2183 dmu = fric_coefs(i,2)-fric_coefs(i,4)
2184 vv2 = (vv - fric_coefs(i,6))**2
2185 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
2186 ENDIF
2187 xmu(i) = max(xmu(i),em30)
2188 ENDDO
2189
2190#include "vectorize.inc"
2191 DO k=1,nforth ! Orthotropic friction couples
2192 i = indexorth(k)
2193c
2194 v2 = vx(i)*dir1(i,1) +vy(i)*dir1(i,2)+vz(i)*dir1(i,3)
2195 vv = max(em30,v2)
2196 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
2197 dmu = fric_coefs(i,3)-fric_coefs(i,1)
2198 vv1 = vv / fric_coefs(i,5)
2199 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
2200 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
2201 dmu = fric_coefs(i,4)-fric_coefs(i,3)
2202 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
2203 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
2204 ELSE
2205 dmu = fric_coefs(i,2)-fric_coefs(i,4)
2206 vv2 = (vv - fric_coefs(i,6))**2
2207 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
2208 ENDIF
2209
2210 v2 = vx(i)*dir2(i,1) +vy(i)*dir2(i,2)+vz(i)*dir2(i,3)
2211 vv = max(em30,v2)
2212 IF(vv>=0.AND.vv<=fric_coefs2(i,5)) THEN
2213 dmu = fric_coefs2(i,3)-fric_coefs2(i,1)
2214 vv1 = vv / fric_coefs2(i,5)
2215 xmu2(i) = fric_coefs2(i,1)+ dmu*vv1*(two-vv1)
2216 ELSEIF(vv>fric_coefs2(i,5).AND.vv<fric_coefs2(i,6)) THEN
2217 dmu = fric_coefs2(i,4)-fric_coefs2(i,3)
2218 vv1 = (vv - fric_coefs2(i,5))/(fric_coefs2(i,6)-fric_coefs2(i,5))
2219 xmu2(i) = fric_coefs2(i,3)+ dmu * (three-two*vv1)*vv1**2
2220 ELSE
2221 dmu = fric_coefs2(i,2)-fric_coefs2(i,4)
2222 vv2 = (vv - fric_coefs2(i,6))**2
2223 xmu2(i) = fric_coefs2(i,2) - dmu / (one + dmu*vv2)
2224 ENDIF
2225 xmu(i) = max(xmu(i),em30)
2226 xmu2(i) = max(xmu2(i),em30)
2227 ENDDO
2228
2229#include "vectorize.inc"
2230 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
2231 i = indexorth(k)
2232 IF(xmu(i)<=em30) THEN
2233 xmu(i) = em30
2234 dir1(i,1:3) = zero
2235 an(k) = zero
2236 ELSE
2237 an(k)=xmu(i)**2
2238 an(k)=one/an(k)
2239 ENDIF
2240 IF(xmu2(i)<=em30) THEN
2241 xmu2(i) = em30
2242 dir2(i,1:3) = zero
2243 bn(k) = zero
2244 ELSE
2245 bn(k)=xmu2(i)**2
2246 bn(k)=one/bn(k)
2247 ENDIF
2248 ENDDO
2249
2250 ELSEIF(mfrot==4)THEN
2251C--- Exponential decay model
2252#include "vectorize.inc"
2253 DO k=1,nfisot ! isotropic friction couples
2254 i = indexisot(k)
2255 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
2256 v2 = (vx(i) - n1(i)*aa)**2
2257 . + (vy(i) - n2(i)*aa)**2
2258 . + (vz(i) - n3(i)*aa)**2
2259 vv = sqrt(max(em30,v2))
2260 xmu(i) = fric_coefs(i,1)
2261 . + (fricc(i)-fric_coefs(i,1))*exp(-fric_coefs(i,2)*vv)
2262 xmu(i) = max(xmu(i),em30)
2263 ENDDO
2264c
2265#include "vectorize.inc"
2266 DO k=1,nforth ! Orthotropic friction couples
2267 i = indexorth(k)
2268c
2269 v2 = vx(i)*dir1(i,1) +vy(i)*dir1(i,2)+vz(i)*dir1(i,3)
2270 vv = max(em30,v2)
2271 xmu(i) = fricc(i)
2272 . + (fric_coefs(i,1)-fricc(i))*exp(-fric_coefs(i,2)*vv)
2273c
2274 v2 = vx(i)*dir2(i,1) +vy(i)*dir2(i,2)+vz(i)*dir2(i,3)
2275 vv = max(em30,v2)
2276 xmu2(i) = fric_coefs2(i,1)
2277 . + (fricc2(i)-fric_coefs2(i,1))*exp(-fric_coefs2(i,2)*vv)
2278
2279 xmu(i) = max(xmu(i),em30)
2280 xmu2(i) = max(xmu2(i),em30)
2281 ENDDO
2282
2283#include "vectorize.inc"
2284 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
2285 i = indexorth(k)
2286 IF(xmu(i)<=em30) THEN
2287 xmu(i) = em30
2288 dir1(i,1:3) = zero
2289 an(k) = zero
2290 ELSE
2291 an(k)=xmu(i)**2
2292 an(k)=one/an(k)
2293 ENDIF
2294 IF(xmu2(i)<=em30) THEN
2295 xmu2(i) = em30
2296 dir2(i,1:3) = zero
2297 bn(k) = zero
2298 ELSE
2299 bn(k)=xmu2(i)**2
2300 bn(k)=one/bn(k)
2301 ENDIF
2302 ENDDO
2303
2304 ENDIF
2305
2306 ENDIF ! IORTHFRIC
2307
2308C------------------
2309C TANGENT FORCE CALCULATION
2310C------------------
2311 fsav4 = zero
2312 fsav5 = zero
2313 fsav6 = zero
2314 fsav12= zero
2315 fsav13= zero
2316 fsav14= zero
2317 fsav15= zero
2318
2319 IF(iorthfric ==0 ) THEN
2320C++ Isotropic Friction
2321C---------------------------------
2322 IF (ifq>=10) THEN
2323C---------------------------------
2324C INCREMENTAL (STIFFNESS) FORMULATION
2325C---------------------------------
2326 IF (ifq==13) THEN
2327 alpha = max(one,alpha0*dt12)
2328 ELSE
2329 alpha = alpha0
2330 ENDIF
2331 IF (inconv==1) THEN
2332 DO i=1,jlt-jlt_tied
2333 fx = stif0(i)*vx(i)*dt12
2334 fy = stif0(i)*vy(i)*dt12
2335 fz = stif0(i)*vz(i)*dt12
2336 fx = cand_fx(index(i)) + alpha*fx
2337 fy = cand_fy(index(i)) + alpha*fy
2338 fz = cand_fz(index(i)) + alpha*fz
2339 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2340 fx = fx - ftn*n1(i)
2341 fy = fy - ftn*n2(i)
2342 fz = fz - ftn*n3(i)
2343 ft = fx*fx + fy*fy + fz*fz
2344 ft = max(ft,em30)
2345 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2346 beta = min(one,xmu(i)*sqrt(fn/ft))
2347 fxt(i) = fx * beta
2348 fyt(i) = fy * beta
2349 fzt(i) = fz * beta
2350 cand_fx(index(i)) = fxt(i)
2351 cand_fy(index(i)) = fyt(i)
2352 cand_fz(index(i)) = fzt(i)
2353 ifpen(index(i)) = 1
2354C------- total force
2355 fxi(i)=fxi(i) + fxt(i)
2356 fyi(i)=fyi(i) + fyt(i)
2357 fzi(i)=fzi(i) + fzt(i)
2358C---------------------------------
2359C CONTACT ENERGY CALCULATION
2360C---------------------------------
2361 IF( intth > 0 .AND.beta/=zero) THEN
2362 efrict(i) = (fx-fxt(i))*fxt(i) + (fy-fyt(i))*fyt(i) +
2363 . (fz-fzt(i))*fzt(i)
2364 efrict(i) = efrict(i)/stif0(i) ! FRICTIONAL ENERGY
2365 qfrict = qfrict + efrict(i)
2366 ENDIF
2367 efric_l(i) = dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
2368 econvt = econvt + efric_l(i)
2369
2370 ENDDO
2371
2372
2373C--------implicit non converge---
2374 ELSE
2375 DO i=1,jlt-jlt_tied
2376 fx = stif0(i)*vx(i)*dt12
2377 fy = stif0(i)*vy(i)*dt12
2378 fz = stif0(i)*vz(i)*dt12
2379 fx = cand_fx(index(i)) + alpha*fx
2380 fy = cand_fy(index(i)) + alpha*fy
2381 fz = cand_fz(index(i)) + alpha*fz
2382 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2383 fx = fx - ftn*n1(i)
2384 fy = fy - ftn*n2(i)
2385 fz = fz - ftn*n3(i)
2386 ft = fx*fx + fy*fy + fz*fz
2387 ft = max(ft,em30)
2388 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2389 beta = min(one,xmu(i)*sqrt(fn/ft))
2390 fxt(i) = fx * beta
2391 fyt(i) = fy * beta
2392 fzt(i) = fz * beta
2393 fxi(i)=fxi(i) + fxt(i)
2394 fyi(i)=fyi(i) + fyt(i)
2395 fzi(i)=fzi(i) + fzt(i)
2396 ifpen(index(i)) = 1
2397 ENDDO
2398 ENDIF
2399C---------------------------------
2400C TOTAL (VISCOUS) FORMULATION + FRICTION FILTERING
2401C---------------------------------
2402 ELSEIF (ifq>0) THEN
2403 IF (ifq==3) THEN
2404 alpha = max(one,alpha0*dt12)
2405 ELSE
2406 alpha = alpha0
2407 ENDIF
2408 alphi = one - alpha
2409 DO i=1,jlt-jlt_tied
2410 vnx = n1(i)*vn(i)
2411 vny = n2(i)*vn(i)
2412 vnz = n3(i)*vn(i)
2413 vx(i) = vx(i) - vnx
2414 vy(i) = vy(i) - vny
2415 vz(i) = vz(i) - vnz
2416 v2 = vx(i)**2 + vy(i)**2 + vz(i)**2
2417 vis2(i) = viscffric(i) * vis2(i)
2418 fm2 = (xmu(i)*fni(i))**2
2419 f2 = vis2(i) * v2
2420 a2 = min(f2,fm2) / max(em30,f2)
2421 aa = sqrt(a2 * vis2(i))
2422 fx = aa * vx(i)
2423 fy = aa * vy(i)
2424 fz = aa * vz(i)
2425 fxt(i) = alpha*fx + alphi*cand_fx(index(i))
2426 fyt(i) = alpha*fy + alphi*cand_fy(index(i))
2427 fzt(i) = alpha*fz + alphi*cand_fz(index(i))
2428 cand_fx(index(i)) = fxt(i)
2429 cand_fy(index(i)) = fyt(i)
2430 cand_fz(index(i)) = fzt(i)
2431 ifpen(index(i)) = 1
2432C------- total force
2433 fxi(i) = fxi(i) + fxt(i)
2434 fyi(i) = fyi(i) + fyt(i)
2435 fzi(i) = fzi(i) + fzt(i)
2436C---------------------------------
2437C CONTACT ENERGY CALCULATION
2438C---------------------------------
2439 efrict(i) = zero
2440 IF( intth > 0) THEN
2441 efrict(i) = dt1*(fxt(i)*vx(i) + fyt(i)*vy(i) + fzt(i)*vz(i)) ! FRICTIONAL ENERGY
2442 qfrict = qfrict + efrict(i)
2443 ENDIF
2444 efric_l(i) = dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
2445 econvt = econvt + efric_l(i)
2446 ENDDO
2447 ELSE
2448C---------------------------------
2449C TOTAL (VISCOUS) FORMULATION / NO FRICTION FILTERING
2450C---------------------------------
2451 DO i=1,jlt-jlt_tied
2452 vnx = n1(i)*vn(i)
2453 vny = n2(i)*vn(i)
2454 vnz = n3(i)*vn(i)
2455 vx(i) = vx(i) - vnx
2456 vy(i) = vy(i) - vny
2457 vz(i) = vz(i) - vnz
2458 v2 = vx(i)**2 + vy(i)**2 + vz(i)**2
2459 vis2(i) = viscffric(i) * vis2(i)
2460 fm2 = (xmu(i)*fni(i))**2
2461 f2 = vis2(i) * v2
2462 a2 = min(f2,fm2) / max(em30,f2)
2463 aa = sqrt(a2 * vis2(i))
2464 fxt(i) = aa * vx(i)
2465 fyt(i) = aa * vy(i)
2466 fzt(i) = aa * vz(i)
2467C------- total force
2468 fxi(i)=fxi(i) + fxt(i)
2469 fyi(i)=fyi(i) + fyt(i)
2470 fzi(i)=fzi(i) + fzt(i)
2471C---------------------------------
2472C CONTACT ENERGY CALCULATION
2473C---------------------------------
2474 efrict(i) = zero
2475 IF( intth > 0) THEN
2476 efrict(i) = aa * v2 * dt1 ! FRICTIONAL ENERGY
2477 qfrict = qfrict + efrict(i)
2478 ENDIF
2479 efric_l(i) = aa * v2 * dt1
2480 econvt = econvt + efric_l(i)
2481 ENDDO
2482 ENDIF
2483
2484
2485 ELSE
2486C++ Orthotropic Friction
2487
2488C---------------------------------
2489 IF (ifq>=10) THEN
2490C---------------------------------
2491C INCREMENTAL (STIFFNESS) FORMULATION
2492C---------------------------------
2493 IF (ifq==13) THEN
2494 alpha = max(one,alpha0*dt12)
2495 ELSE
2496 alpha = alpha0
2497 ENDIF
2498 IF (inconv==1) THEN
2499#include "vectorize.inc"
2500 DO k=1,nfisot ! isotropic friction couples
2501 i = indexisot(k)
2502 fx = stif0(i)*vx(i)*dt12
2503 fy = stif0(i)*vy(i)*dt12
2504 fz = stif0(i)*vz(i)*dt12
2505 fx = cand_fx(index(i)) + alpha*fx
2506 fy = cand_fy(index(i)) + alpha*fy
2507 fz = cand_fz(index(i)) + alpha*fz
2508 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2509 fx = fx - ftn*n1(i)
2510 fy = fy - ftn*n2(i)
2511 fz = fz - ftn*n3(i)
2512 ft = fx*fx + fy*fy + fz*fz
2513 ft = max(ft,em30)
2514 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2515
2516 beta = min(one,xmu(i)*sqrt(fn/ft))
2517
2518 fxt(i) = fx * beta
2519 fyt(i) = fy * beta
2520 fzt(i) = fz * beta
2521
2522 cand_fx(index(i)) = fxt(i)
2523 cand_fy(index(i)) = fyt(i)
2524 cand_fz(index(i)) = fzt(i)
2525 ifpen(index(i)) = 1
2526C------- total force
2527 fxi(i)=fxi(i) + fxt(i)
2528 fyi(i)=fyi(i) + fyt(i)
2529 fzi(i)=fzi(i) + fzt(i)
2530
2531C---------------------------------
2532C CONTACT ENERGY CALCULATION
2533C---------------------------------
2534 efrict(i) = zero
2535 IF( intth > 0 .AND.beta/=zero) THEN
2536 efrict(i) = (fx-fxt(i))*fxt(i) + (fy-fyt(i))*fyt(i) +
2537 . (fz-fzt(i))*fzt(i)
2538 efrict(i) = efrict(i)/stif0(i) ! FRICTIONAL ENERGY
2539 qfrict = qfrict + efrict(i)
2540 ENDIF
2541 efric_l(i)= dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
2542 econvt = econvt + efric_l(i)
2543 ENDDO
2544
2545#include "vectorize.inc"
2546 DO k=1,nforth ! Orthotropic friction couples
2547 i = indexorth(k)
2548 fx = stif0(i)*vx(i)*dt12
2549 fy = stif0(i)*vy(i)*dt12
2550 fz = stif0(i)*vz(i)*dt12
2551 fx = cand_fx(index(i)) + alpha*fx
2552 fy = cand_fy(index(i)) + alpha*fy
2553 fz = cand_fz(index(i)) + alpha*fz
2554
2555 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2556
2557 fx = fx - ftn*n1(i)
2558 fy = fy - ftn*n2(i)
2559 fz = fz - ftn*n3(i)
2560
2561 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2562 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2563
2564 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2565 ft = max(ft,em30)
2566 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2567
2568 beta = fn/ft
2569
2570 IF(beta == 0 ) THEN
2571 fxt(i) = zero
2572 fyt(i) = zero
2573 fzt(i) = zero
2574 ELSEIF(beta > 1) THEN ! Inside the ellipse
2575 fxt(i) = fx
2576 fyt(i) = fy
2577 fzt(i) = fz
2578 ELSE ! outside the ellipse
2579
2580! Projection on local tangent of ellipse (outside of ellipse)
2581! ANS = (Fadh-Fproj).n
2582 nep1 =ftt1*an(k)/fn
2583 nep2 =ftt2*bn(k)/fn
2584 nep =nep1*nep1+nep2*nep2
2585 nep =sqrt(nep)
2586
2587 ep=nep1*ftt1+nep2*ftt2
2588
2589 ans=(ep-sqrt(ep))/max(em20,nep)
2590 nep1 =nep1/max(em20,nep)
2591 nep2 =nep2/max(em20,nep)
2592
2593! Projection on ellipse
2594 c11 =ftt1-ans*nep1
2595 c22 =ftt2-ans*nep2
2596
2597 alphaf = atan(c22/c11)
2598
2599 signc = ftt1/max(em20,abs(ftt1))
2600 csa = signc*abs(cos(alphaf))
2601 signc = ftt2/max(em20,abs(ftt2))
2602 sna = signc*abs(sin(alphaf))
2603! Ft computation
2604 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2605 ftt1 = ft * csa
2606 ftt2 = ft * sna
2607
2608 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2609 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2610 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2611
2612 ENDIF
2613 cand_fx(index(i)) = fxt(i)
2614 cand_fy(index(i)) = fyt(i)
2615 cand_fz(index(i)) = fzt(i)
2616 ifpen(index(i)) = 1
2617C------- total force
2618 fxi(i)=fxi(i) + fxt(i)
2619 fyi(i)=fyi(i) + fyt(i)
2620 fzi(i)=fzi(i) + fzt(i)
2621
2622C---------------------------------
2623C CONTACT ENERGY CALCULATION
2624C---------------------------------
2625 efrict(i) = zero
2626 IF( intth > 0 .AND.beta/=zero) THEN
2627 efrict(i) = (fx-fxt(i))*fxt(i) + (fy-fyt(i))*fyt(i) +
2628 . (fz-fzt(i))*fzt(i)
2629 efrict(i) = efrict(i)/stif0(i) ! FRICTIONAL ENERGY
2630 qfrict = qfrict + efrict(i)
2631 ENDIF
2632 efric_l(i)= dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
2633 econvt = econvt + efric_l(i)
2634 ENDDO
2635
2636C--------implicit non converge---
2637 ELSE
2638#include "vectorize.inc"
2639 DO k=1,nfisot ! isotropic friction couples
2640 i = indexisot(k)
2641 fx = stif0(i)*vx(i)*dt12
2642 fy = stif0(i)*vy(i)*dt12
2643 fz = stif0(i)*vz(i)*dt12
2644 fx = cand_fx(index(i)) + alpha*fx
2645 fy = cand_fy(index(i)) + alpha*fy
2646 fz = cand_fz(index(i)) + alpha*fz
2647 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2648 fx = fx - ftn*n1(i)
2649 fy = fy - ftn*n2(i)
2650 fz = fz - ftn*n3(i)
2651 ft = fx*fx + fy*fy + fz*fz
2652 ft = max(ft,em30)
2653 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2654 beta = min(one,xmu(i)*sqrt(fn/ft))
2655 fxt(i) = fx * beta
2656 fyt(i) = fy * beta
2657 fzt(i) = fz * beta
2658 fxi(i)=fxi(i) + fxt(i)
2659 fyi(i)=fyi(i) + fyt(i)
2660 fzi(i)=fzi(i) + fzt(i)
2661 ifpen(index(i)) = 1
2662 ENDDO
2663
2664#include "vectorize.inc"
2665 DO k=1,nforth ! Orthotropic friction couples
2666 i = indexorth(k)
2667 fx = stif0(i)*vx(i)*dt12
2668 fy = stif0(i)*vy(i)*dt12
2669 fz = stif0(i)*vz(i)*dt12
2670 fx = cand_fx(index(i)) + alpha*fx
2671 fy = cand_fy(index(i)) + alpha*fy
2672 fz = cand_fz(index(i)) + alpha*fz
2673
2674 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2675
2676 fx = fx - ftn*n1(i)
2677 fy = fy - ftn*n2(i)
2678 fz = fz - ftn*n3(i)
2679
2680 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2681 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2682 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2683 ft = max(ft,em30)
2684 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2685
2686 beta = fn/ft
2687
2688 IF(beta == 0 ) THEN
2689 fxt(i) = zero
2690 fyt(i) = zero
2691 fzt(i) = zero
2692 ELSEIF(beta > 1) THEN ! Inside the ellipse
2693 fxt(i) = fx
2694 fyt(i) = fy
2695 fzt(i) = fz
2696 ELSE ! outside the ellipse
2697
2698! Projection on local tangent of ellipse (outside of ellipse)
2699! ANS = (Fadh-Fproj).n
2700 nep1 =ftt1*an(k)/fn
2701 nep2 =ftt2*bn(k)/fn
2702 nep =nep1*nep1+nep2*nep2
2703 nep =sqrt(nep)
2704
2705 ep=nep1*ftt1+nep2*ftt2
2706
2707 ans=(ep-sqrt(ep))/max(em20,nep)
2708 nep1 =nep1/max(em20,nep)
2709 nep2 =nep2/max(em20,nep)
2710
2711! Projection on ellipse
2712 c11 =ftt1-ans*nep1
2713 c22 =ftt2-ans*nep2
2714
2715 alphaf = atan(c22/c11)
2716
2717 signc = ftt1/max(em20,abs(ftt1))
2718 csa = signc*abs(cos(alphaf))
2719 signc = ftt2/max(em20,abs(ftt2))
2720 sna = signc*abs(sin(alphaf))
2721! Ft computation
2722 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2723 ftt1 = ft * csa
2724 ftt2 = ft * sna
2725
2726 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2727 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2728 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2729
2730 ENDIF
2731
2732 fxi(i)=fxi(i) + fxt(i)
2733 fyi(i)=fyi(i) + fyt(i)
2734 fzi(i)=fzi(i) + fzt(i)
2735 ifpen(index(i)) = 1
2736 ENDDO
2737
2738
2739 ENDIF
2740C---------------------------------
2741C TOTAL (VISCOUS) FORMULATION + FRICTION FILTERING
2742C---------------------------------
2743 ELSEIF (ifq>0) THEN
2744 IF (ifq==3) THEN
2745 alpha = max(one,alpha0*dt12)
2746 ELSE
2747 alpha = alpha0
2748 ENDIF
2749 alphi = one - alpha
2750#include "vectorize.inc"
2751 DO k=1,nfisot ! isotropic friction couples
2752 i = indexisot(k)
2753 vnx = n1(i)*vn(i)
2754 vny = n2(i)*vn(i)
2755 vnz = n3(i)*vn(i)
2756 vx(i) = vx(i) - vnx
2757 vy(i) = vy(i) - vny
2758 vz(i) = vz(i) - vnz
2759 v2 = vx(i)**2 + vy(i)**2 + vz(i)**2
2760 vis2(i) = viscffric(i) * vis2(i)
2761 fm2 = (xmu(i)*fni(i))**2
2762 f2 = vis2(i) * v2
2763 a2 = min(f2,fm2) / max(em30,f2)
2764 aa = sqrt(a2 * vis2(i))
2765 fx = aa * vx(i)
2766 fy = aa * vy(i)
2767 fz = aa * vz(i)
2768 fxt(i) = alpha*fx + alphi*cand_fx(index(i))
2769 fyt(i) = alpha*fy + alphi*cand_fy(index(i))
2770 fzt(i) = alpha*fz + alphi*cand_fz(index(i))
2771 cand_fx(index(i)) = fxt(i)
2772 cand_fy(index(i)) = fyt(i)
2773 cand_fz(index(i)) = fzt(i)
2774 ifpen(index(i)) = 1
2775C------- total force
2776 fxi(i) = fxi(i) + fxt(i)
2777 fyi(i) = fyi(i) + fyt(i)
2778 fzi(i) = fzi(i) + fzt(i)
2779C---------------------------------
2780C CONTACT ENERGY CALCULATION
2781C---------------------------------
2782 efrict(i) = zero
2783 IF( intth > 0) THEN
2784 efrict(i) = dt1*(fxt(i)*vx(i) + fyt(i)*vy(i) + fzt(i)*vz(i)) ! FRICTIONAL ENERGY
2785 qfrict = qfrict + efrict(i)
2786 ENDIF
2787 efric_l(i)= dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
2788 econvt = econvt + efric_l(i)
2789 ENDDO
2790
2791#include "vectorize.inc"
2792 DO k=1,nforth ! orthotropic friction couples
2793 i = indexorth(k)
2794 vnx = n1(i)*vn(i)
2795 vny = n2(i)*vn(i)
2796 vnz = n3(i)*vn(i)
2797 vx(i) = vx(i) - vnx
2798 vy(i) = vy(i) - vny
2799 vz(i) = vz(i) - vnz
2800
2801 vis2(i) = viscffric(i) * vis2(i)
2802
2803 vis = sqrt(vis2(i))
2804 fx = vis*vx(i)
2805 fy = vis*vy(i)
2806 fz = vis*vz(i)
2807
2808 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2809 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2810 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2811 ft = max(ft,em30)
2812 fn = fni(i)*fni(i)
2813
2814 beta = fn/ft
2815
2816 IF(beta == 0 ) THEN
2817 fxt(i) = zero
2818 fyt(i) = zero
2819 fzt(i) = zero
2820 ELSEIF(beta > 1) THEN ! Inside the ellipse
2821 fxt(i) = fx
2822 fyt(i) = fy
2823 fzt(i) = fz
2824 ELSE ! outside the ellipse
2825
2826! projection on local tangent of ellipse(outside of ellipse)
2827! ANS = (Fadh-Fproj).n
2828 nep1 =ftt1*an(k)/fn
2829 nep2 =ftt2*bn(k)/fn
2830 nep =nep1*nep1+nep2*nep2
2831 nep =sqrt(nep)
2832
2833 ep=nep1*ftt1+nep2*ftt2
2834
2835 ans=(ep-sqrt(ep))/max(em20,nep)
2836 nep1 =nep1/max(em20,nep)
2837 nep2 =nep2/max(em20,nep)
2838
2839! Projection on ellipse
2840 c11 =ftt1-ans*nep1
2841 c22 =ftt2-ans*nep2
2842
2843 alphaf = atan(c22/c11)
2844
2845 signc = ftt1/max(em20,abs(ftt1))
2846 csa = signc*abs(cos(alphaf))
2847 signc = ftt2/max(em20,abs(ftt2))
2848 sna = signc*abs(sin(alphaf))
2849! Ft computation
2850 ft = fni(i) /sqrt( (csa*csa*an(k) + sna*sna*bn(k)))
2851 ftt1 = ft * csa
2852 ftt2 = ft * sna
2853
2854 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2855 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2856 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2857
2858 ENDIF
2859
2860 fxt(i) = alpha*fxt(i) + alphi*cand_fx(index(i))
2861 fyt(i) = alpha*fyt(i) + alphi*cand_fy(index(i))
2862 fzt(i) = alpha*fzt(i) + alphi*cand_fz(index(i))
2863 cand_fx(index(i)) = fxt(i)
2864 cand_fy(index(i)) = fyt(i)
2865 cand_fz(index(i)) = fzt(i)
2866 ifpen(index(i)) = 1
2867C------- total force
2868 fxi(i) = fxi(i) + fxt(i)
2869 fyi(i) = fyi(i) + fyt(i)
2870 fzi(i) = fzi(i) + fzt(i)
2871C---------------------------------
2872C CONTACT ENERGY CALCULATION
2873C---------------------------------
2874 efrict(i) = zero
2875 IF( intth > 0) THEN
2876 efrict(i) = dt1*(fxt(i)*vx(i) + fyt(i)*vy(i) + fzt(i)*vz(i)) ! FRICTIONAL ENERGY
2877 qfrict = qfrict + efrict(i)
2878 ENDIF
2879 efric_l(i)= dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
2880 econvt = econvt + efric_l(i)
2881 ENDDO
2882
2883
2884
2885 ELSE
2886C---------------------------------
2887C TOTAL (VISCOUS) FORMULATION / NO FRICTION FILTERING
2888C---------------------------------
2889#include "vectorize.inc"
2890 DO k=1,nfisot ! isotropic friction couples
2891 i = indexisot(k)
2892 vnx = n1(i)*vn(i)
2893 vny = n2(i)*vn(i)
2894 vnz = n3(i)*vn(i)
2895 vx(i) = vx(i) - vnx
2896 vy(i) = vy(i) - vny
2897 vz(i) = vz(i) - vnz
2898 v2 = vx(i)**2 + vy(i)**2 + vz(i)**2
2899 vis2(i) = viscffric(i) * vis2(i)
2900 fm2 = (xmu(i)*fni(i))**2
2901 f2 = vis2(i) * v2
2902 a2 = min(f2,fm2) / max(em30,f2)
2903 aa = sqrt(a2 * vis2(i))
2904 fxt(i) = aa * vx(i)
2905 fyt(i) = aa * vy(i)
2906 fzt(i) = aa * vz(i)
2907C------- total force
2908 fxi(i) = fxi(i) + fxt(i)
2909 fyi(i) = fyi(i) + fyt(i)
2910 fzi(i) = fzi(i) + fzt(i)
2911C---------------------------------
2912C CONTACT ENERGY CALCULATION
2913C---------------------------------
2914 efrict(i) = zero
2915 IF( intth > 0) THEN
2916 efrict(i) = aa * v2 * dt1 ! FRICTIONAL ENERGY
2917 qfrict = qfrict + efrict(i)
2918 ENDIF
2919 efric_l(i)= aa * v2 * dt1
2920 econvt = econvt + efric_l(i)
2921 ENDDO
2922
2923#include "vectorize.inc"
2924 DO k=1,nforth ! Orthotropic friction couples
2925 i = indexorth(k)
2926 vnx = n1(i)*vn(i)
2927 vny = n2(i)*vn(i)
2928 vnz = n3(i)*vn(i)
2929 vx(i) = vx(i) - vnx
2930 vy(i) = vy(i) - vny
2931 vz(i) = vz(i) - vnz
2932 vis2(i) = viscffric(i) * vis2(i)
2933
2934 vis = sqrt(vis2(i))
2935 fx = vis*vx(i)
2936 fy = vis*vy(i)
2937 fz = vis*vz(i)
2938
2939 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2940 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2941 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2942 ft = max(ft,em30)
2943 fn = fni(i)*fni(i)
2944
2945 beta = fn/ft
2946
2947 IF(beta == 0 ) THEN
2948 fxt(i) = zero
2949 fyt(i) = zero
2950 fzt(i) = zero
2951 ELSEIF(beta > 1) THEN ! Inside the ellipse
2952 fxt(i) = fx
2953 fyt(i) = fy
2954 fzt(i) = fz
2955 ELSE ! outside the ellipse
2956
2957! Projection on local tangent of ellipse (outside of ellipse)
2958! ANS = (Fadh-Fproj).n
2959 nep1 =ftt1*an(k)/fn
2960 nep2 =ftt2*bn(k)/fn
2961 nep =nep1*nep1+nep2*nep2
2962 nep =sqrt(nep)
2963
2964 ep=nep1*ftt1+nep2*ftt2
2965
2966 ans=(ep-sqrt(ep))/max(em20,nep)
2967 nep1 =nep1/max(em20,nep)
2968 nep2 =nep2/max(em20,nep)
2969
2970! Projection on ellipse
2971 c11 =ftt1-ans*nep1
2972 c22 =ftt2-ans*nep2
2973
2974 alphaf = atan(c22/c11)
2975
2976 signc = ftt1/max(em20,abs(ftt1))
2977 csa = signc*abs(cos(alphaf))
2978 signc = ftt2/max(em20,abs(ftt2))
2979 sna = signc*abs(sin(alphaf))
2980! Ft computation
2981 ft = fni(i) /sqrt( (csa*csa*an(k) + sna*sna*bn(k)))
2982 ftt1 = ft * csa
2983 ftt2 = ft * sna
2984
2985 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2986 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2987 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2988
2989 ENDIF
2990
2991
2992C------- total force
2993 fxi(i)=fxi(i) + fxt(i)
2994 fyi(i)=fyi(i) + fyt(i)
2995 fzi(i)=fzi(i) + fzt(i)
2996C---------------------------------
2997C CONTACT ENERGY CALCULATION
2998C---------------------------------
2999 efrict(i) = zero
3000 IF( intth > 0) THEN
3001 efrict(i) = aa * v2 * dt1 ! FRICTIONAL ENERGY
3002 qfrict = qfrict + efrict(i)
3003 ENDIF
3004 efric_l(i)= aa * v2 * dt1
3005 econvt = econvt + efric_l(i)
3006 ENDDO
3007
3008
3009 ENDIF
3010
3011 ENDIF
3012C---------------------------------
3013 IF(itied /= 0)THEN
3014 DO i=jlt-jlt_tied+1,jlt
3015 fxt(i)= t1x(i)*ft1(i) + t2x(i)*ft2(i)
3016 fyt(i)= t1y(i)*ft1(i) + t2y(i)*ft2(i)
3017 fzt(i)= t1z(i)*ft1(i) + t2z(i)*ft2(i)
3018C-------- total force
3019 fxi(i)=fxi(i) + fxt(i)
3020 fyi(i)=fyi(i) + fyt(i)
3021 fzi(i)=fzi(i) + fzt(i)
3022 END DO
3023 END IF
3024C---------------------------------
3025 IF((anim_v(12)+outp_v(12)+h3d_data%N_VECT_PCONT>0.AND.
3026 . ((tt>=output%TANIM .AND. tt<=output%TANIM_STOP).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
3027 . (manim>=4.AND.manim<=15).OR.h3d_data%MH3D/=0))
3028 . .OR.h3d_data%N_VECT_PCONT_MAX>0)THEN
3029 IF (inconv==1) THEN
3030#include "lockon.inc"
3031 DO i=1,jlt
3032 ftcont(1,ix1(i)) =ftcont(1,ix1(i)) + fxt(i)*h1(i)
3033 ftcont(2,ix1(i)) =ftcont(2,ix1(i)) + fyt(i)*h1(i)
3034 ftcont(3,ix1(i)) =ftcont(3,ix1(i)) + fzt(i)*h1(i)
3035 ftcont(1,ix2(i)) =ftcont(1,ix2(i)) + fxt(i)*h2(i)
3036 ftcont(2,ix2(i)) =ftcont(2,ix2(i)) + fyt(i)*h2(i)
3037 ftcont(3,ix2(i)) =ftcont(3,ix2(i)) + fzt(i)*h2(i)
3038 ftcont(1,ix3(i)) =ftcont(1,ix3(i)) + fxt(i)*h3(i)
3039 ftcont(2,ix3(i)) =ftcont(2,ix3(i)) + fyt(i)*h3(i)
3040 ftcont(3,ix3(i)) =ftcont(3,ix3(i)) + fzt(i)*h3(i)
3041 ftcont(1,ix4(i)) =ftcont(1,ix4(i)) + fxt(i)*h4(i)
3042 ftcont(2,ix4(i)) =ftcont(2,ix4(i)) + fyt(i)*h4(i)
3043 ftcont(3,ix4(i)) =ftcont(3,ix4(i)) + fzt(i)*h4(i)
3044 jg = nsvg(i)
3045 IF(jg>0) THEN
3046C In SPMD: Treatment to be redone after reception node Remote if JG <0
3047 ftcont(1,jg)=ftcont(1,jg)- fxt(i)
3048 ftcont(2,jg)=ftcont(2,jg)- fyt(i)
3049 ftcont(3,jg)=ftcont(3,jg)- fzt(i)
3050 ELSE ! cas noeud remote en SPMD
3051 jg = -jg
3052 ftconti(nin)%P(1,jg)=ftconti(nin)%P(1,jg)-fxt(i)
3053 ftconti(nin)%P(2,jg)=ftconti(nin)%P(2,jg)-fyt(i)
3054 ftconti(nin)%P(3,jg)=ftconti(nin)%P(3,jg)-fzt(i)
3055 ENDIF
3056 ENDDO
3057#include "lockoff.inc"
3058 END IF !(INCONV==1) THEN
3059 ENDIF
3060C
3061C---------------------------------
3062 fsav22= zero
3063 fsav23= zero
3064 fsav24= zero
3065 DO i=1,jlt
3066C if type7 sym of type19 opposite sign for tangential force -> forc_sign = -1
3067 impx=forc_sign*fxt(i)*dt12
3068 impy=forc_sign*fyt(i)*dt12
3069 impz=forc_sign*fzt(i)*dt12
3070 fsav4 =fsav4 +impx
3071 fsav5 =fsav5 +impy
3072 fsav6 =fsav6 +impz
3073C if type7 sym of type19 opposite sign for tangential force -> forc_sign = -1
3074 impx=forc_sign*fxi(i)*dt12
3075 impy=forc_sign*fyi(i)*dt12
3076 impz=forc_sign*fzi(i)*dt12
3077 fsav12=fsav12+abs(impx)
3078 fsav13=fsav13+abs(impy)
3079 fsav14=fsav14+abs(impz)
3080 fsav15=fsav15+sqrt(impx*impx+impy*impy+impz*impz)
3081 xp(i)=xi(i)+pene(i)*n1(i)
3082 yp(i)=yi(i)+pene(i)*n2(i)
3083 zp(i)=zi(i)+pene(i)*n3(i)
3084 fsav22=fsav22+yp(i)*impz-zp(i)*impy
3085 fsav23=fsav23+zp(i)*impx-xp(i)*impz
3086 fsav24=fsav24+xp(i)*impy-yp(i)*impx
3087 ENDDO
3088#include "lockon.inc"
3089 fsav(4) = fsav(4) + fsav4
3090 fsav(5) = fsav(5) + fsav5
3091 fsav(6) = fsav(6) + fsav6
3092 fsav(12) = fsav(12) + fsav12
3093 fsav(13) = fsav(13) + fsav13
3094 fsav(14) = fsav(14) + fsav14
3095 fsav(15) = fsav(15) + fsav15
3096 fsav(22) = fsav(22) + fsav22
3097 fsav(23) = fsav(23) + fsav23
3098 fsav(24) = fsav(24) + fsav24
3099 fsav(25) = fsav(25) + (fheats+fheatm)*qfrict
3100 fsav(26) = fsav(26) + econtt ! ECONTTIED part is not added in CE_ELAST for the moment
3101 fsav(27) = fsav(27) + econvt- (fheats+fheatm)*qfrict
3102 fsav(28) = fsav(28) + econtdt
3103#include "lockoff.inc"
3104C
3105 IF(isensint(1)/=0) THEN
3106 DO i=1,jlt
3107 fsavparit(1,4,i+nft) = forc_sign*fxt(i)
3108 fsavparit(1,5,i+nft) = forc_sign*fyt(i)
3109 fsavparit(1,6,i+nft) = forc_sign*fzt(i)
3110 ENDDO
3111 ENDIF
3112C---------------------------------
3113C SORTIES TH PAR SOUS INTERFACE
3114C---------------------------------
3115 IF(nisub/=0)THEN
3116 DO i=1,jlt
3117 nn = nsvg(i)
3118 IF(nn>0)THEN
3119 in=cn_loc(i)
3120 ie=ce_loc(i)
3121 jj =addsubs(in)
3122 kk =addsubm(ie)
3123 DO WHILE(jj<addsubs(in+1))
3124 jsub=lisubs(jj)
3125 itypsub = typsub(jsub)
3126
3127 IF(itypsub == 1 ) THEN ! Defining specific inter
3128
3129 ksub=lisubm(kk)
3130 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
3131 IF(ksub==jsub)THEN
3132C if type7 sym of type19 opposite sign for normal force -> forc_sign = -1
3133 impx=forc_sign*fxt(i)*dt12
3134 impy=forc_sign*fyt(i)*dt12
3135 impz=forc_sign*fzt(i)*dt12
3136C MAIN side :
3137 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
3138 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
3139 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
3140C
3141C if type7 sym of type19 opposite sign for normal force -> forc_sign = -1
3142 impx=forc_sign*fxi(i)*dt12
3143 impy=forc_sign*fyi(i)*dt12
3144 impz=forc_sign*fzi(i)*dt12
3145 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
3146 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
3147 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
3148C
3149 IF(isensint(jsub+1)/=0) THEN
3150 fsavparit(jsub+1,4,i+nft) = forc_sign*fxt(i)
3151 fsavparit(jsub+1,5,i+nft) = forc_sign*fyt(i)
3152 fsavparit(jsub+1,6,i+nft) = forc_sign*fzt(i)
3153 ENDIF
3154C
3155 fsavsub1(15,jsub)= fsavsub1(15,jsub)
3156 . +sqrt(impx*impx+impy*impy+impz*impz)
3157 fsavsub1(22,jsub)=fsavsub1(22,jsub)
3158 . +yp(i)*impz-zp(i)*impy
3159 fsavsub1(23,jsub)=fsavsub1(23,jsub)
3160 . +zp(i)*impx-xp(i)*impz
3161 fsavsub1(24,jsub)=fsavsub1(24,jsub)
3162 . +xp(i)*impy-yp(i)*impx
3163C
3164 END IF
3165
3166 kk=kk+1
3167 ksub=lisubm(kk)
3168 ENDDO
3169 jj=jj+1
3170
3171 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
3172
3173C if type7 sym of type19 opposite sign for normal force -> forc_sign = -1
3174 impx=forc_sign*fxt(i)*dt12
3175 impy=forc_sign*fyt(i)*dt12
3176 impz=forc_sign*fzt(i)*dt12
3177Cl MAIN side :
3178 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
3179 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
3180 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
3181C
3182C if type7 sym of type19 opposite sign for normal force -> forc_sign = -1
3183 impx=forc_sign*fxi(i)*dt12
3184 impy=forc_sign*fyi(i)*dt12
3185 impz=forc_sign*fzi(i)*dt12
3186 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
3187 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
3188 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
3189C
3190 IF(isensint(jsub+1)/=0) THEN
3191 fsavparit(jsub+1,4,i+nft) = forc_sign*fxt(i)
3192 fsavparit(jsub+1,5,i+nft) = forc_sign*fyt(i)
3193 fsavparit(jsub+1,6,i+nft) = forc_sign*fzt(i)
3194 ENDIF
3195C
3196 fsavsub1(15,jsub)= fsavsub1(15,jsub)
3197 . +sqrt(impx*impx+impy*impy+impz*impz)
3198 fsavsub1(22,jsub)=fsavsub1(22,jsub)
3199 . +yp(i)*impz-zp(i)*impy
3200 fsavsub1(23,jsub)=fsavsub1(23,jsub)
3201 . +zp(i)*impx-xp(i)*impz
3202 fsavsub1(24,jsub)=fsavsub1(24,jsub)
3203 . +xp(i)*impy-yp(i)*impx
3204
3205 jj=jj+1
3206
3207 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfaces
3208C
3209C Find if node is on Surface S1 -> S2 or S2 ->S1
3210 iss2 = bitget(inflg_subs(jj),0)
3211 iss1 = bitget(inflg_subs(jj),1)
3212 ksub=lisubm(kk)
3213 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
3214 ims2 = bitget(inflg_subm(kk),0)
3215 ims1 = bitget(inflg_subm(kk),1)
3216 IF(ksub==jsub)THEN
3217! S and M candidates on the same sub_interface
3218 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
3219 . (ims2 == 1 .AND. iss1 == 1))) THEN
3220 kk=kk+1
3221 ksub=lisubm(kk)
3222 cycle
3223 END IF
3224
3225 impx=forc_sign*fxt(i)*dt12
3226 impy=forc_sign*fyt(i)*dt12
3227 impz=forc_sign*fzt(i)*dt12
3228 IF(ims2 > 0) THEN
3229C main side :
3230 fsavsub1(4,jsub)=fsavsub1(4,jsub)-impx
3231 fsavsub1(5,jsub)=fsavsub1(5,jsub)-impy
3232 fsavsub1(6,jsub)=fsavsub1(6,jsub)-impz
3233 ELSE
3234 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
3235 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
3236 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
3237 ENDIF
3238C
3239 impx=forc_sign*fxi(i)*dt12
3240 impy=forc_sign*fyi(i)*dt12
3241 impz=forc_sign*fzi(i)*dt12
3242 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
3243 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
3244 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
3245C
3246 IF(isensint(jsub+1)/=0) THEN
3247 IF(ims2 > 0) THEN
3248 fsavparit(jsub+1,4,i+nft) = -forc_sign*fxt(i)
3249 fsavparit(jsub+1,5,i+nft) = -forc_sign*fyt(i)
3250 fsavparit(jsub+1,6,i+nft) = -forc_sign*fzt(i)
3251 ELSE
3252 fsavparit(jsub+1,4,i+nft) = forc_sign*fxt(i)
3253 fsavparit(jsub+1,5,i+nft) = forc_sign*fyt(i)
3254 fsavparit(jsub+1,6,i+nft) = forc_sign*fzt(i)
3255 ENDIF
3256 ENDIF
3257C
3258 fsavsub1(15,jsub)= fsavsub1(15,jsub)
3259 . +sqrt(impx*impx+impy*impy+impz*impz)
3260 fsavsub1(22,jsub)=fsavsub1(22,jsub)
3261 . +yp(i)*impz-zp(i)*impy
3262 fsavsub1(23,jsub)=fsavsub1(23,jsub)
3263 . +zp(i)*impx-xp(i)*impz
3264 fsavsub1(24,jsub)=fsavsub1(24,jsub)
3265 . +xp(i)*impy-yp(i)*impx
3266 ENDIF
3267 kk=kk+1
3268 ksub=lisubm(kk)
3269 ENDDO
3270 jj=jj+1
3271
3272 ENDIF
3273
3274 END DO
3275 END IF
3276
3277 ie=ce_loc(i)
3278
3279 kk =addsubm(ie) ! address of M
3280 DO WHILE(kk<addsubm(ie+1))
3281! all sub interfaces of S
3282 ksub=lisubm(kk)
3283! KSUB Croissant
3284 itypsub = typsub(ksub)
3285 IF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
3286 impx=-forc_sign*fxt(i)*dt12
3287 impy=-forc_sign*fyt(i)*dt12
3288 impz=-forc_sign*fzt(i)*dt12
3289C main side :
3290 fsavsub1(4,ksub)=fsavsub1(4,ksub)+impx
3291 fsavsub1(5,ksub)=fsavsub1(5,ksub)+impy
3292 fsavsub1(6,ksub)=fsavsub1(6,ksub)+impz
3293C
3294 impx=fxi(i)*dt12
3295 impy=fyi(i)*dt12
3296 impz=fzi(i)*dt12
3297 fsavsub1(12,ksub)=fsavsub1(12,ksub)+abs(impx)
3298 fsavsub1(13,ksub)=fsavsub1(13,ksub)+abs(impy)
3299 fsavsub1(14,ksub)=fsavsub1(14,ksub)+abs(impz)
3300C
3301 IF(isensint(ksub+1)/=0) THEN
3302 fsavparit(ksub+1,4,i+nft) = -forc_sign*fxt(i)
3303 fsavparit(ksub+1,5,i+nft) = -forc_sign*fyt(i)
3304 fsavparit(ksub+1,6,i+nft) = -forc_sign*fzt(i)
3305 ENDIF
3306C
3307 fsavsub1(15,ksub)= fsavsub1(15,ksub)
3308 . +sqrt(impx*impx+impy*impy+impz*impz)
3309 fsavsub1(22,ksub)=fsavsub1(22,ksub)
3310 . +yp(i)*impz-zp(i)*impy
3311 fsavsub1(23,ksub)=fsavsub1(23,ksub)
3312 . +zp(i)*impx-xp(i)*impz
3313 fsavsub1(24,ksub)=fsavsub1(24,ksub)
3314 . +xp(i)*impy-yp(i)*impx
3315
3316 ENDIF
3317 kk=kk+1
3318 ENDDO
3319 END DO
3320 IF(nspmd>1) THEN
3321C loop split because of a PGI bug
3322 DO i=1,jlt
3323 nn = nsvg(i)
3324 IF(nn<0)THEN
3325 nn = -nn
3326 ie=ce_loc(i)
3327 jj =addsubsfi(nin)%P(nn)
3328 kk =addsubm(ie)
3329 DO WHILE(jj<addsubsfi(nin)%P(nn+1))
3330 jsub=lisubsfi(nin)%P(jj)
3331 itypsub = typsub(jsub)
3332
3333 IF(itypsub == 1 ) THEN ! Defining specific inter
3334
3335 ksub=lisubm(kk)
3336 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
3337 IF(ksub==jsub)THEN
3338C if type7 sym of type19 opposite sign for normal force -> forc_sign = -1
3339 impx=forc_sign*fxt(i)*dt12
3340 impy=forc_sign*fyt(i)*dt12
3341 impz=forc_sign*fzt(i)*dt12
3342C MAIN side :
3343 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
3344 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
3345 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
3346C
3347C if type7 sym of type19 opposite sign for normal force -> forc_sign = -1
3348 impx=forc_sign*fxi(i)*dt12
3349 impy=forc_sign*fyi(i)*dt12
3350 impz=forc_sign*fzi(i)*dt12
3351 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
3352 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
3353 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
3354C
3355 IF(isensint(jsub+1)/=0) THEN
3356 fsavparit(jsub+1,4,i+nft) = forc_sign*fxt(i)
3357 fsavparit(jsub+1,5,i+nft) = forc_sign*fyt(i)
3358 fsavparit(jsub+1,6,i+nft) = forc_sign*fzt(i)
3359 ENDIF
3360C
3361 fsavsub1(15,jsub)= fsavsub1(15,jsub)
3362 . +sqrt(impx*impx+impy*impy+impz*impz)
3363 fsavsub1(22,jsub)=fsavsub1(22,jsub)
3364 . +yp(i)*impz-zp(i)*impy
3365 fsavsub1(23,jsub)=fsavsub1(23,jsub)
3366 . +zp(i)*impx-xp(i)*impz
3367 fsavsub1(24,jsub)=fsavsub1(24,jsub)
3368 . +xp(i)*impy-yp(i)*impx
3369C
3370 END IF
3371
3372 kk=kk+1
3373 ksub=lisubm(kk)
3374 ENDDO
3375 jj=jj+1
3376
3377 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
3378
3379C if type7 sym of type19 opposite sign for normal force -> forc_sign = -1
3380 impx=forc_sign*fxt(i)*dt12
3381 impy=forc_sign*fyt(i)*dt12
3382 impz=forc_sign*fzt(i)*dt12
3383Cl MAIN side :
3384 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
3385 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
3386 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
3387C
3388C if type7 sym of type19 opposite sign for normal force -> forc_sign = -1
3389 impx=forc_sign*fxi(i)*dt12
3390 impy=forc_sign*fyi(i)*dt12
3391 impz=forc_sign*fzi(i)*dt12
3392 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
3393 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
3394 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
3395C
3396 IF(isensint(jsub+1)/=0) THEN
3397 fsavparit(jsub+1,4,i+nft) = forc_sign*fxt(i)
3398 fsavparit(jsub+1,5,i+nft) = forc_sign*fyt(i)
3399 fsavparit(jsub+1,6,i+nft) = forc_sign*fzt(i)
3400 ENDIF
3401C
3402 fsavsub1(15,jsub)= fsavsub1(15,jsub)
3403 . +sqrt(impx*impx+impy*impy+impz*impz)
3404 fsavsub1(22,jsub)=fsavsub1(22,jsub)
3405 . +yp(i)*impz-zp(i)*impy
3406 fsavsub1(23,jsub)=fsavsub1(23,jsub)
3407 . +zp(i)*impx-xp(i)*impz
3408 fsavsub1(24,jsub)=fsavsub1(24,jsub)
3409 . +xp(i)*impy-yp(i)*impx
3410
3411 jj=jj+1
3412 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfaces
3413
3414 iss2 = bitget(inflg_subsfi(nin)%P(jj),0)
3415 iss1 = bitget(inflg_subsfi(nin)%P(jj),1)
3416 ksub=lisubm(kk)
3417 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
3418 ims2 = bitget(inflg_subm(kk),0)
3419 ims1 = bitget(inflg_subm(kk),1)
3420 IF(ksub==jsub)THEN
3421 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
3422 . (ims2 == 1 .AND. iss1 == 1))) THEN
3423 kk=kk+1
3424 ksub=lisubm(kk)
3425 cycle
3426 END IF
3427
3428 impx=forc_sign*fxt(i)*dt12
3429 impy=forc_sign*fyt(i)*dt12
3430 impz=forc_sign*fzt(i)*dt12
3431 IF(ims2 > 0) THEN
3432C main side :
3433 fsavsub1(4,jsub)=fsavsub1(4,jsub)-impx
3434 fsavsub1(5,jsub)=fsavsub1(5,jsub)-impy
3435 fsavsub1(6,jsub)=fsavsub1(6,jsub)-impz
3436 ELSE
3437 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
3438 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
3439 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
3440 ENDIF
3441C
3442 impx=fxi(i)*dt12
3443 impy=fyi(i)*dt12
3444 impz=fzi(i)*dt12
3445 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
3446 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
3447 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
3448C
3449 IF(isensint(jsub+1)/=0) THEN
3450 IF(ims2 > 0) THEN
3451 fsavparit(jsub+1,4,i) = forc_sign*fxt(i)
3452 fsavparit(jsub+1,5,i) = forc_sign*fyt(i)
3453 fsavparit(jsub+1,6,i) = forc_sign*fzt(i)
3454 ELSE
3455 fsavparit(jsub+1,4,i) = -forc_sign*fxt(i)
3456 fsavparit(jsub+1,5,i) = -forc_sign*fyt(i)
3457 fsavparit(jsub+1,6,i) = -forc_sign*fzt(i)
3458 ENDIF
3459 ENDIF
3460C
3461 fsavsub1(15,jsub)= fsavsub1(15,jsub)
3462 . +sqrt(impx*impx+impy*impy+impz*impz)
3463 fsavsub1(22,jsub)=fsavsub1(22,jsub)
3464 . +yp(i)*impz-zp(i)*impy
3465 fsavsub1(23,jsub)=fsavsub1(23,jsub)
3466 . +zp(i)*impx-xp(i)*impz
3467 fsavsub1(24,jsub)=fsavsub1(24,jsub)
3468 . +xp(i)*impy-yp(i)*impx
3469 ENDIF
3470 kk=kk+1
3471 ksub=lisubm(kk)
3472 ENDDO
3473 jj=jj+1
3474
3475 ENDIF
3476 END DO
3477 END IF
3478 END DO
3479 END IF
3480#include "lockon.inc"
3481 DO jsub=1,nisub
3482 nsub=lisub(jsub)
3483 DO j=1,15
3484 fsavsub(j,nsub)=fsavsub(j,nsub)+fsavsub1(j,jsub)
3485 END DO
3486 fsavsub(22,nsub)=fsavsub(22,nsub)+fsavsub1(22,jsub)
3487 fsavsub(23,nsub)=fsavsub(23,nsub)+fsavsub1(23,jsub)
3488 fsavsub(24,nsub)=fsavsub(24,nsub)+fsavsub1(24,jsub)
3489 END DO
3490#include "lockoff.inc"
3491 END IF
3492C---------------------------------
3493#include "lockon.inc"
3494 IF (inconv==1) THEN
3495 econtv = econtv + econvt ! Frictional Energy
3496 econt = econt + econtt ! Elastic Energy
3497 econtd = econtd + econtdt ! Damping Energy
3498 econt_cumu = econt_cumu + econttied ! Elastic energy for tied contact : it is cumulated energy
3499 END IF !(INCONV==1) THEN
3500 IF (intth/=0) THEN
3501 qfric = qfric + (fheats+fheatm)*qfrict ! FRICTIONAL HEAT ADDED TO INTERNAL ENERGY
3502 econtv = econtv - (fheats+fheatm)*qfrict ! FRICTIONAL HEAT REMOVED FROM CONTACT ENERGY
3503 ENDIF
3504#include "lockoff.inc"
3505C---------------------------------
3506 IF(interefric >0)THEN
3507 IF (inconv==1) THEN
3508#include "lockon.inc"
3509 DO i=1,jlt
3510 efric_lm = half*efric_l(i)
3511 output%DATA%EFRIC(interefric,ix1(i)) = output%DATA%EFRIC(interefric,ix1(i)) + h1(i)*(efric_lm-fheatm*efrict(i))
3512 output%DATA%EFRIC(interefric,ix2(i)) = output%DATA%EFRIC(interefric,ix2(i)) + h2(i)*(efric_lm-fheatm*efrict(i))
3513 output%DATA%EFRIC(interefric,ix3(i)) = output%DATA%EFRIC(interefric,ix3(i)) + h3(i)*(efric_lm-fheatm*efrict(i))
3514 output%DATA%EFRIC(interefric,ix4(i)) = output%DATA%EFRIC(interefric,ix4(i)) + h4(i)*(efric_lm-fheatm*efrict(i))
3515 jg = nsvg(i)
3516 efric_ls = half*efric_l(i)
3517 IF(jg>0) THEN
3518 output%DATA%EFRIC(interefric,jg) = output%DATA%EFRIC(interefric,jg) + (efric_ls-fheats*efrict(i))
3519 ELSE ! node remote
3520 jg = -jg
3521 efricfi(nin)%P(jg)=efricfi(nin)%P(jg)+ (efric_ls-fheats*efrict(i))
3522 ENDIF
3523 ENDDO
3524#include "lockoff.inc"
3525 END IF !(INCONV==1) THEN
3526 ENDIF
3527C---------------------------------
3528 IF(h3d_data%N_SCAL_CSE_FRIC >0)THEN
3529 IF (inconv==1) THEN
3530#include "lockon.inc"
3531 DO i=1,jlt
3532 efric_lm = half*efric_l(i)
3533 output%DATA%EFRICG(ix1(i)) = output%DATA%EFRICG(ix1(i)) + h1(i)*(efric_lm-fheatm*efrict(i))
3534 output%DATA%EFRICG(ix2(i)) = output%DATA%EFRICG(ix2(i)) + h2(i)*(efric_lm-fheatm*efrict(i))
3535 output%DATA%EFRICG(ix3(i)) = output%DATA%EFRICG(ix3(i)) + h3(i)*(efric_lm-fheatm*efrict(i))
3536 output%DATA%EFRICG(ix4(i)) = output%DATA%EFRICG(ix4(i)) + h4(i)*(efric_lm-fheatm*efrict(i))
3537 jg = nsvg(i)
3538 efric_ls = half*efric_l(i)
3539 IF(jg>0) THEN
3540 output%DATA%EFRICG(jg) = output%DATA%EFRICG(jg) + (efric_ls-fheats*efrict(i))
3541 ELSE ! node remote
3542 jg = -jg
3543 efricgfi(nin)%P(jg)=efricgfi(nin)%P(jg)+ (efric_ls-fheats*efrict(i))
3544 ENDIF
3545 ENDDO
3546#include "lockoff.inc"
3547 END IF !(INCONV==1) THEN
3548 ENDIF
3549C---------------------------------
3550 IF(kdtint==1)THEN
3551 IF( (visc/=zero)
3552 . .AND.(ivis2==0.OR.ivis2==1))THEN
3553 DO i=1,jlt
3554C C (i) = 2.*C (i)
3555C
3556 IF(msi(i)==zero)THEN
3557 ks(i) =zero
3558 cs(i) =zero
3559 stv(i)=zero
3560 ELSE
3561 cx = four*c(i)*c(i)
3562 cy = eight*msi(i)*kt(i)
3563 aux = sqrt(cx+cy)+two*c(i)
3564 stv(i)= kt(i)*aux*aux/max(cy,em30)
3565 aux = two*cf(i)*cf(i)/max(msi(i),em20)
3566 IF(aux>stv(i))THEN
3567 ks(i) =zero
3568 cs(i) =cf(i)
3569 stv(i)=aux
3570 ELSE
3571 ks(i)= kt(i)
3572 cs(i) =c(i)
3573 ENDIF
3574 ENDIF
3575C
3576 j1=ix1(i)
3577 IF(ms(j1)==zero)THEN
3578 k1(i) =zero
3579 c1(i) =zero
3580 st1(i)=zero
3581 ELSE
3582 k1(i)=kt(i)*abs(h1(i))
3583 c1(i)=c(i)*abs(h1(i))
3584 cx =four*c1(i)*c1(i)
3585 cy =eight*ms(j1)*k1(i)
3586 aux = sqrt(cx+cy)+two*c1(i)
3587 st1(i)= k1(i)*aux*aux/max(cy,em30)
3588 cfi = cf(i)*abs(h1(i))
3589 aux = two*cfi*cfi/max(ms(j1),em20)
3590 IF(aux>st1(i))THEN
3591 k1(i) =zero
3592 c1(i) =cfi
3593 st1(i)=aux
3594 ENDIF
3595 ENDIF
3596C
3597 j1=ix2(i)
3598 IF(ms(j1)==zero)THEN
3599 k2(i) =zero
3600 c2(i) =zero
3601 st2(i)=zero
3602 ELSE
3603 k2(i)=kt(i)*abs(h2(i))
3604 c2(i)=c(i)*abs(h2(i))
3605 cx =four*c2(i)*c2(i)
3606 cy =eight*ms(j1)*k2(i)
3607 aux = sqrt(cx+cy)+two*c2(i)
3608 st2(i)= k2(i)*aux*aux/max(cy,em30)
3609 cfi = cf(i)*abs(h2(i))
3610 aux = two*cfi*cfi/max(ms(j1),em20)
3611 IF(aux>st2(i))THEN
3612 k2(i) =zero
3613 c2(i) =cfi
3614 st2(i)=aux
3615 ENDIF
3616 ENDIF
3617C
3618 j1=ix3(i)
3619 IF(ms(j1)==zero)THEN
3620 k3(i) =zero
3621 c3(i) =zero
3622 st3(i)=zero
3623 ELSE
3624 k3(i)=kt(i)*abs(h3(i))
3625 c3(i)=c(i)*abs(h3(i))
3626 cx =four*c3(i)*c3(i)
3627 cy =eight*ms(j1)*k3(i)
3628 aux = sqrt(cx+cy)+two*c3(i)
3629 st3(i)= k3(i)*aux*aux/max(cy,em30)
3630 cfi = cf(i)*abs(h3(i))
3631 aux = two*cfi*cfi/max(ms(j1),em20)
3632 IF(aux>st3(i))THEN
3633 k3(i) =zero
3634 c3(i) =cfi
3635 st3(i)=aux
3636 ENDIF
3637 ENDIF
3638C
3639 j1=ix4(i)
3640 IF(ms(j1)==zero)THEN
3641 k4(i) =zero
3642 c4(i) =zero
3643 st4(i)=zero
3644 ELSE
3645 k4(i)=kt(i)*abs(h4(i))
3646 c4(i)=c(i)*abs(h4(i))
3647 cx =four*c4(i)*c4(i)
3648 cy =eight*ms(j1)*k4(i)
3649 aux = sqrt(cx+cy)+two*c4(i)
3650 st4(i)= k4(i)*aux*aux/max(cy,em30)
3651 cfi = cf(i)*abs(h4(i))
3652 aux = two*cfi*cfi/max(ms(j1),em20)
3653 IF(aux>st4(i))THEN
3654 k4(i) =zero
3655 c4(i) =cfi
3656 st4(i)=aux
3657 ENDIF
3658 ENDIF
3659 ENDDO
3660C
3661 ELSE
3662 DO i=1,jlt
3663 IF( (viscffric(i)/=zero)
3664 . .AND.(ivis2==0.OR.ivis2==1))THEN
3665 IF(msi(i)==zero)THEN
3666 ks(i) =zero
3667 cs(i) =zero
3668 stv(i)=zero
3669 ELSE
3670 cx = four*c(i)*c(i)
3671 cy = eight*msi(i)*kt(i)
3672 aux = sqrt(cx+cy)+two*c(i)
3673 stv(i)= kt(i)*aux*aux/max(cy,em30)
3674 aux = two*cf(i)*cf(i)/max(msi(i),em20)
3675 IF(aux>stv(i))THEN
3676 ks(i) =zero
3677 cs(i) =cf(i)
3678 stv(i)=aux
3679 ELSE
3680 ks(i)= kt(i)
3681 cs(i) =c(i)
3682 ENDIF
3683 ENDIF
3684C
3685 j1=ix1(i)
3686 IF(ms(j1)==zero)THEN
3687 k1(i) =zero
3688 c1(i) =zero
3689 st1(i)=zero
3690 ELSE
3691 k1(i)=kt(i)*abs(h1(i))
3692 c1(i)=c(i)*abs(h1(i))
3693 cx =four*c1(i)*c1(i)
3694 cy =eight*ms(j1)*k1(i)
3695 aux = sqrt(cx+cy)+two*c1(i)
3696 st1(i)= k1(i)*aux*aux/max(cy,em30)
3697 cfi = cf(i)*abs(h1(i))
3698 aux = two*cfi*cfi/max(ms(j1),em20)
3699 IF(aux>st1(i))THEN
3700 k1(i) =zero
3701 c1(i) =cfi
3702 st1(i)=aux
3703 ENDIF
3704 ENDIF
3705C
3706 j1=ix2(i)
3707 IF(ms(j1)==zero)THEN
3708 k2(i) =zero
3709 c2(i) =zero
3710 st2(i)=zero
3711 ELSE
3712 k2(i)=kt(i)*abs(h2(i))
3713 c2(i)=c(i)*abs(h2(i))
3714 cx =four*c2(i)*c2(i)
3715 cy =eight*ms(j1)*k2(i)
3716 aux = sqrt(cx+cy)+two*c2(i)
3717 st2(i)= k2(i)*aux*aux/max(cy,em30)
3718 cfi = cf(i)*abs(h2(i))
3719 aux = two*cfi*cfi/max(ms(j1),em20)
3720 IF(aux>st2(i))THEN
3721 k2(i) =zero
3722 c2(i) =cfi
3723 st2(i)=aux
3724 ENDIF
3725 ENDIF
3726C
3727 j1=ix3(i)
3728 IF(ms(j1)==zero)THEN
3729 k3(i) =zero
3730 c3(i) =zero
3731 st3(i)=zero
3732 ELSE
3733 k3(i)=kt(i)*abs(h3(i))
3734 c3(i)=c(i)*abs(h3(i))
3735 cx =four*c3(i)*c3(i)
3736 cy =eight*ms(j1)*k3(i)
3737 aux = sqrt(cx+cy)+two*c3(i)
3738 st3(i)= k3(i)*aux*aux/max(cy,em30)
3739 cfi = cf(i)*abs(h3(i))
3740 aux = two*cfi*cfi/max(ms(j1),em20)
3741 IF(aux>st3(i))THEN
3742 k3(i) =zero
3743 c3(i) =cfi
3744 st3(i)=aux
3745 ENDIF
3746 ENDIF
3747C
3748 j1=ix4(i)
3749 IF(ms(j1)==zero)THEN
3750 k4(i) =zero
3751 c4(i) =zero
3752 st4(i)=zero
3753 ELSE
3754 k4(i)=kt(i)*abs(h4(i))
3755 c4(i)=c(i)*abs(h4(i))
3756 cx =four*c4(i)*c4(i)
3757 cy =eight*ms(j1)*k4(i)
3758 aux = sqrt(cx+cy)+two*c4(i)
3759 st4(i)= k4(i)*aux*aux/max(cy,em30)
3760 cfi = cf(i)*abs(h4(i))
3761 aux = two*cfi*cfi/max(ms(j1),em20)
3762 IF(aux>st4(i))THEN
3763 k4(i) =zero
3764 c4(i) =cfi
3765 st4(i)=aux
3766 ENDIF
3767 ENDIF
3768 ELSE
3769 ks(i) =stif(i)
3770 cs(i) =zero
3771 stv(i)=ks(i)
3772 k1(i) =stif(i)*abs(h1(i))
3773 c1(i) =zero
3774 st1(i)=k1(i)
3775 k2(i) =stif(i)*abs(h2(i))
3776 c2(i) =zero
3777 st2(i)=k2(i)
3778 k3(i) =stif(i)*abs(h3(i))
3779 c3(i) =zero
3780 st3(i)=k3(i)
3781 k4(i) =stif(i)*abs(h4(i))
3782 c4(i) =zero
3783 st4(i)=k4(i)
3784 ENDIF
3785 ENDDO
3786 ENDIF
3787 ENDIF
3788C
3789C---------------------
3790 IF(itied /= 0)THEN
3791C
3792 dti = 1.e20
3793C
3794 DO i=1,jlt-jlt_tied
3795 dist = gapv(i)-pene(i)
3796 rdist = half*dist / max(em30,-vn(i))
3797 dti = min(rdist,dti)
3798 END DO
3799C
3800 IF (dtmini>zero) THEN
3801 dtm =dtmini
3802 ELSE
3803 dtm =em30
3804 END IF
3805
3806 IF(dti<=dtm)THEN
3807 dti = 1.e20
3808 DO i=1,jlt-jlt_tied
3809 dist =gapv(i)-pene(i)
3810 dti2 = half*dist / max(em30,-vn(i))
3811 IF(dti2<=dtm.AND.cand_f(1,index(i))==zero)THEN
3812 nn = nsvg(i)
3813 IF(nn>0)THEN
3814 ni = itab(nn)
3815 ELSE
3816 ni = itafi(nin)%P(-nn)
3817 ENDIF
3818#include "lockon.inc"
3819 WRITE(iout,'(A,E12.4,A,I10,A,E12.4,A)')
3820 . ' **WARNING MINIMUM TIME STEP ',dti2,
3821 . ' IN INTERFACE ',noint,'(DTMIN=',dtm,')'
3822 WRITE(iout,'(A,I10,A,I10)')' TYING SECONDARY NODE VS MAIN SEGMENT',
3823 . ni,' IN INTERFACE ',noint
3824 WRITE(iout,'(A,4I10)')' MAIN NODES : ',
3825 . itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i))
3826#include "lockoff.inc"
3827 cand_f(1,index(i))= fns(i)*sign(one,side(i)) ! negative vs NX NY NZ
3828 cand_f(2,index(i))= fxt(i)*t1x(i)+fyt(i)*t1y(i)+fzt(i)*t1z(i)
3829 cand_f(3,index(i))= fxt(i)*t2x(i)+fyt(i)*t2y(i)+fzt(i)*t2z(i)
3830 cand_f(4,index(i))= h1(i)
3831 cand_f(5,index(i))= h2(i)
3832 cand_f(6,index(i))= h3(i)
3833 cand_f(7,index(i))= pene(i)
3834 cand_f(8,index(i))= side(i)
3835
3836 ENDIF
3837 ENDDO
3838 ENDIF
3839
3840 ENDIF ! IF(ITIED /= 0)THEN
3841
3842C--------------------------------------------
3843 IF(itied == 0)THEN
3844
3845 IF(idtm==1.OR.idtm==2.OR.
3846 . idtm==5.OR.idtm==6)THEN
3847 dtmi0 = ep20
3848 IF(kdtint==0)THEN
3849 DO i=1,jlt
3850 dtmi(i) = ep20
3851 mas2 = two * msi(i)
3852 IF(mas2>zero.AND.stif(i)>zero.AND.
3853 . irb(kini(i))==0.AND.irb2(kini(i))==0)THEN
3854 dtmi(i) = min(dtmi(i),dtfac1(10)*sqrt(mas2/stif(i)))
3855 ENDIF
3856 mas2 = two* ms(ix1(i))
3857 IF(mas2>zero.AND.h1(i)*stif(i)>zero.AND.
3858 . irb(kinet(ix1(i)))==0.AND.irb2(kinet(ix1(i)))==0)THEN
3859 dtmi(i) = min(dtmi(i),dtfac1(10)*sqrt(mas2/(h1(i)*stif(i))))
3860 ENDIF
3861 mas2 = two * ms(ix2(i))
3862 IF(mas2>zero.AND.h2(i)*stif(i)>zero.AND.
3863 . irb(kinet(ix2(i)))==0.AND.irb2(kinet(ix2(i)))==0)THEN
3864 dtmi(i) = min(dtmi(i),dtfac1(10)*sqrt(mas2/(h2(i)*stif(i))))
3865 ENDIF
3866 mas2 = two* ms(ix3(i))
3867 IF(mas2>zero.AND.h3(i)*stif(i)>zero.AND.
3868 . irb(kinet(ix3(i)))==0.AND.irb2(kinet(ix3(i)))==0)THEN
3869 dtmi(i) = min(dtmi(i),dtfac1(10)*sqrt(mas2/(h3(i)*stif(i))))
3870 ENDIF
3871 mas2 = two * ms(ix4(i))
3872 IF(mas2>zero.AND.h4(i)*stif(i)>zero.AND.
3873 . irb(kinet(ix4(i)))==0.AND.irb2(kinet(ix4(i)))==0)THEN
3874 dtmi(i) = min(dtmi(i),dtfac1(10)*sqrt(mas2/(h4(i)*stif(i))))
3875 ENDIF
3876 dtmi0 = min(dtmi0,dtmi(i))
3877 ENDDO
3878 ELSE
3879 DO i=1,jlt
3880 dtmi(i) = ep20
3881 mas2 = two * msi(i)
3882 IF(mas2>zero.AND.stv(i)>zero.AND.
3883 . irb(kini(i))==0.AND.irb2(kini(i))==0)THEN
3884 dtmi(i) = min(dtmi(i),dtfac1(10)*sqrt(mas2/stv(i)))
3885 ENDIF
3886 mas2 = two * ms(ix1(i))
3887 IF(mas2>zero.AND.st1(i)>zero.AND.
3888 . irb(kinet(ix1(i)))==0.AND.irb2(kinet(ix1(i)))==0)THEN
3889 dtmi(i) = min(dtmi(i),dtfac1(10)*sqrt(mas2/(st1(i))))
3890 ENDIF
3891 mas2 = two * ms(ix2(i))
3892 IF(mas2>zero.AND.st2(i)>zero.AND.
3893 . irb(kinet(ix2(i)))==0.AND.irb2(kinet(ix2(i)))==0)THEN
3894 dtmi(i) = min(dtmi(i),dtfac1(10)*sqrt(mas2/(st2(i))))
3895 ENDIF
3896 mas2 = two * ms(ix3(i))
3897 IF(mas2>zero.AND.st3(i)>zero.AND.
3898 . irb(kinet(ix3(i)))==0.AND.irb2(kinet(ix3(i)))==0)THEN
3899 dtmi(i) = min(dtmi(i),dtfac1(10)*sqrt(mas2/(st3(i))))
3900 ENDIF
3901 mas2 = two * ms(ix4(i))
3902 IF(mas2>zero.AND.st4(i)>zero.AND.
3903 . irb(kinet(ix4(i)))==0.AND.irb2(kinet(ix4(i)))==0)THEN
3904 dtmi(i) = min(dtmi(i),dtfac1(10)*sqrt(mas2/(st4(i))))
3905 ENDIF
3906 dtmi0 = min(dtmi0,dtmi(i))
3907 ENDDO
3908 ENDIF
3909 IF(dtmi0<=dtm)THEN
3910 DO i=1,jlt
3911 IF(dtmi(i)<=dtm)THEN
3912 jg = nsvg(i)
3913 IF(jg>0)THEN
3914 ni = itab(jg)
3915 ELSE
3916 ni = itafi(nin)%P(-jg)
3917 ENDIF
3918 IF(idtm==1)THEN
3919#include "lockon.inc"
3920 WRITE(iout,'(A,E12.4,A,I10,A,E12.4,A)')
3921 . ' **WARNING MINIMUM TIME STEP ',dtmi(i),
3922 . ' IN INTERFACE ',noint,'(DTMIN=',dtm,')'
3923 WRITE(iout,'(A,I10)') ' SECONDARY NODE : ',ni
3924 WRITE(iout,'(A,4I10)')' MAIN NODES : ',
3925 . itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i))
3926#include "lockoff.inc"
3927 tstop = tt
3928 IF ( istamping == 1) THEN
3929 WRITE(istdo,'(A)')'The run encountered a problem in an in
3930 .terface Type 7.'
3931 WRITE(istdo,'(A)')'You may need to check if there is enou
3932 .gh clearance between the tools,'
3933 WRITE(istdo,'(A)')'and that they do not penetrate each ot
3934 .her during their travel'
3935 WRITE(iout, '(A)')'The run encountered a problem in an in
3936 .terface Type 7.'
3937 WRITE(iout, '(A)')'You may need to check if there is enou
3938 .gh clearance between the tools,'
3939 WRITE(iout, '(A)')'and that they do not penetrate each ot
3940 .her during their travel'
3941 ENDIF
3942 ELSEIF(idtm==2)THEN
3943#include "lockon.inc"
3944 WRITE(iout,'(A,E12.4,A,I10,A,E12.4,A)')
3945 . ' **WARNING MINIMUM TIME STEP ',dtmi(i),
3946 . ' IN INTERFACE ',noint,'(DTMIN=',dtm,')'
3947 WRITE(iout,'(A,I10,A,I10)')' DELETE SECONDARY NODE ',
3948 . ni,' FROM INTERFACE ',noint
3949 WRITE(iout,'(A,4I10)')' MAIN NODES : ',
3950 . itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i))
3951 IF(jg>0) THEN
3952 stfn(cn_loc(i)) = -abs(stfn(cn_loc(i)))
3953 ELSE
3954 stifi(nin)%P(-jg) = -abs(stifi(nin)%P(-jg))
3955 ENDIF
3956 IF ( istamping == 1) THEN
3957 WRITE(istdo,'(A)')'The run encountered a problem in an in
3958 .terface Type 7.'
3959 WRITE(istdo,'(A)')'You may need to check if there is enou
3960 .gh clearance between the tools,'
3961 WRITE(istdo,'(A)')'and that they do not penetrate each ot
3962 .her during their travel'
3963 WRITE(iout, '(A)')'The run encountered a problem in an in
3964 .terface Type 7.'
3965 WRITE(iout, '(A)')'You may need to check if there is enou
3966 .gh clearance between the tools,'
3967 WRITE(iout, '(A)')'and that they do not penetrate each ot
3968 .her during their travel'
3969 ENDIF
3970#include "lockoff.inc"
3971 newfront = -1
3972 ELSEIF(idtm==5)THEN
3973#include "lockon.inc"
3974 WRITE(iout,'(A,E12.4,A,I10,A,E12.4,A)')
3975 . ' **WARNING MINIMUM TIME STEP ',dtmi(i),
3976 . ' IN INTERFACE ',noint,'(DTMIN=',dtm,')'
3977 WRITE(iout,'(A,I10)') ' SECONDARY NODE : ',ni
3978 WRITE(iout,'(A,4I10)')' MAIN NODES : ',
3979 . itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i))
3980#include "lockoff.inc"
3981 mstop = 2
3982 IF ( istamping == 1) THEN
3983 WRITE(istdo,'(A)')'The run encountered a problem in an in
3984 .terface Type 7.'
3985 WRITE(istdo,'(A)')'You may need to check if there is enou
3986 .gh clearance between the tools,'
3987 WRITE(istdo,'(A)')'and that they do not penetrate each ot
3988 .her during their travel'
3989 WRITE(iout, '(A)')'The run encountered a problem in an in
3990 .terface Type 7.'
3991 WRITE(iout, '(A)')'You may need to check if there is enou
3992 .gh clearance between the tools,'
3993 WRITE(iout, '(A)')'and that they do not penetrate each ot
3994 .her during their travel'
3995 ENDIF
3996 ELSEIF(idtm==6.AND.ilagm==2)THEN
3997 IF(kinet(jg)+kinet(ix1(i))+kinet(ix2(i))
3998 . +kinet(ix3(i))+kinet(ix4(i))==0 )THEN
3999 cand_n(index(i)) = -iabs(cand_n(index(i)))
4000 stif(i) = 0.
4001 fxi(i) = 0.
4002 fyi(i) = 0.
4003 fzi(i) = 0.
4004 ENDIF
4005 ENDIF
4006 ENDIF
4007 ENDDO
4008 ENDIF
4009 ENDIF
4010
4011 ELSE ! IF(ITIED == 0)THEN
4012
4013 IF (dtmini>zero) THEN
4014 dtm=dtmini
4015 ELSE
4016 dtm=em30
4017 END IF
4018
4019 dtmi0 = ep20
4020 IF(kdtint==0)THEN
4021 DO i=1,jlt-jlt_tied
4022 dtmi(i) = ep20
4023 mas2 = two * msi(i)
4024 IF(mas2>zero.AND.stif(i)>zero.AND.
4025 . irb(kini(i))==0.AND.irb2(kini(i))==0)THEN
4026 dtmi(i) = min(dtmi(i),sqrt(mas2/stif(i)))
4027 ENDIF
4028 mas2 = two* ms(ix1(i))
4029 IF(mas2>zero.AND.h1(i)*stif(i)>zero.AND.
4030 . irb(kinet(ix1(i)))==0.AND.irb2(kinet(ix1(i)))==0)THEN
4031 dtmi(i) = min(dtmi(i),sqrt(mas2/(h1(i)*stif(i))))
4032 ENDIF
4033 mas2 = two * ms(ix2(i))
4034 IF(mas2>zero.AND.h2(i)*stif(i)>zero.AND.
4035 . irb(kinet(ix2(i)))==0.AND.irb2(kinet(ix2(i)))==0)THEN
4036 dtmi(i) = min(dtmi(i),sqrt(mas2/(h2(i)*stif(i))))
4037 ENDIF
4038 mas2 = two* ms(ix3(i))
4039 IF(mas2>zero.AND.h3(i)*stif(i)>zero.AND.
4040 . irb(kinet(ix3(i)))==0.AND.irb2(kinet(ix3(i)))==0)THEN
4041 dtmi(i) = min(dtmi(i),sqrt(mas2/(h3(i)*stif(i))))
4042 ENDIF
4043 mas2 = two * ms(ix4(i))
4044 IF(mas2>zero.AND.h4(i)*stif(i)>zero.AND.
4045 . irb(kinet(ix4(i)))==0.AND.irb2(kinet(ix4(i)))==0)THEN
4046 dtmi(i) = min(dtmi(i),sqrt(mas2/(h4(i)*stif(i))))
4047 ENDIF
4048 dtmi0 = min(dtmi0,dtmi(i))
4049 ENDDO
4050 ELSE
4051 DO i=1,jlt-jlt_tied
4052 dtmi(i) = ep20
4053 mas2 = two * msi(i)
4054 IF(mas2>zero.AND.stv(i)>zero.AND.
4055 . irb(kini(i))==0.AND.irb2(kini(i))==0)THEN
4056 dtmi(i) = min(dtmi(i),sqrt(mas2/stv(i)))
4057 ENDIF
4058 mas2 = two * ms(ix1(i))
4059 IF(mas2>zero.AND.st1(i)>zero.AND.
4060 . irb(kinet(ix1(i)))==0.AND.irb2(kinet(ix1(i)))==0)THEN
4061 dtmi(i) = min(dtmi(i),sqrt(mas2/(st1(i))))
4062 ENDIF
4063 mas2 = two * ms(ix2(i))
4064 IF(mas2>zero.AND.st2(i)>zero.AND.
4065 . irb(kinet(ix2(i)))==0.AND.irb2(kinet(ix2(i)))==0)THEN
4066 dtmi(i) = min(dtmi(i),sqrt(mas2/(st2(i))))
4067 ENDIF
4068 mas2 = two * ms(ix3(i))
4069 IF(mas2>zero.AND.st3(i)>zero.AND.
4070 . irb(kinet(ix3(i)))==0.AND.irb2(kinet(ix3(i)))==0)THEN
4071 dtmi(i) = min(dtmi(i),sqrt(mas2/(st3(i))))
4072 ENDIF
4073 mas2 = two * ms(ix4(i))
4074 IF(mas2>zero.AND.st4(i)>zero.AND.
4075 . irb(kinet(ix4(i)))==0.AND.irb2(kinet(ix4(i)))==0)THEN
4076 dtmi(i) = min(dtmi(i),sqrt(mas2/(st4(i))))
4077 ENDIF
4078 dtmi0 = min(dtmi0,dtmi(i))
4079 ENDDO
4080 ENDIF
4081 IF(dtmi0<=dtm)THEN
4082 DO i=1,jlt-jlt_tied
4083 IF(dtmi(i)<=dtm.AND.cand_f(1,index(i))==zero)THEN
4084 jg = nsvg(i)
4085 IF(jg>0)THEN
4086 ni = itab(jg)
4087 ELSE
4088 ni = itafi(nin)%P(-jg)
4089 ENDIF
4090#include "lockon.inc"
4091 WRITE(iout,'(A,E12.4,A,I10,A,E12.4,A)')
4092 . ' **WARNING MINIMUM TIME STEP ',dtmi(i),
4093 . ' IN INTERFACE ',noint,'(DTMIN=',dtm,')'
4094 WRITE(iout,'(A,I10,A,I10)')' FREEZE SECONDARY NODE VS MAIN SEGMENT',
4095 . ni,' IN INTERFACE ',noint
4096 WRITE(iout,'(A,4I10)')' MAIN NODES : ',
4097 . itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i))
4098#include "lockoff.inc"
4099 cand_f(1,index(i))= fns(i)*sign(one,side(i)) ! negative vs NX NY NZ
4100 cand_f(2,index(i))= fxt(i)*t1x(i)+fyt(i)*t1y(i)+fzt(i)*t1z(i)
4101 cand_f(3,index(i))= fxt(i)*t2x(i)+fyt(i)*t2y(i)+fzt(i)*t2z(i)
4102 cand_f(4,index(i))= h1(i)
4103 cand_f(5,index(i))= h2(i)
4104 cand_f(6,index(i))= h3(i)
4105 cand_f(7,index(i))= pene(i)
4106 cand_f(8,index(i))= side(i)
4107 ENDIF
4108 ENDDO
4109 ENDIF
4110
4111 ENDIF ! IF(ITIED == 0)THEN, ELSE
4112C=======================================================================
4113 DO i=1,jlt
4114 fx1(i)=fxi(i)*h1(i)
4115 fy1(i)=fyi(i)*h1(i)
4116 fz1(i)=fzi(i)*h1(i)
4117C
4118 fx2(i)=fxi(i)*h2(i)
4119 fy2(i)=fyi(i)*h2(i)
4120 fz2(i)=fzi(i)*h2(i)
4121C
4122 fx3(i)=fxi(i)*h3(i)
4123 fy3(i)=fyi(i)*h3(i)
4124 fz3(i)=fzi(i)*h3(i)
4125C
4126 fx4(i)=fxi(i)*h4(i)
4127 fy4(i)=fyi(i)*h4(i)
4128 fz4(i)=fzi(i)*h4(i)
4129 ENDDO
4130C
4131C SPMD: Identification of interf nodes.useful to send
4132 IF (nspmd>1) THEN
4133Ctmp+1 mic only to avoid compiler bug
4134#include "mic_lockon.inc"
4135 DO i = 1,jlt
4136 nn = nsvg(i)
4137 IF(nn<0)THEN
4138C temporary tag of nsvfi a -
4139 nsvfi(nin)%P(-nn) = -abs(nsvfi(nin)%P(-nn))
4140 ENDIF
4141 ENDDO
4142ctmp+1 mic only
4143#include "mic_lockoff.inc"
4144 ENDIF
4145C
4146C-----------------------------------------------------
4147 IF((ibag/=0).OR.(iadm/=0).OR.(idamp_rdof/=0)) THEN
4148 DO i=1,jlt
4149C IF(PENE(I)/=ZERO)THEN
4150C test modified for consistency with SPMD communication (spmd_i7tools)
4151 IF(fxi(i)/=zero.OR.fyi(i)/=zero.OR.fzi(i)/=zero)THEN
4152 jg = nsvg(i)
4153 IF(jg>0) THEN
4154C In SPMD: Treatment to be redone after reception node Remote if JG <0
4155 icontact(jg)=1
4156 ENDIF
4157 icontact(ix1(i))=1
4158 icontact(ix2(i))=1
4159 icontact(ix3(i))=1
4160 icontact(ix4(i))=1
4161 ENDIF
4162 ENDDO
4163 ENDIF
4164C
4165 IF(iadm/=0)THEN
4166 DO i=1,jlt
4167 jg = nsvg(i)
4168#include "lockon.inc"
4169 IF(jg>0) THEN
4170C In SPMD: Treatment to be redone after reception node Remote if JG <0
4171 rcontact(jg)=min(rcontact(jg),rcurvi(i))
4172 END IF
4173 rcontact(ix1(i))=min(rcontact(ix1(i)),rcurvi(i))
4174 rcontact(ix2(i))=min(rcontact(ix2(i)),rcurvi(i))
4175 rcontact(ix3(i))=min(rcontact(ix3(i)),rcurvi(i))
4176 rcontact(ix4(i))=min(rcontact(ix4(i)),rcurvi(i))
4177#include "lockoff.inc"
4178 END DO
4179 END IF
4180 IF(iadm>=2)THEN
4181 DO i=1,jlt
4182 jg = nsvg(i)
4183#include "lockon.inc"
4184 IF(jg>0) THEN
4185C In SPMD: Treatment to be redone after reception node Remote if JG <0
4186 pcontact(jg)=max(pcontact(jg),pene(i)/(padm*gapv(i)))
4187 acontact(jg)=min(acontact(jg),anglmi(i))
4188 END IF
4189#include "lockoff.inc"
4190 END DO
4191 END IF
4192C
4193 IF(ibcc==0) RETURN
4194C
4195 DO 400 i=1,jlt
4196 IF(pene(i)==zero)GOTO 400
4197 ibcm = ibcc / 8
4198 ibcs = ibcc - 8 * ibcm
4199 IF(ibcs>0) THEN
4200 ig=nsvg(i)
4201 IF(ig>0) THEN
4202C In SPMD: Treatment to be redone after reception node Remote if JG <0
4203 CALL ibcoff(ibcs,icodt(ig))
4204 ENDIF
4205 ENDIF
4206 IF(ibcm>0) THEN
4207 ig=ix1(i)
4208 CALL ibcoff(ibcm,icodt(ig))
4209 ig=ix2(i)
4210 CALL ibcoff(ibcm,icodt(ig))
4211 ig=ix3(i)
4212 CALL ibcoff(ibcm,icodt(ig))
4213 ig=ix4(i)
4214 CALL ibcoff(ibcm,icodt(ig))
4215 ENDIF
4216 400 CONTINUE
4217C
4218 RETURN
4219 END
4220
the
end diagonal values have been computed in the(sparse) matrix id.SOL
alpha
#define alpha
Definition eval.h:35
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
i7curv
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)
Definition i7curv.F:443
i7for3
subroutine i7for3(output, jlt, a, v, ibcc, icodt, fsav, gap, fric, ms, visc, viscf, noint, stfn, itab, cn_loc, stiglo, stifn, stif, fskyi, isky, nx1, nx2, nx3, nx4, ny1, ny2, ny3, ny4, nz1, nz2, nz3, nz4, lb1, lb2, lb3, lb4, lc1, lc2, lc3, lc4, p1, p2, p3, p4, fcont, ix1, ix2, ix3, ix4, nsvg, ivis2, neltst, ityptst, dt2t, gapv, inacti, cand_p, index, kinet, newfront, isecin, nstrf, x, irect, ce_loc, mfrot, ifq, cand_fx, cand_fy, cand_fz, alpha0, ifpen, ibag, icontact, viscn, vxi, vyi, vzi, msi, kini, nin, nisub, lisub, addsubs, addsubm, lisubs, lisubm, fsavsub, cand_n, ilagm, icurv, nod_normal, fncont, ftcont, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi, iadm, rcurvi, rcontact, acontact, pcontact, anglmi, padm, intth, temp, tempi, iform, npc, tf, cmaj, dtmini, drad, fheats, fheatm, efrict, qfric, fni, ifric, jtask, h1, h2, h3, h4, ks, kt, k1, k2, k3, k4, c1, c2, c3, c4, cs, c, cf, tint, xfric, fxi, fyi, fzi, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, isensint, fsavparit, nft, sym_flag_type19, h3d_data, fricc, viscffric, fric_coefs, itied, jlt_tied, cand_f, iorthfric, fric_coefs2, fricc2, viscffric2, nforth, nfisot, indexorth, indexisot, dir1, dir2, tagncont, kloadpinter, loadpinter, loadp_hyd_inter, typsub, inflg_subs, inflg_subm, ninloadp, dgaploadint, s_loadpinter, dgaploadp, interefric)
Definition i7for3.F:78
ibcoff
subroutine ibcoff(ibc, icodt)
Definition ibcoff.F:44
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
h3d_mod
Definition h3d_mod.F:311
tri7box::inflg_subsfi
type(int_pointer), dimension(:), allocatable inflg_subsfi
Definition tri7box.F:505
tri7box::fnconti
type(real_pointer2), dimension(:), allocatable fnconti
Definition tri7box.F:510
tri7box::efricgfi
type(real_pointer), dimension(:), allocatable efricgfi
Definition tri7box.F:511
tri7box::stifi
type(real_pointer), dimension(:), allocatable stifi
Definition tri7box.F:449
tri7box::lisubsfi
type(int_pointer), dimension(:), allocatable lisubsfi
Definition tri7box.F:501
tri7box::nsvfi
type(int_pointer), dimension(:), allocatable nsvfi
Definition tri7box.F:431
tri7box::efricfi
type(real_pointer), dimension(:), allocatable efricfi
Definition tri7box.F:511
tri7box::addsubsfi
type(int_pointer), dimension(:), allocatable addsubsfi
Definition tri7box.F:509
tri7box::ftconti
type(real_pointer2), dimension(:), allocatable ftconti
Definition tri7box.F:510
tri7box::itafi
type(int_pointer), dimension(:), allocatable itafi
Definition tri7box.F:440
main
int main(int argc, char *argv[])
Definition CMakeCCompilerId.c:841
h3d_mod::h3d_database
Definition h3d_mod.F:438