diff --git a/Applications/DataExplorer/DataView/DirectConditionGenerator.cpp b/Applications/DataExplorer/DataView/DirectConditionGenerator.cpp
index 4bbc1edf03a509a82ef8e6cff6819041fef88ce7..329b88a0685be8661631a0d02087bced553f7e26 100644
--- a/Applications/DataExplorer/DataView/DirectConditionGenerator.cpp
+++ b/Applications/DataExplorer/DataView/DirectConditionGenerator.cpp
@@ -13,6 +13,7 @@
  */
 
 #include <fstream>
+#include <memory>
 
 // ThirdParty/logog
 #include "logog/include/logog.hpp"
@@ -61,33 +62,53 @@ const std::vector< std::pair<std::size_t,double> >& DirectConditionGenerator::di
 
 const std::vector< std::pair<std::size_t,double> >& DirectConditionGenerator::directWithSurfaceIntegration(MeshLib::Mesh &mesh, const std::string &filename, double scaling)
 {
-	if (_direct_values.empty())
-	{
-		GeoLib::Raster* raster (FileIO::AsciiRasterInterface::readRaster(filename));
-		if (!raster) {
-			ERR("Error in DirectConditionGenerator::directWithSurfaceIntegration() - could not load raster file.");
-			return _direct_values;
-		}
+	if (!_direct_values.empty()) {
+		ERR(
+		    "Error in DirectConditionGenerator::directWithSurfaceIntegration()"
+		    "- Data vector contains outdated values...");
+		return _direct_values;
+	}
 
-		const MathLib::Vector3 dir(0,0,-1);
-		MeshLib::Mesh* sfc_mesh (MeshLib::MeshSurfaceExtraction::getMeshSurface(mesh, dir, true));
-		std::vector<double> node_area_vec =
-			MeshLib::MeshSurfaceExtraction::getSurfaceAreaForNodes(*sfc_mesh);
-		const std::vector<MeshLib::Node*> &surface_nodes (sfc_mesh->getNodes());
-		const std::size_t nNodes(sfc_mesh->getNNodes());
-		const double no_data (raster->getNoDataValue());
-		_direct_values.reserve(nNodes);
-		for (std::size_t i=0; i<nNodes; ++i)
-		{
-			double val (raster->getValueAtPoint(*surface_nodes[i]));
-			val = (val == no_data) ? 0 : ((val*node_area_vec[i])/scaling);
-			_direct_values.push_back (std::pair<std::size_t, double>(surface_nodes[i]->getID(), val));
-		}
+	std::unique_ptr<GeoLib::Raster> raster(
+	    FileIO::AsciiRasterInterface::readRaster(filename));
+	if (!raster) {
+		ERR(
+		    "Error in DirectConditionGenerator::directWithSurfaceIntegration()"
+		    "- could not load raster file.");
+		return _direct_values;
+	}
 
-		delete raster;
+	MathLib::Vector3 const dir(0.0, 0.0, -1.0);
+	double const angle(90);
+	std::string const prop_name("OriginalSubsurfaceNodeIDs");
+	std::unique_ptr<MeshLib::Mesh> surface_mesh(
+	    MeshLib::MeshSurfaceExtraction::getMeshSurface(
+	        mesh, dir, angle, prop_name));
+
+	std::vector<double> node_area_vec =
+		MeshLib::MeshSurfaceExtraction::getSurfaceAreaForNodes(*surface_mesh);
+	const std::vector<MeshLib::Node*> &surface_nodes(surface_mesh->getNodes());
+	const std::size_t nNodes(surface_mesh->getNNodes());
+	const double no_data(raster->getNoDataValue());
+
+	boost::optional<MeshLib::PropertyVector<int> const&> opt_node_id_pv(
+	    surface_mesh->getProperties().getPropertyVector<int>(prop_name));
+	if (!opt_node_id_pv) {
+		ERR(
+		    "Need subsurface node ids, but the property \"%s\" is not "
+		    "available.",
+		    prop_name.c_str());
+		return _direct_values;
+	}
+
+	MeshLib::PropertyVector<int> const& node_id_pv(*opt_node_id_pv);
+	_direct_values.reserve(nNodes);
+	for (std::size_t i=0; i<nNodes; ++i)
+	{
+		double val(raster->getValueAtPoint(*surface_nodes[i]));
+		val = (val == no_data) ? 0 : ((val*node_area_vec[i])/scaling);
+		_direct_values.push_back(std::pair<std::size_t, double>(node_id_pv[i], val));
 	}
-	else
-		std::cout << "Error in DirectConditionGenerator::directWithSurfaceIntegration() - Data vector contains outdated values..." << std::endl;
 
 	return _direct_values;
 }
diff --git a/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp b/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp
index 82239b3fa14e7bc642d782e12d242a55d89d9e29..dfe202cddde8999f66d3973f359e73b738e3919b 100644
--- a/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp
+++ b/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp
@@ -179,10 +179,8 @@ int main (int argc, char* argv[])
 	const MathLib::Vector3 dir(0,0,-1);
 	double const angle(90);
 	std::unique_ptr<MeshLib::Mesh> surface_mesh(
-		MeshLib::MeshSurfaceExtraction::getMeshSurface(
-			*subsurface_mesh, dir, angle, false
-		)
-	);
+	    MeshLib::MeshSurfaceExtraction::getMeshSurface(*subsurface_mesh, dir,
+	                                                   angle));
 	INFO("done.");
 	delete subsurface_mesh;
 	subsurface_mesh = nullptr;
diff --git a/MeshLib/MeshEditing/AddLayerToMesh.cpp b/MeshLib/MeshEditing/AddLayerToMesh.cpp
index 07736e29863bd70dc86723889388b61154d4da39..ac48a2238f391b124c4aab0e7e6ab64dc454817e 100644
--- a/MeshLib/MeshEditing/AddLayerToMesh.cpp
+++ b/MeshLib/MeshEditing/AddLayerToMesh.cpp
@@ -13,9 +13,10 @@
 
 #include "AddLayerToMesh.h"
 
-#include <vector>
 #include <map>
 #include <memory>
+#include <numeric>
+#include <vector>
 
 #include "logog/include/logog.hpp"
 
@@ -29,8 +30,22 @@
 namespace MeshLib
 {
 
+/** Extrudes point, line, triangle or quad elements to its higher dimensional
+ * versions, i.e. line, quad, prism, hexahedron.
+ *
+ * @param subsfc_nodes the nodes the elements are based on
+ * @param sfc_elem the element of the surface that will be extruded
+ * @param sfc_to_subsfc_id_map relation between the surface nodes of the surface
+ * element and the ids of the nodes of the subsurface mesh
+ * @param subsfc_sfc_id_map mapping of the surface nodes of the current mesh
+ * to the surface nodes of the extruded mesh
+ *
+ * @return extruded element (point -> line, line -> quad, tri -> prism, quad ->
+ * hexahedron)
+*/
 MeshLib::Element* extrudeElement(std::vector<MeshLib::Node*> const& subsfc_nodes,
 	MeshLib::Element const& sfc_elem,
+	MeshLib::PropertyVector<std::size_t> const& sfc_to_subsfc_id_map,
 	std::map<std::size_t, std::size_t> const& subsfc_sfc_id_map)
 {
 	if (sfc_elem.getDimension() > 2)
@@ -41,10 +56,11 @@ MeshLib::Element* extrudeElement(std::vector<MeshLib::Node*> const& subsfc_nodes
 
 	for (unsigned j=0; j<nElemNodes; ++j)
 	{
-		new_nodes[j] = subsfc_nodes[sfc_elem.getNode(j)->getID()];
+		std::size_t const subsfc_id(
+		    sfc_to_subsfc_id_map[sfc_elem.getNode(j)->getID()]);
+		new_nodes[j] = subsfc_nodes[subsfc_id];
 		std::size_t new_idx = (nElemNodes==2) ? (3-j) : (nElemNodes+j);
-		new_nodes[new_idx] =
-		    subsfc_nodes[subsfc_sfc_id_map.at(sfc_elem.getNode(j)->getID())];
+		new_nodes[new_idx] = subsfc_nodes[subsfc_sfc_id_map.at(subsfc_id)];
 	}
 	
 	if (sfc_elem.getGeomType() == MeshLib::MeshElemType::LINE)
@@ -80,12 +96,27 @@ MeshLib::Mesh* addLayerToMesh(MeshLib::Mesh const& mesh, double thickness,
 	const MathLib::Vector3 dir(0, 0, flag);
 	double const angle(90);
 	std::unique_ptr<MeshLib::Mesh> sfc_mesh (nullptr);
-	
+
+	std::string const prop_name("OriginalSubsurfaceNodeIDs");
+
 	if (mesh.getDimension() == 3)
-		sfc_mesh.reset(MeshLib::MeshSurfaceExtraction::getMeshSurface(mesh, dir, angle, true));
-	else
+		sfc_mesh.reset(MeshLib::MeshSurfaceExtraction::getMeshSurface(
+		    mesh, dir, angle, prop_name));
+	else {
 		sfc_mesh = (on_top) ? std::unique_ptr<MeshLib::Mesh>(new MeshLib::Mesh(mesh)) :
 		                      std::unique_ptr<MeshLib::Mesh>(MeshLib::createFlippedMesh(mesh));
+		// add property storing node ids
+		boost::optional<MeshLib::PropertyVector<std::size_t>&> pv(
+		    sfc_mesh->getProperties().createNewPropertyVector<std::size_t>(
+		        prop_name, MeshLib::MeshItemType::Node, 1));
+		if (pv) {
+			pv->resize(sfc_mesh->getNNodes());
+			std::iota(pv->begin(), pv->end(), 0);
+		} else {
+			ERR("Could not create and initialize property.");
+			return nullptr;
+		}
+	}
 	INFO("done.");
 
 	// *** add new surface nodes
@@ -99,24 +130,36 @@ MeshLib::Mesh* addLayerToMesh(MeshLib::Mesh const& mesh, double thickness,
 	std::vector<MeshLib::Node*> const& sfc_nodes(sfc_mesh->getNodes());
 	std::size_t const n_sfc_nodes(sfc_nodes.size());
 
+	// fetch subsurface node ids PropertyVector
+	boost::optional<MeshLib::PropertyVector<std::size_t> const&> opt_node_id_pv(
+	    sfc_mesh->getProperties().getPropertyVector<std::size_t>(prop_name));
+	if (!opt_node_id_pv) {
+		ERR(
+		    "Need subsurface node ids, but the property \"%s\" is not "
+		    "available.",
+		    prop_name.c_str());
+		return nullptr;
+	}
+
+	MeshLib::PropertyVector<std::size_t> const& node_id_pv(*opt_node_id_pv);
 	// *** copy sfc nodes to subsfc mesh node
 	std::map<std::size_t, std::size_t> subsfc_sfc_id_map;
 	for (std::size_t k(0); k<n_sfc_nodes; ++k) {
-		std::size_t const subsfc_id(sfc_nodes[k]->getID());
+		std::size_t const subsfc_id(node_id_pv[k]);
 		std::size_t const sfc_id(k+n_subsfc_nodes);
 		subsfc_sfc_id_map.insert(std::make_pair(subsfc_id, sfc_id));
-		MeshLib::Node const& node (*sfc_nodes[k]);
+		MeshLib::Node const& node(*sfc_nodes[k]);
 		subsfc_nodes.push_back(new MeshLib::Node(
 		    node[0], node[1], node[2] - (flag * thickness), sfc_id));
 	}
 
-	// *** insert new top layer elements into subsfc_mesh
+	// *** insert new layer elements into subsfc_mesh
 	std::vector<MeshLib::Element*> const& sfc_elements(sfc_mesh->getElements());
 	std::size_t const n_sfc_elements(sfc_elements.size());
 	for (std::size_t k(0); k<n_sfc_elements; ++k)
-		subsfc_elements.push_back(
-			extrudeElement(subsfc_nodes, *sfc_elements[k], subsfc_sfc_id_map)
-		);
+		subsfc_elements.push_back(extrudeElement(subsfc_nodes, *sfc_elements[k],
+		                                         node_id_pv,
+		                                         subsfc_sfc_id_map));
 
 	auto new_mesh = new MeshLib::Mesh(name, subsfc_nodes, subsfc_elements);
 
diff --git a/MeshLib/MeshSurfaceExtraction.cpp b/MeshLib/MeshSurfaceExtraction.cpp
index a2a18d14917a248bd7abf16417ddcfd40b822ae5..4b7807636cf1087e646d302e9c737fae3a23f93b 100644
--- a/MeshLib/MeshSurfaceExtraction.cpp
+++ b/MeshLib/MeshSurfaceExtraction.cpp
@@ -69,7 +69,9 @@ std::vector<double> MeshSurfaceExtraction::getSurfaceAreaForNodes(const MeshLib:
 	return node_area_vec;
 }
 
-MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, double angle, bool keepOriginalNodeIds)
+MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(
+    const MeshLib::Mesh& mesh, const MathLib::Vector3& dir, double angle,
+    std::string const& subsfc_node_id_backup_prop_name)
 {
 	if (angle < 0 || angle > 90)
 	{
@@ -107,17 +109,23 @@ MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(const MeshLib::Mesh &mesh,
 	}
 
 	std::vector<std::size_t> id_map;
-	if (keepOriginalNodeIds)
+	if (!subsfc_node_id_backup_prop_name.empty())
 	{
 		id_map.reserve(sfc_nodes.size());
 		for (auto node = sfc_nodes.cbegin(); node != sfc_nodes.cend(); ++node)
 			id_map.push_back((*node)->getID());
 	}
 	MeshLib::Mesh* result (new Mesh(mesh.getName()+"-Surface", sfc_nodes, new_elements));
-	if (keepOriginalNodeIds)
-		for (auto node = sfc_nodes.cbegin(); node != sfc_nodes.cend(); ++node)
-			(*node)->setID(id_map[(*node)->getID()]);
-
+	// transmit the original node ids of the subsurface mesh as a property
+	if (!subsfc_node_id_backup_prop_name.empty()) {
+		boost::optional<MeshLib::PropertyVector<std::size_t>&> orig_node_ids(
+		    result->getProperties().createNewPropertyVector<std::size_t>(
+		        subsfc_node_id_backup_prop_name , MeshLib::MeshItemType::Node, 1));
+		if (orig_node_ids) {
+			orig_node_ids->resize(id_map.size());
+			std::copy(id_map.cbegin(), id_map.cend(), orig_node_ids->begin());
+		}
+	}
 	return result;
 }
 
diff --git a/MeshLib/MeshSurfaceExtraction.h b/MeshLib/MeshSurfaceExtraction.h
index e0dcd6484459a56c2e175a201a24d343f7c6e409..f72789c5efbb318c7cb90008d5f844f7707ab7c6 100644
--- a/MeshLib/MeshSurfaceExtraction.h
+++ b/MeshLib/MeshSurfaceExtraction.h
@@ -45,12 +45,21 @@ public:
 	/**
 	 * Returns the 2d-element mesh representing the surface of the given mesh.
 	 * \param mesh                The original mesh
-	 * \param dir                 The direction in which face normals have to point to be considered surface elements
-	 * \param angle               The angle of the allowed deviation from the given direction (0 <= angle <= 90 degrees)
-	 * \param keepOriginalNodeIds If true, ids of mesh nodes are set to ids in original mesh, otherwise node ids are reset (as usual when creating a mesh)
-	 * \return                    A 2D mesh representing the surface in direction dir
+	 * \param dir                 The direction in which face normals have to
+	 * point to be considered surface elements
+	 * \param angle               The angle of the allowed deviation from the
+	 * given direction (0 <= angle <= 90 degrees)
+	 * \param subsfc_node_id_backup_prop_name Name of the property in the
+	 * surface mesh the subsurface node ids are stored to. This is a mapping
+	 * surface node -> subsurface node which is needed from some applications.
+	 * If the string is empty, there isn't a backup created.
+	 * \return A 2D mesh representing the surface in direction dir
 	 */
-	static MeshLib::Mesh* getMeshSurface(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, double angle, bool keepOriginalNodeIds = false);
+	static MeshLib::Mesh* getMeshSurface(
+	    const MeshLib::Mesh& mesh,
+	    const MathLib::Vector3& dir,
+	    double angle,
+	    std::string const& subsfc_node_id_backup_prop_name = "");
 
 	/**
 	 * Returns the boundary of mesh, i.e. lines for 2D meshes and surfaces for 3D meshes.
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index dbf35a623e5c41525e64d4bb5fa8821402c6a412..3c457be50f904cc13543144e9acbab72b095c280 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -33,7 +33,6 @@ class Node : public MathLib::Point3dWithID
 {
 	/* friend classes: */
 	friend class Mesh;
-	friend class MeshSurfaceExtraction;
 	friend class MeshRevision;
 	friend class MeshLayerMapper;