fixed issue where dem nodes could be located below mesh bottom in prism meshes

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++;
 				}
 

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());