dune-localfunctions  2.2.0
raviartthomas2q2dlocalinterpolation.hh
Go to the documentation of this file.
00001 #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS2Q2DLOCALINTERPOLATION_HH
00002 #define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS2Q2DLOCALINTERPOLATION_HH
00003 
00004 #include <vector>
00005 
00006 #include <dune/geometry/quadraturerules.hh>
00007 
00008 namespace Dune 
00009 {
00010 
00019   template<class LB>
00020   class RT2Q2DLocalInterpolation 
00021   {
00022   
00023 public:
00025     RT2Q2DLocalInterpolation ()
00026     {
00027       sign0 = sign1 = sign2 = sign3 = 1.0;
00028     }
00029 
00035     RT2Q2DLocalInterpolation (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(24);
00082       fill(out.begin(), out.end(), 0.0);
00083 
00084       const int qOrder = 6;
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(); it!=rule.end(); ++it)
00088       {
00089         Scalar qPos = it->position();
00090         typename LB::Traits::DomainType localPos;
00091         
00092         localPos[0] = 0.0;
00093         localPos[1] = qPos;
00094         f.evaluate(localPos, y); 
00095         out[0] += (y[0]*n0[0] + y[1]*n0[1])*it->weight()*sign0;    
00096         out[1] += (y[0]*n0[0] + y[1]*n0[1])*(2.0*qPos - 1.0)*it->weight();
00097         out[2] += (y[0]*n0[0] + y[1]*n0[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign0;
00098 
00099         localPos[0] = 1.0;
00100         localPos[1] = qPos;
00101         f.evaluate(localPos, y); 
00102         out[3] += (y[0]*n1[0] + y[1]*n1[1])*it->weight()*sign1;    
00103         out[4] += (y[0]*n1[0] + y[1]*n1[1])*(1.0 - 2.0*qPos)*it->weight();
00104         out[5] += (y[0]*n1[0] + y[1]*n1[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign1;
00105 
00106         localPos[0] = qPos;
00107         localPos[1] = 0.0;
00108         f.evaluate(localPos, y); 
00109         out[6] += (y[0]*n2[0] + y[1]*n2[1])*it->weight()*sign2;    
00110         out[7] += (y[0]*n2[0] + y[1]*n2[1])*(1.0 - 2.0*qPos)*it->weight();
00111         out[8] += (y[0]*n2[0] + y[1]*n2[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign2;
00112 
00113         localPos[0] = qPos;
00114         localPos[1] = 1.0;
00115         f.evaluate(localPos, y); 
00116         out[9]  += (y[0]*n3[0] + y[1]*n3[1])*it->weight()*sign3;
00117         out[10] += (y[0]*n3[0] + y[1]*n3[1])*(2.0*qPos - 1.0)*it->weight();
00118         out[11] += (y[0]*n3[0] + y[1]*n3[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign3;
00119       }
00120 
00121       const QuadratureRule<Vector,2>& rule2 = QuadratureRules<Vector,2>::rule(GeometryType(GeometryType::cube,2), qOrder);
00122       
00123       for (typename QuadratureRule<Vector,2>::const_iterator it = rule2.begin();
00124            it != rule2.end(); ++it)
00125       {
00126         FieldVector<double,2> qPos = it->position();
00127 
00128         f.evaluate(qPos, y); 
00129         out[12] += y[0]*it->weight();
00130         out[13] += y[1]*it->weight();
00131         out[14] += y[0]*qPos[0]*it->weight();
00132         out[15] += y[1]*qPos[0]*it->weight();
00133         out[16] += y[0]*qPos[1]*it->weight();
00134         out[17] += y[1]*qPos[1]*it->weight();
00135         out[18] += y[0]*qPos[0]*qPos[1]*it->weight();
00136         out[19] += y[1]*qPos[0]*qPos[1]*it->weight();
00137         out[20] += y[0]*qPos[1]*qPos[1]*it->weight();
00138         out[21] += y[1]*qPos[0]*qPos[0]*it->weight();
00139         out[22] += y[0]*qPos[0]*qPos[1]*qPos[1]*it->weight();
00140         out[23] += y[1]*qPos[0]*qPos[0]*qPos[1]*it->weight();
00141       }
00142     }
00143 
00144 private:
00145     typename LB::Traits::RangeFieldType sign0, sign1, sign2, sign3;
00146     typename LB::Traits::DomainType n0, n1, n2, n3;
00147   };
00148 }
00149 #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS2Q2DLOCALINTERPOLATION_HH