dune-grid  2.2.0
vtkwriter.hh
Go to the documentation of this file.
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set et ts=4 sw=2 sts=2:
00003 
00004 #ifndef DUNE_VTKWRITER_HH
00005 #define DUNE_VTKWRITER_HH
00006 
00007 #include <cstring>
00008 #include <iostream>
00009 #include <string>
00010 #include <fstream>
00011 #include <sstream>
00012 #include <iomanip>
00013 
00014 #include <vector>
00015 #include <list>
00016 
00017 #include <dune/common/deprecated.hh>
00018 #include <dune/common/exceptions.hh>
00019 #include <dune/common/indent.hh>
00020 #include <dune/common/iteratorfacades.hh>
00021 #include <dune/common/path.hh>
00022 #include <dune/common/shared_ptr.hh>
00023 #include <dune/geometry/referenceelements.hh>
00024 #include <dune/grid/common/mcmgmapper.hh>
00025 #include <dune/grid/common/gridenums.hh>
00026 #include <dune/grid/io/file/vtk/common.hh>
00027 #include <dune/grid/io/file/vtk/dataarraywriter.hh>
00028 #include <dune/grid/io/file/vtk/function.hh>
00029 #include <dune/grid/io/file/vtk/pvtuwriter.hh>
00030 #include <dune/grid/io/file/vtk/streams.hh>
00031 #include <dune/grid/io/file/vtk/vtuwriter.hh>
00032 
00050 namespace Dune
00051 {
00060   template< class GridView >
00061   class VTKWriter {
00062     // extract types
00063     typedef typename GridView::Grid Grid;
00064     typedef typename GridView::ctype DT;
00065     enum { n = GridView::dimension };
00066     enum { w = GridView::dimensionworld };
00067 
00068     typedef typename GridView::template Codim< 0 >::Entity Cell;
00069     typedef typename GridView::template Codim< n >::Entity Vertex;
00070     typedef Cell Entity;
00071     
00072     typedef typename GridView::IndexSet IndexSet;
00073     
00074     static const PartitionIteratorType VTK_Partition = InteriorBorder_Partition;
00075     
00076     typedef typename GridView::template Codim< 0 >
00077       ::template Partition< VTK_Partition >::Iterator
00078       GridCellIterator;
00079     typedef typename GridView::template Codim< n >
00080       ::template Partition< VTK_Partition >::Iterator
00081       GridVertexIterator;
00082     
00083     typedef MultipleCodimMultipleGeomTypeMapper< GridView, MCMGVertexLayout > VertexMapper;
00084 
00085   public:
00086     typedef Dune::VTKFunction< GridView > VTKFunction;
00087     typedef shared_ptr< const VTKFunction > VTKFunctionPtr;
00088 
00089   protected:
00090     typedef typename std::list<VTKFunctionPtr>::const_iterator FunctionIterator;
00091     
00093 
00098     class CellIterator : public GridCellIterator
00099     {
00100     public:
00102       CellIterator(const GridCellIterator & x) : GridCellIterator(x) {};
00105       const FieldVector<DT,n> position() const
00106         {
00107             return GenericReferenceElements<DT,n>::general((*this)->type()).position(0,0);
00108         }
00109     };
00110 
00111     CellIterator cellBegin() const
00112     {
00113       return gridView_.template begin< 0, VTK_Partition >();
00114     }
00115 
00116     CellIterator cellEnd() const
00117     {
00118       return gridView_.template end< 0, VTK_Partition >();
00119     }
00120     
00122 
00136     class VertexIterator :
00137       public ForwardIteratorFacade<VertexIterator, const Entity, const Entity&, int>
00138     {
00139       GridCellIterator git;
00140       GridCellIterator gend;
00141       VTK::DataMode datamode;
00142       // Index of the currently visited corner within the current element.
00143       // NOTE: this is in Dune-numbering, in contrast to CornerIterator.
00144       int cornerIndexDune;
00145       const VertexMapper & vertexmapper;
00146       std::vector<bool> visited;
00147       // in conforming mode, for each vertex id (as obtained by vertexmapper)
00148       // hold its number in the iteration order (VertexIterator)
00149       int offset;
00150     protected:
00151       void basicIncrement ()
00152       {
00153         if( git == gend )
00154           return;
00155         ++cornerIndexDune;
00156         const int numCorners = git->template count< n >();
00157         if( cornerIndexDune == numCorners )
00158         {
00159           offset += numCorners;
00160           cornerIndexDune = 0;
00161 
00162           ++git;
00163           while( (git != gend) && (git->partitionType() != InteriorEntity) )
00164             ++git;
00165         }
00166       }
00167     public:
00168       VertexIterator(const GridCellIterator & x,
00169                      const GridCellIterator & end,
00170                      const VTK::DataMode & dm,
00171                      const VertexMapper & vm) :
00172         git(x), gend(end), datamode(dm), cornerIndexDune(0),
00173         vertexmapper(vm), visited(vm.size(), false),
00174         offset(0)
00175         {
00176           if (datamode == VTK::conforming && git != gend)
00177             visited[vertexmapper.map(*git,cornerIndexDune,n)] = true;
00178         };
00179       void increment ()
00180         {
00181           switch (datamode)
00182           {
00183           case VTK::conforming:
00184             while(visited[vertexmapper.map(*git,cornerIndexDune,n)])
00185             {
00186               basicIncrement();
00187               if (git == gend) return;
00188             }
00189             visited[vertexmapper.map(*git,cornerIndexDune,n)] = true;
00190             break;
00191           case VTK::nonconforming:
00192             basicIncrement();
00193             break;
00194           }
00195        }
00196       bool equals (const VertexIterator & cit) const
00197         {
00198           return git == cit.git
00199             && cornerIndexDune == cit.cornerIndexDune
00200             && datamode == cit.datamode;
00201         }
00202       const Entity& dereference() const
00203         {
00204           return *git;
00205         }
00207       int localindex () const
00208         {
00209           return cornerIndexDune;
00210         }
00212       const FieldVector<DT,n> & position () const
00213         {
00214           return GenericReferenceElements<DT,n>::general(git->type())
00215             .position(cornerIndexDune,n);
00216         }
00217     };    
00218     
00219     VertexIterator vertexBegin () const
00220     {
00221       return VertexIterator( gridView_.template begin< 0, VTK_Partition >(),
00222                              gridView_.template end< 0, VTK_Partition >(),
00223                              datamode, *vertexmapper );
00224     }
00225 
00226     VertexIterator vertexEnd () const
00227     {
00228       return VertexIterator( gridView_.template end< 0, VTK_Partition >(),
00229                              gridView_.template end< 0, VTK_Partition >(),
00230                              datamode, *vertexmapper );
00231     }
00232     
00234 
00248     class CornerIterator :
00249       public ForwardIteratorFacade<CornerIterator, const Entity, const Entity&, int>
00250     {
00251       GridCellIterator git;
00252       GridCellIterator gend;
00253       VTK::DataMode datamode;
00254       // Index of the currently visited corner within the current element.
00255       // NOTE: this is in VTK-numbering, in contrast to VertexIterator.
00256       int cornerIndexVTK;
00257       const VertexMapper & vertexmapper;
00258       // in conforming mode, for each vertex id (as obtained by vertexmapper)
00259       // hold its number in the iteration order of VertexIterator (*not*
00260       // CornerIterator)
00261       const std::vector<int> & number;
00262       // holds the number of corners of all the elements we have seen so far,
00263       // excluding the current element
00264       int offset;
00265 
00266     public:
00267       CornerIterator(const GridCellIterator & x,
00268                      const GridCellIterator & end,
00269                      const VTK::DataMode & dm,
00270                      const VertexMapper & vm,
00271                      const std::vector<int> & num) :
00272         git(x), gend(end), datamode(dm), cornerIndexVTK(0),
00273         vertexmapper(vm),
00274         number(num), offset(0) {};
00275       void increment ()
00276       {
00277         if( git == gend )
00278           return;
00279         ++cornerIndexVTK;
00280         const int numCorners = git->template count< n >();
00281         if( cornerIndexVTK == numCorners )
00282         {
00283           offset += numCorners;
00284           cornerIndexVTK = 0;
00285 
00286           ++git;
00287           while( (git != gend) && (git->partitionType() != InteriorEntity) )
00288             ++git;
00289         }
00290       }
00291       bool equals (const CornerIterator & cit) const
00292         {
00293           return git == cit.git
00294             && cornerIndexVTK == cit.cornerIndexVTK
00295             && datamode == cit.datamode;
00296         }
00297       const Entity& dereference() const
00298         {
00299           return *git;
00300         }
00302 
00306       int id () const
00307         {
00308           switch (datamode)
00309           {
00310           case VTK::conforming:
00311             return
00312               number[vertexmapper.map(*git,VTK::renumber(*git,cornerIndexVTK),
00313                                       n)];
00314           case VTK::nonconforming:
00315             return offset + VTK::renumber(*git,cornerIndexVTK);
00316           default:
00317             DUNE_THROW(IOError,"VTKWriter: unsupported DataMode" << datamode);
00318           }
00319         }
00320     };
00321     
00322     CornerIterator cornerBegin () const
00323     {
00324       return CornerIterator( gridView_.template begin< 0, VTK_Partition >(),
00325                              gridView_.template end< 0, VTK_Partition >(),
00326                              datamode, *vertexmapper, number );
00327     }
00328     
00329     CornerIterator cornerEnd () const
00330     {
00331       return CornerIterator( gridView_.template end< 0, VTK_Partition >(),
00332                              gridView_.template end< 0, VTK_Partition >(),
00333                              datamode, *vertexmapper, number );
00334     }    
00335 
00336   public:
00344     explicit VTKWriter ( const GridView &gridView,
00345                          VTK::DataMode dm = VTK::conforming )
00346     : gridView_( gridView ),
00347       datamode( dm )
00348     { }
00349 
00354     void addCellData (const VTKFunctionPtr & p)
00355       {
00356         celldata.push_back(p);
00357       }
00358 
00364       void addCellData (VTKFunction* p) // DUNE_DEPRECATED
00365       {
00366         celldata.push_back(VTKFunctionPtr(p));
00367       }
00368 
00384     template<class V>
00385     void addCellData (const V& v, const std::string &name, int ncomps = 1)
00386     {
00387       typedef P0VTKFunction<GridView, V> Function;
00388       for (int c=0;c<ncomps;++c) {
00389         std::stringstream compName;
00390         compName << name;
00391         if (ncomps>1)
00392           compName << "[" << c << "]";
00393         VTKFunction* p = new Function(gridView_, v, compName.str(), ncomps, c);
00394         celldata.push_back(VTKFunctionPtr(p));
00395       }
00396     }
00397 
00403       void addVertexData (VTKFunction* p) // DUNE_DEPRECATED
00404       {
00405         vertexdata.push_back(VTKFunctionPtr(p));
00406       }
00407 
00412     void addVertexData (const VTKFunctionPtr & p)
00413       {
00414         vertexdata.push_back(p);
00415       }
00416 
00432     template<class V>
00433     void addVertexData (const V& v, const std::string &name, int ncomps=1)
00434     {
00435       typedef P1VTKFunction<GridView, V> Function;
00436       for (int c=0;c<ncomps;++c) {
00437         std::stringstream compName;
00438         compName << name;
00439         if (ncomps>1)
00440           compName << "[" << c << "]";
00441         VTKFunction* p = new Function(gridView_, v, compName.str(), ncomps, c);
00442         vertexdata.push_back(VTKFunctionPtr(p));
00443       }
00444     }
00445 
00447     void clear ()
00448       {
00449         celldata.clear();
00450         vertexdata.clear();
00451       }
00452 
00454     virtual ~VTKWriter ()
00455       {
00456         this->clear();
00457       }
00458 
00470     std::string write ( const std::string &name,
00471                         VTK::OutputType type = VTK::ascii )
00472     {
00473       return write( name, type, gridView_.comm().rank(), gridView_.comm().size() );
00474     }
00475 
00502     std::string pwrite ( const std::string & name,  const std::string & path, const std::string & extendpath,
00503                          VTK::OutputType type = VTK::ascii )
00504     {
00505       return pwrite( name, path, extendpath, type, gridView_.comm().rank(), gridView_.comm().size() );
00506     }
00507 
00508   protected:
00510 
00521     std::string getParallelPieceName(const std::string& name,
00522                                      const std::string& path,
00523                                      int commRank, int commSize) const
00524     {
00525       std::ostringstream s;
00526       if(path.size() > 0) {
00527         s << path;
00528         if(path[path.size()-1] != '/')
00529           s << '/';
00530       }
00531       s << 's' << std::setw(4) << std::setfill('0') << commSize << '-';
00532       s << 'p' << std::setw(4) << std::setfill('0') << commRank << '-';
00533       s << name;
00534       if(GridView::dimension > 1)
00535         s << ".vtu";
00536       else
00537         s << ".vtp";
00538       return s.str();
00539     }
00540 
00542 
00552     std::string getParallelHeaderName(const std::string& name,
00553                                       const std::string& path,
00554                                       int commSize) const
00555     {
00556       std::ostringstream s;
00557       if(path.size() > 0) {
00558         s << path;
00559         if(path[path.size()-1] != '/')
00560           s << '/';
00561       }
00562       s << 's' << std::setw(4) << std::setfill('0') << commSize << '-';
00563       s << name;
00564       if(GridView::dimension > 1)
00565         s << ".pvtu";
00566       else
00567         s << ".pvtp";
00568       return s.str();
00569     }
00570 
00572 
00584     std::string getSerialPieceName(const std::string& name,
00585                                    const std::string& path) const
00586     {
00587       static const std::string extension =
00588         GridView::dimension == 1 ? ".vtp" : ".vtu";
00589 
00590       return concatPaths(path, name+extension);
00591     }
00592 
00608     std::string write ( const std::string &name,
00609                         VTK::OutputType type,
00610                         const int commRank,
00611                         const int commSize )
00612     {
00613       // in the parallel case, just use pwrite, it has all the necessary
00614       // stuff, so we don't need to reimplement it here.
00615       if(commSize > 1)
00616         return pwrite(name, "", "", type, commRank, commSize);
00617 
00618       // make data mode visible to private functions
00619       outputtype = type;
00620 
00621       // generate filename for process data
00622       std::string pieceName = getSerialPieceName(name, "");
00623 
00624       // write process data
00625       std::ofstream file;
00626       file.open( pieceName.c_str(), std::ios::binary );
00627       if (! file.is_open())
00628         DUNE_THROW(IOError, "Could not write to piece file " << pieceName);
00629       writeDataFile( file );
00630       file.close();
00631 
00632       return pieceName;
00633     }
00634 
00636 
00659     std::string pwrite(const std::string& name, const std::string& path,
00660                        const std::string& extendpath,
00661                        VTK::OutputType ot, const int commRank,
00662                        const int commSize )
00663       {
00664         // make data mode visible to private functions
00665         outputtype=ot;
00666 
00667         // do some magic because paraview can only cope with relative pathes to piece files
00668         std::ofstream file;
00669         std::string piecepath = concatPaths(path, extendpath);
00670         std::string relpiecepath = relativePath(path, piecepath);
00671 
00672         // write this processes .vtu/.vtp piece file
00673         std::string fullname = getParallelPieceName(name, piecepath, commRank,
00674                                                     commSize);
00675         file.open(fullname.c_str(),std::ios::binary);
00676         if (! file.is_open())
00677           DUNE_THROW(IOError, "Could not write to piecefile file " << fullname);
00678         writeDataFile(file);
00679         file.close();
00680         gridView_.comm().barrier();
00681 
00682         // if we are rank 0, write .pvtu/.pvtp parallel header
00683         fullname = getParallelHeaderName(name, path, commSize);
00684         if( commRank  ==0 )
00685         {
00686           file.open(fullname.c_str());
00687           if (! file.is_open())
00688             DUNE_THROW(IOError, "Could not write to parallel file " << fullname);
00689           writeParallelHeader(file,name,relpiecepath, commSize );
00690           file.close();
00691         }
00692         gridView_.comm().barrier();
00693         return fullname;
00694       }
00695 
00696   private:
00698 
00715     void writeParallelHeader(std::ostream& s, const std::string& piecename,
00716                              const std::string& piecepath, const int commSize)
00717       {
00718         VTK::FileType fileType =
00719           (n == 1) ? VTK::polyData : VTK::unstructuredGrid;
00720 
00721         VTK::PVTUWriter writer(s, fileType);
00722 
00723         writer.beginMain();
00724 
00725         // PPointData
00726         {
00727           std::string scalars;
00728           for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end();
00729                ++it)
00730             if ((*it)->ncomps()==1)
00731             {
00732               scalars = (*it)->name();
00733               break;
00734             }
00735           std::string vectors;
00736           for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end();
00737                ++it)
00738             if ((*it)->ncomps()>1)
00739             {
00740               vectors = (*it)->name();
00741               break;
00742             }
00743           writer.beginPointData(scalars, vectors);
00744         }
00745         for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end();
00746              ++it)
00747         {
00748           unsigned writecomps = (*it)->ncomps();
00749           if(writecomps == 2) writecomps = 3;
00750           writer.addArray<float>((*it)->name(), writecomps);
00751         }
00752         writer.endPointData();
00753 
00754         // PCellData
00755         {
00756           std::string scalars;
00757           for (FunctionIterator it=celldata.begin(); it!=celldata.end();
00758                ++it)
00759             if ((*it)->ncomps()==1)
00760             {
00761               scalars = (*it)->name();
00762               break;
00763             }
00764           std::string vectors;
00765           for (FunctionIterator it=celldata.begin(); it!=celldata.end();
00766                ++it)
00767             if ((*it)->ncomps()>1)
00768             {
00769               vectors = (*it)->name();
00770               break;
00771             }
00772           writer.beginCellData(scalars, vectors);
00773         }
00774         for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it) {
00775           unsigned writecomps = (*it)->ncomps();
00776           if(writecomps == 2) writecomps = 3;
00777           writer.addArray<float>((*it)->name(), writecomps);
00778         }
00779         writer.endCellData();
00780 
00781         // PPoints
00782         writer.beginPoints();
00783         writer.addArray<float>("Coordinates", 3);
00784         writer.endPoints();
00785 
00786         // Pieces
00787         for( int i = 0; i < commSize; ++i )
00788         {
00789           const std::string& fullname = getParallelPieceName(piecename,
00790                                                              piecepath, i,
00791                                                              commSize);
00792           writer.addPiece(fullname);
00793         }
00794 
00795         writer.endMain();
00796       }
00797 
00799     void writeDataFile (std::ostream& s)
00800       {
00801         VTK::FileType fileType =
00802           (n == 1) ? VTK::polyData : VTK::unstructuredGrid;
00803 
00804         VTK::VTUWriter writer(s, outputtype, fileType);
00805 
00806         // Grid characteristics
00807         vertexmapper = new VertexMapper( gridView_ );
00808         if (datamode == VTK::conforming)
00809         {
00810           number.resize(vertexmapper->size());
00811           for (std::vector<int>::size_type i=0; i<number.size(); i++) number[i] = -1;
00812         }
00813         countEntities(nvertices, ncells, ncorners);
00814 
00815         writer.beginMain(ncells, nvertices);
00816         writeAllData(writer);
00817         writer.endMain();
00818 
00819         // write appended binary data section
00820         if(writer.beginAppended())
00821           writeAllData(writer);
00822         writer.endAppended();
00823 
00824         delete vertexmapper; number.clear();
00825       }
00826 
00827     void writeAllData(VTK::VTUWriter& writer) {
00828       // PointData
00829       writeVertexData(writer);
00830 
00831       // CellData
00832       writeCellData(writer);
00833 
00834       // Points
00835       writeGridPoints(writer);
00836 
00837       // Cells
00838       writeGridCells(writer);
00839     }
00840 
00841   protected:
00842     std::string getFormatString() const
00843       {
00844           if (outputtype==VTK::ascii)
00845               return "ascii";
00846           if (outputtype==VTK::base64)
00847               return "binary";
00848           if (outputtype==VTK::appendedraw)
00849               return "appended";
00850           if (outputtype==VTK::appendedbase64)
00851               return "appended";
00852           DUNE_THROW(IOError, "VTKWriter: unsupported OutputType" << outputtype);
00853       }
00854       
00855     std::string getTypeString() const
00856       {
00857           if (n==1)
00858               return "PolyData";
00859           else
00860               return "UnstructuredGrid";
00861       }
00862       
00864     virtual void countEntities(int &nvertices, int &ncells, int &ncorners)
00865       {
00866         nvertices = 0;
00867         ncells = 0;
00868         ncorners = 0;
00869         for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
00870         {
00871           ncells++;
00872           // because of the use of vertexmapper->map(), this iteration must be
00873           // in the order of Dune's numbering.
00874           for (int i=0; i<it->template count<n>(); ++i)
00875           {
00876             ncorners++;
00877             if (datamode == VTK::conforming)
00878             {
00879               int alpha = vertexmapper->map(*it,i,n);
00880               if (number[alpha]<0)
00881                 number[alpha] = nvertices++;
00882             }
00883             else
00884             {
00885               nvertices++;
00886             }
00887           }
00888         }
00889       }
00890 
00892     virtual void writeCellData(VTK::VTUWriter& writer)
00893       {
00894         if(celldata.size() == 0)
00895           return;
00896 
00897         std::string scalars = "";
00898         for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it)
00899           if ((*it)->ncomps()==1)
00900           {
00901             scalars = (*it)->name();
00902             break;
00903           }
00904         std::string vectors = "";
00905         for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it)
00906           if ((*it)->ncomps()>1)
00907           {
00908             vectors = (*it)->name();
00909             break;
00910           }
00911 
00912         writer.beginCellData(scalars, vectors);
00913         for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it)
00914         {
00915           // vtk file format: a vector data always should have 3 comps (with
00916           // 3rd comp = 0 in 2D case)
00917           unsigned writecomps = (*it)->ncomps();
00918           if(writecomps == 2) writecomps = 3;
00919           shared_ptr<VTK::DataArrayWriter<float> > p
00920             (writer.makeArrayWriter<float>((*it)->name(), writecomps,
00921                                             ncells));
00922           if(!p->writeIsNoop())
00923             for (CellIterator i=cellBegin(); i!=cellEnd(); ++i)
00924             {
00925               for (int j=0; j<(*it)->ncomps(); j++)
00926                 p->write((*it)->evaluate(j,*i,i.position()));
00927               // vtk file format: a vector data always should have 3 comps
00928               // (with 3rd comp = 0 in 2D case)
00929               for (unsigned j=(*it)->ncomps(); j < writecomps; ++j)
00930                 p->write(0.0);
00931             }
00932         }
00933         writer.endCellData();
00934       }
00935 
00937     virtual void writeVertexData(VTK::VTUWriter& writer)
00938       {
00939         if(vertexdata.size() == 0)
00940           return;
00941 
00942         std::string scalars = "";
00943         for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it)
00944           if ((*it)->ncomps()==1)
00945           {
00946             scalars = (*it)->name();
00947             break;
00948           }
00949         std::string vectors = "";
00950         for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it)
00951           if ((*it)->ncomps()>1)
00952           {
00953             vectors = (*it)->name();
00954             break;
00955           }
00956 
00957         writer.beginPointData(scalars, vectors);
00958         for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it)
00959         {
00960           // vtk file format: a vector data always should have 3 comps (with
00961           // 3rd comp = 0 in 2D case)
00962           unsigned writecomps = (*it)->ncomps();
00963           if(writecomps == 2) writecomps = 3;
00964           shared_ptr<VTK::DataArrayWriter<float> > p
00965             (writer.makeArrayWriter<float>((*it)->name(), writecomps,
00966                                            nvertices));
00967           if(!p->writeIsNoop())
00968             for (VertexIterator vit=vertexBegin(); vit!=vertexEnd(); ++vit)
00969             {
00970               for (int j=0; j<(*it)->ncomps(); j++)
00971                 p->write((*it)->evaluate(j,*vit,vit.position()));
00972               // vtk file format: a vector data always should have 3 comps
00973               // (with 3rd comp = 0 in 2D case)
00974               for (unsigned j=(*it)->ncomps(); j < writecomps; ++j)
00975                 p->write(0.0);
00976             }
00977         }
00978         writer.endPointData();
00979       }
00980 
00982     virtual void writeGridPoints(VTK::VTUWriter& writer)
00983       {
00984         writer.beginPoints();
00985 
00986         shared_ptr<VTK::DataArrayWriter<float> > p
00987           (writer.makeArrayWriter<float>("Coordinates", 3, nvertices));
00988         if(!p->writeIsNoop()) {
00989           VertexIterator vEnd = vertexEnd();
00990           for (VertexIterator vit=vertexBegin(); vit!=vEnd; ++vit)
00991           {
00992             int dimw=w;
00993             for (int j=0; j<std::min(dimw,3); j++)
00994               p->write(vit->geometry().corner(vit.localindex())[j]);
00995             for (int j=std::min(dimw,3); j<3; j++)
00996               p->write(0.0);
00997           }
00998         }
00999         // free the VTK::DataArrayWriter before touching the stream
01000         p.reset();
01001 
01002         writer.endPoints();
01003       }
01004 
01006     virtual void writeGridCells(VTK::VTUWriter& writer)
01007       {
01008         writer.beginCells();
01009 
01010         // connectivity
01011         {
01012           shared_ptr<VTK::DataArrayWriter<int> > p1
01013             (writer.makeArrayWriter<int>("connectivity", 1, ncorners));
01014           if(!p1->writeIsNoop())
01015             for (CornerIterator it=cornerBegin(); it!=cornerEnd(); ++it)
01016               p1->write(it.id());
01017         }
01018 
01019         // offsets
01020         {
01021           shared_ptr<VTK::DataArrayWriter<int> > p2
01022             (writer.makeArrayWriter<int>("offsets", 1, ncells));
01023           if(!p2->writeIsNoop()) {
01024             int offset = 0;
01025             for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
01026             {
01027               offset += it->template count<n>();
01028               p2->write(offset);
01029             }
01030           }
01031         }
01032 
01033         // types
01034         if (n>1)
01035         {
01036           shared_ptr<VTK::DataArrayWriter<unsigned char> > p3
01037             (writer.makeArrayWriter<unsigned char>("types", 1, ncells));
01038           if(!p3->writeIsNoop())
01039             for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
01040             {
01041               int vtktype = VTK::geometryType(it->type());
01042               p3->write(vtktype);
01043             }
01044         }
01045 
01046         writer.endCells();
01047       }
01048 
01049   protected:
01050     // the list of registered functions
01051     std::list<VTKFunctionPtr> celldata;
01052     std::list<VTKFunctionPtr> vertexdata;
01053 
01054   private: 
01055     // the grid
01056     GridView gridView_;
01057 
01058     // temporary grid information
01059   protected:
01060     int ncells;
01061     int nvertices;
01062     int ncorners;
01063   private:
01064     VertexMapper* vertexmapper;
01065     // in conforming mode, for each vertex id (as obtained by vertexmapper)
01066     // hold its number in the iteration order (VertexIterator)
01067     std::vector<int> number;
01068     VTK::DataMode datamode;
01069   protected:
01070     VTK::OutputType outputtype;
01071   };
01072 
01073 }
01074 
01075 #endif