22#ifndef DUMUX_GEOMETRY_DISTANCE_HH
23#define DUMUX_GEOMETRY_DISTANCE_HH
25#include <dune/common/fvector.hh>
26#include <dune/geometry/quadraturerules.hh>
37template<
class Geometry>
38static inline typename Geometry::ctype
40 const Geometry& geometry,
41 std::size_t integrationOrder = 2)
43 typename Geometry::ctype avgDist = 0.0;
44 const auto& quad = Dune::QuadratureRules<typename Geometry::ctype, Geometry::mydimension>::rule(geometry.type(), integrationOrder);
45 for (
const auto& qp : quad)
46 avgDist += (geometry.global(qp.position())-p).two_norm()*qp.weight()*geometry.integrationElement(qp.position());
47 return avgDist/geometry.volume();
55static inline typename Point::value_type
58 const auto ab = b - a;
59 const auto t = (p - a)*ab/ab.two_norm2();
62 return (proj - p).two_norm2();
71template<
class Geometry>
72static inline typename Geometry::ctype
75 static_assert(Geometry::mydimension == 1,
"Geometry has to be a line");
76 const auto& a = geometry.corner(0);
77 const auto& b = geometry.corner(1);
86static inline typename Point::value_type
96template<
class Geometry>
97static inline typename Geometry::ctype
106static inline typename Point::value_type
109 const auto ab = b - a;
110 const auto ap = p - a;
111 const auto t = ap*ab;
114 return ap.two_norm2();
116 const auto lengthSq = ab.two_norm2();
118 return (b - p).two_norm2();
121 proj.axpy(t/lengthSq, ab);
122 return (proj - p).two_norm2();
129template<
class Geometry>
130static inline typename Geometry::ctype
133 static_assert(Geometry::mydimension == 1,
"Geometry has to be a segment");
134 const auto& a = geometry.corner(0);
135 const auto& b = geometry.corner(1);
144static inline typename Point::value_type
152template<
class Geometry>
153static inline typename Geometry::ctype
163static inline typename Point::value_type
166 static_assert(Point::dimension == 3,
"Only works in 3D");
167 const auto ab = b - a;
168 const auto bc = c - b;
169 const auto ca = a - c;
172 const auto ap = p - a;
173 const auto bp = p - b;
174 const auto cp = p - c;
193 const auto tmp =
normal*ap;
202template<
class Geometry>
203static inline typename Geometry::ctype
206 static_assert(Geometry::coorddimension == 3,
"Only works in 3D");
207 static_assert(Geometry::mydimension == 2,
"Geometry has to be a triangle");
208 assert(geometry.corners() == 3);
209 const auto& a = geometry.corner(0);
210 const auto& b = geometry.corner(1);
211 const auto& c = geometry.corner(2);
220static inline typename Point::value_type
228template<
class Geometry>
229static inline typename Geometry::ctype
238template<
class Geometry>
239static inline typename Geometry::ctype
242 if (geometry.corners() == 3)
244 else if (geometry.corners() == 4)
246 const auto& a = geometry.corner(0);
247 const auto& b = geometry.corner(1);
248 const auto& c = geometry.corner(2);
249 const auto& d = geometry.corner(3);
256 DUNE_THROW(Dune::NotImplemented,
"Polygon with " << geometry.corners() <<
" corners not supported");
264template<
class Geometry>
265static inline typename Geometry::ctype
273template<
class Geometry>
274static inline typename Geometry::ctype
276 const typename Geometry::GlobalCoordinate& b,
277 const Geometry& geometry,
278 std::size_t integrationOrder = 2)
280 typename Geometry::ctype avgDist = 0.0;
281 const auto& quad = Dune::QuadratureRules<typename Geometry::ctype, Geometry::mydimension>::rule(geometry.type(), integrationOrder);
282 for (
const auto& qp : quad)
283 avgDist +=
distancePointSegment(geometry.global(qp.position()), a, b)*qp.weight()*geometry.integrationElement(qp.position());
284 return avgDist/geometry.volume();
291template<
class ctype,
int dimWorld>
292static inline ctype
distance(
const Dune::FieldVector<ctype, dimWorld>& a,
293 const Dune::FieldVector<ctype, dimWorld>& b)
294{
return (a-b).two_norm(); }
300template<
class ctype,
int dimWorld>
302 const Dune::FieldVector<ctype, dimWorld>& b)
303{
return (a-b).two_norm2(); }
309template<
class Geo1,
class Geo2,
310 int dimWorld = Geo1::coorddimension,
311 int dim1 = Geo1::mydimension,
int dim2 = Geo2::mydimension>
314 static_assert(Geo1::coorddimension == Geo2::coorddimension,
"Geometries have to have the same coordinate dimensions");
315 static auto distance(
const Geo1& geo1,
const Geo2& geo2)
317 DUNE_THROW(Dune::NotImplemented,
"Geometry distance computation not implemented for dimworld = "
318 << dimWorld <<
", dim1 = " << dim1 <<
", dim2 = " << dim2);
323template<
class Geo1,
class Geo2,
int dimWorld>
326 static_assert(Geo1::coorddimension == Geo2::coorddimension,
"Geometries have to have the same coordinate dimensions");
327 static auto distance(
const Geo1& geo1,
const Geo2& geo2)
332template<
class Geo1,
class Geo2,
int dimWorld>
335 static_assert(Geo1::coorddimension == Geo2::coorddimension,
"Geometries have to have the same coordinate dimensions");
336 static auto distance(
const Geo1& geo1,
const Geo2& geo2)
341template<
class Geo1,
class Geo2,
int dimWorld>
344 static_assert(Geo1::coorddimension == Geo2::coorddimension,
"Geometries have to have the same coordinate dimensions");
345 static auto distance(
const Geo1& geo1,
const Geo2& geo2)
350template<
class Geo1,
class Geo2,
int dimWorld>
353 static_assert(Geo1::coorddimension == Geo2::coorddimension,
"Geometries have to have the same coordinate dimensions");
354 static inline auto distance(
const Geo1& geo1,
const Geo2& geo2)
359template<
class Geo1,
class Geo2,
int dimWorld>
362 static_assert(Geo1::coorddimension == Geo2::coorddimension,
"Geometries have to have the same coordinate dimensions");
363 static inline auto distance(
const Geo1& geo1,
const Geo2& geo2)
373template<
class Geo1,
class Geo2>
381template<
class Geo1,
class Geo2>
382static inline auto distance(
const Geo1& geo1,
const Geo2& geo2)
394template<
class EntitySet,
class ctype,
int dimworld,
395 typename std::enable_if_t<(EntitySet::Entity::Geometry::mydimension > 0),
int> = 0>
399 ctype& minSquaredDistance,
414 else if (tree.
isLeaf(bBox, node))
416 const std::size_t entityIdx = bBox.child1;
418 const auto geometry = tree.
entitySet().entity(entityIdx).geometry();
419 if constexpr (EntitySet::Entity::Geometry::mydimension == 2)
421 else if constexpr (EntitySet::Entity::Geometry::mydimension == 1)
424 DUNE_THROW(Dune::NotImplemented,
"squaredDistance to entity with dim>2");
437 closestEntity(point, tree, bBox.child0, minSquaredDistance, eIdx);
438 closestEntity(point, tree, bBox.child1, minSquaredDistance, eIdx);
447template<
class EntitySet,
class ctype,
int dimworld,
448 typename std::enable_if_t<(EntitySet::Entity::Geometry::mydimension == 0),
int> = 0>
449void closestEntity(
const Dune::FieldVector<ctype, dimworld>& point,
452 ctype& minSquaredDistance,
459 if (tree.
isLeaf(bBox, node))
461 const std::size_t entityIdx = bBox.child1;
462 const auto& p = tree.
entitySet().entity(entityIdx).geometry().corner(0);
483 closestEntity(point, tree, bBox.child0, minSquaredDistance, eIdx);
484 closestEntity(point, tree, bBox.child1, minSquaredDistance, eIdx);
503template<
class EntitySet,
class ctype,
int dimworld>
504std::pair<ctype, std::size_t>
closestEntity(
const Dune::FieldVector<ctype, dimworld>& point,
506 ctype minSquaredDistance = std::numeric_limits<ctype>::max())
508 std::size_t eIdx = 0;
511 return { minSquaredDistance, eIdx };
518template<
class EntitySet,
class ctype,
int dimworld>
521 ctype minSquaredDistance = std::numeric_limits<ctype>::max())
530template<
class EntitySet,
class ctype,
int dimworld>
531ctype
distance(
const Dune::FieldVector<ctype, dimworld>& point,
533 ctype minSquaredDistance = std::numeric_limits<ctype>::max())
Define some often used mathematical functions.
An axis-aligned bounding box volume hierarchy for dune grids.
static ctype distance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest distance between two points.
Definition: distance.hh:292
Vector normal(const Vector &v)
Create a vector normal to the given one (v is expected to be non-zero)
Definition: normal.hh:36
static Geometry::ctype distancePointPolygon(const typename Geometry::GlobalCoordinate &p, const Geometry &geometry)
Compute the shortest distance from a point to a given polygon geometry.
Definition: distance.hh:266
static Point::value_type distancePointSegment(const Point &p, const Point &a, const Point &b)
Compute the distance from a point to the segment connecting the points a and b.
Definition: distance.hh:145
static Point::value_type squaredDistancePointTriangle(const Point &p, const Point &a, const Point &b, const Point &c)
Compute the shortest squared distance from a point to the triangle connecting the points a,...
Definition: distance.hh:164
static Point::value_type squaredDistancePointLine(const Point &p, const Point &a, const Point &b)
Compute the squared distance from a point to a line through the points a and b.
Definition: distance.hh:56
static Point::value_type distancePointTriangle(const Point &p, const Point &a, const Point &b, const Point &c)
Compute the shortest distance from a point to the triangle connecting the points a,...
Definition: distance.hh:221
static Geometry::ctype squaredDistancePointPolygon(const typename Geometry::GlobalCoordinate &p, const Geometry &geometry)
Compute the shortest squared distance from a point to a given polygon geometry.
Definition: distance.hh:240
static Point::value_type squaredDistancePointSegment(const Point &p, const Point &a, const Point &b)
Compute the squared distance from a point to the segment connecting the points a and b.
Definition: distance.hh:107
static ctype squaredDistance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest squared distance between two points.
Definition: distance.hh:301
std::pair< ctype, std::size_t > closestEntity(const Dune::FieldVector< ctype, dimworld > &point, const BoundingBoxTree< EntitySet > &tree, ctype minSquaredDistance=std::numeric_limits< ctype >::max())
Compute the closest entity in an AABB tree to a point (index and shortest squared distance)
Definition: distance.hh:504
static Point::value_type distancePointLine(const Point &p, const Point &a, const Point &b)
Compute the distance from a point to a line through the points a and b.
Definition: distance.hh:87
static Geometry::ctype averageDistanceSegmentGeometry(const typename Geometry::GlobalCoordinate &a, const typename Geometry::GlobalCoordinate &b, const Geometry &geometry, std::size_t integrationOrder=2)
Compute the average distance from a segment to a geometry by integration.
Definition: distance.hh:275
static Geometry::ctype averageDistancePointGeometry(const typename Geometry::GlobalCoordinate &p, const Geometry &geometry, std::size_t integrationOrder=2)
Compute the average distance from a point to a geometry by integration.
Definition: distance.hh:39
void closestEntity(const Dune::FieldVector< ctype, dimworld > &point, const BoundingBoxTree< EntitySet > &tree, std::size_t node, ctype &minSquaredDistance, std::size_t &eIdx)
Compute the closest entity in an AABB tree (index and shortest squared distance) recursively.
Definition: distance.hh:396
Dune::FieldVector< Scalar, 3 > crossProduct(const Dune::FieldVector< Scalar, 3 > &vec1, const Dune::FieldVector< Scalar, 3 > &vec2)
Cross product of two vectors in three-dimensional Euclidean space.
Definition: math.hh:654
constexpr int sign(const ValueType &value) noexcept
Sign or signum function.
Definition: math.hh:641
ctype squaredDistancePointBoundingBox(const Dune::FieldVector< ctype, dimworld > &point, const ctype *b)
Compute squared distance between point and bounding box.
Definition: boundingboxtree.hh:420
An axis-aligned bounding box volume tree implementation.
Definition: boundingboxtree.hh:66
const ctype * getBoundingBoxCoordinates(std::size_t nodeIdx) const
Get an existing bounding box for a given node.
Definition: boundingboxtree.hh:145
bool isLeaf(const BoundingBoxNode &node, std::size_t nodeIdx) const
Definition: boundingboxtree.hh:154
const EntitySet & entitySet() const
the entity set this tree was built with
Definition: boundingboxtree.hh:133
std::size_t numBoundingBoxes() const
Get the number of bounding boxes currently in the tree.
Definition: boundingboxtree.hh:149
const BoundingBoxNode & getBoundingBoxNode(std::size_t nodeIdx) const
Interface to be used by other bounding box trees.
Definition: boundingboxtree.hh:141
Definition: distance.hh:313
static auto distance(const Geo1 &geo1, const Geo2 &geo2)
Definition: distance.hh:315
static auto distance(const Geo1 &geo1, const Geo2 &geo2)
Definition: distance.hh:327
static auto distance(const Geo1 &geo1, const Geo2 &geo2)
Definition: distance.hh:336
static auto distance(const Geo1 &geo1, const Geo2 &geo2)
Definition: distance.hh:345
static auto distance(const Geo1 &geo1, const Geo2 &geo2)
Definition: distance.hh:354
static auto distance(const Geo1 &geo1, const Geo2 &geo2)
Definition: distance.hh:363