2 #ifndef DUNE_PDELAB_DIFFUSIONMIXED_HH
3 #define DUNE_PDELAB_DIFFUSIONMIXED_HH
8 #include<dune/common/exceptions.hh>
9 #include<dune/common/fvector.hh>
10 #include<dune/common/fmatrix.hh>
12 #include<dune/geometry/type.hh>
13 #include<dune/geometry/quadraturerules.hh>
14 #include<dune/geometry/referenceelements.hh>
35 template<
typename PARAM>
63 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
64 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
67 typedef typename LFSU::template Child<0>::Type VelocitySpace;
68 const VelocitySpace& velocityspace = lfsu.template child<0>();
69 typedef typename LFSU::template Child<1>::Type PressureSpace;
70 const PressureSpace& pressurespace = lfsu.template child<1>();
73 typedef typename VelocitySpace::Traits::FiniteElementType::
74 Traits::LocalBasisType::Traits::DomainFieldType DF;
75 typedef typename VelocitySpace::Traits::FiniteElementType::
76 Traits::LocalBasisType::Traits::RangeFieldType RF;
77 typedef typename VelocitySpace::Traits::FiniteElementType::
78 Traits::LocalBasisType::Traits::JacobianType VelocityJacobianType;
79 typedef typename VelocitySpace::Traits::FiniteElementType::
80 Traits::LocalBasisType::Traits::RangeType VelocityRangeType;
81 typedef typename PressureSpace::Traits::FiniteElementType::
82 Traits::LocalBasisType::Traits::RangeType PressureRangeType;
85 const int dim = EG::Geometry::dimension;
88 Dune::FieldVector<DF,dim> pos; pos = 0.0;
89 typename EG::Geometry::JacobianInverseTransposed jac;
90 jac = eg.geometry().jacobianInverseTransposed(pos);
92 RF det = eg.geometry().integrationElement(pos);
96 Dune::GeometryType gt = eg.geometry().type();
97 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
99 typename PARAM::Traits::PermTensorType tensor;
100 tensor = param.A(eg.entity(),localcenter);
104 const Dune::QuadratureRule<DF,dim>& vrule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder_v);
106 for (
const auto& ip : vrule)
109 std::vector<VelocityRangeType> vbasis(velocityspace.size());
110 velocityspace.finiteElement().localBasis().evaluateFunction(ip.position(),vbasis);
113 std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size());
114 for (std::size_t i=0; i<velocityspace.size(); i++)
116 vtransformedbasis[i] = 0.0;
117 jac.umtv(vbasis[i],vtransformedbasis[i]);
121 VelocityRangeType sigma; sigma = 0.0;
122 for (std::size_t i=0; i<velocityspace.size(); i++)
123 sigma.axpy(x(velocityspace,i),vtransformedbasis[i]);
126 VelocityRangeType Kinvsigma;
127 tensor.mv(sigma,Kinvsigma);
130 RF factor = ip.weight() / det;
131 for (std::size_t i=0; i<velocityspace.size(); i++)
132 r.accumulate(velocityspace,i,(Kinvsigma*vtransformedbasis[i])*factor);
136 const Dune::QuadratureRule<DF,dim>& prule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder_p);
138 for (
const auto& ip : prule)
141 std::vector<VelocityJacobianType> vbasis(velocityspace.size());
142 velocityspace.finiteElement().localBasis().evaluateJacobian(ip.position(),vbasis);
143 std::vector<PressureRangeType> pbasis(pressurespace.size());
144 pressurespace.finiteElement().localBasis().evaluateFunction(ip.position(),pbasis);
147 PressureRangeType u; u = 0.0;
148 for (std::size_t i=0; i<pressurespace.size(); i++)
149 u.axpy(x(pressurespace,i),pbasis[i]);
152 typename PARAM::Traits::RangeFieldType a0value = param.c(eg.entity(),ip.position());
155 RF factor = ip.weight();
156 for (std::size_t i=0; i<pressurespace.size(); i++)
157 r.accumulate(pressurespace,i,-a0value*u*pbasis[i]*factor);
160 std::vector<RF> divergence(velocityspace.size(),0.0);
161 for (std::size_t i=0; i<velocityspace.size(); i++)
162 for (
int j=0; j<
dim; j++)
163 divergence[i] += vbasis[i][j][j];
166 for (std::size_t i=0; i<velocityspace.size(); i++)
167 r.accumulate(velocityspace,i,-u*divergence[i]*factor);
170 RF divergencesigma = 0.0;
171 for (std::size_t i=0; i<velocityspace.size(); i++)
172 divergencesigma += x(velocityspace,i)*divergence[i];
175 for (std::size_t i=0; i<pressurespace.size(); i++)
176 r.accumulate(pressurespace,i,-divergencesigma*pbasis[i]*factor);
181 template<
typename EG,
typename LFSV,
typename R>
185 typedef typename LFSV::template Child<1>::Type PressureSpace;
186 const PressureSpace& pressurespace = lfsv.template child<1>();
189 typedef typename PressureSpace::Traits::FiniteElementType::
190 Traits::LocalBasisType::Traits::DomainFieldType DF;
191 typedef typename PressureSpace::Traits::FiniteElementType::
192 Traits::LocalBasisType::Traits::RangeFieldType RF;
193 typedef typename PressureSpace::Traits::FiniteElementType::
194 Traits::LocalBasisType::Traits::RangeType PressureRangeType;
197 const int dim = EG::Geometry::dimension;
200 Dune::GeometryType gt = eg.geometry().type();
201 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder_p);
203 for (
const auto& ip : rule)
206 std::vector<PressureRangeType> pbasis(pressurespace.size());
207 pressurespace.finiteElement().localBasis().evaluateFunction(ip.position(),pbasis);
210 RF y = param.f(eg.entity(),ip.position());
213 RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
214 for (std::size_t i=0; i<pressurespace.size(); i++)
215 r.accumulate(pressurespace,i,y*pbasis[i]*factor);
220 template<
typename IG,
typename LFSV,
typename R>
224 typedef typename LFSV::template Child<0>::Type VelocitySpace;
225 const VelocitySpace& velocityspace = lfsv.template child<0>();
228 typedef typename VelocitySpace::Traits::FiniteElementType::
229 Traits::LocalBasisType::Traits::DomainFieldType DF;
230 typedef typename VelocitySpace::Traits::FiniteElementType::
231 Traits::LocalBasisType::Traits::RangeFieldType RF;
232 typedef typename VelocitySpace::Traits::FiniteElementType::
233 Traits::LocalBasisType::Traits::RangeType VelocityRangeType;
236 const int dim = IG::dimension;
238 auto inside_cell = ig.inside();
241 Dune::FieldVector<DF,dim> pos; pos = 0.0;
242 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
243 jac = inside_cell.geometry().jacobianInverseTransposed(pos);
245 RF det = inside_cell.geometry().integrationElement(pos);
248 Dune::GeometryType gtface = ig.geometryInInside().type();
249 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder_v);
252 for (
const auto& ip : rule)
255 BCType bctype = param.bctype(ig.intersection(),ip.position());
262 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
265 std::vector<VelocityRangeType> vbasis(velocityspace.size());
266 velocityspace.finiteElement().localBasis().evaluateFunction(local,vbasis);
269 std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size());
270 for (std::size_t i=0; i<velocityspace.size(); i++)
272 vtransformedbasis[i] = 0.0;
273 jac.umtv(vbasis[i],vtransformedbasis[i]);
277 RF y = param.g(inside_cell,local);
280 RF factor = ip.weight()*ig.geometry().integrationElement(ip.position())/det;
281 for (std::size_t i=0; i<velocityspace.size(); i++)
282 r.accumulate(velocityspace,i,y*(vtransformedbasis[i]*ig.unitOuterNormal(ip.position()))*factor);
Definition: diffusionmixed.hh:49
Definition: diffusionmixed.hh:50
Definition: diffusionmixed.hh:46
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
const IG & ig
Definition: constraints.hh:147
Type
Definition: convectiondiffusionparameter.hh:67
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
Definition: diffusionmixed.hh:51
Definition: convectiondiffusionparameter.hh:67
Definition: adaptivity.hh:27
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
static const int dim
Definition: adaptivity.hh:83
DiffusionMixed(const PARAM ¶m_, int qorder_v_=2, int qorder_p_=1)
Definition: diffusionmixed.hh:53
Definition: diffusionmixed.hh:36
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:182
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:64
const EG & eg
Definition: constraints.hh:280
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:221