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