32 use element_mod , only : nixs
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48#include "implicit_f.inc"
49
50
51
52#include "com04_c.inc"
53
54
55
57 . x0(3),vn(3),xyz0(3,*),xyz(3,*),d(3,*),al(*)
58 INTEGER IXS(NIXS,*),NUMTOT,NUMEL,NUMCON,NC(5,*),NVOIS(2,*),NA(*),
59 . ITYP,NODCUT
60
61
62
63 INTEGER IARETE(2,12),NEWL,NUMNEW,I,J,K,L,I1,I2,N1,N2,N3,KK,
64 . NODEL,NVC(2,24),NAC(24),NCA(24),NCB(24),LMIN,KB
65
67 . unm,distmin,distmax,dist0,tol,dist(8),
68 . xc(24),yc(24),zc(24),alc(24),tet(24),dis,xm,ym,zm,
69 . x1,y1,z1,v1,xi,yi,zi,vi,xi1,yi1,zi1,csi,ssi,tetmin,
70 . x12,y12,z12,x23,y23,z23,v12,x2,y2,z2
71 DATA iarete/1,2,2,3,3,4,4,1,1,5,2,6,3,7,4,8,5,6,6,7,7,8,8,5/
72 unm = -one
73
74 newl=0
75 numnew=0
76 numcon=0
77 distmin=ep30
78 dist0=vn(1)*x0(1)+vn(2)*x0(2)+vn(3)*x0(3)
79
80 do209 i=1,numels
81 i1=ixs(2,i)
82 i2=ixs(3,i)
83 x1=abs(xyz0(1,i1)-xyz0(1,i2))+abs(xyz0(2,i1)-xyz0(2,i2))
84 & +abs(xyz0(3,i1)-xyz0(3,i2))
85 209 distmin=
min(distmin,abs(x1))
86 tol=em3*distmin
87
88 do100 i=1,numels
89 distmin= ep30
90 distmax=-ep30
91 nodel=0
92 do110 j=1,8
93 k=ixs(j+1,i)
94 xi=xyz0(1,k)
95 yi=xyz0(2,k)
96 zi=xyz0(3,k)
97 IF(ityp==2)THEN
98 xi=xi-d(1,k)
99 yi=yi-d(2,k)
100 zi=zi-d(3,k)
101 ENDIF
102 dist(j)=xi*vn(1)+yi*vn(2)+zi*vn(3)
103 dist(j)=dist(j)-dist0
104 distmin=
min(dist(j),distmin)
105 110 distmax=
max(dist(j),distmax)
106
107 IF(distmin*distmax> zero)goto100
108
109
110
111 DO 120 k=1,12
112 n1=iarete(1,k)
113 n2=iarete(2,k)
114 IF(dist(n1)*dist(n2)> zero) GOTO 120
115 x1 = xyz0(1,ixs(n1+1,i))
116 y1 = xyz0(2,ixs(n1+1,i))
117 z1 = xyz0(3,ixs(n1+1,i))
118 x2 = xyz0(1,ixs(n2+1,i))
119 y2 = xyz0(2,ixs(n2+1,i))
120 z2 = xyz0(3,ixs(n2+1,i))
121 IF(ityp==2)THEN
122 x1 =x1 - d(1,ixs(n1+1,i))
123 y1 =y1 - d(2,ixs(n1+1,i))
124 z1 =z1 - d(3,ixs(n1+1,i))
125 x2 =x2 - d(1,ixs(n2+1,i))
126 y2 =y2 - d(2,ixs(n2+1,i))
127 z2 =z2 - d(3,ixs(n2+1,i))
128 ENDIF
129 IF(abs(dist(n1)-dist(n2))<tol)THEN
130 nodel=nodel+1
131 xc(nodel) = x1
132 yc(nodel) = y1
133 zc(nodel) = z1
134 nvc(1,nodel)=ixs(n1+1,i)
135 nvc(2,nodel)=ixs(n1+1,i)
136 alc(nodel)=1.
137 nodel=nodel+1
138 xc(nodel) = x2
139 yc(nodel) = y2
140 zc(nodel) = z2
141 nvc(1,nodel)=ixs(n2+1,i)
142 nvc(2,nodel)=ixs(n2+1,i)
143 alc(nodel)=1.
144 ELSE
145 nodel=nodel+1
146 alc(nodel)=dist(n1)/(dist(n1)-dist(n2))
147 nvc(1,nodel)=ixs(n1+1,i)
148 nvc(2,nodel)=ixs(n2+1,i)
149 xc(nodel) = alc(nodel) *x2
150 & +(1-alc(nodel))*x1
151 yc(nodel) = alc(nodel) *y2
152 & +(1-alc(nodel))*y1
153 zc(nodel) = alc(nodel) *z2
154 & +(1-alc(nodel))*z1
155 ENDIF
156 120 CONTINUE
157
158
159
160
161 IF(nodel>2)THEN
162 k=1
163 nac(1)=1
164 DO 124 l=2,nodel
165 DO 125 j=1,l-1
166 dis=abs(xc(l)-xc(j))+
167 & abs(yc(l)-yc(j))+
168 & abs(zc(l)-zc(j))
169 IF(dis<=tol)THEN
170 nac(l)=nac(j)
171 GOTO 124
172 ENDIF
173 125 CONTINUE
174 k=k+1
175 nac(l)=k
176 124 CONTINUE
177 DO 126 l=1,nodel
178 xc(nac(l))=xc(l)
179 yc(nac(l))=yc(l)
180 zc(nac(l))=zc(l)
181 alc(nac(l))=alc(l)
182 nvc(1,nac(l))=nvc(1,l)
183 nvc(2,nac(l))=nvc(2,l)
184 126 CONTINUE
185 nodel=k
186 ENDIF
187
188 IF(nodel<=2)goto100
189 IF(nodel>6)goto100
190
191
192
193 xm=zero
194 ym=zero
195 zm=zero
196 DO 130 k=1,nodel
197 xm=xm+xc(k)/float(nodel)
198 ym=ym+yc(k)/float(nodel)
199 130 zm=zm+zc(k)/float(nodel)
200
201 x1=xc(1)-xm
202 y1=yc(1)-ym
203 z1=zc(1)-zm
204 v1=sqrt(x1**2+y1**2+z1**2)
205 IF(v1<tol) goto100
206 x1=x1/v1
207 y1=y1/v1
208 z1=z1/v1
209 tet(1)=zero
210
211 DO 140 k=2,nodel
212 xi=xc(k)-xm
213 yi=yc(k)-ym
214 zi=zc(k)-zm
215 vi=sqrt(xi**2+yi**2+zi**2)
216 IF(vi<tol) goto100
217 xi=xi/vi
218 yi=yi/vi
219 zi=zi/vi
220 csi=x1*xi+y1*yi+z1*zi
223
224 xi1=y1*zi-yi*z1
225 yi1=z1*xi-zi*x1
226 zi1=x1*yi-xi*y1
227 ssi=xi1*vn(1)+yi1*vn(2)+zi1*vn(3)
230 140 tet(k)=atan2(ssi,csi)
231
232 DO 150 k=1,nodel
233 tetmin=ep30
234 DO 151 l=1,nodel
235 IF(tet(l)<tetmin)THEN
236 lmin=l
237 tetmin=tet(l)
238 ENDIF
239 151 CONTINUE
240 nca(k)=lmin
241 tet(lmin)=ep30
242 150 CONTINUE
243
244
245
246
247 kb=0
248 DO 155 k=1,nodel
249 n1=nca(nodel)
250 IF(k>1) n1=nca(k-1)
251 n2=nca(k)
252 n3=nca(1)
253 IF(k<nodel)n3=nca(k+1)
254
255 x12=xc(n2)-xc(n1)
256 y12=yc(n2)-yc(n1)
257 z12=zc(n2)-zc(n1)
258 x23=xc(n3)-xc(n2)
259 y23=yc(n3)-yc(n2)
260 z23=zc(n3)-zc(n2)
261 v12=sqrt(x12**2+y12**2+z12**2)*sqrt(x23**2+y23**2+z23**2)
262 xi1=(y12*z23-y23*z12)/v12
263 yi1=(z12*x23-z23*x12)/v12
264 zi1=(x12*y23-x23*y12)/v12
265 ssi=xi1*vn(1)+yi1*vn(2)+zi1*vn(3)
266 IF(ssi>em30)THEN
267 IF(kb==4)THEN
268 kb=kb+1
269 ncb(kb)=ncb(1)
270 kb=kb+1
271 ncb(kb)=ncb(kb-2)
272 ENDIF
273 kb=kb+1
274 ncb(kb)=nca(k)+numnew
275
276 ENDIF
277 155 CONTINUE
278
279 IF(kb==3.OR.kb==7)THEN
280 kb=kb+1
281 ncb(kb)=ncb(kb-1)
282
283 ENDIF
284 kb=int(kb/4)
285 DO k=1,kb
286 newl=newl+1
287 kk = (k-1)*4+1
288 nc(1,newl)= ncb(kk)
289 nc(2,newl)= ncb(kk+1)
290 nc(3,newl)= ncb(kk+2)
291 nc(4,newl)= ncb(kk+3)
292 nc(5,newl)= i
293 ENDDO
294
295
296
297 IF(ityp==2)THEN
298 do260 k=1,nodel
299 n1=nvc(1,k)
300 n2=nvc(2,k)
301 xc(k)=alc(k)*xyz0(1,n2)+(1-alc(k))*xyz0(1,n1)
302 yc(k)=alc(k)*xyz0(2,n2)+(1-alc(k))*xyz0(2,n1)
303 zc(k)=alc(k)*xyz0(3,n2)+(1-alc(k))*xyz0(3,n1)
304 260 CONTINUE
305 ENDIF
306 DO 270 k=1,nodel
307 xyz(1,numnew+k)=xc(k)
308 xyz(2,numnew+k)=yc(k)
309 xyz(3,numnew+k)=zc(k)
310 al(numnew+k)=alc(k)
311 nvois(1,numnew+k)=nvc(1,k)
312 nvois(2,numnew+k)=nvc(2,k)
313 270 CONTINUE
314
315 numnew=numnew+nodel
316
317 100 CONTINUE
318
319
320
321
322 k=1
323 na(1)=1
324 DO 1240 i=2,numnew
325 DO 1250 j=1,i-1
326 dis=abs(xyz(1,i)-xyz(1,j))+
327 & abs(xyz(2,i)-xyz(2,j))+
328 & abs(xyz(3,i)-xyz(3,j))
329 IF(dis<=tol)THEN
330 na(i)=na(j)
331 GOTO 1240
332 ENDIF
333 1250 CONTINUE
334 k=k+1
335 na(i)=k
336 1240 CONTINUE
337 numtot=k
338 numel=newl
339
340 DO 1260 i=1,numnew
341 al(na(i))=al(i)
342 nvois(1,na(i))=nvois(1,i)
343 nvois(2,na(i))=nvois(2,i)
344 xyz(1,na(i))=xyz(1,i)
345 xyz(2,na(i))=xyz(2,i)
346 xyz(3,na(i))=xyz(3,i)
347 1260 CONTINUE
348 DO k=1,numel
349 nc(1,k)=na(nc(1,k))+nodcut
350 nc(2,k)=na(nc(2,k))+nodcut
351 nc(3,k)=na(nc(3,k))+nodcut
352 nc(4,k)=na(nc(4,k))+nodcut
353 ENDDO
354 IF(numel==0)THEN
355 x1=0.
356 y1=-vn(3)
357 z1= vn(2)
358 v1=sqrt(y1**2+z1**2)
359 IF(v1>em10)THEN
360 v1=ep04*tol/v1
361 y1=y1*v1
362 z1=z1*v1
363 x2= vn(2)*z1-vn(3)*y1
364 y2= -vn(1)*z1
365 z2= vn(1)*y1
366 ELSE
367 x1=zero
368 y1=ep04*tol
369 z1=zero
370 x2=zero
371 y2=zero
372 z2=ep04*tol
373 ENDIF
374 numel=1
375 numtot=1
376 xyz(1,numtot)=x0(1)-x1-x2
377 xyz(2,numtot)=x0(2)-y1-y2
378 xyz(3,numtot)=x0(3)-z1-z2
379 nc(1,1)=numtot+nodcut
380 nvois(1,numtot)=1
381 nvois(2,numtot)=1
382 al(numtot)=zero
383 numtot=numtot+1
384 xyz(1,numtot)=x0(1)+x1-x2
385 xyz(2,numtot)=x0(2)+y1-y2
386 xyz(3,numtot)=x0(3)+z1-z2
387 nc(2,1)=numtot+nodcut
388 nvois(1,numtot)=1
389 nvois(2,numtot)=1
390 al(numtot)=0.
391 numtot=numtot+1
392 xyz(1,numtot)=x0(1)+x1+x2
393 xyz(2,numtot)=x0(2)+y1+y2
394 xyz(3,numtot)=x0(3)+z1+z2
395 nc(3,1)=numtot+nodcut
396 nvois(1,numtot)=1
397 nvois(2,numtot)=1
398 al(numtot)=0.
399 numtot=numtot+1
400 xyz(1,numtot)=x0(1)-x1+x2
401 xyz(2,numtot)=x0(2)-y1+y2
402 xyz(3,numtot)=x0(3)-z1+z2
403 nc(4,1)=numtot+nodcut
404 nvois(1,numtot)=1
405 nvois(2,numtot)=1
406 al(numtot)=zero
407 nc(5,1)=1
408 ENDIF
409 RETURN