Skip to content
Snippets Groups Projects
Commit 8b801e76 authored by Karsten Rink's avatar Karsten Rink
Browse files

moved node construction for new layers to base class

parent 9ff5e9e7
No related branches found
No related tags found
No related merge requests found
...@@ -63,6 +63,21 @@ double LayeredMeshGenerator::calcEpsilon(GeoLib::Raster const& high, GeoLib::Ras ...@@ -63,6 +63,21 @@ double LayeredMeshGenerator::calcEpsilon(GeoLib::Raster const& high, GeoLib::Ras
return ((max-min)*1e-07); return ((max-min)*1e-07);
} }
MeshLib::Node* LayeredMeshGenerator::getNewLayerNode(MeshLib::Node const& dem_node,
MeshLib::Node const& last_layer_node,
GeoLib::Raster const& raster,
std::size_t new_node_id,
double minimum_thickness) const
{
double const elevation = std::min(raster.interpolateValueAtPoint(dem_node), dem_node[2]);
if ((std::abs(elevation - raster.getNoDataValue()) < std::numeric_limits<double>::epsilon()) ||
(elevation - last_layer_node[2] < minimum_thickness))
return new MeshLib::Node(last_layer_node);
else
return new MeshLib::Node(dem_node[0], dem_node[1], elevation, new_node_id);
}
bool LayeredMeshGenerator::allRastersExist(std::vector<std::string> const& raster_paths) const bool LayeredMeshGenerator::allRastersExist(std::vector<std::string> const& raster_paths) const
{ {
for (auto raster = raster_paths.begin(); raster != raster_paths.end(); ++raster) for (auto raster = raster_paths.begin(); raster != raster_paths.end(); ++raster)
......
...@@ -62,6 +62,12 @@ protected: ...@@ -62,6 +62,12 @@ protected:
/// Adds another layer to the subsurface mesh /// Adds another layer to the subsurface mesh
virtual void addLayerToMesh(MeshLib::Mesh const& mesh_layer, unsigned layer_id, GeoLib::Raster const& raster) = 0; virtual void addLayerToMesh(MeshLib::Mesh const& mesh_layer, unsigned layer_id, GeoLib::Raster const& raster) = 0;
MeshLib::Node* getNewLayerNode(MeshLib::Node const& dem_node,
MeshLib::Node const& last_layer_node,
GeoLib::Raster const& raster,
std::size_t new_node_id,
double minimum_thickness = std::numeric_limits<double>::epsilon()) const;
/// Calculates a data-dependent epsilon value /// Calculates a data-dependent epsilon value
double calcEpsilon(GeoLib::Raster const& high, GeoLib::Raster const& low); double calcEpsilon(GeoLib::Raster const& high, GeoLib::Raster const& low);
......
...@@ -72,20 +72,10 @@ void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer ...@@ -72,20 +72,10 @@ void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer
const std::size_t nNodes (dem_mesh.getNNodes()); const std::size_t nNodes (dem_mesh.getNNodes());
const std::vector<MeshLib::Node*> &nodes (dem_mesh.getNodes()); const std::vector<MeshLib::Node*> &nodes (dem_mesh.getNodes());
const std::size_t node_id_offset (_nodes.size()); const std::size_t node_id_offset (_nodes.size());
const std::size_t last_layer_offset (node_id_offset-nNodes); const std::size_t last_layer_node_offset (node_id_offset-nNodes);
const double no_data_value (raster.getNoDataValue());
for (std::size_t i=0; i<nNodes; ++i) for (std::size_t i=0; i<nNodes; ++i)
{ _nodes.push_back(getNewLayerNode(*nodes[i], *_nodes[last_layer_node_offset + i], raster, _nodes.size(), _elevation_epsilon));
// min of dem elevation and layer elevation
double const elevation = std::min(raster.interpolateValueAtPoint(*nodes[i]), (*nodes[i])[2]);
if ((std::abs(elevation - no_data_value) < std::numeric_limits<double>::epsilon()) ||
(elevation - (*_nodes[last_layer_offset+i])[2] < _elevation_epsilon))
_nodes.push_back(new MeshLib::Node(*_nodes[last_layer_offset+i]));
else
_nodes.push_back(new MeshLib::Node((*nodes[i])[0], (*nodes[i])[1], elevation, _nodes.size()));
}
const std::vector<MeshLib::Element*> &layer_elements (dem_mesh.getElements()); const std::vector<MeshLib::Element*> &layer_elements (dem_mesh.getElements());
for (MeshLib::Element* elem : layer_elements) for (MeshLib::Element* elem : layer_elements)
......
...@@ -153,20 +153,10 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay ...@@ -153,20 +153,10 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
std::size_t const nNodes = dem_mesh.getNNodes(); std::size_t const nNodes = dem_mesh.getNNodes();
std::vector<MeshLib::Node*> const& nodes = dem_mesh.getNodes(); std::vector<MeshLib::Node*> const& nodes = dem_mesh.getNodes();
int const last_layer_node_offset = (layer_id-1) * nNodes; int const last_layer_node_offset = (layer_id-1) * nNodes;
double const no_data_value (raster.getNoDataValue());
// add nodes for new layer // add nodes for new layer
for (std::size_t i=0; i<nNodes; ++i) for (std::size_t i=0; i<nNodes; ++i)
{ _nodes.push_back(getNewLayerNode(*nodes[i], *_nodes[last_layer_node_offset + i], raster, _nodes.size()));
// min of dem elevation and layer elevation
double const elevation = std::min(raster.interpolateValueAtPoint(*nodes[i]), (*nodes[i])[2]);
if ((std::abs(elevation - no_data_value) < std::numeric_limits<double>::epsilon()) ||
(elevation - (*_nodes[last_layer_node_offset + i])[2] < std::numeric_limits<double>::epsilon()))
_nodes.push_back(new MeshLib::Node((*nodes[i])[0], (*nodes[i])[1], elevation, _nodes[last_layer_node_offset +i]->getID()));
else
_nodes.push_back(new MeshLib::Node((*nodes[i])[0], (*nodes[i])[1], elevation, (layer_id * nNodes) + i));
}
std::vector<MeshLib::Element*> const& elems = dem_mesh.getElements(); std::vector<MeshLib::Element*> const& elems = dem_mesh.getElements();
std::size_t const nElems (dem_mesh.getNElements()); std::size_t const nElems (dem_mesh.getNElements());
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment