OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
crklayer4n_adv.F File Reference
#include "implicit_f.inc"
#include "units_c.inc"
#include "com_xfem1.inc"

Go to the source code of this file.

Functions/Subroutines

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)
subroutine lsint4 (y1, z1, y2, z2, y, z, fi)

Function/Subroutine Documentation

◆ crklayer4n_adv()

subroutine crklayer4n_adv ( type (elbuf_struct_), dimension(nxel) xfem_str,
integer nel,
integer nft,
integer, dimension(nixc,*) ixc,
integer, dimension(2,*) elcutc,
integer ilay,
integer nlay,
integer, dimension(*) iel_crk,
integer, dimension(*) inod_crk,
integer, dimension(4,*) iadc_crk,
integer, dimension(*) nodenr,
integer, dimension(nlay,nel) elcrkini,
dir1,
dir2,
integer, dimension(2,*) nodedge,
integer, dimension(*) crknodiad,
integer, dimension(*) knod2elc,
type (xfem_edge_), dimension(*) crkedge,
a_i,
xl2,
xl3,
xl4,
yl2,
yl3,
yl4,
integer, dimension(4,*) xedge4n,
integer, dimension(nel) ngl )

Definition at line 36 of file crklayer4n_adv.F.

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.or..or..or. IF (AREA3 < ZERO AREA1 > ONE AREA2 > ONE AREA3 > ONE ) THEN
694 print*,'error : xfem phantom element area: elcrk=',ELCRK
695 ENDIF
696c
697 ENDDO ! IR=1,NEWCRK
698c-----------
699 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
lsint4
subroutine lsint4(y1, z1, y2, z2, y, z, fi)
Definition crklayer4n_adv.F:712
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

◆ lsint4()

subroutine lsint4 ( y1,
z1,
y2,
z2,
y,
z,
fi )

Definition at line 711 of file crklayer4n_adv.F.

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