4#ifndef DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESVELVECFEM_HH
5#define DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESVELVECFEM_HH
11#include <dune/localfunctions/common/interfaceswitch.hh>
12#include <dune/pdelab/localoperator/idefault.hh>
14#include <dune/pdelab/common/quadraturerules.hh>
15#include <dune/pdelab/localoperator/defaultimp.hh>
16#include <dune/pdelab/localoperator/pattern.hh>
17#include <dune/pdelab/localoperator/flags.hh>
18#include <dune/pdelab/localoperator/dgnavierstokesparameter.hh>
19#include <dune/pdelab/localoperator/navierstokesmass.hh>
24#define PBLOCK (- VBLOCK + 1)
29 template<
class Basis,
class Dummy =
void>
30 struct VectorBasisInterfaceSwitch {
32 using DomainLocal =
typename Basis::Traits::DomainLocal;
34 using RangeField =
typename Basis::Traits::RangeField;
36 static const std::size_t dimRange = Basis::Traits::dimRange;
39 template<
typename Geometry>
40 static void jacobian(
const Basis& basis,
const Geometry& geometry,
41 const DomainLocal& xl,
42 std::vector<FieldMatrix<RangeField, dimRange,
45 jac.resize(basis.size());
46 basis.evaluateJacobian(xl, jac);
52 struct VectorBasisInterfaceSwitch<
53 Basis, typename
std::enable_if<
55 std::integral_constant<
57 Basis::Traits::dimDomain
66 using RangeField =
typename Basis::Traits::RangeFieldType;
68 static const std::size_t dimRange = Basis::Traits::dimRange;
71 template<
typename Geometry>
77 jac.resize(basis.size());
81 basis.evaluateJacobian(xl, ljac);
86 for(std::size_t i = 0; i < basis.size(); ++i)
87 for(std::size_t row=0; row < dimRange; ++row)
88 geojac.mv(ljac[i][row], jac[i][row]);
99 template<
typename PRM>
106 using RF =
typename PRM::Traits::RangeField;
109 using Real =
typename InstatBase::RealType;
111 static const bool navier = PRM::assemble_navier;
112 static const bool full_tensor = PRM::assemble_full_tensor;
117 enum { doPatternVolume =
true };
118 enum { doPatternSkeleton =
true };
121 enum { doSkeletonTwoSided =
false };
124 enum { doAlphaVolume =
true };
125 enum { doAlphaSkeleton =
true };
126 enum { doAlphaBoundary =
true };
127 enum { doLambdaVolume =
true };
142 prm(_prm), superintegration_order(_superintegration_order),
147 void preStep (Real , Real dt,
int )
160 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
161 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
163 const unsigned int dim = EG::Geometry::mydimension;
167 using LFSV_V = TypeTree::Child<LFSV,_0>;
168 const auto& lfsv_v =
child(lfsv,
_0);
169 const auto& lfsu_v =
child(lfsu,
_0);
171 using LFSV_P = TypeTree::Child<LFSV,_1>;
172 const auto& lfsv_p =
child(lfsv,
_1);
173 const auto& lfsu_p =
child(lfsu,
_1);
176 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
177 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
178 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
179 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
180 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
181 using RF =
typename BasisSwitch_V::RangeField;
182 using Range_V =
typename BasisSwitch_V::Range;
183 using Range_P =
typename BasisSwitch_P::Range;
184 using size_type =
typename LFSV::Traits::SizeType;
187 auto geo = eg.geometry();
190 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
191 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
192 const int jac_order = geo.type().isSimplex() ? 0 : 1;
193 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
195 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
198 for (
const auto& ip : quadratureRule(geo,qorder))
200 auto local = ip.position();
201 auto mu = prm.mu(eg,local);
202 auto rho = prm.rho(eg,local);
205 std::vector<Range_V> phi_v(lfsv_v.size());
208 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
209 for (size_type i=0; i<lfsu_v.size(); i++)
210 val_u.axpy(x(lfsu_v,i),phi_v[i]);
214 std::vector<Range_P> phi_p(lfsv_p.size());
215 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
219 for (size_type i=0; i<lfsu_p.size(); i++)
220 val_p.axpy(x(lfsu_p,i),phi_p[i]);
223 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
224 VectorBasisSwitch_V::jacobian
225 (FESwitch_V::basis(lfsv_v.finiteElement()), geo, local, jac_phi_v);
228 std::vector<RF> div_phi_v(lfsv_v.size(),0.0);
229 for (size_type i=0; i<lfsv_v.size(); i++)
230 for (size_type d=0; d<dim; d++)
231 div_phi_v[i] += jac_phi_v[i][d][d];
236 for (size_type i=0; i<lfsu_v.size(); i++){
237 jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
238 div_u += x(lfsu_v,i) * div_phi_v[i];
241 auto detj = geo.integrationElement(ip.position());
242 auto weight = ip.weight() * detj;
244 for (size_type i=0; i<lfsv_v.size(); i++) {
249 contraction(jac_u,jac_phi_v[i],dvdu);
250 r.accumulate(lfsv_v, i, dvdu * mu * weight);
255 r.accumulate(lfsv_v, i, - div_phi_v[i] * val_p * weight);
262 Range_V nabla_u_u(0.0);
263 jac_u.mv(val_u,nabla_u_u);
264 r.accumulate(lfsv_v, i, rho * (nabla_u_u*phi_v[i]) * weight);
269 for (size_type i=0; i<lfsv_p.size(); i++) {
273 r.accumulate(lfsv_p, i, - div_u * phi_p[i] * incomp_scaling * weight);
280 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
281 typename LocalMatrix>
282 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
283 LocalMatrix& mat)
const
285 const unsigned int dim = EG::Geometry::mydimension;
289 using LFSV_V = TypeTree::Child<LFSV,_0>;
290 const auto& lfsv_v =
child(lfsv,
_0);
291 const auto& lfsu_v =
child(lfsu,
_0);
293 using LFSV_P = TypeTree::Child<LFSV,_1>;
294 const auto& lfsv_p =
child(lfsv,
_1);
295 const auto& lfsu_p =
child(lfsu,
_1);
298 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
299 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
300 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
301 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
302 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
303 using RF =
typename BasisSwitch_V::RangeField;
304 using Range_V =
typename BasisSwitch_V::Range;
305 using Range_P =
typename BasisSwitch_P::Range;
306 using size_type =
typename LFSV::Traits::SizeType;
309 auto geo = eg.geometry();
312 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
313 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
314 const int jac_order = geo.type().isSimplex() ? 0 : 1;
315 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
317 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
320 for (
const auto& ip : quadratureRule(geo,qorder))
322 auto local = ip.position();
323 auto mu = prm.mu(eg,local);
324 auto rho = prm.rho(eg,local);
327 std::vector<Range_V> phi_v(lfsv_v.size());
330 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
331 for (size_type i=0; i<lfsu_v.size(); i++)
332 val_u.axpy(x(lfsu_v,i),phi_v[i]);
336 std::vector<Range_P> phi_p(lfsv_p.size());
337 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
340 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
341 VectorBasisSwitch_V::jacobian
342 (FESwitch_V::basis(lfsv_v.finiteElement()), geo, local, jac_phi_v);
344 assert(lfsu_v.size() == lfsv_v.size());
346 std::vector<RF> div_phi_v(lfsv_v.size(),0.0);
347 for (size_type i=0; i<lfsv_v.size(); i++)
348 for (size_type d=0; d<dim; d++)
349 div_phi_v[i] += jac_phi_v[i][d][d];
354 for (size_type i=0; i<lfsu_v.size(); i++){
355 jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
359 auto detj = geo.integrationElement(ip.position());
360 auto weight = ip.weight() * detj;
362 for(size_type i=0; i<lfsv_v.size(); i++) {
364 for(size_type j=0; j<lfsu_v.size(); j++) {
369 contraction(jac_phi_v[j],jac_phi_v[i],dvdu);
370 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * dvdu * weight);
377 Range_V nabla_u_phi_v_j(0.0);
378 jac_u.mv(phi_v[j],nabla_u_phi_v_j);
380 Range_V nabla_phi_v_j_u(0.0);
381 jac_phi_v[j].mv(val_u,nabla_phi_v_j_u);
382 mat.accumulate(lfsv_v,i,lfsu_v,j, rho * ((nabla_u_phi_v_j*phi_v[i]) + (nabla_phi_v_j_u*phi_v[i])) * weight);
387 for(size_type j=0; j<lfsv_p.size(); j++) {
392 mat.accumulate(lfsv_v,i,lfsu_p,j, -phi_p[j] * div_phi_v[i] * weight);
393 mat.accumulate(lfsv_p,j,lfsu_v,i, -phi_p[j] * div_phi_v[i] * incomp_scaling * weight);
401 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
402 void alpha_skeleton (
const IG& ig,
403 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
404 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
405 R& r_s, R& r_n)
const
408 const unsigned int dim = IG::Entity::dimension;
409 const unsigned int dimw = IG::coorddimension;
413 using LFSV_V = TypeTree::Child<LFSV,_0>;
414 const auto& lfsv_v_s =
child(lfsv_s,
_0);
415 const auto& lfsu_v_s =
child(lfsu_s,
_0);
416 const auto& lfsv_v_n =
child(lfsv_n,
_0);
417 const auto& lfsu_v_n =
child(lfsu_n,
_0);
419 using LFSV_P = TypeTree::Child<LFSV,_1>;
420 const auto& lfsv_p_s =
child(lfsv_s,
_1);
421 const auto& lfsu_p_s =
child(lfsu_s,
_1);
422 const auto& lfsv_p_n =
child(lfsv_n,
_1);
423 const auto& lfsu_p_n =
child(lfsu_n,
_1);
426 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
427 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
428 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
429 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
430 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
431 using DF =
typename BasisSwitch_V::DomainField;
432 using RF =
typename BasisSwitch_V::RangeField;
433 using Range_V =
typename BasisSwitch_V::Range;
434 using Range_P =
typename BasisSwitch_P::Range;
435 using size_type =
typename LFSV::Traits::SizeType;
438 const auto& cell_inside = ig.inside();
439 const auto& cell_outside = ig.outside();
442 auto geo = ig.geometry();
443 auto geo_inside = cell_inside.geometry();
444 auto geo_outside = cell_outside.geometry();
447 auto geo_in_inside = ig.geometryInInside();
448 auto geo_in_outside = ig.geometryInOutside();
451 const int v_order = FESwitch_V::basis(lfsv_v_s.finiteElement()).order();
452 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
453 const int qorder = 2*v_order + det_jac_order + superintegration_order;
455 const int epsilon = prm.epsilonIPSymmetryFactor();
456 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
458 auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
461 for (
const auto& ip : quadratureRule(geo,qorder))
465 auto local_s = geo_in_inside.global(ip.position());
466 auto local_n = geo_in_outside.global(ip.position());
469 std::vector<Range_V> phi_v_s(lfsv_v_s.size());
470 std::vector<Range_V> phi_v_n(lfsv_v_n.size());
471 FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
472 FESwitch_V::basis(lfsv_v_n.finiteElement()).evaluateFunction(local_n,phi_v_n);
475 std::vector<Range_P> phi_p_s(lfsv_p_s.size());
476 std::vector<Range_P> phi_p_n(lfsv_p_n.size());
477 FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
478 FESwitch_P::basis(lfsv_p_n.finiteElement()).evaluateFunction(local_n,phi_p_n);
481 Range_P val_p_s(0.0);
482 Range_P val_p_n(0.0);
483 for (size_type i=0; i<lfsu_p_s.size(); i++)
484 val_p_s.axpy(x_s(lfsu_p_s,i),phi_p_s[i]);
485 for (size_type i=0; i<lfsu_p_n.size(); i++)
486 val_p_n.axpy(x_n(lfsu_p_n,i),phi_p_n[i]);
489 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
490 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_n(lfsu_v_n.size());
491 VectorBasisSwitch_V::jacobian
492 (FESwitch_V::basis(lfsv_v_s.finiteElement()), geo_inside, local_s, jac_phi_v_s);
493 VectorBasisSwitch_V::jacobian
494 (FESwitch_V::basis(lfsv_v_n.finiteElement()), geo_outside, local_n, jac_phi_v_n);
497 Range_V val_u_s(0.0);
498 Range_V val_u_n(0.0);
501 for (size_type i=0; i<lfsu_v_s.size(); i++){
502 val_u_s.axpy(x_s(lfsu_v_s,i),phi_v_s[i]);
503 jac_u_s.axpy(x_s(lfsu_v_s,i),jac_phi_v_s[i]);
505 for (size_type i=0; i<lfsu_v_n.size(); i++){
506 val_u_n.axpy(x_n(lfsu_v_n,i),phi_v_n[i]);
507 jac_u_n.axpy(x_n(lfsu_v_n,i),jac_phi_v_n[i]);
510 auto normal = ig.unitOuterNormal(ip.position());
511 auto weight = ip.weight()*geo.integrationElement(ip.position());
512 auto mu = prm.mu(ig,ip.position());
514 auto factor = mu * weight;
517 auto jump = val_u_s - val_u_n;
520 auto mean_p = 0.5*(val_p_s + val_p_n);
524 add_compute_flux(jac_u_s,normal,flux_jac_u);
525 add_compute_flux(jac_u_n,normal,flux_jac_u);
529 for (
size_t i=0; i<lfsv_v_s.size(); i++) {
533 r_s.accumulate(lfsv_v_s, i, -(flux_jac_u * phi_v_s[i]) * factor);
539 add_compute_flux(jac_phi_v_s[i],normal,flux_jac_phi);
540 r_s.accumulate(lfsv_v_s, i, epsilon * 0.5 * (flux_jac_phi * jump) * factor);
545 r_s.accumulate(lfsv_v_s,i, penalty_factor * (jump*phi_v_s[i]) * factor);
550 r_s.accumulate(lfsv_v_s,i, mean_p * (phi_v_s[i]*normal) * weight);
554 for (
size_t i=0; i<lfsv_v_n.size(); i++) {
558 r_n.accumulate(lfsv_v_n, i, (flux_jac_u * phi_v_n[i]) * factor);
564 add_compute_flux(jac_phi_v_n[i],normal,flux_jac_phi);
565 r_n.accumulate(lfsv_v_n, i, epsilon * 0.5 * (flux_jac_phi * jump) * factor);
570 r_n.accumulate(lfsv_v_n,i, -penalty_factor * (jump*phi_v_n[i]) * factor);
575 r_n.accumulate(lfsv_v_n,i, -mean_p * (phi_v_n[i]*normal) * weight);
581 for (
size_t i=0; i<lfsv_p_s.size(); i++)
582 r_s.accumulate(lfsv_p_s,i, 0.5*phi_p_s[i] * (jump*normal) * incomp_scaling * weight);
583 for (
size_t i=0; i<lfsv_p_n.size(); i++)
584 r_n.accumulate(lfsv_p_n,i, 0.5*phi_p_n[i] * (jump*normal) * incomp_scaling * weight);
590 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
591 typename LocalMatrix>
592 void jacobian_skeleton (
const IG& ig,
593 const LFSU& lfsu_s,
const X&,
const LFSV& lfsv_s,
594 const LFSU& lfsu_n,
const X&,
const LFSV& lfsv_n,
595 LocalMatrix& mat_ss, LocalMatrix& mat_sn,
596 LocalMatrix& mat_ns, LocalMatrix& mat_nn)
const
599 const unsigned int dim = IG::Entity::dimension;
600 const unsigned int dimw = IG::coorddimension;
604 using LFSV_V = TypeTree::Child<LFSV,_0>;
605 const auto& lfsv_v_s =
child(lfsv_s,
_0);
606 const auto& lfsu_v_s =
child(lfsu_s,
_0);
607 const auto& lfsv_v_n =
child(lfsv_n,
_0);
608 const auto& lfsu_v_n =
child(lfsu_n,
_0);
610 using LFSV_P = TypeTree::Child<LFSV,_1>;
611 const auto& lfsv_p_s =
child(lfsv_s,
_1);
612 const auto& lfsu_p_s =
child(lfsu_s,
_1);
613 const auto& lfsv_p_n =
child(lfsv_n,
_1);
614 const auto& lfsu_p_n =
child(lfsu_n,
_1);
617 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
618 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
619 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
620 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
621 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
622 using DF =
typename BasisSwitch_V::DomainField;
623 using RF =
typename BasisSwitch_V::RangeField;
624 using Range_V =
typename BasisSwitch_V::Range;
625 using Range_P =
typename BasisSwitch_P::Range;
626 using size_type =
typename LFSV::Traits::SizeType;
629 auto const& cell_inside = ig.inside();
630 auto const& cell_outside = ig.outside();
633 auto geo = ig.geometry();
634 auto geo_inside = cell_inside.geometry();
635 auto geo_outside = cell_outside.geometry();
638 auto geo_in_inside = ig.geometryInInside();
639 auto geo_in_outside = ig.geometryInOutside();
642 const int v_order = FESwitch_V::basis(lfsv_v_s.finiteElement()).order();
643 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
644 const int qorder = 2*v_order + det_jac_order + superintegration_order;
646 auto epsilon = prm.epsilonIPSymmetryFactor();
647 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
649 auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
652 for (
const auto& ip : quadratureRule(geo,qorder))
656 auto local_s = geo_in_inside.global(ip.position());
657 auto local_n = geo_in_outside.global(ip.position());
660 std::vector<Range_V> phi_v_s(lfsv_v_s.size());
661 std::vector<Range_V> phi_v_n(lfsv_v_n.size());
662 FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
663 FESwitch_V::basis(lfsv_v_n.finiteElement()).evaluateFunction(local_n,phi_v_n);
666 std::vector<Range_P> phi_p_s(lfsv_p_s.size());
667 std::vector<Range_P> phi_p_n(lfsv_p_n.size());
668 FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
669 FESwitch_P::basis(lfsv_p_n.finiteElement()).evaluateFunction(local_n,phi_p_n);
672 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
673 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_n(lfsu_v_n.size());
674 VectorBasisSwitch_V::jacobian
675 (FESwitch_V::basis(lfsv_v_s.finiteElement()), geo_inside, local_s, jac_phi_v_s);
676 VectorBasisSwitch_V::jacobian
677 (FESwitch_V::basis(lfsv_v_n.finiteElement()), geo_outside, local_n, jac_phi_v_n);
679 auto normal = ig.unitOuterNormal(ip.position());
680 auto weight = ip.weight()*geo.integrationElement(ip.position());
681 auto mu = prm.mu(ig,ip.position());
683 auto factor = mu * weight;
688 for(size_type i=0; i<lfsv_v_s.size(); i++) {
692 add_compute_flux(jac_phi_v_s[i],normal,flux_jac_phi_i);
699 for(size_type j=0; j<lfsu_v_s.size(); j++) {
701 add_compute_flux(jac_phi_v_s[j],normal,flux_jac_phi_j);
703 mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, -0.5 * (flux_jac_phi_j*phi_v_s[i]) * factor);
704 mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, epsilon * 0.5 * (flux_jac_phi_i*phi_v_s[j]) * factor);
705 mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, penalty_factor * (phi_v_s[j]*phi_v_s[i]) * factor);
708 for(size_type j=0; j<lfsu_v_n.size(); j++) {
710 add_compute_flux(jac_phi_v_n[j],normal,flux_jac_phi_j);
712 mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -0.5 * (flux_jac_phi_j*phi_v_s[i]) * factor);
713 mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -epsilon * 0.5 * (flux_jac_phi_i*phi_v_n[j]) * factor);
714 mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -penalty_factor * (phi_v_n[j]*phi_v_s[i]) * factor);
720 for(size_type j=0; j<lfsu_p_s.size(); j++) {
721 mat_ss.accumulate(lfsv_v_s,i,lfsu_p_s,j, 0.5*phi_p_s[j] * (phi_v_s[i]*normal) * weight);
724 for(size_type j=0; j<lfsu_p_n.size(); j++) {
725 mat_sn.accumulate(lfsv_v_s,i,lfsu_p_n,j, 0.5*phi_p_n[j] * (phi_v_s[i]*normal) * weight);
732 for(size_type i=0; i<lfsv_v_n.size(); i++) {
736 add_compute_flux(jac_phi_v_n[i],normal,flux_jac_phi_i);
743 for(size_type j=0; j<lfsu_v_s.size(); j++) {
745 add_compute_flux(jac_phi_v_s[j],normal,flux_jac_phi_j);
747 mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, 0.5 * (flux_jac_phi_j*phi_v_n[i]) * factor);
748 mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, epsilon * 0.5 * (flux_jac_phi_i*phi_v_s[j]) * factor);
749 mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, -penalty_factor * (phi_v_s[j]*phi_v_n[i]) * factor);
752 for(size_type j=0; j<lfsu_v_n.size(); j++) {
754 add_compute_flux(jac_phi_v_n[j],normal,flux_jac_phi_j);
756 mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, 0.5 * (flux_jac_phi_j*phi_v_n[i]) * factor);
757 mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, -epsilon * 0.5 * (flux_jac_phi_i*phi_v_n[j]) * factor);
758 mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, penalty_factor * (phi_v_n[j]*phi_v_n[i]) * factor);
764 for(size_type j=0; j<lfsu_p_s.size(); j++) {
765 mat_ns.accumulate(lfsv_v_n,i,lfsu_p_s,j, -0.5*phi_p_s[j] * (phi_v_n[i]*normal) * weight);
768 for(size_type j=0; j<lfsu_p_n.size(); j++) {
769 mat_nn.accumulate(lfsv_v_n,i,lfsu_p_n,j, -0.5*phi_p_n[j] * (phi_v_n[i]*normal) * weight);
776 for(size_type i=0; i<lfsv_p_s.size(); i++) {
777 for(size_type j=0; j<lfsu_v_s.size(); j++)
778 mat_ss.accumulate(lfsv_p_s,i,lfsu_v_s,j, 0.5*phi_p_s[i] * (phi_v_s[j]*normal) * incomp_scaling * weight);
780 for(size_type j=0; j<lfsu_v_n.size(); j++)
781 mat_sn.accumulate(lfsv_p_s,i,lfsu_v_n,j, -0.5*phi_p_s[i] * (phi_v_n[j]*normal) * incomp_scaling * weight);
784 for(size_type i=0; i<lfsv_p_n.size(); i++) {
785 for(size_type j=0; j<lfsu_v_s.size(); j++)
786 mat_ns.accumulate(lfsv_p_n,i,lfsu_v_s,j, 0.5*phi_p_n[i] * (phi_v_s[j]*normal) * incomp_scaling * weight);
788 for(size_type j=0; j<lfsu_v_n.size(); j++)
789 mat_nn.accumulate(lfsv_p_n,i,lfsu_v_n,j, -0.5*phi_p_n[i] * (phi_v_n[j]*normal) * incomp_scaling * weight);
796 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
797 void alpha_boundary (
const IG& ig,
798 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
802 const unsigned int dim = IG::Entity::dimension;
803 const unsigned int dimw = IG::coorddimension;
807 using LFSV_V = TypeTree::Child<LFSV,_0>;
808 const auto& lfsv_v =
child(lfsv,
_0);
809 const auto& lfsu_v =
child(lfsu,
_0);
811 using LFSV_P = TypeTree::Child<LFSV,_1>;
812 const auto& lfsv_p =
child(lfsv,
_1);
813 const auto& lfsu_p =
child(lfsu,
_1);
816 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
817 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
818 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
819 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
820 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
821 using DF =
typename BasisSwitch_V::DomainField;
822 using RF =
typename BasisSwitch_V::RangeField;
823 using Range_V =
typename BasisSwitch_V::Range;
824 using Range_P =
typename BasisSwitch_P::Range;
825 using size_type =
typename LFSV::Traits::SizeType;
828 const auto& cell_inside = ig.inside();
831 auto geo = ig.geometry();
832 auto geo_inside = cell_inside.geometry();
835 auto geo_in_inside = ig.geometryInInside();
838 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
839 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
840 const int qorder = 2*v_order + det_jac_order + superintegration_order;
842 auto epsilon = prm.epsilonIPSymmetryFactor();
843 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
845 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
848 for (
const auto& ip : quadratureRule(geo,qorder))
851 auto local = geo_in_inside.global(ip.position());
854 std::vector<Range_V> phi_v(lfsv_v.size());
855 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
858 std::vector<Range_P> phi_p(lfsv_p.size());
859 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
862 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
863 VectorBasisSwitch_V::jacobian
864 (FESwitch_V::basis(lfsv_v.finiteElement()), geo_inside, local, jac_phi_v);
868 for (size_type i=0; i<lfsu_p.size(); i++)
869 val_p.axpy(x(lfsu_p,i),phi_p[i]);
874 for (size_type i=0; i<lfsu_v.size(); i++){
875 val_u.axpy(x(lfsu_v,i),phi_v[i]);
876 jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
879 auto normal = ig.unitOuterNormal(ip.position());
880 auto weight = ip.weight()*geo.integrationElement(ip.position());
881 auto mu = prm.mu(ig,ip.position());
884 auto bctype(prm.bctype(ig,ip.position()));
886 if (bctype == BC::VelocityDirichlet) {
888 auto u0(prm.g(cell_inside,local));
889 auto jump = val_u - u0;
893 add_compute_flux(jac_u,normal,flux_jac_u);
895 for (
size_t i=0; i<lfsv_v.size(); i++) {
899 r.accumulate(lfsv_v,i, -mu * (flux_jac_u * phi_v[i]) * weight);
905 add_compute_flux(jac_phi_v[i],normal,flux_jac_phi);
906 r.accumulate(lfsv_v,i, mu * epsilon * (flux_jac_phi*jump) * weight);
911 r.accumulate(lfsv_v,i, mu * (jump*phi_v[i]) * penalty_factor * weight);
916 r.accumulate(lfsv_v,i, val_p * (phi_v[i]*normal) * weight);
922 for(size_type i=0; i<lfsv_p.size(); i++) {
923 r.accumulate(lfsv_p,i, phi_p[i] * (jump*normal) * incomp_scaling * weight);
927 if (bctype == BC::StressNeumann) {
928 auto stress(prm.j(ig,ip.position(),normal));
930 for(size_type i=0; i<lfsv_v.size(); i++) {
931 r.accumulate(lfsv_v,i, (stress*phi_v[i]) * weight);
939 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
940 typename LocalMatrix>
941 void jacobian_boundary (
const IG& ig,
942 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
943 LocalMatrix& mat)
const
946 const unsigned int dim = IG::Entity::dimension;
947 const unsigned int dimw = IG::coorddimension;
951 using LFSV_V = TypeTree::Child<LFSV,_0>;
952 const auto& lfsv_v =
child(lfsv,
_0);
953 const auto& lfsu_v =
child(lfsu,
_0);
955 using LFSV_P = TypeTree::Child<LFSV,_1>;
956 const auto& lfsv_p =
child(lfsv,
_1);
957 const auto& lfsu_p =
child(lfsu,
_1);
960 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
961 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
962 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
963 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
964 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
965 using DF =
typename BasisSwitch_V::DomainField;
966 using RF =
typename BasisSwitch_V::RangeField;
967 using Range_V =
typename BasisSwitch_V::Range;
968 using Range_P =
typename BasisSwitch_P::Range;
969 using size_type =
typename LFSV::Traits::SizeType;
972 const auto& cell_inside = ig.inside();
975 auto geo = ig.geometry();
976 auto geo_inside = cell_inside.geometry();
979 auto geo_in_inside = ig.geometryInInside();
982 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
983 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
984 const int qorder = 2*v_order + det_jac_order + superintegration_order;
986 auto epsilon = prm.epsilonIPSymmetryFactor();
987 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
989 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
992 for (
const auto& ip : quadratureRule(geo,qorder))
995 auto local = geo_in_inside.global(ip.position());
998 std::vector<Range_V> phi_v(lfsv_v.size());
999 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1002 std::vector<Range_P> phi_p(lfsv_p.size());
1003 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1006 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
1007 VectorBasisSwitch_V::jacobian
1008 (FESwitch_V::basis(lfsv_v.finiteElement()), geo_inside, local, jac_phi_v);
1010 auto normal = ig.unitOuterNormal(ip.position());
1011 auto weight = ip.weight()*geo.integrationElement(ip.position());
1012 auto mu = prm.mu(ig,ip.position());
1015 auto bctype(prm.bctype(ig,ip.position()));
1017 if (bctype == BC::VelocityDirichlet) {
1019 for(size_type i=0; i<lfsv_v.size(); i++) {
1022 add_compute_flux(jac_phi_v[i],normal,flux_jac_phi_i);
1024 for(size_type j=0; j<lfsu_v.size(); j++) {
1030 add_compute_flux(jac_phi_v[j],normal,flux_jac_phi_j);
1032 mat.accumulate(lfsv_v,i,lfsu_v,j, -mu * (flux_jac_phi_j*phi_v[i]) * weight);
1033 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * epsilon * (flux_jac_phi_i*phi_v[j]) *weight);
1038 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (phi_v[j]*phi_v[i]) * penalty_factor * weight);
1044 for(size_type j=0; j<lfsu_p.size(); j++) {
1045 mat.accumulate(lfsv_v,i,lfsu_p,j, phi_p[j] * (phi_v[i]*normal) * weight);
1052 for(size_type i=0; i<lfsv_p.size(); i++) {
1053 for(size_type j=0; j<lfsu_v.size(); j++) {
1054 mat.accumulate(lfsv_p,i,lfsu_v,j, phi_p[i] * (phi_v[j]*normal) * incomp_scaling * weight);
1065 template<
typename EG,
typename LFSV,
typename R>
1066 void lambda_volume (
const EG& eg,
const LFSV& lfsv, R& r)
const
1068 const unsigned int dim = EG::Geometry::mydimension;
1072 using LFSV_V = TypeTree::Child<LFSV,_0>;
1073 const auto& lfsv_v =
child(lfsv,
_0);
1076 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
1077 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
1078 using Range_V =
typename BasisSwitch_V::Range;
1079 using size_type =
typename LFSV::Traits::SizeType;
1082 const auto& cell = eg.entity();
1085 auto geo = eg.geometry();
1088 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1089 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
1090 const int qorder = 2*v_order + det_jac_order + superintegration_order;
1093 for (
const auto& ip : quadratureRule(geo,qorder))
1095 auto local = ip.position();
1099 std::vector<Range_V> phi_v(lfsv_v.size());
1100 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1102 auto weight = ip.weight() * geo.integrationElement(ip.position());
1105 auto fval(prm.f(cell,local));
1110 for(size_type i=0; i<lfsv_v.size(); i++)
1111 r.accumulate(lfsv_v,i, -(fval*phi_v[i]) * weight);
1118 template<
class M,
class RF>
1119 static void contraction(
const M & a,
const M & b, RF & v)
1122 const int an = a.N();
1123 const int am = a.M();
1124 for(
int r=0; r<an; ++r)
1125 for(
int c=0; c<am; ++c)
1126 v += a[r][c] * b[r][c];
1129 template<
class DU,
class R>
1130 static void add_compute_flux(
const DU & du,
const R & n, R & result)
1132 const int N = du.N();
1133 const int M = du.M();
1134 for(
int r=0; r<N; ++r)
1135 for(
int c=0; c<M; ++c)
1136 result[r] += du[r][c] * n[c];
1140 const int superintegration_order;
A dense n x m matrix.
Definition: fmatrix.hh:69
vector space out of a tensor product of fields.
Definition: fvector.hh:96
Wrapper class for geometries.
Definition: geometry.hh:67
Implementation::JacobianInverseTransposed JacobianInverseTransposed
type of jacobian inverse transposed
Definition: geometry.hh:115
@ coorddimension
Definition: geometry.hh:92
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
Return inverse of transposed of Jacobian.
Definition: geometry.hh:265
A local operator for solving the Navier-Stokes equations using a DG discretization and a vector-value...
Definition: dgnavierstokesvelvecfem.hh:104
DGNavierStokesVelVecFEM(PRM &_prm, int _superintegration_order=0)
Constructor.
Definition: dgnavierstokesvelvecfem.hh:141
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:51
constexpr index_constant< 1 > _1
Compile time index with value 1.
Definition: indices.hh:54
Namespace with predefined compile time indices for the range [0,19].
Definition: indices.hh:49
Dune namespace.
Definition: alignedallocator.hh:14
A hierarchical structure of string parameters.
Definition: stokesparameter.hh:32
static void jacobian(const Basis &basis, const Geometry &geometry, const DomainLocal &xl, std::vector< FieldMatrix< RangeField, dimRange, Geometry::coorddimension > > &jac)
Compute global jacobian matrix for vector valued bases.
Definition: dgnavierstokesvelvecfem.hh:72
typename Basis::Traits::DomainType DomainLocal
export vector type of the local coordinates
Definition: dgnavierstokesvelvecfem.hh:64
typename Basis::Traits::RangeFieldType RangeField
export field type of the values
Definition: dgnavierstokesvelvecfem.hh:66