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>
13#include <dune/pdelab/localoperator/pattern.hh>
14#include <dune/pdelab/localoperator/flags.hh>
15#include <dune/pdelab/localoperator/idefault.hh>
17#include <dune/pdelab/common/quadraturerules.hh>
29 template<
typename PRM>
37 enum { doPatternVolume =
true };
40 enum { doAlphaVolume =
true };
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
51 const auto& lfsv_pfs_v =
child(lfsv,
_0);
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,
64 const auto& lfsv_pfs_v =
child(lfsv,
_0);
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,
79 typename LFSU::Traits::FiniteElementType>;
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());
103 for (
const auto& ip : quadratureRule(geo,qorder))
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,
129 typename LFSU::Traits::FiniteElementType>;
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>
186 enum { doPatternVolume =
true };
189 enum { doAlphaVolume =
true };
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
202 const auto& lfsv_v =
child(lfsv,
_0);
203 const auto& lfsu_v =
child(lfsu,
_0);
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());
227 for (
const auto& ip : quadratureRule(geo,qorder))
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>
250 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
256 const auto& lfsv_v =
child(lfsv,
_0);
257 const auto& lfsu_v =
child(lfsu,
_0);
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());
280 for (
const auto& ip : quadratureRule(geo,qorder))
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;
sparsity pattern generator
Definition: pattern.hh:14
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
Default flags for all local operators.
Definition: flags.hh:19
A local operator for the mass term corresponding to the instationary part in the Navier-Stokes equati...
Definition: navierstokesmass.hh:34
A local operator for the mass term corresponding to the instationary part in the Navier-Stokes equati...
Definition: navierstokesmass.hh:183
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:52
std::size_t degree(const Node &node)
Returns the degree of node as run time information.
Definition: nodeinterface.hh:79
Namespace with predefined compile time indices for the range [0,19].
Definition: indices.hh:50
Dune namespace.
Definition: alignedallocator.hh:13
Switch for uniform treatment of local and global basis classes.
Definition: interfaceswitch.hh:154
Switch for uniform treatment of finite element with either the local or the global interface.
Definition: interfaceswitch.hh:30