dune-localfunctions  2.2.0
raviartthomas1q3dlocalinterpolation.hh
Go to the documentation of this file.
00001 #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1Q3DLOCALINTERPOLATION_HH
00002 #define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1Q3DLOCALINTERPOLATION_HH
00003 
00004 #include <vector>
00005 
00006 #include <dune/geometry/quadraturerules.hh>
00007 
00008 namespace Dune 
00009 {
00018   template<class LB>
00019   class RT1Q3DLocalInterpolation
00020   {
00021 
00022 public:
00024     RT1Q3DLocalInterpolation ()
00025     {
00026       sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0;
00027     }
00028 
00034     RT1Q3DLocalInterpolation (unsigned int s)
00035     {
00036       sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 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       if (s & 8)
00050       {
00051         sign3 = -1.0;
00052       }
00053       if (s & 16)
00054       {
00055         sign4 = -1.0;
00056       }
00057       if (s & 32)
00058       {
00059         sign5 = -1.0;
00060       }
00061 
00062       n0[0] = -1.0;
00063       n0[1] =  0.0;
00064       n0[2] =  0.0;
00065       n1[0] =  1.0;
00066       n1[1] =  0.0;
00067       n1[2] =  0.0;
00068       n2[0] =  0.0;
00069       n2[1] = -1.0;
00070       n2[2] =  0.0;
00071       n3[0] =  0.0;
00072       n3[1] =  1.0;
00073       n3[2] =  0.0;
00074       n4[0] =  0.0;
00075       n4[1] =  0.0;
00076       n4[2] = -1.0;
00077       n5[0] =  0.0;
00078       n5[1] =  0.0;
00079       n5[2] =  1.0;
00080     }
00081 
00090     template<class F, class C>
00091     void interpolate (const F& f, std::vector<C>& out) const
00092     {
00093       // f gives v*outer normal at a point on the edge!
00094       typedef typename LB::Traits::RangeFieldType Scalar;
00095       typedef typename LB::Traits::DomainFieldType Vector;
00096       typename F::Traits::RangeType y;
00097 
00098       out.resize(36);
00099       fill(out.begin(), out.end(), 0.0);
00100 
00101       const int qOrder = 3;
00102       const QuadratureRule<Scalar,2>& rule1 = QuadratureRules<Scalar,2>::rule(GeometryType(GeometryType::cube,2), qOrder);
00103 
00104       for (typename QuadratureRule<Scalar,2>::const_iterator it = rule1.begin();
00105            it != rule1.end(); ++it)
00106       {
00107         Dune::FieldVector<Scalar,2> qPos = it->position();
00108         typename LB::Traits::DomainType localPos;
00109 
00110         localPos[0] = 0.0;
00111         localPos[1] = qPos[0];
00112         localPos[2] = qPos[1];
00113         f.evaluate(localPos, y);
00114         out[0] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*it->weight()*sign0;
00115         out[6] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[0] - 1.0)*it->weight();
00116         out[12] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[1] - 1.0)*it->weight();
00117         out[18] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*it->weight();
00118         
00119         localPos[0] = 1.0;
00120         localPos[1] = qPos[0];
00121         localPos[2] = qPos[1];
00122         f.evaluate(localPos, y);
00123         out[1] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*it->weight()*sign1;
00124         out[7] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[0])*it->weight();
00125         out[13] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[1])*it->weight();
00126         out[19] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*it->weight();
00127         
00128         localPos[0] = qPos[0];
00129         localPos[1] = 0.0;
00130         localPos[2] = qPos[1];
00131         f.evaluate(localPos, y);
00132         out[2] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*it->weight()*sign2;
00133         out[8] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(1.0 - 2.0*qPos[0])*it->weight();
00134         out[14] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(2.0*qPos[1] - 1.0)*it->weight();
00135         out[20] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*it->weight();
00136         
00137         localPos[0] = qPos[0];
00138         localPos[1] = 1.0;
00139         localPos[2] = qPos[1];
00140         f.evaluate(localPos, y);
00141         out[3] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*it->weight()*sign3;
00142         out[9] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(2.0*qPos[0] - 1.0)*it->weight();
00143         out[15] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(1.0 - 2.0*qPos[1])*it->weight();
00144         out[21] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*it->weight();
00145         
00146         localPos[0] = qPos[0];
00147         localPos[1] = qPos[1];
00148         localPos[2] = 0.0;
00149         f.evaluate(localPos, y);
00150         out[4] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*it->weight()*sign4;
00151         out[10] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[0])*it->weight();
00152         out[16] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[1])*it->weight();
00153         out[22] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*it->weight();
00154         
00155         localPos[0] = qPos[0];
00156         localPos[1] = qPos[1];
00157         localPos[2] = 1.0;
00158         f.evaluate(localPos, y);
00159         out[5] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*it->weight()*sign5;
00160         out[11] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[0] - 1.0)*it->weight();
00161         out[17] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[1] - 1.0)*it->weight();
00162         out[23] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*it->weight();
00163       }
00164       
00165       const QuadratureRule<Vector,3>& rule2 = QuadratureRules<Vector,3>::rule(GeometryType(GeometryType::cube,3), qOrder);
00166       for (typename QuadratureRule<Vector,3>::const_iterator it = rule2.begin();
00167            it != rule2.end(); ++it)
00168       {
00169         FieldVector<double,3> qPos = it->position();
00170 
00171         f.evaluate(qPos, y);
00172         out[24] += y[0]*it->weight();
00173         out[25] += y[1]*it->weight();
00174         out[26] += y[2]*it->weight();
00175         out[27] += y[0]*qPos[1]*it->weight();
00176         out[28] += y[0]*qPos[2]*it->weight();
00177         out[29] += y[1]*qPos[0]*it->weight();
00178         out[30] += y[1]*qPos[2]*it->weight();
00179         out[31] += y[2]*qPos[0]*it->weight();
00180         out[32] += y[2]*qPos[1]*it->weight();
00181         out[33] += y[0]*qPos[1]*qPos[2]*it->weight();
00182         out[34] += y[1]*qPos[0]*qPos[2]*it->weight();
00183         out[35] += y[2]*qPos[0]*qPos[1]*it->weight();
00184       }
00185     }
00186 
00187 private:
00188     typename LB::Traits::RangeFieldType sign0, sign1, sign2, sign3, sign4, sign5;
00189     typename LB::Traits::DomainType n0, n1, n2, n3, n4, n5;
00190   };
00191 }
00192 #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1Q3DLOCALINTERPOLATION_HH