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