OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
crklen4n_adv.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!|| crklen4n_adv ../engine/source/elements/xfem/crklen4n_adv.F
25!||--- called by ------------------------------------------------------
26!|| cforc3 ../engine/source/elements/shell/coque/cforc3.F
27!|| czforc3 ../engine/source/elements/shell/coquez/czforc3.F
28!||--- calls -----------------------------------------------------
29!|| arret ../engine/source/system/arret.F
30!||--- uses -----------------------------------------------------
31!|| crackxfem_mod ../engine/share/modules/crackxfem_mod.F
32!|| element_mod ../common_source/modules/elements/element_mod.F90
33!||====================================================================
34 SUBROUTINE crklen4n_adv(
35 . NEL ,NFT ,ILAY ,NLAY ,IXC ,
36 . CRKLEN ,ELCRKINI ,IEL_CRK ,DIR1 ,DIR2 ,
37 . NODEDGE ,CRKEDGE ,XEDGE4N ,NGL ,XL2 ,
38 . XL3 ,XL4 ,YL2 ,YL3 ,YL4 ,
39 . ALDT )
40C-----------------------------------------------
41C crack advancement, shells 4N
42C-----------------------------------------------
43C M o d u l e s
44C-----------------------------------------------
45 USE crackxfem_mod
46 use element_mod , only : nixc
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51C-----------------------------------------------
52C C o m m o n B l o c k s
53C-----------------------------------------------
54#include "units_c.inc"
55C-----------------------------------------------
56C D u m m y A r g u m e n t s
57C-----------------------------------------------
58 INTEGER NEL,NFT,ILAY,NLAY
59 INTEGER IXC(NIXC,*),NGL(NEL),IEL_CRK(*),ELCRKINI(NLAY,NEL),
60 . NODEDGE(2,*),XEDGE4N(4,*)
61 my_real DIR1(NLAY,NEL),DIR2(NLAY,NEL),CRKLEN(NEL),ALDT(NEL),
62 . XL2(NEL),YL2(NEL),XL3(NEL),YL3(NEL),XL4(NEL),YL4(NEL)
63 TYPE (XFEM_EDGE_) , DIMENSION(*) :: CRKEDGE
64C-----------------------------------------------
65C L o c a l V a r i a b l e s
66C-----------------------------------------------
67 INTEGER I,K,R,IR,P1,P2,NEWCRK,IED,OK,ELCRK,
68 . fac,iedge,icut,nod1,nod2
69c
70 INTEGER JCT(NEL),EDGEL(4,NEL),TIP(NEL),ECUT(2,NEL),dd(4),d(8),KPERM(8)
71c
72 my_real
73 . xin(2,nel),yin(2,nel),len(4,nel),
74 . xxl(4,nel),yyl(4,nel),beta0(4,nel)
75 my_real
76 . xint,yint,cross,acd,bcd,dlx,dly,
77 . m12,mm,xint0,yint0,dir11,dir22,
78 . beta,bmin,bmax
79C----------
80 DATA d/1,2,2,3,4,3,1,4/
81 DATA dd/2,3,4,1/
82 DATA kperm/1,2,3,4,1,2,3,4/
83 parameter(bmin = 0.01, bmax = 0.99)
84C=======================================================================
85 newcrk = 0
86 DO i=1,nel
87 jct(i) = 0
88 IF (elcrkini(ilay,i) == 5) THEN ! crack propagation
89 newcrk = newcrk + 1
90 jct(newcrk) = i
91 elcrkini(ilay,i) = 2 ! reset for propagation
92 ELSEIF (elcrkini(ilay,i) == -5) THEN ! crack initialization
93 crklen(i) = aldt(i)
94 elcrkini(ilay,i) = 0 ! reset for initialization
95 ENDIF
96 ENDDO
97 IF (newcrk == 0) RETURN
98C---
99 DO ir=1,newcrk
100 i = jct(ir)
101 tip(i) = 0
102 ecut(1:2,i) = 0
103 edgel(1:4,i) = 0
104 beta0(1:4,i) = zero
105 xin(1,i) = zero ! first inters point in local skew
106 yin(1,i) = zero
107 xin(2,i) = zero ! second inters point in local skew
108 yin(2,i) = zero
109c
110 xxl(1,i) = zero
111 yyl(1,i) = zero
112 xxl(2,i) = xl2(i)
113 yyl(2,i) = yl2(i)
114 xxl(3,i) = xl3(i)
115 yyl(3,i) = yl3(i)
116 xxl(4,i) = xl4(i)
117 yyl(4,i) = yl4(i)
118c
119 len(1,i) = xl2(i)*xl2(i) + yl2(i)*yl2(i)
120 len(2,i) = (xl3(i)-xl2(i))*(xl3(i)-xl2(i))+
121 . (yl3(i)-yl2(i))*(yl3(i)-yl2(i))
122 len(3,i) = (xl4(i)-xl3(i))*(xl4(i)-xl3(i))+
123 . (yl4(i)-yl3(i))*(yl4(i)-yl3(i))
124 len(4,i) = xl4(i)*xl4(i) + yl4(i)*yl4(i)
125 ENDDO
126C------------------------------------------------
127c First intersection (already cut edge)
128C------------------------------------------------
129 DO ir=1,newcrk ! loop over elems (layers) with advancing cracks
130 i = jct(ir)
131 elcrk = iel_crk(i+nft) ! DYE SYS XFEM
132 ok = 0
133 icut = 0
134 ied = 0
135 DO k=1,4 ! edges
136 iedge = xedge4n(k,elcrk)
137 IF (iedge > 0) THEN
138 icut = crkedge(ilay)%ICUTEDGE(iedge)
139 IF (icut == 1) THEN
140 nod1 = nodedge(1,iedge) ! noeud std
141 nod2 = nodedge(2,iedge)
142 IF (nod1 == ixc(k+1,i) .and. nod2 == ixc(dd(k)+1,i)) THEN
143 p1 = k
144 p2 = dd(k)
145 ELSE IF (nod2 == ixc(k+1,i) .and. nod1 == ixc(dd(k)+1,i)) THEN
146 p1 = dd(k)
147 p2 = k
148 ENDIF
149 ok = 1
150 ied = k
151 ecut(1,i)= k
152 EXIT
153 ENDIF ! IF (ICUT == 1) THEN
154 ENDIF
155 ENDDO ! DO K=1,4
156C---
157 IF (ok == 1) THEN ! edge found
158 beta = crkedge(ilay)%RATIO(iedge)
159 xin(1,i) = xxl(p1,i) + beta*(xxl(p2,i) - xxl(p1,i))
160 yin(1,i) = yyl(p1,i) + beta*(yyl(p2,i) - yyl(p1,i))
161c
162 edgel(ied,i) = 1 ! Local: First Edge Cup
163 iedge = xedge4n(ied,elcrk) ! N edge element sys xfem
164 tip(i) = crkedge(ilay)%EDGETIP(1,iedge) ! 1 or 2 , start or end of crack
165 ELSE
166 WRITE(iout,*) 'ERROR IN ADVANCING CRACK --- CHECK CRACK TIP'
167 CALL arret(2)
168 ENDIF
169C---
170 END DO ! DO IR=1,NEWCRK
171C--------------------------------------------------
172c Search for second intersection (new cut edge)
173C--------------------------------------------------
174 DO ir=1,newcrk
175 i = jct(ir)
176 elcrk = iel_crk(i+nft)
177 r = ecut(1,i)
178 xint0 = xin(1,i)
179 yint0 = yin(1,i)
180 dir11 =-dir2(ilay,i)
181 dir22 = dir1(ilay,i)
182C---
183 IF (dir11 == zero) THEN
184 DO k=1,3
185 r = kperm(ecut(1,i) + k)
186 iedge = xedge4n(r,elcrk)
187 nod1 = nodedge(1,iedge)
188 nod2 = nodedge(2,iedge)
189 IF (nod1 == ixc(r+1,i) .and. nod2 == ixc(dd(r)+1,i))THEN
190 p1 = r
191 p2 = dd(r)
192 ELSE IF (nod2 == ixc(r+1,i).and.nod1 == ixc(dd(r)+1,i))THEN
193 p1 = dd(r)
194 p2 = r
195 ENDIF
196 dlx = xxl(p2,i) - xxl(p1,i)
197 IF (dlx /= zero) THEN
198 dly = yyl(p2,i) - yyl(p1,i)
199 m12 = dly / dlx
200 xint = xint0
201 yint = yyl(p1,i) + m12*(xint-xxl(p1,i))
202 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero .and.
203 . (yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
204 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
205 beta = sqrt(cross / len(r,i))
206 IF (beta > bmax .OR. beta < bmin) THEN
207 beta = max(beta, bmin)
208 beta = min(beta, bmax)
209 yint = yyl(p1,i) + beta*(yyl(p2,i)-yyl(p1,i))
210 ENDIF
211C
212 ecut(2,i) = r
213 xin(2,i) = xint
214 yin(2,i) = yint
215 edgel(r,i) = 2
216 beta0(r,i) = beta
217 EXIT
218 ENDIF
219 ENDIF
220 ENDDO
221c
222 ELSEIF (dir22 == zero) THEN
223 DO k=1,3
224 r = kperm(ecut(1,i) + k)
225 iedge = xedge4n(r,elcrk)
226 nod1 = nodedge(1,iedge)
227 nod2 = nodedge(2,iedge)
228 IF (nod1 == ixc(r+1,i) .and. nod2 == ixc(dd(r)+1,i)) THEN
229 p1 = r
230 p2 = dd(r)
231 ELSE IF (nod2 == ixc(r+1,i).and.nod1 == ixc(dd(r)+1,i)) THEN
232 p1 = dd(r)
233 p2 = r
234 ENDIF
235 dly = yyl(p2,i) - yyl(p1,i)
236 IF (dly /= zero) THEN
237 dlx = xxl(p2,i) - xxl(p1,i)
238 m12 = dlx / dly
239 yint = yint0
240 xint = xxl(p1,i) + m12*(yint-yyl(p1,i))
241 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero .and.
242 . (yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
243 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
244 beta = sqrt(cross / len(r,i))
245 IF (beta > bmax .OR. beta < bmin) THEN
246 beta = max(beta, bmin)
247 beta = min(beta, bmax)
248 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
249 ENDIF
250C
251 ecut(2,i) = r
252 xin(2,i) = xint
253 yin(2,i) = yint
254 edgel(r,i) = 2
255 beta0(r,i) = beta
256 EXIT
257 ENDIF
258 ENDIF
259 ENDDO
260c
261 ELSEIF (dir11 /= zero .and. dir22 /= zero) THEN
262 DO k=1,3
263 r = kperm(ecut(1,i) + k)
264 iedge = xedge4n(r,elcrk)
265 nod1 = nodedge(1,iedge)
266 nod2 = nodedge(2,iedge)
267 IF (nod1 == ixc(r+1,i) .and. nod2 == ixc(dd(r)+1,i)) THEN
268 p1 = r
269 p2 = dd(r)
270 ELSE IF (nod2 == ixc(r+1,i).and.nod1 == ixc(dd(r)+1,i)) THEN
271 p1 = dd(r)
272 p2 = r
273 ENDIF
274C
275 dlx = xxl(p2,i) - xxl(p1,i)
276 dly = yyl(p2,i) - yyl(p1,i)
277 mm = dir22/dir11
278 IF (dlx == zero) THEN
279 xint = xxl(p1,i)
280 yint = yint0 + mm*(xint-xint0)
281 IF ((yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
282 cross = (yyl(p1,i) - yint)**2
283 beta = sqrt(cross / len(r,i))
284 IF (beta > bmax .OR. beta < bmin) THEN
285 beta = max(beta, bmin)
286 beta = min(beta, bmax)
287 yint = yyl(p1,i) + beta*(yyl(p2,i)-yyl(p1,i))
288 ENDIF
289 ecut(2,i) = r
290 xin(2,i) = xint
291 yin(2,i) = yint
292 edgel(r,i) = 2
293 beta0(r,i) = beta
294 EXIT
295 ENDIF
296 ELSEIF (dly == zero) THEN
297 yint = yyl(p1,i)
298 xint = xint0 + (yint0-yyl(p1,i)) / mm
299 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero) THEN
300 cross = (xxl(p1,i) - xint)**2
301 beta = sqrt(cross / len(r,i))
302 IF (beta > bmax .OR. beta < bmin) THEN
303 beta = max(beta, bmin)
304 beta = min(beta, bmax)
305 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
306 ENDIF
307 ecut(2,i) = r
308 xin(2,i) = xint
309 yin(2,i) = yint
310 edgel(r,i) = 2
311 beta0(r,i) = beta
312 EXIT
313 ENDIF
314 ELSE
315 m12 = dly / dlx
316 IF (mm /= m12) THEN
317 xint = (yint0-yyl(p1,i) + m12*xxl(p1,i) - mm*xint0)/(m12-mm)
318 yint = yint0 + mm*(xint-xint0)
319 acd = (yint-yyl(p1,i))*(xint0 - xxl(p1,i))
320 . - (xint-xxl(p1,i))*(yint0 - yyl(p1,i))
321 bcd = (yint-yyl(p2,i))*(xint0 - xxl(p2,i))
322 . - (xint-xxl(p2,i))*(yint0 - yyl(p2,i))
323 IF (acd*bcd <= zero) THEN
324 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
325 beta = sqrt(cross / len(r,i))
326 IF (beta > bmax .OR. beta < bmin) THEN
327 beta = max(beta, bmin)
328 beta = min(beta, bmax)
329 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
330 yint = yyl(p1,i) + beta*(yyl(p2,i)-yyl(p1,i))
331 ENDIF
332 ecut(2,i) = r
333 xin(2,i) = xint
334 yin(2,i) = yint
335 edgel(r,i) = 2
336 beta0(r,i) = beta
337 EXIT
338 ENDIF
339 ENDIF
340 ENDIF
341 ENDDO
342 ENDIF
343 ENDDO
344C-----------------------------------------------------------------------
345C check for getting both intersections
346C-----------------------------------------------------------------------
347 DO ir=1,newcrk
348 i = jct(ir)
349 fac = 0
350 DO r=1,4
351 IF (edgel(r,i)==1 .or. edgel(r,i)==2) fac=fac+1
352 ENDDO
353 IF (fac /= 2) THEN
354 WRITE(iout,*) 'ERROR IN ADVANCING CRACK. NO CUT EDGES'
355 CALL arret(2)
356 ENDIF
357 crklen(i) = sqrt((xin(2,i) - xin(1,i))**2 + (yin(2,i) - yin(1,i))**2)
358 ENDDO
359c-----------
360 RETURN
361 END
crklen4n_adv
subroutine crklen4n_adv(nel, nft, ilay, nlay, ixc, crklen, elcrkini, iel_crk, dir1, dir2, nodedge, crkedge, xedge4n, ngl, xl2, xl3, xl4, yl2, yl3, yl4, aldt)
Definition crklen4n_adv.F:40
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
crackxfem_mod
Definition crackxfem_mod.F:121
arret
subroutine arret(nn)
Definition arret.F:86