OGS
GeoLib::SurfaceGrid Class Referencefinal

Detailed Description

Definition at line 29 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
 
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 25 of file SurfaceGrid.cpp.

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

Member Function Documentation

◆ getGridCellCoordinates()

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

Definition at line 221 of file SurfaceGrid.cpp.

223 {
224  auto const& min_point{getMinPoint()};
225  std::array<std::size_t, 3> coords{
226  {static_cast<std::size_t>((p[0] - min_point[0]) *
228  static_cast<std::size_t>((p[1] - min_point[1]) *
230  static_cast<std::size_t>((p[2] - min_point[2]) *
231  _inverse_step_sizes[2])}};
232 
233  if (coords[0] >= _n_steps[0] || coords[1] >= _n_steps[1] ||
234  coords[2] >= _n_steps[2])
235  {
236  DBUG(
237  "Computed indices ({:d},{:d},{:d}), max grid cell indices "
238  "({:d},{:d},{:d})",
239  coords[0], coords[1], coords[2], _n_steps[0], _n_steps[1],
240  _n_steps[2]);
241  return std::optional<std::array<std::size_t, 3>>();
242  }
243  return std::optional<std::array<std::size_t, 3>>(coords);
244 }
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:27

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 246 of file SurfaceGrid.cpp.

248 {
249  std::optional<std::array<std::size_t, 3>> optional_c(
251  if (!optional_c)
252  {
253  return false;
254  }
255  std::array<std::size_t, 3> c(optional_c.value());
256 
257  std::size_t const grid_cell_idx(c[0] + c[1] * _n_steps[0] +
258  c[2] * _n_steps[0] * _n_steps[1]);
259  std::vector<Triangle const*> const& triangles(
260  _triangles_in_grid_box[grid_cell_idx]);
261  auto const it = std::find_if(triangles.begin(), triangles.end(),
262  [eps, pnt](auto const* triangle)
263  { return triangle->containsPoint(pnt, eps); });
264  return it != triangles.end();
265 }
std::optional< std::array< std::size_t, 3 > > getGridCellCoordinates(MathLib::Point3d const &p) const

References _n_steps, _triangles_in_grid_box, MaterialPropertyLib::c, and getGridCellCoordinates().

◆ sortTriangleInGridCells()

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

Definition at line 167 of file SurfaceGrid.cpp.

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

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 146 of file SurfaceGrid.cpp.

147 {
148  for (std::size_t l(0); l < sfc->getNumberOfTriangles(); l++)
149  {
150  if (!sortTriangleInGridCells((*sfc)[l]))
151  {
152  Point const& p0(*((*sfc)[l]->getPoint(0)));
153  Point const& p1(*((*sfc)[l]->getPoint(1)));
154  Point const& p2(*((*sfc)[l]->getPoint(2)));
155  auto const& min{getMinPoint()};
156  auto const& max{getMaxPoint()};
157  OGS_FATAL(
158  "Sorting triangle {:d} [({:f},{:f},{:f}), ({:f},{:f},{:f}), "
159  "({:f},{:f},{:f}) into grid. Bounding box is [{:f}, {:f}] x "
160  "[{:f}, {:f}] x [{:f}, {:f}].",
161  l, p0[0], p0[1], p0[2], p1[0], p1[1], p1[2], p2[0], p2[1],
162  p2[2], min[0], max[0], min[1], max[1], min[2], max[2]);
163  }
164  }
165 }
#define OGS_FATAL(...)
Definition: Error.h:26
bool sortTriangleInGridCells(GeoLib::Triangle const *const triangle)
TemplateElement< PointRule1 > Point
Definition: Point.h:20

References GeoLib::AABB::getMaxPoint(), GeoLib::AABB::getMinPoint(), 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 41 of file SurfaceGrid.h.

Referenced by getGridCellCoordinates().

◆ _n_steps

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

Definition at line 42 of file SurfaceGrid.h.

Referenced by getGridCellCoordinates(), isPointInSurface(), and sortTriangleInGridCells().

◆ _step_sizes

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

Definition at line 40 of file SurfaceGrid.h.

◆ _triangles_in_grid_box

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

Definition at line 43 of file SurfaceGrid.h.

Referenced by isPointInSurface(), and sortTriangleInGridCells().


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