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