2 #ifndef DUNE_PDELAB_CONVECTIONDIFFUSIONFEM_HH
3 #define DUNE_PDELAB_CONVECTIONDIFFUSIONFEM_HH
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/fvector.hh>
9 #include<dune/geometry/type.hh>
10 #include<dune/geometry/referenceelements.hh>
11 #include<dune/geometry/quadraturerules.hh>
38 template<
typename T,
typename FiniteElementMap>
57 : param(param_), intorderadd(intorderadd_)
62 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
63 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
66 typedef typename LFSU::Traits::FiniteElementType::
67 Traits::LocalBasisType::Traits::DomainFieldType DF;
68 typedef typename LFSU::Traits::FiniteElementType::
69 Traits::LocalBasisType::Traits::RangeFieldType RF;
70 typedef typename LFSU::Traits::FiniteElementType::
71 Traits::LocalBasisType::Traits::JacobianType JacobianType;
72 typedef typename LFSU::Traits::FiniteElementType::
73 Traits::LocalBasisType::Traits::RangeType RangeType;
75 typedef typename LFSU::Traits::SizeType size_type;
78 const int dim = EG::Entity::dimension;
81 Dune::GeometryType gt = eg.geometry().type();
82 const int intorder = intorderadd+2*lfsu.finiteElement().localBasis().order();
83 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
86 typename T::Traits::PermTensorType tensor;
87 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
88 tensor = param.A(eg.entity(),localcenter);
91 for (
const auto& ip : rule)
96 const std::vector<RangeType>& phi = cache.evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
100 for (size_type i=0; i<lfsu.size(); i++)
101 u += x(lfsu,i)*phi[i];
106 const std::vector<JacobianType>& js = cache.evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
109 const typename EG::Geometry::JacobianInverseTransposed jac =
110 eg.geometry().jacobianInverseTransposed(ip.position());
111 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
112 for (size_type i=0; i<lfsu.size(); i++)
113 jac.mv(js[i][0],gradphi[i]);
116 Dune::FieldVector<RF,dim> gradu(0.0);
117 for (size_type i=0; i<lfsu.size(); i++)
118 gradu.axpy(x(lfsu,i),gradphi[i]);
121 Dune::FieldVector<RF,dim> Agradu(0.0);
122 tensor.umv(gradu,Agradu);
125 typename T::Traits::RangeType b = param.b(eg.entity(),ip.position());
126 typename T::Traits::RangeFieldType c = param.c(eg.entity(),ip.position());
127 typename T::Traits::RangeFieldType f = param.f(eg.entity(),ip.position());
130 RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
131 for (size_type i=0; i<lfsu.size(); i++)
132 r.accumulate(lfsu,i,( Agradu*gradphi[i] - u*(b*gradphi[i]) + (c*u-f)*phi[i] )*factor);
137 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
142 typedef typename LFSU::Traits::FiniteElementType::
143 Traits::LocalBasisType::Traits::DomainFieldType DF;
144 typedef typename LFSU::Traits::FiniteElementType::
145 Traits::LocalBasisType::Traits::RangeFieldType RF;
146 typedef typename LFSU::Traits::FiniteElementType::
147 Traits::LocalBasisType::Traits::JacobianType JacobianType;
148 typedef typename LFSU::Traits::FiniteElementType::
149 Traits::LocalBasisType::Traits::RangeType RangeType;
150 typedef typename LFSU::Traits::SizeType size_type;
153 const int dim = EG::Entity::dimension;
156 Dune::GeometryType gt = eg.geometry().type();
157 const int intorder = intorderadd+2*lfsu.finiteElement().localBasis().order();
158 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
161 typename T::Traits::PermTensorType tensor;
162 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
163 tensor = param.A(eg.entity(),localcenter);
166 for (
const auto& ip : rule)
171 const std::vector<JacobianType>& js = cache.evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
174 const typename EG::Geometry::JacobianInverseTransposed jac
175 = eg.geometry().jacobianInverseTransposed(ip.position());
176 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
177 std::vector<Dune::FieldVector<RF,dim> > Agradphi(lfsu.size());
178 for (size_type i=0; i<lfsu.size(); i++)
180 jac.mv(js[i][0],gradphi[i]);
181 tensor.mv(gradphi[i],Agradphi[i]);
187 const std::vector<RangeType>& phi = cache.evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
190 typename T::Traits::RangeType b = param.b(eg.entity(),ip.position());
191 typename T::Traits::RangeFieldType c = param.c(eg.entity(),ip.position());
194 RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
195 for (size_type j=0; j<lfsu.size(); j++)
196 for (size_type i=0; i<lfsu.size(); i++)
197 mat.accumulate(lfsu,i,lfsu,j,( Agradphi[j]*gradphi[i]-phi[j]*(b*gradphi[i])+c*phi[j]*phi[i] )*factor);
202 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
204 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
208 typedef typename LFSV::Traits::FiniteElementType::
209 Traits::LocalBasisType::Traits::DomainFieldType DF;
210 typedef typename LFSV::Traits::FiniteElementType::
211 Traits::LocalBasisType::Traits::RangeFieldType RF;
212 typedef typename LFSV::Traits::FiniteElementType::
213 Traits::LocalBasisType::Traits::RangeType RangeType;
215 typedef typename LFSV::Traits::SizeType size_type;
218 const int dim = IG::dimension;
221 auto inside_cell = ig.inside();
224 Dune::GeometryType gtface = ig.geometryInInside().type();
225 Dune::FieldVector<DF,dim-1> facecenterlocal = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
227 bctype = param.bctype(ig.intersection(),facecenterlocal);
233 const int intorder = intorderadd+2*lfsu_s.finiteElement().localBasis().order();
234 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
237 for (
const auto& ip : rule)
240 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
245 const std::vector<RangeType>& phi = cache.evaluateFunction(local,lfsu_s.finiteElement().localBasis());
250 typename T::Traits::RangeFieldType j = param.j(ig.intersection(),ip.position());
253 RF factor = ip.weight()*ig.geometry().integrationElement(ip.position());
254 for (size_type i=0; i<lfsu_s.size(); i++)
255 r_s.accumulate(lfsu_s,i,j*phi[i]*factor);
262 for (size_type i=0; i<lfsu_s.size(); i++)
263 u += x_s(lfsu_s,i)*phi[i];
266 typename T::Traits::RangeType b = param.b(inside_cell,local);
267 const Dune::FieldVector<DF,dim> n = ig.unitOuterNormal(ip.position());
270 typename T::Traits::RangeFieldType o = param.o(ig.intersection(),ip.position());
273 RF factor = ip.weight()*ig.geometry().integrationElement(ip.position());
274 for (size_type i=0; i<lfsu_s.size(); i++)
275 r_s.accumulate(lfsu_s,i,( (b*n)*u + o)*phi[i]*factor);
281 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
283 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
287 typedef typename LFSV::Traits::FiniteElementType::
288 Traits::LocalBasisType::Traits::DomainFieldType DF;
289 typedef typename LFSV::Traits::FiniteElementType::
290 Traits::LocalBasisType::Traits::RangeFieldType RF;
291 typedef typename LFSV::Traits::FiniteElementType::
292 Traits::LocalBasisType::Traits::RangeType RangeType;
294 typedef typename LFSV::Traits::SizeType size_type;
297 const int dim = IG::dimension;
300 auto inside_cell = ig.inside();
303 Dune::GeometryType gtface = ig.geometryInInside().type();
304 Dune::FieldVector<DF,dim-1> facecenterlocal = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
306 bctype = param.bctype(ig.intersection(),facecenterlocal);
313 const int intorder = intorderadd+2*lfsu_s.finiteElement().localBasis().order();
314 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
317 for (
const auto& ip : rule)
320 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
325 const std::vector<RangeType>& phi = cache.evaluateFunction(local,lfsu_s.finiteElement().localBasis());
328 typename T::Traits::RangeType b = param.b(inside_cell,local);
329 const Dune::FieldVector<DF,dim> n = ig.unitOuterNormal(ip.position());
332 RF factor = ip.weight()*ig.geometry().integrationElement(ip.position());
333 for (size_type j=0; j<lfsu_s.size(); j++)
334 for (size_type i=0; i<lfsu_s.size(); i++)
335 mat_s.accumulate(lfsu_s,i,lfsu_s,j,(b*n)*phi[j]*phi[i]*factor);
349 typedef typename FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
374 enum { dim = T::Traits::GridViewType::dimension };
376 typedef typename T::Traits::RangeFieldType Real;
395 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
396 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
399 typedef typename LFSU::Traits::FiniteElementType::
400 Traits::LocalBasisType::Traits::DomainFieldType DF;
401 typedef typename LFSU::Traits::FiniteElementType::
402 Traits::LocalBasisType::Traits::RangeFieldType RF;
403 typedef typename LFSU::Traits::FiniteElementType::
404 Traits::LocalBasisType::Traits::RangeType RangeType;
405 typedef typename LFSU::Traits::SizeType size_type;
408 const int dim = EG::Geometry::dimension;
409 const int intorder = 2*lfsu.finiteElement().localBasis().order();
412 Dune::GeometryType gt = eg.geometry().type();
413 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
417 for (
const auto& ip : rule)
420 std::vector<RangeType> phi(lfsu.size());
421 lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
425 for (size_type i=0; i<lfsu.size(); i++)
426 u += x(lfsu,i)*phi[i];
429 typename T::Traits::RangeFieldType c = param.c(eg.entity(),ip.position());
432 typename T::Traits::RangeFieldType f = param.f(eg.entity(),ip.position());
435 RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
436 sum += (f*f-c*c*u*u)*factor;
440 DF h_T = diameter(eg.geometry());
441 r.accumulate(lfsv,0,h_T*h_T*sum);
447 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
449 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
450 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
451 R& r_s, R& r_n)
const
454 typedef typename LFSU::Traits::FiniteElementType::
455 Traits::LocalBasisType::Traits::DomainFieldType DF;
456 typedef typename LFSU::Traits::FiniteElementType::
457 Traits::LocalBasisType::Traits::RangeFieldType RF;
458 typedef typename LFSU::Traits::FiniteElementType::
459 Traits::LocalBasisType::Traits::JacobianType JacobianType;
460 typedef typename LFSU::Traits::SizeType size_type;
463 const int dim = IG::dimension;
466 auto inside_cell = ig.inside();
467 auto outside_cell = ig.outside();
470 const Dune::FieldVector<DF,dim>&
471 inside_local = Dune::ReferenceElements<DF,dim>::general(inside_cell.type()).position(0,0);
472 const Dune::FieldVector<DF,dim>&
473 outside_local = Dune::ReferenceElements<DF,dim>::general(outside_cell.type()).position(0,0);
474 typename T::Traits::PermTensorType A_s, A_n;
475 A_s = param.A(inside_cell,inside_local);
476 A_n = param.A(outside_cell,outside_local);
479 const int intorder = 2*lfsu_s.finiteElement().localBasis().order();
480 Dune::GeometryType gtface = ig.geometryInInside().type();
481 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
484 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
487 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
488 Dune::FieldVector<RF,dim> An_F_s;
490 Dune::FieldVector<RF,dim> An_F_n;
495 for (
const auto& ip : rule)
498 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(ip.position());
499 Dune::FieldVector<DF,dim> iplocal_n = ig.geometryInOutside().global(ip.position());
502 std::vector<JacobianType> gradphi_s(lfsu_s.size());
503 lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
504 std::vector<JacobianType> gradphi_n(lfsu_n.size());
505 lfsu_n.finiteElement().localBasis().evaluateJacobian(iplocal_n,gradphi_n);
508 jac = inside_cell.geometry().jacobianInverseTransposed(iplocal_s);
509 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
510 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
511 jac = outside_cell.geometry().jacobianInverseTransposed(iplocal_n);
512 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
513 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
516 Dune::FieldVector<RF,dim> gradu_s(0.0);
517 for (size_type i=0; i<lfsu_s.size(); i++)
518 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
519 Dune::FieldVector<RF,dim> gradu_n(0.0);
520 for (size_type i=0; i<lfsu_n.size(); i++)
521 gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
524 RF factor = ip.weight() * ig.geometry().integrationElement(ip.position());
525 RF jump = (An_F_s*gradu_s)-(An_F_n*gradu_n);
526 sum += 0.25*jump*jump*factor;
531 DF h_T = std::max(diameter(inside_cell.geometry()),diameter(outside_cell.geometry()));
532 r_s.accumulate(lfsv_s,0,h_T*sum);
533 r_n.accumulate(lfsv_n,0,h_T*sum);
539 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
541 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
545 typedef typename LFSU::Traits::FiniteElementType::
546 Traits::LocalBasisType::Traits::DomainFieldType DF;
547 typedef typename LFSU::Traits::FiniteElementType::
548 Traits::LocalBasisType::Traits::RangeFieldType RF;
549 typedef typename LFSU::Traits::FiniteElementType::
550 Traits::LocalBasisType::Traits::JacobianType JacobianType;
551 typedef typename LFSU::Traits::SizeType size_type;
554 const int dim = IG::dimension;
556 auto inside_cell = ig.inside();
558 const Dune::FieldVector<DF,dim>&
559 inside_local = Dune::ReferenceElements<DF,dim>::general(inside_cell.type()).position(0,0);
560 typename T::Traits::PermTensorType A_s;
561 A_s = param.A(inside_cell,inside_local);
562 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
563 Dune::FieldVector<RF,dim> An_F_s;
567 const int intorder = 2*lfsu_s.finiteElement().localBasis().order();
568 Dune::GeometryType gtface = ig.geometryInInside().type();
569 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
572 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
575 const Dune::FieldVector<DF,dim-1>
576 face_local = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
577 BCType bctype = param.bctype(ig.intersection(),face_local);
583 for (
const auto& ip : rule)
586 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(ip.position());
589 std::vector<JacobianType> gradphi_s(lfsu_s.size());
590 lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
593 jac = inside_cell.geometry().jacobianInverseTransposed(iplocal_s);
594 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
595 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
598 Dune::FieldVector<RF,dim> gradu_s(0.0);
599 for (size_type i=0; i<lfsu_s.size(); i++)
600 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
603 RF j = param.j(ig.intersection(),ip.position());
606 RF factor = ip.weight() * ig.geometry().integrationElement(ip.position());
607 RF jump = j+(An_F_s*gradu_s);
608 sum += jump*jump*factor;
613 DF h_T = diameter(inside_cell.geometry());
614 r_s.accumulate(lfsv_s,0,h_T*sum);
621 typename GEO::ctype diameter (
const GEO& geo)
const
623 typedef typename GEO::ctype DF;
625 const int dim = GEO::coorddimension;
626 for (
int i=0; i<geo.corners(); i++)
628 Dune::FieldVector<DF,dim> xi = geo.corner(i);
629 for (
int j=i+1; j<geo.corners(); j++)
631 Dune::FieldVector<DF,dim> xj = geo.corner(j);
633 hmax = std::max(hmax,xj.two_norm());
662 enum { dim = T::Traits::GridViewType::dimension };
664 typedef typename T::Traits::RangeFieldType Real;
678 : param(param_), time(time_), dt(dt_), cmax(0)
682 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
683 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
686 typedef typename LFSU::Traits::FiniteElementType::
687 Traits::LocalBasisType::Traits::DomainFieldType DF;
688 typedef typename LFSU::Traits::FiniteElementType::
689 Traits::LocalBasisType::Traits::RangeFieldType RF;
690 typedef typename LFSU::Traits::FiniteElementType::
691 Traits::LocalBasisType::Traits::RangeType RangeType;
692 typedef typename LFSU::Traits::SizeType size_type;
695 const int dim = EG::Geometry::dimension;
696 const int intorder = 2*lfsu.finiteElement().localBasis().order();
699 Dune::GeometryType gt = eg.geometry().type();
700 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
707 for (
const auto& ip : rule)
710 std::vector<RangeType> phi(lfsu.size());
711 lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
715 for (size_type i=0; i<lfsu.size(); i++)
716 u += x(lfsu,i)*phi[i];
719 RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
724 typename T::Traits::RangeFieldType f_down = param.f(eg.entity(),ip.position());
725 param.setTime(time+0.5*dt);
726 typename T::Traits::RangeFieldType f_mid = param.f(eg.entity(),ip.position());
727 param.setTime(time+dt);
728 typename T::Traits::RangeFieldType f_up = param.f(eg.entity(),ip.position());
729 RF f_average = (1.0/6.0)*f_down + (2.0/3.0)*f_mid + (1.0/6.0)*f_up;
732 fsum_down += (f_down-f_average)*(f_down-f_average)*factor;
733 fsum_mid += (f_mid-f_average)*(f_mid-f_average)*factor;
734 fsum_up += (f_up-f_average)*(f_up-f_average)*factor;
738 DF h_T = diameter(eg.geometry());
739 r.accumulate(lfsv,0,(h_T*h_T/dt)*sum);
740 r.accumulate(lfsv,0,h_T*h_T * dt*((1.0/6.0)*fsum_down+(2.0/3.0)*fsum_mid+(1.0/6.0)*fsum_up) );
760 typename GEO::ctype diameter (
const GEO& geo)
const
762 typedef typename GEO::ctype DF;
764 const int dim = GEO::coorddimension;
765 for (
int i=0; i<geo.corners(); i++)
767 Dune::FieldVector<DF,dim> xi = geo.corner(i);
768 for (
int j=i+1; j<geo.corners(); j++)
770 Dune::FieldVector<DF,dim> xj = geo.corner(j);
772 hmax = std::max(hmax,xj.two_norm());
781 template<
typename T,
typename EG>
788 template<
typename X,
typename Y>
791 y[0] = t.f(eg.entity(),x);
819 enum { dim = T::Traits::GridViewType::dimension };
821 typedef typename T::Traits::RangeFieldType Real;
836 : param(param_), time(time_), dt(dt_)
840 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
841 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
844 typedef typename LFSU::Traits::FiniteElementType::
845 Traits::LocalBasisType::Traits::DomainFieldType DF;
846 typedef typename LFSU::Traits::FiniteElementType::
847 Traits::LocalBasisType::Traits::RangeFieldType RF;
848 typedef typename LFSU::Traits::FiniteElementType::
849 Traits::LocalBasisType::Traits::RangeType RangeType;
850 typedef typename LFSU::Traits::SizeType size_type;
851 typedef typename LFSU::Traits::FiniteElementType::
852 Traits::LocalBasisType::Traits::JacobianType JacobianType;
855 const int dim = EG::Geometry::dimension;
856 const int intorder = 2*lfsu.finiteElement().localBasis().order();
859 Dune::GeometryType gt = eg.geometry().type();
860 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
864 std::vector<RF> f_up, f_down, f_mid;
866 lfsu.finiteElement().localInterpolation().interpolate(f_adapter,f_down);
867 param.setTime(time+0.5*dt);
868 lfsu.finiteElement().localInterpolation().interpolate(f_adapter,f_mid);
869 param.setTime(time+dt);
870 lfsu.finiteElement().localInterpolation().interpolate(f_adapter,f_up);
875 RF fsum_grad_up(0.0);
876 RF fsum_grad_mid(0.0);
877 RF fsum_grad_down(0.0);
878 for (
const auto& ip : rule)
881 std::vector<RangeType> phi(lfsu.size());
882 lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
886 for (size_type i=0; i<lfsu.size(); i++)
887 u += x(lfsu,i)*phi[i];
890 RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
894 std::vector<JacobianType> js(lfsu.size());
895 lfsu.finiteElement().localBasis().evaluateJacobian(ip.position(),js);
898 const typename EG::Geometry::JacobianInverseTransposed jac =
899 eg.geometry().jacobianInverseTransposed(ip.position());
900 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
901 for (size_type i=0; i<lfsu.size(); i++)
902 jac.mv(js[i][0],gradphi[i]);
905 Dune::FieldVector<RF,dim> gradu(0.0);
906 for (size_type i=0; i<lfsu.size(); i++)
907 gradu.axpy(x(lfsu,i),gradphi[i]);
910 sum_grad += (gradu*gradu)*factor;
913 Dune::FieldVector<RF,dim> gradf_down(0.0);
914 for (size_type i=0; i<lfsu.size(); i++) gradf_down.axpy(f_down[i],gradphi[i]);
915 Dune::FieldVector<RF,dim> gradf_mid(0.0);
916 for (size_type i=0; i<lfsu.size(); i++) gradf_mid.axpy(f_mid[i],gradphi[i]);
917 Dune::FieldVector<RF,dim> gradf_up(0.0);
918 for (size_type i=0; i<lfsu.size(); i++) gradf_up.axpy(f_up[i],gradphi[i]);
919 Dune::FieldVector<RF,dim> gradf_average(0.0);
920 for (size_type i=0; i<lfsu.size(); i++)
921 gradf_average.axpy((1.0/6.0)*f_down[i]+(2.0/3.0)*f_mid[i]+(1.0/6.0)*f_up[i],gradphi[i]);
924 gradf_down -= gradf_average;
925 fsum_grad_down += (gradf_down*gradf_down)*factor;
926 gradf_mid -= gradf_average;
927 fsum_grad_mid += (gradf_mid*gradf_mid)*factor;
928 gradf_up -= gradf_average;
929 fsum_grad_up += (gradf_up*gradf_up)*factor;
933 DF h_T = diameter(eg.geometry());
934 r.accumulate(lfsv,0,dt * sum_grad);
935 r.accumulate(lfsv,0,dt*dt * dt*((1.0/6.0)*fsum_grad_down+(2.0/3.0)*fsum_grad_mid+(1.0/6.0)*fsum_grad_up));
940 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
942 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
946 typedef typename LFSU::Traits::FiniteElementType::
947 Traits::LocalBasisType::Traits::DomainFieldType DF;
948 typedef typename LFSU::Traits::FiniteElementType::
949 Traits::LocalBasisType::Traits::RangeFieldType RF;
952 const int dim = IG::dimension;
955 const int intorder = 2*lfsu_s.finiteElement().localBasis().order();
956 Dune::GeometryType gtface = ig.geometryInInside().type();
957 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
960 const Dune::FieldVector<DF,dim-1>
961 face_local = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
962 BCType bctype = param.bctype(ig.intersection(),face_local);
969 for (
const auto& ip : rule)
973 RF j_down = param.j(ig.intersection(),ip.position());
974 param.setTime(time+0.5*dt);
975 RF j_mid = param.j(ig.intersection(),ip.position());
976 param.setTime(time+dt);
977 RF j_up = param.j(ig.intersection(),ip.position());
980 RF factor = ip.weight() * ig.geometry().integrationElement(ip.position());
981 sum_down += (j_down-j_mid)*(j_down-j_mid)*factor;
982 sum_up += (j_up-j_mid)*(j_up-j_mid)*factor;
987 auto inside_cell = ig.inside();
988 DF h_T = diameter(inside_cell.geometry());
989 r_s.accumulate(lfsv_s,0,(h_T+dt*dt/h_T)*dt*0.5*(sum_down+sum_up));
998 typename GEO::ctype diameter (
const GEO& geo)
const
1000 typedef typename GEO::ctype DF;
1002 const int dim = GEO::coorddimension;
1003 for (
int i=0; i<geo.corners(); i++)
1005 Dune::FieldVector<DF,dim> xi = geo.corner(i);
1006 for (
int j=i+1; j<geo.corners(); j++)
1008 Dune::FieldVector<DF,dim> xj = geo.corner(j);
1010 hmax = std::max(hmax,xj.two_norm());
ConvectionDiffusionFEMResidualEstimator(T ¶m_)
constructor: pass parameter object
Definition: convectiondiffusionfem.hh:390
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusionfem.hh:941
Definition: convectiondiffusionfem.hh:382
Definition: convectiondiffusionfem.hh:50
void setTime(double t)
set time in parameter class
Definition: convectiondiffusionfem.hh:341
Definition: convectiondiffusionfem.hh:386
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusionfem.hh:540
Definition: convectiondiffusionfem.hh:669
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusionfem.hh:203
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
Definition: convectiondiffusionfem.hh:53
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusionfem.hh:138
Definition: convectiondiffusionfem.hh:673
const IG & ig
Definition: constraints.hh:147
Type
Definition: convectiondiffusionparameter.hh:67
Definition: convectiondiffusionfem.hh:830
ConvectionDiffusionTemporalResidualEstimator1(T ¶m_, double time_, double dt_)
constructor: pass parameter object
Definition: convectiondiffusionfem.hh:677
Default flags for all local operators.
Definition: flags.hh:18
void clearCmax()
Definition: convectiondiffusionfem.hh:743
Definition: convectiondiffusionparameter.hh:67
Definition: convectiondiffusionfem.hh:658
sparsity pattern generator
Definition: pattern.hh:13
Definition: convectiondiffusionparameter.hh:67
ConvectionDiffusionTemporalResidualEstimator(T ¶m_, double time_, double dt_)
constructor: pass parameter object
Definition: convectiondiffusionfem.hh:835
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:63
Definition: convectiondiffusionfem.hh:827
ConvectionDiffusionFEM(T ¶m_, int intorderadd_=0)
Definition: convectiondiffusionfem.hh:56
Definition: convectiondiffusionparameter.hh:67
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: convectiondiffusionfem.hh:448
Definition: convectiondiffusionfem.hh:831
Definition: convectiondiffusionfem.hh:39
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:396
Definition: convectiondiffusionfem.hh:385
void jacobian_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, M &mat_s) const
Definition: convectiondiffusionfem.hh:282
Definition: adaptivity.hh:27
Definition: convectiondiffusionfem.hh:815
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
Definition: convectiondiffusionfem.hh:387
double getCmax() const
Definition: convectiondiffusionfem.hh:748
static const int dim
Definition: adaptivity.hh:83
Definition: convectiondiffusionfem.hh:826
Definition: convectiondiffusionfem.hh:371
Definition: convectiondiffusionfem.hh:782
Definition: convectiondiffusionfem.hh:381
Definition: convectiondiffusionfem.hh:670
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:15
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:841
Definition: convectiondiffusionfem.hh:54
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:683
void evaluate(const X &x, Y &y) const
Definition: convectiondiffusionfem.hh:789
CD_RHS_LocalAdapter(const T &t_, const EG &eg_)
Definition: convectiondiffusionfem.hh:785
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538