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>
11 #include<dune/common/static_assert.hh>
13 #include<dune/geometry/type.hh>
14 #include<dune/geometry/quadraturerules.hh>
15 #include<dune/geometry/referenceelements.hh>
36 template<
typename K,
typename A0,
typename F,
typename B,
typename G>
51 DiffusionMixed (
const K& k_,
const A0& a0_,
const F& f_,
const B& bctype_,
const G& g_,
int qorder_v_=2,
int qorder_p_=1)
52 : k(k_), a0(a0_), f(f_), bctype(bctype_), g(g_), qorder_v(qorder_v_), qorder_p(qorder_p_)
56 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
57 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
60 typedef typename LFSU::template Child<0>::Type VelocitySpace;
61 const VelocitySpace& velocityspace = lfsu.template child<0>();
62 typedef typename LFSU::template Child<1>::Type PressureSpace;
63 const PressureSpace& pressurespace = lfsu.template child<1>();
66 typedef typename VelocitySpace::Traits::FiniteElementType::
67 Traits::LocalBasisType::Traits::DomainFieldType DF;
68 typedef typename VelocitySpace::Traits::FiniteElementType::
69 Traits::LocalBasisType::Traits::RangeFieldType RF;
70 typedef typename VelocitySpace::Traits::FiniteElementType::
71 Traits::LocalBasisType::Traits::JacobianType VelocityJacobianType;
72 typedef typename VelocitySpace::Traits::FiniteElementType::
73 Traits::LocalBasisType::Traits::RangeType VelocityRangeType;
74 typedef typename PressureSpace::Traits::FiniteElementType::
75 Traits::LocalBasisType::Traits::RangeType PressureRangeType;
78 const int dim = EG::Geometry::dimension;
81 Dune::FieldVector<DF,dim> pos; pos = 0.0;
82 typename EG::Geometry::JacobianInverseTransposed jac;
83 jac = eg.geometry().jacobianInverseTransposed(pos);
85 RF det = eg.geometry().integrationElement(pos);
88 typename K::Traits::RangeType tensor;
89 Dune::GeometryType gt = eg.geometry().type();
90 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
91 k.evaluate(eg.entity(),localcenter,tensor);
95 const Dune::QuadratureRule<DF,dim>& vrule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder_v);
96 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=vrule.begin(); it!=vrule.end(); ++it)
99 std::vector<VelocityRangeType> vbasis(velocityspace.size());
100 velocityspace.finiteElement().localBasis().evaluateFunction(it->position(),vbasis);
103 std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size());
104 for (std::size_t i=0; i<velocityspace.size(); i++)
106 vtransformedbasis[i] = 0.0;
107 jac.umtv(vbasis[i],vtransformedbasis[i]);
111 VelocityRangeType sigma; sigma = 0.0;
112 for (std::size_t i=0; i<velocityspace.size(); i++)
113 sigma.axpy(x(velocityspace,i),vtransformedbasis[i]);
116 VelocityRangeType Kinvsigma;
117 tensor.mv(sigma,Kinvsigma);
120 RF factor = it->weight() / det;
121 for (std::size_t i=0; i<velocityspace.size(); i++)
122 r.accumulate(velocityspace,i,(Kinvsigma*vtransformedbasis[i])*factor);
126 const Dune::QuadratureRule<DF,dim>& prule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder_p);
127 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=prule.begin(); it!=prule.end(); ++it)
130 std::vector<VelocityJacobianType> vbasis(velocityspace.size());
131 velocityspace.finiteElement().localBasis().evaluateJacobian(it->position(),vbasis);
132 std::vector<PressureRangeType> pbasis(pressurespace.size());
133 pressurespace.finiteElement().localBasis().evaluateFunction(it->position(),pbasis);
136 PressureRangeType u; u = 0.0;
137 for (std::size_t i=0; i<pressurespace.size(); i++)
138 u.axpy(x(pressurespace,i),pbasis[i]);
141 typename A0::Traits::RangeType a0value;
142 a0.evaluate(eg.entity(),it->position(),a0value);
145 RF factor = it->weight();
146 for (std::size_t i=0; i<pressurespace.size(); i++)
147 r.accumulate(pressurespace,i,-a0value*u*pbasis[i]*factor);
150 std::vector<RF> divergence(velocityspace.size(),0.0);
151 for (std::size_t i=0; i<velocityspace.size(); i++)
152 for (
int j=0; j<
dim; j++)
153 divergence[i] += vbasis[i][j][j];
156 for (std::size_t i=0; i<velocityspace.size(); i++)
157 r.accumulate(velocityspace,i,-u*divergence[i]*factor);
160 RF divergencesigma = 0.0;
161 for (std::size_t i=0; i<velocityspace.size(); i++)
162 divergencesigma += x(velocityspace,i)*divergence[i];
165 for (std::size_t i=0; i<pressurespace.size(); i++)
166 r.accumulate(pressurespace,i,-divergencesigma*pbasis[i]*factor);
171 template<
typename EG,
typename LFSV,
typename R>
175 typedef typename LFSV::template Child<1>::Type PressureSpace;
176 const PressureSpace& pressurespace = lfsv.template child<1>();
179 typedef typename PressureSpace::Traits::FiniteElementType::
180 Traits::LocalBasisType::Traits::DomainFieldType DF;
181 typedef typename PressureSpace::Traits::FiniteElementType::
182 Traits::LocalBasisType::Traits::RangeFieldType RF;
183 typedef typename PressureSpace::Traits::FiniteElementType::
184 Traits::LocalBasisType::Traits::RangeType PressureRangeType;
187 const int dim = EG::Geometry::dimension;
190 Dune::GeometryType gt = eg.geometry().type();
191 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder_p);
194 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
197 std::vector<PressureRangeType> pbasis(pressurespace.size());
198 pressurespace.finiteElement().localBasis().evaluateFunction(it->position(),pbasis);
201 typename F::Traits::RangeType y;
202 f.evaluate(eg.entity(),it->position(),y);
205 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
206 for (std::size_t i=0; i<pressurespace.size(); i++)
207 r.accumulate(pressurespace,i,y*pbasis[i]*factor);
212 template<
typename IG,
typename LFSV,
typename R>
216 typedef typename LFSV::template Child<0>::Type VelocitySpace;
217 const VelocitySpace& velocityspace = lfsv.template child<0>();
220 typedef typename VelocitySpace::Traits::FiniteElementType::
221 Traits::LocalBasisType::Traits::DomainFieldType DF;
222 typedef typename VelocitySpace::Traits::FiniteElementType::
223 Traits::LocalBasisType::Traits::RangeFieldType RF;
224 typedef typename VelocitySpace::Traits::FiniteElementType::
225 Traits::LocalBasisType::Traits::RangeType VelocityRangeType;
228 const int dim = IG::dimension;
231 Dune::FieldVector<DF,dim> pos; pos = 0.0;
232 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
233 jac = ig.inside()->geometry().jacobianInverseTransposed(pos);
235 RF det = ig.inside()->geometry().integrationElement(pos);
238 Dune::GeometryType gtface = ig.geometryInInside().type();
239 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder_v);
242 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
246 if( bctype.isNeumann( ig,it->position() ) )
250 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
253 std::vector<VelocityRangeType> vbasis(velocityspace.size());
254 velocityspace.finiteElement().localBasis().evaluateFunction(local,vbasis);
257 std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size());
258 for (std::size_t i=0; i<velocityspace.size(); i++)
260 vtransformedbasis[i] = 0.0;
261 jac.umtv(vbasis[i],vtransformedbasis[i]);
265 typename G::Traits::RangeType y;
266 g.evaluate(*(ig.inside()),local,y);
269 RF factor = it->weight()*ig.geometry().integrationElement(it->position())/det;
270 for (std::size_t i=0; i<velocityspace.size(); i++)
271 r.accumulate(velocityspace,i,y*(vtransformedbasis[i]*ig.unitOuterNormal(it->position()))*factor);
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:213
Definition: diffusionmixed.hh:37
Default flags for all local operators.
Definition: flags.hh:18
Definition: diffusionmixed.hh:49
static const int dim
Definition: adaptivity.hh:82
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:172
const IG & ig
Definition: common/constraints.hh:146
DiffusionMixed(const K &k_, const A0 &a0_, const F &f_, const B &bctype_, const G &g_, int qorder_v_=2, int qorder_p_=1)
Definition: diffusionmixed.hh:51
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:57
sparsity pattern generator
Definition: pattern.hh:14
Definition: diffusionmixed.hh:47
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
const EG & eg
Definition: common/constraints.hh:277
Definition: diffusionmixed.hh:44
Definition: diffusionmixed.hh:48