OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ini_inimap1d.F File Reference
#include "implicit_f.inc"
#include "param_c.inc"
#include "scr17_c.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "mvsiz_p.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine ini_inimap1d (inimap1d, elbuf_tab, ipart, iparg, iparts, ipartq, xgrid, vel, ixs, ixq, ixtg, pm, ipm, bufmat, multi_fvm, pld, npc, igrbric, igrquad, igrsh3n, npts, mat_param, snpc, stf)

Function/Subroutine Documentation

◆ ini_inimap1d()

subroutine ini_inimap1d ( type(inimap1d_struct), dimension(ninimap1d), intent(inout), target inimap1d,
type(elbuf_struct_), dimension(ngroup), intent(inout), target elbuf_tab,
integer, dimension(lipart1, *), intent(in) ipart,
integer, dimension(nparg, ngroup), intent(in) iparg,
integer, dimension(*), intent(in) iparts,
integer, dimension(*), intent(in) ipartq,
dimension(3, *), intent(in) xgrid,
dimension(3, *), intent(inout) vel,
integer, dimension(nixs, numels), intent(in) ixs,
integer, dimension(nixq, numelq), intent(in) ixq,
integer, dimension(nixtg, numeltg), intent(in) ixtg,
dimension(npropm, nummat), intent(in) pm,
integer, dimension(npropmi, *), intent(in) ipm,
dimension(*), intent(in) bufmat,
type(multi_fvm_struct), intent(inout) multi_fvm,
dimension(2, npts/2), target pld,
integer, dimension(*), intent(in) npc,
type (group_), dimension(ngrbric) igrbric,
type (group_), dimension(ngrquad) igrquad,
type (group_), dimension(ngrsh3n) igrsh3n,
integer, intent(in) npts,
type(matparam_struct_), dimension(nummat), intent(in) mat_param,
integer, intent(in) snpc,
integer, intent(in) stf )
Parameters
[in]stfarray sizes

Definition at line 35 of file ini_inimap1d.F.

40C-----------------------------------------------
41C M o d u l e s
42C-----------------------------------------------
43 USE elbufdef_mod
44 USE inimap1d_mod
45 USE multi_fvm_mod
46 USE groupdef_mod
47 USE message_mod
48 USE alefvm_mod , only:alefvm_param
49 USE matparam_def_mod
50 USE multimat_param_mod , ONLY : m51_n0phas, m51_nvphas
51 USE multi_solve_eint_mod , ONLY : multi_solve_eint
52C-----------------------------------------------
53C I m p l i c i t T y p e s
54C-----------------------------------------------
55#include "implicit_f.inc"
56C-----------------------------------------------
57C C o m m o n B l o c k s
58C-----------------------------------------------
59#include "param_c.inc"
60#include "scr17_c.inc"
61! NGROUP
62#include "com01_c.inc"
63! NFUNCT
64#include "com04_c.inc"
65! MVSIZ
66#include "mvsiz_p.inc"
67C-----------------------------------------------
68C D u m m y A r g u m e n t s
69C-----------------------------------------------
70 INTEGER,INTENT(IN) :: SNPC, STF !< array sizes
71 TYPE(INIMAP1D_STRUCT), TARGET, DIMENSION(NINIMAP1D), INTENT(INOUT) :: INIMAP1D
72 TYPE(ELBUF_STRUCT_), DIMENSION(NGROUP), INTENT(INOUT), TARGET :: ELBUF_TAB
73 INTEGER, INTENT(IN) :: IPART(LIPART1, *), NPC(*)
74 INTEGER, INTENT(IN) :: IPARTS(*), IPARTQ(*),IPM(NPROPMI, *),IPARG(NPARG, NGROUP),
75 . IXS(NIXS, NUMELS), IXQ(NIXQ, NUMELQ), IXTG(NIXTG, NUMELTG),NPTS
76 my_real, INTENT(IN) :: xgrid(3, *), pm(npropm, nummat), bufmat(*)
77 my_real, INTENT(IN), TARGET :: pld(2, npts/2)
78 my_real, INTENT(INOUT) :: vel(3, *)
79 TYPE(MULTI_FVM_STRUCT), INTENT(INOUT) :: MULTI_FVM
80 TYPE(MATPARAM_STRUCT_) ,DIMENSION(NUMMAT), INTENT(IN) :: MAT_PARAM
81C-----------------------------------------------
82 TYPE (GROUP_) , DIMENSION(NGRBRIC) :: IGRBRIC
83 TYPE (GROUP_) , DIMENSION(NGRQUAD) :: IGRQUAD
84 TYPE (GROUP_) , DIMENSION(NGRSH3N) :: IGRSH3N
85C-----------------------------------------------
86C L o c a l V a r i a b l e s
87C-----------------------------------------------
88 my_real, TARGET :: pld_inv(npts/2,2)
89 INTEGER :: JJ, KK, IAD, NELEM, NG, NFT, NEL, IFIRST, ILAST, IAD1, IAD2, IMID, FUNC,
90 . FUNC1, FUNC2, MATID, MLW
91 LOGICAL :: CONT
92 INTEGER :: II, IFUNC, NPT, FIRST, LAST, ICOOR, INODE, NODEID, IND, SHIFT,
93 . NODE1, NODE2, NODE3, NODE4, NODE5, NODE6, NODE7, NODE8,
94 . ELEMID
95 my_real :: x0, y0, z0, radius, xc, yc, zc, value1, value2, nx, ny, nz, VALUE,
96 . xnode, ynode, znode, rad1, rad2, x1, y1, z1, xp, yp, zp
97 INTEGER, DIMENSION(:), ALLOCATABLE :: ELEM_LIST
98 INTEGER :: LOCAL_ELEM_LIST(MVSIZ)
99 my_real, DIMENSION(:, :), ALLOCATABLE :: pres, eint
100 my_real :: tt, res, rho
101 TYPE(G_BUFEL_), POINTER :: GBUF
102 INTEGER :: NBNODE, NODELIST(8), ITY, ISOLNOD,
103 . GRBRICID,GRQUADID,GRSH3NID, NBMAT, NBMAT_TARGET, IMAT, KPHASE, NUVAR, IADBUF, NPAR, IFORM,
104 . M51_SUBMAT_ID(4),TAG_MAT(NUMMAT)
105 my_real, POINTER, DIMENSION(:) :: ptrx, ptry
106
107 TYPE pseudo_pointer_array
108 my_real, POINTER, DIMENSION(:) :: temp
109 END TYPE
110 TYPE (PSEUDO_POINTER_ARRAY) :: SUBMAT(21)
111
112 CHARACTER*2 :: Str1, Str2
113 my_real :: vfrac(mvsiz,21),densfrac
114
115 TYPE(BUF_EOS_), POINTER :: EBUF
116 INTEGER NVAREOS
117 INTEGER GMID !< global material identifier (multimaterial law)
118 INTEGER :: NVARTMP_EOS
119
120C-----------------------------------------------
121C S o u r c e L i n e s
122C-----------------------------------------------
123
124 DO kk=1,npts/2
125 pld_inv(kk,1)=pld(1,kk)
126 pld_inv(kk,2)=pld(2,kk)
127 ENDDO
128
129 IF (n2d == 0) THEN
130 ALLOCATE(elem_list(numels))
131 elem_list(1:numels)=0
132 ELSE
133 ALLOCATE(elem_list(numelq + numeltg))
134 elem_list(1:numelq+numeltg)=0
135 ENDIF
136
137 tag_mat(1:nummat)=0
138
139 grbricid = 0
140 grquadid = 0
141 grsh3nid = 0
142
143 DO kk = 1, ninimap1d
144 IF(.NOT.inimap1d(kk)%CORRECTLY_READ)cycle
145 nbmat = inimap1d(kk)%NBMAT
146 IF(nbmat==0)cycle !something went wrong with reader
147 ALLOCATE(pres(mvsiz, nbmat), eint(mvsiz, nbmat))
148 IF (inimap1d(kk)%GRBRICID /= 0) THEN
149 grbricid = inimap1d(kk)%GRBRICID
150 ELSEIF (inimap1d(kk)%GRQUADID /= 0) THEN
151 grquadid = inimap1d(kk)%GRQUADID
152 ELSEIF (inimap1d(kk)%GRSH3NID /= 0) THEN
153 grsh3nid = inimap1d(kk)%GRSH3NID
154 ENDIF
155 x0 = xgrid(1, inimap1d(kk)%NODEID1)
156 y0 = xgrid(2, inimap1d(kk)%NODEID1)
157 z0 = xgrid(3, inimap1d(kk)%NODEID1)
158 IF (inimap1d(kk)%NODEID2 > 0) THEN
159 x1 = xgrid(1, inimap1d(kk)%NODEID2)
160 y1 = xgrid(2, inimap1d(kk)%NODEID2)
161 z1 = xgrid(3, inimap1d(kk)%NODEID2)
162 ENDIF
163 IF (grbricid > 0) THEN
164 nelem = igrbric(grbricid)%NENTITY
165 DO ii = 1, nelem
166 elem_list(ii) = igrbric(grbricid)%ENTITY(ii)
167 ENDDO
168 ELSEIF (grquadid > 0) THEN
169 nelem = igrquad(grquadid)%NENTITY
170 DO ii = 1, nelem
171 elem_list(ii) = igrquad(grquadid)%ENTITY(ii)
172 ENDDO
173 ELSEIF (grsh3nid > 0) THEN
174 nelem = igrsh3n(grsh3nid)%NENTITY
175 DO ii = 1, nelem
176 elem_list(ii) = igrsh3n(grsh3nid)%ENTITY(ii)
177 ENDDO
178 ENDIF
179C If NSPMD != 1 the ids are not ordered
180 CALL quicksort_i(elem_list(1:nelem), 1, nelem)
181 DO ng = 1, ngroup
182 mlw = iparg(1, ng)
183 nel = iparg(2, ng)
184 nft = iparg(3, ng)
185 ity = iparg(5, ng)
186 SELECT CASE(mlw)
187 CASE (3, 4, 6, 49, 51, 151)
188 CONTINUE
189 CASE DEFAULT
190 cycle
191 END SELECT
192 matid=0
193 IF (ity == 1) THEN
194 matid = ixs(1, nft + 1)
195 ELSE IF (ity == 2) THEN
196 matid = ixq(1, nft + 1)
197 ELSE IF (ity == 7) THEN
198 matid = ixtg(1, nft + 1)
199 ELSE
200 cycle ! not compatible elem type (beam, truss, spring, shell ...)
201 ENDIF
202 m51_submat_id(1:4) = 0
203 IF (mlw == 51) THEN
204 iform = ipm(62, matid)
205 IF(iform>=2.AND.iform<=6)cycle
206 m51_submat_id(1:4) = ipm(51:54, matid)
207 jj = minloc(m51_submat_id,1)
208 nbmat_target = 4
209 IF(m51_submat_id(jj)==0)nbmat_target=max(1,jj-1)
210 ELSEIF (mlw == 151)THEN
211 nbmat_target = multi_fvm%NBMAT
212 ELSE
213 nbmat_target = 1
214 ENDIF
215c
216 nuvar = ipm(8, matid)
217 iadbuf = ipm(7, matid)
218 isolnod = iparg(28, ng)
219 gbuf => elbuf_tab(ng)%GBUF
220C Find first index of elem_list such as corresponding elem lies in the group
221 ifirst = 1
222 shift = 0
223 DO WHILE(elem_list(ifirst + shift) < 1 + nft .AND. ifirst + shift <= nelem)
224 shift = shift + 1
225 ENDDO
226 IF (shift > nelem) THEN
227 cycle
228 ELSE
229 ifirst = ifirst + shift
230 ENDIF
231
232 ilast = ifirst
233 shift = 0
234 DO WHILE(elem_list(ilast + shift) <= nel + nft .AND. ilast + shift < nelem)
235 shift = shift + 1
236 IF( shift >= nel)EXIT
237 ENDDO
238C IFIRST + SHIFT is the first element which ID is greater than NFT + NEL
239 IF (ilast + shift < nelem) THEN
240 shift = shift - 1
241 ilast = ilast + shift
242 ELSE
243 ilast = nelem
244 ENDIF
245
246 !IF some elems are targeted in this current group NG, check material law
247 IF(ifirst>0 .AND. ifirst<=ilast )THEN
248 SELECT CASE(mlw)
249 CASE (3, 4, 6, 49, 51, 151)
250 CONTINUE
251 CASE DEFAULT
252 !WRITE(IOUT, '(A, I10)') "**WARNING : /INIMAP1D OPTION NOT COMPATIBLE WITH MATERIAL LAW", MLW
253 cycle
254 END SELECT
255 IF (mlw == 51) THEN
256 iform = ipm(62, matid)
257 IF (iform /= 12 .AND. inimap1d(kk)%FORMULATION == 1) THEN
258 CALL ancmsg(msgid = 1732, msgtype=msgerror, anmode=aninfo,i1 = ipm(1, matid), i2 = inimap1d(kk)%ID)
259 ENDIF
260 ENDIF
261 ENDIF
262
263C Check submaterial consistency
264 IF(tag_mat(matid)==0)THEN
265 IF(nbmat /= nbmat_target)THEN
266 WRITE(str1,'(i2)')nbmat
267 WRITE(str2,'(i2)')nbmat_target
268 CALL ancmsg(msgid = 2048, msgtype=msgerror, anmode=aninfo,
269 . c1="INIMAP1D",
270 . c2="NUMBER OF SUBMATERIAL(S) FROM TARGET MATERIAL LAW:"//str1,
271 . c3=" IS DIFFERENT FROM SUBMATERIAL(S) NUMBER IN TARGET DOMAIN:"//str2,
272 . c4="INIMAP1D",
273 . i1 = ipm(1, matid),
274 . i2 = inimap1d(kk)%ID)
275 tag_mat(matid) = 1 !display this error material once per material law. Not for each cell using this material law.
276 nbmat=min(nbmat,nbmat_target) !in order to finish starter without memory corruption
277 ENDIF
278 ENDIF
279
280 DO ii = ifirst, ilast
281 elemid = elem_list(ii)
282 IF (elemid < 1 + nft .OR. elemid > nel + nft) THEN
283 ENDIF
284 IF (ity == 1 .AND. isolnod == 8) THEN
285C HEXA
286 node1 = ixs(2, elemid)
287 node2 = ixs(3, elemid)
288 node3 = ixs(4, elemid)
289 node4 = ixs(5, elemid)
290 node5 = ixs(6, elemid)
291 node6 = ixs(7, elemid)
292 node7 = ixs(8, elemid)
293 node8 = ixs(9, elemid)
294 xc = one_over_8 * (xgrid(1, node1) + xgrid(1, node2) + xgrid(1, node3) + xgrid(1, node4)
295 . + xgrid(1, node5) + xgrid(1, node6) + xgrid(1, node7) + xgrid(1, node8))
296 yc = one_over_8 * (xgrid(2, node1) + xgrid(2, node2) + xgrid(2, node3) + xgrid(2, node4)
297 . + xgrid(2, node5) + xgrid(2, node6) + xgrid(2, node7) + xgrid(2, node8))
298 zc = one_over_8 * (xgrid(3, node1) + xgrid(3, node2) + xgrid(3, node3) + xgrid(3, node4)
299 . + xgrid(3, node5) + xgrid(3, node6) + xgrid(3, node7) + xgrid(3, node8))
300 nbnode = 8
301 nodelist(1) = 2
302 nodelist(2) = 3
303 nodelist(3) = 4
304 nodelist(4) = 5
305 nodelist(5) = 6
306 nodelist(6) = 7
307 nodelist(7) = 8
308 nodelist(8) = 9
309 ELSEIF (ity == 1 .AND. isolnod == 4) THEN
310C TETRA
311 node1 = ixs(2, elemid)
312 node2 = ixs(4, elemid)
313 node3 = ixs(7, elemid)
314 node4 = ixs(6, elemid)
315 xc = fourth * (xgrid(1, node1) + xgrid(1, node2) + xgrid(1, node3) + xgrid(1, node4))
316 yc = fourth * (xgrid(2, node1) + xgrid(2, node2) + xgrid(2, node3) + xgrid(2, node4))
317 zc = fourth * (xgrid(3, node1) + xgrid(3, node2) + xgrid(3, node3) + xgrid(3, node4))
318 nbnode = 4
319 nodelist(1) = 2
320 nodelist(2) = 4
321 nodelist(3) = 7
322 nodelist(4) = 6
323 ELSEIF (ity == 2) THEN
324C QUAD
325 node1 = ixq(2, elemid)
326 node2 = ixq(3, elemid)
327 node3 = ixq(4, elemid)
328 node4 = ixq(5, elemid)
329 xc = fourth * (xgrid(1, node1) + xgrid(1, node2) + xgrid(1, node3) + xgrid(1, node4))
330 yc = fourth * (xgrid(2, node1) + xgrid(2, node2) + xgrid(2, node3) + xgrid(2, node4))
331 zc = fourth * (xgrid(3, node1) + xgrid(3, node2) + xgrid(3, node3) + xgrid(3, node4))
332 nbnode = 4
333 nodelist(1) = 2
334 nodelist(2) = 3
335 nodelist(3) = 4
336 nodelist(4) = 5
337 ELSEIF (ity == 7) THEN
338C TRIANGLES
339 node1 = ixtg(2, elemid)
340 node2 = ixtg(3, elemid)
341 node3 = ixtg(4, elemid)
342 xc = third * (xgrid(1, node1) + xgrid(1, node2) + xgrid(1, node3))
343 yc = third * (xgrid(2, node1) + xgrid(2, node2) + xgrid(2, node3))
344 zc = third * (xgrid(3, node1) + xgrid(3, node2) + xgrid(3, node3))
345 nbnode = 3
346 nodelist(1) = 2
347 nodelist(2) = 3
348 nodelist(3) = 4
349 ENDIF
350 IF (inimap1d(kk)%PROJ == 3) THEN
351C SPHERICAL
352 radius = sqrt((xc - x0) * (xc - x0) + (yc - y0) * (yc - y0) + (zc - z0) * (zc - z0))
353 ELSEIF (inimap1d(kk)%PROJ == 1) THEN
354C PLANAR
355 radius = (xc - x0) * inimap1d(kk)%NX +
356 . (yc - y0) * inimap1d(kk)%NY +
357 . (zc - z0) * inimap1d(kk)%NZ
358 ELSE
359C CYLINDRICAL PROJ
360 tt = (xc - x0) * (x1 - x0) + (yc - y0) * (y1 - y0) + (zc - z0) * (z1 - z0)
361 tt = tt / ((x1 - x0) * (x1 - x0) +
362 . (y1 - y0) * (y1 - y0) +
363 . (z1 - z0) * (z1 - z0))
364 xp = x0 + tt * (x1 - x0)
365 yp = y0 + tt * (y1 - y0)
366 zp = z0 + tt * (z1 - z0)
367 radius = sqrt((xc - xp) * (xc - xp) + (yc - yp) * (yc - yp) + (zc - zp) * (zc - zp))
368 ENDIF
369
370 DO imat = 1, nbmat
371C Alpha
372 ifunc = inimap1d(kk)%FUNC_ALPHA(imat)
373 IF (ifunc /= 0) THEN
374 IF(ifunc>0)THEN
375 !--function_identifier
376 npt = (npc(ifunc + 1) - npc(ifunc)) / 2
377 first = 1 + (npc(ifunc) - 1) / 2
378 ptrx => pld_inv(1:npts/2,1) !NPTS : 2*nb of all function points NPT:for current function
379 ptry => pld_inv(1:npts/2,2)
380 ELSE
381 !--data from file
382 npt = inimap1d(kk)%NUM_CENTROIDS
383 first = 1
384 ptrx => inimap1d(kk)%X(1:npt)
385 ptry => inimap1d(kk)%SUBMAT(imat)%VFRAC(1:npt)
386 ENDIF
387 shift = 0
388 DO WHILE (ptrx(first + shift) < radius )
389 shift = shift + 1
390 IF(shift >= npt)EXIT
391 ENDDO
392 ind = first + shift
393 IF (shift == 0) THEN
394 res = ptry(first)
395 IF (mlw == 51) THEN
396 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
397 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (1 + kphase - 1)) = res
398 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (23 + kphase - 1)) = res
399 ELSE
400 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%VOL(elemid - nft) = res * gbuf%VOL(elemid - nft)
401 ENDIF
402 ELSEIF (shift < npt) THEN
403 rad1 = ptrx(ind - 1)
404 rad2 = ptrx(ind)
405 value1 = ptry(ind - 1)
406 value2 = ptry(ind)
407 res = (value1 + (value2 - value1) / (rad2 - rad1) * (radius - rad1))
408 IF (mlw == 51) THEN
409 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
410 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (1 + kphase - 1)) = res
411 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (23 + kphase - 1)) = res
412 ELSE
413 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%VOL(elemid - nft) = res * gbuf%VOL(elemid - nft)
414 ENDIF
415 ELSE
416 res = ptry(first + npt - 1)
417 IF (mlw == 51) THEN
418 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
419 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (1 + kphase - 1)) = res
420 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (23 + kphase - 1)) = res
421 ELSE
422 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%VOL(elemid - nft) = res * gbuf%VOL(elemid - nft)
423 ENDIF
424 ENDIF
425 ELSE
426 res = one
427 IF (mlw == 51) THEN
428 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
429 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (1 + kphase - 1)) = res
430 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (23 + kphase - 1)) = res
431 ELSE
432 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%VOL(elemid - nft) = res * gbuf%VOL(elemid - nft)
433 ENDIF
434 ENDIF
435 vfrac(elemid-nft,imat) = res
436C RHO
437 ifunc = inimap1d(kk)%FUNC_RHO(imat)
438 IF(ifunc /= 0)THEN
439 IF(ifunc > 0)THEN
440 !--function_identifier
441 npt = (npc(ifunc + 1) - npc(ifunc)) / 2
442 first = 1 + (npc(ifunc) - 1) / 2
443 ptrx => pld_inv(1:npts/2,1) !NPTS : 2*nb of all function points NPT:for current function
444 ptry => pld_inv(1:npts/2,2)
445 ELSE
446 !--data from file
447 npt = inimap1d(kk)%NUM_CENTROIDS
448 first = 1
449 ptrx => inimap1d(kk)%X(1:npt)
450 ptry => inimap1d(kk)%SUBMAT(imat)%RHO(1:npt)
451 ENDIF
452 shift = 0
453 DO WHILE (ptrx(first + shift) < radius )
454 shift = shift + 1
455 IF( shift >= npt)EXIT
456 ENDDO
457 ind = first + shift
458 IF (shift == 0) THEN
459 res = ptry(first) * inimap1d(kk)%FAC_RHO(imat)
460 IF (mlw == 51) THEN
461 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
462 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (9 + kphase - 1)) = res
463 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (12 + kphase - 1)) = res
464 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (20 + kphase - 1)) = res
465 ELSE
466 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%RHO(elemid - nft) = res
467 ENDIF
468 ELSEIF (shift < npt) THEN
469 rad1 = ptrx(ind - 1)
470 rad2 = ptrx(ind)
471 value1 = ptry(ind - 1)
472 value2 = ptry(ind)
473 res = (value1 + (value2 - value1) / (rad2 - rad1) * (radius - rad1)) * inimap1d(kk)%FAC_RHO(imat)
474 IF (mlw == 51) THEN
475 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
476 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (9 + kphase - 1)) = res
477 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (12 + kphase - 1)) = res
478 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (20 + kphase - 1)) = res
479 ELSE
480 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%RHO(elemid - nft) = res
481 ENDIF
482 ELSE
483 res = ptry(first + npt - 1) * inimap1d(kk)%FAC_RHO(imat)
484 IF (mlw == 51) THEN
485 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
486 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (9 + kphase - 1)) = res
487 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (12 + kphase - 1)) = res
488 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (20 + kphase - 1)) = res
489 ELSE
490 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%RHO(elemid - nft) = res
491 ENDIF
492 ENDIF
493 ELSE
494 res = inimap1d(kk)%FAC_RHO(imat)
495 IF (mlw == 51) THEN
496 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
497 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (9 + kphase - 1)) = res
498 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (12 + kphase - 1)) = res
499 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (20 + kphase - 1)) = res
500 ELSE
501 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%RHO(elemid - nft) = res
502 ENDIF
503 ENDIF
504
505 IF (inimap1d(kk)%FORMULATION == 2) THEN
506C EINT
507 !not managed when reading from file (FORMULATION=1)
508 ifunc = inimap1d(kk)%FUNC_ENER(imat)
509 IF (ifunc > 0) THEN
510 npt = (npc(ifunc + 1) - npc(ifunc)) / 2
511 first = 1 + (npc(ifunc) - 1) / 2
512 shift = 0
513 DO WHILE (ptrx(first + shift) < radius)
514 shift = shift + 1
515 IF( shift >= npt)EXIT
516 ENDDO
517 ind = first + shift
518 IF (shift == 0) THEN
519 res = ptry(first) * inimap1d(kk)%FAC_PRES_ENER(imat)
520 IF (mlw == 51) THEN
521 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
522 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (8 + kphase - 1)) = res
523 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (21 + kphase - 1)) = res
524 ELSE
525 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%EINT(elemid - nft) = res
526 ENDIF
527 ELSEIF (shift < npt) THEN
528 rad1 = ptrx(ind - 1)
529 rad2 = ptrx(ind)
530 value1 = ptry(ind - 1)
531 value2 = ptry(ind)
532 res = (value1 + (value2 - value1) / (rad2 - rad1) * (radius - rad1)) * inimap1d(kk)%FAC_PRES_ENER(imat)
533 IF (mlw == 51) THEN
534 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
535 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (8 + kphase - 1)) = res
536 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (21 + kphase - 1)) = res
537 ELSE
538 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%EINT(elemid - nft) = res
539 ENDIF
540 ELSE
541 res = ptry(first + npt - 1) * inimap1d(kk)%FAC_PRES_ENER(imat)
542 IF (mlw == 51) THEN
543 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
544 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (8 + kphase - 1)) = res
545 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (21 + kphase - 1)) = res
546 ELSE
547 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%EINT(elemid - nft) = res
548 ENDIF
549 ENDIF
550 ELSE
551 res = inimap1d(kk)%FAC_PRES_ENER(imat)
552 IF (mlw == 51) THEN
553 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
554 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (8 + kphase - 1)) = res
555 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (21 + kphase - 1)) = res
556 ELSE
557 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%EINT(elemid - nft) = res
558 ENDIF
559 ENDIF
560 eint(elemid-nft,imat) = res
561 ELSE IF (inimap1d(kk)%FORMULATION == 1) THEN
562C PRES
563 ifunc = inimap1d(kk)%FUNC_PRES(imat)
564 IF(ifunc /= 0)THEN
565 IF(ifunc > 0)THEN
566 !--function_identifier
567 npt = (npc(ifunc + 1) - npc(ifunc)) / 2
568 first = 1 + (npc(ifunc) - 1) / 2
569 ptrx => pld_inv(1:npts/2,1) !NPTS : 2*nb of all function points NPT:for current function
570 ptry => pld_inv(1:npts/2,2)
571 ELSE
572 !--data from file
573 npt = inimap1d(kk)%NUM_CENTROIDS
574 first = 1
575 ptrx => inimap1d(kk)%X(1:npt)
576 ptry => inimap1d(kk)%SUBMAT(imat)%PRES(1:npt)
577 ENDIF
578 shift = 0
579 DO WHILE (ptrx(first + shift) < radius)
580 shift = shift + 1
581 IF( shift >= npt)EXIT
582 ENDDO
583 ind = first + shift
584 IF (shift == 0) THEN
585 pres(elemid - nft, imat) = ptry(first) * inimap1d(kk)%FAC_PRES_ENER(imat)
586 ELSEIF (shift < npt) THEN
587 rad1 = ptrx(ind - 1)
588 rad2 = ptrx(ind)
589 value1 = ptry(ind - 1)
590 value2 = ptry(ind)
591 pres(elemid - nft, imat) = (value1 + (value2 - value1) / (rad2 - rad1) * (radius - rad1)) *
592 . * inimap1d(kk)%FAC_PRES_ENER(imat)
593 ELSE
594 value1 = ptry(first + npt - 1)
595 pres(elemid - nft, imat) = value1 * inimap1d(kk)%FAC_PRES_ENER(imat)
596 ENDIF
597 ELSE
598 pres(elemid - nft, imat) = inimap1d(kk)%FAC_PRES_ENER(imat)
599 ENDIF
600
601
602 IF(multi_fvm%IS_USED)THEN
603 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%SIG(nel*(1-1)+elemid-nft)=-pres(elemid - nft, imat)
604 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%SIG(nel*(2-1)+elemid-nft)=-pres(elemid - nft, imat)
605 elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%SIG(nel*(3-1)+elemid-nft)=-pres(elemid - nft, imat)
606 elbuf_tab(ng)%GBUF%SIG(nel*(1-1)+elemid-nft)=-pres(elemid - nft, imat) !same pressure for all phases
607 elbuf_tab(ng)%GBUF%SIG(nel*(2-1)+elemid-nft)=-pres(elemid - nft, imat)
608 elbuf_tab(ng)%GBUF%SIG(nel*(3-1)+elemid-nft)=-pres(elemid - nft, imat)
609 ELSE
610 elbuf_tab(ng)%GBUF%SIG(nel*(1-1)+elemid-nft)=-pres(elemid - nft, 1)
611 elbuf_tab(ng)%GBUF%SIG(nel*(2-1)+elemid-nft)=-pres(elemid - nft, 1)
612 elbuf_tab(ng)%GBUF%SIG(nel*(3-1)+elemid-nft)=-pres(elemid - nft, 1)
613 IF(mlw==51)THEN
614 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
615 !ELBUF_TAB(NG)%BUFLY(1)%MAT(1,1,1)%VAR(ELEMID - NFT + NEL * (2 + KPHASE - 1)) = -PRES(ELEMID - NFT, 1)
616 !ELBUF_TAB(NG)%BUFLY(1)%MAT(1,1,1)%VAR(ELEMID - NFT + NEL * (3 + KPHASE - 1)) = -PRES(ELEMID - NFT, 1)
617 !ELBUF_TAB(NG)%BUFLY(1)%MAT(1,1,1)%VAR(ELEMID - NFT + NEL * (4 + KPHASE - 1)) = -PRES(ELEMID - NFT, 1)
618 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (18+ kphase - 1)) = pres(elemid - nft, 1)
619 ENDIF
620 ENDIF
621
622 ENDIF
623
624 ENDDO
625
626C Burn Fraction / Detonation Times
627 IF(elbuf_tab(ng)%GBUF%G_TB > 0)elbuf_tab(ng)%GBUF%TB(elemid - nft) = 1e20
628 IF (mlw == 51) THEN
629 kphase = m51_n0phas - 1
630 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (kphase - 1)) = one
631 IF(elbuf_tab(ng)%GBUF%G_BFRAC > 0)elbuf_tab(ng)%GBUF%BFRAC(elemid - nft) = one
632 kphase = m51_n0phas + (4 - 1) * m51_nvphas ! JWL is submat 4
633 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(1 + nel * (15 + kphase - 1)) = 1e20
634 ELSE
635 DO imat=1,nbmat
636 IF(elbuf_tab(ng)%BUFLY(imat)%L_BFRAC > 0)elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%BFRAC(elemid - nft)=one
637 ENDDO
638 IF(elbuf_tab(ng)%GBUF%G_BFRAC > 0)elbuf_tab(ng)%GBUF%BFRAC(elemid - nft) = one
639 ENDIF
640
641
642 IF (.NOT. multi_fvm%IS_USED .AND. alefvm_param%IEnabled == 0) THEN
643C VEL
644 DO inode = 1, nbnode
645 IF (ity == 1) THEN
646 nodeid = ixs(nodelist(inode), elemid)
647 ELSEIF (ity == 2) THEN
648 nodeid = ixq(nodelist(inode), elemid)
649 ELSEIF (ity == 7) THEN
650 nodeid = ixtg(nodelist(inode), elemid)
651 ENDIF
652 IF (inimap1d(kk)%TAGNODE(nodeid) == 0) THEN
653
654 xnode = xgrid(1, nodeid)
655 ynode = xgrid(2, nodeid)
656 znode = xgrid(3, nodeid)
657 IF (inimap1d(kk)%PROJ == 3) THEN
658C SPHERICAL PROJ
659 radius = sqrt((xnode - x0) * (xnode - x0) +
660 . (ynode - y0) * (ynode - y0) +
661 . (znode - z0) * (znode - z0))
662 ELSE IF (inimap1d(kk)%PROJ == 1) THEN
663C PLANAR PROJ
664 radius = (xnode - x0) * inimap1d(kk)%NX +
665 . (ynode - y0) * inimap1d(kk)%NY +
666 . (znode - z0) * inimap1d(kk)%NZ
667 ELSE
668C CYLINDRICAL PROJ
669 tt = (xnode - x0) * (x1 - x0) +
670 . (ynode - y0) * (y1 - y0) +
671 . (znode - z0) * (z1 - z0)
672 tt = tt / ((x1 - x0) * (x1 - x0) +
673 . (y1 - y0) * (y1 - y0) +
674 . (z1 - z0) * (z1 - z0))
675 xp = x0 + tt * (x1 - x0)
676 yp = y0 + tt * (y1 - y0)
677 zp = z0 + tt * (z1 - z0)
678 radius = sqrt((xnode - xp) * (xnode - xp) +
679 . (ynode - yp) * (ynode - yp) +
680 . (znode - zp) * (znode - zp))
681 ENDIF
682 IF (abs(radius) < em10
683 . .AND. inimap1d(kk)%PROJ /= 1) THEN
684 vel(1:3, nodeid) = zero
685 ELSE
686 IF (inimap1d(kk)%PROJ == 3) THEN
687C SPHERICAL
688 nx = (xnode - x0) / radius
689 ny = (ynode - y0) / radius
690 nz = (znode - z0) / radius
691 ELSE IF (inimap1d(kk)%PROJ == 1) THEN
692C PLANAR
693 nx = inimap1d(kk)%NX
694 ny = inimap1d(kk)%NY
695 nz = inimap1d(kk)%NZ
696 ELSE
697C CYLINDRICAL
698 nx = (xnode - xp) / radius
699 ny = (ynode - yp) / radius
700 nz = (znode - zp) / radius
701 ENDIF
702 ifunc = inimap1d(kk)%FUNC_VEL
703 IF(ifunc /= 0)THEN
704 IF(ifunc > 0)THEN
705 !--function_identifier
706 npt = (npc(ifunc + 1) - npc(ifunc)) / 2
707 first = 1 + (npc(ifunc) - 1) / 2
708 ptrx => pld_inv(1:npts/2,1) !NPTS : 2*nb of all function points NPT:for current function
709 ptry => pld_inv(1:npts/2,2)
710 ELSE
711 !--data from file
712 npt = inimap1d(kk)%NUM_NODE_VEL
713 first = 1
714 ptrx => inimap1d(kk)%X_VEL(1:npt)
715 ptry => inimap1d(kk)%VEL(1:npt)
716 ENDIF
717 shift = 0
718 DO WHILE (ptrx(first + shift) < radius )
719 shift = shift + 1
720 IF( shift >= npt)EXIT
721 ENDDO
722 ind = first + shift
723 IF (shift == 0) THEN
724 VALUE = ptry(first)
725 ELSEIF (shift < npt) THEN
726 rad1 = ptrx(ind - 1)
727 rad2 = ptrx(ind)
728 value1 = ptry(ind - 1)
729 value2 = ptry(ind)
730 VALUE = value1 + (value2 - value1) / (rad2 - rad1) * (radius - rad1)
731 ELSE
732 value1 = ptry(first + npt - 1)
733 VALUE = value1
734 ENDIF
735 VALUE = VALUE * inimap1d(kk)%FAC_VEL
736 ELSE
737 VALUE = zero
738 ENDIF
739 vel(1, nodeid) = VALUE * nx
740 vel(2, nodeid) = VALUE * ny
741 vel(3, nodeid) = VALUE * nz
742 ENDIF
743 inimap1d(kk)%TAGNODE(nodeid) = 1
744 ENDIF
745 ENDDO ! INODE
746 ELSE
747
748C Initialize cell centred velocities
749 ifunc = inimap1d(kk)%FUNC_VEL
750 IF(ifunc /= 0)THEN
751 IF(ifunc > 0)THEN
752 !--function_identifier
753 npt = (npc(ifunc + 1) - npc(ifunc)) / 2
754 first = 1 + (npc(ifunc) - 1) / 2
755 ptrx => pld_inv(1:npts/2,1) !NPTS : 2*nb of all function points NPT:for current function
756 ptry => pld_inv(1:npts/2,2)
757 ELSE
758 !--data from file
759 npt = inimap1d(kk)%NUM_NODE_VEL
760 first = 1
761 ptrx => inimap1d(kk)%X_VEL(1:npt)
762 ptry => inimap1d(kk)%VEL(1:npt)
763 ENDIF
764 shift = 0
765 DO WHILE (ptrx(first + shift) < radius )
766 shift = shift + 1
767 IF( shift >= npt)EXIT
768 ENDDO
769 ind = first + shift
770
771 IF (shift == 0) THEN
772 VALUE = ptry(first)
773 ELSEIF (shift < npt) THEN
774 rad1 = ptrx(ind - 1)
775 rad2 = ptrx(ind)
776 value1 = ptry(ind - 1)
777 value2 = ptry(ind)
778 VALUE = value1 + (value2 - value1) / (rad2 - rad1) * (radius - rad1)
779 ELSE
780 value1 = ptry(first + npt - 1)
781 VALUE = value1
782 ENDIF
783 nx = zero
784 ny = zero
785 nz = zero
786 IF (inimap1d(kk)%PROJ == 3) THEN
787C SPHERICAL
788 IF (abs(radius) > em10) THEN
789 nx = (xc - x0) / radius
790 ny = (yc - y0) / radius
791 nz = (zc - z0) / radius
792 ENDIF
793 ELSE IF (inimap1d(kk)%PROJ == 1) THEN
794C PLANAR
795 nx = inimap1d(kk)%NX
796 ny = inimap1d(kk)%NY
797 nz = inimap1d(kk)%NZ
798 ELSE
799C CYLINDRICAL
800 IF (abs(radius) > em10) THEN
801 nx = (xc - xp) / radius
802 ny = (yc - yp) / radius
803 nz = (zc - zp) / radius
804 ENDIF
805 ENDIF
806 VALUE = VALUE * inimap1d(kk)%FAC_VEL
807 ELSE
808 VALUE = zero
809 ENDIF
810 IF (multi_fvm%IS_USED) THEN
811 multi_fvm%VEL(1, elemid) = VALUE * nx
812 multi_fvm%VEL(2, elemid) = VALUE * ny
813 multi_fvm%VEL(3, elemid) = VALUE * nz
814 ELSE IF (alefvm_param%IEnabled > 0) THEN
815 gbuf%MOM(elemid - nft + 0 * nel) = VALUE * nx * gbuf%RHO(elemid - nft)
816 gbuf%MOM(elemid - nft + 1 * nel) = VALUE * ny * gbuf%RHO(elemid - nft)
817 gbuf%MOM(elemid - nft + 2 * nel) = VALUE * nz * gbuf%RHO(elemid - nft)
818 ENDIF
819 ENDIF ! MULTI_FVM
820 ENDDO ! II = IFIRST, ILAST
821
822 IF (inimap1d(kk)%FORMULATION == 2) THEN
823 DO ii = ifirst, ilast
824 elemid = elem_list(ii)
825 gbuf%EINT(elemid - nft) = gbuf%EINT(elemid - nft) * gbuf%RHO(elemid - nft)
826 gbuf%TEMP(elemid - nft) = 300*(gbuf%RHO(elemid - nft) )**0.4
827 ENDDO
828
829 IF (mlw /= 51 .AND. mlw /= 151) THEN
830 submat(1)%TEMP(1:nel) => elbuf_tab(ng)%GBUF%TEMP(1:nel)
831 ELSE IF (mlw == 151) THEN
832 DO imat=1,nbmat
833 submat(imat)%TEMP(1:nel) => elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%TEMP(1:nel)
834 ENDDO
835 ELSE
836 DO imat=1,nbmat
837 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
838 submat(imat)%TEMP(1:nel) => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(1+nel*(16+kphase-1) : nel*(16+kphase ))
839 ENDDO
840 ENDIF
841
842 ELSE IF (inimap1d(kk)%FORMULATION == 1) THEN
843C Compute EINT from pressure
844 IF (mlw /= 151 .AND. mlw /= 51) THEN
845 IF (ity == 1) THEN
846 matid = ixs(1, 1 + nft)
847 ELSEIF (ity == 2) THEN
848 matid = ixq(1, 1 + nft)
849 ELSEIF (ity == 7) THEN
850 matid = ixtg(1, 1 + nft)
851 ENDIF
852 nvartmp_eos = elbuf_tab(ng)%BUFLY(1)%NVARTMP_EOS
853 submat(1)%TEMP(1:nel) => elbuf_tab(ng)%GBUF%TEMP(1:nel)
854 ebuf => elbuf_tab(ng)%BUFLY(1)%EOS(1,1,1)
855 nvareos = elbuf_tab(ng)%BUFLY(1)%NVAR_EOS
856 CALL multi_solve_eint(matid, nft, nel, pres(:, 1), gbuf%EINT, gbuf%RHO, ifirst, ilast, elem_list,
857 . ipm, pm, bufmat, mlw, submat(1)%TEMP,snpc,stf,npc,pld, ebuf%VAR, nvareos,
858 . mat_param(matid),nvartmp_eos,ebuf%VARTMP,nummat,npropmi,npropm)
859 ELSE IF (mlw == 151) THEN
860 DO imat = 1, nbmat
861 IF (ity == 1) THEN
862 gmid = ixs(1, 1 + nft)
863 matid = mat_param(gmid)%MULTIMAT%MID(imat)
864 ELSEIF (ity == 2) THEN
865 gmid = ixq(1, 1 + nft)
866 matid = mat_param(gmid)%MULTIMAT%MID(imat)
867 ELSEIF (ity == 7) THEN
868 gmid = ixtg(1, 1 + nft)
869 matid = mat_param(gmid)%MULTIMAT%MID(imat)
870 ENDIF
871 nvartmp_eos = elbuf_tab(ng)%BUFLY(imat)%NVARTMP_EOS
872 submat(imat)%TEMP(1:nel) => elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%TEMP(1:nel)
873 ebuf => elbuf_tab(ng)%BUFLY(imat)%EOS(1,1,1)
874 nvareos = elbuf_tab(ng)%BUFLY(imat)%NVAR_EOS
875 CALL multi_solve_eint(matid, nft, nel, pres(:, imat), elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%EINT,
876 . elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%RHO, ifirst, ilast, elem_list,
877 . ipm, pm, bufmat, mlw,submat(imat)%TEMP,snpc,stf,npc,pld,ebuf%VAR,nvareos,
878 . mat_param(matid),nvartmp_eos,ebuf%VARTMP,nummat,npropmi,npropm)
879 ENDDO
880 ELSE
881 DO imat = 1, nbmat
882 IF (ity == 1) THEN
883 gmid = ixs(1, 1 + nft)
884 matid = mat_param(gmid)%MULTIMAT%MID(imat)
885 ELSEIF (ity == 2) THEN
886 gmid = ixq(1, 1 + nft)
887 matid = mat_param(gmid)%MULTIMAT%MID(imat)
888 ELSEIF (ity == 7) THEN
889 gmid = ixtg(1, 1 + nft)
890 matid = mat_param(gmid)%MULTIMAT%MID(imat)
891 ENDIF
892 nvartmp_eos = elbuf_tab(ng)%BUFLY(1)%NVARTMP_EOS
893 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
894 submat(imat)%TEMP(1:nel) => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(1+nel*(16+kphase-1) : nel*(16+kphase ))
895 ebuf => elbuf_tab(ng)%BUFLY(1)%EOS(1,1,1)
896 nvareos = elbuf_tab(ng)%BUFLY(1)%NVAR_EOS
897 CALL multi_solve_eint(matid, nft, nel, pres(:, imat), eint(:, imat),
898 . elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(1 + nel * (9 + kphase - 1):nel + nel * (9 + kphase - 1)),
899 . ifirst, ilast, elem_list,
900 . ipm, pm, bufmat, mlw, submat(imat)%TEMP,snpc,stf,npc,pld,ebuf%VAR,nvareos,
901 . mat_param(matid),nvartmp_eos,ebuf%VARTMP,nummat,npropmi,npropm)
902 DO ii = ifirst, ilast
903 elemid = elem_list(ii)
904 res = eint(elemid - nft, imat)
905 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
906 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (8 + kphase - 1)) = res
907 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(elemid - nft + nel * (21 + kphase - 1)) = res
908 ENDDO
909 ENDDO
910 ENDIF
911
912 IF(mlw == 151)THEN
913 IF(nbmat > 1)THEN
914 DO ii = 1,nel
915 elbuf_tab(ng)%GBUF%TEMP(ii) = zero
916 DO imat=1,nbmat
917 densfrac = elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)%RHO(ii) / elbuf_tab(ng)%GBUF%RHO(ii)
918 elbuf_tab(ng)%GBUF%TEMP(ii) = elbuf_tab(ng)%GBUF%TEMP(ii) + vfrac(ii,imat)*densfrac * submat(imat)%TEMP(ii)
919 ENDDO
920 ENDDO
921 ENDIF
922 ELSEIF(mlw==51)THEN
923 DO ii = 1,nel
924 elbuf_tab(ng)%GBUF%TEMP(ii) = zero
925 DO imat=1,nbmat
926 kphase = m51_n0phas + (m51_submat_id(imat) - 1) * m51_nvphas
927 densfrac = elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(ii + nel * (12 + kphase - 1)) / elbuf_tab(ng)%GBUF%RHO(ii)
928 elbuf_tab(ng)%GBUF%TEMP(ii) = elbuf_tab(ng)%GBUF%TEMP(ii) + vfrac(ii,imat)*densfrac * submat(imat)%TEMP(ii)
929 ENDDO
930 elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR(1:nel)= elbuf_tab(ng)%GBUF%TEMP(1:nel) !global temperature law51 buffer
931 ENDDO
932 ENDIF
933
934 ENDIF
935 ENDDO ! NG = 1, NGROUP
936 DEALLOCATE(pres, eint)
937 ENDDO
938C ===================
939C Memory deallocation
940C ===================
941 DEALLOCATE(elem_list)
942
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
alefvm_mod
Definition alefvm_mod.F:88
alefvm_mod::alefvm_param
type(alefvm_param_), target alefvm_param
Definition alefvm_mod.F:121
groupdef_mod
Definition groupdef_mod.F:662
inimap1d_mod
Definition inimap1d_mod.F:35
message_mod
Definition message_mod.F:1249
quicksort_i
recursive subroutine quicksort_i(a, first, last)
Definition quicksort.F:92
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889