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