OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pblast_1.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_1 ../engine/source/loads/pblast/pblast_1.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!|| omp_get_thread_num ../engine/source/engine/openmp_stub.F90
31!|| pblast_parameters__free_air ../common_source/modules/loads/pblast_mod.f90
32!||--- uses -----------------------------------------------------
33!|| groupdef_mod ../common_source/modules/groupdef_mod.F
34!|| h3d_inc_mod ../engine/share/modules/h3d_inc_mod.f
35!|| h3d_mod ../engine/share/modules/h3d_mod.F
36!|| names_and_titles_mod ../common_source/modules/names_and_titles_mod.F
37!|| output_mod ../common_source/modules/output/output_mod.F90
38!|| pblast_mod ../common_source/modules/loads/pblast_mod.F90
39!|| th_mod ../engine/share/modules/th_mod.F
40!|| th_surf_mod ../common_source/modules/interfaces/th_surf_mod.F
41!||====================================================================
42 SUBROUTINE pblast_1(PBLAST ,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(friedlander_params_) :: FRIEDLANDER_PARAMS
94 TYPE(pblast_) :: 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,
99 . ID, ITA_SHIFT,IMODEL,NN(4),NS,KSURF,NDT
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
102 my_real :: zx,zy,zz,norm,nx,ny,nz,area,
103 . cos_theta, alpha_inci, alpha_refl, p_inci, p_refl_,p_inci_, p_refl,z,
104 . i_inci,i_refl,i_inci_,i_refl_, dt_0, t_a,dt_0_,
105 . wave_refl,wave_inci, w13,
106 . dnorm, xdet,ydet,zdet,tdet,wtnt,pmin,t_stop,dx,dy,dz,p,
107 . fac_m_bb, fac_l_bb, fac_t_bb, fac_p_bb, fac_i_bb, ta_first, tt_star, z1_,
108 . decay_inci,decay_refl,zeta,
109 . cst_255_div_ln_z1_on_zn, log10_, npt, ff(3), ta_inf_loc,
110 . surf_patch
111
112 INTEGER, EXTERNAL :: OMP_GET_THREAD_NUM
113
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 (FREE AIR 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_ = half
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 t_stop = fac(13,nl)
190 is = iloadp(02,nl)
191 iz_update = iloadp(06,nl)
192 abac_id = iloadp(07,nl)
193 id = iloadp(08,nl) !user_id
194 ita_shift = iloadp(09,nl)
195 ndt = iloadp(10,nl)
196 imodel = iloadp(11,nl)
197 isiz_seg = iloadp(01,nl)/4
198 ierr1 = 0
199 w13 = (wtnt*fac_m_bb)**third ! '*FAC_M' g->work unit '/FAC_M' : WORK_UNIT -> g
200 z = zero
201 ta_inf_loc = ep20
202
203 is_decay_to_be_computed = .false.
204 IF(imodel == 2)is_decay_to_be_computed=.true.
205 nwarn = 0
206
207 !---------------------------------------------
208 ! LOOP ON SEGMENTS (4N or 3N)
209 !---------------------------------------------
210
211!$OMP DO SCHEDULE(GUIDED,MVSIZ)
212 DO i = 1,isiz_seg
213 n1=lloadp(iloadp(4,nl)+4*(i-1))
214 n2=lloadp(iloadp(4,nl)+4*(i-1)+1)
215 n3=lloadp(iloadp(4,nl)+4*(i-1)+2)
216 n4=lloadp(iloadp(4,nl)+4*(i-1)+3)
217 nn(1)=n1
218 nn(2)=n2
219 nn(3)=n3
220 nn(4)=n4
221
222 IF(n4==0 .OR. n3==n4 )THEN
223 !3-NODE-SEGMENT
224 pblast%PBLAST_TAB(il)%NPt(i) = three
225 npt = three
226 !Segment Centroid
227 zx = x(1,n1)+x(1,n2)+x(1,n3)
228 zy = x(2,n1)+x(2,n2)+x(2,n3)
229 zz = x(3,n1)+x(3,n2)+x(3,n3)
230 zx = zx*third
231 zy = zy*third
232 zz = zz*third
233 !Normal vector : (NX,NY,NZ) = 2*S*n where |n|=1.0
234 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))
235 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))
236 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))
237 !NORM = 2*S
238 norm = sqrt(nx*nx+ny*ny+nz*nz)
239 ELSE
240 !4-NODE-SEGMENT
241 pblast%PBLAST_TAB(il)%NPt(i) = four
242 npt = four
243 !Face Centroid
244 zx = x(1,n1)+x(1,n2)+x(1,n3)+x(1,n4)
245 zy = x(2,n1)+x(2,n2)+x(2,n3)+x(2,n4)
246 zz = x(3,n1)+x(3,n2)+x(3,n3)+x(3,n4)
247 zx = zx*fourth
248 zy = zy*fourth
249 zz = zz*fourth
250 !Normal vector (NX,NY,NZ) = 2*S*n where |n|=1.0
251 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))
252 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))
253 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))
254 !NORM = 2*S
255 norm = sqrt(nx*nx+ny*ny+nz*nz)
256 ENDIF
257
258 pblast%PBLAST_TAB(il)%N(1,i) = n1
259 pblast%PBLAST_TAB(il)%N(2,i) = n2
260 pblast%PBLAST_TAB(il)%N(3,i) = n3
261 pblast%PBLAST_TAB(il)%N(4,i) = n4
262
263 !--------------------------------------------!
264 ! Update Wave parameters !
265 ! (experimental) !
266 ! If requested. Otherwise use Starter param. !
267 ! Default : do not update:use Starter param. !
268 !--------------------------------------------!
269 IF(iz_update == 2)THEN
270
271 dtmin_loc = ep20
272
273 !Dist
274 dx = (xdet - zx)*fac_l_bb ! => working unit to cm
275 dy = (ydet - zy)*fac_l_bb ! => ... to cm
276 dz = (zdet - zz)*fac_l_bb ! => ... to cm
277 dnorm = sqrt(dx*dx+dy*dy+dz*dz)
278
279 !Angle from detonation point
280 cos_theta = dx*nx + dy*ny + dz*nz
281 cos_theta = cos_theta/max(em20,norm*dnorm)
282
283 !scaled distance
284 z = dnorm / w13 !in abac unit ID g,cm,mus
285 !display warning message if out of bounds
286 ! precondition for SUBROUTINE PBLAST_PARAMETERS__FREE_AIR
287 IF(z <= 0.5 .AND. pblast%PBLAST_TAB(il)%TAGMSG(i) == 0 )THEN
288 write(iout, fmt='(A,I0,A)') "Warning : /LOAD/PBLAST id=",id," R/W**(1/3) < 0.5 cm/g**(1/3)"
289 write(istdo,fmt='(A,I0,A)') "Warning : /LOAD/PBLAST id=",id," R/W**(1/3) < 0.5 cm/g**(1/3)"
290 write(iout, fmt='(A)') " Radial Distance R too close to the charge"
291 write(istdo,fmt='(A)') " Radial Distance R too close to the charge"
292
293 if (n4 == 0 .OR. n3 == n4)then
294 write(iout, fmt='(a,3i11)') " -> Segment nodes : ", itab(n1),itab(n2),itab(n3)
295 write(istdo,fmt='(a,3i11)') " -> Segment nodes : ", itab(n1),itab(n2),itab(n3)
296 else
297 write(iout, fmt='(a,4i11)') " -> Segment nodes : ", itab(n1),itab(n2),itab(n3),itab(n4)
298 write(istdo,fmt='(a,4i11)') " -> Segment nodes : ", itab(n1),itab(n2),itab(n3),itab(n4)
299 endif
300 PBLAST%PBLAST_TAB(IL)%TAGMSG(I) = 1
301
302 ELSEIF(Z > 400. AND. PBLAST%PBLAST_TAB(IL)%TAGMSG(I) == 0)THEN
303 write(iout, fmt='(a,i0,a)') "Warning : /LOAD/PBLAST id=",ID," R/W**(1/3) > 400. cm/g**(1/3)"
304 write(istdo,fmt='(a,i0,a)') "Warning : /LOAD/PBLAST id=",ID," R/W**(1/3) > 400. cm/g**(1/3)"
305 WRITE(iout, FMT='(a)') " Radial Distance R too far from the charge"
306 WRITE(istdo,FMT='(a)') " Radial Distance R too far from the charge"
307.OR. if (N4 == 0 N3 == N4)then
308 write(iout, fmt='(a,3i11)') " -> Segment nodes : ", itab(n1),itab(n2),itab(n3)
309 write(istdo,fmt='(a,3i11)') " -> Segment nodes : ", itab(n1),itab(n2),itab(n3)
310 else
311 write(iout, fmt='(a,4i11)') " -> Segment nodes : ", itab(n1),itab(n2),itab(n3),itab(n4)
312 write(istdo,fmt='(a,4i11)') " -> Segment nodes : ", itab(n1),itab(n2),itab(n3),itab(n4)
313 endif
314 PBLAST%PBLAST_TAB(IL)%TAGMSG(I) = 1
315 ENDIF
316
317 !------------------------------------------------------------------!
318 CALL PBLAST_PARAMETERS__FREE_AIR(PBLAST,Z, W13, TDET,
319 + FAC_P_bb, FAC_I_bb, FAC_T_bb,
320 + IS_DECAY_TO_BE_COMPUTED,
321 + FRIEDLANDER_PARAMS,NWARN)
322
323 P_inci = FRIEDLANDER_PARAMS%P_inci
324 P_inci_ = FRIEDLANDER_PARAMS%P_inci_
325 P_refl = FRIEDLANDER_PARAMS%P_refl
326 P_refl_ = FRIEDLANDER_PARAMS%P_refl_
327 I_inci = FRIEDLANDER_PARAMS%I_inci
328 I_inci_ = FRIEDLANDER_PARAMS%I_inci_
329 I_refl = FRIEDLANDER_PARAMS%I_refl
330 I_refl_ = FRIEDLANDER_PARAMS%I_refl_
331 T_A = FRIEDLANDER_PARAMS%T_A
332 DT_0 = FRIEDLANDER_PARAMS%DT_0
333 DT_0_ = FRIEDLANDER_PARAMS%DT_0_
334 decay_inci = FRIEDLANDER_PARAMS%decay_inci
335 decay_refl = FRIEDLANDER_PARAMS%decay_refl
336 !------------------------------------------------------------------!
337
338 TA_INF_LOC = MIN(TA_INF_LOC, T_A)
339
340 !update wave parameters
341 PBLAST%PBLAST_TAB(IL)%cos_theta(I) = cos_theta
342 PBLAST%PBLAST_TAB(IL)%P_inci(I) = P_inci
343 PBLAST%PBLAST_TAB(IL)%P_refl(I) = P_refl
344 PBLAST%PBLAST_TAB(IL)%ta(I) = T_A
345 PBLAST%PBLAST_TAB(IL)%t0(I) = DT_0
346 PBLAST%PBLAST_TAB(IL)%decay_inci(I) = decay_inci
347 PBLAST%PBLAST_TAB(IL)%decay_refl(I) = decay_refl
348
349 DTMIN_LOC = MIN(DTMIN_LOC,DT_0/NDT)
350 IS_UPDATED = .TRUE.
351
352 ELSE! => IZ_UPDATE=1
353
354 !use wave parameters from Starter
355 Z=ZERO
356 cos_theta = PBLAST%PBLAST_TAB(IL)%cos_theta(I)
357 P_inci = PBLAST%PBLAST_TAB(IL)%P_inci(I)
358 P_refl = PBLAST%PBLAST_TAB(IL)%P_refl(I)
359 T_A = PBLAST%PBLAST_TAB(IL)%ta(I)
360 DT_0 = PBLAST%PBLAST_TAB(IL)%t0(I)
361 decay_inci = PBLAST%PBLAST_TAB(IL)%decay_inci(I)
362 decay_refl = PBLAST%PBLAST_TAB(IL)%decay_refl(I)
363 DTMIN_LOC = PBLAST%PBLAST_TAB(IL)%DTMIN
364
365 ENDIF !IF(IZ_UPDATE == 2)
366
367 !Coefficients for wave superimposition
368 !PressureLoad = Reflected_Pressure * cos2X + IncidentPressure * (1 + cos2X -2 cosX)
369 IF(cos_theta<=ZERO)THEN
370 !Surface not facing the point of explosion
371 alpha_refl = ZERO
372 alpha_inci = ONE
373 ELSE
374 alpha_refl = cos_theta**2 ! cos**2 a
375 alpha_inci = ONE + cos_theta - TWO * alpha_refl ! 1 + cos a -2 cos**2 a
376 ENDIF
377
378 !Building pressure waves from Friedlander model. (Modified model can bu introduced later if needed)
379 WAVE_INCI = ZERO
380 WAVE_REFL = ZERO
381 IF(TT_STAR>=T_A)THEN
382 WAVE_INCI = P_inci*(ONE-(TT_STAR-T_A)/DT_0)*exp(-DECAY_inci*(TT_STAR-T_A)/DT_0)
383 WAVE_REFL = P_refl*(ONE-(TT_STAR-T_A)/DT_0)*exp(-DECAY_refl*(TT_STAR-T_A)/DT_0)
384 ELSE
385 WAVE_INCI = ZERO
386 WAVE_REFL = ZERO
387 ENDIF
388 P = alpha_refl * WAVE_REFL + alpha_inci * WAVE_INCI
389 P = MAX(P,PMIN)
390 PBLAST%PBLAST_TAB(IL)%PRES(I) = P
391
392 !!Expand Pressure load to nodes
393 ! FF is nodal force which applied on each node N1,N2,N3, and also N4 if relevant
394 ! FF = FF_elem / NPT = Pload.S.n / NPT where n is the unitary normal vector
395 ! NX,NY,NZ = 2S.n (in all cases:quadrangles & triangles)
396 SURF_PATCH = HALF*SQRT(NX*NX+NY*NY+NZ*NZ) / NPT
397 FF(1) = -P * HALF*NX / NPT ! -P*S/NPT . nx
398 FF(2) = -P * HALF*NY / NPT ! -P*S/NPT . ny
399 FF(3) = -P * HALF*NZ / NPT ! -P*S/NPT . nz
400 !storing force for one node of the current face (for assembly below)
401 PBLAST%PBLAST_TAB(IL)%FX(I) = FF(1)
402 PBLAST%PBLAST_TAB(IL)%FY(I) = FF(2)
403 PBLAST%PBLAST_TAB(IL)%FZ(I) = FF(3)
404 PBLAST%PBLAST_TAB(IL)%SURF_PATCH(I) = SURF_PATCH
405
406 !External Force work
407 ! on a given node : DW = <F,V>*dt
408 ! for this current 4-node or 3-node face : DW = sum( <F_k,V_k>*dt k=1,NPT) where F_k=Fel/NPT
409 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)))))
410
411C----- /TH/SURF -------
412 IF(TH_SURF%LOADP_FLAG > 0 ) THEN
413 NSEGPL = NSEGPL + 1
414 AREA = SURF_PATCH * NPT
415 DO NS=TH_SURF%LOADP_KSEGS(NSEGPL) +1,TH_SURF%LOADP_KSEGS(NSEGPL+1)
416 KSURF = TH_SURF%LOADP_SEGS(NS)
417 th_surf%channels(4,KSURF)= th_surf%channels(4,KSURF) + AREA*P ! mean pressure
418 th_surf%channels(5,KSURF)= th_surf%channels(5,KSURF) + AREA ! suface where pressure is applied
419 ENDDO
420 ENDIF
421
422 ENDDO!next I
423!$OMP END DO
424
425.AND. IF(IMODEL == 2 NWARN > 0)THEN
426 MSGOUT1=''
427 WRITE(MSGOUT1,FMT='(i0,a)') NWARN,
428 . ' segment(s) has excessive positive impulse regarding the peak pressure and positive duration.'
429 MSGOUT2=''
430 MSGOUT2='a triangular waveform will be used instead to maximize the impulse. defining a pmin VALUE is strongly recommended'
431 write(IOUT , FMT='(a,i10,/a,/a)') "Updated parameters for /LOAD/PBLAST id=", ID, MSGOUT1, MSGOUT2
432 write(ISTDO, FMT='(a,i10,/a,/a)') "Updated parameters for /LOAD/PBLAST id=", ID, MSGOUT1, MSGOUT2
433 ENDIF
434
435 CALL MY_BARRIER
436 IF(IS_UPDATED)THEN
437#include "lockon.inc"
438 !---arrival time
439 ZETA = FAC(07,NL)
440 FAC(07,NL) = MIN(TA_INF_LOC, FAC(07,NL)) !smp min value
441 !---time step
442 DTMIN_LOC = (ONE+EM06)*(FAC(07,NL) - TT) ! go directly to trigger time
443 DTMIN_LOC=MAX(PBLAST%PBLAST_TAB(IL)%DTMIN, DTMIN_LOC)
444 !---no update on next cycle
445 IZ_UPDATE = 1 !update done
446 ILOADP(06,NL) = IZ_UPDATE
447#include "lockoff.inc"
448!$OMP SINGLE
449 write(IOUT,FMT='(a,i10,a,e16.8,a,e16.8)') "Updated parameters for /LOAD/PBLAST id=",
450 . ID,' previous first arrival time :',ZETA,
451 . ' is now updated to :',FAC(07,NL)
452 write(ISTDO,FMT='(a,i10,a,e16.8,a,e16.8)') "Updated parameters for /LOAD/PBLAST id=",
453 . ID,' previous first arrival time :',ZETA,
454 . ' is now updated to :',FAC(07,NL)
455!$OMP END SINGLE
456 ENDIF
457
458 !-------------------------------------------------------------------!
459 ! FORCE ASSEMBLY !
460 ! /PARITH/OFF : F directly added in A(1:3,1:NUMNOD). !
461 ! /PARITH/ON : F added FSKY & and automatically treated later !
462 !-------------------------------------------------------------------!
463 ! SPMD/SMP Parith/OFF
464 IF(IPARIT==0) THEN
465!$OMP SINGLE
466 DO I = 1,ISIZ_SEG
467 N1=LLOADP(ILOADP(4,NL)+4*(I-1))
468 N2=LLOADP(ILOADP(4,NL)+4*(I-1)+1)
469 N3=LLOADP(ILOADP(4,NL)+4*(I-1)+2)
470 N4=LLOADP(ILOADP(4,NL)+4*(I-1)+3)
471 A(1,N1)=A(1,N1)+PBLAST%PBLAST_TAB(IL)%FX(I)
472 A(2,N1)=A(2,N1)+PBLAST%PBLAST_TAB(IL)%FY(I)
473 A(3,N1)=A(3,N1)+PBLAST%PBLAST_TAB(IL)%FZ(I)
474 A(1,N2)=A(1,N2)+PBLAST%PBLAST_TAB(IL)%FX(I)
475 A(2,N2)=A(2,N2)+PBLAST%PBLAST_TAB(IL)%FY(I)
476 A(3,N2)=A(3,N2)+PBLAST%PBLAST_TAB(IL)%FZ(I)
477 A(1,N3)=A(1,N3)+PBLAST%PBLAST_TAB(IL)%FX(I)
478 A(2,N3)=A(2,N3)+PBLAST%PBLAST_TAB(IL)%FY(I)
479 A(3,N3)=A(3,N3)+PBLAST%PBLAST_TAB(IL)%FZ(I)
480 IF(PBLAST%PBLAST_TAB(IL)%NPt(I) == FOUR)THEN
481 A(1,N4)=A(1,N4)+PBLAST%PBLAST_TAB(IL)%FX(I)
482 A(2,N4)=A(2,N4)+PBLAST%PBLAST_TAB(IL)%FY(I)
483 A(3,N4)=A(3,N4)+PBLAST%PBLAST_TAB(IL)%FZ(I)
484 ENDIF
485 ENDDO
486!$OMP END SINGLE
487 ELSE
488!$OMP DO SCHEDULE(GUIDED,MVSIZ)
489 DO I = 1,ISIZ_SEG
490 IAD =IADC(ILOADP(4,NL)+4*(I-1))
491 FSKY(1,IAD) =PBLAST%PBLAST_TAB(IL)%FX(I)
492 FSKY(2,IAD) =PBLAST%PBLAST_TAB(IL)%FY(I)
493 FSKY(3,IAD) =PBLAST%PBLAST_TAB(IL)%FZ(I)
494 IAD =IADC(ILOADP(4,NL)+4*(I-1)+1)
495 FSKY(1,IAD) =PBLAST%PBLAST_TAB(IL)%FX(I)
496 FSKY(2,IAD) =PBLAST%PBLAST_TAB(IL)%FY(I)
497 FSKY(3,IAD) =PBLAST%PBLAST_TAB(IL)%FZ(I)
498 IAD =IADC(ILOADP(4,NL)+4*(I-1)+2)
499 FSKY(1,IAD) =PBLAST%PBLAST_TAB(IL)%FX(I)
500 FSKY(2,IAD) =PBLAST%PBLAST_TAB(IL)%FY(I)
501 FSKY(3,IAD) =PBLAST%PBLAST_TAB(IL)%FZ(I)
502 IF(PBLAST%PBLAST_TAB(IL)%NPt(I) == FOUR)THEN
503 IAD =IADC(ILOADP(4,NL)+4*(I-1)+3)
504 FSKY(1,IAD) =PBLAST%PBLAST_TAB(IL)%FX(I)
505 FSKY(2,IAD) =PBLAST%PBLAST_TAB(IL)%FY(I)
506 FSKY(3,IAD) =PBLAST%PBLAST_TAB(IL)%FZ(I)
507 ENDIF
508 ENDDO
509!$OMP END DO
510 ENDIF !IPARIT
511
512
513 !-------------------------------------------!
514 ! ANIMATION FILE /ANIM/VECT/FEXT !
515 ! H3D FILE /H3D/NODA/FEXT !
516 !-------------------------------------------!
517!$OMP SINGLE
518 IF(IANIM_OR_H3D > 0) THEN
519 DO I = 1,ISIZ_SEG
520 N1=PBLAST%PBLAST_TAB(IL)%N(1,I)
521 N2=PBLAST%PBLAST_TAB(IL)%N(2,I)
522 N3=PBLAST%PBLAST_TAB(IL)%N(3,I)
523 N4=PBLAST%PBLAST_TAB(IL)%N(4,I)
524 FEXT(1,N1) = FEXT(1,N1)+PBLAST%PBLAST_TAB(IL)%FX(I)
525 FEXT(2,N1) = FEXT(2,N1)+PBLAST%PBLAST_TAB(IL)%FY(I)
526 FEXT(3,N1) = FEXT(3,N1)+PBLAST%PBLAST_TAB(IL)%FZ(I)
527 FEXT(1,N2) = FEXT(1,N2)+PBLAST%PBLAST_TAB(IL)%FX(I)
528 FEXT(2,N2) = FEXT(2,N2)+PBLAST%PBLAST_TAB(IL)%FY(I)
529 FEXT(3,N2) = FEXT(3,N2)+PBLAST%PBLAST_TAB(IL)%FZ(I)
530 FEXT(1,N3) = FEXT(1,N3)+PBLAST%PBLAST_TAB(IL)%FX(I)
531 FEXT(2,N3) = FEXT(2,N3)+PBLAST%PBLAST_TAB(IL)%FY(I)
532 FEXT(3,N3) = FEXT(3,N3)+PBLAST%PBLAST_TAB(IL)%FZ(I)
533 IF(PBLAST%PBLAST_TAB(IL)%NPt(I)==FOUR)THEN
534 FEXT(1,N4) = FEXT(1,N4)+PBLAST%PBLAST_TAB(IL)%FX(I)
535 FEXT(2,N4) = FEXT(2,N4)+PBLAST%PBLAST_TAB(IL)%FY(I)
536 FEXT(3,N4) = FEXT(3,N4)+PBLAST%PBLAST_TAB(IL)%FZ(I)
537 ENDIF
538 ENDDO
539 ENDIF
540.OR..OR. IF(TH_HAS_NODA_PEXT > 0 ANIM_HAS_NODA_PEXT > 0 H3D_HAS_NODA_PEXT > 0) THEN
541 DO I = 1,ISIZ_SEG
542 N1 = PBLAST%PBLAST_TAB(IL)%N(1,I)
543 N2 = PBLAST%PBLAST_TAB(IL)%N(2,I)
544 N3 = PBLAST%PBLAST_TAB(IL)%N(3,I)
545 N4 = PBLAST%PBLAST_TAB(IL)%N(4,I)
546 SURF_PATCH = PBLAST%PBLAST_TAB(IL)%SURF_PATCH(I)
547 NODA_SURF(N1) = NODA_SURF(N1) + SURF_PATCH
548 NODA_SURF(N2) = NODA_SURF(N2) + SURF_PATCH
549 NODA_SURF(N3) = NODA_SURF(N3) + SURF_PATCH
550 P = PBLAST%PBLAST_TAB(IL)%PRES(I) * SURF_PATCH
551 NODA_PEXT(N1) = NODA_PEXT(N1) + P
552 NODA_PEXT(N2) = NODA_PEXT(N2) + P
553 NODA_PEXT(N3) = NODA_PEXT(N3) + P
554 IF(PBLAST%PBLAST_TAB(IL)%NPT(I) == FOUR)THEN
555 NODA_SURF(N4) = NODA_SURF(N4) + SURF_PATCH
556 NODA_PEXT(N4) = NODA_PEXT(N4) + P
557 ENDIF
558 ENDDO
559 ENDIF
560!$OMP END SINGLE
561
562 RETURN
563
564C-----------------------------------------------
565 IF (IERR1/=0) THEN
566 WRITE(IOUT,*)' ** error in memory allocation - pblast loading'
567 WRITE(ISTDO,*)' ** error in memory allocation - pblast loading'
568 CALL ARRET(2)
569 END IF
570C-----------------------------------------------
571 END SUBROUTINE PBLAST_1
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
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_1
subroutine pblast_1(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_1.F:46
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