2#ifndef DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONDG_HH
3#define DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONDG_HH
8#include<dune/geometry/referenceelements.hh>
10#include<dune/pdelab/common/quadraturerules.hh>
11#include<dune/pdelab/common/function.hh>
12#include<dune/pdelab/localoperator/pattern.hh>
13#include<dune/pdelab/localoperator/flags.hh>
14#include<dune/pdelab/localoperator/idefault.hh>
15#include<dune/pdelab/localoperator/defaultimp.hh>
16#include<dune/pdelab/finiteelement/localbasiscache.hh>
18#include"convectiondiffusionparameter.hh"
29 struct ConvectionDiffusionDGMethod
31 enum Type { NIPG, SIPG, IIPG };
34 struct ConvectionDiffusionDGWeights
36 enum Type { weightsOn, weightsOff };
54 template<
typename T,
typename FiniteElementMap>
61 enum { dim = T::Traits::GridViewType::dimension };
63 using Real =
typename T::Traits::RangeFieldType;
64 using BCType =
typename ConvectionDiffusionBoundaryConditions::Type;
68 enum { doPatternVolume =
true };
69 enum { doPatternSkeleton =
true };
72 enum { doAlphaVolume =
true };
73 enum { doAlphaSkeleton =
true };
74 enum { doAlphaBoundary =
true };
75 enum { doLambdaVolume =
true };
86 ConvectionDiffusionDGMethod::Type method_=ConvectionDiffusionDGMethod::SIPG,
87 ConvectionDiffusionDGWeights::Type weights_=ConvectionDiffusionDGWeights::weightsOn,
95 , intorderadd(intorderadd_)
96 , quadrature_factor(2)
100 if (method==ConvectionDiffusionDGMethod::SIPG) theta = -1.0;
101 if (method==ConvectionDiffusionDGMethod::IIPG) theta = 0.0;
105 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
106 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
109 using RF =
typename LFSU::Traits::FiniteElementType::
110 Traits::LocalBasisType::Traits::RangeFieldType;
111 using size_type =
typename LFSU::Traits::SizeType;
114 const int dim = EG::Entity::dimension;
115 const int order =
std::max(lfsu.finiteElement().localBasis().order(),
116 lfsv.finiteElement().localBasis().order());
119 const auto& cell = eg.entity();
122 auto geo = eg.geometry();
126 auto localcenter = ref_el.position(0,0);
127 auto A = param.A(cell,localcenter);
130 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
131 std::vector<Dune::FieldVector<RF,dim> > gradpsi(lfsv.size());
136 typename EG::Geometry::JacobianInverseTransposed jac;
139 auto intorder = intorderadd + quadrature_factor * order;
140 for (
const auto& ip : quadratureRule(geo,intorder))
143 if (!Impl::permeabilityIsConstantPerCell<T>(param))
145 A = param.A(cell, ip.position());
149 auto& phi = cache[order].evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
150 auto& psi = cache[order].evaluateFunction(ip.position(),lfsv.finiteElement().localBasis());
154 for (size_type i=0; i<lfsu.size(); i++)
155 u += x(lfsu,i)*phi[i];
158 auto& js = cache[order].evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
159 auto& js_v = cache[order].evaluateJacobian(ip.position(),lfsv.finiteElement().localBasis());
162 jac = geo.jacobianInverseTransposed(ip.position());
163 for (size_type i=0; i<lfsu.size(); i++)
164 jac.mv(js[i][0],gradphi[i]);
166 for (size_type i=0; i<lfsv.size(); i++)
167 jac.mv(js_v[i][0],gradpsi[i]);
171 for (size_type i=0; i<lfsu.size(); i++)
172 gradu.axpy(x(lfsu,i),gradphi[i]);
178 auto b = param.b(cell,ip.position());
181 auto c = param.c(cell,ip.position());
184 RF factor = ip.weight() * geo.integrationElement(ip.position());
185 for (size_type i=0; i<lfsv.size(); i++)
186 r.accumulate(lfsv,i,( Agradu*gradpsi[i] - u*(b*gradpsi[i]) + c*u*psi[i] )*factor);
191 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
192 void jacobian_apply_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, Y& y)
const
194 alpha_volume(eg,lfsu,x,lfsv,y);
198 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
199 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
203 using RF =
typename LFSU::Traits::FiniteElementType::
204 Traits::LocalBasisType::Traits::RangeFieldType;
205 using size_type =
typename LFSU::Traits::SizeType;
208 const int dim = EG::Entity::dimension;
209 const int order =
std::max(lfsu.finiteElement().localBasis().order(),
210 lfsv.finiteElement().localBasis().order());
213 const auto& cell = eg.entity();
216 auto geo = eg.geometry();
220 auto localcenter = ref_el.position(0,0);
221 auto A = param.A(cell,localcenter);
224 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
225 std::vector<Dune::FieldVector<RF,dim> > Agradphi(lfsu.size());
228 typename EG::Geometry::JacobianInverseTransposed jac;
231 auto intorder = intorderadd + quadrature_factor * order;
232 for (
const auto& ip : quadratureRule(geo,intorder))
235 if (!Impl::permeabilityIsConstantPerCell<T>(param))
237 A = param.A(cell, ip.position());
241 auto& phi = cache[order].evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
244 auto& js = cache[order].evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
247 jac = geo.jacobianInverseTransposed(ip.position());
248 for (size_type i=0; i<lfsu.size(); i++)
250 jac.mv(js[i][0],gradphi[i]);
251 A.mv(gradphi[i],Agradphi[i]);
255 auto b = param.b(cell,ip.position());
258 auto c = param.c(cell,ip.position());
261 auto factor = ip.weight() * geo.integrationElement(ip.position());
262 for (size_type j=0; j<lfsu.size(); j++)
263 for (size_type i=0; i<lfsu.size(); i++)
264 mat.accumulate(lfsu,i,lfsu,j,( Agradphi[j]*gradphi[i] - phi[j]*(b*gradphi[i]) + c*phi[j]*phi[i] )*factor);
270 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
271 void alpha_skeleton (
const IG& ig,
272 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
273 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
274 R& r_s, R& r_n)
const
277 using RF =
typename LFSV::Traits::FiniteElementType::
278 Traits::LocalBasisType::Traits::RangeFieldType;
279 using size_type =
typename LFSV::Traits::SizeType;
282 const int dim = IG::Entity::dimension;
284 std::max(lfsu_s.finiteElement().localBasis().order(),
285 lfsu_n.finiteElement().localBasis().order()),
286 std::max(lfsv_s.finiteElement().localBasis().order(),
287 lfsv_n.finiteElement().localBasis().order())
291 const auto& cell_inside = ig.inside();
292 const auto& cell_outside = ig.outside();
295 auto geo = ig.geometry();
296 auto geo_inside = cell_inside.geometry();
297 auto geo_outside = cell_outside.geometry();
300 auto geo_in_inside = ig.geometryInInside();
301 auto geo_in_outside = ig.geometryInOutside();
306 auto local_inside = ref_el_inside.position(0,0);
307 auto local_outside = ref_el_outside.position(0,0);
308 auto A_s = param.A(cell_inside,local_inside);
309 auto A_n = param.A(cell_outside,local_outside);
313 auto h_F =
std::min(geo_inside.volume(),geo_outside.volume())/geo.volume();
316 auto n_F = ig.centerUnitOuterNormal();
325 RF harmonic_average(0.0);
326 if (weights==ConvectionDiffusionDGWeights::weightsOn)
328 RF delta_s = (An_F_s*n_F);
329 RF delta_n = (An_F_n*n_F);
330 omega_s = delta_n/(delta_s+delta_n+1e-20);
331 omega_n = delta_s/(delta_s+delta_n+1e-20);
332 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1e-20);
336 omega_s = omega_n = 0.5;
337 harmonic_average = 1.0;
341 auto order_s = lfsu_s.finiteElement().localBasis().order();
342 auto order_n = lfsu_n.finiteElement().localBasis().order();
346 auto penalty_factor = (alpha/h_F) * harmonic_average *
degree*(
degree+dim-1);
349 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
350 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
351 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
352 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_n(lfsv_n.size());
357 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
360 auto intorder = intorderadd+quadrature_factor*order;
361 for (
const auto& ip : quadratureRule(geo,intorder))
364 auto n_F_local = ig.unitOuterNormal(ip.position());
367 if (!Impl::permeabilityIsConstantPerCell<T>(param))
369 A_s = param.A(cell_inside,geo_in_inside.global(ip.position()));
370 A_n = param.A(cell_outside,geo_in_outside.global(ip.position()));
371 A_s.mv(n_F_local,An_F_s);
372 A_n.mv(n_F_local,An_F_n);
373 if (weights==ConvectionDiffusionDGWeights::weightsOn)
375 RF delta_s = (An_F_s*n_F_local);
376 RF delta_n = (An_F_n*n_F_local);
377 omega_s = delta_n/(delta_s+delta_n+1e-20);
378 omega_n = delta_s/(delta_s+delta_n+1e-20);
379 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1e-20);
380 penalty_factor = (alpha/h_F) * harmonic_average *
degree*(
degree+dim-1);
386 auto iplocal_s = geo_in_inside.global(ip.position());
387 auto iplocal_n = geo_in_outside.global(ip.position());
390 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
391 auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
392 auto& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
393 auto& psi_n = cache[order_n].evaluateFunction(iplocal_n,lfsv_n.finiteElement().localBasis());
397 for (size_type i=0; i<lfsu_s.size(); i++)
398 u_s += x_s(lfsu_s,i)*phi_s[i];
400 for (size_type i=0; i<lfsu_n.size(); i++)
401 u_n += x_n(lfsu_n,i)*phi_n[i];
404 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
405 auto& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
406 auto& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
407 auto& gradpsi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsv_n.finiteElement().localBasis());
410 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
411 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
412 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
413 jac = geo_outside.jacobianInverseTransposed(iplocal_n);
414 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
415 for (size_type i=0; i<lfsv_n.size(); i++) jac.mv(gradpsi_n[i][0],tgradpsi_n[i]);
419 for (size_type i=0; i<lfsu_s.size(); i++)
420 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
422 for (size_type i=0; i<lfsu_n.size(); i++)
423 gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
426 auto b = param.b(cell_inside,iplocal_s);
427 auto normalflux = b*n_F_local;
428 RF omegaup_s, omegaup_n;
441 auto factor = ip.weight()*geo.integrationElement(ip.position());
444 auto term1 = (omegaup_s*u_s + omegaup_n*u_n) * normalflux *factor;
445 for (size_type i=0; i<lfsv_s.size(); i++)
446 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
447 for (size_type i=0; i<lfsv_n.size(); i++)
448 r_n.accumulate(lfsv_n,i,-term1 * psi_n[i]);
451 auto term2 = -(omega_s*(An_F_s*gradu_s) + omega_n*(An_F_n*gradu_n)) * factor;
452 for (size_type i=0; i<lfsv_s.size(); i++)
453 r_s.accumulate(lfsv_s,i,term2 * psi_s[i]);
454 for (size_type i=0; i<lfsv_n.size(); i++)
455 r_n.accumulate(lfsv_n,i,-term2 * psi_n[i]);
458 auto term3 = (u_s-u_n) * factor;
459 for (size_type i=0; i<lfsv_s.size(); i++)
460 r_s.accumulate(lfsv_s,i,term3 * theta * omega_s * (An_F_s*tgradpsi_s[i]));
461 for (size_type i=0; i<lfsv_n.size(); i++)
462 r_n.accumulate(lfsv_n,i,term3 * theta * omega_n * (An_F_n*tgradpsi_n[i]));
465 auto term4 = penalty_factor * (u_s-u_n) * factor;
466 for (size_type i=0; i<lfsv_s.size(); i++)
467 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
468 for (size_type i=0; i<lfsv_n.size(); i++)
469 r_n.accumulate(lfsv_n,i,-term4 * psi_n[i]);
474 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
475 void jacobian_apply_skeleton (
const IG& ig,
476 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
477 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
478 Y& y_s, Y& y_n)
const
480 alpha_skeleton(ig,lfsu_s,x_s,lfsv_s,lfsu_n,x_n,lfsv_n,y_s,y_n);
483 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
484 void jacobian_skeleton (
const IG& ig,
485 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
486 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
487 M& mat_ss, M& mat_sn,
488 M& mat_ns, M& mat_nn)
const
491 using RF =
typename LFSV::Traits::FiniteElementType::
492 Traits::LocalBasisType::Traits::RangeFieldType;
493 using size_type =
typename LFSV::Traits::SizeType;
496 const int dim = IG::Entity::dimension;
498 std::max(lfsu_s.finiteElement().localBasis().order(),
499 lfsu_n.finiteElement().localBasis().order()),
500 std::max(lfsv_s.finiteElement().localBasis().order(),
501 lfsv_n.finiteElement().localBasis().order())
505 const auto& cell_inside = ig.inside();
506 const auto& cell_outside = ig.outside();
509 auto geo = ig.geometry();
510 auto geo_inside = cell_inside.geometry();
511 auto geo_outside = cell_outside.geometry();
514 auto geo_in_inside = ig.geometryInInside();
515 auto geo_in_outside = ig.geometryInOutside();
520 auto local_inside = ref_el_inside.position(0,0);
521 auto local_outside = ref_el_outside.position(0,0);
522 auto A_s = param.A(cell_inside,local_inside);
523 auto A_n = param.A(cell_outside,local_outside);
527 auto h_F =
std::min(geo_inside.volume(),geo_outside.volume())/geo.volume();
530 auto n_F = ig.centerUnitOuterNormal();
539 RF harmonic_average(0.0);
540 if (weights==ConvectionDiffusionDGWeights::weightsOn)
542 RF delta_s = (An_F_s*n_F);
543 RF delta_n = (An_F_n*n_F);
544 omega_s = delta_n/(delta_s+delta_n+1e-20);
545 omega_n = delta_s/(delta_s+delta_n+1e-20);
546 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1e-20);
550 omega_s = omega_n = 0.5;
551 harmonic_average = 1.0;
555 auto order_s = lfsu_s.finiteElement().localBasis().order();
556 auto order_n = lfsu_n.finiteElement().localBasis().order();
560 auto penalty_factor = (alpha/h_F) * harmonic_average *
degree*(
degree+dim-1);
563 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
564 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
567 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
570 const int intorder = intorderadd+quadrature_factor*order;
571 for (
const auto& ip : quadratureRule(geo,intorder))
574 auto n_F_local = ig.unitOuterNormal(ip.position());
577 if (!Impl::permeabilityIsConstantPerCell<T>(param))
579 A_s = param.A(cell_inside,geo_in_inside.global(ip.position()));
580 A_n = param.A(cell_outside,geo_in_outside.global(ip.position()));
581 A_s.mv(n_F_local,An_F_s);
582 A_n.mv(n_F_local,An_F_n);
583 if (weights==ConvectionDiffusionDGWeights::weightsOn)
585 RF delta_s = (An_F_s*n_F_local);
586 RF delta_n = (An_F_n*n_F_local);
587 omega_s = delta_n/(delta_s+delta_n+1e-20);
588 omega_n = delta_s/(delta_s+delta_n+1e-20);
589 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1e-20);
590 penalty_factor = (alpha/h_F) * harmonic_average *
degree*(
degree+dim-1);
595 auto iplocal_s = geo_in_inside.global(ip.position());
596 auto iplocal_n = geo_in_outside.global(ip.position());
599 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
600 auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
603 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
604 auto& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
607 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
608 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
609 jac = geo_outside.jacobianInverseTransposed(iplocal_n);
610 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
613 auto b = param.b(cell_inside,iplocal_s);
614 auto normalflux = b*n_F_local;
615 RF omegaup_s, omegaup_n;
628 auto factor = ip.weight() * geo.integrationElement(ip.position());
629 auto ipfactor = penalty_factor * factor;
632 for (size_type j=0; j<lfsu_s.size(); j++) {
633 auto temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
634 for (size_type i=0; i<lfsu_s.size(); i++) {
635 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux *factor * phi_s[i]);
636 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,temp1 * phi_s[i]);
637 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
638 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * ipfactor * phi_s[i]);
641 for (size_type j=0; j<lfsu_n.size(); j++) {
642 auto temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
643 for (size_type i=0; i<lfsu_s.size(); i++) {
644 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,omegaup_n * phi_n[j] * normalflux *factor * phi_s[i]);
645 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,temp1 * phi_s[i]);
646 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
647 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * ipfactor * phi_s[i]);
650 for (size_type j=0; j<lfsu_s.size(); j++) {
651 auto temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
652 for (size_type i=0; i<lfsu_n.size(); i++) {
653 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-omegaup_s * phi_s[j] * normalflux *factor * phi_n[i]);
654 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-temp1 * phi_n[i]);
655 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,phi_s[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
656 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-phi_s[j] * ipfactor * phi_n[i]);
659 for (size_type j=0; j<lfsu_n.size(); j++) {
660 auto temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
661 for (size_type i=0; i<lfsu_n.size(); i++) {
662 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-omegaup_n * phi_n[j] * normalflux *factor * phi_n[i]);
663 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-temp1 * phi_n[i]);
664 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
665 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,phi_n[j] * ipfactor * phi_n[i]);
683 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
684 void residual_boundary_integral (
const IG& ig,
685 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
686 R& r_s,
bool jacobian_apply=
false)
const
689 using RF =
typename LFSV::Traits::FiniteElementType::
690 Traits::LocalBasisType::Traits::RangeFieldType;
691 using size_type =
typename LFSV::Traits::SizeType;
694 const int dim = IG::Entity::dimension;
696 lfsu_s.finiteElement().localBasis().order(),
697 lfsv_s.finiteElement().localBasis().order()
701 const auto& cell_inside = ig.inside();
704 auto geo = ig.geometry();
705 auto geo_inside = cell_inside.geometry();
708 auto geo_in_inside = ig.geometryInInside();
712 auto local_inside = ref_el_inside.position(0,0);
713 auto A_s = param.A(cell_inside,local_inside);
717 auto h_F = geo_inside.volume()/geo.volume();
720 auto n_F = ig.centerUnitOuterNormal();
724 if (weights==ConvectionDiffusionDGWeights::weightsOn)
725 harmonic_average = An_F_s*n_F;
727 harmonic_average = 1.0;
730 auto order_s = lfsu_s.finiteElement().localBasis().order();
734 auto penalty_factor = (alpha/h_F) * harmonic_average *
degree*(
degree+dim-1);
737 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
738 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
742 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
745 auto intorder = intorderadd+quadrature_factor*order;
746 for (
const auto& ip : quadratureRule(geo,intorder))
749 auto n_F_local = ig.unitOuterNormal(ip.position());
752 if (!Impl::permeabilityIsConstantPerCell<T>(param))
754 A_s = param.A(cell_inside,geo_in_inside.global(ip.position()));
755 A_s.mv(n_F_local,An_F_s);
756 if (weights==ConvectionDiffusionDGWeights::weightsOn)
758 harmonic_average = An_F_s*n_F_local;
759 penalty_factor = (alpha/h_F) * harmonic_average *
degree*(
degree+dim-1);
763 auto bctype = param.bctype(ig.intersection(),ip.position());
765 if (bctype == ConvectionDiffusionBoundaryConditions::None)
769 auto iplocal_s = geo_in_inside.global(ip.position());
772 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
773 auto& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
776 RF factor = ip.weight() * geo.integrationElement(ip.position());
778 if (bctype == ConvectionDiffusionBoundaryConditions::Neumann)
780 if (not jacobian_apply){
782 auto j = param.j(ig.intersection(),ip.position());
785 for (size_type i=0; i<lfsv_s.size(); i++)
786 r_s.accumulate(lfsv_s,i,j * psi_s[i] * factor);
793 for (size_type i=0; i<lfsu_s.size(); i++)
794 u_s += x_s(lfsu_s,i)*phi_s[i];
797 auto b = param.b(cell_inside,iplocal_s);
798 auto normalflux = b*n_F_local;
800 if (bctype == ConvectionDiffusionBoundaryConditions::Outflow)
802 if (normalflux<-1e-30)
804 "Outflow boundary condition on inflow! [b("
805 << geo.global(ip.position()) <<
") = "
809 auto term1 = u_s * normalflux *factor;
810 for (size_type i=0; i<lfsv_s.size(); i++)
811 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
813 if (not jacobian_apply){
815 auto o = param.o(ig.intersection(),ip.position());
818 for (size_type i=0; i<lfsv_s.size(); i++)
819 r_s.accumulate(lfsv_s,i,o * psi_s[i] * factor);
825 assert (bctype == ConvectionDiffusionBoundaryConditions::Dirichlet);
826 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
827 auto& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
830 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
831 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
832 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
836 for (size_type i=0; i<lfsu_s.size(); i++)
837 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
840 auto g = param.g(cell_inside,iplocal_s);
847 RF omegaup_s, omegaup_n;
860 auto term1 = (omegaup_s*u_s + omegaup_n*g) * normalflux *factor;
861 for (size_type i=0; i<lfsv_s.size(); i++)
862 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
865 auto term2 = (An_F_s*gradu_s) * factor;
866 for (size_type i=0; i<lfsv_s.size(); i++)
867 r_s.accumulate(lfsv_s,i,-term2 * psi_s[i]);
870 auto term3 = (u_s-g) * factor;
871 for (size_type i=0; i<lfsv_s.size(); i++)
872 r_s.accumulate(lfsv_s,i,term3 * theta * (An_F_s*tgradpsi_s[i]));
875 auto term4 = penalty_factor * (u_s-g) * factor;
876 for (size_type i=0; i<lfsv_s.size(); i++)
877 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
883 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
884 void alpha_boundary (
const IG& ig,
885 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
888 residual_boundary_integral(ig, lfsu_s, x_s, lfsv_s, r_s);
892 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
893 void jacobian_apply_boundary (
const IG& ig,
894 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
898 residual_boundary_integral(ig, lfsu_s, x_s, lfsv_s, y_s,
true);
901 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
902 void jacobian_boundary (
const IG& ig,
903 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
907 using RF =
typename LFSV::Traits::FiniteElementType::
908 Traits::LocalBasisType::Traits::RangeFieldType;
909 using size_type =
typename LFSV::Traits::SizeType;
912 const int dim = IG::Entity::dimension;
914 lfsu_s.finiteElement().localBasis().order(),
915 lfsv_s.finiteElement().localBasis().order()
919 const auto& cell_inside = ig.inside();
922 auto geo = ig.geometry();
923 auto geo_inside = cell_inside.geometry();
926 auto geo_in_inside = ig.geometryInInside();
930 auto local_inside = ref_el_inside.position(0,0);
931 auto A_s = param.A(cell_inside,local_inside);
935 auto h_F = geo_inside.volume()/geo.volume();
938 auto n_F = ig.centerUnitOuterNormal();
942 if (weights==ConvectionDiffusionDGWeights::weightsOn)
943 harmonic_average = An_F_s*n_F;
945 harmonic_average = 1.0;
948 auto order_s = lfsu_s.finiteElement().localBasis().order();
952 auto penalty_factor = (alpha/h_F) * harmonic_average *
degree*(
degree+dim-1);
955 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
958 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
961 auto intorder = intorderadd+quadrature_factor*order;
962 for (
const auto& ip : quadratureRule(geo,intorder))
965 auto n_F_local = ig.unitOuterNormal(ip.position());
968 if (!Impl::permeabilityIsConstantPerCell<T>(param))
970 A_s = param.A(cell_inside,geo_in_inside.global(ip.position()));
971 A_s.mv(n_F_local,An_F_s);
972 if (weights==ConvectionDiffusionDGWeights::weightsOn)
974 harmonic_average = An_F_s*n_F_local;
975 penalty_factor = (alpha/h_F) * harmonic_average *
degree*(
degree+dim-1);
979 auto bctype = param.bctype(ig.intersection(),ip.position());
981 if (bctype == ConvectionDiffusionBoundaryConditions::None ||
982 bctype == ConvectionDiffusionBoundaryConditions::Neumann)
986 auto iplocal_s = geo_in_inside.global(ip.position());
989 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
992 auto factor = ip.weight() * geo.integrationElement(ip.position());
995 auto b = param.b(cell_inside,iplocal_s);
996 auto normalflux = b*n_F_local;
998 if (bctype == ConvectionDiffusionBoundaryConditions::Outflow)
1000 if (normalflux<-1e-30)
1002 "Outflow boundary condition on inflow! [b("
1003 << geo.global(ip.position()) <<
") = "
1004 << b <<
")" << n_F_local <<
" " << normalflux);
1007 for (size_type j=0; j<lfsu_s.size(); j++)
1008 for (size_type i=0; i<lfsu_s.size(); i++)
1009 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * normalflux * factor * phi_s[i]);
1015 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
1018 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
1019 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
1022 RF omegaup_s = normalflux>=0.0 ? 1.0 : 0.0;
1025 for (size_type j=0; j<lfsu_s.size(); j++)
1026 for (size_type i=0; i<lfsu_s.size(); i++)
1027 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux * factor * phi_s[i]);
1030 for (size_type j=0; j<lfsu_s.size(); j++)
1031 for (size_type i=0; i<lfsu_s.size(); i++)
1032 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,-(An_F_s*tgradphi_s[j]) * factor * phi_s[i]);
1035 for (size_type j=0; j<lfsu_s.size(); j++)
1036 for (size_type i=0; i<lfsu_s.size(); i++)
1037 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * (An_F_s*tgradphi_s[i]));
1040 for (size_type j=0; j<lfsu_s.size(); j++)
1041 for (size_type i=0; i<lfsu_s.size(); i++)
1042 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,penalty_factor * phi_s[j] * phi_s[i] * factor);
1047 template<
typename EG,
typename LFSV,
typename R>
1048 void lambda_volume (
const EG& eg,
const LFSV& lfsv, R& r)
const
1051 using size_type =
typename LFSV::Traits::SizeType;
1054 const auto& cell = eg.entity();
1057 auto geo = eg.geometry();
1060 auto order = lfsv.finiteElement().localBasis().order();
1061 auto intorder = intorderadd + 2 * order;
1062 for (
const auto& ip : quadratureRule(geo,intorder))
1065 auto& phi = cache[order].evaluateFunction(ip.position(),lfsv.finiteElement().localBasis());
1068 auto f = param.f(cell,ip.position());
1071 auto factor = ip.weight() * geo.integrationElement(ip.position());
1072 for (size_type i=0; i<lfsv.size(); i++)
1073 r.accumulate(lfsv,i,-f*phi[i]*factor);
1086 ConvectionDiffusionDGMethod::Type method;
1087 ConvectionDiffusionDGWeights::Type weights;
1090 int quadrature_factor;
1093 using LocalBasisType =
typename FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType;
1105 std::vector<Cache> cache;
1108 void element_size (
const GEO& geo,
typename GEO::ctype& hmin,
typename GEO::ctype hmax)
const
1110 using DF =
typename GEO::ctype;
1113 const int dim = GEO::coorddimension;
1124 for (
int i=0; i<Dune::ReferenceElements<DF,dim>::general(
gt).size(dim-1); i++)
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:641
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
Definition: convectiondiffusiondg.hh:60
void setTime(Real t)
set time in parameter class
Definition: convectiondiffusiondg.hh:1078
ConvectionDiffusionDG(T ¶m_, ConvectionDiffusionDGMethod::Type method_=ConvectionDiffusionDGMethod::SIPG, ConvectionDiffusionDGWeights::Type weights_=ConvectionDiffusionDGWeights::weightsOn, Real alpha_=1.0, int intorderadd_=0)
constructor: pass parameter object and define DG-method
Definition: convectiondiffusiondg.hh:85
sparsity pattern generator
Definition: pattern.hh:30
sparsity pattern generator
Definition: pattern.hh:14
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:18
Default flags for all local operators.
Definition: flags.hh:19
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:158
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
std::size_t degree(const Node &node)
Returns the degree of node as run time information.
Definition: nodeinterface.hh:79
Dune namespace.
Definition: alignedallocator.hh:13
Class providing access to the singletons of the reference elements.
Definition: referenceelements.hh:128