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>
37 template<
typename PARAM>
65 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
66 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
69 using VelocitySpace =
typename LFSU::template Child<0>::Type;
70 using PressureSpace =
typename LFSU::template Child<1>::Type;
71 using DF =
typename VelocitySpace::Traits::FiniteElementType::
72 Traits::LocalBasisType::Traits::DomainFieldType;
73 using RF =
typename VelocitySpace::Traits::FiniteElementType::
74 Traits::LocalBasisType::Traits::RangeFieldType;
75 using VelocityJacobianType =
typename VelocitySpace::Traits::FiniteElementType::
76 Traits::LocalBasisType::Traits::JacobianType;
77 using VelocityRangeType =
typename VelocitySpace::Traits::FiniteElementType::
78 Traits::LocalBasisType::Traits::RangeType;
79 using PressureRangeType =
typename PressureSpace::Traits::FiniteElementType::
80 Traits::LocalBasisType::Traits::RangeType;
83 using namespace Indices;
84 const auto& velocityspace = child(lfsu,_0);
85 const auto& pressurespace = lfsu.template child<1>();
88 const int dim = EG::Geometry::mydimension;
91 const auto& cell = eg.entity();
94 auto geo = eg.geometry();
97 Dune::FieldVector<DF,dim> pos;
99 auto jac = geo.jacobianInverseTransposed(pos);
101 auto det = geo.integrationElement(pos);
104 auto ref_el = referenceElement(geo);
105 auto localcenter = ref_el.position(0,0);
106 auto tensor = param.A(cell,localcenter);
110 std::vector<VelocityRangeType> vbasis(velocityspace.size());
111 std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size());
112 VelocityRangeType sigma;
113 VelocityRangeType Kinvsigma;
114 std::vector<VelocityJacobianType> vjacbasis(velocityspace.size());
115 std::vector<PressureRangeType> pbasis(pressurespace.size());
116 std::vector<RF> divergence(velocityspace.size(),0.0);
124 velocityspace.finiteElement().localBasis().evaluateFunction(ip.position(),vbasis);
127 for (std::size_t i=0; i<velocityspace.size(); i++)
129 vtransformedbasis[i] = 0.0;
130 jac.umtv(vbasis[i],vtransformedbasis[i]);
135 for (std::size_t i=0; i<velocityspace.size(); i++)
136 sigma.axpy(x(velocityspace,i),vtransformedbasis[i]);
139 tensor.mv(sigma,Kinvsigma);
142 auto factor = ip.weight() / det;
143 for (std::size_t i=0; i<velocityspace.size(); i++)
144 r.accumulate(velocityspace,i,(Kinvsigma*vtransformedbasis[i])*factor);
152 velocityspace.finiteElement().localBasis().evaluateJacobian(ip.position(),vjacbasis);
153 pressurespace.finiteElement().localBasis().evaluateFunction(ip.position(),pbasis);
158 for (std::size_t i=0; i<pressurespace.size(); i++)
159 u.axpy(x(pressurespace,i),pbasis[i]);
162 auto a0value = param.c(cell,ip.position());
165 RF factor = ip.weight();
166 for (std::size_t i=0; i<pressurespace.size(); i++)
167 r.accumulate(pressurespace,i,-a0value*u*pbasis[i]*factor);
170 for (std::size_t i=0; i<velocityspace.size(); i++){
172 for (
int j=0; j<
dim; j++)
173 divergence[i] += vjacbasis[i][j][j];
177 for (std::size_t i=0; i<velocityspace.size(); i++)
178 r.accumulate(velocityspace,i,-u*divergence[i]*factor);
181 RF divergencesigma = 0.0;
182 for (std::size_t i=0; i<velocityspace.size(); i++)
183 divergencesigma += x(velocityspace,i)*divergence[i];
186 for (std::size_t i=0; i<pressurespace.size(); i++)
187 r.accumulate(pressurespace,i,-divergencesigma*pbasis[i]*factor);
192 template<
typename EG,
typename LFSV,
typename R>
196 using PressureSpace =
typename LFSV::template Child<1>::Type;
197 using PressureRangeType =
typename PressureSpace::Traits::FiniteElementType::
198 Traits::LocalBasisType::Traits::RangeType;
201 using namespace Indices;
202 const auto& pressurespace = child(lfsv,_1);
205 const auto& cell = eg.entity();
208 auto geo = eg.geometry();
211 std::vector<PressureRangeType> pbasis(pressurespace.size());
217 pressurespace.finiteElement().localBasis().evaluateFunction(ip.position(),pbasis);
220 auto y = param.f(cell,ip.position());
223 auto factor = ip.weight() * geo.integrationElement(ip.position());
224 for (std::size_t i=0; i<pressurespace.size(); i++)
225 r.accumulate(pressurespace,i,y*pbasis[i]*factor);
230 template<
typename IG,
typename LFSV,
typename R>
234 using VelocitySpace =
typename LFSV::template Child<0>::Type;
235 using DF =
typename VelocitySpace::Traits::FiniteElementType::
236 Traits::LocalBasisType::Traits::DomainFieldType;
237 using VelocityRangeType =
typename VelocitySpace::Traits::FiniteElementType::
238 Traits::LocalBasisType::Traits::RangeType;
241 using namespace Indices;
242 const auto& velocityspace = child(lfsv,_0);
245 const int dim = IG::Entity::dimension;
248 const auto& cell_inside =
ig.inside();
251 auto geo =
ig.geometry();
252 auto geo_inside = cell_inside.geometry();
255 auto geo_in_inside =
ig.geometryInInside();
258 Dune::FieldVector<DF,dim> pos;
260 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
261 jac = geo_inside.jacobianInverseTransposed(pos);
263 auto det = geo_inside.integrationElement(pos);
266 std::vector<VelocityRangeType> vbasis(velocityspace.size());
267 std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size());
273 auto bctype = param.bctype(
ig.intersection(),ip.position());
280 auto local = geo_in_inside.global(ip.position());
283 velocityspace.finiteElement().localBasis().evaluateFunction(local,vbasis);
286 for (std::size_t i=0; i<velocityspace.size(); i++)
288 vtransformedbasis[i] = 0.0;
289 jac.umtv(vbasis[i],vtransformedbasis[i]);
293 auto y = param.g(cell_inside,local);
296 auto factor = ip.weight()*geo.integrationElement(ip.position())/det;
297 for (std::size_t i=0; i<velocityspace.size(); i++)
298 r.accumulate(velocityspace,i,y*(vtransformedbasis[i]*
ig.unitOuterNormal(ip.position()))*factor);
312 #endif // DUNE_PDELAB_LOCALOPERATOR_DIFFUSIONMIXED_HH