OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sinit22_fvm.F File Reference
#include "implicit_f.inc"
#include "param_c.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "com08_c.inc"
#include "task_c.inc"
#include "vect01_c.inc"
#include "inter22.inc"
#include "mvsiz_p.inc"
#include "comlock.inc"
#include "subvolumes.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sinit22_fvm (ixs, elbuf_tab, iparg, itab, itask, bufbric, nbric_l, x, ale_connectivity, v, nv46, veul, igrnod, ipari, igrtruss, ixt, bufmat, ipm)

Function/Subroutine Documentation

◆ sinit22_fvm()

subroutine sinit22_fvm ( integer, dimension(nixs,*), intent(in) ixs,
type(elbuf_struct_), dimension(ngroup), target elbuf_tab,
integer, dimension(nparg,*), intent(in) iparg,
integer, dimension(*), intent(in) itab,
integer, intent(in) itask,
integer, dimension(*), intent(in) bufbric,
integer nbric_l,
dimension(3,*), intent(inout) x,
type(t_ale_connectivity), intent(in) ale_connectivity,
dimension(3,*), intent(in) v,
integer, intent(in) nv46,
dimension(lveul,*), intent(in) veul,
type (group_), dimension(ngrnod) igrnod,
integer, dimension(*), intent(in) ipari,
type (group_), dimension(ngrtrus) igrtruss,
integer, dimension(nixt,*), intent(in) ixt,
dimension(*), target bufmat,
integer, dimension(npropmi,*), intent(in) ipm )

Definition at line 48 of file sinit22_fvm.F.

53C-----------------------------------------------
54C D e s c r i p t i o n
55C-----------------------------------------------
56C Interface Type22 (/INTER/TYPE22) is an FSI coupling method based on cut cell method.
57C This experimental cut cell method is not completed, abandoned, and is not an official option.
58C
59
60! This subroutine is called for FSI interface 22
61! It handles intersected bricks when lagrangian
62! surface moved in this new cycle.
63!S T E P S
64 !------------------------------------------------------------!
65 ! @00. OUTPUT CUT CELL BUFFER LIST !
66 ! @01. CINEATIC TIME STEP !
67 ! @02. INIT/PARAMETERS !
68 ! @03. INIT/MULTITHREADING !
69 ! @04. INIT/ALLOCATION !
70 ! @05. RESET TAG22 !
71 ! @06. COMPUTE UNEXTENDED main CELL VOLUMES !
72 ! @07. BUILDING BIJECTION APPLICATION for cut cells !
73 ! @08. ADJACENT BRICKS DATA IN CUT CELL BUFFER !
74 ! @09. AREA and VOLUME RATIO UPDATE for cut cells !
75 ! @10. BUILD CELL 9 !
76 ! @11. BUILD UNCUT CELL IN CUT CELL BUFFER !
77 ! @12. STORING BRICK REFERENCE VOLUME (UNCUT) !
78 ! @13. POLY 9 : MARK IF NO NODE ON A FACE !
79 ! @14. FILL MATERIAL BUFFER : %bakMAT%... !
80 ! @15. FINDING & STORING ADJACENT CELLS !
81 ! @16. BUILD POLYGONAL TRACE FOR POLY9 !
82 ! @17. PRE-TREATMENT FOR DBLE SECONDARY WITH BOTH V=0 !
83 ! @18. TAGGING SECONDARY CUT CELL !
84 ! @19. LINKING SECONDARY CELL TO A main ONE !
85 ! @20. LIST ALL SECONDARY CELLS LINKED TO A GIVEN main CELL !
86 ! @21. INDIRECT LINKING !
87 ! @22. INDIRECT LINKING !
88 ! @23. TAG SECONDARY CELLS WHICH ARE BASED ON FORMER main NODE !
89 ! @24. VOLUME MERGING BETWEEN SECONDARY CELL & ITS main ONE !
90 ! @25. MATERIAL MERGING (IF main HAS TO BE DEPORTED) !
91 ! @26. DEMERGING CRITERIA !
92 ! @27. MATERIAL DEMERGING !
93 ! @28. COMPUTE VOLD_Vnew_Cell for newly cell in buffer !
94 ! @29. COMPUTE VOLUME CHANGE FOR EACH CELL !
95 ! @30. LINK SWITCHING !
96 ! @31. LEVEL2 NEIGHBORHOOD (FOR FLUXES & CONVECTION) !
97 ! @32. COUNT POINTS ON A GIVEN POLYHEDRON FACE !
98 ! @33. LOCATE WHERE WAS main FOR NEXT CYCLE !
99 ! @34. SET IIAD22(1:NUMELS) = %ELBUG_TAB(NG)%TAG22(I) !
100 ! @35. SET MLW !
101 ! @36. RESET OLD SECONDARY LIST CONNECTIVITY !
102 ! @37. SET main STRONG NODE !
103 ! @38. DRAW SUPERCELL CONNEXION WITH 1D TRUSSES !
104 ! @39. COMPUTE DVOL PREDICTION FOR EACH CELL !
105 ! @40. BUILD SUPERCELL DVOL (PREDICTION) !
106 ! @41. FIX SUPERCELL DVOL (CORRECTION) !
107 ! @42. TIME ZERO VOLUME INIT (NOT DONE IN STARTER) !
108 ! @43. DEBUG OUTPUT !
109 ! @44. DEBUG - WRITE CUT CELL BUFFER !
110 ! @45. SUPERCELL CENTERS !
111 ! @46. MARK SUPER-CELL CENTERS WITH ORPHAN NODES !
112 ! @47. UNCUT CELLS + POLY 9 : CENTERS !
113 ! @48. MARK CELL CENTERS WITH ORPHAN NODES !
114 ! @49. MARK FACE CENTERS WITH ORPHAN NODES !
115 ! @50. DEALLOCATE !
116 !------------------------------------------------------------!
117
118!
119! BUFFERS :
120! BRICK_LIST : BUT CELL BUFFER(NIN,IB) : NIN interface id, IB:cut cell id
121C-----------------------------------------------
122C N o t e s F o r L a t e r
123C-----------------------------------------------
124C -1- Reset tag22 only for brick in previous buffer (check performances)
125C -2- remove pointers
126C -3- put the allocatables in a module
127C-----------------------------------------------
128C M o d u l e s
129C-----------------------------------------------
130 USE initbuf_mod
131 USE i22bufbric_mod
132 USE i22tri_mod
133 USE elbufdef_mod
134 USE alefvm_mod
135 USE groupdef_mod
136 USE ale_connectivity_mod
137 USE multimat_param_mod , ONLY : m51_n0phas, m51_nvphas
138 use element_mod , only : nixs,nixt
139C-----------------------------------------------
140C I m p l i c i t T y p e s
141C-----------------------------------------------
142#include "implicit_f.inc"
143C-----------------------------------------------
144C C o m m o n B l o c k s
145C-----------------------------------------------
146#include "param_c.inc"
147#include "com01_c.inc"
148#include "com04_c.inc"
149#include "com08_c.inc"
150#include "task_c.inc"
151#include "vect01_c.inc"
152#include "inter22.inc"
153#include "mvsiz_p.inc"
154#include "comlock.inc"
155#include "subvolumes.inc"
156C-----------------------------------------------
157C D u m m y A r g u m e n t s
158C-----------------------------------------------
159 INTEGER,INTENT(IN) :: IXS(NIXS,*) ,IPARG(NPARG,*),ITAB(*) ,NV46 ,BUFBRIC(*),
160 . IPARI(*) ,IXT(NIXT,*) ,ITASK ,IPM(NPROPMI,*)
161 TYPE(ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
162 my_real,INTENT(IN) :: v(3,*), veul(lveul,*)
163 my_real,INTENT(IN),TARGET :: bufmat(*)
164 my_real,INTENT(INOUT) :: x(3,*)
165!
166 TYPE (GROUP_) , DIMENSION(NGRNOD) :: IGRNOD
167 TYPE (GROUP_) , DIMENSION(NGRTRUS) :: IGRTRUSS
168 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECTIVITY
169C-----------------------------------------------
170C L o c a l V a r i a b l e s
171C-----------------------------------------------
172 TYPE(L_BUFEL_) ,POINTER :: LBUF1,LBUF2
173 INTEGER :: I,J,K0,K1,K2,JV,IDV, NBCUT, NBCUTv, NBCUTprev, NEL,NG,NFL,NBF,NBL,NBL1,NBRIC_L,NIN
174 INTEGER :: brickID,tNB,NTAG,IV,NGv,IAD22,NCELLv,ICV,IGR
175 INTEGER :: IB,IBV,IBv2,IBv_i,IBo,ICELL,ICELL2,MCELL,NCELL,MNOD,ID,ITAG(66)
176 INTEGER :: IPOS, LLT_, LLT_o,LLT_v, IDLOCv, IPOSf, IPOSiadj,ICELLv,ICELLv2
177 INTEGER :: INODES(8),INOD,INOD2,INODE,NNODES, NNODES2, ADD, ADD0, ITRIMAT, K, KV, Ko
178 LOGICAL :: lDONE,lStillTruss,lFOUND,lCYCLE,lTARGET,lStillNode, lCOND1, lCOND2
179 my_real :: volcell, volcell51(trimat), volsecnd51(trimat),volsecnd,volcell51_old(trimat)
180 my_real :: adjmain_vol(6), adjmain_face(6), adjmain_centroid(3,6),n(3,6), fac
181 INTEGER :: IDadj_MCELLv(6), IVadj_MCELLv(6), IBadj_MCELLv(6)
182 my_real,DIMENSION(:,:),ALLOCATABLE :: f
183 my_real, POINTER :: pvar, pvarv, pvaro
184 my_real, DIMENSION(:), POINTER :: pvar3
185 TYPE(G_BUFEL_) ,POINTER :: GBUF, GBUFv, GBUFo
186 TYPE(L_BUFEL_) ,POINTER :: LBUF
187 TYPE(POLY_ENTITY),DIMENSION(:), POINTER :: pIsMain ,pIsMainV
188 my_real,DIMENSION(:) , POINTER :: pfullface
189 TYPE(NODE_ENTITY),DIMENSION(:) , POINTER :: pNodWasMain
190 !INTEGER,DIMENSION(:,:), POINTER :: pNAdjCELL
191 INTEGER,DIMENSION(:,:), POINTER :: pAdjBRICK, pAdjBRICKv
192 TYPE(NODE_ENTITY),DIMENSION(:) , POINTER :: pWhereWasMain
193 TYPE(NODE_ENTITY), DIMENSION(:), POINTER :: pWhichCellNod,pWhichCellNodv
194 INTEGER , POINTER :: pMainID
195 my_real,DIMENSION(:) , POINTER :: uparam !=> BRICK_ENTITY%POLY()%Vnew
196 TYPE(POLY_ENTITY), DIMENSION(:), POINTER :: pSUBVOLv
197 my_real , POINTER :: psubvold
198 TYPE(POLY_ENTITY), DIMENSION(:), POINTER :: pSUBVOL
199 CHARACTER*6 :: char1
200 LOGICAL :: bool1, bool2
201 CHARACTER*10 :: debugMAINSECND
202 CHARACTER*10 ,ALLOCATABLE :: debugMAINSECNDv(:,:,:)
203 INTEGER,ALLOCATABLE :: IsMainV(:,:,:)
204 TYPE(BUF_MAT_) ,POINTER :: MBUF,MBUFv, MBUFo
205 INTEGER :: ICRIT_MAT_DEPORTATION, ICRIT_DEMERGE
206! INTEGER,ALLOCATABLE,DIMENSION(:,:,:):: TARGET_DATA !index1:facet - index2:adjacent_cell_number - index3:ib(cut_cell_id)
207 INTEGER,ALLOCATABLE,DIMENSION(:,:,:):: ORIGIN_DATA !index1:facet - index2:adjacent_cell_number - index3:ib(cut_cell_id),J
208 INTEGER,ALLOCATABLE,DIMENSION(:) :: Ntarget,Norigin
209 INTEGER :: MTN_,ITAR, IORIG, MAIN_TAG(6,9) !index1:Face !index2:icv=cell number
210 INTEGER :: IPLA, ISOLNOD, FM, IBm,IBmCur,IBMo,IBMold, IDM, IE, IN, MT
211 INTEGER :: ICELLTAG(9), ICELLTAG_OLD(9),IC, NAdjCELL,NAdjCELLv, IFV, FV, FV2
212 INTEGER :: IVv,MCELLv, MCELLvv, NGvv,IDLOCvv, IBvv,IFVv
213 INTEGER :: GET_UNIQUE_MAIN_CELL, LINK_WITH_UNIQUE_MAIN_CELL,IADJ,NADJ,NADJv,IE_M
214 INTEGER :: NsecndNOD, NP, NP_NODE, JJ, NINTP(6,9), NNOD(6,9), NN, SECid, IDLOC, NewMnod(8), MLW
215 INTEGER :: ITASKv, NTRUS, NPOINT, IRES(2), NewInBuffer
216 INTEGER,ALLOCATABLE,DIMENSION(:,:) :: DESTROY_DATA
217 INTEGER :: I22LOOP,IT,IUNLINK,ITASKB, J1,J2, NTAR, INod_W, IDEMERGE, INod_W_old, ISECND, IAD0,II
218 INTEGER :: FV_old, NumSECND, NumSECND2, IC1, IC2, ITARGET, NumTarget, WasCut, LID, IFAC, IEDG
219 INTEGER :: Cod1, Cod2, Cod3, Icompl, Poly9woNodes, Poly9woNodesV, ISGN(2), ICODE, OLD_ICODE
220 INTEGER :: mainID,NC(8),MOLD,MNEW, IB2,NUM,IADBUF,NUPARAM,NGm,IDLOCm,ICELLm, NP_(9),NODE_ID
221 LOGICAL :: debug_outp,debug_outp2, lSKIP
222
223 my_real,ALLOCATABLE, DIMENSION(:,:) :: vol51, vol51v !index1:IB index2:ITRIMAT
224 my_real,ALLOCATABLE, DIMENSION(:) :: uvar,uvar_adv !(M51_N0PHAS+TRIMAT*M51_NVPHAS+I22LAW37) !NUVAR=M51_N0PHAS+TRIMAT*M51_NVPHAS
225 my_real :: vfrac(4),som, som_(4), somi, adjface(nv46,nadj_f), delta, ppoint(3)
226 my_real :: sumface, vectmp(3), volorig(24), pointtmp(3), point0(3), cut_point(3)
227 my_real :: cut_vel(3)
228 my_real :: eint, eintv, rho,rho_u(3),mom(3),rhov,sigv(6), sig(6),volv,vol,vol_m,vol_s
229 my_real :: var, var_,var__, var_6(6), var_vf(4), vold_phase(4)
230 my_real :: tmp(6),dxmin,vmax,ratiof,ratiofv,ratio,ratiov,uncutvol,uncutvolv
231 my_real :: vuncut, vi,vj, m(9,9) !index1:polyhedra in current conformation, index2:polyhedra in previous confromation
232 my_real :: vsum(3) , n_(3), vnew, vold, dvol_numeric, dvol_predic
233 my_real :: sgn, dvi, dvii, face, face9, m_tot, m_liq, m_toto,m_liqo, rho10, rho20, mfrac
234 my_real :: center(3,8)
235
236 my_real :: det, df11, df12, df21, df22, f1, f2, p1, p2, drho1, drho2, error, tol, r1, psh, pmin
237 my_real :: c1, gam, p0, p, rho1, rho2, mas, mas1, mas2, ssp, ssp1, ssp2, rhoc2_1, rhoc2_2, rhoc2
238 my_real :: alp, alpo, beta, volcellold, vel_m(3), surf_s, norm_s(3), adv, norm, vm, vs
239 my_real :: rho_adv, eint_adv, sig_adv(6), mom_adv(3), zm(3),zs(3), volratio
240 INTEGER :: ITER, NITER, LLTo, LLTm, NGo, IADV
241 INTEGER, ALLOCATABLE, DIMENSION(:) :: ORDER, VALUE
242
243 INTEGER :: NC1,NC2,NC3,NC4,NC5,NC6,NC7,NC8
244 INTEGER :: IAD2, IAD3, LGTH2, LGTH3
245
246 INTERFACE
247 FUNCTION i22aera(NPTS, P, C)
248 my_real :: i22aera(3)
249 INTEGER, INTENT(IN) :: NPTS
250 my_real, INTENT(IN) :: p(3,npts)
251 my_real, INTENT(INOUT) :: c(3)
252 END FUNCTION
253 END INTERFACE
254
255 TYPE pointer_array_r
256 my_real,DIMENSION(:), POINTER :: p
257 END TYPE
258
259 TYPE pointer_array_i
260 INTEGER,DIMENSION(:), POINTER :: p
261 END TYPE
262
263 TYPE(POLY_ENTITY), DIMENSION(:), POINTER :: pFACE ! WARNING : all pointers begin with 1 !!! index 0 impossible pFACE(1:7,:)=>%Face_Cell(0:6,:) <
264 TYPE(POINTER_ARRAY_R) :: pFACEv(9)
265 TYPE(POINTER_ARRAY_I) :: pListNodID(9)
266 TYPE(POINTER_ARRAY_I) :: pListNodIDv(9)
267 !TYPE(POINTER_ARRAY_I) :: pAdjCELL
268
269 INTEGER ICF(4,6)
270 DATA icf/1,4,3,2,3,4,8,7,5,6,7,8,1,2,6,5,2,3,7,6,1,5,8,4/
271
272 my_real ::
273 . aj7(mvsiz), aj8(mvsiz) , aj9(mvsiz),
274 . aji1, aji2, aji3,
275 . aji4, aji5, aji6,
276 . aji7, aji8, aji9,
277 . x17 , x28 , x35 , x46,
278 . y17 , y28 , y35 , y46,
279 . z17 , z28 , z35 , z46,
280 . jac_59_68(mvsiz), jac_67_49(mvsiz), jac_48_57(mvsiz),
281 . aj12, aj45, aj78,
282 . a17 , a28 ,
283 . b17 , b28 ,
284 . c17 , c28 ,
285 . x1(mvsiz),x2(mvsiz),x3(mvsiz),x4(mvsiz),x5(mvsiz),x6(mvsiz),x7(mvsiz),x8(mvsiz),
286 . y1(mvsiz),y2(mvsiz),y3(mvsiz),y4(mvsiz),y5(mvsiz),y6(mvsiz),y7(mvsiz),y8(mvsiz),
287 . z1(mvsiz),z2(mvsiz),z3(mvsiz),z4(mvsiz),z5(mvsiz),z6(mvsiz),z7(mvsiz),z8(mvsiz),
288 . dett(mvsiz),
289 . aj1(mvsiz),aj2(mvsiz),aj3(mvsiz),
290 . aj4(mvsiz),aj5(mvsiz),aj6(mvsiz),hxp,hyp,hzp,
291 . x1_,x2_,x3_,x4_,x5_,x6_,x7_,x8_,
292 . y1_,y2_,y3_,y4_,y5_,y6_,y7_,y8_,
293 . z1_,z2_,z3_,z4_,z5_,z6_,z7_,z8_
294
295C-----------------------------------------------
296C P r e c o n d i t i o n s
297C-----------------------------------------------
298 !IF(INT22 == 0) RETURN !already checked if this subroutine is called
299C-----------------------------------------------
300C S o u r c e L i n e s
301C-----------------------------------------------
302 !------------------------------------------------------------!
303 ! @00. OUTPUT CUT CELL BUFFER LIST !
304 !------------------------------------------------------------!
305 IF(ibug22_ccbuflist /= 0)THEN
306 nin = 1
307 IF(itask==0)THEN
308 write (*,fmt='(A, 1000I7)')"CUT CELL BUFFER : ", ixs(11,brick_list(nin,1:nb)%id)
309 ENDIF
310 ENDIF
311
312 !------------------------------------------------------------!
313 ! @01. CINEATIC TIME STEP !
314 !------------------------------------------------------------!
315 idt_int22 = 1 !kinematic time step by default. Will be replaced by 0 if exist i such as DTIX(I) < DT22_MIN (mqviscb.F)
316 dxmin = minval(dx22min_l(0:nthread-1))
317 vmax = maxval(v22max_l(0:nthread-1))
318 IF(vmax == zero)THEN
319 dt22_min = ep30
320 ELSE
321 dt22_min = dxmin/ncross22 / vmax
322 ENDIF
323 if(ibug22_dtmin==1)print *, "dt22_min = ", dt22_min
324 !!sIF(DT1==ZERO .AND. ITASK==0) print *, "INTER22, Tolerance ratio ",RATIO22
325
326 !------------------------------------------------------------!
327 ! @02. INIT/PARAMETERS !
328 !------------------------------------------------------------!
329 iadv = 0 ! 0: advection step off 1: advection step on !When demerging a new main, advection is done from old main cell (possibly relocated due to a merging process) to new one.
330 nin = 1
331 dbvol = zero
332 dbmass = zero
333 i22loop = 0 !infinite loop detected
334 imergel(itask) = 0
335 IF(trimat<=0)THEN
336 ALLOCATE(uvar(i22law37))
337 ALLOCATE(uvar_adv(i22law37))
338 ELSE
339 ALLOCATE(uvar(m51_n0phas+trimat*m51_nvphas))
340 ALLOCATE(uvar_adv(m51_n0phas+trimat*m51_nvphas))
341 ENDIF
342
343 !------------------------------------------------------------!
344 ! @03. INIT/MULTITHREADING !
345 !------------------------------------------------------------!
346 nbf = 1+itask*nb/nthread
347 nbl = (itask+1)*nb/nthread
348 nbl = min(nbl,nb)
349 tnb = nbl-nbf+1
350 i22_degenerated = 0
351
352 !!!!!!!!!!!!!!!!!!!!!!DEBUG OUTPUT!!!!!!!!!!!!!!!!!!!!!!!
353 debug_outp = .false.
354 debug_outp2= .false.
355
356 do ib=nbf,nbl
357 ie = brick_list(nin,ib)%id
358 if (ie == ibug22_sinit .or. ibug22_sinit == -1)then
359 debug_outp = .true.
360 exit
361 endif
362 enddo
363
364 if(itask==0.AND.debug_outp)then
365 print *, ""
366 print *, " |----------i22sinit.F-----------|"
367 print *, " | INITIALIZATION SUBROUTINE |"
368 print *, " |-------------------------------|"
369 char1 = ' '
370 end if
371 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
372
373 !------------------------------------------------------------!
374 ! @04. INIT/ALLOCATION !
375 !------------------------------------------------------------!
376 ALLOCATE (debugmainsecndv(nbl-nbf+1,6,9))
377! ALLOCATE (TARGET_DATA(tNB,9,1:2)) !index1:IB, index2:1:Ntarget<=9, index3:1=localface_2=adjacent_cell
378 ALLOCATE (origin_data(nbf:nbl,9,1:3)) !index1:IB, index2:1:Ntarget<=9, index3:1=localface_2=adjacent_cell
379 ALLOCATE (ismainv(nbf:nbl,6,9))
380 ALLOCATE (f(6,nbf:nbl))
381 ALLOCATE (vol51(nbf:nbl,trimat),vol51v(nbf:nbl,trimat))
382c ALLOCATE (Ntarget(NBF:NBL))
383 ALLOCATE (norigin(nbf:nbl))
384 !-------------------!
385 CALL my_barrier()
386 !-------------------!
387 !------------------------------------------------------------!
388 ! @05. RESET TAG22 !
389 !------------------------------------------------------------!
390 !TAG22(I) is adress in cut cell buffer "BRICK_LIST(NIN,TAG22(I))"
391 !for interface number NIN, and element I from current group.
392 !TAG22==0 means elem I is not in cut cell buffer.
393 DO ng=itask+1,ngroup,nthread
394 IF(iparg(8,ng) /= 1) THEN
395 CALL initbuf(
396 1 iparg ,ng ,
397 2 mtn ,nel ,nft ,iad ,ity ,
398 3 npt ,jale ,ismstr ,jeul ,jtur ,
399 4 jthe ,jlag ,jmult ,jhbe ,jivf ,
400 5 nvaux ,jpor ,jcvt ,jclose ,ipla ,
401 6 irep ,iint ,igtyp ,israt ,isrot ,
402 7 icsen ,isorth ,isorthg ,ifailure,jsms )
403 IF(jlag /= 1 .AND. ity<=2) THEN
404 IF (mtn /= 0 .AND. iparg(64,ng)==0) THEN
405 lft = 1
406 llt = nel
407 isolnod = iparg(28,ng)
408 IF (ity == 1 .AND. isolnod /= 4) THEN
409 gbuf => elbuf_tab(ng)%GBUF
410 gbuf%TAG22(lft:llt) = 0
411 ENDIF
412 ENDIF
413 ENDIF
414 ENDIF
415 ENDDO
416 !------------------------------------------------------------!
417 ! @06. COMPUTE UNEXTENDED main CELL VOLUMES !
418 !------------------------------------------------------------!
419 DO ng=itask+1,ngroup,nthread
420 IF(iparg(8,ng) /= 1) THEN
421 CALL initbuf(
422 1 iparg ,ng ,
423 2 mtn ,nel ,nft ,iad ,ity ,
424 3 npt ,jale ,ismstr ,jeul ,jtur ,
425 4 jthe ,jlag ,jmult ,jhbe ,jivf ,
426 5 nvaux ,jpor ,jcvt ,jclose ,ipla ,
427 6 irep ,iint ,igtyp ,israt ,isrot ,
428 7 icsen ,isorth ,isorthg ,ifailure,jsms )
429 IF(jlag /= 1 .AND. ity<=2) THEN
430 IF (mtn /= 0 .AND. iparg(64,ng)==0) THEN
431 lft = 1
432 llt = nel
433 isolnod = iparg(28,ng)
434 IF (ity == 1 .AND. isolnod /= 4) THEN
435 gbuf => elbuf_tab(ng)%GBUF
436 IF(jeul/=0)THEN
437 IF(integ8==0)THEN
438 gbuf%VOL(lft:llt) = veul(32,nft+lft:nft+llt)
439 ELSEIF(integ8==1)THEN
440 gbuf%VOL(lft:llt) = veul(52,nft+lft:nft+llt)
441 endif!INTEG8
442 ELSE
443 DO i=lft,llt
444 ii = i+nft
445 x1(i) = x(1,ixs(2,ii)) ; y1(i) = x(2,ixs(2,ii)) ; z1(i) = x(3,ixs(2,ii)) ;
446 x2(i) = x(1,ixs(3,ii)) ; y2(i) = x(2,ixs(3,ii)) ; z2(i) = x(3,ixs(3,ii)) ;
447 x3(i) = x(1,ixs(4,ii)) ; y3(i) = x(2,ixs(4,ii)) ; z3(i) = x(3,ixs(4,ii)) ;
448 x4(i) = x(1,ixs(5,ii)) ; y4(i) = x(2,ixs(5,ii)) ; z4(i) = x(3,ixs(5,ii)) ;
449 x5(i) = x(1,ixs(6,ii)) ; y5(i) = x(2,ixs(6,ii)) ; z5(i) = x(3,ixs(6,ii)) ;
450 x6(i) = x(1,ixs(7,ii)) ; y6(i) = x(2,ixs(7,ii)) ; z6(i) = x(3,ixs(7,ii)) ;
451 x7(i) = x(1,ixs(8,ii)) ; y7(i) = x(2,ixs(8,ii)) ; z7(i) = x(3,ixs(8,ii)) ;
452 x8(i) = x(1,ixs(9,ii)) ; y8(i) = x(2,ixs(9,ii)) ; z8(i) = x(3,ixs(9,ii)) ;
453 ENDDO
454 DO i=lft,llt
455 x17=x7(i)-x1(i)
456 x28=x8(i)-x2(i)
457 x35=x5(i)-x3(i)
458 x46=x6(i)-x4(i)
459 y17=y7(i)-y1(i)
460 y28=y8(i)-y2(i)
461 y35=y5(i)-y3(i)
462 y46=y6(i)-y4(i)
463 z17=z7(i)-z1(i)
464 z28=z8(i)-z2(i)
465 z35=z5(i)-z3(i)
466 z46=z6(i)-z4(i)
467 aj1(i)=x17+x28-x35-x46
468 aj2(i)=y17+y28-y35-y46
469 aj3(i)=z17+z28-z35-z46
470 a17=x17+x46
471 a28=x28+x35
472 b17=y17+y46
473 b28=y28+y35
474 c17=z17+z46
475 c28=z28+z35
476 aj4(i)=a17+a28
477 aj5(i)=b17+b28
478 aj6(i)=c17+c28
479 aj7(i)=a17-a28
480 aj8(i)=b17-b28
481 aj9(i)=c17-c28
482 ENDDO
483 DO i=lft,llt
484 jac_59_68(i)=aj5(i)*aj9(i)-aj6(i)*aj8(i)
485 jac_67_49(i)=aj6(i)*aj7(i)-aj4(i)*aj9(i)
486 jac_48_57(i)=aj4(i)*aj8(i)-aj5(i)*aj7(i)
487 ENDDO
488 DO i=lft,llt
489 dett(i)=one_over_64*(aj1(i)*jac_59_68(i)+aj2(i)*jac_67_49(i)+aj3(i)*jac_48_57(i))
490 ENDDO
491 gbuf%VOL(lft:llt) = dett(lft:llt)
492 endif!JEUL/=0
493 ENDIF
494 ENDIF
495 ENDIF
496 ENDIF
497 ENDDO
498 !------------------------------------------------------------!
499 ! @07. BUILDING BIJECTION APPLICATION for cut cells !
500 ! (local array) <--> (global buffer) !
501 ! NBRIC <--> (I,NG) !
502 !------------------------------------------------------------!
503 !-------------------!
504 CALL my_barrier() !need %TAG22 to be init to 0. (initialisation is multithreading segmenting 1:NG, not 1:NB
505 !-------------------!
506 DO ib=nbf,nbl !1,NB
507 idb(ib) = brick_list(nin,ib)%ID
508 brick_list(nin,ib)%ITASK = itask
509 ENDDO
510 DO ib=nbf,nbl !1,NB
511 !recuperer NG,I
512 !From IB get I and NG
513 ldone = .false.
514 DO ng=1,ngroup
515 nel=iparg(2,ng)
516 nfl=iparg(3,ng)
517 IF((idb(ib)>nfl).AND.(idb(ib)<=nfl+nel))THEN
518 IF(iparg(11,ng)==1) jeul = 1
519 IF(iparg(7,ng)==1) jale = 1
520 idlocb(ib) = idb(ib) - nfl
521 ngb(ib) = ng
522 brick_list(nin,ib)%NG = ngb(ib)
523 brick_list(nin,ib)%IDLOC = idlocb(ib)
524 nelb(ib) = nel
525 ldone =.true.
526 gbuf => elbuf_tab(ngb(ib))%GBUF
527 nbcut = brick_list(nin,ib)%NBCUT
528 gbuf%TAG22(idlocb(ib)) = ib ! les tags ont deja ete remis a zero
529 EXIT !next I
530 ELSE
531 cycle !next NG
532 ENDIF
533 if (.NOT.(ldone))then
534 write( *,*) "int 22 : error in group sorting"
535 stop 2201
536 end if
537 ENDDO !next NG
538 ENDDO !next i
539 !-------------------!
540 CALL my_barrier()
541 !-------------------!
542 !------------------------------------------------------------!
543 ! @08. ADJACENT BRICKS DATA IN CUT CELL BUFFER !
544 !------------------------------------------------------------!
545 DO ib=nbf,nbl
546 brickid = idb(ib)
547 iad2 = ale_connectivity%ee_connect%iad_connect(brickid)
548 lgth2 = ale_connectivity%ee_connect%iad_connect(brickid+1) -
549 . ale_connectivity%ee_connect%iad_connect(brickid)
550 DO j=1,lgth2
551 idv = ale_connectivity%ee_connect%connected(iad2 + j - 1)
552 brick_list(nin,ib)%Adjacent_Brick(j,1) = idv !IV
553 brick_list(nin,ib)%Adjacent_Brick(j,2) = 0 !NGV
554 brick_list(nin,ib)%Adjacent_Brick(j,3) = 0 !IDLOCV
555 brick_list(nin,ib)%Adjacent_Brick(j,4) = 0 !IBV
556 brick_list(nin,ib)%Adjacent_Brick(j,5) = 0 !IFV
557 IF(idv>0)THEN
558 iad3 = ale_connectivity%ee_connect%iad_connect(idv)
559 lgth3 = ale_connectivity%ee_connect%iad_connect(idv+1) -
560 . ale_connectivity%ee_connect%iad_connect(idv)
561 DO jv=1,lgth3
562 IF(ale_connectivity%ee_connect%connected(iad3 + jv - 1)==brickid)THEN
563 brick_list(nin,ib)%Adjacent_Brick(j,5) = jv
564 EXIT
565 ENDIF
566 enddo!next JV
567 CALL get_group_id(idv,ngv,idlocv,iparg)
568 gbufv => elbuf_tab(ngv)%GBUF
569 iad22 = gbufv%TAG22(idlocv)
570 brick_list(nin,ib)%Adjacent_Brick(j,2) = ngv
571 brick_list(nin,ib)%Adjacent_Brick(j,3) = idlocv
572 brick_list(nin,ib)%Adjacent_Brick(j,4) = iad22
573 IF (iad22==0 .AND. brick_list(nin,ib)%NBCUT>0)THEN
574 print *, "**error : inter22 : Lagrangian Surface seems to"
575 print *, " reach eulerian boundary domain. "
576 print *, " Check Surface location and GRBRICK definition"
577 print *, " near related Brick_ID =", ixs(11,brick_list(nin,ib)%id)
578 !print *, "ITASK=", ITASK
579
580 stop 2202 !means that a real cut cell has an adjacent brick which is not in cut cell buffer : NOT EXPECTED unless if /INTER/TYPE22 GrBrick is not correctly defined.
581 ENDIF
582 ELSEIF(idv==0)THEN
583 brick_list(nin,ib)%Adjacent_Brick(j,2) = 0
584 brick_list(nin,ib)%Adjacent_Brick(j,3) = 0
585 brick_list(nin,ib)%Adjacent_Brick(j,4) = 0
586 ELSEIF(idv<0)THEN
587 write (*,*) "unavailable in SPMD"
588 stop 2203
589 ENDIF
590 ENDDO
591 ENDDO
592
593 !------------------------------------------------------------!
594 ! @09. AREA and VOLUME RATIO UPDATE for cut cells !
595 !------------------------------------------------------------!
596 DO ib=nbf,nbl
597 brickid = idb(ib)
598 IF(i22_aleul==2)THEN
599 n(1:3,1) = (/ veul(14,brickid) , veul(20,brickid) , veul(26,brickid) /)
600 n(1:3,2) = (/ veul(15,brickid) , veul(21,brickid) , veul(27,brickid) /)
601 n(1:3,3) = (/ veul(16,brickid) , veul(22,brickid) , veul(28,brickid) /)
602 n(1:3,4) = (/ veul(17,brickid) , veul(23,brickid) , veul(29,brickid) /)
603 n(1:3,5) = (/ veul(18,brickid) , veul(24,brickid) , veul(30,brickid) /)
604 n(1:3,6) = (/ veul(19,brickid) , veul(25,brickid) , veul(31,brickid) /)
605 brick_list(nin,ib)%N(1,1:3) = n(1:3,1) / sqrt(sum(n(1:3,1)*n(1:3,1)))
606 brick_list(nin,ib)%N(2,1:3) = n(1:3,2) / sqrt(sum(n(1:3,2)*n(1:3,2)))
607 brick_list(nin,ib)%N(3,1:3) = n(1:3,3) / sqrt(sum(n(1:3,3)*n(1:3,3)))
608 brick_list(nin,ib)%N(4,1:3) = n(1:3,4) / sqrt(sum(n(1:3,4)*n(1:3,4)))
609 brick_list(nin,ib)%N(5,1:3) = n(1:3,5) / sqrt(sum(n(1:3,5)*n(1:3,5)))
610 brick_list(nin,ib)%N(6,1:3) = n(1:3,6) / sqrt(sum(n(1:3,6)*n(1:3,6)))
611 ELSEIF(i22_aleul==1)THEN
612 ii = brickid
613 !---8 local node numbers NC1 TO NC8 for solid element I ---!
614 nc1=ixs(2,ii)
615 nc2=ixs(3,ii)
616 nc3=ixs(4,ii)
617 nc4=ixs(5,ii)
618 nc5=ixs(6,ii)
619 nc6=ixs(7,ii)
620 nc7=ixs(8,ii)
621 nc8=ixs(9,ii)
622 !
623 !---Coordinates of the 8 nodes
624 x1_=x(1,nc1)
625 y1_=x(2,nc1)
626 z1_=x(3,nc1)
627 !
628 x2_=x(1,nc2)
629 y2_=x(2,nc2)
630 z2_=x(3,nc2)
631 !
632 x3_=x(1,nc3)
633 y3_=x(2,nc3)
634 z3_=x(3,nc3)
635 !
636 x4_=x(1,nc4)
637 y4_=x(2,nc4)
638 z4_=x(3,nc4)
639 !
640 x5_=x(1,nc5)
641 y5_=x(2,nc5)
642 z5_=x(3,nc5)
643 !
644 x6_=x(1,nc6)
645 y6_=x(2,nc6)
646 z6_=x(3,nc6)
647 !
648 x7_=x(1,nc7)
649 y7_=x(2,nc7)
650 z7_=x(3,nc7)
651 !
652 x8_=x(1,nc8)
653 y8_=x(2,nc8)
654 z8_=x(3,nc8)
655 !
656 ! Face-1
657 n(1,1)=(y3_-y1_)*(z2_-z4_) - (z3_-z1_)*(y2_-y4_)
658 n(2,1)=(z3_-z1_)*(x2_-x4_) - (x3_-x1_)*(z2_-z4_)
659 n(3,1)=(x3_-x1_)*(y2_-y4_) - (y3_-y1_)*(x2_-x4_)
660 ! Face-2
661 n(1,2)=(y7_-y4_)*(z3_-z8_) - (z7_-z4_)*(y3_-y8_)
662 n(2,2)=(z7_-z4_)*(x3_-x8_) - (x7_-x4_)*(z3_-z8_)
663 n(3,2)=(x7_-x4_)*(y3_-y8_) - (y7_-y4_)*(x3_-x8_)
664 ! Face-3
665 n(1,3)=(y6_-y8_)*(z7_-z5_) - (z6_-z8_)*(y7_-y5_)
666 n(2,3)=(z6_-z8_)*(x7_-x5_) - (x6_-x8_)*(z7_-z5_)
667 n(3,3)=(x6_-x8_)*(y7_-y5_) - (y6_-y8_)*(x7_-x5_)
668 ! Face-4
669 n(1,4)=(y2_-y5_)*(z6_-z1_) - (z2_-z5_)*(y6_-y1_)
670 n(2,4)=(z2_-z5_)*(x6_-x1_) - (x2_-x5_)*(z6_-z1_)
671 n(3,4)=(x2_-x5_)*(y6_-y1_) - (y2_-y5_)*(x6_-x1_)
672 ! Face-5
673 n(1,5)=(y7_-y2_)*(z6_-z3_) - (z7_-z2_)*(y6_-y3_)
674 n(2,5)=(z7_-z2_)*(x6_-x3_) - (x7_-x2_)*(z6_-z3_)
675 n(3,5)=(x7_-x2_)*(y6_-y3_) - (y7_-y2_)*(x6_-x3_)
676 ! Face-6
677 n(1,6)=(y8_-y1_)*(z4_-z5_) - (z8_-z1_)*(y4_-y5_)
678 n(2,6)=(z8_-z1_)*(x4_-x5_) - (x8_-x1_)*(z4_-z5_)
679 n(3,6)=(x8_-x1_)*(y4_-y5_) - (y8_-y1_)*(x4_-x5_)
680 !
681 brick_list(nin,ib)%N(1,1:3) = n(1:3,1) / sqrt(sum(n(1:3,1)*n(1:3,1)))
682 brick_list(nin,ib)%N(2,1:3) = n(1:3,2) / sqrt(sum(n(1:3,2)*n(1:3,2)))
683 brick_list(nin,ib)%N(3,1:3) = n(1:3,3) / sqrt(sum(n(1:3,3)*n(1:3,3)))
684 brick_list(nin,ib)%N(4,1:3) = n(1:3,4) / sqrt(sum(n(1:3,4)*n(1:3,4)))
685 brick_list(nin,ib)%N(5,1:3) = n(1:3,5) / sqrt(sum(n(1:3,5)*n(1:3,5)))
686 brick_list(nin,ib)%N(6,1:3) = n(1:3,6) / sqrt(sum(n(1:3,6)*n(1:3,6)))
687 ENDIF
688
689 DO j=1,nv46
690 ! compute Brick Face New Area !suplicated eflux3/aflux3, to optimize
691 f(j,ib) = half * sqrt( sum( n(:,j) * n(:,j) ) ) ! N is 2Sn
692 ENDDO
693 pface(1:9) => brick_list(nin,ib)%POLY(1:9)
694! pFACE(2)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(2)%FACE(1:6)%Surf
695! pFACE(3)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(3)%FACE(1:6)%Surf
696! pFACE(4)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(4)%FACE(1:6)%Surf
697! pFACE(5)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(5)%FACE(1:6)%Surf
698! pFACE(6)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(6)%FACE(1:6)%Surf
699! pFACE(7)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(7)%FACE(1:6)%Surf
700! pFACE(8)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(8)%FACE(1:6)%Surf
701! pFACE(9)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(9)%FACE(1:6)%Surf
702 pfullface => brick_list(nin,ib)%Face_BRICK(1:6)
703 pfullface(1:6) = f(1:6,ib)
704 IF(brick_list(nin,ib)%NBCUT==0)THEN
705 pface(1)%FACE(1:6)%Surf = f(1:6,ib)
706 cycle
707 ENDIF
708 ncell = brick_list(nin,ib)%NBCUT !cell 9 not yet taken into account.
709 DO icell=1,ncell
710 !FACES UPDATE
711 IF(brick_list(nin,ib)%POLY(icell)%Vnew<zero)THEN
712 DO j=1,nv46
713 pface(icell)%FACE(j)%Surf = f(j,ib) + pface(icell)%FACE(j)%Surf
714 if(pface(icell)%FACE(j)%Surf<zero)then
715 write (*,*) "**error : inter22 negative cell face"
716 endif
717 enddo! next J
718 ENDIF
719 !VOLUME UPDATE
720 vol = brick_list(nin,ib)%POLY(icell)%Vnew
721 IF(vol<zero)THEN
722 gbuf => elbuf_tab(ngb(ib))%GBUF
723 brick_list(nin,ib)%POLY(icell)%Vnew = gbuf%VOL(idlocb(ib)) + vol
724 ENDIF
725 volratio = brick_list(nin,ib)%POLY(icell)%Vnew / elbuf_tab(ngb(ib))%GBUF%VOL(idlocb(ib))
726 ! THRESHOLD VOLUME : Smax = O(V**(2/3)), dh= EM03 * O(V**(1/3)) => Veps = EM04*O(V) => VOLratio = Veps/V = EM04 (EM04 => dh = 0.01 % * average_edge)
727 IF(abs(volratio) <= em04)THEN
728 i22_degenerated = 1
729 ENDIF
730 enddo!next ICELL
731 ENDDO !next IB (cut cell)
732
733
734
735
736 !------------------------------------------------------------!
737 ! @10. BUILD CELL 9 !
738 ! which is complementary one : !
739 ! (faces,volumes,brick nodes and num points) !
740 !------------------------------------------------------------!
741 DO ib=nbf,nbl
742 plistnodid(1)%p(1:8) => brick_list(nin,ib)%POLY(1)%ListNodID(1:8)
743 plistnodid(2)%p(1:8) => brick_list(nin,ib)%POLY(2)%ListNodID(1:8)
744 plistnodid(3)%p(1:8) => brick_list(nin,ib)%POLY(3)%ListNodID(1:8)
745 plistnodid(4)%p(1:8) => brick_list(nin,ib)%POLY(4)%ListNodID(1:8)
746 plistnodid(5)%p(1:8) => brick_list(nin,ib)%POLY(5)%ListNodID(1:8)
747 plistnodid(6)%p(1:8) => brick_list(nin,ib)%POLY(6)%ListNodID(1:8)
748 plistnodid(7)%p(1:8) => brick_list(nin,ib)%POLY(7)%ListNodID(1:8)
749 plistnodid(8)%p(1:8) => brick_list(nin,ib)%POLY(8)%ListNodID(1:8)
750 plistnodid(9)%p(1:8) => brick_list(nin,ib)%POLY(9)%ListNodID(1:8)
751 psubvol(1:9) => brick_list(nin,ib)%POLY(1:9)
752 pwhichcellnod(1:8) => brick_list(nin,ib)%NODE(1:8)
753 pface(1:9) => brick_list(nin,ib)%POLY(1:9)!%FACE(1:6)%Surf
754! pFACE(2)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(2)%FACE(1:6)%Surf
755! pFACE(3)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(3)%FACE(1:6)%Surf
756! pFACE(4)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(4)%FACE(1:6)%Surf
757! pface(5)%p(1:6) => brick_list(nin,ib)%POLY(5)%FACE(1:6)%Surf
758! pFACE(6)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(6)%FACE(1:6)%Surf
759! pFACE(7)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(7)%FACE(1:6)%Surf
760! pFACE(8)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(8)%FACE(1:6)%Surf
761! pface(9)%p(1:6) => brick_list(nin,ib)%POLY(9)%FACE(1:6)%Surf
762 id = brick_list(nin,ib)%ID
763 ncell = brick_list(nin,ib)%NBCUT
764 cod1 = iabs(brick_list(nin,ib)%SecID_CELL(1))
765 cod2 = iabs(brick_list(nin,ib)%SecID_CELL(2))
766 IF(ncell/=0)THEN !if brick cutting.
767 !---VOLUME AND CONNECTIVITY
768 brick_list(nin,ib)%POLY(9)%FACE(1)%Surf = f(1,ib)-sum( (/ (brick_list(nin,ib)%POLY(k)%FACE(1)%Surf,k=1,ncell)/) )
769 brick_list(nin,ib)%POLY(9)%FACE(2)%Surf = f(2,ib)-sum( (/ (brick_list(nin,ib)%POLY(k)%FACE(2)%Surf,k=1,ncell)/) )
770 brick_list(nin,ib)%POLY(9)%FACE(3)%Surf = f(3,ib)-sum( (/ (brick_list(nin,ib)%POLY(k)%FACE(3)%Surf,k=1,ncell)/) )
771 brick_list(nin,ib)%POLY(9)%FACE(4)%Surf = f(4,ib)-sum( (/ (brick_list(nin,ib)%POLY(k)%FACE(4)%Surf,k=1,ncell)/) )
772 brick_list(nin,ib)%POLY(9)%FACE(5)%Surf = f(5,ib)-sum( (/ (brick_list(nin,ib)%POLY(k)%FACE(5)%Surf,k=1,ncell)/) )
773 brick_list(nin,ib)%POLY(9)%FACE(6)%Surf = f(6,ib)-sum( (/ (brick_list(nin,ib)%POLY(k)%FACE(6)%Surf,k=1,ncell)/) )
774 brick_list(nin,ib)%POLY(9)%FACE0%Surf = sum( (/ (brick_list(nin,ib)%POLY(k)%FACE0%Surf,k=1,ncell) /) ) !sum of cut surfaces
775 brick_list(nin,ib)%POLY(9)%Vnew = -sum(brick_list(nin,ib)%POLY(1:ncell)%Vnew)
776 !print *, "pSUBVOL(9)=", pSUBVOL(9)
777 brick_list(nin,ib)%POLY(9)%NumNOD = 0 !should noty be needed
778 brick_list(nin,ib)%POLY(9)%NumPOINT = 0 !should noty be needed
779 mnod = 0
780 ppoint(1:3) = zero
781 !saving brick nodes used by remaining polyhedron (complementary one)
782 DO j=1,8
783 IF(pwhichcellnod(j)%WhichCell/=0) cycle !a tagged brick node is already used by a polyhedron
784 pwhichcellnod(j)%WhichCell = 9
785 mnod = mnod + 1
786 plistnodid(9)%p(mnod) = j !IXS(J+1,ID)
787 ppoint(1) = ppoint(1) + x(1,ixs(1+j,id))
788 ppoint(2) = ppoint(2) + x(2,ixs(1+j,id))
789 ppoint(3) = ppoint(3) + x(3,ixs(1+j,id))
790 ENDDO
791 brick_list(nin,ib)%POLY(9)%NumNOD = mnod
792 brick_list(nin,ib)%POLY(9)%NumPOINT =
793 . sum(brick_list(nin,ib)%POLY(1:ncell)%NumPOINT ) - sum(brick_list(nin,ib)%POLY(1:ncell)%NumNOD) + mnod
794 gbuf => elbuf_tab(ngb(ib))%GBUF
795 brick_list(nin,ib)%POLY(9)%Vnew = gbuf%VOL(idlocb(ib)) + brick_list(nin,ib)%POLY(9)%Vnew
796 uncutvol = gbuf%VOL(idlocb(ib))
797 IF(abs(brick_list(nin,ib)%POLY(9)%Vnew /uncutvol) <em04
798 . .OR. abs(brick_list(nin,ib)%POLY(1)%Vnew /uncutvol) <em04)THEN
799 i22_degenerated = 1
800 ENDIF
801 !---AND ALSO CELL CENTER
802 !MNOD brick nodes above + K1 intersection points below.
803 npoint = brick_list(nin,ib)%POLY(9)%NumPOINT
804 k1 = 0
805 DO j=1,12
806 nbcut = brick_list(nin,ib)%EDGE(j)%NBCUT
807 DO i=1,nbcut
808 pointtmp(1:3) = brick_list(nin,ib)%EDGE(j)%CUTPOINT(1:3,i)
809 ppoint(1) = ppoint(1) + pointtmp(1)
810 ppoint(2) = ppoint(2) + pointtmp(2)
811 ppoint(3) = ppoint(3) + pointtmp(3)
812 k1 = k1 + 1
813 ENDDO
814 ENDDO
815 IF(npoint/=0)THEN
816 ppoint(1) = ppoint(1) / (k1+mnod) !NPOINT
817 ppoint(2) = ppoint(2) / (k1+mnod) !NPOINT
818 ppoint(3) = ppoint(3) / (k1+mnod) !NPOINT
819 brick_list(nin,ib)%POLY(9)%CellCENTER(1:3) = ppoint(1:3)
820 ENDIF
821 endif!NCELL/=0
822 enddo!next IB
823
824 !------------------------------------------------------------!
825 ! @11. BUILD UNCUT CELL IN CUT CELL BUFFER !
826 !------------------------------------------------------------!
827 DO ib=nbf,nbl
828 IF(brick_list(nin,ib)%NBCUT == 0)THEN
829 psubvol(1:9) => brick_list(nin,ib)%POLY(1:9)
830 psubvol(1:9)%Vnew = 0
831 gbuf => elbuf_tab(ngb(ib))%GBUF
832 psubvol(1)%Vnew = gbuf%VOL(idlocb(ib))
833 brick_list(nin,ib)%Vnew_SCell = gbuf%VOL(idlocb(ib))
834 ENDIF
835 ENDDO
836
837 !------------------------------------------------------------!
838 ! @12. STORING BRICK REFERENCE VOLUME (UNCUT) !
839 !------------------------------------------------------------!
840 DO ib=nbf,nbl
841 brick_list(nin,ib)%UncutVol = elbuf_tab(ngb(ib))%GBUF%VOL(idlocb(ib))
842 ENDDO
843
844 9 CONTINUE
845
846
847 !------------------------------------------------------------!
848 ! @13. POLY 9 : MARK IF NO NODE ON A FACE !
849 !------------------------------------------------------------!
850 !Poly9woNodes
851 DO ib=nbf,nbl
852 ncell = brick_list(nin,ib)%NBCUT
853 brick_list(nin,ib)%Poly9woNodes(1:6,1:2) = 0
854 IF(ncell<= 1) cycle
855 DO j=1,6
856 face = brick_list(nin,ib)%POLY(9)%FACE(j)%Surf
857 IF(face ==zero) cycle
858 !loop on nodes (1,2,3,4,5,6,7,8) composing face J
859 lfound = .false.
860 DO k=1,4
861 inod = int22_buf%nodFACE(j,k)
862 icell = brick_list(nin,ib)%NODE(inod)%WhichCell
863 IF(icell == 9) lfound = .true.
864 ENDDO
865 IF(lfound)cycle ! CELL 9 IS SHARING A NODE SO ITS ADJACENT CELL WAS FOUND ABOVE (See 14. FINDING & STORING ADJACENT CELLS)
866 !Here Cell 9 has a trace on face J but has no nodes on this face. An adjacent celle must be find (otherwise internal forces issue, missing area)
867 brick_list(nin,ib)%Poly9woNodes(j,1) = 1
868 enddo!next J
869 !UPDATE FLUX FOR THIS CASE : NO NODES => NO FLUXES (CLOSED SURFACE HYPOTHESIS)
870 enddo! next IB
871
872
873
874 DO ib=nbf,nbl
875 brick_list(nin,ib)%IsMergeTarget = 0
876 !print *, "0, ismergetarget", ixs(11,brick_list(nin,ib)%id)
877 ENDDO
878
879
880 !------------------------------------------------------------!
881 ! @14. FILL MATERIAL BUFFER : %bakMAT%... !
882 !------------------------------------------------------------!
883 DO ib=nbf,nbl
884 ncell = brick_list(nin,ib)%NBCUT
885 mcell = brick_list(nin,ib)%MainID
886 id = brick_list(nin,ib)%ID
887 ng = brick_list(nin,ib)%NG
888 idloc = brick_list(nin,ib)%IDLOC
889 mlw = brick_list(nin,ib)%MLW
890 gbuf => elbuf_tab(ng)%GBUF
891 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
892 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
893 llt_ = iparg(2,ng)
894 icell = 0
895 IF(mcell==0)THEN
896 brick_list(nin,ib)%bakMAT%RHO = zero
897 brick_list(nin,ib)%bakMAT%rhoE = zero
898 brick_list(nin,ib)%bakMAT%rhoU(1) = zero
899 brick_list(nin,ib)%bakMAT%rhoU(2) = zero
900 brick_list(nin,ib)%bakMAT%rhoU(3) = zero
901 brick_list(nin,ib)%bakMAT%ssp = zero
902 brick_list(nin,ib)%bakMAT%SIG(1) = zero
903 brick_list(nin,ib)%bakMAT%SIG(2) = zero
904 brick_list(nin,ib)%bakMAT%SIG(3) = zero
905 brick_list(nin,ib)%bakMAT%SIG(4) = zero
906 brick_list(nin,ib)%bakMAT%SIG(5) = zero
907 brick_list(nin,ib)%bakMAT%SIG(6) = zero
908 !------MULTI-MATERIAL-----------!
909 IF(i22law37+i22law51 == 0) cycle
910 brick_list(nin,ib)%bakMAT%UVAR(1) = zero
911 brick_list(nin,ib)%bakMAT%UVAR(2) = zero
912 brick_list(nin,ib)%bakMAT%UVAR(3) = zero
913 brick_list(nin,ib)%bakMAT%UVAR(4) = zero
914 brick_list(nin,ib)%bakMAT%UVAR(5) = zero
915 IF(i22law51 == 0) cycle
916 DO k=6, i22law51
917 brick_list(nin,ib)%bakMAT%UVAR(k) = zero
918 enddo!next K
919 ELSE
920 brick_list(nin,ib)%bakMAT%RHO = gbuf%RHO(idloc)
921 brick_list(nin,ib)%bakMAT%rhoE = gbuf%EINT(idloc)
922 brick_list(nin,ib)%bakMAT%rhoU(1) = gbuf%MOM(llt_*(1-1) +idloc)
923 brick_list(nin,ib)%bakMAT%rhoU(2) = gbuf%MOM(llt_*(2-1) +idloc)
924 brick_list(nin,ib)%bakMAT%rhoU(3) = gbuf%MOM(llt_*(3-1) +idloc)
925 brick_list(nin,ib)%bakMAT%ssp = lbuf%SSP(idloc)
926 brick_list(nin,ib)%bakMAT%SIG(1) = gbuf%SIG(llt_*(1-1) +idloc)
927 brick_list(nin,ib)%bakMAT%SIG(2) = gbuf%SIG(llt_*(2-1) +idloc)
928 brick_list(nin,ib)%bakMAT%SIG(3) = gbuf%SIG(llt_*(3-1) +idloc)
929 brick_list(nin,ib)%bakMAT%SIG(4) = gbuf%SIG(llt_*(4-1) +idloc)
930 brick_list(nin,ib)%bakMAT%SIG(5) = gbuf%SIG(llt_*(5-1) +idloc)
931 brick_list(nin,ib)%bakMAT%SIG(6) = gbuf%SIG(llt_*(6-1) +idloc)
932 !------MULTI-MATERIAL-----------!
933 !IF(I22LAW37+I22LAW51 == 0) CYCLE
934 IF(mlw/=37 .AND. mlw/=51)cycle
935 brick_list(nin,ib)%bakMAT%UVAR(1) = mbuf%VAR((1-1)*llt_+idloc)
936 brick_list(nin,ib)%bakMAT%UVAR(2) = mbuf%VAR((2-1)*llt_+idloc)
937 brick_list(nin,ib)%bakMAT%UVAR(3) = mbuf%VAR((3-1)*llt_+idloc)
938 brick_list(nin,ib)%bakMAT%UVAR(4) = mbuf%VAR((4-1)*llt_+idloc)
939 brick_list(nin,ib)%bakMAT%UVAR(5) = mbuf%VAR((5-1)*llt_+idloc)
940 IF(i22law51 == 0) cycle
941 DO k=6, i22law51
942 brick_list(nin,ib)%bakMAT%UVAR(k) = mbuf%VAR((k-1)*llt_+idloc)
943 enddo!next K
944 ENDIF
945 enddo!next IB
946
947 !-------------------!
948 CALL my_barrier()
949 !-------------------!
950
951 !------------------------------------------------------------!
952 ! @15. FINDING & STORING ADJACENT CELLS !
953 !------------------------------------------------------------!
954 DO ib=nbf,nbl
955 pface(1:9) => brick_list(nin,ib)%POLY(1:9)
956! pFACE(2)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(2)%FACE(1:6)%Surf
957! pFACE(3)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(3)%FACE(1:6)%Surf
958! pFACE(4)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(4)%FACE(1:6)%Surf
959! pFACE(5)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(5)%FACE(1:6)%Surf
960! pFACE(6)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(6)%FACE(1:6)%Surf
961! pFACE(7)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(7)%FACE(1:6)%Surf
962! pFACE(8)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(8)%FACE(1:6)%Surf
963! pFACE(9)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(9)%FACE(1:6)%Surf
964 !pNAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%NAdjCell
965 padjbrick => brick_list(nin,ib)%Adjacent_Brick(1:6,1:5)
966 plistnodid(1)%p(1:8) => brick_list(nin,ib)%POLY(1)%ListNodID(1:8)
967 plistnodid(2)%p(1:8) => brick_list(nin,ib)%POLY(2)%ListNodID(1:8)
968 plistnodid(3)%p(1:8) => brick_list(nin,ib)%POLY(3)%ListNodID(1:8)
969 plistnodid(4)%p(1:8) => brick_list(nin,ib)%POLY(4)%ListNodID(1:8)
970 plistnodid(5)%p(1:8) => brick_list(nin,ib)%POLY(5)%ListNodID(1:8)
971 plistnodid(6)%p(1:8) => brick_list(nin,ib)%POLY(6)%ListNodID(1:8)
972 plistnodid(7)%p(1:8) => brick_list(nin,ib)%POLY(7)%ListNodID(1:8)
973 plistnodid(8)%p(1:8) => brick_list(nin,ib)%POLY(8)%ListNodID(1:8)
974 plistnodid(9)%p(1:8) => brick_list(nin,ib)%POLY(9)%ListNodID(1:8)
975 id = brick_list(nin,ib)%ID
976 ncell = brick_list(nin,ib)%NBCUT
977 !pNAdjCELL(1:9,1:6) = 0
978 DO k=1,6
979 brick_list(nin,ib)%POLY(1:9)%FACE(k)%Adjacent_Cell(1) = 0
980 brick_list(nin,ib)%POLY(1:9)%FACE(k)%Adjacent_Cell(2) = 0
981 brick_list(nin,ib)%POLY(1:9)%FACE(k)%Adjacent_Cell(3) = 0
982 brick_list(nin,ib)%POLY(1:9)%FACE(k)%Adjacent_Cell(4) = 0
983 brick_list(nin,ib)%POLY(1:9)%FACE(k)%Adjacent_Cell(5) = 0
984 brick_list(nin,ib)%POLY(1:9)%FACE(k)%NAdjCell = 0
985 ENDDO
986 icell = 0
987 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
988 icell = icell +1
989 IF (icell>ncell .AND. ncell/=0)icell=9
990 DO j=1,nv46 ! loop on brick face
991 IF(pface(icell)%FACE(j)%Surf>zero)THEN ! if polyhedron has a non zero subface on this brick face
992 iv = padjbrick(j,1)
993 IF(iv==0)cycle !if there is no adjacent brick on this face
994 !if there is an adjacent brick on this face
995 iad22 = padjbrick(j,4)
996 IF(iad22 == 0)THEN
997 !adjacent brick is in cut cell buffer : it also uncut
998 brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell = 1
999 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(1) = 1 !necessarily ICELLv=1
1000 ELSEIF(brick_list(nin,iad22)%NBCUT==0)THEN
1001 !adjacent brick is in cut cell buffer but is uncut
1002 brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell = 1
1003 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(1) = 1 !necessarily ICELLv=1
1004 ELSE
1005
1006 pwhichcellnodv(1:8) => brick_list(nin,iad22)%NODE(1:8)
1007 !---brick nodes used by reference cell, must lay on current face J, and are identified by internal id : IXS(1+.,ID)
1008 nnodes = 0
1009 DO k1=1,brick_list(nin,ib)%POLY(icell)%NumNOD !loop on polyhedron nodes !TODO : might be optimized by looping on nodFACE and cycle with pNNODCELL
1010 DO k2=1,4 !loop on face nodes
1011 IF(plistnodid(icell)%p(k1) == int22_buf%nodFACE(j,k2))THEN !polyhedron node on the current face
1012 nnodes = nnodes +1
1013 inodes(nnodes) = ixs(1+int22_buf%nodFACE(j,k2), id)
1014 ENDIF
1015 ENDDO
1016 ENDDO
1017 !---search for local nodes id in adjacent brick.
1018 icelltag(1:9) = 0
1019 jv = padjbrick(j,5)
1020 DO k0=1,4 !loop on adjacent 'local' nodes(ii,j)<->(iv,jv)
1021 k1 = int22_buf%nodFACE(jv,k0)
1022 DO in=1,nnodes
1023 IF(ixs(1+k1,iv)==inodes(in) )THEN
1024 icellv = pwhichcellnodv(k1)%WhichCell
1025 if (icellv==0)then
1026 print *, "ITASK,ICELLv,pWhichCellNodv(K1)",itask,icellv,pwhichcellnodv(k1)%WhichCell
1027 stop 2245
1028 endif
1029 IF(icelltag(icellv)==0)THEN !check if cell is already stored as an adjacent cell, if not store it
1030 icelltag( icellv ) = 1
1031 brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell = brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell + 1
1032 nadjcell = brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell
1033 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(nadjcell) = icellv
1034 ENDIF
1035 ENDIF
1036 enddo!next IN (polyhedron node on face)
1037 enddo!next K0 (adjacent 'local' node)
1038 endif!testing IV & IBV
1039 endif!(pFACE(ICELL)%p(J)>ZERO)
1040
1041 !COMPLEMENTARY POLYHEDRON WITH NO NODES ON A GIVEN FACE
1042 poly9wonodes = brick_list(nin,ib)%Poly9woNodes(j,1)
1043 ! necessarily IBV /= 0
1044 IF(poly9wonodes == 1)THEN
1045 iv = brick_list(nin,ib)%Adjacent_Brick(j,1)
1046 ibv = brick_list(nin,ib)%Adjacent_Brick(j,4)
1047 IF(ibv==0)THEN
1048 brick_list(nin,ib)%Poly9woNodes(j,2) = 0
1049 cycle
1050 ENDIF
1051 jv = brick_list(nin,ib)%Adjacent_Brick(j,5)
1052 poly9wonodesv = brick_list(nin,ibv)%Poly9woNodes(jv,1)
1053 IF(poly9wonodesv==0)THEN
1054 nbcutv = brick_list(nin,ibv)%NBCUT
1055 IF(nbcutv == 0)THEN
1056 !PAS DE VOISIN HYPOTHESE SURFACE FERMEE
1057 !the neighboring polyhedron is polyhedron 9
1058 ! +--------+--+-----------+
1059 ! | \1| !
1060 ! | | |
1061 ! | 9 | 1 |
1062 ! | | |
1063 ! | / | |
1064 ! | / 2 | |
1065 ! +------+----+-----------+
1066 !BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%NAdjCell = BRICK_LIST(NIN,IB)%PoLY(9)%FACE(J)%NAdjCell + 1
1067 !NAdjCell = BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%NAdjCell
1068 !BRICK_LIST(NIN,IB)%Adjacent_Brick(J,9,NadjCell) = 1
1069 brick_list(nin,ib)%Poly9woNodes(j,2) = 0
1070C !DESTRUCTION POLYEDRE VOISIN DETECTE AU DEBUT DU PARAGRAPHE @14.
1071C pNAdjCell(ICELL,J) = 0
1072C pAdjCELL(J,ICELL,1) = 0 !necessarily ICELLv=1
1073C No if not Cell Trapped Error and also we cannot find a neighboring pressure poru contact
1074 ELSE
1075 !the neighboring polyhedron uses the intersection point closest to the vertex node => the neighboring polyhedron systematically becomes polyhedron 9. (drawings required)
1076 ! +--------+--+-----------+ +--------+--+-----------+
1077 ! | \1| ! | \1| 1 !
1078 ! | + ! | +-----------+
1079 ! | 9 | 9 | | 9 | 9 |
1080 ! | +-----------+ | | |
1081 ! | / | 1 | | / | |
1082 ! | / 2 | | | / 2 | |
1083 ! +------+----+-----------+ +------+----+-----------+
1084 !autre hypothese => ... = 1
1085 brick_list(nin,ib)%POLY(9)%FACE(j)%NAdjCell = brick_list(nin,ib)%POLY(9)%FACE(j)%NAdjCell + 1
1086 nadjcell = brick_list(nin,ib)%POLY(9)%FACE(j)%NAdjCell
1087 brick_list(nin,ib)%POLY(9)%FACE(j)%Adjacent_Cell(nadjcell) = 9
1088 brick_list(nin,ib)%Poly9woNodes(j,2) = 9
1089 ENDIF
1090 ELSEIF(poly9wonodesv/=0)THEN
1091 !the neighboring polyhedron is polyhedron 9
1092 ! +--------+--+-----------+
1093 ! | \1| 2 !
1094 ! | +-----------+
1095 ! | 9 | 9 |
1096 ! | +-----------+
1097 ! | / | 1 |
1098 ! | / 2 | |
1099 ! +------+----+-----------+
1100
1101 ! SIMPLY CONNECTING poly_9 <-> polyV_9
1102 brick_list(nin,ib)%POLY(9)%FACE(j)%NAdjCell = brick_list(nin,ib)%POLY(9)%FACE(j)%NAdjCell + 1
1103 nadjcell = brick_list(nin,ib)%POLY(9)%FACE(j)%NAdjCell
1104 brick_list(nin,ib)%POLY(9)%FACE(j)%Adjacent_Cell(nadjcell) = 9
1105 brick_list(nin,ib)%Poly9woNodes(j,2) = 9
1106 endif!(Poly9woNodesV==0)
1107 ENDIF
1108 enddo!next J (face)
1109 enddo!next ICELL
1110 enddo!next IB
1111
1112
1113 !------------------------------------------------------------!
1114 ! @16. BUILD POLYGONAL TRACE FOR POLY9 !
1115 !------------------------------------------------------------!
1116 !CUT CELL BUFFER REMINDER (check common source for complete & update source)
1117 !my_real :: N(3) !normal vector
1118 !my_real :: B(3) !basis
1119 !my_real :: SCUT(1) !Cut Area
1120 !my_real :: SEFF(1) !Efficace Section = Wetted Surface : Cut Aera - Holes (computed by suming lagrangian projection)
1121 !my_real :: P(3,8) !Cut Surface Polygon (coordinates)
1122 !INTEGER :: NP !Cut Surface Polygon (number of points) ! NP<=6 | NP<=8 for additional Scut (closed surface hypothesis)
1123 !BRICK_LIST(NIN,IB)%Edge(J)%NODE(1:2) = EDGE_LIST(NIN,K)%NODE(1:2)
1124 !BRICK_LIST(NIN,IB)%Edge(J)%NBCUT = EDGE_LIST(NIN,K)%NBCUT
1125 !BRICK_LIST(NIN,IB)%Edge(J)%CUTSHELL(1:2) = EDGE_LIST(NIN,K)%CUTSHELL(1:2)
1126 !BRICK_LIST(NIN,IB)%Edge(J)%CUTCOOR(1:2) = EDGE_LIST(NIN,K)%CUTCOOR(1:2)
1127 !BRICK_LIST(NIN,IB)%Edge(J)%CUTPOINT(1:3,1) = EDGE_LIST(NIN,K)%CUTPOINT(1:3,1)
1128 !BRICK_LIST(NIN,IB)%Edge(J)%CUTPOINT(1:3,2) = EDGE_LIST(NIN,K)%CUTPOINT(1:3,2)
1129 !BRICK_LIST(NIN,IB)%Edge(J)%CUTVEL(1:3,1) = EDGE_LIST(NIN,K)%CUTVEL(1:3,1)
1130 !BRICK_LIST(NIN,IB)%Edge(J)%CUTVEL(1:3,2) = EDGE_LIST(NIN,K)%CUTVEL(1:3,2)
1131 !BRICK_LIST(NIN,IB)%Edge(J)%VECTOR(1:3) = EDGE_LIST(NIN,K)%VECTOR(1:3)
1132 !BRICK_LIST(NIN,IB)%Edge(J)%LEN = EDGE_LIST(NIN,K)%LEN
1133 !
1134 !
1135 ! +--------+--+ +-------+---+
1136 ! | \1| | | \ |
1137 ! | + Check WhichCellNode for each node of the face | | +
1138 ! | 9 | if digits 1,2, and 9 tagged then a polygon has => | | P |
1139 ! | + to be computed | | +
1140 ! | / | | | /|
1141 ! | / 2 | | |/ |
1142 ! +------+----+ +------+----+
1143 !
1144 !
1145 !must be counterclockwise from outward ( clockwise from point inside the brick )
1146 DO ib=nbf,nbl
1147 ie = brick_list(nin,ib)%ID
1148 brick_list(nin,ib)%ClosedSurf(1:6) = 0
1149 DO j=1,6
1150 icode = 0
1151 iv = brick_list(nin,ib)%Adjacent_Brick(j,1)
1152 IF(iv==0)cycle
1153 DO k=1,4
1154 inod = int22_buf%nodFACE(j,k)
1155 icell = brick_list(nin,ib)%NODE(inod)%WhichCell
1156 icode = ibset(icode,icell)
1157 ENDDO
1158
1159 !THIS IS USED TO MARK CLOSURE CASES (CLOSURE HYPOTHESIS)
1160 IF(icode == 518 .OR. icode == 6) THEN
1161 face9 = brick_list(nin,ib)%POLY(9)%FACE(j)%Surf
1162 ibv = brick_list(nin,ib)%Adjacent_Brick(j,4)
1163 !necessarily IBV>0 otherwise we cannot have ICODE=518 (in other words: brick IB is intersected so the neighbor is in the buffer)
1164 nbcutv = brick_list(nin,ibv)%NBCUT
1165 IF(nbcutv == 0 .AND. face9>zero) THEN
1166 !IMPORTANT DISTINCTION BECAUSE WE NEED THE NEIGHBOR FOR CONTACT FORCES BUT OTHERWISE IT MUST BE IGNORED FOR FINT AND CONVECTION
1167 brick_list(nin,ib)%ClosedSurf(j) = 1
1168C print *, "brick_id=",ixs(11,iv)," face=",J, " has a closed surface"
1169C print *, "sinit22_fvm marked closure hypothesis."
1170 ENDIF
1171 ENDIF
1172
1173 poly9wonodes = brick_list(nin,ib)%Poly9woNodes(j,1)
1174 icellv = brick_list(nin,ib)%Poly9woNodes(j,2)
1175 IF(poly9wonodes /=0 .AND. icellv /= 0)cycle
1176 IF(poly9wonodes==0 .AND. icode /= 518)cycle
1177 ! poly 1,2,9 on this current face => ICODE = b1000000110 = 518 (bits 1,2,9)
1178 np = 0
1179 ppoint(1:3) = zero
1180 cut_vel(1:3) = zero
1181 nbcutprev = 0
1182 DO inod=4,1,-1
1183 !downside
1184 pointtmp(1:3) = x(1:3,ixs(1+int22_buf%NodFace(j,inod),ie))
1185 ppoint(1) = ppoint(1) + pointtmp(1)
1186 ppoint(2) = ppoint(2) + pointtmp(2)
1187 ppoint(3) = ppoint(3) + pointtmp(3)
1188 !polygon point
1189 iedg = int22_buf%iLeftEdge(j,inod)
1190 isgn(1:2) = (/1,2/)
1191 IF(iedg < 0) isgn(1:2) = (/2,1/)
1192 iedg = iabs(iedg)
1193 nbcut = brick_list(nin,ib)%Edge(iedg)%NBCUT
1194 IF(nbcut == 0)THEN
1195 cycle
1196 ELSEIF(nbcut == 1)THEN
1197 np = np + 1
1198 cut_point(1:3) = brick_list(nin,ib)%Edge(iedg)%CUTPOINT(1:3,1)
1199 brick_list(nin,ib)%PCUT(8+j)%P(1:3,np) = cut_point(1:3)
1200 cut_vel(1:3) = cut_vel(1:3) + brick_list(nin,ib)%Edge(iedg)%CUTVEL(1:3,1)
1201! IF(NBCUTprev==1)THEN
1202! ! Add the top point to polygon 9
1203! NP = NP + 1
1204! NP _NODE = NP_NODE + 1
1205! CUT_Point(1:3) = X(1:3,IXS(1+iEDGE(ISGN(1),iEDG) ,IE))
1206! brick_list(nin,ib)%PCUT(8+j)%P(1:3,np) = cut_point(1:3)
1207! !do not increment CUTVEL
1208! print *, "pause-22- "
1209! pause
1210! ENDIF
1211 ELSE
1212 np = np + 1
1213 cut_point(1:3) = brick_list(nin,ib)%Edge(iedg)%CUTPOINT(1:3,isgn(1))
1214 brick_list(nin,ib)%PCUT(8+j)%P(1:3,np) = cut_point(1:3)
1215 cut_vel(1:3) = cut_vel(1:3) + brick_list(nin,ib)%Edge(iedg)%CUTVEL(1:3,isgn(1))
1216 np = np + 1
1217 cut_point(1:3) = brick_list(nin,ib)%Edge(iedg)%CUTPOINT(1:3,isgn(2))
1218 cut_vel(1:3) = cut_vel(1:3) + brick_list(nin,ib)%Edge(iedg)%CUTVEL(1:3,isgn(2))
1219 brick_list(nin,ib)%PCUT(8+j)%P(1:3,np) = cut_point(1:3)
1220 ENDIF
1221 nbcutprev = nbcut
1222 enddo!next INOD
1223 brick_list(nin,ib)%PCUT(8+j)%NP = np
1224 brick_list(nin,ib)%PCUT(8+j)%B(1) = fourth * ppoint(1)
1225 brick_list(nin,ib)%PCUT(8+j)%B(2) = fourth * ppoint(2)
1226 brick_list(nin,ib)%PCUT(8+j)%B(3) = fourth * ppoint(3)
1227 brick_list(nin,ib)%PCUT(8+j)%N(1:3) = brick_list(nin,ib)%N(j,1:3) * brick_list(nin,ib)%POLY(9)%FACE(j)%Surf
1228 brick_list(nin,ib)%PCUT(8+j)%SEFF = zero
1229 !BRICK_LIST(NIN,IB)%POLY(9)%FACE0%Vel(1:3) = CUT_Vel(1:3) / (NP*ONE)
1230 brick_list(nin,ib)%PCUT(8+j)%Vel(1:3) = cut_vel(1:3) / (np*one)
1231 ppoint(1:3) = brick_list(nin,ib)%PCUT(8+j)%B(1:3)
1232 vectmp(1:3) = i22aera(np, brick_list(nin,ib)%PCUT(8+j)%P(1:3,1:np), ppoint(1:3) )
1233 face = sqrt(sum(vectmp(1:3)*vectmp(1:3)))
1234 brick_list(nin,ib)%PCUT(8+j)%SCUT = face
1235 enddo!next J
1236 enddo!next IB
1237
1238
1239 !-------------------!
1240 CALL my_barrier()
1241 !-------------------!
1242
1243 !------------------------------------------------------------!
1244 ! @17. PRE-TREATMENT FOR DBLE SECND WITH BOTH V=0 !
1245 !------------------------------------------------------------!
1246 nbl1 = nbl
1247 IF(i22_degenerated == 1)THEN
1248 ALLOCATE (destroy_data(7,2*nb)) !on each thread
1249 ELSE
1250 nbl1 = 0
1251 ENDIF
1252 !-------------------!
1253 CALL my_barrier()
1254 !-------------------!
1255 ipos = 0
1256 DO ib=nbf,nbl1
1257 ncell = brick_list(nin,ib)%NBCUT
1258 icell = 0
1259 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
1260 icell = icell +1
1261 IF (icell>ncell .AND. ncell/=0)icell=9
1262 vol = brick_list(nin,ib)%POLY(icell)%Vnew
1263 uncutvol = brick_list(nin,ib)%UncutVol
1264 ratio = vol / uncutvol
1265 IF(abs(ratio)>em04)cycle !Vol=1.01e-05, UNCUTVOL=10 => ratio > EM06. So EM04 is retained. THRESHOLD VOLUME : Smax = O(V**(2/3)), dh= EM03 * O(V**(1/3)) => Veps = EM04*O(V) => ratio=Veps/V = EM04
1266 idloc = maxloc(brick_list(nin,ib)%POLY(icell)%FACE(1:6)%SURF,1)
1267 ratiof = brick_list(nin,ib)%POLY(icell)%FACE(idloc)%SURF / brick_list(nin,ib)%Face_Brick(idloc)
1268 !volume degenere
1269 !We are looking for a neighboring correspondence
1270 j = idloc
1271 nadjcell = brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell
1272 DO iadj = 1,nadjcell
1273 icellv = brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(iadj)
1274 ibv = brick_list(nin,ib)%Adjacent_Brick(j,4)
1275 jv = brick_list(nin,ib)%Adjacent_Brick(j,5)
1276 itaskv = brick_list(nin,ibv)%ITASK
1277 volv = brick_list(nin,ibv)%POLY(icellv)%Vnew
1278 uncutvolv = brick_list(nin,ibv)%UncutVOL
1279 ratiov = volv/uncutvolv
1280 !vol=1.01e-05, uncutvol=10 => ratio > em06. so em04 is retained. threshold volume : smax = o(v**(2/3)), dh= em03 * o(v**(1/3)) => veps = em04*o(v) => ratio=veps/v = em04
1281 IF(abs(ratiov)>em04 .OR. ibv>ib )cycle
1282 !IF( (ITASKv==ITASK.AND.IB>IBv) .OR. (ITASKv<ITASK) )THEN
1283 ipos = ipos + 1
1284 destroy_data(1,ipos) = nin
1285 destroy_data(2,ipos) = ib
1286 destroy_data(3,ipos) = icell
1287 destroy_data(4,ipos) = icellv
1288 destroy_data(5,ipos) = ibv
1289 destroy_data(6,ipos) = j
1290 destroy_data(7,ipos) = jv
1291 ! ENDIF
1292 enddo!next ICV
1293 !J was IDLOC
1294 enddo!next I
1295 enddo!next IB
1296
1297 !detach from the previous loop and launch after barrier to avoid multi-threading issue
1298 if(itask==0 .AND. ibug22_destroy==1 .and. ipos>0)then
1299 print *, ""
1300 print *, " |------i22_destroy_cell.F-------|"
1301 print *, " | INITIALIZATION SUBROUTINE |"
1302 print *, " |-------------------------------|"
1303 print *, ""
1304 endif
1305
1306 !-------------------!
1307 CALL my_barrier()
1308 !-------------------!
1309
1310 DO i=1,ipos
1311 ib = destroy_data(2,i)
1312 icell = destroy_data(3,i)
1313 icellv = destroy_data(4,i)
1314 ibv = destroy_data(5,i)
1315 j = destroy_data(6,i)
1316 jv = destroy_data(7,i)
1317 !print *, "destroy", ixs(11,brick_list(nin,ib)%id),ICELL
1318 CALL destroy_cell( nin, ib,icell,icellv,ibv,j,jv, ixs, itask)
1319 ENDDO
1320
1321 IF(i22_degenerated == 1)THEN ! this is activated if SECONDARY/SECONDARY adjacency with both 0-volumes
1322 !-------------------!
1323 CALL my_barrier()
1324 !-------------------!
1325
1326 i22loop = i22loop + 1
1327 if(i22loop >= 2)then
1328 print *, "**error : inter22, unexpected situation."
1329 endif
1330 i22_degenerated = 0
1331 GOTO 9 !can be optimized later by updating adjacency in destroy_cell.F or looping only on bricks related to destroyed cell. (need to reach finish line)
1332 ENDIF
1333
1334
1335 !------------------------------------------------------------!
1336 ! @18. TAGGING SECONDARY CUT CELL !
1337 !------------------------------------------------------------!
1338 DO ib=nbf,nbl
1339 ng = brick_list(nin,ib)%NG
1340 pmainid => brick_list(nin,ib)%MainID
1341 pismain(1:9) => brick_list(nin,ib)%POLY(1:9)
1342 pismain(1:9)%IsMain = 0
1343 pmainid = 0
1344 IF(brick_list(nin,ib)%NBCUT==0)THEN
1345 pismain(1)%IsMain = 1
1346 pmainid = 1
1347 cycle
1348 ENDIF
1349 gbuf => elbuf_tab(ngb(ib))%GBUF
1350 vol = brick_list(nin,ib)%UncutVol
1351 ncell = brick_list(nin,ib)%NBCUT
1352 icell = 0
1353 !tag the potential mains
1354 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
1355 icell = icell +1
1356 IF (icell>ncell .AND. ncell/=0)icell=9
1357 !if lower than 50% then it is a SECONDARY.
1358 volcell = brick_list(nin,ib)%POLY(icell)%Vnew
1359 fac = volcell / vol
1360 IF(fac>critmerge22)THEN
1361 !print *, "CritMerge22=", CritMerge22
1362 pismain(icell)%IsMain = 1
1363 ENDIF
1364 enddo!next ICELL
1365 k = sum(pismain(1:9)%IsMain)
1366 !In case of several potential main :Unicity with Reproductibility (indepandant from numbering)
1367 IF(k==1)THEN
1368 pmainid = maxloc(pismain(1:9)%IsMain,1)
1369 ELSEIF(k>=2)THEN
1370 mcell = get_unique_main_cell(nin,ib,k)
1371 pismain(1:9)%IsMain = 0
1372 pismain(mcell)%IsMain = 1
1373 pmainid = mcell
1374 !In case of no cell satisfy criteria, choose cell 9 to avoid problem with partial cut
1375 ELSEIF(k==0)THEN
1376 pismain(9)%IsMain = 1
1377 pmainid = 9
1378 ENDIF
1379 enddo!next IB
1380
1381 !-------------------!
1382 CALL my_barrier()
1383 !-------------------!
1384 !------------------------------------------------------------!
1385 ! @19. LINKING SECONDARY CELL TO A main ONE !
1386 !------------------------------------------------------------!
1387 n_unlinked_l(itask) = 0
1388 !if no main cell is found, it becomes itself a main cell (example : INTERF/INT_22/1BRICK-EOS/CEL_only)
1389 DO ib=nbf,nbl
1390 ie = brick_list(nin,ib)%ID
1391 ncell = brick_list(nin,ib)%NBCUT
1392 padjbrick => brick_list(nin,ib)%Adjacent_Brick(1:6,1:5)
1393 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
1394 !pNAdjCell => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%NAdjCell
1395 pismain(1:9) => brick_list(nin,ib)%POLY(1:9)
1396 pmainid => brick_list(nin,ib)%MainID
1397 mcell = pmainid
1398 icell = 0
1399 brick_list(nin,ib)%POLY(1:9)%WhereIsMain(1) = 0
1400 brick_list(nin,ib)%POLY(1:9)%WhereIsMain(2) = 0
1401 brick_list(nin,ib)%POLY(1:9)%WhereIsMain(3) = 0
1402 brick_list(nin,ib)%POLY(1:9)%WhereIsMain(4) = 0
1403 IF(ncell == 0) THEN
1404 brick_list(nin,ib)%POLY(mcell)%WhereIsMain(3) = ie
1405 brick_list(nin,ib)%POLY(mcell)%WhereIsMain(4) = ib
1406 cycle
1407 ENDIF
1408 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
1409 icell = icell +1
1410 IF (icell>ncell .AND. ncell/=0)icell=9
1411 !--------------------------------------!
1412 !if cell needs to be linked to a main one
1413 IF(icell == mcell)THEN
1414 brick_list(nin,ib)%POLY(icell)%WhereIsMain(1:2) = 0
1415 brick_list(nin,ib)%POLY(icell)%WhereIsMain(3) = ie
1416 brick_list(nin,ib)%POLY(icell)%WhereIsMain(4) = ib
1417 ELSE
1418 adjmain_face(1:6) = zero
1419 adjmain_centroid(1:3,1:6) = zero
1420 idadj_mcellv(1:6) = 0
1421 ivadj_mcellv(1:6) = 0
1422 ibadj_mcellv(1:6) = 0
1423 !Loop on face to list adjacent main cells
1424 IF(sum(brick_list(nin,ib)%POLY(icell)%FACE(1:6)%NAdjCell) == 0)THEN
1425 print *, "**error : inter22 - Cell trapped. Check Lagrangian Surface surrounding BRICK ID:",
1426 . ixs(11,brick_list(nin,ib)%ID) ; stop 2206
1427 cycle
1428 ENDIF
1429 DO j=1,nv46
1430 nadjcell = brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell
1431 IF(nadjcell == 0)cycle
1432 iv = padjbrick(j,1)
1433 ibv = padjbrick(j,4)
1434 IF(ibv==0)THEN
1435 adjmain_face(j) = brick_list(nin,ib)%POLY(icell)%FACE(j)%Surf
1436 idadj_mcellv(j) = 1
1437 ivadj_mcellv(j) = iv
1438 ibadj_mcellv(j) = 0
1439 adjmain_centroid(1,j) = sum(x(1,ixs(2:9,iv))) / eight
1440 adjmain_centroid(2,j) = sum(x(2,ixs(2:9,iv))) / eight
1441 adjmain_centroid(3,j) = sum(x(3,ixs(2:9,iv))) / eight
1442 ELSE
1443 DO k=1,nadjcell
1444 icellv = brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(k)
1445 jv = padjbrick(j,5)
1446 IF(brick_list(nin,ibv)%POLY(icellv)%IsMain == 1)THEN
1447 adjmain_face(j) = min(brick_list(nin,ib)%POLY(icell)%FACE(j)%Surf,
1448 . brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%Surf)
1449 adjmain_centroid(1:3,j) = brick_list(nin,ibv)%POLY(icellv)%CellCenter(1:3)
1450 idadj_mcellv(j) = icellv
1451 ivadj_mcellv(j) = iv
1452 ibadj_mcellv(j) = ibv
1453 EXIT !no more main cell out this face (only one)
1454 ENDIF
1455 enddo!next K
1456 ENDIF
1457 ENDDO
1458 sumface = sum(adjmain_face(1:6))
1459 IF(sumface==zero)THEN
1460c print *, "**Warning inter22 : no main cell to link with cell from BRICK ID:",IXS(11,BRICK_LIST(NIN,IB)%ID)
1461 n_unlinked_l(itask) = n_unlinked_l(itask) + 1
1462 unlinked_cells_l(itask,1,n_unlinked_l(itask)) = ib
1463 unlinked_cells_l(itask,2,n_unlinked_l(itask)) = icell
1464 cycle !number of adjacent cell on given face J : no adjacent cell on face J means no need to search out there
1465 ENDIF
1466! adjMAIN_face(1:6) = adjMAIN_face(1:6) / sumFACE
1467 !IPOS = LINK_WITH_UNIQUE_MAIN_CELL(adjMAIN_face, adjMAIN_centroid)
1468 ipos = maxloc(adjmain_face(1:6),1)
1469 brick_list(nin,ib)%POLY(icell)%WhereIsMain(1) = ipos
1470 brick_list(nin,ib)%POLY(icell)%WhereIsMain(2) = idadj_mcellv(ipos)
1471 brick_list(nin,ib)%POLY(icell)%WhereIsMain(3) = ivadj_mcellv(ipos)
1472 brick_list(nin,ib)%POLY(icell)%WhereIsMain(4) = ibadj_mcellv(ipos)
1473 ENDIF
1474 enddo!next ICELL
1475 ENDDO
1476 !-----------------------------------------------------------!
1477 CALL my_barrier() !
1478 !-----------------------------------------------------------!
1479 if((ibug22_undirectlink==1.OR.ibug22_undirectlink==-1) .AND. n_unlinked_l(itask)>0)
1480 . print *, " **SINIT** : UNLINKED CELL SYNTHESIS "
1481 DO i=1,n_unlinked_l(itask)
1482 ib = unlinked_cells_l(itask,1,i)
1483 icell = unlinked_cells_l(itask,2,i)
1484 if(ibug22_undirectlink==1.OR.ibug22_undirectlink==-1)write(*,fmt='(I10,A1,I1)') ixs(11,brick_list(nin,ib)%ID),".",icell
1485 ENDDO
1486 if(ibug22_undirectlink==1 .AND. n_unlinked_l(itask)>0)print *, ""
1487
1488
1489 !------------------------------------------------------------!
1490 ! @20. LIST ALL SECONDARY CELLS LINKED TO A GIVEN main CELL !
1491 !------------------------------------------------------------!
1492 DO ib=nbf,nbl
1493 nbcut = brick_list(nin,ib)%NBCUT
1494 mcell = brick_list(nin,ib)%MainID
1495 IF(mcell == 0) cycle
1496 !pNAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%NAdjCell
1497 ie = brick_list(nin,ib)%ID
1498 nsecndnod = 0
1499 numsecnd = 0
1500 brick_list(nin,ib)%SecndList%Num = 0
1501 brick_list(nin,ib)%SecndList%VOL_Unmerged = brick_list(nin,ib)%POLY(mcell)%Vnew
1502 DO j=1,nv46
1503 iv = brick_list(nin,ib)%Adjacent_Brick(j,1)
1504 ngv = brick_list(nin,ib)%Adjacent_Brick(j,2)
1505 idlocv = brick_list(nin,ib)%Adjacent_Brick(j,3)
1506 ibv = brick_list(nin,ib)%Adjacent_Brick(j,4) !necessarily different from zero because a cut cell (NBCUT>0) has all its neighbor in cut cell buffer (TZINF increased)
1507 ifv = brick_list(nin,ib)%Adjacent_Brick(j,5)
1508 nadj = brick_list(nin,ib)%POLY(mcell)%FACE(j)%NAdjCell
1509 IF(ibv==0)cycle
1510 DO iadj = 1,nadj
1511 icellv = brick_list(nin,ib)%POLY(mcell)%FACE(j)%Adjacent_Cell(iadj)
1512 ie_m = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(3)
1513 IF(ie_m/=ie)cycle
1514 !secnd cell IBV.ICELLv is added into the list since it is linked to IB.MCELL
1515 numsecnd = numsecnd + 1
1516 nsecndnod = nsecndnod + brick_list(nin,ibv)%POLY(icellv)%NumNOD
1517 brick_list(nin,ib)%SecndList%FM(numsecnd) = j
1518 brick_list(nin,ib)%SecndList%FV(numsecnd) = ifv
1519 brick_list(nin,ib)%SecndList%IV(numsecnd) = iv
1520 brick_list(nin,ib)%SecndList%IBV(numsecnd) = ibv
1521 brick_list(nin,ib)%SecndList%ICELLv(numsecnd) = icellv
1522 brick_list(nin,ib)%SecndList%VOL(numsecnd) = brick_list(nin,ibv)%POLY(icellv)%vnew
1523 brick_list(nin,ib)%SecndList%NumNOD_cell(numsecnd) = brick_list(nin,ibv)%POLY(icellv)%NumNOD
1524 brick_list(nin,ib)%SecndList%ListNodID(numsecnd,1:8) = brick_list(nin,ibv)%POLY(icellv)%ListNodID(1:8)
1525 brick_list(nin,ib)%SecndList%SURF_v(numsecnd) = brick_list(nin,ibv)%POLY(icellv)%FACE(ifv)%SURF
1526 enddo!next IADJ
1527 brick_list(nin,ib)%SecndList%Num = numsecnd
1528 enddo!next J
1529 brick_list(nin,ib)%SecndList%NumSecndNodes = nsecndnod
1530 ENDDO
1531
1532 !------------------------------------------------------------!
1533 ! @21. INDIRECT LINKING !
1534 ! TAG RELATED CELL AND STORE ITS INDIRECT LINK !
1535 !------------------------------------------------------------!
1536 !CREATE THE LINKS IN THE CUT CELL BUFFER
1537 DO iunlink=1,n_unlinked_l(itask)
1538 ib = unlinked_cells_l(itask,1,iunlink)
1539 icell = unlinked_cells_l(itask,2,iunlink)
1540 adjface(:,:) = zero
1541 ipos = 0
1542 DO j=1,nv46
1543 nadj = brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell
1544 ibv = brick_list(nin,ib)%Adjacent_Brick(j,4)
1545
1546 DO iadj=1,nadj
1547 jv = brick_list(nin,ib)%Adjacent_Brick(j,5)
1548 icellv = brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(iadj)
1549 ipos = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(1)
1550 IF(ipos==0)cycle
1551 IF(ibv/=0)THEN
1552 adjface(j,iadj) = min(brick_list(nin,ib)%POLY(icell)%FACE(j)%Surf,
1553 . brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%Surf)
1554 ELSE
1555 adjface(j,iadj) = brick_list(nin,ib)%POLY(icell)%FACE(j)%Surf
1556 ENDIF
1557 ENDDO
1558 ENDDO
1559 ires(1:2) = maxloc(adjface)
1560 iposf = ires(1) !face number
1561 iposiadj = ires(2) !adjacent cell number
1562 ibv = brick_list(nin,ib)%Adjacent_Brick(iposf,4)
1563 icellv = brick_list(nin,ib)%POLY(icell)%FACE(iposf)%Adjacent_Cell(iposiadj)
1564 ipos = 0
1565 IF(icellv/=0)
1566 . ipos = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(1)
1567 IF(ipos==0 .OR. adjface(ires(1),ires(2))==zero)THEN
1568 print *, "***error : inter22 unable to treat cell ",ixs(11,brick_list(nin,ib)%ID)
1569 print *, " Cell seems trapped with no adjacent cells"
1570 stop 2207
1571 ENDIF
1572 IF(ipos>=10)THEN
1573 print *, "***error : inter22 unable to treat cell ",ixs(11,brick_list(nin,ib)%ID)
1574 print *, " Cell seems trapped with no adjacent cells"
1575 print *, " Fluid mesh in interaction area should be refined"
1576 stop 2208
1577 ENDIF
1578 j = iposf*10 + ipos
1579 brick_list(nin,ib)%POLY(icell)%WhereIsMain(1) = j
1580 brick_list(nin,ib)%POLY(icell)%WhereIsMain(2) = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(2)
1581 brick_list(nin,ib)%POLY(icell)%WhereIsMain(3) = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(3)
1582 brick_list(nin,ib)%POLY(icell)%WhereIsMain(4) = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(4)
1583 IF((ibug22_undirectlink==1.OR.ibug22_undirectlink==-1) .AND. n_unlinked_l(itask)>0)THEN
1584 write(*,fmt='(A,I10,A1,I1,A,I2,I2,A1,I2,A2)') "unlinked cell:", ixs(11,brick_list(nin,ib)%ID),
1585 . ".",icell, " is now linked to faces ", iposf, ipos,"(",j," )"
1586 ENDIF
1587 ENDDO
1588
1589 !------------------------------------------------------------!
1590 ! @22. INDIRECT LINKING !
1591 ! UPDATE main CELL BUFFER !
1592 !------------------------------------------------------------!
1593 !need to wait all thread finishing to link before updating their supercell volume and their buffer
1594 !-------------------!
1595 CALL my_barrier()
1596 !-------------------!
1597
1598 !MERGE VOLUMES AND UPDATE MAIN BUFFER (secondary main may be on another thread: OpenMP processing)
1599 DO it=0,nthread-1
1600 DO iunlink=1,n_unlinked_l(it)
1601 ib = unlinked_cells_l(it,1,iunlink)
1602 icell = unlinked_cells_l(it,2,iunlink)
1603 ibm = brick_list(nin,ib)%POLY(icell)%WhereIsMain(4)
1604 itaskb = brick_list(nin,ibm)%ITASK
1605
1606! write (*,FMT='(A,I5,A,I5,I5,A,I5)')
1607! . "traitement ib=", ib, "itask,itaskb=", itask, itaskb," N_UNLINKED_L(IT)",N_UNLINKED_L(IT)
1608 IF(itaskb/=itask)cycle
1609 !Link Unlinked secnd cell to main IB,ICELL
1610 volcell = brick_list(nin,ib)%POLY(icell)%Vnew
1611 id = brick_list(nin,ib)%ID
1612 j = brick_list(nin,ib)%POLY(icell)%WhereIsMain(1)
1613 j1 = j/10
1614 j2 = mod(j,10)
1615 ibv = brick_list(nin,ib)%Adjacent_Brick(j1,4) !cell directly linked to a main
1616 IF(ibv==0)THEN
1617 print *, "inter22 :Error lagrangian surface is escaping eulerian domain."
1618 !with CALL ARRET() : truss might be with len=0 => time step issue.
1619 stop
1620 ENDIF
1621 fm = brick_list(nin,ibv)%Adjacent_Brick(j2,5)
1622 numsecnd = brick_list(nin,ibm)%SecndList%Num
1623 nsecndnod = brick_list(nin,ibm)%SecndList%NumSecndNodes
1624 numsecnd = numsecnd + 1
1625 nsecndnod = nsecndnod + brick_list(nin,ib)%POLY(icell)%NumNOD
1626 brick_list(nin,ibm)%SecndList%Num = numsecnd
1627 brick_list(nin,ibm)%SecndList%NumSecndNodes = nsecndnod
1628 brick_list(nin,ibm)%SecndList%FM(numsecnd) = fm
1629 brick_list(nin,ibm)%SecndList%FV(numsecnd) = j !J>10
1630 brick_list(nin,ibm)%SecndList%IV(numsecnd) = id
1631 brick_list(nin,ibm)%SecndList%IBV(numsecnd) = ib
1632 brick_list(nin,ibm)%SecndList%ICELLv(numsecnd) = icell
1633 brick_list(nin,ibm)%SecndList%VOL(numsecnd) = volcell
1634 brick_list(nin,ibm)%SecndList%NumNOD_cell(numsecnd) = brick_list(nin,ib)%POLY(icell)%NumNOD
1635 brick_list(nin,ibm)%SecndList%ListNodID(numsecnd,1:8) = brick_list(nin,ib)%POLY(icell)%ListNodID(1:8)
1636 brick_list(nin,ibm)%SecndList%SURF_v(numsecnd) = brick_list(nin,ib)%POLY(icell)%FACE(j1)%SURF
1637 enddo!next IUNLINK
1638 enddo!next IT
1639
1640 !increase the main volume
1641 !done later
1642
1643
1644 !------------------------------------------------------------!
1645 ! @23. TAG SECND CELLS WHICH ARE BASED ON FORMER main NODE !
1646 !------------------------------------------------------------!
1647! TARGET_DATA(:,:,:) = 0 !tag main cells linked to these secnd ones.
1648! Ntarget(NBF:NBL) = 0
1649 DO ib=nbf,nbl
1650 brick_list(nin,ib)%Ntarget = 0
1651 nbcut = brick_list(nin,ib)%NBCUT
1652 ncell = nbcut
1653 icell = 0
1654 pnodwasmain(1:8) => brick_list(nin,ib)%NODE(1:8)!%NodWasMain
1655 plistnodid(1)%p(1:8) => brick_list(nin,ib)%POLY(1)%ListNodID(1:8)
1656 plistnodid(2)%p(1:8) => brick_list(nin,ib)%POLY(2)%ListNodID(1:8)
1657 plistnodid(3)%p(1:8) => brick_list(nin,ib)%POLY(3)%ListNodID(1:8)
1658 plistnodid(4)%p(1:8) => brick_list(nin,ib)%POLY(4)%ListNodID(1:8)
1659 plistnodid(5)%p(1:8) => brick_list(nin,ib)%POLY(5)%ListNodID(1:8)
1660 plistnodid(6)%p(1:8) => brick_list(nin,ib)%POLY(6)%ListNodID(1:8)
1661 plistnodid(7)%p(1:8) => brick_list(nin,ib)%POLY(7)%ListNodID(1:8)
1662 plistnodid(8)%p(1:8) => brick_list(nin,ib)%POLY(8)%ListNodID(1:8)
1663 plistnodid(9)%p(1:8) => brick_list(nin,ib)%POLY(9)%ListNodID(1:8)
1664 mcell = brick_list(nin,ib)%MainID
1665 ntag = 0
1666 main_tag(:,:) = 0
1667
1668 brick_list(nin,ib)%MergeTarget(1:3,1:5) = 0
1669
1670 inod_w = brick_list(nin,ib)%OldMainStrongNode
1671 IF( inod_w*dt1 > zero)THEN
1672 IF (brick_list(nin,ib)%NODE(inod_w)%WhichCell==mcell)ncell = -1
1673 ENDIF
1674 !Finding secnd cell inside current brick element, which are based on former main cell nodes
1675 !------LOOP--------!
1676 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
1677 icell = icell +1
1678 IF (icell>ncell .AND. ncell/=0)icell=9
1679 IF(icell==mcell)cycle
1680 mnod = brick_list(nin,ib)%POLY(icell)%NumNOD
1681 CALL weighting_cell_nodes(nin,ib,icell,inod_w, idemerge)
1682 IF(sum(pnodwasmain(plistnodid(icell)%p(1:mnod))%NodWasMain)==0)cycle !if none of the polyhedra nodes is based on a former main cell node : nothing to do
1683 IF(pnodwasmain(inod_w)%NodWasMain<=0) cycle !(negative in debug if IB new in the buffer)
1684 !ICELL has at least one node from the previous main cell
1685 ! Lets Get Its Main Cell
1686
1687 !Continuity Criteria To avoid unexpected merging.
1688 !Consistent only with kinematic time step "several cycle to fully cross an element".
1689 vol_s = brick_list(nin,ib)%POLY(icell)%vnew
1690 uncutvol = brick_list(nin,ib)%UncutVol
1691 var = vol_s / (uncutvol)
1692 IF (var < critdvol22)cycle
1693 inod_w = brick_list(nin,ib)%OldMainStrongNode
1694 IF(inod_w==0)cycle !nothing to merge if no main or no cut cell in previous cycle.
1695
1696 ipos = brick_list(nin,ib)%POLY(icell)%WhereIsMain(1)
1697 icv = brick_list(nin,ib)%POLY(icell)%WhereIsMain(2)
1698 ibm = brick_list(nin,ib)%POLY(icell)%WhereIsMain(4)
1699 ntag = ntag + 1 !number of secnd cells wet by former main cell
1700c Target_data (ib, ntag, 1) = iPos! Facet locating the largest main
1701c TARGET_DATA(IB,NTAG,2) = ICV !cell number of this main
1702 brick_list(nin,ib)%MergeTarget(1,ntag) = ipos
1703 brick_list(nin,ib)%MergeTarget(2,ntag) = icv
1704 brick_list(nin,ib)%MergeTarget(3,ntag) = ibm
1705 enddo! next ICELL
1706 !----END-LOOP------!
1707c Ntarget(IB) = NTAG ! number of secnd cell whose nodes were wet by the former main cell.= number of main to do material deportation (if one main per secnd cell)
1708 brick_list(nin,ib)%NTarget = ntag
1709 ! Ntarget always lower or equal to 9 (max number of secnd cell if one main per secnd)
1710 ! enumeration of the target
1711 !Warning : what happen if 2 secnd cells have the same main cell. Duplicates are not possible by tagging main target
1712 enddo!next IB
1713
1714 !------------------------------------------------------------!
1715 ! @24. VOLUME MERGING BETWEEN SECND CELL & ITS main ONE !
1716 ! UNCONDITIONAL !
1717 ! for law51 old volume need also to be merged !
1718 !------------------------------------------------------------!
1719 !loop on intersected bricks.
1720 DO ib=nbf,nbl
1721 ncell = brick_list(nin,ib)%NBCUT
1722 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
1723 !pNAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%NAdjCell
1724 padjbrick => brick_list(nin,ib)%Adjacent_Brick(1:6,1:5)
1725 psubvol(1:9) => brick_list(nin,ib)%POLY(1:9)
1726 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
1727 mcell = brick_list(nin,ib)%MainID
1728 id = brick_list(nin,ib)%ID
1729 ng = brick_list(nin,ib)%NG
1730 idloc = brick_list(nin,ib)%IDLOC
1731 icell = 0
1732 brick_list(nin,ib)%Vnew_SCell = brick_list(nin,ib)%POLY(mcell)%Vnew
1733 !------LOOP--------!
1734 numsecnd = brick_list(nin,ib)%SecndList%Num
1735 DO ic=1,numsecnd
1736 icellv = brick_list(nin,ib)%SecndList%ICELLv(ic)
1737 ibv = brick_list(nin,ib)%SecndList%IBv(ic)
1738 ngv = brick_list(nin,ibv)%NG
1739 idlocv = brick_list(nin,ibv)%IDLOC
1740 iv = brick_list(nin,ibv)%ID
1741 !---MONO-MATERIAL---!
1742 psubvolv(1:9) => brick_list(nin,ibv)%POLY(1:9)
1743 volv = psubvolv(icellv)%Vnew
1744 !pSUBVOL(MCELL) = pSUBVOL(MCELL) + VOLv !no longer cumulated here
1745 brick_list(nin,ib)%Vnew_SCell = brick_list(nin,ib)%Vnew_SCell + volv
1746 enddo!next IC
1747 !initialization of V1OLD
1748 IF(dt1==zero)THEN
1749 !---MULTI-MATERIAL----!
1750 mtn_ = iparg(1,ngb(ib))
1751 IF(mtn_==51)THEN
1752 llt_ = iparg(2,ngb(ib))
1753 !in sigesp51.F:
1754 !V1OLD = V1OLD - TIMESTEP*DDVOL1
1755 !Volume change correspond to DDVOL=DV/DT = f(fluxes). Gird deformation influences relative material velocity ans also fluxes
1756 DO itrimat=1,trimat
1757 ipos = 1
1758 k1 = m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1 ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
1759 k = k1 * llt_
1760 vfrac(itrimat) = mbuf%VAR(k+idlocb(ib)) !backup for material merging
1761 ipos = 11
1762 k1 = m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1 ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
1763 k = k1 * llt_
1764 var = brick_list(nin,ib)%Vnew_SCell*vfrac(itrimat)
1765 mbuf%VAR(k+idlocb(ib)) = var
1766 enddo!next ITRIMAT
1767 endif!(MTN_==51)
1768 endif!TIME==ZERO
1769 enddo!next IB
1770
1771
1772 !------------------------------------------------------------!
1773 ! @25. MATERIAL MERGING (IF main HAS TO BE DEPORTED) !
1774 !------------------------------------------------------------!
1775 ! if Ntarget > 0 then material deporation has to be done from former main cell buffer (GBUF)
1776 ! to the N=Ntarget target main cells.
1777 nbl1 = nbl
1778 IF(dt1==zero)nbl1 = 0
1779 if(ibug22_merge/=0)then
1780 if(itask==0)allocate(tmp22array(7,nb))
1781 call my_barrier
1782 tmp22array(1:7,nbf:nbl1)=zero
1783 endif
1784 DO ib=nbf,nbl1
1785 nbcut = brick_list(nin,ib)%NBCUT
1786 ncell = nbcut
1787 icell = 0
1788 gbuf => elbuf_tab(ngb(ib))%GBUF
1789 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
1790 padjbrick => brick_list(nin,ib)%Adjacent_Brick(1:6,1:5)
1791 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
1792 idloc = brick_list(nin,ib)%IDLOC
1793 vol = zero
1794 volv = zero
1795 ntar = brick_list(nin,ib)%Ntarget
1796 llt_ = iparg(2,ngb(ib))
1797 !------LOOP--------!
1798 !loop on a main target to distribute former main data
1799 DO itar = 1,ntar
1800 j = brick_list(nin,ib)%MergeTarget(1,itar)
1801 icv = brick_list(nin,ib)%MergeTarget(2,itar)
1802 IF(j>=10)THEN
1803 j1 = j/10
1804 j2 = mod(j,10)
1805 ibv_i = brick_list(nin,ib )%Adjacent_Brick(j1,4)
1806 ngv = brick_list(nin,ibv_i)%Adjacent_Brick(j2,2)
1807 idlocv = brick_list(nin,ibv_i)%Adjacent_Brick(j2,3)
1808 ibv = brick_list(nin,ibv_i)%Adjacent_Brick(j2,4)
1809 ELSE
1810 ngv = brick_list(nin,ib )%Adjacent_Brick(j,2)
1811 idlocv = brick_list(nin,ib )%Adjacent_Brick(j,3)
1812 ibv = brick_list(nin,ib )%Adjacent_Brick(j,4) !TARGET CELL DATA IS STORED IN EINT_L(...,IBv), RHOL, UVARL...
1813 ENDIF
1814 gbufv => elbuf_tab(ngv)%GBUF
1815 ratio = one/ntar !in the loop but avoid conditional tests "IF(Ntar>0)" for NBF:NBL if not applicable
1816 itaskv = brick_list(nin,ibv)%ITASK
1817 imergel(itaskv) = 1
1818 brick_list(nin,ibv)%IsMergeTarget = 1
1819 !print *, "set isMergetarget ", ixs(11,brick_list(nin,ibv)%id)
1820 !---debug---!
1821 cod1=ixs(11,brick_list(nin,ib)%id)
1822 cod2=ixs(11,brick_list(nin,ibv)%id)
1823 if (ibug22_merge==-1)then
1824 !!!write(*,FMT='(A,I8,A,I8)') " MERGING : id=", Cod1, " to idv:", Cod2
1825 tmp22array(1,ib)=cod1
1826 tmp22array(2,ib)=cod2
1827 endif
1828 !-----------!
1829 !---------------------!
1830 ! MATERIAL MERGING ! !only if main cell loose all its nodes, main (which became secnd) contents is sent to the adjacent main.
1831 !---------------------!
1832 ! example:
1833 ! EINTv : target main on which secnd cell is linked
1834 ! EINT : old main cell buffer
1835 !-----------------------MONO-MATERIAL------------------------!
1836 vol = brick_list(nin,ib)%Vold_SCell
1837 supercellvol_l(itask,0,ibv) = supercellvol_l(itask,0,ibv) + ratio * vol
1838 !------ENERGIE--------!
1839 eint = ratio * vol * brick_list(nin,ib)%bakMAT%rhoE
1840 eint_l(itask,ibv) = eint_l(itask,ibv) + eint
1841 !------MASS----------!
1842 if (ibug22_merge==-1)then
1843 !!!write(*,FMT='(A,E30.16)') " RATIO=", RATIO
1844 !!!write(*,FMT='(A,E30.16)') " +rho=", GBUF%RHO(IDLOCB(IB))
1845 !!!write(*,FMT='(A,E30.16)') " +rhoUx=",GBUF%MOM(LLT_*(1-1) + IDLOCB(IB))
1846 !!!write(*,FMT='(A,E30.16)') " +Vol=", VOL
1847 tmp22array(3,ib)=ratio
1848 tmp22array(4,ib)=gbuf%RHO(idlocb(ib))
1849 tmp22array(5,ib)=gbuf%MOM(llt_*(1-1) + idlocb(ib))
1850 tmp22array(6,ib)=vol
1851 endif
1852 rho = ratio * vol * brick_list(nin,ib)%bakMAT%RHO
1853 rho_l(itask,ibv) = rho_l(itask,ibv) + rho
1854 !------TENSOR--------!
1855 DO j=1,nv46
1856 sig(j) = ratio * vol * brick_list(nin,ib)%bakMAT%SIG(j)
1857 sig_l(itask,j,ibv) = sig_l(itask,j,ibv) + sig(j)
1858 ENDDO
1859 !----MOMENTUM--------!
1860 !momentum is (rho.Ux, rho.Uy, rho.Uz)
1861 mom(1) = ratio * vol * brick_list(nin,ib)%bakMAT%rhoU(1)
1862 mom_l(1,itask,ibv) = mom_l(1,itask,ibv) + mom(1)
1863 mom(2) = ratio * vol * brick_list(nin,ib)%bakMAT%rhoU(2)
1864 mom_l(2,itask,ibv) = mom_l(2,itask,ibv) + mom(2)
1865 mom(3) = ratio * vol * brick_list(nin,ib)%bakMAT%rhoU(3)
1866 mom_l(3,itask,ibv) = mom_l(3,itask,ibv) + mom(3)
1867 !-----------------------MULTI-MATERIAL-----------------------!
1868 mtn_ = iparg(1,ngb(ib))
1869 IF(mtn_==37)THEN
1870 !UVAR(I,1) : massic percentage of liquid * global density (rho1*V1/V : it needs to give liquid mass multiplying by element volume in aleconve.F)
1871 !UVAR(I,2) : density of gas
1872 !UVAR(I,3) : density of liquid
1873 !UVAR(I,4) : volumetric fraction of liquid
1874 !UVAR(I,5) : volumetric fraction of gas
1875 llt_ = iparg(2,ngb(ib))
1876c print *, "merge-adding ratio uvar5=" , BRICK_LIST(NIN,IB)%bakMAT%UVAR(5)
1877c print *, "merge-adding ratio uvar4=" , BRICK_LIST(NIN,IB)%bakMAT%UVAR(4)
1878c print *, "merge-adding ratio uvar3=" , BRICK_LIST(NIN,IB)%bakMAT%UVAR(3)
1879c print *, "merge-adding ratio uvar2=" , BRICK_LIST(NIN,IB)%bakMAT%UVAR(2)
1880c print *, "merge-adding ratio uvar1=" , BRICK_LIST(NIN,IB)%bakMAT%UVAR(1)
1881
1882 uvarl(itask,ibv,5) = uvarl(itask,ibv,5)+ratio*vol* brick_list(nin,ib)%bakMAT%UVAR(5)
1883 uvarl(itask,ibv,4) = uvarl(itask,ibv,4)+ratio*vol* brick_list(nin,ib)%bakMAT%UVAR(4)
1884 uvarl(itask,ibv,3) = uvarl(itask,ibv,3)+ratio*vol* brick_list(nin,ib)%bakMAT%UVAR(3)*brick_list(nin,ib)%bakMAT%UVAR(4)
1885 uvarl(itask,ibv,2) = uvarl(itask,ibv,2)+ratio*vol* brick_list(nin,ib)%bakMAT%UVAR(2)*brick_list(nin,ib)%bakMAT%UVAR(5)
1886 uvarl(itask,ibv,1) = uvarl(itask,ibv,1)+ratio*vol* brick_list(nin,ib)%bakMAT%UVAR(1)
1887 supercellvol_l(itask,1,ibv) = supercellvol_l(itask,1,ibv) + uvarl(itask,ibv,4)
1888 supercellvol_l(itask,2,ibv) = supercellvol_l(itask,2,ibv) + uvarl(itask,ibv,5)
1889
1890 ELSEIF(mtn_==51)THEN
1891 llt_ = iparg(2,ngb(ib))
1892 !in sigesp51.F:
1893 !V1OLD = V1OLD - TIMESTEP*DDVOL1
1894 !Volume change correspond to DDVOL=DV/DT = f(fluxes). Gird deformation influences relative material velocity ans also fluxes
1895 DO itrimat=1,trimat
1896 ipos = 11
1897 k = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
1898 vol51(ib,itrimat) = brick_list(nin,ib)%bakMAT%UVAR(k) !backup for material merging
1899 var = ratio * vol51(ib,itrimat)
1900 IF(var/=zero) supercellvol_l(itask,itrimat,ibv) = supercellvol_l(itask,itrimat,ibv) + var
1901 enddo!next ITRIMAT
1902 uvarl(itask,ibv,1) = uvarl(itask,ibv,1) + brick_list(nin,ib)%bakMAT%UVAR(1)
1903 uvarl(itask,ibv,2) = uvarl(itask,ibv,2) + brick_list(nin,ib)%bakMAT%UVAR(2)
1904 uvarl(itask,ibv,3) = uvarl(itask,ibv,3) + brick_list(nin,ib)%bakMAT%UVAR(3)
1905c print *, "adding 1 " , BRICK_LIST(NIN,IB)%bakMAT%UVAR(1), UVARL(ITASK,IBv,1)
1906c print *, "adding 2 " , BRICK_LIST(NIN,IB)%bakMAT%UVAR(2), UVARL(ITASK,IBv,2)
1907c print *, "adding 3 " , BRICK_LIST(NIN,IB)%bakMAT%UVAR(3), UVARL(ITASK,IBv,3)
1908 DO itrimat=1,trimat
1909 var = ratio*vol51(ib,itrimat)
1910 IF(var==zero)cycle
1911 DO ipos = 1,m51_nvphas
1912 IF(ipos==11)cycle !volume already treated merged during geomeetrical merging.
1913 k = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
1914 pvar => brick_list(nin,ib)%bakMAT%UVAR(k)
1915 uvarl(itask,ibv,k) = uvarl(itask,ibv,k) + pvar * var
1916c print *, "adding j " , BRICK_LIST(NIN,IB)%bakMAT%UVAR(K) ,VAR, UVARL(ITASK,IBv,K)
1917 enddo!next IPOS
1918 enddo!next ITRIMAT
1919 endif!(MTN_==51)
1920 !------------------------------------------------------------!
1921 enddo! next ITAR
1922 !----END-LOOP------!
1923 enddo!next IB
1924
1925 !--------------!
1926 CALL my_barrier
1927 !--------------!
1928
1929
1930 if(ibug22_merge/=0 .and.dt1>zero)then
1931 if(itask==0)then
1932 do ib=1,nb
1933 cod1=nint(tmp22array(1,ib))
1934 cod2=nint(tmp22array(2,ib))
1935 if(cod1<=0)cycle
1936 write(*,fmt='(A,I8,A,I8)') " MERGING : id=", cod1, " to idv:", cod2
1937 write(*,fmt='(A,E30.16)') " RATIO=", tmp22array(3,ib)
1938 write(*,fmt='(A,E30.16)') " +rho=", tmp22array(4,ib)
1939 write(*,fmt='(A,E30.16)') " +rhoUx=", tmp22array(5,ib)
1940 write(*,fmt='(A,E30.16)') " +Vol=", tmp22array(6,ib)
1941 tmp22array(:,ib) = zero
1942 enddo
1943 endif
1944 call my_barrier !otherwise TMP22ARRAY IS MODIFIED
1945 endif
1946
1947
1948 !OPENMP STACKING
1949
1950! IF(IMERGEL(ITASK)==0)THEN
1951! NBL1=0
1952! ELSE
1953! IF(TRIMAT<0)THEN
1954! ALLOCATE(UVAR (I22LAW37))
1955! ELSE
1956! ALLOCATE(UVAR (M51_N0PHAS+TRIMAT*M51_NVPHAS))
1957! ENDIF
1958! ENDIF
1959
1960 DO ib=nbf,nbl1
1961 nbcut = brick_list(nin,ib)%NBCUT
1962 ncell = nbcut
1963 icell = 0
1964 gbuf => elbuf_tab(ngb(ib))%GBUF
1965 lbuf => elbuf_tab(ngb(ib))%BUFLY(1)%LBUF(1,1,1)
1966 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
1967 padjbrick => brick_list(nin,ib)%Adjacent_Brick(1:6,1:5)
1968 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
1969 mlw = brick_list(nin,ib)%MLW
1970 som = sum(supercellvol_l(0:nthread-1,0,ib))
1971 llt_ = iparg(2,ngb(ib))
1972 IF(mlw==37 .OR. mlw==51)THEN
1973 som_(1) = sum(supercellvol_l(0:nthread-1,1,ib))
1974 som_(2) = sum(supercellvol_l(0:nthread-1,2,ib))
1975 som_(3) = sum(supercellvol_l(0:nthread-1,3,ib))
1976 som_(4) = sum(supercellvol_l(0:nthread-1,4,ib))
1977 ENDIF
1978 IF(som==zero)cycle
1979 ng = brick_list(nin,ib)%NG
1980 idloc = brick_list(nin,ib)%IDLOC
1981 brickid = brick_list(nin,ib)%ID
1982
1983 IF(brick_list(nin,ib)%Ntarget>0)THEN
1984 delta = zero
1985 ELSE
1986 delta = one
1987 ENDIF
1988
1989 if (ibug22_merge==-1)then
1990 !!!write(*,FMT='(A,E30.16)') " TARGET=", IXS(11,brick_list(nin,ib)%id)
1991 !!!write(*,FMT='(A,E30.16)') " DELTA=", DELTA
1992 !!!write(*,FMT='(A,E30.16)') " rho=", BRICK_LIST(NIN,IB)%bakMAT%RHO
1993 !!!write(*,FMT='(A,3E30.16)') " +rhoUx=", BRICK_LIST(NIN,IB)%bakMAT%rhoU(1:3)
1994 !!!write(*,FMT='(A,E30.16)') " +Vol=", BRICK_LIST(NIN,IB)%Vold_SCell
1995 tmp22array(1,ib) = ixs(11,brick_list(nin,ib)%id)
1996 tmp22array(2,ib) = delta
1997 tmp22array(3,ib) = brick_list(nin,ib)%bakMAT%RHO
1998 tmp22array(4:6,ib) = brick_list(nin,ib)%bakMAT%rhoU(1:3)
1999 tmp22array(7,ib) = brick_list(nin,ib)%Vold_SCell
2000 endif
2001
2002 som = som + delta * brick_list(nin,ib)%Vold_SCell
2003 supercellvol_l(0:nthread-1,0,ib) = zero
2004
2005 eint = delta*brick_list(nin,ib)%bakMAT%rhoE*brick_list(nin,ib)%Vold_SCell +sum(eint_l(0:nthread-1,ib))
2006 eint = eint / som
2007 eint_l(0:nthread-1,ib) = zero !to avoid full reinit
2008
2009 rho = delta*brick_list(nin,ib)%bakMAT%RHO*brick_list(nin,ib)%Vold_SCell +sum(rho_l(0:nthread-1,ib))
2010 rho = rho / som
2011 rho_l(0:nthread-1,ib) = zero !to avoid full reinit
2012
2013 !MomX
2014 mom(1) = delta*brick_list(nin,ib)%bakMAT%rhoU(1)*brick_list(nin,ib)%Vold_SCell+sum(mom_l(1,0:nthread-1,ib))
2015 mom(1) = mom(1) / som
2016 mom_l(1,0:nthread-1,ib)= zero !to avoid full reinit
2017 !MomY
2018 mom(2) = delta*brick_list(nin,ib)%bakMAT%rhoU(2)*brick_list(nin,ib)%Vold_SCell+sum(mom_l(2,0:nthread-1,ib))
2019 mom(2) = mom(2) / som
2020 mom_l(2,0:nthread-1,ib)= zero !to avoid full reinit
2021 !MomZ
2022 mom(3) = delta*brick_list(nin,ib)%bakMAT%rhoU(3)*brick_list(nin,ib)%Vold_SCell+sum(mom_l(3,0:nthread-1,ib))
2023 mom(3) = mom(3) / som
2024 mom_l(3,0:nthread-1,ib)= zero !to avoid full reinit
2025
2026 DO j=1,nv46
2027 sig(j) = delta*brick_list(nin,ib)%bakMAT%SIG(j)*brick_list(nin,ib)%Vold_SCell+sum(sig_l(0:nthread-1,j,ib))
2028 sig(j) = sig(j) / som
2029 sig_l(0:nthread-1,j,ib)= zero !to avoid full reinit
2030 ENDDO
2031
2032 gbuf%EINT(idloc) = eint
2033 gbuf%RHO(idloc) = rho ; alefvm_buffer%FCELL(4, brick_list(nin,ib)%ID ) = rho
2034 gbuf%MOM(llt_*(1-1) + idloc) = mom(1) ; alefvm_buffer%FCELL(1, brick_list(nin,ib)%ID ) = mom(1)
2035 gbuf%MOM(llt_*(2-1) + idloc) = mom(2) ; alefvm_buffer%FCELL(2, brick_list(nin,ib)%ID ) = mom(2)
2036 gbuf%MOM(llt_*(3-1) + idloc) = mom(3) ; alefvm_buffer%FCELL(3, brick_list(nin,ib)%ID ) = mom(3)
2037 !Fcell (5,) see conditional block below
2038
2039 DO j=1,nv46
2040 gbuf%SIG(llt_*(j-1)+idloc) = sig(j)
2041 ENDDO
2042
2043 ! Cod1=ixs(11,brick_list(nin,ib)%id)
2044 ! if (IBUG22_merge==Cod1 .or. IBUG22_merge==-1)then
2045 ! write(*,FMT='(A,I8,A,I8)') " AFTER MERGING : id=", Cod1, " is going to be updated"
2046 ! write(*,FMT='(A,E30.16)') " rho=", RHO
2047 ! write(*,FMT='(A,E30.16)') " rho.Ux=", MOM(1)
2048 ! print *,"-------------------------"
2049 ! endif
2050
2051 mtn_ = iparg(1,ngb(ib))
2052 IF(mtn_==37)THEN
2053 !UVAR(I,1) : massic percentage of liquid * global density (rho1*V1/V : it needs to give liquid mass multiplying by element volume in aleconve.F)
2054 !UVAR(I,2) : density of gas
2055 !UVAR(I,3) : density of liquid
2056 !UVAR(I,4) : volumetric fraction of liquid
2057 !UVAR(I,5) : volumetric fraction of gas
2058 llt_ = iparg(2,ngb(ib))
2059
2060 som_(1) = som_(1) + delta*brick_list(nin,ib)%bakMAT%UVAR(4)*brick_list(nin,ib)%Vold_SCell
2061 som_(2) = som_(2) + delta*brick_list(nin,ib)%bakMAT%UVAR(5)*brick_list(nin,ib)%Vold_SCell
2062
2063
2064 uvar(5) = delta*brick_list(nin,ib)%bakMAT%UVAR(5)*brick_list(nin,ib)%Vold_SCell + sum(uvarl(0:nthread-1,ib,5))
2065 uvar(4) = delta*brick_list(nin,ib)%bakMAT%UVAR(4)*brick_list(nin,ib)%Vold_SCell + sum(uvarl(0:nthread-1,ib,4))
2066 uvar(3) = delta*brick_list(nin,ib)%bakMAT%UVAR(3)*brick_list(nin,ib)%bakMAT%UVAR(4)*brick_list(nin,ib)%Vold_SCell
2067 uvar(3) = uvar(3) + sum(uvarl(0:nthread-1,ib,3))
2068 uvar(2) = delta*brick_list(nin,ib)%bakMAT%UVAR(2)*brick_list(nin,ib)%bakMAT%UVAR(5)*brick_list(nin,ib)%Vold_SCell
2069 uvar(2) = uvar(2) + sum(uvarl(0:nthread-1,ib,2))
2070 uvar(1) = delta*brick_list(nin,ib)%bakMAT%UVAR(1)*brick_list(nin,ib)%Vold_SCell + sum(uvarl(0:nthread-1,ib,1))
2071 uvarl(0:nthread-1,ib,1:5) = zero
2072
2073c print *, "merge uvar-5", BRICK_LIST(NIN,IB)%bakMAT%UVAR(5),UVAR(1) / SOM
2074c print *, "merge uvar-4", BRICK_LIST(NIN,IB)%bakMAT%UVAR(4),UVAR(2) / SOM
2075c print *, "merge uvar-3", BRICK_LIST(NIN,IB)%bakMAT%UVAR(3),UVAR(3) / SOM
2076c print *, "merge uvar-2", BRICK_LIST(NIN,IB)%bakMAT%UVAR(2),UVAR(4) / SOM
2077c print *, "merge uvar-1", BRICK_LIST(NIN,IB)%bakMAT%UVAR(1),UVAR(5) / SOM
2078c print *, ""
2079
2080 mbuf%VAR((0*llt_ + idloc)) = uvar(1) / som !som=vol
2081 IF(som_(2)>em20)THEN
2082 mbuf%VAR((1*llt_ + idloc)) = uvar(2) / som_(2) !SOM2 = VOL_G =sum (vol_g , cell)
2083 mbuf%VAR((4*llt_ + idloc)) = uvar(5) / som
2084 ENDIF
2085 IF(som_(1)>em20)THEN
2086 mbuf%VAR((2*llt_ + idloc)) = uvar(3) / som_(1) !SOM3 = VOL_L =sum (vol_l , cell)
2087 mbuf%VAR((3*llt_ + idloc)) = uvar(4) / som
2088 ENDIF
2089 supercellvol_l(0:nthread-1,1:4,ib) = zero
2090
2091
2092 !APPEL SIGEPS37 POUR EQUILIBRE PRESSION
2093 brickid = brick_list(nin,ib)%ID
2094 idloc = idlocb(ib)
2095 ng = ngb(ib)
2096 vol = som
2097 CALL sigeps37_single_cell (
2098 1 elbuf_tab, ixs, bufmat, iparg, ipm,
2099 2 idloc , ng , brickid, vol
2100 . )
2101
2102
2103
2104 ELSEIF(mtn_==51)THEN
2105 llt_ = iparg(2,ngb(ib))
2106 uvar(1) = delta*brick_list(nin,ib)%bakMAT%UVAR(1) + sum(uvarl(0:nthread-1,ib,1))
2107 uvar(2) = delta*brick_list(nin,ib)%bakMAT%UVAR(2) + sum(uvarl(0:nthread-1,ib,2))
2108 uvar(3) = delta*brick_list(nin,ib)%bakMAT%UVAR(3) + sum(uvarl(0:nthread-1,ib,3))
2109 uvarl(0:nthread-1,ib,1:3) = zero
2110 mbuf%VAR((0*llt_ + idloc))= uvar(1) / som
2111 mbuf%VAR((1*llt_ + idloc))= uvar(2) / som
2112 mbuf%VAR((2*llt_ + idloc))= uvar(3) / som
2113 DO itrimat=1,trimat
2114 somi = sum(supercellvol_l(0:nthread-1,itrimat,ib))
2115 supercellvol_l(0:nthread-1,itrimat,ib) = zero
2116 ipos = 11
2117 k = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos)
2118 somi = somi + delta*brick_list(nin,ib)%bakMAT%UVAR(k)
2119 vol51(ib,itrimat) = brick_list(nin,ib)%bakMAT%UVAR(k)
2120 uvar(:) = zero
2121 IF(somi>zero)THEN
2122 DO ipos=1,m51_nvphas
2123 IF(ipos==11)THEN
2124 k1 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1) ! ici UVAR(1:M51_N0PHAS+TRIMAT*M51_NVPHAS)
2125 k = k1 * llt_
2126 mbuf%VAR(k+idlocb(ib)) = somi
2127 cycle
2128 ELSEIF(ipos==1)THEN
2129 k1 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1) ! ici UVAR(1:M51_N0PHAS+TRIMAT*M51_NVPHAS)
2130 k = k1 * llt_
2131 mbuf%VAR(k+idlocb(ib)) = somi / som
2132 cycle
2133 ENDIF
2134 k = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos)
2135 k1 = (k-1)*llt_
2136 pvar => mbuf%VAR(k1+idlocb(ib))
2137 uvar(k) = sum(uvarl(0:nthread-1,ib,k)) + pvar * vol51(ib,itrimat) * delta
2138 uvarl(0:nthread-1,ib,k) = zero
2139 uvar(k) = uvar(k) / somi
2140 IF(ipos/=1)pvar = uvar(k)
2141 enddo!next IPOS
2142 ELSE
2143 !vfrac_i=0 Vol_i=0
2144 ipos = 1
2145 k1 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1)
2146 k = k1 * llt_
2147 mbuf%VAR(k+idlocb(ib)) = zero
2148 ipos = 11
2149 k1 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1)
2150 k = k1 * llt_
2151 mbuf%VAR(k+idlocb(ib)) = zero
2152 endif!(SOMi>ZERO)
2153 enddo!next ITRIMAT
2154
2155 endif! (MTN_==51)
2156
2157 brick_list(nin,ib)%Vold_Scell= som
2158
2159
2160 enddo!next IB
2161
2162
2163
2164
2165 if(ibug22_merge/=0 )then
2166 call my_barrier
2167 if(itask==0)then
2168 if(dt1>zero)then
2169 do ib=1,nb
2170 if (tmp22array(1,ib)==zero )cycle
2171 write(*,fmt='(A,E30.16)') " TARGET=", nint(tmp22array(1,ib))
2172 write(*,fmt='(A,E30.16)') " DELTA=", tmp22array(2,ib)
2173 write(*,fmt='(A,E30.16)') " rho=", tmp22array(3,ib)
2174 write(*,fmt='(A,3E30.16)') " +rhoUx=", tmp22array(4,ib)
2175 write(*,fmt='(A,E30.16)') " +Vol=", tmp22array(5,ib)
2176 enddo
2177 endif!if(DT1>ZERO)
2178 !print *, "deallocate"
2179 deallocate(tmp22array)
2180 endif
2181 endif
2182
2183 !------------------------------------------------------------!
2184 ! @26. DEMERGING CRITERIA !
2185 ! (IF main HAS TO BE DEPORTED THEN !
2186 ! DEMERGING WILL OCCURS IN NEW main CELL) !
2187 ! STEP 1/2 - listing origin cell !
2188 !------------------------------------------------------------!
2189 ! if Ntar > 0 then material deporation has to be done from former main cell buffer (GBUF)
2190 ! to the N=Ntar target main cells.
2191 nbl1 = nbl
2192 IF(dt1==zero)nbl1 = 0
2193 origin_data(:,:,:) = 0 !tag main cells linked to these secnd ones.
2194 norigin(nbf:nbl) = 0
2195 DO ib=nbf,nbl1
2196 nbcut = brick_list(nin,ib)%NBCUT
2197 ncell = nbcut
2198 icell = 0
2199 mcell = brick_list(nin,ib)%mainID
2200 gbuf => elbuf_tab(ngb(ib))%GBUF
2201 pnodwasmain(1:8) => brick_list(nin,ib)%NODE(1:8)!%NodWasMain
2202 pwherewasmain(1:8) => brick_list(nin,ib)%NODE(1:8)!%WhereWasMain
2203 brickid = brick_list(nin,ib)%ID
2204 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
2205 ntag = 0
2206 mtn_ = iparg(1,ngb(ib))
2207 newmnod(1:8) = 0
2208 itag(:) = 0
2209 IF(mcell == 0)cycle
2210 mnod = brick_list(nin,ib)%POLY(mcell)%NumNOD
2211ccc print *, "demerging criteria for id=", ixs(11,brick_list(nin,ib)%id)
2212 !------CHECK NODES--------!
2213 DO k=1,mnod
2214 inod = brick_list(nin,ib)%POLY(mcell)%ListNodID(k)
2215 IF(pnodwasmain(inod)%NodWasMain==1)cycle !node not concerned
2216 !Here node was secnd and just became main
2217 newmnod(inod) = 1
2218 ENDDO
2219
2220 !Continuity Criteria To avoid unexpected demerging.
2221 !Consistent only with kinematic time step "several cycle to fully cross an element".
2222 vol_m = brick_list(nin,ib)%secndlist%VOL_unmerged !BRICK_LIST(NIN,IB)%POLY(MCELL)vnew
2223 uncutvol = brick_list(nin,ib)%UncutVol
2224 var = vol_m / uncutvol
2225ccc print *, " VAR > 1-CritDVol22 : CYCLE", VAR > 1-CritDVol22, VAR, VOL_M
2226ccc print *, " MNOD==0) : CYCLE", MNOD==0
2227!not correct with sharp edges
2228! IF (VAR > 1-CritDVol22)CYCLE
2229 IF(mnod==0) cycle
2230 CALL weighting_cell_nodes(nin,ib,mcell,inod_w,idemerge)
2231
2232 ! IF(pNodWasMain(BRICK_LIST(NIN,IB)%OldMainStrongNode) == 1) IDEMERGE = 0
2233
2234 inod_w_old = brick_list(nin,ib)%OldMainStrongNode
2235ccc print *, " INOD_W_old==0 : CYCLE", INOD_W_old==0
2236 IF(inod_w_old<=0)cycle
2237 IF(brick_list(nin,ib)%NODE(inod_w_old)%WhichCell == mcell) idemerge = 0
2238ccc print *, " IDEMERGE= : ", IDEMERGE
2239
2240 IF( idemerge==1 ) THEN !(negative in debug if IB new in the buffer)
2241! IF(pNodWasMain(INod_W)<=0 ) THEN !(negatif en debug si IB nouveau dans le buffer)
2242! IF(SUM(NewMnod)>=MNOD/2 .AND. MNOD/=0)THEN !DEMERGE criteria : all main nodes were secndo ones
2243 !get related main cell location (through its face)
2244 DO k=1,mnod
2245 inod = brick_list(nin,ib)%POLY(mcell)%ListNodID(k)
2246 j = pwherewasmain(inod)%WhereWasMain
2247 IF(j==0)cycle !means node was already main
2248 IF(itag(j)==1)cycle !list without doublon
2249 itag(j) = 1
2250 !Check if main cell was moved (merged)
2251 !Get Main Cut Cell Address
2252 IF(j>=10)THEN
2253 j1 = j/10
2254 j2 = mod(j,10)
2255 ibv_i = brick_list(nin,ib )%Adjacent_Brick(j1,4)
2256 ibv = brick_list(nin,ibv_i)%Adjacent_Brick(j2,4)
2257 ELSE
2258 ibv = brick_list(nin,ib )%Adjacent_Brick(j,4)
2259 ENDIF
2260 ntar = brick_list(nin,ibv)%Ntarget
2261 IF(ntar == 0)THEN
2262 ntag = ntag + 1
2263 origin_data(ib,ntag,1) = ibv
2264 origin_data(ib,ntag,2) = j
2265 origin_data(ib,ntag,3) = ibv
2266 ELSE
2267 DO itar=1,ntar
2268 ntag = ntag + 1
2269 origin_data(ib,ntag,1) = brick_list(nin,ibv)%MergeTarget(3,itar)
2270 origin_data(ib,ntag,2) = j
2271 origin_data(ib,ntag,3) = ibv
2272 ENDDO
2273 ENDIF
2274 ENDDO
2275 norigin(ib) = ntag
2276 IF(ntag==0)THEN
2277 print *, "**error : inter22, topology issue."
2278 ENDIF
2279 ENDIF
2280 enddo! next IB
2281
2282 !--------------!
2283 CALL my_barrier
2284 !--------------!
2285
2286 !------------------------------------------------------------!
2287 ! @27. MATERIAL DEMERGING !
2288 ! STEP 2 : filling main cell with data from origin cells !
2289 !------------------------------------------------------------!
2290 ! if Ntar > 0 then material deporation has to be done from former main cell buffer (GBUF)
2291 ! to the N=Ntar target main cells.
2292 nbl1 = nbl
2293 IF(dt1==zero)nbl1 = 0
2294
2295 if(ibug22_merge/=0)then
2296 if(itask==0)ALLOCATE (tmp22array(8,nb))
2297 in22 = 0
2298 call my_barrier
2299 tmp22array(1:8,nbf:nbl1)=zero
2300 endif
2301
2302
2303 DO ib=nbf,nbl1
2304 IF(norigin(ib)==0)cycle
2305 nbcut = brick_list(nin,ib)%NBCUT
2306 IF(nbcut==0)cycle
2307 ncell = nbcut
2308 icell = 0
2309 mcell = brick_list(nin,ib)%mainID
2310 gbuf => elbuf_tab(ngb(ib))%GBUF
2311 lbuf => elbuf_tab(ngb(ib))%BUFLY(1)%LBUF(1,1,1)
2312 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
2313 padjbrick => brick_list(nin,ib)%Adjacent_Brick(1:6,1:3)
2314 pnodwasmain(1:8) => brick_list(nin,ib)%NODE(1:8)!%NodWasMain
2315 pwherewasmain(1:8) => brick_list(nin,ib)%NODE(1:8)!%WhereWasMain
2316 psubvol(1:9) => brick_list(nin,ib)%POLY(1:9)
2317 brickid = brick_list(nin,ib)%ID
2318 idloc = brick_list(nin,ib)%IDLOC
2319 mtn_ = iparg(1,ngb(ib))
2320 llt_ = iparg(2,ngb(ib))
2321 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
2322 ntag = 0
2323 volorig(1:24) = zero
2324 volcell = zero
2325 volcell51(:) = zero
2326 volsecnd51(:) = zero
2327 vol = zero
2328 rho = zero
2329 eint = zero
2330 sig(1:6) = zero
2331 mom(1:3) = zero
2332 ssp = zero
2333 uvar(:) = zero
2334 rho_adv = zero
2335 eint_adv = zero
2336 sig_adv(1:6) = zero
2337 mom_adv(1:3) = zero
2338 uvar_adv(:) = zero
2339
2340 !-----------------------------------------------!
2341 ! INIT CELL CONTENT IF IT CAME FROM A MERGE !
2342 !-----------------------------------------------!
2343 IF(brick_list(nin,ib)%IsMergeTarget==1)THEN
2344 ! pause
2345 !!!print *, "ismergeTarget , Norigin", ixs(11,brickID), Norigin(IB), BRICK_LIST(NIN,IB)%Ntarget
2346 ! This cell received material from a merge.it also has to be taken into account
2347 volcell = brick_list(nin,ib)%Vold_Scell
2348 rho = volcell* gbuf%RHO(idloc)
2349 eint = volcell* gbuf%EINT(idloc)
2350 mom(1) = volcell* gbuf%MOM(llt_*(1-1) + idloc)
2351 mom(2) = volcell* gbuf%MOM(llt_*(2-1) + idloc)
2352 mom(3) = volcell* gbuf%MOM(llt_*(3-1) + idloc)
2353 sig(1) = volcell* gbuf%SIG(llt_*(1-1)+idloc)
2354 sig(2) = volcell* gbuf%SIG(llt_*(2-1)+idloc)
2355 sig(3) = volcell* gbuf%SIG(llt_*(3-1)+idloc)
2356 sig(4) = volcell* gbuf%SIG(llt_*(4-1)+idloc)
2357 sig(5) = volcell* gbuf%SIG(llt_*(5-1)+idloc)
2358 sig(6) = volcell* gbuf%SIG(llt_*(6-1)+idloc)
2359 ssp = volcell* lbuf%SSP(idloc)
2360 volcell = zero
2361 IF(mtn_==37)THEN
2362 uvar(5) = volcell*mbuf%VAR(((5-1)*llt_ + idloc))
2363 uvar(4) = volcell*mbuf%VAR(((4-1)*llt_ + idloc))
2364 uvar(3) = volcell*mbuf%VAR(((3-1)*llt_ + idloc))*mbuf%VAR(((4-1)*llt_ + idloc))
2365 uvar(2) = volcell*mbuf%VAR(((2-1)*llt_ + idloc))*mbuf%VAR(((5-1)*llt_ + idloc))
2366 uvar(1) = volcell*mbuf%VAR(((1-1)*llt_ + idloc))
2367 ELSEIF(mtn_==51)THEN
2368 ipos = 11
2369 volsecnd51(1) = mbuf%VAR(idloc + ((m51_n0phas + (1-1)*m51_nvphas )+ipos-1)*llt_)
2370 volsecnd51(2) = mbuf%VAR(idloc + ((m51_n0phas + (2-1)*m51_nvphas )+ipos-1)*llt_)
2371 volsecnd51(3) = mbuf%VAR(idloc + ((m51_n0phas + (3-1)*m51_nvphas )+ipos-1)*llt_)
2372 volsecnd51(4) = mbuf%VAR(idloc + ((m51_n0phas + (4-1)*m51_nvphas )+ipos-1)*llt_)
2373 DO itrimat=1,trimat
2374 DO ipos = 1 , m51_nvphas
2375 IF(ipos==11) cycle
2376 k = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2377 k1 = k * llt_
2378 uvar(k+1) = volsecnd51(itrimat) * mbuf%VAR(k1+idloc)
2379 enddo!next IPOS
2380 enddo!next ITRIMAT
2381 volsecnd51(1:4) = zero
2382 ENDIF
2383 ENDIF
2384
2385 ! Velocity of main is velocity wherever inside the supercell. It will be used to make an advection step get pressure gradient which did not occur on previous step due to homogeneity.
2386 vel_m(1) = gbuf%MOM(3*(idloc-1)+1) / gbuf%RHO(idloc)
2387 vel_m(2) = gbuf%MOM(3*(idloc-1)+2) / gbuf%RHO(idloc)
2388 vel_m(3) = gbuf%MOM(3*(idloc-1)+3) / gbuf%RHO(idloc)
2389
2390 ntar = brick_list(nin,ib)%Ntarget
2391 !-----------------------------------------------!
2392 ! LOOP ON ORIGIN CELLS !
2393 !-----------------------------------------------!
2394 DO iorig = 1, norigin(ib)
2395 !IORIG data regarding IB
2396 ibm = origin_data(ib,iorig,1) !real target
2397 ibo = origin_data(ib,iorig,3) !original brick before knowing it was burged below
2398 iv = brick_list(nin,ibm)%ID
2399 ngv = brick_list(nin,ibm)%NG
2400 idlocv = brick_list(nin,ibm)%IDLOC
2401 j = origin_data(ib,iorig,2)
2402 IF(ibm==ib)cycle
2403 IF(j>=10)THEN
2404 j1 = j/10
2405 j2 = mod(j,10)
2406 ibv = brick_list(nin,ib )%Adjacent_Brick(j1,4)
2407 jv = brick_list(nin,ibv)%Adjacent_Brick(j2,5)
2408 ELSE
2409 jv = brick_list(nin,ib)%Adjacent_Brick(j,5)
2410 ENDIF
2411
2412 numsecnd = old_secndlist(nin,ibo)%Num !this is list from previous cycle
2413 !-----------------------------------------------!
2414 ! GET SECND CONNECTED TO IT !
2415 !-----------------------------------------------!
2416 DO k=1,numsecnd
2417 IF (old_secndlist(nin,ibo)%IBv(k) /= ib)cycle
2418 IF (old_secndlist(nin,ibo)%FM(k) /= jv)cycle !next secnd : this one is not connecter through main face id %FM
2419
2420 !---debug---!
2421 cod1 = ixs(11,brick_list(nin,ib)%id)
2422 cod2 = ixs(11,brick_list(nin,ibo)%id)
2423 if (ibug22_merge==-1)then
2424#include "lockon.inc"
2425 if(ibm==ibo)then
2426 ! write(*,FMT='(A,I8,A,I8)') "DEMERGING : id=", Cod1,
2427 ! . " from idv:", ixs(11,brick_list(nin,ibo)%id)
2428
2429 in22 = in22 + 1
2430 tmp22array(1,in22)=cod1
2431 tmp22array(2,in22)=ixs(11,brick_list(nin,ibo)%id)
2432 tmp22array(3,in22)=0
2433
2434 else
2435 ! write(*,FMT='(A,I8,A,I8,A,I8)') "DEMERGING : id=", Cod1,
2436 ! . " from idv:", ixs(11,brick_list(nin,ibo)%id)," moved to", ixs(11,brick_list(nin,ibm)%id)
2437
2438 in22 = in22 + 1
2439 tmp22array(1,in22)=cod1
2440 tmp22array(2,in22)=ixs(11,brick_list(nin,ibo)%id)
2441 tmp22array(3,in22)=ixs(11,brick_list(nin,ibm)%id)
2442
2443 endif
2444 ! write (*,FMT='(I10,A,F30.16,A,F30.16)')
2445 ! . Cod1," Vold=",brick_list(nin,ib)%Vold_scell, " Vnew=",brick_list(nin,ib)%Vnew_scell
2446 ! write (*,FMT='(I10,A,F30.16,A,F30.16)')
2447 ! . cod2," Vold=",brick_list(nin,ibo)%Vold_scell," Vnew=",brick_list(nin,ibo)%Vnew_scell
2448
2449 tmp22array(4,in22)= zero !brick_list(nin,ib)%Vold_scell !value not consistent, old main was merge elsewhere so scell has change of side since a demerge process is emerging a new scell
2450 tmp22array(5,in22)= brick_list(nin,ib)%Vnew_scell
2451 tmp22array(6,in22)= zero !brick_list(nin,ibo)%Vold_scell
2452 tmp22array(7,in22)= brick_list(nin,ibo)%Vnew_scell
2453 tmp22array(8,in22)=tmp22array(8,in22)+1
2454#include "lockoff.inc"
2455 endif
2456 !-----------!
2457
2458 !here secnd is merged on the correct face, and also to the correct target main (IB)
2459 !!!!!!!!!!!!! VOLsecnd = OLD_SecndList(NIN,IBo)%VOL(K)
2460 !-----------------------------------------------!
2461 ! START TO CALCULATE SECND NEW QUANTITIES !
2462 !-----------------------------------------------!
2463 volsecnd = old_secndlist(nin,ibo)%VOL(k)
2464 !!! VolsecndOLD = OLD_SecndList(NIN,IBo)%VOL(K)
2465 volcell = volcell + volsecnd !cumulative salve volume from previous cycle = old volume for new main cell
2466 !!! VOLcellOLD = VOLcell + VOLsecndOLD
2467 gbuf => elbuf_tab(ngv)%GBUF
2468 lbuf => elbuf_tab(ngv)%BUFLY(1)%LBUF(1,1,1)
2469 llt_v = iparg(2,ngv)
2470 !secnd cell surf
2471 surf_s = old_secndlist(nin,ibo)%SURF_V(k)
2472 fv = old_secndlist(nin,ibo)%FV(k)
2473 IF(fv<10)THEN
2474 norm_s(1:3) = brick_list(nin,ib)%N(fv,1:3)
2475 ELSE
2476 !FV = FV/10
2477 !NORM_S(1:3) = BRICK_LIST(NIN,IB)%N(FV,1:3)
2478 !There are several possible facets for indirect neighbors.It would be necessary to SAVAGARGER J> 10 for each path, then advection
2479 !or change the normal by taking the vector centroid_M-centroid_S
2480 zm(1:3) = brick_list(nin,ibo)%SCellCenter(1:3)
2481 zs(1:3) = brick_list(nin,ib)%POLY(mcell)%CellCenter(1:3)
2482 norm_s(1) = zm(1)-zs(1)
2483 norm_s(2) = zm(2)-zs(2)
2484 norm_s(3) = zm(3)-zs(3)
2485 norm = sqrt(norm_s(1)*norm_s(1)+norm_s(2)*norm_s(2)+norm_s(3)*norm_s(3))
2486 norm_s(1) = norm_s(1) / norm
2487 norm_s(2) = norm_s(2) / norm
2488 norm_s(3) = norm_s(3) / norm
2489 ENDIF
2490
2491
2492 !----------------------------------!
2493 ! ADVECTION STEP IN DEMERGED CELL !
2494 !----------------------------------!
2495 adv = zero
2496 IF(iadv==1)THEN
2497 !MCELL = BRICK_LIST(NIN,IB)%MainID
2498 vm = one ! BRICK_LIST(NIN,IBm)%POLY(MCELL)%Vnew
2499 adv = -half* half* dt1*vm*surf_s* (vel_m(1)*norm_s(1)+vel_m(2)*norm_s(2)+vel_m(3)*norm_s(3))
2500 rho_adv = rho_adv + adv * gbuf%RHO(idlocv)
2501 eint_adv = eint + adv * gbuf%EINT(idlocv)
2502 mom_adv(1) = mom_adv(1) + adv * gbuf%MOM(3*(idlocv-1)+1)
2503 mom_adv(2) = mom_adv(2) + adv * gbuf%MOM(3*(idlocv-1)+2)
2504 mom_adv(3) = mom_adv(3) + adv * gbuf%MOM(3*(idlocv-1)+3)
2505 sig_adv(1) = sig_adv(1) + adv * gbuf%SIG(llt_v*(1-1)+idlocv)
2506 sig_adv(2) = sig_adv(2) + adv * gbuf%SIG(llt_v*(2-1)+idlocv)
2507 sig_adv(3) = sig_adv(3) + adv * gbuf%SIG(llt_v*(3-1)+idlocv)
2508 sig_adv(4) = sig_adv(4) + adv * gbuf%SIG(llt_v*(4-1)+idlocv)
2509 sig_adv(5) = sig_adv(5) + adv * gbuf%SIG(llt_v*(5-1)+idlocv)
2510 sig_adv(6) = sig_adv(6) + adv * gbuf%SIG(llt_v*(6-1)+idlocv)
2511 ENDIF
2512 !-------------------------------!
2513 ! MONO-MATERIAL CUMULATION !
2514 !-------------------------------!
2515 rho = rho + volsecnd * gbuf%RHO(idlocv)
2516 eint = eint + volsecnd * gbuf%EINT(idlocv)
2517 mom(1) = mom(1) + volsecnd * gbuf%MOM(llt_v*(1-1)+idlocv)
2518 mom(2) = mom(2) + volsecnd * gbuf%MOM(llt_v*(2-1)+idlocv)
2519 mom(3) = mom(3) + volsecnd * gbuf%MOM(llt_v*(3-1)+idlocv)
2520 sig(1) = sig(1) + volsecnd * gbuf%SIG(llt_v*(1-1)+idlocv)
2521 sig(2) = sig(2) + volsecnd * gbuf%SIG(llt_v*(2-1)+idlocv)
2522 sig(3) = sig(3) + volsecnd * gbuf%SIG(llt_v*(3-1)+idlocv)
2523 sig(4) = sig(4) + volsecnd * gbuf%SIG(llt_v*(4-1)+idlocv)
2524 sig(5) = sig(5) + volsecnd * gbuf%SIG(llt_v*(5-1)+idlocv)
2525 sig(6) = sig(6) + volsecnd * gbuf%SIG(llt_v*(6-1)+idlocv)
2526 ssp = ssp + volsecnd * lbuf%SSP(idlocv)
2527 !-------------------------------!
2528 ! MULTI-MATERIAL CUMULULATION !
2529 !-------------------------------!
2530 IF(mtn_==37)THEN
2531 !UVAR(I,1) : massic percentage of liquid * global density (rho1*V1/V : it needs to give liquid mass multiplying by element volume in aleconve.F)
2532 !UVAR(I,2) : density of gas
2533 !UVAR(I,3) : density of liquid
2534 !UVAR(I,4) : volumetric fraction of liquid
2535 !UVAR(I,5) : volumetric fraction of gas
2536 mbufv => elbuf_tab(ngv)%BUFLY(1)%MAT(1,1,1)
2537 llt_v = iparg(2,ngv)
2538 uvar(5) = uvar(5) +volsecnd*mbufv%VAR(((5-1)*llt_v + idlocv))
2539 uvar(4) = uvar(4) +volsecnd*mbufv%VAR(((4-1)*llt_v + idlocv))
2540 uvar(3) = uvar(3) +volsecnd*mbufv%VAR(((3-1)*llt_v + idlocv))*mbufv%VAR(((4-1)*llt_v + idlocv))
2541 uvar(2) = uvar(2) +volsecnd*mbufv%VAR(((2-1)*llt_v + idlocv))*mbufv%VAR(((5-1)*llt_v + idlocv))
2542 uvar(1) = uvar(1) +volsecnd*mbufv%VAR(((1-1)*llt_v + idlocv))
2543 uvar_adv(1) = uvar_adv(1) +adv * mbufv%VAR(((1-1)*llt_v + idlocv))
2544
2545 ELSEIF(mtn_==51)THEN
2546 mbufv => elbuf_tab(ngv)%BUFLY(1)%MAT(1,1,1)
2547 llt_ = iparg(2,ngb(ib))
2548 llt_v = iparg(2,ngv)
2549 !in sigesp51.F:
2550 !V1OLD = V1OLD - TIMESTEP*DDVOL1
2551 !Volume change correspond to DDVOL=DV/DT = f(fluxes). Gird deformation influences relative material velocity ans also fluxes
2552 DO itrimat=1,trimat
2553 !get,store Volumetric Fraction
2554 ipos = 1
2555 k1 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2556 kv = k1 * llt_v
2557 vfrac(itrimat) = mbufv%VAR(kv+idlocv)
2558 volsecnd51(itrimat) = vfrac(itrimat) * volsecnd !submaterial volume in secnd cell
2559 volcell51(itrimat) = volcell51(itrimat) + volsecnd51(itrimat) !submaterial volume in ALL secnd cell
2560 enddo!next ITRIMAT
2561 uvar(1) = uvar(1) + mbufv%VAR((0*llt_v + idlocv)) * volsecnd
2562 uvar(2) = uvar(2) + mbufv%VAR((1*llt_v + idlocv)) * volsecnd
2563 uvar(3) = uvar(3) + mbufv%VAR((2*llt_v + idlocv)) * volsecnd
2564 uvar_adv(1) = uvar_adv(1) + mbufv%VAR((0*llt_v + idlocv)) * adv
2565 uvar_adv(2) = uvar_adv(2) + mbufv%VAR((1*llt_v + idlocv)) * adv
2566 uvar_adv(3) = uvar_adv(3) + mbufv%VAR((2*llt_v + idlocv)) * adv
2567 IF(iadv==1)print *, "to do compute/get vfrac on face" !; pause
2568 DO itrimat=1,trimat
2569 DO ipos = 1 , m51_nvphas
2570 ko = ((m51_n0phas + (itrimat-1)*m51_nvphas ))
2571 k1 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2572 kv = k1 * llt_v
2573 IF(ipos==11)THEN
2574 uvar(k1+1) = uvar(k1+1) + volsecnd51(itrimat)
2575 ELSE
2576 uvar(k1+1) = uvar(k1+1) + volsecnd51(itrimat) * mbufv%VAR(kv+idlocv)
2577 uvar_adv(k1+1) = uvar_adv(k1+1) + adv * mbufv%VAR(kv+idlocv)
2578 IF(iadv==1)print *, "to do compute/get vfrac on face" !; pause
2579 ENDIF
2580 enddo!next IPOS
2581 enddo!next ITRIMAT
2582 endif!(MTN_==51)
2583 !---------------------------
2584#include "lockon.inc"
2585 !lock method : to be update later can be removed
2586 !SUBSTRACTING SECND VOLUME FROM SUPERCELL TO DEMERGE
2587 !--MONOMAT
2588 brick_list(nin,ibm)%Vold_Scell = brick_list(nin,ibm)%Vold_Scell - volsecnd
2589 !--MULTIMAT
2590 IF(mtn_==51)THEN
2591 DO itrimat=1,trimat
2592 !remove submaterial volume in detached secnd cell
2593 ipos = 11
2594 k1 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2595 kv = k1 * llt_v
2596 mbufv%VAR(kv+idlocv) = mbufv%VAR(kv+idlocv) - volsecnd51(itrimat)
2597 enddo!next ITRIMAT
2598 ENDIF
2599#include "lockoff.inc"
2600 enddo!next K
2601 enddo!next IORIG
2602
2603 IF(volcell==zero)cycle
2604
2605 !-----------------------------------------------!
2606 ! GET VALUE FROM CUMULATIVE ONE !
2607 !---MONOMAT-------------------------------------!
2608 rho = (rho +iadv*rho_adv )/ volcell
2609 eint = (eint +iadv*eint_adv )/ volcell
2610 mom(1) = (mom(1) +iadv*mom_adv(1))/ volcell
2611 mom(2) = (mom(2) +iadv*mom_adv(2))/ volcell
2612 mom(3) = (mom(3) +iadv*mom_adv(3))/ volcell
2613 sig(:) = (sig(:) +iadv*sig_adv(:))/ volcell
2614 ssp = ssp / volcell ! -> run EOS
2615 !---MULTIMAT -----------------------------------!
2616 IF(mtn_==37)THEN
2617 uvar(1) = (uvar(1)+ iadv*uvar_adv(1)) / volcell
2618 uvar(2) = uvar(2) / volcell
2619 uvar(3) = uvar(3) / volcell
2620 uvar(4) = uvar(4) / volcell
2621 uvar(5) = uvar(5) / volcell
2622 !call sigeps37_single_cell before storing in MBUF%VAR
2623 ELSEIF(mtn_==51)THEN
2624 uvar(1) = (uvar(1) + iadv*uvar_adv(1)) / volcell !sum of full volume because UVAR(1:3) are homogenized data, not submaterial one.
2625 uvar(2) = (uvar(2) + iadv*uvar_adv(2)) / volcell
2626 uvar(3) = (uvar(3) + iadv*uvar_adv(3)) / volcell
2627 DO itrimat=1,trimat
2628 ko = m51_n0phas + (itrimat-1)*m51_nvphas
2629 IF(volcell51(itrimat)/=zero)THEN
2630 uvar(ko+01:ko+10) = (uvar(ko+01:ko+10) + iadv*uvar_adv(ko+01:ko+10)) / volcell51(itrimat) !submaterial data are averaged with sum of submaterial secnd volume
2631 uvar(ko+12:ko+m51_nvphas) = (uvar(ko+12:ko+m51_nvphas) + iadv*uvar_adv(ko+12:ko+m51_nvphas)) / volcell51(itrimat) !except ipos=11 (volume)
2632 ELSE
2633 uvar(ko+01) = zero
2634 uvar(ko+11) = zero
2635 ENDIF
2636 ENDDO
2637 endif!(MTN_==51)
2638 !-----------------------------------------------!
2639
2640 !-----------------------------------------------!
2641 ! STORE VALUES IN DEMERGED CELLS !
2642 !---MONOMAT-------------------------------------!
2643 ng = brick_list(nin,ib)%NG
2644 idloc = brick_list(nin,ib)%IDLOC
2645 llt_ = iparg(2,ng)
2646 gbuf => elbuf_tab(ng)%GBUF
2647 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
2648 gbuf%RHO(idloc) = rho
2649 gbuf%EINT(idloc) = eint
2650 gbuf%MOM(llt_*(1-1)+idloc) = mom(1)
2651 gbuf%MOM(llt_*(2-1)+idloc) = mom(2)
2652 gbuf%MOM(llt_*(3-1)+idloc) = mom(3)
2653 gbuf%SIG(llt_*(1-1)+idloc) = sig(1)
2654 gbuf%SIG(llt_*(2-1)+idloc) = sig(2)
2655 gbuf%SIG(llt_*(3-1)+idloc) = sig(3)
2656 gbuf%SIG(llt_*(4-1)+idloc) = sig(4)
2657 gbuf%SIG(llt_*(5-1)+idloc) = sig(5)
2658 gbuf%SIG(llt_*(6-1)+idloc) = sig(6)
2659 lbuf%SSP(idloc) = ssp
2660 alefvm_buffer%FCELL(4, brick_list(nin,ib)%ID ) = rho
2661 alefvm_buffer%FCELL(1, brick_list(nin,ib)%ID ) = mom(1)
2662 alefvm_buffer%FCELL(2, brick_list(nin,ib)%ID ) = mom(2)
2663 alefvm_buffer%FCELL(3, brick_list(nin,ib)%ID ) = mom(3)
2664 alefvm_buffer%FCELL(6, brick_list(nin,ib)%ID ) = -third*(sig(1)+sig(2)+sig(3))
2665 alefvm_buffer%FCELL(5, brick_list(nin,ib)%ID ) = ssp
2666
2667 !---MUKLTIMAT-------------------------------------!
2668 IF(mtn_==37)THEN
2669 llt_ = iparg(2,ngb(ib))
2670 idloc = idlocb(ib)
2671 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
2672 mt = ixs(1,brickid)
2673 iadbuf = ipm(7,mt)
2674 rho10 = bufmat(iadbuf-1+11)
2675 rho20 = bufmat(iadbuf-1+12)
2676c print *, ""
2677c print *, "demerge-5", MBUF%VAR((5-1)*LLT_+IDLOC) , UVAR(5)
2678c print *, "demerge-4", MBUF%VAR((4-1)*LLT_+IDLOC) , UVAR(4)
2679c print *, "demerge-3", MBUF%VAR((3-1)*LLT_+IDLOC) , UVAR(3)
2680c print *, "demerge-2", MBUF%VAR((2-1)*LLT_+IDLOC) , UVAR(2)
2681c print *, "demerge-1", MBUF%VAR((1-1)*LLT_+IDLOC) , UVAR(1)
2682 IF(uvar(3) == zero) uvar(3) = rho10
2683 IF(uvar(2) == zero) uvar(2) = rho20
2684 mbuf%VAR((5-1)*llt_+idloc) = uvar(5)
2685 mbuf%VAR((4-1)*llt_+idloc) = uvar(4)
2686 mbuf%VAR((3-1)*llt_+idloc) = uvar(3)
2687 mbuf%VAR((2-1)*llt_+idloc) = uvar(2)
2688 mbuf%VAR((1-1)*llt_+idloc) = uvar(1)
2689 ELSEIF(mtn_==51)THEN
2690 llt_ = iparg(2,ngb(ib))
2691 idloc = idlocb(ib)
2692 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
2693 DO itrimat=1,trimat
2694 IF(volcell51(itrimat)/=zero)THEN
2695 !get,store Volumetric Fraction
2696 DO ipos=1,m51_nvphas
2697 k0 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2698 k1 = k0 * llt_
2699 mbuf%VAR(k1+idloc) = uvar(m51_n0phas+(itrimat-1)*m51_nvphas+ipos)
2700 ENDDO
2701
2702 ipos=1
2703 k0 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2704 k1 = k0 * llt_
2705 mbuf%VAR(k1+idloc) = volcell51(itrimat)/volcell
2706 ipos=11
2707 k0 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2708 k1 = k0 * llt_
2709 mbuf%VAR(k1+idloc) = volcell51(itrimat)
2710 ELSE
2711 ipos=1
2712 k0 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2713 k1 = k0 * llt_
2714 mbuf%VAR(k1+idloc) = zero
2715 ipos=11
2716 k0 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2717 k1 = k0 * llt_
2718 mbuf%VAR(k1+idloc) = zero
2719 ENDIF
2720 ENDDO
2721 endif!(MTN_==51)
2722 !-----------------------------------------------!
2723
2724 !-----------------------------------------------!
2725 ! NEW main CELL UPDATE !
2726 !-----------------------------------------------!
2727 brick_list(nin,ib)%Vold_SCell = volcell
2728
2729 !APPEL SIGEPS37 POUR EQUILIBRE PRESSION
2730 brickid = brick_list(nin,ib)%ID
2731 idloc = idlocb(ib)
2732 ng = ngb(ib)
2733 vol = volcell
2734 mtn_ = iparg(1,ng)
2735 IF(mtn_==37)THEN
2736 CALL sigeps37_single_cell (
2737 1 elbuf_tab, ixs, bufmat, iparg, ipm,
2738 2 idloc , ng , brickid, vol
2739 . )
2740 ENDIF
2741
2742 !
2743 !VOLcell = sum (VOLslacve) : each volsecnd is removed from
2744 !its former main cell and cumulative value is set in new
2745 !main cell to demerge.
2746 !-----------------------------------------------------!
2747 ! ADVECTION STEP UPDATE ORIGIN CELL (PREVIOUS MAIN) !
2748 !-----------------------------------------------------!
2749 IF(iadv==1)THEN
2750 ng = brick_list(nin,ibm)%NG
2751 idloc = brick_list(nin,ibm)%IDLOC
2752 gbuf => elbuf_tab(ng)%GBUF
2753 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
2754 gbuf%RHO(idloc) = gbuf%RHO(idloc) - rho_adv/volcell
2755 gbuf%EINT(idloc) = gbuf%EINT(idloc) - eint_adv/volcell
2756 gbuf%MOM(3*(idloc-1) +1) = gbuf%MOM(3*(idloc-1) +1) - mom_adv(1)/volcell
2757 gbuf%MOM(3*(idloc-1) +2) = gbuf%MOM(3*(idloc-1) +2) - mom_adv(2)/volcell
2758 gbuf%MOM(3*(idloc-1) +3) = gbuf%MOM(3*(idloc-1) +3) - mom_adv(3)/volcell
2759 gbuf%SIG(llt_*(1-1) +idloc) = gbuf%SIG(llt_*(1-1) +idloc) - sig_adv(1)/volcell
2760 gbuf%SIG(llt_*(2-1) +idloc) = gbuf%SIG(llt_*(2-1) +idloc) - sig_adv(2)/volcell
2761 gbuf%SIG(llt_*(3-1) +idloc) = gbuf%SIG(llt_*(3-1) +idloc) - sig_adv(3)/volcell
2762 gbuf%SIG(llt_*(4-1) +idloc) = gbuf%SIG(llt_*(4-1) +idloc) - sig_adv(4)/volcell
2763 gbuf%SIG(llt_*(5-1) +idloc) = gbuf%SIG(llt_*(5-1) +idloc) - sig_adv(5)/volcell
2764 gbuf%SIG(llt_*(6-1) +idloc) = gbuf%SIG(llt_*(6-1) +idloc) - sig_adv(6)/volcell
2765 mtn_ = iparg(1,ng)
2766 IF(mtn_==37)THEN
2767 llt_ = iparg(2,ng)
2768 idloc = idlocb(ibm)
2769 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
2770 brickid = brick_list(nin,ibm)%ID
2771 mt = ixs(1,brickid)
2772 iadbuf = ipm(7,mt)
2773 rho10 = bufmat(iadbuf-1+11)
2774 rho20 = bufmat(iadbuf-1+12)
2775 IF(uvar(3) == zero) uvar(3) = rho10
2776 IF(uvar(2) == zero) uvar(2) = rho20
2777 !UVAR(I,1) : massic percentage of liquid * global density (rho1*V1/V : it needs to give liquid mass multiplying by element volume in aleconve.F)
2778 !UVAR(I,2) : density of gas
2779 !UVAR(I,3) : density of liquid
2780 !UVAR(I,4) : volumetric fraction of liquid
2781 !UVAR(I,5) : volumetric fraction of gas
2782 mbuf%VAR((1-1)*llt_+idloc) = uvar(1) - uvar_adv(1)/volcell
2783
2784 !APPEL SIGEPS37 POUR EQUILIBRE PRESSION
2785 brickid = brick_list(nin,ibm)%ID
2786 idloc = idlocb(ibm)
2787 ng = ngb(ibm)
2788 vol = brick_list(nin,ibm)%Vold_Scell
2789 mtn_ = iparg(1,ng)
2790 IF(mtn_==37)THEN
2791 CALL sigeps37_single_cell (
2792 1 elbuf_tab, ixs, bufmat, iparg, ipm,
2793 2 idloc , ng , brickid, vol
2794 . )
2795 ENDIF
2796
2797 ELSEIF(mtn_==51)THEN
2798 ! print *, "todo" ; pause
2799 ENDIF
2800 ENDIF
2801
2802
2803 !-----------------------------------------------!
2804 enddo!next IB
2805
2806 if (ibug22_merge==-1)then
2807 call my_barrier
2808 if(itask==0 .and. dt1>zero)then
2809 allocate(order(in22),value(in22))
2810 DO ib=1,in22
2811 value(ib)=tmp22array(1,ib)
2812 ENDDO
2813 order=0
2814 !CALL QUICKSORT_I2 !(ORDER,IN22,value)
2815 do ib2=1,in22
2816 ib=order(ib2)
2817 cod1=tmp22array(1,ib)
2818 cod2=tmp22array(2,ib)
2819 cod3=tmp22array(3,ib)
2820 !if(cod1==0)cycle
2821 if(cod3==0)then
2822 write(*,fmt='(A,I8,A,I8)') "DEMERGING : id=", cod1, " from idv:", cod2
2823 else
2824 write(*,fmt='(A,I8,A,I8,A,I8)') "DEMERGING : id=", cod1, " from idv:", cod2," moved to", cod3
2825 endif
2826 write (*,fmt='(I10,A,F30.16,A,F30.16)') cod1," Vold=", tmp22array(4,ib), " Vnew=", tmp22array(5,ib)
2827 write (*,fmt='(I10,A,F30.16,A,F30.16)') cod2," Vold=", tmp22array(6,ib), " Vnew=", tmp22array(7,ib)
2828 write (*,fmt='(A,F30.16)') " Number of access =", tmp22array(8,ib)
2829 enddo
2830 endif
2831 call my_barrier
2832 if(itask==0)then
2833 if(allocated(tmp22array))deallocate(tmp22array)
2834 if(dt1>zero .and. allocated(order) )deallocate(order)
2835 endif
2836 endif
2837
2838
2839 !------------------------------------------------------------!
2840 ! @28. COMPUTE VOLD_Vnew_Cell for newly cell in buffer !
2841 !------------------------------------------------------------!
2842 !OLD_Vnew_Cell(1:9) are old cell volumes for a cut cell.
2843 !It is taken from previous cycle.
2844 !If brick was not cut on previous cycle then it must be initialized here
2845 IF(dt1==zero)nbl1 = 0
2846 DO ib=nbf,nbl1
2847 nbcut = brick_list(nin,ib)%NBCUT
2848 wascut = brick_list(nin,ib)%WasCut
2849 newinbuffer = brick_list(nin,ib)%NewInBuffer
2850CCC IF(WASCUT==0.and. NBCUT==0)THEN
2851 IF( newinbuffer == 1 )THEN
2852
2853 mcell = 1
2854 brick_list(nin,ib)%POLY(mcell)%OLD_Vnew = brick_list(nin,ib)%POLY(mcell)%Vnew
2855 brick_list(nin,ib)%POLY(2:9)%OLD_Vnew = zero
2856! print *, "new brick in cut cell buffer:", IXS(11,brick_list(nin,ib)%id)
2857! print *, "new brick in cut cell buffer: BRICK_LIST(NIN,IB)%POLY(MCELL)%OLD_Vnew=",BRICK_LIST(NIN,IB)%POLY(MCELL)%Vnew
2858
2859 ng = brick_list(nin,ib)%NG
2860 idloc = brick_list(nin,ib)%IDLOC
2861 id = brick_list(nin,ib)%ID
2862 gbuf => elbuf_tab(ng)%GBUF
2863 brick_list(nin,ib)%Vold_Scell = gbuf%VOL(idloc)
2864
2865c print *, "NEW IN BUFFER : ", ixs(11,id)
2866c print *, " vold_scell : ", brick_list(nin,ib)%vold_scell
2867c pause
2868
2869 ENDIF
2870 enddo!next IB
2871
2872 !------------------------------------------------------------!
2873 ! @29. COMPUTE VOLUME CHANGE FOR EACH CELL !
2874 !------------------------------------------------------------!
2875 IF(ibug22_tracking>0 .OR. ibug22_tracking==-1.AND.itask==0)print *, "=== TOPOLOGICAL TRACKING ==="
2876 !topology is changing so this is not trivial and must be taken into account
2877 !example : a single secnd hexa can be split into {1 penta + 2 tetra} so
2878 !it is not possible to just make difference between two polyhedra volumes.
2879 nbl1 = nbl
2880 IF(dt1==zero)THEN
2881 nbl1 = 0
2882 DO ib=nbf,nbl
2883 brick_list(nin,ib)%POLY(1:9)%DVOL(1) = zero
2884 brick_list(nin,ib)%DVOL = zero
2885 ENDDO
2886 ENDIF
2887 DO ib=nbf,nbl1
2888 id = brick_list(nin,ib)%ID
2889 mcell = brick_list(nin,ib)%mainID
2890 ncell = brick_list(nin,ib)%NBCUT
2891 icell = 0
2892 m(:,:) = zero
2893 wascut = brick_list(nin,ib)%WasCut
2894 mtn_ = iparg(1,ngb(ib))
2895
2896 IF(wascut>0)THEN !Inevitably M (9, JJ)> Zero, out of the conditional test
2897 !filling relation matrix
2898 DO inod = 1,8
2899 jj = brick_list(nin,ib)%NODE(inod)%WhichCell
2900 ii = brick_list(nin,ib)%NODE(inod)%OLD_WhichCell
2901 IF(ii <0) ii = 1 !just being cut on this cycle so whichcell node was not set on previous cycle
2902 m(ii,jj) = brick_list(nin,ib)%POLY(jj)%Vnew !instead of setting kronecker symbol, volume is directly stored. V=0 <=> delta_ij=0
2903 enddo!next INOD
2904 IF(nbcut>2 .AND. m(9,9)==zero)m(9,9) = brick_list(nin,ib)%POLY(jj)%Vnew
2905 !cas particulier cellule 9 et cellule i (sigma_col>1)
2906 !the 2 cannot be linked because they are not on the same side of the surface.
2907 icell = 9
2908 jj = icell
2909 var = zero
2910 DO ii=1,9
2911 IF(m(ii,jj)>zero)var=var+one
2912 enddo!next II
2913 IF(ncell==0)THEN
2914 ii=maxloc(brick_list(nin,ib)%POLY(1:9)%OLD_Vnew,1)
2915 DO i=1,9
2916 IF (i==ii) cycle
2917 m(i,jj) = zero
2918 ENDDO
2919 ELSE
2920 IF(var>two )THEN
2921 m(1:8,jj) = zero !Several Polyedres disappear, we no longer count them.
2922 ELSEIF(var==two )THEN
2923 !keep the largest volume
2924cc DO II=1,8 ! 1st value is between 1,8 ; last one is at index 9
2925cc VAR = BRICK_LIST(NIN,IB)%POLY(II)%OLD_Vnew
2926cc IF(VAR>ZERO)EXIT
2927cc ENDDO
2928cc IF(VAR>BRICK_LIST(NIN,IB)%POLY(9)%OLD_Vnew)THEN
2929cc M(9,JJ) = ZERO
2930cc ELSE
2931cc M(II,JJ) = ZERO
2932cc ENDIF
2933
2934 !Correction trial April 20, 2010
2935 var = brick_list(nin,ib)%POLY(1)%OLD_Vnew
2936 var_ = brick_list(nin,ib)%POLY(9)%OLD_Vnew
2937 lcond1 = var<var_
2938 var = brick_list(nin,ib)%POLY(1)%Vnew
2939 var_ = brick_list(nin,ib)%POLY(9)%Vnew
2940 lcond2 = var<var_
2941
2942 m(9,9) = zero
2943
2944! IF(lCOND1/=lCOND2)THEN
2945! VAR=M(9,1)
2946! M(9,1) = M(1,9)
2947! M(1,9) = VAR
2948! ENDIF
2949
2950
2951
2952 ENDIF
2953 endif!NCELL==0
2954
2955 ELSEIF(wascut == 0)THEN
2956 mcell = brick_list(nin,ib)%MainID
2957 m(1:9,1:9) = zero
2958 m(1,mcell) = brick_list(nin,ib)%POLY(mcell)%Vnew
2959
2960 ENDIF
2961
2962c IF(IBUG22_TRACKING>0 .OR. IBUG22_TRACKING==-1)then
2963c !display matrix
2964c if (ncell > 0)then
2965c write(*,FMT=('(A,9I15 )')) " ", 1,2,3,4,5,6,7,8,9
2966c do II=1,9
2967c write(*,FMT=('(A,9F15.5 )')) " |", M(II,1:9 )
2968c enddo
2969c endif
2970c ENDIF
2971
2972
2973
2974 !calculating old_volume for each cell
2975 icell = 0
2976 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
2977 icell = icell + 1
2978 IF (icell>ncell .AND. ncell/=0)icell=9
2979 var = zero
2980 var_ = zero
2981 jj = icell
2982 vfrac(1:4) = zero
2983 var_vf(1:4) = zero
2984 vold_phase(1:4) = zero
2985 vj = brick_list(nin,ib)%POLY(jj)%Vnew
2986 DO ii=1,9
2987 IF(m(ii,jj)==zero)cycle
2988 vi = brick_list(nin,ib)%POLY(ii)%OLD_Vnew
2989 !Vi : Vprev, Vj: Vnew for ratio
2990 var = var + vi *vj / sum(m(ii,:)) !this is not a sum of kronecker symbol since volume were stored instead
2991 ENDDO
2992 brick_list(nin,ib)%POLY(icell)%Vold = var
2993
2994c IF(IBUG22_TRACKING==-1 .OR. IBUG22_TRACKING==IXS(11,brick_list(nin,ib)%id))then
2995c write(*,FMT=('(A,I8,A1,I1 )')) " brickID =", ixs(11,id),".",ICELL
2996c write(*,FMT=('(A,I8 )')) " NBCUT =", brick_list(nin,ib)%nbcut
2997c write(*,FMT=('(A,I8 )')) " WASCUT =", brick_list(nin,ib)%wascut
2998c write(*,FMT=('(A,F30.16 )')) " VOLD =", BRICK_LIST(NIN,IB)%POLY(ICELL)%Vold
2999c write(*,FMT=('(A,F30.16 )')) " VNEW =", BRICK_LIST(NIN,IB)%POLY(ICELL)%Vnew
3000c var_ = zero
3001c var__ = zero
3002c icellTag(1:9) = 0
3003c icellTag_old(1:9) = 0
3004c do ii=1,9
3005c if(M(ii,ICELL)>ZERO)then
3006c var_=var_+ONE
3007c icellTag(ii) = 1
3008c endif
3009c enddo
3010c if(var_==ONE)then
3011c var__=zero
3012c do ii=1,9
3013c K = MAXLOC(icellTag,1)
3014c if( M( K ,ii) > ZERO )then
3015c var__=var__+ONE
3016c icelltag_old(ii) = 1
3017c endif
3018c enddo
3019c endif
3020c if ( var_ >= TWO )then
3021c SOM = ZERO
3022c VFRAC(1:4) = ZERO
3023c DO k=1,9
3024c if(icelltag(k)==0)CYCLE
3025c Vi = BRICK_LIST(NIN,IB)%POLY(K)%OLD_Vnew
3026c SOM = SOM + Vi
3027c enddo
3028c write(*,FMT=('(A,I8)')) " merge from poly:", IcellTag(1:9)
3029c elseif (var_ == ONE .AND. var__==ONE)then
3030c write(*,FMT=('(A,I8)')) " cell is stable (same number as previous cycle)", Icelltag(1)
3031c elseif (var__ >= TWO)then
3032c K = MAXLOC(icellTag,1)
3033c write(*,FMT=('(A,I8)')) " cell is a fragment of a previous polyhedron ",K !only one in the column (var_) and on the line (var__)
3034c endif
3035c ENDIF
3036
3037 enddo!next ICELL
3038 ENDDO !next IB
3039
3040
3041
3042 !------------------------------------------------------------!
3043 ! @30. LINK SWITCHING !
3044 ! link has change then transportation must occur !
3045 !------------------------------------------------------------!
3046 nbl1 = nbl
3047 in = 0
3048 debug_outp2 = .false.
3049 if(ibug22_link_switch/=0)then
3050 if(itask==0)ALLOCATE (tmp22array(6,nb))
3051 call my_barrier
3052 tmp22array(1:6,nbf:nbl)=zero
3053 do ib=nbf,nbl
3054 ie = brick_list(nin,ib)%id
3055 if(ixs(11,ie)==ibug22_link_switch .OR. ibug22_link_switch==-1)then
3056 debug_outp2 = .true.
3057 exit
3058 endif
3059 enddo
3060 endif!(IBUG22_LINK_SWITCH/=0)
3061 IF(dt1==zero)nbl1=0
3062 DO ib=nbf,nbl1
3063 ncell = brick_list(nin,ib)%NBCUT
3064 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
3065 !pNAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%NAdjCell
3066 padjbrick => brick_list(nin,ib)%Adjacent_Brick(1:6,1:5)
3067 psubvol(1:9) => brick_list(nin,ib)%POLY(1:9)
3068 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
3069 mcell = brick_list(nin,ib)%MainID
3070 id = brick_list(nin,ib)%ID
3071 ng = brick_list(nin,ib)%NG
3072 idloc = brick_list(nin,ib)%IDLOC
3073 gbuf => elbuf_tab(ngb(ib))%GBUF
3074 lbuf => elbuf_tab(ngb(ib))%BUFLY(1)%LBUF(1,1,1)
3075 icell = 0
3076c print *, "check switch in brickID=", ixs(11,brick_list(nin,ib)%id)
3077 !------LOOP--------!
3078 numsecnd = brick_list(nin,ib)%SecndList%Num
3079 !BRICK_LIST(NIN,IB)%SecndList%ISWITCH(1:NumSECND) = 0
3080 DO ic=1,numsecnd
3081 icellv = brick_list(nin,ib)%SecndList%ICELLv(ic)
3082 ibv = brick_list(nin,ib)%SecndList%IBv(ic)
3083 fm = brick_list(nin,ib)%SecndList%FM(ic)
3084 fv = brick_list(nin,ib)%SecndList%FV(ic)
3085 brickid = brick_list(nin,ib)%ID
3086 ibv = ibv
3087 ngv = brick_list(nin,ibv)%NG
3088 idlocv = brick_list(nin,ibv)%IDLOC
3089 iv = brick_list(nin,ibv)%ID
3090 itaskv = brick_list(nin,ibv)%ITASK
3091 ! loop over the nodes of the adjacent secondary cell
3092 nnodes = brick_list(nin,ibv)%POLY(icellv)%NumNOD
3093 j1 = 0
3094 j2 = 0
3095 DO k1=1,nnodes !Nodes attached to IB
3096 inod = brick_list(nin,ibv)%POLY(icellv)%ListNodID(k1)
3097 fv_old = brick_list(nin,ibv)%NODE(inod)%WhereWasMain !FACE
3098 IF (fv_old == 0)cycle !the node was wet by main cell : MERGE PROCESS APPLIES (see related chapter)
3099 !OLD_M
3100 IF (fv_old == fv)cycle ! NO LINK CHANGE : this is still linked to the same direction
3101 !else FV_old /= FV : THERE WAS A SWITCH
3102 ! GET PREVIOUS FACE TO GET PREVIOUS main CELL
3103 IF(fv_old<=nv46)THEN
3104 ibmo = brick_list(nin,ibv)%Adjacent_Brick(fv_old,4)
3105 ELSE
3106 j = fv_old
3107 j1 = j/10
3108 j2 = mod(j,10)
3109 ibv = brick_list(nin,ibv)%Adjacent_Brick(j1,4)
3110 ibmo = brick_list(nin,ibv)%Adjacent_Brick(j2,4)
3111 ENDIF
3112
3113 ibm = ib
3114 lltm = iparg(2,ngb(ib))
3115 ! IF(IBm==IBmo)CYCLE : impossible it meansFV_old=FV
3116 lcycle = .false.
3117 numtarget = brick_list(nin,ibmo)%Ntarget
3118 IF(numtarget>=2)print *, "**inter22 - Warning Multiple targets",ixs(11,brick_list(nin,ibv)%ID), icellv
3119 DO itarget=1,numtarget
3120 IF (brick_list(nin,ibmo)%MergeTarget(3,itarget)==ibm)THEN
3121 lcycle = .true.
3122 EXIT
3123 ENDIF
3124 ibmo = brick_list(nin,ibmo)%MergeTarget(3,itarget) !In all rigor it would be necessary to take up in all the Target, here we take the first
3125 print *, "**inter22 : multiple targets", ixs(11,brick_list(nin,ibv)%ID), icellv
3126 ENDDO
3127 if(lcycle)cycle
3128 lcycle = .false.
3129 !remove salve cell volume from IBMo(Vold) and add it in IBo(Vold)
3130 !addition is treated now since only current thread is handling this cell IBm \in [NBF,NBL]
3131 !search corresponding secnd cell in IBmo%old_secnds_lit
3132 lfound = .false. !if found old salve cell attached to INOD
3133 numsecnd2 = old_secndlist(nin,ibmo)%Num
3134 ngo = brick_list(nin,ibmo)%NG
3135 llto = iparg(2,ngo)
3136 DO ic2 = 1, numsecnd2
3137 !Get Previous Second Cell
3138 fv2 = old_secndlist(nin,ibmo)%FV(ic2)
3139 icellv2 = old_secndlist(nin,ibmo)%ICELLv(ic2)
3140 ibv2 = old_secndlist(nin,ibmo)%IBv(ic2)
3141 !If the secondary cell was linked by a face other than FV_old then it was not the correct main cell and we pass over
3142 IF(fv2 /= fv_old )cycle
3143 !the secondary cell is attached to the main cell that interests us IBmo
3144 !if we arrive here, it means that the secondary cell (0, in part, or in whole) is changing its link
3145 !we check if one of its nodes corresponds to the INOD node of the new secondary cell in IBv
3146 !we loop over the vertex nodes of the old cell to see which vertex node has switched
3147 volcell = zero
3148 volsecnd = zero
3149 nnodes2 = old_secndlist(nin,ibmo)%NumNOD_Cell(ic2)
3150 DO k2=1,nnodes2
3151 inod2 = old_secndlist(nin,ibmo)%ListNodID(ic2,k2)
3152 icell2 = brick_list(nin,ibv2)%NODE(inod2)%WhichCell
3153 ibmcur = brick_list(nin,ibv2)%POLY(icell2)%WhereIsMain(4)
3154 IF(ibmcur == ib) volsecnd = volsecnd + one/nnodes*brick_list(nin,ibv2)%POLY(icell2)%Vnew !Example Brick 20 becomes main in case 15.14: the Penta in 41 has two nodes, if we complete 2 times we take its volume twice to inject it in Bric 30.
3155 volcell = volcell + one/nnodes*brick_list(nin,ibv2)%POLY(icell2)%Vnew
3156 ENDDO
3157 if (volsecnd == zero)cycle
3158 !BRICK_LIST(NIN,IB)%SecndList%ISWITCH(IC) = 1
3159 !print *, ""
3160 !print *, "VolCell =", VolCell
3161 !print *, "VolSECND =", VolSECND
3162 ratio = volsecnd/volcell !is the extraction ratio in the volume of the old secondary cell.
3163 !here IC2(cycle-1) corresponds to IC(cycle+0)
3164 volcell = ratio * old_secndlist(nin,ibmo)%VOL(ic2)
3165 vold = brick_list(nin,ibm)%Vold_SCell
3166 gbufo => elbuf_tab(ngb(ibmo))%GBUF
3167 !---A.MATERIAL MERGING---!
3168 !Ratio is <1 in case 1 secnd split into several secnd cells which different main links, then collect only relevant quantity not all from old secnd cell.
3169 eint = volcell * gbufo%EINT(idlocb(ibmo))
3170 rho = volcell * gbufo%RHO(idlocb(ibmo))
3171 mom(1) = volcell * gbufo%MOM(llto*(1-1) + idlocb(ibmo))
3172 mom(2) = volcell * gbufo%MOM(llto*(2-1) + idlocb(ibmo))
3173 mom(3) = volcell * gbufo%MOM(llto*(3-1) + idlocb(ibmo))
3174 DO j=1,nv46
3175 sig(j) = volcell * gbufo%SIG(llto*(j-1)+idlocb(ibmo))
3176 ENDDO
3177 gbuf%EINT(idlocb(ibm)) = (gbuf%EINT(idlocb(ibm)) * vold + eint) / (vold+volcell)
3178 gbuf%RHO(idlocb(ibm)) = (gbuf%RHO(idlocb(ibm)) * vold + rho) / (vold+volcell)
3179 gbuf%MOM(lltm*(1-1) + idlocb(ibm)) = (gbuf%MOM(lltm*(1-1) + idlocb(ibm)) * vold + mom(1)) / (vold+volcell)
3180 gbuf%MOM(lltm*(2-1) + idlocb(ibm)) = (gbuf%MOM(lltm*(2-1) + idlocb(ibm)) * vold + mom(2)) / (vold+volcell)
3181 gbuf%MOM(lltm*(3-1) + idlocb(ibm)) = (gbuf%MOM(lltm*(3-1) + idlocb(ibm)) * vold + mom(3)) / (vold+volcell)
3182 alefvm_buffer%FCELL(4, brick_list(nin,ibm)%id )= gbuf%RHO(idlocb(ibm))
3183 alefvm_buffer%FCELL(1, brick_list(nin,ibm)%id )= gbuf%MOM(lltm*(1-1) + idlocb(ibm))
3184 alefvm_buffer%FCELL(2, brick_list(nin,ibm)%id )= gbuf%MOM(lltm*(2-1) + idlocb(ibm))
3185 alefvm_buffer%FCELL(3, brick_list(nin,ibm)%id )= gbuf%MOM(lltm*(3-1) + idlocb(ibm))
3186
3187 !FCELL(5, brick_list(nin,ibm)%id) = ! see mtn 37 conditional block below
3188 !FCELL(6, brick_list(nin,ibm)%id) = ! idem
3189
3190 DO j=1,nv46
3191 gbuf%SIG(lltm*(j-1)+idlocb(ibm)) = (gbuf%SIG(lltm*(j-1)+idlocb(ibm)) * vold + sig(j)) / (vold+volcell)
3192 ENDDO
3193
3194 if(debug_outp2)then
3195 !!!write (*,FMT='(A)') " === LINK SWITCH ==="
3196 in = in + 1
3197 endif
3198
3199 if(ibug22_link_switch==ixs(11,brick_list(nin,ib)%id) .or. ibug22_link_switch==-1 )then
3200 ! print *, "brick target =", ixs(11,brick_list(nin,ib)%id)
3201 ! print *, "brick origin =", ixs(11,brick_list(nin,ibmo)%id)
3202 ! print *, "adding",VOLcell ,'to',BRICK_LIST(NIN,IBm)%Vold_SCell
3203 ! print *, "updated target -old volume- =",BRICK_LIST(NIN,IBm)%Vold_SCell +VOLcell
3204 tmp22array(1,ibm)= ixs(11,brick_list(nin,ib)%id)
3205 tmp22array(2,ibm)= ixs(11,brick_list(nin,ibmo)%id)
3206 tmp22array(3,ibm)= ixs(11,brick_list(nin,ibm)%id)
3207 tmp22array(4,ibm)= brick_list(nin,ibm)%Vold_SCell
3208 tmp22array(5,ibm)= brick_list(nin,ibm)%Vold_SCell +volcell
3209 tmp22array(6,ibm)= volcell
3210 endif
3211
3212 !GBUF%VOL(IDLOCB(IBm)) = GBUF%VOL(IDLOCB(IBm)) + VOLcell
3213 brick_list(nin,ibm)%Vold_SCell = brick_list(nin,ibm)%Vold_SCell + volcell
3214
3215 mtn_ = iparg(1,ngb(ibmo))
3216 volcell51(1:trimat) = zero
3217 IF(mtn_==37)THEN
3218 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
3219 mt = ixs(1,brickid)
3220 iadbuf = ipm(7,mt)
3221 rho10 = bufmat(iadbuf-1+11)
3222 rho20 = bufmat(iadbuf-1+12)
3223 !UVAR(I,1) : massic percentage of liquid * global density (rho1*V1/V : it needs to give liquid mass multiplying by element volume in aleconve.F)
3224 !UVAR(I,2) : density of gas
3225 !UVAR(I,3) : density of liquid
3226 !UVAR(I,4) : volumetric fraction of liquid
3227 !UVAR(I,5) : volumetric fraction of gas
3228 mbuf => elbuf_tab(ngb(ibm))%BUFLY(1)%MAT(1,1,1)
3229 mbufo => elbuf_tab(ngb(ibmo))%BUFLY(1)%MAT(1,1,1)
3230 llt_ = iparg(2,ngb(ibm))
3231 llt_o = iparg(2,ngb(ibmo))
3232
3233 ! UVAR(5) - alpha_gas
3234 pvar => mbuf%VAR((5-1)*llt_ +idlocb(ibm))
3235 pvaro => mbufo%VAR((5-1)*llt_o+idlocb(ibmo))
3236C print *, "link switch 5", pVAR, (pVAR * VOLD + VOLcell * pVARo)/ (VOLD+VOLcell)
3237 pvar = (pvar * vold + volcell * pvaro)/ (vold+volcell)
3238 pvar = max(pvar,zero)
3239
3240 ! Uvar (4) - alpha_liquid
3241 pvar => mbuf%VAR((4-1)*llt_ +idlocb(ibm))
3242 pvaro => mbufo%VAR((4-1)*llt_o+idlocb(ibmo))
3243C print *, "link switch 4", pVAR, (pVAR * VOLD + VOLcell * pVARo)/ (VOLD+VOLcell)
3244 pvar = (pvar * vold + volcell * pvaro)/ (vold+volcell)
3245 pvar = max(pvar,zero)
3246
3247 ! Uvar (3) - Rho_liquid
3248 pvar => mbuf%VAR((3-1)*llt_ +idlocb(ibm))
3249 pvaro => mbufo%VAR((3-1)*llt_o+idlocb(ibmo))
3250 alp = mbuf%VAR((4-1)*llt_ +idlocb(ibm))
3251 alpo = mbufo%VAR((4-1)*llt_o+idlocb(ibmo))
3252C print *, "link switch 3", pVAR, (pVAR*ALP*VOLD+VOLcell*pVARo*ALPo /(VOLD+VOLcell)
3253 pvar = (pvar*alp*vold+volcell*pvaro*alpo)
3254 IF(pvar>zero)pvar=pvar/(alp*vold+alpo*volcell)
3255 pvar = max(pvar,zero)
3256 !cts14june2017!
3257 IF( pvar == zero) pvar = rho10
3258
3259 ! UVAR(2) - rho_gas
3260 pvar => mbuf%VAR((2-1)*llt_ +idlocb(ibm))
3261 pvaro => mbufo%VAR((2-1)*llt_o+idlocb(ibmo))
3262 alp = mbuf%VAR((5-1)*llt_ +idlocb(ibm))
3263 alpo = mbufo%VAR((5-1)*llt_o+idlocb(ibmo))
3264C print *, "link switch 2", pVAR, (pVAR*ALP*VOLD+VOLcell*pVARo*ALPo /(VOLD+VOLcell)
3265 pvar = (pvar*alp*vold+volcell*pvaro*alpo)
3266 IF(pvar>zero)pvar=pvar/(alp*vold+alpo*volcell)
3267 pvar = max(pvar,zero)
3268 !cts14june2017!
3269 IF( pvar == zero) pvar = rho20
3270
3271 ! Uvar (1) - Rho1.v1/V = Rho1.alpha_liquid
3272 pvar => mbuf%VAR((1-1)*llt_ +idlocb(ibm))
3273 pvaro => mbufo%VAR((1-1)*llt_o+idlocb(ibmo))
3274C print *, "link switch 1", pVAR, (pVAR * VOLD + VOLcell * pVARo)/ (VOLD+VOLcell)
3275 pvar = (pvar * vold + volcell * pvaro)/ (vold+volcell)
3276 pvar = max(pvar,zero)
3277
3278 !GBUF%RHO(IDLOCB(IBm)) = (GBUF%RHO(IDLOCB(IBm)) * VOLD + RHO) / (VOLD+VOLcell)
3279
3280
3281
3282 !APPEL SIGEPS37 POUR EQUILIBRE PRESSION
3283 brickid = brick_list(nin,ibm)%ID
3284 idloc = idlocb(ibm)
3285 ng = ngb(ibm)
3286 vol = vold+volcell
3287 CALL sigeps37_single_cell (
3288 1 elbuf_tab, ixs, bufmat, iparg, ipm,
3289 2 idloc , ng , brickid, vol
3290 . )
3291
3292
3293 ELSEIF(mtn_==51)THEN
3294 mbuf => elbuf_tab(ngb(ibm))%BUFLY(1)%MAT(1,1,1)
3295 llt_ = iparg(2,ngb(ibm))
3296 mbufo => elbuf_tab(ngb(ibmo))%BUFLY(1)%MAT(1,1,1)
3297 llt_o = iparg(2,ngb(ibmo))
3298
3299 DO itrimat=1,trimat
3300 !IBm_old - volume fraction!
3301 ipos = 1
3302 k0 = m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1 ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
3303 k = k0 * llt_o
3304 vfrac(itrimat) = mbufo%VAR(k+idlocb(ibmo))
3305 enddo!next ITRIMAT
3306
3307 DO itrimat=1,trimat
3308 ipos = 11
3309 k0 = m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1 ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
3310 k = k0 * llt_
3311 volcell51(itrimat) = volcell * vfrac(itrimat)
3312 mbuf%VAR(k+idlocb(ibm)) = mbuf%VAR(k+idlocb(ibm)) + volcell51(itrimat)
3313 enddo!next ITRIMAT
3314
3315 DO itrimat=1,trimat
3316 !IBm - volume fraction!
3317 ipos = 1
3318 k0 = m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1 ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
3319 k = k0 * llt_
3320 vfrac(itrimat) = mbuf%VAR(k+idlocb(ibm))
3321 IF(volcell51(itrimat)<=zero)cycle !unplug if nothiong to add (otherwise div by zero)
3322 DO ipos = 1,m51_nvphas
3323 IF(ipos==11)cycle !volume already treated
3324 k2 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
3325 k = llt_ * k2
3326 ko = llt_o * k2
3327 pvar => mbuf%VAR(k+idlocb(ibm))
3328 pvar = ( pvar * vfrac(itrimat)*vold + volcell51(itrimat)*mbufo%VAR(ko+idlocb(ibmo)) )
3329 pvar = pvar / (vfrac(itrimat)*vold + volcell51(itrimat))
3330 enddo!next IPOS
3331 enddo!next ITRIMAT
3332
3333 endif! law37 | law51
3334 !---THREAD TAG---!
3335C print *, " --- set VOLD_L += VOLcell, IB, ID ", VOLcell , IBMo, ixs(11,brick_list(nin,ibmo)%ID)
3336 vold_l(itask,0,ibmo) = vold_l(itask,0,ibmo) + volcell
3337 vold_l(itask,1:trimat,ibmo) = vold_l(itask,1:trimat,ibmo) + volcell51(1:trimat)
3338 enddo!next IC2
3339 EXIT !next IC (since this one is treated)
3340 enddo!next K1
3341 enddo!next IC
3342 enddo!next IB
3343
3344 CALL my_barrier
3345
3346 if(ibug22_link_switch/=0 .AND. itask==0 .and. dt1>zero)THEN
3347 write (*,fmt='(A)') " === LINK SWITCH ==="
3348 DO ib=1,nb
3349 IF(tmp22array(1,ib)==0)cycle
3350 print *, "brick target =", tmp22array(1,ib)
3351 print *, "brick origin =", tmp22array(2,ib)
3352 print *, "brick main =", tmp22array(3,ib)
3353 print *, "adding",tmp22array(6,ib) ,'to', tmp22array(4,ib)
3354 print *, "updated target -old volume- =", tmp22array(5,ib)
3355 ENDDO
3356 endif
3357 call my_barrier
3358
3359 !print
3360 if(itask==0)then
3361 if(ibug22_link_switch/=0)then
3362 DO ib=1,nb
3363 DO it = 0,nthread-1
3364 IF (vold_l(it,0,ib) == zero)cycle
3365 print *, ""
3366 print *, " brick ID =", ixs(11,brick_list(nin,ib)%id)
3367 print *, " removing ",vold_l(it,0,ib),'from',brick_list(nin,ib)%Vold_SCell
3368 print *, " new origin volume =", ixs(11,brick_list(nin,ib)%id)
3369 print *, " %vold, %vnew " , brick_list(nin,ib)%Vold_SCell- vold_l(it,0,ib), brick_list(nin,ib)%Vnew_SCell
3370 print *, ""
3371 ENDDO
3372 ENDDO
3373 endif
3374 endif
3375 call my_barrier
3376
3377 !---B.OLD VOLUME UPDATE---!
3378 !substract only on current thread
3379 DO ib=nbf,nbl1
3380 !GBUF => ELBUF_TAB(NGB(IB))%GBUF
3381 DO it = 0,nthread-1
3382 IF (vold_l(it,0,ib) == zero)cycle
3383 if(ibug22_link_switch==ixs(11,brick_list(nin,ib)%id) .or. ibug22_link_switch==-1 )then
3384 !print *, ""
3385 !print *, " brick ID =", ixs(11,brick_list(nin,ib)%id)
3386 !print *, " removing ",VOLD_L(IT,0,IB),'from',BRICK_LIST(NIN,IB)%Vold_SCell
3387 !print *, " new origin volume =", ixs(11,brick_list(nin,ib)%id)
3388 !print *, " %vold, %vnew " , BRICK_LIST(NIN,IB)%Vold_SCell- VOLD_L(IT,0,IB), BRICK_LIST(NIN,IB)%Vnew_SCell
3389 !print *, ""
3390 endif
3391
3392 brick_list(nin,ib)%Vold_SCell = brick_list(nin,ib)%Vold_SCell - vold_l(it,0,ib)
3393
3394 !!!!!!!!!!!!BRICK_LIST(NIN,IB)%Vold_SCell = BRICK_LIST(NIN,IB)%Vold_SCell - VOLD_L(IT,0,IB)
3395 mtn_ = iparg(1,ngb(ib))
3396 IF(mtn_==51)THEN
3397 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
3398 llt_ = iparg(2,ngb(ib))
3399 DO itrimat=1,trimat
3400 ipos = 11
3401 k0 = m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1 ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
3402 k = k0 * llt_
3403 mbuf%VAR(k+idlocb(ib)) = mbuf%VAR(k+idlocb(ib)) - vold_l(it,itrimat,ib)
3404 enddo!next ITRIMAT
3405 endif!(MTN_==51)
3406
3407
3408ccc if(IBUG22_LINK_SWITCH==ixs(11,brick_list(nin,ib)%id) .or. IBUG22_LINK_SWITCH==-1 )then
3409ccc print *, " --reset VOLD_L = ZERO, IB, ID ", ZERO , IB, ixs(11,brick_list(nin,ib)%ID)
3410ccc endif
3411 vold_l(it,0:max(0,trimat),ib) = zero !reset VOLD_L for next cycle
3412 enddo!next IT
3413 enddo!next IB
3414
3415
3416
3417
3418 if(ibug22_link_switch/=0)THEN
3419 call my_barrier
3420 if(itask==0)DEALLOCATE (tmp22array)
3421 endif
3422
3423 !------------------------------------------------------------!
3424 ! @32. LEVEL2 NEIGHBORHOOD (FOR FLUXES & CONVECTION) !
3425 !------------------------------------------------------------!
3426 DO ib=nbf,nbl
3427 nbcut = brick_list(nin,ib)%NBCUT
3428 pismain(1:9) => brick_list(nin,ib)%POLY(1:9)
3429 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
3430 padjbrick => brick_list(nin,ib)%Adjacent_Brick(1:6,1:5)
3431 pmainid => brick_list(nin,ib)%MainID
3432 ncell = nbcut
3433 icell = 0
3434 ie = brick_list(nin,ib)%ID
3435 !we handle the main cells to avoid lockon/lockoff (by avoiding loops on the secondarys)
3436 mcell = pmainid
3437 IF(mcell == 0)cycle ! next hexahedron cell
3438 brick_list(nin,ib)%ADJ_ELEMS%Num = 0
3439 DO j=1,6
3440 iv = padjbrick(j,1)
3441 ngv = padjbrick(j,2)
3442 idlocv = padjbrick(j,3)
3443 ibv = padjbrick(j,4)
3444 IF(ibv>0)THEN !if adjacent brick to face J is in cut cell buffer
3445 nbcutv = brick_list(nin,ibv)%NBCUT
3446 IF(nbcutv==0)cycle
3447 nadjcell = brick_list(nin,ib)%POLY(mcell)%FACE(j)%NAdjCell
3448 DO ic=1,nadjcell
3449 icellv = brick_list(nin,ib)%POLY(mcell)%FACE(j)%Adjacent_Cell(ic)
3450 mcellv = brick_list(nin,ibv)%MainID
3451 IF(icellv == mcellv)cycle
3452 !this adjacent cell is a secnd one, use it as bridge to a new adjacent main cell.
3453 !pWhereIsMainV => BRICK_LIST(NIN,IBv)%POLY(1:9)%WhereIsMain(1:4)
3454 jv = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(1) ! pWhereIsMainV(ICELLv,1)
3455 idm = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(3) ! pWhereIsMainV(ICELLv,3)
3456 IF(ie/=idm)THEN !secnd cell merged to the main IB
3457 !take all main neighbors of the secondary cell except JV face
3458 !add only one neighbor the one with which the secondary cell is merged
3459 brick_list(nin,ib)%ADJ_ELEMS%Num = brick_list(nin,ib)%ADJ_ELEMS%Num + 1
3460 k = brick_list(nin,ib)%ADJ_ELEMS%Num
3461 brick_list(nin,ib)%ADJ_ELEMS%IV(k) = idm
3462 brick_list(nin,ib)%ADJ_ELEMS%IB(k) = ibm
3463 brick_list(nin,ib)%ADJ_ELEMS%ICELL(k) = icellv
3464 brick_list(nin,ib)%ADJ_ELEMS%SecndFACE(k) = jv !to retrieve cut face area
3465 else! !secnd cell merged to main /= IB
3466 padjbrickv => brick_list(nin,ibv)%Adjacent_Brick(1:6,1:5)
3467 DO j2=1,6 !loop on IV.icellv (secnd) adjacent bricks
3468 IF (j2==jv)cycle !behind JV there is IB itself, so skip it
3469 ivv = padjbrick(j2,1)
3470 ngvv = padjbrick(j2,2)
3471 idlocvv = padjbrick(j2,3)
3472 ibvv = padjbrick(j2,4)
3473 ifvv = padjbrick(j2,5)
3474 IF(ivv == 0)cycle
3475
3476 IF(ibvv > 0)THEN
3477 mcellvv = brick_list(nin,ibvv)%MainID
3478 !if main IVV.MCELvv has connexion with IV.ICELLv then add it
3479 nadjcellv = brick_list(nin,ibvv)%POLY(mcellvv)%FACE(ifvv)%NAdjCell
3480 DO k=1,nadjcellv
3481 icv = brick_list(nin,ibvv)%POLY(mcellvv)%FACE(ifvv)%Adjacent_Cell(k)
3482 IF(icv/=icellv)cycle
3483 brick_list(nin,ib)%ADJ_ELEMS%Num = brick_list(nin,ib)%ADJ_ELEMS%Num + 1
3484 k2 = brick_list(nin,ib)%ADJ_ELEMS%Num
3485 brick_list(nin,ib)%ADJ_ELEMS%IV(k2) = ivv
3486 brick_list(nin,ib)%ADJ_ELEMS%IB(k2) = ibvv
3487 brick_list(nin,ib)%ADJ_ELEMS%ICELL(k2) = icellv !=ICv
3488 brick_list(nin,ib)%ADJ_ELEMS%SecndFACE(k2) = j2 !to retrieve cut face area
3489 EXIT
3490 enddo!next K
3491 ELSEIF(ibvv == 0)THEN
3492 brick_list(nin,ib)%ADJ_ELEMS%Num = brick_list(nin,ib)%ADJ_ELEMS%Num + 1
3493 k2 = brick_list(nin,ib)%ADJ_ELEMS%Num
3494 brick_list(nin,ib)%ADJ_ELEMS%IV(k2) = ivv
3495 brick_list(nin,ib)%ADJ_ELEMS%IB(k2) = ibvv
3496 brick_list(nin,ib)%ADJ_ELEMS%ICELL(k2) = 1
3497 brick_list(nin,ib)%ADJ_ELEMS%SecndFACE(k2) = j2 !to retrieve cut face area
3498 ENDIF
3499 enddo!next J2
3500 endif!(IE==IDM)
3501 enddo!next IC
3502 ENDIF
3503 enddo!next J
3504 ENDDO
3505
3506
3507
3508 !------------------------------------------------------------!
3509 ! @32. COUNT POINTS ON A GIVEN POLYHEDRON FACE !
3510 !------------------------------------------------------------!
3511 DO ib=nbf,nbl
3512 nbcut = brick_list(nin,ib)%NBCUT
3513 mcell = brick_list(nin,ib)%MainID
3514 ncell = brick_list(nin,ib)%NBCUT
3515 icell = 0
3516 nintp(1:6,1:9) = 0
3517 nnod(1:6,1:9) = 0
3518 IF(nbcut == 0) THEN
3519 brick_list(nin,ib)%POLY(1)%FACE(1:6)%NumPOINT=(/4,4,4,4,4,4/)
3520 cycle
3521 ENDIF
3522 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
3523 icell = icell + 1
3524 IF (icell>ncell .AND. ncell/=0)THEN
3525 icell=9
3526 cycle
3527 ENDIF
3528 secid = brick_list(nin,ib)%SECID_Cell(icell)
3529 brick_list(nin,ib)%POLY(icell)%FACE(1:6)%NumPOINT = 0
3530 DO j=1,6
3531 jj = gface(j,secid)
3532 IF(jj==0)EXIT
3533 np = gnpt(j,secid)
3534 nn = gnnod(j,secid)
3535 brick_list(nin,ib)%POLY(icell)%FACE(jj)%NumPOINT = np
3536 nintp(jj,icell) = np-nn !only intersection point
3537 nnod(jj,icell) = nn
3538 ENDDO
3539 enddo!next ICELL
3540 IF(icell==9)THEN
3541 DO j=1,6
3542 !sommets
3543 nn = sum( nnod(j,1:ncell) )
3544 nn = 4 - nn !complete
3545 !intersection points
3546 np = sum( nintp(j,1:ncell) )
3547 !Total point number = summits + intersection PT
3548 brick_list(nin,ib)%POLY(9)%FACE(j)%NumPOINT = np + nn
3549 ENDDO
3550 ENDIF
3551 ENDDO
3552
3553 !------------------------------------------------------------!
3554 ! @33. LOCATE WHERE WAS main FOR NEXT CYCLE !
3555 ! & TAG ITS BRICK NODES (update in i22intersect.F) !
3556 !------------------------------------------------------------!
3557 DO ib=nbf,nbl
3558 pwherewasmain(1:8) => brick_list(nin,ib)%NODE(1:8)!%WhereWasMain
3559 !pWhereIsMain => BRICK_LIST(NIN,IB)%POLY(1:9)%WhereIsMain(1:4)
3560 plistnodid(1)%p(1:8) => brick_list(nin,ib)%POLY(1)%ListNodID(1:8)
3561 plistnodid(2)%p(1:8) => brick_list(nin,ib)%POLY(2)%ListNodID(1:8)
3562 plistnodid(3)%p(1:8) => brick_list(nin,ib)%POLY(3)%ListNodID(1:8)
3563 plistnodid(4)%p(1:8) => brick_list(nin,ib)%POLY(4)%ListNodID(1:8)
3564 plistnodid(5)%p(1:8) => brick_list(nin,ib)%POLY(5)%ListNodID(1:8)
3565 plistnodid(6)%p(1:8) => brick_list(nin,ib)%POLY(6)%ListNodID(1:8)
3566 plistnodid(7)%p(1:8) => brick_list(nin,ib)%POLY(7)%ListNodID(1:8)
3567 plistnodid(8)%p(1:8) => brick_list(nin,ib)%POLY(8)%ListNodID(1:8)
3568 plistnodid(9)%p(1:8) => brick_list(nin,ib)%POLY(9)%ListNodID(1:8)
3569 pnodwasmain(1:8) => brick_list(nin,ib)%NODE(1:8)!%NodWasMain
3570 ncell = brick_list(nin,ib)%NBCUT
3571 icell = 0
3572 mcell = brick_list(nin,ib)%MainID
3573 ! nodes wet by main cell
3574 mnod = brick_list(nin,ib)%POLY(mcell)%NumNOD
3575 pnodwasmain(1:8)%NodWasMain = 0
3576 pwherewasmain(1:8)%WhereWasMain = 0
3577 DO j=1,mnod
3578 pnodwasmain(plistnodid(mcell)%p(j))%NodWasMain = 1
3579 ENDDO
3580 ! locate main cell for nodes of a given secnd cell
3581 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
3582 icell = icell +1
3583 IF (icell>ncell .AND. ncell/=0)icell=9
3584 IF(icell == mcell)cycle
3585 ipos = brick_list(nin,ib)%POLY(icell)%WhereIsMain(1) !pWhereIsMain(ICELL,1)
3586 mnod = brick_list(nin,ib)%POLY(icell)%NumNOD
3587 !All polyhedron nodes have the same main cell : face=IPOS
3588 DO j=1,mnod
3589 pwherewasmain(plistnodid(icell)%p(j))%WhereWasMain = ipos
3590 enddo!next J
3591 enddo!next ICELL
3592 ENDDO
3593
3594
3595
3596 !------------------------------------------------------------!
3597 ! @34. SET IIAD22(1:NUMELS) = %ELBUG_TAB(NG)%TAG22(I) !
3598 !------------------------------------------------------------!
3599 DO ib=nbf,nbl
3600 ie = brick_list(nin,ib)%ID
3601 iiad22(nin, ie) = ib
3602 ENDDO
3603
3604 !------------------------------------------------------------!
3605 ! @35. SET MLW !
3606 !------------------------------------------------------------!
3607 DO ib=nbf,nbl
3608 mlw = iparg(1,ngb(ib))
3609 brick_list(nin,ib)%MLW = mlw
3610 ENDDO
3611
3612 !------------------------------------------------------------!
3613 ! @36. RESET OLD SECND LIST CONNECTIVITY !
3614 ! WILL BE DEFINED IN I22_GET_PRV_DATA() !
3615 !------------------------------------------------------------!
3616 !RESET ONLY NON ZERO VALUES
3617 DO ib=nbf,nbl
3618 numsecnd = old_secndlist(nin,ib)%Num
3619 IF (numsecnd==0)cycle
3620 old_secndlist(nin,ib)%Num = 0
3621 old_secndlist(nin,ib)%NumSecndNodes = 0
3622 DO j=1,24
3623 old_secndlist(nin,ib)%FM(j) = 0
3624 old_secndlist(nin,ib)%FV(j) = 0
3625 old_secndlist(nin,ib)%IV(j) = 0
3626 old_secndlist(nin,ib)%IBV(j) = 0
3627 old_secndlist(nin,ib)%ICELLv(j) = 0
3628 old_secndlist(nin,ib)%VOL(j) = zero
3629 old_secndlist(nin,ib)%NumNOD_Cell(j) = 0
3630 old_secndlist(nin,ib)%ListNodID(j,1:8) = 0
3631 ENDDO
3632 enddo!next IB
3633
3634 !------------------------------------------------------------!
3635 ! @37. SET main STRONG NODE !
3636 !------------------------------------------------------------!
3637 DO ib=nbf,nbl
3638 mcell = brick_list(nin,ib)%mainID
3639 CALL weighting_cell_nodes(nin,ib,mcell,inod_w, idemerge)
3640 IF(mcell==0)THEN
3641 brick_list(nin,ib)%OldMainStrongNode = 0
3642 ELSE
3643 brick_list(nin,ib)%OldMainStrongNode = inod_w
3644 ENDIF
3645 ENDDO
3646
3647 !------------------------------------------------------------!
3648 ! @38. DRAW SUPERCELL CONNEXION WITH 1D TRUSSES !
3649 !------------------------------------------------------------!
3650 lstilltruss = .true.
3651 igr = ipari(52) !reserved GRTRUS ---> truss group identifier
3652 ntrus = ipari(53 )!reserved GRTRUS ---> truss group NEL
3653 IF( itask==0 .AND. ntrus/=0 )THEN
3654 ii = 1 ! first node of group node
3655 DO ib=nbf,nbl
3656 mcell = brick_list(nin,ib)%mainID
3657 IF (mcell==0 ) cycle
3658 IF (.NOT.lstilltruss) cycle
3659 point0(1:3) = brick_list(nin,ib)%POLY(mcell)%CellCenter(1:3)
3660
3661 !LOOP ON SECND CELLS TO DRAW CONNEXION
3662 numsecnd = brick_list(nin,ib)%SecndList%Num
3663 DO isecnd=1,numsecnd
3664 ibv = brick_list(nin,ib)%SecndList%IBV(isecnd)
3665 icellv = brick_list(nin,ib)%SecndList%ICELLv(isecnd)
3666 pointtmp(1:3) = brick_list(nin,ibv)%POLY(icellv)%CellCenter(1:3)
3667 !check if group has still one truss
3668 IF (ii>ntrus) THEN
3669 lstilltruss=.false.
3670 EXIT
3671 print *, "** Warning inter22 : no more truss in group to mark cell links"
3672 ENDIF
3673 x(1:3,ixt(2,igrtruss(igr)%ENTITY(ii))) = point0(1:3)
3674 x(1:3,ixt(3,igrtruss(igr)%ENTITY(ii))) = pointtmp(1:3)
3675 !write (*,*) "set TRUS_id=", IXT(NIXT,IGRTRUSS(IGR)%ENTITY(II)) ," corr=(", Point0(1:3), ")-(" ,PointTMP(1:3) ,")"
3676 if(ibug22_truss==1)then
3677 write (*,fmt='(A,I10,A,I10,A,I10)') "set TRUS_id=", ixt(nixt,igrtruss(igr)%ENTITY(ii)) ,
3678 . " main=", ixs(11,brick_list(nin,ib)%id) ," secnd=", ixs(11,brick_list(nin,ibv)%id)
3679 endif
3680 ii = ii + 1 !next truss
3681 ENDDO
3682 ENDDO !next IB
3683 DO ii = 1,ntrus
3684 !print *, "reset TRUS_id=", IXT(NIXT,IGRTRUSS(IGR)%ENTITY(II))
3685 x(1:3,ixt(2,igrtruss(igr)%ENTITY(ii))) = (/zero, zero, zero/)
3686 x(1:3,ixt(3,igrtruss(igr)%ENTITY(ii))) = (/ one, zero, zero/)
3687 ENDDO
3688 ENDIF !ITASK==0
3689
3690
3691 !------------------------------------------------------------!
3692 ! @39. COMPUTE DVOL PREDICTION FOR EACH CELL !
3693 !------------------------------------------------------------!
3694 DO ib=nbf,nbl
3695 nbcut = brick_list(nin,ib)%NBCUT
3696 icell = 0
3697 ncell = nbcut
3698 IF(nbcut == 0)THEN
3699 icell = 1
3700 brick_list(nin,ib)%POLY(icell)%DVOL(1) = zero
3701 cycle
3702 ENDIF
3703 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
3704 icell = icell +1
3705 IF (icell>ncell .AND. ncell/=0)icell=9
3706 IF(icell == 9)THEN
3707 brick_list(nin,ib)%POLY(9)%DVOL(1) = -sum(brick_list(nin,ib)%POLY(1:nbcut)%DVOL(1))
3708 ELSE
3709 ! V_face
3710 vsum(1) = brick_list(nin,ib)%PCUT(icell)%Vel(1)
3711 vsum(2) = brick_list(nin,ib)%PCUT(icell)%Vel(2)
3712 vsum(3) = brick_list(nin,ib)%PCUT(icell)%Vel(3)
3713 !VSUM(1) = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%VEL(1)
3714 !VSUM(2) = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%VEL(2)
3715 !VSUM(3) = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%VEL(3)
3716 !not unitary : NN= 2Sn
3717 n_(1) = brick_list(nin,ib)%PCUT(icell)%N(1)
3718 n_(2) = brick_list(nin,ib)%PCUT(icell)%N(2)
3719 n_(3) = brick_list(nin,ib)%PCUT(icell)%N(3)
3720 dvol_predic = dt1*(vsum(1)*n_(1) + vsum(2)*n_(2) + vsum(3)*n_(3))
3721 brick_list(nin,ib)%POLY(icell)%DVOL(1) = dvol_predic !for next CYCLE
3722 ENDIF
3723 enddo!next ICELL
3724 enddo!next IB
3725 !!!!!!!!
3726 !output!
3727 !!!!!!!!
3728 ! IF(IBUG22_PREDICTION/=0)THEN
3729 ! print *, "==PREDICTION DVOL"
3730 ! DO IB=NBF,NBL
3731 ! IF(IBUG22_PREDICTION>0 .AND. IBUG22_PREDICTION/=IXS(11,brick_list(nin,ib)%id))cycle
3732 ! NCELL = BRICK_LIST(NIN,IB)%NBCUT
3733 ! ICELL = 0
3734 ! DO WHILE (ICELL<=NCELL) ! loop on polyhedron {1:NCELL} U {9}
3735 ! ICELL = ICELL +1
3736 ! IF (ICELL>NCELL .AND. NCELL/=0)ICELL = 9
3737 ! write(*,FMT=('(A,I8,A,I1,A,F30.16)'))"brickID=", ixs(11,brick_list(nin,ib)%id),".",ICELL,"=",
3738 ! . BRICK_LIST(NIN,IB)%POLY(ICELL)%DVOL(1)
3739 ! ENDDO!next ICELL
3740 ! ENDDO!next IB
3741 ! ENDIF
3742 !------------------------------------------------------------!
3743 ! @40. BUILD SUPERCELL DVOL (PREDICTION) !
3744 !------------------------------------------------------------!
3745 !--------------!
3746 CALL my_barrier
3747 !--------------!
3748 !THEN COMPUTE SUPERCELL GEOMETRICAL DVOL
3749 DO ib=nbf,nbl
3750 ncell = brick_list(nin,ib)%NBCUT
3751 mcell = brick_list(nin,ib)%MainID
3752 brick_list(nin,ib)%DVOL = brick_list(nin,ib)%POLY(mcell)%DVOL(1)
3753 !------LOOP--------!
3754 numsecnd = brick_list(nin,ib)%SecndList%Num
3755 DO ic=1,numsecnd
3756 icellv = brick_list(nin,ib)%SecndList%ICELLv(ic)
3757 ibv = brick_list(nin,ib)%SecndList%IBv(ic)
3758 ngv = brick_list(nin,ibv)%NG
3759 idlocv = brick_list(nin,ibv)%IDLOC
3760 iv = brick_list(nin,ibv)%ID !'0' PREDICTION with old conformation
3761 brick_list(nin,ib)%DVOL = brick_list(nin,ib)%DVOL + brick_list(nin,ibv)%POLY(icellv)%DVOL(1) !'1' PREDICTION with new conformation
3762 enddo!next IC
3763 enddo!next IB
3764
3765
3766 !Please note that Treaty identical in law 51 for v1old, v2old ..
3767
3768!!!!!!!!!!!!!!!!!!!!!!!! attention en modifiant ici les vold_scell, on perturbe le critre predictif ci-dessous.
3769
3770
3771
3772
3773 !------------------------------------------------------------!
3774 ! @41. FIX SUPERCELL DVOL (CORRECTION) !
3775 !------------------------------------------------------------!
3776
3777 lskip = .false.
3778
3779c if(IBUG22_dvol /= 0)then
3780c print *, "+---Sinit22_fvm : DVOL majoration"
3781c endif
3782 if(dt1==zero)then
3783 do ib=nbf,nbl
3784 brick_list(nin,ib)%Vold_SCell = brick_list(nin,ib)%Vnew_SCell
3785 enddo
3786 endif
3787CCC RATIO22 = 1.02D00
3788
3789 if(.NOT.lskip)then
3790 DO ib=nbf,1*nbl
3791 icode = brick_list(nin,ib)%ICODE
3792 old_icode = brick_list(nin,ib)%OLD_ICODE
3793 if (brick_list(nin,ib)%Vold_SCell==-111.d00) brick_list(nin,ib)%Vold_SCell = brick_list(nin,ib)%Vnew_SCell !never init.
3794 vnew = brick_list(nin,ib)%Vnew_SCell
3795 vold = brick_list(nin,ib)%Vold_SCell
3796 dvol_numeric = vnew-vold
3797 dvol_predic = brick_list(nin,ib)%dvol
3798 uncutvol = brick_list(nin,ib)%UncutVOL
3799 !IF GEOMETRICAL VALUE IS RISEN ABOVE ITS PREDICTION : CORRECT
3800c if((IBUG22_prediction ==-1 .OR. IBUG22_prediction==ixs(11,brick_list(nin,ib)%id)) )then
3801c print *, "+"
3802c print *, "+------elem_id =",ixs(11,brick_list(nin,ib)%id)
3803c print *, "+--------old_icode =",OLD_ICODE
3804c print *, "+--------icode =",ICODE
3805c print *, "+--------dvol_prediction =",brick_list(nin,ib)%dvol
3806c print *, "+--------dvol_numeric =",VNEW-VOLD
3807c print *, "+--------vnew =",VNEW
3808c print *, "+--------ratio22 =",RATIO22
3809c print *, "+--------DVOL_NUM/V =",ABS(DVOL_NUMERIC)/UncutVOL
3810c endif
3811
3812 IF(abs(dvol_numeric) > ratio22*abs(dvol_predic) .AND. dvol_predic /= zero .AND. ratio22 < 1000 )THEN
3813 IF((icode /= old_icode ) )THEN
3814! IF( (ICODE /= OLD_ICODE ) .OR. (ICODE==OLD_ICODE>ZERO .AND. ABS(DVOL_NUMERIC/DVOL_PREDIC)>1.5d00 ))THEN
3815c if((IBUG22_prediction ==-1 .OR. IBUG22_prediction==ixs(11,brick_list(nin,ib)%id)) )then
3816c print *, "+"
3817c print *, "+------elem_id =",ixs(11,brick_list(nin,ib)%id)
3818c print *, "+--------old_icode =",OLD_ICODE
3819c print *, "+--------icode =",ICODE
3820c print *, "+--------dvol_prediction =",brick_list(nin,ib)%dvol
3821c print *, "+--------dvol_numeric =",VNEW-VOLD
3822c print *, "+--------vnew =",VNEW
3823c endif
3824 brick_list(nin,ib)%Vold_SCell = brick_list(nin,ib)%Vnew_SCell - dvol_predic
3825 ng = brick_list(nin,ib)%NG
3826 idloc = brick_list(nin,ib)%IDLOC
3827 gbuf => elbuf_tab(ng)%GBUF
3828C write (*,FMT='(A,F30.16)' ) " UncutVOL =", brick_list(nin,ib)%uncutvol
3829C write (*,FMT='(A,F30.16)' ) " ABS(DVOL_NUMERIC)/UncutVOL =",ABS(DVOL_NUMERIC)/UncutVOL
3830 !################################################################!
3831 !# TRAITEMENT MULTIMATERIAU LAW51 #!
3832 !################################################################!
3833 mlw = brick_list(nin,ib)%MLW
3834 ng = brick_list(nin,ib)%NG
3835 idloc = brick_list(nin,ib)%IDLOC
3836 IF(mlw==51)THEN
3837 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
3838 llt = iparg(2,ng)
3839 DO itrimat=1,trimat
3840 !IPOS=01 : fractions volumiques
3841 !IPos = 11: old volumes of each phase
3842 ipos = 1
3843 k = llt * (m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1)
3844 vfrac(itrimat) = mbuf%VAR(k+idloc)
3845 ipos = 11
3846 k = llt * (m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1)
3847 volcell51_old(itrimat) = mbuf%VAR(k+idloc)
3848 mbuf%VAR(k+idloc) = max(zero, (vnew-dvol_predic)*vfrac(itrimat) )
3849 volcell51(itrimat) = mbuf%VAR(k+idloc)
3850 enddo!next ITRIMAT
3851 ENDIF
3852 if(ibug22_prediction ==-1 .OR. ibug22_prediction==ixs(11,brick_list(nin,ib)%id))then
3853c print *, "+--------vold =",VOLD ,"->", BRICK_LIST(NIN,IB)%Vold_SCell
3854! print *, "+--------vold_1 =",VolCELL51_OLD(1) ,"->", VolCELL51(1)
3855! print *, "+--------vold_2 =",VolCELL51_OLD(2) ,"->", VolCELL51(2)
3856! print *, "+--------vold_3 =",VolCELL51_OLD(3) ,"->", VolCELL51(3)
3857! print *, "+--------vold_4 =",VolCELL51_OLD(4) ,"->", VolCELL51(4)
3858 endif
3859 endif!IF(ICODE /= OLD_ICODE .AND. OLD_ICODE /= 0)
3860 endif!IF(ABS(DVOL_NUMERIC) > RATIO22*ABS(DVOL_PREDIC) .AND. RATIO22 < 1000 .AND. ABS(DVOL_NUMERIC)/UncutVOL>0.05)
3861 enddo!next IB
3862 endif
3863
3864 IF( (ibug22_prediction ==-1 ))THEN
3865 call my_barrier
3866 IF(itask==0)THEN
3867 DO ib=1,nb
3868 icode = brick_list(nin,ib)%ICODE
3869 old_icode = brick_list(nin,ib)%OLD_ICODE
3870 vnew = brick_list(nin,ib)%Vnew_SCell
3871 vold = brick_list(nin,ib)%Vold_SCell
3872 dvol_numeric = vnew-vold
3873 dvol_predic = brick_list(nin,ib)%dvol
3874 uncutvol = brick_list(nin,ib)%UncutVOL
3875 print *, "+------elem_id =",ixs(11,brick_list(nin,ib)%id)
3876 print *, "+--------old_icode =",old_icode
3877 print *, "+--------icode =",icode
3878 print *, "+--------dvol_prediction =",brick_list(nin,ib)%dvol
3879 print *, "+--------dvol_numeric =",vnew-vold
3880 print *, "+--------vnew =",vnew
3881 print *, "+--------vold =",vold ,"->", brick_list(nin,ib)%Vold_SCell
3882 ENDDO
3883 ENDIF
3884 ENDIF
3885
3886
3887 !------------------------------------------------------------!
3888 ! @42. TIME ZERO VOLUME INIT (NOT DONE IN STARTER) !
3889 ! needed for convection subroutines !
3890 ! sforc3>voln22 not called before convection !
3891 !------------------------------------------------------------!
3892 IF(dt1==zero)THEN
3893 DO ib=nbf,nbl
3894 psubvol(1:9) => brick_list(nin,ib)%POLY(1:9)
3895 psubvold => brick_list(nin,ib)%Vold_SCell
3896 mcell = brick_list(nin,ib)%MainID
3897 IF(mcell==0)cycle
3898 psubvold = brick_list(nin,ib)%Vnew_Scell
3899 ENDDO
3900 ENDIF
3901 DO ib=nbf,nbl
3902 psubvol(1:9) => brick_list(nin,ib)%POLY(1:9)
3903 psubvold => brick_list(nin,ib)%Vold_SCell
3904 mcell = brick_list(nin,ib)%MainID
3905 !IF(MCELL==0)CYCLE
3906 gbuf => elbuf_tab(ngb(ib))%GBUF
3907 IF(psubvold>zero)gbuf%VOL(idlocb(ib)) = psubvold
3908C write (*,FMT='(A,I10,A,F30.16,A,F30.16)')"brickID=", IXS(11,brick_list(nin,ib)%id),
3909C . " vold=",brick_list(nin,ib)%vold_scell," vnew=",brick_list(nin,ib)%vnew_scell
3910 ENDDO
3911
3912
3913 !------------------------------------------------------------!
3914 ! @43. DEBUG OUTPUT !
3915 !------------------------------------------------------------!
3916 if (debug_outp .AND. ibug22_sinit/=0)then
3917 do ib=nbf,nbl
3918 psubvol(1:9) => brick_list(nin,ib)%POLY(1:9)
3919 plistnodid(1)%p(1:8) => brick_list(nin,ib)%POLY(1)%ListNodID(1:8)
3920 plistnodid(2)%p(1:8) => brick_list(nin,ib)%POLY(2)%ListNodID(1:8)
3921 plistnodid(3)%p(1:8) => brick_list(nin,ib)%POLY(3)%ListNodID(1:8)
3922 plistnodid(4)%p(1:8) => brick_list(nin,ib)%POLY(4)%ListNodID(1:8)
3923 plistnodid(5)%p(1:8) => brick_list(nin,ib)%POLY(5)%ListNodID(1:8)
3924 plistnodid(6)%p(1:8) => brick_list(nin,ib)%POLY(6)%ListNodID(1:8)
3925 plistnodid(7)%p(1:8) => brick_list(nin,ib)%POLY(7)%ListNodID(1:8)
3926 plistnodid(8)%p(1:8) => brick_list(nin,ib)%POLY(8)%ListNodID(1:8)
3927 plistnodid(9)%p(1:8) => brick_list(nin,ib)%POLY(9)%ListNodID(1:8)
3928 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
3929 !pNAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%NAdjCell
3930 padjbrick => brick_list(nin,ib)%Adjacent_Brick(1:6,1:5)
3931 pismain(1:9) => brick_list(nin,ib)%POLY(1:9)
3932 pmainid => brick_list(nin,ib)%MainID
3933 id = brick_list(nin,ib)%ID
3934 icell = 0
3935 mcell = mcell
3936 ncell = brick_list(nin,ib)%NBCUT
3937 gbuf => elbuf_tab(ngb(ib))%GBUF
3938 vol = gbuf%VOL(idlocb(ib))
3939 brickid = idb(ib)
3940 if(ibug22_sinit/=ixs(11,brickid) .AND. ibug22_sinit/=-1) cycle
3941 write (*,fmt='(A,I12)') "+=== BRICK ID===", ixs(11,id)
3942 if(ncell==0)then
3943 write (*,fmt='(A )') "| uncut: "
3944 write (*,fmt='(A,1F30.20)') "| volume: ", vol
3945 write (*,fmt='(A,1F30.20)') "| ext. volume: ", brick_list(nin,ib)%Vnew_Scell
3946 write (*,fmt='(A,1F30.20)') "| masse: ", gbuf%VOL(idlocb(ib)) * gbuf%RHO(idlocb(ib))
3947 if(brick_list(nin,ib)%SecndList%Num>0)then
3948 do j=1,brick_list(nin,ib)%SecndList%Num
3949 write (*,fmt='(A )') "| secnd list : "
3950 write (*,fmt='(A,I10 )') "| + J : ", brick_list(nin,ib)%SecndList%FM(j)
3951 write (*,fmt='(A,I10 )') "| + IB : ", brick_list(nin,ib)%SecndList%IBv(j)
3952 write (*,fmt='(A,I10 )') "| + brickID : ", ixs(11,brick_list(nin,ib)%SecndList%IV(j))
3953 enddo
3954 endif
3955 cycle
3956 endif
3957 write (*,fmt='(A,1F30.20)') "| volume: ", vol
3958 write (*,fmt='(A,6F30.20)') "| faces: ", f(1:6,ib)
3959 write (*,fmt='(A,1F30.20)') "| masse: ", gbuf%VOL(idlocb(ib)) * gbuf%RHO(idlocb(ib))
3960 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
3961 icell = icell +1
3962 IF (icell>ncell .AND. ncell/=0)icell=9
3963 debugmainsecnd = '.........|'
3964 mnod = brick_list(nin,ib)%POLY(icell)%NumNOD
3965 write (*,fmt='(A )') "|"
3966 if(icell/=9)then
3967 write (*,fmt='(A,I1,A,A,A1,I2,A,A6)')
3968 . "+== ICELL= ", icell , ", SecType=", brick_list(nin,ib)%SECTYPE(icell) ,
3969 . "(", brick_list(nin,ib)%SecID_Cell(icell) , ") - ", char1
3970 else
3971 write (*,fmt='(A,A6)') "+== REMAINING POLYHEDRON - ", char1
3972 endif
3973
3974 write (*,fmt='(A )') "| |"
3975 write (*,fmt='(A,I1)') "| +===Main/Secnd=", pismain(icell)%IsMain
3976 write (*,fmt='(A,F30.20)') "| +======SUVOLUMES=", psubvol(icell)%Vnew
3977 write (*,fmt='(A,6F30.20)') "| +=======SUBFACES=", brick_list(nin,ib)%POLY(icell)%FACE(1:6)%Surf
3978 write (*,fmt='(A,F30.20)') "| +=======CUT AERA=", brick_list(nin,ib)%PCUT(icell)%SCUT(1)
3979 write (*,fmt='(A,A,I2)') "| +======NUM POINT=", " ",brick_list(nin,ib)%POLY(icell)%NumPOINT
3980 write (*,fmt='(A,A,I1,A,8I12)') "| +======NODE LIST=",
3981 . " (",mnod,") ", plistnodid(icell)%p(1:mnod)
3982 write (*,fmt='(A,A,8I12)') "| | radIDs=",
3983 . " ", ixs(1+plistnodid(icell)%p(1:mnod),id)
3984 write (*,fmt='(A,A,8I12)') "| | userIDs=",
3985 . " ", itab(ixs(1+plistnodid(icell)%p(1:mnod),id))
3986 iad2 = ale_connectivity%ee_connect%iad_connect(brickid)
3987 lgth2 = ale_connectivity%ee_connect%iad_connect(brickid+1) -
3988 . ale_connectivity%ee_connect%iad_connect(brickid)
3989 If(sum(iabs(ale_connectivity%ee_connect%connected(iad2:iad2 + lgth2 - 1)))/=0)then
3990 write (*,fmt='(A,6I10)') "| +===Adj Brick(i)=", padjbrick(1:6,1)
3991 DO j=1,6
3992 IF( padjbrick(j,1)/=0 )THEN
3993 write (*,fmt='(A,6I10)') "| +===Adj Brick(u)=", ixs(11,padjbrick(j,1))
3994 ENDIF
3995 ENDDO
3996 DO j=1,6
3997 nadjcell = brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell
3998 IF(nadjcell/=0)THEN
3999 write (*,fmt='(A,I1,A,5I3)') "| +======Adj Cells, face=",j, " :",
4000 . brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(1:nadjcell)
4001 ENDIF
4002 ENDDO
4003 ENDIF
4004 enddo!next icell
4005 write (*,fmt='(A )') " "
4006 enddo!next IB
4007 if(itask==0.AND.debug_outp)then
4008 write (*,fmt='(A )') " ----sini22_end----"
4009 write (*,fmt='(A )') " "
4010 endif
4011 endif! (IBUG22_sinit/=0)
4012
4013 !------------------------------------------------------------!
4014 ! @44. DEBUG - WRITE CUT CELL BUFFER !
4015 !------------------------------------------------------------!
4016 ! CALL WRITE_CUT_CELL_BUFFER() !post/debug
4017
4018 !------------------------------------------------------------!
4019 ! @45. SUPERCELL CENTERS !
4020 !------------------------------------------------------------!
4021 DO ib=nbf,nbl
4022 ncell = brick_list(nin,ib)%NBCUT
4023 numsecnd = brick_list(nin,ib)%SecndList%Num
4024 mainid = brick_list(nin,ib)%mainID
4025 IF (mainid ==0) cycle
4026 volcell = brick_list(nin,ib)%POLY(mainid)%Vnew
4027 vol = volcell !cumul
4028 ppoint(1:3) = brick_list(nin,ib)%POLY(mainid)%CellCenter(1:3) * volcell
4029 DO ic=1,numsecnd
4030 icellv = brick_list(nin,ib)%SecndList%ICELLv(ic)
4031 ibv = brick_list(nin,ib)%SecndList%IBv(ic)
4032 volcell = brick_list(nin,ibv)%POLY(icellv)%Vnew
4033 point0(1:3) = brick_list(nin,ibv)%POLY(icellv)%CellCenter(1:3)
4034 ppoint(1) = ppoint(1) + volcell * point0(1)
4035 ppoint(2) = ppoint(2) + volcell * point0(2)
4036 ppoint(3) = ppoint(3) + volcell * point0(3)
4037 vol = vol + volcell
4038 enddo!next IC
4039 ppoint(1) = ppoint(1) / vol
4040 ppoint(2) = ppoint(2) / vol
4041 ppoint(3) = ppoint(3) / vol
4042 brick_list(nin,ib)%SCellCenter(1:3) = ppoint(1:3)
4043 enddo!next IB
4044
4045 !------------------------------------------------------------!
4046 ! @46. MARK SUPER-CELL CENTERS WITH ORPHAN NODES !
4047 !------------------------------------------------------------!
4048 lstillnode = .true.
4049 IF(ipari(70)/=0)THEN
4050 nnodes = igrnod(ipari(70))%NENTITY
4051 !!!WRITE(*,*)(ITAB(IBUFSSG(J)),J=IAD0,IAD0+NNODES-1)
4052 IF( itask==0 .AND. nnodes/=0 )THEN
4053 ii = 1 ! first node of group node
4054 DO ib=nbf,nbl
4055 mcell = brick_list(nin,ib)%mainID
4056 IF (mcell==0 ) cycle
4057 IF (.NOT.lstillnode) cycle
4058 point0(1:3) = brick_list(nin,ib)%SCellCenter(1:3)
4059 IF (ii>nnodes) THEN
4060 lstillnode = .false.
4061 print *, "** Warning inter22 : no more node in group to mark cell center. Last one was" ,
4062 . itab(igrnod(ipari(70))%ENTITY(nnodes))
4063 EXIT
4064 ENDIF
4065 x(1:3,igrnod(ipari(70))%ENTITY(ii)) = point0(1:3)
4066 if(ibug22_orphannodes == 1)then
4067 write (*,fmt='(A,I10,A,I10,A,I10)')"set orphan_node_id=",itab(igrnod(ipari(70))%ENTITY(ii))
4068 endif
4069 ii = ii + 1 !next node
4070 ENDDO !next IB
4071 DO ii = 1, nnodes
4072 !print *, "reset SPHCEL_id=", IXT(NIXT,IGRNOD(IPARI(70))%ENTITY(II))
4073 x(1:3,igrnod(ipari(70))%ENTITY(ii)) = (/zero, zero, zero/)
4074 ENDDO
4075 ENDIF !ITASK==0
4076 ENDIF
4077
4078 !------------------------------------------------------------!
4079 ! @47. UNCUT CELLS + POLY 9 : CENTERS !
4080 !------------------------------------------------------------!
4081 !FOR UNCUT BRICKS, Centers for polyehedra faces computed in i22subvol
4082 DO ib=nbf,nbl
4083 ncell = brick_list(nin,ib)%NBCUT
4084 IF(ncell==0)THEN
4085 ie = brick_list(nin,ib)%ID
4086 nc(1:8) = ixs(2:9,ie)
4087 DO j=1,6
4088 brick_list(nin,ib)%POLY(1)%FACE(j)%Center(1) = fourth * sum( x(1, nc(icf(1:4,j)) ) )
4089 brick_list(nin,ib)%POLY(1)%FACE(j)%Center(2) = fourth * sum( x(2, nc(icf(1:4,j)) ) )
4090 brick_list(nin,ib)%POLY(1)%FACE(j)%Center(3) = fourth * sum( x(3, nc(icf(1:4,j)) ) )
4091 ENDDO
4092 ELSE
4093 !simplification, car approximation precendente invalide lorsque limit(FACE9) = 0
4094 ie = brick_list(nin,ib)%ID
4095 nc(1:8) = ixs(2:9,ie)
4096 !
4097 !ICELL 9 only since polyhedra face centers were computed in i22subvol.
4098 ! Poly9 is complementary polyhedron, it is built without graph but deduced by boolean operation from full brick.
4099 ! its face centers must be computed to display face center (velocity post-treatment)
4100 ! face center are only known for cut polyhedra (A & C below), center for poly9 must be deduced (poly B below)
4101 !
4102 ! N4 N3
4103 ! !------------+-------!
4104 ! ! | I3 !
4105 ! ! | !
4106 ! ! | ! Ni : nodes
4107 ! ! | ! Ii : intersection points
4108 ! ! I4 B | C ! A,C : cut polyhedra
4109 ! !\ | ! B : complemntary polyhedron
4110 ! ! \ | ! Ca,Cb,Cc : face centers for each A,B,C polyhedra
4111 ! ! \ | ! Npa,Npb,Npc : number of points for polygon on face
4112 ! ! A \ I1 | I2 !
4113 ! !----+-------|-------!
4114 ! N1 N2
4115 !
4116 ! Ca = (N1+I1+I4)/3
4117 ! Cb = (N4+I1+I2+I3+I4)/5
4118 ! Cc = (N2+N3+I2+I3)/4
4119 ! 3Ca + 5Cb + 4Cc = Sum(Ni) + 2 Sum(Ii)
4120 ! => poly9 C9 = [-Npa.Ca -Npc.Cc + Sum(Ni) +2.Sum(Ii)] / Np9
4121 !
4122 !
4123 Do j=1,6
4124 face = brick_list(nin,ib)%POLY(9)%FACE(j)%Surf
4125 np_(9) = brick_list(nin,ib)%POLY(9)%FACE(j)%NumPOINT
4126 center(:,:) = zero
4127 IF(abs(face)<=em10 .OR. np_(9)==0)THEN
4128 brick_list(nin,ib)%POLY(9)%FACE(j)%Center(1:3) = zero
4129 ELSE
4130 ncell = brick_list(nin,ib)%NBCUT
4131 DO icell = 1 ,ncell
4132 np_(icell) = brick_list(nin,ib)%POLY(icell)%FACE(j)%NumPOINT
4133 center(1:3,icell) = brick_list(nin,ib)%POLY(icell)%FACE(j)%Center(1:3)
4134 ENDDO
4135 !NP_ is number of points (intersection + nodes)
4136 ! pPoint is sum of nodes
4137 !CUT_point is sum of intersection points
4138 ! Center : are face center
4139 ! Point0 : is center of poly 9
4140 np_(ncell+1:9) = 0
4141 np_(9) = brick_list(nin,ib)%POLY(9)%FACE(j)%NumPOINT
4142 ppoint(1) = sum( x(1, nc(icf(1:4,j)) ) )
4143 ppoint(2) = sum( x(2, nc(icf(1:4,j)) ) )
4144 ppoint(3) = sum( x(3, nc(icf(1:4,j)) ) )
4145 cut_point(1:3) = brick_list(nin,ib)%POLY(9)%FACE(j)%Center(1:3) !stored there in i22subvol.F it is SUM(cutpoints on the face J)
4146 point0(1) = ppoint(1) + two*cut_point(1)
4147 . - np_(1)*center(1,1)- np_(2)*center(1,2)- np_(3)*center(1,3)- np_(4)*center(1,4)
4148 . - np_(5)*center(1,5)- np_(6)*center(1,6)- np_(7)*center(1,7)- np_(8)*center(1,8)
4149 point0(2) = ppoint(2) + two*cut_point(2)
4150 . - np_(1)*center(2,1)- np_(2)*center(2,2)- np_(3)*center(2,3)- np_(4)*center(2,4)
4151 . - np_(5)*center(2,5)- np_(6)*center(2,6)- np_(7)*center(2,7)- np_(8)*center(2,8)
4152 point0(3) = ppoint(3) + two*cut_point(3)
4153 . - np_(1)*center(3,1)- np_(2)*center(3,2)- np_(3)*center(3,3)- np_(4)*center(3,4)
4154 . - np_(5)*center(3,5)- np_(6)*center(3,6)- np_(7)*center(3,7)- np_(8)*center(3,8)
4155 point0(1:3) = point0(1:3) / np_(9)
4156 brick_list(nin,ib)%POLY(9)%FACE(j)%Center(1:3) = point0(1:3)
4157 !print *, "poly9 center, ixs(11,IE), face J ", ixs(11,IE), J
4158 !print *, BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%Center(1:3)
4159 ENDIF
4160 enddo!next J
4161 ENDIF
4162 enddo!next IB
4163
4164
4165
4166 !------------------------------------------------------------!
4167 ! @48. MARK CELL CENTERS WITH ORPHAN NODES !
4168 !------------------------------------------------------------!
4169 lstillnode = .true.
4170 IF(ipari(81)/=0)THEN
4171 nnodes = igrnod(ipari(81))%NENTITY
4172 !!!WRITE(*,*)(ITAB(IBUFSSG(J)),J=IAD0,IAD0+NNODES-1)
4173 IF( itask==0 .AND. nnodes/=0 )THEN
4174 ii = 1 ! first node of group node
4175 DO ib=nbf,nbl
4176 icell = 0
4177 ncell = brick_list(nin,ib)%NBCUT
4178 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
4179 icell = icell +1
4180 IF (icell>ncell .AND. ncell/=0)icell=9
4181 IF (.NOT.lstillnode) cycle
4182 point0(1:3) = brick_list(nin,ib)%POLY(icell)%CellCenter(1:3)
4183 IF(ii>nnodes)THEN
4184 lstillnode = .false.
4185 print *, "** Warning inter22 : no more node in group to mark cell center",
4186 . itab(igrnod(ipari(81))%ENTITY(nnodes))
4187 EXIT
4188 ENDIF
4189 brick_list(nin,ib)%POLY(icell)%ID_FREE_NODE = igrnod(ipari(81))%ENTITY(ii)
4190 x(1:3,igrnod(ipari(81))%ENTITY(ii)) = point0(1:3)
4191 if(ibug22_orphannodes == 1)then
4192 write (*,fmt='(A,I10,A,I10,A,I10)')"set orphan_node_id=",
4193 . itab(igrnod(ipari(81))%ENTITY(ii)),"in brick_id=",ixs(11,brick_list(nin,ib)%id)
4194 endif
4195 ii = ii + 1 !next node
4196 ENDDO
4197 ENDDO !next IB
4198 DO ii = 1,nnodes
4199 !print *, "reset SPHCEL_id=", IXT(NIXT,IGRNOD(IPARI(81))%ENTITY(II))
4200 x(1:3,igrnod(ipari(81))%ENTITY(ii)) = zero
4201 ENDDO
4202 ENDIF !ITASK==0
4203 ENDIF
4204
4205
4206 !------------------------------------------------------------!
4207 ! @49. MARK FACE CENTERS WITH ORPHAN NODES !
4208 !------------------------------------------------------------!
4209 IF(ipari(82)/=0)THEN
4210 nnodes = igrnod(ipari(82))%NENTITY
4211 IF(nnodes==0)RETURN
4212 ii = 1 ! first node of group node
4213 DO ib=nbf,nbl
4214 ie = brick_list(nin,ib)%ID
4215 icell = 0
4216 ncell = brick_list(nin,ib)%NBCUT
4217 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
4218 icell = icell +1
4219 IF (icell>ncell .AND. ncell/=0)icell=9
4220 IF(.NOT.lstillnode) cycle
4221 DO j=1, 6
4222 IF (ii>nnodes) THEN
4223 lstillnode = .false.
4224 print *, "** Warning inter22 : no more node in group to mark face centers." ,
4225 . itab(igrnod(ipari(82))%ENTITY(nnodes))
4226 EXIT
4227 ENDIF
4228 node_id = igrnod(ipari(82))%ENTITY(ii)
4229 x(1:3,node_id) = brick_list(nin,ib)%POLY(icell)%FACE(j)%Center(1:3)
4230 ii = ii + 1 !next node
4231 ENDDO
4232 ENDDO !next ICELL
4233 enddo!next IB
4234 ENDIF
4235
4236
4237
4238 !print *, "SINIT22_FVM : NB=", NBL
4239 !do ib=nbf,nbl
4240 ! id = brick_list(nin,ib)%id
4241 ! print *, "ixs =", ixs(11,id)
4242 ! print *, "nbcut = ", brick_list(1,ib)%nbcut
4243 ! print *, "sectype1 =", brick_list(1,ib)%sectype(1)
4244 ! print *, "sectype2 =", brick_list(1,ib)%sectype(2)
4245 ! print *, "---------------------"
4246 !enddo
4247
4248
4249 !------------------------------------------------------------!
4250 ! @50. DEALLOCATE !
4251 !------------------------------------------------------------!
4252 IF(ALLOCATED(debugmainsecndv))DEALLOCATE (debugmainsecndv)
4253 IF(ALLOCATED(ismainv))DEALLOCATE (ismainv)
4254 IF(ALLOCATED(f))DEALLOCATE (f)
4255 IF(ALLOCATED(vol51))DEALLOCATE (vol51,vol51v)
4256 IF(ALLOCATED(origin_data))DEALLOCATE (origin_data)
4257 IF(ALLOCATED(destroy_data))DEALLOCATE (destroy_data)
4258 IF(ALLOCATED(norigin))DEALLOCATE (norigin)
4259
4260
4261 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
destroy_cell
subroutine destroy_cell(nin, ib, icell_target, icellv, ibv, j, jv, ixs, itask)
Definition destroy_cell.F:33
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
get_group_id
subroutine get_group_id(ii, ng, ig, iparg)
Definition get_group_id.F:30
get_unique_main_cell
integer function get_unique_main_cell(nin, ib, k)
Definition get_unique_master_cell.F:35
i22aera
function i22aera(npts, p, c)
Definition i22subvol.F:2382
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
id
initmumps id
Definition multiplerhs_example.m:4
ale_connectivity_mod
Definition ale_connectivity_mod.F:225
alefvm_mod
Definition alefvm_mod.F:88
alefvm_mod::alefvm_buffer
type(alefvm_buffer_), target alefvm_buffer
Definition alefvm_mod.F:120
groupdef_mod
Definition groupdef_mod.F:668
i22bufbric_mod
Definition cut-cell-search_mod.F:95
i22bufbric_mod::brick_list
type(brick_entity), dimension(:,:), allocatable, target brick_list
Definition cut-cell-search_mod.F:105
i22bufbric_mod::int22_buf
type(int22_buf_t) int22_buf
Definition cut-cell-search_mod.F:130
i22tri_mod
Definition cut-cell-search_mod.F:212
i22tri_mod::ibug22_truss
integer ibug22_truss
Definition cut-cell-search_mod.F:372
i22tri_mod::nb
integer nb
Definition cut-cell-search_mod.F:487
i22tri_mod::ibug22_ccbuflist
integer ibug22_ccbuflist
Definition cut-cell-search_mod.F:428
i22tri_mod::ibug22_orphannodes
integer ibug22_orphannodes
Definition cut-cell-search_mod.F:373
i22tri_mod::ibug22_destroy
integer ibug22_destroy
Definition cut-cell-search_mod.F:347
i22tri_mod::ibug22_dtmin
integer ibug22_dtmin
Definition cut-cell-search_mod.F:383
i22tri_mod::ibug22_undirectlink
integer ibug22_undirectlink
Definition cut-cell-search_mod.F:357
i22tri_mod::in22
integer in22
Definition cut-cell-search_mod.F:520
i22tri_mod::ibug22_link_switch
integer ibug22_link_switch
Definition cut-cell-search_mod.F:408
i22tri_mod::ibug22_prediction
integer ibug22_prediction
Definition cut-cell-search_mod.F:418
i22tri_mod::ibug22_merge
integer ibug22_merge
Definition cut-cell-search_mod.F:352
i22tri_mod::ibug22_tracking
integer ibug22_tracking
Definition cut-cell-search_mod.F:413
i22tri_mod::ibug22_sinit
integer ibug22_sinit
Definition cut-cell-search_mod.F:316
initbuf_mod
Definition initbuf.F:160
initbuf_mod::initbuf
subroutine initbuf(iparg, ng, mtn, llt, nft, iad, ity, npt, jale, ismstr, jeul, jtur, jthe, jlag, jmult, jhbe, jivf, mid, jpor, jcvt, jclose, jpla, irep, iint, igtyp, israt, isrot, icsen, isorth, isorthg, ifailure, jsms)
Definition initbuf.F:261
sigeps37_single_cell
subroutine sigeps37_single_cell(elbuf_tab, ixs, bufmat, iparg, ipm, idloc, ng, brickid, vol)
Definition sigeps37_single_cell.F:37
my_barrier
subroutine my_barrier
Definition machine.F:31
weighting_cell_nodes
subroutine weighting_cell_nodes(nin, ib, icell, ires, idemerge)
Definition weighting_cell_nodes.F:32