test_1d3d_intersection.cc 17.4 KB
Newer Older
1
2
3
4
5
6
7
#include <config.h>

#include <iostream>
#include <algorithm>

#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
8
#include <dune/geometry/type.hh>
9
10
#include <dune/geometry/multilineargeometry.hh>

11
#include <dumux/common/math.hh>
12

13
14
#include "transformation.hh"
#include <dumux/geometry/geometryintersection.hh>
15
16
17
18
19
20
21
22
23
24

#ifndef DOXYGEN
template<int dimworld = 3>
Dune::MultiLinearGeometry<double, 1, dimworld>
makeLine(std::initializer_list<Dune::FieldVector<double, dimworld>>&& c)
{
    return {Dune::GeometryTypes::line, c};
}

template<int dimworld = 3>
25
bool testIntersection(const Dune::MultiLinearGeometry<double, dimworld, dimworld>& polyhedron,
26
                      const Dune::MultiLinearGeometry<double, 1, dimworld>& line,
27
                      bool foundExpected, bool verbose)
28
29
30
{
    using Test = Dumux::GeometryIntersection<Dune::MultiLinearGeometry<double,dimworld,dimworld>,
                                             Dune::MultiLinearGeometry<double,1,dimworld> >;
31
    typename Test::Intersection intersection;
32
    bool found = Test::intersection(polyhedron, line, intersection);
33
    if (!found && foundExpected)
34
35
36
    {
        std::cerr << "  Failed detecting intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
    }
37
    else if (found && !foundExpected)
38
39
40
41
42
43
44
45
46
47
    {
        std::cerr << "  Found false positive: intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
    }
    if (verbose)
    {
        if (found && foundExpected)
            std::cout << "  Found intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
        else if (!found && !foundExpected)
            std::cout << "  No intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
    }
48
49
50
51
    return (found == foundExpected);
}
#endif

52
53
54
55
namespace Dumux {

template<class Transformation>
void runIntersectionTest(std::vector<bool>& returns, const Transformation& transform, bool verbose)
56
{
57
58
    using Points = std::vector<Dune::FieldVector<double, 3>>;
    using Geo = Dune::MultiLinearGeometry<double, 3, 3>;
59

60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
    // test tetrahedron-line intersections
    if (verbose) std::cout << "\n  -- Test tetrahedron-line intersections" << std::endl;

    auto cornersTetrahedron = Points({{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}});
    std::transform(cornersTetrahedron.begin(), cornersTetrahedron.end(), cornersTetrahedron.begin(),
                   [&](const auto& p) { return transform(p); });
    const auto tetrahedron = Geo(Dune::GeometryTypes::tetrahedron, cornersTetrahedron);

    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 0.0, 0.0})}, {transform({1.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({1.0, 0.0, 0.0})}, {transform({0.0, 1.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 1.0, 0.0})}, {transform({0.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 0.0, 0.0})}, {transform({0.0, 0.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({1.0, 0.0, 0.0})}, {transform({0.0, 0.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 1.0, 0.0})}, {transform({0.0, 0.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 0.0, 0.5})}, {transform({0.5, 0.0, 0.5})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 0.0, 0.5})}, {transform({0.0, 0.5, 0.5})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.5, 0.0, 0.5})}, {transform({0.0, 0.5, 0.5})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 0.0, 1.0})}, {transform({0.5, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 0.0, 1.0})}, {transform({0.0, 0.5, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 0.0, 1.0})}, {transform({0.5, 0.5, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 0.0, 0.0})}, {transform({1.0, 1.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.25, 0.25, 0.0})}, {transform({0.25, 0.25, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({-1.0, 0.25, 0.5})}, {transform({1.0, 0.25, 0.5})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({1.0, 1.0, 1.0})}, {transform({-1.0, -1.0, -1.0})}}), true, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({1.5, 0.0, 0.5})}, {transform({0.0, 1.5, 0.5})}}), false, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 0.0, 0.0})}, {transform({0.0, 0.0, -1.0})}}), false, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({1.0, 1.0, 0.0})}, {transform({0.0, 0.0, 2.0})}}), false, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({1.0, 0.0, 0.1})}, {transform({0.0, 1.0, 0.1})}}), false, verbose));
    returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 0.0, -0.1})}, {transform({1.0, 1.0, -0.1})}}), false, verbose));

    // test hexahedron-line intersections
    if (verbose) std::cout << "\n  -- Test hexahedron-line intersections" << std::endl;

    auto cornersHexahedron = Points({{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {1.0, 1.0, 0.0},
                                     {0.0, 0.0, 1.0}, {1.0, 0.0, 1.0}, {0.0, 1.0, 1.0}, {1.0, 1.0, 1.0}});
    std::transform(cornersHexahedron.begin(), cornersHexahedron.end(), cornersHexahedron.begin(),
                   [&](const auto& p) { return transform(p); });
    auto hexahedron = Geo(Dune::GeometryTypes::hexahedron, cornersHexahedron);

    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 0.0, 0.0})}, {transform({1.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 0.0, 0.0})}, {transform({0.0, 1.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 0.0, 0.0})}, {transform({0.0, 0.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.0, 0.0, 0.0})}, {transform({0.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.0, 0.0, 0.0})}, {transform({1.0, 0.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.0, 0.0, 0.0})}, {transform({1.0, 1.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 1.0, 0.0})}, {transform({0.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 1.0, 0.0})}, {transform({0.0, 1.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 1.0, 0.0})}, {transform({1.0, 1.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.0, 1.0, 0.0})}, {transform({1.0, 1.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.0, 1.0, 0.0})}, {transform({0.0, 1.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.0, 1.0, 0.0})}, {transform({1.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 0.0, 1.0})}, {transform({1.0, 0.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 0.0, 1.0})}, {transform({0.0, 1.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 0.0, 1.0})}, {transform({0.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.0, 0.0, 1.0})}, {transform({0.0, 0.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.0, 0.0, 1.0})}, {transform({1.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.0, 0.0, 1.0})}, {transform({1.0, 1.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 1.0, 1.0})}, {transform({0.0, 0.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 1.0, 1.0})}, {transform({0.0, 1.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 1.0, 1.0})}, {transform({1.0, 1.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.0, 1.0, 1.0})}, {transform({1.0, 1.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.0, 1.0, 1.0})}, {transform({0.0, 1.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.0, 1.0, 1.0})}, {transform({1.0, 0.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 0.0, 0.0})}, {transform({1.0, 1.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.5, 0.5, 0.5})}, {transform({0.5, 0.5, -2.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.5, 0.0, 0.5})}, {transform({0.5, 1.0, 0.5})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 0.5, 0.5})}, {transform({1.0, 0.5, 0.5})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.5, 0.5, 0.0})}, {transform({0.5, 0.5, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 0.0, 2.0})}, {transform({1.0, 1.0, 2.0})}}), false, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 0.0, 1.1})}, {transform({1.0, 1.0, 1.1})}}), false, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.1, 1.1, 0.0})}, {transform({1.1, 1.1, 1.0})}}), false, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.1, 0.0, 0.0})}, {transform({1.1, 1.0, 1.0})}}), false, verbose));
    returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, -0.1, 0.0})}, {transform({1.0, -0.1, 0.0})}}), false, verbose));

    // test pyramid-line intersections
    if (verbose) std::cout << "\n  -- Test pyramid-line intersections" << std::endl;

    auto cornersPyramid = Points({{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {1.0, 1.0, 0.0}, {0.5, 0.5, 1.0}});
    std::transform(cornersPyramid.begin(), cornersPyramid.end(), cornersPyramid.begin(),
                   [&](const auto& p) { return transform(p); });
    auto pyramid = Geo(Dune::GeometryTypes::pyramid, cornersPyramid);

    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.0, 0.0, 0.0})}, {transform({1.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({1.0, 0.0, 0.0})}, {transform({1.0, 1.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({1.0, 1.0, 0.0})}, {transform({0.0, 1.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.0, 1.0, 0.0})}, {transform({0.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.5, 0.5, 1.0})}, {transform({0.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.5, 0.5, 1.0})}, {transform({1.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.5, 0.5, 1.0})}, {transform({0.0, 1.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.5, 0.5, 1.0})}, {transform({1.0, 1.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.5, 0.5, 1.0})}, {transform({0.5, 0.5, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.25, 0.25, 0.5})}, {transform({0.75, 0.25, 0.5})}}), true, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.75, 0.25, 0.5})}, {transform({0.75, 0.75, 0.5})}}), true, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.75, 0.75, 0.5})}, {transform({0.25, 0.75, 0.5})}}), true, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.25, 0.75, 0.5})}, {transform({0.25, 0.25, 0.5})}}), true, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.0, 0.0, 1.0})}, {transform({1.0, 0.0, 1.0})}}), false, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.0, 0.0, 1.0})}, {transform({0.0, 1.0, 1.0})}}), false, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.0, 0.0, -0.1})}, {transform({1.0, 1.0, -0.1})}}), false, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.0, 1.1, 0.0})}, {transform({1.0, 1.1, 0.0})}}), false, verbose));
    returns.push_back(testIntersection(pyramid, makeLine({{transform({0.4, 0.0, 1.0})}, {transform({0.4, 1.0, 1.0})}}), false, verbose));

    // test prism-line intersections
    if (verbose) std::cout << "\n  -- Test prism-line intersections" << std::endl;

    auto cornersPrism = Points({{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0},
                                {0.0, 0.0, 1.0}, {1.0, 0.0, 1.0}, {0.0, 1.0, 1.0}});
    std::transform(cornersPrism.begin(), cornersPrism.end(), cornersPrism.begin(),
                   [&](const auto& p) { return transform(p); });
    auto prism = Geo(Dune::GeometryTypes::prism, cornersPrism);

    returns.push_back(testIntersection(prism, makeLine({{transform({0.0, 0.0, 0.0})}, {transform({1.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({1.0, 0.0, 0.0})}, {transform({0.0, 1.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({0.0, 1.0, 0.0})}, {transform({0.0, 0.0, 0.0})}}), true, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({0.0, 0.0, 1.0})}, {transform({1.0, 0.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({1.0, 0.0, 1.0})}, {transform({0.0, 1.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({0.0, 1.0, 1.0})}, {transform({0.0, 0.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({0.0, 0.0, 0.0})}, {transform({0.0, 0.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({1.0, 0.0, 0.0})}, {transform({1.0, 0.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({0.0, 1.0, 0.0})}, {transform({0.0, 1.0, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({0.25, 0.25, 0.0})}, {transform({0.25, 0.25, 1.0})}}), true, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({0.0, 0.0, 0.5})}, {transform({1.0, 0.0, 0.5})}}), true, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({1.0, 0.0, 0.5})}, {transform({0.0, 1.0, 0.5})}}), true, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({0.0, 1.0, 0.5})}, {transform({0.0, 0.0, 0.5})}}), true, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({1.0, 1.0, 0.0})}, {transform({1.0, 1.0, 1.0})}}), false, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({2.0, 0.0, 0.5})}, {transform({0.0, 2.0, 0.5})}}), false, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({1.1, 0.0, 0.0})}, {transform({1.1, 0.0, 1.0})}}), false, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({-0.1, 0.0, 1.0})}, {transform({-0.1, 1.0, 1.0})}}), false, verbose));
    returns.push_back(testIntersection(prism, makeLine({{transform({1.0, 0.0, 1.1})}, {transform({0.0, 1.0, 1.1})}}), false, verbose));
}
189

190
} // end namespace Dumux
191

192
193
194
int main (int argc, char *argv[])
{
    using namespace Dumux;
195
196
197

    // collect returns to determine exit code
    std::vector<bool> returns;
198
    constexpr bool verbose = false;
199

200
201
202
203
204
205
    using Vec = Dune::FieldVector<double, 3>;
    for (const double scaling : {1.0, 1e3, 1e12, 1e-12})
        for (const double translation : {0.0, 1.0})
            for (const double angle : {0.0, 0.2*M_PI, 0.5*M_PI, 0.567576567*M_PI, M_PI})
                for (const auto& rotAxis : {Vec(std::sqrt(3.0)/3.0), Vec({std::sqrt(2.0)/2.0, std::sqrt(2.0)/2.0, 0.0})})
                    runIntersectionTest(returns, make3DTransformation<double>(scaling, Vec(translation), rotAxis, angle, true), verbose);
206

207
    // determine the exit code
208
    if (std::any_of(returns.begin(), returns.end(), std::logical_not<bool>{}))
209
210
        return 1;

211
212
213
    std::cout << "\n++++++++++++++++++++++\n"
              << "All tests passed!"
              << "\n++++++++++++++++++++++" << std::endl;
214
215
216

    return 0;
}