Added a ParaView Data writer, i.e. pvd writer

dune/vtk/pvdwriter.hh

+#pragma once
+
+#include <iosfwd>
+#include <string>
+#include <vector>
+#include <tuple>
+
+#include <dune/vtk/vtktypes.hh>
+#include <dune/vtk/vtkwriterinterface.hh>
+
+namespace Dune
+{
+  /// File-Writer for ParaView .pvd files
+  template <class VtkWriter>
+  class PvdWriter
+  {
+    using Self = PvdWriter;
+
+    static_assert(IsVtkWriter<VtkWriter>::value, "Writer must implement the VtkWriterInterface.");
+
+  public:
+    /// Constructor, creates a VtkWriter with constructor arguments forwarded
+    template <class... Args, disableCopyMove<Self,Args...> = 0>
+    explicit PvdWriter (Args&&... args)
+      : vtkWriter_{std::forward<Args>(args)...}
+    {}
+
+    /// Write the attached data to the file
+    void write (double time, std::string const& fn)
+    {
+      write(time, fn, Vtk::BINARY);
+    }
+
+    /// Write the attached data to the file with \ref Vtk::FormatTypes and \ref Vtk::DataTypes
+    void write (double time,
+                std::string const& fn,
+                Vtk::FormatTypes format,
+                Vtk::DataTypes datatype = Vtk::FLOAT32);
+
+    /// Attach point data to the writer, \see VtkFunction for possible arguments
+    template <class Function, class... Args>
+    PvdWriter& 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>
+    PvdWriter& addCellData (Function const& fct, Args&&... args)
+    {
+      vtkWriter_.addCellData(fct, std::forward<Args>(args)...);
+      return *this;
+    }
+
+  protected:
+    /// Write a series of vtk files in a .pvd ParaView Data file
+    void writeFile (double time, std::string const& filename) const;
+
+  protected:
+    VtkWriter vtkWriter_;
+    std::vector<std::pair<double, std::string>> timeSeries_;
+    Vtk::FormatTypes format_;
+    Vtk::DataTypes datatype_;
+  };
+
+} // end namespace Dune
+
+#include "pvdwriter.impl.hh"

dune/vtk/pvdwriter.impl.hh

+#pragma once
+
+#include <iomanip>
+
+#include <dune/vtk/utility/filesystem.hh>
+#include <dune/vtk/utility/string.hh>
+
+namespace Dune {
+
+template <class W>
+void PvdWriter<W>
+  ::write (double time, std::string const& fn, Vtk::FormatTypes format, Vtk::DataTypes datatype)
+{
+  format_ = format;
+  datatype_ = datatype;
+
+#ifndef HAVE_ZLIB
+  if (format_ == Vtk::COMPRESSED)
+    format_ = Vtk::BINARY;
+#endif
+
+  auto p = filesystem::path(fn);
+  auto name = p.stem();
+  p.remove_filename();
+  p /= name.string();
+
+  std::string ext = "." + vtkWriter_.getFileExtension();
+  std::string filename = p.string() + "_t" + std::to_string(timeSeries_.size());
+
+  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)
+    ext = ".p" + vtkWriter_.getFileExtension();
+#endif
+
+  timeSeries_.emplace_back(time, filename + ext);
+  vtkWriter_.write(filename + ext, format_, datatype_);
+
+  if (rank == 0)
+    writeFile(time, p.string() + ".pvd");
+}
+
+
+template <class W>
+void PvdWriter<W>
+  ::writeFile (double time, std::string const& filename) const
+{
+  std::ofstream out(filename, std::ios_base::ate | std::ios::binary);
+  assert(out.is_open());
+
+  if (datatype_ == Vtk::FLOAT32)
+    out << std::setprecision(std::numeric_limits<float>::digits10+2);
+  else
+    out << std::setprecision(std::numeric_limits<double>::digits10+2);
+
+  out << "<?xml version=\"1.0\"?>\n";
+  out << "<VTKFile"
+      << " type=\"Collection\""
+      << " version=\"0.1\""
+      << (format_ != Vtk::ASCII ? " byte_order=\"" << vtkWriter_.getEndian() << "\"" : "")
+      << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\"" : "")
+      << ">\n";
+
+  out << "<Collection>\n";
+
+  // Write all timesteps
+  for (auto const& timestep : timeSeries_) {
+    out << "<DataSet"
+        << " timestep=\"" << timestep.first << "\""
+        << " part=\"0\""
+        << " file=\"" << timestep.second << "\""
+        << " />\n";
+  }
+
+  out << "</Collection>\n";
+  out << "</VTKFile>";
+}
+
+} // end namespace Dune

dune/vtk/vtkwriterinterface.hh

@@ -5,19 +5,23 @@
 #include <string>
 #include <vector>
 
-#include <dune/common/std/optional.hh>
-
 #include <dune/vtk/filewriter.hh>
 #include <dune/vtk/vtkfunction.hh>
 #include <dune/vtk/vtktypes.hh>
 
 namespace Dune
 {
-  /// File-Writer for Vtk .vtu files
+  // forward declaration
+  template <class VtkWriter>
+  class PvdWriter;
+
+  /// Interface for file writers for the Vtk XML file formats
   template <class GridView, class DataCollector>
   class VtkWriterInterface
       : public FileWriter
   {
+    template <class> friend class PvdWriter;
+
   protected:
     static constexpr int dimension = GridView::dimension;
 
@@ -62,7 +66,7 @@ namespace Dune
       return *this;
     }
 
-  protected:
+  private:
     /// Write a serial VTK file in Unstructured format
     virtual void writeSerialFile (std::string const& filename) const = 0;
 
@@ -74,6 +78,7 @@ namespace Dune
     /// Return the file extension of the serial file (not including the dot)
     virtual std::string fileExtension () const = 0;
 
+  protected:
     // 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`.
@@ -114,10 +119,15 @@ namespace Dune
       return (reinterpret_cast<char*>(&i)[1] == 1 ? "BigEndian" : "LittleEndian");
     }
 
+    // provide accessor to \ref fileExtension virtual method
+    std::string getFileExtension () const
+    {
+      return fileExtension();
+    }
+
   protected:
     mutable DataCollector dataCollector_;
 
-    std::string filename_;
     Vtk::FormatTypes format_;
     Vtk::DataTypes datatype_;
 
@@ -129,6 +139,17 @@ namespace Dune
     int compression_level = -1; // in [0,9], -1 ... use default value
   };
 
+
+  template <class Writer>
+  struct IsVtkWriter
+  {
+    template <class GV, class DC>
+    static std::uint16_t test(VtkWriterInterface<GV,DC> const&);
+    static std::uint8_t  test(...); // fall-back overload
+
+    static constexpr bool value = sizeof(decltype(test(std::declval<Writer>()))) > sizeof(std::uint8_t);
+  };
+
 } // end namespace Dune
 
 #include "vtkwriterinterface.impl.hh"

src/CMakeLists.txt

@@ -27,6 +27,10 @@ if (dune-functions_FOUND)
   add_executable("geometrygrid" geometrygrid.cc)
   target_link_dune_default_libraries("geometrygrid")
   target_link_libraries("geometrygrid" dunevtk)
+
+  add_executable("pvdwriter" pvdwriter.cc)
+  target_link_dune_default_libraries("pvdwriter")
+  target_link_libraries("pvdwriter" dunevtk)
 endif ()
 
 if (dune-polygongrid_FOUND)

src/pvdwriter.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <iostream>
+#include <vector>
+
+#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
+#include <dune/common/exceptions.hh> // We use exceptions
+#include <dune/common/filledarray.hh>
+
+#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dune/vtk/writers/vtkunstructuredgridwriter.hh>
+#include <dune/vtk/pvdwriter.hh>
+
+using namespace Dune;
+using namespace Dune::Functions;
+
+template <class GridView>
+void write (std::string prefix, GridView const& gridView)
+{
+  using namespace BasisFactory;
+  auto basis = makeBasis(gridView, lagrange<1>());
+
+  FieldVector<double,GridView::dimension> c;
+  if (GridView::dimension > 0) c[0] = 11.0;
+  if (GridView::dimension > 1) c[1] = 7.0;
+  if (GridView::dimension > 2) c[2] = 3.0;
+
+  // write analytic function
+  auto p1Analytic = makeAnalyticGridViewFunction([&c](auto const& x) { return c.dot(x); }, gridView);
+
+  using Writer = VtkUnstructuredGridWriter<GridView>;
+  PvdWriter<Writer> pvdWriter(gridView);
+
+  pvdWriter.addPointData(p1Analytic, "p1");
+  pvdWriter.addCellData(p1Analytic, "p0");
+  for (double t = 0.0; t < 10.0; t += 1.0) {
+    pvdWriter.write(t, prefix + "_" + std::to_string(GridView::dimension) + "d_ascii.vtu",
+      Vtk::ASCII, Vtk::FLOAT32);
+  }
+}
+
+template <int I>
+using int_ = std::integral_constant<int,I>;
+
+int main (int argc, char** argv)
+{
+  Dune::MPIHelper::instance(argc, argv);
+
+  // Test PvdWriter for YaspGrid
+  using GridType = YaspGrid<2>;
+  FieldVector<double,2> upperRight; upperRight = 1.0;
+  auto numElements = filledArray<2,int>(4);
+  GridType grid(upperRight, numElements);
+  write("yasp", grid.leafGridView());
+}
\ No newline at end of file