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>
24 template<
typename PRM>
38 : p(p_), superintegration_order(superintegration_order_)
42 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
43 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
45 typedef typename LFSV::template Child<0>::Type LFSV_PFS_V;
46 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<0>();
48 for(
unsigned int i=0; i<LFSV_PFS_V::CHILDREN; ++i)
50 scalar_alpha_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),r);
55 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
59 typedef typename LFSV::template Child<0>::Type LFSV_PFS_V;
60 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<0>();
62 for(
unsigned int i=0; i<LFSV_PFS_V::CHILDREN; ++i)
64 scalar_jacobian_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),mat);
69 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
70 void scalar_alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
75 typedef FiniteElementInterfaceSwitch<
76 typename LFSU::Traits::FiniteElementType
78 typedef BasisInterfaceSwitch<
79 typename FESwitch::Basis
83 typedef typename BasisSwitch::DomainField DF;
84 typedef typename BasisSwitch::RangeField RF;
85 typedef typename BasisSwitch::Range RangeType;
87 typedef typename LFSU::Traits::SizeType size_type;
90 const int dim = EG::Geometry::dimension;
93 Dune::GeometryType gt = eg.geometry().type();
94 const int v_order = FESwitch::basis(lfsu.finiteElement()).order();
95 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
96 const int qorder = 2*v_order + det_jac_order + superintegration_order;
97 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
100 for (
const auto& ip : rule)
103 std::vector<RangeType> phi(lfsu.size());
104 FESwitch::basis(lfsu.finiteElement()).evaluateFunction(ip.position(),phi);
106 RF rho = p.rho(eg,ip.position());
109 for (size_type i=0; i<lfsu.size(); i++)
110 u += x(lfsu,i)*phi[i];
113 RF factor = ip.weight() * rho * eg.geometry().integrationElement(ip.position());
115 for (size_type i=0; i<lfsu.size(); i++)
116 r.accumulate(lfsv,i, u*phi[i]*factor);
120 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
121 void scalar_jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
126 typedef FiniteElementInterfaceSwitch<
127 typename LFSU::Traits::FiniteElementType
129 typedef BasisInterfaceSwitch<
130 typename FESwitch::Basis
134 typedef typename BasisSwitch::DomainField DF;
135 typedef typename BasisSwitch::RangeField RF;
136 typedef typename BasisSwitch::Range RangeType;
137 typedef typename LFSU::Traits::SizeType size_type;
140 const int dim = EG::Geometry::dimension;
143 Dune::GeometryType gt = eg.geometry().type();
144 const int v_order = FESwitch::basis(lfsu.finiteElement()).order();
145 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
146 const int qorder = 2*v_order + det_jac_order + superintegration_order;
147 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
150 for (
const auto& ip : rule)
153 std::vector<RangeType> phi(lfsu.size());
154 FESwitch::basis(lfsu.finiteElement()).evaluateFunction(ip.position(),phi);
157 RF rho = p.rho(eg,ip.position());
158 RF factor = ip.weight() * rho * eg.geometry().integrationElement(ip.position());
159 for (size_type j=0; j<lfsu.size(); j++)
160 for (size_type i=0; i<lfsu.size(); i++)
161 mat.accumulate(lfsv,i,lfsu,j, phi[j]*phi[i]*factor);
166 const int superintegration_order;
177 template<
typename PRM>
191 : p(p_), superintegration_order(superintegration_order_)
195 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
196 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
199 const int dim = EG::Geometry::dimension;
202 typedef typename LFSV::template Child<0>::Type LFSV_V;
203 const LFSV_V& lfsv_v = lfsv.template child<0>();
204 const LFSV_V& lfsu_v = lfsu.template child<0>();
207 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
208 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
209 typedef typename BasisSwitch_V::DomainField DF;
210 typedef typename BasisSwitch_V::RangeField RF;
211 typedef typename BasisSwitch_V::Range Range_V;
212 typedef typename LFSV::Traits::SizeType size_type;
215 Dune::GeometryType gt = eg.geometry().type();
216 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
217 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
218 const int qorder = 2*v_order + det_jac_order + superintegration_order;
219 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
222 for (
const auto& ip : rule)
224 const Dune::FieldVector<DF,dim> local = ip.position();
225 const RF rho = p.rho(eg,local);
228 std::vector<Range_V> phi_v(lfsv_v.size());
229 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
233 for(size_type i=0; i<lfsu_v.size(); i++)
234 u.axpy(x(lfsu_v,i),phi_v[i]);
236 const RF factor = ip.weight() * rho * eg.geometry().integrationElement(ip.position());
238 for(size_type i=0; i<lfsv_v.size(); i++)
239 r.accumulate(lfsv_v,i, (u*phi_v[i]) * factor);
245 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
250 const int dim = EG::Geometry::dimension;
253 typedef typename LFSV::template Child<0>::Type LFSV_V;
254 const LFSV_V& lfsv_v = lfsv.template child<0>();
255 const LFSV_V& lfsu_v = lfsu.template child<0>();
258 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
259 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
260 typedef typename BasisSwitch_V::DomainField DF;
261 typedef typename BasisSwitch_V::RangeField RF;
262 typedef typename BasisSwitch_V::Range Range_V;
263 typedef typename LFSV::Traits::SizeType size_type;
266 Dune::GeometryType gt = eg.geometry().type();
267 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
268 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
269 const int qorder = 2*v_order + det_jac_order + superintegration_order;
270 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
273 for (
const auto& ip : rule)
275 const Dune::FieldVector<DF,dim> local = ip.position();
276 const RF rho = p.rho(eg,local);
279 std::vector<Range_V> phi_v(lfsv_v.size());
280 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
282 const RF factor = ip.weight() * rho * eg.geometry().integrationElement(ip.position());
284 for(size_type i=0; i<lfsv_v.size(); i++)
285 for(size_type j=0; j<lfsu_v.size(); j++)
286 mat.accumulate(lfsv_v,i,lfsu_v,j, (phi_v[j]*phi_v[i]) * factor);
292 const int superintegration_order;
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: navierstokesmass.hh:196
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: navierstokesmass.hh:56
NavierStokesVelVecMass(const PRM &p_, int superintegration_order_=0)
Definition: navierstokesmass.hh:190
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: navierstokesmass.hh:246
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: navierstokesmass.hh:43
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
NavierStokesMass(const PRM &p_, int superintegration_order_=0)
Definition: navierstokesmass.hh:37
Definition: navierstokesmass.hh:32
Definition: adaptivity.hh:27
A local operator for the mass term corresponding to the instationary part in the Navier-Stokes equati...
Definition: navierstokesmass.hh:25
static const int dim
Definition: adaptivity.hh:83
Definition: navierstokesmass.hh:185
A local operator for the mass term corresponding to the instationary part in the Navier-Stokes equati...
Definition: navierstokesmass.hh:178
Definition: navierstokesmass.hh:35
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Definition: navierstokesmass.hh:188