From e81a6bd859127263a0fe875ef2cf1176089e2622 Mon Sep 17 00:00:00 2001 From: Kilian Date: Tue, 8 Oct 2019 18:01:39 +0200 Subject: [PATCH 1/2] [common][geometricentityset] Add entity set for Dune geometries --- dumux/common/geometry/geometricentityset.hh | 132 +++++++++++++++++++- 1 file changed, 129 insertions(+), 3 deletions(-) diff --git a/dumux/common/geometry/geometricentityset.hh b/dumux/common/geometry/geometricentityset.hh index 0d75999619..012ace5fcc 100644 --- a/dumux/common/geometry/geometricentityset.hh +++ b/dumux/common/geometry/geometricentityset.hh @@ -21,18 +21,19 @@ * \note This can be used e.g. to contruct a bounding box volume hierarchy of a grid * It defines the minimum requirement for such a set */ -#ifndef DUMUX_GRIDVIEW_GEOMETRIC_ENTITY_SET_HH -#define DUMUX_GRIDVIEW_GEOMETRIC_ENTITY_SET_HH +#ifndef DUMUX_GEOMETRIC_ENTITY_SET_HH +#define DUMUX_GEOMETRIC_ENTITY_SET_HH #include #include +#include #include namespace Dumux { /*! * \ingroup Geometry - * \brief An interface for a set of geometric entities + * \brief An interface for a set of geometric entities based on a GridView * \note This can be used e.g. to contruct a bounding box volume hierarchy of a grid * It defines the minimum requirement for such a set */ @@ -108,6 +109,131 @@ private: }; +/*! + * \ingroup Geometry + * \brief An interface for a set of geometric entities + * \note This can be used e.g. to contruct a bounding box volume hierarchy of a grid + * It defines the minimum requirement for such a set + */ +template +class GeometriesEntitySet +{ + /*! + * \brief Wrapper to turn a geometry into a geometric entity + */ + class EntityWrapper + { + public: + using Geometry = GeoType; + + /*! + * \brief Constructor + */ + EntityWrapper(const Geometry& geo, const std::size_t index) : geo_(geo), index_(index) {} + + /*! + * \brief Constructor + */ + EntityWrapper(Geometry&& geo, const std::size_t index) : geo_(std::move(geo)), index_(index) {} + + /*! + * \brief Returns the geometry + */ + const Geometry& geometry() const + { return geo_; } + + /*! + * \brief Returns the index of the geometry + */ + std::size_t index() const + { return index_; } + + private: + Geometry geo_; + std::size_t index_; + }; + +public: + using Entity = EntityWrapper; + + /*! + * \brief Constructor for initializer_list + */ + GeometriesEntitySet(std::initializer_list&& geometries) + { + std::size_t index = 0; + // note: std::initializer_list::begin() returns const T*, + // thus no moving will be performed and only the copying ctor of + // EntityWrapper can be called + for (auto&& g : geometries) + entities_.emplace_back(g, index++); + } + + /*! + * \brief Constructor for a vector of geometries + */ + GeometriesEntitySet(const std::vector& geometries) + { + std::size_t index = 0; + for (auto&& g : geometries) + entities_.emplace_back(g, index++); + } + + /*! + * \brief Constructor for a vector of geometries + */ + GeometriesEntitySet(std::vector&& geometries) + { + std::size_t index = 0; + for (auto&& g : geometries) + entities_.emplace_back(std::move(g), index++); + } + + /*! + * \brief The world dimension of the entity set + */ + enum { dimensionworld = Entity::Geometry::coorddimension }; + + /*! + * \brief the coordinate type + */ + using ctype = typename Entity::Geometry::ctype; + + /*! + * \brief the number of entities in this set + */ + decltype(auto) size() const + { return entities_.size(); } + + /*! + * \brief begin iterator to enable range-based for iteration + */ + decltype(auto) begin() const + { return entities_.begin(); } + + /*! + * \brief end iterator to enable range-based for iteration + */ + decltype(auto) end() const + { return entities_.end(); } + + /*! + * \brief get an entities index + */ + template + std::size_t index(const Entity& e) const + { return e.index(); } + + /*! + * \brief get an entity from an index + */ + Entity entity(std::size_t index) const + { return entities_[index]; } + +private: + std::vector entities_; +}; + } // end namespace Dumux #endif -- GitLab From 6c80d14bc2f4a2f2e39e7304e9ed7c3b99b65cb4 Mon Sep 17 00:00:00 2001 From: Kilian Date: Mon, 14 Oct 2019 15:16:01 +0200 Subject: [PATCH 2/2] [test][bBoxTree] Add test for intersection with geometries --- test/common/boundingboxtree/test_bboxtree.cc | 126 ++++++++++++++++++- 1 file changed, 124 insertions(+), 2 deletions(-) diff --git a/test/common/boundingboxtree/test_bboxtree.cc b/test/common/boundingboxtree/test_bboxtree.cc index 007eca0625..9c6066f7f2 100644 --- a/test/common/boundingboxtree/test_bboxtree.cc +++ b/test/common/boundingboxtree/test_bboxtree.cc @@ -83,8 +83,8 @@ public: return 0; } - template - int intersectTree(const Dumux::BoundingBoxTree& otherTree, + template + int intersectTree(const Dumux::BoundingBoxTree>& otherTree, const OtherGridView& otherGridView, std::size_t expectedUniqueIntersections, bool checkTotalIntersections = false, @@ -142,6 +142,37 @@ public: return 0; } + template + int intersectTree(const Dumux::BoundingBoxTree>& otherTree, + const ExpectedIntersectingElements& expectedIntersectingElements) + { + Dune::Timer timer; + const auto intersections = intersectingEntities(otherTree, *tree_); + std::cout << "Computed " << intersections.size() << " tree intersections in " << timer.elapsed() << std::endl; + + std::unordered_map> sortedResults; + + for (const auto& i : intersections) + { + sortedResults[i.first()].push_back(i.second()); + std::sort(sortedResults[i.first()].begin(), sortedResults[i.first()].end()); + sortedResults[i.first()].erase(std::unique(sortedResults[i.first()].begin(), sortedResults[i.first()].end()), sortedResults[i.first()].end()); + } + + for (int i = 0; i < expectedIntersectingElements.size(); ++i) + { + if (expectedIntersectingElements.at(i) != sortedResults[i]) + { + std::cout << "Error for geometry type " << otherTree.entitySet().entity(i).geometry().type() << std::endl; + for (int j = 0; j < expectedIntersectingElements.at(i).size(); ++j) + if (expectedIntersectingElements.at(i)[j] != sortedResults[i][j]) + std::cout << "expected index " << expectedIntersectingElements.at(i)[j] << ", got " << sortedResults[i][j] << std::endl; + return 1; + } + } + return 0; + } + private: std::shared_ptr tree_; }; @@ -198,6 +229,97 @@ int main (int argc, char *argv[]) try Dune::AffineGeometry geometry(Dune::GeometryTypes::simplex(2), corners); returns.push_back(test.intersectGeometry(geometry, 2145)); // (33*33/2 - 33/2)*4 + 33 #endif + + // test intersection of grid with 1D geometries (lines) + if (dimWorld > 1) + { + using GeometryType = Dune::MultiLinearGeometry; + const std::vector cornersDiagonal{lowerLeft, upperRight}; + std::vector cornersVertical{GlobalPosition(0.4*scaling), GlobalPosition(0.4*scaling)}; + cornersVertical[1][dimWorld-1] = 0.6*scaling; + GeometryType diagonalLine(Dune::GeometryTypes::line, cornersDiagonal); + GeometryType verticalLine(Dune::GeometryTypes::line, cornersVertical); + std::vector geometries{std::move(diagonalLine), std::move(verticalLine)}; + using GeometriesEntitySet = Dumux::GeometriesEntitySet; + GeometriesEntitySet entitySet(std::move(geometries)); + Dumux::BoundingBoxTree geometriesTree(std::make_shared(entitySet)); + std::unordered_map> expectedIntersectingElements; + + if (dimWorld == 2) + { + expectedIntersectingElements[0] = {0, 34, 68, 102, 136, 170, 204, 238, 272, 306, 340, 374, 408, + 442, 476, 510, 544, 578, 612, 646, 680, 714, 748, 782, 816, + 850, 884, 918, 952, 986, 1020, 1054, 1088}; + expectedIntersectingElements[1] = {442, 475, 508, 541, 574, 607, 640}; + returns.push_back(test.intersectTree(geometriesTree, expectedIntersectingElements)); + } + else + { + expectedIntersectingElements[0] = {0, 1123, 2246, 3369, 4492, 5615, 6738, 7861, 8984, 10107, 11230, + 12353, 13476, 14599, 15722, 16845, 17968, 19091, 20214, 21337, 22460, + 23583, 24706, 25829, 26952, 28075, 29198, 30321, 31444, 32567, + 33690, 34813, 35936}; + expectedIntersectingElements[1] = {14599 , 15688, 16777, 17866, 18955, 20044, 21133}; + returns.push_back(test.intersectTree(geometriesTree, expectedIntersectingElements)); + } + } + + // test intersection with 2D geometry (triangle) + if (dimWorld > 1) + { + using GeometryType = Dune::MultiLinearGeometry; + const GlobalPosition corner0 = lowerLeft; + GlobalPosition corner1 = upperRight; corner1*= 0.1; + GlobalPosition corner2 = corner1; corner2[dimWorld-1] = lowerLeft[dimWorld-1]; + std::vector cornersTriangle{corner0, corner1, corner2}; + const GeometryType triangle(Dune::GeometryTypes::simplex(2), cornersTriangle); + using GeometriesEntitySet = Dumux::GeometriesEntitySet; + GeometriesEntitySet entitySet({triangle}); + Dumux::BoundingBoxTree geometriesTree(std::make_shared(entitySet)); + std::unordered_map> expectedIntersectingElements; + + if (dimWorld == 2) + { + expectedIntersectingElements[0] = {0, 1, 2, 3, 34, 35, 36, 68, 69, 102}; + returns.push_back(test.intersectTree(geometriesTree, expectedIntersectingElements)); + } + else + { + expectedIntersectingElements[0] = {0, 34, 68, 102, 1123, 1157, 1191, 2246, 2280, 3369}; + returns.push_back(test.intersectTree(geometriesTree, expectedIntersectingElements)); + } + } + + // test intersection with 2D geometry (quadrilateral) + if (dimWorld > 1) + { + using GeometryType = Dune::AxisAlignedCubeGeometry; + GlobalPosition lowerLeftCube = upperRight; lowerLeftCube *= 0.4; + GlobalPosition upperRightCube = lowerLeftCube; upperRightCube[0] += 0.2*scaling; upperRightCube[1] += 0.2*scaling; + GeometryType cube(lowerLeftCube, upperRightCube); + using GeometriesEntitySet = Dumux::GeometriesEntitySet; + GeometriesEntitySet entitySet(std::vector{cube}); + Dumux::BoundingBoxTree geometriesTree(std::make_shared(entitySet)); + std::unordered_map> expectedIntersectingElements; + + if (dimWorld == 2) + { + expectedIntersectingElements[0] = {442, 443, 444, 445, 446, 447, 448, 475, 476, 477, 478, 479, 480, 481, + 508, 509, 510, 511, 512, 513, 514, 541, 542, 543, 544, 545, 546, 547, + 574, 575, 576, 577, 578, 579, 580, 607, 608, 609, 610, 611, 612, 613, + 640, 641, 642, 643, 644, 645, 646}; + returns.push_back(test.intersectTree(geometriesTree, expectedIntersectingElements)); + } + else + { + expectedIntersectingElements[0] = {14599, 14600, 14601, 14602, 14603, 14604, 14605, 14632, 14633, 14634, + 14635, 14636, 14637, 14638, 14665, 14666, 14667, 14668, 14669,14670, + 14671, 14698, 14699, 14700, 14701, 14702, 14703, 14704, 14731, 14732, + 14733, 14734, 14735, 14736, 14737, 14764, 14765, 14766, 14767,14768, + 14769, 14770, 14797, 14798, 14799, 14800, 14801, 14802, 14803}; + returns.push_back(test.intersectTree(geometriesTree, expectedIntersectingElements)); + } + } } #if HAVE_DUNE_FOAMGRID && WORLD_DIMENSION == 3 -- GitLab