OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
crklayer3n_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!|| crklayer3n_adv ../engine/source/elements/xfem/crklayer3n_adv.F
25!||--- called by ------------------------------------------------------
26!|| c3forc3 ../engine/source/elements/sh3n/coque3n/c3forc3.F
27!||--- calls -----------------------------------------------------
28!|| arret ../engine/source/system/arret.F
29!|| lsint4 ../engine/source/elements/xfem/crklayer4n_adv.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 crklayer3n_adv(
35 . NEL ,NFT ,ILAY ,NLAY ,IXTG ,
36 . ELCUTC ,ELCRKINI ,IEL_CRKTG,INOD_CRK ,IAD_CRKTG,
37 . NODENR ,DIR1 ,DIR2 ,NODEDGE ,CRKNODIAD,
38 . KNOD2ELC ,CRKEDGE ,XEDGE3N ,NGL ,AREA ,
39 . XL2 ,XL3 ,YL2 ,YL3 )
40c-----------------------------------------------
41C crack advancement, shells 3N
42c-----------------------------------------------
43C M o d u l e s
44C-----------------------------------------------
45 USE crackxfem_mod
46 use element_mod , only : nixtg
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"
55#include "com_xfem1.inc"
56C-----------------------------------------------
57C D u m m y A r g u m e n t s
58C-----------------------------------------------
59 INTEGER NEL,NFT,ILAY,NLAY
60 INTEGER IXTG(NIXTG,*),NGL(NEL),IEL_CRKTG(*),
61 . INOD_CRK(*),NODENR(*),IAD_CRKTG(3,*),ELCRKINI(NLAY,*),
62 . ELCUTC(2,*),NODEDGE(2,*),CRKNODIAD(*),KNOD2ELC(*),XEDGE3N(3,*)
63 my_real, DIMENSION(NEL) :: AREA,XL2,YL2,XL3,YL3
64 my_real dir1(nlay,nel),dir2(nlay,nel)
65 TYPE (XFEM_EDGE_) , DIMENSION(*) :: CRKEDGE
66C-----------------------------------------------
67C L o c a l V a r i a b l e s
68C-----------------------------------------------
69 INTEGER I,K,IR,p1,p2,
70 . ienr1,ienr2,newcrk,ied,ied1,ied2,fac,ok,icrk,pp1,pp2,pp3,
71 . nod1,nod2,ifi1,ifi2,ie1,ie2,ie10,ie20,elcrk,elcrktg,
72 . iedge,icut,sigbeta,itri,nm,np,nx1,nx2,nx3,iad1,iad2,iad3
73 INTEGER JCT(NEL),ELFISS(NEL),EDGEL(3,NEL),
74 . ienr0(3,nel),tip(nel),iadc(3),
75 . dd(3),d(6),dx(6),isign(3),n(3),ienr(3),nn(3),inv(2)
76C----------
77 my_real, DIMENSION(NEL) :: xl1,yl1
78 my_real, DIMENSION(2,NEL) :: xin,yin
79 my_real, DIMENSION(3,NEL) :: xxl,yyl,len,fit
80 my_real beta0(3,nel),xn(3),yn(3),xmi(2),ymi(2)
81 my_real beta,xint,yint,fi,bmin,bmax,
82 . x10,y10,x20,y20,m12,mm,cross1,cross12,
83 . xint0,yint0,dir11,dir22,x1,y1,x2,y2,x3,y3,area1,area2,area3
84c------------------------
85 DATA d/1,2,2,3,1,3/
86 DATA dd/2,3,1/
87 DATA dx/1,2,3,1,2,3/
88 parameter(bmin = 0.01, bmax = 0.99)
89C=======================================================================
90 newcrk = 0
91 DO i=1,nel
92 jct(i) = 0
93 IF (elcrkini(ilay,i) == 1) THEN
94 newcrk = newcrk + 1
95 jct(newcrk) = i
96 ENDIF
97 ENDDO
98 IF (newcrk == 0) RETURN
99C---
100 pp1 = nxel*(ilay-1) + 1
101 pp2 = pp1 + 1
102 pp3 = pp1 + 2
103C---
104 DO i=1,nel
105 beta0(1:3,i) = zero
106 elfiss(i)= 0
107 tip(i) = 0
108C
109 edgel(1,i)=0
110 edgel(2,i)=0
111 edgel(3,i)=0
112C
113 ienr0(1,i)=0
114 ienr0(2,i)=0
115 ienr0(3,i)=0
116 xin(1,i) = zero ! first inters point in local skew
117 yin(1,i) = zero
118 xin(2,i) = zero ! second inters point in local skew
119 yin(2,i) = zero
120 ENDDO
121C
122 inv(1) = 2
123 inv(2) = 1
124C
125C---
126C advance crack inside uncut elements in the layer
127C---
128 DO ir=1,newcrk ! loop over elements with advancing cracks
129 i = jct(ir)
130 elcrktg = iel_crktg(i+nft)
131 ok = 0
132 icut = 0
133 ied = 0
134 DO k=1,3
135 iedge = xedge3n(k,elcrktg)
136 icut = crkedge(ilay)%ICUTEDGE(iedge)
137 nod1 = nodedge(1,iedge)
138 nod2 = nodedge(2,iedge)
139 IF (nod1 == ixtg(k+1,i) .and. nod2 == ixtg(dd(k)+1,i)) THEN
140 p1 = k
141 p2 = dd(k)
142 ELSE IF (nod2 == ixtg(k+1,i) .and. nod1 == ixtg(dd(k)+1,i)) THEN
143 p1 = dd(k)
144 p2 = k
145 ENDIF
146 IF (icut == 1) THEN
147 ok = ok + 1
148 ied = k
149c tag
150 icrk = crkedge(ilay)%EDGEICRK(iedge)
151 ienr1 = crkedge(ilay)%EDGEENR(1,iedge)
152 ienr2 = crkedge(ilay)%EDGEENR(2,iedge)
153 ifi1 = crkedge(ilay)%EDGEIFI(1,iedge)
154 ifi2 = crkedge(ilay)%EDGEIFI(2,iedge)
155c fill
156 elfiss(i) = icrk
157 ienr0(p1,i)= ienr1
158 ienr0(p2,i)= ienr2
159 EXIT
160 ENDIF ! IF(ICUT == 1)THEN
161 ENDDO ! DO K=1,4
162C---
163 IF (ok /= 1) THEN
164 WRITE(iout,*) 'ERROR IN ADVANCING CRACK --- CHECK CRACK TIP'
165 CALL arret(2)
166 ENDIF
167C---
168 edgel(ied,i) = 1
169 iedge = xedge3n(ied,elcrktg)
170 tip(i) = crkedge(ilay)%EDGETIP(1,iedge)
171C---
172 END DO ! DO IR=1,NEWCRK
173C----------------------------------------------------------------------
174c local coords
175C----------------------------------------------------------------------
176 DO i=1,nel
177 xl1(i) = zero
178 yl1(i) = zero
179 xxl(1,i) = xl1(i)
180 yyl(1,i) = yl1(i)
181 xxl(2,i) = xl2(i)
182 yyl(2,i) = yl2(i)
183 xxl(3,i) = xl3(i)
184 yyl(3,i) = yl3(i)
185 ENDDO
186c
187 DO i=1,nel
188 len(1,i) = (xl2(i)-xl1(i))*(xl2(i)-xl1(i))
189 . + (yl2(i)-yl1(i))*(yl2(i)-yl1(i))
190 len(2,i) = (xl3(i)-xl2(i))*(xl3(i)-xl2(i))
191 . + (yl3(i)-yl2(i))*(yl3(i)-yl2(i))
192 len(3,i) = (xl1(i)-xl3(i))*(xl1(i)-xl3(i))
193 . + (yl1(i)-yl3(i))*(yl1(i)-yl3(i))
194 ENDDO
195C------------------------------------------------
196C - intersections -
197C------------------------------------------------
198 DO ir=1,newcrk
199 i=jct(ir)
200 elcrktg = iel_crktg(i+nft)
201 elcrk = elcrktg + ecrkxfec
202 ied1 = 0
203 ied2 = 0
204 DO k=1,3
205 IF(edgel(k,i) > 0)THEN
206 ied1 = edgel(k,i)
207 ied2 = inv(ied1)
208 EXIT
209 END IF
210 END DO
211 DO k=1,3
212 iedge = xedge3n(k,elcrktg)
213 IF (iedge > 0 .and. edgel(k,i) == 1) THEN
214 icut = crkedge(ilay)%ICUTEDGE(iedge)
215 IF (icut == 1) THEN
216 beta = crkedge(ilay)%RATIO(iedge)
217C
218 IF (beta > one .or. beta == zero) THEN
219 WRITE(*,*) 'ERROR NEGATIV BETA, NO INTERSECTION!'
220 CALL arret(2)
221 ENDIF
222C
223 nod1 = nodedge(1,iedge)
224 nod2 = nodedge(2,iedge)
225 IF (nod1 == ixtg(k+1,i) .and. nod2 == ixtg(dd(k)+1,i)) THEN
226 p1 = k
227 p2 = dd(k)
228 ELSEIF (nod2 == ixtg(k+1,i).and.nod1==ixtg(dd(k)+1,i)) THEN
229 p1 = dd(k)
230 p2 = k
231 ENDIF
232 x10 = xxl(p1,i)
233 y10 = yyl(p1,i)
234 x20 = xxl(p2,i)
235 y20 = yyl(p2,i)
236C
237 xint = x10+beta*(x20-x10)
238 yint = y10+beta*(y20-y10)
239 xin(ied1,i) = xint
240 yin(ied1,i) = yint
241 ENDIF
242 ENDIF
243 ENDDO
244C---
245 IF (ied1 == 0 .or. ied2 == 0) GOTO 130
246 xint0 = xin(ied1,i)
247 yint0 = yin(ied1,i)
248C---
249 dir11 = -dir2(ilay,i)
250 dir22 = dir1(ilay,i)
251C---
252 IF (dir11 == zero) THEN
253 DO 140 k=1,3
254 xint = zero
255 yint = zero
256 elcrktg = iel_crktg(i+nft)
257 elcrk = elcrktg + ecrkxfec
258 iedge = xedge3n(k,elcrktg)
259 nod1 = nodedge(1,iedge)
260 nod2 = nodedge(2,iedge)
261 IF(nod1 == ixtg(k+1,i) .and. nod2 == ixtg(dd(k)+1,i))THEN
262 p1 = k
263 p2 = dd(k)
264 ELSE IF(nod2 == ixtg(k+1,i).and.nod1==ixtg(dd(k)+1,i))THEN
265 p1 = dd(k)
266 p2 = k
267 ENDIF
268C
269 IF (edgel(k,i) == ied1) GOTO 140
270 IF (xxl(p1,i) == xxl(p2,i)) GOTO 140 ! no inters (parallel to dir 2)
271 m12 = xxl(p2,i)-xxl(p1,i)
272 m12 = (yyl(p2,i)-yyl(p1,i))/m12
273 xint = xint0
274 yint = yyl(p1,i)+m12*(xint-xxl(p1,i))
275 cross12 = (xint-xxl(p1,i))*(xint-xxl(p2,i))+
276 . (yint-yyl(p1,i))*(yint-yyl(p2,i))
277 IF (cross12 > zero) GOTO 140
278c
279 cross1 = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
280 beta = sqrt(cross1 / len(k,i))
281 beta = max(beta, bmin)
282 beta = min(beta, bmax)
283 beta0(k,i) = beta
284C
285 xin(ied2,i) = xint
286 yin(ied2,i) = yint
287 edgel(k,i) = ied2
288 EXIT
289 140 CONTINUE
290 ELSEIF(dir22 == zero)THEN
291 DO 150 k=1,3
292 xint = zero
293 yint = zero
294 elcrktg = iel_crktg(i+nft)
295 elcrk = elcrktg + ecrkxfec
296 iedge = xedge3n(k,elcrktg)
297 nod1 = nodedge(1,iedge)
298 nod2 = nodedge(2,iedge)
299 IF(nod1 == ixtg(k+1,i) .and. nod2 == ixtg(dd(k)+1,i))THEN
300 p1 = k
301 p2 = dd(k)
302 ELSE IF(nod2 == ixtg(k+1,i).and.nod1==ixtg(dd(k)+1,i))THEN
303 p1 = dd(k)
304 p2 = k
305 ENDIF
306C
307 IF (edgel(k,i) == ied1) GOTO 150
308 IF (yyl(p1,i) == yyl(p2,i)) GOTO 150 ! no inters (parallel to dir 2)
309 m12 = yyl(p2,i)-yyl(p1,i)
310 m12 = (xxl(p2,i)-xxl(p1,i))/m12
311 yint = yint0
312 xint = xxl(p1,i)+m12*(yint-yyl(p1,i))
313 cross12 = (xint-xxl(p1,i))*(xint-xxl(p2,i))+
314 . (yint-yyl(p1,i))*(yint-yyl(p2,i))
315 IF (cross12 > zero) GOTO 150
316C
317 cross1 = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
318 beta = sqrt(cross1 / len(k,i))
319 beta = max(beta, bmin)
320 beta = min(beta, bmax)
321 beta0(k,i) = beta
322C
323 xin(ied2,i) = xint
324 yin(ied2,i) = yint
325 edgel(k,i) = ied2
326 EXIT
327 150 CONTINUE
328 ELSEIF(dir11 /= zero .AND. dir22 /= zero)THEN
329 DO 160 k=1,3
330 xint = zero
331 yint = zero
332 elcrktg = iel_crktg(i+nft)
333 elcrk = elcrktg + ecrkxfec
334 iedge = xedge3n(k,elcrktg)
335 nod1 = nodedge(1,iedge)
336 nod2 = nodedge(2,iedge)
337 IF (nod1 == ixtg(k+1,i) .and. nod2 == ixtg(dd(k)+1,i))THEN
338 p1 = k
339 p2 = dd(k)
340 ELSE IF (nod2 == ixtg(k+1,i).and.nod1==ixtg(dd(k)+1,i))THEN
341 p1 = dd(k)
342 p2 = k
343 ENDIF
344C
345 IF (edgel(k,i) == ied1) GOTO 160
346 IF (xxl(p1,i) == xxl(p2,i)) THEN
347 mm = dir22/dir11
348 xint = xxl(p1,i) ! or = XXL(p2,I)
349 yint = yint0+mm*(xint-xint0)
350 cross12 = (xint-xxl(p1,i))*(xint-xxl(p2,i))+
351 . (yint-yyl(p1,i))*(yint-yyl(p2,i))
352 IF (cross12 > zero) GOTO 160
353C
354 cross1 = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
355 beta = sqrt(cross1 / len(k,i))
356 beta = max(beta, bmin)
357 beta = min(beta, bmax)
358 beta0(k,i) = beta
359C
360 xin(ied2,i) = xint
361 yin(ied2,i) = yint
362 edgel(k,i) = ied2
363 EXIT
364 ELSE
365 mm = dir22/dir11
366 m12 = xxl(p2,i)-xxl(p1,i)
367 m12 = (yyl(p2,i)-yyl(p1,i))/m12
368 IF (mm == m12) GOTO 160 ! no inters (parallel to dir 2)
369 xint = (yint0-yyl(p1,i)+m12*xxl(p1,i)-mm*xint0)/(m12-mm)
370 yint = yint0+mm*(xint-xint0)
371 cross12 = (xint-xxl(p1,i))*(xint-xxl(p2,i))+
372 . (yint-yyl(p1,i))*(yint-yyl(p2,i))
373 IF (cross12 > zero) GOTO 160
374C
375 cross1 = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
376 beta = sqrt(cross1 / len(k,i))
377 beta = max(beta, bmin)
378 beta = min(beta, bmax)
379 beta0(k,i) = beta
380C
381 xin(ied2,i) = xint
382 yin(ied2,i) = yint
383 edgel(k,i) = ied2
384 EXIT
385 ENDIF
386 160 CONTINUE
387 ENDIF
388 130 CONTINUE
389 ENDDO
390C----------------------------------------------------------------------
391C check for getting both intersections
392C
393 DO ir=1,newcrk
394 i = jct(ir)
395 fac = 0
396 DO k=1,3
397 IF (edgel(k,i)==1 .or. edgel(k,i)==2) fac=fac+1
398 ENDDO
399 IF (fac /= 2) THEN
400 WRITE(iout,*) 'ERROR IN ADVANCING CRACK.NO CUT EDGES'
401 CALL arret(2)
402 ENDIF
403 ENDDO
404C---
405 DO ir=1,newcrk
406 i = jct(ir)
407 elcrktg = iel_crktg(i+nft)
408 DO k=1,3
409 ied = edgel(k,i)
410 IF (ied > 0) THEN
411 crkedge(ilay)%IEDGETG(k,elcrktg) = ied ! tag cut edges of each phantom
412 ENDIF
413 ENDDO
414 ENDDO
415C----------------------------------------------------------------------
416C SIGN DISTANCE FOR NEW CRACKED LAYER
417C----------------------------------------------------------------------
418 DO ir=1,newcrk
419 i = jct(ir)
420 fit(1,i) = zero
421 fit(2,i) = zero
422 fit(3,i) = zero
423 xn(1) = xl1(i)
424 yn(1) = yl1(i)
425 xn(2) = xl2(i)
426 yn(2) = yl2(i)
427 xn(3) = xl3(i)
428 yn(3) = yl3(i)
429c
430 DO k=1,3
431 p1 = k
432 p2 = dd(k)
433 ied = edgel(k,i)
434 IF (ied > 0) THEN
435 xmi(ied) = half*(xn(p1)+xn(p2))
436 ymi(ied) = half*(yn(p1)+yn(p2))
437 ENDIF
438 ENDDO
439 DO k=1,3
440 fi=zero
441 CALL lsint4(xmi(1),ymi(1),xmi(2),ymi(2),xn(k),yn(k),fi )
442 IF (fit(k,i) == zero) fit(k,i)=fi
443 ENDDO
444 ENDDO
445C
446 DO ir=1,newcrk
447 i = jct(ir)
448 elcrktg = iel_crktg(i+nft)
449C
450 DO k=1,3
451 ied = edgel(k,i)
452 IF (ied == 2) THEN ! only for second intersection (tip)
453 iedge = xedge3n(k,elcrktg)
454c
455 icut = crkedge(ilay)%ICUTEDGE(iedge)
456 IF (icut > 0) THEN ! already cut before => edge connecting two cracks
457 crkedge(ilay)%ICUTEDGE(iedge) = 3
458 ELSE ! new edge cut
459 crkedge(ilay)%ICUTEDGE(iedge) = 2
460 crkedge(ilay)%RATIO(iedge) = beta0(k,i)
461 ENDIF
462 ENDIF
463 ENDDO
464 ENDDO
465C-----------------------
466c Fill New Cut Phantom Element
467C-----------------------
468 DO ir=1,newcrk
469 i = jct(ir)
470 elcrktg = iel_crktg(i+nft)
471 elcrk = elcrktg + ecrkxfec
472C
473 iadc(1) = iad_crktg(1,elcrktg)
474 iadc(2) = iad_crktg(2,elcrktg)
475 iadc(3) = iad_crktg(3,elcrktg)
476C
477 n(1) = ixtg(2,i)
478 n(2) = ixtg(3,i)
479 n(3) = ixtg(4,i)
480C
481 nn(1) = inod_crk(n(1))
482 nn(2) = inod_crk(n(2))
483 nn(3) = inod_crk(n(3))
484C
485 icrk = elfiss(i)
486C
487 elcutc(1,i) = 2
488 numelcrk = numelcrk + 1
489C
490 isign(1) = int(sign(one,fit(1,i)))
491 isign(2) = int(sign(one,fit(2,i)))
492 isign(3) = int(sign(one,fit(3,i)))
493C
494 IF (fit(1,i) == zero) isign(1) = 0
495 IF (fit(2,i) == zero) isign(2) = 0
496 IF (fit(3,i) == zero) isign(3) = 0
497c
498 itri = 0
499 nm = 0
500 np = 0
501 DO k=1,3
502 IF (isign(k) > 0) THEN
503 itri = itri + 1
504 np = k
505 ELSEIF (isign(k) < 0) THEN
506 nm = k
507 ENDIF
508 ENDDO
509 IF (itri == 1) THEN
510 itri = -1
511 nx1 = np
512 ELSEIF (itri == 2) THEN
513 itri = 1
514 nx1 = nm
515 ENDIF
516 nx2 = dx(nx1+1)
517 nx3 = dx(nx2+1)
518 xfem_phantom(ilay)%ITRI(1,elcrk) = itri
519 xfem_phantom(ilay)%ITRI(2,elcrk) = nx1
520c--------------------------------------------
521c activate enriched nodes
522c--------------------------------------------
523 DO k=1,3
524 IF(ienr0(k,i) /= 0)THEN
525 ienr(k) = ienr0(k,i)
526 ELSE
527 ienr(k) = crknodiad(iadc(k)) + knod2elc(nn(k))*(ilay-1)
528 ENDIF
529 ENDDO
530C
531 sigbeta = 0
532 DO k=1,3
533 ied = edgel(k,i)
534 IF (ied == 2) THEN ! only new cut edge
535 iedge = xedge3n(k,elcrktg)
536 nod1 = nodedge(1,iedge)
537 nod2 = nodedge(2,iedge)
538 ie10 = crkedge(ilay)%EDGEENR(1,iedge)
539 ie20 = crkedge(ilay)%EDGEENR(2,iedge)
540 IF (nod1 == ixtg(k+1,i) .and. nod2 == ixtg(dd(k)+1,i)) THEN
541 ie1 = ienr(k)
542 ie2 = ienr(dd(k))
543 ifi1 = isign(k)
544 ifi2 = isign(dd(k))
545 sigbeta = iedge
546 ELSE IF (nod2 == ixtg(k+1,i).and.nod1 == ixtg(dd(k)+1,i)) THEN
547 ie1 = ienr(dd(k))
548 ie2 = ienr(k)
549 ifi1 = isign(dd(k))
550 ifi2 = isign(k)
551 sigbeta = -iedge
552 END IF
553 crkedge(ilay)%EDGEENR(1,iedge) = max(ie1,ie10)
554 crkedge(ilay)%EDGEENR(2,iedge) = max(ie2,ie20)
555 IF (crkedge(ilay)%EDGEICRK(iedge) == 0)
556 . crkedge(ilay)%EDGEICRK(iedge) = icrk
557 ENDIF ! IF(IED == 2)THEN
558 ENDDO
559C--------------------------------------------
560 crkedge(ilay)%LAYCUT(elcrk) = 1
561 xfem_phantom(ilay)%ELCUT(elcrk) = icrk
562C--------------------------------------------
563 DO k=1,3
564 ied = edgel(k,i)
565 iedge = xedge3n(k,elcrktg)
566 IF (ied > 0) THEN
567 crkedge(ilay)%EDGETIP(1,iedge) = tip(i)
568 crkedge(ilay)%EDGETIP(2,iedge) =
569 . crkedge(ilay)%EDGETIP(2,iedge) + 1
570 ENDIF
571 ENDDO
572C
573 xfem_phantom(ilay)%IFI(iadc(1)) = isign(1)
574 xfem_phantom(ilay)%IFI(iadc(2)) = isign(2)
575 xfem_phantom(ilay)%IFI(iadc(3)) = isign(3)
576C------------------
577c IXEL = 1 => positif element within ILAY (ILEV = PP1)
578C------------------
579C
580 crklvset(pp1)%ENR0(1,iadc(1)) = -ienr(1)
581 crklvset(pp1)%ENR0(1,iadc(2)) = -ienr(2)
582 crklvset(pp1)%ENR0(1,iadc(3)) = -ienr(3)
583C
584 IF (isign(1) > 0) crklvset(pp1)%ENR0(1,iadc(1)) = 0
585 IF (isign(2) > 0) crklvset(pp1)%ENR0(1,iadc(2)) = 0
586 IF (isign(3) > 0) crklvset(pp1)%ENR0(1,iadc(3)) = 0
587C------------------
588c IXEL = 2 => negatif element within ILAY (ILEV = PP2)
589C------------------
590 crklvset(pp2)%ENR0(1,iadc(1)) = -ienr(1)
591 crklvset(pp2)%ENR0(1,iadc(2)) = -ienr(2)
592 crklvset(pp2)%ENR0(1,iadc(3)) = -ienr(3)
593C
594 IF (isign(1) < 0) crklvset(pp2)%ENR0(1,iadc(1)) = 0
595 IF (isign(2) < 0) crklvset(pp2)%ENR0(1,iadc(2)) = 0
596 IF (isign(3) < 0) crklvset(pp2)%ENR0(1,iadc(3)) = 0
597C------------------
598c IXEL = 3 => Third phantom element (ILEV = PP3)
599C------------------
600 IF (itri < 0) THEN ! sign ILEV3 = ILEV2 < 0
601 ie2 = xedge3n(nx3,elcrktg) ! tip edge
602 IF (crkedge(ilay)%ICUTEDGE(ie2) > 1) THEN
603 sigbeta = -sigbeta
604 iad1 = iadc(nx1)
605 iad2 = iadc(nx2)
606 iad3 = iadc(nx3)
607 nod1 = iad3
608 nod2 = iad1
609 crklvset(pp3)%ENR0(1,iad1) = abs(crklvset(pp2)%ENR0(1,iad1))
610 crklvset(pp3)%ENR0(1,iad2) = crklvset(pp2)%ENR0(1,iad2)
611 crklvset(pp3)%ENR0(1,iad3) = crklvset(pp2)%ENR0(1,iad3)
612 crklvset(pp2)%ENR0(1,iad1) = -crknodiad(iad1) - knod2elc(nn(nx1))*(ilay-1)
613c
614c-- AREA factors for each phantom
615 x1 = xxl(nx1,i)
616 y1 = yyl(nx1,i)
617 ied = crkedge(ilay)%IEDGETG(nx1,elcrktg)
618 x2 = xin(ied,i)
619 y2 = yin(ied,i)
620 ied = crkedge(ilay)%IEDGETG(nx3,elcrktg)
621 x3 = xin(ied,i)
622 y3 = yin(ied,i)
623 area1 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
624 area1 = area1 / area(i)
625 x1 = xxl(nx2,i)
626 y1 = yyl(nx2,i)
627 x2 = xxl(nx3,i)
628 y2 = yyl(nx3,i)
629 area2 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
630 area2 = area2 / area(i)
631 area3 = one - area1 - area2
632c
633 ELSE
634c print*,' IMPOSSIBLE CASE'
635 ENDIF
636c
637 ELSEIF (itri > 0) THEN ! sign ILEV3 = ILEV1 > 0
638c
639 ie1 = xedge3n(nx1,elcrktg) ! tip edge
640 IF (crkedge(ilay)%ICUTEDGE(ie1) > 1) THEN
641 iad1 = iadc(nx1)
642 iad2 = iadc(nx2)
643 iad3 = iadc(nx3)
644 crklvset(pp3)%ENR0(1,iad1) = abs(crklvset(pp1)%ENR0(1,iad1))
645 crklvset(pp3)%ENR0(1,iad2) = crklvset(pp1)%ENR0(1,iad2)
646 crklvset(pp3)%ENR0(1,iad3) = crklvset(pp1)%ENR0(1,iad3)
647 crklvset(pp1)%ENR0(1,iad1) = -crknodiad(iad1) - knod2elc(nn(nx1))*(ilay-1)
648c
649c-- AREA factors for each phantom
650 ied1= crkedge(ilay)%IEDGETG(nx1,elcrktg)
651 ied2= crkedge(ilay)%IEDGETG(nx3,elcrktg)
652 x1 = xin(ied1,i)
653 y1 = yin(ied1,i)
654 x2 = xxl(nx2,i)
655 y2 = yyl(nx2,i)
656 x3 = xxl(nx3,i)
657 y3 = yyl(nx3,i)
658 area1 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
659 area1 = area1 / area(i)
660 x1 = xxl(nx1,i)
661 y1 = yyl(nx1,i)
662 x2 = xin(ied1,i)
663 y2 = yin(ied1,i)
664 x3 = xin(ied2,i)
665 y3 = yin(ied2,i)
666 area2 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
667 area2 = area2 / area(i)
668 area3 = one - area1 - area2
669 ELSE
670c print*,' IMPOSSIBLE CASE'
671 ENDIF
672 ENDIF ! ITRI
673c----
674 crklvset(pp1)%AREA(elcrk) = area1
675 crklvset(pp2)%AREA(elcrk) = area2
676 crklvset(pp3)%AREA(elcrk) = area3
677c
678 IF (area3 < zero .or. area1 > one .or. area2 > one .or. area3 > one ) THEN
679 print*,'ERROR : XFEM PHANTOM ELEMENT AREA: ELCRK=',elcrk
680 ENDIF
681c
682 ENDDO ! IR=1,NEWCRK
683C-----------
684 RETURN
685 END
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
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