hierarchicalprismshapefunctions.hh

00001 #ifndef DUNE_HIERARCHICAL_PRISM_SHAPEFUNCTIONS_HH
00002 #define DUNE_HIERARCHICAL_PRISM_SHAPEFUNCTIONS_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 
00024   /***********************************************************
00025    * Hierarchical P2 shape functions for prism
00026    ***********************************************************/
00027 
00033   /*
00034         polynomial of the form: ( a[0] + b[0]*x + b[1]*y + b[2]z) * (a[1] + c[0]*x + c[1]*y + c[2]*z)
00035         derrivative of the form:( aa[0][dir] + bb[0][dir]*x + bb[1][dir]*y + bb[2][dir]*z)
00036         dir=0 for x, dir=1 for y, dir=2 for z
00037   */
00038   template<typename C,typename T>
00039   class P2HierarchicalPrismShapeFunction
00040   {
00041   public:
00042 
00043         enum {dim=3};
00044         enum {comps=1};
00045         enum {m=18}; // number of basis functions
00046 
00047         typedef C CoordType;
00048         typedef T ResultType;
00049         typedef P2HierarchicalPrismShapeFunction ImplementationType;
00050   
00051         // i'th shape function
00052         P2HierarchicalPrismShapeFunction(int i)
00053         {
00054             number=i;
00055 
00056             // A renumbering of the shape functions might simplify this constructor
00057 
00058             if (i<3 || (i>=12 && i<15))
00059                 codim_=3;
00060             else if (i>=9 && i<12)
00061                 codim_=1;
00062             else
00063                 codim_=2;
00064 
00065             static const int entities[] = {0,1,2,    // vertices of the lower triangle
00066                                            1,2,0,    // edges of the lower triangle
00067                                            3, 4, 5,  // vertical edges
00068                                            2,3,1,    // quadrilateral sides
00069                                            3,4,5,    // vertices of the upper triangle
00070                                            7,8,6};   // edges of the upper triangle
00071 
00072             entity_ = entities[i];
00073 
00074             if (i<6) {
00075 
00076                 a1d = 1.0; 
00077                 b1d = -3.0;
00078                 c1d = 2.0;
00079                 pos[2] = 0;
00080             
00081             } else if (i<12) {
00082                 
00083                 a1d = 0.0; 
00084                 b1d = 4.0;
00085                 c1d = -4.0;
00086                 pos[2] = 0.5;
00087                         
00088             } else {
00089                 a1d = 0.0; 
00090                 b1d = -1.0;
00091                 c1d = 2.0;
00092                 pos[2] = 1.0;
00093             }
00094 
00095             switch(i) {
00096 
00097             case 0:
00098             case 6:
00099             case 12:
00100 
00101                 //--interpolation point associated with shape fn
00102                 pos[0]=0.0;
00103                 pos[1]=0.0;
00104                 //--
00105                 coeff=2;
00106                 a[0]=1.0;
00107                 a[1]=0.5;
00108                 b[0]=-1.0;
00109                 b[1]=-1.0;
00110                 c[0]=-1.0;
00111                 c[1]=-1.0;
00112                 //x derivative
00113                 aa[0][0]=-3;
00114                 bb[0][0]=4;
00115                 bb[1][0]=4;
00116                 //y derivative
00117                 aa[0][1]=-3;
00118                 bb[0][1]=4;
00119                 bb[1][1]=4;
00120                 break;
00121             case 1:
00122             case 7:
00123             case 13:
00124 
00125                 //--interpolation point associated with shape fn
00126                 pos[0]=1.0;
00127                 pos[1]=0.0;
00128                 //--
00129                 coeff=2;
00130                 a[0]=0.0;
00131                 a[1]=-0.5;
00132                 b[0]=1.0;
00133                 b[1]=0.0;
00134                 c[0]=1.0;
00135                 c[1]=0.0;
00136                 //x derivative
00137                 aa[0][0]=-1;
00138                 bb[0][0]=4;
00139                 bb[1][0]=0;
00140                 //y derivative
00141                 aa[0][1]=0;
00142                 bb[0][1]=0;
00143                 bb[1][1]=0;
00144                 break;
00145             case 2:
00146             case 8:
00147             case 14:
00148                 //--interpolation point associated with shape fn
00149                 pos[0]=0.0;
00150                 pos[1]=1.0;
00151                 //--
00152                 coeff=2;
00153                 a[0]=0.0;
00154                 a[1]=-0.5;
00155                 b[0]=0.0;
00156                 b[1]=1.0;
00157                 c[0]=0.0;
00158                 c[1]=1.0;
00159                 //x derivative
00160                 aa[0][0]=0;
00161                 bb[0][0]=0;
00162                 bb[1][0]=0;
00163                 //y derivative
00164                 aa[0][1]=-1;
00165                 bb[0][1]=0;
00166                 bb[1][1]=4;
00167                 break;
00168             case 3:
00169             case 9:
00170             case 15:
00171                 //--interpolation point associated with shape fn
00172                 pos[0]=0.5;
00173                 pos[1]=0.5;
00174                 //--
00175                 coeff=4;
00176                 a[0]=0.0;
00177                 a[1]=0.0;
00178                 b[0]=1.0;
00179                 b[1]=0.0;
00180                 c[0]=0.0;
00181                 c[1]=1.0;
00182                 //x derivative
00183                 aa[0][0]=0;
00184                 bb[0][0]=0;
00185                 bb[1][0]=4;
00186                 //y derivative
00187                 aa[0][1]=0;
00188                 bb[0][1]=4;
00189                 bb[1][1]=0;
00190                 break;
00191             case 4:
00192             case 10:
00193             case 16:
00194                 //--interpolation point associated with shape fn
00195                 pos[0]=0.0;
00196                 pos[1]=0.5;
00197                 //--
00198                 coeff=4;
00199                 a[0]=0.0;
00200                 a[1]=1.0;
00201                 b[0]=0.0;
00202                 b[1]=1.0;
00203                 c[0]=-1.0;
00204                 c[1]=-1.0;
00205                 //x derivative
00206                 aa[0][0]=0;
00207                 bb[0][0]=0;
00208                 bb[1][0]=-4;
00209                 //y derivative
00210                 aa[0][1]=4;
00211                 bb[0][1]=-4;
00212                 bb[1][1]=-8;
00213                 break;
00214             case 5:
00215             case 11:
00216             case 17:
00217                 //--interpolation point associated with shape fn
00218                 pos[0]=0.5;
00219                 pos[1]=0.0;
00220                 //--
00221                 coeff=4;
00222                 a[0]=0.0;
00223                 a[1]=1.0;
00224                 b[0]=1.0;
00225                 b[1]=0.0;
00226                 c[0]=-1.0;
00227                 c[1]=-1.0;
00228                 //x derivative
00229                 aa[0][0]=4;
00230                 bb[0][0]=-8;
00231                 bb[1][0]=-4;
00232                 //y derivative
00233                 aa[0][1]=0;
00234                 bb[0][1]=-4;
00235                 bb[1][1]=0;
00236                 break;
00237             default:
00238                 DUNE_THROW(RangeError, "wrong no of shape fns in Prism?");
00239                 break;
00240             }
00241         }
00243         P2HierarchicalPrismShapeFunction ()
00244     {}
00245   
00247         ResultType evaluateFunction (int comp, const FieldVector<CoordType,3>& x) const
00248         {
00249             // Part I: Evaluate a P2 function on the triangle
00250             ResultType phi1=a[0];
00251             ResultType phi2=a[1];
00252             
00253             for (int j=0;j<2;++j)
00254                 {
00255                     phi1 +=b[j]*x[j];
00256                     phi2 +=c[j]*x[j];
00257                 }
00258 
00259             // Part II: Evaluate a P2 function on the segment, multiply the results
00260             ResultType phi3 = a1d + x[2]*b1d + x[2]*x[2]*c1d;
00261 
00262             return coeff * phi1*phi2*phi3;
00263    
00264    
00265         } 
00266   
00267 
00268         ResultType evaluateDerivative (int comp, int dir, const FieldVector<CoordType,3>& x) const
00269   
00270         {
00271             if (dir==2) {
00272                 ResultType phi1=a[0];
00273                 ResultType phi2=a[1];
00274                 
00275                 for (int j=0;j<2;++j) {
00276                     phi1 +=b[j]*x[j];
00277                     phi2 +=c[j]*x[j];
00278                 }
00279                 
00280                 ResultType phi = coeff * phi1 * phi2;
00281                 return phi * (b1d + 2*c1d*x[2]);
00282 
00283             } 
00284             
00285             ResultType deriv=aa[0][dir];
00286             
00287             for (int j=0;j<2;++j)
00288                 deriv +=bb[j][dir]*x[j];
00289             
00290             return deriv * (a1d + x[2]*b1d + x[2]*x[2]*c1d);
00291 
00292         }
00293 
00294   
00295   
00297         int localindex (int comp) const
00298         {
00299           return number;
00300         }
00301 
00303         int codim () const
00304         {
00305           return codim_; 
00306         }
00307 
00309         int entity () const
00310         {
00311           return entity_;
00312         }
00313 
00315         int entityindex () const
00316         {
00317           return 0;
00318         }
00319 
00321         const FieldVector<CoordType,dim>& position () const
00322         {
00323           return pos;
00324         }
00325 
00326   private:
00327       int number,coeff,entity_,codim_;
00328       ResultType a[2], b[2], c[2],aa[2][2], bb[2][2];
00329 
00330       // Coefficients for the 1d-P2 shape functions
00331       ResultType a1d, b1d, c1d;
00332       FieldVector<CoordType,3> pos;
00333   
00334         
00335   };
00336 
00337 
00338   template<typename C, typename T, typename S>
00339   class P2HierarchicalPrismShapeFunctionSet
00340   {
00341   public:
00342 
00343         // compile time sizes
00344         enum { dim=3 };    // maps from R^d
00345         enum { comps=1 };      // to R^1
00346 
00347         enum { m=18 }; // total number of basis functions
00348 
00349         // export types
00350         typedef C CoordType;
00351         typedef T ResultType;
00352         typedef S value_type;
00353         typedef typename S::ImplementationType Imp; // Imp is either S or derived from S
00354 
00356         P2HierarchicalPrismShapeFunctionSet ()
00357         {
00358             for (int i=0; i<m; i++)
00359                 sf[i] = Imp(i); // assignment of derived class objects defined in wrapper
00360 
00361         }
00362 
00364         int size () const
00365         {
00366           return m;
00367         }
00368 
00370         int size (int entity, int codim) const
00371         {
00372           if (codim==dim) 
00373               return 6; 
00374           else if (codim==2)
00375               return 9;
00376           else if (codim==1)
00377               return 3;
00378 
00379           return 0;
00380         }
00381         
00383         const value_type& operator[] (int i) const
00384         {
00385           return sf[i]; // ok derived class reference goes for base class reference
00386         }
00387 
00389         int order () const
00390         {
00391           return 2;
00392         }
00393 
00395         GeometryType type () const
00396         {
00397       static GeometryType prism(GeometryType::prism, dim);
00398           return prism;
00399         }
00400 
00401   private:
00402         S sf[m];
00403   };
00404 
00406   template<typename C, typename T>
00407   class P2HierarchicalPrismShapeFunctionSetContainer
00408   {
00409   public:
00410         // compile time sizes
00411         enum { dim=3 };       
00412         enum { comps=1 };
00413         enum { maxsize=18 };
00414 
00415         // exported types
00416         typedef C CoordType; 
00417         typedef T ResultType;
00418         typedef P2HierarchicalPrismShapeFunctionSet<C,T,P2HierarchicalPrismShapeFunction<C,T> > value_type;
00419 
00420         const value_type& operator() (GeometryType type, int order) const
00421         {
00422           if(type.isPrism()) return p2prism;
00423           DUNE_THROW(NotImplemented, "type not implemented yet");
00424         }
00425 
00426   private:
00427         value_type p2prism;
00428   };
00429 
00432 }
00433 #endif