2 #ifndef DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONDG_HH 3 #define DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONDG_HH 5 #include<dune/common/exceptions.hh> 6 #include<dune/common/fvector.hh> 8 #include<dune/geometry/referenceelements.hh> 54 template<
typename T,
typename FiniteElementMap>
62 enum {
dim = T::Traits::GridViewType::dimension };
64 using Real =
typename T::Traits::RangeFieldType;
69 enum { doPatternVolume =
true };
70 enum { doPatternSkeleton =
true };
73 enum { doAlphaVolume =
true };
74 enum { doAlphaSkeleton =
true };
75 enum { doAlphaBoundary =
true };
76 enum { doLambdaVolume =
true };
93 param(param_), method(method_), weights(weights_),
94 alpha(alpha_), intorderadd(intorderadd_), quadrature_factor(2),
103 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
104 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 107 using RF =
typename LFSU::Traits::FiniteElementType::
108 Traits::LocalBasisType::Traits::RangeFieldType;
109 using size_type =
typename LFSU::Traits::SizeType;
112 const int dim = EG::Entity::dimension;
113 const int order = std::max(lfsu.finiteElement().localBasis().order(),
114 lfsv.finiteElement().localBasis().order());
117 const auto& cell = eg.entity();
120 auto geo = eg.geometry();
124 auto localcenter = ref_el.position(0,0);
125 auto A = param.A(cell,localcenter);
128 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
129 std::vector<Dune::FieldVector<RF,dim> > gradpsi(lfsv.size());
130 Dune::FieldVector<RF,dim> gradu(0.0);
131 Dune::FieldVector<RF,dim> Agradu(0.0);
134 typename EG::Geometry::JacobianInverseTransposed jac;
137 auto intorder = intorderadd + quadrature_factor * order;
141 auto& phi = cache[order].evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
142 auto& psi = cache[order].evaluateFunction(ip.position(),lfsv.finiteElement().localBasis());
146 for (size_type i=0; i<lfsu.size(); i++)
147 u += x(lfsu,i)*phi[i];
150 auto& js = cache[order].evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
151 auto& js_v = cache[order].evaluateJacobian(ip.position(),lfsv.finiteElement().localBasis());
154 jac = geo.jacobianInverseTransposed(ip.position());
155 for (size_type i=0; i<lfsu.size(); i++)
156 jac.mv(js[i][0],gradphi[i]);
158 for (size_type i=0; i<lfsv.size(); i++)
159 jac.mv(js_v[i][0],gradpsi[i]);
163 for (size_type i=0; i<lfsu.size(); i++)
164 gradu.axpy(x(lfsu,i),gradphi[i]);
170 auto b = param.b(cell,ip.position());
173 auto c = param.c(cell,ip.position());
176 RF factor = ip.weight() * geo.integrationElement(ip.position());
177 for (size_type i=0; i<lfsv.size(); i++)
178 r.accumulate(lfsv,i,( Agradu*gradpsi[i] - u*(b*gradpsi[i]) + c*u*psi[i] )*factor);
183 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
186 alpha_volume(eg,lfsu,x,lfsv,y);
190 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
195 using RF =
typename LFSU::Traits::FiniteElementType::
196 Traits::LocalBasisType::Traits::RangeFieldType;
197 using size_type =
typename LFSU::Traits::SizeType;
200 const int dim = EG::Entity::dimension;
201 const int order = std::max(lfsu.finiteElement().localBasis().order(),
202 lfsv.finiteElement().localBasis().order());
205 const auto& cell = eg.entity();
208 auto geo = eg.geometry();
212 auto localcenter = ref_el.position(0,0);
213 auto A = param.A(cell,localcenter);
216 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
217 std::vector<Dune::FieldVector<RF,dim> > Agradphi(lfsu.size());
220 typename EG::Geometry::JacobianInverseTransposed jac;
223 auto intorder = intorderadd + quadrature_factor * order;
227 auto& phi = cache[order].evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
230 auto& js = cache[order].evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
233 jac = geo.jacobianInverseTransposed(ip.position());
234 for (size_type i=0; i<lfsu.size(); i++)
236 jac.mv(js[i][0],gradphi[i]);
237 A.mv(gradphi[i],Agradphi[i]);
241 auto b = param.b(cell,ip.position());
244 auto c = param.c(cell,ip.position());
247 auto factor = ip.weight() * geo.integrationElement(ip.position());
248 for (size_type j=0; j<lfsu.size(); j++)
249 for (size_type i=0; i<lfsu.size(); i++)
250 mat.accumulate(lfsu,i,lfsu,j,( Agradphi[j]*gradphi[i] - phi[j]*(b*gradphi[i]) + c*phi[j]*phi[i] )*factor);
256 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
258 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
259 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
260 R& r_s, R& r_n)
const 263 using RF =
typename LFSV::Traits::FiniteElementType::
264 Traits::LocalBasisType::Traits::RangeFieldType;
265 using size_type =
typename LFSV::Traits::SizeType;
268 const int dim = IG::dimension;
269 const int order = std::max(
270 std::max(lfsu_s.finiteElement().localBasis().order(),
271 lfsu_n.finiteElement().localBasis().order()),
272 std::max(lfsv_s.finiteElement().localBasis().order(),
273 lfsv_n.finiteElement().localBasis().order())
277 const auto& cell_inside = ig.inside();
278 const auto& cell_outside = ig.outside();
281 auto geo = ig.geometry();
282 auto geo_inside = cell_inside.geometry();
283 auto geo_outside = cell_outside.geometry();
286 auto geo_in_inside = ig.geometryInInside();
287 auto geo_in_outside = ig.geometryInOutside();
292 auto local_inside = ref_el_inside.position(0,0);
293 auto local_outside = ref_el_outside.position(0,0);
294 auto A_s = param.A(cell_inside,local_inside);
295 auto A_n = param.A(cell_outside,local_outside);
299 auto h_F = std::min(geo_inside.volume(),geo_outside.volume())/geo.volume();
302 auto n_F = ig.centerUnitOuterNormal();
303 Dune::FieldVector<RF,dim> An_F_s;
305 Dune::FieldVector<RF,dim> An_F_n;
311 RF harmonic_average(0.0);
314 RF delta_s = (An_F_s*n_F);
315 RF delta_n = (An_F_n*n_F);
316 omega_s = delta_n/(delta_s+delta_n+1
e-20);
317 omega_n = delta_s/(delta_s+delta_n+1
e-20);
318 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-20);
322 omega_s = omega_n = 0.5;
323 harmonic_average = 1.0;
327 auto order_s = lfsu_s.finiteElement().localBasis().order();
328 auto order_n = lfsu_n.finiteElement().localBasis().order();
329 auto degree = std::max( order_s, order_n );
332 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
335 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
336 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
337 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
338 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_n(lfsv_n.size());
339 Dune::FieldVector<RF,dim> gradu_s(0.0);
340 Dune::FieldVector<RF,dim> gradu_n(0.0);
343 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
346 auto intorder = intorderadd+quadrature_factor*order;
350 auto n_F_local = ig.unitOuterNormal(ip.position());
353 auto iplocal_s = geo_in_inside.global(ip.position());
354 auto iplocal_n = geo_in_outside.global(ip.position());
357 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
358 auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
359 auto& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
360 auto& psi_n = cache[order_n].evaluateFunction(iplocal_n,lfsv_n.finiteElement().localBasis());
364 for (size_type i=0; i<lfsu_s.size(); i++)
365 u_s += x_s(lfsu_s,i)*phi_s[i];
367 for (size_type i=0; i<lfsu_n.size(); i++)
368 u_n += x_n(lfsu_n,i)*phi_n[i];
371 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
372 auto& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
373 auto& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
374 auto& gradpsi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsv_n.finiteElement().localBasis());
377 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
378 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
379 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
380 jac = geo_outside.jacobianInverseTransposed(iplocal_n);
381 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
382 for (size_type i=0; i<lfsv_n.size(); i++) jac.mv(gradpsi_n[i][0],tgradpsi_n[i]);
386 for (size_type i=0; i<lfsu_s.size(); i++)
387 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
389 for (size_type i=0; i<lfsu_n.size(); i++)
390 gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
393 auto b = param.b(cell_inside,iplocal_s);
394 auto normalflux = b*n_F_local;
395 RF omegaup_s, omegaup_n;
408 auto factor = ip.weight()*geo.integrationElement(ip.position());
411 auto term1 = (omegaup_s*u_s + omegaup_n*u_n) * normalflux *factor;
412 for (size_type i=0; i<lfsv_s.size(); i++)
413 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
414 for (size_type i=0; i<lfsv_n.size(); i++)
415 r_n.accumulate(lfsv_n,i,-term1 * psi_n[i]);
418 auto term2 = -(omega_s*(An_F_s*gradu_s) + omega_n*(An_F_n*gradu_n)) * factor;
419 for (size_type i=0; i<lfsv_s.size(); i++)
420 r_s.accumulate(lfsv_s,i,term2 * psi_s[i]);
421 for (size_type i=0; i<lfsv_n.size(); i++)
422 r_n.accumulate(lfsv_n,i,-term2 * psi_n[i]);
425 auto term3 = (u_s-u_n) * factor;
426 for (size_type i=0; i<lfsv_s.size(); i++)
427 r_s.accumulate(lfsv_s,i,term3 * theta * omega_s * (An_F_s*tgradpsi_s[i]));
428 for (size_type i=0; i<lfsv_n.size(); i++)
429 r_n.accumulate(lfsv_n,i,term3 * theta * omega_n * (An_F_n*tgradpsi_n[i]));
432 auto term4 = penalty_factor * (u_s-u_n) * factor;
433 for (size_type i=0; i<lfsv_s.size(); i++)
434 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
435 for (size_type i=0; i<lfsv_n.size(); i++)
436 r_n.accumulate(lfsv_n,i,-term4 * psi_n[i]);
441 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
443 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
444 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
445 Y& y_s, Y& y_n)
const 447 alpha_skeleton(ig,lfsu_s,x_s,lfsv_s,lfsu_n,x_n,lfsv_n,y_s,y_n);
450 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
452 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
453 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
454 M& mat_ss, M& mat_sn,
455 M& mat_ns, M& mat_nn)
const 458 using RF =
typename LFSV::Traits::FiniteElementType::
459 Traits::LocalBasisType::Traits::RangeFieldType;
460 using size_type =
typename LFSV::Traits::SizeType;
463 const int dim = IG::dimension;
464 const int order = std::max(
465 std::max(lfsu_s.finiteElement().localBasis().order(),
466 lfsu_n.finiteElement().localBasis().order()),
467 std::max(lfsv_s.finiteElement().localBasis().order(),
468 lfsv_n.finiteElement().localBasis().order())
472 const auto& cell_inside = ig.inside();
473 const auto& cell_outside = ig.outside();
476 auto geo = ig.geometry();
477 auto geo_inside = cell_inside.geometry();
478 auto geo_outside = cell_outside.geometry();
481 auto geo_in_inside = ig.geometryInInside();
482 auto geo_in_outside = ig.geometryInOutside();
487 auto local_inside = ref_el_inside.position(0,0);
488 auto local_outside = ref_el_outside.position(0,0);
489 auto A_s = param.A(cell_inside,local_inside);
490 auto A_n = param.A(cell_outside,local_outside);
494 auto h_F = std::min(geo_inside.volume(),geo_outside.volume())/geo.volume();
497 auto n_F = ig.centerUnitOuterNormal();
498 Dune::FieldVector<RF,dim> An_F_s;
500 Dune::FieldVector<RF,dim> An_F_n;
506 RF harmonic_average(0.0);
509 RF delta_s = (An_F_s*n_F);
510 RF delta_n = (An_F_n*n_F);
511 omega_s = delta_n/(delta_s+delta_n+1
e-20);
512 omega_n = delta_s/(delta_s+delta_n+1
e-20);
513 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-20);
517 omega_s = omega_n = 0.5;
518 harmonic_average = 1.0;
522 auto order_s = lfsu_s.finiteElement().localBasis().order();
523 auto order_n = lfsu_n.finiteElement().localBasis().order();
524 auto degree = std::max( order_s, order_n );
527 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
530 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
531 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
534 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
537 const int intorder = intorderadd+quadrature_factor*order;
541 auto n_F_local = ig.unitOuterNormal(ip.position());
544 auto iplocal_s = geo_in_inside.global(ip.position());
545 auto iplocal_n = geo_in_outside.global(ip.position());
548 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
549 auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
552 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
553 auto& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
556 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
557 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
558 jac = geo_outside.jacobianInverseTransposed(iplocal_n);
559 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
562 auto b = param.b(cell_inside,iplocal_s);
563 auto normalflux = b*n_F_local;
564 RF omegaup_s, omegaup_n;
577 auto factor = ip.weight() * geo.integrationElement(ip.position());
578 auto ipfactor = penalty_factor * factor;
581 for (size_type j=0; j<lfsu_s.size(); j++) {
582 auto temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
583 for (size_type i=0; i<lfsu_s.size(); i++) {
584 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux *factor * phi_s[i]);
585 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,temp1 * phi_s[i]);
586 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
587 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * ipfactor * phi_s[i]);
590 for (size_type j=0; j<lfsu_n.size(); j++) {
591 auto temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
592 for (size_type i=0; i<lfsu_s.size(); i++) {
593 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,omegaup_n * phi_n[j] * normalflux *factor * phi_s[i]);
594 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,temp1 * phi_s[i]);
595 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
596 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * ipfactor * phi_s[i]);
599 for (size_type j=0; j<lfsu_s.size(); j++) {
600 auto temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
601 for (size_type i=0; i<lfsu_n.size(); i++) {
602 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-omegaup_s * phi_s[j] * normalflux *factor * phi_n[i]);
603 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-temp1 * phi_n[i]);
604 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,phi_s[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
605 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-phi_s[j] * ipfactor * phi_n[i]);
608 for (size_type j=0; j<lfsu_n.size(); j++) {
609 auto temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
610 for (size_type i=0; i<lfsu_n.size(); i++) {
611 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-omegaup_n * phi_n[j] * normalflux *factor * phi_n[i]);
612 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-temp1 * phi_n[i]);
613 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
614 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,phi_n[j] * ipfactor * phi_n[i]);
622 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
624 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
628 using RF =
typename LFSV::Traits::FiniteElementType::
629 Traits::LocalBasisType::Traits::RangeFieldType;
630 using size_type =
typename LFSV::Traits::SizeType;
633 const int dim = IG::dimension;
634 const int order = std::max(
635 lfsu_s.finiteElement().localBasis().order(),
636 lfsv_s.finiteElement().localBasis().order()
640 const auto& cell_inside = ig.inside();
643 auto geo = ig.geometry();
644 auto geo_inside = cell_inside.geometry();
647 auto geo_in_inside = ig.geometryInInside();
651 auto local_inside = ref_el_inside.position(0,0);
652 auto A_s = param.A(cell_inside,local_inside);
656 auto h_F = geo_inside.volume()/geo.volume();
659 auto n_F = ig.centerUnitOuterNormal();
660 Dune::FieldVector<RF,dim> An_F_s;
664 harmonic_average = An_F_s*n_F;
666 harmonic_average = 1.0;
669 auto order_s = lfsu_s.finiteElement().localBasis().order();
670 auto degree = order_s;
673 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
676 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
677 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
678 Dune::FieldVector<RF,dim> gradu_s(0.0);
681 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
684 auto intorder = intorderadd+quadrature_factor*order;
687 auto bctype = param.bctype(ig.intersection(),ip.position());
693 auto iplocal_s = geo_in_inside.global(ip.position());
696 auto n_F_local = ig.unitOuterNormal(ip.position());
699 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
700 auto& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
703 RF factor = ip.weight() * geo.integrationElement(ip.position());
708 auto j = param.j(ig.intersection(),ip.position());
711 for (size_type i=0; i<lfsv_s.size(); i++)
712 r_s.accumulate(lfsv_s,i,j * psi_s[i] * factor);
719 for (size_type i=0; i<lfsu_s.size(); i++)
720 u_s += x_s(lfsu_s,i)*phi_s[i];
723 auto b = param.b(cell_inside,iplocal_s);
724 auto normalflux = b*n_F_local;
728 if (normalflux<-1
e-30)
729 DUNE_THROW(Dune::Exception,
730 "Outflow boundary condition on inflow! [b(" 731 << geo.global(ip.position()) <<
") = " 735 auto term1 = u_s * normalflux *factor;
736 for (size_type i=0; i<lfsv_s.size(); i++)
737 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
740 auto o = param.o(ig.intersection(),ip.position());
743 for (size_type i=0; i<lfsv_s.size(); i++)
744 r_s.accumulate(lfsv_s,i,o * psi_s[i] * factor);
751 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
752 auto& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
755 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
756 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
757 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
761 for (size_type i=0; i<lfsu_s.size(); i++)
762 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
765 auto g = param.g(cell_inside,iplocal_s);
768 RF omegaup_s, omegaup_n;
781 auto term1 = (omegaup_s*u_s + omegaup_n*g) * normalflux *factor;
782 for (size_type i=0; i<lfsv_s.size(); i++)
783 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
786 auto term2 = (An_F_s*gradu_s) * factor;
787 for (size_type i=0; i<lfsv_s.size(); i++)
788 r_s.accumulate(lfsv_s,i,-term2 * psi_s[i]);
791 auto term3 = (u_s-g) * factor;
792 for (size_type i=0; i<lfsv_s.size(); i++)
793 r_s.accumulate(lfsv_s,i,term3 * theta * (An_F_s*tgradpsi_s[i]));
796 auto term4 = penalty_factor * (u_s-g) * factor;
797 for (size_type i=0; i<lfsv_s.size(); i++)
798 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
802 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
804 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
808 using RF =
typename LFSV::Traits::FiniteElementType::
809 Traits::LocalBasisType::Traits::RangeFieldType;
810 using size_type =
typename LFSV::Traits::SizeType;
813 const int dim = IG::dimension;
814 const int order = std::max(
815 lfsu_s.finiteElement().localBasis().order(),
816 lfsv_s.finiteElement().localBasis().order()
820 const auto& cell_inside = ig.inside();
823 auto geo = ig.geometry();
824 auto geo_inside = cell_inside.geometry();
827 auto geo_in_inside = ig.geometryInInside();
831 auto local_inside = ref_el_inside.position(0,0);
832 auto A_s = param.A(cell_inside,local_inside);
836 auto h_F = geo_inside.volume()/geo.volume();
839 auto n_F = ig.centerUnitOuterNormal();
840 Dune::FieldVector<RF,dim> An_F_s;
844 harmonic_average = An_F_s*n_F;
846 harmonic_average = 1.0;
849 auto order_s = lfsu_s.finiteElement().localBasis().order();
850 auto degree = order_s;
853 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
856 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
859 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
862 auto intorder = intorderadd+quadrature_factor*order;
865 auto bctype = param.bctype(ig.intersection(),ip.position());
872 auto iplocal_s = geo_in_inside.global(ip.position());
875 auto n_F_local = ig.unitOuterNormal(ip.position());
878 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
881 auto factor = ip.weight() * geo.integrationElement(ip.position());
884 auto b = param.b(cell_inside,iplocal_s);
885 auto normalflux = b*n_F_local;
889 if (normalflux<-1
e-30)
890 DUNE_THROW(Dune::Exception,
891 "Outflow boundary condition on inflow! [b(" 892 << geo.global(ip.position()) <<
") = " 893 << b <<
")" << n_F_local <<
" " << normalflux);
896 for (size_type j=0; j<lfsu_s.size(); j++)
897 for (size_type i=0; i<lfsu_s.size(); i++)
898 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * normalflux * factor * phi_s[i]);
904 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
907 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
908 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
911 RF omegaup_s, omegaup_n;
924 for (size_type j=0; j<lfsu_s.size(); j++)
925 for (size_type i=0; i<lfsu_s.size(); i++)
926 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux * factor * phi_s[i]);
929 for (size_type j=0; j<lfsu_s.size(); j++)
930 for (size_type i=0; i<lfsu_s.size(); i++)
931 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,-(An_F_s*tgradphi_s[j]) * factor * phi_s[i]);
934 for (size_type j=0; j<lfsu_s.size(); j++)
935 for (size_type i=0; i<lfsu_s.size(); i++)
936 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * (An_F_s*tgradphi_s[i]));
939 for (size_type j=0; j<lfsu_s.size(); j++)
940 for (size_type i=0; i<lfsu_s.size(); i++)
941 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,penalty_factor * phi_s[j] * phi_s[i] * factor);
946 template<
typename EG,
typename LFSV,
typename R>
950 using size_type =
typename LFSV::Traits::SizeType;
953 const auto& cell = eg.entity();
956 auto geo = eg.geometry();
959 auto order = lfsv.finiteElement().localBasis().order();
960 auto intorder = intorderadd + 2 * order;
964 auto& phi = cache[order].evaluateFunction(ip.position(),lfsv.finiteElement().localBasis());
967 auto f = param.f(cell,ip.position());
970 auto factor = ip.weight() * geo.integrationElement(ip.position());
971 for (size_type i=0; i<lfsv.size(); i++)
972 r.accumulate(lfsv,i,-f*phi[i]*factor);
989 int quadrature_factor;
992 using LocalBasisType =
typename FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType;
1004 std::vector<Cache> cache;
1007 void element_size (
const GEO& geo,
typename GEO::ctype& hmin,
typename GEO::ctype hmax)
const 1009 using DF =
typename GEO::ctype;
1012 const int dim = GEO::coorddimension;
1015 Dune::FieldVector<DF,dim> x = geo.corner(0);
1017 hmin = hmax = x.two_norm();
1022 Dune::GeometryType gt = geo.type();
1023 for (
int i=0; i<Dune::ReferenceElements<DF,dim>::general(gt).size(dim-1); i++)
1025 Dune::FieldVector<DF,dim> x = geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,0,dim));
1026 x -= geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,1,dim));
1027 hmin = std::min(hmin,x.two_norm());
1028 hmax = std::max(hmax,x.two_norm());
1036 #endif // DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONDG_HH const IG & ig
Definition: constraints.hh:148
Definition: convectiondiffusiondg.hh:30
static const int dim
Definition: adaptivity.hh:83
const Entity & e
Definition: localfunctionspace.hh:111
void jacobian_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, M &mat_ss, M &mat_sn, M &mat_ns, M &mat_nn) const
Definition: convectiondiffusiondg.hh:451
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericaljacobianapply.hh:285
Definition: convectiondiffusiondg.hh:32
Definition: convectiondiffusiondg.hh:35
Definition: convectiondiffusionparameter.hh:65
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Definition: convectiondiffusiondg.hh:184
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusiondg.hh:623
Definition: convectiondiffusiondg.hh:32
Definition: convectiondiffusiondg.hh:37
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104
Definition: convectiondiffusionparameter.hh:65
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
ConvectionDiffusionDG(T ¶m_, ConvectionDiffusionDGMethod::Type method_=ConvectionDiffusionDGMethod::NIPG, ConvectionDiffusionDGWeights::Type weights_=ConvectionDiffusionDGWeights::weightsOff, Real alpha_=0.0, int intorderadd_=0)
constructor: pass parameter object and define DG-method
Definition: convectiondiffusiondg.hh:86
Type
Definition: convectiondiffusionparameter.hh:65
Definition: convectiondiffusionparameter.hh:65
Definition: convectiondiffusiondg.hh:32
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: convectiondiffusiondg.hh:257
Type
Definition: convectiondiffusiondg.hh:37
void jacobian_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, M &mat_ss) const
Definition: convectiondiffusiondg.hh:803
sparsity pattern generator
Definition: pattern.hh:29
Definition: convectiondiffusiondg.hh:55
Definition: convectiondiffusionparameter.hh:65
sparsity pattern generator
Definition: pattern.hh:13
Type
Definition: convectiondiffusiondg.hh:32
QuadratureRuleWrapper< QuadratureRule< typename Geometry::ctype, Geometry::mydimension > > quadratureRule(const Geometry &geo, std::size_t order, QuadratureType::Enum quadrature_type=QuadratureType::GaussLegendre)
Returns a quadrature rule for the given geometry.
Definition: quadraturerules.hh:117
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:17
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusiondg.hh:191
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:947
Default flags for all local operators.
Definition: flags.hh:18
void jacobian_apply_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, Y &y_s, Y &y_n) const
Definition: convectiondiffusiondg.hh:442
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:104
Definition: convectiondiffusiondg.hh:37
void setTime(Real t)
set time in parameter class
Definition: convectiondiffusiondg.hh:977