added structured grid writer

dune/vtk/datacollectors/structureddatacollector.hh

+#pragma once
+
+#include <dune/vtk/datacollector.hh>
+
+namespace Dune { namespace experimental
+{
+
+/// Implementation of \ref DataCollector for linear cells, with continuous data.
+template <class GridView>
+class StructuredDataCollector
+    : public DataCollectorInterface<GridView, StructuredDataCollector<GridView>>
+{
+  enum { dim = GridView::dimension };
+
+  using Self = StructuredDataCollector;
+  using Super = DataCollectorInterface<GridView, Self>;
+  using Super::gridView_;
+  using ctype = typename GridView::ctype;
+
+  struct CoordLess
+  {
+    template <class T, int N>
+    bool operator() (FieldVector<T,N> const& lhs, FieldVector<T,N> const& rhs) const
+    {
+      for (int i = N-1; i >= 0; --i) {
+        if (std::abs(lhs[i] - rhs[i]) < std::numeric_limits<T>::epsilon())
+          continue;
+        return lhs[i] < rhs[i];
+      }
+      return false;
+    }
+  };
+
+public:
+  StructuredDataCollector (GridView const& gridView)
+    : Super(gridView)
+  {}
+
+  /// Return number of grid vertices
+  std::uint64_t numPointsImpl () const
+  {
+    return gridView_.size(dim);
+  }
+
+  std::array<int, 6> const& extent () const
+  {
+    return extent_;
+  }
+
+  auto const& origin () const
+  {
+    return origin_;
+  }
+
+  auto const& spacing () const
+  {
+    return spacing_;
+  }
+
+  void updateImpl ()
+  {
+    // TODO: extract this information from the grid
+    int extent = GridView::dimension == 1 ? gridView_.size(dim) :
+                 GridView::dimension == 2 ? isqrt(gridView_.size(dim)) :
+                 GridView::dimension == 3 ? icbrt(gridView_.size(dim)) : 0;
+    for (int i = 0; i < 3; ++i) {
+      if (GridView::dimension > i) {
+        extent_[2*i] = 0;
+        extent_[2*i+1] = extent-1;
+        spacing_[i] = gridView_.grid().domainSize()[i] / (extent-1);
+      } else {
+        extent_[2*i] = extent_[2*i+1] = 0;
+        spacing_[i] = 0;
+      }
+
+      origin_[i] = 0;
+    }
+  }
+
+  /// Return the coordinates of all grid vertices in the order given by the indexSet
+  template <class T>
+  std::vector<T> pointsImpl () const
+  {
+    std::vector<T> data(this->numPoints() * 3);
+    auto const& indexSet = gridView_.indexSet();
+    for (auto const& vertex : vertices(gridView_)) {
+      std::size_t idx = 3 * indexSet.index(vertex);
+      auto v = vertex.geometry().center();
+      for (std::size_t j = 0; j < v.size(); ++j)
+        data[idx + j] = T(v[j]);
+      for (std::size_t j = v.size(); j < 3u; ++j)
+        data[idx + j] = T(0);
+    }
+    return data;
+  }
+
+  /// Evaluate the `fct` at the corners of the elements
+  template <class T, class GlobalFunction>
+  std::vector<T> pointDataImpl (GlobalFunction const& fct) const
+  {
+    std::vector<T> data(this->numPoints() * fct.ncomps());
+    auto const& indexSet = gridView_.indexSet();
+    auto localFct = localFunction(fct);
+    for (auto const& e : elements(gridView_)) {
+      localFct.bind(e);
+      Vtk::CellType cellType{e.type()};
+      auto refElem = referenceElement(e.geometry());
+      for (int j = 0; j < e.subEntities(dim); ++j) {
+        int k = cellType.permutation(j);
+        std::size_t idx = fct.ncomps() * indexSet.subIndex(e,k,dim);
+        for (int comp = 0; comp < fct.ncomps(); ++comp)
+          data[idx + comp] = T(localFct.evaluate(comp, refElem.position(k,dim)));
+      }
+      localFct.unbind();
+    }
+    return data;
+  }
+
+private:
+  static constexpr std::uint32_t isqrt (std::uint64_t x)
+  {
+    std::uint32_t c = 0x8000;
+    std::uint32_t g = 0x8000;
+
+    while (true) {
+      if (g*g > x)
+        g ^= c;
+      c >>= 1;
+      if (c == 0)
+        return g;
+      g |= c;
+    }
+  }
+
+  static constexpr std::uint32_t icbrt (std::uint64_t x)
+  {
+    std::uint32_t y = 0;
+    for (int s = 63; s >= 0; s -= 3) {
+      y += y;
+      std::uint64_t b = 3*y*(std::uint64_t(y) + 1) + 1;
+      if ((x >> s) >= b) {
+        x -= b << s;
+        y++;
+      }
+    }
+    return y;
+  }
+
+private:
+  std::array<int, 6> extent_;
+  FieldVector<ctype,3> spacing_;
+  FieldVector<ctype,3> origin_;
+  std::vector<std::size_t> indexMap_;
+};
+
+}} // end namespace Dune::experimental

dune/vtk/vtkimagedatawriter.hh

+#pragma once
+
+#include <array>
+#include <iosfwd>
+#include <map>
+
+#include "datacollector.hh"
+#include "filewriter.hh"
+#include "vtkfunction.hh"
+#include "vtktypes.hh"
+#include "vtkwriter.hh"
+#include "datacollectors/structureddatacollector.hh"
+
+namespace Dune { namespace experimental
+{
+  /// File-Writer for VTK .vtu files
+  template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
+  class VtkImageDataWriter
+      : public VtkWriter<GridView, DataCollector>
+  {
+    static constexpr int dimension = GridView::dimension;
+
+    using Super = VtkWriter<GridView, DataCollector>;
+    using pos_type = typename Super::pos_type;
+
+  public:
+    /// Constructor, stores the gridView
+    VtkImageDataWriter (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;
+
+    virtual std::string fileExtension () const override
+    {
+      return "vti";
+    }
+
+  private:
+    using Super::dataCollector_;
+    using Super::format_;
+    using Super::datatype_;
+
+    // attached data
+    using Super::pointData_;
+    using Super::cellData_;
+  };
+
+}} // end namespace Dune::experimental
+
+#include "vtkimagedatawriter.impl.hh"

dune/vtk/vtkimagedatawriter.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 "utility/enum.hh"
+#include "utility/filesystem.hh"
+#include "utility/string.hh"
+
+namespace Dune { namespace experimental {
+
+template <class GV, class DC>
+void VtkImageDataWriter<GV,DC>
+  ::writeSerialFile (std::string const& filename) const
+{
+  std::ofstream out(filename, std::ios_base::ate | std::ios::binary);
+  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);
+  }
+
+  dataCollector_.update();
+
+  std::vector<pos_type> offsets; // pos => offset
+  out << "<VTKFile type=\"ImageData\" version=\"1.0\" "
+      << "byte_order=\"" << this->getEndian() << "\" header_type=\"UInt64\""
+      << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\">\n" : ">\n");
+  out << "<ImageData WholeExtent=\"";
+  auto const& extent = dataCollector_.extent();
+  for (int i = 0; i < 3; ++i) {
+    out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
+  }
+  out << "\" Origin=\"";
+  auto const& origin = dataCollector_.origin();
+  for (int i = 0; i < 3; ++i) {
+    out << (i == 0 ? "" : " ") << origin[i];
+  }
+  out << "\" Spacing=\"";
+  auto const& spacing = dataCollector_.spacing();
+  for (int i = 0; i < 3; ++i) {
+    out << (i == 0 ? "" : " ") << spacing[i];
+  }
+  out << "\">\n";
+  out << "<Piece Extent=\"";
+  for (int i = 0; i < 3; ++i) {
+    out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
+  }
+  out << "\">\n";
+
+  { // Write data associated with grid points
+    auto scalar = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+    auto vector = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+    out << "<PointData" << (scalar != pointData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+                        << (vector != pointData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+                        << ">\n";
+    for (auto const& v : pointData_)
+      this->writeData(out, offsets, v, Super::POINT_DATA);
+    out << "</PointData>\n";
+  }
+
+  { // Write data associated with grid cells
+    auto scalar = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+    auto vector = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+    out << "<CellData" << (scalar != cellData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+                       << (vector != cellData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+                       << ">\n";
+    for (auto const& v : cellData_)
+      this->writeData(out, offsets, v, Super::CELL_DATA);
+    out << "</CellData>\n";
+  }
+  out << "</Piece>\n";
+  out << "</ImageData>\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 (datatype_ == 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 (datatype_ == 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) );
+    }
+    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]);
+    }
+  }
+}
+
+
+template <class GV, class DC>
+void VtkImageDataWriter<GV,DC>
+  ::writeParallelFile (std::string const& pfilename, int size) const
+{
+  std::string filename = pfilename + ".p" + this->fileExtension();
+  std::ofstream out(filename, std::ios_base::ate | std::ios::binary);
+
+  out << "<VTKFile type=\"PImageData\" version=\"1.0\" "
+      << "byte_order=\"" << this->getEndian() << "\" header_type=\"UInt64\""
+      << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\">\n" : ">\n");
+  out << "<PImageData GhostLevel=\"0\" WholeExtent=\"";
+  auto const& extent = dataCollector_.extent();
+  for (int i = 0; i < 3; ++i) {
+    out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
+  }
+  out << "\" Origin=\"";
+  auto const& origin = dataCollector_.origin();
+  for (int i = 0; i < 3; ++i) {
+    out << (i == 0 ? "" : " ") << origin[i];
+  }
+  out << "\" Spacing=\"";
+  auto const& spacing = dataCollector_.spacing();
+  for (int i = 0; i < 3; ++i) {
+    out << (i == 0 ? "" : " ") << spacing[i];
+  }
+  out << "\">\n";
+
+  { // Write data associated with grid points
+    auto scalar = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+    auto vector = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+    out << "<PPointData" << (scalar != pointData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+                         << (vector != pointData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+                         << ">\n";
+    for (auto const& v : pointData_) {
+      out << "<PDataArray Name=\"" << v.name() << "\" type=\"" << Vtk::Map::from_datatype[datatype_] << "\""
+          << " NumberOfComponents=\"" << v.ncomps() << "\" />\n";
+    }
+    out << "</PPointData>\n";
+  }
+
+  { // Write data associated with grid cells
+    auto scalar = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+    auto vector = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+    out << "<PCellData" << (scalar != cellData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+                        << (vector != cellData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+                        << ">\n";
+    for (auto const& v : cellData_) {
+      out << "<PDataArray Name=\"" << v.name() << "\" type=\"" << Vtk::Map::from_datatype[datatype_] << "\""
+          << " NumberOfComponents=\"" << v.ncomps() << "\" />\n";
+    }
+    out << "</PCellData>\n";
+  }
+
+  // Write piece file references
+  for (int i = 0; i < size; ++i) {
+    out << "<Piece Source=\"" << pfilename << "_p" << std::to_string(i) << "." << this->fileExtension() << "\" Extent=\"";
+    for (int i = 0; i < 3; ++i) {
+      out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
+    }
+    out << "\" />\n";
+  }
+
+  out << "</PUnstructuredGrid>\n";
+  out << "</VTKFile>";
+}
+
+}} // end namespace Dune::experimental

dune/vtk/vtkstructuredgridwriter.hh

+#pragma once
+
+#include <array>
+#include <iosfwd>
+#include <map>
+
+#include "datacollector.hh"
+#include "filewriter.hh"
+#include "vtkfunction.hh"
+#include "vtktypes.hh"
+#include "vtkwriter.hh"
+#include "datacollectors/structureddatacollector.hh"
+
+namespace Dune { namespace experimental
+{
+  /// File-Writer for VTK .vtu files
+  template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
+  class VtkStructuredGridWriter
+      : public VtkWriter<GridView, DataCollector>
+  {
+    static constexpr int dimension = GridView::dimension;
+
+    using Super = VtkWriter<GridView, DataCollector>;
+    using pos_type = typename Super::pos_type;
+
+  public:
+    /// Constructor, stores the gridView
+    VtkStructuredGridWriter (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;
+
+    virtual std::string fileExtension () const override
+    {
+      return "vts";
+    }
+
+  private:
+    using Super::dataCollector_;
+    using Super::format_;
+    using Super::datatype_;
+
+    // attached data
+    using Super::pointData_;
+    using Super::cellData_;
+  };
+
+}} // end namespace Dune::experimental
+
+#include "vtkstructuredgridwriter.impl.hh"

dune/vtk/vtkstructuredgridwriter.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 "utility/enum.hh"
+#include "utility/filesystem.hh"
+#include "utility/string.hh"
+
+namespace Dune { namespace experimental {
+
+template <class GV, class DC>
+void VtkStructuredGridWriter<GV,DC>
+  ::writeSerialFile (std::string const& filename) const
+{
+  std::ofstream out(filename, std::ios_base::ate | std::ios::binary);
+  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);
+  }
+
+  dataCollector_.update();
+
+  std::vector<pos_type> offsets; // pos => offset
+  out << "<VTKFile type=\"StructuredGrid\" version=\"1.0\" "
+      << "byte_order=\"" << this->getEndian() << "\" header_type=\"UInt64\""
+      << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\">\n" : ">\n");
+  out << "<StructuredGrid WholeExtent=\"";
+  for (int i = 0; i < 3; ++i) {
+    out << (i == 0 ? "" : " ") << dataCollector_.extent()[2*i];
+    out << " " << dataCollector_.extent()[2*i+1];
+  }
+  out << "\">\n";
+  out << "<Piece NumberOfPoints=\"" << dataCollector_.numPoints() << "\" "
+      << "NumberOfCells=\"" << dataCollector_.numCells() << "\" Extent=\"";
+  for (int i = 0; i < 3; ++i) {
+    out << (i == 0 ? "" : " ") << dataCollector_.extent()[2*i];
+    out << " " << dataCollector_.extent()[2*i+1];
+  }
+  out << "\">\n";
+
+  { // Write data associated with grid points
+    auto scalar = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+    auto vector = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+    out << "<PointData" << (scalar != pointData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+                        << (vector != pointData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+                        << ">\n";
+    for (auto const& v : pointData_)
+      this->writeData(out, offsets, v, Super::POINT_DATA);
+    out << "</PointData>\n";
+  }
+
+  { // Write data associated with grid cells
+    auto scalar = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+    auto vector = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+    out << "<CellData" << (scalar != cellData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+                       << (vector != cellData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+                       << ">\n";
+    for (auto const& v : cellData_)
+      this->writeData(out, offsets, v, Super::CELL_DATA);
+    out << "</CellData>\n";
+  }
+
+  // Write point coordinates
+  out << "<Points>\n";
+  this->writePoints(out, offsets);
+  out << "</Points>\n";
+
+  out << "</Piece>\n";
+  out << "</StructuredGrid>\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 (datatype_ == 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 (datatype_ == 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)
+      blocks.push_back( this->template writePointsAppended<float>(out) );
+    else
+      blocks.push_back( this->template writePointsAppended<double>(out) );
+    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]);
+    }
+  }
+}
+
+
+template <class GV, class DC>
+void VtkStructuredGridWriter<GV,DC>
+  ::writeParallelFile (std::string const& pfilename, int size) const
+{
+  std::string filename = pfilename + ".p" + this->fileExtension();
+  std::ofstream out(filename, std::ios_base::ate | std::ios::binary);