geometryintersection.hh 17.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/*****************************************************************************
 *   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    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (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 Common
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/exceptions.hh>
27
#include <dune/common/promotiontraits.hh>
28
#include <dune/geometry/referenceelements.hh>
29
#include <dune/geometry/multilineargeometry.hh>
30
31

#include <dumux/common/math.hh>
32
33
#include <dumux/common/geometry/intersectspointgeometry.hh>
#include <dumux/common/geometry/grahamconvexhull.hh>
34

35
namespace Dumux {
36

37
/*!
38
 * \ingroup Common
39
40
41
42
 * \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.
 */
43
44
45
46
47
template
<class Geometry1, class Geometry2,
  int dimworld = Geometry1::coorddimension,
  int dim1 = Geometry1::mydimension,
  int dim2 = Geometry2::mydimension>
48
class GeometryIntersection
49
50
{
public:
51
    //! Determine if the two geometries intersect and compute the intersection corners
52
    template<class IntersectionType>
53
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& intersection)
54
    {
55
56
        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 = "
57
58
59
60
                                          << dimworld << ", dim1 = " << dim1 << ", dim2 = " << dim2);
    }
};

61
//! polyhedron--segment intersection in 3d space
62
template <class Geometry1, class Geometry2>
63
class GeometryIntersection<Geometry1, Geometry2, 3, 3, 1>
64
{
65
66
67
68
69
70
    enum { dimworld = 3 };
    enum { dim1 = 3 };
    enum { dim2 = 1 };

public:
    using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
71
72
    using Point = Dune::FieldVector<ctype, dimworld>;
    using IntersectionType = std::array<std::vector<Point>, 1>;
73

74
75
76
private:
    static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
    using ReferenceElements = typename Dune::ReferenceElements<ctype, dim1>;
77
78
79
80

public:
    /*!
     *  \brief Colliding segment and convex polyhedron
81
82
     *  \note Algorithm based on the one from "Real-Time Collision Detection" by Christer Ericson,
     *        published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc.
83
84
     *        Basis is the theorem that for any two non-intersecting convex polyhedrons
     *        a separating plane exists.
85
86
     * \param geo1/geo2 The geometries to intersect
     * \param intersection If the geometries collide intersection holds the corner points of
87
     *        the intersection object in global coordinates.
88
     */
89
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& intersection)
90
    {
91
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
92
93
94
95
96
97
98
99

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

        // The initial interval is the whole segment
        // afterward we start clipping the interval
        // by the planes decribed by the facet
100
101
        ctype tfirst = 0.0;
        ctype tlast = 1.0;
102

103
        const std::vector<std::vector<int>> facets = [&]()
104
        {
105
106
107
108
109
110
111
112
113
114
115
116
117
            std::vector<std::vector<int>> facets;
            // 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
                    facets = {{2, 0, 6, 4}, {1, 3, 5, 7}, {0, 1, 4, 5},
                              {3, 2, 7, 6}, {1, 0, 3, 2}, {4, 5, 6, 7}};
                    break;
                case 4: // tetrahedron
                    facets = {{1, 2, 0}, {0, 1, 3}, {0, 3, 2}, {1, 2, 3}};
                    break;
                default:
118
119
                    DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
                                   << geo1.type() << ", "<< geo1.corners() << " corners.");
120
121
122
123
            }

            return facets;
        }();
124
125
126
127
128
129
130
131

        for (const auto& f : facets)
        {
            // compute normal vector by cross product
            const auto v0 = geo1.corner(f[1]) - geo1.corner(f[0]);
            const auto v1 = geo1.corner(f[2]) - geo1.corner(f[0]);
            const auto eps = eps_*v0.two_norm();

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

135
136
            const ctype denom = n*d;
            const ctype dist = n*(a-geo1.corner(f[0]));
137
138
139

            // if denominator is zero the segment in parallel to
            // the plane. If the distance is positive there is no intersection
140
141
            using std::abs;
            if (abs(denom) < eps)
142
143
144
145
146
            {
                if (dist > eps)
                    return false;
            }
            else // not parallel: compute line-plane intersection
147
            {
148
                const ctype t = -dist / denom;
149
                // if entering half space cut tfirst if t is larger
150
151
                using std::signbit;
                if (signbit(denom))
152
                {
153
154
                    if (t > tfirst)
                        tfirst = t;
155
                }
156
157
158
159
160
161
162
163
164
165
166
                // 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 - eps_)
                    return false;

167
168
            }
        }
169
170
        // If we made it until here an intersection exists. We also export
        // the intersections geometry now s(t) = a + t(b-a) in [tfirst, tlast]
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
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
        intersection = {std::vector<Point>({geo2.global(tfirst), geo2.global(tlast)})};
        return true;
    }
};

//! polyhedron--polygon intersection in 3d space
template <class Geometry1, class Geometry2>
class GeometryIntersection<Geometry1, Geometry2, 3, 3, 2>
{
    enum { dimworld = 3 };
    enum { dim1 = 3 };
    enum { dim2 = 2 };

public:
    using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
    using Point = Dune::FieldVector<ctype, dimworld>;
    using IntersectionType = std::vector<std::array<Point, 3>>;

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 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
     * \param intersection A triangulation of the intersection polygon
     */
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& 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(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());

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

            using PolySegTest = GeometryIntersection<Geometry2, SegGeometry>;
242
243
244
            typename PolySegTest::IntersectionType polySegIntersection;
            if (PolySegTest::template intersection<2>(geo2, segGeo, polySegIntersection))
                points.emplace_back(polySegIntersection[0]);
245
246
247
248
249
250
251
252
253
254
255
256
257
258
        }

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

260
261
262
263
264
265
266
                    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));
267

268
269
270
271
272
273
274
275
276
                    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));
277

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

                using PolySegTest = GeometryIntersection<PolyhedronFaceGeometry, SegGeometry>;
281
282
283
                typename PolySegTest::IntersectionType polySegIntersection;
                if (PolySegTest::template intersection<2>(faceGeo, segGeo, polySegIntersection))
                    points.emplace_back(polySegIntersection[0]);
284
285
286
287
288
289
290
291
292
293
            }
        }

        // 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
        {
294
            using std::abs;
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
            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;

        // compute convex hull
        const auto convexHull = grahamConvexHull2d3d(points);
310
        assert(!convexHull.empty());
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

        // the intersections are the triangulation of the convex hull polygon
        intersection = triangulateConvexHull(convexHull);
        return true;
    }
};

//! polygon--segment intersection in 3d space
template <class Geometry1, class Geometry2>
class GeometryIntersection<Geometry1, Geometry2, 3, 2, 1>
{
    enum { dimworld = 3 };
    enum { dim1 = 2 };
    enum { dim2 = 1 };

public:
    using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
    using Point = Dune::FieldVector<ctype, dimworld>;
    using IntersectionType = std::vector<Point>;

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

public:
    /*!
     *  \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
341
     * \param is If the geometries collide, is holds the corner points of
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
     *        the intersection object in global coordinates.
     */
    template<int dimIntersection>
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& is)
    {
        if (dimIntersection != 2)
            DUNE_THROW(Dune::NotImplemented, "Only simplex intersections are currently implemented!");

        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)
            return intersection<dimIntersection>(a, b, c, p, q, is);

        else if (geo1.corners() == 4)
        {
            const auto d = geo1.corner(3);
            if (intersection<dimIntersection>(a, b, d, p, q, is)) return true;
            else if (intersection<dimIntersection>(a, d, c, p, q, is)) return true;
            else return false;
        }

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

    // triangle--segment intersection with points as input
    template<int dimIntersection>
    static bool intersection(const Point& a, const Point& b, const Point& c,
                             const Point& p, const Point& q,
                             IntersectionType& is)
    {
        if (dimIntersection != 2)
            DUNE_THROW(Dune::NotImplemented, "Only simplex intersections are currently implemented!");

        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;
        const auto eps = eps_*ab.two_norm2()*qp.two_norm();
        using std::abs;
        if (abs(denom) < eps)
            return false;

        // 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);
        is = {ip};
423
        return true;
424
425
426
427
428
    }
};

} // end namespace Dumux

429
# endif