OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
detcord.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "vect01_c.inc"
#include "scr11_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine detcord (detonator_cord, x, xc, yc, zc, vdet, vdet2, alt, bt, tb, has_detonator, iopt)

Function/Subroutine Documentation

◆ detcord()

subroutine detcord ( type(detonator_cord_struct_) detonator_cord,
x,
xc,
yc,
zc,
vdet,
vdet2,
alt,
bt,
tb,
logical, dimension(mvsiz) has_detonator,
integer iopt )

Definition at line 33 of file detcord.F.

34C-----------------------------------------------
35C D e s c r i p t i o n
36C-----------------------------------------------
37 !IOPT=0 : def=3
38 !IOPT=1 : piecewise linear - multiple segments (experimental / osbslete)
39 !IOPT=2 : instantaneous - multiple segments (experimental / obsolete)
40 !IOPT=3 : Centripetal-Catmull-Rom SPLINE interpolation + projection along neutral fiber
41C-----------------------------------------------
42C M o d u l e s
43C-----------------------------------------------
44 USE detonators_mod
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48#include "implicit_f.inc"
49C-----------------------------------------------
50C C o m m o n B l o c k s
51C-----------------------------------------------
52#include "mvsiz_p.inc"
53#include "vect01_c.inc"
54#include "scr11_c.inc"
55C-----------------------------------------------
56C D u m m y A r g u m e n t s
57C-----------------------------------------------
58 INTEGER :: IOPT
59 my_real :: x(3,*),xc(mvsiz),yc(mvsiz),zc(mvsiz),bt(mvsiz),vdet,vdet2,alt,tb(mvsiz)
60 TYPE(DETONATOR_CORD_STRUCT_)::DETONATOR_CORD
61 LOGICAL HAS_DETONATOR(MVSIZ)
62C-----------------------------------------------
63C L o c a l V a r i a b l e s
64C-----------------------------------------------
65 INTEGER :: J, I1,I2, II,K, iTdet, i, NPTS, NSEG,IDIST(MVSIZ),NDETCORD,NNOD
66 my_real :: ddmx, dtos, xlp2, ylp2, xlp1 , zlp2 ,ylp1,
67 . zlp1, xl0 , yl0 , zl0 , xl1 , yl1 ,zl1 , xl2 ,yl2 ,zl2,
68 . ps1 , ps2 , dl1 , dl2 , dl(mvsiz), s1 ,s2 , s3, tdetc, tp1, tp2, th, th1,th2,
69 . tc1,tc2,tc0,tc,p1p2,dh,alpha,knots(4),zp(3),
70 . zh1(3),dist1,t1,zh2(3),dist2,t2,dist(mvsiz),xx,yy,zz,d,local_pt(4,3),c(3),t,dd,
71 . len,len1,len2
72
73 type spline_path
74 INTEGER :: NUM_POINTS
75 my_real,ALLOCATABLE,DIMENSION(:,:) :: control_point
76 my_real,ALLOCATABLE,DIMENSION(:) :: length
77 my_real,ALLOCATABLE,DIMENSION(:) :: cumulative_length
78 my_real,ALLOCATABLE,DIMENSION(:,:) :: knots
79 end type spline_path
80
81 type(SPLINE_PATH),TARGET :: USER_SPLINE_PATH
82
83 my_real, POINTER, DIMENSION(:,:) :: ptr
84
85C-----------------------------------------------
86C S o u r c e L i n e s
87C-----------------------------------------------
88 nnod = detonator_cord%NUMNOD
89 dtos = ep20
90
91 IF(iopt == 2)THEN
92 !INSTANTANEOUS IGNITION OF THE WHOLE CORD
93 DO j=1,nnod-1
94 !first node
95 ii = detonator_cord%NODES(j)
96 xlp1 = x(1,ii)
97 ylp1 = x(2,ii)
98 zlp1 = x(3,ii)
99 !second node
100 ii = detonator_cord%NODES(j+1)
101 xlp2 = x(1,ii)
102 ylp2 = x(2,ii)
103 zlp2 = x(3,ii)
104 !
105 xl0 = (xlp1-xlp2)
106 yl0 = (ylp1-ylp2)
107 zl0 = (zlp1-zlp2)
108 !loop on elems to compute burning time from line [P1,P2]
109 DO i=lft,llt
110 xl1 = (xc(i)-xlp1)
111 yl1 = (yc(i)-ylp1)
112 zl1 = (zc(i)-zlp1)
113 xl2 = (xc(i)-xlp2)
114 yl2 = (yc(i)-ylp2)
115 zl2 = (zc(i)-zlp2)
116 ps1 = xl1*xl0+yl1*yl0+zl1*zl0
117 ps2 = xl2*xl0+yl2*yl0+zl2*zl0
118 IF(ps1*ps2 > zero)THEN
119 !projection point not on the segment [P1,P2]
120 dl1 = sqrt(xl1**2+yl1**2+zl1**2)
121 dl2 = sqrt(xl2**2+yl2**2+zl2**2)
122 dl(i)=min(dl1,dl2)
123 ELSE
124 !projection along the segment [P1,P2]
125 s1 = yl1*zl0 - zl1*yl0
126 s2 = zl1*xl0 - xl1*zl0
127 s3 = xl1*yl0 - yl1*xl0
128 dl(i)=sqrt((s1**2+s2**2+s3**2)/(xl0**2+yl0**2+zl0**2))
129 ENDIF
130 bt(i) =alt+dl(i)/vdet
131 IF(bt(i) < abs(tb(i))) tb(i)=-bt(i)
132 END DO !next I
133 END DO !J=1,NNOD
134
135 ELSEIF(iopt == 1)THEN
136 ! DETONATING CORD HAS ITS OWN DETONATION VELOCITY ALONG THE PATH
137 DO j=1,nnod-1
138 !first node
139 ii = detonator_cord%NODES(j)
140 xlp1 = x(1,ii)
141 ylp1 = x(2,ii)
142 zlp1 = x(3,ii)
143 tp1 = detonator_cord%TDET_PATH(j)
144 !second node
145 ii = detonator_cord%NODES(j+1)
146 xlp2 = x(1,ii)
147 ylp2 = x(2,ii)
148 zlp2 = x(3,ii)
149 tp2 = detonator_cord%TDET_PATH(j+1)
150 !loop on elems to compute burning time from line [P1,P2]
151 DO i=lft,llt
152 tdetc = zero
153 xl0 = (xlp1-xlp2)
154 yl0 = (ylp1-ylp2)
155 zl0 = (zlp1-zlp2)
156 xl1 = (xc(i)-xlp1)
157 yl1 = (yc(i)-ylp1)
158 zl1 = (zc(i)-zlp1)
159 xl2 = (xc(i)-xlp2)
160 yl2 = (yc(i)-ylp2)
161 zl2 = (zc(i)-zlp2)
162 ps1 = xl1*xl0+yl1*yl0+zl1*zl0
163 ps2 = xl2*xl0+yl2*yl0+zl2*zl0
164 dl1 = sqrt(xl1**2+yl1**2+zl1**2)
165 dl2 = sqrt(xl2**2+yl2**2+zl2**2)
166 tc1 = tp1 + dl1 /vdet
167 tc2 = tp2 + dl2 /vdet
168 tc = min(tc1,tc2)
169 IF(ps1*ps2 <= zero)THEN
170 s1 = yl1*zl0 - zl1*yl0
171 s2 = zl1*xl0 - xl1*zl0
172 s3 = xl1*yl0 - yl1*xl0
173 !
174 p1p2 = (xl0**2+yl0**2+zl0**2)
175 dl2 = (s1**2+s2**2+s3**2)/p1p2
176 dl(i) = sqrt(dl2)
177 p1p2 = sqrt(p1p2)
178 dh = sqrt((xl1**2+yl1**2+zl1**2)-dl2) !P1H
179 !
180 th1 = tp1 + dh/vdet2
181 th2 = tp2 + (p1p2-dh)/vdet2
182 th = min(th1,th2)
183 tc0 = th + dl(i)/vdet
184 tc = min(tc,tc0)
185 ENDIF
186 bt(i) = tc
187 IF(bt(i) < abs(tb(i))) tb(i)=-bt(i)
188 END DO !next I
189 END DO !J=1,NNOD
190
191 ELSEIF(iopt == 3)THEN
192
193 !CATMULL-ROM PARAMETER
194 ! ALPHA = 0.0 : uniform
195 ! ALPHA = 0.5 : centripetal
196 ! ALPHA = 1.0 : chordal
197 alpha = half
198
199 IF(vdet2 == zero)vdet2=vdet
200
201 !initialize burning times by following the neutral fiber.
202 npts = nnod !number of point composing the cord path
203 nseg = nnod-1 !number of segment composing the cord path
204 ALLOCATE(user_spline_path%control_point(0:npts+1,3))
205 ALLOCATE(user_spline_path%length(npts-1))
206 ALLOCATE(user_spline_path%cumulative_length(npts-1))
207 ALLOCATE(user_spline_path%knots(npts-1,4) )
208 !PT(0 ,1:3)= will be defined by energy minimum
209 !PT(NPTS+1,1:3)= will be defined by energy minimum
210
211 !===CONTROL POINTS===!
212 DO j=1,nnod
213 ii = detonator_cord%NODES(j)
214 xlp1 = x(1,ii)
215 ylp1 = x(2,ii)
216 zlp1 = x(3,ii)
217 user_spline_path%control_point(j,1:3)=(/xlp1, ylp1, zlp1/)
218 ENDDO
219
220 !===END POINTS===!
221 ptr=>user_spline_path%control_point(0:npts+1,1:3)
222 !ADDING FIRST END POINTS - BENDING ENERGY MINIMUM
223 ptr(1+0,1)=half*(five*ptr(1+1,1)-four*ptr(1+2,1)+ptr(1+3,1))
224 ptr(1+0,2)=half*(five*ptr(1+1,2)-four*ptr(1+2,2)+ptr(1+3,2))
225 ptr(1+0,3)=half*(five*ptr(1+1,3)-four*ptr(1+2,3)+ptr(1+3,3))
226 !ADDING LAST END POINTS - BENDING ENERGY MINIMUM
227 ptr(1+npts+1,1)=half*(1.*ptr(1+npts-2,1)-four*ptr(1+npts-1,1)+five*ptr(1+npts,1))
228 ptr(1+npts+1,2)=half*(1.*ptr(1+npts-2,2)-four*ptr(1+npts-1,2)+five*ptr(1+npts,2))
229 ptr(1+npts+1,3)=half*(1.*ptr(1+npts-2,3)-four*ptr(1+npts-1,3)+five*ptr(1+npts,3))
230
231 !DEFINE KNOTS
232 DO j=1,nseg
233 CALL cr_spline_knots(ptr(j ,1), knots, alpha)
234 user_spline_path%knots(j,1:4)=knots(1:4)
235 ENDDO
236
237 !COMPUTE DISTANCE FROM CENTROID Z TO CONTROL POINT ptr(1+:)
238 ! idist(J) is the shortest point in PT from Z(J)
239 DO j=lft,llt
240 dist(j) = 1e20
241 idist(j) = 0
242 DO i=1,npts
243 xx = xc(j)-ptr(1+i,1)
244 yy = yc(j)-ptr(1+i,2)
245 zz = zc(j)-ptr(1+i,3)
246 dd = sqrt(xx*xx+yy*yy+zz*zz)
247 IF(dd < dist(j))THEN
248 dist(j) = dd
249 idist(j) = i
250 ENDIF
251 ENDDO
252 ENDDO
253
254 !SPLINE LENGTHS
255 DO i=1,nseg
256 local_pt(1,1:3)=ptr(1+i-1,1:3)
257 local_pt(2,1:3)=ptr(1+i+0,1:3)
258 local_pt(3,1:3)=ptr(1+i+1,1:3)
259 local_pt(4,1:3)=ptr(1+i+2,1:3)
260 t=1.0
261 CALL cr_spline_length(local_pt,alpha,t,len)
262 user_spline_path%length(i) = len
263 IF(i > 1)THEN
264 user_spline_path%cumulative_length(i) = user_spline_path%cumulative_length(i-1) + len
265 ELSE
266 user_spline_path%cumulative_length(1) = len
267 ENDIF
268 !print *, "spline-i=", I
269 !print *, "len =", len
270 ENDDO
271
272 !TEST POINT PROJECTIONS
273 DO j=lft,llt
274 !nearest point : IDIST(J)
275 i=idist(j)
276 i=min(i,nseg)
277 k=i
278 t=half
279 local_pt(1,1:3)=ptr(1+i-1,1:3)
280 local_pt(2,1:3)=ptr(1+i+0,1:3)
281 local_pt(3,1:3)=ptr(1+i+1,1:3)
282 local_pt(4,1:3)=ptr(1+i+2,1:3)
283 zp(1:3) =(/xc(j),yc(j),zc(j)/)
284 CALL cr_spline_point_proj(local_pt,zp(1:3),alpha,zh1,dist1,t1)
285 CALL cr_spline_length(local_pt,alpha,t1,len1)
286 c(1:3)=zh1(1:3)
287 t=t1 ! position [0,1]
288 len=len1
289 !two adjacent spline on given node IDIST(J), unless for first spline I==1
290 IF(i >= 2)THEN
291 i=i-1
292 local_pt(1,1:3)=ptr(1+i-1,1:3)
293 local_pt(2,1:3)=ptr(1+i+0,1:3)
294 local_pt(3,1:3)=ptr(1+i+1,1:3)
295 local_pt(4,1:3)=ptr(1+i+2,1:3)
296 CALL cr_spline_point_proj(local_pt,zp(1:3),alpha,zh2,dist2,t2)
297 CALL cr_spline_length(local_pt,alpha,t2,len2)
298 IF(dist2 < dist1)THEN
299 c(1:3)=zh2(1:3)
300 t=t2
301 k=i
302 len = len2
303 ENDIF
304 ENDIF
305 !print *, "point=", ZP(1:3)
306 !print *, " ->", C(1:3)
307 !!CALL CR_SPLINE_LENGTH(LOCAL_PT,ALPHA,T,LEN)
308 IF(k > 1)len=len+user_spline_path%cumulative_length(k-1)
309 bt(j) = detonator_cord%TDET_PATH(1) + len/vdet2
310 IF(bt(j) < abs(tb(j))) tb(j)=-bt(j)
311 ENDDO
312
313 ENDIF
314
315 DO i=lft,llt
316 has_detonator(i) = .true.
317 ENDDO
318
319 dto=dtos
320
321C-----------------------------------------------
322 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
cr_spline_knots
subroutine cr_spline_knots(pts, knots, alpha)
Definition cr_spline_knots.F:31
cr_spline_length
subroutine cr_spline_length(local_pt, alpha, t, len)
Definition cr_spline_length.F:32
cr_spline_point_proj
subroutine cr_spline_point_proj(local_pt, z, alpha, zh, h, t)
Definition cr_spline_point_proj.F:32
the
end diagonal values have been computed in the(sparse) matrix id.SOL
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
detonators_mod
Definition detonators_mod.F:79