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