2 #ifndef DUNE_PDELAB_STOKESDG_HH
3 #define DUNE_PDELAB_STOKESDG_HH
5 #include <dune/common/exceptions.hh>
6 #include <dune/common/fvector.hh>
7 #include <dune/common/static_assert.hh>
8 #include <dune/common/parametertreeparser.hh>
10 #include <dune/geometry/quadraturerules.hh>
12 #include <dune/localfunctions/common/interfaceswitch.hh>
23 #define PBLOCK (- VBLOCK + 1)
34 template<
typename PRM,
bool full_tensor = true>
45 typedef typename PRM::Traits::RangeField RF;
78 StokesDG (PRM & _prm,
int _superintegration_order=0) :
92 template<
typename EG,
typename LFSV,
typename R>
96 static const unsigned int dim = EG::Geometry::dimension;
100 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
102 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
103 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
106 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
109 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
110 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
111 const unsigned int vsize = lfsv_v.size();
112 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
113 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
114 const unsigned int psize = lfsv_p.size();
117 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
118 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
119 typedef typename BasisSwitch_V::DomainField DF;
120 typedef typename BasisSwitch_V::Range RT;
121 typedef typename BasisSwitch_V::RangeField RF;
122 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
123 typedef typename LFSV::Traits::SizeType size_type;
126 Dune::GeometryType gt = eg.geometry().type();
127 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
128 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
132 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
135 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
137 const Dune::FieldVector<DF,dim> local = it->position();
141 std::vector<RT> phi_v(vsize);
142 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
145 std::vector<RT> phi_p(psize);
146 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
148 const RF weight = it->weight() * eg.geometry().integrationElement(it->position());
151 typename PRM::Traits::VelocityRange fval(
prm.f(eg,local));
156 const RF factor = weight;
157 for (
unsigned int d=0; d<
dim; d++)
159 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
162 for (size_type i=0; i<vsize; i++)
164 RF val = phi_v[i]*factor;
165 r.accumulate(lfsv_v,i, fval[d] * val);
169 const RF g2 =
prm.g2(eg,it->position());
172 for (
size_t i=0; i<lfsv_p.size(); i++)
174 r.accumulate(lfsv_p,i, g2 * phi_p[i] * factor);
182 template<
typename IG,
typename LFSV,
typename R>
186 static const unsigned int dim = IG::Geometry::dimension;
190 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
192 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
193 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
196 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
199 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
200 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
201 const unsigned int vsize = lfsv_v.size();
202 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
203 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
204 const unsigned int psize = lfsv_p.size();
207 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
208 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
209 typedef typename BasisSwitch_V::DomainField DF;
210 typedef typename BasisSwitch_V::Range RT;
211 typedef typename BasisSwitch_V::RangeField RF;
212 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
215 Dune::GeometryType gtface = ig.geometry().type();
216 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
217 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
218 const int jac_order = gtface.isSimplex() ? 0 : 1;
220 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
222 const int epsilon =
prm.epsilonIPSymmetryFactor();
223 const RF incomp_scaling =
prm.incompressibilityScaling(
current_dt);
226 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
229 Dune::FieldVector<DF,dim-1> flocal = it->position();
230 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(flocal);
233 const RF penalty_factor =
prm.getFaceIP(ig,flocal);
236 std::vector<RT> phi_v(vsize);
237 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
239 std::vector<RT> phi_p(psize);
240 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
242 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
243 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
244 ig.inside()->geometry(), local, grad_phi_v);
246 const Dune::FieldVector<DF,dim> normal = ig.unitOuterNormal(it->position());
247 const RF weight = it->weight()*ig.geometry().integrationElement(it->position());
248 const RF mu =
prm.mu(ig,flocal);
251 typename PRM::Traits::BoundaryCondition::Type bctype(
prm.bctype(ig,flocal));
255 typename PRM::Traits::VelocityRange u0(
prm.g(ig,flocal));
260 RF factor = mu * weight;
261 for (
unsigned int i=0;i<vsize;++i)
263 const RF val = (grad_phi_v[i][0]*normal) * factor;
264 for (
unsigned int d=0;d<
dim;++d)
266 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
267 r.accumulate(lfsv_v,i, - epsilon * val * u0[d]);
272 for (
unsigned int dd=0;dd<
dim;++dd)
274 RF Tval = (grad_phi_v[i][0][d]*normal[dd]) * factor;
275 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
276 r.accumulate(lfsv_v_dd,i, - epsilon * Tval * u0[d]);
284 factor = penalty_factor * weight;
285 for (
unsigned int i=0;i<vsize;++i)
287 const RF val = phi_v[i] * factor;
288 for (
unsigned int d=0;d<
dim;++d)
290 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
291 r.accumulate(lfsv_v,i, -val * u0[d] );
297 for (
unsigned int i=0;i<psize;++i)
299 RF val = phi_p[i]*(u0 * normal) * weight;
300 r.accumulate(lfsv_p,i, - val * incomp_scaling);
305 typename PRM::Traits::VelocityRange stress(
prm.j(ig,flocal,normal));
311 for (
unsigned int i=0;i<vsize;++i)
313 for (
unsigned int d=0;d<
dim;++d)
315 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
316 RF val = stress[d]*phi_v[i] * weight;
317 r.accumulate(lfsv_v,i, val);
325 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
331 static const unsigned int dim = EG::Geometry::dimension;
335 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
336 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
337 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
339 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
342 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
343 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
344 const unsigned int vsize = lfsv_v.size();
345 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
346 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
347 const unsigned int psize = lfsv_p.size();
350 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
351 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
352 typedef typename BasisSwitch_V::DomainField DF;
353 typedef typename BasisSwitch_V::Range RT;
354 typedef typename BasisSwitch_V::RangeField RF;
355 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
356 typedef typename LFSV::Traits::SizeType size_type;
359 Dune::GeometryType gt = eg.geometry().type();
360 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
361 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
362 const int jac_order = gt.isSimplex() ? 0 : 1;
364 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
366 const RF incomp_scaling =
prm.incompressibilityScaling(
current_dt);
369 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
371 const Dune::FieldVector<DF,dim> local = it->position();
372 const RF mu =
prm.mu(eg,local);
375 std::vector<RT> phi_p(psize);
376 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
379 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
380 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
381 eg.geometry(), local, grad_phi_v);
383 const RF detj = eg.geometry().integrationElement(it->position());
384 const RF weight = it->weight() * detj;
389 const RF factor = mu * weight;
390 for (size_type j=0; j<vsize; j++)
392 for (size_type i=0; i<vsize; i++)
395 RF val = (grad_phi_v[j][0]*grad_phi_v[i][0])*factor;
397 for (
unsigned int d=0; d<
dim; d++)
399 const LFSV_V& lfsv_v_d = lfsv_pfs_v.child(d);
400 mat.accumulate(lfsv_v_d,i,lfsv_v_d,j, val);
404 for (
unsigned int dd=0; dd<
dim; dd++){
405 RF Tval = (grad_phi_v[j][0][d]*grad_phi_v[i][0][dd])*factor;
406 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
407 mat.accumulate(lfsv_v_d,i,lfsv_v_dd,j, Tval);
419 for (size_type j=0; j<psize; j++)
421 RF val = -1.0 * phi_p[j]*weight;
422 for (size_type i=0; i<vsize; i++)
424 for (
unsigned int d=0; d<
dim; d++)
426 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
427 mat.accumulate(lfsv_p,j,lfsv_v,i, val*grad_phi_v[i][0][d] * incomp_scaling);
428 mat.accumulate(lfsv_v,i,lfsv_p,j, val*grad_phi_v[i][0][d]);
436 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
439 const LFSU& lfsu_s,
const X&,
const LFSV& lfsv_s,
440 const LFSU& lfsu_n,
const X&,
const LFSV& lfsv_n,
445 static const unsigned int dim = IG::Geometry::dimension;
446 static const unsigned int dimw = IG::Geometry::dimensionworld;
450 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
452 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
453 const LFSV_PFS_V& lfsv_s_pfs_v = lfsv_s.template child<VBLOCK>();
454 const LFSV_PFS_V& lfsv_n_pfs_v = lfsv_n.template child<VBLOCK>();
456 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
459 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
460 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.template child<0>();
461 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.template child<0>();
462 const unsigned int vsize_s = lfsv_s_v.size();
463 const unsigned int vsize_n = lfsv_n_v.size();
464 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
465 const LFSV_P& lfsv_s_p = lfsv_s.template child<PBLOCK>();
466 const LFSV_P& lfsv_n_p = lfsv_n.template child<PBLOCK>();
467 const unsigned int psize_s = lfsv_s_p.size();
468 const unsigned int psize_n = lfsv_n_p.size();
471 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
472 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
473 typedef typename BasisSwitch_V::DomainField DF;
474 typedef typename BasisSwitch_V::Range RT;
475 typedef typename BasisSwitch_V::RangeField RF;
476 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
479 Dune::GeometryType gtface = ig.geometry().type();
480 const int v_order = FESwitch_V::basis(lfsv_s_v.finiteElement()).order();
481 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
483 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
485 const int epsilon =
prm.epsilonIPSymmetryFactor();
486 const RF incomp_scaling =
prm.incompressibilityScaling(
current_dt);
489 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
493 Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(it->position());
494 Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(it->position());
496 const RF penalty_factor =
prm.getFaceIP(ig,it->position());
499 std::vector<RT> phi_v_s(vsize_s);
500 std::vector<RT> phi_v_n(vsize_n);
501 FESwitch_V::basis(lfsv_s_v.finiteElement()).evaluateFunction(local_s,phi_v_s);
502 FESwitch_V::basis(lfsv_n_v.finiteElement()).evaluateFunction(local_n,phi_v_n);
504 std::vector<RT> phi_p_s(psize_s);
505 std::vector<RT> phi_p_n(psize_n);
506 FESwitch_P::basis(lfsv_s_p.finiteElement()).evaluateFunction(local_s,phi_p_s);
507 FESwitch_P::basis(lfsv_n_p.finiteElement()).evaluateFunction(local_n,phi_p_n);
510 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_s(vsize_s);
511 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_s_v.finiteElement()),
512 ig.inside()->geometry(), local_s, grad_phi_v_s);
514 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_n(vsize_n);
515 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_n_v.finiteElement()),
516 ig.outside()->geometry(), local_n, grad_phi_v_n);
518 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(it->position());
519 const RF weight = it->weight()*ig.geometry().integrationElement(it->position());
520 const RF mu =
prm.mu(ig,it->position());
525 assert(vsize_s == vsize_n);
526 RF factor = mu * weight;
527 for (
unsigned int i=0;i<vsize_s;++i)
529 for (
unsigned int j=0;j<vsize_s;++j)
531 RF val = (0.5*(grad_phi_v_s[i][0]*normal)*phi_v_s[j]) * factor;
532 for (
unsigned int d=0;d<
dim;++d)
534 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
535 mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v,i, - val);
536 mat_ss.accumulate(lfsv_s_v,i,lfsv_s_v,j, epsilon*val );
541 for (
unsigned int dd=0;dd<
dim;++dd)
543 RF Tval = (0.5*(grad_phi_v_s[i][0][d]*normal[dd])*phi_v_s[j]) * factor;
544 const LFSV_V& lfsv_s_v_dd = lfsv_s_pfs_v.child(dd);
545 mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v_dd,i, - Tval);
546 mat_ss.accumulate(lfsv_s_v_dd,i,lfsv_s_v,j, epsilon*Tval );
551 for (
unsigned int j=0;j<vsize_n;++j)
554 RF val = (-0.5*(grad_phi_v_s[i][0]*normal)*phi_v_n[j]) * factor;
555 for (
unsigned int d=0;d<
dim;++d)
557 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
558 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
559 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v,i,- val);
560 mat_sn.accumulate(lfsv_s_v,i,lfsv_n_v,j, epsilon*val);
565 for (
unsigned int dd=0;dd<
dim;++dd)
567 RF Tval = (-0.5*(grad_phi_v_s[i][0][d]*normal[dd])*phi_v_n[j]) * factor;
568 const LFSV_V& lfsv_s_v_dd = lfsv_s_pfs_v.child(dd);
569 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v_dd,i,- Tval);
570 mat_sn.accumulate(lfsv_s_v_dd,i,lfsv_n_v,j, epsilon*Tval);
576 for (
unsigned int i=0;i<vsize_n;++i)
578 for (
unsigned int j=0;j<vsize_s;++j)
580 RF val = (0.5*(grad_phi_v_n[i][0]*normal)*phi_v_s[j]) * factor;
581 for (
unsigned int d=0;d<
dim;++d)
583 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
584 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
585 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v,i, - val);
586 mat_ns.accumulate(lfsv_n_v,i,lfsv_s_v,j, epsilon*val );
591 for (
unsigned int dd=0;dd<
dim;++dd)
593 RF Tval = (0.5*(grad_phi_v_n[i][0][d]*normal[dd])*phi_v_s[j]) * factor;
594 const LFSV_V& lfsv_n_v_dd = lfsv_n_pfs_v.child(dd);
595 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v_dd,i, - Tval);
596 mat_ns.accumulate(lfsv_n_v_dd,i,lfsv_s_v,j, epsilon*Tval );
601 for (
unsigned int j=0;j<vsize_n;++j)
604 RF val = (-0.5*(grad_phi_v_n[i][0]*normal)*phi_v_n[j]) * factor;
605 for (
unsigned int d=0;d<
dim;++d)
607 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
608 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v,i,- val);
609 mat_nn.accumulate(lfsv_n_v,i,lfsv_n_v,j, epsilon*val);
614 for (
unsigned int dd=0;dd<
dim;++dd)
616 RF Tval = (-0.5*(grad_phi_v_n[i][0][d]*normal[dd])*phi_v_n[j]) * factor;
617 const LFSV_V& lfsv_n_v_dd = lfsv_n_pfs_v.child(dd);
618 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v_dd,i,- Tval);
619 mat_nn.accumulate(lfsv_n_v_dd,i,lfsv_n_v,j, epsilon*Tval);
628 factor = penalty_factor * weight;
629 for (
unsigned int i=0;i<vsize_s;++i)
631 for (
unsigned int j=0;j<vsize_s;++j)
633 RF val = phi_v_s[i]*phi_v_s[j] * factor;
634 for (
unsigned int d=0;d<
dim;++d)
636 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
637 mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v,i, val);
640 for (
unsigned int j=0;j<vsize_n;++j)
642 RF val = phi_v_s[i]*phi_v_n[j] * factor;
643 for (
unsigned int d=0;d<
dim;++d)
645 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
646 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
647 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v,i, - val);
651 for (
unsigned int i=0;i<vsize_n;++i)
653 for (
unsigned int j=0;j<vsize_s;++j)
655 RF val = phi_v_n[i]*phi_v_s[j] * factor;
656 for (
unsigned int d=0;d<
dim;++d)
658 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
659 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
660 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v,i, - val);
663 for (
unsigned int j=0;j<vsize_n;++j)
665 RF val = phi_v_n[i]*phi_v_n[j] * factor;
666 for (
unsigned int d=0;d<
dim;++d)
668 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
669 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v,i, val);
677 for (
unsigned int i=0;i<vsize_s;++i)
679 for (
unsigned int j=0;j<psize_s;++j)
681 for (
unsigned int d=0;d<
dim;++d)
683 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
684 RF val = 0.5*(phi_p_s[j]*normal[d]*phi_v_s[i]) * weight;
685 mat_ss.accumulate(lfsv_s_v,i,lfsv_s_p,j, val);
686 mat_ss.accumulate(lfsv_s_p,j,lfsv_s_v,i, val * incomp_scaling);
689 for (
unsigned int j=0;j<psize_n;++j)
691 for (
unsigned int d=0;d<
dim;++d)
693 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
694 RF val = 0.5*(phi_p_n[j]*normal[d]*phi_v_s[i]) * weight;
695 mat_sn.accumulate(lfsv_s_v,i,lfsv_n_p,j, val);
696 mat_ns.accumulate(lfsv_n_p,j,lfsv_s_v,i, val * incomp_scaling);
700 for (
unsigned int i=0;i<vsize_n;++i)
702 for (
unsigned int j=0;j<psize_s;++j)
704 for (
unsigned int d=0;d<
dim;++d)
706 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
709 RF val = -0.5*(phi_p_s[j]*normal[d]*phi_v_n[i]) * weight;
710 mat_ns.accumulate(lfsv_n_v,i,lfsv_s_p,j, val);
711 mat_sn.accumulate(lfsv_s_p,j,lfsv_n_v,i, val * incomp_scaling);
714 for (
unsigned int j=0;j<psize_n;++j)
716 for (
unsigned int d=0;d<
dim;++d)
718 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
721 RF val = -0.5*(phi_p_n[j]*normal[d]*phi_v_n[i]) * weight;
722 mat_nn.accumulate(lfsv_n_v,i,lfsv_n_p,j, val);
723 mat_nn.accumulate(lfsv_n_p,j,lfsv_n_v,i, val * incomp_scaling);
731 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
734 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
738 static const unsigned int dim = IG::Geometry::dimension;
739 static const unsigned int dimw = IG::Geometry::dimensionworld;
743 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
745 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
746 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
748 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
751 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
752 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
753 const unsigned int vsize = lfsv_v.size();
754 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
755 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
756 const unsigned int psize = lfsv_p.size();
759 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
760 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
761 typedef typename BasisSwitch_V::DomainField DF;
762 typedef typename BasisSwitch_V::Range RT;
763 typedef typename BasisSwitch_V::RangeField RF;
764 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
767 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
768 Dune::GeometryType gtface = ig.geometry().type();
769 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-1);
771 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
774 typename PRM::Traits::BoundaryCondition::Type bctype(
prm.bctype(ig,rule.begin()->position()));
776 const int epsilon =
prm.epsilonIPSymmetryFactor();
777 const RF incomp_scaling =
prm.incompressibilityScaling(
current_dt);
780 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
783 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
785 const RF penalty_factor =
prm.getFaceIP(ig,it->position() );
788 std::vector<RT> phi_v(vsize);
789 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
791 std::vector<RT> phi_p(psize);
792 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
794 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
795 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
796 ig.inside()->geometry(), local, grad_phi_v);
798 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(it->position());
799 const RF weight = it->weight()*ig.geometry().integrationElement(it->position());
800 const RF mu =
prm.mu(ig,it->position());
803 RF slip_factor = 0.0;
810 slip_factor = BoundarySlipSwitch::boundarySlip(
prm,ig,it->position());
815 const RF factor = weight * (1.0 - slip_factor);
820 for (
unsigned int i=0;i<vsize;++i)
822 for (
unsigned int j=0;j<vsize;++j)
824 RF val = ((grad_phi_v[j][0]*normal)*phi_v[i]) * factor * mu;
825 for (
unsigned int d=0;d<
dim;++d)
827 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
828 mat.accumulate(lfsv_v,i,lfsv_v,j, - val);
829 mat.accumulate(lfsv_v,j,lfsv_v,i, epsilon*val);
834 for (
unsigned int dd=0;dd<
dim;++dd)
836 RF Tval = ((grad_phi_v[j][0][d]*normal[dd])*phi_v[i]) * factor * mu;
837 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
838 mat.accumulate(lfsv_v,i,lfsv_v_dd,j, - Tval);
839 mat.accumulate(lfsv_v_dd,j,lfsv_v,i, epsilon*Tval);
849 for (
unsigned int i=0;i<vsize;++i)
851 for (
unsigned int j=0;j<psize;++j)
853 for (
unsigned int d=0;d<
dim;++d)
855 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
856 RF val = (phi_p[j]*normal[d]*phi_v[i]) * weight;
857 mat.accumulate(lfsv_p,j,lfsv_v,i, val * incomp_scaling);
858 mat.accumulate(lfsv_v,i,lfsv_p,j, val);
865 const RF p_factor = penalty_factor * factor;
866 for (
unsigned int i=0;i<vsize;++i)
868 for (
unsigned int j=0;j<vsize;++j)
870 RF val = phi_v[i]*phi_v[j] * p_factor;
871 for (
unsigned int d=0;d<
dim;++d)
873 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
874 mat.accumulate(lfsv_v,j,lfsv_v,i, val);
882 const RF factor = weight * (slip_factor);
888 for (
unsigned int i=0;i<vsize;++i)
890 for (
unsigned int j=0;j<vsize;++j)
898 RF val = ten_sum * ((grad_phi_v[j][0]*normal)*phi_v[i]) * factor * mu;
899 for (
unsigned int d=0;d<
dim;++d)
901 const LFSV_V& lfsv_v_d = lfsv_pfs_v.child(d);
903 for (
unsigned int dd=0;dd<
dim;++dd)
905 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
907 mat.accumulate(lfsv_v_dd,i,lfsv_v_d,j, -val*normal[d]*normal[dd]);
908 mat.accumulate(lfsv_v_d,j,lfsv_v_dd,i, epsilon*val*normal[d]*normal[dd]);
917 const RF p_factor = penalty_factor * factor;
918 for (
unsigned int i=0;i<vsize;++i)
920 for (
unsigned int j=0;j<vsize;++j)
922 RF val = phi_v[i]*phi_v[j] * p_factor;
923 for (
unsigned int d=0;d<
dim;++d)
925 const LFSV_V& lfsv_v_d = lfsv_pfs_v.child(d);
926 for (
unsigned int dd=0;dd<
dim;++dd)
928 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
929 mat.accumulate(lfsv_v_d,j,lfsv_v_dd,i, val*normal[d]*normal[dd]);
954 template<
typename PRM,
bool full_tensor = true>
961 typedef typename PRM::Traits::RangeField
RF;
974 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
983 static const unsigned int dim = EG::Geometry::dimension;
987 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
988 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
989 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
991 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
994 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
995 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
996 const unsigned int vsize = lfsv_v.size();
999 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
1000 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
1001 typedef typename BasisSwitch_V::DomainField DF;
1002 typedef typename BasisSwitch_V::Range RT;
1003 typedef typename BasisSwitch_V::RangeField
RF;
1006 Dune::GeometryType gt = eg.geometry().type();
1007 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1008 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
1009 const int jac_order = gt.isSimplex() ? 0 : 1;
1011 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
1014 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
1016 const Dune::FieldVector<DF,dim> local = it->position();
1019 const RF rho =
prm.rho(eg,local);
1020 if(rho == 0)
continue;
1023 std::vector<RT> phi_v(vsize);
1024 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1027 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
1028 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1029 eg.geometry(), local, grad_phi_v);
1031 const RF weight = it->weight() * eg.geometry().integrationElement(it->position());
1034 Dune::FieldVector<RF,dim> vu(0.0);
1035 for(
unsigned int d=0; d<
dim; ++d){
1036 const LFSV_V & lfsu_v = lfsv_pfs_v.child(d);
1037 for (
size_t i=0; i<lfsu_v.size(); i++)
1038 vu[d] += x(lfsu_v,i) * phi_v[i];
1041 for(
unsigned int dv=0; dv<
dim; ++dv){
1042 const LFSV_V & lfsv_v = lfsv_pfs_v.child(dv);
1045 Dune::FieldVector<RF,dim> gradu(0.0);
1046 for (
size_t i=0; i<lfsv_v.size(); i++)
1047 gradu.axpy(x(lfsv_v,i),grad_phi_v[i][0]);
1050 for(
unsigned int du=0; du <
dim; ++du){
1051 const LFSV_V & lfsu_v = lfsv_pfs_v.child(du);
1053 for (
size_t i=0; i<vsize; i++)
1054 for(
size_t j=0; j<vsize; j++)
1055 mat.accumulate(lfsv_v,i,lfsu_v,j,
1056 rho * phi_v[j] * gradu[du] * phi_v[i] * weight);
1059 const LFSV_V & lfsu_v = lfsv_pfs_v.child(dv);
1060 for(
size_t j=0; j<vsize; j++){
1061 const Dune::FieldVector<RF,dim> du(grad_phi_v[j][0]);
1062 for (
size_t i=0; i<vsize; i++){
1063 mat.accumulate(lfsv_v,i,lfsu_v,j, rho * (vu * du) * phi_v[i] * weight);
1072 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
1073 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
1079 static const unsigned int dim = EG::Geometry::dimension;
1083 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for StokesDG");
1084 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
1085 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
1087 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesDG");
1090 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
1091 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
1092 const unsigned int vsize = lfsv_v.size();
1095 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
1096 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
1097 typedef typename BasisSwitch_V::DomainField DF;
1098 typedef typename BasisSwitch_V::Range RT;
1099 typedef typename BasisSwitch_V::RangeField
RF;
1102 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1103 Dune::GeometryType gt = eg.geometry().type();
1104 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
1105 const int jac_order = gt.isSimplex() ? 0 : 1;
1107 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
1110 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
1112 const Dune::FieldVector<DF,dim> local = it->position();
1115 const RF rho =
prm.rho(eg,local);
1116 if(rho == 0)
continue;
1119 std::vector<RT> phi_v(vsize);
1120 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1123 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
1124 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1125 eg.geometry(), local, grad_phi_v);
1127 const RF weight = it->weight() * eg.geometry().integrationElement(it->position());
1130 Dune::FieldVector<RF,dim> vu(0.0);
1131 for(
unsigned int d=0; d<
dim; ++d){
1132 const LFSV_V & lfsu_v = lfsv_pfs_v.child(d);
1133 for (
size_t i=0; i<lfsu_v.size(); i++)
1134 vu[d] += x(lfsu_v,i) * phi_v[i];
1137 for(
unsigned int d=0; d<
dim; ++d){
1138 const LFSV_V & lfsu_v = lfsv_pfs_v.child(d);
1141 Dune::FieldVector<RF,dim> gradu(0.0);
1142 for (
size_t i=0; i<lfsu_v.size(); i++)
1143 gradu.axpy(x(lfsu_v,i),grad_phi_v[i][0]);
1146 const RF u_nabla_u = vu * gradu;
1148 for (
size_t i=0; i<vsize; i++)
1149 r.accumulate(lfsu_v,i, rho * u_nabla_u * phi_v[i] * weight);
1161 template<
typename PRM>
1169 typedef typename PRM::Traits::RangeField RF;
1196 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
1202 static const unsigned int dim = EG::Geometry::dimension;
1205 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
1206 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
1208 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for StokesMassDG");
1211 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
1212 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
1213 const unsigned int vsize = lfsv_v.size();
1216 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
1217 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
1218 typedef typename BasisSwitch_V::DomainField DF;
1219 typedef typename BasisSwitch_V::Range RT;
1220 typedef typename BasisSwitch_V::RangeField RF;
1221 typedef typename LFSV::Traits::SizeType size_type;
1224 Dune::GeometryType gt = eg.geometry().type();
1225 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1226 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
1228 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
1231 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
1233 const Dune::FieldVector<DF,dim> local = it->position();
1236 std::vector<RT> psi_v(vsize);
1237 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,psi_v);
1239 const RF rho =
prm.rho(eg,local);
1240 const RF weight = it->weight() * eg.geometry().integrationElement(it->position());
1245 const RF factor = rho * weight;
1246 for (size_type j=0; j<vsize; j++)
1248 for (size_type i=0; i<vsize; i++)
1250 const RF val = (psi_v[j]*psi_v[i])*factor;
1252 for (
unsigned int d=0; d<
dim; d++)
1254 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
1255 mat.accumulate(lfsv_v,i,lfsv_v,j, val);
Definition: stokesdg.hh:61
Compile-time switch for the boundary slip factor.
Definition: stokesdgparameter.hh:195
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:438
Definition: stokesdg.hh:53
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: stokesdg.hh:183
Definition: stokesdg.hh:64
sparsity pattern generator
Definition: pattern.hh:30
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: stokesdg.hh:1073
int superintegration_order
Definition: stokesdg.hh:1264
Real current_dt
Definition: stokesdg.hh:942
Default flags for all local operators.
Definition: flags.hh:18
A local operator for solving the stokes equation using a DG discretization.
Definition: stokesdg.hh:35
int superintegration_order
Definition: stokesdg.hh:941
Definition: stokesdg.hh:62
A local operator for solving the stokes equation using a DG discretization.
Definition: stokesdg.hh:1162
static const int dim
Definition: adaptivity.hh:82
A local operator for solving the navier stokes equation using a DG discretization.
Definition: stokesdg.hh:955
Implement alpha_volume() based on jacobian_volume()
Definition: defaultimp.hh:611
StokesDG(PRM &_prm, int _superintegration_order=0)
Constructor.
Definition: stokesdg.hh:78
double RealType
Definition: idefault.hh:92
Definition: stokesparameter.hh:32
Definition: stokesdg.hh:54
Definition: stokesdg.hh:57
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: stokesdg.hh:976
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
compute
Definition: defaultimp.hh:619
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: stokesdg.hh:327
Definition: stokesdg.hh:63
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: stokesdg.hh:93
StokesBoundaryCondition BC
Boundary condition indicator type.
Definition: stokesdg.hh:959
A dense matrix for storing data associated with the degrees of freedom of a pair of LocalFunctionSpac...
Definition: common/localmatrix.hh:184
const IG & ig
Definition: common/constraints.hh:146
Definition: stokesdg.hh:1174
StokesDG< PRM, full_tensor > StokesLocalOperator
Definition: stokesdg.hh:963
Definition: stokesdg.hh:1177
void preStep(RealType time, RealType dt, int)
Definition: stokesdg.hh:84
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Implement alpha_skeleton() based on jacobian_skeleton()
Definition: defaultimp.hh:650
void jacobian_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: stokesdg.hh:733
sparsity pattern generator
Definition: pattern.hh:14
Definition: stokesparameter.hh:36
PRM & prm
Definition: stokesdg.hh:940
Definition: stokesdg.hh:60
const EG & eg
Definition: common/constraints.hh:277
Definition: stokesparameter.hh:37
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: stokesdg.hh:1198
Definition: stokesparameter.hh:35
PRM & prm
Definition: stokesdg.hh:1263
PRM::Traits::RangeField RF
Common range field type.
Definition: stokesdg.hh:961
NavierStokesDG(PRM &prm_, int superintegration_order_=0)
Definition: stokesdg.hh:970
Implement alpha_boundary() based on jacobian_boundary()
Definition: defaultimp.hh:704
StokesMassDG(PRM &_prm, int _superintegration_order=0)
Constructor.
Definition: stokesdg.hh:1191