p2function.hh

00001 #ifndef DUNE_P2_FUNCTION_HH
00002 #define DUNE_P2_FUNCTION_HH
00003 
00004 #include <dune/istl/bvector.hh>
00005 #include <dune/disc/functions/functions.hh>
00006 #include <dune/disc/functions/globalp2index.hh>
00007 #include <dune/disc/shapefunctions/lagrangeshapefunctions.hh>
00008 
00015 namespace Dune
00016 {
00017 
00031   template<class GridType, class RT, typename IS, int m=1>
00032   class P2Function : virtual public ElementwiseCInfinityFunction<GridType,RT,m>,
00033                      virtual public H1Function<typename GridType::ctype,RT,GridType::dimension,m>,
00034                      virtual public C0GridFunction<GridType,RT,m>
00035   {
00037         typedef typename GridType::ctype DT;
00038 
00040         enum {dim=GridType::dimension};
00041 
00043         typedef typename GridType::template Codim<0>::Entity Entity;
00044 
00045   public:
00046       typedef Dune::BlockVector<Dune::FieldVector<RT,m> > RepresentationType;
00047 
00049       P2Function (const GridType& grid,  const IS& indexset)
00050           : grid_(grid), indexSet_(indexset), mapper_(grid,indexset)
00051         {
00052             data_.resize(mapper_.size());
00053         }
00054 
00056         /* A DifferentiableFunction can say how many derivatives exist
00057           and can be safely evaluated.
00058           \todo This is taken from the P1Function, but I find it strange.
00059           P2Functions are locally C^inf but globally only C^0!
00060          */
00061         virtual int order () const
00062         {
00063           return 2;
00064         }
00065 
00067 
00074         virtual RT eval (int comp, const Dune::FieldVector<DT,dim>& x) const
00075         {
00076           DUNE_THROW(NotImplemented, "global eval not implemented yet");
00077           return 0;
00078         }
00079 
00081 
00086         virtual void evalall (const Dune::FieldVector<DT,dim>& x, Dune::FieldVector<RT,m>& y) const
00087         {
00088           DUNE_THROW(NotImplemented, "global eval not implemented yet");
00089         }
00090 
00092 
00099         virtual RT derivative (int comp, const Dune::FieldVector<int,dim>& d, const Dune::FieldVector<DT,dim>& x) const
00100         {
00101           DUNE_THROW(NotImplemented, "global derivative not implemented yet");
00102         }
00103 
00105 
00111         virtual RT evallocal (int comp, const Entity& e, const Dune::FieldVector<DT,dim>& xi) const
00112         {
00113           RT value=0;
00114           Dune::GeometryType gt = e.type(); // extract type of element
00115           for (int i=0; i<Dune::LagrangeShapeFunctions<DT,RT,dim>::general(gt,2).size(); ++i) {
00116               int index = GlobalP2Index<GridType,dim>::map(e, 
00117                                                            mapper_, 
00118                                                            &Dune::LagrangeShapeFunctions<DT,RT,dim>::general(gt,2)[i]);
00119                 value += Dune::LagrangeShapeFunctions<DT,RT,dim>::general(gt,2)[i].evaluateFunction(0,xi)*(data_)[index][comp];
00120           }
00121           return value;
00122         }
00123 
00125 
00130       virtual void evalalllocal (const Entity& e, const Dune::FieldVector<DT,dim>& xi, 
00131                                  Dune::FieldVector<RT,m>& y) const {
00132 
00133           Dune::GeometryType gt = e.type(); // extract type of element
00134           y = 0;
00135           for (int i=0; i<Dune::LagrangeShapeFunctions<DT,RT,dim>::general(gt,2).size(); ++i) {
00136 
00137               RT basefuncvalue = Dune::LagrangeShapeFunctions<DT,RT,dim>::general(gt,2)[i].evaluateFunction(0,xi);
00138               int index = GlobalP2Index<GridType,dim>::map(e, 
00139                                                            mapper_, 
00140                                                            &Dune::LagrangeShapeFunctions<DT,RT,dim>::general(gt,2)[i]);
00141 
00142               for (int c=0; c<m; c++)
00143                   y[c] += basefuncvalue * data_[index][c];
00144           }
00145       }
00146 
00148 
00156         virtual RT derivativelocal (int comp, const Dune::FieldVector<int,dim>& d, 
00157                                     const Entity& e, const Dune::FieldVector<DT,dim>& xi) const
00158         {
00159             DUNE_THROW(NotImplemented, "global derivative not implemented yet");
00160         }
00161 
00163 
00167       const RepresentationType& operator* () const
00168         {
00169           return data_;
00170         }
00171 
00173 
00177         RepresentationType& operator* ()
00178         {
00179           return data_;
00180         }
00181 
00182   private:
00183       const GridType& grid_;
00184 
00185       const IS& indexSet_;
00186 
00187       RepresentationType data_;
00188 
00189       Dune::MultipleCodimMultipleGeomTypeMapper<GridType, IS, P2Layout> mapper_;
00190   };
00191 
00198   template<class G, class RT, int m=1>
00199   class LeafP2Function : public P2Function<G,RT,typename G::template Codim<0>::LeafIndexSet,m>
00200   {
00201   public:
00204         LeafP2Function (const G& grid)
00205           : P2Function<G,RT,typename G::template Codim<0>::LeafIndexSet,m>(grid,grid.leafIndexSet())
00206         {}
00207   };
00208 
00209 
00216   template<class G, class RT, int m=1>
00217   class LevelP2Function : public P2Function<G,RT,typename G::template Codim<0>::LevelIndexSet,m>
00218   {
00219   public:
00223         LevelP2Function (const G& grid, int level, bool extendoverlap=false)
00224           : P2Function<G,RT,typename G::template Codim<0>::LevelIndexSet,m>(grid,grid.levelIndexSet(level))
00225         {}
00226   };
00227 
00228 }  // end namespace Dune
00229 #endif