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

Go to the source code of this file.

Functions/Subroutines

subroutine mafillpcomp (nef, lakonf, ipnei, neifa, neiel, vfa, area, advfa, xlet, cosa, volume, au, ad, jq, irow, ap, ielfa, ifabou, xle, b, xxn, neq, nzs, hfa, gradpel, bp, xxi, neij, xlen, cosb, ielmatf, mi, a1, a2, a3, velo, veloo, dtimef, shcon, ntmat_, vel, nactdohinv, xrlfa, flux, nefa, nefb, iau6, xxicn, gamma)
 

Function/Subroutine Documentation

◆ mafillpcomp()

subroutine mafillpcomp ( integer, intent(in)  nef,
character*8, dimension(*), intent(in)  lakonf,
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(*), intent(in)  area,
real*8, dimension(*), intent(in)  advfa,
real*8, dimension(*), intent(in)  xlet,
real*8, dimension(*), intent(in)  cosa,
real*8, dimension(*), intent(in)  volume,
real*8, dimension(nzs), intent(inout)  au,
real*8, dimension(neq), intent(inout)  ad,
integer, dimension(*), intent(in)  jq,
integer, dimension(*), intent(in)  irow,
real*8, dimension(*), intent(inout)  ap,
integer, dimension(4,*), intent(in)  ielfa,
integer, dimension(*), intent(in)  ifabou,
real*8, dimension(*), intent(in)  xle,
real*8, dimension(neq), intent(inout)  b,
real*8, dimension(3,*), intent(in)  xxn,
integer  neq,
integer, intent(in)  nzs,
real*8, dimension(3,*), intent(in)  hfa,
real*8, dimension(3,*), intent(in)  gradpel,
real*8, dimension(*), intent(inout)  bp,
real*8, dimension(3,*), intent(in)  xxi,
integer, dimension(*), intent(in)  neij,
real*8, dimension(*), intent(in)  xlen,
real*8, dimension(*), intent(in)  cosb,
integer, dimension(mi(3),*), intent(in)  ielmatf,
integer, dimension(*), intent(in)  mi,
real*8, intent(in)  a1,
real*8, intent(in)  a2,
real*8, intent(in)  a3,
real*8, dimension(nef,0:7), intent(in)  velo,
real*8, dimension(nef,0:7), intent(in)  veloo,
real*8, intent(in)  dtimef,
real*8, dimension(0:3,ntmat_,*), intent(in)  shcon,
integer, intent(in)  ntmat_,
real*8, dimension(nef,0:7), intent(in)  vel,
integer, dimension(*), intent(in)  nactdohinv,
real*8, dimension(3,*), intent(in)  xrlfa,
real*8, dimension(*), intent(in)  flux,
integer  nefa,
integer  nefb,
integer, dimension(6,*)  iau6,
real*8, dimension(3,*)  xxicn,
real*8, dimension(*), intent(in)  gamma 
)
25 !
26 ! filling the lhs and rhs to calculate p
27 !
28  implicit none
29 !
30  character*8 lakonf(*)
31 !
32  integer i,nef,indexf,ipnei(*),j,neifa(*),iel3,k,
33  & neiel(*),iel,ifa,irow(*),ielfa(4,*),compressible,
34  & ifabou(*),neq,jq(*),iel2,indexb,knownflux,indexf2,
35  & j2,neij(*),nzs,imat,nefa,nefb,iau6(6,*),
36  & mi(*),ielmatf(mi(3),*),ntmat_,nactdohinv(*),knownpressure
37 !
38  real*8 coef,vfa(0:7,*),volume(*),area(*),advfa(*),xlet(*),
39  & cosa(*),ad(neq),au(nzs),xle(*),xxn(3,*),ap(*),b(neq),cosb(*),
40  & hfa(3,*),gradpel(3,*),bp(*),xxi(3,*),xlen(*),r,a1,a2,a3,
41  & xflux,constant,velo(nef,0:7),veloo(nef,0:7),dtimef,
42  & shcon(0:3,ntmat_,*),vel(nef,0:7),dd,convec,fluxisobar,
43  & coef1,coef3,xrlfa(3,*),coef2,gamma(*),coefp,coefn,xmach,
44  & flux(*),bp_ifa,aa(8,8),xxicn(3,*)
45 !
46  intent(in) nef,lakonf,ipnei,neifa,neiel,vfa,area,
47  & advfa,xlet,cosa,volume,jq,irow,ielfa,ifabou,xle,
48  & xxn,nzs,hfa,gradpel,xxi,neij,
49  & xlen,cosb,ielmatf,mi,a1,a2,a3,velo,veloo,dtimef,shcon,
50  & ntmat_,vel,nactdohinv,xrlfa,flux,gamma
51 !
52  intent(inout) ad,au,b,ap,bp
53 !
54  do i=nefa,nefb
55  imat=ielmatf(1,i)
56  r=shcon(3,1,imat)
57  indexf=ipnei(i)
58  do j=1,ipnei(i+1)-ipnei(i)
59  knownflux=0
60  knownpressure=0
61  convec=0
62 !
63 ! diffusion
64 !
65  indexf=indexf+1
66  ifa=neifa(indexf)
67  iel=neiel(indexf)
68  if(iel.ne.0) then
69  coef=vfa(5,ifa)*(volume(i)+volume(iel))*area(ifa)/
70  & (advfa(ifa)*2.d0*xlet(indexf)*cosb(indexf))
71  ad(i)=ad(i)+coef
72  au(indexf)=au(indexf)-coef
73  b(i)=b(i)+coef*(vel(iel,4)-vel(i,4))
74  convec=coef*(vel(iel,4)-vel(i,4))
75 !
76 ! correction for non-orthogonal meshes
77 !
78  j2=neij(indexf)
79  indexf2=ipnei(iel)+j2
80  bp_ifa=((gradpel(1,iel)*xxicn(1,indexf2)+
81  & gradpel(2,iel)*xxicn(2,indexf2)+
82  & gradpel(3,iel)*xxicn(3,indexf2))
83  & *xle(indexf2)
84  & -(gradpel(1,i)*xxicn(1,indexf)+
85  & gradpel(2,i)*xxicn(2,indexf)+
86  & gradpel(3,i)*xxicn(3,indexf))
87  & *xle(indexf))
88  b(i)=b(i)+coef*bp_ifa
89  convec=convec+coef*bp_ifa
90 !
91 ! following line is correct if the temperature
92 ! changes from one pressure correction iteration to
93 ! the next
94 !
95  convec=-convec/(vfa(5,ifa)*r*vfa(0,ifa))
96  else
97  iel2=ielfa(2,ifa)
98  if(iel2.lt.0) then
99  if((ifabou(-iel2+1).ne.0).and.
100  & (ifabou(-iel2+2).ne.0).and.
101  & (ifabou(-iel2+3).ne.0)) then
102 !
103 ! all velocity components given
104 !
105  knownflux=1
106  elseif(ifabou(-iel2+5).lt.0) then
107 !
108 ! sliding conditions
109 !
110  knownflux=2
111  elseif(ifabou(-iel2+4).ne.0) then
112  knownpressure=1
113 !
114 ! pressure given (only if not all velocity
115 ! components are given)
116 !
117  coef=vfa(5,ifa)*volume(i)*area(ifa)/
118  & (advfa(ifa)*xle(indexf)*cosa(indexf))
119  ad(i)=ad(i)+coef
120  b(i)=b(i)+coef*vfa(4,ifa)
121  & -coef*vel(i,4)
122  convec=coef*vfa(4,ifa)-coef*vel(i,4)
123 !
124 ! correction for non-orthogonal meshes
125 !
126  bp_ifa=(-(gradpel(1,i)*xxicn(1,indexf)+
127  & gradpel(2,i)*xxicn(2,indexf)+
128  & gradpel(3,i)*xxicn(3,indexf))
129  & *xle(indexf))
130  b(i)=b(i)+coef*bp_ifa
131  convec=convec+coef*bp_ifa
132 !
133 ! following line is correct if the temperature
134 ! changes from one pressure correction iteration to
135 ! the next
136 !
137  convec=-convec/(vfa(5,ifa)*r*vfa(0,ifa))
138  endif
139  endif
140  endif
141 !
142 ! save coefficients for correctvfa.f
143 !
144  if((iel.eq.0).or.(i.lt.iel)) then
145  ap(ifa)=coef
146  bp(ifa)=bp_ifa
147  endif
148 !
149 ! convection
150 !
151 ! flux
152 !
153  if(knownflux.eq.1) then
154  xflux=flux(indexf)
155  elseif(knownflux.eq.2) then
156  xflux=0.d0
157  endif
158 !
159 ! flux based on constant pressure
160 !
161  if(knownflux.eq.0) then
162  fluxisobar=area(ifa)*
163  & (hfa(1,ifa)*xxn(1,indexf)+
164  & hfa(2,ifa)*xxn(2,indexf)+
165  & hfa(3,ifa)*xxn(3,indexf))
166  endif
167 !
168 ! rhs
169 !
170  if(knownflux.eq.0) then
171  b(i)=b(i)-vfa(5,ifa)*fluxisobar
172  elseif(knownflux.eq.1) then
173  b(i)=b(i)-xflux
174  endif
175 !
176  if(knownflux.eq.0) then
177 !
178 ! following line leads to oscillations in the solution
179 ! (only for subsonic and transonic solutions)
180 !
181  coef=fluxisobar/(r*vfa(0,ifa))+convec
182  elseif(knownflux.eq.1) then
183  coef=xflux/(r*vfa(0,ifa)*vfa(5,ifa))
184  else
185  coef=0.d0
186  endif
187 !
188  if(coef.ge.0.d0) then
189 !
190 ! outflowing flux
191 !
192  ad(i)=ad(i)+coef
193 c
194 c retarded central difference
195 c b(i)=b(i)-gamma(ifa)*(vfa(4,ifa)-vel(i,4))*coef
196 c
197  else
198  if(iel.gt.0) then
199 !
200 ! incoming flux from neighboring element
201 !
202  au(indexf)=au(indexf)+coef
203 !
204 c
205 c retarded central difference
206 c b(i)=b(i)-gamma(ifa)*(vfa(4,ifa)-vel(iel,4))*coef
207 c
208  elseif(knownpressure.eq.0) then
209  ad(i)=ad(i)+coef
210  endif
211  endif
212 !
213  enddo
214 !
215 ! transient term
216 !
217 c a1=1.d0/dtimef
218 c a2=-1.d0/dtimef
219 c a3=0.d0/dtimef
220 c constant=volume(i)/(r*vel(i,0))
221 c b(i)=b(i)-
222 c & (a1*vel(i,5)+a2*velo(i,5)+a3*veloo(i,5))*volume(i)
223  b(i)=b(i)-
224  & (vel(i,5)-velo(i,5))*volume(i)/dtimef
225 c constant=a1*constant
226  constant=volume(i)/(r*vel(i,0)*dtimef)
227  ad(i)=ad(i)+constant
228 !
229  enddo
230 !
231  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)