diff --git a/dune/vtk/utility/uid.hh b/dune/vtk/utility/uid.hh new file mode 100644 index 0000000000000000000000000000000000000000..22c2bb7e7bf4eb72c4aa2062cea325b39c428825 --- /dev/null +++ b/dune/vtk/utility/uid.hh @@ -0,0 +1,22 @@ +#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 diff --git a/dune/vtk/vtktimeserieswriter.hh b/dune/vtk/vtktimeserieswriter.hh new file mode 100644 index 0000000000000000000000000000000000000000..bb4d646942fa7dbe1af1908ee6d97970934c3daf --- /dev/null +++ b/dune/vtk/vtktimeserieswriter.hh @@ -0,0 +1,93 @@ +#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" diff --git a/dune/vtk/vtktimeserieswriter.impl.hh b/dune/vtk/vtktimeserieswriter.impl.hh new file mode 100644 index 0000000000000000000000000000000000000000..37b35527c9b6d9df66ed7171667470925e8b9a18 --- /dev/null +++ b/dune/vtk/vtktimeserieswriter.impl.hh @@ -0,0 +1,101 @@ +#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 diff --git a/dune/vtk/vtkwriterinterface.hh b/dune/vtk/vtkwriterinterface.hh index dab897bea478d54ee8b6344da729bb992c32cb23..d5ae604fe58169ef75cbb228b2810e880da69bf6 100644 --- a/dune/vtk/vtkwriterinterface.hh +++ b/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_; diff --git a/dune/vtk/vtkwriterinterface.impl.hh b/dune/vtk/vtkwriterinterface.impl.hh index 26883e356b091c2b8d8900381969a66b4284cffc..945dc4bbcd2ce56da815ed354ec4ccfa74ac5dcd 100644 --- a/dune/vtk/vtkwriterinterface.impl.hh +++ b/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(); diff --git a/dune/vtk/writers/vtkimagedatawriter.hh b/dune/vtk/writers/vtkimagedatawriter.hh index d3b2af6fec89b8e9ce56efd560073f3917aab817..cc6f04e772562b96bc617e8d16a9b3da6c61ede4 100644 --- a/dune/vtk/writers/vtkimagedatawriter.hh +++ b/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_; diff --git a/dune/vtk/writers/vtkimagedatawriter.impl.hh b/dune/vtk/writers/vtkimagedatawriter.impl.hh index 235b4ef89186a2e3f772c9a4b3d4c7246111d2c6..b976b428d85686ca7ff1c69f40b429290d594ac8 100644 --- a/dune/vtk/writers/vtkimagedatawriter.impl.hh +++ b/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]); - } - } } diff --git a/dune/vtk/writers/vtkrectilineargridwriter.hh b/dune/vtk/writers/vtkrectilineargridwriter.hh index b40fdabe21dfd13be82f992a8c31c42e0fd27408..050e899bd9fc9dc3c2c9d518bce388c13b7586be 100644 --- a/dune/vtk/writers/vtkrectilineargridwriter.hh +++ b/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; + private: using Super::dataCollector_; using Super::format_; diff --git a/dune/vtk/writers/vtkrectilineargridwriter.impl.hh b/dune/vtk/writers/vtkrectilineargridwriter.impl.hh index 23b8f09a179abc27519e709217a687fcf8dab722..48a650ee3646ff6e7de7a23840fd560e4dfd180c 100644 --- a/dune/vtk/writers/vtkrectilineargridwriter.impl.hh +++ b/dune/vtk/writers/vtkrectilineargridwriter.impl.hh @@ -47,6 +47,11 @@ void VtkRectilinearGridWriter<GV,DC> out << "<Piece Extent=\"" << join(extent.begin(), extent.end()) << "\">\n"; }); + // Write point coordinates for x, y, and z ordinate + out << "<Coordinates>\n"; + writeCoordinates(out, offsets); + out << "</Coordinates>\n"; + // Write data associated with grid points out << "<PointData" << this->getNames(pointData_) << ">\n"; for (auto const& v : pointData_) @@ -59,53 +64,11 @@ void VtkRectilinearGridWriter<GV,DC> 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"; - } - + 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]); - } - } } @@ -131,6 +94,13 @@ void VtkRectilinearGridWriter<GV,DC> << " WholeExtent=\"" << join(wholeExtent.begin(), wholeExtent.end()) << "\"" << ">\n"; + // Write point coordinates for x, y, and z ordinate + out << "<PCoordinates>\n"; + out << "<PDataArray Name=\"x\" type=\"" << to_string(datatype_) << "\" />\n"; + out << "<PDataArray Name=\"y\" type=\"" << to_string(datatype_) << "\" />\n"; + out << "<PDataArray Name=\"z\" type=\"" << to_string(datatype_) << "\" />\n"; + out << "</PCoordinates>\n"; + // Write data associated with grid points out << "<PPointData" << this->getNames(pointData_) << ">\n"; for (auto const& v : pointData_) { @@ -153,13 +123,6 @@ void VtkRectilinearGridWriter<GV,DC> } out << "</PCellData>\n"; - // Write point coordinates for x, y, and z ordinate - out << "<PCoordinates>\n"; - out << "<PDataArray Name=\"x\" type=\"" << to_string(datatype_) << "\" />\n"; - out << "<PDataArray Name=\"y\" type=\"" << to_string(datatype_) << "\" />\n"; - out << "<PDataArray Name=\"z\" type=\"" << to_string(datatype_) << "\" />\n"; - out << "</PCoordinates>\n"; - // Write piece file references dataCollector_.writePieces([&out,pfilename,ext=this->fileExtension()](int p, auto const& extent, bool write_extent) { @@ -177,13 +140,17 @@ void VtkRectilinearGridWriter<GV,DC> template <class GV, class DC> void VtkRectilinearGridWriter<GV,DC> - ::writeCoordinates (std::ofstream& out, std::vector<pos_type>& offsets) const + ::writeCoordinates (std::ofstream& out, std::vector<pos_type>& offsets, + Std::optional<std::size_t> timestep) const { std::string names = "xyz"; if (format_ == Vtk::ASCII) { auto coordinates = dataCollector_.template coordinates<double>(); for (std::size_t d = 0; d < 3; ++d) { - out << "<DataArray type=\"" << to_string(datatype_) << "\" Name=\"" << names[d] << "\" format=\"ascii\">\n"; + out << "<DataArray type=\"" << to_string(datatype_) << "\" Name=\"" << names[d] << "\" format=\"ascii\""; + if (timestep) + out << " TimeStep=\"" << *timestep << "\""; + out << ">\n"; std::size_t i = 0; for (auto const& c : coordinates[d]) out << c << (++i % 6 != 0 ? ' ' : '\n'); @@ -193,6 +160,8 @@ void VtkRectilinearGridWriter<GV,DC> else { // Vtk::APPENDED format for (std::size_t j = 0; j < 3; ++j) { out << "<DataArray type=\"" << to_string(datatype_) << "\" Name=\"" << names[j] << "\" format=\"appended\""; + if (timestep) + out << " TimeStep=\"" << *timestep << "\""; out << " offset="; offsets.push_back(out.tellp()); out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' '); @@ -203,20 +172,23 @@ void VtkRectilinearGridWriter<GV,DC> template <class GV, class DC> - template <class T> -std::array<std::uint64_t,3> VtkRectilinearGridWriter<GV,DC> - ::writeCoordinatesAppended (std::ofstream& out) const +void VtkRectilinearGridWriter<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"); - auto coordinates = dataCollector_.template coordinates<T>(); - - // write conncetivity, offsets, and types - std::uint64_t bs0 = this->writeAppended(out, coordinates[0]); - std::uint64_t bs1 = this->writeAppended(out, coordinates[1]); - std::uint64_t bs2 = this->writeAppended(out, coordinates[2]); - - return {bs0, bs1, bs2}; + // 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])); + } } } // end namespace Dune diff --git a/dune/vtk/writers/vtkstructuredgridwriter.hh b/dune/vtk/writers/vtkstructuredgridwriter.hh index beb9d564e3347ab0bbf22b10379fae47ff530c3e..f6ebc6a5c2b70afe66f02ffbcffef8e1ebf2850c 100644 --- a/dune/vtk/writers/vtkstructuredgridwriter.hh +++ b/dune/vtk/writers/vtkstructuredgridwriter.hh @@ -29,8 +29,10 @@ namespace Dune public: /// Constructor, stores the gridView - VtkStructuredGridWriter (GridView const& gridView) - : Super(gridView) + VtkStructuredGridWriter (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 "vts"; } + virtual void writeGridAppended (std::ofstream& out, std::vector<std::uint64_t>& blocks) const override; + private: using Super::dataCollector_; using Super::format_; diff --git a/dune/vtk/writers/vtkstructuredgridwriter.impl.hh b/dune/vtk/writers/vtkstructuredgridwriter.impl.hh index 61c69cb8f9936f9f784eb7b6e39283cf970a8b8b..1dcd7241622c62d7b0649dd828a31fa17bfe1c13 100644 --- a/dune/vtk/writers/vtkstructuredgridwriter.impl.hh +++ b/dune/vtk/writers/vtkstructuredgridwriter.impl.hh @@ -45,6 +45,11 @@ void VtkStructuredGridWriter<GV,DC> out << "<Piece Extent=\"" << join(extent.begin(), extent.end()) << "\">\n"; }); + // Write point coordinates + out << "<Points>\n"; + this->writePoints(out, offsets); + out << "</Points>\n"; + // Write data associated with grid points out << "<PointData" << this->getNames(pointData_) << ">\n"; for (auto const& v : pointData_) @@ -57,50 +62,11 @@ void VtkStructuredGridWriter<GV,DC> 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 (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) - blocks.push_back( this->template writePointsAppended<float>(out) ); - else - blocks.push_back( this->template writePointsAppended<double>(out) ); - 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]); - } - } } @@ -126,6 +92,14 @@ void VtkStructuredGridWriter<GV,DC> << " WholeExtent=\"" << join(wholeExtent.begin(), wholeExtent.end()) << "\"" << ">\n"; + // Write points + out << "<PPoints>\n"; + out << "<PDataArray" + << " type=\"" << to_string(datatype_) << "\"" + << " NumberOfComponents=\"3\"" + << " />\n"; + out << "</PPoints>\n"; + // Write data associated with grid points out << "<PPointData" << this->getNames(pointData_) << ">\n"; for (auto const& v : pointData_) { @@ -148,14 +122,6 @@ void VtkStructuredGridWriter<GV,DC> } out << "</PCellData>\n"; - // Write points - out << "<PPoints>\n"; - out << "<PDataArray" - << " type=\"" << to_string(datatype_) << "\"" - << " NumberOfComponents=\"3\"" - << " />\n"; - out << "</PPoints>\n"; - // Write piece file references dataCollector_.writePieces([&out,pfilename,ext=this->fileExtension()](int p, auto const& extent, bool write_extent) { @@ -170,4 +136,21 @@ void VtkStructuredGridWriter<GV,DC> out << "</VTKFile>"; } + +template <class GV, class DC> +void VtkStructuredGridWriter<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"); + + // write points + if (datatype_ == Vtk::FLOAT32) { + auto points = dataCollector_.template points<float>(); + blocks.push_back(this->writeValuesAppended(out, points)); + } else { + auto points = dataCollector_.template points<double>(); + blocks.push_back(this->writeValuesAppended(out, points)); + } +} + } // end namespace Dune diff --git a/dune/vtk/writers/vtkunstructuredgridwriter.hh b/dune/vtk/writers/vtkunstructuredgridwriter.hh index 221cb3931b624837b0d7635a5cf285976c49b948..2429214253eda0c0986ad2b935426e39a697a504 100644 --- a/dune/vtk/writers/vtkunstructuredgridwriter.hh +++ b/dune/vtk/writers/vtkunstructuredgridwriter.hh @@ -22,6 +22,8 @@ namespace Dune class VtkUnstructuredGridWriter : public VtkWriterInterface<GridView, DataCollector> { + template <class> friend class VtkTimeseriesWriter; + static constexpr int dimension = GridView::dimension; using Super = VtkWriterInterface<GridView, DataCollector>; @@ -29,8 +31,10 @@ namespace Dune public: /// Constructor, stores the gridView - VtkUnstructuredGridWriter (GridView const& gridView) - : Super(gridView) + VtkUnstructuredGridWriter (GridView const& gridView, + Vtk::FormatTypes format = Vtk::BINARY, + Vtk::DataTypes datatype = Vtk::FLOAT32) + : Super(gridView, format, datatype) {} private: @@ -42,20 +46,37 @@ namespace Dune /// for [i] in [0,...,size). virtual void writeParallelFile (std::string const& pfilename, int size) const override; + /// Write a series of timesteps in one file + /** + * \param filename The name of the output file + * \param filenameMesh The name of a file where the mesh is stored. Must exist. + * \param timesteps A vector of pairs (timestep, filename) where the filename indicates + * a file where the data of the timestep is stored. + * \param blocks A list of block sizes of the binary data stored in the files. + * Order: (points, cells, pointdata[0], celldata[0], pointdata[1], celldata[1],...) + **/ + void writeTimeseriesSerialFile (std::string const& filename, + std::string const& filenameMesh, + std::vector<std::pair<double, std::string>> const& timesteps, + std::vector<std::uint64_t> const& blocks) const; + + /// Write parallel VTK file for series of timesteps + void writeTimeseriesParallelFile (std::string const& pfilename, int size, + std::vector<std::pair<double, std::string>> const& timesteps) const; + virtual std::string fileExtension () const override { return "vtu"; } + virtual void writeGridAppended (std::ofstream& out, std::vector<std::uint64_t>& blocks) const override; + // 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 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; + std::vector<pos_type>& offsets, + Std::optional<std::size_t> timestep = {}) const; private: using Super::dataCollector_; diff --git a/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh b/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh index 42769b600cc169ad1936a5fe9ac2d4c2b5c19bae..9bc44820d8741cb2729ab202fce3a3a6146d911b 100644 --- a/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh +++ b/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh @@ -45,6 +45,16 @@ void VtkUnstructuredGridWriter<GV,DC> << " NumberOfCells=\"" << dataCollector_.numCells() << "\"" << ">\n"; + // Write point coordinates + out << "<Points>\n"; + this->writePoints(out, offsets); + out << "</Points>\n"; + + // Write element connectivity, types and offsets + out << "<Cells>\n"; + writeCells(out, offsets); + out << "</Cells>\n"; + // Write data associated with grid points out << "<PointData" << this->getNames(pointData_) << ">\n"; for (auto const& v : pointData_) @@ -57,58 +67,11 @@ void VtkUnstructuredGridWriter<GV,DC> this->writeData(out, offsets, v, Super::CELL_DATA); out << "</CellData>\n"; - // Write point coordinates - out << "<Points>\n"; - this->writePoints(out, offsets); - out << "</Points>\n"; - - // Write element connectivity, types and offsets - out << "<Cells>\n"; - writeCells(out, offsets); - out << "</Cells>\n"; - out << "</Piece>\n"; out << "</UnstructuredGrid>\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) - blocks.push_back( this->template writePointsAppended<float>(out) ); - else - blocks.push_back( this->template writePointsAppended<double>(out) ); - - auto bs = writeCellsAppended(out); - blocks.insert(blocks.end(), bs.begin(), bs.end()); - 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]); - } - } } @@ -130,6 +93,14 @@ void VtkUnstructuredGridWriter<GV,DC> out << "<PUnstructuredGrid GhostLevel=\"0\">\n"; + // Write points + out << "<PPoints>\n"; + out << "<PDataArray" + << " type=\"" << to_string(datatype_) << "\"" + << " NumberOfComponents=\"3\"" + << " />\n"; + out << "</PPoints>\n"; + // Write data associated with grid points out << "<PPointData" << this->getNames(pointData_) << ">\n"; for (auto const& v : pointData_) { @@ -152,6 +123,159 @@ void VtkUnstructuredGridWriter<GV,DC> } out << "</PCellData>\n"; + // Write piece file references + for (int p = 0; p < size; ++p) { + std::string piece_source = pfilename + "_p" + std::to_string(p) + "." + this->fileExtension(); + out << "<Piece Source=\"" << piece_source << "\" />\n"; + } + + out << "</PUnstructuredGrid>\n"; + out << "</VTKFile>"; +} + + +template <class GV, class DC> +void VtkUnstructuredGridWriter<GV,DC> + ::writeTimeseriesSerialFile (std::string const& filename, + std::string const& filenameMesh, + std::vector<std::pair<double, std::string>> const& timesteps, + std::vector<std::uint64_t> const& blocks) const +{ + std::ofstream out(filename, std::ios_base::ate | std::ios::binary); + assert(out.is_open()); + + assert(is_a(format_, Vtk::APPENDED)); + 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<std::vector<pos_type>> offsets(timesteps.size()); // pos => offset + out << "<VTKFile" + << " type=\"UnstructuredGrid\"" + << " version=\"1.0\"" + << " byte_order=\"" << this->getEndian() << "\"" + << " header_type=\"UInt64\"" + << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\"" : "") + << ">\n"; + + out << "<UnstructuredGrid" + << " TimeValues=\""; + { + std::size_t i = 0; + for (auto const& timestep : timesteps) + out << timestep.first << (++i % 6 != 0 ? ' ' : '\n'); + } + out << "\">\n"; + + out << "<Piece" + << " NumberOfPoints=\"" << dataCollector_.numPoints() << "\"" + << " NumberOfCells=\"" << dataCollector_.numCells() << "\"" + << ">\n"; + + // Write point coordinates + out << "<Points>\n"; + for (std::size_t i = 0; i < timesteps.size(); ++i) + this->writePoints(out, offsets[i], i); + out << "</Points>\n"; + + // Write element connectivity, types and offsets + out << "<Cells>\n"; + for (std::size_t i = 0; i < timesteps.size(); ++i) + writeCells(out, offsets[i], i); + out << "</Cells>\n"; + + const std::size_t shift = offsets[0].size(); // number of blocks to write the grid + + // Write data associated with grid points + out << "<PointData" << this->getNames(pointData_) << ">\n"; + for (std::size_t i = 0; i < timesteps.size(); ++i) { + for (auto const& v : pointData_) + this->writeData(out, offsets[i], v, Super::POINT_DATA, i); + } + out << "</PointData>\n"; + + // Write data associated with grid cells + out << "<CellData" << this->getNames(cellData_) << ">\n"; + for (std::size_t i = 0; i < timesteps.size(); ++i) { + for (auto const& v : cellData_) + this->writeData(out, offsets[i], v, Super::CELL_DATA, i); + } + out << "</CellData>\n"; + + out << "</Piece>\n"; + out << "</UnstructuredGrid>\n"; + + out << "<AppendedData encoding=\"raw\">\n_"; + pos_type appended_pos = out.tellp(); + + { // write grid (points, cells) + std::ifstream file_mesh(filenameMesh, std::ios_base::in | std::ios_base::binary); + out << file_mesh.rdbuf(); + assert( std::uint64_t(out.tellp()) == std::accumulate(blocks.begin(), std::next(blocks.begin(),shift), std::uint64_t(appended_pos)) ); + } + + // write point-data and cell-data + for (auto const& timestep : timesteps) { + std::ifstream file(timestep.second, std::ios_base::in | std::ios_base::binary); + out << file.rdbuf(); + } + out << "</AppendedData>\n"; + + out << "</VTKFile>"; + + // write correct offsets in file. + pos_type offset = 0; + for (std::size_t i = 0; i < timesteps.size(); ++i) { + offset = 0; + auto const& off = offsets[i]; + + // write mesh data offsets + for (std::size_t j = 0; j < shift; ++j) { + out.seekp(off[j]); + out << '"' << offset << '"'; + offset += pos_type(blocks[j]); + } + } + + std::size_t j = shift; + for (std::size_t i = 0; i < timesteps.size(); ++i) { + auto const& off = offsets[i]; + + for (std::size_t k = shift; k < off.size(); ++k) { + out.seekp(off[k]); + out << '"' << offset << '"'; + offset += pos_type(blocks[j++]); + } + } +} + + +template <class GV, class DC> +void VtkUnstructuredGridWriter<GV,DC> + ::writeTimeseriesParallelFile (std::string const& pfilename, int size, std::vector<std::pair<double, std::string>> const& timesteps) const +{ + std::string filename = pfilename + ".pvtu"; + std::ofstream out(filename, std::ios_base::ate | std::ios::binary); + assert(out.is_open()); + + out << "<VTKFile" + << " type=\"PUnstructuredGrid\"" + << " version=\"1.0\"" + << " byte_order=\"" << this->getEndian() << "\"" + << " header_type=\"UInt64\"" + << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\"" : "") + << ">\n"; + + out << "<PUnstructuredGrid GhostLevel=\"0\"" + << " TimeValues=\""; + { + std::size_t i = 0; + for (auto const& timestep : timesteps) + out << timestep.first << (++i % 6 != 0 ? ' ' : '\n'); + } + out << "\">\n"; + // Write points out << "<PPoints>\n"; out << "<PDataArray" @@ -160,6 +284,34 @@ void VtkUnstructuredGridWriter<GV,DC> << " />\n"; out << "</PPoints>\n"; + // Write data associated with grid points + out << "<PPointData" << this->getNames(pointData_) << ">\n"; + for (std::size_t i = 0; i < timesteps.size(); ++i) { + for (auto const& v : pointData_) { + out << "<PDataArray" + << " Name=\"" << v.name() << "\"" + << " type=\"" << to_string(v.type()) << "\"" + << " NumberOfComponents=\"" << v.ncomps() << "\"" + << " TimeStep=\"" << i << "\"" + << " />\n"; + } + } + out << "</PPointData>\n"; + + // Write data associated with grid cells + out << "<PCellData" << this->getNames(cellData_) << ">\n"; + for (std::size_t i = 0; i < timesteps.size(); ++i) { + for (auto const& v : cellData_) { + out << "<PDataArray" + << " Name=\"" << v.name() << "\"" + << " type=\"" << to_string(v.type()) << "\"" + << " NumberOfComponents=\"" << v.ncomps() << "\"" + << " TimeStep=\"" << i << "\"" + << " />\n"; + } + } + out << "</PCellData>\n"; + // Write piece file references for (int p = 0; p < size; ++p) { std::string piece_source = pfilename + "_p" + std::to_string(p) + "." + this->fileExtension(); @@ -173,23 +325,33 @@ void VtkUnstructuredGridWriter<GV,DC> template <class GV, class DC> void VtkUnstructuredGridWriter<GV,DC> - ::writeCells (std::ofstream& out, std::vector<pos_type>& offsets) const + ::writeCells (std::ofstream& out, std::vector<pos_type>& offsets, + Std::optional<std::size_t> timestep) const { if (format_ == Vtk::ASCII) { auto cells = dataCollector_.cells(); - out << "<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">\n"; + out << "<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\""; + if (timestep) + out << " TimeStep=\"" << *timestep << "\""; + out << ">\n"; std::size_t i = 0; for (auto const& c : cells.connectivity) out << c << (++i % 6 != 0 ? ' ' : '\n'); out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n"; - out << "<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">\n"; + out << "<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\""; + if (timestep) + out << " TimeStep=\"" << *timestep << "\""; + out << ">\n"; i = 0; for (auto const& o : cells.offsets) out << o << (++i % 6 != 0 ? ' ' : '\n'); out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n"; - out << "<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n"; + out << "<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\""; + if (timestep) + out << " TimeStep=\"" << *timestep << "\""; + out << ">\n"; i = 0; for (auto const& t : cells.types) out << int(t) << (++i % 6 != 0 ? ' ' : '\n'); @@ -197,18 +359,24 @@ void VtkUnstructuredGridWriter<GV,DC> } else { // Vtk::APPENDED format out << "<DataArray type=\"Int64\" Name=\"connectivity\" format=\"appended\""; + if (timestep) + out << " TimeStep=\"" << *timestep << "\""; out << " offset="; offsets.push_back(out.tellp()); out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' '); out << "/>\n"; out << "<DataArray type=\"Int64\" Name=\"offsets\" format=\"appended\""; + if (timestep) + out << " TimeStep=\"" << *timestep << "\""; out << " offset="; offsets.push_back(out.tellp()); out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' '); out << "/>\n"; out << "<DataArray type=\"UInt8\" Name=\"types\" format=\"appended\""; + if (timestep) + out << " TimeStep=\"" << *timestep << "\""; out << " offset="; offsets.push_back(out.tellp()); out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' '); @@ -218,20 +386,25 @@ void VtkUnstructuredGridWriter<GV,DC> template <class GV, class DC> -std::array<std::uint64_t,3> VtkUnstructuredGridWriter<GV,DC> - ::writeCellsAppended (std::ofstream& out) const +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"); - auto cells = dataCollector_.cells(); + // write points + if (datatype_ == Vtk::FLOAT32) { + auto points = dataCollector_.template points<float>(); + blocks.push_back(this->writeValuesAppended(out, points)); + } else { + auto points = dataCollector_.template points<double>(); + blocks.push_back(this->writeValuesAppended(out, points)); + } // write conncetivity, offsets, and types - std::uint64_t bs0 = this->writeAppended(out, cells.connectivity); - std::uint64_t bs1 = this->writeAppended(out, cells.offsets); - std::uint64_t bs2 = this->writeAppended(out, cells.types); - - return {bs0, bs1, bs2}; + 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)); } - } // end namespace Dune diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 46958c1bfe2f4151d9f5308adc860ff1fc2c170d..3b6012ab3b883a66622327ecc4ca254bc37ebe6f 100644 --- a/src/CMakeLists.txt +++ b/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("timeserieswriter" timeserieswriter.cc) + target_link_dune_default_libraries("timeserieswriter") + target_link_libraries("timeserieswriter" dunevtk) endif () if (dune-polygongrid_FOUND) diff --git a/src/benchmark.cc b/src/benchmark.cc index e1d731adcda8575d8534d1058f32d9927f86e91e..2439c286d192f354e84f0842bf85f5177bb46134 100644 --- a/src/benchmark.cc +++ b/src/benchmark.cc @@ -61,11 +61,10 @@ template <class GridView> void writer_new (GridView const& gridView) { Timer t; - VtkUnstructuredGridWriter<GridView> vtkWriter(gridView); for (auto const& test_case : test_cases_new) { t.reset(); - vtkWriter.write("writer_new_" + std::get<0>(test_case) + ".vtu", - std::get<1>(test_case), std::get<2>(test_case)); + VtkUnstructuredGridWriter<GridView> vtkWriter(gridView, std::get<1>(test_case), std::get<2>(test_case)); + vtkWriter.write("writer_new_" + std::get<0>(test_case) + ".vtu"); std::cout << " time (writer_new_" + std::get<0>(test_case) + ") = " << t.elapsed() << "\n"; } } diff --git a/src/datacollector.cc b/src/datacollector.cc index 5864a4f013561604623d31c180613863cde51807..482f45d24cd0ca2776061cdf43ff0cf48ca6814a 100644 --- a/src/datacollector.cc +++ b/src/datacollector.cc @@ -32,14 +32,13 @@ using namespace Dune::Functions; template <class DataCollector, class GridView, class Fct1, class Fct2> void write_dc (std::string prefix, GridView const& gridView, Fct1 const& fct1, Fct2 const& fct2) { - VtkUnstructuredGridWriter<GridView, DataCollector> vtkWriter(gridView); + VtkUnstructuredGridWriter<GridView, DataCollector> vtkWriter(gridView, Vtk::ASCII, Vtk::FLOAT32); vtkWriter.addPointData(fct1, "p1"); vtkWriter.addCellData(fct1, "p0"); vtkWriter.addPointData(fct2, "q1"); vtkWriter.addCellData(fct2, "q0"); - vtkWriter.write(prefix + "_" + std::to_string(GridView::dimension) + "d_ascii.vtu", - Vtk::ASCII, Vtk::FLOAT32); + vtkWriter.write(prefix + "_" + std::to_string(GridView::dimension) + "d_ascii.vtu"); } template <class GridView> diff --git a/src/geometrygrid.cc b/src/geometrygrid.cc index 8b56ffd7e49b9f02677c9f514624ce623e147e48..1fe36b37d9cb650d6d858203072fb92d99abc1b5 100644 --- a/src/geometrygrid.cc +++ b/src/geometrygrid.cc @@ -53,10 +53,10 @@ void write (std::string prefix, GridView const& gridView) auto p1Analytic = makeAnalyticGridViewFunction([&c](auto const& x) { return c.dot(x); }, gridView); - VtkWriter<GridView> vtkWriter(gridView); + VtkWriter<GridView> vtkWriter(gridView, Vtk::ASCII, Vtk::FLOAT32); vtkWriter.addPointData(p1Analytic, "q1"); vtkWriter.addCellData(p1Analytic, "q0"); - vtkWriter.write(prefix + "_ascii.vtu", Vtk::ASCII, Vtk::FLOAT32); + vtkWriter.write(prefix + "_ascii.vtu"); } int main (int argc, char** argv) diff --git a/src/legacyvtkwriter.cc b/src/legacyvtkwriter.cc index b001fce11fc11b0a02dec28d0833698393952fd7..c28688d0e09c6acfeb278d162d1e580711d2398e 100644 --- a/src/legacyvtkwriter.cc +++ b/src/legacyvtkwriter.cc @@ -38,7 +38,7 @@ int main(int argc, char** argv) std::shared_ptr<VTKFunction<GridView> const> p1FctWrapped(new P1Function(gridView, p1function, "p1")); using Writer = VtkUnstructuredGridWriter<GridView>; - Writer vtkWriter(gridView); + Writer vtkWriter(gridView, Vtk::ASCII); vtkWriter.addPointData(p1FctWrapped); - vtkWriter.write("test_ascii_float32.vtu", Vtk::ASCII); + vtkWriter.write("test_ascii_float32.vtu"); } diff --git a/src/polygongrid.cc b/src/polygongrid.cc index fd24fd3ed3ef824650b3c801e7e0a1ca76511660..626c339c5f600461ebed7bc8a4d9a890c02a0ac3 100644 --- a/src/polygongrid.cc +++ b/src/polygongrid.cc @@ -57,7 +57,7 @@ int main(int argc, char** argv) GridView gridView = gridPtr->leafGridView(); using Writer = VtkUnstructuredGridWriter<GridView>; - Writer vtkWriter(gridView); + Writer vtkWriter(gridView, Vtk::ASCII); auto p1Analytic = makeAnalyticGridViewFunction([](auto const& x) { return std::sin(10*x[0])*std::cos(10*x[1]); }, gridView); @@ -65,10 +65,5 @@ int main(int argc, char** argv) vtkWriter.addPointData(p1Analytic, "p1"); vtkWriter.addCellData(p1Analytic, "p0"); - vtkWriter.write("poly_ascii_float32.vtu", Vtk::ASCII); - vtkWriter.write("poly_binary_float32.vtu", Vtk::BINARY); - vtkWriter.write("poly_compressed_float32.vtu", Vtk::COMPRESSED); - vtkWriter.write("poly_ascii_float64.vtu", Vtk::ASCII, Vtk::FLOAT64); - vtkWriter.write("poly_binary_float64.vtu", Vtk::BINARY, Vtk::FLOAT64); - vtkWriter.write("poly_compressed_float64.vtu", Vtk::COMPRESSED, Vtk::FLOAT64); + vtkWriter.write("poly_ascii_float32.vtu"); } diff --git a/src/structuredgridwriter.cc b/src/structuredgridwriter.cc index 349244a50ef50b4fdbf8fddeaf02f771a808ef55..e28067ab3ab51ec5baba377bf6b0cf5b37a69b67 100644 --- a/src/structuredgridwriter.cc +++ b/src/structuredgridwriter.cc @@ -45,30 +45,30 @@ void write(std::string prefix, GridView const& gridView) { using Writer = VtkImageDataWriter<GridView>; - Writer vtkWriter(gridView); + Writer vtkWriter(gridView, Vtk::ASCII, Vtk::FLOAT32); vtkWriter.addPointData(fct2, "analytic"); - vtkWriter.write(prefix + "id_ascii_float32.vti", Vtk::ASCII, Vtk::FLOAT32); + vtkWriter.write(prefix + "id_ascii_float32.vti"); } { using Writer = VtkRectilinearGridWriter<GridView>; - Writer vtkWriter(gridView); + Writer vtkWriter(gridView, Vtk::ASCII, Vtk::FLOAT32); vtkWriter.addPointData(fct2, "analytic"); - vtkWriter.write(prefix + "rg_ascii_float32.vtr", Vtk::ASCII, Vtk::FLOAT32); + vtkWriter.write(prefix + "rg_ascii_float32.vtr"); } { using Writer = VtkStructuredGridWriter<GridView>; - Writer vtkWriter(gridView); + Writer vtkWriter(gridView, Vtk::ASCII, Vtk::FLOAT32); vtkWriter.addPointData(fct2, "analytic"); - vtkWriter.write(prefix + "sg_ascii_float32.vts", Vtk::ASCII, Vtk::FLOAT32); + vtkWriter.write(prefix + "sg_ascii_float32.vts"); } { using Writer = VtkUnstructuredGridWriter<GridView>; - Writer vtkWriter(gridView); + Writer vtkWriter(gridView, Vtk::ASCII, Vtk::FLOAT32); vtkWriter.addPointData(fct2, "analytic"); - vtkWriter.write(prefix + "ug_ascii_float32.vts", Vtk::ASCII, Vtk::FLOAT32); + vtkWriter.write(prefix + "ug_ascii_float32.vts"); } } diff --git a/src/test/mixed_element_test.cc b/src/test/mixed_element_test.cc index 6b15247c0d03c25a4efde4c113164b9cb2b01d65..6dc61ef49783918d2c199ceb0ca767e38620e428 100644 --- a/src/test/mixed_element_test.cc +++ b/src/test/mixed_element_test.cc @@ -69,10 +69,9 @@ static TestCases test_cases = { template <class GridView> void writer_test (GridView const& gridView) { - VtkUnstructuredGridWriter<GridView> vtkWriter(gridView); for (auto const& test_case : test_cases) { - vtkWriter.write("/tmp/reader_writer_test_" + std::get<0>(test_case) + ".vtu", - std::get<1>(test_case), std::get<2>(test_case)); + VtkUnstructuredGridWriter<GridView> vtkWriter(gridView, std::get<1>(test_case), std::get<2>(test_case)); + vtkWriter.write("/tmp/reader_writer_test_" + std::get<0>(test_case) + ".vtu"); } } @@ -81,9 +80,9 @@ void reader_test (Test& test) { for (auto const& test_case : test_cases) { auto grid = VtkReader<Grid>::read("/tmp/reader_writer_test_" + std::get<0>(test_case) + ".vtu"); - VtkUnstructuredGridWriter<typename Grid::LeafGridView> vtkWriter(grid->leafGridView()); - vtkWriter.write("/tmp/reader_writer_test_" + std::get<0>(test_case) + "_2.vtu", + VtkUnstructuredGridWriter<typename Grid::LeafGridView> vtkWriter(grid->leafGridView(), std::get<1>(test_case), std::get<2>(test_case)); + vtkWriter.write("/tmp/reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"); test.check(compare_files("/tmp/reader_writer_test_" + std::get<0>(test_case) + ".vtu", "/tmp/reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"), std::get<0>(test_case)); } diff --git a/src/test/reader_writer_test.cc b/src/test/reader_writer_test.cc index d34cbd86b1574d2b0051108c5bc39a6c5a5897dc..da3e1a1585d31b52176d03f95cf31acda5e8eebb 100644 --- a/src/test/reader_writer_test.cc +++ b/src/test/reader_writer_test.cc @@ -69,10 +69,9 @@ static TestCases test_cases = { template <class GridView> void writer_test (GridView const& gridView) { - VtkUnstructuredGridWriter<GridView> vtkWriter(gridView); for (auto const& test_case : test_cases) { - vtkWriter.write("/tmp/reader_writer_test_" + std::get<0>(test_case) + ".vtu", - std::get<1>(test_case), std::get<2>(test_case)); + VtkUnstructuredGridWriter<GridView> vtkWriter(gridView, std::get<1>(test_case), std::get<2>(test_case)); + vtkWriter.write("/tmp/reader_writer_test_" + std::get<0>(test_case) + ".vtu"); } } @@ -81,9 +80,9 @@ void reader_test (Test& test) { for (auto const& test_case : test_cases) { auto grid = VtkReader<Grid>::read("/tmp/reader_writer_test_" + std::get<0>(test_case) + ".vtu"); - VtkUnstructuredGridWriter<typename Grid::LeafGridView> vtkWriter(grid->leafGridView()); - vtkWriter.write("/tmp/reader_writer_test_" + std::get<0>(test_case) + "_2.vtu", + VtkUnstructuredGridWriter<typename Grid::LeafGridView> vtkWriter(grid->leafGridView(), std::get<1>(test_case), std::get<2>(test_case)); + vtkWriter.write("/tmp/reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"); test.check(compare_files("/tmp/reader_writer_test_" + std::get<0>(test_case) + ".vtu", "/tmp/reader_writer_test_" + std::get<0>(test_case) + "_2.vtu")); } diff --git a/src/timeserieswriter.cc b/src/timeserieswriter.cc new file mode 100644 index 0000000000000000000000000000000000000000..e999b222e2923bb777405f2965b235c5df93ceba --- /dev/null +++ b/src/timeserieswriter.cc @@ -0,0 +1,55 @@ +// -*- 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/gridfunctions/analyticgridviewfunction.hh> +#include <dune/grid/yaspgrid.hh> + +#include <dune/vtk/vtktimeserieswriter.hh> +#include <dune/vtk/writers/vtkunstructuredgridwriter.hh> + +using namespace Dune; +using namespace Dune::Functions; + +template <class GridView> +void write (std::string prefix, GridView const& gridView) +{ + FieldVector<double,GridView::dimension> c{11.0, 7.0, 3.0}; + auto p1Analytic = makeAnalyticGridViewFunction([&c](auto const& x) -> float { return c.dot(x); }, gridView); + + using Writer = VtkUnstructuredGridWriter<GridView>; + VtkTimeseriesWriter<Writer> seriesWriter(gridView, Vtk::BINARY, Vtk::FLOAT32); + seriesWriter.addPointData(p1Analytic, "q1"); + seriesWriter.addCellData(p1Analytic, "q0"); + std::string filename = prefix + "_" + std::to_string(GridView::dimension) + "d_binary32.vtu"; + for (double t = 0.0; t < 5; t += 0.5) { + seriesWriter.writeTimestep(t, filename); + } + seriesWriter.write(filename); + + Writer vtkWriter(gridView, Vtk::BINARY, Vtk::FLOAT32); + vtkWriter.addPointData(p1Analytic, "q1"); + vtkWriter.addCellData(p1Analytic, "q0"); + vtkWriter.write(filename); +} + +int main (int argc, char** argv) +{ + Dune::MPIHelper::instance(argc, argv); + + using GridType = YaspGrid<3>; + FieldVector<double,3> upperRight; upperRight = 1.0; + auto numElements = filledArray<3,int>(8); + GridType grid(upperRight, numElements); + write("yasp", grid.leafGridView()); +} \ No newline at end of file diff --git a/src/vtkreader.cc b/src/vtkreader.cc index b02fceeb59096b2808125b79ab597ab238752751..d82e96ba9848bf87c0bfaa65851b18dd08074e6f 100644 --- a/src/vtkreader.cc +++ b/src/vtkreader.cc @@ -37,42 +37,46 @@ int main(int argc, char** argv) GridView gridView = grid.leafGridView(); - VtkUnstructuredGridWriter<GridView> vtkWriter(gridView); - vtkWriter.write("test_ascii_float32.vtu", Vtk::ASCII); - vtkWriter.write("test_binary_float32.vtu", Vtk::BINARY); - vtkWriter.write("test_compressed_float64.vtu", Vtk::COMPRESSED, Vtk::FLOAT64); + VtkUnstructuredGridWriter<GridView> vtkWriter1(gridView, Vtk::ASCII); + vtkWriter1.write("test_ascii_float32.vtu"); + + VtkUnstructuredGridWriter<GridView> vtkWriter2(gridView, Vtk::BINARY); + vtkWriter2.write("test_binary_float32.vtu"); + + VtkUnstructuredGridWriter<GridView> vtkWriter3(gridView, Vtk::COMPRESSED, Vtk::FLOAT64); + vtkWriter3.write("test_compressed_float64.vtu"); } { auto gridPtr = VtkReader<GridType>::read("test_ascii_float32.vtu"); auto& grid = *gridPtr; - VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView()); - vtkWriter.write("test_ascii_float32_2.vtu", Vtk::ASCII); + VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView(), Vtk::ASCII); + vtkWriter.write("test_ascii_float32_2.vtu"); } { auto gridPtr = VtkReader<GridType>::read("test_binary_float32.vtu"); auto& grid = *gridPtr; - VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView()); - vtkWriter.write("test_ascii_float32_3.vtu", Vtk::ASCII); + VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView(), Vtk::ASCII); + vtkWriter.write("test_ascii_float32_3.vtu"); } { auto gridPtr = VtkReader<GridType>::read("test_compressed_float64.vtu"); auto& grid = *gridPtr; - VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView()); - vtkWriter.write("test_ascii_float64_3.vtu", Vtk::ASCII, Vtk::FLOAT64); + VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView(), Vtk::ASCII, Vtk::FLOAT64); + vtkWriter.write("test_ascii_float64_3.vtu"); } { auto gridPtr = VtkReader<GridType,ConnectedGridCreator>::read("test_ascii_float32.vtu"); auto& grid = *gridPtr; - VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView()); - vtkWriter.write("test_ascii_float32_4.vtu", Vtk::ASCII); + VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView(), Vtk::ASCII); + vtkWriter.write("test_ascii_float32_4.vtu"); } if (filesystem::exists("paraview_3d.vtu")) { @@ -81,8 +85,8 @@ int main(int argc, char** argv) auto gridPtr = VtkReader<GridType3d>::read("paraview_3d.vtu"); auto& grid = *gridPtr; - VtkUnstructuredGridWriter<GridView3d> vtkWriter(grid.leafGridView()); - vtkWriter.write("paraview_3d_ascii.vtu", Vtk::ASCII, Vtk::FLOAT64); + VtkUnstructuredGridWriter<GridView3d> vtkWriter(grid.leafGridView(), Vtk::ASCII, Vtk::FLOAT64); + vtkWriter.write("paraview_3d_ascii.vtu"); } if (filesystem::exists("paraview_2d.vtu")) { @@ -92,7 +96,7 @@ int main(int argc, char** argv) auto gridPtr = VtkReader<GridType2d>::read("paraview_2d.vtu"); auto& grid = *gridPtr; - VtkUnstructuredGridWriter<GridView2d> vtkWriter(grid.leafGridView()); - vtkWriter.write("paraview_2d_ascii.vtu", Vtk::ASCII, Vtk::FLOAT64); + VtkUnstructuredGridWriter<GridView2d> vtkWriter(grid.leafGridView(), Vtk::ASCII, Vtk::FLOAT64); + vtkWriter.write("paraview_2d_ascii.vtu"); } } diff --git a/src/vtkwriter.cc b/src/vtkwriter.cc index b15b5ff0fdb5d48c2a597710ac8d7d6d9d7e9e6d..d6819bcb8b39437258b825a5713af09d08f19663 100644 --- a/src/vtkwriter.cc +++ b/src/vtkwriter.cc @@ -54,14 +54,13 @@ void write (std::string prefix, GridView const& gridView) // write analytic function auto p1Analytic = makeAnalyticGridViewFunction([&c](auto const& x) { return c.dot(x); }, gridView); - VtkWriter<GridView> vtkWriter(gridView); - vtkWriter.addPointData(p1Interpol, "p1"); - vtkWriter.addCellData(p1Interpol, "p0"); - vtkWriter.addPointData(p1Analytic, "q1"); - vtkWriter.addCellData(p1Analytic, "q0"); for (auto const& test_case : test_cases) { - vtkWriter.write(prefix + "_" + std::to_string(GridView::dimension) + "d_" + std::get<0>(test_case) + ".vtu", - std::get<1>(test_case), std::get<2>(test_case)); + VtkWriter<GridView> vtkWriter(gridView, std::get<1>(test_case), std::get<2>(test_case)); + vtkWriter.addPointData(p1Interpol, "p1"); + vtkWriter.addCellData(p1Interpol, "p0"); + vtkWriter.addPointData(p1Analytic, "q1"); + vtkWriter.addCellData(p1Analytic, "q0"); + vtkWriter.write(prefix + "_" + std::to_string(GridView::dimension) + "d_" + std::get<0>(test_case) + ".vtu"); } }