Loading [MathJax]/extensions/tex2jax.js
OGS
VoxelGridFromLayeredMeshes.cpp
Go to the documentation of this file.
1
11
12#include <range/v3/algorithm/fill.hpp>
13
14#include "MeshLib/Mesh.h"
19
20static std::string mat_name = "MaterialIDs";
21
22// returns the AABB of all mesh nodes of layers read so far
23void adjustExtent(std::pair<MathLib::Point3d, MathLib::Point3d>& extent,
24 MeshLib::Mesh const& mesh)
25{
26 auto const& nodes = mesh.getNodes();
27 GeoLib::AABB aabb(nodes.cbegin(), nodes.cend());
28 for (std::size_t i = 0; i < 3; ++i)
29 {
30 extent.first[i] = std::min(extent.first[i], aabb.getMinPoint()[i]);
31 extent.second[i] = std::max(extent.second[i], aabb.getMaxPoint()[i]);
32 }
33}
34
35// creates a voxel grid of the AABB of all layers
36std::unique_ptr<MeshLib::Mesh> generateInitialMesh(
37 std::pair<MathLib::Point3d, MathLib::Point3d>& extent,
38 std::array<double, 3> const& res)
39{
40 INFO("Creating initial mesh...");
41 std::array<double, 3> mesh_range{{extent.second[0] - extent.first[0],
42 extent.second[1] - extent.first[1],
43 extent.second[2] - extent.first[2]}};
44 std::array<std::size_t, 3> const n_cells{
45 {static_cast<std::size_t>(std::ceil(mesh_range[0] / res[0])),
46 static_cast<std::size_t>(std::ceil(mesh_range[1] / res[1])),
47 static_cast<std::size_t>(std::ceil(mesh_range[2] / res[2]))}};
48 for (std::size_t i = 0; i < 3; ++i)
49 {
50 double const ext_range = n_cells[i] * res[i];
51 double const offset = (ext_range - mesh_range[i]) / 2.0;
52 mesh_range[i] = ext_range;
53 extent.first[i] -= offset;
54 extent.second[i] += offset;
55 }
56 std::unique_ptr<MeshLib::Mesh> mesh(
58 mesh_range[0], mesh_range[1], mesh_range[2], n_cells[0], n_cells[1],
59 n_cells[2], extent.first));
60 auto mat_id = mesh->getProperties().createNewPropertyVector<int>(
61 mat_name, MeshLib::MeshItemType::Cell, mesh->getNumberOfElements(), 1);
62 if (!mat_id)
63 {
64 return nullptr;
65 }
66 ranges::fill(*mat_id, -1);
67 return mesh;
68}
69
70// returns the element the given node is projected on (or nullptr otherwise)
72 MeshLib::MeshElementGrid const& grid,
73 MathLib::Point3d const& node,
74 double const max_edge)
75{
76 constexpr double max_val = std::numeric_limits<double>::max();
77 MathLib::Point3d const min_vol{
78 {node[0] - max_edge, node[1] - max_edge, -max_val}};
79 MathLib::Point3d const max_vol{
80 {node[0] + max_edge, node[1] + max_edge, max_val}};
81 auto const& intersection_candidates =
82 grid.getElementsInVolume(min_vol, max_vol);
84 intersection_candidates, node);
85}
86
87// casts vote if the given nodes belongs to lower layer, upper layer or no layer
88// at all
89void voteMatId(MathLib::Point3d const& node,
90 MeshLib::MeshElementGrid const& grid,
91 double const max_edge,
92 std::size_t& nullptr_cnt,
93 std::size_t& upper_layer_cnt,
94 std::size_t& lower_layer_cnt)
95{
96 auto const& proj_elem = getProjectedElement(grid, node, max_edge);
97 if (proj_elem == nullptr)
98 {
99 nullptr_cnt++;
100 return;
101 }
102 if (node[2] >
104 {
105 upper_layer_cnt++;
106 return;
107 }
108 lower_layer_cnt++;
109}
110
111// sets material IDs for all elements depending on the layers they are located
112// between
114 std::vector<MeshLib::Mesh const*> const& layers,
115 bool const dilate)
116{
117 INFO("Setting material properties...");
118 std::size_t const n_layers = layers.size();
119 auto const& elems = mesh.getElements();
120 std::size_t const n_elems = mesh.getNumberOfElements();
121 auto mat_ids = mesh.getProperties().getPropertyVector<int>(mat_name);
122 std::vector<bool> is_set(n_elems, false);
123 for (int i = n_layers - 1; i >= 0; --i)
124 {
125 INFO("-> Layer {:d}", n_layers - i - 1);
126 MeshLib::MeshElementGrid const grid(*layers[i]);
127 auto const edgeLengths = minMaxEdgeLength(layers[i]->getElements());
128 double const max_edge = edgeLengths.second;
129 for (std::size_t j = 0; j < n_elems; ++j)
130 {
131 if (is_set[j])
132 {
133 continue;
134 }
135
136 std::size_t nullptr_cnt(0);
137 std::size_t upper_layer_cnt(0);
138 std::size_t lower_layer_cnt(0);
139
140 auto const& node = MeshLib::getCenterOfGravity(*elems[j]);
141 voteMatId(node, grid, max_edge, nullptr_cnt, upper_layer_cnt,
142 lower_layer_cnt);
143 if (nullptr_cnt)
144 {
145 // if no element was found at centre point, vote via corners
146 for (std::size_t k = 0; k < 8; ++k)
147 {
148 MeshLib::Node const& n = *elems[j]->getNode(k);
149 voteMatId(n, grid, max_edge, nullptr_cnt, upper_layer_cnt,
150 lower_layer_cnt);
151 }
152
153 // If the "dilate"-param is set, a mat ID will be assigned if at
154 // least one node was voting for a specific layer. Without the
155 // "dilate"-param, an absolute majority is needed. In case of a
156 // tie, the lower layer will be favoured.
157 if ((upper_layer_cnt == 0 && lower_layer_cnt == 0) ||
158 (!dilate && nullptr_cnt >= upper_layer_cnt &&
159 nullptr_cnt >= lower_layer_cnt))
160 {
161 continue;
162 }
163 if (upper_layer_cnt > lower_layer_cnt)
164 {
165 (*mat_ids)[j] = n_layers - i - 1;
166 }
167 else
168 {
169 is_set[j] = true;
170 }
171 continue;
172 }
173 if (upper_layer_cnt)
174 {
175 (*mat_ids)[j] = n_layers - i - 1;
176 }
177 else
178 {
179 is_set[j] = true;
180 }
181 }
182 }
183 // set all elements above uppermost layer back to -1 so they are
184 // subsequently cut
185 std::replace(mat_ids->begin(), mat_ids->end(),
186 static_cast<int>(n_layers - 1), -1);
187}
188
189// Removes all elements from mesh that have not been marked as being located
190// between two layers. If all elements remain unmarked, a nullptr is returned.
191std::vector<std::size_t> markSpecificElements(MeshLib::Mesh const& mesh,
192 int const mat_id)
193{
194 std::vector<std::size_t> marked_elems;
195 auto const mat_ids = *MeshLib::materialIDs(mesh);
196 std::size_t const n_elems = mat_ids.size();
197 for (std::size_t i = 0; i < n_elems; ++i)
198 {
199 if (mat_ids[i] == mat_id)
200 {
201 marked_elems.push_back(i);
202 }
203 }
204 return marked_elems;
205}
206
207// Creates a VoxelGrid after extending the AABB for each layer.
208std::unique_ptr<MeshLib::Mesh> MeshToolsLib::MeshGenerators::
210 std::pair<MathLib::Point3d, MathLib::Point3d>& extent,
211 std::vector<MeshLib::Mesh const*> const& layers,
212 std::array<double, 3> const cellsize,
213 bool const dilate)
214{
215 for (auto const& layer : layers)
216 {
217 adjustExtent(extent, *layer);
218 }
219
220 std::unique_ptr<MeshLib::Mesh> mesh(generateInitialMesh(extent, cellsize));
221 if (mesh == nullptr)
222 {
223 return nullptr;
224 }
225 setMaterialIDs(*mesh, layers, dilate);
226 auto const marked_elements = markSpecificElements(*mesh, -1);
227 if (marked_elements.size() == mesh->getNumberOfElements())
228 {
229 return nullptr;
230 }
231 std::unique_ptr<MeshLib::Mesh> new_mesh(
232 MeshToolsLib::removeElements(*mesh, marked_elements, "mesh"));
233 return new_mesh;
234}
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
Definition of the Mesh class.
static std::string mat_name
std::vector< std::size_t > markSpecificElements(MeshLib::Mesh const &mesh, int const mat_id)
void voteMatId(MathLib::Point3d const &node, MeshLib::MeshElementGrid const &grid, double const max_edge, std::size_t &nullptr_cnt, std::size_t &upper_layer_cnt, std::size_t &lower_layer_cnt)
MeshLib::Element const * getProjectedElement(MeshLib::MeshElementGrid const &grid, MathLib::Point3d const &node, double const max_edge)
void setMaterialIDs(MeshLib::Mesh &mesh, std::vector< MeshLib::Mesh const * > const &layers, bool const dilate)
std::unique_ptr< MeshLib::Mesh > generateInitialMesh(std::pair< MathLib::Point3d, MathLib::Point3d > &extent, std::array< double, 3 > const &res)
void adjustExtent(std::pair< MathLib::Point3d, MathLib::Point3d > &extent, MeshLib::Mesh const &mesh)
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Definition AABB.h:56
Eigen::Vector3d const & getMaxPoint() const
Definition AABB.h:187
Eigen::Vector3d const & getMinPoint() const
Definition AABB.h:180
std::vector< MeshLib::Element const * > getElementsInVolume(POINT const &min, POINT const &max) const
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:108
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:111
Properties & getProperties()
Definition Mesh.h:136
std::size_t getNumberOfElements() const
Get the number of elements.
Definition Mesh.h:99
PropertyVector< T > const * getPropertyVector(std::string_view name) const
PropertyVector< int > const * materialIDs(Mesh const &mesh)
Definition Mesh.cpp:269
MathLib::Point3d getCenterOfGravity(Element const &element)
Calculates the center of gravity for the mesh element.
Definition Element.cpp:142
MeshLib::Mesh * generateRegularHexMesh(const BaseLib::ISubdivision &div_x, const BaseLib::ISubdivision &div_y, const BaseLib::ISubdivision &div_z, MathLib::Point3d const &origin=MathLib::ORIGIN, std::string const &mesh_name="mesh")
std::unique_ptr< MeshLib::Mesh > createVoxelFromLayeredMesh(std::pair< MathLib::Point3d, MathLib::Point3d > &extent, std::vector< MeshLib::Mesh const * > const &layers, std::array< double, 3 > const cellsize, bool const dilate)
double getElevation(MeshLib::Element const &element, MathLib::Point3d const &node)
MeshLib::Element const * getProjectedElement(std::vector< const MeshLib::Element * > const &elements, MathLib::Point3d const &node)
MeshLib::Mesh * removeElements(const MeshLib::Mesh &mesh, const std::vector< std::size_t > &removed_element_ids, const std::string &new_mesh_name)