OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25for3.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "com06_c.inc"
#include "com08_c.inc"
#include "scr05_c.inc"
#include "scr07_c.inc"
#include "scr11_c.inc"
#include "scr14_c.inc"
#include "scr16_c.inc"
#include "scr18_c.inc"
#include "sms_c.inc"
#include "param_c.inc"
#include "impl1_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i25for3 (jlt, a, v, ibcc, icodt, fsav, ms, visc, viscf, noint, stfn, itab, cn_loc, stiglo, stifn, stif, inacti, index, n1, n2, n3, h1, h2, h3, h4, fcont, pene, nrtm, ix1, ix2, ix3, ix4, nsvg, ivis2, neltst, ityptst, dt2t, kinet, newfront, isecin, nstrf, secfcum, x, irect, ce_loc, mfrot, ifq, secnd_fr, alpha0, ibag, icontact, irtlm, viscn, vxi, vyi, vzi, msi, kini, nin, nisub, lisub, addsubs, addsubm, lisubs, lisubm, inflg_subs, inflg_subm, fsavsub, ilagm, icurv, fncont, ftcont, nsn, xx, yy, zz, xi, yi, zi, anglmi, padm, iadm, rcurvi, rcontact, acontact, pcontact, mskyi_sms, iskyi_sms, nsms, cand_n_n, pene_old, stif_old, mbinflg, ilev, igsti, kmin, intply, nm1, nm2, nm3, msegtyp, jtask, isensint, fsavparit, h3d_data, fricc, viscffric, fric_coefs, gapv, viscfluid, sigmaxadh, viscadhfact, if_adh, areas, base_adh, iorthfric, fric_coefs2, fricc2, viscffric2, nforth, nfisot, indexorth, indexisot, dir1, dir2, apinch, stifpinch, fni, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, fxi, fyi, fzi, intth, drad, fheats, fheatm, qfric, efrict, tagncont, kloadpinter, loadpinter, loadp_hyd_inter, typsub, ncfit, ninloadp, dgaploadint, s_loadpinter, dist, dgaploadpmax, interefric, intcarea, parameters)
subroutine i25sms2 (jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, nin, noint, mskyi_sms, iskyi_sms, nsms, kt, c, cf, dtmini, dti)

Function/Subroutine Documentation

◆ i25for3()

subroutine i25for3 ( integer jlt,
a,
v,
integer ibcc,
integer, dimension(*) icodt,
fsav,
ms,
visc,
viscf,
integer noint,
stfn,
integer, dimension(*) itab,
integer, dimension(mvsiz) cn_loc,
stiglo,
stifn,
stif,
integer inacti,
integer, dimension(mvsiz) index,
n1,
n2,
n3,
h1,
h2,
h3,
h4,
fcont,
pene,
integer nrtm,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) nsvg,
integer ivis2,
integer neltst,
integer ityptst,
dt2t,
integer, dimension(*) kinet,
integer newfront,
integer isecin,
integer, dimension(*) nstrf,
secfcum,
x,
integer, dimension(4,*) irect,
integer, dimension(mvsiz) ce_loc,
integer mfrot,
integer ifq,
secnd_fr,
alpha0,
integer ibag,
integer, dimension(*) icontact,
integer, dimension(4,nsn) irtlm,
viscn,
vxi,
vyi,
vzi,
msi,
integer, dimension(*) kini,
integer nin,
integer nisub,
integer, dimension(*) lisub,
integer, dimension(*) addsubs,
integer, dimension(*) addsubm,
integer, dimension(*) lisubs,
integer, dimension(*) lisubm,
integer, dimension(*) inflg_subs,
integer, dimension(*) inflg_subm,
fsavsub,
integer ilagm,
integer, dimension(3) icurv,
fncont,
ftcont,
integer nsn,
xx,
yy,
zz,
xi,
yi,
zi,
anglmi,
padm,
integer iadm,
rcurvi,
rcontact,
acontact,
pcontact,
mskyi_sms,
integer, dimension(*) iskyi_sms,
integer, dimension(*) nsms,
integer, dimension(*) cand_n_n,
pene_old,
stif_old,
integer, dimension(*) mbinflg,
integer ilev,
integer igsti,
kmin,
integer intply,
nm1,
nm2,
nm3,
integer, dimension(*) msegtyp,
integer jtask,
integer, dimension(*) isensint,
fsavparit,
type(h3d_database) h3d_data,
fricc,
viscffric,
fric_coefs,
gapv,
viscfluid,
sigmaxadh,
viscadhfact,
integer, dimension(mvsiz) if_adh,
areas,
base_adh,
integer iorthfric,
fric_coefs2,
fricc2,
viscffric2,
integer nforth,
integer nfisot,
integer, dimension(mvsiz) indexorth,
integer, dimension(mvsiz) indexisot,
dir1,
dir2,
apinch,
stifpinch,
fni,
fx1,
fy1,
fz1,
fx2,
fy2,
fz2,
fx3,
fy3,
fz3,
fx4,
fy4,
fz4,
fxi,
fyi,
fzi,
integer intth,
drad,
fheats,
fheatm,
qfric,
efrict,
integer, dimension(nloadp_hyd_inter,numnod) tagncont,
integer, dimension(ninter+1), intent(in) kloadpinter,
integer, dimension(s_loadpinter), intent(in) loadpinter,
integer, dimension(nloadp_hyd), intent(in) loadp_hyd_inter,
integer, dimension(*) typsub,
integer ncfit,
integer, intent(in) ninloadp,
dimension(s_loadpinter), intent(in) dgaploadint,
integer, intent(in) s_loadpinter,
dimension(mvsiz), intent(inout) dist,
intent(in) dgaploadpmax,
integer, intent(in) interefric,
integer, intent(in) intcarea,
type (parameters_), intent(in) parameters )

Definition at line 40 of file i25for3.F.

74C-----------------------------------------------
75C M o d u l e s
76C-----------------------------------------------
77 USE tri7box
78 USE plyxfem_mod
79 USE h3d_mod
80 USE pinchtype_mod
81 USE anim_mod
82 USE outputs_mod
83 USE parameters_mod
84C-----------------------------------------------
85C I m p l i c i t T y p e s
86C-----------------------------------------------
87#include "implicit_f.inc"
88#include "comlock.inc"
89C-----------------------------------------------
90C G l o b a l P a r a m e t e r s
91C-----------------------------------------------
92#include "mvsiz_p.inc"
93C-----------------------------------------------
94C C o m m o n B l o c k s
95C-----------------------------------------------
96#include "com01_c.inc"
97#include "com04_c.inc"
98#include "com06_c.inc"
99#include "com08_c.inc"
100#include "scr05_c.inc"
101#include "scr07_c.inc"
102#include "scr11_c.inc"
103#include "scr14_c.inc"
104#include "scr16_c.inc"
105#include "scr18_c.inc"
106#include "sms_c.inc"
107#include "param_c.inc"
108#include "impl1_c.inc"
109C-----------------------------------------------
110C D u m m y A r g u m e n t s
111C-----------------------------------------------
112 INTEGER NELTST,ITYPTST,JLT,IBCC,IBCM,IBCS,IVIS2,INACTI,IBAG,NIN,
113 . NSN,IORTHFRIC,NFORTH ,NFISOT,
114 . ICODT(*), ITAB(*), KINET(*),IRTLM(4,NSN),
115 . MFROT, IFQ, NOINT,NEWFRONT,ISECIN, NSTRF(*),
116 . IRECT(4,*),ICONTACT(*),
117 . KINI(*),MBINFLG(*),ILEV,JTASK,
118 . ISET, IADM,INTTH,IFORM,NRTM,
119 . CAND_N_N(*),INTPLY,MSEGTYP(*),TAGNCONT(NLOADP_HYD_INTER,NUMNOD)
120 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
121 . CN_LOC(MVSIZ),CE_LOC(MVSIZ),INDEX(MVSIZ),NSVG(MVSIZ),
122 . NISUB, LISUB(*), ADDSUBS(*), ADDSUBM(*), LISUBS(*), LISUBM(*),
123 . INFLG_SUBS(*), INFLG_SUBM(*), ILAGM, ICURV(3),
124 . ISKYI_SMS(*), NSMS(*),IGSTI,
125 . ISENSINT(*), IF_ADH(MVSIZ),
126 . INDEXORTH(MVSIZ) ,INDEXISOT(MVSIZ),TYPSUB(*),NCFIT
127 INTEGER , INTENT(IN) :: NINLOADP,S_LOADPINTER
128 INTEGER , INTENT(IN) :: KLOADPINTER(NINTER+1),LOADPINTER(S_LOADPINTER),
129 . LOADP_HYD_INTER(NLOADP_HYD)
130 INTEGER , INTENT(IN) :: INTEREFRIC, INTCAREA
131 my_real
132 . stiglo, x(3,*),
133 . a(3,*), ms(*), v(3,*), fsav(*),fcont(3,*),
134 . secnd_fr(6,*),alpha0,
135 . visc,viscf,vis,dt2t,stfn(*),stifn(*),
136 . fsavsub(nthvki,*), fncont(3,*),ftcont(3,*),
137 . mskyi_sms(*),kmin, viscfluid, sigmaxadh, viscadhfact,
138 . apinch(3,*),stifpinch(*)
139 my_real
140 . fxi(mvsiz), fyi(mvsiz), fzi(mvsiz), fni(mvsiz),
141 . fx1(mvsiz), fx2(mvsiz), fx3(mvsiz), fx4(mvsiz),
142 . fy1(mvsiz), fy2(mvsiz), fy3(mvsiz), fy4(mvsiz),
143 . fz1(mvsiz), fz2(mvsiz), fz3(mvsiz), fz4(mvsiz),
144 . la(mvsiz), lb(mvsiz) , lc(mvsiz) ,stif(mvsiz),
145 . pene(mvsiz),
146 . secfcum(7,numnod,nsect),
147 . tmp(mvsiz),
148 . stifsav(mvsiz), viscn(*),
149 . vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),msi(mvsiz),
150 . xx(mvsiz,5),yy(mvsiz,5),zz(mvsiz,5),
151 . xi(mvsiz),yi(mvsiz),zi(mvsiz),
152 . n1(mvsiz), n2(mvsiz), n3(mvsiz),
153 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
154 . nm1(mvsiz), nm2(mvsiz), nm3(mvsiz),
155 . fsavparit(nisub+1,11,*),fric_coefs(mvsiz,10), fricc(mvsiz),
156 . viscffric(mvsiz), gapv(mvsiz), areas(mvsiz), base_adh(mvsiz),
157 . fric_coefs2(mvsiz,10),fricc2(mvsiz),viscffric2(mvsiz), dir1(mvsiz,3),
158 . dir2(mvsiz,3), efrict(mvsiz) ,
159 . drad, fheats, fheatm, qfric
160 my_real
161 . rcurvi(*), rcontact(*), acontact(*),
162 . pcontact(*),padm, anglmi(*),rstif,
163 . pene_old(5,nsn),stif_old(2,nsn)
164 my_real , INTENT(IN) :: dgaploadpmax
165 my_real , INTENT(IN) :: dgaploadint(s_loadpinter)
166 my_real , INTENT(INOUT) :: dist(mvsiz)
167 TYPE(H3D_DATABASE) :: H3D_DATA
168 TYPE (PARAMETERS_) ,INTENT(IN):: PARAMETERS
169C-----------------------------------------------
170C L o c a l V a r i a b l e s
171C-----------------------------------------------
172 INTEGER I, J1, IG, J, JG , K0,NBINTER,K1S,K,IL,IE, NN, NI,
173 . NA1,NA2,N,ITAG,NN1,NN2,NN3,NN4,
174 . IGRN,ISS1,ISS2,IMS1,IMS2,PP,PPL,ITYPSUB
175 my_real
176 . fxr(mvsiz), fyr(mvsiz), fzr(mvsiz)
177 my_real
178 . fxt(mvsiz),fyt(mvsiz),fzt(mvsiz),
179 . vis2(mvsiz), dtmi(mvsiz), xmu(mvsiz),stif0(mvsiz), ff(mvsiz),
180 . vx(mvsiz), vy(mvsiz), vz(mvsiz), vn(mvsiz),
181 . xp(mvsiz), yp(mvsiz), zp(mvsiz), efric_l(mvsiz),
182 . vnx, vny, vnz, aa, crit,s2,
183 . v2, fm2, dt1inv, visca, fac,alphi,alpha,beta,
184 . fx, fy, fz, f2, mas2, m2sk, dtmi0,dti,ft,fn,fmax,ftn,
185 . facm1, econtt, econvt, h0, econtdt,
186 . d1,d2,d3,d4,a1,a2,a3,a4,efric_ls,efric_lm,
187 . fsav1, fsav2, fsav3, fsav4, fsav5, fsav6, fsav7, fsav8,
188 . fsav9, fsav10, fsav11, fsav12, fsav13, fsav14, fsav15,
189 . fsav22, fsav23, fsav24, fsav29, ffo,
190 . e10, h0d, s2d, sum,
191 . la1d,la2d,la3d,la4d,t1,t1d,t2,t2d,ffd,visd,facd,d1d,
192 . p1s(mvsiz),p2s(mvsiz),p3s(mvsiz),p4s(mvsiz),
193 . d2d,d3d,d4d,vnxd,vnyd,vnzd,v2d,fm2d,f2d,aad,fxd,fyd,fzd,
194 . a1d,a2d,a3d,a4d,vv,ax1,ax2,ay1,ay2,az1,az2,ax,ay,az,
195 . area,p,vv1,vv2,v21,dmu, dti2,h00 ,a0x,a0y,a0z,rx,ry,rz,
196 . anx,any,anz,aan,aax,aay,aaz ,rr,rs,aaa,bbb,stfr,visr ,tm,ts,
197 . vn1(3),vn2(3),vn3(3),vn4(4),
198 . csa ,sna ,alphaf ,ftt1 ,ftt2 ,qfrict,
199 . an(nforth) ,bn(nforth) ,nep1 ,nep2 ,nep ,c11 ,c22 ,ans ,ep,signc
200 my_real
201 . prec
202 my_real
203 . st1(mvsiz),st2(mvsiz),st3(mvsiz),st4(mvsiz),stv(mvsiz),
204 . kt(mvsiz),c(mvsiz),cf(mvsiz),dpene(mvsiz),
205 . ks(mvsiz),k1(mvsiz),k2(mvsiz),k3(mvsiz),k4(mvsiz),
206 . cs(mvsiz),c1(mvsiz),c2(mvsiz),c3(mvsiz),c4(mvsiz),
207 . cx,cy,cfi,aux, fx6(6,mvsiz), fy6(6,mvsiz), fz6(6,mvsiz),
208 . dtmini,stif0_imp(mvsiz),penn,tiny,xmu2(mvsiz)
209 INTEGER JSUB,KSUB,JJ,KK,IN,NSUB
210 INTEGER ISEGTYP(MVSIZ),IXSS(4),NOD1,NOD2,NS
211 my_real
212 . fsavsub1(25,nisub),impx,impy,impz,vnm(mvsiz),sfac,stmin,
213 . fxs(4),fys(4),fzs(4)
214 my_real
215 . stifadh, factor , gapp, dgapload
216C
217 INTEGER BITGET
218 EXTERNAL bitget
219C-----------------------------------------------
220 dtmini=zero
221 IF (iresp==1) THEN
222 prec = fiveem4
223 tiny = em15
224 ELSE
225 prec = em10
226 tiny = em30
227 ENDIF
228 IF(dt1>zero)THEN
229 dt1inv = one/dt1
230 ELSE
231 dt1inv =zero
232 ENDIF
233 econtt = zero
234 econvt = zero
235 econtdt = zero
236 qfrict = zero
237 DO i=1,jlt
238 stif0(i) = stif(i)
239 ENDDO
240
241 efrict(1:jlt) = zero
242 efric_l(1:jlt) = zero
243C-------------------------------------------
244C PENE Offset removed to i24dst3.F
245C-------------------------------------------
246C
247 DO i=1,jlt
248 IF(pene(i) == zero.AND.intth == 0.AND.(ninloadp==0.OR.dgaploadpmax==zero))THEN
249 h1(i) = zero
250 h2(i) = zero
251 h3(i) = zero
252 h4(i) = zero
253 fni(i)= zero
254 ELSEIF(pene(i) == zero)THEN
255 fni(i)= zero
256C
257 fxi(i)=zero
258 fyi(i)=zero
259 fzi(i)=zero
260C
261 fx1(i)=zero
262 fy1(i)=zero
263 fz1(i)=zero
264C
265 fx2(i)=zero
266 fy2(i)=zero
267 fz2(i)=zero
268C
269 fx3(i)=zero
270 fy3(i)=zero
271 fz3(i)=zero
272C
273 fx4(i)=zero
274 fy4(i)=zero
275 fz4(i)=zero
276 ELSE
277 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
278 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
279 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
280 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
281 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
282 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
283 vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
284
285 ENDIF
286 ENDDO
287C-------------------------------------------
288C stifness reduction to avoid positive
289C force jump and energy generation
290C-------------------------------------------
291 DO i=1,jlt
292C
293 IF(pene(i) == zero) cycle
294C
295 dpene(i) = max(zero,-vn(i)*dt1)
296 jg = nsvg(i)
297 n = cand_n_n(i)
298 IF(jg > 0)THEN
299 IF(tt /= zero)THEN
300 IF(pene(i) > pene_old(2,n) .AND.
301 . dpene(i) < pene(i)-pene_old(2,n))THEN
302 stif(i) = (stif_old(2,n)*pene_old(2,n)+ stif(i)*dpene(i))
303 . /pene(i)
304 ELSE
305 stif(i) = stif_old(2,n)
306 ENDIF
307 pene_old(1,n) = pene(i)
308 stif_old(1,n) = stif(i)
309 ELSE
310 pene_old(1,n) = max(pene_old(1,n),pene(i))
311 stif_old(1,n) = max(stif_old(1,n),stif(i))
312 ENDIF
313c if(itab(nsvg(i))==7444)
314c . print *,'i25for3 stif',itab(jg),pene(i),pene_old(2,n),pene(i)-pene_old(2,n),dpene(i),stif_old(2,n),stif(i)
315 ELSE
316 jg = -jg
317 IF(tt /= zero)THEN
318 IF(pene(i) > pene_oldfi(nin)%P(2,jg) .and.
319 . dpene(i) < pene(i)-pene_oldfi(nin)%P(2,jg))THEN
320 stif(i) = (stif_oldfi(nin)%P(2,jg)*pene_oldfi(nin)%P(2,jg)
321 . + stif(i)*dpene(i)) / pene(i)
322 ELSE
323 stif(i) = stif_oldfi(nin)%P(2,jg)
324 ENDIF
325 pene_oldfi(nin)%P(1,jg) = pene(i)
326 stif_oldfi(nin)%P(1,jg) = stif(i)
327 ELSE
328 pene_oldfi(nin)%P(1,jg) =
329 * max(pene_oldfi(nin)%P(1,jg),pene(i))
330 stif_oldfi(nin)%P(1,jg) =
331 * max(stif_oldfi(nin)%P(1,jg),stif(i))
332 ENDIF
333 ENDIF
334c if(nsvg(i)>0.and.3827132.or.
335c . itab(ix1(i))==3827132.or.itab(ix2(i))==3827132.or.itab(ix3(i))==3827132.or.itab(ix4(i))==3827132)
336c . print *,'i25for3 stif',STIF_OLD(2,N),pene_old(2,n),stif_old(1,n),pene_old(1,n),pene(i),pene_old(5,n),
337c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i))
338 ENDDO
339C-------------------------------------------
340 IF (ncfit >0) THEN
341 DO i=1,jlt
342 IF(pene(i) == zero)cycle
343 stif(i) = stif0(i)
344 IF(stiglo<=zero) stif(i) = half*stif(i)
345 fni(i)= -stif(i) * pene(i)
346 ENDDO
347 ELSEIF(ivis2/=-1)THEN ! NO ADHESION CASE
348 DO i=1,jlt
349 IF(pene(i) == zero)cycle
350 IF(stiglo<=zero)THEN
351 stif(i) = half*stif(i)
352 ELSEIF(stif(i)/=zero)THEN
353 stif(i) = stiglo
354 ENDIF
355 fni(i)= -stif(i) * pene(i)
356 ENDDO
357 ELSE ! ADHESION CASE
358C BASE_ADH = base where the adhesion spring is connected
359C IF_ADH = 1 if adhesion spring exists else 0
360C PENE = measured from the outer boundary of adhesion zone (i.e. 2*(real_gap))
361 DO i=1,jlt
362 IF(pene(i) == zero) cycle
363 IF(stiglo<=zero)THEN
364 stif(i) = half*stif(i)
365 ELSEIF(stif(i)/=zero)THEN
366 stif(i) = stiglo
367 ENDIF
368 stifadh = sigmaxadh*areas(i)/max(em30,base_adh(i))
369 IF(pene(i) < base_adh(i))THEN ! INSIDE ADHESION ZONE
370 fni(i) = stifadh*if_adh(i)*(base_adh(i)-pene(i))
371 ELSE ! INSIDE PENETRATION ZONE
372 fni(i) = -stif(i)*(pene(i)-base_adh(i))
373 ENDIF
374 ENDDO
375 ENDIF
376C---------------------------------
377C Contact Energy (elastic)
378C---------------------------------
379 DO i=1,jlt
380 IF(pene(i) == zero)cycle
381 econtt = econtt+half*stif(i)*pene(i)**2
382 END DO
383C---------------------------------
384C DAMPING + FRIC
385C---------------------------------
386 ff(1:jlt)=zero
387C
388 IF(visc/=zero)THEN
389 IF(ivis2==6)THEN
390C---------------------------------
391C VISC QUAD TYPE V227
392C---------------------------------
393 DO i=1,jlt
394 vis2(i) = two * stif(i) * msi(i)
395 ENDDO
396C---------------------------------
397 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
398 DO i=1,jlt
399 IF(pene(i) == zero)cycle
400 fac = stif(i) / max(em30,stif(i))
401 vis = sqrt(vis2(i))
402 ff(i) = fac * visc * vis
403 stif(i) = stif0(i) + ff(i) * dt1inv
404 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
405 ff(i) = ff(i) * vn(i)
406 fni(i) = fni(i) + ff(i)
407 ENDDO
408 ELSE
409 DO i=1,jlt
410 IF(pene(i) == zero)cycle
411 fac = stif(i) / max(em30,stif(i))
412 vis = sqrt(vis2(i))
413 c(i)= fac * visc * vis
414 kt(i)= stif0(i)
415 stif(i) = stif0(i) + c(i) * dt1inv
416 ff(i) = c(i) * vn(i)
417 fni(i) = fni(i) + ff(i)
418 cf(i) = fac*sqrt(viscffric(i))*vis
419 stif(i) = max(stif(i) ,cf(i)*dt1inv)
420 ENDDO
421 ENDIF
422 ELSEIF(ivis2==1)THEN
423C---------------------------------
424C VISC QUAD TYPE V227
425C---------------------------------
426 DO i=1,jlt
427 IF(pene(i) == zero)cycle
428 mas2 = ms(ix1(i))*h1(i)
429 . + ms(ix2(i))*h2(i)
430 . + ms(ix3(i))*h3(i)
431 . + ms(ix4(i))*h4(i)
432 vis2(i) = two* stif(i) * msi(i) * mas2 /
433 . max(em30,msi(i)+mas2)
434 ENDDO
435C---------------------------------
436 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
437 DO i=1,jlt
438 IF(pene(i) == zero)cycle
439 fac = stif(i) / max(em30,stif(i))
440 vis = sqrt(vis2(i))
441 ff(i) = fac * visc * vis
442 stif(i) = stif0(i) + ff(i) * dt1inv
443 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
444 ff(i) = ff(i) * vn(i)
445 fni(i) = fni(i) + ff(i)
446 ENDDO
447 ELSE
448 DO i=1,jlt
449 IF(pene(i) == zero)cycle
450 fac = stif(i) / max(em30,stif(i))
451 vis = sqrt(vis2(i))
452 c(i)= fac * visc * vis
453 kt(i)= stif0(i)
454 stif(i) = stif(i) + c(i) * dt1inv
455 ff(i) = c(i) * vn(i)
456 fni(i) = fni(i) + ff(i)
457 cf(i) = fac*sqrt(viscffric(i))*vis
458 stif(i) = max(stif(i) ,cf(i)*dt1inv)
459 ENDDO
460 ENDIF
461 ELSEIF(ivis2==2)THEN
462C---------------------------------
463C VISC QUAD TYPE
464C---------------------------------
465 DO i=1,jlt
466 IF(pene(i) == zero)cycle
467 vis2(i) = two* stif(i) * msi(i)
468 fac = stif(i) / max(em30,stif(i))
469 vis = sqrt(vis2(i))
470 ff(i) = fac * visc * vis
471 stif(i) = stif0(i) + two * ff(i) * dt1inv
472 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
473 ff(i) = ff(i) * vn(i)
474 fni(i) = fni(i) + ff(i)
475 ENDDO
476 ELSEIF(ivis2==3)THEN
477C---------------------------------
478C VISC QUAD = 0
479C---------------------------------
480 DO i=1,jlt
481 IF(pene(i) == zero)cycle
482 vis2(i) = two * stif(i) * msi(i)
483 fac = stif(i) / max(em30,stif(i))
484 vis = sqrt(vis2(i))
485 ff(i) = fac * visc * vis
486 stif(i) = stif0(i) + two* ff(i) * dt1inv
487 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
488 ff(i) = ff(i) * vn(i)
489 fni(i) = fni(i) + ff(i)
490 ENDDO
491 ELSEIF(ivis2==4)THEN
492C---------------------------------
493C VISC = 0
494C---------------------------------
495 DO i=1,jlt
496 IF(pene(i) == zero)cycle
497 fac = stif(i) / max(em30,stif(i))
498 vis2(i) = two* stif(i) * msi(i)
499 vis = sqrt(vis2(i))
500 stif(i) = max(stif0(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
501 ENDDO
502 ELSEIF(ivis2==5)THEN
503C---------------------------------
504C VISC = 2M/dt => pour visc < 1 , stable : dt < 2M/visc = dt
505C M=M1*M2/M1+M2 pour visc = 1 , choc elastique
506C pour visc = 0.5 , choc mou
507C---------------------------------
508 DO i=1,jlt
509 IF(pene(i) == zero)cycle
510 fac = zero ! missing proper definition of FAC for IVIS2=5 (unused)
511 mas2 = ms(ix1(i))*h1(i)
512 . + ms(ix2(i))*h2(i)
513 . + ms(ix3(i))*h3(i)
514 . + ms(ix4(i))*h4(i)
515 vis2(i) = two* stif(i) * msi(i)
516 vis = 2. * visc * dt1inv * msi(i) * mas2 /
517 . max(em30,msi(i)+mas2)
518 stif(i) = max(stif0(i) ,fac*sqrt(viscffric(i)*vis2(i))*dt1inv)
519 ff(i) = min(zero,ff(i)-fni(i))
520 fni(i) = fni(i)+ff(i)
521 ENDDO
522 ELSEIF(ivis2==-1)THEN
523C---------------------------------
524C ADHESION MODEL, NORMAL DAMPING WITH VISC + ADHESION SPRING
525C---------------------------------
526 DO i=1,jlt
527 IF(pene(i) == zero)cycle
528 mas2 = ms(ix1(i))*h1(i)
529 . + ms(ix2(i))*h2(i)
530 . + ms(ix3(i))*h3(i)
531 . + ms(ix4(i))*h4(i)
532 vis2(i) = two* stif(i) * msi(i) * mas2 /
533 . max(em30,msi(i)+mas2)
534 ENDDO
535C---------------------------------
536 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
537 DO i=1,jlt
538 IF(pene(i) == zero)cycle
539 fac = stif(i) / max(em30,stif(i))
540 vis = sqrt(vis2(i))
541 stifadh = sigmaxadh*areas(i)/max(em30,base_adh(i))
542 ff(i) = fac * visc * vis
543 stif(i) = stif0(i) + ff(i) * dt1inv
544 stif(i) = max(stif(i), stifadh + ff(i) * dt1inv) ! NORMAL DIRECTION
545 stif(i) = max(stif(i),fac*viscadhfact*viscfluid*two/gapv(i)*areas(i)*dt1inv) ! TANGENTIAL DIRECTION
546 ff(i) = ff(i) * vn(i)
547 fni(i) = fni(i) + ff(i)
548 ENDDO
549 ELSE
550 DO i=1,jlt
551 IF(pene(i) == zero)cycle
552 fac = stif(i) / max(em30,stif(i))
553 vis = sqrt(vis2(i))
554 stifadh = sigmaxadh*areas(i)/max(em30,base_adh(i))
555 c(i) = fac * visc * vis
556 kt(i) = stif0(i)
557 stif(i) = stif(i) + c(i) * dt1inv
558 ff(i) = c(i) * vn(i)
559 fni(i) = fni(i) + ff(i)
560 cf(i) = fac*viscadhfact*viscfluid*two/gapv(i)*areas(i) ! tangential direction
561 stif(i) = max(stif(i), stifadh + c(i) * dt1inv) ! NORMAL DIRECTION
562 stif(i) = max(stif(i) ,cf(i)*dt1inv) ! TANGENTIAL DIRECTION
563 ENDDO
564 ENDIF
565 ELSE
566 ENDIF
567 ELSE
568 DO i=1,jlt
569
570 IF(viscffric(i)/=zero ) THEN
571
572 IF(ivis2==6)THEN
573C---------------------------------
574C VISC QUAD TYPE V227
575C---------------------------------
576 vis2(i) = two * stif(i) * msi(i)
577C---------------------------------
578C---------------------------------
579 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
580 IF(pene(i) == zero)cycle
581 fac = stif(i) / max(em30,stif(i))
582 vis = sqrt(vis2(i))
583 ff(i) = fac * visc * vis
584 stif(i) = stif0(i) + ff(i) * dt1inv
585 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
586 ff(i) = ff(i) * vn(i)
587 fni(i) = fni(i) + ff(i)
588 ELSE
589 IF(pene(i) == zero)cycle
590 fac = stif(i) / max(em30,stif(i))
591 vis = sqrt(vis2(i))
592 c(i)= fac * visc * vis
593 kt(i)= stif0(i)
594 stif(i) = stif0(i) + c(i) * dt1inv
595 ff(i) = c(i) * vn(i)
596 fni(i) = fni(i) + ff(i)
597 cf(i) = fac*sqrt(viscffric(i))*vis
598 stif(i) = max(stif(i) ,cf(i)*dt1inv)
599 ENDIF
600 ELSEIF(ivis2==1)THEN
601C---------------------------------
602C VISC QUAD TYPE V227
603C---------------------------------
604 IF(pene(i) == zero)cycle
605 mas2 = ms(ix1(i))*h1(i)
606 . + ms(ix2(i))*h2(i)
607 . + ms(ix3(i))*h3(i)
608 . + ms(ix4(i))*h4(i)
609 vis2(i) = two* stif(i) * msi(i) * mas2 /
610 . max(em30,msi(i)+mas2)
611C---------------------------------
612 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
613 IF(pene(i) == zero)cycle
614 fac = stif(i) / max(em30,stif(i))
615 vis = sqrt(vis2(i))
616 ff(i) = fac * visc * vis
617 stif(i) = stif0(i) + ff(i) * dt1inv
618 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
619 ff(i) = ff(i) * vn(i)
620 fni(i) = fni(i) + ff(i)
621 ELSE
622 IF(pene(i) == zero)cycle
623 fac = stif(i) / max(em30,stif(i))
624 vis = sqrt(vis2(i))
625 c(i)= fac * visc * vis
626 kt(i)= stif0(i)
627 stif(i) = stif(i) + c(i) * dt1inv
628 ff(i) = c(i) * vn(i)
629 fni(i) = fni(i) + ff(i)
630 cf(i) = fac*sqrt(viscffric(i))*vis
631 stif(i) = max(stif(i) ,cf(i)*dt1inv)
632 ENDIF
633 ELSEIF(ivis2==2)THEN
634C---------------------------------
635C VISC QUAD TYPE
636C---------------------------------
637 IF(pene(i) == zero)cycle
638 vis2(i) = two* stif(i) * msi(i)
639 fac = stif(i) / max(em30,stif(i))
640 vis = sqrt(vis2(i))
641 ff(i) = fac * visc * vis
642 stif(i) = stif0(i) + two * ff(i) * dt1inv
643 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
644 ff(i) = ff(i) * vn(i)
645 fni(i) = fni(i) + ff(i)
646 ELSEIF(ivis2==3)THEN
647C---------------------------------
648C VISC QUAD = 0
649C---------------------------------
650 IF(pene(i) == zero)cycle
651 vis2(i) = two * stif(i) * msi(i)
652 fac = stif(i) / max(em30,stif(i))
653 vis = sqrt(vis2(i))
654 ff(i) = fac * visc * vis
655 stif(i) = stif0(i) + two* ff(i) * dt1inv
656 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
657 ff(i) = ff(i) * vn(i)
658 fni(i) = fni(i) + ff(i)
659 ELSEIF(ivis2==4)THEN
660C---------------------------------
661C VISC = 0
662C---------------------------------
663 IF(pene(i) == zero)cycle
664 fac = stif(i) / max(em30,stif(i))
665 vis2(i) = two* stif(i) * msi(i)
666 vis = sqrt(vis2(i))
667 stif(i) = max(stif0(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
668 ELSEIF(ivis2==5)THEN
669C---------------------------------
670C VISC = 2M/dt => pour visc < 1 , stable : dt < 2M/visc = dt
671C M=M1*M2/M1+M2 pour visc = 1 , choc elastique
672C pour visc = 0.5 , choc mou
673C---------------------------------
674 IF(pene(i) == zero)cycle
675 mas2 = ms(ix1(i))*h1(i)
676 . + ms(ix2(i))*h2(i)
677 . + ms(ix3(i))*h3(i)
678 . + ms(ix4(i))*h4(i)
679 vis2(i) = two* stif(i) * msi(i)
680 vis = 2. * visc * dt1inv * msi(i) * mas2 /
681 . max(em30,msi(i)+mas2)
682 stif(i) = max(stif0(i) ,fac*sqrt(viscffric(i)*vis2(i))*dt1inv)
683 ff(i) = vis * vn(i)
684 ff(i) = min(zero,ff(i)-fni(i))
685 fni(i) = fni(i)+ff(i)
686 ELSE
687 ENDIF
688
689
690 ELSE
691cbb initialisation du tableau VIS2 pour eviter des problemes
692 vis2(i) = zero
693 ENDIF
694 ENDDO
695 ENDIF
696C---------------------------------
697C Energy absorbed (damping)
698C---------------------------------
699 DO i=1,jlt
700 IF(pene(i) == zero)cycle
701 jg = nsvg(i)
702 IF(jg > 0)THEN
703 n = cand_n_n(i)
704 econtdt = econtdt+pene_old(3,n)*vn(i)*dt1
705 pene_old(3,n) = half*ff(i)
706 econtdt = econtdt+pene_old(3,n)*vn(i)*dt1
707 ELSE
708 jg = -jg
709 econtdt = econtdt+pene_oldfi(nin)%P(3,jg)*vn(i)*dt1
710 pene_oldfi(nin)%P(3,jg) = half*ff(i)
711 econtdt = econtdt+pene_oldfi(nin)%P(3,jg)*vn(i)*dt1
712 END IF
713 END DO
714C---------------------------------
715C SAUVEGARDE DE L'IMPULSION NORMALE
716C---------------------------------
717 fsav1 = zero
718 fsav2 = zero
719 fsav3 = zero
720 fsav8 = zero
721 fsav9 = zero
722 fsav10= zero
723 fsav11= zero
724 IF (ilev==2) THEN
725 DO i=1,jlt
726 IF(pene(i) == zero)cycle
727 ie=ce_loc(i)
728 ims2 = bitget(mbinflg(ie),1)
729 fxi(i)=n1(i)*fni(i)
730 fyi(i)=n2(i)*fni(i)
731 fzi(i)=n3(i)*fni(i)
732 impx=fxi(i)*dt12
733 impy=fyi(i)*dt12
734 impz=fzi(i)*dt12
735 fsav8 =fsav8 +abs(impx)
736 fsav9 =fsav9 +abs(impy)
737 fsav10=fsav10+abs(impz)
738 IF (ims2 > 0 ) THEN
739 fsav1 =fsav1 -impx
740 fsav2 =fsav2 -impy
741 fsav3 =fsav3 -impz
742 fsav11=fsav11-fni(i)*dt12
743 ELSE
744 fsav1 =fsav1 +impx
745 fsav2 =fsav2 +impy
746 fsav3 =fsav3 +impz
747 fsav11=fsav11+fni(i)*dt12
748 END IF
749 IF(isensint(1)/=0) THEN
750 IF (ims2 >0 ) THEN
751 fsavparit(1,1,i) = -fxi(i)
752 fsavparit(1,2,i) = -fyi(i)
753 fsavparit(1,3,i) = -fzi(i)
754 ELSE
755 fsavparit(1,1,i) = fxi(i)
756 fsavparit(1,2,i) = fyi(i)
757 fsavparit(1,3,i) = fzi(i)
758 END IF
759 ENDIF
760 ENDDO
761 ELSE
762 DO i=1,jlt
763 IF(pene(i) == zero)cycle
764 fxi(i)=n1(i)*fni(i)
765 fyi(i)=n2(i)*fni(i)
766 fzi(i)=n3(i)*fni(i)
767 impx=fxi(i)*dt12
768 impy=fyi(i)*dt12
769 impz=fzi(i)*dt12
770 fsav1 =fsav1 +impx
771 fsav2 =fsav2 +impy
772 fsav3 =fsav3 +impz
773 fsav8 =fsav8 +abs(impx)
774 fsav9 =fsav9 +abs(impy)
775 fsav10=fsav10+abs(impz)
776 fsav11=fsav11+fni(i)*dt12
777 IF(isensint(1)/=0) THEN
778 fsavparit(1,1,i) = fxi(i)
779 fsavparit(1,2,i) = fyi(i)
780 fsavparit(1,3,i) = fzi(i)
781 ENDIF
782 ENDDO
783 END IF !(ILEV==2) THEN
784#include "lockon.inc"
785 fsav(1)=fsav(1)+fsav1
786 fsav(2)=fsav(2)+fsav2
787 fsav(3)=fsav(3)+fsav3
788 fsav(8)=fsav(8)+fsav8
789 fsav(9)=fsav(9)+fsav9
790 fsav(10)=fsav(10)+fsav10
791 fsav(11)=fsav(11)+fsav11
792#include "lockoff.inc"
793C
794C---------------------------------
795 IF((anim_v(12)+outp_v(12)+h3d_data%N_VECT_PCONT>0.AND.
796 . ((tt>=tanim .AND. tt<=tanim_stop).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
797 . (manim>=4.AND.manim<=15).OR.h3d_data%MH3D/=0))
798 . .OR.h3d_data%N_VECT_PCONT_MAX>0)THEN
799 IF (inconv==1) THEN
800#include "lockon.inc"
801 DO i=1,jlt
802 IF(pene(i) == zero)cycle
803 fncont(1,ix1(i)) =fncont(1,ix1(i)) + fxi(i)*h1(i)
804 fncont(2,ix1(i)) =fncont(2,ix1(i)) + fyi(i)*h1(i)
805 fncont(3,ix1(i)) =fncont(3,ix1(i)) + fzi(i)*h1(i)
806 fncont(1,ix2(i)) =fncont(1,ix2(i)) + fxi(i)*h2(i)
807 fncont(2,ix2(i)) =fncont(2,ix2(i)) + fyi(i)*h2(i)
808 fncont(3,ix2(i)) =fncont(3,ix2(i)) + fzi(i)*h2(i)
809 fncont(1,ix3(i)) =fncont(1,ix3(i)) + fxi(i)*h3(i)
810 fncont(2,ix3(i)) =fncont(2,ix3(i)) + fyi(i)*h3(i)
811 fncont(3,ix3(i)) =fncont(3,ix3(i)) + fzi(i)*h3(i)
812 fncont(1,ix4(i)) =fncont(1,ix4(i)) + fxi(i)*h4(i)
813 fncont(2,ix4(i)) =fncont(2,ix4(i)) + fyi(i)*h4(i)
814 fncont(3,ix4(i)) =fncont(3,ix4(i)) + fzi(i)*h4(i)
815 jg = nsvg(i)
816 IF(jg>0) THEN
817C en SPMD : traitement a refaire apres reception noeud remote si JG < 0
818 fncont(1,jg)=fncont(1,jg)- fxi(i)
819 fncont(2,jg)=fncont(2,jg)- fyi(i)
820 fncont(3,jg)=fncont(3,jg)- fzi(i)
821 ELSE ! cas noeud remote en SPMD
822 jg = -jg
823 fnconti(nin)%P(1,jg)=fnconti(nin)%P(1,jg)-fxi(i)
824 fnconti(nin)%P(2,jg)=fnconti(nin)%P(2,jg)-fyi(i)
825 fnconti(nin)%P(3,jg)=fnconti(nin)%P(3,jg)-fzi(i)
826 ENDIF
827 ENDDO
828#include "lockoff.inc"
829 END IF !(INCONV==1) THEN
830 ENDIF
831C
832C---------------------------------
833C SORTIES TH PAR SOUS INTERFACE
834C---------------------------------
835 IF(nisub/=0)THEN
836 DO jsub=1,nisub
837 DO j=1,25
838 fsavsub1(j,jsub)=zero
839 END DO
840 ENDDO
841 DO i=1,jlt
842 IF(pene(i) == zero)cycle
843 nn = nsvg(i)
844 IF(nn>0)THEN
845 in=cn_loc(i)
846
847 IF (msegtyp(ce_loc(i)) < 0) THEN
848 ie= - msegtyp(ce_loc(i))
849 ELSE
850 ie = ce_loc(i)
851 ENDIF
852 IF(ie > nrtm) ie=ie-nrtm
853
854 jj =addsubs(in)
855 kk =addsubm(ie)
856 DO WHILE(jj<addsubs(in+1))
857 jsub=lisubs(jj)
858
859 itypsub = typsub(jsub)
860 IF(itypsub == 1 ) THEN ! Defining specific inter
861
862 iss1 = bitget(inflg_subs(jj),0)
863 iss2 = bitget(inflg_subs(jj),1)
864 igrn = bitget(inflg_subs(jj),2)
865 ksub=lisubm(kk)
866 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
867 ims1 = bitget(inflg_subm(kk),0)
868 ims2 = bitget(inflg_subm(kk),1)
869 IF(ksub==jsub)THEN
870 IF(.NOT.(((ims1 == 1 .OR. ims2 == 1) .AND. igrn == 1).OR.
871 . (ims1 == 1 .AND. iss2 == 1).OR.
872 . (ims2 == 1 .AND. iss1 == 1))) THEN
873 kk=kk+1
874 ksub=lisubm(kk)
875 cycle
876 END IF
877 impx=fxi(i)*dt12
878 impy=fyi(i)*dt12
879 impz=fzi(i)*dt12
880
881 IF(ims2 > 0)THEN
882 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
883 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
884 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
885 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
886 ELSE
887 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
888 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
889 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
890 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
891 END IF
892C
893 IF(isensint(jsub+1)/=0) THEN
894 IF(ims2 > 0)THEN
895 fsavparit(jsub+1,1,i) = -fxi(i)
896 fsavparit(jsub+1,2,i) = -fyi(i)
897 fsavparit(jsub+1,3,i) = -fzi(i)
898 ELSE
899 fsavparit(jsub+1,1,i) = fxi(i)
900 fsavparit(jsub+1,2,i) = fyi(i)
901 fsavparit(jsub+1,3,i) = fzi(i)
902 END IF
903 ENDIF
904C
905 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
906 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
907 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
908C
909 END IF
910 kk=kk+1
911 ksub=lisubm(kk)
912 ENDDO
913 jj=jj+1
914
915 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
916
917 impx=fxi(i)*dt12
918 impy=fyi(i)*dt12
919 impz=fzi(i)*dt12
920
921 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
922 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
923 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
924
925 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
926 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
927 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
928
929 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
930
931 IF(isensint(jsub+1)/=0) THEN
932 fsavparit(jsub+1,1,i) = fxi(i)
933 fsavparit(jsub+1,2,i) = fyi(i)
934 fsavparit(jsub+1,3,i) = fzi(i)
935 ENDIF
936
937 jj=jj+1
938
939 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfs
940
941 iss2 = bitget(inflg_subs(jj),0)
942 iss1 = bitget(inflg_subs(jj),1)
943 ksub=lisubm(kk)
944 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
945 ims2 = bitget(inflg_subm(kk),0)
946 ims1 = bitget(inflg_subm(kk),1)
947 IF(ksub==jsub)THEN
948 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
949 . (ims2 == 1 .AND. iss1 == 1))) THEN
950 kk=kk+1
951 ksub=lisubm(kk)
952 cycle
953 END IF
954
955 impx=fxi(i)*dt12
956 impy=fyi(i)*dt12
957 impz=fzi(i)*dt12
958
959
960 IF(ims2 > 0)THEN
961 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
962 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
963 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
964 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
965 ELSE
966 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
967 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
968 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
969 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
970 END IF
971
972 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
973 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
974 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
975
976
977 IF(isensint(jsub+1)/=0) THEN
978 IF(ims2 > 0)THEN
979 fsavparit(jsub+1,1,i) = -fxi(i)
980 fsavparit(jsub+1,2,i) = -fyi(i)
981 fsavparit(jsub+1,3,i) = -fzi(i)
982 ELSE
983 fsavparit(jsub+1,1,i) = fxi(i)
984 fsavparit(jsub+1,2,i) = fyi(i)
985 fsavparit(jsub+1,3,i) = fzi(i)
986 END IF
987 ENDIF
988 ENDIF
989 kk=kk+1
990 ksub=lisubm(kk)
991 ENDDO
992 jj=jj+1
993
994 ENDIF
995
996
997 END DO
998 END IF
999
1000 IF (msegtyp(ce_loc(i)) < 0) THEN
1001 ie= - msegtyp(ce_loc(i))
1002 ELSE
1003 ie = ce_loc(i)
1004 ENDIF
1005 IF(ie > nrtm) ie=ie-nrtm
1006
1007 kk =addsubm(ie)
1008 DO WHILE(kk<addsubm(ie+1))
1009 ksub=lisubm(kk)
1010
1011 itypsub = typsub(ksub)
1012 IF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
1013
1014 impx=-fxi(i)*dt12
1015 impy=-fyi(i)*dt12
1016 impz=-fzi(i)*dt12
1017
1018 fsavsub1(1,ksub)=fsavsub1(1,ksub)+impx
1019 fsavsub1(2,ksub)=fsavsub1(2,ksub)+impy
1020 fsavsub1(3,ksub)=fsavsub1(3,ksub)+impz
1021
1022 fsavsub1(8,ksub) =fsavsub1(8,ksub) +abs(impx)
1023 fsavsub1(9,ksub) =fsavsub1(9,ksub) +abs(impy)
1024 fsavsub1(10,ksub)=fsavsub1(10,ksub)+abs(impz)
1025
1026 fsavsub1(11,ksub)=fsavsub1(11,ksub)-fni(i)*dt12
1027
1028 IF(isensint(ksub+1)/=0) THEN
1029 fsavparit(ksub+1,1,i) = -fxi(i)
1030 fsavparit(ksub+1,2,i) = -fyi(i)
1031 fsavparit(ksub+1,3,i) = -fzi(i)
1032 ENDIF
1033
1034
1035 ENDIF
1036 kk=kk+1
1037 ENDDO
1038
1039 END DO
1040 IF(nspmd>1) THEN
1041 DO i=1,jlt
1042 IF(pene(i) == zero)cycle
1043 nn = nsvg(i)
1044 IF(nn<0)THEN
1045 nn = -nn
1046
1047 IF (msegtyp(ce_loc(i)) < 0) THEN
1048 ie= - msegtyp(ce_loc(i))
1049 ELSE
1050 ie = ce_loc(i)
1051 ENDIF
1052 IF(ie > nrtm) ie=ie-nrtm
1053
1054 jj =addsubsfi(nin)%P(nn)
1055 kk =addsubm(ie)
1056 DO WHILE(jj<addsubsfi(nin)%P(nn+1))
1057 jsub=lisubsfi(nin)%P(jj)
1058 itypsub = typsub(jsub)
1059 IF(itypsub == 1 ) THEN
1060
1061 iss1 = bitget(inflg_subsfi(nin)%P(jj),0)
1062 iss2 = bitget(inflg_subsfi(nin)%P(jj),1)
1063 igrn = bitget(inflg_subsfi(nin)%P(jj),2)
1064 ksub=lisubm(kk)
1065 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1066 ims1 = bitget(inflg_subm(kk),0)
1067 ims2 = bitget(inflg_subm(kk),1)
1068 IF(ksub==jsub)THEN
1069 IF(.NOT.(((ims1 == 1 .OR. ims2 == 1) .AND. igrn == 1).OR.
1070 . (ims1 == 1 .AND. iss2 == 1).OR.
1071 . (ims2 == 1 .AND. iss1 == 1))) THEN
1072 kk=kk+1
1073 ksub=lisubm(kk)
1074 cycle
1075 END IF
1076 impx=fxi(i)*dt12
1077 impy=fyi(i)*dt12
1078 impz=fzi(i)*dt12
1079
1080 IF(ims2 > 0)THEN
1081 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
1082 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
1083 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
1084 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
1085 ELSE
1086 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1087 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1088 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1089 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1090 END IF
1091C
1092 IF(isensint(jsub+1)/=0) THEN
1093 IF(ims2 > 0)THEN
1094 fsavparit(jsub+1,1,i) = -fxi(i)
1095 fsavparit(jsub+1,2,i) = -fyi(i)
1096 fsavparit(jsub+1,3,i) = -fzi(i)
1097 ELSE
1098 fsavparit(jsub+1,1,i) = fxi(i)
1099 fsavparit(jsub+1,2,i) = fyi(i)
1100 fsavparit(jsub+1,3,i) = fzi(i)
1101 END IF
1102 ENDIF
1103C
1104 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1105 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1106 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1107C
1108 ENDIF
1109 kk=kk+1
1110 ksub=lisubm(kk)
1111 ENDDO
1112 jj=jj+1
1113
1114 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
1115
1116 impx=fxi(i)*dt12
1117 impy=fyi(i)*dt12
1118 impz=fzi(i)*dt12
1119
1120 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1121 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1122 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1123
1124 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1125 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1126 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1127
1128 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1129
1130 IF(isensint(jsub+1)/=0) THEN
1131 fsavparit(jsub+1,1,i) = fxi(i)
1132 fsavparit(jsub+1,2,i) = fyi(i)
1133 fsavparit(jsub+1,3,i) = fzi(i)
1134 ENDIF
1135
1136 jj=jj+1
1137
1138 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfaces
1139
1140 iss2 = bitget(inflg_subsfi(nin)%P(jj),0)
1141 iss1 = bitget(inflg_subsfi(nin)%P(jj),1)
1142 ksub=lisubm(kk)
1143 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1144 ims2 = bitget(inflg_subm(kk),0)
1145 ims1 = bitget(inflg_subm(kk),1)
1146 IF(ksub==jsub)THEN
1147 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
1148 . (ims2 == 1 .AND. iss1 == 1))) THEN
1149 kk=kk+1
1150 ksub=lisubm(kk)
1151 cycle
1152 END IF
1153
1154 impx=fxi(i)*dt12
1155 impy=fyi(i)*dt12
1156 impz=fzi(i)*dt12
1157
1158 IF(ims2 > 0)THEN
1159 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
1160 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
1161 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
1162 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
1163 ELSE
1164 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1165 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1166 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1167 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1168 END IF
1169C
1170 IF(isensint(jsub+1)/=0) THEN
1171 IF(ims2 > 0)THEN
1172 fsavparit(jsub+1,1,i) = -fxi(i)
1173 fsavparit(jsub+1,2,i) = -fyi(i)
1174 fsavparit(jsub+1,3,i) = -fzi(i)
1175 ELSE
1176 fsavparit(jsub+1,1,i) = fxi(i)
1177 fsavparit(jsub+1,2,i) = fyi(i)
1178 fsavparit(jsub+1,3,i) = fzi(i)
1179 END IF
1180 ENDIF
1181
1182 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1183 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1184 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1185
1186 ENDIF
1187 kk=kk+1
1188 ksub=lisubm(kk)
1189 ENDDO
1190 jj=jj+1
1191
1192 ENDIF
1193
1194 END DO
1195 END IF
1196 END DO
1197 END IF
1198 END IF
1199
1200C------------For /LOAD/PRESSURE tag nodes in contact-------------
1201 IF(ninloadp > 0) THEN
1202 DO k = kloadpinter(nin)+1, kloadpinter(nin+1)
1203 pp = loadpinter(k)
1204 ppl = loadp_hyd_inter(pp)
1205 dgapload = dgaploadint(k)
1206 DO i=1,jlt
1207 gapp= gapv(i) + dgapload
1208 IF(pene(i) > zero .OR.(dist(i) <= gapp)) THEN
1209 tagncont(ppl,ix1(i)) = 1
1210 tagncont(ppl,ix2(i)) = 1
1211 tagncont(ppl,ix3(i)) = 1
1212 tagncont(ppl,ix4(i)) = 1
1213 jg = nsvg(i)
1214 IF(jg>0) THEN
1215C SPMD : do same after reception of forces for remote nodes
1216 tagncont(ppl,jg) = 1
1217 ENDIF
1218 ENDIF
1219 ENDDO
1220 ENDDO
1221
1222 ENDIF
1223
1224C---------------------------------
1225C NEW FRICTION MODELS
1226C---------------------------------
1227
1228C Friction coefficient computation
1229C++++++++++++++++++++++++++++++++++++++++
1230
1231 IF(ivis2/=-1) THEN ! friction calculation not needed for adhesion case
1232
1233 IF(iorthfric == 0) THEN
1234C++ Isotropic Friction
1235
1236 IF (mfrot==0) THEN
1237C--- Coulomb friction
1238 DO i=1,jlt
1239 xmu(i) = fricc(i)
1240 ENDDO
1241 ELSEIF (mfrot==1) THEN
1242C--- Viscous friction
1243 DO i=1,jlt
1244 ie=ce_loc(i)
1245 IF(pene(i) == 0) THEN
1246 xmu(i) = zero
1247 cycle
1248 ENDIF
1249 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1250 v2 = (vx(i) - n1(i)*aa)**2
1251 . + (vy(i) - n2(i)*aa)**2
1252 . + (vz(i) - n3(i)*aa)**2
1253 vv = sqrt(max(em30,v2))
1254 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1255 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1256 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1257 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1258 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1259 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1260 ax = ay1*az2 - az1*ay2
1261 ay = az1*ax2 - ax1*az2
1262 az = ax1*ay2 - ay1*ax2
1263 area = half*sqrt(ax*ax+ay*ay+az*az)
1264 p = -fni(i)/area
1265 xmu(i) = fricc(i) + (fric_coefs(i,1) + fric_coefs(i,4)*p ) * p
1266 . +(fric_coefs(i,2) + fric_coefs(i,3)*p) * vv + fric_coefs(i,5)*v2
1267 xmu(i) = max(xmu(i),em30)
1268 ENDDO
1269 ELSEIF(mfrot==2)THEN
1270C--- Darmstad LAW
1271 DO i=1,jlt
1272 ie=ce_loc(i)
1273 IF(pene(i) == 0) THEN
1274 xmu(i) = zero
1275 cycle
1276 ENDIF
1277 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1278 v2 = (vx(i) - n1(i)*aa)**2
1279 . + (vy(i) - n2(i)*aa)**2
1280 . + (vz(i) - n3(i)*aa)**2
1281 vv = sqrt(max(em30,v2))
1282 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1283 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1284 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1285 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1286 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1287 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1288 ax = ay1*az2 - az1*ay2
1289 ay = az1*ax2 - ax1*az2
1290 az = ax1*ay2 - ay1*ax2
1291 area = half*sqrt(ax*ax+ay*ay+az*az)
1292 p = -fni(i)/area
1293 xmu(i) = fricc(i)
1294 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
1295 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
1296 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
1297 xmu(i) = max(xmu(i),em30)
1298 ENDDO
1299 ELSEIF (mfrot==3) THEN
1300C--- Renard LAW
1301 DO i=1,jlt
1302 IF(pene(i) == 0) THEN
1303 xmu(i) = zero
1304 cycle
1305 ENDIF
1306 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1307 v2 = (vx(i) - n1(i)*aa)**2
1308 . + (vy(i) - n2(i)*aa)**2
1309 . + (vz(i) - n3(i)*aa)**2
1310 vv = sqrt(max(em30,v2))
1311 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
1312 dmu = fric_coefs(i,3)-fric_coefs(i,1)
1313 vv1 = vv / fric_coefs(i,5)
1314 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
1315 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
1316 dmu = fric_coefs(i,4)-fric_coefs(i,3)
1317 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
1318 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
1319 ELSE
1320 dmu = fric_coefs(i,2)-fric_coefs(i,4)
1321 vv2 = (vv - fric_coefs(i,6))**2
1322 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
1323 ENDIF
1324 xmu(i) = max(xmu(i),em30)
1325 ENDDO
1326 ELSEIF(mfrot==4)THEN
1327C--- Exponential decay model
1328 DO i=1,jlt
1329 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1330 v2 = (vx(i) - n1(i)*aa)**2
1331 . + (vy(i) - n2(i)*aa)**2
1332 . + (vz(i) - n3(i)*aa)**2
1333 vv = sqrt(max(em30,v2))
1334 xmu(i) = fric_coefs(i,1)
1335 . + (fricc(i)-fric_coefs(i,1))*exp(-fric_coefs(i,2)*vv)
1336 xmu(i) = max(xmu(i),em30)
1337 ENDDO
1338 ENDIF
1339
1340 ELSE
1341C++ Orthotropic Friction
1342
1343 IF (mfrot==0) THEN
1344C--- Coulomb friction
1345#include "vectorize.inc"
1346 DO k=1,nfisot
1347 i = indexisot(k)
1348 xmu(i) = fricc(i)
1349 ENDDO
1350#include "vectorize.inc"
1351 DO k=1,nforth
1352 i = indexorth(k)
1353 xmu(i) = fricc(i)
1354 xmu2(i) = fricc2(i)
1355 IF(xmu(i)<=em30) THEN
1356 xmu(i) = em30
1357 dir1(i,1:3) = zero
1358 an(k) = zero
1359 ELSE
1360 an(k)=xmu(i)**2
1361 an(k)=one/an(k)
1362 ENDIF
1363 IF(xmu2(i)<=em30) THEN
1364 xmu2(i) = em30
1365 dir2(i,1:3) = zero
1366 bn(k) = zero
1367 ELSE
1368 bn(k)=xmu2(i)**2
1369 bn(k)=one/bn(k)
1370 ENDIF
1371 ENDDO
1372 ELSEIF (mfrot==1) THEN
1373C--- Viscous friction
1374#include "vectorize.inc"
1375 DO k=1,nfisot ! isotropic friction couples
1376 i = indexisot(k)
1377 IF(pene(i) == 0) THEN
1378 xmu(i) = zero
1379 cycle
1380 ENDIF
1381 ie=ce_loc(i)
1382 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1383 v2 = (vx(i) - n1(i)*aa)**2
1384 . + (vy(i) - n2(i)*aa)**2
1385 . + (vz(i) - n3(i)*aa)**2
1386 vv = sqrt(max(em30,v2))
1387
1388 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1389 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1390 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1391 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1392 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1393 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1394 ax = ay1*az2 - az1*ay2
1395 ay = az1*ax2 - ax1*az2
1396 az = ax1*ay2 - ay1*ax2
1397 area = half*sqrt(ax*ax+ay*ay+az*az)
1398 p = -fni(i)/area
1399 xmu(i) = fricc(i) + (fric_coefs(i,1) + fric_coefs(i,4)*p ) * p
1400 . +(fric_coefs(i,2) + fric_coefs(i,3)*p) * vv + fric_coefs(i,5)*v2
1401
1402 xmu(i) = max(xmu(i),em30)
1403 ENDDO
1404
1405#include "vectorize.inc"
1406 DO k=1,nforth ! Orthotropic friction couples
1407 i = indexorth(k)
1408 ie=ce_loc(i)
1409 IF(pene(i) == 0) THEN
1410 xmu(i) = zero
1411 xmu2(i) = zero
1412 cycle
1413 ENDIF
1414
1415
1416c
1417 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1418 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1419 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1420 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1421 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1422 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1423 ax = ay1*az2 - az1*ay2
1424 ay = az1*ax2 - ax1*az2
1425 az = ax1*ay2 - ay1*ax2
1426 area = half*sqrt(ax*ax+ay*ay+az*az)
1427 p = -fni(i)/area
1428c
1429 v2 = vx(i)*dir1(k,1) +vy(i)*dir1(k,2)+vz(i)*dir1(k,3)
1430 vv = max(em30,v2)
1431 xmu(i) = fricc(i) + (fric_coefs(i,1) + fric_coefs(i,4)*p ) * p
1432 . +(fric_coefs(i,2) + fric_coefs(i,3)*p) * vv + fric_coefs(i,5)*vv*vv
1433
1434 v2 = vx(i)*dir2(k,1) +vy(i)*dir2(k,2)+vz(i)*dir2(k,3)
1435 vv = max(em30,v2)
1436 xmu2(i) = fricc2(i) + (fric_coefs2(i,1) + fric_coefs2(i,4)*p ) * p
1437 . +(fric_coefs2(i,2) + fric_coefs2(i,3)*p) * vv + fric_coefs2(i,5)*vv*vv
1438
1439 xmu(i) = max(xmu(i),em30)
1440 xmu2(i) = max(xmu2(i),em30)
1441 ENDDO
1442
1443#include "vectorize.inc"
1444 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
1445 i = indexorth(k)
1446 IF(xmu(i)<=em30) THEN
1447 xmu(i) = em30
1448 dir1(i,1:3) = zero
1449 an(k) = zero
1450 ELSE
1451 an(k)=xmu(i)**2
1452 an(k)=one/an(k)
1453 ENDIF
1454 IF(xmu2(i)<=em30) THEN
1455 xmu2(i) = em30
1456 dir2(i,1:3) = zero
1457 bn(k) = zero
1458 ELSE
1459 bn(k)=xmu2(i)**2
1460 bn(k)=one/bn(k)
1461 ENDIF
1462 ENDDO
1463
1464 ELSEIF(mfrot==2)THEN
1465C--- Loi Darmstad
1466#include "vectorize.inc"
1467 DO k=1,nfisot ! isotropic friction couples
1468 i = indexisot(k)
1469 IF(pene(i) == 0) THEN
1470 xmu(i) = zero
1471 cycle
1472 ENDIF
1473
1474 ie=ce_loc(i)
1475 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1476 v2 = (vx(i) - n1(i)*aa)**2
1477 . + (vy(i) - n2(i)*aa)**2
1478 . + (vz(i) - n3(i)*aa)**2
1479 vv = sqrt(max(em30,v2))
1480 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1481 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1482 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1483 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1484 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1485 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1486 ax = ay1*az2 - az1*ay2
1487 ay = az1*ax2 - ax1*az2
1488 az = ax1*ay2 - ay1*ax2
1489 area = half*sqrt(ax*ax+ay*ay+az*az)
1490 p = -fni(i)/area
1491 xmu(i) = fricc(i)
1492 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
1493 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
1494 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
1495 xmu(i) = max(xmu(i),em30)
1496 ENDDO
1497c
1498#include "vectorize.inc"
1499 DO k=1,nforth ! Orthotropic friction couples
1500 i = indexorth(k)
1501 ie=ce_loc(i)
1502 IF(pene(i) == 0) THEN
1503 xmu(i) = zero
1504 xmu2(i) = zero
1505 cycle
1506 ENDIF
1507
1508c
1509 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1510 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1511 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1512 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1513 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1514 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1515 ax = ay1*az2 - az1*ay2
1516 ay = az1*ax2 - ax1*az2
1517 az = ax1*ay2 - ay1*ax2
1518 area = half*sqrt(ax*ax+ay*ay+az*az)
1519 p = -fni(i)/area
1520c
1521 v2 = vx(i)*dir1(k,1) +vy(i)*dir1(k,2)+vz(i)*dir1(k,3)
1522 vv = max(em30,v2)
1523 xmu(i) = fricc(i)
1524 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
1525 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
1526 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
1527c
1528 v2 = vx(i)*dir2(k,1) +vy(i)*dir2(k,2)+vz(i)*dir2(k,3)
1529 vv = max(em30,v2)
1530 xmu2(i) = fricc2(i)
1531 . + fric_coefs2(i,1)*exp(fric_coefs2(i,2)*vv)*p*p
1532 . + fric_coefs2(i,3)*exp(fric_coefs2(i,4)*vv)*p
1533 . + fric_coefs2(i,5)*exp(fric_coefs2(i,6)*vv)
1534
1535 xmu(i) = max(xmu(i),em30)
1536 xmu2(i) = max(xmu2(i),em30)
1537 ENDDO
1538
1539#include "vectorize.inc"
1540 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
1541 i = indexorth(k)
1542 IF(xmu(i)<=em30) THEN
1543 xmu(i) = em30
1544 dir1(i,1:3) = zero
1545 an(k) = zero
1546 ELSE
1547 an(k)=xmu(i)**2
1548 an(k)=one/an(k)
1549 ENDIF
1550 IF(xmu2(i)<=em30) THEN
1551 xmu2(i) = em30
1552 dir2(i,1:3) = zero
1553 bn(k) = zero
1554 ELSE
1555 bn(k)=xmu2(i)**2
1556 bn(k)=one/bn(k)
1557 ENDIF
1558 ENDDO
1559
1560 ELSEIF (mfrot==3) THEN
1561C--- Loi Renard
1562#include "vectorize.inc"
1563 DO k=1,nfisot ! isotropic friction couples
1564 i = indexisot(k)
1565 IF(pene(i) == 0) THEN
1566 xmu(i) = zero
1567 cycle
1568 ENDIF
1569
1570 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1571 v2 = (vx(i) - n1(i)*aa)**2
1572 . + (vy(i) - n2(i)*aa)**2
1573 . + (vz(i) - n3(i)*aa)**2
1574 vv = sqrt(max(em30,v2))
1575 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
1576 dmu = fric_coefs(i,3)-fric_coefs(i,1)
1577 vv1 = vv / fric_coefs(i,5)
1578 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
1579 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
1580 dmu = fric_coefs(i,4)-fric_coefs(i,3)
1581 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
1582 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
1583 ELSE
1584 dmu = fric_coefs(i,2)-fric_coefs(i,4)
1585 vv2 = (vv - fric_coefs(i,6))**2
1586 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
1587 ENDIF
1588 xmu(i) = max(xmu(i),em30)
1589 ENDDO
1590
1591#include "vectorize.inc"
1592 DO k=1,nforth ! Orthotropic friction couples
1593 i = indexorth(k)
1594 IF(pene(i) == 0) THEN
1595 xmu(i) = zero
1596 xmu2(i) = zero
1597 cycle
1598 ENDIF
1599
1600c
1601 v2 = vx(i)*dir1(k,1) +vy(i)*dir1(k,2)+vz(i)*dir1(k,3)
1602 vv = max(em30,v2)
1603 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
1604 dmu = fric_coefs(i,3)-fric_coefs(i,1)
1605 vv1 = vv / fric_coefs(i,5)
1606 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
1607 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
1608 dmu = fric_coefs(i,4)-fric_coefs(i,3)
1609 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
1610 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
1611 ELSE
1612 dmu = fric_coefs(i,2)-fric_coefs(i,4)
1613 vv2 = (vv - fric_coefs(i,6))**2
1614 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
1615 ENDIF
1616
1617 v2 = vx(i)*dir2(k,1) +vy(i)*dir2(k,2)+vz(i)*dir2(k,3)
1618 vv = max(em30,v2)
1619 IF(vv>=0.AND.vv<=fric_coefs2(i,5)) THEN
1620 dmu = fric_coefs2(i,3)-fric_coefs2(i,1)
1621 vv1 = vv / fric_coefs2(i,5)
1622 xmu2(i) = fric_coefs2(i,1)+ dmu*vv1*(two-vv1)
1623 ELSEIF(vv>fric_coefs2(i,5).AND.vv<fric_coefs2(i,6)) THEN
1624 dmu = fric_coefs2(i,4)-fric_coefs2(i,3)
1625 vv1 = (vv - fric_coefs2(i,5))/(fric_coefs2(i,6)-fric_coefs2(i,5))
1626 xmu2(i) = fric_coefs2(i,3)+ dmu * (three-two*vv1)*vv1**2
1627 ELSE
1628 dmu = fric_coefs2(i,2)-fric_coefs2(i,4)
1629 vv2 = (vv - fric_coefs2(i,6))**2
1630 xmu2(i) = fric_coefs2(i,2) - dmu / (one + dmu*vv2)
1631 ENDIF
1632
1633 xmu(i) = max(xmu(i),em30)
1634 xmu2(i) = max(xmu2(i),em30)
1635
1636 ENDDO
1637
1638#include "vectorize.inc"
1639 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
1640 i = indexorth(k)
1641 IF(xmu(i)<=em30) THEN
1642 xmu(i) = em30
1643 dir1(i,1:3) = zero
1644 an(k) = zero
1645 ELSE
1646 an(k)=xmu(i)**2
1647 an(k)=one/an(k)
1648 ENDIF
1649 IF(xmu2(i)<=em30) THEN
1650 xmu2(i) = em30
1651 dir2(i,1:3) = zero
1652 bn(k) = zero
1653 ELSE
1654 bn(k)=xmu2(i)**2
1655 bn(k)=one/bn(k)
1656 ENDIF
1657 ENDDO
1658
1659 ELSEIF(mfrot==4)THEN
1660C--- Exponential decay model
1661#include "vectorize.inc"
1662 DO k=1,nfisot ! isotropic friction couples
1663 i = indexisot(k)
1664 IF(pene(i) == 0) THEN
1665 xmu(i) = zero
1666 xmu2(i) = zero
1667 cycle
1668 ENDIF
1669 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1670 v2 = (vx(i) - n1(i)*aa)**2
1671 . + (vy(i) - n2(i)*aa)**2
1672 . + (vz(i) - n3(i)*aa)**2
1673 vv = sqrt(max(em30,v2))
1674 xmu(i) = fric_coefs(i,1)
1675 . + (fricc(i)-fric_coefs(i,1))*exp(-fric_coefs(i,2)*vv)
1676 xmu(i) = max(xmu(i),em30)
1677 ENDDO
1678c
1679#include "vectorize.inc"
1680 DO k=1,nforth ! Orthotropic friction couples
1681 i = indexorth(k)
1682 IF(pene(i) == 0) THEN
1683 xmu(i) = zero
1684 xmu2(i) = zero
1685 cycle
1686 ENDIF
1687 ie=ce_loc(i)
1688c
1689 v2 = vx(i)*dir1(i,1) +vy(i)*dir1(i,2)+vz(i)*dir1(i,3)
1690 vv = max(em30,v2)
1691 xmu(i) = fricc(i)
1692 . + (fric_coefs(i,1)-fricc(i))*exp(-fric_coefs(i,2)*vv)
1693c
1694 v2 = vx(i)*dir2(i,1) +vy(i)*dir2(i,2)+vz(i)*dir2(i,3)
1695 vv = max(em30,v2)
1696 xmu2(i) = fric_coefs2(i,1)
1697 . + (fricc2(i)-fric_coefs2(i,1))*exp(-fric_coefs2(i,2)*vv)
1698 xmu(i) = max(xmu(i),em30)
1699 xmu2(i) = max(xmu2(i),em30)
1700 ENDDO
1701
1702#include "vectorize.inc"
1703 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
1704 i = indexorth(k)
1705 IF(xmu(i)<=em30) THEN
1706 xmu(i) = em30
1707 dir1(i,1:3) = zero
1708 an(k) = zero
1709 ELSE
1710 an(k)=xmu(i)**2
1711 an(k)=one/an(k)
1712 ENDIF
1713 IF(xmu2(i)<=em30) THEN
1714 xmu2(i) = em30
1715 dir2(i,1:3) = zero
1716 bn(k) = zero
1717 ELSE
1718 bn(k)=xmu2(i)**2
1719 bn(k)=one/bn(k)
1720 ENDIF
1721 ENDDO
1722 ENDIF
1723
1724 ENDIF ! IORTHFRIC
1725 ENDIF
1726C------------------
1727C TANGENT FORCE CALCULATION
1728C------------------
1729 fxt(1:jlt)=zero
1730 fyt(1:jlt)=zero
1731 fzt(1:jlt)=zero
1732C
1733 IF(ivis2==-1)THEN ! ADHESION CASE
1734 DO i=1,jlt
1735 IF(pene(i)==zero)cycle
1736 vnx = n1(i)*vn(i)
1737 vny = n2(i)*vn(i)
1738 vnz = n3(i)*vn(i)
1739 vx(i) = vx(i) - vnx
1740 vy(i) = vy(i) - vny
1741 vz(i) = vz(i) - vnz
1742C SHEAR STRESS = VISCOSITY * DU/DY - NEWTONIAN SHEAR LAW
1743 factor = viscadhfact*viscfluid*two/gapv(i)*areas(i)
1744 fxt(i) = factor*vx(i)
1745 fyt(i) = factor*vy(i)
1746 fzt(i) = factor*vz(i)
1747C------- total force and energy
1748 fxi(i) = fxi(i) + fxt(i)
1749 fyi(i) = fyi(i) + fyt(i)
1750 fzi(i) = fzi(i) + fzt(i)
1751 econtdt = econtdt + dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i)) ! Tangential contact energy
1752 ENDDO
1753 ELSE ! NO ADHESION CASE
1754 IF (ifq /= 0) THEN
1755C---------------------------------
1756C INCREMENTAL (STIFFNESS) FORMULATION
1757C---------------------------------
1758 IF (ifq==13) THEN
1759 alpha = max(one,alpha0*dt12)
1760 ELSE
1761 alpha = alpha0
1762 ENDIF
1763
1764 IF(iorthfric ==0 ) THEN
1765C++ Isotropic Friction
1766
1767 IF (inconv==1) THEN
1768 DO i=1,jlt
1769 IF(pene(i) == zero)cycle
1770 fx = stif0(i)*vx(i)*dt12
1771 fy = stif0(i)*vy(i)*dt12
1772 fz = stif0(i)*vz(i)*dt12
1773 jg = nsvg(i)
1774 IF(jg>0) THEN
1775 n = cand_n_n(i)
1776c SECND_FR(4:6,N) is the old friction force
1777c if(itab(jg)==31774.or.itab(jg)==6992.or.itab(jg)==7106)
1778c . print *,'slav_fr nat',itab(jg),SECND_FR(4:6,N),PENE_OLD(1:5,N),STIF_OLD(1:2,N)
1779 fx = secnd_fr(4,n) + alpha*fx
1780 fy = secnd_fr(5,n) + alpha*fy
1781 fz = secnd_fr(6,n) + alpha*fz
1782 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
1783 fx = fx - ftn*n1(i)
1784 fy = fy - ftn*n2(i)
1785 fz = fz - ftn*n3(i)
1786 ft = fx*fx + fy*fy + fz*fz
1787 ft = max(ft,tiny)
1788 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
1789 beta = min(one,xmu(i)*sqrt(fn/ft))
1790
1791 fxt(i) = fx * beta
1792 fyt(i) = fy * beta
1793 fzt(i) = fz * beta
1794C
1795 secnd_fr(1,n) = fxt(i)
1796 secnd_fr(2,n) = fyt(i)
1797 secnd_fr(3,n) = fzt(i)
1798
1799C
1800c if(itab(jg)==7444)
1801c . print *,'i25for3 natif',itab(jg),pene(i),fxt(i),fyt(i),fzt(i),fxi(i),fyi(i),fzi(i)
1802 ELSE ! cas noeud remote en SPMD
1803 jg = -jg
1804c SECND_FRFI(NIN)%P(4:6,JG) is the old friction force
1805c if(itafi(nin)%p(jg)==31774.or.itafi(nin)%p(jg)==6992.or.itafi(nin)%p(jg)==7106)
1806c . print *,'slav_fr rem',itafi(nin)%p(jg),SECND_FRFI(NIN)%P(4:6,JG),PENE_OLDFI(nin)%p(1:5,JG),STIF_OLDFI(nin)%p(1:2,JG)
1807 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
1808 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
1809 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
1810 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
1811 fx = fx - ftn*n1(i)
1812 fy = fy - ftn*n2(i)
1813 fz = fz - ftn*n3(i)
1814 ft = fx*fx + fy*fy + fz*fz
1815 ft = max(ft,tiny)
1816 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
1817 beta = min(one,xmu(i)*sqrt(fn/ft))
1818 fxt(i) = fx * beta
1819 fyt(i) = fy * beta
1820 fzt(i) = fz * beta
1821C
1822 secnd_frfi(nin)%P(1,jg) = fxt(i)
1823 secnd_frfi(nin)%P(2,jg) = fyt(i)
1824 secnd_frfi(nin)%P(3,jg) = fzt(i)
1825C
1826c if(itafi(nin)%p(jg)==27715)
1827c . print *,'i25for3 remote',itafi(nin)%p(jg),pene(i),fxt(i),fyt(i),fzt(i),fxi(i),fyi(i),fzi(i)
1828 ENDIF
1829C------- total force
1830 fxi(i)=fxi(i) + fxt(i)
1831 fyi(i)=fyi(i) + fyt(i)
1832 fzi(i)=fzi(i) + fzt(i)
1833
1834 IF( intth > 0 .AND.beta/=zero) THEN
1835 efrict(i) = (fx-fxt(i))*fxt(i) + (fy-fyt(i))*fyt(i) +
1836 . (fz-fzt(i))*fzt(i)
1837 efrict(i) = efrict(i)/stif0(i) ! FRICTIONAL ENERGY
1838 qfrict = qfrict + efrict(i)
1839 ENDIF
1840 efric_l(i) = dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
1841 econvt = econvt + efric_l(i)
1842
1843 ENDDO
1844C--------implicit non converge---
1845 ELSE
1846 DO i=1,jlt
1847 IF(pene(i) == zero)cycle
1848 fx = stif0(i)*vx(i)*dt12
1849 fy = stif0(i)*vy(i)*dt12
1850 fz = stif0(i)*vz(i)*dt12
1851 jg = nsvg(i)
1852 n = cand_n_n(i)
1853 IF(jg>0) THEN
1854c SECND_FR(1:3,N) is the old friction force
1855 fx = secnd_fr(4,n) + alpha*fx
1856 fy = secnd_fr(5,n) + alpha*fy
1857 fz = secnd_fr(6,n) + alpha*fz
1858 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
1859 fx = fx - ftn*n1(i)
1860 fy = fy - ftn*n2(i)
1861 fz = fz - ftn*n3(i)
1862 ft = fx*fx + fy*fy + fz*fz
1863 ft = max(ft,tiny)
1864 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
1865 beta = min(one,xmu(i)*sqrt(fn/ft))
1866 fxt(i) = fx * beta
1867 fyt(i) = fy * beta
1868 fzt(i) = fz * beta
1869 fxi(i)=fxi(i) + fxt(i)
1870 fyi(i)=fyi(i) + fyt(i)
1871 fzi(i)=fzi(i) + fzt(i)
1872 ELSE ! cas noeud remote en SPMD
1873 jg = -jg
1874 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
1875 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
1876 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
1877 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
1878 fx = fx - ftn*n1(i)
1879 fy = fy - ftn*n2(i)
1880 fz = fz - ftn*n3(i)
1881 ft = fx*fx + fy*fy + fz*fz
1882 ft = max(ft,tiny)
1883 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
1884 beta = min(one,xmu(i)*sqrt(fn/ft))
1885 fxt(i) = fx * beta
1886 fyt(i) = fy * beta
1887 fzt(i) = fz * beta
1888 fxi(i)=fxi(i) + fxt(i)
1889 fyi(i)=fyi(i) + fyt(i)
1890 fzi(i)=fzi(i) + fzt(i)
1891 ENDIF
1892C IFPEN(INDEX(I)) = 1
1893 ENDDO
1894 ENDIF
1895
1896
1897 ELSE
1898C++ Orthotropic Friction
1899
1900 IF (inconv==1) THEN
1901#include "vectorize.inc"
1902 DO k=1,nfisot ! isotropic friction couples
1903 i = indexisot(k)
1904 IF(pene(i) == zero)cycle
1905 fx = stif0(i)*vx(i)*dt12
1906 fy = stif0(i)*vy(i)*dt12
1907 fz = stif0(i)*vz(i)*dt12
1908 jg = nsvg(i)
1909 IF(jg>0) THEN
1910 n = cand_n_n(i)
1911c SECND_FR(4:6,N) is the old friction force
1912c if(itab(jg)==31774.or.itab(jg)==6992.or.itab(jg)==7106)
1913c . print *,'slav_fr nat',itab(jg),SECND_FR(4:6,N),PENE_OLD(1:5,N),STIF_OLD(1:2,N)
1914 fx = secnd_fr(4,n) + alpha*fx
1915 fy = secnd_fr(5,n) + alpha*fy
1916 fz = secnd_fr(6,n) + alpha*fz
1917 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
1918 fx = fx - ftn*n1(i)
1919 fy = fy - ftn*n2(i)
1920 fz = fz - ftn*n3(i)
1921 ft = fx*fx + fy*fy + fz*fz
1922 ft = max(ft,tiny)
1923 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
1924 beta = min(one,xmu(i)*sqrt(fn/ft))
1925 fxt(i) = fx * beta
1926 fyt(i) = fy * beta
1927 fzt(i) = fz * beta
1928C
1929 secnd_fr(1,n) = fxt(i)
1930 secnd_fr(2,n) = fyt(i)
1931 secnd_fr(3,n) = fzt(i)
1932C
1933c if(itab(jg)==7444)
1934c . print *,'i25for3 natif',itab(jg),pene(i),fxt(i),fyt(i),fzt(i),fxi(i),fyi(i),fzi(i)
1935 ELSE ! cas noeud remote en SPMD
1936 jg = -jg
1937c SECND_FRFI(NIN)%P(4:6,JG) is the old friction force
1938c if(itafi(nin)%p(jg)==31774.or.itafi(nin)%p(jg)==6992.or.itafi(nin)%p(jg)==7106)
1939c . print *,'slav_fr rem',itafi(nin)%p(jg),SECND_FRFI(NIN)%P(4:6,JG),PENE_OLDFI(nin)%p(1:5,JG),STIF_OLDFI(nin)%p(1:2,JG)
1940 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
1941 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
1942 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
1943 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
1944 fx = fx - ftn*n1(i)
1945 fy = fy - ftn*n2(i)
1946 fz = fz - ftn*n3(i)
1947 ft = fx*fx + fy*fy + fz*fz
1948 ft = max(ft,tiny)
1949 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
1950 beta = min(one,xmu(i)*sqrt(fn/ft))
1951 fxt(i) = fx * beta
1952 fyt(i) = fy * beta
1953 fzt(i) = fz * beta
1954C
1955 secnd_frfi(nin)%P(1,jg) = fxt(i)
1956 secnd_frfi(nin)%P(2,jg) = fyt(i)
1957 secnd_frfi(nin)%P(3,jg) = fzt(i)
1958C
1959c if(itafi(nin)%p(jg)==27715)
1960c . print *,'i25for3 remote',itafi(nin)%p(jg),pene(i),fxt(i),fyt(i),fzt(i),fxi(i),fyi(i),fzi(i)
1961 ENDIF
1962C------- total force
1963 fxi(i)=fxi(i) + fxt(i)
1964 fyi(i)=fyi(i) + fyt(i)
1965 fzi(i)=fzi(i) + fzt(i)
1966
1967 IF( intth > 0 .AND.beta/=zero) THEN
1968 efrict(i) = (fx-fxt(i))*fxt(i) + (fy-fyt(i))*fyt(i) +
1969 . (fz-fzt(i))*fzt(i)
1970 efrict(i) = efrict(i)/stif0(i) ! FRICTIONAL ENERGY
1971 qfrict = qfrict + efrict(i)
1972 ENDIF
1973 efric_l(i) = dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
1974 econvt = econvt + efric_l(i)
1975
1976 ENDDO
1977
1978#include "vectorize.inc"
1979 DO k=1,nforth ! Orthotropic friction couples
1980 i = indexorth(k)
1981 IF(pene(i) == zero)cycle
1982 fx = stif0(i)*vx(i)*dt12
1983 fy = stif0(i)*vy(i)*dt12
1984 fz = stif0(i)*vz(i)*dt12
1985 jg = nsvg(i)
1986 IF(jg>0) THEN
1987 n = cand_n_n(i)
1988c SECND_FR(4:6,N) is the old friction force
1989c if(itab(jg)==31774.or.itab(jg)==6992.or.itab(jg)==7106)
1990c . print *,'slav_fr nat',itab(jg),SECND_FR(4:6,N),PENE_OLD(1:5,N),STIF_OLD(1:2,N)
1991 fx = secnd_fr(4,n) + alpha*fx
1992 fy = secnd_fr(5,n) + alpha*fy
1993 fz = secnd_fr(6,n) + alpha*fz
1994 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
1995 fx = fx - ftn*n1(i)
1996 fy = fy - ftn*n2(i)
1997 fz = fz - ftn*n3(i)
1998
1999 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2000 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2001 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2002 ft = max(ft,em30)
2003 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2004
2005 beta = fn/ft
2006
2007 IF(beta == 0 ) THEN
2008 fxt(i) = zero
2009 fyt(i) = zero
2010 fzt(i) = zero
2011 ELSEIF(beta > 1) THEN ! Inside the ellipse
2012 fxt(i) = fx
2013 fyt(i) = fy
2014 fzt(i) = fz
2015 ELSE ! outside the ellipse
2016
2017! Projection on local tangent of ellipse (outside of ellipse)
2018! ANS = (Fadh-Fproj).n
2019 nep1 =ftt1*an(k)/fn
2020 nep2 =ftt2*bn(k)/fn
2021 nep =nep1*nep1+nep2*nep2
2022 nep =sqrt(nep)
2023
2024 ep=nep1*ftt1+nep2*ftt2
2025
2026 ans=(ep-sqrt(ep))/max(em20,nep)
2027 nep1 =nep1/max(em20,nep)
2028 nep2 =nep2/max(em20,nep)
2029
2030! Projection on ellipse
2031 c11 =ftt1-ans*nep1
2032 c22 =ftt2-ans*nep2
2033
2034 alphaf = atan(c22/c11)
2035
2036 signc = ftt1/max(em20,abs(ftt1))
2037 csa = signc*abs(cos(alphaf))
2038 signc = ftt2/max(em20,abs(ftt2))
2039 sna = signc*abs(sin(alphaf))
2040! Ft computation
2041 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2042 ftt1 = ft * csa
2043 ftt2 = ft * sna
2044
2045 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2046 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2047 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2048
2049 ENDIF
2050C
2051 secnd_fr(1,n) = fxt(i)
2052 secnd_fr(2,n) = fyt(i)
2053 secnd_fr(3,n) = fzt(i)
2054C
2055c if(itab(jg)==7444)
2056c . print *,'i25for3 natif',itab(jg),pene(i),fxt(i),fyt(i),fzt(i),fxi(i),fyi(i),fzi(i)
2057 ELSE ! cas noeud remote en SPMD
2058 jg = -jg
2059c SECND_FRFI(NIN)%P(4:6,JG) is the old friction force
2060c if(itafi(nin)%p(jg)==31774.or.itafi(nin)%p(jg)==6992.or.itafi(nin)%p(jg)==7106)
2061c . print *,'slav_fr rem',itafi(nin)%p(jg),SECND_FRFI(NIN)%P(4:6,JG),PENE_OLDFI(nin)%p(1:5,JG),STIF_OLDFI(nin)%p(1:2,JG)
2062 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2063 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2064 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2065 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2066 fx = fx - ftn*n1(i)
2067 fy = fy - ftn*n2(i)
2068 fz = fz - ftn*n3(i)
2069
2070 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2071 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2072 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2073 ft = max(ft,em30)
2074 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2075
2076 beta = fn/ft
2077
2078 IF(beta == 0 ) THEN
2079 fxt(i) = zero
2080 fyt(i) = zero
2081 fzt(i) = zero
2082 ELSEIF(beta > 1) THEN ! Inside the ellipse
2083 fxt(i) = fx
2084 fyt(i) = fy
2085 fzt(i) = fz
2086 ELSE ! outside the ellipse
2087
2088! Projection on local tangent of ellipse (outside of ellipse)
2089! ANS = (Fadh-Fproj).n
2090 nep1 =ftt1*an(k)/fn
2091 nep2 =ftt2*bn(k)/fn
2092 nep =nep1*nep1+nep2*nep2
2093 nep =sqrt(nep)
2094
2095 ep=nep1*ftt1+nep2*ftt2
2096
2097 ans=(ep-sqrt(ep))/max(em20,nep)
2098 nep1 =nep1/max(em20,nep)
2099 nep2 =nep2/max(em20,nep)
2100
2101! Projection on ellipse
2102 c11 =ftt1-ans*nep1
2103 c22 =ftt2-ans*nep2
2104
2105 alphaf = atan(c22/c11)
2106
2107 signc = ftt1/max(em20,abs(ftt1))
2108 csa = signc*abs(cos(alphaf))
2109 signc = ftt2/max(em20,abs(ftt2))
2110 sna = signc*abs(sin(alphaf))
2111! Ft computation
2112 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2113 ftt1 = ft * csa
2114 ftt2 = ft * sna
2115
2116 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2117 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2118 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2119
2120 ENDIF
2121C
2122 secnd_frfi(nin)%P(1,jg) = fxt(i)
2123 secnd_frfi(nin)%P(2,jg) = fyt(i)
2124 secnd_frfi(nin)%P(3,jg) = fzt(i)
2125C
2126c if(itafi(nin)%p(jg)==27715)
2127c . print *,'i25for3 remote',itafi(nin)%p(jg),pene(i),fxt(i),fyt(i),fzt(i),fxi(i),fyi(i),fzi(i)
2128 ENDIF
2129C------- total force
2130 fxi(i)=fxi(i) + fxt(i)
2131 fyi(i)=fyi(i) + fyt(i)
2132 fzi(i)=fzi(i) + fzt(i)
2133
2134 IF( intth > 0 .AND.beta/=zero) THEN
2135 efrict(i) = (fx-fxt(i))*fxt(i) + (fy-fyt(i))*fyt(i) +
2136 . (fz-fzt(i))*fzt(i)
2137 efrict(i) = efrict(i)/stif0(i) ! FRICTIONAL ENERGY
2138 qfrict = qfrict + efrict(i)
2139 ENDIF
2140 efric_l(i) = dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
2141 econvt = econvt + efric_l(i)
2142
2143 ENDDO
2144C--------implicit non converge---
2145 ELSE
2146#include "vectorize.inc"
2147 DO k=1,nfisot ! isotropic friction couples
2148 i = indexisot(k)
2149 IF(pene(i) == zero)cycle
2150 fx = stif0(i)*vx(i)*dt12
2151 fy = stif0(i)*vy(i)*dt12
2152 fz = stif0(i)*vz(i)*dt12
2153 jg = nsvg(i)
2154 n = cand_n_n(i)
2155 IF(jg>0) THEN
2156c SECND_FR(1:3,N) is the old friction force
2157 fx = secnd_fr(4,n) + alpha*fx
2158 fy = secnd_fr(5,n) + alpha*fy
2159 fz = secnd_fr(6,n) + alpha*fz
2160 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2161 fx = fx - ftn*n1(i)
2162 fy = fy - ftn*n2(i)
2163 fz = fz - ftn*n3(i)
2164 ft = fx*fx + fy*fy + fz*fz
2165 ft = max(ft,tiny)
2166 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2167 beta = min(one,xmu(i)*sqrt(fn/ft))
2168 fxt(i) = fx * beta
2169 fyt(i) = fy * beta
2170 fzt(i) = fz * beta
2171 fxi(i)=fxi(i) + fxt(i)
2172 fyi(i)=fyi(i) + fyt(i)
2173 fzi(i)=fzi(i) + fzt(i)
2174 ELSE ! cas noeud remote en SPMD
2175 jg = -jg
2176 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2177 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2178 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2179 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2180 fx = fx - ftn*n1(i)
2181 fy = fy - ftn*n2(i)
2182 fz = fz - ftn*n3(i)
2183 ft = fx*fx + fy*fy + fz*fz
2184 ft = max(ft,tiny)
2185 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2186 beta = min(one,xmu(i)*sqrt(fn/ft))
2187 fxt(i) = fx * beta
2188 fyt(i) = fy * beta
2189 fzt(i) = fz * beta
2190 fxi(i)=fxi(i) + fxt(i)
2191 fyi(i)=fyi(i) + fyt(i)
2192 fzi(i)=fzi(i) + fzt(i)
2193 ENDIF
2194
2195 ENDDO
2196
2197#include "vectorize.inc"
2198 DO k=1,nforth ! Orthotropic friction couples
2199 i = indexorth(k)
2200 IF(pene(i) == zero)cycle
2201 fx = stif0(i)*vx(i)*dt12
2202 fy = stif0(i)*vy(i)*dt12
2203 fz = stif0(i)*vz(i)*dt12
2204 jg = nsvg(i)
2205 n = cand_n_n(i)
2206 IF(jg>0) THEN
2207c SECND_FR(1:3,N) is the old friction force
2208 fx = secnd_fr(4,n) + alpha*fx
2209 fy = secnd_fr(5,n) + alpha*fy
2210 fz = secnd_fr(6,n) + alpha*fz
2211 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2212 fx = fx - ftn*n1(i)
2213 fy = fy - ftn*n2(i)
2214 fz = fz - ftn*n3(i)
2215 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2216 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2217 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2218 ft = max(ft,em30)
2219 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2220
2221 beta = fn/ft
2222
2223 IF(beta == 0 ) THEN
2224 fxt(i) = zero
2225 fyt(i) = zero
2226 fzt(i) = zero
2227 ELSEIF(beta > 1) THEN ! Inside the ellipse
2228 fxt(i) = fx
2229 fyt(i) = fy
2230 fzt(i) = fz
2231 ELSE ! outside the ellipse
2232
2233! Projection on local tangent of ellipse (outside of ellipse)
2234! ANS = (Fadh-Fproj).n
2235 nep1 =ftt1*an(k)/fn
2236 nep2 =ftt2*bn(k)/fn
2237 nep =nep1*nep1+nep2*nep2
2238 nep =sqrt(nep)
2239
2240 ep=nep1*ftt1+nep2*ftt2
2241
2242 ans=(ep-sqrt(ep))/max(em20,nep)
2243 nep1 =nep1/max(em20,nep)
2244 nep2 =nep2/max(em20,nep)
2245
2246! Projection on ellipse
2247 c11 =ftt1-ans*nep1
2248 c22 =ftt2-ans*nep2
2249
2250 alphaf = atan(c22/c11)
2251
2252 signc = ftt1/max(em20,abs(ftt1))
2253 csa = signc*abs(cos(alphaf))
2254 signc = ftt2/max(em20,abs(ftt2))
2255 sna = signc*abs(sin(alphaf))
2256! Ft computation
2257 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2258 ftt1 = ft * csa
2259 ftt2 = ft * sna
2260
2261 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2262 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2263 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2264
2265 ENDIF
2266 fxi(i)=fxi(i) + fxt(i)
2267 fyi(i)=fyi(i) + fyt(i)
2268 fzi(i)=fzi(i) + fzt(i)
2269 ELSE ! cas noeud remote en SPMD
2270 jg = -jg
2271 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2272 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2273 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2274 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2275 fx = fx - ftn*n1(i)
2276 fy = fy - ftn*n2(i)
2277 fz = fz - ftn*n3(i)
2278 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2279 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2280 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2281 ft = max(ft,em30)
2282 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2283
2284 beta = fn/ft
2285
2286 IF(beta == 0 ) THEN
2287 fxt(i) = zero
2288 fyt(i) = zero
2289 fzt(i) = zero
2290 ELSEIF(beta > 1) THEN ! Inside the ellipse
2291 fxt(i) = fx
2292 fyt(i) = fy
2293 fzt(i) = fz
2294 ELSE ! outside the ellipse
2295
2296! Projection on local tangent of ellipse (outside of ellipse)
2297! ANS = (Fadh-Fproj).n
2298 nep1 =ftt1*an(k)/fn
2299 nep2 =ftt2*bn(k)/fn
2300 nep =nep1*nep1+nep2*nep2
2301 nep =sqrt(nep)
2302
2303 ep=nep1*ftt1+nep2*ftt2
2304
2305 ans=(ep-sqrt(ep))/max(em20,nep)
2306 nep1 =nep1/max(em20,nep)
2307 nep2 =nep2/max(em20,nep)
2308
2309! Projection on ellipse
2310 c11 =ftt1-ans*nep1
2311 c22 =ftt2-ans*nep2
2312
2313 alphaf = atan(c22/c11)
2314
2315 signc = ftt1/max(em20,abs(ftt1))
2316 csa = signc*abs(cos(alphaf))
2317 signc = ftt2/max(em20,abs(ftt2))
2318 sna = signc*abs(sin(alphaf))
2319! Ft computation
2320 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2321 ftt1 = ft * csa
2322 ftt2 = ft * sna
2323
2324 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2325 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2326 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2327
2328 ENDIF
2329 fxi(i)=fxi(i) + fxt(i)
2330 fyi(i)=fyi(i) + fyt(i)
2331 fzi(i)=fzi(i) + fzt(i)
2332 ENDIF
2333 ENDDO
2334 ENDIF
2335
2336
2337
2338 ENDIF
2339
2340 ENDIF
2341 ENDIF
2342C
2343C---------------------------------
2344 IF((anim_v(12)+outp_v(12)+h3d_data%N_VECT_PCONT>0.AND.
2345 . ((tt>=tanim .AND. tt<=tanim_stop).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
2346 . (manim>=4.AND.manim<=15).OR.h3d_data%MH3D/=0))
2347 . .OR.h3d_data%N_VECT_PCONT_MAX>0)THEN
2348 IF (inconv==1) THEN
2349#include "lockon.inc"
2350 DO i=1,jlt
2351 IF(pene(i) == zero)cycle
2352 ftcont(1,ix1(i)) =ftcont(1,ix1(i)) + fxt(i)*h1(i)
2353 ftcont(2,ix1(i)) =ftcont(2,ix1(i)) + fyt(i)*h1(i)
2354 ftcont(3,ix1(i)) =ftcont(3,ix1(i)) + fzt(i)*h1(i)
2355 ftcont(1,ix2(i)) =ftcont(1,ix2(i)) + fxt(i)*h2(i)
2356 ftcont(2,ix2(i)) =ftcont(2,ix2(i)) + fyt(i)*h2(i)
2357 ftcont(3,ix2(i)) =ftcont(3,ix2(i)) + fzt(i)*h2(i)
2358 ftcont(1,ix3(i)) =ftcont(1,ix3(i)) + fxt(i)*h3(i)
2359 ftcont(2,ix3(i)) =ftcont(2,ix3(i)) + fyt(i)*h3(i)
2360 ftcont(3,ix3(i)) =ftcont(3,ix3(i)) + fzt(i)*h3(i)
2361 ftcont(1,ix4(i)) =ftcont(1,ix4(i)) + fxt(i)*h4(i)
2362 ftcont(2,ix4(i)) =ftcont(2,ix4(i)) + fyt(i)*h4(i)
2363 ftcont(3,ix4(i)) =ftcont(3,ix4(i)) + fzt(i)*h4(i)
2364 jg = nsvg(i)
2365 IF(jg>0) THEN
2366C en SPMD : traitement a refaire apres reception noeud remote si JG < 0
2367 ftcont(1,jg)=ftcont(1,jg)- fxt(i)
2368 ftcont(2,jg)=ftcont(2,jg)- fyt(i)
2369 ftcont(3,jg)=ftcont(3,jg)- fzt(i)
2370 ELSE ! cas noeud remote en SPMD
2371 jg = -jg
2372 ftconti(nin)%P(1,jg)=ftconti(nin)%P(1,jg)-fxt(i)
2373 ftconti(nin)%P(2,jg)=ftconti(nin)%P(2,jg)-fyt(i)
2374 ftconti(nin)%P(3,jg)=ftconti(nin)%P(3,jg)-fzt(i)
2375 ENDIF
2376 ENDDO
2377#include "lockoff.inc"
2378 END IF !(INCONV==1) THEN
2379 ENDIF
2380C
2381C---------------------------------
2382 fsav4 = zero
2383 fsav5 = zero
2384 fsav6 = zero
2385 fsav12= zero
2386 fsav13= zero
2387 fsav14= zero
2388 fsav15= zero
2389 fsav22= zero
2390 fsav23= zero
2391 fsav24= zero
2392 fsav29= zero
2393 DO i=1,jlt
2394 IF(pene(i) == zero)cycle
2395 impx=fxt(i)*dt12
2396 impy=fyt(i)*dt12
2397 impz=fzt(i)*dt12
2398 fsav4 =fsav4 +impx
2399 fsav5 =fsav5 +impy
2400 fsav6 =fsav6 +impz
2401 impx=fxi(i)*dt12
2402 impy=fyi(i)*dt12
2403 impz=fzi(i)*dt12
2404 fsav12=fsav12+abs(impx)
2405 fsav13=fsav13+abs(impy)
2406 fsav14=fsav14+abs(impz)
2407 fsav15=fsav15+sqrt(impx*impx+impy*impy+impz*impz)
2408 xp(i)=xi(i)+pene(i)*n1(i)
2409 yp(i)=yi(i)+pene(i)*n2(i)
2410 zp(i)=zi(i)+pene(i)*n3(i)
2411 fsav22=fsav22+yp(i)*impz-zp(i)*impy
2412 fsav23=fsav23+zp(i)*impx-xp(i)*impz
2413 fsav24=fsav24+xp(i)*impy-yp(i)*impx
2414 ENDDO
2415
2416 IF(intcarea > 0) THEN
2417 DO i=1,jlt
2418 IF(pene(i) == zero)cycle
2419 jg = nsvg(i)
2420 IF(jg>0) THEN ! Only local nodes in i25for3, remote nodes will be taken into account later
2421 fsav29=fsav29+parameters%INTAREAN(jg)
2422 ENDIF
2423 ENDDO
2424 ENDIF
2425
2426#include "lockon.inc"
2427 fsav(4) = fsav(4) + fsav4
2428 fsav(5) = fsav(5) + fsav5
2429 fsav(6) = fsav(6) + fsav6
2430 fsav(12) = fsav(12) + fsav12
2431 fsav(13) = fsav(13) + fsav13
2432 fsav(14) = fsav(14) + fsav14
2433 fsav(15) = fsav(15) + fsav15
2434 fsav(22) = fsav(22) + fsav22
2435 fsav(23) = fsav(23) + fsav23
2436 fsav(24) = fsav(24) + fsav24
2437 fsav(25) = fsav(25) + (fheats+fheatm)*qfrict
2438 fsav(26) = fsav(26) + econtt
2439 fsav(27) = fsav(27) + econvt - (fheats+fheatm)*qfrict
2440 fsav(28) = fsav(28) + econtdt
2441 fsav(29) = fsav(29) + fsav29
2442#include "lockoff.inc"
2443C
2444 IF(isensint(1)/=0) THEN
2445 DO i=1,jlt
2446 IF(pene(i) == zero)cycle
2447 fsavparit(1,4,i) = fxt(i)
2448 fsavparit(1,5,i) = fyt(i)
2449 fsavparit(1,6,i) = fzt(i)
2450 ENDDO
2451 ENDIF
2452C
2453C---------------------------------
2454C SORTIES TH PAR SOUS INTERFACE
2455C---------------------------------
2456 IF(nisub/=0)THEN
2457 DO i=1,jlt
2458 IF(pene(i) == zero)cycle
2459 nn = nsvg(i)
2460 IF(nn>0)THEN
2461 in=cn_loc(i)
2462
2463 IF (msegtyp(ce_loc(i)) < 0) THEN
2464 ie= - msegtyp(ce_loc(i))
2465 ELSE
2466 ie = ce_loc(i)
2467 ENDIF
2468 IF(ie > nrtm) ie=ie-nrtm
2469
2470 jj =addsubs(in) ! Addresse du S
2471 kk =addsubm(ie) ! Addresse du M
2472
2473 DO WHILE(jj<addsubs(in+1))
2474! all sub interfaces of S
2475 jsub=lisubs(jj)
2476! JSUB Croissant
2477 itypsub = typsub(jsub)
2478 IF(itypsub == 1 ) THEN
2479C
2480C Find if node is on Surface S1 S2 ou GRNOD
2481 iss1 = bitget(inflg_subs(jj),0)
2482 iss2 = bitget(inflg_subs(jj),1)
2483 igrn = bitget(inflg_subs(jj),2)
2484 ksub=lisubm(kk)
2485 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
2486 ims1 = bitget(inflg_subm(kk),0)
2487 ims2 = bitget(inflg_subm(kk),1)
2488 IF(ksub==jsub)THEN
2489! S and M candidates on the same sub_interface
2490 IF(.NOT.(((ims1 == 1 .OR. ims2 == 1) .AND. igrn == 1).OR.
2491 . (ims1 == 1 .AND. iss2 == 1).OR.
2492 . (ims2 == 1 .AND. iss1 == 1))) THEN
2493 kk=kk+1
2494 ksub=lisubm(kk)
2495 cycle
2496 END IF
2497 impx=fxt(i)*dt12
2498 impy=fyt(i)*dt12
2499 impz=fzt(i)*dt12
2500C main side :
2501 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
2502 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
2503 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
2504C
2505 impx=fxi(i)*dt12
2506 impy=fyi(i)*dt12
2507 impz=fzi(i)*dt12
2508 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
2509 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
2510 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
2511C
2512 IF(isensint(jsub+1)/=0) THEN
2513 fsavparit(jsub+1,4,i) = fxt(i)
2514 fsavparit(jsub+1,5,i) = fyt(i)
2515 fsavparit(jsub+1,6,i) = fzt(i)
2516 ENDIF
2517C
2518 fsavsub1(15,jsub)= fsavsub1(15,jsub)
2519 . +sqrt(impx*impx+impy*impy+impz*impz)
2520 fsavsub1(22,jsub)=fsavsub1(22,jsub)
2521 . +yp(i)*impz-zp(i)*impy
2522 fsavsub1(23,jsub)=fsavsub1(23,jsub)
2523 . +zp(i)*impx-xp(i)*impz
2524 fsavsub1(24,jsub)=fsavsub1(24,jsub)
2525 . +xp(i)*impy-yp(i)*impx
2526C
2527 IF(parameters%INTCAREA > 0) fsavsub1(25,jsub) = fsavsub1(25,jsub) + parameters%INTAREAN(nn)
2528C
2529 ENDIF
2530 kk=kk+1
2531 ksub=lisubm(kk)
2532 ENDDO
2533 jj=jj+1
2534
2535 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
2536
2537 impx=fxt(i)*dt12
2538 impy=fyt(i)*dt12
2539 impz=fzt(i)*dt12
2540C main side :
2541 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
2542 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
2543 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
2544C
2545 impx=fxi(i)*dt12
2546 impy=fyi(i)*dt12
2547 impz=fzi(i)*dt12
2548 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
2549 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
2550 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
2551C
2552 IF(isensint(jsub+1)/=0) THEN
2553 fsavparit(jsub+1,4,i) = fxt(i)
2554 fsavparit(jsub+1,5,i) = fyt(i)
2555 fsavparit(jsub+1,6,i) = fzt(i)
2556 ENDIF
2557C
2558 fsavsub1(15,jsub)= fsavsub1(15,jsub)
2559 . +sqrt(impx*impx+impy*impy+impz*impz)
2560 fsavsub1(22,jsub)=fsavsub1(22,jsub)
2561 . +yp(i)*impz-zp(i)*impy
2562 fsavsub1(23,jsub)=fsavsub1(23,jsub)
2563 . +zp(i)*impx-xp(i)*impz
2564 fsavsub1(24,jsub)=fsavsub1(24,jsub)
2565 . +xp(i)*impy-yp(i)*impx
2566C
2567 IF(parameters%INTCAREA > 0) fsavsub1(25,jsub) = fsavsub1(25,jsub) + parameters%INTAREAN(nn)
2568C
2569 jj=jj+1
2570
2571 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfaces
2572
2573C
2574C Find if node is on Surface S1 S2 ou GRNOD
2575 iss2 = bitget(inflg_subs(jj),0)
2576 iss1 = bitget(inflg_subs(jj),1)
2577 ksub=lisubm(kk)
2578 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
2579 ims2 = bitget(inflg_subm(kk),0)
2580 ims1 = bitget(inflg_subm(kk),1)
2581 IF(ksub==jsub)THEN
2582! S and M candidates on the same sub_interface
2583 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
2584 . (ims2 == 1 .AND. iss1 == 1))) THEN
2585 kk=kk+1
2586 ksub=lisubm(kk)
2587 cycle
2588 END IF
2589
2590 impx=fxt(i)*dt12
2591 impy=fyt(i)*dt12
2592 impz=fzt(i)*dt12
2593
2594 IF(ims2 > 0)THEN
2595C main side :
2596 fsavsub1(4,jsub)=fsavsub1(4,jsub)-impx
2597 fsavsub1(5,jsub)=fsavsub1(5,jsub)-impy
2598 fsavsub1(6,jsub)=fsavsub1(6,jsub)-impz
2599 ELSE
2600 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
2601 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
2602 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
2603 ENDIF
2604C
2605 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
2606 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
2607 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
2608C
2609 IF(isensint(jsub+1)/=0) THEN
2610 IF(ims2 > 0)THEN
2611 fsavparit(jsub+1,4,i) = -fxt(i)
2612 fsavparit(jsub+1,5,i) = -fyt(i)
2613 fsavparit(jsub+1,6,i) = -fzt(i)
2614 ELSE
2615 fsavparit(jsub+1,4,i) = fxt(i)
2616 fsavparit(jsub+1,5,i) = fyt(i)
2617 fsavparit(jsub+1,6,i) = fzt(i)
2618 ENDIF
2619 ENDIF
2620C
2621 fsavsub1(15,jsub)= fsavsub1(15,jsub)
2622 . +sqrt(impx*impx+impy*impy+impz*impz)
2623 fsavsub1(22,jsub)=fsavsub1(22,jsub)
2624 . +yp(i)*impz-zp(i)*impy
2625 fsavsub1(23,jsub)=fsavsub1(23,jsub)
2626 . +zp(i)*impx-xp(i)*impz
2627 fsavsub1(24,jsub)=fsavsub1(24,jsub)
2628 . +xp(i)*impy-yp(i)*impx
2629C
2630 IF(parameters%INTCAREA > 0) fsavsub1(25,jsub) = fsavsub1(25,jsub) + parameters%INTAREAN(nn)
2631C
2632 ENDIF
2633 kk=kk+1
2634 ksub=lisubm(kk)
2635 ENDDO
2636 jj=jj+1
2637
2638 ENDIF
2639 END DO
2640 END IF
2641
2642
2643 IF (msegtyp(ce_loc(i)) < 0) THEN
2644 ie= - msegtyp(ce_loc(i))
2645 ELSE
2646 ie = ce_loc(i)
2647 ENDIF
2648 IF(ie > nrtm) ie=ie-nrtm
2649
2650 kk =addsubm(ie) ! Addresse du M
2651 DO WHILE(kk<addsubm(ie+1))
2652! all sub interfaces of S
2653 ksub=lisubm(kk)
2654! KSUB Croissant
2655 itypsub = typsub(ksub)
2656 IF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
2657 impx=-fxt(i)*dt12
2658 impy=-fyt(i)*dt12
2659 impz=-fzt(i)*dt12
2660C main side :
2661 fsavsub1(4,ksub)=fsavsub1(4,ksub)+impx
2662 fsavsub1(5,ksub)=fsavsub1(5,ksub)+impy
2663 fsavsub1(6,ksub)=fsavsub1(6,ksub)+impz
2664C
2665 impx=fxi(i)*dt12
2666 impy=fyi(i)*dt12
2667 impz=fzi(i)*dt12
2668 fsavsub1(12,ksub)=fsavsub1(12,ksub)+abs(impx)
2669 fsavsub1(13,ksub)=fsavsub1(13,ksub)+abs(impy)
2670 fsavsub1(14,ksub)=fsavsub1(14,ksub)+abs(impz)
2671C
2672 IF(isensint(ksub+1)/=0) THEN
2673 fsavparit(ksub+1,4,i) = -fxt(i)
2674 fsavparit(ksub+1,5,i) = -fyt(i)
2675 fsavparit(ksub+1,6,i) = -fzt(i)
2676 ENDIF
2677C
2678 fsavsub1(15,ksub)= fsavsub1(15,ksub)
2679 . +sqrt(impx*impx+impy*impy+impz*impz)
2680 fsavsub1(22,ksub)=fsavsub1(22,ksub)
2681 . +yp(i)*impz-zp(i)*impy
2682 fsavsub1(23,ksub)=fsavsub1(23,ksub)
2683 . +zp(i)*impx-xp(i)*impz
2684 fsavsub1(24,ksub)=fsavsub1(24,ksub)
2685 . +xp(i)*impy-yp(i)*impx
2686
2687 IF(parameters%INTCAREA > 0) THEN
2688 IF(nn > 0) THEN
2689 fsavsub1(25,ksub) = fsavsub1(25,ksub) + parameters%INTAREAN(nn)
2690 ELSE
2691 fsavsub1(25,ksub) = fsavsub1(25,ksub) + intareanfi(nin)%P(-nn)
2692 ENDIF
2693 ENDIF
2694
2695 ENDIF
2696 kk=kk+1
2697 ENDDO
2698
2699
2700 END DO
2701C
2702 IF(nspmd>1) THEN
2703 DO i=1,jlt
2704 IF(pene(i) == zero)cycle
2705 nn = nsvg(i)
2706 IF(nn<0)THEN
2707 nn = -nn
2708
2709 IF (msegtyp(ce_loc(i)) < 0) THEN
2710 ie= - msegtyp(ce_loc(i))
2711 ELSE
2712 ie = ce_loc(i)
2713 ENDIF
2714 IF(ie > nrtm) ie=ie-nrtm
2715
2716 jj =addsubsfi(nin)%P(nn)
2717 kk =addsubm(ie)
2718 DO WHILE(jj<addsubsfi(nin)%P(nn+1))
2719 jsub=lisubsfi(nin)%P(jj)
2720 itypsub = typsub(jsub)
2721 IF(itypsub == 1 ) THEN
2722 iss1 = bitget(inflg_subsfi(nin)%P(jj),0)
2723 iss2 = bitget(inflg_subsfi(nin)%P(jj),1)
2724 igrn = bitget(inflg_subsfi(nin)%P(jj),2)
2725 ksub=lisubm(kk)
2726 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
2727 ims1 = bitget(inflg_subm(kk),0)
2728 ims2 = bitget(inflg_subm(kk),1)
2729 IF(ksub==jsub)THEN
2730 IF(.NOT.(((ims1 == 1 .OR. ims2 == 1) .AND. igrn == 1).OR.
2731 . (ims1 == 1 .AND. iss2 == 1).OR.
2732 . (ims2 == 1 .AND. iss1 == 1))) THEN
2733 kk=kk+1
2734 ksub=lisubm(kk)
2735 cycle
2736 END IF
2737 impx=fxt(i)*dt12
2738 impy=fyt(i)*dt12
2739 impz=fzt(i)*dt12
2740C main side :
2741 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
2742 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
2743 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
2744C
2745 impx=fxi(i)*dt12
2746 impy=fyi(i)*dt12
2747 impz=fzi(i)*dt12
2748 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
2749 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
2750 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
2751C
2752 IF(isensint(jsub+1)/=0) THEN
2753 fsavparit(jsub+1,4,i) = fxt(i)
2754 fsavparit(jsub+1,5,i) = fyt(i)
2755 fsavparit(jsub+1,6,i) = fzt(i)
2756 ENDIF
2757C
2758 fsavsub1(15,jsub)= fsavsub1(15,jsub)
2759 . +sqrt(impx*impx+impy*impy+impz*impz)
2760 fsavsub1(22,jsub)=fsavsub1(22,jsub)
2761 . +yp(i)*impz-zp(i)*impy
2762 fsavsub1(23,jsub)=fsavsub1(23,jsub)
2763 . +zp(i)*impx-xp(i)*impz
2764 fsavsub1(24,jsub)=fsavsub1(24,jsub)
2765 . +xp(i)*impy-yp(i)*impx
2766
2767 IF(parameters%INTCAREA > 0) fsavsub1(25,jsub) = fsavsub1(25,jsub) + intareanfi(nin)%P(nn)
2768
2769 ENDIF
2770 kk=kk+1
2771 ksub=lisubm(kk)
2772 ENDDO
2773 jj=jj+1
2774
2775 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
2776
2777 impx=fxt(i)*dt12
2778 impy=fyt(i)*dt12
2779 impz=fzt(i)*dt12
2780C main side :
2781 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
2782 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
2783 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
2784C
2785 impx=fxi(i)*dt12
2786 impy=fyi(i)*dt12
2787 impz=fzi(i)*dt12
2788 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
2789 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
2790 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
2791C
2792 IF(isensint(jsub+1)/=0) THEN
2793 fsavparit(jsub+1,4,i) = fxt(i)
2794 fsavparit(jsub+1,5,i) = fyt(i)
2795 fsavparit(jsub+1,6,i) = fzt(i)
2796 ENDIF
2797C
2798 fsavsub1(15,jsub)= fsavsub1(15,jsub)
2799 . +sqrt(impx*impx+impy*impy+impz*impz)
2800 fsavsub1(22,jsub)=fsavsub1(22,jsub)
2801 . +yp(i)*impz-zp(i)*impy
2802 fsavsub1(23,jsub)=fsavsub1(23,jsub)
2803 . +zp(i)*impx-xp(i)*impz
2804 fsavsub1(24,jsub)=fsavsub1(24,jsub)
2805 . +xp(i)*impy-yp(i)*impx
2806
2807 IF(parameters%INTCAREA > 0) fsavsub1(25,jsub) = fsavsub1(25,jsub) + intareanfi(nin)%P(nn)
2808
2809 jj=jj+1
2810
2811 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfaces
2812
2813 iss2 = bitget(inflg_subsfi(nin)%P(jj),0)
2814 iss1 = bitget(inflg_subsfi(nin)%P(jj),1)
2815 igrn = bitget(inflg_subsfi(nin)%P(jj),2)
2816 ksub=lisubm(kk)
2817 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
2818 ims2 = bitget(inflg_subm(kk),0)
2819 ims1 = bitget(inflg_subm(kk),1)
2820 IF(ksub==jsub)THEN
2821 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
2822 . (ims2 == 1 .AND. iss1 == 1))) THEN
2823 kk=kk+1
2824 ksub=lisubm(kk)
2825 cycle
2826 END IF
2827
2828 impx=fxt(i)*dt12
2829 impy=fyt(i)*dt12
2830 impz=fzt(i)*dt12
2831 IF(ims2 > 0)THEN
2832C main side :
2833 fsavsub1(4,jsub)=fsavsub1(4,jsub)-impx
2834 fsavsub1(5,jsub)=fsavsub1(5,jsub)-impy
2835 fsavsub1(6,jsub)=fsavsub1(6,jsub)-impz
2836 ELSE
2837 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
2838 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
2839 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
2840 ENDIF
2841C
2842 impx=fxi(i)*dt12
2843 impy=fyi(i)*dt12
2844 impz=fzi(i)*dt12
2845
2846 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
2847 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
2848 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
2849C
2850 IF(isensint(jsub+1)/=0) THEN
2851 IF(ims2 > 0)THEN
2852 fsavparit(jsub+1,4,i) = -fxt(i)
2853 fsavparit(jsub+1,5,i) = -fyt(i)
2854 fsavparit(jsub+1,6,i) = -fzt(i)
2855 ELSE
2856 fsavparit(jsub+1,4,i) = fxt(i)
2857 fsavparit(jsub+1,5,i) = fyt(i)
2858 fsavparit(jsub+1,6,i) = fzt(i)
2859 ENDIF
2860 ENDIF
2861C
2862 fsavsub1(15,jsub)= fsavsub1(15,jsub)
2863 . +sqrt(impx*impx+impy*impy+impz*impz)
2864 fsavsub1(22,jsub)=fsavsub1(22,jsub)
2865 . +yp(i)*impz-zp(i)*impy
2866 fsavsub1(23,jsub)=fsavsub1(23,jsub)
2867 . +zp(i)*impx-xp(i)*impz
2868 fsavsub1(24,jsub)=fsavsub1(24,jsub)
2869 . +xp(i)*impy-yp(i)*impx
2870
2871 IF(parameters%INTCAREA > 0) fsavsub1(25,jsub) = fsavsub1(25,jsub) + intareanfi(nin)%P(nn)
2872 ENDIF
2873 kk=kk+1
2874 ksub=lisubm(kk)
2875 ENDDO
2876 jj=jj+1
2877
2878 ENDIF
2879
2880 END DO
2881
2882 END IF
2883 END DO
2884 END IF
2885#include "lockon.inc"
2886 DO jsub=1,nisub
2887 nsub=lisub(jsub)
2888 DO j=1,15
2889 fsavsub(j,nsub)=fsavsub(j,nsub)+fsavsub1(j,jsub)
2890 END DO
2891 fsavsub(22,nsub)=fsavsub(22,nsub)+fsavsub1(22,jsub)
2892 fsavsub(23,nsub)=fsavsub(23,nsub)+fsavsub1(23,jsub)
2893 fsavsub(24,nsub)=fsavsub(24,nsub)+fsavsub1(24,jsub)
2894 IF(parameters%INTCAREA > 0) fsavsub(29,nsub)=fsavsub(29,nsub)+fsavsub1(25,jsub)
2895 END DO
2896#include "lockoff.inc"
2897 END IF
2898C---------------------------------
2899#include "lockon.inc"
2900 IF (inconv==1) THEN
2901 econtv = econtv + econvt ! Frictional Energy
2902 econt = econt + econtt ! Elastic Energy
2903 econtd = econtd + econtdt ! Damping Energy
2904 IF (intth > 0) THEN
2905 qfric = qfric + (fheats+fheatm)*qfrict ! FRICTIONAL HEAT ADDED TO INTERNAL ENERGY
2906 econtv = econtv - (fheats+fheatm)*qfrict ! FRICTIONAL HEAT REMOVED FROM CONTACT ENERGY
2907 ENDIF
2908 END IF !(INCONV==1) THEN
2909#include "lockoff.inc"
2910C---------------------------------
2911C---------------------------------
2912 IF(interefric >0)THEN
2913 IF (inconv==1) THEN
2914#include "lockon.inc"
2915 DO i=1,jlt
2916 IF(pene(i) == zero)cycle
2917 efric_lm = half*efric_l(i)
2918 efric(interefric,ix1(i)) = efric(interefric,ix1(i)) + h1(i)*(efric_lm-fheatm*efrict(i))
2919 efric(interefric,ix2(i)) = efric(interefric,ix2(i)) + h2(i)*(efric_lm-fheatm*efrict(i))
2920 efric(interefric,ix3(i)) = efric(interefric,ix3(i)) + h3(i)*(efric_lm-fheatm*efrict(i))
2921 efric(interefric,ix4(i)) = efric(interefric,ix4(i)) + h4(i)*(efric_lm-fheatm*efrict(i))
2922 jg = nsvg(i)
2923 efric_ls = half*efric_l(i)
2924 IF(jg>0) THEN
2925 efric(interefric,jg) = efric(interefric,jg) + (efric_ls-fheats*efrict(i))
2926 ELSE ! node remote
2927 jg = -jg
2928 efricfi(nin)%P(jg)=efricfi(nin)%P(jg)+ (efric_ls-fheats*efrict(i))
2929 ENDIF
2930 ENDDO
2931#include "lockoff.inc"
2932 END IF !(INCONV==1) THEN
2933 ENDIF
2934C---------------------------------
2935 IF(h3d_data%N_SCAL_CSE_FRIC >0)THEN
2936 IF (inconv==1) THEN
2937#include "lockon.inc"
2938 DO i=1,jlt
2939 IF(pene(i) == zero)cycle
2940 efric_lm = half*efric_l(i)
2941 efricg(ix1(i)) = efricg(ix1(i)) + h1(i)*(efric_lm-fheatm*efrict(i))
2942 efricg(ix2(i)) = efricg(ix2(i)) + h2(i)*(efric_lm-fheatm*efrict(i))
2943 efricg(ix3(i)) = efricg(ix3(i)) + h3(i)*(efric_lm-fheatm*efrict(i))
2944 efricg(ix4(i)) = efricg(ix4(i)) + h4(i)*(efric_lm-fheatm*efrict(i))
2945 jg = nsvg(i)
2946 efric_ls = half*efric_l(i)
2947 IF(jg>0) THEN
2948 efricg(jg) = efricg(jg) + (efric_ls-fheats*efrict(i))
2949 ELSE ! node remote
2950 jg = -jg
2951 efricgfi(nin)%P(jg)=efricgfi(nin)%P(jg)+ (efric_ls-fheats*efrict(i))
2952 ENDIF
2953 ENDDO
2954#include "lockoff.inc"
2955 END IF !(INCONV==1) THEN
2956 ENDIF
2957C---------------------------------
2958 IF(kdtint==1)THEN
2959 IF( (visc/=zero)
2960 . .AND.(ivis2==6.OR.ivis2==1))THEN
2961 DO i=1,jlt
2962 IF(pene(i) == zero)cycle
2963C C(I)=2.*C(I)
2964C
2965 IF(msi(i)==zero)THEN
2966 ks(i) =zero
2967 cs(i) =zero
2968 stv(i)=zero
2969 ELSE
2970 cx = four*c(i)*c(i)
2971 cy = eight*msi(i)*kt(i)
2972 aux = sqrt(cx+cy)+two*c(i)
2973 stv(i)= kt(i)*aux*aux/max(cy,em30)
2974 aux = two*cf(i)*cf(i)/max(msi(i),em20)
2975 IF(aux>stv(i))THEN
2976 ks(i) =zero
2977 cs(i) =cf(i)
2978 stv(i)=aux
2979 ELSE
2980 ks(i)= kt(i)
2981 cs(i) =c(i)
2982 ENDIF
2983 ENDIF
2984C
2985 j1=ix1(i)
2986 IF(ms(j1)==zero)THEN
2987 k1(i) =zero
2988 c1(i) =zero
2989 st1(i)=zero
2990 ELSE
2991 k1(i)=kt(i)*abs(h1(i))
2992 c1(i)=c(i)*abs(h1(i))
2993 cx =four*c1(i)*c1(i)
2994 cy =eight*ms(j1)*k1(i)
2995 aux = sqrt(cx+cy)+two*c1(i)
2996 st1(i)= k1(i)*aux*aux/max(cy,em30)
2997 cfi = cf(i)*abs(h1(i))
2998 aux = two*cfi*cfi/max(ms(j1),em20)
2999 IF(aux>st1(i))THEN
3000 k1(i) =zero
3001 c1(i) =cfi
3002 st1(i)=aux
3003 ENDIF
3004 ENDIF
3005C
3006 j1=ix2(i)
3007 IF(ms(j1)==zero)THEN
3008 k2(i) =zero
3009 c2(i) =zero
3010 st2(i)=zero
3011 ELSE
3012 k2(i)=kt(i)*abs(h2(i))
3013 c2(i)=c(i)*abs(h2(i))
3014 cx =four*c2(i)*c2(i)
3015 cy =eight*ms(j1)*k2(i)
3016 aux = sqrt(cx+cy)+two*c2(i)
3017 st2(i)= k2(i)*aux*aux/max(cy,em30)
3018 cfi = cf(i)*abs(h2(i))
3019 aux = two*cfi*cfi/max(ms(j1),em20)
3020 IF(aux>st2(i))THEN
3021 k2(i) =zero
3022 c2(i) =cfi
3023 st2(i)=aux
3024 ENDIF
3025 ENDIF
3026C
3027 j1=ix3(i)
3028 IF(ms(j1)==zero)THEN
3029 k3(i) =zero
3030 c3(i) =zero
3031 st3(i)=zero
3032 ELSE
3033 k3(i)=kt(i)*abs(h3(i))
3034 c3(i)=c(i)*abs(h3(i))
3035 cx =four*c3(i)*c3(i)
3036 cy =eight*ms(j1)*k3(i)
3037 aux = sqrt(cx+cy)+two*c3(i)
3038 st3(i)= k3(i)*aux*aux/max(cy,em30)
3039 cfi = cf(i)*abs(h3(i))
3040 aux = two*cfi*cfi/max(ms(j1),em20)
3041 IF(aux>st3(i))THEN
3042 k3(i) =zero
3043 c3(i) =cfi
3044 st3(i)=aux
3045 ENDIF
3046 ENDIF
3047C
3048 j1=ix4(i)
3049 IF(ms(j1)==zero)THEN
3050 k4(i) =zero
3051 c4(i) =zero
3052 st4(i)=zero
3053 ELSE
3054 k4(i)=kt(i)*abs(h4(i))
3055 c4(i)=c(i)*abs(h4(i))
3056 cx =four*c4(i)*c4(i)
3057 cy =eight*ms(j1)*k4(i)
3058 aux = sqrt(cx+cy)+two*c4(i)
3059 st4(i)= k4(i)*aux*aux/max(cy,em30)
3060 cfi = cf(i)*abs(h4(i))
3061 aux = two*cfi*cfi/max(ms(j1),em20)
3062 IF(aux>st4(i))THEN
3063 k4(i) =zero
3064 c4(i) =cfi
3065 st4(i)=aux
3066 ENDIF
3067 ENDIF
3068 ENDDO
3069C
3070 ELSE
3071 DO i=1,jlt
3072 IF(viscffric(i)/=zero) THEN
3073 IF(pene(i) == zero)cycle
3074C C(I)=2.*C(I)
3075C
3076 IF(msi(i)==zero)THEN
3077 ks(i) =zero
3078 cs(i) =zero
3079 stv(i)=zero
3080 ELSE
3081 cx = four*c(i)*c(i)
3082 cy = eight*msi(i)*kt(i)
3083 aux = sqrt(cx+cy)+two*c(i)
3084 stv(i)= kt(i)*aux*aux/max(cy,em30)
3085 aux = two*cf(i)*cf(i)/max(msi(i),em20)
3086 IF(aux>stv(i))THEN
3087 ks(i) =zero
3088 cs(i) =cf(i)
3089 stv(i)=aux
3090 ELSE
3091 ks(i)= kt(i)
3092 cs(i) =c(i)
3093 ENDIF
3094 ENDIF
3095C
3096 j1=ix1(i)
3097 IF(ms(j1)==zero)THEN
3098 k1(i) =zero
3099 c1(i) =zero
3100 st1(i)=zero
3101 ELSE
3102 k1(i)=kt(i)*abs(h1(i))
3103 c1(i)=c(i)*abs(h1(i))
3104 cx =four*c1(i)*c1(i)
3105 cy =eight*ms(j1)*k1(i)
3106 aux = sqrt(cx+cy)+two*c1(i)
3107 st1(i)= k1(i)*aux*aux/max(cy,em30)
3108 cfi = cf(i)*abs(h1(i))
3109 aux = two*cfi*cfi/max(ms(j1),em20)
3110 IF(aux>st1(i))THEN
3111 k1(i) =zero
3112 c1(i) =cfi
3113 st1(i)=aux
3114 ENDIF
3115 ENDIF
3116C
3117 j1=ix2(i)
3118 IF(ms(j1)==zero)THEN
3119 k2(i) =zero
3120 c2(i) =zero
3121 st2(i)=zero
3122 ELSE
3123 k2(i)=kt(i)*abs(h2(i))
3124 c2(i)=c(i)*abs(h2(i))
3125 cx =four*c2(i)*c2(i)
3126 cy =eight*ms(j1)*k2(i)
3127 aux = sqrt(cx+cy)+two*c2(i)
3128 st2(i)= k2(i)*aux*aux/max(cy,em30)
3129 cfi = cf(i)*abs(h2(i))
3130 aux = two*cfi*cfi/max(ms(j1),em20)
3131 IF(aux>st2(i))THEN
3132 k2(i) =zero
3133 c2(i) =cfi
3134 st2(i)=aux
3135 ENDIF
3136 ENDIF
3137C
3138 j1=ix3(i)
3139 IF(ms(j1)==zero)THEN
3140 k3(i) =zero
3141 c3(i) =zero
3142 st3(i)=zero
3143 ELSE
3144 k3(i)=kt(i)*abs(h3(i))
3145 c3(i)=c(i)*abs(h3(i))
3146 cx =four*c3(i)*c3(i)
3147 cy =eight*ms(j1)*k3(i)
3148 aux = sqrt(cx+cy)+two*c3(i)
3149 st3(i)= k3(i)*aux*aux/max(cy,em30)
3150 cfi = cf(i)*abs(h3(i))
3151 aux = two*cfi*cfi/max(ms(j1),em20)
3152 IF(aux>st3(i))THEN
3153 k3(i) =zero
3154 c3(i) =cfi
3155 st3(i)=aux
3156 ENDIF
3157 ENDIF
3158C
3159 j1=ix4(i)
3160 IF(ms(j1)==zero)THEN
3161 k4(i) =zero
3162 c4(i) =zero
3163 st4(i)=zero
3164 ELSE
3165 k4(i)=kt(i)*abs(h4(i))
3166 c4(i)=c(i)*abs(h4(i))
3167 cx =four*c4(i)*c4(i)
3168 cy =eight*ms(j1)*k4(i)
3169 aux = sqrt(cx+cy)+two*c4(i)
3170 st4(i)= k4(i)*aux*aux/max(cy,em30)
3171 cfi = cf(i)*abs(h4(i))
3172 aux = two*cfi*cfi/max(ms(j1),em20)
3173 IF(aux>st4(i))THEN
3174 k4(i) =zero
3175 c4(i) =cfi
3176 st4(i)=aux
3177 ENDIF
3178 ENDIF
3179
3180
3181 ELSE
3182 IF(pene(i) == zero)cycle
3183 ks(i) =stif(i)
3184 cs(i) =zero
3185 stv(i)=ks(i)
3186 k1(i) =stif(i)*abs(h1(i))
3187 c1(i) =zero
3188 st1(i)=k1(i)
3189 k2(i) =stif(i)*abs(h2(i))
3190 c2(i) =zero
3191 st2(i)=k2(i)
3192 k3(i) =stif(i)*abs(h3(i))
3193 c3(i) =zero
3194 st3(i)=k3(i)
3195 k4(i) =stif(i)*abs(h4(i))
3196 c4(i) =zero
3197 st4(i)=k4(i)
3198 ENDIF
3199 ENDDO
3200 ENDIF
3201 ENDIF
3202C
3203C=======================================================================
3204C---------------------------------
3205 itag = 0
3206 DO i=1,jlt
3207 IF(pene(i) == zero)cycle
3208!!
3209 fx1(i)=fxi(i)*h1(i)
3210 fy1(i)=fyi(i)*h1(i)
3211 fz1(i)=fzi(i)*h1(i)
3212C
3213 fx2(i)=fxi(i)*h2(i)
3214 fy2(i)=fyi(i)*h2(i)
3215 fz2(i)=fzi(i)*h2(i)
3216C
3217 fx3(i)=fxi(i)*h3(i)
3218 fy3(i)=fyi(i)*h3(i)
3219 fz3(i)=fzi(i)*h3(i)
3220C
3221 fx4(i)=fxi(i)*h4(i)
3222 fy4(i)=fyi(i)*h4(i)
3223 fz4(i)=fzi(i)*h4(i)
3224C
3225 ENDDO
3226C
3227C
3228C spmd : identification des noeuds interf. utiles a envoyer
3229 IF (nspmd>1) THEN
3230ctmp+1 mic only
3231#include "mic_lockon.inc"
3232 DO i = 1,jlt
3233 nn = nsvg(i)
3234 IF(nn<0 .AND. h1(i)+h2(i)+h3(i)+h4(i)/=0)THEN
3235C tag temporaire de NSVFI a -
3236 nsvfi(nin)%P(-nn) = -abs(nsvfi(nin)%P(-nn))
3237 ENDIF
3238 ENDDO
3239ctmp+1 mic only
3240#include "mic_lockoff.inc"
3241 ENDIF
3242C
3243 IF(idtmins==2.OR.idtmins_int/=0)THEN
3244 dti=dt2t
3245 CALL i25sms2(jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
3246 2 nsvg ,h1 ,h2 ,h3 ,h4 ,stif ,
3247 3 nin ,noint ,mskyi_sms, iskyi_sms,nsms ,
3248 4 kt ,c ,cf ,dtmini,dti )
3249 IF(dti<dt2t)THEN
3250 dt2t = dti
3251 neltst = noint
3252 ityptst = 10
3253 ENDIF
3254 ENDIF
3255C
3256 IF(idtmins_int/=0)THEN
3257 stif(1:jlt)=zero
3258 END IF
3259C
3260
3261
3262
3263C
3264 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT>0)THEN
3265 IF (inconv==1) THEN
3266#include "lockon.inc"
3267 DO i=1,jlt
3268 IF(pene(i) == zero)cycle
3269 fcont(1,ix1(i)) =fcont(1,ix1(i)) + fx1(i)
3270 fcont(2,ix1(i)) =fcont(2,ix1(i)) + fy1(i)
3271 fcont(3,ix1(i)) =fcont(3,ix1(i)) + fz1(i)
3272 fcont(1,ix2(i)) =fcont(1,ix2(i)) + fx2(i)
3273 fcont(2,ix2(i)) =fcont(2,ix2(i)) + fy2(i)
3274 fcont(3,ix2(i)) =fcont(3,ix2(i)) + fz2(i)
3275 fcont(1,ix3(i)) =fcont(1,ix3(i)) + fx3(i)
3276 fcont(2,ix3(i)) =fcont(2,ix3(i)) + fy3(i)
3277 fcont(3,ix3(i)) =fcont(3,ix3(i)) + fz3(i)
3278 fcont(1,ix4(i)) =fcont(1,ix4(i)) + fx4(i)
3279 fcont(2,ix4(i)) =fcont(2,ix4(i)) + fy4(i)
3280 fcont(3,ix4(i)) =fcont(3,ix4(i)) + fz4(i)
3281 jg = nsvg(i)
3282 IF(jg>0) THEN
3283C en SPMD : traitement a refaire apres reception noeud remote si JG < 0
3284 fcont(1,jg)=fcont(1,jg)- fxi(i)
3285 fcont(2,jg)=fcont(2,jg)- fyi(i)
3286 fcont(3,jg)=fcont(3,jg)- fzi(i)
3287 ENDIF
3288 ENDDO
3289#include "lockoff.inc"
3290 END IF !(INCONV==1) THEN
3291 ENDIF
3292C-----------------------------------------------------
3293 IF(isecin>0.AND.inconv==1)THEN
3294 k0=nstrf(25)
3295 IF(nstrf(1)+nstrf(2)/=0)THEN
3296 DO i=1,nsect
3297 nbinter=nstrf(k0+14)
3298 k1s=k0+30
3299 DO j=1,nbinter
3300 IF(nstrf(k1s)==noint)THEN
3301 IF(isecut/=0)THEN
3302#include "lockon.inc"
3303 DO k=1,jlt
3304 IF(pene(k) == zero)cycle
3305C attention aux signes pour le cumul des efforts
3306C a rendre conforme avec CFORC3
3307 IF(secfcum(4,ix1(k),i)==1.)THEN
3308 secfcum(1,ix1(k),i)=secfcum(1,ix1(k),i)-fx1(k)
3309 secfcum(2,ix1(k),i)=secfcum(2,ix1(k),i)-fy1(k)
3310 secfcum(3,ix1(k),i)=secfcum(3,ix1(k),i)-fz1(k)
3311 ENDIF
3312 IF(secfcum(4,ix2(k),i)==1.)THEN
3313 secfcum(1,ix2(k),i)=secfcum(1,ix2(k),i)-fx2(k)
3314 secfcum(2,ix2(k),i)=secfcum(2,ix2(k),i)-fy2(k)
3315 secfcum(3,ix2(k),i)=secfcum(3,ix2(k),i)-fz2(k)
3316 ENDIF
3317 IF(secfcum(4,ix3(k),i)==1.)THEN
3318 secfcum(1,ix3(k),i)=secfcum(1,ix3(k),i)-fx3(k)
3319 secfcum(2,ix3(k),i)=secfcum(2,ix3(k),i)-fy3(k)
3320 secfcum(3,ix3(k),i)=secfcum(3,ix3(k),i)-fz3(k)
3321 ENDIF
3322 IF(secfcum(4,ix4(k),i)==1.)THEN
3323 secfcum(1,ix4(k),i)=secfcum(1,ix4(k),i)-fx4(k)
3324 secfcum(2,ix4(k),i)=secfcum(2,ix4(k),i)-fy4(k)
3325 secfcum(3,ix4(k),i)=secfcum(3,ix4(k),i)-fz4(k)
3326 ENDIF
3327 jg = nsvg(k)
3328 IF(jg>0) THEN
3329C en SPMD : traitement a refaire apres reception noeud remote si JG < 0
3330 IF(secfcum(4,jg,i)==1.)THEN
3331 secfcum(1,jg,i)=secfcum(1,jg,i)+fxi(k)
3332 secfcum(2,jg,i)=secfcum(2,jg,i)+fyi(k)
3333 secfcum(3,jg,i)=secfcum(3,jg,i)+fzi(k)
3334 ENDIF
3335 ENDIF
3336 ENDDO
3337#include "lockoff.inc"
3338 ENDIF
3339C +fsav(section)
3340 ENDIF
3341 k1s=k1s+1
3342 ENDDO
3343 k0=nstrf(k0+24)
3344 ENDDO
3345 ENDIF
3346 ENDIF
3347C-----------------------------------------------------
3348 IF(ibag/=0.OR.iadm/=0)THEN
3349 DO i=1,jlt
3350 IF(pene(i) == zero)cycle
3351C test modifie pour coherence avec communication SPMD (spmd_i7tools)
3352c IF(FXI(I)/=ZERO.OR.FYI(I)/=ZERO.OR.FZI(I)/=ZERO)THEN
3353 jg = nsvg(i)
3354 IF(jg>0) THEN
3355C en SPMD : traitement a refaire apres reception noeud remote si JG < 0
3356 icontact(jg)=1
3357 ENDIF
3358 icontact(ix1(i))=1
3359 icontact(ix2(i))=1
3360 icontact(ix3(i))=1
3361 icontact(ix4(i))=1
3362 ENDDO
3363 ENDIF
3364C
3365 IF(iadm/=0)THEN
3366 END IF
3367 IF(iadm>=2)THEN
3368 END IF
3369C
3370 IF(ibcc==0) RETURN
3371C
3372 DO i=1,jlt
3373 IF(pene(i) == zero)cycle
3374 ibcm = ibcc / 8
3375 ibcs = ibcc - 8 * ibcm
3376 IF(ibcs>0) THEN
3377 ig=nsvg(i)
3378C---------obsolate option, no need to update
3379 IF(ig>0.AND.ig<=numnod) THEN
3380C en SPMD : traitement a refaire apres reception noeud remote si JG < 0
3381 CALL ibcoff(ibcs,icodt(ig))
3382 ENDIF
3383 ENDIF
3384 IF(ibcm>0) THEN
3385 ig=ix1(i)
3386 CALL ibcoff(ibcm,icodt(ig))
3387 ig=ix2(i)
3388 CALL ibcoff(ibcm,icodt(ig))
3389 ig=ix3(i)
3390 CALL ibcoff(ibcm,icodt(ig))
3391 ig=ix4(i)
3392 CALL ibcoff(ibcm,icodt(ig))
3393 ENDIF
3394 ENDDO
3395C
3396 RETURN
bitget
integer function bitget(i, n)
Definition bitget.F:37
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
i25sms2
subroutine i25sms2(jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, nin, noint, mskyi_sms, iskyi_sms, nsms, kt, c, cf, dtmini, dti)
Definition i25for3.F:3414
ibcoff
subroutine ibcoff(ibc, icodt)
Definition ibcoff.F:44
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
anim_mod
Definition anim_mod.F:96
h3d_mod
Definition h3d_mod.F:308
outputs_mod
Definition outputs_mod.F:59
parameters_mod
Definition parameters_mod.F:49
pinchtype_mod
Definition pinchtype_mod.F:69
plyxfem_mod
Definition plyxfem_mod.F:31
tri7box::stif_oldfi
type(real_pointer2), dimension(:), allocatable stif_oldfi
Definition tri7box.F:545
tri7box::secnd_frfi
type(real_pointer2), dimension(:), allocatable secnd_frfi
Definition tri7box.F:543
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::lisubsfi
type(int_pointer), dimension(:), allocatable lisubsfi
Definition tri7box.F:501
tri7box::nsvfi
type(int_pointer), dimension(:), allocatable nsvfi
Definition tri7box.F:431
tri7box::intareanfi
type(real_pointer), dimension(:), allocatable intareanfi
Definition tri7box.F:554
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::pene_oldfi
type(real_pointer2), dimension(:), allocatable pene_oldfi
Definition tri7box.F:544

◆ i25sms2()

subroutine i25sms2 ( integer jlt,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) nsvg,
h1,
h2,
h3,
h4,
stif,
integer nin,
integer noint,
mskyi_sms,
integer, dimension(lskyi_sms,*) iskyi_sms,
integer, dimension(*) nsms,
kt,
c,
cf,
dtmini,
dti )

Definition at line 3410 of file i25for3.F.

3414C-----------------------------------------------
3415C M o d u l e s
3416C-----------------------------------------------
3417 USE tri7box
3418 USE message_mod
3419C-----------------------------------------------
3420C I m p l i c i t T y p e s
3421C-----------------------------------------------
3422#include "implicit_f.inc"
3423#include "comlock.inc"
3424C-----------------------------------------------
3425C G l o b a l P a r a m e t e r s
3426C-----------------------------------------------
3427#include "mvsiz_p.inc"
3428C-----------------------------------------------
3429C C o m m o n B l o c k s
3430C-----------------------------------------------
3431#include "parit_c.inc"
3432#include "task_c.inc"
3433#include "sms_c.inc"
3434C-----------------------------------------------
3435C D u m m y A r g u m e n t s
3436C-----------------------------------------------
3437 INTEGER JLT,NIN,NOINT,
3438 . IX1(MVSIZ),IX2(MVSIZ),IX3(MVSIZ),IX4(MVSIZ),NSVG(MVSIZ),
3439 . NSMS(*), ISKYI_SMS(LSKYI_SMS,*)
3440 my_real
3441 . h1(mvsiz),h2(mvsiz),h3(mvsiz),h4(mvsiz),stif(mvsiz),
3442 . mskyi_sms(*), kt(mvsiz), c(mvsiz), cf(mvsiz), dtmini, dti
3443C-----------------------------------------------
3444C L o c a l V a r i a b l e s
3445C-----------------------------------------------
3446 INTEGER I, IG, NISKYL1, NISKYL, NN,JG,IXSS(4),NFIC,J,NS
3447 my_real
3448 . hh(mvsiz),mas, dts,dtm_int,fici,fics(4)
3449C
3450C
3451 niskyl1 = 0
3452 DO i=1,jlt
3453 hh(i)=h1(i)+h2(i)+h3(i)+h4(i)
3454 IF(nsms(i)==0.OR.stif(i)==zero.OR.hh(i)==zero) cycle
3455 IF (h1(i)/=zero) niskyl1 = niskyl1 + 1
3456 IF (h2(i)/=zero) niskyl1 = niskyl1 + 1
3457 IF (h3(i)/=zero) niskyl1 = niskyl1 + 1
3458 IF (h4(i)/=zero) niskyl1 = niskyl1 + 1
3459 ENDDO
3460#include "lockon.inc"
3461 niskyl = nisky_sms
3462 nisky_sms = nisky_sms + niskyl1
3463#include "lockoff.inc"
3464C
3465 IF (niskyl+niskyl1 > lskyi_sms) THEN
3466 CALL ancmsg(msgid=26,anmode=aninfo)
3467 CALL arret(2)
3468 ENDIF
3469C
3470 IF (dtmini>zero) THEN
3471 dtm_int=dtmini
3472 ELSE
3473 dtm_int=dtmins_int
3474 ENDIF
3475
3476 DO i=1,jlt
3477 IF(nsms(i)==0.OR.stif(i)==zero.OR.hh(i)==zero) cycle
3478C
3479 IF(nsms(i)>0)THEN
3480 dts = dtmins/dtfacs
3481 dti=min(dti,dtmins)
3482 ELSE
3483 dts = dtm_int/dtfacs_int
3484 dti=min(dti,dtm_int)
3485 END IF
3486
3487 mas= dts * ( dts * kt(i) + c(i) )
3488 mas = half * max( mas, dts * cf(i) )
3489C
3490 ig =nsvg(i)
3491 IF(ig > 0)THEN
3492 IF(h1(i)/=zero)THEN
3493 niskyl=niskyl+1
3494 mskyi_sms(niskyl)=abs(h1(i))*mas
3495 iskyi_sms(niskyl,1)=ig
3496 iskyi_sms(niskyl,2)=ix1(i)
3497 iskyi_sms(niskyl,3)=ispmd+1
3498 END IF
3499 IF(h2(i)/=zero)THEN
3500 niskyl=niskyl+1
3501 mskyi_sms(niskyl)=abs(h2(i))*mas
3502 iskyi_sms(niskyl,1)=ig
3503 iskyi_sms(niskyl,2)=ix2(i)
3504 iskyi_sms(niskyl,3)=ispmd+1
3505 END IF
3506 IF(h3(i)/=zero)THEN
3507 niskyl=niskyl+1
3508 mskyi_sms(niskyl)=abs(h3(i))*mas
3509 iskyi_sms(niskyl,1)=ig
3510 iskyi_sms(niskyl,2)=ix3(i)
3511 iskyi_sms(niskyl,3)=ispmd+1
3512 END IF
3513 IF(h4(i)/=zero)THEN
3514 niskyl=niskyl+1
3515 mskyi_sms(niskyl)=abs(h4(i))*mas
3516 iskyi_sms(niskyl,1)=ig
3517 iskyi_sms(niskyl,2)=ix4(i)
3518 iskyi_sms(niskyl,3)=ispmd+1
3519 END IF
3520 ELSE
3521 nn = -ig
3522 IF(h1(i)/=zero)THEN
3523 niskyl=niskyl+1
3524 mskyi_sms(niskyl)=abs(h1(i))*mas
3525 iskyi_sms(niskyl,1)=nodamsfi(nin)%P(nn)
3526 iskyi_sms(niskyl,2)=ix1(i)
3527 iskyi_sms(niskyl,3)=procamsfi(nin)%P(nn)
3528 END IF
3529 IF(h2(i)/=zero)THEN
3530 niskyl=niskyl+1
3531 mskyi_sms(niskyl)=abs(h2(i))*mas
3532 iskyi_sms(niskyl,1)=nodamsfi(nin)%P(nn)
3533 iskyi_sms(niskyl,2)=ix2(i)
3534 iskyi_sms(niskyl,3)=procamsfi(nin)%P(nn)
3535 END IF
3536 IF(h3(i)/=zero)THEN
3537 niskyl=niskyl+1
3538 mskyi_sms(niskyl)=abs(h3(i))*mas
3539 iskyi_sms(niskyl,1)=nodamsfi(nin)%P(nn)
3540 iskyi_sms(niskyl,2)=ix3(i)
3541 iskyi_sms(niskyl,3)=procamsfi(nin)%P(nn)
3542 END IF
3543 IF(h4(i)/=zero)THEN
3544 niskyl=niskyl+1
3545 mskyi_sms(niskyl)=abs(h4(i))*mas
3546 iskyi_sms(niskyl,1)=nodamsfi(nin)%P(nn)
3547 iskyi_sms(niskyl,2)=ix4(i)
3548 iskyi_sms(niskyl,3)=procamsfi(nin)%P(nn)
3549 END IF
3550 END IF
3551 ENDDO
3552C
3553 RETURN
message_mod
Definition message_mod.F:1249
tri7box::nodamsfi
type(int_pointer), dimension(:), allocatable nodamsfi
Definition tri7box.F:440
tri7box::procamsfi
type(int_pointer), dimension(:), allocatable procamsfi
Definition tri7box.F:440
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
arret
subroutine arret(nn)
Definition arret.F:87