OGS
LayeredVolume.cpp
Go to the documentation of this file.
1 
15 #include "LayeredVolume.h"
16 
17 #include "GeoLib/Raster.h"
18 #include "MeshLib/Elements/Quad.h"
19 #include "MeshLib/Elements/Tri.h"
24 #include "MeshLib/Properties.h"
25 
27  const MeshLib::Mesh& mesh,
28  const std::vector<GeoLib::Raster const*>& rasters,
29  double minimum_thickness,
30  double noDataReplacementValue)
31 {
32  if (mesh.getDimension() != 2)
33  {
34  return false;
35  }
36 
37  _elevation_epsilon = calcEpsilon(*rasters[0], *rasters.back());
38  if (_elevation_epsilon <= 0)
39  {
40  return false;
41  }
42 
43  // remove line elements, only tri + quad remain
44  MeshLib::ElementSearch ex(mesh);
46  std::unique_ptr<MeshLib::Mesh> top(
47  removeElements(mesh, ex.getSearchedElementIDs(), "MeshLayer"));
48  if (top == nullptr)
49  {
50  top = std::make_unique<MeshLib::Mesh>(mesh);
51  }
52 
54  *top, *rasters.back(), noDataReplacementValue))
55  {
56  return false;
57  }
58 
59  std::unique_ptr<MeshLib::Mesh> bottom(new MeshLib::Mesh(*top));
60  if (!MeshLib::MeshLayerMapper::layerMapping(*bottom, *rasters[0], 0))
61  {
62  return false;
63  }
64 
65  this->_minimum_thickness = minimum_thickness;
66  _nodes = MeshLib::copyNodeVector(bottom->getNodes());
67  _elements = MeshLib::copyElementVector(bottom->getElements(), _nodes);
68  if (!_materials.empty())
69  {
70  ERR("The materials vector is not empty.");
71  return false;
72  }
73  _materials.resize(_elements.size(), 0);
74 
75  // map each layer and attach to subsurface mesh
76  const std::size_t nRasters(rasters.size());
77  for (std::size_t i = 1; i < nRasters; ++i)
78  {
79  this->addLayerToMesh(*top, i, *rasters[i]);
80  }
81 
82  // close boundaries between layers
83  this->addLayerBoundaries(*top, nRasters);
84  this->removeCongruentElements(nRasters, top->getNumberOfElements());
85 
86  return true;
87 }
88 
90  unsigned layer_id,
91  GeoLib::Raster const& raster)
92 {
93  const std::size_t nNodes(dem_mesh.getNumberOfNodes());
94  const std::vector<MeshLib::Node*>& nodes(dem_mesh.getNodes());
95  const std::size_t node_id_offset(_nodes.size());
96  const std::size_t last_layer_node_offset(node_id_offset - nNodes);
97 
98  for (std::size_t i = 0; i < nNodes; ++i)
99  {
100  _nodes.push_back(getNewLayerNode(*nodes[i],
101  *_nodes[last_layer_node_offset + i],
102  raster,
103  _nodes.size()));
104  }
105 
106  const std::vector<MeshLib::Element*>& layer_elements(
107  dem_mesh.getElements());
108  for (MeshLib::Element* elem : layer_elements)
109  {
110  if (elem->getGeomType() == MeshLib::MeshElemType::TRIANGLE)
111  {
112  std::array<MeshLib::Node*, 3> tri_nodes = {
113  {_nodes[node_id_offset + getNodeIndex(*elem, 0)],
114  _nodes[node_id_offset + getNodeIndex(*elem, 1)],
115  _nodes[node_id_offset + getNodeIndex(*elem, 2)]}};
116  _elements.push_back(new MeshLib::Tri(tri_nodes));
117  _materials.push_back(layer_id);
118  }
119  else if (elem->getGeomType() == MeshLib::MeshElemType::QUAD)
120  {
121  std::array<MeshLib::Node*, 4> quad_nodes = {
122  {_nodes[node_id_offset + getNodeIndex(*elem, 0)],
123  _nodes[node_id_offset + getNodeIndex(*elem, 1)],
124  _nodes[node_id_offset + getNodeIndex(*elem, 2)],
125  _nodes[node_id_offset + getNodeIndex(*elem, 3)]}};
126  _elements.push_back(new MeshLib::Quad(quad_nodes));
127  _materials.push_back(layer_id);
128  }
129  }
130 }
131 
133  std::size_t nLayers)
134 {
135  const unsigned nLayerBoundaries(nLayers - 1);
136  const std::size_t nNodes(layer.getNumberOfNodes());
137  const std::vector<MeshLib::Element*>& layer_elements(layer.getElements());
138  for (MeshLib::Element* elem : layer_elements)
139  {
140  const std::size_t nElemNodes(elem->getNumberOfBaseNodes());
141  for (unsigned i = 0; i < nElemNodes; ++i)
142  {
143  if (elem->getNeighbor(i) == nullptr)
144  {
145  for (unsigned j = 0; j < nLayerBoundaries; ++j)
146  {
147  const std::size_t offset(j * nNodes);
148  MeshLib::Node* n0 = _nodes[offset + getNodeIndex(*elem, i)];
149  MeshLib::Node* n1 =
150  _nodes[offset +
151  getNodeIndex(*elem, (i + 1) % nElemNodes)];
152  MeshLib::Node* n2 =
153  _nodes[offset + nNodes +
154  getNodeIndex(*elem, (i + 1) % nElemNodes)];
155  MeshLib::Node* n3 =
156  _nodes[offset + nNodes + getNodeIndex(*elem, i)];
157 
158  auto const v0 =
159  Eigen::Map<Eigen::Vector3d const>(n0->getCoords());
160  auto const v1 =
161  Eigen::Map<Eigen::Vector3d const>(n1->getCoords());
162  auto const v2 =
163  Eigen::Map<Eigen::Vector3d const>(n2->getCoords());
164  auto const v3 =
165  Eigen::Map<Eigen::Vector3d const>(n3->getCoords());
166  double const eps = std::numeric_limits<double>::epsilon();
167 
168  if ((v2 - v1).norm() > eps)
169  {
170  const std::array tri_nodes = {n0, n2, n1};
171  _elements.push_back(new MeshLib::Tri(tri_nodes));
172  _materials.push_back(nLayers + j);
173  }
174  if ((v3 - v0).norm() > eps)
175  {
176  const std::array tri_nodes = {n0, n3, n2};
177  _elements.push_back(new MeshLib::Tri(tri_nodes));
178  _materials.push_back(nLayers + j);
179  }
180  }
181  }
182  }
183  }
184 }
185 
186 void LayeredVolume::removeCongruentElements(std::size_t nLayers,
187  std::size_t nElementsPerLayer)
188 {
189  for (std::size_t i = nLayers - 1; i > 0; --i)
190  {
191  const std::size_t lower_offset((i - 1) * nElementsPerLayer);
192  const std::size_t upper_offset(i * nElementsPerLayer);
193  for (std::size_t j = 0; j < nElementsPerLayer; ++j)
194  {
195  MeshLib::Element const* const high(_elements[upper_offset + j]);
196  MeshLib::Element* const low(_elements[lower_offset + j]);
197 
198  unsigned count(0);
199  const std::size_t nElemNodes(low->getNumberOfBaseNodes());
200  for (std::size_t k = 0; k < nElemNodes; ++k)
201  {
202  if (getNodeIndex(*high, k) == getNodeIndex(*low, k))
203  {
204  low->setNode(k, _nodes[getNodeIndex(*high, k)]);
205  count++;
206  }
207  }
208 
209  if (count == nElemNodes)
210  {
211  delete _elements[upper_offset + j];
212  // mark element and material entries for deletion
213  _elements[upper_offset + j] = nullptr;
214  _materials[upper_offset + j] = -1;
215  }
216  else
217  {
219  _attribute_points.emplace_back(
220  attr[0],
221  attr[1],
222  (attr[2] + MeshLib::getCenterOfGravity(*low)[2]) / 2.0,
223  _materials[lower_offset + j]);
224  }
225  }
226  }
227  // delete marked entries
228  auto elem_vec_end =
229  std::remove(_elements.begin(), _elements.end(), nullptr);
230  _elements.erase(elem_vec_end, _elements.end());
231  auto mat_vec_end = std::remove(_materials.begin(), _materials.end(), -1);
232  _materials.erase(mat_vec_end, _materials.end());
233 }
Definition of Duplicate functions.
Definition of the LayeredVolume class.
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42
Definition of the MeshLayerMapper class.
Definition of the class Properties that implements a container of properties.
Definition of the Quad class.
Definition of the GeoLib::Raster class.
Definition of the Tri class.
Class Raster is used for managing raster data.
Definition: Raster.h:42
std::vector< int > _materials
static double calcEpsilon(GeoLib::Raster const &low, GeoLib::Raster const &high)
Calculates a data-dependent epsilon value.
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
void removeCongruentElements(std::size_t nLayers, std::size_t nElementsPerLayer)
Removes duplicate 2D elements (possible due to outcroppings)
std::vector< MeshLib::Node > _attribute_points
Definition: LayeredVolume.h:80
void addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer_id, GeoLib::Raster const &raster) override
Adds another layer to the subsurface mesh.
bool createRasterLayers(const MeshLib::Mesh &mesh, const std::vector< GeoLib::Raster const * > &rasters, double minimum_thickness, double noDataReplacementValue=0.0) override
void addLayerBoundaries(const MeshLib::Mesh &layer, std::size_t nLayers)
Creates boundary surfaces between the mapped layers to make the volumes watertight.
const T * getCoords() const
Definition: TemplatePoint.h:75
Element search class.
Definition: ElementSearch.h:28
std::size_t searchByElementType(MeshElemType eleType)
Marks all elements of the given element type.
const std::vector< std::size_t > & getSearchedElementIDs() const
return marked elements
Definition: ElementSearch.h:33
virtual unsigned getNumberOfBaseNodes() const =0
virtual void setNode(unsigned idx, Node *node)=0
static bool layerMapping(MeshLib::Mesh const &mesh, const GeoLib::Raster &raster, double nodata_replacement=0.0, bool const ignore_nodata=false)
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
double norm(MatrixOrVector const &x, MathLib::VecNormType type)
Definition: LinAlg.h:88
MeshLib::Mesh * removeElements(const MeshLib::Mesh &mesh, const std::vector< std::size_t > &removed_element_ids, const std::string &new_mesh_name)
std::vector< Node * > copyNodeVector(const std::vector< Node * > &nodes)
Creates a deep copy of a Node vector.
std::size_t getNodeIndex(Element const &element, unsigned const idx)
Definition: Element.cpp:225
std::vector< Element * > copyElementVector(std::vector< Element * > const &elements, std::vector< Node * > const &new_nodes, std::vector< std::size_t > const *const node_id_map)
MeshLib::Node getCenterOfGravity(Element const &element)
Calculates the center of gravity for the mesh element.
Definition: Element.cpp:126