dune-localfunctions
2.2.0
|
00001 // -*- tab-width: 4; indent-tabs-mode: nil -*- 00002 #ifndef DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH 00003 #define DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH 00004 00009 #include <vector> 00010 #include <dune/common/fvector.hh> 00011 #include <dune/common/fmatrix.hh> 00012 00013 #include <dune/localfunctions/common/localbasis.hh> 00014 #include <dune/localfunctions/common/localkey.hh> 00015 00016 namespace Dune 00017 { 00018 template<class D, class R, int dim> 00019 class HierarchicalSimplexP2WithElementBubbleLocalBasis 00020 { 00021 public: 00022 HierarchicalSimplexP2WithElementBubbleLocalBasis() 00023 { 00024 DUNE_THROW(Dune::NotImplemented,"HierarchicalSimplexP2LocalBasis not implemented for dim > 3."); 00025 } 00026 }; 00027 00042 template<class D, class R> 00043 class HierarchicalSimplexP2WithElementBubbleLocalBasis<D,R,1> 00044 { 00045 public: 00047 typedef LocalBasisTraits<D,1,Dune::FieldVector<D,1>,R,1,Dune::FieldVector<R,1>, 00048 Dune::FieldMatrix<R,1,1> > Traits; 00049 00051 unsigned int size () const 00052 { 00053 return 3; 00054 } 00055 00057 inline void evaluateFunction (const typename Traits::DomainType& in, 00058 std::vector<typename Traits::RangeType>& out) const 00059 { 00060 out.resize(3); 00061 00062 out[0] = 1-in[0]; 00063 out[1] = in[0]; 00064 out[2] = 1-4*(in[0]-0.5)*(in[0]-0.5); 00065 } 00066 00068 inline void 00069 evaluateJacobian (const typename Traits::DomainType& in, // position 00070 std::vector<typename Traits::JacobianType>& out) const // return value 00071 { 00072 out.resize(3); 00073 00074 out[0][0][0] = -1; 00075 out[1][0][0] = 1; 00076 out[2][0][0] = 4-8*in[0]; 00077 } 00078 00081 unsigned int order () const 00082 { 00083 return 2; 00084 } 00085 00086 }; 00087 00108 template<class D, class R> 00109 class HierarchicalSimplexP2WithElementBubbleLocalBasis<D,R,2> 00110 { 00111 public: 00113 typedef LocalBasisTraits<D,2,Dune::FieldVector<D,2>,R,1,Dune::FieldVector<R,1>, 00114 Dune::FieldMatrix<R,1,2> > Traits; 00115 00117 unsigned int size () const 00118 { 00119 return 7; 00120 } 00121 00123 inline void evaluateFunction (const typename Traits::DomainType& in, 00124 std::vector<typename Traits::RangeType>& out) const 00125 { 00126 out.resize(7); 00127 00128 out[0] = 1 - in[0] - in[1]; 00129 out[1] = 4*in[0]*(1-in[0]-in[1]); 00130 out[2] = in[0]; 00131 out[3] = 4*in[1]*(1-in[0]-in[1]); 00132 out[4] = 4*in[0]*in[1]; 00133 out[5] = in[1]; 00134 out[6] = 27*in[0]*in[1]*(1-in[0]-in[1]); 00135 00136 } 00137 00139 inline void 00140 evaluateJacobian (const typename Traits::DomainType& in, // position 00141 std::vector<typename Traits::JacobianType>& out) const // return value 00142 { 00143 out.resize(7); 00144 00145 out[0][0][0] = -1; out[0][0][1] = -1; 00146 out[1][0][0] = 4-8*in[0]-4*in[1]; out[1][0][1] = -4*in[0]; 00147 out[2][0][0] = 1; out[2][0][1] = 0; 00148 out[3][0][0] = -4*in[1]; out[3][0][1] = 4-4*in[0]-8*in[1]; 00149 out[4][0][0] = 4*in[1]; out[4][0][1] = 4*in[0]; 00150 out[5][0][0] = 0; out[5][0][1] = 1; 00151 00152 // Cubic bubble 00153 out[6][0][0] = 27 * in[1] * (1 - 2*in[0] - in[1]); 00154 out[6][0][1] = 27 * in[0] * (1 - 2*in[1] - in[0]); 00155 00156 } 00157 00160 unsigned int order () const 00161 { 00162 return 3; 00163 } 00164 00165 }; 00166 00191 template<class D, class R> 00192 class HierarchicalSimplexP2WithElementBubbleLocalBasis<D,R,3> 00193 { 00194 public: 00196 typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>, 00197 Dune::FieldMatrix<R,1,3> > Traits; 00198 00200 unsigned int size () const 00201 { 00202 return 11; 00203 } 00204 00206 void evaluateFunction (const typename Traits::DomainType& in, 00207 std::vector<typename Traits::RangeType>& out) const 00208 { 00209 out.resize(10); 00210 00211 out[0] = 1 - in[0] - in[1] - in[2]; 00212 out[1] = 4 * in[0] * (1 - in[0] - in[1] - in[2]); 00213 out[2] = in[0]; 00214 out[3] = 4 * in[1] * (1 - in[0] - in[1] - in[2]); 00215 out[4] = 4 * in[0] * in[1]; 00216 out[5] = in[1]; 00217 out[6] = 4 * in[2] * (1 - in[0] - in[1] - in[2]); 00218 out[7] = 4 * in[0] * in[2]; 00219 out[8] = 4 * in[1] * in[2]; 00220 out[9] = in[2]; 00221 00222 // quartic element bubble 00223 out[10] = 81*in[0]*in[1]*in[2]*(1-in[0]-in[1]-in[2]); 00224 } 00225 00227 void evaluateJacobian (const typename Traits::DomainType& in, // position 00228 std::vector<typename Traits::JacobianType>& out) const // return value 00229 { 00230 out.resize(10); 00231 00232 out[0][0][0] = -1; out[0][0][1] = -1; out[0][0][2] = -1; 00233 out[1][0][0] = 4-8*in[0]-4*in[1]-4*in[2]; out[1][0][1] = -4*in[0]; out[1][0][2] = -4*in[0]; 00234 out[2][0][0] = 1; out[2][0][1] = 0; out[2][0][2] = 0; 00235 out[3][0][0] = -4*in[1]; out[3][0][1] = 4-4*in[0]-8*in[1]-4*in[2]; out[3][0][2] = -4*in[1]; 00236 out[4][0][0] = 4*in[1]; out[4][0][1] = 4*in[0]; out[4][0][2] = 0; 00237 out[5][0][0] = 0; out[5][0][1] = 1; out[5][0][2] = 0; 00238 out[6][0][0] = -4*in[2]; out[6][0][1] = -4*in[2]; out[6][0][2] = 4-4*in[0]-4*in[1]-8*in[2]; 00239 out[7][0][0] = 4*in[2]; out[7][0][1] = 0; out[7][0][2] = 4*in[0]; 00240 out[8][0][0] = 0; out[8][0][1] = 4*in[2]; out[8][0][2] = 4*in[1]; 00241 out[9][0][0] = 0; out[9][0][1] = 0; out[9][0][2] = 1; 00242 00243 out[10][0][0] = 81 * in[1] * in[2] * (1 - 2*in[0] - in[1] - in[2]); 00244 out[10][0][1] = 81 * in[0] * in[2] * (1 - in[0] - 2*in[1] - in[2]); 00245 out[10][0][2] = 81 * in[0] * in[1] * (1 - in[0] - in[1] - 2*in[2]); 00246 } 00247 00248 00251 unsigned int order () const 00252 { 00253 return 4; 00254 } 00255 00256 }; 00257 00258 00284 template <int dim> 00285 class HierarchicalSimplexP2WithElementBubbleLocalCoefficients 00286 { 00287 // The binomial coefficient: dim+1 over 1 00288 static const int numVertices = dim+1; 00289 00290 // The binomial coefficient: dim+1 over 2 00291 static const int numEdges = (dim+1)*dim / 2; 00292 00293 public: 00295 HierarchicalSimplexP2WithElementBubbleLocalCoefficients () 00296 : li(numVertices+numEdges + 1) 00297 { 00298 if (dim!=2) 00299 DUNE_THROW(NotImplemented, "only for 2d"); 00300 00301 li[0] = Dune::LocalKey(0,2,0); // Vertex (0,0) 00302 li[1] = Dune::LocalKey(0,1,0); // Edge (0.5, 0) 00303 li[2] = Dune::LocalKey(1,2,0); // Vertex (1,0) 00304 li[3] = Dune::LocalKey(1,1,0); // Edge (0, 0.5) 00305 li[4] = Dune::LocalKey(2,1,0); // Edge (0.5, 0.5) 00306 li[5] = Dune::LocalKey(2,2,0); // Vertex (0,1) 00307 li[6] = Dune::LocalKey(0,0,0); // Element (1/3, 1/3) 00308 } 00309 00311 size_t size () const 00312 { 00313 return numVertices+numEdges + 1; 00314 } 00315 00317 const Dune::LocalKey& localKey (size_t i) const 00318 { 00319 return li[i]; 00320 } 00321 00322 private: 00323 std::vector<Dune::LocalKey> li; 00324 }; 00325 00326 template<class LB> 00327 class HierarchicalSimplexP2WithElementBubbleLocalInterpolation 00328 { 00329 public: 00330 00332 template<typename F, typename C> 00333 void interpolate (const F& f, std::vector<C>& out) const 00334 { 00335 typename LB::Traits::DomainType x; 00336 typename LB::Traits::RangeType y; 00337 00338 out.resize(7); 00339 00340 // vertices 00341 x[0] = 0.0; x[1] = 0.0; f.evaluate(x,y); out[0] = y; 00342 x[0] = 1.0; x[1] = 0.0; f.evaluate(x,y); out[2] = y; 00343 x[0] = 0.0; x[1] = 1.0; f.evaluate(x,y); out[5] = y; 00344 00345 // edge bubbles 00346 x[0] = 0.5; x[1] = 0.0; f.evaluate(x,y); 00347 out[1] = y - out[0]*(1-x[0]) - out[2]*x[0]; 00348 00349 x[0] = 0.0; x[1] = 0.5; f.evaluate(x,y); 00350 out[3] = y - out[0]*(1-x[1]) - out[5]*x[1]; 00351 00352 x[0] = 0.5; x[1] = 0.5; f.evaluate(x,y); 00353 out[4] = y - out[2]*x[0] - out[5]*x[1]; 00354 00355 // element bubble 00356 x[0] = 1.0/3; x[1] = 1.0/3; f.evaluate(x,y); 00357 00359 HierarchicalSimplexP2WithElementBubbleLocalBasis<double,double,2> shapeFunctions; 00360 std::vector<typename LB::Traits::RangeType> sfValues; 00361 shapeFunctions.evaluateFunction(x, sfValues); 00362 00363 out[6] = y; 00364 for (int i=0; i<6; i++) 00365 out[6] -= out[i]*sfValues[i]; 00366 00367 } 00368 00369 }; 00370 00371 00372 } 00373 #endif