DUNE PDELab (git)

taylorhoodnavierstokes.hh
1// -*- tab-width: 2; indent-tabs-mode: nil -*-
2#ifndef DUNE_PDELAB_LOCALOPERATOR_TAYLORHOODNAVIERSTOKES_HH
3#define DUNE_PDELAB_LOCALOPERATOR_TAYLORHOODNAVIERSTOKES_HH
4
5#include<vector>
6
9
11#include<dune/geometry/referenceelements.hh>
12
13#include<dune/pdelab/common/quadraturerules.hh>
15
16#include"defaultimp.hh"
17#include"pattern.hh"
18#include"idefault.hh"
19#include"flags.hh"
20#include"l2.hh"
21#include"stokesparameter.hh"
22#include"navierstokesmass.hh"
23
24namespace Dune {
25 namespace PDELab {
26
30
51 template<typename P>
53 public FullVolumePattern,
55 public InstationaryLocalOperatorDefaultMethods<typename P::Traits::RangeField>
56 {
57 public:
60
62 using Real = typename InstatBase::RealType;
63
64 static const bool navier = P::assemble_navier;
65 static const bool full_tensor = P::assemble_full_tensor;
66
67 // pattern assembly flags
68 enum { doPatternVolume = true };
69
70 // residual assembly flags
71 enum { doAlphaVolume = true };
72 enum { doLambdaVolume = true };
73 enum { doLambdaBoundary = true };
74
75 using PhysicalParameters = P;
76
77 TaylorHoodNavierStokes (PhysicalParameters& p, int superintegration_order_ = 0)
78
79 : _p(p)
80 , superintegration_order(superintegration_order_)
81 {}
82
83 // set time in parameter class
84 void setTime(Real t)
85 {
87 _p.setTime(t);
88 }
89
90 // volume integral depending on test and ansatz functions
91 template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
92 void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
93 {
94 // define types
95 using namespace Indices;
96 using LFSU_V_PFS = TypeTree::Child<LFSU,_0>;
97 using LFSU_V = TypeTree::Child<LFSU_V_PFS,_0>;
98 using LFSU_P = TypeTree::Child<LFSU,_1>;
99 using RF = typename LFSU_V::Traits::FiniteElementType::
100 Traits::LocalBasisType::Traits::RangeFieldType;
101 using RT_V = typename LFSU_V::Traits::FiniteElementType::
102 Traits::LocalBasisType::Traits::RangeType;
103 using JacobianType_V = typename LFSU_V::Traits::FiniteElementType::
104 Traits::LocalBasisType::Traits::JacobianType;
105 using RT_P = typename LFSU_P::Traits::FiniteElementType::
106 Traits::LocalBasisType::Traits::RangeType;
107
108 // extract local function spaces
109 const auto& lfsu_v_pfs = child(lfsu,_0);
110 const unsigned int vsize = lfsu_v_pfs.child(0).size();
111 const auto& lfsu_p = child(lfsu,_1);
112 const unsigned int psize = lfsu_p.size();
113
114 // dimensions
115 const int dim = EG::Geometry::mydimension;
116
117 // get geometry
118 auto geo = eg.geometry();
119
120 // determine quadrature order
121 const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
122 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
123 const int jac_order = geo.type().isSimplex() ? 0 : 1;
124 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
125
126 // Initialize vectors outside for loop
127 typename EG::Geometry::JacobianInverseTransposed jac;
128 std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
129 std::vector<RT_P> psi(psize);
131 std::vector<RT_V> phi(vsize);
133
134 // loop over quadrature points
135 for (const auto& ip : quadratureRule(geo,qorder))
136 {
137 // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
138 std::vector<JacobianType_V> js(vsize);
139 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
140
141 // transform gradient to real element
142 jac = geo.jacobianInverseTransposed(ip.position());
143 for (size_t i=0; i<vsize; i++)
144 {
145 gradphi[i] = 0.0;
146 jac.umv(js[i][0],gradphi[i]);
147 }
148
149 // evaluate basis functions
150 lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
151
152 // compute u (if Navier term enabled)
153 if(navier)
154 {
155 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
156
157 for(int d=0; d<dim; ++d)
158 {
159 vu[d] = 0.0;
160 const auto& lfsu_v = lfsu_v_pfs.child(d);
161 for (size_t i=0; i<lfsu_v.size(); i++)
162 vu[d] += x(lfsu_v,i) * phi[i];
163 }
164 }
165
166 // Compute velocity jacobian
167 for(int d=0; d<dim; ++d){
168 jacu[d] = 0.0;
169 const auto& lfsu_v = lfsu_v_pfs.child(d);
170 for (size_t i=0; i<lfsu_v.size(); i++)
171 jacu[d].axpy(x(lfsu_v,i),gradphi[i]);
172 }
173
174 // compute pressure
175 RT_P func_p(0.0);
176 for (size_t i=0; i<lfsu_p.size(); i++)
177 func_p += x(lfsu_p,i) * psi[i];
178
179 // Viscosity and density
180 const auto mu = _p.mu(eg,ip.position());
181 const auto rho = _p.rho(eg,ip.position());
182
183 // geometric weight
184 const auto factor = ip.weight() * geo.integrationElement(ip.position());
185
186 for(int d=0; d<dim; ++d){
187
188 const auto& lfsu_v = lfsu_v_pfs.child(d);
189
190 //compute u * grad u_d
191 const auto u_nabla_u = vu * jacu[d];
192
193 for (size_t i=0; i<vsize; i++){
194
195 // integrate mu * grad u * grad phi_i
196 r.accumulate(lfsu_v,i, mu * (jacu[d] * gradphi[i]) * factor);
197
198 // integrate (grad u)^T * grad phi_i
199 if (full_tensor)
200 for(int dd=0; dd<dim; ++dd)
201 r.accumulate(lfsu_v,i, mu * (jacu[dd][d] * gradphi[i][dd]) * factor);
202
203 // integrate div phi_i * p
204 r.accumulate(lfsu_v,i,- (func_p * gradphi[i][d]) * factor);
205
206 // integrate u * grad u * phi_i
207 if(navier)
208 r.accumulate(lfsu_v,i, rho * u_nabla_u * phi[i] * factor);
209 }
210
211 }
212
213 // integrate div u * psi_i
214 for (size_t i=0; i<psize; i++)
215 for(int d=0; d<dim; ++d)
216 // divergence of u is the trace of the velocity jacobian
217 r.accumulate(lfsu_p,i, -1.0 * jacu[d][d] * psi[i] * factor);
218
219 }
220 }
221
222
223 // volume integral depending on test functions
224 template<typename EG, typename LFSV, typename R>
225 void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
226 {
227 // define types
228 using namespace Indices;
229 using LFSV_V_PFS = TypeTree::Child<LFSV,_0>;
230 using LFSV_V = TypeTree::Child<LFSV_V_PFS,_0>;
231 using LFSV_P = TypeTree::Child<LFSV,_1>;
232 using RT_V = typename LFSV_V::Traits::FiniteElementType::
233 Traits::LocalBasisType::Traits::RangeType;
234 using RT_P = typename LFSV_P::Traits::FiniteElementType::
235 Traits::LocalBasisType::Traits::RangeType;
236
237 // extract local function spaces
238 const auto& lfsv_v_pfs = child(lfsv,_0);
239 const unsigned int vsize = lfsv_v_pfs.child(0).size();
240 const auto& lfsv_p = child(lfsv,_1);
241 const unsigned int psize = lfsv_p.size();
242
243 // dimensions
244 const int dim = EG::Geometry::mydimension;
245
246 // get cell
247 const auto& cell = eg.entity();
248
249 // get geometry
250 auto geo = eg.geometry();
251
252 // determine quadrature order
253 const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
254 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
255 const int qorder = 2*v_order + det_jac_order + superintegration_order;
256
257 // Initialize vectors outside for loop
258 std::vector<RT_V> phi(vsize);
259 std::vector<RT_P> psi(psize);
260
261 // loop over quadrature points
262 for (const auto& ip : quadratureRule(geo,qorder))
263 {
264 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
265
266 lfsv_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
267
268 // forcing term
269 const auto f1 = _p.f(cell,ip.position());
270
271 // geometric weight
272 const auto factor = ip.weight() * geo.integrationElement(ip.position());
273
274 for(int d=0; d<dim; ++d){
275
276 const auto& lfsv_v = lfsv_v_pfs.child(d);
277
278 for (size_t i=0; i<vsize; i++)
279 {
280 // integrate f1 * phi_i
281 r.accumulate(lfsv_v,i, -f1[d]*phi[i] * factor);
282 }
283
284 }
285
286 const auto g2 = _p.g2(eg,ip.position());
287
288 // integrate div u * psi_i
289 for (size_t i=0; i<psize; i++)
290 {
291 r.accumulate(lfsv_p,i, g2 * psi[i] * factor);
292 }
293
294 }
295 }
296
297
298 // residual of boundary term
299 template<typename IG, typename LFSV, typename R>
300 void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r) const
301 {
302 // define types
303 using namespace Indices;
304 using LFSV_V_PFS = TypeTree::Child<LFSV,_0>;
305 using LFSV_V = TypeTree::Child<LFSV_V_PFS,_0>;
306 using RT_V = typename LFSV_V::Traits::FiniteElementType::
307 Traits::LocalBasisType::Traits::RangeType;
308
309 // extract local velocity function spaces
310 const auto& lfsv_v_pfs = child(lfsv,_0);
311 const unsigned int vsize = lfsv_v_pfs.child(0).size();
312
313 // dimensions
314 static const unsigned int dim = IG::Entity::dimension;
315
316 // get geometry
317 auto geo = ig.geometry();
318
319 // Get geometry of intersection in local coordinates of inside cell
320 auto geo_in_inside = ig.geometryInInside();
321
322 // determine quadrature order
323 const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
324 const int det_jac_order = geo_in_inside.type().isSimplex() ? 0 : (dim-2);
325 const int jac_order = geo_in_inside.type().isSimplex() ? 0 : 1;
326 const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
327
328 // Initialize vectors outside for loop
329 std::vector<RT_V> phi(vsize);
330
331 // loop over quadrature points and integrate normal flux
332 for (const auto& ip : quadratureRule(geo,qorder))
333 {
334 // evaluate boundary condition type
335 auto bctype = _p.bctype(ig,ip.position());
336
337 // skip rest if we are on Dirichlet boundary
338 if (bctype != BC::StressNeumann)
339 continue;
340
341 // position of quadrature point in local coordinates of element
342 auto local = geo_in_inside.global(ip.position());
343
344 // evaluate basis functions
345 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
346
347 const auto factor = ip.weight() * geo.integrationElement(ip.position());
348 const auto normal = ig.unitOuterNormal(ip.position());
349
350 // evaluate flux boundary condition
351 const auto neumann_stress = _p.j(ig,ip.position(),normal);
352
353 for(unsigned int d=0; d<dim; ++d)
354 {
355
356 const auto& lfsv_v = lfsv_v_pfs.child(d);
357
358 for (size_t i=0; i<vsize; i++)
359 {
360 r.accumulate(lfsv_v,i, neumann_stress[d] * phi[i] * factor);
361 }
362
363 }
364 }
365 }
366
367
368 template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
369 void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
370 M& mat) const
371 {
372 // define types
373 using namespace Indices;
374 using LFSU_V_PFS = TypeTree::Child<LFSU,_0>;
375 using LFSU_V = TypeTree::Child<LFSU_V_PFS,_0>;
376 using LFSU_P = TypeTree::Child<LFSU,_1>;
377 using RF = typename LFSU_V::Traits::FiniteElementType::
378 Traits::LocalBasisType::Traits::RangeFieldType;
379 using RT_V = typename LFSU_V::Traits::FiniteElementType::
380 Traits::LocalBasisType::Traits::RangeType;
381 using JacobianType_V = typename LFSU_V::Traits::FiniteElementType::
382 Traits::LocalBasisType::Traits::JacobianType;
383 using RT_P = typename LFSU_P::Traits::FiniteElementType::
384 Traits::LocalBasisType::Traits::RangeType;
385
386 // extract local function spaces
387 const auto& lfsu_v_pfs = child(lfsu,_0);
388 const unsigned int vsize = lfsu_v_pfs.child(0).size();
389 const auto& lfsu_p = child(lfsu,_1);
390 const unsigned int psize = lfsu_p.size();
391
392 // dimensions
393 const int dim = EG::Geometry::mydimension;
394
395 // get geometry
396 auto geo = eg.geometry();
397
398 // determine quadrature order
399 const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
400 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
401 const int jac_order = geo.type().isSimplex() ? 0 : 1;
402 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
403
404 // Initialize vectors outside for loop
405 typename EG::Geometry::JacobianInverseTransposed jac;
406 std::vector<JacobianType_V> js(vsize);
407 std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
408 std::vector<RT_P> psi(psize);
409 std::vector<RT_V> phi(vsize);
411 Dune::FieldVector<RF,dim> gradu_d(0.0);
412
413 // loop over quadrature points
414 for (const auto& ip : quadratureRule(geo,qorder))
415 {
416 // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
417 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
418
419 // transform gradient to real element
420 jac = geo.jacobianInverseTransposed(ip.position());
421 for (size_t i=0; i<vsize; i++)
422 {
423 gradphi[i] = 0.0;
424 jac.umv(js[i][0],gradphi[i]);
425 }
426
427 // evaluate basis functions
428 lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
429
430 // compute u (if Navier term enabled)
431 if(navier){
432 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
433
434 for(int d = 0; d < dim; ++d){
435 vu[d] = 0.0;
436 const auto& lfsv_v = lfsu_v_pfs.child(d);
437 for(size_t l = 0; l < vsize; ++l)
438 vu[d] += x(lfsv_v,l) * phi[l];
439 }
440 }
441
442 // Viscosity and density
443 const auto mu = _p.mu(eg,ip.position());
444 const auto rho = _p.rho(eg,ip.position());
445
446 const auto factor = ip.weight() * geo.integrationElement(ip.position());
447
448 for(int d=0; d<dim; ++d){
449
450 const auto& lfsv_v = lfsu_v_pfs.child(d);
451 const auto& lfsu_v = lfsv_v;
452
453 // Derivatives of d-th velocity component
454 gradu_d = 0.0;
455 if(navier)
456 for(size_t l =0; l < vsize; ++l)
457 gradu_d.axpy(x(lfsv_v,l), gradphi[l]);
458
459 for (size_t i=0; i<lfsv_v.size(); i++){
460
461 // integrate grad phi_u_i * grad phi_v_i (viscous force)
462 for (size_t j=0; j<lfsv_v.size(); j++){
463 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[i] * gradphi[j]) * factor);
464
465 // integrate (grad phi_u_i)^T * grad phi_v_i (viscous force)
466 if(full_tensor)
467 for(int dd=0; dd<dim; ++dd){
468 const auto& lfsu_v = lfsu_v_pfs.child(dd);
469 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[j][d] * gradphi[i][dd]) * factor);
470 }
471
472 }
473
474 // integrate grad_d phi_v_d * p_u (pressure force)
475 for (size_t j=0; j<psize; j++)
476 mat.accumulate(lfsv_v,i,lfsu_p,j, - (gradphi[i][d] * psi[j]) * factor);
477
478 if(navier){
479 for(int k =0; k < dim; ++k){
480 const auto& lfsu_v = lfsu_v_pfs.child(k);
481
482 const auto pre_factor = factor * rho * gradu_d[k] * phi[i];
483
484 for(size_t j=0; j< lfsu_v.size(); ++j)
485 mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * phi[j]);
486 } // k
487
488 const auto pre_factor = factor * rho * phi[i];
489 for(size_t j=0; j< lfsu_v.size(); ++j)
490 mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * (vu * gradphi[j]));
491 }
492
493 }
494
495 for (size_t i=0; i<psize; i++){
496 for (size_t j=0; j<lfsu_v.size(); j++)
497 mat.accumulate(lfsu_p,i,lfsu_v,j, - (gradphi[j][d] * psi[i]) * factor);
498 }
499 } // d
500 } // it
501 }
502
503 private:
504 P& _p;
505 const int superintegration_order;
506 };
507
509 } // namespace PDELab
510} // namespace Dune
511
512#endif // DUNE_PDELAB_LOCALOPERATOR_TAYLORHOODNAVIERSTOKES_HH
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
sparsity pattern generator
Definition: pattern.hh:14
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104
Default flags for all local operators.
Definition: flags.hh:19
A local operator for the Navier-Stokes equations.
Definition: taylorhoodnavierstokes.hh:56
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:52
constexpr index_constant< 1 > _1
Compile time index with value 1.
Definition: indices.hh:55
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
Definition: stokesparameter.hh:32
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)