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