From 247219385cabdf1c58d1515a10ba7c4e4de41b28 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Thu, 23 Aug 2018 16:45:22 +0200
Subject: [PATCH] added structured grid writer
---
.../datacollectors/structureddatacollector.hh | 156 +++++++++++++++
dune/vtk/vtkimagedatawriter.hh | 58 ++++++
dune/vtk/vtkimagedatawriter.impl.hh | 183 ++++++++++++++++++
dune/vtk/vtkstructuredgridwriter.hh | 58 ++++++
dune/vtk/vtkstructuredgridwriter.impl.hh | 181 +++++++++++++++++
src/CMakeLists.txt | 4 +
src/structuredgridwriter.cc | 99 ++++++++++
7 files changed, 739 insertions(+)
create mode 100644 dune/vtk/datacollectors/structureddatacollector.hh
create mode 100644 dune/vtk/vtkimagedatawriter.hh
create mode 100644 dune/vtk/vtkimagedatawriter.impl.hh
create mode 100644 dune/vtk/vtkstructuredgridwriter.hh
create mode 100644 dune/vtk/vtkstructuredgridwriter.impl.hh
create mode 100644 src/structuredgridwriter.cc
diff --git a/dune/vtk/datacollectors/structureddatacollector.hh b/dune/vtk/datacollectors/structureddatacollector.hh
new file mode 100644
index 0000000..80b54dd
--- /dev/null
+++ b/dune/vtk/datacollectors/structureddatacollector.hh
@@ -0,0 +1,156 @@
+#pragma once
+
+#include <dune/vtk/datacollector.hh>
+
+namespace Dune { namespace experimental
+{
+
+/// Implementation of \ref DataCollector for linear cells, with continuous data.
+template <class GridView>
+class StructuredDataCollector
+ : public DataCollectorInterface<GridView, StructuredDataCollector<GridView>>
+{
+ enum { dim = GridView::dimension };
+
+ using Self = StructuredDataCollector;
+ using Super = DataCollectorInterface<GridView, Self>;
+ using Super::gridView_;
+ using ctype = typename GridView::ctype;
+
+ struct CoordLess
+ {
+ template <class T, int N>
+ bool operator() (FieldVector<T,N> const& lhs, FieldVector<T,N> const& rhs) const
+ {
+ for (int i = N-1; i >= 0; --i) {
+ if (std::abs(lhs[i] - rhs[i]) < std::numeric_limits<T>::epsilon())
+ continue;
+ return lhs[i] < rhs[i];
+ }
+ return false;
+ }
+ };
+
+public:
+ StructuredDataCollector (GridView const& gridView)
+ : Super(gridView)
+ {}
+
+ /// Return number of grid vertices
+ std::uint64_t numPointsImpl () const
+ {
+ return gridView_.size(dim);
+ }
+
+ std::array<int, 6> const& extent () const
+ {
+ return extent_;
+ }
+
+ auto const& origin () const
+ {
+ return origin_;
+ }
+
+ auto const& spacing () const
+ {
+ return spacing_;
+ }
+
+ void updateImpl ()
+ {
+ // TODO: extract this information from the grid
+ int extent = GridView::dimension == 1 ? gridView_.size(dim) :
+ GridView::dimension == 2 ? isqrt(gridView_.size(dim)) :
+ GridView::dimension == 3 ? icbrt(gridView_.size(dim)) : 0;
+ for (int i = 0; i < 3; ++i) {
+ if (GridView::dimension > i) {
+ extent_[2*i] = 0;
+ extent_[2*i+1] = extent-1;
+ spacing_[i] = gridView_.grid().domainSize()[i] / (extent-1);
+ } else {
+ extent_[2*i] = extent_[2*i+1] = 0;
+ spacing_[i] = 0;
+ }
+
+ origin_[i] = 0;
+ }
+ }
+
+ /// Return the coordinates of all grid vertices in the order given by the indexSet
+ template <class T>
+ std::vector<T> pointsImpl () const
+ {
+ std::vector<T> data(this->numPoints() * 3);
+ auto const& indexSet = gridView_.indexSet();
+ for (auto const& vertex : vertices(gridView_)) {
+ std::size_t idx = 3 * indexSet.index(vertex);
+ auto v = vertex.geometry().center();
+ for (std::size_t j = 0; j < v.size(); ++j)
+ data[idx + j] = T(v[j]);
+ for (std::size_t j = v.size(); j < 3u; ++j)
+ data[idx + j] = T(0);
+ }
+ return data;
+ }
+
+ /// Evaluate the `fct` at the corners of the elements
+ template <class T, class GlobalFunction>
+ std::vector<T> pointDataImpl (GlobalFunction const& fct) const
+ {
+ std::vector<T> data(this->numPoints() * fct.ncomps());
+ auto const& indexSet = gridView_.indexSet();
+ auto localFct = localFunction(fct);
+ for (auto const& e : elements(gridView_)) {
+ localFct.bind(e);
+ Vtk::CellType cellType{e.type()};
+ auto refElem = referenceElement(e.geometry());
+ for (int j = 0; j < e.subEntities(dim); ++j) {
+ int k = cellType.permutation(j);
+ std::size_t idx = fct.ncomps() * indexSet.subIndex(e,k,dim);
+ for (int comp = 0; comp < fct.ncomps(); ++comp)
+ data[idx + comp] = T(localFct.evaluate(comp, refElem.position(k,dim)));
+ }
+ localFct.unbind();
+ }
+ return data;
+ }
+
+private:
+ static constexpr std::uint32_t isqrt (std::uint64_t x)
+ {
+ std::uint32_t c = 0x8000;
+ std::uint32_t g = 0x8000;
+
+ while (true) {
+ if (g*g > x)
+ g ^= c;
+ c >>= 1;
+ if (c == 0)
+ return g;
+ g |= c;
+ }
+ }
+
+ static constexpr std::uint32_t icbrt (std::uint64_t x)
+ {
+ std::uint32_t y = 0;
+ for (int s = 63; s >= 0; s -= 3) {
+ y += y;
+ std::uint64_t b = 3*y*(std::uint64_t(y) + 1) + 1;
+ if ((x >> s) >= b) {
+ x -= b << s;
+ y++;
+ }
+ }
+ return y;
+ }
+
+private:
+ std::array<int, 6> extent_;
+ FieldVector<ctype,3> spacing_;
+ FieldVector<ctype,3> origin_;
+ std::vector<std::size_t> indexMap_;
+};
+
+}} // end namespace Dune::experimental
diff --git a/dune/vtk/vtkimagedatawriter.hh b/dune/vtk/vtkimagedatawriter.hh
new file mode 100644
index 0000000..b1e4345
--- /dev/null
+++ b/dune/vtk/vtkimagedatawriter.hh
@@ -0,0 +1,58 @@
+#pragma once
+
+#include <array>
+#include <iosfwd>
+#include <map>
+
+#include "datacollector.hh"
+#include "filewriter.hh"
+#include "vtkfunction.hh"
+#include "vtktypes.hh"
+#include "vtkwriter.hh"
+#include "datacollectors/structureddatacollector.hh"
+
+namespace Dune { namespace experimental
+{
+ /// File-Writer for VTK .vtu files
+ template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
+ class VtkImageDataWriter
+ : public VtkWriter<GridView, DataCollector>
+ {
+ static constexpr int dimension = GridView::dimension;
+
+ using Super = VtkWriter<GridView, DataCollector>;
+ using pos_type = typename Super::pos_type;
+
+ public:
+ /// Constructor, stores the gridView
+ VtkImageDataWriter (GridView const& gridView)
+ : Super(gridView)
+ {}
+
+ private:
+ /// Write a serial VTK file in Unstructured format
+ virtual void writeSerialFile (std::string const& filename) const override;
+
+ /// Write a parallel VTK file `pfilename.pvtu` in Unstructured format,
+ /// with `size` the number of pieces and serial files given by `pfilename_p[i].vtu`
+ /// for [i] in [0,...,size).
+ virtual void writeParallelFile (std::string const& pfilename, int size) const override;
+
+ virtual std::string fileExtension () const override
+ {
+ return "vti";
+ }
+
+ private:
+ using Super::dataCollector_;
+ using Super::format_;
+ using Super::datatype_;
+
+ // attached data
+ using Super::pointData_;
+ using Super::cellData_;
+ };
+
+}} // end namespace Dune::experimental
+
+#include "vtkimagedatawriter.impl.hh"
diff --git a/dune/vtk/vtkimagedatawriter.impl.hh b/dune/vtk/vtkimagedatawriter.impl.hh
new file mode 100644
index 0000000..0f0ac77
--- /dev/null
+++ b/dune/vtk/vtkimagedatawriter.impl.hh
@@ -0,0 +1,183 @@
+#pragma once
+
+#include <iomanip>
+#include <iostream>
+#include <iterator>
+#include <fstream>
+#include <sstream>
+#include <string>
+
+#include <dune/geometry/referenceelements.hh>
+#include <dune/geometry/type.hh>
+
+#include "utility/enum.hh"
+#include "utility/filesystem.hh"
+#include "utility/string.hh"
+
+namespace Dune { namespace experimental {
+
+template <class GV, class DC>
+void VtkImageDataWriter<GV,DC>
+ ::writeSerialFile (std::string const& filename) const
+{
+ std::ofstream out(filename, std::ios_base::ate | std::ios::binary);
+ if (format_ == Vtk::ASCII) {
+ if (datatype_ == Vtk::FLOAT32)
+ out << std::setprecision(std::numeric_limits<float>::digits10+2);
+ else
+ out << std::setprecision(std::numeric_limits<double>::digits10+2);
+ }
+
+ dataCollector_.update();
+
+ std::vector<pos_type> offsets; // pos => offset
+ out << "<VTKFile type=\"ImageData\" version=\"1.0\" "
+ << "byte_order=\"" << this->getEndian() << "\" header_type=\"UInt64\""
+ << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\">\n" : ">\n");
+ out << "<ImageData WholeExtent=\"";
+ auto const& extent = dataCollector_.extent();
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
+ }
+ out << "\" Origin=\"";
+ auto const& origin = dataCollector_.origin();
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << origin[i];
+ }
+ out << "\" Spacing=\"";
+ auto const& spacing = dataCollector_.spacing();
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << spacing[i];
+ }
+ out << "\">\n";
+ out << "<Piece Extent=\"";
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
+ }
+ out << "\">\n";
+
+ { // Write data associated with grid points
+ auto scalar = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+ auto vector = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+ out << "<PointData" << (scalar != pointData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+ << (vector != pointData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+ << ">\n";
+ for (auto const& v : pointData_)
+ this->writeData(out, offsets, v, Super::POINT_DATA);
+ out << "</PointData>\n";
+ }
+
+ { // Write data associated with grid cells
+ auto scalar = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+ auto vector = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+ out << "<CellData" << (scalar != cellData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+ << (vector != cellData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+ << ">\n";
+ for (auto const& v : cellData_)
+ this->writeData(out, offsets, v, Super::CELL_DATA);
+ out << "</CellData>\n";
+ }
+ 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 (datatype_ == 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 (datatype_ == 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";
+ }
+
+ 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]);
+ }
+ }
+}
+
+
+template <class GV, class DC>
+void VtkImageDataWriter<GV,DC>
+ ::writeParallelFile (std::string const& pfilename, int size) const
+{
+ std::string filename = pfilename + ".p" + this->fileExtension();
+ std::ofstream out(filename, std::ios_base::ate | std::ios::binary);
+
+ out << "<VTKFile type=\"PImageData\" version=\"1.0\" "
+ << "byte_order=\"" << this->getEndian() << "\" header_type=\"UInt64\""
+ << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\">\n" : ">\n");
+ out << "<PImageData GhostLevel=\"0\" WholeExtent=\"";
+ auto const& extent = dataCollector_.extent();
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
+ }
+ out << "\" Origin=\"";
+ auto const& origin = dataCollector_.origin();
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << origin[i];
+ }
+ out << "\" Spacing=\"";
+ auto const& spacing = dataCollector_.spacing();
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << spacing[i];
+ }
+ out << "\">\n";
+
+ { // Write data associated with grid points
+ auto scalar = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+ auto vector = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+ out << "<PPointData" << (scalar != pointData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+ << (vector != pointData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+ << ">\n";
+ for (auto const& v : pointData_) {
+ out << "<PDataArray Name=\"" << v.name() << "\" type=\"" << Vtk::Map::from_datatype[datatype_] << "\""
+ << " NumberOfComponents=\"" << v.ncomps() << "\" />\n";
+ }
+ out << "</PPointData>\n";
+ }
+
+ { // Write data associated with grid cells
+ auto scalar = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+ auto vector = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+ out << "<PCellData" << (scalar != cellData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+ << (vector != cellData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+ << ">\n";
+ for (auto const& v : cellData_) {
+ out << "<PDataArray Name=\"" << v.name() << "\" type=\"" << Vtk::Map::from_datatype[datatype_] << "\""
+ << " NumberOfComponents=\"" << v.ncomps() << "\" />\n";
+ }
+ out << "</PCellData>\n";
+ }
+
+ // Write piece file references
+ for (int i = 0; i < size; ++i) {
+ out << "<Piece Source=\"" << pfilename << "_p" << std::to_string(i) << "." << this->fileExtension() << "\" Extent=\"";
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
+ }
+ out << "\" />\n";
+ }
+
+ out << "</PUnstructuredGrid>\n";
+ out << "</VTKFile>";
+}
+
+}} // end namespace Dune::experimental
diff --git a/dune/vtk/vtkstructuredgridwriter.hh b/dune/vtk/vtkstructuredgridwriter.hh
new file mode 100644
index 0000000..8c3aaa2
--- /dev/null
+++ b/dune/vtk/vtkstructuredgridwriter.hh
@@ -0,0 +1,58 @@
+#pragma once
+
+#include <array>
+#include <iosfwd>
+#include <map>
+
+#include "datacollector.hh"
+#include "filewriter.hh"
+#include "vtkfunction.hh"
+#include "vtktypes.hh"
+#include "vtkwriter.hh"
+#include "datacollectors/structureddatacollector.hh"
+
+namespace Dune { namespace experimental
+{
+ /// File-Writer for VTK .vtu files
+ template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
+ class VtkStructuredGridWriter
+ : public VtkWriter<GridView, DataCollector>
+ {
+ static constexpr int dimension = GridView::dimension;
+
+ using Super = VtkWriter<GridView, DataCollector>;
+ using pos_type = typename Super::pos_type;
+
+ public:
+ /// Constructor, stores the gridView
+ VtkStructuredGridWriter (GridView const& gridView)
+ : Super(gridView)
+ {}
+
+ private:
+ /// Write a serial VTK file in Unstructured format
+ virtual void writeSerialFile (std::string const& filename) const override;
+
+ /// Write a parallel VTK file `pfilename.pvtu` in Unstructured format,
+ /// with `size` the number of pieces and serial files given by `pfilename_p[i].vtu`
+ /// for [i] in [0,...,size).
+ virtual void writeParallelFile (std::string const& pfilename, int size) const override;
+
+ virtual std::string fileExtension () const override
+ {
+ return "vts";
+ }
+
+ private:
+ using Super::dataCollector_;
+ using Super::format_;
+ using Super::datatype_;
+
+ // attached data
+ using Super::pointData_;
+ using Super::cellData_;
+ };
+
+}} // end namespace Dune::experimental
+
+#include "vtkstructuredgridwriter.impl.hh"
diff --git a/dune/vtk/vtkstructuredgridwriter.impl.hh b/dune/vtk/vtkstructuredgridwriter.impl.hh
new file mode 100644
index 0000000..7900e33
--- /dev/null
+++ b/dune/vtk/vtkstructuredgridwriter.impl.hh
@@ -0,0 +1,181 @@
+#pragma once
+
+#include <iomanip>
+#include <iostream>
+#include <iterator>
+#include <fstream>
+#include <sstream>
+#include <string>
+
+#include <dune/geometry/referenceelements.hh>
+#include <dune/geometry/type.hh>
+
+#include "utility/enum.hh"
+#include "utility/filesystem.hh"
+#include "utility/string.hh"
+
+namespace Dune { namespace experimental {
+
+template <class GV, class DC>
+void VtkStructuredGridWriter<GV,DC>
+ ::writeSerialFile (std::string const& filename) const
+{
+ std::ofstream out(filename, std::ios_base::ate | std::ios::binary);
+ if (format_ == Vtk::ASCII) {
+ if (datatype_ == Vtk::FLOAT32)
+ out << std::setprecision(std::numeric_limits<float>::digits10+2);
+ else
+ out << std::setprecision(std::numeric_limits<double>::digits10+2);
+ }
+
+ dataCollector_.update();
+
+ std::vector<pos_type> offsets; // pos => offset
+ out << "<VTKFile type=\"StructuredGrid\" version=\"1.0\" "
+ << "byte_order=\"" << this->getEndian() << "\" header_type=\"UInt64\""
+ << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\">\n" : ">\n");
+ out << "<StructuredGrid WholeExtent=\"";
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << dataCollector_.extent()[2*i];
+ out << " " << dataCollector_.extent()[2*i+1];
+ }
+ out << "\">\n";
+ out << "<Piece NumberOfPoints=\"" << dataCollector_.numPoints() << "\" "
+ << "NumberOfCells=\"" << dataCollector_.numCells() << "\" Extent=\"";
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << dataCollector_.extent()[2*i];
+ out << " " << dataCollector_.extent()[2*i+1];
+ }
+ out << "\">\n";
+
+ { // Write data associated with grid points
+ auto scalar = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+ auto vector = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+ out << "<PointData" << (scalar != pointData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+ << (vector != pointData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+ << ">\n";
+ for (auto const& v : pointData_)
+ this->writeData(out, offsets, v, Super::POINT_DATA);
+ out << "</PointData>\n";
+ }
+
+ { // Write data associated with grid cells
+ auto scalar = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+ auto vector = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+ out << "<CellData" << (scalar != cellData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+ << (vector != cellData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+ << ">\n";
+ for (auto const& v : cellData_)
+ 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 (datatype_ == 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 (datatype_ == 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";
+ }
+
+ 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]);
+ }
+ }
+}
+
+
+template <class GV, class DC>
+void VtkStructuredGridWriter<GV,DC>
+ ::writeParallelFile (std::string const& pfilename, int size) const
+{
+ std::string filename = pfilename + ".p" + this->fileExtension();
+ std::ofstream out(filename, std::ios_base::ate | std::ios::binary);
+
+ out << "<VTKFile type=\"PStructuredGrid\" version=\"1.0\" "
+ << "byte_order=\"" << this->getEndian() << "\" header_type=\"UInt64\""
+ << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\">\n" : ">\n");
+ out << "<PStructuredGrid GhostLevel=\"0\" WholeExtent=\"";
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << dataCollector_.extent()[2*i];
+ out << " " << dataCollector_.extent()[2*i+1];
+ } // TODO: WholeExtent is probably a global information. needs synchronization
+ out << "\">\n";
+
+ { // Write data associated with grid points
+ auto scalar = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+ auto vector = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+ out << "<PPointData" << (scalar != pointData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+ << (vector != pointData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+ << ">\n";
+ for (auto const& v : pointData_) {
+ out << "<PDataArray Name=\"" << v.name() << "\" type=\"" << Vtk::Map::from_datatype[datatype_] << "\""
+ << " NumberOfComponents=\"" << v.ncomps() << "\" />\n";
+ }
+ out << "</PPointData>\n";
+ }
+
+ { // Write data associated with grid cells
+ auto scalar = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() == 1; });
+ auto vector = std::find_if(cellData_.begin(), cellData_.end(), [](auto const& v) { return v.ncomps() > 1; });
+ out << "<PCellData" << (scalar != cellData_.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
+ << (vector != cellData_.end() ? " Vectors=\"" + vector->name() + "\"" : "")
+ << ">\n";
+ for (auto const& v : cellData_) {
+ out << "<PDataArray Name=\"" << v.name() << "\" type=\"" << Vtk::Map::from_datatype[datatype_] << "\""
+ << " NumberOfComponents=\"" << v.ncomps() << "\" />\n";
+ }
+ out << "</PCellData>\n";
+ }
+
+ // Write points
+ out << "<PPoints>\n";
+ out << "<PDataArray type=\"" << Vtk::Map::from_datatype[datatype_] << "\" NumberOfComponents=\"3\" />\n";
+ out << "</PPoints>\n";
+
+ // Write piece file references
+ for (int i = 0; i < size; ++i) {
+ out << "<Piece Source=\"" << pfilename << "_p" << std::to_string(i) << "." << this->fileExtension() << "\" Extent=\"";
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << dataCollector_.extent()[2*i];
+ out << " " << dataCollector_.extent()[2*i+1];
+ } // TODO: need extent of piece
+ out << "\" />\n";
+ }
+
+ out << "</PUnstructuredGrid>\n";
+ out << "</VTKFile>";
+}
+
+}} // end namespace Dune::experimental
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index a1d6e57..0978279 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -18,4 +18,8 @@ add_executable("quadratic" quadratic.cc)
target_link_dune_default_libraries("quadratic")
target_link_libraries("quadratic" dunevtk)
+add_executable("structuredgridwriter" structuredgridwriter.cc)
+target_link_dune_default_libraries("structuredgridwriter")
+target_link_libraries("structuredgridwriter" dunevtk)
+
add_subdirectory(test)
\ No newline at end of file
diff --git a/src/structuredgridwriter.cc b/src/structuredgridwriter.cc
new file mode 100644
index 0000000..8d08855
--- /dev/null
+++ b/src/structuredgridwriter.cc
@@ -0,0 +1,99 @@
+// -*- 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/functionspacebases/interpolate.hh>
+#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
+
+#include <dune/vtk/vtkstructuredgridwriter.hh>
+#include <dune/vtk/vtkimagedatawriter.hh>
+
+using namespace Dune;
+using namespace Dune::experimental;
+using namespace Dune::Functions;
+
+#define GRID_TYPE 1
+
+int main(int argc, char** argv)
+{
+ Dune::MPIHelper::instance(argc, argv);
+
+ const int dim = 3;
+
+ using GridType = YaspGrid<dim>;
+ FieldVector<double,dim> upperRight; upperRight = 1.0;
+ auto numElements = filledArray<dim,int>(8);
+ GridType grid(upperRight,numElements);
+
+ using GridView = typename GridType::LeafGridView;
+ GridView gridView = grid.leafGridView();
+
+ using namespace BasisFactory;
+ auto basis = makeBasis(gridView, lagrange<1>());
+
+ std::vector<double> p1function(basis.dimension());
+
+ interpolate(basis, p1function, [](auto const& x) {
+ return 100*x[0] + 10*x[1] + 1*x[2];
+ });
+
+ // write discrete global-basis function
+ auto p1FctWrapped = makeDiscreteGlobalBasisFunction<double>(basis, p1function);
+
+ {
+ using Writer = VtkStructuredGridWriter<GridView>;
+ Writer vtkWriter(gridView);
+ vtkWriter.addPointData(p1FctWrapped, "p1");
+ vtkWriter.addCellData(p1FctWrapped, "p0");
+
+ // write analytic function
+ auto p1Analytic = makeAnalyticGridViewFunction([](auto const& x) {
+ return std::sin(10*x[0])*std::cos(10*x[1])+std::sin(10*x[2]);
+ }, gridView);
+
+ vtkWriter.addPointData(p1Analytic, "analytic");
+
+ vtkWriter.write("sg_ascii_float32.vts", Vtk::ASCII);
+ vtkWriter.write("sg_binary_float32.vts", Vtk::BINARY);
+ vtkWriter.write("sg_compressed_float32.vts", Vtk::COMPRESSED);
+ vtkWriter.write("sg_ascii_float64.vts", Vtk::ASCII, Vtk::FLOAT64);
+ vtkWriter.write("sg_binary_float64.vts", Vtk::BINARY, Vtk::FLOAT64);
+ vtkWriter.write("sg_compressed_float64.vts", Vtk::COMPRESSED, Vtk::FLOAT64);
+ }
+
+ {
+ using Writer = VtkImageDataWriter<GridView>;
+ Writer vtkWriter(gridView);
+ vtkWriter.addPointData(p1FctWrapped, "p1");
+ vtkWriter.addCellData(p1FctWrapped, "p0");
+
+ // write analytic function
+ auto p1Analytic = makeAnalyticGridViewFunction([](auto const& x) {
+ return std::sin(10*x[0])*std::cos(10*x[1])+std::sin(10*x[2]);
+ }, gridView);
+
+ vtkWriter.addPointData(p1Analytic, "analytic");
+
+ vtkWriter.write("id_ascii_float32.vti", Vtk::ASCII);
+ vtkWriter.write("id_binary_float32.vti", Vtk::BINARY);
+ vtkWriter.write("id_compressed_float32.vti", Vtk::COMPRESSED);
+ vtkWriter.write("id_ascii_float64.vti", Vtk::ASCII, Vtk::FLOAT64);
+ vtkWriter.write("id_binary_float64.vti", Vtk::BINARY, Vtk::FLOAT64);
+ vtkWriter.write("id_compressed_float64.vti", Vtk::COMPRESSED, Vtk::FLOAT64);
+ }
+}
--
GitLab