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>
61 enum { dim = T::Traits::GridViewType::dimension };
63 using Real =
typename T::Traits::RangeFieldType;
95 , intorderadd(intorderadd_)
96 , quadrature_factor(2)
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();
125 auto ref_el = referenceElement(geo);
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());
132 Dune::FieldVector<RF,dim> gradu(0.0);
133 Dune::FieldVector<RF,dim> Agradu(0.0);
136 typename EG::Geometry::JacobianInverseTransposed jac;
139 auto intorder = intorderadd + quadrature_factor * order;
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>
198 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
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();
219 auto ref_el = referenceElement(geo);
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;
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>
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;
283 const int order = std::max(
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();
304 auto ref_el_inside = referenceElement(geo_inside);
305 auto ref_el_outside = referenceElement(geo_outside);
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();
317 Dune::FieldVector<RF,dim> An_F_s;
319 Dune::FieldVector<RF,dim> An_F_n;
325 RF harmonic_average(0.0);
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+1
e-20);
331 omega_n = delta_s/(delta_s+delta_n+1
e-20);
332 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-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();
343 auto degree = std::max( order_s, order_n );
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());
353 Dune::FieldVector<RF,dim> gradu_s(0.0);
354 Dune::FieldVector<RF,dim> gradu_n(0.0);
357 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
360 auto intorder = intorderadd+quadrature_factor*order;
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);
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+1
e-20);
378 omega_n = delta_s/(delta_s+delta_n+1
e-20);
379 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-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>
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
483 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
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;
497 const int order = std::max(
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();
518 auto ref_el_inside = referenceElement(geo_inside);
519 auto ref_el_outside = referenceElement(geo_outside);
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();
531 Dune::FieldVector<RF,dim> An_F_s;
533 Dune::FieldVector<RF,dim> An_F_n;
539 RF harmonic_average(0.0);
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+1
e-20);
545 omega_n = delta_s/(delta_s+delta_n+1
e-20);
546 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-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();
557 auto degree = std::max( order_s, order_n );
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;
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);
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+1
e-20);
588 omega_n = delta_s/(delta_s+delta_n+1
e-20);
589 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-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>
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;
695 const int order = std::max(
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();
711 auto ref_el_inside = referenceElement(geo_inside);
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();
721 Dune::FieldVector<RF,dim> An_F_s;
725 harmonic_average = An_F_s*n_F;
727 harmonic_average = 1.0;
730 auto order_s = lfsu_s.finiteElement().localBasis().order();
731 auto degree = order_s;
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());
739 Dune::FieldVector<RF,dim> gradu_s(0.0);
742 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
745 auto intorder = intorderadd+quadrature_factor*order;
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);
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());
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());
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;
802 if (normalflux<-1
e-30)
803 DUNE_THROW(Dune::Exception,
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);
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>
885 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
892 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
894 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
901 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
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;
913 const int order = std::max(
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();
929 auto ref_el_inside = referenceElement(geo_inside);
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();
939 Dune::FieldVector<RF,dim> An_F_s;
943 harmonic_average = An_F_s*n_F;
945 harmonic_average = 1.0;
948 auto order_s = lfsu_s.finiteElement().localBasis().order();
949 auto degree = order_s;
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;
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);
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());
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;
1000 if (normalflux<-1
e-30)
1001 DUNE_THROW(Dune::Exception,
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>
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;
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);
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;
1116 Dune::FieldVector<DF,dim> x = geo.corner(0);
1118 hmin = hmax = x.two_norm();
1123 Dune::GeometryType gt = geo.type();
1124 for (
int i=0; i<Dune::ReferenceElements<DF,dim>::general(gt).size(dim-1); i++)
1126 Dune::FieldVector<DF,dim> x = geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,0,dim));
1127 x -= geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,1,dim));
1128 hmin = std::min(hmin,x.two_norm());
1129 hmax = std::max(hmax,x.two_norm());
1137 #endif // DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONDG_HH