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