4 #ifndef DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKES_HH 5 #define DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKES_HH 7 #include <dune/common/exceptions.hh> 8 #include <dune/common/fvector.hh> 9 #include <dune/common/parametertreeparser.hh> 11 #include <dune/geometry/quadraturerules.hh> 13 #include <dune/localfunctions/common/interfaceswitch.hh> 32 template<
typename PRM>
39 using RF =
typename PRM::Traits::RangeField;
44 static const bool navier = PRM::assemble_navier;
45 static const bool full_tensor = PRM::assemble_full_tensor;
76 prm(_prm), superintegration_order(_superintegration_order),
94 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
95 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 98 const unsigned int dim = EG::Geometry::mydimension;
101 using namespace Indices;
102 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
103 const auto& lfsv_pfs_v = child(lfsv,_0);
106 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
107 const auto& lfsv_v = child(lfsv_pfs_v,_0);
108 const unsigned int vsize = lfsv_v.size();
109 using LFSV_P = TypeTree::Child<LFSV,_1>;
110 const auto& lfsv_p = child(lfsv,_1);
111 const unsigned int psize = lfsv_p.size();
114 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
115 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
116 using RT =
typename BasisSwitch_V::Range;
117 using RF =
typename BasisSwitch_V::RangeField;
118 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
119 using size_type =
typename LFSV::Traits::SizeType;
122 auto geo = eg.geometry();
125 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
126 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
127 const int jac_order = geo.type().isSimplex() ? 0 : 1;
128 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
130 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
133 std::vector<RT> phi_v(vsize);
134 Dune::FieldVector<RF,dim> vu(0.0);
135 std::vector<RT> phi_p(psize);
136 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
137 Dune::FieldMatrix<RF,dim,dim> jacu(0.0);
142 auto local = ip.position();
143 auto mu = prm.mu(eg,local);
144 auto rho = prm.rho(eg,local);
148 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
150 for(
unsigned int d=0; d<
dim; ++d) {
152 const auto& lfsu_v = lfsv_pfs_v.child(d);
153 for(size_type i=0; i<vsize; i++)
154 vu[d] += x(lfsu_v,i) * phi_v[i];
159 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
163 for(size_type i=0; i<psize; i++)
164 p += x(lfsv_p,i) * phi_p[i];
167 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
168 geo, local, grad_phi_v);
171 for(
unsigned int d = 0; d<
dim; ++d) {
173 const auto& lfsv_v = lfsv_pfs_v.child(d);
174 for(size_type i=0; i<vsize; i++)
175 jacu[d].axpy(x(lfsv_v,i), grad_phi_v[i][0]);
178 auto detj = geo.integrationElement(ip.position());
179 auto weight = ip.weight() * detj;
181 for(
unsigned int d = 0; d<
dim; ++d) {
182 const auto& lfsv_v = lfsv_pfs_v.child(d);
184 for(size_type i=0; i<vsize; i++) {
188 r.accumulate(lfsv_v,i, mu * (jacu[d]*grad_phi_v[i][0]) * weight);
192 for(
unsigned int dd = 0; dd<
dim; ++dd)
193 r.accumulate(lfsv_v,i, mu * jacu[dd][d] * grad_phi_v[i][0][dd] * weight);
198 r.accumulate(lfsv_v,i, -p * grad_phi_v[i][0][d] * weight);
205 auto u_nabla_u = vu * jacu[d];
207 r.accumulate(lfsv_v,i, rho * u_nabla_u * phi_v[i] * weight);
216 for(size_type i=0; i<psize; i++)
217 for(
unsigned int d = 0; d <
dim; ++d)
219 r.accumulate(lfsv_p,i, -jacu[d][d] * phi_p[i] * incomp_scaling * weight);
225 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
231 const unsigned int dim = EG::Geometry::mydimension;
234 using namespace Indices;
235 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
236 const auto& lfsv_pfs_v = child(lfsv,_0);
239 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
240 const auto& lfsv_v = child(lfsv_pfs_v,_0);
241 const unsigned int vsize = lfsv_v.size();
242 using LFSV_P = TypeTree::Child<LFSV,_1>;
243 const auto& lfsv_p = child(lfsv,_1);
244 const unsigned int psize = lfsv_p.size();
247 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
248 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
249 using RT =
typename BasisSwitch_V::Range;
250 using RF =
typename BasisSwitch_V::RangeField;
251 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
252 using size_type =
typename LFSV::Traits::SizeType;
255 auto geo = eg.geometry();
258 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
259 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
260 const int jac_order = geo.type().isSimplex() ? 0 : 1;
261 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
263 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
266 std::vector<RT> phi_v(vsize);
267 Dune::FieldVector<RF,dim> vu(0.0);
268 std::vector<RT> phi_p(psize);
269 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
270 Dune::FieldVector<RF,dim> gradu_dv(0.0);
275 auto local = ip.position();
276 auto mu = prm.mu(eg,local);
277 auto rho = prm.rho(eg,local);
281 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
283 for(
unsigned int d=0; d<
dim; ++d) {
285 const auto& lfsu_v = lfsv_pfs_v.child(d);
286 for(size_type i=0; i<vsize; i++)
287 vu[d] += x(lfsu_v,i) * phi_v[i];
292 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
295 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
296 geo, local, grad_phi_v);
298 auto detj = geo.integrationElement(ip.position());
299 auto weight = ip.weight() * detj;
301 for(
unsigned int dv = 0; dv<
dim; ++dv) {
302 const LFSV_V& lfsv_v = lfsv_pfs_v.child(dv);
307 for(size_type l=0; l<vsize; ++l)
308 gradu_dv.axpy(x(lfsv_v,l), grad_phi_v[l][0]);
310 for(size_type i=0; i<vsize; i++) {
312 for(size_type j=0; j<vsize; j++) {
316 mat.accumulate(lfsv_v,i,lfsv_v,j, mu * (grad_phi_v[j][0]*grad_phi_v[i][0]) * weight);
320 for(
unsigned int du = 0; du<
dim; ++du) {
321 const auto& lfsu_v = lfsv_pfs_v.child(du);
322 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (grad_phi_v[j][0][dv]*grad_phi_v[i][0][du]) * weight);
330 for(size_type j=0; j<psize; j++) {
331 mat.accumulate(lfsv_p,j,lfsv_v,i, -phi_p[j] * grad_phi_v[i][0][dv] * incomp_scaling * weight);
332 mat.accumulate(lfsv_v,i,lfsv_p,j, -phi_p[j] * grad_phi_v[i][0][dv] * weight);
341 for(size_type j=0; j<vsize; j++)
342 mat.accumulate(lfsv_v,i,lfsv_v,j, rho * (vu * grad_phi_v[j][0]) * phi_v[i] * weight);
345 for(
unsigned int du = 0; du <
dim; ++du) {
346 const auto& lfsu_v = lfsv_pfs_v.child(du);
347 for(size_type j=0; j<vsize; j++)
348 mat.accumulate(lfsv_v,i,lfsu_v,j, rho * phi_v[j] * gradu_dv[du] * phi_v[i] * weight);
361 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
363 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
364 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
365 R& r_s, R& r_n)
const 368 const unsigned int dim = IG::dimension;
371 using namespace Indices;
372 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
373 const auto& lfsv_s_pfs_v = child(lfsv_s,_0);
374 const auto& lfsv_n_pfs_v = child(lfsv_n,_0);
377 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
378 const auto& lfsv_s_v = child(lfsv_s_pfs_v,_0);
379 const auto& lfsv_n_v = child(lfsv_n_pfs_v,_0);
380 const unsigned int vsize_s = lfsv_s_v.size();
381 const unsigned int vsize_n = lfsv_n_v.size();
382 using LFSV_P = TypeTree::Child<LFSV,_1>;
383 const auto& lfsv_s_p = child(lfsv_s,_1);
384 const auto& lfsv_n_p = child(lfsv_n,_1);
385 const unsigned int psize_s = lfsv_s_p.size();
386 const unsigned int psize_n = lfsv_n_p.size();
389 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
390 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
391 using RT =
typename BasisSwitch_V::Range;
392 using RF =
typename BasisSwitch_V::RangeField;
393 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
396 const auto& cell_inside = ig.inside();
397 const auto& cell_outside = ig.outside();
400 auto geo = ig.geometry();
401 auto geo_inside = cell_inside.geometry();
402 auto geo_outside = cell_outside.geometry();
405 auto geo_in_inside = ig.geometryInInside();
406 auto geo_in_outside = ig.geometryInOutside();
409 const int v_order = FESwitch_V::basis(lfsv_s_v.finiteElement()).order();
410 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
411 const int qorder = 2*v_order + det_jac_order + superintegration_order;
413 const int epsilon = prm.epsilonIPSymmetryFactor();
414 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
417 std::vector<RT> phi_v_s(vsize_s);
418 std::vector<RT> phi_v_n(vsize_n);
419 Dune::FieldVector<RF,dim> u_s(0.0), u_n(0.0);
420 std::vector<RT> phi_p_s(psize_s);
421 std::vector<RT> phi_p_n(psize_n);
422 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_s(vsize_s);
423 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_n(vsize_n);
424 Dune::FieldMatrix<RF,dim,dim> jacu_s(0.0), jacu_n(0.0);
426 auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
433 auto local_s = geo_in_inside.global(ip.position());
434 auto local_n = geo_in_outside.global(ip.position());
437 FESwitch_V::basis(lfsv_s_v.finiteElement()).evaluateFunction(local_s,phi_v_s);
438 FESwitch_V::basis(lfsv_n_v.finiteElement()).evaluateFunction(local_n,phi_v_n);
441 assert(vsize_s == vsize_n);
442 for(
unsigned int d=0; d<
dim; ++d) {
445 const auto& lfsv_s_v = lfsv_s_pfs_v.child(d);
446 const auto& lfsv_n_v = lfsv_n_pfs_v.child(d);
447 for(
unsigned int i=0; i<vsize_s; i++) {
448 u_s[d] += x_s(lfsv_s_v,i) * phi_v_s[i];
449 u_n[d] += x_n(lfsv_n_v,i) * phi_v_n[i];
454 FESwitch_P::basis(lfsv_s_p.finiteElement()).evaluateFunction(local_s,phi_p_s);
455 FESwitch_P::basis(lfsv_n_p.finiteElement()).evaluateFunction(local_n,phi_p_n);
458 assert(psize_s == psize_n);
459 RF p_s(0.0), p_n(0.0);
460 for(
unsigned int i=0; i<psize_s; i++) {
461 p_s += x_s(lfsv_s_p,i) * phi_p_s[i];
462 p_n += x_n(lfsv_n_p,i) * phi_p_n[i];
466 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_s_v.finiteElement()),
467 geo_inside, local_s, grad_phi_v_s);
469 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_n_v.finiteElement()),
470 geo_outside, local_n, grad_phi_v_n);
473 for(
unsigned int d=0; d<
dim; ++d) {
476 const auto& lfsv_s_v = lfsv_s_pfs_v.child(d);
477 const auto& lfsv_n_v = lfsv_n_pfs_v.child(d);
478 for(
unsigned int i=0; i<vsize_s; i++) {
479 jacu_s[d].axpy(x_s(lfsv_s_v,i), grad_phi_v_s[i][0]);
480 jacu_n[d].axpy(x_n(lfsv_n_v,i), grad_phi_v_n[i][0]);
484 auto normal = ig.unitOuterNormal(ip.position());
485 auto weight = ip.weight()*geo.integrationElement(ip.position());
486 auto mu = prm.mu(ig,ip.position());
488 auto factor = mu * weight;
490 for(
unsigned int d=0; d<
dim; ++d) {
491 const auto& lfsv_s_v = lfsv_s_pfs_v.child(d);
492 const auto& lfsv_n_v = lfsv_n_pfs_v.child(d);
497 auto val = 0.5 * ((jacu_s[d] * normal) + (jacu_n[d] * normal)) * factor;
498 for(
unsigned int i=0; i<vsize_s; i++) {
499 r_s.accumulate(lfsv_s_v,i, -val * phi_v_s[i]);
500 r_n.accumulate(lfsv_n_v,i, val * phi_v_n[i]);
503 for(
unsigned int dd=0; dd<
dim; ++dd) {
504 auto Tval = 0.5 * (jacu_s[dd][d] + jacu_n[dd][d]) * normal[dd] * factor;
505 r_s.accumulate(lfsv_s_v,i, -Tval * phi_v_s[i]);
506 r_n.accumulate(lfsv_n_v,i, Tval * phi_v_n[i]);
514 auto jumpu_d = u_s[d] - u_n[d];
515 for(
unsigned int i=0; i<vsize_s; i++) {
516 r_s.accumulate(lfsv_s_v,i, epsilon * 0.5 * (grad_phi_v_s[i][0] * normal) * jumpu_d * factor);
517 r_n.accumulate(lfsv_n_v,i, epsilon * 0.5 * (grad_phi_v_n[i][0] * normal) * jumpu_d * factor);
520 for(
unsigned int dd=0; dd<
dim; ++dd) {
521 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);
522 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);
530 for(
unsigned int i=0; i<vsize_s; i++) {
531 r_s.accumulate(lfsv_s_v,i, penalty_factor * jumpu_d * phi_v_s[i] * factor);
532 r_n.accumulate(lfsv_n_v,i, -penalty_factor * jumpu_d * phi_v_n[i] * factor);
538 auto mean_p = 0.5*(p_s + p_n);
539 for(
unsigned int i=0; i<vsize_s; i++) {
540 r_s.accumulate(lfsv_s_v,i, mean_p * phi_v_s[i] * normal[d] * weight);
541 r_n.accumulate(lfsv_n_v,i, -mean_p * phi_v_n[i] * normal[d] * weight);
548 auto jumpu_n = (u_s*normal) - (u_n*normal);
549 for(
unsigned int i=0; i<psize_s; i++) {
550 r_s.accumulate(lfsv_s_p,i, 0.5 * phi_p_s[i] * jumpu_n * incomp_scaling * weight);
551 r_n.accumulate(lfsv_n_p,i, 0.5 * phi_p_n[i] * jumpu_n * incomp_scaling * weight);
558 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
561 const LFSU& lfsu_s,
const X&,
const LFSV& lfsv_s,
562 const LFSU& lfsu_n,
const X&,
const LFSV& lfsv_n,
567 const unsigned int dim = IG::dimension;
570 using namespace Indices;
571 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
572 const auto& lfsv_s_pfs_v = child(lfsv_s,_0);
573 const auto& lfsv_n_pfs_v = child(lfsv_n,_0);
576 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
577 const auto& lfsv_s_v = child(lfsv_s_pfs_v,_0);
578 const auto& lfsv_n_v = child(lfsv_n_pfs_v,_0);
579 const unsigned int vsize_s = lfsv_s_v.size();
580 const unsigned int vsize_n = lfsv_n_v.size();
581 using LFSV_P = TypeTree::Child<LFSV,_1>;
582 const auto& lfsv_s_p = child(lfsv_s,_1);
583 const auto& lfsv_n_p = child(lfsv_n,_1);
584 const unsigned int psize_s = lfsv_s_p.size();
585 const unsigned int psize_n = lfsv_n_p.size();
588 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
589 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
590 using RT =
typename BasisSwitch_V::Range;
591 using RF =
typename BasisSwitch_V::RangeField;
592 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
595 auto const& cell_inside = ig.inside();
596 auto const& cell_outside = ig.outside();
599 auto geo = ig.geometry();
600 auto geo_inside = cell_inside.geometry();
601 auto geo_outside = cell_outside.geometry();
604 auto geo_in_inside = ig.geometryInInside();
605 auto geo_in_outside = ig.geometryInOutside();
608 const int v_order = FESwitch_V::basis(lfsv_s_v.finiteElement()).order();
609 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
610 const int qorder = 2*v_order + det_jac_order + superintegration_order;
612 const int epsilon = prm.epsilonIPSymmetryFactor();
613 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
616 std::vector<RT> phi_v_s(vsize_s);
617 std::vector<RT> phi_v_n(vsize_n);
618 std::vector<RT> phi_p_s(psize_s);
619 std::vector<RT> phi_p_n(psize_n);
620 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_s(vsize_s);
621 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_n(vsize_n);
623 auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
630 auto local_s = geo_in_inside.global(ip.position());
631 auto local_n = geo_in_outside.global(ip.position());
634 FESwitch_V::basis(lfsv_s_v.finiteElement()).evaluateFunction(local_s,phi_v_s);
635 FESwitch_V::basis(lfsv_n_v.finiteElement()).evaluateFunction(local_n,phi_v_n);
637 FESwitch_P::basis(lfsv_s_p.finiteElement()).evaluateFunction(local_s,phi_p_s);
638 FESwitch_P::basis(lfsv_n_p.finiteElement()).evaluateFunction(local_n,phi_p_n);
641 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_s_v.finiteElement()),
642 geo_inside, local_s, grad_phi_v_s);
644 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_n_v.finiteElement()),
645 geo_outside, local_n, grad_phi_v_n);
647 auto normal = ig.unitOuterNormal(ip.position());
648 auto weight = ip.weight()*geo.integrationElement(ip.position());
649 auto mu = prm.mu(ig,ip.position());
651 assert(vsize_s == vsize_n);
652 auto factor = mu * weight;
654 for(
unsigned int d = 0; d <
dim; ++d) {
655 const auto& lfsv_s_v = lfsv_s_pfs_v.child(d);
656 const auto& lfsv_n_v = lfsv_n_pfs_v.child(d);
662 for(
unsigned int i=0; i<vsize_s; ++i) {
664 for(
unsigned int j=0; j<vsize_s; ++j) {
665 auto val = (0.5*(grad_phi_v_s[i][0]*normal)*phi_v_s[j]) * factor;
666 mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v,i, -val);
667 mat_ss.accumulate(lfsv_s_v,i,lfsv_s_v,j, epsilon * val);
668 mat_ss.accumulate(lfsv_s_v,i,lfsv_s_v,j, phi_v_s[i] * phi_v_s[j] * penalty_factor * factor);
672 for(
unsigned int dd = 0; dd <
dim; ++dd) {
673 auto Tval = (0.5*(grad_phi_v_s[i][0][d]*normal[dd])*phi_v_s[j]) * factor;
674 const auto& lfsv_s_v_dd = lfsv_s_pfs_v.child(dd);
675 mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v_dd,i, - Tval);
676 mat_ss.accumulate(lfsv_s_v_dd,i,lfsv_s_v,j, epsilon*Tval );
681 for(
unsigned int j=0; j<vsize_n; ++j) {
683 auto val = (-0.5*(grad_phi_v_s[i][0]*normal)*phi_v_n[j]) * factor;
684 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v,i,- val);
685 mat_sn.accumulate(lfsv_s_v,i,lfsv_n_v,j, epsilon*val);
686 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v,i, -phi_v_s[i] * phi_v_n[j] * penalty_factor * factor);
690 for (
unsigned int dd=0;dd<
dim;++dd) {
691 auto Tval = (-0.5*(grad_phi_v_s[i][0][d]*normal[dd])*phi_v_n[j]) * factor;
692 const auto& lfsv_s_v_dd = lfsv_s_pfs_v.child(dd);
693 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v_dd,i,- Tval);
694 mat_sn.accumulate(lfsv_s_v_dd,i,lfsv_n_v,j, epsilon*Tval);
699 for(
unsigned int j=0; j<vsize_s; ++j) {
700 auto val = (0.5*(grad_phi_v_n[i][0]*normal)*phi_v_s[j]) * factor;
701 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v,i, - val);
702 mat_ns.accumulate(lfsv_n_v,i,lfsv_s_v,j, epsilon*val );
703 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v,i, -phi_v_n[i] * phi_v_s[j] * penalty_factor * factor);
707 for (
unsigned int dd=0;dd<
dim;++dd) {
708 auto Tval = (0.5*(grad_phi_v_n[i][0][d]*normal[dd])*phi_v_s[j]) * factor;
709 const auto& lfsv_n_v_dd = lfsv_n_pfs_v.child(dd);
710 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v_dd,i, - Tval);
711 mat_ns.accumulate(lfsv_n_v_dd,i,lfsv_s_v,j, epsilon*Tval );
716 for(
unsigned int j=0; j<vsize_n; ++j) {
718 auto val = (-0.5*(grad_phi_v_n[i][0]*normal)*phi_v_n[j]) * factor;
719 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v,i, - val);
720 mat_nn.accumulate(lfsv_n_v,i,lfsv_n_v,j, epsilon*val);
721 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v,i, phi_v_n[i] * phi_v_n[j] * penalty_factor * factor);
725 for (
unsigned int dd=0;dd<
dim;++dd) {
726 auto Tval = (-0.5*(grad_phi_v_n[i][0][d]*normal[dd])*phi_v_n[j]) * factor;
727 const auto& lfsv_n_v_dd = lfsv_n_pfs_v.child(dd);
728 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v_dd,i,- Tval);
729 mat_nn.accumulate(lfsv_n_v_dd,i,lfsv_n_v,j, epsilon*Tval);
738 for(
unsigned int j=0; j<psize_s; ++j) {
739 auto val = 0.5*(phi_p_s[j]*normal[d]*phi_v_s[i]) * weight;
740 mat_ss.accumulate(lfsv_s_v,i,lfsv_s_p,j, val);
741 mat_ss.accumulate(lfsv_s_p,j,lfsv_s_v,i, val * incomp_scaling);
744 for(
unsigned int j=0; j<psize_n; ++j) {
745 auto val = 0.5*(phi_p_n[j]*normal[d]*phi_v_s[i]) * weight;
746 mat_sn.accumulate(lfsv_s_v,i,lfsv_n_p,j, val);
747 mat_ns.accumulate(lfsv_n_p,j,lfsv_s_v,i, val * incomp_scaling);
750 for (
unsigned int j=0; j<psize_s;++j) {
752 auto val = -0.5*(phi_p_s[j]*normal[d]*phi_v_n[i]) * weight;
753 mat_ns.accumulate(lfsv_n_v,i,lfsv_s_p,j, val);
754 mat_sn.accumulate(lfsv_s_p,j,lfsv_n_v,i, val * incomp_scaling);
757 for (
unsigned int j=0; j<psize_n;++j) {
759 auto val = -0.5*(phi_p_n[j]*normal[d]*phi_v_n[i]) * weight;
760 mat_nn.accumulate(lfsv_n_v,i,lfsv_n_p,j, val);
761 mat_nn.accumulate(lfsv_n_p,j,lfsv_n_v,i, val * incomp_scaling);
770 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
772 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
776 const unsigned int dim = IG::dimension;
779 using namespace Indices;
780 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
781 const auto& lfsv_pfs_v = child(lfsv,_0);
784 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
785 const auto& lfsv_v = child(lfsv_pfs_v,_0);
786 const unsigned int vsize = lfsv_v.size();
787 using LFSV_P = TypeTree::Child<LFSV,_1>;
788 const auto& lfsv_p = child(lfsv,_1);
789 const unsigned int psize = lfsv_p.size();
792 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
793 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
794 using RT =
typename BasisSwitch_V::Range;
795 using RF =
typename BasisSwitch_V::RangeField;
796 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
799 const auto& cell_inside = ig.inside();
802 auto geo = ig.geometry();
803 auto geo_inside = cell_inside.geometry();
806 auto geo_in_inside = ig.geometryInInside();
809 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
810 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
811 const int qorder = 2*v_order + det_jac_order + superintegration_order;
813 auto epsilon = prm.epsilonIPSymmetryFactor();
814 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
817 std::vector<RT> phi_v(vsize);
818 Dune::FieldVector<RF,dim> u(0.0);
819 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
820 Dune::FieldMatrix<RF,dim,dim> jacu(0.0);
822 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
828 auto local = geo_in_inside.global(ip.position());
831 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
834 for(
unsigned int d=0; d<
dim; ++d) {
836 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
837 for(
unsigned int i=0; i<vsize; i++)
838 u[d] += x(lfsv_v,i) * phi_v[i];
842 std::vector<RT> phi_p(psize);
843 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
847 for(
unsigned int i=0; i<psize; i++)
848 p += x(lfsv_p,i) * phi_p[i];
851 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
852 geo_inside, local, grad_phi_v);
855 for(
unsigned int d=0; d<
dim; ++d) {
857 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
858 for(
unsigned int i=0; i<vsize; i++)
859 jacu[d].axpy(x(lfsv_v,i), grad_phi_v[i][0]);
862 auto normal = ig.unitOuterNormal(ip.position());
863 auto weight = ip.weight()*geo.integrationElement(ip.position());
864 auto mu = prm.mu(ig,ip.position());
867 auto bctype(prm.bctype(ig,ip.position()));
870 RF slip_factor = 0.0;
877 slip_factor = BoundarySlipSwitch::boundarySlip(prm,ig,ip.position());
883 auto factor = weight * (1.0 - slip_factor);
885 for(
unsigned int d = 0; d <
dim; ++d) {
886 const auto& lfsv_v = lfsv_pfs_v.child(d);
888 auto val = (jacu[d] * normal) * factor * mu;
889 for(
unsigned int i=0; i<vsize; i++) {
893 r.accumulate(lfsv_v,i, -val * phi_v[i]);
894 r.accumulate(lfsv_v,i, epsilon * mu * (grad_phi_v[i][0] * normal) * u[d] * factor);
897 for(
unsigned int dd=0; dd<
dim; ++dd) {
898 r.accumulate(lfsv_v,i, -mu * jacu[dd][d] * normal[dd] * phi_v[i] * factor);
899 r.accumulate(lfsv_v,i, epsilon * mu * grad_phi_v[i][0][dd] * u[dd] * normal[d] * factor);
906 r.accumulate(lfsv_v,i, u[d] * phi_v[i] * mu * penalty_factor * factor);
911 r.accumulate(lfsv_v,i, p * phi_v[i] * normal[d] * weight);
918 for(
unsigned int i=0; i<psize; i++) {
919 r.accumulate(lfsv_p,i, phi_p[i] * (u * normal) * incomp_scaling * weight);
925 auto factor = weight * (slip_factor);
931 for(
unsigned int d = 0; d <
dim; ++d) {
932 const auto& lfsv_v = lfsv_pfs_v.child(d);
934 for(
unsigned int i=0; i<vsize; i++) {
938 for(
unsigned int dd = 0; dd <
dim; ++dd)
939 r.accumulate(lfsv_v,i, -ten_sum * mu * (jacu[dd] * normal) * normal[dd] *phi_v[i] * normal[d] * factor);
940 r.accumulate(lfsv_v,i, epsilon * ten_sum * mu * (u * normal) * (grad_phi_v[i][0] * normal) * normal[d] * factor);
945 r.accumulate(lfsv_v,i, mu * penalty_factor * (u * normal) * phi_v[i] * normal[d] * factor);
954 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
957 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
961 const unsigned int dim = IG::dimension;
964 using namespace Indices;
965 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
966 const auto& lfsv_pfs_v = child(lfsv,_0);
969 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
970 const auto& lfsv_v = child(lfsv_pfs_v,_0);
971 const unsigned int vsize = lfsv_v.size();
972 using LFSV_P = TypeTree::Child<LFSV,_1>;
973 const auto& lfsv_p = child(lfsv,_1);
974 const unsigned int psize = lfsv_p.size();
977 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
978 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
979 using RT =
typename BasisSwitch_V::Range;
980 using RF =
typename BasisSwitch_V::RangeField;
981 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
984 const auto& cell_inside = ig.inside();
987 auto geo = ig.geometry();
988 auto geo_inside = cell_inside.geometry();
991 auto geo_in_inside = ig.geometryInInside();
994 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
995 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
996 const int qorder = 2*v_order + det_jac_order + superintegration_order;
998 auto epsilon = prm.epsilonIPSymmetryFactor();
999 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
1002 std::vector<RT> phi_v(vsize);
1003 std::vector<RT> phi_p(psize);
1004 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
1006 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
1012 auto local = geo_in_inside.global(ip.position());
1015 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1017 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1019 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1020 geo_inside, local, grad_phi_v);
1022 auto normal = ig.unitOuterNormal(ip.position());
1023 auto weight = ip.weight()*geo.integrationElement(ip.position());
1024 auto mu = prm.mu(ig,ip.position());
1027 auto bctype(prm.bctype(ig,ip.position()));
1030 RF slip_factor = 0.0;
1037 slip_factor = BoundarySlipSwitch::boundarySlip(prm,ig,ip.position());
1043 auto factor = weight * (1.0 - slip_factor);
1045 for(
unsigned int d = 0; d <
dim; ++d) {
1046 const auto& lfsv_v = lfsv_pfs_v.child(d);
1048 for(
unsigned int i=0; i<vsize; i++) {
1050 for(
unsigned int j=0; j<vsize; j++) {
1054 auto val = ((grad_phi_v[j][0]*normal)*phi_v[i]) * factor * mu;
1055 mat.accumulate(lfsv_v,i,lfsv_v,j, - val);
1056 mat.accumulate(lfsv_v,j,lfsv_v,i, epsilon*val);
1060 for(
unsigned int dd = 0; dd <
dim; ++dd) {
1061 const auto& lfsv_v_dd = lfsv_pfs_v.child(dd);
1062 auto Tval = ((grad_phi_v[j][0][d]*normal[dd])*phi_v[i]) * factor * mu;
1063 mat.accumulate(lfsv_v,i,lfsv_v_dd,j, - Tval);
1064 mat.accumulate(lfsv_v_dd,j,lfsv_v,i, epsilon*Tval);
1070 mat.accumulate(lfsv_v,j,lfsv_v,i, phi_v[i] * phi_v[j] * mu * penalty_factor * factor);
1077 for(
unsigned int j=0; j<psize; j++) {
1078 mat.accumulate(lfsv_p,j,lfsv_v,i, (phi_p[j]*normal[d]*phi_v[i]) * weight * incomp_scaling);
1079 mat.accumulate(lfsv_v,i,lfsv_p,j, (phi_p[j]*normal[d]*phi_v[i]) * weight);
1088 auto factor = weight * (slip_factor);
1094 for (
unsigned int i=0;i<vsize;++i)
1096 for (
unsigned int j=0;j<vsize;++j)
1104 auto val = ten_sum * ((grad_phi_v[j][0]*normal)*phi_v[i]) * factor * mu;
1105 for (
unsigned int d=0;d<
dim;++d)
1107 const auto& lfsv_v_d = lfsv_pfs_v.child(d);
1109 for (
unsigned int dd=0;dd<
dim;++dd)
1111 const auto& lfsv_v_dd = lfsv_pfs_v.child(dd);
1113 mat.accumulate(lfsv_v_dd,i,lfsv_v_d,j, -val*normal[d]*normal[dd]);
1114 mat.accumulate(lfsv_v_d,j,lfsv_v_dd,i, epsilon*val*normal[d]*normal[dd]);
1123 auto p_factor = mu * penalty_factor * factor;
1124 for (
unsigned int i=0;i<vsize;++i)
1126 for (
unsigned int j=0;j<vsize;++j)
1128 auto val = phi_v[i]*phi_v[j] * p_factor;
1129 for (
unsigned int d=0;d<
dim;++d)
1131 const auto& lfsv_v_d = lfsv_pfs_v.child(d);
1132 for (
unsigned int dd=0;dd<
dim;++dd)
1134 const auto& lfsv_v_dd = lfsv_pfs_v.child(dd);
1135 mat.accumulate(lfsv_v_d,j,lfsv_v_dd,i, val*normal[d]*normal[dd]);
1147 template<
typename EG,
typename LFSV,
typename R>
1151 static const unsigned int dim = EG::Geometry::mydimension;
1154 using namespace Indices;
1155 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
1156 const auto& lfsv_pfs_v = child(lfsv,_0);
1159 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
1160 const auto& lfsv_v = child(lfsv_pfs_v,_0);
1161 const unsigned int vsize = lfsv_v.size();
1162 using LFSV_P = TypeTree::Child<LFSV,_1>;
1163 const auto& lfsv_p = child(lfsv,_1);
1164 const unsigned int psize = lfsv_p.size();
1167 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
1168 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
1169 using RT =
typename BasisSwitch_V::Range;
1170 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
1171 using size_type =
typename LFSV::Traits::SizeType;
1174 const auto& cell = eg.entity();
1177 auto geo = eg.geometry();
1180 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1181 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
1182 const int qorder = 2*v_order + det_jac_order + superintegration_order;
1185 std::vector<RT> phi_v(vsize);
1186 std::vector<RT> phi_p(psize);
1191 auto local = ip.position();
1195 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1198 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1200 auto weight = ip.weight() * geo.integrationElement(ip.position());
1203 auto fval(prm.f(cell,local));
1208 auto factor = weight;
1209 for (
unsigned int d=0; d<
dim; d++) {
1210 const auto& lfsv_v = lfsv_pfs_v.child(d);
1213 for (size_type i=0; i<vsize; i++) {
1214 auto val = phi_v[i]*factor;
1215 r.accumulate(lfsv_v,i, -fval[d] * val);
1219 auto g2 = prm.g2(eg,ip.position());
1222 for (
size_t i=0; i<lfsv_p.size(); i++) {
1223 r.accumulate(lfsv_p,i, g2 * phi_p[i] * factor);
1231 template<
typename IG,
typename LFSV,
typename R>
1235 static const unsigned int dim = IG::dimension;
1238 using namespace Indices;
1239 using LFSV_PFS_V = TypeTree::Child<LFSV,_0>;
1240 const auto& lfsv_pfs_v = child(lfsv,_0);
1243 using LFSV_V = TypeTree::Child<LFSV_PFS_V,_0>;
1244 const auto& lfsv_v = child(lfsv_pfs_v,_0);
1245 const unsigned int vsize = lfsv_v.size();
1246 using LFSV_P = TypeTree::Child<LFSV,_1>;
1247 const auto& lfsv_p = child(lfsv,_1);
1248 const unsigned int psize = lfsv_p.size();
1251 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
1252 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
1253 using DF =
typename BasisSwitch_V::DomainField;
1254 using RT =
typename BasisSwitch_V::Range;
1255 using RF =
typename BasisSwitch_V::RangeField;
1256 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
1259 const auto& cell_inside = ig.inside();
1262 auto geo = ig.geometry();
1263 auto geo_inside = cell_inside.geometry();
1266 auto geo_in_inside = ig.geometryInInside();
1269 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1270 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
1271 const int jac_order = geo.type().isSimplex() ? 0 : 1;
1272 const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
1274 auto epsilon = prm.epsilonIPSymmetryFactor();
1275 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
1278 std::vector<RT> phi_v(vsize);
1279 std::vector<RT> phi_p(psize);
1280 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
1282 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
1288 Dune::FieldVector<DF,dim-1> flocal = ip.position();
1289 Dune::FieldVector<DF,dim> local = geo_in_inside.global(flocal);
1293 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1295 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1298 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1299 geo_inside, local, grad_phi_v);
1301 auto normal = ig.unitOuterNormal(ip.position());
1302 auto weight = ip.weight()*geo.integrationElement(ip.position());
1303 auto mu = prm.mu(ig,flocal);
1306 auto bctype(prm.bctype(ig,flocal));
1310 auto u0(prm.g(cell_inside,local));
1312 auto factor = mu * weight;
1313 for(
unsigned int d = 0; d <
dim; ++d) {
1314 const auto& lfsv_v = lfsv_pfs_v.child(d);
1316 for(
unsigned int i=0; i<vsize; i++) {
1320 r.accumulate(lfsv_v,i, -epsilon * (grad_phi_v[i][0] * normal) * factor * u0[d]);
1324 for(
unsigned int dd = 0; dd <
dim; ++dd) {
1325 const auto& lfsv_v_dd = lfsv_pfs_v.child(dd);
1326 auto Tval = (grad_phi_v[i][0][d]*normal[dd]) * factor;
1327 r.accumulate(lfsv_v_dd,i, -epsilon * Tval * u0[d]);
1333 r.accumulate(lfsv_v,i, -phi_v[i] * penalty_factor * u0[d] * factor);
1341 for (
unsigned int i=0;i<psize;++i)
1343 auto val = phi_p[i]*(u0 * normal) * weight;
1344 r.accumulate(lfsv_p,i, - val * incomp_scaling);
1349 auto stress(prm.j(ig,flocal,normal));
1355 for(
unsigned int d = 0; d <
dim; ++d) {
1356 const auto& lfsv_v = lfsv_pfs_v.child(d);
1358 for(
unsigned int i=0; i<vsize; i++)
1359 r.accumulate(lfsv_v,i, stress[d] * phi_v[i] * weight);
1367 const int superintegration_order;
1373 #endif // DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKES_HH Definition: dgnavierstokes.hh:58
const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
Definition: dgnavierstokes.hh:50
Compile-time switch for the boundary slip factor.
Definition: dgnavierstokesparameter.hh:175
void alpha_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, R &r_s, R &r_n) const
Definition: dgnavierstokes.hh:362
void jacobian_skeleton(const IG &ig, const LFSU &lfsu_s, const X &, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &, const LFSV &lfsv_n, LocalMatrix &mat_ss, LocalMatrix &mat_sn, LocalMatrix &mat_ns, LocalMatrix &mat_nn) const
Definition: dgnavierstokes.hh:560
Definition: stokesparameter.hh:36
Definition: stokesparameter.hh:35
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: dgnavierstokes.hh:227
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
const P & p
Definition: constraints.hh:147
Definition: dgnavierstokes.hh:54
void jacobian_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: dgnavierstokes.hh:956
Definition: dgnavierstokes.hh:57
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: dgnavierstokes.hh:1148
sparsity pattern generator
Definition: pattern.hh:29
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: dgnavierstokes.hh:1232
void preStep(Real, Real dt, int)
Definition: dgnavierstokes.hh:81
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: dgnavierstokes.hh:95
Definition: dgnavierstokes.hh:51
DGNavierStokes(PRM &_prm, int _superintegration_order=0)
Constructor.
Definition: dgnavierstokes.hh:75
Definition: stokesparameter.hh:32
Definition: dgnavierstokes.hh:59
A dense matrix for storing data associated with the degrees of freedom of a pair of LocalFunctionSpac...
Definition: localmatrix.hh:184
sparsity pattern generator
Definition: pattern.hh:13
QuadratureRuleWrapper< QuadratureRule< typename Geometry::ctype, Geometry::mydimension > > quadratureRule(const Geometry &geo, std::size_t order, QuadratureType::Enum quadrature_type=QuadratureType::GaussLegendre)
Returns a quadrature rule for the given geometry.
Definition: quadraturerules.hh:117
void alpha_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: dgnavierstokes.hh:771
Definition: dgnavierstokes.hh:60
Default flags for all local operators.
Definition: flags.hh:18
A local operator for solving the Navier-Stokes equations using a DG discretization.
Definition: dgnavierstokes.hh:33
void setTime(Real t)
Definition: dgnavierstokes.hh:87
Definition: dgnavierstokes.hh:61
Definition: stokesparameter.hh:37
R RealType
Definition: idefault.hh:92