hierarchicalsimplexshapefunctions.hh

00001 #ifndef DUNE_HIERARCHICAL_SIMPLEX_SHAPEFUNCTIONS_HH
00002 #define DUNE_HIERARCHICAL_SIMPLEX_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 
00009 #include"hierarchicalcubeshapefunctions.hh"
00010 
00017 namespace Dune
00018 {
00026   /***********************************************************
00027    * The interface for  Simplex shape functions of arbitrary
00028    * order and element type 
00029    * now support only triangle and tetrahedron
00030    ***********************************************************/
00031 
00032 
00033     /***********************************************************
00034      * P2 shape functions for the simplex 
00035      ***********************************************************/
00041     /*  specialization for 1D 
00042         polynomials of the form coeff[0]*x^2 + coeff[1]*x + coeff[2]
00043     */
00044   
00045 
00046 template<typename C, typename T, int dim>
00047 class P2HierarchicalSimplexShapeFunction;
00048 
00049 template<typename C, typename T>
00050 class P2HierarchicalSimplexShapeFunction<C,T,1>
00051   : public P2HierarchicalCubeShapeFunction<C,T,1>
00052 {
00053 public:
00054   enum {dim=1};
00055   enum {comps=1};
00056   enum {m=3}; 
00057   typedef C CoordType;
00058   typedef T ResultType;
00059   typedef P2HierarchicalSimplexShapeFunction ImplementationType;
00060   P2HierarchicalSimplexShapeFunction(int i,int en,int co)
00061     : P2HierarchicalCubeShapeFunction<C,T,1>(i,en,co,FieldVector<int,1>(i))
00062     {}
00063 
00064   P2HierarchicalSimplexShapeFunction()
00065     {}
00066 };
00067 
00068 //specialization 1D
00069 //-------------------------------------------------
00070 template<typename C, typename T, int d, typename S>
00071 class P2HierarchicalSimplexShapeFunctionSet;
00072 
00073 template<typename C, typename T,typename S>
00074 class P2HierarchicalSimplexShapeFunctionSet<C,T,1,S>
00075 
00076 {
00077 
00078 public:
00079   enum {dim=1};
00080   enum {comps=1};
00081   enum {m=3}; 
00082 
00083   typedef C CoordType;
00084   typedef T ResultType;
00085   typedef S value_type;
00086   typedef typename S::ImplementationType Imp; 
00088   P2HierarchicalSimplexShapeFunctionSet()
00089   {
00090     ReferenceSimplex<C,1> simpline;
00091     int j=0;
00092  for (int c=2; c>=0; --c)
00093    for(int e=0;e<simpline.size(c);++e)
00094      {
00095       sf[j] = Imp(j,e,c);
00096       j++;
00097      }
00098   }
00099   
00101   int size () const
00102   {
00103     return m;
00104   }
00106   int size (int entity, int codim) const
00107   {
00108     return 1;
00109   }
00111   const value_type& operator[] (int i) const
00112   {
00113     return sf[i]; // ok derived class reference goes for base class reference
00114   }
00115 
00117   int order () const
00118   {
00119     return 2;
00120   }
00121 
00123   GeometryType type () const
00124   {
00125     static GeometryType simplex(GeometryType::simplex, dim);
00126     return simplex;
00127   }
00128 private:
00129   S sf[m];
00130 };
00131 
00133  // specialization for line
00134   template<typename C, typename T, int d>
00135   class P2HierarchicalSimplexShapeFunctionSetContainer;
00136 
00137 template<typename C,typename T>
00138 class P2HierarchicalSimplexShapeFunctionSetContainer<C,T,1>
00139   {
00140 
00141   public:
00142 
00143     enum {dim=1};
00144     enum{comps=1};
00145     enum {maxsize=3};
00146 
00147     // exported types
00148     typedef C CoordType; 
00149     typedef T ResultType;
00150     typedef P2HierarchicalSimplexShapeFunctionSet<C,T,dim,P2HierarchicalSimplexShapeFunction<C,T,dim> > value_type;
00151 
00152     const value_type& operator() (GeometryType type, int order) const
00153     {
00154       if (type.isSimplex()) return p2simplex;
00155       DUNE_THROW(NotImplemented,"type not yet implemented");
00156         }
00157   private:
00158     value_type p2simplex;
00159   };
00160 
00161 /*  specialization for 2D 
00162     A class for piecewise quadratic shape functions in a simplex geometry for triangles
00163  */
00164   /*
00165 polynomial of the form: coeff * ( a[0] + b[0]*x + b[1]*y )(a[1] + c[0]*x + c[1]*y )
00166 derrivative of the form:( aa[0][dir] + bb[0][dir]*x + bb[1][dir]*y)
00167 dir=0 for x, dir=1 for y
00168        */
00169 
00170 template<typename C, typename T, int dim>
00171 class P2HierarchicalSimplexShapeFunction;
00172 
00173 template<typename C, typename T>
00174 class P2HierarchicalSimplexShapeFunction<C,T,2>
00175 {
00176 
00177 public:
00178 
00179   enum {dim=2};
00180   enum{comps=1};
00181   enum {m=6}; // no of basis functions=6 
00182 
00183   typedef C CoordType;
00184   typedef T ResultType;
00185   typedef P2HierarchicalSimplexShapeFunction ImplementationType;
00186 
00187   P2HierarchicalSimplexShapeFunction(int i,int en,int co)
00188     {
00189       number=i;
00190       ent = en;
00191       cod = co;
00192     
00193       switch(i)
00194         {
00195         case 0:
00196           //--interpolation point associated with shape fn
00197           pos[0]=0.0;
00198           pos[1]=0.0;
00199           //--
00200           coeff=1;
00201           a[0]= 1.0;
00202           a[1]= 1.0;
00203           b[0]=-1.0;
00204           b[1]=-1.0;
00205           c[0]= 0.0;
00206           c[1]= 0.0;
00207           //x derivative
00208           aa[0][0]=-1;
00209           bb[0][0]=0;
00210           bb[1][0]=0;
00211           //y derivative
00212           aa[0][1]=-1;
00213           bb[0][1]=0;
00214           bb[1][1]=0;
00215           break;
00216         case 1:
00217           //--interpolation point associated with shape fn
00218           pos[0]=1.0;
00219           pos[1]=0.0;
00220           //--
00221           coeff=1;
00222           a[0]=0.0;
00223           a[1]=1;
00224           b[0]=1.0;
00225           b[1]=0.0;
00226           c[0]=0.0;
00227           c[1]=0.0;
00228           //x derivative
00229           aa[0][0]=1;
00230           bb[0][0]=0;
00231           bb[1][0]=0;
00232           //y derivative
00233           aa[0][1]=0;
00234           bb[0][1]=0;
00235           bb[1][1]=0;
00236           break;
00237         case 2:
00238           //--interpolation point associated with shape fn
00239           pos[0]=0.0;
00240           pos[1]=1.0;
00241           //--
00242           coeff=1;
00243           a[0]=0.0;
00244           a[1]=1.0;
00245           b[0]=0.0;
00246           b[1]=1.0;
00247           c[0]=0.0;
00248           c[1]=0.0;
00249           //x derivative
00250           aa[0][0]=0;
00251           bb[0][0]=0;
00252           bb[1][0]=0;
00253           //y derivative
00254           aa[0][1]=1;
00255           bb[0][1]=0;
00256           bb[1][1]=0;
00257           break;
00258         case 3:
00259           //--interpolation point associated with shape fn
00260           pos[0]=0.5;
00261           pos[1]=0.5;
00262           //--
00263           coeff=4;
00264           a[0]=0.0;
00265           a[1]=0.0;
00266           b[0]=1.0;
00267           b[1]=0.0;
00268           c[0]=0.0;
00269           c[1]=1.0;
00270          //x derivative
00271           aa[0][0]=0;
00272           bb[0][0]=0;
00273           bb[1][0]=4;
00274           //y derivative
00275           aa[0][1]=0;
00276           bb[0][1]=4;
00277           bb[1][1]=0;
00278           break;
00279         case 4:
00280           //--interpolation point associated with shape fn
00281           pos[0]=0.0;
00282           pos[1]=0.5;
00283           //--
00284           coeff=4;
00285           a[0]=0.0;
00286           a[1]=1.0;
00287           b[0]=0.0;
00288           b[1]=1.0;
00289           c[0]=-1.0;
00290           c[1]=-1.0;
00291           //x derivative
00292           aa[0][0]=0;
00293           bb[0][0]=0;
00294           bb[1][0]=-4;
00295           //y derivative
00296           aa[0][1]=4;
00297           bb[0][1]=-4;
00298           bb[1][1]=-8;
00299           break;
00300         case 5:
00301           //--interpolation point associated with shape fn
00302           pos[0]=0.5;
00303           pos[1]=0.0;
00304           //--
00305           coeff=4;
00306           a[0]=0.0;
00307           a[1]=1.0;
00308           b[0]=1.0;
00309           b[1]=0.0;
00310           c[0]=-1.0;
00311           c[1]=-1.0;
00312           //x derivative
00313           aa[0][0]=4;
00314           bb[0][0]=-8;
00315           bb[1][0]=-4;
00316           //y derivative
00317           aa[0][1]=0;
00318           bb[0][1]=-4;
00319           bb[1][1]=0;
00320           break;
00321         default:
00322           DUNE_THROW(RangeError, "wrong no of shape fns? check simplexshapefunctins.hh");
00323           break;
00324         }
00325     }
00326 
00327  P2HierarchicalSimplexShapeFunction()
00328   {}
00329 
00331   ResultType evaluateFunction (int comp, const FieldVector<CoordType,2>& x) const
00332   {
00333    
00334     ResultType phi1=a[0];
00335     ResultType phi2=a[1];
00336 
00337     for (int j=0;j<2;++j)
00338       {
00339         phi1 +=b[j]*x[j];
00340         phi2 +=c[j]*x[j];
00341       }
00342     ResultType phi=coeff*phi1*phi2;
00343     return phi;
00344   } 
00345   
00346 
00347   ResultType evaluateDerivative (int comp, int dir, const FieldVector<CoordType,2>& x) const
00348   
00349   {
00350 
00351     ResultType deriv=aa[0][dir];
00352    
00353     for (int j=0;j<2;++j)
00354       {
00355         deriv +=bb[j][dir]*x[j];
00356       }
00357         
00358     return deriv;
00359   }
00360 
00361   
00362   
00364   int localindex (int comp) const
00365   {
00366     return number;
00367   }
00368 
00370   int codim () const
00371   {
00372     return cod;
00373   }
00374 
00376   int entity () const
00377   {
00378     return ent;
00379   }
00380 
00382   int entityindex () const
00383   {
00384     return 0;
00385   }
00386 
00388   const FieldVector<CoordType,2>& position () const
00389   {
00390     return pos;
00391   }
00392 
00393 private:
00394   int number,coeff,ent,cod;
00395   ResultType a[2], b[2], c[2],aa[2][2], bb[2][2];
00396   FieldVector<CoordType,2> pos;
00397   
00398 };
00399 
00400 
00401 //specialization 2D
00402 //-------------------------------------------------
00403 template<typename C, typename T, int d, typename S>
00404 class P2HierarchicalSimplexShapeFunctionSet;
00405 
00406 template<typename C, typename T,typename S>
00407 class P2HierarchicalSimplexShapeFunctionSet<C,T,2,S>
00408 
00409 {
00410 
00411 public:
00412   enum {dim=2};
00413   enum {comps=1};
00414   enum {m=6}; 
00415 
00416   typedef C CoordType;
00417   typedef T ResultType;
00418   typedef S value_type;
00419   typedef typename S::ImplementationType Imp; // Imp is either S or derived from S
00421   P2HierarchicalSimplexShapeFunctionSet()
00422   {
00423     ReferenceSimplex<C,2> triang;
00424     int j=0;
00425  for (int c=2; c>=1; --c)
00426    for(int e=0;e<triang.size(c);++e)
00427      {
00428       sf[j] = Imp(j,e,c);
00429       j++;
00430      }
00431   }
00432  
00433   
00435   int size () const
00436   {
00437     return m;
00438   }
00440   int size (int entity, int codim) const
00441   {
00442     return 1;
00443   }
00445   const value_type& operator[] (int i) const
00446   {
00447     return sf[i]; // ok derived class reference goes for base class reference
00448   }
00449 
00451   int order () const
00452   {
00453     return 2;
00454   }
00455 
00457   GeometryType type () const
00458   {
00459     static GeometryType simplex(GeometryType::simplex, dim);
00460     return simplex;
00461   }
00462 private:
00463   S sf[m];
00464 };
00465 
00467  // specialization for triangle
00468   template<typename C, typename T, int d>
00469   class P2HierarchicalSimplexShapeFunctionSetContainer;
00470 
00471 template<typename C,typename T>
00472 class P2HierarchicalSimplexShapeFunctionSetContainer<C,T,2>
00473   {
00474 
00475   public:
00476 
00477     enum {dim=2};
00478     enum{comps=1};
00479     enum {maxsize=6};
00480 
00481     // exported types
00482     typedef C CoordType; 
00483     typedef T ResultType;
00484     typedef P2HierarchicalSimplexShapeFunctionSet<C,T,dim,P2HierarchicalSimplexShapeFunction<C,T,dim> > value_type;
00485 
00486     const value_type& operator() (GeometryType type, int order) const
00487     {
00488       if (type.isSimplex()) return p2simplex;
00489       DUNE_THROW(NotImplemented,"type not yet implemented");
00490         }
00491   private:
00492     value_type p2simplex;
00493   };
00494 
00495 /*  specialization for 3D 
00496     A class for piecewise quadratic shape functions in a simplex geometry for tetrahedron
00497     Evaluation is done at the quadrature points
00498  */
00499 template<typename C, typename T, int dim>
00500 class P2HierarchicalSimplexShapeFunction;
00501 
00502  template<typename C, typename T>
00503  class P2HierarchicalSimplexShapeFunction<C,T,3>
00504 {
00505 public:
00506   enum {dim=3};
00507   enum{comps=1};
00508   enum {m=10}; // no of basis functions=10  
00509 
00510   typedef C CoordType;
00511   typedef T ResultType;
00512   typedef P2HierarchicalSimplexShapeFunction ImplementationType;
00513 
00514 
00515 P2HierarchicalSimplexShapeFunction(int i,int en,int co)
00516     {
00517       number=i;
00518       ent = en;
00519       cod = co;
00520     
00521       switch(i)
00522         {
00523         case 0:
00524           //--interpolation point associated with shape fn
00525           pos[0]=0.0;
00526           pos[1]=0.0;
00527           pos[2]=0.0;
00528           //--
00529           coeff=1;
00530           a[0]=1.0;
00531           a[1]=1.0;
00532           b[0]=-1.0;
00533           b[1]=-1.0;
00534           b[2]=-1.0;
00535           c[0]= 0.0;
00536           c[1]= 0.0;
00537           c[2]= 0.0;
00538           //x derivative
00539           aa[0][0]=-1.0;
00540           bb[0][0]=0.0;
00541           bb[1][0]=0.0;
00542           bb[2][0]=0.0;
00543           //y derivative
00544           aa[0][1]=-1.0;
00545           bb[0][1]=0.0;
00546           bb[1][1]=0.0;
00547           bb[2][1]=0.0;
00548           // z derivative
00549           aa[0][2]=-1.0;
00550           bb[0][2]=0.0;
00551           bb[1][2]=0.0;
00552           bb[2][2]=0.0;
00553           break;
00554         case 1:
00555           //--interpolation point associated with shape fn
00556           pos[0]=1.0;
00557           pos[1]=0.0;
00558           pos[2]=0.0;
00559           //--
00560           coeff=1;
00561           a[0]=0.0;
00562           a[1]=1.0;
00563           b[0]=1.0;
00564           b[1]=0.0;
00565           b[2]=0.0;
00566           c[0]=0.0;
00567           c[1]=0.0;
00568           c[2]=0.0;
00569            //x derivative
00570           aa[0][0]=1.0;
00571           bb[0][0]=0.0;
00572           bb[1][0]=0.0;
00573           bb[2][0]=0.0;
00574           //y derivative
00575           aa[0][1]=0.0;
00576           bb[0][1]=0.0;
00577           bb[1][1]=0.0;
00578           bb[2][1]=0.0;
00579           // z derivative
00580           aa[0][2]=0.0;
00581           bb[0][2]=0.0;
00582           bb[1][2]=0.0;
00583           bb[2][2]=0.0;
00584           break;
00585         case 2:
00586           //--interpolation point associated with shape fn
00587           pos[0]=0.0;
00588           pos[1]=1.0;
00589           pos[2]=0.0;
00590           //--
00591           coeff=1;
00592           a[0]=0.0;
00593           a[1]=1.0;
00594           b[0]=0.0;
00595           b[1]=1.0;
00596           b[2]=0.0;
00597           c[0]=0.0;
00598           c[1]=0.0;
00599           c[2]=0.0;
00600           //x derivative
00601           aa[0][0]=0.0;
00602           bb[0][0]=0.0;
00603           bb[1][0]=0.0;
00604           bb[2][0]=0.0;
00605           //y derivative
00606           aa[0][1]=1.0;
00607           bb[0][1]=0.0;
00608           bb[1][1]=0.0;
00609           bb[2][1]=0.0;
00610           // z derivative
00611           aa[0][2]=0.0;
00612           bb[0][2]=0.0;
00613           bb[1][2]=0.0;
00614           bb[2][2]=0.0;
00615           break;
00616         case 3:
00617           //--interpolation point associated with shape fn
00618           pos[0]=0.0;
00619           pos[1]=0.0;
00620           pos[2]=1.0;
00621           //--
00622           coeff=1;
00623           a[0]=0.0;
00624           a[1]=1.0;
00625           b[0]=0.0;
00626           b[1]=0.0;
00627           b[2]=1.0;
00628           c[0]=0.0;
00629           c[1]=0.0;
00630           c[2]=0.0;
00631           //x derivative
00632           aa[0][0]=0.0;
00633           bb[0][0]=0.0;
00634           bb[1][0]=0.0;
00635           bb[2][0]=0.0;
00636           //y derivative
00637           aa[0][1]=0.0;
00638           bb[0][1]=0.0;
00639           bb[1][1]=0.0;
00640           bb[2][1]=0.0;
00641           // z derivative
00642           aa[0][2]=1.0;
00643           bb[0][2]=0.0;
00644           bb[1][2]=0.0;
00645           bb[2][2]=0.0;
00646           break;
00647         case 4:
00648           //--interpolation point associated with shape fn
00649           pos[0]=0.5;
00650           pos[1]=0.0;
00651           pos[2]=0.0;
00652           //--
00653           coeff=4;
00654           a[0]=0.0;
00655           a[1]=1.0;
00656           b[0]=1.0;
00657           b[1]=0.0;
00658           b[2]=0.0;
00659           c[0]=-1.0;
00660           c[1]=-1.0;
00661           c[2]=-1.0;
00662           //x derivative
00663           aa[0][0]=4.0;
00664           bb[0][0]=-8.0;
00665           bb[1][0]=-4.0;
00666           bb[2][0]=-4.0;
00667           //y derivative
00668           aa[0][1]=0.0;
00669           bb[0][1]=-4.0;
00670           bb[1][1]=0.0;
00671           bb[2][1]=0.0;
00672           // z derivative
00673           aa[0][2]=0.0;
00674           bb[0][2]=-4.0;
00675           bb[1][2]=0.0;
00676           bb[2][2]=0.0;
00677           break;
00678         case 5:
00679           //--interpolation point associated with shape fn
00680           pos[0]=0.5;
00681           pos[1]=0.5;
00682           pos[2]=0.0;
00683           //--
00684           coeff=4;
00685           a[0]=0.0;
00686           a[1]=0.0;
00687           b[0]=1.0;
00688           b[1]=0.0;
00689           b[2]=0.0;
00690           c[0]=0.0;
00691           c[1]=1.0;
00692           c[2]=0.0;
00693           //x derivative
00694           aa[0][0]=0.0;
00695           bb[0][0]=0.0;
00696           bb[1][0]=4.0;
00697           bb[2][0]=0.0;
00698           //y derivative
00699           aa[0][1]=0.0;
00700           bb[0][1]=4.0;
00701           bb[1][1]=0.0;
00702           bb[2][1]=0.0;
00703           // z derivative
00704           aa[0][2]=0.0;
00705           bb[0][2]=0.0;
00706           bb[1][2]=0.0;
00707           bb[2][2]=0.0;
00708           break;
00709         case 6:
00710           //--interpolation point associated with shape fn
00711           pos[0]=0.0;
00712           pos[1]=0.5;
00713           pos[2]=0.0;
00714           //--
00715           coeff=4;
00716           a[0]=0.0;
00717           a[1]=1.0;
00718           b[0]=0.0;
00719           b[1]=1.0;
00720           b[2]=0.0;
00721           c[0]=-1.0;
00722           c[1]=-1.0;
00723           c[2]=-1.0;
00724          //x derivative
00725           aa[0][0]=0.0;
00726           bb[0][0]=0.0;
00727           bb[1][0]=-4.0;
00728           bb[2][0]=0.0;
00729           //y derivative
00730           aa[0][1]=4.0;
00731           bb[0][1]=-4.0;
00732           bb[1][1]=-8.0;
00733           bb[2][1]=-4.0;
00734           // z derivative
00735           aa[0][2]=0.0;
00736           bb[0][2]=0.0;
00737           bb[1][2]=-4.0;
00738           bb[2][2]=0.0;
00739           break;
00740         case 7:
00741           //--interpolation point associated with shape fn
00742           pos[0]=0.0;
00743           pos[1]=0.0;
00744           pos[2]=0.5;
00745           //--
00746           coeff=4;
00747           a[0]=0.0;
00748           a[1]=1.0;
00749           b[0]=0.0;
00750           b[1]=0.0;
00751           b[2]=1.0;
00752           c[0]=-1.0;
00753           c[1]=-1.0;
00754           c[2]=-1.0;
00755           //x derivative
00756           aa[0][0]=0.0;
00757           bb[0][0]=0.0;
00758           bb[1][0]=0.0;
00759           bb[2][0]=-4.0;
00760           //y derivative
00761           aa[0][1]=0.0;
00762           bb[0][1]=0.0;
00763           bb[1][1]=0.0;
00764           bb[2][1]=-4.0;
00765           // z derivative
00766           aa[0][2]=4.0;
00767           bb[0][2]=-4.0;
00768           bb[1][2]=-4.0;
00769           bb[2][2]=-8.0;
00770           break;
00771         case 8:
00772           //--interpolation point associated with shape fn
00773           pos[0]=0.5;
00774           pos[1]=0.0;
00775           pos[2]=0.5;
00776           //--
00777           coeff=4;
00778           a[0]=0.0;
00779           a[1]=0.0;
00780           b[0]=1.0;
00781           b[1]=0.0;
00782           b[2]=0.0;
00783           c[0]=0.0;
00784           c[1]=0.0;
00785           c[2]=1.0;
00786          //x derivative
00787           aa[0][0]=0.0;
00788           bb[0][0]=0.0;
00789           bb[1][0]=0.0;
00790           bb[2][0]=4.0;
00791           //y derivative
00792           aa[0][1]=0.0;
00793           bb[0][1]=0.0;
00794           bb[1][1]=0.0;
00795           bb[2][1]=0.0;
00796           // z derivative
00797           aa[0][2]=0.0;
00798           bb[0][2]=4.0;
00799           bb[1][2]=0.0;
00800           bb[2][2]=0.0;
00801           break;
00802         case 9:
00803           //--interpolation point associated with shape fn
00804           pos[0]=0.0;
00805           pos[1]=0.5;
00806           pos[2]=0.5;
00807           //--
00808           coeff=4;
00809           a[0]=0.0;
00810           a[1]=0.0;
00811           b[0]=0.0;
00812           b[1]=1.0;
00813           b[2]=0.0;
00814           c[0]=0.0;
00815           c[1]=0.0;
00816           c[2]=1.0;
00817          //x derivative
00818           aa[0][0]=0.0;
00819           bb[0][0]=0.0;
00820           bb[1][0]=0.0;
00821           bb[2][0]=0.0;
00822           //y derivative
00823           aa[0][1]=0.0;
00824           bb[0][1]=0.0;
00825           bb[1][1]=0.0;
00826           bb[2][1]=4.0;
00827           // z derivative
00828           aa[0][2]=0.0;
00829           bb[0][2]=0.0;
00830           bb[1][2]=4.0;
00831           bb[2][2]=0.0;
00832           break;
00833         default:
00834           DUNE_THROW(RangeError, "wrong no of shape fns? check out please");
00835           break;
00836         }
00837     }
00838 
00839 
00840 // default constructor
00841   P2HierarchicalSimplexShapeFunction()
00842     {
00843 
00844     }
00845 
00847   ResultType evaluateFunction (int comp, const FieldVector<CoordType,3>& x) const
00848   {
00849    
00850     ResultType phi1=a[0];
00851     ResultType phi2=a[1];
00852 
00853     for (int j=0;j<3;++j)
00854       {
00855         phi1 +=b[j]*x[j];
00856         phi2 +=c[j]*x[j];
00857       }
00858     ResultType phi=coeff*phi1*phi2;
00859     return phi;
00860   } 
00861   
00862 
00863   ResultType evaluateDerivative (int comp, int dir, const FieldVector<CoordType,3>& x) const
00864   
00865   {
00866 
00867     ResultType deriv=aa[0][dir];
00868    
00869     for (int j=0;j<3;++j)
00870       {
00871         deriv +=bb[j][dir]*x[j];
00872       }
00873         
00874     return deriv;
00875   
00876   }
00877 
00878    
00880   int localindex (int comp) const
00881   {
00882     return number;
00883   }
00884 
00886   int codim () const
00887   {
00888     return cod;
00889   }
00890 
00892   int entity () const
00893   {
00894     return ent;
00895   }
00896 
00898   int entityindex () const
00899   {
00900     return 0;
00901   }
00902 
00904   const FieldVector<CoordType,dim>& position () const
00905   {
00906     return pos;
00907   }
00908 
00909 private:
00910 
00911  int number,coeff,ent,cod;
00912  ResultType a[2], b[3], c[3],aa[3][3], bb[3][3];
00913   FieldVector<CoordType,3> pos;
00914   
00915         
00916 };
00917 
00918 
00919 //specialization 3D
00920 //-----------------------
00921 template<typename C, typename T, int d, typename S>
00922 class P2HierarchicalSimplexShapeFunctionSet;
00923 
00924 template<typename C, typename T,typename S>
00925 class P2HierarchicalSimplexShapeFunctionSet<C,T,3,S>
00926 
00927 {
00928 
00929 public:
00930   enum {dim=3};
00931   enum {comps=1};
00932   enum {m=10}; 
00933 
00934   typedef C CoordType;
00935   typedef T ResultType;
00936   typedef S value_type;
00937   typedef typename S::ImplementationType Imp; // Imp is either S or derived from S
00939   P2HierarchicalSimplexShapeFunctionSet()
00940   {
00941     ReferenceSimplex<C,3> tetrahed;
00942     int j=0;
00943  for (int c=3; c>=2; --c)
00944    for(int e=0;e<tetrahed.size(c);++e)
00945      {
00946       sf[j] = Imp(j,e,c);
00947       j++;
00948      }
00949   }
00950  
00951   
00953   int size () const
00954   {
00955     return m;
00956   }
00958   int size (int entity, int codim) const
00959   {
00960     return 1;
00961   }
00963   const value_type& operator[] (int i) const
00964   {
00965     return sf[i]; // ok derived class reference goes for base class reference
00966   }
00967 
00969   int order () const
00970   {
00971     return 2;
00972   }
00973 
00975   GeometryType type () const
00976   {
00977     static GeometryType simplex(GeometryType::simplex, dim);
00978     return simplex;
00979   }
00980 private:
00981   S sf[m];
00982 };
00983 
00984 
00986 
00987 // specialization for tetrahedron
00988   template<typename C, typename T, int d>
00989   class P2HierarchicalSimplexShapeFunctionSetContainer;
00990 
00991 template<typename C,typename T>
00992 class P2HierarchicalSimplexShapeFunctionSetContainer<C,T,3>
00993   {
00994   public:
00995 
00996     enum {dim=3};
00997     enum{comps=1};
00998     enum {maxsize=10};
00999 
01000     // exported types
01001     typedef C CoordType; 
01002     typedef T ResultType;
01003     typedef P2HierarchicalSimplexShapeFunctionSet<C,T,dim,P2HierarchicalSimplexShapeFunction<C,T,dim> > value_type;
01004 
01005     const value_type& operator() (GeometryType type, int order) const
01006     {
01007       if (type.isSimplex()) return p2simplex;
01008       DUNE_THROW(NotImplemented,"type not yet implemented");
01009         }
01010   private:
01011     value_type p2simplex;
01012   };
01013 
01014 }
01015 #endif