added RectilinearGridWriter

dune/vtk/datacollectors/structureddatacollector.hh

@@ -30,6 +30,7 @@ public:
   StructuredDataCollectorInterface (GridView const& gridView)
     : Super(gridView)
     , defaultDataCollector_(gridView)
+    , ghostLevel_(gridView.overlapSize(0))
   {}
 
   /// Return number of grid vertices
@@ -70,6 +71,21 @@ public:
     this->asDerived().writePiecesImpl(writer);
   }
 
+  int ghostLevel () const
+  {
+    return this->asDerived().ghostLevelImpl();
+  }
+
+  /// Return the coordinates along the ordinates x, y, and z
+  template <class T>
+  std::array<std::vector<T>, 3> coordinates () const
+  {
+    return this->asDerived().template coordinatesImpl<T>();
+  }
+
+
+public:
+
   /// Return the coordinates of all grid vertices in the order given by the indexSet
   template <class T>
   std::vector<T> pointsImpl () const
@@ -84,8 +100,39 @@ public:
     return defaultDataCollector_.template pointData<T>(fct);
   }
 
+
+private:
+
+  int ghostLevelImpl () const
+  {
+    return ghostLevel_;
+  }
+
+  template <class T>
+  std::array<std::vector<T>, 3> coordinatesImpl () const
+  {
+    auto origin = this->origin();
+    auto spacing = this->spacing();
+
+    std::array<std::vector<T>, 3> ordinates{};
+    writeLocalPiece([&ordinates,&origin,&spacing](auto const& extent) {
+      for (std::size_t d = 0; d < GridView::dimension; ++d) {
+        auto s = extent[2*d+1] - extent[2*d] + 1;
+        ordinates[d].resize(s);
+        for (std::size_t i = 0; i < s; ++i)
+          ordinates[d][i] = origin[d] + (extent[2*d] + i)*spacing[d];
+      }
+    });
+
+    for (std::size_t d = GridView::dimension; d < 3; ++d)
+      ordinates[d].resize(1, T(0));
+
+    return ordinates;
+  }
+
 private:
   DefaultDataCollector<GridView> defaultDataCollector_;
+  int ghostLevel_;
 };
 
 }} // end namespace Dune::experimental

dune/vtk/datacollectors/yaspdatacollector.hh

@@ -113,6 +113,29 @@ public:
     }
   }
 
+  template <class T>
+  std::array<std::vector<T>, 3> coordinatesImpl () const
+  {
+    auto it = gridView_.grid().begin(level_);
+    auto const& coords = it->coords;
+
+    std::array<std::vector<T>, 3> ordinates{};
+    writeLocalPieceImpl([&ordinates,&coords](auto const& extent)
+    {
+      for (std::size_t d = 0; d < dim; ++d) {
+        auto s = extent[2*d+1] - extent[2*d] + 1;
+        ordinates[d].resize(s);
+        for (std::size_t i = 0; i < s; ++i)
+          ordinates[d][i] = coords.coordinate(d, extent[2*d] + i);
+      }
+    });
+
+    for (std::size_t d = dim; d < 3; ++d)
+      ordinates[d].resize(1, T(0));
+
+    return ordinates;
+  }
+
 private:
   std::array<int, 6> wholeExtent_;
   FieldVector<ctype,3> spacing_;

dune/vtk/vtktypes.cc

@@ -21,14 +21,14 @@ std::string to_string (FormatTypes type)
 std::string to_string (DataTypes type)
 {
   switch (type) {
-    case INT8: return "Int8";
-    case UINT8: return "UInt8";
-    case INT16: return "Int16";
-    case UINT16: return "UInt16";
-    case INT32: return "Int32";
-    case UINT32: return "UInt32";
-    case INT64: return "Int64";
-    case UINT64: return "UInt64";
+    case INT8:    return "Int8";
+    case UINT8:   return "UInt8";
+    case INT16:   return "Int16";
+    case UINT16:  return "UInt16";
+    case INT32:   return "Int32";
+    case UINT32:  return "UInt32";
+    case INT64:   return "Int64";
+    case UINT64:  return "UInt64";
     case FLOAT32: return "Float32";
     case FLOAT64: return "Float64";
     default:

dune/vtk/vtkwriter.hh

 #pragma once
 
-#include <array>
-#include <iosfwd>
-#include <map>
+#include <dune/vtk/writers/vtkimagedatawriter.hh>
+#include <dune/vtk/writers/vtkrectilineargridwriter.hh>
+#include <dune/vtk/writers/vtkstructuredgridwriter.hh>
+#include <dune/vtk/writers/vtkunstructuredgridwriter.hh>
 
-#include <dune/common/std/optional.hh>
+#if HAVE_DUNE_SPGRID
+#include <dune/grid/spgrid.hh>
+#include <dune/vtk/datacollectors/spdatacollector.hh>
+#endif
 
-#include "datacollector.hh"
-#include "filewriter.hh"
-#include "vtkfunction.hh"
-#include "vtktypes.hh"
+#include <dune/grid/geometrygrid.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/vtk/datacollectors/yaspdatacollector.hh>
 
 namespace Dune { namespace experimental
 {
-  /// File-Writer for Vtk .vtu files
-  template <class GridView, class DataCollector>
-  class VtkWriter
-      : public FileWriter
+  namespace Impl
   {
-  protected:
-    static constexpr int dimension = GridView::dimension;
-
-    using GlobalFunction = VTKFunction<GridView>;
-    using LocalFunction = VTKLocalFunction<GridView>;
-    using pos_type = typename std::ostream::pos_type;
-
-    enum PositionTypes {
-      POINT_DATA,
-      CELL_DATA
-    };
-
-  public:
-    /// Constructor, stores the gridView
-    VtkWriter (GridView const& gridView)
-      : dataCollector_(gridView)
-    {}
-
-    /// Write the attached data to the file
-    virtual void write (std::string const& fn) override
-    {
-      write(fn, Vtk::BINARY);
-    }
-
-    /// Write the attached data to the file with \ref Vtk::FormatTypes and \ref Vtk::DataTypes
-    void write (std::string const& fn,
-                Vtk::FormatTypes format,
-                Vtk::DataTypes datatype = Vtk::FLOAT32);
-
-    /// Attach point data to the writer
-    template <class GridViewFunction>
-    VtkWriter& addPointData (GridViewFunction const& gridViewFct,
-                              std::string const& name = {},
-                              int ncomps = 1)
-    {
-      pointData_.emplace_back(gridViewFct, name, ncomps);
-      return *this;
-    }
-
-    /// Attach cell data to the writer
-    template <class GridViewFunction>
-    VtkWriter& addCellData (GridViewFunction const& gridViewFct,
-                            std::string const& name = {},
-                            int ncomps = 1)
-    {
-      cellData_.emplace_back(gridViewFct, name, ncomps);
-      return *this;
-    }
-
-  protected:
-    /// Write a serial VTK file in Unstructured format
-    virtual void writeSerialFile (std::string const& filename) const = 0;
-
-    /// Write a parallel VTK file `pfilename.pvtu` in Unstructured format,
-    /// with `size` the number of pieces and serial files given by `pfilename_p[i].vtu`
-    /// for [i] in [0,...,size).
-    virtual void writeParallelFile (std::string const& pfilename, int size) const = 0;
-
-    /// Return the file extension of the serial file (not including the dot)
-    virtual std::string fileExtension () const = 0;
-
-    // Write the point or cell values given by the grid function `fct` to the
-    // output stream `out`. In case of binary format, stores the streampos of XML
-    // attributes "offset" in the vector `offsets`.
-    void writeData (std::ofstream& out,
-                    std::vector<pos_type>& offsets,
-                    GlobalFunction const& fct,
-                    PositionTypes type) const;
-
-    // Write the coordinates of the vertices to the output stream `out`. In case
-    // of binary format, stores the streampos of XML attributes "offset" in the
-    // vector `offsets`.
-    void writePoints (std::ofstream& out,
-                      std::vector<pos_type>& offsets) const;
-
-    // Write the element connectivity to the output stream `out`. In case
-    // of binary format, stores the streampos of XML attributes "offset" in the
-    // vector `offsets`.
-    void writeCells (std::ofstream& oust,
-                     std::vector<pos_type>& offsets) const;
-
-    // Collect point or cell data (depending on \ref PositionTypes) and pass
-    // the resulting vector to \ref writeAppended.
-    template <class T>
-    std::uint64_t writeDataAppended (std::ofstream& out,
-                                     GlobalFunction const& localFct,
-                                     PositionTypes type) const;
-
-    // Collect point positions and pass the resulting vector to \ref writeAppended.
-    template <class T>
-    std::uint64_t writePointsAppended (std::ofstream& out) const;
-
-    // Collect element connectivity, offsets and element types, and pass the
-    // resulting vectors to \ref writeAppended.
-    std::array<std::uint64_t,3> writeCellsAppended (std::ofstream& out) const;
-
-    // Write the `values` in blocks (possibly compressed) to the output
-    // stream `out`. Return the written block size.
-    template <class T>
-    std::uint64_t writeAppended (std::ofstream& out, std::vector<T> const& values) const;
-
-
-    Std::optional<std::string> getScalarName (std::vector<GlobalFunction> const& data) const
+    // The default writer assumes an unstructured grid
+    template <class GridView, class Grid>
+    struct VtkWriterImpl
     {
-      auto scalar = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.ncomps() == 1; });
-      return scalar != data.end() ? Std::optional<std::string>{scalar->name()} : Std::optional<std::string>{};
-    }
+      using type = VtkUnstructuredGridWriter<GridView>;
+    };
 
-    Std::optional<std::string> getVectorName (std::vector<GlobalFunction> const& data) const
+#if HAVE_DUNE_SPGRID
+    // A structured grid with constant spacing in x, y, and z direction.
+    template <class GridView, class ct, int dim, template <int> class Ref, class Comm>
+    struct VtkWriterImpl<GridView, SPGrid<ct,dim,Ref,Comm>>
     {
-      auto vector = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.ncomps() == 3; });
-      return vector != data.end() ? Std::optional<std::string>{vector->name()} : Std::optional<std::string>{};
-    }
+      using type = VtkImageDataWriter<GridView, SPDataCollector<GridView>>;
+    };
+#endif
 
-    Std::optional<std::string> getTensorName (std::vector<GlobalFunction> const& data) const
+    // A structured grid with constant spacing in x, y, and z direction.
+    template <class GridView, int dim, class Coordinates>
+    struct VtkWriterImpl<GridView, YaspGrid<dim,Coordinates>>
     {
-      auto tensor = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.ncomps() == 9; });
-      return tensor != data.end() ? Std::optional<std::string>{tensor->name()} : Std::optional<std::string>{};
-    }
+      using type = VtkImageDataWriter<GridView, YaspDataCollector<GridView>>;
+    };
 
-    std::string getNames (std::vector<GlobalFunction> const& data) const
+    // A structured grid with coordinates in x, y, and z direction with arbitrary spacing
+    template <class GridView, int dim, class ct>
+    struct VtkWriterImpl<GridView, YaspGrid<dim,TensorProductCoordinates<ct,dim>>>
     {
-      auto n1 = getScalarName(data);
-      auto n2 = getVectorName(data);
-      auto n3 = getScalarName(data);
-      return (n1 ? " Scalars=\"" + *n1 + "\"" : "")
-           + (n2 ? " Vectors=\"" + *n2 + "\"" : "")
-           + (n3 ? " Tensors=\"" + *n3 + "\"" : "");
-    }
+      using type = VtkRectilinearGridWriter<GridView, YaspDataCollector<GridView>>;
+    };
 
-    // Returns endianness
-    std::string getEndian () const
+    // A transformed structured grid has structured connectivity but unstructured point
+    // coordinates.
+    template <class GridView, int dim, class Coordinates, class CoordFunction, class Allocator>
+    struct VtkWriterImpl<GridView, GeometryGrid<YaspGrid<dim,Coordinates>, CoordFunction, Allocator>>
     {
-      short i = 1;
-      return (reinterpret_cast<char*>(&i)[1] == 1 ? "BigEndian" : "LittleEndian");
-    }
-
-  protected:
-    mutable DataCollector dataCollector_;
+      using type = VtkStructuredGridWriter<GridView, YaspDataCollector<GridView>>;
+    };
 
-    std::string filename_;
-    Vtk::FormatTypes format_;
-    Vtk::DataTypes datatype_;
+  } // end namespace Impl
 
-    // attached data
-    std::vector<GlobalFunction> pointData_;
-    std::vector<GlobalFunction> cellData_;
 
-    std::size_t const block_size = 1024*32;
-    int compression_level = -1; // in [0,9], -1 ... use default value
-  };
+  /// Default choice for several grid types, uses the default data-collector.
+  template <class GridView>
+  using VtkWriter = Impl::VtkWriterImpl<GridView, typename GridView::Grid>;
 
 }} // end namespace Dune::experimental
-
-#include "vtkwriter.impl.hh"

dune/vtk/vtkimagedatawriter.hh → dune/vtk/writers/vtkimagedatawriter.hh

@@ -4,23 +4,24 @@
 #include <iosfwd>
 #include <map>
 
-#include "datacollector.hh"
-#include "filewriter.hh"
-#include "vtkfunction.hh"
-#include "vtktypes.hh"
-#include "vtkwriter.hh"
-#include "datacollectors/structureddatacollector.hh"
+#include <dune/vtk/datacollector.hh>
+#include <dune/vtk/filewriter.hh>
+#include <dune/vtk/vtkfunction.hh>
+#include <dune/vtk/vtktypes.hh>
+#include <dune/vtk/datacollectors/structureddatacollector.hh>
+
+#include "vtkwriterinterface.hh"
 
 namespace Dune { namespace experimental
 {
   /// File-Writer for VTK .vtu files
   template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
   class VtkImageDataWriter
-      : public VtkWriter<GridView, DataCollector>
+      : public VtkWriterInterface<GridView, DataCollector>
   {
     static constexpr int dimension = GridView::dimension;
 
-    using Super = VtkWriter<GridView, DataCollector>;
+    using Super = VtkWriterInterface<GridView, DataCollector>;
     using pos_type = typename Super::pos_type;
 
   public:

dune/vtk/vtkimagedatawriter.impl.hh → dune/vtk/writers/vtkimagedatawriter.impl.hh

@@ -10,9 +10,9 @@
 #include <dune/geometry/referenceelements.hh>
 #include <dune/geometry/type.hh>
 
-#include "utility/enum.hh"
-#include "utility/filesystem.hh"
-#include "utility/string.hh"
+#include <dune/vtk/utility/enum.hh>
+#include <dune/vtk/utility/filesystem.hh>
+#include <dune/vtk/utility/string.hh>
 
 namespace Dune { namespace experimental {
 
@@ -21,6 +21,8 @@ void VtkImageDataWriter<GV,DC>
   ::writeSerialFile (std::string const& filename) const
 {
   std::ofstream out(filename, std::ios_base::ate | std::ios::binary);
+  assert(out.is_open());
+
   if (format_ == Vtk::ASCII) {
     if (datatype_ == Vtk::FLOAT32)
       out << std::setprecision(std::numeric_limits<float>::digits10+2);
@@ -105,6 +107,7 @@ void VtkImageDataWriter<GV,DC>
 {
   std::string filename = pfilename + ".p" + this->fileExtension();
   std::ofstream out(filename, std::ios_base::ate | std::ios::binary);
+  assert(out.is_open());
 
   out << "<VTKFile"
       << " type=\"StructuredGrid\""
@@ -118,7 +121,7 @@ void VtkImageDataWriter<GV,DC>
   auto const& origin = dataCollector_.origin();
   auto const& spacing = dataCollector_.spacing();
   out << "<PImageData"
-      << " GhostLevel=\"0\""
+      << " GhostLevel=\"" << dataCollector_.ghostLevel() << "\""
       << " WholeExtent=\"" << join(wholeExtent.begin(), wholeExtent.end()) << "\""
       << " Origin=\"" << join(origin.begin(), origin.end()) << "\""
       << " Spacing=\"" << join(spacing.begin(), spacing.end()) << "\""

dune/vtk/writers/vtkrectilineargridwriter.hh

+#pragma once
+
+#include <array>
+#include <iosfwd>
+#include <map>
+
+#include <dune/vtk/datacollector.hh>
+#include <dune/vtk/filewriter.hh>
+#include <dune/vtk/vtkfunction.hh>
+#include <dune/vtk/vtktypes.hh>
+#include <dune/vtk/datacollectors/structureddatacollector.hh>
+
+#include "vtkwriterinterface.hh"
+
+namespace Dune { namespace experimental
+{
+  /// File-Writer for VTK .vtu files
+  template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
+  class VtkRectilinearGridWriter
+      : public VtkWriterInterface<GridView, DataCollector>
+  {
+    static constexpr int dimension = GridView::dimension;
+
+    using Super = VtkWriterInterface<GridView, DataCollector>;
+    using pos_type = typename Super::pos_type;
+
+  public:
+    /// Constructor, stores the gridView
+    VtkRectilinearGridWriter (GridView const& gridView)
+      : Super(gridView)
+    {}
+
+  private:
+    /// Write a serial VTK file in Unstructured format
+    virtual void writeSerialFile (std::string const& filename) const override;
+
+    /// Write a parallel VTK file `pfilename.pvtu` in Unstructured format,
+    /// with `size` the number of pieces and serial files given by `pfilename_p[i].vtu`
+    /// for [i] in [0,...,size).
+    virtual void writeParallelFile (std::string const& pfilename, int size) const override;
+
+    void writeCoordinates (std::ofstream& out, std::vector<pos_type>& offsets) const;
+
+    template <class T>
+    std::array<std::uint64_t, 3> writeCoordinatesAppended (std::ofstream& out) const;
+
+    virtual std::string fileExtension () const override
+    {
+      return "vtr";
+    }
+
+  private:
+    using Super::dataCollector_;
+    using Super::format_;
+    using Super::datatype_;
+
+    // attached data
+    using Super::pointData_;
+    using Super::cellData_;
+  };
+
+}} // end namespace Dune::experimental
+
+#include "vtkrectilineargridwriter.impl.hh"

dune/vtk/writers/vtkrectilineargridwriter.impl.hh

+#pragma once
+
+#include <iomanip>
+#include <iostream>
+#include <iterator>
+#include <fstream>
+#include <sstream>
+#include <string>
+
+#include <dune/geometry/referenceelements.hh>
+#include <dune/geometry/type.hh>
+
+#include <dune/vtk/utility/enum.hh>
+#include <dune/vtk/utility/filesystem.hh>
+#include <dune/vtk/utility/string.hh>
+
+namespace Dune { namespace experimental {
+
+template <class GV, class DC>
+void VtkRectilinearGridWriter<GV,DC>
+  ::writeSerialFile (std::string const& filename) const
+{
+  std::ofstream out(filename, std::ios_base::ate | std::ios::binary);
+  assert(out.is_open());
+  if (format_ == Vtk::ASCII) {
+    if (datatype_ == Vtk::FLOAT32)
+      out << std::setprecision(std::numeric_limits<float>::digits10+2);
+    else
+      out << std::setprecision(std::numeric_limits<double>::digits10+2);
+  }
+
+  std::vector<pos_type> offsets; // pos => offset
+  out << "<VTKFile"
+      << " type=\"RectilinearGrid\""
+      << " version=\"1.0\""
+      << " byte_order=\"" << this->getEndian() << "\""
+      << " header_type=\"UInt64\""
+      << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\"" : "")
+      << ">\n";
+
+  auto const& wholeExtent = dataCollector_.wholeExtent();
+  auto const& origin = dataCollector_.origin();
+  auto const& spacing = dataCollector_.spacing();
+  out << "<RectilinearGrid"
+      << " WholeExtent=\"" << join(wholeExtent.begin(), wholeExtent.end()) << "\""
+      << ">\n";
+
+  dataCollector_.writeLocalPiece([&out](auto const& extent) {
+    out << "<Piece Extent=\"" << join(extent.begin(), extent.end()) << "\">\n";
+  });
+
+  // Write data associated with grid points
+  out << "<PointData" << this->getNames(pointData_) << ">\n";
+  for (auto const& v : pointData_)
+    this->writeData(out, offsets, v, Super::POINT_DATA);
+  out << "</PointData>\n";
+
+  // Write data associated with grid cells
+  out << "<CellData" << this->getNames(cellData_) << ">\n";
+  for (auto const& v : cellData_)
+    this->writeData(out, offsets, v, Super::CELL_DATA);
+  out << "</CellData>\n";
+
+  // Write point coordinates for x, y, and z ordinate
+  out << "<Coordinates>\n";
+  writeCoordinates(out, offsets);
+  out << "</Coordinates>\n";
+
+  out << "</Piece>\n";
+  out << "</RectilinearGrid>\n";
+
+  std::vector<std::uint64_t> blocks; // size of i'th appended block
+  pos_type appended_pos = 0;
+  if (is_a(format_, Vtk::APPENDED)) {
+    out << "<AppendedData encoding=\"raw\">\n_";
+    appended_pos = out.tellp();
+    for (auto const& v : pointData_) {
+      if (v.type() == Vtk::FLOAT32)
+        blocks.push_back( this->template writeDataAppended<float>(out, v, Super::POINT_DATA) );
+      else
+        blocks.push_back( this->template writeDataAppended<double>(out, v, Super::POINT_DATA) );
+    }
+    for (auto const& v : cellData_) {
+      if (v.type() == Vtk::FLOAT32)
+        blocks.push_back( this->template writeDataAppended<float>(out, v, Super::CELL_DATA) );
+      else
+        blocks.push_back( this->template writeDataAppended<double>(out, v, Super::CELL_DATA) );
+    }
+
+    if (datatype_ == Vtk::FLOAT32) {
+      auto bs = writeCoordinatesAppended<float>(out);
+      blocks.insert(blocks.end(), bs.begin(), bs.end());
+    } else {
+      auto bs = writeCoordinatesAppended<double>(out);
+      blocks.insert(blocks.end(), bs.begin(), bs.end());
+    }
+    out << "</AppendedData>\n";
+  }
+
+  out << "</VTKFile>";
+
+  // fillin offset values and block sizes
+  if (is_a(format_, Vtk::APPENDED)) {
+    pos_type offset = 0;
+    for (std::size_t i = 0; i < offsets.size(); ++i) {
+      out.seekp(offsets[i]);
+      out << '"' << offset << '"';
+      offset += pos_type(blocks[i]);
+    }
+  }
+}
+