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