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