CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
mafillv.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine mafillv (nef, ipnei, neifa, neiel, vfa, xxn, area, auv, adv, jq, irow, nzs, bv, vel, cosa, umfa, xlet, xle, gradvfa, xxi, body, volume, ielfa, lakonf, ifabou, nbody, dtimef, velo, veloo, sel, xrlfa, gamma, xxj, nactdohinv, a1, a2, a3, flux, nefa, nefb, icyclic, c, ifatie, iau6, xxni, xxnj, iturbulent, gradvel)
 

Function/Subroutine Documentation

◆ mafillv()

subroutine mafillv ( integer, intent(in)  nef,
integer, dimension(*), intent(in)  ipnei,
integer, dimension(*), intent(in)  neifa,
integer, dimension(*), intent(in)  neiel,
real*8, dimension(0:7,*), intent(in)  vfa,
real*8, dimension(3,*), intent(in)  xxn,
real*8, dimension(*), intent(in)  area,
real*8, dimension(*), intent(inout)  auv,
real*8, dimension(*), intent(inout)  adv,
integer, dimension(*), intent(in)  jq,
integer, dimension(*), intent(in)  irow,
integer, intent(in)  nzs,
real*8, dimension(nef,3), intent(inout)  bv,
real*8, dimension(nef,0:7), intent(in)  vel,
real*8, dimension(*), intent(in)  cosa,
real*8, dimension(*), intent(in)  umfa,
real*8, dimension(*), intent(in)  xlet,
real*8, dimension(*), intent(in)  xle,
real*8, dimension(3,3,*), intent(in)  gradvfa,
real*8, dimension(3,*), intent(in)  xxi,
real*8, dimension(0:3,*), intent(in)  body,
real*8, dimension(*), intent(in)  volume,
integer, dimension(4,*), intent(in)  ielfa,
character*8, dimension(*), intent(in)  lakonf,
integer, dimension(*), intent(in)  ifabou,
integer, intent(in)  nbody,
real*8, intent(in)  dtimef,
real*8, dimension(nef,0:7), intent(in)  velo,
real*8, dimension(nef,0:7), intent(in)  veloo,
real*8, dimension(3,*), intent(inout)  sel,
real*8, dimension(3,*), intent(in)  xrlfa,
real*8, dimension(*), intent(in)  gamma,
real*8, dimension(3,*), intent(in)  xxj,
integer, dimension(*), intent(in)  nactdohinv,
real*8, intent(in)  a1,
real*8, intent(in)  a2,
real*8, intent(in)  a3,
real*8, dimension(*), intent(in)  flux,
integer  nefa,
integer  nefb,
integer  icyclic,
real*8, dimension(3,3)  c,
integer, dimension(*)  ifatie,
integer, dimension(6,*)  iau6,
real*8, dimension(3,*)  xxni,
real*8, dimension(3,*)  xxnj,
integer  iturbulent,
real*8, dimension(3,3,*), intent(in)  gradvel 
)
25 !
26  implicit none
27 !
28  character*8 lakonf(*)
29 !
30  integer i,nef,indexf,ipnei(*),j,ifa,iel,neifa(*),icyclic,
31  & neiel(*),jq(*),irow(*),nzs,iwall,ielfa(4,*),nefa,nefb,
32  & ipointer,ifabou(*),nbody,k,indexb,nactdohinv(*),
33  & ifatie(*),iau6(6,*),iturbulent
34 !
35  real*8 xflux,vfa(0:7,*),xxn(3,*),area(*),auv(*),adv(*),bv(nef,3),
36  & vel(nef,0:7),cosa(*),umfa(*),xlet(*),xle(*),coef,gradvfa(3,3,*),
37  & xxi(3,*),body(0:3,*),volume(*),coef2,dtimef,velo(nef,0:7),
38  & veloo(nef,0:7),rhovel,constant,sel(3,*),xrlfa(3,*),gamma(*),
39  & xxj(3,*),a1,a2,a3,flux(*),c(3,3),xxni(3,*),xxnj(3,*),difcoef,
40  & xl1,xl2,aa,bb,gradvel(3,3,*)
41 !
42  intent(in) nef,ipnei,neifa,neiel,vfa,xxn,area,
43  & jq,irow,nzs,vel,cosa,umfa,xlet,xle,gradvfa,xxi,
44  & body,volume,ielfa,lakonf,ifabou,nbody,
45  & dtimef,velo,veloo,xrlfa,gamma,xxj,nactdohinv,a1,
46  & a2,a3,flux,gradvel
47 !
48  intent(inout) adv,auv,bv,sel
49 !
50  do i=nefa,nefb
51  do indexf=ipnei(i)+1,ipnei(i+1)
52 !
53 ! convection
54 !
55  ifa=neifa(indexf)
56  iel=neiel(indexf)
57  xflux=flux(indexf)
58 !
59 ! vfa sometimes undefined?
60 !
61  if(xflux.ge.0.d0) then
62 !
63 ! outflowing xflux
64 !
65  if(iel.eq.0) then
66  adv(i)=adv(i)+xflux
67  else
68  adv(i)=adv(i)+xflux
69  do k=1,3
70 c cd
71 c bv(i,k)=bv(i,k)-(vfa(k,ifa)-vel(i,k))*xflux
72 c retarded gamma
73  bv(i,k)=bv(i,k)-gamma(ifa)*(vfa(k,ifa)
74  & -vel(i,k))*xflux
75 c improved smart
76 c aa=gamma(ifa)
77 c bb=0.d0
78 c if(aa.lt.0.5d0) then
79 c bb=2.d0/3.d0
80 c elseif(aa.lt.0.8d0) then
81 c bb=3.d0/8.d0
82 c endif
83 c bv(i,k)=bv(i,k)-(aa*vel(i,k)
84 c & +(1.d0-aa)*vel(iel,k)
85 c & -2.d0*(1.d0-aa-bb)*
86 c & (gradvel(k,1,i)*xxj(1,indexf)+
87 c & gradvel(k,2,i)*xxj(2,indexf)+
88 c & gradvel(k,3,i)*xxj(3,indexf))*xlet(indexf)
89 c & -vel(i,k))
90  enddo
91  endif
92  else
93  if(iel.gt.0) then
94 !
95 ! incoming flux from neighboring element
96 !
97  if((icyclic.eq.0).or.(ifatie(ifa).eq.0)) then
98  auv(indexf)=auv(indexf)+xflux
99  do k=1,3
100 c cd
101 c bv(i,k)=bv(i,k)
102 c & -(vfa(k,ifa)-vel(iel,k))*xflux
103 c retarded gamma
104  bv(i,k)=bv(i,k)-gamma(ifa)*
105  & (vfa(k,ifa)-vel(iel,k))*xflux
106 c improved smart
107 c aa=gamma(ifa)
108 c bb=0.d0
109 c if(aa.lt.0.5d0) then
110 c bb=2.d0/3.d0
111 c elseif(aa.lt.0.8d0) then
112 c bb=3.d0/8.d0
113 c endif
114 c bv(i,k)=bv(i,k)-(aa*vel(iel,k)
115 c & +(1.d0-aa)*vel(i,k)
116 c & +2.d0*(1.d0-aa-bb)*
117 c & (gradvel(k,1,iel)*xxj(1,indexf)+
118 c & gradvel(k,2,iel)*xxj(2,indexf)+
119 c & gradvel(k,3,iel)*xxj(3,indexf))*xlet(indexf)
120 c & -vel(iel,k))
121  enddo
122  else
123  do k=1,3
124  bv(i,k)=bv(i,k)-vfa(k,ifa)*xflux
125  enddo
126  endif
127  else
128 !
129 ! incoming flux through boundary
130 !
131  if(ielfa(2,ifa).lt.0) then
132  indexb=-ielfa(2,ifa)
133  if(((ifabou(indexb+1).ne.0).and.
134  & (ifabou(indexb+2).ne.0).and.
135  & (ifabou(indexb+3).ne.0)).or.
136  & (dabs(xflux).lt.1.d-10)) then
137  do k=1,3
138  bv(i,k)=bv(i,k)-vfa(k,ifa)*xflux
139  enddo
140  else
141 c write(*,*) '*ERROR in mafillv: not all'
142 c write(*,*) ' components of an incoming'
143 c write(*,*) ' flux through face ',
144 c & indexf-ipnei(i)
145 c write(*,*)' of element ',nactdohinv(i),
146 c & ' are given',
147 c &vfa(1,ifa),vfa(2,ifa),vfa(3,ifa)
148  endif
149  else
150 c write(*,*) '*ERROR in mafillv: not all'
151 c write(*,*) ' components of an incoming'
152 c write(*,*) ' flux through face ',
153 c & indexf-ipnei(i)
154 c write(*,*)' of element ',nactdohinv(i),
155 c & ' are given',
156 c &vfa(1,ifa),vfa(2,ifa),vfa(3,ifa)
157  endif
158  endif
159  endif
160 !
161 ! diffusion
162 !
163  if(iturbulent.eq.0) then
164  difcoef=umfa(ifa)
165  else
166  difcoef=umfa(ifa)+vfa(5,ifa)*vfa(6,ifa)/vfa(7,ifa)
167  endif
168 !
169  if(iel.ne.0) then
170 !
171 ! neighboring element
172 !
173  coef=difcoef*area(ifa)/xlet(indexf)
174  adv(i)=adv(i)+coef
175  if((icyclic.eq.0).or.(ifatie(ifa).eq.0)) then
176  auv(indexf)=auv(indexf)-coef
177  elseif(ifatie(ifa).gt.0) then
178 !
179 ! for cyclic symmetry the term is retarded, since
180 ! otherwise the x-, y- and z- components are linked
181 ! (i.e. the x-, y- and z- momentum equations cannot
182 ! be solved separately any more)
183 !
184  do k=1,3
185  bv(i,k)=bv(i,k)+
186  & (c(k,1)*vel(iel,1)+c(k,2)*vel(iel,2)+
187  & c(k,3)*vel(iel,3))*coef
188  enddo
189  else
190  do k=1,3
191  bv(i,k)=bv(i,k)+
192  & (c(1,k)*vel(iel,1)+c(2,k)*vel(iel,2)+
193  & c(3,k)*vel(iel,3))*coef
194  enddo
195  endif
196 !
197 ! correction for non-orthogonal grid
198 !
199  do k=1,3
200  bv(i,k)=bv(i,k)+difcoef*area(ifa)*
201  & (gradvfa(k,1,ifa)*xxnj(1,indexf)+
202  & gradvfa(k,2,ifa)*xxnj(2,indexf)+
203  & gradvfa(k,3,ifa)*xxnj(3,indexf))
204  enddo
205  else
206 !
207 ! boundary; check whether wall (specified by user),
208 ! outlet (no velocity boundary conditions or
209 ! none of those
210 !
211  iwall=0
212  ipointer=abs(ielfa(2,ifa))
213  if(ipointer.gt.0) then
214  iwall=ifabou(ipointer+5)
215  endif
216 c if(iwall.lt.1) then
217  if(iwall.eq.0) then
218 !
219 ! external face, but no wall
220 !
221 c if((ifabou(ipointer+1).ne.0).or.
222 c & (ifabou(ipointer+2).ne.0).or.
223 c & (ifabou(ipointer+3).ne.0)) then
224 c
225  if(ielfa(3,ifa).gt.0) then
226 !
227 ! no outlet: face velocity fixed
228 !
229  coef=difcoef*area(ifa)/xle(indexf)
230  adv(i)=adv(i)+coef
231  do k=1,3
232  bv(i,k)=bv(i,k)+coef*vfa(k,ifa)
233  enddo
234  else
235 !
236 ! outlet: no diffusion
237 !
238  endif
239 !
240 ! correction for non-orthogonal grid
241 !
242  do k=1,3
243  bv(i,k)=bv(i,k)+difcoef*area(ifa)*
244  & (gradvfa(k,1,ifa)*xxni(1,indexf)+
245  & gradvfa(k,2,ifa)*xxni(2,indexf)+
246  & gradvfa(k,3,ifa)*xxni(3,indexf))
247  enddo
248 c else
249  elseif(iwall.gt.0) then
250 !
251 ! wall
252 !
253  coef=difcoef*area(ifa)/(xle(indexf)*cosa(indexf))
254  adv(i)=adv(i)+coef
255 !
256 ! correction for non-orthogonal grid and nonzero
257 ! wall velocity
258 !
259  coef2=((vel(i,1)-vfa(1,ifa))*xxn(1,indexf)+
260  & (vel(i,2)-vfa(2,ifa))*xxn(2,indexf)+
261  & (vel(i,3)-vfa(3,ifa))*xxn(3,indexf))*coef
262  do k=1,3
263  bv(i,k)=bv(i,k)+coef*vfa(k,ifa)+
264  & coef2*xxn(k,indexf)
265  enddo
266  else
267 !
268 ! sliding conditions (no shear stress)
269 !
270  coef=difcoef*area(ifa)/(xle(indexf)*cosa(indexf))
271 !
272 ! correction for non-orthogonal grid
273 !
274  coef2=((vel(i,1)-vfa(1,ifa))*xxn(1,indexf)+
275  & (vel(i,2)-vfa(2,ifa))*xxn(2,indexf)+
276  & (vel(i,3)-vfa(3,ifa))*xxn(3,indexf))*coef
277  do k=1,3
278  bv(i,k)=bv(i,k)-coef2*xxn(k,indexf)
279  enddo
280  endif
281  endif
282 !
283 ! pressure
284 !
285  enddo
286 !
287 ! body force
288 !
289  rhovel=vel(i,5)*volume(i)
290 !
291  if(nbody.gt.0) then
292  do k=1,3
293  bv(i,k)=bv(i,k)+rhovel*body(k,i)
294  enddo
295  endif
296 !
297 ! transient term
298 !
299  constant=rhovel
300  do k=1,3
301  bv(i,k)=bv(i,k)-(a2*velo(i,k)+a3*veloo(i,k))*constant
302  enddo
303  constant=a1*constant
304  adv(i)=adv(i)+constant
305 !
306 ! copying b into sel (rhs without pressure)
307 !
308  do j=1,3
309  sel(j,i)=bv(i,j)
310  enddo
311 !
312 ! pressure contribution to b
313 !
314  if(iturbulent.eq.0) then
315  do indexf=ipnei(i)+1,ipnei(i+1)
316  ifa=neifa(indexf)
317  do k=1,3
318  bv(i,k)=bv(i,k)
319  & -vfa(4,ifa)*xxn(k,indexf)*area(ifa)
320  enddo
321  enddo
322  else
323  do indexf=ipnei(i)+1,ipnei(i+1)
324  ifa=neifa(indexf)
325  do k=1,3
326  bv(i,k)=bv(i,k)
327  & -(vfa(4,ifa)+2.d0*vfa(5,ifa)*vfa(6,ifa)/3.d0)
328  & *xxn(k,indexf)*area(ifa)
329  enddo
330  enddo
331  endif
332 !
333  enddo
334 !
335  return
flux
subroutine flux(node1, node2, nodem, nelem, lakon, kon, ipkon, nactdog, identity, ielprop, prop, kflag, v, xflow, f, nodef, idirf, df, cp, R, rho, physcon, g, co, dvi, numf, vold, set, shcon, nshcon, rhcon, nrhcon, ntmat_, mi, ider, ttime, time, iaxial)
Definition: flux.f:24
Hosted by OpenAircraft.com, (Michigan UAV, LLC)