vtkwriter.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_VTKWRITER_HH
#define DUNE_VTKWRITER_HH
 
#include <cstring>
#include <iostream>
#include <string>
#include <fstream>
#include <sstream>
#include <iomanip>
 
#include <vector>
#include <list>
 
#include <dune/common/exceptions.hh>
#include <dune/common/indent.hh>
#include <dune/common/iteratorfacades.hh>
#include <dune/common/path.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/gridenums.hh>
#include <dune/grid/io/file/vtk/common.hh>
#include <dune/grid/io/file/vtk/dataarraywriter.hh>
#include <dune/grid/io/file/vtk/function.hh>
#include <dune/grid/io/file/vtk/pvtuwriter.hh>
#include <dune/grid/io/file/vtk/streams.hh>
#include <dune/grid/io/file/vtk/vtuwriter.hh>
 
namespace Dune
{
  template< class GridView >
  class VTKWriter {
    // extract types
    typedef typename GridView::Grid Grid;
    typedef typename GridView::ctype DT;
    enum { n = GridView::dimension };
    enum { w = GridView::dimensionworld };
 
    typedef typename GridView::template Codim< 0 >::Entity Cell;
    typedef typename GridView::template Codim< n >::Entity Vertex;
    typedef Cell Entity;
 
    typedef typename GridView::IndexSet IndexSet;
 
    static const PartitionIteratorType VTK_Partition = InteriorBorder_Partition;
    //static const PartitionIteratorType VTK_Partition = All_Partition;
 
    typedef typename GridView::template Codim< 0 >
    ::template Partition< VTK_Partition >::Iterator
    GridCellIterator;
    typedef typename GridView::template Codim< n >
    ::template Partition< VTK_Partition >::Iterator
    GridVertexIterator;
 
    typedef MultipleCodimMultipleGeomTypeMapper< GridView, MCMGVertexLayout > VertexMapper;
 
    // return true if entity should be skipped in Vertex and Corner iterator
    static bool skipEntity( const PartitionType entityType )
    {
      switch( VTK_Partition )
      {
        // for All_Partition no entity has to be skipped
        case All_Partition:             return false;
        case InteriorBorder_Partition:  return ( entityType != InteriorEntity );
        default: DUNE_THROW(NotImplemented,"Add check for this partition type");
      }
      return false ;
    }
 
  public:
    typedef Dune::VTKFunction< GridView > VTKFunction;
    typedef shared_ptr< const VTKFunction > VTKFunctionPtr;
 
  protected:
    typedef typename std::list<VTKFunctionPtr>::const_iterator FunctionIterator;
 
 
    class CellIterator : public GridCellIterator
    {
    public:
      CellIterator(const GridCellIterator & x) : GridCellIterator(x) {}
      const FieldVector<DT,n> position() const
      {
        return ReferenceElements<DT,n>::general((*this)->type()).position(0,0);
      }
    };
 
    CellIterator cellBegin() const
    {
      return gridView_.template begin< 0, VTK_Partition >();
    }
 
    CellIterator cellEnd() const
    {
      return gridView_.template end< 0, VTK_Partition >();
    }
 
 
    class VertexIterator :
      public ForwardIteratorFacade<VertexIterator, const Entity, const Entity&, int>
    {
      GridCellIterator git;
      GridCellIterator gend;
      VTK::DataMode datamode;
      // Index of the currently visited corner within the current element.
      // NOTE: this is in Dune-numbering, in contrast to CornerIterator.
      int cornerIndexDune;
      const VertexMapper & vertexmapper;
      std::vector<bool> visited;
      // in conforming mode, for each vertex id (as obtained by vertexmapper)
      // hold its number in the iteration order (VertexIterator)
      int offset;
    protected:
      void basicIncrement ()
      {
        if( git == gend )
          return;
        ++cornerIndexDune;
        const int numCorners = git->template count< n >();
        if( cornerIndexDune == numCorners )
        {
          offset += numCorners;
          cornerIndexDune = 0;
 
          ++git;
          while( (git != gend) && skipEntity( git->partitionType() ) )
            ++git;
        }
      }
    public:
      VertexIterator(const GridCellIterator & x,
                     const GridCellIterator & end,
                     const VTK::DataMode & dm,
                     const VertexMapper & vm) :
        git(x), gend(end), datamode(dm), cornerIndexDune(0),
        vertexmapper(vm), visited(vm.size(), false),
        offset(0)
      {
        if (datamode == VTK::conforming && git != gend)
          visited[vertexmapper.map(*git,cornerIndexDune,n)] = true;
      }
      void increment ()
      {
        switch (datamode)
        {
        case VTK::conforming :
          while(visited[vertexmapper.map(*git,cornerIndexDune,n)])
          {
            basicIncrement();
            if (git == gend) return;
          }
          visited[vertexmapper.map(*git,cornerIndexDune,n)] = true;
          break;
        case VTK::nonconforming :
          basicIncrement();
          break;
        }
      }
      bool equals (const VertexIterator & cit) const
      {
        return git == cit.git
               && cornerIndexDune == cit.cornerIndexDune
               && datamode == cit.datamode;
      }
      const Entity& dereference() const
      {
        return *git;
      }
      int localindex () const
      {
        return cornerIndexDune;
      }
      const FieldVector<DT,n> & position () const
      {
        return ReferenceElements<DT,n>::general(git->type())
               .position(cornerIndexDune,n);
      }
    };
 
    VertexIterator vertexBegin () const
    {
      return VertexIterator( gridView_.template begin< 0, VTK_Partition >(),
                             gridView_.template end< 0, VTK_Partition >(),
                             datamode, *vertexmapper );
    }
 
    VertexIterator vertexEnd () const
    {
      return VertexIterator( gridView_.template end< 0, VTK_Partition >(),
                             gridView_.template end< 0, VTK_Partition >(),
                             datamode, *vertexmapper );
    }
 
 
    class CornerIterator :
      public ForwardIteratorFacade<CornerIterator, const Entity, const Entity&, int>
    {
      GridCellIterator git;
      GridCellIterator gend;
      VTK::DataMode datamode;
      // Index of the currently visited corner within the current element.
      // NOTE: this is in VTK-numbering, in contrast to VertexIterator.
      int cornerIndexVTK;
      const VertexMapper & vertexmapper;
      // in conforming mode, for each vertex id (as obtained by vertexmapper)
      // hold its number in the iteration order of VertexIterator (*not*
      // CornerIterator)
      const std::vector<int> & number;
      // holds the number of corners of all the elements we have seen so far,
      // excluding the current element
      int offset;
 
    public:
      CornerIterator(const GridCellIterator & x,
                     const GridCellIterator & end,
                     const VTK::DataMode & dm,
                     const VertexMapper & vm,
                     const std::vector<int> & num) :
        git(x), gend(end), datamode(dm), cornerIndexVTK(0),
        vertexmapper(vm),
        number(num), offset(0) {}
      void increment ()
      {
        if( git == gend )
          return;
        ++cornerIndexVTK;
        const int numCorners = git->template count< n >();
        if( cornerIndexVTK == numCorners )
        {
          offset += numCorners;
          cornerIndexVTK = 0;
 
          ++git;
          while( (git != gend) && skipEntity( git->partitionType() ) )
            ++git;
        }
      }
      bool equals (const CornerIterator & cit) const
      {
        return git == cit.git
               && cornerIndexVTK == cit.cornerIndexVTK
               && datamode == cit.datamode;
      }
      const Entity& dereference() const
      {
        return *git;
      }
 
      int id () const
      {
        switch (datamode)
        {
        case VTK::conforming :
          return
            number[vertexmapper.map(*git,VTK::renumber(*git,cornerIndexVTK),
                                    n)];
        case VTK::nonconforming :
          return offset + VTK::renumber(*git,cornerIndexVTK);
        default :
          DUNE_THROW(IOError,"VTKWriter: unsupported DataMode" << datamode);
        }
      }
    };
 
    CornerIterator cornerBegin () const
    {
      return CornerIterator( gridView_.template begin< 0, VTK_Partition >(),
                             gridView_.template end< 0, VTK_Partition >(),
                             datamode, *vertexmapper, number );
    }
 
    CornerIterator cornerEnd () const
    {
      return CornerIterator( gridView_.template end< 0, VTK_Partition >(),
                             gridView_.template end< 0, VTK_Partition >(),
                             datamode, *vertexmapper, number );
    }
 
  public:
    explicit VTKWriter ( const GridView &gridView,
                         VTK::DataMode dm = VTK::conforming )
      : gridView_( gridView ),
        datamode( dm )
    { }
 
    void addCellData (const VTKFunctionPtr & p)
    {
      celldata.push_back(p);
    }
 
    void addCellData (VTKFunction* p)
    {
      celldata.push_back(VTKFunctionPtr(p));
    }
 
    template<class V>
    void addCellData (const V& v, const std::string &name, int ncomps = 1)
    {
      typedef P0VTKFunction<GridView, V> Function;
      for (int c=0; c<ncomps; ++c) {
        std::stringstream compName;
        compName << name;
        if (ncomps>1)
          compName << "[" << c << "]";
        VTKFunction* p = new Function(gridView_, v, compName.str(), ncomps, c);
        celldata.push_back(VTKFunctionPtr(p));
      }
    }
 
    void addVertexData (VTKFunction* p)
    {
      vertexdata.push_back(VTKFunctionPtr(p));
    }
 
    void addVertexData (const VTKFunctionPtr & p)
    {
      vertexdata.push_back(p);
    }
 
    template<class V>
    void addVertexData (const V& v, const std::string &name, int ncomps=1)
    {
      typedef P1VTKFunction<GridView, V> Function;
      for (int c=0; c<ncomps; ++c) {
        std::stringstream compName;
        compName << name;
        if (ncomps>1)
          compName << "[" << c << "]";
        VTKFunction* p = new Function(gridView_, v, compName.str(), ncomps, c);
        vertexdata.push_back(VTKFunctionPtr(p));
      }
    }
 
    void clear ()
    {
      celldata.clear();
      vertexdata.clear();
    }
 
    virtual ~VTKWriter ()
    {
      this->clear();
    }
 
    std::string write ( const std::string &name,
                        VTK::OutputType type = VTK::ascii )
    {
      return write( name, type, gridView_.comm().rank(), gridView_.comm().size() );
    }
 
    std::string pwrite ( const std::string & name,  const std::string & path, const std::string & extendpath,
                         VTK::OutputType type = VTK::ascii )
    {
      return pwrite( name, path, extendpath, type, gridView_.comm().rank(), gridView_.comm().size() );
    }
 
  protected:
 
    std::string getParallelPieceName(const std::string& name,
                                     const std::string& path,
                                     int commRank, int commSize) const
    {
      std::ostringstream s;
      if(path.size() > 0) {
        s << path;
        if(path[path.size()-1] != '/')
          s << '/';
      }
      s << 's' << std::setw(4) << std::setfill('0') << commSize << '-';
      s << 'p' << std::setw(4) << std::setfill('0') << commRank << '-';
      s << name;
      if(GridView::dimension > 1)
        s << ".vtu";
      else
        s << ".vtp";
      return s.str();
    }
 
 
    std::string getParallelHeaderName(const std::string& name,
                                      const std::string& path,
                                      int commSize) const
    {
      std::ostringstream s;
      if(path.size() > 0) {
        s << path;
        if(path[path.size()-1] != '/')
          s << '/';
      }
      s << 's' << std::setw(4) << std::setfill('0') << commSize << '-';
      s << name;
      if(GridView::dimension > 1)
        s << ".pvtu";
      else
        s << ".pvtp";
      return s.str();
    }
 
 
    std::string getSerialPieceName(const std::string& name,
                                   const std::string& path) const
    {
      static const std::string extension =
        GridView::dimension == 1 ? ".vtp" : ".vtu";
 
      return concatPaths(path, name+extension);
    }
 
    std::string write ( const std::string &name,
                        VTK::OutputType type,
                        const int commRank,
                        const int commSize )
    {
      // in the parallel case, just use pwrite, it has all the necessary
      // stuff, so we don't need to reimplement it here.
      if(commSize > 1)
        return pwrite(name, "", "", type, commRank, commSize);
 
      // make data mode visible to private functions
      outputtype = type;
 
      // generate filename for process data
      std::string pieceName = getSerialPieceName(name, "");
 
      // write process data
      std::ofstream file;
      file.exceptions(std::ios_base::badbit | std::ios_base::failbit |
                      std::ios_base::eofbit);
      file.open( pieceName.c_str(), std::ios::binary );
      if (! file.is_open())
        DUNE_THROW(IOError, "Could not write to piece file " << pieceName);
      writeDataFile( file );
      file.close();
 
      return pieceName;
    }
 
 
    std::string pwrite(const std::string& name, const std::string& path,
                       const std::string& extendpath,
                       VTK::OutputType ot, const int commRank,
                       const int commSize )
    {
      // make data mode visible to private functions
      outputtype=ot;
 
      // do some magic because paraview can only cope with relative pathes to piece files
      std::ofstream file;
      file.exceptions(std::ios_base::badbit | std::ios_base::failbit |
                      std::ios_base::eofbit);
      std::string piecepath = concatPaths(path, extendpath);
      std::string relpiecepath = relativePath(path, piecepath);
 
      // write this processes .vtu/.vtp piece file
      std::string fullname = getParallelPieceName(name, piecepath, commRank,
                                                  commSize);
      file.open(fullname.c_str(),std::ios::binary);
      if (! file.is_open())
        DUNE_THROW(IOError, "Could not write to piecefile file " << fullname);
      writeDataFile(file);
      file.close();
      gridView_.comm().barrier();
 
      // if we are rank 0, write .pvtu/.pvtp parallel header
      fullname = getParallelHeaderName(name, path, commSize);
      if( commRank  ==0 )
      {
        file.open(fullname.c_str());
        if (! file.is_open())
          DUNE_THROW(IOError, "Could not write to parallel file " << fullname);
        writeParallelHeader(file,name,relpiecepath, commSize );
        file.close();
      }
      gridView_.comm().barrier();
      return fullname;
    }
 
  private:
 
    void writeParallelHeader(std::ostream& s, const std::string& piecename,
                             const std::string& piecepath, const int commSize)
    {
      VTK::FileType fileType =
        (n == 1) ? VTK::polyData : VTK::unstructuredGrid;
 
      VTK::PVTUWriter writer(s, fileType);
 
      writer.beginMain();
 
      // PPointData
      {
        std::string scalars;
        for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end();
             ++it)
          if ((*it)->ncomps()==1)
          {
            scalars = (*it)->name();
            break;
          }
        std::string vectors;
        for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end();
             ++it)
          if ((*it)->ncomps()>1)
          {
            vectors = (*it)->name();
            break;
          }
        writer.beginPointData(scalars, vectors);
      }
      for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end();
           ++it)
      {
        unsigned writecomps = (*it)->ncomps();
        if(writecomps == 2) writecomps = 3;
        writer.addArray<float>((*it)->name(), writecomps);
      }
      writer.endPointData();
 
      // PCellData
      {
        std::string scalars;
        for (FunctionIterator it=celldata.begin(); it!=celldata.end();
             ++it)
          if ((*it)->ncomps()==1)
          {
            scalars = (*it)->name();
            break;
          }
        std::string vectors;
        for (FunctionIterator it=celldata.begin(); it!=celldata.end();
             ++it)
          if ((*it)->ncomps()>1)
          {
            vectors = (*it)->name();
            break;
          }
        writer.beginCellData(scalars, vectors);
      }
      for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it) {
        unsigned writecomps = (*it)->ncomps();
        if(writecomps == 2) writecomps = 3;
        writer.addArray<float>((*it)->name(), writecomps);
      }
      writer.endCellData();
 
      // PPoints
      writer.beginPoints();
      writer.addArray<float>("Coordinates", 3);
      writer.endPoints();
 
      // Pieces
      for( int i = 0; i < commSize; ++i )
      {
        const std::string& fullname = getParallelPieceName(piecename,
                                                           piecepath, i,
                                                           commSize);
        writer.addPiece(fullname);
      }
 
      writer.endMain();
    }
 
    void writeDataFile (std::ostream& s)
    {
      VTK::FileType fileType =
        (n == 1) ? VTK::polyData : VTK::unstructuredGrid;
 
      VTK::VTUWriter writer(s, outputtype, fileType);
 
      // Grid characteristics
      vertexmapper = new VertexMapper( gridView_ );
      if (datamode == VTK::conforming)
      {
        number.resize(vertexmapper->size());
        for (std::vector<int>::size_type i=0; i<number.size(); i++) number[i] = -1;
      }
      countEntities(nvertices, ncells, ncorners);
 
      writer.beginMain(ncells, nvertices);
      writeAllData(writer);
      writer.endMain();
 
      // write appended binary data section
      if(writer.beginAppended())
        writeAllData(writer);
      writer.endAppended();
 
      delete vertexmapper; number.clear();
    }
 
    void writeAllData(VTK::VTUWriter& writer) {
      // PointData
      writeVertexData(writer);
 
      // CellData
      writeCellData(writer);
 
      // Points
      writeGridPoints(writer);
 
      // Cells
      writeGridCells(writer);
    }
 
  protected:
    std::string getFormatString() const
    {
      if (outputtype==VTK::ascii)
        return "ascii";
      if (outputtype==VTK::base64)
        return "binary";
      if (outputtype==VTK::appendedraw)
        return "appended";
      if (outputtype==VTK::appendedbase64)
        return "appended";
      DUNE_THROW(IOError, "VTKWriter: unsupported OutputType" << outputtype);
    }
 
    std::string getTypeString() const
    {
      if (n==1)
        return "PolyData";
      else
        return "UnstructuredGrid";
    }
 
    virtual void countEntities(int &nvertices, int &ncells, int &ncorners)
    {
      nvertices = 0;
      ncells = 0;
      ncorners = 0;
      for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
      {
        ncells++;
        // because of the use of vertexmapper->map(), this iteration must be
        // in the order of Dune's numbering.
        for (int i=0; i<it->template count<n>(); ++i)
        {
          ncorners++;
          if (datamode == VTK::conforming)
          {
            int alpha = vertexmapper->map(*it,i,n);
            if (number[alpha]<0)
              number[alpha] = nvertices++;
          }
          else
          {
            nvertices++;
          }
        }
      }
    }
 
    virtual void writeCellData(VTK::VTUWriter& writer)
    {
      if(celldata.size() == 0)
        return;
 
      std::string scalars = "";
      for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it)
        if ((*it)->ncomps()==1)
        {
          scalars = (*it)->name();
          break;
        }
      std::string vectors = "";
      for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it)
        if ((*it)->ncomps()>1)
        {
          vectors = (*it)->name();
          break;
        }
 
      writer.beginCellData(scalars, vectors);
      for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it)
      {
        // vtk file format: a vector data always should have 3 comps (with
        // 3rd comp = 0 in 2D case)
        unsigned writecomps = (*it)->ncomps();
        if(writecomps == 2) writecomps = 3;
        shared_ptr<VTK::DataArrayWriter<float> > p
          (writer.makeArrayWriter<float>((*it)->name(), writecomps,
                                         ncells));
        if(!p->writeIsNoop())
          for (CellIterator i=cellBegin(); i!=cellEnd(); ++i)
          {
            for (int j=0; j<(*it)->ncomps(); j++)
              p->write((*it)->evaluate(j,*i,i.position()));
            // vtk file format: a vector data always should have 3 comps
            // (with 3rd comp = 0 in 2D case)
            for (unsigned j=(*it)->ncomps(); j < writecomps; ++j)
              p->write(0.0);
          }
      }
      writer.endCellData();
    }
 
    virtual void writeVertexData(VTK::VTUWriter& writer)
    {
      if(vertexdata.size() == 0)
        return;
 
      std::string scalars = "";
      for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it)
        if ((*it)->ncomps()==1)
        {
          scalars = (*it)->name();
          break;
        }
      std::string vectors = "";
      for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it)
        if ((*it)->ncomps()>1)
        {
          vectors = (*it)->name();
          break;
        }
 
      writer.beginPointData(scalars, vectors);
      for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it)
      {
        // vtk file format: a vector data always should have 3 comps (with
        // 3rd comp = 0 in 2D case)
        unsigned writecomps = (*it)->ncomps();
        if(writecomps == 2) writecomps = 3;
        shared_ptr<VTK::DataArrayWriter<float> > p
          (writer.makeArrayWriter<float>((*it)->name(), writecomps,
                                         nvertices));
        if(!p->writeIsNoop())
          for (VertexIterator vit=vertexBegin(); vit!=vertexEnd(); ++vit)
          {
            for (int j=0; j<(*it)->ncomps(); j++)
              p->write((*it)->evaluate(j,*vit,vit.position()));
            // vtk file format: a vector data always should have 3 comps
            // (with 3rd comp = 0 in 2D case)
            for (unsigned j=(*it)->ncomps(); j < writecomps; ++j)
              p->write(0.0);
          }
      }
      writer.endPointData();
    }
 
    virtual void writeGridPoints(VTK::VTUWriter& writer)
    {
      writer.beginPoints();
 
      shared_ptr<VTK::DataArrayWriter<float> > p
        (writer.makeArrayWriter<float>("Coordinates", 3, nvertices));
      if(!p->writeIsNoop()) {
        VertexIterator vEnd = vertexEnd();
        for (VertexIterator vit=vertexBegin(); vit!=vEnd; ++vit)
        {
          int dimw=w;
          for (int j=0; j<std::min(dimw,3); j++)
            p->write(vit->geometry().corner(vit.localindex())[j]);
          for (int j=std::min(dimw,3); j<3; j++)
            p->write(0.0);
        }
      }
      // free the VTK::DataArrayWriter before touching the stream
      p.reset();
 
      writer.endPoints();
    }
 
    virtual void writeGridCells(VTK::VTUWriter& writer)
    {
      writer.beginCells();
 
      // connectivity
      {
        shared_ptr<VTK::DataArrayWriter<int> > p1
          (writer.makeArrayWriter<int>("connectivity", 1, ncorners));
        if(!p1->writeIsNoop())
          for (CornerIterator it=cornerBegin(); it!=cornerEnd(); ++it)
            p1->write(it.id());
      }
 
      // offsets
      {
        shared_ptr<VTK::DataArrayWriter<int> > p2
          (writer.makeArrayWriter<int>("offsets", 1, ncells));
        if(!p2->writeIsNoop()) {
          int offset = 0;
          for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
          {
            offset += it->template count<n>();
            p2->write(offset);
          }
        }
      }
 
      // types
      if (n>1)
      {
        shared_ptr<VTK::DataArrayWriter<unsigned char> > p3
          (writer.makeArrayWriter<unsigned char>("types", 1, ncells));
        if(!p3->writeIsNoop())
          for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
          {
            int vtktype = VTK::geometryType(it->type());
            p3->write(vtktype);
          }
      }
 
      writer.endCells();
    }
 
  protected:
    // the list of registered functions
    std::list<VTKFunctionPtr> celldata;
    std::list<VTKFunctionPtr> vertexdata;
 
    // the grid
    GridView gridView_;
 
    // temporary grid information
    int ncells;
    int nvertices;
    int ncorners;
  private:
    VertexMapper* vertexmapper;
    // in conforming mode, for each vertex id (as obtained by vertexmapper)
    // hold its number in the iteration order (VertexIterator)
    std::vector<int> number;
    VTK::DataMode datamode;
  protected:
    VTK::OutputType outputtype;
  };
 
}
 
#endif