4#ifndef DUNE_PDELAB_LOCALOPERATOR_NAVIERSTOKESMASS_HH
5#define DUNE_PDELAB_LOCALOPERATOR_NAVIERSTOKESMASS_HH
7#include <dune/geometry/referenceelements.hh>
8#include <dune/localfunctions/common/interfaceswitch.hh>
10#include <dune/typetree/compositenode.hh>
11#include <dune/typetree/utility.hh>
29 template<
typename PRM>
43 : p(p_), superintegration_order(superintegration_order_)
47 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
48 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
50 using namespace Indices;
51 const auto& lfsv_pfs_v = child(lfsv,_0);
52 for(
unsigned int i=0; i<TypeTree::degree(lfsv_pfs_v); ++i)
54 scalar_alpha_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),r);
59 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
60 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
63 using namespace Indices;
64 const auto& lfsv_pfs_v = child(lfsv,_0);
65 for(
unsigned int i=0; i<TypeTree::degree(lfsv_pfs_v); ++i)
67 scalar_jacobian_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),mat);
72 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
73 void scalar_alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
78 using FESwitch = FiniteElementInterfaceSwitch<
79 typename LFSU::Traits::FiniteElementType>;
80 using BasisSwitch = BasisInterfaceSwitch<
81 typename FESwitch::Basis>;
84 using RF =
typename BasisSwitch::RangeField;
85 using RangeType =
typename BasisSwitch::Range;
86 using size_type =
typename LFSU::Traits::SizeType;
89 const int dim = EG::Geometry::mydimension;
92 auto geo = eg.geometry();
95 const int v_order = FESwitch::basis(lfsu.finiteElement()).order();
96 const int det_jac_order = geo.type().isSimplex() ? 0 : (
dim-1);
97 const int qorder = 2*v_order + det_jac_order + superintegration_order;
100 std::vector<RangeType> phi(lfsu.size());
106 FESwitch::basis(lfsu.finiteElement()).evaluateFunction(ip.position(),phi);
108 auto rho = p.rho(eg,ip.position());
111 for (size_type i=0; i<lfsu.size(); i++)
112 u += x(lfsu,i)*phi[i];
115 auto factor = ip.weight() * rho * geo.integrationElement(ip.position());
117 for (size_type i=0; i<lfsu.size(); i++)
118 r.accumulate(lfsv,i, u*phi[i]*factor);
122 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
123 void scalar_jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
128 using FESwitch = FiniteElementInterfaceSwitch<
129 typename LFSU::Traits::FiniteElementType>;
130 using BasisSwitch = BasisInterfaceSwitch<
131 typename FESwitch::Basis>;
134 using RangeType =
typename BasisSwitch::Range;
135 using size_type =
typename LFSU::Traits::SizeType;
138 const int dim = EG::Geometry::mydimension;
141 auto geo = eg.geometry();
144 const int v_order = FESwitch::basis(lfsu.finiteElement()).order();
145 const int det_jac_order = geo.type().isSimplex() ? 0 : (
dim-1);
146 const int qorder = 2*v_order + det_jac_order + superintegration_order;
149 std::vector<RangeType> phi(lfsu.size());
152 for (
const auto& ip : quadratureRule(geo,qorder))
155 FESwitch::basis(lfsu.finiteElement()).evaluateFunction(ip.position(),phi);
158 auto rho = p.rho(eg,ip.position());
159 auto factor = ip.weight() * rho * geo.integrationElement(ip.position());
160 for (size_type j=0; j<lfsu.size(); j++)
161 for (size_type i=0; i<lfsu.size(); i++)
162 mat.accumulate(lfsv,i,lfsu,j, phi[j]*phi[i]*factor);
167 const int superintegration_order;
178 template<
typename PRM>
192 : p(p_), superintegration_order(superintegration_order_)
196 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
197 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
200 using namespace Indices;
201 using LFSV_V = TypeTree::Child<LFSV,_0>;
202 const auto& lfsv_v = child(lfsv,_0);
203 const auto& lfsu_v = child(lfsu,_0);
206 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
207 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
208 using Range_V =
typename BasisSwitch_V::Range;
209 using size_type =
typename LFSV::Traits::SizeType;
212 const int dim = EG::Geometry::mydimension;
215 auto geo = eg.geometry();
218 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
219 const int det_jac_order = geo.type().isSimplex() ? 0 : (
dim-1);
220 const int qorder = 2*v_order + det_jac_order + superintegration_order;
223 std::vector<Range_V> phi_v(lfsv_v.size());
229 auto local = ip.position();
230 auto rho = p.rho(eg,local);
233 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
237 for(size_type i=0; i<lfsu_v.size(); i++)
238 u.axpy(x(lfsu_v,i),phi_v[i]);
240 auto factor = ip.weight() * rho * geo.integrationElement(ip.position());
242 for(size_type i=0; i<lfsv_v.size(); i++)
243 r.accumulate(lfsv_v,i, (u*phi_v[i]) * factor);
249 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
254 using namespace Indices;
255 using LFSV_V = TypeTree::Child<LFSV,_0>;
256 const auto& lfsv_v = child(lfsv,_0);
257 const auto& lfsu_v = child(lfsu,_0);
260 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
261 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
262 using Range_V =
typename BasisSwitch_V::Range;
263 using size_type =
typename LFSV::Traits::SizeType;
266 const int dim = EG::Geometry::mydimension;
269 auto geo = eg.geometry();
272 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
273 const int det_jac_order = geo.type().isSimplex() ? 0 : (
dim-1);
274 const int qorder = 2*v_order + det_jac_order + superintegration_order;
277 std::vector<Range_V> phi_v(lfsv_v.size());
282 auto local = ip.position();
283 auto rho = p.rho(eg,local);
286 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
288 auto factor = ip.weight() * rho * geo.integrationElement(ip.position());
290 for(size_type i=0; i<lfsv_v.size(); i++)
291 for(size_type j=0; j<lfsu_v.size(); j++)
292 mat.accumulate(lfsv_v,i,lfsu_v,j, (phi_v[j]*phi_v[i]) * factor);
298 const int superintegration_order;
static const int dim
Definition adaptivity.hh:84
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
Default flags for all local operators.
Definition flags.hh:19
Default class for additional methods in instationary local operators.
Definition idefault.hh:90
A local operator for the mass term corresponding to the instationary part in the Navier-Stokes equati...
Definition navierstokesmass.hh:34
@ doAlphaVolume
Definition navierstokesmass.hh:40
@ doPatternVolume
Definition navierstokesmass.hh:37
NavierStokesMass(const PRM &p_, int superintegration_order_=0)
Definition navierstokesmass.hh:42
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition navierstokesmass.hh:48
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition navierstokesmass.hh:60
A local operator for the mass term corresponding to the instationary part in the Navier-Stokes equati...
Definition navierstokesmass.hh:183
@ doPatternVolume
Definition navierstokesmass.hh:186
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition navierstokesmass.hh:197
@ doAlphaVolume
Definition navierstokesmass.hh:189
NavierStokesVelVecMass(const PRM &p_, int superintegration_order_=0)
Definition navierstokesmass.hh:191
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition navierstokesmass.hh:250
sparsity pattern generator
Definition pattern.hh:14