geometrycollision.hh 6.95 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/*****************************************************************************
 *   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
 * \brief A class for collision detection of two geometries
20
 *        and computation of intersection corners
21
 */
22
23
#ifndef DUMUX_GEOMETRY_COLLISION_HH
#define DUMUX_GEOMETRY_COLLISION_HH
24
25

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

#include <dumux/common/math.hh>

namespace Dumux
{

34
35
36
37
38
/*!
 * \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.
 */
39
40
41
42
43
44
45
template
<class Geometry1, class Geometry2,
  int dimworld = Geometry1::coorddimension,
  int dim1 = Geometry1::mydimension,
  int dim2 = Geometry2::mydimension>
class GeometryCollision
{
46
47
48
    using Scalar = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
    using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;

49
public:
50
51
    //! Determine if the two geometries intersect and compute the intersection corners
    static bool collide(const Geometry1& geo1, const Geometry2& geo2, std::vector<GlobalPosition>& intersection)
52
53
    {
        static_assert(dimworld == Geometry2::coorddimension, "Can only collide geometries of same coordinate dimension");
54
        DUNE_THROW(Dune::NotImplemented, "Geometry collision detection with intersection computation not implemented for dimworld = "
55
56
57
58
59
60
61
62
63
64
65
                                          << dimworld << ", dim1 = " << dim1 << ", dim2 = " << dim2);
    }
};

//! Geometry collision detection with 3d and 1d geometry in 3d space
template <class Geometry1, class Geometry2>
class GeometryCollision<Geometry1, Geometry2, 3, 3, 1>
{
    static const int dimworld = 3;
    static const int dim1 = 3;
    static const int dim2 = 1;
66
67
    static const int dimis = dim2; // the intersection dimension
    using Scalar = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
68
    using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim1>;
69
    using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
70
71
72
73
74
75
76
77
78

    static constexpr Scalar eps_ = 1.5e-7; // base epsilon for floating point comparisons

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

        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
        Scalar tfirst = 0;
        Scalar tlast = 1;

97
98
        std::vector<std::vector<int>> facets;
        // sort facet corner so that normal n = (p1-p0)x(p2-p0) always points outwards
99
100
        switch (geo1.corners())
        {
101
            // todo compute intersection geometries
102
            case 8: // hexahedron
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
                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:
                DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type " << geo1.type() << ", "<< geo1.corners() << " corners.");
        }

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

120
            auto n = crossProduct(v0, v1);
121
122
123
124
125
126
127
128
129
130
131
132
133
            n /= n.two_norm();

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

            // if denominator is zero the segment in parallel to
            // the plane. If the distance is positive there is no intersection
            if (std::abs(denom) < eps)
            {
                if (dist > eps)
                    return false;
            }
            else // not parallel: compute line-plane intersection
134
            {
135
136
137
                const Scalar t = -dist / denom;
                // if entering half space cut tfirst if t is larger
                if (std::signbit(denom))
138
                {
139
140
                    if (t > tfirst)
                        tfirst = t;
141
                }
142
143
144
145
146
147
148
149
150
151
152
                // 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;

153
154
            }
        }
155
156
157
158
        // 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]
        intersection = {geo2.global(tfirst), geo2.global(tlast)};
        return true;
159
160
161
162
163
    }
};

} // end namespace Dumux

164
# endif