p0function.hh

00001 // $Id: p0function.hh 520 2008-10-07 22:06:21Z sander $
00002 
00003 #ifndef DUNE_P0FUNCTION_HH
00004 #define DUNE_P0FUNCTION_HH
00005 
00006 //C++ includes
00007 #include<new>
00008 #include<iostream>
00009 #include<vector>
00010 #include<list>
00011 #include <map>
00012 
00013 // Dune includes
00014 #include<dune/common/fvector.hh>
00015 #include<dune/common/exceptions.hh>
00016 #include<dune/grid/common/grid.hh>
00017 #include<dune/grid/common/mcmgmapper.hh>
00018 #include<dune/grid/common/universalmapper.hh>
00019 #include<dune/grid/io/file/vtk/vtkwriter.hh>
00020 #include<dune/istl/bvector.hh>
00021 #include<dune/istl/operators.hh>
00022 #include<dune/istl/bcrsmatrix.hh>
00023 #include<dune/disc/shapefunctions/lagrangeshapefunctions.hh>
00024 
00025 // same directory includes
00026 #include"functions.hh"
00027 
00033 namespace Dune
00034 {
00044   // forward declaration
00045   //  template<class G, class RT> class P0FEFunctionManager;
00046 
00047 
00049 
00058   template<class GV, class RT, int m=1>
00059   class P0Function :
00060     virtual public GridFunction<typename GV::Grid,RT,m>,
00061     virtual public GridFunctionGlobalEvalDefault<GV,RT,m>,
00062     virtual public FunctionDefault<typename GV::Grid::ctype,RT,GV::Grid::dimension,m>,
00063     virtual public L2Function<typename GV::Grid::ctype,RT,GV::Grid::dimension,m>
00064   {
00066         typedef typename GV::Grid G;
00067 
00069         typedef typename G::ctype DT;
00070 
00072         enum {n=G::dimension};
00073 
00075         typedef typename G::template Codim<0>::Entity Entity;
00076 
00078         template<int dim>
00079         struct P0Layout
00080         {
00081           bool contains (Dune::GeometryType gt)
00082           {
00083               return gt.dim() == dim;
00084           }
00085         }; 
00086 
00088         P0Function (const P0Function&);
00089 
00090   public:
00091         typedef BlockVector<FieldVector<RT,m> > RepresentationType;
00092 
00094         P0Function (const GV& gv_)
00095           : GridFunctionGlobalEvalDefault<GV,RT,m>(gv_),
00096             gv(gv_), mapper_(gv_.grid(),gv_.indexSet())
00097         {
00098           oldcoeff = 0;
00099           try {
00100                 coeff = new RepresentationType(mapper_.size());
00101           }
00102           catch (std::bad_alloc) {
00103                 std::cerr << "not enough memory in P0Function" << std::endl;
00104                 throw; // rethrow exception
00105           }
00106 //        std::cout << "making  function with " << mapper_.size() << " components" << std::endl;
00107         }
00108 
00110         ~P0Function ()
00111         {
00112           delete coeff;
00113           if (oldcoeff!=0) delete oldcoeff;
00114         }
00115 
00117 
00123         virtual RT evallocal (int comp, const Entity& e, const Dune::FieldVector<DT,n>& xi) const
00124         {
00125           return (*coeff)[mapper_.map(e)][comp];
00126         }
00127 
00129 
00134         virtual void evalalllocal (const Entity& e, const Dune::FieldVector<DT,G::dimension>& xi, 
00135                                                            Dune::FieldVector<RT,m>& y) const
00136         {
00137           for (int c=0; c<m; c++)
00138                 y[c] = (*coeff)[mapper_.map(e)][c];
00139         }
00140 
00141 
00143 
00149         void interpolate (const C0GridFunction<G,RT,m>& u)
00150         {
00151           typedef typename GV::template Codim<0>::Iterator Iterator;
00152 
00153           Iterator eendit = gv.template end<0>();
00154           for (Iterator it = gv.template begin<0>(); it!=eendit; ++it)
00155                 {
00156                   Dune::GeometryType gt = it->type();
00157                   for (int c=0; c<m; c++)
00158                         (*coeff)[mapper_.map(*it)][c] = 
00159                           u.evallocal(c,*it,Dune::ReferenceElements<DT,n>::general(gt).position(0,0));
00160                 }
00161         }
00162 
00164 
00168         const RepresentationType& operator* () const
00169         {
00170           return (*coeff);
00171         }
00172 
00174 
00178         RepresentationType& operator* ()
00179         {
00180           return (*coeff);
00181         }
00182 
00187         void preAdapt ()
00188         {
00189           DUNE_THROW(NotImplemented, "preAdapt()");
00190         }
00191 
00200         void postAdapt ()
00201         {
00202           DUNE_THROW(NotImplemented, "postAdapt()");
00203         }
00204  
00207         const MultipleCodimMultipleGeomTypeMapper<G,typename GV::IndexSet,P0Layout>& mapper () const
00208         {
00209           return mapper_;
00210         }
00211 
00213         void vtkout (VTKWriter<GV>& vtkwriter, std::string s) const
00214         {
00215           typename VTKWriter<GV>::VTKFunction *p = new VTKGridFunctionWrapper<typename GV::IndexSet,RT,m>(*this,s);
00216           vtkwriter.addCellData(p);
00217         }
00218         
00219   private:
00220         // The GridView.  Don't use a reference here since Grid::levelView
00221         // returns a temporary
00222         const GV gv;
00223 
00224         // we need a mapper
00225         MultipleCodimMultipleGeomTypeMapper<G,typename GV::IndexSet,P0Layout> mapper_;
00226 
00227         // and a dynamically allocated vector
00228         RepresentationType* coeff;
00229 
00230         // saved pointer in update phase
00231         RepresentationType* oldcoeff;
00232   };
00233 
00234 
00237   template<class G, class RT, int m=1>
00238   class LeafP0Function : public P0Function<typename G::LeafGridView,RT,m>
00239   {
00240   public:
00241         LeafP0Function (const G& grid) 
00242           : GridFunctionGlobalEvalDefault<typename G::LeafGridView,RT,m>(grid.leafView())
00243           , P0Function<typename G::LeafGridView,RT,m>(grid.leafView())
00244         {}
00245   };
00246 
00247 
00250   template<class G, class RT, int m=1>
00251   class LevelP0Function : public P0Function<typename G::LevelGridView,RT,m>
00252   {
00253   public:
00254     LevelP0Function (const G& grid, int level) 
00255       : GridFunctionGlobalEvalDefault<typename G::LevelGridView,RT,m>(grid.levelView(level))
00256       , P0Function<typename G::LevelGridView,RT,m>(grid.levelView(level))
00257     {}
00258   };
00259 
00260 
00261 
00267   template<class G, class T>
00268   class LeafP0FunctionWrapper : virtual public GridFunction<G,typename T::value_type,1>,
00269                                            virtual public L2Function<typename G::ctype,typename T::value_type,G::dimension,1>
00270   {
00272         typedef typename G::ctype DT;
00273         typedef typename T::value_type RT;
00274         typedef typename G::template Codim<0>::LeafIndexSet IS;
00275         enum {m=1};
00276 
00278         enum {n=G::dimension};
00279 
00281         typedef typename G::template Codim<0>::Entity Entity;
00282 
00284         template<int dim>
00285         struct P0Layout
00286         {
00287           bool contains (Dune::GeometryType gt)
00288           {
00289               return gt.dim() == dim;
00290           }
00291         }; 
00292 
00294         LeafP0FunctionWrapper (const LeafP0FunctionWrapper&);
00295         LeafP0FunctionWrapper& operator= (const LeafP0FunctionWrapper&);
00296 
00297   public:
00298         typedef T RepresentationType;
00299 
00301         LeafP0FunctionWrapper (const G& g, const T& v) 
00302           : grid(g), is(g.leafIndexSet()), mapper(g,g.leafIndexSet()), coeff(v)
00303         {
00304           if ((unsigned int)mapper.size()!=v.size())
00305                 DUNE_THROW(MathError,"LeafP0FunctionWrapper: size of vector does not match grid size");
00306 
00307         }
00308 
00310 
00316         virtual RT eval (int comp, const Dune::FieldVector<DT,n>& x) const
00317         {
00318           DUNE_THROW(NotImplemented, "global eval not implemented yet");
00319           return 0;
00320         }
00321 
00323 
00327         virtual void evalall (const Dune::FieldVector<DT,n>& x, Dune::FieldVector<RT,m>& y) const
00328         {
00329           DUNE_THROW(NotImplemented, "global eval not implemented yet");
00330         }
00331 
00333 
00339         virtual RT evallocal (int comp, const Entity& e, const Dune::FieldVector<DT,n>& xi) const
00340         {
00341           return coeff[mapper.map(e)];
00342         }
00343 
00345 
00350         virtual void evalalllocal (const Entity& e, const Dune::FieldVector<DT,G::dimension>& xi, 
00351                                                            Dune::FieldVector<RT,m>& y) const
00352         {
00353           y[0] = coeff[mapper.map(e)];
00354         }
00355         
00356   private:
00357         // a reference to the grid
00358         const G& grid;
00359 
00360         // reference to index set on the grid (might be level or leaf)
00361         const IS& is;
00362 
00363         // we need a mapper
00364         MultipleCodimMultipleGeomTypeMapper<G,IS,P0Layout> mapper;
00365 
00366         // and a dynamically allocated vector
00367         const T& coeff;
00368   };
00369 
00370 
00373 }
00374 #endif