37
38
39
40
41
43
44
45
46#include "implicit_f.inc"
47
48
49
50#include "units_c.inc"
51#include "com_xfem1.inc"
52
53
54
55 INTEGER NEL,NFT,ILAY,NLAY
56 INTEGER IXTG(NIXTG,*),NGL(NEL),IEL_CRKTG(*),ELCRKINI(NLAY,*),
57 . NODEDGE(2,*),XEDGE3N(3,*)
59 TYPE (XFEM_EDGE_) , DIMENSION(*) :: CRKEDGE
60 my_real,
DIMENSION(NEL) :: xl2,yl2
61
62
63
64 INTEGER I,J,K,IR,p1,p2,NEWCRK,IED,IED1,IED2,FAC,OK,ICRK,
65 . NOD1,NOD2,ELCRK,ELCRKTG,IEDGE,ICUT
66 INTEGER JCT(NEL),EDGEL(3,NEL),ELTIP(NEL),TIP(NEL)
67 INTEGER DD(3),D(6),ISIGN(3),N(3),IENR(3),NN(3),INV(2)
68
69 my_real,
DIMENSION(NEL) :: xl1,yl1
70 my_real,
DIMENSION(2,NEL) :: xin,yin
71 my_real,
DIMENSION(3,NEL) :: xxl,yyl,len
72 my_real beta0(3,nel),xn(3),yn(3),zn(3),xmi(2),ymi(2)
73 my_real beta,xint,yint,bmin,bmax,x10,y10,z10,x20,y20,z20,
74 . m12,mm,cross1,cross12,xint0,yint0,dir11,dir22
75
76 DATA d/1,2,2,3,1,3/
77 DATA dd/2,3,1/
78 parameter(bmin = 0.01, bmax = 0.99)
79
80 newcrk = 0
81 DO i=1,nel
82 jct(i) = 0
83 IF (elcrkini(ilay,i) == 5) THEN
84 newcrk = newcrk + 1
85 jct(newcrk) = i
86 elcrkini(ilay,i) = 2
87 ELSEIF (elcrkini(ilay,i) == -5) THEN
88 crklen(i) = aldt(i)
89 elcrkini(ilay,i) = -1
90 ENDIF
91 ENDDO
92 IF (newcrk == 0) RETURN
93
94 DO i=1,nel
95 beta0(1:3,i) = zero
96 tip(i) = 0
97 edgel(1,i)=0
98 edgel(2,i)=0
99 edgel(3,i)=0
100 xin(1,i) = zero
101 yin(1,i) = zero
102 xin(2,i) = zero
103 yin(2,i) = zero
104 ENDDO
105
106 inv(1) = 2
107 inv(2) = 1
108
109
110
111
112 DO ir=1,newcrk
113 i = jct(ir)
114 elcrktg = iel_crktg(i+nft)
115 ok = 0
116 icut = 0
117 ied = 0
118 DO k=1,3
119 iedge = xedge3n(k,elcrktg)
120 icut = crkedge(ilay)%ICUTEDGE(iedge)
121 nod1 = nodedge(1,iedge)
122 nod2 = nodedge(2,iedge)
123 IF (nod1 == ixtg(k+1,i) .and. nod2 == ixtg(dd(k)+1,i)) THEN
124 p1 = k
125 p2 = dd(k)
126 ELSE IF (nod2 == ixtg(k+1,i) .and. nod1 == ixtg(dd(k)+1,i)) THEN
127 p1 = dd(k)
128 p2 = k
129 ENDIF
130 IF (icut == 1) THEN
131 ok = ok + 1
132 ied = k
133
134 icrk = crkedge(ilay)%EDGEICRK(iedge)
135 EXIT
136 ENDIF
137 ENDDO
138
139 IF (ok /= 1) THEN
140 WRITE(iout,*) 'ERROR IN ADVANCING CRACK --- CHECK CRACK TIP'
142 ENDIF
143
144 edgel(ied,i) = 1
145 iedge = xedge3n(ied,elcrktg)
146 tip(i) = crkedge(ilay)%EDGETIP(1,iedge)
147
148 END DO
149
150
151
152 DO i=1,nel
153 xl1(i) = zero
154 yl1(i) = zero
155 xxl(1,i)=xl1(i)
156 yyl(1,i)=yl1(i)
157 xxl(2,i)=xl2(i)
158 yyl(2,i)=yl2(i)
159 xxl(3,i)=xl3(i)
160 yyl(3,i)=yl3(i)
161 ENDDO
162
163 DO i=1,nel
164 len(1,i) = (xl2(i)-xl1(i))*(xl2(i)-xl1(i))
165 . + (yl2(i)-yl1(i))*(yl2(i)-yl1(i))
166 len(2,i) = (xl3(i)-xl2(i))*(xl3(i)-xl2(i))
167 . + (yl3(i)-yl2(i))*(yl3(i)-yl2(i))
168 len(3,i) = (xl1(i)-xl3(i))*(xl1(i)-xl3(i))
169 . + (yl1(i)-yl3(i))*(yl1(i)-yl3(i))
170 ENDDO
171
172
173
174 DO ir=1,newcrk
175 i=jct(ir)
176 elcrktg = iel_crktg(i+nft)
177 elcrk = elcrktg + ecrkxfec
178 ied1 = 0
179 ied2 = 0
180 DO k=1,3
181 IF(edgel(k,i) > 0)THEN
182 ied1 = edgel(k,i)
183 ied2 = inv(ied1)
184 EXIT
185 END IF
186 END DO
187 DO k=1,3
188 iedge = xedge3n(k,elcrktg)
189 IF (iedge > 0 .and. edgel(k,i) == 1) THEN
190 icut = crkedge(ilay)%ICUTEDGE(iedge)
191 IF (icut == 1) THEN
192 beta = crkedge(ilay)%RATIO(iedge)
193
194 IF (beta > one .or. beta == zero) THEN
195 WRITE(*,*) 'ERROR NEGATIV BETA, NO INTERSECTION!'
197 ENDIF
198
199 nod1 = nodedge(1,iedge)
200 nod2 = nodedge(2,iedge)
201 IF (nod1 == ixtg(k+1,i) .and. nod2 == ixtg(dd(k)+1,i)) THEN
202 p1 = k
203 p2 = dd(k)
204 ELSEIF (nod2 == ixtg(k+1,i).and.nod1==ixtg(dd(k)+1,i)) THEN
205 p1 = dd(k)
206 p2 = k
207 ENDIF
208 x10 = xxl(p1,i)
209 y10 = yyl(p1,i)
210 x20 = xxl(p2,i)
211 y20 = yyl(p2,i)
212
213 xint = x10+beta*(x20-x10)
214 yint = y10+beta*(y20-y10)
215 xin(ied1,i) = xint
216 yin(ied1,i) = yint
217 ENDIF
218 ENDIF
219 ENDDO
220
221 IF (ied1 == 0 .or. ied2 == 0) GOTO 130
222 xint0 = xin(ied1,i)
223 yint0 = yin(ied1,i)
224
225 dir11 = -dir2(ilay,i)
226 dir22 = dir1(ilay,i)
227
228 IF (dir11 == zero) THEN
229 DO 140 k=1,3
230 xint = zero
231 yint = zero
232 elcrktg = iel_crktg(i+nft)
233 elcrk = elcrktg + ecrkxfec
234 iedge = xedge3n(k,elcrktg)
235 nod1 = nodedge(1,iedge)
236 nod2 = nodedge(2,iedge)
237 IF(nod1 == ixtg(k+1,i) .and. nod2 == ixtg(dd(k)+1,i))THEN
238 p1 = k
239 p2 = dd(k)
240 ELSE IF(nod2 == ixtg(k+1,i).and.nod1==ixtg(dd(k)+1,i))THEN
241 p1 = dd(k)
242 p2 = k
243 ENDIF
244
245 IF (edgel(k,i) == ied1) GOTO 140
246 IF (xxl(p1,i) == xxl(p2,i)) GOTO 140
247 m12 = xxl(p2,i)-xxl(p1,i)
248 m12 = (yyl(p2,i)-yyl(p1,i))/m12
249 xint = xint0
250 yint = yyl(p1,i)+m12*(xint-xxl(p1,i))
251 cross12 = (xint-xxl(p1,i))*(xint-xxl(p2,i))+
252 . (yint-yyl(p1,i))*(yint-yyl(p2,i))
253 IF (cross12 > zero) GOTO 140
254
255 cross1 = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
256 beta = sqrt(cross1 / len(k,i))
257 beta =
max(beta, bmin)
258 beta =
min(beta, bmax)
259 beta0(k,i) = beta
260
261 xin(ied2,i) = xint
262 yin(ied2,i) = yint
263 edgel(k,i) = ied2
264 EXIT
265 140 CONTINUE
266 ELSEIF(dir22 == zero)THEN
267 DO 150 k=1,3
268 xint = zero
269 yint = zero
270 elcrktg = iel_crktg(i+nft)
271 elcrk = elcrktg + ecrkxfec
272 iedge = xedge3n(k,elcrktg)
273 nod1 = nodedge(1,iedge)
274 nod2 = nodedge(2,iedge)
275 IF(nod1 == ixtg(k+1,i) .and. nod2 == ixtg(dd(k)+1,i))THEN
276 p1 = k
277 p2 = dd(k)
278 ELSE IF(nod2 == ixtg(k+1,i).and.nod1==ixtg(dd(k)+1,i))THEN
279 p1 = dd(k)
280 p2 = k
281 ENDIF
282
283 IF (edgel(k,i) == ied1) GOTO 150
284 IF (yyl(p1,i) == yyl(p2,i)) GOTO 150
285 m12 = yyl(p2,i)-yyl(p1,i)
286 m12 = (xxl(p2,i)-xxl(p1,i))/m12
287 yint = yint0
288 xint = xxl(p1,i)+m12*(yint-yyl(p1,i))
289 cross12 = (xint-xxl(p1,i))*(xint-xxl(p2,i))+
290 . (yint-yyl(p1,i))*(yint-yyl(p2,i))
291 IF (cross12 > zero) GOTO 150
292
293 cross1 = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
294 beta = sqrt(cross1 / len(k,i))
295 beta =
max(beta, bmin)
296 beta =
min(beta, bmax)
297 beta0(k,i) = beta
298
299 xin(ied2,i) = xint
300 yin(ied2,i) = yint
301 edgel(k,i) = ied2
302 EXIT
303 150 CONTINUE
304 ELSEIF(dir11 /= zero .AND. dir22 /= zero)THEN
305 DO 160 k=1,3
306 xint = zero
307 yint = zero
308 elcrktg = iel_crktg(i+nft)
309 elcrk = elcrktg + ecrkxfec
310 iedge = xedge3n(k,elcrktg)
311 nod1 = nodedge(1,iedge)
312 nod2 = nodedge(2,iedge)
313 IF (nod1 == ixtg(k+1,i) .and. nod2 == ixtg(dd(k)+1,i))THEN
314 p1 = k
315 p2 = dd(k)
316 ELSE IF (nod2 == ixtg(k+1,i).and.nod1==ixtg(dd(k)+1,i))THEN
317 p1 = dd(k)
318 p2 = k
319 ENDIF
320
321 IF (edgel(k,i) == ied1) GOTO 160
322 IF (xxl(p1,i) == xxl(p2,i)) THEN
323 mm = dir22/dir11
324 xint = xxl(p1,i)
325 yint = yint0+mm*(xint-xint0)
326 cross12 = (xint-xxl(p1,i))*(xint-xxl(p2,i))+
327 . (yint-yyl(p1,i))*(yint-yyl(p2,i))
328 IF (cross12 > zero) GOTO 160
329
330 cross1 = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
331 beta = sqrt(cross1 / len(k,i))
332 beta =
max(beta, bmin)
333 beta =
min(beta, bmax)
334 beta0(k,i) = beta
335
336 xin(ied2,i) = xint
337 yin(ied2,i) = yint
338 edgel(k,i) = ied2
339 EXIT
340 ELSE
341 mm = dir22/dir11
342 m12 = xxl(p2,i)-xxl(p1,i)
343 m12 = (yyl(p2,i)-yyl(p1,i))/m12
344 IF (mm == m12) GOTO 160
345 xint = (yint0-yyl(p1,i)+m12*xxl(p1,i)-mm*xint0)/(m12-mm)
346 yint = yint0+mm*(xint-xint0)
347 cross12 = (xint-xxl(p1,i))*(xint-xxl(p2,i
348 . (yint-yyl(p1,i))*(yint-yyl(p2,i))
349 IF (cross12 > zero) GOTO 160
350
351 cross1 = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
352 beta = sqrt(cross1 / len(k,i))
353 beta =
max(beta, bmin)
354 beta =
min(beta, bmax)
355 beta0(k,i) = beta
356
357 xin(ied2,i) = xint
358 yin(ied2,i) = yint
359 edgel(k,i) = ied2
360 EXIT
361 ENDIF
362 160 CONTINUE
363 ENDIF
364 130 CONTINUE
365 ENDDO
366
367
368
369 DO ir=1,newcrk
370 i = jct(ir)
371 fac = 0
372 DO k=1,3
373 IF (edgel(k,i)==1 .or. edgel(k,i)==2) fac=fac+1
374 ENDDO
375 IF (fac /= 2) THEN
376 WRITE(iout,*) 'ERROR IN ADVANCING CRACK.NO CUT EDGES'
378 ENDIF
379 crklen(i) = sqrt((xin(2,i) - xin(1,i))**2 + (yin(2,i) - yin(1,i))**2)
380 ENDDO
381
382 RETURN