| 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-FileCopyrightText: 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 Functionality to triangulate point clouds | ||
| 11 | * \note Most of the implemented algorithms currently do not scale for large point clouds | ||
| 12 | * and are only meant to be used with a small number of points | ||
| 13 | */ | ||
| 14 | #ifndef DUMUX_GEOMETRY_TRIANGULATION_HH | ||
| 15 | #define DUMUX_GEOMETRY_TRIANGULATION_HH | ||
| 16 | |||
| 17 | #include <vector> | ||
| 18 | #include <array> | ||
| 19 | #include <algorithm> | ||
| 20 | #include <type_traits> | ||
| 21 | #include <tuple> | ||
| 22 | |||
| 23 | #include <dune/common/exceptions.hh> | ||
| 24 | #include <dune/common/fvector.hh> | ||
| 25 | |||
| 26 | #include <dumux/common/math.hh> | ||
| 27 | #include <dumux/geometry/intersectspointsimplex.hh> | ||
| 28 | #include <dumux/geometry/grahamconvexhull.hh> | ||
| 29 | |||
| 30 | namespace Dumux { | ||
| 31 | namespace TriangulationPolicy { | ||
| 32 | |||
| 33 | //! Policy that expects a point cloud that represents a convex | ||
| 34 | //! hull. Inserts the mid point and connects it to the points of the hull. | ||
| 35 | struct MidPointPolicy {}; | ||
| 36 | |||
| 37 | //! Policy that first finds the convex hull | ||
| 38 | //! and then uses the mid point policy to construct a triangulation | ||
| 39 | struct ConvexHullPolicy {}; | ||
| 40 | |||
| 41 | //! Delaunay-type triangulations. | ||
| 42 | struct DelaunayPolicy {}; | ||
| 43 | |||
| 44 | #ifndef DOXYGEN | ||
| 45 | namespace Detail { | ||
| 46 | using DefaultDimPolicies = std::tuple<DelaunayPolicy, MidPointPolicy, ConvexHullPolicy>; | ||
| 47 | } // end namespace Detail | ||
| 48 | #endif | ||
| 49 | |||
| 50 | //! Default policy for a given dimension | ||
| 51 | template<int dim, int dimWorld> | ||
| 52 | using DefaultPolicy = std::tuple_element_t<dim-1, Detail::DefaultDimPolicies>; | ||
| 53 | |||
| 54 | } // end namespace TriangulationPolicy | ||
| 55 | |||
| 56 | /*! | ||
| 57 | * \ingroup Geometry | ||
| 58 | * \brief The default data type to store triangulations | ||
| 59 | * \note Stores each simplex separate and without connectivity information | ||
| 60 | * \note We usually use this for sub-triangulation to allow for quadrature rules and interpolation on intersection geometries | ||
| 61 | * This is neither meant to be used to store large amounts of data nor as a mesh-like object | ||
| 62 | */ | ||
| 63 | template<int dim, int dimWorld, class ctype> | ||
| 64 | using Triangulation = std::vector< std::array<Dune::FieldVector<ctype, dimWorld>, dim+1> >; | ||
| 65 | |||
| 66 | /*! | ||
| 67 | * \ingroup Geometry | ||
| 68 | * \brief Triangulate area given points of a convex hull (1d) | ||
| 69 | * \tparam dim Specifies the dimension of the resulting triangulation | ||
| 70 | * \tparam dimWorld The dimension of the coordinate space | ||
| 71 | * \tparam Policy Specifies the algorithm to be used for triangulation | ||
| 72 | * | ||
| 73 | * \note this specialization is for 1d discretization using segments | ||
| 74 | * \note sorts the points and connects them via segments | ||
| 75 | */ | ||
| 76 | template< int dim, int dimWorld, class Policy = TriangulationPolicy::DefaultPolicy<dim, dimWorld>, | ||
| 77 | class RandomAccessContainer, | ||
| 78 | std::enable_if_t< std::is_same_v<Policy, TriangulationPolicy::DelaunayPolicy> | ||
| 79 | && dim == 1, int> = 0 > | ||
| 80 | inline Triangulation<dim, dimWorld, typename RandomAccessContainer::value_type::value_type> | ||
| 81 | triangulate(const RandomAccessContainer& points) | ||
| 82 | { | ||
| 83 | using ctype = typename RandomAccessContainer::value_type::value_type; | ||
| 84 | using Point = Dune::FieldVector<ctype, dimWorld>; | ||
| 85 | |||
| 86 | static_assert(std::is_same_v<typename RandomAccessContainer::value_type, Point>, | ||
| 87 | "Triangulation expects Dune::FieldVector as point type"); | ||
| 88 | |||
| 89 | if (points.size() == 2) | ||
| 90 | return Triangulation<dim, dimWorld, ctype>({ {points[0], points[1]} }); | ||
| 91 | |||
| 92 | //! \todo sort points and create polyline | ||
| 93 | assert(points.size() > 1); | ||
| 94 | DUNE_THROW(Dune::NotImplemented, "1d triangulation for point cloud size > 2"); | ||
| 95 | } | ||
| 96 | |||
| 97 | /*! | ||
| 98 | * \ingroup Geometry | ||
| 99 | * \brief Triangulate area given points of a convex hull (2d) | ||
| 100 | * \tparam dim Specifies the dimension of the resulting triangulation | ||
| 101 | * \tparam dimWorld The dimension of the coordinate space | ||
| 102 | * \tparam Policy Specifies the algorithm to be used for triangulation | ||
| 103 | * | ||
| 104 | * \note this specialization is for 2d triangulations using mid point policy | ||
| 105 | * \note Assumes all points of the convex hull are coplanar | ||
| 106 | * \note This inserts a mid point and connects all corners with that point to triangles | ||
| 107 | * \note Assumes points are given as a ordered sequence representing the polyline forming the convex hull | ||
| 108 | */ | ||
| 109 | template< int dim, int dimWorld, class Policy = TriangulationPolicy::DefaultPolicy<dim, dimWorld>, | ||
| 110 | class RandomAccessContainer, | ||
| 111 | std::enable_if_t< std::is_same_v<Policy, TriangulationPolicy::MidPointPolicy> | ||
| 112 | && dim == 2, int> = 0 > | ||
| 113 | inline Triangulation<dim, dimWorld, typename RandomAccessContainer::value_type::value_type> | ||
| 114 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 110812 times.
|
110812 | triangulate(const RandomAccessContainer& convexHullPoints) |
| 115 | { | ||
| 116 | using ctype = typename RandomAccessContainer::value_type::value_type; | ||
| 117 | using Point = Dune::FieldVector<ctype, dimWorld>; | ||
| 118 | using Triangle = std::array<Point, 3>; | ||
| 119 | |||
| 120 | static_assert(std::is_same_v<typename RandomAccessContainer::value_type, Point>, | ||
| 121 | "Triangulation expects Dune::FieldVector as point type"); | ||
| 122 | |||
| 123 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 110812 times.
|
110812 | if (convexHullPoints.size() < 3) |
| 124 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Try to triangulate point cloud with less than 3 points!"); | |
| 125 | |||
| 126 |
2/2✓ Branch 0 taken 3911 times.
✓ Branch 1 taken 106901 times.
|
110812 | if (convexHullPoints.size() == 3) |
| 127 | 3911 | return std::vector<Triangle>(1, {convexHullPoints[0], convexHullPoints[1], convexHullPoints[2]}); | |
| 128 | |||
| 129 | 106901 | Point midPoint(0.0); | |
| 130 |
2/2✓ Branch 0 taken 430380 times.
✓ Branch 1 taken 106901 times.
|
537281 | for (const auto& p : convexHullPoints) |
| 131 | 430380 | midPoint += p; | |
| 132 | 106901 | midPoint /= convexHullPoints.size(); | |
| 133 | |||
| 134 | 106901 | std::vector<Triangle> triangulation; | |
| 135 |
1/2✓ Branch 1 taken 106901 times.
✗ Branch 2 not taken.
|
106901 | triangulation.reserve(convexHullPoints.size()); |
| 136 | |||
| 137 |
2/2✓ Branch 0 taken 323479 times.
✓ Branch 1 taken 106901 times.
|
430380 | for (std::size_t i = 0; i < convexHullPoints.size()-1; ++i) |
| 138 |
1/2✓ Branch 1 taken 323479 times.
✗ Branch 2 not taken.
|
323479 | triangulation.emplace_back(Triangle{midPoint, convexHullPoints[i], convexHullPoints[i+1]}); |
| 139 | |||
| 140 |
1/4✓ Branch 1 taken 106901 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
106901 | triangulation.emplace_back(Triangle{midPoint, convexHullPoints[convexHullPoints.size()-1], convexHullPoints[0]}); |
| 141 | |||
| 142 | 106901 | return triangulation; | |
| 143 | 106901 | } | |
| 144 | |||
| 145 | /*! | ||
| 146 | * \ingroup Geometry | ||
| 147 | * \brief Triangulate area given points (2d) | ||
| 148 | * \tparam dim Specifies the dimension of the resulting triangulation | ||
| 149 | * \tparam dimWorld The dimension of the coordinate space | ||
| 150 | * \tparam Policy Specifies the algorithm to be used for triangulation | ||
| 151 | * | ||
| 152 | * \note this specialization is for 2d triangulations using the convex hull policy | ||
| 153 | * That means we first construct the convex hull of the points and then | ||
| 154 | * triangulate the convex hull using the midpoint policy | ||
| 155 | * \note Assumes points are unique and not all colinear (will throw a Dune::InvalidStateException) | ||
| 156 | */ | ||
| 157 | template< int dim, int dimWorld, class Policy = TriangulationPolicy::DefaultPolicy<dim, dimWorld>, | ||
| 158 | class RandomAccessContainer, | ||
| 159 | std::enable_if_t< std::is_same_v<Policy, TriangulationPolicy::ConvexHullPolicy> | ||
| 160 | && dim == 2, int> = 0 > | ||
| 161 | inline Triangulation<dim, dimWorld, typename RandomAccessContainer::value_type::value_type> | ||
| 162 | 18383 | triangulate(const RandomAccessContainer& points) | |
| 163 | { | ||
| 164 | 18383 | const auto convexHullPoints = grahamConvexHull<2>(points); | |
| 165 |
1/2✓ Branch 1 taken 18383 times.
✗ Branch 2 not taken.
|
18383 | return triangulate<dim, dimWorld, TriangulationPolicy::MidPointPolicy>(convexHullPoints); |
| 166 | 18383 | } | |
| 167 | |||
| 168 | /*! | ||
| 169 | * \ingroup Geometry | ||
| 170 | * \brief Triangulate volume given a point cloud (3d) | ||
| 171 | * \tparam dim Specifies the dimension of the resulting triangulation | ||
| 172 | * \tparam dimWorld The dimension of the coordinate space | ||
| 173 | * \tparam Policy Specifies the algorithm to be used for triangulation | ||
| 174 | * | ||
| 175 | * \note this specialization is for 3d triangulations using the convex hull policy | ||
| 176 | * \note Assumes points are unique and not all coplanar | ||
| 177 | */ | ||
| 178 | template< int dim, int dimWorld, class Policy = TriangulationPolicy::DefaultPolicy<dim, dimWorld>, | ||
| 179 | class RandomAccessContainer, | ||
| 180 | std::enable_if_t< std::is_same_v<Policy, TriangulationPolicy::ConvexHullPolicy> | ||
| 181 | && dim == 3, int> = 0 > | ||
| 182 | inline Triangulation<dim, dimWorld, typename RandomAccessContainer::value_type::value_type> | ||
| 183 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3215 times.
|
3215 | triangulate(const RandomAccessContainer& points) |
| 184 | { | ||
| 185 | using ctype = typename RandomAccessContainer::value_type::value_type; | ||
| 186 | using Point = Dune::FieldVector<ctype, dimWorld>; | ||
| 187 | using Tetrahedron = std::array<Point, 4>; | ||
| 188 | |||
| 189 | static_assert(std::is_same_v<typename RandomAccessContainer::value_type, Point>, | ||
| 190 | "Triangulation expects Dune::FieldVector as point type"); | ||
| 191 | |||
| 192 | 3215 | const auto numPoints = points.size(); | |
| 193 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3215 times.
|
3215 | if (numPoints < 4) |
| 194 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Trying to create 3D triangulation of point cloud with less than 4 points!"); | |
| 195 | |||
| 196 |
2/2✓ Branch 0 taken 94 times.
✓ Branch 1 taken 3121 times.
|
3215 | if (numPoints == 4) |
| 197 | 94 | return std::vector<Tetrahedron>(1, {points[0], points[1], points[2], points[3]}); | |
| 198 | |||
| 199 | // compute the mid point of the point cloud (not the midpoint of the convex hull but this | ||
| 200 | // should not matter too much for the applications we have in mind here) | ||
| 201 | 3121 | Point midPoint(0.0); | |
| 202 | 3121 | Point lowerLeft(1e100); | |
| 203 | 3121 | Point upperRight(-1e100); | |
| 204 |
2/2✓ Branch 0 taken 25020 times.
✓ Branch 1 taken 3121 times.
|
28141 | for (const auto& p : points) |
| 205 | { | ||
| 206 | 100080 | midPoint += p; | |
| 207 |
2/2✓ Branch 0 taken 75060 times.
✓ Branch 1 taken 25020 times.
|
100080 | for (int i = 0; i < dimWorld; ++i) |
| 208 | { | ||
| 209 | using std::max; using std::min; | ||
| 210 |
4/4✓ Branch 0 taken 38353 times.
✓ Branch 1 taken 36707 times.
✓ Branch 2 taken 16664 times.
✓ Branch 3 taken 58396 times.
|
113413 | lowerLeft[i] = min(p[i], lowerLeft[i]); |
| 211 |
2/2✓ Branch 0 taken 16664 times.
✓ Branch 1 taken 58396 times.
|
91724 | upperRight[i] = max(p[i], upperRight[i]); |
| 212 | } | ||
| 213 | } | ||
| 214 | 3121 | midPoint /= numPoints; | |
| 215 | |||
| 216 | 3121 | auto magnitude = 0.0; | |
| 217 | using std::max; | ||
| 218 |
2/2✓ Branch 0 taken 9363 times.
✓ Branch 1 taken 3121 times.
|
12484 | for (int i = 0; i < dimWorld; ++i) |
| 219 |
2/2✓ Branch 0 taken 1519 times.
✓ Branch 1 taken 7844 times.
|
10882 | magnitude = max(upperRight[i] - lowerLeft[i], magnitude); |
| 220 | 3121 | const auto eps = 1e-7*magnitude; | |
| 221 | 3121 | const auto eps2 = eps*magnitude; | |
| 222 | |||
| 223 | // reserve memory conservatively to avoid reallocation | ||
| 224 | 3121 | std::vector<Tetrahedron> triangulation; | |
| 225 |
1/2✓ Branch 1 taken 3121 times.
✗ Branch 2 not taken.
|
3121 | triangulation.reserve(numPoints); |
| 226 | |||
| 227 | // make a buffer for storing coplanar points and indices | ||
| 228 | 3121 | std::vector<Point> coplanarPointBuffer; | |
| 229 |
3/6✓ Branch 0 taken 3105 times.
✓ Branch 1 taken 16 times.
✓ Branch 3 taken 3121 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
6226 | coplanarPointBuffer.reserve(std::min<std::size_t>(12, numPoints-1)); |
| 230 | |||
| 231 | // remember coplanar cluster planes | ||
| 232 | // coplanar clusters are uniquely identified by a point and a normal (plane) | ||
| 233 | // we only want to add each cluster once (when handling the first triangle in the cluster) | ||
| 234 | 3121 | std::vector<std::pair<int, Point>> coplanarClusters; | |
| 235 |
1/2✓ Branch 1 taken 3121 times.
✗ Branch 2 not taken.
|
3121 | coplanarClusters.reserve(numPoints/3); |
| 236 | |||
| 237 | // brute force algorithm: Try all possible triangles and check | ||
| 238 | // if they are triangles of the convex hull. This is achieved by checking if all | ||
| 239 | // other points are in the same half-space (side) | ||
| 240 | // the algorithm is O(n^4), so only do this for very small point clouds | ||
| 241 |
2/2✓ Branch 0 taken 25020 times.
✓ Branch 1 taken 3121 times.
|
28141 | for (int i = 0; i < numPoints; ++i) |
| 242 | { | ||
| 243 |
2/2✓ Branch 0 taken 88792 times.
✓ Branch 1 taken 25020 times.
|
113812 | for (int j = i+1; j < numPoints; ++j) |
| 244 | { | ||
| 245 |
2/2✓ Branch 0 taken 185889 times.
✓ Branch 1 taken 88792 times.
|
274681 | for (int k = j+1; k < numPoints; ++k) |
| 246 | { | ||
| 247 | 185889 | const auto pointI = points[i]; | |
| 248 | 185889 | const auto ab = points[j] - pointI; | |
| 249 | 185889 | const auto ac = points[k] - pointI; | |
| 250 |
2/2✓ Branch 1 taken 49095 times.
✓ Branch 2 taken 136794 times.
|
185889 | const auto normal = crossProduct(ab, ac); |
| 251 | |||
| 252 | // clear list of coplanar points w.r.t to triangle ijk | ||
| 253 | 185889 | coplanarPointBuffer.clear(); | |
| 254 | |||
| 255 | // check if this triangle ijk is admissible which means | ||
| 256 | // it is on the convex hull (all other points in the cloud are in the same half-space/side) | ||
| 257 | 557667 | const bool isAdmissible = [&]() | |
| 258 | { | ||
| 259 | // check if points are colinear and we can't form a triangle | ||
| 260 | // if so, skip this triangle | ||
| 261 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 185889 times.
|
371778 | if (normal.two_norm2() < eps2*eps2) |
| 262 | return false; | ||
| 263 | |||
| 264 | int marker = 0; // 0 means undecided side (or coplanar) | ||
| 265 |
2/2✓ Branch 0 taken 914190 times.
✓ Branch 1 taken 18895 times.
|
933085 | for (int m = 0; m < numPoints; ++m) |
| 266 | { | ||
| 267 |
6/6✓ Branch 0 taken 752910 times.
✓ Branch 1 taken 161280 times.
✓ Branch 2 taken 654013 times.
✓ Branch 3 taken 98897 times.
✓ Branch 4 taken 616197 times.
✓ Branch 5 taken 37816 times.
|
914190 | if (m != i && m != j && m != k) |
| 268 | { | ||
| 269 | // check scalar product with surface normal to decide side | ||
| 270 | 616197 | const auto ad = points[m] - pointI; | |
| 271 |
2/2✓ Branch 0 taken 496852 times.
✓ Branch 1 taken 119345 times.
|
1232394 | const auto sp = normal*ad; |
| 272 | |||
| 273 | // if the sign changes wrt the previous sign, the triangle is not part of the convex hull | ||
| 274 | using std::abs; using std::signbit; | ||
| 275 | 616197 | const bool coplanar = abs(sp) < eps2*magnitude; | |
| 276 |
4/4✓ Branch 0 taken 496852 times.
✓ Branch 1 taken 119345 times.
✓ Branch 2 taken 245017 times.
✓ Branch 3 taken 251835 times.
|
616197 | int newMarker = coplanar ? 0 : signbit(sp) ? -1 : 1; |
| 277 | |||
| 278 | // make decision for a side as soon as the next marker is != 0 | ||
| 279 | // keep track of the previous marker | ||
| 280 |
2/2✓ Branch 0 taken 163784 times.
✓ Branch 1 taken 452413 times.
|
616197 | if (marker == 0 && newMarker != 0) |
| 281 | 163784 | marker = newMarker; | |
| 282 | |||
| 283 | // if marker flips, not all points are on one side | ||
| 284 | // zero marker (undecided side / coplanar point) shouldn't abort the process | ||
| 285 |
2/2✓ Branch 0 taken 518890 times.
✓ Branch 1 taken 97307 times.
|
616197 | if (newMarker != 0 && marker != newMarker) |
| 286 | 166994 | return false; | |
| 287 | |||
| 288 | // handle possible coplanar points | ||
| 289 |
2/2✓ Branch 0 taken 119345 times.
✓ Branch 1 taken 399545 times.
|
518890 | if (coplanar) |
| 290 | { | ||
| 291 | using std::abs; | ||
| 292 |
4/4✓ Branch 0 taken 97597 times.
✓ Branch 1 taken 21748 times.
✓ Branch 2 taken 69687 times.
✓ Branch 3 taken 27910 times.
|
119345 | if (m < k && std::find_if( |
| 293 | coplanarClusters.begin(), coplanarClusters.end(), | ||
| 294 |
6/6✓ Branch 0 taken 198921 times.
✓ Branch 1 taken 158284 times.
✓ Branch 2 taken 137024 times.
✓ Branch 3 taken 220181 times.
✓ Branch 4 taken 67337 times.
✓ Branch 5 taken 69687 times.
|
693150 | [=](const auto& c){ return c.first == std::min(m, i) && abs(ab*c.second) < eps2*magnitude; } |
| 295 |
2/2✓ Branch 0 taken 69687 times.
✓ Branch 1 taken 27910 times.
|
97597 | ) != coplanarClusters.end()) |
| 296 | { | ||
| 297 | // this cluster has already been handled | ||
| 298 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 69675 times.
|
69687 | coplanarPointBuffer.clear(); |
| 299 | 69687 | return false; | |
| 300 | } | ||
| 301 | else | ||
| 302 | 49658 | coplanarPointBuffer.push_back(points[m]); | |
| 303 | } | ||
| 304 | } | ||
| 305 | } | ||
| 306 | |||
| 307 | // if there are coplanar points complete the cluster with i, j, and k | ||
| 308 | // and store the cluster information for future lookup | ||
| 309 |
2/2✓ Branch 0 taken 18383 times.
✓ Branch 1 taken 512 times.
|
18895 | if (!coplanarPointBuffer.empty()) |
| 310 | { | ||
| 311 | 18383 | coplanarPointBuffer.insert(coplanarPointBuffer.end(), { points[i], points[j], points[k] }); | |
| 312 | 18383 | coplanarClusters.emplace_back(std::make_pair(i, normal)); | |
| 313 | } | ||
| 314 | |||
| 315 | // we require that not all points are coplanar, so | ||
| 316 | // there will be always at least one non-coplanar point | ||
| 317 | // to check on which side the other points are. | ||
| 318 | // Hence, once we get here, the triangle or coplanar point cluster is part of the convex hull | ||
| 319 | return true; | ||
| 320 |
1/2✓ Branch 1 taken 185889 times.
✗ Branch 2 not taken.
|
185889 | }(); |
| 321 | |||
| 322 |
2/2✓ Branch 0 taken 18895 times.
✓ Branch 1 taken 166994 times.
|
185889 | if (isAdmissible) |
| 323 | { | ||
| 324 | // check if we have a cluster of coplanar points forming on of the | ||
| 325 | // faces of the convex hull, if yes, compute (2d) convex hull first and triangulate | ||
| 326 |
2/2✓ Branch 0 taken 18383 times.
✓ Branch 1 taken 512 times.
|
18895 | if (!coplanarPointBuffer.empty()) |
| 327 | { | ||
| 328 |
1/2✓ Branch 1 taken 18383 times.
✗ Branch 2 not taken.
|
18383 | const auto triangles = triangulate<2, 3, TriangulationPolicy::ConvexHullPolicy>(coplanarPointBuffer); |
| 329 |
2/2✓ Branch 0 taken 73804 times.
✓ Branch 1 taken 18383 times.
|
92187 | for (const auto& triangle : triangles) |
| 330 | { | ||
| 331 | 147608 | const auto ab = triangle[1] - triangle[0]; | |
| 332 | 73804 | const auto ac = triangle[2] - triangle[0]; | |
| 333 | 73804 | const auto normal = crossProduct(ab, ac); | |
| 334 | 147608 | const auto am = midPoint - triangle[0]; | |
| 335 |
2/2✓ Branch 0 taken 37057 times.
✓ Branch 1 taken 36747 times.
|
73804 | const auto sp = normal*am; |
| 336 | using std::signbit; | ||
| 337 | 73804 | const bool isBelow = signbit(sp); | |
| 338 |
2/2✓ Branch 0 taken 37057 times.
✓ Branch 1 taken 36747 times.
|
73804 | if (isBelow) |
| 339 | 37057 | triangulation.emplace_back(Tetrahedron{ | |
| 340 |
1/2✓ Branch 1 taken 37057 times.
✗ Branch 2 not taken.
|
37057 | triangle[0], triangle[2], triangle[1], midPoint |
| 341 | }); | ||
| 342 | else | ||
| 343 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
36747 | triangulation.emplace_back(Tetrahedron{ |
| 344 |
1/2✓ Branch 1 taken 36747 times.
✗ Branch 2 not taken.
|
36747 | triangle[0], triangle[1], triangle[2], midPoint |
| 345 | }); | ||
| 346 | } | ||
| 347 | 18383 | } | |
| 348 | else | ||
| 349 | { | ||
| 350 | 512 | const auto am = midPoint - pointI; | |
| 351 |
2/2✓ Branch 0 taken 257 times.
✓ Branch 1 taken 255 times.
|
512 | const auto sp = normal*am; |
| 352 | using std::signbit; | ||
| 353 | 512 | const bool isBelow = signbit(sp); | |
| 354 |
2/2✓ Branch 0 taken 257 times.
✓ Branch 1 taken 255 times.
|
512 | if (isBelow) |
| 355 | 257 | triangulation.emplace_back(Tetrahedron{ | |
| 356 |
1/2✓ Branch 1 taken 257 times.
✗ Branch 2 not taken.
|
257 | pointI, points[k], points[j], midPoint |
| 357 | }); | ||
| 358 | else | ||
| 359 | 255 | triangulation.emplace_back(Tetrahedron{ | |
| 360 |
1/2✓ Branch 1 taken 255 times.
✗ Branch 2 not taken.
|
255 | pointI, points[j], points[k], midPoint |
| 361 | }); | ||
| 362 | } | ||
| 363 | } | ||
| 364 | } | ||
| 365 | } | ||
| 366 | } | ||
| 367 | |||
| 368 | // sanity check: if points are not coplanar, then using the mid point policy, we get at least 4 tetrahedrons | ||
| 369 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3121 times.
|
3121 | if (triangulation.size() < 4) |
| 370 |
0/20✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
|
3121 | DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the triangulation!"); |
| 371 | |||
| 372 | 3121 | return triangulation; | |
| 373 | 6242 | } | |
| 374 | |||
| 375 | } // end namespace Dumux | ||
| 376 | |||
| 377 | #endif | ||
| 378 |