From 8998ec6e193a4a90884e96a0bd8dc364029240ff Mon Sep 17 00:00:00 2001 From: Simon Praetorius <simon.praetorius@tu-dresden.de> Date: Thu, 30 Aug 2018 17:43:25 +0200 Subject: [PATCH] Initial version of timeseries writer --- dune/vtk/vtktimeserieswriter.hh | 67 +++++++ dune/vtk/vtktimeserieswriter.impl.hh | 84 +++++++++ dune/vtk/vtkwriterinterface.hh | 38 ++-- dune/vtk/vtkwriterinterface.impl.hh | 16 +- dune/vtk/writers/vtkimagedatawriter.hh | 6 +- dune/vtk/writers/vtkrectilineargridwriter.hh | 6 +- dune/vtk/writers/vtkstructuredgridwriter.hh | 6 +- dune/vtk/writers/vtkunstructuredgridwriter.hh | 14 +- .../writers/vtkunstructuredgridwriter.impl.hh | 164 ++++++++++++++++++ src/CMakeLists.txt | 4 + src/timeserieswriter.cc | 59 +++++++ 11 files changed, 431 insertions(+), 33 deletions(-) create mode 100644 dune/vtk/vtktimeserieswriter.hh create mode 100644 dune/vtk/vtktimeserieswriter.impl.hh create mode 100644 src/timeserieswriter.cc diff --git a/dune/vtk/vtktimeserieswriter.hh b/dune/vtk/vtktimeserieswriter.hh new file mode 100644 index 0000000..78feb8f --- /dev/null +++ b/dune/vtk/vtktimeserieswriter.hh @@ -0,0 +1,67 @@ +#pragma once + +#include <iosfwd> +#include <map> +#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 + template <class VtkWriter> + class VtkTimeseriesWriter + { + protected: + using Self = VtkTimeseriesWriter; + using pos_type = typename std::ostream::pos_type; + + public: + /// Constructor, stores the gridView + template <class... Args, disableCopyMove<Self, Args...> = 0> + VtkTimeseriesWriter (Args&&... args) + : vtkWriter_{std::forward<Args>(args)...} + {} + + /// Write the attached data to the file with \ref Vtk::FormatTypes and \ref Vtk::DataTypes + void writeTimestep (double time, std::string const& fn); + + // NOTE: requires a aforging call to \ref write + void write (std::string const& fn); + + /// 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_; + + bool initialized_ = false; + + // block size of attached data + std::vector<std::uint64_t> blocksize_; + + 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 0000000..cc20ead --- /dev/null +++ b/dune/vtk/vtktimeserieswriter.impl.hh @@ -0,0 +1,84 @@ +#pragma once + +#include <algorithm> +#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 p = filesystem::path(fn); + auto name = p.stem(); + p.remove_filename(); + p /= name.string(); + + vtkWriter_.dataCollector_.update(); + + if (!initialized_) { + // write points and cells only once + filenameMesh_ = p.string() + ".mesh.data"; + std::ofstream file_mesh(filenameMesh_, std::ios_base::ate | std::ios::binary); + + if (vtkWriter_.datatype_ == Vtk::FLOAT32) + blocksize_.push_back( vtkWriter_.template writePointsAppended<float>(file_mesh) ); + else + blocksize_.push_back( vtkWriter_.template writePointsAppended<double>(file_mesh) ); + + auto bs = vtkWriter_.writeCellsAppended(file_mesh); + blocksize_.insert(blocksize_.end(), bs.begin(), bs.end()); + initialized_ = true; + + std::cout << "blocksize = [" << join(blocksize_.begin(), blocksize_.end()) << "]\n"; + } + + std::string filenameData = p.string() + "_t" + std::to_string(timesteps_.size()) + ".data"; + std::ofstream file_data(filenameData, std::ios_base::ate | std::ios::binary); + + for (auto const& v : vtkWriter_.pointData_) { + if (v.type() == Vtk::FLOAT32) + blocksize_.push_back( vtkWriter_.template writeDataAppended<float>(file_data, v, W::POINT_DATA) ); + else + blocksize_.push_back( vtkWriter_.template writeDataAppended<double>(file_data, v, W::POINT_DATA) ); + } + for (auto const& v : vtkWriter_.cellData_) { + if (v.type() == Vtk::FLOAT32) + blocksize_.push_back( vtkWriter_.template writeDataAppended<float>(file_data, v, W::CELL_DATA) ); + else + blocksize_.push_back( vtkWriter_.template writeDataAppended<double>(file_data, v, W::CELL_DATA) ); + } + + timesteps_.emplace_back(time, filenameData); +} + + +template <class W> +void VtkTimeseriesWriter<W> + ::write (std::string const& fn) +{ + auto p = filesystem::path(fn); + auto name = p.stem(); + p.remove_filename(); + p /= name.string(); + std::string filename = p.string() + "." + vtkWriter_.getFileExtension(); + vtkWriter_.writeTimeseriesFile(filename, filenameMesh_, timesteps_, blocksize_); +} + +} // end namespace Dune diff --git a/dune/vtk/vtkwriterinterface.hh b/dune/vtk/vtkwriterinterface.hh index dab897b..221a96d 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> @@ -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 26883e3..d36dafc 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); diff --git a/dune/vtk/writers/vtkimagedatawriter.hh b/dune/vtk/writers/vtkimagedatawriter.hh index d3b2af6..0db3f86 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: diff --git a/dune/vtk/writers/vtkrectilineargridwriter.hh b/dune/vtk/writers/vtkrectilineargridwriter.hh index b40fdab..f7df76b 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: diff --git a/dune/vtk/writers/vtkstructuredgridwriter.hh b/dune/vtk/writers/vtkstructuredgridwriter.hh index beb9d56..b27d890 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: diff --git a/dune/vtk/writers/vtkunstructuredgridwriter.hh b/dune/vtk/writers/vtkunstructuredgridwriter.hh index 221cb39..f513923 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,6 +46,12 @@ 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 + void writeTimeseriesFile (std::string const& filename, + std::string const& filenameMesh, + std::vector<std::pair<double, std::string>> const& timesteps, + std::vector<std::uint64_t> const& blocksize) const; + virtual std::string fileExtension () const override { return "vtu"; diff --git a/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh b/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh index 42769b6..5abec03 100644 --- a/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh +++ b/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh @@ -171,6 +171,170 @@ void VtkUnstructuredGridWriter<GV,DC> } +template <class GV, class DC> +void VtkUnstructuredGridWriter<GV,DC> + ::writeTimeseriesFile (std::string const& filename, + std::string const& filenameMesh, + std::vector<std::pair<double, std::string>> const& timesteps, + std::vector<std::uint64_t> const& blocksize) 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=\"\n"; + for (auto const& timestep : timesteps) + out << timestep.first << "\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) { + out << "<DataArray" + << " type=\"" << to_string(datatype_) << "\"" + << " NumberOfComponents=\"3\"" + << " TimeStep=\"" << i << "\"" + << " format=\"appended\"" + << " offset="; + offsets[i].push_back(out.tellp()); + out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' '); + out << "/>\n"; + } + out << "</Points>\n"; + + // Write element connectivity, types and offsets + out << "<Cells>\n"; + for (std::size_t i = 0; i < timesteps.size(); ++i) { + out << "<DataArray type=\"Int64\" Name=\"connectivity\" format=\"appended\"" + << " TimeStep=\"" << i << "\"" + << " offset="; + offsets[i].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\"" + << " TimeStep=\"" << i << "\"" + << " offset="; + offsets[i].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\"" + << " TimeStep=\"" << i << "\"" + << " offset="; + offsets[i].push_back(out.tellp()); + out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' '); + out << "/>\n"; + } + out << "</Cells>\n"; + + const std::size_t shift = offsets[0].size(); // number of blocks to write the grid + std::cout << "shift = " << shift << "\n"; + + // 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_) { + out << "<DataArray" + << " Name=\"" << v.name() << "\"" + << " TimeStep=\"" << i << "\"" + << " type=\"" << to_string(datatype_) << "\"" + << " NumberOfComponents=\"" << v.ncomps() << "\"" + << " format=\"appended\"" + << " offset="; + offsets[i].push_back(out.tellp()); + out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' '); + out << "/>\n"; + } + } + 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_) { + out << "<DataArray" + << " Name=\"" << v.name() << "\"" + << " TimeStep=\"" << i << "\"" + << " type=\"" << to_string(datatype_) << "\"" + << " NumberOfComponents=\"" << v.ncomps() << "\"" + << " format=\"appended\"" + << " offset="; + offsets[i].push_back(out.tellp()); + out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' '); + out << "/>\n"; + } + } + out << "</CellData>\n"; + + out << "</Piece>\n"; + out << "</UnstructuredGrid>\n"; + + std::vector<std::uint64_t> blocks; // size of i'th appended block + pos_type appended_pos = 0; + out << "<AppendedData encoding=\"raw\">\n_"; + appended_pos = out.tellp(); + + std::ifstream file_mesh(filenameMesh, std::ios_base::in | std::ios_base::binary); + out << file_mesh.rdbuf(); + file_mesh.close(); + assert( std::uint64_t(out.tellp()) == std::accumulate(blocksize.begin(), std::next(blocksize.begin(),shift), std::uint64_t(appended_pos)) ); + + 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(blocksize[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(blocksize[j++]); + } + } +} + + template <class GV, class DC> void VtkUnstructuredGridWriter<GV,DC> ::writeCells (std::ofstream& out, std::vector<pos_type>& offsets) const diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 46958c1..3b6012a 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/timeserieswriter.cc b/src/timeserieswriter.cc new file mode 100644 index 0000000..35f2752 --- /dev/null +++ b/src/timeserieswriter.cc @@ -0,0 +1,59 @@ +// -*- 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/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) +{ + 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>; + VtkTimeseriesWriter<Writer> vtkWriter(gridView, Vtk::BINARY, Vtk::FLOAT32); + vtkWriter.addPointData(p1Analytic, "q1"); + vtkWriter.addCellData(p1Analytic, "q0"); + for (double t = 0.0; t < 5; t += 0.5) { + vtkWriter.writeTimestep(t, prefix + "_" + std::to_string(GridView::dimension) + "d_ascii.vtu"); + } + vtkWriter.write(prefix + "_" + std::to_string(GridView::dimension) + "d_ascii.vtu"); +} + +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>(4); + GridType grid(upperRight, numElements); + write("yasp", grid.leafGridView()); +} \ No newline at end of file -- GitLab