pyramidshapefunctions.hh

00001 #ifndef DUNE_PYRAMIDSHAPEFUNCTIONS_HH
00002 #define DUNE_PYRAMIDSHAPEFUNCTIONS_HH
00003 
00004 #include<iostream>
00005 #include<dune/common/fvector.hh>
00006 #include<dune/common/exceptions.hh>
00007 #include<dune/common/geometrytype.hh>
00008 
00015 namespace Dune
00016 {
00017 
00023   /***********************************************************
00024    * The interface for  pyramid shape functions of arbitrary
00025    * order 
00026    ***********************************************************/
00027 
00028 
00029   /***********************************************************
00030    * P0 shape functions for pyramid
00031    ***********************************************************/
00032 
00037   template<typename C, typename T>
00038   class P0PyramidShapeFunction
00039   {
00040   public:
00041 
00042         // compile time sizes
00043         enum { dim=3 };    // maps from R^d
00044         enum { comps=1 };      // to R^1
00045 
00046         enum { m=1 }; // total number of basis functions
00047 
00048         // export types
00049         typedef C CoordType;
00050         typedef T ResultType;
00051         typedef P0PyramidShapeFunction ImplementationType;
00052 
00054         P0PyramidShapeFunction ()
00055          {         
00056              pos[0] = 2/5.0;
00057              pos[1] = 2/5.0;
00058              pos[2] = 1/5.0;
00059            
00060          }
00061 
00062 
00064         ResultType evaluateFunction (int comp, const FieldVector<CoordType,3>& x) const
00065         {
00066           return 1;
00067         }
00068 
00070         ResultType evaluateDerivative (int comp, int dir, const FieldVector<CoordType,3>& x) const
00071         {
00072           return 0;
00073         }
00074 
00076         int localindex (int comp) const
00077         {
00078           return 0;
00079         }
00080 
00082         int codim () const
00083         {
00084           return 0;
00085         }
00086 
00088         int entity () const
00089         {
00090           return 0;
00091         }
00092 
00094         int entityindex () const
00095         {
00096           return 0;
00097         }
00098 
00100         const FieldVector<CoordType,dim>& position () const
00101         {
00102           return pos;
00103         }
00104         
00105   private:
00106         FieldVector<CoordType,3> pos;
00107   };
00108 
00109 
00110   template<typename C, typename T, typename S>
00111   class P0PyramidShapeFunctionSet
00112   {
00113   public:
00114 
00115         // compile time sizes
00116         enum { dim=3 };    // maps from R^d
00117         enum { comps=1 };      // to R^1
00118 
00119         enum { m=1 }; // total number of basis functions
00120 
00121         // export types
00122         typedef C CoordType;
00123         typedef T ResultType;
00124         typedef S value_type;
00125         typedef typename S::ImplementationType Imp; // Imp is either S or derived from S
00126 
00128         P0PyramidShapeFunctionSet ()
00129         {
00130           sf = Imp(); // assignment of derived class objects defined in wrapper
00131         }
00132 
00134         int size () const
00135         {
00136           return 1;
00137         }
00138 
00140         int size (int entity, int codim) const
00141         {
00142           if (codim==0) return 1; else return 0;
00143         }
00144         
00146         const value_type& operator[] (int i) const
00147         {
00148           return sf; // ok derived class reference goes for base class reference
00149         }
00150 
00152         int order () const
00153         {
00154           return 0;
00155         }
00156 
00158         GeometryType type () const
00159         {
00160       static GeometryType pyramid(GeometryType::pyramid, dim);
00161           return pyramid;
00162         }
00163 
00164   private:
00165         S sf;
00166   };
00167 
00168 
00169 /***********************************************************
00170  * P1 shape functions for pyramid
00171  ***********************************************************/
00172 
00177   template<typename C, typename T>
00178   class P1PyramidShapeFunction
00179   {
00180   public:
00181 
00182         // compile time sizes
00183         enum { dim=3 };    // maps from R^d
00184         enum { comps=1 };      // to R^1
00185 
00186         enum { m=5 }; // total number of basis functions
00187 
00188         // export types
00189         typedef C CoordType;
00190         typedef T ResultType;
00191         typedef P1PyramidShapeFunction ImplementationType;
00192 
00194         P1PyramidShapeFunction (int i)
00195         {          
00196           number =i;
00197 
00198           switch (i)
00199                 {
00200                 case 0:
00201                   pos[0]=0.0;
00202                   pos[1]=0.0;
00203                   pos[2]=0.0;
00204                   break;
00205                 case 1:
00206                   pos[0]=1.0;
00207                   pos[1]=0.0;
00208                   pos[2]=0.0;
00209                   break;
00210                 case 2:
00211                   pos[0]=1.0;
00212                   pos[1]=1.0;
00213                   pos[2]=0.0;
00214                   break;
00215                 case 3:
00216                   pos[0]=0.0;
00217                   pos[1]=1.0;
00218                   pos[2]=0.0;
00219                   break;
00220                 case 4:
00221                   pos[0]=0.0;
00222                   pos[1]=0.0;
00223                   pos[2]=1.0;
00224                   break;
00225                 default:
00226                   DUNE_THROW(RangeError, "wrong no of shape fns in Pyramid?");
00227                   break;
00228                 }
00229            
00230         }
00231 
00233         P1PyramidShapeFunction ()
00234     {}
00235 
00237         ResultType evaluateFunction (int comp, const FieldVector<CoordType,3>& x) const
00238         {
00239           ResultType phi;
00240           switch(number)
00241             {
00242             case 0:
00243               if(x[0] > x[1])
00244                         {
00245                           phi=((1-x[0])*(1-x[1])-x[2]*(1-x[1]));
00246                           return phi;
00247                         }
00248               else
00249                         {
00250                           phi=((1-x[0])*(1-x[1])-x[2]*(1-x[0]));
00251                           return phi;
00252                         }
00253               break;
00254 
00255             case 1:
00256               if(x[0] > x[1])
00257                         {
00258                           phi=(x[0]*(1-x[1])-x[2]*x[1]);
00259                           return phi;
00260                         }
00261               else
00262                         {
00263                           phi=(x[0]*(1-x[1])-x[2]*x[0]);
00264                           return phi;
00265                         }
00266               break;
00267 
00268             case 2:
00269               if(x[0] > x[1])
00270                         {
00271                           phi=(x[0]*x[1]+x[2]*x[1]);
00272                           return phi;
00273                         }
00274               else
00275                         {
00276                           phi=(x[0]*x[1]+x[2]*x[0]);
00277                           return phi;
00278                         }
00279               break;
00280 
00281             case 3:
00282               if(x[0] > x[1])
00283                         {
00284                           phi=((1-x[0])*x[1]-x[2]*x[1]);
00285                           return phi;
00286                         }
00287               else
00288                         { 
00289                           phi=((1-x[0])*x[1]-x[2]*x[0]);
00290                           return phi;
00291                         }
00292               break;
00293             case 4:
00294               
00295               phi=x[2];
00296               return phi;
00297         
00298               break;
00299             default:
00300               DUNE_THROW(RangeError, "wrong no of shape fns in Pyramid?");
00301               break;
00302             }
00303           return phi;
00304           
00305         }
00306 
00308     ResultType evaluateDerivative (int comp, int dir, const FieldVector<CoordType,3>& x) const
00309     {
00310       ResultType dphi;
00311       switch(number)
00312                 {
00313                 case 0:
00314                   if(x[0] > x[1])
00315                         {
00316                           switch(dir)
00317                                 {
00318                                 case 0:
00319                                   dphi=-1.0+x[1];
00320                                   return dphi;
00321                                   break;
00322                                 case 1:
00323                                   dphi=-1.0+x[0]+x[2];
00324                                   return dphi;
00325                                   break;
00326                                 case 2:
00327                                   dphi=-1.0+x[1];
00328                                   return dphi;
00329                                   break;
00330                                 }
00331 
00332                         }
00333                   else
00334                         {
00335                   
00336                           switch(dir)
00337                                 {
00338                                 case 0:
00339                                   dphi=-1.0+x[1]+x[2];
00340                                   return dphi;
00341                                   break;
00342                                 case 1:
00343                                   dphi=-1.0+x[0];
00344                                   return dphi;
00345                                   break;
00346                                 case 2:
00347                                   dphi=-1.0+x[0];
00348                                   return dphi;
00349                                   break;
00350                                 }
00351 
00352                         }
00353                   break;
00354 
00355                 case 1:
00356                   if(x[0] > x[1])
00357                         {
00358                  
00359                           switch(dir)
00360                                 {
00361                                 case 0:
00362                                   dphi=1.0-x[1];
00363                                   return dphi;
00364                                   break;
00365                                 case 1:
00366                                   dphi=-x[0]-x[2];
00367                                   return dphi;
00368                                   break;
00369                                 case 2:
00370                                   dphi=-x[1];
00371                                   return dphi;
00372                                   break;
00373                                 }
00374                         }
00375                   else
00376                         {
00377                  
00378                           switch(dir)
00379                                 {
00380                                 case 0:
00381                                   dphi=1.0-x[1]-x[2];
00382                                   return dphi;
00383                                   break;
00384                                 case 1:
00385                                   dphi=-x[0];
00386                                   return dphi;
00387                                   break;
00388                                 case 2:
00389                                   dphi=-x[0];
00390                                   return dphi;
00391                                   break;
00392                                 }
00393 
00394                         }
00395                   break;
00396 
00397                 case 2:
00398                   if(x[0] > x[1])
00399                         {
00400                           switch(dir)
00401                                 {
00402                                 case 0:
00403                                   dphi=x[1];
00404                                   return dphi;
00405                                   break;
00406                                 case 1:
00407                                   dphi=x[0]+x[2];
00408                                   return dphi;
00409                                   break;
00410                                 case 2:
00411                                   dphi=x[1];
00412                                   return dphi;
00413                                   break;
00414                                 }
00415 
00416                         }
00417                   else
00418                         {
00419                   
00420                           switch(dir)
00421                                 {
00422                                 case 0:
00423                                   dphi=x[1]+x[2];
00424                                   return dphi;
00425                                   break;
00426                                 case 1:
00427                                   dphi=x[0];
00428                                   return dphi;
00429                                   break;
00430                                 case 2:
00431                                   dphi=x[0];
00432                                   return dphi;
00433                                   break;
00434                                 }
00435 
00436                         }
00437                   break;
00438 
00439                 case 3:
00440                   if(x[0] > x[1])
00441                         {
00442                  
00443                           switch(dir)
00444                                 {
00445                                 case 0:
00446                                   dphi=-x[1];
00447                                   return dphi;
00448                                   break;
00449                                 case 1:
00450                                   dphi=1.0-x[0]-x[2];
00451                                   return dphi;
00452                                   break;
00453                                 case 2:
00454                                   dphi=-x[1];
00455                                   return dphi;
00456                                   break;
00457                                 }
00458 
00459                         }
00460                   else
00461                         { 
00462                  
00463                           switch(dir)
00464                                 {
00465                                 case 0:
00466                                   dphi=-x[1]-x[2];
00467                                   return dphi;
00468                                   break;
00469                                 case 1:
00470                                   dphi=1.0-x[0];
00471                                   return dphi;
00472                                   break;
00473                                 case 2:
00474                                   dphi=-x[0];
00475                                   return dphi;
00476                                   break;
00477                                 }
00478                         }
00479                   break;
00480                 case 4:
00481                   switch(dir)
00482                         {
00483                         case 0:
00484                           dphi=0;
00485                           return dphi;
00486                           break;
00487                         case 1:
00488                           dphi=0;
00489                           return dphi;
00490                           break;
00491                         case 2:
00492                           dphi=1.0;
00493                           return dphi;
00494                           break;
00495                         }
00496               
00497                   break;
00498                 default:
00499           DUNE_THROW(RangeError, "wrong no of shape fns in Pyramid?");
00500                   break;
00501                 }
00502       return 0;
00503     }
00504 
00506         int localindex (int comp) const
00507         {
00508           return number;
00509         }
00510 
00512         int codim () const
00513         {
00514           return dim;
00515         }
00516 
00518         int entity () const
00519         {
00520           return number;
00521         }
00522 
00524         int entityindex () const
00525         {
00526           return 0;
00527         }
00528 
00530         const FieldVector<CoordType,dim>& position () const
00531         {
00532           return pos;
00533         }
00534         
00535   private:
00536     int number;
00537         FieldVector<CoordType,3> pos;
00538   };
00539 
00540 
00541   template<typename C, typename T, typename S>
00542   class P1PyramidShapeFunctionSet
00543   {
00544   public:
00545 
00546         // compile time sizes
00547         enum { dim=3 };    // maps from R^d
00548         enum { comps=1 };      // to R^1
00549 
00550         enum { m=5 }; // total number of basis functions
00551 
00552         // export types
00553         typedef C CoordType;
00554         typedef T ResultType;
00555         typedef S value_type;
00556         typedef typename S::ImplementationType Imp; // Imp is either S or derived from S
00557 
00559         P1PyramidShapeFunctionSet ()
00560           {
00561                 for (int i=0; i<m; i++)
00562                   sf[i] = Imp(i); // assignment of derived class objects defined in wrapper
00563           }
00564 
00566         int size () const
00567           {
00568                 return m;
00569           }
00570 
00572         int size (int entity, int codim) const
00573           {
00574                 if (codim==dim) return m; else return 0;
00575           }
00576         
00578         const value_type& operator[] (int i) const
00579           {
00580                 return sf[i]; // ok derived class reference goes for base class reference
00581           }
00582 
00584         int order () const
00585           {
00586                 return 1;
00587           }
00588 
00590         GeometryType type () const
00591           {
00592         static GeometryType pyramid(GeometryType::pyramid, dim);
00593                 return pyramid;
00594           }
00595 
00596   private:
00597         S sf[m];
00598   };
00599 
00600 
00601 
00602 
00604     template<typename C, typename T>
00605   class P0PyramidShapeFunctionSetContainer
00606   {
00607   public:
00608         // compile time sizes
00609         enum { dim=3 };       
00610         enum { comps=1 };
00611         enum { maxsize=1 };
00612 
00613         // exported types
00614         typedef C CoordType; 
00615         typedef T ResultType;
00616         typedef P0PyramidShapeFunctionSet<C,T,P0PyramidShapeFunction<C,T> > value_type;
00617 
00618         const value_type& operator() (GeometryType type, int order) const
00619         {
00620 
00621           if(type.isPyramid()) return p0pyramid;
00622           DUNE_THROW(NotImplemented, "type not implemented yet");
00623         }
00624   private:
00625         value_type p0pyramid;
00626   };
00627  
00629 
00630   template<typename C, typename T>
00631   class P1PyramidShapeFunctionSetContainer
00632   {
00633   public:
00634         // compile time sizes
00635         enum { dim=3 };       
00636         enum { comps=1 };
00637         enum { maxsize=5 };
00638 
00639         // exported types
00640         typedef C CoordType; 
00641         typedef T ResultType;
00642         typedef P1PyramidShapeFunctionSet<C,T,P1PyramidShapeFunction<C,T> > value_type;
00643 
00644         const value_type& operator() (GeometryType type, int order) const
00645         {
00646 
00647           if(type.isPyramid()) return p1pyramid;
00648           DUNE_THROW(NotImplemented, "type not implemented yet");
00649         }
00650   private:
00651         value_type p1pyramid;
00652   };
00654 }
00655 #endif