4 #ifndef DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKES_HH
5 #define DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKES_HH
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>
25 #define PBLOCK (- VBLOCK + 1)
35 template<
typename PRM>
42 typedef typename PRM::Traits::RangeField RF;
47 static const bool navier = PRM::assemble_navier;
48 static const bool full_tensor = PRM::assemble_full_tensor;
79 prm(_prm), superintegration_order(_superintegration_order),
97 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
98 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
101 const unsigned int dim = EG::Geometry::dimension;
105 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokes");
106 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
107 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
109 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for DGNavierStokes");
112 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
113 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
114 const unsigned int vsize = lfsv_v.size();
115 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
116 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
117 const unsigned int psize = lfsv_p.size();
120 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
121 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
122 typedef typename BasisSwitch_V::DomainField DF;
123 typedef typename BasisSwitch_V::Range RT;
124 typedef typename BasisSwitch_V::RangeField RF;
125 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
126 typedef typename LFSV::Traits::SizeType size_type;
129 Dune::GeometryType gt = eg.geometry().type();
130 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
131 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
132 const int jac_order = gt.isSimplex() ? 0 : 1;
133 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
134 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
136 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
139 for (
const auto& ip : rule)
141 const Dune::FieldVector<DF,dim> local = ip.position();
142 const RF mu = prm.mu(eg,local);
143 const RF rho = prm.rho(eg,local);
146 std::vector<RT> phi_v(vsize);
147 Dune::FieldVector<RF,dim> vu(0.0);
149 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
151 for(
unsigned int d=0; d<
dim; ++d) {
152 const LFSV_V& lfsu_v = lfsv_pfs_v.child(d);
153 for(size_type i=0; i<vsize; i++)
154 vu[d] += x(lfsu_v,i) * phi_v[i];
159 std::vector<RT> phi_p(psize);
160 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
164 for(size_type i=0; i<psize; i++)
165 p += x(lfsv_p,i) * phi_p[i];
168 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
169 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
170 eg.geometry(), local, grad_phi_v);
173 Dune::FieldMatrix<RF,dim,dim> jacu(0.0);
174 for(
unsigned int d = 0; d<
dim; ++d) {
175 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
176 for(size_type i=0; i<vsize; i++)
177 jacu[d].axpy(x(lfsv_v,i), grad_phi_v[i][0]);
180 const RF detj = eg.geometry().integrationElement(ip.position());
181 const RF weight = ip.weight() * detj;
183 for(
unsigned int d = 0; d<
dim; ++d) {
184 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
186 for(size_type i=0; i<vsize; i++) {
190 r.accumulate(lfsv_v,i, mu * (jacu[d]*grad_phi_v[i][0]) * weight);
194 for(
unsigned int dd = 0; dd<
dim; ++dd)
195 r.accumulate(lfsv_v,i, mu * jacu[dd][d] * grad_phi_v[i][0][dd] * weight);
200 r.accumulate(lfsv_v,i, -p * grad_phi_v[i][0][d] * weight);
207 const RF u_nabla_u = vu * jacu[d];
209 r.accumulate(lfsv_v,i, rho * u_nabla_u * phi_v[i] * weight);
218 for(size_type i=0; i<psize; i++)
219 for(
unsigned int d = 0; d <
dim; ++d)
221 r.accumulate(lfsv_p,i, -jacu[d][d] * phi_p[i] * incomp_scaling * weight);
227 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
233 const unsigned int dim = EG::Geometry::dimension;
237 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokes");
238 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
239 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
241 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for DGNavierStokes");
244 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
245 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
246 const unsigned int vsize = lfsv_v.size();
247 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
248 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
249 const unsigned int psize = lfsv_p.size();
252 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
253 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
254 typedef typename BasisSwitch_V::DomainField DF;
255 typedef typename BasisSwitch_V::Range RT;
256 typedef typename BasisSwitch_V::RangeField RF;
257 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
258 typedef typename LFSV::Traits::SizeType size_type;
261 Dune::GeometryType gt = eg.geometry().type();
262 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
263 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
264 const int jac_order = gt.isSimplex() ? 0 : 1;
265 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
266 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
268 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
271 for (
const auto& ip : rule)
273 const Dune::FieldVector<DF,dim> local = ip.position();
274 const RF mu = prm.mu(eg,local);
275 const RF rho = prm.rho(eg,local);
278 std::vector<RT> phi_v(vsize);
279 Dune::FieldVector<RF,dim> vu(0.0);
281 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
283 for(
unsigned int d=0; d<
dim; ++d) {
284 const LFSV_V& lfsu_v = lfsv_pfs_v.child(d);
285 for(size_type i=0; i<vsize; i++)
286 vu[d] += x(lfsu_v,i) * phi_v[i];
291 std::vector<RT> phi_p(psize);
292 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
295 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
296 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
297 eg.geometry(), local, grad_phi_v);
299 const RF detj = eg.geometry().integrationElement(ip.position());
300 const RF weight = ip.weight() * detj;
302 for(
unsigned int dv = 0; dv<
dim; ++dv) {
303 const LFSV_V& lfsv_v = lfsv_pfs_v.child(dv);
306 Dune::FieldVector<RF,dim> gradu_dv(0.0);
308 for(size_type l=0; l<vsize; ++l)
309 gradu_dv.axpy(x(lfsv_v,l), grad_phi_v[l][0]);
311 for(size_type i=0; i<vsize; i++) {
313 for(size_type j=0; j<vsize; j++) {
317 mat.accumulate(lfsv_v,i,lfsv_v,j, mu * (grad_phi_v[j][0]*grad_phi_v[i][0]) * weight);
321 for(
unsigned int du = 0; du<
dim; ++du) {
322 const LFSV_V& lfsu_v = lfsv_pfs_v.child(du);
323 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (grad_phi_v[j][0][dv]*grad_phi_v[i][0][du]) * weight);
331 for(size_type j=0; j<psize; j++) {
332 mat.accumulate(lfsv_p,j,lfsv_v,i, -phi_p[j] * grad_phi_v[i][0][dv] * incomp_scaling * weight);
333 mat.accumulate(lfsv_v,i,lfsv_p,j, -phi_p[j] * grad_phi_v[i][0][dv] * weight);
342 for(size_type j=0; j<vsize; j++)
343 mat.accumulate(lfsv_v,i,lfsv_v,j, rho * (vu * grad_phi_v[j][0]) * phi_v[i] * weight);
346 for(
unsigned int du = 0; du <
dim; ++du) {
347 const LFSV_V& lfsu_v = lfsv_pfs_v.child(du);
348 for(size_type j=0; j<vsize; j++)
349 mat.accumulate(lfsv_v,i,lfsu_v,j, rho * phi_v[j] * gradu_dv[du] * phi_v[i] * weight);
362 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
364 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
365 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
366 R& r_s, R& r_n)
const
369 const unsigned int dim = IG::Geometry::dimension;
370 const unsigned int dimw = IG::Geometry::dimensionworld;
374 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokes");
376 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
377 const LFSV_PFS_V& lfsv_s_pfs_v = lfsv_s.template child<VBLOCK>();
378 const LFSV_PFS_V& lfsv_n_pfs_v = lfsv_n.template child<VBLOCK>();
380 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for DGNavierStokes");
383 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
384 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.template child<0>();
385 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.template child<0>();
386 const unsigned int vsize_s = lfsv_s_v.size();
387 const unsigned int vsize_n = lfsv_n_v.size();
388 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
389 const LFSV_P& lfsv_s_p = lfsv_s.template child<PBLOCK>();
390 const LFSV_P& lfsv_n_p = lfsv_n.template child<PBLOCK>();
391 const unsigned int psize_s = lfsv_s_p.size();
392 const unsigned int psize_n = lfsv_n_p.size();
395 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
396 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
397 typedef typename BasisSwitch_V::DomainField DF;
398 typedef typename BasisSwitch_V::Range RT;
399 typedef typename BasisSwitch_V::RangeField RF;
400 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
403 auto inside_cell = ig.inside();
404 auto outside_cell = ig.outside();
407 Dune::GeometryType gtface = ig.geometry().type();
408 const int v_order = FESwitch_V::basis(lfsv_s_v.finiteElement()).order();
409 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
410 const int qorder = 2*v_order + det_jac_order + superintegration_order;
411 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
413 const int epsilon = prm.epsilonIPSymmetryFactor();
414 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
417 for (
const auto& ip : rule)
421 Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(ip.position());
422 Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(ip.position());
424 const RF penalty_factor = prm.getFaceIP(ig,ip.position());
427 std::vector<RT> phi_v_s(vsize_s);
428 std::vector<RT> phi_v_n(vsize_n);
429 FESwitch_V::basis(lfsv_s_v.finiteElement()).evaluateFunction(local_s,phi_v_s);
430 FESwitch_V::basis(lfsv_n_v.finiteElement()).evaluateFunction(local_n,phi_v_n);
433 assert(vsize_s == vsize_n);
434 Dune::FieldVector<RF,dim> u_s(0.0), u_n(0.0);
435 for(
unsigned int d=0; d<
dim; ++d) {
436 const LFSV_V & lfsv_s_v = lfsv_s_pfs_v.child(d);
437 const LFSV_V & lfsv_n_v = lfsv_n_pfs_v.child(d);
438 for(
unsigned int i=0; i<vsize_s; i++) {
439 u_s[d] += x_s(lfsv_s_v,i) * phi_v_s[i];
440 u_n[d] += x_n(lfsv_n_v,i) * phi_v_n[i];
445 std::vector<RT> phi_p_s(psize_s);
446 std::vector<RT> phi_p_n(psize_n);
447 FESwitch_P::basis(lfsv_s_p.finiteElement()).evaluateFunction(local_s,phi_p_s);
448 FESwitch_P::basis(lfsv_n_p.finiteElement()).evaluateFunction(local_n,phi_p_n);
451 assert(psize_s == psize_n);
452 RF p_s(0.0), p_n(0.0);
453 for(
unsigned int i=0; i<psize_s; i++) {
454 p_s += x_s(lfsv_s_p,i) * phi_p_s[i];
455 p_n += x_n(lfsv_n_p,i) * phi_p_n[i];
459 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_s(vsize_s);
460 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_s_v.finiteElement()),
461 inside_cell.geometry(), local_s, grad_phi_v_s);
463 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_n(vsize_n);
464 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_n_v.finiteElement()),
465 outside_cell.geometry(), local_n, grad_phi_v_n);
468 Dune::FieldMatrix<RF,dim,dim> jacu_s(0.0), jacu_n(0.0);
469 for(
unsigned int d=0; d<
dim; ++d) {
470 const LFSV_V & lfsv_s_v = lfsv_s_pfs_v.child(d);
471 const LFSV_V & lfsv_n_v = lfsv_n_pfs_v.child(d);
472 for(
unsigned int i=0; i<vsize_s; i++) {
473 jacu_s[d].axpy(x_s(lfsv_s_v,i), grad_phi_v_s[i][0]);
474 jacu_n[d].axpy(x_n(lfsv_n_v,i), grad_phi_v_n[i][0]);
478 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
479 const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
480 const RF mu = prm.mu(ig,ip.position());
482 const RF factor = mu * weight;
484 for(
unsigned int d=0; d<
dim; ++d) {
485 const LFSV_V & lfsv_s_v = lfsv_s_pfs_v.child(d);
486 const LFSV_V & lfsv_n_v = lfsv_n_pfs_v.child(d);
491 RF val = 0.5 * ((jacu_s[d] * normal) + (jacu_n[d] * normal)) * factor;
492 for(
unsigned int i=0; i<vsize_s; i++) {
493 r_s.accumulate(lfsv_s_v,i, -val * phi_v_s[i]);
494 r_n.accumulate(lfsv_n_v,i, val * phi_v_n[i]);
497 for(
unsigned int dd=0; dd<
dim; ++dd) {
498 RF Tval = 0.5 * (jacu_s[dd][d] + jacu_n[dd][d]) * normal[dd] * factor;
499 r_s.accumulate(lfsv_s_v,i, -Tval * phi_v_s[i]);
500 r_n.accumulate(lfsv_n_v,i, Tval * phi_v_n[i]);
508 RF jumpu_d = u_s[d] - u_n[d];
509 for(
unsigned int i=0; i<vsize_s; i++) {
510 r_s.accumulate(lfsv_s_v,i, epsilon * 0.5 * (grad_phi_v_s[i][0] * normal) * jumpu_d * factor);
511 r_n.accumulate(lfsv_n_v,i, epsilon * 0.5 * (grad_phi_v_n[i][0] * normal) * jumpu_d * factor);
514 for(
unsigned int dd=0; dd<
dim; ++dd) {
515 r_s.accumulate(lfsv_s_v,i, epsilon * 0.5 * grad_phi_v_s[i][0][dd] * (u_s[dd] - u_n[dd]) * normal[d] * factor);
516 r_n.accumulate(lfsv_n_v,i, epsilon * 0.5 * grad_phi_v_n[i][0][dd] * (u_s[dd] - u_n[dd]) * normal[d] * factor);
524 for(
unsigned int i=0; i<vsize_s; i++) {
525 r_s.accumulate(lfsv_s_v,i, penalty_factor * jumpu_d * phi_v_s[i] * weight);
526 r_n.accumulate(lfsv_n_v,i, -penalty_factor * jumpu_d * phi_v_n[i] * weight);
532 RF mean_p = 0.5*(p_s + p_n);
533 for(
unsigned int i=0; i<vsize_s; i++) {
534 r_s.accumulate(lfsv_s_v,i, mean_p * phi_v_s[i] * normal[d] * weight);
535 r_n.accumulate(lfsv_n_v,i, -mean_p * phi_v_n[i] * normal[d] * weight);
542 RF jumpu_n = (u_s*normal) - (u_n*normal);
543 for(
unsigned int i=0; i<psize_s; i++) {
544 r_s.accumulate(lfsv_s_p,i, 0.5 * phi_p_s[i] * jumpu_n * incomp_scaling * weight);
545 r_n.accumulate(lfsv_n_p,i, 0.5 * phi_p_n[i] * jumpu_n * incomp_scaling * weight);
552 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
555 const LFSU& lfsu_s,
const X&,
const LFSV& lfsv_s,
556 const LFSU& lfsu_n,
const X&,
const LFSV& lfsv_n,
561 const unsigned int dim = IG::Geometry::dimension;
562 const unsigned int dimw = IG::Geometry::dimensionworld;
566 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokes");
568 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
569 const LFSV_PFS_V& lfsv_s_pfs_v = lfsv_s.template child<VBLOCK>();
570 const LFSV_PFS_V& lfsv_n_pfs_v = lfsv_n.template child<VBLOCK>();
572 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for DGNavierStokes");
575 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
576 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.template child<0>();
577 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.template child<0>();
578 const unsigned int vsize_s = lfsv_s_v.size();
579 const unsigned int vsize_n = lfsv_n_v.size();
580 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
581 const LFSV_P& lfsv_s_p = lfsv_s.template child<PBLOCK>();
582 const LFSV_P& lfsv_n_p = lfsv_n.template child<PBLOCK>();
583 const unsigned int psize_s = lfsv_s_p.size();
584 const unsigned int psize_n = lfsv_n_p.size();
587 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
588 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
589 typedef typename BasisSwitch_V::DomainField DF;
590 typedef typename BasisSwitch_V::Range RT;
591 typedef typename BasisSwitch_V::RangeField RF;
592 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
595 auto inside_cell = ig.inside();
596 auto outside_cell = ig.outside();
599 Dune::GeometryType gtface = ig.geometry().type();
600 const int v_order = FESwitch_V::basis(lfsv_s_v.finiteElement()).order();
601 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
602 const int qorder = 2*v_order + det_jac_order + superintegration_order;
603 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
605 const int epsilon = prm.epsilonIPSymmetryFactor();
606 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
609 for (
const auto& ip : rule)
613 Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(ip.position());
614 Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(ip.position());
616 const RF penalty_factor = prm.getFaceIP(ig,ip.position());
619 std::vector<RT> phi_v_s(vsize_s);
620 std::vector<RT> phi_v_n(vsize_n);
621 FESwitch_V::basis(lfsv_s_v.finiteElement()).evaluateFunction(local_s,phi_v_s);
622 FESwitch_V::basis(lfsv_n_v.finiteElement()).evaluateFunction(local_n,phi_v_n);
624 std::vector<RT> phi_p_s(psize_s);
625 std::vector<RT> phi_p_n(psize_n);
626 FESwitch_P::basis(lfsv_s_p.finiteElement()).evaluateFunction(local_s,phi_p_s);
627 FESwitch_P::basis(lfsv_n_p.finiteElement()).evaluateFunction(local_n,phi_p_n);
630 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_s(vsize_s);
631 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_s_v.finiteElement()),
632 inside_cell.geometry(), local_s, grad_phi_v_s);
634 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_n(vsize_n);
635 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_n_v.finiteElement()),
636 outside_cell.geometry(), local_n, grad_phi_v_n);
638 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
639 const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
640 const RF mu = prm.mu(ig,ip.position());
642 assert(vsize_s == vsize_n);
643 const RF factor = mu * weight;
645 for(
unsigned int d = 0; d <
dim; ++d) {
646 const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
647 const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
653 for(
unsigned int i=0; i<vsize_s; ++i) {
655 for(
unsigned int j=0; j<vsize_s; ++j) {
656 RF val = (0.5*(grad_phi_v_s[i][0]*normal)*phi_v_s[j]) * factor;
657 mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v,i, -val);
658 mat_ss.accumulate(lfsv_s_v,i,lfsv_s_v,j, epsilon * val);
659 mat_ss.accumulate(lfsv_s_v,i,lfsv_s_v,j, phi_v_s[i] * phi_v_s[j] * penalty_factor * weight);
663 for(
unsigned int dd = 0; dd <
dim; ++dd) {
664 RF Tval = (0.5*(grad_phi_v_s[i][0][d]*normal[dd])*phi_v_s[j]) * factor;
665 const LFSV_V& lfsv_s_v_dd = lfsv_s_pfs_v.child(dd);
666 mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v_dd,i, - Tval);
667 mat_ss.accumulate(lfsv_s_v_dd,i,lfsv_s_v,j, epsilon*Tval );
672 for(
unsigned int j=0; j<vsize_n; ++j) {
674 RF val = (-0.5*(grad_phi_v_s[i][0]*normal)*phi_v_n[j]) * factor;
675 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v,i,- val);
676 mat_sn.accumulate(lfsv_s_v,i,lfsv_n_v,j, epsilon*val);
677 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v,i, -phi_v_s[i] * phi_v_n[j] * penalty_factor * weight);
681 for (
unsigned int dd=0;dd<
dim;++dd) {
682 RF Tval = (-0.5*(grad_phi_v_s[i][0][d]*normal[dd])*phi_v_n[j]) * factor;
683 const LFSV_V& lfsv_s_v_dd = lfsv_s_pfs_v.child(dd);
684 mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v_dd,i,- Tval);
685 mat_sn.accumulate(lfsv_s_v_dd,i,lfsv_n_v,j, epsilon*Tval);
690 for(
unsigned int j=0; j<vsize_s; ++j) {
691 RF val = (0.5*(grad_phi_v_n[i][0]*normal)*phi_v_s[j]) * factor;
692 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v,i, - val);
693 mat_ns.accumulate(lfsv_n_v,i,lfsv_s_v,j, epsilon*val );
694 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v,i, -phi_v_n[i] * phi_v_s[j] * penalty_factor * weight);
698 for (
unsigned int dd=0;dd<
dim;++dd) {
699 RF Tval = (0.5*(grad_phi_v_n[i][0][d]*normal[dd])*phi_v_s[j]) * factor;
700 const LFSV_V& lfsv_n_v_dd = lfsv_n_pfs_v.child(dd);
701 mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v_dd,i, - Tval);
702 mat_ns.accumulate(lfsv_n_v_dd,i,lfsv_s_v,j, epsilon*Tval );
707 for(
unsigned int j=0; j<vsize_n; ++j) {
709 RF val = (-0.5*(grad_phi_v_n[i][0]*normal)*phi_v_n[j]) * factor;
710 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v,i, - val);
711 mat_nn.accumulate(lfsv_n_v,i,lfsv_n_v,j, epsilon*val);
712 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v,i, phi_v_n[i] * phi_v_n[j] * penalty_factor * weight);
716 for (
unsigned int dd=0;dd<
dim;++dd) {
717 RF Tval = (-0.5*(grad_phi_v_n[i][0][d]*normal[dd])*phi_v_n[j]) * factor;
718 const LFSV_V& lfsv_n_v_dd = lfsv_n_pfs_v.child(dd);
719 mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v_dd,i,- Tval);
720 mat_nn.accumulate(lfsv_n_v_dd,i,lfsv_n_v,j, epsilon*Tval);
729 for(
unsigned int j=0; j<psize_s; ++j) {
730 RF val = 0.5*(phi_p_s[j]*normal[d]*phi_v_s[i]) * weight;
731 mat_ss.accumulate(lfsv_s_v,i,lfsv_s_p,j, val);
732 mat_ss.accumulate(lfsv_s_p,j,lfsv_s_v,i, val * incomp_scaling);
735 for(
unsigned int j=0; j<psize_n; ++j) {
736 RF val = 0.5*(phi_p_n[j]*normal[d]*phi_v_s[i]) * weight;
737 mat_sn.accumulate(lfsv_s_v,i,lfsv_n_p,j, val);
738 mat_ns.accumulate(lfsv_n_p,j,lfsv_s_v,i, val * incomp_scaling);
741 for (
unsigned int j=0; j<psize_s;++j) {
743 RF val = -0.5*(phi_p_s[j]*normal[d]*phi_v_n[i]) * weight;
744 mat_ns.accumulate(lfsv_n_v,i,lfsv_s_p,j, val);
745 mat_sn.accumulate(lfsv_s_p,j,lfsv_n_v,i, val * incomp_scaling);
748 for (
unsigned int j=0; j<psize_n;++j) {
750 RF val = -0.5*(phi_p_n[j]*normal[d]*phi_v_n[i]) * weight;
751 mat_nn.accumulate(lfsv_n_v,i,lfsv_n_p,j, val);
752 mat_nn.accumulate(lfsv_n_p,j,lfsv_n_v,i, val * incomp_scaling);
761 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
763 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
767 const unsigned int dim = IG::Geometry::dimension;
768 const unsigned int dimw = IG::Geometry::dimensionworld;
772 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokes");
774 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
775 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
777 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for DGNavierStokes");
780 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
781 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
782 const unsigned int vsize = lfsv_v.size();
783 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
784 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
785 const unsigned int psize = lfsv_p.size();
788 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
789 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
790 typedef typename BasisSwitch_V::DomainField DF;
791 typedef typename BasisSwitch_V::Range RT;
792 typedef typename BasisSwitch_V::RangeField RF;
793 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
796 auto inside_cell = ig.inside();
799 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
800 Dune::GeometryType gtface = ig.geometry().type();
801 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-1);
802 const int qorder = 2*v_order + det_jac_order + superintegration_order;
803 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
805 const int epsilon = prm.epsilonIPSymmetryFactor();
806 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
809 for (
const auto& ip : rule)
812 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
814 const RF penalty_factor = prm.getFaceIP(ig,ip.position() );
817 std::vector<RT> phi_v(vsize);
818 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
821 Dune::FieldVector<RF,dim> u(0.0);
822 for(
unsigned int d=0; d<
dim; ++d) {
823 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
824 for(
unsigned int i=0; i<vsize; i++)
825 u[d] += x(lfsv_v,i) * phi_v[i];
829 std::vector<RT> phi_p(psize);
830 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
834 for(
unsigned int i=0; i<psize; i++)
835 p += x(lfsv_p,i) * phi_p[i];
838 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
839 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
840 inside_cell.geometry(), local, grad_phi_v);
843 Dune::FieldMatrix<RF,dim,dim> jacu(0.0);
844 for(
unsigned int d=0; d<
dim; ++d) {
845 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
846 for(
unsigned int i=0; i<vsize; i++)
847 jacu[d].axpy(x(lfsv_v,i), grad_phi_v[i][0]);
850 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
851 const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
852 const RF mu = prm.mu(ig,ip.position());
855 typename PRM::Traits::BoundaryCondition::Type bctype(prm.bctype(ig,ip.position()));
858 RF slip_factor = 0.0;
865 slip_factor = BoundarySlipSwitch::boundarySlip(prm,ig,ip.position());
871 const RF factor = weight * (1.0 - slip_factor);
873 for(
unsigned int d = 0; d <
dim; ++d) {
874 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
876 RF val = (jacu[d] * normal) * factor * mu;
877 for(
unsigned int i=0; i<vsize; i++) {
881 r.accumulate(lfsv_v,i, -val * phi_v[i]);
882 r.accumulate(lfsv_v,i, epsilon * mu * (grad_phi_v[i][0] * normal) * u[d] * factor);
885 for(
unsigned int dd=0; dd<
dim; ++dd) {
886 r.accumulate(lfsv_v,i, -mu * jacu[dd][d] * normal[dd] * phi_v[i] * factor);
887 r.accumulate(lfsv_v,i, epsilon * mu * grad_phi_v[i][0][dd] * u[dd] * normal[d] * factor);
894 r.accumulate(lfsv_v,i, u[d] * phi_v[i] * penalty_factor * factor);
899 r.accumulate(lfsv_v,i, p * phi_v[i] * normal[d] * weight);
906 for(
unsigned int i=0; i<psize; i++) {
907 r.accumulate(lfsv_p,i, phi_p[i] * (u * normal) * incomp_scaling * weight);
913 const RF factor = weight * (slip_factor);
919 for(
unsigned int d = 0; d <
dim; ++d) {
920 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
922 for(
unsigned int i=0; i<vsize; i++) {
926 for(
unsigned int dd = 0; dd <
dim; ++dd)
927 r.accumulate(lfsv_v,i, -ten_sum * mu * (jacu[dd] * normal) * normal[dd] *phi_v[i] * normal[d] * factor);
928 r.accumulate(lfsv_v,i, epsilon * ten_sum * mu * (u * normal) * (grad_phi_v[i][0] * normal) * normal[d] * factor);
933 r.accumulate(lfsv_v,i, penalty_factor * (u * normal) * phi_v[i] * normal[d] * factor);
942 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
945 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
949 const unsigned int dim = IG::Geometry::dimension;
950 const unsigned int dimw = IG::Geometry::dimensionworld;
954 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokes");
956 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
957 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
959 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for DGNavierStokes");
962 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
963 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
964 const unsigned int vsize = lfsv_v.size();
965 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
966 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
967 const unsigned int psize = lfsv_p.size();
970 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
971 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
972 typedef typename BasisSwitch_V::DomainField DF;
973 typedef typename BasisSwitch_V::Range RT;
974 typedef typename BasisSwitch_V::RangeField RF;
975 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
978 auto inside_cell = ig.inside();
981 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
982 Dune::GeometryType gtface = ig.geometry().type();
983 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-1);
984 const int qorder = 2*v_order + det_jac_order + superintegration_order;
985 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
987 const int epsilon = prm.epsilonIPSymmetryFactor();
988 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
991 for (
const auto& ip : rule)
994 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
996 const RF penalty_factor = prm.getFaceIP(ig,ip.position() );
999 std::vector<RT> phi_v(vsize);
1000 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1002 std::vector<RT> phi_p(psize);
1003 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1005 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
1006 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1007 inside_cell.geometry(), local, grad_phi_v);
1009 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
1010 const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
1011 const RF mu = prm.mu(ig,ip.position());
1014 typename PRM::Traits::BoundaryCondition::Type bctype(prm.bctype(ig,ip.position()));
1017 RF slip_factor = 0.0;
1024 slip_factor = BoundarySlipSwitch::boundarySlip(prm,ig,ip.position());
1030 const RF factor = weight * (1.0 - slip_factor);
1032 for(
unsigned int d = 0; d <
dim; ++d) {
1033 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
1035 for(
unsigned int i=0; i<vsize; i++) {
1037 for(
unsigned int j=0; j<vsize; j++) {
1041 RF val = ((grad_phi_v[j][0]*normal)*phi_v[i]) * factor * mu;
1042 mat.accumulate(lfsv_v,i,lfsv_v,j, - val);
1043 mat.accumulate(lfsv_v,j,lfsv_v,i, epsilon*val);
1047 for(
unsigned int dd = 0; dd <
dim; ++dd) {
1048 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
1049 RF Tval = ((grad_phi_v[j][0][d]*normal[dd])*phi_v[i]) * factor * mu;
1050 mat.accumulate(lfsv_v,i,lfsv_v_dd,j, - Tval);
1051 mat.accumulate(lfsv_v_dd,j,lfsv_v,i, epsilon*Tval);
1057 mat.accumulate(lfsv_v,j,lfsv_v,i, phi_v[i] * phi_v[j] * penalty_factor * factor);
1064 for(
unsigned int j=0; j<psize; j++) {
1065 mat.accumulate(lfsv_p,j,lfsv_v,i, (phi_p[j]*normal[d]*phi_v[i]) * weight * incomp_scaling);
1066 mat.accumulate(lfsv_v,i,lfsv_p,j, (phi_p[j]*normal[d]*phi_v[i]) * weight);
1075 const RF factor = weight * (slip_factor);
1081 for (
unsigned int i=0;i<vsize;++i)
1083 for (
unsigned int j=0;j<vsize;++j)
1091 RF val = ten_sum * ((grad_phi_v[j][0]*normal)*phi_v[i]) * factor * mu;
1092 for (
unsigned int d=0;d<
dim;++d)
1094 const LFSV_V& lfsv_v_d = lfsv_pfs_v.child(d);
1096 for (
unsigned int dd=0;dd<
dim;++dd)
1098 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
1100 mat.accumulate(lfsv_v_dd,i,lfsv_v_d,j, -val*normal[d]*normal[dd]);
1101 mat.accumulate(lfsv_v_d,j,lfsv_v_dd,i, epsilon*val*normal[d]*normal[dd]);
1110 const RF p_factor = penalty_factor * factor;
1111 for (
unsigned int i=0;i<vsize;++i)
1113 for (
unsigned int j=0;j<vsize;++j)
1115 RF val = phi_v[i]*phi_v[j] * p_factor;
1116 for (
unsigned int d=0;d<
dim;++d)
1118 const LFSV_V& lfsv_v_d = lfsv_pfs_v.child(d);
1119 for (
unsigned int dd=0;dd<
dim;++dd)
1121 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
1122 mat.accumulate(lfsv_v_d,j,lfsv_v_dd,i, val*normal[d]*normal[dd]);
1134 template<
typename EG,
typename LFSV,
typename R>
1138 static const unsigned int dim = EG::Geometry::dimension;
1142 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokes");
1144 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
1145 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
1148 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for DGNavierStokes");
1151 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
1152 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
1153 const unsigned int vsize = lfsv_v.size();
1154 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
1155 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
1156 const unsigned int psize = lfsv_p.size();
1159 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
1160 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
1161 typedef typename BasisSwitch_V::DomainField DF;
1162 typedef typename BasisSwitch_V::Range RT;
1163 typedef typename BasisSwitch_V::RangeField RF;
1164 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
1165 typedef typename LFSV::Traits::SizeType size_type;
1168 Dune::GeometryType gt = eg.geometry().type();
1169 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1170 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
1171 const int qorder = 2*v_order + det_jac_order + superintegration_order;
1173 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
1176 for (
const auto& ip : rule)
1178 const Dune::FieldVector<DF,dim> local = ip.position();
1182 std::vector<RT> phi_v(vsize);
1183 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1186 std::vector<RT> phi_p(psize);
1187 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1189 const RF weight = ip.weight() * eg.geometry().integrationElement(ip.position());
1192 typename PRM::Traits::VelocityRange fval(prm.f(eg,local));
1197 const RF factor = weight;
1198 for (
unsigned int d=0; d<
dim; d++) {
1199 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
1202 for (size_type i=0; i<vsize; i++) {
1203 RF val = phi_v[i]*factor;
1204 r.accumulate(lfsv_v,i, -fval[d] * val);
1208 const RF g2 = prm.g2(eg,ip.position());
1211 for (
size_t i=0; i<lfsv_p.size(); i++) {
1212 r.accumulate(lfsv_p,i, g2 * phi_p[i] * factor);
1220 template<
typename IG,
typename LFSV,
typename R>
1224 static const unsigned int dim = IG::Geometry::dimension;
1228 ((LFSV::CHILDREN == 2),
"You seem to use the wrong function space for DGNavierStokes");
1230 typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
1231 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
1234 ((LFSV_PFS_V::CHILDREN == dim),
"You seem to use the wrong function space for DGNavierStokes");
1237 typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
1238 const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
1239 const unsigned int vsize = lfsv_v.size();
1240 typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
1241 const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
1242 const unsigned int psize = lfsv_p.size();
1245 typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
1246 typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
1247 typedef typename BasisSwitch_V::DomainField DF;
1248 typedef typename BasisSwitch_V::Range RT;
1249 typedef typename BasisSwitch_V::RangeField RF;
1250 typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
1253 auto inside_cell = ig.inside();
1256 Dune::GeometryType gtface = ig.geometry().type();
1257 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1258 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
1259 const int jac_order = gtface.isSimplex() ? 0 : 1;
1260 const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
1261 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
1263 const int epsilon = prm.epsilonIPSymmetryFactor();
1264 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
1267 for (
const auto& ip : rule)
1270 Dune::FieldVector<DF,dim-1> flocal = ip.position();
1271 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(flocal);
1274 const RF penalty_factor = prm.getFaceIP(ig,flocal);
1277 std::vector<RT> phi_v(vsize);
1278 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1280 std::vector<RT> phi_p(psize);
1281 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1283 std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
1284 BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1285 inside_cell.geometry(), local, grad_phi_v);
1287 const Dune::FieldVector<DF,dim> normal = ig.unitOuterNormal(ip.position());
1288 const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
1289 const RF mu = prm.mu(ig,flocal);
1292 typename PRM::Traits::BoundaryCondition::Type bctype(prm.bctype(ig,flocal));
1296 typename PRM::Traits::VelocityRange u0(prm.g(inside_cell,local));
1298 RF factor = mu * weight;
1299 for(
unsigned int d = 0; d <
dim; ++d) {
1300 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
1302 for(
unsigned int i=0; i<vsize; i++) {
1306 r.accumulate(lfsv_v,i, -epsilon * (grad_phi_v[i][0] * normal) * factor * u0[d]);
1310 for(
unsigned int dd = 0; dd <
dim; ++dd) {
1311 const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
1312 RF Tval = (grad_phi_v[i][0][d]*normal[dd]) * factor;
1313 r.accumulate(lfsv_v_dd,i, -epsilon * Tval * u0[d]);
1319 r.accumulate(lfsv_v,i, -phi_v[i] * penalty_factor * u0[d] * weight);
1327 for (
unsigned int i=0;i<psize;++i)
1329 RF val = phi_p[i]*(u0 * normal) * weight;
1330 r.accumulate(lfsv_p,i, - val * incomp_scaling);
1335 typename PRM::Traits::VelocityRange stress(prm.j(ig,flocal,normal));
1341 for(
unsigned int d = 0; d <
dim; ++d) {
1342 const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
1344 for(
unsigned int i=0; i<vsize; i++)
1345 r.accumulate(lfsv_v,i, stress[d] * phi_v[i] * weight);
1353 const int superintegration_order;
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: dgnavierstokes.hh:1135
void alpha_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: dgnavierstokes.hh:762
Definition: dgnavierstokes.hh:61
void setTime(doublet_)
set time for subsequent evaluation
Definition: idefault.hh:104
Definition: stokesparameter.hh:36
Definition: stokesparameter.hh:37
A dense matrix for storing data associated with the degrees of freedom of a pair of LocalFunctionSpac...
Definition: localmatrix.hh:184
Definition: stokesparameter.hh:35
DGNavierStokes(PRM &_prm, int _superintegration_order=0)
Constructor.
Definition: dgnavierstokes.hh:78
Definition: dgnavierstokes.hh:57
void setTime(Real t)
Definition: dgnavierstokes.hh:90
const IG & ig
Definition: constraints.hh:147
Definition: dgnavierstokes.hh:60
Definition: dgnavierstokes.hh:62
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
Definition: stokesparameter.hh:32
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: dgnavierstokes.hh:1221
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: dgnavierstokes.hh:554
Definition: dgnavierstokes.hh:54
Definition: adaptivity.hh:27
static const int dim
Definition: adaptivity.hh:83
Definition: dgnavierstokes.hh:63
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: dgnavierstokes.hh:363
const P & p
Definition: constraints.hh:146
void jacobian_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: dgnavierstokes.hh:944
sparsity pattern generator
Definition: pattern.hh:29
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: dgnavierstokes.hh:229
double RealType
Definition: idefault.hh:92
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: dgnavierstokes.hh:98
Definition: dgnavierstokes.hh:64
A local operator for solving the Navier-Stokes equations using a DG discretization.
Definition: dgnavierstokes.hh:36
Definition: dgnavierstokes.hh:53
Compile-time switch for the boundary slip factor.
Definition: dgnavierstokesparameter.hh:159
void preStep(RealType, RealType dt, int)
Definition: dgnavierstokes.hh:84
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89