diff --git a/MeshLib/MeshGenerators/LayeredVolume.cpp b/MeshLib/MeshGenerators/LayeredVolume.cpp
index 3263dcb0343c30bc69accbfb9605e2c0e53c5409..7dca13f4f08e49e3f1e68531fd7bffaae445c70c 100644
--- a/MeshLib/MeshGenerators/LayeredVolume.cpp
+++ b/MeshLib/MeshGenerators/LayeredVolume.cpp
@@ -69,9 +69,9 @@ void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned lay
for (std::size_t i=0; i<nNodes; ++i)
{
- if (layer_id > 0 &&
- ((*layer_nodes[i])[2] == no_data_value ||
- (*_nodes[last_layer_offset+i])[2]-(*layer_nodes[i])[2] < _elevation_epsilon))
+ if ((layer_id > 0) &&
+ ((std::abs((*layer_nodes[i])[2] - no_data_value) < std::numeric_limits<double>::epsilon()) ||
+ ((*_nodes[last_layer_offset+i])[2]-(*layer_nodes[i])[2] < _elevation_epsilon)))
_nodes.push_back(new MeshLib::Node(*_nodes[last_layer_offset+i]));
else
_nodes.push_back(new MeshLib::Node(layer_nodes[i]->getCoords(), _nodes.size()));
@@ -83,16 +83,16 @@ void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned lay
if (elem->getGeomType() == MeshElemType::TRIANGLE)
{
std::array<MeshLib::Node*,3> tri_nodes = {{ _nodes[node_id_offset+elem->getNodeIndex(0)],
- _nodes[node_id_offset+elem->getNodeIndex(1)],
- _nodes[node_id_offset+elem->getNodeIndex(2)] }};
+ _nodes[node_id_offset+elem->getNodeIndex(1)],
+ _nodes[node_id_offset+elem->getNodeIndex(2)] }};
_elements.push_back(new MeshLib::Tri(tri_nodes, layer_id+1));
}
else if (elem->getGeomType() == MeshElemType::QUAD)
{
std::array<MeshLib::Node*,4> quad_nodes = {{ _nodes[node_id_offset+elem->getNodeIndex(0)],
- _nodes[node_id_offset+elem->getNodeIndex(1)],
- _nodes[node_id_offset+elem->getNodeIndex(2)],
- _nodes[node_id_offset+elem->getNodeIndex(3)] }};
+ _nodes[node_id_offset+elem->getNodeIndex(1)],
+ _nodes[node_id_offset+elem->getNodeIndex(2)],
+ _nodes[node_id_offset+elem->getNodeIndex(3)] }};
_elements.push_back(new MeshLib::Quad(quad_nodes, layer_id+1));
}
}
@@ -138,15 +138,15 @@ void LayeredVolume::removeCongruentElements(std::size_t nLayers, std::size_t nEl
const std::size_t lower_offset ( i * nElementsPerLayer);
for (std::size_t j=0; j<nElementsPerLayer; ++j)
{
- MeshLib::Element const*const high (_elements[upper_offset+j]);
- MeshLib::Element *low (_elements[lower_offset+j]);
+ MeshLib::Element *const high (_elements[upper_offset+j]);
+ MeshLib::Element const*const low (_elements[lower_offset+j]);
unsigned count(0);
const std::size_t nElemNodes (low->getNNodes());
for (std::size_t k=0; k<nElemNodes; ++k)
if (high->getNodeIndex(k) == low->getNodeIndex(k))
{
- low->setNode(k, _nodes[high->getNodeIndex(k)]);
+ high->setNode(k, _nodes[low->getNodeIndex(k)]);
count++;
}
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.cpp b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
index c56b846d7089d364e21ead757fad6eb349914d5a..656ca3b7bff3a0e78632f8fcc8b6549ef26f9ce9 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.cpp
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
@@ -108,30 +108,38 @@ MeshLib::Mesh* MeshLayerMapper::createStaticLayers(MeshLib::Mesh const& mesh, st
bool MeshLayerMapper::createRasterLayers(MeshLib::Mesh const& mesh, std::vector<GeoLib::Raster const*> const& rasters, double noDataReplacementValue)
{
- const std::size_t nLayers(rasters.size());
- if (nLayers < 1 || mesh.getDimension() != 2)
- {
- ERR("MeshLayerMapper::createRasterLayers(): A 2D mesh and at least two rasters required as input.");
- return nullptr;
- }
+ const std::size_t nLayers(rasters.size());
+ if (nLayers < 1 || mesh.getDimension() != 2)
+ {
+ ERR("MeshLayerMapper::createRasterLayers(): A 2D mesh and at least two rasters required as input.");
+ return nullptr;
+ }
- MeshLib::Mesh* dem_mesh (new MeshLib::Mesh(mesh));
- if (layerMapping(*dem_mesh, *rasters.back(), 0))
- {
- std::size_t const nNodes = mesh.getNNodes();
- _nodes.reserve(nLayers * nNodes);
+ MeshLib::Mesh* top (new MeshLib::Mesh(mesh));
+ if (!layerMapping(*top, *rasters.back(), 0))
+ return false;
+
+ MeshLib::Mesh* bottom (new MeshLib::Mesh(mesh));
+ if (!layerMapping(*bottom, *rasters[0], 0))
+ return false;
- // number of triangles in the original mesh
- std::size_t const nElems (std::count_if(mesh.getElements().begin(), mesh.getElements().end(),
- [](MeshLib::Element const* elem) { return (elem->getGeomType() == MeshElemType::TRIANGLE);}));
- _elements.reserve(nElems * (nLayers-1));
+ std::size_t const nNodes = mesh.getNNodes();
+ _nodes.reserve(nLayers * nNodes);
- for (std::size_t i=0; i<nLayers; ++i)
- addLayerToMesh(*dem_mesh, i, *rasters[i]);
+ // number of triangles in the original mesh
+ std::size_t const nElems (std::count_if(mesh.getElements().begin(), mesh.getElements().end(),
+ [](MeshLib::Element const* elem) { return (elem->getGeomType() == MeshElemType::TRIANGLE);}));
+ _elements.reserve(nElems * (nLayers-1));
- return true;
- }
- return false;
+ // bottom layer
+ std::vector<MeshLib::Node*> const& nodes = bottom->getNodes();
+ std::copy(nodes.cbegin(), nodes.cend(), std::back_inserter(_nodes));
+
+ // the other layers
+ for (std::size_t i=1; i<nLayers; ++i)
+ addLayerToMesh(*top, i, *rasters[i]);
+
+ return true;
}
void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer_id, GeoLib::Raster const& raster)
@@ -147,17 +155,13 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
// min of dem elevation and layer elevation
double const elevation = std::min(raster.interpolateValueAtPoint(*nodes[i]), (*nodes[i])[2]);
- if ((layer_id > 0) &&
- ((std::abs(elevation - no_data_value) < std::numeric_limits<double>::epsilon()) ||
- (elevation - (*_nodes[last_layer_node_offset + i])[2] < std::numeric_limits<double>::epsilon())))
+ 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));
}
- if (layer_id == 0)
- return;
-
std::vector<MeshLib::Element*> const& elems = dem_mesh.getElements();
std::size_t const nElems (dem_mesh.getNElements());