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