4 #ifndef DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESVELVECFEM_HH
5 #define DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESVELVECFEM_HH
7 #include <dune/common/exceptions.hh>
8 #include <dune/common/fvector.hh>
9 #include <dune/common/parametertree.hh>
11 #include <dune/geometry/quadraturerules.hh>
13 #include <dune/localfunctions/common/interfaceswitch.hh>
25 #define PBLOCK (- VBLOCK + 1)
30 template<
class Basis,
class Dummy =
void>
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,
44 Geometry::coorddimension> >& jac)
46 jac.resize(basis.size());
47 basis.evaluateJacobian(xl, jac);
54 Basis, typename enable_if<
58 Basis::Traits::dimDomain
67 typedef typename Basis::Traits::RangeFieldType
RangeField;
69 static const std::size_t
dimRange = Basis::Traits::dimRange;
72 template<
typename Geometry>
73 static void jacobian(
const Basis& basis,
const Geometry& geometry,
74 const DomainLocal& xl,
75 std::vector<FieldMatrix<RangeField, dimRange,
76 Geometry::coorddimension> >& jac)
78 jac.resize(basis.size());
80 std::vector<FieldMatrix<
82 basis.evaluateJacobian(xl, ljac);
84 const typename Geometry::JacobianInverseTransposed& geojac =
85 geometry.jacobianInverseTransposed(xl);
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 typedef typename PRM::Traits::RangeField RF;
112 static const bool navier = PRM::assemble_navier;
113 static const bool full_tensor = PRM::assemble_full_tensor;
143 prm(_prm), superintegration_order(_superintegration_order),
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::dimension;
168 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokesVelVecFEM");
170 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
171 const LFSV_V& lfsv_v = lfsv.template child<VBLOCK>();
172 const LFSV_V& lfsu_v = lfsu.template child<VBLOCK>();
174 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
175 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
176 const LFSV_P& lfsu_p = lfsu.template child<PBLOCK>();
179 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
180 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
182 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
183 typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
184 typedef typename BasisSwitch_V::DomainField DF;
185 typedef typename BasisSwitch_V::RangeField RF;
186 typedef typename BasisSwitch_V::Range Range_V;
187 typedef typename BasisSwitch_P::Range Range_P;
188 typedef typename LFSV::Traits::SizeType size_type;
191 Dune::GeometryType gt = eg.geometry().type();
192 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
193 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
194 const int jac_order = gt.isSimplex() ? 0 : 1;
195 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
196 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
198 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
201 for (
const auto& ip : rule)
203 const Dune::FieldVector<DF,dim> local = ip.position();
204 const RF mu = prm.mu(eg,local);
205 const RF rho = prm.rho(eg,local);
208 std::vector<Range_V> phi_v(lfsv_v.size());
211 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
212 for (size_type i=0; i<lfsu_v.size(); i++)
213 val_u.axpy(x(lfsu_v,i),phi_v[i]);
217 std::vector<Range_P> phi_p(lfsv_p.size());
218 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
222 for (size_type i=0; i<lfsu_p.size(); i++)
223 val_p.axpy(x(lfsu_p,i),phi_p[i]);
226 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
227 VectorBasisSwitch_V::jacobian
228 (FESwitch_V::basis(lfsv_v.finiteElement()), eg.geometry(), local, jac_phi_v);
231 std::vector<RF> div_phi_v(lfsv_v.size(),0.0);
232 for (size_type i=0; i<lfsv_v.size(); i++)
233 for (size_type d=0; d<
dim; d++)
234 div_phi_v[i] += jac_phi_v[i][d][d];
237 Dune::FieldMatrix<RF,dim,dim> jac_u(0.0);
239 for (size_type i=0; i<lfsu_v.size(); i++){
240 jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
241 div_u += x(lfsu_v,i) * div_phi_v[i];
244 const RF detj = eg.geometry().integrationElement(ip.position());
245 const RF weight = ip.weight() * detj;
247 for (size_type i=0; i<lfsv_v.size(); i++) {
251 RF dvdu(0); contraction(jac_u,jac_phi_v[i],dvdu);
252 r.accumulate(lfsv_v, i, dvdu * mu * weight);
257 r.accumulate(lfsv_v, i, - div_phi_v[i] * val_p * weight);
264 Range_V nabla_u_u(0.0);
265 jac_u.mv(val_u,nabla_u_u);
266 r.accumulate(lfsv_v, i, rho * (nabla_u_u*phi_v[i]) * weight);
271 for (size_type i=0; i<lfsv_p.size(); i++) {
275 r.accumulate(lfsv_p, i, - div_u * phi_p[i] * incomp_scaling * weight);
282 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
287 const unsigned int dim = EG::Geometry::dimension;
291 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokesVelVecFEM");
293 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
294 const LFSV_V& lfsv_v = lfsv.template child<VBLOCK>();
295 const LFSV_V& lfsu_v = lfsu.template child<VBLOCK>();
297 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
298 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
299 const LFSV_P& lfsu_p = lfsu.template child<PBLOCK>();
302 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
303 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
305 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
306 typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
307 typedef typename BasisSwitch_V::DomainField DF;
308 typedef typename BasisSwitch_V::RangeField RF;
309 typedef typename BasisSwitch_V::Range Range_V;
310 typedef typename BasisSwitch_P::Range Range_P;
311 typedef typename LFSV::Traits::SizeType size_type;
314 Dune::GeometryType gt = eg.geometry().type();
315 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
316 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
317 const int jac_order = gt.isSimplex() ? 0 : 1;
318 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
319 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
321 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
324 for (
const auto& ip : rule)
326 const Dune::FieldVector<DF,dim> local = ip.position();
327 const RF mu = prm.mu(eg,local);
328 const RF rho = prm.rho(eg,local);
331 std::vector<Range_V> phi_v(lfsv_v.size());
334 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
335 for (size_type i=0; i<lfsu_v.size(); i++)
336 val_u.axpy(x(lfsu_v,i),phi_v[i]);
340 std::vector<Range_P> phi_p(lfsv_p.size());
341 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
344 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
345 VectorBasisSwitch_V::jacobian
346 (FESwitch_V::basis(lfsv_v.finiteElement()), eg.geometry(), local, jac_phi_v);
348 assert(lfsu_v.size() == lfsv_v.size());
350 std::vector<RF> div_phi_v(lfsv_v.size(),0.0);
351 for (size_type i=0; i<lfsv_v.size(); i++)
352 for (size_type d=0; d<
dim; d++)
353 div_phi_v[i] += jac_phi_v[i][d][d];
356 Dune::FieldMatrix<RF,dim,dim> jac_u(0.0);
358 for (size_type i=0; i<lfsu_v.size(); i++){
359 jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
363 const RF detj = eg.geometry().integrationElement(ip.position());
364 const RF weight = ip.weight() * detj;
366 for(size_type i=0; i<lfsv_v.size(); i++) {
368 for(size_type j=0; j<lfsu_v.size(); j++) {
372 RF dvdu(0.0); contraction(jac_phi_v[j],jac_phi_v[i],dvdu);
373 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * dvdu * weight);
380 Range_V nabla_u_phi_v_j(0.0);
381 jac_u.mv(phi_v[j],nabla_u_phi_v_j);
383 Range_V nabla_phi_v_j_u(0.0);
384 jac_phi_v[j].mv(val_u,nabla_phi_v_j_u);
385 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);
390 for(size_type j=0; j<lfsv_p.size(); j++) {
395 mat.accumulate(lfsv_v,i,lfsu_p,j, -phi_p[j] * div_phi_v[i] * weight);
396 mat.accumulate(lfsv_p,j,lfsu_v,i, -phi_p[j] * div_phi_v[i] * incomp_scaling * weight);
404 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
406 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
407 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
408 R& r_s, R& r_n)
const
411 const unsigned int dim = IG::Geometry::dimension;
412 const unsigned int dimw = IG::Geometry::dimensionworld;
416 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokesVelVecFEM");
418 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
419 const LFSV_V& lfsv_v_s = lfsv_s.template child<VBLOCK>();
420 const LFSV_V& lfsu_v_s = lfsu_s.template child<VBLOCK>();
421 const LFSV_V& lfsv_v_n = lfsv_n.template child<VBLOCK>();
422 const LFSV_V& lfsu_v_n = lfsu_n.template child<VBLOCK>();
424 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
425 const LFSV_P& lfsv_p_s = lfsv_s.template child<PBLOCK>();
426 const LFSV_P& lfsu_p_s = lfsu_s.template child<PBLOCK>();
427 const LFSV_P& lfsv_p_n = lfsv_n.template child<PBLOCK>();
428 const LFSV_P& lfsu_p_n = lfsu_n.template child<PBLOCK>();
431 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
432 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
434 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
435 typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
436 typedef typename BasisSwitch_V::DomainField DF;
437 typedef typename BasisSwitch_V::RangeField RF;
438 typedef typename BasisSwitch_V::Range Range_V;
439 typedef typename BasisSwitch_P::Range Range_P;
440 typedef typename LFSV::Traits::SizeType size_type;
443 auto inside_cell = ig.inside();
444 auto outside_cell = ig.outside();
447 Dune::GeometryType gtface = ig.geometry().type();
448 const int v_order = FESwitch_V::basis(lfsv_v_s.finiteElement()).order();
449 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
450 const int qorder = 2*v_order + det_jac_order + superintegration_order;
451 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
453 const int epsilon = prm.epsilonIPSymmetryFactor();
454 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
457 for (
const auto& ip : rule)
461 Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(ip.position());
462 Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(ip.position());
464 const RF penalty_factor = prm.getFaceIP(ig,ip.position());
467 std::vector<Range_V> phi_v_s(lfsv_v_s.size());
468 std::vector<Range_V> phi_v_n(lfsv_v_n.size());
469 FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
470 FESwitch_V::basis(lfsv_v_n.finiteElement()).evaluateFunction(local_n,phi_v_n);
473 std::vector<Range_P> phi_p_s(lfsv_p_s.size());
474 std::vector<Range_P> phi_p_n(lfsv_p_n.size());
475 FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
476 FESwitch_P::basis(lfsv_p_n.finiteElement()).evaluateFunction(local_n,phi_p_n);
479 Range_P val_p_s(0.0);
480 Range_P val_p_n(0.0);
481 for (size_type i=0; i<lfsu_p_s.size(); i++)
482 val_p_s.axpy(x_s(lfsu_p_s,i),phi_p_s[i]);
483 for (size_type i=0; i<lfsu_p_n.size(); i++)
484 val_p_n.axpy(x_n(lfsu_p_n,i),phi_p_n[i]);
487 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
488 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_n(lfsu_v_n.size());
489 VectorBasisSwitch_V::jacobian
490 (FESwitch_V::basis(lfsv_v_s.finiteElement()), inside_cell.geometry(), local_s, jac_phi_v_s);
491 VectorBasisSwitch_V::jacobian
492 (FESwitch_V::basis(lfsv_v_n.finiteElement()), outside_cell.geometry(), local_n, jac_phi_v_n);
495 Range_V val_u_s(0.0);
496 Range_V val_u_n(0.0);
497 Dune::FieldMatrix<RF,dim,dim> jac_u_s(0.0);
498 Dune::FieldMatrix<RF,dim,dim> jac_u_n(0.0);
499 for (size_type i=0; i<lfsu_v_s.size(); i++){
500 val_u_s.axpy(x_s(lfsu_v_s,i),phi_v_s[i]);
501 jac_u_s.axpy(x_s(lfsu_v_s,i),jac_phi_v_s[i]);
503 for (size_type i=0; i<lfsu_v_n.size(); i++){
504 val_u_n.axpy(x_n(lfsu_v_n,i),phi_v_n[i]);
505 jac_u_n.axpy(x_n(lfsu_v_n,i),jac_phi_v_n[i]);
508 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
509 const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
510 const RF mu = prm.mu(ig,ip.position());
512 const RF factor = mu * weight;
515 const Dune::FieldVector<DF,dimw> jump = val_u_s - val_u_n;
518 const RF mean_p = 0.5*(val_p_s + val_p_n);
521 Dune::FieldVector<DF,dimw> flux_jac_u(0.0);
522 add_compute_flux(jac_u_s,normal,flux_jac_u);
523 add_compute_flux(jac_u_n,normal,flux_jac_u);
527 for (
size_t i=0; i<lfsv_v_s.size(); i++) {
531 r_s.accumulate(lfsv_v_s, i, -(flux_jac_u * phi_v_s[i]) * factor);
536 Dune::FieldVector<DF,dimw> flux_jac_phi(0.0);
537 add_compute_flux(jac_phi_v_s[i],normal,flux_jac_phi);
538 r_s.accumulate(lfsv_v_s, i, epsilon * 0.5 * (flux_jac_phi * jump) * factor);
543 r_s.accumulate(lfsv_v_s,i, penalty_factor * (jump*phi_v_s[i]) * weight);
548 r_s.accumulate(lfsv_v_s,i, mean_p * (phi_v_s[i]*normal) * weight);
552 for (
size_t i=0; i<lfsv_v_n.size(); i++) {
556 r_n.accumulate(lfsv_v_n, i, (flux_jac_u * phi_v_n[i]) * factor);
561 Dune::FieldVector<DF,dimw> flux_jac_phi(0.0);
562 add_compute_flux(jac_phi_v_n[i],normal,flux_jac_phi);
563 r_n.accumulate(lfsv_v_n, i, epsilon * 0.5 * (flux_jac_phi * jump) * factor);
568 r_n.accumulate(lfsv_v_n,i, -penalty_factor * (jump*phi_v_n[i]) * weight);
573 r_n.accumulate(lfsv_v_n,i, -mean_p * (phi_v_n[i]*normal) * weight);
579 for (
size_t i=0; i<lfsv_p_s.size(); i++)
580 r_s.accumulate(lfsv_p_s,i, 0.5*phi_p_s[i] * (jump*normal) * incomp_scaling * weight);
581 for (
size_t i=0; i<lfsv_p_n.size(); i++)
582 r_n.accumulate(lfsv_p_n,i, 0.5*phi_p_n[i] * (jump*normal) * incomp_scaling * weight);
588 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
591 const LFSU& lfsu_s,
const X&,
const LFSV& lfsv_s,
592 const LFSU& lfsu_n,
const X&,
const LFSV& lfsv_n,
597 const unsigned int dim = IG::Geometry::dimension;
598 const unsigned int dimw = IG::Geometry::dimensionworld;
602 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokesVelVecFEM");
604 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
605 const LFSV_V& lfsv_v_s = lfsv_s.template child<VBLOCK>();
606 const LFSV_V& lfsu_v_s = lfsu_s.template child<VBLOCK>();
607 const LFSV_V& lfsv_v_n = lfsv_n.template child<VBLOCK>();
608 const LFSV_V& lfsu_v_n = lfsu_n.template child<VBLOCK>();
610 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
611 const LFSV_P& lfsv_p_s = lfsv_s.template child<PBLOCK>();
612 const LFSV_P& lfsu_p_s = lfsu_s.template child<PBLOCK>();
613 const LFSV_P& lfsv_p_n = lfsv_n.template child<PBLOCK>();
614 const LFSV_P& lfsu_p_n = lfsu_n.template child<PBLOCK>();
617 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
618 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
620 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
621 typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
622 typedef typename BasisSwitch_V::DomainField DF;
623 typedef typename BasisSwitch_V::RangeField RF;
624 typedef typename BasisSwitch_V::Range Range_V;
625 typedef typename BasisSwitch_P::Range Range_P;
626 typedef typename LFSV::Traits::SizeType size_type;
629 auto inside_cell = ig.inside();
630 auto outside_cell = ig.outside();
633 Dune::GeometryType gtface = ig.geometry().type();
634 const int v_order = FESwitch_V::basis(lfsv_v_s.finiteElement()).order();
635 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
636 const int qorder = 2*v_order + det_jac_order + superintegration_order;
637 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
639 const int epsilon = prm.epsilonIPSymmetryFactor();
640 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
643 for (
const auto& ip : rule)
647 Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(ip.position());
648 Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(ip.position());
650 const RF penalty_factor = prm.getFaceIP(ig,ip.position());
653 std::vector<Range_V> phi_v_s(lfsv_v_s.size());
654 std::vector<Range_V> phi_v_n(lfsv_v_n.size());
655 FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
656 FESwitch_V::basis(lfsv_v_n.finiteElement()).evaluateFunction(local_n,phi_v_n);
659 std::vector<Range_P> phi_p_s(lfsv_p_s.size());
660 std::vector<Range_P> phi_p_n(lfsv_p_n.size());
661 FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
662 FESwitch_P::basis(lfsv_p_n.finiteElement()).evaluateFunction(local_n,phi_p_n);
665 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
666 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_n(lfsu_v_n.size());
667 VectorBasisSwitch_V::jacobian
668 (FESwitch_V::basis(lfsv_v_s.finiteElement()), inside_cell.geometry(), local_s, jac_phi_v_s);
669 VectorBasisSwitch_V::jacobian
670 (FESwitch_V::basis(lfsv_v_n.finiteElement()), outside_cell.geometry(), local_n, jac_phi_v_n);
672 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
673 const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
674 const RF mu = prm.mu(ig,ip.position());
676 const RF factor = mu * weight;
681 for(size_type i=0; i<lfsv_v_s.size(); i++) {
684 Dune::FieldVector<DF,dimw> flux_jac_phi_i(0.0);
685 add_compute_flux(jac_phi_v_s[i],normal,flux_jac_phi_i);
692 for(size_type j=0; j<lfsu_v_s.size(); j++) {
693 Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
694 add_compute_flux(jac_phi_v_s[j],normal,flux_jac_phi_j);
696 mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, -0.5 * (flux_jac_phi_j*phi_v_s[i]) * factor);
697 mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, epsilon * 0.5 * (flux_jac_phi_i*phi_v_s[j]) * factor);
698 mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, penalty_factor * (phi_v_s[j]*phi_v_s[i]) * weight);
701 for(size_type j=0; j<lfsu_v_n.size(); j++) {
702 Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
703 add_compute_flux(jac_phi_v_n[j],normal,flux_jac_phi_j);
705 mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -0.5 * (flux_jac_phi_j*phi_v_s[i]) * factor);
706 mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -epsilon * 0.5 * (flux_jac_phi_i*phi_v_n[j]) * factor);
707 mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -penalty_factor * (phi_v_n[j]*phi_v_s[i]) * weight);
713 for(size_type j=0; j<lfsu_p_s.size(); j++) {
714 mat_ss.accumulate(lfsv_v_s,i,lfsu_p_s,j, 0.5*phi_p_s[j] * (phi_v_s[i]*normal) * weight);
717 for(size_type j=0; j<lfsu_p_n.size(); j++) {
718 mat_sn.accumulate(lfsv_v_s,i,lfsu_p_n,j, 0.5*phi_p_n[j] * (phi_v_s[i]*normal) * weight);
725 for(size_type i=0; i<lfsv_v_n.size(); i++) {
728 Dune::FieldVector<DF,dimw> flux_jac_phi_i(0.0);
729 add_compute_flux(jac_phi_v_n[i],normal,flux_jac_phi_i);
736 for(size_type j=0; j<lfsu_v_s.size(); j++) {
737 Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
738 add_compute_flux(jac_phi_v_s[j],normal,flux_jac_phi_j);
740 mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, 0.5 * (flux_jac_phi_j*phi_v_n[i]) * factor);
741 mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, epsilon * 0.5 * (flux_jac_phi_i*phi_v_s[j]) * factor);
742 mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, -penalty_factor * (phi_v_s[j]*phi_v_n[i]) * weight);
745 for(size_type j=0; j<lfsu_v_n.size(); j++) {
746 Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
747 add_compute_flux(jac_phi_v_n[j],normal,flux_jac_phi_j);
749 mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, 0.5 * (flux_jac_phi_j*phi_v_n[i]) * factor);
750 mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, -epsilon * 0.5 * (flux_jac_phi_i*phi_v_n[j]) * factor);
751 mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, penalty_factor * (phi_v_n[j]*phi_v_n[i]) * weight);
757 for(size_type j=0; j<lfsu_p_s.size(); j++) {
758 mat_ns.accumulate(lfsv_v_n,i,lfsu_p_s,j, -0.5*phi_p_s[j] * (phi_v_n[i]*normal) * weight);
761 for(size_type j=0; j<lfsu_p_n.size(); j++) {
762 mat_nn.accumulate(lfsv_v_n,i,lfsu_p_n,j, -0.5*phi_p_n[j] * (phi_v_n[i]*normal) * weight);
769 for(size_type i=0; i<lfsv_p_s.size(); i++) {
770 for(size_type j=0; j<lfsu_v_s.size(); j++)
771 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);
773 for(size_type j=0; j<lfsu_v_n.size(); j++)
774 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);
777 for(size_type i=0; i<lfsv_p_n.size(); i++) {
778 for(size_type j=0; j<lfsu_v_s.size(); j++)
779 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);
781 for(size_type j=0; j<lfsu_v_n.size(); j++)
782 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);
789 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
791 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
795 const unsigned int dim = IG::Geometry::dimension;
796 const unsigned int dimw = IG::Geometry::dimensionworld;
800 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokesVelVecFEM");
802 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
803 const LFSV_V& lfsv_v = lfsv.template child<VBLOCK>();
804 const LFSV_V& lfsu_v = lfsu.template child<VBLOCK>();
806 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
807 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
808 const LFSV_P& lfsu_p = lfsu.template child<PBLOCK>();
811 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
812 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
814 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
815 typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
816 typedef typename BasisSwitch_V::DomainField DF;
817 typedef typename BasisSwitch_V::RangeField RF;
818 typedef typename BasisSwitch_V::Range Range_V;
819 typedef typename BasisSwitch_P::Range Range_P;
820 typedef typename LFSV::Traits::SizeType size_type;
823 auto inside_cell = ig.inside();
826 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
827 Dune::GeometryType gtface = ig.geometry().type();
828 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-1);
829 const int qorder = 2*v_order + det_jac_order + superintegration_order;
830 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
832 const int epsilon = prm.epsilonIPSymmetryFactor();
833 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
836 for (
const auto& ip : rule)
839 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
841 const RF penalty_factor = prm.getFaceIP(ig,ip.position() );
844 std::vector<Range_V> phi_v(lfsv_v.size());
845 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
848 std::vector<Range_P> phi_p(lfsv_p.size());
849 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
852 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
853 VectorBasisSwitch_V::jacobian
854 (FESwitch_V::basis(lfsv_v.finiteElement()), inside_cell.geometry(), local, jac_phi_v);
858 for (size_type i=0; i<lfsu_p.size(); i++)
859 val_p.axpy(x(lfsu_p,i),phi_p[i]);
863 Dune::FieldMatrix<RF,dim,dim> jac_u(0.0);
864 for (size_type i=0; i<lfsu_v.size(); i++){
865 val_u.axpy(x(lfsu_v,i),phi_v[i]);
866 jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
869 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
870 const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
871 const RF mu = prm.mu(ig,ip.position());
874 typename PRM::Traits::BoundaryCondition::Type bctype(prm.bctype(ig,ip.position()));
878 typename PRM::Traits::VelocityRange u0(prm.g(inside_cell,local));
879 const Dune::FieldVector<DF,dimw> jump = val_u - u0;
882 Dune::FieldVector<DF,dimw> flux_jac_u(0.0);
883 add_compute_flux(jac_u,normal,flux_jac_u);
885 for (
size_t i=0; i<lfsv_v.size(); i++) {
889 r.accumulate(lfsv_v,i, -mu * (flux_jac_u * phi_v[i]) * weight);
894 Dune::FieldVector<DF,dimw> flux_jac_phi(0.0);
895 add_compute_flux(jac_phi_v[i],normal,flux_jac_phi);
896 r.accumulate(lfsv_v,i, mu * epsilon * (flux_jac_phi*jump) * weight);
901 r.accumulate(lfsv_v,i, (jump*phi_v[i]) * penalty_factor * weight);
906 r.accumulate(lfsv_v,i, val_p * (phi_v[i]*normal) * weight);
912 for(size_type i=0; i<lfsv_p.size(); i++) {
913 r.accumulate(lfsv_p,i, phi_p[i] * (jump*normal) * incomp_scaling * weight);
918 typename PRM::Traits::VelocityRange stress(prm.j(ig,ip.position(),normal));
920 for(size_type i=0; i<lfsv_v.size(); i++) {
921 r.accumulate(lfsv_v,i, (stress*phi_v[i]) * weight);
929 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
932 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
936 const unsigned int dim = IG::Geometry::dimension;
937 const unsigned int dimw = IG::Geometry::dimensionworld;
941 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokesVelVecFEM");
943 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
944 const LFSV_V& lfsv_v = lfsv.template child<VBLOCK>();
945 const LFSV_V& lfsu_v = lfsu.template child<VBLOCK>();
947 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
948 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
949 const LFSV_P& lfsu_p = lfsu.template child<PBLOCK>();
952 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
953 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
955 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
956 typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
957 typedef typename BasisSwitch_V::DomainField DF;
958 typedef typename BasisSwitch_V::RangeField RF;
959 typedef typename BasisSwitch_V::Range Range_V;
960 typedef typename BasisSwitch_P::Range Range_P;
961 typedef typename LFSV::Traits::SizeType size_type;
964 auto inside_cell = ig.inside();
967 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
968 Dune::GeometryType gtface = ig.geometry().type();
969 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-1);
970 const int qorder = 2*v_order + det_jac_order + superintegration_order;
971 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
973 const int epsilon = prm.epsilonIPSymmetryFactor();
974 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
977 for (
const auto& ip : rule)
980 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
982 const RF penalty_factor = prm.getFaceIP(ig,ip.position() );
985 std::vector<Range_V> phi_v(lfsv_v.size());
986 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
989 std::vector<Range_P> phi_p(lfsv_p.size());
990 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
993 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
994 VectorBasisSwitch_V::jacobian
995 (FESwitch_V::basis(lfsv_v.finiteElement()), inside_cell.geometry(), local, jac_phi_v);
997 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
998 const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
999 const RF mu = prm.mu(ig,ip.position());
1002 typename PRM::Traits::BoundaryCondition::Type bctype(prm.bctype(ig,ip.position()));
1006 for(size_type i=0; i<lfsv_v.size(); i++) {
1008 Dune::FieldVector<DF,dimw> flux_jac_phi_i(0.0);
1009 add_compute_flux(jac_phi_v[i],normal,flux_jac_phi_i);
1011 for(size_type j=0; j<lfsu_v.size(); j++) {
1016 Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
1017 add_compute_flux(jac_phi_v[j],normal,flux_jac_phi_j);
1019 mat.accumulate(lfsv_v,i,lfsu_v,j, -mu * (flux_jac_phi_j*phi_v[i]) * weight);
1020 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * epsilon * (flux_jac_phi_i*phi_v[j]) *weight);
1025 mat.accumulate(lfsv_v,i,lfsu_v,j, (phi_v[j]*phi_v[i]) * penalty_factor * weight);
1031 for(size_type j=0; j<lfsu_p.size(); j++) {
1032 mat.accumulate(lfsv_v,i,lfsu_p,j, phi_p[j] * (phi_v[i]*normal) * weight);
1039 for(size_type i=0; i<lfsv_p.size(); i++) {
1040 for(size_type j=0; j<lfsu_v.size(); j++) {
1041 mat.accumulate(lfsv_p,i,lfsu_v,j, phi_p[i] * (phi_v[j]*normal) * incomp_scaling * weight);
1052 template<
typename EG,
typename LFSV,
typename R>
1055 const unsigned int dim = EG::Geometry::dimension;
1059 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokesVelVecFEM");
1061 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
1062 const LFSV_V& lfsv_v = lfsv.template child<VBLOCK>();
1065 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
1066 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
1068 typedef typename BasisSwitch_V::DomainField DF;
1069 typedef typename BasisSwitch_V::RangeField RF;
1070 typedef typename BasisSwitch_V::Range Range_V;
1071 typedef typename LFSV::Traits::SizeType size_type;
1074 Dune::GeometryType gt = eg.geometry().type();
1075 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1076 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
1077 const int qorder = 2*v_order + det_jac_order + superintegration_order;
1079 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
1082 for (
const auto& ip : rule)
1084 const Dune::FieldVector<DF,dim> local = ip.position();
1088 std::vector<Range_V> phi_v(lfsv_v.size());
1089 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1091 const RF weight = ip.weight() * eg.geometry().integrationElement(ip.position());
1094 typename PRM::Traits::VelocityRange fval(prm.f(eg,local));
1099 for(size_type i=0; i<lfsv_v.size(); i++)
1100 r.accumulate(lfsv_v,i, -(fval*phi_v[i]) * weight);
1107 template<
class M,
class RF>
1108 static void contraction(
const M & a,
const M & b, RF & v)
1111 const int an = a.N();
1112 const int am = a.M();
1113 for(
int r=0; r<an; ++r)
1114 for(
int c=0; c<am; ++c)
1115 v += a[r][c] * b[r][c];
1118 template<
class DU,
class R>
1119 static void add_compute_flux(
const DU & du,
const R & n, R & result)
1121 const int N = du.N();
1122 const int M = du.M();
1123 for(
int r=0; r<N; ++r)
1124 for(
int c=0; c<M; ++c)
1125 result[r] += du[r][c] * n[c];
1129 const int superintegration_order;
Basis::Traits::RangeFieldType RangeField
export field type of the values
Definition: dgnavierstokesvelvecfem.hh:67
Basis::Traits::DomainType DomainLocal
export vector type of the local coordinates
Definition: dgnavierstokesvelvecfem.hh:65
Definition: dgnavierstokesvelvecfem.hh:31
void setTime(doublet_)
set time for subsequent evaluation
Definition: idefault.hh:104
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: dgnavierstokesvelvecfem.hh:162
Definition: stokesparameter.hh:36
A dense matrix for storing data associated with the degrees of freedom of a pair of LocalFunctionSpac...
Definition: localmatrix.hh:184
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: dgnavierstokesvelvecfem.hh:284
Definition: stokesparameter.hh:35
static const std::size_t dimRange
export dimension of the values
Definition: dgnavierstokesvelvecfem.hh:37
A local operator for solving the stokes equation using a DG discretization and a vector-valued finite...
Definition: dgnavierstokesvelvecfem.hh:101
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: dgnavierstokesvelvecfem.hh:1053
const IG & ig
Definition: constraints.hh:147
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
Definition: dgnavierstokesvelvecfem.hh:122
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:41
Basis::Traits::RangeField RangeField
export field type of the values
Definition: dgnavierstokesvelvecfem.hh:35
Default flags for all local operators.
Definition: flags.hh:18
DGNavierStokesVelVecFEM(PRM &_prm, int _superintegration_order=0)
Constructor.
Definition: dgnavierstokesvelvecfem.hh:142
sparsity pattern generator
Definition: pattern.hh:13
Definition: stokesparameter.hh:32
void preStep(RealType, RealType dt, int)
Definition: dgnavierstokesvelvecfem.hh:148
void setTime(Real t)
Definition: dgnavierstokesvelvecfem.hh:154
Definition: dgnavierstokesvelvecfem.hh:128
Definition: dgnavierstokesvelvecfem.hh:127
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: dgnavierstokesvelvecfem.hh:405
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
Definition: adaptivity.hh:27
Definition: dgnavierstokesvelvecfem.hh:125
static const int dim
Definition: adaptivity.hh:83
void alpha_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: dgnavierstokesvelvecfem.hh:790
sparsity pattern generator
Definition: pattern.hh:29
Definition: dgnavierstokesvelvecfem.hh:119
double RealType
Definition: idefault.hh:92
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: dgnavierstokesvelvecfem.hh:590
Basis::Traits::DomainLocal DomainLocal
export vector type of the local coordinates
Definition: dgnavierstokesvelvecfem.hh:33
Definition: dgnavierstokesvelvecfem.hh:118
const EG & eg
Definition: constraints.hh:280
void jacobian_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: dgnavierstokesvelvecfem.hh:931
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Definition: dgnavierstokesvelvecfem.hh:126