31
32
33
35
36
37
38#include "implicit_f.inc"
39
40
41
42 INTEGER IVOLU(*), (*), NN
44 . rvolu(*), x(3,*)
45
46
47
48 INTEGER I, II, NBX, NBY, NBRIC, NNB, IFV, NNB1, J, K, NV, NALL,
49 . , TNC(4), TNV(5,6), KK, L, N1, N2, N3, N4, N5, N6, N7,
50 . N8
52 . vx3, vy3, vz3, vx1, vy1, vz1,
norm, ss, vx2, vy2, vz2,
53 . x0, y0, z0, zlmin, zlmax, xx, yy, zz, zl, lx, ly, dx, dy,
54 . xlc, ylc, zlc1, zlc2, xxg, yyg, zzg, xc, yc, zc, l1, l2,
55 . l3, x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4,
56 . nx1, ny1, nz1, nx2, ny2, nz2, area1, area2, nx, ny, nz,
57 . rr, xcf, ycf, zcf, vx, vy, vz
58 INTEGER, DIMENSION(:), ALLOCATABLE :: ITAG
59C
60 INTEGER FAC(4,6)
61 DATA fac /1,4,3,2,
62 . 5,6,7,8,
63 . 1,2,6,5,
64 . 2,3,7,6,
65 . 3,4,8,7,
66 . 4,1,5,8/
67
68 vx3=rvolu(35)
69 vy3=rvolu(36)
70 vz3=rvolu(37)
71 vx1=rvolu(38)
72 vy1=rvolu(39)
73 vz1=rvolu(40)
74
75 norm=sqrt(vx3**2+vy3**2+vz3**2)
79 ss=vx3*vx1+vy3*vy1+vz3*vz1
80 vx1=vx1-ss*vx3
81 vy1=vy1-ss*vy3
82 vz1=vz1-ss*vz3
83 norm=sqrt(vx1**2+vy1**2+vz1**2)
87 vx2=vy3*vz1-vz3*vy1
88 vy2=vz3*vx1-vx3*vz1
89 vz2=vx3*vy1-vy3*vx1
90
91 x0=rvolu(41)
92 y0=rvolu(42)
93 z0=rvolu(43)
94
95 zlmin=ep30
96 zlmax=-ep30
97 DO i=1,nn
98 ii=ibuf(i)
99 xx=x(1,ii)
100 yy=x(2,ii)
101 zz=x(3,ii)
102 zl=(xx-x0)*vx3+(yy-y0)*vy3+(zz-z0)*vz3
105 ENDDO
106
107 ifv=ivolu(45)
108 nbx=ivolu(54)
109 nby=ivolu(55)
110
111 lx=rvolu(44)
112 ly=rvolu(45)
113 dx=two*lx/nbx
114 dy=two*ly/nby
115 zlc1=zlmin-
max(lx,ly)
116 zlc2=zlmax+
max(lx,ly)
117 nnb=0
118 nnb1=(nbx+1)*(nby+1)
119 DO i=1,nby+1
120 ylc=-ly+(i-1)*dy
121 DO j=1,nbx+1
122 xlc=-lx+(j-1)*dx
123 nnb=nnb+1
126 fvdata(ifv)%XB(3,nnb)=zlc1
127 fvdata(ifv)%XB(1,nnb1+nnb)=xlc
128 fvdata(ifv)%XB(2,nnb1+nnb)=ylc
129 fvdata(ifv)%XB(3,nnb1+nnb)=zlc2
130 ENDDO
131 ENDDO
132
133 nnb=(nbx+1)*(nby+1)*2
134 DO i=1,nnb
138 xxg=x0+xx*vx1+yy*vx2+zz*vx3
139 yyg=y0+xx*vy1+yy*vy2+zz*vy3
140 zzg=z0+xx*vz1+yy*vz2+zz*vz3
144 ENDDO
145
146 nbric=0
147 DO i=1,nby
148 DO j=1,nbx
149 nbric=nbric+1
150 fvdata(ifv)%BRIC(1,nbric)=(i-1)*(nbx+1)+j
151 fvdata(ifv)%BRIC(2,nbric)=(i-1)*(nbx+1)+j+1
152 fvdata(ifv)%BRIC(3,nbric)=i*(nbx+1)+j+1
153 fvdata(ifv)%BRIC(4,nbric)=i*(nbx+1)+j
154 fvdata(ifv)%BRIC(5,nbric)=nnb1+(i-1)*(nbx+1)+j
155 fvdata(ifv)%BRIC(6,nbric)=nnb1+(i-1)*(nbx+1)+j+1
156 fvdata(ifv)%BRIC(7,nbric)=nnb1+i*(nbx+1)+j+1
157 fvdata(ifv)%BRIC(8,nbric)=nnb1+i*(nbx+1)+j
158 ENDDO
159 ENDDO
160
161 DO i=1,nbric
162 DO j=1,6
164 DO k=1,4
165 fvdata(ifv)%SFAC(j,k,i)=zero
166 ENDDO
167 ENDDO
168 DO j=1,7
169 fvdata(ifv)%TBRIC(6+j,i)=0
170 ENDDO
171 ENDDO
172
173 ALLOCATE(itag(nnb))
174 DO i=1,nnb
175 itag(i)=0
176 ENDDO
177
178 DO i=1,nbric
179 DO j=1,8
180 itag(
fvdata(ifv)%BRIC(j,i))=1
181 ENDDO
182 nv=0
183 DO j=1,nbric
184 IF (j==i) cycle
185 nall=0
186 nc=0
187 DO k=1,8
188 IF (itag(
fvdata(ifv)%BRIC(k,j))==1)
THEN
189 nc=nc+1
190 tnc(nc)=
fvdata(ifv)%BRIC(k,j)
191 ENDIF
192 nall=nall+itag(
fvdata(ifv)%BRIC(k,j))
193 ENDDO
194 IF (nall/=4) cycle
195 nv=nv+1
196 tnv(1,nv)=j
197 tnv(2,nv)=tnc(1)
198 tnv(3,nv)=tnc(2)
199 tnv(4,nv)=tnc(3)
200 tnv(5,nv)=tnc(4)
201 ENDDO
202 DO j=1,8
203 itag(
fvdata(ifv)%BRIC(j,i))=0
204 ENDDO
205 DO j=1,6
206 DO k=1,4
207 kk=fac(k,j)
208 itag(
fvdata(ifv)%BRIC(kk,i))=1
209 ENDDO
210 DO k=1,nv
211 nall=1
212 DO l=1,4
213 nall=nall*itag(tnv(1+l,k))
214 ENDDO
215 IF (nall==1)
fvdata(ifv)%TBRIC(j,i)=tnv(1,k)
216 ENDDO
217 DO k=1,4
218 kk=fac(k,j)
219 itag(
fvdata(ifv)%BRIC(kk,i))=0
220 ENDDO
221 ENDDO
222 ENDDO
223 DEALLOCATE(itag)
224
225 DO i=1,nbric
234 xc=one_over_8*(
fvdata(ifv)%XB(1,n1)+
fvdata(ifv)%XB(1,n2)+
238 yc=one_over_8*(
fvdata(ifv)%XB(2,n1)+
fvdata(ifv)%XB(2,n2)+
242 zc=one_over_8*(
fvdata(ifv)%XB(3,n1)+
fvdata(ifv)%XB(3,n2)+
246
247 IF (i==1) THEN
257 l1=sqrt(vx1**2+vy1**2+vz1**2)
258 l2=sqrt(vx2**2+vy2**2+vz2**2)
259 l3=sqrt(vx3**2+vy3**2+vz3**2)
261 ENDIF
262
263 DO j=1,6
264 n1=
fvdata(ifv)%BRIC(fac(1,j),i)
265 n2=
fvdata(ifv)%BRIC(fac(2,j),i)
266 n3=
fvdata(ifv)%BRIC(fac(3,j),i)
267 n4=
fvdata(ifv)%BRIC(fac(4,j),i)
280
281 vx1=x2-x1
282 vy1=y2-y1
283 vz1=z2-z1
284 vx2=x3-x1
285 vy2=y3-y1
286 vz2=z3-z1
287 vx3=x4-x1
288 vy3=y4-y1
289 vz3=z4-z1
290 nx1=vy1*vz2-vz1*vy2
291 ny1=vz1*vx2-vx1*vz2
292 nz1=vx1*vy2-vy1*vx2
293 nx2=vy2*vz3-vz2*vy3
294 ny2=vz2*vx3-vx2*vz3
295 nz2=vx2*vy3-vy2*vx3
296 area1=half*sqrt(nx1**2+ny1**2+nz1**2)
297 area2=half*sqrt(nx2**2+ny2**2+nz2**2)
298 fvdata(ifv)%SFAC(j,1,i)=area1+area2
299
300 nx=half*(nx1+nx2)
301 ny=half*(ny1+ny2)
302 nz=half*(nz1+nz2)
303 rr=sqrt(nx**2+ny**2+nz**2)
304 xcf=fourth*(x1+x2+x3+x4)
305 ycf=fourth*(y1+y2+y3+y4)
306 zcf=fourth*(z1+z2+z3+z4)
307 vx=xc-xcf
308 vy=yc-ycf
309 vz=zc-zcf
310 ss=vx*nx+vy*ny+vz*nz
311 IF(rr == zero) cycle
312 IF (ss<=zero) THEN
313 fvdata(ifv)%SFAC(j,2,i)=nx/rr
314 fvdata(ifv)%SFAC(j,3,i)=ny/rr
315 fvdata(ifv)%SFAC(j,4,i)=nz/rr
316 ELSE
317 fvdata(ifv)%SFAC(j,2,i)=-nx/rr
318 fvdata(ifv)%SFAC(j,3,i)=-ny/rr
319 fvdata(ifv)%SFAC(j,4,i)=-nz/rr
320 ENDIF
321 ENDDO
322 ENDDO
323
324 RETURN
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
type(fvbag_data), dimension(:), allocatable fvdata