dune-localfunctions
2.2.0
|
00001 #ifndef DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI12DLOCALINTERPOLATION_HH 00002 #define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI12DLOCALINTERPOLATION_HH 00003 00004 #include <vector> 00005 00006 #include <dune/geometry/quadraturerules.hh> 00007 00008 namespace Dune 00009 { 00018 template<class LB> 00019 class BDM12DLocalInterpolation 00020 { 00021 00022 public: 00024 BDM12DLocalInterpolation () 00025 { 00026 sign0 = sign1 = sign2 = 1.0; 00027 } 00028 00034 BDM12DLocalInterpolation (unsigned int s) 00035 { 00036 sign0 = sign1 = sign2 = 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 00050 n0[0] = 0.0; 00051 n0[1] = -1.0; 00052 n1[0] = -1.0; 00053 n1[1] = 0.0; 00054 n2[0] = 1.0/sqrt(2.0); 00055 n2[1] = 1.0/sqrt(2.0); 00056 c0 = 0.5*n0[0] - 1.0*n0[1]; 00057 c1 = -1.0*n1[0] + 0.5*n1[1]; 00058 c2 = 0.5*n2[0] + 0.5*n2[1]; 00059 } 00060 00069 template<typename F, typename C> 00070 void interpolate (const F& f, std::vector<C>& out) const 00071 { 00072 // f gives v*outer normal at a point on the edge! 00073 typedef typename LB::Traits::RangeFieldType Scalar; 00074 typename F::Traits::RangeType y; 00075 00076 out.resize(6); 00077 fill(out.begin(), out.end(), 0.0); 00078 00079 const int qOrder = 4; 00080 const Dune::QuadratureRule<Scalar,1>& rule = Dune::QuadratureRules<Scalar,1>::rule(Dune::GeometryType(Dune::GeometryType::simplex,1), qOrder); 00081 00082 for (typename Dune::QuadratureRule<Scalar,1>::const_iterator it=rule.begin(); it!=rule.end(); ++it) 00083 { 00084 Scalar qPos = it->position(); 00085 typename LB::Traits::DomainType localPos; 00086 00087 localPos[0] = qPos; 00088 localPos[1] = 0.0; 00089 f.evaluate(localPos, y); 00090 out[0] += (y[0]*n0[0] + y[1]*n0[1])*it->weight()*sign0/c0; 00091 out[3] += (y[0]*n0[0] + y[1]*n0[1])*(2.0*qPos - 1.0)*it->weight()/c0; 00092 00093 localPos[0] = 0.0; 00094 localPos[1] = qPos; 00095 f.evaluate(localPos, y); 00096 out[1] += (y[0]*n1[0] + y[1]*n1[1])*it->weight()*sign1/c1; 00097 out[4] += (y[0]*n1[0] + y[1]*n1[1])*(1.0 - 2.0*qPos)*it->weight()/c1; 00098 00099 localPos[0] = 1.0 - qPos; 00100 localPos[1] = qPos; 00101 f.evaluate(localPos, y); 00102 out[2] += (y[0]*n2[0] + y[1]*n2[1])*it->weight()*sign2/c2; 00103 out[5] += (y[0]*n2[0] + y[1]*n2[1])*(2.0*qPos - 1.0)*it->weight()/c2; 00104 } 00105 } 00106 00107 private: 00108 typename LB::Traits::RangeFieldType sign0,sign1,sign2; 00109 typename LB::Traits::DomainType n0,n1,n2; 00110 typename LB::Traits::RangeFieldType c0,c1,c2; 00111 }; 00112 } 00113 00114 #endif // DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI12DLOCALINTERPOLATION_HH