dune-localfunctions
2.2.0
|
00001 #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1Q3DLOCALINTERPOLATION_HH 00002 #define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1Q3DLOCALINTERPOLATION_HH 00003 00004 #include <vector> 00005 00006 #include <dune/geometry/quadraturerules.hh> 00007 00008 namespace Dune 00009 { 00018 template<class LB> 00019 class RT1Q3DLocalInterpolation 00020 { 00021 00022 public: 00024 RT1Q3DLocalInterpolation () 00025 { 00026 sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0; 00027 } 00028 00034 RT1Q3DLocalInterpolation (unsigned int s) 00035 { 00036 sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 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 if (s & 8) 00050 { 00051 sign3 = -1.0; 00052 } 00053 if (s & 16) 00054 { 00055 sign4 = -1.0; 00056 } 00057 if (s & 32) 00058 { 00059 sign5 = -1.0; 00060 } 00061 00062 n0[0] = -1.0; 00063 n0[1] = 0.0; 00064 n0[2] = 0.0; 00065 n1[0] = 1.0; 00066 n1[1] = 0.0; 00067 n1[2] = 0.0; 00068 n2[0] = 0.0; 00069 n2[1] = -1.0; 00070 n2[2] = 0.0; 00071 n3[0] = 0.0; 00072 n3[1] = 1.0; 00073 n3[2] = 0.0; 00074 n4[0] = 0.0; 00075 n4[1] = 0.0; 00076 n4[2] = -1.0; 00077 n5[0] = 0.0; 00078 n5[1] = 0.0; 00079 n5[2] = 1.0; 00080 } 00081 00090 template<class F, class C> 00091 void interpolate (const F& f, std::vector<C>& out) const 00092 { 00093 // f gives v*outer normal at a point on the edge! 00094 typedef typename LB::Traits::RangeFieldType Scalar; 00095 typedef typename LB::Traits::DomainFieldType Vector; 00096 typename F::Traits::RangeType y; 00097 00098 out.resize(36); 00099 fill(out.begin(), out.end(), 0.0); 00100 00101 const int qOrder = 3; 00102 const QuadratureRule<Scalar,2>& rule1 = QuadratureRules<Scalar,2>::rule(GeometryType(GeometryType::cube,2), qOrder); 00103 00104 for (typename QuadratureRule<Scalar,2>::const_iterator it = rule1.begin(); 00105 it != rule1.end(); ++it) 00106 { 00107 Dune::FieldVector<Scalar,2> qPos = it->position(); 00108 typename LB::Traits::DomainType localPos; 00109 00110 localPos[0] = 0.0; 00111 localPos[1] = qPos[0]; 00112 localPos[2] = qPos[1]; 00113 f.evaluate(localPos, y); 00114 out[0] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*it->weight()*sign0; 00115 out[6] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[0] - 1.0)*it->weight(); 00116 out[12] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[1] - 1.0)*it->weight(); 00117 out[18] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*it->weight(); 00118 00119 localPos[0] = 1.0; 00120 localPos[1] = qPos[0]; 00121 localPos[2] = qPos[1]; 00122 f.evaluate(localPos, y); 00123 out[1] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*it->weight()*sign1; 00124 out[7] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[0])*it->weight(); 00125 out[13] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[1])*it->weight(); 00126 out[19] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*it->weight(); 00127 00128 localPos[0] = qPos[0]; 00129 localPos[1] = 0.0; 00130 localPos[2] = qPos[1]; 00131 f.evaluate(localPos, y); 00132 out[2] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*it->weight()*sign2; 00133 out[8] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(1.0 - 2.0*qPos[0])*it->weight(); 00134 out[14] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(2.0*qPos[1] - 1.0)*it->weight(); 00135 out[20] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*it->weight(); 00136 00137 localPos[0] = qPos[0]; 00138 localPos[1] = 1.0; 00139 localPos[2] = qPos[1]; 00140 f.evaluate(localPos, y); 00141 out[3] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*it->weight()*sign3; 00142 out[9] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(2.0*qPos[0] - 1.0)*it->weight(); 00143 out[15] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(1.0 - 2.0*qPos[1])*it->weight(); 00144 out[21] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*it->weight(); 00145 00146 localPos[0] = qPos[0]; 00147 localPos[1] = qPos[1]; 00148 localPos[2] = 0.0; 00149 f.evaluate(localPos, y); 00150 out[4] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*it->weight()*sign4; 00151 out[10] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[0])*it->weight(); 00152 out[16] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[1])*it->weight(); 00153 out[22] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*it->weight(); 00154 00155 localPos[0] = qPos[0]; 00156 localPos[1] = qPos[1]; 00157 localPos[2] = 1.0; 00158 f.evaluate(localPos, y); 00159 out[5] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*it->weight()*sign5; 00160 out[11] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[0] - 1.0)*it->weight(); 00161 out[17] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[1] - 1.0)*it->weight(); 00162 out[23] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*it->weight(); 00163 } 00164 00165 const QuadratureRule<Vector,3>& rule2 = QuadratureRules<Vector,3>::rule(GeometryType(GeometryType::cube,3), qOrder); 00166 for (typename QuadratureRule<Vector,3>::const_iterator it = rule2.begin(); 00167 it != rule2.end(); ++it) 00168 { 00169 FieldVector<double,3> qPos = it->position(); 00170 00171 f.evaluate(qPos, y); 00172 out[24] += y[0]*it->weight(); 00173 out[25] += y[1]*it->weight(); 00174 out[26] += y[2]*it->weight(); 00175 out[27] += y[0]*qPos[1]*it->weight(); 00176 out[28] += y[0]*qPos[2]*it->weight(); 00177 out[29] += y[1]*qPos[0]*it->weight(); 00178 out[30] += y[1]*qPos[2]*it->weight(); 00179 out[31] += y[2]*qPos[0]*it->weight(); 00180 out[32] += y[2]*qPos[1]*it->weight(); 00181 out[33] += y[0]*qPos[1]*qPos[2]*it->weight(); 00182 out[34] += y[1]*qPos[0]*qPos[2]*it->weight(); 00183 out[35] += y[2]*qPos[0]*qPos[1]*it->weight(); 00184 } 00185 } 00186 00187 private: 00188 typename LB::Traits::RangeFieldType sign0, sign1, sign2, sign3, sign4, sign5; 00189 typename LB::Traits::DomainType n0, n1, n2, n3, n4, n5; 00190 }; 00191 } 00192 #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1Q3DLOCALINTERPOLATION_HH