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

redesign of readers and grid creators, parallel traversal restricted to interor+border

parent 938c9f1e
*.vtu
*.pvtu
*.vti
*.pvti
*.vtr
*.pvtr
*.vtp
*.pvtp
*.vts
*.pvts
build*/
*.log
*.tmp
\ No newline at end of file
......@@ -5,10 +5,21 @@
namespace Dune {
template <class GridView, class Derived>
template <class GridView, class Derived, class Partition>
class DataCollectorInterface
{
public:
/// The partitionset to collect data from
static constexpr auto partition = Partition{};
/// The dimension of the grid
enum { dim = GridView::dimension };
/// The dimension of the world
enum { dow = GridView::dimensionworld };
public:
/// Store a copy of the GridView
DataCollectorInterface (GridView const& gridView)
: gridView_(gridView)
{}
......@@ -25,7 +36,13 @@ public:
return asDerived().ghostLevelImpl();
}
/// Return the number of points in the grid
/// \brief Return the number of cells in (this partition of the) grid
std::uint64_t numCells () const
{
return asDerived().numCellsImpl();
}
/// Return the number of points in (this partition of the) grid
std::uint64_t numPoints () const
{
return asDerived().numPointsImpl();
......@@ -59,8 +76,9 @@ public:
return asDerived().template pointDataImpl<T>(fct);
}
/// Return a flat vector of function values evaluated at the cells. \see pointData.
/// Return a flat vector of function values evaluated at the cells in the order of traversal.
/**
* \see pointData.
* Note: Cells might be descibred explicitly by connectivity, offsets, and types, e.g. in
* an UnstructuredGrid, or might be described implicitly by the grid type, e.g. in
* StructuredGrid.
......@@ -95,7 +113,7 @@ public: // default implementations
return gridView_.overlapSize(0);
}
// Evaluate `fct` in center of cell
// Evaluate `fct` in center of cell.
template <class T, class VtkFunction>
std::vector<T> cellDataImpl (VtkFunction const& fct) const;
......
......@@ -6,20 +6,20 @@
namespace Dune {
template <class GridView, class Derived>
template <class GV, class D, class P>
template <class T, class VtkFunction>
std::vector<T> DataCollectorInterface<GridView, Derived>
std::vector<T> DataCollectorInterface<GV,D,P>
::cellDataImpl (VtkFunction const& fct) const
{
std::vector<T> data(gridView_.size(0) * fct.ncomps());
MultipleCodimMultipleGeomTypeMapper<GridView> mapper(gridView_, mcmgElementLayout());
std::vector<T> data;
data.reserve(this->numCells() * fct.ncomps());
auto localFct = localFunction(fct);
for (auto const& e : elements(gridView_, Partitions::all)) {
for (auto const& e : elements(gridView_, partition)) {
localFct.bind(e);
auto refElem = referenceElement<T,GridView::dimension>(e.type());
std::size_t idx = fct.ncomps() * mapper.index(e);
auto refElem = referenceElement<T,dim>(e.type());
for (int comp = 0; comp < fct.ncomps(); ++comp)
data[idx + comp] = T(localFct.evaluate(comp, refElem.position(0,0)));
data.emplace_back(localFct.evaluate(comp, refElem.position(0,0)));
localFct.unbind();
}
return data;
......
......@@ -9,39 +9,57 @@ namespace Dune
{
/// Implementation of \ref DataCollector for linear cells, with continuous data.
template <class GridView>
template <class GridView, class Partition>
class ContinuousDataCollector
: public UnstructuredDataCollectorInterface<GridView, ContinuousDataCollector<GridView>>
: public UnstructuredDataCollectorInterface<GridView, ContinuousDataCollector<GridView,Partition>, Partition>
{
enum { dim = GridView::dimension };
using Self = ContinuousDataCollector;
using Super = UnstructuredDataCollectorInterface<GridView, Self>;
using Super = UnstructuredDataCollectorInterface<GridView, Self, Partition>;
public:
using Super::dim;
using Super::partition;
public:
ContinuousDataCollector (GridView const& gridView)
: Super(gridView)
{}
/// Collect the vertex indices
void updateImpl ()
{
numPoints_ = 0;
indexMap_.resize(gridView_.size(dim));
auto const& indexSet = gridView_.indexSet();
for (auto const& vertex : vertices(gridView_, partition))
indexMap_[indexSet.index(vertex)] = std::int64_t(numPoints_++);
if (gridView_.comm().size() > 1) {
auto&& e = elements(gridView_, partition);
numCells_ = std::distance(std::begin(e), std::end(e));
} else {
numCells_ = gridView_.size(0);
}
}
/// Return number of grid vertices
std::uint64_t numPointsImpl () const
{
return gridView_.size(dim);
return numPoints_;
}
/// 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(gridView_.size(dim) * 3);
auto const& indexSet = gridView_.indexSet();
for (auto const& vertex : vertices(gridView_, Partitions::all)) {
std::size_t idx = 3 * indexSet.index(vertex);
std::vector<T> data;
data.reserve(numPoints_ * 3);
for (auto const& vertex : vertices(gridView_, partition)) {
auto v = vertex.geometry().center();
for (std::size_t j = 0; j < v.size(); ++j)
data[idx + j] = T(v[j]);
data.emplace_back(v[j]);
for (std::size_t j = v.size(); j < 3u; ++j)
data[idx + j] = T(0);
data.emplace_back(0);
}
return data;
}
......@@ -49,11 +67,12 @@ public:
/// Return a vector of global unique ids of the points
std::vector<std::uint64_t> pointIdsImpl () const
{
std::vector<std::uint64_t> data(gridView_.size(dim));
std::vector<std::uint64_t> data;
data.reserve(numPoints_);
GlobalIndexSet<GridView> globalIndexSet(gridView_, dim);
auto const& indexSet = gridView_.indexSet();
for (auto const& vertex : vertices(gridView_, Partitions::all)) {
data[indexSet.index(vertex)] = std::uint64_t(globalIndexSet.index(vertex));
for (auto const& vertex : vertices(gridView_, partition)) {
data.emplace_back(globalIndexSet.index(vertex));
}
return data;
}
......@@ -61,7 +80,7 @@ public:
/// Return number of grid cells
std::uint64_t numCellsImpl () const
{
return gridView_.size(0);
return numCells_;
}
/// Return the types, offsets and connectivity of the cells, using the same connectivity as
......@@ -76,15 +95,15 @@ public:
});
Cells cells;
cells.connectivity.reserve(gridView_.size(0) * maxVertices);
cells.offsets.reserve(gridView_.size(0));
cells.types.reserve(gridView_.size(0));
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_, Partitions::all)) {
for (auto const& c : elements(gridView_, partition)) {
Vtk::CellType cellType(c.type());
for (unsigned int j = 0; j < c.subEntities(dim); ++j)
cells.connectivity.push_back( std::int64_t(indexSet.subIndex(c,cellType.permutation(j),dim)) );
cells.connectivity.emplace_back(indexMap_[indexSet.subIndex(c,cellType.permutation(j),dim)]);
cells.offsets.push_back(old_o += c.subEntities(dim));
cells.types.push_back(cellType.type());
}
......@@ -95,15 +114,15 @@ public:
template <class T, class GlobalFunction>
std::vector<T> pointDataImpl (GlobalFunction const& fct) const
{
std::vector<T> data(gridView_.size(dim) * fct.ncomps());
std::vector<T> data(numPoints_ * fct.ncomps());
auto const& indexSet = gridView_.indexSet();
auto localFct = localFunction(fct);
for (auto const& e : elements(gridView_, Partitions::all)) {
for (auto const& e : elements(gridView_, partition)) {
localFct.bind(e);
Vtk::CellType cellType{e.type()};
auto refElem = referenceElement(e.geometry());
for (unsigned int j = 0; j < e.subEntities(dim); ++j) {
std::size_t idx = fct.ncomps() * indexSet.subIndex(e,cellType.permutation(j),dim);
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)));
}
......@@ -114,6 +133,9 @@ public:
protected:
using Super::gridView_;
std::uint64_t numPoints_ = 0;
std::uint64_t numCells_ = 0;
std::vector<std::int64_t> indexMap_;
};
} // end namespace Dune
......@@ -6,14 +6,16 @@ namespace Dune
{
/// Implementation of \ref DataCollector for linear cells, with discontinuous data.
template <class GridView>
template <class GridView, class Partition>
class DiscontinuousDataCollector
: public UnstructuredDataCollectorInterface<GridView, DiscontinuousDataCollector<GridView>>
: public UnstructuredDataCollectorInterface<GridView, DiscontinuousDataCollector<GridView,Partition>, Partition>
{
enum { dim = GridView::dimension };
using Self = DiscontinuousDataCollector;
using Super = UnstructuredDataCollectorInterface<GridView, Self>;
using Super = UnstructuredDataCollectorInterface<GridView, Self, Partition>;
public:
using Super::dim;
using Super::partition;
public:
DiscontinuousDataCollector (GridView const& gridView)
......@@ -24,10 +26,12 @@ public:
void updateImpl ()
{
numPoints_ = 0;
numCells_ = 0;
indexMap_.resize(gridView_.size(dim));
std::int64_t vertex_idx = 0;
auto const& indexSet = gridView_.indexSet();
for (auto const& c : elements(gridView_, Partitions::interior)) {
for (auto const& c : elements(gridView_, partition)) {
numCells_++;
numPoints_ += c.subEntities(dim);
for (unsigned int i = 0; i < c.subEntities(dim); ++i)
indexMap_[indexSet.subIndex(c, i, dim)] = vertex_idx++;
......@@ -46,7 +50,7 @@ public:
{
std::vector<T> data(numPoints_ * 3);
auto const& indexSet = gridView_.indexSet();
for (auto const& element : elements(gridView_, Partitions::interior)) {
for (auto const& element : elements(gridView_, partition)) {
for (unsigned 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);
......@@ -62,7 +66,7 @@ public:
/// Return number of grid cells
std::uint64_t numCellsImpl () const
{
return gridView_.size(0);
return numCells_;
}
/// Connect the corners of each cell. The leads to a global discontinuous grid
......@@ -70,12 +74,12 @@ public:
{
Cells cells;
cells.connectivity.reserve(numPoints_);
cells.offsets.reserve(gridView_.size(0));
cells.types.reserve(gridView_.size(0));
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_, Partitions::interior)) {
for (auto const& c : elements(gridView_, partition)) {
Vtk::CellType cellType(c.type());
for (unsigned int j = 0; j < c.subEntities(dim); ++j) {
std::int64_t vertex_idx = indexMap_[indexSet.subIndex(c,cellType.permutation(j),dim)];
......@@ -95,7 +99,7 @@ public:
std::vector<T> data(numPoints_ * fct.ncomps());
auto const& indexSet = gridView_.indexSet();
auto localFct = localFunction(fct);
for (auto const& e : elements(gridView_, Partitions::interior)) {
for (auto const& e : elements(gridView_, partition)) {
localFct.bind(e);
Vtk::CellType cellType{e.type()};
auto refElem = referenceElement(e.geometry());
......@@ -111,6 +115,7 @@ public:
protected:
using Super::gridView_;
std::uint64_t numCells_ = 0;
std::uint64_t numPoints_ = 0;
std::vector<std::int64_t> indexMap_;
};
......
......@@ -8,12 +8,14 @@ namespace Dune
/// Implementation of \ref DataCollector for quadratic cells, with continuous data.
template <class GridView>
class QuadraticDataCollector
: public UnstructuredDataCollectorInterface<GridView, QuadraticDataCollector<GridView>>
: public UnstructuredDataCollectorInterface<GridView, QuadraticDataCollector<GridView>, Partitions::All>
{
enum { dim = GridView::dimension };
using Self = QuadraticDataCollector;
using Super = UnstructuredDataCollectorInterface<GridView, Self>;
using Super = UnstructuredDataCollectorInterface<GridView, Self, Partitions::All>;
public:
using Super::dim;
using Super::partition; // NOTE: quadratic data-collector currently implemented for the All partition only
public:
QuadraticDataCollector (GridView const& gridView)
......@@ -36,7 +38,7 @@ public:
{
std::vector<T> data(this->numPoints() * 3);
auto const& indexSet = gridView_.indexSet();
for (auto const& element : elements(gridView_, Partitions::interior)) {
for (auto const& element : elements(gridView_, partition)) {
auto geometry = element.geometry();
auto refElem = referenceElement<T,dim>(element.type());
......@@ -82,7 +84,7 @@ public:
std::int64_t old_o = 0;
auto const& indexSet = gridView_.indexSet();
for (auto const& c : elements(gridView_, Partitions::interior)) {
for (auto const& c : elements(gridView_, partition)) {
Vtk::CellType cellType(c.type(), Vtk::QUADRATIC);
for (unsigned int j = 0; j < c.subEntities(dim); ++j) {
int k = cellType.permutation(j);
......@@ -107,7 +109,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_, Partitions::interior)) {
for (auto const& e : elements(gridView_, partition)) {
localFct.bind(e);
Vtk::CellType cellType{e.type(), Vtk::QUADRATIC};
auto refElem = referenceElement(e.geometry());
......
......@@ -15,12 +15,14 @@ template <class GridView>
class SPDataCollector
: public StructuredDataCollectorInterface<GridView, SPDataCollector<GridView>>
{
enum { dim = GridView::dimension };
using Self = SPDataCollector;
using Super = StructuredDataCollectorInterface<GridView, Self>;
using ctype = typename GridView::ctype;
public:
using Super::dim;
using Super::partition;
public:
SPDataCollector (GridView const& gridView)
: Super(gridView)
......
......@@ -19,6 +19,10 @@ protected:
using SubDataCollector = ContinuousDataCollector<GridView>;
using ctype = typename GridView::ctype;
public:
using Super::dim;
using Super::partition;
public:
StructuredDataCollectorInterface (GridView const& gridView)
: Super(gridView)
......@@ -103,6 +107,16 @@ public: // default implementation:
#endif
}
/// Return number of grid cells
std::uint64_t numCellsImpl () const
{
auto extent = this->extent();
std::uint64_t num = 1;
for (int d = 0; d < dim; ++d)
num *= extent[2*d+1] - extent[2*d];
return num;
}
/// Return number of grid vertices
std::uint64_t numPointsImpl () const
{
......@@ -192,14 +206,14 @@ public: // default implementation:
auto extent = this->extent();
std::array<std::vector<T>, 3> ordinates{};
for (int d = 0; d < GridView::dimension; ++d) {
for (int d = 0; d < dim; ++d) {
auto s = extent[2*d+1] - extent[2*d] + 1;
ordinates[d].resize(s);
for (int i = 0; i < s; ++i)
ordinates[d][i] = origin[d] + (extent[2*d] + i)*spacing[d];
}
for (int d = GridView::dimension; d < 3; ++d)
for (int d = dim; d < 3; ++d)
ordinates[d].resize(1, T(0));
return ordinates;
......
......@@ -14,23 +14,21 @@ struct Cells
std::vector<std::int64_t> connectivity;
};
template <class GridView, class Derived>
template <class GridView, class Derived, class Partition>
class UnstructuredDataCollectorInterface
: public DataCollectorInterface<GridView, Derived>
: public DataCollectorInterface<GridView, Derived, Partition>
{
using Super = DataCollectorInterface<GridView, Derived>;
using Super = DataCollectorInterface<GridView, Derived, Partition>;
public:
using Super::dim;
using Super::partition;
public:
UnstructuredDataCollectorInterface (GridView const& gridView)
: Super(gridView)
{}
/// \brief Return the number of cells in the grid
std::uint64_t numCells () const
{
return this->asDerived().numCellsImpl();
}
/// \brief Return cell types, offsets, and connectivity. \see Cells
Cells cells () const
{
......
......@@ -11,12 +11,14 @@ template <class GridView>
class YaspDataCollector
: public StructuredDataCollectorInterface<GridView, YaspDataCollector<GridView>>
{
enum { dim = GridView::dimension };
using Self = YaspDataCollector;
using Super = StructuredDataCollectorInterface<GridView, Self>;
using ctype = typename GridView::ctype;
public:
using Super::dim;
using Super::partition;
public:
YaspDataCollector (GridView const& gridView)
: Super(gridView)
......
......@@ -4,18 +4,18 @@ namespace Dune
{
// forward declaration of all classes in dune-vtk
template <class GridView, class Derived>
template <class GridView, class Derived, class Partition = Partitions::InteriorBorder>
class DataCollectorInterface;
// @{ datacollectors
template <class GridView, class Derived>
template <class GridView, class Derived, class Partition = Partitions::InteriorBorder>
class UnstructuredDataCollectorInterface;
// @{ unstructured-datacollectors
template <class GridView>
template <class GridView, class Partition = Partitions::InteriorBorder>
class ContinuousDataCollector;
template <class GridView>
template <class GridView, class Partition = Partitions::InteriorBorder>
class DiscontinuousDataCollector;
template <class GridView>
......
......@@ -21,18 +21,7 @@ public:
public:
GridCreatorInterface (GridFactory<Grid>& factory)
: factory_(&factory)
{
#if DUNE_VERSION_LT(DUNE_GRID,2,7)
// old GridFactory implementation does not provide access to its collective communication
#if HAVE_MPI
MPI_Comm_rank(MPI_COMM_WORLD, &rank_);
MPI_Comm_size(MPI_COMM_WORLD, &numRanks_);
#endif
#else
rank_ = factory.comm().rank();
numRanks_ = factory.comm().size();
#endif
}
{}
/// Insert all points as vertices into the factory
void insertVertices (std::vector<GlobalCoordinate> const& points,
......@@ -61,16 +50,10 @@ public:
return *factory_;
}
/// Return the MPI_Comm_rank of the factory, (or of MPI_COMM_WORLD)
int rank () const
{
return rank_;
}
/// Return the MPI_Comm_size of the factory, (or of MPI_COMM_WORLD)
int size () const
/// Return the mpi collective communicator
auto comm () const
{
return numRanks_;
return MPIHelper::getCollectiveCommunication();
}
protected: // cast to derived type
......@@ -107,8 +90,6 @@ public: // default implementations
protected:
GridFactory<Grid>* factory_;
int rank_ = 0;
int numRanks_ = 1;
};
} // end namespace Dune
......@@ -76,14 +76,13 @@ namespace Dune
vtk_cell.push_back(new_idx);
}
if (cellType.noPermutation())
if (cellType.noPermutation()) {
factory().insertElement(type,vtk_cell);
else {
} else {
// apply index permutation
std::vector<unsigned int> cell(nNodes);
for (int j = 0; j < nNodes; ++j)
cell[j] = vtk_cell[cellType.permutation(j)];
factory().insertElement(type,cell);
}
}
......
......@@ -29,20 +29,19 @@ namespace Dune
: Super(factory)
{}
using Super::factory;
void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
std::vector<std::uint64_t> const& point_ids)
{
assert(point_ids.size() == points.size());
for (std::size_t i = 0; i < points.