| 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-FileCopyrightInfo: 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 | 754373 | 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 | 754373 | const auto a = geo2.corner(0); | |
| 135 | 754373 | const auto b = geo2.corner(1); | |
| 136 | 754373 | 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 | 754373 | tfirst = 0.0; | |
| 142 | 754373 | tlast = 1.0; | |
| 143 | |||
| 144 |
1/2✓ Branch 1 taken 359140 times.
✗ Branch 2 not taken.
|
1508746 | const auto facets = getFacetCornerIndices(geo1); |
| 145 |
4/4✓ Branch 0 taken 1577669 times.
✓ Branch 1 taken 124140 times.
✓ Branch 2 taken 1577669 times.
✓ Branch 3 taken 124140 times.
|
3860328 | for (const auto& f : facets) |
| 146 | { | ||
| 147 |
1/2✓ Branch 1 taken 25123 times.
✗ Branch 2 not taken.
|
2215647 | const auto n = computeNormal(f); |
| 148 | |||
| 149 |
1/2✓ Branch 1 taken 25123 times.
✗ Branch 2 not taken.
|
3092591 | const auto c0 = geo1.corner(f[0]); |
| 150 | 2215647 | const ctype denom = n*d; | |
| 151 | 2232652 | 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.
|
7974719 | const auto edge1 = geo1.corner(f[1]) - geo1.corner(f[0]); |
| 155 | 2215647 | 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 |
4/4✓ Branch 0 taken 382776 times.
✓ Branch 1 taken 1194893 times.
✓ Branch 2 taken 382776 times.
✓ Branch 3 taken 1194893 times.
|
4431294 | if (abs(denom) < eps) |
| 161 | { | ||
| 162 |
2/2✓ Branch 0 taken 320704 times.
✓ Branch 1 taken 62072 times.
|
422989 | if (dist > eps) |
| 163 | 618438 | return false; | |
| 164 | } | ||
| 165 | else // not parallel: compute line-line intersection | ||
| 166 | { | ||
| 167 | 1792658 | const ctype t = -dist / denom; | |
| 168 | // if entering half space cut tfirst if t is larger | ||
| 169 | using std::signbit; | ||
| 170 |
4/4✓ Branch 0 taken 671793 times.
✓ Branch 1 taken 523100 times.
✓ Branch 2 taken 671793 times.
✓ Branch 3 taken 523100 times.
|
3585316 | if (signbit(denom)) |
| 171 | { | ||
| 172 |
2/2✓ Branch 0 taken 395915 times.
✓ Branch 1 taken 275878 times.
|
1012979 | if (t > tfirst) |
| 173 | 637726 | tfirst = t; | |
| 174 | } | ||
| 175 | // if exiting half space cut tlast if t is smaller | ||
| 176 | else | ||
| 177 | { | ||
| 178 |
2/2✓ Branch 0 taken 371542 times.
✓ Branch 1 taken 151558 times.
|
779679 | if (t < tlast) |
| 179 | 605347 | 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 821553 times.
✓ Branch 1 taken 373340 times.
|
1792658 | if (tfirst > tlast - baseEps) |
| 184 | return false; | ||
| 185 | } | ||
| 186 | } | ||
| 187 | |||
| 188 | // If we made it until here an intersection exists. | ||
| 189 | return true; | ||
| 190 | } | ||
| 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 | 1896112 | 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 | 1896112 | const auto v1 = geo1.corner(1) - geo1.corner(0); | |
| 255 | 1896112 | const auto v2 = geo2.corner(1) - geo2.corner(0); | |
| 256 | 1896112 | const auto ac = geo2.corner(0) - geo1.corner(0); | |
| 257 | |||
| 258 | 5688336 | auto n2 = Point({-1.0*v2[1], v2[0]}); | |
| 259 | 3792224 | n2 /= n2.two_norm(); | |
| 260 | |||
| 261 | // compute distance of first corner on geo2 to line1 | ||
| 262 | const auto dist12 = n2*ac; | ||
| 263 | |||
| 264 | // first check parallel segments | ||
| 265 | using std::abs; | ||
| 266 | using std::sqrt; | ||
| 267 | |||
| 268 | 1896112 | const auto v1Norm2 = v1.two_norm2(); | |
| 269 | 1896112 | const auto eps = eps_*sqrt(v1Norm2); | |
| 270 | 1896112 | const auto eps2 = eps_*v1Norm2; | |
| 271 | |||
| 272 | 1896112 | const auto sp = n2*v1; | |
| 273 |
4/4✓ Branch 0 taken 947872 times.
✓ Branch 1 taken 948240 times.
✓ Branch 2 taken 947872 times.
✓ Branch 3 taken 948240 times.
|
3792224 | if (abs(sp) < eps) |
| 274 | { | ||
| 275 |
4/4✓ Branch 0 taken 213648 times.
✓ Branch 1 taken 734224 times.
✓ Branch 2 taken 213648 times.
✓ Branch 3 taken 734224 times.
|
1895744 | 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 213640 times.
|
213648 | 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 63892 times.
✓ Branch 3 taken 149748 times.
|
640920 | if ( (geo2.corner(1) - geo1.corner(0)).two_norm2() < eps2 ) |
| 290 | 63892 | { intersection = geo2.corner(1); return true; } | |
| 291 | |||
| 292 |
2/2✓ Branch 2 taken 63896 times.
✓ Branch 3 taken 85852 times.
|
449244 | if ( (geo2.corner(0) - geo1.corner(1)).two_norm2() < eps2 ) |
| 293 | 63896 | { intersection = geo2.corner(0); return true; } | |
| 294 | } | ||
| 295 | |||
| 296 | // no intersection | ||
| 297 | 85856 | return false; | |
| 298 | } | ||
| 299 | |||
| 300 | // intersection point on v1 in local coords | ||
| 301 | 948240 | const auto t1 = dist12 / sp; | |
| 302 | |||
| 303 | // check if the local coords are valid | ||
| 304 |
4/4✓ Branch 0 taken 796348 times.
✓ Branch 1 taken 151892 times.
✓ Branch 2 taken 644404 times.
✓ Branch 3 taken 151944 times.
|
948240 | if (t1 < -1.0*eps_ || t1 > 1.0 + eps_) |
| 305 | return false; | ||
| 306 | |||
| 307 | // compute global coordinates | ||
| 308 | 1288808 | auto isPoint = geo1.global(t1); | |
| 309 | |||
| 310 | // check if point is in bounding box of v2 | ||
| 311 | 644404 | const auto c = geo2.corner(0); | |
| 312 | 644404 | const auto d = geo2.corner(1); | |
| 313 | |||
| 314 | using std::min; using std::max; | ||
| 315 |
24/24✓ Branch 0 taken 60 times.
✓ Branch 1 taken 644344 times.
✓ Branch 2 taken 60 times.
✓ Branch 3 taken 644344 times.
✓ Branch 4 taken 60 times.
✓ Branch 5 taken 644344 times.
✓ Branch 6 taken 8 times.
✓ Branch 7 taken 644396 times.
✓ Branch 8 taken 8 times.
✓ Branch 9 taken 644396 times.
✓ Branch 10 taken 8 times.
✓ Branch 11 taken 644396 times.
✓ Branch 12 taken 322208 times.
✓ Branch 13 taken 322196 times.
✓ Branch 14 taken 322208 times.
✓ Branch 15 taken 322196 times.
✓ Branch 16 taken 322208 times.
✓ Branch 17 taken 322196 times.
✓ Branch 18 taken 322236 times.
✓ Branch 19 taken 322168 times.
✓ Branch 20 taken 322236 times.
✓ Branch 21 taken 322168 times.
✓ Branch 22 taken 322236 times.
✓ Branch 23 taken 322168 times.
|
2900040 | 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 |
4/4✓ Branch 0 taken 302220 times.
✓ Branch 1 taken 342184 times.
✓ Branch 2 taken 302220 times.
✓ Branch 3 taken 342184 times.
|
1288808 | if ( intersectsPointBoundingBox(isPoint, bBox.data()) ) |
| 318 | { | ||
| 319 | 302220 | intersection = std::move(isPoint); | |
| 320 | 302220 | 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 | 68 | 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 | 68 | const auto& a = geo1.corner(0); | |
| 338 | 68 | const auto& b = geo1.corner(1); | |
| 339 | 68 | const auto ab = b - a; | |
| 340 | |||
| 341 | 68 | const auto& p = geo2.corner(0); | |
| 342 | 68 | const auto& q = geo2.corner(1); | |
| 343 | 68 | const auto pq = q - p; | |
| 344 | 204 | Dune::FieldVector<ctype, dimworld> n2{-pq[1], pq[0]}; | |
| 345 | |||
| 346 | using std::max; | ||
| 347 | 68 | const auto abNorm2 = ab.two_norm2(); | |
| 348 | 68 | const auto pqNorm2 = pq.two_norm2(); | |
| 349 |
2/2✓ Branch 0 taken 20 times.
✓ Branch 1 taken 48 times.
|
68 | const auto eps2 = eps_*max(abNorm2, pqNorm2); |
| 350 | |||
| 351 | // non-parallel segments do not intersect in a segment. | ||
| 352 | using std::abs; | ||
| 353 |
4/4✓ Branch 0 taken 64 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 64 times.
✓ Branch 3 taken 4 times.
|
204 | if (abs(n2*ab) > eps2) |
| 354 | return false; | ||
| 355 | |||
| 356 | // check if the segments lie on the same line | ||
| 357 | 64 | const auto ap = p - a; | |
| 358 |
4/4✓ Branch 0 taken 60 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 60 times.
✓ Branch 3 taken 4 times.
|
192 | if (abs(ap*n2) > eps2) |
| 359 | return false; | ||
| 360 | |||
| 361 | // compute scaled local coordinates of corner 1/2 of segment2 on segment1 | ||
| 362 | 60 | auto t1 = ab*ap; | |
| 363 | 120 | auto t2 = ab*(q - a); | |
| 364 | |||
| 365 | using std::swap; | ||
| 366 |
2/2✓ Branch 0 taken 33 times.
✓ Branch 1 taken 27 times.
|
60 | if (t1 > t2) |
| 367 | 33 | swap(t1, t2); | |
| 368 | |||
| 369 | using std::clamp; | ||
| 370 |
2/2✓ Branch 0 taken 44 times.
✓ Branch 1 taken 16 times.
|
60 | t1 = clamp(t1, 0.0, abNorm2); |
| 371 |
2/2✓ Branch 0 taken 56 times.
✓ Branch 1 taken 4 times.
|
60 | t2 = clamp(t2, 0.0, abNorm2); |
| 372 | |||
| 373 |
4/4✓ Branch 0 taken 28 times.
✓ Branch 1 taken 32 times.
✓ Branch 2 taken 28 times.
✓ Branch 3 taken 32 times.
|
120 | if (abs(t2-t1) < eps2) |
| 374 | return false; | ||
| 375 | |||
| 376 | 56 | intersection = Intersection({geo1.global(t1/abNorm2), geo1.global(t2/abNorm2)}); | |
| 377 | 28 | 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 | 320 | 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/4✓ Branch 1 taken 306 times.
✓ Branch 2 taken 14 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
320 | 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 | 1224 | intersection = Intersection({geo2.global(tfirst), geo2.global(tlast)}); | |
| 419 | 306 | 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 | 8 | 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 | 4800 | static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast) | |
| 457 | { | ||
| 458 | // lambda to obtain the facet corners on geo1 | ||
| 459 | 4752 | auto getFacetCorners = [] (const Geometry1& geo1) | |
| 460 | { | ||
| 461 |
4/10✓ Branch 0 taken 13 times.
✓ Branch 1 taken 4715 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 13 times.
✓ Branch 6 taken 11 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
4752 | std::vector< std::array<int, 2> > facetCorners; |
| 462 |
8/12✓ Branch 0 taken 13 times.
✓ Branch 1 taken 4715 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13 times.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 13 times.
✓ Branch 7 taken 11 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 13 times.
✓ Branch 10 taken 11 times.
✗ Branch 11 not taken.
|
4800 | switch (geo1.corners()) |
| 463 | { | ||
| 464 | 26 | case 4: // quadrilateral | |
| 465 |
2/4✓ Branch 1 taken 4717 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
|
4730 | facetCorners = {{0, 2}, {3, 1}, {1, 0}, {2, 3}}; |
| 466 | break; | ||
| 467 | 22 | case 3: // triangle | |
| 468 |
2/4✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
|
22 | facetCorners = {{1, 0}, {0, 2}, {2, 1}}; |
| 469 | break; | ||
| 470 | ✗ | default: | |
| 471 | ✗ | DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type " | |
| 472 | << geo1.type() << " with "<< geo1.corners() << " corners."); | ||
| 473 | } | ||
| 474 | |||
| 475 | 14256 | return facetCorners; | |
| 476 | }; | ||
| 477 | |||
| 478 | // lambda to obtain the normal on a facet | ||
| 479 | 4800 | const auto center1 = geo1.center(); | |
| 480 | 55311 | auto computeNormal = [&geo1, ¢er1] (const std::array<int, 2>& facetCorners) | |
| 481 | { | ||
| 482 | 33674 | const auto c0 = geo1.corner(facetCorners[0]); | |
| 483 | 33674 | const auto c1 = geo1.corner(facetCorners[1]); | |
| 484 | 16837 | const auto edge = c1 - c0; | |
| 485 | |||
| 486 | 50511 | Dune::FieldVector<ctype, dimworld> n({-1.0*edge[1], edge[0]}); | |
| 487 | 33674 | 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 16753 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 84 times.
|
50511 | if ( n*(center1-c0) > 0.0 ) |
| 492 | n *= -1.0; | ||
| 493 | |||
| 494 | 16837 | return n; | |
| 495 | }; | ||
| 496 | |||
| 497 | 4800 | 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 | 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.
|
4406 | 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 | 196032 | 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 |
3/6✓ Branch 1 taken 195772 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 195772 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 195772 times.
✗ Branch 7 not taken.
|
784128 | std::vector<Point> points; points.reserve(6); |
| 565 | |||
| 566 | // add polygon1 corners that are inside polygon2 | ||
| 567 |
4/4✓ Branch 0 taken 783048 times.
✓ Branch 1 taken 195772 times.
✓ Branch 2 taken 783048 times.
✓ Branch 3 taken 195772 times.
|
1175892 | for (int i = 0; i < geo1.corners(); ++i) |
| 568 |
4/5✓ Branch 1 taken 134324 times.
✓ Branch 2 taken 648724 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 44 times.
✓ Branch 5 taken 76 times.
|
1567736 | if (intersectsPointGeometry(geo1.corner(i), geo2)) |
| 569 |
2/3✓ Branch 1 taken 134324 times.
✓ Branch 2 taken 44 times.
✗ Branch 3 not taken.
|
268736 | points.emplace_back(geo1.corner(i)); |
| 570 | |||
| 571 |
2/2✓ Branch 0 taken 195192 times.
✓ Branch 1 taken 580 times.
|
196032 | const auto numPoints1 = points.size(); |
| 572 |
4/4✓ Branch 0 taken 195192 times.
✓ Branch 1 taken 580 times.
✓ Branch 2 taken 195192 times.
✓ Branch 3 taken 580 times.
|
391868 | if (numPoints1 != geo1.corners()) |
| 573 | { | ||
| 574 | // add polygon2 corners that are inside polygon1 | ||
| 575 |
4/4✓ Branch 0 taken 780744 times.
✓ Branch 1 taken 195192 times.
✓ Branch 2 taken 780744 times.
✓ Branch 3 taken 195192 times.
|
1172688 | for (int i = 0; i < geo2.corners(); ++i) |
| 576 |
4/5✓ Branch 1 taken 510524 times.
✓ Branch 2 taken 270220 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 48 times.
✓ Branch 5 taken 72 times.
|
1563448 | if (intersectsPointGeometry(geo2.corner(i), geo1)) |
| 577 |
2/3✓ Branch 1 taken 510524 times.
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
|
1022496 | points.emplace_back(geo2.corner(i)); |
| 578 | |||
| 579 |
2/4✓ Branch 0 taken 195192 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 195192 times.
✗ Branch 3 not taken.
|
390904 | if (points.empty()) |
| 580 | return false; | ||
| 581 | |||
| 582 |
6/6✓ Branch 0 taken 118392 times.
✓ Branch 1 taken 76800 times.
✓ Branch 2 taken 118392 times.
✓ Branch 3 taken 76800 times.
✓ Branch 4 taken 118392 times.
✓ Branch 5 taken 76800 times.
|
586320 | if (points.size() - numPoints1 != geo2.corners()) |
| 583 | { | ||
| 584 |
1/2✓ Branch 1 taken 118356 times.
✗ Branch 2 not taken.
|
118528 | const auto refElement1 = referenceElement(geo1); |
| 585 |
1/2✓ Branch 1 taken 118356 times.
✗ Branch 2 not taken.
|
118528 | 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 |
4/4✓ Branch 0 taken 473532 times.
✓ Branch 1 taken 118392 times.
✓ Branch 2 taken 473532 times.
✓ Branch 3 taken 118392 times.
|
1066600 | for (int i = 0; i < refElement1.size(dim1-1); ++i) |
| 591 | { | ||
| 592 | 474036 | const auto localEdgeGeom1 = refElement1.template geometry<dim1-1>(i); | |
| 593 |
6/9✓ Branch 2 taken 473532 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 473424 times.
✓ Branch 5 taken 108 times.
✓ Branch 6 taken 473424 times.
✓ Branch 7 taken 108 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 108 times.
✗ Branch 10 not taken.
|
2370180 | const auto edge1 = SegGeometry( Dune::GeometryTypes::line, |
| 594 |
1/4✓ Branch 3 taken 473424 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
1895688 | std::vector<Point>( {geo1.global(localEdgeGeom1.corner(0)), |
| 595 | geo1.global(localEdgeGeom1.corner(1))} )); | ||
| 596 | |||
| 597 |
4/4✓ Branch 0 taken 1894056 times.
✓ Branch 1 taken 473532 times.
✓ Branch 2 taken 1894056 times.
✓ Branch 3 taken 473532 times.
|
4266180 | for (int j = 0; j < refElement2.size(dim2-1); ++j) |
| 598 | { | ||
| 599 | 1896072 | const auto localEdgeGeom2 = refElement2.template geometry<dim2-1>(j); | |
| 600 |
6/11✓ Branch 2 taken 1894056 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1893696 times.
✓ Branch 5 taken 360 times.
✓ Branch 6 taken 1893696 times.
✓ Branch 7 taken 360 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 360 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
9480360 | const auto edge2 = SegGeometry( Dune::GeometryTypes::line, |
| 601 |
1/5✓ Branch 3 taken 1893696 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
7583568 | std::vector<Point>( {geo2.global(localEdgeGeom2.corner(0)), |
| 602 | geo2.global(localEdgeGeom2.corner(1))} )); | ||
| 603 | |||
| 604 | using EdgeTest = GeometryIntersection<SegGeometry, SegGeometry, PointPolicy>; | ||
| 605 | 1896072 | typename EdgeTest::Intersection edgeIntersection; | |
| 606 |
2/2✓ Branch 1 taken 429668 times.
✓ Branch 2 taken 1464388 times.
|
1896072 | if (EdgeTest::intersection(edge1, edge2, edgeIntersection)) |
| 607 |
1/2✓ Branch 1 taken 429668 times.
✗ Branch 2 not taken.
|
429988 | points.emplace_back(edgeIntersection); |
| 608 | } | ||
| 609 | } | ||
| 610 | } | ||
| 611 | } | ||
| 612 | |||
| 613 |
2/4✓ Branch 0 taken 195772 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 195772 times.
✗ Branch 3 not taken.
|
392040 | if (points.empty()) |
| 614 | return false; | ||
| 615 | |||
| 616 | // remove duplicates | ||
| 617 | 783896 | const auto eps = (geo1.corner(0) - geo1.corner(1)).two_norm()*eps_; | |
| 618 | 588060 | std::sort(points.begin(), points.end(), [&eps](const auto& a, const auto& b) -> bool | |
| 619 | { | ||
| 620 | using std::abs; | ||
| 621 |
48/208✗ 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 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 not taken.
✗ Branch 23 not taken.
✓ Branch 24 taken 315638 times.
✓ Branch 25 taken 855127 times.
✓ Branch 26 taken 315638 times.
✓ Branch 27 taken 855127 times.
✓ Branch 28 taken 315638 times.
✓ Branch 29 taken 855127 times.
✓ Branch 30 taken 315638 times.
✓ Branch 31 taken 855127 times.
✗ 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 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✓ Branch 48 taken 368 times.
✓ Branch 49 taken 364 times.
✓ Branch 50 taken 368 times.
✓ Branch 51 taken 364 times.
✓ Branch 52 taken 368 times.
✓ Branch 53 taken 364 times.
✓ Branch 54 taken 368 times.
✓ Branch 55 taken 364 times.
✗ 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 not taken.
✗ Branch 71 not taken.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
✗ Branch 94 not taken.
✗ Branch 95 not taken.
✓ Branch 96 taken 308332 times.
✓ Branch 97 taken 570504 times.
✓ Branch 98 taken 308332 times.
✓ Branch 99 taken 570504 times.
✓ Branch 100 taken 308332 times.
✓ Branch 101 taken 570504 times.
✓ Branch 102 taken 308332 times.
✓ Branch 103 taken 570504 times.
✗ Branch 104 not taken.
✗ Branch 105 not taken.
✗ Branch 106 not taken.
✗ Branch 107 not taken.
✗ 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 taken 136 times.
✓ Branch 121 taken 192 times.
✓ Branch 122 taken 136 times.
✓ Branch 123 taken 192 times.
✓ Branch 124 taken 136 times.
✓ Branch 125 taken 192 times.
✓ Branch 126 taken 136 times.
✓ Branch 127 taken 192 times.
✗ 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 not taken.
✗ Branch 141 not taken.
✗ Branch 142 not taken.
✗ Branch 143 not taken.
✗ Branch 144 not taken.
✗ Branch 145 not taken.
✗ Branch 146 not taken.
✗ Branch 147 not taken.
✗ Branch 148 not taken.
✗ Branch 149 not taken.
✗ Branch 150 not taken.
✗ Branch 151 not taken.
✗ Branch 152 not taken.
✗ Branch 153 not taken.
✗ Branch 154 not taken.
✗ Branch 155 not taken.
✗ 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 not taken.
✗ Branch 166 not taken.
✗ Branch 167 not taken.
✗ Branch 168 not taken.
✗ Branch 169 not taken.
✗ Branch 170 not taken.
✗ Branch 171 not taken.
✗ Branch 172 not taken.
✗ Branch 173 not taken.
✗ Branch 174 not taken.
✗ Branch 175 not taken.
✗ Branch 176 not taken.
✗ Branch 177 not taken.
✗ Branch 178 not taken.
✗ Branch 179 not taken.
✗ Branch 180 not taken.
✗ Branch 181 not taken.
✗ Branch 182 not taken.
✗ Branch 183 not taken.
✗ Branch 184 not taken.
✗ Branch 185 not taken.
✗ Branch 186 not taken.
✗ Branch 187 not taken.
✗ Branch 188 not taken.
✗ Branch 189 not taken.
✗ Branch 190 not taken.
✗ Branch 191 not taken.
✓ Branch 192 taken 372 times.
✓ Branch 193 taken 216 times.
✓ Branch 194 taken 372 times.
✓ Branch 195 taken 216 times.
✓ Branch 196 taken 372 times.
✓ Branch 197 taken 216 times.
✓ Branch 198 taken 372 times.
✓ Branch 199 taken 216 times.
✓ Branch 200 taken 112 times.
✓ Branch 201 taken 100 times.
✓ Branch 202 taken 112 times.
✓ Branch 203 taken 100 times.
✓ Branch 204 taken 112 times.
✓ Branch 205 taken 100 times.
✓ Branch 206 taken 112 times.
✓ Branch 207 taken 100 times.
|
8205844 | return (abs(a[0]-b[0]) > eps ? a[0] < b[0] : a[1] < b[1]); |
| 622 | }); | ||
| 623 | |||
| 624 | 2346532 | auto removeIt = std::unique(points.begin(), points.end(), [&eps](const auto& a, const auto&b) | |
| 625 | { | ||
| 626 | 1759272 | return (b-a).two_norm() < eps; | |
| 627 | }); | ||
| 628 | |||
| 629 |
8/8✓ Branch 0 taken 150872 times.
✓ Branch 1 taken 44900 times.
✓ Branch 2 taken 150872 times.
✓ Branch 3 taken 44900 times.
✓ Branch 4 taken 150872 times.
✓ Branch 5 taken 44900 times.
✓ Branch 6 taken 150872 times.
✓ Branch 7 taken 44900 times.
|
784080 | points.erase(removeIt, points.end()); |
| 630 | |||
| 631 | // return false if we don't have at least three unique points | ||
| 632 |
4/4✓ Branch 0 taken 78584 times.
✓ Branch 1 taken 117188 times.
✓ Branch 2 taken 78584 times.
✓ Branch 3 taken 117188 times.
|
392040 | if (points.size() < 3) |
| 633 | return false; | ||
| 634 | |||
| 635 | // intersection polygon is convex hull of above points | ||
| 636 |
2/4✓ Branch 1 taken 78584 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 78584 times.
|
78820 | intersection = grahamConvexHull<2>(points); |
| 637 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 78584 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 78584 times.
|
157640 | assert(!intersection.empty()); |
| 638 | return true; | ||
| 639 | } | ||
| 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 | 132832 | 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 | static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast) | ||
| 754 | { | ||
| 755 | // lambda to obtain facet corners on geo1 | ||
| 756 | 346620 | auto getFacetCorners = [] (const Geometry1& geo1) | |
| 757 | { | ||
| 758 |
5/7✓ Branch 0 taken 19955 times.
✓ Branch 1 taken 166905 times.
✓ Branch 2 taken 6400 times.
✓ Branch 3 taken 4720 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2432 times.
✗ Branch 6 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 |
9/10✓ Branch 0 taken 19955 times.
✓ Branch 1 taken 166905 times.
✓ Branch 2 taken 6400 times.
✓ Branch 3 taken 4720 times.
✓ Branch 4 taken 2432 times.
✓ Branch 5 taken 19955 times.
✓ Branch 6 taken 9040 times.
✓ Branch 7 taken 6400 times.
✓ Branch 8 taken 4720 times.
✗ Branch 9 not taken.
|
240527 | switch (geo1.corners()) |
| 761 | { | ||
| 762 | // todo compute intersection geometries | ||
| 763 | 19955 | case 8: // hexahedron | |
| 764 |
32/72✓ Branch 1 taken 177820 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 177820 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 177820 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 177820 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 177820 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 177820 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 177820 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 177820 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 177820 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 177820 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 177820 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 177820 times.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✓ Branch 38 taken 177820 times.
✓ Branch 39 taken 1066920 times.
✓ Branch 40 taken 177820 times.
✓ Branch 41 taken 1066920 times.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✓ Branch 54 taken 2432 times.
✗ Branch 55 not taken.
✓ Branch 57 taken 2432 times.
✗ Branch 58 not taken.
✓ Branch 60 taken 2432 times.
✗ Branch 61 not taken.
✓ Branch 63 taken 2432 times.
✗ Branch 64 not taken.
✓ Branch 66 taken 2432 times.
✗ Branch 67 not taken.
✓ Branch 69 taken 2432 times.
✗ Branch 70 not taken.
✓ Branch 72 taken 2432 times.
✗ Branch 73 not taken.
✓ Branch 75 taken 2432 times.
✗ Branch 76 not taken.
✓ Branch 78 taken 2432 times.
✗ Branch 79 not taken.
✓ Branch 81 taken 2432 times.
✗ Branch 82 not taken.
✓ Branch 84 taken 2432 times.
✗ Branch 85 not taken.
✓ Branch 87 taken 2432 times.
✗ Branch 88 not taken.
✗ Branch 90 not taken.
✓ Branch 91 taken 2432 times.
✓ Branch 92 taken 14592 times.
✓ Branch 93 taken 2432 times.
✓ Branch 94 taken 14592 times.
✗ Branch 95 not taken.
✗ Branch 102 not taken.
✗ Branch 103 not taken.
✗ Branch 104 not taken.
✗ Branch 105 not taken.
|
1602313 | 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 | 9040 | case 6: // prism | |
| 768 |
14/30✓ Branch 1 taken 9040 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9040 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9040 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 9040 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 9040 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 9040 times.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✓ Branch 32 taken 9040 times.
✓ Branch 33 taken 45200 times.
✓ Branch 34 taken 9040 times.
✓ Branch 35 taken 45200 times.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
|
63280 | facetCorners = {{0, 2, 1}, {3, 4, 5}, {0, 1, 3, 4}, {0, 3, 2, 5}, {1, 2, 4, 5}}; |
| 769 | 9040 | break; | |
| 770 | 6400 | case 5: // pyramid | |
| 771 |
14/30✓ Branch 1 taken 6400 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6400 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6400 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6400 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 6400 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 6400 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 6400 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 6400 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 6400 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 6400 times.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✓ Branch 32 taken 6400 times.
✓ Branch 33 taken 32000 times.
✓ Branch 34 taken 6400 times.
✓ Branch 35 taken 32000 times.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
|
44800 | facetCorners = {{0, 2, 1, 3}, {0, 1, 4}, {1, 3, 4}, {2, 4, 3}, {0, 4, 2}}; |
| 772 | 6400 | break; | |
| 773 | 4720 | case 4: // tetrahedron | |
| 774 |
12/26✓ Branch 1 taken 4720 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4720 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4720 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4720 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4720 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 4720 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 4720 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 4720 times.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 4720 times.
✓ Branch 27 taken 18880 times.
✓ Branch 28 taken 4720 times.
✓ Branch 29 taken 18880 times.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
|
28320 | facetCorners = {{1, 0, 2}, {0, 1, 3}, {0, 3, 2}, {1, 2, 3}}; |
| 775 | 4720 | break; | |
| 776 | ✗ | default: | |
| 777 | ✗ | DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type " | |
| 778 | << geo1.type() << ", "<< geo1.corners() << " corners."); | ||
| 779 | } | ||
| 780 | |||
| 781 | 400824 | return facetCorners; | |
| 782 | }; | ||
| 783 | |||
| 784 | // lambda to obtain the normal on a facet | ||
| 785 | 1613980 | 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 | 1633582 | 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 | 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 | 49 | const auto b1Norm2 = b1.two_norm2(); | |
| 872 | |||
| 873 | using std::max; | ||
| 874 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 1 times.
|
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 47 times.
✓ Branch 2 taken 2 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 | 47 | auto d1 = geo2.corner(0) - geo1.corner(0); | |
| 884 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 37 times.
|
94 | 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 |
4/4✓ Branch 0 taken 43 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 43 times.
✓ Branch 3 taken 4 times.
|
141 | if (abs(d1*n2) > eps3) |
| 892 | return false; | ||
| 893 | |||
| 894 | // the candidate intersection points | ||
| 895 |
3/6✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 43 times.
✗ Branch 7 not taken.
|
129 | std::vector<Point> points; points.reserve(8); |
| 896 | |||
| 897 | // add polygon1 corners that are inside polygon2 | ||
| 898 |
4/4✓ Branch 0 taken 130 times.
✓ Branch 1 taken 43 times.
✓ Branch 2 taken 130 times.
✓ Branch 3 taken 43 times.
|
216 | 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 |
2/2✓ Branch 0 taken 38 times.
✓ Branch 1 taken 5 times.
|
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 |
4/4✓ Branch 0 taken 126 times.
✓ Branch 1 taken 38 times.
✓ Branch 2 taken 126 times.
✓ Branch 3 taken 38 times.
|
202 | 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 |
2/4✓ Branch 0 taken 38 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 38 times.
✗ Branch 3 not taken.
|
76 | 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 |
4/4✓ Branch 0 taken 114 times.
✓ Branch 1 taken 38 times.
✓ Branch 2 taken 114 times.
✓ Branch 3 taken 38 times.
|
266 | for (int i = 0; i < referenceElement1.size(dim1-1); ++i) |
| 921 | { | ||
| 922 | 114 | const auto localEdgeGeom1 = referenceElement1.template geometry<dim1-1>(i); | |
| 923 |
4/8✓ Branch 2 taken 114 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 114 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 114 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 114 times.
✗ Branch 10 not taken.
|
570 | const auto edge1 = SegGeometry( |
| 924 | Dune::GeometryTypes::line, | ||
| 925 |
0/2✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
228 | std::vector<Point>({ |
| 926 | geo1.global(localEdgeGeom1.corner(0)), | ||
| 927 | geo1.global(localEdgeGeom1.corner(1)) | ||
| 928 | }) | ||
| 929 | ); | ||
| 930 | |||
| 931 |
4/4✓ Branch 0 taken 378 times.
✓ Branch 1 taken 114 times.
✓ Branch 2 taken 378 times.
✓ Branch 3 taken 114 times.
|
870 | for (int j = 0; j < referenceElement2.size(dim2-1); ++j) |
| 932 | { | ||
| 933 | 378 | const auto localEdgeGeom2 = referenceElement2.template geometry<dim2-1>(j); | |
| 934 |
4/10✓ Branch 2 taken 378 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 378 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 378 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 378 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
1890 | const auto edge2 = SegGeometry( |
| 935 | Dune::GeometryTypes::line, | ||
| 936 |
0/2✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
756 | std::vector<Point>({ |
| 937 | geo2.global(localEdgeGeom2.corner(0)), | ||
| 938 | 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 |
2/4✓ Branch 0 taken 43 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 43 times.
✗ Branch 3 not taken.
|
86 | if (points.empty()) |
| 952 | return false; | ||
| 953 | |||
| 954 | // remove duplicates | ||
| 955 | 43 | const auto eps = maxNorm*eps_; | |
| 956 | 129 | std::sort(points.begin(), points.end(), [eps] (const auto& a, const auto& b) -> bool | |
| 957 | { | ||
| 958 | using std::abs; | ||
| 959 |
16/104✗ 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 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 not taken.
✗ Branch 23 not taken.
✓ Branch 24 taken 60 times.
✓ Branch 25 taken 248 times.
✓ Branch 26 taken 60 times.
✓ Branch 27 taken 248 times.
✓ Branch 28 taken 60 times.
✓ Branch 29 taken 248 times.
✓ Branch 30 taken 60 times.
✓ Branch 31 taken 248 times.
✗ 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 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ 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 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 not taken.
✗ Branch 71 not taken.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
✗ Branch 94 not taken.
✗ Branch 95 not taken.
✓ Branch 96 taken 71 times.
✓ Branch 97 taken 140 times.
✓ Branch 98 taken 71 times.
✓ Branch 99 taken 140 times.
✓ Branch 100 taken 71 times.
✓ Branch 101 taken 140 times.
✓ Branch 102 taken 71 times.
✓ Branch 103 taken 140 times.
|
2076 | return (abs(a[0]-b[0]) > eps ? a[0] < b[0] |
| 960 |
12/78✗ 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 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 taken 84 times.
✓ Branch 19 taken 164 times.
✓ Branch 20 taken 84 times.
✓ Branch 21 taken 164 times.
✓ Branch 22 taken 84 times.
✓ Branch 23 taken 164 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 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 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ 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 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 not taken.
✗ Branch 71 not taken.
✓ Branch 72 taken 48 times.
✓ Branch 73 taken 92 times.
✓ Branch 74 taken 48 times.
✓ Branch 75 taken 92 times.
✓ Branch 76 taken 48 times.
✓ Branch 77 taken 92 times.
|
1164 | : (abs(a[1]-b[1]) > eps ? a[1] < b[1] |
| 961 |
4/26✗ 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 taken 84 times.
✓ Branch 7 taken 164 times.
✗ 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 not taken.
✗ Branch 23 not taken.
✓ Branch 24 taken 48 times.
✓ Branch 25 taken 92 times.
|
388 | : (a[2] < b[2]))); |
| 962 | }); | ||
| 963 | |||
| 964 | 43 | const auto squaredEps = eps*eps; | |
| 965 | 215 | 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 | points.end() | ||
| 969 | ); | ||
| 970 | |||
| 971 | // return false if we don't have at least three unique points | ||
| 972 |
4/4✓ Branch 0 taken 21 times.
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 22 times.
|
86 | if (points.size() < 3) |
| 973 | return false; | ||
| 974 | |||
| 975 | // intersection polygon is convex hull of above points | ||
| 976 |
2/4✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 21 times.
|
21 | intersection = grahamConvexHull<2>(points); |
| 977 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 21 times.
|
42 | assert(!intersection.empty()); |
| 978 | return true; | ||
| 979 | } | ||
| 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 | 42239 | 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 |
3/6✓ Branch 1 taken 34631 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 34631 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 34631 times.
✗ Branch 7 not taken.
|
168956 | std::vector<Point> points; points.reserve(10); |
| 1051 | |||
| 1052 | // add 3d geometry corners that are inside the 2d geometry | ||
| 1053 |
4/4✓ Branch 0 taken 249456 times.
✓ Branch 1 taken 34631 times.
✓ Branch 2 taken 249456 times.
✓ Branch 3 taken 34631 times.
|
394798 | for (int i = 0; i < geo1.corners(); ++i) |
| 1054 |
5/5✓ Branch 1 taken 215168 times.
✓ Branch 2 taken 34288 times.
✓ Branch 3 taken 10816 times.
✓ Branch 4 taken 209982 times.
✓ Branch 5 taken 19634 times.
|
592888 | if (intersectsPointGeometry(geo1.corner(i), geo2)) |
| 1055 |
2/3✓ Branch 1 taken 13304 times.
✓ Branch 2 taken 8118 times.
✗ Branch 3 not taken.
|
56358 | points.emplace_back(geo1.corner(i)); |
| 1056 | |||
| 1057 | // add 2d geometry corners that are inside the 3d geometry | ||
| 1058 |
4/4✓ Branch 0 taken 105037 times.
✓ Branch 1 taken 34631 times.
✓ Branch 2 taken 100525 times.
✓ Branch 3 taken 33503 times.
|
211211 | for (int i = 0; i < geo2.corners(); ++i) |
| 1059 |
5/5✓ Branch 1 taken 3014 times.
✓ Branch 2 taken 88675 times.
✓ Branch 3 taken 13348 times.
✓ Branch 4 taken 16 times.
✓ Branch 5 taken 60 times.
|
233808 | if (intersectsPointGeometry(geo2.corner(i), geo1)) |
| 1060 |
2/3✓ Branch 1 taken 3014 times.
✓ Branch 2 taken 8130 times.
✗ Branch 3 not taken.
|
14786 | 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 34611 times.
✗ Branch 2 not taken.
|
42239 | const auto refElement1 = referenceElement(geo1); |
| 1067 |
1/2✓ Branch 1 taken 23101 times.
✗ Branch 2 not taken.
|
42239 | 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 |
4/4✓ Branch 0 taken 374184 times.
✓ Branch 1 taken 34631 times.
✓ Branch 2 taken 374184 times.
✓ Branch 3 taken 34631 times.
|
973199 | for (int i = 0; i < refElement1.size(dim1-1); ++i) |
| 1074 | { | ||
| 1075 | 465480 | const auto localEdgeGeom = refElement1.template geometry<dim1-1>(i); | |
| 1076 | 930960 | const auto p = geo1.global(localEdgeGeom.corner(0)); | |
| 1077 | 930960 | const auto q = geo1.global(localEdgeGeom.corner(1)); | |
| 1078 |
6/16✓ Branch 1 taken 374184 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 374184 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 374184 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 374184 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 374184 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 374184 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
2792880 | const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q}); |
| 1079 | |||
| 1080 | using PolySegTest = GeometryIntersection<Geometry2, SegGeometry, PointPolicy>; | ||
| 1081 | 465480 | typename PolySegTest::Intersection polySegIntersection; | |
| 1082 |
3/4✓ Branch 1 taken 374184 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 51273 times.
✓ Branch 4 taken 322911 times.
|
465480 | if (PolySegTest::intersection(geo2, segGeo, polySegIntersection)) |
| 1083 |
1/2✓ Branch 1 taken 51273 times.
✗ Branch 2 not taken.
|
65225 | points.emplace_back(std::move(polySegIntersection)); |
| 1084 | } | ||
| 1085 | |||
| 1086 | // add intersection points of all polygon faces (codim 1) with the polyhedron faces | ||
| 1087 |
4/4✓ Branch 0 taken 193990 times.
✓ Branch 1 taken 34631 times.
✓ Branch 2 taken 193990 times.
✓ Branch 3 taken 34631 times.
|
521515 | for (int i = 0; i < refElement1.size(1); ++i) |
| 1088 | { | ||
| 1089 |
2/6✓ Branch 1 taken 193990 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 193990 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
764562 | const auto faceGeo = [&]() |
| 1090 | { | ||
| 1091 | 581970 | const auto localFaceGeo = refElement1.template geometry<1>(i); | |
| 1092 |
4/4✓ Branch 0 taken 120750 times.
✓ Branch 1 taken 27592 times.
✓ Branch 2 taken 120750 times.
✓ Branch 3 taken 27592 times.
|
387980 | if (localFaceGeo.corners() == 4) |
| 1093 | { | ||
| 1094 | 332796 | const auto a = geo1.global(localFaceGeo.corner(0)); | |
| 1095 | 332796 | const auto b = geo1.global(localFaceGeo.corner(1)); | |
| 1096 | 332796 | const auto c = geo1.global(localFaceGeo.corner(2)); | |
| 1097 | 332796 | const auto d = geo1.global(localFaceGeo.corner(3)); | |
| 1098 | |||
| 1099 |
4/10✓ Branch 1 taken 120750 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 120750 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 120750 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 120750 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
831990 | return PolyhedronFaceGeometry(Dune::GeometryTypes::cube(2), std::vector<Point>{a, b, c, d}); |
| 1100 | } | ||
| 1101 | else | ||
| 1102 | { | ||
| 1103 | 55184 | const auto a = geo1.global(localFaceGeo.corner(0)); | |
| 1104 | 55184 | const auto b = geo1.global(localFaceGeo.corner(1)); | |
| 1105 | 55184 | const auto c = geo1.global(localFaceGeo.corner(2)); | |
| 1106 | |||
| 1107 |
4/10✓ Branch 1 taken 27592 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27592 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 27592 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 27592 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
137960 | return PolyhedronFaceGeometry(Dune::GeometryTypes::simplex(2), std::vector<Point>{a, b, c}); |
| 1108 | } | ||
| 1109 |
2/2✓ Branch 0 taken 120750 times.
✓ Branch 1 taken 27592 times.
|
193990 | }(); |
| 1110 | |||
| 1111 |
4/4✓ Branch 0 taken 588834 times.
✓ Branch 1 taken 193990 times.
✓ Branch 2 taken 588834 times.
✓ Branch 3 taken 193990 times.
|
1693546 | for (int j = 0; j < refElement2.size(1); ++j) |
| 1112 | { | ||
| 1113 | 726954 | const auto localEdgeGeom = refElement2.template geometry<1>(j); | |
| 1114 | 726954 | const auto p = geo2.global(localEdgeGeom.corner(0)); | |
| 1115 | 726954 | const auto q = geo2.global(localEdgeGeom.corner(1)); | |
| 1116 | |||
| 1117 |
6/16✓ Branch 1 taken 588834 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 588834 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 588834 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 588834 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 588834 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 588834 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
4361724 | const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q}); |
| 1118 | |||
| 1119 | using PolySegTest = GeometryIntersection<PolyhedronFaceGeometry, SegGeometry, PointPolicy>; | ||
| 1120 | 726954 | typename PolySegTest::Intersection polySegIntersection; | |
| 1121 |
3/4✓ Branch 1 taken 588834 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 65786 times.
✓ Branch 4 taken 523048 times.
|
726954 | if (PolySegTest::intersection(faceGeo, segGeo, polySegIntersection)) |
| 1122 |
1/2✓ Branch 1 taken 65786 times.
✗ Branch 2 not taken.
|
71258 | points.emplace_back(std::move(polySegIntersection)); |
| 1123 | } | ||
| 1124 | } | ||
| 1125 | |||
| 1126 | // return if no intersection points were found | ||
| 1127 |
4/4✓ Branch 0 taken 19844 times.
✓ Branch 1 taken 14787 times.
✓ Branch 2 taken 19844 times.
✓ Branch 3 taken 14787 times.
|
84478 | if (points.empty()) return false; |
| 1128 | |||
| 1129 | // remove duplicates | ||
| 1130 | 58816 | const auto eps = (geo1.corner(0) - geo1.corner(1)).two_norm()*eps_; | |
| 1131 | 1264318 | const auto notEqual = [eps] (auto a, auto b) { using std::abs; return abs(b-a) > eps; }; | |
| 1132 | 71196 | std::sort(points.begin(), points.end(), [notEqual](const auto& a, const auto& b) -> bool | |
| 1133 | { | ||
| 1134 |
160/384✗ 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 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 74023 times.
✓ Branch 17 taken 122487 times.
✓ Branch 18 taken 74023 times.
✓ Branch 19 taken 122487 times.
✓ Branch 20 taken 74023 times.
✓ Branch 21 taken 122487 times.
✓ Branch 22 taken 74023 times.
✓ Branch 23 taken 122487 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 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 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ 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 64 taken 368 times.
✓ Branch 65 taken 364 times.
✓ Branch 66 taken 368 times.
✓ Branch 67 taken 364 times.
✓ Branch 68 taken 368 times.
✓ Branch 69 taken 364 times.
✓ Branch 70 taken 368 times.
✓ Branch 71 taken 364 times.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✓ Branch 88 taken 33178 times.
✓ Branch 89 taken 69931 times.
✓ Branch 90 taken 33178 times.
✓ Branch 91 taken 69931 times.
✓ Branch 92 taken 33178 times.
✓ Branch 93 taken 69931 times.
✓ Branch 94 taken 33178 times.
✓ Branch 95 taken 69931 times.
✗ Branch 96 not taken.
✗ 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 not taken.
✗ Branch 105 not taken.
✗ Branch 106 not taken.
✗ Branch 107 not taken.
✗ 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 taken 720 times.
✓ Branch 137 taken 1774 times.
✓ Branch 138 taken 720 times.
✓ Branch 139 taken 1774 times.
✓ Branch 140 taken 720 times.
✓ Branch 141 taken 1774 times.
✓ Branch 142 taken 720 times.
✓ Branch 143 taken 1774 times.
✓ Branch 144 taken 78 times.
✓ Branch 145 taken 69 times.
✓ Branch 146 taken 78 times.
✓ Branch 147 taken 69 times.
✓ Branch 148 taken 78 times.
✓ Branch 149 taken 69 times.
✓ Branch 150 taken 78 times.
✓ Branch 151 taken 69 times.
✓ Branch 152 taken 35 times.
✓ Branch 153 taken 23 times.
✓ Branch 154 taken 35 times.
✓ Branch 155 taken 23 times.
✓ Branch 156 taken 35 times.
✓ Branch 157 taken 23 times.
✓ Branch 158 taken 35 times.
✓ Branch 159 taken 23 times.
✓ Branch 160 taken 1 times.
✓ Branch 161 taken 9 times.
✓ Branch 162 taken 1 times.
✓ Branch 163 taken 9 times.
✓ Branch 164 taken 1 times.
✓ Branch 165 taken 9 times.
✓ Branch 166 taken 1 times.
✓ Branch 167 taken 9 times.
✓ Branch 168 taken 4 times.
✓ Branch 169 taken 4 times.
✓ Branch 170 taken 4 times.
✓ Branch 171 taken 4 times.
✓ Branch 172 taken 4 times.
✓ Branch 173 taken 4 times.
✓ Branch 174 taken 4 times.
✓ Branch 175 taken 4 times.
✓ Branch 176 taken 4 times.
✓ Branch 177 taken 2 times.
✓ Branch 178 taken 4 times.
✓ Branch 179 taken 2 times.
✓ Branch 180 taken 4 times.
✓ Branch 181 taken 2 times.
✓ Branch 182 taken 4 times.
✓ Branch 183 taken 2 times.
✓ Branch 184 taken 1 times.
✓ Branch 185 taken 1 times.
✓ Branch 186 taken 1 times.
✓ Branch 187 taken 1 times.
✓ Branch 188 taken 1 times.
✓ Branch 189 taken 1 times.
✓ Branch 190 taken 1 times.
✓ Branch 191 taken 1 times.
✓ Branch 192 taken 1 times.
✗ Branch 193 not taken.
✓ Branch 194 taken 1 times.
✗ Branch 195 not taken.
✓ Branch 196 taken 1 times.
✗ Branch 197 not taken.
✓ Branch 198 taken 1 times.
✗ Branch 199 not taken.
✗ Branch 200 not taken.
✗ Branch 201 not taken.
✗ Branch 202 not taken.
✗ Branch 203 not taken.
✗ Branch 204 not taken.
✗ Branch 205 not taken.
✗ Branch 206 not taken.
✗ Branch 207 not taken.
✓ Branch 208 taken 18140 times.
✓ Branch 209 taken 21376 times.
✓ Branch 210 taken 18140 times.
✓ Branch 211 taken 21376 times.
✓ Branch 212 taken 18140 times.
✓ Branch 213 taken 21376 times.
✓ Branch 214 taken 18140 times.
✓ Branch 215 taken 21376 times.
✗ Branch 216 not taken.
✗ Branch 217 not taken.
✗ Branch 218 not taken.
✗ Branch 219 not taken.
✗ Branch 220 not taken.
✗ Branch 221 not taken.
✗ Branch 222 not taken.
✗ Branch 223 not taken.
✗ Branch 224 not taken.
✗ Branch 225 not taken.
✗ Branch 226 not taken.
✗ Branch 227 not taken.
✗ Branch 228 not taken.
✗ Branch 229 not taken.
✗ Branch 230 not taken.
✗ Branch 231 not taken.
✗ Branch 232 not taken.
✗ Branch 233 not taken.
✗ Branch 234 not taken.
✗ Branch 235 not taken.
✗ Branch 236 not taken.
✗ Branch 237 not taken.
✗ Branch 238 not taken.
✗ Branch 239 not taken.
✗ Branch 240 not taken.
✗ Branch 241 not taken.
✗ Branch 242 not taken.
✗ Branch 243 not taken.
✗ Branch 244 not taken.
✗ Branch 245 not taken.
✗ Branch 246 not taken.
✗ Branch 247 not taken.
✗ Branch 248 not taken.
✗ Branch 249 not taken.
✗ Branch 250 not taken.
✗ Branch 251 not taken.
✗ Branch 252 not taken.
✗ Branch 253 not taken.
✗ Branch 254 not taken.
✗ Branch 255 not taken.
✗ Branch 256 not taken.
✗ Branch 257 not taken.
✗ Branch 258 not taken.
✗ Branch 259 not taken.
✗ Branch 260 not taken.
✗ Branch 261 not taken.
✗ Branch 262 not taken.
✗ Branch 263 not taken.
✗ Branch 264 not taken.
✗ Branch 265 not taken.
✗ Branch 266 not taken.
✗ Branch 267 not taken.
✗ Branch 268 not taken.
✗ Branch 269 not taken.
✗ Branch 270 not taken.
✗ Branch 271 not taken.
✗ Branch 272 not taken.
✗ Branch 273 not taken.
✗ Branch 274 not taken.
✗ Branch 275 not taken.
✗ Branch 276 not taken.
✗ Branch 277 not taken.
✗ Branch 278 not taken.
✗ Branch 279 not taken.
✓ Branch 280 taken 6399 times.
✓ Branch 281 taken 7477 times.
✓ Branch 282 taken 6399 times.
✓ Branch 283 taken 7477 times.
✓ Branch 284 taken 6399 times.
✓ Branch 285 taken 7477 times.
✓ Branch 286 taken 6399 times.
✓ Branch 287 taken 7477 times.
✓ Branch 288 taken 19 times.
✓ Branch 289 taken 32 times.
✓ Branch 290 taken 19 times.
✓ Branch 291 taken 32 times.
✓ Branch 292 taken 19 times.
✓ Branch 293 taken 32 times.
✓ Branch 294 taken 19 times.
✓ Branch 295 taken 32 times.
✓ Branch 296 taken 1 times.
✓ Branch 297 taken 28 times.
✓ Branch 298 taken 1 times.
✓ Branch 299 taken 28 times.
✓ Branch 300 taken 1 times.
✓ Branch 301 taken 28 times.
✓ Branch 302 taken 1 times.
✓ Branch 303 taken 28 times.
✗ Branch 304 not taken.
✓ Branch 305 taken 4 times.
✗ Branch 306 not taken.
✓ Branch 307 taken 4 times.
✗ Branch 308 not taken.
✓ Branch 309 taken 4 times.
✗ Branch 310 not taken.
✓ Branch 311 taken 4 times.
✗ Branch 312 not taken.
✗ Branch 313 not taken.
✗ Branch 314 not taken.
✗ Branch 315 not taken.
✗ Branch 316 not taken.
✗ Branch 317 not taken.
✗ Branch 318 not taken.
✗ Branch 319 not taken.
✗ Branch 320 not taken.
✗ Branch 321 not taken.
✗ Branch 322 not taken.
✗ Branch 323 not taken.
✗ Branch 324 not taken.
✗ Branch 325 not taken.
✗ Branch 326 not taken.
✗ Branch 327 not taken.
✗ Branch 328 not taken.
✓ Branch 329 taken 4 times.
✗ Branch 330 not taken.
✓ Branch 331 taken 4 times.
✗ Branch 332 not taken.
✓ Branch 333 taken 4 times.
✗ Branch 334 not taken.
✓ Branch 335 taken 4 times.
✗ Branch 336 not taken.
✓ Branch 337 taken 4 times.
✗ Branch 338 not taken.
✓ Branch 339 taken 4 times.
✗ Branch 340 not taken.
✓ Branch 341 taken 4 times.
✗ Branch 342 not taken.
✓ Branch 343 taken 4 times.
✗ Branch 344 not taken.
✗ Branch 345 not taken.
✗ Branch 346 not taken.
✗ Branch 347 not taken.
✗ Branch 348 not taken.
✗ Branch 349 not taken.
✗ Branch 350 not taken.
✗ Branch 351 not taken.
✓ Branch 352 taken 372 times.
✓ Branch 353 taken 216 times.
✓ Branch 354 taken 372 times.
✓ Branch 355 taken 216 times.
✓ Branch 356 taken 372 times.
✓ Branch 357 taken 216 times.
✓ Branch 358 taken 372 times.
✓ Branch 359 taken 216 times.
✓ Branch 360 taken 206 times.
✓ Branch 361 taken 866 times.
✓ Branch 362 taken 206 times.
✓ Branch 363 taken 866 times.
✓ Branch 364 taken 206 times.
✓ Branch 365 taken 866 times.
✓ Branch 366 taken 206 times.
✓ Branch 367 taken 866 times.
✓ Branch 368 taken 9500 times.
✓ Branch 369 taken 8172 times.
✓ Branch 370 taken 9500 times.
✓ Branch 371 taken 8172 times.
✓ Branch 372 taken 9500 times.
✓ Branch 373 taken 8172 times.
✓ Branch 374 taken 9500 times.
✓ Branch 375 taken 8172 times.
✓ Branch 376 taken 3035 times.
✓ Branch 377 taken 4261 times.
✓ Branch 378 taken 3035 times.
✓ Branch 379 taken 4261 times.
✓ Branch 380 taken 3035 times.
✓ Branch 381 taken 4261 times.
✓ Branch 382 taken 3035 times.
✓ Branch 383 taken 4261 times.
|
1532756 | return (notEqual(a[0], b[0]) ? a[0] < b[0] |
| 1135 |
87/288✗ 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 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 13759 times.
✓ Branch 13 taken 108728 times.
✓ Branch 14 taken 13759 times.
✓ Branch 15 taken 108728 times.
✓ Branch 16 taken 13759 times.
✓ Branch 17 taken 108728 times.
✗ 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 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 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✓ Branch 48 taken 364 times.
✗ Branch 49 not taken.
✓ Branch 50 taken 364 times.
✗ Branch 51 not taken.
✓ Branch 52 taken 364 times.
✗ 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 taken 10648 times.
✓ Branch 67 taken 59283 times.
✓ Branch 68 taken 10648 times.
✓ Branch 69 taken 59283 times.
✓ Branch 70 taken 10648 times.
✓ Branch 71 taken 59283 times.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
✗ Branch 94 not taken.
✗ Branch 95 not taken.
✗ Branch 96 not taken.
✗ 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 taken 1774 times.
✗ Branch 104 not taken.
✓ Branch 105 taken 1774 times.
✗ Branch 106 not taken.
✓ Branch 107 taken 1774 times.
✗ Branch 108 not taken.
✓ Branch 109 taken 69 times.
✗ Branch 110 not taken.
✓ Branch 111 taken 69 times.
✗ Branch 112 not taken.
✓ Branch 113 taken 69 times.
✗ Branch 114 not taken.
✓ Branch 115 taken 23 times.
✗ Branch 116 not taken.
✓ Branch 117 taken 23 times.
✗ Branch 118 not taken.
✓ Branch 119 taken 23 times.
✗ Branch 120 not taken.
✓ Branch 121 taken 9 times.
✗ Branch 122 not taken.
✓ Branch 123 taken 9 times.
✗ Branch 124 not taken.
✓ Branch 125 taken 9 times.
✗ Branch 126 not taken.
✓ Branch 127 taken 4 times.
✗ Branch 128 not taken.
✓ Branch 129 taken 4 times.
✗ Branch 130 not taken.
✓ Branch 131 taken 4 times.
✗ Branch 132 not taken.
✓ Branch 133 taken 2 times.
✗ Branch 134 not taken.
✓ Branch 135 taken 2 times.
✗ Branch 136 not taken.
✓ Branch 137 taken 2 times.
✗ Branch 138 not taken.
✓ Branch 139 taken 1 times.
✗ Branch 140 not taken.
✓ Branch 141 taken 1 times.
✗ Branch 142 not taken.
✓ Branch 143 taken 1 times.
✗ Branch 144 not taken.
✗ Branch 145 not taken.
✗ Branch 146 not taken.
✗ Branch 147 not taken.
✗ Branch 148 not taken.
✗ Branch 149 not taken.
✗ Branch 150 not taken.
✗ Branch 151 not taken.
✗ Branch 152 not taken.
✗ Branch 153 not taken.
✗ Branch 154 not taken.
✗ Branch 155 not taken.
✓ Branch 156 taken 10304 times.
✓ Branch 157 taken 11072 times.
✓ Branch 158 taken 10304 times.
✓ Branch 159 taken 11072 times.
✓ Branch 160 taken 10304 times.
✓ Branch 161 taken 11072 times.
✗ Branch 162 not taken.
✗ Branch 163 not taken.
✗ Branch 164 not taken.
✗ Branch 165 not taken.
✗ Branch 166 not taken.
✗ Branch 167 not taken.
✗ Branch 168 not taken.
✗ Branch 169 not taken.
✗ Branch 170 not taken.
✗ Branch 171 not taken.
✗ Branch 172 not taken.
✗ Branch 173 not taken.
✗ Branch 174 not taken.
✗ Branch 175 not taken.
✗ Branch 176 not taken.
✗ Branch 177 not taken.
✗ Branch 178 not taken.
✗ Branch 179 not taken.
✗ Branch 180 not taken.
✗ Branch 181 not taken.
✗ Branch 182 not taken.
✗ Branch 183 not taken.
✗ Branch 184 not taken.
✗ Branch 185 not taken.
✗ Branch 186 not taken.
✗ Branch 187 not taken.
✗ Branch 188 not taken.
✗ Branch 189 not taken.
✗ Branch 190 not taken.
✗ Branch 191 not taken.
✗ Branch 192 not taken.
✗ Branch 193 not taken.
✗ Branch 194 not taken.
✗ Branch 195 not taken.
✗ Branch 196 not taken.
✗ Branch 197 not taken.
✗ Branch 198 not taken.
✗ Branch 199 not taken.
✗ Branch 200 not taken.
✗ Branch 201 not taken.
✗ Branch 202 not taken.
✗ Branch 203 not taken.
✗ Branch 204 not taken.
✗ Branch 205 not taken.
✗ Branch 206 not taken.
✗ Branch 207 not taken.
✗ Branch 208 not taken.
✗ Branch 209 not taken.
✓ Branch 210 taken 1861 times.
✓ Branch 211 taken 5616 times.
✓ Branch 212 taken 1861 times.
✓ Branch 213 taken 5616 times.
✓ Branch 214 taken 1861 times.
✓ Branch 215 taken 5616 times.
✓ Branch 216 taken 23 times.
✓ Branch 217 taken 9 times.
✓ Branch 218 taken 23 times.
✓ Branch 219 taken 9 times.
✓ Branch 220 taken 23 times.
✓ Branch 221 taken 9 times.
✓ Branch 222 taken 21 times.
✓ Branch 223 taken 7 times.
✓ Branch 224 taken 21 times.
✓ Branch 225 taken 7 times.
✓ Branch 226 taken 21 times.
✓ Branch 227 taken 7 times.
✗ Branch 228 not taken.
✓ Branch 229 taken 4 times.
✗ Branch 230 not taken.
✓ Branch 231 taken 4 times.
✗ Branch 232 not taken.
✓ Branch 233 taken 4 times.
✗ Branch 234 not taken.
✗ Branch 235 not taken.
✗ Branch 236 not taken.
✗ Branch 237 not taken.
✗ Branch 238 not taken.
✗ Branch 239 not taken.
✗ Branch 240 not taken.
✗ Branch 241 not taken.
✗ Branch 242 not taken.
✗ Branch 243 not taken.
✗ Branch 244 not taken.
✗ Branch 245 not taken.
✗ Branch 246 not taken.
✓ Branch 247 taken 4 times.
✗ Branch 248 not taken.
✓ Branch 249 taken 4 times.
✗ Branch 250 not taken.
✓ Branch 251 taken 4 times.
✗ Branch 252 not taken.
✓ Branch 253 taken 4 times.
✗ Branch 254 not taken.
✓ Branch 255 taken 4 times.
✗ Branch 256 not taken.
✓ Branch 257 taken 4 times.
✗ Branch 258 not taken.
✗ Branch 259 not taken.
✗ Branch 260 not taken.
✗ Branch 261 not taken.
✗ Branch 262 not taken.
✗ Branch 263 not taken.
✓ Branch 264 taken 216 times.
✗ Branch 265 not taken.
✓ Branch 266 taken 216 times.
✗ Branch 267 not taken.
✓ Branch 268 taken 216 times.
✗ Branch 269 not taken.
✗ Branch 270 not taken.
✓ Branch 271 taken 866 times.
✗ Branch 272 not taken.
✓ Branch 273 taken 866 times.
✗ Branch 274 not taken.
✓ Branch 275 taken 866 times.
✓ Branch 276 taken 4756 times.
✓ Branch 277 taken 3416 times.
✓ Branch 278 taken 4756 times.
✓ Branch 279 taken 3416 times.
✓ Branch 280 taken 4756 times.
✓ Branch 281 taken 3416 times.
✓ Branch 282 taken 818 times.
✓ Branch 283 taken 3443 times.
✓ Branch 284 taken 818 times.
✓ Branch 285 taken 3443 times.
✓ Branch 286 taken 818 times.
✓ Branch 287 taken 3443 times.
|
711312 | : (notEqual(a[1], b[1]) ? a[1] < b[1] |
| 1136 |
29/96✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 13759 times.
✓ Branch 5 taken 108728 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 10648 times.
✓ Branch 23 taken 59283 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.
|
237104 | : (a[2] < b[2]))); |
| 1137 | }); | ||
| 1138 | |||
| 1139 | 23732 | const auto squaredEps = eps*eps; | |
| 1140 | 145332 | points.erase(std::unique( | |
| 1141 | points.begin(), points.end(), | ||
| 1142 | 232890 | [squaredEps] (const auto& a, const auto&b) { return (b-a).two_norm2() < squaredEps; }), | |
| 1143 | points.end() | ||
| 1144 | ); | ||
| 1145 | |||
| 1146 | // return false if we don't have more than three unique points | ||
| 1147 |
4/4✓ Branch 0 taken 13016 times.
✓ Branch 1 taken 6828 times.
✓ Branch 2 taken 13016 times.
✓ Branch 3 taken 6828 times.
|
47464 | if (points.size() < 3) return false; |
| 1148 | |||
| 1149 | // intersection polygon is convex hull of above points | ||
| 1150 |
2/4✓ Branch 1 taken 13016 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 13016 times.
|
15932 | intersection = grahamConvexHull<2>(points); |
| 1151 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 13016 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13016 times.
|
31864 | assert(!intersection.empty()); |
| 1152 | |||
| 1153 | return true; | ||
| 1154 | } | ||
| 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 | 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 | 12782 | 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 |
1/2✓ Branch 1 taken 124 times.
✗ Branch 2 not taken.
|
12782 | const auto refElement1 = referenceElement(geo1); |
| 1242 |
1/2✓ Branch 1 taken 124 times.
✗ Branch 2 not taken.
|
12782 | const auto refElement2 = referenceElement(geo2); |
| 1243 | |||
| 1244 | // the candidate intersection points | ||
| 1245 |
2/6✓ Branch 1 taken 6453 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6453 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
38346 | std::vector<Point> points; |
| 1246 |
3/6✓ Branch 1 taken 6453 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6453 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6453 times.
✗ Branch 8 not taken.
|
38346 | points.reserve(refElement1.size(2) + refElement2.size(2)); |
| 1247 | |||
| 1248 | // add corners inside the other geometry | ||
| 1249 | 25936 | const auto addPointIntersections = [&](const auto& g1, const auto& g2) | |
| 1250 | { | ||
| 1251 |
12/24✓ Branch 0 taken 5448 times.
✓ Branch 1 taken 681 times.
✓ Branch 2 taken 5448 times.
✓ Branch 3 taken 681 times.
✓ Branch 4 taken 1732 times.
✓ Branch 5 taken 433 times.
✓ Branch 6 taken 1732 times.
✓ Branch 7 taken 433 times.
✓ Branch 8 taken 94336 times.
✓ Branch 9 taken 11792 times.
✓ Branch 10 taken 94336 times.
✓ Branch 11 taken 11792 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.
|
127328 | for (int i = 0; i < g1.corners(); ++i) |
| 1252 |
7/14✓ Branch 1 taken 106 times.
✓ Branch 2 taken 4010 times.
✓ Branch 3 taken 1332 times.
✓ Branch 5 taken 762 times.
✓ Branch 6 taken 970 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 28008 times.
✓ Branch 9 taken 66328 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
199316 | if (const auto& c = g1.corner(i); intersectsPointGeometry(c, g2)) |
| 1253 |
0/8✗ Branch 2 not taken.
✗ Branch 3 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.
|
29528 | points.emplace_back(c); |
| 1254 | }; | ||
| 1255 | |||
| 1256 |
1/2✓ Branch 1 taken 6453 times.
✗ Branch 2 not taken.
|
12782 | addPointIntersections(geo1, geo2); |
| 1257 |
1/2✓ Branch 1 taken 6453 times.
✗ Branch 2 not taken.
|
12782 | 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 | 25564 | const auto addEdgeIntersections = [&](const auto& g1, const auto& g2, const auto& ref1, const auto& ref2) | |
| 1268 | { | ||
| 1269 |
12/12✓ Branch 0 taken 4086 times.
✓ Branch 1 taken 681 times.
✓ Branch 2 taken 4086 times.
✓ Branch 3 taken 681 times.
✓ Branch 4 taken 1732 times.
✓ Branch 5 taken 433 times.
✓ Branch 6 taken 1732 times.
✓ Branch 7 taken 433 times.
✓ Branch 8 taken 70752 times.
✓ Branch 9 taken 11792 times.
✓ Branch 10 taken 70752 times.
✓ Branch 11 taken 11792 times.
|
166046 | for (int i = 0; i < ref1.size(1); ++i) |
| 1270 | { | ||
| 1271 |
6/12✓ Branch 1 taken 4086 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2598 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1732 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1732 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 70752 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 70752 times.
✗ Branch 14 not taken.
|
304792 | const auto faceGeo = [&]() |
| 1272 | { | ||
| 1273 | 229710 | const auto localFaceGeo = ref1.template geometry<1>(i); | |
| 1274 |
2/4✓ Branch 0 taken 1488 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1488 times.
✗ Branch 3 not taken.
|
153140 | if (localFaceGeo.corners() == 4) |
| 1275 | { | ||
| 1276 | 149676 | const auto a = g1.global(localFaceGeo.corner(0)); | |
| 1277 | 149676 | const auto b = g1.global(localFaceGeo.corner(1)); | |
| 1278 | 149676 | const auto c = g1.global(localFaceGeo.corner(2)); | |
| 1279 | 149676 | const auto d = g1.global(localFaceGeo.corner(3)); | |
| 1280 | |||
| 1281 | return PolyhedronFaceGeometry( | ||
| 1282 | Dune::GeometryTypes::cube(2), std::vector<Point>{a, b, c, d} | ||
| 1283 |
4/10✓ Branch 1 taken 1488 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1488 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1488 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1488 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
374190 | ); |
| 1284 | } | ||
| 1285 | else | ||
| 1286 | { | ||
| 1287 | 3464 | const auto a = g1.global(localFaceGeo.corner(0)); | |
| 1288 | 3464 | const auto b = g1.global(localFaceGeo.corner(1)); | |
| 1289 | 3464 | const auto c = g1.global(localFaceGeo.corner(2)); | |
| 1290 | |||
| 1291 | return PolyhedronFaceGeometry( | ||
| 1292 | Dune::GeometryTypes::simplex(2), std::vector<Point>{a, b, c} | ||
| 1293 |
0/10✗ 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 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
8660 | ); |
| 1294 | } | ||
| 1295 |
1/2✓ Branch 0 taken 1488 times.
✗ Branch 1 not taken.
|
76570 | }(); |
| 1296 | |||
| 1297 |
12/12✓ Branch 0 taken 33444 times.
✓ Branch 1 taken 4086 times.
✓ Branch 2 taken 33444 times.
✓ Branch 3 taken 4086 times.
✓ Branch 4 taken 20784 times.
✓ Branch 5 taken 1732 times.
✓ Branch 6 taken 20784 times.
✓ Branch 7 taken 1732 times.
✓ Branch 8 taken 849024 times.
✓ Branch 9 taken 70752 times.
✓ Branch 10 taken 849024 times.
✓ Branch 11 taken 70752 times.
|
1883074 | for (int j = 0; j < ref2.size(2); ++j) |
| 1298 | { | ||
| 1299 | 903252 | const auto localEdgeGeom = ref2.template geometry<2>(j); | |
| 1300 | 1806504 | const auto p = g2.global(localEdgeGeom.corner(0)); | |
| 1301 | 1806504 | const auto q = g2.global(localEdgeGeom.corner(1)); | |
| 1302 | |||
| 1303 |
18/48✓ Branch 1 taken 33444 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 33444 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 33444 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 33444 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 33444 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 33444 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 21 taken 20784 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 20784 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 20784 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 20784 times.
✗ Branch 30 not taken.
✓ Branch 32 taken 20784 times.
✗ Branch 33 not taken.
✓ Branch 34 taken 20784 times.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✓ Branch 41 taken 849024 times.
✗ Branch 42 not taken.
✓ Branch 44 taken 849024 times.
✗ Branch 45 not taken.
✓ Branch 47 taken 849024 times.
✗ Branch 48 not taken.
✓ Branch 49 taken 849024 times.
✗ Branch 50 not taken.
✓ Branch 52 taken 849024 times.
✗ Branch 53 not taken.
✓ Branch 54 taken 849024 times.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
|
5419512 | const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q}); |
| 1304 | |||
| 1305 | using PolySegTest = GeometryIntersection<PolyhedronFaceGeometry, SegGeometry, PointPolicy>; | ||
| 1306 | 903252 | typename PolySegTest::Intersection polySegIntersection; | |
| 1307 |
9/12✓ Branch 1 taken 33444 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5797 times.
✓ Branch 4 taken 27647 times.
✓ Branch 6 taken 20784 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1234 times.
✓ Branch 9 taken 19550 times.
✓ Branch 11 taken 849024 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 109584 times.
✓ Branch 14 taken 739440 times.
|
903252 | if (PolySegTest::intersection(faceGeo, segGeo, polySegIntersection)) |
| 1308 |
3/6✓ Branch 1 taken 5797 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1234 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 109584 times.
✗ Branch 8 not taken.
|
116615 | points.emplace_back(std::move(polySegIntersection)); |
| 1309 | } | ||
| 1310 | } | ||
| 1311 | }; | ||
| 1312 | |||
| 1313 |
1/2✓ Branch 1 taken 6453 times.
✗ Branch 2 not taken.
|
12782 | addEdgeIntersections(geo1, geo2, refElement1, refElement2); |
| 1314 |
1/2✓ Branch 1 taken 6453 times.
✗ Branch 2 not taken.
|
12782 | addEdgeIntersections(geo2, geo1, refElement2, refElement1); |
| 1315 | |||
| 1316 | // return if no intersection points were found | ||
| 1317 |
4/4✓ Branch 0 taken 6407 times.
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 6407 times.
✓ Branch 3 taken 46 times.
|
25564 | if (points.empty()) |
| 1318 | return false; | ||
| 1319 | |||
| 1320 | // remove duplicates | ||
| 1321 | 36274 | const auto norm = (geo1.corner(0) - geo1.corner(1)).two_norm(); | |
| 1322 | 12690 | const auto eps = norm*eps_; | |
| 1323 | 2239566 | const auto notEqual = [eps] (auto a, auto b) { using std::abs; return abs(b-a) > eps; }; | |
| 1324 | 38070 | std::sort(points.begin(), points.end(), [notEqual](const auto& a, const auto& b) -> bool | |
| 1325 | { | ||
| 1326 |
160/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 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 2280 times.
✓ Branch 17 taken 8243 times.
✓ Branch 18 taken 2280 times.
✓ Branch 19 taken 8243 times.
✓ Branch 20 taken 2280 times.
✓ Branch 21 taken 8243 times.
✓ Branch 22 taken 2280 times.
✓ Branch 23 taken 8243 times.
✓ Branch 24 taken 1945 times.
✓ Branch 25 taken 3041 times.
✓ Branch 26 taken 1945 times.
✓ Branch 27 taken 3041 times.
✓ Branch 28 taken 1945 times.
✓ Branch 29 taken 3041 times.
✓ Branch 30 taken 1945 times.
✓ Branch 31 taken 3041 times.
✓ Branch 32 taken 3927 times.
✓ Branch 33 taken 8929 times.
✓ Branch 34 taken 3927 times.
✓ Branch 35 taken 8929 times.
✓ Branch 36 taken 3927 times.
✓ Branch 37 taken 8929 times.
✓ Branch 38 taken 3927 times.
✓ Branch 39 taken 8929 times.
✓ Branch 40 taken 224 times.
✓ Branch 41 taken 1161 times.
✓ Branch 42 taken 224 times.
✓ Branch 43 taken 1161 times.
✓ Branch 44 taken 224 times.
✓ Branch 45 taken 1161 times.
✓ Branch 46 taken 224 times.
✓ Branch 47 taken 1161 times.
✓ Branch 48 taken 122 times.
✓ Branch 49 taken 1107 times.
✓ Branch 50 taken 122 times.
✓ Branch 51 taken 1107 times.
✓ Branch 52 taken 122 times.
✓ Branch 53 taken 1107 times.
✓ Branch 54 taken 122 times.
✓ Branch 55 taken 1107 times.
✓ Branch 56 taken 25 times.
✓ Branch 57 taken 126 times.
✓ Branch 58 taken 25 times.
✓ Branch 59 taken 126 times.
✓ Branch 60 taken 25 times.
✓ Branch 61 taken 126 times.
✓ Branch 62 taken 25 times.
✓ Branch 63 taken 126 times.
✓ Branch 64 taken 55 times.
✓ Branch 65 taken 155 times.
✓ Branch 66 taken 55 times.
✓ Branch 67 taken 155 times.
✓ Branch 68 taken 55 times.
✓ Branch 69 taken 155 times.
✓ Branch 70 taken 55 times.
✓ Branch 71 taken 155 times.
✓ Branch 72 taken 22 times.
✓ Branch 73 taken 108 times.
✓ Branch 74 taken 22 times.
✓ Branch 75 taken 108 times.
✓ Branch 76 taken 22 times.
✓ Branch 77 taken 108 times.
✓ Branch 78 taken 22 times.
✓ Branch 79 taken 108 times.
✓ Branch 80 taken 16 times.
✓ Branch 81 taken 80 times.
✓ Branch 82 taken 16 times.
✓ Branch 83 taken 80 times.
✓ Branch 84 taken 16 times.
✓ Branch 85 taken 80 times.
✓ Branch 86 taken 16 times.
✓ Branch 87 taken 80 times.
✓ Branch 88 taken 475 times.
✓ Branch 89 taken 1183 times.
✓ Branch 90 taken 475 times.
✓ Branch 91 taken 1183 times.
✓ Branch 92 taken 475 times.
✓ Branch 93 taken 1183 times.
✓ Branch 94 taken 475 times.
✓ Branch 95 taken 1183 times.
✗ Branch 96 not taken.
✗ 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 73198 times.
✓ Branch 105 taken 224692 times.
✓ Branch 106 taken 73198 times.
✓ Branch 107 taken 224692 times.
✓ Branch 108 taken 73198 times.
✓ Branch 109 taken 224692 times.
✓ Branch 110 taken 73198 times.
✓ Branch 111 taken 224692 times.
✓ Branch 112 taken 20280 times.
✓ Branch 113 taken 69291 times.
✓ Branch 114 taken 20280 times.
✓ Branch 115 taken 69291 times.
✓ Branch 116 taken 20280 times.
✓ Branch 117 taken 69291 times.
✓ Branch 118 taken 20280 times.
✓ Branch 119 taken 69291 times.
✓ Branch 120 taken 30500 times.
✓ Branch 121 taken 71604 times.
✓ Branch 122 taken 30500 times.
✓ Branch 123 taken 71604 times.
✓ Branch 124 taken 30500 times.
✓ Branch 125 taken 71604 times.
✓ Branch 126 taken 30500 times.
✓ Branch 127 taken 71604 times.
✓ Branch 128 taken 3065 times.
✓ Branch 129 taken 5142 times.
✓ Branch 130 taken 3065 times.
✓ Branch 131 taken 5142 times.
✓ Branch 132 taken 3065 times.
✓ Branch 133 taken 5142 times.
✓ Branch 134 taken 3065 times.
✓ Branch 135 taken 5142 times.
✓ Branch 136 taken 1745 times.
✓ Branch 137 taken 4305 times.
✓ Branch 138 taken 1745 times.
✓ Branch 139 taken 4305 times.
✓ Branch 140 taken 1745 times.
✓ Branch 141 taken 4305 times.
✓ Branch 142 taken 1745 times.
✓ Branch 143 taken 4305 times.
✓ Branch 144 taken 1570 times.
✓ Branch 145 taken 3730 times.
✓ Branch 146 taken 1570 times.
✓ Branch 147 taken 3730 times.
✓ Branch 148 taken 1570 times.
✓ Branch 149 taken 3730 times.
✓ Branch 150 taken 1570 times.
✓ Branch 151 taken 3730 times.
✓ Branch 152 taken 625 times.
✓ Branch 153 taken 1532 times.
✓ Branch 154 taken 625 times.
✓ Branch 155 taken 1532 times.
✓ Branch 156 taken 625 times.
✓ Branch 157 taken 1532 times.
✓ Branch 158 taken 625 times.
✓ Branch 159 taken 1532 times.
✓ Branch 160 taken 125 times.
✓ Branch 161 taken 637 times.
✓ Branch 162 taken 125 times.
✓ Branch 163 taken 637 times.
✓ Branch 164 taken 125 times.
✓ Branch 165 taken 637 times.
✓ Branch 166 taken 125 times.
✓ Branch 167 taken 637 times.
✗ Branch 168 not taken.
✗ Branch 169 not taken.
✗ Branch 170 not taken.
✗ Branch 171 not taken.
✗ Branch 172 not taken.
✗ Branch 173 not taken.
✗ Branch 174 not taken.
✗ Branch 175 not taken.
✓ Branch 176 taken 1398 times.
✓ Branch 177 taken 2111 times.
✓ Branch 178 taken 1398 times.
✓ Branch 179 taken 2111 times.
✓ Branch 180 taken 1398 times.
✓ Branch 181 taken 2111 times.
✓ Branch 182 taken 1398 times.
✓ Branch 183 taken 2111 times.
✓ Branch 184 taken 18369 times.
✓ Branch 185 taken 69559 times.
✓ Branch 186 taken 18369 times.
✓ Branch 187 taken 69559 times.
✓ Branch 188 taken 18369 times.
✓ Branch 189 taken 69559 times.
✓ Branch 190 taken 18369 times.
✓ Branch 191 taken 69559 times.
|
2546808 | return (notEqual(a[0], b[0]) ? a[0] < b[0] |
| 1327 |
117/144✗ 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 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 2464 times.
✓ Branch 13 taken 5779 times.
✓ Branch 14 taken 2464 times.
✓ Branch 15 taken 5779 times.
✓ Branch 16 taken 2464 times.
✓ Branch 17 taken 5779 times.
✓ Branch 18 taken 1098 times.
✓ Branch 19 taken 1943 times.
✓ Branch 20 taken 1098 times.
✓ Branch 21 taken 1943 times.
✓ Branch 22 taken 1098 times.
✓ Branch 23 taken 1943 times.
✓ Branch 24 taken 2238 times.
✓ Branch 25 taken 6691 times.
✓ Branch 26 taken 2238 times.
✓ Branch 27 taken 6691 times.
✓ Branch 28 taken 2238 times.
✓ Branch 29 taken 6691 times.
✓ Branch 30 taken 98 times.
✓ Branch 31 taken 1063 times.
✓ Branch 32 taken 98 times.
✓ Branch 33 taken 1063 times.
✓ Branch 34 taken 98 times.
✓ Branch 35 taken 1063 times.
✓ Branch 36 taken 67 times.
✓ Branch 37 taken 1040 times.
✓ Branch 38 taken 67 times.
✓ Branch 39 taken 1040 times.
✓ Branch 40 taken 67 times.
✓ Branch 41 taken 1040 times.
✓ Branch 42 taken 19 times.
✓ Branch 43 taken 107 times.
✓ Branch 44 taken 19 times.
✓ Branch 45 taken 107 times.
✓ Branch 46 taken 19 times.
✓ Branch 47 taken 107 times.
✓ Branch 48 taken 62 times.
✓ Branch 49 taken 93 times.
✓ Branch 50 taken 62 times.
✓ Branch 51 taken 93 times.
✓ Branch 52 taken 62 times.
✓ Branch 53 taken 93 times.
✓ Branch 54 taken 18 times.
✓ Branch 55 taken 90 times.
✓ Branch 56 taken 18 times.
✓ Branch 57 taken 90 times.
✓ Branch 58 taken 18 times.
✓ Branch 59 taken 90 times.
✓ Branch 60 taken 4 times.
✓ Branch 61 taken 76 times.
✓ Branch 62 taken 4 times.
✓ Branch 63 taken 76 times.
✓ Branch 64 taken 4 times.
✓ Branch 65 taken 76 times.
✓ Branch 66 taken 200 times.
✓ Branch 67 taken 983 times.
✓ Branch 68 taken 200 times.
✓ Branch 69 taken 983 times.
✓ Branch 70 taken 200 times.
✓ Branch 71 taken 983 times.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✓ Branch 78 taken 60960 times.
✓ Branch 79 taken 163732 times.
✓ Branch 80 taken 60960 times.
✓ Branch 81 taken 163732 times.
✓ Branch 82 taken 60960 times.
✓ Branch 83 taken 163732 times.
✓ Branch 84 taken 16279 times.
✓ Branch 85 taken 53012 times.
✓ Branch 86 taken 16279 times.
✓ Branch 87 taken 53012 times.
✓ Branch 88 taken 16279 times.
✓ Branch 89 taken 53012 times.
✓ Branch 90 taken 13436 times.
✓ Branch 91 taken 58168 times.
✓ Branch 92 taken 13436 times.
✓ Branch 93 taken 58168 times.
✓ Branch 94 taken 13436 times.
✓ Branch 95 taken 58168 times.
✓ Branch 96 taken 2140 times.
✓ Branch 97 taken 3002 times.
✓ Branch 98 taken 2140 times.
✓ Branch 99 taken 3002 times.
✓ Branch 100 taken 2140 times.
✓ Branch 101 taken 3002 times.
✓ Branch 102 taken 520 times.
✓ Branch 103 taken 3785 times.
✓ Branch 104 taken 520 times.
✓ Branch 105 taken 3785 times.
✓ Branch 106 taken 520 times.
✓ Branch 107 taken 3785 times.
✓ Branch 108 taken 1442 times.
✓ Branch 109 taken 2288 times.
✓ Branch 110 taken 1442 times.
✓ Branch 111 taken 2288 times.
✓ Branch 112 taken 1442 times.
✓ Branch 113 taken 2288 times.
✓ Branch 114 taken 125 times.
✓ Branch 115 taken 1407 times.
✓ Branch 116 taken 125 times.
✓ Branch 117 taken 1407 times.
✓ Branch 118 taken 125 times.
✓ Branch 119 taken 1407 times.
✗ Branch 120 not taken.
✓ Branch 121 taken 637 times.
✗ Branch 122 not taken.
✓ Branch 123 taken 637 times.
✗ Branch 124 not taken.
✓ Branch 125 taken 637 times.
✗ 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 taken 406 times.
✓ Branch 133 taken 1705 times.
✓ Branch 134 taken 406 times.
✓ Branch 135 taken 1705 times.
✓ Branch 136 taken 406 times.
✓ Branch 137 taken 1705 times.
✓ Branch 138 taken 17014 times.
✓ Branch 139 taken 52545 times.
✓ Branch 140 taken 17014 times.
✓ Branch 141 taken 52545 times.
✓ Branch 142 taken 17014 times.
✓ Branch 143 taken 52545 times.
|
1430208 | : (notEqual(a[1], b[1]) ? a[1] < b[1] |
| 1328 |
39/48✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2464 times.
✓ Branch 5 taken 5779 times.
✓ Branch 6 taken 1098 times.
✓ Branch 7 taken 1943 times.
✓ Branch 8 taken 2238 times.
✓ Branch 9 taken 6691 times.
✓ Branch 10 taken 98 times.
✓ Branch 11 taken 1063 times.
✓ Branch 12 taken 67 times.
✓ Branch 13 taken 1040 times.
✓ Branch 14 taken 19 times.
✓ Branch 15 taken 107 times.
✓ Branch 16 taken 62 times.
✓ Branch 17 taken 93 times.
✓ Branch 18 taken 18 times.
✓ Branch 19 taken 90 times.
✓ Branch 20 taken 4 times.
✓ Branch 21 taken 76 times.
✓ Branch 22 taken 200 times.
✓ Branch 23 taken 983 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.
|
476736 | : (a[2] < b[2]))); |
| 1329 | }); | ||
| 1330 | |||
| 1331 | 12690 | const auto squaredEps = eps*eps; | |
| 1332 | 199142 | points.erase(std::unique( | |
| 1333 | points.begin(), points.end(), | ||
| 1334 | 143780 | [squaredEps] (const auto& a, const auto&b) { return (b-a).two_norm2() < squaredEps; }), | |
| 1335 | points.end() | ||
| 1336 | ); | ||
| 1337 | |||
| 1338 | // return false if we don't have more than four unique points (dim+1) | ||
| 1339 |
4/4✓ Branch 0 taken 3747 times.
✓ Branch 1 taken 2660 times.
✓ Branch 2 taken 3747 times.
✓ Branch 3 taken 2660 times.
|
25380 | 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 | 3747 | const bool coplanar = [&] | |
| 1344 | { | ||
| 1345 | 11080 | const auto p0 = points[0]; | |
| 1346 | 18735 | 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 | 3747 | const auto epsCoplanar = normal.two_norm()*norm*eps_; | |
| 1350 |
4/4✓ Branch 0 taken 196 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 196 times.
✓ Branch 3 taken 8 times.
|
7333 | for (int i = 3; i < points.size(); ++i) |
| 1351 | { | ||
| 1352 | 9820 | const auto ad = points[i] - p0; | |
| 1353 | using std::abs; | ||
| 1354 |
4/4✓ Branch 0 taken 108 times.
✓ Branch 1 taken 88 times.
✓ Branch 2 taken 108 times.
✓ Branch 3 taken 88 times.
|
9820 | if (abs(normal*ad) > epsCoplanar) |
| 1355 | 1324 | return false; | |
| 1356 | } | ||
| 1357 | |||
| 1358 | return true; | ||
| 1359 | 7378 | }(); | |
| 1360 | |||
| 1361 |
2/2✓ Branch 0 taken 1324 times.
✓ Branch 1 taken 2423 times.
|
7378 | 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 1324 times.
✗ Branch 2 not taken.
|
2540 | intersection = points; |
| 1366 | |||
| 1367 | return true; | ||
| 1368 | } | ||
| 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 | 354449 | 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 73770 times.
✓ Branch 2 taken 280618 times.
|
354449 | 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 | 295172 | intersection = Intersection({geo2.global(tfirst), geo2.global(tlast)}); | |
| 1424 | 73813 | 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 | 2705188 | 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 | 2705188 | const auto p = geo2.corner(0); | |
| 1445 | 2705188 | const auto q = geo2.corner(1); | |
| 1446 | |||
| 1447 | 2705188 | const auto a = geo1.corner(0); | |
| 1448 | 2705188 | const auto b = geo1.corner(1); | |
| 1449 | 2705188 | const auto c = geo1.corner(2); | |
| 1450 | |||
| 1451 |
4/4✓ Branch 0 taken 464047 times.
✓ Branch 1 taken 1388731 times.
✓ Branch 2 taken 464047 times.
✓ Branch 3 taken 1388731 times.
|
5383304 | if (geo1.corners() == 3) |
| 1452 | 783135 | return intersect_<Policy>(a, b, c, p, q, is); | |
| 1453 | |||
| 1454 |
2/4✓ Branch 0 taken 1388731 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1388731 times.
✗ Branch 3 not taken.
|
3817034 | else if (geo1.corners() == 4) |
| 1455 | { | ||
| 1456 | 1922053 | Intersection is1, is2; | |
| 1457 | bool hasSegment1, hasSegment2; | ||
| 1458 | |||
| 1459 | 1922053 | const auto d = geo1.corner(3); | |
| 1460 | 1922053 | const bool intersects1 = intersect_<Policy>(a, b, d, p, q, is1, hasSegment1); | |
| 1461 | 1922053 | const bool intersects2 = intersect_<Policy>(a, d, c, p, q, is2, hasSegment2); | |
| 1462 | |||
| 1463 |
4/4✓ Branch 0 taken 1374151 times.
✓ Branch 1 taken 28116 times.
✓ Branch 2 taken 1346312 times.
✓ Branch 3 taken 27839 times.
|
1922053 | if (hasSegment1 || hasSegment2) |
| 1464 | return false; | ||
| 1465 | |||
| 1466 |
2/2✓ Branch 0 taken 97613 times.
✓ Branch 1 taken 1248699 times.
|
1861312 | if (intersects1) { is = is1; return true; } |
| 1467 |
2/2✓ Branch 0 taken 48264 times.
✓ Branch 1 taken 1200435 times.
|
1746745 | 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 | 354449 | static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast) | |
| 1483 | { | ||
| 1484 | // lambda to obtain the facet corners on geo1 | ||
| 1485 | 354388 | auto getFacetCorners = [] (const Geometry1& geo1) | |
| 1486 | { | ||
| 1487 |
3/10✗ Branch 0 not taken.
✓ Branch 1 taken 354364 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 13 times.
✓ Branch 6 taken 11 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
354388 | std::vector< std::array<int, 2> > facetCorners; |
| 1488 |
6/12✗ Branch 0 not taken.
✓ Branch 1 taken 354364 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 354364 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 13 times.
✓ Branch 7 taken 11 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 13 times.
✓ Branch 10 taken 11 times.
✗ Branch 11 not taken.
|
708776 | 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 | break; | ||
| 1493 | 354375 | case 3: // triangle | |
| 1494 |
2/4✓ Branch 1 taken 354364 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
|
354375 | facetCorners = {{1, 0}, {0, 2}, {2, 1}}; |
| 1495 | 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 | 1063164 | return facetCorners; | |
| 1502 | }; | ||
| 1503 | |||
| 1504 | 354449 | const auto center1 = geo1.center(); | |
| 1505 | 1063299 | const auto normal1 = crossProduct(geo1.corner(1) - geo1.corner(0), | |
| 1506 | 1417700 | geo1.corner(2) - geo1.corner(0)); | |
| 1507 | |||
| 1508 | // lambda to obtain the normal on a facet | ||
| 1509 | 3330613 | auto computeNormal = [¢er1, &normal1, &geo1] (const std::array<int, 2>& facetCorners) | |
| 1510 | { | ||
| 1511 | 1488082 | const auto c0 = geo1.corner(facetCorners[0]); | |
| 1512 | 1488082 | const auto c1 = geo1.corner(facetCorners[1]); | |
| 1513 | 744041 | const auto edge = c1 - c0; | |
| 1514 | |||
| 1515 | 744041 | Dune::FieldVector<ctype, dimworld> n = crossProduct(edge, normal1); | |
| 1516 | 1488082 | 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 743957 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 84 times.
✗ Branch 3 not taken.
|
1488082 | if ( n*(center1-c0) > 0.0 ) |
| 1521 | n *= -1.0; | ||
| 1522 | |||
| 1523 | 744041 | return n; | |
| 1524 | }; | ||
| 1525 | |||
| 1526 | 354449 | 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 | 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 | 3268581 | 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 | 3268581 | hasSegmentIs = false; | |
| 1551 | |||
| 1552 | 3268581 | const auto ab = b - a; | |
| 1553 | 3268581 | const auto ac = c - a; | |
| 1554 | 3268581 | const auto qp = p - q; | |
| 1555 | |||
| 1556 | // compute the triangle normal that defines the triangle plane | ||
| 1557 | 3268581 | 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 | 3268581 | const auto denom = normal*qp; | |
| 1562 | 3268581 | const auto abnorm2 = ab.two_norm2(); | |
| 1563 | 3268581 | const auto eps = eps_*abnorm2*qp.two_norm(); | |
| 1564 | using std::abs; | ||
| 1565 |
4/4✓ Branch 0 taken 1157013 times.
✓ Branch 1 taken 752908 times.
✓ Branch 2 taken 1157013 times.
✓ Branch 3 taken 752908 times.
|
6537162 | if (abs(denom) < eps) |
| 1566 | { | ||
| 1567 | 1373701 | const auto pa = a - p; | |
| 1568 | 1373701 | const auto denom2 = normal*pa; | |
| 1569 |
2/2✓ Branch 0 taken 292120 times.
✓ Branch 1 taken 864893 times.
|
4121103 | 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 | 354364 | SegmentIntersectionType segmentIs; | |
| 1580 | |||
| 1581 | 354364 | Triangle triangle(Dune::GeometryTypes::simplex(2), std::array<Point, 3>({a, b, c})); | |
| 1582 | 354364 | Segment segment(Dune::GeometryTypes::simplex(1), std::array<Point, 2>({p, q})); | |
| 1583 |
2/2✓ Branch 1 taken 58012 times.
✓ Branch 2 taken 234108 times.
|
354364 | if (SegmentIntersectionAlgorithm::intersection(triangle, segment, segmentIs)) |
| 1584 | { | ||
| 1585 | 73750 | hasSegmentIs = true; | |
| 1586 | 73750 | 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.
|
989344 | for (const auto& ip : {p, q}) |
| 1592 | { | ||
| 1593 | 516658 | if ( intersectsPointSimplex(ip, a, b) | |
| 1594 |
2/2✓ Branch 1 taken 344380 times.
✓ Branch 2 taken 25067 times.
|
461569 | || 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.
|
953112 | || intersectsPointSimplex(ip, c, a) ) |
| 1596 | { | ||
| 1597 | 88542 | is = ip; | |
| 1598 | 88542 | return true; | |
| 1599 | } | ||
| 1600 | } | ||
| 1601 | |||
| 1602 | 192072 | 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 | 1894880 | const auto ap = p - a; | |
| 1608 | 1894880 | const auto t = (ap*normal)/denom; | |
| 1609 |
2/2✓ Branch 0 taken 597679 times.
✓ Branch 1 taken 155229 times.
|
1894880 | if (t < 0.0 - eps_) return false; |
| 1610 |
2/2✓ Branch 0 taken 449497 times.
✓ Branch 1 taken 148182 times.
|
1436115 | 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 | 987776 | const auto e = crossProduct(qp, ap); | |
| 1615 | 987776 | const auto v = (ac*e)/denom; | |
| 1616 |
4/4✓ Branch 0 taken 333450 times.
✓ Branch 1 taken 116047 times.
✓ Branch 2 taken 262013 times.
✓ Branch 3 taken 71437 times.
|
987776 | if (v < -eps_ || v > 1.0 + eps_) return false; |
| 1617 | 423202 | const auto w = -(ab*e)/denom; | |
| 1618 |
4/4✓ Branch 0 taken 201694 times.
✓ Branch 1 taken 60319 times.
✓ Branch 2 taken 165820 times.
✓ Branch 3 taken 35874 times.
|
423202 | 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 | 225203 | const auto u = 1.0 - v - w; | |
| 1623 | |||
| 1624 | 225203 | Point ip(0.0); | |
| 1625 | ip.axpy(u, a); | ||
| 1626 | ip.axpy(v, b); | ||
| 1627 | 225203 | ip.axpy(w, c); | |
| 1628 | 225203 | is = ip; | |
| 1629 | 225203 | 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 | 418 | 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 1 taken 100 times.
✓ Branch 2 taken 318 times.
|
836 | if ( n.two_norm2() < eps2*v1Norm2 ) |
| 1700 | { | ||
| 1701 | // check if they lie on the same line | ||
| 1702 |
2/2✓ Branch 1 taken 51 times.
✓ Branch 2 taken 49 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.
|
10 | 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.
|
84 | 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.
|
76 | 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 | 636 | const auto t2 = - 1.0*(ac*v1Normal) / (v2*v1Normal); | |
| 1734 | |||
| 1735 | // check if the local coords are valid | ||
| 1736 |
4/4✓ Branch 0 taken 254 times.
✓ Branch 1 taken 64 times.
✓ Branch 2 taken 214 times.
✓ Branch 3 taken 40 times.
|
318 | if (t2 < -1.0*eps_ || t2 > 1.0 + eps_) |
| 1737 | return false; | ||
| 1738 | |||
| 1739 |
2/2✓ Branch 3 taken 162 times.
✓ Branch 4 taken 52 times.
|
428 | 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 | 478 | const auto pqNorm2 = pq.two_norm2(); | |
| 1767 | |||
| 1768 | using std::max; | ||
| 1769 |
2/2✓ Branch 0 taken 35 times.
✓ Branch 1 taken 443 times.
|
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 474 times.
✓ Branch 2 taken 4 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 470 times.
✓ Branch 2 taken 4 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 334 times.
✓ Branch 1 taken 136 times.
|
470 | tp = clamp(tp, 0.0, abNorm2); |
| 1796 |
2/2✓ Branch 0 taken 466 times.
✓ Branch 1 taken 4 times.
|
470 | tq = clamp(tq, 0.0, abNorm2); |
| 1797 | |||
| 1798 |
4/4✓ Branch 0 taken 338 times.
✓ Branch 1 taken 132 times.
✓ Branch 2 taken 338 times.
✓ Branch 3 taken 132 times.
|
940 | if (abs(tp-tq) < eps2) |
| 1799 | return false; | ||
| 1800 | |||
| 1801 | 676 | 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 |