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