OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i22for3.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "com06_c.inc"
#include "com08_c.inc"
#include "scr07_c.inc"
#include "scr11_c.inc"
#include "scr14_c.inc"
#include "scr16_c.inc"
#include "parit_c.inc"
#include "param_c.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ i22for3()

subroutine i22for3 ( integer jlt,
a,
v,
integer ibc,
integer, dimension(*) icodt,
fsav,
gap,
fric,
ms,
visc,
viscf,
integer noint,
stfn,
integer, dimension(*) itab,
integer, dimension(mvsiz) cb_loc,
stiglo,
stifn,
stif,
fskyi,
integer, dimension(*) 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,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) nsvg,
integer ivis2,
integer neltst,
integer ityptst,
dt2t,
integer intth,
gapv,
integer inacti,
cand_p,
integer, dimension(mvsiz) index,
integer niskyfi,
integer, dimension(*) kinet,
integer newfront,
integer isecin,
integer, dimension(*) nstrf,
secfcum,
x,
integer, dimension(4,*) irect,
integer, dimension(mvsiz) ce_loc,
integer mfrot,
integer ifq,
frot_p,
cand_fx,
cand_fy,
cand_fz,
alpha0,
integer, dimension(*) ifpen,
integer ibag,
integer, dimension(*) icontact,
viscn,
vxi,
vyi,
vzi,
msi,
integer, dimension(*) kini,
integer nin,
integer nisub,
integer, dimension(*) lisub,
integer, dimension(*) addsubs,
integer, dimension(*) addsubm,
integer, dimension(*) lisubs,
integer, dimension(*) lisubm,
fsavsub,
integer, dimension(*) cand_n,
integer ilagm,
integer, dimension(3) icurv,
pres,
fncont,
ms0,
n_scut,
n_surf,
cog,
integer, dimension(*) cand_e,
swet,
type (elbuf_struct_), dimension(ngroup), target elbuf_tab,
x1,
x2,
x3,
x4,
y1,
y2,
y3,
y4,
z1,
z2,
z3,
z4,
integer, dimension(nixs,*) ixs,
integer nv46,
delta,
integer, dimension(*) isensint,
fsavparit,
integer, dimension(nparg,*) iparg,
type(h3d_database) h3d_data )

Definition at line 39 of file i22for3.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
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:128
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
i22aera
function i22aera(npts, p, c)
Definition i22subvol.F:2380
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
anim_mod
Definition anim_mod.F:96
h3d_mod
Definition h3d_mod.F:308
i22bufbric_mod
Definition cut-cell-search_mod.F:95
i22bufbric_mod::brick_list
type(brick_entity), dimension(:,:), allocatable, target brick_list
Definition cut-cell-search_mod.F:105
i22tri_mod
Definition cut-cell-search_mod.F:212
i22tri_mod::nbcut_max
integer nbcut_max
Definition cut-cell-search_mod.F:517
i22tri_mod::ibug22_fcont
integer ibug22_fcont
Definition cut-cell-search_mod.F:336