OGS 6.1.0-1721-g6382411ad
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  dim[k] = true;
40  }
41  return dim;
42  };
43 
44  MathLib::Point3d const& min_pnt(_aabb.getMinPoint());
45  MathLib::Point3d 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], max_pnt[2] - min_pnt[2] }};
50 
51  const std::size_t n_eles(sfc_mesh.getNumberOfElements());
52  const std::size_t n_eles_per_cell(100);
53 
54  // *** condition: n_eles / n_cells < n_eles_per_cell
55  // where n_cells = _n_steps[0] * _n_steps[1] * _n_steps[2]
56  // *** with _n_steps[0] = ceil(pow(n_eles*delta[0]*delta[0]/(n_eles_per_cell*delta[1]*delta[2]), 1/3.)));
57  // _n_steps[1] = _n_steps[0] * delta[1]/delta[0],
58  // _n_steps[2] = _n_steps[0] * delta[2]/delta[0]
59  auto sc_ceil = [](double v){
60  return static_cast<std::size_t>(std::ceil(v));
61  };
62 
63  switch (dim.count()) {
64  case 3: // 3d case
65  _n_steps[0] = sc_ceil(std::cbrt(
66  n_eles*delta[0]*delta[0]/(n_eles_per_cell*delta[1]*delta[2])));
67  _n_steps[1] = sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0));
68  _n_steps[2] = sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0));
69  break;
70  case 2: // 2d cases
71  if (dim[0] && dim[2]) { // 2d case: xz plane, y = const
72  _n_steps[0] = sc_ceil(std::sqrt(n_eles*delta[0]/(n_eles_per_cell*delta[2])));
73  _n_steps[2] = sc_ceil(_n_steps[0]*delta[2]/delta[0]);
74  }
75  else if (dim[0] && dim[1]) { // 2d case: xy plane, z = const
76  _n_steps[0] = sc_ceil(std::sqrt(n_eles*delta[0]/(n_eles_per_cell*delta[1])));
77  _n_steps[1] = sc_ceil(_n_steps[0] * delta[1]/delta[0]);
78  }
79  else if (dim[1] && dim[2]) { // 2d case: yz plane, x = const
80  _n_steps[1] = sc_ceil(std::sqrt(n_eles*delta[1]/(n_eles_per_cell*delta[2])));
81  _n_steps[2] = sc_ceil(n_eles * delta[2] / (n_eles_per_cell*delta[1]));
82  }
83  break;
84  case 1: // 1d cases
85  for (std::size_t k(0); k<3; ++k) {
86  if (dim[k]) {
87  _n_steps[k] = sc_ceil(static_cast<double>(n_eles)/n_eles_per_cell);
88  }
89  }
90  }
91 
92  // some frequently used expressions to fill the vector of elements per grid
93  // cell
94  for (std::size_t k(0); k<3; k++) {
95  _step_sizes[k] = delta[k] / _n_steps[k];
96  _inverse_step_sizes[k] = 1.0 / _step_sizes[k];
97  }
98 
100  sortElementsInGridCells(sfc_mesh);
101 }
102 
104 {
105  return _aabb.getMinPoint();
106 }
107 
109 {
110  return _aabb.getMaxPoint();
111 }
112 
114 {
115  for (auto const element : sfc_mesh.getElements()) {
116  if (! sortElementInGridCells(*element)) {
117  OGS_FATAL("Sorting element (id=%d) into mesh element grid.",
118  element->getID());
119  }
120  }
121 }
122 
124 {
125  std::array<std::size_t,3> min;
126  std::array<std::size_t,3> max;
127  std::pair<bool, std::array<std::size_t, 3>> c(
128  getGridCellCoordinates(*(static_cast<MathLib::Point3d const*>(element.getNode(0)))));
129  if (c.first) {
130  min = c.second;
131  max = min;
132  } else {
133  return false;
134  }
135 
136  std::vector<std::array<std::size_t,3>> coord_vecs(element.getNumberOfNodes());
137  for (std::size_t k(1); k<element.getNumberOfNodes(); ++k) {
138  // compute coordinates of the grid for each node of the element
139  c = getGridCellCoordinates(*(static_cast<MathLib::Point3d const*>(element.getNode(k))));
140  if (!c.first)
141  return false;
142 
143  for (std::size_t j(0); j < 3; ++j)
144  {
145  if (min[j] > c.second[j])
146  min[j] = c.second[j];
147  if (max[j] < c.second[j])
148  max[j] = c.second[j];
149  }
150  }
151 
152  const std::size_t n_plane(_n_steps[0]*_n_steps[1]);
153 
154  // If a node of an element is almost equal to the upper right point of the
155  // AABB the grid cell coordinates computed by getGridCellCoordintes() could
156  // be to large (due to numerical errors). The following lines ensure that
157  // the grid cell coordinates are in the valid range.
158  for (std::size_t k(0); k<3; ++k)
159  max[k] = std::min(_n_steps[k]-1, max[k]);
160 
161  // insert the element into the grid cells
162  for (std::size_t i(min[0]); i<=max[0]; i++) {
163  for (std::size_t j(min[1]); j<=max[1]; j++) {
164  for (std::size_t k(min[2]); k<=max[2]; k++) {
165  _elements_in_grid_box[i+j*_n_steps[0]+k*n_plane].push_back(&element);
166  }
167  }
168  }
169 
170  return true;
171 }
172 
173 std::pair<bool, std::array<std::size_t, 3>>
175 {
176  bool valid(true);
177  std::array<std::size_t, 3> coords;
178 
179  for (std::size_t k(0); k<3; ++k) {
180  const double d(p[k]-_aabb.getMinPoint()[k]);
181  if (d < 0.0) {
182  valid = false;
183  coords[k] = 0;
184  } else if (_aabb.getMaxPoint()[k] <= p[k]) {
185  valid = false;
186  coords[k] = _n_steps[k]-1;
187  } else {
188  coords[k] = static_cast<std::size_t>(d * _inverse_step_sizes[k]);
189  }
190  }
191 
192  return std::make_pair(valid, coords);
193 }
194 
195 #ifndef NDEBUG
197  GeoLib::GEOObjects& geometries,
198  std::string& geometry_name)
199 {
200  std::vector<std::string> cell_names;
201 
202  auto addPoints = [&geometries](
203  MathLib::Point3d const& p, std::array<double, 3> const& d,
204  std::array<std::size_t, 3> const& c, std::string& name) {
205  auto pnts = std::make_unique<std::vector<GeoLib::Point*>>();
206  pnts->push_back(new GeoLib::Point(p[0]+c[0]*d[0], p[1]+c[1]*d[1], p[2]+c[2]*d[2]));
207  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]));
208  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]));
209  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]));
210  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]));
211  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]));
212  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]));
213  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]));
214  std::array<double,3> ulps; // unit in the last place
215  double const towards(std::numeric_limits<double>::max());
216  ulps[0] = std::nextafter(d[0], towards)-d[0];
217  ulps[1] = std::nextafter(d[1], towards)-d[1];
218  ulps[2] = std::nextafter(d[2], towards)-d[2];
219  double const tolerance(std::min(std::min(ulps[0], ulps[1]), ulps[2]));
220  geometries.addPointVec(std::move(pnts), name, nullptr, tolerance);
221  };
222 
223  for (std::size_t i(0); i<grid._n_steps[0]; ++i) {
224  for (std::size_t j(0); j<grid._n_steps[1]; ++j) {
225  for (std::size_t k(0); k<grid._n_steps[2]; ++k) {
226  cell_names.emplace_back("Grid-"+std::to_string(i)+"-"
227  +std::to_string(j)+"-"+std::to_string(k));
228  addPoints(grid._aabb.getMinPoint(), grid._step_sizes,
229  {{i, j, k}}, cell_names.back());
230  auto plys = std::make_unique<std::vector<GeoLib::Polyline*>>();
231  auto & points = *geometries.getPointVec(cell_names.back());
232 
233  auto* ply_bottom(new GeoLib::Polyline(points));
234  for (std::size_t l(0); l < 4; ++l)
235  ply_bottom->addPoint(l);
236  ply_bottom->addPoint(0); // close to bottom surface
237  plys->push_back(ply_bottom);
238 
239  auto* ply_top(new GeoLib::Polyline(points));
240  for (std::size_t l(4); l<8; ++l)
241  ply_top->addPoint(l);
242  ply_top->addPoint(4); // close to top surface
243  plys->push_back(ply_top);
244 
245  auto* ply_04(new GeoLib::Polyline(points));
246  ply_04->addPoint(0);
247  ply_04->addPoint(4);
248  plys->push_back(ply_04);
249 
250  auto* ply_15(new GeoLib::Polyline(points));
251  ply_15->addPoint(1);
252  ply_15->addPoint(5);
253  plys->push_back(ply_15);
254 
255  auto* ply_26(new GeoLib::Polyline(points));
256  ply_26->addPoint(2);
257  ply_26->addPoint(6);
258  plys->push_back(ply_26);
259 
260  auto* ply_37(new GeoLib::Polyline(points));
261  ply_37->addPoint(3);
262  ply_37->addPoint(7);
263  plys->push_back(ply_37);
264 
265  geometries.addPolylineVec(std::move(plys), cell_names.back(),
266  nullptr);
267  }
268  }
269  }
270  if (geometries.mergeGeometries(cell_names, geometry_name) == 2)
271  geometry_name = cell_names.front();
272 }
273 #endif // NDEBUG
274 
275 } // end namespace MeshLib
MathLib::Point3d const & getMaxPoint() const
Definition: AABB.h:161
Container class for geometric objects.
Definition: GEOObjects.h:62
void addPointVec(std::unique_ptr< std::vector< Point *>> points, std::string &name, std::unique_ptr< std::map< std::string, std::size_t >> pnt_names=nullptr, double eps=sqrt(std::numeric_limits< double >::epsilon()))
Definition: GEOObjects.cpp:39
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
MathLib::Point3d const & getMinPoint() const
Definition: AABB.h:154
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:124
MathLib::Point3d const & getMaxPoint() const
void sortElementsInGridCells(MeshLib::Mesh const &sfc_mesh)
Interface for heuristic search length strategy.
Definition: ProjectData.h:28
const std::vector< Point * > * getPointVec(const std::string &name) const
Definition: GEOObjects.cpp:57
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:71
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:375
MeshElementGrid(MeshLib::Mesh const &mesh)
const Node * getNode(unsigned i) const
Definition: Element.cpp:144