dune-localfunctions
2.2.0
|
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