39
41 USE elbufdef_mod
42
43
44
45#include "implicit_f.inc"
46
47
48
49#include "mvsiz_p.inc"
50
51
52
53#include "com_xfem1.inc"
54#include "units_c.inc"
55
56
57
58 INTEGER NFT,JFT,JLT,NXLAY
59 INTEGER IXC(NIXC,*),INOD_CRK(*),KNOD2ELC(*),IADC_CRK(4,*),
60 . IEL_CRK(*),ELCUTC(2,*),NODEDGE(2,*),CRKNODIAD(*),(4,*)
62 TYPE (ELBUF_STRUCT_) :: ELBUF_STR
63 TYPE (ELBUF_STRUCT_), DIMENSION(NXEL) :: XFEM_STR
64 TYPE (XFEM_EDGE_) , DIMENSION(*) :: CRKEDGE
65
66
67
68 INTEGER I,J,K,ELCRK,ILEV,ELCUT,pp1,pp2,pp3,IADC(4),IENR0(4),
69 . IENR(4),IED,IEDGE,r,IE10,IE20,IE1,IE2,NOD1,NOD2,N(4),NX(4),
70 . DD(4),ISIGN1,ISIGN2,ISIGN3,ISIGN4,IAD1,IAD2,IAD3,IAD4,
71 . ISIGN0(NXEL,4),p1,p2,LAYCUT,ICUTEDGE,IBOUNDEDGE,
72 . NTAG(4),EDGEENR(4),ENR(4),
73 . ,ITRI,NX1,NX2,NX3,NX4,NM,NP
75 . x1g(mvsiz),x2g(mvsiz),x3g(mvsiz),x4g(mvsiz),
76 . y1g(mvsiz),y2g(mvsiz),y3g(mvsiz),y4g(mvsiz),
77 . z1g(mvsiz),z2g(mvsiz),z3g(mvsiz),z4g(mvsiz),
area(mvsiz),
78 . lxyz0(3),rx(mvsiz),ry(mvsiz),rz(mvsiz),
79 . sx(mvsiz),sy(mvsiz),sz(mvsiz),r11(mvsiz),r12(mvsiz),
80 . r13(mvsiz),r21(mvsiz),r22(mvsiz),r23(mvsiz),r31(mvsiz),
81 . r32(mvsiz),r33(mvsiz),xl1(mvsiz),yl1(mvsiz),xl2(mvsiz),
82 . yl2(mvsiz),xl3(mvsiz),yl3(mvsiz),xl4(mvsiz),yl4(mvsiz),
83 . fit(4,mvsiz),offg(mvsiz),xin(2,mvsiz),yin(2,mvsiz),
84 . xxl(4,mvsiz),yyl(4,mvsiz),xn(4),yn(4),dx(8),xm(2),ym(2)
85 my_real xxx,yyy,zzz,fi,beta,x10,y10,z10,x20,y20,z20,
86 . x1,y1,x2,y2,x3,y3,x4,y4,area1,area2,area3
87 DATA dd/2,3,4,1/
88 DATA dx/1,2,3,4,1,2,3,4/
89
90
91
92 x1 = zero
93 x2 = zero
94 x3 = zero
95 x4 = zero
96 y1 = zero
97 y2 = zero
98 y3 = zero
99 y4 = zero
100 area1 = zero
101 area2 = zero
102 area3 = zero
103 nx1 = 0
104 nx2 = 0
105 nx3 = 0
106 nx4 = 0
107 nm = 0
108 np = 0
109 p1 = 0
110 p2 = 0
111 pp1 = 0
112 pp2 = 0
113 DO i=jft,jlt
114 x1g(i)=x(1,ixc(2,i+nft))
115 y1g(i)=x(2,ixc(2,i+nft))
116 z1g(i)=x(3,ixc(2,i+nft))
117 x2g(i)=x(1,ixc(3,i+nft))
118 y2g(i)=x(2,ixc(3,i+nft))
119 z2g(i)=x(3,ixc(3,i+nft))
120 x3g(i)=x(1,ixc(4,i+nft))
121 y3g(i)=x(2,ixc(4,i+nft))
122 z3g(i)=x(3,ixc(4,i+nft))
123 x4g(i)=x(1,ixc(5,i+nft))
124 y4g(i)=x(2,ixc(5,i+nft))
125 z4g(i)=x(3,ixc(5,i+nft))
126 ENDDO
127
128
129
130 DO i=jft,jlt
131 rx(i) = x2g(i)+x3g(i)-x1g(i)-x4g(i)
132 sx(i) = x3g(i)+x4g(i)-x1g(i)-x2g(i)
133 ry(i) = y2g(i)+y3g(i)-y1g(i)-y4g(i)
134 sy(i) = y3g(i)+y4g(i)-y1g(i)-y2g(i)
135 rz(i) = z2g(i)+z3g(i)-z1g(i)-z4g(i)
136 sz(i) = z3g(i)+z4g(i)-z1g(i)-z2g(i)
137 offg(i) = elbuf_str%GBUF%OFF(i)
138 ENDDO
139 k = 0
141 . rx, ry, rz,
142 . sx, sy, sz,
143 . r11,r12,r13,r21,r22,r23,r31,r32,r33,
area,offg )
144
145
146
147 DO i=jft,jlt
148 lxyz0(1)=fourth*(x3g(i)+x4g(i)+x1g(i)+x2g(i))
149 lxyz0(2)=fourth*(y3g(i)+y4g(i)+y1g(i)+y2g(i))
150 lxyz0(3)=fourth*(z3g(i)+z4g(i)+z1g(i)+z2g(i))
151 xxx = x1g(i)-lxyz0(1)
152 yyy = y1g(i)-lxyz0(2)
153 zzz = z1g(i)-lxyz0(3)
154 xl1(i)=r11(i)*xxx+r21(i)*yyy+r31(i)*zzz
155 yl1(i)=r12(i)*xxx+r22(i)*yyy+r32(i)*zzz
156 xxx = x2g(i)-lxyz0(1)
157 yyy = y2g(i)-lxyz0(2)
158 zzz = z2g(i)-lxyz0(3)
159 xl2(i)=r11(i)*xxx+r21(i)*yyy+r31(i)*zzz
160 yl2(i)=r12(i)*xxx+r22(i)*yyy+r32(i)*zzz
161 xxx = x3g(i)-lxyz0(1)
162 yyy = y3g(i)-lxyz0(2)
163 zzz = z3g(i)-lxyz0(3)
164 xl3(i)=r11(i)*xxx+r21(i)*yyy+r31(i)*zzz
165 yl3(i)=r12(i)*xxx+r22(i)*yyy+r32(i)*zzz
166 xxx = x4g(i)-lxyz0(1)
167 yyy = y4g(i)-lxyz0(2)
168 zzz = z4g(i)-lxyz0(3)
169 xl4(i)=r11(i)*xxx+r21(i)*yyy+r31(i)*zzz
170 yl4(i)=r12(i)*xxx+r22(i)*yyy+r32(i)*zzz
172 ENDDO
173
174
175
176 DO ilay=1,nxlay
177
178 pp1 = nxel*(ilay-1)+1
179 pp2 = pp1 + 1
180 pp3 = pp2 + 1
181
182 DO i=jft,jlt
183 fit(1,i)=zero
184 fit(2,i)=zero
185 fit(3,i)=zero
186 fit(4,i)=zero
187 ENDDO
188
189 DO i=jft,jlt
190 elcrk = iel_crk(i+nft)
191 laycut = crkedge(ilay)%LAYCUT(elcrk)
192 IF (laycut /= 0) THEN
193 xn(1) = xl1(i)
194 yn(1) = yl1(i)
195 xn(2) = xl2(i)
196 yn(2) = yl2(i)
197 xn(3) = xl3(i)
198 yn(3) = yl3(i)
199 xn(4) = xl4(i)
200 yn(4) = yl4(i)
201 xxl(1,i)= xl1(i)
202 yyl(1,i)= yl1(i)
203 xxl(2,i)= xl2(i)
204 yyl(2,i)= yl2(i)
205 xxl(3,i)= xl3(i)
206 yyl(3,i)= yl3(i)
207 xxl(4,i)= xl4(i)
208 yyl(4,i)= yl4(i)
209 DO k=1,4
210 p1 = k
211 p2 = dd(k)
212 ied = crkedge(ilay)%IEDGEC(k,elcrk)
213 IF (ied > 0) THEN
214 iedge = xedge4n(k,elcrk)
215 beta = crkedge(ilay)%RATIO(iedge)
216 nod1 = nodedge(1,iedge)
217 nod2 = nodedge(2,iedge)
218 IF (nod1 == ixc(k+1,i+nft) .and. nod2 == ixc(dd(k)+1,i+nft)) THEN
219 p1 = k
220 p2 = dd(k)
221 ELSEIF (nod2 == ixc(k+1,i+nft).and.nod1 == ixc(dd(k)+1,i+nft))THEN
222 p1 = dd(k)
223 p2 = k
224 ENDIF
225 xin(ied,i) = xn(p1) + beta*(xn(p2) - xn(p1))
226 yin(ied,i) = yn(p1) + beta*(yn(p2) - yn(p1))
227 xm(ied) = half*(xn(p1)+xn(p2))
228 ym(ied) = half*(yn(p1)+yn(p2))
229 ENDIF
230 ENDDO
231
232 DO k=1,4
233 fi=zero
234 CALL lsint4(xm(1),ym(1),xm(2),ym(2),xn(k),yn(k),fi )
235 IF (fit(k,i)==zero) fit(k,i) = fi
236 ENDDO
237 ENDIF
238 ENDDO
239
240 IF (ilay == 1) THEN
241 DO i=jft,jlt
242 elcrk = iel_crk(i+nft)
243 elcut = crkedge(ilay)%LAYCUT(elcrk)
244 IF (elcut /= 0) THEN
245 elcutc(1,i+nft) = 1
246 elcutc(2,i+nft) = 1
247 ENDIF
248 ENDDO
249 ENDIF
250
251
252
253 DO i=jft,jlt
254 elcrk = iel_crk(i+nft)
255 laycut = crkedge(ilay)%LAYCUT(elcrk)
256 IF (laycut /= 0) THEN
257
258 iadc(1) = iadc_crk(1,elcrk)
259 iadc(2) = iadc_crk(2,elcrk)
260 iadc(3) = iadc_crk(3,elcrk)
261 iadc(4) = iadc_crk(4,elcrk)
262
263 ienr0(1) = crknodiad(iadc(1))
264 ienr0(2) = crknodiad(iadc(2))
265 ienr0(3) = crknodiad(iadc(3))
266 ienr0(4) = crknodiad(iadc(4))
267
268 n(1) = ixc(2,i+nft)
269 n(2) = ixc(3,i+nft)
270 n(3) = ixc(4,i+nft)
271 n(4) = ixc(5,i+nft)
272
273 nx(1) = inod_crk(n(1))
274 nx(2) = inod_crk(n(2))
275 nx(3) = inod_crk(n(3))
276 nx(4) = inod_crk(n(4))
277
278 ienr(1) = ienr0(1) + knod2elc(nx(1))*(ilay-1)
279 ienr(2) = ienr0(2) + knod2elc(nx(2))*(ilay-1)
280 ienr(3) = ienr0(3) + knod2elc(nx(3))*(ilay-1)
281 ienr(4) = ienr0(4) + knod2elc(nx(4))*(ilay-1)
282
283 ntag(1:4) = 0
284 edgeenr(1:4) = 0
285 enr(1:4) = 0
286
287 DO r=1,4
288 ied = crkedge(ilay)%IEDGEC(r,elcrk)
289 IF (ied > 0) THEN
290 ntag(r) = ntag(r) + 1
291 ntag(dd(r)) = ntag(dd(r)) + 1
292
293 iedge = xedge4n(r,elcrk)
294 nod1 = nodedge(1,iedge)
295 nod2 = nodedge(2,iedge)
296 ie10 = crkedge(ilay)%EDGEENR(1,iedge)
297 ie20 = crkedge(ilay)%EDGEENR(2,iedge)
298 IF (nod1 == n(r) .and. nod2 == n(dd(r))) THEN
299 p1 = r
300 p2 = dd(r)
301 ELSEIF (nod2 == n(r) .and. nod1 == n(dd(r))) THEN
302 p1 = dd(r)
303 p2 = r
304 ENDIF
305 edgeenr(p1) = ie10
306 edgeenr(p2) = ie20
307 ENDIF
308 ENDDO
309
310 DO r=1,4
311 IF(ntag(r) > 0)THEN
312 enr(r) = edgeenr(r)
313 ELSE
314 enr(r) = ienr(r)
315 ENDIF
316 ENDDO
317
318 DO r=1,4
319 IF (ienr(r) > ienrnod) THEN
320 WRITE(iout,*) 'ERROR CRACK INITIATION,ENRICHMENT NODE EXCEEDED'
322 ENDIF
323 ENDDO
324
325 isign1 = int(sign(one,fit(1,i)))
326 isign2 = int(sign(one,fit(2,i)))
327 isign3 = int(sign(one,fit(3,i)))
328 isign4 = int(sign(one,fit(4,i)))
329
330 IF (fit(1,i) == zero) isign1 = 0
331 IF (fit(2,i) == zero) isign2 = 0
332 IF (fit(3,i) == zero) isign3 = 0
333 IF (fit(4,i) == zero) isign4 = 0
334
335 DO j=1,nxel
336 isign0(j,1) = isign1
337 isign0(j,2) = isign2
338 isign0(j,3) = isign3
339 isign0(j,4) = isign4
340 ENDDO
341
342 DO k=1,4
343 ied = crkedge(ilay)%IEDGEC(k,elcrk)
344 IF (ied > 0) THEN
345 iedge = xedge4n(k,elcrk)
346 nod1 = nodedge(1,iedge)
347 nod2 = nodedge(2,iedge)
348 IF (nod1 == n(k) .and. nod2 == n(dd(k))) THEN
349 p1 = k
350 p2 = dd(k)
351 ELSEIF (nod2 == n(k) .and. nod1 == n(dd(k))) THEN
352 p1 = dd(k)
353 p2 = k
354 ENDIF
355 icutedge = crkedge(ilay)%ICUTEDGE(iedge)
356 iboundedge = crkedge(ilay)%IBORDEDGE(iedge)
357 IF (icutedge == 2 .AND. iboundedge == 0) THEN
358 enr(p1) = -enr(p1)
359 enr(p2) = -enr(p2)
360 ENDIF
361 ENDIF
362 ENDDO
363
364 itri = 0
365 nm = 0
366 np = 0
367 DO k=1,4
368 IF (isign0(1,k) > 0) THEN
369 itri = itri + 1
370 np = k
371 ELSEIF (isign0(1,k) < 0) THEN
372 nm = k
373 ENDIF
374 ENDDO
375
376 IF (itri == 1) THEN
377 itri = -1
378 nx1 = np
379 ELSEIF (itri == 3) THEN
380 itri = 1
381 nx1 = nm
382 ELSEIF (itri == 2) THEN
383 itri = 0
384 nx1 = np
385 IF (np > 0 .and. isign0(1,np-1) > 0) THEN
386 nx1 = np-1
387 ELSE
388 nx1 = np
389 ENDIF
390 ENDIF
391
392 nx2 = dx(nx1+1)
393 nx3 = dx(nx1+2)
394 nx4 = dx(nx1+3)
395 iad1 = iadc(nx1)
396 iad2 = iadc(nx2)
397 iad3 = iadc(nx3)
398 iad4 = iadc(nx4)
401
402
403
404 IF (itri < 0) THEN
405
406 ied = crkedge(ilay)%IEDGEC(nx4,elcrk)
407 x1 = xin(ied,i)
408 y1 = yin(ied,i)
409
410 x2 = xxl(nx3,i)
411 y2 = yyl(nx3,i)
412 x3 = xxl(nx4,i)
413 y3 = yyl(nx4,i)
414 area2 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
415
416 ied = crkedge(ilay)%IEDGEC(nx1,elcrk)
417 x2 = xin(ied,i)
418 y2 = yin(ied,i)
419 x3 = xxl(nx1,i)
420 y3 = yyl(nx1,i)
421 area1 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
422 area1 = area1 /
area(i)
423 area2 = area2 /
area(i)
424 area3 = one - area1 - area2
425 ELSEIF (itri > 0) THEN
426
427 ied = crkedge(ilay)%IEDGEC(nx1,elcrk)
428 x1 = xin(ied,i)
429 y1 = yin(ied,i)
430
431 x2 = xxl(nx2,i)
432 y2 = yyl(nx2,i)
433 x3 = xxl(nx3,i)
434 y3 = yyl(nx3,i)
435 area1 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
436
437 ied = crkedge(ilay)%IEDGEC(nx4,elcrk)
438 x2 = xin(ied,i)
439 y2 = yin(ied,i)
440 x3 = xxl(nx1,i)
441 y3 = yyl(nx1,i)
442 area2 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
443
444 area1 = area1 /
area(i)
445 area2 = area2 /
area(i)
446 area3 = one - area1 - area2
450 ELSE
451 x1 = xxl(nx1,i)
452 y1 = yyl(nx1,i)
453 x2 = xxl(nx2,i)
454 y2 = yyl(nx2,i)
455 ied = crkedge(ilay)%IEDGEC(nx2,elcrk)
456 IF (ied > 0) THEN
457 x3 = xin(ied,i)
458 y3 = yin(ied,i)
459 ELSE
460
461 ENDIF
462 ied = crkedge(ilay)%IEDGEC(nx4,elcrk)
463 IF (ied > 0) THEN
464 x4 = xin(ied,i)
465 y4 = yin(ied,i)
466 ELSE
467
468 ENDIF
469 area1 = half*abs(x1*y2 - x2*y1 + x2*y3 - x3*y2 +
470 . x3*y4 - x4*y3 + x4*y1 - x1*y4)
471 area1 = area1 /
area(i)
472 area2 = one - area1
473 area3 = zero
474 ENDIF
475
479
480
481 ilev = pp1
482
483
484 IF (itri == 0) THEN
485 crklvset(ilev)%ENR0(1,iadc(1)) = abs(enr(1))
486 crklvset(ilev)%ENR0(1,iadc(2)) = enr(2)
487 crklvset(ilev)%ENR0(1,iadc(3)) = enr(3)
488 crklvset(ilev)%ENR0(1,iadc(4)) = enr(4)
489 ELSE
490 crklvset(ilev)%ENR0(1,iadc(1)) = enr(1)
491 crklvset(ilev)%ENR0(1,iadc(2)) = enr(2)
492 crklvset(ilev)%ENR0(1,iadc(3)) = enr(3)
493 crklvset(ilev)%ENR0(1,iadc(4)) = enr(4)
494 ENDIF
495
496 IF(isign0(1,1) > 0)
crklvset(ilev)%ENR0(1,iadc(1)) = 0
497 IF(isign0(1,2) > 0)
crklvset(ilev)%ENR0(1,iadc(2)) = 0
498 IF(isign0(1,3) > 0)
crklvset(ilev)%ENR0(1,iadc(3)) = 0
499 IF(isign0(1,4) > 0)
crklvset(ilev)%ENR0(1,iadc(4)) = 0
500
505
510
511
512 ilev = pp2
513
514
515 crklvset(ilev)%ENR0(1,iadc(1)) = enr(1)
516 crklvset(ilev)%ENR0(1,iadc(2)) = enr(2)
517 crklvset(ilev)%ENR0(1,iadc(3)) = enr(3)
518 crklvset(ilev)%ENR0(1,iadc(4)) = enr(4)
519
520 IF(isign0(2,1) < 0)
crklvset(ilev)%ENR0(1,iadc(1)) = 0
521 IF(isign0(2,2) < 0)
crklvset(ilev)%ENR0
522 IF(isign0(2,3) < 0)
crklvset(ilev)%ENR0(1,iadc(3)) = 0
523 IF(isign0(2,4) < 0)
crklvset(ilev)%ENR0(1,iadc(4)) = 0
524
529
530
531 ilev = pp3
532
533 IF (itri < 0) THEN
534 ie1 = xedge4n(nx2,elcrk)
535 ie2 = xedge4n(nx4,elcrk)
536 IF (crkedge(ilay)%ICUTEDGE(ie2) == 2) THEN
541
542 ELSEIF (crkedge(ilay)%ICUTEDGE(ie1) == 2) THEN
543
544 ENDIF
545
546 ELSEIF (itri > 0) THEN
547
548 ie1 = xedge4n(nx1,elcrk)
549 ie2 = xedge4n(nx4,elcrk)
550 IF (crkedge(ilay)%ICUTEDGE(ie1) == 2) THEN
555 crklvset(pp1)%ENR0(1,iadc(nx1)) = -crknodiad(iadc(nx1)) - knod2elc(nx(nx1))*(ilay-1)
556
557 ELSEIF (crkedge(ilay)%ICUTEDGE(ie2) == 2) THEN
558
559 ENDIF
560 ELSEIF (itri == 0) THEN
561 xfem_str(nxel)%GBUF%OFF(i) = zero
562 xfem_str(nxel)%BUFLY(ilay)%LBUF(1,1,1)%OFF(i) = 0
563 ENDIF
568
569 ENDIF
570 ENDDO
571 ENDDO
572
573 RETURN
subroutine clskew3(jft, jlt, irep, rx, ry, rz, sx, sy, sz, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z, det)
subroutine lsint4(y1, z1, y2, z2, y, z, fi)
subroutine area(d1, x, x2, y, y2, eint, stif0)
type(xfem_phantom_), dimension(:), allocatable xfem_phantom
type(xfem_lvset_), dimension(:), allocatable crklvset