OGS
GeoLib::SurfaceGrid Class Referencefinal

Detailed Description

Definition at line 21 of file SurfaceGrid.h.

#include <SurfaceGrid.h>

Inheritance diagram for GeoLib::SurfaceGrid:
[legend]
Collaboration diagram for GeoLib::SurfaceGrid:
[legend]

Public Member Functions

 SurfaceGrid (GeoLib::Surface const *const sfc)
bool isPointInSurface (MathLib::Point3d const &pnt, double eps=std::numeric_limits< double >::epsilon()) const
Public Member Functions inherited from GeoLib::AABB
template<typename PNT_TYPE>
 AABB (std::vector< PNT_TYPE * > const &pnts, std::vector< std::size_t > const &ids)
template<typename InputIterator>
 AABB (InputIterator first, InputIterator last)
template<typename PNT_TYPE>
bool update (PNT_TYPE const &p)
template<typename T>
bool containsPoint (T const &pnt, double eps) const
template<typename T>
bool containsPointXY (T const &pnt) const
MinMaxPoints getMinMaxPoints () const
Eigen::Vector3d const & getMinPoint () const
Eigen::Vector3d const & getMaxPoint () const
bool containsAABB (AABB const &other_aabb) const

Private Member Functions

void sortTrianglesInGridCells (GeoLib::Surface const *const sfc)
bool sortTriangleInGridCells (GeoLib::Triangle const *const triangle)
std::optional< std::array< std::size_t, 3 > > getGridCellCoordinates (MathLib::Point3d const &p) const

Private Attributes

std::array< double, 3 > _step_sizes {}
std::array< double, 3 > _inverse_step_sizes {}
std::array< std::size_t, 3 > _n_steps
std::vector< std::vector< GeoLib::Triangle const * > > _triangles_in_grid_box

Constructor & Destructor Documentation

◆ SurfaceGrid()

GeoLib::SurfaceGrid::SurfaceGrid ( GeoLib::Surface const *const sfc)
explicit

Definition at line 17 of file SurfaceGrid.cpp.

18 : AABB(sfc->getAABB()), _n_steps({{1, 1, 1}})
19{
20 auto [min_point, max_point] = getMinMaxPoints();
21 // enlarge the bounding, such that the points with maximal coordinates
22 // fits into the grid
23 for (std::size_t k(0); k < 3; ++k)
24 {
25 max_point[k] += std::abs(max_point[k]) * 1e-6;
26 if (std::abs(max_point[k]) < std::numeric_limits<double>::epsilon())
27 {
28 max_point[k] = (max_point[k] - min_point[k]) * (1.0 + 1e-6);
29 }
30 }
31
32 Eigen::Vector3d delta = max_point - min_point;
33
34 if (delta[0] < std::numeric_limits<double>::epsilon())
35 {
36 const double max_delta(std::max(delta[1], delta[2]));
37 min_point[0] -= max_delta * 0.5e-3;
38 max_point[0] += max_delta * 0.5e-3;
39 delta[0] = max_point[0] - min_point[0];
40 }
41
42 if (delta[1] < std::numeric_limits<double>::epsilon())
43 {
44 const double max_delta(std::max(delta[0], delta[2]));
45 min_point[1] -= max_delta * 0.5e-3;
46 max_point[1] += max_delta * 0.5e-3;
47 delta[1] = max_point[1] - min_point[1];
48 }
49
50 if (delta[2] < std::numeric_limits<double>::epsilon())
51 {
52 const double max_delta(std::max(delta[0], delta[1]));
53 min_point[2] -= max_delta * 0.5e-3;
54 max_point[2] += max_delta * 0.5e-3;
55 delta[2] = max_point[2] - min_point[2];
56 }
57
58 update(min_point);
59 update(max_point);
60
61 const std::size_t n_tris(sfc->getNumberOfTriangles());
62 const std::size_t n_tris_per_cell(5);
63
64 Eigen::Matrix<bool, 3, 1> dim =
65 delta.array() >= std::numeric_limits<double>::epsilon();
66
67 // *** condition: n_tris / n_cells < n_tris_per_cell
68 // where n_cells = _n_steps[0] * _n_steps[1] * _n_steps[2]
69 // *** with _n_steps[0] =
70 // ceil(pow(n_tris*delta[0]*delta[0]/(n_tris_per_cell*delta[1]*delta[2]),
71 // 1/3.)));
72 // _n_steps[1] = _n_steps[0] * delta[1]/delta[0],
73 // _n_steps[2] = _n_steps[0] * delta[2]/delta[0]
74 auto sc_ceil = [](double v)
75 { return static_cast<std::size_t>(std::ceil(v)); };
76 switch (dim.count())
77 {
78 case 3: // 3d case
79 _n_steps[0] =
80 sc_ceil(std::cbrt(n_tris * delta[0] * delta[0] /
81 (n_tris_per_cell * delta[1] * delta[2])));
82 _n_steps[1] =
83 sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0));
84 _n_steps[2] =
85 sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0));
86 break;
87 case 2: // 2d cases
88 if (dim[0] && dim[2])
89 { // 2d case: xz plane, y = const
90 _n_steps[0] = sc_ceil(std::sqrt(n_tris * delta[0] /
91 (n_tris_per_cell * delta[2])));
92 _n_steps[2] = sc_ceil(_n_steps[0] * delta[2] / delta[0]);
93 }
94 else if (dim[0] && dim[1])
95 { // 2d case: xy plane, z = const
96 _n_steps[0] = sc_ceil(std::sqrt(n_tris * delta[0] /
97 (n_tris_per_cell * delta[1])));
98 _n_steps[1] = sc_ceil(_n_steps[0] * delta[1] / delta[0]);
99 }
100 else if (dim[1] && dim[2])
101 { // 2d case: yz plane, x = const
102 _n_steps[1] = sc_ceil(std::sqrt(n_tris * delta[1] /
103 (n_tris_per_cell * delta[2])));
104 _n_steps[2] =
105 sc_ceil(n_tris * delta[2] / (n_tris_per_cell * delta[1]));
106 }
107 break;
108 case 1: // 1d cases
109 for (std::size_t k(0); k < 3; ++k)
110 {
111 if (dim[k])
112 {
113 _n_steps[k] =
114 sc_ceil(static_cast<double>(n_tris) / n_tris_per_cell);
115 }
116 }
117 }
118
119 // some frequently used expressions to fill the grid vectors
120 for (std::size_t k(0); k < 3; k++)
121 {
122 _step_sizes[k] = delta[k] / _n_steps[k];
123 if (delta[k] > std::numeric_limits<double>::epsilon())
124 {
125 _inverse_step_sizes[k] = 1.0 / _step_sizes[k];
126 }
127 else
128 {
129 _inverse_step_sizes[k] = 0;
130 }
131 }
132
135}
AABB(std::vector< PNT_TYPE * > const &pnts, std::vector< std::size_t > const &ids)
Definition AABB.h:53
MinMaxPoints getMinMaxPoints() const
Definition AABB.h:163
bool update(PNT_TYPE const &p)
Definition AABB.h:97
void sortTrianglesInGridCells(GeoLib::Surface const *const sfc)
std::vector< std::vector< GeoLib::Triangle const * > > _triangles_in_grid_box
Definition SurfaceGrid.h:37
std::array< std::size_t, 3 > _n_steps
Definition SurfaceGrid.h:36
std::array< double, 3 > _inverse_step_sizes
Definition SurfaceGrid.h:35
std::array< double, 3 > _step_sizes
Definition SurfaceGrid.h:34

References GeoLib::AABB::AABB(), and _n_steps.

Member Function Documentation

◆ getGridCellCoordinates()

std::optional< std::array< std::size_t, 3 > > GeoLib::SurfaceGrid::getGridCellCoordinates ( MathLib::Point3d const & p) const
private

Definition at line 211 of file SurfaceGrid.cpp.

213{
214 auto const& min_point{getMinPoint()};
215 std::array<std::size_t, 3> coords{
216 {static_cast<std::size_t>((p[0] - min_point[0]) *
218 static_cast<std::size_t>((p[1] - min_point[1]) *
220 static_cast<std::size_t>((p[2] - min_point[2]) *
222
223 if (coords[0] >= _n_steps[0] || coords[1] >= _n_steps[1] ||
224 coords[2] >= _n_steps[2])
225 {
226 DBUG(
227 "Computed indices ({:d},{:d},{:d}), max grid cell indices "
228 "({:d},{:d},{:d})",
229 coords[0], coords[1], coords[2], _n_steps[0], _n_steps[1],
230 _n_steps[2]);
231 return {};
232 }
233 return coords;
234}
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:22
Eigen::Vector3d const & getMinPoint() const
Definition AABB.h:169
constexpr ranges::views::view_closure coords
Definition Mesh.h:223

References _inverse_step_sizes, _n_steps, DBUG(), and GeoLib::AABB::getMinPoint().

Referenced by isPointInSurface(), and sortTriangleInGridCells().

◆ isPointInSurface()

bool GeoLib::SurfaceGrid::isPointInSurface ( MathLib::Point3d const & pnt,
double eps = std::numeric_limits<double>::epsilon() ) const

Definition at line 236 of file SurfaceGrid.cpp.

238{
239 auto const grid_cell_coordinates = getGridCellCoordinates(pnt);
240 if (!grid_cell_coordinates)
241 {
242 return false;
243 }
244
245 std::size_t const grid_cell_index =
246 (*grid_cell_coordinates)[0] +
247 (*grid_cell_coordinates)[1] * _n_steps[0] +
248 (*grid_cell_coordinates)[2] * _n_steps[0] * _n_steps[1];
249
250 std::vector<Triangle const*> const& triangles(
251 _triangles_in_grid_box[grid_cell_index]);
252 return std::any_of(triangles.begin(), triangles.end(),
253 [eps, pnt](auto const* triangle)
254 { return triangle->containsPoint(pnt, eps); });
255}
std::optional< std::array< std::size_t, 3 > > getGridCellCoordinates(MathLib::Point3d const &p) const

References _n_steps, _triangles_in_grid_box, and getGridCellCoordinates().

◆ sortTriangleInGridCells()

bool GeoLib::SurfaceGrid::sortTriangleInGridCells ( GeoLib::Triangle const *const triangle)
private

Definition at line 157 of file SurfaceGrid.cpp.

158{
159 // compute grid coordinates for each triangle point
160 std::optional<std::array<std::size_t, 3> const> c_p0(
161 getGridCellCoordinates(*(triangle->getPoint(0))));
162 if (!c_p0)
163 {
164 return false;
165 }
166 std::optional<std::array<std::size_t, 3> const> c_p1(
167 getGridCellCoordinates(*(triangle->getPoint(1))));
168 if (!c_p1)
169 {
170 return false;
171 }
172 std::optional<std::array<std::size_t, 3> const> c_p2(
173 getGridCellCoordinates(*(triangle->getPoint(2))));
174 if (!c_p2)
175 {
176 return false;
177 }
178
179 // determine interval in grid (grid cells the triangle will be inserted)
180 std::size_t const i_min(
181 std::min(std::min((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0]));
182 std::size_t const i_max(
183 std::max(std::max((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0]));
184 std::size_t const j_min(
185 std::min(std::min((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1]));
186 std::size_t const j_max(
187 std::max(std::max((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1]));
188 std::size_t const k_min(
189 std::min(std::min((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2]));
190 std::size_t const k_max(
191 std::max(std::max((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2]));
192
193 const std::size_t n_plane(_n_steps[0] * _n_steps[1]);
194
195 // insert the triangle into the grid cells
196 for (std::size_t i(i_min); i <= i_max; i++)
197 {
198 for (std::size_t j(j_min); j <= j_max; j++)
199 {
200 for (std::size_t k(k_min); k <= k_max; k++)
201 {
202 _triangles_in_grid_box[i + j * _n_steps[0] + k * n_plane]
203 .push_back(triangle);
204 }
205 }
206 }
207
208 return true;
209}

References _n_steps, _triangles_in_grid_box, getGridCellCoordinates(), and GeoLib::Triangle::getPoint().

Referenced by sortTrianglesInGridCells().

◆ sortTrianglesInGridCells()

void GeoLib::SurfaceGrid::sortTrianglesInGridCells ( GeoLib::Surface const *const sfc)
private

Definition at line 137 of file SurfaceGrid.cpp.

138{
139 for (std::size_t l(0); l < sfc->getNumberOfTriangles(); l++)
140 {
141 if (!sortTriangleInGridCells((*sfc)[l]))
142 {
143 Point const& p0(*((*sfc)[l]->getPoint(0)));
144 Point const& p1(*((*sfc)[l]->getPoint(1)));
145 Point const& p2(*((*sfc)[l]->getPoint(2)));
146 auto const [min, max] = getMinMaxPoints();
147 OGS_FATAL(
148 "Sorting triangle {:d} [({:f},{:f},{:f}), ({:f},{:f},{:f}), "
149 "({:f},{:f},{:f}) into grid. Bounding box is [{:f}, {:f}] x "
150 "[{:f}, {:f}] x [{:f}, {:f}].",
151 l, p0[0], p0[1], p0[2], p1[0], p1[1], p1[2], p2[0], p2[1],
152 p2[2], min[0], max[0], min[1], max[1], min[2], max[2]);
153 }
154 }
155}
#define OGS_FATAL(...)
Definition Error.h:19
bool sortTriangleInGridCells(GeoLib::Triangle const *const triangle)
TemplateElement< PointRule1 > Point

References GeoLib::AABB::getMinMaxPoints(), GeoLib::Surface::getNumberOfTriangles(), OGS_FATAL, and sortTriangleInGridCells().

Member Data Documentation

◆ _inverse_step_sizes

std::array<double, 3> GeoLib::SurfaceGrid::_inverse_step_sizes {}
private

Definition at line 35 of file SurfaceGrid.h.

35{};

Referenced by getGridCellCoordinates().

◆ _n_steps

std::array<std::size_t, 3> GeoLib::SurfaceGrid::_n_steps
private

◆ _step_sizes

std::array<double, 3> GeoLib::SurfaceGrid::_step_sizes {}
private

Definition at line 34 of file SurfaceGrid.h.

34{};

◆ _triangles_in_grid_box

std::vector<std::vector<GeoLib::Triangle const*> > GeoLib::SurfaceGrid::_triangles_in_grid_box
private

Definition at line 37 of file SurfaceGrid.h.

Referenced by isPointInSurface(), and sortTriangleInGridCells().


The documentation for this class was generated from the following files: