Commit d14e1487 authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/geo-ent-set' into 'master'

[common][geometricentityset] Add entity set for Dune geometries

Closes #737

See merge request !1741
parents 9550ab25 6c80d14b
......@@ -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 <memory>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/geometry/multilineargeometry.hh>
#include <dumux/common/entitymap.hh>
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 GeoType>
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<typename Entity::Geometry>&& 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<typename Entity::Geometry>& 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<typename Entity::Geometry>&& 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<class Entity>
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<Entity> entities_;
};
} // end namespace Dumux
#endif
......@@ -83,8 +83,8 @@ public:
return 0;
}
template <class OtherEntitySet, class OtherGridView>
int intersectTree(const Dumux::BoundingBoxTree<OtherEntitySet>& otherTree,
template <class OtherGridView>
int intersectTree(const Dumux::BoundingBoxTree<GridViewGeometricEntitySet<OtherGridView>>& otherTree,
const OtherGridView& otherGridView,
std::size_t expectedUniqueIntersections,
bool checkTotalIntersections = false,
......@@ -142,6 +142,37 @@ public:
return 0;
}
template <class GeometryType, class ExpectedIntersectingElements>
int intersectTree(const Dumux::BoundingBoxTree<GeometriesEntitySet<GeometryType>>& 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<std::size_t, std::vector<std::size_t>> 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<BoundingBoxTree> tree_;
};
......@@ -198,6 +229,97 @@ int main (int argc, char *argv[]) try
Dune::AffineGeometry<double, 2, WORLD_DIMENSION> 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<Scalar, 1, dimWorld>;
const std::vector<GlobalPosition> cornersDiagonal{lowerLeft, upperRight};
std::vector<GlobalPosition> 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<GeometryType> geometries{std::move(diagonalLine), std::move(verticalLine)};
using GeometriesEntitySet = Dumux::GeometriesEntitySet<GeometryType>;
GeometriesEntitySet entitySet(std::move(geometries));
Dumux::BoundingBoxTree<GeometriesEntitySet> geometriesTree(std::make_shared<GeometriesEntitySet>(entitySet));
std::unordered_map<std::size_t, std::vector<std::size_t>> 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<Scalar, 2, dimWorld>;
const GlobalPosition corner0 = lowerLeft;
GlobalPosition corner1 = upperRight; corner1*= 0.1;
GlobalPosition corner2 = corner1; corner2[dimWorld-1] = lowerLeft[dimWorld-1];
std::vector<GlobalPosition> cornersTriangle{corner0, corner1, corner2};
const GeometryType triangle(Dune::GeometryTypes::simplex(2), cornersTriangle);
using GeometriesEntitySet = Dumux::GeometriesEntitySet<GeometryType>;
GeometriesEntitySet entitySet({triangle});
Dumux::BoundingBoxTree<GeometriesEntitySet> geometriesTree(std::make_shared<GeometriesEntitySet>(entitySet));
std::unordered_map<std::size_t, std::vector<std::size_t>> 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<Scalar, 2, dimWorld>;
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<GeometryType>;
GeometriesEntitySet entitySet(std::vector<GeometryType>{cube});
Dumux::BoundingBoxTree<GeometriesEntitySet> geometriesTree(std::make_shared<GeometriesEntitySet>(entitySet));
std::unordered_map<std::size_t, std::vector<std::size_t>> 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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment