Merge branch 'feature/timeseries' into 'master'

Feature/timeseries

See merge request spraetor/dune-vtk!4

dune/vtk/utility/uid.hh

+#pragma once
+
+#include <cstdlib>
+#include <cstring>
+#include <ctime>
+#include <string>
+
+namespace Dune
+{
+  inline std::string uid (std::size_t len = 8)
+  {
+    static const auto digits = "0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ";
+    static const int N = std::strlen(digits);
+
+    std::string id(len,' ');
+    for (std::size_t i = 0; i < len; ++i)
+      id[i] = digits[std::rand()%N];
+
+    return id;
+  }
+
+} // end namespace Dune

dune/vtk/vtktimeserieswriter.hh

+#pragma once
+
+#include <string>
+#include <tuple>
+#include <vector>
+
+#include <dune/vtk/filewriter.hh>
+#include <dune/vtk/vtktypes.hh>
+#include <dune/vtk/utility/filesystem.hh>
+#include <dune/vtk/utility/uid.hh>
+
+namespace Dune
+{
+  /// File-Writer for Vtk timeseries .vtu files
+  /**
+   * \tparam VtkWriter  Type of a FileWriter derived from \ref VtkWriterInterface that
+   *                    additionally supports writeTimeseriesSerialFile() and writeTimeseriesParallelFile(),
+   *                    e.g. \ref VtkUnstructuredGridWriter.
+   **/
+  template <class VtkWriter>
+  class VtkTimeseriesWriter
+      : public FileWriter
+  {
+  protected:
+    using Self = VtkTimeseriesWriter;
+    using pos_type = typename std::ostream::pos_type;
+
+    template <class W>
+    using HasWriteTimeseriesSerialFile = decltype(&W::writeTimeseriesSerialFile);
+    static_assert(Std::is_detected<HasWriteTimeseriesSerialFile, VtkWriter>::value, "");
+
+    template <class W>
+    using HasWriteTimeseriesParallelFile = decltype(&W::writeTimeseriesParallelFile);
+    static_assert(Std::is_detected<HasWriteTimeseriesParallelFile, VtkWriter>::value, "");
+
+  public:
+    /// Constructor, stores the gridView
+    template <class... Args, disableCopyMove<Self, Args...> = 0>
+    VtkTimeseriesWriter (Args&&... args)
+      : vtkWriter_{std::forward<Args>(args)...}
+    {
+      assert(vtkWriter_.format_ != Vtk::ASCII && "Timeseries writer requires APPENDED mode");
+      std::srand(std::time(nullptr));
+      // put temporary file to /tmp directory
+      tmpDir_ = filesystem::path("/tmp/vtktimeserieswriter_" + uid(10) + "/");
+      assert( filesystem::exists("/tmp") );
+      filesystem::create_directories(tmpDir_);
+    }
+
+    ~VtkTimeseriesWriter ()
+    {
+      std::remove(tmpDir_.string().c_str());
+    }
+
+    /// Write the attached data to the file with \ref Vtk::FormatTypes and \ref Vtk::DataTypes
+    void writeTimestep (double time, std::string const& fn);
+
+    /// Writes all timesteps to single timeseries file.
+    // NOTE: requires an aforging call to \ref writeTimestep
+    virtual void write (std::string const& fn) override;
+
+    /// Attach point data to the writer, \see VtkFunction for possible arguments
+    template <class Function, class... Args>
+    VtkTimeseriesWriter& addPointData (Function const& fct, Args&&... args)
+    {
+      vtkWriter_.addPointData(fct, std::forward<Args>(args)...);
+      return *this;
+    }
+
+    /// Attach cell data to the writer, \see VtkFunction for possible arguments
+    template <class Function, class... Args>
+    VtkTimeseriesWriter& addCellData (Function const& fct, Args&&... args)
+    {
+      vtkWriter_.addCellData(fct, std::forward<Args>(args)...);
+      return *this;
+    }
+
+  protected:
+    VtkWriter vtkWriter_;
+    filesystem::path tmpDir_;
+
+    bool initialized_ = false;
+
+    // block size of attached data
+    std::vector<std::uint64_t> blocks_;
+
+    std::string filenameMesh_;
+    std::vector<std::pair<double, std::string>> timesteps_;
+  };
+
+} // end namespace Dune
+
+#include "vtktimeserieswriter.impl.hh"

dune/vtk/vtktimeserieswriter.impl.hh

+#pragma once
+
+#include <algorithm>
+#include <cstdio>
+#include <iomanip>
+#include <iostream>
+#include <iterator>
+#include <fstream>
+#include <sstream>
+#include <string>
+
+#ifdef HAVE_ZLIB
+#include <zlib.h>
+#endif
+
+#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 {
+
+template <class W>
+void VtkTimeseriesWriter<W>
+  ::writeTimestep (double time, std::string const& fn)
+{
+  auto name = filesystem::path(fn).stem();
+  auto tmp = tmpDir_;
+  tmp /= name.string();
+
+  vtkWriter_.dataCollector_.update();
+
+  std::string filenameBase = tmp.string();
+
+  int rank = 0;
+  int num_ranks = 1;
+  #ifdef HAVE_MPI
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    MPI_Comm_size(MPI_COMM_WORLD, &num_ranks);
+    if (num_ranks > 1)
+      filenameBase = tmp.string() + "_p" + std::to_string(rank);
+  #endif
+
+  if (!initialized_) {
+    // write points and cells only once
+    filenameMesh_ = filenameBase + ".mesh.vtkdata";
+    std::ofstream out(filenameMesh_, std::ios_base::ate | std::ios::binary);
+    vtkWriter_.writeGridAppended(out, blocks_);
+    initialized_ = true;
+  }
+
+  std::string filenameData = filenameBase + "_t" + std::to_string(timesteps_.size()) + ".vtkdata";
+  std::ofstream out(filenameData, std::ios_base::ate | std::ios::binary);
+  vtkWriter_.writeDataAppended(out, blocks_);
+  timesteps_.emplace_back(time, filenameData);
+}
+
+
+template <class W>
+void VtkTimeseriesWriter<W>
+  ::write (std::string const& fn)
+{
+  assert( initialized_ );
+  assert( timesteps_.size() < 1000 ); // maximal number of allowed timesteps in timeseries file
+
+  auto p = filesystem::path(fn);
+  auto name = p.stem();
+  p.remove_filename();
+  p /= name.string();
+
+  std::string filenameBase = p.string() + "_ts";
+
+  int rank = 0;
+  int num_ranks = 1;
+#ifdef HAVE_MPI
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    MPI_Comm_size(MPI_COMM_WORLD, &num_ranks);
+    if (num_ranks > 1)
+      filenameBase = p.string() + "_ts_p" + std::to_string(rank);
+#endif
+
+  std::string filename = filenameBase + "." + vtkWriter_.getFileExtension();
+  vtkWriter_.writeTimeseriesSerialFile(filename, filenameMesh_, timesteps_, blocks_);
+
+#ifdef HAVE_MPI
+  if (num_ranks > 1 && rank == 0)
+    vtkWriter_.writeTimeseriesParallelFile(p.string() + "_ts", num_ranks, timesteps_);
+#endif
+
+  // remove all temporary data files
+  int ec = std::remove(filenameMesh_.c_str());
+  assert(ec == 0);
+  for (auto const& timestep : timesteps_) {
+    ec = std::remove(timestep.second.c_str());
+    assert(ec == 0);
+  }
+}
+
+} // end namespace Dune

dune/vtk/vtkwriterinterface.hh

@@ -5,6 +5,10 @@
 #include <string>
 #include <vector>
 
+#ifdef HAVE_ZLIB
+#include <zlib.h>
+#endif
+
 #include <dune/common/std/optional.hh>
 
 #include <dune/vtk/filewriter.hh>
@@ -13,11 +17,17 @@
 
 namespace Dune
 {
+  // forward declaration
+  template <class VtkWriter>
+  class VtkTimeseriesWriter;
+
   /// File-Writer for Vtk .vtu files
   template <class GridView, class DataCollector>
   class VtkWriterInterface
       : public FileWriter
   {
+    template <class> friend class VtkTimeseriesWriter;
+
   protected:
     static constexpr int dimension = GridView::dimension;
 
@@ -31,20 +41,23 @@ namespace Dune
 
   public:
     /// Constructor, stores the gridView
-    VtkWriterInterface (GridView const& gridView)
+    VtkWriterInterface (GridView const& gridView,
+                        Vtk::FormatTypes format = Vtk::BINARY,
+                        Vtk::DataTypes datatype = Vtk::FLOAT32)
       : dataCollector_(gridView)
-    {}
-
-    /// Write the attached data to the file
-    virtual void write (std::string const& fn) override
+      , format_(format)
+      , datatype_(datatype)
     {
-      write(fn, Vtk::BINARY);
+#ifndef HAVE_ZLIB
+      if (format_ == Vtk::COMPRESSED) {
+        std::cout << "Dune is compiled without compression. Falling back to BINARY VTK output!\n";
+        format_ = Vtk::BINARY;
+      }
+#endif
     }
 
-    /// 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);
+    /// Write the attached data to the file
+    virtual void write (std::string const& fn) override;
 
     /// Attach point data to the writer, \see VtkFunction for possible arguments
     template <class Function, class... Args>
@@ -74,35 +87,35 @@ namespace Dune
     /// Return the file extension of the serial file (not including the dot)
     virtual std::string fileExtension () const = 0;
 
+    /// Write points and cells in raw/compressed format to output stream
+    virtual void writeGridAppended (std::ofstream& out, std::vector<std::uint64_t>& blocks) 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,
                     VtkFunction const& fct,
-                    PositionTypes type) const;
+                    PositionTypes type,
+                    Std::optional<std::size_t> timestep = {}) 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,
-                                     VtkFunction const& fct,
-                                     PositionTypes type) const;
+    // Write points-data and cell-data in raw/compressed format to output stream
+    void writeDataAppended (std::ofstream& out, std::vector<std::uint64_t>& blocks) 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;
+                      std::vector<pos_type>& offsets,
+                      Std::optional<std::size_t> timestep = {}) const;
 
-    // Collect point positions and pass the resulting vector to \ref writeAppended.
-    template <class T>
-    std::uint64_t writePointsAppended (std::ofstream& out) const;
+    // Write Appended section and fillin offset values to XML attributes
+    void writeAppended (std::ofstream& out, std::vector<pos_type> const& offsets) 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::uint64_t writeValuesAppended (std::ofstream& out, std::vector<T> const& values) const;
 
     /// Return PointData/CellData attributes for the name of the first scalar/vector/tensor DataArray
     std::string getNames (std::vector<VtkFunction> const& data) const;
@@ -114,6 +127,11 @@ namespace Dune
       return (reinterpret_cast<char*>(&i)[1] == 1 ? "BigEndian" : "LittleEndian");
     }
 
+    std::string getFileExtension () const
+    {
+      return fileExtension();
+    }
+
   protected:
     mutable DataCollector dataCollector_;
 

dune/vtk/vtkwriterinterface.impl.hh

@@ -8,10 +8,6 @@
 #include <sstream>
 #include <string>
 
-#ifdef HAVE_ZLIB
-#include <zlib.h>
-#endif
-
 #include <dune/geometry/referenceelements.hh>
 #include <dune/geometry/type.hh>
 
@@ -23,18 +19,8 @@ namespace Dune {
 
 template <class GV, class DC>
 void VtkWriterInterface<GV,DC>
-  ::write (std::string const& fn, Vtk::FormatTypes format, Vtk::DataTypes datatype)
+  ::write (std::string const& fn)
 {
-  format_ = format;
-  datatype_ = datatype;
-
-#ifndef HAVE_ZLIB
-  if (format_ == Vtk::COMPRESSED) {
-    std::cout << "Dune is compiled without compression. Falling back to BINARY VTK output!\n";
-    format_ = Vtk::BINARY;
-  }
-#endif
-
   dataCollector_.update();
 
   auto p = filesystem::path(fn);
@@ -44,21 +30,21 @@ void VtkWriterInterface<GV,DC>
 
   std::string filename = p.string() + "." + fileExtension();
 
+  int rank = 0;
+  int num_ranks = 1;
+
 #ifdef HAVE_MPI
-  int rank = -1;
-  int num_ranks = -1;
   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
   MPI_Comm_size(MPI_COMM_WORLD, &num_ranks);
-  if (num_ranks > 1) {
+  if (num_ranks > 1)
     filename = p.string() + "_p" + std::to_string(rank) + "." + fileExtension();
+#endif
 
-    writeSerialFile(filename);
-    if (rank == 0) {
-      writeParallelFile(p.string(), num_ranks);
-    }
-  } else {
-    writeSerialFile(filename);
-  }
+  writeSerialFile(filename);
+
+#ifdef HAVE_MPI
+  if (num_ranks > 1 && rank == 0)
+    writeParallelFile(p.string(), num_ranks);
 #endif
 }
 
@@ -66,12 +52,16 @@ void VtkWriterInterface<GV,DC>
 template <class GV, class DC>
 void VtkWriterInterface<GV,DC>
   ::writeData (std::ofstream& out, std::vector<pos_type>& offsets,
-               VtkFunction const& fct, PositionTypes type) const
+               VtkFunction const& fct, PositionTypes type,
+               Std::optional<std::size_t> timestep) const
 {
   out << "<DataArray Name=\"" << fct.name() << "\" type=\"" << to_string(fct.type()) << "\""
-      << " NumberOfComponents=\"" << fct.ncomps() << "\" format=\"" << (format_ == Vtk::ASCII ? "ascii\">\n" : "appended\"");
+      << " NumberOfComponents=\"" << fct.ncomps() << "\" format=\"" << (format_ == Vtk::ASCII ? "ascii\"" : "appended\"");
+  if (timestep)
+    out << " TimeStep=\"" << *timestep << "\"";
 
   if (format_ == Vtk::ASCII) {
+    out << ">\n";
     std::size_t i = 0;
     if (type == POINT_DATA) {
       auto data = dataCollector_.template pointData<double>(fct);
@@ -94,12 +84,16 @@ void VtkWriterInterface<GV,DC>
 
 template <class GV, class DC>
 void VtkWriterInterface<GV,DC>
-  ::writePoints (std::ofstream& out, std::vector<pos_type>& offsets) const
+  ::writePoints (std::ofstream& out, std::vector<pos_type>& offsets,
+                Std::optional<std::size_t> timestep) const
 {
   out << "<DataArray type=\"" << to_string(datatype_) << "\""
-      << " NumberOfComponents=\"3\" format=\"" << (format_ == Vtk::ASCII ? "ascii\">\n" : "appended\"");
+      << " NumberOfComponents=\"3\" format=\"" << (format_ == Vtk::ASCII ? "ascii\"" : "appended\"");
+  if (timestep)
+    out << " TimeStep=\"" << *timestep << "\"";
 
   if (format_ == Vtk::ASCII) {
+    out << ">\n";
     auto points = dataCollector_.template points<double>();
     std::size_t i = 0;
     for (auto const& v : points)
@@ -113,33 +107,44 @@ void VtkWriterInterface<GV,DC>
   }
 }
 
-
 template <class GV, class DC>
-  template <class T>
-std::uint64_t VtkWriterInterface<GV,DC>
-  ::writeDataAppended (std::ofstream& out, VtkFunction const& fct, PositionTypes type) const
+void VtkWriterInterface<GV,DC>
+  ::writeDataAppended (std::ofstream& out, std::vector<std::uint64_t>& blocks) const
 {
-  assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
-
-  if (type == POINT_DATA) {
-    auto data = dataCollector_.template pointData<T>(fct);
-    return this->writeAppended(out, data);
-  } else {
-    auto data = dataCollector_.template cellData<T>(fct);
-    return this->writeAppended(out, data);
+  for (auto const& v : pointData_) {
+    blocks.push_back( v.type() == Vtk::FLOAT32
+      ? this->writeValuesAppended(out, dataCollector_.template pointData<float>(v))
+      : this->writeValuesAppended(out, dataCollector_.template pointData<double>(v)));
+  }
+  for (auto const& v : cellData_) {
+    blocks.push_back( v.type() == Vtk::FLOAT32
+      ? this->writeValuesAppended(out, dataCollector_.template cellData<float>(v))
+      : this->writeValuesAppended(out, dataCollector_.template cellData<double>(v)));
   }
 }
 
 
 template <class GV, class DC>
-  template <class T>
-std::uint64_t VtkWriterInterface<GV,DC>
-  ::writePointsAppended (std::ofstream& out) const
+void VtkWriterInterface<GV,DC>
+  ::writeAppended (std::ofstream& out, std::vector<pos_type> const& offsets) const
 {
-  assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
+  if (is_a(format_, Vtk::APPENDED)) {
+    out << "<AppendedData encoding=\"raw\">\n_";
+    std::vector<std::uint64_t> blocks;
+    writeGridAppended(out, blocks);
+    writeDataAppended(out, blocks);
+    out << "</AppendedData>\n";
+    pos_type appended_pos = out.tellp();
+
+    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]);
+    }
 
-  auto points = dataCollector_.template points<T>();
-  return this->writeAppended(out, points);
+    out.seekp(appended_pos);
+  }
 }
 
 
@@ -187,7 +192,7 @@ std::uint64_t writeCompressed (unsigned char const* buffer, unsigned char* buffe
 template <class GV, class DC>
   template <class T>
 std::uint64_t VtkWriterInterface<GV,DC>
-  ::writeAppended (std::ofstream& out, std::vector<T> const& values) const
+  ::writeValuesAppended (std::ofstream& out, std::vector<T> const& values) const
 {
   assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
   pos_type begin_pos = out.tellp();

dune/vtk/writers/vtkimagedatawriter.hh

@@ -29,8 +29,10 @@ namespace Dune
 
   public:
     /// Constructor, stores the gridView
-    VtkImageDataWriter (GridView const& gridView)
-      : Super(gridView)
+    VtkImageDataWriter (GridView const& gridView,
+                        Vtk::FormatTypes format = Vtk::BINARY,
+                        Vtk::DataTypes datatype = Vtk::FLOAT32)
+      : Super(gridView, format, datatype)
     {}
 
   private:
@@ -47,6 +49,8 @@ namespace Dune
       return "vti";
     }
 
+    virtual void writeGridAppended (std::ofstream& /*out*/, std::vector<std::uint64_t>& /*blocks*/) const override {}
+
   private:
     using Super::dataCollector_;
     using Super::format_;

dune/vtk/writers/vtkimagedatawriter.impl.hh

@@ -67,37 +67,8 @@ void VtkImageDataWriter<GV,DC>
   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 (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) );
-    }
-    out << "</AppendedData>\n";
-  }
-
+  this->writeAppended(out, offsets);
   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]);
-    }
-  }
 }
 
 

dune/vtk/writers/vtkrectilineargridwriter.hh

@@ -29,8 +29,10 @@ namespace Dune
 
   public:
     /// Constructor, stores the gridView
-    VtkRectilinearGridWriter (GridView const& gridView)
-      : Super(gridView)
+    VtkRectilinearGridWriter (GridView const& gridView,
+                              Vtk::FormatTypes format = Vtk::BINARY,
+                              Vtk::DataTypes datatype = Vtk::FLOAT32)
+      : Super(gridView, format, datatype)
     {}
 
   private:
@@ -42,7 +44,8 @@ namespace Dune
     /// 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;
+    void writeCoordinates (std::ofstream& out, std::vector<pos_type>& offsets,
+                           Std::optional<std::size_t> timestep = {}) const;
 
     template <class T>
     std::array<std::uint64_t, 3> writeCoordinatesAppended (std::ofstream& out) const;
@@ -52,6 +55,8 @@ namespace Dune
       return "vtr";
     }
 
+    virtual void writeGridAppended (std::ofstream& out, std::vector<std::uint64_t>& blocks) const override;