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