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> 53 template<
typename T,
typename FiniteElementMap>
61 enum {
dim = T::Traits::GridViewType::dimension };
63 using Real =
typename T::Traits::RangeFieldType;
68 enum { doPatternVolume =
true };
69 enum { doPatternSkeleton =
true };
72 enum { doAlphaVolume =
true };
73 enum { doAlphaSkeleton =
true };
74 enum { doAlphaBoundary =
true };
75 enum { doLambdaVolume =
true };
92 param(param_), method(method_), weights(weights_),
93 alpha(alpha_), intorderadd(intorderadd_), quadrature_factor(2),
102 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
103 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 106 using RF =
typename LFSU::Traits::FiniteElementType::
107 Traits::LocalBasisType::Traits::RangeFieldType;
108 using size_type =
typename LFSU::Traits::SizeType;
111 const int dim = EG::Entity::dimension;
112 const int order = std::max(lfsu.finiteElement().localBasis().order(),
113 lfsv.finiteElement().localBasis().order());
116 const auto& cell = eg.entity();
119 auto geo = eg.geometry();
122 auto ref_el = referenceElement(geo);
123 auto localcenter = ref_el.position(0,0);
124 auto A = param.A(cell,localcenter);
127 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
128 std::vector<Dune::FieldVector<RF,dim> > gradpsi(lfsv.size());
129 Dune::FieldVector<RF,dim> gradu(0.0);
130 Dune::FieldVector<RF,dim> Agradu(0.0);
133 typename EG::Geometry::JacobianInverseTransposed jac;
136 auto intorder = intorderadd + quadrature_factor * order;
140 auto& phi = cache[order].evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
141 auto& psi = cache[order].evaluateFunction(ip.position(),lfsv.finiteElement().localBasis());
145 for (size_type i=0; i<lfsu.size(); i++)
146 u += x(lfsu,i)*phi[i];
149 auto& js = cache[order].evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
150 auto& js_v = cache[order].evaluateJacobian(ip.position(),lfsv.finiteElement().localBasis());
153 jac = geo.jacobianInverseTransposed(ip.position());
154 for (size_type i=0; i<lfsu.size(); i++)
155 jac.mv(js[i][0],gradphi[i]);
157 for (size_type i=0; i<lfsv.size(); i++)
158 jac.mv(js_v[i][0],gradpsi[i]);
162 for (size_type i=0; i<lfsu.size(); i++)
163 gradu.axpy(x(lfsu,i),gradphi[i]);
169 auto b = param.b(cell,ip.position());
172 auto c = param.c(cell,ip.position());
175 RF factor = ip.weight() * geo.integrationElement(ip.position());
176 for (size_type i=0; i<lfsv.size(); i++)
177 r.accumulate(lfsv,i,( Agradu*gradpsi[i] - u*(b*gradpsi[i]) + c*u*psi[i] )*factor);
182 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
185 alpha_volume(eg,lfsu,x,lfsv,y);
189 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
194 using RF =
typename LFSU::Traits::FiniteElementType::
195 Traits::LocalBasisType::Traits::RangeFieldType;
196 using size_type =
typename LFSU::Traits::SizeType;
199 const int dim = EG::Entity::dimension;
200 const int order = std::max(lfsu.finiteElement().localBasis().order(),
201 lfsv.finiteElement().localBasis().order());
204 const auto& cell = eg.entity();
207 auto geo = eg.geometry();
210 auto ref_el = referenceElement(geo);
211 auto localcenter = ref_el.position(0,0);
212 auto A = param.A(cell,localcenter);
215 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
216 std::vector<Dune::FieldVector<RF,dim> > Agradphi(lfsu.size());
219 typename EG::Geometry::JacobianInverseTransposed jac;
222 auto intorder = intorderadd + quadrature_factor * order;
226 auto& phi = cache[order].evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
229 auto& js = cache[order].evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
232 jac = geo.jacobianInverseTransposed(ip.position());
233 for (size_type i=0; i<lfsu.size(); i++)
235 jac.mv(js[i][0],gradphi[i]);
236 A.mv(gradphi[i],Agradphi[i]);
240 auto b = param.b(cell,ip.position());
243 auto c = param.c(cell,ip.position());
246 auto factor = ip.weight() * geo.integrationElement(ip.position());
247 for (size_type j=0; j<lfsu.size(); j++)
248 for (size_type i=0; i<lfsu.size(); i++)
249 mat.accumulate(lfsu,i,lfsu,j,( Agradphi[j]*gradphi[i] - phi[j]*(b*gradphi[i]) + c*phi[j]*phi[i] )*factor);
255 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
257 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
258 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
259 R& r_s, R& r_n)
const 262 using RF =
typename LFSV::Traits::FiniteElementType::
263 Traits::LocalBasisType::Traits::RangeFieldType;
264 using size_type =
typename LFSV::Traits::SizeType;
267 const int dim = IG::Entity::dimension;
268 const int order = std::max(
269 std::max(lfsu_s.finiteElement().localBasis().order(),
270 lfsu_n.finiteElement().localBasis().order()),
271 std::max(lfsv_s.finiteElement().localBasis().order(),
272 lfsv_n.finiteElement().localBasis().order())
276 const auto& cell_inside = ig.inside();
277 const auto& cell_outside = ig.outside();
280 auto geo = ig.geometry();
281 auto geo_inside = cell_inside.geometry();
282 auto geo_outside = cell_outside.geometry();
285 auto geo_in_inside = ig.geometryInInside();
286 auto geo_in_outside = ig.geometryInOutside();
289 auto ref_el_inside = referenceElement(geo_inside);
290 auto ref_el_outside = referenceElement(geo_outside);
291 auto local_inside = ref_el_inside.position(0,0);
292 auto local_outside = ref_el_outside.position(0,0);
293 auto A_s = param.A(cell_inside,local_inside);
294 auto A_n = param.A(cell_outside,local_outside);
298 auto h_F = std::min(geo_inside.volume(),geo_outside.volume())/geo.volume();
301 auto n_F = ig.centerUnitOuterNormal();
302 Dune::FieldVector<RF,dim> An_F_s;
304 Dune::FieldVector<RF,dim> An_F_n;
310 RF harmonic_average(0.0);
313 RF delta_s = (An_F_s*n_F);
314 RF delta_n = (An_F_n*n_F);
315 omega_s = delta_n/(delta_s+delta_n+1
e-20);
316 omega_n = delta_s/(delta_s+delta_n+1
e-20);
317 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-20);
321 omega_s = omega_n = 0.5;
322 harmonic_average = 1.0;
326 auto order_s = lfsu_s.finiteElement().localBasis().order();
327 auto order_n = lfsu_n.finiteElement().localBasis().order();
328 auto degree = std::max( order_s, order_n );
331 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
334 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
335 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
336 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
337 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_n(lfsv_n.size());
338 Dune::FieldVector<RF,dim> gradu_s(0.0);
339 Dune::FieldVector<RF,dim> gradu_n(0.0);
342 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
345 auto intorder = intorderadd+quadrature_factor*order;
349 auto n_F_local = ig.unitOuterNormal(ip.position());
352 auto iplocal_s = geo_in_inside.global(ip.position());
353 auto iplocal_n = geo_in_outside.global(ip.position());
356 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
357 auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
358 auto& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
359 auto& psi_n = cache[order_n].evaluateFunction(iplocal_n,lfsv_n.finiteElement().localBasis());
363 for (size_type i=0; i<lfsu_s.size(); i++)
364 u_s += x_s(lfsu_s,i)*phi_s[i];
366 for (size_type i=0; i<lfsu_n.size(); i++)
367 u_n += x_n(lfsu_n,i)*phi_n[i];
370 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
371 auto& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
372 auto& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
373 auto& gradpsi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsv_n.finiteElement().localBasis());
376 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
377 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
378 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
379 jac = geo_outside.jacobianInverseTransposed(iplocal_n);
380 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
381 for (size_type i=0; i<lfsv_n.size(); i++) jac.mv(gradpsi_n[i][0],tgradpsi_n[i]);
385 for (size_type i=0; i<lfsu_s.size(); i++)
386 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
388 for (size_type i=0; i<lfsu_n.size(); i++)
389 gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
392 auto b = param.b(cell_inside,iplocal_s);
393 auto normalflux = b*n_F_local;
394 RF omegaup_s, omegaup_n;
407 auto factor = ip.weight()*geo.integrationElement(ip.position());
410 auto term1 = (omegaup_s*u_s + omegaup_n*u_n) * normalflux *factor;
411 for (size_type i=0; i<lfsv_s.size(); i++)
412 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
413 for (size_type i=0; i<lfsv_n.size(); i++)
414 r_n.accumulate(lfsv_n,i,-term1 * psi_n[i]);
417 auto term2 = -(omega_s*(An_F_s*gradu_s) + omega_n*(An_F_n*gradu_n)) * factor;
418 for (size_type i=0; i<lfsv_s.size(); i++)
419 r_s.accumulate(lfsv_s,i,term2 * psi_s[i]);
420 for (size_type i=0; i<lfsv_n.size(); i++)
421 r_n.accumulate(lfsv_n,i,-term2 * psi_n[i]);
424 auto term3 = (u_s-u_n) * factor;
425 for (size_type i=0; i<lfsv_s.size(); i++)
426 r_s.accumulate(lfsv_s,i,term3 * theta * omega_s * (An_F_s*tgradpsi_s[i]));
427 for (size_type i=0; i<lfsv_n.size(); i++)
428 r_n.accumulate(lfsv_n,i,term3 * theta * omega_n * (An_F_n*tgradpsi_n[i]));
431 auto term4 = penalty_factor * (u_s-u_n) * factor;
432 for (size_type i=0; i<lfsv_s.size(); i++)
433 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
434 for (size_type i=0; i<lfsv_n.size(); i++)
435 r_n.accumulate(lfsv_n,i,-term4 * psi_n[i]);
440 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
442 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
443 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
444 Y& y_s, Y& y_n)
const 446 alpha_skeleton(ig,lfsu_s,x_s,lfsv_s,lfsu_n,x_n,lfsv_n,y_s,y_n);
449 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
451 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
452 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
453 M& mat_ss, M& mat_sn,
454 M& mat_ns, M& mat_nn)
const 457 using RF =
typename LFSV::Traits::FiniteElementType::
458 Traits::LocalBasisType::Traits::RangeFieldType;
459 using size_type =
typename LFSV::Traits::SizeType;
462 const int dim = IG::Entity::dimension;
463 const int order = std::max(
464 std::max(lfsu_s.finiteElement().localBasis().order(),
465 lfsu_n.finiteElement().localBasis().order()),
466 std::max(lfsv_s.finiteElement().localBasis().order(),
467 lfsv_n.finiteElement().localBasis().order())
471 const auto& cell_inside = ig.inside();
472 const auto& cell_outside = ig.outside();
475 auto geo = ig.geometry();
476 auto geo_inside = cell_inside.geometry();
477 auto geo_outside = cell_outside.geometry();
480 auto geo_in_inside = ig.geometryInInside();
481 auto geo_in_outside = ig.geometryInOutside();
484 auto ref_el_inside = referenceElement(geo_inside);
485 auto ref_el_outside = referenceElement(geo_outside);
486 auto local_inside = ref_el_inside.position(0,0);
487 auto local_outside = ref_el_outside.position(0,0);
488 auto A_s = param.A(cell_inside,local_inside);
489 auto A_n = param.A(cell_outside,local_outside);
493 auto h_F = std::min(geo_inside.volume(),geo_outside.volume())/geo.volume();
496 auto n_F = ig.centerUnitOuterNormal();
497 Dune::FieldVector<RF,dim> An_F_s;
499 Dune::FieldVector<RF,dim> An_F_n;
505 RF harmonic_average(0.0);
508 RF delta_s = (An_F_s*n_F);
509 RF delta_n = (An_F_n*n_F);
510 omega_s = delta_n/(delta_s+delta_n+1
e-20);
511 omega_n = delta_s/(delta_s+delta_n+1
e-20);
512 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-20);
516 omega_s = omega_n = 0.5;
517 harmonic_average = 1.0;
521 auto order_s = lfsu_s.finiteElement().localBasis().order();
522 auto order_n = lfsu_n.finiteElement().localBasis().order();
523 auto degree = std::max( order_s, order_n );
526 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
529 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
530 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
533 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
536 const int intorder = intorderadd+quadrature_factor*order;
540 auto n_F_local = ig.unitOuterNormal(ip.position());
543 auto iplocal_s = geo_in_inside.global(ip.position());
544 auto iplocal_n = geo_in_outside.global(ip.position());
547 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
548 auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
551 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
552 auto& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
555 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
556 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
557 jac = geo_outside.jacobianInverseTransposed(iplocal_n);
558 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
561 auto b = param.b(cell_inside,iplocal_s);
562 auto normalflux = b*n_F_local;
563 RF omegaup_s, omegaup_n;
576 auto factor = ip.weight() * geo.integrationElement(ip.position());
577 auto ipfactor = penalty_factor * factor;
580 for (size_type j=0; j<lfsu_s.size(); j++) {
581 auto temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
582 for (size_type i=0; i<lfsu_s.size(); i++) {
583 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux *factor * phi_s[i]);
584 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,temp1 * phi_s[i]);
585 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
586 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * ipfactor * phi_s[i]);
589 for (size_type j=0; j<lfsu_n.size(); j++) {
590 auto temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
591 for (size_type i=0; i<lfsu_s.size(); i++) {
592 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,omegaup_n * phi_n[j] * normalflux *factor * phi_s[i]);
593 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,temp1 * phi_s[i]);
594 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
595 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * ipfactor * phi_s[i]);
598 for (size_type j=0; j<lfsu_s.size(); j++) {
599 auto temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
600 for (size_type i=0; i<lfsu_n.size(); i++) {
601 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-omegaup_s * phi_s[j] * normalflux *factor * phi_n[i]);
602 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-temp1 * phi_n[i]);
603 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,phi_s[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
604 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-phi_s[j] * ipfactor * phi_n[i]);
607 for (size_type j=0; j<lfsu_n.size(); j++) {
608 auto temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
609 for (size_type i=0; i<lfsu_n.size(); i++) {
610 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-omegaup_n * phi_n[j] * normalflux *factor * phi_n[i]);
611 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-temp1 * phi_n[i]);
612 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
613 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,phi_n[j] * ipfactor * phi_n[i]);
621 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
623 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
627 using RF =
typename LFSV::Traits::FiniteElementType::
628 Traits::LocalBasisType::Traits::RangeFieldType;
629 using size_type =
typename LFSV::Traits::SizeType;
632 const int dim = IG::Entity::dimension;
633 const int order = std::max(
634 lfsu_s.finiteElement().localBasis().order(),
635 lfsv_s.finiteElement().localBasis().order()
639 const auto& cell_inside = ig.inside();
642 auto geo = ig.geometry();
643 auto geo_inside = cell_inside.geometry();
646 auto geo_in_inside = ig.geometryInInside();
649 auto ref_el_inside = referenceElement(geo_inside);
650 auto local_inside = ref_el_inside.position(0,0);
651 auto A_s = param.A(cell_inside,local_inside);
655 auto h_F = geo_inside.volume()/geo.volume();
658 auto n_F = ig.centerUnitOuterNormal();
659 Dune::FieldVector<RF,dim> An_F_s;
663 harmonic_average = An_F_s*n_F;
665 harmonic_average = 1.0;
668 auto order_s = lfsu_s.finiteElement().localBasis().order();
669 auto degree = order_s;
672 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
675 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
676 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
677 Dune::FieldVector<RF,dim> gradu_s(0.0);
680 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
683 auto intorder = intorderadd+quadrature_factor*order;
686 auto bctype = param.bctype(ig.intersection(),ip.position());
692 auto iplocal_s = geo_in_inside.global(ip.position());
695 auto n_F_local = ig.unitOuterNormal(ip.position());
698 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
699 auto& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
702 RF factor = ip.weight() * geo.integrationElement(ip.position());
707 auto j = param.j(ig.intersection(),ip.position());
710 for (size_type i=0; i<lfsv_s.size(); i++)
711 r_s.accumulate(lfsv_s,i,j * psi_s[i] * factor);
718 for (size_type i=0; i<lfsu_s.size(); i++)
719 u_s += x_s(lfsu_s,i)*phi_s[i];
722 auto b = param.b(cell_inside,iplocal_s);
723 auto normalflux = b*n_F_local;
727 if (normalflux<-1
e-30)
728 DUNE_THROW(Dune::Exception,
729 "Outflow boundary condition on inflow! [b(" 730 << geo.global(ip.position()) <<
") = " 734 auto term1 = u_s * normalflux *factor;
735 for (size_type i=0; i<lfsv_s.size(); i++)
736 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
739 auto o = param.o(ig.intersection(),ip.position());
742 for (size_type i=0; i<lfsv_s.size(); i++)
743 r_s.accumulate(lfsv_s,i,o * psi_s[i] * factor);
750 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
751 auto& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
754 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
755 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
756 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
760 for (size_type i=0; i<lfsu_s.size(); i++)
761 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
764 auto g = param.g(cell_inside,iplocal_s);
767 RF omegaup_s, omegaup_n;
780 auto term1 = (omegaup_s*u_s + omegaup_n*g) * normalflux *factor;
781 for (size_type i=0; i<lfsv_s.size(); i++)
782 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
785 auto term2 = (An_F_s*gradu_s) * factor;
786 for (size_type i=0; i<lfsv_s.size(); i++)
787 r_s.accumulate(lfsv_s,i,-term2 * psi_s[i]);
790 auto term3 = (u_s-g) * factor;
791 for (size_type i=0; i<lfsv_s.size(); i++)
792 r_s.accumulate(lfsv_s,i,term3 * theta * (An_F_s*tgradpsi_s[i]));
795 auto term4 = penalty_factor * (u_s-g) * factor;
796 for (size_type i=0; i<lfsv_s.size(); i++)
797 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
801 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
803 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
807 using RF =
typename LFSV::Traits::FiniteElementType::
808 Traits::LocalBasisType::Traits::RangeFieldType;
809 using size_type =
typename LFSV::Traits::SizeType;
812 const int dim = IG::Entity::dimension;
813 const int order = std::max(
814 lfsu_s.finiteElement().localBasis().order(),
815 lfsv_s.finiteElement().localBasis().order()
819 const auto& cell_inside = ig.inside();
822 auto geo = ig.geometry();
823 auto geo_inside = cell_inside.geometry();
826 auto geo_in_inside = ig.geometryInInside();
829 auto ref_el_inside = referenceElement(geo_inside);
830 auto local_inside = ref_el_inside.position(0,0);
831 auto A_s = param.A(cell_inside,local_inside);
835 auto h_F = geo_inside.volume()/geo.volume();
838 auto n_F = ig.centerUnitOuterNormal();
839 Dune::FieldVector<RF,dim> An_F_s;
843 harmonic_average = An_F_s*n_F;
845 harmonic_average = 1.0;
848 auto order_s = lfsu_s.finiteElement().localBasis().order();
849 auto degree = order_s;
852 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
855 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
858 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
861 auto intorder = intorderadd+quadrature_factor*order;
864 auto bctype = param.bctype(ig.intersection(),ip.position());
871 auto iplocal_s = geo_in_inside.global(ip.position());
874 auto n_F_local = ig.unitOuterNormal(ip.position());
877 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
880 auto factor = ip.weight() * geo.integrationElement(ip.position());
883 auto b = param.b(cell_inside,iplocal_s);
884 auto normalflux = b*n_F_local;
888 if (normalflux<-1
e-30)
889 DUNE_THROW(Dune::Exception,
890 "Outflow boundary condition on inflow! [b(" 891 << geo.global(ip.position()) <<
") = " 892 << b <<
")" << n_F_local <<
" " << normalflux);
895 for (size_type j=0; j<lfsu_s.size(); j++)
896 for (size_type i=0; i<lfsu_s.size(); i++)
897 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * normalflux * factor * phi_s[i]);
903 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
906 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
907 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
910 RF omegaup_s = normalflux>=0.0 ? 1.0 : 0.0;
913 for (size_type j=0; j<lfsu_s.size(); j++)
914 for (size_type i=0; i<lfsu_s.size(); i++)
915 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux * factor * phi_s[i]);
918 for (size_type j=0; j<lfsu_s.size(); j++)
919 for (size_type i=0; i<lfsu_s.size(); i++)
920 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,-(An_F_s*tgradphi_s[j]) * factor * phi_s[i]);
923 for (size_type j=0; j<lfsu_s.size(); j++)
924 for (size_type i=0; i<lfsu_s.size(); i++)
925 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * (An_F_s*tgradphi_s[i]));
928 for (size_type j=0; j<lfsu_s.size(); j++)
929 for (size_type i=0; i<lfsu_s.size(); i++)
930 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,penalty_factor * phi_s[j] * phi_s[i] * factor);
935 template<
typename EG,
typename LFSV,
typename R>
939 using size_type =
typename LFSV::Traits::SizeType;
942 const auto& cell = eg.entity();
945 auto geo = eg.geometry();
948 auto order = lfsv.finiteElement().localBasis().order();
949 auto intorder = intorderadd + 2 * order;
953 auto& phi = cache[order].evaluateFunction(ip.position(),lfsv.finiteElement().localBasis());
956 auto f = param.f(cell,ip.position());
959 auto factor = ip.weight() * geo.integrationElement(ip.position());
960 for (size_type i=0; i<lfsv.size(); i++)
961 r.accumulate(lfsv,i,-f*phi[i]*factor);
978 int quadrature_factor;
981 using LocalBasisType =
typename FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType;
993 std::vector<Cache> cache;
996 void element_size (
const GEO& geo,
typename GEO::ctype& hmin,
typename GEO::ctype hmax)
const 998 using DF =
typename GEO::ctype;
1001 const int dim = GEO::coorddimension;
1004 Dune::FieldVector<DF,dim> x = geo.corner(0);
1006 hmin = hmax = x.two_norm();
1011 Dune::GeometryType gt = geo.type();
1012 for (
int i=0; i<Dune::ReferenceElements<DF,dim>::general(gt).size(dim-1); i++)
1014 Dune::FieldVector<DF,dim> x = geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,0,dim));
1015 x -= geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,1,dim));
1016 hmin = std::min(hmin,x.two_norm());
1017 hmax = std::max(hmax,x.two_norm());
1025 #endif // DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONDG_HH
Definition: convectiondiffusiondg.hh:31
static const int dim
Definition: adaptivity.hh:84
const Entity & e
Definition: localfunctionspace.hh:120
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:29
Definition: convectiondiffusiondg.hh:29
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104
Definition: convectiondiffusiondg.hh:36
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
sparsity pattern generator
Definition: pattern.hh:13
Definition: convectiondiffusiondg.hh:31
Definition: convectiondiffusionparameter.hh:65
Definition: convectiondiffusiondg.hh:31
Definition: convectiondiffusionparameter.hh:65
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:441
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:17
Definition: convectiondiffusiondg.hh:36
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusiondg.hh:190
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
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:103
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Definition: convectiondiffusiondg.hh:183
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:450
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:802
void setTime(Real t)
set time in parameter class
Definition: convectiondiffusiondg.hh:966
Definition: convectiondiffusionparameter.hh:65
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:256
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericaljacobianapply.hh:285
const IG & ig
Definition: constraints.hh:148
Definition: convectiondiffusiondg.hh:54
Type
Definition: convectiondiffusiondg.hh:36
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:936
Type
Definition: convectiondiffusiondg.hh:31
Type
Definition: convectiondiffusionparameter.hh:65
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:622
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:85
Definition: convectiondiffusionparameter.hh:65
Definition: convectiondiffusiondg.hh:34