dune-pdelab  2.4-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 
13 
14 namespace Dune {
15  namespace PDELab {
16 
24  template<typename PRM>
26  public FullVolumePattern ,
29  {
30  public:
31  // pattern assembly flags
32  enum { doPatternVolume = true };
33 
34  // residual assembly flags
35  enum { doAlphaVolume = true };
36 
37  NavierStokesMass (const PRM & p_, int superintegration_order_ = 0)
38  : p(p_), superintegration_order(superintegration_order_)
39  {}
40 
41  // volume integral depending on test and ansatz functions
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
44  {
45  typedef typename LFSV::template Child<0>::Type LFSV_PFS_V;
46  const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<0>();
47 
48  for(unsigned int i=0; i<LFSV_PFS_V::CHILDREN; ++i)
49  {
50  scalar_alpha_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),r);
51  }
52  }
53 
54  // jacobian of volume term
55  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
56  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
57  M& mat) const
58  {
59  typedef typename LFSV::template Child<0>::Type LFSV_PFS_V;
60  const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<0>();
61 
62  for(unsigned int i=0; i<LFSV_PFS_V::CHILDREN; ++i)
63  {
64  scalar_jacobian_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),mat);
65  }
66  }
67 
68  private:
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,
71  R& r) const
72  {
73 
74  // Switches between local and global interface
75  typedef FiniteElementInterfaceSwitch<
76  typename LFSU::Traits::FiniteElementType
77  > FESwitch;
78  typedef BasisInterfaceSwitch<
79  typename FESwitch::Basis
80  > BasisSwitch;
81 
82  // domain and range field type
83  typedef typename BasisSwitch::DomainField DF;
84  typedef typename BasisSwitch::RangeField RF;
85  typedef typename BasisSwitch::Range RangeType;
86 
87  typedef typename LFSU::Traits::SizeType size_type;
88 
89  // dimensions
90  const int dim = EG::Geometry::dimension;
91 
92  // select quadrature rule
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);
98 
99  // loop over quadrature points
100  for (const auto& ip : rule)
101  {
102  // evaluate basis functions
103  std::vector<RangeType> phi(lfsu.size());
104  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(ip.position(),phi);
105 
106  RF rho = p.rho(eg,ip.position());
107  // evaluate u
108  RF u=0.0;
109  for (size_type i=0; i<lfsu.size(); i++)
110  u += x(lfsu,i)*phi[i];
111 
112  // u*phi_i
113  RF factor = ip.weight() * rho * eg.geometry().integrationElement(ip.position());
114 
115  for (size_type i=0; i<lfsu.size(); i++)
116  r.accumulate(lfsv,i, u*phi[i]*factor);
117  }
118  }
119 
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,
122  M& mat) const
123  {
124 
125  // Switches between local and global interface
126  typedef FiniteElementInterfaceSwitch<
127  typename LFSU::Traits::FiniteElementType
128  > FESwitch;
129  typedef BasisInterfaceSwitch<
130  typename FESwitch::Basis
131  > BasisSwitch;
132 
133  // domain and range field type
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;
138 
139  // dimensions
140  const int dim = EG::Geometry::dimension;
141 
142  // select quadrature rule
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);
148 
149  // loop over quadrature points
150  for (const auto& ip : rule)
151  {
152  // evaluate basis functions
153  std::vector<RangeType> phi(lfsu.size());
154  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(ip.position(),phi);
155 
156  // integrate phi_j*phi_i
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);
162  }
163  }
164 
165  const PRM& p;
166  const int superintegration_order;
167  }; // end class NavierStokesMass
168 
177  template<typename PRM>
179  public FullVolumePattern ,
182  {
183  public:
184  // pattern assembly flags
185  enum { doPatternVolume = true };
186 
187  // residual assembly flags
188  enum { doAlphaVolume = true };
189 
190  NavierStokesVelVecMass (const PRM & p_, int superintegration_order_ = 0)
191  : p(p_), superintegration_order(superintegration_order_)
192  {}
193 
194  // volume integral depending on test and ansatz functions
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
197  {
198  // dimensions
199  const int dim = EG::Geometry::dimension;
200 
201  // subspaces
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>();
205 
206  // domain and range field type
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;
213 
214  // select quadrature rule
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);
220 
221  // loop over quadrature points
222  for (const auto& ip : rule)
223  {
224  const Dune::FieldVector<DF,dim> local = ip.position();
225  const RF rho = p.rho(eg,local);
226 
227  // compute basis functions
228  std::vector<Range_V> phi_v(lfsv_v.size());
229  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
230 
231  // compute u
232  Range_V u(0.0);
233  for(size_type i=0; i<lfsu_v.size(); i++)
234  u.axpy(x(lfsu_v,i),phi_v[i]);
235 
236  const RF factor = ip.weight() * rho * eg.geometry().integrationElement(ip.position());
237 
238  for(size_type i=0; i<lfsv_v.size(); i++)
239  r.accumulate(lfsv_v,i, (u*phi_v[i]) * factor);
240 
241  } // end loop quadrature points
242  } // end alpha_volume
243 
244  // jacobian of volume term
245  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
246  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
247  M& mat) const
248  {
249  // dimensions
250  const int dim = EG::Geometry::dimension;
251 
252  // subspaces
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>();
256 
257  // domain and range field type
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;
264 
265  // select quadrature rule
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);
271 
272  // loop over quadrature points
273  for (const auto& ip : rule)
274  {
275  const Dune::FieldVector<DF,dim> local = ip.position();
276  const RF rho = p.rho(eg,local);
277 
278  // compute basis functions
279  std::vector<Range_V> phi_v(lfsv_v.size());
280  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
281 
282  const RF factor = ip.weight() * rho * eg.geometry().integrationElement(ip.position());
283 
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);
287  } // end loop quadrature points
288  } // end jacobian_volume
289 
290  private :
291  const PRM& p;
292  const int superintegration_order;
293  }; // end class NavierStokesVelVecMass
294 
295  } // end namespace PDELab
296 } // end namespace Dune
297 #endif
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
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