From 5408049536bd88773d6a26b5b6e5335aca031cce Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Sun, 2 Sep 2018 20:30:23 +0200
Subject: [PATCH] corrected writer of vector/tensor values data
---
dune/vtk/defaultvtkfunction.hh | 32 +++++--
dune/vtk/vtkfunction.hh | 58 ++++++++-----
dune/vtk/vtktypes.hh | 3 +-
dune/vtk/vtkwriterinterface.hh | 9 +-
dune/vtk/vtkwriterinterface.impl.hh | 67 +++++++++++----
dune/vtk/writers/vtkimagedatawriter.impl.hh | 16 +---
.../writers/vtkrectilineargridwriter.impl.hh | 22 +----
.../writers/vtkstructuredgridwriter.impl.hh | 26 ++----
.../writers/vtkunstructuredgridwriter.impl.hh | 65 +++-----------
src/CMakeLists.txt | 4 +
src/vectorwriter.cc | 85 +++++++++++++++++++
src/vtkreader.cc | 2 +-
12 files changed, 234 insertions(+), 155 deletions(-)
create mode 100644 src/vectorwriter.cc
diff --git a/dune/vtk/defaultvtkfunction.hh b/dune/vtk/defaultvtkfunction.hh
index d32f8a4..796144f 100644
--- a/dune/vtk/defaultvtkfunction.hh
+++ b/dune/vtk/defaultvtkfunction.hh
@@ -4,6 +4,13 @@
namespace Dune
{
+ template <class T, int N>
+ class FieldVector;
+
+ template <class T, int N, int M>
+ class FieldMatrix;
+
+
/// Type erasure for dune-functions LocalFunction interface
template <class GridView, class LocalFunction>
class LocalFunctionWrapper final
@@ -17,9 +24,6 @@ namespace Dune
template <class F, class D>
using Range = std::decay_t<decltype(std::declval<F>()(std::declval<D>()))>;
- template <class F, class D>
- using VectorValued = decltype(std::declval<Range<F,D>>()[0u]);
-
public:
template <class LocalFct, disableCopyMove<Self, LocalFct> = 0>
LocalFunctionWrapper (LocalFct&& localFct)
@@ -38,22 +42,32 @@ namespace Dune
virtual double evaluate (int comp, LocalCoordinate const& xi) const override
{
- return evaluateImpl(comp, xi, Std::is_detected<VectorValued,LocalFunction,LocalCoordinate>{});
+ return evaluateImpl(comp, localFct_(xi));
}
private:
// Evaluate a component of a vector valued data
- double evaluateImpl (int comp, LocalCoordinate const& xi, std::true_type) const
+ template <class T, int N, int M>
+ double evaluateImpl (int comp, FieldMatrix<T,N,M> const& mat) const
+ {
+ int r = comp / 3;
+ int c = comp % 3;
+ return r < N && c < M ? mat[r][c] : 0.0;
+ }
+
+ // Evaluate a component of a vector valued data
+ template <class T, int N>
+ double evaluateImpl (int comp, FieldVector<T,N> const& vec) const
{
- auto y = localFct_(xi);
- return comp < y.size() ? y[comp] : 0.0;
+ return comp < N ? vec[comp] : 0.0;
}
// Return the scalar values
- double evaluateImpl (int comp, LocalCoordinate const& xi, std::false_type) const
+ template <class T>
+ double evaluateImpl (int comp, T const& value) const
{
assert(comp == 0);
- return localFct_(xi);
+ return value;
}
private:
diff --git a/dune/vtk/vtkfunction.hh b/dune/vtk/vtkfunction.hh
index 8fa14a7..148a065 100644
--- a/dune/vtk/vtkfunction.hh
+++ b/dune/vtk/vtkfunction.hh
@@ -2,6 +2,7 @@
#include <type_traits>
+#include <dune/common/std/optional.hh>
#include <dune/common/std/type_traits.hh>
#include "vtklocalfunction.hh"
@@ -9,17 +10,41 @@
namespace Dune
{
+ template <class T, int N>
+ class FieldVector;
+
+ template <class T, int N, int M>
+ class FieldMatrix;
+
+ namespace Impl
+ {
+ template <class T, class = void>
+ struct SizeImpl
+ : std::integral_constant<int, 1> {};
+
+ template <class T, int N>
+ struct SizeImpl<FieldVector<T,N>>
+ : std::integral_constant<int, N> {};
+
+ template <class T, int N, int M>
+ struct SizeImpl<FieldMatrix<T,N,M>>
+ : std::integral_constant<int, N*M> {};
+ }
+
+ template <class T>
+ constexpr int Size = Impl::SizeImpl<std::decay_t<T>>::value;
+
+
template <class GridView>
class VtkFunction
{
template <class F>
using HasLocalFunction = decltype(localFunction(std::declval<F>()));
- template <class F>
- using Domain = typename std::decay_t<F>::EntitySet::GlobalCoordinate;
+ using Domain = typename GridView::template Codim<0>::Entity::Geometry::LocalCoordinate;
template <class F>
- using Range = std::decay_t<decltype(std::declval<F>()(std::declval<Domain<F>>()))>;
+ using Range = std::decay_t<decltype(std::declval<F>()(std::declval<Domain>()))>;
public:
/// Constructor VtkFunction from legacy VTKFunction
@@ -32,25 +57,18 @@ namespace Dune
/// Construct VtkFunction from dune-functions GridFunction with Signature
template <class F,
- std::enable_if_t<Std::is_detected<HasLocalFunction,F>::value, int> = 0,
- std::enable_if_t<Std::is_detected<Range,F>::value,int> = 0>
- VtkFunction (F&& fct, std::string name, int ncomps = 1, Vtk::DataTypes type = Vtk::Map::type<Range<F>>)
+ std::enable_if_t<Std::is_detected<HasLocalFunction,F>::value, int> = 0>
+ VtkFunction (F&& fct, std::string name,
+ Std::optional<int> ncomps = {},
+ Std::optional<Vtk::DataTypes> type = {})
: localFct_(localFunction(std::forward<F>(fct)))
, name_(std::move(name))
- , ncomps_(ncomps)
- , type_(type)
- {}
+ {
+ using R = Range<decltype(localFunction(std::forward<F>(fct)))>;
- /// Construct VtkFunction from dune-functions GridFunction without Signature
- template <class F,
- std::enable_if_t<Std::is_detected<HasLocalFunction,F>::value, int> = 0,
- std::enable_if_t<not Std::is_detected<Range,F>::value,int> = 0>
- VtkFunction (F const& fct, std::string name, int ncomps = 1, Vtk::DataTypes type = Vtk::FLOAT32)
- : localFct_(localFunction(std::forward<F>(fct)))
- , name_(std::move(name))
- , ncomps_(ncomps)
- , type_(type)
- {}
+ ncomps_ = ncomps ? *ncomps : Size<R>;
+ type_ = type ? *type : Vtk::Map::type<R>;
+ }
VtkFunction () = default;
@@ -69,7 +87,7 @@ namespace Dune
/// Return the number of components of the Range
int ncomps () const
{
- return ncomps_;
+ return ncomps_ > 3 ? 9 : ncomps_ > 1 ? 3 : 1; // tensor, vector, scalar
}
/// Return the VTK Datatype associated with the functions range type
diff --git a/dune/vtk/vtktypes.hh b/dune/vtk/vtktypes.hh
index 250d5a8..5e50c6a 100644
--- a/dune/vtk/vtktypes.hh
+++ b/dune/vtk/vtktypes.hh
@@ -5,6 +5,7 @@
#include <string>
#include <vector>
+#include <dune/common/ftraits.hh>
#include <dune/geometry/type.hh>
namespace Dune
@@ -79,7 +80,7 @@ namespace Dune
static constexpr DataTypes typeImpl (Type<long double>) { return FLOAT64; }
template <class T>
- static constexpr DataTypes type = typeImpl(Type<T>{});
+ static constexpr DataTypes type = typeImpl(Type<typename FieldTraits<T>::field_type>{});
};
diff --git a/dune/vtk/vtkwriterinterface.hh b/dune/vtk/vtkwriterinterface.hh
index c32a062..9f1b386 100644
--- a/dune/vtk/vtkwriterinterface.hh
+++ b/dune/vtk/vtkwriterinterface.hh
@@ -5,10 +5,6 @@
#include <string>
#include <vector>
-#ifdef HAVE_ZLIB
-#include <zlib.h>
-#endif
-
#include <dune/common/std/optional.hh>
#include <dune/vtk/filewriter.hh>
#include <dune/vtk/vtkfunction.hh>
@@ -121,6 +117,11 @@ namespace Dune
template <class T>
std::uint64_t writeValuesAppended (std::ofstream& out, std::vector<T> const& values) const;
+ template <class T>
+ void writeValuesAscii (std::ofstream& out, std::vector<T> const& values) const;
+
+ void writeHeader (std::ofstream& out, std::string const& type) const;
+
/// Return PointData/CellData attributes for the name of the first scalar/vector/tensor DataArray
std::string getNames (std::vector<VtkFunction> const& data) const;
diff --git a/dune/vtk/vtkwriterinterface.impl.hh b/dune/vtk/vtkwriterinterface.impl.hh
index 4242cf6..8a6b9fa 100644
--- a/dune/vtk/vtkwriterinterface.impl.hh
+++ b/dune/vtk/vtkwriterinterface.impl.hh
@@ -8,6 +8,10 @@
#include <sstream>
#include <string>
+#ifdef HAVE_ZLIB
+#include <zlib.h>
+#endif
+
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
@@ -82,17 +86,11 @@ void VtkWriterInterface<GV,DC>
if (format_ == Vtk::ASCII) {
out << ">\n";
- std::size_t i = 0;
- 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 = dataCollector_.template cellData<double>(fct);
- for (auto const& v : data)
- out << v << (++i % 6 != 0 ? ' ' : '\n');
- }
- out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n";
+ if (type == POINT_DATA)
+ writeValuesAscii(out, dataCollector_.template pointData<double>(fct));
+ else
+ writeValuesAscii(out, dataCollector_.template cellData<double>(fct));
+ out << "</DataArray>\n";
} else {
out << " offset=";
offsets.push_back(out.tellp());
@@ -114,11 +112,8 @@ void VtkWriterInterface<GV,DC>
if (format_ == Vtk::ASCII) {
out << ">\n";
- auto points = dataCollector_.template points<double>();
- std::size_t i = 0;
- for (auto const& v : points)
- out << v << (++i % 6 != 0 ? ' ' : '\n');
- out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n";
+ writeValuesAscii(out, dataCollector_.template points<double>());
+ out << "</DataArray>\n";
} else {
out << " offset=";
offsets.push_back(out.tellp());
@@ -168,6 +163,46 @@ void VtkWriterInterface<GV,DC>
}
+namespace Impl {
+
+ template <class T, std::enable_if_t<(sizeof(T)>1), int> = 0>
+ T const& printable (T const& t) { return t; }
+
+ std::int16_t printable (std::int8_t c) { return std::int16_t(c); }
+ std::uint16_t printable (std::uint8_t c) { return std::uint16_t(c); }
+
+} // end namespace Impl
+
+
+template <class GV, class DC>
+ template <class T>
+void VtkWriterInterface<GV,DC>
+ ::writeValuesAscii (std::ofstream& out, std::vector<T> const& values) const
+{
+ assert(is_a(format_, Vtk::ASCII) && "Function should by called only in ascii mode!\n");
+ std::size_t i = 0;
+ for (auto const& v : values)
+ out << Impl::printable(v) << (++i % 6 != 0 ? ' ' : '\n');
+ if (i % 6 != 0)
+ out << '\n';
+}
+
+template <class GV, class DC>
+void VtkWriterInterface<GV,DC>
+ ::writeHeader (std::ofstream& out, std::string const& type) const
+{
+ out << "<VTKFile"
+ << " type=\"" << type << "\""
+ << " version=\"1.0\""
+ << " header_type=\"UInt64\"";
+ if (format_ != Vtk::ASCII)
+ out << " byte_order=\"" << getEndian() << "\"";
+ if (format_ == Vtk::COMPRESSED)
+ out << " compressor=\"vtkZLibDataCompressor\"";
+ out << ">\n";
+}
+
+
namespace Impl {
template <class T>
diff --git a/dune/vtk/writers/vtkimagedatawriter.impl.hh b/dune/vtk/writers/vtkimagedatawriter.impl.hh
index f37691c..0eb3c38 100644
--- a/dune/vtk/writers/vtkimagedatawriter.impl.hh
+++ b/dune/vtk/writers/vtkimagedatawriter.impl.hh
@@ -21,13 +21,7 @@ void VtkImageDataWriter<GV,DC>
::writeSerialFile (std::ofstream& out) const
{
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";
+ this->writeHeader(out, "ImageData");
auto const& wholeExtent = dataCollector_.wholeExtent();
auto const& origin = dataCollector_.origin();
@@ -66,13 +60,7 @@ template <class GV, class DC>
void VtkImageDataWriter<GV,DC>
::writeParallelFile (std::ofstream& out, std::string const& pfilename, int /*size*/) const
{
- out << "<VTKFile"
- << " type=\"PImageData\""
- << " version=\"1.0\""
- << " byte_order=\"" << this->getEndian() << "\""
- << " header_type=\"UInt64\""
- << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\"" : "")
- << ">\n";
+ this->writeHeader(out, "PImageData");
auto const& wholeExtent = dataCollector_.wholeExtent();
auto const& origin = dataCollector_.origin();
diff --git a/dune/vtk/writers/vtkrectilineargridwriter.impl.hh b/dune/vtk/writers/vtkrectilineargridwriter.impl.hh
index 3ac7235..8bad8f4 100644
--- a/dune/vtk/writers/vtkrectilineargridwriter.impl.hh
+++ b/dune/vtk/writers/vtkrectilineargridwriter.impl.hh
@@ -21,13 +21,7 @@ void VtkRectilinearGridWriter<GV,DC>
::writeSerialFile (std::ofstream& out) const
{
std::vector<pos_type> offsets; // pos => offset
- out << "<VTKFile"
- << " type=\"RectilinearGrid\""
- << " version=\"1.0\""
- << " byte_order=\"" << this->getEndian() << "\""
- << " header_type=\"UInt64\""
- << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\"" : "")
- << ">\n";
+ this->writeHeader(out, "RectilinearGrid");
auto const& wholeExtent = dataCollector_.wholeExtent();
out << "<RectilinearGrid"
@@ -67,13 +61,7 @@ template <class GV, class DC>
void VtkRectilinearGridWriter<GV,DC>
::writeParallelFile (std::ofstream& out, std::string const& pfilename, int /*size*/) const
{
- out << "<VTKFile"
- << " type=\"PRectilinearGrid\""
- << " version=\"1.0\""
- << " byte_order=\"" << this->getEndian() << "\""
- << " header_type=\"UInt64\""
- << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\"" : "")
- << ">\n";
+ this->writeHeader(out, "PRectilinearGrid");
auto const& wholeExtent = dataCollector_.wholeExtent();
out << "<PRectilinearGrid"
@@ -138,10 +126,8 @@ void VtkRectilinearGridWriter<GV,DC>
if (timestep)
out << " TimeStep=\"" << *timestep << "\"";
out << ">\n";
- std::size_t i = 0;
- for (auto const& c : coordinates[d])
- out << c << (++i % 6 != 0 ? ' ' : '\n');
- out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n";
+ this->writeValuesAscii(out, coordinates[d]);
+ out << "</DataArray>\n";
}
}
else { // Vtk::APPENDED format
diff --git a/dune/vtk/writers/vtkstructuredgridwriter.impl.hh b/dune/vtk/writers/vtkstructuredgridwriter.impl.hh
index 30e6323..baf7d7e 100644
--- a/dune/vtk/writers/vtkstructuredgridwriter.impl.hh
+++ b/dune/vtk/writers/vtkstructuredgridwriter.impl.hh
@@ -21,13 +21,7 @@ void VtkStructuredGridWriter<GV,DC>
::writeSerialFile (std::ofstream& out) const
{
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";
+ this->writeHeader(out, "StructuredGrid");
auto const& wholeExtent = dataCollector_.wholeExtent();
out << "<StructuredGrid WholeExtent=\"" << join(wholeExtent.begin(), wholeExtent.end()) << "\">\n";
@@ -65,13 +59,7 @@ template <class GV, class DC>
void VtkStructuredGridWriter<GV,DC>
::writeParallelFile (std::ofstream& out, std::string const& pfilename, int /*size*/) const
{
- out << "<VTKFile"
- << " type=\"PStructuredGrid\""
- << " version=\"1.0\""
- << " byte_order=\"" << this->getEndian() << "\""
- << " header_type=\"UInt64\""
- << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\"" : "")
- << ">\n";
+ this->writeHeader(out, "PStructuredGrid");
auto const& wholeExtent = dataCollector_.wholeExtent();
out << "<PStructuredGrid"
@@ -131,13 +119,9 @@ void VtkStructuredGridWriter<GV,DC>
assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
// write points
- if (datatype_ == Vtk::FLOAT32) {
- auto points = dataCollector_.template points<float>();
- blocks.push_back(this->writeValuesAppended(out, points));
- } else {
- auto points = dataCollector_.template points<double>();
- blocks.push_back(this->writeValuesAppended(out, points));
- }
+ blocks.push_back( datatype_ == Vtk::FLOAT32
+ ? this->writeValuesAppended(out, dataCollector_.template points<float>())
+ : this->writeValuesAppended(out, dataCollector_.template points<double>()) );
}
} // end namespace Dune
diff --git a/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh b/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh
index 0a4e8b7..d915751 100644
--- a/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh
+++ b/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh
@@ -21,15 +21,9 @@ void VtkUnstructuredGridWriter<GV,DC>
::writeSerialFile (std::ofstream& out) const
{
std::vector<pos_type> offsets; // pos => offset
- out << "<VTKFile"
- << " type=\"UnstructuredGrid\""
- << " version=\"1.0\""
- << " byte_order=\"" << this->getEndian() << "\""
- << " header_type=\"UInt64\""
- << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\"" : "")
- << ">\n";
-
+ this->writeHeader(out, "UnstructuredGrid");
out << "<UnstructuredGrid>\n";
+
out << "<Piece"
<< " NumberOfPoints=\"" << dataCollector_.numPoints() << "\""
<< " NumberOfCells=\"" << dataCollector_.numCells() << "\""
@@ -69,14 +63,7 @@ template <class GV, class DC>
void VtkUnstructuredGridWriter<GV,DC>
::writeParallelFile (std::ofstream& out, std::string const& pfilename, int size) const
{
- out << "<VTKFile"
- << " type=\"PUnstructuredGrid\""
- << " version=\"1.0\""
- << " byte_order=\"" << this->getEndian() << "\""
- << " header_type=\"UInt64\""
- << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\"" : "")
- << ">\n";
-
+ this->writeHeader(out, "PUnstructuredGrid");
out << "<PUnstructuredGrid GhostLevel=\"0\">\n";
// Write points
@@ -130,14 +117,7 @@ void VtkUnstructuredGridWriter<GV,DC>
assert(is_a(format_, Vtk::APPENDED));
std::vector<std::vector<pos_type>> offsets(timesteps.size()); // pos => offset
- out << "<VTKFile"
- << " type=\"UnstructuredGrid\""
- << " version=\"1.0\""
- << " byte_order=\"" << this->getEndian() << "\""
- << " header_type=\"UInt64\""
- << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\"" : "")
- << ">\n";
-
+ this->writeHeader(out, "UnstructuredGrid");
out << "<UnstructuredGrid"
<< " TimeValues=\"";
{
@@ -237,14 +217,7 @@ void VtkUnstructuredGridWriter<GV,DC>
int size,
std::vector<std::pair<double, std::string>> const& timesteps) const
{
- out << "<VTKFile"
- << " type=\"PUnstructuredGrid\""
- << " version=\"1.0\""
- << " byte_order=\"" << this->getEndian() << "\""
- << " header_type=\"UInt64\""
- << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\"" : "")
- << ">\n";
-
+ this->writeHeader(out, "PUnstructuredGrid");
out << "<PUnstructuredGrid GhostLevel=\"0\""
<< " TimeValues=\"";
{
@@ -312,28 +285,22 @@ void VtkUnstructuredGridWriter<GV,DC>
if (timestep)
out << " TimeStep=\"" << *timestep << "\"";
out << ">\n";
- std::size_t i = 0;
- for (auto const& c : cells.connectivity)
- out << c << (++i % 6 != 0 ? ' ' : '\n');
- out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n";
+ this->writeValuesAscii(out, cells.connectivity);
+ out << "</DataArray>\n";
out << "<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\"";
if (timestep)
out << " TimeStep=\"" << *timestep << "\"";
out << ">\n";
- i = 0;
- for (auto const& o : cells.offsets)
- out << o << (++i % 6 != 0 ? ' ' : '\n');
- out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n";
+ this->writeValuesAscii(out, cells.offsets);
+ out << "</DataArray>\n";
out << "<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\"";
if (timestep)
out << " TimeStep=\"" << *timestep << "\"";
out << ">\n";
- i = 0;
- for (auto const& t : cells.types)
- out << int(t) << (++i % 6 != 0 ? ' ' : '\n');
- out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n";
+ this->writeValuesAscii(out, cells.types);
+ out << "</DataArray>\n";
}
else { // Vtk::APPENDED format
out << "<DataArray type=\"Int64\" Name=\"connectivity\" format=\"appended\"";
@@ -370,13 +337,9 @@ void VtkUnstructuredGridWriter<GV,DC>
assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
// write points
- if (datatype_ == Vtk::FLOAT32) {
- auto points = dataCollector_.template points<float>();
- blocks.push_back(this->writeValuesAppended(out, points));
- } else {
- auto points = dataCollector_.template points<double>();
- blocks.push_back(this->writeValuesAppended(out, points));
- }
+ blocks.push_back( datatype_ == Vtk::FLOAT32
+ ? this->writeValuesAppended(out, dataCollector_.template points<float>())
+ : this->writeValuesAppended(out, dataCollector_.template points<double>()) );
// write conncetivity, offsets, and types
auto cells = dataCollector_.cells();
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 82056a0..dcab3e6 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -35,6 +35,10 @@ if (dune-functions_FOUND)
add_executable("pvdwriter" pvdwriter.cc)
target_link_dune_default_libraries("pvdwriter")
target_link_libraries("pvdwriter" dunevtk)
+
+ add_executable("vectorwriter" vectorwriter.cc)
+ target_link_dune_default_libraries("vectorwriter")
+ target_link_libraries("vectorwriter" dunevtk)
endif ()
if (dune-polygongrid_FOUND)
diff --git a/src/vectorwriter.cc b/src/vectorwriter.cc
new file mode 100644
index 0000000..c7eac50
--- /dev/null
+++ b/src/vectorwriter.cc
@@ -0,0 +1,85 @@
+// -*- 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 <set>
+#include <vector>
+
+#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
+#include <dune/common/exceptions.hh> // We use exceptions
+
+#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/yaspgrid.hh>
+
+#include <dune/vtk/vtkwriter.hh>
+
+using namespace Dune;
+using namespace Dune::Functions;
+
+using TestCases = std::set<std::tuple<std::string,Vtk::FormatTypes,Vtk::DataTypes>>;
+static TestCases test_cases = {
+ {"ascii32", Vtk::ASCII, Vtk::FLOAT32},
+ {"bin32", Vtk::BINARY, Vtk::FLOAT32},
+ {"zlib32", Vtk::COMPRESSED, Vtk::FLOAT32},
+ {"ascii64", Vtk::ASCII, Vtk::FLOAT64},
+ {"bin64", Vtk::BINARY, Vtk::FLOAT64},
+ {"zlib64", Vtk::COMPRESSED, Vtk::FLOAT64},
+};
+
+template <class GridView>
+void write (std::string prefix, GridView const& gridView)
+{
+ using namespace BasisFactory;
+ auto basis = makeBasis(gridView, lagrange<1>());
+
+ FieldVector<double,GridView::dimensionworld> c;
+ if (GridView::dimensionworld > 0) c[0] = 11.0;
+ if (GridView::dimensionworld > 1) c[1] = 7.0;
+ if (GridView::dimensionworld > 2) c[2] = 3.0;
+
+ auto vec = [&c](auto const& x) {
+ FieldVector<double,GridView::dimensionworld> result; result = c.dot(x);
+ return result;
+ };
+ auto mat = [&c](auto const& x) {
+ FieldMatrix<double,GridView::dimensionworld,GridView::dimensionworld> result;
+ for (int i = 0; i < GridView::dimensionworld; ++i)
+ result[i][i] = c.dot(x);
+ return result;
+ };
+ auto vec_fct = makeAnalyticGridViewFunction(vec, gridView);
+ auto mat_fct = makeAnalyticGridViewFunction(mat, gridView);
+
+ for (auto const& test_case : test_cases) {
+ VtkWriter<GridView> vtkWriter(gridView, std::get<1>(test_case), std::get<2>(test_case));
+ vtkWriter.addPointData(vec_fct, "vec");
+ vtkWriter.addPointData(mat_fct, "mat");
+ vtkWriter.write(prefix + "_" + std::to_string(GridView::dimensionworld) + "d_" + std::get<0>(test_case) + ".vtu");
+ }
+}
+
+template <int I>
+using int_ = std::integral_constant<int,I>;
+
+int main (int argc, char** argv)
+{
+ Dune::MPIHelper::instance(argc, argv);
+
+ // Test VtkWriter for YaspGrid
+ Hybrid::forEach(std::make_tuple(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>(8);
+ GridType grid(upperRight, numElements, 0, 0);
+ write("yasp_vec", grid.leafGridView());
+ });
+}
\ No newline at end of file
diff --git a/src/vtkreader.cc b/src/vtkreader.cc
index 8453d99..d82e96b 100644
--- a/src/vtkreader.cc
+++ b/src/vtkreader.cc
@@ -31,7 +31,7 @@ int main(int argc, char** argv)
{
FieldVector<double,dim> lowerLeft; lowerLeft = 0.0;
FieldVector<double,dim> upperRight; upperRight = 1.0;
- auto numElements = filledArray<dim,unsigned int>(8);
+ auto numElements = filledArray<dim,unsigned int>(4);
auto gridPtr = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft, upperRight, numElements);
auto& grid = *gridPtr;
--
GitLab