Commit 8998ec6e authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Initial version of timeseries writer

parent cecca335
#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"
#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
......@@ -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_;
......
......@@ -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);
......
......@@ -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:
......
......@@ -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:
......
......@@ -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:
......
......@@ -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";
......
......@@ -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
......
......@@ -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)
......
// -*- 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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment