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