40
41
42
43
44
46 use element_mod , only : nixc
47
48
49
50#include "implicit_f.inc"
51
52
53
54#include "units_c.inc"
55
56
57
58 INTEGER NEL,NFT,ILAY,NLAY
59 INTEGER IXC(NIXC,*),NGL(NEL),IEL_CRK(*),ELCRKINI(NLAY,NEL),
60 . NODEDGE(2,*),XEDGE4N(4,*)
61 my_real dir1(nlay,nel),dir2(nlay,nel),crklen(nel),aldt(nel),
62 . xl2(nel),yl2(nel),xl3(nel),yl3(nel),xl4(nel),yl4(nel)
63 TYPE (XFEM_EDGE_) , DIMENSION(*) :: CRKEDGE
64
65
66
67 INTEGER I,K,R,IR,P1,P2,NEWCRK,IED,OK,ELCRK,
68 . FAC,IEDGE,ICUT,NOD1,NOD2
69
70 INTEGER JCT(NEL),EDGEL(4,NEL),TIP(NEL),ECUT(2,NEL),dd(4),d(8),KPERM(8)
71
73 . xin(2,nel),yin(2,nel),len(4,nel),
74 . xxl(4,nel),yyl(4,nel),beta0(4,nel)
76 . xint,yint,cross,acd,bcd,dlx,dly,
77 . m12,mm,xint0,yint0,dir11,dir22,
78 . beta,bmin,bmax
79
80 DATA d/1,2,2,3,4,3,1,4/
81 DATA dd/2,3,4,1/
82 DATA kperm/1,2,3,4,1,2,3,4/
83 parameter(bmin = 0.01, bmax = 0.99)
84
85 newcrk = 0
86 DO i=1,nel
87 jct(i) = 0
88 IF (elcrkini(ilay,i) == 5) THEN
89 newcrk = newcrk + 1
90 jct(newcrk) = i
91 elcrkini(ilay,i) = 2
92 ELSEIF (elcrkini(ilay,i) == -5) THEN
93 crklen(i) = aldt(i)
94 elcrkini(ilay,i) = 0
95 ENDIF
96 ENDDO
97 IF (newcrk == 0) RETURN
98
99 DO ir=1,newcrk
100 i = jct(ir)
101 tip(i) = 0
102 ecut(1:2,i) = 0
103 edgel(1:4,i) = 0
104 beta0(1:4,i) = zero
105 xin(1,i) = zero
106 yin(1,i) = zero
107 xin(2,i) = zero
108 yin(2,i) = zero
109
110 xxl(1,i) = zero
111 yyl(1,i) = zero
112 xxl(2,i) = xl2(i)
113 yyl(2,i) = yl2(i)
114 xxl(3,i) = xl3(i)
115 yyl(3,i) = yl3(i)
116 xxl(4,i) = xl4(i)
117 yyl(4,i) = yl4(i)
118
119 len(1,i) = xl2(i)*xl2(i) + yl2(i)*yl2(i)
120 len(2,i) = (xl3(i)-xl2(i))*(xl3(i)-xl2(i))+
121 . (yl3(i)-yl2(i))*(yl3(i)-yl2(i))
122 len(3,i) = (xl4(i)-xl3(i))*(xl4(i)-xl3(i))+
123 . (yl4(i)-yl3(i))*(yl4(i)-yl3(i))
124 len(4,i) = xl4(i)*xl4(i) + yl4(i)*yl4(i)
125 ENDDO
126
127
128
129 DO ir=1,newcrk
130 i = jct(ir)
131 elcrk = iel_crk(i+nft)
132 ok = 0
133 icut = 0
134 ied = 0
135 DO k=1,4
136 iedge = xedge4n(k,elcrk)
137 IF (iedge > 0) THEN
138 icut = crkedge(ilay)%ICUTEDGE(iedge)
139 IF (icut == 1) THEN
140 nod1 = nodedge(1,iedge)
141 nod2 = nodedge(2,iedge)
142 IF (nod1 == ixc(k+1,i) .and. nod2 == ixc(dd(k)+1,i)) THEN
143 p1 = k
144 p2 = dd(k)
145 ELSE IF (nod2 == ixc(k+1,i) .and. nod1 == ixc(dd(k)+1,i)) THEN
146 p1 = dd(k)
147 p2 = k
148 ENDIF
149 ok = 1
150 ied = k
151 ecut(1,i)= k
152 EXIT
153 ENDIF
154 ENDIF
155 ENDDO
156
157 IF (ok == 1) THEN
158 beta = crkedge(ilay)%RATIO(iedge)
159 xin(1,i) = xxl(p1,i) + beta*(xxl(p2,i) - xxl(p1,i))
160 yin(1,i) = yyl(p1,i) + beta*(yyl(p2,i) - yyl(p1,i))
161
162 edgel(ied,i) = 1
163 iedge = xedge4n(ied,elcrk)
164 tip(i) = crkedge(ilay)%EDGETIP(1,iedge)
165 ELSE
166 WRITE(iout,*) 'ERROR IN ADVANCING CRACK --- CHECK CRACK TIP'
168 ENDIF
169
170 END DO
171
172
173
174 DO ir=1,newcrk
175 i = jct(ir)
176 elcrk = iel_crk(i+nft)
177 r = ecut(1,i)
178 xint0 = xin(1,i)
179 yint0 = yin(1,i)
180 dir11 =-dir2(ilay,i)
181 dir22 = dir1(ilay,i)
182
183 IF (dir11 == zero) THEN
184 DO k=1,3
185 r = kperm(ecut(1,i) + k)
186 iedge = xedge4n(r,elcrk)
187 nod1 = nodedge(1,iedge)
188 nod2 = nodedge(2,iedge)
189 IF (nod1 == ixc(r+1,i) .and. nod2 == ixc(dd(r)+1,i))THEN
190 p1 = r
191 p2 = dd(r)
192 ELSE IF (nod2 == ixc(r+1,i).and.nod1 == ixc(dd(r)+1,i))THEN
193 p1 = dd(r)
194 p2 = r
195 ENDIF
196 dlx = xxl(p2,i) - xxl(p1,i)
197 IF (dlx /= zero) THEN
198 dly = yyl(p2,i) - yyl(p1,i)
199 m12 = dly / dlx
200 xint = xint0
201 yint = yyl(p1,i) + m12*(xint-xxl(p1,i))
202 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero .and.
203 . (yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
204 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
205 beta = sqrt(cross / len(r,i))
206 IF (beta > bmax .OR. beta < bmin) THEN
207 beta =
max(beta, bmin)
208 beta =
min(beta, bmax)
209 yint = yyl(p1,i) + beta*(yyl(p2,i)-yyl(p1,i))
210 ENDIF
211
212 ecut(2,i) = r
213 xin(2,i) = xint
214 yin(2,i) = yint
215 edgel(r,i) = 2
216 beta0(r,i) = beta
217 EXIT
218 ENDIF
219 ENDIF
220 ENDDO
221
222 ELSEIF (dir22 == zero) THEN
223 DO k=1,3
224 r = kperm(ecut(1,i) + k)
225 iedge = xedge4n(r,elcrk)
226 nod1 = nodedge(1,iedge)
227 nod2 = nodedge(2,iedge)
228 IF (nod1 == ixc(r+1,i) .and. nod2 == ixc(dd(r)+1,i)) THEN
229 p1 = r
230 p2 = dd(r)
231 ELSE IF (nod2 == ixc(r+1,i).and.nod1 == ixc(dd(r)+1,i)) THEN
232 p1 = dd(r)
233 p2 = r
234 ENDIF
235 dly = yyl(p2,i) - yyl(p1,i)
236 IF (dly /= zero) THEN
237 dlx = xxl(p2,i) - xxl(p1,i)
238 m12 = dlx / dly
239 yint = yint0
240 xint = xxl(p1,i) + m12*(yint-yyl(p1,i))
241 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero .and.
242 . (yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
243 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
244 beta = sqrt(cross / len(r,i))
245 IF (beta > bmax .OR. beta < bmin) THEN
246 beta =
max(beta, bmin)
247 beta =
min(beta, bmax)
248 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
249 ENDIF
250
251 ecut(2,i) = r
252 xin(2,i) = xint
253 yin(2,i) = yint
254 edgel(r,i) = 2
255 beta0(r,i) = beta
256 EXIT
257 ENDIF
258 ENDIF
259 ENDDO
260
261 ELSEIF (dir11 /= zero .and. dir22 /= zero) THEN
262 DO k=1,3
263 r = kperm(ecut(1,i) + k)
264 iedge = xedge4n(r,elcrk)
265 nod1 = nodedge(1,iedge)
266 nod2 = nodedge(2,iedge)
267 IF (nod1 == ixc(r+1,i) .and. nod2 == ixc(dd(r)+1,i)) THEN
268 p1 = r
269 p2 = dd(r)
270 ELSE IF (nod2 == ixc(r+1,i).and.nod1 == ixc(dd(r)+1,i)) THEN
271 p1 = dd(r)
272 p2 = r
273 ENDIF
274
275 dlx = xxl(p2,i) - xxl(p1,i)
276 dly = yyl(p2,i) - yyl(p1,i)
277 mm = dir22/dir11
278 IF (dlx == zero) THEN
279 xint = xxl(p1,i)
280 yint = yint0 + mm*(xint-xint0)
281 IF ((yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
282 cross = (yyl(p1,i) - yint)**2
283 beta = sqrt(cross / len(r,i))
284 IF (beta > bmax .OR. beta < bmin) THEN
285 beta =
max(beta, bmin)
286 beta =
min(beta, bmax)
287 yint = yyl(p1,i) + beta*(yyl(p2,i)-yyl(p1,i))
288 ENDIF
289 ecut(2,i) = r
290 xin(2,i) = xint
291 yin(2,i) = yint
292 edgel(r,i) = 2
293 beta0(r,i) = beta
294 EXIT
295 ENDIF
296 ELSEIF (dly == zero) THEN
297 yint = yyl(p1,i)
298 xint = xint0 + (yint0-yyl(p1,i)) / mm
299 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero) THEN
300 cross = (xxl(p1,i) - xint)**2
301 beta = sqrt(cross / len(r,i))
302 IF (beta > bmax .OR. beta < bmin) THEN
303 beta =
max(beta, bmin)
304 beta =
min(beta, bmax)
305 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
306 ENDIF
307 ecut(2,i) = r
308 xin(2,i) = xint
309 yin(2,i) = yint
310 edgel(r,i) = 2
311 beta0(r,i) = beta
312 EXIT
313 ENDIF
314 ELSE
315 m12 = dly / dlx
316 IF (mm /= m12) THEN
317 xint = (yint0-yyl(p1,i) + m12*xxl(p1,i) - mm*xint0)/(m12-mm)
318 yint = yint0 + mm*(xint-xint0)
319 acd = (yint-yyl(p1,i))*(xint0 - xxl(p1,i))
320 . - (xint-xxl(p1,i))*(yint0 - yyl(p1,i))
321 bcd = (yint-yyl(p2,i))*(xint0 - xxl(p2,i))
322 . - (xint-xxl(p2,i))*(yint0 - yyl(p2,i))
323 IF (acd*bcd <= zero) THEN
324 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
325 beta = sqrt(cross / len(r,i))
326 IF (beta > bmax .OR. beta < bmin) THEN
327 beta =
max(beta, bmin)
328 beta =
min(beta, bmax)
329 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
330 yint = yyl(p1,i) + beta*(yyl(p2,i)-yyl(p1,i))
331 ENDIF
332 ecut(2,i) = r
333 xin(2,i) = xint
334 yin(2,i) = yint
335 edgel(r,i) = 2
336 beta0(r,i) = beta
337 EXIT
338 ENDIF
339 ENDIF
340 ENDIF
341 ENDDO
342 ENDIF
343 ENDDO
344
345
346
347 DO ir=1,newcrk
348 i = jct(ir)
349 fac = 0
350 DO r=1,4
351 IF (edgel(r,i)==1 .or. edgel(r,i)==2) fac=fac+1
352 ENDDO
353 IF (fac /= 2) THEN
354 WRITE(iout,*) 'ERROR IN ADVANCING CRACK. NO CUT EDGES'
356 ENDIF
357 crklen(i) = sqrt((xin(2,i) - xin(1,i))
358 ENDDO
359
360 RETURN