Implement writer with flexible header type

Add method to set compressor type

README.md

@@ -43,11 +43,11 @@ See also the `src/` directory for more examples.
 
 ## Comparison with Dune::VTKWriter
 In Dune-Grid there is a VTK writer available, that is a bit different from the
-proposed one. A comparions:
+proposed one. A comparison:
 
 | **Property**       | **Dune-Grid** | **Dune-Vtk** |
 | ------------------ | :-----------: | :----------: |
-| VTK version        | 0.1           | 1.0          |
+| VTK version        | 0.1           | 0.1/1.0      |
 | UnstructuredGrid   | **x**         | **x**        |
 | PolyData           | (1d)          | -            |
 | StructuredGrid     | -             | **x**        |

dune/vtk/function.hh

@@ -58,7 +58,7 @@ namespace Dune
         : localFct_(fct)
         , name_(fct->name())
         , ncomps_(fct->ncomps())
-        , type_(type ? *type : Vtk::FLOAT64)
+        , type_(type ? *type : Vtk::FLOAT32)
       {}
 
       /// Construct VtkFunction from dune-functions GridFunction with Signature

dune/vtk/types.cc

@@ -39,6 +39,18 @@ std::string to_string (DataTypes type)
   }
 }
 
+std::string to_string (CompressorTypes type)
+{
+  switch (type) {
+    case ZLIB: return "vtkZLibDataCompressor";
+    case LZ4:  return "vtkLZ4DataCompressor";
+    case LZMA: return "vtkLZMADataCompressor";
+    default:
+      DUNE_THROW(RangeError, "CompressorTypes not found.");
+      std::abort();
+  }
+}
+
 GeometryType to_geometry (std::uint8_t cell)
 {
   switch (cell) {

dune/vtk/types.hh

@@ -39,6 +39,7 @@ namespace Dune
       LZ4,
       LZMA
     };
+    std::string to_string (CompressorTypes);
 
     enum CellParametrization {
       LINEAR,

dune/vtk/vtkreader.impl.hh

@@ -304,30 +304,6 @@ void VtkReader<Grid,Creator,Field>::readParallelFileFromStream (std::ifstream& i
 }
 
 
-template <class Grid, class Creator, class Field>
-  template <class Caller>
-void VtkReader<Grid,Creator,Field>::readAppendedCaller(Vtk::DataTypes floatType, Vtk::DataTypes headerType, Caller caller)
-{
-  switch (headerType) {
-    case Vtk::UINT32:
-      switch (floatType) {
-        case Vtk::FLOAT32: caller(MetaType<float>{},  MetaType<std::uint32_t>{}); break;
-        case Vtk::FLOAT64: caller(MetaType<double>{}, MetaType<std::uint32_t>{}); break;
-        default: VTK_ASSERT_MSG(false, "Unsupported FloatingPoint type"); break;
-      }
-      break;
-    case Vtk::UINT64:
-      switch (floatType) {
-        case Vtk::FLOAT32: caller(MetaType<float>{},  MetaType<std::uint64_t>{}); break;
-        case Vtk::FLOAT64: caller(MetaType<double>{}, MetaType<std::uint64_t>{}); break;
-        default: VTK_ASSERT_MSG(false, "Unsupported FloatingPoint type"); break;
-      }
-      break;
-    default: VTK_ASSERT_MSG(false, "Unsupported Header Integer type"); break;
-  }
-}
-
-
 // @{ implementation detail
 /**
  * Read ASCII data from `input` stream into vector `values`

dune/vtk/vtkwriterinterface.hh

@@ -50,29 +50,34 @@ namespace Dune
      *
      * \param gridView  Implementation of Dune::GridView
      * \param format    Format of the VTK file, either Vtk::BINARY, Vtk::ASCII, or Vtk::COMPRESSED
-     * \param datatype  Output datatype used for coordinates and other global floating point values
+     * \param datatype  Data type of a single component of the point coordinates [Vtk::FLOAT32]
+     * \param headertype  Integer type used in binary data headers [Vtk::UINT32]
      **/
     VtkWriterInterface (GridView const& gridView,
-                        Vtk::FormatTypes format = Vtk::BINARY,
-                        Vtk::DataTypes datatype = Vtk::FLOAT32)
-      : VtkWriterInterface(std::make_shared<DataCollector>(gridView), format, datatype)
+                        Vtk::FormatTypes format = Vtk::FormatTypes::BINARY,
+                        Vtk::DataTypes datatype = Vtk::DataTypes::FLOAT32,
+                        Vtk::DataTypes headertype = Vtk::DataTypes::UINT32)
+      : VtkWriterInterface(std::make_shared<DataCollector>(gridView), format, datatype, headertype)
     {}
 
     /// \brief Constructor, wraps the passed DataCollector in a non-destroying shared_ptr
     VtkWriterInterface (DataCollector& dataCollector,
-                        Vtk::FormatTypes format = Vtk::BINARY,
-                        Vtk::DataTypes datatype = Vtk::FLOAT32)
-      : VtkWriterInterface(stackobject_to_shared_ptr(dataCollector), format, datatype)
+                        Vtk::FormatTypes format = Vtk::FormatTypes::BINARY,
+                        Vtk::DataTypes datatype = Vtk::DataTypes::FLOAT32,
+                        Vtk::DataTypes headertype = Vtk::DataTypes::UINT32)
+      : VtkWriterInterface(stackobject_to_shared_ptr(dataCollector), format, datatype, headertype)
     {}
 
     /// \brief Constructor, stores the passed DataCollector
     VtkWriterInterface (std::shared_ptr<DataCollector> dataCollector,
-                        Vtk::FormatTypes format = Vtk::BINARY,
-                        Vtk::DataTypes datatype = Vtk::FLOAT32)
+                        Vtk::FormatTypes format = Vtk::FormatTypes::BINARY,
+                        Vtk::DataTypes datatype = Vtk::DataTypes::FLOAT32,
+                        Vtk::DataTypes headertype = Vtk::DataTypes::UINT32)
       : dataCollector_(std::move(dataCollector))
     {
       setFormat(format);
       setDatatype(datatype);
+      setHeadertype(headertype);
     }
 
 
@@ -120,20 +125,43 @@ namespace Dune
     {
       format_ = format;
 
-#if !HAVE_VTK_ZLIB
-      if (format_ == Vtk::COMPRESSED) {
+      if (format_ == Vtk::FormatTypes::COMPRESSED) {
+#if HAVE_VTK_ZLIB
+        compressor_ = Vtk::CompressorTypes::ZLIB;
+#else
         std::cout << "Dune is compiled without compression. Falling back to BINARY VTK output!\n";
-        format_ = Vtk::BINARY;
-      }
+        format_ = Vtk::FormatTypes::BINARY;
 #endif
+      } else {
+        compressor_ = Vtk::CompressorTypes::NONE;
+      }
+
     }
 
-    // Sets the global datatype used for coordinates and other global float values
+    /// Sets the global datatype used for coordinates and other global float values
     void setDatatype (Vtk::DataTypes datatype)
     {
       datatype_ = datatype;
     }
 
+    /// Sets the integer type used in binary data headers
+    void setHeadertype (Vtk::DataTypes datatype)
+    {
+      headertype_ = datatype;
+    }
+
+    /// Sets the compressor type used in binary data headers, Additionally a compression
+    /// level can be passed with level = -1 means: default compression level. Level must be in [0-9]
+    void setCompressor (Vtk::CompressorTypes compressor, int level = -1)
+    {
+      compressor_ = compressor;
+      compression_level = level;
+      VTK_ASSERT(level >= -1 && level <= 9);
+
+      if (compressor_ != Vtk::CompressorTypes::NONE)
+        format_ = Vtk::FormatTypes::COMPRESSED;
+    }
+
   private:
     /// Write a serial VTK file in Unstructured format
     virtual void writeSerialFile (std::ofstream& out) const = 0;
@@ -174,8 +202,8 @@ namespace Dune
 
     // Write the `values` in blocks (possibly compressed) to the output
     // stream `out`. Return the written block size.
-    template <class T>
-    std::uint64_t writeValuesAppended (std::ofstream& out, std::vector<T> const& values) const;
+    template <class HeaderType, class FloatType>
+    std::uint64_t writeValuesAppended (std::ofstream& out, std::vector<FloatType> const& values) const;
 
     // Write the `values` in a space and newline separated list of ascii representations.
     // The precision is controlled by the datatype and numerical_limits::digits10.
@@ -224,6 +252,8 @@ namespace Dune
 
     Vtk::FormatTypes format_;
     Vtk::DataTypes datatype_;
+    Vtk::DataTypes headertype_;
+    Vtk::CompressorTypes compressor_ = Vtk::CompressorTypes::NONE;
 
     // attached data
     std::vector<VtkFunction> pointData_;

dune/vtk/vtkwriterinterface.impl.hh

@@ -17,6 +17,7 @@
 #include <dune/geometry/type.hh>
 
 #include <dune/vtk/utility/enum.hh>
+#include <dune/vtk/utility/errors.hh>
 #include <dune/vtk/utility/filesystem.hh>
 #include <dune/vtk/utility/string.hh>
 
@@ -125,19 +126,27 @@ void VtkWriterInterface<GV,DC>
   }
 }
 
+
 template <class GV, class DC>
 void VtkWriterInterface<GV,DC>
   ::writeDataAppended (std::ofstream& out, std::vector<std::uint64_t>& blocks) const
 {
   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)));
+    Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(v.type(), headertype_,
+    [&](auto f, auto h) {
+      using F = typename decltype(f)::type;
+      using H = typename decltype(h)::type;
+      blocks.push_back(this->template writeValuesAppended<H>(out, dataCollector_->template pointData<F>(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)));
+    Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(v.type(), headertype_,
+    [&](auto f, auto h) {
+      using F = typename decltype(f)::type;
+      using H = typename decltype(h)::type;
+      blocks.push_back(this->template writeValuesAppended<H>(out, dataCollector_->template cellData<F>(v)));
+    });
   }
 }
 
@@ -146,7 +155,7 @@ template <class GV, class DC>
 void VtkWriterInterface<GV,DC>
   ::writeAppended (std::ofstream& out, std::vector<pos_type> const& offsets) const
 {
-  if (is_a(format_, Vtk::APPENDED)) {
+  if (is_a(format_, Vtk::FormatTypes::APPENDED)) {
     out << "<AppendedData encoding=\"raw\">\n_";
     std::vector<std::uint64_t> blocks;
     writeGridAppended(out, blocks);
@@ -197,11 +206,11 @@ void VtkWriterInterface<GV,DC>
   out << "<VTKFile"
       << " type=\"" << type << "\""
       << " version=\"1.0\""
-      << " header_type=\"UInt64\"";
+      << " header_type=\"" << to_string(headertype_) << "\"";
   if (format_ != Vtk::ASCII)
     out << " byte_order=\"" << getEndian() << "\"";
-  if (format_ == Vtk::COMPRESSED)
-    out << " compressor=\"vtkZLibDataCompressor\"";
+  if (compressor_ != Vtk::NONE)
+    out << " compressor=\"" << to_string(compressor_) << "\"";
   out << ">\n";
 }
 
@@ -218,10 +227,8 @@ std::uint64_t writeValuesToBuffer (std::size_t max_num_values, unsigned char* bu
   return bs;
 }
 
-
-template <class OStream>
-std::uint64_t writeCompressed (unsigned char const* buffer, unsigned char* buffer_out,
-                               std::uint64_t bs, std::uint64_t cbs, int level, OStream& outb)
+inline std::uint64_t compressBuffer_zlib (unsigned char const* buffer, unsigned char* buffer_out,
+                                          std::uint64_t bs, std::uint64_t cbs, int level)
 {
 #if HAVE_VTK_ZLIB
   uLongf uncompressed_space = uLongf(bs);
@@ -233,8 +240,6 @@ std::uint64_t writeCompressed (unsigned char const* buffer, unsigned char* buffe
   if (compress2(out, &compressed_space, in, uncompressed_space, level) != Z_OK) {
     std::cerr << "Zlib error while compressing data.\n";
     std::abort();
-  } else {
-    outb.write((char*)out, compressed_space);
   }
 
   return compressed_space;
@@ -245,55 +250,97 @@ std::uint64_t writeCompressed (unsigned char const* buffer, unsigned char* buffe
 #endif
 }
 
+inline std::uint64_t compressBuffer_lz4 (unsigned char const* /* buffer */, unsigned char* /* buffer_out */,
+                                         std::uint64_t /* bs */, std::uint64_t /* cbs */, int /* level */)
+{
+#if HAVE_VTK_LZ4
+  std::cerr << "LZ4 Compression not yet implemented" << std::endl;
+  std::abort();
+  return 0;
+#else
+  std::cerr << "LZ4 Compression not supported. Provide the LZ4 package to CMake." << std::endl;
+  std::abort();
+  return 0;
+#endif
+}
+
+inline std::uint64_t compressBuffer_lzma (unsigned char const* /* buffer */, unsigned char* /* buffer_out */,
+                                          std::uint64_t /* bs */, std::uint64_t /* cbs */, int /* level */)
+{
+#if HAVE_VTK_LZMA
+  std::cerr << "LZMA Compression not yet implemented" << std::endl;
+  std::abort();
+  return 0;
+#else
+  std::cerr << "LZMA Compression not supported. Provide the LZMA package to CMake." << std::endl;
+  std::abort();
+  return 0;
+#endif
+}
+
 } // end namespace Impl
 
 template <class GV, class DC>
-  template <class T>
+  template <class HeaderType, class FloatType>
 std::uint64_t VtkWriterInterface<GV,DC>
-  ::writeValuesAppended (std::ofstream& out, std::vector<T> const& values) const
+  ::writeValuesAppended (std::ofstream& out, std::vector<FloatType> const& values) const
 {
   assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
   pos_type begin_pos = out.tellp();
 
-  std::uint64_t size = values.size() * sizeof(T);
+  HeaderType size = values.size() * sizeof(FloatType);
 
-  std::uint64_t num_full_blocks = size / block_size;
-  std::uint64_t last_block_size = size % block_size;
-  std::uint64_t num_blocks = num_full_blocks + (last_block_size > 0 ? 1 : 0);
+  HeaderType num_full_blocks = size / block_size;
+  HeaderType last_block_size = size % block_size;
+  HeaderType num_blocks = num_full_blocks + (last_block_size > 0 ? 1 : 0);
 
   // write block-size(s)
-  std::uint64_t zero = 0;
-  if (format_ == Vtk::COMPRESSED) {
-    out.write((char*)&num_blocks, sizeof(std::uint64_t));
-    out.write((char*)&block_size, sizeof(std::uint64_t));
-    out.write((char*)&last_block_size, sizeof(std::uint64_t));
-    for (std::uint64_t i = 0; i < num_blocks; ++i)
-      out.write((char*)&zero, sizeof(std::uint64_t));
+  HeaderType zero = 0;
+  if (compressor_ != Vtk::NONE) {
+    out.write((char*)&num_blocks, sizeof(HeaderType));
+    out.write((char*)&block_size, sizeof(HeaderType));
+    out.write((char*)&last_block_size, sizeof(HeaderType));
+    for (HeaderType i = 0; i < num_blocks; ++i)
+      out.write((char*)&zero, sizeof(HeaderType));
   } else {
-    out.write((char*)&size, sizeof(std::uint64_t));
+    out.write((char*)&size, sizeof(HeaderType));
   }
 
-  std::uint64_t compressed_block_size = block_size + (block_size + 999)/1000 + 12;
+  HeaderType compressed_block_size = block_size + (block_size + 999)/1000 + 12;
   std::vector<unsigned char> buffer(block_size);
   std::vector<unsigned char> buffer_out;
-  std::size_t num_values = block_size / sizeof(T);
+  if (compressor_ != Vtk::NONE)
+    buffer_out.resize(std::size_t(compressed_block_size));
 
-  std::vector<std::uint64_t> cbs(std::size_t(num_blocks), 0); // compressed block sizes
+  std::size_t num_values = block_size / sizeof(FloatType);
+
+  std::vector<HeaderType> cbs(std::size_t(num_blocks), 0); // compressed block sizes
   for (std::size_t i = 0; i < std::size_t(num_blocks); ++i) {
-    std::uint64_t bs = Impl::writeValuesToBuffer<T>(num_values, buffer.data(), values, i*num_values);
+    HeaderType bs = Impl::writeValuesToBuffer<FloatType>(num_values, buffer.data(), values, i*num_values);
 
-    if (format_ == Vtk::COMPRESSED) {
-      buffer_out.resize(std::size_t(compressed_block_size));
-      cbs[i] = Impl::writeCompressed(buffer.data(), buffer_out.data(), bs,
-                                     compressed_block_size, compression_level, out);
-    } else
+    switch (compressor_) {
+      case Vtk::NONE:
         out.write((char*)buffer.data(), bs);
+        break;
+      case Vtk::ZLIB:
+        cbs[i] = Impl::compressBuffer_zlib(buffer.data(), buffer_out.data(), bs, compressed_block_size, compression_level);
+        out.write((char*)buffer_out.data(), cbs[i]);
+        break;
+      case Vtk::LZ4:
+        cbs[i] = Impl::compressBuffer_zlib(buffer.data(), buffer_out.data(), bs, compressed_block_size, compression_level);
+        out.write((char*)buffer_out.data(), cbs[i]);
+        break;
+      case Vtk::LZMA:
+        cbs[i] = Impl::compressBuffer_zlib(buffer.data(), buffer_out.data(), bs, compressed_block_size, compression_level);
+        out.write((char*)buffer_out.data(), cbs[i]);
+        break;
+    }
   }
 
   pos_type end_pos = out.tellp();
-  if (format_ == Vtk::COMPRESSED) {
-    out.seekp(begin_pos + std::streamoff(3*sizeof(std::uint64_t)));
-    out.write((char*)cbs.data(), std::streamsize(num_blocks*sizeof(std::uint64_t)));
+  if (compressor_ != Vtk::NONE) {
+    out.seekp(begin_pos + std::streamoff(3*sizeof(HeaderType)));
+    out.write((char*)cbs.data(), std::streamsize(num_blocks*sizeof(HeaderType)));
     out.seekp(end_pos);
   }
 

dune/vtk/writers/vtkimagedatawriter.hh

@@ -51,6 +51,7 @@ namespace Dune
     using Super::dataCollector_;
     using Super::format_;
     using Super::datatype_;
+    using Super::headertype_;
 
     // attached data
     using Super::pointData_;

dune/vtk/writers/vtkrectilineargridwriter.hh

@@ -57,6 +57,7 @@ namespace Dune
     using Super::dataCollector_;
     using Super::format_;
     using Super::datatype_;
+    using Super::headertype_;
 
     // attached data
     using Super::pointData_;

dune/vtk/writers/vtkrectilineargridwriter.impl.hh

@@ -151,17 +151,15 @@ void VtkRectilinearGridWriter<GV,DC>
   assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
 
   // write coordinates along axis
-  if (datatype_ == Vtk::FLOAT32) {
-    auto coordinates = dataCollector_->template coordinates<float>();
-    blocks.push_back(this->writeValuesAppended(out, coordinates[0]));
-    blocks.push_back(this->writeValuesAppended(out, coordinates[1]));
-    blocks.push_back(this->writeValuesAppended(out, coordinates[2]));
-  } else {
-    auto coordinates = dataCollector_->template coordinates<double>();
-    blocks.push_back(this->writeValuesAppended(out, coordinates[0]));
-    blocks.push_back(this->writeValuesAppended(out, coordinates[1]));
-    blocks.push_back(this->writeValuesAppended(out, coordinates[2]));
-  }
+  Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(datatype_, headertype_,
+  [&](auto f, auto h) {
+    using F = typename decltype(f)::type;
+    using H = typename decltype(h)::type;
+    auto coordinates = dataCollector_->template coordinates<F>();
+    blocks.push_back(this->template writeValuesAppended<H>(out, coordinates[0]));
+    blocks.push_back(this->template writeValuesAppended<H>(out, coordinates[1]));
+    blocks.push_back(this->template writeValuesAppended<H>(out, coordinates[2]));
+  });
 }
 
 } // end namespace Dune

dune/vtk/writers/vtkstructuredgridwriter.hh

@@ -51,6 +51,7 @@ namespace Dune
     using Super::dataCollector_;
     using Super::format_;
     using Super::datatype_;
+    using Super::headertype_;
 
     // attached data
     using Super::pointData_;

dune/vtk/writers/vtkstructuredgridwriter.impl.hh

@@ -119,9 +119,12 @@ void VtkStructuredGridWriter<GV,DC>
   assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
 
   // write points
-  blocks.push_back( datatype_ == Vtk::FLOAT32
-    ? this->writeValuesAppended(out, dataCollector_->template points<float>())
-    : this->writeValuesAppended(out, dataCollector_->template points<double>()) );
+  Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(datatype_, headertype_,
+  [&](auto f, auto h) {
+    using F = typename decltype(f)::type;
+    using H = typename decltype(h)::type;
+    blocks.push_back(this->template writeValuesAppended<H>(out, dataCollector_->template points<F>()));
+  });
 }
 
 } // end namespace Dune

dune/vtk/writers/vtkunstructuredgridwriter.hh

@@ -83,6 +83,7 @@ namespace Dune
     using Super::dataCollector_;
     using Super::format_;
     using Super::datatype_;
+    using Super::headertype_;
 
     // attached data
     using Super::pointData_;

dune/vtk/writers/vtkunstructuredgridwriter.impl.hh

@@ -333,6 +333,7 @@ void VtkUnstructuredGridWriter<GV,DC>
   }
 }
 
+
 template <class GV, class DC>
 void VtkUnstructuredGridWriter<GV,DC>
   ::writePointIds (std::ofstream& out,
@@ -362,26 +363,32 @@ void VtkUnstructuredGridWriter<GV,DC>
   }
 }
 
+
 template <class GV, class DC>
 void VtkUnstructuredGridWriter<GV,DC>
   ::writeGridAppended (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");
 
+  Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(datatype_, headertype_,
+  [&](auto f, auto h) {
+    using F = typename decltype(f)::type;
+    using H = typename decltype(h)::type;
+
     // write points
-  blocks.push_back( datatype_ == Vtk::FLOAT32
-    ? this->writeValuesAppended(out, dataCollector_->template points<float>())
-    : this->writeValuesAppended(out, dataCollector_->template points<double>()) );
+    blocks.push_back(this->template writeValuesAppended<H>(out, dataCollector_->template points<F>()));
 
-  // write conncetivity, offsets, and types
+    // write connectivity, offsets, and types
     auto cells = dataCollector_->cells();
-  blocks.push_back(this->writeValuesAppended(out, cells.connectivity));
-  blocks.push_back(this->writeValuesAppended(out, cells.offsets));
-  blocks.push_back(this->writeValuesAppended(out, cells.types));
+    blocks.push_back(this->template writeValuesAppended<H>(out, cells.connectivity));
+    blocks.push_back(this->template writeValuesAppended<H>(out, cells.offsets));
+    blocks.push_back(this->template writeValuesAppended<H>(out, cells.types));
 
+    // optionally, write global point IDs
     auto ids = dataCollector_->pointIds();
     if (!ids.empty())
-    blocks.push_back(this->writeValuesAppended(out, ids));
+      blocks.push_back(this->template writeValuesAppended<H>(out, ids));
+  });
 }
 
 } // end namespace Dune