dune-localfunctions  2.2.0
brezzidouglasmarini12dlocalinterpolation.hh
Go to the documentation of this file.
00001 #ifndef DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI12DLOCALINTERPOLATION_HH
00002 #define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI12DLOCALINTERPOLATION_HH
00003 
00004 #include <vector>
00005 
00006 #include <dune/geometry/quadraturerules.hh>
00007 
00008 namespace Dune 
00009 {
00018   template<class LB>
00019   class BDM12DLocalInterpolation 
00020   {
00021 
00022 public:
00024     BDM12DLocalInterpolation ()
00025     {
00026         sign0 = sign1 = sign2 = 1.0;
00027     }
00028 
00034     BDM12DLocalInterpolation (unsigned int s)
00035     {
00036       sign0 = sign1 = sign2 = 1.0;
00037       if (s & 1)
00038       {
00039         sign0 = -1.0;
00040       }
00041       if (s & 2)
00042       {
00043         sign1 = -1.0;
00044       }
00045       if (s & 4)
00046       {
00047         sign2 = -1.0;
00048       }
00049       
00050       n0[0] =  0.0;
00051       n0[1] = -1.0;
00052       n1[0] = -1.0;
00053       n1[1] =  0.0;
00054       n2[0] =  1.0/sqrt(2.0);
00055       n2[1] =  1.0/sqrt(2.0);
00056       c0 =  0.5*n0[0] - 1.0*n0[1];
00057       c1 = -1.0*n1[0] + 0.5*n1[1];
00058       c2 =  0.5*n2[0] + 0.5*n2[1];
00059     }
00060 
00069     template<typename F, typename C>
00070     void interpolate (const F& f, std::vector<C>& out) const
00071     {
00072       // f gives v*outer normal at a point on the edge!
00073       typedef typename LB::Traits::RangeFieldType Scalar;
00074       typename F::Traits::RangeType y;
00075 
00076       out.resize(6);
00077       fill(out.begin(), out.end(), 0.0);
00078 
00079       const int qOrder = 4;
00080       const Dune::QuadratureRule<Scalar,1>& rule = Dune::QuadratureRules<Scalar,1>::rule(Dune::GeometryType(Dune::GeometryType::simplex,1), qOrder);
00081 
00082       for (typename Dune::QuadratureRule<Scalar,1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
00083       {
00084         Scalar qPos = it->position();
00085         typename LB::Traits::DomainType localPos;
00086         
00087         localPos[0] = qPos;
00088         localPos[1] = 0.0;
00089         f.evaluate(localPos, y); 
00090         out[0] += (y[0]*n0[0] + y[1]*n0[1])*it->weight()*sign0/c0;    
00091         out[3] += (y[0]*n0[0] + y[1]*n0[1])*(2.0*qPos - 1.0)*it->weight()/c0;    
00092 
00093         localPos[0] = 0.0;
00094         localPos[1] = qPos;
00095         f.evaluate(localPos, y); 
00096         out[1] += (y[0]*n1[0] + y[1]*n1[1])*it->weight()*sign1/c1;    
00097         out[4] += (y[0]*n1[0] + y[1]*n1[1])*(1.0 - 2.0*qPos)*it->weight()/c1;    
00098 
00099         localPos[0] = 1.0 - qPos;
00100         localPos[1] = qPos;
00101         f.evaluate(localPos, y); 
00102         out[2] += (y[0]*n2[0] + y[1]*n2[1])*it->weight()*sign2/c2;    
00103         out[5] += (y[0]*n2[0] + y[1]*n2[1])*(2.0*qPos - 1.0)*it->weight()/c2;    
00104       }
00105     }
00106 
00107 private:
00108     typename LB::Traits::RangeFieldType sign0,sign1,sign2;
00109     typename LB::Traits::DomainType n0,n1,n2;
00110     typename LB::Traits::RangeFieldType c0,c1,c2;
00111   };
00112 }
00113 
00114 #endif // DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI12DLOCALINTERPOLATION_HH