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