36
37
38
40
41
42
43#include "implicit_f.inc"
44
45
46
47#include "com04_c.inc"
48#include "units_c.inc"
49#include "comlock.inc"
50
51
52
53 INTEGER ,NROT
54 INTEGER ,DIMENSION(NIXTG,NEL),INTENT(IN) :: IXTG
55 INTEGER ,DIMENSION(NUMNOD) ,INTENT(IN) :: ITAB
56 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL,FWAVE_EL
57
58 my_real ,
DIMENSION(NEL,NROT) ,
INTENT(IN) :: dir_a
59 my_real ,
DIMENSION(NEL,2) ,
INTENT(IN) :: crkdir
60 my_real ,
DIMENSION(NEL) ,
INTENT(IN) :: xl2,xl3,yl2,yl3
61 TYPE (FAILWAVE_STR_) :: FAILWAVE
62
63
64
65 INTEGER I,II,K,N1,N2,N3,INTERSECTION,INT1,INT2,INT3,INT4,INT5,INT6,
66 . LEVEL,NEWCRK1,NEWCRK2,
67 . NCURR,MAXLEV,FOUND_EXISTING,IDEBUG
68 INTEGER ,DIMENSION(NEL) :: INDX1,INDX2
69 INTEGER ,DIMENSION(3) :: IDF1,IDF2,NOD_ID,NOD_NN
70
71 my_real :: dx1,dx2,dy1,dy2,dir11,dir22,cosx,sinx,cosy,siny,
72 . lmax,xm,ym,x1,y1,x2,y2,x3,y3,x4,y4,x5
73 . rx,ry,rat1,rat2
74 INTEGER SEG_INTERSECT
76
77
78
79 idebug = 0
80
81
82 SELECT CASE (failwave%WAVE_MOD)
83
84 CASE (1)
85
86 DO i=1,nel
87 IF (fwave_el(i) < 0) THEN
88 n1 = failwave%IDXI(ixtg(2,i))
89 n2 = failwave%IDXI(ixtg(3,i))
90 n3 = failwave%IDXI(ixtg(4,i))
91 failwave%FWAVE_NOD_STACK(1,n1,1) = 1
92 failwave%FWAVE_NOD_STACK(1,n2,1) = 1
93 failwave%FWAVE_NOD_STACK(1,n3,1) = 1
94 failwave%MAXLEV_STACK(n1) = 1
95 failwave%MAXLEV_STACK(n2) = 1
96 failwave%MAXLEV_STACK(n3) = 1
97 ENDIF
98 ENDDO
99
100 CASE (2,3)
101
102 newcrk1 = 0
103 newcrk2 = 0
104 DO i=1,nel
105 IF (fwave_el(i) == -1) THEN
106 newcrk1 = newcrk1 + 1
107 indx1(newcrk1) = i
108 ELSEIF (fwave_el(i) == -2) THEN
109 newcrk2 = newcrk2 + 1
110 indx2(newcrk2) = i
111 ELSEIF (fwave_el(i) == -3) THEN
112 newcrk1 = newcrk1 + 1
113 newcrk2 = newcrk2 + 1
114 indx1(newcrk1) = i
115 indx2(newcrk2) = i
116 ENDIF
117 ENDDO
118
119 IF (newcrk1 + newcrk2 > 0) THEN
120
121 rat1 = half * tan(pi/six)
122 rat2 = one - rat1
123
124
125
126 DO ii=1,newcrk1
127 i = indx1(ii)
128 n1 = ixtg(2,i)
129 n2 = ixtg(3,i)
130 n3 = ixtg(4,i)
131 nod_nn(1) = failwave%IDXI(n1)
132 nod_nn(2) = failwave%IDXI(n2)
133 nod_nn(3) = failwave%IDXI(n3)
134 nod_id(1) = itab(n1)
135 nod_id(2) = itab(n2)
136 nod_id(3) = itab(n3)
137 idf1(:) = 0
138 idf2(:) = 0
139
140 idebug = 0
141
142
143 IF (nrot == 0) THEN
144 dir11 = -crkdir(i,2)
145 dir22 = crkdir(i,1)
146 ELSE
147 cosx = dir_a(i,1)
148 sinx = dir_a(i,2)
149 cosy =-crkdir(i,2)
150 siny = crkdir(i,1)
151 dir11 = cosx*cosy - sinx*siny
152 dir22 = cosx*siny + sinx*cosy
153 ENDIF
154
155 xm = (xl2(i) + xl3(i)) * third
156 ym = (yl2(i) + yl3(i)) * third
157 lmax =
max(xl2(i)*2 + yl2(i)**2, xl3(i)**2 + yl3(i))
158 lmax = sqrt(lmax) * two
159
160 dx1 = xm - dir11 * lmax
161 dy1 = ym - dir22 * lmax
162 dx2 = xm + dir11 * lmax
163 dy2 = ym + dir22 * lmax
164
165 x1 = xl2(i)*rat1
166 y1 = yl2(i)*rat1
167 x2 = xl2(i)*rat2
168 y2 = yl2(i)*rat2
169 rx = xl3(i) - xl2(i)
170 ry = yl3(i) - yl2(i)
171 x3 = xl2(i) + rx * rat1
172 y3 = yl2(i) + ry * rat1
173 x4 = xl2(i) + rx * rat2
174 y4 = yl2(i) + ry * rat2
175 rx =-xl3(i)
176 ry =-yl3(i)
177 x5 = xl3(i) + rx * rat1
178 y5 = yl3(i) + ry * rat1
179 x6 = xl3(i) + rx * rat2
180 y6 = yl3(i) + ry * rat2
181
182 int1 =
seg_intersect(x6,y6,x1,y1,dx1,dy1,dx2,dy2,xint,yint,idebug)
183 int2 =
seg_intersect(x2,y2,x3,y3,dx1,dy1,dx2,dy2,xint,yint,idebug)
184 int3 =
seg_intersect(x4,y4,x5,y5,dx1,dy1,dx2,dy2,xint,yint,idebug)
185 int4 =
seg_intersect(x1,y1,x2,y2,dx1,dy1,dx2,dy2,xint,yint,idebug)
186 int5 =
seg_intersect(x3,y3,x4,y4,dx1,dy1,dx2,dy2,xint,yint,idebug)
187 int6 =
seg_intersect(x5,y5,x6,y6,dx1,dy1,dx2,dy2,xint,yint,idebug)
188
189 intersection = 0
190 IF (int1 == 1) THEN
191 idf1(1) = nod_id(3)
192 idf2(1) = nod_id(2)
193 idf1(2) = nod_id(3)
194 idf1(3) = nod_id(2)
195 intersection = 1
196 ELSE IF (int2 == 1) THEN
197 idf1(1) = nod_id(3)
198 idf1(3) = nod_id(1)
199 idf1(2) = nod_id(1)
200 idf2(2) = nod_id(3)
201 intersection = 1
202 ELSE IF (int3 == 1) THEN
203 idf1(1) = nod_id(2)
204 idf1(2) = nod_id(1)
205 idf1(3) = nod_id(2)
206 idf2(3) = nod_id(1)
207 intersection = 1
208 ELSE IF (int4 == 1) THEN
209 idf1(1) = nod_id(2)
210 idf1(2) = nod_id(1)
211 idf1(3) = nod_id(2)
212 idf2(3) = nod_id(1)
213 intersection = 1
214 ELSE IF (int5 == 1) THEN
215 idf1(1) = nod_id(3)
216 idf2(1) = nod_id(2)
217 idf1(2) = nod_id(3)
218 idf1(3) = nod_id(2)
219 intersection = 1
220 ELSE IF (int6 == 1) THEN
221 idf1(1) = nod_id(3)
222 idf1(3) = nod_id(1)
223 idf1(2) = nod_id(1)
224 idf2(2) = nod_id(3)
225 intersection = 1
226 END IF
227
228 IF (intersection == 1) THEN
229 DO k=1,3
230 ncurr = nod_nn(k)
231
232
233 failwave%MAXLEV_STACK(ncurr) = failwave%MAXLEV_STACK(ncurr) + 1
234 maxlev = failwave%MAXLEV_STACK(ncurr)
235
236
237 IF (maxlev > failwave%SIZE) THEN
238#include "lockon.inc"
239 WRITE(iout,*) 'ERROR IN FAILWAVE PROPAGATION: ELEMENT =',ngl(i),
240 . 'LEVEL=',maxlev
241#include "lockoff.inc"
242 maxlev = failwave%SIZE
243 failwave%MAXLEV_STACK(ncurr) = maxlev
244 ENDIF
245 failwave%FWAVE_NOD_STACK(1,ncurr,maxlev) = idf1(k)
246 failwave%FWAVE_NOD_STACK(2,ncurr,maxlev) = idf2(k)
247 END DO
248
249 ELSE
250
251#include "lockon.inc"
252 WRITE(iout,*) 'ERROR ADVANCING CRACK IN ELEMENT, 1st DIR ',ngl(i)
253#include "lockoff.inc"
254 ENDIF
255
256 ENDDO
257
258
259
260 DO ii=1,newcrk2
261 i = indx2(ii)
262 n1 = ixtg(2,i)
263 n2 = ixtg(3,i)
264 n3 = ixtg(4,i)
265 nod_nn(1) = failwave%IDXI(n1)
266 nod_nn(2) = failwave%IDXI(n2)
267 nod_nn(3) = failwave%IDXI(n3)
268 nod_id(1) = itab(n1)
269 nod_id(2) = itab(n2)
270 nod_id(3) = itab(n3)
271 idf1(:) = 0
272 idf2(:) = 0
273
274 IF (nrot == 0) THEN
275 dir11 = crkdir(i,1)
276 dir22 = crkdir(i,2)
277 ELSE
278 cosx = dir_a(i,1)
279 sinx = dir_a(i,2)
280 cosy = crkdir(i,1)
281 siny = crkdir(i,2)
282 dir11 = cosx*cosy - sinx*siny
283 dir22 = cosx*siny + sinx*cosy
284 ENDIF
285
286 xm = (xl2(i) + xl3(i)) * third
287 ym = (yl2(i) + yl3(i)) * third
288 lmax =
max(xl2(i)*2 + yl2(i)**2, xl3(i)**2 + yl3(i))
289 lmax = sqrt(lmax) * two
290
291 dx1 = xm - dir11 * lmax
292 dy1 = ym - dir22 * lmax
293 dx2 = xm + dir11 * lmax
294 dy2 = ym + dir22 * lmax
295
296 x1 = xl2(i)*rat1
297 y1 = yl2(i)*rat1
298 x2 = xl2(i)*rat2
299 y2 = yl2(i)*rat2
300 rx = xl3(i) - xl2(i)
301 ry = yl3(i) - yl2(i)
302 x3 = xl2(i) + rx * rat1
303 y3 = yl2(i) + ry * rat1
304 x4 = xl2(i) + rx * rat2
305 y4 = yl2(i) + ry * rat2
306 rx =-xl3(i)
307 ry =-yl3(i)
308 x5 = xl3(i) + rx * rat1
309 y5 = yl3(i) + ry * rat1
310 x6 = xl3(i) + rx * rat2
311 y6 = yl3(i) + ry * rat2
312
313 int1 =
seg_intersect(x6,y6,x1,y1,dx1,dy1,dx2,dy2,xint,yint,idebug)
314 int2 =
seg_intersect(x2,y2,x3,y3,dx1,dy1,dx2,dy2,xint,yint,idebug)
315 int3 =
seg_intersect(x4,y4,x5,y5,dx1,dy1,dx2,dy2,xint,yint,idebug)
316 int4 =
seg_intersect(x1,y1,x2,y2,dx1,dy1,dx2,dy2,xint,yint,idebug)
317 int5 =
seg_intersect(x3,y3,x4,y4,dx1,dy1,dx2,dy2,xint,yint,idebug)
318 int6 =
seg_intersect(x5,y5,x6,y6,dx1,dy1,dx2,dy2,xint,yint,idebug)
319
320 intersection = 0
321 IF (int1 == 1) THEN
322 idf1(1) = nod_id(3)
323 idf2(1) = nod_id(2)
324 idf1(2) = nod_id(3)
325 idf1(3) = nod_id(2)
326 intersection = 1
327 ELSE IF (int2 == 1) THEN
328 idf1(1) = nod_id(3)
329 idf1(3) = nod_id(1)
330 idf1(2) = nod_id(1)
331 idf2(2) = nod_id(3)
332 intersection = 1
333 ELSE IF (int3 == 1) THEN
334 idf1(1) = nod_id(2)
335 idf1(2) = nod_id(1)
336 idf1(3) = nod_id(2)
337 idf2(3) = nod_id(1)
338 intersection = 1
339 ELSE IF (int4 == 1) THEN
340 idf1(1) = nod_id(2)
341 idf1(2) = nod_id(1)
342 idf1(3) = nod_id(2)
343 idf2(3) = nod_id(1)
344 intersection = 1
345 ELSE IF (int5 == 1) THEN
346 idf1(1) = nod_id(3)
347 idf2(1) = nod_id(2)
348 idf1(2) = nod_id(3)
349 idf1(3) = nod_id(2)
350 intersection = 1
351 ELSE IF (int6 == 1) THEN
352 idf1(1) = nod_id(3)
353 idf1(3) = nod_id(1)
354 idf1(2) = nod_id(1)
355 idf2(2) = nod_id(3)
356 intersection = 1
357 END IF
358
359 IF (intersection == 1) THEN
360 DO k=1,3
361 ncurr = nod_nn(k)
362
363
364 failwave%MAXLEV_STACK(ncurr) = failwave%MAXLEV_STACK(ncurr) + 1
365 maxlev = failwave%MAXLEV_STACK(ncurr)
366
367
368 IF (maxlev > failwave%SIZE) THEN
369#include "lockon.inc"
370 WRITE(iout,*) 'ERROR IN FAILWAVE PROPAGATION: ELEMENT =',ngl(i),
371 . 'LEVEL=',maxlev
372#include "lockoff.inc"
373 maxlev = failwave%SIZE
374 failwave%MAXLEV_STACK(ncurr) = maxlev
375 ENDIF
376 failwave%FWAVE_NOD_STACK(1,ncurr,maxlev) = idf1(k)
377 failwave%FWAVE_NOD_STACK(2,ncurr,maxlev) = idf2(k)
378 END DO
379
380 ELSE
381
382#include "lockon.inc"
383 WRITE(iout,*) 'ERROR ADVANCING CRACK IN ELEMENT, 2nd DIR ',ngl(i)
384#include "lockoff.inc"
385 ENDIF
386
387 ENDDO
388 END IF
389
390
391 END SELECT
392
393 RETURN
integer function seg_intersect(x1, y1, x2, y2, x3, y3, x4, y4, xint, yint, idebug)