4#ifndef DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKES_HH
5#define DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKES_HH
11#include <dune/geometry/referenceelements.hh>
13#include <dune/localfunctions/common/interfaceswitch.hh>
14#include <dune/pdelab/localoperator/idefault.hh>
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>
31 template<
typename PRM>
38 using RF =
typename PRM::Traits::RangeField;
41 using Real =
typename InstatBase::RealType;
43 static const bool navier = PRM::assemble_navier;
44 static const bool full_tensor = PRM::assemble_full_tensor;
49 enum { doPatternVolume =
true };
50 enum { doPatternSkeleton =
true };
53 enum { doSkeletonTwoSided =
false };
56 enum { doAlphaVolume =
true };
57 enum { doAlphaSkeleton =
true };
58 enum { doAlphaBoundary =
true };
59 enum { doLambdaVolume =
true };
60 enum { doLambdaBoundary =
true };
75 prm(_prm), superintegration_order(_superintegration_order),
80 void preStep (Real , Real dt,
int )
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
97 const unsigned int dim = EG::Geometry::mydimension;
101 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
102 const auto& lfsv_pfs_v =
child(lfsv,
_0);
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();
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;
121 auto geo = eg.geometry();
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;
129 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
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);
139 for (
const auto& ip : quadratureRule(geo,qorder))
141 auto local = ip.position();
142 auto mu = prm.mu(eg,local);
143 auto rho = prm.rho(eg,local);
147 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
149 for(
unsigned int d=0; d<dim; ++d) {
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];
158 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
162 for(size_type i=0; i<psize; i++)
163 p += x(lfsv_p,i) * phi_p[i];
166 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
167 geo, local, grad_phi_v);
170 for(
unsigned int d = 0; d<dim; ++d) {
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]);
177 auto detj = geo.integrationElement(ip.position());
178 auto weight = ip.weight() * detj;
180 for(
unsigned int d = 0; d<dim; ++d) {
181 const auto& lfsv_v = lfsv_pfs_v.child(d);
183 for(size_type i=0; i<vsize; i++) {
187 r.accumulate(lfsv_v,i, mu * (jacu[d]*grad_phi_v[i][0]) * weight);
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);
197 r.accumulate(lfsv_v,i, -p * grad_phi_v[i][0][d] * weight);
204 auto u_nabla_u = vu * jacu[d];
206 r.accumulate(lfsv_v,i, rho * u_nabla_u * phi_v[i] * weight);
215 for(size_type i=0; i<psize; i++)
216 for(
unsigned int d = 0; d < dim; ++d)
218 r.accumulate(lfsv_p,i, -jacu[d][d] * phi_p[i] * incomp_scaling * weight);
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
230 const unsigned int dim = EG::Geometry::mydimension;
234 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
235 const auto& lfsv_pfs_v =
child(lfsv,
_0);
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();
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;
254 auto geo = eg.geometry();
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;
262 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
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);
272 for (
const auto& ip : quadratureRule(geo,qorder))
274 auto local = ip.position();
275 auto mu = prm.mu(eg,local);
276 auto rho = prm.rho(eg,local);
280 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
282 for(
unsigned int d=0; d<dim; ++d) {
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];
291 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
294 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
295 geo, local, grad_phi_v);
297 auto detj = geo.integrationElement(ip.position());
298 auto weight = ip.weight() * detj;
300 for(
unsigned int dv = 0; dv<dim; ++dv) {
301 const LFSV_V& lfsv_v = lfsv_pfs_v.child(dv);
306 for(size_type l=0; l<vsize; ++l)
307 gradu_dv.axpy(x(lfsv_v,l), grad_phi_v[l][0]);
309 for(size_type i=0; i<vsize; i++) {
311 for(size_type j=0; j<vsize; j++) {
315 mat.accumulate(lfsv_v,i,lfsv_v,j, mu * (grad_phi_v[j][0]*grad_phi_v[i][0]) * weight);
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);
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);
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);
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);
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
367 const unsigned int dim = IG::Entity::dimension;
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);
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();
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 >;
395 const auto& cell_inside = ig.inside();
396 const auto& cell_outside = ig.outside();
399 auto geo = ig.geometry();
400 auto geo_inside = cell_inside.geometry();
401 auto geo_outside = cell_outside.geometry();
404 auto geo_in_inside = ig.geometryInInside();
405 auto geo_in_outside = ig.geometryInOutside();
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;
412 const int epsilon = prm.epsilonIPSymmetryFactor();
413 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
416 std::vector<RT> phi_v_s(vsize_s);
417 std::vector<RT> phi_v_n(vsize_n);
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);
425 auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
428 for (
const auto& ip : quadratureRule(geo,qorder))
432 auto local_s = geo_in_inside.global(ip.position());
433 auto local_n = geo_in_outside.global(ip.position());
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);
440 assert(vsize_s == vsize_n);
441 for(
unsigned int d=0; d<dim; ++d) {
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];
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);
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];
465 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_s_v.finiteElement()),
466 geo_inside, local_s, grad_phi_v_s);
468 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_n_v.finiteElement()),
469 geo_outside, local_n, grad_phi_v_n);
472 for(
unsigned int d=0; d<dim; ++d) {
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]);
483 auto normal = ig.unitOuterNormal(ip.position());
484 auto weight = ip.weight()*geo.integrationElement(ip.position());
485 auto mu = prm.mu(ig,ip.position());
487 auto factor = mu * weight;
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);
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]);
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]);
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);
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);
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);
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);
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);
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
566 const unsigned int dim = IG::Entity::dimension;
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);
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();
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 >;
594 auto const& cell_inside = ig.inside();
595 auto const& cell_outside = ig.outside();
598 auto geo = ig.geometry();
599 auto geo_inside = cell_inside.geometry();
600 auto geo_outside = cell_outside.geometry();
603 auto geo_in_inside = ig.geometryInInside();
604 auto geo_in_outside = ig.geometryInOutside();
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;
611 const int epsilon = prm.epsilonIPSymmetryFactor();
612 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
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);
622 auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
625 for (
const auto& ip : quadratureRule(geo,qorder))
629 auto local_s = geo_in_inside.global(ip.position());
630 auto local_n = geo_in_outside.global(ip.position());
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);
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);
640 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_s_v.finiteElement()),
641 geo_inside, local_s, grad_phi_v_s);
643 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_n_v.finiteElement()),
644 geo_outside, local_n, grad_phi_v_n);
646 auto normal = ig.unitOuterNormal(ip.position());
647 auto weight = ip.weight()*geo.integrationElement(ip.position());
648 auto mu = prm.mu(ig,ip.position());
650 assert(vsize_s == vsize_n);
651 auto factor = mu * weight;
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);
661 for(
unsigned int i=0; i<vsize_s; ++i) {
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);
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 );
680 for(
unsigned int j=0; j<vsize_n; ++j) {
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);
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);
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);
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 );
715 for(
unsigned int j=0; j<vsize_n; ++j) {
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);
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);
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);
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);
749 for (
unsigned int j=0; j<psize_s;++j) {
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);
756 for (
unsigned int j=0; j<psize_n;++j) {
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);
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,
775 const unsigned int dim = IG::Entity::dimension;
779 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
780 const auto& lfsv_pfs_v =
child(lfsv,
_0);
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();
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 >;
798 const auto& cell_inside = ig.inside();
801 auto geo = ig.geometry();
802 auto geo_inside = cell_inside.geometry();
805 auto geo_in_inside = ig.geometryInInside();
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;
812 auto epsilon = prm.epsilonIPSymmetryFactor();
813 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
816 std::vector<RT> phi_v(vsize);
818 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
821 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
824 for (
const auto& ip : quadratureRule(geo,qorder))
827 auto local = geo_in_inside.global(ip.position());
830 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
833 for(
unsigned int d=0; d<dim; ++d) {
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];
841 std::vector<RT> phi_p(psize);
842 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
846 for(
unsigned int i=0; i<psize; i++)
847 p += x(lfsv_p,i) * phi_p[i];
850 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
851 geo_inside, local, grad_phi_v);
854 for(
unsigned int d=0; d<dim; ++d) {
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]);
861 auto normal = ig.unitOuterNormal(ip.position());
862 auto weight = ip.weight()*geo.integrationElement(ip.position());
863 auto mu = prm.mu(ig,ip.position());
866 auto bctype(prm.bctype(ig,ip.position()));
869 RF slip_factor = 0.0;
870 using BoundarySlipSwitch = NavierStokesDGImp::VariableBoundarySlipSwitch<PRM>;
871 if (bctype == BC::SlipVelocity)
876 slip_factor = BoundarySlipSwitch::boundarySlip(prm,ig,ip.position());
879 if (bctype == BC::VelocityDirichlet || bctype == BC::SlipVelocity)
882 auto factor = weight * (1.0 - slip_factor);
884 for(
unsigned int d = 0; d < dim; ++d) {
885 const auto& lfsv_v = lfsv_pfs_v.child(d);
887 auto val = (jacu[d] * normal) * factor * mu;
888 for(
unsigned int i=0; i<vsize; i++) {
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);
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);
905 r.accumulate(lfsv_v,i, u[d] * phi_v[i] * mu * penalty_factor * factor);
910 r.accumulate(lfsv_v,i, p * phi_v[i] * normal[d] * weight);
917 for(
unsigned int i=0; i<psize; i++) {
918 r.accumulate(lfsv_p,i, phi_p[i] * (u * normal) * incomp_scaling * weight);
922 if(bctype == BC::SlipVelocity)
924 auto factor = weight * (slip_factor);
930 for(
unsigned int d = 0; d < dim; ++d) {
931 const auto& lfsv_v = lfsv_pfs_v.child(d);
933 for(
unsigned int i=0; i<vsize; i++) {
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);
944 r.accumulate(lfsv_v,i, mu * penalty_factor * (u * normal) * phi_v[i] * normal[d] * factor);
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
960 const unsigned int dim = IG::Entity::dimension;
964 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
965 const auto& lfsv_pfs_v =
child(lfsv,
_0);
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();
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 >;
983 const auto& cell_inside = ig.inside();
986 auto geo = ig.geometry();
987 auto geo_inside = cell_inside.geometry();
990 auto geo_in_inside = ig.geometryInInside();
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;
997 auto epsilon = prm.epsilonIPSymmetryFactor();
998 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
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);
1005 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
1008 for (
const auto& ip : quadratureRule(geo,qorder))
1011 auto local = geo_in_inside.global(ip.position());
1014 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1016 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1018 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1019 geo_inside, local, grad_phi_v);
1021 auto normal = ig.unitOuterNormal(ip.position());
1022 auto weight = ip.weight()*geo.integrationElement(ip.position());
1023 auto mu = prm.mu(ig,ip.position());
1026 auto bctype(prm.bctype(ig,ip.position()));
1029 RF slip_factor = 0.0;
1030 using BoundarySlipSwitch = NavierStokesDGImp::VariableBoundarySlipSwitch<PRM>;
1031 if (bctype == BC::SlipVelocity)
1036 slip_factor = BoundarySlipSwitch::boundarySlip(prm,ig,ip.position());
1039 if (bctype == BC::VelocityDirichlet || bctype == BC::SlipVelocity)
1042 auto factor = weight * (1.0 - slip_factor);
1044 for(
unsigned int d = 0; d < dim; ++d) {
1045 const auto& lfsv_v = lfsv_pfs_v.child(d);
1047 for(
unsigned int i=0; i<vsize; i++) {
1049 for(
unsigned int j=0; j<vsize; j++) {
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);
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);
1069 mat.accumulate(lfsv_v,j,lfsv_v,i, phi_v[i] * phi_v[j] * mu * penalty_factor * factor);
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);
1085 if (bctype == BC::SlipVelocity)
1087 auto factor = weight * (slip_factor);
1093 for (
unsigned int i=0;i<vsize;++i)
1095 for (
unsigned int j=0;j<vsize;++j)
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)
1106 const auto& lfsv_v_d = lfsv_pfs_v.child(d);
1108 for (
unsigned int dd=0;dd<dim;++dd)
1110 const auto& lfsv_v_dd = lfsv_pfs_v.child(dd);
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]);
1122 auto p_factor = mu * penalty_factor * factor;
1123 for (
unsigned int i=0;i<vsize;++i)
1125 for (
unsigned int j=0;j<vsize;++j)
1127 auto val = phi_v[i]*phi_v[j] * p_factor;
1128 for (
unsigned int d=0;d<dim;++d)
1130 const auto& lfsv_v_d = lfsv_pfs_v.child(d);
1131 for (
unsigned int dd=0;dd<dim;++dd)
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]);
1146 template<
typename EG,
typename LFSV,
typename R>
1147 void lambda_volume (
const EG& eg,
const LFSV& lfsv, R& r)
const
1150 static const unsigned int dim = EG::Geometry::mydimension;
1154 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
1155 const auto& lfsv_pfs_v =
child(lfsv,
_0);
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();
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;
1173 const auto& cell = eg.entity();
1176 auto geo = eg.geometry();
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;
1184 std::vector<RT> phi_v(vsize);
1185 std::vector<RT> phi_p(psize);
1188 for (
const auto& ip : quadratureRule(geo,qorder))
1190 auto local = ip.position();
1194 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1197 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1199 auto weight = ip.weight() * geo.integrationElement(ip.position());
1202 auto fval(prm.f(cell,local));
1207 auto factor = weight;
1208 for (
unsigned int d=0; d<dim; d++) {
1209 const auto& lfsv_v = lfsv_pfs_v.child(d);
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);
1218 auto g2 = prm.g2(eg,ip.position());
1221 for (
size_t i=0; i<lfsv_p.size(); i++) {
1222 r.accumulate(lfsv_p,i, g2 * phi_p[i] * factor);
1230 template<
typename IG,
typename LFSV,
typename R>
1231 void lambda_boundary (
const IG& ig,
const LFSV& lfsv, R& r)
const
1234 static const unsigned int dim = IG::Entity::dimension;
1238 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
1239 const auto& lfsv_pfs_v =
child(lfsv,
_0);
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();
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 >;
1258 const auto& cell_inside = ig.inside();
1261 auto geo = ig.geometry();
1262 auto geo_inside = cell_inside.geometry();
1265 auto geo_in_inside = ig.geometryInInside();
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;
1273 auto epsilon = prm.epsilonIPSymmetryFactor();
1274 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
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);
1281 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
1284 for (
const auto& ip : quadratureRule(geo,qorder))
1292 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1294 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1297 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1298 geo_inside, local, grad_phi_v);
1300 auto normal = ig.unitOuterNormal(ip.position());
1301 auto weight = ip.weight()*geo.integrationElement(ip.position());
1302 auto mu = prm.mu(ig,flocal);
1305 auto bctype(prm.bctype(ig,flocal));
1307 if (bctype == BC::VelocityDirichlet)
1309 auto u0(prm.g(cell_inside,local));
1311 auto factor = mu * weight;
1312 for(
unsigned int d = 0; d < dim; ++d) {
1313 const auto& lfsv_v = lfsv_pfs_v.child(d);
1315 for(
unsigned int i=0; i<vsize; i++) {
1319 r.accumulate(lfsv_v,i, -epsilon * (grad_phi_v[i][0] * normal) * factor * u0[d]);
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]);
1332 r.accumulate(lfsv_v,i, -phi_v[i] * penalty_factor * u0[d] * factor);
1340 for (
unsigned int i=0;i<psize;++i)
1342 auto val = phi_p[i]*(u0 * normal) * weight;
1343 r.accumulate(lfsv_p,i, - val * incomp_scaling);
1346 if (bctype == BC::StressNeumann)
1348 auto stress(prm.j(ig,flocal,normal));
1354 for(
unsigned int d = 0; d < dim; ++d) {
1355 const auto& lfsv_v = lfsv_pfs_v.child(d);
1357 for(
unsigned int i=0; i<vsize; i++)
1358 r.accumulate(lfsv_v,i, stress[d] * phi_v[i] * weight);
1366 const int superintegration_order;
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
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