geometryintersection.hh 53.5 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
27
#include <tuple>

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

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

38
namespace Dumux {
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
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;
};

79
80
81
82
83
84
85
86
87
88
89
90
91
//! Policy structure for polyhedron-like intersections
template<class ct, int dw>
struct PolyhedronPolicy
{
    using ctype = ct;
    using Point = Dune::FieldVector<ctype, dw>;
    using Polyhedron = std::vector<Point>;

    static constexpr int dimWorld = dw;
    static constexpr int dimIntersection = 3;
    using Intersection = Polyhedron;
};

92
//! default policy chooser class
93
template<class Geometry1, class Geometry2>
94
95
96
class DefaultPolicyChooser
{
    static constexpr int dimworld = Geometry1::coorddimension;
97
    static constexpr int isDim = std::min( int(Geometry1::mydimension), int(Geometry2::mydimension) );
98
99
100
101
    static_assert(dimworld == int(Geometry2::coorddimension),
                  "Geometries must have the same coordinate dimension!");
    static_assert(int(Geometry1::mydimension) <= 3 && int(Geometry2::mydimension) <= 3,
                  "Geometries must have dimension 3 or less.");
102
    using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
103
104
105
106
107

    using DefaultPolicies = std::tuple<PointPolicy<ctype, dimworld>,
                                       SegmentPolicy<ctype, dimworld>,
                                       PolygonPolicy<ctype, dimworld>,
                                       PolyhedronPolicy<ctype, dimworld>>;
108
public:
109
    using type = std::tuple_element_t<isDim, DefaultPolicies>;
110
111
112
113
114
115
116
};

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

} // end namespace IntersectionPolicy
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
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
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

204
/*!
205
 * \ingroup Geometry
206
207
208
209
 * \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.
 */
210
211
template
<class Geometry1, class Geometry2,
212
213
214
215
 class Policy = IntersectionPolicy::DefaultPolicy<Geometry1, Geometry2>,
 int dimworld = Geometry1::coorddimension,
 int dim1 = Geometry1::mydimension,
 int dim2 = Geometry2::mydimension>
216
class GeometryIntersection
217
{
218
219
    static constexpr int dimWorld = Policy::dimWorld;

220
public:
221
222
223
224
225
226
    using ctype = typename Policy::ctype;
    using Point = typename Policy::Point;
    using Intersection = typename Policy::Intersection;

    //! Determine if the two geometries intersect and compute the intersection geometry
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
227
    {
228
229
        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 = "
230
231
232
233
                                          << dimworld << ", dim1 = " << dim1 << ", dim2 = " << dim2);
    }
};

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
338
339
340
341
342
343
344
345
346
347
348
349
350
/*!
 * \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;

    /*!
     * \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");
    }
};

351
352
/*!
 * \ingroup Geometry
Dennis Gläser's avatar
Dennis Gläser committed
353
 * \brief A class for polygon--segment intersection in 2d space
354
 */
355
356
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 2, 1>
357
358
359
360
361
362
{
    enum { dimworld = 2 };
    enum { dim1 = 2 };
    enum { dim2 = 1 };

public:
363
364
365
366
    using ctype = typename Policy::ctype;
    using Point = typename Policy::Point;
    using Intersection = typename Policy::Intersection;

367
368
369
370
371
372
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
373
     * \brief Colliding segment and convex polygon
374
375
376
     * \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.
377
     * \note This overload is used when segment-like intersections are seeked.
378
     */
379
380
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
381
    {
382
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
383

384
385
        ctype tfirst, tlast;
        if (intersect_(geo1, geo2, tfirst, tlast))
386
        {
387
388
389
390
391
            // 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;
        }
392

393
394
        return false;
    }
395

396
397
398
399
400
401
402
403
404
405
406
    /*!
     * \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");
407

408
409
        ctype tfirst, tlast;
        if (!intersect_(geo1, geo2, tfirst, tlast))
410
        {
411
            // if start & end point are same, the intersection is a point
412
            using std::abs;
413
            if ( tfirst > tlast - eps_ )
414
            {
415
416
                intersection = Intersection(geo2.global(tfirst));
                return true;
417
418
            }
        }
419
420

        return false;
421
    }
422
423
424
425
426
427
428

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)
     {
429
430
         // lambda to obtain the facet corners on geo1
         auto getFacetCorners = [] (const Geometry1& geo1)
431
         {
432
             std::vector< std::array<int, 2> > facetCorners;
433
434
435
             switch (geo1.corners())
             {
                 case 4: // quadrilateral
436
                     facetCorners = {{0, 2}, {3, 1}, {1, 0}, {2, 3}};
437
438
                     break;
                 case 3: // triangle
439
                     facetCorners = {{1, 0}, {0, 2}, {2, 1}};
440
441
442
443
444
445
                     break;
                 default:
                     DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
                                    << geo1.type() << " with "<< geo1.corners() << " corners.");
             }

446
447
             return facetCorners;
         };
448

449
         // lambda to obtain the normal on a facet
450
         const auto center1 = geo1.center();
451
         auto computeNormal = [&geo1, &center1] (const std::array<int, 2>& facetCorners)
452
         {
453
454
             const auto c0 = geo1.corner(facetCorners[0]);
             const auto c1 = geo1.corner(facetCorners[1]);
455
456
457
458
459
460
461
462
463
464
             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;

465
466
             return n;
         };
467

468
         return Impl::computeSegmentIntersection(geo1, geo2, eps_, tfirst, tlast, getFacetCorners, computeNormal);
469
     }
470
471
};

472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
/*!
 * \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);
    }
};

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
/*!
 * \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;

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));

544
545
        const auto numPoints1 = points.size();
        if (numPoints1 != geo1.corners())
546
        {
547
548
549
550
551
552
553
            // 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;
554

555
            if (points.size() - numPoints1 != geo2.corners())
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
                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);
                    }
                }
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
            }
        }

        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
609
        intersection = grahamConvexHull<2>(points);
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
        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");
    }
};

643
644
645
646
/*!
 * \ingroup Geometry
 * \brief A class for polyhedron--segment intersection in 3d space
 */
647
648
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 3, 1>
649
{
650
651
652
653
654
    enum { dimworld = 3 };
    enum { dim1 = 3 };
    enum { dim2 = 1 };

public:
655
656
657
658
    using ctype = typename Policy::ctype;
    using Point = typename Policy::Point;
    using Intersection = typename Policy::Intersection;

659
660
661
private:
    static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
    using ReferenceElements = typename Dune::ReferenceElements<ctype, dim1>;
662
663
664
665

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

680
681
        ctype tfirst, tlast;
        if (intersect_(geo1, geo2, tfirst, tlast))
682
        {
683
684
685
686
687
            // 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;
        }
688

689
690
        return false;
    }
691

692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
    /*!
     *  \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");
707

708
709
710
711
        ctype tfirst, tlast;
        if (intersect_(geo1, geo2, tfirst, tlast))
        {
            // if start & end point are same, the intersection is a point
712
            using std::abs;
713
            if ( abs(tfirst - tlast) < eps_ )
714
            {
715
716
                intersection = Intersection(geo2.global(tfirst));
                return true;
717
718
            }
        }
719
720

        return false;
721
    }
722
723
724
725
726
727
728

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)
     {
729
730
         // lambda to obtain facet corners on geo1
         auto getFacetCorners = [] (const Geometry1& geo1)
731
         {
732
             std::vector< std::vector<int> > facetCorners;
733
734
735
736
737
             // 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
738
                     facetCorners = {{2, 0, 6, 4}, {1, 3, 5, 7}, {0, 1, 4, 5},
739
740
741
                               {3, 2, 7, 6}, {1, 0, 3, 2}, {4, 5, 6, 7}};
                     break;
                 case 4: // tetrahedron
742
                     facetCorners = {{1, 0, 2}, {0, 1, 3}, {0, 3, 2}, {1, 2, 3}};
743
744
745
746
747
748
                     break;
                 default:
                     DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
                                    << geo1.type() << ", "<< geo1.corners() << " corners.");
             }

749
750
             return facetCorners;
         };
751

752
753
         // lambda to obtain the normal on a facet
         auto computeNormal = [&geo1] (const std::vector<int>& facetCorners)
754
         {
755
756
             const auto v0 = geo1.corner(facetCorners[1]) - geo1.corner(facetCorners[0]);
             const auto v1 = geo1.corner(facetCorners[2]) - geo1.corner(facetCorners[0]);
757
758
759
760

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

761
762
             return n;
         };
763

764
         return Impl::computeSegmentIntersection(geo1, geo2, eps_, tfirst, tlast, getFacetCorners, computeNormal);
765
     }
766
767
};

768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
/*!
 * \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);
    }
};

792
793
794
795
/*!
 * \ingroup Geometry
 * \brief A class for polyhedron--polygon intersection in 3d space
 */
796
797
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 3, 2>
798
799
800
801
802
803
{
    enum { dimworld = 3 };
    enum { dim1 = 3 };
    enum { dim2 = 2 };

public:
804
805
806
807
    using ctype = typename Policy::ctype;
    using Point = typename Policy::Point;
    using Intersection = typename Policy::Intersection;

808
809
810
811
812
813
814
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:
    /*!
815
     * \brief Colliding polygon and convex polyhedron
816
817
818
819
820
821
822
823
824
     * \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
825
     * \param intersection Container to store the corner points of the polygon (as convex hull)
826
     * \note This overload is used when polygon like intersections are seeked
827
     */
828
829
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 2, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
830
    {
831
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852

        // 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());

853
854
855
        // intersection policy for point-like intersections (used later)
        using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>;

856
857
858
859
860
861
862
863
        // 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});

864
865
866
867
            using PolySegTest = GeometryIntersection<Geometry2, SegGeometry, PointPolicy>;
            typename PolySegTest::Intersection polySegIntersection;
            if (PolySegTest::intersection(geo2, segGeo, polySegIntersection))
                points.emplace_back(polySegIntersection);
868
869
870
871
872
873
874
875
876
877
878
879
880
881
        }

        // 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));
882

883
884
885
886
887
888
889
                    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));
890

891
892
893
894
895
896
897
898
899
                    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));
900

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

903
904
905
906
                using PolySegTest = GeometryIntersection<PolyhedronFaceGeometry, SegGeometry, PointPolicy>;
                typename PolySegTest::Intersection polySegIntersection;
                if (PolySegTest::intersection(faceGeo, segGeo, polySegIntersection))
                    points.emplace_back(polySegIntersection);
907
908
909
910
911
912
913
914
915
916
            }
        }

        // 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
        {
917
            using std::abs;
918
919
920
921
922
923
924
925
926
927
928
929
930
            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;

931
        // intersection polygon is convex hull of above points
932
        intersection = grahamConvexHull<2>(points);
933
        assert(!intersection.empty());
934
935
936

        return true;
    }
937
938
939
940
941
942
943
944
945
946
947
948

    /*!
     * \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
949
     * \param intersection Container to store the intersection result
950
951
952
953
954
955
956
957
     * \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");
    }
958
959
};

960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
/*!
 * \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);
    }
};

984
985
986
987
/*!
 * \ingroup Geometry
 * \brief A class for polygon--segment intersection in 3d space
 */
988
989
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 2, 1>
990
991
992
993
994
995
{
    enum { dimworld = 3 };
    enum { dim1 = 2 };
    enum { dim2 = 1 };

public:
996
997
998
999
    using ctype = typename Policy::ctype;
    using Point = typename Policy::Point;
    using Intersection = typename Policy::Intersection;

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

public:
1005
1006
1007
1008
1009
    /*!
     *  \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. (Chapter 5.3.6)
     * \param geo1/geo2 The geometries to intersect
1010
     * \param intersection If the geometries collide, is holds the corner points of
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
     *        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 == 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");

        ctype tfirst, tlast;
        if (intersect_(geo1, geo2, tfirst, tlast))
        {
            // 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;
        }

        return false;
    }

1031
1032
1033
1034
1035
    /*!
     *  \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. (Chapter 5.3.6)
     * \param geo1/geo2 The geometries to intersect
1036
     * \param is If the geometries collide, is holds the corner points of
1037
     *        the intersection object in global coordinates.
1038
     * \note This overload is used when point-like intersections are seeked
1039
     */
1040
1041
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& is)
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
    {
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");

        const auto p = geo2.corner(0);
        const auto q = geo2.corner(1);

        const auto a = geo1.corner(0);
        const auto b = geo1.corner(1);
        const auto c = geo1.corner(2);

        if (geo1.corners() == 3)
1053
            return intersect_<Policy>(a, b, c, p, q, is);
1054
1055
1056

        else if (geo1.corners() == 4)
        {
1057
1058
1059
            Intersection is1, is2;
            bool hasSegment1, hasSegment2;

1060
            const auto d = geo1.corner(3);
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
            const bool intersects1 = intersect_<Policy>(a, b, d, p, q, is1, hasSegment1);
            const bool intersects2 = intersect_<Policy>(a, d, c, p, q, is2, hasSegment2);

            if (hasSegment1 || hasSegment2)
                return false;

            if (intersects1) { is = is1; return true; }
            if (intersects2) { is = is2; return true; }

            return false;
1071
1072
1073
1074
1075
1076
1077
        }

        else
            DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
                          << geo1.type() << ", "<< geo1.corners() << " corners.");
    }

1078
private:
1079
1080
1081
1082
1083
1084
    /*!
     * \brief Obtain local coordinates of start/end point of the intersecting segment.
     */
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
    static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast)
    {
1085
1086
        // lambda to obtain the facet corners on geo1
        auto getFacetCorners = [] (const Geometry1& geo1)
1087
        {
1088
            std::vector< std::array<int, 2> > facetCorners;
1089
1090
1091
            switch (geo1.corners())
            {
                case 4: // quadrilateral
1092
                    facetCorners = {{0, 2}, {3, 1}, {1, 0}, {2, 3}};
1093
1094
                    break;
                case 3: // triangle
1095
                    facetCorners = {{1, 0}, {0, 2}, {2, 1}};
1096
1097
1098
1099
1100
1101
                    break;
                default:
                    DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
                                   << geo1.type() << " with "<< geo1.corners() << " corners.");
            }

1102
1103
            return facetCorners;
        };
1104
1105
1106
1107
1108

        const auto center1 = geo1.center();
        const auto normal1 = crossProduct(geo1.corner(1) - geo1.corner(0),
                                          geo1.corner(2) - geo1.corner(0));

1109
1110
        // lambda to obtain the normal on a facet
        auto computeNormal = [&center1, &normal1, &geo1] (const std::array<int, 2>& facetCorners)
1111
        {
1112
1113
            const auto c0 = geo1.corner(facetCorners[0]);
            const auto c1 = geo1.corner(facetCorners[1]);
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
            const auto edge = c1 - c0;

            Dune::FieldVector<ctype, dimworld> n = crossProduct(edge, normal1);
            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;

1124
1125
            return n;
        };
1126

1127
        return Impl::computeSegmentIntersection(geo1, geo2, eps_, tfirst, tlast, getFacetCorners, computeNormal);
1128
1129
1130
1131
1132
    }

    /*!
     * \brief triangle--segment point-like intersection with points as input.
     */
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
    static bool intersect_(const Point& a, const Point& b, const Point& c,
                           const Point& p, const Point& q,
                           Intersection& is)
    {
        bool hasSegment;
        return intersect_(a, b, c, p, q, is, hasSegment);
    }

    /*!
     * \brief triangle--segment point-like intersection with points as input. Also
     *        stores if a segment-like intersection was found in the provided boolean.
     */
    template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
    static bool intersect_(const Point& a, const Point& b, const Point& c,
                           const Point& p, const Point& q,
                           Intersection& is, bool& hasSegmentIs)
    {
        hasSegmentIs = false;
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162

        const auto ab = b - a;
        const auto ac = c - a;
        const auto qp = p - q;

        // compute the triangle normal that defines the triangle plane
        const auto normal = crossProduct(ab, ac);

        // compute the denominator
        // if denom is 0 the segment is parallel and we can return
        const auto denom = normal*qp;
1163
1164
        const auto abnorm2 = ab.two_norm2();
        const auto eps = eps_*abnorm2*qp.two_norm();
1165
1166
        using std::abs;
        if (abs(denom) < eps)
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
        {
            const auto pa = a - p;
            const auto denom2 = normal*pa;
            if (abs(denom2) > eps_*pa.two_norm()*abnorm2)
                return false;

            // if we get here, we are in-plane. Check if a
            // segment-like intersection with geometry 1 exists.
            using SegmentPolicy = typename IntersectionPolicy::SegmentPolicy<ctype, dimworld>;
            using Triangle = Dune::AffineGeometry<ctype, 2, dimworld>;
            using Segment = Dune::AffineGeometry<ctype, 1, dimworld>;
            using SegmentIntersectionAlgorithm = GeometryIntersection<Triangle, Segment, SegmentPolicy>;
            using SegmentIntersectionType = typename SegmentIntersectionAlgorithm::Intersection;
            SegmentIntersectionType segmentIs;

            Triangle triangle(Dune::GeometryTypes::simplex(2), std::array<Point, 3>({a, b, c}));
            Segment segment(Dune::GeometryTypes::simplex(1), std::array<Point, 2>({p, q}));
            if (SegmentIntersectionAlgorithm::intersection(triangle, segment, segmentIs))
            {
                hasSegmentIs = true;
                return false;
            }

            // We are in-plane. A point-like
            // intersection can only be on the edges
            for (const auto& ip : {p, q})
            {
                if ( intersectsPointSimplex(ip, a, b)
                     || intersectsPointSimplex(ip, b, c)
                     || intersectsPointSimplex(ip, c, a) )
                {
                    is = ip;
                    return true;
                }
            }

1203
            return false;
1204
        }
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228

        // compute intersection t value of pq with plane of triangle.
        // a segment intersects if and only if 0 <= t <= 1.
        const auto ap = p - a;
        const auto t = (ap*normal)/denom;
        if (t < 0.0 - eps_) return false;
        if (t > 1.0 + eps_) return false;

        // compute the barycentric coordinates and check if the intersection point
        // is inside the bounds of the triangle
        const auto e = crossProduct(qp, ap);
        const auto v = (ac*e)/denom;
        if (v < -eps_ || v > 1.0 + eps_) return false;
        const auto w = -(ab*e)/denom;
        if (w < -eps_ || v + w > 1.0 + eps_) return false;

        // Now we are sure there is an intersection points
        // Perform delayed division compute the last barycentric coordinate component
        const auto u = 1.0 - v - w;

        Point ip(0.0);
        ip.axpy(u, a);
        ip.axpy(v, b);
        ip.axpy(w, c);
1229
        is = ip;
1230
        return true;
1231
1232
1233
    }
};

1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
/*!
 * \ingroup Geometry
 * \brief A class for segment--polygon intersection in 3d space
 */
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 1, 2>
: public GeometryIntersection<Geometry2, Geometry1, Policy, 3, 2, 1>
{
    using Base = GeometryIntersection<Geometry2, Geometry1, Policy, 3, 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 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);
    }
};

1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
/*!
 * \ingroup Geometry
 * \brief A class for segment--segment intersection in 3d space
 */
template <class Geometry1, class Geometry2, class Policy>
class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 1, 1>
{
    enum { dimworld = 3 };
    enum { dim1 = 1 };
    enum { dim2 = 1 };

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

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

public:
    /*!
     * \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 == 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, "segment-segment intersection detection for point-like intersections");
    }

    /*!
     *  \brief Colliding two segments in 3D
     * \param geo1/geo2 The geometries to intersect
Melanie Lipp's avatar
Melanie Lipp committed
1296
     * \param intersection If the geometries collide, is holds the corner points of
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
     *        the intersection object in global coordinates.
     * \note This overload is used when segment-like intersections are seeked
     */
    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");

        const auto& a = geo1.corner(0);
        const auto& b = geo1.corner(1);
        const auto& p = geo2.corner(0);
        const auto& q = geo2.corner(1);

        const auto ab = b-a;
        const auto pq = q-p;
        const auto abNorm = ab.two_norm();
        const auto pqNorm = pq.two_norm();

        using std::max;
        const auto maxNorm = max(abNorm, pqNorm);
        const auto eps2 = eps_*maxNorm*maxNorm;

        // if the segment are not parallel there is no segment intersection
        using std::abs;
        if (!(abs(abs(ab*pq) - abNorm*pqNorm) < eps2))
            return false;

        const auto ap = (p-a);
        const auto aq = (q-a);

        // points have to be colinear
        if (!(crossProduct(ap, aq).two_norm() < eps2))
            return false;

        // scaled local coordinates
        // we save the division for the normalization for now
        // and do it in the very end if we find an intersection
        auto tp = ap*ab;
        auto tq = aq*ab;
        const auto abnorm2 = abNorm*abNorm;

        // make sure they are sorted
        using std::swap;
        if (tp > tq)
            swap(tp, tq);

        using std::min; using std::max;
1344
1345
        tp = min(abnorm2, max(0.0, tp));
        tq = max(0.0, min(abnorm2, tq));
1346

1347
        if (abs(tp-tq) < eps2)
1348
1349
1350
1351
1352
1353
1354
            return false;

        intersection = Intersection({geo1.global(tp/abnorm2), geo1.global(tq/abnorm2)});
        return true;
    }
};

1355
1356
} // end namespace Dumux

1357
# endif