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/.gitlab-ci.yml b/.gitlab-ci.yml index e473bb7be28cec6476e136823ad690b9dae0e64e..9d718d9758574d77b3771b16609766418c18e431 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -3,18 +3,26 @@ before_script: - duneci-install-module https://gitlab.dune-project.org/martin.nolte/dune-polygongrid.git - duneci-install-module https://gitlab.dune-project.org/extensions/dune-spgrid.git -debian:10 gcc-8-17: - image: registry.dune-project.org/docker/ci/dune:2.6-debian-10-gcc-8-17 +dune-2.6 gcc-6-14: + image: registry.dune-project.org/docker/ci/dune-pdelab-deps:2.6-debian-9-gcc-6-14 script: duneci-standard-test -debian:10 clang-6-libcpp-17: - image: registry.dune-project.org/docker/ci/dune:2.6-debian-10-clang-6-libcpp-17 +dune-2.6 gcc-8-17: + image: registry.dune-project.org/docker/ci/dune-pdelab-deps:2.6-debian-10-gcc-8-17 script: duneci-standard-test -debian:9 gcc-6-14: - image: registry.dune-project.org/docker/ci/dune:2.6-debian-9-gcc-6-14 +dune-2.6 clang-7-17: + image: registry.dune-project.org/docker/ci/dune-pdelab-deps:2.6-debian-10-clang-7-libcpp-17 script: duneci-standard-test -ubuntu:18.04 clang-6-17: - image: registry.dune-project.org/docker/ci/dune:2.6-ubuntu-18.04-clang-6-17 +dune-git gcc-7-14: + image: registry.dune-project.org/docker/ci/dune-pdelab-deps:git-debian-10-gcc-7-14 script: duneci-standard-test + +dune-git gcc-8-17: + image: registry.dune-project.org/docker/ci/dune-pdelab-deps:git-debian-10-gcc-8-noassert-17 + script: duneci-standard-test + +dune-git clang-7-17: + image: registry.dune-project.org/docker/ci/dune-pdelab-deps:git-debian-10-clang-7-libcpp-17 + script: duneci-standard-test \ No newline at end of file diff --git a/dune.module b/dune.module index bbfcce42ffbd535a14b3fd5a3bc3f54d782ae2d1..4c1f8ee87949d93700f079752e9e5b0bafaa798a 100644 --- a/dune.module +++ b/dune.module @@ -4,7 +4,7 @@ #Name of the module Module: dune-vtk -Version: 0.1 +Version: 0.2 Maintainer: simon.praetorius@tu-dresden.de #depending on Depends: dune-common (>= 2.6) dune-geometry (>= 2.6) dune-grid diff --git a/dune/vtk/CMakeLists.txt b/dune/vtk/CMakeLists.txt index 7c6d8a5a9bc9352d3851e9b8d6324e9fb9df4433..64ee710adc312727d0bffe819454aad69540c1a0 100644 --- a/dune/vtk/CMakeLists.txt +++ b/dune/vtk/CMakeLists.txt @@ -8,7 +8,6 @@ install(FILES defaultvtkfunction.hh filereader.hh filewriter.hh - gridcreator.hh legacyvtkfunction.hh pvdwriter.hh pvdwriter.impl.hh @@ -26,5 +25,6 @@ install(FILES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vtkwriter) add_subdirectory(datacollectors) +add_subdirectory(gridcreators) add_subdirectory(utility) add_subdirectory(writers) diff --git a/dune/vtk/datacollectorinterface.hh b/dune/vtk/datacollectorinterface.hh index 8d967f3afcd22d853a96e80c3ebbaa38ba2b7b83..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) {} @@ -19,13 +30,19 @@ public: asDerived().updateImpl(); } - /// Return the number of overlapping elements + /// Return the number of ghost elements int ghostLevel () const { 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 36c4c692751a79b46810cc974479ce022a8728b2..192af3871be2478cec4cff1df8fadd48c90d24c2 100644 --- a/dune/vtk/datacollectors/continuousdatacollector.hh +++ b/dune/vtk/datacollectors/continuousdatacollector.hh @@ -3,43 +3,75 @@ #include <numeric> #include "unstructureddatacollector.hh" +#include <dune/grid/utility/globalindexset.hh> + 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; + } + + /// Return a vector of global unique ids of the points + std::vector<std::uint64_t> pointIdsImpl () const + { + std::vector<std::uint64_t> data; + data.reserve(numPoints_); + GlobalIndexSet<GridView> globalIndexSet(gridView_, dim); + for (auto const& vertex : vertices(gridView_, partition)) { + data.emplace_back(globalIndexSet.index(vertex)); } return data; } @@ -47,7 +79,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 @@ -62,15 +94,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()); } @@ -81,15 +113,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))); } @@ -100,6 +132,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 01c67603f2b8c78def753a7792bad9863c1f4173..05a9f07f9370bede3cb37a2a622ed8161c29052f 100644 --- a/dune/vtk/datacollectors/structureddatacollector.hh +++ b/dune/vtk/datacollectors/structureddatacollector.hh @@ -9,18 +9,6 @@ namespace Dune { - -namespace Impl -{ - // Should be specialized for concrete structured grid - template <class GridView, class Grid> - struct StructuredDataCollectorImpl; -} - -template <class GridView> -using StructuredDataCollector = typename Impl::StructuredDataCollectorImpl<GridView, typename GridView::Grid>::type; - - /// The Interface for structured data-collectors template <class GridView, class Derived> class StructuredDataCollectorInterface @@ -31,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) @@ -115,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 { @@ -204,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 e9a55fb1fcb378ddeeada6f3c9b6c927bbb42808..802cf7c5c107f9bcf35d21c381d5bf119f71b62d 100644 --- a/dune/vtk/datacollectors/unstructureddatacollector.hh +++ b/dune/vtk/datacollectors/unstructureddatacollector.hh @@ -14,29 +14,39 @@ 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 { return this->asDerived().cellsImpl(); } + std::vector<std::uint64_t> pointIds () const + { + return this->asDerived().pointIdsImpl(); + } + +protected: + // default implementation + std::vector<std::uint64_t> pointIdsImpl () const + { + return {}; + } + protected: using Super::gridView_; }; 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/filereader.hh b/dune/vtk/filereader.hh index abc9d0c7a66ce946bbaca9d51f6de83bbcebc38f..ff7e36ccad09cdde75d4f4a9883f781d27eb4ff7 100644 --- a/dune/vtk/filereader.hh +++ b/dune/vtk/filereader.hh @@ -4,6 +4,8 @@ #include <string> #include <utility> +#include <dune/vtk/forward.hh> + namespace Dune { template <class Grid, class FilerReaderImp> diff --git a/dune/vtk/filewriter.hh b/dune/vtk/filewriter.hh index 8873e75644fffecabdb1fc4b20026d8943c7e4ae..376018506398c6b41afcd10380821db52cbf37b7 100644 --- a/dune/vtk/filewriter.hh +++ b/dune/vtk/filewriter.hh @@ -3,6 +3,8 @@ #include <string> #include <dune/common/std/optional.hh> +#include <dune/vtk/forward.hh> + namespace Dune { class FileWriter diff --git a/dune/vtk/forward.hh b/dune/vtk/forward.hh new file mode 100644 index 0000000000000000000000000000000000000000..3ade39f6b098dad6e5c9df4f7363571b6a2dc2d7 --- /dev/null +++ b/dune/vtk/forward.hh @@ -0,0 +1,104 @@ +#pragma once + +namespace Dune +{ + // forward declaration of all classes in dune-vtk + + template <class GridView, class Derived, class Partition = Partitions::InteriorBorder> + class DataCollectorInterface; + + // @{ datacollectors + template <class GridView, class Derived, class Partition = Partitions::InteriorBorder> + class UnstructuredDataCollectorInterface; + + // @{ unstructured-datacollectors + template <class GridView, class Partition = Partitions::InteriorBorder> + class ContinuousDataCollector; + + template <class GridView, class Partition = Partitions::InteriorBorder> + class DiscontinuousDataCollector; + + template <class GridView> + class QuadraticDataCollector; + // @} unstructured-datacollectors + + template <class GridView, class Derived> + class StructuredDataCollectorInterface; + + namespace Impl + { + // Should be specialized for concrete structured grid + template <class GridView, class Grid> + struct StructuredDataCollectorImpl; + } + + template <class GridView> + using StructuredDataCollector = typename Impl::StructuredDataCollectorImpl<GridView, typename GridView::Grid>::type; + + // @{ structured-datacollectors + template <class GridView> + class SPDataCollector; + + template <class GridView> + class YaspDataCollector; + // @} structured-datacollectors + + // @} datacollectors + + template <class Grid, class Derived> + class GridCreatorInterface; + + template <class GridCreator, class Derived> + struct DerivedGridCreator; + + // @{ gridcreators + template <class Grid> + struct ContinuousGridCreator; + + template <class Grid> + struct DiscontinuousGridCreator; + + template <class Grid> + struct ParallelGridCreator; + + template <class Grid> + struct SerialGridCreator; + // @} gridcreators + + + template <class Grid, class FilerReaderImp> + class FileReader; + + template <class Grid, class GridCreator = ContinuousGridCreator<Grid>> + class VtkReader; + + + class FileWriter; + + // @{ filewriters + template <class VtkWriter> + class PvdWriter; + + template <class VtkWriter> + class VtkTimeseriesWriter; + + template <class GridView, class DataCollector> + class VtkWriterInterface; + + // @{ vtkwriters + template <class GridView, class DataCollector = StructuredDataCollector<GridView>> + class VtkImageDataWriter; + + template <class GridView, class DataCollector = StructuredDataCollector<GridView>> + class VtkRectilinearGridWriter; + + template <class GridView, class DataCollector = StructuredDataCollector<GridView>> + class VtkStructuredGridWriter; + + template <class GridView, class DataCollector = ContinuousDataCollector<GridView>> + class VtkUnstructuredGridWriter; + // @} vtkwriters + + // @} filewriters + +} // end namespace Dune diff --git a/dune/vtk/gridcreator.hh b/dune/vtk/gridcreator.hh deleted file mode 100644 index 4fec1291e6a7f9d3f86d86386b30716056485627..0000000000000000000000000000000000000000 --- a/dune/vtk/gridcreator.hh +++ /dev/null @@ -1,119 +0,0 @@ -#pragma once - -#include <cassert> -#include <cstdint> -#include <limits> -#include <vector> - -#include <dune/common/exceptions.hh> -#include <dune/grid/common/gridfactory.hh> - -#include "vtktypes.hh" - -namespace Dune -{ - // Create a grid where the input points and connectivity is already - // connected correctly. - struct DefaultGridCreator - { - template <class Grid, class Coord> - static void create (GridFactory<Grid>& factory, - std::vector<Coord> const& points, - std::vector<std::uint8_t> const& types, - std::vector<std::int64_t> const& offsets, - std::vector<std::int64_t> const& connectivity) - { - for (auto const& p : points) - factory.insertVertex(p); - - std::size_t idx = 0; - for (std::size_t i = 0; i < types.size(); ++i) { - auto type = Vtk::to_geometry(types[i]); - Vtk::CellType cellType{type}; - auto refElem = referenceElement<double,Grid::dimension>(type); - - int nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]); - assert(nNodes == refElem.size(Grid::dimension)); - std::vector<unsigned int> vtk_cell; vtk_cell.reserve(nNodes); - for (int j = 0; j < nNodes; ++j) - vtk_cell.push_back( connectivity[idx++] ); - - if (cellType.noPermutation()) - factory.insertElement(type,vtk_cell); - 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); - } - } - } - }; - - - // Create a grid where the input points are not connected and the connectivity - // describes separated elements. - struct ConnectedGridCreator - { - struct CoordLess - { - template <class T, int N> - bool operator() (FieldVector<T,N> const& lhs, FieldVector<T,N> const& rhs) const - { - for (int i = 0; i < N; ++i) { - if (std::abs(lhs[i] - rhs[i]) < std::numeric_limits<T>::epsilon()) - continue; - return lhs[i] < rhs[i]; - } - return false; - } - }; - - template <class Grid, class Coord> - static void create (GridFactory<Grid>& factory, - std::vector<Coord> const& points, - std::vector<std::uint8_t> const& types, - std::vector<std::int64_t> const& offsets, - std::vector<std::int64_t> const& connectivity) - { - std::size_t idx = 0; - std::map<Coord, std::size_t, CoordLess> unique_points; - for (auto const& p : points) { - auto b = unique_points.emplace(std::make_pair(p,idx)); - if (b.second) { - factory.insertVertex(p); - ++idx; - } - } - - idx = 0; - for (std::size_t i = 0; i < types.size(); ++i) { - auto type = Vtk::to_geometry(types[i]); - Vtk::CellType cellType{type}; - - int nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]); - assert(nNodes > 0); - std::vector<unsigned int> vtk_cell; vtk_cell.reserve(nNodes); - for (int j = 0; j < nNodes; ++j) { - std::size_t v_j = connectivity[idx++]; - std::size_t new_idx = unique_points[points[v_j]]; - vtk_cell.push_back(new_idx); - } - - if (cellType.noPermutation()) - factory.insertElement(type,vtk_cell); - 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); - } - } - } - }; - -} // end namespace Dune diff --git a/dune/vtk/gridcreatorinterface.hh b/dune/vtk/gridcreatorinterface.hh new file mode 100644 index 0000000000000000000000000000000000000000..56fb6f0352c1b1b3b04c936e682617bc46bea85a --- /dev/null +++ b/dune/vtk/gridcreatorinterface.hh @@ -0,0 +1,95 @@ +#pragma once + +#include <cstdint> +#include <string> +#include <vector> + +#include <dune/common/version.hh> +#include <dune/grid/common/gridfactory.hh> + +#include <dune/vtk/forward.hh> + +namespace Dune { + +template <class G, class Derived> +class GridCreatorInterface +{ +public: + using Grid = G; + using GlobalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::GlobalCoordinate; + +public: + GridCreatorInterface (GridFactory<Grid>& factory) + : factory_(&factory) + {} + + /// Insert all points as vertices into the factory + void insertVertices (std::vector<GlobalCoordinate> const& points, + std::vector<std::uint64_t> const& point_ids) + { + asDerived().insertVerticesImpl(points, point_ids); + } + + /// Create elements based on type and connectivity description + void insertElements (std::vector<std::uint8_t> const& types, + std::vector<std::int64_t> const& offsets, + std::vector<std::int64_t> const& connectivity) + { + asDerived().insertElementsImpl(types, offsets, connectivity); + } + + /// Insert part of a grid stored in file into factory + void insertPieces (std::vector<std::string> const& pieces) + { + asDerived().insertPiecesImpl(pieces); + } + + /// Return the associated GridFactory + GridFactory<Grid>& factory () + { + return *factory_; + } + + /// Return the mpi collective communicator + auto comm () const + { + return MPIHelper::getCollectiveCommunication(); + } + +protected: // cast to derived type + + Derived& asDerived () + { + return static_cast<Derived&>(*this); + } + + const Derived& asDerived () const + { + return static_cast<const Derived&>(*this); + } + +public: // default implementations + + void insertVerticesImpl (std::vector<GlobalCoordinate> const&, + std::vector<std::uint64_t> const&) + { + /* do nothing */ + } + + void insertElementsImpl (std::vector<std::uint8_t> const&, + std::vector<std::int64_t> const&, + std::vector<std::int64_t> const&) + { + /* do nothing */ + } + + void insertPiecesImpl (std::vector<std::string> const&) + { + /* do nothing */; + } + +protected: + GridFactory<Grid>* factory_; +}; + +} // end namespace Dune diff --git a/dune/vtk/gridcreators/CMakeLists.txt b/dune/vtk/gridcreators/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..299b4d3351ee9b10b178018221801d1ed696363d --- /dev/null +++ b/dune/vtk/gridcreators/CMakeLists.txt @@ -0,0 +1,9 @@ +#install headers +install(FILES + common.hh + continuousgridcreator.hh + derivedgridcreator.hh + discontinuousgridcreator.hh + parallelgridcreator.hh + serialgridcreator.hh + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vtkwriter/gridcreators) diff --git a/dune/vtk/gridcreators/common.hh b/dune/vtk/gridcreators/common.hh new file mode 100644 index 0000000000000000000000000000000000000000..e47eb4dcd50ee57d37b7bd746972ab52c694b00f --- /dev/null +++ b/dune/vtk/gridcreators/common.hh @@ -0,0 +1,22 @@ +#pragma once + +#include <type_traits> + +namespace Dune +{ + template <class Factory, class... Args> + using HasInsertVertex = decltype( std::declval<Factory>().insertVertex(std::declval<Args>()...) ); + + namespace Impl + { + template <class GF, class = void> + struct VertexIdType { using type = unsigned int; }; + + template <class GF> + struct VertexIdType<GF, typename GF::VertexId> { using type = typename GF::VertexId; }; + } + + template <class GF> + using VertexId_t = typename Impl::VertexIdType<GF>::type; + +} // end namespace Dune diff --git a/dune/vtk/gridcreators/continuousgridcreator.hh b/dune/vtk/gridcreators/continuousgridcreator.hh new file mode 100644 index 0000000000000000000000000000000000000000..03270a268d7fd78f1d91c74fdeb11da896ededa1 --- /dev/null +++ b/dune/vtk/gridcreators/continuousgridcreator.hh @@ -0,0 +1,69 @@ +#pragma once + +#include <cassert> +#include <cstdint> +#include <limits> +#include <vector> + +#include <dune/common/exceptions.hh> +#include <dune/common/hybridutilities.hh> +#include <dune/grid/common/gridfactory.hh> + +#include <dune/vtk/vtktypes.hh> +#include <dune/vtk/gridcreatorinterface.hh> + +namespace Dune +{ + // Create a grid where the input points and connectivity is already + // connected correctly. + template <class Grid> + struct ContinuousGridCreator + : public GridCreatorInterface<Grid, ContinuousGridCreator<Grid>> + { + using Super = GridCreatorInterface<Grid, ContinuousGridCreator<Grid>>; + using GlobalCoordinate = typename Super::GlobalCoordinate; + + ContinuousGridCreator (GridFactory<Grid>& factory) + : Super(factory) + {} + + using Super::factory; + + void insertVerticesImpl (std::vector<GlobalCoordinate> const& points, + std::vector<std::uint64_t> const& /*point_ids*/) + { + for (auto const& p : points) + factory().insertVertex(p); + } + + void insertElementsImpl (std::vector<std::uint8_t> const& types, + std::vector<std::int64_t> const& offsets, + std::vector<std::int64_t> const& connectivity) + { + std::size_t idx = 0; + for (std::size_t i = 0; i < types.size(); ++i) { + auto type = Vtk::to_geometry(types[i]); + Vtk::CellType cellType{type}; + auto refElem = referenceElement<double,Grid::dimension>(type); + + int nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]); + assert(nNodes == refElem.size(Grid::dimension)); + std::vector<unsigned int> vtk_cell; vtk_cell.reserve(nNodes); + for (int j = 0; j < nNodes; ++j) + vtk_cell.push_back( connectivity[idx++] ); + + if (cellType.noPermutation()) + factory().insertElement(type,vtk_cell); + 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); + } + } + } + }; + +} // end namespace Dune diff --git a/dune/vtk/gridcreators/derivedgridcreator.hh b/dune/vtk/gridcreators/derivedgridcreator.hh new file mode 100644 index 0000000000000000000000000000000000000000..91fa8463f028872f9160a8cce9638945615399c3 --- /dev/null +++ b/dune/vtk/gridcreators/derivedgridcreator.hh @@ -0,0 +1,51 @@ +#pragma once + +#include <cstdint> +#include <string> +#include <vector> + +#include <dune/grid/common/gridfactory.hh> +#include <dune/vtk/forward.hh> +#include <dune/vtk/gridcreatorinterface.hh> +#include <dune/vtk/gridcreators/common.hh> +#include <dune/vtk/gridcreators/continuousgridcreator.hh> + +namespace Dune +{ + template <class GridCreator, class Derived> + struct DerivedGridCreator + : public GridCreatorInterface<typename GridCreator::Grid, Derived> + { + using Self = DerivedGridCreator; + using Super = GridCreatorInterface<typename GridCreator::Grid, Derived>; + using Grid = typename GridCreator::Grid; + using GlobalCoordinate = typename Super::GlobalCoordinate; + + DerivedGridCreator (GridFactory<Grid>& factory) + : Super(factory) + , gridCreator_(factory) + {} + + void insertVerticesImpl (std::vector<GlobalCoordinate> const& points, + std::vector<std::uint64_t> const& point_ids) + { + gridCreator_.insertVertices(points, point_ids); + } + + void insertElementsImpl (std::vector<std::uint8_t> const& types, + std::vector<std::int64_t> const& offsets, + std::vector<std::int64_t> const& connectivity) + { + gridCreator_.insertElements(types, offsets, connectivity); + } + + void insertPiecesImpl (std::vector<std::string> const& pieces) + { + gridCreator_.insertPieces(pieces); + } + + private: + GridCreator gridCreator_; + }; + +} // end namespace Dune diff --git a/dune/vtk/gridcreators/discontinuousgridcreator.hh b/dune/vtk/gridcreators/discontinuousgridcreator.hh new file mode 100644 index 0000000000000000000000000000000000000000..5160b00e92a29521b9c42ac3814ce8b2901da63c --- /dev/null +++ b/dune/vtk/gridcreators/discontinuousgridcreator.hh @@ -0,0 +1,96 @@ +#pragma once + +#include <cassert> +#include <cstdint> +#include <limits> +#include <vector> + +#include <dune/common/exceptions.hh> +#include <dune/common/hybridutilities.hh> +#include <dune/grid/common/gridfactory.hh> + +#include <dune/vtk/vtktypes.hh> +#include <dune/vtk/gridcreatorinterface.hh> +namespace Dune +{ + // Create a grid where the input points are not connected and the connectivity + // describes separated elements. + template <class Grid> + struct DiscontinuousGridCreator + : public GridCreatorInterface<Grid, DiscontinuousGridCreator<Grid>> + { + using Super = GridCreatorInterface<Grid, DiscontinuousGridCreator<Grid>>; + using GlobalCoordinate = typename Super::GlobalCoordinate; + + struct CoordLess + { + template <class T, int N> + bool operator() (FieldVector<T,N> const& lhs, FieldVector<T,N> const& rhs) const + { + for (int i = 0; i < N; ++i) { + if (std::abs(lhs[i] - rhs[i]) < std::numeric_limits<T>::epsilon()) + continue; + return lhs[i] < rhs[i]; + } + return false; + } + }; + + DiscontinuousGridCreator (GridFactory<Grid>& factory) + : Super(factory) + {} + + using Super::factory; + void insertVerticesImpl (std::vector<GlobalCoordinate> const& points, + std::vector<std::uint64_t> const& /*point_ids*/) + { + points_ = &points; + uniquePoints_.clear(); + std::size_t idx = 0; + + for (auto const& p : points) { + auto b = uniquePoints_.emplace(std::make_pair(p,idx)); + if (b.second) { + factory().insertVertex(p); + ++idx; + } + } + } + + void insertElementsImpl (std::vector<std::uint8_t> const& types, + std::vector<std::int64_t> const& offsets, + std::vector<std::int64_t> const& connectivity) + { + assert(points_ != nullptr); + std::size_t idx = 0; + for (std::size_t i = 0; i < types.size(); ++i) { + auto type = Vtk::to_geometry(types[i]); + Vtk::CellType cellType{type}; + + int nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]); + assert(nNodes > 0); + std::vector<unsigned int> vtk_cell; vtk_cell.reserve(nNodes); + for (int j = 0; j < nNodes; ++j) { + std::size_t v_j = connectivity[idx++]; + std::size_t new_idx = uniquePoints_[(*points_)[v_j]]; + vtk_cell.push_back(new_idx); + } + + if (cellType.noPermutation()) { + factory().insertElement(type,vtk_cell); + } 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); + } + } + } + + private: + std::vector<GlobalCoordinate> const* points_ = nullptr; + std::map<GlobalCoordinate, std::size_t, CoordLess> uniquePoints_; + }; + +} // end namespace Dune diff --git a/dune/vtk/gridcreators/parallelgridcreator.hh b/dune/vtk/gridcreators/parallelgridcreator.hh new file mode 100644 index 0000000000000000000000000000000000000000..fc2a08c09bb7ca1d469597830bba5646c9a9ee1e --- /dev/null +++ b/dune/vtk/gridcreators/parallelgridcreator.hh @@ -0,0 +1,49 @@ +#pragma once + +#include <cstdint> +#include <string> +#include <vector> + +#include <dune/grid/common/gridfactory.hh> +#include <dune/vtk/forward.hh> +#include <dune/vtk/gridcreatorinterface.hh> +#include <dune/vtk/gridcreators/common.hh> +#include <dune/vtk/gridcreators/derivedgridcreator.hh> +#include <dune/vtk/gridcreators/continuousgridcreator.hh> + +namespace Dune +{ + template <class Grid> + struct ParallelGridCreator + : public DerivedGridCreator<ContinuousGridCreator<Grid>, ParallelGridCreator<Grid>> + { + using Self = ParallelGridCreator; + using Super = DerivedGridCreator<ContinuousGridCreator<Grid>, Self>; + using GlobalCoordinate = typename Super::GlobalCoordinate; + using VertexId = VertexId_t<GridFactory<Grid>>; + + // The GridFactory must support insertion of global vertex IDs + static_assert(Std::is_detected<HasInsertVertex, GridFactory<Grid>, GlobalCoordinate, VertexId>{}, ""); + + ParallelGridCreator (GridFactory<Grid>& factory) + : 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) + this->factory().insertVertex(points[i], VertexId(point_ids[i])); + } + + void insertPiecesImpl (std::vector<std::string> const& pieces) + { + if (int(pieces.size()) == this->comm().size()) { + VtkReader<Grid, Self> pieceReader(this->factory()); + pieceReader.readFromFile(pieces[this->comm().rank()], true); + } + } + }; + +} // end namespace Dune diff --git a/dune/vtk/gridcreators/serialgridcreator.hh b/dune/vtk/gridcreators/serialgridcreator.hh new file mode 100644 index 0000000000000000000000000000000000000000..75aeb13b76b5ccac3bc5ac94f36d39009e281c40 --- /dev/null +++ b/dune/vtk/gridcreators/serialgridcreator.hh @@ -0,0 +1,73 @@ +#pragma once + +#include <cstdint> +#include <string> +#include <vector> + +#include <dune/grid/common/gridfactory.hh> +#include <dune/vtk/forward.hh> +#include <dune/vtk/gridcreatorinterface.hh> +#include <dune/vtk/gridcreators/discontinuousgridcreator.hh> + +namespace Dune +{ + template <class Grid> + struct SerialGridCreator + : public GridCreatorInterface<Grid, SerialGridCreator<Grid>> + { + using Self = SerialGridCreator; + using Super = GridCreatorInterface<Grid, Self>; + using GlobalCoordinate = typename Super::GlobalCoordinate; + + SerialGridCreator (GridFactory<Grid>& factory) + : Super(factory) + {} + + void insertVerticesImpl (std::vector<GlobalCoordinate> const& points, + std::vector<std::uint64_t> const& /*point_ids*/) + { + shift_.push_back(points_.size()); + points_.reserve(points_.size() + points.size()); + points_.insert(points_.end(), points.begin(), points.end()); + } + + void insertElementsImpl (std::vector<std::uint8_t> const& types, + std::vector<std::int64_t> const& offsets, + std::vector<std::int64_t> const& connectivity) + { + types_.reserve(types_.size() + types.size()); + types_.insert(types_.end(), types.begin(), types.end()); + + offsets_.reserve(offsets_.size() + offsets.size()); + std::transform(offsets.begin(), offsets.end(), std::back_inserter(offsets_), + [shift=offsets_.empty() ? 0 : offsets_.back()](std::int64_t o) { return o + shift; }); + + connectivity_.reserve(connectivity_.size() + connectivity.size()); + std::transform(connectivity.begin(), connectivity.end(), std::back_inserter(connectivity_), + [shift=shift_.back()](std::int64_t idx) { return idx + shift; }); + } + + void insertPiecesImpl (std::vector<std::string> const& pieces) + { + if (this->comm().rank() == 0) { + VtkReader<Grid, Self> pieceReader(*this); + for (std::string const& piece : pieces) { + pieceReader.readFromFile(piece, false); + pieceReader.createGrid(false); + } + + DiscontinuousGridCreator<Grid> creator(this->factory()); + creator.insertVertices(points_, {}); + creator.insertElements(types_, offsets_, connectivity_); + } + } + + private: + std::vector<GlobalCoordinate> points_; + std::vector<std::uint8_t> types_; + std::vector<std::int64_t> offsets_; + std::vector<std::int64_t> connectivity_; + std::vector<std::int64_t> shift_; + }; + +} // end namespace Dune diff --git a/dune/vtk/pvdwriter.hh b/dune/vtk/pvdwriter.hh index 46343bd01d10381c5b30da2379958de38f8f6ac8..8fe864e6386a3ed4717cfcd258a3fa11d9ceb661 100644 --- a/dune/vtk/pvdwriter.hh +++ b/dune/vtk/pvdwriter.hh @@ -7,6 +7,7 @@ #include <dune/vtk/vtktypes.hh> #include <dune/vtk/filewriter.hh> +#include <dune/vtk/forward.hh> namespace Dune { 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 3b90c76706f5e0036f07a374cfd8ea78b84b4160..d49199dd3061f645ef83e7270294e57f155834ce 100644 --- a/dune/vtk/vtkreader.hh +++ b/dune/vtk/vtkreader.hh @@ -3,20 +3,23 @@ #include <iosfwd> #include <map> -#include "filereader.hh" -#include "gridcreator.hh" -#include "vtktypes.hh" +#include <dune/vtk/filereader.hh> +#include <dune/vtk/forward.hh> +#include <dune/vtk/vtktypes.hh> + +// default GridCreator +#include <dune/vtk/gridcreators/continuousgridcreator.hh> namespace Dune { - /// File-Reader for Vtk .vtu files + /// File-Reader for Vtk unstructured .vtu files /** * Reads .vtu files and constructs a grid from the cells stored in the file * Additionally, stored data can be read. * * Assumption on the file structure: Each XML tag must be on a separate line. **/ - template <class Grid, class GridCreator = DefaultGridCreator> + template <class Grid, class GridCreator> class VtkReader : public FileReader<Grid, VtkReader<Grid, GridCreator>> { @@ -37,13 +40,26 @@ namespace Dune using GlobalCoordinate = typename Entity::Geometry::GlobalCoordinate; public: - /// Constructor. Stores a pointer to the GridFactory. - VtkReader (GridFactory<Grid>& factory) - : factory_(&factory) + /// Constructor. Creates a new GridCreator with the passed factory + template <class... Args> + explicit VtkReader (Args&&... args) + : creatorStorage_(std::make_unique<GridCreator>(std::forward<Args>(args)...)) + , creator_(*creatorStorage_) + {} + + /// Constructor. Stores a references to the passed creator + VtkReader (GridCreator& creator) + : creator_(creator) {} - /// Read the grid from file with `filename` into the GridFactory `factory` - void readFromFile (std::string const& filename); + /// Read the grid from file with `filename` into the GridFactory \ref factory_ + void readFromFile (std::string const& filename, bool create = true); + + /// Read the grid from and input stream into the GridFactory \ref factory_ + void readSerialFileFromStream (std::ifstream& input, bool create = true); + + /// Read the grid from and input stream into the GridFactory \ref factory_ + void readParallelFileFromStream (std::ifstream& input, int rank, int size, bool create = true); /// Implementation of \ref FileReader interface static void readFactoryImpl (GridFactory<Grid>& factory, std::string const& filename) @@ -52,6 +68,15 @@ namespace Dune reader.readFromFile(filename); } + /// Construct a grid using the GridCreator + void createGrid(bool insertPieces = true); + + /// Return the filenames of parallel pieces + std::vector<std::string> const& pieces () const + { + return pieces_; + } + private: // Read values stored on the cells with name `name` template <class T> @@ -69,7 +94,7 @@ namespace Dune } // Read vertex coordinates from `input` stream and store in into `factory` - Sections readPoints (std::ifstream& input); + Sections readPoints (std::ifstream& input, std::string name); template <class T> void readPointsAppended (std::ifstream& input); @@ -101,32 +126,40 @@ namespace Dune // Read attributes from current xml tag std::map<std::string, std::string> parseXml(std::string const& line, bool& closed); - // Construct a grid using the GridFactory `factory` and the read vectors - // \ref vec_points, \ref vec_types, \ref vec_offsets, and \ref vec_connectivity, - // by passing to \ref GridCreator. - void createGrid() const; + // clear all vectors + void clear (); + + auto comm () const + { + return MPIHelper::getCollectiveCommunication(); + } private: - GridFactory<Grid>* factory_; + std::unique_ptr<GridCreator> creatorStorage_ = nullptr; + GridCreator& creator_; /// Data format, i.e. ASCII, BINARY or COMPRESSED. Read from xml attributes. Vtk::FormatTypes format_; // Temporary data to construct the grid elements std::vector<GlobalCoordinate> vec_points; - std::vector<std::uint8_t> vec_types; //< VTK cell type ID - std::vector<std::int64_t> vec_offsets; //< offset of vertices of cell + std::vector<std::uint64_t> vec_point_ids; //< Global unique vertex ID + std::vector<std::uint8_t> vec_types; //< VTK cell type ID + std::vector<std::int64_t> vec_offsets; //< offset of vertices of cell std::vector<std::int64_t> vec_connectivity; //< vertex indices of cell - std::size_t numberOfCells_; //< Number of cells in the grid - std::size_t numberOfPoints_; // Number of vertices in the grid + std::size_t numberOfCells_ = 0; //< Number of cells in the grid + std::size_t numberOfPoints_ = 0; //< Number of vertices in the grid // offset information for appended data // map Name -> {DataType,NumberOfComponents,Offset} std::map<std::string, DataArrayAttributes> dataArray_; + // vector of filenames of parallel pieces + std::vector<std::string> pieces_; + /// Offset of beginning of appended data - std::uint64_t offset0_; + std::uint64_t offset0_ = 0; }; } // end namespace Dune diff --git a/dune/vtk/vtkreader.impl.hh b/dune/vtk/vtkreader.impl.hh index 0d19e794e142e3053353ddddaee435ec5e9a64be..3fd29ba1c1503ab5d882f9b302416ee44166a27d 100644 --- a/dune/vtk/vtkreader.impl.hh +++ b/dune/vtk/vtkreader.impl.hh @@ -3,11 +3,12 @@ #include <iterator> #include <string> -#ifdef HAVE_VTK_ZLIB +#if HAVE_VTK_ZLIB #include <zlib.h> #endif #include <dune/common/classname.hh> +#include <dune/common/version.hh> #include "utility/filesystem.hh" #include "utility/string.hh" @@ -15,7 +16,7 @@ namespace Dune { template <class Grid, class Creator> -void VtkReader<Grid,Creator>::readFromFile (std::string const& filename) +void VtkReader<Grid,Creator>::readFromFile (std::string const& filename, bool create) { // check whether file exists! if (!filesystem::exists(filename)) @@ -24,6 +25,22 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename) std::ifstream input(filename, std::ios_base::in | std::ios_base::binary); assert(input.is_open()); + std::string ext = filesystem::path(filename).extension().string(); + if (ext == ".vtu") { + readSerialFileFromStream(input, create); + pieces_.push_back(filename); + } else if (ext == ".pvtu") { + readParallelFileFromStream(input, comm().rank(), comm().size(), create); + } else { + DUNE_THROW(IOError, "File has unknown file-extension '" << ext << "'. Allowed are only '.vtu' and '.pvtu'."); + } +} + + +template <class Grid, class Creator> +void VtkReader<Grid,Creator>::readSerialFileFromStream (std::ifstream& input, bool create) +{ + clear(); std::string compressor = ""; std::string data_name = "", data_format = ""; Vtk::DataTypes data_type = Vtk::UNKNOWN; @@ -50,9 +67,6 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename) compressor = attr["compressor"]; assert(compressor == "vtkZLibDataCompressor"); // only ZLib compression supported } - - // std::cout << "<VTKFile type='" << attr["type"] << "' version='" << attr["version"] << "' header_type='" << attr["header_type"] << "' byte_order='" << attr["byte_order"] << "' compressor='" << attr["compressor"] << "'>\n"; - section = VTK_FILE; } else if (isSection(line, "/VTKFile", section, VTK_FILE)) @@ -121,8 +135,6 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename) // Store attributes of DataArray dataArray_[data_name] = {data_type, data_components, data_offset}; - // std::cout << "<DataArray type='" << attr["type"] << "' Name='" << data_name << "' NumberOfComponents='" << attr["NumberOfComponents"] << "' format='" << data_format << "' offset='" << attr["offset"] << "' " << (closed ? "/" : "") << ">\n"; - // Skip section in appended mode if (data_format == "appended") { if (!closed) { @@ -187,7 +199,7 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename) } break; case POINTS_DATA_ARRAY: - section = readPoints(input); + section = readPoints(input, data_name); break; case CD_DATA_ARRAY: if (data_type == Vtk::FLOAT32) { @@ -213,7 +225,59 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename) if (section != NO_SECTION) DUNE_THROW(IOError, "VTK-File is incomplete. It must end with </VTKFile>!"); - createGrid(); + if (create) + createGrid(); +} + + +template <class Grid, class Creator> +void VtkReader<Grid,Creator>::readParallelFileFromStream (std::ifstream& input, int commRank, int commSize, bool create) +{ + clear(); + + Sections section = NO_SECTION; + for (std::string line; std::getline(input, line); ) { + ltrim(line); + + if (isSection(line, "VTKFile", section)) { + bool closed = false; + auto attr = parseXml(line, closed); + + if (!attr["type"].empty()) + assert(attr["type"] == "PUnstructuredGrid"); + if (!attr["version"].empty()) + assert(std::stod(attr["version"]) == 1.0); + if (!attr["byte_order"].empty()) + assert(attr["byte_order"] == "LittleEndian"); + if (!attr["header_type"].empty()) + assert(attr["header_type"] == "UInt64"); + if (!attr["compressor"].empty()) + assert(attr["compressor"] == "vtkZLibDataCompressor"); // only ZLib compression supported + section = VTK_FILE; + } + else if (isSection(line, "/VTKFile", section, VTK_FILE)) + section = NO_SECTION; + else if (isSection(line, "PUnstructuredGrid", section, VTK_FILE)) + section = UNSTRUCTURED_GRID; + else if (isSection(line, "/PUnstructuredGrid", section, UNSTRUCTURED_GRID)) + section = VTK_FILE; + else if (isSection(line, "Piece", section, UNSTRUCTURED_GRID)) { + bool closed = false; + auto attr = parseXml(line, closed); + + assert(attr.count("Source") > 0); + pieces_.push_back(attr["Source"]); + } + + if (section == NO_SECTION) + break; + } + + if (section != NO_SECTION) + DUNE_THROW(IOError, "VTK-File is incomplete. It must end with </VTKFile>!"); + + if (create) + createGrid(); } @@ -250,14 +314,16 @@ Sections readDataArray (IStream& input, std::vector<T>& values, std::size_t max_ template <class Grid, class Creator> typename VtkReader<Grid,Creator>::Sections -VtkReader<Grid,Creator>::readPoints (std::ifstream& input) +VtkReader<Grid,Creator>::readPoints (std::ifstream& input, std::string name) { using T = typename GlobalCoordinate::value_type; assert(numberOfPoints_ > 0); assert(dataArray_["points"].components == 3u); + Sections sec; + std::vector<T> point_values; - auto sec = readDataArray(input, point_values, 3*numberOfPoints_, POINTS_DATA_ARRAY, POINTS); + sec = readDataArray(input, point_values, 3*numberOfPoints_, POINTS_DATA_ARRAY, POINTS); assert(sec == POINTS); assert(point_values.size() == 3*numberOfPoints_); @@ -282,7 +348,6 @@ void VtkReader<Grid,Creator>::readPointsAppended (std::ifstream& input) { assert(numberOfPoints_ > 0); assert(dataArray_["points"].components == 3u); - std::vector<T> point_values; readAppended(input, point_values, dataArray_["points"].offset); assert(point_values.size() == 3*numberOfPoints_); @@ -321,6 +386,9 @@ VtkReader<Grid,Creator>::readCells (std::ifstream& input, std::string name) else max_size = numberOfCells_ * max_vertices; sec = readDataArray(input, vec_connectivity, max_size, CELLS_DATA_ARRAY, CELLS); + } else if (name == "global_point_ids") { + sec = readDataArray(input, vec_point_ids, numberOfPoints_, CELLS_DATA_ARRAY, CELLS); + assert(vec_point_ids.size() == numberOfPoints_); } return sec; @@ -346,6 +414,13 @@ void VtkReader<Grid,Creator>::readCellsAppended (std::ifstream& input) assert(connectivity_data.type == Vtk::INT64); readAppended(input, vec_connectivity, connectivity_data.offset); assert(vec_connectivity.size() == std::size_t(vec_offsets.back())); + + if (dataArray_.count("global_point_ids") > 0) { + auto point_id_data = dataArray_["global_point_ids"]; + assert(point_id_data.type == Vtk::UINT64); + readAppended(input, vec_point_ids, point_id_data.offset); + assert(vec_point_ids.size() == numberOfPoints_); + } } @@ -361,7 +436,7 @@ template <class T, class IStream> void read_compressed (T* buffer, unsigned char* buffer_in, std::uint64_t bs, std::uint64_t cbs, IStream& input) { -#ifdef HAVE_VTK_ZLIB +#if HAVE_VTK_ZLIB uLongf uncompressed_space = uLongf(bs); uLongf compressed_space = uLongf(cbs); @@ -436,13 +511,18 @@ void VtkReader<Grid,Creator>::readAppended (std::ifstream& input, std::vector<T> template <class Grid, class Creator> -void VtkReader<Grid,Creator>::createGrid () const +void VtkReader<Grid,Creator>::createGrid (bool insertPieces) { assert(vec_points.size() == numberOfPoints_); assert(vec_types.size() == numberOfCells_); assert(vec_offsets.size() == numberOfCells_); - Creator::create(*factory_, vec_points, vec_types, vec_offsets, vec_connectivity); + if (!vec_points.empty()) + creator_.insertVertices(vec_points, vec_point_ids); + if (!vec_types.empty()) + creator_.insertElements(vec_types, vec_offsets, vec_connectivity); + if (insertPieces) + creator_.insertPieces(pieces_); } // Assume input already read the line <AppendedData ...> @@ -513,4 +593,21 @@ std::map<std::string, std::string> VtkReader<Grid,Creator>::parseXml (std::strin return attr; } + +template <class Grid, class Creator> +void VtkReader<Grid,Creator>::clear () +{ + vec_points.clear(); + vec_point_ids.clear(); + vec_types.clear(); + vec_offsets.clear(); + vec_connectivity.clear(); + dataArray_.clear(); + pieces_.clear(); + + numberOfCells_ = 0; + numberOfPoints_ = 0; + offset0_ = 0; +} + } // end namespace Dune diff --git a/dune/vtk/vtktimeserieswriter.hh b/dune/vtk/vtktimeserieswriter.hh index 587467731857c948eeb8b2b9c69ee83c15b38402..a2b022e4038e7a799329327b3b6d0e43bac6892a 100644 --- a/dune/vtk/vtktimeserieswriter.hh +++ b/dune/vtk/vtktimeserieswriter.hh @@ -5,6 +5,7 @@ #include <vector> #include <dune/vtk/filewriter.hh> +#include <dune/vtk/forward.hh> #include <dune/vtk/vtktypes.hh> #include <dune/vtk/utility/filesystem.hh> #include <dune/vtk/utility/uid.hh> 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/vtkwriter.hh b/dune/vtk/vtkwriter.hh index 91befcfed7d93ce92232634d0813638239c2eb67..ff4414594356d293a7956060cb4832a1781d8219 100644 --- a/dune/vtk/vtkwriter.hh +++ b/dune/vtk/vtkwriter.hh @@ -12,6 +12,7 @@ #include <dune/grid/geometrygrid.hh> #include <dune/grid/yaspgrid.hh> +#include <dune/vtk/forward.hh> #include <dune/vtk/datacollectors/yaspdatacollector.hh> namespace Dune diff --git a/dune/vtk/vtkwriterinterface.hh b/dune/vtk/vtkwriterinterface.hh index 674bde17324d9608df505f3f68863fecb6f62d0a..b753f415848d6e10b31f168f6da0c5a6124c6033 100644 --- a/dune/vtk/vtkwriterinterface.hh +++ b/dune/vtk/vtkwriterinterface.hh @@ -7,18 +7,12 @@ #include <dune/common/std/optional.hh> #include <dune/vtk/filewriter.hh> +#include <dune/vtk/forward.hh> #include <dune/vtk/vtkfunction.hh> #include <dune/vtk/vtktypes.hh> namespace Dune { - // forward declaration - template <class VtkWriter> - class VtkTimeseriesWriter; - - template <class VtkWriter> - class PvdWriter; - /// Interface for file writers for the Vtk XML file formats template <class GridView, class DataCollector> class VtkWriterInterface @@ -31,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 { @@ -44,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) { @@ -154,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/dune/vtk/writers/vtkimagedatawriter.hh b/dune/vtk/writers/vtkimagedatawriter.hh index 50db5b7e7d38877f1afbda07637cb36bb76a2b02..e5260e3c512b037858529a4ca5353982f4628416 100644 --- a/dune/vtk/writers/vtkimagedatawriter.hh +++ b/dune/vtk/writers/vtkimagedatawriter.hh @@ -5,6 +5,7 @@ #include <map> #include <dune/vtk/filewriter.hh> +#include <dune/vtk/forward.hh> #include <dune/vtk/vtkfunction.hh> #include <dune/vtk/vtktypes.hh> #include <dune/vtk/datacollectors/structureddatacollector.hh> @@ -18,7 +19,7 @@ namespace Dune * Requirement: * - DataCollector must be a model of \ref StructuredDataCollector **/ - template <class GridView, class DataCollector = StructuredDataCollector<GridView>> + template <class GridView, class DataCollector> class VtkImageDataWriter : public VtkWriterInterface<GridView, DataCollector> { diff --git a/dune/vtk/writers/vtkrectilineargridwriter.hh b/dune/vtk/writers/vtkrectilineargridwriter.hh index f7c0317c197abb745f56379bc8c1f7367c5b7252..5bb2da2864496d65f8c4c02facbac20f187c9125 100644 --- a/dune/vtk/writers/vtkrectilineargridwriter.hh +++ b/dune/vtk/writers/vtkrectilineargridwriter.hh @@ -5,6 +5,7 @@ #include <map> #include <dune/vtk/filewriter.hh> +#include <dune/vtk/forward.hh> #include <dune/vtk/vtkfunction.hh> #include <dune/vtk/vtktypes.hh> #include <dune/vtk/datacollectors/structureddatacollector.hh> @@ -18,7 +19,7 @@ namespace Dune * Requirement: * - DataCollector must be a model of \ref StructuredDataCollector **/ - template <class GridView, class DataCollector = StructuredDataCollector<GridView>> + template <class GridView, class DataCollector> class VtkRectilinearGridWriter : public VtkWriterInterface<GridView, DataCollector> { diff --git a/dune/vtk/writers/vtkstructuredgridwriter.hh b/dune/vtk/writers/vtkstructuredgridwriter.hh index 522d4723c4c9409d36f2f26fabb32771f7a25674..d5c240f1926c413cbc7b7df4c979dc7a4b389db2 100644 --- a/dune/vtk/writers/vtkstructuredgridwriter.hh +++ b/dune/vtk/writers/vtkstructuredgridwriter.hh @@ -5,6 +5,7 @@ #include <map> #include <dune/vtk/filewriter.hh> +#include <dune/vtk/forward.hh> #include <dune/vtk/vtkfunction.hh> #include <dune/vtk/vtktypes.hh> #include <dune/vtk/datacollectors/structureddatacollector.hh> @@ -18,7 +19,7 @@ namespace Dune * Requirement: * - DataCollector must be a model of \ref StructuredDataCollector **/ - template <class GridView, class DataCollector = StructuredDataCollector<GridView>> + template <class GridView, class DataCollector> class VtkStructuredGridWriter : public VtkWriterInterface<GridView, DataCollector> { diff --git a/dune/vtk/writers/vtkunstructuredgridwriter.hh b/dune/vtk/writers/vtkunstructuredgridwriter.hh index fb2003fc1d02ba74f9af238b0e4a6f200848233d..ee88bf7552ef1f44f005435a6c7523c2b84cc6b9 100644 --- a/dune/vtk/writers/vtkunstructuredgridwriter.hh +++ b/dune/vtk/writers/vtkunstructuredgridwriter.hh @@ -5,6 +5,7 @@ #include <map> #include <dune/vtk/filewriter.hh> +#include <dune/vtk/forward.hh> #include <dune/vtk/vtkfunction.hh> #include <dune/vtk/vtktypes.hh> #include <dune/vtk/datacollectors/continuousdatacollector.hh> @@ -18,7 +19,7 @@ namespace Dune * Requirement: * - DataCollector must be a model of \ref DataCollector **/ - template <class GridView, class DataCollector = ContinuousDataCollector<GridView>> + template <class GridView, class DataCollector> class VtkUnstructuredGridWriter : public VtkWriterInterface<GridView, DataCollector> { @@ -75,10 +76,14 @@ namespace Dune // Write the element connectivity to the output stream `out`. In case // of binary format, stores the streampos of XML attributes "offset" in the // vector `offsets`. - void writeCells (std::ofstream& oust, + void writeCells (std::ofstream& out, std::vector<pos_type>& offsets, Std::optional<std::size_t> timestep = {}) const; + void writePointIds (std::ofstream& out, + std::vector<pos_type>& offsets, + Std::optional<std::size_t> timestep = {}) const; + private: using Super::dataCollector_; using Super::format_; diff --git a/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh b/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh index d91575187906daf8d243109363d725d608d361c4..4c4b36e5155be11cd4af1c3cc50f259fd0f67f4d 100644 --- a/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh +++ b/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh @@ -37,6 +37,7 @@ void VtkUnstructuredGridWriter<GV,DC> // Write element connectivity, types and offsets out << "<Cells>\n"; writeCells(out, offsets); + writePointIds(out, offsets); out << "</Cells>\n"; // Write data associated with grid points @@ -134,14 +135,17 @@ void VtkUnstructuredGridWriter<GV,DC> // Write point coordinates out << "<Points>\n"; - for (std::size_t i = 0; i < timesteps.size(); ++i) + for (std::size_t i = 0; i < timesteps.size(); ++i) { this->writePoints(out, offsets[i], i); + } out << "</Points>\n"; // Write element connectivity, types and offsets out << "<Cells>\n"; - for (std::size_t i = 0; i < timesteps.size(); ++i) + for (std::size_t i = 0; i < timesteps.size(); ++i) { writeCells(out, offsets[i], i); + writePointIds(out, offsets[i], i); + } out << "</Cells>\n"; const std::size_t shift = offsets[0].size(); // number of blocks to write the grid @@ -329,6 +333,34 @@ void VtkUnstructuredGridWriter<GV,DC> } } +template <class GV, class DC> +void VtkUnstructuredGridWriter<GV,DC> + ::writePointIds (std::ofstream& out, + std::vector<pos_type>& offsets, + Std::optional<std::size_t> timestep) const +{ + auto ids = dataCollector_.pointIds(); + if (ids.empty()) + return; + + if (format_ == Vtk::ASCII) { + out << "<DataArray type=\"UInt64\" Name=\"global_point_ids\" format=\"ascii\""; + if (timestep) + out << " TimeStep=\"" << *timestep << "\""; + out << ">\n"; + this->writeValuesAscii(out, ids); + out << "</DataArray>\n"; + } + else { // Vtk::APPENDED format + out << "<DataArray type=\"UInt64\" Name=\"global_point_ids\" format=\"appended\""; + if (timestep) + out << " TimeStep=\"" << *timestep << "\""; + out << " offset="; + offsets.push_back(out.tellp()); + out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' '); + out << "/>\n"; + } +} template <class GV, class DC> void VtkUnstructuredGridWriter<GV,DC> @@ -346,6 +378,10 @@ void VtkUnstructuredGridWriter<GV,DC> blocks.push_back(this->writeValuesAppended(out, cells.connectivity)); blocks.push_back(this->writeValuesAppended(out, cells.offsets)); blocks.push_back(this->writeValuesAppended(out, cells.types)); + + auto ids = dataCollector_.pointIds(); + if (!ids.empty()) + blocks.push_back(this->writeValuesAppended(out, ids)); } } // end namespace Dune diff --git a/src/test/CMakeLists.txt b/src/test/CMakeLists.txt index 77409fec3c8b51e4b73abb957d624162ff4ce476..86a7f5ffc9fc24f93e4f94656353a139d211aa60 100644 --- a/src/test/CMakeLists.txt +++ b/src/test/CMakeLists.txt @@ -1,6 +1,9 @@ dune_add_test(SOURCES reader_writer_test.cc LINK_LIBRARIES dunevtk) +dune_add_test(SOURCES parallel_reader_writer_test.cc + LINK_LIBRARIES dunevtk) + dune_add_test(SOURCES mixed_element_test.cc LINK_LIBRARIES dunevtk CMAKE_GUARD HAVE_UG) \ No newline at end of file diff --git a/src/test/parallel_reader_writer_test.cc b/src/test/parallel_reader_writer_test.cc new file mode 100644 index 0000000000000000000000000000000000000000..a15261b5cb88ffb614c6f990cd2a0ea989493f29 --- /dev/null +++ b/src/test/parallel_reader_writer_test.cc @@ -0,0 +1,199 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: + +#ifdef HAVE_CONFIG_H +# include "config.h" +#endif + +#include <cstring> +#include <iostream> +#include <vector> + +#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI +#include <dune/common/filledarray.hh> +#include <dune/common/std/type_traits.hh> +#include <dune/common/test/testsuite.hh> + +#include <dune/grid/uggrid.hh> +#include <dune/grid/yaspgrid.hh> +#include <dune/grid/utility/structuredgridfactory.hh> + +#if HAVE_DUNE_ALUGRID +#include <dune/alugrid/grid.hh> +#endif + +#include <dune/vtk/vtkreader.hh> +#include <dune/vtk/writers/vtkunstructuredgridwriter.hh> +#include <dune/vtk/gridcreators/parallelgridcreator.hh> +#include <dune/vtk/gridcreators/serialgridcreator.hh> + +using namespace Dune; + +// see https://stackoverflow.com/questions/6163611/compare-two-files +bool compare_files (std::string const& fn1, std::string const& fn2) +{ + std::ifstream in1(fn1, std::ios::binary); + std::ifstream in2(fn2, std::ios::binary); + if (!in1 || !in2) { + std::cout << "can not find file " << fn1 << " or file " << fn2 << "\n"; + return false; + } + + std::ifstream::pos_type size1 = in1.seekg(0, std::ifstream::end).tellg(); + in1.seekg(0, std::ifstream::beg); + + std::ifstream::pos_type size2 = in2.seekg(0, std::ifstream::end).tellg(); + in2.seekg(0, std::ifstream::beg); + + if (size1 != size2) + return false; + + static const std::size_t BLOCKSIZE = 4096; + std::size_t remaining = size1; + + while (remaining) { + char buffer1[BLOCKSIZE], buffer2[BLOCKSIZE]; + std::size_t size = std::min(BLOCKSIZE, remaining); + + in1.read(buffer1, size); + in2.read(buffer2, size); + + if (0 != std::memcmp(buffer1, buffer2, size)) + return false; + + remaining -= size; + } + + return true; +} + + +template <class GF> +using HasParallelGridFactoryImpl = decltype(std::declval<GF>().createGrid(true,true,std::string(""),true)); + +template <class G> +using HasParallelGridFactory = Std::is_detected<HasParallelGridFactoryImpl, GridFactory<G>>; + + +template <class Test> +void compare (Test& test, filesystem::path const& dir, filesystem::path const& name) +{ + test.check(compare_files(dir.string() + '/' + name.string() + ".vtu", + dir.string() + '/' + name.string() + "_2.vtu")); +} + +template <class GridView> +void writer_test (GridView const& gridView, std::string base_name) +{ + VtkUnstructuredGridWriter<GridView> vtkWriter(gridView, Vtk::ASCII, Vtk::FLOAT32); + vtkWriter.write(base_name + ".vtu"); +} + +template <class Grid, class Creator> +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"; + + // 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(); + writer_test(gridPtr->leafGridView(), base_name); + } + + mpi.getCollectiveCommunication().barrier(); // need a barrier between write and read + + // Step 2: read the grid from file1 and write it back to file2 + { GridFactory<Grid> factory; + 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 + template <int dim> + using ALUGridType = Dune::ALUGrid<dim, dim, Dune::simplex, Dune::conforming>; +#endif + +template <int I> +using int_ = std::integral_constant<int,I>; + +int main (int argc, char** argv) +{ + auto& mpi = Dune::MPIHelper::instance(argc, argv); + + TestSuite test{}; + +#if HAVE_UG + reader_writer_test<UGGrid<2>, SerialGridCreator<UGGrid<2>>>(mpi, test, "UGGrid<2>"); + reader_writer_test<UGGrid<3>, SerialGridCreator<UGGrid<3>>>(mpi, test, "UGGrid<3>"); +#endif + +#if HAVE_DUNE_ALUGRID + // Test VtkWriter for ALUGrid. + reader_writer_test<ALUGridType<2>, SerialGridCreator<ALUGridType<2>>>(mpi, test, "ALUGridType<2>"); + reader_writer_test<ALUGridType<2>, ParallelGridCreator<ALUGridType<2>>>(mpi, test, "ALUGridType<2, Parallel>", false); + + reader_writer_test<ALUGridType<3>, SerialGridCreator<ALUGridType<3>>>(mpi, test, "ALUGridType<3>"); + #if DUNE_VERSION_LT(DUNE_GRID,2,7) + // Currently the 2.7 branch is not working, due to a new bisection compatibility check in 3d + reader_writer_test<ALUGridType<3>, ParallelGridCreator<ALUGridType<3>>>(mpi, test, "ALUGridType<3, Parallel>", false); + #endif +#endif + + return test.exit(); +} diff --git a/src/test/reader_writer_test.cc b/src/test/reader_writer_test.cc index 2421af8caf7106dc3890286c01e1b3ec907d22a8..f0051fc40e5c43cb129bb4b4e88f170e69e9231c 100644 --- a/src/test/reader_writer_test.cc +++ b/src/test/reader_writer_test.cc @@ -17,8 +17,13 @@ #include <dune/grid/yaspgrid.hh> #include <dune/grid/utility/structuredgridfactory.hh> +#if HAVE_DUNE_ALUGRID +#include <dune/alugrid/grid.hh> +#endif + #include <dune/vtk/vtkreader.hh> #include <dune/vtk/writers/vtkunstructuredgridwriter.hh> +#include <dune/vtk/gridcreators/continuousgridcreator.hh> using namespace Dune; @@ -27,6 +32,10 @@ bool compare_files (std::string const& fn1, std::string const& fn2) { std::ifstream in1(fn1, std::ios::binary); std::ifstream in2(fn2, std::ios::binary); + if (!in1 || !in2) { + std::cout << "can not find file " << fn1 << " or file " << fn2 << "\n"; + return false; + } std::ifstream::pos_type size1 = in1.seekg(0, std::ifstream::end).tellg(); in1.seekg(0, std::ifstream::beg); @@ -66,6 +75,13 @@ static TestCases test_cases = { {"zlib64", Vtk::COMPRESSED, Vtk::FLOAT64}, }; +template <class Test> +void compare (Test& test, filesystem::path const& dir, filesystem::path const& name) +{ + test.check(compare_files(dir.string() + '/' + name.string() + ".vtu", + dir.string() + '/' + name.string() + "_2.vtu")); +} + template <class GridView> void writer_test (GridView const& gridView) { @@ -76,15 +92,35 @@ void writer_test (GridView const& gridView) } template <class Grid, class Test> -void reader_test (Test& test) +void reader_test (MPIHelper& mpi, Test& test) { + std::string ext = ".vtu"; for (auto const& test_case : test_cases) { - auto grid = VtkReader<Grid>::read("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"); - test.check(compare_files("reader_writer_test_" + std::get<0>(test_case) + ".vtu", - "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{factory.createGrid()}; + pieces1 = reader.pieces(); + + 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 = reader2.pieces(); + } + + 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])); } } @@ -94,13 +130,17 @@ using int_ = std::integral_constant<int,I>; int main (int argc, char** argv) { - Dune::MPIHelper::instance(argc, argv); + auto& mpi = Dune::MPIHelper::instance(argc, argv); + if (mpi.size() > 1) { + std::cout << "Parallel VtkReader not yet supported\n"; + return 0; + } TestSuite test{}; -#ifdef HAVE_UG +#if HAVE_UG // Test VtkWriter for UGGrid - Hybrid::forEach(std::make_tuple(int_<2>{}, int_<3>{}), [&test](auto dim) + Hybrid::forEach(std::make_tuple(int_<2>{}, int_<3>{}), [&test,&mpi](auto dim) { using GridType = UGGrid<dim.value>; { @@ -108,11 +148,31 @@ int main (int argc, char** argv) FieldVector<double,dim.value> upperRight; upperRight = 1.0; auto numElements = filledArray<dim.value,unsigned int>(4); auto gridPtr = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft, upperRight, numElements); + gridPtr->loadBalance(); + + writer_test(gridPtr->leafGridView()); + } + + reader_test<GridType>(mpi,test); + }); +#endif + +#if DUNE_VERSION_LT(DUNE_GRID,2,7) && HAVE_DUNE_ALUGRID + // Test VtkWriter for ALUGrid. Currently the 2.7 branch is not working. + Hybrid::forEach(std::make_tuple(int_<2>{}, int_<3>{}), [&test,&mpi](auto dim) + { + using GridType = Dune::ALUGrid<dim.value, dim.value, Dune::simplex, Dune::conforming>; + { + FieldVector<double,dim.value> lowerLeft; lowerLeft = 0.0; + FieldVector<double,dim.value> upperRight; upperRight = 1.0; + auto numElements = filledArray<dim.value,unsigned int>(4); + auto gridPtr = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft, upperRight, numElements); + gridPtr->loadBalance(); writer_test(gridPtr->leafGridView()); } - reader_test<GridType>(test); + reader_test<GridType>(mpi,test); }); #endif @@ -121,7 +181,7 @@ int main (int argc, char** argv) { using GridType = YaspGrid<dim.value>; FieldVector<double,dim.value> upperRight; upperRight = 1.0; - auto numElements = filledArray<dim.value,int>(4); + auto numElements = filledArray<dim.value,int>(8); GridType grid(upperRight, numElements); writer_test(grid.leafGridView()); }); diff --git a/src/vtkreader.cc b/src/vtkreader.cc index d82e96ba9848bf87c0bfaa65851b18dd08074e6f..277e47a0b81b6d5ea77b82d3e25e1c3aafb41064 100644 --- a/src/vtkreader.cc +++ b/src/vtkreader.cc @@ -16,6 +16,8 @@ #include <dune/grid/utility/structuredgridfactory.hh> #include <dune/vtk/vtkreader.hh> +#include <dune/vtk/gridcreators/continuousgridcreator.hh> +#include <dune/vtk/gridcreators/discontinuousgridcreator.hh> #include <dune/vtk/writers/vtkunstructuredgridwriter.hh> using namespace Dune; @@ -72,7 +74,7 @@ int main(int argc, char** argv) } { - auto gridPtr = VtkReader<GridType,ConnectedGridCreator>::read("test_ascii_float32.vtu"); + auto gridPtr = VtkReader<GridType,DiscontinuousGridCreator<GridType>>::read("test_ascii_float32.vtu"); auto& grid = *gridPtr; VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView(), Vtk::ASCII); diff --git a/src/vtkwriter.cc b/src/vtkwriter.cc index d44255f424eb2236ee03e1213293a6b110e42fef..435c3ff1a63bb463c45c9ddf70a3227a328b8809 100644 --- a/src/vtkwriter.cc +++ b/src/vtkwriter.cc @@ -38,7 +38,7 @@ static TestCases test_cases = { template <class GridView> void write (std::string prefix, GridView const& gridView) { -#if DUNE_VERSION_LT(DUNE_FUNCTIONS, 2, 7) +#if DUNE_VERSION_LT(DUNE_FUNCTIONS, 2, 6) using namespace BasisBuilder; #else using namespace BasisFactory;