Commit be9d72e6 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Implementation and test of timeseries writer

parent 8998ec6e
#pragma once
#include <iosfwd>
#include <map>
#include <string>
#include <tuple>
#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
......@@ -16,6 +12,7 @@ namespace Dune
/// File-Writer for Vtk .vtu files
template <class VtkWriter>
class VtkTimeseriesWriter
: public FileWriter
{
protected:
using Self = VtkTimeseriesWriter;
......@@ -26,13 +23,16 @@ namespace Dune
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");
}
/// 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);
/// 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>
......
#pragma once
#include <algorithm>
#include <cstdio>
#include <iomanip>
#include <iostream>
#include <iterator>
......@@ -34,7 +35,7 @@ void VtkTimeseriesWriter<W>
if (!initialized_) {
// write points and cells only once
filenameMesh_ = p.string() + ".mesh.data";
filenameMesh_ = p.string() + ".mesh.vtkdata";
std::ofstream file_mesh(filenameMesh_, std::ios_base::ate | std::ios::binary);
if (vtkWriter_.datatype_ == Vtk::FLOAT32)
......@@ -45,11 +46,9 @@ void VtkTimeseriesWriter<W>
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::string filenameData = p.string() + "_t" + std::to_string(timesteps_.size()) + ".vtkdata";
std::ofstream file_data(filenameData, std::ios_base::ate | std::ios::binary);
for (auto const& v : vtkWriter_.pointData_) {
......@@ -73,12 +72,23 @@ 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 filename = p.string() + "." + vtkWriter_.getFileExtension();
std::string filename = p.string() + ".ts." + vtkWriter_.getFileExtension();
vtkWriter_.writeTimeseriesFile(filename, filenameMesh_, timesteps_, blocksize_);
// 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
......@@ -197,9 +197,12 @@ void VtkUnstructuredGridWriter<GV,DC>
<< ">\n";
out << "<UnstructuredGrid"
<< " TimeValues=\"\n";
<< " TimeValues=\"";
{
std::size_t i = 0;
for (auto const& timestep : timesteps)
out << timestep.first << "\n";
out << timestep.first << (++i % 6 != 0 ? ' ' : '\n');
}
out << "\">\n";
out << "<Piece"
......@@ -249,7 +252,6 @@ void VtkUnstructuredGridWriter<GV,DC>
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";
......
......@@ -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";
}
}
......
......@@ -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>
......
......@@ -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)
......
......@@ -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");
}
......@@ -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");
}
......@@ -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");
}
}
......
......@@ -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));
}
......
......@@ -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"));
}
......
......@@ -12,8 +12,6 @@
#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>
......@@ -26,25 +24,23 @@ 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);
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> vtkWriter(gridView, Vtk::BINARY, Vtk::FLOAT32);
vtkWriter.addPointData(p1Analytic, "q1");
vtkWriter.addCellData(p1Analytic, "q0");
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) {
vtkWriter.writeTimestep(t, prefix + "_" + std::to_string(GridView::dimension) + "d_ascii.vtu");
seriesWriter.writeTimestep(t, filename);
}
vtkWriter.write(prefix + "_" + std::to_string(GridView::dimension) + "d_ascii.vtu");
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)
......
......@@ -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");
}
}
......@@ -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);
for (auto const& test_case : test_cases) {
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");
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.write(prefix + "_" + std::to_string(GridView::dimension) + "d_" + std::get<0>(test_case) + ".vtu");
}
}
......
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