OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m5in3.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| m5in3 ../starter/source/initial_conditions/detonation/m5in3.F
25!||--- called by ------------------------------------------------------
26!|| m105init ../starter/source/materials/mat/mat105/m105init.F
27!|| m51init ../starter/source/materials/mat/mat051/m51init.F
28!|| m97init ../starter/source/materials/mat/mat097/m97init.F
29!|| matini ../starter/source/materials/mat_share/matini.F
30!|| multifluid_init3 ../starter/source/multifluid/multifluid_init3.F
31!|| multifluid_init3t ../starter/source/multifluid/multifluid_init3t.F
32!||--- calls -----------------------------------------------------
33!|| detcord ../starter/source/initial_conditions/detonation/detcord.F
34!||--- uses -----------------------------------------------------
35!|| detonators_mod ../starter/share/modules1/detonators_mod.F
36!|| message_mod ../starter/share/message_module/message_mod.F
37!||====================================================================
38 SUBROUTINE m5in3(PM,MAT,M151_ID,DETONATORS,TB,IPARG,X,IX,NIX)
39C-----------------------------------------------
40C M o d u l e s
41C-----------------------------------------------
42 USE message_mod
43 USE detonators_mod
44 USE constant_mod , ONLY : ep21
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48#include "implicit_f.inc"
49C-----------------------------------------------
50C G l o b a l P a r a m e t e r s
51C-----------------------------------------------
52#include "mvsiz_p.inc"
53C-----------------------------------------------
54C C o m m o n B l o c k s
55C-----------------------------------------------
56#include "param_c.inc"
57#include "scr03_c.inc"
58#include "units_c.inc"
59#include "vect01_c.inc"
60#include "scr11_c.inc"
61#include "com04_c.inc"
62C-----------------------------------------------
63C D u m m y A r g u m e n t s
64C-----------------------------------------------
65 INTEGER :: MAT(*),IPARG(NPARG),NIX,IX(NIX,*)
66 INTEGER,INTENT(IN) :: M151_ID
67 my_real :: pm(npropm,nummat), tb(*), x(3, *)
68 TYPE(detonators_struct_)::DETONATORS
69C-----------------------------------------------
70C L o c a l V a r i a b l e s
71C-----------------------------------------------
72 INTEGER I, N, MTL,IOPT, NPE2
73 INTEGER NDETPS, NDETSG, NECRAN, NDETPL, NDETCORD, NDETPS_NO_SHADOW, NDETPS_SHADOW
74 LOGICAL HAS_DETONATOR(MVSIZ)
75 my_real
76 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz),
77 . x5(mvsiz), x6(mvsiz), x7(mvsiz), x8(mvsiz),
78 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
79 . y5(mvsiz), y6(mvsiz), y7(mvsiz), y8(mvsiz),
80 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz),
81 . z5(mvsiz), z6(mvsiz), z7(mvsiz), z8(mvsiz),
82 . xc(mvsiz), yc(mvsiz), zc(mvsiz),
83 . bt(mvsiz), dl(mvsiz),
84 . alt , xlp, ylp, zlp, xlp1, ylp1, zlp1,
85 . xlp2,ylp2, zlp2, xl0, yl0, zl0, xl1,
86 . yl1 , zl1, xl2, yl2, zl2, ps1, ps2,
87 . dl1, dl2, s1, s2, s3,
88 . nx, ny, nz, nn, vdet,vdet2
89 INTEGER :: NODE1, NODE2, NODE3, NODE4, NODE5, NODE6, NODE7, NODE8, II
90 INTEGER :: GRNOD_ID, INOD, NNOD, NODE_ID
91 INTEGER :: I_SHADOW_FLAG
92C=======================================================================
93C
94 ndetps = detonators%N_DET_POINT
95 ndetsg = detonators%N_DET_LINE
96 necran = detonators%N_DET_WAVE_SHAPER
97 ndetpl = detonators%N_DET_PLANE
98 ndetcord = detonators%N_DET_CORD
99
100 ndetps = detonators%N_DET_POINT
101 ndetsg = detonators%N_DET_LINE
102 necran = detonators%N_DET_WAVE_SHAPER
103 ndetpl = detonators%N_DET_PLANE
104 ndetcord = detonators%N_DET_CORD
105
106 ! numbering /DFS/DETPOINT without shadowing option (then arrival time deduced from radial distance)
107 ndetps_no_shadow = 0
108 ndetps_shadow = 0
109 DO i = 1, detonators%N_DET_POINT
110 i_shadow_flag = detonators%POINT(i)%SHADOW
111 IF(i_shadow_flag == 0)THEN
112 ndetps_no_shadow = ndetps_no_shadow + 1
113 ELSE
114 ndetps_shadow = ndetps_shadow + 1
115 ENDIF
116 ENDDO
117
118!------------------------------------!
119 !treat only detonators without shadowing option
120 IF(detonators%N_DET - ndetps_shadow > 0) THEN
121
122 IF (jsph == 1) THEN
123C SPH -> only X1, Y1, Z1, have to be used
124 DO i = lft, llt
125 ii = i + nft
126 node1 = ix(3, ii)
127 x1(i) = x(1, node1)
128 x2(i) = zero
129 x3(i) = zero
130 x4(i) = zero
131 x5(i) = zero
132 x6(i) = zero
133 x7(i) = zero
134 x8(i) = zero
135 y1(i) = x(2, node1)
136 y2(i) = zero
137 y3(i) = zero
138 y4(i) = zero
139 y5(i) = zero
140 y6(i) = zero
141 y7(i) = zero
142 y8(i) = zero
143 z1(i) = x(3, node1)
144 z2(i) = zero
145 z3(i) = zero
146 z4(i) = zero
147 z5(i) = zero
148 z6(i) = zero
149 z7(i) = zero
150 z8(i) = zero
151 ENDDO
152 ELSE
153 DO i = lft, llt
154 ii = i + nft
155 node1 = ix(2, ii)
156 node2 = ix(3, ii)
157 node3 = ix(4, ii)
158 node4 = ix(5, ii)
159 node5 = ix(6, ii)
160 node6 = ix(7, ii)
161 node7 = ix(8, ii)
162 node8 = ix(9, ii)
163 x1(i) = x(1, node1)
164 x2(i) = x(1, node2)
165 x3(i) = x(1, node3)
166 x4(i) = x(1, node4)
167 x5(i) = x(1, node5)
168 x6(i) = x(1, node6)
169 x7(i) = x(1, node7)
170 x8(i) = x(1, node8)
171 y1(i) = x(2, node1)
172 y2(i) = x(2, node2)
173 y3(i) = x(2, node3)
174 y4(i) = x(2, node4)
175 y5(i) = x(2, node5)
176 y6(i) = x(2, node6)
177 y7(i) = x(2, node7)
178 y8(i) = x(2, node8)
179 z1(i) = x(3, node1)
180 z2(i) = x(3, node2)
181 z3(i) = x(3, node3)
182 z4(i) = x(3, node4)
183 z5(i) = x(3, node5)
184 z6(i) = x(3, node6)
185 z7(i) = x(3, node7)
186 z8(i) = x(3, node8)
187 ENDDO
188 ENDIF
189
190 !--------------------------------------------------!
191 ! INITIALIZATION BEFORE LIGHTING TIME COMPUTATION !
192 !--------------------------------------------------!
193 DO i=lft,llt
194 tb(i) = ep21
195 has_detonator(i)=.false. ! will be removed -1 for each detonator associated to this elem
196 IF(jsph == 0)THEN
197 IF(iparg(28) == 4) THEN
198 xc(i)=fourth*(x1(i)+x2(i)+x3(i)+x4(i))
199 yc(i)=fourth*(y1(i)+y2(i)+y3(i)+y4(i))
200 zc(i)=fourth*(z1(i)+z2(i)+z3(i)+z4(i))
201 ELSE
202 xc(i)=one_over_8*(x1(i)+x2(i)+x3(i)+x4(i)+x5(i)+x6(i)+x7(i)+x8(i))
203 yc(i)=one_over_8*(y1(i)+y2(i)+y3(i)+y4(i)+y5(i)+y6(i)+y7(i)+y8(i))
204 zc(i)=one_over_8*(z1(i)+z2(i)+z3(i)+z4(i)+z5(i)+z6(i)+z7(i)+z8(i))
205 ENDIF
206 ELSE
207 xc(i)=x1(i)
208 yc(i)=y1(i)
209 zc(i)=z1(i)
210 ENDIF
211 END DO
212
213 !---------------------------------!
214 ! LIGHTING TIME FOR /DFS/DETPOINT !
215 !---------------------------------!
216 IF(ndetps /= 0) THEN
217 DO i=lft,llt
218 DO n=1,ndetps
219 i_shadow_flag = detonators%POINT(n)%SHADOW
220 IF(i_shadow_flag /= 0)cycle ! solve Eikonal equation instead
221 mtl=detonators%POINT(n)%MAT
222 grnod_id=detonators%POINT(n)%GRNOD_ID
223 IF(mtl == 0 .OR. mtl == mat(i) .OR. mtl == m151_id) THEN
224 !--- SINGLE NODE ---!
225 IF(grnod_id == 0)THEN
226 alt=detonators%POINT(n)%TDET
227 xlp=detonators%POINT(n)%XDET
228 ylp=detonators%POINT(n)%YDET
229 zlp=detonators%POINT(n)%ZDET
230 dl(i) =(xc(i)-xlp)**2+(yc(i)-ylp)**2+(zc(i)-zlp)**2
231 dl(i)=sqrt(dl(i))
232 bt(i) =alt+dl(i)/pm(38,mat(i))
233 IF(bt(i) < abs(tb(i))) tb(i)=-bt(i)
234 has_detonator(i)=.true.
235 !--- GROUP OF NODE ---!
236 ELSE
237 nnod = detonators%POINT(n)%NNOD
238 alt=detonators%POINT(n)%TDET
239 has_detonator(i)=.true.
240 DO inod=1,nnod
241 node_id=detonators%POINT(n)%NODLIST(inod)
242 xlp=x(1,node_id)
243 ylp=x(2,node_id)
244 zlp=x(3,node_id)
245 dl(i) =(xc(i)-xlp)**2+(yc(i)-ylp)**2+(zc(i)-zlp)**2
246 dl(i)=sqrt(dl(i))
247 bt(i) =alt+dl(i)/pm(38,mat(i))
248 IF(bt(i) < abs(tb(i))) tb(i)=-bt(i)
249 ENDDO !next INOD
250 ENDIF
251 END IF !(MTL == 0 .OR. MTL == MAT(I) .OR. MTL == M151_ID)
252 END DO
253 END DO
254 ENDIF
255
256 !---------------------------------!
257 ! LIGHTING TIME FOR /DFS/DETLINE !
258 !---------------------------------!
259 IF(ndetsg /= 0) THEN
260 DO n=1,ndetsg
261 alt=detonators%LINE(n)%TDET
262 mtl=detonators%LINE(n)%MAT
263 xlp1=detonators%LINE(n)%XDET_1
264 ylp1=detonators%LINE(n)%YDET_1
265 zlp1=detonators%LINE(n)%ZDET_1
266 xlp2=detonators%LINE(n)%XDET_2
267 ylp2=detonators%LINE(n)%YDET_2
268 zlp2=detonators%LINE(n)%ZDET_2
269 DO i=lft,llt
270 IF(mtl == 0 .OR. mtl == mat(i) .OR. mtl == m151_id) THEN
271 xl0 =(xlp1-xlp2)
272 yl0 =(ylp1-ylp2)
273 zl0 =(zlp1-zlp2)
274 xl1 =(xc(i)-xlp1)
275 yl1 =(yc(i)-ylp1)
276 zl1 =(zc(i)-zlp1)
277 xl2 =(xc(i)-xlp2)
278 yl2 =(yc(i)-ylp2)
279 zl2 =(zc(i)-zlp2)
280 ps1 =xl1*xl0+yl1*yl0+zl1*zl0
281 ps2 =xl2*xl0+yl2*yl0+zl2*zl0
282 IF(ps1*ps2 > zero)THEN
283 dl1 =sqrt(xl1**2+yl1**2+zl1**2)
284 dl2 =sqrt(xl2**2+yl2**2+zl2**2)
285 dl(i)=min(dl1,dl2)
286 ELSE
287 s1 =yl1*zl0 - zl1*yl0
288 s2 =zl1*xl0 - xl1*zl0
289 s3 =xl1*yl0 - yl1*xl0
290 dl(i)=sqrt((s1**2+s2**2+s3**2)/(xl0**2+yl0**2+zl0**2))
291 ENDIF
292 bt(i) =alt+dl(i)/pm(38,mat(i))
293 IF(bt(i) < abs(tb(i))) tb(i)=-bt(i)
294 has_detonator(i)=.true.
295 END IF
296 END DO
297 END DO
298 ENDIF
299
300 !---------------------------------!
301 ! LIGHTING TIME FOR /DFS/DETPLAN !
302 !---------------------------------!
303 IF(ndetpl /= 0) THEN
304 DO n=1,ndetpl
305 alt=detonators%PLANE(n)%TDET
306 mtl=detonators%PLANE(n)%MAT
307 xlp=detonators%PLANE(n)%XDET
308 ylp=detonators%PLANE(n)%YDET
309 zlp=detonators%PLANE(n)%ZDET
310 nx=detonators%PLANE(n)%NX
311 ny=detonators%PLANE(n)%NY
312 nz=detonators%PLANE(n)%NZ
313 nn=sqrt(nx**2+ny**2+nz**2)
314 nn=max(nn,em20)
315 dl1=xlp*nx + ylp*ny + zlp*nz
316 dl1 = dl1/nn
317 DO i=lft,llt
318 IF(mtl == 0 .OR. mtl == mat(i) .OR. mtl == m151_id) THEN
319 ! DL = (OC.ON) / ||ON||
320 ! C: centroid
321 ! P: detonation base
322 ! N: detonation vector
323 dl(i) = (xc(i)-xlp)*nx + (yc(i)-ylp)*ny + (zc(i)-zlp)*nz
324 dl(i) = abs(dl(i))
325 dl(i) = dl(i)/nn
326 bt(i) =alt+dl(i)/pm(38,mat(i))
327 IF(bt(i) < abs(tb(i))) tb(i)=-bt(i)
328 has_detonator(i)=.true.
329 END IF
330 END DO
331 END DO
332 ENDIF
333
334 !---------------------------------!
335 ! LIGHTING TIME FOR /DFS/DETCORD !
336 !---------------------------------!
337 vdet = zero
338 IF(ndetcord /= 0) THEN
339 DO n=1,ndetcord
340 alt = detonators%CORD(n)%TDET
341 mtl = detonators%CORD(n)%MAT
342 vdet2 = detonators%CORD(n)%VDET
343 iopt = detonators%CORD(n)%IOPT
344 npe2 = detonators%CORD(n)%NUMNOD
345 dto0 = alt
346 vdto=pm(38,mat(1))
347 IF(vdet == zero)vdet=pm(38,mat(1)) !optional detonation velocity
348 IF(mtl /= mat(1) .AND. mtl /= 0 .AND. mtl /= m151_id) cycle
349 dto0 = alt
350 CALL detcord(detonators%CORD(n),x,xc,yc,zc,vdto,vdet2,alt,bt,tb,has_detonator,iopt)
351 END do! next N
352 ENDIF
353
354 !--------------------------------------!
355 END IF !-> LIGHTING TIME COMPUTATION COMPLETE !
356 !--------------------------------------!
357
358 RETURN
359 END
360
my_real
#define my_real
Definition cppsort.cpp:32
detcord
subroutine detcord(detonator_cord, x, xc, yc, zc, vdet, vdet2, alt, bt, tb, has_detonator, iopt)
Definition detcord.F:34
m5in3
subroutine m5in3(pm, mat, m151_id, detonators, tb, iparg, x, ix, nix)
Definition m5in3.F:39
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
constant_mod
Definition constant_mod.F:294
detonators_mod
Definition detonators_mod.F:79
message_mod
Definition message_mod.F:1249
detonators_mod::detonators_struct_
Definition detonators_mod.F:154