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>
36 MeshLib::Mesh const& mesh, std::vector<float>
const& layer_thickness_vector,
37 std::string
const& mesh_name)
39 std::vector<float> thickness;
40 for (std::size_t i = 0; i < layer_thickness_vector.size(); ++i)
42 if (layer_thickness_vector[i] > std::numeric_limits<float>::epsilon())
44 thickness.push_back(layer_thickness_vector[i]);
48 WARN(
"Ignoring layer {:d} with thickness {:f}.", i,
49 layer_thickness_vector[i]);
53 const std::size_t nLayers(thickness.size());
56 ERR(
"MeshLayerMapper::createStaticLayers(): A 2D mesh with nLayers > 0 "
57 "is required as input.");
63 const std::size_t nElems(
66 { return (elem->getDimension() == 2); }));
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);
79 ERR(
"Could not create PropertyVector object 'MaterialIDs'.");
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) |
91 for (
unsigned layer_id = 0; layer_id <= nLayers; ++layer_id)
94 unsigned node_offset(nNodes * layer_id);
97 z_offset += thickness[layer_id - 1];
100 std::transform(nodes.cbegin(), nodes.cend(),
101 new_nodes.begin() + node_offset,
103 return new MeshLib::Node((*node)[0], (*node)[1],
104 (*node)[2] - z_offset);
114 node_offset -= nNodes;
116 for (
unsigned i = 0; i < nOrgElems; ++i)
127 for (
unsigned j = 0; j < nElemNodes; ++j)
129 const unsigned node_id =
131 e_nodes[j] = new_nodes[node_id + nNodes];
132 e_nodes[j + nElemNodes] = new_nodes[node_id];
146 OGS_FATAL(
"MeshLayerMapper: Unknown element type to extrude.");
156 std::vector<GeoLib::Raster const*>
const& rasters,
157 double minimum_thickness,
158 double noDataReplacementValue)
160 const std::size_t nLayers(rasters.size());
163 ERR(
"MeshLayerMapper::createRasterLayers(): A 2D mesh and at least two "
164 "rasters required as input.");
168 auto top = std::make_unique<MeshLib::Mesh>(mesh);
169 if (!
layerMapping(*top, *rasters.back(), noDataReplacementValue))
174 auto bottom = std::make_unique<MeshLib::Mesh>(mesh);
182 _nodes.reserve(nLayers * nNodes);
185 std::size_t
const nElems(std::count_if(
188 { return (elem->getGeomType() == MeshLib::MeshElemType::TRIANGLE); }));
189 _elements.reserve(nElems * (nLayers - 1));
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); });
198 for (std::size_t i = 0; i < nLayers - 1; ++i)
210 const unsigned pyramid_base[3][4] = {
217 std::vector<MeshLib::Node*>
const& top_nodes = dem_mesh.
getNodes();
218 int const last_layer_node_offset = layer_id * nNodes;
221 for (std::size_t i = 0; i < nNodes; ++i)
224 *
_nodes[last_layer_node_offset + i],
228 std::vector<MeshLib::Element*>
const& elems = dem_mesh.
getElements();
231 for (std::size_t i = 0; i < nElems; ++i)
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)
244 _nodes[
_nodes[last_layer_node_offset + getNodeIndex(*elem, j)]
246 new_elem_nodes[node_counter] =
247 (
_nodes[last_layer_node_offset + getNodeIndex(*elem, j) +
249 if (new_elem_nodes[j]->getID() !=
250 new_elem_nodes[node_counter]->getID())
260 switch (node_counter)
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];
280 std::array<MeshLib::Node*, 4> tet_nodes{};
281 std::copy(new_elem_nodes.begin(),
282 new_elem_nodes.begin() + node_counter,
296 double const nodata_replacement,
297 bool const ignore_nodata)
301 ERR(
"MeshLayerMapper::layerMapping() - requires 2D mesh");
306 const double x0(header.
origin[0]);
307 const double y0(header.
origin[1]);
310 const std::pair<double, double> xDim(
311 x0, x0 + header.
n_cols * delta);
312 const std::pair<double, double> yDim(
313 y0, y0 + header.
n_rows * delta);
316 const std::vector<MeshLib::Node*>& nodes = mesh.
getNodes();
317 for (
unsigned i = 0; i < nNodes; ++i)
322 (*nodes[i])[2] = nodata_replacement;
327 constexpr double eps = std::numeric_limits<double>::epsilon();
328 if (std::abs(elevation - header.
no_data) < eps)
334 elevation = nodata_replacement;
340 (*nodes[i])[2] = elevation;
351 ERR(
"MshLayerMapper::mapToStaticValue() - requires 2D mesh");
355 std::vector<MeshLib::Node*>
const& nodes(mesh.
getNodes());
Definition of the Hex class.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition of the MeshLayerMapper class.
Definition of the class Properties that implements a container of properties.
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.
double getValueAtPoint(const MathLib::Point3d &pnt) const
bool isPntOnRaster(MathLib::Point3d const &pnt) const
double interpolateValueAtPoint(const MathLib::Point3d &pnt) const
RasterHeader const & getHeader() const
Returns the complete header information.
double _minimum_thickness
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.
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
std::size_t getNumberOfNodes() const
Get the number of nodes.
std::size_t getNumberOfElements() const
Get the number of elements.
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
PropertyVector< T > * createNewPropertyVector(std::string_view name, MeshItemType mesh_item_type, std::size_t n_components=1)