OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
crklayer4n_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!|| crklayer4n_adv ../engine/source/elements/xfem/crklayer4n_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!|| lsint4 ../engine/source/elements/xfem/crklayer4n_adv.F
31!||--- uses -----------------------------------------------------
32!|| crackxfem_mod ../engine/share/modules/crackxfem_mod.F
33!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
34!|| element_mod ../common_source/modules/elements/element_mod.F90
35!||====================================================================
36 SUBROUTINE crklayer4n_adv(
37 . XFEM_STR ,NEL ,NFT ,IXC ,ELCUTC ,
38 . ILAY ,NLAY ,IEL_CRK ,INOD_CRK ,
39 . IADC_CRK ,NODENR ,ELCRKINI ,DIR1 ,DIR2 ,
40 . NODEDGE ,CRKNODIAD ,KNOD2ELC ,CRKEDGE ,A_I ,
41 . XL2 ,XL3 ,XL4 ,YL2 ,YL3 ,
42 . YL4 ,XEDGE4N ,NGL )
43C-----------------------------------------------
44C crack advancement, shells 4N
45C-----------------------------------------------
46C M o d u l e s
47C-----------------------------------------------
48 USE elbufdef_mod
49 USE crackxfem_mod
50 use element_mod , only : nixc
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 "units_c.inc"
59#include "com_xfem1.inc"
60C-----------------------------------------------
61C D u m m y A r g u m e n t s
62C-----------------------------------------------
63 INTEGER NEL,NFT,IXC(NIXC,*),ILAY,NLAY,NGL(NEL),IEL_CRK(*),
64 . INOD_CRK(*),NODENR(*),IADC_CRK(4,*),ELCRKINI(NLAY,NEL),
65 . ELCUTC(2,*),NODEDGE(2,*),CRKNODIAD(*),KNOD2ELC(*),XEDGE4N(4,*)
66 my_real DIR1(NLAY,NEL),DIR2(NLAY,NEL),A_I(NEL),
67 . XL2(NEL),YL2(NEL),XL3(NEL),YL3(NEL),XL4(NEL),YL4(NEL)
68
69 TYPE (ELBUF_STRUCT_), DIMENSION(NXEL) :: XFEM_STR
70 TYPE (XFEM_EDGE_) , DIMENSION(*) :: CRKEDGE
71C-----------------------------------------------
72C L o c a l V a r i a b l e s
73C-----------------------------------------------
74 INTEGER I, K, II, R, IR, P1, P2, PP1, PP2, PP3, NEWCRK, IED, OK,
75 . icrk,elcrk,nx1,nx2,nx3,nx4,nm,np,fac,ifi1,ifi2,iedge,
76 . icut,sigbeta,nod1,nod2,ie1,ie2,ie10,ie20,itri,iad1,iad2,iad3,iad4
77c
78 INTEGER ECUT(2,NEL),EDGEL(4,NEL),IENR0(4,NEL),dd(4),d(8),KPERM(8),
79 . N(4),NN(4),IADC(4),IENR(4),ISIGN(4)
80 INTEGER , DIMENSION(NEL) :: JCT,ELFISS,TIP
81c
82 my_real, DIMENSION(2) :: XMI,YMI
83 my_real, DIMENSION(2,NEL) :: XIN,YIN
84 my_real, DIMENSION(4,NEL) :: len,xxl,yyl,fit,beta0
85 my_real
86 . xint,yint,fi,cross,acd,bcd,dlx,dly,m12,mm,xint0,yint0,beta,
87 . bmin,bmax,dir11,dir22,x1,y1,x2,y2,x3,y3,x4,y4,area1,area2,area3
88C----------
89 DATA d/1,2,2,3,4,3,1,4/
90 DATA dd/2,3,4,1/
91 DATA kperm/1,2,3,4,1,2,3,4/
92 parameter(bmin = 0.01, bmax = 0.99)
93C=======================================================================
94 newcrk = 0
95 DO i=1,nel
96 jct(i) = 0
97 IF (elcrkini(ilay,i) == 1) THEN ! crack propagation
98 newcrk = newcrk + 1
99 jct(newcrk) = i
100 ENDIF
101 ENDDO
102 IF (newcrk == 0) RETURN
103C---
104 ii = nxel*(ilay-1)
105 pp1 = ii + 1
106 pp2 = ii + 2
107 pp3 = ii + 3
108C---
109 DO ir=1,newcrk
110 i = jct(ir)
111 elfiss(i)= 0
112 tip(i) = 0
113 ecut(1:2,i) = 0
114 edgel(1:4,i) = 0
115 ienr0(1:4,i) = 0
116 beta0(1:4,i) = zero
117 xin(1,i) = zero ! first inters point in local skew
118 yin(1,i) = zero
119 xin(2,i) = zero ! second inters point in local skew
120 yin(2,i) = zero
121c
122 xxl(1,i) = zero
123 yyl(1,i) = zero
124 xxl(2,i) = xl2(i)
125 yyl(2,i) = yl2(i)
126 xxl(3,i) = xl3(i)
127 yyl(3,i) = yl3(i)
128 xxl(4,i) = xl4(i)
129 yyl(4,i) = yl4(i)
130c
131 len(1,i) = xl2(i)*xl2(i) + yl2(i)*yl2(i)
132 len(2,i) = (xl3(i)-xl2(i))*(xl3(i)-xl2(i))+
133 . (yl3(i)-yl2(i))*(yl3(i)-yl2(i))
134 len(3,i) = (xl4(i)-xl3(i))*(xl4(i)-xl3(i))+
135 . (yl4(i)-yl3(i))*(yl4(i)-yl3(i))
136 len(4,i) = xl4(i)*xl4(i) + yl4(i)*yl4(i)
137 ENDDO
138C------------------------------------------------
139c First intersection (already cut edge)
140C------------------------------------------------
141 DO ir=1,newcrk ! loop over elems (layers) with advancing cracks
142 i = jct(ir)
143 elcrk = iel_crk(i+nft) ! DYE SYS XFEM
144 ok = 0
145 icut = 0
146 ied = 0
147 DO k=1,4 ! edges
148 iedge = xedge4n(k,elcrk)
149 IF (iedge > 0) THEN
150 icut = crkedge(ilay)%ICUTEDGE(iedge)
151 IF (icut == 1) THEN
152 nod1 = nodedge(1,iedge) ! noeud std
153 nod2 = nodedge(2,iedge)
154 IF (nod1 == ixc(k+1,i) .and. nod2 == ixc(dd(k)+1,i)) THEN
155 p1 = k
156 p2 = dd(k)
157 ELSE IF (nod2 == ixc(k+1,i) .and. nod1 == ixc(dd(k)+1,i)) THEN
158 p1 = dd(k)
159 p2 = k
160 ENDIF
161 ok = 1
162 ied = k
163 ecut(1,i)= k
164 EXIT
165 ENDIF ! IF (ICUT == 1) THEN
166 ENDIF
167 ENDDO ! DO K=1,4
168C---
169 IF (ok == 1) THEN ! edge found
170 beta = crkedge(ilay)%RATIO(iedge)
171 xin(1,i) = xxl(p1,i) + beta*(xxl(p2,i) - xxl(p1,i))
172 yin(1,i) = yyl(p1,i) + beta*(yyl(p2,i) - yyl(p1,i))
173c
174 elfiss(i) = crkedge(ilay)%EDGEICRK(iedge) ! ID of crack that propagates
175 ienr0(p1,i) = crkedge(ilay)%EDGEENR(1,iedge) ! enrichissement 1 noeud
176 ienr0(p2,i) = crkedge(ilay)%EDGEENR(2,iedge) ! enrichissement 2 noeud
177 edgel(ied,i) = 1 ! Local: First Edge Cup
178 iedge = xedge4n(ied,elcrk) ! N edge element sys xfem
179 tip(i) = crkedge(ilay)%EDGETIP(1,iedge) ! 1 or 2 , start or end of crack
180 ELSE
181 WRITE(iout,*) 'ERROR IN ADVANCING CRACK --- CHECK CRACK TIP'
182 CALL arret(2)
183 ENDIF
184C---
185 END DO ! DO IR=1,NEWCRK
186C--------------------------------------------------
187c Search for second intersection (new cut edge)
188C--------------------------------------------------
189 DO ir=1,newcrk
190 i = jct(ir)
191 elcrk = iel_crk(i+nft)
192 r = ecut(1,i)
193 xint0 = xin(1,i)
194 yint0 = yin(1,i)
195 dir11 =-dir2(ilay,i)
196 dir22 = dir1(ilay,i)
197C---
198 IF (dir11 == zero) THEN
199 DO k=1,3
200 r = kperm(ecut(1,i) + k)
201 iedge = xedge4n(r,elcrk)
202 nod1 = nodedge(1,iedge)
203 nod2 = nodedge(2,iedge)
204 IF (nod1 == ixc(r+1,i) .and. nod2 == ixc(dd(r)+1,i))THEN
205 p1 = r
206 p2 = dd(r)
207 ELSE IF (nod2 == ixc(r+1,i).and.nod1 == ixc(dd(r)+1,i))THEN
208 p1 = dd(r)
209 p2 = r
210 ENDIF
211 dlx = xxl(p2,i) - xxl(p1,i)
212 IF (dlx /= zero) THEN
213 dly = yyl(p2,i) - yyl(p1,i)
214 m12 = dly / dlx
215 xint = xint0
216 yint = yyl(p1,i) + m12*(xint-xxl(p1,i))
217 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero .and.
218 . (yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
219 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
220 beta = sqrt(cross / len(r,i))
221 IF (beta > bmax .OR. beta < bmin) THEN
222 beta = max(beta, bmin)
223 beta = min(beta, bmax)
224 yint = yyl(p1,i) + beta*(yyl(p2,i)-yyl(p1,i))
225 ENDIF
226C
227 ecut(2,i) = r
228 xin(2,i) = xint
229 yin(2,i) = yint
230 edgel(r,i) = 2
231 beta0(r,i) = beta
232 EXIT
233 ENDIF
234 ENDIF
235 ENDDO
236c
237 ELSEIF (dir22 == zero) THEN
238 DO k=1,3
239 r = kperm(ecut(1,i) + k)
240 iedge = xedge4n(r,elcrk)
241 nod1 = nodedge(1,iedge)
242 nod2 = nodedge(2,iedge)
243 IF (nod1 == ixc(r+1,i) .and. nod2 == ixc(dd(r)+1,i)) THEN
244 p1 = r
245 p2 = dd(r)
246 ELSE IF (nod2 == ixc(r+1,i).and.nod1 == ixc(dd(r)+1,i)) THEN
247 p1 = dd(r)
248 p2 = r
249 ENDIF
250 dly = yyl(p2,i) - yyl(p1,i)
251 IF (dly /= zero) THEN
252 dlx = xxl(p2,i) - xxl(p1,i)
253 m12 = dlx / dly
254 yint = yint0
255 xint = xxl(p1,i) + m12*(yint-yyl(p1,i))
256 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero .and.
257 . (yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
258 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
259 beta = sqrt(cross / len(r,i))
260 IF (beta > bmax .OR. beta < bmin) THEN
261 beta = max(beta, bmin)
262 beta = min(beta, bmax)
263 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
264 ENDIF
265C
266 ecut(2,i) = r
267 xin(2,i) = xint
268 yin(2,i) = yint
269 edgel(r,i) = 2
270 beta0(r,i) = beta
271 EXIT
272 ENDIF
273 ENDIF
274 ENDDO
275c
276 ELSEIF (dir11 /= zero .and. dir22 /= zero) THEN
277 DO k=1,3
278 r = kperm(ecut(1,i) + k)
279 iedge = xedge4n(r,elcrk)
280 nod1 = nodedge(1,iedge)
281 nod2 = nodedge(2,iedge)
282 IF (nod1 == ixc(r+1,i) .and. nod2 == ixc(dd(r)+1,i)) THEN
283 p1 = r
284 p2 = dd(r)
285 ELSE IF (nod2 == ixc(r+1,i).and.nod1 == ixc(dd(r)+1,i)) THEN
286 p1 = dd(r)
287 p2 = r
288 ENDIF
289C
290 dlx = xxl(p2,i) - xxl(p1,i)
291 dly = yyl(p2,i) - yyl(p1,i)
292 mm = dir22/dir11
293 IF (dlx == zero) THEN
294 xint = xxl(p1,i)
295 yint = yint0 + mm*(xint-xint0)
296 IF ((yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
297 cross = (yyl(p1,i) - yint)**2
298 beta = sqrt(cross / len(r,i))
299 IF (beta > bmax .OR. beta < bmin) THEN
300 beta = max(beta, bmin)
301 beta = min(beta, bmax)
302 yint = yyl(p1,i) + beta*(yyl(p2,i)-yyl(p1,i))
303 ENDIF
304 ecut(2,i) = r
305 xin(2,i) = xint
306 yin(2,i) = yint
307 edgel(r,i) = 2
308 beta0(r,i) = beta
309 EXIT
310 ENDIF
311 ELSEIF (dly == zero) THEN
312 yint = yyl(p1,i)
313 xint = xint0 + (yint0-yyl(p1,i)) / mm
314 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero) THEN
315 cross = (xxl(p1,i) - xint)**2
316 beta = sqrt(cross / len(r,i))
317 IF (beta > bmax .OR. beta < bmin) THEN
318 beta = max(beta, bmin)
319 beta = min(beta, bmax)
320 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
321 ENDIF
322 ecut(2,i) = r
323 xin(2,i) = xint
324 yin(2,i) = yint
325 edgel(r,i) = 2
326 beta0(r,i) = beta
327 EXIT
328 ENDIF
329 ELSE
330 m12 = dly / dlx
331 IF (mm /= m12) THEN
332 xint = (yint0-yyl(p1,i) + m12*xxl(p1,i) - mm*xint0)/(m12-mm)
333 yint = yint0 + mm*(xint-xint0)
334 acd = (yint-yyl(p1,i))*(xint0 - xxl(p1,i))
335 . - (xint-xxl(p1,i))*(yint0 - yyl(p1,i))
336 bcd = (yint-yyl(p2,i))*(xint0 - xxl(p2,i))
337 . - (xint-xxl(p2,i))*(yint0 - yyl(p2,i))
338 IF (acd*bcd <= zero) THEN
339 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
340 beta = sqrt(cross / len(r,i))
341 IF (beta > bmax .OR. beta < bmin) THEN
342 beta = max(beta, bmin)
343 beta = min(beta, bmax)
344 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
345 yint = yyl(p1,i) + beta*(yyl(p2,i)-yyl(p1,i))
346 ENDIF
347 ecut(2,i) = r
348 xin(2,i) = xint
349 yin(2,i) = yint
350 edgel(r,i) = 2
351 beta0(r,i) = beta
352 EXIT
353 ENDIF
354 ENDIF
355 ENDIF
356 ENDDO
357 ENDIF
358 ENDDO
359C-----------------------------------------------------------------------
360C check for getting both intersections
361C-----------------------------------------------------------------------
362 DO ir=1,newcrk
363 i=jct(ir)
364 fac = 0
365 DO r=1,4
366 IF (edgel(r,i)==1 .or. edgel(r,i)==2) fac=fac+1
367 ENDDO
368 IF (fac /= 2) THEN
369 WRITE(iout,*) 'ERROR IN ADVANCING CRACK. NO CUT EDGES'
370 CALL arret(2)
371 ENDIF
372 ENDDO
373c tag cut edge numbers of each layer
374 DO ir=1,newcrk
375 i = jct(ir)
376 elcrk = iel_crk(i+nft)
377 DO k=1,4
378 ied = edgel(k,i)
379 crkedge(ilay)%IEDGEC(k,elcrk) = ied ! IED = 0,1,2
380 ENDDO
381 ENDDO
382C----------------------------------------------------------------------
383C SIGN DISTANCE FOR NEW CRACKED LAYER
384C----------------------------------------------------------------------
385 DO ir=1,newcrk
386 i = jct(ir)
387 fit(1,i) = zero
388 fit(2,i) = zero
389 fit(3,i) = zero
390 fit(4,i) = zero
391c
392 DO k=1,4
393 p1 = k
394 p2 = dd(k)
395 ied = edgel(k,i)
396 IF (ied > 0) THEN
397 xmi(ied) = half*(xxl(p1,i)+xxl(p2,i))
398 ymi(ied) = half*(yyl(p1,i)+yyl(p2,i))
399 ENDIF
400 ENDDO
401 DO k=1,4
402 fi = zero
403 CALL lsint4(xmi(1),ymi(1),xmi(2),ymi(2),xxl(k,i),yyl(k,i),fi)
404 IF (fit(k,i) == zero) fit(k,i)=fi
405 ENDDO
406 ENDDO
407C-------------------
408c Loop over newly cut elements
409C-------------------
410 DO ir=1,newcrk
411 i = jct(ir)
412 elcrk = iel_crk(i+nft)
413C
414 k = ecut(2,i) ! second intersection (tip)
415 iedge = xedge4n(k,elcrk)
416 icut = crkedge(ilay)%ICUTEDGE(iedge)
417 IF (icut > 0) THEN ! already cut before => edge connecting two cracks
418 crkedge(ilay)%ICUTEDGE(iedge) = 3 ! 2 cracks on the same edge
419 ELSE
420 crkedge(ilay)%ICUTEDGE(iedge) = 2 ! edge cut
421 crkedge(ilay)%RATIO(iedge) = beta0(k,i)
422 ENDIF
423 ENDDO ! IR=1,NEWCRK
424C-----------------------
425C FILL new cut layer - main loop
426C-----------------------
427 DO ir=1,newcrk
428 i = jct(ir)
429 elcrk = iel_crk(i+nft) ! n elem sys xfem
430 elcutc(1,i) = 2 ! FLAG of the std cut element (= 1.2)
431 numelcrk = numelcrk + 1 ! + 1 cut element (for debug/animation by element)
432C
433 iadc(1) = iadc_crk(1,elcrk)
434 iadc(2) = iadc_crk(2,elcrk)
435 iadc(3) = iadc_crk(3,elcrk)
436 iadc(4) = iadc_crk(4,elcrk)
437c
438c if(ispmd==6.and.elcrk==1) print*,' ELCRK:,IAD=',ELCRK,IADC(1),IADC(2),IADC(3),IADC(4)
439c if(ispmd==3.and.elcrk==1220) print*,' ELCRK:,IAD=',ELCRK,IADC(1),IADC(2),IADC(3),IADC(4)
440c if(ispmd==1.and.elcrk==595) print*,' ELCRK:,IAD=',ELCRK,IADC(1),IADC(2),IADC(3),IADC(4)
441c if(ispmd==1.and.elcrk==636) print*,' ELCRK:,IAD=',ELCRK,IADC(1),IADC(2),IADC(3),IADC(4)
442c if(nspmd==1.and.(elcrk==4059.or.elcrk==3959.or.elcrk==3960.or.elcrk==3860))
443c
444c
445 n(1) = ixc(2,i) ! n node sys std
446 n(2) = ixc(3,i)
447 n(3) = ixc(4,i)
448 n(4) = ixc(5,i)
449C
450 nn(1) = inod_crk(n(1)) ! n node sys xfem
451 nn(2) = inod_crk(n(2))
452 nn(3) = inod_crk(n(3))
453 nn(4) = inod_crk(n(4))
454C
455 icrk = elfiss(i) ! ID of crack that propagates
456C
457 isign(1) = int(sign(one,fit(1,i)))
458 isign(2) = int(sign(one,fit(2,i)))
459 isign(3) = int(sign(one,fit(3,i)))
460 isign(4) = int(sign(one,fit(4,i)))
461C
462 IF (fit(1,i) == zero) isign(1) = 0
463 IF (fit(2,i) == zero) isign(2) = 0
464 IF (fit(3,i) == zero) isign(3) = 0
465 IF (fit(4,i) == zero) isign(4) = 0
466c--------------------------------------------
467c activate enriched nodes
468c--------------------------------------------
469 DO k=1,4
470 IF (ienr0(k,i) == 0) THEN
471 ienr(k) = crknodiad(iadc(k)) + knod2elc(nn(k))*(ilay-1) ! default => node std
472 ELSE
473 ienr(k) = ienr0(k,i) ! enriched node
474 ENDIF
475 ENDDO
476C
477 sigbeta = 0
478 k = ecut(2,i) ! only new cut edge (tip)
479 iedge = xedge4n(k,elcrk)
480 nod1 = nodedge(1,iedge)
481 nod2 = nodedge(2,iedge)
482 ie10 = crkedge(ilay)%EDGEENR(1,iedge)
483 ie20 = crkedge(ilay)%EDGEENR(2,iedge)
484 IF (nod1 == ixc(k+1,i) .and. nod2 == ixc(dd(k)+1,i)) THEN
485 ie1 = ienr(k)
486 ie2 = ienr(dd(k))
487 ifi1 = isign(k)
488 ifi2 = isign(dd(k))
489 sigbeta = iedge
490 ELSE IF (nod2 == ixc(k+1,i) .and. nod1 == ixc(dd(k)+1,i)) THEN
491 ie1 = ienr(dd(k))
492 ie2 = ienr(k)
493 ifi1 = isign(dd(k))
494 ifi2 = isign(k)
495 sigbeta = -iedge
496 END IF
497 crkedge(ilay)%EDGEENR(1,iedge) = max(ie1,ie10)
498 crkedge(ilay)%EDGEENR(2,iedge) = max(ie2,ie20)
499 IF (crkedge(ilay)%EDGEICRK(iedge) == 0)
500 . crkedge(ilay)%EDGEICRK(iedge) = icrk
501c--------------------------
502 itri = 0
503 nm = 0
504 np = 0
505 DO k=1,4
506 ied = edgel(k,i)
507 IF (ied > 0) THEN
508 iedge = xedge4n(k,elcrk)
509 crkedge(ilay)%EDGETIP(1,iedge) = tip(i)
510 crkedge(ilay)%EDGETIP(2,iedge) =
511 . crkedge(ilay)%EDGETIP(2,iedge) + 1
512 ENDIF
513 IF (isign(k) > 0) THEN
514 itri = itri + 1
515 np = k
516 ELSEIF (isign(k) < 0) THEN
517 nm = k
518 ENDIF
519 ENDDO
520c---
521 IF (itri == 1) THEN
522 itri = -1
523 nx1 = np
524 ELSEIF (itri == 3) THEN
525 itri = 1
526 nx1 = nm
527 ELSEIF (itri == 2) THEN
528 itri = 0
529 IF (np > 1 .and. isign(np-1) > 0) THEN
530 nx1 = np-1
531 ELSE
532 nx1 = np
533 ENDIF
534 ENDIF
535c
536c print*,' Avancement: IEL,ELCRK,proc=',I+NFT,ELCRK,ispmd
537c print*,' IAD=',IADC(1),IADC(2),IADC(3),IADC(4)
538c print*,' NSX=',NN(1),NN(2),NN(3),NN(4),ITRI
539c---
540 nx2 = kperm(nx1+1)
541 nx3 = kperm(nx1+2)
542 nx4 = kperm(nx1+3)
543C
544 crkedge(ilay)%LAYCUT(elcrk) = 1
545 xfem_phantom(ilay)%ITRI(1,elcrk) = itri
546 xfem_phantom(ilay)%ITRI(2,elcrk) = nx1
547 xfem_phantom(ilay)%ELCUT(elcrk) = icrk
548 xfem_phantom(ilay)%IFI(iadc(1)) = isign(1)
549 xfem_phantom(ilay)%IFI(iadc(2)) = isign(2)
550 xfem_phantom(ilay)%IFI(iadc(3)) = isign(3)
551 xfem_phantom(ilay)%IFI(iadc(4)) = isign(4)
552C------------------
553C IXEL = 1 => positif element within ILAY (ILEV = PP1)
554C------------------
555 crklvset(pp1)%ENR0(1,iadc(1)) = -ienr(1)
556 crklvset(pp1)%ENR0(1,iadc(2)) = -ienr(2)
557 crklvset(pp1)%ENR0(1,iadc(3)) = -ienr(3)
558 crklvset(pp1)%ENR0(1,iadc(4)) = -ienr(4)
559C
560 IF (isign(1) > 0) crklvset(pp1)%ENR0(1,iadc(1)) = 0
561 IF (isign(2) > 0) crklvset(pp1)%ENR0(1,iadc(2)) = 0
562 IF (isign(3) > 0) crklvset(pp1)%ENR0(1,iadc(3)) = 0
563 IF (isign(4) > 0) crklvset(pp1)%ENR0(1,iadc(4)) = 0
564C------------------
565C IXEL = 2 => negatif element within ILAY (ILEV = PP2)
566C------------------
567 crklvset(pp2)%ENR0(1,iadc(1)) = -ienr(1)
568 crklvset(pp2)%ENR0(1,iadc(2)) = -ienr(2)
569 crklvset(pp2)%ENR0(1,iadc(3)) = -ienr(3)
570 crklvset(pp2)%ENR0(1,iadc(4)) = -ienr(4)
571C
572 IF (isign(1) < 0) crklvset(pp2)%ENR0(1,iadc(1)) = 0
573 IF (isign(2) < 0) crklvset(pp2)%ENR0(1,iadc(2)) = 0
574 IF (isign(3) < 0) crklvset(pp2)%ENR0(1,iadc(3)) = 0
575 IF (isign(4) < 0) crklvset(pp2)%ENR0(1,iadc(4)) = 0
576C------------------
577C IXEL = 3 => Third phantom element (ILEV = PP3)
578C------------------
579 IF (itri == 0) THEN
580 x1 = xxl(nx1,i)
581 y1 = yyl(nx1,i)
582 x2 = xxl(nx2,i)
583 y2 = yyl(nx2,i)
584 ied = edgel(nx2,i)
585 IF (ied > 0) THEN
586 x3 = xin(ied,i)
587 y3 = yin(ied,i)
588 ELSE
589 ENDIF
590 ied = edgel(nx4,i)
591 IF (ied > 0) THEN
592 x4 = xin(ied,i)
593 y4 = yin(ied,i)
594 ELSE
595 ENDIF
596 area1 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
597 .
598 area1 = area1 * a_i(i)
599 area2 = one - area1
600 area3 = zero
601c
602 xfem_str(3)%GBUF%OFF(i) = zero
603c XFEM_STR(3)%BUFLY(ILAY)%OFF(I) = 0
604 xfem_str(3)%BUFLY(ilay)%LBUF(1,1,1)%OFF(i) = zero
605c
606 ELSEIF (itri < 0) THEN ! cut en 3, IPH3 = IPH2 < 0
607 ie1 = xedge4n(nx2,elcrk)
608 ie2 = xedge4n(nx4,elcrk)
609c
610 IF (crkedge(ilay)%ICUTEDGE(ie2) > 1) THEN
611 sigbeta = -sigbeta
612 iad1 = iadc(nx1)
613 iad2 = iadc(nx2)
614 iad3 = iadc(nx3)
615 iad4 = iadc(nx4)
616 nod1 = iad4
617 nod2 = iad1
618 crklvset(pp3)%ENR0(1,iad1) = abs(crklvset(pp2)%ENR0(1,iad1))
619 crklvset(pp3)%ENR0(1,iad2) = crklvset(pp2)%ENR0(1,iad2)
620 crklvset(pp3)%ENR0(1,iad3) = crklvset(pp2)%ENR0(1,iad3)
621 crklvset(pp3)%ENR0(1,iad4) = crklvset(pp2)%ENR0(1,iad4)
622 crklvset(pp2)%ENR0(1,iad1) = -crknodiad(iad1) - knod2elc(nn(nx1))*(ilay-1)
623c
624c-- IXEL=2 AREA factor
625 ied = crkedge(ilay)%IEDGEC(nx4,elcrk)
626 x1 = xin(ied,i)
627 y1 = yin(ied,i)
628 x2 = xxl(nx3,i)
629 y2 = yyl(nx3,i)
630 x3 = xxl(nx4,i)
631 y3 = yyl(nx4,i)
632 area2 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
633c-- IXEL=1 AREA factor
634 ied = crkedge(ilay)%IEDGEC(nx1,elcrk)
635 x2 = xin(ied,i)
636 y2 = yin(ied,i)
637 x3 = xxl(nx1,i)
638 y3 = yyl(nx1,i)
639 area1 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
640c
641 area1 = area1 * a_i(i)
642 area2 = area2 * a_i(i)
643 area3 = one - area1 - area2
644 ELSEIF (crkedge(ilay)%ICUTEDGE(ie1) > 1) THEN
645 print*,' IMPOSSIBLE CASE'
646 ENDIF
647C
648 ELSEIF (itri > 0) THEN ! cut en 3, IPH3 = IPH1
649 ie1 = xedge4n(nx1,elcrk)
650 ie2 = xedge4n(nx4,elcrk)
651c
652 IF (crkedge(ilay)%ICUTEDGE(ie1) > 1) THEN
653 iad1 = iadc(nx1)
654 iad2 = iadc(nx2)
655 iad3 = iadc(nx3)
656 iad4 = iadc(nx4)
657 nod1 = iad2
658 nod2 = iad1
659 crklvset(pp3)%ENR0(1,iad1) = abs(crklvset(pp1)%ENR0(1,iad1))
660 crklvset(pp3)%ENR0(1,iad2) = crklvset(pp1)%ENR0(1,iad2)
661 crklvset(pp3)%ENR0(1,iad3) = crklvset(pp1)%ENR0(1,iad3)
662 crklvset(pp3)%ENR0(1,iad4) = crklvset(pp1)%ENR0(1,iad4)
663 crklvset(pp1)%ENR0(1,iad1) = -crknodiad(iad1) - knod2elc(nn(nx1))*(ilay-1)
664c-- IXEL=1 AREA factor
665 ied = crkedge(ilay)%IEDGEC(nx1,elcrk)
666 x1 = xin(ied,i)
667 y1 = yin(ied,i)
668 x2 = xxl(nx2,i)
669 y2 = yyl(nx2,i)
670 x3 = xxl(nx3,i)
671 y3 = yyl(nx3,i)
672 area1 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
673c-- IXEL=2 AREA factor
674 ied = crkedge(ilay)%IEDGEC(nx4,elcrk)
675 x2 = xin(ied,i)
676 y2 = yin(ied,i)
677 x3 = xxl(nx1,i)
678 y3 = yyl(nx1,i)
679 area2 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
680c
681 area1 = area1 * a_i(i)
682 area2 = area2 * a_i(i)
683 area3 = one - area1 - area2
684 ELSEIF (crkedge(ilay)%ICUTEDGE(ie2) > 1) THEN
685 print*,' IMPOSSIBLE CASE'
686 ENDIF
687 ENDIF ! ITRI
688c----
689 crklvset(pp1)%AREA(elcrk) = area1
690 crklvset(pp2)%AREA(elcrk) = area2
691 crklvset(pp3)%AREA(elcrk) = area3
692c
693 IF (area3 < zero .or. area1 > one .or. area2 > one .or. area3 > one ) THEN
694 print*,'ERROR : XFEM PHANTOM ELEMENT AREA: ELCRK=',elcrk
695 ENDIF
696c
697 ENDDO ! IR=1,NEWCRK
698c-----------
699 RETURN
700 END
701!||====================================================================
702!|| lsint4 ../engine/source/elements/xfem/crklayer4n_adv.F
703!||--- called by ------------------------------------------------------
704!|| crklayer3n_adv ../engine/source/elements/xfem/crklayer3n_adv.f
705!|| crklayer3n_ini ../engine/source/elements/xfem/crklayer3n_ini.F
706!|| crklayer4n_adv ../engine/source/elements/xfem/crklayer4n_adv.F
707!|| crklayer4n_ini ../engine/source/elements/xfem/crklayer4n_ini.F
708!|| enrichc_ini ../engine/source/elements/xfem/enrichc_ini.F
709!|| enrichtg_ini ../engine/source/elements/xfem/enrichtg_ini.F
710!||====================================================================
711 SUBROUTINE lsint4(Y1, Z1, Y2, Z2, Y, Z, FI)
712C-----------------------------------------------
713C I m p l i c i t T y p e s
714C-----------------------------------------------
715#include "implicit_f.inc"
716C-----------------------------------------------
717C D u m m y A r g u m e n t s
718C-----------------------------------------------
719 my_real
720 . y1,z1,y2,z2,y,z,fi
721C-----------------------------------------------
722C L o c a l V a r i a b l e s
723C-----------------------------------------------
724 my_real
725 . area,ab
726C-----------------------------------------------
727 fi = zero
728 area = ((y2*z-y*z2)-(y1*z-y*z1)+(y1*z2-z1*y2))
729 ab = (y2-y1)**2 + (z2-z1)**2
730 IF (ab > zero) fi = area/sqrt(ab)
731C-----------
732 RETURN
733 END
my_real
#define my_real
Definition cppsort.cpp:32
crklayer3n_adv
subroutine crklayer3n_adv(nel, nft, ilay, nlay, ixtg, elcutc, elcrkini, iel_crktg, inod_crk, iad_crktg, nodenr, dir1, dir2, nodedge, crknodiad, knod2elc, crkedge, xedge3n, ngl, area, xl2, xl3, yl2, yl3)
Definition crklayer3n_adv.F:40
lsint4
subroutine lsint4(y1, z1, y2, z2, y, z, fi)
Definition crklayer4n_adv.F:712
crklayer4n_adv
subroutine crklayer4n_adv(xfem_str, nel, nft, ixc, elcutc, ilay, nlay, iel_crk, inod_crk, iadc_crk, nodenr, elcrkini, dir1, dir2, nodedge, crknodiad, knod2elc, crkedge, a_i, xl2, xl3, xl4, yl2, yl3, yl4, xedge4n, ngl)
Definition crklayer4n_adv.F:43
czforc3
subroutine czforc3(timers, elbuf_str, jft, jlt, nft, npt, itab, mtn, ipri, ithk, neltst, istrain, ipla, dt1, dt2t, pm, geo, partsav, ixc, ityptst, bufmat, tf, npf, iadc, failwave, x, dr, v, vr, f, m, stifn, stifr, fsky, tani, indxof, ismstr, group_param, ipartc, thke, nvc, iofc, ihbe, f11, f12, f13, f14, f21, f22, f23, f24, f31, f32, f33, f34, m11, m12, m13, m14, m21, m22, m23, m24, m31, m32, m33, m34, kfts, fzero, igeo, ipm, ifailure, itask, jthe, temp, fthe, fthesky, iexpan, gresav, grth, igrth, xedge4n, msc, dmelc, jsms, table, iparg, mat_elem, ixfem, knod2elc, sensors, elcutc, inod_crk, iel_crk, nodenr, iadc_crk, nodedge, crknodiad, condn, condnsky, stack, isubstack, xfem_str, crkedge, drape_sh4n, nel, nloc_dmg, indx_drape, igre, jtur, dt, ncycle, snpc, stf, glob_therm, idel7nok, userl_avail, maxfunc, sbufmat, ipart, lipart1)
Definition czforc3.F:116
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
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
crackxfem_mod::xfem_phantom
type(xfem_phantom_), dimension(:), allocatable xfem_phantom
Definition crackxfem_mod.F:139
crackxfem_mod::crklvset
type(xfem_lvset_), dimension(:), allocatable crklvset
Definition crackxfem_mod.F:146
arret
subroutine arret(nn)
Definition arret.F:86