From 390d3b97eb8d4e6641cc2ad2af31e36afc74f005 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Fri, 24 Aug 2018 18:04:40 +0200
Subject: [PATCH] Structured grid writer and data collector extended to support
YaspGrid and SPGrid
---
dune.module | 1 +
dune/vtk/datacollector.hh | 2 +-
.../datacollectors/continuousdatacollector.hh | 6 +-
.../discontinuousdatacollector.hh | 8 +-
.../datacollectors/quadraticdatacollector.hh | 6 +-
dune/vtk/datacollectors/spdatacollector.hh | 102 +++++++++++++
.../datacollectors/structureddatacollector.hh | 140 ++++--------------
dune/vtk/datacollectors/yaspdatacollector.hh | 124 ++++++++++++++++
dune/vtk/vtkimagedatawriter.hh | 2 +-
dune/vtk/vtkimagedatawriter.impl.hh | 47 +++---
dune/vtk/vtkstructuredgridwriter.hh | 3 +-
dune/vtk/vtkstructuredgridwriter.impl.hh | 50 ++++---
dune/vtk/vtkunstructuredgridwriter.impl.hh | 2 -
dune/vtk/vtkwriter.impl.hh | 2 +
src/structuredgridwriter.cc | 124 ++++++++++------
15 files changed, 403 insertions(+), 216 deletions(-)
create mode 100644 dune/vtk/datacollectors/spdatacollector.hh
create mode 100644 dune/vtk/datacollectors/yaspdatacollector.hh
diff --git a/dune.module b/dune.module
index d2413eb..0016381 100644
--- a/dune.module
+++ b/dune.module
@@ -8,3 +8,4 @@ Version: 0.1
Maintainer: simon.praetorius@tu-dresden.de
#depending on
Depends: dune-grid dune-functions
+Suggests: dune-spgrid
diff --git a/dune/vtk/datacollector.hh b/dune/vtk/datacollector.hh
index bfbae79..a5bd9fc 100644
--- a/dune/vtk/datacollector.hh
+++ b/dune/vtk/datacollector.hh
@@ -113,7 +113,7 @@ protected: // default implementations
std::vector<T> data(numCells() * fct.ncomps());
auto const& indexSet = gridView_.indexSet();
auto localFct = localFunction(fct);
- for (auto const& e : elements(gridView_)) {
+ for (auto const& e : elements(gridView_, Partitions::all)) {
localFct.bind(e);
auto geometry = e.geometry();
std::size_t idx = fct.ncomps() * indexSet.index(e);
diff --git a/dune/vtk/datacollectors/continuousdatacollector.hh b/dune/vtk/datacollectors/continuousdatacollector.hh
index dea7840..340f558 100644
--- a/dune/vtk/datacollectors/continuousdatacollector.hh
+++ b/dune/vtk/datacollectors/continuousdatacollector.hh
@@ -33,7 +33,7 @@ public:
{
std::vector<T> data(this->numPoints() * 3);
auto const& indexSet = gridView_.indexSet();
- for (auto const& vertex : vertices(gridView_)) {
+ for (auto const& vertex : vertices(gridView_, Partitions::all)) {
std::size_t idx = 3 * indexSet.index(vertex);
auto v = vertex.geometry().center();
for (std::size_t j = 0; j < v.size(); ++j)
@@ -61,7 +61,7 @@ public:
cells.types.reserve(this->numCells());
std::int64_t old_o = 0;
- for (auto const& c : elements(gridView_)) {
+ for (auto const& c : elements(gridView_, Partitions::all)) {
Vtk::CellType cellType(c.type());
for (int j = 0; j < c.subEntities(dim); ++j)
cells.connectivity.push_back( std::int64_t(indexSet.subIndex(c,cellType.permutation(j),dim)) );
@@ -79,7 +79,7 @@ public:
std::vector<T> data(this->numPoints() * fct.ncomps());
auto const& indexSet = gridView_.indexSet();
auto localFct = localFunction(fct);
- for (auto const& e : elements(gridView_)) {
+ for (auto const& e : elements(gridView_, Partitions::all)) {
localFct.bind(e);
Vtk::CellType cellType{e.type()};
auto refElem = referenceElement(e.geometry());
diff --git a/dune/vtk/datacollectors/discontinuousdatacollector.hh b/dune/vtk/datacollectors/discontinuousdatacollector.hh
index 80ecfe0..88d60a3 100644
--- a/dune/vtk/datacollectors/discontinuousdatacollector.hh
+++ b/dune/vtk/datacollectors/discontinuousdatacollector.hh
@@ -30,7 +30,7 @@ public:
indexMap_.resize(gridView_.size(dim));
std::int64_t vertex_idx = 0;
auto const& indexSet = gridView_.indexSet();
- for (auto const& c : elements(gridView_)) {
+ for (auto const& c : elements(gridView_, Partitions::interior)) {
numPoints_ += c.subEntities(dim);
for (int i = 0; i < c.subEntities(dim); ++i)
indexMap_[indexSet.subIndex(c, i, dim)] = vertex_idx++;
@@ -49,7 +49,7 @@ public:
{
std::vector<T> data(this->numPoints() * 3);
auto const& indexSet = gridView_.indexSet();
- for (auto const& element : elements(gridView_)) {
+ for (auto const& element : elements(gridView_, Partitions::interior)) {
for (int i = 0; i < element.subEntities(dim); ++i) {
std::size_t idx = 3 * indexMap_[indexSet.subIndex(element, i, dim)];
auto v = element.geometry().corner(i);
@@ -72,7 +72,7 @@ public:
std::int64_t old_o = 0;
auto const& indexSet = gridView_.indexSet();
- for (auto const& c : elements(gridView_)) {
+ for (auto const& c : elements(gridView_, Partitions::interior)) {
Vtk::CellType cellType(c.type());
for (int j = 0; j < c.subEntities(dim); ++j) {
std::int64_t vertex_idx = indexMap_[indexSet.subIndex(c,cellType.permutation(j),dim)];
@@ -92,7 +92,7 @@ public:
std::vector<T> data(this->numPoints() * fct.ncomps());
auto const& indexSet = gridView_.indexSet();
auto localFct = localFunction(fct);
- for (auto const& e : elements(gridView_)) {
+ for (auto const& e : elements(gridView_, Partitions::interior)) {
localFct.bind(e);
Vtk::CellType cellType{e.type()};
auto refElem = referenceElement(e.geometry());
diff --git a/dune/vtk/datacollectors/quadraticdatacollector.hh b/dune/vtk/datacollectors/quadraticdatacollector.hh
index 7308158..df3aeb3 100644
--- a/dune/vtk/datacollectors/quadraticdatacollector.hh
+++ b/dune/vtk/datacollectors/quadraticdatacollector.hh
@@ -37,7 +37,7 @@ public:
{
std::vector<T> data(this->numPoints() * 3);
auto const& indexSet = gridView_.indexSet();
- for (auto const& element : elements(gridView_)) {
+ for (auto const& element : elements(gridView_, Partitions::interior)) {
auto geometry = element.geometry();
auto refElem = referenceElement<T,dim>(element.type());
@@ -77,7 +77,7 @@ public:
std::int64_t old_o = 0;
auto const& indexSet = gridView_.indexSet();
- for (auto const& c : elements(gridView_)) {
+ for (auto const& c : elements(gridView_, Partitions::interior)) {
Vtk::CellType cellType(c.type(), Vtk::QUADRATIC);
for (int j = 0; j < c.subEntities(dim); ++j) {
int k = cellType.permutation(j);
@@ -103,7 +103,7 @@ public:
std::vector<T> data(this->numPoints() * fct.ncomps());
auto const& indexSet = gridView_.indexSet();
auto localFct = localFunction(fct);
- for (auto const& e : elements(gridView_)) {
+ for (auto const& e : elements(gridView_, Partitions::interior)) {
localFct.bind(e);
Vtk::CellType cellType{e.type(), Vtk::QUADRATIC};
auto refElem = referenceElement(e.geometry());
diff --git a/dune/vtk/datacollectors/spdatacollector.hh b/dune/vtk/datacollectors/spdatacollector.hh
new file mode 100644
index 0000000..5ec8d53
--- /dev/null
+++ b/dune/vtk/datacollectors/spdatacollector.hh
@@ -0,0 +1,102 @@
+#pragma once
+
+#ifdef HAVE_DUNE_SPGRID
+#include <dune/grid/spgrid.hh>
+#endif
+#include <dune/vtk/datacollectors/structureddatacollector.hh>
+
+namespace Dune { namespace experimental
+{
+#ifdef HAVE_DUNE_SPGRID
+
+// Specialization for SPGrid
+template <class GridView>
+class SPDataCollector
+ : public StructuredDataCollectorInterface<GridView, SPDataCollector<GridView>>
+{
+ enum { dim = GridView::dimension };
+
+ using Self = SPDataCollector;
+ using Super = StructuredDataCollectorInterface<GridView, Self>;
+ using Super::gridView_;
+ using ctype = typename GridView::ctype;
+
+public:
+ SPDataCollector (GridView const& gridView)
+ : Super(gridView)
+ , wholeExtent_(filledArray<6,int>(0))
+ , origin_(0.0)
+ , spacing_(0.0)
+ {}
+
+ std::array<int, 6> const& wholeExtentImpl () const
+ {
+ return wholeExtent_;
+ }
+
+ auto const& originImpl () const
+ {
+ return origin_;
+ }
+
+ auto const& spacingImpl () const
+ {
+ return spacing_;
+ }
+
+ void updateImpl ()
+ {
+ Super::updateImpl();
+
+ auto const& begin = gridView_.impl().gridLevel().globalMesh().begin();
+ auto const& end = gridView_.impl().gridLevel().globalMesh().end();
+ auto const& cube = gridView_.grid().domain().cube();
+
+ for (int i = 0; i < dim; ++i) {
+ wholeExtent_[2*i] = begin[i];
+ wholeExtent_[2*i+1] = end[i];
+ spacing_[i] = cube.width()[i] / (end[i] - begin[i]);
+ origin_[i] = cube.origin()[i];
+ }
+ }
+
+ template <class Writer>
+ void writeLocalPieceImpl (Writer const& writer) const
+ {
+ auto extent = filledArray<6,int>(0);
+ auto const& localMesh = gridView_.impl().gridLevel().localMesh();
+ for (int i = 0; i < dim; ++i) {
+ extent[2*i] = localMesh.begin()[i];
+ extent[2*i+1] = localMesh.end()[i];
+ }
+ writer(extent);
+ }
+
+ template <class Writer>
+ void writePiecesImpl (Writer const& writer) const
+ {
+ auto extent = filledArray<6,int>(0);
+ auto const& partitionList = gridView_.impl().gridLevel().decomposition_;
+ // for (auto const& part : partitionList)
+ for (std::size_t p = 0; p < partitionList.size(); ++p)
+ {
+ auto const& part = partitionList[p];
+ for (int i = 0; i < dim; ++i) {
+ extent[2*i] = part.begin()[i];
+ extent[2*i+1] = part.begin()[i] + part.width(i);
+ }
+ writer(p, extent, true);
+ }
+ }
+
+private:
+ std::array<int, 6> wholeExtent_;
+ std::array<int, 6> extent_;
+ FieldVector<ctype,3> spacing_;
+ FieldVector<ctype,3> origin_;
+ std::vector<std::size_t> indexMap_;
+};
+
+#endif
+
+}} // end namespace Dune::experimental
diff --git a/dune/vtk/datacollectors/structureddatacollector.hh b/dune/vtk/datacollectors/structureddatacollector.hh
index 80b54dd..667696f 100644
--- a/dune/vtk/datacollectors/structureddatacollector.hh
+++ b/dune/vtk/datacollectors/structureddatacollector.hh
@@ -1,156 +1,78 @@
#pragma once
-#include <dune/vtk/datacollector.hh>
+#include <dune/vtk/datacollectors/continuousdatacollector.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>>
+template <class GridView, class Derived>
+class StructuredDataCollectorInterface
+ : public DataCollectorInterface<GridView, Derived>
{
- enum { dim = GridView::dimension };
-
- using Self = StructuredDataCollector;
- using Super = DataCollectorInterface<GridView, Self>;
+protected:
+ using Super = DataCollectorInterface<GridView, Derived>;
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)
+ StructuredDataCollectorInterface (GridView const& gridView)
: Super(gridView)
+ , defaultDataCollector_(gridView)
{}
/// Return number of grid vertices
std::uint64_t numPointsImpl () const
{
- return gridView_.size(dim);
+ return gridView_.size(GridView::dimension);
}
- std::array<int, 6> const& extent () const
+ void updateImpl ()
{
- return extent_;
+ defaultDataCollector_.update();
}
- auto const& origin () const
+ std::array<int, 6> const& wholeExtent () const
{
- return origin_;
+ return this->asDerived().wholeExtentImpl();
}
- auto const& spacing () const
+ FieldVector<ctype, 3> const& origin () const
{
- return spacing_;
+ return this->asDerived().originImpl();
}
- void updateImpl ()
+ FieldVector<ctype, 3> const& spacing () const
{
- // 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;
- }
+ return this->asDerived().spacingImpl();
+ }
- origin_[i] = 0;
- }
+ template <class Writer>
+ void writeLocalPiece (Writer const& writer) const
+ {
+ this->asDerived().writeLocalPieceImpl(writer);
+ }
+
+ template <class Writer>
+ void writePieces (Writer const& writer) const
+ {
+ this->asDerived().writePiecesImpl(writer);
}
/// 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;
+ return defaultDataCollector_.template points<T>();
}
/// 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;
+ return defaultDataCollector_.template pointData<T>(fct);
}
private:
- std::array<int, 6> extent_;
- FieldVector<ctype,3> spacing_;
- FieldVector<ctype,3> origin_;
- std::vector<std::size_t> indexMap_;
+ DefaultDataCollector<GridView> defaultDataCollector_;
};
}} // end namespace Dune::experimental
diff --git a/dune/vtk/datacollectors/yaspdatacollector.hh b/dune/vtk/datacollectors/yaspdatacollector.hh
new file mode 100644
index 0000000..e557941
--- /dev/null
+++ b/dune/vtk/datacollectors/yaspdatacollector.hh
@@ -0,0 +1,124 @@
+#pragma once
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/vtk/datacollectors/structureddatacollector.hh>
+
+namespace Dune { namespace experimental
+{
+// Specialization for YaspGrid
+template <class GridView>
+class YaspDataCollector
+ : public StructuredDataCollectorInterface<GridView, YaspDataCollector<GridView>>
+{
+ enum { dim = GridView::dimension };
+
+ using Self = YaspDataCollector;
+ using Super = StructuredDataCollectorInterface<GridView, Self>;
+ using Super::gridView_;
+ using ctype = typename GridView::ctype;
+
+public:
+ YaspDataCollector (GridView const& gridView)
+ : Super(gridView)
+ , wholeExtent_(filledArray<6,int>(0))
+ , origin_(0.0)
+ , spacing_(0.0)
+ , level_(gridView.template begin<0,Interior_Partition>()->level())
+ {}
+
+ std::array<int, 6> const& wholeExtentImpl () const
+ {
+ return wholeExtent_;
+ }
+
+ auto const& originImpl () const
+ {
+ return origin_;
+ }
+
+ auto const& spacingImpl () const
+ {
+ return spacing_;
+ }
+
+ void updateImpl ()
+ {
+ Super::updateImpl();
+
+ for (int i = 0; i < dim; ++i) {
+ wholeExtent_[2*i] = 0;
+ wholeExtent_[2*i+1] = gridView_.grid().levelSize(level_,i);
+ }
+
+ auto it = gridView_.grid().begin(level_);
+ initGeometry(it->coords);
+ }
+
+ template <class Coords>
+ void initGeometry (Coords const& coords) { DUNE_THROW(NotImplemented, "Coordinate-Type not implemented!"); }
+
+ void initGeometry (EquidistantCoordinates<ctype,dim> const& coords)
+ {
+ for (int i = 0; i < dim; ++i) {
+ spacing_[i] = coords.meshsize(i,0);
+ origin_[i] = 0;
+ }
+ }
+
+ void initGeometry (EquidistantOffsetCoordinates<ctype,dim> const& coords)
+ {
+ for (int i = 0; i < dim; ++i) {
+ spacing_[i] = coords.meshsize(i,0);
+ origin_[i] = coords.origin(i);
+ }
+ }
+
+ void initGeometry (TensorProductCoordinates<ctype,dim> const& coords)
+ {
+ for (int i = 0; i < dim; ++i) {
+ spacing_[i] = coords.meshsize(i,0); // is not constant, but also not used.
+ origin_[i] = coords.coordinate(i,0);
+ }
+ }
+
+ template <class Writer>
+ void writeLocalPieceImpl (Writer const& writer) const
+ {
+ auto const& gl = *gridView_.grid().begin(level_);
+ auto const& g = gl.interior[0];
+ auto extent = filledArray<6,int>(0);
+ auto const& gc = *g.dataBegin();
+ for (int i = 0; i < dim; ++i) {
+ extent[2*i] = gc.min(i);
+ extent[2*i+1] = gc.max(i)+1;
+ }
+ writer(extent);
+ }
+
+ template <class Writer>
+ void writePiecesImpl (Writer const& writer) const
+ {
+ auto extent = filledArray<6,int>(0);
+ // for (auto const& part : gridView_.gridLevel().template partition<InteriorEntity>())
+ // {
+ // for (int i = 0; i < dim; ++i) {
+ // extent[2*i] = part.begin()[i];
+ // extent[2*i+1] = part.end()[i]-1;
+ // }
+
+ int num_ranks = -1;
+ MPI_Comm_size(MPI_COMM_WORLD, &num_ranks);
+ for (int p = 0; p < num_ranks; ++p) {
+ writer(p, extent, false);
+ }
+ }
+
+private:
+ std::array<int, 6> wholeExtent_;
+ FieldVector<ctype,3> spacing_;
+ FieldVector<ctype,3> origin_;
+ std::vector<std::size_t> indexMap_;
+ int level_;
+};
+
+}} // end namespace Dune::experimental
diff --git a/dune/vtk/vtkimagedatawriter.hh b/dune/vtk/vtkimagedatawriter.hh
index b1e4345..55a1df5 100644
--- a/dune/vtk/vtkimagedatawriter.hh
+++ b/dune/vtk/vtkimagedatawriter.hh
@@ -14,7 +14,7 @@
namespace Dune { namespace experimental
{
/// File-Writer for VTK .vtu files
- template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
+ template <class GridView, class DataCollector>
class VtkImageDataWriter
: public VtkWriter<GridView, DataCollector>
{
diff --git a/dune/vtk/vtkimagedatawriter.impl.hh b/dune/vtk/vtkimagedatawriter.impl.hh
index 0f0ac77..b392047 100644
--- a/dune/vtk/vtkimagedatawriter.impl.hh
+++ b/dune/vtk/vtkimagedatawriter.impl.hh
@@ -28,16 +28,14 @@ void VtkImageDataWriter<GV,DC>
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();
+ auto const& wholeExtent = dataCollector_.wholeExtent();
for (int i = 0; i < 3; ++i) {
- out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
+ out << (i == 0 ? "" : " ") << wholeExtent[2*i] << " " << wholeExtent[2*i+1];
}
out << "\" Origin=\"";
auto const& origin = dataCollector_.origin();
@@ -50,11 +48,15 @@ void VtkImageDataWriter<GV,DC>
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";
+
+ dataCollector_.writeLocalPiece([&out](std::array<int,6> const& extent)
+ {
+ 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; });
@@ -125,9 +127,9 @@ void VtkImageDataWriter<GV,DC>
<< "byte_order=\"" << this->getEndian() << "\" header_type=\"UInt64\""
<< (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\">\n" : ">\n");
out << "<PImageData GhostLevel=\"0\" WholeExtent=\"";
- auto const& extent = dataCollector_.extent();
+ auto const& wholeExtent = dataCollector_.wholeExtent();
for (int i = 0; i < 3; ++i) {
- out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
+ out << (i == 0 ? "" : " ") << wholeExtent[2*i] << " " << wholeExtent[2*i+1];
}
out << "\" Origin=\"";
auto const& origin = dataCollector_.origin();
@@ -168,15 +170,20 @@ void VtkImageDataWriter<GV,DC>
}
// 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";
+ dataCollector_.writePieces([&out,pfilename,ext=this->fileExtension()](int p, std::array<int,6> const& extent, bool write_extent)
+ {
+ out << "<Piece Source=\"" << pfilename << "_p" << std::to_string(p) << "." << ext << "\"";
+ if (write_extent) {
+ out << " Extent=\"";
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
+ }
+ out << "\"";
+ }
+ out << " />\n";
+ });
+
+ out << "</PImageData>\n";
out << "</VTKFile>";
}
diff --git a/dune/vtk/vtkstructuredgridwriter.hh b/dune/vtk/vtkstructuredgridwriter.hh
index 8c3aaa2..a3c3b68 100644
--- a/dune/vtk/vtkstructuredgridwriter.hh
+++ b/dune/vtk/vtkstructuredgridwriter.hh
@@ -9,12 +9,11 @@
#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>>
+ template <class GridView, class DataCollector>
class VtkStructuredGridWriter
: public VtkWriter<GridView, DataCollector>
{
diff --git a/dune/vtk/vtkstructuredgridwriter.impl.hh b/dune/vtk/vtkstructuredgridwriter.impl.hh
index 7900e33..d7f6294 100644
--- a/dune/vtk/vtkstructuredgridwriter.impl.hh
+++ b/dune/vtk/vtkstructuredgridwriter.impl.hh
@@ -28,26 +28,26 @@ void VtkStructuredGridWriter<GV,DC>
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=\"";
+ auto const& wholeExtent = dataCollector_.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 << (i == 0 ? "" : " ") << wholeExtent[2*i] << " " << wholeExtent[2*i+1];
}
out << "\">\n";
+ dataCollector_.writeLocalPiece([&out](std::array<int,6> const& extent)
+ {
+ 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; });
@@ -127,10 +127,10 @@ void VtkStructuredGridWriter<GV,DC>
<< "byte_order=\"" << this->getEndian() << "\" header_type=\"UInt64\""
<< (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\">\n" : ">\n");
out << "<PStructuredGrid GhostLevel=\"0\" WholeExtent=\"";
+ auto const& wholeExtent = dataCollector_.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 << (i == 0 ? "" : " ") << wholeExtent[2*i] << " " << wholeExtent[2*i+1];
+ }
out << "\">\n";
{ // Write data associated with grid points
@@ -165,16 +165,20 @@ void VtkStructuredGridWriter<GV,DC>
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";
- }
+ dataCollector_.writePieces([&out,pfilename,ext=this->fileExtension()](int p, std::array<int,6> const& extent, bool write_extent)
+ {
+ out << "<Piece Source=\"" << pfilename << "_p" << std::to_string(p) << "." << ext << "\"";
+ if (write_extent) {
+ out << " Extent=\"";
+ for (int i = 0; i < 3; ++i) {
+ out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
+ }
+ out << "\"";
+ }
+ out << " />\n";
+ });
- out << "</PUnstructuredGrid>\n";
+ out << "</PStructuredGrid>\n";
out << "</VTKFile>";
}
diff --git a/dune/vtk/vtkunstructuredgridwriter.impl.hh b/dune/vtk/vtkunstructuredgridwriter.impl.hh
index 4b619c3..226f004 100644
--- a/dune/vtk/vtkunstructuredgridwriter.impl.hh
+++ b/dune/vtk/vtkunstructuredgridwriter.impl.hh
@@ -28,8 +28,6 @@ void VtkUnstructuredGridWriter<GV,DC>
out << std::setprecision(std::numeric_limits<double>::digits10+2);
}
- dataCollector_.update();
-
std::vector<pos_type> offsets; // pos => offset
out << "<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" "
<< "byte_order=\"" << this->getEndian() << "\" header_type=\"UInt64\""
diff --git a/dune/vtk/vtkwriter.impl.hh b/dune/vtk/vtkwriter.impl.hh
index f70936a..184e477 100644
--- a/dune/vtk/vtkwriter.impl.hh
+++ b/dune/vtk/vtkwriter.impl.hh
@@ -34,6 +34,8 @@ void VtkWriter<GV,DC>
}
#endif
+ dataCollector_.update();
+
auto p = filesystem::path(fn);
auto name = p.stem();
p.remove_filename();
diff --git a/src/structuredgridwriter.cc b/src/structuredgridwriter.cc
index 8d08855..5fbcb7d 100644
--- a/src/structuredgridwriter.cc
+++ b/src/structuredgridwriter.cc
@@ -22,78 +22,106 @@
#include <dune/vtk/vtkstructuredgridwriter.hh>
#include <dune/vtk/vtkimagedatawriter.hh>
+#include <dune/vtk/datacollectors/yaspdatacollector.hh>
+#include <dune/vtk/datacollectors/spdatacollector.hh>
using namespace Dune;
using namespace Dune::experimental;
using namespace Dune::Functions;
-#define GRID_TYPE 1
-
-int main(int argc, char** argv)
+namespace Impl_
{
- Dune::MPIHelper::instance(argc, argv);
+ template <class GridView, class Grid>
+ struct StructuredDataCollector;
- const int dim = 3;
+#ifdef HAVE_DUNE_SPGRID
+ template<class GridView, class ct, int dim, template< int > class Ref, class Comm>
+ struct StructuredDataCollector<GridView, SPGrid<ct,dim,Ref,Comm>>
+ {
+ using type = SPDataCollector<GridView>;
+ };
+#endif
+
+ template<class GridView, int dim, class Coordinates>
+ struct StructuredDataCollector<GridView, YaspGrid<dim,Coordinates>>
+ {
+ using type = YaspDataCollector<GridView>;
+ };
+}
+
+template <class GridView>
+using StructuredDataCollector = typename Impl_::StructuredDataCollector<GridView, typename GridView::Grid>::type;
- 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();
+template <int dim>
+using int_ = std::integral_constant<int,dim>;
+template <class GridView>
+void write(std::string prefix, GridView const& gridView)
+{
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);
+ auto fct1 = makeDiscreteGlobalBasisFunction<double>(basis, p1function);
+ auto fct2 = makeAnalyticGridViewFunction([](auto const& x) {
+ return std::sin(10*x[0])*std::cos(10*x[1])+std::sin(10*x[2]);
+ }, gridView);
{
- using Writer = VtkStructuredGridWriter<GridView>;
+ using Writer = VtkStructuredGridWriter<GridView, StructuredDataCollector<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);
+ vtkWriter.addPointData(fct1, "p1");
+ vtkWriter.addCellData(fct1, "p0");
+ vtkWriter.addPointData(fct2, "analytic");
+ vtkWriter.write(prefix + "sg_ascii_float32.vts", Vtk::ASCII);
}
{
- using Writer = VtkImageDataWriter<GridView>;
+ using Writer = VtkImageDataWriter<GridView, StructuredDataCollector<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);
+ vtkWriter.addPointData(fct1, "p1");
+ vtkWriter.addCellData(fct1, "p0");
+ vtkWriter.addPointData(fct2, "analytic");
+ vtkWriter.write(prefix + "id_ascii_float32.vti", Vtk::ASCII);
}
}
+
+template <int dim>
+void write_yaspgrid(std::integral_constant<int,dim>)
+{
+ using GridType = YaspGrid<dim>;
+ FieldVector<double,dim> upperRight; upperRight = 1.0;
+ auto numElements = filledArray<dim,int>(8);
+ GridType grid(upperRight,numElements,0,0);
+
+ write("yasp_" + std::to_string(dim) + "d_", grid.leafGridView());
+}
+
+template <int dim>
+void write_spgrid(std::integral_constant<int,dim>)
+{
+#ifdef HAVE_DUNE_SPGRID
+ using GridType = SPGrid<double,dim, SPIsotropicRefinement>;
+ FieldVector<double,dim> upperRight; upperRight = 1.0;
+ auto numElements = filledArray<dim,int>(8);
+ GridType grid(SPDomain<double,dim>::unitCube(),numElements);
+
+ write("sp_" + std::to_string(dim) + "d_", grid.leafGridView());
+#endif
+}
+
+int main(int argc, char** argv)
+{
+ Dune::MPIHelper::instance(argc, argv);
+
+ Hybrid::forEach(std::make_tuple(int_<1>{}, int_<2>{}, int_<3>{}), [](auto const dim)
+ {
+ write_yaspgrid(dim);
+ write_spgrid(dim);
+ });
+}
--
GitLab