OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
crklayer4n_adv.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| crklayer4n_adv ../engine/source/elements/xfem/crklayer4n_adv.F
25!||--- called by ------------------------------------------------------
26!|| cforc3 ../engine/source/elements/shell/coque/cforc3.F
27!|| czforc3 ../engine/source/elements/shell/coquez/czforc3.F
28!||--- calls -----------------------------------------------------
29!|| arret ../engine/source/system/arret.F
30!|| lsint4 ../engine/source/elements/xfem/crklayer4n_adv.F
31!||--- uses -----------------------------------------------------
32!|| crackxfem_mod ../engine/share/modules/crackxfem_mod.F
33!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
34!||====================================================================
35 SUBROUTINE crklayer4n_adv(
36 . XFEM_STR ,NEL ,NFT ,IXC ,ELCUTC ,
37 . ILAY ,NLAY ,IEL_CRK ,INOD_CRK ,
38 . IADC_CRK ,NODENR ,ELCRKINI ,DIR1 ,DIR2 ,
39 . NODEDGE ,CRKNODIAD ,KNOD2ELC ,CRKEDGE ,A_I ,
40 . XL2 ,XL3 ,XL4 ,YL2 ,YL3 ,
41 . YL4 ,XEDGE4N ,NGL )
42C-----------------------------------------------
43C crack advancement, shells 4N
44C-----------------------------------------------
45C M o d u l e s
46C-----------------------------------------------
47 USE elbufdef_mod
48 USE crackxfem_mod
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53C-----------------------------------------------
54C C o m m o n B l o c k s
55C-----------------------------------------------
56#include "units_c.inc"
57#include "com_xfem1.inc"
58C-----------------------------------------------
59C D u m m y A r g u m e n t s
60C-----------------------------------------------
61 INTEGER NEL,NFT,IXC(NIXC,*),ILAY,NLAY,NGL(NEL),IEL_CRK(*),
62 . INOD_CRK(*),NODENR(*),IADC_CRK(4,*),ELCRKINI(NLAY,NEL),
63 . ELCUTC(2,*),NODEDGE(2,*),CRKNODIAD(*),KNOD2ELC(*),XEDGE4N(4,*)
64 my_real DIR1(NLAY,NEL),DIR2(NLAY,NEL),A_I(NEL),
65 . XL2(NEL),YL2(NEL),XL3(NEL),YL3(NEL),XL4(NEL),YL4(NEL)
66 TYPE (ELBUF_STRUCT_), TARGET :: ELBUF_STR
67 TYPE (ELBUF_STRUCT_), DIMENSION(NXEL) :: XFEM_STR
68 TYPE (XFEM_EDGE_) , DIMENSION(*) :: CRKEDGE
69C-----------------------------------------------
70C L o c a l V a r i a b l e s
71C-----------------------------------------------
72 INTEGER I,J,K,II,R,IR,P1,P2,PP1,PP2,PP3,NEWCRK,IED,OK,
73 . icrk,ielcrk,elcrk,nx1,nx2,nx3,nx4,nm,np,fac,ifi1,ifi2,iedge,
74 . icut,sigbeta,nod1,nod2,ie1,ie2,ie10,ie20,itri,iad1,iad2,iad3,iad4
75c
76 INTEGER ECUT(2,NEL),EDGEL(4,NEL),IENR0(4,NEL),dd(4),d(8),KPERM(8),
77 . N(4),NN(4),IADC(4),IENR(4),ENR(4),ISIGN(4),EN(2,4)
78 INTEGER , DIMENSION(NEL) :: JCT,ELFISS,TIP
79c
80 my_real, DIMENSION(2) :: XMI,YMI
81 my_real, DIMENSION(2,NEL) :: XIN,YIN
82 my_real, DIMENSION(4,NEL) :: len,xxl,yyl,fit,beta0
83 my_real
84 . xint,yint,fi,cross,acd,bcd,dlx,dly,d12,m12,mm,xint0,yint0,beta,
85 . bmin,bmax,dir11,dir22,x1,y1,x2,y2,x3,y3,x4,y4,area1,area2,area3
86C----------
87 DATA d/1,2,2,3,4,3,1,4/
88 DATA dd/2,3,4,1/
89 DATA kperm/1,2,3,4,1,2,3,4/
90 parameter(bmin = 0.01, bmax = 0.99)
91C=======================================================================
92 newcrk = 0
93 DO i=1,nel
94 jct(i) = 0
95 IF (elcrkini(ilay,i) == 1) THEN ! avancement de fissure
96 newcrk = newcrk + 1
97 jct(newcrk) = i
98 ENDIF
99 ENDDO
100 IF (newcrk == 0) RETURN
101C---
102 ii = nxel*(ilay-1)
103 pp1 = ii + 1
104 pp2 = ii + 2
105 pp3 = ii + 3
106C---
107 DO ir=1,newcrk
108 i = jct(ir)
109 elfiss(i)= 0
110 tip(i) = 0
111 ecut(1:2,i) = 0
112 edgel(1:4,i) = 0
113 ienr0(1:4,i) = 0
114 beta0(1:4,i) = zero
115 xin(1,i) = zero ! first inters point in local skew
116 yin(1,i) = zero
117 xin(2,i) = zero ! second inters point in local skew
118 yin(2,i) = zero
119c
120 xxl(1,i) = zero
121 yyl(1,i) = zero
122 xxl(2,i) = xl2(i)
123 yyl(2,i) = yl2(i)
124 xxl(3,i) = xl3(i)
125 yyl(3,i) = yl3(i)
126 xxl(4,i) = xl4(i)
127 yyl(4,i) = yl4(i)
128c
129 len(1,i) = xl2(i)*xl2(i) + yl2(i)*yl2(i)
130 len(2,i) = (xl3(i)-xl2(i))*(xl3(i)-xl2(i))+
131 . (yl3(i)-yl2(i))*(yl3(i)-yl2(i))
132 len(3,i) = (xl4(i)-xl3(i))*(xl4(i)-xl3(i))+
133 . (yl4(i)-yl3(i))*(yl4(i)-yl3(i))
134 len(4,i) = xl4(i)*xl4(i) + yl4(i)*yl4(i)
135 ENDDO
136C------------------------------------------------
137c First intersection (already cut edge)
138C------------------------------------------------
139 DO ir=1,newcrk ! loop over elems (layers) with advancing cracks
140 i = jct(ir)
141 elcrk = iel_crk(i+nft) ! N element sys xfem
142 ok = 0
143 icut = 0
144 ied = 0
145 DO k=1,4 ! edges
146 iedge = xedge4n(k,elcrk)
147 IF (iedge > 0) THEN
148 icut = crkedge(ilay)%ICUTEDGE(iedge)
149 IF (icut == 1) THEN
150 nod1 = nodedge(1,iedge) ! noeud std
151 nod2 = nodedge(2,iedge)
152 IF (nod1 == ixc(k+1,i) .and. nod2 == ixc(dd(k)+1,i)) THEN
153 p1 = k
154 p2 = dd(k)
155 ELSE IF (nod2 == ixc(k+1,i) .and. nod1 == ixc(dd(k)+1,i)) THEN
156 p1 = dd(k)
157 p2 = k
158 ENDIF
159 ok = 1
160 ied = k
161 ecut(1,i)= k
162 EXIT
163 ENDIF ! IF (ICUT == 1) THEN
164 ENDIF
165 ENDDO ! DO K=1,4
166C---
167 IF (ok == 1) THEN ! edge found
168 beta = crkedge(ilay)%RATIO(iedge)
169 xin(1,i) = xxl(p1,i) + beta*(xxl(p2,i) - xxl(p1,i))
170 yin(1,i) = yyl(p1,i) + beta*(yyl(p2,i) - yyl(p1,i))
171c
172 elfiss(i) = crkedge(ilay)%EDGEICRK(iedge) ! Id fissure qui avance
173 ienr0(p1,i) = crkedge(ilay)%EDGEENR(1,iedge) ! enrichissement 1 noeud
174 ienr0(p2,i) = crkedge(ilay)%EDGEENR(2,iedge) ! enrichissement 2 noeud
175 edgel(ied,i) = 1 ! local : premier edge coupe
176 iedge = xedge4n(ied,elcrk) ! N edge element sys xfem
177 tip(i) = crkedge(ilay)%EDGETIP(1,iedge) ! 1 ou 2 , debut ou fin de fissure
178 ELSE
179 WRITE(iout,*) 'ERROR IN ADVANCING CRACK --- CHECK CRACK TIP'
180 CALL arret(2)
181 ENDIF
182C---
183 END DO ! DO IR=1,NEWCRK
184C--------------------------------------------------
185c Search for second intersection (new cut edge)
186C--------------------------------------------------
187 DO ir=1,newcrk
188 i = jct(ir)
189 elcrk = iel_crk(i+nft)
190 r = ecut(1,i)
191 xint0 = xin(1,i)
192 yint0 = yin(1,i)
193 dir11 =-dir2(ilay,i)
194 dir22 = dir1(ilay,i)
195C---
196 IF (dir11 == zero) THEN
197 DO k=1,3
198 r = kperm(ecut(1,i) + k)
199 iedge = xedge4n(r,elcrk)
200 nod1 = nodedge(1,iedge)
201 nod2 = nodedge(2,iedge)
202 IF (nod1 == ixc(r+1,i) .and. nod2 == ixc(dd(r)+1,i))THEN
203 p1 = r
204 p2 = dd(r)
205 ELSE IF (nod2 == ixc(r+1,i).and.nod1 == ixc(dd(r)+1,i))THEN
206 p1 = dd(r)
207 p2 = r
208 ENDIF
209 dlx = xxl(p2,i) - xxl(p1,i)
210 IF (dlx /= zero) THEN
211 dly = yyl(p2,i) - yyl(p1,i)
212 m12 = dly / dlx
213 xint = xint0
214 yint = yyl(p1,i) + m12*(xint-xxl(p1,i))
215 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero .and.
216 . (yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
217 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
218 beta = sqrt(cross / len(r,i))
219 IF (beta > bmax .OR. beta < bmin) THEN
220 beta = max(beta, bmin)
221 beta = min(beta, bmax)
222 yint = yyl(p1,i) + beta*(yyl(p2,i)-yyl(p1,i))
223 ENDIF
224C
225 ecut(2,i) = r
226 xin(2,i) = xint
227 yin(2,i) = yint
228 edgel(r,i) = 2
229 beta0(r,i) = beta
230 EXIT
231 ENDIF
232 ENDIF
233 ENDDO
234c
235 ELSEIF (dir22 == zero) THEN
236 DO k=1,3
237 r = kperm(ecut(1,i) + k)
238 iedge = xedge4n(r,elcrk)
239 nod1 = nodedge(1,iedge)
240 nod2 = nodedge(2,iedge)
241 IF (nod1 == ixc(r+1,i) .and. nod2 == ixc(dd(r)+1,i)) THEN
242 p1 = r
243 p2 = dd(r)
244 ELSE IF (nod2 == ixc(r+1,i).and.nod1 == ixc(dd(r)+1,i)) THEN
245 p1 = dd(r)
246 p2 = r
247 ENDIF
248 dly = yyl(p2,i) - yyl(p1,i)
249 IF (dly /= zero) THEN
250 dlx = xxl(p2,i) - xxl(p1,i)
251 m12 = dlx / dly
252 yint = yint0
253 xint = xxl(p1,i) + m12*(yint-yyl(p1,i))
254 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero .and.
255 . (yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
256 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
257 beta = sqrt(cross / len(r,i))
258 IF (beta > bmax .OR. beta < bmin) THEN
259 beta = max(beta, bmin)
260 beta = min(beta, bmax)
261 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
262 ENDIF
263C
264 ecut(2,i) = r
265 xin(2,i) = xint
266 yin(2,i) = yint
267 edgel(r,i) = 2
268 beta0(r,i) = beta
269 EXIT
270 ENDIF
271 ENDIF
272 ENDDO
273c
274 ELSEIF (dir11 /= zero .and. dir22 /= zero) THEN
275 DO k=1,3
276 r = kperm(ecut(1,i) + k)
277 iedge = xedge4n(r,elcrk)
278 nod1 = nodedge(1,iedge)
279 nod2 = nodedge(2,iedge)
280 IF (nod1 == ixc(r+1,i) .and. nod2 == ixc(dd(r)+1,i)) THEN
281 p1 = r
282 p2 = dd(r)
283 ELSE IF (nod2 == ixc(r+1,i).and.nod1 == ixc(dd(r)+1,i)) THEN
284 p1 = dd(r)
285 p2 = r
286 ENDIF
287C
288 dlx = xxl(p2,i) - xxl(p1,i)
289 dly = yyl(p2,i) - yyl(p1,i)
290 mm = dir22/dir11
291 IF (dlx == zero) THEN
292 xint = xxl(p1,i)
293 yint = yint0 + mm*(xint-xint0)
294 IF ((yint-yyl(p1,i))*(yint-yyl(p2,i)) <= zero) THEN
295 cross = (yyl(p1,i) - yint)**2
296 beta = sqrt(cross / len(r,i))
297 IF (beta > bmax .OR. beta < bmin) THEN
298 beta = max(beta, bmin)
299 beta = min(beta, bmax)
300 yint = yyl(p1,i) + beta*(yyl(p2,i)-yyl(p1,i))
301 ENDIF
302 ecut(2,i) = r
303 xin(2,i) = xint
304 yin(2,i) = yint
305 edgel(r,i) = 2
306 beta0(r,i) = beta
307 EXIT
308 ENDIF
309 ELSEIF (dly == zero) THEN
310 yint = yyl(p1,i)
311 xint = xint0 + (yint0-yyl(p1,i)) / mm
312 IF ((xint-xxl(p1,i))*(xint-xxl(p2,i)) <= zero) THEN
313 cross = (xxl(p1,i) - xint)**2
314 beta = sqrt(cross / len(r,i))
315 IF (beta > bmax .OR. beta < bmin) THEN
316 beta = max(beta, bmin)
317 beta = min(beta, bmax)
318 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
319 ENDIF
320 ecut(2,i) = r
321 xin(2,i) = xint
322 yin(2,i) = yint
323 edgel(r,i) = 2
324 beta0(r,i) = beta
325 EXIT
326 ENDIF
327 ELSE
328 m12 = dly / dlx
329 IF (mm /= m12) THEN
330 xint = (yint0-yyl(p1,i) + m12*xxl(p1,i) - mm*xint0)/(m12-mm)
331 yint = yint0 + mm*(xint-xint0)
332 acd = (yint-yyl(p1,i))*(xint0 - xxl(p1,i))
333 . - (xint-xxl(p1,i))*(yint0 - yyl(p1,i))
334 bcd = (yint-yyl(p2,i))*(xint0 - xxl(p2,i))
335 . - (xint-xxl(p2,i))*(yint0 - yyl(p2,i))
336 IF (acd*bcd <= zero) THEN
337 cross = (xxl(p1,i) - xint)**2 + (yyl(p1,i) - yint)**2
338 beta = sqrt(cross / len(r,i))
339 IF (beta > bmax .OR. beta < bmin) THEN
340 beta = max(beta, bmin)
341 beta = min(beta, bmax)
342 xint = xxl(p1,i) + beta*(xxl(p2,i)-xxl(p1,i))
343 yint = yyl(p1,i) + beta*(yyl(p2,i)-yyl(p1,i))
344 ENDIF
345 ecut(2,i) = r
346 xin(2,i) = xint
347 yin(2,i) = yint
348 edgel(r,i) = 2
349 beta0(r,i) = beta
350 EXIT
351 ENDIF
352 ENDIF
353 ENDIF
354 ENDDO
355 ENDIF
356 ENDDO
357C-----------------------------------------------------------------------
358C check for getting both intersections
359C-----------------------------------------------------------------------
360 DO ir=1,newcrk
361 i=jct(ir)
362 fac = 0
363 DO r=1,4
364 IF (edgel(r,i)==1 .or. edgel(r,i)==2) fac=fac+1
365 ENDDO
366 IF (fac /= 2) THEN
367 WRITE(iout,*) 'ERROR IN ADVANCING CRACK. NO CUT EDGES'
368 CALL arret(2)
369 ENDIF
370 ENDDO
371c tag cut edge numbers of each layer
372 DO ir=1,newcrk
373 i = jct(ir)
374 elcrk = iel_crk(i+nft)
375 DO k=1,4
376 ied = edgel(k,i)
377 crkedge(ilay)%IEDGEC(k,elcrk) = ied ! IED = 0,1,2
378 ENDDO
379 ENDDO
380C----------------------------------------------------------------------
381C SIGN DISTANCE FOR NEW CRACKED LAYER
382C----------------------------------------------------------------------
383 DO ir=1,newcrk
384 i = jct(ir)
385 fit(1,i) = zero
386 fit(2,i) = zero
387 fit(3,i) = zero
388 fit(4,i) = zero
389c
390 DO k=1,4
391 p1 = k
392 p2 = dd(k)
393 ied = edgel(k,i)
394 IF (ied > 0) THEN
395 xmi(ied) = half*(xxl(p1,i)+xxl(p2,i))
396 ymi(ied) = half*(yyl(p1,i)+yyl(p2,i))
397 ENDIF
398 ENDDO
399 DO k=1,4
400 fi = zero
401 CALL lsint4(xmi(1),ymi(1),xmi(2),ymi(2),xxl(k,i),yyl(k,i),fi)
402 IF (fit(k,i) == zero) fit(k,i)=fi
403 ENDDO
404 ENDDO
405C-------------------
406c Loop over newly cut elements
407C-------------------
408 DO ir=1,newcrk
409 i = jct(ir)
410 elcrk = iel_crk(i+nft)
411C
412 k = ecut(2,i) ! second intersection (tip)
413 iedge = xedge4n(k,elcrk)
414 icut = crkedge(ilay)%ICUTEDGE(iedge)
415 IF (icut > 0) THEN ! already cut before => edge connecting two cracks
416 crkedge(ilay)%ICUTEDGE(iedge) = 3 ! 2 cracks sur le meme edge
417 ELSE
418 crkedge(ilay)%ICUTEDGE(iedge) = 2 ! edge cut
419 crkedge(ilay)%RATIO(iedge) = beta0(k,i)
420 ENDIF
421 ENDDO ! IR=1,NEWCRK
422C-----------------------
423C FILL new cut layer - main loop
424C-----------------------
425 DO ir=1,newcrk
426 i = jct(ir)
427 elcrk = iel_crk(i+nft) ! n elem sys xfem
428 elcutc(1,i) = 2 ! flag de l'element std coupe (= 1,2)
429 numelcrk = numelcrk + 1 ! + 1 element coupe (pour debug/anim par element)
430C
431 iadc(1) = iadc_crk(1,elcrk)
432 iadc(2) = iadc_crk(2,elcrk)
433 iadc(3) = iadc_crk(3,elcrk)
434 iadc(4) = iadc_crk(4,elcrk)
435c
436c if(ispmd==6.and.elcrk==1) print*,' ELCRK:,IAD=',ELCRK,IADC(1),IADC(2),IADC(3),IADC(4)
437c if(ispmd==3.and.elcrk==1220) print*,' ELCRK:,IAD=',ELCRK,IADC(1),IADC(2),IADC(3),IADC(4)
438c if(ispmd==1.and.elcrk==595) print*,' ELCRK:,IAD=',ELCRK,IADC(1),IADC(2),IADC(3),IADC(4)
439c if(ispmd==1.and.elcrk==636) print*,' ELCRK:,IAD=',ELCRK,IADC(1),IADC(2),IADC(3),IADC(4)
440c if(nspmd==1.and.(elcrk==4059.or.elcrk==3959.or.elcrk==3960.or.elcrk==3860))
441c
442c
443 n(1) = ixc(2,i) ! n node sys std
444 n(2) = ixc(3,i)
445 n(3) = ixc(4,i)
446 n(4) = ixc(5,i)
447C
448 nn(1) = inod_crk(n(1)) ! n node sys xfem
449 nn(2) = inod_crk(n(2))
450 nn(3) = inod_crk(n(3))
451 nn(4) = inod_crk(n(4))
452C
453 icrk = elfiss(i) ! Id fissure qui avance
454C
455 isign(1) = int(sign(one,fit(1,i)))
456 isign(2) = int(sign(one,fit(2,i)))
457 isign(3) = int(sign(one,fit(3,i)))
458 isign(4) = int(sign(one,fit(4,i)))
459C
460 IF (fit(1,i) == zero) isign(1) = 0
461 IF (fit(2,i) == zero) isign(2) = 0
462 IF (fit(3,i) == zero) isign(3) = 0
463 IF (fit(4,i) == zero) isign(4) = 0
464c--------------------------------------------
465c activate enriched nodes
466c--------------------------------------------
467 DO k=1,4
468 IF (ienr0(k,i) == 0) THEN
469 ienr(k) = crknodiad(iadc(k)) + knod2elc(nn(k))*(ilay-1) ! default => node std
470 ELSE
471 ienr(k) = ienr0(k,i) ! enriched node
472 ENDIF
473 ENDDO
474C
475 sigbeta = 0
476 k = ecut(2,i) ! only new cut edge (tip)
477 iedge = xedge4n(k,elcrk)
478 nod1 = nodedge(1,iedge)
479 nod2 = nodedge(2,iedge)
480 ie10 = crkedge(ilay)%EDGEENR(1,iedge)
481 ie20 = crkedge(ilay)%EDGEENR(2,iedge)
482 IF (nod1 == ixc(k+1,i) .and. nod2 == ixc(dd(k)+1,i)) THEN
483 ie1 = ienr(k)
484 ie2 = ienr(dd(k))
485 ifi1 = isign(k)
486 ifi2 = isign(dd(k))
487 sigbeta = iedge
488 ELSE IF (nod2 == ixc(k+1,i) .and. nod1 == ixc(dd(k)+1,i)) THEN
489 ie1 = ienr(dd(k))
490 ie2 = ienr(k)
491 ifi1 = isign(dd(k))
492 ifi2 = isign(k)
493 sigbeta = -iedge
494 END IF
495 crkedge(ilay)%EDGEENR(1,iedge) = max(ie1,ie10)
496 crkedge(ilay)%EDGEENR(2,iedge) = max(ie2,ie20)
497 IF (crkedge(ilay)%EDGEICRK(iedge) == 0)
498 . crkedge(ilay)%EDGEICRK(iedge) = icrk
499c--------------------------
500 itri = 0
501 nm = 0
502 np = 0
503 DO k=1,4
504 ied = edgel(k,i)
505 IF (ied > 0) THEN
506 iedge = xedge4n(k,elcrk)
507 crkedge(ilay)%EDGETIP(1,iedge) = tip(i)
508 crkedge(ilay)%EDGETIP(2,iedge) =
509 . crkedge(ilay)%EDGETIP(2,iedge) + 1
510 ENDIF
511 IF (isign(k) > 0) THEN
512 itri = itri + 1
513 np = k
514 ELSEIF (isign(k) < 0) THEN
515 nm = k
516 ENDIF
517 ENDDO
518c---
519 IF (itri == 1) THEN
520 itri = -1
521 nx1 = np
522 ELSEIF (itri == 3) THEN
523 itri = 1
524 nx1 = nm
525 ELSEIF (itri == 2) THEN
526 itri = 0
527 IF (np > 1 .and. isign(np-1) > 0) THEN
528 nx1 = np-1
529 ELSE
530 nx1 = np
531 ENDIF
532 ENDIF
533c
534c print*,' Avancement: IEL,ELCRK,proc=',I+NFT,ELCRK,ispmd
535c print*,' IAD=',IADC(1),IADC(2),IADC(3),IADC(4)
536c print*,' NSX=',NN(1),NN(2),NN(3),NN(4),ITRI
537c---
538 nx2 = kperm(nx1+1)
539 nx3 = kperm(nx1+2)
540 nx4 = kperm(nx1+3)
541C
542 crkedge(ilay)%LAYCUT(elcrk) = 1
543 xfem_phantom(ilay)%ITRI(1,elcrk) = itri
544 xfem_phantom(ilay)%ITRI(2,elcrk) = nx1
545 xfem_phantom(ilay)%ELCUT(elcrk) = icrk
546 xfem_phantom(ilay)%IFI(iadc(1)) = isign(1)
547 xfem_phantom(ilay)%IFI(iadc(2)) = isign(2)
548 xfem_phantom(ilay)%IFI(iadc(3)) = isign(3)
549 xfem_phantom(ilay)%IFI(iadc(4)) = isign(4)
550C------------------
551C IXEL = 1 => positif element within ILAY (ILEV = PP1)
552C------------------
553 crklvset(pp1)%ENR0(1,iadc(1)) = -ienr(1)
554 crklvset(pp1)%ENR0(1,iadc(2)) = -ienr(2)
555 crklvset(pp1)%ENR0(1,iadc(3)) = -ienr(3)
556 crklvset(pp1)%ENR0(1,iadc(4)) = -ienr(4)
557C
558 IF (isign(1) > 0) crklvset(pp1)%ENR0(1,iadc(1)) = 0
559 IF (isign(2) > 0) crklvset(pp1)%ENR0(1,iadc(2)) = 0
560 IF (isign(3) > 0) crklvset(pp1)%ENR0(1,iadc(3)) = 0
561 IF (isign(4) > 0) crklvset(pp1)%ENR0(1,iadc(4)) = 0
562C------------------
563C IXEL = 2 => negatif element within ILAY (ILEV = PP2)
564C------------------
565 crklvset(pp2)%ENR0(1,iadc(1)) = -ienr(1)
566 crklvset(pp2)%ENR0(1,iadc(2)) = -ienr(2)
567 crklvset(pp2)%ENR0(1,iadc(3)) = -ienr(3)
568 crklvset(pp2)%ENR0(1,iadc(4)) = -ienr(4)
569C
570 IF (isign(1) < 0) crklvset(pp2)%ENR0(1,iadc(1)) = 0
571 IF (isign(2) < 0) crklvset(pp2)%ENR0(1,iadc(2)) = 0
572 IF (isign(3) < 0) crklvset(pp2)%ENR0(1,iadc(3)) = 0
573 IF (isign(4) < 0) crklvset(pp2)%ENR0(1,iadc(4)) = 0
574C------------------
575C IXEL = 3 => Third phantom element (ILEV = PP3)
576C------------------
577 IF (itri == 0) THEN
578 x1 = xxl(nx1,i)
579 y1 = yyl(nx1,i)
580 x2 = xxl(nx2,i)
581 y2 = yyl(nx2,i)
582 ied = edgel(nx2,i)
583 IF (ied > 0) THEN
584 x3 = xin(ied,i)
585 y3 = yin(ied,i)
586 ELSE
587 ENDIF
588 ied = edgel(nx4,i)
589 IF (ied > 0) THEN
590 x4 = xin(ied,i)
591 y4 = yin(ied,i)
592 ELSE
593 ENDIF
594 area1 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
595 .
596 area1 = area1 * a_i(i)
597 area2 = one - area1
598 area3 = zero
599c
600 xfem_str(3)%GBUF%OFF(i) = zero
601c XFEM_STR(3)%BUFLY(ILAY)%OFF(I) = 0
602 xfem_str(3)%BUFLY(ilay)%LBUF(1,1,1)%OFF(i) = zero
603c
604 ELSEIF (itri < 0) THEN ! cut en 3, IPH3 = IPH2 < 0
605 ie1 = xedge4n(nx2,elcrk)
606 ie2 = xedge4n(nx4,elcrk)
607c
608 IF (crkedge(ilay)%ICUTEDGE(ie2) > 1) THEN
609 sigbeta = -sigbeta
610 iad1 = iadc(nx1)
611 iad2 = iadc(nx2)
612 iad3 = iadc(nx3)
613 iad4 = iadc(nx4)
614 nod1 = iad4
615 nod2 = iad1
616 crklvset(pp3)%ENR0(1,iad1) = abs(crklvset(pp2)%ENR0(1,iad1))
617 crklvset(pp3)%ENR0(1,iad2) = crklvset(pp2)%ENR0(1,iad2)
618 crklvset(pp3)%ENR0(1,iad3) = crklvset(pp2)%ENR0(1,iad3)
619 crklvset(pp3)%ENR0(1,iad4) = crklvset(pp2)%ENR0(1,iad4)
620 crklvset(pp2)%ENR0(1,iad1) = -crknodiad(iad1) - knod2elc(nn(nx1))*(ilay-1)
621c
622c-- IXEL=2 AREA factor
623 ied = crkedge(ilay)%IEDGEC(nx4,elcrk)
624 x1 = xin(ied,i)
625 y1 = yin(ied,i)
626 x2 = xxl(nx3,i)
627 y2 = yyl(nx3,i)
628 x3 = xxl(nx4,i)
629 y3 = yyl(nx4,i)
630 area2 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
631c-- IXEL=1 AREA factor
632 ied = crkedge(ilay)%IEDGEC(nx1,elcrk)
633 x2 = xin(ied,i)
634 y2 = yin(ied,i)
635 x3 = xxl(nx1,i)
636 y3 = yyl(nx1,i)
637 area1 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
638c
639 area1 = area1 * a_i(i)
640 area2 = area2 * a_i(i)
641 area3 = one - area1 - area2
642 ELSEIF (crkedge(ilay)%ICUTEDGE(ie1) > 1) THEN
643 print*,' IMPOSSIBLE CASE'
644 ENDIF
645C
646 ELSEIF (itri > 0) THEN ! cut en 3, IPH3 = IPH1
647 ie1 = xedge4n(nx1,elcrk)
648 ie2 = xedge4n(nx4,elcrk)
649c
650 IF (crkedge(ilay)%ICUTEDGE(ie1) > 1) THEN
651 iad1 = iadc(nx1)
652 iad2 = iadc(nx2)
653 iad3 = iadc(nx3)
654 iad4 = iadc(nx4)
655 nod1 = iad2
656 nod2 = iad1
657 crklvset(pp3)%ENR0(1,iad1) = abs(crklvset(pp1)%ENR0(1,iad1))
658 crklvset(pp3)%ENR0(1,iad2) = crklvset(pp1)%ENR0(1,iad2)
659 crklvset(pp3)%ENR0(1,iad3) = crklvset(pp1)%ENR0(1,iad3)
660 crklvset(pp3)%ENR0(1,iad4) = crklvset(pp1)%ENR0(1,iad4)
661 crklvset(pp1)%ENR0(1,iad1) = -crknodiad(iad1) - knod2elc(nn(nx1))*(ilay-1)
662c-- IXEL=1 AREA factor
663 ied = crkedge(ilay)%IEDGEC(nx1,elcrk)
664 x1 = xin(ied,i)
665 y1 = yin(ied,i)
666 x2 = xxl(nx2,i)
667 y2 = yyl(nx2,i)
668 x3 = xxl(nx3,i)
669 y3 = yyl(nx3,i)
670 area1 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
671c-- IXEL=2 AREA factor
672 ied = crkedge(ilay)%IEDGEC(nx4,elcrk)
673 x2 = xin(ied,i)
674 y2 = yin(ied,i)
675 x3 = xxl(nx1,i)
676 y3 = yyl(nx1,i)
677 area2 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
678c
679 area1 = area1 * a_i(i)
680 area2 = area2 * a_i(i)
681 area3 = one - area1 - area2
682 ELSEIF (crkedge(ilay)%ICUTEDGE(ie2) > 1) THEN
683 print*,' IMPOSSIBLE CASE'
684 ENDIF
685 ENDIF ! ITRI
686c----
687 crklvset(pp1)%AREA(elcrk) = area1
688 crklvset(pp2)%AREA(elcrk) = area2
689 crklvset(pp3)%AREA(elcrk) = area3
690c
691 IF (area3 < zero .or. area1 > one .or. area2 > one .or. area3 > one ) THEN
692 print*,'ERROR : XFEM PHANTOM ELEMENT AREA: ELCRK=',elcrk
693 ENDIF
694c
695 ENDDO ! IR=1,NEWCRK
696c-----------
697 RETURN
698 END
699!||====================================================================
700!|| lsint4 ../engine/source/elements/xfem/crklayer4n_adv.F
701!||--- called by ------------------------------------------------------
702!|| crklayer3n_adv ../engine/source/elements/xfem/crklayer3n_adv.F
703!|| crklayer3n_ini ../engine/source/elements/xfem/crklayer3n_ini.F
704!|| crklayer4n_adv ../engine/source/elements/xfem/crklayer4n_adv.F
705!|| crklayer4n_ini ../engine/source/elements/xfem/crklayer4n_ini.f
706!|| enrichc_ini ../engine/source/elements/xfem/enrichc_ini.F
707!|| enrichtg_ini ../engine/source/elements/xfem/enrichtg_ini.F
708!||====================================================================
709 SUBROUTINE lsint4(Y1, Z1, Y2, Z2, Y, Z, FI)
710C-----------------------------------------------
711C I m p l i c i t T y p e s
712C-----------------------------------------------
713#include "implicit_f.inc"
714C-----------------------------------------------
715C D u m m y A r g u m e n t s
716C-----------------------------------------------
717 my_real
718 . y1,z1,y2,z2,y,z,fi
719C-----------------------------------------------
720C L o c a l V a r i a b l e s
721C-----------------------------------------------
722 my_real
723 . area,ab
724C-----------------------------------------------
725 fi = zero
726 area = ((y2*z-y*z2)-(y1*z-y*z1)+(y1*z2-z1*y2))
727 ab = (y2-y1)**2 + (z2-z1)**2
728 IF (ab > zero) fi = area/sqrt(ab)
729C-----------
730 RETURN
731 END
my_real
#define my_real
Definition cppsort.cpp:32
lsint4
subroutine lsint4(y1, z1, y2, z2, y, z, fi)
Definition crklayer4n_adv.F:710
crklayer4n_adv
subroutine crklayer4n_adv(xfem_str, nel, nft, ixc, elcutc, ilay, nlay, iel_crk, inod_crk, iadc_crk, nodenr, elcrkini, dir1, dir2, nodedge, crknodiad, knod2elc, crkedge, a_i, xl2, xl3, xl4, yl2, yl3, yl4, xedge4n, ngl)
Definition crklayer4n_adv.F:42
crklayer4n_ini
subroutine crklayer4n_ini(xfem_str, nel, nft, ixc, elcutc, ilay, nlay, iel_crk, inod_crk, iadc_crk, nodenr, elcrkini, dir1, dir2, nodedge, crknodiad, knod2elc, crkedge, a_i, xl2, xl3, xl4, yl2, yl3, yl4, xedge4n, ngl)
Definition crklayer4n_ini.F:42
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
crackxfem_mod
Definition crackxfem_mod.F:121
crackxfem_mod::xfem_phantom
type(xfem_phantom_), dimension(:), allocatable xfem_phantom
Definition crackxfem_mod.F:139
crackxfem_mod::crklvset
type(xfem_lvset_), dimension(:), allocatable crklvset
Definition crackxfem_mod.F:146
arret
subroutine arret(nn)
Definition arret.F:87