dune-pdelab  2.5-dev
navierstokesmass.hh
Go to the documentation of this file.
1 // -*- tab-width: 2; indent-tabs-mode: nil -*-
2 // vi: set et ts=2 sw=2 sts=2:
3 
4 #ifndef DUNE_PDELAB_LOCALOPERATOR_NAVIERSTOKESMASS_HH
5 #define DUNE_PDELAB_LOCALOPERATOR_NAVIERSTOKESMASS_HH
6 
7 #include <dune/geometry/quadraturerules.hh>
8 #include <dune/localfunctions/common/interfaceswitch.hh>
9 
10 #include <dune/typetree/compositenode.hh>
11 #include <dune/typetree/utility.hh>
12 
16 
19 
20 namespace Dune {
21  namespace PDELab {
22 
30  template<typename PRM>
32  public FullVolumePattern ,
34  public InstationaryLocalOperatorDefaultMethods<typename PRM::Traits::RangeField>
35  {
36  public:
37  // pattern assembly flags
38  enum { doPatternVolume = true };
39 
40  // residual assembly flags
41  enum { doAlphaVolume = true };
42 
43  NavierStokesMass (const PRM & p_, int superintegration_order_ = 0)
44  : p(p_), superintegration_order(superintegration_order_)
45  {}
46 
47  // volume integral depending on test and ansatz functions
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
50  {
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)
54  {
55  scalar_alpha_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),r);
56  }
57  }
58 
59  // jacobian of volume term
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,
62  M& mat) const
63  {
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)
67  {
68  scalar_jacobian_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),mat);
69  }
70  }
71 
72  private:
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,
75  R& r) const
76  {
77 
78  // Switches between local and global interface
79  using FESwitch = FiniteElementInterfaceSwitch<
80  typename LFSU::Traits::FiniteElementType>;
81  using BasisSwitch = BasisInterfaceSwitch<
82  typename FESwitch::Basis>;
83 
84  // Define types
85  using RF = typename BasisSwitch::RangeField;
86  using RangeType = typename BasisSwitch::Range;
87  using size_type = typename LFSU::Traits::SizeType;
88 
89  // Dimensions
90  const int dim = EG::Geometry::mydimension;
91 
92  // Get geometry
93  auto geo = eg.geometry();
94 
95  // Determine quadrature order
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;
99 
100  // Initialize vectors outside for loop
101  std::vector<RangeType> phi(lfsu.size());
102 
103  // Loop over quadrature points
104  for (const auto& ip : quadratureRule(geo,qorder))
105  {
106  // evaluate basis functions
107  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(ip.position(),phi);
108 
109  auto rho = p.rho(eg,ip.position());
110  // evaluate u
111  RF u=0.0;
112  for (size_type i=0; i<lfsu.size(); i++)
113  u += x(lfsu,i)*phi[i];
114 
115  // u*phi_i
116  auto factor = ip.weight() * rho * geo.integrationElement(ip.position());
117 
118  for (size_type i=0; i<lfsu.size(); i++)
119  r.accumulate(lfsv,i, u*phi[i]*factor);
120  }
121  }
122 
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,
125  M& mat) const
126  {
127 
128  // Switches between local and global interface
129  using FESwitch = FiniteElementInterfaceSwitch<
130  typename LFSU::Traits::FiniteElementType>;
131  using BasisSwitch = BasisInterfaceSwitch<
132  typename FESwitch::Basis>;
133 
134  // Define types
135  using RangeType = typename BasisSwitch::Range;
136  using size_type = typename LFSU::Traits::SizeType;
137 
138  // Dimensions
139  const int dim = EG::Geometry::mydimension;
140 
141  // Get geometry
142  auto geo = eg.geometry();
143 
144  // Determine quadrature order
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;
148 
149  // Initialize vectors outside for loop
150  std::vector<RangeType> phi(lfsu.size());
151 
152  // Loop over quadrature points
153  for (const auto& ip : quadratureRule(geo,qorder))
154  {
155  // evaluate basis functions
156  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(ip.position(),phi);
157 
158  // integrate phi_j*phi_i
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);
164  }
165  }
166 
167  const PRM& p;
168  const int superintegration_order;
169  }; // end class NavierStokesMass
170 
179  template<typename PRM>
181  public FullVolumePattern ,
184  {
185  public:
186  // pattern assembly flags
187  enum { doPatternVolume = true };
188 
189  // residual assembly flags
190  enum { doAlphaVolume = true };
191 
192  NavierStokesVelVecMass (const PRM & p_, int superintegration_order_ = 0)
193  : p(p_), superintegration_order(superintegration_order_)
194  {}
195 
196  // volume integral depending on test and ansatz functions
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
199  {
200  // subspaces
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);
205 
206  // Define types
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;
211 
212  // dimensions
213  const int dim = EG::Geometry::mydimension;
214 
215  // Get geometry
216  auto geo = eg.geometry();
217 
218  // Determine quadrature order
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;
222 
223  // Initialize vectors outside for loop
224  std::vector<Range_V> phi_v(lfsv_v.size());
225  Range_V u(0.0);
226 
227  // loop over quadrature points
228  for (const auto& ip : quadratureRule(geo,qorder))
229  {
230  auto local = ip.position();
231  auto rho = p.rho(eg,local);
232 
233  // compute basis functions
234  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
235 
236  // compute u
237  u = 0.0;
238  for(size_type i=0; i<lfsu_v.size(); i++)
239  u.axpy(x(lfsu_v,i),phi_v[i]);
240 
241  auto factor = ip.weight() * rho * geo.integrationElement(ip.position());
242 
243  for(size_type i=0; i<lfsv_v.size(); i++)
244  r.accumulate(lfsv_v,i, (u*phi_v[i]) * factor);
245 
246  } // end loop quadrature points
247  } // end alpha_volume
248 
249  // jacobian of volume term
250  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
251  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
252  M& mat) const
253  {
254  // subspaces
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);
259 
260  // Define types
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;
265 
266  // dimensions
267  const int dim = EG::Geometry::mydimension;
268 
269  // Get geometry
270  auto geo = eg.geometry();
271 
272  // Determine quadrature order
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;
276 
277  // Initialize vectors outside for loop
278  std::vector<Range_V> phi_v(lfsv_v.size());
279 
280  // Loop over quadrature points
281  for (const auto& ip : quadratureRule(geo,qorder))
282  {
283  auto local = ip.position();
284  auto rho = p.rho(eg,local);
285 
286  // compute basis functions
287  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
288 
289  auto factor = ip.weight() * rho * geo.integrationElement(ip.position());
290 
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);
294  } // end loop quadrature points
295  } // end jacobian_volume
296 
297  private :
298  const PRM& p;
299  const int superintegration_order;
300  }; // end class NavierStokesVelVecMass
301 
302  } // end namespace PDELab
303 } // end namespace Dune
304 #endif // DUNE_PDELAB_LOCALOPERATOR_NAVIERSTOKESMASS_HH
static const int dim
Definition: adaptivity.hh:83
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
NavierStokesVelVecMass(const PRM &p_, int superintegration_order_=0)
Definition: navierstokesmass.hh:192
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
NavierStokesMass(const PRM &p_, int superintegration_order_=0)
Definition: navierstokesmass.hh:43
sparsity pattern generator
Definition: pattern.hh:13
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
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
Definition: navierstokesmass.hh:41
Default flags for all local operators.
Definition: flags.hh:18
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: navierstokesmass.hh:61
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: navierstokesmass.hh:251
const P & p
Definition: constraints.hh:147
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: navierstokesmass.hh:198