OGS
MeshElementGrid.cpp
Go to the documentation of this file.
1 /*
2  * \brief Definition of the class MeshElementGrid.
3  *
4  * \file
6  * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
8  * See accompanying file LICENSE.txt or
9  */
10
11 #include "MeshElementGrid.h"
12
13 #include <algorithm>
14 #include <bitset>
15 #include <cmath>
16 #include <memory>
17
19 #include "MeshLib/Mesh.h"
20 #include "MeshLib/Node.h"
21
22 namespace MeshLib
23 {
25  : _aabb{mesh.getNodes().cbegin(), mesh.getNodes().cend()},
26  _n_steps({{1, 1, 1}})
27 {
28  auto getDimensions = [](auto const& min, auto const& max)
29  {
30  std::bitset<3> dim; // all bits are set to zero.
31  for (std::size_t k(0); k < 3; ++k)
32  {
33  double const tolerance(
34  std::nexttoward(max[k], std::numeric_limits<double>::max()) -
35  max[k]);
36  if (std::abs(max[k] - min[k]) > tolerance)
37  {
38  dim[k] = true;
39  }
40  }
41  return dim;
42  };
43
44  auto const& min_pnt(_aabb.getMinPoint());
45  auto const& max_pnt(_aabb.getMaxPoint());
46  auto const dim = getDimensions(min_pnt, max_pnt);
47
48  std::array<double, 3> delta{{max_pnt[0] - min_pnt[0],
49  max_pnt[1] - min_pnt[1],
50  max_pnt[2] - min_pnt[2]}};
51
52  const std::size_t n_eles(mesh.getNumberOfElements());
53  const std::size_t n_eles_per_cell(100);
54
55  // *** condition: n_eles / n_cells < n_eles_per_cell
56  // where n_cells = _n_steps[0] * _n_steps[1] * _n_steps[2]
57  // *** with _n_steps[0] =
58  // ceil(pow(n_eles*delta[0]*delta[0]/(n_eles_per_cell*delta[1]*delta[2]),
59  // 1/3.)));
60  // _n_steps[1] = _n_steps[0] * delta[1]/delta[0],
61  // _n_steps[2] = _n_steps[0] * delta[2]/delta[0]
62  auto sc_ceil = [](double v)
63  { return static_cast<std::size_t>(std::ceil(v)); };
64
65  switch (dim.count())
66  {
67  case 3: // 3d case
68  _n_steps[0] =
69  sc_ceil(std::cbrt(n_eles * delta[0] * delta[0] /
70  (n_eles_per_cell * delta[1] * delta[2])));
71  _n_steps[1] =
72  sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0));
73  _n_steps[2] =
74  sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0));
75  break;
76  case 2: // 2d cases
77  if (dim[0] && dim[2])
78  { // 2d case: xz plane, y = const
79  _n_steps[0] = sc_ceil(std::sqrt(n_eles * delta[0] /
80  (n_eles_per_cell * delta[2])));
81  _n_steps[2] = sc_ceil(_n_steps[0] * delta[2] / delta[0]);
82  }
83  else if (dim[0] && dim[1])
84  { // 2d case: xy plane, z = const
85  _n_steps[0] = sc_ceil(std::sqrt(n_eles * delta[0] /
86  (n_eles_per_cell * delta[1])));
87  _n_steps[1] = sc_ceil(_n_steps[0] * delta[1] / delta[0]);
88  }
89  else if (dim[1] && dim[2])
90  { // 2d case: yz plane, x = const
91  _n_steps[1] = sc_ceil(std::sqrt(n_eles * delta[1] /
92  (n_eles_per_cell * delta[2])));
93  _n_steps[2] =
94  sc_ceil(n_eles * delta[2] / (n_eles_per_cell * delta[1]));
95  }
96  break;
97  case 1: // 1d cases
98  for (std::size_t k(0); k < 3; ++k)
99  {
100  if (dim[k])
101  {
102  _n_steps[k] =
103  sc_ceil(static_cast<double>(n_eles) / n_eles_per_cell);
104  }
105  }
106  }
107
108  // some frequently used expressions to fill the vector of elements per grid
109  // cell
110  for (std::size_t k(0); k < 3; k++)
111  {
112  _step_sizes[k] = delta[k] / _n_steps[k];
113  _inverse_step_sizes[k] = 1.0 / _step_sizes[k];
114  }
115
116  _elements_in_grid_box.resize(_n_steps[0] * _n_steps[1] * _n_steps[2]);
117  sortElementsInGridCells(mesh);
118 }
119
120 Eigen::Vector3d const& MeshElementGrid::getMinPoint() const
121 {
122  return _aabb.getMinPoint();
123 }
124
125 Eigen::Vector3d const& MeshElementGrid::getMaxPoint() const
126 {
127  return _aabb.getMaxPoint();
128 }
129
131 {
132  for (auto const element : mesh.getElements())
133  {
134  if (!sortElementInGridCells(*element))
135  {
136  OGS_FATAL("Sorting element (id={:d}) into mesh element grid.",
137  element->getID());
138  }
139  }
140 }
141
143 {
144  std::array<std::size_t, 3> min{};
145  std::array<std::size_t, 3> max{};
146  std::pair<bool, std::array<std::size_t, 3>> c(getGridCellCoordinates(
147  *(static_cast<MathLib::Point3d const*>(element.getNode(0)))));
148  if (c.first)
149  {
150  min = c.second;
151  max = min;
152  }
153  else
154  {
155  return false;
156  }
157
158  for (std::size_t k(1); k < element.getNumberOfNodes(); ++k)
159  {
160  // compute coordinates of the grid for each node of the element
162  *(static_cast<MathLib::Point3d const*>(element.getNode(k))));
163  if (!c.first)
164  {
165  return false;
166  }
167
168  for (std::size_t j(0); j < 3; ++j)
169  {
170  if (min[j] > c.second[j])
171  {
172  min[j] = c.second[j];
173  }
174  if (max[j] < c.second[j])
175  {
176  max[j] = c.second[j];
177  }
178  }
179  }
180
181  const std::size_t n_plane(_n_steps[0] * _n_steps[1]);
182
183  // If a node of an element is almost equal to the upper right point of the
184  // AABB the grid cell coordinates computed by getGridCellCoordintes() could
185  // be to large (due to numerical errors). The following lines ensure that
186  // the grid cell coordinates are in the valid range.
187  for (std::size_t k(0); k < 3; ++k)
188  {
189  max[k] = std::min(_n_steps[k] - 1, max[k]);
190  }
191
192  // insert the element into the grid cells
193  for (std::size_t i(min[0]); i <= max[0]; i++)
194  {
195  for (std::size_t j(min[1]); j <= max[1]; j++)
196  {
197  for (std::size_t k(min[2]); k <= max[2]; k++)
198  {
199  _elements_in_grid_box[i + j * _n_steps[0] + k * n_plane]
200  .push_back(&element);
201  }
202  }
203  }
204
205  return true;
206 }
207
208 std::pair<bool, std::array<std::size_t, 3>>
210 {
211  bool valid(true);
212  std::array<std::size_t, 3> coords{};
213
214  for (std::size_t k(0); k < 3; ++k)
215  {
216  const double d(p[k] - _aabb.getMinPoint()[k]);
217  if (d < 0.0)
218  {
219  valid = false;
220  coords[k] = 0;
221  }
222  else if (_aabb.getMaxPoint()[k] <= p[k])
223  {
224  valid = false;
225  coords[k] = _n_steps[k] - 1;
226  }
227  else
228  {
229  coords[k] = static_cast<std::size_t>(d * _inverse_step_sizes[k]);
230  }
231  }
232
233  return std::make_pair(valid, coords);
234 }
235 } // namespace MeshLib
Definition of the Element class.
#define OGS_FATAL(...)
Definition: Error.h:26
std::vector< std::size_t > getNodes(GeoLib::Point const &pnt, std::vector< MeshLib::Node * > const &nodes, MeshLib::PropertyVector< int > const &mat_ids, std::pair< int, int > const &mat_limits, std::pair< double, double > const &elevation_limits, MeshLib::Mesh const &mesh)
Definition of the Mesh class.
Definition of the Node class.
std::pair< std::size_t, std::size_t > getDimensions(std::string const &line)
Returns raster dimensions from the "Zone"-description.
Eigen::Vector3d const & getMinPoint() const
Definition: AABB.h:174
Eigen::Vector3d const & getMaxPoint() const
Definition: AABB.h:181
virtual const Node * getNode(unsigned idx) const =0
virtual unsigned getNumberOfNodes() const =0
void sortElementsInGridCells(MeshLib::Mesh const &mesh)
std::vector< std::vector< MeshLib::Element const * > > _elements_in_grid_box
MeshElementGrid(MeshLib::Mesh const &mesh)
std::pair< bool, std::array< std::size_t, 3 > > getGridCellCoordinates(MathLib::Point3d const &p) const
bool sortElementInGridCells(MeshLib::Element const &element)
Eigen::Vector3d const & getMinPoint() const
std::array< double, 3 > _inverse_step_sizes
std::array< std::size_t, 3 > _n_steps
Eigen::Vector3d const & getMaxPoint() const
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:98
constexpr std::array c
static double const tolerance