4#ifndef DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESVELVECFEM_HH
5#define DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESVELVECFEM_HH
12#include <dune/localfunctions/common/interfaceswitch.hh>
13#include <dune/pdelab/localoperator/idefault.hh>
15#include <dune/pdelab/common/quadraturerules.hh>
16#include <dune/pdelab/localoperator/defaultimp.hh>
17#include <dune/pdelab/localoperator/pattern.hh>
18#include <dune/pdelab/localoperator/flags.hh>
19#include <dune/pdelab/localoperator/dgnavierstokesparameter.hh>
20#include <dune/pdelab/localoperator/navierstokesmass.hh>
25#define PBLOCK (- VBLOCK + 1)
30 template<
class Basis,
class Dummy =
void>
31 struct VectorBasisInterfaceSwitch {
33 using DomainLocal =
typename Basis::Traits::DomainLocal;
35 using RangeField =
typename Basis::Traits::RangeField;
37 static const std::size_t dimRange = Basis::Traits::dimRange;
40 template<
typename Geometry>
41 static void jacobian(
const Basis& basis,
const Geometry& geometry,
42 const DomainLocal& xl,
43 std::vector<FieldMatrix<RangeField, dimRange,
46 jac.resize(basis.size());
47 basis.evaluateJacobian(xl, jac);
53 struct VectorBasisInterfaceSwitch<
54 Basis, typename
std::enable_if<
56 std::integral_constant<
58 Basis::Traits::dimDomain
67 using RangeField =
typename Basis::Traits::RangeFieldType;
69 static const std::size_t dimRange = Basis::Traits::dimRange;
72 template<
typename Geometry>
78 jac.resize(basis.size());
82 basis.evaluateJacobian(xl, ljac);
87 for(std::size_t i = 0; i < basis.size(); ++i)
88 for(std::size_t row=0; row < dimRange; ++row)
89 geojac.mv(ljac[i][row], jac[i][row]);
100 template<
typename PRM>
107 using RF =
typename PRM::Traits::RangeField;
110 using Real =
typename InstatBase::RealType;
112 static const bool navier = PRM::assemble_navier;
113 static const bool full_tensor = PRM::assemble_full_tensor;
118 enum { doPatternVolume =
true };
119 enum { doPatternSkeleton =
true };
122 enum { doSkeletonTwoSided =
false };
125 enum { doAlphaVolume =
true };
126 enum { doAlphaSkeleton =
true };
127 enum { doAlphaBoundary =
true };
128 enum { doLambdaVolume =
true };
143 prm(_prm), superintegration_order(_superintegration_order),
148 void preStep (Real , Real dt,
int )
161 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
162 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
164 const unsigned int dim = EG::Geometry::mydimension;
168 using LFSV_V = TypeTree::Child<LFSV,_0>;
169 const auto& lfsv_v =
child(lfsv,
_0);
170 const auto& lfsu_v =
child(lfsu,
_0);
172 using LFSV_P = TypeTree::Child<LFSV,_1>;
173 const auto& lfsv_p =
child(lfsv,
_1);
174 const auto& lfsu_p =
child(lfsu,
_1);
177 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
178 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
179 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
180 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
181 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
182 using RF =
typename BasisSwitch_V::RangeField;
183 using Range_V =
typename BasisSwitch_V::Range;
184 using Range_P =
typename BasisSwitch_P::Range;
185 using size_type =
typename LFSV::Traits::SizeType;
188 auto geo = eg.geometry();
191 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
192 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
193 const int jac_order = geo.type().isSimplex() ? 0 : 1;
194 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
196 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
199 for (
const auto& ip : quadratureRule(geo,qorder))
201 auto local = ip.position();
202 auto mu = prm.mu(eg,local);
203 auto rho = prm.rho(eg,local);
206 std::vector<Range_V> phi_v(lfsv_v.size());
209 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
210 for (size_type i=0; i<lfsu_v.size(); i++)
211 val_u.axpy(x(lfsu_v,i),phi_v[i]);
215 std::vector<Range_P> phi_p(lfsv_p.size());
216 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
220 for (size_type i=0; i<lfsu_p.size(); i++)
221 val_p.axpy(x(lfsu_p,i),phi_p[i]);
224 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
225 VectorBasisSwitch_V::jacobian
226 (FESwitch_V::basis(lfsv_v.finiteElement()), geo, local, jac_phi_v);
229 std::vector<RF> div_phi_v(lfsv_v.size(),0.0);
230 for (size_type i=0; i<lfsv_v.size(); i++)
231 for (size_type d=0; d<dim; d++)
232 div_phi_v[i] += jac_phi_v[i][d][d];
237 for (size_type i=0; i<lfsu_v.size(); i++){
238 jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
239 div_u += x(lfsu_v,i) * div_phi_v[i];
242 auto detj = geo.integrationElement(ip.position());
243 auto weight = ip.weight() * detj;
245 for (size_type i=0; i<lfsv_v.size(); i++) {
250 contraction(jac_u,jac_phi_v[i],dvdu);
251 r.accumulate(lfsv_v, i, dvdu * mu * weight);
256 r.accumulate(lfsv_v, i, - div_phi_v[i] * val_p * weight);
263 Range_V nabla_u_u(0.0);
264 jac_u.mv(val_u,nabla_u_u);
265 r.accumulate(lfsv_v, i, rho * (nabla_u_u*phi_v[i]) * weight);
270 for (size_type i=0; i<lfsv_p.size(); i++) {
274 r.accumulate(lfsv_p, i, - div_u * phi_p[i] * incomp_scaling * weight);
281 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
282 typename LocalMatrix>
283 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
284 LocalMatrix& mat)
const
286 const unsigned int dim = EG::Geometry::mydimension;
290 using LFSV_V = TypeTree::Child<LFSV,_0>;
291 const auto& lfsv_v =
child(lfsv,
_0);
292 const auto& lfsu_v =
child(lfsu,
_0);
294 using LFSV_P = TypeTree::Child<LFSV,_1>;
295 const auto& lfsv_p =
child(lfsv,
_1);
296 const auto& lfsu_p =
child(lfsu,
_1);
299 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
300 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
301 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
302 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
303 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
304 using RF =
typename BasisSwitch_V::RangeField;
305 using Range_V =
typename BasisSwitch_V::Range;
306 using Range_P =
typename BasisSwitch_P::Range;
307 using size_type =
typename LFSV::Traits::SizeType;
310 auto geo = eg.geometry();
313 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
314 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
315 const int jac_order = geo.type().isSimplex() ? 0 : 1;
316 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
318 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
321 for (
const auto& ip : quadratureRule(geo,qorder))
323 auto local = ip.position();
324 auto mu = prm.mu(eg,local);
325 auto rho = prm.rho(eg,local);
328 std::vector<Range_V> phi_v(lfsv_v.size());
331 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
332 for (size_type i=0; i<lfsu_v.size(); i++)
333 val_u.axpy(x(lfsu_v,i),phi_v[i]);
337 std::vector<Range_P> phi_p(lfsv_p.size());
338 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
341 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
342 VectorBasisSwitch_V::jacobian
343 (FESwitch_V::basis(lfsv_v.finiteElement()), geo, local, jac_phi_v);
345 assert(lfsu_v.size() == lfsv_v.size());
347 std::vector<RF> div_phi_v(lfsv_v.size(),0.0);
348 for (size_type i=0; i<lfsv_v.size(); i++)
349 for (size_type d=0; d<dim; d++)
350 div_phi_v[i] += jac_phi_v[i][d][d];
355 for (size_type i=0; i<lfsu_v.size(); i++){
356 jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
360 auto detj = geo.integrationElement(ip.position());
361 auto weight = ip.weight() * detj;
363 for(size_type i=0; i<lfsv_v.size(); i++) {
365 for(size_type j=0; j<lfsu_v.size(); j++) {
370 contraction(jac_phi_v[j],jac_phi_v[i],dvdu);
371 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * dvdu * weight);
378 Range_V nabla_u_phi_v_j(0.0);
379 jac_u.mv(phi_v[j],nabla_u_phi_v_j);
381 Range_V nabla_phi_v_j_u(0.0);
382 jac_phi_v[j].mv(val_u,nabla_phi_v_j_u);
383 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);
388 for(size_type j=0; j<lfsv_p.size(); j++) {
393 mat.accumulate(lfsv_v,i,lfsu_p,j, -phi_p[j] * div_phi_v[i] * weight);
394 mat.accumulate(lfsv_p,j,lfsu_v,i, -phi_p[j] * div_phi_v[i] * incomp_scaling * weight);
402 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
403 void alpha_skeleton (
const IG& ig,
404 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
405 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
406 R& r_s, R& r_n)
const
409 const unsigned int dim = IG::Entity::dimension;
410 const unsigned int dimw = IG::coorddimension;
414 using LFSV_V = TypeTree::Child<LFSV,_0>;
415 const auto& lfsv_v_s =
child(lfsv_s,
_0);
416 const auto& lfsu_v_s =
child(lfsu_s,
_0);
417 const auto& lfsv_v_n =
child(lfsv_n,
_0);
418 const auto& lfsu_v_n =
child(lfsu_n,
_0);
420 using LFSV_P = TypeTree::Child<LFSV,_1>;
421 const auto& lfsv_p_s =
child(lfsv_s,
_1);
422 const auto& lfsu_p_s =
child(lfsu_s,
_1);
423 const auto& lfsv_p_n =
child(lfsv_n,
_1);
424 const auto& lfsu_p_n =
child(lfsu_n,
_1);
427 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
428 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
429 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
430 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
431 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
432 using DF =
typename BasisSwitch_V::DomainField;
433 using RF =
typename BasisSwitch_V::RangeField;
434 using Range_V =
typename BasisSwitch_V::Range;
435 using Range_P =
typename BasisSwitch_P::Range;
436 using size_type =
typename LFSV::Traits::SizeType;
439 const auto& cell_inside = ig.inside();
440 const auto& cell_outside = ig.outside();
443 auto geo = ig.geometry();
444 auto geo_inside = cell_inside.geometry();
445 auto geo_outside = cell_outside.geometry();
448 auto geo_in_inside = ig.geometryInInside();
449 auto geo_in_outside = ig.geometryInOutside();
452 const int v_order = FESwitch_V::basis(lfsv_v_s.finiteElement()).order();
453 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
454 const int qorder = 2*v_order + det_jac_order + superintegration_order;
456 const int epsilon = prm.epsilonIPSymmetryFactor();
457 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
459 auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
462 for (
const auto& ip : quadratureRule(geo,qorder))
466 auto local_s = geo_in_inside.global(ip.position());
467 auto local_n = geo_in_outside.global(ip.position());
470 std::vector<Range_V> phi_v_s(lfsv_v_s.size());
471 std::vector<Range_V> phi_v_n(lfsv_v_n.size());
472 FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
473 FESwitch_V::basis(lfsv_v_n.finiteElement()).evaluateFunction(local_n,phi_v_n);
476 std::vector<Range_P> phi_p_s(lfsv_p_s.size());
477 std::vector<Range_P> phi_p_n(lfsv_p_n.size());
478 FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
479 FESwitch_P::basis(lfsv_p_n.finiteElement()).evaluateFunction(local_n,phi_p_n);
482 Range_P val_p_s(0.0);
483 Range_P val_p_n(0.0);
484 for (size_type i=0; i<lfsu_p_s.size(); i++)
485 val_p_s.axpy(x_s(lfsu_p_s,i),phi_p_s[i]);
486 for (size_type i=0; i<lfsu_p_n.size(); i++)
487 val_p_n.axpy(x_n(lfsu_p_n,i),phi_p_n[i]);
490 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
491 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_n(lfsu_v_n.size());
492 VectorBasisSwitch_V::jacobian
493 (FESwitch_V::basis(lfsv_v_s.finiteElement()), geo_inside, local_s, jac_phi_v_s);
494 VectorBasisSwitch_V::jacobian
495 (FESwitch_V::basis(lfsv_v_n.finiteElement()), geo_outside, local_n, jac_phi_v_n);
498 Range_V val_u_s(0.0);
499 Range_V val_u_n(0.0);
502 for (size_type i=0; i<lfsu_v_s.size(); i++){
503 val_u_s.axpy(x_s(lfsu_v_s,i),phi_v_s[i]);
504 jac_u_s.axpy(x_s(lfsu_v_s,i),jac_phi_v_s[i]);
506 for (size_type i=0; i<lfsu_v_n.size(); i++){
507 val_u_n.axpy(x_n(lfsu_v_n,i),phi_v_n[i]);
508 jac_u_n.axpy(x_n(lfsu_v_n,i),jac_phi_v_n[i]);
511 auto normal = ig.unitOuterNormal(ip.position());
512 auto weight = ip.weight()*geo.integrationElement(ip.position());
513 auto mu = prm.mu(ig,ip.position());
515 auto factor = mu * weight;
518 auto jump = val_u_s - val_u_n;
521 auto mean_p = 0.5*(val_p_s + val_p_n);
525 add_compute_flux(jac_u_s,normal,flux_jac_u);
526 add_compute_flux(jac_u_n,normal,flux_jac_u);
530 for (
size_t i=0; i<lfsv_v_s.size(); i++) {
534 r_s.accumulate(lfsv_v_s, i, -(flux_jac_u * phi_v_s[i]) * factor);
540 add_compute_flux(jac_phi_v_s[i],normal,flux_jac_phi);
541 r_s.accumulate(lfsv_v_s, i, epsilon * 0.5 * (flux_jac_phi * jump) * factor);
546 r_s.accumulate(lfsv_v_s,i, penalty_factor * (jump*phi_v_s[i]) * factor);
551 r_s.accumulate(lfsv_v_s,i, mean_p * (phi_v_s[i]*normal) * weight);
555 for (
size_t i=0; i<lfsv_v_n.size(); i++) {
559 r_n.accumulate(lfsv_v_n, i, (flux_jac_u * phi_v_n[i]) * factor);
565 add_compute_flux(jac_phi_v_n[i],normal,flux_jac_phi);
566 r_n.accumulate(lfsv_v_n, i, epsilon * 0.5 * (flux_jac_phi * jump) * factor);
571 r_n.accumulate(lfsv_v_n,i, -penalty_factor * (jump*phi_v_n[i]) * factor);
576 r_n.accumulate(lfsv_v_n,i, -mean_p * (phi_v_n[i]*normal) * weight);
582 for (
size_t i=0; i<lfsv_p_s.size(); i++)
583 r_s.accumulate(lfsv_p_s,i, 0.5*phi_p_s[i] * (jump*normal) * incomp_scaling * weight);
584 for (
size_t i=0; i<lfsv_p_n.size(); i++)
585 r_n.accumulate(lfsv_p_n,i, 0.5*phi_p_n[i] * (jump*normal) * incomp_scaling * weight);
591 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
592 typename LocalMatrix>
593 void jacobian_skeleton (
const IG& ig,
594 const LFSU& lfsu_s,
const X&,
const LFSV& lfsv_s,
595 const LFSU& lfsu_n,
const X&,
const LFSV& lfsv_n,
596 LocalMatrix& mat_ss, LocalMatrix& mat_sn,
597 LocalMatrix& mat_ns, LocalMatrix& mat_nn)
const
600 const unsigned int dim = IG::Entity::dimension;
601 const unsigned int dimw = IG::coorddimension;
605 using LFSV_V = TypeTree::Child<LFSV,_0>;
606 const auto& lfsv_v_s =
child(lfsv_s,
_0);
607 const auto& lfsu_v_s =
child(lfsu_s,
_0);
608 const auto& lfsv_v_n =
child(lfsv_n,
_0);
609 const auto& lfsu_v_n =
child(lfsu_n,
_0);
611 using LFSV_P = TypeTree::Child<LFSV,_1>;
612 const auto& lfsv_p_s =
child(lfsv_s,
_1);
613 const auto& lfsu_p_s =
child(lfsu_s,
_1);
614 const auto& lfsv_p_n =
child(lfsv_n,
_1);
615 const auto& lfsu_p_n =
child(lfsu_n,
_1);
618 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
619 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
620 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
621 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
622 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
623 using DF =
typename BasisSwitch_V::DomainField;
624 using RF =
typename BasisSwitch_V::RangeField;
625 using Range_V =
typename BasisSwitch_V::Range;
626 using Range_P =
typename BasisSwitch_P::Range;
627 using size_type =
typename LFSV::Traits::SizeType;
630 auto const& cell_inside = ig.inside();
631 auto const& cell_outside = ig.outside();
634 auto geo = ig.geometry();
635 auto geo_inside = cell_inside.geometry();
636 auto geo_outside = cell_outside.geometry();
639 auto geo_in_inside = ig.geometryInInside();
640 auto geo_in_outside = ig.geometryInOutside();
643 const int v_order = FESwitch_V::basis(lfsv_v_s.finiteElement()).order();
644 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
645 const int qorder = 2*v_order + det_jac_order + superintegration_order;
647 auto epsilon = prm.epsilonIPSymmetryFactor();
648 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
650 auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
653 for (
const auto& ip : quadratureRule(geo,qorder))
657 auto local_s = geo_in_inside.global(ip.position());
658 auto local_n = geo_in_outside.global(ip.position());
661 std::vector<Range_V> phi_v_s(lfsv_v_s.size());
662 std::vector<Range_V> phi_v_n(lfsv_v_n.size());
663 FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
664 FESwitch_V::basis(lfsv_v_n.finiteElement()).evaluateFunction(local_n,phi_v_n);
667 std::vector<Range_P> phi_p_s(lfsv_p_s.size());
668 std::vector<Range_P> phi_p_n(lfsv_p_n.size());
669 FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
670 FESwitch_P::basis(lfsv_p_n.finiteElement()).evaluateFunction(local_n,phi_p_n);
673 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
674 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_n(lfsu_v_n.size());
675 VectorBasisSwitch_V::jacobian
676 (FESwitch_V::basis(lfsv_v_s.finiteElement()), geo_inside, local_s, jac_phi_v_s);
677 VectorBasisSwitch_V::jacobian
678 (FESwitch_V::basis(lfsv_v_n.finiteElement()), geo_outside, local_n, jac_phi_v_n);
680 auto normal = ig.unitOuterNormal(ip.position());
681 auto weight = ip.weight()*geo.integrationElement(ip.position());
682 auto mu = prm.mu(ig,ip.position());
684 auto factor = mu * weight;
689 for(size_type i=0; i<lfsv_v_s.size(); i++) {
693 add_compute_flux(jac_phi_v_s[i],normal,flux_jac_phi_i);
700 for(size_type j=0; j<lfsu_v_s.size(); j++) {
702 add_compute_flux(jac_phi_v_s[j],normal,flux_jac_phi_j);
704 mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, -0.5 * (flux_jac_phi_j*phi_v_s[i]) * factor);
705 mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, epsilon * 0.5 * (flux_jac_phi_i*phi_v_s[j]) * factor);
706 mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, penalty_factor * (phi_v_s[j]*phi_v_s[i]) * factor);
709 for(size_type j=0; j<lfsu_v_n.size(); j++) {
711 add_compute_flux(jac_phi_v_n[j],normal,flux_jac_phi_j);
713 mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -0.5 * (flux_jac_phi_j*phi_v_s[i]) * factor);
714 mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -epsilon * 0.5 * (flux_jac_phi_i*phi_v_n[j]) * factor);
715 mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -penalty_factor * (phi_v_n[j]*phi_v_s[i]) * factor);
721 for(size_type j=0; j<lfsu_p_s.size(); j++) {
722 mat_ss.accumulate(lfsv_v_s,i,lfsu_p_s,j, 0.5*phi_p_s[j] * (phi_v_s[i]*normal) * weight);
725 for(size_type j=0; j<lfsu_p_n.size(); j++) {
726 mat_sn.accumulate(lfsv_v_s,i,lfsu_p_n,j, 0.5*phi_p_n[j] * (phi_v_s[i]*normal) * weight);
733 for(size_type i=0; i<lfsv_v_n.size(); i++) {
737 add_compute_flux(jac_phi_v_n[i],normal,flux_jac_phi_i);
744 for(size_type j=0; j<lfsu_v_s.size(); j++) {
746 add_compute_flux(jac_phi_v_s[j],normal,flux_jac_phi_j);
748 mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, 0.5 * (flux_jac_phi_j*phi_v_n[i]) * factor);
749 mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, epsilon * 0.5 * (flux_jac_phi_i*phi_v_s[j]) * factor);
750 mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, -penalty_factor * (phi_v_s[j]*phi_v_n[i]) * factor);
753 for(size_type j=0; j<lfsu_v_n.size(); j++) {
755 add_compute_flux(jac_phi_v_n[j],normal,flux_jac_phi_j);
757 mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, 0.5 * (flux_jac_phi_j*phi_v_n[i]) * factor);
758 mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, -epsilon * 0.5 * (flux_jac_phi_i*phi_v_n[j]) * factor);
759 mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, penalty_factor * (phi_v_n[j]*phi_v_n[i]) * factor);
765 for(size_type j=0; j<lfsu_p_s.size(); j++) {
766 mat_ns.accumulate(lfsv_v_n,i,lfsu_p_s,j, -0.5*phi_p_s[j] * (phi_v_n[i]*normal) * weight);
769 for(size_type j=0; j<lfsu_p_n.size(); j++) {
770 mat_nn.accumulate(lfsv_v_n,i,lfsu_p_n,j, -0.5*phi_p_n[j] * (phi_v_n[i]*normal) * weight);
777 for(size_type i=0; i<lfsv_p_s.size(); i++) {
778 for(size_type j=0; j<lfsu_v_s.size(); j++)
779 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);
781 for(size_type j=0; j<lfsu_v_n.size(); j++)
782 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);
785 for(size_type i=0; i<lfsv_p_n.size(); i++) {
786 for(size_type j=0; j<lfsu_v_s.size(); j++)
787 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);
789 for(size_type j=0; j<lfsu_v_n.size(); j++)
790 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);
797 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
798 void alpha_boundary (
const IG& ig,
799 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
803 const unsigned int dim = IG::Entity::dimension;
804 const unsigned int dimw = IG::coorddimension;
808 using LFSV_V = TypeTree::Child<LFSV,_0>;
809 const auto& lfsv_v =
child(lfsv,
_0);
810 const auto& lfsu_v =
child(lfsu,
_0);
812 using LFSV_P = TypeTree::Child<LFSV,_1>;
813 const auto& lfsv_p =
child(lfsv,
_1);
814 const auto& lfsu_p =
child(lfsu,
_1);
817 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
818 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
819 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
820 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
821 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
822 using DF =
typename BasisSwitch_V::DomainField;
823 using RF =
typename BasisSwitch_V::RangeField;
824 using Range_V =
typename BasisSwitch_V::Range;
825 using Range_P =
typename BasisSwitch_P::Range;
826 using size_type =
typename LFSV::Traits::SizeType;
829 const auto& cell_inside = ig.inside();
832 auto geo = ig.geometry();
833 auto geo_inside = cell_inside.geometry();
836 auto geo_in_inside = ig.geometryInInside();
839 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
840 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
841 const int qorder = 2*v_order + det_jac_order + superintegration_order;
843 auto epsilon = prm.epsilonIPSymmetryFactor();
844 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
846 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
849 for (
const auto& ip : quadratureRule(geo,qorder))
852 auto local = geo_in_inside.global(ip.position());
855 std::vector<Range_V> phi_v(lfsv_v.size());
856 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
859 std::vector<Range_P> phi_p(lfsv_p.size());
860 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
863 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
864 VectorBasisSwitch_V::jacobian
865 (FESwitch_V::basis(lfsv_v.finiteElement()), geo_inside, local, jac_phi_v);
869 for (size_type i=0; i<lfsu_p.size(); i++)
870 val_p.axpy(x(lfsu_p,i),phi_p[i]);
875 for (size_type i=0; i<lfsu_v.size(); i++){
876 val_u.axpy(x(lfsu_v,i),phi_v[i]);
877 jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
880 auto normal = ig.unitOuterNormal(ip.position());
881 auto weight = ip.weight()*geo.integrationElement(ip.position());
882 auto mu = prm.mu(ig,ip.position());
885 auto bctype(prm.bctype(ig,ip.position()));
887 if (bctype == BC::VelocityDirichlet) {
889 auto u0(prm.g(cell_inside,local));
890 auto jump = val_u - u0;
894 add_compute_flux(jac_u,normal,flux_jac_u);
896 for (
size_t i=0; i<lfsv_v.size(); i++) {
900 r.accumulate(lfsv_v,i, -mu * (flux_jac_u * phi_v[i]) * weight);
906 add_compute_flux(jac_phi_v[i],normal,flux_jac_phi);
907 r.accumulate(lfsv_v,i, mu * epsilon * (flux_jac_phi*jump) * weight);
912 r.accumulate(lfsv_v,i, mu * (jump*phi_v[i]) * penalty_factor * weight);
917 r.accumulate(lfsv_v,i, val_p * (phi_v[i]*normal) * weight);
923 for(size_type i=0; i<lfsv_p.size(); i++) {
924 r.accumulate(lfsv_p,i, phi_p[i] * (jump*normal) * incomp_scaling * weight);
928 if (bctype == BC::StressNeumann) {
929 auto stress(prm.j(ig,ip.position(),normal));
931 for(size_type i=0; i<lfsv_v.size(); i++) {
932 r.accumulate(lfsv_v,i, (stress*phi_v[i]) * weight);
940 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
941 typename LocalMatrix>
942 void jacobian_boundary (
const IG& ig,
943 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
944 LocalMatrix& mat)
const
947 const unsigned int dim = IG::Entity::dimension;
948 const unsigned int dimw = IG::coorddimension;
952 using LFSV_V = TypeTree::Child<LFSV,_0>;
953 const auto& lfsv_v =
child(lfsv,
_0);
954 const auto& lfsu_v =
child(lfsu,
_0);
956 using LFSV_P = TypeTree::Child<LFSV,_1>;
957 const auto& lfsv_p =
child(lfsv,
_1);
958 const auto& lfsu_p =
child(lfsu,
_1);
961 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
962 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
963 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
964 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
965 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
966 using DF =
typename BasisSwitch_V::DomainField;
967 using RF =
typename BasisSwitch_V::RangeField;
968 using Range_V =
typename BasisSwitch_V::Range;
969 using Range_P =
typename BasisSwitch_P::Range;
970 using size_type =
typename LFSV::Traits::SizeType;
973 const auto& cell_inside = ig.inside();
976 auto geo = ig.geometry();
977 auto geo_inside = cell_inside.geometry();
980 auto geo_in_inside = ig.geometryInInside();
983 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
984 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
985 const int qorder = 2*v_order + det_jac_order + superintegration_order;
987 auto epsilon = prm.epsilonIPSymmetryFactor();
988 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
990 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
993 for (
const auto& ip : quadratureRule(geo,qorder))
996 auto local = geo_in_inside.global(ip.position());
999 std::vector<Range_V> phi_v(lfsv_v.size());
1000 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1003 std::vector<Range_P> phi_p(lfsv_p.size());
1004 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1007 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
1008 VectorBasisSwitch_V::jacobian
1009 (FESwitch_V::basis(lfsv_v.finiteElement()), geo_inside, local, jac_phi_v);
1011 auto normal = ig.unitOuterNormal(ip.position());
1012 auto weight = ip.weight()*geo.integrationElement(ip.position());
1013 auto mu = prm.mu(ig,ip.position());
1016 auto bctype(prm.bctype(ig,ip.position()));
1018 if (bctype == BC::VelocityDirichlet) {
1020 for(size_type i=0; i<lfsv_v.size(); i++) {
1023 add_compute_flux(jac_phi_v[i],normal,flux_jac_phi_i);
1025 for(size_type j=0; j<lfsu_v.size(); j++) {
1031 add_compute_flux(jac_phi_v[j],normal,flux_jac_phi_j);
1033 mat.accumulate(lfsv_v,i,lfsu_v,j, -mu * (flux_jac_phi_j*phi_v[i]) * weight);
1034 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * epsilon * (flux_jac_phi_i*phi_v[j]) *weight);
1039 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (phi_v[j]*phi_v[i]) * penalty_factor * weight);
1045 for(size_type j=0; j<lfsu_p.size(); j++) {
1046 mat.accumulate(lfsv_v,i,lfsu_p,j, phi_p[j] * (phi_v[i]*normal) * weight);
1053 for(size_type i=0; i<lfsv_p.size(); i++) {
1054 for(size_type j=0; j<lfsu_v.size(); j++) {
1055 mat.accumulate(lfsv_p,i,lfsu_v,j, phi_p[i] * (phi_v[j]*normal) * incomp_scaling * weight);
1066 template<
typename EG,
typename LFSV,
typename R>
1067 void lambda_volume (
const EG& eg,
const LFSV& lfsv, R& r)
const
1069 const unsigned int dim = EG::Geometry::mydimension;
1073 using LFSV_V = TypeTree::Child<LFSV,_0>;
1074 const auto& lfsv_v =
child(lfsv,
_0);
1077 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
1078 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
1079 using Range_V =
typename BasisSwitch_V::Range;
1080 using size_type =
typename LFSV::Traits::SizeType;
1083 const auto& cell = eg.entity();
1086 auto geo = eg.geometry();
1089 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1090 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
1091 const int qorder = 2*v_order + det_jac_order + superintegration_order;
1094 for (
const auto& ip : quadratureRule(geo,qorder))
1096 auto local = ip.position();
1100 std::vector<Range_V> phi_v(lfsv_v.size());
1101 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1103 auto weight = ip.weight() * geo.integrationElement(ip.position());
1106 auto fval(prm.f(cell,local));
1111 for(size_type i=0; i<lfsv_v.size(); i++)
1112 r.accumulate(lfsv_v,i, -(fval*phi_v[i]) * weight);
1119 template<
class M,
class RF>
1120 static void contraction(
const M & a,
const M & b, RF & v)
1123 const int an = a.N();
1124 const int am = a.M();
1125 for(
int r=0; r<an; ++r)
1126 for(
int c=0; c<am; ++c)
1127 v += a[r][c] * b[r][c];
1130 template<
class DU,
class R>
1131 static void add_compute_flux(
const DU & du,
const R & n, R & result)
1133 const int N = du.N();
1134 const int M = du.M();
1135 for(
int r=0; r<N; ++r)
1136 for(
int c=0; c<M; ++c)
1137 result[r] += du[r][c] * n[c];
1141 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
Wrapper class for geometries.
Definition: geometry.hh:71
Implementation::JacobianInverseTransposed JacobianInverseTransposed
type of jacobian inverse transposed
Definition: geometry.hh:120
static constexpr int coorddimension
dimension of embedding coordinate system
Definition: geometry.hh:97
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
Return inverse of transposed of Jacobian.
Definition: geometry.hh:328
A local operator for solving the Navier-Stokes equations using a DG discretization and a vector-value...
Definition: dgnavierstokesvelvecfem.hh:105
DGNavierStokesVelVecFEM(PRM &_prm, int _superintegration_order=0)
Constructor.
Definition: dgnavierstokesvelvecfem.hh:142
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.
Traits for type conversions and type information.
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
A hierarchical structure of string parameters.
Definition: stokesparameter.hh:32
typename Basis::Traits::RangeFieldType RangeField
export field type of the values
Definition: dgnavierstokesvelvecfem.hh:67
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:73
typename Basis::Traits::DomainType DomainLocal
export vector type of the local coordinates
Definition: dgnavierstokesvelvecfem.hh:65