OGS 6.2.1-97-g73d1aeda3
MeshElementGrid.cpp
Go to the documentation of this file.
1 /*
2  * \brief Definition of the class MeshElementGrid.
3  *
4  * \copyright
5  * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
6  * Distributed under a Modified BSD License.
7  * See accompanying file LICENSE.txt or
8  */
9 
10 #include "MeshElementGrid.h"
11 
12 #include <algorithm>
13 #include <bitset>
14 #include <cmath>
15 #include <memory>
16 
17 #include <logog/include/logog.hpp>
18 
19 #include "../Mesh.h"
20 #include "../Node.h"
21 #include "../Elements/Element.h"
22 
23 #include "GeoLib/GEOObjects.h"
24 
25 namespace MeshLib {
26 
28  _aabb{sfc_mesh.getNodes().cbegin(), sfc_mesh.getNodes().cend()},
29  _n_steps({{1,1,1}})
30 {
31  auto getDimensions =
32  [](MathLib::Point3d const& min, MathLib::Point3d const& max)
33  {
34  std::bitset<3> dim; // all bits are set to zero.
35  for (std::size_t k(0); k < 3; ++k) {
36  double const tolerance(
37  std::nexttoward(max[k],std::numeric_limits<double>::max())-max[k]);
38  if (std::abs(max[k] - min[k]) > tolerance)
39  {
40  dim[k] = true;
41  }
42  }
43  return dim;
44  };
45 
46  MathLib::Point3d const& min_pnt(_aabb.getMinPoint());
47  MathLib::Point3d const& max_pnt(_aabb.getMaxPoint());
48  auto const dim = getDimensions(min_pnt, max_pnt);
49 
50  std::array<double, 3> delta{{ max_pnt[0] - min_pnt[0],
51  max_pnt[1] - min_pnt[1], max_pnt[2] - min_pnt[2] }};
52 
53  const std::size_t n_eles(sfc_mesh.getNumberOfElements());
54  const std::size_t n_eles_per_cell(100);
55 
56  // *** condition: n_eles / n_cells < n_eles_per_cell
57  // where n_cells = _n_steps[0] * _n_steps[1] * _n_steps[2]
58  // *** with _n_steps[0] = ceil(pow(n_eles*delta[0]*delta[0]/(n_eles_per_cell*delta[1]*delta[2]), 1/3.)));
59  // _n_steps[1] = _n_steps[0] * delta[1]/delta[0],
60  // _n_steps[2] = _n_steps[0] * delta[2]/delta[0]
61  auto sc_ceil = [](double v){
62  return static_cast<std::size_t>(std::ceil(v));
63  };
64 
65  switch (dim.count()) {
66  case 3: // 3d case
67  _n_steps[0] = sc_ceil(std::cbrt(
68  n_eles*delta[0]*delta[0]/(n_eles_per_cell*delta[1]*delta[2])));
69  _n_steps[1] = sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0));
70  _n_steps[2] = sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0));
71  break;
72  case 2: // 2d cases
73  if (dim[0] && dim[2]) { // 2d case: xz plane, y = const
74  _n_steps[0] = sc_ceil(std::sqrt(n_eles*delta[0]/(n_eles_per_cell*delta[2])));
75  _n_steps[2] = sc_ceil(_n_steps[0]*delta[2]/delta[0]);
76  }
77  else if (dim[0] && dim[1]) { // 2d case: xy plane, z = const
78  _n_steps[0] = sc_ceil(std::sqrt(n_eles*delta[0]/(n_eles_per_cell*delta[1])));
79  _n_steps[1] = sc_ceil(_n_steps[0] * delta[1]/delta[0]);
80  }
81  else if (dim[1] && dim[2]) { // 2d case: yz plane, x = const
82  _n_steps[1] = sc_ceil(std::sqrt(n_eles*delta[1]/(n_eles_per_cell*delta[2])));
83  _n_steps[2] = sc_ceil(n_eles * delta[2] / (n_eles_per_cell*delta[1]));
84  }
85  break;
86  case 1: // 1d cases
87  for (std::size_t k(0); k<3; ++k) {
88  if (dim[k]) {
89  _n_steps[k] = sc_ceil(static_cast<double>(n_eles)/n_eles_per_cell);
90  }
91  }
92  }
93 
94  // some frequently used expressions to fill the vector of elements per grid
95  // cell
96  for (std::size_t k(0); k<3; k++) {
97  _step_sizes[k] = delta[k] / _n_steps[k];
98  _inverse_step_sizes[k] = 1.0 / _step_sizes[k];
99  }
100 
102  sortElementsInGridCells(sfc_mesh);
103 }
104 
106 {
107  return _aabb.getMinPoint();
108 }
109 
111 {
112  return _aabb.getMaxPoint();
113 }
114 
116 {
117  for (auto const element : sfc_mesh.getElements()) {
118  if (! sortElementInGridCells(*element)) {
119  OGS_FATAL("Sorting element (id=%d) into mesh element grid.",
120  element->getID());
121  }
122  }
123 }
124 
126 {
127  std::array<std::size_t,3> min;
128  std::array<std::size_t,3> max;
129  std::pair<bool, std::array<std::size_t, 3>> c(
130  getGridCellCoordinates(*(static_cast<MathLib::Point3d const*>(element.getNode(0)))));
131  if (c.first) {
132  min = c.second;
133  max = min;
134  } else {
135  return false;
136  }
137 
138  std::vector<std::array<std::size_t,3>> coord_vecs(element.getNumberOfNodes());
139  for (std::size_t k(1); k<element.getNumberOfNodes(); ++k) {
140  // compute coordinates of the grid for each node of the element
141  c = getGridCellCoordinates(*(static_cast<MathLib::Point3d const*>(element.getNode(k))));
142  if (!c.first)
143  {
144  return false;
145  }
146 
147  for (std::size_t j(0); j < 3; ++j)
148  {
149  if (min[j] > c.second[j])
150  {
151  min[j] = c.second[j];
152  }
153  if (max[j] < c.second[j])
154  {
155  max[j] = c.second[j];
156  }
157  }
158  }
159 
160  const std::size_t n_plane(_n_steps[0]*_n_steps[1]);
161 
162  // If a node of an element is almost equal to the upper right point of the
163  // AABB the grid cell coordinates computed by getGridCellCoordintes() could
164  // be to large (due to numerical errors). The following lines ensure that
165  // the grid cell coordinates are in the valid range.
166  for (std::size_t k(0); k < 3; ++k)
167  {
168  max[k] = std::min(_n_steps[k] - 1, max[k]);
169  }
170 
171  // insert the element into the grid cells
172  for (std::size_t i(min[0]); i<=max[0]; i++) {
173  for (std::size_t j(min[1]); j<=max[1]; j++) {
174  for (std::size_t k(min[2]); k<=max[2]; k++) {
175  _elements_in_grid_box[i+j*_n_steps[0]+k*n_plane].push_back(&element);
176  }
177  }
178  }
179 
180  return true;
181 }
182 
183 std::pair<bool, std::array<std::size_t, 3>>
185 {
186  bool valid(true);
187  std::array<std::size_t, 3> coords;
188 
189  for (std::size_t k(0); k<3; ++k) {
190  const double d(p[k]-_aabb.getMinPoint()[k]);
191  if (d < 0.0) {
192  valid = false;
193  coords[k] = 0;
194  } else if (_aabb.getMaxPoint()[k] <= p[k]) {
195  valid = false;
196  coords[k] = _n_steps[k]-1;
197  } else {
198  coords[k] = static_cast<std::size_t>(d * _inverse_step_sizes[k]);
199  }
200  }
201 
202  return std::make_pair(valid, coords);
203 }
204 
205 #ifndef NDEBUG
207  GeoLib::GEOObjects& geometries,
208  std::string& geometry_name)
209 {
210  std::vector<std::string> cell_names;
211 
212  auto addPoints = [&geometries](
213  MathLib::Point3d const& p, std::array<double, 3> const& d,
214  std::array<std::size_t, 3> const& c, std::string& name) {
215  auto pnts = std::make_unique<std::vector<GeoLib::Point*>>();
216  pnts->push_back(new GeoLib::Point(p[0]+c[0]*d[0], p[1]+c[1]*d[1], p[2]+c[2]*d[2]));
217  pnts->push_back(new GeoLib::Point(p[0]+c[0]*d[0], p[1]+(c[1]+1)*d[1], p[2]+c[2]*d[2]));
218  pnts->push_back(new GeoLib::Point(p[0]+(c[0]+1)*d[0], p[1]+(c[1]+1)*d[1], p[2]+c[2]*d[2]));
219  pnts->push_back(new GeoLib::Point(p[0]+(c[0]+1)*d[0], p[1]+c[1]*d[1], p[2]+c[2]*d[2]));
220  pnts->push_back(new GeoLib::Point(p[0]+c[0]*d[0], p[1]+c[1]*d[1], p[2]+(c[2]+1)*d[2]));
221  pnts->push_back(new GeoLib::Point(p[0]+c[0]*d[0], p[1]+(c[1]+1)*d[1], p[2]+(c[2]+1)*d[2]));
222  pnts->push_back(new GeoLib::Point(p[0]+(c[0]+1)*d[0], p[1]+(c[1]+1)*d[1], p[2]+(c[2]+1)*d[2]));
223  pnts->push_back(new GeoLib::Point(p[0]+(c[0]+1)*d[0], p[1]+c[1]*d[1], p[2]+(c[2]+1)*d[2]));
224  std::array<double,3> ulps; // unit in the last place
225  double const towards(std::numeric_limits<double>::max());
226  ulps[0] = std::nextafter(d[0], towards)-d[0];
227  ulps[1] = std::nextafter(d[1], towards)-d[1];
228  ulps[2] = std::nextafter(d[2], towards)-d[2];
229  double const tolerance(std::min(std::min(ulps[0], ulps[1]), ulps[2]));
230  geometries.addPointVec(std::move(pnts), name, nullptr, tolerance);
231  };
232 
233  for (std::size_t i(0); i<grid._n_steps[0]; ++i) {
234  for (std::size_t j(0); j<grid._n_steps[1]; ++j) {
235  for (std::size_t k(0); k<grid._n_steps[2]; ++k) {
236  cell_names.emplace_back("Grid-"+std::to_string(i)+"-"
237  +std::to_string(j)+"-"+std::to_string(k));
238  addPoints(grid._aabb.getMinPoint(), grid._step_sizes,
239  {{i, j, k}}, cell_names.back());
240  auto plys = std::make_unique<std::vector<GeoLib::Polyline*>>();
241  auto & points = *geometries.getPointVec(cell_names.back());
242 
243  auto* ply_bottom(new GeoLib::Polyline(points));
244  for (std::size_t l(0); l < 4; ++l)
245  {
246  ply_bottom->addPoint(l);
247  }
248  ply_bottom->addPoint(0); // close to bottom surface
249  plys->push_back(ply_bottom);
250 
251  auto* ply_top(new GeoLib::Polyline(points));
252  for (std::size_t l(4); l < 8; ++l)
253  {
254  ply_top->addPoint(l);
255  }
256  ply_top->addPoint(4); // close to top surface
257  plys->push_back(ply_top);
258 
259  auto* ply_04(new GeoLib::Polyline(points));
260  ply_04->addPoint(0);
261  ply_04->addPoint(4);
262  plys->push_back(ply_04);
263 
264  auto* ply_15(new GeoLib::Polyline(points));
265  ply_15->addPoint(1);
266  ply_15->addPoint(5);
267  plys->push_back(ply_15);
268 
269  auto* ply_26(new GeoLib::Polyline(points));
270  ply_26->addPoint(2);
271  ply_26->addPoint(6);
272  plys->push_back(ply_26);
273 
274  auto* ply_37(new GeoLib::Polyline(points));
275  ply_37->addPoint(3);
276  ply_37->addPoint(7);
277  plys->push_back(ply_37);
278 
279  geometries.addPolylineVec(std::move(plys), cell_names.back(),
280  nullptr);
281  }
282  }
283  }
284  if (geometries.mergeGeometries(cell_names, geometry_name) == 2)
285  {
286  geometry_name = cell_names.front();
287  }
288 }
289 #endif // NDEBUG
290 
291 } // end namespace MeshLib
MathLib::Point3d const & getMaxPoint() const
Definition: AABB.h:173
Container class for geometric objects.
Definition: GEOObjects.h:62
std::array< std::size_t, 3 > _n_steps
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition: Mesh.h:105
std::array< double, 3 > _inverse_step_sizes
void addPointVec(std::unique_ptr< std::vector< Point *>> points, std::string &name, std::unique_ptr< std::map< std::string, std::size_t >> pnt_id_name_map=nullptr, double eps=sqrt(std::numeric_limits< double >::epsilon()))
Definition: GEOObjects.cpp:45
MathLib::Point3d const & getMinPoint() const
Definition: AABB.h:166
std::vector< std::vector< MeshLib::Element const * > > _elements_in_grid_box
void addPolylineVec(std::unique_ptr< std::vector< Polyline *>> lines, const std::string &name, std::unique_ptr< std::map< std::string, std::size_t >> ply_names=nullptr)
Definition: GEOObjects.cpp:137
MathLib::Point3d const & getMaxPoint() const
void sortElementsInGridCells(MeshLib::Mesh const &sfc_mesh)
Interface for heuristic search length strategy.
Definition: ProjectData.h:30
const std::vector< Point * > * getPointVec(const std::string &name) const
Definition: GEOObjects.cpp:63
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
Definition: Polyline.h:50
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:108
std::pair< std::size_t, std::size_t > getDimensions(std::string const &line)
Returns raster dimensions from the "Zone"-description.
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
virtual unsigned getNumberOfNodes() const =0
Returns the number of all nodes including both linear and nonlinear nodes.
std::array< double, 3 > _step_sizes
bool sortElementInGridCells(MeshLib::Element const &element)
MathLib::Point3d const & getMinPoint() const
std::pair< bool, std::array< std::size_t, 3 > > getGridCellCoordinates(MathLib::Point3d const &p) const
friend void getGridGeometry(MeshLib::MeshElementGrid const &grid, GeoLib::GEOObjects &geometries, std::string &geometry_name)
int mergeGeometries(std::vector< std::string > const &names, std::string &merged_geo_name)
Definition: GEOObjects.cpp:443
MeshElementGrid(MeshLib::Mesh const &mesh)
const Node * getNode(unsigned i) const
Definition: Element.cpp:158