OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i24for3.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!|| i24for3 ../engine/source/interfaces/int24/i24for3.f
25!||--- called by ------------------------------------------------------
26!|| i24mainf ../engine/source/interfaces/int24/i24main.F
27!||--- calls -----------------------------------------------------
28!|| bitget ../engine/source/interfaces/intsort/i20sto.F
29!|| i24_save_sub ../engine/source/interfaces/int24/i24_save_sub.F
30!|| i24ass0 ../engine/source/interfaces/int24/i24for3.F
31!|| i24ass2 ../engine/source/interfaces/int24/i24for3.F
32!|| i24for1_fic ../engine/source/interfaces/int24/i24for3e.F
33!|| i24for1_ficr ../engine/source/interfaces/int24/i24for3e.F
34!|| i24intarea_fic ../engine/source/interfaces/int24/i24intarea_fic.F90
35!|| i24sms2 ../engine/source/interfaces/int24/i24for3.F
36!|| ibcoff ../engine/source/interfaces/interf/ibcoff.F
37!||--- uses -----------------------------------------------------
38!|| h3d_mod ../engine/share/modules/h3d_mod.F
39!|| i24intarea_fic_mod ../engine/source/interfaces/int24/i24intarea_fic.F90
40!|| output_mod ../common_source/modules/output/output_mod.F90
41!|| parameters_mod ../common_source/modules/interfaces/parameters_mod.F
42!|| plyxfem_mod ../engine/share/modules/plyxfem_mod.f
43!|| tri7box ../engine/share/modules/tri7box.F
44!||====================================================================
45 SUBROUTINE i24for3(OUTPUT, JLT ,A ,V ,IBCC ,ICODT ,
46 2 FSAV ,GAP ,FRIC ,MS ,VISC ,
47 3 VISCF ,NOINT ,STFN ,ITAB ,CN_LOC ,
48 4 STIGLO ,STIFN ,STIF ,FSKYI ,ISKY ,
49 5 N1 ,N2 ,N3 ,H1 ,H2 ,
50 6 H3 ,H4 ,FCONT ,PENE ,
51 7 IX1 ,IX2 ,IX3 ,IX4 ,NSVG ,
52 8 IVIS2 ,NELTST ,ITYPTST ,DT2T ,SUBTRIA ,
53 9 GAPV ,INACTI ,INDEX ,NISKYFI ,
54 A KINET ,NEWFRONT ,ISECIN ,NSTRF ,SECFCUM ,
55 B X ,IRECT ,CE_LOC ,MFROT ,IFQ ,
56 C FROT_P ,SECND_FR ,ALPHA0 ,
57 D IBAG ,ICONTACT ,IRTLM ,
58 E VISCN ,VXI ,VYI ,VZI ,MSI ,
59 F KINI ,NIN ,NISUB ,LISUB ,ADDSUBS ,
60 G ADDSUBM ,LISUBS ,LISUBM ,FSAVSUB ,CAND_N ,
61 H ILAGM ,ICURV ,FNCONT ,FTCONT ,NSN ,
62 I X1 ,X2 ,X3 ,X4 ,Y1 ,
63 J Y2 ,Y3 ,Y4 ,Z1 ,Z2 ,
64 K Z3 ,Z4 ,XI ,YI ,ZI ,
65 L IADM ,RCURVI ,RCONTACT ,ACONTACT ,PCONTACT ,
66 M ANGLMI ,PADM ,INTTH ,PHI , FTHE ,
67 N FTHESKYI ,TEMP ,TEMPI ,RSTIF ,IFORM ,
68 O MSKYI_SMS ,ISKYI_SMS ,NSMS ,CAND_N_N ,PENE_OLD ,
69 P STIF_OLD ,MBINFLG ,ILEV ,IGSTI ,KMIN ,
70 P INTPLY ,IPLY ,INOD_PXFEM,NM1 ,NM2 ,
71 R NM3 ,NREBOU ,IRTSE ,NSNE ,IS2SE ,
72 S IS2PT ,MSEGTYP ,JTASK ,ISENSINT ,
73 U FSAVPARIT ,NFT ,H3D_DATA,FRICC ,VISCFFRIC ,
74 V FRIC_COEFS ,T2MAIN_SMS,INTNITSCHE,FORNEQSI,IORTHFRIC ,
75 W FRIC_COEFS2 ,FRICC2 ,VISCFFRIC2,NFORTH ,NFISOT ,
76 X INDEXORTH , INDEXISOT,DIR1 ,DIR2 ,T2FAC_SMS ,
77 Y F_PFIT ,TAGNCONT ,KLOADPINTER,LOADPINTER,LOADP_HYD_INTER,
78 Z TYPSUB ,INFLG_SUBS,INFLG_SUBM ,NINLOADP ,DGAPLOADINT,
79 1 S_LOADPINTER,DIST ,IXX ,INTEREFRIC,INTCAREA,
80 2 PARAMETERS ,PENREF ,KMAX ,S_ADDSUBM, S_LISUBM,S_TYPSUB,
81 3 NISUBMAX,I_STOK,NRTM ,NRTSE ,NSNR)
82C-----------------------------------------------
83C M o d u l e s
84C-----------------------------------------------
85 USE tri7box
86 USE plyxfem_mod
87 USE h3d_mod
88 USE output_mod
89 USE parameters_mod
90 USE i24intarea_fic_mod , ONLY : i24intarea_fic
91C-----------------------------------------------
92C I m p l i c i t T y p e s
93C-----------------------------------------------
94#include "implicit_f.inc"
95#include "comlock.inc"
96C-----------------------------------------------
97C G l o b a l P a r a m e t e r s
98C-----------------------------------------------
99#include "mvsiz_p.inc"
100C-----------------------------------------------
101C C o m m o n B l o c k s
102C-----------------------------------------------
103#include "com01_c.inc"
104#include "com04_c.inc"
105#include "com06_c.inc"
106#include "com08_c.inc"
107#include "scr05_c.inc"
108#include "scr07_c.inc"
109#include "scr11_c.inc"
110#include "scr14_c.inc"
111#include "scr16_c.inc"
112#include "scr18_c.inc"
113#include "sms_c.inc"
114#include "parit_c.inc"
115#include "param_c.inc"
116#include "impl1_c.inc"
117C-----------------------------------------------
118C D u m m y A r g u m e n t s
119C-----------------------------------------------
120 TYPE(output_), INTENT(INOUT) :: OUTPUT
121 INTEGER S_ADDSUBM ! Size of ADDSUBM
122 INTEGER S_LISUBM ! Size of LISUBM
123 INTEGER S_TYPSUB ! Size of TYPSUB
124 INTEGER NISUBMAX ! Size of ISENSINT
125 INTEGER I_STOK ! Number of contact candidates / Size of FSAVPARIT
126 INTEGER NELTST,ITYPTST,JLT,IBCC,IBCM,IBCS,IVIS2,INACTI,IBAG,NIN,
127 . INTNITSCHE,IORTHFRIC,NFORTH ,NFISOT,
128 . NSN,ICODT(*), ITAB(*), ISKY(*), KINET(*),IRTLM(2,NSN),
129 . MFROT, IFQ, NOINT,NEWFRONT,ISECIN, NSTRF(*),
130 . IRECT(4,*),ICONTACT(*), CAND_N(*),
131 . KINI(*),SUBTRIA(MVSIZ),MBINFLG(*),ILEV,
132 . NISKYFI,IADM,INTTH,IFORM,CAND_N_N(*),INTPLY,
133 . MSEGTYP(*),TAGNCONT(NLOADP_HYD_INTER,NUMNOD),
134 . TYPSUB(S_TYPSUB),JTASK
135 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
136 . CN_LOC(MVSIZ),CE_LOC(MVSIZ),INDEX(MVSIZ),NSVG(MVSIZ),
137 . NISUB, LISUB(*), ADDSUBS(*), ADDSUBM(S_ADDSUBM), LISUBS(*),
138 . LISUBM(*),ILAGM,ICURV(3), NREBOU,
139 . ISKYI_SMS(*), NSMS(*),IGSTI,IPLY(4,*),INOD_PXFEM(*),
140 . IRTSE(5,*),NSNE ,IS2SE(2,*),IS2PT(*), ISENSINT(NISUBMAX+1),NFT,T2MAIN_SMS(6,*),
141 . INDEXORTH(MVSIZ) ,INDEXISOT(MVSIZ),INFLG_SUBS(*), INFLG_SUBM(*)
142 INTEGER , INTENT(IN) :: NINLOADP,S_LOADPINTER
143 INTEGER , INTENT(IN) :: KLOADPINTER(NINTER+1),LOADPINTER(S_LOADPINTER),
144 . LOADP_HYD_INTER(NLOADP_HYD)
145 INTEGER , INTENT(IN) :: IXX(MVSIZ,13)
146 INTEGER , INTENT(IN) :: INTEREFRIC, INTCAREA
147 INTEGER , INTENT(IN) :: NRTM, NRTSE, NSNR
148 my_real
149 . STIGLO,FROT_P(*), X(3,*),
150 . A(3,*), MS(*), V(3,*), FSAV(*),FCONT(3,*),
151 . SECND_FR(6,*),ALPHA0,
152 . GAP, FRIC,VISC,VISCF,VIS,DT2T,STFN(*),STIFN(*),
153 . FSKYI(LSKYI,NFSKYI),FSAVSUB(NTHVKI,*), FNCONT(3,*),FTCONT(3,*),
154 . MSKYI_SMS(*),KMIN
155 my_real
156 . stif(mvsiz),
157 . gapv(mvsiz), pene(mvsiz),
158 . secfcum(7,numnod,nsect),
159 .
160 . viscn(*),
161 . vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),msi(mvsiz),
162 . x1(mvsiz),y1(mvsiz),z1(mvsiz),
163 . x2(mvsiz),y2(mvsiz),z2(mvsiz),
164 . x3(mvsiz),y3(mvsiz),z3(mvsiz),
165 . x4(mvsiz),y4(mvsiz),z4(mvsiz),
166 . xi(mvsiz),yi(mvsiz),zi(mvsiz),
167 . n1(mvsiz), n2(mvsiz), n3(mvsiz),
168 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
169 . nm1(mvsiz), nm2(mvsiz), nm3(mvsiz),
170 . phi(mvsiz), fthe(*),ftheskyi(lskyi),temp(*), tempi(mvsiz),
171 . fsavparit(nisub+1,11,*),
172 . fricc(mvsiz),viscffric(mvsiz),fric_coefs(mvsiz,10),forneqsi(mvsiz,3),
173 . fric_coefs2(mvsiz,10),fricc2(mvsiz),viscffric2(mvsiz),
174 . dir1(mvsiz,3),dir2(mvsiz,3),t2fac_sms(*)
175
176
177 my_real
178 . rcurvi(*), rcontact(*), acontact(*),
179 . pcontact(*),padm, anglmi(*),rstif,
180 . pene_old(5,nsn),stif_old(2,nsn),f_pfit
181 my_real , INTENT(IN) :: dgaploadint(s_loadpinter)
182 my_real , INTENT(INOUT) :: dist(mvsiz)
183 my_real, DIMENSION(MVSIZ), INTENT(IN) :: penref
184 my_real, INTENT(IN) :: kmax
185
186 TYPE(h3d_database) :: H3D_DATA
187 TYPE (PARAMETERS_) ,INTENT(IN):: PARAMETERS
188C-----------------------------------------------
189C L o c a l V a r i a b l e s
190C-----------------------------------------------
191 INTEGER I, J1, IG, J, JG, K0, NBINTER, K1S, K, IE, NN,
192 . N,IMS2,ITAG,NN1,NN2,NN3,
193 . NN4,IFIT,PP,PPL,ITYPSUB,ISS1,ISS2,IMS1,IGRN
194 my_real
195 . FXI(MVSIZ), FYI(MVSIZ), FZI(MVSIZ), FNI(MVSIZ),
196 . FXT(MVSIZ),FYT(MVSIZ),FZT(MVSIZ),
197 . FX1(MVSIZ), FX2(MVSIZ), FX3(MVSIZ), FX4(MVSIZ),
198 . FY1(MVSIZ), FY2(MVSIZ), FY3(MVSIZ), FY4(MVSIZ),
199 . FZ1(MVSIZ), FZ2(MVSIZ), FZ3(MVSIZ), FZ4(MVSIZ),
200 . VIS2(MVSIZ), XMU(MVSIZ),STIF0(MVSIZ),
201 . VX(MVSIZ), VY(MVSIZ), VZ(MVSIZ), VN(MVSIZ),
202 . XP(MVSIZ), YP(MVSIZ), ZP(MVSIZ),EFRIC_L(MVSIZ),
203 . AA,
204 . V2, DT1INV, FAC,FF,ALPHA,BETA,
205 . FX, FY, FZ, MAS2,DTI,FT,FN,FMAX,FTN,
206 . FACM1, ECONTT, ECONVT, ECONTDT,
207 . FSAV1, FSAV2, FSAV3, FSAV4, FSAV5, FSAV6, FSAV8,
208 . FSAV9, FSAV10, FSAV11, FSAV12, FSAV13, FSAV14, FSAV15,
209 . FSAV22, FSAV23, FSAV24, FFO,FSAV29,
210 . VV,AX1,AX2,AY1,AY2,AZ1,AZ2,AX,AY,AZ,
211 . AREA,P,VV1,VV2,DMU ,
212 . VN1(3),VN2(3),VN3(3),VN4(4),
213 . csa ,sna ,alphaf ,ftt1 ,ftt2 ,
214 . an(nforth) ,bn(nforth) ,nep1 ,nep2 ,nep ,c11 ,c22 ,ans ,ep ,
215 . signc,dgapload,gapp,efric_ls,efric_lm
216 my_real
217 . prec,facv
218 my_real
219 . st1(mvsiz),st2(mvsiz),st3(mvsiz),st4(mvsiz),stv(mvsiz),
220 . kt(mvsiz),c(mvsiz),cf(mvsiz),dpene(mvsiz),
221 . ks(mvsiz),k1(mvsiz),k2(mvsiz),k3(mvsiz),k4(mvsiz),
222 . cs(mvsiz),c1(mvsiz),c2(mvsiz),c3(mvsiz),c4(mvsiz),
223 . cx,cy,cfi,aux,
224 . phi1(mvsiz), phi2(mvsiz), phi3(mvsiz), phi4(mvsiz),
225 . dtmini,stif0_imp(mvsiz),penn,tiny,pendr,
226 . fnnit1,fnnit2,fnnit3,fnnitsche,fnitxt(mvsiz),fnityt(mvsiz),fnitzt(mvsiz),
227 . xmu2(mvsiz),facn,fnon,stifkt(mvsiz),fackt,arean_fic
228 INTEGER JSUB, KSUB, JJ, KK, IN, NSUB, NCFIT
229
230 my_real
231 . fsavsub1(25,nisub),impx,impy,impz,vnm(mvsiz),sfac,stmin,
232 . hh,finc,ddp,prec1
233
234 DOUBLE PRECISION RP1, RP2
235 INTEGER TAGIP(MVSIZ)
236C
237 INTEGER BITGET
238 EXTERNAL BITGET
239C-----------------------------------------------
240 DTMINI=zero
241 fac = zero
242 IF (iresp==1) THEN
243 prec = fiveem4
244 tiny = em15
245 prec1= em05
246 ELSE
247 prec = em10
248 tiny = em30
249 prec1= zero
250 ENDIF
251 IF(dt1>zero)THEN
252 dt1inv = one/dt1
253 ELSE
254 dt1inv =zero
255 ENDIF
256 econtt = zero
257 econvt = zero
258 econtdt = zero
259 DO i=1,jlt
260 stif0(i) = stif(i)
261 stifkt(i) = stif(i)
262 ENDDO
263 efric_l(1:jlt) = zero
264C---------------------
265C COURBURE FIXE
266C---------------------
267 IF(icurv(1)==1)THEN
268 ELSEIF(icurv(1)==2)THEN
269 ELSEIF(icurv(1) == 3)THEN
270 ENDIF
271 ncfit = icurv(2)
272 IF (inacti==-1.AND.impl_s==0.AND.ncfit>0) THEN
273 ifit = 1
274 ELSE
275 ifit = 0
276 END IF
277 fnon = ep02
278 fmax = ep03
279 IF (igsti==-1) fmax = kmax
280C-------------------------------------------
281C PENE Offset removed to i24dst3.F
282C-------------------------------------------
283C
284 IF(intply == 0) THEN
285 DO i=1,jlt
286 IF(pene(i) == zero)THEN
287 h1(i) = zero
288 h2(i) = zero
289 h3(i) = zero
290 h4(i) = zero
291 fni(i)= zero
292C
293 fxi(i)=zero
294 fyi(i)=zero
295 fzi(i)=zero
296C
297 fx1(i)=zero
298 fy1(i)=zero
299 fz1(i)=zero
300C
301 fx2(i)=zero
302 fy2(i)=zero
303 fz2(i)=zero
304C
305 fx3(i)=zero
306 fy3(i)=zero
307 fz3(i)=zero
308C
309 fx4(i)=zero
310 fy4(i)=zero
311 fz4(i)=zero
312 ELSE
313 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
314 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
315 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
316 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
317 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
318 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
319 vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
320 vnm(i) = nm1(i)*vx(i) + nm2(i)*vy(i) + nm3(i)*vz(i)
321 ENDIF
322 ENDDO
323 ELSE
324 DO i=1,jlt
325 IF(pene(i) == zero)THEN
326 h1(i) = zero
327 h2(i) = zero
328 h3(i) = zero
329 h4(i) = zero
330 fni(i)= zero
331C
332 fxi(i)=zero
333 fyi(i)=zero
334 fzi(i)=zero
335C
336 fx1(i)=zero
337 fy1(i)=zero
338 fz1(i)=zero
339C
340 fx2(i)=zero
341 fy2(i)=zero
342 fz2(i)=zero
343C
344 fx3(i)=zero
345 fy3(i)=zero
346 fz3(i)=zero
347C
348 fx4(i)=zero
349 fy4(i)=zero
350 fz4(i)=zero
351 ELSE
352 jj = iply(1,i)
353 nn1 =inod_pxfem(ix1(i))
354 nn2 =inod_pxfem(ix2(i))
355 nn3 =inod_pxfem(ix3(i))
356 nn4 =inod_pxfem(ix4(i))
357C
358 vn1(1) = v(1,ix1(i))
359 vn1(2) = v(2,ix1(i))
360 vn1(3) = v(3,ix1(i))
361C
362 vn2(1) = v(1,ix2(i))
363 vn2(2) = v(2,ix2(i))
364 vn2(3) = v(3,ix2(i))
365C
366 vn3(1) = v(1,ix3(i))
367 vn3(2) = v(2,ix3(i))
368 vn3(3) = v(3,ix3(i))
369C
370 vn4(1) = v(1,ix4(i))
371 vn4(2) = v(2,ix4(i))
372 vn4(3) = v(3,ix4(i))
373C
374 IF(iplyxfem /= 2)THEN
375C
376C possible energy creation after full delamination
377C - seen with IPLYXFEM=1 -
378 jj = iply(1,i)
379 IF(nn1 > 0 .AND. jj > 0) THEN
380 vn1(1) = vn1(1) + ply(jj)%V(1,nn1)
381 vn1(2) = vn1(2) + ply(jj)%V(2,nn1)
382 vn1(3) = vn1(3) + ply(jj)%V(3,nn1)
383 ENDIF
384 jj = iply(2,i)
385 IF(nn2 > 0 .AND. jj > 0) THEN
386 vn2(1) = vn2(1) + ply(jj)%V(1,nn2)
387 vn2(2) = vn2(2) + ply(jj)%V(2,nn2)
388 vn2(3) = vn2(3) + ply(jj)%V(3,nn2)
389 ENDIF
390 jj = iply(3,i)
391 IF(nn3 > 0 .AND. jj > 0) THEN
392 vn3(1) = vn3(1) + ply(jj)%V(1,nn3)
393 vn3(2) = vn3(2) + ply(jj)%V(2,nn3)
394 vn3(3) = vn3(3) + ply(jj)%V(3,nn3)
395 ENDIF
396 jj = iply(4,i)
397 IF(nn4 > 0 .AND. jj > 0) THEN
398 vn4(1) = vn4(1) + ply(jj)%V(1,nn4)
399 vn4(2) = vn4(2) + ply(jj)%V(2,nn4)
400 vn4(3) = vn4(3) + ply(jj)%V(3,nn4)
401 ENDIF
402 END IF
403
404 vx(i) = vxi(i) - h1(i)*vn1(1) - h2(i)*vn2(1)
405 . - h3(i)*vn3(1) - h4(i)*vn4(1)
406 vy(i) = vyi(i) - h1(i)*vn1(2) - h2(i)*vn2(2)
407 . - h3(i)*vn3(2) - h4(i)*vn4(2)
408 vz(i) = vzi(i) - h1(i)*vn1(3) - h2(i)*vn2(3)
409 . - h3(i)*vn3(3) - h4(i)*vn4(3)
410 vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
411 vnm(i) = nm1(i)*vx(i) + nm2(i)*vy(i) + nm3(i)*vz(i)
412 ENDIF
413 ENDDO
414
415 ENDIF
416C-------------------------------------------
417C stifness reduction to avoid positive
418C force jump and energy generation
419C-------------------------------------------
420 IF (ifit >0.OR.(inacti==-1.AND.impl_s==0)) THEN
421 DO i=1,jlt
422 IF(pene(i) == zero)cycle
423 stif(i) = stif0(i)
424 IF(stiglo<=zero) stif(i) = half*stif(i)
425C---- nonlinear stif
426 IF(penref(i) > zero) THEN
427 pendr = (pene(i)/penref(i))**2
428 fackt = min(fmax,(one+three*fnon*pendr))
429 facn = min(fmax,(one+fnon*pendr))
430 stifkt(i)= stifkt(i)*fackt
431 stif(i)= stif(i)*facn
432 END IF
433 fni(i)= -stif(i) * pene(i)
434 ENDDO
435 ELSEIF (igsti==-1.AND.impl_s==0) THEN
436c----- first version w/o force jump treatment
437 DO i=1,jlt
438 IF(pene(i) == zero)cycle
439 stif(i) = stif0(i)
440C---- nonlinear stif
441 IF(penref(i) > zero) THEN
442 pendr = (pene(i)/penref(i))**2
443 facn = min(fmax,(one+fnon*pendr))
444 stif(i)= min(kmax,stif(i)*facn)
445 fackt = three*(stif(i)/stif0(i)-one)
446 stifkt(i)= max(stif(i),stifkt(i)*fackt)
447 END IF
448 fni(i)= -stif(i) * pene(i)
449 ENDDO
450 ELSEIF(impl_s>0.AND.igsti==6)THEN
451C-------0.001*STIF initially (first contact), special case(rebound) 0.0001*STIF but >=KMIN------
452C--- NREBOU<0 rebound in this cycle, : -1 -> at first contact; -2 -> after first contact
453 IF (nrebou == -1) THEN
454 sfac = em04
455 ELSE
456 sfac = em03
457 END IF
458 DO i=1,jlt
459 IF(pene(i) == zero)cycle
460 jg = nsvg(i)
461 n = cand_n_n(i)
462 IF(jg > 0)THEN
463 IF (stif_old(1,n)==zero.OR.nrebou==-1) THEN
464 stif_old(1,n) = sfac*stif0(i)
465 ELSEIF (nrebou <-1) THEN
466 stmin = em04*stif0(i)
467 stif_old(1,n) = max(stmin,em01*stif_old(1,n))
468 END IF
469 stif_old(1,n) = max(kmin,stif_old(1,n))
470 stif(i) = stif_old(1,n)
471 ELSE
472 jg = -jg
473 IF (stif_oldfi(nin)%P(1,jg)==zero.OR.nrebou==-1) THEN
474 stif_oldfi(nin)%P(1,jg)= sfac*stif0(i)
475 ELSEIF (nrebou <-1) THEN
476 stmin = em04*stif0(i)
477 stif_oldfi(nin)%P(1,jg)=
478 + max(stmin,em01*stif_oldfi(nin)%P(1,jg))
479 END IF
480 stif_oldfi(nin)%P(1,jg)= max(kmin,stif_oldfi(nin)%P(1,jg))
481 stif(i) = stif_oldfi(nin)%P(1,jg)
482 END IF
483 END DO
484 IF(inconv >=0)THEN
485C---------if multi-contact in the same groupe---
486 DO i=1,jlt
487 IF(pene(i) == zero)cycle
488 jg = nsvg(i)
489 n = cand_n_n(i)
490 IF(jg > 0)THEN
491 stif0_imp(i) = stif_old(1,n)
492 ELSE
493 stif0_imp(i) = stif_oldfi(nin)%P(1,-jg)
494 END IF
495 END DO
496 DO i=1,jlt
497 IF(pene(i) == zero)cycle
498 jg = nsvg(i)
499 n = cand_n_n(i)
500 IF(jg > 0)THEN
501C---------if multi-contact in the same groupe---
502 IF(pene_old(2,n)/= zero.OR.pene_old(5,n)/= zero)THEN
503 IF(pene(i) > pene_old(2,n) )THEN
504 facm1 = pene(i)/max(em20,pene_old(2,n))
505 IF (stif0_imp(i)/stif0(i) <em02) THEN
506 facm1 = onep2*facm1
507 penn = three_half
508 ELSE
509 penn = onep05
510 END IF
511 facm1=min(penn,facm1)
512 stif(i) = facm1*stif0_imp(i)
513C-----------never higher than STIF0----------------
514 stif(i) = min(stif(i),stif0(i))
515 ENDIF
516C------------quand PENE_OLD(2,N)=PENE_OLD(1,N)
517 IF(inconv ==1 )THEN
518 pene_old(1,n) = pene(i)
519 stif_old(1,n) = stif(i)
520 ENDIF
521 ELSEIF(inconv ==1 )THEN
522c new impact, possible multi impact
523#include "lockon.inc"
524 pene_old(1,n) = max(pene_old(1,n),pene(i))
525 stif_old(1,n) = min(stif_old(1,n),stif(i))
526#include "lockoff.inc"
527 ENDIF
528 ELSE
529 jg = -jg
530 IF(pene_oldfi(nin)%P(2,jg)/= zero.OR.
531 . pene_oldfi(nin)%P(5,jg)/= zero)THEN
532 IF(pene(i) > pene_oldfi(nin)%P(2,jg) )THEN
533 facm1 = pene(i)/max(em20,pene_oldfi(nin)%P(2,jg))
534 IF (stif0_imp(i)/stif0(i) <em02) THEN
535 facm1 = onep2*facm1
536 penn = three_half
537 ELSE
538 penn = onep05
539 END IF
540 facm1=min(penn,facm1)
541 stif(i) = facm1*stif0_imp(i)
542 stif(i) = min(stif(i),stif0(i))
543 ENDIF
544 IF(inconv ==1 )THEN
545 pene_oldfi(nin)%P(1,jg) = pene(i)
546 stif_oldfi(nin)%P(1,jg) = stif(i)
547 END IF !(INCONV ==1 )THEN
548 ELSEIF(inconv ==1 )THEN
549#include "lockon.inc"
550 pene_oldfi(nin)%P(1,jg) =
551 * max(pene_oldfi(nin)%P(1,jg),pene(i))
552 stif_oldfi(nin)%P(1,jg) =
553 * min(stif_oldfi(nin)%P(1,jg),stif(i))
554#include "lockoff.inc"
555 ENDIF
556 ENDIF
557 ENDDO
558C-------to pass for implicit
559 DO i=1,jlt
560 IF(pene(i) == zero)cycle
561 jg = nsvg(i)
562 n = cand_n_n(i)
563 IF(jg > 0)THEN
564 stif_old(2,n) = stif(i)
565 ELSE
566 stif_oldfi(nin)%P(2,-jg) = stif(i)
567 END IF
568 END DO
569 END IF !(INCONV>=0)THEN
570 DO i=1,jlt
571 IF(pene(i) == zero)cycle
572 IF(stiglo<=zero)THEN
573 stif(i) = half*stif(i)
574ccc ECONTT = ECONTT + HALF*STIF(I)*PENE(I)**2
575 ELSEIF(stif(i)/=zero)THEN
576 stif(i) = stiglo
577ccc ECONTT = ECONTT + HALF*STIGLO*PENE(I)**2
578 ENDIF
579 fni(i)= -stif(i) * pene(i)
580 ENDDO
581
582C------explicit or implicit+Istif/=6
583 ELSE
584 DO i=1,jlt
585 IF(pene(i) == zero)cycle
586 dpene(i) = max(zero,-vnm(i)*dt1)
587 jg = nsvg(i)
588 n = cand_n_n(i)
589 IF(jg > 0)THEN
590 IF(tt /= zero.AND.(pene_old(2,n)/= zero.OR.pene_old(5,n)/= zero))THEN
591 ddp = pene(i)-pene_old(2,n)-dpene(i)
592 IF(pene(i) > pene_old(2,n) .AND.ddp>zero)THEN
593 rp1 = pene_old(2,n)/pene(i)
594 rp2 = dpene(i)/pene(i)
595 stif(i) = stif_old(2,n)*rp1+ stif(i)*rp2
596 ELSEIF(.NOT.(pene_old(2,n)== zero.AND.stif_old(2,n)==zero))THEN
597C-----------first contact STIF_OLD(2,N) should not be zero
598 stif(i) = stif_old(2,n)
599 END IF
600 IF(inconv ==1 )THEN
601 pene_old(1,n) = pene(i)
602 stif_old(1,n) = stif(i)
603 ENDIF
604 ELSEIF(inconv ==1 )THEN
605c new impact, possible multi impact
606#include "lockon.inc"
607 pene_old(1,n) = max(pene_old(1,n),pene(i))
608 stif_old(1,n) = max(stif_old(1,n),stif(i))
609#include "lockoff.inc"
610 ENDIF
611 ELSE
612 jg = -jg
613 IF(tt /= zero.AND.(pene_oldfi(nin)%P(2,jg)/= zero.OR.
614 . pene_oldfi(nin)%P(5,jg)/= zero))THEN
615 ddp = pene(i)-pene_oldfi(nin)%P(2,jg)-dpene(i)
616 IF(pene(i) > pene_oldfi(nin)%P(2,jg) .AND.ddp>zero)THEN
617 rp1 = pene_oldfi(nin)%P(2,jg)/pene(i)
618 rp2 = dpene(i)/pene(i)
619 stif(i) = stif_oldfi(nin)%P(2,jg)*rp1+ stif(i)*rp2
620 ELSEIF(.NOT.(pene_oldfi(nin)%P(2,jg)== zero.AND.
621 . stif_oldfi(nin)%P(2,jg)==zero))THEN
622 stif(i) = stif_oldfi(nin)%P(2,jg)
623 END IF
624 IF(inconv ==1 )THEN
625 pene_oldfi(nin)%P(1,jg) = pene(i)
626 stif_oldfi(nin)%P(1,jg) = stif(i)
627 ENDIF
628 ELSEIF(inconv ==1 )THEN
629#include "lockon.inc"
630 pene_oldfi(nin)%P(1,jg) =
631 * max(pene_oldfi(nin)%P(1,jg),pene(i))
632 stif_oldfi(nin)%P(1,jg) =
633 * max(stif_oldfi(nin)%P(1,jg),stif(i))
634#include "lockoff.inc"
635 ENDIF
636 ENDIF
637 ENDDO
638C-------------------------------------------
639 DO i=1,jlt
640 IF(pene(i) == zero)cycle
641 IF(stiglo<=zero)THEN
642 stif(i) = half*stif(i)
643 ELSEIF(stif(i)/=zero)THEN
644 stif(i) = stiglo
645 ENDIF
646 fni(i)= -stif(i) * pene(i)
647
648 ENDDO
649C-------------------------------------------
650 ENDIF !IF(IMPL_S>0.AND.IGSTI==6)
651C-------------------------------------------
652C Contact energy (elastic)
653C-------------------------------------------
654 DO i=1,jlt
655 IF(pene(i) == zero)cycle
656 econtt = econtt+half*stif(i)*pene(i)**2
657 ENDDO
658C
659 IF(intnitsche > 0 ) THEN
660 DO i=1,jlt
661 IF(pene(i) == zero)cycle
662 fnnit1 = forneqsi(i,1)*n1(i)
663 fnnit2 = forneqsi(i,2)*n2(i)
664 fnnit3 = forneqsi(i,3)*n3(i)
665 fnnitsche = fnnit1 + fnnit2 + fnnit3
666C
667 fni(i) = min(zero,fni(i) - half*fnnitsche)
668C
669 ENDDO
670 ENDIF
671C---------------------------------
672C DAMPING + FRIC
673C---------------------------------
674C goto 999
675 IF(visc/=zero)THEN
676 IF(ivis2==0)THEN
677C---------------------------------
678C VISC QUAD TYPE V227
679C---------------------------------
680 DO i=1,jlt
681 vis2(i) = two * stif(i) * msi(i)
682 ENDDO
683C---------------------------------
684C---------------------------------
685 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
686 DO i=1,jlt
687 IF(pene(i) == zero)cycle
688 fac = stif(i) / max(em30,stif(i))
689 vis = sqrt(vis2(i))
690 ff = fac * visc * vis
691 stif(i) = stif0(i) + ff * dt1inv
692 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
693 ffo = ff
694 ff = ff * vn(i)
695C ECONVT = ECONVT + FF * VN(I) * DT1
696 fni(i) = fni(i) + ff
697 ENDDO
698 ELSE
699 DO i=1,jlt
700 IF(pene(i) == zero)cycle
701 fac = stif(i) / max(em30,stif(i))
702 vis = sqrt(vis2(i))
703 c(i)= fac * visc * vis
704 kt(i)= stif0(i)
705 stif(i) = stif0(i) + c(i) * dt1inv
706 ff = c(i) * vn(i)
707C ECONVT = ECONVT + FF * VN(I) * DT1
708 fni(i) = fni(i) + ff
709 cf(i) = fac*sqrt(viscffric(i))*vis
710 stif(i) = max(stif(i) ,cf(i)*dt1inv)
711 ENDDO
712 ENDIF
713 ELSEIF(ivis2==1.OR.ivis2==2)THEN
714C---------------------------------
715C VISC QUAD TYPE V227
716C---------------------------------
717 IF(ivis2==1)THEN
718 facv = two
719 ELSE
720 facv = four
721 END IF
722 DO i=1,jlt
723 IF(pene(i) == zero)cycle
724 mas2 = ms(ix1(i))*h1(i)
725 . + ms(ix2(i))*h2(i)
726 . + ms(ix3(i))*h3(i)
727 . + ms(ix4(i))*h4(i)
728 vis2(i) = facv* stif(i) * msi(i) * mas2 /
729 . max(em30,msi(i)+mas2)
730 ENDDO
731C---------------------------------
732C---------------------------------
733 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
734 DO i=1,jlt
735 IF(pene(i) == zero)cycle
736 fac = stif(i) / max(em30,stif(i))
737 vis = sqrt(vis2(i))
738 ff = fac * visc * vis
739 stif(i) = stif0(i) + ff * dt1inv
740 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
741 ffo = ff
742 ff = ff * vn(i)
743C ECONTDT = ECONTDT + FF * VN(I) * DT1
744 fni(i) = fni(i) + ff
745 ENDDO
746 ELSE
747 DO i=1,jlt
748 IF(pene(i) == zero)cycle
749 fac = stif(i) / max(em30,stif(i))
750 vis = sqrt(vis2(i))
751 c(i)= fac * visc * vis
752 kt(i)= stif0(i)
753 stif(i) = stif(i) + c(i) * dt1inv
754 ff = c(i) * vn(i)
755C ECONTDT = ECONTDT + FF * VN(I) * DT1
756 fni(i) = fni(i) + ff
757 cf(i) = fac*sqrt(viscffric(i))*vis
758 stif(i) = max(stif(i) ,cf(i)*dt1inv)
759 ENDDO
760 ENDIF
761 ELSEIF(ivis2==3)THEN
762C---------------------------------
763C VISC QUAD = 0
764C---------------------------------
765 DO i=1,jlt
766 IF(pene(i) == zero)cycle
767 vis2(i) = two * stif(i) * msi(i)
768 fac = stif(i) / max(em30,stif(i))
769 vis = sqrt(vis2(i))
770 ff = fac * visc * vis
771 stif(i) = stif0(i) + two* ff * dt1inv
772 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
773 ff = ff * vn(i)
774C ECONTDT = ECONTDT+ FF * VN(I) * DT1
775 fni(i) = fni(i) + ff
776 ENDDO
777 ELSEIF(ivis2==4)THEN
778C---------------------------------
779C VISC = 0
780C---------------------------------
781 DO i=1,jlt
782 IF(pene(i) == zero)cycle
783 fac = stif(i) / max(em30,stif(i))
784 vis2(i) = two* stif(i) * msi(i)
785 vis = sqrt(vis2(i))
786 stif(i) = max(stif0(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
787 ENDDO
788 ELSEIF(ivis2==5)THEN
789C---------------------------------
790C Visc = 2m/dt => For visc <1, stable: Dt <2m/visc = Dt
791C M = m1*m2/m1+m2 for visc = 1, elastic shock
792C For visc = 0.5, elastic collision
793C---------------------------------
794 DO i=1,jlt
795 IF(pene(i) == zero)cycle
796 mas2 = ms(ix1(i))*h1(i)
797 . + ms(ix2(i))*h2(i)
798 . + ms(ix3(i))*h3(i)
799 . + ms(ix4(i))*h4(i)
800 vis2(i) = two* stif(i) * msi(i)
801 vis = 2. * visc * dt1inv * msi(i) * mas2 /
802 . max(em30,msi(i)+mas2)
803 stif(i) = max(stif0(i) ,fac*sqrt(viscffric(i)*vis2(i))*dt1inv)
804 ff = vis * vn(i)
805 econtdt = econtdt + min(zero,ff-fni(i)) * vn(i) * dt1
806 fni(i) = min(fni(i),ff)
807 ENDDO
808 ELSE
809 ENDIF
810 ELSE
811 DO i=1,jlt
812 IF(viscffric(i)/=zero)THEN
813
814 IF(ivis2==0)THEN
815C---------------------------------
816C VISC QUAD TYPE V227
817C---------------------------------
818 vis2(i) = two * stif(i) * msi(i)
819C---------------------------------
820C---------------------------------
821 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
822 IF(pene(i) == zero)cycle
823 fac = stif(i) / max(em30,stif(i))
824 vis = sqrt(vis2(i))
825 ff = fac * visc * vis
826 stif(i) = stif0(i) + ff * dt1inv
827 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
828 ffo = ff
829 ff = ff * vn(i)
830C ECONTDT = ECONTDT + FF * VN(I) * DT1
831 fni(i) = fni(i) + ff
832 ELSE
833 IF(pene(i) == zero)cycle
834 fac = stif(i) / max(em30,stif(i))
835 vis = sqrt(vis2(i))
836 c(i)= fac * visc * vis
837 kt(i)= stif0(i)
838 stif(i) = stif0(i) + c(i) * dt1inv
839 ff = c(i) * vn(i)
840C ECONTDT = ECONTDT + FF * VN(I) * DT1
841 fni(i) = fni(i) + ff
842 cf(i) = fac*sqrt(viscffric(i))*vis
843 stif(i) = max(stif(i) ,cf(i)*dt1inv)
844 ENDIF
845 ELSEIF(ivis2==1.OR.ivis2==2)THEN
846C---------------------------------
847C VISC QUAD TYPE V227
848C---------------------------------
849 IF(ivis2==1)THEN
850 facv = two
851 ELSE
852 facv = four
853 END IF
854 IF(pene(i) == zero)cycle
855 mas2 = ms(ix1(i))*h1(i)
856 . + ms(ix2(i))*h2(i)
857 . + ms(ix3(i))*h3(i)
858 . + ms(ix4(i))*h4(i)
859 vis2(i) = facv* stif(i) * msi(i) * mas2 /
860 . max(em30,msi(i)+mas2)
861C---------------------------------
862C---------------------------------
863 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
864 IF(pene(i) == zero)cycle
865 fac = stif(i) / max(em30,stif(i))
866 vis = sqrt(vis2(i))
867 ff = fac * visc * vis
868 stif(i) = stif0(i) + ff * dt1inv
869 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
870 ffo = ff
871 ff = ff * vn(i)
872C ECONTDT = ECONTDT + FF * VN(I) * DT1
873 fni(i) = fni(i) + ff
874 ELSE
875 IF(pene(i) == zero)cycle
876 fac = stif(i) / max(em30,stif(i))
877 vis = sqrt(vis2(i))
878 c(i)= fac * visc * vis
879 kt(i)= stif0(i)
880 stif(i) = stif(i) + c(i) * dt1inv
881 ff = c(i) * vn(i)
882C ECONTDT = ECONTDT + FF * VN(I) * DT1
883 fni(i) = fni(i) + ff
884 cf(i) = fac*sqrt(viscffric(i))*vis
885 stif(i) = max(stif(i) ,cf(i)*dt1inv)
886 ENDIF
887 ELSEIF(ivis2==2)THEN
888C---------------------------------
889C VISC QUAD TYPE
890C---------------------------------
891 IF(pene(i) == zero)cycle
892 vis2(i) = two* stif(i) * msi(i)
893 fac = stif(i) / max(em30,stif(i))
894 vis = sqrt(vis2(i))
895 ff = fac * visc * vis
896 stif(i) = stif0(i) + two * ff * dt1inv
897 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
898 ff = ff * vn(i)
899C ECONTDT = ECONTDT + FF * VN(I) * DT1
900 fni(i) = fni(i) + ff
901 ELSEIF(ivis2==3)THEN
902C---------------------------------
903C VISC QUAD = 0
904C---------------------------------
905 IF(pene(i) == zero)cycle
906 vis2(i) = two * stif(i) * msi(i)
907 fac = stif(i) / max(em30,stif(i))
908 vis = sqrt(vis2(i))
909 ff = fac * visc * vis
910 stif(i) = stif0(i) + two* ff * dt1inv
911 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
912 ff = ff * vn(i)
913C ECONTDT = ECONTDT + FF * VN(I) * DT1
914 fni(i) = fni(i) + ff
915 ELSEIF(ivis2==4)THEN
916C---------------------------------
917C VISC = 0
918C---------------------------------
919 IF(pene(i) == zero)cycle
920 fac = stif(i) / max(em30,stif(i))
921 vis2(i) = two* stif(i) * msi(i)
922 vis = sqrt(vis2(i))
923 stif(i) = max(stif0(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
924 ELSEIF(ivis2==5)THEN
925C---------------------------------
926C Visc = 2m/dt => For visc <1, stable: Dt <2m/visc = Dt
927C M = m1*m2/m1+m2 for visc = 1, elastic shock
928C For visc = 0.5, elastic collision
929C---------------------------------
930 IF(pene(i) == zero)cycle
931 mas2 = ms(ix1(i))*h1(i)
932 . + ms(ix2(i))*h2(i)
933 . + ms(ix3(i))*h3(i)
934 . + ms(ix4(i))*h4(i)
935 vis2(i) = two* stif(i) * msi(i)
936 vis = 2. * visc * dt1inv * msi(i) * mas2 /
937 . max(em30,msi(i)+mas2)
938 stif(i) = max(stif0(i) ,fac*sqrt(viscffric(i)*vis2(i))*dt1inv)
939 ff = vis * vn(i)
940 econtdt = econtdt + min(zero,ff-fni(i)) * vn(i) * dt1
941 fni(i) = min(fni(i),ff)
942 ELSE
943 ENDIF
944
945
946 ELSE
947cbb initialization of the vis2 array to avoid problems
948 vis2(i) = zero
949 ENDIF
950 ENDDO
951 ENDIF
952
953C---------------------------------
954C SAUVEGARDE DE L'IMPULSION NORMALE
955C---------------------------------
956 fsav1 = zero
957 fsav2 = zero
958 fsav3 = zero
959 fsav8 = zero
960 fsav9 = zero
961 fsav10= zero
962 fsav11= zero
963 IF (ilev==2) THEN
964 DO i=1,jlt
965 IF(pene(i) == zero)cycle
966 ie=ce_loc(i)
967 ims2 = bitget(mbinflg(ie),1)
968 fxi(i)=n1(i)*fni(i)
969 fyi(i)=n2(i)*fni(i)
970 fzi(i)=n3(i)*fni(i)
971 impx=fxi(i)*dt12
972 impy=fyi(i)*dt12
973 impz=fzi(i)*dt12
974 fsav8 =fsav8 +abs(impx)
975 fsav9 =fsav9 +abs(impy)
976 fsav10=fsav10+abs(impz)
977 IF (ims2 >0 ) THEN
978 fsav1 =fsav1 -impx
979 fsav2 =fsav2 -impy
980 fsav3 =fsav3 -impz
981 fsav11=fsav11-fni(i)*dt12
982 ELSE
983 fsav1 =fsav1 +impx
984 fsav2 =fsav2 +impy
985 fsav3 =fsav3 +impz
986 fsav11=fsav11+fni(i)*dt12
987 END IF
988 ENDDO
989 ELSE
990 DO i=1,jlt
991 IF(pene(i) == zero)cycle
992 fxi(i)=n1(i)*fni(i)
993 fyi(i)=n2(i)*fni(i)
994 fzi(i)=n3(i)*fni(i)
995 impx=fxi(i)*dt12
996 impy=fyi(i)*dt12
997 impz=fzi(i)*dt12
998 fsav1 =fsav1 +impx
999 fsav2 =fsav2 +impy
1000 fsav3 =fsav3 +impz
1001 fsav8 =fsav8 +abs(impx)
1002 fsav9 =fsav9 +abs(impy)
1003 fsav10=fsav10+abs(impz)
1004 fsav11=fsav11+fni(i)*dt12
1005 ENDDO
1006 END IF !(ILEV==2) THEN
1007#include "lockon.inc"
1008 fsav(1)=fsav(1)+fsav1
1009 fsav(2)=fsav(2)+fsav2
1010 fsav(3)=fsav(3)+fsav3
1011 fsav(8)=fsav(8)+fsav8
1012 fsav(9)=fsav(9)+fsav9
1013 fsav(10)=fsav(10)+fsav10
1014 fsav(11)=fsav(11)+fsav11
1015#include "lockoff.inc"
1016C
1017 IF(isensint(1)/=0) THEN
1018 DO i=1,jlt
1019 IF(pene(i) == zero)cycle
1020 fsavparit(1,1,i+nft) = fxi(i)
1021 fsavparit(1,2,i+nft) = fyi(i)
1022 fsavparit(1,3,i+nft) = fzi(i)
1023 ENDDO
1024 ENDIF
1025C---------------------------------
1026 IF((anim_v(12)+outp_v(12)+h3d_data%N_VECT_PCONT>0.AND.
1027 . ((tt>=output%TANIM.AND.tt<=output%TANIM_STOP).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
1028 . (manim>=4.AND.manim<=15).OR.h3d_data%MH3D/=0))
1029 . .OR.h3d_data%N_VECT_PCONT_MAX>0)THEN
1030 IF (inconv==1) THEN
1031#include "lockon.inc"
1032 DO i=1,jlt
1033 IF(pene(i) == zero)cycle
1034 fncont(1,ix1(i)) =fncont(1,ix1(i)) + fxi(i)*h1(i)
1035 fncont(2,ix1(i)) =fncont(2,ix1(i)) + fyi(i)*h1(i)
1036 fncont(3,ix1(i)) =fncont(3,ix1(i)) + fzi(i)*h1(i)
1037 fncont(1,ix2(i)) =fncont(1,ix2(i)) + fxi(i)*h2(i)
1038 fncont(2,ix2(i)) =fncont(2,ix2(i)) + fyi(i)*h2(i)
1039 fncont(3,ix2(i)) =fncont(3,ix2(i)) + fzi(i)*h2(i)
1040 fncont(1,ix3(i)) =fncont(1,ix3(i)) + fxi(i)*h3(i)
1041 fncont(2,ix3(i)) =fncont(2,ix3(i)) + fyi(i)*h3(i)
1042 fncont(3,ix3(i)) =fncont(3,ix3(i)) + fzi(i)*h3(i)
1043 fncont(1,ix4(i)) =fncont(1,ix4(i)) + fxi(i)*h4(i)
1044 fncont(2,ix4(i)) =fncont(2,ix4(i)) + fyi(i)*h4(i)
1045 fncont(3,ix4(i)) =fncont(3,ix4(i)) + fzi(i)*h4(i)
1046 jg = nsvg(i)
1047 IF(jg>0) THEN
1048C In SPMD: Treatment to be redone after reception node Remote if JG <0
1049 IF (jg >numnod) THEN
1050 ig = jg - numnod
1051 CALL i24for1_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
1052 + ig ,fxi(i) ,fyi(i) ,fzi(i) ,fncont ,
1053 + -1 )
1054 ELSE
1055 fncont(1,jg)=fncont(1,jg)- fxi(i)
1056 fncont(2,jg)=fncont(2,jg)- fyi(i)
1057 fncont(3,jg)=fncont(3,jg)- fzi(i)
1058 END IF
1059 ELSE ! cas noeud remote en SPMD
1060 jg = -jg
1061 IF(isedge_fi(nin)%P(jg)==1)THEN
1062 CALL i24for1_ficr(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
1063 + jg ,fxi(i) ,fyi(i) ,fzi(i) ,fnconti(nin)%P(1,1) ,
1064 + -1 ,nin )
1065 ELSE
1066 fnconti(nin)%P(1,jg)=fnconti(nin)%P(1,jg)-fxi(i)
1067 fnconti(nin)%P(2,jg)=fnconti(nin)%P(2,jg)-fyi(i)
1068 fnconti(nin)%P(3,jg)=fnconti(nin)%P(3,jg)-fzi(i)
1069 ENDIF
1070 ENDIF
1071 ENDDO
1072#include "lockoff.inc"
1073 END IF !(INCONV==1) THEN
1074 ENDIF
1075C
1076C---------------------------------
1077C SORTIES TH PAR SOUS INTERFACE
1078C---------------------------------
1079 IF(nisub/=0)THEN
1080 DO jsub=1,nisub
1081 DO j=1,25
1082 fsavsub1(j,jsub)=zero
1083 END DO
1084 ENDDO
1085 DO i=1,jlt
1086 IF(pene(i) == zero)cycle
1087 nn = nsvg(i)
1088 IF(nn>0)THEN
1089 in=cn_loc(i)
1090 IF (msegtyp(ce_loc(i)) < 0) THEN
1091 ie= - msegtyp(ce_loc(i))
1092 ELSE
1093 ie = ce_loc(i)
1094 ENDIF
1095 IF(ie > nrtm) ie=ie-nrtm
1096 jj =addsubs(in)
1097 kk =addsubm(ie)
1098 DO WHILE(jj<addsubs(in+1))
1099
1100 jsub=lisubs(jj)
1101 itypsub = typsub(jsub)
1102
1103 IF(itypsub == 1 ) THEN ! Defining specific inter
1104 iss1 = bitget(inflg_subs(jj),0)
1105 iss2 = bitget(inflg_subs(jj),1)
1106 igrn = bitget(inflg_subs(jj),2)
1107 ksub=lisubm(kk)
1108 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1109 ims1 = bitget(inflg_subm(kk),0)
1110 ims2 = bitget(inflg_subm(kk),1)
1111 IF(ksub==jsub)THEN
1112 IF(.NOT.(((ims1 == 1 .OR. ims2 == 1) .AND. igrn == 1).OR.
1113 . (ims1 == 1 .AND. iss2 == 1).OR.
1114 . (ims2 == 1 .AND. iss1 == 1))) THEN
1115 kk=kk+1
1116 ksub=lisubm(kk)
1117 cycle
1118 END IF
1119C
1120 impx=fxi(i)*dt12
1121 impy=fyi(i)*dt12
1122 impz=fzi(i)*dt12
1123 IF(ims2 > 0)THEN
1124 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
1125 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
1126 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
1127 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
1128 ELSE
1129 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1130 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1131 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1132 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1133 END IF
1134C
1135 IF(isensint(jsub+1)/=0) THEN
1136 IF(ims2 > 0)THEN
1137 fsavparit(jsub+1,1,i) = -fxi(i)
1138 fsavparit(jsub+1,2,i) = -fyi(i)
1139 fsavparit(jsub+1,3,i) = -fzi(i)
1140 ELSE
1141 fsavparit(jsub+1,1,i) = fxi(i)
1142 fsavparit(jsub+1,2,i) = fyi(i)
1143 fsavparit(jsub+1,3,i) = fzi(i)
1144 END IF
1145 ENDIF
1146C
1147 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1148 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1149 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1150C
1151 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1152
1153 IF(parameters%INTCAREA > 0) THEN
1154 IF(nn <=numnod) THEN
1155 fsavsub1(25,jsub) = fsavsub1(25,jsub) + parameters%INTAREAN(nn)
1156 ELSE
1157 ig = nn - numnod
1158 CALL i24intarea_fic(irtse ,nsne ,is2se ,is2pt ,ig ,
1159 + nrtse , numnod ,parameters%INTAREAN, arean_fic)
1160 fsavsub1(25,jsub) = fsavsub1(25,jsub) + arean_fic
1161 ENDIF
1162 ENDIF
1163C
1164 END IF
1165
1166 kk=kk+1
1167 ksub=lisubm(kk)
1168 ENDDO
1169 jj=jj+1
1170
1171 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surf : secondary side
1172
1173 impx=fxi(i)*dt12
1174 impy=fyi(i)*dt12
1175 impz=fzi(i)*dt12
1176
1177 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1178 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1179 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1180
1181 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1182 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1183 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1184
1185 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1186
1187 IF(isensint(jsub+1)/=0) THEN
1188 fsavparit(jsub+1,1,i+nft) = fxi(i)
1189 fsavparit(jsub+1,2,i+nft) = fyi(i)
1190 fsavparit(jsub+1,3,i+nft) = fzi(i)
1191 ENDIF
1192
1193 IF(parameters%INTCAREA > 0) THEN
1194 IF(nn <=numnod) THEN
1195 fsavsub1(25,jsub) = fsavsub1(25,jsub) + parameters%INTAREAN(nn)
1196 ELSE
1197 ig = nn - numnod
1198 CALL i24intarea_fic(irtse ,nsne ,is2se ,is2pt ,ig ,
1199 + nrtse , numnod , parameters%INTAREAN, arean_fic)
1200 fsavsub1(25,jsub) = fsavsub1(25,jsub) + arean_fic
1201 ENDIF
1202 ENDIF
1203
1204 jj=jj+1
1205
1206 ELSEIF(itypsub == 3) THEN !Inter =0 : collecting forces from all inter with 2 surfs
1207
1208
1209
1210 iss2 = bitget(inflg_subs(jj),0)
1211 iss1 = bitget(inflg_subs(jj),1)
1212 ksub=lisubm(kk)
1213
1214 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1215 ims2 = bitget(inflg_subm(kk),0)
1216 ims1 = bitget(inflg_subm(kk),1)
1217
1218 IF(ksub==jsub)THEN
1219 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
1220 . (ims2 == 1 .AND. iss1 == 1))) THEN
1221 kk=kk+1
1222 ksub=lisubm(kk)
1223 cycle
1224 END IF
1225
1226 impx=fxi(i)*dt12
1227 impy=fyi(i)*dt12
1228 impz=fzi(i)*dt12
1229
1230 IF(ims2 > 0)THEN
1231 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
1232 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
1233 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
1234 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
1235 ELSE
1236 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1237 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1238 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1239 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1240 END IF
1241C
1242 IF(isensint(jsub+1)/=0) THEN
1243 IF(ims2 > 0)THEN
1244 fsavparit(jsub+1,1,i+nft) = -fxi(i)
1245 fsavparit(jsub+1,2,i+nft) = -fyi(i)
1246 fsavparit(jsub+1,3,i+nft) = -fzi(i)
1247 ELSE
1248 fsavparit(jsub+1,1,i+nft) = fxi(i)
1249 fsavparit(jsub+1,2,i+nft) = fyi(i)
1250 fsavparit(jsub+1,3,i+nft) = fzi(i)
1251 END IF
1252 ENDIF
1253C
1254 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1255 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1256 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1257C
1258 IF(parameters%INTCAREA > 0) THEN
1259 IF(nn <=numnod) THEN
1260 fsavsub1(25,jsub) = fsavsub1(25,jsub) + parameters%INTAREAN(nn)
1261 ELSE
1262 ig = nn - numnod
1263 CALL i24intarea_fic(irtse ,nsne ,is2se ,is2pt ,ig ,
1264 + nrtse , numnod , parameters%INTAREAN, arean_fic)
1265 fsavsub1(25,jsub) = fsavsub1(25,jsub) + arean_fic
1266 ENDIF
1267 ENDIF
1268 END IF
1269
1270 kk=kk+1
1271 ksub=lisubm(kk)
1272 ENDDO
1273 jj=jj+1
1274 ENDIF
1275
1276 END DO
1277 END IF
1278
1279
1280 IF (msegtyp(ce_loc(i)) < 0) THEN
1281 ie= - msegtyp(ce_loc(i))
1282 ELSE
1283 ie = ce_loc(i)
1284 ENDIF
1285 IF(ie > nrtm) ie=ie-nrtm
1286
1287 ! Compiler issue workaround,
1288 ! Loop packed in subroutine to avoid optimization breaking the code.
1289 CALL i24_save_sub(numnod,mvsiz,nisub,s_addsubm,s_lisubm,s_typsub,nisubmax,i_stok,
1290 * ie,itypsub,nin,i,nn,nft,
1291 * addsubm,lisubm,typsub,
1292 * parameters%INTAREAN,parameters%INTCAREA,isensint,
1293 * fxi,fyi,fzi,fni,dt12,
1294 * fsavsub1,fsavparit ,nrtse,
1295 * irtse,nsne,is2se ,is2pt,nsnr)
1296
1297 END DO
1298
1299
1300 IF(nspmd>1) THEN
1301 DO i=1,jlt
1302 IF(pene(i) == zero)cycle
1303 nn = nsvg(i)
1304 IF(nn<0)THEN
1305 nn = -nn
1306 IF (msegtyp(ce_loc(i)) < 0) THEN
1307 ie= - msegtyp(ce_loc(i))
1308 ELSE
1309 ie = ce_loc(i)
1310 ENDIF
1311 jj =addsubsfi(nin)%P(nn)
1312 kk =addsubm(ie)
1313 DO WHILE(jj<addsubsfi(nin)%P(nn+1))
1314 jsub=lisubsfi(nin)%P(jj)
1315 itypsub = typsub(jsub)
1316
1317 IF(itypsub == 1 ) THEN ! Defining specific inter
1318
1319 iss1 = bitget(inflg_subsfi(nin)%P(jj),0)
1320 iss2 = bitget(inflg_subsfi(nin)%P(jj),1)
1321 igrn = bitget(inflg_subsfi(nin)%P(jj),2)
1322 ksub=lisubm(kk)
1323 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1324 ims1 = bitget(inflg_subm(kk),0)
1325 ims2 = bitget(inflg_subm(kk),1)
1326 IF(ksub==jsub)THEN
1327 IF(.NOT.(((ims1 == 1 .OR. ims2 == 1) .AND. igrn == 1).OR.
1328 . (ims1 == 1 .AND. iss2 == 1).OR.
1329 . (ims2 == 1 .AND. iss1 == 1))) THEN
1330 kk=kk+1
1331 ksub=lisubm(kk)
1332 cycle
1333 END IF
1334 impx=fxi(i)*dt12
1335 impy=fyi(i)*dt12
1336 impz=fzi(i)*dt12
1337C
1338 IF(ims2 > 0)THEN
1339 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
1340 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
1341 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
1342 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
1343 ELSE
1344 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1345 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1346 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1347 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1348 END IF
1349C
1350 IF(isensint(jsub+1)/=0) THEN
1351 IF(ims2 > 0)THEN
1352 fsavparit(jsub+1,1,i) = -fxi(i)
1353 fsavparit(jsub+1,2,i) = -fyi(i)
1354 fsavparit(jsub+1,3,i) = -fzi(i)
1355 ELSE
1356 fsavparit(jsub+1,1,i) = fxi(i)
1357 fsavparit(jsub+1,2,i) = fyi(i)
1358 fsavparit(jsub+1,3,i) = fzi(i)
1359 END IF
1360 ENDIF
1361C
1362 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1363 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1364 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1365C
1366 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1367C
1368 IF(parameters%INTCAREA > 0) THEN
1369 IF(isedge_fi(nin)%P(nn)==1)THEN
1370 CALL i24intarea_fic(irtse_fi(nin)%P ,nsnr ,is2se_fi(nin)%P ,is2pt_fi(nin)%P ,nn ,
1371 + nsnr , nsnr , intareanfi(nin)%P, arean_fic)
1372 fsavsub1(25,jsub) = fsavsub1(25,jsub) + arean_fic
1373 ELSE ! cas noeud remote en SPMD non edge
1374 fsavsub1(25,jsub) = fsavsub1(25,jsub) + intareanfi(nin)%P(nn)
1375 ENDIF
1376 ENDIF
1377 END IF
1378
1379 kk=kk+1
1380 ksub=lisubm(kk)
1381 ENDDO
1382 jj=jj+1
1383
1384 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surf : secondary side
1385
1386 impx=fxi(i)*dt12
1387 impy=fyi(i)*dt12
1388 impz=fzi(i)*dt12
1389
1390 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1391 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1392 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1393
1394 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1395 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1396 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1397
1398 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1399
1400 IF(isensint(jsub+1)/=0) THEN
1401 fsavparit(jsub+1,1,i+nft) = fxi(i)
1402 fsavparit(jsub+1,2,i+nft) = fyi(i)
1403 fsavparit(jsub+1,3,i+nft) = fzi(i)
1404 ENDIF
1405
1406 IF(parameters%INTCAREA > 0) THEN
1407 IF(isedge_fi(nin)%P(nn)==1)THEN
1408 CALL i24intarea_fic(irtse_fi(nin)%P ,nsnr ,is2se_fi(nin)%P ,is2pt_fi(nin)%P ,nn ,
1409 + nsnr , nsnr , intareanfi(nin)%P, arean_fic)
1410 fsavsub1(25,jsub) = fsavsub1(25,jsub) + arean_fic
1411 ELSE ! cas noeud remote en SPMD non edge
1412 fsavsub1(25,jsub) = fsavsub1(25,jsub) + intareanfi(nin)%P(nn)
1413 ENDIF
1414 ENDIF
1415
1416 jj=jj+1
1417
1418 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfaces
1419
1420 iss2 = bitget(inflg_subsfi(nin)%P(jj),0)
1421 iss1 = bitget(inflg_subsfi(nin)%P(jj),1)
1422 ksub=lisubm(kk)
1423 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1424 ims2 = bitget(inflg_subm(kk),0)
1425 ims1 = bitget(inflg_subm(kk),1)
1426 IF(ksub==jsub)THEN
1427 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
1428 . (ims2 == 1 .AND. iss1 == 1))) THEN
1429 kk=kk+1
1430 ksub=lisubm(kk)
1431 cycle
1432 END IF
1433
1434 impx=fxi(i)*dt12
1435 impy=fyi(i)*dt12
1436 impz=fzi(i)*dt12
1437
1438 IF(ims2 > 0)THEN
1439 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
1440 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
1441 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
1442 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
1443 ELSE
1444 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1445 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1446 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1447 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1448 END IF
1449C
1450 IF(isensint(jsub+1)/=0) THEN
1451 IF(ims2 > 0)THEN
1452 fsavparit(jsub+1,1,i+nft) = -fxi(i)
1453 fsavparit(jsub+1,2,i+nft) = -fyi(i)
1454 fsavparit(jsub+1,3,i+nft) = -fzi(i)
1455 ELSE
1456 fsavparit(jsub+1,1,i+nft) = fxi(i)
1457 fsavparit(jsub+1,2,i+nft) = fyi(i)
1458 fsavparit(jsub+1,3,i+nft) = fzi(i)
1459 END IF
1460 ENDIF
1461
1462 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1463 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1464 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1465
1466 IF(parameters%INTCAREA > 0) THEN
1467 IF(isedge_fi(nin)%P(nn)==1)THEN
1468 CALL i24intarea_fic(irtse_fi(nin)%P ,nsnr ,is2se_fi(nin)%P ,is2pt_fi(nin)%P ,nn ,
1469 + nsnr , nsnr , intareanfi(nin)%P, arean_fic)
1470 fsavsub1(25,jsub) = fsavsub1(25,jsub) + arean_fic
1471 ELSE ! cas noeud remote en SPMD non edge
1472 fsavsub1(25,jsub) = fsavsub1(25,jsub) + intareanfi(nin)%P(nn)
1473 ENDIF
1474 ENDIF
1475
1476 ENDIF
1477 kk=kk+1
1478 ksub=lisubm(kk)
1479 ENDDO
1480 jj=jj+1
1481
1482 ENDIF
1483
1484 END DO
1485 END IF
1486 END DO
1487 END IF
1488 END IF
1489
1490C---------------------------------
1491C NEW FRICTION MODELS
1492C---------------------------------
1493
1494C Friction coefficient computation
1495C++++++++++++++++++++++++++++++++++++++++
1496
1497 IF(iorthfric == 0) THEN
1498C++ Isotropic Friction
1499
1500 IF (mfrot==0) THEN
1501C--- Coulomb friction
1502 DO i=1,jlt
1503 IF(pene(i) == 0) THEN
1504 xmu(i) = zero
1505 cycle
1506 ENDIF
1507 xmu(i) = fricc(i)
1508 ENDDO
1509 ELSEIF (mfrot==1) THEN
1510C--- Viscous friction
1511 DO i=1,jlt
1512 IF(pene(i) == 0) THEN
1513 xmu(i) = zero
1514 cycle
1515 ENDIF
1516 ie=ce_loc(i)
1517 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1518 v2 = (vx(i) - n1(i)*aa)**2
1519 . + (vy(i) - n2(i)*aa)**2
1520 . + (vz(i) - n3(i)*aa)**2
1521 vv = sqrt(max(em30,v2))
1522 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1523 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1524 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1525 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1526 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1527 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1528 ax = ay1*az2 - az1*ay2
1529 ay = az1*ax2 - ax1*az2
1530 az = ax1*ay2 - ay1*ax2
1531 area = half*sqrt(ax*ax+ay*ay+az*az)
1532 p = -fni(i)/area
1533 xmu(i) = fricc(i) + (fric_coefs(i,1) + fric_coefs(i,4)*p ) * p
1534 . +(fric_coefs(i,2) + fric_coefs(i,3)*p) * vv + fric_coefs(i,5)*v2
1535 xmu(i) = max(xmu(i),em30)
1536 ENDDO
1537 ELSEIF(mfrot==2)THEN
1538C--- Darmstad LAW
1539 DO i=1,jlt
1540 IF(pene(i) == 0) THEN
1541 xmu(i) = zero
1542 cycle
1543 ENDIF
1544 ie=ce_loc(i)
1545 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1546 v2 = (vx(i) - n1(i)*aa)**2
1547 . + (vy(i) - n2(i)*aa)**2
1548 . + (vz(i) - n3(i)*aa)**2
1549 vv = sqrt(max(em30,v2))
1550 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1551 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1552 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1553 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1554 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1555 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1556 ax = ay1*az2 - az1*ay2
1557 ay = az1*ax2 - ax1*az2
1558 az = ax1*ay2 - ay1*ax2
1559 area = half*sqrt(ax*ax+ay*ay+az*az)
1560 p = -fni(i)/area
1561 xmu(i) = fricc(i)
1562 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
1563 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
1564 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
1565 xmu(i) = max(xmu(i),em30)
1566 ENDDO
1567 ELSEIF (mfrot==3) THEN
1568C--- Renard LAW
1569 DO i=1,jlt
1570 IF(pene(i) == 0) THEN
1571 xmu(i) = zero
1572 cycle
1573 ENDIF
1574 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1575 v2 = (vx(i) - n1(i)*aa)**2
1576 . + (vy(i) - n2(i)*aa)**2
1577 . + (vz(i) - n3(i)*aa)**2
1578 vv = sqrt(max(em30,v2))
1579 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
1580 dmu = fric_coefs(i,3)-fric_coefs(i,1)
1581 vv1 = vv / fric_coefs(i,5)
1582 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
1583 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
1584 dmu = fric_coefs(i,4)-fric_coefs(i,3)
1585 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
1586 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
1587 ELSE
1588 dmu = fric_coefs(i,2)-fric_coefs(i,4)
1589 vv2 = (vv - fric_coefs(i,6))**2
1590 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
1591 ENDIF
1592 xmu(i) = max(xmu(i),em30)
1593 ENDDO
1594 ELSEIF(mfrot==4)THEN
1595C--- Exponential decay model
1596 DO i=1,jlt
1597 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1598 v2 = (vx(i) - n1(i)*aa)**2
1599 . + (vy(i) - n2(i)*aa)**2
1600 . + (vz(i) - n3(i)*aa)**2
1601 vv = sqrt(max(em30,v2))
1602 xmu(i) = fric_coefs(i,1)
1603 . + (fricc(i)-fric_coefs(i,1))*exp(-fric_coefs(i,2)*vv)
1604 xmu(i) = max(xmu(i),em30)
1605 ENDDO
1606 ENDIF
1607
1608 ELSE
1609C++ Orthotropic Friction
1610
1611 IF (mfrot==0) THEN
1612C--- Coulomb friction
1613#include "vectorize.inc"
1614 DO k=1,nfisot
1615 i = indexisot(k)
1616 xmu(i) = fricc(i)
1617 ENDDO
1618
1619#include "vectorize.inc"
1620 DO k=1,nforth
1621 i = indexorth(k)
1622 xmu(i) = fricc(i)
1623 xmu2(i) = fricc2(i)
1624 IF(xmu(i)<=em30) THEN
1625 xmu(i) = em30
1626 dir1(i,1:3) = zero
1627 an(k) = zero
1628 ELSE
1629 an(k)=xmu(i)**2
1630 an(k)=one/an(k)
1631 ENDIF
1632 IF(xmu2(i)<=em30) THEN
1633 xmu2(i) = em30
1634 dir2(i,1:3) = zero
1635 bn(k) = zero
1636 ELSE
1637 bn(k)=xmu2(i)**2
1638 bn(k)=one/bn(k)
1639 ENDIF
1640 ENDDO
1641 ELSEIF (mfrot==1) THEN
1642C--- Viscous friction
1643#include "vectorize.inc"
1644 DO k=1,nfisot ! isotropic friction couples
1645 i = indexisot(k)
1646 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1647 v2 = (vx(i) - n1(i)*aa)**2
1648 . + (vy(i) - n2(i)*aa)**2
1649 . + (vz(i) - n3(i)*aa)**2
1650 vv = sqrt(max(em30,v2))
1651
1652 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1653 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1654 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1655 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1656 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1657 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1658 ax = ay1*az2 - az1*ay2
1659 ay = az1*ax2 - ax1*az2
1660 az = ax1*ay2 - ay1*ax2
1661 area = half*sqrt(ax*ax+ay*ay+az*az)
1662 p = -fni(i)/area
1663 xmu(i) = fricc(i) + (fric_coefs(i,1) + fric_coefs(i,4)*p ) * p
1664 . +(fric_coefs(i,2) + fric_coefs(i,3)*p) * vv + fric_coefs(i,5)*v2
1665 xmu(i) = max(xmu(i),em30)
1666 ENDDO
1667
1668#include "vectorize.inc"
1669 DO k=1,nforth ! Orthotropic friction couples
1670 i = indexorth(k)
1671 ie=ce_loc(i)
1672c
1673 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1674 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1675 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1676 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1677 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1678 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1679 ax = ay1*az2 - az1*ay2
1680 ay = az1*ax2 - ax1*az2
1681 az = ax1*ay2 - ay1*ax2
1682 area = half*sqrt(ax*ax+ay*ay+az*az)
1683 p = -fni(i)/area
1684c
1685 v2 = vx(i)*dir1(k,1) +vy(i)*dir1(k,2)+vz(i)*dir1(k,3)
1686 vv = max(em30,v2)
1687 xmu(i) = fricc(i) + (fric_coefs(i,1) + fric_coefs(i,4)*p ) * p
1688 . +(fric_coefs(i,2) + fric_coefs(i,3)*p) * vv + fric_coefs(i,5)*vv*vv
1689
1690 v2 = vx(i)*dir2(k,1) +vy(i)*dir2(k,2)+vz(i)*dir2(k,3)
1691 vv = max(em30,v2)
1692 xmu2(i) = fricc2(i) + (fric_coefs2(i,1) + fric_coefs2(i,4)*p ) * p
1693 . +(fric_coefs2(i,2) + fric_coefs2(i,3)*p) * vv + fric_coefs2(i,5)*vv*vv
1694
1695 xmu(i) = max(xmu(i),em30)
1696 xmu2(i) = max(xmu2(i),em30)
1697
1698 ENDDO
1699
1700#include "vectorize.inc"
1701 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
1702 i = indexorth(k)
1703 IF(xmu(i)<=em30) THEN
1704 xmu(i) = em30
1705 dir1(i,1:3) = zero
1706 an(k) = zero
1707 ELSE
1708 an(k)=xmu(i)**2
1709 an(k)=one/an(k)
1710 ENDIF
1711 IF(xmu2(i)<=em30) THEN
1712 xmu2(i) = em30
1713 dir2(i,1:3) = zero
1714 bn(k) = zero
1715 ELSE
1716 bn(k)=xmu2(i)**2
1717 bn(k)=one/bn(k)
1718 ENDIF
1719 ENDDO
1720
1721
1722 ELSEIF(mfrot==2)THEN
1723C--- Loi Darmstad
1724#include "vectorize.inc"
1725 DO k=1,nfisot ! isotropic friction couples
1726 i = indexisot(k)
1727 ie=ce_loc(i)
1728 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1729 v2 = (vx(i) - n1(i)*aa)**2
1730 . + (vy(i) - n2(i)*aa)**2
1731 . + (vz(i) - n3(i)*aa)**2
1732 vv = sqrt(max(em30,v2))
1733 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1734 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1735 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1736 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1737 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1738 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1739 ax = ay1*az2 - az1*ay2
1740 ay = az1*ax2 - ax1*az2
1741 az = ax1*ay2 - ay1*ax2
1742 area = half*sqrt(ax*ax+ay*ay+az*az)
1743 p = -fni(i)/area
1744 xmu(i) = fricc(i)
1745 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
1746 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
1747 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
1748 xmu(i) = max(xmu(i),em30)
1749 ENDDO
1750c
1751#include "vectorize.inc"
1752 DO k=1,nforth ! Orthotropic friction couples
1753 i = indexorth(k)
1754 ie=ce_loc(i)
1755c
1756 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1757 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1758 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1759 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1760 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1761 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1762 ax = ay1*az2 - az1*ay2
1763 ay = az1*ax2 - ax1*az2
1764 az = ax1*ay2 - ay1*ax2
1765 area = half*sqrt(ax*ax+ay*ay+az*az)
1766 p = -fni(i)/area
1767c
1768 v2 = vx(i)*dir1(k,1) +vy(i)*dir1(k,2)+vz(i)*dir1(k,3)
1769 vv = max(em30,v2)
1770 xmu(i) = fricc(i)
1771 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
1772 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
1773 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
1774c
1775 v2 = vx(i)*dir2(k,1) +vy(i)*dir2(k,2)+vz(i)*dir2(k,3)
1776 vv = max(em30,v2)
1777 xmu2(i) = fricc2(i)
1778 . + fric_coefs2(i,1)*exp(fric_coefs2(i,2)*vv)*p*p
1779 . + fric_coefs2(i,3)*exp(fric_coefs2(i,4)*vv)*p
1780 . + fric_coefs2(i,5)*exp(fric_coefs2(i,6)*vv)
1781
1782 xmu(i) = max(xmu(i),em30)
1783 xmu2(i) = max(xmu2(i),em30)
1784
1785 ENDDO
1786
1787#include "vectorize.inc"
1788 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
1789 i = indexorth(k)
1790 IF(xmu(i)<=em30) THEN
1791 xmu(i) = em30
1792 dir1(i,1:3) = zero
1793 an(k) = zero
1794 ELSE
1795 an(k)=xmu(i)**2
1796 an(k)=one/an(k)
1797 ENDIF
1798 IF(xmu2(i)<=em30) THEN
1799 xmu2(i) = em30
1800 dir2(i,1:3) = zero
1801 bn(k) = zero
1802 ELSE
1803 bn(k)=xmu2(i)**2
1804 bn(k)=one/bn(k)
1805 ENDIF
1806 ENDDO
1807
1808 ELSEIF (mfrot==3) THEN
1809C--- Loi Renard
1810#include "vectorize.inc"
1811 DO k=1,nfisot ! isotropic friction couples
1812 i = indexisot(k)
1813 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1814 v2 = (vx(i) - n1(i)*aa)**2
1815 . + (vy(i) - n2(i)*aa)**2
1816 . + (vz(i) - n3(i)*aa)**2
1817 vv = sqrt(max(em30,v2))
1818 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
1819 dmu = fric_coefs(i,3)-fric_coefs(i,1)
1820 vv1 = vv / fric_coefs(i,5)
1821 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
1822 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
1823 dmu = fric_coefs(i,4)-fric_coefs(i,3)
1824 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
1825 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
1826 ELSE
1827 dmu = fric_coefs(i,2)-fric_coefs(i,4)
1828 vv2 = (vv - fric_coefs(i,6))**2
1829 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
1830 ENDIF
1831 xmu(i) = max(xmu(i),em30)
1832 ENDDO
1833
1834#include "vectorize.inc"
1835 DO k=1,nforth ! Orthotropic friction couples
1836 i = indexorth(k)
1837c
1838 v2 = vx(i)*dir1(k,1) +vy(i)*dir1(k,2)+vz(i)*dir1(k,3)
1839 vv = max(em30,v2)
1840 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
1841 dmu = fric_coefs(i,3)-fric_coefs(i,1)
1842 vv1 = vv / fric_coefs(i,5)
1843 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
1844 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
1845 dmu = fric_coefs(i,4)-fric_coefs(i,3)
1846 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
1847 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
1848 ELSE
1849 dmu = fric_coefs(i,2)-fric_coefs(i,4)
1850 vv2 = (vv - fric_coefs(i,6))**2
1851 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
1852 ENDIF
1853
1854 v2 = vx(i)*dir2(k,1) +vy(i)*dir2(k,2)+vz(i)*dir2(k,3)
1855 vv = max(em30,v2)
1856 IF(vv>=0.AND.vv<=fric_coefs2(i,5)) THEN
1857 dmu = fric_coefs2(i,3)-fric_coefs2(i,1)
1858 vv1 = vv / fric_coefs2(i,5)
1859 xmu2(i) = fric_coefs2(i,1)+ dmu*vv1*(two-vv1)
1860 ELSEIF(vv>fric_coefs2(i,5).AND.vv<fric_coefs2(i,6)) THEN
1861 dmu = fric_coefs2(i,4)-fric_coefs2(i,3)
1862 vv1 = (vv - fric_coefs2(i,5))/(fric_coefs2(i,6)-fric_coefs2(i,5))
1863 xmu2(i) = fric_coefs2(i,3)+ dmu * (three-two*vv1)*vv1**2
1864 ELSE
1865 dmu = fric_coefs2(i,2)-fric_coefs2(i,4)
1866 vv2 = (vv - fric_coefs2(i,6))**2
1867 xmu2(i) = fric_coefs2(i,2) - dmu / (one + dmu*vv2)
1868 ENDIF
1869 xmu(i) = max(xmu(i),em30)
1870 xmu2(i) = max(xmu2(i),em30)
1871 ENDDO
1872
1873#include "vectorize.inc"
1874 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
1875 i = indexorth(k)
1876 IF(xmu(i)<=em30) THEN
1877 xmu(i) = em30
1878 dir1(i,1:3) = zero
1879 an(k) = zero
1880 ELSE
1881 an(k)=xmu(i)**2
1882 an(k)=one/an(k)
1883 ENDIF
1884 IF(xmu2(i)<=em30) THEN
1885 xmu2(i) = em30
1886 dir2(i,1:3) = zero
1887 bn(k) = zero
1888 ELSE
1889 bn(k)=xmu2(i)**2
1890 bn(k)=one/bn(k)
1891 ENDIF
1892 ENDDO
1893
1894
1895 ELSEIF(mfrot==4)THEN
1896C--- Exponential decay model
1897#include "vectorize.inc"
1898 DO k=1,nfisot ! isotropic friction couples
1899 i = indexisot(k)
1900 IF(pene(i) == 0) THEN
1901 xmu(i) = zero
1902 xmu2(i) = zero
1903 cycle
1904 ENDIF
1905 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1906 v2 = (vx(i) - n1(i)*aa)**2
1907 . + (vy(i) - n2(i)*aa)**2
1908 . + (vz(i) - n3(i)*aa)**2
1909 vv = sqrt(max(em30,v2))
1910 xmu(i) = fric_coefs(i,1)
1911 . + (fricc(i)-fric_coefs(i,1))*exp(-fric_coefs(i,2)*vv)
1912 xmu(i) = max(xmu(i),em30)
1913 ENDDO
1914c
1915#include "vectorize.inc"
1916 DO k=1,nforth ! Orthotropic friction couples
1917 i = indexorth(k)
1918 IF(pene(i) == 0) THEN
1919 xmu(i) = zero
1920 xmu2(i) = zero
1921 cycle
1922 ENDIF
1923c
1924 v2 = vx(i)*dir1(i,1) +vy(i)*dir1(i,2)+vz(i)*dir1(i,3)
1925 vv = max(em30,v2)
1926 xmu(i) = fricc(i)
1927 . + (fric_coefs(i,1)-fricc(i))*exp(-fric_coefs(i,2)*vv)
1928c
1929 v2 = vx(i)*dir2(i,1) +vy(i)*dir2(i,2)+vz(i)*dir2(i,3)
1930 vv = max(em30,v2)
1931 xmu2(i) = fric_coefs2(i,1)
1932 . + (fricc2(i)-fric_coefs2(i,1))*exp(-fric_coefs2(i,2)*vv)
1933 xmu(i) = max(xmu(i),em30)
1934 xmu2(i) = max(xmu2(i),em30)
1935 ENDDO
1936
1937#include "vectorize.inc"
1938 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
1939 i = indexorth(k)
1940 IF(xmu(i)<=em30) THEN
1941 xmu(i) = em30
1942 dir1(i,1:3) = zero
1943 an(k) = zero
1944 ELSE
1945 an(k)=xmu(i)**2
1946 an(k)=one/an(k)
1947 ENDIF
1948 IF(xmu2(i)<=em30) THEN
1949 xmu2(i) = em30
1950 dir2(i,1:3) = zero
1951 bn(k) = zero
1952 ELSE
1953 bn(k)=xmu2(i)**2
1954 bn(k)=one/bn(k)
1955 ENDIF
1956 ENDDO
1957
1958 ENDIF
1959
1960 ENDIF ! IORTHFRIC
1961C------------------
1962C TANGENT FORCE CALCULATION
1963C------------------
1964 fsav4 = zero
1965 fsav5 = zero
1966 fsav6 = zero
1967 fsav12= zero
1968 fsav13= zero
1969 fsav14= zero
1970 fsav15= zero
1971 IF (ifq /= 0) THEN
1972C---------------------------------
1973C INCREMENTAL (STIFFNESS) FORMULATION
1974C---------------------------------
1975 IF (ifq==13) THEN
1976 alpha = max(one,alpha0*dt12)
1977 ELSE
1978 alpha = alpha0
1979 ENDIF
1980
1981 IF(intnitsche > 0 ) THEN
1982
1983 DO i=1,jlt
1984 IF(pene(i) == zero)cycle
1985 ftn = forneqsi(i,1)*n1(i) + forneqsi(i,2)*n2(i) + forneqsi(i,3)*n3(i)
1986 fnitxt(i) = half*(forneqsi(i,1) - ftn*n1(i))
1987 fnityt(i) = half*(forneqsi(i,2) - ftn*n2(i))
1988 fnitzt(i) = half*(forneqsi(i,3) - ftn*n3(i))
1989 ENDDO
1990 ELSE
1991 DO i=1,jlt
1992 fnitxt(i) = zero
1993 fnityt(i) = zero
1994 fnitzt(i) = zero
1995 ENDDO
1996 ENDIF
1997
1998 IF(iorthfric ==0 ) THEN
1999C++ Isotropic Friction
2000
2001 IF (inconv==1) THEN
2002 DO i=1,jlt
2003 IF(pene(i) == zero)cycle
2004 fx = stif0(i)*vx(i)*dt12
2005 fy = stif0(i)*vy(i)*dt12
2006 fz = stif0(i)*vz(i)*dt12
2007 jg = nsvg(i)
2008 IF(jg>0) THEN
2009 n = cand_n_n(i)
2010c SECND_FR(4:6,N) is the old friction force
2011 fx = secnd_fr(4,n) + alpha*fx
2012 fy = secnd_fr(5,n) + alpha*fy
2013 fz = secnd_fr(6,n) + alpha*fz
2014 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2015 fx = fx - ftn*n1(i)
2016 fy = fy - ftn*n2(i)
2017 fz = fz - ftn*n3(i)
2018 fx = fx + fnitxt(i)
2019 fy = fy + fnityt(i)
2020 fz = fz + fnitzt(i)
2021 ft = fx*fx + fy*fy + fz*fz
2022 ft = max(ft,tiny)
2023 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2024 beta = min(one,xmu(i)*sqrt(fn/ft))
2025 fxt(i) = fx * beta
2026 fyt(i) = fy * beta
2027 fzt(i) = fz * beta
2028
2029#include "lockon.inc"
2030 IF (abs(fxt(i)- fnitxt(i))>abs(secnd_fr(1,n))) secnd_fr(1,n) = fxt(i)- fnitxt(i)
2031 IF (abs(fyt(i)- fnityt(i))>abs(secnd_fr(2,n))) secnd_fr(2,n) = fyt(i)- fnityt(i)
2032 IF (abs(fzt(i)- fnitzt(i))>abs(secnd_fr(3,n))) secnd_fr(3,n) = fzt(i)- fnitzt(i)
2033C case equal abs but opposite sign
2034 IF ((fxt(i)- fnitxt(i))==-secnd_fr(1,n)) secnd_fr(1,n) = abs(fxt(i)- fnitxt(i))
2035 IF ((fyt(i)- fnityt(i))==-secnd_fr(2,n)) secnd_fr(2,n) = abs(fyt(i)- fnityt(i))
2036 IF ((fzt(i)- fnitzt(i))==-secnd_fr(3,n)) secnd_fr(3,n) = abs(fzt(i)- fnitzt(i))
2037
2038#include "lockoff.inc"
2039 ELSE ! cas noeud remote en SPMD
2040 jg = -jg
2041c SECND_FRFI(NIN)%P(4:6,JG) is the old friction force
2042 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2043 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2044 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2045 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2046 fx = fx - ftn*n1(i)
2047 fy = fy - ftn*n2(i)
2048 fz = fz - ftn*n3(i)
2049 fx = fx + fnitxt(i)
2050 fy = fy + fnityt(i)
2051 fz = fz + fnitzt(i)
2052 ft = fx*fx + fy*fy + fz*fz
2053 ft = max(ft,tiny)
2054 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2055 beta = min(one,xmu(i)*sqrt(fn/ft))
2056 fxt(i) = fx * beta
2057 fyt(i) = fy * beta
2058 fzt(i) = fz * beta
2059#include "lockon.inc"
2060 IF ( abs(fxt(i)- fnitxt(i)) > abs(secnd_frfi(nin)%P(1,jg)) )
2061 * secnd_frfi(nin)%P(1,jg) = fxt(i)- fnitxt(i)
2062 IF ( abs(fyt(i)- fnityt(i)) > abs(secnd_frfi(nin)%P(2,jg)) )
2063 * secnd_frfi(nin)%P(2,jg) = fyt(i)- fnityt(i)
2064 IF ( abs(fzt(i)- fnitzt(i)) > abs(secnd_frfi(nin)%P(3,jg)) )
2065 * secnd_frfi(nin)%P(3,jg) = fzt(i)- fnitzt(i)
2066C
2067 IF ((fxt(i)- fnitxt(i))== - secnd_frfi(nin)%P(1,jg) )
2068 * secnd_frfi(nin)%P(1,jg) = abs(fxt(i)- fnitxt(i))
2069 IF ((fyt(i)- fnityt(i))== - secnd_frfi(nin)%P(2,jg) )
2070 * secnd_frfi(nin)%P(2,jg) = abs(fyt(i)- fnityt(i))
2071 IF ((fzt(i)- fnitzt(i))== - secnd_frfi(nin)%P(3,jg) )
2072 * secnd_frfi(nin)%P(3,jg) = abs(fzt(i)- fnitzt(i))
2073#include "lockoff.inc"
2074 ENDIF
2075C------- total force
2076 fxi(i)=fxi(i) + fxt(i)
2077 fyi(i)=fyi(i) + fyt(i)
2078 fzi(i)=fzi(i) + fzt(i)
2079 efric_l(i) = dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
2080 econvt = econvt + efric_l(i)
2081 ENDDO
2082C--------implicit non converge---
2083 ELSE
2084 DO i=1,jlt
2085 IF(pene(i) == zero)cycle
2086 fx = stif0(i)*vx(i)*dt12
2087 fy = stif0(i)*vy(i)*dt12
2088 fz = stif0(i)*vz(i)*dt12
2089 jg = nsvg(i)
2090 n = cand_n_n(i)
2091 IF(jg>0) THEN
2092c SECND_FR(4:6,N) is the old friction force
2093 fx = secnd_fr(4,n) + alpha*fx
2094 fy = secnd_fr(5,n) + alpha*fy
2095 fz = secnd_fr(6,n) + alpha*fz
2096 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2097 fx = fx - ftn*n1(i)
2098 fy = fy - ftn*n2(i)
2099 fz = fz - ftn*n3(i)
2100 fx = fx + fnitxt(i)
2101 fy = fy + fnityt(i)
2102 fz = fz + fnitzt(i)
2103 ft = fx*fx + fy*fy + fz*fz
2104 ft = max(ft,tiny)
2105 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2106 beta = min(one,xmu(i)*sqrt(fn/ft))
2107 fxt(i) = fx * beta
2108 fyt(i) = fy * beta
2109 fzt(i) = fz * beta
2110 fxi(i)=fxi(i) + fxt(i)
2111 fyi(i)=fyi(i) + fyt(i)
2112 fzi(i)=fzi(i) + fzt(i)
2113 ELSE ! cas noeud remote en SPMD
2114 jg = -jg
2115 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2116 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2117 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2118 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2119 fx = fx - ftn*n1(i)
2120 fy = fy - ftn*n2(i)
2121 fz = fz - ftn*n3(i)
2122 fx = fx + fnitxt(i)
2123 fy = fy + fnityt(i)
2124 fz = fz + fnitzt(i)
2125 ft = fx*fx + fy*fy + fz*fz
2126 ft = max(ft,tiny)
2127 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2128 beta = min(one,xmu(i)*sqrt(fn/ft))
2129 fxt(i) = fx * beta
2130 fyt(i) = fy * beta
2131 fzt(i) = fz * beta
2132 fxi(i)=fxi(i) + fxt(i)
2133 fyi(i)=fyi(i) + fyt(i)
2134 fzi(i)=fzi(i) + fzt(i)
2135 ENDIF
2136C IFPEN(INDEX(I)) = 1
2137 ENDDO
2138 ENDIF
2139
2140 ELSE
2141C++ Orthotropic Friction
2142
2143 IF (inconv==1) THEN
2144#include "vectorize.inc"
2145 DO k=1,nfisot ! isotropic friction couples
2146 i = indexisot(k)
2147 IF(pene(i) == zero)cycle
2148 fx = stif0(i)*vx(i)*dt12
2149 fy = stif0(i)*vy(i)*dt12
2150 fz = stif0(i)*vz(i)*dt12
2151 jg = nsvg(i)
2152 IF(jg>0) THEN
2153 n = cand_n_n(i)
2154c SECND_FR(4:6,N) is the old friction force
2155 fx = secnd_fr(4,n) + alpha*fx
2156 fy = secnd_fr(5,n) + alpha*fy
2157 fz = secnd_fr(6,n) + alpha*fz
2158 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2159 fx = fx - ftn*n1(i)
2160 fy = fy - ftn*n2(i)
2161 fz = fz - ftn*n3(i)
2162 fx = fx + fnitxt(i)
2163 fy = fy + fnityt(i)
2164 fz = fz + fnitzt(i)
2165 ft = fx*fx + fy*fy + fz*fz
2166 ft = max(ft,tiny)
2167 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2168 beta = min(one,xmu(i)*sqrt(fn/ft))
2169 fxt(i) = fx * beta
2170 fyt(i) = fy * beta
2171 fzt(i) = fz * beta
2172#include "lockon.inc"
2173 IF (abs(fxt(i)- fnitxt(i))>abs(secnd_fr(1,n))) secnd_fr(1,n) = fxt(i)- fnitxt(i)
2174 IF (abs(fyt(i)- fnityt(i))>abs(secnd_fr(2,n))) secnd_fr(2,n) = fyt(i)- fnityt(i)
2175 IF (abs(fzt(i)- fnitzt(i))>abs(secnd_fr(3,n))) secnd_fr(3,n) = fzt(i)- fnitzt(i)
2176C case equal abs but opposite sign
2177 IF ((fxt(i)- fnitxt(i))==-secnd_fr(1,n)) secnd_fr(1,n) = abs(fxt(i)- fnitxt(i))
2178 IF ((fyt(i)- fnityt(i))==-secnd_fr(2,n)) secnd_fr(2,n) = abs(fyt(i)- fnityt(i))
2179 IF ((fzt(i)- fnitzt(i))==-secnd_fr(3,n)) secnd_fr(3,n) = abs(fzt(i)- fnitzt(i))
2180#include "lockoff.inc"
2181 ELSE ! cas noeud remote en SPMD
2182 jg = -jg
2183c SECND_FRFI(NIN)%P(4:6,JG) is the old friction force
2184 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2185 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2186 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2187 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2188 fx = fx - ftn*n1(i)
2189 fy = fy - ftn*n2(i)
2190 fz = fz - ftn*n3(i)
2191 fx = fx + fnitxt(i)
2192 fy = fy + fnityt(i)
2193 fz = fz + fnitzt(i)
2194 ft = fx*fx + fy*fy + fz*fz
2195 ft = max(ft,tiny)
2196 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2197 beta = min(one,xmu(i)*sqrt(fn/ft))
2198 fxt(i) = fx * beta
2199 fyt(i) = fy * beta
2200 fzt(i) = fz * beta
2201#include "lockon.inc"
2202 IF ( abs(fxt(i)- fnitxt(i)) > abs(secnd_frfi(nin)%P(1,jg)) )
2203 * secnd_frfi(nin)%P(1,jg) = fxt(i)- fnitxt(i)
2204 IF ( abs(fyt(i)- fnityt(i)) > abs(secnd_frfi(nin)%P(2,jg)) )
2205 * secnd_frfi(nin)%P(2,jg) = fyt(i)- fnityt(i)
2206 IF ( abs(fzt(i)- fnitzt(i)) > abs(secnd_frfi(nin)%P(3,jg)) )
2207 * secnd_frfi(nin)%P(3,jg) = fzt(i)- fnitzt(i)
2208C
2209 IF ((fxt(i)- fnitxt(i))== - secnd_frfi(nin)%P(1,jg) )
2210 * secnd_frfi(nin)%P(1,jg) = abs(fxt(i)- fnitxt(i))
2211 IF ((fyt(i)- fnityt(i))== - secnd_frfi(nin)%P(2,jg) )
2212 * secnd_frfi(nin)%P(2,jg) = abs(fyt(i)- fnityt(i))
2213 IF ((fzt(i)- fnitzt(i))== - secnd_frfi(nin)%P(3,jg) )
2214 * secnd_frfi(nin)%P(3,jg) = abs(fzt(i)- fnitzt(i))
2215#include "lockoff.inc"
2216 ENDIF
2217C------- total force
2218 fxi(i)=fxi(i) + fxt(i)
2219 fyi(i)=fyi(i) + fyt(i)
2220 fzi(i)=fzi(i) + fzt(i)
2221 efric_l(i) = dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
2222 econvt = econvt + efric_l(i)
2223 ENDDO
2224
2225#include "vectorize.inc"
2226 DO k=1,nforth ! Orthotropic friction couples
2227 i = indexorth(k)
2228 IF(pene(i) == zero)cycle
2229 fx = stif0(i)*vx(i)*dt12
2230 fy = stif0(i)*vy(i)*dt12
2231 fz = stif0(i)*vz(i)*dt12
2232 jg = nsvg(i)
2233 IF(jg>0) THEN
2234 n = cand_n_n(i)
2235c SECND_FR(4:6,N) is the old friction force
2236 fx = secnd_fr(4,n) + alpha*fx
2237 fy = secnd_fr(5,n) + alpha*fy
2238 fz = secnd_fr(6,n) + alpha*fz
2239 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2240 fx = fx - ftn*n1(i)
2241 fy = fy - ftn*n2(i)
2242 fz = fz - ftn*n3(i)
2243
2244 fx = fx + fnitxt(i)
2245 fy = fy + fnityt(i)
2246 fz = fz + fnitzt(i)
2247
2248 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2249 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2250 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2251 ft = max(ft,em30)
2252 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2253
2254 beta = fn/ft
2255
2256 IF(beta == 0 ) THEN
2257 fxt(i) = zero
2258 fyt(i) = zero
2259 fzt(i) = zero
2260 ELSEIF(beta > 1) THEN ! Inside the ellipse
2261 fxt(i) = fx
2262 fyt(i) = fy
2263 fzt(i) = fz
2264 ELSE ! outside the ellipse
2265
2266! Projection on local tangent of ellipse (outside of ellipse)
2267! ANS = (Fadh-Fproj).n
2268 nep1 =ftt1*an(k)/fn
2269 nep2 =ftt2*bn(k)/fn
2270 nep =nep1*nep1+nep2*nep2
2271 nep =sqrt(nep)
2272 ep=nep1*ftt1+nep2*ftt2
2273
2274 ans=(ep-sqrt(ep))/max(em20,nep)
2275 nep1 =nep1/max(em20,nep)
2276 nep2 =nep2/max(em20,nep)
2277! Projection on ellipse
2278 c11 =ftt1-ans*nep1
2279 c22 =ftt2-ans*nep2
2280 alphaf = atan(c22/c11)
2281 signc = ftt1/max(em20,abs(ftt1))
2282 csa = signc*abs(cos(alphaf))
2283 signc = ftt2/max(em20,abs(ftt2))
2284 sna = signc*abs(sin(alphaf))
2285! Ft computation
2286 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2287 ftt1 = ft * csa
2288 ftt2 = ft * sna
2289
2290 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2291 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2292 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2293
2294 ENDIF
2295
2296#include "lockon.inc"
2297 IF (abs(fxt(i))>abs(secnd_fr(1,n))) secnd_fr(1,n) = fxt(i)
2298 IF (abs(fyt(i))>abs(secnd_fr(2,n))) secnd_fr(2,n) = fyt(i)
2299 IF (abs(fzt(i))>abs(secnd_fr(3,n))) secnd_fr(3,n) = fzt(i)
2300C case equal abs but opposite sign
2301 IF (fxt(i)==-secnd_fr(1,n)) secnd_fr(1,n) = abs(fxt(i))
2302 IF (fyt(i)==-secnd_fr(2,n)) secnd_fr(2,n) = abs(fyt(i))
2303 IF (fzt(i)==-secnd_fr(3,n)) secnd_fr(3,n) = abs(fzt(i))
2304#include "lockoff.inc"
2305 ELSE ! cas noeud remote en SPMD
2306 jg = -jg
2307c SECND_FRFI(NIN)%P(4:6,JG) is the old friction force
2308 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2309 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2310 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2311 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2312 fx = fx - ftn*n1(i)
2313 fy = fy - ftn*n2(i)
2314 fz = fz - ftn*n3(i)
2315 fx = fx + fnitxt(i)
2316 fy = fy + fnityt(i)
2317 fz = fz + fnitzt(i)
2318 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2319 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2320 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2321 ft = max(ft,em30)
2322 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2323
2324 beta = fn/ft
2325
2326 IF(beta == 0 ) THEN
2327 fxt(i) = zero
2328 fyt(i) = zero
2329 fzt(i) = zero
2330 ELSEIF(beta > 1) THEN ! Inside the ellipse
2331 fxt(i) = fx
2332 fyt(i) = fy
2333 fzt(i) = fz
2334 ELSE ! outside the ellipse
2335
2336! Projection on local tangent of ellipse (outside of ellipse)
2337! ANS = (Fadh-Fproj).n
2338 nep1 =ftt1*an(k)/fn
2339 nep2 =ftt2*bn(k)/fn
2340 nep =nep1*nep1+nep2*nep2
2341 nep =sqrt(nep)
2342
2343 ep=nep1*ftt1+nep2*ftt2
2344
2345 ans=(ep-sqrt(ep))/max(em20,nep)
2346 nep1 =nep1/max(em20,nep)
2347 nep2 =nep2/max(em20,nep)
2348
2349! Projection on ellipse
2350 c11 =ftt1-ans*nep1
2351 c22 =ftt2-ans*nep2
2352
2353 alphaf = atan(c22/c11)
2354
2355 signc = ftt1/max(em20,abs(ftt1))
2356 csa = signc*abs(cos(alphaf))
2357 signc = ftt2/max(em20,abs(ftt2))
2358 sna = signc*abs(sin(alphaf))
2359! Ft computation
2360 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2361 ftt1 = ft * csa
2362 ftt2 = ft * sna
2363
2364 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2365 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2366 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2367
2368 ENDIF
2369#include "lockon.inc"
2370 IF ( abs(fxt(i)) > abs(secnd_frfi(nin)%P(1,jg)) )
2371 * secnd_frfi(nin)%P(1,jg) = fxt(i)
2372 IF ( abs(fyt(i)) > abs(secnd_frfi(nin)%P(2,jg)) )
2373 * secnd_frfi(nin)%P(2,jg) = fyt(i)
2374 IF ( abs(fzt(i)) > abs(secnd_frfi(nin)%P(3,jg)) )
2375 * secnd_frfi(nin)%P(3,jg) = fzt(i)
2376C
2377 IF (fxt(i)== - secnd_frfi(nin)%P(1,jg) )
2378 * secnd_frfi(nin)%P(1,jg) = abs(fxt(i))
2379 IF (fyt(i)== - secnd_frfi(nin)%P(2,jg) )
2380 * secnd_frfi(nin)%P(2,jg) = abs(fyt(i))
2381 IF (fzt(i)== - secnd_frfi(nin)%P(3,jg) )
2382 * secnd_frfi(nin)%P(3,jg) = abs(fzt(i))
2383#include "lockoff.inc"
2384 ENDIF
2385C------- total force
2386 fxi(i)=fxi(i) + fxt(i)
2387 fyi(i)=fyi(i) + fyt(i)
2388 fzi(i)=fzi(i) + fzt(i)
2389 efric_l(i) = dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
2390 econvt = econvt + efric_l(i)
2391 ENDDO
2392
2393
2394C--------implicit non converge---
2395 ELSE
2396#include "vectorize.inc"
2397 DO k=1,nfisot ! isotropic friction couples
2398 i = indexisot(k)
2399 IF(pene(i) == zero)cycle
2400 fx = stif0(i)*vx(i)*dt12
2401 fy = stif0(i)*vy(i)*dt12
2402 fz = stif0(i)*vz(i)*dt12
2403 jg = nsvg(i)
2404 n = cand_n_n(i)
2405 IF(jg>0) THEN
2406c SECND_FR(4:6,N) is the old friction force
2407 fx = secnd_fr(4,n) + alpha*fx
2408 fy = secnd_fr(5,n) + alpha*fy
2409 fz = secnd_fr(6,n) + alpha*fz
2410 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2411 fx = fx - ftn*n1(i)
2412 fy = fy - ftn*n2(i)
2413 fz = fz - ftn*n3(i)
2414 fx = fx + fnitxt(i)
2415 fy = fy + fnityt(i)
2416 fz = fz + fnitzt(i)
2417 ft = fx*fx + fy*fy + fz*fz
2418 ft = max(ft,tiny)
2419 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2420 beta = min(one,xmu(i)*sqrt(fn/ft))
2421 fxt(i) = fx * beta
2422 fyt(i) = fy * beta
2423 fzt(i) = fz * beta
2424 fxi(i)=fxi(i) + fxt(i)
2425 fyi(i)=fyi(i) + fyt(i)
2426 fzi(i)=fzi(i) + fzt(i)
2427 ELSE ! cas noeud remote en SPMD
2428 jg = -jg
2429 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2430 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2431 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2432 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2433 fx = fx - ftn*n1(i)
2434 fy = fy - ftn*n2(i)
2435 fz = fz - ftn*n3(i)
2436 fx = fx + fnitxt(i)
2437 fy = fy + fnityt(i)
2438 fz = fz + fnitzt(i)
2439 ft = fx*fx + fy*fy + fz*fz
2440 ft = max(ft,tiny)
2441 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2442 beta = min(one,xmu(i)*sqrt(fn/ft))
2443 fxt(i) = fx * beta
2444 fyt(i) = fy * beta
2445 fzt(i) = fz * beta
2446 fxi(i)=fxi(i) + fxt(i)
2447 fyi(i)=fyi(i) + fyt(i)
2448 fzi(i)=fzi(i) + fzt(i)
2449 ENDIF
2450C IFPEN(INDEX(I)) = 1
2451 ENDDO
2452
2453#include "vectorize.inc"
2454 DO k=1,nforth ! Orthotropic friction couples
2455 i = indexorth(k)
2456 IF(pene(i) == zero)cycle
2457 fx = stif0(i)*vx(i)*dt12
2458 fy = stif0(i)*vy(i)*dt12
2459 fz = stif0(i)*vz(i)*dt12
2460 jg = nsvg(i)
2461 n = cand_n_n(i)
2462 IF(jg>0) THEN
2463c SECND_FR(4:6,N) is the old friction force
2464 fx = secnd_fr(4,n) + alpha*fx
2465 fy = secnd_fr(5,n) + alpha*fy
2466 fz = secnd_fr(6,n) + alpha*fz
2467 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2468 fx = fx - ftn*n1(i)
2469 fy = fy - ftn*n2(i)
2470 fz = fz - ftn*n3(i)
2471 fx = fx + fnitxt(i)
2472 fy = fy + fnityt(i)
2473 fz = fz + fnitzt(i)
2474
2475 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2476 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2477 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2478 ft = max(ft,em30)
2479 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2480
2481 beta = fn/ft
2482
2483 IF(beta == 0 ) THEN
2484 fxt(i) = zero
2485 fyt(i) = zero
2486 fzt(i) = zero
2487 ELSEIF(beta > 1) THEN ! Inside the ellipse
2488 fxt(i) = fx
2489 fyt(i) = fy
2490 fzt(i) = fz
2491 ELSE ! outside the ellipse
2492
2493! Projection on local tangent of ellipse (outside of ellipse)
2494! ANS = (Fadh-Fproj).n
2495 nep1 =ftt1*an(k)/fn
2496 nep2 =ftt2*bn(k)/fn
2497 nep =nep1*nep1+nep2*nep2
2498 nep =sqrt(nep)
2499
2500 ep=nep1*ftt1+nep2*ftt2
2501
2502 ans=(ep-sqrt(ep))/max(em20,nep)
2503 nep1 =nep1/max(em20,nep)
2504 nep2 =nep2/max(em20,nep)
2505
2506! Projection on ellipse
2507 c11 =ftt1-ans*nep1
2508 c22 =ftt2-ans*nep2
2509
2510 alphaf = atan(c22/c11)
2511
2512 signc = ftt1/max(em20,abs(ftt1))
2513 csa = signc*abs(cos(alphaf))
2514 signc = ftt2/max(em20,abs(ftt2))
2515 sna = signc*abs(sin(alphaf))
2516! Ft computation
2517 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2518 ftt1 = ft * csa
2519 ftt2 = ft * sna
2520
2521 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2522 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2523 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2524
2525 ENDIF
2526 fxi(i)=fxi(i) + fxt(i)
2527 fyi(i)=fyi(i) + fyt(i)
2528 fzi(i)=fzi(i) + fzt(i)
2529 ELSE ! cas noeud remote en SPMD
2530 jg = -jg
2531 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2532 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2533 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2534 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2535 fx = fx - ftn*n1(i)
2536 fy = fy - ftn*n2(i)
2537 fz = fz - ftn*n3(i)
2538 fx = fx + fnitxt(i)
2539 fy = fy + fnityt(i)
2540 fz = fz + fnitzt(i)
2541
2542 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2543 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2544 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2545 ft = max(ft,em30)
2546 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2547
2548 beta = fn/ft
2549
2550 IF(beta == 0 ) THEN
2551 fxt(i) = zero
2552 fyt(i) = zero
2553 fzt(i) = zero
2554 ELSEIF(beta > 1) THEN ! Inside the ellipse
2555 fxt(i) = fx
2556 fyt(i) = fy
2557 fzt(i) = fz
2558 ELSE ! outside the ellipse
2559
2560! projection on local tangent of ellipse(outside of ellipse)
2561! ANS = (Fadh-Fproj).n
2562 nep1 =ftt1*an(k)/fn
2563 nep2 =ftt2*bn(k)/fn
2564 nep =nep1*nep1+nep2*nep2
2565 nep =sqrt(nep)
2566
2567 ep=nep1*ftt1+nep2*ftt2
2568
2569 ans=(ep-sqrt(ep))/max(em20,nep)
2570 nep1 =nep1/max(em20,nep)
2571 nep2 =nep2/max(em20,nep)
2572
2573! projection on ellipse
2574 c11 =ftt1-ans*nep1
2575 c22 =ftt2-ans*nep2
2576
2577 alphaf = atan(c22/c11)
2578
2579 signc = ftt1/max(em20,abs(ftt1))
2580 csa = signc*abs(cos(alphaf))
2581 signc = ftt2/max(em20,abs(ftt2))
2582 sna = signc*abs(sin(alphaf))
2583! Ft computation
2584 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2585 ftt1 = ft * csa
2586 ftt2 = ft * sna
2587
2588 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2589 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2590 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2591
2592 ENDIF
2593
2594 fxi(i)=fxi(i) + fxt(i)
2595 fyi(i)=fyi(i) + fyt(i)
2596 fzi(i)=fzi(i) + fzt(i)
2597 ENDIF
2598C IFPEN(INDEX(I)) = 1
2599 ENDDO
2600
2601 ENDIF !imp
2602
2603 ENDIF ! iorth
2604
2605
2606
2607 ENDIF !ifq
2608C
2609C---------------------------------
2610 IF((anim_v(12)+outp_v(12)+h3d_data%N_VECT_PCONT>0.AND.
2611 . ((tt>=output%TANIM.AND.tt<=output%TANIM_STOP).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
2612 . (manim>=4.AND.manim<=15).OR.h3d_data%MH3D/=0))
2613 . .OR.h3d_data%N_VECT_PCONT_MAX>0)THEN
2614 IF (inconv==1) THEN
2615#include "lockon.inc"
2616 DO i=1,jlt
2617 IF(pene(i) == zero)cycle
2618 ftcont(1,ix1(i)) =ftcont(1,ix1(i)) + fxt(i)*h1(i)
2619 ftcont(2,ix1(i)) =ftcont(2,ix1(i)) + fyt(i)*h1(i)
2620 ftcont(3,ix1(i)) =ftcont(3,ix1(i)) + fzt(i)*h1(i)
2621 ftcont(1,ix2(i)) =ftcont(1,ix2(i)) + fxt(i)*h2(i)
2622 ftcont(2,ix2(i)) =ftcont(2,ix2(i)) + fyt(i)*h2(i)
2623 ftcont(3,ix2(i)) =ftcont(3,ix2(i)) + fzt(i)*h2(i)
2624 ftcont(1,ix3(i)) =ftcont(1,ix3(i)) + fxt(i)*h3(i)
2625 ftcont(2,ix3(i)) =ftcont(2,ix3(i)) + fyt(i)*h3(i)
2626 ftcont(3,ix3(i)) =ftcont(3,ix3(i)) + fzt(i)*h3(i)
2627 ftcont(1,ix4(i)) =ftcont(1,ix4(i)) + fxt(i)*h4(i)
2628 ftcont(2,ix4(i)) =ftcont(2,ix4(i)) + fyt(i)*h4(i)
2629 ftcont(3,ix4(i)) =ftcont(3,ix4(i)) + fzt(i)*h4(i)
2630 jg = nsvg(i)
2631 IF(jg>0) THEN
2632C In SPMD: Treatment to be redone after reception node Remote if JG <0
2633 IF (jg >numnod) THEN
2634 ig = jg - numnod
2635 CALL i24for1_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
2636 + ig ,fxt(i) ,fyt(i) ,fzt(i) ,ftcont ,
2637 + -1 )
2638 ELSE
2639 ftcont(1,jg)=ftcont(1,jg)- fxt(i)
2640 ftcont(2,jg)=ftcont(2,jg)- fyt(i)
2641 ftcont(3,jg)=ftcont(3,jg)- fzt(i)
2642 ENDIF
2643 ELSE
2644 jg=-jg
2645 IF(isedge_fi(nin)%P(jg)==1)THEN
2646 CALL i24for1_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
2647 + jg ,fxt(i) ,fyt(i) ,fzt(i) ,ftconti(nin)%P(1,1) ,
2648 + -1 )
2649 ELSE ! cas noeud remote en SPMD non edge
2650 ftconti(nin)%P(1,jg)=ftconti(nin)%P(1,jg)-fxt(i)
2651 ftconti(nin)%P(2,jg)=ftconti(nin)%P(2,jg)-fyt(i)
2652 ftconti(nin)%P(3,jg)=ftconti(nin)%P(3,jg)-fzt(i)
2653 ENDIF
2654 END IF
2655 ENDDO
2656#include "lockoff.inc"
2657 END IF !(INCONV==1) THEN
2658 ENDIF
2659C
2660C---------------------------------
2661 fsav22= zero
2662 fsav23= zero
2663 fsav24= zero
2664 fsav29= zero
2665 DO i=1,jlt
2666 IF(pene(i) == zero)cycle
2667 impx=fxt(i)*dt12
2668 impy=fyt(i)*dt12
2669 impz=fzt(i)*dt12
2670 fsav4 =fsav4 +impx
2671 fsav5 =fsav5 +impy
2672 fsav6 =fsav6 +impz
2673 impx=fxi(i)*dt12
2674 impy=fyi(i)*dt12
2675 impz=fzi(i)*dt12
2676 fsav12=fsav12+abs(impx)
2677 fsav13=fsav13+abs(impy)
2678 fsav14=fsav14+abs(impz)
2679 fsav15=fsav15+sqrt(impx*impx+impy*impy+impz*impz)
2680 xp(i)=xi(i)+pene(i)*n1(i)
2681 yp(i)=yi(i)+pene(i)*n2(i)
2682 zp(i)=zi(i)+pene(i)*n3(i)
2683 fsav22=fsav22+yp(i)*impz-zp(i)*impy
2684 fsav23=fsav23+zp(i)*impx-xp(i)*impz
2685 fsav24=fsav24+xp(i)*impy-yp(i)*impx
2686 ENDDO
2687
2688 IF(intcarea > 0) THEN
2689 DO i=1,jlt
2690 IF(pene(i) == zero)cycle
2691 jg = nsvg(i)
2692 IF(jg>0) THEN ! Only local nodes in i25for3, remote nodes will be taken into account later
2693 IF(jg <= numnod) THEN
2694 fsav29=fsav29+parameters%INTAREAN(jg)
2695 ELSE
2696 ig = jg - numnod
2697 CALL i24intarea_fic(irtse ,nsne ,is2se ,is2pt ,ig ,
2698 + nrtse , numnod , parameters%INTAREAN, arean_fic)
2699 fsav29=fsav29 + arean_fic
2700 ENDIF
2701 ELSE
2702 jg=-jg
2703 IF(isedge_fi(nin)%P(jg)==1)THEN
2704 CALL i24intarea_fic(irtse_fi(nin)%P ,nsnr ,is2se_fi(nin)%P ,is2pt_fi(nin)%P ,jg ,
2705 + nsnr , nsnr ,intareanfi(nin)%P, arean_fic)
2706 fsav29=fsav29 + arean_fic
2707 ELSE ! cas noeud remote en SPMD non edge
2708 fsav29=fsav29 + intareanfi(nin)%P(jg)
2709 ENDIF
2710 ENDIF
2711 ENDDO
2712 ENDIF
2713#include "lockon.inc"
2714 fsav(4) = fsav(4) + fsav4
2715 fsav(5) = fsav(5) + fsav5
2716 fsav(6) = fsav(6) + fsav6
2717 fsav(12) = fsav(12) + fsav12
2718 fsav(13) = fsav(13) + fsav13
2719 fsav(14) = fsav(14) + fsav14
2720 fsav(15) = fsav(15) + fsav15
2721 fsav(22) = fsav(22) + fsav22
2722 fsav(23) = fsav(23) + fsav23
2723 fsav(24) = fsav(24) + fsav24
2724 fsav(26) = fsav(26) + econtt
2725 fsav(27) = fsav(27) + econvt
2726 fsav(28) = fsav(28) + econtdt
2727 fsav(29) = fsav(29) + fsav29
2728#include "lockoff.inc"
2729C
2730 IF(isensint(1)/=0) THEN
2731 DO i=1,jlt
2732 fsavparit(1,4,i+nft) = fxi(i)
2733 fsavparit(1,5,i+nft) = fyi(i)
2734 fsavparit(1,6,i+nft) = fzi(i)
2735 ENDDO
2736 ENDIF
2737C---------------------------------
2738C SORTIES TH PAR SOUS INTERFACE
2739C---------------------------------
2740 IF(nisub/=0)THEN
2741 DO i=1,jlt
2742 IF(pene(i) == zero)cycle
2743 nn = nsvg(i)
2744 IF(nn>0)THEN
2745 in=cn_loc(i)
2746 IF (msegtyp(ce_loc(i)) < 0) THEN
2747 ie= - msegtyp(ce_loc(i))
2748 ELSE
2749 ie = ce_loc(i)
2750 ENDIF
2751 jj =addsubs(in)
2752 kk =addsubm(ie)
2753 DO WHILE(jj<addsubs(in+1))
2754 jsub=lisubs(jj)
2755 itypsub = typsub(jsub)
2756
2757 IF(itypsub == 1 ) THEN
2758 iss1 = bitget(inflg_subs(jj),0)
2759 iss2 = bitget(inflg_subs(jj),1)
2760 igrn = bitget(inflg_subs(jj),2)
2761 ksub=lisubm(kk)
2762 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
2763 ims1 = bitget(inflg_subm(kk),0)
2764 ims2 = bitget(inflg_subm(kk),1)
2765 IF(ksub==jsub)THEN
2766! S and M candidates on the same sub_interface
2767 IF(.NOT.(((ims1 == 1 .OR. ims2 == 1) .AND. igrn == 1).OR.
2768 . (ims1 == 1 .AND. iss2 == 1).OR.
2769 . (ims2 == 1 .AND. iss1 == 1))) THEN
2770 kk=kk+1
2771 ksub=lisubm(kk)
2772 cycle
2773 END IF
2774 impx=fxt(i)*dt12
2775 impy=fyt(i)*dt12
2776 impz=fzt(i)*dt12
2777C main side :
2778 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
2779 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
2780 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
2781C
2782 impx=fxi(i)*dt12
2783 impy=fyi(i)*dt12
2784 impz=fzi(i)*dt12
2785 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
2786 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
2787 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
2788C
2789 IF(isensint(jsub+1)/=0) THEN
2790 fsavparit(jsub+1,4,i+nft) = fxt(i)
2791 fsavparit(jsub+1,5,i+nft) = fyt(i)
2792 fsavparit(jsub+1,6,i+nft) = fzt(i)
2793 ENDIF
2794C
2795 fsavsub1(15,jsub)= fsavsub1(15,jsub)
2796 . +sqrt(impx*impx+impy*impy+impz*impz)
2797 fsavsub1(22,jsub)=fsavsub1(22,jsub)
2798 . +yp(i)*impz-zp(i)*impy
2799 fsavsub1(23,jsub)=fsavsub1(23,jsub)
2800 . +zp(i)*impx-xp(i)*impz
2801 fsavsub1(24,jsub)=fsavsub1(24,jsub)
2802 . +xp(i)*impy-yp(i)*impx
2803C
2804 END IF
2805
2806 kk=kk+1
2807 ksub=lisubm(kk)
2808 ENDDO
2809 jj=jj+1
2810
2811 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface : second side
2812
2813 impx=fxt(i)*dt12
2814 impy=fyt(i)*dt12
2815 impz=fzt(i)*dt12
2816C main side :
2817 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
2818 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
2819 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
2820C
2821 impx=fxi(i)*dt12
2822 impy=fyi(i)*dt12
2823 impz=fzi(i)*dt12
2824 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
2825 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
2826 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
2827C
2828 IF(isensint(jsub+1)/=0) THEN
2829 fsavparit(jsub+1,4,i+nft) = fxt(i)
2830 fsavparit(jsub+1,5,i+nft) = fyt(i)
2831 fsavparit(jsub+1,6,i+nft) = fzt(i)
2832 ENDIF
2833C
2834 fsavsub1(15,jsub)= fsavsub1(15,jsub)
2835 . +sqrt(impx*impx+impy*impy+impz*impz)
2836 fsavsub1(22,jsub)=fsavsub1(22,jsub)
2837 . +yp(i)*impz-zp(i)*impy
2838 fsavsub1(23,jsub)=fsavsub1(23,jsub)
2839 . +zp(i)*impx-xp(i)*impz
2840 fsavsub1(24,jsub)=fsavsub1(24,jsub)
2841 . +xp(i)*impy-yp(i)*impx
2842
2843 jj=jj+1
2844
2845
2846 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfaces
2847C
2848C Find if node is on Surface S1 -> S2 or S2 ->S1
2849 iss2 = bitget(inflg_subs(jj),0)
2850 iss1 = bitget(inflg_subs(jj),1)
2851 ksub=lisubm(kk)
2852 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
2853 ims2 = bitget(inflg_subm(kk),0)
2854 ims1 = bitget(inflg_subm(kk),1)
2855 IF(ksub==jsub)THEN
2856! S and M candidates on the same sub_interface
2857 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
2858 . (ims2 == 1 .AND. iss1 == 1))) THEN
2859 kk=kk+1
2860 ksub=lisubm(kk)
2861 cycle
2862 END IF
2863
2864 impx=fxt(i)*dt12
2865 impy=fyt(i)*dt12
2866 impz=fzt(i)*dt12
2867 IF(ims2 > 0) THEN
2868C main side :
2869 fsavsub1(4,jsub)=fsavsub1(4,jsub)-impx
2870 fsavsub1(5,jsub)=fsavsub1(5,jsub)-impy
2871 fsavsub1(6,jsub)=fsavsub1(6,jsub)-impz
2872 ELSE
2873 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
2874 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
2875 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
2876 ENDIF
2877C
2878 impx=fxi(i)*dt12
2879 impy=fyi(i)*dt12
2880 impz=fzi(i)*dt12
2881 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
2882 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
2883 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
2884C
2885 IF(isensint(jsub+1)/=0) THEN
2886 IF(ims2 > 0) THEN
2887 fsavparit(jsub+1,4,i+nft) = -fxt(i)
2888 fsavparit(jsub+1,5,i+nft) = -fyt(i)
2889 fsavparit(jsub+1,6,i+nft) = -fzt(i)
2890 ELSE
2891 fsavparit(jsub+1,4,i+nft) = fxt(i)
2892 fsavparit(jsub+1,5,i+nft) = fyt(i)
2893 fsavparit(jsub+1,6,i+nft) = fzt(i)
2894 ENDIF
2895 ENDIF
2896C
2897 fsavsub1(15,jsub)= fsavsub1(15,jsub)
2898 . +sqrt(impx*impx+impy*impy+impz*impz)
2899 fsavsub1(22,jsub)=fsavsub1(22,jsub)
2900 . +yp(i)*impz-zp(i)*impy
2901 fsavsub1(23,jsub)=fsavsub1(23,jsub)
2902 . +zp(i)*impx-xp(i)*impz
2903 fsavsub1(24,jsub)=fsavsub1(24,jsub)
2904 . +xp(i)*impy-yp(i)*impx
2905C
2906 ENDIF
2907 kk=kk+1
2908 ksub=lisubm(kk)
2909 ENDDO
2910 jj=jj+1
2911
2912
2913 ENDIF
2914
2915 END DO
2916 END IF
2917
2918 IF (msegtyp(ce_loc(i)) < 0) THEN
2919 ie= - msegtyp(ce_loc(i))
2920 ELSE
2921 ie = ce_loc(i)
2922 ENDIF
2923 IF(ie > nrtm) ie=ie-nrtm
2924
2925 kk =addsubm(ie) ! address of m
2926 DO WHILE(kk<addsubm(ie+1))
2927! all sub interfaces of s
2928 ksub=lisubm(kk)
2929! KSUB Croissant
2930 itypsub = typsub(ksub)
2931 IF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
2932 impx=-fxt(i)*dt12
2933 impy=-fyt(i)*dt12
2934 impz=-fzt(i)*dt12
2935C main side :
2936 fsavsub1(4,ksub)=fsavsub1(4,ksub)+impx
2937 fsavsub1(5,ksub)=fsavsub1(5,ksub)+impy
2938 fsavsub1(6,ksub)=fsavsub1(6,ksub)+impz
2939C
2940 impx=fxi(i)*dt12
2941 impy=fyi(i)*dt12
2942 impz=fzi(i)*dt12
2943 fsavsub1(12,ksub)=fsavsub1(12,ksub)+abs(impx)
2944 fsavsub1(13,ksub)=fsavsub1(13,ksub)+abs(impy)
2945 fsavsub1(14,ksub)=fsavsub1(14,ksub)+abs(impz)
2946C
2947 IF(isensint(ksub+1)/=0) THEN
2948 fsavparit(ksub+1,4,i+nft) = -fxt(i)
2949 fsavparit(ksub+1,5,i+nft) = -fyt(i)
2950 fsavparit(ksub+1,6,i+nft) = -fzt(i)
2951 ENDIF
2952C
2953 fsavsub1(15,ksub)= fsavsub1(15,ksub)
2954 . +sqrt(impx*impx+impy*impy+impz*impz)
2955 fsavsub1(22,ksub)=fsavsub1(22,ksub)
2956 . +yp(i)*impz-zp(i)*impy
2957 fsavsub1(23,ksub)=fsavsub1(23,ksub)
2958 . +zp(i)*impx-xp(i)*impz
2959 fsavsub1(24,ksub)=fsavsub1(24,ksub)
2960 . +xp(i)*impy-yp(i)*impx
2961
2962 ENDIF
2963 kk=kk+1
2964 ENDDO
2965
2966 END DO
2967
2968
2969 IF(nspmd>1) THEN
2970C loop split because of a PGI bug
2971 DO i=1,jlt
2972 IF(pene(i) == zero)cycle
2973 nn = nsvg(i)
2974 IF(nn<0)THEN
2975 nn = -nn
2976 IF (msegtyp(ce_loc(i)) < 0) THEN
2977 ie= - msegtyp(ce_loc(i))
2978 ELSE
2979 ie = ce_loc(i)
2980 ENDIF
2981 jj =addsubsfi(nin)%P(nn)
2982 kk =addsubm(ie)
2983 DO WHILE(jj<addsubsfi(nin)%P(nn+1))
2984 jsub=lisubsfi(nin)%P(jj)
2985 itypsub = typsub(jsub)
2986
2987 IF(itypsub == 1 ) THEN
2988 iss1 = bitget(inflg_subsfi(nin)%P(jj),0)
2989 iss2 = bitget(inflg_subsfi(nin)%P(jj),1)
2990 igrn = bitget(inflg_subsfi(nin)%P(jj),2)
2991 ksub=lisubm(kk)
2992 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
2993 ims1 = bitget(inflg_subm(kk),0)
2994 ims2 = bitget(inflg_subm(kk),1)
2995 IF(ksub==jsub)THEN
2996 IF(.NOT.(((ims1 == 1 .OR. ims2 == 1) .AND. igrn == 1).OR.
2997 . (ims1 == 1 .AND. iss2 == 1).OR.
2998 . (ims2 == 1 .AND. iss1 == 1))) THEN
2999 kk=kk+1
3000 ksub=lisubm(kk)
3001 cycle
3002 END IF
3003 impx=fxt(i)*dt12
3004 impy=fyt(i)*dt12
3005 impz=fzt(i)*dt12
3006C main side :
3007 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
3008 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
3009 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
3010C
3011 impx=fxi(i)*dt12
3012 impy=fyi(i)*dt12
3013 impz=fzi(i)*dt12
3014 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
3015 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
3016 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
3017C
3018 IF(isensint(jsub+1)/=0) THEN
3019 fsavparit(jsub+1,4,i+nft) = fxt(i)
3020 fsavparit(jsub+1,5,i+nft) = fyt(i)
3021 fsavparit(jsub+1,6,i+nft) = fzt(i)
3022 ENDIF
3023C
3024 fsavsub1(15,jsub)= fsavsub1(15,jsub)
3025 . +sqrt(impx*impx+impy*impy+impz*impz)
3026 fsavsub1(22,jsub)=fsavsub1(22,jsub)
3027 . +yp(i)*impz-zp(i)*impy
3028 fsavsub1(23,jsub)=fsavsub1(23,jsub)
3029 . +zp(i)*impx-xp(i)*impz
3030 fsavsub1(24,jsub)=fsavsub1(24,jsub)
3031 . +xp(i)*impy-yp(i)*impx
3032
3033C
3034 END IF
3035
3036 kk=kk+1
3037 ksub=lisubm(kk)
3038 ENDDO
3039 jj=jj+1
3040
3041 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
3042
3043 impx=fxt(i)*dt12
3044 impy=fyt(i)*dt12
3045 impz=fzt(i)*dt12
3046C main side :
3047 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
3048 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
3049 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
3050C
3051 impx=fxi(i)*dt12
3052 impy=fyi(i)*dt12
3053 impz=fzi(i)*dt12
3054 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
3055
3056 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
3057 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
3058C
3059 IF(isensint(jsub+1)/=0) THEN
3060 fsavparit(jsub+1,4,i+nft) = fxt(i)
3061 fsavparit(jsub+1,5,i+nft) = fyt(i)
3062 fsavparit(jsub+1,6,i+nft) = fzt(i)
3063 ENDIF
3064C
3065 fsavsub1(15,jsub)= fsavsub1(15,jsub)
3066 . +sqrt(impx*impx+impy*impy+impz*impz)
3067 fsavsub1(22,jsub)=fsavsub1(22,jsub)
3068 . +yp(i)*impz-zp(i)*impy
3069 fsavsub1(23,jsub)=fsavsub1(23,jsub)
3070 . +zp(i)*impx-xp(i)*impz
3071 fsavsub1(24,jsub)=fsavsub1(24,jsub)
3072 . +xp(i)*impy-yp(i)*impx
3073
3074 jj=jj+1
3075
3076 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfaces
3077
3078 iss2 = bitget(inflg_subsfi(nin)%P(jj),0)
3079 iss1 = bitget(inflg_subsfi(nin)%P(jj),1)
3080 ksub=lisubm(kk)
3081 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
3082 ims2 = bitget(inflg_subm(kk),0)
3083 ims1 = bitget(inflg_subm(kk),1)
3084 IF(ksub==jsub)THEN
3085 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
3086 . (ims2 == 1 .AND. iss1 == 1))) THEN
3087 kk=kk+1
3088 ksub=lisubm(kk)
3089 cycle
3090 END IF
3091
3092 impx=fxi(i)*dt12
3093 impy=fyi(i)*dt12
3094 impz=fzi(i)*dt12
3095
3096 IF(ims2 > 0)THEN
3097 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
3098 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
3099 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
3100 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
3101 ELSE
3102 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
3103 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
3104 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
3105 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
3106 END IF
3107C
3108 IF(isensint(jsub+1)/=0) THEN
3109 IF(ims2 > 0 ) THEN
3110 fsavparit(jsub+1,4,i+nft) = -fxt(i)
3111 fsavparit(jsub+1,5,i+nft) = -fyt(i)
3112 fsavparit(jsub+1,6,i+nft) = -fzt(i)
3113 ELSE
3114 fsavparit(jsub+1,4,i+nft) = fxt(i)
3115 fsavparit(jsub+1,5,i+nft) = fyt(i)
3116 fsavparit(jsub+1,6,i+nft) = fzt(i)
3117 ENDIF
3118 ENDIF
3119
3120 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
3121 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
3122 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
3123
3124 ENDIF
3125 kk=kk+1
3126 ksub=lisubm(kk)
3127 ENDDO
3128 jj=jj+1
3129
3130 ENDIF
3131 END DO
3132 END IF
3133 END DO
3134 END IF
3135#include "lockon.inc"
3136 DO jsub=1,nisub
3137 nsub=lisub(jsub)
3138 DO j=1,15
3139 fsavsub(j,nsub)=fsavsub(j,nsub)+fsavsub1(j,jsub)
3140 END DO
3141 fsavsub(22,nsub)=fsavsub(22,nsub)+fsavsub1(22,jsub)
3142 fsavsub(23,nsub)=fsavsub(23,nsub)+fsavsub1(23,jsub)
3143 fsavsub(24,nsub)=fsavsub(24,nsub)+fsavsub1(24,jsub)
3144 IF(parameters%INTCAREA > 0) fsavsub(29,nsub)=fsavsub(29,nsub)+fsavsub1(25,jsub)
3145 END DO
3146#include "lockoff.inc"
3147 END IF
3148C---------------------------------
3149#include "lockon.inc"
3150 IF (inconv==1) THEN
3151 econtv = econtv + econvt ! frictional energy
3152 econt = econt + econtt ! Elastic Energy
3153 econtd = econtd + econtdt ! Damping Energy
3154 END IF !(INCONV==1) THEN
3155#include "lockoff.inc"
3156C---------------------------------
3157C---------------------------------
3158 IF(interefric >0)THEN
3159 IF (inconv==1) THEN
3160#include "lockon.inc"
3161 DO i=1,jlt
3162 IF(pene(i) == zero)cycle
3163 efric_lm = half*efric_l(i)
3164 output%DATA%EFRIC(interefric,ix1(i)) = output%DATA%EFRIC(interefric,ix1(i)) + h1(i)*efric_lm
3165 output%DATA%EFRIC(interefric,ix2(i)) = output%DATA%EFRIC(interefric,ix2(i)) + h2(i)*efric_lm
3166 output%DATA%EFRIC(interefric,ix3(i)) = output%DATA%EFRIC(interefric,ix3(i)) + h3(i)*efric_lm
3167 output%DATA%EFRIC(interefric,ix4(i)) = output%DATA%EFRIC(interefric,ix4(i)) + h4(i)*efric_lm
3168 jg = nsvg(i)
3169 efric_ls = half*efric_l(i)
3170 IF(jg>0) THEN
3171 output%DATA%EFRIC(interefric,jg) = output%DATA%EFRIC(interefric,jg) + efric_ls
3172 ELSE ! node remote
3173 jg = -jg
3174 efricfi(nin)%P(jg)=efricfi(nin)%P(jg)+ efric_ls
3175 ENDIF
3176 ENDDO
3177#include "lockoff.inc"
3178 END IF !(INCONV==1) THEN
3179 ENDIF
3180C---------------------------------
3181 IF(h3d_data%N_SCAL_CSE_FRIC >0)THEN
3182 IF (inconv==1) THEN
3183#include "lockon.inc"
3184 DO i=1,jlt
3185 IF(pene(i) == zero)cycle
3186 efric_lm = half*efric_l(i)
3187 output%DATA%EFRICG(ix1(i)) = output%DATA%EFRICG(ix1(i)) + h1(i)*efric_lm
3188 output%DATA%EFRICG(ix2(i)) = output%DATA%EFRICG(ix2(i)) + h2(i)*efric_lm
3189 output%DATA%EFRICG(ix3(i)) = output%DATA%EFRICG(ix3(i)) + h3(i)*efric_lm
3190 output%DATA%EFRICG(ix4(i)) = output%DATA%EFRICG(ix4(i)) + h4(i)*efric_lm
3191 jg = nsvg(i)
3192 efric_ls = half*efric_l(i)
3193 IF(jg>0) THEN
3194 output%DATA%EFRICG(jg) = output%DATA%EFRICG(jg) + efric_ls
3195 ELSE ! node remote
3196 jg = -jg
3197 efricgfi(nin)%P(jg)=efricgfi(nin)%P(jg)+ efric_ls
3198 ENDIF
3199 ENDDO
3200#include "lockoff.inc"
3201 END IF !(INCONV==1) THEN
3202 ENDIF
3203C---------------------------------
3204 IF(kdtint==1)THEN
3205 IF( (visc/=zero)
3206 . .AND.(ivis2==0.OR.ivis2==1))THEN
3207 DO i=1,jlt
3208 IF(pene(i) == zero)cycle
3209C C (i) = 2.*C (i)
3210C
3211 IF(msi(i)==zero)THEN
3212 ks(i) =zero
3213 cs(i) =zero
3214 stv(i)=zero
3215 ELSE
3216 cx = four*c(i)*c(i)
3217 cy = eight*msi(i)*kt(i)
3218 aux = sqrt(cx+cy)+two*c(i)
3219 stv(i)= kt(i)*aux*aux/max(cy,em30)
3220 aux = two*cf(i)*cf(i)/max(msi(i),em20)
3221 IF(aux>stv(i))THEN
3222 ks(i) =zero
3223 cs(i) =cf(i)
3224 stv(i)=aux
3225 ELSE
3226 ks(i)= kt(i)
3227 cs(i) =c(i)
3228 ENDIF
3229 ENDIF
3230C
3231 j1=ix1(i)
3232 IF(ms(j1)==zero)THEN
3233 k1(i) =zero
3234 c1(i) =zero
3235 st1(i)=zero
3236 ELSE
3237 k1(i)=kt(i)*abs(h1(i))
3238 c1(i)=c(i)*abs(h1(i))
3239 cx =four*c1(i)*c1(i)
3240 cy =eight*ms(j1)*k1(i)
3241 aux = sqrt(cx+cy)+two*c1(i)
3242 st1(i)= k1(i)*aux*aux/max(cy,em30)
3243 cfi = cf(i)*abs(h1(i))
3244 aux = two*cfi*cfi/max(ms(j1),em20)
3245 IF(aux>st1(i))THEN
3246 k1(i) =zero
3247 c1(i) =cfi
3248 st1(i)=aux
3249 ENDIF
3250 ENDIF
3251C
3252 j1=ix2(i)
3253 IF(ms(j1)==zero)THEN
3254 k2(i) =zero
3255 c2(i) =zero
3256 st2(i)=zero
3257 ELSE
3258 k2(i)=kt(i)*abs(h2(i))
3259 c2(i)=c(i)*abs(h2(i))
3260 cx =four*c2(i)*c2(i)
3261 cy =eight*ms(j1)*k2(i)
3262 aux = sqrt(cx+cy)+two*c2(i)
3263 st2(i)= k2(i)*aux*aux/max(cy,em30)
3264 cfi = cf(i)*abs(h2(i))
3265 aux = two*cfi*cfi/max(ms(j1),em20)
3266 IF(aux>st2(i))THEN
3267 k2(i) =zero
3268 c2(i) =cfi
3269 st2(i)=aux
3270 ENDIF
3271 ENDIF
3272C
3273 j1=ix3(i)
3274 IF(ms(j1)==zero)THEN
3275 k3(i) =zero
3276 c3(i) =zero
3277 st3(i)=zero
3278 ELSE
3279 k3(i)=kt(i)*abs(h3(i))
3280 c3(i)=c(i)*abs(h3(i))
3281 cx =four*c3(i)*c3(i)
3282 cy =eight*ms(j1)*k3(i)
3283 aux = sqrt(cx+cy)+two*c3(i)
3284 st3(i)= k3(i)*aux*aux/max(cy,em30)
3285 cfi = cf(i)*abs(h3(i))
3286 aux = two*cfi*cfi/max(ms(j1),em20)
3287 IF(aux>st3(i))THEN
3288 k3(i) =zero
3289 c3(i) =cfi
3290 st3(i)=aux
3291 ENDIF
3292 ENDIF
3293C
3294 j1=ix4(i)
3295 IF(ms(j1)==zero)THEN
3296 k4(i) =zero
3297 c4(i) =zero
3298 st4(i)=zero
3299 ELSE
3300 k4(i)=kt(i)*abs(h4(i))
3301 c4(i)=c(i)*abs(h4(i))
3302 cx =four*c4(i)*c4(i)
3303 cy =eight*ms(j1)*k4(i)
3304 aux = sqrt(cx+cy)+two*c4(i)
3305 st4(i)= k4(i)*aux*aux/max(cy,em30)
3306 cfi = cf(i)*abs(h4(i))
3307 aux = two*cfi*cfi/max(ms(j1),em20)
3308 IF(aux>st4(i))THEN
3309 k4(i) =zero
3310 c4(i) =cfi
3311 st4(i)=aux
3312 ENDIF
3313 ENDIF
3314 ENDDO
3315C
3316 ELSE
3317 DO i=1,jlt
3318 IF(viscffric(i)/=zero
3319 . .AND.(ivis2==0.OR.ivis2==1))THEN
3320 IF(pene(i) == zero)cycle
3321C C (i) = 2.*C (i)
3322C
3323 IF(msi(i)==zero)THEN
3324 ks(i) =zero
3325 cs(i) =zero
3326 stv(i)=zero
3327 ELSE
3328 cx = four*c(i)*c(i)
3329 cy = eight*msi(i)*kt(i)
3330 aux = sqrt(cx+cy)+two*c(i)
3331 stv(i)= kt(i)*aux*aux/max(cy,em30)
3332 aux = two*cf(i)*cf(i)/max(msi(i),em20)
3333 IF(aux>stv(i))THEN
3334 ks(i) =zero
3335 cs(i) =cf(i)
3336 stv(i)=aux
3337 ELSE
3338 ks(i)= kt(i)
3339 cs(i) =c(i)
3340 ENDIF
3341 ENDIF
3342C
3343 j1=ix1(i)
3344 IF(ms(j1)==zero)THEN
3345 k1(i) =zero
3346 c1(i) =zero
3347 st1(i)=zero
3348 ELSE
3349 k1(i)=kt(i)*abs(h1(i))
3350 c1(i)=c(i)*abs(h1(i))
3351 cx =four*c1(i)*c1(i)
3352 cy =eight*ms(j1)*k1(i)
3353 aux = sqrt(cx+cy)+two*c1(i)
3354 st1(i)= k1(i)*aux*aux/max(cy,em30)
3355 cfi = cf(i)*abs(h1(i))
3356 aux = two*cfi*cfi/max(ms(j1),em20)
3357 IF(aux>st1(i))THEN
3358 k1(i) =zero
3359 c1(i) =cfi
3360 st1(i)=aux
3361 ENDIF
3362 ENDIF
3363C
3364 j1=ix2(i)
3365 IF(ms(j1)==zero)THEN
3366 k2(i) =zero
3367 c2(i) =zero
3368 st2(i)=zero
3369 ELSE
3370 k2(i)=kt(i)*abs(h2(i))
3371 c2(i)=c(i)*abs(h2(i))
3372 cx =four*c2(i)*c2(i)
3373 cy =eight*ms(j1)*k2(i)
3374 aux = sqrt(cx+cy)+two*c2(i)
3375 st2(i)= k2(i)*aux*aux/max(cy,em30)
3376 cfi = cf(i)*abs(h2(i))
3377 aux = two*cfi*cfi/max(ms(j1),em20)
3378 IF(aux>st2(i))THEN
3379 k2(i) =zero
3380 c2(i) =cfi
3381 st2(i)=aux
3382 ENDIF
3383 ENDIF
3384C
3385 j1=ix3(i)
3386 IF(ms(j1)==zero)THEN
3387 k3(i) =zero
3388 c3(i) =zero
3389 st3(i)=zero
3390 ELSE
3391 k3(i)=kt(i)*abs(h3(i))
3392 c3(i)=c(i)*abs(h3(i))
3393 cx =four*c3(i)*c3(i)
3394 cy =eight*ms(j1)*k3(i)
3395 aux = sqrt(cx+cy)+two*c3(i)
3396 st3(i)= k3(i)*aux*aux/max(cy,em30)
3397 cfi = cf(i)*abs(h3(i))
3398 aux = two*cfi*cfi/max(ms(j1),em20)
3399 IF(aux>st3(i))THEN
3400 k3(i) =zero
3401 c3(i) =cfi
3402 st3(i)=aux
3403 ENDIF
3404 ENDIF
3405C
3406 j1=ix4(i)
3407 IF(ms(j1)==zero)THEN
3408 k4(i) =zero
3409 c4(i) =zero
3410 st4(i)=zero
3411 ELSE
3412 k4(i)=kt(i)*abs(h4(i))
3413 c4(i)=c(i)*abs(h4(i))
3414 cx =four*c4(i)*c4(i)
3415 cy =eight*ms(j1)*k4(i)
3416 aux = sqrt(cx+cy)+two*c4(i)
3417 st4(i)= k4(i)*aux*aux/max(cy,em30)
3418 cfi = cf(i)*abs(h4(i))
3419 aux = two*cfi*cfi/max(ms(j1),em20)
3420 IF(aux>st4(i))THEN
3421 k4(i) =zero
3422 c4(i) =cfi
3423 st4(i)=aux
3424 ENDIF
3425 ENDIF
3426 ELSE
3427 IF(pene(i) == zero)cycle
3428 ks(i) =stif(i)
3429 cs(i) =zero
3430 stv(i)=ks(i)
3431 k1(i) =stif(i)*abs(h1(i))
3432 c1(i) =zero
3433 st1(i)=k1(i)
3434 k2(i) =stif(i)*abs(h2(i))
3435 c2(i) =zero
3436 st2(i)=k2(i)
3437 k3(i) =stif(i)*abs(h3(i))
3438 c3(i) =zero
3439 st3(i)=k3(i)
3440 k4(i) =stif(i)*abs(h4(i))
3441 c4(i) =zero
3442 st4(i)=k4(i)
3443 ENDIF
3444 ENDDO
3445 ENDIF
3446 ENDIF
3447C
3448C=======================================================================
3449C---------------------------------
3450 itag = 0
3451 DO i=1,jlt
3452 IF(pene(i) == zero)cycle
3453!!
3454 fx1(i)=fxi(i)*h1(i)
3455 fy1(i)=fyi(i)*h1(i)
3456 fz1(i)=fzi(i)*h1(i)
3457C
3458 fx2(i)=fxi(i)*h2(i)
3459 fy2(i)=fyi(i)*h2(i)
3460 fz2(i)=fzi(i)*h2(i)
3461C
3462 fx3(i)=fxi(i)*h3(i)
3463 fy3(i)=fyi(i)*h3(i)
3464 fz3(i)=fzi(i)*h3(i)
3465C
3466 fx4(i)=fxi(i)*h4(i)
3467 fy4(i)=fyi(i)*h4(i)
3468 fz4(i)=fzi(i)*h4(i)
3469C
3470 phi1(i) = zero
3471 phi2(i) = zero
3472 phi3(i) = zero
3473 phi4(i) = zero
3474 ENDDO
3475
3476
3477 IF(idtmins==2.OR.idtmins_int/=0)THEN
3478 dti=dt2t
3479 CALL i24sms2(jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
3480 2 nsvg ,h1 ,h2 ,h3 ,h4 ,stif ,
3481 3 nin ,noint ,mskyi_sms, iskyi_sms,nsms ,
3482 4 kt ,c ,cf ,dtmini,dti ,
3483 5 irtse ,nsne ,is2se,is2pt ,t2main_sms ,
3484 6 t2fac_sms)
3485 IF(dti<dt2t)THEN
3486 dt2t = dti
3487 neltst = noint
3488 ityptst = 10
3489 ENDIF
3490 ENDIF
3491C
3492 IF(idtmins_int/=0)THEN
3493 stif(1:jlt)=zero
3494 END IF
3495C------------For /LOAD/PRESSURE tag nodes in contact-------------
3496 tagip(1:jlt) = 0
3497 IF(ninloadp > 0) THEN
3498 DO k = kloadpinter(nin)+1, kloadpinter(nin+1)
3499 pp = loadpinter(k)
3500 ppl = loadp_hyd_inter(pp)
3501 dgapload = dgaploadint(k)
3502 DO i=1,jlt
3503 gapp= gapv(i) + dgapload
3504 IF((pene(i)/=zero).OR.(dist(i) <= gapp.AND.dist(i) >= zero)) THEN
3505 tagip(i) = 1
3506 IF(pene(i)==zero) THEN
3507 ix1(i) = ixx(i,1)
3508 ix2(i) = ixx(i,2)
3509 ix3(i) = ixx(i,3)
3510 ix4(i) = ixx(i,4)
3511 ENDIF
3512 tagncont(ppl,ixx(i,1)) = 1
3513 tagncont(ppl,ixx(i,2)) = 1
3514 tagncont(ppl,ixx(i,3)) = 1
3515 tagncont(ppl,ixx(i,4)) = 1
3516 jg = nsvg(i)
3517 IF(jg>0) THEN
3518C SPMD : do same after reception of forces for remote nodes
3519 IF (jg <= numnod) THEN
3520 tagncont(ppl,jg) = 1
3521 ELSE
3522 ig = jg - numnod
3523 IF (ig > 0) THEN
3524 IF (is2se(1,ig) > 0) THEN
3525 ie = is2se(1,ig)
3526 ELSEIF(is2se(2,ig) > 0) THEN
3527 ie = is2se(2,ig)
3528 END IF
3529 DO j =1,4
3530 tagncont(ppl,irtse(j,ie)) = 1
3531 END DO
3532 ENDIF
3533 ENDIF
3534 ENDIF
3535 ENDIF
3536 ENDDO
3537 ENDDO
3538
3539 ENDIF
3540C
3541C
3542C SPMD: Identification of interf nodes.useful to send
3543 IF (nspmd>1) THEN
3544ctmp+1 mic only
3545#include "mic_lockon.inc"
3546 DO i = 1,jlt
3547 nn = nsvg(i)
3548 IF(nn<0)THEN
3549C temporary tag of nsvfi a -
3550 hh = h1(i)+h2(i)+h3(i)+h4(i)
3551 IF(hh /= zero.OR.tagip(i)==1)THEN
3552 IF(isedge_fi(nin)%P(-nn)==0)THEN
3553 nsvfi(nin)%P(-nn) = -abs(nsvfi(nin)%P(-nn))
3554 ELSE
3555 j1=irtse_fi(nin)%P(1,-nn)
3556 nsvfi(nin)%P(j1) = -abs(nsvfi(nin)%P(j1))
3557 j1=irtse_fi(nin)%P(2,-nn)
3558 nsvfi(nin)%P(j1) = -abs(nsvfi(nin)%P(j1))
3559 j1=irtse_fi(nin)%P(3,-nn)
3560 nsvfi(nin)%P(j1) = -abs(nsvfi(nin)%P(j1))
3561 j1=irtse_fi(nin)%P(4,-nn)
3562 nsvfi(nin)%P(j1) = -abs(nsvfi(nin)%P(j1))
3563 ENDIF
3564 ENDIF
3565 ENDIF
3566 ENDDO
3567ctmp+1 mic only
3568#include "mic_lockoff.inc"
3569 ENDIF
3570
3571C
3572 IF (impl_s ==0.AND.(inacti==-1.OR.igsti==-1)) stif(1:jlt)=max(stif(1:jlt),stifkt(1:jlt))
3573 IF(iparit==3)THEN
3574 stop 789
3575 ELSEIF(iparit==0)THEN
3576 CALL i24ass0(jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
3577 2 nsvg ,h1 ,h2 ,h3 ,h4 ,stif ,
3578 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,fz2 ,
3579 4 fx3 ,fy3 ,fz3 ,fx4 ,fy4 ,fz4 ,
3580 5 fxi ,fyi ,fzi ,a ,stifn ,nin ,
3581 6 intth ,phi ,fthe ,phi1 , phi2 ,phi3 ,
3582 7 phi4 ,intply ,iply ,inod_pxfem ,
3583 8 irtse ,nsne ,is2se ,is2pt,jtask )
3584
3585 ELSE
3586 CALL i24ass2(jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
3587 2 nsvg ,h1 ,h2 ,h3 ,h4 ,stif ,
3588 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,fz2 ,
3589 4 fx3 ,fy3 ,fz3 ,fx4 ,fy4 ,fz4 ,
3590 5 fxi ,fyi ,fzi ,fskyi,isky ,niskyfi,
3591 6 nin ,noint ,intth,phi ,ftheskyi ,phi1,
3592 7 phi2 ,phi3 , phi4 ,itab ,intply,iply ,
3593 8 inod_pxfem,irtse,nsne,is2se,is2pt,tagip )
3594 ENDIF
3595C
3596 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT>0)THEN
3597 IF (inconv==1) THEN
3598#include "lockon.inc"
3599 DO i=1,jlt
3600 IF(pene(i) == zero)cycle
3601 fcont(1,ix1(i)) =fcont(1,ix1(i)) + fx1(i)
3602 fcont(2,ix1(i)) =fcont(2,ix1(i)) + fy1(i)
3603 fcont(3,ix1(i)) =fcont(3,ix1(i)) + fz1(i)
3604 fcont(1,ix2(i)) =fcont(1,ix2(i)) + fx2(i)
3605 fcont(2,ix2(i)) =fcont(2,ix2(i)) + fy2(i)
3606 fcont(3,ix2(i)) =fcont(3,ix2(i)) + fz2(i)
3607 fcont(1,ix3(i)) =fcont(1,ix3(i)) + fx3(i)
3608 fcont(2,ix3(i)) =fcont(2,ix3(i)) + fy3(i)
3609 fcont(3,ix3(i)) =fcont(3,ix3(i)) + fz3(i)
3610 fcont(1,ix4(i)) =fcont(1,ix4(i)) + fx4(i)
3611 fcont(2,ix4(i)) =fcont(2,ix4(i)) + fy4(i)
3612 fcont(3,ix4(i)) =fcont(3,ix4(i)) + fz4(i)
3613 jg = nsvg(i)
3614 IF(jg>0) THEN
3615C In SPMD: Treatment to be redone after reception node Remote if JG <0
3616 IF (jg >numnod) THEN
3617 ig = jg - numnod
3618 CALL i24for1_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
3619 + ig ,fxi(i) ,fyi(i) ,fzi(i) ,fcont ,
3620 + -1 )
3621 ELSE
3622 fcont(1,jg)=fcont(1,jg)- fxi(i)
3623 fcont(2,jg)=fcont(2,jg)- fyi(i)
3624 fcont(3,jg)=fcont(3,jg)- fzi(i)
3625 ENDIF
3626
3627 ELSE
3628! On QA test ../miniqa/INTERF/INT_24/implicit_Int24_Ipen0/data/CONVI7TU_0001.rad 2 x 2 JG = -1
3629C IF(ISEDGE_FI(NIN)%P(JG)==1)THEN
3630C ELSE
3631C ENDIF
3632
3633 ENDIF
3634 ENDDO
3635#include "lockoff.inc"
3636 END IF !(INCONV==1) THEN
3637 ENDIF
3638C-----------------------------------------------------
3639 IF(isecin>0.AND.inconv==1)THEN
3640 k0=nstrf(25)
3641 IF(nstrf(1)+nstrf(2)/=0)THEN
3642 DO i=1,nsect
3643 nbinter=nstrf(k0+14)
3644 k1s=k0+30
3645 DO j=1,nbinter
3646 IF(nstrf(k1s)==noint)THEN
3647 IF(isecut/=0)THEN
3648#include "lockon.inc"
3649 DO k=1,jlt
3650 IF(pene(k) == zero)cycle
3651C pay attention to the signs for the accumulation of forces
3652C To make it conform with CFORC3
3653 IF(secfcum(4,ix1(k),i)==1.)THEN
3654 secfcum(1,ix1(k),i)=secfcum(1,ix1(k),i)-fx1(k)
3655 secfcum(2,ix1(k),i)=secfcum(2,ix1(k),i)-fy1(k)
3656 secfcum(3,ix1(k),i)=secfcum(3,ix1(k),i)-fz1(k)
3657 ENDIF
3658 IF(secfcum(4,ix2(k),i)==1.)THEN
3659 secfcum(1,ix2(k),i)=secfcum(1,ix2(k),i)-fx2(k)
3660 secfcum(2,ix2(k),i)=secfcum(2,ix2(k),i)-fy2(k)
3661 secfcum(3,ix2(k),i)=secfcum(3,ix2(k),i)-fz2(k)
3662 ENDIF
3663 IF(secfcum(4,ix3(k),i)==1.)THEN
3664 secfcum(1,ix3(k),i)=secfcum(1,ix3(k),i)-fx3(k)
3665 secfcum(2,ix3(k),i)=secfcum(2,ix3(k),i)-fy3(k)
3666 secfcum(3,ix3(k),i)=secfcum(3,ix3(k),i)-fz3(k)
3667 ENDIF
3668 IF(secfcum(4,ix4(k),i)==1.)THEN
3669 secfcum(1,ix4(k),i)=secfcum(1,ix4(k),i)-fx4(k)
3670 secfcum(2,ix4(k),i)=secfcum(2,ix4(k),i)-fy4(k)
3671 secfcum(3,ix4(k),i)=secfcum(3,ix4(k),i)-fz4(k)
3672 ENDIF
3673 jg = nsvg(k)
3674 IF(jg>0) THEN
3675C In SPMD: Treatment to be redone after reception node Remote if JG <0
3676 IF (jg >numnod) THEN
3677 ig = jg - numnod
3678 IF(secfcum(4,jg,i)==1.)THEN
3679 CALL i24for1_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
3680 + ig ,fxi(i) ,fyi(i) ,fzi(i) ,secfcum(1,1,i),
3681 + 1 )
3682 ENDIF
3683 ELSE
3684 IF(secfcum(4,jg,i)==1.)THEN
3685 secfcum(1,jg,i)=secfcum(1,jg,i)+fxi(k)
3686 secfcum(2,jg,i)=secfcum(2,jg,i)+fyi(k)
3687 secfcum(3,jg,i)=secfcum(3,jg,i)+fzi(k)
3688 ENDIF
3689 ENDIF
3690 ELSE
3691C Node Remote
3692 jg=-jg
3693 IF(isedge_fi(nin)%P(jg)==1)THEN
3694 IF(secfcum(4,jg,i)==1.)THEN
3695 CALL i24for1_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
3696 + jg ,fxi(i) ,fyi(i) ,fzi(i) ,secfcum(1,1,i),
3697 + 1 )
3698 ENDIF
3699 ELSE
3700 IF(secfcum(4,jg,i)==1.)THEN
3701 secfcum(1,jg,i)=secfcum(1,jg,i)+fxi(k)
3702 secfcum(2,jg,i)=secfcum(2,jg,i)+fyi(k)
3703 secfcum(3,jg,i)=secfcum(3,jg,i)+fzi(k)
3704 ENDIF
3705 ENDIF
3706 ENDIF
3707 ENDDO
3708#include "lockoff.inc"
3709 ENDIF
3710C +fsav(section)
3711 ENDIF
3712 k1s=k1s+1
3713 ENDDO
3714 k0=nstrf(k0+24)
3715 ENDDO
3716 ENDIF
3717 ENDIF
3718C-----------------------------------------------------
3719 IF(ibag/=0.OR.iadm/=0)THEN
3720 DO i=1,jlt
3721 IF(pene(i) == zero)cycle
3722C test modified for consistency with spmd communication (spmd_i7tools)
3723c IF(FXI(I)/=ZERO.OR.FYI(I)/=ZERO.OR.FZI(I)/=ZERO)THEN
3724 jg = nsvg(i)
3725 IF(jg>0) THEN
3726C In SPMD: Treatment to be redone after reception node Remote if JG <0
3727 icontact(jg)=1
3728 ENDIF
3729 icontact(ix1(i))=1
3730 icontact(ix2(i))=1
3731 icontact(ix3(i))=1
3732 icontact(ix4(i))=1
3733 ENDDO
3734 ENDIF
3735C
3736 IF(iadm/=0)THEN
3737 END IF
3738 IF(iadm>=2)THEN
3739 END IF
3740C
3741 IF(ibcc==0) RETURN
3742C
3743 DO i=1,jlt
3744 IF(pene(i) == zero)cycle
3745 ibcm = ibcc / 8
3746 ibcs = ibcc - 8 * ibcm
3747 IF(ibcs>0) THEN
3748 ig=nsvg(i)
3749C---------obsolate option, no need to update
3750 IF(ig>0.AND.ig<=numnod) THEN
3751C In SPMD: Treatment to be redone after reception node Remote if JG <0
3752 CALL ibcoff(ibcs,icodt(ig))
3753 ENDIF
3754 ENDIF
3755 IF(ibcm>0) THEN
3756 ig=ix1(i)
3757 CALL ibcoff(ibcm,icodt(ig))
3758 ig=ix2(i)
3759 CALL ibcoff(ibcm,icodt(ig))
3760 ig=ix3(i)
3761 CALL ibcoff(ibcm,icodt(ig))
3762 ig=ix4(i)
3763 CALL ibcoff(ibcm,icodt(ig))
3764 ENDIF
3765 ENDDO
3766C
3767 RETURN
3768 END
3769!||====================================================================
3770!|| i24ass2 ../engine/source/interfaces/int24/i24for3.F
3771!||--- called by ------------------------------------------------------
3772!|| i24for3 ../engine/source/interfaces/int24/i24for3.f
3773!||--- calls -----------------------------------------------------
3774!|| ancmsg ../engine/source/output/message/message.F
3775!|| arret ../engine/source/system/arret.F
3776!|| i24forc_fic ../engine/source/interfaces/int24/i24for3e.F
3777!||--- uses -----------------------------------------------------
3778!|| debug_mod ../engine/share/modules/debug_mod.F
3779!|| message_mod ../engine/share/message_module/message_mod.F
3780!|| plyxfem_mod ../engine/share/modules/plyxfem_mod.F
3781!|| tri7box ../engine/share/modules/tri7box.F
3782!||====================================================================
3783 SUBROUTINE i24ass2(JLT ,IX1 ,IX2 ,IX3 ,IX4 ,
3784 2 NSVG ,H1 ,H2 ,H3 ,H4 ,STIF ,
3785 3 FX1 ,FY1 ,FZ1 ,FX2 ,FY2 ,FZ2 ,
3786 4 FX3 ,FY3 ,FZ3 ,FX4 ,FY4 ,FZ4 ,
3787 5 FXI ,FYI ,FZI ,FSKYI,ISKY ,NISKYFI,
3788 6 NIN ,NOINT ,INTTH,PHI ,FTHESKYI,PHI1,
3789 7 PHI2 ,PHI3 , PHI4 ,ITAB,INTPLY ,IPLY ,
3790 8 INOD ,IRTSE ,NSNE ,IS2SE ,IS2PT,TAGIP)
3791C-----------------------------------------------
3792C M o d u l e s
3793C-----------------------------------------------
3794 USE tri7box
3795 USE message_mod
3796 USE plyxfem_mod
3797 USE debug_mod
3798C-----------------------------------------------
3799C I m p l i c i t T y p e s
3800C-----------------------------------------------
3801#include "implicit_f.inc"
3802#include "comlock.inc"
3803C-----------------------------------------------
3804C G l o b a l P a r a m e t e r s
3805C-----------------------------------------------
3806#include "mvsiz_p.inc"
3807C-----------------------------------------------
3808C C o m m o n B l o c k s
3809C-----------------------------------------------
3810#include "parit_c.inc"
3811#include "com04_c.inc"
3812C-----------------------------------------------
3813C D u m m y A r g u m e n t s
3814C-----------------------------------------------
3815 INTEGER JLT,NISKYFI,NIN,NOINT,INTTH,
3816 . ISKY(*),ITAB(*),
3817 . IX1(MVSIZ),IX2(MVSIZ),IX3(MVSIZ),IX4(MVSIZ),NSVG(MVSIZ),
3818 . INTPLY,IPLY(4,*),INOD(*),
3819 . IRTSE(5,*),NSNE ,IS2SE(2,*),IS2PT(*)
3820 INTEGER , INTENT(INOUT) :: TAGIP(MVSIZ)
3821 my_real
3822 . H1(MVSIZ),H2(MVSIZ),H3(MVSIZ),H4(MVSIZ),STIF(MVSIZ),
3823 . FX1(MVSIZ),FY1(MVSIZ),FZ1(MVSIZ),
3824 . FX2(MVSIZ),FY2(MVSIZ),FZ2(MVSIZ),
3825 . FX3(MVSIZ),FY3(MVSIZ),FZ3(MVSIZ),
3826 . fx4(mvsiz),fy4(mvsiz),fz4(mvsiz),
3827 . fxi(mvsiz),fyi(mvsiz),fzi(mvsiz),
3828 . fskyi(lskyi,nfskyi),ftheskyi(lskyi),phi(mvsiz),
3829 . phi1(*),phi2(*) ,phi3(*) ,phi4(*)
3830C-----------------------------------------------
3831C L o c a l V a r i a b l e s
3832C-----------------------------------------------
3833 my_real
3834 . hh(mvsiz),fici(5),fics(4,4),ficsth(4,5)
3835 INTEGER I, J1, IG, NISKYL1, NISKYL, IGP, IGM, NISKYFIL, J
3836 INTEGER JG, IXSS(4), NFIC
3837
3838C
3839 niskyl1 = 0
3840 tagip(1:jlt) = 0
3841 DO i = 1, jlt
3842 hh(i)=h1(i)+h2(i)+h3(i)+h4(i)
3843 IF(hh(i)==zero.AND.tagip(i)==0)cycle
3844 IF (h1(i)/=zero.OR.tagip(i)==1) niskyl1 = niskyl1 + 1
3845 IF (h2(i)/=zero.OR.tagip(i)==1) niskyl1 = niskyl1 + 1
3846 IF (h3(i)/=zero.OR.tagip(i)==1) niskyl1 = niskyl1 + 1
3847 IF (h4(i)/=zero.OR.tagip(i)==1) niskyl1 = niskyl1 + 1
3848 ENDDO
3849C
3850C Precalcul impact local / remote
3851C
3852 igp = 0
3853 igm = 0
3854 DO i=1,jlt
3855 IF(hh(i)==zero.AND.tagip(i)==0)cycle
3856 ig =nsvg(i)
3857 IF(ig>0) THEN
3858 igp = igp+1
3859 IF (ig > numnod) igp = igp+3
3860 ELSE
3861 igm = igm+1
3862 IF(isedge_fi(nin)%P(-ig)==1)igm = igm+3
3863 ENDIF
3864 ENDDO
3865C
3866#include "lockon.inc"
3867 niskyl = nisky
3868 nisky = nisky + niskyl1 + igp
3869 niskyfil = niskyfi
3870 niskyfi = niskyfi + igm
3871#include "lockoff.inc"
3872C
3873 IF (niskyl+niskyl1+igp > lskyi) THEN
3874 CALL ancmsg(msgid=26,anmode=aninfo)
3875 CALL arret(2)
3876 ENDIF
3877 IF (niskyfil+igm > nlskyfi(nin)) THEN
3878 CALL ancmsg(msgid=26,anmode=aninfo)
3879 CALL arret(2)
3880 ENDIF
3881C
3882 IF(intth == 0 ) THEN
3883 IF(intply == 0 )THEN
3884 DO i=1,jlt
3885 IF (h1(i)/=zero.OR.tagip(i)==1) THEN
3886 niskyl = niskyl + 1
3887 fskyi(niskyl,1)=fx1(i)
3888 fskyi(niskyl,2)=fy1(i)
3889 fskyi(niskyl,3)=fz1(i)
3890 fskyi(niskyl,4)=stif(i)*abs(h1(i))
3891 isky(niskyl) = ix1(i)
3892 ENDIF
3893 ENDDO
3894 DO i=1,jlt
3895 IF (h2(i)/=zero.OR.tagip(i)==1) THEN
3896 niskyl = niskyl + 1
3897 fskyi(niskyl,1)=fx2(i)
3898 fskyi(niskyl,2)=fy2(i)
3899 fskyi(niskyl,3)=fz2(i)
3900 fskyi(niskyl,4)=stif(i)*abs(h2(i))
3901 isky(niskyl) = ix2(i)
3902 ENDIF
3903 ENDDO
3904 DO i=1,jlt
3905 IF (h3(i)/=zero.OR.tagip(i)==1) THEN
3906 niskyl = niskyl + 1
3907 fskyi(niskyl,1)=fx3(i)
3908 fskyi(niskyl,2)=fy3(i)
3909 fskyi(niskyl,3)=fz3(i)
3910 fskyi(niskyl,4)=stif(i)*abs(h3(i))
3911 isky(niskyl) = ix3(i)
3912 ENDIF
3913 ENDDO
3914 DO i=1,jlt
3915 IF (h4(i)/=zero.OR.tagip(i)==1) THEN
3916 niskyl = niskyl + 1
3917 fskyi(niskyl,1)=fx4(i)
3918 fskyi(niskyl,2)=fy4(i)
3919 fskyi(niskyl,3)=fz4(i)
3920 fskyi(niskyl,4)=stif(i)*abs(h4(i))
3921 isky(niskyl) = ix4(i)
3922 ENDIF
3923 ENDDO
3924C
3925 DO i=1,jlt
3926 IF(hh(i)==zero.AND.tagip(i)==0)cycle
3927 ig =nsvg(i)
3928 IF(ig>0) THEN
3929 IF (ig >numnod) THEN
3930 jg = ig - numnod
3931 nfic=1
3932 fici(nfic)=fxi(i)
3933 nfic = nfic + 1
3934 fici(nfic)=fyi(i)
3935 nfic = nfic + 1
3936 fici(nfic)=fzi(i)
3937 nfic = nfic + 1
3938 fici(nfic)=stif(i)
3939 CALL i24forc_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
3940 + jg ,nfic ,fici ,fics ,ixss )
3941 DO j = 1,4
3942 j1 = ixss(j)
3943 niskyl = niskyl + 1
3944 fskyi(niskyl,1)=-fics(j,1)
3945 fskyi(niskyl,2)=-fics(j,2)
3946 fskyi(niskyl,3)=-fics(j,3)
3947 fskyi(niskyl,4)= fics(j,4)
3948 isky(niskyl) = j1
3949 END DO
3950 ELSE
3951 niskyl = niskyl + 1
3952 fskyi(niskyl,1)=-fxi(i)
3953 fskyi(niskyl,2)=-fyi(i)
3954 fskyi(niskyl,3)=-fzi(i)
3955 fskyi(niskyl,4)= stif(i)
3956 isky(niskyl) = ig
3957 END IF !(IG >NUMNOD) THEN
3958 ELSE
3959 ig = -ig
3960 IF(isedge_fi(nin)%P(ig)==1)THEN
3961 nfic=1
3962 fici(nfic)=fxi(i)
3963 nfic = nfic + 1
3964 fici(nfic)=fyi(i)
3965 nfic = nfic + 1
3966 fici(nfic)=fzi(i)
3967 nfic = nfic + 1
3968 fici(nfic)=stif(i)
3969 CALL i24forc_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) , is2pt_fi(nin)%P(1) ,
3970 + ig ,nfic ,fici ,fics ,ixss )
3971
3972 DO j = 1,4
3973 j1 = ixss(j)
3974 niskyfil = niskyfil + 1
3975 fskyfi(nin)%P(1,niskyfil)=-fics(j,1)
3976 fskyfi(nin)%P(2,niskyfil)=-fics(j,2)
3977 fskyfi(nin)%P(3,niskyfil)=-fics(j,3)
3978 fskyfi(nin)%P(4,niskyfil)= fics(j,4)
3979 iskyfi(nin)%P(niskyfil) = j1
3980 END DO
3981
3982 ELSE
3983 niskyfil = niskyfil + 1
3984 fskyfi(nin)%P(1,niskyfil)=-fxi(i)
3985 fskyfi(nin)%P(2,niskyfil)=-fyi(i)
3986 fskyfi(nin)%P(3,niskyfil)=-fzi(i)
3987 fskyfi(nin)%P(4,niskyfil)= stif(i)
3988 iskyfi(nin)%P(niskyfil) = ig
3989 ENDIF
3990 ENDIF
3991 ENDDO
3992 ELSE
3993C with plyxfem formulation s
3994 DO i=1,jlt
3995 IF (h1(i)/=zero.OR.tagip(i)==1) THEN
3996 niskyl = niskyl + 1
3997 fskyi(niskyl,1)=fx1(i)
3998 fskyi(niskyl,2)=fy1(i)
3999 fskyi(niskyl,3)=fz1(i)
4000 fskyi(niskyl,4)=stif(i)*abs(h1(i))
4001 isky(niskyl) = ix1(i)
4002C
4003 plyskyi%FSKYI(niskyl,1)=fx1(i)
4004 plyskyi%FSKYI(niskyl,2)=fy1(i)
4005 plyskyi%FSKYI(niskyl,3)=fz1(i)
4006 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h1(i))
4007 plyskyi%FSKYI(niskyl,5)=iply(1,i)
4008 ENDIF
4009 ENDDO
4010 DO i=1,jlt
4011 IF (h2(i)/=zero.OR.tagip(i)==1) THEN
4012 niskyl = niskyl + 1
4013 fskyi(niskyl,1)=fx2(i)
4014 fskyi(niskyl,2)=fy2(i)
4015 fskyi(niskyl,3)=fz2(i)
4016 fskyi(niskyl,4)=stif(i)*abs(h2(i))
4017 isky(niskyl) = ix2(i)
4018C
4019 plyskyi%FSKYI(niskyl,1)=fx2(i)
4020 plyskyi%FSKYI(niskyl,2)=fy2(i)
4021 plyskyi%FSKYI(niskyl,3)=fz2(i)
4022 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h2(i))
4023 plyskyi%FSKYI(niskyl,5)=iply(2,i)
4024 ENDIF
4025 ENDDO
4026 DO i=1,jlt
4027 IF (h3(i)/=zero.OR.tagip(i)==1) THEN
4028 niskyl = niskyl + 1
4029 fskyi(niskyl,1)=fx3(i)
4030 fskyi(niskyl,2)=fy3(i)
4031 fskyi(niskyl,3)=fz3(i)
4032 fskyi(niskyl,4)=stif(i)*abs(h3(i))
4033 isky(niskyl) = ix3(i)
4034C
4035 plyskyi%FSKYI(niskyl,1)=fx3(i)
4036 plyskyi%FSKYI(niskyl,2)=fy3(i)
4037 plyskyi%FSKYI(niskyl,3)=fz3(i)
4038 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h3(i))
4039 plyskyi%FSKYI(niskyl,5)=iply(3,i)
4040 ENDIF
4041 ENDDO
4042 DO i=1,jlt
4043 IF (h4(i)/=zero.OR.tagip(i)==1) THEN
4044 niskyl = niskyl + 1
4045 fskyi(niskyl,1)=fx4(i)
4046 fskyi(niskyl,2)=fy4(i)
4047 fskyi(niskyl,3)=fz4(i)
4048 fskyi(niskyl,4)=stif(i)*abs(h4(i))
4049 isky(niskyl) = ix4(i)
4050C
4051 plyskyi%FSKYI(niskyl,1)=fx4(i)
4052 plyskyi%FSKYI(niskyl,2)=fy4(i)
4053 plyskyi%FSKYI(niskyl,3)=fz4(i)
4054 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h4(i))
4055 plyskyi%FSKYI(niskyl,5)=iply(4,i)
4056 ENDIF
4057 ENDDO
4058C
4059 DO i=1,jlt
4060 IF(hh(i)==zero.AND.tagip(i)==0)cycle
4061 ig =nsvg(i)
4062 IF(ig>0) THEN
4063 IF (ig >numnod) THEN
4064 jg = ig - numnod
4065 nfic=1
4066 fici(nfic)=fxi(i)
4067 nfic = nfic + 1
4068 fici(nfic)=fyi(i)
4069 nfic = nfic + 1
4070 fici(nfic)=fzi(i)
4071 nfic = nfic + 1
4072 fici(nfic)=stif(i)
4073 CALL i24forc_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
4074 + jg ,nfic ,fici ,fics ,ixss )
4075 DO j = 1,4
4076 j1 = ixss(j)
4077 niskyl = niskyl + 1
4078 fskyi(niskyl,1)=-fics(j,1)
4079 fskyi(niskyl,2)=-fics(j,2)
4080 fskyi(niskyl,3)=-fics(j,3)
4081 fskyi(niskyl,4)= fics(j,4)
4082 isky(niskyl) = j1
4083 END DO
4084 ELSE
4085 niskyl = niskyl + 1
4086 fskyi(niskyl,1)=-fxi(i)
4087 fskyi(niskyl,2)=-fyi(i)
4088 fskyi(niskyl,3)=-fzi(i)
4089 fskyi(niskyl,4)= stif(i)
4090 isky(niskyl) = ig
4091 END IF !(IG >NUMNOD) THEN
4092C
4093 ELSE
4094 ig = -ig
4095 IF(isedge_fi(nin)%P(ig)==1)THEN
4096 jg = ig
4097 nfic=1
4098 fici(nfic)=fxi(i)
4099 nfic = nfic + 1
4100 fici(nfic)=fyi(i)
4101 nfic = nfic + 1
4102 fici(nfic)=fzi(i)
4103 nfic = nfic + 1
4104 fici(nfic)=stif(i)
4105 CALL i24forc_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1),
4106 + jg ,nfic ,fici ,fics ,ixss )
4107 DO j = 1,4
4108 j1 = ixss(j)
4109 niskyfil = niskyfil + 1
4110 fskyfi(nin)%P(1,niskyfil)=-fics(j,1)
4111 fskyfi(nin)%P(2,niskyfil)=-fics(j,2)
4112 fskyfi(nin)%P(3,niskyfil)=-fics(j,3)
4113 fskyfi(nin)%P(4,niskyfil)= fics(j,4)
4114 iskyfi(nin)%P(niskyfil) = j1
4115 END DO
4116 ELSE
4117 niskyfil = niskyfil + 1
4118 fskyfi(nin)%P(1,niskyfil)=-fxi(i)
4119 fskyfi(nin)%P(2,niskyfil)=-fyi(i)
4120 fskyfi(nin)%P(3,niskyfil)=-fzi(i)
4121 fskyfi(nin)%P(4,niskyfil)= stif(i)
4122 iskyfi(nin)%P(niskyfil) = ig
4123 ENDIF
4124C must be be done for second/plyxfem node if it must be considered -----
4125 ENDIF
4126 ENDDO
4127 ENDIF
4128C Thermique
4129 ELSE
4130 IF(intply == 0) THEN
4131 DO i=1,jlt
4132 IF (h1(i)/=0.) THEN
4133 niskyl = niskyl + 1
4134 fskyi(niskyl,1)=fx1(i)
4135 fskyi(niskyl,2)=fy1(i)
4136 fskyi(niskyl,3)=fz1(i)
4137 fskyi(niskyl,4)=stif(i)*abs(h1(i))
4138 isky(niskyl) = ix1(i)
4139 ftheskyi(niskyl) = phi1(i)
4140 ENDIF
4141 ENDDO
4142 DO i=1,jlt
4143 IF (h2(i)/=zero) THEN
4144 niskyl = niskyl + 1
4145 fskyi(niskyl,1)=fx2(i)
4146 fskyi(niskyl,2)=fy2(i)
4147 fskyi(niskyl,3)=fz2(i)
4148 fskyi(niskyl,4)=stif(i)*abs(h2(i))
4149 isky(niskyl) = ix2(i)
4150 ftheskyi(niskyl) = phi2(i)
4151 ENDIF
4152 ENDDO
4153 DO i=1,jlt
4154 IF (h3(i)/=zero) THEN
4155 niskyl = niskyl + 1
4156 fskyi(niskyl,1)=fx3(i)
4157 fskyi(niskyl,2)=fy3(i)
4158 fskyi(niskyl,3)=fz3(i)
4159 fskyi(niskyl,4)=stif(i)*abs(h3(i))
4160 isky(niskyl) = ix3(i)
4161 ftheskyi(niskyl) = phi3(i)
4162 ENDIF
4163 ENDDO
4164 DO i=1,jlt
4165 IF (h4(i)/=zero) THEN
4166 niskyl = niskyl + 1
4167 fskyi(niskyl,1)=fx4(i)
4168 fskyi(niskyl,2)=fy4(i)
4169 fskyi(niskyl,3)=fz4(i)
4170 fskyi(niskyl,4)=stif(i)*abs(h4(i))
4171 isky(niskyl) = ix4(i)
4172 ftheskyi(niskyl) = phi4(i)
4173 ENDIF
4174 ENDDO
4175C
4176 DO i=1,jlt
4177 IF(hh(i)==zero)cycle
4178 ig =nsvg(i)
4179 IF(ig>0) THEN
4180 IF (ig >numnod) THEN
4181 jg = ig - numnod
4182 nfic=1
4183 fici(nfic)=fxi(i)
4184 nfic = nfic + 1
4185 fici(nfic)=fyi(i)
4186 nfic = nfic + 1
4187 fici(nfic)=fzi(i)
4188 nfic = nfic + 1
4189 fici(nfic)=stif(i)
4190 nfic = nfic + 1
4191 fici(nfic)=phi(i)
4192 CALL i24forc_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
4193 + jg ,nfic ,fici ,ficsth ,ixss )
4194 DO j = 1,4
4195 j1 = ixss(j)
4196 niskyl = niskyl + 1
4197 fskyi(niskyl,1)=-ficsth(j,1)
4198 fskyi(niskyl,2)=-ficsth(j,2)
4199 fskyi(niskyl,3)=-ficsth(j,3)
4200 fskyi(niskyl,4)= ficsth(j,4)
4201 ftheskyi(niskyl)=ficsth(j,5)
4202 isky(niskyl) = j1
4203 END DO
4204 ELSE
4205 niskyl = niskyl + 1
4206 fskyi(niskyl,1)=-fxi(i)
4207 fskyi(niskyl,2)=-fyi(i)
4208 fskyi(niskyl,3)=-fzi(i)
4209 fskyi(niskyl,4)= stif(i)
4210 isky(niskyl) = ig
4211 ftheskyi(niskyl)=phi(i)
4212 END IF !(IG >NUMNOD) THEN
4213 ELSE
4214 ig = -ig
4215 IF(isedge_fi(nin)%P(ig)==1)THEN
4216 jg = ig
4217 nfic=1
4218 fici(nfic)=fxi(i)
4219 nfic = nfic + 1
4220 fici(nfic)=fyi(i)
4221 nfic = nfic + 1
4222 fici(nfic)=fzi(i)
4223 nfic = nfic + 1
4224 fici(nfic)=stif(i)
4225 nfic = nfic + 1
4226 fici(nfic)=phi(i)
4227 CALL i24forc_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
4228 + jg ,nfic ,fici ,ficsth ,ixss )
4229
4230 DO j = 1,4
4231 j1 = ixss(j)
4232 niskyfil = niskyfil + 1
4233 fskyfi(nin)%P(1,niskyfil)=-ficsth(j,1)
4234 fskyfi(nin)%P(2,niskyfil)=-ficsth(j,2)
4235 fskyfi(nin)%P(3,niskyfil)=-ficsth(j,3)
4236 fskyfi(nin)%P(4,niskyfil)= ficsth(j,4)
4237 ftheskyfi(nin)%P(niskyfil)=ficsth(j,5)
4238 isky(niskyl) = j1
4239 END DO
4240 ELSE
4241 niskyfil = niskyfil + 1
4242 fskyfi(nin)%P(1,niskyfil)=-fxi(i)
4243 fskyfi(nin)%P(2,niskyfil)=-fyi(i)
4244 fskyfi(nin)%P(3,niskyfil)=-fzi(i)
4245 fskyfi(nin)%P(4,niskyfil)= stif(i)
4246 iskyfi(nin)%P(niskyfil) = ig
4247 ftheskyfi(nin)%P(niskyfil)=phi(i)
4248 ENDIF
4249 ENDIF
4250 ENDDO
4251C with plyxfem formulation
4252 ELSE
4253 DO i=1,jlt
4254 IF (h1(i)/=0.) THEN
4255 niskyl = niskyl + 1
4256 fskyi(niskyl,1)=fx1(i)
4257 fskyi(niskyl,2)=fy1(i)
4258 fskyi(niskyl,3)=fz1(i)
4259 fskyi(niskyl,4)=stif(i)*abs(h1(i))
4260 isky(niskyl) = ix1(i)
4261 ftheskyi(niskyl) = phi1(i)
4262C
4263 plyskyi%FSKYI(niskyl,1)=fx1(i)
4264 plyskyi%FSKYI(niskyl,2)=fy1(i)
4265 plyskyi%FSKYI(niskyl,3)=fz1(i)
4266 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h1(i))
4267 plyskyi%FSKYI(niskyl,5)=iply(1,i)
4268 ENDIF
4269 ENDDO
4270 DO i=1,jlt
4271 IF (h2(i)/=zero) THEN
4272 niskyl = niskyl + 1
4273 fskyi(niskyl,1)=fx2(i)
4274 fskyi(niskyl,2)=fy2(i)
4275 fskyi(niskyl,3)=fz2(i)
4276 fskyi(niskyl,4)=stif(i)*abs(h2(i))
4277 isky(niskyl) = ix2(i)
4278 ftheskyi(niskyl) = phi2(i)
4279C
4280 plyskyi%FSKYI(niskyl,1)=fx2(i)
4281 plyskyi%FSKYI(niskyl,2)=fy2(i)
4282 plyskyi%FSKYI(niskyl,3)=fz2(i)
4283 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h2(i))
4284 plyskyi%FSKYI(niskyl,5)=iply(2,i)
4285 ENDIF
4286 ENDDO
4287 DO i=1,jlt
4288 IF (h3(i)/=zero) THEN
4289 niskyl = niskyl + 1
4290 fskyi(niskyl,1)=fx3(i)
4291 fskyi(niskyl,2)=fy3(i)
4292 fskyi(niskyl,3)=fz3(i)
4293 fskyi(niskyl,4)=stif(i)*abs(h3(i))
4294 isky(niskyl) = ix3(i)
4295 ftheskyi(niskyl) = phi3(i)
4296C
4297 plyskyi%FSKYI(niskyl,1)=fx3(i)
4298 plyskyi%FSKYI(niskyl,2)=fy3(i)
4299 plyskyi%FSKYI(niskyl,3)=fz3(i)
4300 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h3(i))
4301 plyskyi%FSKYI(niskyl,5)=iply(3,i)
4302 ENDIF
4303 ENDDO
4304 DO i=1,jlt
4305 IF (h4(i)/=zero) THEN
4306 niskyl = niskyl + 1
4307 fskyi(niskyl,1)=fx4(i)
4308 fskyi(niskyl,2)=fy4(i)
4309 fskyi(niskyl,3)=fz4(i)
4310 fskyi(niskyl,4)=stif(i)*abs(h4(i))
4311 isky(niskyl) = ix4(i)
4312 ftheskyi(niskyl) = phi4(i)
4313C
4314 plyskyi%FSKYI(niskyl,1)=fx4(i)
4315 plyskyi%FSKYI(niskyl,2)=fy4(i)
4316 plyskyi%FSKYI(niskyl,3)=fz4(i)
4317 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h4(i))
4318 plyskyi%FSKYI(niskyl,5)=iply(4,i)
4319 ENDIF
4320 ENDDO
4321C
4322 DO i=1,jlt
4323 IF(hh(i)==zero)cycle
4324 ig =nsvg(i)
4325 IF(ig>0) THEN
4326 IF (ig >numnod) THEN
4327 jg = ig - numnod
4328 nfic=1
4329 fici(nfic)=fxi(i)
4330 nfic = nfic + 1
4331 fici(nfic)=fyi(i)
4332 nfic = nfic + 1
4333 fici(nfic)=fzi(i)
4334 nfic = nfic + 1
4335 fici(nfic)=stif(i)
4336 nfic = nfic + 1
4337 fici(nfic)=phi(i)
4338 CALL i24forc_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
4339 + jg ,nfic ,fici ,ficsth ,ixss )
4340 DO j = 1,4
4341 j1 = ixss(j)
4342 niskyl = niskyl + 1
4343 fskyi(niskyl,1)=-ficsth(j,1)
4344 fskyi(niskyl,2)=-ficsth(j,2)
4345 fskyi(niskyl,3)=-ficsth(j,3)
4346 fskyi(niskyl,4)= ficsth(j,4)
4347 ftheskyi(niskyl)=ficsth(j,5)
4348 isky(niskyl) = j1
4349 END DO
4350 ELSE
4351 niskyl = niskyl + 1
4352 fskyi(niskyl,1)=-fxi(i)
4353 fskyi(niskyl,2)=-fyi(i)
4354 fskyi(niskyl,3)=-fzi(i)
4355 fskyi(niskyl,4)= stif(i)
4356 isky(niskyl) = ig
4357 ftheskyi(niskyl)=phi(i)
4358 END IF !(IG >NUMNOD)
4359CC the second node is not plyxfem
4360cc PLYSKYI(1)%FSKYI(NISKYL,1)=-FXI(I)
4361cc PLYSKYI(1)%FSKYI(NISKYL,2)=-FYI(I)
4362cc PLYSKYI(1)%FSKYI(NISKYL,3)=-FZI(I)
4363cc PLYSKYI(I)%FSKYI(NISKYL,4)=STIF(I)
4364 ELSE
4365 ig = -ig
4366 IF(isedge_fi(nin)%P(ig)==1)THEN
4367 jg = ig
4368 nfic=1
4369 fici(nfic)=fxi(i)
4370 nfic = nfic + 1
4371 fici(nfic)=fyi(i)
4372 nfic = nfic + 1
4373 fici(nfic)=fzi(i)
4374 nfic = nfic + 1
4375 fici(nfic)=stif(i)
4376 nfic = nfic + 1
4377 fici(nfic)=phi(i)
4378 CALL i24forc_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
4379 + jg ,nfic ,fici ,ficsth ,ixss )
4380 DO j = 1,4
4381 j1 = ixss(j)
4382 niskyfil = niskyfil + 1
4383 fskyfi(nin)%P(1,niskyfil)=-ficsth(j,1)
4384 fskyfi(nin)%P(2,niskyfil)=-ficsth(j,2)
4385 fskyfi(nin)%P(3,niskyfil)=-ficsth(j,3)
4386 fskyfi(nin)%P(4,niskyfil)= ficsth(j,4)
4387 ftheskyfi(nin)%P(niskyfil)=ficsth(j,5)
4388 iskyfi(nin)%P(niskyfil) = j1
4389 END DO
4390
4391 ELSE
4392 niskyfil = niskyfil + 1
4393 fskyfi(nin)%P(1,niskyfil)=-fxi(i)
4394 fskyfi(nin)%P(2,niskyfil)=-fyi(i)
4395 fskyfi(nin)%P(3,niskyfil)=-fzi(i)
4396 fskyfi(nin)%P(4,niskyfil)= stif(i)
4397 iskyfi(nin)%P(niskyfil) = ig
4398 ftheskyfi(nin)%P(niskyfil)=phi(i)
4399 ENDIF
4400 ENDIF
4401 ENDDO
4402 ENDIF
4403 ENDIF
4404C
4405 RETURN
4406 END
4407!||====================================================================
4408!|| i24ass0 ../engine/source/interfaces/int24/i24for3.F
4409!||--- called by ------------------------------------------------------
4410!|| i24for3 ../engine/source/interfaces/int24/i24for3.F
4411!||--- calls -----------------------------------------------------
4412!|| i24forc_fic ../engine/source/interfaces/int24/i24for3e.F
4413!||--- uses -----------------------------------------------------
4414!|| plyxfem_mod ../engine/share/modules/plyxfem_mod.F
4415!|| tri7box ../engine/share/modules/tri7box.F
4416!||====================================================================
4417 SUBROUTINE i24ass0(JLT ,IX1 ,IX2 ,IX3 ,IX4 ,
4418 2 NSVG ,H1 ,H2 ,H3 ,H4 ,STIF ,
4419 3 FX1 ,FY1 ,FZ1 ,FX2 ,FY2 ,FZ2 ,
4420 4 FX3 ,FY3 ,FZ3 ,FX4 ,FY4 ,FZ4 ,
4421 5 FXI ,FYI ,FZI ,A ,STIFN ,NIN ,
4422 6 INTTH ,PHI ,FTHE ,PHI1 , PHI2 ,PHI3 ,
4423 7 PHI4 ,INTPLY, IPLY,INOD,IRTSE ,NSNE ,
4424 8 IS2SE ,IS2PT ,JTASK )
4425C-----------------------------------------------
4426C M o d u l e s
4427C-----------------------------------------------
4428 USE tri7box
4429 USE plyxfem_mod
4430C-----------------------------------------------
4431C I m p l i c i t T y p e s
4432C-----------------------------------------------
4433#include "implicit_f.inc"
4434C-----------------------------------------------
4435C G l o b a l P a r a m e t e r s
4436C-----------------------------------------------
4437#include "mvsiz_p.inc"
4438C-----------------------------------------------
4439C C o m m o n B l o c k s
4440C-----------------------------------------------
4441#include "com04_c.inc"
4442C-----------------------------------------------
4443C D u m m y A r g u m e n t s
4444C-----------------------------------------------
4445 INTEGER JLT, NIN,INTTH,
4446 . ix1(mvsiz),ix2(mvsiz),ix3(mvsiz),ix4(mvsiz),nsvg(mvsiz),
4447 . intply,iply(4,*),inod(*),irtse(5,*),nsne,is2se(2,*),
4448 . is2pt(*) ,jtask
4449 my_real
4450 . h1(mvsiz),h2(mvsiz),h3(mvsiz),h4(mvsiz),stif(mvsiz),
4451 . fx1(mvsiz),fy1(mvsiz),fz1(mvsiz),
4452 . fx2(mvsiz),fy2(mvsiz),fz2(mvsiz),
4453 . fx3(mvsiz),fy3(mvsiz),fz3(mvsiz),
4454 . fx4(mvsiz),fy4(mvsiz),fz4(mvsiz),
4455 . fxi(mvsiz),fyi(mvsiz),fzi(mvsiz),
4456 . a(3,*), stifn(*),phi(*), fthe(*),
4457 . phi1(*), phi2(*), phi3(*), phi4(*)
4458C-----------------------------------------------
4459C L o c a l V a r i a b l e s
4460C-----------------------------------------------
4461 my_real
4462 . hh(mvsiz),fici(5),fics(4,4),ficsth(4,5)
4463 INTEGER I, J1, IG,ILY,NN,JG,IXSS(4),NFIC,J,ISHIFT,NODFI
4464
4465C
4466 IF(INTTH == 0) then
4467 DO i=1,jlt
4468 hh(i)=h1(i)+h2(i)+h3(i)+h4(i)
4469 IF(hh(i)==zero)cycle
4470 j1=ix1(i)
4471 a(1,j1)=a(1,j1)+fx1(i)
4472 a(2,j1)=a(2,j1)+fy1(i)
4473 a(3,j1)=a(3,j1)+fz1(i)
4474 stifn(j1) = stifn(j1) + stif(i)*abs(h1(i))
4475C
4476 j1=ix2(i)
4477 a(1,j1)=a(1,j1)+fx2(i)
4478 a(2,j1)=a(2,j1)+fy2(i)
4479 a(3,j1)=a(3,j1)+fz2(i)
4480 stifn(j1) = stifn(j1) + stif(i)*abs(h2(i))
4481C
4482 j1=ix3(i)
4483 a(1,j1)=a(1,j1)+fx3(i)
4484 a(2,j1)=a(2,j1)+fy3(i)
4485 a(3,j1)=a(3,j1)+fz3(i)
4486 stifn(j1) = stifn(j1) + stif(i)*abs(h3(i))
4487C
4488 j1=ix4(i)
4489 a(1,j1)=a(1,j1)+fx4(i)
4490 a(2,j1)=a(2,j1)+fy4(i)
4491 a(3,j1)=a(3,j1)+fz4(i)
4492 stifn(j1) = stifn(j1) + stif(i)*abs(h4(i))
4493 ENDDO
4494 ELSE
4495 DO i=1,jlt
4496 hh(i)=h1(i)+h2(i)+h3(i)+h4(i)
4497 IF(hh(i)==zero)cycle
4498 j1=ix1(i)
4499 a(1,j1)=a(1,j1)+fx1(i)
4500 a(2,j1)=a(2,j1)+fy1(i)
4501 a(3,j1)=a(3,j1)+fz1(i)
4502 stifn(j1) = stifn(j1) + stif(i)*abs(h1(i))
4503 fthe(j1) = fthe(j1) + phi1(i)
4504C
4505 j1=ix2(i)
4506 a(1,j1)=a(1,j1)+fx2(i)
4507 a(2,j1)=a(2,j1)+fy2(i)
4508 a(3,j1)=a(3,j1)+fz2(i)
4509 stifn(j1) = stifn(j1) + stif(i)*abs(h2(i))
4510 fthe(j1) = fthe(j1) + phi2(i)
4511C
4512 j1=ix3(i)
4513 a(1,j1)=a(1,j1)+fx3(i)
4514 a(2,j1)=a(2,j1)+fy3(i)
4515 a(3,j1)=a(3,j1)+fz3(i)
4516 stifn(j1) = stifn(j1) + stif(i)*abs(h3(i))
4517 fthe(j1) = fthe(j1) + phi3(i)
4518C
4519 j1=ix4(i)
4520 a(1,j1)=a(1,j1)+fx4(i)
4521 a(2,j1)=a(2,j1)+fy4(i)
4522 a(3,j1)=a(3,j1)+fz4(i)
4523 stifn(j1) = stifn(j1) + stif(i)*abs(h4(i))
4524 fthe(j1) = fthe(j1) + phi4(i)
4525 ENDDO
4526 ENDIF
4527 IF(intply > 0) THEN
4528 DO i=1,jlt
4529 hh(i)=h1(i)+h2(i)+h3(i)+h4(i)
4530 IF(hh(i)==zero)cycle
4531 j1=ix1(i)
4532 nn = inod(j1)
4533 ily = iply(1,i)
4534 IF(inod(j1) > 0 .AND. ily > 0 ) THEN
4535 ply(ily)%A(1,nn)=ply(ily)%A(1,nn) + fx1(i)
4536 ply(ily)%A(2,nn)=ply(ily)%A(2,nn) + fy1(i)
4537 ply(ily)%A(3,nn)=ply(ily)%A(3,nn) + fz1(i)
4538 ply(ily)%A(4,nn)=ply(ily)%A(4,nn) + stif(i)*abs(h1(i))
4539 ENDIF
4540C
4541 j1=ix2(i)
4542 nn = inod(j1)
4543 ily = iply(2,i)
4544 IF(inod(j1) > 0 .AND. ily > 0) THEN
4545 ply(ily)%A(1,nn)=ply(ily)%A(1,nn) +fx2(i)
4546 ply(ily)%A(2,nn)=ply(ily)%A(2,nn) +fy2(i)
4547 ply(ily)%A(3,nn)=ply(ily)%A(3,nn) +fz2(i)
4548 ply(ily)%A(4,nn)=ply(ily)%A(4,nn) + stif(i)*abs(h2(i))
4549 ENDIF
4550C
4551 j1=ix3(i)
4552 nn = inod(j1)
4553 ily = iply(3,i)
4554 IF(inod(j1) > 0 .AND. ily > 0) THEN
4555 ply(ily)%A(1,nn)=ply(ily)%A(1,nn) +fx3(i)
4556 ply(ily)%A(2,nn)=ply(ily)%A(2,nn) +fy3(i)
4557 ply(ily)%A(3,nn)=ply(ily)%A(3,nn) +fz3(i)
4558 ply(ily)%A(4,nn)=ply(ily)%A(4,nn) + stif(i)*abs(h3(i))
4559 ENDIF
4560C
4561 j1=ix4(i)
4562 nn = inod(j1)
4563 ily = iply(4,i)
4564 IF(inod(j1) > 0 .AND. ily > 0) THEN
4565 ply(ily)%A(1,nn)=ply(ily)%A(1,nn) +fx4(i)
4566 ply(ily)%A(2,nn)=ply(ily)%A(2,nn) +fy4(i)
4567 ply(ily)%A(3,nn)=ply(ily)%A(3,nn) +fz4(i)
4568 ply(ily)%A(4,nn)=ply(ily)%A(4,nn) + stif(i)*abs(h4(i))
4569 ENDIF
4570 ENDDO
4571 ENDIF
4572C
4573 nodfi = nlskyfi(nin)
4574 ishift = nodfi*(jtask-1)
4575 IF(intth == 0 ) THEN
4576 DO i=1,jlt
4577 IF(hh(i)==zero)cycle
4578 ig=nsvg(i)
4579 IF(ig>0)THEN
4580 IF (ig >numnod) THEN
4581 jg = ig - numnod
4582 nfic=1
4583 fici(nfic)=fxi(i)
4584 nfic = nfic + 1
4585 fici(nfic)=fyi(i)
4586 nfic = nfic + 1
4587 fici(nfic)=fzi(i)
4588 nfic = nfic + 1
4589 fici(nfic)=stif(i)
4590 CALL i24forc_fic( 3,irtse ,nsne ,is2se ,is2pt ,
4591 + jg ,nfic ,fici ,fics ,ixss )
4592 DO j = 1,4
4593 j1 = ixss(j)
4594 a(1,j1)=a(1,j1)-fics(j,1)
4595 a(2,j1)=a(2,j1)-fics(j,2)
4596 a(3,j1)=a(3,j1)-fics(j,3)
4597 stifn(j1) = stifn(j1) + fics(j,4)
4598 END DO
4599 ELSE
4600 a(1,ig)=a(1,ig)-fxi(i)
4601 a(2,ig)=a(2,ig)-fyi(i)
4602 a(3,ig)=a(3,ig)-fzi(i)
4603 stifn(ig) = stifn(ig) + stif(i)
4604 END IF !(IG >NUMNOD) THEN
4605 ELSE
4606 ig = -ig
4607 IF(isedge_fi(nin)%P(ig)==1)THEN
4608 jg = ig
4609 nfic=1
4610 fici(nfic)=fxi(i)
4611 nfic = nfic + 1
4612 fici(nfic)=fyi(i)
4613 nfic = nfic + 1
4614 fici(nfic)=fzi(i)
4615 nfic = nfic + 1
4616 fici(nfic)=stif(i)
4617 CALL i24forc_fic( 3,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
4618 + jg ,nfic ,fici ,fics ,ixss )
4619 DO j = 1,4
4620 j1 = ixss(j)
4621 afi(nin)%P(1,j1+ishift)=afi(nin)%P(1,j1+ishift)-fics(j,1)
4622 afi(nin)%P(2,j1+ishift)=afi(nin)%P(2,j1+ishift)-fics(j,2)
4623 afi(nin)%P(3,j1+ishift)=afi(nin)%P(3,j1+ishift)-fics(j,3)
4624 stnfi(nin)%P(j1+ishift)=stnfi(nin)%P(j1+ishift) + fics(j,4)
4625 END DO
4626
4627 ELSE
4628 afi(nin)%P(1,ig+ishift)=afi(nin)%P(1,ig+ishift)-fxi(i)
4629 afi(nin)%P(2,ig+ishift)=afi(nin)%P(2,ig+ishift)-fyi(i)
4630 afi(nin)%P(3,ig+ishift)=afi(nin)%P(3,ig+ishift)-fzi(i)
4631 stnfi(nin)%P(ig+ishift)=stnfi(nin)%P(ig+ishift) + stif(i)
4632 ENDIF
4633 ENDIF
4634 ENDDO
4635C
4636 ELSE
4637 DO i=1,jlt
4638 IF(hh(i)==zero)cycle
4639 ig=nsvg(i)
4640 IF(ig>0)THEN
4641 IF (ig >numnod) THEN
4642 jg = ig - numnod
4643 nfic=1
4644 fici(nfic)=fxi(i)
4645 nfic = nfic + 1
4646 fici(nfic)=fyi(i)
4647 nfic = nfic + 1
4648 fici(nfic)=fzi(i)
4649 nfic = nfic + 1
4650 fici(nfic)=stif(i)
4651 nfic = nfic + 1
4652 fici(nfic)=phi(i)
4653 CALL i24forc_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
4654 + jg ,nfic ,fici ,ficsth ,ixss )
4655 DO j = 1,4
4656 j1 = ixss(j)
4657 a(1,j1)=a(1,j1)-ficsth(j,1)
4658 a(2,j1)=a(2,j1)-ficsth(j,2)
4659 a(3,j1)=a(3,j1)-ficsth(j,3)
4660 stifn(j1) = stifn(j1) + ficsth(j,4)
4661 fthe(j1)=fthe(j1) + ficsth(j,5)
4662 END DO
4663 ELSE
4664 a(1,ig)=a(1,ig)-fxi(i)
4665 a(2,ig)=a(2,ig)-fyi(i)
4666 a(3,ig)=a(3,ig)-fzi(i)
4667 stifn(ig) = stifn(ig) + stif(i)
4668 fthe(ig)=fthe(ig) + phi(i)
4669 END IF !(IG >NUMNOD) THEN
4670 ELSE
4671 ig = -ig
4672 IF(isedge_fi(nin)%P(ig)==1)THEN
4673 jg = ig
4674 nfic=1
4675 fici(nfic)=fxi(i)
4676 nfic = nfic + 1
4677 fici(nfic)=fyi(i)
4678 nfic = nfic + 1
4679 fici(nfic)=fzi(i)
4680 nfic = nfic + 1
4681 fici(nfic)=stif(i)
4682 nfic = nfic + 1
4683 fici(nfic)=phi(i)
4684 CALL i24forc_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
4685 + jg ,nfic ,fici ,ficsth ,ixss )
4686 DO j = 1,4
4687 j1 = ixss(j)
4688 afi(nin)%P(1,j1+ishift)=afi(nin)%P(1,j1+ishift)-ficsth(j,1)
4689 afi(nin)%P(2,j1+ishift)=afi(nin)%P(2,j1+ishift)-ficsth(j,2)
4690 afi(nin)%P(3,j1+ishift)=afi(nin)%P(3,j1+ishift)-ficsth(j,3)
4691
4692 stnfi(nin)%P(j1+ishift)=stnfi(nin)%P(j1+ishift) + ficsth(j,4)
4693 fthefi(nin)%P(j1+ishift)= fthefi(nin)%P(j1+ishift) + ficsth(j,5)
4694 END DO
4695
4696 ELSE
4697 afi(nin)%P(1,ig+ishift)=afi(nin)%P(1,ig+ishift)-fxi(i)
4698 afi(nin)%P(2,ig+ishift)=afi(nin)%P(2,ig+ishift)-fyi(i)
4699 afi(nin)%P(3,ig+ishift)=afi(nin)%P(3,ig+ishift)-fzi(i)
4700 stnfi(nin)%P(ig+ishift)=stnfi(nin)%P(ig+ishift) + stif(i)
4701 fthefi(nin)%P(ig+ishift)= fthefi(nin)%P(ig+ishift) + phi(i)
4702 ENDIF
4703 ENDIF
4704 ENDDO
4705 ENDIF
4706
4707cc IF(INTPLY > 0) THEN
4708cc DO I=1,JLT
4709cc IF(HH(I)==ZERO)CYCLE
4710cc IG=NSVG(I)
4711cc NN = INOD(IG)
4712cc IF(IG>0 .AND. NN > 0)THEN
4713cc ILY = IPLY(5,I)
4714cc PLY(ILY)%A(1,NN) = PLY(ILY)%A(1,NN) - FXI(I)
4715cc PLY(ILY)%A(2,NN) = PLY(ILY)%A(2,NN) - FYI(I)
4716cc PLY(ILY)%A(3,NN) = PLY(ILY)%A(3,NN) - FZI(I)
4717cc PLY(ILY)%A(4,NN) = PLY(ILY)%A(4,NN) + STIF(I)
4718cc ELSE
4719C have the SPMD appearance
4720cc IG = -IG
4721cc AFI(NIN)%P(1,IG)=AFI(NIN)%P(1,IG)-FXI(I)
4722cc AFI(NIN)%P(2,IG)=AFI(NIN)%P(2,IG)-FYI(I)
4723cc AFI(NIN)%P(3,IG)=AFI(NIN)%P(3,IG)-FZI(I)
4724cc ENDIF
4725cc ENDDO
4726cc ENDIF
4727C
4728 RETURN
4729 END
4730C
4731!||====================================================================
4732!|| i24sms2 ../engine/source/interfaces/int24/i24for3.F
4733!||--- called by ------------------------------------------------------
4734!|| i24for3 ../engine/source/interfaces/int24/i24for3.F
4735!||--- calls -----------------------------------------------------
4736!|| ancmsg ../engine/source/output/message/message.F
4737!|| arret ../engine/source/system/arret.F
4738!|| i24forc_fic ../engine/source/interfaces/int24/i24for3e.f
4739!||--- uses -----------------------------------------------------
4740!|| message_mod ../engine/share/message_module/message_mod.F
4741!|| tri7box ../engine/share/modules/tri7box.F
4742!||====================================================================
4743 SUBROUTINE i24sms2(JLT ,IX1 ,IX2 ,IX3 ,IX4 ,
4744 2 NSVG ,H1 ,H2 ,H3 ,H4 ,STIF ,
4745 3 NIN ,NOINT ,MSKYI_SMS ,ISKYI_SMS,NSMS ,
4746 4 KT ,C ,CF ,DTMINI,DTI ,
4747 5 IRTSE ,NSNE ,IS2SE,IS2PT ,T2MAIN_SMS,
4748 6 T2FAC_SMS)
4749C-----------------------------------------------
4750C M o d u l e s
4751C-----------------------------------------------
4752 USE tri7box
4753 USE message_mod
4754C-----------------------------------------------
4755C I m p l i c i t T y p e s
4756C-----------------------------------------------
4757#include "implicit_f.inc"
4758#include "comlock.inc"
4759C-----------------------------------------------
4760C G l o b a l P a r a m e t e r s
4761C-----------------------------------------------
4762#include "mvsiz_p.inc"
4763C-----------------------------------------------
4764C C o m m o n B l o c k s
4765C-----------------------------------------------
4766#include "parit_c.inc"
4767#include "task_c.inc"
4768#include "sms_c.inc"
4769#include "com04_c.inc"
4770C-----------------------------------------------
4771C D u m m y A r g u m e n t s
4772C-----------------------------------------------
4773 INTEGER JLT,NIN,NOINT,
4774 . IX1(MVSIZ),IX2(MVSIZ),IX3(MVSIZ),IX4(MVSIZ),NSVG(MVSIZ),
4775 . NSMS(*), ISKYI_SMS(LSKYI_SMS,*),IRTSE(5,*),NSNE,IS2SE(2,*),
4776 . IS2PT(*),T2MAIN_SMS(6,*)
4777 my_real
4778 . H1(MVSIZ),H2(MVSIZ),H3(MVSIZ),H4(MVSIZ),STIF(MVSIZ),
4779 . MSKYI_SMS(*), KT(MVSIZ), C(MVSIZ), CF(MVSIZ), DTMINI, DTI,T2FAC_SMS(*)
4780C-----------------------------------------------
4781C L o c a l V a r i a b l e s
4782C-----------------------------------------------
4783 INTEGER I, IG, NISKYL1, NISKYL, NN, IXSS(4), J, NS, NL, K, KK
4784 my_real
4785 . hh(mvsiz),mas, dts,dtm_int,fici,fics(4),fac
4786
4787C-----------------------------------------------
4788C IF contact on nodes of TYPE2 additional connections must be created (T2MAIN_SMS(1) > 1)
4789C-----------------------------------------------
4790C
4791 niskyl1 = 0
4792C
4793 DO i=1,jlt
4794
4795 hh(i)=h1(i)+h2(i)+h3(i)+h4(i)
4796 IF(nsms(i)==0.OR.stif(i)==zero.OR.hh(i)==zero) cycle
4797C
4798 ig =nsvg(i)
4799 IF(ig > 0) THEN
4800 IF (ig > numnod) THEN
4801 nl = 0
4802 ns = ig - numnod
4803 fici = one
4804 CALL i24forc_fic(3,irtse,nsne,is2se,is2pt,ns,1,fici,fics,ixss)
4805 DO j=1,4
4806 nl = nl + t2main_sms(1,ixss(j))*ceiling(fics(j))
4807 ENDDO
4808 ELSE
4809 nl = t2main_sms(1,ig)
4810 ENDIF
4811 ELSE
4812 IF(isedge_fi(nin)%P(-ig)==1) THEN
4813 nl = 0
4814 fici = one
4815 CALL i24forc_fic(3,irtse_fi(nin)%P(1,1),nsne,is2se_fi(nin)%P(1,1), is2pt_fi(nin)%P(1),-ig,1,fici,fics,ixss)
4816 DO j=1,4
4817 nl = nl + t2main_sms_fi(nin)%P(1,ixss(j))*ceiling(fics(j))
4818 ENDDO
4819 ELSE
4820 nl = t2main_sms_fi(nin)%P(1,-ig)
4821 ENDIF
4822 ENDIF
4823C
4824 IF (h1(i)/=zero) niskyl1 = niskyl1 + nl*t2main_sms(1,ix1(i))
4825 IF (h2(i)/=zero) niskyl1 = niskyl1 + nl*t2main_sms(1,ix2(i))
4826 IF (h3(i)/=zero) niskyl1 = niskyl1 + nl*t2main_sms(1,ix3(i))
4827 IF (h4(i)/=zero) niskyl1 = niskyl1 + nl*t2main_sms(1,ix4(i))
4828C
4829 ENDDO
4830C
4831#include "lockon.inc"
4832 niskyl = nisky_sms
4833 nisky_sms = nisky_sms + niskyl1
4834#include "lockoff.inc"
4835C
4836 IF (niskyl+niskyl1 > lskyi_sms) THEN
4837 CALL ancmsg(msgid=26,anmode=aninfo)
4838 CALL arret(2)
4839 ENDIF
4840C
4841 IF (dtmini>zero) THEN
4842 dtm_int=dtmini
4843 ELSE
4844 dtm_int=dtmins_int
4845 ENDIF
4846
4847 DO i=1,jlt
4848 IF(nsms(i)==0.OR.stif(i)==zero.OR.hh(i)==zero) cycle
4849C
4850 IF(nsms(i)>0)THEN
4851 dts = dtmins/dtfacs
4852 dti=min(dti,dtmins)
4853 ELSE
4854 dts = dtm_int/dtfacs_int
4855 dti=min(dti,dtm_int)
4856 END IF
4857
4858 mas= dts * ( dts * kt(i) + c(i) )
4859 mas = half * max( mas, dts * cf(i) )
4860C
4861C-----------------------------------------------
4862C FAC = 1.0 for all cases
4863C FAC = 0.5 for type2 crossed connections - T2MAIN(i)*T2MAIN(j)=16
4864C-----------------------------------------------
4865C
4866 ig =nsvg(i)
4867 IF(ig > 0)THEN
4868 IF (ig > numnod) THEN
4869C-----------------------------------------------
4870C-- Edge 2 Edge contact
4871C-----------------------------------------------
4872 ns = ig - numnod
4873 fici = mas
4874 CALL i24forc_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
4875 + ns ,1 ,fici ,fics ,ixss )
4876C
4877 DO j = 1, 4
4878 IF (fics(j) > zero) THEN
4879 IF(h1(i)/=zero)THEN
4880 fac = t2fac_sms(ixss(j))
4881 IF (t2main_sms(1,ixss(j))*t2main_sms(1,ix1(i))==16) fac = half*max(t2fac_sms(ixss(j)),t2fac_sms(ix1(i)))
4882 DO k=1,t2main_sms(1,ixss(j))
4883 DO kk=1,t2main_sms(1,ix1(i))
4884 niskyl=niskyl+1
4885 mskyi_sms(niskyl)=abs(h1(i))*fics(j)*fac
4886 iskyi_sms(niskyl,1)=t2main_sms(k+1,ixss(j))
4887 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix1(i))
4888 iskyi_sms(niskyl,3)=ispmd+1
4889 ENDDO
4890 ENDDO
4891 ENDIF
4892 IF(h2(i)/=zero)THEN
4893 fac = t2fac_sms(ixss(j))
4894 IF (t2main_sms(1,ixss(j))*t2main_sms(1,ix2(i))==16) fac = half*max(t2fac_sms(ixss(j)),t2fac_sms(ix2(i)))
4895 DO k=1,t2main_sms(1,ixss(j))
4896 DO kk=1,t2main_sms(1,ix2(i))
4897 niskyl=niskyl+1
4898 mskyi_sms(niskyl)=abs(h2(i))*fics(j)*fac
4899 iskyi_sms(niskyl,1)=t2main_sms(k+1,ixss(j))
4900 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix2(i))
4901 iskyi_sms(niskyl,3)=ispmd+1
4902 ENDDO
4903 ENDDO
4904 ENDIF
4905 IF(h3(i)/=zero)THEN
4906 fac = t2fac_sms(ixss(j))
4907 IF (t2main_sms(1,ixss(j))*t2main_sms(1,ix3(i))==16) fac = half*max(t2fac_sms(ixss(j)),t2fac_sms(ix3(i)))
4908 DO k=1,t2main_sms(1,ixss(j))
4909 DO kk=1,t2main_sms(1,ix3(i))
4910 niskyl=niskyl+1
4911 mskyi_sms(niskyl)=abs(h3(i))*fics(j)*fac
4912 iskyi_sms(niskyl,1)=t2main_sms(k+1,ixss(j))
4913 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix3(i))
4914 iskyi_sms(niskyl,3)=ispmd+1
4915 ENDDO
4916 ENDDO
4917 ENDIF
4918 IF(h4(i)/=zero)THEN
4919 fac = t2fac_sms(ixss(j))
4920 IF (t2main_sms(1,ixss(j))*t2main_sms(1,ix4(i))==16) fac = half*max(t2fac_sms(ixss(j)),t2fac_sms(ix4(i)))
4921 DO k=1,t2main_sms(1,ixss(j))
4922 DO kk=1,t2main_sms(1,ix4(i))
4923 niskyl=niskyl+1
4924 mskyi_sms(niskyl)=abs(h4(i))*fics(j)*fac
4925 iskyi_sms(niskyl,1)=t2main_sms(k+1,ixss(j))
4926 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix4(i))
4927 iskyi_sms(niskyl,3)=ispmd+1
4928 ENDDO
4929 ENDDO
4930 ENDIF
4931 ENDIF
4932 END DO
4933 ELSE
4934C-----------------------------------------------
4935C-- Main / node contact
4936C-----------------------------------------------
4937 IF(h1(i)/=zero)THEN
4938 fac = t2fac_sms(ig)
4939 IF (t2main_sms(1,ig)*t2main_sms(1,ix1(i))==16) fac = half*max(t2fac_sms(ig),t2fac_sms(ix1(i)))
4940 DO k=1,t2main_sms(1,ig)
4941 DO kk=1,t2main_sms(1,ix1(i))
4942 niskyl=niskyl+1
4943 mskyi_sms(niskyl)=abs(h1(i))*mas*fac
4944 iskyi_sms(niskyl,1)=t2main_sms(k+1,ig)
4945 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix1(i))
4946 iskyi_sms(niskyl,3)=ispmd+1
4947 ENDDO
4948 ENDDO
4949 END IF
4950 IF(h2(i)/=zero)THEN
4951 fac = t2fac_sms(ig)
4952 IF (t2main_sms(1,ig)*t2main_sms(1,ix2(i))==16) fac = half*max(t2fac_sms(ig),t2fac_sms(ix2(i)))
4953 DO k=1,t2main_sms(1,ig)
4954 DO kk=1,t2main_sms(1,ix2(i))
4955 niskyl=niskyl+1
4956 mskyi_sms(niskyl)=abs(h2(i))*mas*fac
4957 iskyi_sms(niskyl,1)=t2main_sms(k+1,ig)
4958 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix2(i))
4959 iskyi_sms(niskyl,3)=ispmd+1
4960 ENDDO
4961 ENDDO
4962 END IF
4963 IF(h3(i)/=zero)THEN
4964 fac = t2fac_sms(ig)
4965 IF (t2main_sms(1,ig)*t2main_sms(1,ix3(i))==16) fac = half*max(t2fac_sms(ig),t2fac_sms(ix3(i)))
4966 DO k=1,t2main_sms(1,ig)
4967 DO kk=1,t2main_sms(1,ix3(i))
4968 niskyl=niskyl+1
4969 mskyi_sms(niskyl)=abs(h3(i))*mas*fac
4970 iskyi_sms(niskyl,1)=t2main_sms(k+1,ig)
4971 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix3(i))
4972 iskyi_sms(niskyl,3)=ispmd+1
4973 ENDDO
4974 ENDDO
4975 END IF
4976 IF(h4(i)/=zero)THEN
4977 fac = t2fac_sms(ig)
4978 IF (t2main_sms(1,ig)*t2main_sms(1,ix4(i))==16) fac = half*max(t2fac_sms(ig),t2fac_sms(ix4(i)))
4979 DO k=1,t2main_sms(1,ig)
4980 DO kk=1,t2main_sms(1,ix4(i))
4981 niskyl=niskyl+1
4982 mskyi_sms(niskyl)=abs(h4(i))*mas*fac
4983 iskyi_sms(niskyl,1)=t2main_sms(k+1,ig)
4984 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix4(i))
4985 iskyi_sms(niskyl,3)=ispmd+1
4986 ENDDO
4987 ENDDO
4988 END IF
4989 END IF !(IG > NUMNOD) THNE
4990C
4991C-----------------------------------------------
4992 ELSE
4993C-----------------------------------------------
4994C
4995 nn = -ig
4996 IF (isedge_fi(nin)%P(nn)==1)THEN
4997C-----------------------------------------------
4998C-- E2E remote contact
4999C-----------------------------------------------
5000 fici = mas
5001 CALL i24forc_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) , is2pt_fi(nin)%P(1) ,
5002 + nn ,1 ,fici ,fics ,ixss )
5003C
5004 DO j = 1, 4
5005 IF (fics(j) > zero) THEN
5006 IF(h1(i)/=zero)THEN
5007 fac = t2fac_sms_fi(nin)%P(ixss(j))
5008 IF (t2main_sms_fi(nin)%P(1,ixss(j))*t2main_sms(1,ix1(i))==16)
5009 . fac = half*max(t2fac_sms_fi(nin)%P(ixss(j)),t2fac_sms(ix1(i)))
5010 DO k=1,t2main_sms_fi(nin)%P(1,ixss(j))
5011 DO kk=1,t2main_sms(1,ix1(i))
5012 niskyl=niskyl+1
5013 mskyi_sms(niskyl)=abs(h1(i))*fics(j)*fac
5014 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,ixss(j))
5015 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix1(i))
5016 iskyi_sms(niskyl,3)=procamsfi(nin)%P(ixss(j))
5017 ENDDO
5018 ENDDO
5019 END IF
5020 IF(h2(i)/=zero)THEN
5021 fac = t2fac_sms_fi(nin)%P(ixss(j))
5022 IF (t2main_sms_fi(nin)%P(1,ixss(j))*t2main_sms(1,ix2(i))==16)
5023 . fac = half*max(t2fac_sms_fi(nin)%P(ixss(j)),t2fac_sms(ix2(i)))
5024 DO k=1,t2main_sms_fi(nin)%P(1,ixss(j))
5025 DO kk=1,t2main_sms(1,ix2(i))
5026 niskyl=niskyl+1
5027 mskyi_sms(niskyl)=abs(h2(i))*fics(j)*fac
5028 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,ixss(j))
5029 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix2(i))
5030 iskyi_sms(niskyl,3)=procamsfi(nin)%P(ixss(j))
5031 ENDDO
5032 ENDDO
5033 END IF
5034 IF(h3(i)/=zero)THEN
5035 fac = t2fac_sms_fi(nin)%P(ixss(j))
5036 IF (t2main_sms_fi(nin)%P(1,ixss(j))*t2main_sms(1,ix3(i))==16)
5037 . fac = half*max(t2fac_sms_fi(nin)%P(ixss(j)),t2fac_sms(ix3(i)))
5038 DO k=1,t2main_sms_fi(nin)%P(1,ixss(j))
5039 DO kk=1,t2main_sms(1,ix3(i))
5040 niskyl=niskyl+1
5041 mskyi_sms(niskyl)=abs(h3(i))*fics(j)*fac
5042 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,ixss(j))
5043 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix3(i))
5044 iskyi_sms(niskyl,3)=procamsfi(nin)%P(ixss(j))
5045 ENDDO
5046 ENDDO
5047 END IF
5048 IF(h4(i)/=zero)THEN
5049 fac = t2fac_sms_fi(nin)%P(ixss(j))
5050 IF (t2main_sms_fi(nin)%P(1,ixss(j))*t2main_sms(1,ix4(i))==16)
5051 . fac = half*max(t2fac_sms_fi(nin)%P(ixss(j)),t2fac_sms(ix4(i)))
5052 DO k=1,t2main_sms_fi(nin)%P(1,ixss(j))
5053 DO kk=1,t2main_sms(1,ix4(i))
5054 niskyl=niskyl+1
5055 mskyi_sms(niskyl)=abs(h4(i))*fics(j)*fac
5056 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,ixss(j))
5057 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix4(i))
5058 iskyi_sms(niskyl,3)=procamsfi(nin)%P(ixss(j))
5059 ENDDO
5060 ENDDO
5061 END IF
5062 ENDIF
5063 ENDDO
5064 ELSE
5065C-----------------------------------------------
5066C-- Main / node remote contact
5067C-----------------------------------------------
5068 IF(h1(i)/=zero)THEN
5069 fac = t2fac_sms_fi(nin)%P(nn)
5070 IF (t2main_sms_fi(nin)%P(1,nn)*t2main_sms(1,ix1(i))==16)
5071 . fac = half*max(t2fac_sms_fi(nin)%P(nn),t2fac_sms(ix1(i)))
5072 DO k=1,t2main_sms_fi(nin)%P(1,nn)
5073 DO kk=1,t2main_sms(1,ix1(i))
5074 niskyl=niskyl+1
5075 mskyi_sms(niskyl)=abs(h1(i))*mas*fac
5076 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,nn)
5077 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix1(i))
5078 iskyi_sms(niskyl,3)=procamsfi(nin)%P(nn)
5079 ENDDO
5080 ENDDO
5081 END IF
5082 IF(h2(i)/=zero)THEN
5083 fac = t2fac_sms_fi(nin)%P(nn)
5084 IF (t2main_sms_fi(nin)%P(1,nn)*t2main_sms(1,ix2(i))==16)
5085 . fac = half*max(t2fac_sms_fi(nin)%P(nn),t2fac_sms(ix2(i)))
5086 DO k=1,t2main_sms_fi(nin)%P(1,nn)
5087 DO kk=1,t2main_sms(1,ix2(i))
5088 niskyl=niskyl+1
5089 mskyi_sms(niskyl)=abs(h2(i))*mas*fac
5090 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,nn)
5091 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix2(i))
5092 iskyi_sms(niskyl,3)=procamsfi(nin)%P(nn)
5093 ENDDO
5094 ENDDO
5095 END IF
5096 IF(h3(i)/=zero)THEN
5097 fac = t2fac_sms_fi(nin)%P(nn)
5098 IF (t2main_sms_fi(nin)%P(1,nn)*t2main_sms(1,ix3(i))==16)
5099 . fac = half*max(t2fac_sms_fi(nin)%P(nn),t2fac_sms(ix3(i)))
5100 DO k=1,t2main_sms_fi(nin)%P(1,nn)
5101 DO kk=1,t2main_sms(1,ix3(i))
5102 niskyl=niskyl+1
5103 mskyi_sms(niskyl)=abs(h3(i))*mas*fac
5104 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,nn)
5105 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix3(i))
5106 iskyi_sms(niskyl,3)=procamsfi(nin)%P(nn)
5107 ENDDO
5108 ENDDO
5109 END IF
5110 IF(h4(i)/=zero)THEN
5111 fac = t2fac_sms_fi(nin)%P(nn)
5112 IF (t2main_sms_fi(nin)%P(1,nn)*t2main_sms(1,ix4(i))==16)
5113 . fac = half*max(t2fac_sms_fi(nin)%P(nn),t2fac_sms(ix4(i)))
5114 DO k=1,t2main_sms_fi(nin)%P(1,nn)
5115 DO kk=1,t2main_sms(1,ix4(i))
5116 niskyl=niskyl+1
5117 mskyi_sms(niskyl)=abs(h4(i))*mas*fac
5118 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,nn)
5119 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix4(i))
5120 iskyi_sms(niskyl,3)=procamsfi(nin)%P(nn)
5121 ENDDO
5122 ENDDO
5123 END IF
5124 ENDIF
5125 END IF
5126 ENDDO
5127C
5128 RETURN
5129 END
5130C
5131
my_real
#define my_real
Definition cppsort.cpp:32
i24_save_sub
subroutine i24_save_sub(numnod, mvsiz, nisub, s_addsubm, s_lisubm, s_typsub, nisubmax, i_stok, ie, itypsub, nin, i, nn, nft, addsubm, lisubm, typsub, intarean, intcarea, isensint, fxi, fyi, fzi, fni, dt12, fsavsub1, fsavparit, nrtse, irtse, nsne, is2se, is2pt, nsnr)
Definition i24_save_sub.F:40
i24for3
subroutine i24for3(output, jlt, a, v, ibcc, icodt, fsav, gap, fric, ms, visc, viscf, noint, stfn, itab, cn_loc, stiglo, stifn, stif, fskyi, isky, n1, n2, n3, h1, h2, h3, h4, fcont, pene, ix1, ix2, ix3, ix4, nsvg, ivis2, neltst, ityptst, dt2t, subtria, gapv, inacti, index, niskyfi, kinet, newfront, isecin, nstrf, secfcum, x, irect, ce_loc, mfrot, ifq, frot_p, secnd_fr, alpha0, ibag, icontact, irtlm, viscn, vxi, vyi, vzi, msi, kini, nin, nisub, lisub, addsubs, addsubm, lisubs, lisubm, fsavsub, cand_n, ilagm, icurv, fncont, ftcont, nsn, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi, iadm, rcurvi, rcontact, acontact, pcontact, anglmi, padm, intth, phi, fthe, ftheskyi, temp, tempi, rstif, iform, mskyi_sms, iskyi_sms, nsms, cand_n_n, pene_old, stif_old, mbinflg, ilev, igsti, kmin, intply, iply, inod_pxfem, nm1, nm2, nm3, nrebou, irtse, nsne, is2se, is2pt, msegtyp, jtask, isensint, fsavparit, nft, h3d_data, fricc, viscffric, fric_coefs, t2main_sms, intnitsche, forneqsi, iorthfric, fric_coefs2, fricc2, viscffric2, nforth, nfisot, indexorth, indexisot, dir1, dir2, t2fac_sms, f_pfit, tagncont, kloadpinter, loadpinter, loadp_hyd_inter, typsub, inflg_subs, inflg_subm, ninloadp, dgaploadint, s_loadpinter, dist, ixx, interefric, intcarea, parameters, penref, kmax, s_addsubm, s_lisubm, s_typsub, nisubmax, i_stok, nrtm, nrtse, nsnr)
Definition i24for3.F:82
i24sms2
subroutine i24sms2(jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, nin, noint, mskyi_sms, iskyi_sms, nsms, kt, c, cf, dtmini, dti, irtse, nsne, is2se, is2pt, t2main_sms, t2fac_sms)
Definition i24for3.F:4749
i24ass2
subroutine i24ass2(jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, fxi, fyi, fzi, fskyi, isky, niskyfi, nin, noint, intth, phi, ftheskyi, phi1, phi2, phi3, phi4, itab, intply, iply, inod, irtse, nsne, is2se, is2pt, tagip)
Definition i24for3.F:3791
i24ass0
subroutine i24ass0(jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, fxi, fyi, fzi, a, stifn, nin, intth, phi, fthe, phi1, phi2, phi3, phi4, intply, iply, inod, irtse, nsne, is2se, is2pt, jtask)
Definition i24for3.F:4425
i24for3e
subroutine i24for3e()
Definition i24for3e.F:27
i24for1_fic
subroutine i24for1_fic(npt, irtse, nsne, is2se, is2pt, ns, fxi, fyi, fzi, for1, inega)
Definition i24for3e.F:148
i24forc_fic
subroutine i24forc_fic(npt, irtse, nsne, is2se, is2pt, ns, nfic, fici, fics, ixss)
Definition i24for3e.F:39
i24for1_ficr
subroutine i24for1_ficr(npt, irtse, nsne, is2se, is2pt, ns, fxi, fyi, fzi, for1, inega, ni)
Definition i24for3e.F:285
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
debug_mod
Definition debug_mod.F:47
h3d_mod
Definition h3d_mod.F:311
message_mod
Definition message_mod.F:1257
parameters_mod
Definition parameters_mod.F:49
plyxfem_mod
Definition plyxfem_mod.F:31
plyxfem_mod::ply
type(ply_data), dimension(:), allocatable ply
Definition plyxfem_mod.F:92
plyxfem_mod::plyskyi
type(ply_data), allocatable plyskyi
Definition plyxfem_mod.F:93
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::is2pt_fi
type(int_pointer), dimension(:), allocatable is2pt_fi
Definition tri7box.F:537
tri7box::inflg_subsfi
type(int_pointer), dimension(:), allocatable inflg_subsfi
Definition tri7box.F:505
tri7box::ftheskyfi
type(real_pointer), dimension(:), allocatable ftheskyfi
Definition tri7box.F:449
tri7box::fnconti
type(real_pointer2), dimension(:), allocatable fnconti
Definition tri7box.F:510
tri7box::is2se_fi
type(int_pointer2), dimension(:), allocatable is2se_fi
Definition tri7box.F:536
tri7box::efricgfi
type(real_pointer), dimension(:), allocatable efricgfi
Definition tri7box.F:511
tri7box::stnfi
type(real_pointer), dimension(:), allocatable stnfi
Definition tri7box.F:449
tri7box::lisubsfi
type(int_pointer), dimension(:), allocatable lisubsfi
Definition tri7box.F:501
tri7box::irtse_fi
type(int_pointer2), dimension(:), allocatable irtse_fi
Definition tri7box.F:535
tri7box::afi
type(real_pointer2), dimension(:), allocatable afi
Definition tri7box.F:459
tri7box::fskyfi
type(real_pointer2), dimension(:), allocatable fskyfi
Definition tri7box.F:459
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::iskyfi
type(int_pointer), dimension(:), allocatable iskyfi
Definition tri7box.F:480
tri7box::isedge_fi
type(int_pointer), dimension(:), allocatable isedge_fi
Definition tri7box.F:540
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::t2fac_sms_fi
type(real_pointer), dimension(:), allocatable t2fac_sms_fi
Definition tri7box.F:557
tri7box::nlskyfi
integer, dimension(:), allocatable nlskyfi
Definition tri7box.F:512
tri7box::fthefi
type(real_pointer), dimension(:), allocatable fthefi
Definition tri7box.F:449
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
tri7box::t2main_sms_fi
type(int_pointer2), dimension(:), allocatable t2main_sms_fi
Definition tri7box.F:558
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