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>
64 enum {
dim = T::Traits::GridViewType::dimension };
66 using Real =
typename T::Traits::RangeFieldType;
71 enum { doPatternVolume =
true };
72 enum { doPatternSkeleton =
true };
75 enum { doAlphaVolume =
true };
76 enum { doAlphaSkeleton =
true };
77 enum { doAlphaBoundary =
true };
78 enum { doLambdaVolume =
true };
97 param(param_), method(method_), weights(weights_),
98 alpha(alpha_), intorderadd(intorderadd_), quadrature_factor(2),
107 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
108 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 111 using RF =
typename LFSU::Traits::FiniteElementType::
112 Traits::LocalBasisType::Traits::RangeFieldType;
113 using size_type =
typename LFSU::Traits::SizeType;
116 const int dim = EG::Entity::dimension;
117 const int order = std::max(lfsu.finiteElement().localBasis().order(),
118 lfsv.finiteElement().localBasis().order());
121 const auto& cell = eg.entity();
124 auto geo = eg.geometry();
128 auto localcenter = ref_el.position(0,0);
129 auto A = param.A(cell,localcenter);
132 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
133 std::vector<Dune::FieldVector<RF,dim> > gradpsi(lfsv.size());
134 Dune::FieldVector<RF,dim> gradu(0.0);
135 Dune::FieldVector<RF,dim> Agradu(0.0);
138 typename EG::Geometry::JacobianInverseTransposed jac;
141 auto intorder = intorderadd + quadrature_factor * order;
145 auto& phi = cache[order].evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
146 auto& psi = cache[order].evaluateFunction(ip.position(),lfsv.finiteElement().localBasis());
150 for (size_type i=0; i<lfsu.size(); i++)
151 u += x(lfsu,i)*phi[i];
154 auto& js = cache[order].evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
155 auto& js_v = cache[order].evaluateJacobian(ip.position(),lfsv.finiteElement().localBasis());
158 jac = geo.jacobianInverseTransposed(ip.position());
159 for (size_type i=0; i<lfsu.size(); i++)
160 jac.mv(js[i][0],gradphi[i]);
162 for (size_type i=0; i<lfsv.size(); i++)
163 jac.mv(js_v[i][0],gradpsi[i]);
167 for (size_type i=0; i<lfsu.size(); i++)
168 gradu.axpy(x(lfsu,i),gradphi[i]);
174 auto b = param.b(cell,ip.position());
177 auto c = param.c(cell,ip.position());
180 RF factor = ip.weight() * geo.integrationElement(ip.position());
181 for (size_type i=0; i<lfsv.size(); i++)
182 r.accumulate(lfsv,i,( Agradu*gradpsi[i] - u*(b*gradpsi[i]) + c*u*psi[i] )*factor);
187 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
192 using RF =
typename LFSU::Traits::FiniteElementType::
193 Traits::LocalBasisType::Traits::RangeFieldType;
194 using size_type =
typename LFSU::Traits::SizeType;
197 const int dim = EG::Entity::dimension;
198 const int order = std::max(lfsu.finiteElement().localBasis().order(),
199 lfsv.finiteElement().localBasis().order());
202 const auto& cell = eg.entity();
205 auto geo = eg.geometry();
209 auto localcenter = ref_el.position(0,0);
210 auto A = param.A(cell,localcenter);
213 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
214 std::vector<Dune::FieldVector<RF,dim> > Agradphi(lfsu.size());
217 typename EG::Geometry::JacobianInverseTransposed jac;
220 auto intorder = intorderadd + quadrature_factor * order;
224 auto& phi = cache[order].evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
227 auto& js = cache[order].evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
230 jac = geo.jacobianInverseTransposed(ip.position());
231 for (size_type i=0; i<lfsu.size(); i++)
233 jac.mv(js[i][0],gradphi[i]);
234 A.mv(gradphi[i],Agradphi[i]);
238 auto b = param.b(cell,ip.position());
241 auto c = param.c(cell,ip.position());
244 auto factor = ip.weight() * geo.integrationElement(ip.position());
245 for (size_type j=0; j<lfsu.size(); j++)
246 for (size_type i=0; i<lfsu.size(); i++)
247 mat.accumulate(lfsu,i,lfsu,j,( Agradphi[j]*gradphi[i] - phi[j]*(b*gradphi[i]) + c*phi[j]*phi[i] )*factor);
253 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
255 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
256 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
257 R& r_s, R& r_n)
const 260 using RF =
typename LFSV::Traits::FiniteElementType::
261 Traits::LocalBasisType::Traits::RangeFieldType;
262 using size_type =
typename LFSV::Traits::SizeType;
265 const int dim = IG::dimension;
266 const int order = std::max(
267 std::max(lfsu_s.finiteElement().localBasis().order(),
268 lfsu_n.finiteElement().localBasis().order()),
269 std::max(lfsv_s.finiteElement().localBasis().order(),
270 lfsv_n.finiteElement().localBasis().order())
274 const auto& cell_inside = ig.inside();
275 const auto& cell_outside = ig.outside();
278 auto geo = ig.geometry();
279 auto geo_inside = cell_inside.geometry();
280 auto geo_outside = cell_outside.geometry();
283 auto geo_in_inside = ig.geometryInInside();
284 auto geo_in_outside = ig.geometryInOutside();
289 auto local_inside = ref_el_inside.position(0,0);
290 auto local_outside = ref_el_outside.position(0,0);
291 auto A_s = param.A(cell_inside,local_inside);
292 auto A_n = param.A(cell_outside,local_outside);
296 auto h_F = std::min(geo_inside.volume(),geo_outside.volume())/geo.volume();
299 auto n_F = ig.centerUnitOuterNormal();
300 Dune::FieldVector<RF,dim> An_F_s;
302 Dune::FieldVector<RF,dim> An_F_n;
308 RF harmonic_average(0.0);
311 RF delta_s = (An_F_s*n_F);
312 RF delta_n = (An_F_n*n_F);
313 omega_s = delta_n/(delta_s+delta_n+1
e-20);
314 omega_n = delta_s/(delta_s+delta_n+1
e-20);
315 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-20);
319 omega_s = omega_n = 0.5;
320 harmonic_average = 1.0;
324 auto order_s = lfsu_s.finiteElement().localBasis().order();
325 auto order_n = lfsu_n.finiteElement().localBasis().order();
326 auto degree = std::max( order_s, order_n );
329 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
332 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
333 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
334 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
335 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_n(lfsv_n.size());
336 Dune::FieldVector<RF,dim> gradu_s(0.0);
337 Dune::FieldVector<RF,dim> gradu_n(0.0);
340 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
343 auto intorder = intorderadd+quadrature_factor*order;
347 auto n_F_local = ig.unitOuterNormal(ip.position());
350 auto iplocal_s = geo_in_inside.global(ip.position());
351 auto iplocal_n = geo_in_outside.global(ip.position());
354 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
355 auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
356 auto& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
357 auto& psi_n = cache[order_n].evaluateFunction(iplocal_n,lfsv_n.finiteElement().localBasis());
361 for (size_type i=0; i<lfsu_s.size(); i++)
362 u_s += x_s(lfsu_s,i)*phi_s[i];
364 for (size_type i=0; i<lfsu_n.size(); i++)
365 u_n += x_n(lfsu_n,i)*phi_n[i];
368 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
369 auto& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
370 auto& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
371 auto& gradpsi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsv_n.finiteElement().localBasis());
374 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
375 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
376 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
377 jac = geo_outside.jacobianInverseTransposed(iplocal_n);
378 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
379 for (size_type i=0; i<lfsv_n.size(); i++) jac.mv(gradpsi_n[i][0],tgradpsi_n[i]);
383 for (size_type i=0; i<lfsu_s.size(); i++)
384 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
386 for (size_type i=0; i<lfsu_n.size(); i++)
387 gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
390 auto b = param.b(cell_inside,iplocal_s);
391 auto normalflux = b*n_F_local;
392 RF omegaup_s, omegaup_n;
405 auto factor = ip.weight()*geo.integrationElement(ip.position());
408 auto term1 = (omegaup_s*u_s + omegaup_n*u_n) * normalflux *factor;
409 for (size_type i=0; i<lfsv_s.size(); i++)
410 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
411 for (size_type i=0; i<lfsv_n.size(); i++)
412 r_n.accumulate(lfsv_n,i,-term1 * psi_n[i]);
415 auto term2 = -(omega_s*(An_F_s*gradu_s) + omega_n*(An_F_n*gradu_n)) * factor;
416 for (size_type i=0; i<lfsv_s.size(); i++)
417 r_s.accumulate(lfsv_s,i,term2 * psi_s[i]);
418 for (size_type i=0; i<lfsv_n.size(); i++)
419 r_n.accumulate(lfsv_n,i,-term2 * psi_n[i]);
422 auto term3 = (u_s-u_n) * factor;
423 for (size_type i=0; i<lfsv_s.size(); i++)
424 r_s.accumulate(lfsv_s,i,term3 * theta * omega_s * (An_F_s*tgradpsi_s[i]));
425 for (size_type i=0; i<lfsv_n.size(); i++)
426 r_n.accumulate(lfsv_n,i,term3 * theta * omega_n * (An_F_n*tgradpsi_n[i]));
429 auto term4 = penalty_factor * (u_s-u_n) * factor;
430 for (size_type i=0; i<lfsv_s.size(); i++)
431 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
432 for (size_type i=0; i<lfsv_n.size(); i++)
433 r_n.accumulate(lfsv_n,i,-term4 * psi_n[i]);
437 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
439 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
440 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
441 M& mat_ss, M& mat_sn,
442 M& mat_ns, M& mat_nn)
const 445 using RF =
typename LFSV::Traits::FiniteElementType::
446 Traits::LocalBasisType::Traits::RangeFieldType;
447 using size_type =
typename LFSV::Traits::SizeType;
450 const int dim = IG::dimension;
451 const int order = std::max(
452 std::max(lfsu_s.finiteElement().localBasis().order(),
453 lfsu_n.finiteElement().localBasis().order()),
454 std::max(lfsv_s.finiteElement().localBasis().order(),
455 lfsv_n.finiteElement().localBasis().order())
459 const auto& cell_inside = ig.inside();
460 const auto& cell_outside = ig.outside();
463 auto geo = ig.geometry();
464 auto geo_inside = cell_inside.geometry();
465 auto geo_outside = cell_outside.geometry();
468 auto geo_in_inside = ig.geometryInInside();
469 auto geo_in_outside = ig.geometryInOutside();
474 auto local_inside = ref_el_inside.position(0,0);
475 auto local_outside = ref_el_outside.position(0,0);
476 auto A_s = param.A(cell_inside,local_inside);
477 auto A_n = param.A(cell_outside,local_outside);
481 auto h_F = std::min(geo_inside.volume(),geo_outside.volume())/geo.volume();
484 auto n_F = ig.centerUnitOuterNormal();
485 Dune::FieldVector<RF,dim> An_F_s;
487 Dune::FieldVector<RF,dim> An_F_n;
493 RF harmonic_average(0.0);
496 RF delta_s = (An_F_s*n_F);
497 RF delta_n = (An_F_n*n_F);
498 omega_s = delta_n/(delta_s+delta_n+1
e-20);
499 omega_n = delta_s/(delta_s+delta_n+1
e-20);
500 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-20);
504 omega_s = omega_n = 0.5;
505 harmonic_average = 1.0;
509 auto order_s = lfsu_s.finiteElement().localBasis().order();
510 auto order_n = lfsu_n.finiteElement().localBasis().order();
511 auto degree = std::max( order_s, order_n );
514 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
517 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
518 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
521 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
524 const int intorder = intorderadd+quadrature_factor*order;
528 auto n_F_local = ig.unitOuterNormal(ip.position());
531 auto iplocal_s = geo_in_inside.global(ip.position());
532 auto iplocal_n = geo_in_outside.global(ip.position());
535 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
536 auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
539 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
540 auto& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
543 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
544 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
545 jac = geo_outside.jacobianInverseTransposed(iplocal_n);
546 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
549 auto b = param.b(cell_inside,iplocal_s);
550 auto normalflux = b*n_F_local;
551 RF omegaup_s, omegaup_n;
564 auto factor = ip.weight() * geo.integrationElement(ip.position());
565 auto ipfactor = penalty_factor * factor;
568 for (size_type j=0; j<lfsu_s.size(); j++) {
569 auto temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
570 for (size_type i=0; i<lfsu_s.size(); i++) {
571 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux *factor * phi_s[i]);
572 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,temp1 * phi_s[i]);
573 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
574 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * ipfactor * phi_s[i]);
577 for (size_type j=0; j<lfsu_n.size(); j++) {
578 auto temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
579 for (size_type i=0; i<lfsu_s.size(); i++) {
580 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,omegaup_n * phi_n[j] * normalflux *factor * phi_s[i]);
581 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,temp1 * phi_s[i]);
582 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
583 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * ipfactor * phi_s[i]);
586 for (size_type j=0; j<lfsu_s.size(); j++) {
587 auto temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
588 for (size_type i=0; i<lfsu_n.size(); i++) {
589 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-omegaup_s * phi_s[j] * normalflux *factor * phi_n[i]);
590 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-temp1 * phi_n[i]);
591 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,phi_s[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
592 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-phi_s[j] * ipfactor * phi_n[i]);
595 for (size_type j=0; j<lfsu_n.size(); j++) {
596 auto temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
597 for (size_type i=0; i<lfsu_n.size(); i++) {
598 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-omegaup_n * phi_n[j] * normalflux *factor * phi_n[i]);
599 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-temp1 * phi_n[i]);
600 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
601 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,phi_n[j] * ipfactor * phi_n[i]);
609 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
611 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
615 using RF =
typename LFSV::Traits::FiniteElementType::
616 Traits::LocalBasisType::Traits::RangeFieldType;
617 using size_type =
typename LFSV::Traits::SizeType;
620 const int dim = IG::dimension;
621 const int order = std::max(
622 lfsu_s.finiteElement().localBasis().order(),
623 lfsv_s.finiteElement().localBasis().order()
627 const auto& cell_inside = ig.inside();
630 auto geo = ig.geometry();
631 auto geo_inside = cell_inside.geometry();
634 auto geo_in_inside = ig.geometryInInside();
638 auto local_inside = ref_el_inside.position(0,0);
639 auto A_s = param.A(cell_inside,local_inside);
643 auto h_F = geo_inside.volume()/geo.volume();
646 auto n_F = ig.centerUnitOuterNormal();
647 Dune::FieldVector<RF,dim> An_F_s;
651 harmonic_average = An_F_s*n_F;
653 harmonic_average = 1.0;
656 auto order_s = lfsu_s.finiteElement().localBasis().order();
657 auto degree = order_s;
660 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
663 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
664 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
665 Dune::FieldVector<RF,dim> gradu_s(0.0);
668 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
671 auto intorder = intorderadd+quadrature_factor*order;
674 auto bctype = param.bctype(ig.intersection(),ip.position());
680 auto iplocal_s = geo_in_inside.global(ip.position());
683 auto n_F_local = ig.unitOuterNormal(ip.position());
686 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
687 auto& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
690 RF factor = ip.weight() * geo.integrationElement(ip.position());
695 auto j = param.j(ig.intersection(),ip.position());
698 for (size_type i=0; i<lfsv_s.size(); i++)
699 r_s.accumulate(lfsv_s,i,j * psi_s[i] * factor);
706 for (size_type i=0; i<lfsu_s.size(); i++)
707 u_s += x_s(lfsu_s,i)*phi_s[i];
710 auto b = param.b(cell_inside,iplocal_s);
711 auto normalflux = b*n_F_local;
715 if (normalflux<-1
e-30)
716 DUNE_THROW(Dune::Exception,
717 "Outflow boundary condition on inflow! [b(" 718 << geo.global(ip.position()) <<
") = " 722 auto term1 = u_s * normalflux *factor;
723 for (size_type i=0; i<lfsv_s.size(); i++)
724 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
727 auto o = param.o(ig.intersection(),ip.position());
730 for (size_type i=0; i<lfsv_s.size(); i++)
731 r_s.accumulate(lfsv_s,i,o * psi_s[i] * factor);
738 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
739 auto& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
742 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
743 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
744 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
748 for (size_type i=0; i<lfsu_s.size(); i++)
749 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
752 auto g = param.g(cell_inside,iplocal_s);
755 RF omegaup_s, omegaup_n;
768 auto term1 = (omegaup_s*u_s + omegaup_n*g) * normalflux *factor;
769 for (size_type i=0; i<lfsv_s.size(); i++)
770 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
773 auto term2 = (An_F_s*gradu_s) * factor;
774 for (size_type i=0; i<lfsv_s.size(); i++)
775 r_s.accumulate(lfsv_s,i,-term2 * psi_s[i]);
778 auto term3 = (u_s-g) * factor;
779 for (size_type i=0; i<lfsv_s.size(); i++)
780 r_s.accumulate(lfsv_s,i,term3 * theta * (An_F_s*tgradpsi_s[i]));
783 auto term4 = penalty_factor * (u_s-g) * factor;
784 for (size_type i=0; i<lfsv_s.size(); i++)
785 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
789 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
791 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
795 using RF =
typename LFSV::Traits::FiniteElementType::
796 Traits::LocalBasisType::Traits::RangeFieldType;
797 using size_type =
typename LFSV::Traits::SizeType;
800 const int dim = IG::dimension;
801 const int order = std::max(
802 lfsu_s.finiteElement().localBasis().order(),
803 lfsv_s.finiteElement().localBasis().order()
807 const auto& cell_inside = ig.inside();
810 auto geo = ig.geometry();
811 auto geo_inside = cell_inside.geometry();
814 auto geo_in_inside = ig.geometryInInside();
818 auto local_inside = ref_el_inside.position(0,0);
819 auto A_s = param.A(cell_inside,local_inside);
823 auto h_F = geo_inside.volume()/geo.volume();
826 auto n_F = ig.centerUnitOuterNormal();
827 Dune::FieldVector<RF,dim> An_F_s;
831 harmonic_average = An_F_s*n_F;
833 harmonic_average = 1.0;
836 auto order_s = lfsu_s.finiteElement().localBasis().order();
837 auto degree = order_s;
840 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
843 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
846 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
849 auto intorder = intorderadd+quadrature_factor*order;
852 auto bctype = param.bctype(ig.intersection(),ip.position());
859 auto iplocal_s = geo_in_inside.global(ip.position());
862 auto n_F_local = ig.unitOuterNormal(ip.position());
865 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
868 auto factor = ip.weight() * geo.integrationElement(ip.position());
871 auto b = param.b(cell_inside,iplocal_s);
872 auto normalflux = b*n_F_local;
876 if (normalflux<-1
e-30)
877 DUNE_THROW(Dune::Exception,
878 "Outflow boundary condition on inflow! [b(" 879 << geo.global(ip.position()) <<
") = " 880 << b <<
")" << n_F_local <<
" " << normalflux);
883 for (size_type j=0; j<lfsu_s.size(); j++)
884 for (size_type i=0; i<lfsu_s.size(); i++)
885 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * normalflux * factor * phi_s[i]);
891 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
894 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
895 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
898 RF omegaup_s, omegaup_n;
911 for (size_type j=0; j<lfsu_s.size(); j++)
912 for (size_type i=0; i<lfsu_s.size(); i++)
913 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux * factor * phi_s[i]);
916 for (size_type j=0; j<lfsu_s.size(); j++)
917 for (size_type i=0; i<lfsu_s.size(); i++)
918 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,-(An_F_s*tgradphi_s[j]) * factor * phi_s[i]);
921 for (size_type j=0; j<lfsu_s.size(); j++)
922 for (size_type i=0; i<lfsu_s.size(); i++)
923 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * (An_F_s*tgradphi_s[i]));
926 for (size_type j=0; j<lfsu_s.size(); j++)
927 for (size_type i=0; i<lfsu_s.size(); i++)
928 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,penalty_factor * phi_s[j] * phi_s[i] * factor);
933 template<
typename EG,
typename LFSV,
typename R>
937 using size_type =
typename LFSV::Traits::SizeType;
940 const auto& cell = eg.entity();
943 auto geo = eg.geometry();
946 auto order = lfsv.finiteElement().localBasis().order();
947 auto intorder = intorderadd + 2 * order;
951 auto& phi = cache[order].evaluateFunction(ip.position(),lfsv.finiteElement().localBasis());
954 auto f = param.f(cell,ip.position());
957 auto factor = ip.weight() * geo.integrationElement(ip.position());
958 for (size_type i=0; i<lfsv.size(); i++)
959 r.accumulate(lfsv,i,-f*phi[i]*factor);
976 int quadrature_factor;
979 using LocalBasisType =
typename FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType;
991 std::vector<Cache> cache;
994 void element_size (
const GEO& geo,
typename GEO::ctype& hmin,
typename GEO::ctype hmax)
const 996 using DF =
typename GEO::ctype;
999 const int dim = GEO::coorddimension;
1002 Dune::FieldVector<DF,dim> x = geo.corner(0);
1004 hmin = hmax = x.two_norm();
1009 Dune::GeometryType gt = geo.type();
1010 for (
int i=0; i<Dune::ReferenceElements<DF,dim>::general(gt).size(dim-1); i++)
1012 Dune::FieldVector<DF,dim> x = geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,0,dim));
1013 x -= geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,1,dim));
1014 hmin = std::min(hmin,x.two_norm());
1015 hmax = std::max(hmax,x.two_norm());
1023 #endif // DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONDG_HH 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
Definition: convectiondiffusiondg.hh:35
const IG & ig
Definition: constraints.hh:148
Definition: convectiondiffusionparameter.hh:65
static const int dim
Definition: adaptivity.hh:83
Type
Definition: convectiondiffusionparameter.hh:65
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
sparsity pattern generator
Definition: pattern.hh:29
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:610
Definition: convectiondiffusionparameter.hh:65
Definition: convectiondiffusiondg.hh:32
Type
Definition: convectiondiffusiondg.hh:32
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericaljacobianapply.hh:285
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:88
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:108
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume() ...
Definition: numericaljacobianapply.hh:32
Type
Definition: convectiondiffusiondg.hh:37
Definition: convectiondiffusiondg.hh:37
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton() ...
Definition: numericaljacobianapply.hh:180
Definition: convectiondiffusiondg.hh:32
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:17
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusiondg.hh:188
Definition: convectiondiffusiondg.hh:30
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:934
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
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:438
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:254
sparsity pattern generator
Definition: pattern.hh:13
const Entity & e
Definition: localfunctionspace.hh:111
void setTime(Real t)
set time in parameter class
Definition: convectiondiffusiondg.hh:964
Definition: convectiondiffusiondg.hh:32
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:790
Definition: convectiondiffusiondg.hh:37
Definition: convectiondiffusiondg.hh:55
Definition: convectiondiffusionparameter.hh:65
Definition: convectiondiffusionparameter.hh:65
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
Default flags for all local operators.
Definition: flags.hh:18