| 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 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 | 98884 | 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 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 98884 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 98884 times.
|
197768 | 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 94973 times.
|
98884 | if (convexHullPoints.size() == 3) |
| 127 | 7822 | return std::vector<Triangle>(1, {convexHullPoints[0], convexHullPoints[1], convexHullPoints[2]}); | |
| 128 | |||
| 129 | 94973 | Point midPoint(0.0); | |
| 130 |
4/4✓ Branch 0 taken 382670 times.
✓ Branch 1 taken 94973 times.
✓ Branch 2 taken 382670 times.
✓ Branch 3 taken 94973 times.
|
955286 | for (const auto& p : convexHullPoints) |
| 131 | 382670 | midPoint += p; | |
| 132 | 94973 | midPoint /= convexHullPoints.size(); | |
| 133 | |||
| 134 |
1/2✓ Branch 1 taken 94973 times.
✗ Branch 2 not taken.
|
94973 | std::vector<Triangle> triangulation; |
| 135 |
1/2✓ Branch 1 taken 94973 times.
✗ Branch 2 not taken.
|
94973 | triangulation.reserve(convexHullPoints.size()); |
| 136 | |||
| 137 |
4/4✓ Branch 0 taken 287697 times.
✓ Branch 1 taken 94973 times.
✓ Branch 2 taken 287697 times.
✓ Branch 3 taken 94973 times.
|
765340 | for (std::size_t i = 0; i < convexHullPoints.size()-1; ++i) |
| 138 |
3/6✓ Branch 1 taken 287697 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 287697 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 287697 times.
✗ Branch 8 not taken.
|
863091 | triangulation.emplace_back(Triangle{midPoint, convexHullPoints[i], convexHullPoints[i+1]}); |
| 139 | |||
| 140 |
2/6✓ Branch 1 taken 94973 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 94973 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
189946 | triangulation.emplace_back(Triangle{midPoint, convexHullPoints[convexHullPoints.size()-1], convexHullPoints[0]}); |
| 141 | |||
| 142 | 189946 | return triangulation; | |
| 143 | } | ||
| 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 | 7067 | triangulate(const RandomAccessContainer& points) | |
| 163 | { | ||
| 164 |
1/4✓ Branch 1 taken 7067 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
14134 | const auto convexHullPoints = grahamConvexHull<2>(points); |
| 165 |
1/2✓ Branch 1 taken 7067 times.
✗ Branch 2 not taken.
|
14134 | return triangulate<dim, dimWorld, TriangulationPolicy::MidPointPolicy>(convexHullPoints); |
| 166 | } | ||
| 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 | 1329 | 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 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1329 times.
|
1329 | const auto numPoints = points.size(); |
| 193 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1329 times.
|
1329 | 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 1235 times.
|
1329 | if (numPoints == 4) |
| 197 | 188 | 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 | 1235 | Point midPoint(0.0); | |
| 202 | 1235 | Point lowerLeft(1e100); | |
| 203 | 1235 | Point upperRight(-1e100); | |
| 204 |
4/4✓ Branch 0 taken 9932 times.
✓ Branch 1 taken 1235 times.
✓ Branch 2 taken 9932 times.
✓ Branch 3 taken 1235 times.
|
12402 | for (const auto& p : points) |
| 205 | { | ||
| 206 | midPoint += p; | ||
| 207 |
2/2✓ Branch 0 taken 29796 times.
✓ Branch 1 taken 9932 times.
|
39728 | for (int i = 0; i < dimWorld; ++i) |
| 208 | { | ||
| 209 | using std::max; using std::min; | ||
| 210 |
8/8✓ Branch 0 taken 15698 times.
✓ Branch 1 taken 14098 times.
✓ Branch 2 taken 15698 times.
✓ Branch 3 taken 14098 times.
✓ Branch 4 taken 15698 times.
✓ Branch 5 taken 14098 times.
✓ Branch 6 taken 7146 times.
✓ Branch 7 taken 22650 times.
|
105086 | lowerLeft[i] = min(p[i], lowerLeft[i]); |
| 211 |
6/6✓ Branch 0 taken 7146 times.
✓ Branch 1 taken 22650 times.
✓ Branch 2 taken 7146 times.
✓ Branch 3 taken 22650 times.
✓ Branch 4 taken 7146 times.
✓ Branch 5 taken 22650 times.
|
103680 | upperRight[i] = max(p[i], upperRight[i]); |
| 212 | } | ||
| 213 | } | ||
| 214 | 1235 | midPoint /= numPoints; | |
| 215 | |||
| 216 | 1235 | auto magnitude = 0.0; | |
| 217 | using std::max; | ||
| 218 |
2/2✓ Branch 0 taken 3705 times.
✓ Branch 1 taken 1235 times.
|
4940 | for (int i = 0; i < dimWorld; ++i) |
| 219 |
6/6✓ Branch 0 taken 1491 times.
✓ Branch 1 taken 2214 times.
✓ Branch 2 taken 1491 times.
✓ Branch 3 taken 2214 times.
✓ Branch 4 taken 1491 times.
✓ Branch 5 taken 2214 times.
|
12606 | magnitude = max(upperRight[i] - lowerLeft[i], magnitude); |
| 220 | 1235 | const auto eps = 1e-7*magnitude; | |
| 221 | 1235 | const auto eps2 = eps*magnitude; | |
| 222 | |||
| 223 | // reserve memory conservatively to avoid reallocation | ||
| 224 |
1/4✓ Branch 1 taken 1235 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
1235 | std::vector<Tetrahedron> triangulation; |
| 225 |
1/2✓ Branch 1 taken 1235 times.
✗ Branch 2 not taken.
|
1235 | triangulation.reserve(numPoints); |
| 226 | |||
| 227 | // make a buffer for storing coplanar points and indices | ||
| 228 |
3/4✓ Branch 0 taken 1219 times.
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1235 times.
|
2470 | std::vector<Point> coplanarPointBuffer; |
| 229 |
3/6✓ Branch 0 taken 1219 times.
✓ Branch 1 taken 16 times.
✓ Branch 3 taken 1235 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
2454 | 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 |
2/6✓ Branch 1 taken 1235 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1235 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
2470 | std::vector<std::pair<int, Point>> coplanarClusters; |
| 235 |
1/2✓ Branch 1 taken 1235 times.
✗ Branch 2 not taken.
|
1235 | 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 9932 times.
✓ Branch 1 taken 1235 times.
|
11167 | for (int i = 0; i < numPoints; ++i) |
| 242 | { | ||
| 243 |
2/2✓ Branch 0 taken 35984 times.
✓ Branch 1 taken 9932 times.
|
45916 | for (int j = i+1; j < numPoints; ++j) |
| 244 | { | ||
| 245 |
2/2✓ Branch 0 taken 80273 times.
✓ Branch 1 taken 35984 times.
|
116257 | for (int k = j+1; k < numPoints; ++k) |
| 246 | { | ||
| 247 | 80273 | const auto pointI = points[i]; | |
| 248 | 160546 | const auto ab = points[j] - pointI; | |
| 249 | 160546 | const auto ac = points[k] - pointI; | |
| 250 | 80273 | const auto normal = crossProduct(ab, ac); | |
| 251 | |||
| 252 | // clear list of coplanar points w.r.t to triangle ijk | ||
| 253 |
2/2✓ Branch 0 taken 18919 times.
✓ Branch 1 taken 61354 times.
|
80273 | 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 | 80273 | const bool isAdmissible = [&]() | |
| 258 | { | ||
| 259 | // check if points are colinear and we can't form a triangle | ||
| 260 | // if so, skip this triangle | ||
| 261 | 80273 | 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 408640 times.
✓ Branch 1 taken 7577 times.
|
416217 | for (int m = 0; m < numPoints; ++m) |
| 266 | { | ||
| 267 | 1047795 | if (m != i && m != j && m != k) | |
| 268 | { | ||
| 269 | // check scalar product with surface normal to decide side | ||
| 270 | 568380 | const auto ad = points[m] - pointI; | |
| 271 | 284190 | 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 |
2/2✓ Branch 0 taken 236513 times.
✓ Branch 1 taken 47677 times.
|
284190 | const bool coplanar = abs(sp) < eps2*magnitude; |
| 276 |
4/4✓ Branch 0 taken 116711 times.
✓ Branch 1 taken 119802 times.
✓ Branch 2 taken 116711 times.
✓ Branch 3 taken 119802 times.
|
473026 | 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 71370 times.
✓ Branch 1 taken 212820 times.
|
284190 | if (marker == 0 && newMarker != 0) |
| 281 | 71370 | 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 239689 times.
✓ Branch 1 taken 44501 times.
|
284190 | if (newMarker != 0 && marker != newMarker) |
| 286 | 72696 | return false; | |
| 287 | |||
| 288 | // handle possible coplanar points | ||
| 289 |
2/2✓ Branch 0 taken 47677 times.
✓ Branch 1 taken 192012 times.
|
239689 | if (coplanar) |
| 290 | { | ||
| 291 | using std::abs; | ||
| 292 |
2/2✓ Branch 0 taken 39131 times.
✓ Branch 1 taken 8546 times.
|
47677 | if (m < k && std::find_if( |
| 293 | 39131 | coplanarClusters.begin(), coplanarClusters.end(), | |
| 294 |
8/8✓ Branch 0 taken 78217 times.
✓ Branch 1 taken 63984 times.
✓ Branch 2 taken 54040 times.
✓ Branch 3 taken 88161 times.
✓ Branch 4 taken 25845 times.
✓ Branch 5 taken 28195 times.
✓ Branch 6 taken 25845 times.
✓ Branch 7 taken 28195 times.
|
328498 | [=](const auto& c){ return c.first == std::min(m, i) && abs(ab*c.second) < eps2*magnitude; } |
| 295 |
4/4✓ Branch 2 taken 28195 times.
✓ Branch 3 taken 10936 times.
✓ Branch 4 taken 28195 times.
✓ Branch 5 taken 10936 times.
|
78262 | ) != coplanarClusters.end()) |
| 296 | { | ||
| 297 | // this cluster has already been handled | ||
| 298 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 28183 times.
|
28195 | coplanarPointBuffer.clear(); |
| 299 | 28195 | return false; | |
| 300 | } | ||
| 301 | else | ||
| 302 | 38964 | 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 |
4/4✓ Branch 0 taken 7067 times.
✓ Branch 1 taken 510 times.
✓ Branch 2 taken 7067 times.
✓ Branch 3 taken 510 times.
|
15154 | if (!coplanarPointBuffer.empty()) |
| 310 | { | ||
| 311 | 21201 | coplanarPointBuffer.insert(coplanarPointBuffer.end(), { points[i], points[j], points[k] }); | |
| 312 | 7067 | 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 |
10/12✓ Branch 1 taken 80273 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 80273 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 339789 times.
✓ Branch 6 taken 68851 times.
✓ Branch 7 taken 299366 times.
✓ Branch 8 taken 40423 times.
✓ Branch 9 taken 284190 times.
✓ Branch 10 taken 15176 times.
✓ Branch 11 taken 236513 times.
✓ Branch 12 taken 47677 times.
|
874577 | }(); |
| 321 | |||
| 322 |
2/2✓ Branch 0 taken 7577 times.
✓ Branch 1 taken 72696 times.
|
80273 | 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 |
4/4✓ Branch 0 taken 7067 times.
✓ Branch 1 taken 510 times.
✓ Branch 2 taken 7067 times.
✓ Branch 3 taken 510 times.
|
15154 | if (!coplanarPointBuffer.empty()) |
| 327 | { | ||
| 328 |
2/6✓ Branch 1 taken 7067 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7067 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
21201 | const auto triangles = triangulate<2, 3, TriangulationPolicy::ConvexHullPolicy>(coplanarPointBuffer); |
| 329 |
4/4✓ Branch 0 taken 28540 times.
✓ Branch 1 taken 7067 times.
✓ Branch 2 taken 28540 times.
✓ Branch 3 taken 7067 times.
|
49741 | for (const auto& triangle : triangles) |
| 330 | { | ||
| 331 | 85620 | const auto ab = triangle[1] - triangle[0]; | |
| 332 | 85620 | const auto ac = triangle[2] - triangle[0]; | |
| 333 | 28540 | const auto normal = crossProduct(ab, ac); | |
| 334 | 57080 | const auto am = midPoint - triangle[0]; | |
| 335 | 28540 | const auto sp = normal*am; | |
| 336 | using std::signbit; | ||
| 337 |
2/2✓ Branch 0 taken 14425 times.
✓ Branch 1 taken 14115 times.
|
28540 | const bool isBelow = signbit(sp); |
| 338 |
2/2✓ Branch 0 taken 14425 times.
✓ Branch 1 taken 14115 times.
|
28540 | if (isBelow) |
| 339 | 14425 | triangulation.emplace_back(Tetrahedron{ | |
| 340 |
4/8✓ Branch 1 taken 14425 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14425 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14425 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 14425 times.
✗ Branch 11 not taken.
|
57700 | triangle[0], triangle[2], triangle[1], midPoint |
| 341 | }); | ||
| 342 | else | ||
| 343 | 14115 | triangulation.emplace_back(Tetrahedron{ | |
| 344 |
4/8✓ Branch 1 taken 14115 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14115 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14115 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 14115 times.
✗ Branch 11 not taken.
|
56460 | triangle[0], triangle[1], triangle[2], midPoint |
| 345 | }); | ||
| 346 | } | ||
| 347 | } | ||
| 348 | else | ||
| 349 | { | ||
| 350 | 510 | const auto am = midPoint - pointI; | |
| 351 | 510 | const auto sp = normal*am; | |
| 352 | using std::signbit; | ||
| 353 |
2/2✓ Branch 0 taken 256 times.
✓ Branch 1 taken 254 times.
|
510 | const bool isBelow = signbit(sp); |
| 354 |
2/2✓ Branch 0 taken 256 times.
✓ Branch 1 taken 254 times.
|
510 | if (isBelow) |
| 355 | 256 | triangulation.emplace_back(Tetrahedron{ | |
| 356 |
3/6✓ Branch 1 taken 256 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 256 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 256 times.
✗ Branch 8 not taken.
|
768 | pointI, points[k], points[j], midPoint |
| 357 | }); | ||
| 358 | else | ||
| 359 | 254 | triangulation.emplace_back(Tetrahedron{ | |
| 360 |
3/6✓ Branch 1 taken 254 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 254 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 254 times.
✗ Branch 8 not taken.
|
762 | 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 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1235 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1235 times.
|
2470 | if (triangulation.size() < 4) |
| 370 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the triangulation!"); | |
| 371 | |||
| 372 |
2/4✓ Branch 0 taken 1235 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1235 times.
✗ Branch 3 not taken.
|
2470 | return triangulation; |
| 373 | } | ||
| 374 | |||
| 375 | } // end namespace Dumux | ||
| 376 | |||
| 377 | #endif | ||
| 378 |