Commit 34d21cab authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

cleanup of source code, removed several warnings

parent 4803f702
......@@ -63,7 +63,7 @@ public:
std::int64_t old_o = 0;
for (auto const& c : elements(gridView_, Partitions::all)) {
Vtk::CellType cellType(c.type());
for (int j = 0; j < c.subEntities(dim); ++j)
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.offsets.push_back(old_o += c.subEntities(dim));
cells.types.push_back(cellType.type());
......@@ -83,7 +83,7 @@ public:
localFct.bind(e);
Vtk::CellType cellType{e.type()};
auto refElem = referenceElement(e.geometry());
for (int j = 0; j < e.subEntities(dim); ++j) {
for (unsigned 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)));
......
......@@ -32,7 +32,7 @@ public:
auto const& indexSet = gridView_.indexSet();
for (auto const& c : elements(gridView_, Partitions::interior)) {
numPoints_ += c.subEntities(dim);
for (int i = 0; i < c.subEntities(dim); ++i)
for (unsigned int i = 0; i < c.subEntities(dim); ++i)
indexMap_[indexSet.subIndex(c, i, dim)] = vertex_idx++;
}
}
......@@ -50,7 +50,7 @@ public:
std::vector<T> data(this->numPoints() * 3);
auto const& indexSet = gridView_.indexSet();
for (auto const& element : elements(gridView_, Partitions::interior)) {
for (int i = 0; i < element.subEntities(dim); ++i) {
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);
for (std::size_t j = 0; j < v.size(); ++j)
......@@ -74,7 +74,7 @@ public:
auto const& indexSet = gridView_.indexSet();
for (auto const& c : elements(gridView_, Partitions::interior)) {
Vtk::CellType cellType(c.type());
for (int j = 0; j < c.subEntities(dim); ++j) {
for (unsigned 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);
}
......@@ -96,7 +96,7 @@ public:
localFct.bind(e);
Vtk::CellType cellType{e.type()};
auto refElem = referenceElement(e.geometry());
for (int j = 0; j < e.subEntities(dim); ++j) {
for (unsigned 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)));
......
......@@ -42,7 +42,7 @@ public:
auto refElem = referenceElement<T,dim>(element.type());
// vertices
for (int i = 0; i < element.subEntities(dim); ++i) {
for (unsigned 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)
......@@ -51,7 +51,7 @@ public:
data[idx + j] = T(0);
}
// edge centers
for (int i = 0; i < element.subEntities(dim-1); ++i) {
for (unsigned 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)
......@@ -79,12 +79,12 @@ public:
auto const& indexSet = gridView_.indexSet();
for (auto const& c : elements(gridView_, Partitions::interior)) {
Vtk::CellType cellType(c.type(), Vtk::QUADRATIC);
for (int j = 0; j < c.subEntities(dim); ++j) {
for (unsigned 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) {
for (unsigned 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);
......@@ -107,13 +107,13 @@ public:
localFct.bind(e);
Vtk::CellType cellType{e.type(), Vtk::QUADRATIC};
auto refElem = referenceElement(e.geometry());
for (int j = 0; j < e.subEntities(dim); ++j) {
for (unsigned 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) {
for (unsigned 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)
......
......@@ -25,6 +25,7 @@ public:
SPDataCollector (GridView const& gridView)
: Super(gridView)
, wholeExtent_(filledArray<6,int>(0))
, extent_(filledArray<6,int>(0))
, origin_(0.0)
, spacing_(0.0)
{}
......@@ -34,6 +35,11 @@ public:
return wholeExtent_;
}
std::array<int, 6> const& extentImpl () const
{
return extent_;
}
auto const& originImpl () const
{
return origin_;
......@@ -48,6 +54,7 @@ public:
{
Super::updateImpl();
auto const& localMesh = gridView_.impl().gridLevel().localMesh();
auto const& begin = gridView_.impl().gridLevel().globalMesh().begin();
auto const& end = gridView_.impl().gridLevel().globalMesh().end();
auto const& cube = gridView_.grid().domain().cube();
......@@ -55,40 +62,13 @@ public:
for (int i = 0; i < dim; ++i) {
wholeExtent_[2*i] = begin[i];
wholeExtent_[2*i+1] = end[i];
extent_[2*i] = localMesh.begin()[i];
extent_[2*i+1] = localMesh.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_;
......@@ -99,8 +79,8 @@ private:
namespace Impl
{
template<class GridView, class ct, int dim, template< int > class Ref, class Comm>
struct StructuredDataCollector<GridView, SPGrid<ct,dim,Ref,Comm>>
template <class GridView, class ct, int dim, template< int > class Ref, class Comm>
struct StructuredDataCollectorImpl<GridView, SPGrid<ct,dim,Ref,Comm>>
{
using type = SPDataCollector<GridView>;
};
......
......@@ -9,11 +9,11 @@ namespace Impl
{
// Should be specialized for concrete structured grid
template <class GridView, class Grid>
struct StructuredDataCollector;
struct StructuredDataCollectorImpl;
}
template <class GridView>
using StructuredDataCollector = typename Impl::StructuredDataCollector<GridView, typename GridView::Grid>::type;
using StructuredDataCollector = typename Impl::StructuredDataCollectorImpl<GridView, typename GridView::Grid>::type;
/// The Interface for structured data-collectors
......@@ -33,67 +33,94 @@ public:
, ghostLevel_(gridView.overlapSize(0))
{}
/// Return number of grid vertices
std::uint64_t numPointsImpl () const
{
return gridView_.size(GridView::dimension);
}
void updateImpl ()
/// Sequence of Index pairs [begin, end) for the cells in each direction
std::array<int, 6> wholeExtent () const
{
defaultDataCollector_.update();
return this->asDerived().wholeExtentImpl();
}
std::array<int, 6> const& wholeExtent () const
/// Sequence of Index pairs [begin, end) for the cells in each direction of the local partition
std::array<int, 6> extent () const
{
return this->asDerived().wholeExtentImpl();
return this->asDerived().extentImpl();
}
FieldVector<ctype, 3> const& origin () const
/// Lower left corner of the grid
FieldVector<ctype, 3> origin () const
{
return this->asDerived().originImpl();
}
FieldVector<ctype, 3> const& spacing () const
/// Constant grid spacing in each coordinate direction
FieldVector<ctype, 3> spacing () const
{
return this->asDerived().spacingImpl();
}
/// Call the `writer` with extent
template <class Writer>
void writeLocalPiece (Writer const& writer) const
{
this->asDerived().writeLocalPieceImpl(writer);
}
/// Call the `writer` with piece number and piece extent
template <class Writer>
void writePieces (Writer const& writer) const
{
this->asDerived().writePiecesImpl(writer);
}
/// Return the number of overlapping elements
int ghostLevel () const
{
return this->asDerived().ghostLevelImpl();
}
/// Return the coordinates along the ordinates x, y, and z
/// The coordinates defines point coordinates for an extent by specifying the ordinate along each axis.
template <class T>
std::array<std::vector<T>, 3> coordinates () const
{
return this->asDerived().template coordinatesImpl<T>();
}
public:
/// Return number of grid vertices
std::uint64_t numPointsImpl () const
{
return gridView_.size(GridView::dimension);
}
/// \copyref DefaultDataCollector::update.
void updateImpl ()
{
defaultDataCollector_.update();
/// Return the coordinates of all grid vertices in the order given by the indexSet
#if HAVE_MPI
int rank = -1;
int num_ranks = -1;
MPI_Comm_rank(gridView_.comm(), &rank);
MPI_Comm_size(gridView_.comm(), &num_ranks);
if (rank == 0) {
extents_.resize(num_ranks);
requests_.resize(num_ranks, MPI_REQUEST_NULL);
for (int i = 1; i < num_ranks; ++i)
MPI_Irecv(extents_[i].data(), extents_[i].size(), MPI_INT, i, /*tag=*/6, gridView_.comm(), &requests_[i]);
}
sendRequest_ = MPI_REQUEST_NULL;
#endif
}
/// \copydoc DefaultDataCollector::points.
template <class T>
std::vector<T> pointsImpl () const
{
return defaultDataCollector_.template points<T>();
}
/// Evaluate the `fct` at the corners of the elements
/// \copydoc DefaultDataCollector::pointData
template <class T, class GlobalFunction>
std::vector<T> pointDataImpl (GlobalFunction const& fct) const
{
......@@ -101,30 +128,79 @@ public:
}
private:
/// Default implementation for \ref writeLocalPiece. Calculates the extent and communicates it to
/// rank 0.
template <class Writer>
void writeLocalPieceImpl (Writer const& writer) const
{
auto&& extent = this->extent();
#if HAVE_MPI
int sendFlag = 0;
MPI_Status sendStatus;
MPI_Test(&sendRequest_, &sendFlag, &sendStatus);
if (sendFlag) {
int rank = -1;
MPI_Comm_rank(gridView_.comm(), &rank);
if (rank != 0) {
MPI_Isend(extent.data(), extent.size(), MPI_INT, 0, /*tag=*/6, gridView_.comm(), &sendRequest_);
} else {
extents_[0] = extent;
}
}
#endif
writer(extent);
}
/// Receive extent from all ranks and call the `writer` with the rank's extent vector
template <class Writer>
void writePiecesImpl (Writer const& writer) const
{
writer(0, extents_[0], true);
#if HAVE_MPI
int num_ranks = -1;
MPI_Comm_size(gridView_.comm(), &num_ranks);
for (int p = 1; p < num_ranks; ++p) {
int idx = -1;
MPI_Status status;
MPI_Waitany(num_ranks, requests_.data(), &idx, &status);
if (idx != MPI_UNDEFINED) {
assert(idx == status.MPI_SOURCE && status.MPI_TAG == 6);
writer(idx, extents_[idx], true);
}
}
#endif
}
/// Return the \ref GridView::overlapSize
int ghostLevelImpl () const
{
return ghostLevel_;
}
/// Ordinate along each axis with constant \ref spacing from the \ref origin
template <class T>
std::array<std::vector<T>, 3> coordinatesImpl () const
{
auto origin = this->origin();
auto spacing = this->spacing();
auto extent = this->extent();
std::array<std::vector<T>, 3> ordinates{};
writeLocalPiece([&ordinates,&origin,&spacing](auto const& extent) {
for (std::size_t d = 0; d < GridView::dimension; ++d) {
auto s = extent[2*d+1] - extent[2*d] + 1;
ordinates[d].resize(s);
for (std::size_t i = 0; i < s; ++i)
ordinates[d][i] = origin[d] + (extent[2*d] + i)*spacing[d];
}
});
for (std::size_t d = GridView::dimension; d < 3; ++d)
for (int d = 0; d < GridView::dimension; ++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)
ordinates[d].resize(1, T(0));
return ordinates;
......@@ -133,6 +209,12 @@ private:
private:
DefaultDataCollector<GridView> defaultDataCollector_;
int ghostLevel_;
#if HAVE_MPI
mutable std::vector<std::array<int,6>> extents_;
mutable std::vector<MPI_Request> requests_;
mutable MPI_Request sendRequest_ = MPI_REQUEST_NULL;
#endif
};
}} // end namespace Dune::experimental
......@@ -21,6 +21,7 @@ public:
YaspDataCollector (GridView const& gridView)
: Super(gridView)
, wholeExtent_(filledArray<6,int>(0))
, extent_(filledArray<6,int>(0))
, origin_(0.0)
, spacing_(0.0)
, level_(gridView.template begin<0,All_Partition>()->level())
......@@ -31,6 +32,11 @@ public:
return wholeExtent_;
}
std::array<int, 6> const& extentImpl () const
{
return extent_;
}
auto const& originImpl () const
{
return origin_;
......@@ -44,34 +50,24 @@ public:
void updateImpl ()
{
Super::updateImpl();
localPieceCalled_ = false;
for (int i = 0; i < dim; ++i) {
wholeExtent_[2*i] = 0;
wholeExtent_[2*i+1] = gridView_.grid().levelSize(level_,i);
}
auto const& gl = *gridView_.grid().begin(level_);
auto const& g = gl.interior[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;
}
auto it = gridView_.grid().begin(level_);
initGeometry(it->coords);
#if HAVE_MPI
int rank = -1;
int num_ranks = -1;
MPI_Comm_rank(gridView_.comm(), &rank);
MPI_Comm_size(gridView_.comm(), &num_ranks);
if (rank == 0) {
extents_.resize(num_ranks);
requests_.resize(num_ranks);
for (int i = 1; i < num_ranks; ++i)
MPI_Irecv(extents_[i].data(), extents_[i].size(), MPI_INT, i, /*tag=*/6, gridView_.comm(), &requests_[i]);
}
#endif
}
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) {
......@@ -96,51 +92,8 @@ public:
}
}
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;
}
#if HAVE_MPI
if (!localPieceCalled_) {
int rank = -1;
MPI_Comm_rank(gridView_.comm(), &rank);
if (rank != 0) {
MPI_Isend(extent.data(), extent.size(), MPI_INT, 0, /*tag=*/6, gridView_.comm(), &sendRequest_);
} else {
extents_[0] = extent;
}
localPieceCalled_ = true;
}
#endif
writer(extent);
}
template <class Writer>
void writePiecesImpl (Writer const& writer) const
{
assert(localPieceCalled_);
#if HAVE_MPI
int num_ranks = -1;
MPI_Comm_size(gridView_.comm(), &num_ranks);
std::vector<MPI_Status> status(num_ranks-1);
MPI_Waitall(num_ranks-1, requests_.data()+1, status.data());
for (int p = 0; p < num_ranks; ++p) {
writer(p, extents_[p], true);
}
#endif
}
/// Extract the ordinates from the Coordinates object of the current level
template <class T>
std::array<std::vector<T>, 3> coordinatesImpl () const
{
......@@ -148,17 +101,14 @@ public:
auto const& coords = it->coords;
std::array<std::vector<T>, 3> ordinates{};
writeLocalPieceImpl([&ordinates,&coords](auto const& extent)
{
for (std::size_t d = 0; d < dim; ++d) {
auto s = extent[2*d+1] - extent[2*d] + 1;
ordinates[d].resize(s);
for (std::size_t i = 0; i < s; ++i)
ordinates[d][i] = coords.coordinate(d, extent[2*d] + i);
}
});
for (std::size_t d = dim; d < 3; ++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] = coords.coordinate(d, extent_[2*d] + i);
}
for (int d = dim; d < 3; ++d)
ordinates[d].resize(1, T(0));
return ordinates;
......@@ -166,21 +116,17 @@ public:
private:
std::array<int, 6> wholeExtent_;
FieldVector<ctype,3> spacing_;
std::array<int, 6> extent_;
FieldVector<ctype,3> origin_;
FieldVector<ctype,3> spacing_;
std::vector<std::size_t> indexMap_;
int level_;
mutable std::vector<std::array<int,6>> extents_;
mutable std::vector<MPI_Request> requests_;
mutable MPI_Request sendRequest_;
mutable bool localPieceCalled_ = false;
};
namespace Impl
{
template<class GridView, int dim, class Coordinates>
struct StructuredDataCollector<GridView, YaspGrid<dim,Coordinates>>
template <class GridView, int dim, class Coordinates>
struct StructuredDataCollectorImpl<GridView, YaspGrid<dim,Coordinates>>
{
using type = YaspDataCollector<GridView>;
};
......
......@@ -32,10 +32,10 @@ namespace Dune { namespace experimental
Vtk::CellType cellType{type};
auto refElem = referenceElement<double,Grid::dimension>(type);
std::size_t nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]);