OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
crklayer3n_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 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)

Function/Subroutine Documentation

◆ crklayer3n_adv()

subroutine crklayer3n_adv ( integer nel,
integer nft,
integer ilay,
integer nlay,
integer, dimension(nixtg,*) ixtg,
integer, dimension(2,*) elcutc,
integer, dimension(nlay,*) elcrkini,
integer, dimension(*) iel_crktg,
integer, dimension(*) inod_crk,
integer, dimension(3,*) iad_crktg,
integer, dimension(*) nodenr,
dir1,
dir2,
integer, dimension(2,*) nodedge,
integer, dimension(*) crknodiad,
integer, dimension(*) knod2elc,
type (xfem_edge_), dimension(*) crkedge,
integer, dimension(3,*) xedge3n,
integer, dimension(nel) ngl,
dimension(nel) area,
dimension(nel) xl2,
dimension(nel) xl3,
dimension(nel) yl2,
dimension(nel) yl3 )

Definition at line 33 of file crklayer3n_adv.F.

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 IF (iedge > 0 .and. edgel(k,i) == 1) THEN
212 icut = crkedge(ilay)%ICUTEDGE(iedge)
213 IF (icut == 1) THEN
214 beta = crkedge(ilay)%RATIO(iedge)
215C
216 IF (beta > one .or. 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
my_real
#define my_real
Definition cppsort.cpp:32
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