OGS
MeshLayerMapper.cpp
Go to the documentation of this file.
1
15#include "MeshLayerMapper.h"
16
17#include <algorithm>
18#include <range/v3/view/iota.hpp>
19#include <range/v3/view/join.hpp>
20#include <range/v3/view/repeat_n.hpp>
21#include <range/v3/view/transform.hpp>
22
23#include "BaseLib/Logging.h"
24#include "GeoLib/Raster.h"
25#include "MathLib/MathTools.h"
30#include "MeshLib/Properties.h"
32
33namespace MeshToolsLib
34{
36 MeshLib::Mesh const& mesh, std::vector<float> const& layer_thickness_vector,
37 std::string const& mesh_name)
38{
39 std::vector<float> thickness;
40 for (std::size_t i = 0; i < layer_thickness_vector.size(); ++i)
41 {
42 if (layer_thickness_vector[i] > std::numeric_limits<float>::epsilon())
43 {
44 thickness.push_back(layer_thickness_vector[i]);
45 }
46 else
47 {
48 WARN("Ignoring layer {:d} with thickness {:f}.", i,
49 layer_thickness_vector[i]);
50 }
51 }
52
53 const std::size_t nLayers(thickness.size());
54 if (nLayers < 1 || mesh.getDimension() != 2)
55 {
56 ERR("MeshLayerMapper::createStaticLayers(): A 2D mesh with nLayers > 0 "
57 "is required as input.");
58 return nullptr;
59 }
60
61 const std::size_t nNodes = mesh.getNumberOfNodes();
62 // count number of 2d elements in the original mesh
63 const std::size_t nElems(
64 std::count_if(mesh.getElements().begin(), mesh.getElements().end(),
65 [](MeshLib::Element const* elem)
66 { return (elem->getDimension() == 2); }));
67
68 const std::size_t nOrgElems(mesh.getNumberOfElements());
69 const std::vector<MeshLib::Node*>& nodes = mesh.getNodes();
70 const std::vector<MeshLib::Element*>& elems = mesh.getElements();
71 std::vector<MeshLib::Node*> new_nodes(nNodes + (nLayers * nNodes));
72 std::vector<MeshLib::Element*> new_elems;
73 new_elems.reserve(nElems * nLayers);
74 MeshLib::Properties properties;
75 auto* const materials = properties.createNewPropertyVector<int>(
76 "MaterialIDs", MeshLib::MeshItemType::Cell);
77 if (!materials)
78 {
79 ERR("Could not create PropertyVector object 'MaterialIDs'.");
80 return nullptr;
81 }
82
83 auto layer_to_material_id = [&](unsigned const layer_id)
84 { return ranges::views::repeat_n(nLayers - layer_id, nOrgElems); };
85 materials->assign(ranges::views::iota(0u, nLayers + 1) |
86 ranges::views::transform(layer_to_material_id) |
87 ranges::views::join);
88
89 double z_offset(0.0);
90
91 for (unsigned layer_id = 0; layer_id <= nLayers; ++layer_id)
92 {
93 // add nodes for new layer
94 unsigned node_offset(nNodes * layer_id);
95 if (layer_id > 0)
96 {
97 z_offset += thickness[layer_id - 1];
98 }
99
100 std::transform(nodes.cbegin(), nodes.cend(),
101 new_nodes.begin() + node_offset,
102 [&z_offset](MeshLib::Node* node) {
103 return new MeshLib::Node((*node)[0], (*node)[1],
104 (*node)[2] - z_offset);
105 });
106
107 // starting with 2nd layer create prism or hex elements connecting the
108 // last layer with the current one
109 if (layer_id == 0)
110 {
111 continue;
112 }
113
114 node_offset -= nNodes;
115
116 for (unsigned i = 0; i < nOrgElems; ++i)
117 {
118 const MeshLib::Element* sfc_elem(elems[i]);
119 if (sfc_elem->getDimension() < 2)
120 { // ignore line-elements
121 continue;
122 }
123
124 const unsigned nElemNodes(sfc_elem->getNumberOfBaseNodes());
125 auto** e_nodes = new MeshLib::Node*[2 * nElemNodes];
126
127 for (unsigned j = 0; j < nElemNodes; ++j)
128 {
129 const unsigned node_id =
130 sfc_elem->getNode(j)->getID() + node_offset;
131 e_nodes[j] = new_nodes[node_id + nNodes];
132 e_nodes[j + nElemNodes] = new_nodes[node_id];
133 }
135 {
136 // extrude triangles to prism
137 new_elems.push_back(new MeshLib::Prism(e_nodes));
138 }
139 else if (sfc_elem->getGeomType() == MeshLib::MeshElemType::QUAD)
140 {
141 // extrude quads to hexes
142 new_elems.push_back(new MeshLib::Hex(e_nodes));
143 }
144 else
145 {
146 OGS_FATAL("MeshLayerMapper: Unknown element type to extrude.");
147 }
148 }
149 }
150 return new MeshLib::Mesh(mesh_name, new_nodes, new_elems,
151 true /* compute_element_neighbors */, properties);
152}
153
155 MeshLib::Mesh const& mesh,
156 std::vector<GeoLib::Raster const*> const& rasters,
157 double minimum_thickness,
158 double noDataReplacementValue)
159{
160 const std::size_t nLayers(rasters.size());
161 if (nLayers < 2 || mesh.getDimension() != 2)
162 {
163 ERR("MeshLayerMapper::createRasterLayers(): A 2D mesh and at least two "
164 "rasters required as input.");
165 return false;
166 }
167
168 auto top = std::make_unique<MeshLib::Mesh>(mesh);
169 if (!layerMapping(*top, *rasters.back(), noDataReplacementValue))
170 {
171 return false;
172 }
173
174 auto bottom = std::make_unique<MeshLib::Mesh>(mesh);
175 if (!layerMapping(*bottom, *rasters[0], 0))
176 {
177 return false;
178 }
179
180 this->_minimum_thickness = minimum_thickness;
181 std::size_t const nNodes = mesh.getNumberOfNodes();
182 _nodes.reserve(nLayers * nNodes);
183
184 // number of triangles in the original mesh
185 std::size_t const nElems(std::count_if(
186 mesh.getElements().begin(), mesh.getElements().end(),
187 [](MeshLib::Element const* elem)
188 { return (elem->getGeomType() == MeshLib::MeshElemType::TRIANGLE); }));
189 _elements.reserve(nElems * (nLayers - 1));
190 _materials.reserve(nElems * (nLayers - 1));
191
192 // add bottom layer
193 std::vector<MeshLib::Node*> const& nodes = bottom->getNodes();
194 std::transform(nodes.begin(), nodes.end(), std::back_inserter(_nodes),
195 [](auto const* node) { return new MeshLib::Node(*node); });
196
197 // add the other layers
198 for (std::size_t i = 0; i < nLayers - 1; ++i)
199 {
200 addLayerToMesh(*top, i, *rasters[i + 1]);
201 }
202
203 return true;
204}
205
207 unsigned layer_id,
208 GeoLib::Raster const& raster)
209{
210 const unsigned pyramid_base[3][4] = {
211 {1, 3, 4, 2}, // Point 4 missing
212 {2, 4, 3, 0}, // Point 5 missing
213 {0, 3, 4, 1}, // Point 6 missing
214 };
215
216 std::size_t const nNodes = dem_mesh.getNumberOfNodes();
217 std::vector<MeshLib::Node*> const& top_nodes = dem_mesh.getNodes();
218 int const last_layer_node_offset = layer_id * nNodes;
219
220 // add nodes for new layer
221 for (std::size_t i = 0; i < nNodes; ++i)
222 {
223 _nodes.push_back(getNewLayerNode(*top_nodes[i],
224 *_nodes[last_layer_node_offset + i],
225 raster, _nodes.size()));
226 }
227
228 std::vector<MeshLib::Element*> const& elems = dem_mesh.getElements();
229 std::size_t const nElems(dem_mesh.getNumberOfElements());
230
231 for (std::size_t i = 0; i < nElems; ++i)
232 {
233 MeshLib::Element* elem(elems[i]);
235 {
236 continue;
237 }
238 unsigned node_counter(3);
239 unsigned missing_idx(0);
240 std::array<MeshLib::Node*, 6> new_elem_nodes{};
241 for (unsigned j = 0; j < 3; ++j)
242 {
243 new_elem_nodes[j] =
244 _nodes[_nodes[last_layer_node_offset + getNodeIndex(*elem, j)]
245 ->getID()];
246 new_elem_nodes[node_counter] =
247 (_nodes[last_layer_node_offset + getNodeIndex(*elem, j) +
248 nNodes]);
249 if (new_elem_nodes[j]->getID() !=
250 new_elem_nodes[node_counter]->getID())
251 {
252 node_counter++;
253 }
254 else
255 {
256 missing_idx = j;
257 }
258 }
259
260 switch (node_counter)
261 {
262 case 6:
263 _elements.push_back(new MeshLib::Prism(new_elem_nodes));
264 _materials.push_back(layer_id);
265 break;
266 case 5:
267 {
268 std::array<MeshLib::Node*, 5> pyramid_nodes{};
269 pyramid_nodes[0] = new_elem_nodes[pyramid_base[missing_idx][0]];
270 pyramid_nodes[1] = new_elem_nodes[pyramid_base[missing_idx][1]];
271 pyramid_nodes[2] = new_elem_nodes[pyramid_base[missing_idx][2]];
272 pyramid_nodes[3] = new_elem_nodes[pyramid_base[missing_idx][3]];
273 pyramid_nodes[4] = new_elem_nodes[missing_idx];
274 _elements.push_back(new MeshLib::Pyramid(pyramid_nodes));
275 _materials.push_back(layer_id);
276 break;
277 }
278 case 4:
279 {
280 std::array<MeshLib::Node*, 4> tet_nodes{};
281 std::copy(new_elem_nodes.begin(),
282 new_elem_nodes.begin() + node_counter,
283 tet_nodes.begin());
284 _elements.push_back(new MeshLib::Tet(tet_nodes));
285 _materials.push_back(layer_id);
286 break;
287 }
288 default:
289 continue;
290 }
291 }
292}
293
295 GeoLib::Raster const& raster,
296 double const nodata_replacement,
297 bool const ignore_nodata)
298{
299 if (mesh.getDimension() != 2)
300 {
301 ERR("MeshLayerMapper::layerMapping() - requires 2D mesh");
302 return false;
303 }
304
305 GeoLib::RasterHeader const& header(raster.getHeader());
306 const double x0(header.origin[0]);
307 const double y0(header.origin[1]);
308 const double delta(header.cell_size);
309
310 const std::pair<double, double> xDim(
311 x0, x0 + header.n_cols * delta); // extension in x-dimension
312 const std::pair<double, double> yDim(
313 y0, y0 + header.n_rows * delta); // extension in y-dimension
314
315 const std::size_t nNodes(mesh.getNumberOfNodes());
316 const std::vector<MeshLib::Node*>& nodes = mesh.getNodes();
317 for (unsigned i = 0; i < nNodes; ++i)
318 {
319 if (!ignore_nodata && !raster.isPntOnRaster(*nodes[i]))
320 {
321 // use either default value or elevation from layer above
322 (*nodes[i])[2] = nodata_replacement;
323 continue;
324 }
325
326 double elevation(raster.getValueAtPoint(*nodes[i]));
327 constexpr double eps = std::numeric_limits<double>::epsilon();
328 if (std::abs(elevation - header.no_data) < eps)
329 {
330 if (ignore_nodata)
331 {
332 continue;
333 }
334 elevation = nodata_replacement;
335 }
336 else
337 {
338 elevation = raster.interpolateValueAtPoint(*nodes[i]);
339 }
340 (*nodes[i])[2] = elevation;
341 }
342
343 return true;
344}
345
347 double const value)
348{
349 if (mesh.getDimension() != 2)
350 {
351 ERR("MshLayerMapper::mapToStaticValue() - requires 2D mesh");
352 return false;
353 }
354
355 std::vector<MeshLib::Node*> const& nodes(mesh.getNodes());
356 for (MeshLib::Node* node : nodes)
357 {
358 (*node)[2] = value;
359 }
360 return true;
361}
362
363} // namespace MeshToolsLib
#define OGS_FATAL(...)
Definition Error.h:26
Definition of the Hex class.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
Definition of the MeshLayerMapper class.
Definition of the class Properties that implements a container of properties.
Definition of the MeshSurfaceExtraction class.
Definition of the Prism class.
Definition of the Pyramid class.
Definition of the GeoLib::Raster class.
Definition of the Tet class.
Class Raster is used for managing raster data.
Definition Raster.h:49
double getValueAtPoint(const MathLib::Point3d &pnt) const
Definition Raster.cpp:56
bool isPntOnRaster(MathLib::Point3d const &pnt) const
Definition Raster.cpp:151
double interpolateValueAtPoint(const MathLib::Point3d &pnt) const
Definition Raster.cpp:85
RasterHeader const & getHeader() const
Returns the complete header information.
Definition Raster.h:85
std::vector< int > _materials
MeshLib::Node * getNewLayerNode(MeshLib::Node const &dem_node, MeshLib::Node const &last_layer_node, GeoLib::Raster const &raster, std::size_t new_node_id) const
std::vector< MeshLib::Node * > _nodes
std::vector< MeshLib::Element * > _elements
std::size_t getID() const
virtual MeshElemType getGeomType() const =0
virtual unsigned getNumberOfBaseNodes() const =0
virtual const Node * getNode(unsigned idx) const =0
virtual constexpr unsigned getDimension() const =0
Get dimension of the mesh element.
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
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition Mesh.h:90
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition Mesh.h:102
std::size_t getNumberOfElements() const
Get the number of elements.
Definition Mesh.h:99
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
Definition Properties.h:33
PropertyVector< T > * createNewPropertyVector(std::string_view name, MeshItemType mesh_item_type, std::size_t n_components=1)
Definition Properties.h:17
bool createRasterLayers(MeshLib::Mesh const &mesh, std::vector< GeoLib::Raster const * > const &rasters, double minimum_thickness, double noDataReplacementValue=0.0) override
static MeshLib::Mesh * createStaticLayers(MeshLib::Mesh const &mesh, std::vector< float > const &layer_thickness_vector, std::string const &mesh_name="SubsurfaceMesh")
void addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer_id, GeoLib::Raster const &raster) override
Adds another layer to a subsurface mesh.
static bool layerMapping(MeshLib::Mesh const &mesh, const GeoLib::Raster &raster, double nodata_replacement=0.0, bool const ignore_nodata=false)
static bool mapToStaticValue(MeshLib::Mesh const &mesh, double value)
Maps the elevation of all mesh nodes to the specified static value.
Contains the relevant information when storing a geoscientific raster data.
Definition Raster.h:28
MathLib::Point3d origin
Definition Raster.h:32
std::size_t n_cols
Definition Raster.h:29
std::size_t n_rows
Definition Raster.h:30