Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
/*
* \brief Declaration of the MeshElementGrid class.
*
* \copyright
* Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/
#ifndef MESHELEMENTGRID_H_
#define MESHELEMENTGRID_H_
#include <array>
#include <limits>
#include <vector>
#include "GeoLib/AABB.h"
#include "MathLib/Point3d.h"
namespace GeoLib {
class GEOObjects;
}
namespace MeshLib {
// forward declarations
class Mesh;
class Element;
/// MeshElementGrid implements a grid data structure supporting search
/// operations and covers a given mesh. It consists of grid cells that are all
/// of the same size. Grid cells contain pointers to intersecting mesh elements.
/// @attention The user has to ensure the validity of the pointers while the
/// MeshElementGrid instance lives.
class MeshElementGrid final {
#ifndef NDEBUG
friend void getGridGeometry(MeshLib::MeshElementGrid const& grid,
GeoLib::GEOObjects& geometries,
std::string& geometry_name);
#endif
public:
/// Constructs a grid. Grid cells contains intersecting mesh elements.
/// @param mesh the MeshLib::Mesh instance the grid will be constructed from
explicit MeshElementGrid(MeshLib::Mesh const& mesh);
/// Fill and return a vector containing elements of all grid cells that have
/// a non-empty intersection with the box that is defined by min and max.
/// @param min min point of the box
/// @param max max point of the box
/// @return a (possible empty) vector of elements
template <typename POINT>
std::vector<MeshLib::Element const*> getElementsInVolume(
POINT const& min, POINT const& max) const
{
auto const min_coords(getGridCellCoordinates(min));
auto const max_coords(getGridCellCoordinates(max));
if (!min_coords.first) {
WARN(
"MeshElementGrid::getElementsInVolume: Min point (%f,%f,%f) "
"outside of MeshElementGrid [%f,%f) x [%f,%f) x [%f,%f).",
min[0], min[1], min[2], _aabb.getMinPoint()[0],
_aabb.getMaxPoint()[0], _aabb.getMinPoint()[1],
_aabb.getMaxPoint()[1], _aabb.getMinPoint()[2],
_aabb.getMaxPoint()[2]);
}
if (!max_coords.first) {
WARN(
"MeshElementGrid::getElementsInVolume: Max point (%f,%f,%f) "
"outside of MeshElementGrid [%f,%f) x [%f,%f) x [%f,%f).",
max[0], max[1], max[2], _aabb.getMinPoint()[0],
_aabb.getMaxPoint()[0], _aabb.getMinPoint()[1],
_aabb.getMaxPoint()[1], _aabb.getMinPoint()[2],
_aabb.getMaxPoint()[2]);
}
std::vector<MeshLib::Element const*> elements_vec;
const std::size_t n_plane(_n_steps[0]*_n_steps[1]);
for (std::size_t i(min_coords.second[0]); i<=max_coords.second[0]; i++) {
for (std::size_t j(min_coords.second[1]); j<=max_coords.second[1]; j++) {
for (std::size_t k(min_coords.second[2]); k<=max_coords.second[2]; k++) {
std::size_t idx(i+j*_n_steps[0]+k*n_plane);
std::copy(_elements_in_grid_box[idx].begin(),
_elements_in_grid_box[idx].end(),
elements_vec.end());
}
}
}
return elements_vec;
}
private:
void sortElementsInGridCells(MeshLib::Mesh const& sfc_mesh);
bool sortElementInGridCells(MeshLib::Element const& element);
GeoLib::AABB _aabb;
/// Computes the grid cell coordinates for given point. The first element of
/// the returned pair (bool) is true if the point is within the grid, else
/// false.
std::pair<bool, std::array<std::size_t,3>>
getGridCellCoordinates(MathLib::Point3d const& p) const;
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<MeshLib::Element const*>> _elements_in_grid_box;
};
#ifndef NDEBUG
/// Transfers the grid cells to a geometry. The grid-geometry can be viewed
/// using the DataExplorer. The visualization can be used to explain the
/// algorithm or to find bugs in the algorithm.
void getGridGeometry(MeshLib::MeshElementGrid const& grid,
GeoLib::GEOObjects& geometries,
std::string& geometry_name);
#endif
} // end namespace MeshLib
#endif /* MESHELEMENTGRID_H_ */