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