From f781632fb65cf5c88ae5a4558583f26842350b10 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Wed, 22 Aug 2018 17:42:50 +0200
Subject: [PATCH] added datacollector and quadratic vtk cell types
---
dune/vtk/datacollector.hh | 419 +++++++++++++++++++++++++++++++++
dune/vtk/filereader.hh | 5 +-
dune/vtk/filewriter.hh | 4 +-
dune/vtk/gridcreator.hh | 12 +-
dune/vtk/vtkfunction.hh | 38 ++-
dune/vtk/vtkreader.hh | 4 +-
dune/vtk/vtkreader.impl.hh | 6 +-
dune/vtk/vtktypes.cc | 146 ++++++++----
dune/vtk/vtktypes.hh | 42 +++-
dune/vtk/vtkwriter.hh | 41 ++--
dune/vtk/vtkwriter.impl.hh | 244 +++++++------------
src/CMakeLists.txt | 12 +
src/benchmark.cc | 114 +++++++++
src/datacollector.cc | 115 +++++++++
src/quadratic.cc | 84 +++++++
src/test/mixed_element_test.cc | 3 +-
src/test/reader_writer_test.cc | 1 +
src/vtkreader.cc | 1 +
src/vtkwriter.cc | 5 +-
19 files changed, 1015 insertions(+), 281 deletions(-)
create mode 100644 dune/vtk/datacollector.hh
create mode 100644 src/benchmark.cc
create mode 100644 src/datacollector.cc
create mode 100644 src/quadratic.cc
diff --git a/dune/vtk/datacollector.hh b/dune/vtk/datacollector.hh
new file mode 100644
index 0000000..328a86a
--- /dev/null
+++ b/dune/vtk/datacollector.hh
@@ -0,0 +1,419 @@
+#pragma once
+
+#include <algorithm>
+#include <cstdint>
+#include <vector>
+
+#include "vtktypes.hh"
+
+namespace Dune { namespace experimental {
+
+struct Cells
+{
+ std::vector<std::uint8_t> types;
+ std::vector<std::int64_t> offsets;
+ std::vector<std::int64_t> connectivity;
+};
+
+class DataCollectorInterface
+{
+private:
+ DataCollectorInterface () = default;
+
+public:
+ /// \brief Return the number of points in the grid
+ std::uint64_t numPoints () const;
+
+ /// \brief Return the number of cells in the grid
+ std::uint64_t numCells () const;
+
+ /// Update the DataCollector on the current GridView
+ void update ();
+
+ /// \brief Return a flat vector of point coordinates
+ /**
+ * All coordinates are extended to 3 components and concatenated.
+ * [p0_x, p0_y, p0_z, p1_x, p1_y, p1_z, ...]
+ * If the GridView::dimensionworld < 3, the remaining components are
+ * set to 0
+ **/
+ template <class T>
+ std::vector<T> points () const;
+
+ /// \brief Return cell types, offsets, and connectivity. \see Cells
+ Cells cells () const;
+
+ /// \brief Return a flat vector of function values evaluated at the grid points.
+ /**
+ * In case of a vector valued function, flat the vector entries:
+ * [fct(p0)_0, fct(p0)_1, fct(p0)_2, fct(p1)_0, ...]
+ * where the vector dimension must be 3 (possible extended by 0s)
+ * In case of tensor valued function, flat the tensor row-wise:
+ * [fct(p0)_00, fct(p0)_01, fct(p0)_02, fct(p0)_10, fct(p0)_11, fct(p0)_12, fct(p0)_20...]
+ * where the tensor dimension must be 3x3 (possible extended by 0s)
+ **/
+ template <class T, class GlobalFunction>
+ std::vector<T> pointData (GlobalFunction const& fct) const;
+
+ /// \brief Return a flat vector of function values evaluated at the grid cells. \see pointData.
+ template <class T, class GlobalFunction>
+ std::vector<T> cellData (GlobalFunction const& fct) const;
+};
+
+
+/// Implementation of \ref DataCollector for linear cells, with continuous data.
+template <class GridView>
+class DefaultDataCollector
+{
+ enum { dim = GridView::dimension };
+
+public:
+ DefaultDataCollector (GridView const& gridView)
+ : gridView_(gridView)
+ {}
+
+ void update () {}
+
+ std::uint64_t numPoints () const
+ {
+ return gridView_.size(dim);
+ }
+
+ std::uint64_t numCells () const
+ {
+ return gridView_.size(0);
+ }
+
+ template <class T>
+ std::vector<T> points () const
+ {
+ std::vector<T> data(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;
+ }
+
+ Cells cells () const
+ {
+ auto const& indexSet = gridView_.indexSet();
+ auto types = indexSet.types(0);
+ int maxVertices = std::accumulate(types.begin(), types.end(), 1, [](int m, GeometryType t) {
+ auto refElem = referenceElement<double,dim>(t);
+ return std::max(m, refElem.size(dim));
+ });
+
+ Cells cells;
+ cells.connectivity.reserve(numCells() * maxVertices);
+ cells.offsets.reserve(numCells());
+ cells.types.reserve(numCells());
+
+ std::int64_t old_o = 0;
+ for (auto const& c : elements(gridView_)) {
+ 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)) );
+ cells.offsets.push_back(old_o += c.subEntities(dim));
+ cells.types.push_back(cellType.type());
+ }
+
+ return cells;
+ }
+
+ template <class T, class GlobalFunction>
+ std::vector<T> pointData (GlobalFunction const& fct) const
+ {
+ std::vector<T> data(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) {
+ std::size_t idx = fct.ncomps() * indexSet.subIndex(e,cellType.permutation(j),dim);
+ for (int comp = 0; comp < fct.ncomps(); ++comp)
+ data[idx + comp] = T(localFct.evaluate(comp, refElem.position(cellType.permutation(j),dim)));
+ }
+ localFct.unbind();
+ }
+ return data;
+ }
+
+ template <class T, class GlobalFunction>
+ std::vector<T> cellData (GlobalFunction const& fct) const
+ {
+ std::vector<T> data(numCells() * fct.ncomps());
+ auto const& indexSet = gridView_.indexSet();
+ auto localFct = localFunction(fct);
+ for (auto const& e : elements(gridView_)) {
+ localFct.bind(e);
+ auto refElem = referenceElement(e.geometry());
+ std::size_t idx = fct.ncomps() * indexSet.index(e);
+ for (int comp = 0; comp < fct.ncomps(); ++comp)
+ data[idx + comp] = T(localFct.evaluate(comp, refElem.position(0,0)));
+ localFct.unbind();
+ }
+ return data;
+ }
+
+private:
+ GridView gridView_;
+};
+
+
+/// Implementation of \ref DataCollector for linear cells, with discontinuous data.
+template <class GridView>
+class DiscontinuousDataCollector
+{
+ enum { dim = GridView::dimension };
+
+public:
+ DiscontinuousDataCollector (GridView const& gridView)
+ : gridView_(gridView)
+ {}
+
+ void update ()
+ {
+ numPoints_ = 0;
+ indexMap_.resize(gridView_.size(dim));
+ std::int64_t vertex_idx = 0;
+ auto const& indexSet = gridView_.indexSet();
+ for (auto const& c : elements(gridView_)) {
+ numPoints_ += c.subEntities(dim);
+ for (int i = 0; i < c.subEntities(dim); ++i)
+ indexMap_[indexSet.subIndex(c, i, dim)] = vertex_idx++;
+ }
+ }
+
+ std::uint64_t numPoints () const
+ {
+ return numPoints_;
+ }
+
+ std::uint64_t numCells () const
+ {
+ return gridView_.size(0);
+ }
+
+ template <class T>
+ std::vector<T> points () const
+ {
+ std::vector<T> data(numPoints() * 3);
+ auto const& indexSet = gridView_.indexSet();
+ for (auto const& element : elements(gridView_)) {
+ 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);
+ 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;
+ }
+
+ Cells cells () const
+ {
+ Cells cells;
+ cells.connectivity.reserve(numPoints());
+ cells.offsets.reserve(numCells());
+ cells.types.reserve(numCells());
+
+ std::int64_t old_o = 0;
+ auto const& indexSet = gridView_.indexSet();
+ for (auto const& c : elements(gridView_)) {
+ 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)];
+ cells.connectivity.push_back(vertex_idx);
+ }
+ cells.offsets.push_back(old_o += c.subEntities(dim));
+ cells.types.push_back(cellType.type());
+ }
+
+ return cells;
+ }
+
+ template <class T, class GlobalFunction>
+ std::vector<T> pointData (GlobalFunction const& fct) const
+ {
+ std::vector<T> data(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) {
+ std::size_t idx = fct.ncomps() * indexMap_[indexSet.subIndex(e, cellType.permutation(j), dim)];
+ for (int comp = 0; comp < fct.ncomps(); ++comp)
+ data[idx + comp] = T(localFct.evaluate(comp, refElem.position(cellType.permutation(j),dim)));
+ }
+ localFct.unbind();
+ }
+ return data;
+ }
+
+ template <class T, class GlobalFunction>
+ std::vector<T> cellData (GlobalFunction const& fct) const
+ {
+ std::vector<T> data(numCells() * fct.ncomps());
+ auto const& indexSet = gridView_.indexSet();
+ auto localFct = localFunction(fct);
+ for (auto const& e : elements(gridView_)) {
+ localFct.bind(e);
+ auto refElem = referenceElement(e.geometry());
+ std::size_t idx = fct.ncomps() * indexSet.index(e);
+ for (int comp = 0; comp < fct.ncomps(); ++comp)
+ data[idx + comp] = T(localFct.evaluate(comp, refElem.position(0,0)));
+ localFct.unbind();
+ }
+ return data;
+ }
+
+private:
+ GridView gridView_;
+ std::uint64_t numPoints_ = 0;
+ std::vector<std::int64_t> indexMap_;
+};
+
+
+/// Implementation of \ref DataCollector for quadratic cells, with continuous data.
+template <class GridView>
+class QuadraticDataCollector
+{
+ enum { dim = GridView::dimension };
+
+public:
+ QuadraticDataCollector (GridView const& gridView)
+ : gridView_(gridView)
+ {}
+
+ void update () {}
+
+ std::uint64_t numPoints () const
+ {
+ return gridView_.size(dim) + gridView_.size(dim-1);
+ }
+
+ std::uint64_t numCells () const
+ {
+ return gridView_.size(0);
+ }
+
+ template <class T>
+ std::vector<T> points () const
+ {
+ std::vector<T> data(numPoints() * 3);
+ auto const& indexSet = gridView_.indexSet();
+ for (auto const& element : elements(gridView_)) {
+ auto geometry = element.geometry();
+ auto refElem = referenceElement<T,dim>(element.type());
+
+ // vertices
+ for (int i = 0; i < element.subEntities(dim); ++i) {
+ std::size_t idx = 3 * indexSet.subIndex(element, i, dim);
+ auto v = geometry.global(refElem.position(i,dim));
+ 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);
+ }
+ // edge centers
+ for (int i = 0; i < element.subEntities(dim-1); ++i) {
+ std::size_t idx = 3 * (indexSet.subIndex(element, i, dim-1) + gridView_.size(dim));
+ auto v = geometry.global(refElem.position(i,dim-1));
+ 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;
+ }
+
+ Cells cells () const
+ {
+ Cells cells;
+ cells.connectivity.reserve(numPoints());
+ cells.offsets.reserve(numCells());
+ cells.types.reserve(numCells());
+
+ std::int64_t old_o = 0;
+ auto const& indexSet = gridView_.indexSet();
+ for (auto const& c : elements(gridView_)) {
+ Vtk::CellType cellType(c.type(), Vtk::QUADRATIC);
+ for (int j = 0; j < c.subEntities(dim); ++j) {
+ int k = cellType.permutation(j);
+ std::int64_t point_idx = indexSet.subIndex(c,k,dim);
+ cells.connectivity.push_back(point_idx);
+ }
+ for (int j = 0; j < c.subEntities(dim-1); ++j) {
+ int k = cellType.permutation(c.subEntities(dim) + j);
+ std::int64_t point_idx = (indexSet.subIndex(c,k,dim-1) + gridView_.size(dim));
+ cells.connectivity.push_back(point_idx);
+ }
+ cells.offsets.push_back(old_o += c.subEntities(dim)+c.subEntities(dim-1));
+ cells.types.push_back(cellType.type());
+ }
+
+ return cells;
+ }
+
+ template <class T, class GlobalFunction>
+ std::vector<T> pointData (GlobalFunction const& fct) const
+ {
+ std::vector<T> data(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(), Vtk::QUADRATIC};
+ 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)));
+ }
+ for (int j = 0; j < e.subEntities(dim-1); ++j) {
+ int k = cellType.permutation(e.subEntities(dim) + j);
+ std::size_t idx = fct.ncomps() * (indexSet.subIndex(e, k, dim-1) + gridView_.size(dim));
+ for (int comp = 0; comp < fct.ncomps(); ++comp)
+ data[idx + comp] = T(localFct.evaluate(comp, refElem.position(k, dim-1)));
+ }
+ localFct.unbind();
+ }
+ return data;
+ }
+
+ template <class T, class GlobalFunction>
+ std::vector<T> cellData (GlobalFunction const& fct) const
+ {
+ std::vector<T> data(numCells() * fct.ncomps());
+ auto const& indexSet = gridView_.indexSet();
+ auto localFct = localFunction(fct);
+ for (auto const& e : elements(gridView_)) {
+ localFct.bind(e);
+ auto refElem = referenceElement(e.geometry());
+ std::size_t idx = fct.ncomps() * indexSet.index(e);
+ for (int comp = 0; comp < fct.ncomps(); ++comp)
+ data[idx + comp] = T(localFct.evaluate(comp, refElem.position(0,0)));
+ localFct.unbind();
+ }
+ return data;
+ }
+
+private:
+ GridView gridView_;
+};
+
+}} // end namespace Dune::experimental
diff --git a/dune/vtk/filereader.hh b/dune/vtk/filereader.hh
index 5e2c313..aad0764 100644
--- a/dune/vtk/filereader.hh
+++ b/dune/vtk/filereader.hh
@@ -4,7 +4,7 @@
#include <string>
#include <utility>
-namespace Dune
+namespace Dune { namespace experimental
{
template <class Grid, class FilerReaderImp>
class FileReader
@@ -68,5 +68,4 @@ namespace Dune
}
};
-
-} // end namespace Dune
+}} // end namespace Dune::experimental
diff --git a/dune/vtk/filewriter.hh b/dune/vtk/filewriter.hh
index ffe21e0..f184f34 100644
--- a/dune/vtk/filewriter.hh
+++ b/dune/vtk/filewriter.hh
@@ -2,7 +2,7 @@
#include <string>
-namespace Dune
+namespace Dune { namespace experimental
{
class FileWriter
{
@@ -14,4 +14,4 @@ namespace Dune
virtual void write (std::string const& filename) = 0;
};
-} // end namespace Dune
+}} // end namespace Dune::experimental
diff --git a/dune/vtk/gridcreator.hh b/dune/vtk/gridcreator.hh
index 6b969a7..7399bdc 100644
--- a/dune/vtk/gridcreator.hh
+++ b/dune/vtk/gridcreator.hh
@@ -10,7 +10,7 @@
#include "vtktypes.hh"
-namespace Dune
+namespace Dune { namespace experimental
{
// Create a grid where the input points and connectivity is already
// connected correctly.
@@ -28,9 +28,9 @@ namespace Dune
std::size_t idx = 0;
for (std::size_t i = 0; i < types.size(); ++i) {
- if (Vtk::Map::type.count(types[i]) == 0)
+ if (Vtk::Map::from_type.count(types[i]) == 0)
DUNE_THROW(Exception, "Unknown ElementType: " << types[i]);
- auto type = Vtk::Map::type[types[i]];
+ auto type = Vtk::Map::from_type[types[i]];
Vtk::CellType cellType{type};
auto refElem = referenceElement<double,Grid::dimension>(type);
@@ -92,9 +92,9 @@ namespace Dune
idx = 0;
for (std::size_t i = 0; i < types.size(); ++i) {
- if (Vtk::Map::type.count(types[i]) == 0)
+ if (Vtk::Map::from_type.count(types[i]) == 0)
DUNE_THROW(Exception, "Unknown ElementType: " << types[i]);
- auto type = Vtk::Map::type[types[i]];
+ auto type = Vtk::Map::from_type[types[i]];
Vtk::CellType cellType{type};
std::size_t nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]);
@@ -120,4 +120,4 @@ namespace Dune
}
};
-} // end namespace Dune
+}} // end namespace Dune::experimental
diff --git a/dune/vtk/vtkfunction.hh b/dune/vtk/vtkfunction.hh
index 1c6820b..3220e0d 100644
--- a/dune/vtk/vtkfunction.hh
+++ b/dune/vtk/vtkfunction.hh
@@ -5,7 +5,7 @@
#include <dune/common/std/type_traits.hh>
#include <dune/functions/common/typeerasure.hh>
-namespace Dune
+namespace Dune { namespace experimental
{
/// An abstract base class for LocalFunctions
template <class GridView>
@@ -69,7 +69,7 @@ namespace Dune
double evaluateImpl (int comp, LocalCoordinate const& xi, std::true_type) const
{
auto y = this->get()(xi);
- return y[comp];
+ return comp < y.size() ? y[comp] : 0.0;
}
// Return the scalar values
@@ -128,7 +128,7 @@ namespace Dune
{
public:
/// Create a local function
- virtual VTKLocalFunction<GridView> makelocalFunction () const = 0;
+ virtual VTKLocalFunction<GridView> makeLocalFunction () const = 0;
/// Virtual destructor
virtual ~VTKFunctionInterface () = default;
@@ -143,12 +143,9 @@ namespace Dune
{
public:
using Wrapper::Wrapper;
- using Function = typename Wrapper::Wrapped;
- using Interface = VTKFunctionInterface<GridView>;
-
- virtual VTKLocalFunction<GridView> makelocalFunction () const override
+ virtual VTKLocalFunction<GridView> makeLocalFunction () const override
{
- return {localFunction(this->get())};
+ return VTKLocalFunction<GridView>{localFunction(this->get())};
}
};
};
@@ -163,16 +160,17 @@ namespace Dune
VTKFunctionImpl<GridView>::template Model>;
template <class F>
- using IsGridFunction = decltype(localFunction(std::declval<F>()));
+ using HasLocalFunction = decltype(localFunction(std::declval<F>()));
public:
- template <class F, disableCopyMove<VTKFunction, F> = 0>
+ template <class F>
VTKFunction (F&& f, std::string name, int ncomps = 1)
: Super(std::forward<F>(f))
, name_(std::move(name))
- , ncomps_(ncomps)
+ , ncomps_(ncomps > 1 ? 3 : 1)
{
- static_assert(Std::is_detected<IsGridFunction,F>::value,
+ assert(1 <= ncomps && ncomps <= 3);
+ static_assert(Std::is_detected<HasLocalFunction,F>::value,
"Requires A GridFunction to be passed to the VTKFunction.");
}
@@ -181,13 +179,7 @@ namespace Dune
/// Create a LocalFunction
friend VTKLocalFunction<GridView> localFunction (VTKFunction const& self)
{
- return self.asInterface().makelocalFunction();
- }
-
- /// Return the number of components of the Range
- int ncomps () const
- {
- return ncomps_; // TODO: detect automatically. Is this always possible?
+ return self.asInterface().makeLocalFunction();
}
/// Return a name associated with the function
@@ -196,10 +188,16 @@ namespace Dune
return name_;
}
+ /// Return the number of components of the Range
+ int ncomps () const
+ {
+ return ncomps_;
+ }
+
private:
std::string name_;
int ncomps_ = 1;
};
-} // end namespace Dune
+}} // end namespace Dune::experimental
diff --git a/dune/vtk/vtkreader.hh b/dune/vtk/vtkreader.hh
index 3b90c76..2dc6beb 100644
--- a/dune/vtk/vtkreader.hh
+++ b/dune/vtk/vtkreader.hh
@@ -7,7 +7,7 @@
#include "gridcreator.hh"
#include "vtktypes.hh"
-namespace Dune
+namespace Dune { namespace experimental
{
/// File-Reader for Vtk .vtu files
/**
@@ -129,6 +129,6 @@ namespace Dune
std::uint64_t offset0_;
};
-} // end namespace Dune
+}} // end namespace Dune::experimental
#include "vtkreader.impl.hh"
diff --git a/dune/vtk/vtkreader.impl.hh b/dune/vtk/vtkreader.impl.hh
index 819528c..0c735cb 100644
--- a/dune/vtk/vtkreader.impl.hh
+++ b/dune/vtk/vtkreader.impl.hh
@@ -12,7 +12,7 @@
#include "utility/filesystem.hh"
#include "utility/string.hh"
-namespace Dune {
+namespace Dune { namespace experimental {
template <class Grid, class Creator>
void VtkReader<Grid,Creator>::readFromFile (std::string const& filename)
@@ -91,7 +91,7 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename)
bool closed = false;
auto attr = parseXml(line, closed);
- data_type = Vtk::Map::datatype[attr["type"]];
+ data_type = Vtk::Map::to_datatype[attr["type"]];
if (!attr["Name"].empty())
data_name = to_lower(attr["Name"]);
@@ -512,4 +512,4 @@ std::map<std::string, std::string> VtkReader<Grid,Creator>::parseXml (std::strin
return attr;
}
-} // end namespace Dune
+}} // end namespace Dune::experimental
diff --git a/dune/vtk/vtktypes.cc b/dune/vtk/vtktypes.cc
index e0d08ec..feae227 100644
--- a/dune/vtk/vtktypes.cc
+++ b/dune/vtk/vtktypes.cc
@@ -2,21 +2,21 @@
#include <iostream>
-namespace Dune {
+namespace Dune { namespace experimental {
namespace Vtk {
-std::map<std::uint8_t, GeometryType> Map::type = {
- { 1, GeometryTypes::vertex },
- { 3, GeometryTypes::line },
- { 5, GeometryTypes::triangle },
- { 9, GeometryTypes::quadrilateral },
- {10, GeometryTypes::tetrahedron },
- {12, GeometryTypes::hexahedron },
- {13, GeometryTypes::prism },
- {14, GeometryTypes::pyramid },
+std::map<std::uint8_t, GeometryType> Map::from_type = {
+ {VERTEX, GeometryTypes::vertex },
+ {LINE, GeometryTypes::line },
+ {TRIANGLE, GeometryTypes::triangle },
+ {QUAD, GeometryTypes::quadrilateral },
+ {TETRA, GeometryTypes::tetrahedron },
+ {HEXAHEDRON, GeometryTypes::hexahedron },
+ {WEDGE, GeometryTypes::prism },
+ {PYRAMID, GeometryTypes::pyramid },
};
-std::map<std::string, DataTypes> Map::datatype = {
+std::map<std::string, DataTypes> Map::to_datatype = {
{"Int8", INT8},
{"UInt8", UINT8},
{"Int16", INT16},
@@ -29,49 +29,93 @@ std::map<std::string, DataTypes> Map::datatype = {
{"Float64", FLOAT64}
};
+std::map<DataTypes, std::string> Map::from_datatype = {
+ {INT8, "Int8"},
+ {UINT8, "UInt8"},
+ {INT16, "Int16"},
+ {UINT16, "UInt16"},
+ {INT32, "Int32"},
+ {UINT32, "UInt32"},
+ {INT64, "Int64"},
+ {UINT64, "UInt64"},
+ {FLOAT32, "Float32"},
+ {FLOAT64, "Float64"}
+};
-CellType::CellType (GeometryType const& t)
+CellType::CellType (GeometryType const& t, CellParametrization parametrization)
: noPermutation_(true)
{
- if (t.isVertex()) {
- type_ = 1;
- permutation_ = {0};
- }
- else if (t.isLine()) {
- type_ = 3;
- permutation_ = {0,1};
- }
- else if (t.isTriangle()) {
- type_ = 5;
- permutation_ = {0,1,2};
- }
- else if (t.isQuadrilateral()) {
- type_ = 9;
- permutation_ = {0,1,2,3};
- }
- else if (t.isTetrahedron()) {
- type_ = 10;
- permutation_ = {0,1,2,3};
- }
- else if (t.isHexahedron()) {
- type_ = 12;
- permutation_ = {0,1,3,2,4,5,7,6};
- noPermutation_ = false;
- }
- else if (t.isPrism()) {
- type_ = 13;
- permutation_ = {0,2,1,3,5,4};
- noPermutation_ = false;
- }
- else if (t.isPyramid()) {
- type_ = 14;
- permutation_ = {0,1,3,2,4};
- noPermutation_ = false;
- }
- else {
- std::cerr << "Geometry Type not supported by VTK!\n";
- std::abort();
+ if (parametrization == LINEAR) {
+ if (t.isVertex()) {
+ type_ = VERTEX;
+ permutation_ = {0};
+ }
+ else if (t.isLine()) {
+ type_ = LINE;
+ permutation_ = {0,1};
+ }
+ else if (t.isTriangle()) {
+ type_ = TRIANGLE;
+ permutation_ = {0,1,2};
+ }
+ else if (t.isQuadrilateral()) {
+ type_ = QUAD;
+ permutation_ = {0,1,3,2};
+ noPermutation_ = false;
+ }
+ else if (t.isTetrahedron()) {
+ type_ = TETRA;
+ permutation_ = {0,1,2,3};
+ }
+ else if (t.isHexahedron()) {
+ type_ = HEXAHEDRON;
+ permutation_ = {0,1,3,2,4,5,7,6};
+ noPermutation_ = false;
+ }
+ else if (t.isPrism()) {
+ type_ = WEDGE;
+ permutation_ = {0,2,1,3,5,4};
+ noPermutation_ = false;
+ }
+ else if (t.isPyramid()) {
+ type_ = PYRAMID;
+ permutation_ = {0,1,3,2,4};
+ noPermutation_ = false;
+ }
+ else {
+ std::cerr << "Geometry Type not supported by VTK!\n";
+ std::abort();
+ }
+ } else if (parametrization == QUADRATIC) {
+ if (t.isLine()) {
+ type_ = QUADRATIC_EDGE;
+ permutation_ = {0,1, 0};
+ }
+ else if (t.isTriangle()) {
+ type_ = QUADRATIC_TRIANGLE;
+ permutation_ = {0,1,2, 0,2,1};
+ noPermutation_ = false;
+ }
+ else if (t.isQuadrilateral()) {
+ type_ = QUADRATIC_QUAD;
+ permutation_ = {0,1,3,2, 3,1,0,2}; // maybe need inverse permutation as well
+ noPermutation_ = false;
+ }
+ else if (t.isTetrahedron()) {
+ type_ = QUADRATIC_TETRA;
+ permutation_ = {0,1,2,3, 0,2,1,3,4,5};
+ noPermutation_ = false;
+ }
+ else if (t.isHexahedron()) {
+ type_ = QUADRATIC_HEXAHEDRON;
+ permutation_ = {0,1,3,2,4,5,7,6, 6,5,7,4,10,9,11,8,0,1,3,2}; // maybe need inverse permutation as well
+ noPermutation_ = false;
+ }
+ else {
+ std::cerr << "Geometry Type not supported by VTK!\n";
+ std::abort();
+ }
}
}
-}} // end namespace Vtk::Dune
+}}} // end namespace Dune::experimental::Vtk
diff --git a/dune/vtk/vtktypes.hh b/dune/vtk/vtktypes.hh
index 8376e7e..68342ce 100644
--- a/dune/vtk/vtktypes.hh
+++ b/dune/vtk/vtktypes.hh
@@ -7,7 +7,7 @@
#include <dune/geometry/type.hh>
-namespace Dune
+namespace Dune { namespace experimental
{
namespace Vtk
{
@@ -28,20 +28,40 @@ namespace Dune
FLOAT64 = 64
};
- enum PositionTypes {
- VERTEX_DATA,
- CELL_DATA
+ enum CellParametrization {
+ LINEAR,
+ QUADRATIC
};
- enum ContinuityTypes {
- CONFORMING,
- NONCONFORMING
+ enum CellTypes : std::uint8_t {
+ // Linear VTK cell types
+ VERTEX = 1,
+ POLY_VERTEX = 2, // not supported
+ LINE = 3,
+ POLY_LINE = 4, // not supported
+ TRIANGLE = 5,
+ TRIANGLE_STRIP = 6, // not supported
+ POLYGON = 7, // not supported
+ PIXEL = 8, // not supported
+ QUAD = 9,
+ TETRA = 10,
+ VOXEL = 11, // not supported
+ HEXAHEDRON = 12,
+ WEDGE = 13,
+ PYRAMID = 14,
+ // Quadratic VTK cell types
+ QUADRATIC_EDGE = 21,
+ QUADRATIC_TRIANGLE = 22,
+ QUADRATIC_QUAD = 23,
+ QUADRATIC_TETRA = 24,
+ QUADRATIC_HEXAHEDRON = 25
};
struct Map
{
- static std::map<std::uint8_t, GeometryType> type; // VTK Cell type -> Dune::GeometryType
- static std::map<std::string, DataTypes> datatype; // String -> DataTypes
+ static std::map<std::uint8_t, GeometryType> from_type; // VTK Cell type -> Dune::GeometryType
+ static std::map<std::string, DataTypes> to_datatype; // String -> DataTypes
+ static std::map<DataTypes, std::string> from_datatype; // DataTypes -> String
};
@@ -49,7 +69,7 @@ namespace Dune
class CellType
{
public:
- CellType (GeometryType const& t);
+ CellType (GeometryType const& t, CellParametrization = LINEAR);
/// Return VTK Cell type
std::uint8_t type () const
@@ -75,4 +95,4 @@ namespace Dune
};
} // end namespace Vtk
-} // end namespace Dune
+}} // end namespace Dune::experimental
diff --git a/dune/vtk/vtkwriter.hh b/dune/vtk/vtkwriter.hh
index 05a809b..d3dcca1 100644
--- a/dune/vtk/vtkwriter.hh
+++ b/dune/vtk/vtkwriter.hh
@@ -4,14 +4,15 @@
#include <iosfwd>
#include <map>
+#include "datacollector.hh"
#include "filewriter.hh"
#include "vtkfunction.hh"
#include "vtktypes.hh"
-namespace Dune
+namespace Dune { namespace experimental
{
/// File-Writer for Vtk .vtu files
- template <class GridView>
+ template <class GridView, class DataCollector = DefaultDataCollector<GridView>>
class VtkWriter
: public FileWriter
{
@@ -21,10 +22,15 @@ namespace Dune
using LocalFunction = VTKLocalFunction<GridView>;
using pos_type = typename std::ostream::pos_type;
+ enum PositionTypes {
+ POINT_DATA,
+ CELL_DATA
+ };
+
public:
/// Constructor, stores the gridView
VtkWriter (GridView const& gridView)
- : gridView_(gridView)
+ : dataCollector_(gridView)
{}
/// Write the attached data to the file
@@ -38,13 +44,13 @@ namespace Dune
Vtk::FormatTypes format,
Vtk::DataTypes datatype = Vtk::FLOAT32);
- /// Attach vertex data to the writer
+ /// Attach point data to the writer
template <class GridViewFunction>
- VtkWriter& addVertexData (GridViewFunction const& gridViewFct,
+ VtkWriter& addPointData (GridViewFunction const& gridViewFct,
std::string const& name = {},
int ncomps = 1)
{
- vertexData_.emplace_back(gridViewFct, name, ncomps);
+ pointData_.emplace_back(gridViewFct, name, ncomps);
return *this;
}
@@ -59,17 +65,17 @@ namespace Dune
}
private:
- // Write \ref vertexData_ and \ref cellData_ with set \ref format_ and
+ // Write \ref pointData_ and \ref cellData_ with set \ref format_ and
// \ref datatype_ to file given by filename
void writeImpl (std::string const& filename) const;
- // Write the vertex or cell values given by the grid function `fct` to the
+ // Write the point or cell values given by the grid function `fct` to the
// output stream `out`. In case of binary format, stores the streampos of XML
// attributes "offset" in the vector `offsets`.
void writeData (std::ofstream& out,
std::vector<pos_type>& offsets,
GlobalFunction const& fct,
- Vtk::PositionTypes type) const;
+ PositionTypes type) const;
// Write the coordinates of the vertices to the output stream `out`. In case
// of binary format, stores the streampos of XML attributes "offset" in the
@@ -89,14 +95,14 @@ namespace Dune
std::uint64_t writeAppended (std::ofstream& out,
std::vector<T> const& values) const;
- // Collect vertex or cell data (depending on \ref PositionTypes) and pass
+ // Collect point or cell data (depending on \ref PositionTypes) and pass
// the resulting vector to \ref writeAppended.
template <class T>
std::uint64_t writeDataAppended (std::ofstream& out,
GlobalFunction const& localFct,
- Vtk::PositionTypes type) const;
+ PositionTypes type) const;
- // Collect vertex positions and pass the resulting vector to \ref writeAppended.
+ // Collect point positions and pass the resulting vector to \ref writeAppended.
template <class T>
std::uint64_t writePointsAppended (std::ofstream& out) const;
@@ -111,26 +117,21 @@ namespace Dune
return (reinterpret_cast<char*>(&i)[1] == 1 ? "BigEndian" : "LittleEndian");
}
- Vtk::CellType getType (GeometryType const& t) const
- {
- return {t};
- }
-
private:
- GridView gridView_;
+ mutable DataCollector dataCollector_;
std::string filename_;
Vtk::FormatTypes format_;
Vtk::DataTypes datatype_;
// attached data
- std::vector<GlobalFunction> vertexData_;
+ std::vector<GlobalFunction> pointData_;
std::vector<GlobalFunction> cellData_;
std::size_t const block_size = 1024*32;
int compression_level = -1; // in [0,9], -1 ... use default value
};
-} // end namespace Dune
+}} // end namespace Dune::experimental
#include "vtkwriter.impl.hh"
diff --git a/dune/vtk/vtkwriter.impl.hh b/dune/vtk/vtkwriter.impl.hh
index 942c2bc..e445274 100644
--- a/dune/vtk/vtkwriter.impl.hh
+++ b/dune/vtk/vtkwriter.impl.hh
@@ -18,10 +18,10 @@
#include "utility/filesystem.hh"
#include "utility/string.hh"
-namespace Dune {
+namespace Dune { namespace experimental {
-template <class GridView>
-void VtkWriter<GridView>::write (std::string const& fn, Vtk::FormatTypes format, Vtk::DataTypes datatype)
+template <class GV, class DC>
+void VtkWriter<GV,DC>::write (std::string const& fn, Vtk::FormatTypes format, Vtk::DataTypes datatype)
{
format_ = format;
datatype_ = datatype;
@@ -52,44 +52,59 @@ void VtkWriter<GridView>::write (std::string const& fn, Vtk::FormatTypes format,
}
-template <class GridView>
-void VtkWriter<GridView>::writeImpl (std::string const& filename) const
+template <class GV, class DC>
+void VtkWriter<GV,DC>::writeImpl (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+1);
+ out << std::setprecision(std::numeric_limits<float>::digits10+2);
else
- out << std::setprecision(std::numeric_limits<double>::digits10+1);
+ 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=\"" << getEndian() << "\" header_type=\"UInt64\""
<< (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\">\n" : ">\n");
out << "<UnstructuredGrid>\n";
- out << "<Piece NumberOfPoints=\"" << gridView_.size(dimension) << "\" "
- << "NumberOfCells=\"" << gridView_.size(0) << "\">\n";
-
- // TODO: detect vector and scalars by ncomps() property of data
- out << "<PointData" << (!vertexData_.empty() ? " Scalars=\"" + vertexData_.front().name() + "\"" : "") << ">\n";
- for (auto const& v : vertexData_)
- writeData(out, offsets, v, Vtk::VERTEX_DATA);
- out << "</PointData>\n";
+ out << "<Piece NumberOfPoints=\"" << dataCollector_.numPoints() << "\" "
+ << "NumberOfCells=\"" << dataCollector_.numCells() << "\">\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_)
+ writeData(out, offsets, v, POINT_DATA);
+ out << "</PointData>\n";
+ }
- // TODO: detect vector and scalars by ncomps() property of data
- out << "<CellData" << (!cellData_.empty() ? " Scalars=\"" + cellData_.front().name() + "\"" : "") << ">\n";
- for (auto const& v : cellData_)
- writeData(out, offsets, v, Vtk::CELL_DATA);
- out << "</CellData>\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_)
+ writeData(out, offsets, v, CELL_DATA);
+ out << "</CellData>\n";
+ }
+ // Write point coordinates
out << "<Points>\n";
writePoints(out, offsets);
out << "</Points>\n";
+ // Write element connectivity, types and offsets
out << "<Cells>\n";
writeCells(out, offsets);
out << "</Cells>\n";
+
out << "</Piece>\n";
out << "</UnstructuredGrid>\n";
@@ -98,17 +113,17 @@ void VtkWriter<GridView>::writeImpl (std::string const& filename) const
if (is_a(format_, Vtk::APPENDED)) {
out << "<AppendedData encoding=\"raw\">\n_";
appended_pos = out.tellp();
- for (auto const& v : vertexData_) {
+ for (auto const& v : pointData_) {
if (datatype_ == Vtk::FLOAT32)
- blocks.push_back( writeDataAppended<float>(out, v, Vtk::VERTEX_DATA) );
+ blocks.push_back( writeDataAppended<float>(out, v, POINT_DATA) );
else
- blocks.push_back( writeDataAppended<double>(out, v, Vtk::VERTEX_DATA) );
+ blocks.push_back( writeDataAppended<double>(out, v, POINT_DATA) );
}
for (auto const& v : cellData_) {
if (datatype_ == Vtk::FLOAT32)
- blocks.push_back( writeDataAppended<float>(out, v, Vtk::CELL_DATA) );
+ blocks.push_back( writeDataAppended<float>(out, v, CELL_DATA) );
else
- blocks.push_back( writeDataAppended<double>(out, v, Vtk::CELL_DATA) );
+ blocks.push_back( writeDataAppended<double>(out, v, CELL_DATA) );
}
if (datatype_ == Vtk::FLOAT32)
blocks.push_back( writePointsAppended<float>(out) );
@@ -133,83 +148,21 @@ void VtkWriter<GridView>::writeImpl (std::string const& filename) const
}
-// @{ implementation details
-// TODO: allow to overload these function
-
-template <class T, class GridView>
-std::vector<T> getPoints (GridView const& gridView)
-{
- const int dim = GridView::dimension;
- std::vector<T> data(gridView.size(dim) * 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;
-}
-
-template <class T, class GridView, class GlobalFunction>
-std::vector<T> getVertexData (GridView const& gridView, GlobalFunction const& fct)
-{
- const int dim = GridView::dimension;
- std::vector<T> data(gridView.size(dim) * 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) {
- std::size_t idx = fct.ncomps() * indexSet.subIndex(e,cellType.permutation(j),dim);
- for (int comp = 0; comp < fct.ncomps(); ++comp)
- data[idx + comp] = T(localFct.evaluate(comp, refElem.position(cellType.permutation(j),dim)));
- }
- localFct.unbind();
- }
- return data;
-}
-
-template <class T, class GridView, class GlobalFunction>
-std::vector<T> getCellData (GridView const& gridView, GlobalFunction const& fct)
-{
- const int dim = GridView::dimension;
- std::vector<T> data(gridView.size(0) * fct.ncomps());
- auto const& indexSet = gridView.indexSet();
- auto localFct = localFunction(fct);
- for (auto const& e : elements(gridView)) {
- localFct.bind(e);
- auto refElem = referenceElement(e.geometry());
- std::size_t idx = fct.ncomps() * indexSet.index(e);
- for (int comp = 0; comp < fct.ncomps(); ++comp)
- data[idx + comp] = T(localFct.evaluate(comp, refElem.position(0,0)));
- localFct.unbind();
- }
- return data;
-}
-
-// @}
-
-
-template <class GridView>
-void VtkWriter<GridView>::writeData (std::ofstream& out, std::vector<pos_type>& offsets,
- GlobalFunction const& fct, Vtk::PositionTypes type) const
+template <class GV, class DC>
+void VtkWriter<GV,DC>::writeData (std::ofstream& out, std::vector<pos_type>& offsets,
+ GlobalFunction const& fct, PositionTypes type) const
{
- out << "<DataArray Name=\"" << fct.name() << "\" type=\"" << (datatype_ == Vtk::FLOAT32 ? "Float32" : "Float64") << "\""
+ out << "<DataArray Name=\"" << fct.name() << "\" type=\"" << Vtk::Map::from_datatype[datatype_] << "\""
<< " NumberOfComponents=\"" << fct.ncomps() << "\" format=\"" << (format_ == Vtk::ASCII ? "ascii\">\n" : "appended\"");
if (format_ == Vtk::ASCII) {
std::size_t i = 0;
- if (type == Vtk::VERTEX_DATA) {
- auto data = getVertexData<double>(gridView_, fct);
+ if (type == POINT_DATA) {
+ auto data = dataCollector_.template pointData<double>(fct);
for (auto const& v : data)
out << v << (++i % 6 != 0 ? ' ' : '\n');
} else {
- auto data = getCellData<double>(gridView_, fct);
+ auto data = dataCollector_.template cellData<double>(fct);
for (auto const& v : data)
out << v << (++i % 6 != 0 ? ' ' : '\n');
}
@@ -223,14 +176,14 @@ void VtkWriter<GridView>::writeData (std::ofstream& out, std::vector<pos_type>&
}
-template <class GridView>
-void VtkWriter<GridView>::writePoints (std::ofstream& out, std::vector<pos_type>& offsets) const
+template <class GV, class DC>
+void VtkWriter<GV,DC>::writePoints (std::ofstream& out, std::vector<pos_type>& offsets) const
{
- out << "<DataArray type=\"" << (datatype_ == Vtk::FLOAT32 ? "Float32" : "Float64") << "\""
+ out << "<DataArray type=\"" << Vtk::Map::from_datatype[datatype_] << "\""
<< " NumberOfComponents=\"3\" format=\"" << (format_ == Vtk::ASCII ? "ascii\">\n" : "appended\"");
if (format_ == Vtk::ASCII) {
- auto points = getPoints<double>(gridView_);
+ auto points = dataCollector_.template points<double>();
std::size_t i = 0;
for (auto const& v : points)
out << v << (++i % 6 != 0 ? ' ' : '\n');
@@ -244,31 +197,27 @@ void VtkWriter<GridView>::writePoints (std::ofstream& out, std::vector<pos_type>
}
-template <class GridView>
-void VtkWriter<GridView>::writeCells (std::ofstream& out, std::vector<pos_type>& offsets) const
+template <class GV, class DC>
+void VtkWriter<GV,DC>::writeCells (std::ofstream& out, std::vector<pos_type>& offsets) const
{
- auto const& indexSet = gridView_.indexSet();
if (format_ == Vtk::ASCII) {
+ auto cells = dataCollector_.cells();
out << "<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">\n";
std::size_t i = 0;
- for (auto const& c : elements(gridView_)) {
- auto cellType = getType(c.type());
- for (int j = 0; j < c.subEntities(dimension); ++j)
- out << indexSet.subIndex(c,cellType.permutation(j),dimension) << (++i % 6 != 0 ? ' ' : '\n');
- }
+ for (auto const& c : cells.connectivity)
+ out << c << (++i % 6 != 0 ? ' ' : '\n');
out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n";
out << "<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">\n";
i = 0;
- std::size_t old_o = 0;
- for (auto const& c : elements(gridView_))
- out << (old_o += c.subEntities(dimension)) << (++i % 6 != 0 ? ' ' : '\n');
+ for (auto const& o : cells.offsets)
+ out << o << (++i % 6 != 0 ? ' ' : '\n');
out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n";
out << "<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n";
i = 0;
- for (auto const& c : elements(gridView_))
- out << int(getType(c.type()).type()) << (++i % 6 != 0 ? ' ' : '\n');
+ for (auto const& t : cells.types)
+ out << int(t) << (++i % 6 != 0 ? ' ' : '\n');
out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n";
}
else { // Vtk::APPENDED format
@@ -293,7 +242,7 @@ void VtkWriter<GridView>::writeCells (std::ofstream& out, std::vector<pos_type>&
}
-// @{ implementation details
+namespace Impl {
template <class T>
std::uint64_t writeValuesToBuffer (std::size_t max_num_values, unsigned char* buffer,
@@ -332,11 +281,12 @@ std::uint64_t writeCompressed (unsigned char const* buffer, unsigned char* buffe
#endif
}
-// @}
+} // end namespace Impl
+
-template <class GridView>
+template <class GV, class DC>
template <class T>
-std::uint64_t VtkWriter<GridView>::writeAppended (std::ofstream& out, std::vector<T> const& values) const
+std::uint64_t VtkWriter<GV,DC>::writeAppended (std::ofstream& out, std::vector<T> const& values) const
{
assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
pos_type begin_pos = out.tellp();
@@ -366,12 +316,12 @@ std::uint64_t VtkWriter<GridView>::writeAppended (std::ofstream& out, std::vecto
std::vector<std::uint64_t> cbs(std::size_t(num_blocks), 0); // compressed block sizes
for (std::size_t i = 0; i < std::size_t(num_blocks); ++i) {
- std::uint64_t bs = writeValuesToBuffer<T>(num_values, buffer.data(), values, i*num_values);
+ std::uint64_t bs = Impl::writeValuesToBuffer<T>(num_values, buffer.data(), values, i*num_values);
if (format_ == Vtk::COMPRESSED) {
buffer_out.resize(std::size_t(compressed_block_size));
- cbs[i] = writeCompressed(buffer.data(), buffer_out.data(), bs,
- compressed_block_size, compression_level, out);
+ cbs[i] = Impl::writeCompressed(buffer.data(), buffer_out.data(), bs,
+ compressed_block_size, compression_level, out);
} else
out.write((char*)buffer.data(), bs);
}
@@ -387,72 +337,46 @@ std::uint64_t VtkWriter<GridView>::writeAppended (std::ofstream& out, std::vecto
}
-template <class GridView>
+template <class GV, class DC>
template <class T>
-std::uint64_t VtkWriter<GridView>::writeDataAppended (std::ofstream& out, GlobalFunction const& fct, Vtk::PositionTypes type) const
+std::uint64_t VtkWriter<GV,DC>::writeDataAppended (std::ofstream& out, GlobalFunction const& fct, PositionTypes type) const
{
assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
- if (type == Vtk::VERTEX_DATA) {
- auto data = getVertexData<T>(gridView_, fct);
+ if (type == POINT_DATA) {
+ auto data = dataCollector_.template pointData<T>(fct);
return writeAppended(out, data);
} else {
- auto data = getCellData<T>(gridView_, fct);
+ auto data = dataCollector_.template cellData<T>(fct);
return writeAppended(out, data);
}
}
-template <class GridView>
+template <class GV, class DC>
template <class T>
-std::uint64_t VtkWriter<GridView>::writePointsAppended (std::ofstream& out) const
+std::uint64_t VtkWriter<GV,DC>::writePointsAppended (std::ofstream& out) const
{
assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
- auto points = getPoints<T>(gridView_);
+ auto points = dataCollector_.template points<T>();
return writeAppended(out, points);
}
-template <class GridView>
-std::array<std::uint64_t,3> VtkWriter<GridView>::writeCellsAppended (std::ofstream& out) const
+template <class GV, class DC>
+std::array<std::uint64_t,3> VtkWriter<GV,DC>::writeCellsAppended (std::ofstream& out) const
{
assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
- auto const& indexSet = gridView_.indexSet();
- auto types = indexSet.types(0);
- int maxVertices = std::accumulate(types.begin(), types.end(), 1, [](int m, GeometryType t) {
- auto refElem = referenceElement<double,dimension>(t);
- return std::max(m, refElem.size(dimension));
- });
-
- std::vector<std::int64_t> connectivity;
- std::vector<std::int64_t> cell_offsets;
- std::vector<std::uint8_t> cell_types;
-
- connectivity.reserve(gridView_.size(0) * maxVertices);
- cell_offsets.reserve(gridView_.size(0));
- cell_types.reserve(gridView_.size(0));
-
- std::int64_t old_o = 0;
- for (auto const& c : elements(gridView_)) {
- auto cellType = getType(c.type());
- for (int j = 0; j < c.subEntities(dimension); ++j)
- connectivity.push_back( std::int64_t(indexSet.subIndex(c,cellType.permutation(j),dimension)) );
- cell_offsets.push_back(old_o += c.subEntities(dimension));
- cell_types.push_back(cellType.type());
- }
-
- // write conncetivity
- std::uint64_t bs0 = writeAppended(out, connectivity);
-
- // write offsets
- std::uint64_t bs1 = writeAppended(out, cell_offsets);
+ auto cells = dataCollector_.cells();
- // write cell types
- std::uint64_t bs2 = writeAppended(out, cell_types);
+ // write conncetivity, offsets, and types
+ std::uint64_t bs0 = writeAppended(out, cells.connectivity);
+ std::uint64_t bs1 = writeAppended(out, cells.offsets);
+ std::uint64_t bs2 = writeAppended(out, cells.types);
return {bs0, bs1, bs2};
}
-} // end namespace Dec
+}} // end namespace Dune::experimental
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 497e78f..a1d6e57 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -6,4 +6,16 @@ add_executable("vtkreader" vtkreader.cc)
target_link_dune_default_libraries("vtkreader")
target_link_libraries("vtkreader" dunevtk)
+add_executable("benchmark" benchmark.cc)
+target_link_dune_default_libraries("benchmark")
+target_link_libraries("benchmark" dunevtk)
+
+add_executable("datacollector" datacollector.cc)
+target_link_dune_default_libraries("datacollector")
+target_link_libraries("datacollector" dunevtk)
+
+add_executable("quadratic" quadratic.cc)
+target_link_dune_default_libraries("quadratic")
+target_link_libraries("quadratic" dunevtk)
+
add_subdirectory(test)
\ No newline at end of file
diff --git a/src/benchmark.cc b/src/benchmark.cc
new file mode 100644
index 0000000..8bcd419
--- /dev/null
+++ b/src/benchmark.cc
@@ -0,0 +1,114 @@
+// -*- 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 <cstring>
+#include <iostream>
+#include <vector>
+
+#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
+#include <dune/common/filledarray.hh>
+#include <dune/common/timer.hh>
+#include <dune/common/test/testsuite.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/io/file/vtk.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#include <dune/vtk/vtkwriter.hh>
+
+using namespace Dune;
+
+using TestCasesOld = std::set<std::tuple<std::string,VTK::OutputType,VTK::DataMode>>;
+static TestCasesOld test_cases_old = {
+ {"asciiC", VTK::ascii, VTK::conforming},
+// {"asciiNC", VTK::ascii, VTK::nonconforming},
+// {"base64C", VTK::appendedbase64, VTK::conforming},
+// {"base64NC", VTK::appendedbase64, VTK::nonconforming},
+ {"binC", VTK::appendedraw, VTK::conforming},
+// {"binaryNC", VTK::appendedraw, VTK::nonconforming}
+};
+
+
+using TestCasesNew = std::set<std::tuple<std::string,experimental::Vtk::FormatTypes,experimental::Vtk::DataTypes>>;
+static TestCasesNew test_cases_new = {
+ {"ascii32", experimental::Vtk::ASCII, experimental::Vtk::FLOAT32},
+ {"bin32", experimental::Vtk::BINARY, experimental::Vtk::FLOAT32},
+ // {"zlib32", experimental::Vtk::COMPRESSED, experimental::Vtk::FLOAT32},
+ // {"ascii64", experimental::Vtk::ASCII, experimental::Vtk::FLOAT64},
+ // {"bin64", experimental::Vtk::BINARY, experimental::Vtk::FLOAT64},
+ // {"zlib64", experimental::Vtk::COMPRESSED, experimental::Vtk::FLOAT64}
+};
+
+template <class GridView>
+void writer_old (GridView const& gridView)
+{
+ Timer t;
+ for (auto const& test_case : test_cases_old) {
+ t.reset();
+ VTKWriter<GridView> vtkWriter(gridView, std::get<2>(test_case));
+ vtkWriter.write("/tmp/writer_old_" + std::get<0>(test_case) + ".vtu",
+ std::get<1>(test_case));
+ std::cout << " time (writer_old_" + std::get<0>(test_case) + ") = " << t.elapsed() << "\n";
+ }
+}
+
+template <class GridView>
+void writer_new (GridView const& gridView)
+{
+ Timer t;
+ experimental::VtkWriter<GridView> vtkWriter(gridView);
+ for (auto const& test_case : test_cases_new) {
+ t.reset();
+ vtkWriter.write("/tmp/writer_new_" + std::get<0>(test_case) + ".vtu",
+ std::get<1>(test_case), std::get<2>(test_case));
+ std::cout << " time (writer_new_" + std::get<0>(test_case) + ") = " << t.elapsed() << "\n";
+ }
+}
+
+template <int I>
+using int_ = std::integral_constant<int,I>;
+
+int main (int argc, char** argv)
+{
+ Dune::MPIHelper::instance(argc, argv);
+
+ TestSuite test{};
+
+#ifdef HAVE_UG
+ // Test VtkWriter for UGGrid
+ std::cout << "UGGrid\n";
+ Hybrid::forEach(std::make_tuple(int_<2>{}, int_<3>{}), [&test](auto dim)
+ {
+ using GridType = UGGrid<dim.value>;
+ FieldVector<double,dim.value> lowerLeft; lowerLeft = 0.0;
+ FieldVector<double,dim.value> upperRight; upperRight = 1.0;
+ auto numElements = filledArray<dim.value,unsigned int>(10);
+ auto gridPtr = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft, upperRight, numElements);
+
+ std::cout << "DIMENSION " << dim.value << "\n";
+ writer_old(gridPtr->leafGridView());
+ writer_new(gridPtr->leafGridView());
+ });
+#endif
+
+ // Test VtkWriter for YaspGrid
+ std::cout << "YaspGrid\n";
+ Hybrid::forEach(std::make_tuple(int_<1>{}, int_<2>{}, int_<3>{}), [](auto dim)
+ {
+ using GridType = YaspGrid<dim.value>;
+ FieldVector<double,dim.value> upperRight; upperRight = 1.0;
+ auto numElements = filledArray<dim.value,int>(10);
+ GridType grid(upperRight, numElements);
+
+ std::cout << "DIMENSION " << dim.value << "\n";
+ writer_old(grid.leafGridView());
+ writer_new(grid.leafGridView());
+ });
+
+ return test.exit();
+}
diff --git a/src/datacollector.cc b/src/datacollector.cc
new file mode 100644
index 0000000..2bdd846
--- /dev/null
+++ b/src/datacollector.cc
@@ -0,0 +1,115 @@
+// -*- 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/vtkwriter.hh>
+#include <dune/vtk/vtkreader.hh>
+
+using namespace Dune;
+using namespace Dune::experimental;
+using namespace Dune::Functions;
+
+#define GRID_TYPE 2
+static const int dim = 3;
+#if GRID_TYPE == 1
+ using GridType = YaspGrid<dim>;
+#elif GRID_TYPE == 2
+ using GridType = UGGrid<dim>;
+#endif
+using GridViewType = typename GridType::LeafGridView;
+
+void write()
+{
+#if GRID_TYPE == 1
+ FieldVector<double,dim> upperRight; upperRight = 1.0;
+ auto numElements = filledArray<dim,int>(4);
+ GridType grid(upperRight,numElements);
+#elif GRID_TYPE == 2
+ FieldVector<double,dim> lowerLeft; lowerLeft = 0.0;
+ FieldVector<double,dim> upperRight; upperRight = 1.0;
+ auto numElements = filledArray<dim,unsigned int>(4);
+ auto gridPtr = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft, upperRight, numElements);
+ auto& grid = *gridPtr;
+#endif
+
+ GridViewType 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);
+ // 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);
+
+
+ { // Default DataCollector
+ using Writer = VtkWriter<GridViewType>;
+ Writer vtkWriter(gridView);
+ vtkWriter.addPointData(p1FctWrapped, "p1");
+ vtkWriter.addCellData(p1FctWrapped, "p0");
+ vtkWriter.addPointData(p1Analytic, "analytic");
+
+ vtkWriter.write("c_ascii_float32.vtu", Vtk::ASCII);
+ vtkWriter.write("c_binary_float32.vtu", Vtk::BINARY);
+ vtkWriter.write("c_compressed_float32.vtu", Vtk::COMPRESSED);
+ vtkWriter.write("c_ascii_float64.vtu", Vtk::ASCII, Vtk::FLOAT64);
+ vtkWriter.write("c_binary_float64.vtu", Vtk::BINARY, Vtk::FLOAT64);
+ vtkWriter.write("c_compressed_float64.vtu", Vtk::COMPRESSED, Vtk::FLOAT64);
+ }
+
+ { // Discontinuous DataCollector
+ using Writer = VtkWriter<GridViewType, DiscontinuousDataCollector<GridViewType>>;
+ Writer vtkWriter(gridView);
+ vtkWriter.addPointData(p1FctWrapped, "p1");
+ vtkWriter.addCellData(p1FctWrapped, "p0");
+ vtkWriter.addPointData(p1Analytic, "analytic");
+
+ vtkWriter.write("dc_ascii_float32.vtu", Vtk::ASCII);
+ vtkWriter.write("dc_binary_float32.vtu", Vtk::BINARY);
+ vtkWriter.write("dc_compressed_float32.vtu", Vtk::COMPRESSED);
+ vtkWriter.write("dc_ascii_float64.vtu", Vtk::ASCII, Vtk::FLOAT64);
+ vtkWriter.write("dc_binary_float64.vtu", Vtk::BINARY, Vtk::FLOAT64);
+ vtkWriter.write("dc_compressed_float64.vtu", Vtk::COMPRESSED, Vtk::FLOAT64);
+ }
+}
+
+int main(int argc, char** argv)
+{
+ Dune::MPIHelper::instance(argc, argv);
+
+ write();
+ {
+ auto gridPtr = VtkReader<GridType, ConnectedGridCreator>::read("dc_ascii_float32.vtu");
+ auto& grid = *gridPtr;
+
+ VtkWriter<GridViewType> vtkWriter(grid.leafGridView());
+ vtkWriter.write("c_ascii_float32_2.vtu", Vtk::ASCII);
+ }
+}
diff --git a/src/quadratic.cc b/src/quadratic.cc
new file mode 100644
index 0000000..0015eb4
--- /dev/null
+++ b/src/quadratic.cc
@@ -0,0 +1,84 @@
+// -*- 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/vtkwriter.hh>
+
+using namespace Dune;
+using namespace Dune::experimental;
+using namespace Dune::Functions;
+
+#define GRID_TYPE 2
+
+int main(int argc, char** argv)
+{
+ Dune::MPIHelper::instance(argc, argv);
+
+ const int dim = 3;
+
+#if GRID_TYPE == 1
+ using GridType = YaspGrid<dim>;
+ FieldVector<double,dim> upperRight; upperRight = 1.0;
+ auto numElements = filledArray<dim,int>(4);
+ GridType grid(upperRight,numElements);
+#elif GRID_TYPE == 2
+ using GridType = UGGrid<dim>;
+ FieldVector<double,dim> lowerLeft; lowerLeft = 0.0;
+ FieldVector<double,dim> upperRight; upperRight = 1.0;
+ auto numElements = filledArray<dim,unsigned int>(4);
+ auto gridPtr = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft, upperRight, numElements);
+ auto& grid = *gridPtr;
+#endif
+
+ 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 = VtkWriter<GridView, QuadraticDataCollector<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("p2_ascii_float32.vtu", Vtk::ASCII);
+ vtkWriter.write("p2_binary_float32.vtu", Vtk::BINARY);
+ vtkWriter.write("p2_compressed_float32.vtu", Vtk::COMPRESSED);
+ vtkWriter.write("p2_ascii_float64.vtu", Vtk::ASCII, Vtk::FLOAT64);
+ vtkWriter.write("p2_binary_float64.vtu", Vtk::BINARY, Vtk::FLOAT64);
+ vtkWriter.write("p2_compressed_float64.vtu", Vtk::COMPRESSED, Vtk::FLOAT64);
+}
diff --git a/src/test/mixed_element_test.cc b/src/test/mixed_element_test.cc
index 1d92c5e..5321d9d 100644
--- a/src/test/mixed_element_test.cc
+++ b/src/test/mixed_element_test.cc
@@ -21,6 +21,7 @@
#include <dune/vtk/vtkwriter.hh>
using namespace Dune;
+using namespace Dune::experimental;
// see https://stackoverflow.com/questions/6163611/compare-two-files
bool compare_files (std::string const& fn1, std::string const& fn2)
@@ -85,7 +86,7 @@ void reader_test (Test& test)
vtkWriter.write("/tmp/reader_writer_test_" + std::get<0>(test_case) + "_2.vtu",
std::get<1>(test_case), std::get<2>(test_case));
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"));
+ "/tmp/reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"), std::get<0>(test_case));
}
}
diff --git a/src/test/reader_writer_test.cc b/src/test/reader_writer_test.cc
index 0bcec08..deb2a98 100644
--- a/src/test/reader_writer_test.cc
+++ b/src/test/reader_writer_test.cc
@@ -21,6 +21,7 @@
#include <dune/vtk/vtkwriter.hh>
using namespace Dune;
+using namespace Dune::experimental;
// see https://stackoverflow.com/questions/6163611/compare-two-files
bool compare_files (std::string const& fn1, std::string const& fn2)
diff --git a/src/vtkreader.cc b/src/vtkreader.cc
index a6002b0..12d0732 100644
--- a/src/vtkreader.cc
+++ b/src/vtkreader.cc
@@ -19,6 +19,7 @@
#include <dune/vtk/vtkwriter.hh>
using namespace Dune;
+using namespace Dune::experimental;
int main(int argc, char** argv)
{
diff --git a/src/vtkwriter.cc b/src/vtkwriter.cc
index 9275963..7a4a44d 100644
--- a/src/vtkwriter.cc
+++ b/src/vtkwriter.cc
@@ -23,6 +23,7 @@
#include <dune/vtk/vtkwriter.hh>
using namespace Dune;
+using namespace Dune::experimental;
using namespace Dune::Functions;
#define GRID_TYPE 2
@@ -64,7 +65,7 @@ int main(int argc, char** argv)
using Writer = VtkWriter<GridView>;
Writer vtkWriter(gridView);
- vtkWriter.addVertexData(p1FctWrapped, "p1");
+ vtkWriter.addPointData(p1FctWrapped, "p1");
vtkWriter.addCellData(p1FctWrapped, "p0");
// write analytic function
@@ -72,7 +73,7 @@ int main(int argc, char** argv)
return std::sin(10*x[0])*std::cos(10*x[1])+std::sin(10*x[2]);
}, gridView);
- vtkWriter.addVertexData(p1Analytic, "analytic");
+ vtkWriter.addPointData(p1Analytic, "analytic");
vtkWriter.write("test_ascii_float32.vtu", Vtk::ASCII);
vtkWriter.write("test_binary_float32.vtu", Vtk::BINARY);
--
GitLab