Commit 37414d10 authored by Timo Koch's avatar Timo Koch
Browse files

[bboxtree] Add method to intersection two AABB trees

We can collide two AABB trees efficiently now. Thus,
AABB trees of two grids can be intersected. This returns
a pair of vectors. They have the same length. Each entry
contains the element index of the one grid in the first
vector and the intersecting element of the second grid in
the second vector.

This implements primitive tests for intersection of two AABB trees
only for 3d with 1d grids for now. If you want to intersect
grids of other dimensions implement the method collides in the
geometry collision class added with this commit.
parent efc01f0a
......@@ -30,6 +30,7 @@
#include <dune/common/timer.hh>
#include <dune/common/exceptions.hh>
#include <dumux/common/math.hh>
#include <dumux/common/geometrycollision.hh>
namespace Dumux {
......@@ -825,6 +826,28 @@ public:
return entities;
}
// Compute all intersections between entities and a point
template<class OtherGridView>
std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
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;
// Call the recursive find function to find candidates
computeCollisions_(otherTree,
this->numBoundingBoxes_() - 1,
otherTree.numBoundingBoxes_() -1,
myEntities,
otherEntities);
return std::make_pair(myEntities, otherEntities);
}
// Get an element from a global element index
Element entity(unsigned int eIdx) const
{ return indexToElementMap_->entity(eIdx); }
......@@ -930,6 +953,76 @@ private:
}
}
// Compute collisions with other bounding box tree recursively
template <class OtherGridView>
void computeCollisions_(const Dumux::BoundingBoxTree<OtherGridView>& treeB,
unsigned int nodeA,
unsigned int nodeB,
std::vector<unsigned int>& entitiesA,
std::vector<unsigned int>& entitiesB) const
{
// get alias
const auto& treeA = *this;
// Get the bounding box for the current node
const auto& bBoxA = treeA.getBoundingBox_(nodeA);
const auto& bBoxB = treeB.getBoundingBox_(nodeB);
// if the two bounding boxes don't collide we can stop searching
if (!BoundingBoxTreeHelper<dimworld>::
boundingBoxInBoundingBox(treeA.getBoundingBoxCoordinates_(nodeA),
treeB.getBoundingBoxCoordinates_(nodeB)))
return;
// Check if we have a leaf in treeA or treeB
const bool isLeafA = treeA.isLeaf_(bBoxA, nodeA);
const bool isLeafB = treeB.isLeaf_(bBoxB, nodeB);
// If both boxes are leaves and collide add them
if (isLeafA && isLeafB)
{
const unsigned int eIdxA = bBoxA.child_1;
const unsigned int eIdxB = bBoxB.child_1;
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);
}
}
// 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);
}
// 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);
}
// 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);
}
else
{
computeCollisions_(treeB, nodeA, bBoxB.child_0, entitiesA, entitiesB);
computeCollisions_(treeB, nodeA, bBoxB.child_1, entitiesA, entitiesB);
}
}
// Add a new bounding box to the tree
inline unsigned int addBoundingBox_(const BoundingBox& bBox,
const double* b)
......
/*****************************************************************************
* 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 A class for collision detection of two geometries
* and computation of intersections
*/
#ifndef DUMUX_GEOMETRY_COLLISION_AND_INTERSECTION_HH
#define DUMUX_GEOMETRY_COLLISION_AND_INTERSECTION_HH
#include <dune/common/exceptions.hh>
#include <dune/geometry/referenceelements.hh>
#include <dumux/common/math.hh>
namespace Dumux
{
template
<class Geometry1, class Geometry2,
int dimworld = Geometry1::coorddimension,
int dim1 = Geometry1::mydimension,
int dim2 = Geometry2::mydimension>
class GeometryCollision
{
public:
static bool collide(const Geometry1& geo1, const Geometry2& geo2)
{
static_assert(dimworld == Geometry2::coorddimension, "Can only collide geometries of same coordinate dimension");
DUNE_THROW(Dune::NotImplemented, "Geometry collision detection not implemented for dimworld = "
<< dimworld << ", dim1 = " << dim1 << ", dim2 = " << dim2);
}
};
//! Geometry collision detection with 3d and 1d geometry in 3d space
template <class Geometry1, class Geometry2>
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;
using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim1>;
static constexpr Scalar eps_ = 1.5e-7; // base epsilon for floating point comparisons
public:
/*!
* \brief Colliding segment and convex polyhedron
* \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.
*/
static bool collide(const Geometry1& geo1, const Geometry2& geo2)
{
static_assert(dimworld == Geometry2::coorddimension, "Can only collide geometries of same coordinate dimension");
const auto a = geo2.corner(0);
const auto b = geo2.corner(1);
const auto d = b - a;
// The initial interval is the whole segment
// afterward we start clipping the interval
// by the planes decribed by the facet
Scalar tfirst = 0;
Scalar tlast = 1;
switch (geo1.corners())
{
case 8: // hexahedron
{
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)
{
// 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 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;
}
default:
DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type " << geo1.type());
}
}
};
} // end namespace Dumux
# endif
\ No newline at end of file
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