hierarchicalcubeshapefunctions.hh

00001 // $Id: cubeshapefunctions.hh 336 2006-05-03 13:09:05Z oliver $
00002 
00003 #ifndef DUNE_HIERARCHICAL_CUBE_SHAPEFUNCTIONS_HH
00004 #define DUNE_HIERARCHICAL_CUBE_SHAPEFUNCTIONS_HH
00005 
00006 #include<iostream>
00007 #include<dune/common/fvector.hh>
00008 #include<dune/common/exceptions.hh>
00009 #include<dune/common/misc.hh>
00010 #include<dune/common/geometrytype.hh>
00011 #include<dune/grid/common/referenceelements.hh>
00012 
00018 namespace Dune
00019 {
00030   /***********************************************************
00031    * P2 shape functions for the cube of any dimension
00032    ***********************************************************/
00033 
00044   template<typename C, typename T, int d>
00045   class P2HierarchicalCubeShapeFunction
00046   {
00047   public:
00048 
00049         // compile time sizes
00050         enum { dim=d };    // maps from R^d
00051         enum { comps=1 };      // to R^1
00052 
00053         enum { m=Power_m_p<3,d>::power }; // total number of basis functions
00054 
00055         // export types
00056         typedef C CoordType;
00057         typedef T ResultType;
00058         typedef P2HierarchicalCubeShapeFunction ImplementationType;
00059 
00062         P2HierarchicalCubeShapeFunction (int no, int en, int co, const FieldVector<int,dim>& ipos) // make it the i'th shape function
00063         {
00064           num = no;
00065           ent = en;
00066           cod = co;
00067 
00068           // Find out if this is a bubble function
00069           bool isBubble = false;
00070           for (int j=0; j<dim; j++)
00071               isBubble |= ipos[j]==1;
00072 
00073           if (isBubble) {
00074 
00075               // Quadratic bubble functions
00076               for (int j=0; j<dim; j++) {
00077                   
00078                   if (ipos[j]==0) {
00079                       
00080                       a[j] = 1.0; 
00081                       b[j] = -3.0;
00082                       c[j] = 2.0;
00083                       pos[j] = 0;
00084                   }
00085                   if (ipos[j]==1) {
00086                       
00087                       a[j] = 0.0; 
00088                       b[j] = 4.0;
00089                       c[j] = -4.0;
00090                       pos[j] = 0.5;
00091                   }
00092                   if (ipos[j]==2) {
00093                       
00094                       a[j] = 0.0; 
00095                       b[j] = -1.0;
00096                       c[j] = 2.0;
00097                       pos[j] = 1.0;
00098                   }
00099               }
00100           } else {
00101               
00102               // P1 shape functions
00103               for (int j=0; j<dim; j++) {
00104 
00105                   if (ipos[j]==0) {
00106                       
00107                       a[j] = 1.0; 
00108                       b[j] = -1.0;
00109                       c[j] = 0.0;
00110                       pos[j] = 0;
00111                   }
00112                   if (ipos[j]==2) {
00113                       a[j] = 0.0; 
00114                       b[j] = 1.0;
00115                       c[j] = 0.0;
00116                       pos[j] = 1.0;
00117                   }
00118               }
00119 
00120           }
00121 
00122         }
00123 
00125         P2HierarchicalCubeShapeFunction ()
00126         {       }
00127 
00129         ResultType evaluateFunction (int comp, const FieldVector<CoordType,d>& x) const
00130         {
00131           ResultType phi = a[0]+x[0]*b[0] +x[0]*x[0]*c[0];
00132           for (int j=1; j<dim; j++) phi *= a[j]+x[j]*b[j]+x[j]*x[j]*c[j];
00133           return phi;
00134         }
00135 
00137         ResultType evaluateDerivative (int comp, int dir, const FieldVector<CoordType,d>& x) const
00138         {
00139           ResultType deriv = b[dir]+2*c[dir]*x[dir];
00140           for (int j=0; j<dim; j++)
00141                 if (j!=dir) deriv *= a[j]+x[j]*b[j]+x[j]*x[j]*c[j];
00142           return deriv;
00143         }
00144 
00146         int localindex (int comp) const
00147         {
00148           return num;
00149         }
00150 
00152         int codim () const
00153         {
00154           return cod;
00155         }
00156 
00158         int entity () const
00159         {
00160           return ent;
00161         }
00162 
00164         int entityindex () const
00165         {
00166           return 0;
00167         }
00168 
00170         const FieldVector<CoordType,dim>& position () const
00171         {
00172           return pos;
00173         }
00174         
00175   private:
00176         int num;
00177         int ent;
00178         int cod;
00179         ResultType a[dim], b[dim], c[dim]; // store coefficients for this shape function
00180         FieldVector<CoordType,d> pos;
00181   };
00182 
00183 
00188   template<typename C, typename T, int d, typename S>
00189   class P2HierarchicalCubeShapeFunctionSet
00190   {
00191   public:
00192 
00193         // compile time sizes
00194         enum { dim=d };    // maps from R^d
00195         enum { comps=1 };      // to R^1
00196 
00197         enum { m=Power_m_p<3,d>::power }; // total number of basis functions
00198 
00199         // export types
00200         typedef C CoordType;
00201         typedef T ResultType;
00202         typedef S value_type;
00203         typedef typename S::ImplementationType Imp; // Imp is either S or derived from S
00204 
00206         P2HierarchicalCubeShapeFunctionSet ()
00207         {
00208           ReferenceCube<C,d> cube;
00209           int i=0;
00210           for (int c=0; c<=dim; c++)
00211                 for (int e=0; e<cube.size(c); e++)
00212                   {
00213                         sf[i] = Imp(i,e,c,cube.iposition(e,c));
00214                         i++;
00215                   }
00216         }
00217 
00219         int size () const
00220         {
00221           return m;
00222         }
00223 
00225         int size (int entity, int codim) const
00226         {
00227           return 1;
00228         }
00229         
00231         const value_type& operator[] (int i) const
00232         {
00233           return sf[i]; // ok derived class reference goes for base class reference
00234         }
00235 
00237         int order () const
00238         {
00239           return 2;
00240         }
00241 
00243         GeometryType type () const
00244         {
00245       static GeometryType cube(GeometryType::cube, dim);
00246           return cube;
00247         }
00248 
00249   private:
00250         S sf[m];
00251   };
00252 
00253 
00255   template<typename C, typename T, int d>
00256   class P2HierarchicalCubeShapeFunctionSetContainer
00257   {
00258   public:
00259         // compile time sizes
00260         enum { dim=d };       
00261         enum { comps=1 };
00262         enum { maxsize=Power_m_p<3,d>::power };
00263 
00264         // exported types
00265         typedef C CoordType; 
00266         typedef T ResultType;
00267         typedef P2HierarchicalCubeShapeFunctionSet<C,T,d,P2HierarchicalCubeShapeFunction<C,T,d> > value_type;
00268 
00269         const value_type& operator() (GeometryType type, int order) const
00270         {
00271       if (type.isCube()) return p2cube;
00272           DUNE_THROW(NotImplemented, "type not implemented yet");
00273         }
00274   private:
00275         value_type p2cube;
00276   };
00277  
00278 
00279 
00281 }
00282 #endif