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

Go to the source code of this file.

Functions/Subroutines

subroutine mafillvcomp (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)
 

Function/Subroutine Documentation

◆ mafillvcomp()

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