| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup Geometry | ||
| 10 | * \brief A class for collision detection of two geometries | ||
| 11 | * and computation of intersection corners | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_GEOMETRY_INTERSECTION_HH | ||
| 14 | #define DUMUX_GEOMETRY_INTERSECTION_HH | ||
| 15 | |||
| 16 | #include <tuple> | ||
| 17 | |||
| 18 | #include <dune/common/exceptions.hh> | ||
| 19 | #include <dune/common/promotiontraits.hh> | ||
| 20 | #include <dune/geometry/multilineargeometry.hh> | ||
| 21 | #include <dune/geometry/referenceelements.hh> | ||
| 22 | |||
| 23 | #include <dumux/common/math.hh> | ||
| 24 | #include <dumux/geometry/intersectspointgeometry.hh> | ||
| 25 | #include <dumux/geometry/grahamconvexhull.hh> | ||
| 26 | #include <dumux/geometry/boundingboxtree.hh> | ||
| 27 | |||
| 28 | namespace Dumux { | ||
| 29 | namespace IntersectionPolicy { | ||
| 30 | |||
| 31 | //! Policy structure for point-like intersections | ||
| 32 | template<class ct, int dw> | ||
| 33 | struct PointPolicy | ||
| 34 | { | ||
| 35 | using ctype = ct; | ||
| 36 | using Point = Dune::FieldVector<ctype, dw>; | ||
| 37 | |||
| 38 | static constexpr int dimWorld = dw; | ||
| 39 | static constexpr int dimIntersection = 0; | ||
| 40 | using Intersection = Point; | ||
| 41 | }; | ||
| 42 | |||
| 43 | //! Policy structure for segment-like intersections | ||
| 44 | template<class ct, int dw> | ||
| 45 | struct SegmentPolicy | ||
| 46 | { | ||
| 47 | using ctype = ct; | ||
| 48 | using Point = Dune::FieldVector<ctype, dw>; | ||
| 49 | using Segment = std::array<Point, 2>; | ||
| 50 | |||
| 51 | static constexpr int dimWorld = dw; | ||
| 52 | static constexpr int dimIntersection = 1; | ||
| 53 | using Intersection = Segment; | ||
| 54 | }; | ||
| 55 | |||
| 56 | //! Policy structure for polygon-like intersections | ||
| 57 | template<class ct, int dw> | ||
| 58 | struct PolygonPolicy | ||
| 59 | { | ||
| 60 | using ctype = ct; | ||
| 61 | using Point = Dune::FieldVector<ctype, dw>; | ||
| 62 | using Polygon = std::vector<Point>; | ||
| 63 | |||
| 64 | static constexpr int dimWorld = dw; | ||
| 65 | static constexpr int dimIntersection = 2; | ||
| 66 | using Intersection = Polygon; | ||
| 67 | }; | ||
| 68 | |||
| 69 | //! Policy structure for polyhedron-like intersections | ||
| 70 | template<class ct, int dw> | ||
| 71 | struct PolyhedronPolicy | ||
| 72 | { | ||
| 73 | using ctype = ct; | ||
| 74 | using Point = Dune::FieldVector<ctype, dw>; | ||
| 75 | using Polyhedron = std::vector<Point>; | ||
| 76 | |||
| 77 | static constexpr int dimWorld = dw; | ||
| 78 | static constexpr int dimIntersection = 3; | ||
| 79 | using Intersection = Polyhedron; | ||
| 80 | }; | ||
| 81 | |||
| 82 | //! default policy chooser class | ||
| 83 | template<class Geometry1, class Geometry2> | ||
| 84 | class DefaultPolicyChooser | ||
| 85 | { | ||
| 86 | static constexpr int dimworld = Geometry1::coorddimension; | ||
| 87 | static constexpr int isDim = std::min( int(Geometry1::mydimension), int(Geometry2::mydimension) ); | ||
| 88 | static_assert(dimworld == int(Geometry2::coorddimension), | ||
| 89 | "Geometries must have the same coordinate dimension!"); | ||
| 90 | static_assert(int(Geometry1::mydimension) <= 3 && int(Geometry2::mydimension) <= 3, | ||
| 91 | "Geometries must have dimension 3 or less."); | ||
| 92 | using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType; | ||
| 93 | |||
| 94 | using DefaultPolicies = std::tuple<PointPolicy<ctype, dimworld>, | ||
| 95 | SegmentPolicy<ctype, dimworld>, | ||
| 96 | PolygonPolicy<ctype, dimworld>, | ||
| 97 | PolyhedronPolicy<ctype, dimworld>>; | ||
| 98 | public: | ||
| 99 | using type = std::tuple_element_t<isDim, DefaultPolicies>; | ||
| 100 | }; | ||
| 101 | |||
| 102 | //! Helper alias to define the default intersection policy | ||
| 103 | template<class Geometry1, class Geometry2> | ||
| 104 | using DefaultPolicy = typename DefaultPolicyChooser<Geometry1, Geometry2>::type; | ||
| 105 | |||
| 106 | } // end namespace IntersectionPolicy | ||
| 107 | |||
| 108 | namespace Detail { | ||
| 109 | |||
| 110 | /*! | ||
| 111 | * \ingroup Geometry | ||
| 112 | * \brief Algorithm to find segment-like intersections of a polygon/polyhedron with a | ||
| 113 | * segment. The result is stored in the form of the local coordinates tfirst | ||
| 114 | * and tlast on the segment geo1. | ||
| 115 | * \param geo1 the first geometry | ||
| 116 | * \param geo2 the second geometry (dimGeo2 < dimGeo1) | ||
| 117 | * \param baseEps the base epsilon used for floating point comparisons | ||
| 118 | * \param tfirst stores the local coordinate of beginning of intersection segment | ||
| 119 | * \param tlast stores the local coordinate of the end of intersection segment | ||
| 120 | * \param getFacetCornerIndices Function to obtain the facet corner indices on geo1 | ||
| 121 | * \param computeNormal Function to obtain the normal vector on a facet on geo1 | ||
| 122 | */ | ||
| 123 | template< class Geo1, class Geo2, class ctype, | ||
| 124 | class GetFacetCornerIndices, class ComputeNormalFunction > | ||
| 125 | 787323 | bool computeSegmentIntersection(const Geo1& geo1, const Geo2& geo2, ctype baseEps, | |
| 126 | ctype& tfirst, ctype& tlast, | ||
| 127 | const GetFacetCornerIndices& getFacetCornerIndices, | ||
| 128 | const ComputeNormalFunction& computeNormal) | ||
| 129 | { | ||
| 130 | static_assert(int(Geo2::mydimension) == 1, "Geometry2 must be a segment"); | ||
| 131 | static_assert(int(Geo1::mydimension) > int(Geo2::mydimension), | ||
| 132 | "Dimension of Geometry1 must be higher than that of Geometry2"); | ||
| 133 | |||
| 134 | 787323 | const auto a = geo2.corner(0); | |
| 135 | 787323 | const auto b = geo2.corner(1); | |
| 136 | 787323 | const auto d = b - a; | |
| 137 | |||
| 138 | // The initial interval is the whole segment. | ||
| 139 | // Afterwards we start clipping the interval | ||
| 140 | // at the intersections with the facets of geo1 | ||
| 141 | 787323 | tfirst = 0.0; | |
| 142 | 787323 | tlast = 1.0; | |
| 143 | |||
| 144 | 787323 | const auto facets = getFacetCornerIndices(geo1); | |
| 145 |
2/2✓ Branch 0 taken 1651073 times.
✓ Branch 1 taken 133802 times.
|
2434648 | for (const auto& f : facets) |
| 146 | { | ||
| 147 |
1/2✓ Branch 1 taken 25123 times.
✗ Branch 2 not taken.
|
2289051 | const auto n = computeNormal(f); |
| 148 | |||
| 149 |
4/5✓ Branch 1 taken 629111 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 75369 times.
✓ Branch 4 taken 25123 times.
✓ Branch 0 taken 1811964 times.
|
6882741 | const auto c0 = geo1.corner(f[0]); |
| 150 | 2289051 | const ctype denom = n*d; | |
| 151 |
1/2✓ Branch 1 taken 25123 times.
✗ Branch 2 not taken.
|
2314174 | const ctype dist = n*(a-c0); |
| 152 | |||
| 153 | // use first edge of the facet to scale eps | ||
| 154 |
2/4✓ Branch 1 taken 25123 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 25123 times.
✗ Branch 5 not taken.
|
6862368 | const auto edge1 = geo1.corner(f[1]) - geo1.corner(f[0]); |
| 155 |
2/2✓ Branch 0 taken 404804 times.
✓ Branch 1 taken 1246269 times.
|
2289051 | const ctype eps = baseEps*edge1.two_norm(); |
| 156 | |||
| 157 | // if denominator is zero the segment in parallel to the edge. | ||
| 158 | // If the distance is positive there is no intersection | ||
| 159 | using std::abs; | ||
| 160 |
2/2✓ Branch 0 taken 404804 times.
✓ Branch 1 taken 1246269 times.
|
2289051 | if (abs(denom) < eps) |
| 161 | { | ||
| 162 |
2/2✓ Branch 0 taken 335436 times.
✓ Branch 1 taken 69368 times.
|
445017 | if (dist > eps) |
| 163 | 641726 | return false; | |
| 164 | } | ||
| 165 | else // not parallel: compute line-line intersection | ||
| 166 | { | ||
| 167 |
2/2✓ Branch 0 taken 699855 times.
✓ Branch 1 taken 546414 times.
|
1844034 | const ctype t = -dist / denom; |
| 168 | // if entering half space cut tfirst if t is larger | ||
| 169 | using std::signbit; | ||
| 170 |
2/2✓ Branch 0 taken 699855 times.
✓ Branch 1 taken 546414 times.
|
1844034 | if (signbit(denom)) |
| 171 | { | ||
| 172 |
2/2✓ Branch 0 taken 412265 times.
✓ Branch 1 taken 287590 times.
|
1041041 | if (t > tfirst) |
| 173 | 654076 | tfirst = t; | |
| 174 | } | ||
| 175 | // if exiting half space cut tlast if t is smaller | ||
| 176 | else | ||
| 177 | { | ||
| 178 |
2/2✓ Branch 0 taken 384524 times.
✓ Branch 1 taken 161890 times.
|
802993 | if (t < tlast) |
| 179 | 618329 | tlast = t; | |
| 180 | } | ||
| 181 | // there is no intersection of the interval is empty | ||
| 182 | // use unscaled epsilon since t is in local coordinates | ||
| 183 |
2/2✓ Branch 0 taken 856937 times.
✓ Branch 1 taken 389332 times.
|
1844034 | if (tfirst > tlast - baseEps) |
| 184 | return false; | ||
| 185 | } | ||
| 186 | } | ||
| 187 | |||
| 188 | // If we made it until here an intersection exists. | ||
| 189 | return true; | ||
| 190 | 787323 | } | |
| 191 | |||
| 192 | } // end namespace Detail | ||
| 193 | |||
| 194 | /*! | ||
| 195 | * \ingroup Geometry | ||
| 196 | * \brief A class for geometry collision detection and intersection calculation | ||
| 197 | * The class can be specialized for combinations of dimworld, dim1, dim2, where | ||
| 198 | * dimworld is the world dimension embedding a grid of dim1 and a grid of dim2. | ||
| 199 | */ | ||
| 200 | template | ||
| 201 | <class Geometry1, class Geometry2, | ||
| 202 | class Policy = IntersectionPolicy::DefaultPolicy<Geometry1, Geometry2>, | ||
| 203 | int dimworld = Geometry1::coorddimension, | ||
| 204 | int dim1 = Geometry1::mydimension, | ||
| 205 | int dim2 = Geometry2::mydimension> | ||
| 206 | class GeometryIntersection | ||
| 207 | { | ||
| 208 | static constexpr int dimWorld = Policy::dimWorld; | ||
| 209 | |||
| 210 | public: | ||
| 211 | using ctype = typename Policy::ctype; | ||
| 212 | using Point = typename Policy::Point; | ||
| 213 | using Intersection = typename Policy::Intersection; | ||
| 214 | |||
| 215 | //! Determine if the two geometries intersect and compute the intersection geometry | ||
| 216 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | ||
| 217 | { | ||
| 218 | static_assert(dimworld == Geometry2::coorddimension, "Can only intersect geometries of same coordinate dimension"); | ||
| 219 | DUNE_THROW(Dune::NotImplemented, "Geometry intersection detection with intersection computation not implemented for dimworld = " | ||
| 220 | << dimworld << ", dim1 = " << dim1 << ", dim2 = " << dim2); | ||
| 221 | } | ||
| 222 | }; | ||
| 223 | |||
| 224 | /*! | ||
| 225 | * \ingroup Geometry | ||
| 226 | * \brief A class for segment--segment intersection in 2d space | ||
| 227 | */ | ||
| 228 | template <class Geometry1, class Geometry2, class Policy> | ||
| 229 | class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 1, 1> | ||
| 230 | { | ||
| 231 | enum { dimworld = 2 }; | ||
| 232 | enum { dim1 = 1 }; | ||
| 233 | enum { dim2 = 1 }; | ||
| 234 | |||
| 235 | // base epsilon for floating point comparisons | ||
| 236 | static constexpr typename Policy::ctype eps_ = 1.5e-7; | ||
| 237 | |||
| 238 | public: | ||
| 239 | using ctype = typename Policy::ctype; | ||
| 240 | using Point = typename Policy::Point; | ||
| 241 | using Intersection = typename Policy::Intersection; | ||
| 242 | |||
| 243 | /*! | ||
| 244 | * \brief Colliding two segments | ||
| 245 | * \param geo1/geo2 The geometries to intersect | ||
| 246 | * \param intersection The intersection point | ||
| 247 | * \note This overload is used when point-like intersections are seeked | ||
| 248 | */ | ||
| 249 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0> | ||
| 250 | 1898288 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | |
| 251 | { | ||
| 252 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 253 | |||
| 254 | 1898288 | const auto v1 = geo1.corner(1) - geo1.corner(0); | |
| 255 | 1898288 | const auto v2 = geo2.corner(1) - geo2.corner(0); | |
| 256 | 3796576 | const auto ac = geo2.corner(0) - geo1.corner(0); | |
| 257 | |||
| 258 | 1898288 | auto n2 = Point({-1.0*v2[1], v2[0]}); | |
| 259 |
2/2✓ Branch 0 taken 3796576 times.
✓ Branch 1 taken 1898288 times.
|
7593152 | n2 /= n2.two_norm(); |
| 260 | |||
| 261 | // compute distance of first corner on geo2 to line1 | ||
| 262 |
2/2✓ Branch 0 taken 3796576 times.
✓ Branch 1 taken 1898288 times.
|
5694864 | const auto dist12 = n2*ac; |
| 263 | |||
| 264 | // first check parallel segments | ||
| 265 | using std::abs; | ||
| 266 | using std::sqrt; | ||
| 267 | |||
| 268 | 1898288 | const auto v1Norm2 = v1.two_norm2(); | |
| 269 | 1898288 | const auto eps = eps_*sqrt(v1Norm2); | |
| 270 | 1898288 | const auto eps2 = eps_*v1Norm2; | |
| 271 | |||
| 272 |
2/2✓ Branch 0 taken 948960 times.
✓ Branch 1 taken 949328 times.
|
1898288 | const auto sp = n2*v1; |
| 273 |
2/2✓ Branch 0 taken 948960 times.
✓ Branch 1 taken 949328 times.
|
1898288 | if (abs(sp) < eps) |
| 274 | { | ||
| 275 |
2/2✓ Branch 0 taken 735232 times.
✓ Branch 1 taken 213728 times.
|
948960 | if (abs(dist12) > eps) |
| 276 | return false; | ||
| 277 | |||
| 278 | // point intersection can only be one of the corners | ||
| 279 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 213720 times.
|
213728 | if ( (v1*v2) < 0.0 ) |
| 280 | { | ||
| 281 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
|
8 | if ( ac.two_norm2() < eps2 ) |
| 282 | 4 | { intersection = geo2.corner(0); return true; } | |
| 283 | |||
| 284 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
|
12 | if ( (geo2.corner(1) - geo1.corner(1)).two_norm2() < eps2 ) |
| 285 | ✗ | { intersection = geo2.corner(1); return true; } | |
| 286 | } | ||
| 287 | else | ||
| 288 | { | ||
| 289 |
2/2✓ Branch 2 taken 63896 times.
✓ Branch 3 taken 149824 times.
|
641160 | if ( (geo2.corner(1) - geo1.corner(0)).two_norm2() < eps2 ) |
| 290 | 63896 | { intersection = geo2.corner(1); return true; } | |
| 291 | |||
| 292 |
2/2✓ Branch 2 taken 63900 times.
✓ Branch 3 taken 85924 times.
|
449472 | if ( (geo2.corner(0) - geo1.corner(1)).two_norm2() < eps2 ) |
| 293 | 63900 | { intersection = geo2.corner(0); return true; } | |
| 294 | } | ||
| 295 | |||
| 296 | // no intersection | ||
| 297 | 85928 | return false; | |
| 298 | } | ||
| 299 | |||
| 300 | // intersection point on v1 in local coords | ||
| 301 | 949328 | const auto t1 = dist12 / sp; | |
| 302 | |||
| 303 | // check if the local coords are valid | ||
| 304 |
4/4✓ Branch 0 taken 152092 times.
✓ Branch 1 taken 797236 times.
✓ Branch 2 taken 152144 times.
✓ Branch 3 taken 645092 times.
|
949328 | if (t1 < -1.0*eps_ || t1 > 1.0 + eps_) |
| 305 | return false; | ||
| 306 | |||
| 307 | // compute global coordinates | ||
| 308 | 645092 | auto isPoint = geo1.global(t1); | |
| 309 | |||
| 310 | // check if point is in bounding box of v2 | ||
| 311 | 645092 | const auto c = geo2.corner(0); | |
| 312 |
2/2✓ Branch 1 taken 60 times.
✓ Branch 2 taken 645032 times.
|
645092 | const auto d = geo2.corner(1); |
| 313 | |||
| 314 | using std::min; using std::max; | ||
| 315 |
10/10✓ Branch 0 taken 60 times.
✓ Branch 1 taken 645032 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 645084 times.
✓ Branch 4 taken 322552 times.
✓ Branch 5 taken 322540 times.
✓ Branch 6 taken 322580 times.
✓ Branch 7 taken 322512 times.
✓ Branch 8 taken 342544 times.
✓ Branch 9 taken 302548 times.
|
1290292 | std::array<ctype, 4> bBox({ min(c[0], d[0]), min(c[1], d[1]), max(c[0], d[0]), max(c[1], d[1]) }); |
| 316 | |||
| 317 |
2/2✓ Branch 0 taken 342544 times.
✓ Branch 1 taken 302548 times.
|
645092 | if ( intersectsPointBoundingBox(isPoint, bBox.data()) ) |
| 318 | { | ||
| 319 | 302548 | intersection = std::move(isPoint); | |
| 320 | 302548 | return true; | |
| 321 | } | ||
| 322 | |||
| 323 | return false; | ||
| 324 | } | ||
| 325 | |||
| 326 | /*! | ||
| 327 | * \brief Colliding two segments | ||
| 328 | * \param geo1/geo2 The geometries to intersect | ||
| 329 | * \param intersection Container to store corners of intersection segment | ||
| 330 | * \note this overload is used when segment-like intersections are seeked | ||
| 331 | */ | ||
| 332 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0> | ||
| 333 | 70 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | |
| 334 | { | ||
| 335 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 336 | |||
| 337 | 70 | const auto& a = geo1.corner(0); | |
| 338 | 70 | const auto& b = geo1.corner(1); | |
| 339 | 70 | const auto ab = b - a; | |
| 340 | |||
| 341 | 70 | const auto& p = geo2.corner(0); | |
| 342 | 70 | const auto& q = geo2.corner(1); | |
| 343 | 70 | const auto pq = q - p; | |
| 344 | 70 | Dune::FieldVector<ctype, dimworld> n2{-pq[1], pq[0]}; | |
| 345 | |||
| 346 | using std::max; | ||
| 347 | 70 | const auto abNorm2 = ab.two_norm2(); | |
| 348 |
2/2✓ Branch 0 taken 20 times.
✓ Branch 1 taken 50 times.
|
70 | const auto pqNorm2 = pq.two_norm2(); |
| 349 | 70 | const auto eps2 = eps_*max(abNorm2, pqNorm2); | |
| 350 | |||
| 351 | // non-parallel segments do not intersect in a segment. | ||
| 352 | using std::abs; | ||
| 353 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 66 times.
|
70 | if (abs(n2*ab) > eps2) |
| 354 | return false; | ||
| 355 | |||
| 356 | // check if the segments lie on the same line | ||
| 357 | 66 | const auto ap = p - a; | |
| 358 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 62 times.
|
66 | if (abs(ap*n2) > eps2) |
| 359 | return false; | ||
| 360 | |||
| 361 | // compute scaled local coordinates of corner 1/2 of segment2 on segment1 | ||
| 362 | 62 | auto t1 = ab*ap; | |
| 363 | 62 | auto t2 = ab*(q - a); | |
| 364 | |||
| 365 | using std::swap; | ||
| 366 |
2/2✓ Branch 0 taken 33 times.
✓ Branch 1 taken 29 times.
|
62 | if (t1 > t2) |
| 367 | 33 | swap(t1, t2); | |
| 368 | |||
| 369 | using std::clamp; | ||
| 370 |
2/2✓ Branch 0 taken 16 times.
✓ Branch 1 taken 46 times.
|
62 | t1 = clamp(t1, 0.0, abNorm2); |
| 371 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 58 times.
|
62 | t2 = clamp(t2, 0.0, abNorm2); |
| 372 | |||
| 373 |
2/2✓ Branch 0 taken 32 times.
✓ Branch 1 taken 30 times.
|
62 | if (abs(t2-t1) < eps2) |
| 374 | return false; | ||
| 375 | |||
| 376 | 32 | intersection = Intersection({geo1.global(t1/abNorm2), geo1.global(t2/abNorm2)}); | |
| 377 | 30 | return true; | |
| 378 | } | ||
| 379 | }; | ||
| 380 | |||
| 381 | /*! | ||
| 382 | * \ingroup Geometry | ||
| 383 | * \brief A class for polygon--segment intersection in 2d space | ||
| 384 | */ | ||
| 385 | template <class Geometry1, class Geometry2, class Policy> | ||
| 386 | class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 2, 1> | ||
| 387 | { | ||
| 388 | enum { dimworld = 2 }; | ||
| 389 | enum { dim1 = 2 }; | ||
| 390 | enum { dim2 = 1 }; | ||
| 391 | |||
| 392 | public: | ||
| 393 | using ctype = typename Policy::ctype; | ||
| 394 | using Point = typename Policy::Point; | ||
| 395 | using Intersection = typename Policy::Intersection; | ||
| 396 | |||
| 397 | private: | ||
| 398 | static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons | ||
| 399 | |||
| 400 | public: | ||
| 401 | /*! | ||
| 402 | * \brief Colliding segment and convex polygon | ||
| 403 | * \param geo1/geo2 The geometries to intersect | ||
| 404 | * \param intersection If the geometries collide intersection holds the | ||
| 405 | * corner points of the intersection object in global coordinates. | ||
| 406 | * \note This overload is used when segment-like intersections are seeked. | ||
| 407 | */ | ||
| 408 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0> | ||
| 409 | 4750 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | |
| 410 | { | ||
| 411 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 412 | |||
| 413 | ctype tfirst, tlast; | ||
| 414 |
2/2✓ Branch 1 taken 500 times.
✓ Branch 2 taken 4250 times.
|
4750 | if (intersect_(geo1, geo2, tfirst, tlast)) |
| 415 | { | ||
| 416 | // the intersection exists. Export the intersections geometry now: | ||
| 417 | // s(t) = a + t(b-a) in [tfirst, tlast] | ||
| 418 | 500 | intersection = Intersection({geo2.global(tfirst), geo2.global(tlast)}); | |
| 419 | 500 | return true; | |
| 420 | } | ||
| 421 | |||
| 422 | return false; | ||
| 423 | } | ||
| 424 | |||
| 425 | /*! | ||
| 426 | * \brief Colliding segment and convex polygon | ||
| 427 | * \param geo1/geo2 The geometries to intersect | ||
| 428 | * \param intersection If the geometries collide intersection holds the | ||
| 429 | * corner points of the intersection object in global coordinates. | ||
| 430 | * \note This overload is used when point-like intersections are seeked. | ||
| 431 | */ | ||
| 432 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0> | ||
| 433 | 24 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename P::Intersection& intersection) | |
| 434 | { | ||
| 435 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 436 | |||
| 437 | ctype tfirst, tlast; | ||
| 438 |
2/2✓ Branch 1 taken 4 times.
✓ Branch 2 taken 20 times.
|
24 | if (!intersect_(geo1, geo2, tfirst, tlast)) |
| 439 | { | ||
| 440 | // if start & end point are same, the intersection is a point | ||
| 441 | using std::abs; | ||
| 442 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if ( tfirst > tlast - eps_ ) |
| 443 | { | ||
| 444 | 4 | intersection = Intersection(geo2.global(tfirst)); | |
| 445 | 4 | return true; | |
| 446 | } | ||
| 447 | } | ||
| 448 | |||
| 449 | return false; | ||
| 450 | } | ||
| 451 | |||
| 452 | private: | ||
| 453 | /*! | ||
| 454 | * \brief Obtain local coordinates of start/end point of the intersecting segment. | ||
| 455 | */ | ||
| 456 | 4822 | static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast) | |
| 457 | { | ||
| 458 | // lambda to obtain the facet corners on geo1 | ||
| 459 | 4774 | auto getFacetCorners = [] (const Geometry1& geo1) | |
| 460 | { | ||
| 461 |
1/2✓ Branch 1 taken 4726 times.
✗ Branch 2 not taken.
|
4774 | std::vector< std::array<int, 2> > facetCorners; |
| 462 |
4/6✓ Branch 0 taken 13 times.
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13 times.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
|
48 | switch (geo1.corners()) |
| 463 | { | ||
| 464 | 26 | case 4: // quadrilateral | |
| 465 |
2/5✓ Branch 1 taken 4739 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
|
4752 | facetCorners = {{0, 2}, {3, 1}, {1, 0}, {2, 3}}; |
| 466 | 26 | break; | |
| 467 | 22 | case 3: // triangle | |
| 468 |
2/8✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 6 taken 11 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
22 | facetCorners = {{1, 0}, {0, 2}, {2, 1}}; |
| 469 | 22 | break; | |
| 470 | ✗ | default: | |
| 471 |
0/60✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 83 not taken.
✗ Branch 84 not taken.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 93 not taken.
✗ Branch 94 not taken.
✗ Branch 96 not taken.
✗ Branch 97 not taken.
|
4726 | DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type " |
| 472 | << geo1.type() << " with "<< geo1.corners() << " corners."); | ||
| 473 | } | ||
| 474 | |||
| 475 | 4774 | return facetCorners; | |
| 476 | ✗ | }; | |
| 477 | |||
| 478 | // lambda to obtain the normal on a facet | ||
| 479 | 4822 | const auto center1 = geo1.center(); | |
| 480 | 21727 | auto computeNormal = [&geo1, ¢er1] (const std::array<int, 2>& facetCorners) | |
| 481 | { | ||
| 482 | 16905 | const auto c0 = geo1.corner(facetCorners[0]); | |
| 483 | 33642 | const auto c1 = geo1.corner(facetCorners[1]); | |
| 484 | 16905 | const auto edge = c1 - c0; | |
| 485 | |||
| 486 | 16905 | Dune::FieldVector<ctype, dimworld> n({-1.0*edge[1], edge[0]}); | |
| 487 | 16905 | n /= n.two_norm(); | |
| 488 | |||
| 489 | // orientation of the element is not known | ||
| 490 | // make sure normal points outwards of element | ||
| 491 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 16821 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 84 times.
|
50715 | if ( n*(center1-c0) > 0.0 ) |
| 492 | 16905 | n *= -1.0; | |
| 493 | |||
| 494 | 16905 | return n; | |
| 495 | }; | ||
| 496 | |||
| 497 | 4822 | return Detail::computeSegmentIntersection(geo1, geo2, eps_, tfirst, tlast, getFacetCorners, computeNormal); | |
| 498 | } | ||
| 499 | }; | ||
| 500 | |||
| 501 | /*! | ||
| 502 | * \ingroup Geometry | ||
| 503 | * \brief A class for segment--polygon intersection in 2d space | ||
| 504 | */ | ||
| 505 | template <class Geometry1, class Geometry2, class Policy> | ||
| 506 | class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 1, 2> | ||
| 507 | : public GeometryIntersection<Geometry2, Geometry1, Policy, 2, 2, 1> | ||
| 508 | { | ||
| 509 | using Base = GeometryIntersection<Geometry2, Geometry1, Policy, 2, 2, 1>; | ||
| 510 | |||
| 511 | public: | ||
| 512 | /*! | ||
| 513 | * \brief Colliding segment and convex polygon | ||
| 514 | * \param geo1/geo2 The geometries to intersect | ||
| 515 | * \param intersection If the geometries collide intersection holds the | ||
| 516 | * corner points of the intersection object in global coordinates. | ||
| 517 | * \note This forwards to the polygon-segment specialization with swapped arguments. | ||
| 518 | */ | ||
| 519 | template<class P = Policy> | ||
| 520 | 4428 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename Base::Intersection& intersection) | |
| 521 | { | ||
| 522 |
1/2✓ Branch 1 taken 4384 times.
✗ Branch 2 not taken.
|
4428 | return Base::intersection(geo2, geo1, intersection); |
| 523 | } | ||
| 524 | }; | ||
| 525 | |||
| 526 | /*! | ||
| 527 | * \ingroup Geometry | ||
| 528 | * \brief A class for polygon--polygon intersection in 2d space | ||
| 529 | */ | ||
| 530 | template <class Geometry1, class Geometry2, class Policy> | ||
| 531 | class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 2, 2> | ||
| 532 | { | ||
| 533 | enum { dimworld = 2 }; | ||
| 534 | enum { dim1 = 2 }; | ||
| 535 | enum { dim2 = 2 }; | ||
| 536 | |||
| 537 | public: | ||
| 538 | using ctype = typename Policy::ctype; | ||
| 539 | using Point = typename Policy::Point; | ||
| 540 | using Intersection = typename Policy::Intersection; | ||
| 541 | |||
| 542 | private: | ||
| 543 | static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons | ||
| 544 | |||
| 545 | public: | ||
| 546 | /*! | ||
| 547 | * \brief Colliding two polygons | ||
| 548 | * \note First we find the vertex candidates for the intersection region as follows: | ||
| 549 | * Add polygon vertices that are inside the other polygon | ||
| 550 | * Add intersections of polygon edges | ||
| 551 | * Remove duplicate points from the list | ||
| 552 | * Compute the convex hull polygon | ||
| 553 | * Return a triangulation of that polygon as intersection | ||
| 554 | * \param geo1/geo2 The geometries to intersect | ||
| 555 | * \param intersection Container to store the corner points of the polygon (as convex hull) | ||
| 556 | * \note This overload is used when polygon like intersections are seeked | ||
| 557 | */ | ||
| 558 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 2, int> = 0> | ||
| 559 | 196448 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | |
| 560 | { | ||
| 561 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 562 | |||
| 563 | // the candidate intersection points | ||
| 564 |
1/2✓ Branch 1 taken 196188 times.
✗ Branch 2 not taken.
|
196448 | std::vector<Point> points; points.reserve(6); |
| 565 | |||
| 566 | // add polygon1 corners that are inside polygon2 | ||
| 567 |
3/3✓ Branch 0 taken 784336 times.
✓ Branch 1 taken 196396 times.
✓ Branch 2 taken 104 times.
|
982072 | for (int i = 0; i < geo1.corners(); ++i) |
| 568 |
5/6✓ Branch 2 taken 300 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 44 times.
✓ Branch 5 taken 76 times.
✓ Branch 0 taken 134472 times.
✓ Branch 1 taken 649876 times.
|
785624 | if (intersectsPointGeometry(geo1.corner(i), geo2)) |
| 569 |
2/3✓ Branch 1 taken 134472 times.
✓ Branch 2 taken 56 times.
✗ Branch 3 not taken.
|
134556 | points.emplace_back(geo1.corner(i)); |
| 570 | |||
| 571 |
2/2✓ Branch 0 taken 195412 times.
✓ Branch 1 taken 580 times.
|
196448 | const auto numPoints1 = points.size(); |
| 572 |
2/2✓ Branch 0 taken 195608 times.
✓ Branch 1 taken 580 times.
|
196448 | if (numPoints1 != geo1.corners()) |
| 573 | { | ||
| 574 | // add polygon2 corners that are inside polygon1 | ||
| 575 |
3/3✓ Branch 0 taken 782288 times.
✓ Branch 1 taken 195692 times.
✓ Branch 2 taken 36 times.
|
979316 | for (int i = 0; i < geo2.corners(); ++i) |
| 576 |
6/6✓ Branch 2 taken 120 times.
✓ Branch 3 taken 124 times.
✓ Branch 4 taken 180 times.
✓ Branch 5 taken 72 times.
✓ Branch 0 taken 511388 times.
✓ Branch 1 taken 270900 times.
|
783448 | if (intersectsPointGeometry(geo2.corner(i), geo1)) |
| 577 |
2/3✓ Branch 1 taken 511512 times.
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
|
512260 | points.emplace_back(geo2.corner(i)); |
| 578 | |||
| 579 |
2/2✓ Branch 0 taken 195596 times.
✓ Branch 1 taken 12 times.
|
195868 | if (points.empty()) |
| 580 | return false; | ||
| 581 | |||
| 582 |
2/2✓ Branch 0 taken 118664 times.
✓ Branch 1 taken 76932 times.
|
195844 | if (points.size() - numPoints1 != geo2.corners()) |
| 583 | { | ||
| 584 |
1/2✓ Branch 1 taken 118628 times.
✗ Branch 2 not taken.
|
118800 | const auto refElement1 = referenceElement(geo1); |
| 585 | 118800 | const auto refElement2 = referenceElement(geo2); | |
| 586 | |||
| 587 | // add intersections of edges | ||
| 588 | using SegGeometry = Dune::MultiLinearGeometry<ctype, 1, dimworld>; | ||
| 589 | using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>; | ||
| 590 |
2/2✓ Branch 1 taken 474580 times.
✓ Branch 2 taken 118664 times.
|
593884 | for (int i = 0; i < refElement1.size(dim1-1); ++i) |
| 591 | { | ||
| 592 | 475084 | const auto localEdgeGeom1 = refElement1.template geometry<dim1-1>(i); | |
| 593 |
1/2✓ Branch 0 taken 474580 times.
✗ Branch 1 not taken.
|
475084 | const auto edge1 = SegGeometry( Dune::GeometryTypes::line, |
| 594 |
1/2✓ Branch 2 taken 228 times.
✗ Branch 3 not taken.
|
475432 | std::vector<Point>( {geo1.global(localEdgeGeom1.corner(0)), |
| 595 |
1/2✓ Branch 1 taken 474352 times.
✗ Branch 2 not taken.
|
949820 | geo1.global(localEdgeGeom1.corner(1))} )); |
| 596 | |||
| 597 |
2/2✓ Branch 1 taken 1898248 times.
✓ Branch 2 taken 474580 times.
|
2375348 | for (int j = 0; j < refElement2.size(dim2-1); ++j) |
| 598 | { | ||
| 599 | 1900264 | const auto localEdgeGeom2 = refElement2.template geometry<dim2-1>(j); | |
| 600 |
1/2✓ Branch 0 taken 1898248 times.
✗ Branch 1 not taken.
|
1900264 | const auto edge2 = SegGeometry( Dune::GeometryTypes::line, |
| 601 |
1/6✓ Branch 2 taken 360 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
1900624 | std::vector<Point>( {geo2.global(localEdgeGeom2.corner(0)), |
| 602 |
1/2✓ Branch 1 taken 1897888 times.
✗ Branch 2 not taken.
|
3800168 | geo2.global(localEdgeGeom2.corner(1))} )); |
| 603 | |||
| 604 | using EdgeTest = GeometryIntersection<SegGeometry, SegGeometry, PointPolicy>; | ||
| 605 | 1900264 | typename EdgeTest::Intersection edgeIntersection; | |
| 606 |
2/2✓ Branch 1 taken 430324 times.
✓ Branch 2 taken 1467924 times.
|
1900264 | if (EdgeTest::intersection(edge1, edge2, edgeIntersection)) |
| 607 |
1/2✓ Branch 1 taken 430324 times.
✗ Branch 2 not taken.
|
430644 | points.emplace_back(edgeIntersection); |
| 608 | } | ||
| 609 | } | ||
| 610 | } | ||
| 611 | } | ||
| 612 | |||
| 613 |
1/2✓ Branch 0 taken 196176 times.
✗ Branch 1 not taken.
|
196424 | if (points.empty()) |
| 614 | return false; | ||
| 615 | |||
| 616 | // remove duplicates | ||
| 617 | 589128 | const auto eps = (geo1.corner(0) - geo1.corner(1)).two_norm()*eps_; | |
| 618 |
12/48✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 315993 times.
✓ Branch 5 taken 855602 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 368 times.
✓ Branch 9 taken 364 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 308676 times.
✓ Branch 23 taken 570760 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 136 times.
✓ Branch 27 taken 192 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✓ Branch 44 taken 372 times.
✓ Branch 45 taken 216 times.
✓ Branch 46 taken 112 times.
✓ Branch 47 taken 100 times.
|
2249315 | std::sort(points.begin(), points.end(), [&eps](const auto& a, const auto& b) -> bool |
| 619 | { | ||
| 620 | using std::abs; | ||
| 621 |
12/48✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 315993 times.
✓ Branch 5 taken 855602 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 368 times.
✓ Branch 9 taken 364 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 308676 times.
✓ Branch 23 taken 570760 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 136 times.
✓ Branch 27 taken 192 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✓ Branch 44 taken 372 times.
✓ Branch 45 taken 216 times.
✓ Branch 46 taken 112 times.
✓ Branch 47 taken 100 times.
|
2052891 | return (abs(a[0]-b[0]) > eps ? a[0] < b[0] : a[1] < b[1]); |
| 622 | }); | ||
| 623 | |||
| 624 | 1076660 | auto removeIt = std::unique(points.begin(), points.end(), [&eps](const auto& a, const auto&b) | |
| 625 | { | ||
| 626 | 880236 | return (b-a).two_norm() < eps; | |
| 627 | }); | ||
| 628 | |||
| 629 |
2/2✓ Branch 1 taken 78952 times.
✓ Branch 2 taken 117224 times.
|
196424 | points.erase(removeIt, points.end()); |
| 630 | |||
| 631 | // return false if we don't have at least three unique points | ||
| 632 |
2/2✓ Branch 0 taken 78952 times.
✓ Branch 1 taken 117224 times.
|
196424 | if (points.size() < 3) |
| 633 | return false; | ||
| 634 | |||
| 635 | // intersection polygon is convex hull of above points | ||
| 636 |
2/6✗ Branch 1 not taken.
✓ Branch 2 taken 78952 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 78952 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
79188 | intersection = grahamConvexHull<2>(points); |
| 637 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78952 times.
|
79188 | assert(!intersection.empty()); |
| 638 | return true; | ||
| 639 | 392480 | } | |
| 640 | |||
| 641 | /*! | ||
| 642 | * \brief Colliding two polygons | ||
| 643 | * \param geo1/geo2 The geometries to intersect | ||
| 644 | * \param intersection Container to store the corners of intersection segment | ||
| 645 | * \note this overload is used when segment-like intersections are seeked | ||
| 646 | * \todo implement this query | ||
| 647 | */ | ||
| 648 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0> | ||
| 649 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | ||
| 650 | { | ||
| 651 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 652 | DUNE_THROW(Dune::NotImplemented, "Polygon-polygon intersection detection for segment-like intersections"); | ||
| 653 | } | ||
| 654 | |||
| 655 | /*! | ||
| 656 | * \brief Colliding two polygons | ||
| 657 | * \param geo1/geo2 The geometries to intersect | ||
| 658 | * \param intersection The intersection point | ||
| 659 | * \note this overload is used when point-like intersections are seeked | ||
| 660 | * \todo implement this query | ||
| 661 | */ | ||
| 662 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0> | ||
| 663 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | ||
| 664 | { | ||
| 665 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 666 | DUNE_THROW(Dune::NotImplemented, "Polygon-polygon intersection detection for touching points"); | ||
| 667 | } | ||
| 668 | }; | ||
| 669 | |||
| 670 | /*! | ||
| 671 | * \ingroup Geometry | ||
| 672 | * \brief A class for polyhedron--segment intersection in 3d space | ||
| 673 | */ | ||
| 674 | template <class Geometry1, class Geometry2, class Policy> | ||
| 675 | class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 3, 1> | ||
| 676 | { | ||
| 677 | enum { dimworld = 3 }; | ||
| 678 | enum { dim1 = 3 }; | ||
| 679 | enum { dim2 = 1 }; | ||
| 680 | |||
| 681 | public: | ||
| 682 | using ctype = typename Policy::ctype; | ||
| 683 | using Point = typename Policy::Point; | ||
| 684 | using Intersection = typename Policy::Intersection; | ||
| 685 | |||
| 686 | private: | ||
| 687 | static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons | ||
| 688 | |||
| 689 | public: | ||
| 690 | /*! | ||
| 691 | * \brief Colliding segment and convex polyhedron | ||
| 692 | * \note Algorithm based on the one from "Real-Time Collision Detection" by Christer Ericson, | ||
| 693 | * published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc. | ||
| 694 | * Basis is the theorem that for any two non-intersecting convex polyhedrons | ||
| 695 | * a separating plane exists. | ||
| 696 | * \param geo1/geo2 The geometries to intersect | ||
| 697 | * \param intersection If the geometries collide intersection holds the corner points of | ||
| 698 | * the intersection object in global coordinates. | ||
| 699 | * \note This overload is used when segment-like intersections are seeked. | ||
| 700 | */ | ||
| 701 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0> | ||
| 702 | 346620 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | |
| 703 | { | ||
| 704 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 705 | |||
| 706 | ctype tfirst, tlast; | ||
| 707 |
2/2✓ Branch 0 taken 49864 times.
✓ Branch 1 taken 150548 times.
|
346620 | if (intersect_(geo1, geo2, tfirst, tlast)) |
| 708 | { | ||
| 709 | // Intersection exists. Export the intersections geometry now: | ||
| 710 | // s(t) = a + t(b-a) in [tfirst, tlast] | ||
| 711 | 66416 | intersection = Intersection({geo2.global(tfirst), geo2.global(tlast)}); | |
| 712 | 50408 | return true; | |
| 713 | } | ||
| 714 | |||
| 715 | return false; | ||
| 716 | } | ||
| 717 | |||
| 718 | /*! | ||
| 719 | * \brief Colliding segment and convex polyhedron | ||
| 720 | * \note Algorithm based on the one from "Real-Time Collision Detection" by Christer Ericson, | ||
| 721 | * published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc. | ||
| 722 | * Basis is the theorem that for any two non-intersecting convex polyhedrons | ||
| 723 | * a separating plane exists. | ||
| 724 | * \param geo1/geo2 The geometries to intersect | ||
| 725 | * \param intersection If the geometries collide intersection holds the corner points of | ||
| 726 | * the intersection object in global coordinates. | ||
| 727 | * \note This overload is used when point-like intersections are seeked. | ||
| 728 | */ | ||
| 729 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0> | ||
| 730 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | ||
| 731 | { | ||
| 732 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 733 | |||
| 734 | ctype tfirst, tlast; | ||
| 735 | if (intersect_(geo1, geo2, tfirst, tlast)) | ||
| 736 | { | ||
| 737 | // if start & end point are same, the intersection is a point | ||
| 738 | using std::abs; | ||
| 739 | if ( abs(tfirst - tlast) < eps_ ) | ||
| 740 | { | ||
| 741 | intersection = Intersection(geo2.global(tfirst)); | ||
| 742 | return true; | ||
| 743 | } | ||
| 744 | } | ||
| 745 | |||
| 746 | return false; | ||
| 747 | } | ||
| 748 | |||
| 749 | private: | ||
| 750 | /*! | ||
| 751 | * \brief Obtain local coordinates of start/end point of the intersecting segment. | ||
| 752 | */ | ||
| 753 | 200412 | static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast) | |
| 754 | { | ||
| 755 | // lambda to obtain facet corners on geo1 | ||
| 756 |
16/49✓ Branch 1 taken 180252 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 180252 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 180252 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 180252 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 180252 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 9040 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 9040 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 9040 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 9040 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 6400 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 6400 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 6400 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 6400 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 4720 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 4720 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 4720 times.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 15 not taken.
|
693240 | auto getFacetCorners = [] (const Geometry1& geo1) |
| 757 | { | ||
| 758 |
3/6✓ Branch 0 taken 4195 times.
✓ Branch 1 taken 157865 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2432 times.
✗ Branch 5 not taken.
|
200412 | std::vector< std::vector<int> > facetCorners; |
| 759 | // sort facet corner so that normal n = (p1-p0)x(p2-p0) always points outwards | ||
| 760 |
4/5✓ Branch 0 taken 19955 times.
✓ Branch 1 taken 9040 times.
✓ Branch 2 taken 6400 times.
✓ Branch 3 taken 4720 times.
✗ Branch 4 not taken.
|
40115 | switch (geo1.corners()) |
| 761 | { | ||
| 762 | // todo compute intersection geometries | ||
| 763 |
1/2✓ Branch 1 taken 19955 times.
✗ Branch 2 not taken.
|
19955 | case 8: // hexahedron |
| 764 |
9/19✓ Branch 1 taken 177820 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 177820 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1066920 times.
✓ Branch 7 taken 177820 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 13 taken 2432 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2432 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 14592 times.
✓ Branch 19 taken 2432 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 8 taken 94560 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
1622268 | facetCorners = {{2, 0, 6, 4}, {1, 3, 5, 7}, {0, 1, 4, 5}, |
| 765 | {3, 2, 7, 6}, {1, 0, 3, 2}, {4, 5, 6, 7}}; | ||
| 766 | 19955 | break; | |
| 767 |
1/2✓ Branch 1 taken 9040 times.
✗ Branch 2 not taken.
|
9040 | case 6: // prism |
| 768 |
5/12✓ Branch 1 taken 9040 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9040 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 45200 times.
✓ Branch 7 taken 9040 times.
✓ Branch 8 taken 45200 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
72320 | facetCorners = {{0, 2, 1}, {3, 4, 5}, {0, 1, 3, 4}, {0, 3, 2, 5}, {1, 2, 4, 5}}; |
| 769 | 9040 | break; | |
| 770 |
1/2✓ Branch 1 taken 6400 times.
✗ Branch 2 not taken.
|
6400 | case 5: // pyramid |
| 771 |
5/12✓ Branch 1 taken 6400 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6400 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 32000 times.
✓ Branch 7 taken 6400 times.
✓ Branch 8 taken 32000 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
51200 | facetCorners = {{0, 2, 1, 3}, {0, 1, 4}, {1, 3, 4}, {2, 4, 3}, {0, 4, 2}}; |
| 772 | 6400 | break; | |
| 773 |
1/2✓ Branch 1 taken 4720 times.
✗ Branch 2 not taken.
|
4720 | case 4: // tetrahedron |
| 774 |
5/12✓ Branch 1 taken 4720 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4720 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 18880 times.
✓ Branch 7 taken 4720 times.
✓ Branch 8 taken 18880 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
33040 | facetCorners = {{1, 0, 2}, {0, 1, 3}, {0, 3, 2}, {1, 2, 3}}; |
| 775 | 4720 | break; | |
| 776 | ✗ | default: | |
| 777 |
0/38✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 30 not taken.
✗ Branch 33 not taken.
✗ Branch 36 not taken.
✗ Branch 39 not taken.
✗ Branch 42 not taken.
✗ Branch 45 not taken.
|
160297 | DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type " |
| 778 | << geo1.type() << ", "<< geo1.corners() << " corners."); | ||
| 779 | } | ||
| 780 | |||
| 781 | 200412 | return facetCorners; | |
| 782 | ✗ | }; | |
| 783 | |||
| 784 | // lambda to obtain the normal on a facet | ||
| 785 | 1538611 | auto computeNormal = [&geo1] (const std::vector<int>& facetCorners) | |
| 786 | { | ||
| 787 | 2024767 | const auto v0 = geo1.corner(facetCorners[1]) - geo1.corner(facetCorners[0]); | |
| 788 | 2024767 | const auto v1 = geo1.corner(facetCorners[2]) - geo1.corner(facetCorners[0]); | |
| 789 | |||
| 790 | 816791 | auto n = crossProduct(v0, v1); | |
| 791 | 816791 | n /= n.two_norm(); | |
| 792 | |||
| 793 | 816791 | return n; | |
| 794 | }; | ||
| 795 | |||
| 796 | 200412 | return Detail::computeSegmentIntersection(geo1, geo2, eps_, tfirst, tlast, getFacetCorners, computeNormal); | |
| 797 | } | ||
| 798 | }; | ||
| 799 | |||
| 800 | /*! | ||
| 801 | * \ingroup Geometry | ||
| 802 | * \brief A class for segment--polyhedron intersection in 3d space | ||
| 803 | */ | ||
| 804 | template <class Geometry1, class Geometry2, class Policy> | ||
| 805 | class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 1, 3> | ||
| 806 | : public GeometryIntersection<Geometry2, Geometry1, Policy, 3, 3, 1> | ||
| 807 | { | ||
| 808 | using Base = GeometryIntersection<Geometry2, Geometry1, Policy, 3, 3, 1>; | ||
| 809 | public: | ||
| 810 | /*! | ||
| 811 | * \brief Colliding segment and convex polyhedron | ||
| 812 | * \param geo1/geo2 The geometries to intersect | ||
| 813 | * \param intersection If the geometries collide intersection holds the | ||
| 814 | * corner points of the intersection object in global coordinates. | ||
| 815 | * \note This forwards to the polyhedron-segment specialization with swapped arguments. | ||
| 816 | */ | ||
| 817 | template<class P = Policy> | ||
| 818 | 145753 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename Base::Intersection& intersection) | |
| 819 | { | ||
| 820 |
1/2✓ Branch 1 taken 143776 times.
✗ Branch 2 not taken.
|
145753 | return Base::intersection(geo2, geo1, intersection); |
| 821 | } | ||
| 822 | }; | ||
| 823 | |||
| 824 | /*! | ||
| 825 | * \ingroup Geometry | ||
| 826 | * \brief A class for polygon--polygon intersections in 3d space | ||
| 827 | */ | ||
| 828 | template <class Geometry1, class Geometry2, class Policy> | ||
| 829 | class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 2, 2> | ||
| 830 | { | ||
| 831 | enum { dimworld = 3 }; | ||
| 832 | enum { dim1 = 2 }; | ||
| 833 | enum { dim2 = 2 }; | ||
| 834 | |||
| 835 | public: | ||
| 836 | using ctype = typename Policy::ctype; | ||
| 837 | using Point = typename Policy::Point; | ||
| 838 | using Intersection = typename Policy::Intersection; | ||
| 839 | |||
| 840 | private: | ||
| 841 | static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons | ||
| 842 | |||
| 843 | public: | ||
| 844 | /*! | ||
| 845 | * \brief Colliding two polygons | ||
| 846 | * \note First we find the vertex candidates for the intersection region as follows: | ||
| 847 | * Add vertices of first polygon that are inside the second polygon | ||
| 848 | * Add intersections of polygon edges | ||
| 849 | * Remove duplicate points from the list | ||
| 850 | * Compute the convex hull polygon | ||
| 851 | * \param geo1/geo2 The geometries to intersect | ||
| 852 | * \param intersection Container to store the corner points of the polygon (as convex hull) | ||
| 853 | * \note This overload is used when polygon like intersections are seeked | ||
| 854 | */ | ||
| 855 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 2, int> = 0> | ||
| 856 | 49 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | |
| 857 | { | ||
| 858 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 859 | |||
| 860 | // if the geometries are not parallel, there cannot be a polygon intersection | ||
| 861 | 49 | const auto a1 = geo1.corner(1) - geo1.corner(0); | |
| 862 | 49 | const auto b1 = geo1.corner(2) - geo1.corner(0); | |
| 863 | 49 | const auto n1 = crossProduct(a1, b1); | |
| 864 | |||
| 865 | 49 | const auto a2 = geo2.corner(1) - geo2.corner(0); | |
| 866 | 49 | const auto b2 = geo2.corner(2) - geo2.corner(0); | |
| 867 | 49 | const auto n2 = crossProduct(a2, b2); | |
| 868 | |||
| 869 | // compute an epsilon scaled with larger edge length | ||
| 870 | 49 | const auto a1Norm2 = a1.two_norm2(); | |
| 871 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 1 times.
|
49 | const auto b1Norm2 = b1.two_norm2(); |
| 872 | |||
| 873 | using std::max; | ||
| 874 | 49 | const auto maxNorm2 = max(a1Norm2, b1Norm2); | |
| 875 | 49 | const auto eps2 = maxNorm2*eps_; | |
| 876 | |||
| 877 | // early exit if polygons are not parallel | ||
| 878 | using std::abs; | ||
| 879 |
2/2✓ Branch 1 taken 2 times.
✓ Branch 2 taken 47 times.
|
98 | if (crossProduct(n1, n2).two_norm2() > eps2*maxNorm2) |
| 880 | return false; | ||
| 881 | |||
| 882 | // they are parallel, verify that they are on the same plane | ||
| 883 | 94 | auto d1 = geo2.corner(0) - geo1.corner(0); | |
| 884 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 37 times.
|
47 | if (d1.two_norm2() < eps2) |
| 885 | 20 | d1 = geo2.corner(1) - geo1.corner(0); | |
| 886 | |||
| 887 | using std::sqrt; | ||
| 888 | 47 | const auto maxNorm = sqrt(maxNorm2); | |
| 889 | 47 | const auto eps3 = maxNorm*eps2; | |
| 890 | |||
| 891 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 43 times.
|
47 | if (abs(d1*n2) > eps3) |
| 892 | return false; | ||
| 893 | |||
| 894 | // the candidate intersection points | ||
| 895 |
1/2✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
|
43 | std::vector<Point> points; points.reserve(8); |
| 896 | |||
| 897 | // add polygon1 corners that are inside polygon2 | ||
| 898 |
2/2✓ Branch 1 taken 130 times.
✓ Branch 2 taken 43 times.
|
173 | for (int i = 0; i < geo1.corners(); ++i) |
| 899 |
3/4✓ Branch 2 taken 130 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 50 times.
✓ Branch 5 taken 80 times.
|
130 | if (intersectsPointGeometry(geo1.corner(i), geo2)) |
| 900 |
1/2✓ Branch 2 taken 50 times.
✗ Branch 3 not taken.
|
50 | points.emplace_back(geo1.corner(i)); |
| 901 | |||
| 902 |
2/2✓ Branch 0 taken 38 times.
✓ Branch 1 taken 5 times.
|
43 | const auto numPoints1 = points.size(); |
| 903 | 43 | const bool resultIsGeo1 = numPoints1 == geo1.corners(); | |
| 904 |
2/2✓ Branch 0 taken 38 times.
✓ Branch 1 taken 5 times.
|
43 | if (!resultIsGeo1) |
| 905 | { | ||
| 906 | // add polygon2 corners that are inside polygon1 | ||
| 907 |
2/2✓ Branch 1 taken 126 times.
✓ Branch 2 taken 38 times.
|
164 | for (int i = 0; i < geo2.corners(); ++i) |
| 908 |
3/4✓ Branch 2 taken 126 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 50 times.
✓ Branch 5 taken 76 times.
|
126 | if (intersectsPointGeometry(geo2.corner(i), geo1)) |
| 909 |
1/2✓ Branch 2 taken 50 times.
✗ Branch 3 not taken.
|
50 | points.emplace_back(geo2.corner(i)); |
| 910 | |||
| 911 |
1/2✓ Branch 0 taken 38 times.
✗ Branch 1 not taken.
|
38 | const bool resultIsGeo2 = (points.size() - numPoints1) == geo2.corners(); |
| 912 |
1/2✓ Branch 0 taken 38 times.
✗ Branch 1 not taken.
|
38 | if (!resultIsGeo2) |
| 913 | { | ||
| 914 | 38 | const auto referenceElement1 = referenceElement(geo1); | |
| 915 | 38 | const auto referenceElement2 = referenceElement(geo2); | |
| 916 | |||
| 917 | // add intersections of edges | ||
| 918 | using SegGeometry = Dune::MultiLinearGeometry<ctype, 1, dimworld>; | ||
| 919 | using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>; | ||
| 920 |
2/2✓ Branch 1 taken 114 times.
✓ Branch 2 taken 38 times.
|
152 | for (int i = 0; i < referenceElement1.size(dim1-1); ++i) |
| 921 | { | ||
| 922 | 114 | const auto localEdgeGeom1 = referenceElement1.template geometry<dim1-1>(i); | |
| 923 |
1/2✓ Branch 0 taken 114 times.
✗ Branch 1 not taken.
|
114 | const auto edge1 = SegGeometry( |
| 924 | Dune::GeometryTypes::line, | ||
| 925 |
1/2✓ Branch 1 taken 114 times.
✗ Branch 2 not taken.
|
114 | std::vector<Point>({ |
| 926 | 114 | geo1.global(localEdgeGeom1.corner(0)), | |
| 927 | 114 | geo1.global(localEdgeGeom1.corner(1)) | |
| 928 | }) | ||
| 929 | ); | ||
| 930 | |||
| 931 |
2/2✓ Branch 1 taken 378 times.
✓ Branch 2 taken 114 times.
|
492 | for (int j = 0; j < referenceElement2.size(dim2-1); ++j) |
| 932 | { | ||
| 933 | 378 | const auto localEdgeGeom2 = referenceElement2.template geometry<dim2-1>(j); | |
| 934 |
1/2✓ Branch 0 taken 378 times.
✗ Branch 1 not taken.
|
378 | const auto edge2 = SegGeometry( |
| 935 | Dune::GeometryTypes::line, | ||
| 936 |
1/4✓ Branch 1 taken 378 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
378 | std::vector<Point>({ |
| 937 | 378 | geo2.global(localEdgeGeom2.corner(0)), | |
| 938 | 378 | geo2.global(localEdgeGeom2.corner(1)) | |
| 939 | }) | ||
| 940 | ); | ||
| 941 | |||
| 942 | using EdgeTest = GeometryIntersection<SegGeometry, SegGeometry, PointPolicy>; | ||
| 943 | 378 | typename EdgeTest::Intersection edgeIntersection; | |
| 944 |
2/2✓ Branch 1 taken 154 times.
✓ Branch 2 taken 224 times.
|
378 | if (EdgeTest::intersection(edge1, edge2, edgeIntersection)) |
| 945 |
1/2✓ Branch 1 taken 154 times.
✗ Branch 2 not taken.
|
154 | points.emplace_back(edgeIntersection); |
| 946 | } | ||
| 947 | } | ||
| 948 | } | ||
| 949 | } | ||
| 950 | |||
| 951 |
1/2✓ Branch 0 taken 43 times.
✗ Branch 1 not taken.
|
43 | if (points.empty()) |
| 952 | return false; | ||
| 953 | |||
| 954 | // remove duplicates | ||
| 955 | 43 | const auto eps = maxNorm*eps_; | |
| 956 |
4/24✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 60 times.
✓ Branch 5 taken 248 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 71 times.
✓ Branch 23 taken 140 times.
|
562 | std::sort(points.begin(), points.end(), [eps] (const auto& a, const auto& b) -> bool |
| 957 | { | ||
| 958 | using std::abs; | ||
| 959 |
4/24✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 60 times.
✓ Branch 5 taken 248 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 71 times.
✓ Branch 23 taken 140 times.
|
907 | return (abs(a[0]-b[0]) > eps ? a[0] < b[0] |
| 960 |
4/24✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 84 times.
✓ Branch 5 taken 164 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 48 times.
✓ Branch 23 taken 92 times.
|
388 | : (abs(a[1]-b[1]) > eps ? a[1] < b[1] |
| 961 | 256 | : (a[2] < b[2]))); | |
| 962 | }); | ||
| 963 | |||
| 964 | 43 | const auto squaredEps = eps*eps; | |
| 965 |
2/2✓ Branch 2 taken 21 times.
✓ Branch 3 taken 22 times.
|
43 | points.erase(std::unique( |
| 966 | points.begin(), points.end(), | ||
| 967 | 422 | [squaredEps] (const auto& a, const auto&b) { return (b-a).two_norm2() < squaredEps; }), | |
| 968 | 43 | points.end() | |
| 969 | ); | ||
| 970 | |||
| 971 | // return false if we don't have at least three unique points | ||
| 972 |
2/2✓ Branch 0 taken 21 times.
✓ Branch 1 taken 22 times.
|
43 | if (points.size() < 3) |
| 973 | return false; | ||
| 974 | |||
| 975 | // intersection polygon is convex hull of above points | ||
| 976 |
2/6✗ Branch 1 not taken.
✓ Branch 2 taken 21 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
21 | intersection = grahamConvexHull<2>(points); |
| 977 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
|
21 | assert(!intersection.empty()); |
| 978 | return true; | ||
| 979 | 43 | } | |
| 980 | |||
| 981 | /*! | ||
| 982 | * \brief Colliding two polygons | ||
| 983 | * \param geo1/geo2 The geometries to intersect | ||
| 984 | * \param intersection Container to store the corners of intersection segment | ||
| 985 | * \note this overload is used when segment-like intersections are seeked | ||
| 986 | * \todo implement this query | ||
| 987 | */ | ||
| 988 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0> | ||
| 989 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | ||
| 990 | { | ||
| 991 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 992 | DUNE_THROW(Dune::NotImplemented, "Polygon-polygon intersection detection for segment-like intersections"); | ||
| 993 | } | ||
| 994 | |||
| 995 | /*! | ||
| 996 | * \brief Colliding two polygons | ||
| 997 | * \param geo1/geo2 The geometries to intersect | ||
| 998 | * \param intersection The intersection point | ||
| 999 | * \note this overload is used when point-like intersections are seeked | ||
| 1000 | * \todo implement this query | ||
| 1001 | */ | ||
| 1002 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0> | ||
| 1003 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | ||
| 1004 | { | ||
| 1005 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 1006 | DUNE_THROW(Dune::NotImplemented, "Polygon-polygon intersection detection for touching points"); | ||
| 1007 | } | ||
| 1008 | }; | ||
| 1009 | |||
| 1010 | /*! | ||
| 1011 | * \ingroup Geometry | ||
| 1012 | * \brief A class for polyhedron--polygon intersection in 3d space | ||
| 1013 | */ | ||
| 1014 | template <class Geometry1, class Geometry2, class Policy> | ||
| 1015 | class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 3, 2> | ||
| 1016 | { | ||
| 1017 | enum { dimworld = 3 }; | ||
| 1018 | enum { dim1 = 3 }; | ||
| 1019 | enum { dim2 = 2 }; | ||
| 1020 | |||
| 1021 | public: | ||
| 1022 | using ctype = typename Policy::ctype; | ||
| 1023 | using Point = typename Policy::Point; | ||
| 1024 | using Intersection = typename Policy::Intersection; | ||
| 1025 | |||
| 1026 | private: | ||
| 1027 | static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons | ||
| 1028 | |||
| 1029 | public: | ||
| 1030 | /*! | ||
| 1031 | * \brief Colliding polygon and convex polyhedron | ||
| 1032 | * \note First we find the vertex candidates for the intersection region as follows: | ||
| 1033 | * Add triangle vertices that are inside the tetrahedron | ||
| 1034 | * Add tetrahedron vertices that are inside the triangle | ||
| 1035 | * Add all intersection points of tetrahedron edges (codim 2) with the triangle (codim 0) (6*1 tests) | ||
| 1036 | * Add all intersection points of triangle edges (codim 1) with tetrahedron faces (codim 1) (3*4 tests) | ||
| 1037 | * Remove duplicate points from the list | ||
| 1038 | * Compute the convex hull polygon | ||
| 1039 | * Return a triangulation of that polygon as intersection | ||
| 1040 | * \param geo1/geo2 The geometries to intersect | ||
| 1041 | * \param intersection Container to store the corner points of the polygon (as convex hull) | ||
| 1042 | * \note This overload is used when polygon like intersections are seeked | ||
| 1043 | */ | ||
| 1044 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 2, int> = 0> | ||
| 1045 | 42719 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | |
| 1046 | { | ||
| 1047 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 1048 | |||
| 1049 | // the candidate intersection points | ||
| 1050 |
1/2✓ Branch 1 taken 35111 times.
✗ Branch 2 not taken.
|
42719 | std::vector<Point> points; points.reserve(10); |
| 1051 | |||
| 1052 | // add 3d geometry corners that are inside the 2d geometry | ||
| 1053 |
3/3✓ Branch 0 taken 253136 times.
✓ Branch 1 taken 35251 times.
✓ Branch 2 taken 20 times.
|
356879 | for (int i = 0; i < geo1.corners(); ++i) |
| 1054 |
5/5✓ Branch 2 taken 8616 times.
✓ Branch 3 taken 18918 times.
✓ Branch 4 taken 221370 times.
✓ Branch 5 taken 144 times.
✓ Branch 1 taken 244680 times.
|
314160 | if (intersectsPointGeometry(geo1.corner(i), geo2)) |
| 1055 |
2/3✓ Branch 1 taken 23326 times.
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
|
34158 | points.emplace_back(geo1.corner(i)); |
| 1056 | |||
| 1057 | // add 2d geometry corners that are inside the 3d geometry | ||
| 1058 |
3/3✓ Branch 0 taken 6432 times.
✓ Branch 1 taken 102133 times.
✓ Branch 2 taken 33503 times.
|
172696 | for (int i = 0; i < geo2.corners(); ++i) |
| 1059 |
5/5✓ Branch 2 taken 90579 times.
✓ Branch 3 taken 13348 times.
✓ Branch 4 taken 16 times.
✓ Branch 5 taken 60 times.
✓ Branch 1 taken 3030 times.
|
129977 | if (intersectsPointGeometry(geo2.corner(i), geo1)) |
| 1060 |
2/3✓ Branch 1 taken 3030 times.
✓ Branch 2 taken 8130 times.
✗ Branch 3 not taken.
|
11480 | points.emplace_back(geo2.corner(i)); |
| 1061 | |||
| 1062 | // get some geometry types | ||
| 1063 | using PolyhedronFaceGeometry = Dune::MultiLinearGeometry<ctype, 2, dimworld>; | ||
| 1064 | using SegGeometry = Dune::MultiLinearGeometry<ctype, 1, dimworld>; | ||
| 1065 | |||
| 1066 |
1/2✓ Branch 1 taken 23581 times.
✗ Branch 2 not taken.
|
42719 | const auto refElement1 = referenceElement(geo1); |
| 1067 | 42719 | const auto refElement2 = referenceElement(geo2); | |
| 1068 | |||
| 1069 | // intersection policy for point-like intersections (used later) | ||
| 1070 | using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>; | ||
| 1071 | |||
| 1072 | // add intersection points of all polyhedron edges (codim dim-1) with the polygon | ||
| 1073 |
2/2✓ Branch 1 taken 379944 times.
✓ Branch 2 taken 35111 times.
|
513959 | for (int i = 0; i < refElement1.size(dim1-1); ++i) |
| 1074 | { | ||
| 1075 | 471240 | const auto localEdgeGeom = refElement1.template geometry<dim1-1>(i); | |
| 1076 | 900852 | const auto p = geo1.global(localEdgeGeom.corner(0)); | |
| 1077 | 471240 | const auto q = geo1.global(localEdgeGeom.corner(1)); | |
| 1078 |
2/4✓ Branch 1 taken 379944 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 379944 times.
✗ Branch 4 not taken.
|
942480 | const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q}); |
| 1079 | |||
| 1080 | using PolySegTest = GeometryIntersection<Geometry2, SegGeometry, PointPolicy>; | ||
| 1081 | 471240 | typename PolySegTest::Intersection polySegIntersection; | |
| 1082 |
3/4✓ Branch 1 taken 379944 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 53193 times.
✓ Branch 4 taken 326751 times.
|
471240 | if (PolySegTest::intersection(geo2, segGeo, polySegIntersection)) |
| 1083 |
1/2✓ Branch 1 taken 53193 times.
✗ Branch 2 not taken.
|
67145 | points.emplace_back(std::move(polySegIntersection)); |
| 1084 | } | ||
| 1085 | |||
| 1086 | // add intersection points of all polygon faces (codim 1) with the polyhedron faces | ||
| 1087 |
2/2✓ Branch 1 taken 196870 times.
✓ Branch 2 taken 35111 times.
|
285237 | for (int i = 0; i < refElement1.size(1); ++i) |
| 1088 | { | ||
| 1089 |
1/2✓ Branch 1 taken 196870 times.
✗ Branch 2 not taken.
|
288166 | const auto faceGeo = [&]() |
| 1090 | { | ||
| 1091 | 196870 | const auto localFaceGeo = refElement1.template geometry<1>(i); | |
| 1092 |
5/8✓ Branch 0 taken 124806 times.
✓ Branch 1 taken 27592 times.
✓ Branch 2 taken 1536 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 26136 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 16800 times.
✗ Branch 7 not taken.
|
196870 | if (localFaceGeo.corners() == 4) |
| 1093 | { | ||
| 1094 | 169278 | const auto a = geo1.global(localFaceGeo.corner(0)); | |
| 1095 | 169278 | const auto b = geo1.global(localFaceGeo.corner(1)); | |
| 1096 | 169278 | const auto c = geo1.global(localFaceGeo.corner(2)); | |
| 1097 | 169278 | const auto d = geo1.global(localFaceGeo.corner(3)); | |
| 1098 | |||
| 1099 |
12/24✓ Branch 1 taken 124806 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 124806 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 124806 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1536 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 1536 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1536 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 26136 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 26136 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 26136 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 16800 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 16800 times.
✗ Branch 29 not taken.
✓ Branch 30 taken 16800 times.
✗ Branch 31 not taken.
|
507834 | return PolyhedronFaceGeometry(Dune::GeometryTypes::cube(2), std::vector<Point>{a, b, c, d}); |
| 1100 | } | ||
| 1101 | else | ||
| 1102 | { | ||
| 1103 | 27592 | const auto a = geo1.global(localFaceGeo.corner(0)); | |
| 1104 | 27592 | const auto b = geo1.global(localFaceGeo.corner(1)); | |
| 1105 | 27592 | const auto c = geo1.global(localFaceGeo.corner(2)); | |
| 1106 | |||
| 1107 |
3/24✓ Branch 1 taken 27592 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27592 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 27592 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
|
82776 | return PolyhedronFaceGeometry(Dune::GeometryTypes::simplex(2), std::vector<Point>{a, b, c}); |
| 1108 | } | ||
| 1109 | }(); | ||
| 1110 | |||
| 1111 |
2/2✓ Branch 1 taken 600354 times.
✓ Branch 2 taken 196870 times.
|
980992 | for (int j = 0; j < refElement2.size(1); ++j) |
| 1112 | { | ||
| 1113 | 738474 | const auto localEdgeGeom = refElement2.template geometry<1>(j); | |
| 1114 | 738474 | const auto p = geo2.global(localEdgeGeom.corner(0)); | |
| 1115 | 738474 | const auto q = geo2.global(localEdgeGeom.corner(1)); | |
| 1116 | |||
| 1117 |
2/6✓ Branch 1 taken 600354 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 600354 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
1476948 | const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q}); |
| 1118 | |||
| 1119 | using PolySegTest = GeometryIntersection<PolyhedronFaceGeometry, SegGeometry, PointPolicy>; | ||
| 1120 | 738474 | typename PolySegTest::Intersection polySegIntersection; | |
| 1121 |
3/4✓ Branch 1 taken 600354 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 66202 times.
✓ Branch 4 taken 534152 times.
|
738474 | if (PolySegTest::intersection(faceGeo, segGeo, polySegIntersection)) |
| 1122 |
1/2✓ Branch 1 taken 66202 times.
✗ Branch 2 not taken.
|
71674 | points.emplace_back(std::move(polySegIntersection)); |
| 1123 | } | ||
| 1124 | } | ||
| 1125 | |||
| 1126 | // return if no intersection points were found | ||
| 1127 |
2/2✓ Branch 0 taken 20324 times.
✓ Branch 1 taken 14787 times.
|
42719 | if (points.empty()) return false; |
| 1128 | |||
| 1129 | // remove duplicates | ||
| 1130 | 66446 | const auto eps = (geo1.corner(0) - geo1.corner(1)).two_norm()*eps_; | |
| 1131 |
69/192✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 74023 times.
✓ Branch 9 taken 131455 times.
✓ Branch 10 taken 18007 times.
✓ Branch 11 taken 113448 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✓ Branch 32 taken 368 times.
✓ Branch 33 taken 364 times.
✓ Branch 34 taken 364 times.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✓ Branch 44 taken 33178 times.
✓ Branch 45 taken 73723 times.
✓ Branch 46 taken 12784 times.
✓ Branch 47 taken 60939 times.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
✓ Branch 68 taken 720 times.
✓ Branch 69 taken 1774 times.
✗ Branch 70 not taken.
✓ Branch 71 taken 1774 times.
✓ Branch 72 taken 78 times.
✓ Branch 73 taken 69 times.
✗ Branch 74 not taken.
✓ Branch 75 taken 69 times.
✓ Branch 76 taken 35 times.
✓ Branch 77 taken 23 times.
✗ Branch 78 not taken.
✓ Branch 79 taken 23 times.
✓ Branch 80 taken 1 times.
✓ Branch 81 taken 9 times.
✗ Branch 82 not taken.
✓ Branch 83 taken 9 times.
✓ Branch 84 taken 4 times.
✓ Branch 85 taken 4 times.
✗ Branch 86 not taken.
✓ Branch 87 taken 4 times.
✓ Branch 88 taken 4 times.
✓ Branch 89 taken 2 times.
✗ Branch 90 not taken.
✓ Branch 91 taken 2 times.
✓ Branch 92 taken 1 times.
✓ Branch 93 taken 1 times.
✗ Branch 94 not taken.
✓ Branch 95 taken 1 times.
✓ Branch 96 taken 1 times.
✗ Branch 97 not taken.
✗ Branch 98 not taken.
✗ Branch 99 not taken.
✗ Branch 100 not taken.
✗ Branch 101 not taken.
✗ Branch 102 not taken.
✗ Branch 103 not taken.
✓ Branch 104 taken 18140 times.
✓ Branch 105 taken 21376 times.
✓ Branch 106 taken 10304 times.
✓ Branch 107 taken 11072 times.
✗ Branch 108 not taken.
✗ Branch 109 not taken.
✗ Branch 110 not taken.
✗ Branch 111 not taken.
✗ Branch 112 not taken.
✗ Branch 113 not taken.
✗ Branch 114 not taken.
✗ Branch 115 not taken.
✗ Branch 116 not taken.
✗ Branch 117 not taken.
✗ Branch 118 not taken.
✗ Branch 119 not taken.
✗ Branch 120 not taken.
✗ Branch 121 not taken.
✗ Branch 122 not taken.
✗ Branch 123 not taken.
✗ Branch 124 not taken.
✗ Branch 125 not taken.
✗ Branch 126 not taken.
✗ Branch 127 not taken.
✗ Branch 128 not taken.
✗ Branch 129 not taken.
✗ Branch 130 not taken.
✗ Branch 131 not taken.
✗ Branch 132 not taken.
✗ Branch 133 not taken.
✗ Branch 134 not taken.
✗ Branch 135 not taken.
✗ Branch 136 not taken.
✗ Branch 137 not taken.
✗ Branch 138 not taken.
✗ Branch 139 not taken.
✓ Branch 140 taken 6399 times.
✓ Branch 141 taken 7477 times.
✓ Branch 142 taken 1861 times.
✓ Branch 143 taken 5616 times.
✓ Branch 144 taken 19 times.
✓ Branch 145 taken 32 times.
✓ Branch 146 taken 23 times.
✓ Branch 147 taken 9 times.
✓ Branch 148 taken 1 times.
✓ Branch 149 taken 28 times.
✓ Branch 150 taken 21 times.
✓ Branch 151 taken 7 times.
✗ Branch 152 not taken.
✓ Branch 153 taken 4 times.
✗ Branch 154 not taken.
✓ Branch 155 taken 4 times.
✗ Branch 156 not taken.
✗ Branch 157 not taken.
✗ Branch 158 not taken.
✗ Branch 159 not taken.
✗ Branch 160 not taken.
✗ Branch 161 not taken.
✗ Branch 162 not taken.
✗ Branch 163 not taken.
✗ Branch 164 not taken.
✓ Branch 165 taken 4 times.
✗ Branch 166 not taken.
✓ Branch 167 taken 4 times.
✗ Branch 168 not taken.
✓ Branch 169 taken 4 times.
✗ Branch 170 not taken.
✓ Branch 171 taken 4 times.
✗ Branch 172 not taken.
✗ Branch 173 not taken.
✗ Branch 174 not taken.
✗ Branch 175 not taken.
✓ Branch 176 taken 372 times.
✓ Branch 177 taken 216 times.
✓ Branch 178 taken 216 times.
✗ Branch 179 not taken.
✓ Branch 180 taken 206 times.
✓ Branch 181 taken 866 times.
✗ Branch 182 not taken.
✓ Branch 183 taken 866 times.
✓ Branch 184 taken 9500 times.
✓ Branch 185 taken 8172 times.
✓ Branch 186 taken 4756 times.
✓ Branch 187 taken 3416 times.
✓ Branch 188 taken 3035 times.
✓ Branch 189 taken 4261 times.
✓ Branch 190 taken 818 times.
✓ Branch 191 taken 3443 times.
|
403725 | const auto notEqual = [eps] (auto a, auto b) { using std::abs; return abs(b-a) > eps; }; |
| 1132 | 107782 | std::sort(points.begin(), points.end(), [notEqual](const auto& a, const auto& b) -> bool | |
| 1133 | { | ||
| 1134 |
40/96✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 74023 times.
✓ Branch 5 taken 131455 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 368 times.
✓ Branch 17 taken 364 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 33178 times.
✓ Branch 23 taken 73723 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✓ Branch 34 taken 720 times.
✓ Branch 35 taken 1774 times.
✓ Branch 36 taken 78 times.
✓ Branch 37 taken 69 times.
✓ Branch 38 taken 35 times.
✓ Branch 39 taken 23 times.
✓ Branch 40 taken 1 times.
✓ Branch 41 taken 9 times.
✓ Branch 42 taken 4 times.
✓ Branch 43 taken 4 times.
✓ Branch 44 taken 4 times.
✓ Branch 45 taken 2 times.
✓ Branch 46 taken 1 times.
✓ Branch 47 taken 1 times.
✓ Branch 48 taken 1 times.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✓ Branch 52 taken 18140 times.
✓ Branch 53 taken 21376 times.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✓ Branch 70 taken 6399 times.
✓ Branch 71 taken 7477 times.
✓ Branch 72 taken 19 times.
✓ Branch 73 taken 32 times.
✓ Branch 74 taken 1 times.
✓ Branch 75 taken 28 times.
✗ Branch 76 not taken.
✓ Branch 77 taken 4 times.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✓ Branch 83 taken 4 times.
✗ Branch 84 not taken.
✓ Branch 85 taken 4 times.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✓ Branch 88 taken 372 times.
✓ Branch 89 taken 216 times.
✓ Branch 90 taken 206 times.
✓ Branch 91 taken 866 times.
✓ Branch 92 taken 9500 times.
✓ Branch 93 taken 8172 times.
✓ Branch 94 taken 3035 times.
✓ Branch 95 taken 4261 times.
|
645813 | return (notEqual(a[0], b[0]) ? a[0] < b[0] |
| 1135 |
29/96✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 18007 times.
✓ Branch 5 taken 113448 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 364 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 12784 times.
✓ Branch 23 taken 60939 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✓ Branch 35 taken 1774 times.
✗ Branch 36 not taken.
✓ Branch 37 taken 69 times.
✗ Branch 38 not taken.
✓ Branch 39 taken 23 times.
✗ Branch 40 not taken.
✓ Branch 41 taken 9 times.
✗ Branch 42 not taken.
✓ Branch 43 taken 4 times.
✗ Branch 44 not taken.
✓ Branch 45 taken 2 times.
✗ Branch 46 not taken.
✓ Branch 47 taken 1 times.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✓ Branch 52 taken 10304 times.
✓ Branch 53 taken 11072 times.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✓ Branch 70 taken 1861 times.
✓ Branch 71 taken 5616 times.
✓ Branch 72 taken 23 times.
✓ Branch 73 taken 9 times.
✓ Branch 74 taken 21 times.
✓ Branch 75 taken 7 times.
✗ Branch 76 not taken.
✓ Branch 77 taken 4 times.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✓ Branch 83 taken 4 times.
✗ Branch 84 not taken.
✓ Branch 85 taken 4 times.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✓ Branch 88 taken 216 times.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✓ Branch 91 taken 866 times.
✓ Branch 92 taken 4756 times.
✓ Branch 93 taken 3416 times.
✓ Branch 94 taken 818 times.
✓ Branch 95 taken 3443 times.
|
249864 | : (notEqual(a[1], b[1]) ? a[1] < b[1] |
| 1136 | 200710 | : (a[2] < b[2]))); | |
| 1137 | }); | ||
| 1138 | |||
| 1139 | 24212 | const auto squaredEps = eps*eps; | |
| 1140 |
2/2✓ Branch 2 taken 13496 times.
✓ Branch 3 taken 6828 times.
|
24212 | points.erase(std::unique( |
| 1141 | points.begin(), points.end(), | ||
| 1142 | 267146 | [squaredEps] (const auto& a, const auto&b) { return (b-a).two_norm2() < squaredEps; }), | |
| 1143 | 24212 | points.end() | |
| 1144 | ); | ||
| 1145 | |||
| 1146 | // return false if we don't have more than three unique points | ||
| 1147 |
2/2✓ Branch 0 taken 13496 times.
✓ Branch 1 taken 6828 times.
|
24212 | if (points.size() < 3) return false; |
| 1148 | |||
| 1149 | // intersection polygon is convex hull of above points | ||
| 1150 |
2/6✗ Branch 1 not taken.
✓ Branch 2 taken 13496 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 13496 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
16412 | intersection = grahamConvexHull<2>(points); |
| 1151 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 13496 times.
|
16412 | assert(!intersection.empty()); |
| 1152 | |||
| 1153 | return true; | ||
| 1154 | 42719 | } | |
| 1155 | |||
| 1156 | /*! | ||
| 1157 | * \brief Colliding segment and convex polyhedron | ||
| 1158 | * \note First we find the vertex candidates for the intersection region as follows: | ||
| 1159 | * Add triangle vertices that are inside the tetrahedron | ||
| 1160 | * Add tetrahedron vertices that are inside the triangle | ||
| 1161 | * Add all intersection points of tetrahedron edges (codim 2) with the triangle (codim 0) (6*1 tests) | ||
| 1162 | * Add all intersection points of triangle edges (codim 1) with tetrahedron faces (codim 1) (3*4 tests) | ||
| 1163 | * Remove duplicate points from the list | ||
| 1164 | * Compute the convex hull polygon | ||
| 1165 | * Return a triangulation of that polygon as intersection | ||
| 1166 | * \param geo1/geo2 The geometries to intersect | ||
| 1167 | * \param intersection Container to store the intersection result | ||
| 1168 | * \todo implement overloads for segment or point-like intersections | ||
| 1169 | */ | ||
| 1170 | template<class P = Policy, std::enable_if_t<P::dimIntersection != 2, int> = 0> | ||
| 1171 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | ||
| 1172 | { | ||
| 1173 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 1174 | DUNE_THROW(Dune::NotImplemented, "Polyhedron-polygon intersection detection only implemented for polygon-like intersections"); | ||
| 1175 | } | ||
| 1176 | }; | ||
| 1177 | |||
| 1178 | /*! | ||
| 1179 | * \ingroup Geometry | ||
| 1180 | * \brief A class for polygon--polyhedron intersection in 3d space | ||
| 1181 | */ | ||
| 1182 | template <class Geometry1, class Geometry2, class Policy> | ||
| 1183 | class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 2, 3> | ||
| 1184 | : public GeometryIntersection<Geometry2, Geometry1, Policy, 3, 3, 2> | ||
| 1185 | { | ||
| 1186 | using Base = GeometryIntersection<Geometry2, Geometry1, Policy, 3, 3, 2>; | ||
| 1187 | public: | ||
| 1188 | /*! | ||
| 1189 | * \brief Colliding polygon and convex polyhedron | ||
| 1190 | * \param geo1/geo2 The geometries to intersect | ||
| 1191 | * \param intersection If the geometries collide intersection holds the | ||
| 1192 | * corner points of the intersection object in global coordinates. | ||
| 1193 | * \note This forwards to the polyhedron-polygon specialization with swapped arguments. | ||
| 1194 | */ | ||
| 1195 | template<class P = Policy> | ||
| 1196 | 4808 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename Base::Intersection& intersection) | |
| 1197 | { | ||
| 1198 |
3/6✓ Branch 1 taken 196 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 256 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4356 times.
✗ Branch 8 not taken.
|
4808 | return Base::intersection(geo2, geo1, intersection); |
| 1199 | } | ||
| 1200 | }; | ||
| 1201 | |||
| 1202 | /*! | ||
| 1203 | * \ingroup Geometry | ||
| 1204 | * \brief A class for polyhedron--polyhedron intersection in 3d space | ||
| 1205 | */ | ||
| 1206 | template <class Geometry1, class Geometry2, class Policy> | ||
| 1207 | class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 3, 3> | ||
| 1208 | { | ||
| 1209 | enum { dimworld = 3 }; | ||
| 1210 | enum { dim1 = 3 }; | ||
| 1211 | enum { dim2 = 3 }; | ||
| 1212 | |||
| 1213 | public: | ||
| 1214 | using ctype = typename Policy::ctype; | ||
| 1215 | using Point = typename Policy::Point; | ||
| 1216 | using Intersection = typename Policy::Intersection; | ||
| 1217 | |||
| 1218 | private: | ||
| 1219 | static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons | ||
| 1220 | |||
| 1221 | public: | ||
| 1222 | /*! | ||
| 1223 | * \brief Colliding two convex polyhedra | ||
| 1224 | * \note First we find the vertex candidates for the intersection region as follows: | ||
| 1225 | * Add vertices that are inside the other geometry for both geometries | ||
| 1226 | * Add all intersection points of edges (codim 2) with the other tetrahedron's faces triangle | ||
| 1227 | * Remove duplicate points from the list | ||
| 1228 | * Return a triangulation of the polyhedron formed by the convex hull of the point cloud | ||
| 1229 | * \param geo1/geo2 The geometries to intersect | ||
| 1230 | * \param intersection Container to store the corner points of the polygon (as convex hull) | ||
| 1231 | * \note This overload is used when polyhedron-like intersections are seeked | ||
| 1232 | */ | ||
| 1233 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 3, int> = 0> | ||
| 1234 | 14668 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | |
| 1235 | { | ||
| 1236 | static_assert( | ||
| 1237 | int(dimworld) == int(Geometry2::coorddimension), | ||
| 1238 | "Can only collide geometries of same coordinate dimension" | ||
| 1239 | ); | ||
| 1240 | |||
| 1241 | 14668 | const auto refElement1 = referenceElement(geo1); | |
| 1242 | 14668 | const auto refElement2 = referenceElement(geo2); | |
| 1243 | |||
| 1244 | // the candidate intersection points | ||
| 1245 | 14668 | std::vector<Point> points; | |
| 1246 |
1/2✓ Branch 2 taken 8339 times.
✗ Branch 3 not taken.
|
14668 | points.reserve(refElement1.size(2) + refElement2.size(2)); |
| 1247 | |||
| 1248 | // add corners inside the other geometry | ||
| 1249 | 55568 | const auto addPointIntersections = [&](const auto& g1, const auto& g2) | |
| 1250 | { | ||
| 1251 |
6/6✓ Branch 0 taken 18552 times.
✓ Branch 1 taken 4303 times.
✓ Branch 2 taken 1980 times.
✓ Branch 3 taken 15521 times.
✓ Branch 4 taken 96222 times.
✓ Branch 5 taken 11792 times.
|
148370 | for (int i = 0; i < g1.corners(); ++i) |
| 1252 |
8/8✓ Branch 1 taken 15194 times.
✓ Branch 2 taken 4010 times.
✓ Branch 4 taken 762 times.
✓ Branch 5 taken 986 times.
✓ Branch 7 taken 28008 times.
✓ Branch 8 taken 66328 times.
✓ Branch 6 taken 15072 times.
✓ Branch 3 taken 1332 times.
|
131692 | if (const auto& c = g1.corner(i); intersectsPointGeometry(c, g2)) |
| 1253 | 44632 | points.emplace_back(c); | |
| 1254 | }; | ||
| 1255 | |||
| 1256 |
1/2✓ Branch 1 taken 8339 times.
✗ Branch 2 not taken.
|
14668 | addPointIntersections(geo1, geo2); |
| 1257 |
1/2✓ Branch 1 taken 8339 times.
✗ Branch 2 not taken.
|
14668 | addPointIntersections(geo2, geo1); |
| 1258 | |||
| 1259 | // get geometry types for the facets | ||
| 1260 | using PolyhedronFaceGeometry = Dune::MultiLinearGeometry<ctype, 2, dimworld>; | ||
| 1261 | using SegGeometry = Dune::MultiLinearGeometry<ctype, 1, dimworld>; | ||
| 1262 | |||
| 1263 | // intersection policy for point-like intersections (used later) | ||
| 1264 | using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>; | ||
| 1265 | |||
| 1266 | // add intersection points of all edges with the faces of the other polyhedron | ||
| 1267 | 34994 | const auto addEdgeIntersections = [&](const auto& g1, const auto& g2, const auto& ref1, const auto& ref2) | |
| 1268 | { | ||
| 1269 |
6/6✓ Branch 1 taken 15402 times.
✓ Branch 2 taken 2567 times.
✓ Branch 4 taken 13048 times.
✓ Branch 5 taken 2319 times.
✓ Branch 7 taken 70752 times.
✓ Branch 8 taken 11792 times.
|
115880 | for (int i = 0; i < ref1.size(1); ++i) |
| 1270 | { | ||
| 1271 | 219548 | const auto faceGeo = [&]() | |
| 1272 | { | ||
| 1273 | 99202 | const auto localFaceGeo = ref1.template geometry<1>(i); | |
| 1274 |
4/6✓ Branch 0 taken 15402 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 11316 times.
✓ Branch 3 taken 1732 times.
✓ Branch 4 taken 70752 times.
✗ Branch 5 not taken.
|
99202 | if (localFaceGeo.corners() == 4) |
| 1275 | { | ||
| 1276 | 97470 | const auto a = g1.global(localFaceGeo.corner(0)); | |
| 1277 | 97470 | const auto b = g1.global(localFaceGeo.corner(1)); | |
| 1278 | 97470 | const auto c = g1.global(localFaceGeo.corner(2)); | |
| 1279 | 97470 | const auto d = g1.global(localFaceGeo.corner(3)); | |
| 1280 | |||
| 1281 |
3/6✓ Branch 1 taken 15402 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11316 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 70752 times.
✗ Branch 8 not taken.
|
194940 | return PolyhedronFaceGeometry( |
| 1282 | Dune::GeometryTypes::cube(2), std::vector<Point>{a, b, c, d} | ||
| 1283 |
6/12✓ Branch 1 taken 15402 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 15402 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 11316 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 11316 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 70752 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 70752 times.
✗ Branch 14 not taken.
|
194940 | ); |
| 1284 | } | ||
| 1285 | else | ||
| 1286 | { | ||
| 1287 | 1732 | const auto a = g1.global(localFaceGeo.corner(0)); | |
| 1288 | 1732 | const auto b = g1.global(localFaceGeo.corner(1)); | |
| 1289 | 1732 | const auto c = g1.global(localFaceGeo.corner(2)); | |
| 1290 | |||
| 1291 |
1/6✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 4 taken 1732 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
3464 | return PolyhedronFaceGeometry( |
| 1292 | Dune::GeometryTypes::simplex(2), std::vector<Point>{a, b, c} | ||
| 1293 |
2/12✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 6 taken 1732 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1732 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
|
3464 | ); |
| 1294 | } | ||
| 1295 | }(); | ||
| 1296 | |||
| 1297 |
6/6✓ Branch 1 taken 169236 times.
✓ Branch 2 taken 15402 times.
✓ Branch 4 taken 156576 times.
✓ Branch 5 taken 13048 times.
✓ Branch 7 taken 849024 times.
✓ Branch 8 taken 70752 times.
|
1274038 | for (int j = 0; j < ref2.size(2); ++j) |
| 1298 | { | ||
| 1299 | 1174836 | const auto localEdgeGeom = ref2.template geometry<2>(j); | |
| 1300 | 2180436 | const auto p = g2.global(localEdgeGeom.corner(0)); | |
| 1301 | 1174836 | const auto q = g2.global(localEdgeGeom.corner(1)); | |
| 1302 | |||
| 1303 |
6/18✓ Branch 1 taken 169236 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 169236 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 8 taken 156576 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 156576 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 15 taken 849024 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 849024 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
|
2349672 | const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q}); |
| 1304 | |||
| 1305 | using PolySegTest = GeometryIntersection<PolyhedronFaceGeometry, SegGeometry, PointPolicy>; | ||
| 1306 | 1174836 | typename PolySegTest::Intersection polySegIntersection; | |
| 1307 |
9/12✓ Branch 1 taken 169236 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6421 times.
✓ Branch 4 taken 162815 times.
✓ Branch 6 taken 156576 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 7442 times.
✓ Branch 9 taken 149134 times.
✓ Branch 11 taken 849024 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 109584 times.
✓ Branch 14 taken 739440 times.
|
1174836 | if (PolySegTest::intersection(faceGeo, segGeo, polySegIntersection)) |
| 1308 |
3/6✓ Branch 1 taken 6421 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7442 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 109584 times.
✗ Branch 8 not taken.
|
123447 | points.emplace_back(std::move(polySegIntersection)); |
| 1309 | } | ||
| 1310 | } | ||
| 1311 | }; | ||
| 1312 | |||
| 1313 |
1/2✓ Branch 1 taken 8339 times.
✗ Branch 2 not taken.
|
14668 | addEdgeIntersections(geo1, geo2, refElement1, refElement2); |
| 1314 |
1/2✓ Branch 1 taken 8339 times.
✗ Branch 2 not taken.
|
14668 | addEdgeIntersections(geo2, geo1, refElement2, refElement1); |
| 1315 | |||
| 1316 | // return if no intersection points were found | ||
| 1317 |
2/2✓ Branch 0 taken 46 times.
✓ Branch 1 taken 8293 times.
|
14668 | if (points.empty()) |
| 1318 | return false; | ||
| 1319 | |||
| 1320 | // remove duplicates | ||
| 1321 | 40944 | const auto norm = (geo1.corner(0) - geo1.corner(1)).two_norm(); | |
| 1322 | 14576 | const auto eps = norm*eps_; | |
| 1323 |
79/96✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 25784 times.
✓ Branch 9 taken 41571 times.
✓ Branch 10 taken 17208 times.
✓ Branch 11 taken 24363 times.
✓ Branch 12 taken 3395 times.
✓ Branch 13 taken 4965 times.
✓ Branch 14 taken 1994 times.
✓ Branch 15 taken 2971 times.
✓ Branch 16 taken 4989 times.
✓ Branch 17 taken 10435 times.
✓ Branch 18 taken 2902 times.
✓ Branch 19 taken 7533 times.
✓ Branch 20 taken 322 times.
✓ Branch 21 taken 1365 times.
✓ Branch 22 taken 238 times.
✓ Branch 23 taken 1127 times.
✓ Branch 24 taken 132 times.
✓ Branch 25 taken 1249 times.
✓ Branch 26 taken 93 times.
✓ Branch 27 taken 1156 times.
✓ Branch 28 taken 29 times.
✓ Branch 29 taken 196 times.
✓ Branch 30 taken 61 times.
✓ Branch 31 taken 135 times.
✓ Branch 32 taken 153 times.
✓ Branch 33 taken 207 times.
✓ Branch 34 taken 64 times.
✓ Branch 35 taken 143 times.
✓ Branch 36 taken 24 times.
✓ Branch 37 taken 220 times.
✓ Branch 38 taken 44 times.
✓ Branch 39 taken 176 times.
✓ Branch 40 taken 16 times.
✓ Branch 41 taken 80 times.
✓ Branch 42 taken 4 times.
✓ Branch 43 taken 76 times.
✓ Branch 44 taken 10657 times.
✓ Branch 45 taken 10347 times.
✓ Branch 46 taken 5680 times.
✓ Branch 47 taken 4667 times.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✓ Branch 52 taken 73198 times.
✓ Branch 53 taken 224692 times.
✓ Branch 54 taken 60960 times.
✓ Branch 55 taken 163732 times.
✓ Branch 56 taken 20280 times.
✓ Branch 57 taken 69291 times.
✓ Branch 58 taken 16279 times.
✓ Branch 59 taken 53012 times.
✓ Branch 60 taken 30500 times.
✓ Branch 61 taken 71604 times.
✓ Branch 62 taken 13436 times.
✓ Branch 63 taken 58168 times.
✓ Branch 64 taken 3065 times.
✓ Branch 65 taken 5142 times.
✓ Branch 66 taken 2140 times.
✓ Branch 67 taken 3002 times.
✓ Branch 68 taken 1745 times.
✓ Branch 69 taken 4305 times.
✓ Branch 70 taken 520 times.
✓ Branch 71 taken 3785 times.
✓ Branch 72 taken 1570 times.
✓ Branch 73 taken 3730 times.
✓ Branch 74 taken 1442 times.
✓ Branch 75 taken 2288 times.
✓ Branch 76 taken 625 times.
✓ Branch 77 taken 1532 times.
✓ Branch 78 taken 125 times.
✓ Branch 79 taken 1407 times.
✓ Branch 80 taken 125 times.
✓ Branch 81 taken 637 times.
✗ Branch 82 not taken.
✓ Branch 83 taken 637 times.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✓ Branch 88 taken 1398 times.
✓ Branch 89 taken 2111 times.
✓ Branch 90 taken 406 times.
✓ Branch 91 taken 1705 times.
✓ Branch 92 taken 18369 times.
✓ Branch 93 taken 69559 times.
✓ Branch 94 taken 17014 times.
✓ Branch 95 taken 52545 times.
|
732180 | const auto notEqual = [eps] (auto a, auto b) { using std::abs; return abs(b-a) > eps; }; |
| 1324 | 629196 | std::sort(points.begin(), points.end(), [notEqual](const auto& a, const auto& b) -> bool | |
| 1325 | { | ||
| 1326 |
40/48✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 25784 times.
✓ Branch 5 taken 41571 times.
✓ Branch 6 taken 3395 times.
✓ Branch 7 taken 4965 times.
✓ Branch 8 taken 4989 times.
✓ Branch 9 taken 10435 times.
✓ Branch 10 taken 322 times.
✓ Branch 11 taken 1365 times.
✓ Branch 12 taken 132 times.
✓ Branch 13 taken 1249 times.
✓ Branch 14 taken 29 times.
✓ Branch 15 taken 196 times.
✓ Branch 16 taken 153 times.
✓ Branch 17 taken 207 times.
✓ Branch 18 taken 24 times.
✓ Branch 19 taken 220 times.
✓ Branch 20 taken 16 times.
✓ Branch 21 taken 80 times.
✓ Branch 22 taken 10657 times.
✓ Branch 23 taken 10347 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 73198 times.
✓ Branch 27 taken 224692 times.
✓ Branch 28 taken 20280 times.
✓ Branch 29 taken 69291 times.
✓ Branch 30 taken 30500 times.
✓ Branch 31 taken 71604 times.
✓ Branch 32 taken 3065 times.
✓ Branch 33 taken 5142 times.
✓ Branch 34 taken 1745 times.
✓ Branch 35 taken 4305 times.
✓ Branch 36 taken 1570 times.
✓ Branch 37 taken 3730 times.
✓ Branch 38 taken 625 times.
✓ Branch 39 taken 1532 times.
✓ Branch 40 taken 125 times.
✓ Branch 41 taken 637 times.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✓ Branch 44 taken 1398 times.
✓ Branch 45 taken 2111 times.
✓ Branch 46 taken 18369 times.
✓ Branch 47 taken 69559 times.
|
1242852 | return (notEqual(a[0], b[0]) ? a[0] < b[0] |
| 1327 |
39/48✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 17208 times.
✓ Branch 5 taken 24363 times.
✓ Branch 6 taken 1994 times.
✓ Branch 7 taken 2971 times.
✓ Branch 8 taken 2902 times.
✓ Branch 9 taken 7533 times.
✓ Branch 10 taken 238 times.
✓ Branch 11 taken 1127 times.
✓ Branch 12 taken 93 times.
✓ Branch 13 taken 1156 times.
✓ Branch 14 taken 61 times.
✓ Branch 15 taken 135 times.
✓ Branch 16 taken 64 times.
✓ Branch 17 taken 143 times.
✓ Branch 18 taken 44 times.
✓ Branch 19 taken 176 times.
✓ Branch 20 taken 4 times.
✓ Branch 21 taken 76 times.
✓ Branch 22 taken 5680 times.
✓ Branch 23 taken 4667 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 60960 times.
✓ Branch 27 taken 163732 times.
✓ Branch 28 taken 16279 times.
✓ Branch 29 taken 53012 times.
✓ Branch 30 taken 13436 times.
✓ Branch 31 taken 58168 times.
✓ Branch 32 taken 2140 times.
✓ Branch 33 taken 3002 times.
✓ Branch 34 taken 520 times.
✓ Branch 35 taken 3785 times.
✓ Branch 36 taken 1442 times.
✓ Branch 37 taken 2288 times.
✓ Branch 38 taken 125 times.
✓ Branch 39 taken 1407 times.
✗ Branch 40 not taken.
✓ Branch 41 taken 637 times.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✓ Branch 44 taken 406 times.
✓ Branch 45 taken 1705 times.
✓ Branch 46 taken 17014 times.
✓ Branch 47 taken 52545 times.
|
523238 | : (notEqual(a[1], b[1]) ? a[1] < b[1] |
| 1328 | 382628 | : (a[2] < b[2]))); | |
| 1329 | }); | ||
| 1330 | |||
| 1331 | 14576 | const auto squaredEps = eps*eps; | |
| 1332 |
2/2✓ Branch 2 taken 2660 times.
✓ Branch 3 taken 5633 times.
|
14576 | points.erase(std::unique( |
| 1333 | points.begin(), points.end(), | ||
| 1334 | 319572 | [squaredEps] (const auto& a, const auto&b) { return (b-a).two_norm2() < squaredEps; }), | |
| 1335 | 14576 | points.end() | |
| 1336 | ); | ||
| 1337 | |||
| 1338 | // return false if we don't have more than four unique points (dim+1) | ||
| 1339 |
2/2✓ Branch 0 taken 2660 times.
✓ Branch 1 taken 5633 times.
|
14576 | if (points.size() < 4) |
| 1340 | return false; | ||
| 1341 | |||
| 1342 | // check if the points are coplanar (then we don't have a 3-dimensional intersection) | ||
| 1343 | 24161 | const bool coplanar = [&] | |
| 1344 | { | ||
| 1345 | 5633 | const auto p0 = points[0]; | |
| 1346 | 16899 | const auto normal = crossProduct(points[1] - p0, points[2] - p0); | |
| 1347 | // include the normal in eps (instead of norm*norm) since the normal can become very small | ||
| 1348 | // if the first three points are very close together | ||
| 1349 | 5633 | const auto epsCoplanar = normal.two_norm()*norm*eps_; | |
| 1350 |
4/4✓ Branch 0 taken 4154 times.
✓ Branch 1 taken 23 times.
✓ Branch 2 taken 4528 times.
✓ Branch 3 taken 2400 times.
|
11105 | for (int i = 3; i < points.size(); ++i) |
| 1351 | { | ||
| 1352 | 17364 | const auto ad = points[i] - p0; | |
| 1353 | using std::abs; | ||
| 1354 |
4/4✓ Branch 0 taken 2146 times.
✓ Branch 1 taken 2008 times.
✓ Branch 2 taken 1064 times.
✓ Branch 3 taken 3464 times.
|
8682 | if (abs(normal*ad) > epsCoplanar) |
| 1355 | 3210 | return false; | |
| 1356 | } | ||
| 1357 | |||
| 1358 | return true; | ||
| 1359 | 9264 | }(); | |
| 1360 | |||
| 1361 |
2/2✓ Branch 0 taken 2423 times.
✓ Branch 1 taken 3210 times.
|
9264 | if (coplanar) |
| 1362 | return false; | ||
| 1363 | |||
| 1364 | // we return the intersection as point cloud (that can be triangulated later) | ||
| 1365 |
1/2✓ Branch 1 taken 3210 times.
✗ Branch 2 not taken.
|
4426 | intersection = points; |
| 1366 | |||
| 1367 | return true; | ||
| 1368 | 14668 | } | |
| 1369 | |||
| 1370 | /*! | ||
| 1371 | * \brief Colliding segment and convex polyhedron | ||
| 1372 | * \param geo1/geo2 The geometries to intersect | ||
| 1373 | * \param intersection Container to store the intersection result | ||
| 1374 | * \todo implement overloads for polygon-, segment- or point-like intersections | ||
| 1375 | */ | ||
| 1376 | template<class P = Policy, std::enable_if_t<P::dimIntersection != 3, int> = 0> | ||
| 1377 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | ||
| 1378 | { | ||
| 1379 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 1380 | DUNE_THROW(Dune::NotImplemented, "Polyhedron-polygon intersection detection only implemented for polygon-like intersections"); | ||
| 1381 | } | ||
| 1382 | }; | ||
| 1383 | |||
| 1384 | /*! | ||
| 1385 | * \ingroup Geometry | ||
| 1386 | * \brief A class for polygon--segment intersection in 3d space | ||
| 1387 | */ | ||
| 1388 | template <class Geometry1, class Geometry2, class Policy> | ||
| 1389 | class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 2, 1> | ||
| 1390 | { | ||
| 1391 | enum { dimworld = 3 }; | ||
| 1392 | enum { dim1 = 2 }; | ||
| 1393 | enum { dim2 = 1 }; | ||
| 1394 | |||
| 1395 | public: | ||
| 1396 | using ctype = typename Policy::ctype; | ||
| 1397 | using Point = typename Policy::Point; | ||
| 1398 | using Intersection = typename Policy::Intersection; | ||
| 1399 | |||
| 1400 | private: | ||
| 1401 | static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons | ||
| 1402 | |||
| 1403 | public: | ||
| 1404 | /*! | ||
| 1405 | * \brief Colliding segment and convex polyhedron | ||
| 1406 | * \note Algorithm based on the one from "Real-Time Collision Detection" by Christer Ericson, | ||
| 1407 | * published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc. (Chapter 5.3.6) | ||
| 1408 | * \param geo1/geo2 The geometries to intersect | ||
| 1409 | * \param intersection If the geometries collide, is holds the corner points of | ||
| 1410 | * the intersection object in global coordinates. | ||
| 1411 | * \note This overload is used when point-like intersections are seeked | ||
| 1412 | */ | ||
| 1413 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0> | ||
| 1414 | 387377 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | |
| 1415 | { | ||
| 1416 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 1417 | |||
| 1418 | ctype tfirst, tlast; | ||
| 1419 |
2/2✓ Branch 1 taken 83418 times.
✓ Branch 2 taken 303898 times.
|
387377 | if (intersect_(geo1, geo2, tfirst, tlast)) |
| 1420 | { | ||
| 1421 | // the intersection exists. Export the intersections geometry now: | ||
| 1422 | // s(t) = a + t(b-a) in [tfirst, tlast] | ||
| 1423 | 166882 | intersection = Intersection({geo2.global(tfirst), geo2.global(tlast)}); | |
| 1424 | 83461 | return true; | |
| 1425 | } | ||
| 1426 | |||
| 1427 | return false; | ||
| 1428 | } | ||
| 1429 | |||
| 1430 | /*! | ||
| 1431 | * \brief Colliding segment and convex polyhedron | ||
| 1432 | * \note Algorithm based on the one from "Real-Time Collision Detection" by Christer Ericson, | ||
| 1433 | * published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc. (Chapter 5.3.6) | ||
| 1434 | * \param geo1/geo2 The geometries to intersect | ||
| 1435 | * \param is If the geometries collide, is holds the corner points of | ||
| 1436 | * the intersection object in global coordinates. | ||
| 1437 | * \note This overload is used when point-like intersections are seeked | ||
| 1438 | */ | ||
| 1439 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0> | ||
| 1440 | 3282916 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& is) | |
| 1441 | { | ||
| 1442 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 1443 | |||
| 1444 | 3282916 | const auto p = geo2.corner(0); | |
| 1445 | 3282916 | const auto q = geo2.corner(1); | |
| 1446 | |||
| 1447 | 3282916 | const auto a = geo1.corner(0); | |
| 1448 | 3282916 | const auto b = geo1.corner(1); | |
| 1449 | 3282916 | const auto c = geo1.corner(2); | |
| 1450 | |||
| 1451 |
2/2✓ Branch 0 taken 464047 times.
✓ Branch 1 taken 1671835 times.
|
3244324 | if (geo1.corners() == 3) |
| 1452 | 783135 | return intersect_<Policy>(a, b, c, p, q, is); | |
| 1453 | |||
| 1454 |
1/2✓ Branch 0 taken 1671835 times.
✗ Branch 1 not taken.
|
2461189 | else if (geo1.corners() == 4) |
| 1455 | { | ||
| 1456 | 2499781 | Intersection is1, is2; | |
| 1457 | bool hasSegment1, hasSegment2; | ||
| 1458 | |||
| 1459 | 2499781 | const auto d = geo1.corner(3); | |
| 1460 |
1/2✓ Branch 1 taken 19296 times.
✗ Branch 2 not taken.
|
2499781 | const bool intersects1 = intersect_<Policy>(a, b, d, p, q, is1, hasSegment1); |
| 1461 |
1/2✓ Branch 1 taken 19296 times.
✗ Branch 2 not taken.
|
2499781 | const bool intersects2 = intersect_<Policy>(a, d, c, p, q, is2, hasSegment2); |
| 1462 | |||
| 1463 |
4/4✓ Branch 0 taken 1658191 times.
✓ Branch 1 taken 32940 times.
✓ Branch 2 taken 1626008 times.
✓ Branch 3 taken 32183 times.
|
2499781 | if (hasSegment1 || hasSegment2) |
| 1464 | return false; | ||
| 1465 | |||
| 1466 |
2/2✓ Branch 0 taken 102473 times.
✓ Branch 1 taken 1523535 times.
|
2420704 | if (intersects1) { is = is1; return true; } |
| 1467 |
2/2✓ Branch 0 taken 52572 times.
✓ Branch 1 taken 1470963 times.
|
2296417 | if (intersects2) { is = is2; return true; } |
| 1468 | |||
| 1469 | return false; | ||
| 1470 | } | ||
| 1471 | |||
| 1472 | else | ||
| 1473 | ✗ | DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type " | |
| 1474 | << geo1.type() << ", "<< geo1.corners() << " corners."); | ||
| 1475 | } | ||
| 1476 | |||
| 1477 | private: | ||
| 1478 | /*! | ||
| 1479 | * \brief Obtain local coordinates of start/end point of the intersecting segment. | ||
| 1480 | */ | ||
| 1481 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0> | ||
| 1482 | 387377 | static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast) | |
| 1483 | { | ||
| 1484 | // lambda to obtain the facet corners on geo1 | ||
| 1485 | 387316 | auto getFacetCorners = [] (const Geometry1& geo1) | |
| 1486 | { | ||
| 1487 | 387316 | std::vector< std::array<int, 2> > facetCorners; | |
| 1488 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 387292 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13 times.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
|
387316 | switch (geo1.corners()) |
| 1489 | { | ||
| 1490 | 13 | case 4: // quadrilateral | |
| 1491 |
1/4✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
|
13 | facetCorners = {{0, 2}, {3, 1}, {1, 0}, {2, 3}}; |
| 1492 | 13 | break; | |
| 1493 | 387303 | case 3: // triangle | |
| 1494 |
2/8✓ Branch 1 taken 387292 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 6 taken 11 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
387303 | facetCorners = {{1, 0}, {0, 2}, {2, 1}}; |
| 1495 | 387303 | break; | |
| 1496 | ✗ | default: | |
| 1497 | ✗ | DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type " | |
| 1498 | << geo1.type() << " with "<< geo1.corners() << " corners."); | ||
| 1499 | } | ||
| 1500 | |||
| 1501 | 387316 | return facetCorners; | |
| 1502 | ✗ | }; | |
| 1503 | |||
| 1504 | 387377 | const auto center1 = geo1.center(); | |
| 1505 | 387425 | const auto normal1 = crossProduct(geo1.corner(1) - geo1.corner(0), | |
| 1506 | 774754 | geo1.corner(2) - geo1.corner(0)); | |
| 1507 | |||
| 1508 | // lambda to obtain the normal on a facet | ||
| 1509 | 1204754 | auto computeNormal = [¢er1, &normal1, &geo1] (const std::array<int, 2>& facetCorners) | |
| 1510 | { | ||
| 1511 | 817283 | const auto c0 = geo1.corner(facetCorners[0]); | |
| 1512 | 817377 | const auto c1 = geo1.corner(facetCorners[1]); | |
| 1513 | 817377 | const auto edge = c1 - c0; | |
| 1514 | |||
| 1515 | 817377 | Dune::FieldVector<ctype, dimworld> n = crossProduct(edge, normal1); | |
| 1516 | 817377 | n /= n.two_norm(); | |
| 1517 | |||
| 1518 | // orientation of the element is not known | ||
| 1519 | // make sure normal points outwards of element | ||
| 1520 |
2/4✓ Branch 0 taken 817293 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 84 times.
✗ Branch 3 not taken.
|
2452131 | if ( n*(center1-c0) > 0.0 ) |
| 1521 | 817377 | n *= -1.0; | |
| 1522 | |||
| 1523 | 817377 | return n; | |
| 1524 | }; | ||
| 1525 | |||
| 1526 | 387377 | return Detail::computeSegmentIntersection(geo1, geo2, eps_, tfirst, tlast, getFacetCorners, computeNormal); | |
| 1527 | } | ||
| 1528 | |||
| 1529 | /*! | ||
| 1530 | * \brief triangle--segment point-like intersection with points as input. | ||
| 1531 | */ | ||
| 1532 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0> | ||
| 1533 | 464047 | static bool intersect_(const Point& a, const Point& b, const Point& c, | |
| 1534 | const Point& p, const Point& q, | ||
| 1535 | Intersection& is) | ||
| 1536 | { | ||
| 1537 | bool hasSegment; | ||
| 1538 | 464047 | return intersect_(a, b, c, p, q, is, hasSegment); | |
| 1539 | } | ||
| 1540 | |||
| 1541 | /*! | ||
| 1542 | * \brief triangle--segment point-like intersection with points as input. Also | ||
| 1543 | * stores if a segment-like intersection was found in the provided boolean. | ||
| 1544 | */ | ||
| 1545 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0> | ||
| 1546 | 3846309 | static bool intersect_(const Point& a, const Point& b, const Point& c, | |
| 1547 | const Point& p, const Point& q, | ||
| 1548 | Intersection& is, bool& hasSegmentIs) | ||
| 1549 | { | ||
| 1550 | 3846309 | hasSegmentIs = false; | |
| 1551 | |||
| 1552 | 3846309 | const auto ab = b - a; | |
| 1553 | 3846309 | const auto ac = c - a; | |
| 1554 | 3846309 | const auto qp = p - q; | |
| 1555 | |||
| 1556 | // compute the triangle normal that defines the triangle plane | ||
| 1557 | 3846309 | const auto normal = crossProduct(ab, ac); | |
| 1558 | |||
| 1559 | // compute the denominator | ||
| 1560 | // if denom is 0 the segment is parallel and we can return | ||
| 1561 |
2/2✓ Branch 0 taken 5729763 times.
✓ Branch 1 taken 1909921 times.
|
15385236 | const auto denom = normal*qp; |
| 1562 | 3846309 | const auto abnorm2 = ab.two_norm2(); | |
| 1563 |
2/2✓ Branch 0 taken 1157013 times.
✓ Branch 1 taken 752908 times.
|
7692618 | const auto eps = eps_*abnorm2*qp.two_norm(); |
| 1564 | using std::abs; | ||
| 1565 |
2/2✓ Branch 0 taken 1157013 times.
✓ Branch 1 taken 752908 times.
|
3846309 | if (abs(denom) < eps) |
| 1566 | { | ||
| 1567 | 1758853 | const auto pa = a - p; | |
| 1568 | 1758853 | const auto denom2 = normal*pa; | |
| 1569 |
2/2✓ Branch 0 taken 864893 times.
✓ Branch 1 taken 292120 times.
|
3517706 | if (abs(denom2) > eps_*pa.two_norm()*abnorm2) |
| 1570 | return false; | ||
| 1571 | |||
| 1572 | // if we get here, we are in-plane. Check if a | ||
| 1573 | // segment-like intersection with geometry 1 exists. | ||
| 1574 | using SegmentPolicy = typename IntersectionPolicy::SegmentPolicy<ctype, dimworld>; | ||
| 1575 | using Triangle = Dune::AffineGeometry<ctype, 2, dimworld>; | ||
| 1576 | using Segment = Dune::AffineGeometry<ctype, 1, dimworld>; | ||
| 1577 | using SegmentIntersectionAlgorithm = GeometryIntersection<Triangle, Segment, SegmentPolicy>; | ||
| 1578 | using SegmentIntersectionType = typename SegmentIntersectionAlgorithm::Intersection; | ||
| 1579 | 387292 | SegmentIntersectionType segmentIs; | |
| 1580 | |||
| 1581 |
1/2✓ Branch 1 taken 292120 times.
✗ Branch 2 not taken.
|
387292 | Triangle triangle(Dune::GeometryTypes::simplex(2), std::array<Point, 3>({a, b, c})); |
| 1582 |
1/2✓ Branch 1 taken 292120 times.
✗ Branch 2 not taken.
|
387292 | Segment segment(Dune::GeometryTypes::simplex(1), std::array<Point, 2>({p, q})); |
| 1583 |
3/4✓ Branch 1 taken 292120 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 58012 times.
✓ Branch 4 taken 234108 times.
|
387292 | if (SegmentIntersectionAlgorithm::intersection(triangle, segment, segmentIs)) |
| 1584 | { | ||
| 1585 | 83398 | hasSegmentIs = true; | |
| 1586 | 83398 | return false; | |
| 1587 | } | ||
| 1588 | |||
| 1589 | // We are in-plane. A point-like | ||
| 1590 | // intersection can only be on the edges | ||
| 1591 |
2/2✓ Branch 0 taken 424188 times.
✓ Branch 1 taken 146602 times.
|
778354 | for (const auto& ip : {p, q}) |
| 1592 | { | ||
| 1593 | 563146 | if ( intersectsPointSimplex(ip, a, b) | |
| 1594 |
2/2✓ Branch 1 taken 344380 times.
✓ Branch 2 taken 25067 times.
|
507945 | || intersectsPointSimplex(ip, b, c) |
| 1595 |
4/4✓ Branch 0 taken 369447 times.
✓ Branch 1 taken 54741 times.
✓ Branch 3 taken 7698 times.
✓ Branch 4 taken 336682 times.
|
1045944 | || intersectsPointSimplex(ip, c, a) ) |
| 1596 | { | ||
| 1597 | 88686 | is = ip; | |
| 1598 | 88686 | return true; | |
| 1599 | } | ||
| 1600 | } | ||
| 1601 | |||
| 1602 | 215208 | return false; | |
| 1603 | } | ||
| 1604 | |||
| 1605 | // compute intersection t value of pq with plane of triangle. | ||
| 1606 | // a segment intersects if and only if 0 <= t <= 1. | ||
| 1607 | 2087456 | const auto ap = p - a; | |
| 1608 | 2087456 | const auto t = (ap*normal)/denom; | |
| 1609 |
2/2✓ Branch 0 taken 155229 times.
✓ Branch 1 taken 597679 times.
|
2087456 | if (t < 0.0 - eps_) return false; |
| 1610 |
2/2✓ Branch 0 taken 148182 times.
✓ Branch 1 taken 449497 times.
|
1589635 | if (t > 1.0 + eps_) return false; |
| 1611 | |||
| 1612 | // compute the barycentric coordinates and check if the intersection point | ||
| 1613 | // is inside the bounds of the triangle | ||
| 1614 | 1102240 | const auto e = crossProduct(qp, ap); | |
| 1615 | 1102240 | const auto v = (ac*e)/denom; | |
| 1616 |
4/4✓ Branch 0 taken 116047 times.
✓ Branch 1 taken 333450 times.
✓ Branch 2 taken 71437 times.
✓ Branch 3 taken 262013 times.
|
1102240 | if (v < -eps_ || v > 1.0 + eps_) return false; |
| 1617 | 448814 | const auto w = -(ab*e)/denom; | |
| 1618 |
4/4✓ Branch 0 taken 60319 times.
✓ Branch 1 taken 201694 times.
✓ Branch 2 taken 35874 times.
✓ Branch 3 taken 165820 times.
|
448814 | if (w < -eps_ || v + w > 1.0 + eps_) return false; |
| 1619 | |||
| 1620 | // Now we are sure there is an intersection points | ||
| 1621 | // Perform delayed division compute the last barycentric coordinate component | ||
| 1622 | 234923 | const auto u = 1.0 - v - w; | |
| 1623 | |||
| 1624 | 234923 | Point ip(0.0); | |
| 1625 |
2/2✓ Branch 0 taken 497460 times.
✓ Branch 1 taken 165820 times.
|
939692 | ip.axpy(u, a); |
| 1626 |
2/2✓ Branch 0 taken 497460 times.
✓ Branch 1 taken 165820 times.
|
939692 | ip.axpy(v, b); |
| 1627 | 234923 | ip.axpy(w, c); | |
| 1628 | 234923 | is = ip; | |
| 1629 | 234923 | return true; | |
| 1630 | } | ||
| 1631 | }; | ||
| 1632 | |||
| 1633 | /*! | ||
| 1634 | * \ingroup Geometry | ||
| 1635 | * \brief A class for segment--polygon intersection in 3d space | ||
| 1636 | */ | ||
| 1637 | template <class Geometry1, class Geometry2, class Policy> | ||
| 1638 | class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 1, 2> | ||
| 1639 | : public GeometryIntersection<Geometry2, Geometry1, Policy, 3, 2, 1> | ||
| 1640 | { | ||
| 1641 | using Base = GeometryIntersection<Geometry2, Geometry1, Policy, 3, 2, 1>; | ||
| 1642 | public: | ||
| 1643 | /*! | ||
| 1644 | * \brief Colliding segment and convex polygon | ||
| 1645 | * \param geo1/geo2 The geometries to intersect | ||
| 1646 | * \param intersection If the geometries collide intersection holds the | ||
| 1647 | * corner points of the intersection object in global coordinates. | ||
| 1648 | * \note This forwards to the polyhedron-polygon specialization with swapped arguments. | ||
| 1649 | */ | ||
| 1650 | template<class P = Policy> | ||
| 1651 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename Base::Intersection& intersection) | ||
| 1652 | { | ||
| 1653 | return Base::intersection(geo2, geo1, intersection); | ||
| 1654 | } | ||
| 1655 | }; | ||
| 1656 | |||
| 1657 | /*! | ||
| 1658 | * \ingroup Geometry | ||
| 1659 | * \brief A class for segment--segment intersection in 3d space | ||
| 1660 | */ | ||
| 1661 | template <class Geometry1, class Geometry2, class Policy> | ||
| 1662 | class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 1, 1> | ||
| 1663 | { | ||
| 1664 | enum { dimworld = 3 }; | ||
| 1665 | enum { dim1 = 1 }; | ||
| 1666 | enum { dim2 = 1 }; | ||
| 1667 | |||
| 1668 | public: | ||
| 1669 | using ctype = typename Policy::ctype; | ||
| 1670 | using Point = typename Policy::Point; | ||
| 1671 | using Intersection = typename Policy::Intersection; | ||
| 1672 | |||
| 1673 | private: | ||
| 1674 | static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons | ||
| 1675 | |||
| 1676 | public: | ||
| 1677 | /*! | ||
| 1678 | * \brief Colliding two segments | ||
| 1679 | * \param geo1/geo2 The geometries to intersect | ||
| 1680 | * \param intersection The intersection point | ||
| 1681 | * \note This overload is used when point-like intersections are seeked | ||
| 1682 | */ | ||
| 1683 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0> | ||
| 1684 | 418 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | |
| 1685 | { | ||
| 1686 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 1687 | |||
| 1688 | 418 | const auto v1 = geo1.corner(1) - geo1.corner(0); | |
| 1689 | 418 | const auto v2 = geo2.corner(1) - geo2.corner(0); | |
| 1690 | 836 | const auto ac = geo2.corner(0) - geo1.corner(0); | |
| 1691 | |||
| 1692 | 418 | const auto v1Norm2 = v1.two_norm2(); | |
| 1693 | 418 | const auto eps2 = eps_*v1Norm2; | |
| 1694 | |||
| 1695 | 418 | const auto n = crossProduct(v1, v2); | |
| 1696 | |||
| 1697 | // first check if segments are parallel | ||
| 1698 | using std::abs; | ||
| 1699 |
2/2✓ Branch 0 taken 100 times.
✓ Branch 1 taken 318 times.
|
418 | if ( n.two_norm2() < eps2*v1Norm2 ) |
| 1700 | { | ||
| 1701 | // check if they lie on the same line | ||
| 1702 |
2/2✓ Branch 1 taken 49 times.
✓ Branch 2 taken 51 times.
|
200 | if (crossProduct(v1, ac).two_norm2() > eps2) |
| 1703 | return false; | ||
| 1704 | |||
| 1705 | // they lie on the same line, | ||
| 1706 | // if so, point intersection can only be one of the corners | ||
| 1707 | 51 | const auto sp = v1*v2; | |
| 1708 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 42 times.
|
51 | if ( sp < 0.0 ) |
| 1709 | { | ||
| 1710 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 5 times.
|
9 | if ( ac.two_norm2() < eps2 ) |
| 1711 | 4 | { intersection = geo2.corner(0); return true; } | |
| 1712 | |||
| 1713 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 5 times.
|
15 | if ( (geo2.corner(1) - geo1.corner(1)).two_norm2() < eps2 ) |
| 1714 | ✗ | { intersection = geo2.corner(1); return true; } | |
| 1715 | } | ||
| 1716 | else | ||
| 1717 | { | ||
| 1718 |
2/2✓ Branch 2 taken 4 times.
✓ Branch 3 taken 38 times.
|
126 | if ( (geo2.corner(1) - geo1.corner(0)).two_norm2() < eps2 ) |
| 1719 | 4 | { intersection = geo2.corner(1); return true; } | |
| 1720 | |||
| 1721 |
2/2✓ Branch 2 taken 8 times.
✓ Branch 3 taken 30 times.
|
114 | if ( (geo2.corner(0) - geo1.corner(1)).two_norm2() < eps2 ) |
| 1722 | 8 | { intersection = geo2.corner(0); return true; } | |
| 1723 | } | ||
| 1724 | |||
| 1725 | // no intersection | ||
| 1726 | 35 | return false; | |
| 1727 | } | ||
| 1728 | |||
| 1729 | // in-plane normal vector on v1 | ||
| 1730 | 318 | const auto v1Normal = crossProduct(v1, n); | |
| 1731 | |||
| 1732 | // intersection point on v2 in local coords | ||
| 1733 | 318 | const auto t2 = - 1.0*(ac*v1Normal) / (v2*v1Normal); | |
| 1734 | |||
| 1735 | // check if the local coords are valid | ||
| 1736 |
4/4✓ Branch 0 taken 64 times.
✓ Branch 1 taken 254 times.
✓ Branch 2 taken 40 times.
✓ Branch 3 taken 214 times.
|
318 | if (t2 < -1.0*eps_ || t2 > 1.0 + eps_) |
| 1737 | return false; | ||
| 1738 | |||
| 1739 |
2/2✓ Branch 2 taken 162 times.
✓ Branch 3 taken 52 times.
|
214 | if (auto ip = geo2.global(t2); intersectsPointGeometry(ip, geo1)) |
| 1740 | 162 | { intersection = std::move(ip); return true; } | |
| 1741 | |||
| 1742 | 52 | return false; | |
| 1743 | } | ||
| 1744 | |||
| 1745 | /*! | ||
| 1746 | * \brief Colliding two segments in 3D | ||
| 1747 | * \param geo1/geo2 The geometries to intersect | ||
| 1748 | * \param intersection If the geometries collide, is holds the corner points of | ||
| 1749 | * the intersection object in global coordinates. | ||
| 1750 | * \note This overload is used when segment-like intersections are seeked | ||
| 1751 | */ | ||
| 1752 | template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0> | ||
| 1753 | 478 | static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) | |
| 1754 | { | ||
| 1755 | static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); | ||
| 1756 | |||
| 1757 | 478 | const auto& a = geo1.corner(0); | |
| 1758 | 478 | const auto& b = geo1.corner(1); | |
| 1759 | 478 | const auto ab = b-a; | |
| 1760 | |||
| 1761 | 478 | const auto& p = geo2.corner(0); | |
| 1762 | 478 | const auto& q = geo2.corner(1); | |
| 1763 | 478 | const auto pq = q-p; | |
| 1764 | |||
| 1765 | 478 | const auto abNorm2 = ab.two_norm2(); | |
| 1766 |
2/2✓ Branch 0 taken 35 times.
✓ Branch 1 taken 443 times.
|
478 | const auto pqNorm2 = pq.two_norm2(); |
| 1767 | |||
| 1768 | using std::max; | ||
| 1769 | 478 | const auto eps2 = eps_*max(abNorm2, pqNorm2); | |
| 1770 | |||
| 1771 | // if the segment are not parallel there is no segment intersection | ||
| 1772 | using std::abs; | ||
| 1773 |
2/2✓ Branch 1 taken 4 times.
✓ Branch 2 taken 474 times.
|
956 | if (crossProduct(ab, pq).two_norm2() > eps2*eps2) |
| 1774 | return false; | ||
| 1775 | |||
| 1776 | 474 | const auto ap = (p-a); | |
| 1777 | 474 | const auto aq = (q-a); | |
| 1778 | |||
| 1779 | // points have to be colinear | ||
| 1780 |
2/2✓ Branch 1 taken 4 times.
✓ Branch 2 taken 470 times.
|
948 | if (crossProduct(ap, aq).two_norm2() > eps2*eps2) |
| 1781 | return false; | ||
| 1782 | |||
| 1783 | // scaled local coordinates | ||
| 1784 | // we save the division for the normalization for now | ||
| 1785 | // and do it in the very end if we find an intersection | ||
| 1786 | 470 | auto tp = ap*ab; | |
| 1787 | 470 | auto tq = aq*ab; | |
| 1788 | |||
| 1789 | // make sure they are sorted | ||
| 1790 | using std::swap; | ||
| 1791 |
2/2✓ Branch 0 taken 33 times.
✓ Branch 1 taken 437 times.
|
470 | if (tp > tq) |
| 1792 | 33 | swap(tp, tq); | |
| 1793 | |||
| 1794 | using std::clamp; | ||
| 1795 |
2/2✓ Branch 0 taken 136 times.
✓ Branch 1 taken 334 times.
|
470 | tp = clamp(tp, 0.0, abNorm2); |
| 1796 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 466 times.
|
470 | tq = clamp(tq, 0.0, abNorm2); |
| 1797 | |||
| 1798 |
2/2✓ Branch 0 taken 132 times.
✓ Branch 1 taken 338 times.
|
470 | if (abs(tp-tq) < eps2) |
| 1799 | return false; | ||
| 1800 | |||
| 1801 | 338 | intersection = Intersection({geo1.global(tp/abNorm2), geo1.global(tq/abNorm2)}); | |
| 1802 | 338 | return true; | |
| 1803 | } | ||
| 1804 | }; | ||
| 1805 | |||
| 1806 | } // end namespace Dumux | ||
| 1807 | |||
| 1808 | #endif | ||
| 1809 |