dune-localfunctions
2.2.0
|
00001 #ifndef DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1Q2DLOCALINTERPOLATION_HH 00002 #define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1Q2DLOCALINTERPOLATION_HH 00003 00004 #include <vector> 00005 00006 #include <dune/geometry/quadraturerules.hh> 00007 00008 namespace Dune 00009 { 00010 00019 template<class LB> 00020 class BDM1Q2DLocalInterpolation 00021 { 00022 00023 public: 00025 BDM1Q2DLocalInterpolation () 00026 { 00027 sign0 = sign1 = sign2 = sign3 = 1.0; 00028 } 00029 00035 BDM1Q2DLocalInterpolation (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(8); 00082 fill(out.begin(), out.end(), 0.0); 00083 00084 const int qOrder = 4; 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(); 00088 it != rule.end(); ++it) 00089 { 00090 Scalar qPos = it->position(); 00091 typename LB::Traits::DomainType localPos; 00092 00093 localPos[0] = 0.0; 00094 localPos[1] = qPos; 00095 f.evaluate(localPos, y); 00096 out[0] += (y[0]*n0[0] + y[1]*n0[1])*it->weight()*sign0; 00097 out[1] += (y[0]*n0[0] + y[1]*n0[1])*(2.0*qPos - 1.0)*it->weight(); 00098 00099 localPos[0] = 1.0; 00100 localPos[1] = qPos; 00101 f.evaluate(localPos, y); 00102 out[2] += (y[0]*n1[0] + y[1]*n1[1])*it->weight()*sign1; 00103 out[3] += (y[0]*n1[0] + y[1]*n1[1])*(1.0 - 2.0*qPos)*it->weight(); 00104 00105 localPos[0] = qPos; 00106 localPos[1] = 0.0; 00107 f.evaluate(localPos, y); 00108 out[4] += (y[0]*n2[0] + y[1]*n2[1])*it->weight()*sign2; 00109 out[5] += (y[0]*n2[0] + y[1]*n2[1])*(1.0 - 2.0*qPos)*it->weight(); 00110 00111 localPos[0] = qPos; 00112 localPos[1] = 1.0; 00113 f.evaluate(localPos, y); 00114 out[6] += (y[0]*n3[0] + y[1]*n3[1])*it->weight()*sign3; 00115 out[7] += (y[0]*n3[0] + y[1]*n3[1])*(2.0*qPos - 1.0)*it->weight(); 00116 } 00117 } 00118 00119 private: 00120 typename LB::Traits::RangeFieldType sign0, sign1, sign2, sign3; 00121 typename LB::Traits::DomainType n0, n1, n2, n3; 00122 }; 00123 } 00124 #endif // DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1Q2DLOCALINTERPOLATION_HH