40
41
42
43 USE elbufdef_mod
44 USE multi_fvm_mod
48 USE matparam_def_mod, ONLY : matparam_struct_
49 use element_mod , only : nixtg
50
51
52
53#include "implicit_f.inc"
54
55
56
57
58#include "mvsiz_p.inc"
59
60
61
62
63#include "param_c.inc"
64
65#include "vect01_c.inc"
66
67#include "scry_c.inc"
68
69#include "scr17_c.inc"
70
71#include "com04_c.inc"
72
73#include "com01_c.inc"
74
75
76
77 TYPE(ELBUF_STRUCT_), INTENT(IN), TARGET :: ELBUF_STR
78 INTEGER, INTENT(IN) :: NEL, NSIGS, IXTG(NIXTG, *),
79 . IGEO(*), PTSH3N(*), NPF(*), ILOADP(SIZLOADP, *),
80 . IPART(LIPART1, *), IPARTTG(*), IPM(NPROPMI, *)
81 INTEGER, INTENT(INOUT) :: IPARG(*)
82 INTEGER, INTENT(OUT) :: NVC
83 my_real,
INTENT(IN) :: xgrid(3, *), facload(lfacload, *)
84 my_real,
INTENT(INOUT) :: pm(npropm, *)
85 my_real,
INTENT(INOUT) :: geo(*), sigi(nsigs, *), skew(lskew, *),
86 . tf(*), bufmat(*)
87 TYPE(MULTI_FVM_STRUCT) :: MULTI_FVM
88 LOGICAL :: ERROR_THROWN
89 TYPE(DETONATORS_STRUCT_) DETONATORS
90 TYPE(t_ale_connectivity), INTENT(INOUT) :: ALE_CONNECTIVITY
91 TYPE(MATPARAM_STRUCT_) ,DIMENSION(NUMMAT) ,INTENT(IN) :: MAT_PARAM
92
93
94
95 TYPE(L_BUFEL_) ,POINTER :: LBUF
96 TYPE(G_BUFEL_) ,POINTER :: GBUF
97 TYPE(BUF_MAT_) ,POINTER :: MBUF
98 INTEGER :: NLAY, NGL(MVSIZ), MAT(MVSIZ), (MVSIZ), IX2(MVSIZ), IX3(MVSIZ)
100 . y1(mvsiz), z1(mvsiz),
101 . y2(mvsiz), z2(mvsiz),
102 . y3(mvsiz), z3(mvsiz),
103 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
104 . lgth1(mvsiz), lgth2(mvsiz), lgth3(mvsiz),pres,vfrac
106 INTEGER :: II, I, IP, IBID, ILAY, MATLAW
107
108
109
110 IF(n2d>0 .AND. mtn /=151) THEN
112 . msgtype=msgerror,
113 . anmode=aninfo_blind_1)
114 RETURN
115 ENDIF
116
117
118
119 tempel(:) = zero
120
121 gbuf => elbuf_str%GBUF
122
123 nlay = elbuf_str%NLAY
124
125 DO ii = 1, nel
126 i = ii + nft
127 ix1(ii) = ixtg(1 + 1, i)
128 ix2(ii) = ixtg(1 + 2, i)
129 ix3(ii) = ixtg(1 + 3, i)
130 y1(ii) = xgrid(2, ixtg(1 + 1, i))
131 z1(ii) = xgrid(3, ixtg(1 + 1, i))
132 y2(ii) = xgrid(2, ixtg(1 + 2, i))
133 z2(ii) = xgrid(3, ixtg(1 + 2, i))
134 y3(ii) = xgrid(2, ixtg(1 + 3, i))
135 z3(ii) = xgrid(3, ixtg(1 + 3, i))
136 ngl(ii) = ixtg(6, i)
137 nx(ii) = half * ((y2(ii) - y1(ii)) * (z3(ii) - z1(ii)) -
138 . (z2(ii) - z1(ii)) * (y3(ii) - y1(ii)))
139 ny(ii) = zero
140 nz(ii) = zero
141 gbuf%AREA(ii) = sqrt(nx(ii) * nx(ii) + ny(ii) * ny(ii) + nz(ii) * nz(ii))
142 IF (n2d /= 1) THEN
143 gbuf%VOL(ii) = gbuf%AREA(ii)
144 ELSE
145 gbuf%VOL(ii) = (y1(ii) + y2(ii) + y3(ii)) * (
146 . y1(ii) * (z2(ii) - z3(ii)) +
147 . y2(ii) * (z3(ii) - z1(ii)) +
148 . y3(ii) * (z1(ii) - z2(ii))) * one_over_6
149 ENDIF
150
151 lgth1(ii) = sqrt((y2(ii) - y1(ii)) * (y2(ii) - y1(ii)) +
152 . (z2(ii) - z1(ii)) * (z2(ii) - z1(ii)))
153 lgth2(ii) = sqrt((y3(ii) - y2(ii)) * (y3(ii) - y2(ii)) +
154 . (z3(ii) - z2(ii)) * (z3(ii) - z2(ii)))
155 lgth3(ii) = sqrt((y1(ii) - y3(ii)) * (y1(ii) - y3(ii)) +
156 . (z1(ii) - z3(ii)) * (z1(ii) - z3(ii)))
157 gbuf%DELTAX(ii) = gbuf%AREA(ii) /
max(lgth1(ii), lgth2(ii), lgth3(ii))
158 ENDDO
159
160 CALL c3veok3(nvc ,ix1 ,ix2 ,ix3 )
161
162 pm(104,ixtg(1, 1 + nft)) = zero
163
164
165 DO ilay = 1, nlay
166
167 lbuf => elbuf_str%BUFLY(ilay)%LBUF(1,1,1)
168 mbuf => elbuf_str%BUFLY(ilay)%MAT(1,1,1)
169 DO ii = 1, nel
170
171 mat(ii) = mat_param( ixtg(1,ii+nft) )%MULTIMAT%MID(ilay)
172
173 lbuf%VOL(ii) = mat_param( ixtg(1,ii+nft) )%MULTIMAT%VFRAC(ilay
174 ENDDO
175
176 ip = 1
177 ibid = 0
178 CALL matini(pm, ixtg, nixtg, xgrid,
179 . geo, ale_connectivity, detonators, iparg,
180 . sigi, nel, skew, igeo,
181 . ipart,iparttg,
182 . mat, ipm, nsigs, numsh3n, ptsh3n,
183 . ip, ngl, npf, tf, bufmat,
184 . gbuf, lbuf, mbuf, elbuf_str, iloadp,
185 . facload, gbuf%DELTAX,tempel,mat_param )
186
187 vfrac = mat_param( ixtg(1,1+nft) )%MULTIMAT%VFRAC(ilay)
188 pres = pm(104, mat_param( ixtg(1,1+nft) )%MULTIMAT%MID(ilay) )
189 pm(104,ixtg(1, 1 + nft)) = pm(104,ixtg(1, 1 + nft)) + vfrac * pres
190
191 matlaw = ipm(2, mat(1))
192 IF (matlaw == 5) THEN
193
194 IF (.NOT. error_thrown) THEN
195 IF (pm(44, mat(1)) == zero) THEN
196 CALL ancmsg(msgid = 1623, msgtype = msgerror, anmode = aninfo,
197 . i1 = ipm(1, ixtg(1, 1 + nft)), i2 = ipm(1, mat(1)))
198 ENDIF
199 error_thrown = .true.
200 ENDIF
201 CALL m5in2t(pm, mat, ipm(1, ixtg(1,1+nft)), detonators, lbuf%TB, xgrid, ixtg, nixtg)
202 ENDIF
203 ENDDO
204
205 IF (nlay > 1) THEN
206
207
208 gbuf%RHO(1:nel)=zero
209 DO ilay = 1, nlay
210 lbuf => elbuf_str%BUFLY(ilay)%LBUF(1,1,1)
211 DO ii = 1, nel
212 gbuf%RHO(ii) = gbuf%RHO(ii) + lbuf%RHO(ii) * mat_param( ixtg(1,ii+nft) )%MULTIMAT%VFRAC(ilay)
213 ENDDO
214 ENDDO
215
216
217 gbuf%TEMP(1:nel)=zero
218 DO ilay = 1, nlay
219 lbuf => elbuf_str%BUFLY(ilay)%LBUF(1,1,1)
220 DO ii = 1, nel
221 gbuf%TEMP(ii) = gbuf%TEMP(ii) + lbuf%TEMP(ii) *
222 . mat_param( ixtg(1,ii+nft) )%MULTIMAT%VFRAC(ilay)*lbuf%RHO(ii)/gbuf%RHO(ii)
223 ENDDO
224 ENDDO
225
226 ENDIF
subroutine c3veok3(nvc, ix1, ix2, ix3)
subroutine m5in2t(pm, mat, m151_id, detonators, tb, x, ix, nix)
subroutine matini(pm, ix, nix, x, geo, ale_connectivity, detonators, iparg, sigi, nel, skew, igeo, ipart, ipartel, mat, ipm, nsig, nums, pt, ipt, ngl, npf, tf, bufmat, gbuf, lbuf, mbuf, elbuf_str, iloadp, facload, ddeltax, tempel, mat_param)
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)