DUNE PDELab (git)

navierstokesmass.hh
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
13#include <dune/pdelab/localoperator/pattern.hh>
14#include <dune/pdelab/localoperator/flags.hh>
15#include <dune/pdelab/localoperator/idefault.hh>
16
17#include <dune/pdelab/common/quadraturerules.hh>
18
19namespace 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
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
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
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
typename impl::_Child< Node, indices... >::type Child
Template alias for the type of a child node given by a list of child indices.
Definition: childextraction.hh:225
ImplementationDefined child(Node &&node, Indices... indices)
Extracts the child of a node given by a sequence of compile-time and run-time indices.
Definition: childextraction.hh:128
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)