OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fvbag1.F File Reference
#include "implicit_f.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "com06_c.inc"
#include "com08_c.inc"
#include "scr02_c.inc"
#include "scr07_c.inc"
#include "scr18_c.inc"
#include "param_c.inc"
#include "units_c.inc"
#include "task_c.inc"
#include "warn_c.inc"
#include "mvsiz_p.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine fvbag1 (nn, nel, ibuf, elem, njet, ibagjet, rbagjet, nvent, ibaghol, rbaghol, p, rho, tk, u, sspk, x, v, a, nsensor, sensor_tab, fsav, npc, tf, ivolu, rvolu, mpolh, qpolh, epolh, centroid_polh, ppolh, rpolh, gpolh, ssppolh, npolh, ifvnod, rfvnod, ifvtri, ifvpoly, ifvtadr, ifvpolh, ifvpadr, info, nns, nntr, ifv, npolha, dlh, cpapolh, cpbpolh, cpcpolh, rmwpolh, itagel, elsini, icontact, idpolh, elfmass, elfvel, ibufa, elema, tagela, pa, rhoa, tka, ua, brna, nna, ntga, ibpolh, dtpolh, nnt, nelt, xxxa, vvva, ncona, porosity, lgauge, gauge, ityp, igeo, sspka, geo, pm, ipm, tpolh, elfehpy, cpdpolh, cpepolh, cpfpolh, eltg, iparg, mattg, igrouptg, igroupc, elbuf_tab, fext, cfl_coef, pdisp_old, pdisp, h3d_data, itab, wfext, python)

Function/Subroutine Documentation

◆ fvbag1()

subroutine fvbag1 ( integer nn,
integer nel,
integer, dimension(*) ibuf,
integer, dimension(3,*) elem,
integer njet,
integer, dimension(nibjet,*) ibagjet,
rbagjet,
integer nvent,
integer, dimension(nibhol,*) ibaghol,
rbaghol,
p,
rho,
tk,
u,
sspk,
x,
v,
a,
integer, intent(in) nsensor,
type (sensor_str_), dimension(nsensor), intent(in) sensor_tab,
fsav,
integer, dimension(*) npc,
tf,
integer, dimension(*) ivolu,
rvolu,
mpolh,
qpolh,
epolh,
dimension(3,npolh), intent(inout) centroid_polh,
ppolh,
rpolh,
gpolh,
dimension(npolh), intent(inout) ssppolh,
integer, intent(in) npolh,
integer, dimension(3,*) ifvnod,
rfvnod,
integer, dimension(6,*) ifvtri,
integer, dimension(*) ifvpoly,
integer, dimension(*) ifvtadr,
integer, dimension(*) ifvpolh,
integer, dimension(*) ifvpadr,
integer info,
integer nns,
integer nntr,
integer ifv,
integer npolha,
dlh,
cpapolh,
cpbpolh,
cpcpolh,
rmwpolh,
integer, dimension(*) itagel,
elsini,
integer, dimension(*) icontact,
integer, dimension(npolh) idpolh,
elfmass,
elfvel,
integer, dimension(*) ibufa,
integer, dimension(3,*) elema,
integer, dimension(*) tagela,
pa,
rhoa,
tka,
ua,
integer, dimension(8,*) brna,
integer nna,
integer ntga,
integer, dimension(npolh) ibpolh,
dtpolh,
integer nnt,
integer nelt,
xxxa,
vvva,
integer, dimension(16,*) ncona,
porosity,
integer, dimension(3,nbgauge) lgauge,
gauge,
integer ityp,
integer, dimension(npropgi,numgeo) igeo,
sspka,
geo,
pm,
integer, dimension(npropmi,nummat) ipm,
tpolh,
elfehpy,
cpdpolh,
cpepolh,
cpfpolh,
integer, dimension(*) eltg,
integer, dimension(nparg,ngroup) iparg,
integer, dimension(*) mattg,
integer, dimension(numeltg) igrouptg,
integer, dimension(numelc) igroupc,
type(elbuf_struct_), dimension(ngroup) elbuf_tab,
fext,
cfl_coef,
pdisp_old,
pdisp,
type(h3d_database) h3d_data,
integer, dimension(numnod), intent(in) itab,
double precision, intent(inout) wfext,
type(python_) python )

Definition at line 51 of file fvbag1.F.

74C-----------------------------------------------
75C M o d u l e s
76C-----------------------------------------------
77 USE my_alloc_mod
78 USE fvbag_mod
79 USE message_mod
80 USE elbufdef_mod
81 USE fvmbag_meshcontrol_mod
82 USE sensor_mod
83 USE anim_mod
84 USE h3d_mod , only:h3d_database
85 USE cast_mod, only: double_to_my_real
86 USE python_funct_mod, only : python_
87 USE finter_mixed_mod, only: finter_mixed
88C-----------------------------------------------
89C I m p l i c i t T y p e s
90C-----------------------------------------------
91#include "implicit_f.inc"
92C-----------------------------------------------
93C C o m m o n B l o c k s
94C-----------------------------------------------
95#include "com01_c.inc"
96#include "com04_c.inc"
97#include "com06_c.inc"
98#include "com08_c.inc"
99#include "scr02_c.inc"
100#include "scr07_c.inc"
101#include "scr18_c.inc"
102#include "param_c.inc"
103#include "units_c.inc"
104#include "task_c.inc"
105#include "warn_c.inc"
106#include "mvsiz_p.inc"
107C-----------------------------------------------
108C D u m m y A r g u m e n t s
109C-----------------------------------------------
110 INTEGER ,INTENT(IN) :: NSENSOR, NPOLH
111 INTEGER NN, NEL, IBUF(*), ELEM(3,*), NJET, IBAGJET(NIBJET,*),
112 . NVENT, IBAGHOL(NIBHOL,*),
113 . NPC(*), IVOLU(*), IFVNOD(3,*), IFVTRI(6,*),
114 . IFVPOLY(*),IFVTADR(*), IFVPOLH(*), IFVPADR(*), INFO,
115 . NNS, NNTR, IFV, NPOLHA, ITAGEL(*), ICONTACT(*),
116 . IDPOLH(NPOLH), IBUFA(*), ELEMA(3,*), TAGELA(*), BRNA(8,*),
117 . NNA, NTGA, IBPOLH(NPOLH), NNT, NELT, NCONA(16,*),
118 . LGAUGE(3,NBGAUGE), ITYP, IGEO(NPROPGI,NUMGEO)
119 INTEGER ELTG(*)
120 INTEGER IPM(NPROPMI,NUMMAT), IPARG(NPARG,NGROUP)
121 INTEGER MATTG(*), IGROUPTG(NUMELTG), IGROUPC(NUMELC)
122
123 my_real rbagjet(nrbjet,*), rbaghol(nrbhol,*), p(*), rho(*),
124 . tk(*), u(3,*), sspk(*), x(3,numnod), v(3,numnod), a(3,numnod),
125 . fsav(nthvki), tf(*), rvolu(*), mpolh(npolh), qpolh(3,npolh),
126 . epolh(npolh), ppolh(npolh), rpolh(npolh), gpolh(npolh), rfvnod(2,*), dlh,
127 . cpapolh(npolh), cpbpolh(npolh), cpcpolh(npolh), rmwpolh(npolh), elsini(*),
128 . elfmass(*), elfvel(*), pa(*), rhoa(*), tka(*),sspka(*), ua(3,*),
129 . dtpolh(npolh), xxxa(3,*), vvva(3,*), porosity(*),
130 . gauge(llgauge,nbgauge), geo(npropg,numgeo), pm(npropm,nummat), elfehpy(*),
131 . tpolh(npolh), cpdpolh(npolh), cpepolh(npolh), cpfpolh(npolh), fext(3,numnod), cfl_coef,
132 . pdisp_old, pdisp
133 my_real, INTENT(INOUT) :: ssppolh(npolh)
134
135 TYPE(ELBUF_STRUCT_), DIMENSION(NGROUP) :: ELBUF_TAB
136 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) ,INTENT(IN) :: SENSOR_TAB
137 TYPE(H3D_DATABASE) :: H3D_DATA
138 INTEGER,INTENT(IN)::ITAB(NUMNOD)
139 my_real, INTENT(INOUT) :: centroid_polh(3,npolh)
140 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
141 type(python_) :: PYTHON
142C-----------------------------------------------
143C L o c a l V a r i a b l e s
144C-----------------------------------------------
145 INTEGER I, II, IINJ, ISU, IEL, N1, N2, N3, JEL,
146 . NN1, NN2, NN3, J, JJ, K, KK, IMASS, IFLU, ISENS, IVEL,
147 . ITEMP, NV, ILVOUT, NPOLYG,NTRIA, IPRI,
148 . IVENT, IDEF, IPORT, IPORP, IPORA, IPORT1, IPORP1, IPORA1,
149 . ITVENT, PHDT, I1, I2, NNSA, IVDP,ITTF, IDPDEF,
150 . IDTPDEF,PMAIN,N21,N22,ITYPL,
151 . PREPARE_ANIM,IVENTYP,IMETHOD_LC, IEL1, IEL2, IDT_FVMBAG
152 INTEGER TITREVENT(20)
153 INTEGER, DIMENSION(:,:), ALLOCATABLE :: SURF_INT,SURF_EXT
154 INTEGER, DIMENSION(:,:), ALLOCATABLE :: TAGSURF
155 INTEGER, DIMENSION(:), ALLOCATABLE :: ICONT
156 INTEGER, DIMENSION(:), ALLOCATABLE :: ITAGP
157
158 my_real, DIMENSION(:),ALLOCATABLE :: dm
159 my_real, DIMENSION(:),ALLOCATABLE :: nodarea
160 my_real, DIMENSION(:),ALLOCATABLE :: area_max
161 my_real, DIMENSION(:),ALLOCATABLE :: parea
162 my_real, DIMENSION(:),ALLOCATABLE :: pvolu
163 my_real, DIMENSION(:),ALLOCATABLE :: rhoel
164 my_real, DIMENSION(:),ALLOCATABLE :: tkel
165 my_real, DIMENSION(:),ALLOCATABLE :: sspel
166 my_real, DIMENSION(:),ALLOCATABLE :: elarea
167 my_real, DIMENSION(:),ALLOCATABLE :: dmini
168 my_real, DIMENSION(:),ALLOCATABLE :: dmcpa
169 my_real, DIMENSION(:),ALLOCATABLE :: dmcpb
170 my_real, DIMENSION(:),ALLOCATABLE :: dmcpc
171 my_real, DIMENSION(:),ALLOCATABLE :: dmrmw
172 my_real, DIMENSION(:),ALLOCATABLE :: elsout
173 my_real, DIMENSION(:),ALLOCATABLE :: tmp
174 my_real, DIMENSION(:),ALLOCATABLE :: poro
175 my_real, DIMENSION(:),ALLOCATABLE :: vola
176 my_real, DIMENSION(:),ALLOCATABLE :: hmin
177 my_real, DIMENSION(:),ALLOCATABLE :: vel
178 my_real, DIMENSION(:),ALLOCATABLE :: svent
179 my_real, DIMENSION(:),ALLOCATABLE :: dmcpd
180 my_real, DIMENSION(:),ALLOCATABLE :: dmcpe
181 my_real, DIMENSION(:),ALLOCATABLE :: dmcpf
182 my_real, DIMENSION(:),ALLOCATABLE :: pel
183 my_real, DIMENSION(:),ALLOCATABLE :: dls
184 my_real, DIMENSION(:),ALLOCATABLE :: eldmass
185 my_real, DIMENSION(:),ALLOCATABLE :: eldehpy
186 my_real, DIMENSION(:,:),ALLOCATABLE :: pnod
187 my_real, DIMENSION(:,:),ALLOCATABLE :: norm
188 my_real, DIMENSION(:,:),ALLOCATABLE :: pnorm
189 my_real, DIMENSION(:,:),ALLOCATABLE :: pu
190 my_real, DIMENSION(:,:),ALLOCATABLE :: pw
191 my_real, DIMENSION(:,:),ALLOCATABLE :: vnod
192 my_real, DIMENSION(:,:),ALLOCATABLE :: vgrid
193 my_real, DIMENSION(:,:),ALLOCATABLE :: xxx
194 my_real, DIMENSION(:,:),ALLOCATABLE :: vvv
195
196 my_real x1, y1, z1, x2, y2, z2, x3, y3, z3,
197 . x12, y12, z12, x13, y13, z13, nrx, nry, nrz, area2,
198 . ksi, eta,
199 . area,denom,denom_max, nx, ny,
200 . nz, gamai, cpai, cpbi, cpci, rmwi, ti,
201 . cpdi, cpei, cpfi,
202 . rhoi, datainj(6,njet), scalt, fmass, gmass_old,
203 . gmtot_old, fvel, tstart, tsg, gmass, dgmass, injvel,
204 . gmtot, rmwg, cpa, cpb, cpc, ftemp, efac, cpg,
205 . cvg, wsurf(3), gama, fel(nelt),
206 . dq(3,npolh), de(npolh), dmi, dqi(3), dei, rho1, re1, ux1,
207 . uy1, uz1, rho2, re2, ux2, uy2, uz2, vfx, vfy, vfz, ss,
208 . alpha, rhom, rem, ruxm, ruym, ruzm, p1, p2,
209 . uel(3,nelt), area1, area3, ff, temp, temp2, pext, volg, volph,
210 . areap, areael, xx, yy, zz, qa, qb, dtx, al, dd, ad, ssp,
211 . qx, vv, msini, rmw, cv, gama1, gama2,
212 . qvisc(npolh), amtot, pmean, pmean2,
213 . area_vent(nvent), pm_vent(nvent), scalp, scals, exten,
214 . avent, bvent, aini, fport, fporp, fpora, fport1, fporp1,
215 . fpora1, aout, aout1, deri, pdef, pvoltmp,
216 . dtpdefi, dtpdefc, tvent,aoutot, sventot,
217 . x23, y23, z23, x31, y31, z31, l12, l23, l31, h1, h2, h3,
218 . fac1, fac2, fac3, hm1, area_old, pp, vmax, uu,
219 . pcrit, rhoinj, pinj, vx1, vy1, vz1, vx2, vy2, vz2,
220 . vx3, vy3, vz3, vvx, vvy, vvz, fac,
221 . fvdp, rbid, ttf,
222 . cpd, cpe, cpf, temp0, tt1,pstag, volno,tstope,
223 . enint, massflow, wfext0, dmout, dhout,
224 . cpam, cpbm, cpcm, cpdm, cpem, cpfm, rmwm, tmp_vec(3), vmat(3),ssp1,ssp2,sspt,ntria_tot,
225 . pres_av, rho_av, temp_av, cp_av, cv_av, rgas_av, centroid(3)
226 DOUBLE PRECISION :: CP
227
228 LOGICAL :: boolTAGEL,EXIT_NEG_VOL,INJECTION_STARTED
229 LOGICAL :: UP_SWITCH
230C
231 my_real get_u_func
232 EXTERNAL get_u_func
233 CHARACTER*20 VENTTITLE
234 my_real :: dtold, param_lambda
235 INTEGER :: PARAM_L_TYPE, NNODES,nodeID
236
237C-----------------------------------------------
238C C o m m e n t s
239C-----------------------------------------------
240C NEL : number of triangles (external surface only)
241C NELT : number of triangles (both internal & external surfaces)
242C NNS : number of nodes of TETRA mesh
243C NNA : number of 'additional' nodes
244C NNT : total number of nodes
245C NNTR : number of triangles (all polyhedra)
246C NPOLH : number of polyhedron
247C NPOLY : number of polygons
248
249C-----------------------------------------------
250C S o u r c e L i n e s
251C-----------------------------------------------
252 CALL my_alloc(tagsurf,2,nelt+1)
253 CALL my_alloc(icont,nnt)
254 CALL my_alloc(itagp,npolh)
255 !
256 CALL my_alloc(area_max,npolh)
257 CALL my_alloc(dm,npolh)
258 CALL my_alloc(nodarea,nnt)
259 CALL my_alloc(pnod,3,nns)
260 CALL my_alloc(elarea,nelt)
261 CALL my_alloc(norm,3,nelt)
262 CALL my_alloc(parea,nntr)
263 CALL my_alloc(pnorm,3,nntr)
264 CALL my_alloc(pvolu,npolh)
265 CALL my_alloc(rhoel,nelt)
266 CALL my_alloc(tkel,nelt)
267 CALL my_alloc(sspel,nelt)
268 CALL my_alloc(pu,3,npolh)
269 CALL my_alloc(pw,3,npolh)
270 CALL my_alloc(dmini,npolh)
271 CALL my_alloc(dmcpa,npolh)
272 CALL my_alloc(dmcpb,npolh)
273 CALL my_alloc(dmcpc,npolh)
274 CALL my_alloc(dmrmw,npolh)
275 CALL my_alloc(elsout,nelt)
276 CALL my_alloc(tmp,nel)
277 CALL my_alloc(poro,nelt)
278 CALL my_alloc(vola,nna)
279 CALL my_alloc(hmin,nntr)
280 CALL my_alloc(vel,nelt)
281 CALL my_alloc(vnod,3,nns)
282 CALL my_alloc(vgrid,3,nntr)
283 CALL my_alloc(svent,nvent)
284 CALL my_alloc(xxx,3,nnt)
285 CALL my_alloc(dmcpd,npolh)
286 CALL my_alloc(dmcpe,npolh)
287 CALL my_alloc(dmcpf,npolh)
288 CALL my_alloc(pel,nelt)
289 CALL my_alloc(dls,npolh)
290 CALL my_alloc(eldmass,nel)
291 CALL my_alloc(eldehpy,nel)
292 CALL my_alloc(vvv,3,nnt)
293
294
295 volph = zero
296 IF(nbgauge > 0 ) THEN
297 ALLOCATE(fvspmd(ifv)%GGG(3,nnt))
298 ALLOCATE(fvspmd(ifv)%GGA(3,nna))
299 ENDIF
300 ALLOCATE(fvspmd(ifv)%AAA(3,nnt))
301 fvspmd(ifv)%AAA(1:3,1:nnt) = zero
302
303 IF (nspmd == 1) THEN
304C needed for partih/on
305 DO i=1,fvspmd(ifv)%NN_L+fvspmd(ifv)%NNI_L
306 i1=fvspmd(ifv)%IBUF_L(1,i)
307 i2=fvspmd(ifv)%IBUF_L(2,i)
308 xxx(1,i1)=x(1,i2)
309 xxx(2,i1)=x(2,i2)
310 xxx(3,i1)=x(3,i2)
311 ENDDO
312C
313 DO i=1,fvspmd(ifv)%NN_L+fvspmd(ifv)%NNI_L
314 i1=fvspmd(ifv)%IBUF_L(1,i)
315 i2=fvspmd(ifv)%IBUF_L(2,i)
316 vvv(1,i1)=v(1,i2)
317 vvv(2,i1)=v(2,i2)
318 vvv(3,i1)=v(3,i2)
319 ENDDO
320C
321 IF (intbag /= 0.AND.nvent > 0.AND.ivolu(39) /= 0) THEN
322 DO i=1,fvspmd(ifv)%NN_L+fvspmd(ifv)%NNI_L
323 i1=fvspmd(ifv)%IBUF_L(1,i)
324 i2=fvspmd(ifv)%IBUF_L(2,i)
325 icont(i1)=icontact(i2)
326 ENDDO
327 ENDIF
328
329 ELSE
330 CALL spmd_fvb_gath(ifv, x, xxx, rbid, rbid, 4)
331C
332 CALL spmd_fvb_gath(ifv, v, vvv, rbid, rbid, 4)
333C
334 IF (intbag /= 0.AND.nvent > 0.AND.ivolu(39) /= 0)
335 . CALL spmd_fvb_igath(ifv, icontact, icont)
336
337C
338 IF (ispmd/=fvspmd(ifv)%PMAIN-1) GOTO 1000
339 END IF
340C
341
342!$OMP PARALLEL PRIVATE(I,J)
343!$OMP+ PRIVATE(I1,I2,KSI,ETA,N1,N2,N3,IEL,X1,Y1,Z1,FAC)
344!$OMP+ PRIVATE(X2,Y2,Z2,X3,Y3,Z3,VX1,VX2,VX3,VY1,VY2,VY3,VZ1,VZ2,VZ3)
345!$OMP+ PRIVATE(X12,Y12,Z12,X13,Y13,Z13,AREA2,NRX,NRY,NRZ,JJ,AREA,KK,K)
346!$OMP+ PRIVATE(NX,NY,NZ,PVOLTMP)
347
348C---------------------------------------------------
349C Normal Vectors & Surfaces
350C---------------------------------------------------
351
352!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
353 DO i=1,nnt
354 nodarea(i)=zero
355 ENDDO
356!$OMP END DO
357!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
358 DO iel=1,nelt
359 fel(iel)=zero
360 ENDDO
361!$OMP END DO
362!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
363 DO iel=1,nelt
364 n1=elem(1,iel)
365 n2=elem(2,iel)
366 n3=elem(3,iel)
367 x1=xxx(1,n1)
368 x2=xxx(1,n2)
369 x3=xxx(1,n3)
370 y1=xxx(2,n1)
371 y2=xxx(2,n2)
372 y3=xxx(2,n3)
373 z1=xxx(3,n1)
374 z2=xxx(3,n2)
375 z3=xxx(3,n3)
376 x12=x2-x1
377 y12=y2-y1
378 z12=z2-z1
379 x13=x3-x1
380 y13=y3-y1
381 z13=z3-z1
382 nrx=y12*z13-z12*y13
383 nry=z12*x13-x12*z13
384 nrz=x12*y13-y12*x13
385 area2=sqrt(nrx**2+nry**2+nrz**2)
386 elarea(iel)=area2
387 norm(1,iel)=nrx/area2
388 norm(2,iel)=nry/area2
389 norm(3,iel)=nrz/area2
390 IF(iel <= nel) THEN
391 tmp(iel) = one_over_6 * (x1*nrx+y1*nry+z1*nrz)
392 ENDIF
393 pel(iel)=zero
394 ENDDO
395!$OMP END DO
396
397!$OMP SINGLE
398 volg=zero
399 DO iel=1,nelt
400 n1=elem(1,iel)
401 n2=elem(2,iel)
402 n3=elem(3,iel)
403 area2 = elarea(iel)
404 elarea(iel) = half * area2
405 nodarea(n1)=nodarea(n1)+one_over_6*area2
406 nodarea(n2)=nodarea(n2)+one_over_6*area2
407 nodarea(n3)=nodarea(n3)+one_over_6*area2
408 IF(iel <= nel) THEN
409 volg=volg+tmp(iel)
410 ENDIF
411 ENDDO
412C
413 IF (dt1==zero.OR.ivolu(39)==0) THEN
414 DO iel=1,nelt
415 elfmass(iel)=zero
416 elfehpy(iel)=zero
417 ENDDO
418 fsav(1:nthvki)=zero
419 ENDIF
420!$OMP END SINGLE
421C---------------------------------------------------
422C Evolution of FV mesh
423C---------------------------------------------------
424C Nodes
425C
426
427
428!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
429 DO i=1,nns
430 IF (ifvnod(1,i) == 1) THEN
431 iel=ifvnod(2,i)
432 ksi=rfvnod(1,i)
433 eta=rfvnod(2,i)
434C
435 n1=elema(1,iel)
436 n2=elema(2,iel)
437 n3=elema(3,iel)
438C
439 IF (tagela(iel) > 0) THEN
440 x1=xxx(1,n1)
441 y1=xxx(2,n1)
442 z1=xxx(3,n1)
443 x2=xxx(1,n2)
444 y2=xxx(2,n2)
445 z2=xxx(3,n2)
446 x3=xxx(1,n3)
447 y3=xxx(2,n3)
448 z3=xxx(3,n3)
449 vx1=vvv(1,n1)
450 vx2=vvv(1,n2)
451 vx3=vvv(1,n3)
452 vy1=vvv(2,n1)
453 vy2=vvv(2,n2)
454 vy3=vvv(2,n3)
455 vz1=vvv(3,n1)
456 vz2=vvv(3,n2)
457 vz3=vvv(3,n3)
458 ELSEIF (tagela(iel) < 0) THEN
459 x1=xxxa(1,n1)
460 y1=xxxa(2,n1)
461 z1=xxxa(3,n1)
462 x2=xxxa(1,n2)
463 y2=xxxa(2,n2)
464 z2=xxxa(3,n2)
465 x3=xxxa(1,n3)
466 y3=xxxa(2,n3)
467 z3=xxxa(3,n3)
468 vx1=vvva(1,n1)
469 vx2=vvva(1,n2)
470 vx3=vvva(1,n3)
471 vy1=vvva(2,n1)
472 vy2=vvva(2,n2)
473 vy3=vvva(2,n3)
474 vz1=vvva(3,n1)
475 vz2=vvva(3,n2)
476 vz3=vvva(3,n3)
477 ENDIF
478 pnod(1,i)=(one-ksi-eta)*x1+ksi*x2+eta*x3
479 pnod(2,i)=(one-ksi-eta)*y1+ksi*y2+eta*y3
480 pnod(3,i)=(one-ksi-eta)*z1+ksi*z2+eta*z3
481 vnod(1,i)=(one-ksi-eta)*vx1+ksi*vx2+eta*vx3
482 vnod(2,i)=(one-ksi-eta)*vy1+ksi*vy2+eta*vy3
483 vnod(3,i)=(one-ksi-eta)*vz1+ksi*vz2+eta*vz3
484C
485 ELSEIF (ifvnod(1,i) == 2) THEN
486 i2=ifvnod(3,i)
487 pnod(1:3,i) = xxxa(1:3,i2)
488 vnod(1:3,i) = vvva(1:3,i2)
489 ENDIF
490 ENDDO
491!$OMP END DO
492
493C internal node
494!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
495 DO i=1,nns
496 IF (ifvnod(1,i) == 3) THEN
497 i1=ifvnod(2,i)
498 i2=ifvnod(3,i)
499 fac=rfvnod(1,i)
500 pnod(1,i)=fac*pnod(1,i1)+(one-fac)*pnod(1,i2)
501 pnod(2,i)=fac*pnod(2,i1)+(one-fac)*pnod(2,i2)
502 pnod(3,i)=fac*pnod(3,i1)+(one-fac)*pnod(3,i2)
503 vnod(1,i)=fac*vnod(1,i1)+(one-fac)*vnod(1,i2)
504 vnod(2,i)=fac*vnod(2,i1)+(one-fac)*vnod(2,i2)
505 vnod(3,i)=fac*vnod(3,i1)+(one-fac)*vnod(3,i2)
506 ENDIF
507 ENDDO
508!$OMP END DO
509C Normal vector, area, grid velocity in triangle center
510!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
511 DO i=1,nntr
512 n1=ifvtri(1,i)
513 n2=ifvtri(2,i)
514 n3=ifvtri(3,i)
515 x1=pnod(1,n1)
516 y1=pnod(2,n1)
517 z1=pnod(3,n1)
518 x2=pnod(1,n2)
519 y2=pnod(2,n2)
520 z2=pnod(3,n2)
521 x3=pnod(1,n3)
522 y3=pnod(2,n3)
523 z3=pnod(3,n3)
524 x12=x2-x1
525 y12=y2-y1
526 z12=z2-z1
527 x13=x3-x1
528 y13=y3-y1
529 z13=z3-z1
530 nrx=y12*z13-z12*y13
531 nry=z12*x13-x12*z13
532 nrz=x12*y13-y12*x13
533 area2=sqrt(nrx**2+nry**2+nrz**2)
534 parea(i)=half*area2
535 IF (area2 > 0) THEN
536 pnorm(1,i)=nrx/area2
537 pnorm(2,i)=nry/area2
538 pnorm(3,i)=nrz/area2
539 ELSE
540 pnorm(1,i)=zero
541 pnorm(2,i)=zero
542 pnorm(3,i)=zero
543 ENDIF
544 IF(ivolu(39)==0) THEN
545 vgrid(1,i) = zero
546 vgrid(2,i) = zero
547 vgrid(3,i) = zero
548 ELSE
549 vx1=vnod(1,n1)
550 vy1=vnod(2,n1)
551 vz1=vnod(3,n1)
552 vx2=vnod(1,n2)
553 vy2=vnod(2,n2)
554 vz2=vnod(3,n2)
555 vx3=vnod(1,n3)
556 vy3=vnod(2,n3)
557 vz3=vnod(3,n3)
558 vgrid(1,i)=third*(vx1+vx2+vx3)
559 vgrid(2,i)=third*(vy1+vy2+vy3)
560 vgrid(3,i)=third*(vz1+vz2+vz3)
561 ENDIF
562 ENDDO
563!$OMP END DO
564C Polyhedra centers
565
566!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
567 DO i=1,npolh
568 pvolu(i)= zero
569 pvoltmp = zero
570 area_max(i) = zero
571 DO j=ifvpadr(i),ifvpadr(i+1)-1 !polygons composing the poyhedron
572 jj=ifvpolh(j)
573 DO k=ifvtadr(jj), ifvtadr(jj+1)-1 !triangle composing the polygon
574 kk=ifvpoly(k)
575 area=parea(kk) !triangle area7
576 area_max(i)=max(area_max(i), area)
577 iel=ifvtri(4,kk)
578 IF (iel > 0) THEN
579 nx=pnorm(1,kk)
580 ny=pnorm(2,kk)
581 nz=pnorm(3,kk)
582 ELSE
583 i1 =ifvtri(5,kk)
584 i2 =ifvtri(6,kk)
585 IF (i1 == i) THEN
586 nx=pnorm(1,kk)
587 ny=pnorm(2,kk)
588 nz=pnorm(3,kk)
589 ENDIF
590 IF (i2 == i) THEN
591 nx=-pnorm(1,kk)
592 ny=-pnorm(2,kk)
593 nz=-pnorm(3,kk)
594 ENDIF
595 IF (i1 == i.AND.i2 == i) THEN
596 nx=zero
597 ny=zero
598 nz=zero
599 ENDIF
600 ENDIF
601 n1=ifvtri(1,kk)
602 x1=pnod(1,n1)
603 y1=pnod(2,n1)
604 z1=pnod(3,n1)
605 pvoltmp=pvoltmp+third*area*(x1*nx+y1*ny+z1*nz)
606
607 enddo! next K(triangle)
608 enddo! next J (polygon)
609 pvolu(i) = pvoltmp
610 enddo! next I (polyhedron)
611!$OMP END DO
612!$OMP END PARALLEL
613C---------------------------------------------------
614C Verifications
615C---------------------------------------------------
616 areap=zero
617 DO i=1,nntr
618 iel=ifvtri(4,i)
619 IF (iel > 0) areap=areap+parea(i)
620 ENDDO
621 areael=zero
622 DO iel=1,nel
623 areael=areael+elarea(iel)
624 ENDDO
625 rvolu(18)=areael
626 ilvout=ivolu(44)
627 volph=zero
628 exit_neg_vol = .false.
629 DO i=1,npolh
630 volph=volph+pvolu(i)
631 IF(pvolu(i) <= 0) exit_neg_vol = .true.
632 ENDDO
633
634 ipri=mod(ncycle,iabs(ncpri))
635 IF(exit_neg_vol) THEN
636 DO i=1,npolh
637 IF (pvolu(i) <= zero) THEN
638 info=1
639 ierr=ierr+1
640 CALL ancmsg(msgid=185,anmode=aninfo,
641 . i1=idpolh(i),r1=pvolu(i),i3=i,i4=npolh)
642 CALL arret(2)
643 ENDIF
644 ENDDO
645 ENDIF
646
647
648C---------------------------------------------------
649C Calcul des quantites dans les polyhedres
650C---------------------------------------------------
651 IF (dt1 == zero) THEN
652 gamai=rvolu(1)
653 gpolh(1:npolh)=gamai
654 cpai =rvolu(7)
655 cpapolh(1:npolh)=cpai
656 cpbi =rvolu(8)
657 cpbpolh(1:npolh)=cpbi
658 cpci =rvolu(9)
659 cpcpolh(1:npolh)=cpci
660 rmwi =rvolu(10)
661 rmwpolh(1:npolh)=rmwi
662 ti =rvolu(13)
663 tpolh(1:npolh)=ti
664 cpdi =rvolu(56)
665 cpdpolh(1:npolh)=cpdi
666 cpei =rvolu(57)
667 cpepolh(1:npolh)=cpei
668 cpfi =rvolu(58)
669 cpfpolh(1:npolh)=cpfi
670 rhoi =rvolu(62)
671 rpolh(1:npolh)=rhoi
672 ENDIF
673 IF (dt1==zero.OR.ivolu(39)==0) THEN
674 rhoi =rvolu(62)
675 efac =rvolu(66)
676 mpolh(1:npolh)=rhoi*pvolu(1:npolh)
677 epolh(1:npolh)=mpolh(1:npolh)*efac
678 pu(1,1:npolh) = rvolu(67)
679 pu(2,1:npolh) = rvolu(68)
680 pu(3,1:npolh) = rvolu(69)
681 qpolh(1,1:npolh) = mpolh(1:npolh)*pu(1,1:npolh)
682 qpolh(2,1:npolh) = mpolh(1:npolh)*pu(2,1:npolh)
683 qpolh(3,1:npolh) = mpolh(1:npolh)*pu(3,1:npolh)
684 ENDIF
685
686 pext =rvolu(3)
687 datainj(1:6,1:njet)=zero
688 IF(ivolu(39) == 0) GO TO 250
689
690C---------------------------------------------------
691C Injecteurs
692C---------------------------------------------------
693 DO iel=1,nelt
694 IF (itagel(iel) > 0) THEN
695 iinj=itagel(iel)
696 datainj(1,iinj)=datainj(1,iinj)+elarea(iel)
697 ENDIF
698 ENDDO
699C
700 scalt =rvolu(26)
701 IF(ityp == 6) THEN
702 CALL fvinjt6(njet, ibagjet, rbagjet, npc, tf, nsensor,
703 . sensor_tab,scalt, datainj, python )
704 ELSEIF(ityp == 8) THEN
705 CALL fvinjt8(njet, ibagjet, rbagjet, npc, tf, nsensor,
706 . sensor_tab,scalt, igeo, geo, pm, ivolu, datainj,python)
707 ENDIF
708C
709 DO iinj=1,njet
710 isens=ibagjet(4,iinj)
711 fvel=rbagjet(15,iinj)
712 ivel=ibagjet(11,iinj)
713 IF(isens == 0)THEN
714 tstart=zero
715 ELSE
716 tstart=sensor_tab(isens)%TSTART
717 ENDIF
718C ITTF : 1 shift event time (can be removed in a later version)
719C ITTF : 2 shift event time + event function
720C ITTF : 1 or 2 +10 1 flag has been activated and RVOLU(60) is filled
721 ittf=ivolu(17)
722 IF (tt >= tstart.AND.(ittf == 1.OR.ittf == 2.OR.ittf == 3))THEN
723 ittf=ittf+10
724 rvolu(60)=tstart
725 ivolu(17)=ittf
726 END IF
727 IF (tt >= tstart.AND.dt1 > zero)THEN
728 tsg=(tt-tstart)*scalt
729 IF (ivel > 0) THEN
730 injvel=fvel*finter_mixed(python,nfunct,ivel,tsg,npc,tf)
731 ELSE
732 injvel = fvel
733 ENDIF
734 ELSE
735 injvel=zero
736 ENDIF
737 datainj(3,iinj)=injvel
738 ENDDO
739C---------------------------------------------------
740C Vent holes
741C---------------------------------------------------
742 IF (intbag == 0) THEN
743 poro(1:nelt)=zero
744 ELSE
745!$OMP PARALLEL PRIVATE(IEL,N1,N2,N3)
746!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
747 DO iel=1,nelt
748 n1=elem(1,iel)
749 n2=elem(2,iel)
750 n3=elem(3,iel)
751 poro(iel)=zero
752 IF (icont(n1) /= 0) poro(iel)=poro(iel)+third
753 IF (icont(n2) /= 0) poro(iel)=poro(iel)+third
754 IF (icont(n3) /= 0) poro(iel)=poro(iel)+third
755 ENDDO
756!$OMP END DO
757!$OMP END PARALLEL
758 ENDIF
759C
760 pext =rvolu(3)
761 scalt=rvolu(26)
762 scalp=rvolu(27)
763 scals=rvolu(28)
764 ttf =rvolu(60)
765C-----------------------------------------------------------
766C Surface de fuite effective des triangles de l'airbag
767C Correction des scale factors sur la surface des vent holes
768C-----------------------------------------------------------
769 CALL fvvent0(
770 1 elsout ,aoutot ,nvent ,nelt ,ittf ,
771 2 elarea ,elsini ,elem ,itagel ,svent ,
772 3 ibaghol ,rvolu ,rbaghol ,poro ,p ,
773 4 eltg ,iparg ,mattg ,nel ,porosity,
774 5 ipm ,pm ,elbuf_tab,igroupc ,igrouptg)
775C---------------------------------------------------
776C Bilans volumes finis
777C---------------------------------------------------
778!$OMP PARALLEL PRIVATE(I,J,JJ,K,KK,IEL,JEL,AREA,N1,N2,N3)
779!$OMP+ PRIVATE(X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,X12,Y12,Z12,X23,Y23,Z23,Z31,Y31,X31)
780!$OMP+ PRIVATE(L12,L23,L31,H1)
781!$OMP SECTIONS
782!$OMP SECTION
783 fel(1:nelt) = zero
784 elfvel(1:nelt)=zero
785!$OMP SECTION
786 eldmass(1:nel)=zero
787 eldehpy(1:nel)=zero
788!$OMP SECTION
789 itagp(1:npolh)=0
790 dm(1:npolh)=zero
791!$OMP SECTION
792 dq(1,1:npolh)=zero
793 dq(2,1:npolh)=zero
794!$OMP SECTION
795 dq(3,1:npolh)=zero
796 de(1:npolh)=zero
797!$OMP SECTION
798 dmini(1:npolh)=zero
799 dmcpa(1:npolh)=zero
800!$OMP SECTION
801 dmcpb(1:npolh)=zero
802 dmcpc(1:npolh)=zero
803!$OMP SECTION
804 dmrmw(1:npolh)=zero
805 dmcpd(1:npolh)=zero
806!$OMP SECTION
807 dmcpe(1:npolh)=zero
808 dmcpf(1:npolh)=zero
809!$OMP SECTION
810C----------------------
811C Transports
812C----------------------
813 pu(1,1:npolh)=qpolh(1,1:npolh)/mpolh(1:npolh)
814 pu(2,1:npolh)=qpolh(2,1:npolh)/mpolh(1:npolh)
815 pu(3,1:npolh)=qpolh(3,1:npolh)/mpolh(1:npolh)
816!$OMP END SECTIONS
817!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
818 DO i=1,npolh
819 DO j=ifvpadr(i),ifvpadr(i+1)-1
820 jj=ifvpolh(j)
821 DO k=ifvtadr(jj),ifvtadr(jj+1)-1
822 kk=ifvpoly(k)
823 iel=abs(ifvtri(4,kk))
824 jel=ifvtri(4,kk)
825 area=parea(kk)
826 IF(jel == 0) THEN
827C Permeabilite du triangle interne
828 n1=ifvtri(1,kk)
829 n2=ifvtri(2,kk)
830 n3=ifvtri(3,kk)
831 x1=pnod(1,n1)
832 y1=pnod(2,n1)
833 z1=pnod(3,n1)
834 x2=pnod(1,n2)
835 y2=pnod(2,n2)
836 z2=pnod(3,n2)
837 x3=pnod(1,n3)
838 y3=pnod(2,n3)
839 z3=pnod(3,n3)
840 x12=x2-x1
841 y12=y2-y1
842 z12=z2-z1
843 x23=x3-x2
844 y23=y3-y2
845 z23=z3-z2
846 x31=x1-x3
847 y31=y1-y3
848 z31=z1-z3
849 l12=x12**2+y12**2+z12**2
850 l23=x23**2+y23**2+z23**2
851 l31=x31**2+y31**2+z31**2
852 h1=four*area**2/max(l12,l23,l31)
853 hmin(kk)=sqrt(h1)
854 ENDIF
855 ENDDO
856 ENDDO
857 ENDDO
858!$OMP END PARALLEL
859C----------------
860C Calcul des flux
861C----------------
862 DO i=1,npolh
863 DO j=ifvpadr(i),ifvpadr(i+1)-1
864 jj=ifvpolh(j)
865 DO k=ifvtadr(jj),ifvtadr(jj+1)-1
866 kk =ifvpoly(k)
867 iel =abs(ifvtri(4,kk))
868 jel =ifvtri(4,kk)
869 area=parea(kk)
870 IF (area == zero) cycle
871 booltagel = .false.
872 IF(iel > 0)THEN
873 IF(itagel(iel) > 0)booltagel = .true.
874 ENDIF
875 ! IF (JEL > 0.OR.(JEL < 0.AND.ITAGEL(IEL) > 0)) THEN debugger issue in case of IEL=0 (out of bounds)
876 IF (jel > 0.OR.(jel < 0.AND.booltagel)) THEN
877 ii=i
878 nrx=norm(1,iel)
879 nry=norm(2,iel)
880 nrz=norm(3,iel)
881 nv = 0
882 IF(jel < 0) THEN
883 IF (ifvtri(5,kk) == i) THEN
884 nv=ifvtri(6,kk)
885 nx=pnorm(1,kk)
886 ny=pnorm(2,kk)
887 nz=pnorm(3,kk)
888 ELSEIF(ifvtri(6,kk) == i) THEN
889 nv=ifvtri(5,kk)
890 nx=-pnorm(1,kk)
891 ny=-pnorm(2,kk)
892 nz=-pnorm(3,kk)
893 ENDIF
894 !check if adjacent polyhedron already updated
895 IF (nv < i) GOTO 100
896 fac=nx*nrx+ny*nry+nz*nrz
897 IF(fac < 0) ii=nv
898 ENDIF
899C--------------------------------------------
900C Triangle on airbag envelope
901C internal triangle is an injector
902C--------------------------------------------
903 IF (itagel(iel) > 0) THEN
904C Injector
905 iinj=itagel(iel)
906 dmi=datainj(2,iinj)*area/datainj(1,iinj)
907 IF (nv == i .AND. jel < 0) dmi=dmi*half ! le triangle est traite deux fois
908 dqi(1)=-dmi*datainj(3,iinj)*nrx
909 dqi(2)=-dmi*datainj(3,iinj)*nry
910 dqi(3)=-dmi*datainj(3,iinj)*nrz
911C inj. enthalpy DEI=DMI*DATAINJ(4,IINJ)/DATAINJ(5,IINJ)
912 dei=dmi*datainj(4,iinj)
913C
914 dm(ii)=dm(ii)+dmi
915 dq(1,ii)=dq(1,ii)+dqi(1)
916 dq(2,ii)=dq(2,ii)+dqi(2)
917 dq(3,ii)=dq(3,ii)+dqi(3)
918 de(ii)=de(ii)+dei
919C gas mixture
920 dmcpa(ii)=dmcpa(ii)+dmi*rbagjet(2,iinj)
921 dmcpb(ii)=dmcpb(ii)+dmi*rbagjet(3,iinj)
922 dmcpc(ii)=dmcpc(ii)+dmi*rbagjet(4,iinj)
923 dmrmw(ii)=dmrmw(ii)+dmi*rbagjet(1,iinj)
924 IF(ityp == 8) THEN
925 dmcpd(ii)=dmcpd(ii)+dmi*rbagjet(16,iinj)
926 dmcpe(ii)=dmcpe(ii)+dmi*rbagjet(17,iinj)
927 dmcpf(ii)=dmcpf(ii)+dmi*rbagjet(18,iinj)
928 ENDIF
929C
930 elfmass(iel)=elfmass(iel)+dmi*dt1
931 elfehpy(iel)=elfehpy(iel)+dei*dt1
932 elfvel(iel) =-datainj(3,iinj)
933 ELSEIF (itagel(iel) < 0) THEN
934C Vent hole
935 ivent =-itagel(iel)
936 idef =ibaghol(1,ivent)
937 itvent =ibaghol(10,ivent)
938C
939 nx=norm(1,iel)
940 ny=norm(2,iel)
941 nz=norm(3,iel)
942 aout=elsout(iel)*area/elarea(iel)
943 p1=ppolh(i)
944 rho1=rpolh(i)
945 re1=epolh(i)/pvolu(i)
946 gama1=gpolh(i)
947 ux1=pu(1,i)
948 uy1=pu(2,i)
949 uz1=pu(3,i)
950 uu =zero
951C LOCAL VELOCITY
952 IF (itvent == 1) THEN
953 uu=nx*ux1+ny*uy1+nz*uz1
954 rho2=rho1
955 IF(p1 < pext) uu=zero
956C BERNOULLI
957 ELSEIF ((itvent==2.AND.p1>zero).OR.(itvent==4.AND.p1>=pext)) THEN
958 vv=nx*ux1+ny*uy1+nz*uz1
959 vv=max(vv,zero)
960 eta=(gama1-one)/gama1
961 ksi=one/eta
962 pstag=p1*(one+half*eta*rho1*vv*vv/p1)**ksi
963 pcrit=pstag*(two/(gama1+one))**ksi
964 p2 = max(pext,pcrit)
965 rho2 =rho1*(p2/p1)**(one/gama1)
966 uu=two*p1*(one-(p2/p1)**eta)/(rho1*eta)
967 uu=max(uu,zero)
968 uu=sqrt(uu)
969 vmax=half*((p1-pext)*pvolu(i)/(gama1-one))
970 . /max(em20,rho2*aout*dt1*gama1*re1/rho1)
971 vmax=max(vmax,zero)
972 uu=min(uu,vmax)
973C CHEMKIN
974 ELSEIF (itvent == 3) THEN
975 p2=p1-pext
976 ivdp=ibaghol(9,ivent)
977 fvdp=rbaghol(13,ivent)
978 uu=fvdp*get_u_func(ivdp,p2*scalp,deri)
979 IF(p2 < zero) THEN
980 rho2=rvolu(62)
981 ELSE
982 rho2=rho1
983 ENDIF
984C BERNOULLI with flow-in
985 ELSEIF (itvent==4.AND.p1<pext) THEN
986 gamai=rvolu(1)
987 rhoi=rvolu(62)
988 eta=(gamai-one)/gamai
989 pcrit=pext*(two/(gamai+one))**(one/eta)
990 p2 = max(p1,pcrit)
991 rho2 =rhoi*(p2/pext)**(one/gamai)
992 uu=two*pext*(one-(p2/pext)**eta)/(rhoi*eta)
993 uu=max(uu,zero)
994 uu=sqrt(uu)
995 vmax=half*((pext-p1)*pvolu(i)/(gama1-one))
996 . /max(em20,rho2*aout*dt1*rvolu(63))
997 vmax=max(vmax,zero)
998 uu=min(uu,vmax)
999 uu=-uu
1000C GRAEFE
1001 ELSEIF (itvent==5) THEN
1002 p2=p1-pext
1003 rho2=rho1
1004 uu=two*p2/rho2
1005 uu=max(uu,zero)
1006 uu=sqrt(uu)
1007C
1008 ENDIF
1009C
1010 IF (uu > zero.AND.idef == 1) THEN
1011 dmout=uu*rho2*aout
1012 dm(i)=dm(i)-dmout
1013 dq(1,i)=dq(1,i)-dmout*ux1
1014 dq(2,i)=dq(2,i)-dmout*uy1
1015 dq(3,i)=dq(3,i)-dmout*uz1
1016 dhout=dmout*gama1*re1/rho1
1017 de(i)=de(i)-dhout
1018 !DMINI(I)=DMINI(I)-DMOUT
1019 cpai=cpapolh(i)
1020 cpbi=cpbpolh(i)
1021 cpci=cpcpolh(i)
1022 rmwi=rmwpolh(i)
1023 dmcpa(i)=dmcpa(i)-dmout*cpai
1024 dmcpb(i)=dmcpb(i)-dmout*cpbi
1025 dmcpc(i)=dmcpc(i)-dmout*cpci
1026 dmrmw(i)=dmrmw(i)-dmout*rmwi
1027 IF(ityp == 8) THEN
1028 cpdi=cpdpolh(i)
1029 cpei=cpepolh(i)
1030 cpfi=cpfpolh(i)
1031 dmcpd(i)=dmcpd(i)-dmout*cpdi
1032 dmcpe(i)=dmcpe(i)-dmout*cpei
1033 dmcpf(i)=dmcpf(i)-dmout*cpfi
1034 ENDIF
1035 elfmass(iel)=elfmass(iel)-dmout*dt1
1036 elfehpy(iel)=elfehpy(iel)-dhout*dt1
1037 elfvel(iel)=elfvel(iel)+uu*area/elarea(iel)
1038 eldmass(iel)=-dmout
1039 eldehpy(iel)=-dhout
1040 ENDIF
1041C
1042 IF (uu < zero.AND.idef == 1) THEN
1043 dmi=-uu*rho2*aout
1044 dm(i)=dm(i)+dmi
1045 dq(1,i)=dq(1,i)+dmi*ux1
1046 dq(2,i)=dq(2,i)+dmi*uy1
1047 dq(3,i)=dq(3,i)+dmi*uz1
1048 cpai=rvolu(7)
1049 cpbi=rvolu(8)
1050 cpci=rvolu(9)
1051 rmwi=rvolu(10)
1052 dmcpa(i)=dmcpa(i)+dmi*cpai
1053 dmcpb(i)=dmcpb(i)+dmi*cpbi
1054 dmcpc(i)=dmcpc(i)+dmi*cpci
1055 dmrmw(i)=dmrmw(i)+dmi*rmwi
1056 IF(ityp == 8) THEN
1057 cpdi=rvolu(56)
1058 cpei=rvolu(57)
1059 cpfi=rvolu(58)
1060 dmcpd(i)=dmcpd(i)+dmi*cpdi
1061 dmcpe(i)=dmcpe(i)+dmi*cpei
1062 dmcpf(i)=dmcpf(i)+dmi*cpfi
1063 ENDIF
1064 de(i)=de(i)+dmi*rvolu(63)
1065 elfmass(iel)=elfmass(iel)+dmi*dt1
1066 elfehpy(iel)=elfehpy(iel)+dmi*dt1*rvolu(63)
1067 elfvel(iel) =elfvel(iel) +uu*area/elarea(iel)
1068 eldmass(iel)=dmi
1069 eldehpy(iel)=dmi*rvolu(63)
1070 ENDIF
1071 ENDIF
1072 ELSE
1073C----------------------------------------------------------------
1074C internal triangle is not an injector => communication between polyhedra
1075C----------------------------------------------------------------
1076 nv = 0
1077 IF (ifvtri(5,kk) == i) THEN
1078 nv=ifvtri(6,kk)
1079 nx=pnorm(1,kk)
1080 ny=pnorm(2,kk)
1081 nz=pnorm(3,kk)
1082 ELSEIF (ifvtri(6,kk) == i) THEN
1083 nv=ifvtri(5,kk)
1084 nx=-pnorm(1,kk)
1085 ny=-pnorm(2,kk)
1086 nz=-pnorm(3,kk)
1087 ENDIF
1088 !check if adjacent polyhedron already updated
1089 IF (nv <= i) GOTO 100
1090 rho1=rpolh(i)
1091 re1=epolh(i)/pvolu(i)
1092 ux1=pu(1,i)
1093 uy1=pu(2,i)
1094 uz1=pu(3,i)
1095 gama1=gpolh(i)
1096 rho2=rpolh(nv)
1097 re2=epolh(nv)/pvolu(nv)
1098 ux2=pu(1,nv)
1099 uy2=pu(2,nv)
1100 uz2=pu(3,nv)
1101 gama2=gpolh(nv)
1102C
1103 vfx=half*(ux1+ux2)
1104 vfy=half*(uy1+uy2)
1105 vfz=half*(uz1+uz2)
1106 IF(jel == 0) THEN
1107C Permeabilite du triangle interne
1108 hm1=hmin(kk)
1109 h1 = one-rvolu(51)/hm1
1110 IF(h1 > 0) THEN
1111 area=area*h1
1112 ELSE
1113 area = zero
1114 ENDIF
1115
1116C AREA=AREA*MAX(ZERO,ONE-RVOLU(51)/HM1)
1117 ELSE
1118C Triangle appuy sur structure interne => flux fonction de la porosite
1119 IF(porosity(iel-nel) == zero) GOTO 50
1120 area=area*porosity(iel-nel)
1121 ENDIF
1122C
1123 ss=nx*(vfx-vgrid(1,kk))+ny*(vfy-vgrid(2,kk))
1124 . +nz*(vfz-vgrid(3,kk))
1125C Full upwind
1126 IF (ss <= zero) THEN
1127 alpha=zero
1128 ELSE
1129 alpha=one
1130 ENDIF
1131C
1132 rhom=alpha*rho1+(one-alpha)*rho2
1133 ruxm=alpha*rho1*ux1+(one-alpha)*rho2*ux2
1134 ruym=alpha*rho1*uy1+(one-alpha)*rho2*uy2
1135 ruzm=alpha*rho1*uz1+(one-alpha)*rho2*uz2
1136 rem=alpha*gama1*re1+(one-alpha)*gama2*re2
1137C
1138 massflow = rhom*ss*area
1139 dm(i)=dm(i)-massflow
1140 dq(1,i)=dq(1,i)-ruxm*ss*area
1141 dq(2,i)=dq(2,i)-ruym*ss*area
1142 dq(3,i)=dq(3,i)-ruzm*ss*area
1143 de(i)=de(i)-rem*ss*area
1144 dm(nv)=dm(nv)+massflow
1145 dq(1,nv)=dq(1,nv)+ruxm*ss*area
1146 dq(2,nv)=dq(2,nv)+ruym*ss*area
1147 dq(3,nv)=dq(3,nv)+ruzm*ss*area
1148 de(nv)=de(nv)+rem*ss*area
1149 cpam = alpha * cpapolh(i) + (one - alpha) * cpapolh(nv)
1150 cpbm = alpha * cpbpolh(i) + (one - alpha) * cpbpolh(nv)
1151 cpcm = alpha * cpcpolh(i) + (one - alpha) * cpcpolh(nv)
1152 rmwm = alpha * rmwpolh(i) + (one - alpha) * rmwpolh(nv)
1153 dmcpa(i) = dmcpa(i) - massflow * cpam
1154 dmcpa(nv) = dmcpa(nv) + massflow * cpam
1155 dmcpb(i) = dmcpb(i) - massflow * cpbm
1156 dmcpb(nv) = dmcpb(nv) + massflow * cpbm
1157 dmcpc(i) = dmcpc(i) - massflow * cpcm
1158 dmcpc(nv) = dmcpc(nv) + massflow * cpcm
1159 dmrmw(i) = dmrmw(i) - massflow * rmwm
1160 dmrmw(nv) = dmrmw(nv) + massflow * rmwm
1161 IF (ityp == 8) THEN
1162 cpdm = alpha * cpdpolh(i) + (one - alpha) * cpdpolh(nv)
1163 cpem = alpha * cpepolh(i) + (one - alpha) * cpepolh(nv)
1164 cpfm = alpha * cpfpolh(i) + (one - alpha) * cpfpolh(nv)
1165 dmcpd(i) = dmcpd(i) - massflow * cpdm
1166 dmcpd(nv) = dmcpd(nv) + massflow * cpdm
1167 dmcpe(i) = dmcpe(i) - massflow * cpem
1168 dmcpe(nv) = dmcpe(nv) + massflow * cpem
1169 dmcpf(i) = dmcpf(i) - massflow * cpfm
1170 dmcpf(nv) = dmcpf(nv) + massflow * cpfm
1171 ENDIF
1172C
1173 IF(jel < 0) THEN
1174 elfmass(iel)=elfmass(iel)+massflow*dt1
1175 elfvel(iel) =elfvel(iel) +ss*area/elarea(iel)
1176 ENDIF
1177C gas mixture
1178c$$$ IF (SS > ZERO) THEN
1179c$$$ DMINI(I)=DMINI(I)-MASSFLOW
1180c$$$ DMCPA(NV)=DMCPA(NV)+MASSFLOW*CPAPOLH(I)
1181c$$$ DMCPB(NV)=DMCPB(NV)+MASSFLOW*CPBPOLH(I)
1182c$$$ DMCPC(NV)=DMCPC(NV)+MASSFLOW*CPCPOLH(I)
1183c$$$ DMRMW(NV)=DMRMW(NV)+MASSFLOW*RMWPOLH(I)
1184c$$$ ELSE
1185c$$$ DMINI(NV)=DMINI(NV)+MASSFLOW
1186c$$$ DMCPA(I)=DMCPA(I)-MASSFLOW*CPAPOLH(NV)
1187c$$$ DMCPB(I)=DMCPB(I)-MASSFLOW*CPBPOLH(NV)
1188c$$$ DMCPC(I)=DMCPC(I)-MASSFLOW*CPCPOLH(NV)
1189c$$$ DMRMW(I)=DMRMW(I)-MASSFLOW*RMWPOLH(NV)
1190c$$$ ENDIF
1191c$$$ IF(ITYP == 8) THEN
1192c$$$ IF (SS > ZERO) THEN
1193c$$$ DMCPD(NV)=DMCPD(NV)+MASSFLOW*CPDPOLH(I)
1194c$$$ DMCPE(NV)=DMCPE(NV)+MASSFLOW*CPEPOLH(I)
1195c$$$ DMCPF(NV)=DMCPF(NV)+MASSFLOW*CPFPOLH(I)
1196c$$$ ELSE
1197c$$$ DMCPD(I)=DMCPD(I)-MASSFLOW*CPDPOLH(NV)
1198c$$$ DMCPE(I)=DMCPE(I)-MASSFLOW*CPEPOLH(NV)
1199c$$$ DMCPF(I)=DMCPF(I)-MASSFLOW*CPFPOLH(NV)
1200c$$$ ENDIF
1201c$$$ ENDIF
1202 ENDIF
1203 50 ENDDO ! next triangle
1204 100 ENDDO ! next polygon
1205 ENDDO ! next polyhedron
1206C------------------------
1207C
1208
1209!$OMP PARALLEL PRIVATE(I,J,JJ,MSINI,CPA,CPB,CPC,CPD,CPE,CPF,EFAC,TEMP0)
1210!$OMP+ PRIVATE(NTRIA,NTRIA_TOT,NPOLYG,WSURF,SSP1,SSP2,SSPt)
1211!$OMP+ PRIVATE(TEMP,CP,CV,RMW,Denom,Denom_max,VMAT,KK, AREA,IEL1,IEL2)
1212!$OMP+ PRIVATE(AL,DD,AD,GAMA,SSP,QX,VV,TMP_VEC)
1213!$OMP+ PRIVATE(ITYPL, TEMP2)
1214
1215 itypl = ityp
1216!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1217 DO i=1,npolh
1218 itagp(i)=0
1219 msini=mpolh(i)
1220 mpolh(i)=mpolh(i)+dm(i)*dt1
1221 qpolh(1,i)=qpolh(1,i)+dq(1,i)*dt1
1222 qpolh(2,i)=qpolh(2,i)+dq(2,i)*dt1
1223 qpolh(3,i)=qpolh(3,i)+dq(3,i)*dt1
1224 epolh(i)=epolh(i)+de(i)*dt1
1225 dq(1,i)=zero
1226 dq(2,i)=zero
1227 dq(3,i)=zero
1228 de(i)=zero
1229 IF(mpolh(i) <= zero.OR.epolh(i) <= zero) THEN
1230 tpolh(i)=zero
1231 ELSE
1232C gas mixture
1233 msini=msini+dmini(i)*dt1
1234 cpa=(msini*cpapolh(i)+dmcpa(i)*dt1)/mpolh(i)
1235 cpb=(msini*cpbpolh(i)+dmcpb(i)*dt1)/mpolh(i)
1236 cpc=(msini*cpcpolh(i)+dmcpc(i)*dt1)/mpolh(i)
1237 rmw=(msini*rmwpolh(i)+dmrmw(i)*dt1)/mpolh(i)
1238 efac=epolh(i)/mpolh(i)
1239 IF(ityp == 8) THEN
1240 cpd=(msini*cpdpolh(i)+dmcpd(i)*dt1)/mpolh(i)
1241 cpe=(msini*cpepolh(i)+dmcpe(i)*dt1)/mpolh(i)
1242 cpf=(msini*cpfpolh(i)+dmcpf(i)*dt1)/mpolh(i)
1243 ENDIF
1244C-------------------------
1245C Temperature
1246C-------------------------
1247 temp0=rvolu(25)
1248 temp = temp0
1249 CALL fvtemp(itypl , efac , cpa , cpb , cpc ,
1250 . cpd , cpe , cpf , rmw , temp0,
1251 . temp )
1252 cp=cpa+cpb*temp+cpc*temp*temp
1253 cpapolh(i)=cpa
1254 cpbpolh(i)=cpb
1255 cpcpolh(i)=cpc
1256 IF(ityp == 8) THEN
1257 cpdpolh(i)=cpd
1258 cpepolh(i)=cpe
1259 cpfpolh(i)=cpf
1260 cp=cp+cpd*temp**3+cpf*temp**4
1261 temp2 = temp * temp
1262 IF(cpe /= zero .AND. temp2 > zero) THEN
1263 cp=cp+cpe/(temp2)
1264 ENDIF
1265 ENDIF
1266 cv=double_to_my_real(cp)-rmw
1267 rmwpolh(i)=rmw
1268 IF(cp > huge(cv) .OR. cp < -huge(cv) .OR. cv == zero) THEN
1269 ! abs(CP) == Infinity
1270 gpolh(i)=one
1271 ELSE
1272 gpolh(i)=double_to_my_real(cp)/cv
1273 ENDIF
1274 tpolh(i)=temp
1275C---------------------------------------------------
1276C Pressure & Density
1277C---------------------------------------------------
1278 rpolh(i)=mpolh(i)/pvolu(i)
1279 ppolh(i)=rpolh(i)*rmw*temp
1280 ENDIF
1281 ENDDO
1282!$OMP END DO NOWAIT
1283
1284C---------------------------------------------------
1285C TIME STEP & ARTIFICIAL VISCOSITY
1286C---------------------------------------------------
1287 dtx=ep30
1288 phdt=0
1289 qa=rvolu(46)
1290 qb=rvolu(47)
1291
1292
1293 ! /DT/FVMBAG/[ID_DT_OPTION]
1294 ! ID_DT_OPTION = 0 obsolete
1295 ! ID_DT_OPTION = 1 (legacy method), apply to all FVMBAG1
1296 ! ID_DT_OPTION = 2 (2022.3), apply to given FVMBAG1 or FVMBAG2.
1297 ! in this case characteristic length is calculated depending on L_TYPE flag
1298 ! and smoothing is applied : dt <- min(dt, dt_old*(1+lambda))
1299 !--------!-------------------------------------------------------------------!
1300 ! VALUE ! DESCRIPTION !
1301 !------------!--------!-------------------------------------------------------------------!
1302 ! ! 0,1 ! Constant (Starter minimum characteristic length) !
1303 ! L_TYPE ! 2 ! VOLUME^1/3 at current time !
1304 ! ! 3 ! VOLUME / SMAX at current time !
1305 !------------!--------!-------------------------------------------------------------------!
1306
1307
1308 !------------------------------!
1309 idt_fvmbag = fvdata(ifv)%ID_DT_OPTION ! /DT/FVMBAG/[IDT_FVMBAG]
1310
1311 SELECT CASE(idt_fvmbag)
1312
1313 ! /DT/FVMBAG/1
1314 CASE(0,1)
1315 ! ------------
1316 param_l_type = 0
1317
1318!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1319 DO i=1,npolh
1320 IF (ibpolh(i) == 0 .AND.idtmin(52) < 2) THEN
1321 dls(i)=dlh ! constant length calculated during starter
1322 ELSE
1323 dls(i)=pvolu(i)**third
1324 ENDIF
1325 ENDDO
1326!$OMP END DO
1327
1328 ! /DT/FVMBAG/2
1329 CASE(2)
1330 ! ------------
1331 param_l_type = fvdata(ifv)%L_TYPE
1332
1333 SELECT CASE(param_l_type)
1334
1335 ! /DT/FVMBAG/2 +L_TYPE=0,1
1336 CASE(0,1) !dl = starter constant
1337!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1338 DO i=1,npolh
1339 dls(i)=dlh
1340 ENDDO
1341!$OMP END DO
1342
1343 ! /DT/FVMBAG/2 +L_TYPE=2
1344 CASE(2) ! dl=VOL**1/3
1345!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1346 DO i=1,npolh
1347 dls(i)=pvolu(i)**third
1348 ENDDO
1349!$OMP END DO
1350
1351 ! /DT/FVMBAG/2 +L_TYPE=3
1352 CASE(3)
1353!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1354 DO i=1,npolh
1355 dls(i)=pvolu(i)
1356 ENDDO
1357!$OMP END DO
1358
1359 END select!PARAM_L_TYPE
1360
1361
1362 END select!IDT_FVMBAG
1363
1364 !------------------------------!
1365 ! --- TIME STEP ---!
1366 !------------------------------!
1367 SELECT CASE (idt_fvmbag) ! IDT_FVMBAG=FVDATA(IFV)%ID_DT_OPTION
1368
1369 CASE(0,1) !/DT/FVMBAG/0 & /DT/FVMBAG/1
1370!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1371 DO i=1,npolh
1372 IF(mpolh(i)<=zero.OR.epolh(i)<=zero.OR.gpolh(i)<=one) THEN
1373 dtpolh(i)=ep30
1374 qvisc(i) =zero
1375 ELSE
1376 al=pvolu(i)**third
1377 dd=dm(i)/mpolh(i)
1378 ad=max(zero,dd)
1379 gama=gpolh(i)
1380 ssp=sqrt((gama-one)*gama*epolh(i)/mpolh(i))
1381 qvisc(i)=rpolh(i)*ad*al*(qa*qa*al*ad+qb*ssp)
1382 qx=qb*ssp+al*qa*qa*ad
1383 ssp=qx+sqrt(qx*qx+ssp*ssp)
1384 vv=sqrt(pu(1,i)**2+pu(2,i)**2+pu(3,i)**2)
1385 dtpolh(i)=cfl_coef*dls(i)/(ssp+vv)
1386 ssppolh(i)=ssp
1387 ENDIF
1388 ENDDO
1389!$OMP END DO
1390
1391 CASE(2) ! !/DT/FVMBAG/2
1392 param_l_type = fvdata(ifv)%L_TYPE !characteristic length formulation
1393
1394 IF(param_l_type /= 3)THEN !L_TYPE=0/1 & 2
1395!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1396 DO i=1,npolh
1397 IF(mpolh(i) <= zero.OR.epolh(i) <= zero.OR.gpolh(i) <= one) THEN
1398 dtpolh(i)=ep30
1399 qvisc(i) =zero
1400 ELSE
1401 ! 1.calculating Grid Velocity for each polyhedron
1402 pw(1:3,i)=zero
1403 npolyg = ifvpadr(i+1) - ifvpadr(i)
1404 ntria_tot = zero
1405 !loop over polygons
1406 DO j=ifvpadr(i),ifvpadr(i+1)-1
1407 jj=ifvpolh(j)
1408 ntria = ifvtadr(jj+1) - ifvtadr(jj)
1409 wsurf(1:3)=zero
1410 ntria_tot = ntria_tot + ntria
1411 !loop over triangles
1412 DO k=ifvtadr(jj),ifvtadr(jj+1)-1
1413 kk =ifvpoly(k)
1414 IF(kk>0)THEN
1415 wsurf(1)=wsurf(1)+vgrid(1,kk)
1416 wsurf(2)=wsurf(2)+vgrid(2,kk)
1417 wsurf(3)=wsurf(3)+vgrid(3,kk)
1418 ENDIF
1419 enddo! next K (triangle)
1420 pw(1,i)=pw(1,i)+wsurf(1)
1421 pw(2,i)=pw(2,i)+wsurf(2)
1422 pw(3,i)=pw(3,i)+wsurf(3)
1423 enddo! next J (polygon)
1424 !PW:polyhedron mean value over all triangles
1425 pw(1,i)=pw(1,i)/ntria_tot
1426 pw(2,i)=pw(2,i)/ntria_tot
1427 pw(3,i)=pw(3,i)/ntria_tot
1428 ! 2.artificial viscosity
1429 al=pvolu(i)**third
1430 dd=dm(i)/mpolh(i)
1431 ad=max(zero,dd)
1432 gama=gpolh(i)
1433 ssp=sqrt((gama-one)*gama*epolh(i)/mpolh(i))
1434 qvisc(i)=rpolh(i)*ad*al*(qa*qa*al*ad+qb*ssp)
1435 qx=qb*ssp+al*qa*qa*ad
1436 ssp=qx+sqrt(qx*qx+ssp*ssp)
1437 ! 3.relative velocty (PU:material velocity : PW:grid velocity)
1438 tmp_vec(1)=pu(1,i)-pw(1,i)
1439 tmp_vec(2)=pu(2,i)-pw(2,i)
1440 tmp_vec(3)=pu(3,i)-pw(3,i)
1441 vv=sqrt(tmp_vec(1)*tmp_vec(1)+tmp_vec(2)*tmp_vec(2)+tmp_vec(3)*tmp_vec(3))
1442 ! 4.time step
1443 dtpolh(i)=cfl_coef*dls(i)/(ssp+vv)
1444 ssppolh(i)=ssp
1445 ENDIF
1446 ENDDO
1447!$OMP END DO
1448
1449 ELSE !L_TYPE=3
1450
1451!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1452 DO i=1,npolh
1453 IF(mpolh(i) <= zero.OR.epolh(i) <= zero.OR.gpolh(i) <= one) THEN
1454 dtpolh(i)=ep30
1455 qvisc(i) =zero
1456 ELSE
1457 ! 1.calculating Grid Velocity for each polyhedron
1458 pw(1:3,i)=zero
1459 denom_max = zero
1460 DO j=ifvpadr(i),ifvpadr(i+1)-1
1461 jj=ifvpolh(j)
1462 ntria = ifvtadr(jj+1) - ifvtadr(jj)
1463 wsurf(1:3)=zero
1464 DO k=ifvtadr(jj),ifvtadr(jj+1)-1
1465 kk =ifvpoly(k)
1466 area=parea(kk)
1467 iel1=ifvtri(5, kk)
1468 iel2=ifvtri(6, kk)
1469 IF(iel2==0)iel2=i !boundary surface
1470 IF(iel1==0)iel1=i
1471 vmat(1) = half*(pu(1,iel1)+pu(1,iel2))
1472 vmat(2) = half*(pu(2,iel1)+pu(2,iel2))
1473 vmat(3) = half*(pu(3,iel1)+pu(3,iel2))
1474 gama=gpolh(i)
1475 ssp1= sqrt((gama-one)*gama*epolh(iel1)/mpolh(iel1))
1476 ssp2= sqrt((gama-one)*gama*epolh(iel2)/mpolh(iel2))
1477 sspt = half*(ssp1+ssp2)
1478 tmp_vec(1) = (vmat(1)-vgrid(1,kk))*pnorm(1,kk)
1479 tmp_vec(2) = (vmat(2)-vgrid(2,kk))*pnorm(2,kk)
1480 tmp_vec(3) = (vmat(3)-vgrid(3,kk))*pnorm(3,kk)
1481 denom = area*(sqrt(tmp_vec(1)*tmp_vec(1)+tmp_vec(2)*tmp_vec(2)+tmp_vec(3)*tmp_vec(3)) + sspt)
1482 denom_max = max(denom, denom_max)
1483 enddo! next K (triangle)
1484 enddo! next J (polygon)
1485 ! 2.artificial viscosity
1486 al=pvolu(i)**third
1487 dd=dm(i)/mpolh(i)
1488 ad=max(zero,dd)
1489 gama=gpolh(i)
1490 ssp=sqrt((gama-one)*gama*epolh(i)/mpolh(i))
1491 qvisc(i)=rpolh(i)*ad*al*(qa*qa*al*ad+qb*ssp)
1492 qx=qb*ssp+al*qa*qa*ad
1493 ssp=qx+sqrt(qx*qx+ssp*ssp)
1494 ! 3.time step
1495 dtpolh(i)=cfl_coef*dls(i)/(denom_max)
1496 ssppolh(i)=ssp
1497 ENDIF
1498 ENDDO
1499!$OMP END DO
1500
1501 ENDIF
1502
1503
1504 END SELECT
1505
1506
1507!$OMP END PARALLEL
1508
1509 ! FVMBAG minimal time step
1510 DO i=1,npolh
1511 IF (dtpolh(i) < dtx) THEN
1512 dtx=dtpolh(i)
1513 phdt=i
1514 ENDIF
1515 ENDDO
1516
1517 !optional time step growth limitor (/DT/FVMBAG/2, lambda parameter)
1518 IF(ncycle > 0 .AND. idt_fvmbag == 2)THEN
1519 param_lambda = fvdata(ifv)%LAMBDA
1520 IF(param_lambda > zero)THEN
1521 dtold = fvdata(ifv)%DTOLD
1522 IF(dtold > zero)THEN
1523 dtx = min( dtx, dtold*(one+param_lambda))
1524 ENDIF
1525 ENDIF
1526 ENDIF
1527 fvdata(ifv)%DTOLD = dtx
1528
1529 IF (ilvout >= 1.AND.ipri == 0) THEN
1530 WRITE(iout,'(A,I10,A,E12.4,A,I10)')' ** MONITORED VOLUME ID: ',ivolu(1),' TIME STEP:',dtx,' FINITE VOLUME:',idpolh(phdt)
1531 ENDIF
1532
1533 IF (dtx < dt2) THEN
1534 dt2=dtx
1535 nelts=ivolu(1)
1536 itypts=52
1537 ENDIF
1538C---------------------------------------------------
1539C Ibntegration of pressure gradient and surfacic forces
1540C---------------------------------------------------
1541!$OMP PARALLEL PRIVATE(N1,N2,N3,IEL,NX,NY,NZ)
1542!$OMP+ PRIVATE(VX1,VX2,VX3,VY1,VY2,VY3,VZ1,VZ2,VZ3)
1543!$OMP+ PRIVATE(VVX,VVY,VVZ)
1544!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1545 DO iel=1,nelt
1546 n1=elem(1,iel)
1547 n2=elem(2,iel)
1548 n3=elem(3,iel)
1549 vx1=vvv(1,n1)
1550 vy1=vvv(2,n1)
1551 vz1=vvv(3,n1)
1552 vx2=vvv(1,n2)
1553 vy2=vvv(2,n2)
1554 vz2=vvv(3,n2)
1555 vx3=vvv(1,n3)
1556 vy3=vvv(2,n3)
1557 vz3=vvv(3,n3)
1558 vvx=third*(vx1+vx2+vx3)
1559 vvy=third*(vy1+vy2+vy3)
1560 vvz=third*(vz1+vz2+vz3)
1561 nx=norm(1,iel)
1562 ny=norm(2,iel)
1563 nz=norm(3,iel)
1564 vel(iel)=nx*vvx+ny*vvy+nz*vvz
1565 ENDDO
1566!$OMP END DO
1567!$OMP END PARALLEL
1568 DO i=1,npolh
1569 IF (mpolh(i) <= zero.OR.epolh(i) <= zero) cycle
1570 DO j=ifvpadr(i),ifvpadr(i+1)-1
1571 jj=ifvpolh(j)
1572 DO k=ifvtadr(jj),ifvtadr(jj+1)-1
1573 kk=ifvpoly(k)
1574 iel=ifvtri(4,kk)
1575 area=parea(kk)
1576 IF (area == zero) cycle
1577 IF (iel > 0) THEN
1578 nx=norm(1,iel)
1579 ny=norm(2,iel)
1580 nz=norm(3,iel)
1581 IF (itagel(iel) < 0) THEN
1582C Vent hole
1583 ivent =-itagel(iel)
1584 idef =ibaghol(1,ivent)
1585 iventyp=ibaghol(13,ivent)
1586C
1587 IF(idef==1.AND.iventyp==0) THEN
1588 aout=elsout(iel)*area/elarea(iel)
1589 area=area-aout
1590 pp=pext
1591 fel(iel)=fel(iel)+pp*aout
1592 dq(1,i)=dq(1,i)-pp*aout*nx
1593 dq(2,i)=dq(2,i)-pp*aout*ny
1594 dq(3,i)=dq(3,i)-pp*aout*nz
1595 de(i)=de(i)-pp*vel(iel)*aout
1596 ENDIF
1597 ENDIF
1598C Triangle on the airbag envelope
1599 pp=ppolh(i)+qvisc(i)
1600 fel(iel)=fel(iel)+pp*area
1601 dq(1,i)=dq(1,i)-pp*area*nx
1602 dq(2,i)=dq(2,i)-pp*area*ny
1603 dq(3,i)=dq(3,i)-pp*area*nz
1604C Work of pressure forces in the airbag displacement
1605 de(i)=de(i)-pp*vel(iel)*area
1606C Heat Loss or Gain (version 2018)
1607 de(i)=de(i)-rvolu(19)*area*(tpolh(i)-rvolu(25))
1608 ELSEIF(iel == 0) THEN
1609C triangle permeability
1610 area_old=area
1611 area=area*max(zero,one-rvolu(51)/hmin(kk))
1612 area1=area_old-area
1613C internal triangle = communication between polyhedra
1614 nv = 0
1615 IF (ifvtri(5,kk) == i) THEN
1616 nv=ifvtri(6,kk)
1617 nx=pnorm(1,kk)
1618 ny=pnorm(2,kk)
1619 nz=pnorm(3,kk)
1620 ELSEIF (ifvtri(6,kk) == i) THEN
1621 nv=ifvtri(5,kk)
1622 nx=-pnorm(1,kk)
1623 ny=-pnorm(2,kk)
1624 nz=-pnorm(3,kk)
1625 ENDIF
1626 IF (nv < i .OR.mpolh(nv) <= zero.OR.
1627 . epolh(nv) <= zero) GOTO 200
1628 p1=ppolh(i)+qvisc(i)
1629 p2=ppolh(nv)+qvisc(nv)
1630 pmean=half*(p1+p2)
1631C
1632 dq(1,i)=dq(1,i)-pmean*nx*area-p1*nx*area1
1633 dq(2,i)=dq(2,i)-pmean*ny*area-p1*ny*area1
1634 dq(3,i)=dq(3,i)-pmean*nz*area-p1*nz*area1
1635 dq(1,nv)=dq(1,nv)+pmean*nx*area+p2*nx*area1
1636 dq(2,nv)=dq(2,nv)+pmean*ny*area+p2*ny*area1
1637 dq(3,nv)=dq(3,nv)+pmean*nz*area+p2*nz*area1
1638C added FA
1639 ss=nx*vgrid(1,kk)+ny*vgrid(2,kk)+nz*vgrid(3,kk)
1640 de(i) =de(i) -pmean*ss*area-p1*ss*area1
1641 de(nv)=de(nv)+pmean*ss*area+p2*ss*area1
1642 ELSEIF(iel < 0) THEN
1643C internal triangle on internal surface = communication between polyhedra
1644 nv = 0
1645 IF (ifvtri(5,kk) == i) THEN
1646 nv=ifvtri(6,kk)
1647 nx=pnorm(1,kk)
1648 ny=pnorm(2,kk)
1649 nz=pnorm(3,kk)
1650 ELSEIF (ifvtri(6,kk) == i) THEN
1651 nv=ifvtri(5,kk)
1652 nx=-pnorm(1,kk)
1653 ny=-pnorm(2,kk)
1654 nz=-pnorm(3,kk)
1655 ENDIF
1656 IF (nv <= i .OR.mpolh(nv) <= zero.OR.
1657 . epolh(nv) <= zero) GOTO 200
1658 p1=ppolh(i)+qvisc(i)
1659 p2=ppolh(nv)+qvisc(nv)
1660 jel=-iel
1661 pmean=half*(p1+p2)
1662 area_old=area
1663 area=area*porosity(jel-nel)
1664 area1=area_old-area
1665 ss=nx*vgrid(1,kk)+ny*vgrid(2,kk)+nz*vgrid(3,kk)
1666C
1667 nrx=norm(1,jel)
1668 nry=norm(2,jel)
1669 nrz=norm(3,jel)
1670 fac=nx*nrx+ny*nry+nz*nrz
1671C
1672 IF(itagel(jel) > 0) THEN
1673 IF(fac < 0) THEN
1674 fel(jel)=fel(jel)+(p2-p1)*area1
1675 dq(1,i)=dq(1,i)-pmean*nx*area+p1*nrx*area1
1676 dq(2,i)=dq(2,i)-pmean*ny*area+p1*nry*area1
1677 dq(3,i)=dq(3,i)-pmean*nz*area+p1*nrz*area1
1678 de(i)=de(i)-pmean*ss*area+p1*vel(jel)*area1
1679 dq(1,nv)=dq(1,nv)+pmean*nx*area-p2*nrx*area1
1680 dq(2,nv)=dq(2,nv)+pmean*ny*area-p2*nry*area1
1681 dq(3,nv)=dq(3,nv)+pmean*nz*area-p2*nrz*area1
1682 de(nv)=de(nv)+pmean*ss*area-p2*vel(jel)*area1
1683 ELSE
1684 fel(jel)=fel(jel)+(p1-p2)*area1
1685 dq(1,i)=dq(1,i)-pmean*nx*area-p1*nrx*area1
1686 dq(2,i)=dq(2,i)-pmean*ny*area-p1*nry*area1
1687 dq(3,i)=dq(3,i)-pmean*nz*area-p1*nrz*area1
1688 de(i)=de(i)-pmean*ss*area-p1*vel(jel)*area1
1689 dq(1,nv)=dq(1,nv)+pmean*nx*area+p2*nrx*area1
1690 dq(2,nv)=dq(2,nv)+pmean*ny*area+p2*nry*area1
1691 dq(3,nv)=dq(3,nv)+pmean*nz*area+p2*nrz*area1
1692 de(nv)=de(nv)+pmean*ss*area+p2*vel(jel)*area1
1693 ENDIF
1694 ELSE
1695 IF(fac < 0) THEN
1696 fel(jel)=fel(jel)+(p2-p1)*area1
1697 ELSE
1698 fel(jel)=fel(jel)+(p1-p2)*area1
1699 ENDIF
1700 dq(1,i)=dq(1,i)-pmean*nx*area-p1*nx*area1
1701 dq(2,i)=dq(2,i)-pmean*ny*area-p1*ny*area1
1702 dq(3,i)=dq(3,i)-pmean*nz*area-p1*nz*area1
1703 de(i)=de(i)-pmean*ss*area-p1*ss*area1
1704 dq(1,nv)=dq(1,nv)+pmean*nx*area+p2*nx*area1
1705 dq(2,nv)=dq(2,nv)+pmean*ny*area+p2*ny*area1
1706 dq(3,nv)=dq(3,nv)+pmean*nz*area+p2*nz*area1
1707 de(nv)=de(nv)+pmean*ss*area+p2*ss*area1
1708 ENDIF
1709 ENDIF
1710 ENDDO
1711 200 ENDDO
1712 ENDDO
1713
1714C---------------------------------------------------
1715C New Polyhedra quantities
1716C---------------------------------------------------
1717!$OMP PARALLEL PRIVATE(I)
1718!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1719 DO i=1,npolh
1720 qpolh(1,i)=qpolh(1,i)+dq(1,i)*dt1
1721 qpolh(2,i)=qpolh(2,i)+dq(2,i)*dt1
1722 qpolh(3,i)=qpolh(3,i)+dq(3,i)*dt1
1723 epolh(i)=epolh(i)+de(i)*dt1
1724 ENDDO
1725!$OMP END DO
1726!$OMP END PARALLEL
1727
1728C---------------------------------------------------
1729C Impressions
1730C---------------------------------------------------
1731 250 IF (ilvout >= 1 .AND. ipri == 0) THEN
1732 WRITE(iout,*)
1733 WRITE(iout,'(A25,I10,A31)')' ** MONITORED VOLUME ID: ',ivolu(1), ' - FINITE VOLUME MESH UPDATE **'
1734 WRITE(iout,'(A,I8)')' NUMBER OF FINITE VOLUMES : ',npolh
1735 WRITE(iout,'(A,G16.9,A17,G16.9,A1)')' SUM VOLUME FINITE VOLUMES : ',volph,' (VOLUME AIRBAG: ',volg,')'
1736 WRITE(iout,'(A,G16.9,A17,G16.9,A1)')' SUM AREA SURFACE POLYGONS : ',areap,' (AREA AIRBAG : ',areael,')'
1737 WRITE(iout,*)
1738 ENDIF
1739C---------------------------------------------------
1740C TH output
1741C---------------------------------------------------
1742 amtot=zero
1743 pmean=zero
1744 enint=zero
1745 fsav(1:17)=zero
1746 DO i=1,npolh
1747 amtot=amtot+mpolh(i)
1748 enint=enint+epolh(i)
1749 pmean=pmean+ppolh(i)*pvolu(i)
1750 gama=gpolh(i)
1751 rmw=rmwpolh(i)
1752 temp=tpolh(i)
1753 cpa=cpapolh(i)
1754 cpb=cpbpolh(i)
1755 cpc=cpcpolh(i)
1756 cp=cpa+cpb*temp+cpc*temp*temp
1757 IF(ityp == 8) THEN
1758 cpd=cpdpolh(i)
1759 cpe=cpepolh(i)
1760 cpf=cpfpolh(i)
1761 cp=cp+cpd*temp**3+cpf*temp**4
1762 temp2 = temp * temp
1763 IF(temp2 > zero) cp=cp+cpe/(temp*temp)
1764 ENDIF
1765 cv=double_to_my_real(cp)-rmw
1766 fsav(5) =fsav(5) +mpolh(i)*temp
1767 fsav(10)=fsav(10)+mpolh(i)*double_to_my_real(cp)
1768 fsav(11)=fsav(11)+mpolh(i)*cv
1769 fsav(12)=fsav(12)+mpolh(i)*gama
1770 ENDDO
1771C Pressure dispersion
1772 pdisp = zero
1773 pmean2 = pmean / volph
1774 DO i = 1, npolh
1775 pdisp = pdisp + pvolu(i) * (ppolh(i) - pmean2)**2
1776 ENDDO
1777 pdisp = sqrt(pdisp / volph) / pmean2
1778 fsav(19) = pdisp
1779 injection_started = .false.
1780 DO iinj = 1, njet
1781 injection_started = injection_started .OR. (datainj(2, iinj) > zero)
1782 ENDDO
1783 IF (.NOT. injection_started) THEN
1784 pdisp = zero
1785 ENDIF
1786
1787C
1788 fsav(1)=amtot
1789 fsav(2)=volph
1790 pres_av = pmean / volph
1791 fsav(3) = pres_av
1792 fsav(4)=areap
1793 fsav(7)=zero
1794 dmout =zero
1795 dhout =zero
1796 cp_av = fsav(10) / amtot
1797 fsav(10) = cp_av
1798 cv_av = fsav(11) / amtot
1799 fsav(11) = cv_av
1800 fsav(12) = cp_av / cv_av
1801 fsav(14) = npolh
1802 rgas_av = cp_av - cv_av
1803 rho_av = amtot / volph
1804 fsav(5) = fsav(5) / amtot ! sum(mass[i]*T[i])/total_mass
1805
1806 IF(ivolu(39) == 0) THEN
1807 fsav(6) =zero
1808 fsav(13)=zero
1809 ELSE
1810 fsav(6) =aoutot
1811 DO iel=1,nel
1812 IF (itagel(iel) < 0) THEN
1813 ivent=-itagel(iel)
1814 rbaghol(18,ivent)=rbaghol(18,ivent)+elfvel(iel)*elsout(iel)
1815 rbaghol(19,ivent)=rbaghol(19,ivent)-elfmass(iel)
1816 rbaghol(20,ivent)=rbaghol(20,ivent)-elfehpy(iel)
1817 rbaghol(21,ivent)=rbaghol(21,ivent)+eldmass(iel)
1818 rbaghol(22,ivent)=rbaghol(22,ivent)+eldehpy(iel)
1819 ENDIF
1820 ENDDO
1821 DO ivent=1,nvent
1822 sventot=rbaghol(16,ivent)+rbaghol(17,ivent)
1823 IF(sventot <= zero) cycle
1824 rbaghol(18,ivent)=rbaghol(18,ivent)/sventot
1825 ENDDO
1826 DO ivent=1,nvent
1827 fsav(7)=fsav(7)+(rbaghol(16,ivent)+rbaghol(17,ivent))*rbaghol(18,ivent)
1828 dmout =dmout+rbaghol(21,ivent)
1829 dhout =dhout+rbaghol(22,ivent)
1830 ENDDO
1831 fsav(7) =fsav(7) /max(aoutot,em20)
1832 fsav(13)=dtx ! not DTPOLH(PHDT) since DTX may be updated with smoothing factor
1833
1834 IF(ityp == 8) THEN
1835 CALL fvinjt8_1(njet, ibagjet , rbagjet , igeo, geo, pm,
1836 . ivolu, rvolu, dmout, dhout)
1837 ttf =rvolu(60)
1838 IF (ivolu(74) == 0) THEN
1839 up_switch = tt-ttf >= rvolu(70) .OR. npolh <= ivolu(37)
1840 ENDIF
1841 IF (ivolu(74) == 1) THEN
1842C Automatic switch to uniform pressure when dispersion of pressure is low
1843 up_switch = ((pdisp < pdisp_old) .AND. (pdisp < rvolu(73)))
1844 up_switch = up_switch .OR. tt-ttf >= rvolu(70)
1845 ENDIF
1846 IF (ivolu(74) == 2) THEN
1847C Automatic switch to uniform pressure when dispersion of pressure is low
1848 up_switch = ((pdisp < pdisp_old) .AND. (pdisp < rvolu(73)))
1849 up_switch = up_switch .AND. tt-ttf >= rvolu(70)
1850 ENDIF
1851
1852 IF (up_switch) THEN
1853 rvolu(12)=fsav(3) ! pression
1854 rvolu(13)=fsav(5) ! temperature
1855 rvolu(16)=fsav(2) ! volume
1856 rvolu(20)=fsav(1) ! masse
1857 ENDIF
1858 ENDIF
1859
1860 ENDIF
1861
1862
1863
1864 DO iel=1,nelt
1865 IF (itagel(iel) > 0) THEN
1866 fsav(15)=fsav(15)+elfmass(iel)
1867 fsav(16)=fsav(16)+elfehpy(iel)
1868 ENDIF
1869 ENDDO
1870 fsav(17)=enint
1871
1872
1873C---------------------------------------------------
1874C Fluid velocity on nodes
1875C---------------------------------------------------
1876 IF( (tt>=tanim .AND. tt<=tanim_stop) .OR.tt >= toutp.OR. nbgauge > 0 .OR.
1877 . (manim >= 4.AND.manim <= 15)) THEN
1878C Computations for anim, also needed if there are gauges.
1879 prepare_anim = 1
1880 uel(1:3,1:nelt)=zero
1881 rhoel(1:nelt)=zero
1882 tkel(1:nelt) =zero
1883 sspel(1:nelt)=zero
1884 ELSE
1885 prepare_anim = 0
1886 ENDIF
1887
1888!$OMP PARALLEL PRIVATE(I,J,JJ,K,KK,IEL,JEL,II,NV)
1889!$OMP+ PRIVATE(NX,NY,NZ,NRX,NRY,NRZ)
1890!$OMP+ PRIVATE(FAC,AREA)
1891!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1892 DO i=1,npolh
1893 DO j=ifvpadr(i),ifvpadr(i+1)-1
1894 jj=ifvpolh(j)
1895 DO k=ifvtadr(jj),ifvtadr(jj+1)-1
1896 kk=ifvpoly(k)
1897 iel=abs(ifvtri(4,kk))
1898 jel=ifvtri(4,kk)
1899 area=parea(kk)
1900 IF (jel > 0) THEN
1901 pel(iel) =pel(iel) +ppolh(i)*area/elarea(iel)
1902 ELSEIF (jel < 0) THEN
1903 nrx=norm(1,iel)
1904 nry=norm(2,iel)
1905 nrz=norm(3,iel)
1906 IF (ifvtri(5,kk) == i) THEN
1907 nv=ifvtri(6,kk)
1908 nx=pnorm(1,kk)
1909 ny=pnorm(2,kk)
1910 nz=pnorm(3,kk)
1911 ELSEIF(ifvtri(6,kk) == i) THEN
1912 nv=ifvtri(5,kk)
1913 nx=-pnorm(1,kk)
1914 ny=-pnorm(2,kk)
1915 nz=-pnorm(3,kk)
1916 ENDIF
1917 IF (nv >= i) THEN
1918 IF (nv == i) area=area*half
1919 ii=i
1920 fac=nx*nrx+ny*nry+nz*nrz
1921 IF(fac < 0) ii=nv
1922 pel(iel) =pel(iel) +ppolh(ii)*area/elarea(iel)
1923 ENDIF
1924 ENDIF
1925 ENDDO ! next triangle
1926 ENDDO ! next polygon
1927 ENDDO ! next polyhedron
1928!$OMP END DO
1929
1930 IF(prepare_anim == 1) THEN
1931!$OMP SINGLE
1932 DO i=1,npolh
1933 DO j=ifvpadr(i),ifvpadr(i+1)-1
1934 jj=ifvpolh(j)
1935 DO k=ifvtadr(jj),ifvtadr(jj+1)-1
1936 kk=ifvpoly(k)
1937 iel=abs(ifvtri(4,kk))
1938 jel=ifvtri(4,kk)
1939 area=parea(kk)
1940 IF (jel > 0) THEN
1941 uel(1,iel)=uel(1,iel)+pu(1,i)*area/elarea(iel)
1942 uel(2,iel)=uel(2,iel)+pu(2,i)*area/elarea(iel)
1943 uel(3,iel)=uel(3,iel)+pu(3,i)*area/elarea(iel)
1944 rhoel(iel)=rhoel(iel)+rpolh(i)*area/elarea(iel)
1945 tkel(iel) =tkel(iel) +tpolh(i)*area/elarea(iel)
1946 sspel(iel)=sspel(iel)+sqrt(gpolh(i)*ppolh(i)/rpolh(i))*area/elarea(iel)
1947 ELSEIF (jel < 0) THEN
1948 nrx=norm(1,iel)
1949 nry=norm(2,iel)
1950 nrz=norm(3,iel)
1951 IF (ifvtri(5,kk) == i) THEN
1952 nv=ifvtri(6,kk)
1953 nx=pnorm(1,kk)
1954 ny=pnorm(2,kk)
1955 nz=pnorm(3,kk)
1956 ELSEIF(ifvtri(6,kk) == i) THEN
1957 nv=ifvtri(5,kk)
1958 nx=-pnorm(1,kk)
1959 ny=-pnorm(2,kk)
1960 nz=-pnorm(3,kk)
1961 ENDIF
1962 IF (nv < i) GOTO 310
1963 IF (nv == i) area=area*half
1964 ii=i
1965 fac=nx*nrx+ny*nry+nz*nrz
1966 IF(fac < 0) ii=nv
1967 uel(1,iel)=uel(1,iel)+pu(1,ii)*area/elarea(iel)
1968 uel(2,iel)=uel(2,iel)+pu(2,ii)*area/elarea(iel)
1969 uel(3,iel)=uel(3,iel)+pu(3,ii)*area/elarea(iel)
1970 rhoel(iel)=rhoel(iel)+rpolh(ii)*area/elarea(iel)
1971 tkel(iel) =tkel(iel) +tpolh(ii)*area/elarea(iel)
1972 sspel(iel)=sspel(iel)+sqrt(gpolh(ii)*ppolh(ii)/rpolh(ii))*area/elarea(iel)
1973 ENDIF
1974 ENDDO ! next triangle
1975 310 ENDDO ! next polygon
1976 ENDDO ! next polyhedron
1977!$OMP END SINGLE
1978 ENDIF
1979
1980
1981
1982C------------------
1983C VENT HOLES
1984C------------------
1985 IF(prepare_anim == 1) THEN
1986!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1987 DO iel=1,nel
1988 IF(itagel(iel)<0) THEN
1989 uel(1,iel)=elfvel(iel)*norm(1,iel)
1990 uel(2,iel)=elfvel(iel)*norm(2,iel)
1991 uel(3,iel)=elfvel(iel)*norm(3,iel)
1992 ENDIF
1993 ENDDO
1994!$OMP END DO NOWAIT
1995
1996!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1997 DO i=1,nnt
1998 u(1,i)=zero
1999 u(2,i)=zero
2000 u(3,i)=zero
2001 rho(i) =zero
2002 tk(i) =zero
2003 sspk(i) =zero
2004 ENDDO
2005!OMP END DO
2006
2007!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
2008 DO i=1,nna
2009 vola(i)=zero
2010 pa(i)=zero
2011 rhoa(i)=zero
2012 tka(i)=zero
2013 sspka(i)=zero
2014 ua(1,i)=zero
2015 ua(2,i)=zero
2016 ua(3,i)=zero
2017 ENDDO
2018!$OMP ENDDO
2019
2020 ENDIF
2021
2022!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
2023 DO i=1,nnt
2024 p(i)=zero
2025 ENDDO
2026!$OMP END DO
2027!$OMP END PARALLEL
2028
2029 DO iel=1,nelt
2030 n1=elem(1,iel)
2031 n2=elem(2,iel)
2032 n3=elem(3,iel)
2033 area1=nodarea(n1)
2034 area2=nodarea(n2)
2035 area3=nodarea(n3)
2036 fac1=third*elarea(iel)/area1
2037 fac2=third*elarea(iel)/area2
2038 fac3=third*elarea(iel)/area3
2039 p(n1) =p(n1) +fac1*pel(iel)
2040 p(n2) =p(n2) +fac2*pel(iel)
2041 p(n3) =p(n3) +fac3*pel(iel)
2042 ENDDO
2043
2044 IF(prepare_anim == 1) THEN
2045 DO iel=1,nelt
2046 n1=elem(1,iel)
2047 n2=elem(2,iel)
2048 n3=elem(3,iel)
2049 area1=nodarea(n1)
2050 area2=nodarea(n2)
2051 area3=nodarea(n3)
2052 fac1=third*elarea(iel)/area1
2053 fac2=third*elarea(iel)/area2
2054 fac3=third*elarea(iel)/area3
2055 u(1,n1)=u(1,n1)+fac1*uel(1,iel)
2056 u(2,n1)=u(2,n1)+fac1*uel(2,iel)
2057 u(3,n1)=u(3,n1)+fac1*uel(3,iel)
2058 u(1,n2)=u(1,n2)+fac2*uel(1,iel)
2059 u(2,n2)=u(2,n2)+fac2*uel(2,iel)
2060 u(3,n2)=u(3,n2)+fac2*uel(3,iel)
2061 u(1,n3)=u(1,n3)+fac3*uel(1,iel)
2062 u(2,n3)=u(2,n3)+fac3*uel(2,iel)
2063 u(3,n3)=u(3,n3)+fac3*uel(3,iel)
2064 rho(n1)=rho(n1)+fac1*rhoel(iel)
2065 rho(n2)=rho(n2)+fac2*rhoel(iel)
2066 rho(n3)=rho(n3)+fac3*rhoel(iel)
2067 tk(n1)=tk(n1)+fac1*tkel(iel)
2068 tk(n2)=tk(n2)+fac2*tkel(iel)
2069 tk(n3)=tk(n3)+fac3*tkel(iel)
2070 sspk(n1)=sspk(n1)+fac1*sspel(iel)
2071 sspk(n2)=sspk(n2)+fac2*sspel(iel)
2072 sspk(n3)=sspk(n3)+fac3*sspel(iel)
2073 ENDDO
2074
2075 DO i=1,npolh
2076 IF (ibpolh(i) /= 0) THEN
2077 temp=tpolh(i)
2078 ii=iabs(ibpolh(i))
2079 IF(brna(1,ii) == brna(4,ii).AND.brna(5,ii) == brna(8,ii))THEN
2080C Pentahedron
2081 DO k=1,6
2082 j=k
2083 IF(k >= 4) j=j+1
2084 jj=brna(j,ii)
2085 vola(jj)=vola(jj)+one_over_6*pvolu(i)
2086 pa(jj)=pa(jj)+one_over_6*ppolh(i)*pvolu(i)
2087 rhoa(jj)=rhoa(jj)+one_over_6*rpolh(i)*pvolu(i)
2088 tka(jj)=tka(jj)+one_over_6*temp*pvolu(i)
2089 sspka(jj)=sspka(jj)+one_over_6*sqrt(gpolh(i)*ppolh(i)/rpolh(i))*pvolu(i)
2090 ua(1,jj)=ua(1,jj)+one_over_6*pu(1,i)*pvolu(i)
2091 ua(2,jj)=ua(2,jj)+one_over_6*pu(2,i)*pvolu(i)
2092 ua(3,jj)=ua(3,jj)+one_over_6*pu(3,i)*pvolu(i)
2093 ENDDO
2094 ELSEIF(brna(5,ii)==brna(8,ii).AND.
2095 . brna(6,ii)==brna(8,ii).AND.
2096 . brna(7,ii)==brna(8,ii))THEN
2097C Pyramid
2098 DO j=1,5
2099 jj=brna(j,ii)
2100 vola(jj)=vola(jj)+one_fifth*pvolu(i)
2101 pa(jj)=pa(jj)+one_fifth*ppolh(i)*pvolu(i)
2102 rhoa(jj)=rhoa(jj)+one_fifth*rpolh(i)*pvolu(i)
2103 tka(jj)=tka(jj)+one_fifth*temp*pvolu(i)
2104 sspka(jj)=sspka(jj)+one_fifth*sqrt(gpolh(i)*ppolh(i)/rpolh(i))*pvolu(i)
2105 ua(1,jj)=ua(1,jj)+one_fifth*pu(1,i)*pvolu(i)
2106 ua(2,jj)=ua(2,jj)+one_fifth*pu(2,i)*pvolu(i)
2107 ua(3,jj)=ua(3,jj)+one_fifth*pu(3,i)*pvolu(i)
2108 ENDDO
2109 ELSE
2110C Hexahedron or Tetrahedron
2111 DO j=1,8
2112 jj=brna(j,ii)
2113 vola(jj)=vola(jj)+one_over_8*pvolu(i)
2114 pa(jj)=pa(jj)+one_over_8*ppolh(i)*pvolu(i)
2115 rhoa(jj)=rhoa(jj)+one_over_8*rpolh(i)*pvolu(i)
2116 tka(jj)=tka(jj)+one_over_8*temp*pvolu(i)
2117 sspka(jj)=sspka(jj)+one_over_8*sqrt(gpolh(i)*ppolh(i)/rpolh(i))*pvolu(i)
2118 ua(1,jj)=ua(1,jj)+one_over_8*pu(1,i)*pvolu(i)
2119 ua(2,jj)=ua(2,jj)+one_over_8*pu(2,i)*pvolu(i)
2120 ua(3,jj)=ua(3,jj)+one_over_8*pu(3,i)*pvolu(i)
2121 ENDDO
2122 ENDIF
2123 ENDIF
2124 ENDDO
2125
2126 DO i=1,nna
2127 volno=vola(i)
2128 IF(volno > zero) THEN
2129 pa(i) =pa(i)/volno
2130 rhoa(i)=rhoa(i)/volno
2131 tka(i) =tka(i)/volno
2132 sspka(i) =sspka(i)/volno
2133 ua(1,i)=ua(1,i)/volno
2134 ua(2,i)=ua(2,i)/volno
2135 ua(3,i)=ua(3,i)/volno
2136 ENDIF
2137 ENDDO
2138 ENDIF ! prepare_anim
2139
2140
2141C---------------------------------------------------
2142C Potential opening of vent holes
2143C---------------------------------------------------
2144 DO ivent=1,nvent
2145 area_vent(ivent)=zero
2146 pm_vent(ivent)=zero
2147 ENDDO
2148
2149C
2150 DO iel=1,nel
2151 IF (itagel(iel) < 0) THEN
2152 ivent=-itagel(iel)
2153 area_vent(ivent)=area_vent(ivent)+elarea(iel)
2154 pm_vent(ivent)=pm_vent(ivent)+fel(iel)
2155 ENDIF
2156 ENDDO
2157
2158
2159
2160 DO ivent=1,nvent
2161 IF (area_vent(ivent) > zero) THEN
2162 pm_vent(ivent)=pm_vent(ivent)/area_vent(ivent)
2163 ELSE
2164 pm_vent(ivent)=zero
2165 ENDIF
2166C
2167 idef =ibaghol(1,ivent)
2168 idtpdef= ibaghol(11,ivent)
2169 idpdef = ibaghol(12,ivent)
2170C
2171 IF(idef == 2) cycle
2172C
2173 pdef =rbaghol(1,ivent)
2174 dtpdefi=rbaghol(4,ivent)
2175 dtpdefc=rbaghol(5,ivent)
2176 tvent =rbaghol(3,ivent)
2177 tstope =rbaghol(14,ivent)
2178 IF (idtpdef == 0) THEN
2179 IF(pm_vent(ivent) > pdef+pext) THEN
2180 dtpdefc = dtpdefc+dt1
2181 END IF
2182 ELSE IF (idtpdef == 1) THEN
2183 IF (pm_vent(ivent) > pdef+pext) THEN
2184 idpdef = 1
2185 END IF
2186 IF (idpdef == 1) THEN
2187 dtpdefc = dtpdefc+dt1
2188 END IF
2189 END IF
2190 ibaghol(12,ivent) = idpdef
2191 rbaghol(5,ivent) = dtpdefc
2192C
2193 ittf = ivolu(17)
2194 ttf =rvolu(60)
2195 DO k=1,20
2196 titrevent(k)=ibaghol(14+k,ivent)
2197 venttitle(k:k) = achar(titrevent(k))
2198 ENDDO
2199 IF(ittf==11.OR.ittf==12.OR.ittf==13) THEN
2200 IF (idef==0.AND.pm_vent(ivent) > pdef+pext
2201 . .AND.dtpdefc > dtpdefi
2202 . .AND.volph > em3*areap**three_half
2203 . .AND.tt < tstope+ttf) THEN
2204 idef=1
2205 WRITE(iout,'(A)')' ** AIRBAG VENT HOLE MEMBRANE IS DEFLATED **'
2206 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),' VENT HOLE NUMBER',ivent,' ',venttitle
2207 WRITE(istdo,'(2A)')' ** VENT HOLE MEMBRANE IS DEFLATED ',venttitle
2208 ENDIF
2209 IF(idef == 0 .AND. tt > tvent+ttf.AND. tt < tstope+ttf) THEN
2210 idef=1
2211 WRITE(iout,'(A)') ' ** AIRBAG VENTING STARTS **'
2212 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),' VENT HOLE NUMBER',ivent,' ',venttitle
2213 WRITE(istdo,'(2A)') ' ** VENTING STARTS ',venttitle
2214 ENDIF
2215 IF(idef == 1 .AND. tt >= tstope+ttf) THEN
2216 idef=2
2217 WRITE(iout,'(A)') ' ** AIRBAG VENTING STOPS **'
2218 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),' VENT HOLE NUMBER',ivent,' ',venttitle
2219 WRITE(istdo,'(2A)') ' ** VENTING STOPS ',venttitle
2220 ENDIF
2221
2222 ELSE IF(ittf==0) THEN
2223 IF (idef==0.AND.pm_vent(ivent)>pdef+pext
2224 . .AND.dtpdefc>dtpdefi
2225 . .AND.volph>em3*areap**three_half
2226 . .AND.tt<tstope) THEN
2227 idef=1
2228 WRITE(iout,'(A)')' ** AIRBAG VENT HOLE MEMBRANE IS DEFLATED **'
2229 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1), ' VENT HOLE NUMBER',ivent,' ',venttitle
2230 WRITE(istdo,'(2A)')' ** VENT HOLE MEMBRANE IS DEFLATED ',venttitle
2231 ENDIF
2232 IF(idef == 0 .AND. tt > tvent .AND. tt < tstope) THEN
2233 idef=1
2234 WRITE(iout,'(A)') ' ** AIRBAG VENTING STARTS **'
2235 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),' VENT HOLE NUMBER',ivent,' ',venttitle
2236 WRITE(istdo,'(2A)') ' ** VENTING STARTS ',venttitle
2237 ENDIF
2238 IF(idef == 1 .AND. tt >= tstope) THEN
2239 idef=2
2240 WRITE(iout,'(A)') ' ** AIRBAG VENTING STOPS **'
2241 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),
2242 . ' VENT HOLE NUMBER',ivent,' ',venttitle
2243 WRITE(istdo,'(2A)') ' ** VENTING STOPS ',venttitle
2244 ENDIF
2245 END IF
2246 ibaghol(1,ivent)=idef
2247 ENDDO
2248C---------------------------------------------------
2249C Force on Airbag + OUTPUT
2250C---------------------------------------------------
2251 IF(ivolu(39)==1) THEN
2252 wfext0=wfext
2253 DO iel=1,nelt
2254 n1=elem(1,iel)
2255 n2=elem(2,iel)
2256 n3=elem(3,iel)
2257 nrx=norm(1,iel)
2258 nry=norm(2,iel)
2259 nrz=norm(3,iel)
2260 IF(iel <= nel) THEN
2261 ff=fel(iel)-pext*elarea(iel)
2262 ELSE
2263 ff=fel(iel)
2264 ENDIF
2265 fvspmd(ifv)%AAA(1,n1)=fvspmd(ifv)%AAA(1,n1)+third*ff*nrx
2266 fvspmd(ifv)%AAA(2,n1)=fvspmd(ifv)%AAA(2,n1)+third*ff*nry
2267 fvspmd(ifv)%AAA(3,n1)=fvspmd(ifv)%AAA(3,n1)+third*ff*nrz
2268 fvspmd(ifv)%AAA(1,n2)=fvspmd(ifv)%AAA(1,n2)+third*ff*nrx
2269 fvspmd(ifv)%AAA(2,n2)=fvspmd(ifv)%AAA(2,n2)+third*ff*nry
2270 fvspmd(ifv)%AAA(3,n2)=fvspmd(ifv)%AAA(3,n2)+third*ff*nrz
2271 fvspmd(ifv)%AAA(1,n3)=fvspmd(ifv)%AAA(1,n3)+third*ff*nrx
2272 fvspmd(ifv)%AAA(2,n3)=fvspmd(ifv)%AAA(2,n3)+third*ff*nry
2273 fvspmd(ifv)%AAA(3,n3)=fvspmd(ifv)%AAA(3,n3)+third*ff*nrz
2274C
2275 wfext=wfext+ff*vel(iel)*dt1
2276 ENDDO
2277 fsav(18)=fsav(18)+wfext-wfext0
2278 ENDIF
2279
2280
2281
2282 !------------------------------------------------------------!
2283 ! ANIM / H3D : LOCATE POLYHEDRA CENTROIDS !
2284 !------------------------------------------------------------!
2285 !whenever output file is required
2286 IF((tt>=tanim .AND. tt<=tanim_stop) .OR. tt>=toutp .OR.
2287 . (tt>=h3d_data%TH3D .AND. tt<=h3d_data%TH3D_STOP) .OR.
2288 . (manim>=4.AND.manim<=15).OR.h3d_data%MH3D/=0)THEN
2289
2290 !check if IH3D_FLAG ENABLES OUTPUT. OTHERWISE CALCULATION OF CENTROID IS USELESS
2291 IF(ivolu(75)==1)THEN
2292
2293 DO i=1,npolh !loop over polyhedra composing the airbag mesh
2294 centroid(1:3)=zero
2295 denom=zero
2296 DO j=ifvpadr(i),ifvpadr(i+1)-1 !loop over polygons composing the polyhedron
2297 jj=ifvpolh(j)
2298 DO k=ifvtadr(jj),ifvtadr(jj+1)-1 !loop over triangles composing the polygon
2299 kk =ifvpoly(k)
2300 IF(kk>0)THEN
2301 ntria=ntria+1
2302 !triangle nodes
2303 n1=ifvtri(1,kk)
2304 n2=ifvtri(2,kk)
2305 n3=ifvtri(3,kk)
2306 !coordinates
2307 x1=pnod(1,n1)
2308 y1=pnod(2,n1)
2309 z1=pnod(3,n1)
2310 x2=pnod(1,n2)
2311 y2=pnod(2,n2)
2312 z2=pnod(3,n2)
2313 x3=pnod(1,n3)
2314 y3=pnod(2,n3)
2315 z3=pnod(3,n3)
2316 !normal vector
2317 nx= (y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)
2318 ny= (z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)
2319 nz= (x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)
2320 area=half*sqrt(nx*nx + ny*ny + nz*nz)
2321 denom=denom+area
2322 nx=half*nx
2323 ny=half*ny
2324 nz=half*nz
2325 !centroid x,y,z
2326 centroid(1) = centroid(1)+ area*( third*(x1+x2+x3) )
2327 centroid(2) = centroid(2)+ area*( third*(y1+y2+y3) )
2328 centroid(3) = centroid(3)+ area*( third*(z1+z2+z3) )
2329 ENDIF
2330 enddo! next K (triangle)
2331 enddo! next J (polygon)
2332 IF(denom>zero)THEN
2333 centroid(1) = centroid(1)/denom
2334 centroid(2) = centroid(2)/denom
2335 centroid(3) = centroid(3)/denom
2336 ENDIF
2337
2338 centroid_polh(1,i) = centroid(1)
2339 centroid_polh(2,i) = centroid(2)
2340 centroid_polh(3,i) = centroid(3)
2341
2342 enddo! next I (next polyhedron)
2343
2344 ENDIF !IH3D_FLAG
2345 endif! H3D/ANIM CRITERION
2346
2347C----------------------
2348C TH Gauges output
2349C----------------------
2350 IF (nbgauge > 0) THEN
2351C To be done even if gauge nodes are not related to FVMBAG.
2352C otherwise we may to miss some in SPMD
2353 DO i=1,nnt
2354 fvspmd(ifv)%GGG(1,i)=p(i)
2355 fvspmd(ifv)%GGG(2,i)=rho(i)
2356 fvspmd(ifv)%GGG(3,i)=tk(i)
2357 ENDDO
2358
2359 DO i=1,nna
2360 fvspmd(ifv)%GGA(1,i)=pa(i)
2361 fvspmd(ifv)%GGA(2,i)=rhoa(i)
2362 fvspmd(ifv)%GGA(3,i)=tka(i)
2363 ENDDO
2364 ENDIF
2365
2366 1000 CONTINUE
2367 IF (ALLOCATED(nodarea)) DEALLOCATE(nodarea)
2368 IF (ALLOCATED(tagsurf)) DEALLOCATE(tagsurf)
2369 IF (ALLOCATED(icont)) DEALLOCATE(icont)
2370 IF (ALLOCATED(itagp)) DEALLOCATE(itagp)
2371 !
2372 IF (ALLOCATED(area_max)) DEALLOCATE(area_max)
2373 IF (ALLOCATED(dm)) DEALLOCATE(dm)
2374 IF (ALLOCATED(pnod)) DEALLOCATE(pnod)
2375 IF (ALLOCATED(elarea))DEALLOCATE(elarea)
2376 IF (ALLOCATED(norm)) DEALLOCATE(norm)
2377 IF (ALLOCATED(parea)) DEALLOCATE(parea)
2378 IF (ALLOCATED(pnorm)) DEALLOCATE(pnorm)
2379 IF (ALLOCATED(pvolu)) DEALLOCATE(pvolu)
2380 IF (ALLOCATED(rhoel)) DEALLOCATE(rhoel)
2381 IF (ALLOCATED(tkel)) DEALLOCATE(tkel)
2382 IF (ALLOCATED(sspel)) DEALLOCATE(sspel)
2383 IF (ALLOCATED(pu)) DEALLOCATE(pu)
2384 IF (ALLOCATED(pw)) DEALLOCATE(pw)
2385 IF (ALLOCATED(dmini)) DEALLOCATE(dmini)
2386 IF (ALLOCATED(dmcpa)) DEALLOCATE(dmcpa)
2387 IF (ALLOCATED(dmcpb)) DEALLOCATE(dmcpb)
2388 IF (ALLOCATED(dmcpc)) DEALLOCATE(dmcpc)
2389 IF (ALLOCATED(dmrmw)) DEALLOCATE(dmrmw)
2390 IF (ALLOCATED(elsout))DEALLOCATE(elsout)
2391 IF (ALLOCATED(tmp)) DEALLOCATE(tmp)
2392 IF (ALLOCATED(poro)) DEALLOCATE(poro)
2393 IF (ALLOCATED(vola)) DEALLOCATE(vola)
2394 IF (ALLOCATED(hmin)) DEALLOCATE(hmin)
2395 IF (ALLOCATED(vel)) DEALLOCATE(vel)
2396 IF (ALLOCATED(vnod)) DEALLOCATE(vnod)
2397 IF (ALLOCATED(vgrid)) DEALLOCATE(vgrid)
2398 IF (ALLOCATED(svent)) DEALLOCATE(svent)
2399 IF (ALLOCATED(xxx)) DEALLOCATE(xxx)
2400 IF (ALLOCATED(dmcpd)) DEALLOCATE(dmcpd)
2401 IF (ALLOCATED(dmcpe)) DEALLOCATE(dmcpe)
2402 IF (ALLOCATED(dmcpf)) DEALLOCATE(dmcpf)
2403 IF (ALLOCATED(pel)) DEALLOCATE(pel)
2404 IF (ALLOCATED(dls)) DEALLOCATE(dls)
2405 IF (ALLOCATED(eldmass)) DEALLOCATE(eldmass)
2406 IF (ALLOCATED(eldehpy)) DEALLOCATE(eldehpy)
2407 IF (ALLOCATED(vvv)) DEALLOCATE(vvv)
2408 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
alpha
#define alpha
Definition eval.h:35
fvinjt6
subroutine fvinjt6(njet, ibagjet, rbagjet, npc, tf, nsensor, sensor_tab, scalt, datainj, python)
Definition fvinjt6.F:35
fvinjt8
subroutine fvinjt8(njet, ibagjet, rbagjet, npc, tf, nsensor, sensor_tab, scalt, igeo, geo, pm, ivolu, datainj, python)
Definition fvinjt8.F:36
fvinjt8_1
subroutine fvinjt8_1(njet, ibagjet, rbagjet, igeo, geo, pm, ivolu, rvolu, dmout, dhout)
Definition fvinjt8_1.F:32
fvtemp
subroutine fvtemp(ityp, efac_in, cpa_in, cpb_in, cpc_in, cpd_in, cpe_in, cpf_in, rmw_in, temp0_in, temp_in)
Definition fvtemp.F:35
fvvent0
subroutine fvvent0(elsout, aoutot, nvent, nelt, ittf, elarea, elsini, elem, itagel, svent, ibaghol, rvolu, rbaghol, poro, p, eltg, iparg, mattg, nel, porosity, ipm, pm, elbuf_tab, igroupc, igrouptg)
Definition fvvent0.F:45
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
anim_mod
Definition anim_mod.F:96
fvbag_mod
Definition fvbag_mod.F:44
fvbag_mod::fvspmd
type(fvbag_spmd), dimension(:), allocatable fvspmd
Definition fvbag_mod.F:129
fvbag_mod::fvdata
type(fvbag_data), dimension(:), allocatable fvdata
Definition fvbag_mod.F:128
fvmbag_meshcontrol_mod
Definition fvmbag_meshcontrol_mod.F:46
h3d_mod
Definition h3d_mod.F:308
message_mod
Definition message_mod.F:1249
poro
subroutine poro(geo, nodpor, ms, x, v, w, af, am, skew, weight, nporgeo)
Definition poro.F:40
spmd_fvb_gath
subroutine spmd_fvb_gath(ifv, x, xxx, xxxa, xxxsa, ido)
Definition spmd_fvb_gath.F:40
spmd_fvb_igath
subroutine spmd_fvb_igath(ifv, itabs, itabg)
Definition spmd_fvb_igath.F:34
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:889
arret
subroutine arret(nn)
Definition arret.F:87
h3d_mod::h3d_database
Definition h3d_mod.F:435