diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..0558040d18d468008693772154cbf110fcc710cc --- /dev/null +++ b/.gitignore @@ -0,0 +1,15 @@ +*.vtu +*.pvtu +*.vti +*.pvti +*.vtr +*.pvtr +*.vtp +*.pvtp +*.vts +*.pvts + +build*/ + +*.log +*.tmp \ No newline at end of file diff --git a/dune/vtk/datacollectorinterface.hh b/dune/vtk/datacollectorinterface.hh index 79f53924ec1b32c560654af98b37758a9452cecb..f2b896812f96263709c376fff7e2d06dfb633250 100644 --- a/dune/vtk/datacollectorinterface.hh +++ b/dune/vtk/datacollectorinterface.hh @@ -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; diff --git a/dune/vtk/datacollectorinterface.impl.hh b/dune/vtk/datacollectorinterface.impl.hh index abf5e6ad9fbc3fbd3a67f43f316689a9741ad940..eedcfd35154eb6e5380fda8abb0fb9dc0bd683b6 100644 --- a/dune/vtk/datacollectorinterface.impl.hh +++ b/dune/vtk/datacollectorinterface.impl.hh @@ -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; diff --git a/dune/vtk/datacollectors/continuousdatacollector.hh b/dune/vtk/datacollectors/continuousdatacollector.hh index 118dfb958d36f539cccf65eb0d58b9406ad8cba4..62bc577316697b0b238aa0696f79bacf3b684859 100644 --- a/dune/vtk/datacollectors/continuousdatacollector.hh +++ b/dune/vtk/datacollectors/continuousdatacollector.hh @@ -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 diff --git a/dune/vtk/datacollectors/discontinuousdatacollector.hh b/dune/vtk/datacollectors/discontinuousdatacollector.hh index 4a50662f81abc9ec423ee7472f84e3f78a4420ac..25089a77b33e595a5afc09df8ec54ffb664b95b7 100644 --- a/dune/vtk/datacollectors/discontinuousdatacollector.hh +++ b/dune/vtk/datacollectors/discontinuousdatacollector.hh @@ -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_; }; diff --git a/dune/vtk/datacollectors/quadraticdatacollector.hh b/dune/vtk/datacollectors/quadraticdatacollector.hh index 55e98c24d9576e73a8f95ce57308cffc6167148a..62a1497efc340bd5a8fc92cdc68c7e638829807f 100644 --- a/dune/vtk/datacollectors/quadraticdatacollector.hh +++ b/dune/vtk/datacollectors/quadraticdatacollector.hh @@ -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()); diff --git a/dune/vtk/datacollectors/spdatacollector.hh b/dune/vtk/datacollectors/spdatacollector.hh index 3d416d33d9f35b0363adf56c38e46d98391f9201..2af97ec09badbe325ce110c61d0a132dce69e51a 100644 --- a/dune/vtk/datacollectors/spdatacollector.hh +++ b/dune/vtk/datacollectors/spdatacollector.hh @@ -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) diff --git a/dune/vtk/datacollectors/structureddatacollector.hh b/dune/vtk/datacollectors/structureddatacollector.hh index c7c1167c4ae5f00d3d1e51ac03701080984475bf..05a9f07f9370bede3cb37a2a622ed8161c29052f 100644 --- a/dune/vtk/datacollectors/structureddatacollector.hh +++ b/dune/vtk/datacollectors/structureddatacollector.hh @@ -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; diff --git a/dune/vtk/datacollectors/unstructureddatacollector.hh b/dune/vtk/datacollectors/unstructureddatacollector.hh index 9a964dbb7ea4d850198a2e0366c031cf6ab73314..802cf7c5c107f9bcf35d21c381d5bf119f71b62d 100644 --- a/dune/vtk/datacollectors/unstructureddatacollector.hh +++ b/dune/vtk/datacollectors/unstructureddatacollector.hh @@ -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 { diff --git a/dune/vtk/datacollectors/yaspdatacollector.hh b/dune/vtk/datacollectors/yaspdatacollector.hh index 12315d175e602f38023ce5afa8a672b65f459f06..fbfe881d07181081871e80613313e553893d2757 100644 --- a/dune/vtk/datacollectors/yaspdatacollector.hh +++ b/dune/vtk/datacollectors/yaspdatacollector.hh @@ -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) diff --git a/dune/vtk/forward.hh b/dune/vtk/forward.hh index efc8538d98d13c56d479b458c0b6266074f6d44e..3ade39f6b098dad6e5c9df4f7363571b6a2dc2d7 100644 --- a/dune/vtk/forward.hh +++ b/dune/vtk/forward.hh @@ -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> diff --git a/dune/vtk/gridcreatorinterface.hh b/dune/vtk/gridcreatorinterface.hh index 634a248eabc11ecb698f0f2fdf3895454dde1529..56fb6f0352c1b1b3b04c936e682617bc46bea85a 100644 --- a/dune/vtk/gridcreatorinterface.hh +++ b/dune/vtk/gridcreatorinterface.hh @@ -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 diff --git a/dune/vtk/gridcreators/discontinuousgridcreator.hh b/dune/vtk/gridcreators/discontinuousgridcreator.hh index c392d4f3f8016d6c826c7f642f7b6382631c9ed0..5160b00e92a29521b9c42ac3814ce8b2901da63c 100644 --- a/dune/vtk/gridcreators/discontinuousgridcreator.hh +++ b/dune/vtk/gridcreators/discontinuousgridcreator.hh @@ -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); } } diff --git a/dune/vtk/gridcreators/parallelgridcreator.hh b/dune/vtk/gridcreators/parallelgridcreator.hh index 4b69554705c7bbf100944c9248dfd7ba956001ed..fc2a08c09bb7ca1d469597830bba5646c9a9ee1e 100644 --- a/dune/vtk/gridcreators/parallelgridcreator.hh +++ b/dune/vtk/gridcreators/parallelgridcreator.hh @@ -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.size(); ++i) - factory().insertVertex(points[i], VertexId(point_ids[i])); + this->factory().insertVertex(points[i], VertexId(point_ids[i])); } void insertPiecesImpl (std::vector<std::string> const& pieces) { - if (int(pieces.size()) == this->size()) { - VtkReader<Grid, Self> pieceReader(factory()); - pieceReader.readFromFile(pieces[this->rank()], true); + if (int(pieces.size()) == this->comm().size()) { + VtkReader<Grid, Self> pieceReader(this->factory()); + pieceReader.readFromFile(pieces[this->comm().rank()], true); } } }; diff --git a/dune/vtk/gridcreators/serialgridcreator.hh b/dune/vtk/gridcreators/serialgridcreator.hh index c19d33f70914bfbd694d39ab92cdd193c9f36bf8..75aeb13b76b5ccac3bc5ac94f36d39009e281c40 100644 --- a/dune/vtk/gridcreators/serialgridcreator.hh +++ b/dune/vtk/gridcreators/serialgridcreator.hh @@ -49,10 +49,12 @@ namespace Dune void insertPiecesImpl (std::vector<std::string> const& pieces) { - if (this->rank() == 0) { + if (this->comm().rank() == 0) { VtkReader<Grid, Self> pieceReader(*this); - for (std::string const& piece : pieces) - pieceReader.readFromFile(piece); + for (std::string const& piece : pieces) { + pieceReader.readFromFile(piece, false); + pieceReader.createGrid(false); + } DiscontinuousGridCreator<Grid> creator(this->factory()); creator.insertVertices(points_, {}); diff --git a/dune/vtk/pvdwriter.impl.hh b/dune/vtk/pvdwriter.impl.hh index 8841309fe1fab4cf32b65e6768d139f5873f141d..64d5ca6de3a8644df86f32dad93c8ce8d349cd32 100644 --- a/dune/vtk/pvdwriter.impl.hh +++ b/dune/vtk/pvdwriter.impl.hh @@ -25,15 +25,15 @@ void PvdWriter<W> std::string ext = "." + vtkWriter_.getFileExtension(); - int rank = vtkWriter_.rank_; - int numRanks = vtkWriter_.numRanks_; - if (numRanks > 1) + int commRank = vtkWriter_.comm().rank(); + int commSize = vtkWriter_.comm().size(); + if (commSize > 1) ext = ".p" + vtkWriter_.getFileExtension(); timesteps_.emplace_back(time, rel_fn + ext); vtkWriter_.write(seq_fn + ext); - if (rank == 0 && writeCollection) { + if (commRank == 0 && writeCollection) { std::ofstream out(pvd_fn + ".pvd", std::ios_base::ate | std::ios::binary); assert(out.is_open()); @@ -56,8 +56,8 @@ void PvdWriter<W> p.remove_filename(); p /= name.string(); - int rank = vtkWriter_.rank_; - if (rank == 0) { + int commRank = vtkWriter_.comm().rank(); + if (commRank == 0) { std::ofstream out(p.string() + ".pvd", std::ios_base::ate | std::ios::binary); assert(out.is_open()); diff --git a/dune/vtk/vtkreader.hh b/dune/vtk/vtkreader.hh index 82a2b536b8b205cb9b4e3746c34e098da2f02126..d49199dd3061f645ef83e7270294e57f155834ce 100644 --- a/dune/vtk/vtkreader.hh +++ b/dune/vtk/vtkreader.hh @@ -69,7 +69,7 @@ namespace Dune } /// Construct a grid using the GridCreator - void createGrid(); + void createGrid(bool insertPieces = true); /// Return the filenames of parallel pieces std::vector<std::string> const& pieces () const @@ -129,6 +129,11 @@ namespace Dune // clear all vectors void clear (); + auto comm () const + { + return MPIHelper::getCollectiveCommunication(); + } + private: std::unique_ptr<GridCreator> creatorStorage_ = nullptr; GridCreator& creator_; diff --git a/dune/vtk/vtkreader.impl.hh b/dune/vtk/vtkreader.impl.hh index 0f8099ac98ab9f488216408b737294573db46f10..3fd29ba1c1503ab5d882f9b302416ee44166a27d 100644 --- a/dune/vtk/vtkreader.impl.hh +++ b/dune/vtk/vtkreader.impl.hh @@ -28,8 +28,9 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename, bool cr std::string ext = filesystem::path(filename).extension().string(); if (ext == ".vtu") { readSerialFileFromStream(input, create); + pieces_.push_back(filename); } else if (ext == ".pvtu") { - readParallelFileFromStream(input, creator_.rank(), creator_.size(), create); + readParallelFileFromStream(input, comm().rank(), comm().size(), create); } else { DUNE_THROW(IOError, "File has unknown file-extension '" << ext << "'. Allowed are only '.vtu' and '.pvtu'."); } @@ -510,7 +511,7 @@ void VtkReader<Grid,Creator>::readAppended (std::ifstream& input, std::vector<T> template <class Grid, class Creator> -void VtkReader<Grid,Creator>::createGrid () +void VtkReader<Grid,Creator>::createGrid (bool insertPieces) { assert(vec_points.size() == numberOfPoints_); assert(vec_types.size() == numberOfCells_); @@ -520,7 +521,8 @@ void VtkReader<Grid,Creator>::createGrid () creator_.insertVertices(vec_points, vec_point_ids); if (!vec_types.empty()) creator_.insertElements(vec_types, vec_offsets, vec_connectivity); - creator_.insertPieces(pieces_); + if (insertPieces) + creator_.insertPieces(pieces_); } // Assume input already read the line <AppendedData ...> diff --git a/dune/vtk/vtktimeserieswriter.impl.hh b/dune/vtk/vtktimeserieswriter.impl.hh index e2cd4b579655394d3bc25fcc853ea6d8e42797e1..1255a4f560412bbbd5ad3363f3bcfdfd8b17efac 100644 --- a/dune/vtk/vtktimeserieswriter.impl.hh +++ b/dune/vtk/vtktimeserieswriter.impl.hh @@ -49,10 +49,8 @@ void VtkTimeseriesWriter<W> std::string filenameBase = tmp.string(); - int rank = vtkWriter_.rank_; - int numRanks = vtkWriter_.numRanks_; - if (numRanks > 1) - filenameBase = tmp.string() + "_p" + std::to_string(rank); + if (vtkWriter_.comm().size() > 1) + filenameBase = tmp.string() + "_p" + std::to_string(vtkWriter_.comm().rank()); if (!initialized_) { // write points and cells only once @@ -90,10 +88,10 @@ void VtkTimeseriesWriter<W> std::string parallel_fn = data_dir.string() + '/' + name.string() + "_ts"; std::string rel_fn = rel_dir.string() + '/' + name.string() + "_ts"; - int rank = vtkWriter_.rank_; - int numRanks = vtkWriter_.numRanks_; - if (numRanks > 1) - serial_fn += "_p" + std::to_string(rank); + int commRank = vtkWriter_.comm().rank(); + int commSize = vtkWriter_.comm().size(); + if (commSize > 1) + serial_fn += "_p" + std::to_string(commRank); { // write serial file std::ofstream serial_out(serial_fn + "." + vtkWriter_.getFileExtension(), @@ -108,7 +106,7 @@ void VtkTimeseriesWriter<W> vtkWriter_.writeTimeseriesSerialFile(serial_out, filenameMesh_, timesteps_, blocks_); } - if (numRanks > 1 && rank == 0) { + if (commSize > 1 && commRank == 0) { // write parallel file std::ofstream parallel_out(parallel_fn + ".p" + vtkWriter_.getFileExtension(), std::ios_base::ate | std::ios::binary); @@ -119,7 +117,7 @@ void VtkTimeseriesWriter<W> ? std::numeric_limits<float>::digits10+2 : std::numeric_limits<double>::digits10+2); - vtkWriter_.writeTimeseriesParallelFile(parallel_out, rel_fn, numRanks, timesteps_); + vtkWriter_.writeTimeseriesParallelFile(parallel_out, rel_fn, commSize, timesteps_); } } diff --git a/dune/vtk/vtkwriterinterface.hh b/dune/vtk/vtkwriterinterface.hh index 7f6ea2541157f74a935e30935a493d2e99734641..b753f415848d6e10b31f168f6da0c5a6124c6033 100644 --- a/dune/vtk/vtkwriterinterface.hh +++ b/dune/vtk/vtkwriterinterface.hh @@ -25,6 +25,7 @@ namespace Dune static constexpr int dimension = GridView::dimension; using VtkFunction = Dune::VtkFunction<GridView>; + using Communicator = CollectiveCommunication<typename MPIHelper::MPICommunicator>; using pos_type = typename std::ostream::pos_type; enum PositionTypes { @@ -38,8 +39,6 @@ namespace Dune Vtk::FormatTypes format = Vtk::BINARY, Vtk::DataTypes datatype = Vtk::FLOAT32) : dataCollector_(gridView) - , rank_(gridView.comm().rank()) - , numRanks_(gridView.comm().size()) , format_(format) , datatype_(datatype) { @@ -148,13 +147,14 @@ namespace Dune return datatype_; } + auto comm () const + { + return MPIHelper::getCollectiveCommunication(); + } + protected: mutable DataCollector dataCollector_; - // the rank and size of the gridView collective communication - int rank_ = 0; - int numRanks_ = 1; - Vtk::FormatTypes format_; Vtk::DataTypes datatype_; diff --git a/dune/vtk/vtkwriterinterface.impl.hh b/dune/vtk/vtkwriterinterface.impl.hh index 9e35b7af4e9a3bf421c23de462f7702d72553e58..ab08febbefc40218e6e76ec18bb24ca1de041786 100644 --- a/dune/vtk/vtkwriterinterface.impl.hh +++ b/dune/vtk/vtkwriterinterface.impl.hh @@ -39,8 +39,8 @@ void VtkWriterInterface<GV,DC> std::string parallel_fn = data_dir.string() + '/' + name.string(); std::string rel_fn = rel_dir.string() + '/' + name.string(); - if (numRanks_ > 1) - serial_fn += "_p" + std::to_string(rank_); + if (comm().size() > 1) + serial_fn += "_p" + std::to_string(comm().rank()); { // write serial file std::ofstream serial_out(serial_fn + "." + fileExtension(), std::ios_base::ate | std::ios::binary); @@ -54,7 +54,7 @@ void VtkWriterInterface<GV,DC> writeSerialFile(serial_out); } - if (numRanks_ > 1 && rank_ == 0) { + if (comm().size() > 1 && comm().rank() == 0) { // write parallel file std::ofstream parallel_out(parallel_fn + ".p" + fileExtension(), std::ios_base::ate | std::ios::binary); assert(parallel_out.is_open()); @@ -64,7 +64,7 @@ void VtkWriterInterface<GV,DC> ? std::numeric_limits<float>::digits10+2 : std::numeric_limits<double>::digits10+2); - writeParallelFile(parallel_out, rel_fn, numRanks_); + writeParallelFile(parallel_out, rel_fn, comm().size()); } } diff --git a/src/test/parallel_reader_writer_test.cc b/src/test/parallel_reader_writer_test.cc index cb126ba2c1e52e25bff8d4140166fd7d54f44ee9..856c96e53c6e82399ed4ade4e4d657bdc0625e4c 100644 --- a/src/test/parallel_reader_writer_test.cc +++ b/src/test/parallel_reader_writer_test.cc @@ -66,16 +66,13 @@ bool compare_files (std::string const& fn1, std::string const& fn2) return true; } -using TestCase = std::tuple<std::string,Vtk::FormatTypes,Vtk::DataTypes>; -using TestCases = std::set<TestCase>; -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 G> struct HasParallelGridFactory : std::false_type {}; +#if DUNE_VERSION_GT(DUNE_GRID,2,6) && HAVE_DUNE_ALUGRID +template<int dim, int dimworld, Dune::ALUGridElementType elType, Dune::ALUGridRefinementType refineType, class Comm> +struct HasParallelGridFactory<Dune::ALUGrid<dim,dimworld,elType,refineType,Comm>> : std::true_type {}; +#endif + template <class Test> void compare (Test& test, filesystem::path const& dir, filesystem::path const& name) @@ -85,71 +82,87 @@ void compare (Test& test, filesystem::path const& dir, filesystem::path const& n } template <class GridView> -void writer_test (GridView const& gridView) +void writer_test (GridView const& gridView, std::string base_name) { - for (auto const& test_case : test_cases) { - VtkUnstructuredGridWriter<GridView> vtkWriter(gridView, std::get<1>(test_case), std::get<2>(test_case)); - vtkWriter.write("parallel_reader_writer_test_" + std::get<0>(test_case) + ".vtu"); - } + VtkUnstructuredGridWriter<GridView> vtkWriter(gridView, Vtk::ASCII, Vtk::FLOAT32); + vtkWriter.write(base_name + ".vtu"); } -template <class G> struct IsALUGrid : std::false_type {}; -#if DUNE_VERSION_GT(DUNE_GRID,2,6) && HAVE_DUNE_ALUGRID -template<int dim, int dimworld, Dune::ALUGridElementType elType, Dune::ALUGridRefinementType refineType, class Comm> -struct IsALUGrid<Dune::ALUGrid<dim,dimworld,elType,refineType,Comm>> : std::true_type {}; -#endif - template <class Grid, class Creator> -void reader_test (MPIHelper& mpi, TestSuite& test) +void reader_writer_test(MPIHelper& mpi, TestSuite& test, std::string const& testName, bool doLoadBalance = true) { + std::cout << "== " << testName << "\n"; + std::string base_name = "parallel_rw_dim" + std::to_string(Grid::dimension); + std::vector<std::string> pieces1, pieces2; + std::string ext = ".vtu"; if (mpi.size() > 1) ext = ".pvtu"; - TestCase test_case = {"ascii32", Vtk::ASCII, Vtk::FLOAT32}; - - GridFactory<Grid> factory; - VtkReader<Grid, Creator> reader{factory}; - reader.readFromFile("parallel_reader_writer_test_" + std::get<0>(test_case) + ext); - - std::unique_ptr<Grid> grid = Hybrid::ifElse(IsALUGrid<Grid>{}, - [&](auto id) { return id(factory).createGrid(std::true_type{}); }, - [&](auto id) { return id(factory).createGrid(); }); - std::vector<std::string> pieces1 = grid->comm().size() > 1 ? - reader.pieces() : - std::vector<std::string>{"parallel_reader_writer_test_" + std::get<0>(test_case) + ".vtu"}; - - VtkUnstructuredGridWriter<typename Grid::LeafGridView> vtkWriter(grid->leafGridView(), - std::get<1>(test_case), std::get<2>(test_case)); - vtkWriter.write("parallel_reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"); - - GridFactory<Grid> factory2; - VtkReader<Grid, Creator> reader2{factory2}; - reader2.readFromFile("parallel_reader_writer_test_" + std::get<0>(test_case) + "_2" + ext, false); - std::vector<std::string> pieces2 = grid->comm().size() > 1 ? - reader.pieces() : - std::vector<std::string>{"parallel_reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"}; - - test.check(pieces1.size() == pieces2.size(), "pieces1.size == pieces2.size"); - for (std::size_t i = 0; i < pieces1.size(); ++i) - test.check(compare_files(pieces1[i], pieces2[i])); -} + // Step 1: create a new grid and write it to file1 + const int dim = Grid::dimension; + { + FieldVector<double,dim> lowerLeft; lowerLeft = 0.0; + FieldVector<double,dim> upperRight; upperRight = 1.0; + auto numElements = filledArray<dim,unsigned int>(4); + auto gridPtr = StructuredGridFactory<Grid>::createSimplexGrid(lowerLeft, upperRight, numElements); + gridPtr->loadBalance(); + // std::cout << "write1\n"; + writer_test(gridPtr->leafGridView(), base_name); + } + mpi.getCollectiveCommunication().barrier(); // need a barrier between write and read -template <class Grid, class Creator> -void reader_writer_test(MPIHelper& mpi, TestSuite& test, std::string const& testName) -{ - std::cout << "== " << testName << "\n"; - const int dim = Grid::dimension; - FieldVector<double,dim> lowerLeft; lowerLeft = 0.0; - FieldVector<double,dim> upperRight; upperRight = 1.0; - auto numElements = filledArray<dim,unsigned int>(4); - auto gridPtr = StructuredGridFactory<Grid>::createSimplexGrid(lowerLeft, upperRight, numElements); - gridPtr->loadBalance(); - - writer_test(gridPtr->leafGridView()); - MPI_Barrier(MPI_COMM_WORLD); - reader_test<Grid, Creator>(mpi,test); + // Step 2: read the grid from file1 and write it back to file2 + { GridFactory<Grid> factory; + // std::cout << "read1\n"; + VtkReader<Grid, Creator> reader{factory}; + reader.readFromFile(base_name + ext); + + std::unique_ptr<Grid> grid{ Hybrid::ifElse(HasParallelGridFactory<Grid>{}, + [&](auto id) { return id(factory).createGrid(std::true_type{}); }, + [&](auto id) { return id(factory).createGrid(); }) }; + if (doLoadBalance) + grid->loadBalance(); + + writer_test(grid->leafGridView(), base_name + "_1"); + } + + mpi.getCollectiveCommunication().barrier(); + + // Step 3: read the (parallel) file1 to get the piece filenames + { GridFactory<Grid> factory; + VtkReader<Grid, Creator> reader{factory}; + reader.readFromFile(base_name + "_1" + ext); + + std::unique_ptr<Grid> grid{ Hybrid::ifElse(HasParallelGridFactory<Grid>{}, + [&](auto id) { return id(factory).createGrid(std::true_type{}); }, + [&](auto id) { return id(factory).createGrid(); }) }; + if (doLoadBalance) + grid->loadBalance(); + pieces1 = reader.pieces(); + + writer_test(grid->leafGridView(), base_name + "_2"); + } + + mpi.getCollectiveCommunication().barrier(); + + // Step 4: read the (parallel) file2 to get the piece filenames + { GridFactory<Grid> factory; + VtkReader<Grid, Creator> reader{factory}; + reader.readFromFile(base_name + "_2" + ext, false); + + pieces2 = reader.pieces(); + } + + mpi.getCollectiveCommunication().barrier(); + + // Step 4: compare the pieces + if (mpi.rank() == 0) { + test.check(pieces1.size() == pieces2.size(), "pieces1.size == pieces2.size"); + for (std::size_t i = 0; i < pieces1.size(); ++i) + test.check(compare_files(pieces1[i], pieces2[i]), "compare(" + pieces1[i] + ", " + pieces2[i] + ")"); + } } #if HAVE_DUNE_ALUGRID @@ -175,10 +188,10 @@ int main (int argc, char** argv) reader_writer_test<ALUGridType<2>, SerialGridCreator<ALUGridType<2>>>(mpi, test, "ALUGridType<2>"); reader_writer_test<ALUGridType<3>, SerialGridCreator<ALUGridType<3>>>(mpi, test, "ALUGridType<3>"); -// #if DUNE_VERSION_LT(DUNE_GRID,2,7) - reader_writer_test<ALUGridType<2>, ParallelGridCreator<ALUGridType<2>>>(mpi, test, "ALUGridType<2, Parallel>"); - reader_writer_test<ALUGridType<3>, ParallelGridCreator<ALUGridType<3>>>(mpi, test, "ALUGridType<3, Parallel>"); -// #endif + if (HasParallelGridFactory<ALUGridType<2>>{}) + reader_writer_test<ALUGridType<2>, ParallelGridCreator<ALUGridType<2>>>(mpi, test, "ALUGridType<2, Parallel>", false); + if (HasParallelGridFactory<ALUGridType<3>>{}) + reader_writer_test<ALUGridType<3>, ParallelGridCreator<ALUGridType<3>>>(mpi, test, "ALUGridType<3, Parallel>", false); #endif return test.exit(); diff --git a/src/test/reader_writer_test.cc b/src/test/reader_writer_test.cc index 860ff561e2c8aa4a42422a9405821a273f6e5486..99b3b88e2eda706ab719c7966a4c672533135002 100644 --- a/src/test/reader_writer_test.cc +++ b/src/test/reader_writer_test.cc @@ -105,27 +105,31 @@ void reader_test (MPIHelper& mpi, Test& test) ext = ".pvtu"; for (auto const& test_case : test_cases) { - GridFactory<Grid> factory; - VtkReader<Grid> reader{factory}; - reader.readFromFile("reader_writer_test_" + std::get<0>(test_case) + ext); - - std::unique_ptr<Grid> grid = Hybrid::ifElse(IsALUGrid<Grid>{}, - [&](auto id) { return id(factory).createGrid(std::true_type{}); }, - [&](auto id) { return id(factory).createGrid(); }); - std::vector<std::string> pieces1 = grid->comm().size() > 1 ? - reader.pieces() : - std::vector<std::string>{"reader_writer_test_" + std::get<0>(test_case) + ".vtu"}; - - VtkUnstructuredGridWriter<typename Grid::LeafGridView> vtkWriter(grid->leafGridView(), - std::get<1>(test_case), std::get<2>(test_case)); - vtkWriter.write("reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"); - - GridFactory<Grid> factory2; - VtkReader<Grid> reader2{factory2}; - reader2.readFromFile("reader_writer_test_" + std::get<0>(test_case) + "_2" + ext, false); - std::vector<std::string> pieces2 = grid->comm().size() > 1 ? - reader.pieces() : - std::vector<std::string>{"reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"}; + std::vector<std::string> pieces1, pieces2; + + { + GridFactory<Grid> factory; + VtkReader<Grid> reader{factory}; + reader.readFromFile("reader_writer_test_" + std::get<0>(test_case) + ext); + + std::unique_ptr<Grid> grid{ Hybrid::ifElse(IsALUGrid<Grid>{}, + [&](auto id) { return id(factory).createGrid(std::true_type{}); }, + [&](auto id) { return id(factory).createGrid(); }) }; + pieces1 = mpi.size() > 1 ? reader.pieces() : + std::vector<std::string>{"reader_writer_test_" + std::get<0>(test_case) + ".vtu"}; + + VtkUnstructuredGridWriter<typename Grid::LeafGridView> vtkWriter(grid->leafGridView(), + std::get<1>(test_case), std::get<2>(test_case)); + vtkWriter.write("reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"); + } + + { + GridFactory<Grid> factory2; + VtkReader<Grid> reader2{factory2}; + reader2.readFromFile("reader_writer_test_" + std::get<0>(test_case) + "_2" + ext, false); + pieces2 = mpi.size() > 1 ? reader2.pieces() : + std::vector<std::string>{"reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"}; + } test.check(pieces1.size() == pieces2.size(), "pieces1.size == pieces2.size"); for (std::size_t i = 0; i < pieces1.size(); ++i)