OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
enrichtg_ini.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "com_xfem1.inc"
#include "units_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine enrichtg_ini (elbuf_str, ixtg, nft, jft, jlt, nxlay, iad_crktg, iel_crktg, inod_crk, elcutc, nodedge, crknodiad, knod2elc, x, crkedge, xedge3n)

Function/Subroutine Documentation

◆ enrichtg_ini()

subroutine enrichtg_ini ( type (elbuf_struct_) elbuf_str,
integer, dimension(nixtg,*) ixtg,
integer nft,
integer jft,
integer jlt,
integer nxlay,
integer, dimension(3,*) iad_crktg,
integer, dimension(*) iel_crktg,
integer, dimension(*) inod_crk,
integer, dimension(2,*) elcutc,
integer, dimension(2,*) nodedge,
integer, dimension(*) crknodiad,
integer, dimension(*) knod2elc,
x,
type (xfem_edge_), dimension(*) crkedge,
integer, dimension(3,*) xedge3n )

Definition at line 36 of file enrichtg_ini.F.

40C-----------------------------------------------
41 USE crackxfem_mod
42 USE elbufdef_mod
43 use element_mod , only : nixtg
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48C-----------------------------------------------
49C G l o b a l P a r a m e t e r s
50C-----------------------------------------------
51#include "mvsiz_p.inc"
52C-----------------------------------------------
53C C o m m o n B l o c k s
54C-----------------------------------------------
55#include "com_xfem1.inc"
56#include "units_c.inc"
57C-----------------------------------------------
58C D u m m y A r g u m e n t s
59C-----------------------------------------------
60 INTEGER IXTG(NIXTG,*),NFT,JFT,JLT,INOD_CRK(*),KNOD2ELC(*),
61 . IAD_CRKTG(3,*),IEL_CRKTG(*),ELCUTC(2,*),NODEDGE(2,*),
62 . CRKNODIAD(*),XEDGE3N(3,*),NXLAY
63 my_real
64 . x(3,*)
65 TYPE (ELBUF_STRUCT_) :: ELBUF_STR
66 TYPE (XFEM_EDGE_) , DIMENSION(*) :: CRKEDGE
67C-----------------------------------------------
68C L o c a l V a r i a b l e s
69C-----------------------------------------------
70 INTEGER I,K,ELCRK,ELCRKTG,ELCUT,pp1,pp2,pp3,IADC(3),IENR0(3),
71 . IENR(3),IED,IEDGE,IE10,IE20,IE1,IE2,NOD1,NOD2,N(3),DX(6),
72 . NX(3),dd(3),LAYCUT,ISIGN1,ISIGN2,ISIGN3,SIGBETA,
73 . ISIGN0(NXEL,3),p1,p2,
74 . NTAG(3),EDGEENR(3),ENR(3),ILAY,ITRI,NX1,NX2,NX3,NM,NP
75 my_real
76 . x1g(mvsiz),x2g(mvsiz),x3g(mvsiz),y1g(mvsiz),y2g(mvsiz),
77 . y3g(mvsiz),z1g(mvsiz),z2g(mvsiz),z3g(mvsiz),area(mvsiz),
78 . lxyz0(2),rx(mvsiz),ry(mvsiz),rz(mvsiz),
79 . sx(mvsiz),sy(mvsiz),sz(mvsiz),r11(mvsiz),r12(mvsiz),
80 . r13(mvsiz),r21(mvsiz),r22(mvsiz),r23(mvsiz),r31(mvsiz),
81 . r32(mvsiz),r33(mvsiz),xl1(mvsiz),yl1(mvsiz),xl2(mvsiz),
82 . yl2(mvsiz),xl3(mvsiz),yl3(mvsiz),
83 . fit(3,mvsiz),offg(mvsiz),xin(2,mvsiz),yin(2,mvsiz),fi,
84 . xxl(3,mvsiz),yyl(3,mvsiz),xn(3),yn(3),xm(2),ym(2)
85 my_real x1,y1,x2,y2,x3,y3,
86 . beta,area1,area2,area3
87 DATA dd/2,3,1/
88 DATA dx/1,2,3,1,2,3/
89C=======================================================================
90C COORDONNEES (GLOBAL)
91C-----------------------
92 p1 = 0
93 p2 = 0
94 nx1 = 0
95 nx2 = 0
96 area1 = 0
97 area2 = 0
98 area3 = 0
99 DO i=jft,jlt
100 x1g(i)=x(1,ixtg(2,i+nft))
101 y1g(i)=x(2,ixtg(2,i+nft))
102 z1g(i)=x(3,ixtg(2,i+nft))
103 x2g(i)=x(1,ixtg(3,i+nft))
104 y2g(i)=x(2,ixtg(3,i+nft))
105 z2g(i)=x(3,ixtg(3,i+nft))
106 x3g(i)=x(1,ixtg(4,i+nft))
107 y3g(i)=x(2,ixtg(4,i+nft))
108 z3g(i)=x(3,ixtg(4,i+nft))
109 ENDDO
110C----------------------------
111C LOCAL SYSTEM
112C----------------------------
113 DO i=jft,jlt
114 rx(i) = x2g(i)-x1g(i)
115 ry(i) = y2g(i)-y1g(i)
116 rz(i) = z2g(i)-z1g(i)
117 sx(i) = x3g(i)-x1g(i)
118 sy(i) = y3g(i)-y1g(i)
119 sz(i) = z3g(i)-z1g(i)
120 offg(i) = elbuf_str%GBUF%OFF(i)
121 ENDDO
122 k = 0
123 CALL clskew3(jft,jlt,k,
124 . rx, ry, rz,
125 . sx, sy, sz,
126 . r11,r12,r13,r21,r22,r23,r31,r32,r33,area,offg )
127C--------------------------
128C GLOBAL-->LOCAL TRANSFORMATION
129C--------------------------
130 DO i=jft,jlt
131 xl1(i) = zero
132 yl1(i) = zero
133 xl2(i)=r11(i)*rx(i)+r21(i)*ry(i)+r31(i)*rz(i)
134 yl2(i)=r12(i)*rx(i)+r22(i)*ry(i)+r32(i)*rz(i)
135 xl3(i)=r11(i)*sx(i)+r21(i)*sy(i)+r31(i)*sz(i)
136 yl3(i)=r12(i)*sx(i)+r22(i)*sy(i)+r32(i)*sz(i)
137 area(i)=half*area(i)
138 ENDDO
139C
140C LOCAL TO CENTER OF ELEM (possible warped element)
141C
142 DO i=jft,jlt
143 lxyz0(1)=third*(xl1(i)+xl2(i)+xl3(i))
144 lxyz0(2)=third*(yl1(i)+yl2(i)+yl3(i))
145C
146 xl1(i)=xl1(i)-lxyz0(1)
147 yl1(i)=yl1(i)-lxyz0(2)
148 xl2(i)=xl2(i)-lxyz0(1)
149 yl2(i)=yl2(i)-lxyz0(2)
150 xl3(i)=xl3(i)-lxyz0(1)
151 yl3(i)=yl3(i)-lxyz0(2)
152 ENDDO
153c-----------------------------------------------------
154c Loop over layers
155c-----------------------------------------------------
156 DO ilay=1,nxlay
157 pp1 = nxel*(ilay-1)+1
158 pp2 = pp1 + 1
159 pp3 = pp2 + 1
160C
161 DO i=jft,jlt
162 fit(1,i)=zero
163 fit(2,i)=zero
164 fit(3,i)=zero
165 ENDDO
166C
167 DO i=jft,jlt
168 elcrktg = iel_crktg(i+nft)
169 elcrk = elcrktg + ecrkxfec
170 laycut = crkedge(ilay)%LAYCUT(elcrk)
171 IF (laycut /= 0) THEN
172 xn(1) = xl1(i)
173 yn(1) = yl1(i)
174 xn(2) = xl2(i)
175 yn(2) = yl2(i)
176 xn(3) = xl3(i)
177 yn(3) = yl3(i)
178 xxl(1,i)=xl1(i)
179 yyl(1,i)=yl1(i)
180 xxl(2,i)=xl2(i)
181 yyl(2,i)=yl2(i)
182 xxl(3,i)=xl3(i)
183 yyl(3,i)=yl3(i)
184 DO k=1,3
185 p1 = k
186 p2 = dd(k)
187 ied = crkedge(ilay)%IEDGETG(k,elcrktg)
188 IF (ied > 0) THEN
189 iedge = xedge3n(k,elcrktg)
190 beta = crkedge(ilay)%RATIO(iedge)
191 nod1 = nodedge(1,iedge)
192 nod2 = nodedge(2,iedge)
193 IF (nod1 == ixtg(k+1,i+nft) .and. nod2 == ixtg(dd(k)+1,i+nft)) THEN
194 p1 = k
195 p2 = dd(k)
196 ELSEIF (nod2 == ixtg(k+1,i+nft).and.nod1 == ixtg(dd(k)+1,i+nft))THEN
197 p1 = dd(k)
198 p2 = k
199 ENDIF
200 xin(ied,i) = xn(p1) + beta*(xn(p2) - xn(p1))
201 yin(ied,i) = yn(p1) + beta*(yn(p2) - yn(p1))
202 xm(ied) = half*(xn(p1)+xn(p2))
203 ym(ied) = half*(yn(p1)+yn(p2))
204 ENDIF
205 ENDDO
206C
207 DO k=1,3
208 fi = zero
209 CALL lsint4(xm(1),ym(1),xm(2),ym(2),xn(k),yn(k),fi )
210 IF (fit(k,i)==zero) fit(k,i) = fi
211 ENDDO
212 ENDIF
213 ENDDO
214C
215 IF (ilay == 1) THEN
216 DO i=jft,jlt
217 elcrktg = iel_crktg(i+nft)
218 elcrk = elcrktg + ecrkxfec
219 elcut = crkedge(ilay)%LAYCUT(elcrk)
220 IF (elcut /= 0) THEN
221 elcutc(1,i+nft) = 1
222 elcutc(2,i+nft) = 1
223 ENDIF
224 ENDDO
225 ENDIF
226c--------------------------------------------------------------------
227 DO i=jft,jlt
228 elcrktg = iel_crktg(i+nft)
229 elcrk = elcrktg + ecrkxfec
230 laycut = crkedge(ilay)%LAYCUT(elcrk)
231 IF (laycut /= 0) THEN
232c
233 iadc(1) = iad_crktg(1,elcrktg)
234 iadc(2) = iad_crktg(2,elcrktg)
235 iadc(3) = iad_crktg(3,elcrktg)
236C
237 ienr0(1) = crknodiad(iadc(1))
238 ienr0(2) = crknodiad(iadc(2))
239 ienr0(3) = crknodiad(iadc(3))
240C
241 n(1) = ixtg(2,i+nft)
242 n(2) = ixtg(3,i+nft)
243 n(3) = ixtg(4,i+nft)
244C
245 nx(1) = inod_crk(n(1))
246 nx(2) = inod_crk(n(2))
247 nx(3) = inod_crk(n(3))
248C
249 ienr(1) = ienr0(1) + knod2elc(nx(1))*(ilay-1)
250 ienr(2) = ienr0(2) + knod2elc(nx(2))*(ilay-1)
251 ienr(3) = ienr0(3) + knod2elc(nx(3))*(ilay-1)
252C
253 ntag(1:3) = 0
254 edgeenr(1:3) = 0
255 enr(1:3) = 0
256C
257 DO k=1,3
258 ied = crkedge(ilay)%IEDGETG(k,elcrktg)
259 IF (ied > 0) THEN
260 ntag(k) = ntag(k) + 1
261 ntag(dd(k)) = ntag(dd(k)) + 1
262 iedge = xedge3n(k,elcrktg)
263 nod1 = nodedge(1,iedge)
264 nod2 = nodedge(2,iedge)
265 ie10 = crkedge(ilay)%EDGEENR(1,iedge)
266 ie20 = crkedge(ilay)%EDGEENR(2,iedge)
267 IF (nod1 == n(k) .and. nod2 == n(dd(k))) THEN
268 p1 = k
269 p2 = dd(k)
270 ELSEIF (nod2 == n(k) .and. nod1 == n(dd(k))) THEN
271 p1 = dd(k)
272 p2 = k
273 ENDIF
274 edgeenr(p1) = ie10
275 edgeenr(p2) = ie20
276 ENDIF
277 ENDDO
278C
279 DO k=1,3
280 IF (ntag(k) > 0) THEN
281 enr(k) = edgeenr(k)
282 ELSE
283 enr(k) = ienr(k)
284 ENDIF
285 ENDDO
286C
287 DO k=1,3
288 IF (ienr(k) > ienrnod) THEN
289 WRITE(iout,*) 'ERROR CRACK INITIATION,ENRICHMENT NODE EXCEEDED'
290 CALL arret(2)
291 ENDIF
292 ENDDO
293C
294 isign1 = int(sign(one,fit(1,i)))
295 isign2 = int(sign(one,fit(2,i)))
296 isign3 = int(sign(one,fit(3,i)))
297C
298 IF (fit(1,i) == zero) isign1 = 0
299 IF (fit(2,i) == zero) isign2 = 0
300 IF (fit(3,i) == zero) isign3 = 0
301C
302 DO k=1,nxel
303 isign0(k,1) = isign1
304 isign0(k,2) = isign2
305 isign0(k,3) = isign3
306 ENDDO
307c------------------------------------------
308 itri = 0
309 nm = 0
310 np = 0
311 DO k=1,3
312 IF (isign0(1,k) > 0) THEN
313 itri = itri + 1
314 np = k
315 ELSEIF (isign0(1,k) < 0) THEN
316 nm = k
317 ENDIF
318 ENDDO
319 IF (itri == 1) THEN
320 itri = -1
321 nx1 = np
322 ELSEIF (itri == 2) THEN
323 itri = 1
324 nx1 = nm
325 ENDIF
326 nx2 = dx(nx1+1)
327 nx3 = dx(nx2+1)
328 xfem_phantom(ilay)%ITRI(1,elcrk) = itri
329 xfem_phantom(ilay)%ITRI(2,elcrk) = nx1
330C
331 xfem_phantom(ilay)%IFI(iadc(1)) = isign0(1,1)
332 xfem_phantom(ilay)%IFI(iadc(2)) = isign0(1,2)
333 xfem_phantom(ilay)%IFI(iadc(3)) = isign0(1,3)
334C-------------------
335c IXEL = 1 => positif element within ILAY (ILEV = PP1)
336C-------------------
337 crklvset(pp1)%ENR0(1,iadc(1)) = enr(1)
338 crklvset(pp1)%ENR0(1,iadc(2)) = enr(2)
339 crklvset(pp1)%ENR0(1,iadc(3)) = enr(3)
340C
341 IF (isign0(1,1) > 0) crklvset(pp1)%ENR0(1,iadc(1)) = 0
342 IF (isign0(1,2) > 0) crklvset(pp1)%ENR0(1,iadc(2)) = 0
343 IF (isign0(1,3) > 0) crklvset(pp1)%ENR0(1,iadc(3)) = 0
344C
345 crklvset(pp1)%ENR0(2,iadc(1)) = crklvset(pp1)%ENR0(1,iadc(1))
346 crklvset(pp1)%ENR0(2,iadc(2)) = crklvset(pp1)%ENR0(1,iadc(2))
347 crklvset(pp1)%ENR0(2,iadc(3)) = crklvset(pp1)%ENR0(1,iadc(3))
348C-------------------
349c IXEL = 2 => negatif element within ILAY (ILEV = PP2)
350C-------------------
351 crklvset(pp2)%ENR0(1,iadc(1)) = enr(1)
352 crklvset(pp2)%ENR0(1,iadc(2)) = enr(2)
353 crklvset(pp2)%ENR0(1,iadc(3)) = enr(3)
354C
355 IF (isign0(2,1) < 0) crklvset(pp2)%ENR0(1,iadc(1)) = 0
356 IF (isign0(2,2) < 0) crklvset(pp2)%ENR0(1,iadc(2)) = 0
357 IF (isign0(2,3) < 0) crklvset(pp2)%ENR0(1,iadc(3)) = 0
358C
359 crklvset(pp2)%ENR0(2,iadc(1)) = crklvset(pp2)%ENR0(1,iadc(1))
360 crklvset(pp2)%ENR0(2,iadc(2)) = crklvset(pp2)%ENR0(1,iadc(2))
361 crklvset(pp2)%ENR0(2,iadc(3)) = crklvset(pp2)%ENR0(1,iadc(3))
362C-------------------
363c IXEL = 3 => Third phantom element (ILEV = PP3)
364C-------------------
365 IF (itri < 0) THEN ! sign ILEV3 = ILEV2 < 0
366 ie1 = xedge3n(nx1,elcrktg)
367 ie2 = xedge3n(nx3,elcrktg)
368c IF (CRKEDGE(ILAY)%ICUTEDGE(IE2) == 2) THEN ! tip edge
369 crklvset(pp3)%ENR0(1,iadc(nx1)) = abs(crklvset(pp2)%ENR0(1,iadc(nx1)))
370 crklvset(pp3)%ENR0(1,iadc(nx2)) = crklvset(pp2)%ENR0(1,iadc(nx2))
371 crklvset(pp3)%ENR0(1,iadc(nx3)) = crklvset(pp2)%ENR0(1,iadc(nx3))
372 crklvset(pp2)%ENR0(1,iadc(nx1)) = -crknodiad(iadc(nx1)) - knod2elc(nx(nx1))*(ilay-1)
373c
374 beta = crkedge(ilay)%RATIO(ie2)
375 nod1 = nodedge(1,ie2)
376 nod2 = nodedge(2,ie2)
377 k = nx3 ! starting edge node (local)
378 IF (nod1 == ixtg(k+1,i+nft) .and.
379 . nod2 == ixtg(dd(k)+1,i+nft)) THEN
380 p1 = k
381 p2 = dd(k)
382 sigbeta = ie2
383 ELSEIF (nod2 == ixtg(k+1,i+nft) .and.
384 . nod1 == ixtg(dd(k)+1,i+nft)) THEN
385 p1 = dd(k)
386 p2 = k
387 sigbeta = -ie2
388 ENDIF
389c
390c-- AREA factors
391 x1 = xxl(nx1,i)
392 y1 = yyl(nx1,i)
393 ied = crkedge(ilay)%IEDGETG(nx1,elcrktg)
394 x2 = xin(ied,i)
395 y2 = yin(ied,i)
396 ied = crkedge(ilay)%IEDGETG(nx3,elcrktg)
397 x3 = xin(ied,i)
398 y3 = yin(ied,i)
399 area1 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
400 area1 = area1 / area(i)
401 x1 = xxl(nx2,i)
402 y1 = yyl(nx2,i)
403 x2 = xxl(nx3,i)
404 y2 = yyl(nx3,i)
405 area2 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
406 area3 = one - area1 - area2
407c
408 ELSEIF (itri > 0) THEN ! sign ILEV3 = ILEV1 > 0
409c
410 ie1 = xedge3n(nx1,elcrktg) ! tip edge
411c IF (CRKEDGE(ILAY)%ICUTEDGE(IE1) == 2) THEN
412 crklvset(pp3)%ENR0(1,iadc(nx1)) = abs(crklvset(pp1)%ENR0(1,iadc(nx1)))
413 crklvset(pp3)%ENR0(1,iadc(nx2)) = crklvset(pp1)%ENR0(1,iadc(nx2))
414 crklvset(pp3)%ENR0(1,iadc(nx3)) = crklvset(pp1)%ENR0(1,iadc(nx3))
415 crklvset(pp1)%ENR0(1,iadc(nx1)) = -crknodiad(iadc(nx1)) - knod2elc(nx(nx1))*(ilay-1)
416c
417 beta = crkedge(ilay)%RATIO(ie1)
418 nod1 = nodedge(1,ie1)
419 nod2 = nodedge(2,ie1)
420 k = nx1 ! starting edge node (local)
421 IF (nod1 == ixtg(k+1,i+nft) .and.
422 . nod2 == ixtg(dd(k)+1,i+nft)) THEN
423 p1 = k
424 p2 = dd(k)
425 sigbeta = ie1
426 ELSEIF (nod2 == ixtg(k+1,i+nft) .and.
427 . nod1 == ixtg(dd(k)+1,i+nft)) THEN
428 p1 = dd(k)
429 p2 = k
430 sigbeta = -ie1
431 ENDIF
432c
433c-- AREA factors for each phantom
434 x1 = xxl(nx1,i)
435 y1 = yyl(nx1,i)
436 ied = crkedge(ilay)%IEDGETG(nx3,elcrktg)
437 x2 = xin(ied,i)
438 y2 = yin(ied,i)
439 ied = crkedge(ilay)%IEDGETG(nx1,elcrktg)
440 x3 = xin(ied,i)
441 y3 = yin(ied,i)
442 area1 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
443 area1 = area1 / area(i)
444 x1 = xxl(nx2,i)
445 y1 = yyl(nx2,i)
446 x3 = xxl(nx3,i)
447 y3 = yyl(nx3,i)
448 area2 = half*abs((x1-x3)*(y2-y1) - (x1-x2)*(y3-y1))
449 area3 = one - area1 - area2
450c
451 ENDIF ! ITRI
452c
453 crklvset(pp3)%ENR0(2,iadc(1)) = crklvset(pp3)%ENR0(1,iadc(1))
454 crklvset(pp3)%ENR0(2,iadc(2)) = crklvset(pp3)%ENR0(1,iadc(2))
455 crklvset(pp3)%ENR0(2,iadc(3)) = crklvset(pp3)%ENR0(1,iadc(3))
456c
457 crklvset(pp1)%AREA(elcrk) = area1
458 crklvset(pp2)%AREA(elcrk) = area2
459 crklvset(pp3)%AREA(elcrk) = area3
460C-------------------
461 ENDIF ! LAYCUT /= 0
462 ENDDO ! DO I=JFT,JLT
463 ENDDO ! DO ILAY=1,NXLAY
464C-------------
465 RETURN
clskew3
subroutine clskew3(jft, jlt, irep, rx, ry, rz, sx, sy, sz, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z, det)
Definition clskew.F:34
my_real
#define my_real
Definition cppsort.cpp:32
lsint4
subroutine lsint4(y1, z1, y2, z2, y, z, fi)
Definition crklayer4n_adv.F:712
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
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:86