Commit 15d9301b authored by Timo Koch's avatar Timo Koch
Browse files

[bboxtree] Improve geometry collision and return intersection corners

The geoemtry collision noe return a vector of intersection corners.
The collision algorithm with another budning box tree returns
a vector of BoundingBoxTreeIntersection object.

The object has a method first and second to get the intersecting
entity indices as well as a method corners. From that you
could build e.g. a intersection geometry. The intersections are
not unique yet. There might be interesctions with the same
corner points but different elements.
parent cb4a9df1
......@@ -733,6 +733,32 @@ public:
};
};
//! An intersection object resulting from the collision of two bounding box tree primitives
template<class GridView1, class GridView2>
class BoundingBoxTreeIntersection
{
const static int dimworld = GridView1::dimensionworld;
using Scalar = typename GridView1::ctype;
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
public:
BoundingBoxTreeIntersection(unsigned int a, unsigned int b, std::vector<GlobalPosition>&& c)
: a_(a), b_(b), corners_(c) {}
inline unsigned int first() const
{ return a_; }
inline unsigned int second() const
{ return b_; }
inline std::vector<GlobalPosition> corners() const
{ return corners_; }
private:
unsigned int a_, b_;
std::vector<GlobalPosition> corners_;
};
//! An index to element map
template <class GridView>
class IndexToElementMap
......@@ -762,22 +788,19 @@ class BoundingBoxTree
static const int dim = GridView::dimension;
static const int dimworld = GridView::dimensionworld;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename Dune::ReferenceElements<double, dim> ReferenceElements;
typedef typename Dune::ReferenceElement<double, dim> ReferenceElement;
using Element = typename GridView::template Codim<0>::Entity;
using Scalar = typename GridView::ctype;
using GlobalPosition = typename Dune::FieldVector<Scalar, dimworld>;
public:
// Default Constructor
//! Default Constructor
BoundingBoxTree() {}
// Constructor with gridView
//! Constructor with gridView
BoundingBoxTree(const GridView& leafGridView)
{ build(leafGridView); }
// Destructor
~BoundingBoxTree() {}
// Build up bounding box tree for a grid
//! Build up bounding box tree for a grid
void build(const GridView& leafGridView)
{
// clear data if any
......@@ -817,7 +840,7 @@ public:
<< timer.stop() << " seconds." << std::endl;
}
// Compute all intersections between entities and a point
//! Compute all intersections between entities and a point
std::vector<unsigned int> computeEntityCollisions(const Dune::FieldVector<double, dimworld>& point) const
{
// Call the recursive find function to find candidates
......@@ -826,29 +849,27 @@ public:
return entities;
}
// Compute all intersections between entities and a point
//! Compute all intersections between entities and a point
template<class OtherGridView>
std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
std::vector<BoundingBoxTreeIntersection<GridView, OtherGridView>>
computeEntityCollisions(const Dumux::BoundingBoxTree<OtherGridView>& otherTree) const
{
// check if the world dimensions match
static_assert(dimworld == OtherGridView::dimensionworld, "Can only collide bounding box trees of same world dimension");
// Create data structure for return type
std::vector<unsigned int> myEntities;
std::vector<unsigned int> otherEntities;
std::vector<BoundingBoxTreeIntersection<GridView, OtherGridView>> intersections;
// Call the recursive find function to find candidates
computeCollisions_(otherTree,
this->numBoundingBoxes_() - 1,
otherTree.numBoundingBoxes_() -1,
myEntities,
otherEntities);
intersections);
return std::make_pair(myEntities, otherEntities);
return intersections;
}
// Get an element from a global element index
//! Get an element from a global element index
Element entity(unsigned int eIdx) const
{ return indexToElementMap_->entity(eIdx); }
......@@ -958,8 +979,7 @@ private:
void computeCollisions_(const Dumux::BoundingBoxTree<OtherGridView>& treeB,
unsigned int nodeA,
unsigned int nodeB,
std::vector<unsigned int>& entitiesA,
std::vector<unsigned int>& entitiesB) const
std::vector<BoundingBoxTreeIntersection<GridView, OtherGridView>>& intersections) const
{
// get alias
const auto& treeA = *this;
......@@ -986,40 +1006,38 @@ private:
auto geometryA = treeA.entity(eIdxA).geometry();
auto geometryB = treeB.entity(eIdxB).geometry();
if (Dumux::GeometryCollision<decltype(geometryA), decltype(geometryB)>::
collide(geometryA, geometryB))
{
// TODO think about data structure here
entitiesA.push_back(eIdxA);
entitiesB.push_back(eIdxB);
}
using CollisionType = Dumux::GeometryCollision<decltype(geometryA), decltype(geometryB)>;
std::vector<GlobalPosition> intersection;
if (CollisionType::collide(geometryA, geometryB, intersection))
intersections.emplace_back(eIdxA, eIdxB, std::move(intersection));
}
// if we reached the leaf in treeA, just continue in treeB
else if (isLeafA)
{
computeCollisions_(treeB, nodeA, bBoxB.child_0, entitiesA, entitiesB);
computeCollisions_(treeB, nodeA, bBoxB.child_1, entitiesA, entitiesB);
computeCollisions_(treeB, nodeA, bBoxB.child_0, intersections);
computeCollisions_(treeB, nodeA, bBoxB.child_1, intersections);
}
// if we reached the leaf in treeB, just continue in treeA
else if (isLeafB)
{
computeCollisions_(treeB, bBoxA.child_0, nodeB, entitiesA, entitiesB);
computeCollisions_(treeB, bBoxA.child_1, nodeB, entitiesA, entitiesB);
computeCollisions_(treeB, bBoxA.child_0, nodeB, intersections);
computeCollisions_(treeB, bBoxA.child_1, nodeB, intersections);
}
// we know now that both trees didn't reach the leaf yet so
// we continue with the larger tree first (bigger node number)
else if (nodeA > nodeB)
{
computeCollisions_(treeB, bBoxA.child_0, nodeB, entitiesA, entitiesB);
computeCollisions_(treeB, bBoxA.child_1, nodeB, entitiesA, entitiesB);
computeCollisions_(treeB, bBoxA.child_0, nodeB, intersections);
computeCollisions_(treeB, bBoxA.child_1, nodeB, intersections);
}
else
{
computeCollisions_(treeB, nodeA, bBoxB.child_0, entitiesA, entitiesB);
computeCollisions_(treeB, nodeA, bBoxB.child_1, entitiesA, entitiesB);
computeCollisions_(treeB, nodeA, bBoxB.child_0, intersections);
computeCollisions_(treeB, nodeA, bBoxB.child_1, intersections);
}
}
......
......@@ -17,12 +17,13 @@
/*!
* \file
* \brief A class for collision detection of two geometries
* and computation of intersections
* and computation of intersection corners
*/
#ifndef DUMUX_GEOMETRY_COLLISION_AND_INTERSECTION_HH
#define DUMUX_GEOMETRY_COLLISION_AND_INTERSECTION_HH
#ifndef DUMUX_GEOMETRY_COLLISION_HH
#define DUMUX_GEOMETRY_COLLISION_HH
#include <dune/common/exceptions.hh>
#include <dune/common/promotiontraits.hh>
#include <dune/geometry/referenceelements.hh>
#include <dumux/common/math.hh>
......@@ -30,6 +31,7 @@
namespace Dumux
{
//! A class for geometry collision detection and intersection calculation
template
<class Geometry1, class Geometry2,
int dimworld = Geometry1::coorddimension,
......@@ -37,11 +39,15 @@ template
int dim2 = Geometry2::mydimension>
class GeometryCollision
{
using Scalar = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
public:
static bool collide(const Geometry1& geo1, const Geometry2& geo2)
//! Determine if the two geometries intersect and compute the intersection corners
static bool collide(const Geometry1& geo1, const Geometry2& geo2, std::vector<GlobalPosition>& intersection)
{
static_assert(dimworld == Geometry2::coorddimension, "Can only collide geometries of same coordinate dimension");
DUNE_THROW(Dune::NotImplemented, "Geometry collision detection not implemented for dimworld = "
DUNE_THROW(Dune::NotImplemented, "Geometry collision detection with intersection computation not implemented for dimworld = "
<< dimworld << ", dim1 = " << dim1 << ", dim2 = " << dim2);
}
};
......@@ -53,8 +59,10 @@ class GeometryCollision<Geometry1, Geometry2, 3, 3, 1>
static const int dimworld = 3;
static const int dim1 = 3;
static const int dim2 = 1;
using Scalar = typename Geometry1::ctype;
static const int dimis = dim2; // the intersection dimension
using Scalar = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim1>;
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
static constexpr Scalar eps_ = 1.5e-7; // base epsilon for floating point comparisons
......@@ -64,8 +72,10 @@ public:
* \note Algorithm from "Real-Time Collision Detection" by Christer Ericson
* Basis is the theorem that for any two non-intersecting convex polyhedrons
* a separating plane exists.
* \param intersection If the geometries collide intersection holds the corner points of
* the intersection object in global coordinates.
*/
static bool collide(const Geometry1& geo1, const Geometry2& geo2)
static bool collide(const Geometry1& geo1, const Geometry2& geo2, std::vector<GlobalPosition>& intersection)
{
static_assert(dimworld == Geometry2::coorddimension, "Can only collide geometries of same coordinate dimension");
......@@ -79,65 +89,71 @@ public:
Scalar tfirst = 0;
Scalar tlast = 1;
std::vector<std::vector<int>> facets;
// sort facet corner so that normal n = (p1-p0)x(p2-p0) always points outwards
switch (geo1.corners())
{
// todo compute intersection geometries
case 8: // hexahedron
facets = {{2, 0, 6, 4}, {1, 3, 5, 7}, {0, 1, 4, 5},
{3, 2, 7, 6}, {1, 0, 3, 2}, {4, 5, 6, 7}};
break;
case 4: // tetrahedron
facets = {{1, 2, 0}, {0, 1, 3}, {0, 3, 2}, {1, 2, 3}};
break;
default:
DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type " << geo1.type() << ", "<< geo1.corners() << " corners.");
}
for (const auto& f : facets)
{
// compute normal vector by cross product
const auto v0 = geo1.corner(f[1]) - geo1.corner(f[0]);
const auto v1 = geo1.corner(f[2]) - geo1.corner(f[0]);
const auto eps = eps_*v0.two_norm();
auto n = Dumux::crossProduct(v0, v1);
n /= n.two_norm();
const Scalar denom = n*d;
const Scalar dist = n*(a-geo1.corner(f[0]));
// if denominator is zero the segment in parallel to
// the plane. If the distance is positive there is no intersection
if (std::abs(denom) < eps)
{
if (dist > eps)
return false;
}
else // not parallel: compute line-plane intersection
{
std::vector<std::vector<int>> facets = {{2, 0, 6, 4}, {1, 3, 5, 7}, {0, 1, 4, 5},
{3, 2, 7, 6}, {1, 0, 3, 2}, {4, 5, 6, 7}};
for (const auto& f : facets)
const Scalar t = -dist / denom;
// if entering half space cut tfirst if t is larger
if (std::signbit(denom))
{
// compute normal vector by cross product
const auto v0 = geo1.corner(f[1]) - geo1.corner(f[0]);
const auto v1 = geo1.corner(f[2]) - geo1.corner(f[0]);
const auto eps = eps_*v0.two_norm();
auto n = Dumux::crossProduct(v0, v1);
n /= n.two_norm();
const Scalar denom = n*d;
const Scalar dist = n*(a-geo1.corner(f[0]));
// if denominator is zero the segment in parallel to
// the plane. If the distance is positive there is no intersection
if (std::abs(denom) < eps)
{
if (dist > eps)
return false;
}
else // not parallel: compute line-plane intersection
{
const Scalar t = -dist / denom;
// if entering half space cut tfirst if t is larger
if (std::signbit(denom))
{
if (t > tfirst)
tfirst = t;
}
// if exiting half space cut tlast if t is smaller
else
{
if (t < tlast)
tlast = t;
}
// there is no intersection of the interval is empty
// use unscaled epsilon since t is in local coordinates
if (tfirst > tlast - eps_)
return false;
}
if (t > tfirst)
tfirst = t;
}
// If we made it until here an intersection exists. We actually
// also have the intersections geometry now s(t) = a + t(b-a) in [tfirst, tlast]
//std::cout << "intersection geometry: " << geo2.global(tfirst) << " -> " << geo2.global(tlast) << std::endl;
return true;
// if exiting half space cut tlast if t is smaller
else
{
if (t < tlast)
tlast = t;
}
// there is no intersection of the interval is empty
// use unscaled epsilon since t is in local coordinates
if (tfirst > tlast - eps_)
return false;
}
default:
DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type " << geo1.type());
}
// If we made it until here an intersection exists. We also export
// the intersections geometry now s(t) = a + t(b-a) in [tfirst, tlast]
intersection = {geo2.global(tfirst), geo2.global(tlast)};
return true;
}
};
} // end namespace Dumux
# endif
\ No newline at end of file
# endif
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