Commit 390d3b97 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Structured grid writer and data collector extended to support YaspGrid and SPGrid

parent 24721938
......@@ -8,3 +8,4 @@ Version: 0.1
Maintainer: simon.praetorius@tu-dresden.de
#depending on
Depends: dune-grid dune-functions
Suggests: dune-spgrid
......@@ -113,7 +113,7 @@ protected: // default implementations
std::vector<T> data(numCells() * fct.ncomps());
auto const& indexSet = gridView_.indexSet();
auto localFct = localFunction(fct);
for (auto const& e : elements(gridView_)) {
for (auto const& e : elements(gridView_, Partitions::all)) {
localFct.bind(e);
auto geometry = e.geometry();
std::size_t idx = fct.ncomps() * indexSet.index(e);
......
......@@ -33,7 +33,7 @@ public:
{
std::vector<T> data(this->numPoints() * 3);
auto const& indexSet = gridView_.indexSet();
for (auto const& vertex : vertices(gridView_)) {
for (auto const& vertex : vertices(gridView_, Partitions::all)) {
std::size_t idx = 3 * indexSet.index(vertex);
auto v = vertex.geometry().center();
for (std::size_t j = 0; j < v.size(); ++j)
......@@ -61,7 +61,7 @@ public:
cells.types.reserve(this->numCells());
std::int64_t old_o = 0;
for (auto const& c : elements(gridView_)) {
for (auto const& c : elements(gridView_, Partitions::all)) {
Vtk::CellType cellType(c.type());
for (int j = 0; j < c.subEntities(dim); ++j)
cells.connectivity.push_back( std::int64_t(indexSet.subIndex(c,cellType.permutation(j),dim)) );
......@@ -79,7 +79,7 @@ public:
std::vector<T> data(this->numPoints() * fct.ncomps());
auto const& indexSet = gridView_.indexSet();
auto localFct = localFunction(fct);
for (auto const& e : elements(gridView_)) {
for (auto const& e : elements(gridView_, Partitions::all)) {
localFct.bind(e);
Vtk::CellType cellType{e.type()};
auto refElem = referenceElement(e.geometry());
......
......@@ -30,7 +30,7 @@ public:
indexMap_.resize(gridView_.size(dim));
std::int64_t vertex_idx = 0;
auto const& indexSet = gridView_.indexSet();
for (auto const& c : elements(gridView_)) {
for (auto const& c : elements(gridView_, Partitions::interior)) {
numPoints_ += c.subEntities(dim);
for (int i = 0; i < c.subEntities(dim); ++i)
indexMap_[indexSet.subIndex(c, i, dim)] = vertex_idx++;
......@@ -49,7 +49,7 @@ public:
{
std::vector<T> data(this->numPoints() * 3);
auto const& indexSet = gridView_.indexSet();
for (auto const& element : elements(gridView_)) {
for (auto const& element : elements(gridView_, Partitions::interior)) {
for (int i = 0; i < element.subEntities(dim); ++i) {
std::size_t idx = 3 * indexMap_[indexSet.subIndex(element, i, dim)];
auto v = element.geometry().corner(i);
......@@ -72,7 +72,7 @@ public:
std::int64_t old_o = 0;
auto const& indexSet = gridView_.indexSet();
for (auto const& c : elements(gridView_)) {
for (auto const& c : elements(gridView_, Partitions::interior)) {
Vtk::CellType cellType(c.type());
for (int j = 0; j < c.subEntities(dim); ++j) {
std::int64_t vertex_idx = indexMap_[indexSet.subIndex(c,cellType.permutation(j),dim)];
......@@ -92,7 +92,7 @@ public:
std::vector<T> data(this->numPoints() * fct.ncomps());
auto const& indexSet = gridView_.indexSet();
auto localFct = localFunction(fct);
for (auto const& e : elements(gridView_)) {
for (auto const& e : elements(gridView_, Partitions::interior)) {
localFct.bind(e);
Vtk::CellType cellType{e.type()};
auto refElem = referenceElement(e.geometry());
......
......@@ -37,7 +37,7 @@ public:
{
std::vector<T> data(this->numPoints() * 3);
auto const& indexSet = gridView_.indexSet();
for (auto const& element : elements(gridView_)) {
for (auto const& element : elements(gridView_, Partitions::interior)) {
auto geometry = element.geometry();
auto refElem = referenceElement<T,dim>(element.type());
......@@ -77,7 +77,7 @@ public:
std::int64_t old_o = 0;
auto const& indexSet = gridView_.indexSet();
for (auto const& c : elements(gridView_)) {
for (auto const& c : elements(gridView_, Partitions::interior)) {
Vtk::CellType cellType(c.type(), Vtk::QUADRATIC);
for (int j = 0; j < c.subEntities(dim); ++j) {
int k = cellType.permutation(j);
......@@ -103,7 +103,7 @@ public:
std::vector<T> data(this->numPoints() * fct.ncomps());
auto const& indexSet = gridView_.indexSet();
auto localFct = localFunction(fct);
for (auto const& e : elements(gridView_)) {
for (auto const& e : elements(gridView_, Partitions::interior)) {
localFct.bind(e);
Vtk::CellType cellType{e.type(), Vtk::QUADRATIC};
auto refElem = referenceElement(e.geometry());
......
#pragma once
#ifdef HAVE_DUNE_SPGRID
#include <dune/grid/spgrid.hh>
#endif
#include <dune/vtk/datacollectors/structureddatacollector.hh>
namespace Dune { namespace experimental
{
#ifdef HAVE_DUNE_SPGRID
// Specialization for SPGrid
template <class GridView>
class SPDataCollector
: public StructuredDataCollectorInterface<GridView, SPDataCollector<GridView>>
{
enum { dim = GridView::dimension };
using Self = SPDataCollector;
using Super = StructuredDataCollectorInterface<GridView, Self>;
using Super::gridView_;
using ctype = typename GridView::ctype;
public:
SPDataCollector (GridView const& gridView)
: Super(gridView)
, wholeExtent_(filledArray<6,int>(0))
, origin_(0.0)
, spacing_(0.0)
{}
std::array<int, 6> const& wholeExtentImpl () const
{
return wholeExtent_;
}
auto const& originImpl () const
{
return origin_;
}
auto const& spacingImpl () const
{
return spacing_;
}
void updateImpl ()
{
Super::updateImpl();
auto const& begin = gridView_.impl().gridLevel().globalMesh().begin();
auto const& end = gridView_.impl().gridLevel().globalMesh().end();
auto const& cube = gridView_.grid().domain().cube();
for (int i = 0; i < dim; ++i) {
wholeExtent_[2*i] = begin[i];
wholeExtent_[2*i+1] = end[i];
spacing_[i] = cube.width()[i] / (end[i] - begin[i]);
origin_[i] = cube.origin()[i];
}
}
template <class Writer>
void writeLocalPieceImpl (Writer const& writer) const
{
auto extent = filledArray<6,int>(0);
auto const& localMesh = gridView_.impl().gridLevel().localMesh();
for (int i = 0; i < dim; ++i) {
extent[2*i] = localMesh.begin()[i];
extent[2*i+1] = localMesh.end()[i];
}
writer(extent);
}
template <class Writer>
void writePiecesImpl (Writer const& writer) const
{
auto extent = filledArray<6,int>(0);
auto const& partitionList = gridView_.impl().gridLevel().decomposition_;
// for (auto const& part : partitionList)
for (std::size_t p = 0; p < partitionList.size(); ++p)
{
auto const& part = partitionList[p];
for (int i = 0; i < dim; ++i) {
extent[2*i] = part.begin()[i];
extent[2*i+1] = part.begin()[i] + part.width(i);
}
writer(p, extent, true);
}
}
private:
std::array<int, 6> wholeExtent_;
std::array<int, 6> extent_;
FieldVector<ctype,3> spacing_;
FieldVector<ctype,3> origin_;
std::vector<std::size_t> indexMap_;
};
#endif
}} // end namespace Dune::experimental
#pragma once
#include <dune/vtk/datacollector.hh>
#include <dune/vtk/datacollectors/continuousdatacollector.hh>
namespace Dune { namespace experimental
{
/// Implementation of \ref DataCollector for linear cells, with continuous data.
template <class GridView>
class StructuredDataCollector
: public DataCollectorInterface<GridView, StructuredDataCollector<GridView>>
template <class GridView, class Derived>
class StructuredDataCollectorInterface
: public DataCollectorInterface<GridView, Derived>
{
enum { dim = GridView::dimension };
using Self = StructuredDataCollector;
using Super = DataCollectorInterface<GridView, Self>;
protected:
using Super = DataCollectorInterface<GridView, Derived>;
using Super::gridView_;
using ctype = typename GridView::ctype;
struct CoordLess
{
template <class T, int N>
bool operator() (FieldVector<T,N> const& lhs, FieldVector<T,N> const& rhs) const
{
for (int i = N-1; i >= 0; --i) {
if (std::abs(lhs[i] - rhs[i]) < std::numeric_limits<T>::epsilon())
continue;
return lhs[i] < rhs[i];
}
return false;
}
};
public:
StructuredDataCollector (GridView const& gridView)
StructuredDataCollectorInterface (GridView const& gridView)
: Super(gridView)
, defaultDataCollector_(gridView)
{}
/// Return number of grid vertices
std::uint64_t numPointsImpl () const
{
return gridView_.size(dim);
return gridView_.size(GridView::dimension);
}
std::array<int, 6> const& extent () const
void updateImpl ()
{
return extent_;
defaultDataCollector_.update();
}
auto const& origin () const
std::array<int, 6> const& wholeExtent () const
{
return origin_;
return this->asDerived().wholeExtentImpl();
}
auto const& spacing () const
FieldVector<ctype, 3> const& origin () const
{
return spacing_;
return this->asDerived().originImpl();
}
void updateImpl ()
FieldVector<ctype, 3> const& spacing () const
{
// TODO: extract this information from the grid
int extent = GridView::dimension == 1 ? gridView_.size(dim) :
GridView::dimension == 2 ? isqrt(gridView_.size(dim)) :
GridView::dimension == 3 ? icbrt(gridView_.size(dim)) : 0;
for (int i = 0; i < 3; ++i) {
if (GridView::dimension > i) {
extent_[2*i] = 0;
extent_[2*i+1] = extent-1;
spacing_[i] = gridView_.grid().domainSize()[i] / (extent-1);
} else {
extent_[2*i] = extent_[2*i+1] = 0;
spacing_[i] = 0;
}
return this->asDerived().spacingImpl();
}
origin_[i] = 0;
}
template <class Writer>
void writeLocalPiece (Writer const& writer) const
{
this->asDerived().writeLocalPieceImpl(writer);
}
template <class Writer>
void writePieces (Writer const& writer) const
{
this->asDerived().writePiecesImpl(writer);
}
/// Return the coordinates of all grid vertices in the order given by the indexSet
template <class T>
std::vector<T> pointsImpl () const
{
std::vector<T> data(this->numPoints() * 3);
auto const& indexSet = gridView_.indexSet();
for (auto const& vertex : vertices(gridView_)) {
std::size_t idx = 3 * indexSet.index(vertex);
auto v = vertex.geometry().center();
for (std::size_t j = 0; j < v.size(); ++j)
data[idx + j] = T(v[j]);
for (std::size_t j = v.size(); j < 3u; ++j)
data[idx + j] = T(0);
}
return data;
return defaultDataCollector_.template points<T>();
}
/// Evaluate the `fct` at the corners of the elements
template <class T, class GlobalFunction>
std::vector<T> pointDataImpl (GlobalFunction const& fct) const
{
std::vector<T> data(this->numPoints() * fct.ncomps());
auto const& indexSet = gridView_.indexSet();
auto localFct = localFunction(fct);
for (auto const& e : elements(gridView_)) {
localFct.bind(e);
Vtk::CellType cellType{e.type()};
auto refElem = referenceElement(e.geometry());
for (int j = 0; j < e.subEntities(dim); ++j) {
int k = cellType.permutation(j);
std::size_t idx = fct.ncomps() * indexSet.subIndex(e,k,dim);
for (int comp = 0; comp < fct.ncomps(); ++comp)
data[idx + comp] = T(localFct.evaluate(comp, refElem.position(k,dim)));
}
localFct.unbind();
}
return data;
}
private:
static constexpr std::uint32_t isqrt (std::uint64_t x)
{
std::uint32_t c = 0x8000;
std::uint32_t g = 0x8000;
while (true) {
if (g*g > x)
g ^= c;
c >>= 1;
if (c == 0)
return g;
g |= c;
}
}
static constexpr std::uint32_t icbrt (std::uint64_t x)
{
std::uint32_t y = 0;
for (int s = 63; s >= 0; s -= 3) {
y += y;
std::uint64_t b = 3*y*(std::uint64_t(y) + 1) + 1;
if ((x >> s) >= b) {
x -= b << s;
y++;
}
}
return y;
return defaultDataCollector_.template pointData<T>(fct);
}
private:
std::array<int, 6> extent_;
FieldVector<ctype,3> spacing_;
FieldVector<ctype,3> origin_;
std::vector<std::size_t> indexMap_;
DefaultDataCollector<GridView> defaultDataCollector_;
};
}} // end namespace Dune::experimental
#pragma once
#include <dune/grid/yaspgrid.hh>
#include <dune/vtk/datacollectors/structureddatacollector.hh>
namespace Dune { namespace experimental
{
// Specialization for YaspGrid
template <class GridView>
class YaspDataCollector
: public StructuredDataCollectorInterface<GridView, YaspDataCollector<GridView>>
{
enum { dim = GridView::dimension };
using Self = YaspDataCollector;
using Super = StructuredDataCollectorInterface<GridView, Self>;
using Super::gridView_;
using ctype = typename GridView::ctype;
public:
YaspDataCollector (GridView const& gridView)
: Super(gridView)
, wholeExtent_(filledArray<6,int>(0))
, origin_(0.0)
, spacing_(0.0)
, level_(gridView.template begin<0,Interior_Partition>()->level())
{}
std::array<int, 6> const& wholeExtentImpl () const
{
return wholeExtent_;
}
auto const& originImpl () const
{
return origin_;
}
auto const& spacingImpl () const
{
return spacing_;
}
void updateImpl ()
{
Super::updateImpl();
for (int i = 0; i < dim; ++i) {
wholeExtent_[2*i] = 0;
wholeExtent_[2*i+1] = gridView_.grid().levelSize(level_,i);
}
auto it = gridView_.grid().begin(level_);
initGeometry(it->coords);
}
template <class Coords>
void initGeometry (Coords const& coords) { DUNE_THROW(NotImplemented, "Coordinate-Type not implemented!"); }
void initGeometry (EquidistantCoordinates<ctype,dim> const& coords)
{
for (int i = 0; i < dim; ++i) {
spacing_[i] = coords.meshsize(i,0);
origin_[i] = 0;
}
}
void initGeometry (EquidistantOffsetCoordinates<ctype,dim> const& coords)
{
for (int i = 0; i < dim; ++i) {
spacing_[i] = coords.meshsize(i,0);
origin_[i] = coords.origin(i);
}
}
void initGeometry (TensorProductCoordinates<ctype,dim> const& coords)
{
for (int i = 0; i < dim; ++i) {
spacing_[i] = coords.meshsize(i,0); // is not constant, but also not used.
origin_[i] = coords.coordinate(i,0);
}
}
template <class Writer>
void writeLocalPieceImpl (Writer const& writer) const
{
auto const& gl = *gridView_.grid().begin(level_);
auto const& g = gl.interior[0];
auto extent = filledArray<6,int>(0);
auto const& gc = *g.dataBegin();
for (int i = 0; i < dim; ++i) {
extent[2*i] = gc.min(i);
extent[2*i+1] = gc.max(i)+1;
}
writer(extent);
}
template <class Writer>
void writePiecesImpl (Writer const& writer) const
{
auto extent = filledArray<6,int>(0);
// for (auto const& part : gridView_.gridLevel().template partition<InteriorEntity>())
// {
// for (int i = 0; i < dim; ++i) {
// extent[2*i] = part.begin()[i];
// extent[2*i+1] = part.end()[i]-1;
// }
int num_ranks = -1;
MPI_Comm_size(MPI_COMM_WORLD, &num_ranks);
for (int p = 0; p < num_ranks; ++p) {
writer(p, extent, false);
}
}
private:
std::array<int, 6> wholeExtent_;
FieldVector<ctype,3> spacing_;
FieldVector<ctype,3> origin_;
std::vector<std::size_t> indexMap_;
int level_;
};
}} // end namespace Dune::experimental
......@@ -14,7 +14,7 @@
namespace Dune { namespace experimental
{
/// File-Writer for VTK .vtu files
template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
template <class GridView, class DataCollector>
class VtkImageDataWriter
: public VtkWriter<GridView, DataCollector>
{
......
......@@ -28,16 +28,14 @@ void VtkImageDataWriter<GV,DC>
out << std::setprecision(std::numeric_limits<double>::digits10+2);
}
dataCollector_.update();
std::vector<pos_type> offsets; // pos => offset
out << "<VTKFile type=\"ImageData\" version=\"1.0\" "
<< "byte_order=\"" << this->getEndian() << "\" header_type=\"UInt64\""
<< (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\">\n" : ">\n");
out << "<ImageData WholeExtent=\"";
auto const& extent = dataCollector_.extent();
auto const& wholeExtent = dataCollector_.wholeExtent();
for (int i = 0; i < 3; ++i) {
out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
out << (i == 0 ? "" : " ") << wholeExtent[2*i] << " " << wholeExtent[2*i+1];
}
out << "\" Origin=\"";
auto const& origin = dataCollector_.origin();
......@@ -50,11 +48,15 @@ void VtkImageDataWriter<GV,DC>
out << (i == 0 ? "" : " ") << spacing[i];
}
out << "\">\n";
out << "<Piece Extent=\"";
for (int i = 0; i < 3; ++i) {
out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
}
out << "\">\n";
dataCollector_.writeLocalPiece([&out](std::array<int,6> const& extent)
{
out << "<Piece Extent=\"";
for (int i = 0; i < 3; ++i) {
out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
}
out << "\">\n";
});
{ // Write data associated with grid points
auto scalar = std::find_if(pointData_.begin(), pointData_.end(), [](auto const& v) { return v.ncomps() == 1; });
......@@ -125,9 +127,9 @@ void VtkImageDataWriter<GV,DC>
<< "byte_order=\"" << this->getEndian() << "\" header_type=\"UInt64\""
<< (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\">\n" : ">\n");
out << "<PImageData GhostLevel=\"0\" WholeExtent=\"";
auto const& extent = dataCollector_.extent();
auto const& wholeExtent = dataCollector_.wholeExtent();
for (int i = 0; i < 3; ++i) {
out << (i == 0 ? "" : " ") << extent[2*i] << " " << extent[2*i+1];
out << (i == 0 ? "" : " ") << wholeExtent[2*i] << " " << wholeExtent[2*i+1];
}