DUNE PDELab (git)

dgnavierstokes.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_DGNAVIERSTOKES_HH
5#define DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKES_HH
6
10
11#include <dune/geometry/referenceelements.hh>
12
13#include <dune/localfunctions/common/interfaceswitch.hh>
14#include <dune/pdelab/localoperator/idefault.hh>
15
16#include <dune/pdelab/common/quadraturerules.hh>
17#include <dune/pdelab/localoperator/defaultimp.hh>
18#include <dune/pdelab/localoperator/pattern.hh>
19#include <dune/pdelab/localoperator/flags.hh>
20#include <dune/pdelab/localoperator/dgnavierstokesparameter.hh>
21#include <dune/pdelab/localoperator/navierstokesmass.hh>
22
23namespace Dune {
24 namespace PDELab {
25
31 template<typename PRM>
35 public InstationaryLocalOperatorDefaultMethods<typename PRM::Traits::RangeField>
36 {
38 using RF = typename PRM::Traits::RangeField;
39
41 using Real = typename InstatBase::RealType;
42
43 static const bool navier = PRM::assemble_navier;
44 static const bool full_tensor = PRM::assemble_full_tensor;
45
46 public:
47
48 // pattern assembly flags
49 enum { doPatternVolume = true };
50 enum { doPatternSkeleton = true };
51
52 // call the assembler for each face only once
53 enum { doSkeletonTwoSided = false };
54
55 // residual assembly flags
56 enum { doAlphaVolume = true };
57 enum { doAlphaSkeleton = true };
58 enum { doAlphaBoundary = true };
59 enum { doLambdaVolume = true };
60 enum { doLambdaBoundary = true };
61
74 DGNavierStokes (PRM& _prm, int _superintegration_order=0) :
75 prm(_prm), superintegration_order(_superintegration_order),
76 current_dt(1.0)
77 {}
78
79 // Store current dt
80 void preStep (Real , Real dt, int )
81 {
82 current_dt = dt;
83 }
84
85 // set time in parameter class
86 void setTime(Real t)
87 {
89 prm.setTime(t);
90 }
91
92 // volume integral depending on test and ansatz functions
93 template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
94 void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
95 {
96 // dimensions
97 const unsigned int dim = EG::Geometry::mydimension;
98
99 // subspaces
100 using namespace Indices;
101 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
102 const auto& lfsv_pfs_v = child(lfsv,_0);
103
104 // ... we assume all velocity components are the same
105 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
106 const auto& lfsv_v = child(lfsv_pfs_v,_0);
107 const unsigned int vsize = lfsv_v.size();
108 using LFSV_P = TypeTree::Child<LFSV,_1>;
109 const auto& lfsv_p = child(lfsv,_1);
110 const unsigned int psize = lfsv_p.size();
111
112 // domain and range field type
113 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
114 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
115 using RT = typename BasisSwitch_V::Range;
116 using RF = typename BasisSwitch_V::RangeField;
117 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
118 using size_type = typename LFSV::Traits::SizeType;
119
120 // Get geometry
121 auto geo = eg.geometry();
122
123 // Determine quadrature order
124 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
125 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
126 const int jac_order = geo.type().isSimplex() ? 0 : 1;
127 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
128
129 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
130
131 // Initialize vectors outside for loop
132 std::vector<RT> phi_v(vsize);
134 std::vector<RT> phi_p(psize);
135 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
137
138 // loop over quadrature points
139 for (const auto& ip : quadratureRule(geo,qorder))
140 {
141 auto local = ip.position();
142 auto mu = prm.mu(eg,local);
143 auto rho = prm.rho(eg,local);
144
145 // compute u (if Navier term enabled)
146 if(navier) {
147 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
148
149 for(unsigned int d=0; d<dim; ++d) {
150 vu[d] = 0.0;
151 const auto& lfsu_v = lfsv_pfs_v.child(d);
152 for(size_type i=0; i<vsize; i++)
153 vu[d] += x(lfsu_v,i) * phi_v[i];
154 }
155 } // end navier
156
157 // and value of pressure shape functions
158 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
159
160 // compute pressure
161 RF p(0.0);
162 for(size_type i=0; i<psize; i++)
163 p += x(lfsv_p,i) * phi_p[i];
164
165 // compute gradients
166 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
167 geo, local, grad_phi_v);
168
169 // compute velocity jacobian
170 for(unsigned int d = 0; d<dim; ++d) {
171 jacu[d] = 0.0;
172 const auto& lfsv_v = lfsv_pfs_v.child(d);
173 for(size_type i=0; i<vsize; i++)
174 jacu[d].axpy(x(lfsv_v,i), grad_phi_v[i][0]);
175 }
176
177 auto detj = geo.integrationElement(ip.position());
178 auto weight = ip.weight() * detj;
179
180 for(unsigned int d = 0; d<dim; ++d) {
181 const auto& lfsv_v = lfsv_pfs_v.child(d);
182
183 for(size_type i=0; i<vsize; i++) {
184 //================================================//
185 // \int (mu*grad_u*grad_v)
186 //================================================//
187 r.accumulate(lfsv_v,i, mu * (jacu[d]*grad_phi_v[i][0]) * weight);
188
189 // Assemble symmetric part for (grad u)^T
190 if(full_tensor)
191 for(unsigned int dd = 0; dd<dim; ++dd)
192 r.accumulate(lfsv_v,i, mu * jacu[dd][d] * grad_phi_v[i][0][dd] * weight);
193
194 //================================================//
195 // \int -p \nabla\cdot v
196 //================================================//
197 r.accumulate(lfsv_v,i, -p * grad_phi_v[i][0][d] * weight);
198
199 //================================================//
200 // \int \rho ((u\cdot\nabla ) u )\cdot v
201 //================================================//
202 if(navier) {
203 // compute u * grad u_d
204 auto u_nabla_u = vu * jacu[d];
205
206 r.accumulate(lfsv_v,i, rho * u_nabla_u * phi_v[i] * weight);
207 } // end navier
208
209 } // end i
210 } // end d
211
212 //================================================//
213 // \int -q \nabla\cdot u
214 //================================================//
215 for(size_type i=0; i<psize; i++)
216 for(unsigned int d = 0; d < dim; ++d)
217 // divergence of u is the trace of the velocity jacobian
218 r.accumulate(lfsv_p,i, -jacu[d][d] * phi_p[i] * incomp_scaling * weight);
219
220 } // end loop quadrature points
221 } // end alpha_volume
222
223 // jacobian of volume term
224 template<typename EG, typename LFSU, typename X, typename LFSV,
225 typename LocalMatrix>
226 void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
227 LocalMatrix& mat) const
228 {
229 // dimensions
230 const unsigned int dim = EG::Geometry::mydimension;
231
232 // subspaces
233 using namespace Indices;
234 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
235 const auto& lfsv_pfs_v = child(lfsv,_0);
236
237 // ... we assume all velocity components are the same
238 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
239 const auto& lfsv_v = child(lfsv_pfs_v,_0);
240 const unsigned int vsize = lfsv_v.size();
241 using LFSV_P = TypeTree::Child<LFSV,_1>;
242 const auto& lfsv_p = child(lfsv,_1);
243 const unsigned int psize = lfsv_p.size();
244
245 // domain and range field type
246 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
247 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
248 using RT = typename BasisSwitch_V::Range;
249 using RF = typename BasisSwitch_V::RangeField;
250 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
251 using size_type = typename LFSV::Traits::SizeType;
252
253 // Get geometry
254 auto geo = eg.geometry();
255
256 // Determine quadrature order
257 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
258 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
259 const int jac_order = geo.type().isSimplex() ? 0 : 1;
260 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
261
262 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
263
264 // Initialize vectors outside for loop
265 std::vector<RT> phi_v(vsize);
267 std::vector<RT> phi_p(psize);
268 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
269 Dune::FieldVector<RF,dim> gradu_dv(0.0);
270
271 // loop over quadrature points
272 for (const auto& ip : quadratureRule(geo,qorder))
273 {
274 auto local = ip.position();
275 auto mu = prm.mu(eg,local);
276 auto rho = prm.rho(eg,local);
277
278 // compute u (if Navier term enabled)
279 if(navier) {
280 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
281
282 for(unsigned int d=0; d<dim; ++d) {
283 vu[d] = 0.0;
284 const auto& lfsu_v = lfsv_pfs_v.child(d);
285 for(size_type i=0; i<vsize; i++)
286 vu[d] += x(lfsu_v,i) * phi_v[i];
287 }
288 } // end navier
289
290 // and value of pressure shape functions
291 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
292
293 // compute gradients
294 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
295 geo, local, grad_phi_v);
296
297 auto detj = geo.integrationElement(ip.position());
298 auto weight = ip.weight() * detj;
299
300 for(unsigned int dv = 0; dv<dim; ++dv) {
301 const LFSV_V& lfsv_v = lfsv_pfs_v.child(dv);
302
303 // gradient of dv-th velocity component
304 gradu_dv = 0.0;
305 if(navier)
306 for(size_type l=0; l<vsize; ++l)
307 gradu_dv.axpy(x(lfsv_v,l), grad_phi_v[l][0]);
308
309 for(size_type i=0; i<vsize; i++) {
310
311 for(size_type j=0; j<vsize; j++) {
312 //================================================//
313 // \int (mu*grad_u*grad_v)
314 //================================================//
315 mat.accumulate(lfsv_v,i,lfsv_v,j, mu * (grad_phi_v[j][0]*grad_phi_v[i][0]) * weight);
316
317 // Assemble symmetric part for (grad u)^T
318 if(full_tensor)
319 for(unsigned int du = 0; du<dim; ++du) {
320 const auto& lfsu_v = lfsv_pfs_v.child(du);
321 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (grad_phi_v[j][0][dv]*grad_phi_v[i][0][du]) * weight);
322 }
323 }
324
325 //================================================//
326 // - q * div u
327 // - p * div v
328 //================================================//
329 for(size_type j=0; j<psize; j++) {
330 mat.accumulate(lfsv_p,j,lfsv_v,i, -phi_p[j] * grad_phi_v[i][0][dv] * incomp_scaling * weight);
331 mat.accumulate(lfsv_v,i,lfsv_p,j, -phi_p[j] * grad_phi_v[i][0][dv] * weight);
332 }
333
334 //================================================//
335 // \int \rho ((u\cdot\nabla ) u )\cdot v
336 //================================================//
337 if(navier) {
338
339 // block diagonal contribution
340 for(size_type j=0; j<vsize; j++)
341 mat.accumulate(lfsv_v,i,lfsv_v,j, rho * (vu * grad_phi_v[j][0]) * phi_v[i] * weight);
342
343 // remaining contribution
344 for(unsigned int du = 0; du < dim; ++du) {
345 const auto& lfsu_v = lfsv_pfs_v.child(du);
346 for(size_type j=0; j<vsize; j++)
347 mat.accumulate(lfsv_v,i,lfsu_v,j, rho * phi_v[j] * gradu_dv[du] * phi_v[i] * weight);
348 }
349
350 } // end navier
351
352 } // end i
353 } // end dv
354
355 } // end loop quadrature points
356 } // end jacobian_volume
357
358 // skeleton integral depending on test and ansatz functions
359 // each face is only visited ONCE!
360 template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
361 void alpha_skeleton (const IG& ig,
362 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
363 const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
364 R& r_s, R& r_n) const
365 {
366 // dimensions
367 const unsigned int dim = IG::Entity::dimension;
368
369 // subspaces
370 using namespace Indices;
371 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
372 const auto& lfsv_s_pfs_v = child(lfsv_s,_0);
373 const auto& lfsv_n_pfs_v = child(lfsv_n,_0);
374
375 // ... we assume all velocity components are the same
376 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
377 const auto& lfsv_s_v = child(lfsv_s_pfs_v,_0);
378 const auto& lfsv_n_v = child(lfsv_n_pfs_v,_0);
379 const unsigned int vsize_s = lfsv_s_v.size();
380 const unsigned int vsize_n = lfsv_n_v.size();
381 using LFSV_P = TypeTree::Child<LFSV,_1>;
382 const auto& lfsv_s_p = child(lfsv_s,_1);
383 const auto& lfsv_n_p = child(lfsv_n,_1);
384 const unsigned int psize_s = lfsv_s_p.size();
385 const unsigned int psize_n = lfsv_n_p.size();
386
387 // domain and range field type
388 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
389 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
390 using RT = typename BasisSwitch_V::Range;
391 using RF = typename BasisSwitch_V::RangeField;
392 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
393
394 // References to inside and outside cells
395 const auto& cell_inside = ig.inside();
396 const auto& cell_outside = ig.outside();
397
398 // Get geometries
399 auto geo = ig.geometry();
400 auto geo_inside = cell_inside.geometry();
401 auto geo_outside = cell_outside.geometry();
402
403 // Get geometry of intersection in local coordinates of cell_inside and cell_outside
404 auto geo_in_inside = ig.geometryInInside();
405 auto geo_in_outside = ig.geometryInOutside();
406
407 // Determine quadrature order
408 const int v_order = FESwitch_V::basis(lfsv_s_v.finiteElement()).order();
409 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
410 const int qorder = 2*v_order + det_jac_order + superintegration_order;
411
412 const int epsilon = prm.epsilonIPSymmetryFactor();
413 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
414
415 // Initialize vectors outside for loop
416 std::vector<RT> phi_v_s(vsize_s);
417 std::vector<RT> phi_v_n(vsize_n);
418 Dune::FieldVector<RF,dim> u_s(0.0), u_n(0.0);
419 std::vector<RT> phi_p_s(psize_s);
420 std::vector<RT> phi_p_n(psize_n);
421 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_s(vsize_s);
422 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_n(vsize_n);
423 Dune::FieldMatrix<RF,dim,dim> jacu_s(0.0), jacu_n(0.0);
424
425 auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
426
427 // loop over quadrature points and integrate normal flux
428 for (const auto& ip : quadratureRule(geo,qorder))
429 {
430
431 // position of quadrature point in local coordinates of element
432 auto local_s = geo_in_inside.global(ip.position());
433 auto local_n = geo_in_outside.global(ip.position());
434
435 // value of velocity shape functions
436 FESwitch_V::basis(lfsv_s_v.finiteElement()).evaluateFunction(local_s,phi_v_s);
437 FESwitch_V::basis(lfsv_n_v.finiteElement()).evaluateFunction(local_n,phi_v_n);
438
439 // evaluate u
440 assert(vsize_s == vsize_n);
441 for(unsigned int d=0; d<dim; ++d) {
442 u_s[d] = 0.0;
443 u_n[d] = 0.0;
444 const auto& lfsv_s_v = lfsv_s_pfs_v.child(d);
445 const auto& lfsv_n_v = lfsv_n_pfs_v.child(d);
446 for(unsigned int i=0; i<vsize_s; i++) {
447 u_s[d] += x_s(lfsv_s_v,i) * phi_v_s[i];
448 u_n[d] += x_n(lfsv_n_v,i) * phi_v_n[i];
449 }
450 }
451
452 // value of pressure shape functions
453 FESwitch_P::basis(lfsv_s_p.finiteElement()).evaluateFunction(local_s,phi_p_s);
454 FESwitch_P::basis(lfsv_n_p.finiteElement()).evaluateFunction(local_n,phi_p_n);
455
456 // evaluate pressure
457 assert(psize_s == psize_n);
458 RF p_s(0.0), p_n(0.0);
459 for(unsigned int i=0; i<psize_s; i++) {
460 p_s += x_s(lfsv_s_p,i) * phi_p_s[i];
461 p_n += x_n(lfsv_n_p,i) * phi_p_n[i];
462 }
463
464 // compute gradients
465 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_s_v.finiteElement()),
466 geo_inside, local_s, grad_phi_v_s);
467
468 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_n_v.finiteElement()),
469 geo_outside, local_n, grad_phi_v_n);
470
471 // evaluate velocity jacobian
472 for(unsigned int d=0; d<dim; ++d) {
473 jacu_s[d] = 0.0;
474 jacu_n[d] = 0.0;
475 const auto& lfsv_s_v = lfsv_s_pfs_v.child(d);
476 const auto& lfsv_n_v = lfsv_n_pfs_v.child(d);
477 for(unsigned int i=0; i<vsize_s; i++) {
478 jacu_s[d].axpy(x_s(lfsv_s_v,i), grad_phi_v_s[i][0]);
479 jacu_n[d].axpy(x_n(lfsv_n_v,i), grad_phi_v_n[i][0]);
480 }
481 }
482
483 auto normal = ig.unitOuterNormal(ip.position());
484 auto weight = ip.weight()*geo.integrationElement(ip.position());
485 auto mu = prm.mu(ig,ip.position());
486
487 auto factor = mu * weight;
488
489 for(unsigned int d=0; d<dim; ++d) {
490 const auto& lfsv_s_v = lfsv_s_pfs_v.child(d);
491 const auto& lfsv_n_v = lfsv_n_pfs_v.child(d);
492
493 //================================================//
494 // diffusion term
495 //================================================//
496 auto val = 0.5 * ((jacu_s[d] * normal) + (jacu_n[d] * normal)) * factor;
497 for(unsigned int i=0; i<vsize_s; i++) {
498 r_s.accumulate(lfsv_s_v,i, -val * phi_v_s[i]);
499 r_n.accumulate(lfsv_n_v,i, val * phi_v_n[i]);
500
501 if(full_tensor) {
502 for(unsigned int dd=0; dd<dim; ++dd) {
503 auto Tval = 0.5 * (jacu_s[dd][d] + jacu_n[dd][d]) * normal[dd] * factor;
504 r_s.accumulate(lfsv_s_v,i, -Tval * phi_v_s[i]);
505 r_n.accumulate(lfsv_n_v,i, Tval * phi_v_n[i]);
506 }
507 }
508 } // end i
509
510 //================================================//
511 // (non-)symmetric IP term
512 //================================================//
513 auto jumpu_d = u_s[d] - u_n[d];
514 for(unsigned int i=0; i<vsize_s; i++) {
515 r_s.accumulate(lfsv_s_v,i, epsilon * 0.5 * (grad_phi_v_s[i][0] * normal) * jumpu_d * factor);
516 r_n.accumulate(lfsv_n_v,i, epsilon * 0.5 * (grad_phi_v_n[i][0] * normal) * jumpu_d * factor);
517
518 if(full_tensor) {
519 for(unsigned int dd=0; dd<dim; ++dd) {
520 r_s.accumulate(lfsv_s_v,i, epsilon * 0.5 * grad_phi_v_s[i][0][dd] * (u_s[dd] - u_n[dd]) * normal[d] * factor);
521 r_n.accumulate(lfsv_n_v,i, epsilon * 0.5 * grad_phi_v_n[i][0][dd] * (u_s[dd] - u_n[dd]) * normal[d] * factor);
522 }
523 }
524 } // end i
525
526 //================================================//
527 // standard IP term integral
528 //================================================//
529 for(unsigned int i=0; i<vsize_s; i++) {
530 r_s.accumulate(lfsv_s_v,i, penalty_factor * jumpu_d * phi_v_s[i] * factor);
531 r_n.accumulate(lfsv_n_v,i, -penalty_factor * jumpu_d * phi_v_n[i] * factor);
532 } // end i
533
534 //================================================//
535 // pressure-velocity-coupling in momentum equation
536 //================================================//
537 auto mean_p = 0.5*(p_s + p_n);
538 for(unsigned int i=0; i<vsize_s; i++) {
539 r_s.accumulate(lfsv_s_v,i, mean_p * phi_v_s[i] * normal[d] * weight);
540 r_n.accumulate(lfsv_n_v,i, -mean_p * phi_v_n[i] * normal[d] * weight);
541 } // end i
542 } // end d
543
544 //================================================//
545 // incompressibility constraint
546 //================================================//
547 auto jumpu_n = (u_s*normal) - (u_n*normal);
548 for(unsigned int i=0; i<psize_s; i++) {
549 r_s.accumulate(lfsv_s_p,i, 0.5 * phi_p_s[i] * jumpu_n * incomp_scaling * weight);
550 r_n.accumulate(lfsv_n_p,i, 0.5 * phi_p_n[i] * jumpu_n * incomp_scaling * weight);
551 } // end i
552
553 } // end loop quadrature points
554 } // end alpha_skeleton
555
556 // jacobian of skeleton term
557 template<typename IG, typename LFSU, typename X, typename LFSV,
558 typename LocalMatrix>
559 void jacobian_skeleton (const IG& ig,
560 const LFSU& lfsu_s, const X&, const LFSV& lfsv_s,
561 const LFSU& lfsu_n, const X&, const LFSV& lfsv_n,
562 LocalMatrix& mat_ss, LocalMatrix& mat_sn,
563 LocalMatrix& mat_ns, LocalMatrix& mat_nn) const
564 {
565 // dimensions
566 const unsigned int dim = IG::Entity::dimension;
567
568 // subspaces
569 using namespace Indices;
570 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
571 const auto& lfsv_s_pfs_v = child(lfsv_s,_0);
572 const auto& lfsv_n_pfs_v = child(lfsv_n,_0);
573
574 // ... we assume all velocity components are the same
575 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
576 const auto& lfsv_s_v = child(lfsv_s_pfs_v,_0);
577 const auto& lfsv_n_v = child(lfsv_n_pfs_v,_0);
578 const unsigned int vsize_s = lfsv_s_v.size();
579 const unsigned int vsize_n = lfsv_n_v.size();
580 using LFSV_P = TypeTree::Child<LFSV,_1>;
581 const auto& lfsv_s_p = child(lfsv_s,_1);
582 const auto& lfsv_n_p = child(lfsv_n,_1);
583 const unsigned int psize_s = lfsv_s_p.size();
584 const unsigned int psize_n = lfsv_n_p.size();
585
586 // domain and range field type
587 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
588 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
589 using RT = typename BasisSwitch_V::Range;
590 using RF = typename BasisSwitch_V::RangeField;
591 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
592
593 // References to inside and outside cells
594 auto const& cell_inside = ig.inside();
595 auto const& cell_outside = ig.outside();
596
597 // Get geometries
598 auto geo = ig.geometry();
599 auto geo_inside = cell_inside.geometry();
600 auto geo_outside = cell_outside.geometry();
601
602 // Get geometry of intersection in local coordinates of cell_inside and cell_outside
603 auto geo_in_inside = ig.geometryInInside();
604 auto geo_in_outside = ig.geometryInOutside();
605
606 // Determine quadrature order
607 const int v_order = FESwitch_V::basis(lfsv_s_v.finiteElement()).order();
608 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
609 const int qorder = 2*v_order + det_jac_order + superintegration_order;
610
611 const int epsilon = prm.epsilonIPSymmetryFactor();
612 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
613
614 // Initialize vectors outside for loop
615 std::vector<RT> phi_v_s(vsize_s);
616 std::vector<RT> phi_v_n(vsize_n);
617 std::vector<RT> phi_p_s(psize_s);
618 std::vector<RT> phi_p_n(psize_n);
619 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_s(vsize_s);
620 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_n(vsize_n);
621
622 auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
623
624 // loop over quadrature points and integrate normal flux
625 for (const auto& ip : quadratureRule(geo,qorder))
626 {
627
628 // position of quadrature point in local coordinates of element
629 auto local_s = geo_in_inside.global(ip.position());
630 auto local_n = geo_in_outside.global(ip.position());
631
632 // value of velocity shape functions
633 FESwitch_V::basis(lfsv_s_v.finiteElement()).evaluateFunction(local_s,phi_v_s);
634 FESwitch_V::basis(lfsv_n_v.finiteElement()).evaluateFunction(local_n,phi_v_n);
635 // and value of pressure shape functions
636 FESwitch_P::basis(lfsv_s_p.finiteElement()).evaluateFunction(local_s,phi_p_s);
637 FESwitch_P::basis(lfsv_n_p.finiteElement()).evaluateFunction(local_n,phi_p_n);
638
639 // compute gradients
640 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_s_v.finiteElement()),
641 geo_inside, local_s, grad_phi_v_s);
642
643 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_n_v.finiteElement()),
644 geo_outside, local_n, grad_phi_v_n);
645
646 auto normal = ig.unitOuterNormal(ip.position());
647 auto weight = ip.weight()*geo.integrationElement(ip.position());
648 auto mu = prm.mu(ig,ip.position());
649
650 assert(vsize_s == vsize_n);
651 auto factor = mu * weight;
652
653 for(unsigned int d = 0; d < dim; ++d) {
654 const auto& lfsv_s_v = lfsv_s_pfs_v.child(d);
655 const auto& lfsv_n_v = lfsv_n_pfs_v.child(d);
656
657 //================================================//
658 // - (\mu \int < \nabla u > . normal . [v])
659 // \mu \int \frac{\sigma}{|\gamma|^\beta} [u] \cdot [v]
660 //================================================//
661 for(unsigned int i=0; i<vsize_s; ++i) {
662
663 for(unsigned int j=0; j<vsize_s; ++j) {
664 auto val = (0.5*(grad_phi_v_s[i][0]*normal)*phi_v_s[j]) * factor;
665 mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v,i, -val);
666 mat_ss.accumulate(lfsv_s_v,i,lfsv_s_v,j, epsilon * val);
667 mat_ss.accumulate(lfsv_s_v,i,lfsv_s_v,j, phi_v_s[i] * phi_v_s[j] * penalty_factor * factor);
668
669 // Assemble symmetric part for (grad u)^T
670 if(full_tensor) {
671 for(unsigned int dd = 0; dd < dim; ++dd) {
672 auto Tval = (0.5*(grad_phi_v_s[i][0][d]*normal[dd])*phi_v_s[j]) * factor;
673 const auto& lfsv_s_v_dd = lfsv_s_pfs_v.child(dd);
674 mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v_dd,i, - Tval);
675 mat_ss.accumulate(lfsv_s_v_dd,i,lfsv_s_v,j, epsilon*Tval );
676 }
677 }
678 }
679
680 for(unsigned int j=0; j<vsize_n; ++j) {
681 // the normal vector flipped, thus the sign flips
682 auto val = (-0.5*(grad_phi_v_s[i][0]*normal)*phi_v_n[j]) * factor;
683 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v,i,- val);
684 mat_sn.accumulate(lfsv_s_v,i,lfsv_n_v,j, epsilon*val);
685 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v,i, -phi_v_s[i] * phi_v_n[j] * penalty_factor * factor);
686
687 // Assemble symmetric part for (grad u)^T
688 if(full_tensor) {
689 for (unsigned int dd=0;dd<dim;++dd) {
690 auto Tval = (-0.5*(grad_phi_v_s[i][0][d]*normal[dd])*phi_v_n[j]) * factor;
691 const auto& lfsv_s_v_dd = lfsv_s_pfs_v.child(dd);
692 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v_dd,i,- Tval);
693 mat_sn.accumulate(lfsv_s_v_dd,i,lfsv_n_v,j, epsilon*Tval);
694 }
695 }
696 }
697
698 for(unsigned int j=0; j<vsize_s; ++j) {
699 auto val = (0.5*(grad_phi_v_n[i][0]*normal)*phi_v_s[j]) * factor;
700 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v,i, - val);
701 mat_ns.accumulate(lfsv_n_v,i,lfsv_s_v,j, epsilon*val );
702 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v,i, -phi_v_n[i] * phi_v_s[j] * penalty_factor * factor);
703
704 // Assemble symmetric part for (grad u)^T
705 if(full_tensor) {
706 for (unsigned int dd=0;dd<dim;++dd) {
707 auto Tval = (0.5*(grad_phi_v_n[i][0][d]*normal[dd])*phi_v_s[j]) * factor;
708 const auto& lfsv_n_v_dd = lfsv_n_pfs_v.child(dd);
709 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v_dd,i, - Tval);
710 mat_ns.accumulate(lfsv_n_v_dd,i,lfsv_s_v,j, epsilon*Tval );
711 }
712 }
713 }
714
715 for(unsigned int j=0; j<vsize_n; ++j) {
716 // the normal vector flipped, thus the sign flips
717 auto val = (-0.5*(grad_phi_v_n[i][0]*normal)*phi_v_n[j]) * factor;
718 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v,i, - val);
719 mat_nn.accumulate(lfsv_n_v,i,lfsv_n_v,j, epsilon*val);
720 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v,i, phi_v_n[i] * phi_v_n[j] * penalty_factor * factor);
721
722 // Assemble symmetric part for (grad u)^T
723 if(full_tensor) {
724 for (unsigned int dd=0;dd<dim;++dd) {
725 auto Tval = (-0.5*(grad_phi_v_n[i][0][d]*normal[dd])*phi_v_n[j]) * factor;
726 const auto& lfsv_n_v_dd = lfsv_n_pfs_v.child(dd);
727 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v_dd,i,- Tval);
728 mat_nn.accumulate(lfsv_n_v_dd,i,lfsv_n_v,j, epsilon*Tval);
729 }
730 }
731 }
732
733 //================================================//
734 // \int <q> [u] n
735 // \int <p> [v] n
736 //================================================//
737 for(unsigned int j=0; j<psize_s; ++j) {
738 auto val = 0.5*(phi_p_s[j]*normal[d]*phi_v_s[i]) * weight;
739 mat_ss.accumulate(lfsv_s_v,i,lfsv_s_p,j, val);
740 mat_ss.accumulate(lfsv_s_p,j,lfsv_s_v,i, val * incomp_scaling);
741 }
742
743 for(unsigned int j=0; j<psize_n; ++j) {
744 auto val = 0.5*(phi_p_n[j]*normal[d]*phi_v_s[i]) * weight;
745 mat_sn.accumulate(lfsv_s_v,i,lfsv_n_p,j, val);
746 mat_ns.accumulate(lfsv_n_p,j,lfsv_s_v,i, val * incomp_scaling);
747 }
748
749 for (unsigned int j=0; j<psize_s;++j) {
750 // the normal vector flipped, thus the sign flips
751 auto val = -0.5*(phi_p_s[j]*normal[d]*phi_v_n[i]) * weight;
752 mat_ns.accumulate(lfsv_n_v,i,lfsv_s_p,j, val);
753 mat_sn.accumulate(lfsv_s_p,j,lfsv_n_v,i, val * incomp_scaling);
754 }
755
756 for (unsigned int j=0; j<psize_n;++j) {
757 // the normal vector flipped, thus the sign flips
758 auto val = -0.5*(phi_p_n[j]*normal[d]*phi_v_n[i]) * weight;
759 mat_nn.accumulate(lfsv_n_v,i,lfsv_n_p,j, val);
760 mat_nn.accumulate(lfsv_n_p,j,lfsv_n_v,i, val * incomp_scaling);
761 }
762 } // end i
763 } // end d
764
765 } // end loop quadrature points
766 } // end jacobian_skeleton
767
768 // boundary term
769 template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
770 void alpha_boundary (const IG& ig,
771 const LFSU& lfsu, const X& x, const LFSV& lfsv,
772 R& r) const
773 {
774 // dimensions
775 const unsigned int dim = IG::Entity::dimension;
776
777 // subspaces
778 using namespace Indices;
779 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
780 const auto& lfsv_pfs_v = child(lfsv,_0);
781
782 // ... we assume all velocity components are the same
783 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
784 const auto& lfsv_v = child(lfsv_pfs_v,_0);
785 const unsigned int vsize = lfsv_v.size();
786 using LFSV_P = TypeTree::Child<LFSV,_1>;
787 const auto& lfsv_p = child(lfsv,_1);
788 const unsigned int psize = lfsv_p.size();
789
790 // domain and range field type
791 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
792 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
793 using RT = typename BasisSwitch_V::Range;
794 using RF = typename BasisSwitch_V::RangeField;
795 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
796
797 // References to inside cell
798 const auto& cell_inside = ig.inside();
799
800 // Get geometries
801 auto geo = ig.geometry();
802 auto geo_inside = cell_inside.geometry();
803
804 // Get geometry of intersection in local coordinates of cell_inside
805 auto geo_in_inside = ig.geometryInInside();
806
807 // Determine quadrature order
808 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
809 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
810 const int qorder = 2*v_order + det_jac_order + superintegration_order;
811
812 auto epsilon = prm.epsilonIPSymmetryFactor();
813 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
814
815 // Initialize vectors outside for loop
816 std::vector<RT> phi_v(vsize);
818 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
820
821 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
822
823 // loop over quadrature points and integrate normal flux
824 for (const auto& ip : quadratureRule(geo,qorder))
825 {
826 // position of quadrature point in local coordinates of element
827 auto local = geo_in_inside.global(ip.position());
828
829 // value of velocity shape functions
830 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
831
832 // evaluate u
833 for(unsigned int d=0; d<dim; ++d) {
834 u[d] = 0.0;
835 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
836 for(unsigned int i=0; i<vsize; i++)
837 u[d] += x(lfsv_v,i) * phi_v[i];
838 }
839
840 // value of pressure shape functions
841 std::vector<RT> phi_p(psize);
842 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
843
844 // evaluate pressure
845 RF p(0.0);
846 for(unsigned int i=0; i<psize; i++)
847 p += x(lfsv_p,i) * phi_p[i];
848
849 // compute gradients
850 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
851 geo_inside, local, grad_phi_v);
852
853 // evaluate velocity jacobian
854 for(unsigned int d=0; d<dim; ++d) {
855 jacu[d] = 0.0;
856 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
857 for(unsigned int i=0; i<vsize; i++)
858 jacu[d].axpy(x(lfsv_v,i), grad_phi_v[i][0]);
859 }
860
861 auto normal = ig.unitOuterNormal(ip.position());
862 auto weight = ip.weight()*geo.integrationElement(ip.position());
863 auto mu = prm.mu(ig,ip.position());
864
865 // evaluate boundary condition type
866 auto bctype(prm.bctype(ig,ip.position()));
867
868 // Slip factor smoothly switching between slip and no slip conditions.
869 RF slip_factor = 0.0;
870 using BoundarySlipSwitch = NavierStokesDGImp::VariableBoundarySlipSwitch<PRM>;
871 if (bctype == BC::SlipVelocity)
872 // Calls boundarySlip(..) function of parameter
873 // class if available, i.e. if
874 // enable_variable_slip is defined. Otherwise
875 // returns 1.0;
876 slip_factor = BoundarySlipSwitch::boundarySlip(prm,ig,ip.position());
877
878 // velocity boundary condition
879 if (bctype == BC::VelocityDirichlet || bctype == BC::SlipVelocity)
880 {
881 // on BC::VelocityDirichlet: 1.0 - slip_factor = 1.0
882 auto factor = weight * (1.0 - slip_factor);
883
884 for(unsigned int d = 0; d < dim; ++d) {
885 const auto& lfsv_v = lfsv_pfs_v.child(d);
886
887 auto val = (jacu[d] * normal) * factor * mu;
888 for(unsigned int i=0; i<vsize; i++) {
889 //================================================//
890 // - (\mu \int \nabla u. normal . v)
891 //================================================//
892 r.accumulate(lfsv_v,i, -val * phi_v[i]);
893 r.accumulate(lfsv_v,i, epsilon * mu * (grad_phi_v[i][0] * normal) * u[d] * factor);
894
895 if(full_tensor) {
896 for(unsigned int dd=0; dd<dim; ++dd) {
897 r.accumulate(lfsv_v,i, -mu * jacu[dd][d] * normal[dd] * phi_v[i] * factor);
898 r.accumulate(lfsv_v,i, epsilon * mu * grad_phi_v[i][0][dd] * u[dd] * normal[d] * factor);
899 }
900 }
901
902 //================================================//
903 // \mu \int \sigma / |\gamma|^\beta v u
904 //================================================//
905 r.accumulate(lfsv_v,i, u[d] * phi_v[i] * mu * penalty_factor * factor);
906
907 //================================================//
908 // \int p v n
909 //================================================//
910 r.accumulate(lfsv_v,i, p * phi_v[i] * normal[d] * weight);
911 } // end i
912 } // end d
913
914 //================================================//
915 // \int q u n
916 //================================================//
917 for(unsigned int i=0; i<psize; i++) {
918 r.accumulate(lfsv_p,i, phi_p[i] * (u * normal) * incomp_scaling * weight);
919 }
920 } // Velocity Dirichlet
921
922 if(bctype == BC::SlipVelocity)
923 {
924 auto factor = weight * (slip_factor);
925
926 RF ten_sum = 1.0;
927 if(full_tensor)
928 ten_sum = 2.0;
929
930 for(unsigned int d = 0; d < dim; ++d) {
931 const auto& lfsv_v = lfsv_pfs_v.child(d);
932
933 for(unsigned int i=0; i<vsize; i++) {
934 //================================================//
935 // - (\mu \int \nabla u. normal . v)
936 //================================================//
937 for(unsigned int dd = 0; dd < dim; ++dd)
938 r.accumulate(lfsv_v,i, -ten_sum * mu * (jacu[dd] * normal) * normal[dd] *phi_v[i] * normal[d] * factor);
939 r.accumulate(lfsv_v,i, epsilon * ten_sum * mu * (u * normal) * (grad_phi_v[i][0] * normal) * normal[d] * factor);
940
941 //================================================//
942 // \mu \int \sigma / |\gamma|^\beta v u
943 //================================================//
944 r.accumulate(lfsv_v,i, mu * penalty_factor * (u * normal) * phi_v[i] * normal[d] * factor);
945 } // end i
946 } // end d
947
948 } // Slip Velocity
949 } // end loop quadrature points
950 } // end alpha_boundary
951
952 // jacobian of boundary term
953 template<typename IG, typename LFSU, typename X, typename LFSV,
954 typename LocalMatrix>
955 void jacobian_boundary (const IG& ig,
956 const LFSU& lfsu, const X& x, const LFSV& lfsv,
957 LocalMatrix& mat) const
958 {
959 // dimensions
960 const unsigned int dim = IG::Entity::dimension;
961
962 // subspaces
963 using namespace Indices;
964 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
965 const auto& lfsv_pfs_v = child(lfsv,_0);
966
967 // ... we assume all velocity components are the same
968 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
969 const auto& lfsv_v = child(lfsv_pfs_v,_0);
970 const unsigned int vsize = lfsv_v.size();
971 using LFSV_P = TypeTree::Child<LFSV,_1>;
972 const auto& lfsv_p = child(lfsv,_1);
973 const unsigned int psize = lfsv_p.size();
974
975 // domain and range field type
976 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
977 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
978 using RT = typename BasisSwitch_V::Range;
979 using RF = typename BasisSwitch_V::RangeField;
980 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
981
982 // References to inside cell
983 const auto& cell_inside = ig.inside();
984
985 // Get geometries
986 auto geo = ig.geometry();
987 auto geo_inside = cell_inside.geometry();
988
989 // Get geometry of intersection in local coordinates of cell_inside
990 auto geo_in_inside = ig.geometryInInside();
991
992 // Determine quadrature order
993 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
994 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
995 const int qorder = 2*v_order + det_jac_order + superintegration_order;
996
997 auto epsilon = prm.epsilonIPSymmetryFactor();
998 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
999
1000 // Initialize vectors outside for loop
1001 std::vector<RT> phi_v(vsize);
1002 std::vector<RT> phi_p(psize);
1003 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
1004
1005 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
1006
1007 // loop over quadrature points and integrate normal flux
1008 for (const auto& ip : quadratureRule(geo,qorder))
1009 {
1010 // position of quadrature point in local coordinates of element
1011 auto local = geo_in_inside.global(ip.position());
1012
1013 // value of velocity shape functions
1014 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1015 // and value of pressure shape functions
1016 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1017
1018 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1019 geo_inside, local, grad_phi_v);
1020
1021 auto normal = ig.unitOuterNormal(ip.position());
1022 auto weight = ip.weight()*geo.integrationElement(ip.position());
1023 auto mu = prm.mu(ig,ip.position());
1024
1025 // evaluate boundary condition type
1026 auto bctype(prm.bctype(ig,ip.position()));
1027
1028 // Slip factor smoothly switching between slip and no slip conditions.
1029 RF slip_factor = 0.0;
1030 using BoundarySlipSwitch = NavierStokesDGImp::VariableBoundarySlipSwitch<PRM>;
1031 if (bctype == BC::SlipVelocity)
1032 // Calls boundarySlip(..) function of parameter
1033 // class if available, i.e. if
1034 // enable_variable_slip is defined. Otherwise
1035 // returns 1.0;
1036 slip_factor = BoundarySlipSwitch::boundarySlip(prm,ig,ip.position());
1037
1038 // velocity boundary condition
1039 if (bctype == BC::VelocityDirichlet || bctype == BC::SlipVelocity)
1040 {
1041 // on BC::VelocityDirichlet: 1.0 - slip_factor = 1.0
1042 auto factor = weight * (1.0 - slip_factor);
1043
1044 for(unsigned int d = 0; d < dim; ++d) {
1045 const auto& lfsv_v = lfsv_pfs_v.child(d);
1046
1047 for(unsigned int i=0; i<vsize; i++) {
1048
1049 for(unsigned int j=0; j<vsize; j++) {
1050 //================================================//
1051 // - (\mu \int \nabla u. normal . v)
1052 //================================================//
1053 auto val = ((grad_phi_v[j][0]*normal)*phi_v[i]) * factor * mu;
1054 mat.accumulate(lfsv_v,i,lfsv_v,j, - val);
1055 mat.accumulate(lfsv_v,j,lfsv_v,i, epsilon*val);
1056
1057 // Assemble symmetric part for (grad u)^T
1058 if(full_tensor) {
1059 for(unsigned int dd = 0; dd < dim; ++dd) {
1060 const auto& lfsv_v_dd = lfsv_pfs_v.child(dd);
1061 auto Tval = ((grad_phi_v[j][0][d]*normal[dd])*phi_v[i]) * factor * mu;
1062 mat.accumulate(lfsv_v,i,lfsv_v_dd,j, - Tval);
1063 mat.accumulate(lfsv_v_dd,j,lfsv_v,i, epsilon*Tval);
1064 }
1065 }
1066 //================================================//
1067 // \mu \int \sigma / |\gamma|^\beta v u
1068 //================================================//
1069 mat.accumulate(lfsv_v,j,lfsv_v,i, phi_v[i] * phi_v[j] * mu * penalty_factor * factor);
1070 }
1071
1072 //================================================//
1073 // \int q u n
1074 // \int p v n
1075 //================================================//
1076 for(unsigned int j=0; j<psize; j++) {
1077 mat.accumulate(lfsv_p,j,lfsv_v,i, (phi_p[j]*normal[d]*phi_v[i]) * weight * incomp_scaling);
1078 mat.accumulate(lfsv_v,i,lfsv_p,j, (phi_p[j]*normal[d]*phi_v[i]) * weight);
1079 }
1080 } // end i
1081 } // end d
1082
1083 } // Velocity Dirichlet
1084
1085 if (bctype == BC::SlipVelocity)
1086 {
1087 auto factor = weight * (slip_factor);
1088
1089 //================================================//
1090 // - (\mu \int \nabla u. normal . v)
1091 //================================================//
1092
1093 for (unsigned int i=0;i<vsize;++i) // ansatz
1094 {
1095 for (unsigned int j=0;j<vsize;++j) // test
1096 {
1097 RF ten_sum = 1.0;
1098
1099 // Assemble symmetric part for (grad u)^T
1100 if(full_tensor)
1101 ten_sum = 2.0;
1102
1103 auto val = ten_sum * ((grad_phi_v[j][0]*normal)*phi_v[i]) * factor * mu;
1104 for (unsigned int d=0;d<dim;++d)
1105 {
1106 const auto& lfsv_v_d = lfsv_pfs_v.child(d);
1107
1108 for (unsigned int dd=0;dd<dim;++dd)
1109 {
1110 const auto& lfsv_v_dd = lfsv_pfs_v.child(dd);
1111
1112 mat.accumulate(lfsv_v_dd,i,lfsv_v_d,j, -val*normal[d]*normal[dd]);
1113 mat.accumulate(lfsv_v_d,j,lfsv_v_dd,i, epsilon*val*normal[d]*normal[dd]);
1114 }
1115 }
1116 }
1117 }
1118
1119 //================================================//
1120 // \mu \int \sigma / |\gamma|^\beta v u
1121 //================================================//
1122 auto p_factor = mu * penalty_factor * factor;
1123 for (unsigned int i=0;i<vsize;++i)
1124 {
1125 for (unsigned int j=0;j<vsize;++j)
1126 {
1127 auto val = phi_v[i]*phi_v[j] * p_factor;
1128 for (unsigned int d=0;d<dim;++d)
1129 {
1130 const auto& lfsv_v_d = lfsv_pfs_v.child(d);
1131 for (unsigned int dd=0;dd<dim;++dd)
1132 {
1133 const auto& lfsv_v_dd = lfsv_pfs_v.child(dd);
1134 mat.accumulate(lfsv_v_d,j,lfsv_v_dd,i, val*normal[d]*normal[dd]);
1135 }
1136 }
1137 }
1138 }
1139
1140 } // Slip Velocity
1141 } // end loop quadrature points
1142 } // end jacobian_boundary
1143
1144 // volume integral depending only on test functions,
1145 // contains f on the right hand side
1146 template<typename EG, typename LFSV, typename R>
1147 void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
1148 {
1149 // dimensions
1150 static const unsigned int dim = EG::Geometry::mydimension;
1151
1152 // subspaces
1153 using namespace Indices;
1154 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
1155 const auto& lfsv_pfs_v = child(lfsv,_0);
1156
1157 // we assume all velocity components are the same type
1158 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
1159 const auto& lfsv_v = child(lfsv_pfs_v,_0);
1160 const unsigned int vsize = lfsv_v.size();
1161 using LFSV_P = TypeTree::Child<LFSV,_1>;
1162 const auto& lfsv_p = child(lfsv,_1);
1163 const unsigned int psize = lfsv_p.size();
1164
1165 // domain and range field type
1166 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
1167 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
1168 using RT = typename BasisSwitch_V::Range;
1169 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
1170 using size_type = typename LFSV::Traits::SizeType;
1171
1172 // Get cell
1173 const auto& cell = eg.entity();
1174
1175 // Get geometries
1176 auto geo = eg.geometry();
1177
1178 // Determine quadrature order
1179 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1180 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
1181 const int qorder = 2*v_order + det_jac_order + superintegration_order;
1182
1183 // Initialize vectors outside for loop
1184 std::vector<RT> phi_v(vsize);
1185 std::vector<RT> phi_p(psize);
1186
1187 // loop over quadrature points
1188 for (const auto& ip : quadratureRule(geo,qorder))
1189 {
1190 auto local = ip.position();
1191 //const Dune::FieldVector<DF,dimw> global = eg.geometry().global(local);
1192
1193 // values of velocity shape functions
1194 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1195
1196 // values of pressure shape functions
1197 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1198
1199 auto weight = ip.weight() * geo.integrationElement(ip.position());
1200
1201 // evaluate source term
1202 auto fval(prm.f(cell,local));
1203
1204 //================================================//
1205 // \int (f*v)
1206 //================================================//
1207 auto factor = weight;
1208 for (unsigned int d=0; d<dim; d++) {
1209 const auto& lfsv_v = lfsv_pfs_v.child(d);
1210
1211 // and store for each velocity component
1212 for (size_type i=0; i<vsize; i++) {
1213 auto val = phi_v[i]*factor;
1214 r.accumulate(lfsv_v,i, -fval[d] * val);
1215 }
1216 }
1217
1218 auto g2 = prm.g2(eg,ip.position());
1219
1220 // integrate div u * psi_i
1221 for (size_t i=0; i<lfsv_p.size(); i++) {
1222 r.accumulate(lfsv_p,i, g2 * phi_p[i] * factor);
1223 }
1224
1225 } // end loop quadrature points
1226 } // end lambda_volume
1227
1228 // boundary integral independent of ansatz functions,
1229 // Neumann and Dirichlet boundary conditions, DG penalty term's right hand side
1230 template<typename IG, typename LFSV, typename R>
1231 void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r) const
1232 {
1233 // dimensions
1234 static const unsigned int dim = IG::Entity::dimension;
1235
1236 // subspaces
1237 using namespace Indices;
1238 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
1239 const auto& lfsv_pfs_v = child(lfsv,_0);
1240
1241 // ... we assume all velocity components are the same
1242 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
1243 const auto& lfsv_v = child(lfsv_pfs_v,_0);
1244 const unsigned int vsize = lfsv_v.size();
1245 using LFSV_P = TypeTree::Child<LFSV,_1>;
1246 const auto& lfsv_p = child(lfsv,_1);
1247 const unsigned int psize = lfsv_p.size();
1248
1249 // domain and range field type
1250 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
1251 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
1252 using DF = typename BasisSwitch_V::DomainField;
1253 using RT = typename BasisSwitch_V::Range;
1254 using RF = typename BasisSwitch_V::RangeField;
1255 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
1256
1257 // References to inside cell
1258 const auto& cell_inside = ig.inside();
1259
1260 // Get geometries
1261 auto geo = ig.geometry();
1262 auto geo_inside = cell_inside.geometry();
1263
1264 // Get geometry of intersection in local coordinates of cell_inside
1265 auto geo_in_inside = ig.geometryInInside();
1266
1267 // Determine quadrature order
1268 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1269 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
1270 const int jac_order = geo.type().isSimplex() ? 0 : 1;
1271 const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
1272
1273 auto epsilon = prm.epsilonIPSymmetryFactor();
1274 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
1275
1276 // Initialize vectors outside for loop
1277 std::vector<RT> phi_v(vsize);
1278 std::vector<RT> phi_p(psize);
1279 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
1280
1281 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
1282
1283 // loop over quadrature points and integrate normal flux
1284 for (const auto& ip : quadratureRule(geo,qorder))
1285 {
1286 // position of quadrature point in local coordinates of element
1287 Dune::FieldVector<DF,dim-1> flocal = ip.position();
1288 Dune::FieldVector<DF,dim> local = geo_in_inside.global(flocal);
1289 //Dune::FieldVector<DF,dimw> global = ig.geometry().global(flocal);
1290
1291 // value of velocity shape functions
1292 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1293 // and value of pressure shape functions
1294 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1295
1296 // compute gradients
1297 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1298 geo_inside, local, grad_phi_v);
1299
1300 auto normal = ig.unitOuterNormal(ip.position());
1301 auto weight = ip.weight()*geo.integrationElement(ip.position());
1302 auto mu = prm.mu(ig,flocal);
1303
1304 // evaluate boundary condition type
1305 auto bctype(prm.bctype(ig,flocal));
1306
1307 if (bctype == BC::VelocityDirichlet)
1308 {
1309 auto u0(prm.g(cell_inside,local));
1310
1311 auto factor = mu * weight;
1312 for(unsigned int d = 0; d < dim; ++d) {
1313 const auto& lfsv_v = lfsv_pfs_v.child(d);
1314
1315 for(unsigned int i=0; i<vsize; i++) {
1316 //================================================//
1317 // \mu \int \nabla v \cdot u_0 \cdot n
1318 //================================================//
1319 r.accumulate(lfsv_v,i, -epsilon * (grad_phi_v[i][0] * normal) * factor * u0[d]);
1320
1321 // Assemble symmetric part for (grad u)^T
1322 if(full_tensor) {
1323 for(unsigned int dd = 0; dd < dim; ++dd) {
1324 const auto& lfsv_v_dd = lfsv_pfs_v.child(dd);
1325 auto Tval = (grad_phi_v[i][0][d]*normal[dd]) * factor;
1326 r.accumulate(lfsv_v_dd,i, -epsilon * Tval * u0[d]);
1327 }
1328 }
1329 //================================================//
1330 // \int \sigma / |\gamma|^\beta v u_0
1331 //================================================//
1332 r.accumulate(lfsv_v,i, -phi_v[i] * penalty_factor * u0[d] * factor);
1333
1334 } // end i
1335 } // end d
1336
1337 //================================================//
1338 // \int q u_0 n
1339 //================================================//
1340 for (unsigned int i=0;i<psize;++i) // test
1341 {
1342 auto val = phi_p[i]*(u0 * normal) * weight;
1343 r.accumulate(lfsv_p,i, - val * incomp_scaling);
1344 }
1345 } // end BC velocity
1346 if (bctype == BC::StressNeumann)
1347 {
1348 auto stress(prm.j(ig,flocal,normal));
1349
1350 //std::cout << "Pdirichlet\n";
1351 //================================================//
1352 // \int p u n
1353 //================================================//
1354 for(unsigned int d = 0; d < dim; ++d) {
1355 const auto& lfsv_v = lfsv_pfs_v.child(d);
1356
1357 for(unsigned int i=0; i<vsize; i++)
1358 r.accumulate(lfsv_v,i, stress[d] * phi_v[i] * weight);
1359 }
1360 }
1361 } // end loop quadrature points
1362 } // end lambda_boundary
1363
1364 private :
1365 PRM& prm;
1366 const int superintegration_order;
1367 Real current_dt;
1368 }; // end class DGNavierStokes
1369
1370 } // end namespace PDELab
1371} // end namespace Dune
1372#endif // DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKES_HH
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
A local operator for solving the Navier-Stokes equations using a DG discretization.
Definition: dgnavierstokes.hh:36
DGNavierStokes(PRM &_prm, int _superintegration_order=0)
Constructor.
Definition: dgnavierstokes.hh:74
sparsity pattern generator
Definition: pattern.hh:30
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 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
Various parser methods to get data into a ParameterTree object.
Definition: stokesparameter.hh:32
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)