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