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

[bboxtree] Extract intersection algorithms from bounding box tree

* Generalize bounding box tree
* Extract intersection algorithms
* Fix some epsilon issues in interval (1d) and triangle (2d) point intersection
* Adapt point sources to new interface
* Adapt fv grid geometry to new interface -> element map moved to basefvgeometry
* Improve bounding box tree test
parent e1b0a997
This diff is collapsed.
......@@ -18,33 +18,68 @@
*****************************************************************************/
/*!
* \file
* \brief Base class for the finite volume geometry vector for cell-centered TPFA models
* This builds up the sub control volumes and sub control volume faces
* for each element.
* \brief A map from indices to entities using grid entity seeds
*/
#ifndef DUMUX_ELEMENT_MAP_HH
#define DUMUX_ELEMENT_MAP_HH
#ifndef DUMUX_ENTITY_INDEX_MAP_HH
#define DUMUX_ENTITY_INDEX_MAP_HH
namespace Dumux {
//! An index to element map
template <class GridView>
class ElementMap
: public std::vector<typename GridView::Traits::Grid::template Codim<0>::EntitySeed>
/*!
* \brief A map from indices to entities using grid entity seeds
*/
template <class GridView, int codim = 0>
class EntityMap
{
using Grid = typename GridView::Traits::Grid;
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
public:
ElementMap(const GridView& gridView_) : grid_(gridView_.grid()) {}
using Grid = typename GridView::Traits::Grid;
using Entity = typename Grid::template Codim<codim>::Entity;
using EntitySeed = typename Grid::template Codim<codim>::EntitySeed;
//! constructor moving a ready seed list in here
EntityMap(const Grid& grid, std::vector<EntitySeed>&& seeds)
: grid_(grid)
, seeds_(std::move(seeds))
{}
Element element(IndexType eIdx) const
{ return grid_.entity((*this)[eIdx]); }
//! constructor with all entites of codim
template<class Mapper>
EntityMap(const Grid& grid, const Mapper& mapper)
: grid_(grid)
{
update(mapper);
}
//! update the map after the grid changed
void update(std::vector<EntitySeed>&& seeds)
{ seeds_.swap(std::move(seeds)); }
//! update the map after the grid changed
template<class Mapper>
void update(const Mapper& mapper)
{
const auto& gv = grid_.leafGridView();
seeds_.resize(gv.size(codim));
for (const auto& entity : entities(gv, Dune::Codim<codim>()))
seeds_[mapper.index(entity)] = entity.seed();
}
//! get an element from an index i
Entity operator[](std::size_t i) const
{ return grid_.entity(seeds_[i]); }
//! get the size of the map
std::size_t size() const
{ return seeds_.size(); }
private:
const Grid& grid_;
std::vector<EntitySeed> seeds_;
};
template<class GridView>
using ElementMap = EntityMap<GridView, 0>;
} // end namespace Dumux
#endif
\ No newline at end of file
#endif
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief An axis-aligned bounding box volume hierarchy for dune grids
*
* Dumux implementation of an AABB tree
* Inspired by the AABB tree implementation in DOLFIN by Anders Logg which has the
* following license info: DOLFIN is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*/
#ifndef DUMUX_BOUNDINGBOXTREE_HH
#define DUMUX_BOUNDINGBOXTREE_HH
#include <vector>
#include <array>
#include <algorithm>
#include <numeric>
#include <type_traits>
#include <dune/common/promotiontraits.hh>
#include <dune/common/timer.hh>
#include <dune/common/fvector.hh>
namespace Dumux
{
/*!
* \brief An axis-aligned bounding box volume tree implementation
*
* The class constructs a hierarchical structure of bounding box volumes around
* grid entities. This class can be used to efficiently compute intersections
* between a grid and other geometrical object. It only implements the intersection
* of two of such bounding box trees, so that two independent grids can be intersected.
* \tparam GeometricEntitySet has the following requirements
* * export dimensionworld, ctype
* * a size() member function returning the number of entities
* * begin() and end() member function returning at least forward iterators to entities
* * an index() method returning a consecutive index given an entity
* * an entity() method returning an entity given the consecutive index
* * entities have the following requirements:
* * a member function geometry() returning a geometry with the member functions
* * corner() and corners() returning global coordinates and number of corners
*/
template <class GeometricEntitySet>
class BoundingBoxTree
{
enum { dimworld = GeometricEntitySet::dimensionworld };
using ctype = typename GeometricEntitySet::ctype;
/*!
* \brief Bounding box node data structure
* leaf nodes are indicated by setting child0 to
* the node itself and child1 to the index of the entity in the bounding box.
*/
struct BoundingBoxNode
{
std::size_t child0;
std::size_t child1;
};
public:
//! the type of entity set this tree was built with
using EntitySet = GeometricEntitySet;
//! Default Constructor
BoundingBoxTree() = default;
//! Constructor with gridView
BoundingBoxTree(std::shared_ptr<const GeometricEntitySet> set)
{ build(set); }
//! Build up bounding box tree for a grid with leafGridView
void build(std::shared_ptr<const GeometricEntitySet> set)
{
// set the pointer to the entity set
entitySet_ = set;
// clear all internal data
boundingBoxNodes_.clear();
boundingBoxCoordinates_.clear();
// start the timer
Dune::Timer timer;
// Create bounding boxes for all elements
const auto numLeaves = set->size();
// reserve enough space for the nodes and the coordinates
const auto numNodes = 2*numLeaves - 1;
boundingBoxNodes_.reserve(numNodes);
boundingBoxCoordinates_.reserve(numNodes*2*dimworld);
// create a vector for leaf boxes (min and max for all dims)
std::vector<ctype> leafBoxes(2*dimworld*numLeaves);
for (const auto& geometricEntity : *set)
computeEntityBoundingBox_(leafBoxes.data() + 2*dimworld*set->index(geometricEntity), geometricEntity);
// create the leaf partition, the set of available indices (to be sorted)
std::vector<std::size_t> leafPartition(numLeaves);
std::iota(leafPartition.begin(), leafPartition.end(), 0);
// Recursively build the bounding box tree
build_(leafBoxes, leafPartition.begin(), leafPartition.end());
// We are done, log output
std::cout << "Computed bounding box tree with " << numBoundingBoxes()
<< " nodes for " << numLeaves << " grid entites in "
<< timer.stop() << " seconds." << std::endl;
}
//! the entity set this tree was built with
const EntitySet& entitySet() const
{ return *entitySet_; }
/////////////////////////////////////////////////////
//! Interface to be used by other bounding box trees
/////////////////////////////////////////////////////
//! Get an existing bounding box for a given node
const BoundingBoxNode& getBoundingBoxNode(std::size_t nodeIdx) const
{ return boundingBoxNodes_[nodeIdx]; }
//! Get an existing bounding box for a given node
const ctype* getBoundingBoxCoordinates(std::size_t nodeIdx) const
{ return boundingBoxCoordinates_.data() + 2*dimworld*nodeIdx; }
//! Get the number of bounding boxes currently in the tree
std::size_t numBoundingBoxes() const
{ return boundingBoxNodes_.size(); }
//! Check whether a bounding box node is a leaf node
//! Leaf nodes have itself as child0
bool isLeaf(const BoundingBoxNode& node, std::size_t nodeIdx) const
{ return node.child0 == nodeIdx; }
private:
//! vector of bounding box nodes
std::vector<BoundingBoxNode> boundingBoxNodes_;
//! vector of bounding box coordinates
std::vector<ctype> boundingBoxCoordinates_;
//! a pointer to the entity set
std::shared_ptr<const EntitySet> entitySet_;
//! Compute the bounding box of a grid entity
template <class Entity>
void computeEntityBoundingBox_(ctype* b, const Entity& entity) const
{
// get the bounding box coordinates
ctype* xMin = b;
ctype* xMax = b + dimworld;
// get mesh entity data
auto geometry = entity.geometry();
// Get coordinates of first vertex
auto corner = geometry.corner(0);
for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
xMin[dimIdx] = xMax[dimIdx] = corner[dimIdx];
// Compute the min and max over the remaining vertices
for (std::size_t cornerIdx = 1; cornerIdx < geometry.corners(); ++cornerIdx)
{
corner = geometry.corner(cornerIdx);
for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
{
using std::max;
using std::min;
xMin[dimIdx] = min(xMin[dimIdx], corner[dimIdx]);
xMax[dimIdx] = max(xMax[dimIdx], corner[dimIdx]);
}
}
}
//! Build bounding box tree for all entities recursively
std::size_t build_(const std::vector<ctype>& leafBoxes,
const std::vector<std::size_t>::iterator& begin,
const std::vector<std::size_t>::iterator& end)
{
assert(begin < end);
// If we reached the end of the recursion, i.e. only a leaf box is left
if (end - begin == 1)
{
// Get the bounding box coordinates for the leaf
const std::size_t leafNodeIdx = *begin;
const auto beginCoords = leafBoxes.begin() + 2*dimworld*leafNodeIdx;
const auto endCoords = beginCoords + 2*dimworld;
// Store the data in the bounding box
// leaf nodes are indicated by setting child0 to
// the node itself and child1 to the index of the entity in the bounding box.
return addBoundingBox_(BoundingBoxNode{numBoundingBoxes(), leafNodeIdx}, beginCoords, endCoords);
}
// Compute the bounding box of all bounding boxes in the range [begin, end]
const auto bCoords = computeBBoxOfBBoxes_(leafBoxes, begin, end);
// sort bounding boxes along the longest axis
const auto axis = computeLongestAxis_(bCoords);
// nth_element sorts the range to make sure that middle points to the coordinate median in axis direction
// this is the most expensive part of the algorithm
auto middle = begin + (end - begin)/2;
std::nth_element(begin, middle, end, [&leafBoxes, axis](std::size_t i, std::size_t j)
{
const ctype* bi = leafBoxes.data() + 2*dimworld*i;
const ctype* bj = leafBoxes.data() + 2*dimworld*j;
return bi[axis] + bi[axis + dimworld] < bj[axis] + bj[axis + dimworld];
});
// split the bounding boxes into two at the middle iterator and call build recursively, each
// call resulting in a new node of this bounding box, i.e. the root will be added at the end of the process.
return addBoundingBox_(BoundingBoxNode{build_(leafBoxes, begin, middle), build_(leafBoxes, middle, end)},
bCoords.begin(), bCoords.end());
}
//! Add a new bounding box to the tree
template <class Iterator>
std::size_t addBoundingBox_(BoundingBoxNode&& node,
const Iterator& coordBegin,
const Iterator& coordEnd)
{
// Add the bounding box
boundingBoxNodes_.emplace_back(node);
// Add the bounding box's coordinates
boundingBoxCoordinates_.insert(boundingBoxCoordinates_.end(), coordBegin, coordEnd);
// return the index of the new node
return boundingBoxNodes_.size() - 1;
}
//! Compute the bounding box of a vector of bounding boxes
std::array<ctype, 2*dimworld>
computeBBoxOfBBoxes_(const std::vector<ctype>& leafBoxes,
const std::vector<std::size_t>::iterator& begin,
const std::vector<std::size_t>::iterator& end)
{
std::array<ctype, 2*dimworld> bBoxCoords;
// copy the iterator and get coordinates of first box
auto it = begin;
const auto* bFirst = leafBoxes.data() + 2*dimworld*(*it);
for (int coordIdx = 0; coordIdx < 2*dimworld; ++coordIdx)
bBoxCoords[coordIdx] = bFirst[coordIdx];
// Compute min and max with the remaining boxes
for (; it != end; ++it)
{
const auto* b = leafBoxes.data() + 2*dimworld*(*it);
for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx)
if (b[coordIdx] < bBoxCoords[coordIdx])
bBoxCoords[coordIdx] = b[coordIdx];
for (int coordIdx = dimworld; coordIdx < 2*dimworld; ++coordIdx)
if (b[coordIdx] > bBoxCoords[coordIdx])
bBoxCoords[coordIdx] = b[coordIdx];
}
return bBoxCoords;
}
//! Compute the bounding box of a vector of bounding boxes
std::size_t computeLongestAxis_(const std::array<ctype, 2*dimworld>& bCoords)
{
std::array<ctype, dimworld> axisLength;
for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx)
axisLength[coordIdx] = bCoords[dimworld + coordIdx] - bCoords[coordIdx];
return std::distance(axisLength.begin(), std::max_element(axisLength.begin(), axisLength.end()));
}
};
/*!
* \brief Check whether a point is intersectin a bounding box
* \param point The point
* \param b Pointer to bounding box coordinates
*/
template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 3, int> = 0>
inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
{
static constexpr ctype eps_ = 1.0e-7;
const ctype eps0 = eps_*(b[3] - b[0]);
const ctype eps1 = eps_*(b[4] - b[1]);
const ctype eps2 = eps_*(b[5] - b[2]);
return (b[0] - eps0 <= point[0] && point[0] <= b[3] + eps0 &&
b[1] - eps1 <= point[1] && point[1] <= b[4] + eps1 &&
b[2] - eps2 <= point[2] && point[2] <= b[5] + eps2);
}
template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 2, int> = 0>
inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
{
static constexpr ctype eps_ = 1.0e-7;
const ctype eps0 = eps_*(b[2] - b[0]);
const ctype eps1 = eps_*(b[3] - b[1]);
return (b[0] - eps0 <= point[0] && point[0] <= b[2] + eps0 &&
b[1] - eps1 <= point[1] && point[1] <= b[3] + eps1);
}
template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 1, int> = 0>
inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
{
static constexpr ctype eps_ = 1.0e-7;
const ctype eps0 = eps_*(b[1] - b[0]);
return b[0] - eps0 <= point[0] && point[0] <= b[1] + eps0;
}
/*!
* \brief Check whether a bounding box is intersecting another bounding box
* \param a Pointer to first bounding box coordinates
* \param b Pointer to second bounding box coordinates
*/
template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 3, int> = 0>
inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
{
using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
static constexpr ctype eps_ = 1.0e-7;
const ctype eps0 = eps_*(b[3] - b[0]);
const ctype eps1 = eps_*(b[4] - b[1]);
const ctype eps2 = eps_*(b[5] - b[2]);
return (b[0] - eps0 <= a[3] && a[0] <= b[3] + eps0 &&
b[1] - eps1 <= a[4] && a[1] <= b[4] + eps1 &&
b[2] - eps2 <= a[5] && a[2] <= b[5] + eps2);
}
template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 2, int> = 0>
inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
{
using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
static constexpr ctype eps_ = 1.0e-7;
const ctype eps0 = eps_*(b[2] - b[0]);
const ctype eps1 = eps_*(b[3] - b[1]);
return (b[0] - eps0 <= a[2] && a[0] <= b[2] + eps0 &&
b[1] - eps1 <= a[3] && a[1] <= b[3] + eps1);
}
template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 1, int> = 0>
inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
{
using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
static constexpr ctype eps_ = 1.0e-1;
const ctype eps0 = eps_*(b[1] - b[0]);
return b[0] - eps0 <= a[1] && a[0] <= b[1] + eps0;
}
} // end namespace Dumux
#endif
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief A class storing intersections from intersecting two bounding box trees
*/
#ifndef DUMUX_BOUNDING_BOX_TREE_INTERSECTION_HH
#define DUMUX_BOUNDING_BOX_TREE_INTERSECTION_HH
#include <dune/common/fvector.hh>
#include <dune/common/promotiontraits.hh>
namespace Dumux
{
/*!
* \brief An intersection object resulting from the intersection of two bounding box tree primitives
*
* After is has been found that two leaf bounding boxes intersect a primitive test has to be
* performed to see if the actual entities inside the bounding box intersect too. The result
* if such an intersection is found as an object of this class containing the indices of the
* intersecting entities and the corners of the intersection object.
*/
template<class EntitySet0, class EntitySet1>
class BoundingBoxTreeIntersection
{
enum { dimworld = EntitySet0::dimensionworld };
using ctype = typename Dune::PromotionTraits<typename EntitySet0::ctype, typename EntitySet1::ctype>::PromotedType;
using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
public:
explicit BoundingBoxTreeIntersection(std::size_t a,
std::size_t b,
std::vector<GlobalPosition>&& c)
: a_(a)
, b_(b)
, corners_(std::move(c))
{
static_assert(int(EntitySet0::dimensionworld) == int(EntitySet1::dimensionworld),
"Can only store intersections of entity sets with the same world dimension");
}
//! Get the index of the intersecting entity belonging to this grid
std::size_t first() const
{ return a_; }
//! Get the index of the intersecting entity belonging to the other grid
std::size_t second() const
{ return b_; }
//! Get the corners of the intersection geometry
std::vector<GlobalPosition> corners() const
{ return corners_; }
/*!
* \brief Check if the corners of this intersection match with the given corners
* \note This is useful to check if the intersection geometry of two intersections coincide.
*/
bool cornersMatch(const std::vector<GlobalPosition>& otherCorners) const
{
if (otherCorners.size() != corners_.size())
return false;
const auto eps = 1.5e-7*(corners_[1] - corners_[0]).two_norm();
for (int i = 0; i < corners_.size(); ++i)
if ((corners_[i] - otherCorners[i]).two_norm() > eps)
return false;
return true;
}
private:
std::size_t a_, b_; //!< Indices of the intersection elements
std::vector<GlobalPosition> corners_; //!< the corner points of the intersection geometry
};
} // end namespace Dumux
#endif
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \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
*/
#ifndef DUMUX_GRIDVIEW_GEOMETRIC_ENTITY_SET_HH
#define DUMUX_GRIDVIEW_GEOMETRIC_ENTITY_SET_HH
#include <memory>
#include <dune/grid/common/mcmgmapper.hh>
#include <dumux/common/entitymap.hh>
namespace Dumux {
template <class GridView, int codim = 0>
class GridViewGeometricEntitySet
{
using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
using Entity = typename GridView::template Codim<codim>::Entity;