OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i22for3.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!|| i22for3 ../engine/source/interfaces/int22/i22for3.F
25!||--- called by ------------------------------------------------------
26!|| i22mainf ../engine/source/interfaces/int22/i22mainf.F
27!||--- calls -----------------------------------------------------
28!|| i22aera ../engine/source/interfaces/int22/i22subvol.f
29!|| i22ass0 ../engine/source/interfaces/int22/i22assembly.F
30!|| i22ass2 ../engine/source/interfaces/int22/i22assembly.F
31!|| i7ass3 ../engine/source/interfaces/int07/i7ass3.F
32!||--- uses -----------------------------------------------------
33!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
34!|| element_mod ../common_source/modules/elements/element_mod.F90
35!|| h3d_mod ../engine/share/modules/h3d_mod.F
36!|| i22bufbric_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
37!|| i22tri_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
38!|| output_mod ../common_source/modules/output/output_mod.F90
39!||====================================================================
40 SUBROUTINE i22for3(OUTPUT, JLT ,A ,V ,IBC ,ICODT ,
41 2 FSAV ,GAP ,FRIC ,MS ,VISC ,
42 3 VISCF ,NOINT ,STFN ,ITAB ,CB_LOC ,
43 4 STIGLO ,STIFN ,STIF ,FSKYI ,ISKY ,
44 6 NX1 ,NX2 ,NX3 ,NX4 ,NY1 ,
45 7 NY2 ,NY3 ,NY4 ,NZ1 ,NZ2 ,
46 8 NZ3 ,NZ4 ,LB1 ,LB2 ,LB3 ,
47 9 LB4 ,LC1 ,LC2 ,LC3 ,LC4 ,
48 A P1 ,P2 ,P3 ,P4 ,FCONT ,
49 B IX1 ,IX2 ,IX3 ,IX4 ,NSVG ,
50 C IVIS2 ,NELTST ,ITYPTST ,DT2T ,INTTH ,
51 D GAPV ,INACTI ,CAND_P ,INDEX ,NISKYFI,
52 E KINET ,NEWFRONT ,ISECIN ,NSTRF ,SECFCUM,
53 F X ,IRECT ,CE_LOC ,MFROT ,IFQ ,
54 G FROT_P ,CAND_FX ,CAND_FY ,CAND_FZ ,ALPHA0,
55 H IFPEN ,IBAG ,ICONTACT ,
56 J VISCN ,VXI ,VYI ,VZI ,MSI ,
57 K KINI ,NIN ,NISUB ,LISUB ,ADDSUBS,
58 L ADDSUBM ,LISUBS ,LISUBM ,FSAVSUB ,CAND_N ,
59 M ILAGM ,ICURV ,PRES ,FNCONT ,MS0 ,
60 N N_SCUT ,N_SURF ,CoG ,CAND_E ,Swet ,
61 O ELBUF_TAB ,X1 ,X2 ,X3 ,X4 ,
62 3 Y1 ,Y2 ,Y3 ,Y4 ,Z1 ,
63 4 Z2 ,Z3 ,Z4 ,IXS ,NV46 ,
64 5 Delta ,ISENSINT ,FSAVPARIT ,IPARG ,H3D_DATA)
65C-----------------------------------------------
66C D e s c r i p t i o n
67C-----------------------------------------------
68C Interface Type22 (/INTER/TYPE22) is an FSI coupling method based on cut cell method.
69C This experimental cut cell method is not completed, abandoned, and is not an official option.
70C
71C-----------------------------------------------
72C M o d u l e s
73C-----------------------------------------------
74 USE elbufdef_mod
75 USE i22tri_mod
76 USE i22bufbric_mod
77 USE h3d_mod
78 USE output_mod
79 use element_mod , only : nixs
80C-----------------------------------------------
81C I m p l i c i t T y p e s
82C-----------------------------------------------
83#include "implicit_f.inc"
84#include "comlock.inc"
85C-----------------------------------------------
86C G l o b a l P a r a m e t e r s
87C-----------------------------------------------
88#include "mvsiz_p.inc"
89C-----------------------------------------------
90C C o m m o n B l o c k s
91C-----------------------------------------------
92#include "com01_c.inc"
93#include "com04_c.inc"
94#include "com06_c.inc"
95#include "com08_c.inc"
96#include "scr07_c.inc"
97#include "scr11_c.inc"
98#include "scr14_c.inc"
99#include "scr16_c.inc"
100#include "parit_c.inc"
101#include "param_c.inc"
102C-----------------------------------------------
103C D u m m y A r g u m e n t s
104C-----------------------------------------------
105 TYPE(output_), INTENT(INOUT) :: OUTPUT
106 TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
107
108 INTEGER NELTST,ITYPTST,JLT,IBC,IBCM,IBCS,IVIS2,INACTI,IBAG,NIN,
109 . ICODT(*), ITAB(*), ISKY(*), KINET(*),
110 . MFROT, IFQ, NOINT,NEWFRONT,ISECIN, NSTRF(*),
111 . IRECT(4,*),IFPEN(*) ,ICONTACT(*), CAND_N(*),
112 . KINI(*),
113 . ISET, NISKYFI,INTTH,IFORM,NV46
114 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
115 . CB_LOC(MVSIZ),CE_LOC(MVSIZ),INDEX(MVSIZ),NSVG(MVSIZ),
116 . NISUB, LISUB(*), ADDSUBS(*), ADDSUBM(*), LISUBS(*),
117 . LISUBM(*),ILAGM,ICURV(3),CAND_E(*),ixs(nixs,*),
118 . ISENSINT(*), IPARG(NPARG,*)
119 my_real
120 . STIGLO,CAND_P(*),FROT_P(*), X(3,*),MS0(*),
121 . A(3,*), MS(*), V(3,*), FSAV(*),FCONT(3,*),
122 . CAND_FX(*),CAND_FY(*),CAND_FZ(*),ALPHA0,
123 . GAP, FRIC,VISC,VISCF,VIS,DT2T,STFN(*),STIFN(*),
124 . FSKYI(LSKYI,NFSKYI),FSAVSUB(NTHVKI,*),FNCONT(3,*)
125 my_real
126 . NX1(MVSIZ), NX2(MVSIZ), NX3(MVSIZ), NX4(MVSIZ),
127 . NY1(MVSIZ), NY2(MVSIZ), NY3(MVSIZ), NY4(MVSIZ),
128 . NZ1(MVSIZ), NZ2(MVSIZ), NZ3(MVSIZ), NZ4(MVSIZ),
129 . lb1(mvsiz), lb2(mvsiz), lb3(mvsiz), lb4(mvsiz),
130 . lc1(mvsiz), lc2(mvsiz), lc3(mvsiz), lc4(mvsiz),
131 . p1(mvsiz), p2(mvsiz), p3(mvsiz), p4(mvsiz), stif(mvsiz),
132 . gapv(mvsiz),
133 . secfcum(7,numnod,nsect),tmp(mvsiz),
134 . stifsav(mvsiz), viscn(*),
135 . vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),msi(mvsiz),
136 . pres(*), rstif ,
137 . cog(1:3,nbcut_max,mvsiz), n_scut(3,nbcut_max, mvsiz), swet(nbcut_max,mvsiz)
138 my_real ::
139 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz),
140 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
141 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz),
142 . delta(4,nbcut_max,mvsiz),fsavparit(nisub+1,11,*),face
143 TYPE(h3d_database) :: H3D_DATA
144
145C-----------------------------------------------
146C L o c a l V a r i a b l e s
147C-----------------------------------------------
148 TYPE(G_BUFEL_) ,POINTER :: GBUF1, GBUF2
149 INTEGER I, J1, J2, IG, J, JG , K0,NBINTER,K1S,K,IL,IE, NN, NI,NA1,NA2
150 INTEGER IBv, Jv, NBCUTv, IEv, LLT1, LLT2
151 my_real
152 . FXI(MVSIZ), FYI(MVSIZ), FZI(MVSIZ), FNI(MVSIZ),
153 . FXT(MVSIZ),FYT(MVSIZ),FZT(MVSIZ),
154 . FX1(MVSIZ), FX2(MVSIZ), FX3(MVSIZ), FX4(MVSIZ),
155 . FY1(MVSIZ), FY2(MVSIZ), FY3(MVSIZ), FY4(MVSIZ),
156 . FZ1(MVSIZ), FZ2(MVSIZ), FZ3(MVSIZ), FZ4(MVSIZ),
157 . N1(MVSIZ), N2(MVSIZ), N3(MVSIZ), PENE(MVSIZ),
158 . VIS2(MVSIZ), DTMI(MVSIZ), XMU(MVSIZ),STIF0(MVSIZ),
159 . H1(MVSIZ), H2(MVSIZ), H3(MVSIZ), H4(MVSIZ),
160 . VX(MVSIZ), VY(MVSIZ), VZ(MVSIZ), VN(MVSIZ),DIST(MVSIZ),
161 . VNX, VNY, VNZ, AA, CRIT,S2,RDIST,
162 . v2, fm2, visca, fac,ff(3),alphi,alpha,beta,
163 . fx, fy, fz, f2, mas2, m2sk, dtmi0,dti,ft,fn,fmax,ftn,
164 . facm1, econtt, econvt, h0, la1, la2, la3, la4,
165 . d1,d2,d3,d4,a1,a2,a3,a4,
166 . fsav1, fsav2, fsav3, fsav4, fsav5, fsav6, fsav7, fsav8,
167 . fsav9, fsav10, fsav11, fsav12, fsav13, fsav14, fsav15, ffo,
168 . e10, h0d, s2d,
169 . la1d,la2d,la3d,la4d,t1,t1d,t2,t2d,ffd,visd,facd,d1d,
170 . p1s(mvsiz),p2s(mvsiz),p3s(mvsiz),p4s(mvsiz),
171 . d2d,d3d,d4d,vnxd,vnyd,vnzd,v2d,fm2d,f2d,aad,fxd,fyd,fzd,
172 . a1d,a2d,a3d,a4d,vv,ax1,ax2,ay1,ay2,az1,az2,ax,ay,az,
173 . area,p,vv1,vv2,v21,dmu, dti2,h00 ,a0x,a0y,a0z,rx,ry,rz,
174 . anx,any,anz,aan,aax,aay,aaz ,rr,rs,aaa ,tm,ts
175 my_real
176 . surfx,surfy,surfz,n_surf(3,*),m1,m2,mf,theta
177 my_real
178 . st1(mvsiz),st2(mvsiz),st3(mvsiz),st4(mvsiz),stv(mvsiz),
179 . kt(mvsiz),c(mvsiz),cf(mvsiz),
180 . ks(mvsiz),k1(mvsiz),k2(mvsiz),k3(mvsiz),k4(mvsiz),
181 . cs(mvsiz),c1(mvsiz),c2(mvsiz),c3(mvsiz),c4(mvsiz),
182 . cx,cy,cfi,aux,phi1(mvsiz), phi2(mvsiz), phi3(mvsiz),
183 . phi4(mvsiz), n_scut_(3)
184
185 INTEGER JSUB,KSUB,JJ,KK,IN,NSUB,IBID, NG, IB,IV, NBCUT,ICUT, NP_RECT(MVSIZ),
186 . IEM1,IEM2, IBM1, IBM2, JM1,JM2
187 INTEGER NG1,NG2,IDLOC1,IDLOC2,NP
188
189 my_real
190 . FSAVSUB1(15,NISUB),IMPX,IMPY,IMPZ,PP1,PP2,PP3,PP4,BID
191 my_real :: DP(MVSIZ), W(4,MVSIZ), Q, Slag, TMP2(3), Pt1(3),Pt2(3), Pt3(3), Pt4(3),distance(4),Dsum,UN1,UN2,ZC1,ZC2
192
193 INTEGER :: ICELL1, ICELL2, ICELLv, MCELL1, MCELL2, IB1, IB2
194C-----------------------------------------------
195 INTERFACE
196 FUNCTION I22AERA(NPTS,P, C)
197 INTEGER :: NPTS
198 my_real :: P(3,NPTS), I22AERA(3), C(3)
199 END FUNCTION
200 END INTERFACE
201C-----------------------------------------------
202C S o u r c e L i n e s
203C-----------------------------------------------
204 econtt = zero
205 econvt = zero
206 DO i=1,jlt
207 stif0(i) = stif(i)
208 ENDDO
209
210
211 bid = zero
212 ibid = 0
213
214! !---------------------------------!
215! ! PRESSURE CORRECTION !
216! !---------------------------------!
217! IF (IBAG/=0)THEN
218!#include "lockon.inc"
219! DO I=1,JLT
220! PRES(CE_LOC(I)) = PRES(CE_LOC(I)) + FNI(I)
221! ENDDO
222!#include "lockoff.inc"
223! END IF
224
225C---------------------------------
226#include "lockon.inc"
227 econtv = econtv + econvt
228 econt = econt + econtt
229#include "lockoff.inc"
230C---------------------------------
231C
232 !---------------------------------!
233 ! PARAMETER INIT. !
234 !---------------------------------!
235 DO i=1,jlt
236 IF(ix3(i)/=ix4(i))THEN
237 np_rect(i) = 4
238 w(1:4,i) = fourth
239 ELSE
240 np_rect(i) = 3
241 w(1:4,i) = third
242 ENDIF
243 ENDDO
244
245 !---------------------------------!
246 ! CONTACT FORCES. !
247 !---------------------------------!
248
249 DO i=1,jlt
250
251c print *, " I===== : ", I
252c print *, " CELL : ", IXS(11,BRICK_LIST(NIN,CB_LOC(I))%ID)
253c print *, " NBCUT : ", BRICK_LIST(NIN,CB_LOC(I))%NBCUT
254c print *, " SHELL : ", ITAB(IRECT(1:4,iabs(CE_LOC(I))))
255c print *, " NCYCLE : ", NCYCLE
256
257
258 ib = cb_loc(i)
259 ie = brick_list(nin,ib)%ID
260 nbcut = brick_list(nin,ib)%NBCUT
261
262 fx1(i) = zero
263 fy1(i) = zero
264 fz1(i) = zero
265 fx2(i) = zero
266 fy2(i) = zero
267 fz2(i) = zero
268 fx3(i) = zero
269 fy3(i) = zero
270 fz3(i) = zero
271 fx4(i) = zero
272 fy4(i) = zero
273 fz4(i) = zero
274 fxi(i) = zero
275 fyi(i) = zero
276 fzi(i) = zero
277 dp(i) = zero
278 fni(i) = zero
279
280c IF(NBCUT==0) print *, " C Y C L E"
281 IF(nbcut==0) cycle
282
283
284 pt1(1:3) = (/x1(i),y1(i),z1(i)/)
285 pt2(1:3) = (/x2(i),y2(i),z2(i)/)
286 pt3(1:3) = (/x3(i),y3(i),z3(i)/)
287 pt4(1:3) = (/x4(i),y4(i),z4(i)/)
288 tmp(1:3) = i22aera( np_rect(i), (/pt1,pt2,pt3,pt4/) ,tmp2) !to be optimized, almoste calculated in i22cor or i22dst3
289 slag= sqrt(sum(tmp(1:3)*tmp(1:3)))
290
291 DO icut=1,nbcut
292
293 IF(cand_e(i)<0) cycle !ICUT=1,8 do not project shell which do not intersect
294 !
295 ! +------------+------------+
296 ! + | + | + 20 with facette 1 + ( + : intersection)
297 ! + | + | + 20 with facette 2 - ( - : couple sans intersection)
298 ! + 20 | + 30 | + 30 with facette 1 -
299 ! + | + | + 30 with facette 2 +
300 ! + IX | I + I | IX +
301 ! +------------+------------+
302 !
303 ! +--------\----+
304 ! + \I +
305 ! + ---+-\
306 ! + 20 + | <---non intersecting edges but near virtual face (closure hypothesis)
307 ! + __+_/
308 ! + IX /II+
309 ! +---------/---+
310 !
311 ! +--------\----+-------|----+
312 ! + \I + | +
313 ! + ---+-\ | +
314 ! + 20 + | 30 | +
315 ! + __+_/ | +
316 ! + IX /II+ I | IX +
317 ! +---------/---+-------|----+
318 ! |
319 ! | <------ face of this edge must not be projected on virtual cut surface from cell 20
320 ! this case is not treated (lagrangian face too complex or multiple surface, self-impact, etc...)
321 ! Use a simple surface : smooth and without holes
322
323 !-------------------------------!
324 ! GET RELATIVE PRESSURE. !
325 ! FROM BOTH SIDE OF CUT SURFACE !
326 !-------------------------------!
327 ! %P is pressure in cut cell buffer.
328 !Bijection : Scut=Icut <-> (icell1,icell2)
329 icell1 = icut
330 icell2 = 9
331 !Surjection (icell1,icell2) ->> ( (mcell1,ibv1) , (mcell2,ibv2) )
332 iem1 = brick_list(nin,ib)%POLY(icell1)%WhereIsMain(3)
333 iem2 = brick_list(nin,ib)%POLY(icell2)%WhereIsMain(3)
334 jm1 = brick_list(nin,ib)%POLY(icell1)%WhereIsMain(1)
335 jm2 = brick_list(nin,ib)%POLY(icell2)%WhereIsMain(1)
336 IF(iem1==ie)THEN
337 ibm1 = ib
338 ELSE
339 IF(jm1<=nv46)THEN
340 ibm1 = brick_list(nin,ib)%Adjacent_Brick(jm1,4)
341 ELSE
342 j1 = jm1/10
343 j2 = mod(jm1,10)
344 ibm1 = brick_list(nin,ib )%Adjacent_Brick(j1,4)
345 ibm1 = brick_list(nin,ibm1)%Adjacent_Brick(j2,4)
346 ENDIF
347 ENDIF
348 IF(iem2==ie)THEN
349 ibm2 = ib
350 ELSE
351 IF(jm2<=nv46)THEN
352 ibm2 = brick_list(nin,ib)%Adjacent_Brick(jm2,4)
353 ELSE
354 j1 = jm2/10
355 j2 = mod(jm2,10)
356 ibm2 = brick_list(nin,ib )%Adjacent_Brick(j1,4)
357 ibm2 = brick_list(nin,ibm2)%Adjacent_Brick(j2,4)
358 ENDIF
359 ENDIF
360 ng1 = brick_list(nin,ibm1)%NG
361 ng2 = brick_list(nin,ibm2)%NG
362 idloc1 = brick_list(nin,ibm1)%IDLOC
363 idloc2 = brick_list(nin,ibm2)%IDLOC
364 gbuf1 => elbuf_tab(ng1)%GBUF
365 gbuf2 => elbuf_tab(ng2)%GBUF
366 llt1 = iparg(2,ng1)
367 llt2 = iparg(2,ng2)
368 p1(i) = -third * (gbuf1%SIG(idloc1 )+
369 . gbuf1%SIG(idloc1 +llt1 )+
370 . gbuf1%SIG(idloc1 +llt1*2) )
371 m1 = brick_list(nin,ibm1)%MACH
372 un1 = brick_list(nin,ib )%POLY(icell1)%FACE0%U_N(1) !for poly_id=ICELL1 in [1,NBCUT], there is a single cut section
373 zc1 = brick_list(nin,ibm1)%RHOC
374
375 p2(i) = -third * (gbuf2%SIG(idloc2 )+
376 . gbuf2%SIG(idloc2 +llt2 )+
377 . gbuf2%SIG(idloc2 +llt2*2) )
378 un2 = brick_list(nin,ib )%POLY(icell2)%FACE0%U_N(icell1) !icell2=9 => cut surface sharing icell1=icut and icell2=9, is Scut_id=icut=icell1
379 zc2 = brick_list(nin,ibm2)%RHOC
380 m2 = brick_list(nin,ibm2)%MACH
381
382 mf = half*(m1+m2) ! with grid velocity Mi must be calculated consequently in alefvm_stress*
383 theta = min(one,mf)
384
385 p1(i) = p1(i) + zc1*un1*theta
386 p2(i) = p2(i) + zc2*un2*theta
387
388 dp(i) = p1(i) - p2(i)
389
390 !print *, i, dp(i), dp(i)+ZC1*UN1-ZC2*UN2
391 !-------------------------------!
392 q = sum( n_surf(1:3,iabs(cand_e(i))) * n_scut(1:3,icut,i) ) !N_SCUT(1:3,ICUT,I)= S*n(1:3)
393 ffo = dp(i) * swet(icut,i)
394 IF(q<zero) ffo = -ffo !sign change to adjust normal direction (non oriented surface)
395
396 !force at cog
397 ff(:) = ffo * n_surf(:,iabs(cand_e(i)))
398 !rint *, "F_cog=", FF(:)
399
400 !distributed forces using weighting coefficients
401 fx1(i)= fx1(i) + delta(1,icut,i) * ff(1)
402 fy1(i)= fy1(i) + delta(1,icut,i) * ff(2)
403 fz1(i)= fz1(i) + delta(1,icut,i) * ff(3)
404
405 fx2(i)= fx2(i) + delta(2,icut,i) * ff(1)
406 fy2(i)= fy2(i) + delta(2,icut,i) * ff(2)
407 fz2(i)= fz2(i) + delta(2,icut,i) * ff(3)
408
409 fx3(i)= fx3(i) + delta(3,icut,i) * ff(1)
410 fy3(i)= fy3(i) + delta(3,icut,i) * ff(2)
411 fz3(i)= fz3(i) + delta(3,icut,i) * ff(3)
412
413 IF(np_rect(i)==4)THEN
414 fx4(i)= fx4(i) + delta(4,icut,i) * ff(1)
415 fy4(i)= fy4(i) + delta(4,icut,i) * ff(2)
416 fz4(i)= fz4(i) + delta(4,icut,i) * ff(3)
417 ENDIF
418
419C print *, "FF=", FF(1), ITAB(IRECT(1:4,CE_LOC(I)))
420
421 fxi(i)= fxi(i) + ff(1)
422 fyi(i)= fyi(i) + ff(2)
423 fzi(i)= fzi(i) + ff(3)
424
425 if(ibug22_fcont==-1)then
426 print *, "######################################################"
427 print *, "##I22FOR ; candidat I=", i, " ICUT=", icut
428 print *, "######################################################"
429 print *, " JLT : ", i,"/",jlt
430 print *, " CELL : ", ixs(11,brick_list(nin,cb_loc(i))%ID)
431 print *, " SHELL : ", itab(irect(1:4,iabs(ce_loc(i))))
432 print *, " ICUT : ", icut
433 print *, " NCYCLE : ", ncycle
434 print *, " P1 =", p1(i)
435 print *, " P2 =", p2(i)
436 print *, " DP =", dp(i)
437 print *, " Swet(ICUT,I) =", swet(icut,i)
438 print *, " F=DP*Swet =", dp(i) * swet(icut,i)
439 print *, " Slag =", slag
440 print *, " N_SURF() =", n_surf(:,iabs(cand_e(i)))
441 print *, "-----------------"
442 write (*,fmt='(A,4E20.12,A,E20.12)') " DELTA : ", delta(1:4,icut,i), " | SUM=", sum(delta(1:4,icut,i))
443 write(*,fmt='(A,I8,A,3E30.16)')"FX(",itab(irect(1,iabs(ce_loc(i)))),") = ", fx1(i),fy1(i),fz1(i)
444 write(*,fmt='(A,I8,A,3E30.16)')"FX(",itab(irect(2,iabs(ce_loc(i)))),") = ", fx2(i),fy2(i),fz2(i)
445 write(*,fmt='(A,I8,A,3E30.16)')"FX(",itab(irect(3,iabs(ce_loc(i)))),") = ", fx3(i),fy3(i),fz3(i)
446 write(*,fmt='(A,I8,A,3E30.16)')"FX(",itab(irect(4,iabs(ce_loc(i)))),") = ", fx4(i),fy4(i),fz4(i)
447 print *, "-----------------"
448 print *, " FXI(I)=",fxi(i)
449 print *, " FYI(I)=",fyi(i)
450 print *, " FZI(I)=",fzi(i)
451 print *, ""
452 endif
453 enddo!next ICUT
454
455 !previous treatment was for face0
456 !now treating other faces
457
458 DO j=1,6
459 np = brick_list(nin,ib)%PCUT(8+j)%NP
460 IF(np<=0)cycle
461 icell1 = 9
462 face = brick_list(nin,ib)%POLY(9)%FACE(j)%Surf
463 IF(face<=zero) cycle
464 !-------------------------------!
465 ! GET RELATIVE PRESSURE. !
466 ! FROM BOTH SIDE OF CUT SURFACE !
467 !-------------------------------!
468 !Bijection : Scut=Icut <-> (icell1,icell2)
469 !Surjection (icell1,icell2) ->> ( (mcell1,ibv1) , (mcell2,ibv2) )
470 iem1 = brick_list(nin,ib)%POLY(icell1)%WhereIsMain(3)
471 jm1 = brick_list(nin,ib)%POLY(icell1)%WhereIsMain(1)
472 ibv = brick_list(nin,ib)%Adjacent_Brick(j,4)
473 iev = brick_list(nin,ib)%Adjacent_Brick(j,1)
474 jv = brick_list(nin,ib)%Adjacent_Brick(j,5)
475 IF(iev == 0 .OR. ibv == 0) cycle !closed surface is outside the domain
476 nbcutv = brick_list(nin,ibv)%NBCUT
477 IF(nbcutv == 0)THEN
478 icell2 = 1
479 ELSE
480 icell2 = 9
481 ENDIF
482 IF(icell2 == 0) icell2 = 1 !closed surface hypothesis
483 iem2 = brick_list(nin,ibv)%POLY(icell2)%WhereIsMain(3)
484 jm2 = brick_list(nin,ibv)%POLY(icell2)%WhereIsMain(1)
485 IF(iem1==ie)THEN
486 ibm1 = ib
487 ELSE
488 IF(jm1<=nv46)THEN !1<= direct <=6 ; indirect >= 10
489C if (jm1==0)then
490C print *, "i22for3.F:issue"
491C pause
492C endif
493 ibm1 = brick_list(nin,ib)%Adjacent_Brick(jm1,4)
494 ELSE
495 j1 = jm1/10
496 j2 = mod(jm1,10)
497 ibm1 = brick_list(nin,ib )%Adjacent_Brick(j1,4)
498 ibm1 = brick_list(nin,ibm1)%Adjacent_Brick(j2,4)
499 ENDIF
500 ENDIF
501 IF(iem2 == iev)THEN
502 ibm2 = ibv
503 ELSE
504 IF(jm2<=nv46)THEN
505 ibm2 = brick_list(nin,ibv)%Adjacent_Brick(jm2,4)
506 ELSE
507 j1 = jm2/10
508 j2 = mod(jm2,10)
509 ibm2 = brick_list(nin,ib )%Adjacent_Brick(j1,4)
510 ibm2 = brick_list(nin,ibm2)%Adjacent_Brick(j2,4)
511 ENDIF
512 ENDIF
513 ng1 = brick_list(nin,ibm1)%NG
514 ng2 = brick_list(nin,ibm2)%NG
515 idloc1 = brick_list(nin,ibm1)%IDLOC
516 idloc2 = brick_list(nin,ibm2)%IDLOC
517 m1 = brick_list(nin,ibm1)%MACH
518 m2 = brick_list(nin,ibm2)%MACH
519 gbuf1 => elbuf_tab(ng1)%GBUF
520 gbuf2 => elbuf_tab(ng2)%GBUF
521 llt1 = iparg(2,ng1)
522 llt2 = iparg(2,ng2)
523 p1(i) = -third * (gbuf1%SIG(idloc1 )+
524 . gbuf1%SIG(idloc1 +llt1 )+
525 . gbuf1%SIG(idloc1 +llt1*2) )
526 p2(i) = -third * (gbuf2%SIG(idloc2 )+
527 . gbuf2%SIG(idloc2 +llt2 )+
528 . gbuf2%SIG(idloc2 +llt2*2) )
529
530 un1 = brick_list(nin,ib )%POLY(icell1)%FACE(j)%U_N
531 un2 = brick_list(nin,ibv )%POLY(icell2)%FACE(jv)%U_N
532
533 zc1 = brick_list(nin,ibm1)%RHOC
534 zc2 = brick_list(nin,ibm1)%RHOC
535
536 mf = half*(m1+m2)
537 theta = min(one,mf)
538
539 p1(i) = p1(i) + zc1*un1*theta
540 p2(i) = p2(i) + zc2*un2*theta
541 dp(i) = p1(i) - p2(i)
542 !-------------------------------!
543 n_scut_(1:3) = brick_list(nin,ib)%N(j,1:3)*brick_list(nin,ib)%POLY(9)%FACE(j)%Surf
544 q = sum( n_surf(1:3,iabs(cand_e(i))) * n_scut_(1:3) )
545 ffo = dp(i) * swet(8+j,i)
546 IF(q<zero) ffo = -ffo !sign change to adjust normal direction (non oriented surface)
547
548 !force at cog
549 ff(:) = ffo * n_surf(:,iabs(cand_e(i)))
550 !rint *, "F_cog=", FF(:)
551
552 !distributed forces using weighting coefficients
553 fx1(i) = fx1(i) + delta(1,8+j,i) * ff(1)
554 fy1(i) = fy1(i) + delta(1,8+j,i) * ff(2)
555 fz1(i) = fz1(i) + delta(1,8+j,i) * ff(3)
556
557 fx2(i) = fx2(i) + delta(2,8+j,i) * ff(1)
558 fy2(i) = fy2(i) + delta(2,8+j,i) * ff(2)
559 fz2(i) = fz2(i) + delta(2,8+j,i) * ff(3)
560
561 fx3(i) = fx3(i) + delta(3,8+j,i) * ff(1)
562 fy3(i) = fy3(i) + delta(3,8+j,i) * ff(2)
563 fz3(i) = fz3(i) + delta(3,8+j,i) * ff(3)
564
565 IF(np_rect(i)==4)THEN
566 fx4(i) = fx4(i) + delta(4,8+j,i) * ff(1)
567 fy4(i) = fy4(i) + delta(4,8+j,i) * ff(2)
568 fz4(i) = fz4(i) + delta(4,8+j,i) * ff(3)
569 ENDIF
570
571C print *, "FF=", FF(1), ITAB(IRECT(1:4,CE_LOC(I)))
572
573 fxi(i) = fxi(i) + ff(1)
574 fyi(i) = fyi(i) + ff(2)
575 fzi(i) = fzi(i) + ff(3)
576
577 if(ibug22_fcont==-1)then
578 print *, "#################################"
579 print *, "##I22FOR ; facette I=", i
580 print *, "#################################"
581 print *, " CELL : ", ixs(11,brick_list(nin,cb_loc(i))%ID)
582 print *, " SHELL : ", itab(irect(1:4,iabs(ce_loc(i))))
583 print *, " NCYCLE : ", ncycle
584 print *, " P1 =", p1(i)
585 print *, " P2 =", p2(i)
586 print *, " DP =", dp(i)
587 print *, " Swet(8+J,I) =", swet(8+j,i)
588 print *, " F=DP*Swet =", dp(i) * swet(8+j,i)
589 print *, " Slag =", slag
590 print *, " N_SURF() =", n_surf(:,iabs(cand_e(i)))
591 print *, " FN =", ffo
592 print *, "-----------------"
593 write (*,fmt='(A,4E20.12,A,E20.12)') " DELTA : ", delta(1:4,8+j,i), " | SUM=", sum(delta(1:4,8+j,i))
594 write(*,fmt='(A,I8,A,3E30.16)')"FX(",itab(irect(1,iabs(ce_loc(i)))),") = ", fx1(i),fy1(i),fz1(i)
595 write(*,fmt='(A,I8,A,3E30.16)')"FX(",itab(irect(2,iabs(ce_loc(i)))),") = ", fx2(i),fy2(i),fz2(i)
596 write(*,fmt='(A,I8,A,3E30.16)')"FX(",itab(irect(3,iabs(ce_loc(i)))),") = ", fx3(i),fy3(i),fz3(i)
597 write(*,fmt='(A,I8,A,3E30.16)')"FX(",itab(irect(4,iabs(ce_loc(i)))),") = ", fx4(i),fy4(i),fz4(i)
598 print *, "-----------------"
599 print *, " FXI(I)=",fxi(i)
600 print *, " FYI(I)=",fyi(i)
601 print *, " FZI(I)=",fzi(i)
602 print *, ""
603 endif
604 enddo!next J
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626 enddo!next I
627
628C---------------------------------
629C SAUVEGARDE DE L'IMPULSION NORMALE
630C---------------------------------
631 fsav1 = zero
632 fsav2 = zero
633 fsav3 = zero
634 fsav8 = zero
635 fsav9 = zero
636 fsav10 = zero
637 fsav11 = zero
638 DO i=1,jlt
639 impx = fxi(i)*dt12
640 impy = fyi(i)*dt12
641 impz = fzi(i)*dt12
642 fni(i) = sqrt(fxi(i)*fxi(i) + fyi(i)*fyi(i) + fzi(i)*fzi(i))
643 !FNX,FNY,FNZ
644 fsav1 = fsav1 +impx
645 fsav2 = fsav2 +impy
646 fsav3 = fsav3 +impz
647 !|FNX|,|FNY|,|FNZ|
648 fsav8 = fsav8 +abs(impx)
649 fsav9 = fsav9 +abs(impy)
650 fsav10 = fsav10+abs(impz)
651
652 fsav11 = fsav11+fni(i)*dt12
653 ENDDO
654#include "lockon.inc"
655 fsav(1) = fsav(1) + fsav1
656 fsav(2) = fsav(2) + fsav2
657 fsav(3) = fsav(3) + fsav3
658 fsav(8) = fsav(8) + fsav8
659 fsav(9) = fsav(9) + fsav9
660 fsav(10)= fsav(10)+ fsav10
661 fsav(11)= fsav(11)+ fsav11
662#include "lockoff.inc"
663C
664 IF(isensint(1)/=0) THEN
665 DO i=1,jlt
666 fsavparit(1,1,i) = fxi(i)
667 fsavparit(1,2,i) = fyi(i)
668 fsavparit(1,3,i) = fzi(i)
669 ENDDO
670 ENDIF
671
672 !---------------------------------!
673 ! ASSEMBLY PARITH ON/OFF !
674 !---------------------------------!
675 !these array could be removed from inter22 later during source cleaning
676 h1 = zero
677 h2 = zero
678 h3 = zero
679 h4 = zero
680
681 SELECT CASE (iparit)
682 CASE (0)
683c print *, "PARITH=OFF (IPARIT=0)"
684 CALL i22ass0(
685 1 jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
686 2 nsvg ,h1 ,h2 ,h3 ,h4 ,stif ,
687 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,fz2 ,
688 4 fx3 ,fy3 ,fz3 ,fx4 ,fy4 ,fz4 ,
689 5 fxi ,fyi ,fzi ,a ,stifn,nin ,
690 6 intth ,bid ,bid ,bid ,bid ,bid ,
691 7 bid )
692 CASE (3)
693c print *, "PARITH=ON (new) (IPARIT=3)"
694 CALL i7ass3(
695 1 jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
696 2 nsvg ,h1 ,h2 ,h3 ,h4 ,stif ,
697 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,fz2 ,
698 4 fx3 ,fy3 ,fz3 ,fx4 ,fy4 ,fz4 ,
699 5 fxi ,fyi ,fzi ,a ,stifn)
700 CASE DEFAULT
701 CALL i22ass2(
702 1 jlt ,ix1 ,ix2 ,ix3 ,ix4 ,itab ,
703 2 nsvg ,h1 ,h2 ,h3 ,h4 ,stif ,
704 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,fz2 ,
705 4 fx3 ,fy3 ,fz3 ,fx4 ,fy4 ,fz4 ,
706 5 fxi ,fyi ,fzi ,fskyi ,isky ,niskyfi,
707 6 nin ,noint ,intth ,bid ,bid ,bid ,
708 7 bid ,bid ,bid ,cb_loc ,ce_loc ,irect ,
709 8 ixs)
710
711 END SELECT
712
713C
714 !---------------------------------!
715 ! /ANIM/VECT/CONT !
716 !---------------------------------!
717 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT>0)THEN
718#include "lockon.inc"
719 DO i=1,jlt
720 ib = cb_loc(i)
721 nbcut = brick_list(nin,ib)%NBCUT
722 IF(nbcut == 0) cycle
723 fcont(1,ix1(i)) =fcont(1,ix1(i)) + fx1(i)
724 fcont(2,ix1(i)) =fcont(2,ix1(i)) + fy1(i)
725 fcont(3,ix1(i)) =fcont(3,ix1(i)) + fz1(i)
726 fcont(1,ix2(i)) =fcont(1,ix2(i)) + fx2(i)
727 fcont(2,ix2(i)) =fcont(2,ix2(i)) + fy2(i)
728 fcont(3,ix2(i)) =fcont(3,ix2(i)) + fz2(i)
729 fcont(1,ix3(i)) =fcont(1,ix3(i)) + fx3(i)
730 fcont(2,ix3(i)) =fcont(2,ix3(i)) + fy3(i)
731 fcont(3,ix3(i)) =fcont(3,ix3(i)) + fz3(i)
732 fcont(1,ix4(i)) =fcont(1,ix4(i)) + fx4(i)
733 fcont(2,ix4(i)) =fcont(2,ix4(i)) + fy4(i)
734 fcont(3,ix4(i)) =fcont(3,ix4(i)) + fz4(i)
735 if(ibug22_fcont==-1)then
736 print *, "#################################"
737 print *, "##I22FOR ; facette I=", i
738 print *, "## FCONT /ANIM/VECT/CONT #"
739 print *, "#################################"
740 print *, " CELL : ", ixs(11,brick_list(nin,cb_loc(i))%ID)
741 print *, " SHELL : ", itab(irect(1:4,iabs(ce_loc(i))))
742 print *, " NCYCLE : ", ncycle
743 print *, "------------------"
744 write(*,fmt='(A,I8,A,3E30.16)')"FCONT(",itab(irect(1,iabs(ce_loc(i)))),") = ", FCONT(1:3,IX1(I))
745 write(*,FMT='(A,I8,A,3E30.16)')"fcont(",itab(IRECT(2,iabs(CE_LOC(I)))),") = ", FCONT(1:3,IX2(I))
746 write(*,FMT='(A,I8,A,3E30.16)')"fcont(",itab(IRECT(3,iabs(CE_LOC(I)))),") = ", FCONT(1:3,IX3(I))
747 write(*,FMT='(A,I8,A,3E30.16)')"fcont(",itab(IRECT(4,iabs(CE_LOC(I)))),") = ", FCONT(1:3,IX4(I))
748 print *, ""
749 endif
750 ENDDO
751#include "lockoff.inc"
752 ENDIF
753
754
755 ! A TRAITER PLUS TARD
756
757 !---------------------------------!
758 ! /ANIM/VECT/PCONT !
759 !---------------------------------!
760.AND. IF((ANIM_V(12)+OUTP_V(12)+H3D_DATA%N_VECT_PCONT>0
761.AND..OR..OR..AND..OR. . ((TT>=OUTPUT%TANIMTT<=OUTPUT%TANIM_STOP)TT>=TOUTP(TT>=H3D_DATA%TH3DTT<=H3D_DATA%TH3D_STOP)
762.AND..OR. . (MANIM>=4MANIM<=15)H3D_DATA%MH3D/=0))
763.OR. . H3D_DATA%N_VECT_PCONT_MAX>0)THEN
764#include "lockon.inc"
765 DO I=1,JLT
766 FNCONT(1,IX1(I)) =FNCONT(1,IX1(I)) + FXI(I)
767 FNCONT(2,IX1(I)) =FNCONT(2,IX1(I)) + FYI(I)
768 FNCONT(3,IX1(I)) =FNCONT(3,IX1(I)) + FZI(I)
769 FNCONT(1,IX2(I)) =FNCONT(1,IX2(I)) + FXI(I)
770 FNCONT(2,IX2(I)) =FNCONT(2,IX2(I)) + FYI(I)
771 FNCONT(3,IX2(I)) =FNCONT(3,IX2(I)) + FZI(I)
772 FNCONT(1,IX3(I)) =FNCONT(1,IX3(I)) + FXI(I)
773 FNCONT(2,IX3(I)) =FNCONT(2,IX3(I)) + FYI(I)
774 FNCONT(3,IX3(I)) =FNCONT(3,IX3(I)) + FZI(I)
775 FNCONT(1,IX4(I)) =FNCONT(1,IX4(I)) + FXI(I)
776 FNCONT(2,IX4(I)) =FNCONT(2,IX4(I)) + FYI(I)
777 FNCONT(3,IX4(I)) =FNCONT(3,IX4(I)) + FZI(I)
778 if(ibug22_fcont==-1)then
779 print *, "#################################"
780 print *, "##I22FOR ; facette I=", i
781 print *, "## FCONT /ANIM/VECT/CONT #"
782 print *, "#################################"
783 print *, " CELL : ", ixs(11,brick_list(nin,cb_loc(i))%ID)
784 print *, " SHELL : ", itab(irect(1:4,iabs(ce_loc(i))))
785 print *, " NCYCLE : ", ncycle
786 print *, "-------------------"
787 write(*,fmt='(A,I8,A,3E30.16)')"FNCONT(",itab(irect(1,iabs(ce_loc(i)))),") = ", fncont(1:3,ix1(i))
788 write(*,fmt='(A,I8,A,3E30.16)')"FNCONT(",itab(irect(2,iabs(ce_loc(i)))),") = ", fncont(1:3,ix2(i))
789 write(*,fmt='(A,I8,A,3E30.16)')"FNCONT(",itab(irect(3,iabs(ce_loc(i)))),") = ", fncont(1:3,ix3(i))
790 write(*,fmt='(A,I8,A,3E30.16)')"FNCONT(",itab(irect(4,iabs(ce_loc(i)))),") = ", fncont(1:3,ix4(i))
791
792 print *, ""
793 endif
794 ENDDO
795#include "lockoff.inc"
796 ENDIF
797C-----------------------------------------------------
798
799
800
801C-----------------------------------------------------
802C
803 if(ibug22_fcont==-1)then
804 print *, "================================================="
805 print *, ""
806 print *, ""
807 print *, ""
808 endif
809
810 RETURN
811 END
alpha
#define alpha
Definition eval.h:35
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
i22ass2
subroutine i22ass2(jlt, ix1, ix2, ix3, ix4, itab, 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, cb_loc, ce_loc, irect, ixs)
Definition i22assembly.F:129
i22ass0
subroutine i22ass0(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)
Definition i22assembly.F:37
i22for3
subroutine i22for3(output, jlt, a, v, ibc, icodt, fsav, gap, fric, ms, visc, viscf, noint, stfn, itab, cb_loc, stiglo, stifn, stif, fskyi, isky, nx1, nx2, nx3, nx4, ny1, ny2, ny3, ny4, nz1, nz2, nz3, nz4, lb1, lb2, lb3, lb4, lc1, lc2, lc3, lc4, p1, p2, p3, p4, fcont, ix1, ix2, ix3, ix4, nsvg, ivis2, neltst, ityptst, dt2t, intth, gapv, inacti, cand_p, index, niskyfi, kinet, newfront, isecin, nstrf, secfcum, x, irect, ce_loc, mfrot, ifq, frot_p, cand_fx, cand_fy, cand_fz, alpha0, ifpen, ibag, icontact, viscn, vxi, vyi, vzi, msi, kini, nin, nisub, lisub, addsubs, addsubm, lisubs, lisubm, fsavsub, cand_n, ilagm, icurv, pres, fncont, ms0, n_scut, n_surf, cog, cand_e, swet, elbuf_tab, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, ixs, nv46, delta, isensint, fsavparit, iparg, h3d_data)
Definition i22for3.F:65
i22subvol
subroutine i22subvol(x, nin, itask, ipari, itab, ixs, ixtg, v, iparg, elbuf_tab, w, igrsh3n)
Definition i22subvol.F:41
i22aera
function i22aera(npts, p, c)
Definition i22subvol.F:2382
i7ass3
subroutine i7ass3(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, i8a, i8stifn)
Definition i7ass3.F:332
min
#define min(a, b)
Definition macros.h:20
h3d_mod
Definition h3d_mod.F:311
i22bufbric_mod
Definition cut-cell-search_mod.F:95
i22tri_mod
Definition cut-cell-search_mod.F:212
i22tri_mod::nbcut_max
integer nbcut_max
Definition cut-cell-search_mod.F:517
h3d_mod::h3d_database
Definition h3d_mod.F:438