2 #ifndef DUNE_PDELAB_STOKESDG_HH
3 #define DUNE_PDELAB_STOKESDG_HH
5 #warning This file is deprecated, include the header dune/pdelab/localoperator/dgnavierstokes.hh instead!
7 #include <dune/common/exceptions.hh>
8 #include <dune/common/fvector.hh>
9 #include <dune/common/parametertreeparser.hh>
11 #include <dune/geometry/quadraturerules.hh>
13 #include <dune/localfunctions/common/interfaceswitch.hh>
24 #define PBLOCK (- VBLOCK + 1)
35 template<
typename PRM,
bool full_tensor = true>
46 typedef typename PRM::Traits::RangeField RF;
49 typedef typename InstatBase::RealType Real;
54 enum { doPatternVolume =
true };
55 enum { doPatternSkeleton =
true };
58 enum { doSkeletonTwoSided =
false };
61 enum { doAlphaVolume =
true };
62 enum { doAlphaSkeleton =
true };
63 enum { doAlphaBoundary =
true };
64 enum { doLambdaVolume =
true };
65 enum { doLambdaBoundary =
true };
79 DUNE_DEPRECATED_MSG(
"Deprecated in DUNE-PDELab 2.4, use DGNavierStokes instead!")
80 StokesDG (PRM & _prm,
int _superintegration_order=0) :
81 prm(_prm), superintegration_order(_superintegration_order),
94 InstatBase::setTime(t);
100 template<
typename EG,
typename LFSV,
typename R>
104 static const unsigned int dim = EG::Geometry::dimension;
108 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
110 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
111 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
114 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
117 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
118 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
119 const unsigned int vsize = lfsv_v.size();
120 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
121 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
122 const unsigned int psize = lfsv_p.size();
125 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
126 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
127 typedef typename BasisSwitch_V::DomainField DF;
128 typedef typename BasisSwitch_V::Range RT;
129 typedef typename BasisSwitch_V::RangeField RF;
130 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
131 typedef typename LFSV::Traits::SizeType size_type;
134 Dune::GeometryType gt = eg.geometry().type();
135 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
136 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
138 const int qorder = v_order + det_jac_order + superintegration_order;
140 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
143 for (
const auto& ip : rule)
145 const Dune::FieldVector<DF,dim> local = ip.position();
149 std::vector<RT> phi_v(vsize);
150 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
153 std::vector<RT> phi_p(psize);
154 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
156 const RF weight = ip.weight() * eg.geometry().integrationElement(ip.position());
159 typename PRM::Traits::VelocityRange fval(prm.f(eg,local));
164 const RF factor = weight;
165 for (
unsigned int d=0; d<
dim; d++)
167 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
170 for (size_type i=0; i<vsize; i++)
172 RF val = phi_v[i]*factor;
173 r.accumulate(lfsv_v,i, -fval[d] * val);
177 const RF g2 = prm.g2(eg,ip.position());
180 for (
size_t i=0; i<lfsv_p.size(); i++)
182 r.accumulate(lfsv_p,i, g2 * phi_p[i] * factor);
190 template<
typename IG,
typename LFSV,
typename R>
194 static const unsigned int dim = IG::Geometry::dimension;
198 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
200 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
201 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
204 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
207 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
208 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
209 const unsigned int vsize = lfsv_v.size();
210 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
211 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
212 const unsigned int psize = lfsv_p.size();
215 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
216 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
217 typedef typename BasisSwitch_V::DomainField DF;
218 typedef typename BasisSwitch_V::Range RT;
219 typedef typename BasisSwitch_V::RangeField RF;
220 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
223 auto inside_cell = ig.inside();
226 Dune::GeometryType gtface = ig.geometry().type();
227 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
228 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
229 const int jac_order = gtface.isSimplex() ? 0 : 1;
230 const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
231 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
233 const int epsilon = prm.epsilonIPSymmetryFactor();
234 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
237 for (
const auto& ip : rule)
240 Dune::FieldVector<DF,dim-1> flocal = ip.position();
241 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(flocal);
244 const RF penalty_factor = prm.getFaceIP(ig,flocal);
247 std::vector<RT> phi_v(vsize);
248 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
250 std::vector<RT> phi_p(psize);
251 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
253 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
254 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
255 inside_cell.geometry(), local, grad_phi_v);
257 const Dune::FieldVector<DF,dim> normal = ig.unitOuterNormal(ip.position());
258 const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
259 const RF mu = prm.mu(ig,flocal);
262 typename PRM::Traits::BoundaryCondition::Type bctype(prm.bctype(ig,flocal));
264 if (bctype == BC::VelocityDirichlet)
266 typename PRM::Traits::VelocityRange u0(prm.g(inside_cell,local));
271 RF factor = mu * weight;
272 for (
unsigned int i=0;i<vsize;++i)
274 const RF val = (grad_phi_v[i][0]*normal) * factor;
275 for (
unsigned int d=0;d<
dim;++d)
277 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
278 r.accumulate(lfsv_v,i, - epsilon * val * u0[d]);
283 for (
unsigned int dd=0;dd<
dim;++dd)
285 RF Tval = (grad_phi_v[i][0][d]*normal[dd]) * factor;
286 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
287 r.accumulate(lfsv_v_dd,i, - epsilon * Tval * u0[d]);
295 factor = penalty_factor * weight;
296 for (
unsigned int i=0;i<vsize;++i)
298 const RF val = phi_v[i] * factor;
299 for (
unsigned int d=0;d<
dim;++d)
301 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
302 r.accumulate(lfsv_v,i, -val * u0[d] );
308 for (
unsigned int i=0;i<psize;++i)
310 RF val = phi_p[i]*(u0 * normal) * weight;
311 r.accumulate(lfsv_p,i, - val * incomp_scaling);
314 if (bctype == BC::StressNeumann)
316 typename PRM::Traits::VelocityRange stress(prm.j(ig,flocal,normal));
322 for (
unsigned int i=0;i<vsize;++i)
324 for (
unsigned int d=0;d<
dim;++d)
326 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
327 RF val = stress[d]*phi_v[i] * weight;
328 r.accumulate(lfsv_v,i, val);
336 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
342 static const unsigned int dim = EG::Geometry::dimension;
346 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
347 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
348 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
350 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
353 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
354 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
355 const unsigned int vsize = lfsv_v.size();
356 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
357 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
358 const unsigned int psize = lfsv_p.size();
361 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
362 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
363 typedef typename BasisSwitch_V::DomainField DF;
364 typedef typename BasisSwitch_V::Range RT;
365 typedef typename BasisSwitch_V::RangeField RF;
366 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
367 typedef typename LFSV::Traits::SizeType size_type;
370 Dune::GeometryType gt = eg.geometry().type();
371 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
372 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
373 const int jac_order = gt.isSimplex() ? 0 : 1;
374 const int qorder = 2*v_order - 2 + 2*jac_order + det_jac_order + superintegration_order;
375 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
377 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
380 for (
const auto& ip : rule)
382 const Dune::FieldVector<DF,dim> local = ip.position();
383 const RF mu = prm.mu(eg,local);
386 std::vector<RT> phi_p(psize);
387 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
390 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
391 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
392 eg.geometry(), local, grad_phi_v);
394 const RF detj = eg.geometry().integrationElement(ip.position());
395 const RF weight = ip.weight() * detj;
400 const RF factor = mu * weight;
401 for (size_type j=0; j<vsize; j++)
403 for (size_type i=0; i<vsize; i++)
406 RF val = (grad_phi_v[j][0]*grad_phi_v[i][0])*factor;
408 for (
unsigned int d=0; d<
dim; d++)
410 const LFSV_V& lfsv_v_d = lfsv_pfs_v.child(d);
411 mat.accumulate(lfsv_v_d,i,lfsv_v_d,j, val);
415 for (
unsigned int dd=0; dd<
dim; dd++){
416 RF Tval = (grad_phi_v[j][0][d]*grad_phi_v[i][0][dd])*factor;
417 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
418 mat.accumulate(lfsv_v_d,i,lfsv_v_dd,j, Tval);
430 for (size_type j=0; j<psize; j++)
432 RF val = -1.0 * phi_p[j]*weight;
433 for (size_type i=0; i<vsize; i++)
435 for (
unsigned int d=0; d<
dim; d++)
437 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
438 mat.accumulate(lfsv_p,j,lfsv_v,i, val*grad_phi_v[i][0][d] * incomp_scaling);
439 mat.accumulate(lfsv_v,i,lfsv_p,j, val*grad_phi_v[i][0][d]);
447 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
450 const LFSU& lfsu_s,
const X&,
const LFSV& lfsv_s,
451 const LFSU& lfsu_n,
const X&,
const LFSV& lfsv_n,
456 static const unsigned int dim = IG::Geometry::dimension;
457 static const unsigned int dimw = IG::Geometry::dimensionworld;
461 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
463 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
464 const LFSV_PFS_V& lfsv_s_pfs_v = lfsv_s.template child<VBLOCK>();
465 const LFSV_PFS_V& lfsv_n_pfs_v = lfsv_n.template child<VBLOCK>();
467 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
470 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
471 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.template child<0>();
472 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.template child<0>();
473 const unsigned int vsize_s = lfsv_s_v.size();
474 const unsigned int vsize_n = lfsv_n_v.size();
475 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
476 const LFSV_P& lfsv_s_p = lfsv_s.template child<PBLOCK>();
477 const LFSV_P& lfsv_n_p = lfsv_n.template child<PBLOCK>();
478 const unsigned int psize_s = lfsv_s_p.size();
479 const unsigned int psize_n = lfsv_n_p.size();
482 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
483 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
484 typedef typename BasisSwitch_V::DomainField DF;
485 typedef typename BasisSwitch_V::Range RT;
486 typedef typename BasisSwitch_V::RangeField RF;
487 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
490 auto inside_cell = ig.inside();
491 auto outside_cell = ig.outside();
494 Dune::GeometryType gtface = ig.geometry().type();
495 const int v_order = FESwitch_V::basis(lfsv_s_v.finiteElement()).order();
496 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
497 const int qorder = 2*v_order + det_jac_order + superintegration_order;
498 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
500 const int epsilon = prm.epsilonIPSymmetryFactor();
501 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
504 for (
const auto& ip : rule)
508 Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(ip.position());
509 Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(ip.position());
511 const RF penalty_factor = prm.getFaceIP(ig,ip.position());
514 std::vector<RT> phi_v_s(vsize_s);
515 std::vector<RT> phi_v_n(vsize_n);
516 FESwitch_V::basis(lfsv_s_v.finiteElement()).evaluateFunction(local_s,phi_v_s);
517 FESwitch_V::basis(lfsv_n_v.finiteElement()).evaluateFunction(local_n,phi_v_n);
519 std::vector<RT> phi_p_s(psize_s);
520 std::vector<RT> phi_p_n(psize_n);
521 FESwitch_P::basis(lfsv_s_p.finiteElement()).evaluateFunction(local_s,phi_p_s);
522 FESwitch_P::basis(lfsv_n_p.finiteElement()).evaluateFunction(local_n,phi_p_n);
525 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_s(vsize_s);
526 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_s_v.finiteElement()),
527 inside_cell.geometry(), local_s, grad_phi_v_s);
529 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_n(vsize_n);
530 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_n_v.finiteElement()),
531 outside_cell.geometry(), local_n, grad_phi_v_n);
533 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
534 const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
535 const RF mu = prm.mu(ig,ip.position());
540 assert(vsize_s == vsize_n);
541 RF factor = mu * weight;
542 for (
unsigned int i=0;i<vsize_s;++i)
544 for (
unsigned int j=0;j<vsize_s;++j)
546 RF val = (0.5*(grad_phi_v_s[i][0]*normal)*phi_v_s[j]) * factor;
547 for (
unsigned int d=0;d<
dim;++d)
549 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
550 mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v,i, - val);
551 mat_ss.accumulate(lfsv_s_v,i,lfsv_s_v,j, epsilon*val );
556 for (
unsigned int dd=0;dd<
dim;++dd)
558 RF Tval = (0.5*(grad_phi_v_s[i][0][d]*normal[dd])*phi_v_s[j]) * factor;
559 const LFSV_V& lfsv_s_v_dd = lfsv_s_pfs_v.child(dd);
560 mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v_dd,i, - Tval);
561 mat_ss.accumulate(lfsv_s_v_dd,i,lfsv_s_v,j, epsilon*Tval );
566 for (
unsigned int j=0;j<vsize_n;++j)
569 RF val = (-0.5*(grad_phi_v_s[i][0]*normal)*phi_v_n[j]) * factor;
570 for (
unsigned int d=0;d<
dim;++d)
572 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
573 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
574 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v,i,- val);
575 mat_sn.accumulate(lfsv_s_v,i,lfsv_n_v,j, epsilon*val);
580 for (
unsigned int dd=0;dd<
dim;++dd)
582 RF Tval = (-0.5*(grad_phi_v_s[i][0][d]*normal[dd])*phi_v_n[j]) * factor;
583 const LFSV_V& lfsv_s_v_dd = lfsv_s_pfs_v.child(dd);
584 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v_dd,i,- Tval);
585 mat_sn.accumulate(lfsv_s_v_dd,i,lfsv_n_v,j, epsilon*Tval);
591 for (
unsigned int i=0;i<vsize_n;++i)
593 for (
unsigned int j=0;j<vsize_s;++j)
595 RF val = (0.5*(grad_phi_v_n[i][0]*normal)*phi_v_s[j]) * factor;
596 for (
unsigned int d=0;d<
dim;++d)
598 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
599 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
600 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v,i, - val);
601 mat_ns.accumulate(lfsv_n_v,i,lfsv_s_v,j, epsilon*val );
606 for (
unsigned int dd=0;dd<
dim;++dd)
608 RF Tval = (0.5*(grad_phi_v_n[i][0][d]*normal[dd])*phi_v_s[j]) * factor;
609 const LFSV_V& lfsv_n_v_dd = lfsv_n_pfs_v.child(dd);
610 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v_dd,i, - Tval);
611 mat_ns.accumulate(lfsv_n_v_dd,i,lfsv_s_v,j, epsilon*Tval );
616 for (
unsigned int j=0;j<vsize_n;++j)
619 RF val = (-0.5*(grad_phi_v_n[i][0]*normal)*phi_v_n[j]) * factor;
620 for (
unsigned int d=0;d<
dim;++d)
622 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
623 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v,i,- val);
624 mat_nn.accumulate(lfsv_n_v,i,lfsv_n_v,j, epsilon*val);
629 for (
unsigned int dd=0;dd<
dim;++dd)
631 RF Tval = (-0.5*(grad_phi_v_n[i][0][d]*normal[dd])*phi_v_n[j]) * factor;
632 const LFSV_V& lfsv_n_v_dd = lfsv_n_pfs_v.child(dd);
633 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v_dd,i,- Tval);
634 mat_nn.accumulate(lfsv_n_v_dd,i,lfsv_n_v,j, epsilon*Tval);
643 factor = penalty_factor * weight;
644 for (
unsigned int i=0;i<vsize_s;++i)
646 for (
unsigned int j=0;j<vsize_s;++j)
648 RF val = phi_v_s[i]*phi_v_s[j] * factor;
649 for (
unsigned int d=0;d<
dim;++d)
651 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
652 mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v,i, val);
655 for (
unsigned int j=0;j<vsize_n;++j)
657 RF val = phi_v_s[i]*phi_v_n[j] * factor;
658 for (
unsigned int d=0;d<
dim;++d)
660 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
661 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
662 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v,i, - val);
666 for (
unsigned int i=0;i<vsize_n;++i)
668 for (
unsigned int j=0;j<vsize_s;++j)
670 RF val = phi_v_n[i]*phi_v_s[j] * factor;
671 for (
unsigned int d=0;d<
dim;++d)
673 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
674 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
675 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v,i, - val);
678 for (
unsigned int j=0;j<vsize_n;++j)
680 RF val = phi_v_n[i]*phi_v_n[j] * factor;
681 for (
unsigned int d=0;d<
dim;++d)
683 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
684 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v,i, val);
692 for (
unsigned int i=0;i<vsize_s;++i)
694 for (
unsigned int j=0;j<psize_s;++j)
696 for (
unsigned int d=0;d<
dim;++d)
698 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
699 RF val = 0.5*(phi_p_s[j]*normal[d]*phi_v_s[i]) * weight;
700 mat_ss.accumulate(lfsv_s_v,i,lfsv_s_p,j, val);
701 mat_ss.accumulate(lfsv_s_p,j,lfsv_s_v,i, val * incomp_scaling);
704 for (
unsigned int j=0;j<psize_n;++j)
706 for (
unsigned int d=0;d<
dim;++d)
708 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
709 RF val = 0.5*(phi_p_n[j]*normal[d]*phi_v_s[i]) * weight;
710 mat_sn.accumulate(lfsv_s_v,i,lfsv_n_p,j, val);
711 mat_ns.accumulate(lfsv_n_p,j,lfsv_s_v,i, val * incomp_scaling);
715 for (
unsigned int i=0;i<vsize_n;++i)
717 for (
unsigned int j=0;j<psize_s;++j)
719 for (
unsigned int d=0;d<
dim;++d)
721 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
724 RF val = -0.5*(phi_p_s[j]*normal[d]*phi_v_n[i]) * weight;
725 mat_ns.accumulate(lfsv_n_v,i,lfsv_s_p,j, val);
726 mat_sn.accumulate(lfsv_s_p,j,lfsv_n_v,i, val * incomp_scaling);
729 for (
unsigned int j=0;j<psize_n;++j)
731 for (
unsigned int d=0;d<
dim;++d)
733 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
736 RF val = -0.5*(phi_p_n[j]*normal[d]*phi_v_n[i]) * weight;
737 mat_nn.accumulate(lfsv_n_v,i,lfsv_n_p,j, val);
738 mat_nn.accumulate(lfsv_n_p,j,lfsv_n_v,i, val * incomp_scaling);
746 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
749 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
753 static const unsigned int dim = IG::Geometry::dimension;
754 static const unsigned int dimw = IG::Geometry::dimensionworld;
758 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
760 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
761 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
763 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
766 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
767 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
768 const unsigned int vsize = lfsv_v.size();
769 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
770 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
771 const unsigned int psize = lfsv_p.size();
774 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
775 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
776 typedef typename BasisSwitch_V::DomainField DF;
777 typedef typename BasisSwitch_V::Range RT;
778 typedef typename BasisSwitch_V::RangeField RF;
779 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
782 auto inside_cell = ig.inside();
785 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
786 Dune::GeometryType gtface = ig.geometry().type();
787 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-1);
788 const int qorder = 2*v_order + det_jac_order + superintegration_order;
789 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
792 typename PRM::Traits::BoundaryCondition::Type bctype(prm.bctype(ig,rule.begin()->position()));
794 const int epsilon = prm.epsilonIPSymmetryFactor();
795 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
798 for (
const auto& ip : rule)
801 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
803 const RF penalty_factor = prm.getFaceIP(ig,ip.position() );
806 std::vector<RT> phi_v(vsize);
807 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
809 std::vector<RT> phi_p(psize);
810 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
812 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
813 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
814 inside_cell.geometry(), local, grad_phi_v);
816 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
817 const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
818 const RF mu = prm.mu(ig,ip.position());
821 RF slip_factor = 0.0;
823 if (bctype == BC::SlipVelocity)
828 slip_factor = BoundarySlipSwitch::boundarySlip(prm,ig,ip.position());
831 if (bctype == BC::VelocityDirichlet || bctype == BC::SlipVelocity)
833 const RF factor = weight * (1.0 - slip_factor);
838 for (
unsigned int i=0;i<vsize;++i)
840 for (
unsigned int j=0;j<vsize;++j)
842 RF val = ((grad_phi_v[j][0]*normal)*phi_v[i]) * factor * mu;
843 for (
unsigned int d=0;d<
dim;++d)
845 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
846 mat.accumulate(lfsv_v,i,lfsv_v,j, - val);
847 mat.accumulate(lfsv_v,j,lfsv_v,i, epsilon*val);
852 for (
unsigned int dd=0;dd<
dim;++dd)
854 RF Tval = ((grad_phi_v[j][0][d]*normal[dd])*phi_v[i]) * factor * mu;
855 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
856 mat.accumulate(lfsv_v,i,lfsv_v_dd,j, - Tval);
857 mat.accumulate(lfsv_v_dd,j,lfsv_v,i, epsilon*Tval);
867 for (
unsigned int i=0;i<vsize;++i)
869 for (
unsigned int j=0;j<psize;++j)
871 for (
unsigned int d=0;d<
dim;++d)
873 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
874 RF val = (phi_p[j]*normal[d]*phi_v[i]) * weight;
875 mat.accumulate(lfsv_p,j,lfsv_v,i, val * incomp_scaling);
876 mat.accumulate(lfsv_v,i,lfsv_p,j, val);
883 const RF p_factor = penalty_factor * factor;
884 for (
unsigned int i=0;i<vsize;++i)
886 for (
unsigned int j=0;j<vsize;++j)
888 RF val = phi_v[i]*phi_v[j] * p_factor;
889 for (
unsigned int d=0;d<
dim;++d)
891 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
892 mat.accumulate(lfsv_v,j,lfsv_v,i, val);
898 if (bctype == BC::SlipVelocity)
900 const RF factor = weight * (slip_factor);
906 for (
unsigned int i=0;i<vsize;++i)
908 for (
unsigned int j=0;j<vsize;++j)
916 RF val = ten_sum * ((grad_phi_v[j][0]*normal)*phi_v[i]) * factor * mu;
917 for (
unsigned int d=0;d<
dim;++d)
919 const LFSV_V& lfsv_v_d = lfsv_pfs_v.child(d);
921 for (
unsigned int dd=0;dd<
dim;++dd)
923 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
925 mat.accumulate(lfsv_v_dd,i,lfsv_v_d,j, -val*normal[d]*normal[dd]);
926 mat.accumulate(lfsv_v_d,j,lfsv_v_dd,i, epsilon*val*normal[d]*normal[dd]);
935 const RF p_factor = penalty_factor * factor;
936 for (
unsigned int i=0;i<vsize;++i)
938 for (
unsigned int j=0;j<vsize;++j)
940 RF val = phi_v[i]*phi_v[j] * p_factor;
941 for (
unsigned int d=0;d<
dim;++d)
943 const LFSV_V& lfsv_v_d = lfsv_pfs_v.child(d);
944 for (
unsigned int dd=0;dd<
dim;++dd)
946 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
947 mat.accumulate(lfsv_v_d,j,lfsv_v_dd,i, val*normal[d]*normal[dd]);
972 template<
typename PRM,
bool full_tensor = true>
979 typedef typename PRM::Traits::RangeField
RF;
985 using StokesLocalOperator::prm;
986 using StokesLocalOperator::superintegration_order;
988 DUNE_DEPRECATED_MSG(
"Deprecated in DUNE-PDELab 2.4, use DGNavierStokes instead!")
990 : StokesLocalOperator(prm_ ,superintegration_order_)
993 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
999 StokesLocalOperator::jacobian_volume(eg,lfsu,x,lfsv,mat);
1002 static const unsigned int dim = EG::Geometry::dimension;
1006 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
1007 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
1008 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
1010 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
1013 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
1014 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
1015 const unsigned int vsize = lfsv_v.size();
1018 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
1019 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
1020 typedef typename BasisSwitch_V::DomainField DF;
1021 typedef typename BasisSwitch_V::Range RT;
1022 typedef typename BasisSwitch_V::RangeField RF;
1025 Dune::GeometryType gt = eg.geometry().type();
1026 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1027 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
1028 const int jac_order = gt.isSimplex() ? 0 : 1;
1029 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
1030 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
1033 for (
const auto& ip : rule)
1035 const Dune::FieldVector<DF,dim> local = ip.position();
1038 const RF rho = prm.rho(eg,local);
1039 if(rho == 0)
continue;
1042 std::vector<RT> phi_v(vsize);
1043 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1046 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
1047 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1048 eg.geometry(), local, grad_phi_v);
1050 const RF weight = ip.weight() * eg.geometry().integrationElement(ip.position());
1053 Dune::FieldVector<RF,dim> vu(0.0);
1054 for(
unsigned int d=0; d<
dim; ++d){
1055 const LFSV_V & lfsu_v = lfsv_pfs_v.child(d);
1056 for (
size_t i=0; i<lfsu_v.size(); i++)
1057 vu[d] += x(lfsu_v,i) * phi_v[i];
1060 for(
unsigned int dv=0; dv<
dim; ++dv){
1061 const LFSV_V & lfsv_v = lfsv_pfs_v.child(dv);
1064 Dune::FieldVector<RF,dim> gradu(0.0);
1065 for (
size_t i=0; i<lfsv_v.size(); i++)
1066 gradu.axpy(x(lfsv_v,i),grad_phi_v[i][0]);
1069 for(
unsigned int du=0; du <
dim; ++du){
1070 const LFSV_V & lfsu_v = lfsv_pfs_v.child(du);
1072 for (
size_t i=0; i<vsize; i++)
1073 for(
size_t j=0; j<vsize; j++)
1074 mat.accumulate(lfsv_v,i,lfsu_v,j,
1075 rho * phi_v[j] * gradu[du] * phi_v[i] * weight);
1078 const LFSV_V & lfsu_v = lfsv_pfs_v.child(dv);
1079 for(
size_t j=0; j<vsize; j++){
1080 const Dune::FieldVector<RF,dim> du(grad_phi_v[j][0]);
1081 for (
size_t i=0; i<vsize; i++){
1082 mat.accumulate(lfsv_v,i,lfsu_v,j, rho * (vu * du) * phi_v[i] * weight);
1091 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
1092 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
1095 StokesLocalOperator::alpha_volume(eg,lfsu,x,lfsv,r);
1098 static const unsigned int dim = EG::Geometry::dimension;
1102 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
1103 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
1104 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
1106 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
1109 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
1110 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
1111 const unsigned int vsize = lfsv_v.size();
1114 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
1115 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
1116 typedef typename BasisSwitch_V::DomainField DF;
1117 typedef typename BasisSwitch_V::Range RT;
1118 typedef typename BasisSwitch_V::RangeField RF;
1121 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1122 Dune::GeometryType gt = eg.geometry().type();
1123 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
1124 const int jac_order = gt.isSimplex() ? 0 : 1;
1125 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
1126 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
1129 for (
const auto& ip : rule)
1131 const Dune::FieldVector<DF,dim> local = ip.position();
1134 const RF rho = prm.rho(eg,local);
1135 if(rho == 0)
continue;
1138 std::vector<RT> phi_v(vsize);
1139 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1142 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
1143 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1144 eg.geometry(), local, grad_phi_v);
1146 const RF weight = ip.weight() * eg.geometry().integrationElement(ip.position());
1149 Dune::FieldVector<RF,dim> vu(0.0);
1150 for(
unsigned int d=0; d<
dim; ++d){
1151 const LFSV_V & lfsu_v = lfsv_pfs_v.child(d);
1152 for (
size_t i=0; i<lfsu_v.size(); i++)
1153 vu[d] += x(lfsu_v,i) * phi_v[i];
1156 for(
unsigned int d=0; d<
dim; ++d){
1157 const LFSV_V & lfsu_v = lfsv_pfs_v.child(d);
1160 Dune::FieldVector<RF,dim> gradu(0.0);
1161 for (
size_t i=0; i<lfsu_v.size(); i++)
1162 gradu.axpy(x(lfsu_v,i),grad_phi_v[i][0]);
1165 const RF u_nabla_u = vu * gradu;
1167 for (
size_t i=0; i<vsize; i++)
1168 r.accumulate(lfsu_v,i, rho * u_nabla_u * phi_v[i] * weight);
1180 template<
typename PRM>
1188 typedef typename PRM::Traits::RangeField RF;
1193 enum { doPatternVolume =
true };
1196 enum { doAlphaVolume =
true };
1210 DUNE_DEPRECATED_MSG(
"Deprecated in DUNE-PDELab 2.4, use NavierStokesMass instead!")
1212 prm(_prm), superintegration_order(_superintegration_order)
1216 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
1222 static const unsigned int dim = EG::Geometry::dimension;
1225 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
1226 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
1228 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesMassDG");
1231 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
1232 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
1233 const unsigned int vsize = lfsv_v.size();
1236 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
1237 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
1238 typedef typename BasisSwitch_V::DomainField DF;
1239 typedef typename BasisSwitch_V::Range RT;
1240 typedef typename BasisSwitch_V::RangeField RF;
1241 typedef typename LFSV::Traits::SizeType size_type;
1244 Dune::GeometryType gt = eg.geometry().type();
1245 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1246 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
1247 const int qorder = 2*v_order + det_jac_order + superintegration_order;
1248 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
1251 for (
const auto& ip : rule)
1253 const Dune::FieldVector<DF,dim> local = ip.position();
1256 std::vector<RT> psi_v(vsize);
1257 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,psi_v);
1259 const RF rho = prm.rho(eg,local);
1260 const RF weight = ip.weight() * eg.geometry().integrationElement(ip.position());
1265 const RF factor = rho * weight;
1266 for (size_type j=0; j<vsize; j++)
1268 for (size_type i=0; i<vsize; i++)
1270 const RF val = (psi_v[j]*psi_v[i])*factor;
1272 for (
unsigned int d=0; d<
dim; d++)
1274 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
1275 mat.accumulate(lfsv_v,i,lfsv_v,j, val);
void setTime(Real t)
Definition: stokesdg.hh:92
int superintegration_order
Definition: stokesdg.hh:1284
Implement alpha_volume() based on jacobian_volume()
Definition: defaultimp.hh:611
Real current_dt
Definition: stokesdg.hh:960
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: stokesdg.hh:101
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: stokesdg.hh:1092
Implement alpha_boundary() based on jacobian_boundary()
Definition: defaultimp.hh:704
StokesBoundaryCondition BC
Boundary condition indicator type.
Definition: stokesdg.hh:977
A dense matrix for storing data associated with the degrees of freedom of a pair of LocalFunctionSpac...
Definition: localmatrix.hh:184
A local operator for solving the stokes equation using a DG discretization.
Definition: stokesdg.hh:36
const IG & ig
Definition: constraints.hh:147
void jacobian_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: stokesdg.hh:748
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: stokesdg.hh:191
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
Definition: stokesparameter.hh:32
PRM & prm
Definition: stokesdg.hh:958
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: stokesdg.hh:995
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: stokesdg.hh:449
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: stokesdg.hh:338
StokesDG< PRM, full_tensor > StokesLocalOperator
Definition: stokesdg.hh:981
Definition: adaptivity.hh:27
void preStep(RealType, RealType dt, int)
Definition: stokesdg.hh:86
PRM::Traits::RangeField RF
Common range field type.
Definition: stokesdg.hh:979
A local operator for the mass term corresponding to the instationary part in the Navier-Stokes equati...
Definition: navierstokesmass.hh:25
static const int dim
Definition: adaptivity.hh:83
A local operator for solving the navier stokes equation using a DG discretization.
Definition: stokesdg.hh:973
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: stokesdg.hh:1218
sparsity pattern generator
Definition: pattern.hh:29
PRM & prm
Definition: stokesdg.hh:1283
double RealType
Definition: idefault.hh:92
A local operator for solving the stokes equation using a DG discretization.
Definition: stokesdg.hh:1181
A local operator for solving the Navier-Stokes equations using a DG discretization.
Definition: dgnavierstokes.hh:36
Compile-time switch for the boundary slip factor.
Definition: dgnavierstokesparameter.hh:159
Implement alpha_skeleton() based on jacobian_skeleton()
Definition: defaultimp.hh:650
const EG & eg
Definition: constraints.hh:280
int superintegration_order
Definition: stokesdg.hh:959
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89