diff --git a/dune/vtk/datacollectors/continuousdatacollector.hh b/dune/vtk/datacollectors/continuousdatacollector.hh index 340f558a545a39d926eea9ff4762dd5e8e566fb1..82b726645a06b86a65681f525aef33b3937e89c9 100644 --- a/dune/vtk/datacollectors/continuousdatacollector.hh +++ b/dune/vtk/datacollectors/continuousdatacollector.hh @@ -63,7 +63,7 @@ public: std::int64_t old_o = 0; for (auto const& c : elements(gridView_, Partitions::all)) { Vtk::CellType cellType(c.type()); - for (int j = 0; j < c.subEntities(dim); ++j) + for (unsigned int j = 0; j < c.subEntities(dim); ++j) cells.connectivity.push_back( std::int64_t(indexSet.subIndex(c,cellType.permutation(j),dim)) ); cells.offsets.push_back(old_o += c.subEntities(dim)); cells.types.push_back(cellType.type()); @@ -83,7 +83,7 @@ public: localFct.bind(e); Vtk::CellType cellType{e.type()}; auto refElem = referenceElement(e.geometry()); - for (int j = 0; j < e.subEntities(dim); ++j) { + for (unsigned int j = 0; j < e.subEntities(dim); ++j) { std::size_t idx = fct.ncomps() * indexSet.subIndex(e,cellType.permutation(j),dim); for (int comp = 0; comp < fct.ncomps(); ++comp) data[idx + comp] = T(localFct.evaluate(comp, refElem.position(cellType.permutation(j),dim))); diff --git a/dune/vtk/datacollectors/discontinuousdatacollector.hh b/dune/vtk/datacollectors/discontinuousdatacollector.hh index 88d60a3b5aca50bff8d484e67289feec8d843bc3..f922eb97569aedd969ea07bafa2a6a7be1e18a3c 100644 --- a/dune/vtk/datacollectors/discontinuousdatacollector.hh +++ b/dune/vtk/datacollectors/discontinuousdatacollector.hh @@ -32,7 +32,7 @@ public: auto const& indexSet = gridView_.indexSet(); for (auto const& c : elements(gridView_, Partitions::interior)) { numPoints_ += c.subEntities(dim); - for (int i = 0; i < c.subEntities(dim); ++i) + for (unsigned int i = 0; i < c.subEntities(dim); ++i) indexMap_[indexSet.subIndex(c, i, dim)] = vertex_idx++; } } @@ -50,7 +50,7 @@ public: std::vector<T> data(this->numPoints() * 3); auto const& indexSet = gridView_.indexSet(); for (auto const& element : elements(gridView_, Partitions::interior)) { - for (int i = 0; i < element.subEntities(dim); ++i) { + for (unsigned int i = 0; i < element.subEntities(dim); ++i) { std::size_t idx = 3 * indexMap_[indexSet.subIndex(element, i, dim)]; auto v = element.geometry().corner(i); for (std::size_t j = 0; j < v.size(); ++j) @@ -74,7 +74,7 @@ public: auto const& indexSet = gridView_.indexSet(); for (auto const& c : elements(gridView_, Partitions::interior)) { Vtk::CellType cellType(c.type()); - for (int j = 0; j < c.subEntities(dim); ++j) { + for (unsigned int j = 0; j < c.subEntities(dim); ++j) { std::int64_t vertex_idx = indexMap_[indexSet.subIndex(c,cellType.permutation(j),dim)]; cells.connectivity.push_back(vertex_idx); } @@ -96,7 +96,7 @@ public: localFct.bind(e); Vtk::CellType cellType{e.type()}; auto refElem = referenceElement(e.geometry()); - for (int j = 0; j < e.subEntities(dim); ++j) { + for (unsigned int j = 0; j < e.subEntities(dim); ++j) { std::size_t idx = fct.ncomps() * indexMap_[indexSet.subIndex(e, cellType.permutation(j), dim)]; for (int comp = 0; comp < fct.ncomps(); ++comp) data[idx + comp] = T(localFct.evaluate(comp, refElem.position(cellType.permutation(j),dim))); diff --git a/dune/vtk/datacollectors/quadraticdatacollector.hh b/dune/vtk/datacollectors/quadraticdatacollector.hh index df3aeb3db0bd2aeeb7f17057919284a6dae4e3a8..a39042e8701622ddb19224778222e21cfc983f5d 100644 --- a/dune/vtk/datacollectors/quadraticdatacollector.hh +++ b/dune/vtk/datacollectors/quadraticdatacollector.hh @@ -42,7 +42,7 @@ public: auto refElem = referenceElement<T,dim>(element.type()); // vertices - for (int i = 0; i < element.subEntities(dim); ++i) { + for (unsigned int i = 0; i < element.subEntities(dim); ++i) { std::size_t idx = 3 * indexSet.subIndex(element, i, dim); auto v = geometry.global(refElem.position(i,dim)); for (std::size_t j = 0; j < v.size(); ++j) @@ -51,7 +51,7 @@ public: data[idx + j] = T(0); } // edge centers - for (int i = 0; i < element.subEntities(dim-1); ++i) { + for (unsigned int i = 0; i < element.subEntities(dim-1); ++i) { std::size_t idx = 3 * (indexSet.subIndex(element, i, dim-1) + gridView_.size(dim)); auto v = geometry.global(refElem.position(i,dim-1)); for (std::size_t j = 0; j < v.size(); ++j) @@ -79,12 +79,12 @@ public: auto const& indexSet = gridView_.indexSet(); for (auto const& c : elements(gridView_, Partitions::interior)) { Vtk::CellType cellType(c.type(), Vtk::QUADRATIC); - for (int j = 0; j < c.subEntities(dim); ++j) { + for (unsigned int j = 0; j < c.subEntities(dim); ++j) { int k = cellType.permutation(j); std::int64_t point_idx = indexSet.subIndex(c,k,dim); cells.connectivity.push_back(point_idx); } - for (int j = 0; j < c.subEntities(dim-1); ++j) { + for (unsigned int j = 0; j < c.subEntities(dim-1); ++j) { int k = cellType.permutation(c.subEntities(dim) + j); std::int64_t point_idx = (indexSet.subIndex(c,k,dim-1) + gridView_.size(dim)); cells.connectivity.push_back(point_idx); @@ -107,13 +107,13 @@ public: localFct.bind(e); Vtk::CellType cellType{e.type(), Vtk::QUADRATIC}; auto refElem = referenceElement(e.geometry()); - for (int j = 0; j < e.subEntities(dim); ++j) { + for (unsigned int j = 0; j < e.subEntities(dim); ++j) { int k = cellType.permutation(j); std::size_t idx = fct.ncomps() * indexSet.subIndex(e, k, dim); for (int comp = 0; comp < fct.ncomps(); ++comp) data[idx + comp] = T(localFct.evaluate(comp, refElem.position(k, dim))); } - for (int j = 0; j < e.subEntities(dim-1); ++j) { + for (unsigned int j = 0; j < e.subEntities(dim-1); ++j) { int k = cellType.permutation(e.subEntities(dim) + j); std::size_t idx = fct.ncomps() * (indexSet.subIndex(e, k, dim-1) + gridView_.size(dim)); for (int comp = 0; comp < fct.ncomps(); ++comp) diff --git a/dune/vtk/datacollectors/spdatacollector.hh b/dune/vtk/datacollectors/spdatacollector.hh index f70cf553de1ae9fa6df143d828c2f508268c20c5..60545ecac5489e7835970580825776f74f64739c 100644 --- a/dune/vtk/datacollectors/spdatacollector.hh +++ b/dune/vtk/datacollectors/spdatacollector.hh @@ -25,6 +25,7 @@ public: SPDataCollector (GridView const& gridView) : Super(gridView) , wholeExtent_(filledArray<6,int>(0)) + , extent_(filledArray<6,int>(0)) , origin_(0.0) , spacing_(0.0) {} @@ -34,6 +35,11 @@ public: return wholeExtent_; } + std::array<int, 6> const& extentImpl () const + { + return extent_; + } + auto const& originImpl () const { return origin_; @@ -48,6 +54,7 @@ public: { Super::updateImpl(); + auto const& localMesh = gridView_.impl().gridLevel().localMesh(); auto const& begin = gridView_.impl().gridLevel().globalMesh().begin(); auto const& end = gridView_.impl().gridLevel().globalMesh().end(); auto const& cube = gridView_.grid().domain().cube(); @@ -55,40 +62,13 @@ public: for (int i = 0; i < dim; ++i) { wholeExtent_[2*i] = begin[i]; wholeExtent_[2*i+1] = end[i]; + extent_[2*i] = localMesh.begin()[i]; + extent_[2*i+1] = localMesh.end()[i]; spacing_[i] = cube.width()[i] / (end[i] - begin[i]); origin_[i] = cube.origin()[i]; } } - template <class Writer> - void writeLocalPieceImpl (Writer const& writer) const - { - auto extent = filledArray<6,int>(0); - auto const& localMesh = gridView_.impl().gridLevel().localMesh(); - for (int i = 0; i < dim; ++i) { - extent[2*i] = localMesh.begin()[i]; - extent[2*i+1] = localMesh.end()[i]; - } - writer(extent); - } - - template <class Writer> - void writePiecesImpl (Writer const& writer) const - { - auto extent = filledArray<6,int>(0); - auto const& partitionList = gridView_.impl().gridLevel().decomposition_; - // for (auto const& part : partitionList) - for (std::size_t p = 0; p < partitionList.size(); ++p) - { - auto const& part = partitionList[p]; - for (int i = 0; i < dim; ++i) { - extent[2*i] = part.begin()[i]; - extent[2*i+1] = part.begin()[i] + part.width(i); - } - writer(p, extent, true); - } - } - private: std::array<int, 6> wholeExtent_; std::array<int, 6> extent_; @@ -99,8 +79,8 @@ private: namespace Impl { - template<class GridView, class ct, int dim, template< int > class Ref, class Comm> - struct StructuredDataCollector<GridView, SPGrid<ct,dim,Ref,Comm>> + template <class GridView, class ct, int dim, template< int > class Ref, class Comm> + struct StructuredDataCollectorImpl<GridView, SPGrid<ct,dim,Ref,Comm>> { using type = SPDataCollector<GridView>; }; diff --git a/dune/vtk/datacollectors/structureddatacollector.hh b/dune/vtk/datacollectors/structureddatacollector.hh index 6c3124dcc2b43d213b8d4fad3999764c4506f43b..0f9967abb44d9dce71ba074657367bdb8890e8f0 100644 --- a/dune/vtk/datacollectors/structureddatacollector.hh +++ b/dune/vtk/datacollectors/structureddatacollector.hh @@ -9,11 +9,11 @@ namespace Impl { // Should be specialized for concrete structured grid template <class GridView, class Grid> - struct StructuredDataCollector; + struct StructuredDataCollectorImpl; } template <class GridView> -using StructuredDataCollector = typename Impl::StructuredDataCollector<GridView, typename GridView::Grid>::type; +using StructuredDataCollector = typename Impl::StructuredDataCollectorImpl<GridView, typename GridView::Grid>::type; /// The Interface for structured data-collectors @@ -33,67 +33,94 @@ public: , ghostLevel_(gridView.overlapSize(0)) {} - /// Return number of grid vertices - std::uint64_t numPointsImpl () const - { - return gridView_.size(GridView::dimension); - } - - void updateImpl () + /// Sequence of Index pairs [begin, end) for the cells in each direction + std::array<int, 6> wholeExtent () const { - defaultDataCollector_.update(); + return this->asDerived().wholeExtentImpl(); } - std::array<int, 6> const& wholeExtent () const + /// Sequence of Index pairs [begin, end) for the cells in each direction of the local partition + std::array<int, 6> extent () const { - return this->asDerived().wholeExtentImpl(); + return this->asDerived().extentImpl(); } - FieldVector<ctype, 3> const& origin () const + /// Lower left corner of the grid + FieldVector<ctype, 3> origin () const { return this->asDerived().originImpl(); } - FieldVector<ctype, 3> const& spacing () const + /// Constant grid spacing in each coordinate direction + FieldVector<ctype, 3> spacing () const { return this->asDerived().spacingImpl(); } + /// Call the `writer` with extent template <class Writer> void writeLocalPiece (Writer const& writer) const { this->asDerived().writeLocalPieceImpl(writer); } + /// Call the `writer` with piece number and piece extent template <class Writer> void writePieces (Writer const& writer) const { this->asDerived().writePiecesImpl(writer); } + /// Return the number of overlapping elements int ghostLevel () const { return this->asDerived().ghostLevelImpl(); } - /// Return the coordinates along the ordinates x, y, and z + /// The coordinates defines point coordinates for an extent by specifying the ordinate along each axis. template <class T> std::array<std::vector<T>, 3> coordinates () const { return this->asDerived().template coordinatesImpl<T>(); } - public: + /// Return number of grid vertices + std::uint64_t numPointsImpl () const + { + return gridView_.size(GridView::dimension); + } + + /// \copyref DefaultDataCollector::update. + void updateImpl () + { + defaultDataCollector_.update(); - /// Return the coordinates of all grid vertices in the order given by the indexSet +#if HAVE_MPI + int rank = -1; + int num_ranks = -1; + MPI_Comm_rank(gridView_.comm(), &rank); + MPI_Comm_size(gridView_.comm(), &num_ranks); + + if (rank == 0) { + extents_.resize(num_ranks); + requests_.resize(num_ranks, MPI_REQUEST_NULL); + for (int i = 1; i < num_ranks; ++i) + MPI_Irecv(extents_[i].data(), extents_[i].size(), MPI_INT, i, /*tag=*/6, gridView_.comm(), &requests_[i]); + } + + sendRequest_ = MPI_REQUEST_NULL; +#endif + } + + /// \copydoc DefaultDataCollector::points. template <class T> std::vector<T> pointsImpl () const { return defaultDataCollector_.template points<T>(); } - /// Evaluate the `fct` at the corners of the elements + /// \copydoc DefaultDataCollector::pointData template <class T, class GlobalFunction> std::vector<T> pointDataImpl (GlobalFunction const& fct) const { @@ -101,30 +128,79 @@ public: } -private: + /// Default implementation for \ref writeLocalPiece. Calculates the extent and communicates it to + /// rank 0. + template <class Writer> + void writeLocalPieceImpl (Writer const& writer) const + { + auto&& extent = this->extent(); + +#if HAVE_MPI + int sendFlag = 0; + MPI_Status sendStatus; + MPI_Test(&sendRequest_, &sendFlag, &sendStatus); + + if (sendFlag) { + int rank = -1; + MPI_Comm_rank(gridView_.comm(), &rank); + if (rank != 0) { + MPI_Isend(extent.data(), extent.size(), MPI_INT, 0, /*tag=*/6, gridView_.comm(), &sendRequest_); + } else { + extents_[0] = extent; + } + } +#endif + + writer(extent); + } + + /// Receive extent from all ranks and call the `writer` with the rank's extent vector + template <class Writer> + void writePiecesImpl (Writer const& writer) const + { + writer(0, extents_[0], true); + +#if HAVE_MPI + int num_ranks = -1; + MPI_Comm_size(gridView_.comm(), &num_ranks); + for (int p = 1; p < num_ranks; ++p) { + int idx = -1; + MPI_Status status; + MPI_Waitany(num_ranks, requests_.data(), &idx, &status); + if (idx != MPI_UNDEFINED) { + assert(idx == status.MPI_SOURCE && status.MPI_TAG == 6); + writer(idx, extents_[idx], true); + } + } +#endif + } + + + /// Return the \ref GridView::overlapSize int ghostLevelImpl () const { return ghostLevel_; } + + /// Ordinate along each axis with constant \ref spacing from the \ref origin template <class T> std::array<std::vector<T>, 3> coordinatesImpl () const { auto origin = this->origin(); auto spacing = this->spacing(); + auto extent = this->extent(); std::array<std::vector<T>, 3> ordinates{}; - writeLocalPiece([&ordinates,&origin,&spacing](auto const& extent) { - for (std::size_t d = 0; d < GridView::dimension; ++d) { - auto s = extent[2*d+1] - extent[2*d] + 1; - ordinates[d].resize(s); - for (std::size_t i = 0; i < s; ++i) - ordinates[d][i] = origin[d] + (extent[2*d] + i)*spacing[d]; - } - }); - - for (std::size_t d = GridView::dimension; d < 3; ++d) + for (int d = 0; d < GridView::dimension; ++d) { + auto s = extent[2*d+1] - extent[2*d] + 1; + ordinates[d].resize(s); + for (int i = 0; i < s; ++i) + ordinates[d][i] = origin[d] + (extent[2*d] + i)*spacing[d]; + } + + for (int d = GridView::dimension; d < 3; ++d) ordinates[d].resize(1, T(0)); return ordinates; @@ -133,6 +209,12 @@ private: private: DefaultDataCollector<GridView> defaultDataCollector_; int ghostLevel_; + +#if HAVE_MPI + mutable std::vector<std::array<int,6>> extents_; + mutable std::vector<MPI_Request> requests_; + mutable MPI_Request sendRequest_ = MPI_REQUEST_NULL; +#endif }; }} // end namespace Dune::experimental diff --git a/dune/vtk/datacollectors/yaspdatacollector.hh b/dune/vtk/datacollectors/yaspdatacollector.hh index 7130dc52e482e8ca520de2c452de3630bae06015..7f5b352b83ef2914aa81bd87db22f7297d44fce8 100644 --- a/dune/vtk/datacollectors/yaspdatacollector.hh +++ b/dune/vtk/datacollectors/yaspdatacollector.hh @@ -21,6 +21,7 @@ public: YaspDataCollector (GridView const& gridView) : Super(gridView) , wholeExtent_(filledArray<6,int>(0)) + , extent_(filledArray<6,int>(0)) , origin_(0.0) , spacing_(0.0) , level_(gridView.template begin<0,All_Partition>()->level()) @@ -31,6 +32,11 @@ public: return wholeExtent_; } + std::array<int, 6> const& extentImpl () const + { + return extent_; + } + auto const& originImpl () const { return origin_; @@ -44,34 +50,24 @@ public: void updateImpl () { Super::updateImpl(); - localPieceCalled_ = false; for (int i = 0; i < dim; ++i) { wholeExtent_[2*i] = 0; wholeExtent_[2*i+1] = gridView_.grid().levelSize(level_,i); } + auto const& gl = *gridView_.grid().begin(level_); + auto const& g = gl.interior[0]; + auto const& gc = *g.dataBegin(); + for (int i = 0; i < dim; ++i) { + extent_[2*i] = gc.min(i); + extent_[2*i+1] = gc.max(i)+1; + } + auto it = gridView_.grid().begin(level_); initGeometry(it->coords); - -#if HAVE_MPI - int rank = -1; - int num_ranks = -1; - MPI_Comm_rank(gridView_.comm(), &rank); - MPI_Comm_size(gridView_.comm(), &num_ranks); - - if (rank == 0) { - extents_.resize(num_ranks); - requests_.resize(num_ranks); - for (int i = 1; i < num_ranks; ++i) - MPI_Irecv(extents_[i].data(), extents_[i].size(), MPI_INT, i, /*tag=*/6, gridView_.comm(), &requests_[i]); - } -#endif } - template <class Coords> - void initGeometry (Coords const& coords) { DUNE_THROW(NotImplemented, "Coordinate-Type not implemented!"); } - void initGeometry (EquidistantCoordinates<ctype,dim> const& coords) { for (int i = 0; i < dim; ++i) { @@ -96,51 +92,8 @@ public: } } - template <class Writer> - void writeLocalPieceImpl (Writer const& writer) const - { - auto const& gl = *gridView_.grid().begin(level_); - auto const& g = gl.interior[0]; - auto extent = filledArray<6,int>(0); - auto const& gc = *g.dataBegin(); - for (int i = 0; i < dim; ++i) { - extent[2*i] = gc.min(i); - extent[2*i+1] = gc.max(i)+1; - } - -#if HAVE_MPI - if (!localPieceCalled_) { - int rank = -1; - MPI_Comm_rank(gridView_.comm(), &rank); - if (rank != 0) { - MPI_Isend(extent.data(), extent.size(), MPI_INT, 0, /*tag=*/6, gridView_.comm(), &sendRequest_); - } else { - extents_[0] = extent; - } - localPieceCalled_ = true; - } -#endif - - writer(extent); - } - - template <class Writer> - void writePiecesImpl (Writer const& writer) const - { - assert(localPieceCalled_); -#if HAVE_MPI - int num_ranks = -1; - MPI_Comm_size(gridView_.comm(), &num_ranks); - - std::vector<MPI_Status> status(num_ranks-1); - MPI_Waitall(num_ranks-1, requests_.data()+1, status.data()); - - for (int p = 0; p < num_ranks; ++p) { - writer(p, extents_[p], true); - } -#endif - } + /// Extract the ordinates from the Coordinates object of the current level template <class T> std::array<std::vector<T>, 3> coordinatesImpl () const { @@ -148,17 +101,14 @@ public: auto const& coords = it->coords; std::array<std::vector<T>, 3> ordinates{}; - writeLocalPieceImpl([&ordinates,&coords](auto const& extent) - { - for (std::size_t d = 0; d < dim; ++d) { - auto s = extent[2*d+1] - extent[2*d] + 1; - ordinates[d].resize(s); - for (std::size_t i = 0; i < s; ++i) - ordinates[d][i] = coords.coordinate(d, extent[2*d] + i); - } - }); - - for (std::size_t d = dim; d < 3; ++d) + for (int d = 0; d < dim; ++d) { + auto s = extent_[2*d+1] - extent_[2*d] + 1; + ordinates[d].resize(s); + for (int i = 0; i < s; ++i) + ordinates[d][i] = coords.coordinate(d, extent_[2*d] + i); + } + + for (int d = dim; d < 3; ++d) ordinates[d].resize(1, T(0)); return ordinates; @@ -166,21 +116,17 @@ public: private: std::array<int, 6> wholeExtent_; - FieldVector<ctype,3> spacing_; + std::array<int, 6> extent_; FieldVector<ctype,3> origin_; + FieldVector<ctype,3> spacing_; std::vector<std::size_t> indexMap_; int level_; - - mutable std::vector<std::array<int,6>> extents_; - mutable std::vector<MPI_Request> requests_; - mutable MPI_Request sendRequest_; - mutable bool localPieceCalled_ = false; }; namespace Impl { - template<class GridView, int dim, class Coordinates> - struct StructuredDataCollector<GridView, YaspGrid<dim,Coordinates>> + template <class GridView, int dim, class Coordinates> + struct StructuredDataCollectorImpl<GridView, YaspGrid<dim,Coordinates>> { using type = YaspDataCollector<GridView>; }; diff --git a/dune/vtk/gridcreator.hh b/dune/vtk/gridcreator.hh index ec44cad31d224de92439e9d3da22e202e1531ba3..50504b8bb94eef6b537555bf3a61a596ebbdbf0a 100644 --- a/dune/vtk/gridcreator.hh +++ b/dune/vtk/gridcreator.hh @@ -32,10 +32,10 @@ namespace Dune { namespace experimental Vtk::CellType cellType{type}; auto refElem = referenceElement<double,Grid::dimension>(type); - std::size_t nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]); + 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 (std::size_t j = 0; j < nNodes; ++j) + for (int j = 0; j < nNodes; ++j) vtk_cell.push_back( connectivity[idx++] ); if (cellType.noPermutation()) @@ -43,7 +43,7 @@ namespace Dune { namespace experimental else { // apply index permutation std::vector<unsigned int> cell(nNodes); - for (std::size_t j = 0; j < nNodes; ++j) + for (int j = 0; j < nNodes; ++j) cell[j] = vtk_cell[cellType.permutation(j)]; factory.insertElement(type,cell); @@ -93,10 +93,10 @@ namespace Dune { namespace experimental auto type = Vtk::to_geometry(types[i]); Vtk::CellType cellType{type}; - std::size_t nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]); + int nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]); assert(nNodes > 0); std::vector<unsigned int> vtk_cell; vtk_cell.reserve(nNodes); - for (std::size_t j = 0; j < nNodes; ++j) { + 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); @@ -107,7 +107,7 @@ namespace Dune { namespace experimental else { // apply index permutation std::vector<unsigned int> cell(nNodes); - for (std::size_t j = 0; j < nNodes; ++j) + for (int j = 0; j < nNodes; ++j) cell[j] = vtk_cell[cellType.permutation(j)]; factory.insertElement(type,cell); diff --git a/dune/vtk/vtkreader.impl.hh b/dune/vtk/vtkreader.impl.hh index 0c735cbed148cf287ef27dd9d9004786b4256e5d..320fae1c52557c989c9b2057b88b9ca56910b374 100644 --- a/dune/vtk/vtkreader.impl.hh +++ b/dune/vtk/vtkreader.impl.hh @@ -344,7 +344,7 @@ 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() == vec_offsets.back()); + assert(vec_connectivity.size() == std::size_t(vec_offsets.back())); } @@ -368,7 +368,7 @@ void read_compressed (T* buffer, unsigned char* buffer_in, Bytef* uncompressed_buffer = reinterpret_cast<Bytef*>(buffer); input.read((char*)(compressed_buffer), compressed_space); - assert(input.gcount() == compressed_space); + assert(uLongf(input.gcount()) == compressed_space); if (uncompress(uncompressed_buffer, &uncompressed_space, compressed_buffer, compressed_space) != Z_OK) { std::cerr << "Zlib error while uncompressing data.\n"; @@ -429,7 +429,7 @@ void VtkReader<Grid,Creator>::readAppended (std::ifstream& input, std::vector<T> } } else { input.read((char*)(values.data()), size); - assert(input.gcount() == size); + assert(input.gcount() == std::streamsize(size)); } } diff --git a/dune/vtk/writers/vtkimagedatawriter.impl.hh b/dune/vtk/writers/vtkimagedatawriter.impl.hh index c87389fd0d815b3180ef6eb1c31bc6bd01c51250..caff7e56de2b8d9190140463762e3040d651f5d8 100644 --- a/dune/vtk/writers/vtkimagedatawriter.impl.hh +++ b/dune/vtk/writers/vtkimagedatawriter.impl.hh @@ -103,7 +103,7 @@ void VtkImageDataWriter<GV,DC> template <class GV, class DC> void VtkImageDataWriter<GV,DC> - ::writeParallelFile (std::string const& pfilename, int size) const + ::writeParallelFile (std::string const& pfilename, int /*size*/) const { std::string filename = pfilename + ".p" + this->fileExtension(); std::ofstream out(filename, std::ios_base::ate | std::ios::binary); diff --git a/dune/vtk/writers/vtkrectilineargridwriter.impl.hh b/dune/vtk/writers/vtkrectilineargridwriter.impl.hh index 90ef0f87c53662ae1e59a7ed93b7992ab49fec17..6d37be65518708621946bd873e41b6ddc3417fe1 100644 --- a/dune/vtk/writers/vtkrectilineargridwriter.impl.hh +++ b/dune/vtk/writers/vtkrectilineargridwriter.impl.hh @@ -39,8 +39,6 @@ void VtkRectilinearGridWriter<GV,DC> << ">\n"; auto const& wholeExtent = dataCollector_.wholeExtent(); - auto const& origin = dataCollector_.origin(); - auto const& spacing = dataCollector_.spacing(); out << "<RectilinearGrid" << " WholeExtent=\"" << join(wholeExtent.begin(), wholeExtent.end()) << "\"" << ">\n"; @@ -113,7 +111,7 @@ void VtkRectilinearGridWriter<GV,DC> template <class GV, class DC> void VtkRectilinearGridWriter<GV,DC> - ::writeParallelFile (std::string const& pfilename, int size) const + ::writeParallelFile (std::string const& pfilename, int /*size*/) const { std::string filename = pfilename + ".p" + this->fileExtension(); std::ofstream out(filename, std::ios_base::ate | std::ios::binary); @@ -128,8 +126,6 @@ void VtkRectilinearGridWriter<GV,DC> << ">\n"; auto const& wholeExtent = dataCollector_.wholeExtent(); - auto const& origin = dataCollector_.origin(); - auto const& spacing = dataCollector_.spacing(); out << "<PRectilinearGrid" << " GhostLevel=\"" << dataCollector_.ghostLevel() << "\"" << " WholeExtent=\"" << join(wholeExtent.begin(), wholeExtent.end()) << "\"" diff --git a/dune/vtk/writers/vtkstructuredgridwriter.impl.hh b/dune/vtk/writers/vtkstructuredgridwriter.impl.hh index eabb07ffcd55103c2882c919b2012baf42fe423e..1cab9b55dd27d38fc4c32c195fa656a4793c57e9 100644 --- a/dune/vtk/writers/vtkstructuredgridwriter.impl.hh +++ b/dune/vtk/writers/vtkstructuredgridwriter.impl.hh @@ -106,7 +106,7 @@ void VtkStructuredGridWriter<GV,DC> template <class GV, class DC> void VtkStructuredGridWriter<GV,DC> - ::writeParallelFile (std::string const& pfilename, int size) const + ::writeParallelFile (std::string const& pfilename, int /*size*/) const { std::string filename = pfilename + ".p" + this->fileExtension(); std::ofstream out(filename, std::ios_base::ate | std::ios::binary); diff --git a/dune/vtk/writers/vtkunstructuredgridwriter.hh b/dune/vtk/writers/vtkunstructuredgridwriter.hh index 1ec340544e533f423af69c9c3390425bc4e55a7c..44a4aa80af1c7249affd005f323d6783f27334ca 100644 --- a/dune/vtk/writers/vtkunstructuredgridwriter.hh +++ b/dune/vtk/writers/vtkunstructuredgridwriter.hh @@ -44,6 +44,16 @@ namespace Dune { namespace experimental return "vtu"; } + // 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, + std::vector<pos_type>& offsets) const; + + // Collect element connectivity, offsets and element types, and pass the + // resulting vectors to \ref writeAppended. + std::array<std::uint64_t,3> writeCellsAppended (std::ofstream& out) 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 f62c591cce90709e2c1082d5b5179ca7d4e00f1a..4d34751ad9385cca35913a54f6278c8219dbc991 100644 --- a/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh +++ b/dune/vtk/writers/vtkunstructuredgridwriter.impl.hh @@ -64,7 +64,7 @@ void VtkUnstructuredGridWriter<GV,DC> // Write element connectivity, types and offsets out << "<Cells>\n"; - this->writeCells(out, offsets); + writeCells(out, offsets); out << "</Cells>\n"; out << "</Piece>\n"; @@ -92,7 +92,8 @@ void VtkUnstructuredGridWriter<GV,DC> blocks.push_back( this->template writePointsAppended<float>(out) ); else blocks.push_back( this->template writePointsAppended<double>(out) ); - auto bs = this->writeCellsAppended(out); + + auto bs = writeCellsAppended(out); blocks.insert(blocks.end(), bs.begin(), bs.end()); out << "</AppendedData>\n"; } @@ -169,4 +170,68 @@ void VtkUnstructuredGridWriter<GV,DC> out << "</VTKFile>"; } + +template <class GV, class DC> +void VtkUnstructuredGridWriter<GV,DC> + ::writeCells (std::ofstream& out, std::vector<pos_type>& offsets) const +{ + if (format_ == Vtk::ASCII) { + auto cells = dataCollector_.cells(); + out << "<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">\n"; + std::size_t i = 0; + for (auto const& c : cells.connectivity) + out << c << (++i % 6 != 0 ? ' ' : '\n'); + out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n"; + + out << "<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">\n"; + i = 0; + for (auto const& o : cells.offsets) + out << o << (++i % 6 != 0 ? ' ' : '\n'); + out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n"; + + out << "<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n"; + i = 0; + for (auto const& t : cells.types) + out << int(t) << (++i % 6 != 0 ? ' ' : '\n'); + out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n"; + } + else { // Vtk::APPENDED format + out << "<DataArray type=\"Int64\" Name=\"connectivity\" format=\"appended\""; + out << " offset="; + offsets.push_back(out.tellp()); + out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' '); + out << "/>\n"; + + out << "<DataArray type=\"Int64\" Name=\"offsets\" format=\"appended\""; + out << " offset="; + offsets.push_back(out.tellp()); + out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' '); + out << "/>\n"; + + out << "<DataArray type=\"UInt8\" Name=\"types\" format=\"appended\""; + 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> +std::array<std::uint64_t,3> VtkUnstructuredGridWriter<GV,DC> + ::writeCellsAppended (std::ofstream& out) const +{ + assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n"); + + auto cells = dataCollector_.cells(); + + // write conncetivity, offsets, and types + std::uint64_t bs0 = this->writeAppended(out, cells.connectivity); + std::uint64_t bs1 = this->writeAppended(out, cells.offsets); + std::uint64_t bs2 = this->writeAppended(out, cells.types); + + return {bs0, bs1, bs2}; +} + + }} // end namespace Dune::experimental diff --git a/dune/vtk/writers/vtkwriterinterface.hh b/dune/vtk/writers/vtkwriterinterface.hh index 5929666407176e187298ba633536494f5d606660..caf6e43ca84b70c92e483f720970122a878042d5 100644 --- a/dune/vtk/writers/vtkwriterinterface.hh +++ b/dune/vtk/writers/vtkwriterinterface.hh @@ -87,65 +87,37 @@ namespace Dune { namespace experimental GlobalFunction const& fct, PositionTypes type) const; - // Write the coordinates of the vertices to the output stream `out`. In case - // of binary format, stores the streampos of XML attributes "offset" in the - // vector `offsets`. - void writePoints (std::ofstream& out, - std::vector<pos_type>& offsets) const; - - // 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, - std::vector<pos_type>& offsets) const; - // Collect point or cell data (depending on \ref PositionTypes) and pass // the resulting vector to \ref writeAppended. template <class T> std::uint64_t writeDataAppended (std::ofstream& out, - GlobalFunction const& localFct, + GlobalFunction const& fct, PositionTypes type) const; + // Write the coordinates of the vertices to the output stream `out`. In case + // of binary format, stores the streampos of XML attributes "offset" in the + // vector `offsets`. + void writePoints (std::ofstream& out, + std::vector<pos_type>& offsets) const; + // Collect point positions and pass the resulting vector to \ref writeAppended. template <class T> std::uint64_t writePointsAppended (std::ofstream& out) const; - // Collect element connectivity, offsets and element types, and pass the - // resulting vectors to \ref writeAppended. - std::array<std::uint64_t,3> writeCellsAppended (std::ofstream& out) const; - // Write the `values` in blocks (possibly compressed) to the output // stream `out`. Return the written block size. template <class T> std::uint64_t writeAppended (std::ofstream& out, std::vector<T> const& values) const; - - Std::optional<std::string> getScalarName (std::vector<GlobalFunction> const& data) const + /// Return PointData/CellData attributes for the name of the first scalar/vector/tensor DataArray + std::string getNames (std::vector<GlobalFunction> const& data) const { auto scalar = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.ncomps() == 1; }); - return scalar != data.end() ? Std::optional<std::string>{scalar->name()} : Std::optional<std::string>{}; - } - - Std::optional<std::string> getVectorName (std::vector<GlobalFunction> const& data) const - { auto vector = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.ncomps() == 3; }); - return vector != data.end() ? Std::optional<std::string>{vector->name()} : Std::optional<std::string>{}; - } - - Std::optional<std::string> getTensorName (std::vector<GlobalFunction> const& data) const - { auto tensor = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.ncomps() == 9; }); - return tensor != data.end() ? Std::optional<std::string>{tensor->name()} : Std::optional<std::string>{}; - } - - std::string getNames (std::vector<GlobalFunction> const& data) const - { - auto n1 = getScalarName(data); - auto n2 = getVectorName(data); - auto n3 = getTensorName(data); - return (n1 ? " Scalars=\"" + *n1 + "\"" : "") - + (n2 ? " Vectors=\"" + *n2 + "\"" : "") - + (n3 ? " Tensors=\"" + *n3 + "\"" : ""); + return (scalar != data.end() ? " Scalars=\"" + scalar->name() + "\"" : "") + + (vector != data.end() ? " Vectors=\"" + vector->name() + "\"" : "") + + (tensor != data.end() ? " Tensors=\"" + tensor->name() + "\"" : ""); } // Returns endianness diff --git a/dune/vtk/writers/vtkwriterinterface.impl.hh b/dune/vtk/writers/vtkwriterinterface.impl.hh index dc16992960399f31da60ee4b696c94a2f63d23dd..5b39fb25fb068424f56e188213817420fecea52d 100644 --- a/dune/vtk/writers/vtkwriterinterface.impl.hh +++ b/dune/vtk/writers/vtkwriterinterface.impl.hh @@ -113,52 +113,6 @@ void VtkWriterInterface<GV,DC> } -template <class GV, class DC> -void VtkWriterInterface<GV,DC> - ::writeCells (std::ofstream& out, std::vector<pos_type>& offsets) const -{ - if (format_ == Vtk::ASCII) { - auto cells = dataCollector_.cells(); - out << "<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">\n"; - std::size_t i = 0; - for (auto const& c : cells.connectivity) - out << c << (++i % 6 != 0 ? ' ' : '\n'); - out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n"; - - out << "<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">\n"; - i = 0; - for (auto const& o : cells.offsets) - out << o << (++i % 6 != 0 ? ' ' : '\n'); - out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n"; - - out << "<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n"; - i = 0; - for (auto const& t : cells.types) - out << int(t) << (++i % 6 != 0 ? ' ' : '\n'); - out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n"; - } - else { // Vtk::APPENDED format - out << "<DataArray type=\"Int64\" Name=\"connectivity\" format=\"appended\""; - out << " offset="; - offsets.push_back(out.tellp()); - out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' '); - out << "/>\n"; - - out << "<DataArray type=\"Int64\" Name=\"offsets\" format=\"appended\""; - out << " offset="; - offsets.push_back(out.tellp()); - out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' '); - out << "/>\n"; - - out << "<DataArray type=\"UInt8\" Name=\"types\" format=\"appended\""; - 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> template <class T> std::uint64_t VtkWriterInterface<GV,DC> @@ -188,23 +142,6 @@ std::uint64_t VtkWriterInterface<GV,DC> } -template <class GV, class DC> -std::array<std::uint64_t,3> VtkWriterInterface<GV,DC> - ::writeCellsAppended (std::ofstream& out) const -{ - assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n"); - - auto cells = dataCollector_.cells(); - - // write conncetivity, offsets, and types - std::uint64_t bs0 = this->writeAppended(out, cells.connectivity); - std::uint64_t bs1 = this->writeAppended(out, cells.offsets); - std::uint64_t bs2 = this->writeAppended(out, cells.types); - - return {bs0, bs1, bs2}; -} - - namespace Impl { template <class T>