GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/geometry/makegeometry.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 48 59 81.4%
Functions: 4 5 80.0%
Branches: 67 244 27.5%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup Geometry
10 * \brief Create Dune geometries from user-specified points
11 */
12 #ifndef DUMUX_GEOMETRY_MAKE_GEOMETRY_HH
13 #define DUMUX_GEOMETRY_MAKE_GEOMETRY_HH
14
15 #include <vector>
16 #include <array>
17 #include <limits>
18 #include <dune/common/fvector.hh>
19 #include <dune/common/fmatrix.hh>
20 #include <dune/common/exceptions.hh>
21 #include <dune/geometry/multilineargeometry.hh>
22 #include <dumux/common/math.hh>
23 #include <dumux/geometry/grahamconvexhull.hh>
24
25 namespace Dumux {
26
27 /*!
28 * \ingroup Geometry
29 * \brief Checks if four points lie within the same plane
30 */
31 template<class CoordScalar>
32 2568 bool pointsAreCoplanar(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points, const CoordScalar scale)
33 {
34
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2568 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2568 times.
 5136 if (points.size() != 4)
35 ✗ DUNE_THROW(Dune::InvalidStateException, "Check only works for 4 points!");
36
37 // (see "Real-Time Collision Detection" by Christer Ericson)
38 2568 Dune::FieldMatrix<CoordScalar, 4, 4> M;
39
2/2
✓ Branch 0 taken 7704 times.
✓ Branch 1 taken 2568 times.
 10272 for (int i = 0; i < 3; ++i )
40 61632 M[i] = {points[0][i], points[1][i], points[2][i], points[3][i]};
41 2568 M[3] = {1.0*scale, 1.0*scale, 1.0*scale, 1.0*scale};
42
43 using std::abs;
44 2568 return abs(M.determinant()) < 1.5e-7*scale*scale*scale*scale;
45 }
46
47 /*!
48 * \ingroup Geometry
49 * \brief Checks if four points lie within the same plane.
50 */
51 template<class CoordScalar>
52 18 bool pointsAreCoplanar(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points)
53 {
54 18 Dune::FieldVector<CoordScalar, 3> bBoxMin(std::numeric_limits<CoordScalar>::max());
55 18 Dune::FieldVector<CoordScalar, 3> bBoxMax(std::numeric_limits<CoordScalar>::lowest());
56
4/4
✓ Branch 0 taken 72 times.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 72 times.
✓ Branch 3 taken 18 times.
 126 for (const auto& p : points)
57 {
58
2/2
✓ Branch 0 taken 216 times.
✓ Branch 1 taken 72 times.
 288 for (int i=0; i<3; i++)
59 {
60 using std::min;
61 using std::max;
62
8/8
✓ Branch 0 taken 63 times.
✓ Branch 1 taken 153 times.
✓ Branch 2 taken 63 times.
✓ Branch 3 taken 153 times.
✓ Branch 4 taken 63 times.
✓ Branch 5 taken 153 times.
✓ Branch 6 taken 90 times.
✓ Branch 7 taken 126 times.
 711 bBoxMin[i] = min(bBoxMin[i], p[i]);
63
6/6
✓ Branch 0 taken 90 times.
✓ Branch 1 taken 126 times.
✓ Branch 2 taken 90 times.
✓ Branch 3 taken 126 times.
✓ Branch 4 taken 90 times.
✓ Branch 5 taken 126 times.
 828 bBoxMax[i] = max(bBoxMax[i], p[i]);
64 }
65 }
66
67 18 const auto size = (bBoxMax - bBoxMin).two_norm();
68
69 18 return pointsAreCoplanar(points, size);
70 }
71
72 /*!
73 * \ingroup Geometry
74 * \brief Returns a vector of points following the dune ordering.
75 * Convenience method that creates a temporary object in case no array of orientations is desired.
76 *
77 * \param points The user-specified vector of points (potentially in wrong order).
78 */
79 template<class CoordScalar>
80 std::vector<Dune::FieldVector<CoordScalar, 3>> getReorderedPoints(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points)
81 {
82 std::array<int, 4> tmp;
83 return getReorderedPoints(points, tmp);
84 }
85
86 /*!
87 * \ingroup Geometry
88 * \brief Returns a vector of points following the dune ordering.
89 *
90 * \param points The user-specified vector of points (potentially in wrong order).
91 * \param orientations An array of orientations that can be useful for further processing.
92 */
93 template<class CoordScalar>
94 2550 std::vector<Dune::FieldVector<CoordScalar, 3>> getReorderedPoints(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points,
95 std::array<int, 4>& orientations)
96 {
97
2/4
✓ Branch 0 taken 2550 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2550 times.
✗ Branch 3 not taken.
 5100 if(points.size() == 4)
98 {
99 2550 auto& p0 = points[0];
100 2550 auto& p1 = points[1];
101 2550 auto& p2 = points[2];
102 2550 auto& p3 = points[3];
103
104 // check if the points define a proper quadrilateral
105 7650 const auto normal = crossProduct((p1 - p0), (p2 - p0));
106
107 2550 orientations = { getOrientation(p0, p3, p2, normal),
108 2550 getOrientation(p0, p3, p1, normal),
109 2550 getOrientation(p2, p1, p0, normal),
110 2550 getOrientation(p2, p1, p3, normal) };
111
112
113 // check if the points follow the dune ordering (see http://www.dcs.gla.ac.uk/~pat/52233/slides/Geometry1x1.pdf)
114
3/12
✗ Branch 0 not taken.
✓ Branch 1 taken 2550 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2550 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2550 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 7650 const bool diagonalsIntersect = (orientations[0] != orientations[1]) && (orientations[2] != orientations[3]);
115
116 // the points conform with the dune ordering
117 if(diagonalsIntersect)
118 ✗ return points;
119
120 // the points do not conform with the dune ordering, re-order
121 using GlobalPosition = Dune::FieldVector<CoordScalar, 3>;
122
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2550 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2550 times.
 5100 if(!diagonalsIntersect && orientations[0] == 1)
123 5100 return std::vector<GlobalPosition>{p1, p0, p2, p3};
124 ✗ else if(!diagonalsIntersect && orientations[0] == -1)
125 ✗ return std::vector<GlobalPosition>{p3, p1, p0, p2};
126 else
127 ✗ DUNE_THROW(Dune::InvalidStateException, "Could not reorder points");
128 }
129 else
130 ✗ DUNE_THROW(Dune::NotImplemented, "Reorder for " << points.size() << " points.");
131 }
132
133 /*!
134 * \ingroup Geometry
135 * \brief Creates a dune quadrilateral geometry given 4 corner points.
136 *
137 * \tparam CoordScalar The CoordScalar type.
138 * \tparam enableSanityCheck Turn on/off sanity check and reordering of points
139 * \param points The user-specified vector of points (potentially in wrong order).
140 */
141 template<class CoordScalar, bool enableSanityCheck = true>
142 2550 auto makeDuneQuadrilaterial(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points)
143 {
144
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2550 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2550 times.
 5100 if (points.size() != 4)
145 ✗ DUNE_THROW(Dune::InvalidStateException, "A quadrilateral needs 4 corner points!");
146
147 using GlobalPosition = Dune::FieldVector<CoordScalar, 3>;
148 static constexpr auto coordDim = GlobalPosition::dimension;
149 static constexpr auto dim = coordDim-1;
150 using GeometryType = Dune::MultiLinearGeometry<CoordScalar, dim, coordDim>;
151
152 // if no sanity check if desired, use the given points directly to construct the geometry
153 if (!enableSanityCheck)
154 return GeometryType(Dune::GeometryTypes::quadrilateral, points);
155
156 // compute size
157 2550 Dune::FieldVector<CoordScalar, 3> bBoxMin(std::numeric_limits<CoordScalar>::max());
158 2550 Dune::FieldVector<CoordScalar, 3> bBoxMax(std::numeric_limits<CoordScalar>::lowest());
159
4/4
✓ Branch 0 taken 10200 times.
✓ Branch 1 taken 2550 times.
✓ Branch 2 taken 10200 times.
✓ Branch 3 taken 2550 times.
 15300 for (const auto& p : points)
160 {
161
2/2
✓ Branch 0 taken 30600 times.
✓ Branch 1 taken 10200 times.
 40800 for (int i = 0; i < 3; i++)
162 {
163 using std::min;
164 using std::max;
165
8/8
✓ Branch 0 taken 12743 times.
✓ Branch 1 taken 17857 times.
✓ Branch 2 taken 12743 times.
✓ Branch 3 taken 17857 times.
✓ Branch 4 taken 12743 times.
✓ Branch 5 taken 17857 times.
✓ Branch 6 taken 12998 times.
✓ Branch 7 taken 17602 times.
 104543 bBoxMin[i] = min(bBoxMin[i], p[i]);
166
6/6
✓ Branch 0 taken 12998 times.
✓ Branch 1 taken 17602 times.
✓ Branch 2 taken 12998 times.
✓ Branch 3 taken 17602 times.
✓ Branch 4 taken 12998 times.
✓ Branch 5 taken 17602 times.
 117796 bBoxMax[i] = max(bBoxMax[i], p[i]);
167 }
168 }
169
170 2550 const auto size = (bBoxMax - bBoxMin).two_norm();
171
172 // otherwise, perform a number of checks and corrections
173
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2550 times.
 2550 if (!pointsAreCoplanar(points, size))
174 ✗ DUNE_THROW(Dune::InvalidStateException, "Points do not lie within a plane");
175
176
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 5100 auto corners = grahamConvexHull<2>(points);
177
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2550 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2550 times.
 5100 if (corners.size() != 4)
178 ✗ DUNE_THROW(Dune::InvalidStateException, "Points do not span a strictly convex polygon!");
179
180 // make sure points conform with dune ordering
181 std::array<int, 4> orientations;
182
2/4
✓ Branch 1 taken 2550 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2550 times.
 2550 corners = getReorderedPoints(corners, orientations);
183
184
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 2550 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2550 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2550 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2550 times.
 10200 if (std::any_of(orientations.begin(), orientations.end(), [](auto i){ return i == 0; }))
185 ✗ DUNE_THROW(Dune::InvalidStateException, "More than two points lie on the same line.");
186
187
2/6
✓ Branch 1 taken 2550 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2550 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 5100 const auto quadrilateral = GeometryType(Dune::GeometryTypes::quadrilateral, corners);
188
189 2550 const auto eps = 1e-7;
190
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2550 times.
 2550 if (quadrilateral.volume() < eps*size*size)
191 ✗ DUNE_THROW(Dune::InvalidStateException, "Something went wrong, geometry has area of zero");
192
193
2/4
✓ Branch 1 taken 2550 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2550 times.
✗ Branch 4 not taken.
 2550 return quadrilateral;
194 }
195
196
197
198 } // end namespace Dumux
199
200 #endif
201