resolved dem-nodes below bottom layer for tet mesh generation

Conflicts:
	MeshLib/MeshGenerators/LayeredVolume.cpp

MeshLib/MeshGenerators/LayeredVolume.cpp

@@ -18,10 +18,11 @@
 
 #include "Raster.h"
 
-#include "Elements/Tri.h"
-#include "Elements/Quad.h"
-#include "MeshGenerators/MeshLayerMapper.h"
-#include "MeshEditing/ElementExtraction.h"
+#include "MeshLib/Elements/Tri.h"
+#include "MeshLib/Elements/Quad.h"
+#include "MeshLib/MeshEditing/ElementExtraction.h"
+#include "MeshLib/MeshEditing/DuplicateMeshComponents.h"
+#include "MeshLib/MeshGenerators/MeshLayerMapper.h"
 
 
 bool LayeredVolume::createRasterLayers(const MeshLib::Mesh &mesh, const std::vector<GeoLib::Raster const*> &rasters, double noDataReplacementValue)
@@ -32,52 +33,61 @@ bool LayeredVolume::createRasterLayers(const MeshLib::Mesh &mesh, const std::vec
 	_elevation_epsilon = calcEpsilon(*rasters[0], *rasters.back());
 	if (_elevation_epsilon <= 0)
 		return false;
-
+	
 	// remove line elements, only tri + quad remain
 	MeshLib::ElementExtraction ex(mesh);
 	ex.searchByElementType(MeshElemType::LINE);
-	MeshLib::Mesh* mesh_layer (ex.removeMeshElements("MeshLayer"));
-	if (mesh_layer==nullptr)
-		mesh_layer = new MeshLib::Mesh(mesh);
+	MeshLib::Mesh* top (ex.removeMeshElements("MeshLayer"));
+	if (top==nullptr)
+		top = new MeshLib::Mesh(mesh);
 
-	// map each layer and attach to subsurface mesh
-	const std::size_t nRasters (rasters.size());
-	for (size_t i=0; i<nRasters; ++i)
+	if (!MeshLayerMapper::layerMapping(*top, *rasters[0], noDataReplacementValue))
+		return false;
+
+	MeshLib::Mesh* bottom (new MeshLib::Mesh(*top));
+	if (!MeshLayerMapper::layerMapping(*bottom, *rasters.back(), 0))
 	{
-		const double replacement_value = (i==(nRasters-1)) ? noDataReplacementValue : rasters[i]->getNoDataValue();
-		if (!MeshLayerMapper::layerMapping(*mesh_layer, *rasters[i], replacement_value))
-		{
-			this->cleanUpOnError();
-			return false;
-		}
-		this->addLayerToMesh(*mesh_layer, i, *rasters[i]);
+		delete top;
+		return false;
 	}
+
+	_nodes = MeshLib::copyNodeVector(bottom->getNodes());
+	_elements = MeshLib::copyElementVector(bottom->getElements(), _nodes);
+	delete bottom;
+
+	// map each layer and attach to subsurface mesh
+	const std::size_t nRasters (rasters.size());
+	for (int i=nRasters-2; i>=0; --i)
+		this->addLayerToMesh(*top, nRasters-i-1, *rasters[i]);
+
 	// close boundaries between layers
-	this->addLayerBoundaries(*mesh_layer, nRasters);
-	this->removeCongruentElements(nRasters, mesh_layer->getNElements());
-	delete mesh_layer;
+	this->addLayerBoundaries(*top, nRasters);
+	this->removeCongruentElements(nRasters, top->getNElements());
+	delete top;
 	return true;
 }
 
-void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned layer_id, GeoLib::Raster const& raster)
+void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer_id, GeoLib::Raster const& raster)
 {
-	const std::vector<MeshLib::Node*> &layer_nodes (mesh_layer.getNodes());
-	const std::size_t nNodes (layer_nodes.size());
+	const std::size_t nNodes (dem_mesh.getNNodes());
+	const std::vector<MeshLib::Node*> const& nodes (dem_mesh.getNodes());
 	const std::size_t node_id_offset (_nodes.size());
 	const std::size_t last_layer_offset (node_id_offset-nNodes);
 	const double no_data_value (raster.getNoDataValue());
 
 	for (std::size_t i=0; i<nNodes; ++i)
 	{
-		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)))
+		// 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(layer_nodes[i]->getCoords(), _nodes.size()));
+			_nodes.push_back(new MeshLib::Node((*nodes[i])[0], (*nodes[i])[1], elevation, _nodes.size()));
 	}
 
-	const std::vector<MeshLib::Element*> &layer_elements (mesh_layer.getElements());
+	const std::vector<MeshLib::Element*> &layer_elements (dem_mesh.getElements());
 	for (MeshLib::Element* elem : layer_elements)
 	{
 		if (elem->getGeomType() == MeshElemType::TRIANGLE)
@@ -85,7 +95,7 @@ void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned lay
 			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)] }};
-			_elements.push_back(new MeshLib::Tri(tri_nodes, layer_id+1));
+			_elements.push_back(new MeshLib::Tri(tri_nodes, layer_id));
 		}
 		else if (elem->getGeomType() == MeshElemType::QUAD)
 		{
@@ -93,7 +103,7 @@ void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned lay
 			                                             _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));
+			_elements.push_back(new MeshLib::Quad(quad_nodes, layer_id));
 		}
 	}
 }
@@ -132,21 +142,21 @@ void LayeredVolume::addLayerBoundaries(const MeshLib::Mesh &layer, std::size_t n
 
 void LayeredVolume::removeCongruentElements(std::size_t nLayers, std::size_t nElementsPerLayer)
 {
-	for (std::size_t i=1; i<nLayers; ++i)
+	for (std::size_t i=nLayers-1; i>0; --i)
 	{
-		const std::size_t upper_offset ((i-1) * nElementsPerLayer);
-		const std::size_t lower_offset ( i    * nElementsPerLayer);
+		const std::size_t lower_offset ((i-1) * nElementsPerLayer);
+		const std::size_t upper_offset ( i    * nElementsPerLayer);
 		for (std::size_t j=0; j<nElementsPerLayer; ++j)
 		{
-			MeshLib::Element *const high (_elements[upper_offset+j]);
-			MeshLib::Element const*const low  (_elements[lower_offset+j]);
+			MeshLib::Element const*const high (_elements[upper_offset+j]);
+			MeshLib::Element *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))
 				{
-					high->setNode(k, _nodes[low->getNodeIndex(k)]);
+					low->setNode(k, _nodes[high->getNodeIndex(k)]);
 					count++;
 				}