42
43
44
46
47
48
49#include "implicit_f.inc"
50
51
52
53#include "mvsiz_p.inc"
54
55
56
57#include "param_c.inc"
58
59
60
61 INTEGER LEDGE(NLEDGE,*), IRECT(4,*), CAND_M(*), CAND_S(*), ADMSR(4,*),
62 . LBOUND(*), JLT, NEDGE, IEDGE, ITAB(*), IGAP0, IGAP,
63 . N1(MVSIZ), N2(MVSIZ),
64 . M1(MVSIZ), M2(MVSIZ)
65
67 . x(3,*),
68 . xxs1(mvsiz), xxs2(mvsiz), xys1(mvsiz), xys2(mvsiz),
69 . xzs1(mvsiz), xzs2(mvsiz), xxm1(mvsiz), xxm2(mvsiz),
70 . xym1(mvsiz), xym2(mvsiz), xzm1(mvsiz), xzm2(mvsiz),
71 . gape(*),gapve(mvsiz),gap_e_l(nedge),
72 . ex(mvsiz), ey(mvsiz), ez(mvsiz), fx(mvsiz), fy(mvsiz), fz(mvsiz)
73 real*4 edg_bisector(3,4,*), vtx_bisector(3,2,*)
74
75
76
77 INTEGER I ,NN, J, JRM, K, KRM, I1, J1, I2, J2,
78 . IE, JE, SOL_EDGE, SH_EDGE, IM(MVSIZ), IS(MVSIZ)
80 . aaa, dx, dy, dz, dd, nni, ni2, invcos, gape_m(mvsiz), gape_s(mvsiz)
81
82 DO i=1,jlt
83 IF(cand_s(i)<=nedge) THEN
84
85 i1=ledge(1,cand_s(i))
86 j1=ledge(2,cand_s(i))
87 n1(i)=ledge(5,cand_s(i))
88 n2(i)=ledge(6,cand_s(i))
89
90 IF(irect(j1,i1)==n1(i).AND.irect(mod(j1,4)+1,i1)==n2(i))THEN
91 is(i)= 1
92 ELSEIF(irect(j1,i1)==n2(i).AND.irect(mod(j1,4)+1,i1)==n1(i))THEN
93 is(i)=-1
94 ELSE
95 print *,'i25cor3e - internal problem',cand_s(i),n1(i),n2(i),
96 . irect(j1,i1),irect(mod(j1,4)+1,i1)
97 END IF
98
99 i2=ledge(1,cand_m(i))
100 j2=ledge(2,cand_m(i))
101 m1(i)=ledge(5,cand_m(i))
102 m2(i)=ledge(6,cand_m(i))
103
104 IF(irect(j2,i2)==m1(i).AND.irect(mod(j2,4)+1,i2)==m2(i))THEN
105 im(i)= 1
106 ELSEIF(irect(j2,i2)==m2(i).AND.irect(mod(j2,4)+1,i2)==m1(i))THEN
107 im(i)=-1
108 ELSE
109 print *,'i25cor3e - internal problem',cand_m(i),m1(i),m2(i),
110 . irect(j2,i2),irect(mod(j2,4)+1,i2)
111 END IF
112
113 xxs1(i) = x(1,n1(i))
114 xys1(i) = x(2,n1(i))
115 xzs1(i) = x(3,n1(i))
116 xxs2(i) = x(1,n2(i))
117 xys2(i) = x(2,n2(i))
118 xzs2(i) = x(3,n2(i))
119 xxm1(i) = x(1,m1(i))
120 xym1(i) = x(2,m1(i))
121 xzm1(i) = x(3,m1(i))
122 xxm2(i) = x(1,m2(i))
123 xym2(i) = x(2,m2(i))
124 xzm2(i) = x(3,m2(i))
125 END IF
126 END DO
127
128 DO i=1,jlt
129 gape_m(i)=gape(cand_m(i))
130 IF(cand_s(i)<=nedge) THEN
131 gape_s(i)=gape(cand_s(i))
132 ELSE
133 END IF
134 gapve(i)=gape_m(i)+gape_s(i)
135 END DO
136 IF(igap==3) THEN
137 DO i=1,jlt
138 gape_m(i)=
min(gape_m(i),gap_e_l(cand_m(i)))
139 IF(cand_s(i)<=nedge) THEN
140 gape_s(i)=
min(gape_s(i),gap_e_l(cand_s(i)))
141 gapve(i)=
min(gap_e_l(cand_m(i))+gap_e_l(cand_s(i)),gapve(i))
142 ENDIF
143 ENDDO
144 ENDIF
145
146 sol_edge=iedge/10
147 sh_edge =iedge-10*sol_edge
148
149 ex(1:jlt)=zero
150 ey(1:jlt)=zero
151 ez(1:jlt)=zero
152 fx(1:jlt)=zero
153 fy(1:jlt)=zero
154 fz(1:jlt)=zero
155 IF(sh_edge/=0)THEN
156 DO i=1,jlt
157
158 ie=ledge(1,cand_m(i))
159 je=ledge(2,cand_m(i))
160 ex(i) = edg_bisector(1,je,ie)
161 ey(i) = edg_bisector(2,je,ie)
162 ez(i) = edg_bisector(3,je,ie)
163
164 IF(iabs(ledge(7,cand_s(i)))/=1)THEN
165 ie=ledge(1,cand_s(i))
166 je=ledge(2,cand_s(i))
167 fx(i) = edg_bisector(1,je,ie)
168 fy(i) = edg_bisector(2,je,ie)
169 fz(i) = edg_bisector(3,je,ie)
170 END IF
171
172 END DO
173 END IF
174
175 DO i=1,jlt
176
177 IF(ledge(3,cand_m(i))/=0)cycle
178
179 ie=ledge(1,cand_m(i))
180 je=ledge(2,cand_m(i))
181 IF(im(i)==1)THEN
182 i1=admsr(je,ie)
183 i2=admsr(mod(je,4)+1,ie)
184 ELSE
185 i2=admsr(je,ie)
186 i1=admsr(mod(je,4)+1,ie)
187 END IF
188
189 ex(i) = edg_bisector(1,je,ie)
190 ey(i) = edg_bisector(2,je,ie)
191 ez(i) = edg_bisector(3,je,ie)
192
193
194 dx = vtx_bisector(1,1,i1)+vtx_bisector(1,2,i1)
195 dy = vtx_bisector(2,1,i1)+vtx_bisector(2,2,i1)
196 dz = vtx_bisector(3,1,i1)+vtx_bisector(3,2,i1)
197
198 nni = ex(i)*dx + ey(i)*dy + ez(i)*dz
199 ni2 = dx*dx + dy*dy + dz*dz
200
201 IF(nni < zero)THEN
202 dx=dx-two*nni*ex(i)
203 dy=dy-two*nni*ey(i)
204 dz=dz-two*nni*ez(i)
205 nni=-nni
206 END IF
207
208 IF(two*nni*nni < ni2)THEN
209
210 aaa = sqrt(
max(zero,ni2-nni*nni)) - nni
211 dx = dx + aaa*ex(i)
212 dy = dy + aaa*ey(i)
213 dz = dz + aaa*ez(i)
214 ENDIF
215 dd=one/
max(em20,sqrt(dx*dx+dy*dy+dz*dz))
216 dx = dx*dd
217 dy = dy*dd
218 dz = dz*dd
219 invcos = one /
max(em20,ex(i)*dx + ey(i)*dy + ez(i)*dz)
220 dx = dx*invcos
221 dy = dy*invcos
222 dz = dz*invcos
223
224 xxm1(i) = xxm1(i)-gape_m(i)*dx
225 xym1(i) = xym1(i)-gape_m(i)*dy
226 xzm1(i) = xzm1(i)-gape_m(i)*dz
227
228
229 dx = vtx_bisector(1,1,i2)+vtx_bisector(1,2,i2)
230 dy = vtx_bisector(2,1,i2)+vtx_bisector(2,2,i2)
231 dz = vtx_bisector(3,1,i2)+vtx_bisector(3,2,i2)
232
233 nni = ex(i)*dx + ey(i)*dy + ez(i)*dz
234 ni2 = dx*dx + dy*dy + dz*dz
235
236 IF(nni < zero)THEN
237 dx=dx-two*nni*ex(i)
238 dy=dy-two*nni*ey(i)
239 dz=dz-two*nni*ez(i)
240 nni=-nni
241 END IF
242
243 IF(two*nni*nni < ni2)THEN
244
245 aaa = sqrt(
max(zero,ni2-nni*nni)) - nni
246 dx = dx + aaa*ex(i)
247 dy = dy + aaa*ey(i)
248 dz = dz + aaa*ez(i)
249 ENDIF
250 dd=one/
max(em20,sqrt(dx*dx+dy*dy+dz*dz))
251 dx = dx*dd
252 dy = dy*dd
253 dz = dz*dd
254 invcos = one /
max(em20,ex(i)*dx + ey(i)*dy + ez(i)*dz)
255 dx = dx*invcos
256 dy = dy*invcos
257 dz = dz*invcos
258
259 xxm2(i) = xxm2(i)-gape_m(i)*dx
260 xym2(i) = xym2(i)-gape_m(i)*dy
261 xzm2(i) = xzm2(i)-gape_m(i)*dz
262
263 END DO
264
265 IF(igap0/=0) THEN
266 DO i=1,jlt
267
268 IF(ledge(3,cand_s(i))/=0)cycle
269
270 ie=ledge(1,cand_s(i))
271 je=ledge(2,cand_s(i))
272 IF(is(i)==1)THEN
273 i1=admsr(je,ie)
274 i2=admsr(mod(je,4)+1,ie)
275 ELSE
276 i2=admsr(je,ie)
277 i1=admsr(mod(je,4)+1,ie)
278 END IF
279
280 fx(i) = edg_bisector(1,je,ie)
281 fy(i) = edg_bisector(2,je,ie)
282 fz(i) = edg_bisector(3,je,ie)
283
284
285 dx = vtx_bisector(1,1,i1)+vtx_bisector(1,2,i1)
286 dy = vtx_bisector(2,1,i1)+vtx_bisector(2,2,i1)
287 dz = vtx_bisector(3,1,i1)+vtx_bisector(3,2,i1)
288
289 nni = fx(i)*dx + fy(i)*dy + fz(i)*dz
290 ni2 = dx*dx + dy*dy + dz*dz
291
292 IF(nni < zero)THEN
293 dx=dx-two*nni*fx(i)
294 dy=dy-two*nni*fy(i)
295 dz=dz-two*nni*fz(i)
296 nni=-nni
297 END IF
298
299 IF(two*nni*nni < ni2)THEN
300
301 aaa = sqrt(
max(zero,ni2-nni*nni)) - nni
302 dx = dx + aaa*fx(i)
303 dy = dy + aaa*fy(i)
304 dz = dz + aaa*fz(i)
305 ENDIF
306 dd=one/
max(em20,sqrt(dx*dx+dy*dy+dz*dz))
307 dx = dx*dd
308 dy = dy*dd
309 dz = dz*dd
310 invcos = one /
max(em20,fx(i)*dx + fy(i)*dy + fz(i)*dz)
311 dx = dx*invcos
312 dy = dy*invcos
313 dz = dz*invcos
314
315 xxs1(i) = xxs1(i)-gape_s(i)*dx
316 xys1(i) = xys1(i)-gape_s(i)*dy
317 xzs1(i) = xzs1(i)-gape_s(i)*dz
318
319
320 dx = vtx_bisector(1,1,i2)+vtx_bisector(1,2,i2)
321 dy = vtx_bisector(2,1,i2)+vtx_bisector(2,2,i2)
322 dz = vtx_bisector(3,1,i2)+vtx_bisector(3,2,i2)
323
324 nni = fx(i)*dx + fy(i)*dy + fz(i)*dz
325 ni2 = dx*dx + dy*dy + dz*dz
326
327 IF(nni < zero)THEN
328 dx=dx-two*nni*fx(i)
329 dy=dy-two*nni*fy(i)
330 dz=dz-two*nni*fz(i)
331 nni=-nni
332 END IF
333
334 IF(two*nni*nni < ni2)THEN
335
336 aaa = sqrt(
max(zero,ni2-nni*nni)) - nni
337 dx = dx + aaa*fx(i)
338 dy = dy + aaa*fy(i)
339 dz = dz + aaa*fz(i)
340 ENDIF
341 dd=one/
max(em20,sqrt(dx*dx+dy*dy+dz*dz))
342 dx = dx*dd
343 dy = dy*dd
344 dz = dz*dd
345 invcos = one /
max(em20,fx(i)*dx + fy(i)*dy + fz(i)*dz)
346 dx = dx*invcos
347 dy = dy*invcos
348 dz = dz*invcos
349
350 xxs2(i) = xxs2(i)-gape_s(i)*dx
351 xys2(i) = xys2(i)-gape_s(i)*dy
352 xzs2(i) = xzs2(i)-gape_s(i)*dz
353
354 END DO
355 END IF
356
357 RETURN