37
38
39
41 USE elbufdef_mod
42 use element_mod , only : nixtg
43
44
45
46#include "implicit_f.inc"
47
48
49
50#include "mvsiz_p.inc"
51
52
53
54#include "com_xfem1.inc"
55
56
57
58 INTEGER LFT,LLT,NFT,NXLAY
59 INTEGER IELCRKTG(*),EDGETG(3,*),IEDGESH3(3,*),ELCUT(*),
60 . XNOD(2,2),IXTG(NIXTG,*),NODEDGE(2,*),TAGSKYTG(3,*),KNOD2ELC(*),
61 . TAGEDGE(*)
63 . x1l(*),y1l(*),x2l(*),y2l(*),x3l(*),y3l(*),beta0(2)
64
65 TYPE (ELBUF_STRUCT_), TARGET :: ELBUF_STR
66 TYPE (ELBUF_STRUCT_), DIMENSION(NXEL) , TARGET :: XFEM_STR
67 TYPE (XFEM_LVSET_) , DIMENSION(NLEVMAX) :: CRKLVSET
68 TYPE (XFEM_SHELL_) , DIMENSION(NLEVMAX) :: CRKSHELL
69 TYPE (XFEM_EDGE_) , DIMENSION(NXLAYMAX) :: CRKEDGE
70 TYPE (XFEM_PHANTOM_), DIMENSION(NXLAYMAX) :: XFEM_PHANTOM
71
72
73
74 INTEGER I,K,II,IL,ILAY,ELCRK,IED,IEDGE,ICRK,p1,p2,
75 . NOD1,NOD2,JCRK,IXEL
76 INTEGER dd(3),d1(3),d2(3),IFI(2),ILEV(NXEL),N(3),ISIGN0(3),
77 . IENR0(3),IENR(3),NTAG(3)
79 . fit(3,mvsiz),xn(3),yn(3),xmi(2),ymi(2),beta
80 . off_phantom
82 EXTERNAL lsintx
83 DATA dd/2,3,1/
84 DATA d1/2,3,4/
85 DATA d2/3,4,2/
86
87 TYPE(G_BUFEL_) , POINTER :: GBUF
88 TYPE(L_BUFEL_) , POINTER :: LBUF
89
90 p2 = 0
91 DO i=lft,llt
92 xn(1)=x1l(i)
93 yn(1)=y1l(i)
94 xn(2)=x2l(i)
95 yn(2)=y2l(i)
96 xn(3)=x3l(i)
97 yn(3)=y3l(i)
98 IF (elcut(i+nft) > 0) THEN
99 DO k=1,3
100 p1 = k
101 p2 = dd(k)
102 ied = edgetg(k,i+nft)
103 IF (ied > 0) THEN
104 xmi(ied) = half*(xn(p1) + xn(p2))
105 ymi(ied) = half*(yn(p1) + yn(p2))
106 ENDIF
107 ENDDO
108
109 DO k=1,3
110 fit(k,i) = lsintx(xmi(1),ymi(1),xmi(2),ymi(2),xn(k),yn(k))
111 ENDDO
112 ENDIF
113 ENDDO
114
115 DO i=lft,llt
116 elcrk = ielcrktg(i+nft)
117 beta(1,i) = zero
118 beta(2,i) = zero
119 IF(elcut(i+nft) > 0)THEN
120
121 jcrk = elcrk - ecrkxfec
122 DO k=1,3
123 iedge = iedgesh3(k,jcrk)
124 ied = edgetg(k,i+nft)
125 IF (ied > 0) THEN
126 nod1 = nodedge(1,iedge)
127 nod2 = nodedge(2,iedge)
128 IF (nod1 == xnod(ied,1) .and. nod2 == xnod(ied,2)) THEN
129 beta(ied,i) = beta0(ied)
130 ELSE IF (nod2 == xnod(ied,1) .and. nod1 == xnod(ied,2)) THEN
131 beta(ied,i) = one - beta0(ied)
132 END IF
133 ENDIF
134 ENDDO
135 ENDIF
136 ENDDO
137
138 DO ilay=1,nxlay
139 ii = nxel*(ilay-1)
140 DO k=1,nxel
141 ilev(k) = ii + k
142 ENDDO
143 DO i=lft,llt
144 elcrk = ielcrktg(i+nft)
145 jcrk = elcrk - ecrkxfec
146 IF (elcut(i+nft) > 0) THEN
147 icrk = crkshell(ilev(1))%PHANTOMG(elcrk)
148 crklvset(ilev(1))%ELCUT(elcrk) = icrk
149 crklvset(ilev(2))%ELCUT(elcrk) = -icrk
150
151 xfem_phantom(ilay)%ELCUT(elcrk) = icrk
152 crkedge(ilay)%LAYCUT(elcrk) = 2
153
154 n(1) = ixtg(2,i+nft)
155 n(2) = ixtg(3,i+nft)
156 n(3) = ixtg(4,i+nft)
157
158 isign0(1) = int(sign(one,fit(1,i))) * icrk
159 isign0(2) = int(sign(one,fit(2,i))) * icrk
160 isign0(3) = int(sign(one,fit(3,i))) * icrk
161
162 ntag(1:3) = 0
163
164 DO k=1,3
165 ienr0(k) = 0
166 ienr(k) = 0
167 ied = edgetg(k,i+nft)
168 IF (ied > 0) THEN
169 ntag(k) = ntag(k) + 1
170 ntag(dd(k)) = ntag(dd(k)) + 1
171 ENDIF
172 ENDDO
173
174 DO k=1,3
175 ied = edgetg(k,i+nft)
176 iedge = iedgesh3(k,jcrk)
177 IF(ied > 0)THEN
178 nod1 = nodedge(1,iedge)
179 nod2 = nodedge(2,iedge)
180 IF(nod1 == n(k) .and. nod2 == n(dd(k)))THEN
181 p1 = k
182 p2 = dd(k)
183 ELSE IF(nod2 == n(k) .and. nod1 == n(dd(k)))THEN
184 p1 = dd(k)
185 p2 = k
186 END IF
187 IF(ntag(p1) > 0.AND.crkedge(ilay)%EDGEENR(1,iedge) > 0)
188 . ienr0(p1) = crkedge(ilay)%EDGEENR(1,iedge)
189 IF(ntag(p2) > 0.AND.crkedge(ilay)%EDGEENR(2,iedge) > 0)
190 . ienr0(p2) = crkedge(ilay)%EDGEENR(2,iedge)
191 ENDIF
192 ENDDO
193
194 DO k=1,3
195 IF(ienr0(k) /= 0)THEN
196 ienr(k) = ienr0(k)
197 ELSE
198 ienr(k) = tagskytg(k,i+nft)+knod2elc(n(k))*(ilay-1)
199 ENDIF
200 ENDDO
201
202 DO k=1,3
203 ied = edgetg(k,i+nft)
204 iedge = iedgesh3(k,jcrk)
205 IF(ied > 0)THEN
206 DO il=1,nxel
207 crklvset(ilev(il))%EDGETG(k,jcrk) = ied
208 crklvset(ilev(il))%ICUTEDGE(iedge) = 1
209 crklvset(ilev(il))%RATIOEDGE(iedge) = beta(ied,i)
210 ENDDO
211
212 crkedge(ilay)%EDGETIP(1,iedge) =
max(ied,
213 . crkedge(ilay)%EDGETIP(1,iedge))
214 crkedge(ilay)%EDGETIP(2,iedge) =
215 . crkedge(ilay)%EDGETIP(2,iedge) + 1
216
217
218
219 IF(crkedge(ilay)%EDGEICRK(iedge) == 0)
220 . crkedge(ilay)%EDGEICRK(iedge) = icrk
221
222 nod1 = nodedge(1,iedge)
223 nod2 = nodedge(2,iedge)
224 ifi(1:2) = 0
225 p1 = 0
226 p2 = 0
227 IF(nod1 == n(k) .and. nod2 == n(dd(k)))THEN
228 ifi(1) = isign0(k)
229 ifi(2) = isign0(dd(k))
230 p1 = k
231 p2 = dd(k)
232 ELSE IF(nod2 == n(k) .and. nod1 == n(dd(k)))THEN
233 ifi(1) = isign0(dd(k))
234 ifi(2) = isign0(k)
235 p1 = dd(k)
236 p2 = k
237 END IF
238 IF(crkedge(ilay)%EDGEIFI(1,iedge) == 0)
239 . crkedge(ilay)%EDGEIFI(1,iedge) = ifi(1)
240 IF(crkedge(ilay)%EDGEIFI(2,iedge) == 0)
241 . crkedge(ilay)%EDGEIFI(2,iedge) = ifi(2)
242
243 IF(crkedge(ilay)%EDGEENR(1,iedge) == 0)
244 . crkedge(ilay)%EDGEENR(1,iedge) = ienr(p1)
245
246 IF(crkedge(ilay)%EDGEENR(2,iedge) == 0)
247 . crkedge(ilay)%EDGEENR(2,iedge) = ienr(p2)
248 ENDIF
249 ENDDO
250 ENDIF
251 ENDDO
252 ENDDO
253
254
255
256 IF (nxlay > 1)THEN
257
258 DO ixel=1,2
259 DO ilay=1,nxlay
260 lbuf => xfem_str(ixel)%BUFLY(ilay)%LBUF(1,1,1)
261 DO i=lft,llt
262 IF(elcut(i+nft) > 0)THEN
263 off_phantom = lbuf%OFF(i)
264 lbuf%OFF(i) = - off_phantom
265 ENDIF
266 ENDDO
267 ENDDO
268 ENDDO
269 ELSE
270
271 DO ixel=1,2
272 gbuf => xfem_str(ixel)%GBUF
273 DO i=lft,llt
274 IF(elcut(i+nft) > 0)THEN
275 off_phantom = gbuf%OFF(i)
276 gbuf%OFF(i) = - off_phantom
277 ENDIF
278 ENDDO
279 ENDDO
280 ENDIF
281
282
283
284 DO i=lft,llt
285 IF(elcut(i+nft) > 0)THEN
286 elbuf_str%GBUF%OFF(i) = zero
287 ENDIF
288 ENDDO
289
290 DO i=lft,llt
291 elcrk = ielcrktg(i+nft) - ecrkxfec
292 IF (elcut(i+nft) > 0) THEN
293 DO k=1,3
294 ied = edgetg(k,i+nft)
295 iedge = iedgesh3(k,elcrk)
296 IF (ied > 0 .and. iedge > 0) THEN
297 tagedge(iedge) = tagedge(iedge) + 1
298 ENDIF
299 ENDDO
300 ENDIF
301 ENDDO
302
303 RETURN