2 #ifndef DUNE_PDELAB_CONVECTIONDIFFUSIONDG_HH 3 #define DUNE_PDELAB_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 };
90 param(param_), method(method_), weights(weights_),
91 alpha(alpha_), intorderadd(intorderadd_), quadrature_factor(2),
99 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
100 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 103 using RF =
typename LFSU::Traits::FiniteElementType::
104 Traits::LocalBasisType::Traits::RangeFieldType;
105 using size_type =
typename LFSU::Traits::SizeType;
108 const int dim = EG::Entity::dimension;
109 const int order = std::max(lfsu.finiteElement().localBasis().order(),
110 lfsv.finiteElement().localBasis().order());
113 const auto& cell = eg.entity();
116 auto geo = eg.geometry();
120 auto localcenter = ref_el.position(0,0);
121 auto A = param.A(cell,localcenter);
124 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
125 std::vector<Dune::FieldVector<RF,dim> > gradpsi(lfsv.size());
126 Dune::FieldVector<RF,dim> gradu(0.0);
127 Dune::FieldVector<RF,dim> Agradu(0.0);
130 typename EG::Geometry::JacobianInverseTransposed jac;
133 auto intorder = intorderadd + quadrature_factor * order;
137 auto& phi = cache[order].evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
138 auto& psi = cache[order].evaluateFunction(ip.position(),lfsv.finiteElement().localBasis());
142 for (size_type i=0; i<lfsu.size(); i++)
143 u += x(lfsu,i)*phi[i];
146 auto& js = cache[order].evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
147 auto& js_v = cache[order].evaluateJacobian(ip.position(),lfsv.finiteElement().localBasis());
150 jac = geo.jacobianInverseTransposed(ip.position());
151 for (size_type i=0; i<lfsu.size(); i++)
152 jac.mv(js[i][0],gradphi[i]);
154 for (size_type i=0; i<lfsv.size(); i++)
155 jac.mv(js_v[i][0],gradpsi[i]);
159 for (size_type i=0; i<lfsu.size(); i++)
160 gradu.axpy(x(lfsu,i),gradphi[i]);
166 auto b = param.b(cell,ip.position());
169 auto c = param.c(cell,ip.position());
172 RF factor = ip.weight() * geo.integrationElement(ip.position());
173 for (size_type i=0; i<lfsv.size(); i++)
174 r.accumulate(lfsv,i,( Agradu*gradpsi[i] - u*(b*gradpsi[i]) + c*u*psi[i] )*factor);
179 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
184 using RF =
typename LFSU::Traits::FiniteElementType::
185 Traits::LocalBasisType::Traits::RangeFieldType;
186 using size_type =
typename LFSU::Traits::SizeType;
189 const int dim = EG::Entity::dimension;
190 const int order = std::max(lfsu.finiteElement().localBasis().order(),
191 lfsv.finiteElement().localBasis().order());
194 const auto& cell = eg.entity();
197 auto geo = eg.geometry();
201 auto localcenter = ref_el.position(0,0);
202 auto A = param.A(cell,localcenter);
205 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
206 std::vector<Dune::FieldVector<RF,dim> > Agradphi(lfsu.size());
209 typename EG::Geometry::JacobianInverseTransposed jac;
212 auto intorder = intorderadd + quadrature_factor * order;
216 auto& phi = cache[order].evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
219 auto& js = cache[order].evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
222 jac = geo.jacobianInverseTransposed(ip.position());
223 for (size_type i=0; i<lfsu.size(); i++)
225 jac.mv(js[i][0],gradphi[i]);
226 A.mv(gradphi[i],Agradphi[i]);
230 auto b = param.b(cell,ip.position());
233 auto c = param.c(cell,ip.position());
236 auto factor = ip.weight() * geo.integrationElement(ip.position());
237 for (size_type j=0; j<lfsu.size(); j++)
238 for (size_type i=0; i<lfsu.size(); i++)
239 mat.accumulate(lfsu,i,lfsu,j,( Agradphi[j]*gradphi[i] - phi[j]*(b*gradphi[i]) + c*phi[j]*phi[i] )*factor);
245 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
247 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
248 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
249 R& r_s, R& r_n)
const 252 using RF =
typename LFSV::Traits::FiniteElementType::
253 Traits::LocalBasisType::Traits::RangeFieldType;
254 using size_type =
typename LFSV::Traits::SizeType;
257 const int dim = IG::dimension;
258 const int order = std::max(
259 std::max(lfsu_s.finiteElement().localBasis().order(),
260 lfsu_n.finiteElement().localBasis().order()),
261 std::max(lfsv_s.finiteElement().localBasis().order(),
262 lfsv_n.finiteElement().localBasis().order())
266 const auto& cell_inside = ig.inside();
267 const auto& cell_outside = ig.outside();
270 auto geo = ig.geometry();
271 auto geo_inside = cell_inside.geometry();
272 auto geo_outside = cell_outside.geometry();
275 auto geo_in_inside = ig.geometryInInside();
276 auto geo_in_outside = ig.geometryInOutside();
281 auto local_inside = ref_el_inside.position(0,0);
282 auto local_outside = ref_el_outside.position(0,0);
283 auto A_s = param.A(cell_inside,local_inside);
284 auto A_n = param.A(cell_outside,local_outside);
287 auto h_F = std::min(geo_inside.volume(),geo_outside.volume())/geo.volume();
290 auto n_F = ig.centerUnitOuterNormal();
291 Dune::FieldVector<RF,dim> An_F_s;
293 Dune::FieldVector<RF,dim> An_F_n;
299 RF harmonic_average(0.0);
302 RF delta_s = (An_F_s*n_F);
303 RF delta_n = (An_F_n*n_F);
304 omega_s = delta_n/(delta_s+delta_n+1
e-20);
305 omega_n = delta_s/(delta_s+delta_n+1
e-20);
306 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-20);
310 omega_s = omega_n = 0.5;
311 harmonic_average = 1.0;
315 auto order_s = lfsu_s.finiteElement().localBasis().order();
316 auto order_n = lfsu_n.finiteElement().localBasis().order();
317 auto degree = std::max( order_s, order_n );
320 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
323 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
324 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
325 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
326 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_n(lfsv_n.size());
327 Dune::FieldVector<RF,dim> gradu_s(0.0);
328 Dune::FieldVector<RF,dim> gradu_n(0.0);
331 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
334 auto intorder = intorderadd+quadrature_factor*order;
338 auto n_F_local = ig.unitOuterNormal(ip.position());
341 auto iplocal_s = geo_in_inside.global(ip.position());
342 auto iplocal_n = geo_in_outside.global(ip.position());
345 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
346 auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
347 auto& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
348 auto& psi_n = cache[order_n].evaluateFunction(iplocal_n,lfsv_n.finiteElement().localBasis());
352 for (size_type i=0; i<lfsu_s.size(); i++)
353 u_s += x_s(lfsu_s,i)*phi_s[i];
355 for (size_type i=0; i<lfsu_n.size(); i++)
356 u_n += x_n(lfsu_n,i)*phi_n[i];
359 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
360 auto& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
361 auto& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
362 auto& gradpsi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsv_n.finiteElement().localBasis());
365 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
366 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
367 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
368 jac = geo_outside.jacobianInverseTransposed(iplocal_n);
369 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
370 for (size_type i=0; i<lfsv_n.size(); i++) jac.mv(gradpsi_n[i][0],tgradpsi_n[i]);
374 for (size_type i=0; i<lfsu_s.size(); i++)
375 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
377 for (size_type i=0; i<lfsu_n.size(); i++)
378 gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
381 auto b = param.b(cell_inside,iplocal_s);
382 auto normalflux = b*n_F_local;
383 RF omegaup_s, omegaup_n;
396 auto factor = ip.weight()*geo.integrationElement(ip.position());
399 auto term1 = (omegaup_s*u_s + omegaup_n*u_n) * normalflux *factor;
400 for (size_type i=0; i<lfsv_s.size(); i++)
401 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
402 for (size_type i=0; i<lfsv_n.size(); i++)
403 r_n.accumulate(lfsv_n,i,-term1 * psi_n[i]);
406 auto term2 = -(omega_s*(An_F_s*gradu_s) + omega_n*(An_F_n*gradu_n)) * factor;
407 for (size_type i=0; i<lfsv_s.size(); i++)
408 r_s.accumulate(lfsv_s,i,term2 * psi_s[i]);
409 for (size_type i=0; i<lfsv_n.size(); i++)
410 r_n.accumulate(lfsv_n,i,-term2 * psi_n[i]);
413 auto term3 = (u_s-u_n) * factor;
414 for (size_type i=0; i<lfsv_s.size(); i++)
415 r_s.accumulate(lfsv_s,i,term3 * theta * omega_s * (An_F_s*tgradpsi_s[i]));
416 for (size_type i=0; i<lfsv_n.size(); i++)
417 r_n.accumulate(lfsv_n,i,term3 * theta * omega_n * (An_F_n*tgradpsi_n[i]));
420 auto term4 = penalty_factor * (u_s-u_n) * factor;
421 for (size_type i=0; i<lfsv_s.size(); i++)
422 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
423 for (size_type i=0; i<lfsv_n.size(); i++)
424 r_n.accumulate(lfsv_n,i,-term4 * psi_n[i]);
428 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
430 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
431 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
432 M& mat_ss, M& mat_sn,
433 M& mat_ns, M& mat_nn)
const 436 using RF =
typename LFSV::Traits::FiniteElementType::
437 Traits::LocalBasisType::Traits::RangeFieldType;
438 using size_type =
typename LFSV::Traits::SizeType;
441 const int dim = IG::dimension;
442 const int order = std::max(
443 std::max(lfsu_s.finiteElement().localBasis().order(),
444 lfsu_n.finiteElement().localBasis().order()),
445 std::max(lfsv_s.finiteElement().localBasis().order(),
446 lfsv_n.finiteElement().localBasis().order())
450 const auto& cell_inside = ig.inside();
451 const auto& cell_outside = ig.outside();
454 auto geo = ig.geometry();
455 auto geo_inside = cell_inside.geometry();
456 auto geo_outside = cell_outside.geometry();
459 auto geo_in_inside = ig.geometryInInside();
460 auto geo_in_outside = ig.geometryInOutside();
465 auto local_inside = ref_el_inside.position(0,0);
466 auto local_outside = ref_el_outside.position(0,0);
467 auto A_s = param.A(cell_inside,local_inside);
468 auto A_n = param.A(cell_outside,local_outside);
471 auto h_F = std::min(geo_inside.volume(),geo_outside.volume())/geo.volume();
474 auto n_F = ig.centerUnitOuterNormal();
475 Dune::FieldVector<RF,dim> An_F_s;
477 Dune::FieldVector<RF,dim> An_F_n;
483 RF harmonic_average(0.0);
486 RF delta_s = (An_F_s*n_F);
487 RF delta_n = (An_F_n*n_F);
488 omega_s = delta_n/(delta_s+delta_n+1
e-20);
489 omega_n = delta_s/(delta_s+delta_n+1
e-20);
490 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-20);
494 omega_s = omega_n = 0.5;
495 harmonic_average = 1.0;
499 auto order_s = lfsu_s.finiteElement().localBasis().order();
500 auto order_n = lfsu_n.finiteElement().localBasis().order();
501 auto degree = std::max( order_s, order_n );
504 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
507 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
508 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
511 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
514 const int intorder = intorderadd+quadrature_factor*order;
518 auto n_F_local = ig.unitOuterNormal(ip.position());
521 auto iplocal_s = geo_in_inside.global(ip.position());
522 auto iplocal_n = geo_in_outside.global(ip.position());
525 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
526 auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
529 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
530 auto& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
533 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
534 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
535 jac = geo_outside.jacobianInverseTransposed(iplocal_n);
536 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
539 auto b = param.b(cell_inside,iplocal_s);
540 auto normalflux = b*n_F_local;
541 RF omegaup_s, omegaup_n;
554 auto factor = ip.weight() * geo.integrationElement(ip.position());
555 auto ipfactor = penalty_factor * factor;
558 for (size_type j=0; j<lfsu_s.size(); j++) {
559 auto temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
560 for (size_type i=0; i<lfsu_s.size(); i++) {
561 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux *factor * phi_s[i]);
562 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,temp1 * phi_s[i]);
563 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
564 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * ipfactor * phi_s[i]);
567 for (size_type j=0; j<lfsu_n.size(); j++) {
568 auto temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
569 for (size_type i=0; i<lfsu_s.size(); i++) {
570 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,omegaup_n * phi_n[j] * normalflux *factor * phi_s[i]);
571 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,temp1 * phi_s[i]);
572 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
573 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * ipfactor * phi_s[i]);
576 for (size_type j=0; j<lfsu_s.size(); j++) {
577 auto temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
578 for (size_type i=0; i<lfsu_n.size(); i++) {
579 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-omegaup_s * phi_s[j] * normalflux *factor * phi_n[i]);
580 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-temp1 * phi_n[i]);
581 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,phi_s[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
582 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-phi_s[j] * ipfactor * phi_n[i]);
585 for (size_type j=0; j<lfsu_n.size(); j++) {
586 auto temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
587 for (size_type i=0; i<lfsu_n.size(); i++) {
588 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-omegaup_n * phi_n[j] * normalflux *factor * phi_n[i]);
589 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-temp1 * phi_n[i]);
590 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
591 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,phi_n[j] * ipfactor * phi_n[i]);
599 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
601 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
605 using RF =
typename LFSV::Traits::FiniteElementType::
606 Traits::LocalBasisType::Traits::RangeFieldType;
607 using size_type =
typename LFSV::Traits::SizeType;
610 const int dim = IG::dimension;
611 const int order = std::max(
612 lfsu_s.finiteElement().localBasis().order(),
613 lfsv_s.finiteElement().localBasis().order()
617 const auto& cell_inside = ig.inside();
620 auto geo = ig.geometry();
621 auto geo_inside = cell_inside.geometry();
624 auto geo_in_inside = ig.geometryInInside();
628 auto local_inside = ref_el_inside.position(0,0);
629 auto A_s = param.A(cell_inside,local_inside);
632 auto h_F = geo_inside.volume()/geo.volume();
635 auto n_F = ig.centerUnitOuterNormal();
636 Dune::FieldVector<RF,dim> An_F_s;
640 harmonic_average = An_F_s*n_F;
642 harmonic_average = 1.0;
645 auto order_s = lfsu_s.finiteElement().localBasis().order();
646 auto degree = order_s;
649 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
652 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
653 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
654 Dune::FieldVector<RF,dim> gradu_s(0.0);
657 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
660 auto intorder = intorderadd+quadrature_factor*order;
663 auto bctype = param.bctype(ig.intersection(),ip.position());
669 auto iplocal_s = geo_in_inside.global(ip.position());
672 auto n_F_local = ig.unitOuterNormal(ip.position());
675 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
676 auto& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
679 RF factor = ip.weight() * geo.integrationElement(ip.position());
684 auto j = param.j(ig.intersection(),ip.position());
687 for (size_type i=0; i<lfsv_s.size(); i++)
688 r_s.accumulate(lfsv_s,i,j * psi_s[i] * factor);
695 for (size_type i=0; i<lfsu_s.size(); i++)
696 u_s += x_s(lfsu_s,i)*phi_s[i];
699 auto b = param.b(cell_inside,iplocal_s);
700 auto normalflux = b*n_F_local;
704 if (normalflux<-1
e-30)
705 DUNE_THROW(Dune::Exception,
706 "Outflow boundary condition on inflow! [b(" 707 << geo.global(ip.position()) <<
") = " 711 auto term1 = u_s * normalflux *factor;
712 for (size_type i=0; i<lfsv_s.size(); i++)
713 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
716 auto o = param.o(ig.intersection(),ip.position());
719 for (size_type i=0; i<lfsv_s.size(); i++)
720 r_s.accumulate(lfsv_s,i,o * psi_s[i] * factor);
727 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
728 auto& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
731 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
732 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
733 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
737 for (size_type i=0; i<lfsu_s.size(); i++)
738 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
741 auto g = param.g(cell_inside,iplocal_s);
744 RF omegaup_s, omegaup_n;
757 auto term1 = (omegaup_s*u_s + omegaup_n*g) * normalflux *factor;
758 for (size_type i=0; i<lfsv_s.size(); i++)
759 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
762 auto term2 = (An_F_s*gradu_s) * factor;
763 for (size_type i=0; i<lfsv_s.size(); i++)
764 r_s.accumulate(lfsv_s,i,-term2 * psi_s[i]);
767 auto term3 = (u_s-g) * factor;
768 for (size_type i=0; i<lfsv_s.size(); i++)
769 r_s.accumulate(lfsv_s,i,term3 * theta * (An_F_s*tgradpsi_s[i]));
772 auto term4 = penalty_factor * (u_s-g) * factor;
773 for (size_type i=0; i<lfsv_s.size(); i++)
774 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
778 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
780 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
784 using RF =
typename LFSV::Traits::FiniteElementType::
785 Traits::LocalBasisType::Traits::RangeFieldType;
786 using size_type =
typename LFSV::Traits::SizeType;
789 const int dim = IG::dimension;
790 const int order = std::max(
791 lfsu_s.finiteElement().localBasis().order(),
792 lfsv_s.finiteElement().localBasis().order()
796 const auto& cell_inside = ig.inside();
799 auto geo = ig.geometry();
800 auto geo_inside = cell_inside.geometry();
803 auto geo_in_inside = ig.geometryInInside();
807 auto local_inside = ref_el_inside.position(0,0);
808 auto A_s = param.A(cell_inside,local_inside);
811 auto h_F = geo_inside.volume()/geo.volume();
814 auto n_F = ig.centerUnitOuterNormal();
815 Dune::FieldVector<RF,dim> An_F_s;
819 harmonic_average = An_F_s*n_F;
821 harmonic_average = 1.0;
824 auto order_s = lfsu_s.finiteElement().localBasis().order();
825 auto degree = order_s;
828 auto penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
831 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
834 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
837 auto intorder = intorderadd+quadrature_factor*order;
840 auto bctype = param.bctype(ig.intersection(),ip.position());
847 auto iplocal_s = geo_in_inside.global(ip.position());
850 auto n_F_local = ig.unitOuterNormal(ip.position());
853 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
856 auto factor = ip.weight() * geo.integrationElement(ip.position());
859 auto b = param.b(cell_inside,iplocal_s);
860 auto normalflux = b*n_F_local;
864 if (normalflux<-1
e-30)
865 DUNE_THROW(Dune::Exception,
866 "Outflow boundary condition on inflow! [b(" 867 << geo.global(ip.position()) <<
") = " 868 << b <<
")" << n_F_local <<
" " << normalflux);
871 for (size_type j=0; j<lfsu_s.size(); j++)
872 for (size_type i=0; i<lfsu_s.size(); i++)
873 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * normalflux * factor * phi_s[i]);
879 auto& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
882 jac = geo_inside.jacobianInverseTransposed(iplocal_s);
883 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
886 RF omegaup_s, omegaup_n;
899 for (size_type j=0; j<lfsu_s.size(); j++)
900 for (size_type i=0; i<lfsu_s.size(); i++)
901 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux * factor * phi_s[i]);
904 for (size_type j=0; j<lfsu_s.size(); j++)
905 for (size_type i=0; i<lfsu_s.size(); i++)
906 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,-(An_F_s*tgradphi_s[j]) * factor * phi_s[i]);
909 for (size_type j=0; j<lfsu_s.size(); j++)
910 for (size_type i=0; i<lfsu_s.size(); i++)
911 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * (An_F_s*tgradphi_s[i]));
914 for (size_type j=0; j<lfsu_s.size(); j++)
915 for (size_type i=0; i<lfsu_s.size(); i++)
916 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,penalty_factor * phi_s[j] * phi_s[i] * factor);
921 template<
typename EG,
typename LFSV,
typename R>
925 using size_type =
typename LFSV::Traits::SizeType;
928 const auto& cell = eg.entity();
931 auto geo = eg.geometry();
934 auto order = lfsv.finiteElement().localBasis().order();
935 auto intorder = intorderadd + 2 * order;
939 auto& phi = cache[order].evaluateFunction(ip.position(),lfsv.finiteElement().localBasis());
942 auto f = param.f(cell,ip.position());
945 auto factor = ip.weight() * geo.integrationElement(ip.position());
946 for (size_type i=0; i<lfsv.size(); i++)
947 r.accumulate(lfsv,i,-f*phi[i]*factor);
964 int quadrature_factor;
967 using LocalBasisType =
typename FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType;
979 std::vector<Cache> cache;
982 void element_size (
const GEO& geo,
typename GEO::ctype& hmin,
typename GEO::ctype hmax)
const 984 using DF =
typename GEO::ctype;
987 const int dim = GEO::coorddimension;
990 Dune::FieldVector<DF,dim> x = geo.corner(0);
992 hmin = hmax = x.two_norm();
997 Dune::GeometryType gt = geo.type();
998 for (
int i=0; i<Dune::ReferenceElements<DF,dim>::general(gt).size(dim-1); i++)
1000 Dune::FieldVector<DF,dim> x = geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,0,dim));
1001 x -= geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,1,dim));
1002 hmin = std::min(hmin,x.two_norm());
1003 hmax = std::max(hmax,x.two_norm());
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104
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:600
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:429
sparsity pattern generator
Definition: pattern.hh:13
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton() ...
Definition: numericaljacobianapply.hh:180
Definition: convectiondiffusiondg.hh:30
const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
Definition: adaptivity.hh:27
void setTime(double t)
set time in parameter class
Definition: convectiondiffusiondg.hh:952
ConvectionDiffusionDG(T ¶m_, ConvectionDiffusionDGMethod::Type method_=ConvectionDiffusionDGMethod::NIPG, ConvectionDiffusionDGWeights::Type weights_=ConvectionDiffusionDGWeights::weightsOff, Real alpha_=0.0, int intorderadd_=0)
constructor: pass parameter object
Definition: convectiondiffusiondg.hh:81
Definition: convectiondiffusionparameter.hh:65
Definition: convectiondiffusiondg.hh:37
Type
Definition: convectiondiffusionparameter.hh:65
Definition: convectiondiffusiondg.hh:32
sparsity pattern generator
Definition: pattern.hh:29
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericaljacobianapply.hh:285
Default flags for all local operators.
Definition: flags.hh:18
const Entity & e
Definition: localfunctionspace.hh:111
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: convectiondiffusionparameter.hh:65
Definition: convectiondiffusionparameter.hh:65
Definition: convectiondiffusiondg.hh:32
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Definition: convectiondiffusionparameter.hh:65
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusiondg.hh:180
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:779
Type
Definition: convectiondiffusiondg.hh:32
Definition: convectiondiffusiondg.hh:55
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:17
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:246
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
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:100
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:922
Definition: convectiondiffusiondg.hh:35