OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_setoff_npg_c.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!|| fail_setoff_npg_c ../engine/source/materials/fail/fail_setoff_npg_c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!|| usermat_shell ../engine/source/materials/mat_share/usermat_shell.F
28!||--- uses -----------------------------------------------------
29!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
30!|| failwave_mod ../common_source/modules/failwave_mod.F
31!|| mat_elem_mod ../common_source/modules/mat_elem/mat_elem_mod.F90
32!|| stack_mod ../engine/share/modules/stack_mod.F
33!||====================================================================
34 SUBROUTINE fail_setoff_npg_c(
35 . ELBUF_STR,MAT_ELEM ,GEO ,PID ,
36 . NGL ,NEL ,IR ,IS ,
37 . NLAY ,NPTTOT ,THK_LY ,THKLY ,
38 . OFF ,NPG ,STACK ,ISUBSTACK,
39 . IGTYP ,FAILWAVE ,FWAVE_EL ,NLAY_MAX ,
40 . LAYNPT_MAX,NUMGEO ,IPG ,NUMSTACK ,
41 . IGEO ,PRINT_FAIL)
42C-----------------------------------------------
43C M o d u l e s
44C-----------------------------------------------
45 USE mat_elem_mod
46 USE stack_mod
47 USE failwave_mod
48 USE stack_mod
49 USE elbufdef_mod
50C-----------------------------------------------
51C I m p l i c i t T y p e s
52C-----------------------------------------------
53#include "implicit_f.inc"
54C-----------------------------------------------
55C C o m m o n B l o c k s
56C-----------------------------------------------
57#include "param_c.inc"
58#include "com08_c.inc"
59#include "units_c.inc"
60#include "comlock.inc"
61C-----------------------------------------------
62C D u m m y A r g u m e n t s
63C-----------------------------------------------
64 TYPE(elbuf_struct_), INTENT(INOUT), TARGET :: ELBUF_STR
65 my_real, DIMENSION(NPROPG,NUMGEO), INTENT(IN) :: GEO
66 INTEGER, DIMENSION(NPROPGI,NUMGEO),INTENT(IN) :: IGEO
67 INTEGER, INTENT(IN) :: PID,NEL,IR,IS,NLAY,NPTTOT,NPG,IGTYP,
68 . ISUBSTACK,NLAY_MAX,LAYNPT_MAX,NUMGEO,
69 . IPG,NUMSTACK
70 INTEGER, DIMENSION(NEL), INTENT(IN) :: NGL
71 my_real, DIMENSION(NEL,NLAY_MAX*LAYNPT_MAX), INTENT(IN) :: thk_ly
72 my_real, DIMENSION(NPTTOT*NEL), INTENT(IN) :: thkly
73 my_real, DIMENSION(NEL), INTENT(INOUT) :: off
74 TYPE (STACK_PLY), INTENT(IN) :: STACK
75 TYPE (FAILWAVE_STR_), INTENT(IN), TARGET :: FAILWAVE
76 INTEGER, DIMENSION(NEL), INTENT(INOUT) :: FWAVE_EL
77 TYPE (MAT_ELEM_) ,INTENT(INOUT) :: MAT_ELEM
78 LOGICAL, DIMENSION(NEL), INTENT(INOUT) :: PRINT_FAIL
79C-----------------------------------------------
80C L o c a l V a r i a b l e s
81C-----------------------------------------------
82 INTEGER I,II,IEL,IPOS,IL,IFL,IP,IPT,IG,JPG,NPTR,NPTS,NPTT,
83 . countpg,nindxly,ipt_all,nfail,ipweight,ipthkly,
84 . ips,ipr,id_ply,imat
85 my_real :: p_thickg,fail_exp,thfact,norm,dfail,npfail
86 my_real, DIMENSION(NLAY,100) :: pthkf
87 INTEGER, DIMENSION(NEL) :: INDXLY,FAIL_NUM
88 INTEGER, DIMENSION(:), POINTER :: FOFF,LAY_OFF,OFFPG
89 my_real, DIMENSION(NLAY) :: weight,p_thkly
90 TYPE(l_bufel_) ,POINTER :: LBUF
91 CHARACTER(LEN=NCHARTITLE) :: FAIL_NAME
92c-----------------------------------------------------------------------
93c NPTT NUMBER OF INTEGRATION POINTS IN CURRENT LAYER
94c NPTTF NUMBER OF FAILED INTEGRATION POINTS IN THE LAYER
95c NPTTOT NUMBER OF INTEGRATION POINTS IN ALL LAYERS (TOTAL)
96c OFFPG(NEL,NPG) failure flag of PG in each layer 1=alive ,0=dead
97c THK_LY Ratio of layer thickness / element thickness
98c THK Total element thickness
99C=======================================================================
100c
101 imat = 1
102 !=================================================================
103 ! RECOVER PARAMETERS AND INITIALIZATION
104 !=================================================================
105 p_thickg = geo(42,pid)
106 fail_exp = geo(43,pid)
107 IF (igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
108 ipthkly = 1 + 4*nlay ! Address of PTHKLY in STACK%GEO table
109 ipweight = ipthkly + nlay ! Address of WEIGHT in STACK%GEO table
110 ELSE
111 ipthkly = 700 ! Address of PTHKLY in GEO table
112 ipweight = 900 ! Address of WEIGHT in GEO table
113 ENDIF
114 nptr = elbuf_str%NPTR
115 npts = elbuf_str%NPTS
116 jpg = (ipg-1)*nel
117c
118 DO il=1,nlay
119 nfail = elbuf_str%BUFLY(il)%NFAIL
120 imat = elbuf_str%BUFLY(il)%IMAT
121 DO ifl = 1,nfail
122 pthkf(il,ifl) = mat_elem%MAT_PARAM(imat)%FAIL(ifl)%PTHK
123 END DO
124 END DO
125 !=================================================================
126 ! 1 LAYER PROPERTIES - TYPE1/TYPE9
127 !=================================================================
128 IF (nlay == 1) THEN
129c
130 ! Only 1 layer with several integration points
131 il = 1
132c
133 ! material internal identifier number
134 imat = elbuf_str%BUFLY(il)%IMAT
135c
136 ! Check PTHICKFAIL coming from failure criteria
137 nfail = elbuf_str%BUFLY(il)%NFAIL
138 DO ifl = 1,nfail
139 ! -> Percentage of broken thickness
140 IF (pthkf(il,ifl) > zero) THEN
141 pthkf(il,ifl) = min(pthkf(il,ifl),abs(p_thickg))
142 pthkf(il,ifl) = max(min(pthkf(il,ifl),one-em06),em06)
143 ! -> Ratio of broken integration points
144 ELSEIF (pthkf(il,ifl) < zero) THEN
145 pthkf(il,ifl) = max(pthkf(il,ifl),-abs(p_thickg))
146 pthkf(il,ifl) = min(max(pthkf(il,ifl),-one+em6),-em06)
147 ! -> If not defined in the failure criterion, the value of the property is used by default
148 ELSE
149 pthkf(il,ifl) = p_thickg
150 ENDIF
151 ENDDO ! |-> IFL
152c
153 ! Check in-plane Gauss point failure
154 nptt = elbuf_str%BUFLY(il)%NPTT
155 offpg => elbuf_str%BUFLY(il)%OFFPG(jpg+1:jpg+nel)
156 DO iel=1,nel
157 IF (off(iel) == one) THEN
158 DO ifl = 1,nfail
159 thfact = zero
160 npfail = zero
161 DO ipt=1,nptt
162 foff => elbuf_str%BUFLY(il)%FAIL(ir,is,ipt)%FLOC(ifl)%OFF
163 ! Computation of broken fraction of thickness /
164 ! ratio of intg. points
165 IF (foff(iel) < 1) THEN
166 ipos = (ipt-1)*nel + iel
167 thfact = thfact + thkly(ipos)
168 npfail = npfail + one/nptt
169 ENDIF
170 ! Comparison with critical value PTHICKFAIL
171 IF (((thfact >= pthkf(il,ifl)).AND.(pthkf(il,ifl) > zero)).OR.
172 . ((npfail >= abs(pthkf(il,ifl))).AND.(pthkf(il,ifl) < zero))) THEN
173 offpg(iel) = 0
174 fail_num(iel) = ifl
175 ENDIF
176 ENDDO ! |-> IPT
177 ENDDO ! |-> IFL
178 ENDIF
179 ENDDO ! |-> IEL
180c
181 ! Check element failure (when IPG = NPG = 4)
182 IF (ipg == npg) THEN
183 DO iel=1,nel
184 IF (off(iel) == one) THEN
185 countpg = 0
186 DO ig=1,npg
187 jpg = (ig-1)*nel
188 countpg = countpg + elbuf_str%BUFLY(il)%OFFPG(jpg+iel)
189 ENDDO ! |-> IG
190 IF (countpg == 0) THEN
191 off(iel) = four_over_5
192 print_fail(iel) = .false.
193 fail_name = mat_elem%MAT_PARAM(imat)%FAIL(fail_num(iel))%KEYWORD
194#include "lockon.inc"
195 WRITE(iout, 1000) trim(fail_name),ngl(iel)
196 WRITE(istdo,1100) trim(fail_name),ngl(iel),tt
197#include "lockoff.inc"
198 IF (failwave%WAVE_MOD > 0) fwave_el(iel) = -1
199 ENDIF
200 ENDIF
201 ENDDO ! |-> IEL
202 ENDIF
203c
204 !=================================================================
205 ! MULTI LAYER PROPERTIES / 1 INTG POINT - TYPE10/11/16/17/51/52
206 !=================================================================
207 ELSEIF (nlay == npttot) THEN
208c
209 ! Only one integration points in each layer
210 ipt = 1
211c
212 ! Check in-plane Gauss point failure for each layers
213 DO il=1,nlay
214 nfail = elbuf_str%BUFLY(il)%NFAIL
215 lay_off => elbuf_str%BUFLY(il)%OFF
216 offpg =>elbuf_str%BUFLY(il)%OFFPG(jpg+1:jpg+nel)
217 imat = elbuf_str%BUFLY(il)%IMAT
218 DO iel=1,nel
219 IF (off(iel) == one .AND. lay_off(iel) == 1) THEN
220 DO ifl = 1,nfail
221 foff => elbuf_str%BUFLY(il)%FAIL(ir,is,ipt)%FLOC(ifl)%OFF
222 IF (foff(iel) < 1) THEN
223 offpg(iel) = 0
224 fail_num(iel) = ifl
225 ENDIF
226 ENDDO ! |-> IFL
227 ENDIF
228 ENDDO ! |-> IEL
229 ! Check layer failure only if IPG = NPG = 4
230 IF (ipg == npg) THEN
231 nindxly = 0
232 lay_off => elbuf_str%BUFLY(il)%OFF
233 DO iel = 1,nel
234 IF (off(iel) == one) THEN
235 IF (lay_off(iel) == 1) THEN
236 countpg = 0
237 DO ig=1,npg
238 jpg = (ig-1)*nel
239 countpg = countpg + elbuf_str%BUFLY(il)%OFFPG(jpg+iel)
240 ENDDO ! |-> NPG
241 IF (countpg == 0) THEN
242 nindxly = nindxly + 1
243 indxly(nindxly) = iel
244 lay_off(iel) = 0
245 ENDIF
246 ENDIF
247 ENDIF
248 ENDDO ! |-> IEL
249 ! Printing out layer/ply failure message
250 IF (nindxly > 0) THEN
251 ! -> Print out ply failure
252 IF (igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
253 IF (igtyp == 17 .OR. igtyp == 51) THEN
254 id_ply = igeo(1,stack%IGEO(2+il,isubstack))
255 ELSE
256 id_ply = ply_info(1,stack%IGEO(2+il,isubstack)-numstack)
257 ENDIF
258 DO i = 1,nindxly
259 fail_name = mat_elem%MAT_PARAM(imat)%FAIL(fail_num(indxly(i)))%KEYWORD
260#include "lockon.inc"
261 WRITE(iout, 3000) trim(fail_name),id_ply,ngl(indxly(i))
262 WRITE(istdo,3100) trim(fail_name),id_ply,ngl(indxly(i)),tt
263#include "lockoff.inc"
264 ENDDO ! |-> I
265 ! -> Print out layer failure
266 ELSE
267 DO i = 1,nindxly
268 fail_name = mat_elem%MAT_PARAM(imat)%FAIL(fail_num(indxly(i)))%KEYWORD
269#include "lockon.inc"
270 WRITE(iout, 2000) trim(fail_name),il,ngl(indxly(i))
271 WRITE(istdo,2100) trim(fail_name),il,ngl(indxly(i)),tt
272#include "lockoff.inc"
273 ENDDO ! |-> I
274 ENDIF
275 ENDIF
276 ENDIF
277 ENDDO
278c
279 ! Check element failure
280 DO iel=1,nel
281 IF (off(iel) == one) THEN
282 thfact = zero
283 norm = zero
284 npfail = zero
285 DO il=1,nlay
286 IF (igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
287 weight(il) = stack%GEO(ipweight+ il,isubstack)
288 ELSE
289 weight(il) = geo(ipweight + il,pid)
290 ENDIF
291 lay_off => elbuf_str%BUFLY(il)%OFF
292 ipos = (il-1)*nel + iel
293 dfail = thkly(ipos)*weight(il)
294 norm = norm + dfail
295 IF (off(iel) == one .AND. lay_off(iel) == 0) THEN
296 thfact = thfact + thkly(ipos)*weight(il)
297 npfail = npfail + one/nlay
298 ENDIF
299 ENDDO ! |-> IL
300 IF (((thfact >= p_thickg*norm).AND.(p_thickg > zero)).OR.
301 . ((npfail >= abs(p_thickg)).AND.(p_thickg < zero))) THEN
302 off(iel) = four_over_5
303 IF (failwave%WAVE_MOD > 0) fwave_el(iel) = -1
304 ENDIF
305 ENDIF
306 ENDDO ! |-> IEL
307c
308 !=======================================================================
309 ! MULTI LAYER PROPERTIES / SEVERAL INTG. POINTS - TYPE51/TYPE52
310 !=======================================================================
311 ELSE
312c
313 ! Several integration points in each layer
314 ipt_all = 0
315c
316 ! Check PTHICKFAIL coming from failure criteria
317 DO il = 1,nlay
318 nfail = elbuf_str%BUFLY(il)%NFAIL
319 p_thkly(il) = stack%GEO(ipthkly + il,isubstack)
320 DO ifl = 1,nfail
321 ! -> Percentage of broken thickness
322 IF (pthkf(il,ifl) > zero) THEN
323 pthkf(il,ifl) = min(pthkf(il,ifl),abs(p_thkly(il)))
324 pthkf(il,ifl) = max(min(pthkf(il,ifl),one-em06),em06)
325 ! -> Ratio of broken integration points
326 ELSEIF (pthkf(il,ifl) < zero) THEN
327 pthkf(il,ifl) = max(pthkf(il,ifl),-abs(p_thkly(il)))
328 pthkf(il,ifl) = min(max(pthkf(il,ifl),-one+em6),-em06)
329 ! -> If not defined in the failure criterion, the value of the property is used by default
330 ELSE
331 pthkf(il,ifl) = p_thkly(il)
332 ENDIF
333 ENDDO ! |-> IFL
334 ENDDO ! |-> IL
335c
336 ! Check in-plane Gauss point failure for each layers
337 DO il=1,nlay
338 nptt = elbuf_str%BUFLY(il)%NPTT
339 nfail = elbuf_str%BUFLY(il)%NFAIL
340 lay_off => elbuf_str%BUFLY(il)%OFF
341 offpg =>elbuf_str%BUFLY(il)%OFFPG(jpg+1:jpg+nel)
342 weight(il) = stack%GEO(ipweight + il,isubstack)
343 imat = elbuf_str%BUFLY(il)%IMAT
344 DO iel=1,nel
345 IF (off(iel) == one .AND. lay_off(iel) == 1) THEN
346 DO ifl = 1,nfail
347 thfact = zero
348 npfail = zero
349 DO ipt = 1,nptt
350 foff => elbuf_str%BUFLY(il)%FAIL(ir,is,ipt)%FLOC(ifl)%OFF
351 IF (foff(iel) < one) THEN
352 ip = ipt_all + ipt
353 ipos = (ip-1)*nel + iel
354 thfact = thfact + thkly(ipos)/thk_ly(iel,il)
355 npfail = npfail + one/nptt
356 ENDIF
357 IF (((thfact >= pthkf(il,ifl)).AND.(pthkf(il,ifl)>zero)).OR.
358 . ((thfact >= abs(pthkf(il,ifl))).AND.(pthkf(il,ifl)<zero))) THEN
359 offpg(iel) = 0
360 fail_num(iel) = ifl
361 ENDIF
362 ENDDO ! |-> IPT
363 ENDDO ! |-> IFL
364 ENDIF
365 ENDDO ! |-> IEL
366 ipt_all = ipt_all + nptt
367 ENDDO ! |-> IL
368 ! Check layer failure only if IPG = NPG = 4
369 IF (ipg == npg) THEN
370 DO iel=1,nel
371 IF (off(iel) == one) THEN
372 DO il=1,nlay
373 nfail = elbuf_str%BUFLY(il)%NFAIL
374 lay_off => elbuf_str%BUFLY(il)%OFF
375 nindxly = 0
376 IF (lay_off(iel) == 1) THEN
377 countpg = 0
378 DO ig=1,npg
379 jpg = (ig-1)*nel
380 countpg = countpg + elbuf_str%BUFLY(il)%OFFPG(jpg+iel)
381 ENDDO ! |-> IG
382 IF (countpg == 0) THEN ! all Gauss pts failed
383 nindxly = nindxly + 1
384 indxly(nindxly) = iel
385 lay_off(iel) = 0
386 nptt = elbuf_str%BUFLY(il)%NPTT
387 DO ifl = 1,nfail
388 DO ipr=1,nptr
389 DO ips=1,npts
390 DO ipt=1,nptt
391 foff => elbuf_str%BUFLY(il)%FAIL(ipr,ips,ipt)%FLOC(ifl)%OFF
392 foff(iel) = 0
393 ENDDO ! |-> IPT
394 ENDDO ! |-> IPS
395 ENDDO ! |-> IPR
396 ENDDO ! |-> IFL
397 ENDIF
398 ENDIF
399 ! Printing out ply failure message
400 IF (nindxly > 0) THEN
401 IF (igtyp == 51) THEN
402 id_ply = igeo(1,stack%IGEO(2+il,isubstack))
403 ELSE
404 id_ply = ply_info(1,stack%IGEO(2+il,isubstack)-numstack)
405 ENDIF
406 DO i = 1,nindxly
407 fail_name = mat_elem%MAT_PARAM(imat)%FAIL(fail_num(indxly(i)))%KEYWORD
408#include "lockon.inc"
409 WRITE(iout, 3000) trim(fail_name),id_ply,ngl(indxly(i))
410 WRITE(istdo,3100) trim(fail_name),id_ply,ngl(indxly(i)),tt
411#include "lockoff.inc"
412 ENDDO ! |-> I
413 ENDIF
414 ENDDO ! |-> IL
415 ENDIF
416 ENDDO ! |-> IEL
417c
418 ! Check element failure
419 DO iel=1,nel
420 IF (off(iel) == one) THEN
421 thfact = zero
422 norm = zero
423 npfail = zero
424 DO il=1,nlay
425 weight(il) = stack%GEO(ipweight+ il,isubstack)
426 lay_off => elbuf_str%BUFLY(il)%OFF
427 dfail = (thk_ly(iel,il)*weight(il))**fail_exp
428 norm = norm + dfail
429 IF (lay_off(iel) == 0) THEN
430 thfact = thfact + dfail
431 npfail = npfail + one/nlay
432 ENDIF
433 ENDDO ! |-> IL
434 thfact = thfact**(one/fail_exp)
435 norm = norm**(one/fail_exp)
436 IF (((thfact >= p_thickg*norm).AND.(p_thickg > zero)).OR.
437 . ((thfact >= abs(p_thickg)).AND.(p_thickg < zero))) THEN
438 off(iel) = four_over_5
439 IF (failwave%WAVE_MOD > 0) fwave_el(iel) = -1
440 ENDIF
441 ENDIF
442 ENDDO ! |-> IEL
443 ENDIF
444c
445 ENDIF ! IGTYP PROPERTY TYPE
446 !=======================================================================
447c
448 !=======================================================================
449 ! PRINTING OUT FORMATS FOR LAYERS/PLYS FAILURE
450 !=======================================================================
451 1000 FORMAT(1x,'-- RUPTURE (',a,') OF SHELL ELEMENT NUMBER ',i10)
452 1100 FORMAT(1x,'-- RUPTURE (',a,') OF SHELL ELEMENT :',i10,' AT TIME :',g11.4)
453 2000 FORMAT(1x,'-- FAILURE (',a,') OF LAYER',i3, ' ,SHELL ELEMENT NUMBER ',i10)
454 2100 FORMAT(1x,'-- FAILURE (',a,') OF LAYER',i3, ' ,SHELL ELEMENT NUMBER ',i10,
455 . 1x,'AT TIME :',g11.4)
456 3000 FORMAT(1x,'-- FAILURE (',a,') OF PLY ID ',i10, ' ,SHELL ELEMENT NUMBER ',i10)
457 3100 FORMAT(1x,'-- FAILURE (',a,') OF PLY ID ',i10, ' ,SHELL ELEMENT NUMBER ',i10,
458 . 1x,'AT TIME :',g11.4)
459 END
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
fail_setoff_npg_c
subroutine fail_setoff_npg_c(elbuf_str, mat_elem, geo, pid, ngl, nel, ir, is, nlay, npttot, thk_ly, thkly, off, npg, stack, isubstack, igtyp, failwave, fwave_el, nlay_max, laynpt_max, numgeo, ipg, numstack, igeo, print_fail)
Definition fail_setoff_npg_c.F:42
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
failwave_mod
Definition failwave_mod.F:61
stack_mod
Definition stack_mod.F:72
stack_mod::ply_info
integer, dimension(:,:), allocatable ply_info
Definition stack_mod.F:133