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