OGS
MeshLib::MeshLayerMapper Class Referencefinal

Detailed Description

Manipulating and adding prism element layers to an existing 2D mesh.

Definition at line 25 of file MeshLayerMapper.h.

#include <MeshLayerMapper.h>

Inheritance diagram for MeshLib::MeshLayerMapper:
[legend]
Collaboration diagram for MeshLib::MeshLayerMapper:
[legend]

Public Member Functions

 ~MeshLayerMapper () override=default
 
bool createRasterLayers (MeshLib::Mesh const &mesh, std::vector< GeoLib::Raster const * > const &rasters, double minimum_thickness, double noDataReplacementValue=0.0) override
 
- Public Member Functions inherited from LayeredMeshGenerator
virtual bool createLayers (MeshLib::Mesh const &mesh, std::vector< GeoLib::Raster const * > const &rasters, double minimum_thickness, double noDataReplacementValue=0.0) final
 
std::unique_ptr< MeshLib::MeshgetMesh (std::string const &mesh_name) const
 Returns a mesh of the subsurface representation. More...
 

Static Public Member Functions

static MeshLib::MeshcreateStaticLayers (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. More...
 

Private Member Functions

void addLayerToMesh (const MeshLib::Mesh &dem_mesh, unsigned layer_id, GeoLib::Raster const &raster) override
 Adds another layer to a subsurface mesh. More...
 

Additional Inherited Members

- Protected Member Functions inherited from LayeredMeshGenerator
 LayeredMeshGenerator ()=default
 
virtual ~LayeredMeshGenerator ()=default
 
MeshLib::NodegetNewLayerNode (MeshLib::Node const &dem_node, MeshLib::Node const &last_layer_node, GeoLib::Raster const &raster, std::size_t new_node_id) const
 
- Static Protected Member Functions inherited from LayeredMeshGenerator
static double calcEpsilon (GeoLib::Raster const &low, GeoLib::Raster const &high)
 Calculates a data-dependent epsilon value. More...
 
- Protected Attributes inherited from LayeredMeshGenerator
double _elevation_epsilon {0.0001}
 
double _minimum_thickness {std::numeric_limits<float>::epsilon()}
 
std::vector< int > _materials
 
std::vector< MeshLib::Node * > _nodes
 
std::vector< MeshLib::Element * > _elements
 

Constructor & Destructor Documentation

◆ ~MeshLayerMapper()

MeshLib::MeshLayerMapper::~MeshLayerMapper ( )
overridedefault

Member Function Documentation

◆ addLayerToMesh()

void MeshLib::MeshLayerMapper::addLayerToMesh ( const MeshLib::Mesh dem_mesh,
unsigned  layer_id,
GeoLib::Raster const &  raster 
)
overrideprivatevirtual

Adds another layer to a subsurface mesh.

Implements LayeredMeshGenerator.

Definition at line 198 of file MeshLayerMapper.cpp.

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]);
226  if (elem->getGeomType() != MeshLib::MeshElemType::TRIANGLE)
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 }
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::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition: Mesh.h:95
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
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

References LayeredMeshGenerator::_elements, LayeredMeshGenerator::_materials, LayeredMeshGenerator::_nodes, MathLib::LinAlg::copy(), MeshLib::Mesh::getElements(), MeshLib::Element::getGeomType(), LayeredMeshGenerator::getNewLayerNode(), MeshLib::getNodeIndex(), MeshLib::Mesh::getNodes(), MeshLib::Mesh::getNumberOfElements(), MeshLib::Mesh::getNumberOfNodes(), and MeshLib::TRIANGLE.

Referenced by createRasterLayers().

◆ createRasterLayers()

bool MeshLib::MeshLayerMapper::createRasterLayers ( MeshLib::Mesh const &  mesh,
std::vector< GeoLib::Raster const * > const &  rasters,
double  minimum_thickness,
double  noDataReplacementValue = 0.0 
)
overridevirtual

Based on a 2D triangle mesh this method creates a 3D mesh with a given number of prism-layers. Note: While this method would technically also work with quad meshes, this is discouraged as quad elements will most likely not be coplanar after the mapping process which result in invaled mesh elements.

Parameters
meshThe 2D triangle mesh that is the basis for the new 3D prism mesh
rastersContaining all the raster-data for the subsurface layers from bottom to top (starting with the bottom of the oldest layer and ending with the DEM)
minimum_thicknessMinimum thickness of each of the newly created layers (i.e. nodes with a vertical distance smaller than this will be collapsed)
noDataReplacementValueDefault z-coordinate if there are mesh nodes not located on the DEM raster (i.e. raster_paths[0])
Returns
A mesh with the requested number of layers of prism elements (also including Tet- & Pyramid-elements in case of degenerated prisms)

Implements LayeredMeshGenerator.

Definition at line 146 of file MeshLayerMapper.cpp.

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 }
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42
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)

References LayeredMeshGenerator::_elements, LayeredMeshGenerator::_materials, LayeredMeshGenerator::_minimum_thickness, LayeredMeshGenerator::_nodes, addLayerToMesh(), ERR(), MeshLib::Mesh::getDimension(), MeshLib::Mesh::getElements(), MeshLib::Mesh::getNumberOfNodes(), and layerMapping().

◆ createStaticLayers()

MeshLib::Mesh * MeshLib::MeshLayerMapper::createStaticLayers ( MeshLib::Mesh const &  mesh,
std::vector< float > const &  layer_thickness_vector,
std::string const &  mesh_name = "SubsurfaceMesh" 
)
static

Based on a 2D triangle-or quad mesh this method creates a 3D mesh with a given number of prism- or hex-layers

Parameters
meshThe triangle/quad mesh that is the basis for the new prism/hex mesh
layer_thickness_vectorThe size of the vector equals the number of layers of prism/hex elements that will be extruded from the triangle/quad elements of the original mesh. The thickness of the \(i\)-th layer corresponds to the \(i\) entry of the vector.
mesh_nameThe name of the newly created mesh
Returns
A mesh with the requested number of layers of prism/hex elements

Definition at line 31 of file MeshLayerMapper.cpp.

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 }
#define OGS_FATAL(...)
Definition: Error.h:26
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:37
std::size_t getID() const
Definition: Point3dWithID.h:62
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)

References MeshLib::Cell, MeshLib::Properties::createNewPropertyVector(), ERR(), MeshLib::Mesh::getDimension(), MeshLib::Element::getDimension(), MeshLib::Mesh::getElements(), MeshLib::Element::getGeomType(), MathLib::Point3dWithID::getID(), MeshLib::Element::getNode(), MeshLib::Mesh::getNodes(), MeshLib::Element::getNumberOfBaseNodes(), MeshLib::Mesh::getNumberOfElements(), MeshLib::Mesh::getNumberOfNodes(), OGS_FATAL, MeshLib::QUAD, MeshLib::TRIANGLE, and WARN().

Referenced by MeshLayerEditDialog::createPrismMesh(), and MeshLayerEditDialog::createTetMesh().

◆ layerMapping()

bool MeshLib::MeshLayerMapper::layerMapping ( MeshLib::Mesh const &  mesh,
const GeoLib::Raster raster,
double  nodata_replacement = 0.0,
bool const  ignore_nodata = false 
)
static

Maps the elevation of nodes of a given 2D mesh according to the raster. At locations where no information is given, node elevation is set to noDataReplacementValue.

Definition at line 286 of file MeshLayerMapper.cpp.

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 }
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
Contains the relevant information when storing a geoscientific raster data.
Definition: Raster.h:25

References GeoLib::RasterHeader::cell_size, ERR(), MeshLib::Mesh::getDimension(), GeoLib::Raster::getHeader(), MeshLib::Mesh::getNodes(), MeshLib::Mesh::getNumberOfNodes(), GeoLib::Raster::getValueAtPoint(), GeoLib::Raster::interpolateValueAtPoint(), GeoLib::Raster::isPntOnRaster(), GeoLib::RasterHeader::n_cols, GeoLib::RasterHeader::n_rows, GeoLib::RasterHeader::no_data, and GeoLib::RasterHeader::origin.

Referenced by LayeredVolume::createRasterLayers(), createRasterLayers(), and MeshView::openMap2dMeshDialog().

◆ mapToStaticValue()

bool MeshLib::MeshLayerMapper::mapToStaticValue ( MeshLib::Mesh const &  mesh,
double  value 
)
static

Maps the elevation of all mesh nodes to the specified static value.

Definition at line 339 of file MeshLayerMapper.cpp.

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 }

References ERR(), MeshLib::Mesh::getDimension(), and MeshLib::Mesh::getNodes().

Referenced by MeshView::openMap2dMeshDialog().


The documentation for this class was generated from the following files: