31 for (std::size_t k(0); k < 3; ++k)
34 std::nexttoward(max[k], std::numeric_limits<double>::max()) -
36 if (std::abs(max[k] - min[k]) >
tolerance)
44 auto const& min_pnt(_aabb.getMinPoint());
45 auto const& max_pnt(_aabb.getMaxPoint());
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]}};
52 const std::size_t n_eles(mesh.getNumberOfElements());
53 const std::size_t n_eles_per_cell(100);
62 auto sc_ceil = [](
double v)
63 {
return static_cast<std::size_t
>(std::ceil(v)); };
69 sc_ceil(std::cbrt(n_eles * delta[0] * delta[0] /
70 (n_eles_per_cell * delta[1] * delta[2])));
72 sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0));
74 sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0));
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]);
83 else if (dim[0] && dim[1])
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]);
89 else if (dim[1] && dim[2])
91 _n_steps[1] = sc_ceil(std::sqrt(n_eles * delta[1] /
92 (n_eles_per_cell * delta[2])));
94 sc_ceil(n_eles * delta[2] / (n_eles_per_cell * delta[1]));
98 for (std::size_t k(0); k < 3; ++k)
103 sc_ceil(
static_cast<double>(n_eles) / n_eles_per_cell);
110 for (std::size_t k(0); k < 3; k++)
112 _step_sizes[k] = delta[k] / _n_steps[k];
113 _inverse_step_sizes[k] = 1.0 / _step_sizes[k];
116 _elements_in_grid_box.resize(_n_steps[0] * _n_steps[1] * _n_steps[2]);
117 sortElementsInGridCells(mesh);
136 OGS_FATAL(
"Sorting element (id={:d}) into mesh element grid.",
144 std::array<std::size_t, 3> min{};
145 std::array<std::size_t, 3> max{};
168 for (std::size_t j(0); j < 3; ++j)
170 if (min[j] >
c.second[j])
172 min[j] =
c.second[j];
174 if (max[j] <
c.second[j])
176 max[j] =
c.second[j];
187 for (std::size_t k(0); k < 3; ++k)
189 max[k] = std::min(
_n_steps[k] - 1, max[k]);
193 for (std::size_t i(min[0]); i <= max[0]; i++)
195 for (std::size_t j(min[1]); j <= max[1]; j++)
197 for (std::size_t k(min[2]); k <= max[2]; k++)
200 .push_back(&element);
208 std::pair<bool, std::array<std::size_t, 3>>
212 std::array<std::size_t, 3> coords{};
214 for (std::size_t k(0); k < 3; ++k)
233 return std::make_pair(valid, coords);
Definition of the Element class.
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.
Eigen::Vector3d const & getMinPoint() const
Eigen::Vector3d const & getMaxPoint() const
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.
static double const tolerance