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