43
44
45
46#include "implicit_f.inc"
47
48
49
50#include "units_c.inc"
51
52
53
54 INTEGER NNS, NNTR, NPOLH
55 INTEGER IBUF(*), IBUFA(*), ELEMA(3,*), TAGELA(*),
56 . IVOLU(*), IBPOLH(*), IFVNOD(3,*),IFVTRI(6,*),
57 . IFVPOLY(*),IFVTADR(*),IFVPOLH(*), IFVPADR(*)
58 my_real x(3,*), rfvnod(2,*), dlh
59
60
61
62 INTEGER I, J, K, I1, I2, IEL, JJ, KK,
63 . N1, N2, N3, NN1, NN2, NN3,
64 . NADBR
65 INTEGER IDBR1, IDBR2, IDBR3
67 . x1, y1, z1, x2, y2, z2, x3, y3, z3,
68 . nx, ny, nz, area2, ksi, eta,
area, fac,
69 . areapoly, areapolymax,
70 . pnod(3,nns), pvolu(npolh),
71 . parea(nntr), pnorm(3,nntr)
73 . x12, y12, z12, x23, y23, z23, x31, y31, z31,
74 . l12, l23, l31, dl0, dl1, dl2, dll, dlbr1, dlbr2, dlbr3
75
76
77
78 x12 = zero
79 y12 = zero
80 z12 = zero
81 x23 = zero
82 y23 = zero
83 z23 = zero
84 x31 = zero
85 y31 = zero
86 z31 = zero
87 dl0 = dlh
88 dl1 = ep30
89 dl2 = ep30
90 dlbr1 = ep30
91 dlbr2 = ep30
92 dlbr3 = ep30
93 x1 = zero
94 y1 = zero
95 z1 = zero
96 x2 = zero
97 y2 = zero
98 z2 = zero
99 x3 = zero
100 y3 = zero
101 z3 = zero
102
103
104 DO i=1,nns
105 IF (ifvnod(1,i)==1) THEN
106 iel=ifvnod(2,i)
107 ksi=rfvnod(1,i)
108 eta=rfvnod(2,i)
109
110 n1=elema(1,iel)
111 n2=elema(2,iel)
112 n3=elema(3,iel)
113
114 IF (tagela(iel)>0) THEN
115 nn1=ibuf(n1)
116 nn2=ibuf(n2)
117 nn3=ibuf(n3)
118 x1=x(1,nn1)
119 y1=x(2,nn1)
120 z1=x(3,nn1)
121 x2=x(1,nn2)
122 y2=x(2,nn2)
123 z2=x(3,nn2)
124 x3=x(1,nn3)
125 y3=x(2,nn3)
126 z3=x(3,nn3)
127 ELSEIF (tagela(iel)<0) THEN
128 nn1=ibufa(n1)
129 nn2=ibufa(n2)
130 nn3=ibufa(n3)
131 x1=x(1,nn1)
132 y1=x(2,nn1)
133 z1=x(3,nn1)
134 x2=x(1,nn2)
135 y2=x(2,nn2)
136 z2=x(3,nn2)
137 x3=x(1,nn3)
138 y3=x(2,nn3)
139 z3=x(3,nn3)
140 ENDIF
141 pnod(1,i)=(one-ksi-eta)*x1+ksi*x2+eta*x3
142 pnod(2,i)=(one-ksi-eta)*y1+ksi*y2+eta*y3
143 pnod(3,i)=(one-ksi-eta)*z1+ksi*z2+eta*z3
144
145 ELSEIF (ifvnod(1,i)==2) THEN
146 i2=ifvnod(2,i)
147 pnod(1,i)=x(1,i2)
148 pnod(2,i)=x(2,i2)
149 pnod(3,i)=x(3,i2)
150 ENDIF
151 ENDDO
152
153 DO i=1,nns
154 IF (ifvnod(1,i)==3) THEN
155 i1=ifvnod(2,i)
156 i2=ifvnod(3,i)
157 fac=rfvnod(1,i)
158 pnod(1,i)=fac*pnod(1,i1)+(one-fac)*pnod(1,i2)
159 pnod(2,i)=fac*pnod(2,i1)+(one-fac)*pnod(2,i2)
160 pnod(3,i)=fac*pnod(3,i1)+(one-fac)*pnod(3,i2)
161 ENDIF
162 ENDDO
163
164
165
166 DO i=1,nntr
167 n1=ifvtri(1,i)
168 n2=ifvtri(2,i)
169 n3=ifvtri(3,i)
170 CALL fvnormal(pnod,n1,n2,n3,0,nx,ny,nz)
171 area2=sqrt(nx*nx+ny*ny+nz*nz)
172 parea(i)=half*area2
173 IF (area2>0) THEN
174 pnorm(1,i)=nx/area2
175 pnorm(2,i)=ny/area2
176 pnorm(3,i)=nz/area2
177 ELSE
178 pnorm(1,i)=zero
179 pnorm(2,i)=zero
180 pnorm(3,i)=zero
181 ENDIF
182 ENDDO
183
184
185
186 dl0=dlh
187 dl1=ep30
188 dl2=ep30
189 dlbr1=ep30
190 dlbr2=ep30
191 dlbr3=ep30
192 nadbr=0
193 DO i=1,npolh
194 pvolu(i)=zero
195 areapolymax=zero
196
197 DO j=ifvpadr(i),ifvpadr(i+1)-1
198 jj=ifvpolh(j)
199 areapoly=zero
200
201 DO k=ifvtadr(jj), ifvtadr(jj+1)-1
202 kk=ifvpoly(k)
204 areapoly=areapoly+
area
205 iel=ifvtri(4,kk)
206 IF (iel>0) THEN
207 nx=pnorm(1,kk)
208 ny=pnorm(2,kk)
209 nz=pnorm(3,kk)
210 ELSE
211 IF (ifvtri(5,kk)==i) THEN
212 nx=pnorm(1,kk)
213 ny=pnorm(2,kk)
214 nz=pnorm(3,kk)
215 ELSEIF (ifvtri(6,kk)==i) THEN
216 nx=-pnorm(1,kk)
217 ny=-pnorm(2,kk)
218 nz=-pnorm(3,kk)
219 ENDIF
220 ENDIF
221 n1=ifvtri(1,kk)
222 x1=pnod(1,n1)
223 y1=pnod(2,n1)
224 z1=pnod(3,n1)
225 pvolu(i)=pvolu(i)+third*
area*(x1*nx+y1*ny+z1*nz)
226
227 n2=ifvtri(2,kk)
228 x2=pnod(1,n2)
229 y2=pnod(2,n2)
230 z2=pnod(3,n2)
231 n3=ifvtri(3,kk)
232 x3=pnod(1,n3)
233 y3=pnod(2,n3)
234 z3=pnod(3,n3)
235 x12=x2-x1
236 y12=y2-y1
237 z12=z2-z1
238 x23=x3-x2
239 y23=y3-y2
240 z23=z3-z2
241 x31=x1-x3
242 y31=y1-y3
243 z31=z1-z3
244 l12=x12**2+y12**2+z12**2
245 l23=x23**2+y23**2+z23**2
246 l31=x31**2+y31**2+z31**2
247 IF(ibpolh(i)==0) THEN
248 dl2=
min(dl2,l12,l23,l31)
249 ELSE
251 IF(dll < dlbr2) THEN
252 dlbr2=dll
253 idbr2=i
254 ENDIF
255 ENDIF
256 ENDDO
257 areapolymax=
max(areapolymax,areapoly)
258 ENDDO
259
260 IF(ibpolh(i)==0) THEN
261 dl1=
min(dl1,pvolu(i))
262 ELSE
263 nadbr=nadbr+1
264 IF(pvolu(i) < dlbr1) THEN
265 dlbr1=pvolu(i)
266 idbr1=i
267 ENDIF
268 dll=pvolu(i)/areapolymax
269 IF(dll < dlbr3) THEN
270 dlbr3=dll
271 idbr3=i
272 ENDIF
273
274 IF(pvolu(i)<0)THEN
275 WRITE(iout,'(A,E12.4,3I10)') 'NEGATIVE VOLUME',
276 . pvolu(i),i,ibpolh(i),nadbr
277 ENDIF
278 ENDIF
279 ENDDO
280
281 IF(dl1 > zero) THEN
282 dl1=dl1**third
283 ELSE
284 dl1=zero
285 ENDIF
286 dl2=sqrt(dl2)
287 IF(dlbr1 > zero) THEN
288 dlbr1=dlbr1**third
289 ELSE
290 dlbr1=zero
291 ENDIF
292 dlbr2=sqrt(dlbr2)
293 IF(dl0==zero) dlh=dlbr2
294
295
296
297 WRITE(iout,1000) ivolu(1),npolh,npolh-nadbr,dl0,dl1,dl2
298 IF(nadbr > 0) THEN
299 WRITE(iout,1001) nadbr,
300 . dlbr1,idbr1,iabs(ibpolh(idbr1)),
301 . dlbr2,idbr2,iabs(ibpolh(idbr2)),
302 . dlbr3,idbr3,iabs(ibpolh(idbr3))
303 ENDIF
304
3051000 FORMAT(
306 . //' FVMBAG: FINITE VOLUME MINIMUM LENGTH '/
307 . ' -------------------------------------'/
308 . /5x,'VOLUME NUMBER ',i10,
309 . /5x,'TOTAL NUMBER OF FINITE VOLUMES.. . . . . . .=',i10,
310 . /5x,'NUMBER OF POLYHEDRA . . . . .. . . . . . . .=',i10,
311 . /5x,' MINIMUM LENGTH USED FOR TIME STEP. . . .=',1pg20.13,
312 . /5x,' MINIMUM LENGTH BASED ON VOLUME . . . . .=',1pg20.13,
313 . /5x,' MINIMUM LENGTH BASED ON NODAL DISTANCE .=',1pg20.13)
3141001 FORMAT(
315 . 5x,'NUMBER OF ADDITIONAL BRICKS. . . . . . . . .=',i10,
316 . /5x,' MINIMUM LENGTH BASED ON VOLUME . . . . .=',1pg20.13,' (FINITE VOLUME ID ',i10,')',' (BRICK ID ',i10,')',
317 . /5x,' MINIMUM LENGTH BASED ON NODAL DISTANCE .=',1pg20.13,' (FINITE VOLUME ID ',i10,')',' (BRICK ID ',i10,')',
318 . /5x,' MINIMUM LENGTH BASED ON VOLUME/AREA. . .=',1pg20.13,' (FINITE VOLUME ID ',i10,')',' (BRICK ID ',i10,')',/)
319
320 RETURN
subroutine fvnormal(x, n1, n2, n3, n4, nx, ny, nz)
subroutine area(d1, x, x2, y, y2, eint, stif0)