4 #ifndef DUNE_PDELAB_LOCALOPERATOR_NAVIERSTOKESMASS_HH 5 #define DUNE_PDELAB_LOCALOPERATOR_NAVIERSTOKESMASS_HH 7 #include <dune/geometry/quadraturerules.hh> 8 #include <dune/localfunctions/common/interfaceswitch.hh> 10 #include <dune/typetree/compositenode.hh> 11 #include <dune/typetree/utility.hh> 30 template<
typename PRM>
44 : p(p_), superintegration_order(superintegration_order_)
48 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
49 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 51 using namespace Indices;
52 const auto& lfsv_pfs_v = child(lfsv,_0);
53 for(
unsigned int i=0; i<TypeTree::degree(lfsv_pfs_v); ++i)
55 scalar_alpha_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),r);
60 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
61 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
64 using namespace Indices;
65 const auto& lfsv_pfs_v = child(lfsv,_0);
66 for(
unsigned int i=0; i<TypeTree::degree(lfsv_pfs_v); ++i)
68 scalar_jacobian_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),mat);
73 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
74 void scalar_alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
79 using FESwitch = FiniteElementInterfaceSwitch<
80 typename LFSU::Traits::FiniteElementType>;
81 using BasisSwitch = BasisInterfaceSwitch<
82 typename FESwitch::Basis>;
85 using RF =
typename BasisSwitch::RangeField;
86 using RangeType =
typename BasisSwitch::Range;
87 using size_type =
typename LFSU::Traits::SizeType;
90 const int dim = EG::Geometry::mydimension;
93 auto geo = eg.geometry();
96 const int v_order = FESwitch::basis(lfsu.finiteElement()).order();
97 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
98 const int qorder = 2*v_order + det_jac_order + superintegration_order;
101 std::vector<RangeType> phi(lfsu.size());
107 FESwitch::basis(lfsu.finiteElement()).evaluateFunction(ip.position(),phi);
109 auto rho =
p.rho(eg,ip.position());
112 for (size_type i=0; i<lfsu.size(); i++)
113 u += x(lfsu,i)*phi[i];
116 auto factor = ip.weight() * rho * geo.integrationElement(ip.position());
118 for (size_type i=0; i<lfsu.size(); i++)
119 r.accumulate(lfsv,i, u*phi[i]*factor);
123 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
124 void scalar_jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
129 using FESwitch = FiniteElementInterfaceSwitch<
130 typename LFSU::Traits::FiniteElementType>;
131 using BasisSwitch = BasisInterfaceSwitch<
132 typename FESwitch::Basis>;
135 using RangeType =
typename BasisSwitch::Range;
136 using size_type =
typename LFSU::Traits::SizeType;
139 const int dim = EG::Geometry::mydimension;
142 auto geo = eg.geometry();
145 const int v_order = FESwitch::basis(lfsu.finiteElement()).order();
146 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
147 const int qorder = 2*v_order + det_jac_order + superintegration_order;
150 std::vector<RangeType> phi(lfsu.size());
156 FESwitch::basis(lfsu.finiteElement()).evaluateFunction(ip.position(),phi);
159 auto rho =
p.rho(eg,ip.position());
160 auto factor = ip.weight() * rho * geo.integrationElement(ip.position());
161 for (size_type j=0; j<lfsu.size(); j++)
162 for (size_type i=0; i<lfsu.size(); i++)
163 mat.accumulate(lfsv,i,lfsu,j, phi[j]*phi[i]*factor);
168 const int superintegration_order;
179 template<
typename PRM>
193 :
p(p_), superintegration_order(superintegration_order_)
197 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
198 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 201 using namespace Indices;
202 using LFSV_V = TypeTree::Child<LFSV,_0>;
203 const auto& lfsv_v = child(lfsv,_0);
204 const auto& lfsu_v = child(lfsu,_0);
207 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
208 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
209 using Range_V =
typename BasisSwitch_V::Range;
210 using size_type =
typename LFSV::Traits::SizeType;
213 const int dim = EG::Geometry::mydimension;
216 auto geo = eg.geometry();
219 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
220 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
221 const int qorder = 2*v_order + det_jac_order + superintegration_order;
224 std::vector<Range_V> phi_v(lfsv_v.size());
230 auto local = ip.position();
231 auto rho =
p.rho(eg,local);
234 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
238 for(size_type i=0; i<lfsu_v.size(); i++)
239 u.axpy(x(lfsu_v,i),phi_v[i]);
241 auto factor = ip.weight() * rho * geo.integrationElement(ip.position());
243 for(size_type i=0; i<lfsv_v.size(); i++)
244 r.accumulate(lfsv_v,i, (u*phi_v[i]) * factor);
250 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
255 using namespace Indices;
256 using LFSV_V = TypeTree::Child<LFSV,_0>;
257 const auto& lfsv_v = child(lfsv,_0);
258 const auto& lfsu_v = child(lfsu,_0);
261 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
262 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
263 using Range_V =
typename BasisSwitch_V::Range;
264 using size_type =
typename LFSV::Traits::SizeType;
267 const int dim = EG::Geometry::mydimension;
270 auto geo = eg.geometry();
273 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
274 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
275 const int qorder = 2*v_order + det_jac_order + superintegration_order;
278 std::vector<Range_V> phi_v(lfsv_v.size());
283 auto local = ip.position();
284 auto rho =
p.rho(eg,local);
287 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
289 auto factor = ip.weight() * rho * geo.integrationElement(ip.position());
291 for(size_type i=0; i<lfsv_v.size(); i++)
292 for(size_type j=0; j<lfsu_v.size(); j++)
293 mat.accumulate(lfsv_v,i,lfsu_v,j, (phi_v[j]*phi_v[i]) * factor);
299 const int superintegration_order;
304 #endif // DUNE_PDELAB_LOCALOPERATOR_NAVIERSTOKESMASS_HH 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
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: navierstokesmass.hh:49
static const int dim
Definition: adaptivity.hh:83
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: navierstokesmass.hh:198
A local operator for the mass term corresponding to the instationary part in the Navier-Stokes equati...
Definition: navierstokesmass.hh:31
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
NavierStokesMass(const PRM &p_, int superintegration_order_=0)
Definition: navierstokesmass.hh:43
const P & p
Definition: constraints.hh:147
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
sparsity pattern generator
Definition: pattern.hh:13
Definition: navierstokesmass.hh:38
A local operator for the mass term corresponding to the instationary part in the Navier-Stokes equati...
Definition: navierstokesmass.hh:180
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: navierstokesmass.hh:61
NavierStokesVelVecMass(const PRM &p_, int superintegration_order_=0)
Definition: navierstokesmass.hh:192
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: navierstokesmass.hh:251
Default flags for all local operators.
Definition: flags.hh:18
Definition: navierstokesmass.hh:41