OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pblast_2.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!|| pblast_2 ../engine/source/loads/pblast/pblast_2.F
25!||--- called by ------------------------------------------------------
26!|| pblast_load_computation ../engine/source/loads/pblast/pblast.F
27!||--- calls -----------------------------------------------------
28!|| arret ../engine/source/system/arret.F
29!|| my_barrier ../engine/source/system/machine.F
30!|| pblast_parameters__surface_burst ../common_source/modules/loads/pblast_mod.F90
31!||--- uses -----------------------------------------------------
32!|| groupdef_mod ../common_source/modules/groupdef_mod.F
33!|| h3d_inc_mod ../engine/share/modules/h3d_inc_mod.F
34!|| h3d_mod ../engine/share/modules/h3d_mod.F
35!|| names_and_titles_mod ../common_source/modules/names_and_titles_mod.F
36!|| output_mod ../common_source/modules/output/output_mod.F90
37!|| pblast_mod ../common_source/modules/loads/pblast_mod.F90
38!|| th_mod ../engine/share/modules/th_mod.F
39!|| th_surf_mod ../common_source/modules/interfaces/th_surf_mod.F
40!||====================================================================
41 SUBROUTINE pblast_2(PBLAST,
42 * ILOADP ,FAC ,A ,V ,X ,
43 1 IADC ,FSKY ,LLOADP ,FEXT ,NODA_SURF,NODA_PEXT,
44 2 ITAB ,H3D_DATA ,NL ,DTMIN_LOC ,WFEXT_LOC,
45 3 TH_SURF ,NSEGPL)
46C-----------------------------------------------
47C M o d u l e s
48C-----------------------------------------------
49 USE h3d_mod
50 USE pblast_mod
51 USE groupdef_mod
52 USE h3d_inc_mod
53 USE th_surf_mod , ONLY : th_surf_
54 USE names_and_titles_mod , ONLY : ncharline
55 USE th_mod ,ONLY : th_has_noda_pext
56 USE output_mod , ONLY : h3d_has_noda_pext, anim_has_noda_pext
57C-----------------------------------------------
58C I m p l i c i t T y p e s
59C-----------------------------------------------
60#include "implicit_f.inc"
61#include "comlock.inc"
62#include "param_c.inc"
63C-----------------------------------------------
64C C o m m o n B l o c k s
65C-----------------------------------------------
66#include "com04_c.inc"
67#include "com08_c.inc"
68#include "parit_c.inc"
69#include "scr14_c.inc"
70#include "scr16_c.inc"
71#include "mvsiz_p.inc"
72#include "units_c.inc"
73#include "sysunit.inc"
74#include "tabsiz_c.inc"
75C-----------------------------------------------
76C D u m m y A r g u m e n t s
77C-----------------------------------------------
78 INTEGER,INTENT(IN) :: LLOADP(SLLOADP)
79 INTEGER,INTENT(INOUT) :: ILOADP(SIZLOADP,NLOADP)
80 INTEGER,INTENT(IN) :: IADC(*)
81 INTEGER, INTENT(IN) :: ITAB(NUMNOD),NL
82 my_real,INTENT(INOUT) :: dtmin_loc
83 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT_LOC
84 my_real,INTENT(IN) :: v(3,numnod),x(3,numnod)
85 my_real,INTENT(INOUT) :: fac(lfacload,nloadp)
86 my_real,INTENT(INOUT) :: a(3,numnod),fsky(8,sfsky/8)
87 my_real,INTENT(INOUT) :: fext(3,numnod)
88 my_real,INTENT(INOUT) :: noda_surf(numnod)
89 my_real,INTENT(INOUT) :: noda_pext(numnod)
90 TYPE(h3d_database),INTENT(IN) :: H3D_DATA
91 TYPE (TH_SURF_) , INTENT(INOUT) :: TH_SURF
92 INTEGER, INTENT(INOUT) :: NSEGPL
93 TYPE(pblast_),INTENT(INOUT) :: PBLAST
94C-----------------------------------------------
95C L o c a l V a r i a b l e s
96C-----------------------------------------------
97 INTEGER :: N1, N2, N3, N4, IL, IS, IAD, I, IANIM_OR_H3D,IZ_UPDATE,ABAC_ID,ISIZ_SEG,IERR1,ID, ITA_SHIFT,NS,KSURF
98 INTEGER :: NDT,IMODEL,NN(4)
99 INTEGER :: ISHAPE
100 INTEGER :: NWARN !< number of segment for which it is not possible to correlate the positive impulse. It may happen that for a given Pmax and dt0 value even building a triangle shape leads to a lower impulse that the targeted one.
101 my_real :: zx,zy,zz,norm,nx,ny,nz,nnx,nny,nnz,hz,area
102 my_real :: lambda,cos_theta, alpha_inci, alpha_refl, p_inci,p_refl,z
103 my_real :: i_inci,i_refl,dt_0,t_a,wave_refl,wave_inci, w13
104 my_real :: xdet,ydet,zdet,tdet,wtnt,pmin,t_stop,p,fac_m_bb, fac_l_bb, fac_t_bb, fac_p_bb, fac_i_bb, ta_first, tt_star
105 my_real :: z1_
106 my_real :: decay_inci,decay_refl
107 my_real :: cst_255_div_ln_z1_on_zn, log10_, npt, ff(3)
108 my_real :: projz(3), tmp(3), lg, zg, hc
109 my_real :: base_x,base_y,base_z
110 my_real :: ta_inf_loc,zeta
111 my_real :: dnorm, dx, dy, dz
112 my_real :: surf_patch
113 TYPE(FRIEDLANDER_PARAMS_) :: FRIEDLANDER_PARAMS
114 LOGICAL,SAVE :: IS_UPDATED
115 LOGICAL :: IS_DECAY_TO_BE_COMPUTED
116
117 CHARACTER(LEN=NCHARLINE) :: MSGOUT1,MSGOUT2
118
119 DATA cst_255_div_ln_z1_on_zn/-38.147316611455952998/
120 DATA log10_ /2.30258509299405000000/
121
122C-----------------------------------------------
123C D e s c r i p t i o n
124C-----------------------------------------------
125C This subroutines is applying pressure load to a segment submitted to a blast wave (SURFACE BURST model).
126C Preussre time histories are built from "UFC 3-340-02, Dec. 5th 2008" tables which are hardcoded in unit system {g, cm, mus}
127C (table sampling : 256 points ; see pblast_mod.F)
128C
129C T* = T + TA_INF (shift with first arrival time for all pblast option, TA_INF=0. if ita_shift/=2)
130C If request made to update blast profile (iz_update==2) then it is made once on T*=TDET
131C
132C-----------------------------------------------
133C P r e - C o n d i t i o n
134C-----------------------------------------------
135 !--- subroutine not relevant if no /LOAD/PBLAST option
136 IF(pblast%NLOADP_B == 0)RETURN
137
138 !--- time step
139 ta_first = fac(07,nl) ! arrival time for first segment to be loaded
140 il= nl-nloadp_f
141 tt_star = tt + pblast%PBLAST_DT%TA_INF
142 iz_update = iloadp(06,nl)
143 tdet = fac(01,nl)
144 ta_first = fac(07,nl) + pblast%PBLAST_DT%TA_INF ! arrival time for first segment to be loaded
145 IF(iz_update ==1)THEN
146 !if no update during engine
147 dtmin_loc = (one+em06)*(ta_first - tt)
148 dtmin_loc=max(pblast%PBLAST_TAB(il)%DTMIN, dtmin_loc) !go directly to arrival time but avoid too
149 IF(tt_star<ta_first)RETURN ! (nothing to load for now)
150 ELSE
151 IF(tdet > tt)THEN
152 dtmin_loc = (one+em06)*(tdet - tt)
153 dtmin_loc=max(pblast%PBLAST_TAB(il)%DTMIN, dtmin_loc)!go directly to arrival time but avoid too
154 ELSE
155 dtmin_loc = pblast%PBLAST_TAB(il)%DTMIN
156 ENDIF
157 IF(tt_star<tdet)RETURN ! (nothing to update for now)
158 ENDIF
159C-----------------------------------------------,
160C S o u r c e C o d e
161C-----------------------------------------------
162 ianim_or_h3d = anim_v(5)+outp_v(5)+h3d_data%N_VECT_FINT + anim_v(6)+outp_v(6)+h3d_data%N_VECT_FEXT !output
163
164 !Index Bijection
165 z1_ = 0.500000000000000
166
167 !translation from Working unit System to {big bang} unit system
168 fac_m_bb = fac_mass*ep03
169 fac_l_bb = fac_length*ep02
170 fac_t_bb = fac_time*ep06
171 fac_p_bb = fac_m_bb/fac_l_bb/fac_t_bb/fac_t_bb
172 fac_i_bb = fac_p_bb*fac_t_bb
173 fac_i_bb = fac_m_bb/fac_l_bb/fac_t_bb
174
175 is_updated = .false.
176 CALL my_barrier
177
178 !-----------------------------------------------
179 ! FREE AIR BURST
180 !-----------------------------------------------
181 il = nl-nloadp_f
182 tdet = fac(01,nl)
183 xdet = fac(02,nl)
184 ydet = fac(03,nl)
185 zdet = fac(04,nl)
186 wtnt = fac(05,nl)
187 pmin = fac(06,nl)
188 ta_first = fac(07,nl)
189 nnx = fac(08,nl)
190 nny = fac(09,nl)
191 nnz = fac(10,nl)
192 hc = fac(11,nl)
193 t_stop = fac(13,nl)
194 is = iloadp(02,nl)
195 ishape = iloadp(03,nl)
196 iz_update = iloadp(06,nl)
197 abac_id = iloadp(07,nl)
198 id = iloadp(08,nl) !user_id
199 ita_shift = iloadp(09,nl)
200 ndt = iloadp(10,nl)
201 imodel = iloadp(11,nl)
202 isiz_seg = iloadp(01,nl)/4
203 ierr1 = 0
204 w13 = (wtnt*fac_m_bb)**third ! '*FAC_M' g->work unit '/FAC_M' : WORK_UNIT -> g
205 z = zero
206 ta_inf_loc = ep20
207
208 is_decay_to_be_computed = .false.
209 IF(imodel == 2)is_decay_to_be_computed=.true.
210 nwarn = 0
211
212 !---------------------------------------------
213 ! LOOP ON SEGMENTS (4N or 3N)
214 !---------------------------------------------
215
216!$OMP DO SCHEDULE(GUIDED,MVSIZ)
217 DO i = 1,isiz_seg
218 n1=lloadp(iloadp(4,nl)+4*(i-1))
219 n2=lloadp(iloadp(4,nl)+4*(i-1)+1)
220 n3=lloadp(iloadp(4,nl)+4*(i-1)+2)
221 n4=lloadp(iloadp(4,nl)+4*(i-1)+3)
222 nn(1)=n1
223 nn(2)=n2
224 nn(3)=n3
225 nn(4)=n4
226
227 IF(n4 == 0 .OR. n3 == n4 )THEN
228 !3-NODE-SEGMENT
229 pblast%PBLAST_TAB(il)%NPt(i) = three
230 npt = three
231 !Segment Centroid
232 zx = x(1,n1)+x(1,n2)+x(1,n3)
233 zy = x(2,n1)+x(2,n2)+x(2,n3)
234 zz = x(3,n1)+x(3,n2)+x(3,n3)
235 zx = zx*third
236 zy = zy*third
237 zz = zz*third
238 !Normal vector : (NX,NY,NZ) = 2*S*n where |n|=1.0
239 nx = (x(2,n3)-x(2,n1))*(x(3,n3)-x(3,n2)) - (x(3,n3)-x(3,n1))*(x(2,n3)-x(2,n2))
240 ny = (x(3,n3)-x(3,n1))*(x(1,n3)-x(1,n2)) - (x(1,n3)-x(1,n1))*(x(3,n3)-x(3,n2))
241 nz = (x(1,n3)-x(1,n1))*(x(2,n3)-x(2,n2)) - (x(2,n3)-x(2,n1))*(x(1,n3)-x(1,n2))
242 !NORM = 2*S
243 norm = sqrt(nx*nx+ny*ny+nz*nz)
244 ELSE
245 !4-NODE-SEGMENT
246 pblast%PBLAST_TAB(il)%NPt(i) = four
247 npt = four
248 !Face Centroid
249 zx = x(1,n1)+x(1,n2)+x(1,n3)+x(1,n4)
250 zy = x(2,n1)+x(2,n2)+x(2,n3)+x(2,n4)
251 zz = x(3,n1)+x(3,n2)+x(3,n3)+x(3,n4)
252 zx = zx*fourth
253 zy = zy*fourth
254 zz = zz*fourth
255 !Normal vector (NX,NY,NZ) = 2*S*n where |n|=1.0
256 nx = (x(2,n3)-x(2,n1))*(x(3,n4)-x(3,n2)) - (x(3,n3)-x(3,n1))*(x(2,n4)-x(2,n2))
257 ny = (x(3,n3)-x(3,n1))*(x(1,n4)-x(1,n2)) - (x(1,n3)-x(1,n1))*(x(3,n4)-x(3,n2))
258 nz = (x(1,n3)-x(1,n1))*(x(2,n4)-x(2,n2)) - (x(2,n3)-x(2,n1))*(x(1,n4)-x(1,n2))
259 !NORM = 2*S
260 norm = sqrt(nx*nx+ny*ny+nz*nz)
261 ENDIF
262
263 pblast%PBLAST_TAB(il)%N(1,i) = n1
264 pblast%PBLAST_TAB(il)%N(2,i) = n2
265 pblast%PBLAST_TAB(il)%N(3,i) = n3
266 pblast%PBLAST_TAB(il)%N(4,i) = n4
267
268 !Determine Height of centroid (structure face)
269 ! Basis = DETPOINT
270 ! NN is ground vector
271 hz = ( nnx*zx + nny*zy + nnz*zz - nnx*xdet - nny*ydet - nnz*zdet )
272
273 !--------------------------------------------!
274 ! Update Wave parameters !
275 ! (experimental) !
276 ! If requested. Otherwise use Starter param. !
277 ! Default : do not update:use Starter param. !
278 !--------------------------------------------!
279 IF(iz_update == 2)THEN
280
281 dtmin_loc = ep20
282
283 IF(hz >= -em10 .OR. ishape == 2)THEN
284
285 !OVER GROUND
286
287 !Planar Wave : angle with targeted face
288 base_x = xdet
289 base_y = ydet
290 base_z = zdet
291 lambda = (base_x-zx)*nnx + (base_y-zy)*nny + (base_z-zz)*nnz
292
293 projz(1) = zx + lambda*nnx
294 projz(2) = zy + lambda*nny
295 projz(3) = zz + lambda*nnz
296
297 tmp(1) = (projz(1)-xdet)
298 tmp(2) = (projz(2)-ydet)
299 tmp(3) = (projz(3)-zdet)
300 lg = sqrt(tmp(1)*tmp(1)+tmp(2)*tmp(2)+tmp(3)*tmp(3))
301 zg = lg*fac_l_bb/w13 !scaled horizontal distance (ground) {g,cm,mus}
302
303
304 ! scaled distance
305 cos_theta = 0
306 SELECT CASE(ishape)
307 CASE (1,2)
308 ! Dist
309 dx = (xdet - zx)*fac_l_bb
310 dy = (ydet - zy)*fac_l_bb
311 dz = (zdet - zz)*fac_l_bb
312 dnorm = sqrt(dx*dx+dy*dy+dz*dz) ! *FAC_L_bb cm->work unit /FAC_L_bb : WORK_UNIT -> cm
313 z = dnorm / w13 !spherical
314 !Angle from detonation point
315 cos_theta = dx*nx + dy*ny + dz*nz
316 cos_theta = cos_theta/max(em20,norm*dnorm)
317 CASE(3)
318 z = zg !cylindrical - in abac unit id g,cm,mus
319 ! Angle from detonation point
320 cos_theta = (xdet-projz(1))*nx + (ydet-projz(2))*ny + (zdet-projz(3))*nz
321 cos_theta = cos_theta/max(em20,lg*norm)
322 END SELECT
323
324 !checking bounds
325 IF(z > 0.5 .AND. z < 400.) then
326 !Ground Distance inside the interpolation table
327
328 elseif(z <= 0.5 .AND. pblast%PBLAST_TAB(il)%TAGMSG(i) == 0)then
329 write(iout, fmt='(A,I0,A)') "Warning : /LOAD/PBLAST id=",id," Rg/W**(1/3) < 0.5 cm/g**(1/3)"
330 write(istdo,fmt='(A,I0,A)') "Warning : /LOAD/PBLAST id=",id," Rg/W**(1/3) < 0.5 cm/g**(1/3)"
331 write(iout, fmt='(A)') " Horizontal Distance on Ground (Rg) is too close to the charge"
332 write(istdo,fmt='(A)') " Horizontal Distance on Ground (Rg) is too close to the charge"
333 if (n4 == 0 .OR. n3 == n4)then
334 write(iout, fmt='(A,3I11)') " -> Segment nodes : ",itab(n1),itab(n2),itab(n3)
335 write(istdo,fmt='(A,3I11)') " -> Segment nodes : ",itab(n1),itab(n2),itab(n3)
336 else
337 write(iout, fmt='(A,4I11)') " -> Segment nodes : ",itab(n1),itab(n2),itab(n3),itab(n4)
338 write(istdo,fmt='(A,4I11)') " -> Segment nodes : ",itab(n1),itab(n2),itab(n3),itab(n4)
339 endif
340 pblast%PBLAST_TAB(il)%TAGMSG(i) = 1
341
342 elseif(z > 400. .AND. pblast%PBLAST_TAB(il)%TAGMSG(i) == 0)then
343 write(iout, fmt='(A,I0,A)') "Warning : /LOAD/PBLAST id=",id," Rg/W**(1/3) > 400. cm/g**(1/3)"
344 write(istdo,fmt='(a,i0,a)') "warning : /load/pblast id=",id," rg/w**(1/3) > 400. cm/g**(1/3)"
345 write(iout, fmt='(A)') " Horizontal Distance on Ground (Rg) is too far from the charge"
346 write(istdo,fmt='(A)') " Horizontal Distance on Ground (Rg) is too far from the charge"
347 if (n4 == 0 .OR. n3 == n4)then
348 write(iout, fmt='(A,3I11)') " -> Segment nodes : ",itab(n1),itab(n2),itab(n3)
349 write(istdo,fmt='(A,3I11)') " -> Segment nodes : ",itab(n1),itab(n2),itab(n3)
350 else
351 write(iout, fmt='(A,4I11)') " -> Segment nodes : ",itab(n1),itab(n2),itab(n3),itab(n4)
352 write(istdo,fmt='(A,4I11)') " -> Segment nodes : ",itab(n1),itab(n2),itab(n3),itab(n4)
353 endif
354 pblast%PBLAST_TAB(il)%TAGMSG(i) = 1
355
356 endif
357
358 !------------------------------------------------------------------!
359 CALL pblast_parameters__surface_burst( pblast,
360 + z, w13, tdet,
361 + fac_p_bb, fac_i_bb, fac_t_bb,
362 + is_decay_to_be_computed,
363 + friedlander_params, nwarn )
364 p_inci = friedlander_params%P_inci
365 p_refl = friedlander_params%P_refl
366 i_inci = friedlander_params%I_inci
367 i_refl = friedlander_params%I_refl
368 t_a = friedlander_params%T_A
369 dt_0 = friedlander_params%DT_0
370 decay_inci = friedlander_params%decay_inci
371 decay_refl = friedlander_params%decay_refl
372 !------------------------------------------------------------------!
373
374 ta_inf_loc = min(ta_inf_loc, t_a)
375
376 !update wave parameters
377 pblast%PBLAST_TAB(il)%cos_theta(i) = cos_theta
378 pblast%PBLAST_TAB(il)%P_inci(i) = p_inci
379 pblast%PBLAST_TAB(il)%P_refl(i) = p_refl
380 pblast%PBLAST_TAB(il)%ta(i) = t_a
381 pblast%PBLAST_TAB(il)%t0(i) = dt_0
382 pblast%PBLAST_TAB(il)%decay_inci(i) = decay_inci
383 pblast%PBLAST_TAB(il)%decay_refl(i) = decay_refl
384
385 ELSE
386 !nothing to compute underground
387 z=zero
388 cos_theta = zero
389 p_inci = zero
390 p_refl = zero
391 t_a = ep20
392 dt_0 = ep20
393 decay_inci = zero
394 decay_refl = zero
395 ENDIF
396
397 dtmin_loc = min(dtmin_loc,dt_0/ndt)
398 iz_update = 1 !update done
399 iloadp(06,nl) = iz_update
400 is_updated=.true.
401
402 else! => IZ_UPDATE=1
403
404 !use wave parameters from Starter
405 z = zero ! ZG not used here since all wave characteristics are stored.
406 cos_theta = pblast%PBLAST_TAB(il)%cos_theta(i)
407 p_inci = pblast%PBLAST_TAB(il)%P_inci(i)
408 p_refl = pblast%PBLAST_TAB(il)%P_refl(i)
409 t_a = pblast%PBLAST_TAB(il)%ta(i)
410 dt_0 = pblast%PBLAST_TAB(il)%t0(i)
411 decay_inci = pblast%PBLAST_TAB(il)%decay_inci(i)
412 decay_refl = pblast%PBLAST_TAB(il)%decay_refl(i)
413 dtmin_loc = pblast%PBLAST_TAB(il)%DTMIN
414
415 ENDIF !IF(IZ_UPDATE == 2)
416
417 !Coefficients for wave superimposition
418 !PressureLoad = Reflected_Pressure * cos2X + IncidentPressure * (1 + cos2X -2 cosX)
419 IF(cos_theta<=zero)THEN
420 !Surface not facing the point of explosion
421 alpha_refl = zero
422 alpha_inci = one
423 ELSE
424 alpha_refl = cos_theta**2 ! cos**2 a
425 alpha_inci = one + cos_theta - two * alpha_refl ! 1 + cos a -2 cos**2 a
426 ENDIF
427
428 !Building pressure waves from Friedlander model. (Modified model can bu introduced later if needed)
429 wave_inci = zero
430 wave_refl = zero
431 IF(tt_star>=t_a)THEN
432 wave_inci = p_inci*(one-(tt_star-t_a)/dt_0)*exp(-decay_inci*(tt_star-t_a)/dt_0)
433 wave_refl = p_refl*(one-(tt_star-t_a)/dt_0)*exp(-decay_refl*(tt_star-t_a)/dt_0)
434 ELSE
435 wave_inci = zero
436 wave_refl = zero
437 ENDIF
438 p = alpha_refl * wave_refl + alpha_inci * wave_inci
439 p = max(p,pmin)
440 pblast%PBLAST_TAB(il)%PRES(i) = p
441
442 !!Expand Pressure load to nodes
443 ! FF is nodal force which applied on each node N1,N2,N3, and also N4 if relevant
444 ! FF = FF_elem / NPT = Pload.S.n / NPT where n is the unitary normal vector
445 ! NX,NY,NZ = 2S.n (in all cases:quadrangles & triangles)
446 surf_patch = half*sqrt(nx*nx+ny*ny+nz*nz) / npt
447 ff(1) = -p * half*nx / npt ! -P*S/NPT . nx
448 ff(2) = -p * half*ny / npt ! -P*S/NPT . ny
449 ff(3) = -p * half*nz / npt ! -P*S/NPT . nz
450 !storing force for one node of the current face (for assembly below)
451 pblast%PBLAST_TAB(il)%FX(i) = ff(1)
452 pblast%PBLAST_TAB(il)%FY(i) = ff(2)
453 pblast%PBLAST_TAB(il)%FZ(i) = ff(3)
454 pblast%PBLAST_TAB(il)%SURF_PATCH(i) = surf_patch
455
456 !External Force work
457 ! on a given node : DW = <F,V>*dt
458 ! for this current 4-node or 3-node face : DW = sum( <F_k,V_k>*dt k=1,NPT) where F_k=Fel/NPT
459 wfext_loc=wfext_loc+dt1*(ff(1)*sum(v(1,nn(1:nint(npt)))) +ff(2)*sum(v(2,nn(1:nint(npt)))) +ff(3)*sum(v(3,nn(1:nint(npt)))))
460
461C----- /TH/SURF -------
462 IF(th_surf%LOADP_FLAG > 0 ) THEN
463 nsegpl = nsegpl + 1
464 area = surf_patch*npt
465 DO ns=th_surf%LOADP_KSEGS(nsegpl) +1,th_surf%LOADP_KSEGS(nsegpl+1)
466 ksurf = th_surf%LOADP_SEGS(ns)
467 th_surf%channels(4,ksurf)= th_surf%channels(4,ksurf) + area*p ! mean pressure
468 th_surf%channels(5,ksurf)= th_surf%channels(5,ksurf) + area ! surface where pressure is applied
469 ENDDO
470 ENDIF
471
472
473 enddo!next I
474!$OMP END DO
475
476 IF(imodel == 2 .AND. nwarn > 0)THEN
477 msgout1=''
478 WRITE(msgout1,fmt='(I0,A)') nwarn,
479 . ' SEGMENT(S) HAS EXCESSIVE POSITIVE IMPULSE REGARDING THE PEAK PRESSURE AND POSITIVE DURATION.'
480 msgout2=''
481 msgout2='A TRIANGULAR WAVEFORM WILL BE USED INSTEAD TO MAXIMIZE THE IMPULSE. DEFINING A PMIN VALUE IS STRONGLY RECOMMENDED'
482 write(iout , fmt='(A,I10,/A,/A)') "Updated parameters for /LOAD/PBLAST id=", id, msgout1, msgout2
483 write(istdo, fmt='(A,I10,/A,/A)') "Updated parameters for /LOAD/PBLAST id=", id, msgout1, msgout2
484 ENDIF
485
486 CALL my_barrier
487 IF(is_updated)THEN
488#include "lockon.inc"
489 !---arrival time
490 zeta = fac(07,nl)
491 fac(07,nl) = min(ta_inf_loc, fac(07,nl)) !smp min value
492 !---time step
493 dtmin_loc = (one+em06)*(fac(07,nl) - tt) ! go directly to trigger time
494 dtmin_loc=max(pblast%PBLAST_TAB(il)%DTMIN, dtmin_loc)
495 !---no update on next cycle
496 iz_update = 1 !update done
497 iloadp(06,nl) = iz_update
498#include "lockoff.inc"
499!$OMP SINGLE
500 write(*,fmt='(A,I10,A,E16.8,A,E16.8)') "Updated parameters for /LOAD/PBLAST id=",
501 . id,' previous first arrival time :',zeta,
502 . ' is now updated to :',fac(07,nl)
503!$OMP END SINGLE
504 ENDIF
505
506 !-------------------------------------------------------------------!
507 ! FORCE ASSEMBLY !
508 ! /PARITH/OFF : F directly added in A(1:3,1:NUMNOD). !
509 ! /PARITH/ON : F added FSKY & and automatically treated later !
510 !-------------------------------------------------------------------!
511 ! SPMD/SMP Parith/OFF
512 IF(iparit==0) THEN
513!$OMP SINGLE
514 DO i = 1,isiz_seg
515 n1=lloadp(iloadp(4,nl)+4*(i-1))
516 n2=lloadp(iloadp(4,nl)+4*(i-1)+1)
517 n3=lloadp(iloadp(4,nl)+4*(i-1)+2)
518 n4=lloadp(iloadp(4,nl)+4*(i-1)+3)
519 a(1,n1)=a(1,n1)+pblast%PBLAST_TAB(il)%FX(i)
520 a(2,n1)=a(2,n1)+pblast%PBLAST_TAB(il)%FY(i)
521 a(3,n1)=a(3,n1)+pblast%PBLAST_TAB(il)%FZ(i)
522 a(1,n2)=a(1,n2)+pblast%PBLAST_TAB(il)%FX(i)
523 a(2,n2)=a(2,n2)+pblast%PBLAST_TAB(il)%FY(i)
524 a(3,n2)=a(3,n2)+pblast%PBLAST_TAB(il)%FZ(i)
525 a(1,n3)=a(1,n3)+pblast%PBLAST_TAB(il)%FX(i)
526 a(2,n3)=a(2,n3)+pblast%PBLAST_TAB(il)%FY(i)
527 a(3,n3)=a(3,n3)+pblast%PBLAST_TAB(il)%FZ(i)
528 IF(pblast%PBLAST_TAB(il)%NPt(i) == four)THEN
529 a(1,n4)=a(1,n4)+pblast%PBLAST_TAB(il)%FX(i)
530 a(2,n4)=a(2,n4)+pblast%PBLAST_TAB(il)%FY(i)
531 a(3,n4)=a(3,n4)+pblast%PBLAST_TAB(il)%FZ(i)
532 ENDIF
533 ENDDO
534!$OMP END SINGLE
535 ELSE
536!$OMP DO SCHEDULE(GUIDED,MVSIZ)
537 DO i = 1,isiz_seg
538 iad =iadc(iloadp(4,nl)+4*(i-1))
539 fsky(1,iad) =pblast%PBLAST_TAB(il)%FX(i)
540 fsky(2,iad) =pblast%PBLAST_TAB(il)%FY(i)
541 fsky(3,iad) =pblast%PBLAST_TAB(il)%FZ(i)
542 iad =iadc(iloadp(4,nl)+4*(i-1)+1)
543 fsky(1,iad) =pblast%PBLAST_TAB(il)%FX(i)
544 fsky(2,iad) =pblast%PBLAST_TAB(il)%FY(i)
545 fsky(3,iad) =pblast%PBLAST_TAB(il)%FZ(i)
546 iad =iadc(iloadp(4,nl)+4*(i-1)+2)
547 fsky(1,iad) =pblast%PBLAST_TAB(il)%FX(i)
548 fsky(2,iad) =pblast%PBLAST_TAB(il)%FY(i)
549 fsky(3,iad) =pblast%PBLAST_TAB(il)%FZ(i)
550 IF(pblast%PBLAST_TAB(il)%NPt(i) == four)THEN
551 iad =iadc(iloadp(4,nl)+4*(i-1)+3)
552 fsky(1,iad) =pblast%PBLAST_TAB(il)%FX(i)
553 fsky(2,iad) =pblast%PBLAST_TAB(il)%FY(i)
554 fsky(3,iad) =pblast%PBLAST_TAB(il)%FZ(i)
555 ENDIF
556 ENDDO
557!$OMP END DO
558 ENDIF !IPARIT
559
560
561 !-------------------------------------------!
562 ! ANIMATION FILE /ANIM/VECT/FEXT !
563 ! H3D FILE /H3D/NODA/FEXT !
564 !-------------------------------------------!
565!$OMP SINGLE
566 IF(ianim_or_h3d > 0) THEN
567 DO i = 1,isiz_seg
568 n1=pblast%PBLAST_TAB(il)%N(1,i)
569 n2=pblast%PBLAST_TAB(il)%N(2,i)
570 n3=pblast%PBLAST_TAB(il)%N(3,i)
571 n4=pblast%PBLAST_TAB(il)%N(4,i)
572 fext(1,n1) = fext(1,n1)+pblast%PBLAST_TAB(il)%FX(i)
573 fext(2,n1) = fext(2,n1)+pblast%PBLAST_TAB(il)%FY(i)
574 fext(3,n1) = fext(3,n1)+pblast%PBLAST_TAB(il)%FZ(i)
575 fext(1,n2) = fext(1,n2)+pblast%PBLAST_TAB(il)%FX(i)
576 fext(2,n2) = fext(2,n2)+pblast%PBLAST_TAB(il)%FY(i)
577 fext(3,n2) = fext(3,n2)+pblast%PBLAST_TAB(il)%FZ(i)
578 fext(1,n3) = fext(1,n3)+pblast%PBLAST_TAB(il)%FX(i)
579 fext(2,n3) = fext(2,n3)+pblast%PBLAST_TAB(il)%FY(i)
580 fext(3,n3) = fext(3,n3)+pblast%PBLAST_TAB(il)%FZ(i)
581 IF(pblast%PBLAST_TAB(il)%NPt(i)==four)THEN
582 fext(1,n4) = fext(1,n4)+pblast%PBLAST_TAB(il)%FX(i)
583 fext(2,n4) = fext(2,n4)+pblast%PBLAST_TAB(il)%FY(i)
584 fext(3,n4) = fext(3,n4)+pblast%PBLAST_TAB(il)%FZ(i)
585 ENDIF
586 ENDDO
587 ENDIF
588 IF(th_has_noda_pext > 0 .OR. anim_has_noda_pext > 0 .OR. h3d_has_noda_pext > 0) THEN
589 DO i = 1,isiz_seg
590 n1 = pblast%PBLAST_TAB(il)%N(1,i)
591 n2 = pblast%PBLAST_TAB(il)%N(2,i)
592 n3 = pblast%PBLAST_TAB(il)%N(3,i)
593 n4 = pblast%PBLAST_TAB(il)%N(4,i)
594 surf_patch = pblast%PBLAST_TAB(il)%SURF_PATCH(i)
595 noda_surf(n1) = noda_surf(n1) + surf_patch
596 noda_surf(n2) = noda_surf(n2) + surf_patch
597 noda_surf(n3) = noda_surf(n3) + surf_patch
598 p = pblast%PBLAST_TAB(il)%PRES(i) * surf_patch
599 noda_pext(n1) = noda_pext(n1) + p
600 noda_pext(n2) = noda_pext(n2) + p
601 noda_pext(n3) = noda_pext(n3) + p
602 IF(pblast%PBLAST_TAB(il)%NPT(i) == four)THEN
603 noda_surf(n4) = noda_surf(n4) + surf_patch
604 noda_pext(n4) = noda_pext(n4) + p
605 ENDIF
606 ENDDO
607 ENDIF
608!$OMP END SINGLE
609
610 RETURN
611
612C-----------------------------------------------
613 IF (ierr1/=0) THEN
614 WRITE(iout,*)' ** ERROR IN MEMORY ALLOCATION - PBLAST LOADING'
615 WRITE(istdo,*)' ** ERROR IN MEMORY ALLOCATION - PBLAST LOADING'
616 CALL arret(2)
617 END IF
618C-----------------------------------------------
619
620 END SUBROUTINE
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
id
initmumps id
Definition multiplerhs_example.m:4
groupdef_mod
Definition groupdef_mod.F:662
h3d_inc_mod
Definition h3d_inc_mod.F:42
h3d_mod
Definition h3d_mod.F:308
names_and_titles_mod
Definition names_and_titles_mod.F:997
names_and_titles_mod::ncharline
integer, parameter ncharline
Definition names_and_titles_mod.F:1001
th_mod
Definition th_mod.F:42
th_mod::th_has_noda_pext
integer th_has_noda_pext
Definition th_mod.F:121
th_surf_mod
OPTION /TH/SURF outputs of Pressure and Area needed Tabs.
Definition th_surf_mod.F:60
pblast_2
subroutine pblast_2(pblast, iloadp, fac, a, v, x, iadc, fsky, lloadp, fext, noda_surf, noda_pext, itab, h3d_data, nl, dtmin_loc, wfext_loc, th_surf, nsegpl)
Definition pblast_2.F:46
arret
subroutine arret(nn)
Definition arret.F:87
my_barrier
subroutine my_barrier
Definition machine.F:31
h3d_mod::h3d_database
Definition h3d_mod.F:435
th_surf_mod::th_surf_
Definition th_surf_mod.F:75