geometryintersection.hh 51.2 KB
Newer Older
1
2
3
4
5
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
6
 *   the Free Software Foundation, either version 3 of the License, or       *
7
8
9
10
11
12
13
14
15
16
17
18
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
19
 * \ingroup Geometry
20
 * \brief A class for collision detection of two geometries
21
 *        and computation of intersection corners
22
 */
23
24
#ifndef DUMUX_GEOMETRY_INTERSECTION_HH
#define DUMUX_GEOMETRY_INTERSECTION_HH
25

26
#include <dune/common/deprecated.hh>
27
#include <dune/common/exceptions.hh>
28
#include <dune/common/promotiontraits.hh>
29
#include <dune/geometry/referenceelements.hh>
30
#include <dune/geometry/multilineargeometry.hh>
31
32

#include <dumux/common/math.hh>
33
34
#include <dumux/common/geometry/intersectspointgeometry.hh>
#include <dumux/common/geometry/grahamconvexhull.hh>
35
#include <dumux/common/geometry/boundingboxtree.hh>
36

37
namespace Dumux {
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
namespace IntersectionPolicy {

//! Policy structure for point-like intersections
template<class ct, int dw>
struct PointPolicy
{
    using ctype = ct;
    using Point = Dune::FieldVector<ctype, dw>;

    static constexpr int dimWorld = dw;
    static constexpr int dimIntersection = 0;
    using Intersection = Point;
};

//! Policy structure for segment-like intersections
template<class ct, int dw>
struct SegmentPolicy
{
    using ctype = ct;
    using Point = Dune::FieldVector<ctype, dw>;
    using Segment = std::array<Point, 2>;

    static constexpr int dimWorld = dw;
    static constexpr int dimIntersection = 1;
    using Intersection = Segment;
};

//! Policy structure for polygon-like intersections
template<class ct, int dw>
struct PolygonPolicy
{
    using ctype = ct;
    using Point = Dune::FieldVector<ctype, dw>;
    using Polygon = std::vector<Point>;

    static constexpr int dimWorld = dw;
    static constexpr int dimIntersection = 2;
    using Intersection = Polygon;
};

//! default policy chooser class
template<class Geometry1, class Geometry2, int dim2 = Geometry2::mydimension>
class DefaultPolicyChooser
{
    static constexpr int dimworld = Geometry1::coorddimension;
83
    static constexpr int isDim = std::min( int(Geometry1::mydimension), int(Geometry2::mydimension) );
84
    static_assert(dimworld == int(Geometry2::coorddimension), "Geometries must have the same coordinate dimension!");
85
    static_assert(isDim == 1 || isDim == 2, "No chooser class specialization available for given dimensions");
86
87
    using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
public:
88
    using type = std::conditional_t< isDim == 2,
89
90
91
92
93
94
95
96
97
                                     PolygonPolicy<ctype, Geometry1::coorddimension>,
                                     SegmentPolicy<ctype, Geometry1::coorddimension> >;
};

//! Helper alias to define the default intersection policy
template<class Geometry1, class Geometry2>
using DefaultPolicy = typename DefaultPolicyChooser<Geometry1, Geometry2>::type;

} // end namespace IntersectionPolicy
98

99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
namespace Impl {

/*!
 * \ingroup Geometry
 * \brief Algorithm to find segment-like intersections of a polgon/polyhedron with a
 *        segment. The result is stored in the form of the local coordinates tfirst
 *        and tlast on the segment geo1.
 * \param geo1 the first geometry
 * \param geo2 the second geometry (dimGeo2 < dimGeo1)
 * \param baseEps the base epsilon used for floating point comparisons
 * \param tfirst stores the local coordinate of beginning of intersection segment
 * \param tlast stores the local coordinate of the end of intersection segment
 * \param getFacetCornerIndices Function to obtain the facet corner indices on geo1
 * \param computeNormal Function to obtain the normal vector on a facet on geo1
 */
template< class Geo1, class Geo2, class ctype,
          class GetFacetCornerIndices, class ComputeNormalFunction >
bool computeSegmentIntersection(const Geo1& geo1, const Geo2& geo2, ctype baseEps,
                                ctype& tfirst, ctype& tlast,
                                const GetFacetCornerIndices& getFacetCornerIndices,
                                const ComputeNormalFunction& computeNormal)
{
    static_assert(int(Geo2::mydimension) == 1, "Geometry2 must be a segment");
    static_assert(int(Geo1::mydimension) > int(Geo2::mydimension),
                  "Dimension of Geometry1 must be higher than that of Geometry2");

    const auto a = geo2.corner(0);
    const auto b = geo2.corner(1);
    const auto d = b - a;

    // The initial interval is the whole segment.
    // Afterwards we start clipping the interval
    // at the intersections with the facets of geo1
    tfirst = 0.0;
    tlast = 1.0;

    const auto facets = getFacetCornerIndices(geo1);
    for (const auto& f : facets)
    {
        const auto n = computeNormal(f);

        const auto c0 = geo1.corner(f[0]);
        const ctype denom = n*d;
        const ctype dist = n*(a-c0);

        // use first edge of the facet to scale eps
        const auto edge1 = geo1.corner(f[1]) - geo1.corner(f[0]);
        const ctype eps = baseEps*edge1.two_norm();

        // if denominator is zero the segment in parallel to the edge.
        // If the distance is positive there is no intersection
        using std::abs;
        if (abs(denom) < eps)
        {
            if (dist > eps)
                return false;
        }
        else // not parallel: compute line-line intersection
        {
            const ctype t = -dist / denom;
            // if entering half space cut tfirst if t is larger
            using std::signbit;
            if (signbit(denom))
            {
                if (t > tfirst)
                    tfirst = t;
            }
            // if exiting half space cut tlast if t is smaller
            else
            {
                if (t < tlast)
                    tlast = t;
            }
            // there is no intersection of the interval is empty
            // use unscaled epsilon since t is in local coordinates
            if (tfirst > tlast - baseEps)
                return false;
        }
    }

    // If we made it until here an intersection exists.
    return true;
}

} // end namespace Impl

185
/*!
186
 * \ingroup Geometry
187
188
189
190
 * \brief A class for geometry collision detection and intersection calculation
 * The class can be specialized for combinations of dimworld, dim1, dim2, where
 * dimworld is the world dimension embedding a grid of dim1 and a grid of dim2.
 */
191
192
template
<class Geometry1, class Geometry2,
193
194
195
196
 class Policy = IntersectionPolicy::DefaultPolicy<Geometry1, Geometry2>,
 int dimworld = Geometry1::coorddimension,
 int dim1 = Geometry1::mydimension,
 int dim2 = Geometry2::mydimension>
197
class GeometryIntersection
198
{
199
200
    static constexpr int dimWorld = Policy::dimWorld;

201
public:
202
203
204
205
206
207
208
209
210
    using ctype = typename Policy::ctype;
    using Point = typename Policy::Point;
    using Intersection = typename Policy::Intersection;

    //! Deprecated alias, will be removed after 3.1
    using IntersectionType DUNE_DEPRECATED_MSG("Please use Intersection instead") = Intersection;

    //! Determine if the two geometries intersect and compute the intersection geometry
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
211
    {
212
213
        static_assert(dimworld == Geometry2::coorddimension, "Can only intersect geometries of same coordinate dimension");
        DUNE_THROW(Dune::NotImplemented, "Geometry intersection detection with intersection computation not implemented for dimworld = "
214
215
216
217
                                          << dimworld << ", dim1 = " << dim1 << ", dim2 = " << dim2);
    }
};

218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
/*!
 * \ingroup Geometry
 * \brief A class for segment--segment intersection in 2d space
 */
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 1, 1>
{
    enum { dimworld = 2 };
    enum { dim1 = 1 };
    enum { dim2 = 1 };

    // base epsilon for floating point comparisons
    static constexpr typename Policy::ctype eps_ = 1.5e-7;

public:
    using ctype = typename Policy::ctype;
    using Point = typename Policy::Point;
    using Intersection = typename Policy::Intersection;

    //! Deprecated alias, will be removed after 3.1
    using IntersectionType DUNE_DEPRECATED_MSG("Please use Intersection instead") = Intersection;

    /*!
     * \brief Colliding two segments
     * \param geo1/geo2 The geometries to intersect
     * \param intersection The intersection point
     * \note This overload is used when point-like intersections are seeked
     */
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
    {
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");

        const auto v1 = geo1.corner(1) - geo1.corner(0);
        const auto v2 = geo2.corner(1) - geo2.corner(0);
        const auto ac = geo2.corner(0) - geo1.corner(0);

        auto n2 = Point({-1.0*v2[1], v2[0]});
        n2 /= n2.two_norm();

        // compute distance of first corner on geo2 to line1
        const auto dist12 = n2*ac;

        // first check parallel segments
        using std::abs;
        using std::sqrt;

        const auto v1Norm2 = v1.two_norm2();
        const auto eps = eps_*sqrt(v1Norm2);
        const auto eps2 = eps_*v1Norm2;

        const auto sp = n2*v1;
        if (abs(sp) < eps)
        {
            if (abs(dist12) > eps)
                return false;

            // point intersection can only be one of the corners
            if ( (v1*v2) < 0.0 )
            {
                if ( ac.two_norm2() < eps2 )
                { intersection = geo2.corner(0); return true; }

                if ( (geo2.corner(1) - geo1.corner(1)).two_norm2() < eps2 )
                { intersection = geo2.corner(1); return true; }
            }
            else
            {
                if ( (geo2.corner(1) - geo1.corner(0)).two_norm2() < eps2 )
                { intersection = geo2.corner(1); return true; }

                if ( (geo2.corner(0) - geo1.corner(1)).two_norm2() < eps2 )
                { intersection = geo2.corner(0); return true; }
            }

            // no intersection
            return false;
        }

        // intersection point on v1 in local coords
        const auto t1 = dist12 / sp;

        // check if the local coords are valid
        if (t1 < -1.0*eps_ || t1 > 1.0 + eps_)
            return false;

        // compute global coordinates
        auto isPoint = geo1.global(t1);

        // check if point is in bounding box of v2
        const auto c = geo2.corner(0);
        const auto d = geo2.corner(1);

        using std::min; using std::max;
        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]) });

        if ( intersectsPointBoundingBox(isPoint, bBox.data()) )
        {
            intersection = std::move(isPoint);
            return true;
        }

        return false;
    }

    /*!
     * \brief Colliding two segments
     * \param geo1/geo2 The geometries to intersect
     * \param intersection Container to store corners of intersection segment
     * \note this overload is used when segment-like intersections are seeked
     * \todo implement this query
     */
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
    {
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
        DUNE_THROW(Dune::NotImplemented, "segment-segment intersection detection for segment-like intersections");
    }
};

338
339
/*!
 * \ingroup Geometry
Dennis Gläser's avatar
Dennis Gläser committed
340
 * \brief A class for polygon--segment intersection in 2d space
341
 */
342
343
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 2, 1>
344
345
346
347
348
349
{
    enum { dimworld = 2 };
    enum { dim1 = 2 };
    enum { dim2 = 1 };

public:
350
351
352
353
354
355
    using ctype = typename Policy::ctype;
    using Point = typename Policy::Point;
    using Intersection = typename Policy::Intersection;

    //! Deprecated alias, will be removed after 3.1
    using IntersectionType DUNE_DEPRECATED_MSG("Please use Intersection instead") = Intersection;
356
357
358
359
360
361
362

private:
    static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
    using ReferenceElements = typename Dune::ReferenceElements<ctype, dim1>;

public:
    /*!
Dennis Gläser's avatar
Dennis Gläser committed
363
     * \brief Colliding segment and convex polygon
364
365
366
     * \param geo1/geo2 The geometries to intersect
     * \param intersection If the geometries collide intersection holds the
     *        corner points of the intersection object in global coordinates.
367
     * \note This overload is used when segment-like intersections are seeked.
368
     */
369
370
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
371
    {
372
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
373

374
375
        ctype tfirst, tlast;
        if (intersect_(geo1, geo2, tfirst, tlast))
376
        {
377
378
379
380
381
            // the intersection exists. Export the intersections geometry now:
            // s(t) = a + t(b-a) in [tfirst, tlast]
            intersection = Intersection({geo2.global(tfirst), geo2.global(tlast)});
            return true;
        }
382

383
384
        return false;
    }
385

386
387
388
389
390
391
392
393
394
395
396
    /*!
     * \brief Colliding segment and convex polygon
     * \param geo1/geo2 The geometries to intersect
     * \param intersection If the geometries collide intersection holds the
     *        corner points of the intersection object in global coordinates.
     * \note This overload is used when point-like intersections are seeked.
     */
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename P::Intersection& intersection)
    {
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
397

398
399
        ctype tfirst, tlast;
        if (!intersect_(geo1, geo2, tfirst, tlast))
400
        {
401
            // if start & end point are same, the intersection is a point
402
            using std::abs;
403
            if ( tfirst > tlast - eps_ )
404
            {
405
406
                intersection = Intersection(geo2.global(tfirst));
                return true;
407
408
            }
        }
409
410

        return false;
411
    }
412
413
414
415
416
417
418

private:
    /*!
     * \brief Obtain local coordinates of start/end point of the intersecting segment.
     */
     static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast)
     {
419
420
         // lambda to obtain the facet corners on geo1
         auto getFacetCorners = [] (const Geometry1& geo1)
421
         {
422
             std::vector< std::array<int, 2> > facetCorners;
423
424
425
             switch (geo1.corners())
             {
                 case 4: // quadrilateral
426
                     facetCorners = {{0, 2}, {3, 1}, {1, 0}, {2, 3}};
427
428
                     break;
                 case 3: // triangle
429
                     facetCorners = {{1, 0}, {0, 2}, {2, 1}};
430
431
432
433
434
435
                     break;
                 default:
                     DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
                                    << geo1.type() << " with "<< geo1.corners() << " corners.");
             }

436
437
             return facetCorners;
         };
438

439
         // lambda to obtain the normal on a facet
440
         const auto center1 = geo1.center();
441
         auto computeNormal = [&geo1, &center1] (const std::array<int, 2>& facetCorners)
442
         {
443
444
             const auto c0 = geo1.corner(facetCorners[0]);
             const auto c1 = geo1.corner(facetCorners[1]);
445
446
447
448
449
450
451
452
453
454
             const auto edge = c1 - c0;

             Dune::FieldVector<ctype, dimworld> n({-1.0*edge[1], edge[0]});
             n /= n.two_norm();

             // orientation of the element is not known
             // make sure normal points outwards of element
             if ( n*(center1-c0) > 0.0 )
                 n *= -1.0;

455
456
             return n;
         };
457

458
         return Impl::computeSegmentIntersection(geo1, geo2, eps_, tfirst, tlast, getFacetCorners, computeNormal);
459
     }
460
461
};

462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
/*!
 * \ingroup Geometry
 * \brief A class for segment--polygon intersection in 2d space
 */
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 1, 2>
: public GeometryIntersection<Geometry2, Geometry1, Policy, 2, 2, 1>
{
    using Base = GeometryIntersection<Geometry2, Geometry1, Policy, 2, 2, 1>;

public:
    /*!
     * \brief Colliding segment and convex polygon
     * \param geo1/geo2 The geometries to intersect
     * \param intersection If the geometries collide intersection holds the
     *        corner points of the intersection object in global coordinates.
     * \note This forwards to the polygon-segment specialization with swapped arguments.
     */
    template<class P = Policy>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename Base::Intersection& intersection)
    {
        return Base::intersection(geo2, geo1, intersection);
    }
};

487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
/*!
 * \ingroup Geometry
 * \brief A class for polygon--polygon intersection in 2d space
 */
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 2, 2>
{
    enum { dimworld = 2 };
    enum { dim1 = 2 };
    enum { dim2 = 2 };

public:
    using ctype = typename Policy::ctype;
    using Point = typename Policy::Point;
    using Intersection = typename Policy::Intersection;

    //! Deprecated alias, will be removed after 3.1
    using IntersectionType DUNE_DEPRECATED_MSG("Please use Intersection instead") = Intersection;

private:
    static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
    using ReferenceElementsGeo1 = typename Dune::ReferenceElements<ctype, dim1>;
    using ReferenceElementsGeo2 = typename Dune::ReferenceElements<ctype, dim2>;

public:
    /*!
     * \brief Colliding two polygons
     * \note First we find the vertex candidates for the intersection region as follows:
     *       Add polygon vertices that are inside the other polygon
     *       Add intersections of polygon edges
     *       Remove duplicate points from the list
     *       Compute the convex hull polygon
     *       Return a triangulation of that polygon as intersection
     * \param geo1/geo2 The geometries to intersect
     * \param intersection Container to store the corner points of the polygon (as convex hull)
     * \note This overload is used when polygon like intersections are seeked
     */
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 2, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
    {
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");

        // the candidate intersection points
        std::vector<Point> points; points.reserve(6);

        // add polygon1 corners that are inside polygon2
        for (int i = 0; i < geo1.corners(); ++i)
            if (intersectsPointGeometry(geo1.corner(i), geo2))
                points.emplace_back(geo1.corner(i));

        // add polygon2 corners that are inside polygon1
        for (int i = 0; i < geo2.corners(); ++i)
            if (intersectsPointGeometry(geo2.corner(i), geo1))
                points.emplace_back(geo2.corner(i));

        if (points.empty())
            return false;

        const auto referenceElement1 = ReferenceElementsGeo1::general(geo1.type());
        const auto referenceElement2 = ReferenceElementsGeo2::general(geo2.type());

        // add intersections of edges
        using SegGeometry = Dune::MultiLinearGeometry<ctype, 1, dimworld>;
        using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>;
        for (int i = 0; i < referenceElement1.size(dim1-1); ++i)
        {
            const auto localEdgeGeom1 = referenceElement1.template geometry<dim1-1>(i);
            const auto edge1 = SegGeometry( Dune::GeometryTypes::line,
                                            std::vector<Point>( {geo1.global(localEdgeGeom1.corner(0)),
                                                                 geo1.global(localEdgeGeom1.corner(1))} ));

            for (int j = 0; j < referenceElement2.size(dim2-1); ++j)
            {
                const auto localEdgeGeom2 = referenceElement2.template geometry<dim2-1>(j);
                const auto edge2 = SegGeometry( Dune::GeometryTypes::line,
                                                std::vector<Point>( {geo2.global(localEdgeGeom2.corner(0)),
                                                                     geo2.global(localEdgeGeom2.corner(1))} ));

                using EdgeTest = GeometryIntersection<SegGeometry, SegGeometry, PointPolicy>;
                typename EdgeTest::Intersection edgeIntersection;
                if (EdgeTest::intersection(edge1, edge2, edgeIntersection))
                    points.emplace_back(edgeIntersection);
            }
        }

        if (points.empty())
            return false;

        // remove duplicates
        const auto eps = (geo1.corner(0) - geo1.corner(1)).two_norm()*eps_;
        std::sort(points.begin(), points.end(), [&eps](const auto& a, const auto& b) -> bool
        {
            using std::abs;
            return (abs(a[0]-b[0]) > eps ? a[0] < b[0] : a[1] < b[1]);
        });

        auto removeIt = std::unique(points.begin(), points.end(), [&eps](const auto& a, const auto&b)
        {
            return (b-a).two_norm() < eps;
        });

        points.erase(removeIt, points.end());

        // return false if we don't have at least three unique points
        if (points.size() < 3)
            return false;

        // intersection polygon is convex hull of above points
        intersection = grahamConvexHull(points);
        assert(!intersection.empty());
        return true;
    }

    /*!
     * \brief Colliding two polygons
     * \param geo1/geo2 The geometries to intersect
     * \param intersection Container to store the corners of intersection segment
     * \note this overload is used when segment-like intersections are seeked
     * \todo implement this query
     */
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
    {
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
        DUNE_THROW(Dune::NotImplemented, "Polygon-polygon intersection detection for segment-like intersections");
    }

    /*!
     * \brief Colliding two polygons
     * \param geo1/geo2 The geometries to intersect
     * \param intersection The intersection point
     * \note this overload is used when point-like intersections are seeked
     * \todo implement this query
     */
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
    {
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
        DUNE_THROW(Dune::NotImplemented, "Polygon-polygon intersection detection for touching points");
    }
};

629
630
631
632
/*!
 * \ingroup Geometry
 * \brief A class for polyhedron--segment intersection in 3d space
 */
633
634
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 3, 1>
635
{
636
637
638
639
640
    enum { dimworld = 3 };
    enum { dim1 = 3 };
    enum { dim2 = 1 };

public:
641
642
643
644
645
646
    using ctype = typename Policy::ctype;
    using Point = typename Policy::Point;
    using Intersection = typename Policy::Intersection;

    //! Deprecated alias, will be removed after 3.1
    using IntersectionType DUNE_DEPRECATED_MSG("Please use Intersection instead") = Intersection;
647

648
649
650
private:
    static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
    using ReferenceElements = typename Dune::ReferenceElements<ctype, dim1>;
651
652
653
654

public:
    /*!
     *  \brief Colliding segment and convex polyhedron
655
656
     *  \note Algorithm based on the one from "Real-Time Collision Detection" by Christer Ericson,
     *        published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc.
657
658
     *        Basis is the theorem that for any two non-intersecting convex polyhedrons
     *        a separating plane exists.
659
660
     * \param geo1/geo2 The geometries to intersect
     * \param intersection If the geometries collide intersection holds the corner points of
661
     *        the intersection object in global coordinates.
662
     * \note This overload is used when segment-like intersections are seeked.
663
     */
664
665
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
666
    {
667
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
668

669
670
        ctype tfirst, tlast;
        if (intersect_(geo1, geo2, tfirst, tlast))
671
        {
672
673
674
675
676
            // Intersection exists. Export the intersections geometry now:
            // s(t) = a + t(b-a) in [tfirst, tlast]
            intersection = Intersection({geo2.global(tfirst), geo2.global(tlast)});
            return true;
        }
677

678
679
        return false;
    }
680

681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
    /*!
     *  \brief Colliding segment and convex polyhedron
     *  \note Algorithm based on the one from "Real-Time Collision Detection" by Christer Ericson,
     *        published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc.
     *        Basis is the theorem that for any two non-intersecting convex polyhedrons
     *        a separating plane exists.
     * \param geo1/geo2 The geometries to intersect
     * \param intersection If the geometries collide intersection holds the corner points of
     *        the intersection object in global coordinates.
     * \note This overload is used when point-like intersections are seeked.
     */
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
    {
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
696

697
698
699
700
        ctype tfirst, tlast;
        if (intersect_(geo1, geo2, tfirst, tlast))
        {
            // if start & end point are same, the intersection is a point
701
            using std::abs;
702
            if ( abs(tfirst - tlast) < eps_ )
703
            {
704
705
                intersection = Intersection(geo2.global(tfirst));
                return true;
706
707
            }
        }
708
709

        return false;
710
    }
711
712
713
714
715
716
717

private:
    /*!
     * \brief Obtain local coordinates of start/end point of the intersecting segment.
     */
     static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast)
     {
718
719
         // lambda to obtain facet corners on geo1
         auto getFacetCorners = [] (const Geometry1& geo1)
720
         {
721
             std::vector< std::vector<int> > facetCorners;
722
723
724
725
726
             // sort facet corner so that normal n = (p1-p0)x(p2-p0) always points outwards
             switch (geo1.corners())
             {
                 // todo compute intersection geometries
                 case 8: // hexahedron
727
                     facetCorners = {{2, 0, 6, 4}, {1, 3, 5, 7}, {0, 1, 4, 5},
728
729
730
                               {3, 2, 7, 6}, {1, 0, 3, 2}, {4, 5, 6, 7}};
                     break;
                 case 4: // tetrahedron
731
                     facetCorners = {{1, 0, 2}, {0, 1, 3}, {0, 3, 2}, {1, 2, 3}};
732
733
734
735
736
737
                     break;
                 default:
                     DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
                                    << geo1.type() << ", "<< geo1.corners() << " corners.");
             }

738
739
             return facetCorners;
         };
740

741
742
         // lambda to obtain the normal on a facet
         auto computeNormal = [&geo1] (const std::vector<int>& facetCorners)
743
         {
744
745
             const auto v0 = geo1.corner(facetCorners[1]) - geo1.corner(facetCorners[0]);
             const auto v1 = geo1.corner(facetCorners[2]) - geo1.corner(facetCorners[0]);
746
747
748
749

             auto n = crossProduct(v0, v1);
             n /= n.two_norm();

750
751
             return n;
         };
752

753
         return Impl::computeSegmentIntersection(geo1, geo2, eps_, tfirst, tlast, getFacetCorners, computeNormal);
754
     }
755
756
};

757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
/*!
 * \ingroup Geometry
 * \brief A class for segment--polyhedron intersection in 3d space
 */
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 1, 3>
: public GeometryIntersection<Geometry2, Geometry1, Policy, 3, 3, 1>
{
    using Base = GeometryIntersection<Geometry2, Geometry1, Policy, 3, 3, 1>;
public:
    /*!
     * \brief Colliding segment and convex polyhedron
     * \param geo1/geo2 The geometries to intersect
     * \param intersection If the geometries collide intersection holds the
     *        corner points of the intersection object in global coordinates.
     * \note This forwards to the polyhedron-segment specialization with swapped arguments.
     */
    template<class P = Policy>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename Base::Intersection& intersection)
    {
        return Base::intersection(geo2, geo1, intersection);
    }
};

781
782
783
784
/*!
 * \ingroup Geometry
 * \brief A class for polyhedron--polygon intersection in 3d space
 */
785
786
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 3, 2>
787
788
789
790
791
792
{
    enum { dimworld = 3 };
    enum { dim1 = 3 };
    enum { dim2 = 2 };

public:
793
794
795
796
797
798
    using ctype = typename Policy::ctype;
    using Point = typename Policy::Point;
    using Intersection = typename Policy::Intersection;

    //! Deprecated alias, will be removed after 3.1
    using IntersectionType DUNE_DEPRECATED_MSG("Please use Intersection instead") = Intersection;
799
800
801
802
803
804
805
806

private:
    static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
    using ReferenceElementsGeo1 = typename Dune::ReferenceElements<ctype, dim1>;
    using ReferenceElementsGeo2 = typename Dune::ReferenceElements<ctype, dim2>;

public:
    /*!
807
     * \brief Colliding polygon and convex polyhedron
808
809
810
811
812
813
814
815
816
     * \note First we find the vertex candidates for the intersection region as follows:
     *       Add triangle vertices that are inside the tetrahedron
     *       Add tetrahedron vertices that are inside the triangle
     *       Add all intersection points of tetrahedron edges (codim 2) with the triangle (codim 0) (6*1 tests)
     *       Add all intersection points of triangle edges (codim 1) with tetrahedron faces (codim 1) (3*4 tests)
     *       Remove duplicate points from the list
     *       Compute the convex hull polygon
     *       Return a triangulation of that polygon as intersection
     * \param geo1/geo2 The geometries to intersect
817
     * \param intersection Container to store the corner points of the polygon (as convex hull)
818
     * \note This overload is used when polygon like intersections are seeked
819
     */
820
821
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 2, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
822
    {
823
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844

        // the candidate intersection points
        std::vector<Point> points; points.reserve(10);

        // add 3d geometry corners that are inside the 2d geometry
        for (int i = 0; i < geo1.corners(); ++i)
            if (intersectsPointGeometry(geo1.corner(i), geo2))
                points.emplace_back(geo1.corner(i));

        // add 2d geometry corners that are inside the 3d geometry
        for (int i = 0; i < geo2.corners(); ++i)
            if (intersectsPointGeometry(geo2.corner(i), geo1))
                points.emplace_back(geo2.corner(i));

        // get some geometry types
        using PolyhedronFaceGeometry = Dune::MultiLinearGeometry<ctype, 2, dimworld>;
        using SegGeometry = Dune::MultiLinearGeometry<ctype, 1, dimworld>;

        const auto referenceElement1 = ReferenceElementsGeo1::general(geo1.type());
        const auto referenceElement2 = ReferenceElementsGeo2::general(geo2.type());

845
846
847
        // intersection policy for point-like intersections (used later)
        using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>;

848
849
850
851
852
853
854
855
        // add intersection points of all polyhedron edges (codim dim-1) with the polygon
        for (int i = 0; i < referenceElement1.size(dim1-1); ++i)
        {
            const auto localEdgeGeom = referenceElement1.template geometry<dim1-1>(i);
            const auto p = geo1.global(localEdgeGeom.corner(0));
            const auto q = geo1.global(localEdgeGeom.corner(1));
            const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q});

856
857
858
859
            using PolySegTest = GeometryIntersection<Geometry2, SegGeometry, PointPolicy>;
            typename PolySegTest::Intersection polySegIntersection;
            if (PolySegTest::intersection(geo2, segGeo, polySegIntersection))
                points.emplace_back(polySegIntersection);
860
861
862
863
864
865
866
867
868
869
870
871
872
873
        }

        // add intersection points of all polygon faces (codim 1) with the polyhedron faces
        for (int i = 0; i < referenceElement1.size(1); ++i)
        {
            const auto faceGeo = [&]()
            {
                const auto localFaceGeo = referenceElement1.template geometry<1>(i);
                if (localFaceGeo.corners() == 4)
                {
                    const auto a = geo1.global(localFaceGeo.corner(0));
                    const auto b = geo1.global(localFaceGeo.corner(1));
                    const auto c = geo1.global(localFaceGeo.corner(2));
                    const auto d = geo1.global(localFaceGeo.corner(3));
874

875
876
877
878
879
880
881
                    return PolyhedronFaceGeometry(Dune::GeometryTypes::cube(2), std::vector<Point>{a, b, c, d});
                }
                else
                {
                    const auto a = geo1.global(localFaceGeo.corner(0));
                    const auto b = geo1.global(localFaceGeo.corner(1));
                    const auto c = geo1.global(localFaceGeo.corner(2));
882

883
884
885
886
887
888
889
890
891
                    return PolyhedronFaceGeometry(Dune::GeometryTypes::simplex(2), std::vector<Point>{a, b, c});
                }
            }();

            for (int j = 0; j < referenceElement2.size(1); ++j)
            {
                const auto localEdgeGeom = referenceElement2.template geometry<1>(j);
                const auto p = geo2.global(localEdgeGeom.corner(0));
                const auto q = geo2.global(localEdgeGeom.corner(1));
892

893
894
                const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q});

895
896
897
898
                using PolySegTest = GeometryIntersection<PolyhedronFaceGeometry, SegGeometry, PointPolicy>;
                typename PolySegTest::Intersection polySegIntersection;
                if (PolySegTest::intersection(faceGeo, segGeo, polySegIntersection))
                    points.emplace_back(polySegIntersection);
899
900
901
902
903
904
905
906
907
908
            }
        }

        // return if no intersection points were found
        if (points.empty()) return false;

        // remove duplicates
        const auto eps = (geo1.corner(0) - geo1.corner(1)).two_norm()*eps_;
        std::sort(points.begin(), points.end(), [&eps](const auto& a, const auto& b) -> bool
        {
909
            using std::abs;
910
911
912
913
914
915
916
917
918
919
920
921
922
            return (abs(a[0]-b[0]) > eps ? a[0] < b[0] : (abs(a[1]-b[1]) > eps ? a[1] < b[1] : (a[2] < b[2])));
        });

        auto removeIt = std::unique(points.begin(), points.end(), [&eps](const auto& a, const auto&b)
        {
            return (b-a).two_norm() < eps;
        });

        points.erase(removeIt, points.end());

        // return false if we don't have more than three unique points
        if (points.size() < 3) return false;

923
924
925
        // intersection polygon is convex hull of above points
        intersection = grahamConvexHull2d3d(points);
        assert(!intersection.empty());
926
927
928

        return true;
    }
929
930
931
932
933
934
935
936
937
938
939
940

    /*!
     * \brief Colliding segment and convex polyhedron
     * \note First we find the vertex candidates for the intersection region as follows:
     *       Add triangle vertices that are inside the tetrahedron
     *       Add tetrahedron vertices that are inside the triangle
     *       Add all intersection points of tetrahedron edges (codim 2) with the triangle (codim 0) (6*1 tests)
     *       Add all intersection points of triangle edges (codim 1) with tetrahedron faces (codim 1) (3*4 tests)
     *       Remove duplicate points from the list
     *       Compute the convex hull polygon
     *       Return a triangulation of that polygon as intersection
     * \param geo1/geo2 The geometries to intersect
941
     * \param intersection Container to store the intersection result
942
943
944
945
946
947
948
949
     * \todo implement overloads for segment or point-like intersections
     */
    template<class P = Policy, std::enable_if_t<P::dimIntersection != 2, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
    {
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
        DUNE_THROW(Dune::NotImplemented, "Polyhedron-polygon intersection detection only implemented for polygon-like intersections");
    }
950
951
};

952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
/*!
 * \ingroup Geometry
 * \brief A class for polygon--polyhedron intersection in 3d space
 */
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 2, 3>
: public GeometryIntersection<Geometry2, Geometry1, Policy, 3, 3, 2>
{
    using Base = GeometryIntersection<Geometry2, Geometry1, Policy, 3, 3, 2>;
public:
    /*!
     * \brief Colliding polygon and convex polyhedron
     * \param geo1/geo2 The geometries to intersect
     * \param intersection If the geometries collide intersection holds the
     *        corner points of the intersection object in global coordinates.
     * \note This forwards to the polyhedron-polygon specialization with swapped arguments.
     */
    template<class P = Policy>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename Base::Intersection& intersection)
    {
        return Base::intersection(geo2, geo1, intersection);
    }
};

976
977
978
979
/*!
 * \ingroup Geometry
 * \brief A class for polygon--segment intersection in 3d space
 */
980
981
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 2, 1>
982
983
984
985
986
987
{
    enum { dimworld = 3 };
    enum { dim1 = 2 };
    enum { dim2 = 1 };

public:
988
989
990
991
992
993
    using ctype = typename Policy::ctype;
    using Point = typename Policy::Point;
    using Intersection = typename Policy::Intersection;

    //! Deprecated alias, will be removed after 3.1
    using IntersectionType DUNE_DEPRECATED_MSG("Please use Intersection instead") = Intersection;
994
995
996
997
998
999

private:
    static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
    using ReferenceElements = typename Dune::ReferenceElements<ctype, dim1>;

public:
1000
    /*!