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/referenceelements.hh>
36 template<
typename PARAM>
64 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
65 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
68 using VelocitySpace =
typename LFSU::template Child<0>::Type;
69 using PressureSpace =
typename LFSU::template Child<1>::Type;
70 using DF =
typename VelocitySpace::Traits::FiniteElementType::
71 Traits::LocalBasisType::Traits::DomainFieldType;
72 using RF =
typename VelocitySpace::Traits::FiniteElementType::
73 Traits::LocalBasisType::Traits::RangeFieldType;
74 using VelocityJacobianType =
typename VelocitySpace::Traits::FiniteElementType::
75 Traits::LocalBasisType::Traits::JacobianType;
76 using VelocityRangeType =
typename VelocitySpace::Traits::FiniteElementType::
77 Traits::LocalBasisType::Traits::RangeType;
78 using PressureRangeType =
typename PressureSpace::Traits::FiniteElementType::
79 Traits::LocalBasisType::Traits::RangeType;
82 using namespace Indices;
83 const auto& velocityspace = child(lfsu,_0);
84 const auto& pressurespace = lfsu.template child<1>();
87 const int dim = EG::Geometry::mydimension;
90 const auto& cell = eg.entity();
93 auto geo = eg.geometry();
96 Dune::FieldVector<DF,dim> pos;
98 auto jac = geo.jacobianInverseTransposed(pos);
100 auto det = geo.integrationElement(pos);
103 auto ref_el = referenceElement(geo);
104 auto localcenter = ref_el.position(0,0);
105 auto tensor = param.A(cell,localcenter);
109 std::vector<VelocityRangeType> vbasis(velocityspace.size());
110 std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size());
111 VelocityRangeType sigma;
112 VelocityRangeType Kinvsigma;
113 std::vector<VelocityJacobianType> vjacbasis(velocityspace.size());
114 std::vector<PressureRangeType> pbasis(pressurespace.size());
115 std::vector<RF> divergence(velocityspace.size(),0.0);
123 velocityspace.finiteElement().localBasis().evaluateFunction(ip.position(),vbasis);
126 for (std::size_t i=0; i<velocityspace.size(); i++)
128 vtransformedbasis[i] = 0.0;
129 jac.umtv(vbasis[i],vtransformedbasis[i]);
134 for (std::size_t i=0; i<velocityspace.size(); i++)
135 sigma.axpy(x(velocityspace,i),vtransformedbasis[i]);
138 tensor.mv(sigma,Kinvsigma);
141 auto factor = ip.weight() / det;
142 for (std::size_t i=0; i<velocityspace.size(); i++)
143 r.accumulate(velocityspace,i,(Kinvsigma*vtransformedbasis[i])*factor);
151 velocityspace.finiteElement().localBasis().evaluateJacobian(ip.position(),vjacbasis);
152 pressurespace.finiteElement().localBasis().evaluateFunction(ip.position(),pbasis);
157 for (std::size_t i=0; i<pressurespace.size(); i++)
158 u.axpy(x(pressurespace,i),pbasis[i]);
161 auto a0value = param.c(cell,ip.position());
164 RF factor = ip.weight();
165 for (std::size_t i=0; i<pressurespace.size(); i++)
166 r.accumulate(pressurespace,i,-a0value*u*pbasis[i]*factor);
169 for (std::size_t i=0; i<velocityspace.size(); i++){
171 for (
int j=0; j<
dim; j++)
172 divergence[i] += vjacbasis[i][j][j];
176 for (std::size_t i=0; i<velocityspace.size(); i++)
177 r.accumulate(velocityspace,i,-u*divergence[i]*factor);
180 RF divergencesigma = 0.0;
181 for (std::size_t i=0; i<velocityspace.size(); i++)
182 divergencesigma += x(velocityspace,i)*divergence[i];
185 for (std::size_t i=0; i<pressurespace.size(); i++)
186 r.accumulate(pressurespace,i,-divergencesigma*pbasis[i]*factor);
191 template<
typename EG,
typename LFSV,
typename R>
195 using PressureSpace =
typename LFSV::template Child<1>::Type;
196 using PressureRangeType =
typename PressureSpace::Traits::FiniteElementType::
197 Traits::LocalBasisType::Traits::RangeType;
200 using namespace Indices;
201 const auto& pressurespace = child(lfsv,_1);
204 const auto& cell = eg.entity();
207 auto geo = eg.geometry();
210 std::vector<PressureRangeType> pbasis(pressurespace.size());
216 pressurespace.finiteElement().localBasis().evaluateFunction(ip.position(),pbasis);
219 auto y = param.f(cell,ip.position());
222 auto factor = ip.weight() * geo.integrationElement(ip.position());
223 for (std::size_t i=0; i<pressurespace.size(); i++)
224 r.accumulate(pressurespace,i,y*pbasis[i]*factor);
229 template<
typename IG,
typename LFSV,
typename R>
233 using VelocitySpace =
typename LFSV::template Child<0>::Type;
234 using DF =
typename VelocitySpace::Traits::FiniteElementType::
235 Traits::LocalBasisType::Traits::DomainFieldType;
236 using VelocityRangeType =
typename VelocitySpace::Traits::FiniteElementType::
237 Traits::LocalBasisType::Traits::RangeType;
240 using namespace Indices;
241 const auto& velocityspace = child(lfsv,_0);
244 const int dim = IG::Entity::dimension;
247 const auto& cell_inside =
ig.inside();
250 auto geo =
ig.geometry();
251 auto geo_inside = cell_inside.geometry();
254 auto geo_in_inside =
ig.geometryInInside();
257 Dune::FieldVector<DF,dim> pos;
259 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
260 jac = geo_inside.jacobianInverseTransposed(pos);
262 auto det = geo_inside.integrationElement(pos);
265 std::vector<VelocityRangeType> vbasis(velocityspace.size());
266 std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size());
272 auto bctype = param.bctype(
ig.intersection(),ip.position());
279 auto local = geo_in_inside.global(ip.position());
282 velocityspace.finiteElement().localBasis().evaluateFunction(local,vbasis);
285 for (std::size_t i=0; i<velocityspace.size(); i++)
287 vtransformedbasis[i] = 0.0;
288 jac.umtv(vbasis[i],vtransformedbasis[i]);
292 auto y = param.g(cell_inside,local);
295 auto factor = ip.weight()*geo.integrationElement(ip.position())/det;
296 for (std::size_t i=0; i<velocityspace.size(); i++)
297 r.accumulate(velocityspace,i,y*(vtransformedbasis[i]*
ig.unitOuterNormal(ip.position()))*factor);
static const int dim
Definition: adaptivity.hh:84
const IG & ig
Definition: constraints.hh:149
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
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
Type
Definition: convectiondiffusionparameter.hh:113
@ Neumann
Definition: convectiondiffusionparameter.hh:113
Definition: diffusionmixed.hh:41
@ doLambdaBoundary
Definition: diffusionmixed.hh:52
@ doPatternVolume
Definition: diffusionmixed.hh:47
DiffusionMixed(const PARAM ¶m_, int qorder_v_=2, int qorder_p_=1)
Definition: diffusionmixed.hh:54
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:192
@ doLambdaVolume
Definition: diffusionmixed.hh:51
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:230
@ doAlphaVolume
Definition: diffusionmixed.hh:50
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:65
Default flags for all local operators.
Definition: flags.hh:19
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:32
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
sparsity pattern generator
Definition: pattern.hh:14