42
43
44
45
46
47 USE elbufdef_mod
49
50
51
52#include "implicit_f.inc"
53
54
55
56#include "units_c.inc"
57#include "com04_c.inc"
58#include "com_xfem1.inc"
59#include "comlock.inc"
60
61
62
63 INTEGER NEL,NFT,ILAY,IXC(NIXC,*),NLAY,NGL(NEL),IEL_CRK(*),
64 . INOD_CRK(*),NODENR(*),IADC_CRK(4,*),ELCRKINI(NLAY,NEL),
65 . ELCUTC(2,*),NODEDGE(2,*),CRKNODIAD(*),KNOD2ELC(*),XEDGE4N(4,*)
66C
67 my_real,
DIMENSION(NLAY,NEL) :: dir1,dir2
69 TYPE (ELBUF_STRUCT_), DIMENSION(NXEL), TARGET :: XFEM_STR
70 TYPE (XFEM_EDGE_) , DIMENSION(NXLAYMAX) :: CRKEDGE
71
72
73
74 INTEGER I,J,K,R,II,IR,P1,P2,PP1,PP2,PP3,NEWCRK,LAYCUT,FAC,
75 . IEDGE,ICUT,IED,IED1,IED2,ICRK,ELCRK,NOD1,NOD2,IE1,IE2,ITRI,
76 . IFI1,IFI2,IE10,IE20,NP,NX1,NX2,NX3,NX4
77 INTEGER N(4),NN(4),dd(4),d(8),DX(8),
78 . ISIGN(4),IENR0(4),IENR(4),INV(2),ECUT(2,NEL),
79 . IADC(4),JCT(NEL),EDGEL(4,NEL),KPERM(8)
80C
82 . xm(nel),ym(nel),xl2(nel),yl2(nel),xl3(nel),yl3(nel),
83 . xl4(nel),yl4(nel),xin(2,nel),yin(2,nel),xmi(2),ymi(2),
84 . fit(4,nel),xxl(4,nel),yyl(4,nel),len(4,nel),beta0(4,nel)
85
87 . xint,yint,fi,xxx,yyy,zzz,d12,m12,mm,cross,acd,bcd,dlx,dly,
88 . xint0,yint0,dir11,dir22,beta,bmin,bmax,
89 . x1,y1,x2,y2,x3,y3,x4,y4,area1,area2
90
91 DATA d/1,2,2,3,4,3,1,4/
92 DATA dd/2,3,4,1/
93 DATA dx/1,2,3,4,1,2,3,4/
94 DATA kperm/1,2,3,4,1,2,3,4/
95 parameter(bmin = 0.01, bmax = 0.99)
96
97 newcrk = 0
98 DO i=1,nel
99 jct(i) = 0
100 IF (elcrkini(ilay,i) == -1) THEN
101 newcrk = newcrk + 1
102 jct(newcrk) = i
103 ENDIF
104 ENDDO
105 IF (newcrk == 0) RETURN
106
107 ii = nxel*(ilay-1)
108 pp1 = ii + 1
109 pp2 = ii + 2
110 pp3 = ii + 3
111
112 DO ir=1,newcrk
113 i = jct(ir)
114 edgel(1:4,i) = 0
115 ecut(1:2,i) = 0
116 beta0(1:4,i) = zero
117 xin(1,i)=zero
118 yin(1,i)=zero
119 xin(2,i)=zero
120 yin(2,i)=zero
121
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)
130 xm(i) = fourth*(xl2(i)+xl3(i)+xl4(i))
131 ym(i) = fourth*(yl2(i)+yl3(i)+yl4(i))
132
133 len(1,i) = xl2(i)*xl2(i) + yl2(i)*yl2(i)
134 len(2,i) = (xl3(i)-xl2(i))*(xl3(i)-xl2(i))+
135 . (yl3(i)-yl2(i))*(yl3(i)-yl2(i))
136 len(3,i) = (xl4(i)-xl3(i))*(xl4(i)-xl3(i))+
137 . (yl4(i)-yl3(i))*(yl4(i)-yl3(i))
138 len(4,i) = xl4(i)*xl4(i) + yl4(i)*yl4(i)
139 END DO
140
141
142
143 DO ir=1,newcrk
144 i = jct(ir)
145 elcrk = iel_crk(i+nft)
146
147 dir11 = -dir2(ilay,i)
148 dir22 = dir1(ilay,i)
149 fac = 0
150
151 DO k=1,4
152 iedge = xedge4n(k,elcrk)
153 nod1 = nodedge(1,iedge)
154 nod2 = nodedge(2,iedge)
155 IF (nod1 == ixc(k+1,i) .and. nod2 == ixc(dd(k)+1,i)) THEN
156 p1 = k
157 p2 = dd(k)
158 ELSEIF (nod2 == ixc(k+1,i).and.nod1 == ixc(dd(k)+1,i)) THEN
159 p1 = dd(k)
160 p2 = k
161 ENDIF
162
163 IF (dir11 == zero) THEN
164 dlx = xxl(p2,i) - xxl(p1,i)
165 IF (dlx /= zero) THEN
166 dly = yyl(p2,i) - yyl(p1,i)
167 m12 = dly / dlx
168 xint = xm(i)
169 yint = yyl(p1,i) + m12*(xint-xxl(p1,i))
170 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero .and.
171 . (yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
172 fac = 1
173 EXIT
174 ENDIF
175 ENDIF
176
177 ELSEIF (dir22 == zero) THEN
178 dly = yyl(p2,i) - yyl(p1,i)
179 IF (dly /= zero) THEN
180 dlx = xxl(p2,i) - xxl(p1,i)
181 m12 = dlx / dly
182 yint = ym(i)
183 xint = xxl(p1,i) + m12*(yint-yyl(p1,i))
184 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero .and.
185 . (yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
186 fac = 1
187 EXIT
188 ENDIF
189 ENDIF
190
191 ELSEIF (dir11 /= zero .AND. dir22 /= zero) THEN
192 dlx = xxl(p2,i) - xxl(p1,i)
193 dly = yyl(p2,i) - yyl(p1,i)
194 mm = dir22/dir11
195 IF (dlx == zero) THEN
196 xint = xxl(p1,i)
197 yint = ym(i) + mm*(xint-xm(i))
198 IF ((yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
199 fac = 1
200 ENDIF
201 ELSEIF (dly == zero) THEN
202 yint = yyl(p1,i)
203 xint = xm(i) + (ym(i)-yyl(p1,i)) / mm
204 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero) THEN
205 fac = fac + 1
206 EXIT
207 ENDIF
208 ELSE
209 m12 = dly / dlx
210 IF (mm /= m12) THEN
211 xint = (ym(i)-yyl(p1,i) + m12*xxl(p1,i) - mm*xm(i))/(m12-mm)
212 yint = ym(i) + mm*(xint-xm(i))
213 acd = (yint-yyl(p1,i))*(xm(i) - xxl(p1,i))
214 . - (xint-xxl(p1,i))*(ym(i) - yyl(p1,i))
215 bcd = (yint-yyl(p2,i))*(xm(i) - xxl(p2,i))
216 . - (xint-xxl(p2,i))*(ym(i) - yyl(p2,i))
217 IF (acd*bcd <= em3) THEN
218 fac = 1
219 EXIT
220 ENDIF
221 ENDIF
222 ENDIF
223 ENDIF
224
225 ENDDO
226
227 IF (fac == 1) THEN
228 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
229 beta = sqrt(cross / len(k,i))
230 IF (beta > bmax .OR. beta < bmin) THEN
231 beta =
max(beta, bmin)
232 beta =
min(beta, bmax)
233 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
234 yint = yyl(p1,i) + beta*(yyl
235 ENDIF
236
237 ecut(fac,i)= k
238 xin(fac,i) = xint
239 yin(fac,i) = yint
240 edgel(k,i) = fac
241 beta0(k,i) = beta
242 ELSE
243 WRITE(iout,*) 'ERROR IN ADVANCING CRACK --- CHECK CRACK TIP'
245 ENDIF
246
247 ENDDO
248
249
250
251 DO ir=1,newcrk
252 i = jct(ir)
253 elcrk = iel_crk(i+nft)
254 xint0 = xin(1,i)
255 yint0 = yin(1,i)
256 dir11 =-dir2(ilay,i)
257 dir22 = dir1(ilay,i)
258
259 k = ecut(1,i)
260 r = kperm(k + 2)
261 iedge = xedge4n(r,elcrk)
262 nod1 = nodedge(1,iedge)
263 nod2 = nodedge(2,iedge)
264 IF (nod1 == ixc(r+1,i) .and. nod2 == ixc(dd(r)+1,i))THEN
265 p1 = r
266 p2 = dd(r)
267 ELSE IF (nod2 == ixc(r+1,i).and.nod1 == ixc(dd(r)+1,i))THEN
268 p1 = dd(r)
269 p2 = r
270 ENDIF
271
272 IF (dir11 == zero) THEN
273 dlx = xxl(p2,i) - xxl(p1,i)
274 IF (dlx /= zero) THEN
275 dly = yyl(p2,i) - yyl(p1,i)
276 m12 = dly / dlx
277 xint = xm(i)
278 yint = yyl(p1,i) + m12*(xint-xxl(p1,i))
279 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero .and.
280 . (yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
281 fac = 2
282 ENDIF
283 ENDIF
284
285 ELSEIF (dir22 == zero) THEN
286 dly = yyl(p2,i) - yyl(p1,i)
287 IF (dly /= zero) THEN
288 dlx = xxl(p2,i) - xxl(p1,i)
289 m12 = dlx / dly
290 yint = ym(i)
291 xint = xxl(p1,i) + m12*(yint-yyl(p1,i))
292 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero .and.
293 . (yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
294 fac = 2
295 ENDIF
296 ENDIF
297
298 ELSEIF (dir11 /= zero .AND. dir22 /= zero) THEN
299 dlx = xxl(p2,i) - xxl(p1,i)
300 dly = yyl(p2,i) - yyl(p1,i)
301 mm = dir22/dir11
302 IF (dlx == zero) THEN
303 xint = xxl(p1,i)
304 yint = ym(i) + mm*(xint-xm(i))
305 IF ((yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
306 fac = 2
307 ENDIF
308 ELSEIF (dly == zero) THEN
309 yint = yyl(p1,i)
310 xint = xm(i) + (ym(i)-yyl(p1,i)) / mm
311 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero) THEN
312 fac = 2
313 ENDIF
314 ELSE
315 m12 = dly / dlx
316 IF (mm /= m12) THEN
317 xint = (ym(i)-yyl(p1,i) + m12*xxl(p1,i) - mm*xm(i))/(m12-mm)
318 yint = ym(i) + mm*(xint-xm(i))
319 acd = (yint-yyl(p1,i))*(xm(i) - xxl(p1,i))
320 . - (xint-xxl(p1,i))*(ym(i) - yyl(p1,i))
321 bcd = (yint-yyl(p2,i))*(xm(i) - xxl(p2,i))
322 . - (xint-xxl(p2,i))*(ym(i) - yyl(p2,i))
323
324 fac = 2
325
326 ENDIF
327 ENDIF
328 ENDIF
329
330 IF (fac == 2) THEN
331 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
332 beta = sqrt(cross / len(r,i))
333 IF (beta > bmax .OR. beta < bmin) THEN
334 beta =
max(beta, bmin)
335 beta =
min(beta, bmax)
336 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
337 yint = yyl(p1,i) + beta*(yyl(p2,i)-yyl(p1,i))
338 ENDIF
339
340 ecut(2,i) = r
341 xin(2,i) = xint
342 yin(2,i) = yint
343 edgel(r,i)= 2
344 beta0(r,i)= beta
345 ENDIF
346 ENDDO
347
348
349
350 DO ir=1,newcrk
351 i = jct(ir)
352 fac = 0
353 DO j=1,2
354 k = ecut(j,i)
355 IF (edgel(k,i)==1 .or. edgel(k,i)==2) fac=fac+1
356 ENDDO
357 IF (fac /= 2) THEN
358 WRITE(iout,*) 'ERROR IN INITIATION CRACK.NO CUT EDGES'
360 ENDIF
361 ENDDO
362
363
364
365 DO ir=1,newcrk
366 i = jct(ir)
367 elcrk = iel_crk(i+nft)
368 DO j=1,2
369 k = ecut(j,i)
370 crkedge(ilay)%IEDGEC(k,elcrk) = edgel(k,i)
371 ENDDO
372 ENDDO
373
374
375
376 DO ir=1,newcrk
377 i = jct(ir)
378 fit(1,i)=zero
379 fit(2,i)=zero
380 fit(3,i)=zero
381 fit(4,i)=zero
382 DO k=1,4
383 p1 = k
384 p2 = dd(k)
385 ied = edgel(k,i)
386 IF (ied > 0) THEN
387 xmi(ied) = half*(xxl(p1,i)+xxl(p2,i))
388 ymi(ied) = half*(yyl(p1,i)+yyl(p2,i))
389 ENDIF
390 ENDDO
391
392 DO k=1,4
393 fi = zero
394 CALL lsint4(xmi(1),ymi(1),xmi(2),ymi(2),xxl(k,i),yyl(k,i),fi )
395 IF (fit(k,i)==zero) fit(k,i) = fi
396 ENDDO
397 ENDDO
398
399 DO ir=1,newcrk
400 i = jct(ir)
401 elcrk = iel_crk(i+nft)
402 DO j=1,2
403 k = ecut(j,i)
404 iedge = xedge4n(k,elcrk)
405 icut = crkedge(ilay)%ICUTEDGE(iedge)
406 IF (icut > 0) THEN
407 crkedge(ilay)%ICUTEDGE(iedge) = 3
408 ELSE
409 crkedge(ilay)%ICUTEDGE(iedge) = 2
410 crkedge(ilay)%RATIO(iedge) = beta0(k,i)
411 ENDIF
412 ENDDO
413 ENDDO
414
415
416
417 DO ir=1,newcrk
418 i = jct(ir)
419 elcrk = iel_crk(i+nft)
420 elcutc(1,i) = 2
421 numelcrk = numelcrk + 1
422
423 iadc(1) = iadc_crk(1,elcrk)
424 iadc(2) = iadc_crk(2,elcrk)
425 iadc(3) = iadc_crk(3,elcrk)
426 iadc(4) = iadc_crk(4,elcrk)
427
428 n(1) = ixc(2,i)
429 n(2) = ixc(3,i)
430 n(3) = ixc(4,i)
431 n(4) = ixc(5,i)
432 nn(1) = inod_crk(n(1))
433 nn(2) = inod_crk(n(2))
434 nn(3) = inod_crk(n(3))
435 nn(4) = inod_crk(n(4))
436
437 isign(1) = int(sign(one,fit(1,i)))
438 isign(2) = int(sign(one,fit(2,i)))
439 isign(3) = int(sign(one,fit(3,i)))
440 isign(4) = int(sign(one,fit(4,i)))
441
442 IF (fit(1,i) == zero) isign(1) = 0
443 IF (fit(2,i) == zero) isign(2) = 0
444 IF (fit(3,i) == zero) isign(3) = 0
445 IF (fit(4,i) == zero) isign(4) = 0
446
447 icrk =
crkshell(pp1)%CRKSHELLID(elcrk)
448
449 ienr0(1) = crknodiad(iadc(1))
450 ienr0(2) = crknodiad(iadc(2))
451 ienr0(3) = crknodiad(iadc(3))
452 ienr0(4) = crknodiad(iadc(4))
453
454 ienr(1) = ienr0(1) + knod2elc(nn(1))*(ilay-1)
455 ienr(2) = ienr0(2) + knod2elc(nn(2))*(ilay-1)
456 ienr(3) = ienr0(3) + knod2elc(nn(3))*(ilay-1)
457 ienr(4) = ienr0(4) + knod2elc(nn(4))*(ilay-1)
458
459 DO j=1,2
460 k = ecut(j,i)
461 iedge = xedge4n(k,elcrk)
462 nod1 = nodedge(1,iedge)
463 nod2 = nodedge(2,iedge)
464 ie10 = crkedge(ilay)%EDGEENR(1,iedge)
465 ie20 = crkedge(ilay)%EDGEENR(2,iedge)
466 IF (nod1 == n(k) .and. nod2 == n(dd(k))) THEN
467 ie1 = ienr(k)
468 ie2 = ienr(dd(k))
469 ifi1 = isign(k)
470 ifi2 = isign(dd(k))
471 ELSE IF (nod2 == n(k) .and. nod1 == n(dd(k))) THEN
472 ie1 = ienr(dd(k))
473 ie2 = ienr(k)
474 ifi1 = isign(dd(k))
475 ifi2 = isign(k)
476 END IF
477 crkedge(ilay)%EDGEENR(1,iedge) =
max(ie1,ie10)
478 crkedge(ilay)%EDGEENR(2,iedge) =
max(ie2,ie20)
479 IF (crkedge(ilay)%EDGEICRK(iedge) == 0)
480 . crkedge(ilay)%EDGEICRK(iedge)
481 ENDDO
482
483 crkedge(ilay)%LAYCUT(elcrk) = -1
485
486 np = 0
487 DO k=1,4
488 ied = edgel(k,i)
489 iedge = xedge4n(k,elcrk)
490 IF (ied > 0) THEN
491 crkedge(ilay)%EDGETIP(1,iedge) = ied
492 crkedge(ilay)%EDGETIP(2,iedge) =
493 . crkedge(ilay)%EDGETIP(2,iedge) + 1
494 ENDIF
495 IF (isign(k) > 0) np = k
496 ENDDO
497
498 itri = 0
499 nx1 = np
500 IF (np > 0 .and. isign(np-1) > 0) THEN
501 nx1 = np-1
502 ELSE
503 nx1 = np
504 ENDIF
507 nx2 = dx(nx1+1)
508 nx3 = dx(nx1+2)
509 nx4 = dx(nx1+3)
510
511 x1 = xxl(nx1,i)
512 y1 = yyl(nx1
513 x2 = xxl(nx2,i)
514 y2 = yyl(nx2,i)
515 ied = edgel(nx2,i)
516 IF (ied > 0) THEN
517 x3 = xin(ied,i)
518 y3 = yin(ied,i)
519 ELSE
520 print*,' ERROR: K,IED=',k,ied
521 ENDIF
522 ied = edgel(nx4,i)
523 IF (ied > 0) THEN
524 x4 = xin(ied,i)
525 y4 = yin(ied,i)
526 ELSE
527 print*,' ERROR: K,IED=',k,ied
528 ENDIF
529 area1 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
530 area1 = area1 * a_i(i)
531 area2 = one - area1
535
540
541
542
543 crklvset(pp1)%ENR0(1,iadc(1)) = -ienr(1)
544 crklvset(pp1)%ENR0(1,iadc(2)) = -ienr(2)
545 crklvset(pp1)%ENR0(1,iadc(3)) = -ienr(3)
546 crklvset(pp1)%ENR0(1,iadc(4)) = -ienr(4)
547
548 IF (isign(1) > 0)
crklvset(pp1)%ENR0(1,iadc(1)) = 0
549 IF (isign(2) > 0)
crklvset(pp1)%ENR0(1,iadc(2)) = 0
550 IF (isign(3) > 0)
crklvset(pp1)%ENR0(1,iadc(3)) = 0
551 IF (isign(4) > 0)
crklvset(pp1)%ENR0(1,iadc(4)) = 0
552
553
554
555 crklvset(pp2)%ENR0(1,iadc(1)) = -ienr(1)
556 crklvset(pp2)%ENR0(1,iadc(2)) = -ienr(2)
557 crklvset(pp2)%ENR0(1,iadc(3)) = -ienr(3)
558 crklvset(pp2)%ENR0(1,iadc(4)) = -ienr(4)
559
560 IF (isign(1) < 0)
crklvset(pp2)%ENR0(1,iadc(1)) = 0
561 IF (isign(2) < 0)
crklvset(pp2)%ENR0(1,iadc(2)) = 0
562 IF (isign(3) < 0)
crklvset(pp2)%ENR0(1,iadc(3)) = 0
563 IF (isign(4) < 0)
crklvset(pp2)%ENR0(1,iadc(4)) = 0
564
565
566
567 xfem_str(nxel)%GBUF%OFF(i) = zero
568 xfem_str(nxel)%BUFLY(ilay)%LBUF(1,1,1)%OFF(i) = zero
569
570 ENDDO
571
572 nlevset = nlevset + 1
573
574 RETURN
subroutine lsint4(y1, z1, y2, z2, y, z, fi)
type(xfem_phantom_), dimension(:), allocatable xfem_phantom
type(xfem_shell_), dimension(:), allocatable crkshell
type(xfem_lvset_), dimension(:), allocatable crklvset