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

corrected writer of vector/tensor values data

parent 6536d541
......@@ -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:
......
......@@ -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
......
......@@ -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>{});
};
......
......@@ -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;
......
......@@ -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>
......
......@@ -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();
......
......@@ -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
......
......@@ -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
......@@ -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();
......
......@@ -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)
......
// -*- 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)