dune-localfunctions  2.2.0
hierarchicalsimplexp2withelementbubble.hh
Go to the documentation of this file.
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