geometryintersection.hh 7.31 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
29
30
31
32
33
34
#include <dune/geometry/referenceelements.hh>

#include <dumux/common/math.hh>

namespace Dumux
{

35
/*!
36
 * \ingroup Common
37
38
39
40
 * \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.
 */
41
42
43
44
45
template
<class Geometry1, class Geometry2,
  int dimworld = Geometry1::coorddimension,
  int dim1 = Geometry1::mydimension,
  int dim2 = Geometry2::mydimension>
46
class GeometryIntersection
47
48
{
public:
49
50
51
52
    using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
    using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
    using IntersectionType = std::vector<GlobalPosition>;

53
    //! Determine if the two geometries intersect and compute the intersection corners
54
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& intersection)
55
    {
56
57
        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 = "
58
59
60
61
62
63
                                          << dimworld << ", dim1 = " << dim1 << ", dim2 = " << dim2);
    }
};

//! Geometry collision detection with 3d and 1d geometry in 3d space
template <class Geometry1, class Geometry2>
64
class GeometryIntersection<Geometry1, Geometry2, 3, 3, 1>
65
{
66
67
68
69
70
71
72
73
74
    enum { dimworld = 3 };
    enum { dim1 = 3 };
    enum { dim2 = 1 };
    enum { dimis = dim2 }; // the intersection dimension

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

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

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.
86
87
     * \param geo1/geo2 The geometries to intersect
     * \param intersection If the geometries collide intersection holds the corner points of
88
     *        the intersection object in global coordinates.
89
     */
90
    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& intersection)
91
    {
92
        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
93
94
95
96
97
98
99
100

        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
101
102
        ctype tfirst = 0.0;
        ctype tlast = 1.0;
103

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

            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
171
172
        // 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;
173
174
175
176
177
    }
};

} // end namespace Dumux

178
# endif