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/common/static_assert.hh>
10 #include<dune/geometry/type.hh>
11 #include<dune/geometry/referenceelements.hh>
12 #include<dune/geometry/quadraturerules.hh>
39 template<
typename T,
typename FiniteElementMap>
58 : param(param_), intorderadd(intorderadd_)
63 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
64 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
67 typedef typename LFSU::Traits::FiniteElementType::
68 Traits::LocalBasisType::Traits::DomainFieldType DF;
69 typedef typename LFSU::Traits::FiniteElementType::
70 Traits::LocalBasisType::Traits::RangeFieldType RF;
71 typedef typename LFSU::Traits::FiniteElementType::
72 Traits::LocalBasisType::Traits::JacobianType JacobianType;
73 typedef typename LFSU::Traits::FiniteElementType::
74 Traits::LocalBasisType::Traits::RangeType RangeType;
76 typedef typename LFSU::Traits::SizeType size_type;
79 const int dim = EG::Geometry::dimension;
82 Dune::GeometryType gt = eg.geometry().type();
83 const int intorder = intorderadd+2*lfsu.finiteElement().localBasis().order();
84 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
87 typename T::Traits::PermTensorType tensor;
88 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
89 tensor = param.A(eg.entity(),localcenter);
92 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
97 const std::vector<RangeType>& phi = cache.evaluateFunction(it->position(),lfsu.finiteElement().localBasis());
101 for (size_type i=0; i<lfsu.size(); i++)
102 u += x(lfsu,i)*phi[i];
107 const std::vector<JacobianType>& js = cache.evaluateJacobian(it->position(),lfsu.finiteElement().localBasis());
110 const typename EG::Geometry::JacobianInverseTransposed jac =
111 eg.geometry().jacobianInverseTransposed(it->position());
112 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
113 for (size_type i=0; i<lfsu.size(); i++)
114 jac.mv(js[i][0],gradphi[i]);
117 Dune::FieldVector<RF,dim> gradu(0.0);
118 for (size_type i=0; i<lfsu.size(); i++)
119 gradu.axpy(x(lfsu,i),gradphi[i]);
122 Dune::FieldVector<RF,dim> Agradu(0.0);
123 tensor.umv(gradu,Agradu);
126 typename T::Traits::RangeType b = param.b(eg.entity(),it->position());
127 typename T::Traits::RangeFieldType c = param.c(eg.entity(),it->position());
128 typename T::Traits::RangeFieldType
f = param.f(eg.entity(),it->position());
131 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
132 for (size_type i=0; i<lfsu.size(); i++)
133 r.accumulate(lfsu,i,( Agradu*gradphi[i] - u*(b*gradphi[i]) + (c*u-f)*phi[i] )*factor);
138 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
143 typedef typename LFSU::Traits::FiniteElementType::
144 Traits::LocalBasisType::Traits::DomainFieldType DF;
145 typedef typename LFSU::Traits::FiniteElementType::
146 Traits::LocalBasisType::Traits::RangeFieldType RF;
147 typedef typename LFSU::Traits::FiniteElementType::
148 Traits::LocalBasisType::Traits::JacobianType JacobianType;
149 typedef typename LFSU::Traits::FiniteElementType::
150 Traits::LocalBasisType::Traits::RangeType RangeType;
151 typedef typename LFSU::Traits::SizeType size_type;
154 const int dim = EG::Geometry::dimension;
157 Dune::GeometryType gt = eg.geometry().type();
158 const int intorder = intorderadd+2*lfsu.finiteElement().localBasis().order();
159 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
162 typename T::Traits::PermTensorType tensor;
163 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
164 tensor = param.A(eg.entity(),localcenter);
167 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
172 const std::vector<JacobianType>& js = cache.evaluateJacobian(it->position(),lfsu.finiteElement().localBasis());
175 const typename EG::Geometry::JacobianInverseTransposed jac
176 = eg.geometry().jacobianInverseTransposed(it->position());
177 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
178 std::vector<Dune::FieldVector<RF,dim> > Agradphi(lfsu.size());
179 for (size_type i=0; i<lfsu.size(); i++)
181 jac.mv(js[i][0],gradphi[i]);
182 tensor.mv(gradphi[i],Agradphi[i]);
188 const std::vector<RangeType>& phi = cache.evaluateFunction(it->position(),lfsu.finiteElement().localBasis());
191 typename T::Traits::RangeType b = param.b(eg.entity(),it->position());
192 typename T::Traits::RangeFieldType c = param.c(eg.entity(),it->position());
195 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
196 for (size_type j=0; j<lfsu.size(); j++)
197 for (size_type i=0; i<lfsu.size(); i++)
198 mat.accumulate(lfsu,i,lfsu,j,( Agradphi[j]*gradphi[i]-phi[j]*(b*gradphi[i])+c*phi[j]*phi[i] )*factor);
203 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
205 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
209 typedef typename LFSV::Traits::FiniteElementType::
210 Traits::LocalBasisType::Traits::DomainFieldType DF;
211 typedef typename LFSV::Traits::FiniteElementType::
212 Traits::LocalBasisType::Traits::RangeFieldType RF;
213 typedef typename LFSV::Traits::FiniteElementType::
214 Traits::LocalBasisType::Traits::RangeType RangeType;
216 typedef typename LFSV::Traits::SizeType size_type;
219 const int dim = IG::dimension;
222 Dune::GeometryType gtface = ig.geometryInInside().type();
223 Dune::FieldVector<DF,dim-1> facecenterlocal = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
225 bctype = param.bctype(ig.intersection(),facecenterlocal);
231 const int intorder = intorderadd+2*lfsu_s.finiteElement().localBasis().order();
232 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
235 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
238 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
243 const std::vector<RangeType>& phi = cache.evaluateFunction(local,lfsu_s.finiteElement().localBasis());
248 typename T::Traits::RangeFieldType j = param.j(ig.intersection(),it->position());
251 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
252 for (size_type i=0; i<lfsu_s.size(); i++)
253 r_s.accumulate(lfsu_s,i,j*phi[i]*factor);
260 for (size_type i=0; i<lfsu_s.size(); i++)
261 u += x_s(lfsu_s,i)*phi[i];
264 typename T::Traits::RangeType b = param.b(*(ig.inside()),local);
265 const Dune::FieldVector<DF,dim> n = ig.unitOuterNormal(it->position());
268 typename T::Traits::RangeFieldType o = param.o(ig.intersection(),it->position());
271 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
272 for (size_type i=0; i<lfsu_s.size(); i++)
273 r_s.accumulate(lfsu_s,i,( (b*n)*u + o)*phi[i]*factor);
279 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
281 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
285 typedef typename LFSV::Traits::FiniteElementType::
286 Traits::LocalBasisType::Traits::DomainFieldType DF;
287 typedef typename LFSV::Traits::FiniteElementType::
288 Traits::LocalBasisType::Traits::RangeFieldType RF;
289 typedef typename LFSV::Traits::FiniteElementType::
290 Traits::LocalBasisType::Traits::RangeType RangeType;
292 typedef typename LFSV::Traits::SizeType size_type;
295 const int dim = IG::dimension;
298 Dune::GeometryType gtface = ig.geometryInInside().type();
299 Dune::FieldVector<DF,dim-1> facecenterlocal = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
301 bctype = param.bctype(ig.intersection(),facecenterlocal);
308 const int intorder = intorderadd+2*lfsu_s.finiteElement().localBasis().order();
309 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
312 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
315 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
320 const std::vector<RangeType>& phi = cache.evaluateFunction(local,lfsu_s.finiteElement().localBasis());
323 typename T::Traits::RangeType b = param.b(*(ig.inside()),local);
324 const Dune::FieldVector<DF,dim> n = ig.unitOuterNormal(it->position());
327 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
328 for (size_type j=0; j<lfsu_s.size(); j++)
329 for (size_type i=0; i<lfsu_s.size(); i++)
330 mat_s.accumulate(lfsu_s,i,lfsu_s,j,(b*n)*phi[j]*phi[i]*factor);
344 typedef typename FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
369 enum { dim = T::Traits::GridViewType::dimension };
371 typedef typename T::Traits::RangeFieldType Real;
390 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
391 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
394 typedef typename LFSU::Traits::FiniteElementType::
395 Traits::LocalBasisType::Traits::DomainFieldType DF;
396 typedef typename LFSU::Traits::FiniteElementType::
397 Traits::LocalBasisType::Traits::RangeFieldType RF;
398 typedef typename LFSU::Traits::FiniteElementType::
399 Traits::LocalBasisType::Traits::RangeType RangeType;
400 typedef typename LFSU::Traits::SizeType size_type;
403 const int dim = EG::Geometry::dimension;
404 const int intorder = 2*lfsu.finiteElement().localBasis().order();
407 Dune::GeometryType gt = eg.geometry().type();
408 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
412 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
415 std::vector<RangeType> phi(lfsu.size());
416 lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
420 for (size_type i=0; i<lfsu.size(); i++)
421 u += x(lfsu,i)*phi[i];
424 typename T::Traits::RangeFieldType c = param.c(eg.entity(),it->position());
427 typename T::Traits::RangeFieldType
f = param.f(eg.entity(),it->position());
430 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
431 sum += (f*f-c*c*u*u)*factor;
435 DF h_T = diameter(eg.geometry());
436 r.accumulate(lfsv,0,h_T*h_T*sum);
442 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
444 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
445 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
446 R& r_s, R& r_n)
const
449 typedef typename LFSU::Traits::FiniteElementType::
450 Traits::LocalBasisType::Traits::DomainFieldType DF;
451 typedef typename LFSU::Traits::FiniteElementType::
452 Traits::LocalBasisType::Traits::RangeFieldType RF;
453 typedef typename LFSU::Traits::FiniteElementType::
454 Traits::LocalBasisType::Traits::JacobianType JacobianType;
455 typedef typename LFSU::Traits::SizeType size_type;
458 const int dim = IG::dimension;
461 const Dune::FieldVector<DF,dim>&
462 inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
463 const Dune::FieldVector<DF,dim>&
464 outside_local = Dune::ReferenceElements<DF,dim>::general(ig.outside()->type()).position(0,0);
465 typename T::Traits::PermTensorType A_s, A_n;
466 A_s = param.A(*(ig.inside()),inside_local);
467 A_n = param.A(*(ig.outside()),outside_local);
470 const int intorder = 2*lfsu_s.finiteElement().localBasis().order();
471 Dune::GeometryType gtface = ig.geometryInInside().type();
472 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
475 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
478 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
479 Dune::FieldVector<RF,dim> An_F_s;
481 Dune::FieldVector<RF,dim> An_F_n;
486 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
489 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
490 Dune::FieldVector<DF,dim> iplocal_n = ig.geometryInOutside().global(it->position());
493 std::vector<JacobianType> gradphi_s(lfsu_s.size());
494 lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
495 std::vector<JacobianType> gradphi_n(lfsu_n.size());
496 lfsu_n.finiteElement().localBasis().evaluateJacobian(iplocal_n,gradphi_n);
499 jac = ig.inside()->geometry().jacobianInverseTransposed(iplocal_s);
500 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
501 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
502 jac = ig.outside()->geometry().jacobianInverseTransposed(iplocal_n);
503 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
504 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
507 Dune::FieldVector<RF,dim> gradu_s(0.0);
508 for (size_type i=0; i<lfsu_s.size(); i++)
509 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
510 Dune::FieldVector<RF,dim> gradu_n(0.0);
511 for (size_type i=0; i<lfsu_n.size(); i++)
512 gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
515 RF factor = it->weight() * ig.geometry().integrationElement(it->position());
516 RF jump = (An_F_s*gradu_s)-(An_F_n*gradu_n);
517 sum += 0.25*jump*jump*factor;
522 DF h_T = std::max(diameter(ig.inside()->geometry()),diameter(ig.outside()->geometry()));
523 r_s.accumulate(lfsv_s,0,h_T*sum);
524 r_n.accumulate(lfsv_n,0,h_T*sum);
530 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
532 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
536 typedef typename LFSU::Traits::FiniteElementType::
537 Traits::LocalBasisType::Traits::DomainFieldType DF;
538 typedef typename LFSU::Traits::FiniteElementType::
539 Traits::LocalBasisType::Traits::RangeFieldType RF;
540 typedef typename LFSU::Traits::FiniteElementType::
541 Traits::LocalBasisType::Traits::JacobianType JacobianType;
542 typedef typename LFSU::Traits::SizeType size_type;
545 const int dim = IG::dimension;
548 const Dune::FieldVector<DF,dim>&
549 inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
550 typename T::Traits::PermTensorType A_s;
551 A_s = param.A(*(ig.inside()),inside_local);
552 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
553 Dune::FieldVector<RF,dim> An_F_s;
557 const int intorder = 2*lfsu_s.finiteElement().localBasis().order();
558 Dune::GeometryType gtface = ig.geometryInInside().type();
559 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
562 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
565 const Dune::FieldVector<DF,dim-1>
566 face_local = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
567 BCType bctype = param.bctype(ig.intersection(),face_local);
573 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
576 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
579 std::vector<JacobianType> gradphi_s(lfsu_s.size());
580 lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
583 jac = ig.inside()->geometry().jacobianInverseTransposed(iplocal_s);
584 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
585 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
588 Dune::FieldVector<RF,dim> gradu_s(0.0);
589 for (size_type i=0; i<lfsu_s.size(); i++)
590 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
593 RF j = param.j(ig.intersection(),it->position());
596 RF factor = it->weight() * ig.geometry().integrationElement(it->position());
597 RF jump = j+(An_F_s*gradu_s);
598 sum += jump*jump*factor;
603 DF h_T = diameter(ig.inside()->geometry());
604 r_s.accumulate(lfsv_s,0,h_T*sum);
611 typename GEO::ctype diameter (
const GEO& geo)
const
613 typedef typename GEO::ctype DF;
615 const int dim = GEO::coorddimension;
616 for (
int i=0; i<geo.corners(); i++)
618 Dune::FieldVector<DF,dim> xi = geo.corner(i);
619 for (
int j=i+1; j<geo.corners(); j++)
621 Dune::FieldVector<DF,dim> xj = geo.corner(j);
623 hmax = std::max(hmax,xj.two_norm());
652 enum { dim = T::Traits::GridViewType::dimension };
654 typedef typename T::Traits::RangeFieldType Real;
668 : param(param_), time(time_), dt(dt_), cmax(0)
672 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
673 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
676 typedef typename LFSU::Traits::FiniteElementType::
677 Traits::LocalBasisType::Traits::DomainFieldType DF;
678 typedef typename LFSU::Traits::FiniteElementType::
679 Traits::LocalBasisType::Traits::RangeFieldType RF;
680 typedef typename LFSU::Traits::FiniteElementType::
681 Traits::LocalBasisType::Traits::RangeType RangeType;
682 typedef typename LFSU::Traits::SizeType size_type;
685 const int dim = EG::Geometry::dimension;
686 const int intorder = 2*lfsu.finiteElement().localBasis().order();
689 Dune::GeometryType gt = eg.geometry().type();
690 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
697 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
700 std::vector<RangeType> phi(lfsu.size());
701 lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
705 for (size_type i=0; i<lfsu.size(); i++)
706 u += x(lfsu,i)*phi[i];
709 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
714 typename T::Traits::RangeFieldType f_down = param.f(eg.entity(),it->position());
715 param.setTime(time+0.5*dt);
716 typename T::Traits::RangeFieldType f_mid = param.f(eg.entity(),it->position());
717 param.setTime(time+dt);
718 typename T::Traits::RangeFieldType f_up = param.f(eg.entity(),it->position());
719 RF f_average = (1.0/6.0)*f_down + (2.0/3.0)*f_mid + (1.0/6.0)*f_up;
722 fsum_down += (f_down-f_average)*(f_down-f_average)*factor;
723 fsum_mid += (f_mid-f_average)*(f_mid-f_average)*factor;
724 fsum_up += (f_up-f_average)*(f_up-f_average)*factor;
728 DF h_T = diameter(eg.geometry());
729 r.accumulate(lfsv,0,(h_T*h_T/dt)*sum);
730 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) );
750 typename GEO::ctype diameter (
const GEO& geo)
const
752 typedef typename GEO::ctype DF;
754 const int dim = GEO::coorddimension;
755 for (
int i=0; i<geo.corners(); i++)
757 Dune::FieldVector<DF,dim> xi = geo.corner(i);
758 for (
int j=i+1; j<geo.corners(); j++)
760 Dune::FieldVector<DF,dim> xj = geo.corner(j);
762 hmax = std::max(hmax,xj.two_norm());
771 template<
typename T,
typename EG>
778 template<
typename X,
typename Y>
781 y[0] = t.f(eg.entity(),x);
809 enum { dim = T::Traits::GridViewType::dimension };
811 typedef typename T::Traits::RangeFieldType Real;
826 : param(param_), time(time_), dt(dt_)
830 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
831 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
834 typedef typename LFSU::Traits::FiniteElementType::
835 Traits::LocalBasisType::Traits::DomainFieldType DF;
836 typedef typename LFSU::Traits::FiniteElementType::
837 Traits::LocalBasisType::Traits::RangeFieldType RF;
838 typedef typename LFSU::Traits::FiniteElementType::
839 Traits::LocalBasisType::Traits::RangeType RangeType;
840 typedef typename LFSU::Traits::SizeType size_type;
841 typedef typename LFSU::Traits::FiniteElementType::
842 Traits::LocalBasisType::Traits::JacobianType JacobianType;
845 const int dim = EG::Geometry::dimension;
846 const int intorder = 2*lfsu.finiteElement().localBasis().order();
849 Dune::GeometryType gt = eg.geometry().type();
850 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
854 std::vector<RF> f_up, f_down, f_mid;
856 lfsu.finiteElement().localInterpolation().interpolate(f_adapter,f_down);
857 param.setTime(time+0.5*dt);
858 lfsu.finiteElement().localInterpolation().interpolate(f_adapter,f_mid);
859 param.setTime(time+dt);
860 lfsu.finiteElement().localInterpolation().interpolate(f_adapter,f_up);
865 RF fsum_grad_up(0.0);
866 RF fsum_grad_mid(0.0);
867 RF fsum_grad_down(0.0);
868 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
871 std::vector<RangeType> phi(lfsu.size());
872 lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
876 for (size_type i=0; i<lfsu.size(); i++)
877 u += x(lfsu,i)*phi[i];
880 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
884 std::vector<JacobianType> js(lfsu.size());
885 lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
888 const typename EG::Geometry::JacobianInverseTransposed jac =
889 eg.geometry().jacobianInverseTransposed(it->position());
890 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
891 for (size_type i=0; i<lfsu.size(); i++)
892 jac.mv(js[i][0],gradphi[i]);
895 Dune::FieldVector<RF,dim> gradu(0.0);
896 for (size_type i=0; i<lfsu.size(); i++)
897 gradu.axpy(x(lfsu,i),gradphi[i]);
900 sum_grad += (gradu*gradu)*factor;
903 Dune::FieldVector<RF,dim> gradf_down(0.0);
904 for (size_type i=0; i<lfsu.size(); i++) gradf_down.axpy(f_down[i],gradphi[i]);
905 Dune::FieldVector<RF,dim> gradf_mid(0.0);
906 for (size_type i=0; i<lfsu.size(); i++) gradf_mid.axpy(f_mid[i],gradphi[i]);
907 Dune::FieldVector<RF,dim> gradf_up(0.0);
908 for (size_type i=0; i<lfsu.size(); i++) gradf_up.axpy(f_up[i],gradphi[i]);
909 Dune::FieldVector<RF,dim> gradf_average(0.0);
910 for (size_type i=0; i<lfsu.size(); i++)
911 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]);
914 gradf_down -= gradf_average;
915 fsum_grad_down += (gradf_down*gradf_down)*factor;
916 gradf_mid -= gradf_average;
917 fsum_grad_mid += (gradf_mid*gradf_mid)*factor;
918 gradf_up -= gradf_average;
919 fsum_grad_up += (gradf_up*gradf_up)*factor;
923 DF h_T = diameter(eg.geometry());
924 r.accumulate(lfsv,0,dt * sum_grad);
925 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));
930 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
932 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
936 typedef typename LFSU::Traits::FiniteElementType::
937 Traits::LocalBasisType::Traits::DomainFieldType DF;
938 typedef typename LFSU::Traits::FiniteElementType::
939 Traits::LocalBasisType::Traits::RangeFieldType RF;
942 const int dim = IG::dimension;
945 const int intorder = 2*lfsu_s.finiteElement().localBasis().order();
946 Dune::GeometryType gtface = ig.geometryInInside().type();
947 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
950 const Dune::FieldVector<DF,dim-1>
951 face_local = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
952 BCType bctype = param.bctype(ig.intersection(),face_local);
959 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
963 RF j_down = param.j(ig.intersection(),it->position());
964 param.setTime(time+0.5*dt);
965 RF j_mid = param.j(ig.intersection(),it->position());
966 param.setTime(time+dt);
967 RF j_up = param.j(ig.intersection(),it->position());
970 RF factor = it->weight() * ig.geometry().integrationElement(it->position());
971 sum_down += (j_down-j_mid)*(j_down-j_mid)*factor;
972 sum_up += (j_up-j_mid)*(j_up-j_mid)*factor;
977 DF h_T = diameter(ig.inside()->geometry());
978 r_s.accumulate(lfsv_s,0,(h_T+dt*dt/h_T)*dt*0.5*(sum_down+sum_up));
987 typename GEO::ctype diameter (
const GEO& geo)
const
989 typedef typename GEO::ctype DF;
991 const int dim = GEO::coorddimension;
992 for (
int i=0; i<geo.corners(); i++)
994 Dune::FieldVector<DF,dim> xi = geo.corner(i);
995 for (
int j=i+1; j<geo.corners(); j++)
997 Dune::FieldVector<DF,dim> xj = geo.corner(j);
999 hmax = std::max(hmax,xj.two_norm());
Definition: convectiondiffusionparameter.hh:68
Definition: convectiondiffusionfem.hh:51
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:443
ConvectionDiffusionTemporalResidualEstimator1(T ¶m_, double time_, double dt_)
constructor: pass parameter object
Definition: convectiondiffusionfem.hh:667
Definition: convectiondiffusionfem.hh:660
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusionfem.hh:139
Definition: convectiondiffusionfem.hh:380
Default flags for all local operators.
Definition: flags.hh:18
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:391
Definition: convectiondiffusionfem.hh:772
Definition: convectiondiffusionfem.hh:40
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:280
void evaluate(const X &x, Y &y) const
Definition: convectiondiffusionfem.hh:779
Definition: convectiondiffusionfem.hh:382
Definition: convectiondiffusionfem.hh:817
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:64
Definition: convectiondiffusionfem.hh:663
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
static const int dim
Definition: adaptivity.hh:82
Definition: convectiondiffusionfem.hh:376
Definition: convectiondiffusionfem.hh:648
void clearCmax()
Definition: convectiondiffusionfem.hh:733
Definition: convectiondiffusionfem.hh:820
ConvectionDiffusionFEMResidualEstimator(T ¶m_)
constructor: pass parameter object
Definition: convectiondiffusionfem.hh:385
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
const IG & ig
Definition: common/constraints.hh:146
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:204
void setTime(double t)
set time in parameter class
Definition: convectiondiffusionfem.hh:336
ConvectionDiffusionFEM(T ¶m_, int intorderadd_=0)
Definition: convectiondiffusionfem.hh:57
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:931
Definition: convectiondiffusionfem.hh:659
Definition: convectiondiffusionfem.hh:377
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:531
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:16
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
const F & f
Definition: common/constraints.hh:145
Definition: convectiondiffusionfem.hh:805
Definition: convectiondiffusionfem.hh:54
Definition: convectiondiffusionfem.hh:381
sparsity pattern generator
Definition: pattern.hh:14
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
const EG & eg
Definition: common/constraints.hh:277
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:673
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:831
Definition: convectiondiffusionfem.hh:816
Definition: convectiondiffusionfem.hh:55
Definition: convectiondiffusionparameter.hh:68
Definition: convectiondiffusionfem.hh:821
ConvectionDiffusionTemporalResidualEstimator(T ¶m_, double time_, double dt_)
constructor: pass parameter object
Definition: convectiondiffusionfem.hh:825
CD_RHS_LocalAdapter(const T &t_, const EG &eg_)
Definition: convectiondiffusionfem.hh:775
Type
Definition: convectiondiffusionparameter.hh:68
Definition: convectiondiffusionparameter.hh:68
double getCmax() const
Definition: convectiondiffusionfem.hh:738
Definition: convectiondiffusionfem.hh:366