diff --git a/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp b/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp
index d16106a875b0711133a68f293881369d6c8df60b..6303dc9fa47a588983039bcab50bbd107b29662d 100644
--- a/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp
+++ b/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp
@@ -63,6 +63,21 @@ double LayeredMeshGenerator::calcEpsilon(GeoLib::Raster const& high, GeoLib::Ras
     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
 {
     for (auto raster = raster_paths.begin(); raster != raster_paths.end(); ++raster)
diff --git a/MeshLib/MeshGenerators/LayeredMeshGenerator.h b/MeshLib/MeshGenerators/LayeredMeshGenerator.h
index 9cd3f81599642af2a2099426aac00241ae637fb8..5152483da1806fffa55a2489f9f78df0e823ab3f 100644
--- a/MeshLib/MeshGenerators/LayeredMeshGenerator.h
+++ b/MeshLib/MeshGenerators/LayeredMeshGenerator.h
@@ -62,6 +62,12 @@ protected:
     /// Adds another layer to the subsurface mesh
     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
     double calcEpsilon(GeoLib::Raster const& high, GeoLib::Raster const& low);
 
diff --git a/MeshLib/MeshGenerators/LayeredVolume.cpp b/MeshLib/MeshGenerators/LayeredVolume.cpp
index cc36b1e745376fdd773750d84c331fe964de805d..1bbaba483d92040a2fb27d5e8f9f8f9a8f6c0ed8 100644
--- a/MeshLib/MeshGenerators/LayeredVolume.cpp
+++ b/MeshLib/MeshGenerators/LayeredVolume.cpp
@@ -72,20 +72,10 @@ void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer
 	const std::size_t nNodes (dem_mesh.getNNodes());
 	const std::vector<MeshLib::Node*> &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());
+	const std::size_t last_layer_node_offset (node_id_offset-nNodes);
 
 	for (std::size_t i=0; i<nNodes; ++i)
-	{
-		// 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()));
-	}
+		_nodes.push_back(getNewLayerNode(*nodes[i], *_nodes[last_layer_node_offset + i], raster, _nodes.size(), _elevation_epsilon));
 
 	const std::vector<MeshLib::Element*> &layer_elements (dem_mesh.getElements());
 	for (MeshLib::Element* elem : layer_elements)
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.cpp b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
index aaf22cb0d2501871ba3cbf815c6e1dd2c3f03ddb..79ae86356e59c11cf0662a96b74cc27b0831ed41 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.cpp
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
@@ -153,20 +153,10 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
     std::size_t const nNodes = dem_mesh.getNNodes();
     std::vector<MeshLib::Node*> const& nodes = dem_mesh.getNodes();
     int const last_layer_node_offset = (layer_id-1) * nNodes;
-    double const no_data_value (raster.getNoDataValue());
 
     // add nodes for new layer
     for (std::size_t i=0; i<nNodes; ++i)
-    {
-        // 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));
-    }
+        _nodes.push_back(getNewLayerNode(*nodes[i], *_nodes[last_layer_node_offset + i], raster, _nodes.size()));
 
     std::vector<MeshLib::Element*> const& elems = dem_mesh.getElements();
     std::size_t const nElems (dem_mesh.getNElements());