2 #ifndef DUNE_PDELAB_LOCALOPERATOR_DIFFUSIONMIXED_HH 3 #define DUNE_PDELAB_LOCALOPERATOR_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> 38 template<
typename PARAM>
66 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
67 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 70 using VelocitySpace =
typename LFSU::template Child<0>::Type;
71 using PressureSpace =
typename LFSU::template Child<1>::Type;
72 using DF =
typename VelocitySpace::Traits::FiniteElementType::
73 Traits::LocalBasisType::Traits::DomainFieldType;
74 using RF =
typename VelocitySpace::Traits::FiniteElementType::
75 Traits::LocalBasisType::Traits::RangeFieldType;
76 using VelocityJacobianType =
typename VelocitySpace::Traits::FiniteElementType::
77 Traits::LocalBasisType::Traits::JacobianType;
78 using VelocityRangeType =
typename VelocitySpace::Traits::FiniteElementType::
79 Traits::LocalBasisType::Traits::RangeType;
80 using PressureRangeType =
typename PressureSpace::Traits::FiniteElementType::
81 Traits::LocalBasisType::Traits::RangeType;
84 using namespace Indices;
85 const auto& velocityspace = child(lfsu,_0);
86 const auto& pressurespace = lfsu.template child<1>();
89 const int dim = EG::Geometry::mydimension;
92 const auto& cell = eg.entity();
95 auto geo = eg.geometry();
98 Dune::FieldVector<DF,dim> pos;
100 auto jac = geo.jacobianInverseTransposed(pos);
102 auto det = geo.integrationElement(pos);
106 auto localcenter = ref_el.position(0,0);
107 auto tensor = param.A(cell,localcenter);
111 std::vector<VelocityRangeType> vbasis(velocityspace.size());
112 std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size());
113 VelocityRangeType sigma;
114 VelocityRangeType Kinvsigma;
115 std::vector<VelocityJacobianType> vjacbasis(velocityspace.size());
116 std::vector<PressureRangeType> pbasis(pressurespace.size());
117 std::vector<RF> divergence(velocityspace.size(),0.0);
125 velocityspace.finiteElement().localBasis().evaluateFunction(ip.position(),vbasis);
128 for (std::size_t i=0; i<velocityspace.size(); i++)
130 vtransformedbasis[i] = 0.0;
131 jac.umtv(vbasis[i],vtransformedbasis[i]);
136 for (std::size_t i=0; i<velocityspace.size(); i++)
137 sigma.axpy(x(velocityspace,i),vtransformedbasis[i]);
140 tensor.mv(sigma,Kinvsigma);
143 auto factor = ip.weight() / det;
144 for (std::size_t i=0; i<velocityspace.size(); i++)
145 r.accumulate(velocityspace,i,(Kinvsigma*vtransformedbasis[i])*factor);
153 velocityspace.finiteElement().localBasis().evaluateJacobian(ip.position(),vjacbasis);
154 pressurespace.finiteElement().localBasis().evaluateFunction(ip.position(),pbasis);
159 for (std::size_t i=0; i<pressurespace.size(); i++)
160 u.axpy(x(pressurespace,i),pbasis[i]);
163 auto a0value = param.c(cell,ip.position());
166 RF factor = ip.weight();
167 for (std::size_t i=0; i<pressurespace.size(); i++)
168 r.accumulate(pressurespace,i,-a0value*u*pbasis[i]*factor);
171 for (std::size_t i=0; i<velocityspace.size(); i++){
173 for (
int j=0; j<
dim; j++)
174 divergence[i] += vjacbasis[i][j][j];
178 for (std::size_t i=0; i<velocityspace.size(); i++)
179 r.accumulate(velocityspace,i,-u*divergence[i]*factor);
182 RF divergencesigma = 0.0;
183 for (std::size_t i=0; i<velocityspace.size(); i++)
184 divergencesigma += x(velocityspace,i)*divergence[i];
187 for (std::size_t i=0; i<pressurespace.size(); i++)
188 r.accumulate(pressurespace,i,-divergencesigma*pbasis[i]*factor);
193 template<
typename EG,
typename LFSV,
typename R>
197 using PressureSpace =
typename LFSV::template Child<1>::Type;
198 using PressureRangeType =
typename PressureSpace::Traits::FiniteElementType::
199 Traits::LocalBasisType::Traits::RangeType;
202 using namespace Indices;
203 const auto& pressurespace = child(lfsv,_1);
206 const auto& cell = eg.entity();
209 auto geo = eg.geometry();
212 std::vector<PressureRangeType> pbasis(pressurespace.size());
218 pressurespace.finiteElement().localBasis().evaluateFunction(ip.position(),pbasis);
221 auto y = param.f(cell,ip.position());
224 auto factor = ip.weight() * geo.integrationElement(ip.position());
225 for (std::size_t i=0; i<pressurespace.size(); i++)
226 r.accumulate(pressurespace,i,y*pbasis[i]*factor);
231 template<
typename IG,
typename LFSV,
typename R>
235 using VelocitySpace =
typename LFSV::template Child<0>::Type;
236 using DF =
typename VelocitySpace::Traits::FiniteElementType::
237 Traits::LocalBasisType::Traits::DomainFieldType;
238 using VelocityRangeType =
typename VelocitySpace::Traits::FiniteElementType::
239 Traits::LocalBasisType::Traits::RangeType;
242 using namespace Indices;
243 const auto& velocityspace = child(lfsv,_0);
246 const int dim = IG::dimension;
249 const auto& cell_inside = ig.inside();
252 auto geo = ig.geometry();
253 auto geo_inside = cell_inside.geometry();
256 auto geo_in_inside = ig.geometryInInside();
259 Dune::FieldVector<DF,dim> pos;
261 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
262 jac = geo_inside.jacobianInverseTransposed(pos);
264 auto det = geo_inside.integrationElement(pos);
267 std::vector<VelocityRangeType> vbasis(velocityspace.size());
268 std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size());
274 auto bctype = param.bctype(ig.intersection(),ip.position());
281 auto local = geo_in_inside.global(ip.position());
284 velocityspace.finiteElement().localBasis().evaluateFunction(local,vbasis);
287 for (std::size_t i=0; i<velocityspace.size(); i++)
289 vtransformedbasis[i] = 0.0;
290 jac.umtv(vbasis[i],vtransformedbasis[i]);
294 auto y = param.g(cell_inside,local);
297 auto factor = ip.weight()*geo.integrationElement(ip.position())/det;
298 for (std::size_t i=0; i<velocityspace.size(); i++)
299 r.accumulate(velocityspace,i,y*(vtransformedbasis[i]*ig.unitOuterNormal(ip.position()))*factor);
313 #endif // DUNE_PDELAB_LOCALOPERATOR_DIFFUSIONMIXED_HH const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
Definition: diffusionmixed.hh:53
DiffusionMixed(const PARAM ¶m_, int qorder_v_=2, int qorder_p_=1)
Definition: diffusionmixed.hh:56
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:232
Definition: diffusionmixed.hh:49
Definition: diffusionmixed.hh:52
Definition: diffusionmixed.hh:54
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
Type
Definition: convectiondiffusionparameter.hh:65
Definition: convectiondiffusionparameter.hh:65
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:67
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:194
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume() ...
Definition: numericaljacobianapply.hh:32
Definition: diffusionmixed.hh:39
sparsity pattern generator
Definition: pattern.hh:13
QuadratureRuleWrapper< QuadratureRule< typename Geometry::ctype, Geometry::mydimension > > quadratureRule(const Geometry &geo, std::size_t order, QuadratureType::Enum quadrature_type=QuadratureType::GaussLegendre)
Returns a quadrature rule for the given geometry.
Definition: quadraturerules.hh:117
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
Default flags for all local operators.
Definition: flags.hh:18
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:31