dune-localfunctions
2.2.0
|
00001 #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS12DLOCALINTERPOLATION_HH 00002 #define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS12DLOCALINTERPOLATION_HH 00003 00004 #include <vector> 00005 00006 #include <dune/geometry/quadraturerules.hh> 00007 00008 namespace Dune 00009 { 00010 00019 template<class LB> 00020 class RT12DLocalInterpolation 00021 { 00022 00023 public: 00025 RT12DLocalInterpolation () 00026 { 00027 sign0 = sign1 = sign2 = 1.0; 00028 } 00029 00035 RT12DLocalInterpolation (unsigned int s) 00036 { 00037 sign0 = sign1 = sign2 = 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 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 typedef typename LB::Traits::DomainFieldType Vector; 00075 typename F::Traits::RangeType y; 00076 00077 out.resize(8); 00078 fill(out.begin(), out.end(), 0.0); 00079 00080 const int qOrder1 = 4; 00081 const Dune::QuadratureRule<Scalar,1>& rule1 = Dune::QuadratureRules<Scalar,1>::rule(Dune::GeometryType(Dune::GeometryType::simplex,1), qOrder1); 00082 00083 for (typename Dune::QuadratureRule<Scalar,1>::const_iterator it = rule1.begin(); 00084 it != rule1.end(); ++it) 00085 { 00086 Scalar qPos = it->position(); 00087 typename LB::Traits::DomainType localPos; 00088 00089 localPos[0] = qPos; 00090 localPos[1] = 0.0; 00091 f.evaluate(localPos, y); 00092 out[0] += (y[0]*n0[0] + y[1]*n0[1])*it->weight()*sign0/c0; 00093 out[3] += (y[0]*n0[0] + y[1]*n0[1])*(2.0*qPos - 1.0)*it->weight()/c0; 00094 00095 localPos[0] = 0.0; 00096 localPos[1] = qPos; 00097 f.evaluate(localPos, y); 00098 out[1] += (y[0]*n1[0] + y[1]*n1[1])*it->weight()*sign1/c1; 00099 out[4] += (y[0]*n1[0] + y[1]*n1[1])*(1.0 - 2.0*qPos)*it->weight()/c1; 00100 00101 localPos[0] = 1.0 - qPos; 00102 localPos[1] = qPos; 00103 f.evaluate(localPos, y); 00104 out[2] += (y[0]*n2[0] + y[1]*n2[1])*it->weight()*sign2/c2; 00105 out[5] += (y[0]*n2[0] + y[1]*n2[1])*(2.0*qPos - 1.0)*it->weight()/c2; 00106 } 00107 00108 const int qOrder2 = 8; 00109 const Dune::QuadratureRule<Vector,2>& rule2 = Dune::QuadratureRules<Vector,2>::rule(Dune::GeometryType(Dune::GeometryType::simplex,2), qOrder2); 00110 00111 for (typename Dune::QuadratureRule<Vector,2>::const_iterator it = rule2.begin(); 00112 it != rule2.end(); ++it) 00113 { 00114 Dune::FieldVector<double,2> qPos = it->position(); 00115 00116 f.evaluate(qPos, y); 00117 out[6] += y[0]*it->weight(); 00118 out[7] += y[1]*it->weight(); 00119 } 00120 } 00121 00122 private: 00123 typename LB::Traits::RangeFieldType sign0,sign1,sign2; 00124 typename LB::Traits::DomainType n0,n1,n2; 00125 typename LB::Traits::RangeFieldType c0,c1,c2; 00126 }; 00127 } 00128 #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS12DLOCALINTERPOLATION_HH