diff --git a/Applications/DataExplorer/DataView/MeshLayerEditDialog.cpp b/Applications/DataExplorer/DataView/MeshLayerEditDialog.cpp
index 7901fa7e79fbc185a70457115a9d4fa0f450760e..c83cfc3e63e0cef71b8ab92d79a9b55cb5c96325 100644
--- a/Applications/DataExplorer/DataView/MeshLayerEditDialog.cpp
+++ b/Applications/DataExplorer/DataView/MeshLayerEditDialog.cpp
@@ -195,7 +195,6 @@ MeshLib::Mesh* MeshLayerEditDialog::createPrismMesh()
 	const unsigned nLayers = _layerEdit->text().toInt();
 
 	MeshLib::MeshLayerMapper mapper;
-	MeshLib::Mesh* new_mesh (nullptr);
 
 	QTime myTimer0;
 	myTimer0.start();
@@ -207,18 +206,20 @@ MeshLib::Mesh* MeshLayerEditDialog::createPrismMesh()
 		for (int i=nLayers; i>=0; --i)
 			raster_paths.push_back(this->_edits[i]->text().toStdString());
 		if (mapper.createLayers(*_msh, raster_paths, minimum_thickness))
-			new_mesh= mapper.getMesh("SubsurfaceMesh");
+		{
+			INFO("Mesh construction time: %d ms.", myTimer0.elapsed());
+			return mapper.getMesh("SubsurfaceMesh").release();
+		}
+		return nullptr;
 	}
 	else
 	{
 		std::vector<float> layer_thickness;
 		for (unsigned i=0; i<nLayers; ++i)
 			layer_thickness.push_back(this->_edits[i]->text().toFloat());
-		new_mesh = mapper.createStaticLayers(*_msh, layer_thickness);
+		INFO("Mesh construction time: %d ms.", myTimer0.elapsed());
+		return mapper.createStaticLayers(*_msh, layer_thickness);
 	}
-	INFO("Mesh construction time: %d ms.", myTimer0.elapsed());
-
-	return new_mesh;
 }
 
 MeshLib::Mesh* MeshLayerEditDialog::createTetMesh()
@@ -231,7 +232,7 @@ MeshLib::Mesh* MeshLayerEditDialog::createTetMesh()
 		return nullptr;
 
 	const unsigned nLayers = _layerEdit->text().toInt();
-	MeshLib::Mesh* tg_mesh (nullptr);
+	MeshLib::Mesh* tg_mesh(nullptr);
 	QTime myTimer0;
 	myTimer0.start();
 
@@ -244,7 +245,7 @@ MeshLib::Mesh* MeshLayerEditDialog::createTetMesh()
 			raster_paths.push_back(this->_edits[i]->text().toStdString());
 		LayeredVolume lv;
 		if (lv.createLayers(*_msh, raster_paths, minimum_thickness))
-			tg_mesh = lv.getMesh("SubsurfaceMesh");
+			tg_mesh = lv.getMesh("SubsurfaceMesh").release();
 
 		if (tg_mesh)
 		{
diff --git a/Applications/Utils/MeshEdit/appendLinesAlongPolyline.cpp b/Applications/Utils/MeshEdit/appendLinesAlongPolyline.cpp
index 7f5ff9589bff37aaf087e8a9bcc1f2bed4542a45..a330bb5ee28fc46df2cd3e8ab1c32bca25c95344 100644
--- a/Applications/Utils/MeshEdit/appendLinesAlongPolyline.cpp
+++ b/Applications/Utils/MeshEdit/appendLinesAlongPolyline.cpp
@@ -84,16 +84,15 @@ int main (int argc, char* argv[])
 	INFO("Mesh read: %d nodes, %d elements.", mesh->getNNodes(), mesh->getNElements());
 
 	// add line elements
-	MeshLib::Mesh* new_mesh = MeshGeoToolsLib::appendLinesAlongPolylines(*mesh, *ply_vec);
+	std::unique_ptr<MeshLib::Mesh> new_mesh =
+	    MeshGeoToolsLib::appendLinesAlongPolylines(*mesh, *ply_vec);
 	INFO("Mesh created: %d nodes, %d elements.", new_mesh->getNNodes(), new_mesh->getNElements());
 
 	// write into a file
 	FileIO::Legacy::MeshIO meshIO;
-	meshIO.setMesh(new_mesh);
+	meshIO.setMesh(new_mesh.get());
 	meshIO.writeToFile(mesh_out.getValue());
 
-	delete new_mesh;
-
 	delete custom_format;
 	delete logog_cout;
 	LOGOG_SHUTDOWN();
diff --git a/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp b/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp
index 3509e1281dd390a68b45bd1b910e740a81e2dd9a..8d75b1302ad7be80e648f3009a59ae1ffd2b9232 100644
--- a/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp
+++ b/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp
@@ -25,16 +25,25 @@
 
 namespace MeshGeoToolsLib
 {
-
-MeshLib::Mesh* appendLinesAlongPolylines(const MeshLib::Mesh &mesh, const GeoLib::PolylineVec &ply_vec)
+std::unique_ptr<MeshLib::Mesh> appendLinesAlongPolylines(
+    const MeshLib::Mesh& mesh, const GeoLib::PolylineVec& ply_vec)
 {
-
 	// copy existing nodes and elements
 	std::vector<MeshLib::Node*> vec_new_nodes = MeshLib::copyNodeVector(mesh.getNodes());
 	std::vector<MeshLib::Element*> vec_new_eles = MeshLib::copyElementVector(mesh.getElements(), vec_new_nodes);
-	unsigned max_matID = 0;
-	for (auto e : vec_new_eles)
-		max_matID = std::max(max_matID, e->getValue());
+
+	std::vector<int> new_mat_ids;
+	{
+		auto mat_ids = mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+		if (mat_ids) {
+			new_mat_ids.reserve((*mat_ids).size());
+			std::copy((*mat_ids).cbegin(), (*mat_ids).cend(),
+				std::back_inserter(new_mat_ids));
+		}
+	}
+	int max_matID(0);
+	if (!new_mat_ids.empty())
+		max_matID = *(std::max_element(new_mat_ids.cbegin(), new_mat_ids.cend()));
 
 	const std::size_t n_ply (ply_vec.size());
 	// for each polyline
@@ -57,13 +66,25 @@ MeshLib::Mesh* appendLinesAlongPolylines(const MeshLib::Mesh &mesh, const GeoLib
 			std::array<MeshLib::Node*, 2> element_nodes;
 			element_nodes[0] = vec_new_nodes[vec_nodes_on_ply[i]];
 			element_nodes[1] = vec_new_nodes[vec_nodes_on_ply[i+1]];
-			vec_new_eles.push_back(new MeshLib::Line(element_nodes, max_matID+k+1));
+			vec_new_eles.push_back(
+				new MeshLib::Line(element_nodes, vec_new_eles.size()));
+			new_mat_ids.push_back(max_matID+k+1);
 		}
 	}
 
 	// generate a mesh
-	const std::string new_mesh_name = mesh.getName() + "_with_lines";
-	return new MeshLib::Mesh(new_mesh_name, vec_new_nodes, vec_new_eles);
+	const std::string name = mesh.getName() + "_with_lines";
+	std::unique_ptr<MeshLib::Mesh> new_mesh(
+		new MeshLib::Mesh(name, vec_new_nodes, vec_new_eles));
+	auto opt_mat_pv = new_mesh->getProperties().createNewPropertyVector<int>(
+		"MaterialIDs", MeshLib::MeshItemType::Cell);
+	if (opt_mat_pv) {
+		auto & mat_pv = *opt_mat_pv;
+		mat_pv.reserve(new_mat_ids.size());
+		std::copy(new_mat_ids.cbegin(), new_mat_ids.cend(),
+			std::back_inserter(mat_pv));
+	}
+	return new_mesh;
 }
 
 } // MeshGeoToolsLib
diff --git a/MeshGeoToolsLib/AppendLinesAlongPolyline.h b/MeshGeoToolsLib/AppendLinesAlongPolyline.h
index 7fd4159b92607907a2d2a7913456f28e6f050ea7..5213cb558f2e2e901f23ac5fe069a397768abff5 100644
--- a/MeshGeoToolsLib/AppendLinesAlongPolyline.h
+++ b/MeshGeoToolsLib/AppendLinesAlongPolyline.h
@@ -9,6 +9,8 @@
 #ifndef APPENDLINESALONGPOLYLINES_H_
 #define APPENDLINESALONGPOLYLINES_H_
 
+#include <memory>
+
 namespace GeoLib
 {
 class Polyline;
@@ -37,8 +39,8 @@ namespace MeshGeoToolsLib
  * @param ply_vec   polyline vector whose nodes are used to create line elements
  * @return a new mesh which is copied from a given mesh and additionally includes line elements
  */
-MeshLib::Mesh* appendLinesAlongPolylines(const MeshLib::Mesh &mesh, const GeoLib::PolylineVec &ply_vec);
-
+std::unique_ptr<MeshLib::Mesh> appendLinesAlongPolylines(
+    const MeshLib::Mesh& mesh, const GeoLib::PolylineVec& ply_vec);
 }
 
 #endif
diff --git a/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp b/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp
index 37809ccd550b58c905988dd7dc24546340e8d243..b8a37a71b056d34538c98c92389f6943e7a45b66 100644
--- a/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp
+++ b/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp
@@ -29,7 +29,7 @@ BoundaryElementsAtPoint::BoundaryElementsAtPoint(
 	std::array<MeshLib::Node*, 1> const nodes = {{
 	    const_cast<MeshLib::Node*>(_mesh.getNode(node_ids[0]))}};
 
-	_boundary_elements.push_back(new MeshLib::Point{nodes, 0, node_ids[0]});
+	_boundary_elements.push_back(new MeshLib::Point{nodes, node_ids[0]});
 }
 
 BoundaryElementsAtPoint::~BoundaryElementsAtPoint()
diff --git a/MeshLib/ElementStatus.cpp b/MeshLib/ElementStatus.cpp
index 47c67671f396e4fa1c80c74307fe4c87ad5dd0f4..a825faf7d57a6872e2038eea39623a8c7c1bcc7c 100644
--- a/MeshLib/ElementStatus.cpp
+++ b/MeshLib/ElementStatus.cpp
@@ -34,14 +34,19 @@ ElementStatus::ElementStatus(Mesh const* const mesh,
                              std::vector<unsigned> const& vec_inactive_matIDs)
     : ElementStatus(mesh, !vec_inactive_matIDs.empty())
 {
-	const std::size_t nElems (_mesh->getNElements());
-	for (auto material_id : vec_inactive_matIDs) {
-		for (auto e : _mesh->getElements())
-			if (e->getValue() == material_id)
-				this->setElementStatus(e->getID(), false);
+	auto materialIds = mesh->getProperties().getPropertyVector<int>("MaterialIDs");
+	if (materialIds) {
+		for (auto material_id : vec_inactive_matIDs) {
+			for (auto e : _mesh->getElements()) {
+				if ((*materialIds)[e->getID()] == material_id) {
+					setElementStatus(e->getID(), false);
+				}
+			}
+		}
 	}
 
-	_vec_active_eles.reserve(this->getNActiveElements());
+	_vec_active_eles.reserve(getNActiveElements());
+	const std::size_t nElems (_mesh->getNElements());
 	for (std::size_t i=0; i<nElems; ++i)
 		if (_element_status[i])
 			_vec_active_eles.push_back(const_cast<MeshLib::Element*>(_mesh->getElement(i)));
diff --git a/MeshLib/Elements/TemplateElement-impl.h b/MeshLib/Elements/TemplateElement-impl.h
index cba91a51bb326b423673efeec0c369ee5ed735e3..6667db0c115f927db86beb301d235e293a676974 100644
--- a/MeshLib/Elements/TemplateElement-impl.h
+++ b/MeshLib/Elements/TemplateElement-impl.h
@@ -22,6 +22,16 @@ TemplateElement<ELEMENT_RULE>::TemplateElement(Node* nodes[n_all_nodes], unsigne
 	this->_content = ELEMENT_RULE::computeVolume(this->_nodes);
 }
 
+template <class ELEMENT_RULE>
+TemplateElement<ELEMENT_RULE>::TemplateElement(Node* nodes[n_all_nodes], std::size_t id)
+: Element(id)
+{
+	this->_nodes = nodes;
+	this->_neighbors = new Element*[getNNeighbors()];
+	std::fill(this->_neighbors, this->_neighbors + getNNeighbors(), nullptr);
+	this->_content = ELEMENT_RULE::computeVolume(this->_nodes);
+}
+
 template <class ELEMENT_RULE>
 TemplateElement<ELEMENT_RULE>::TemplateElement(std::array<Node*, n_all_nodes> const& nodes, unsigned value, std::size_t id)
 : Element(value, id)
@@ -33,9 +43,20 @@ TemplateElement<ELEMENT_RULE>::TemplateElement(std::array<Node*, n_all_nodes> co
 	this->_content = ELEMENT_RULE::computeVolume(this->_nodes);
 }
 
+template <class ELEMENT_RULE>
+TemplateElement<ELEMENT_RULE>::TemplateElement(std::array<Node*, n_all_nodes> const& nodes, std::size_t id)
+: Element(id)
+{
+	this->_nodes = new Node*[n_all_nodes];
+	std::copy(nodes.begin(), nodes.end(), this->_nodes);
+	this->_neighbors = new Element*[getNNeighbors()];
+	std::fill(this->_neighbors, this->_neighbors + getNNeighbors(), nullptr);
+	this->_content = ELEMENT_RULE::computeVolume(this->_nodes);
+}
+
 template <class ELEMENT_RULE>
 TemplateElement<ELEMENT_RULE>::TemplateElement(const TemplateElement &e)
-: Element(e.getValue(), e.getID())
+: Element(e.getID())
 {
 	this->_nodes = new Node*[n_all_nodes];
 	for (unsigned i=0; i<n_all_nodes; i++)
diff --git a/MeshLib/Elements/TemplateElement.h b/MeshLib/Elements/TemplateElement.h
index 9fe5b5aaf144e5f3b74e827db94569eb2021a609..bb50401f5da49ac0716dbc06705bea48aaeadeac 100644
--- a/MeshLib/Elements/TemplateElement.h
+++ b/MeshLib/Elements/TemplateElement.h
@@ -48,7 +48,15 @@ public:
 	 * @param value  element value, e.g. material ID
 	 * @param id     element id
 	 */
-	TemplateElement(Node* nodes[n_all_nodes], unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
+	TemplateElement(Node* nodes[n_all_nodes], unsigned value, std::size_t id);
+
+	/**
+	 * Constructor with an array of mesh nodes.
+	 *
+	 * @param nodes  an array of pointers of mesh nodes which form this element
+	 * @param id     element id
+	 */
+	TemplateElement(Node* nodes[n_all_nodes], std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/**
 	 * Constructor with an array of mesh nodes
@@ -57,7 +65,15 @@ public:
 	 * @param value  element value, e.g. material ID
 	 * @param id     element id
 	 */
-	TemplateElement(std::array<Node*, n_all_nodes> const& nodes, unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
+	TemplateElement(std::array<Node*, n_all_nodes> const& nodes, unsigned value, std::size_t id);
+
+	/**
+	 * Constructor with an array of mesh nodes
+	 *
+	 * @param nodes  an array of pointers of mesh nodes which form this element
+	 * @param id     element id
+	 */
+	TemplateElement(std::array<Node*, n_all_nodes> const& nodes, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Copy constructor
 	TemplateElement(const TemplateElement &e);
diff --git a/MeshLib/MeshEditing/DuplicateMeshComponents.cpp b/MeshLib/MeshEditing/DuplicateMeshComponents.cpp
index 61f67419dd4878e3744a1c70705c8fd8152e0624..d358c6d6d96e7afb005010f14131b1bd0fae16f9 100644
--- a/MeshLib/MeshEditing/DuplicateMeshComponents.cpp
+++ b/MeshLib/MeshEditing/DuplicateMeshComponents.cpp
@@ -67,7 +67,7 @@ MeshLib::Element* copyElement(MeshLib::Element const*const element, const std::v
 	MeshLib::Node** new_nodes = new MeshLib::Node*[element->getNBaseNodes()];
 	for (unsigned i=0; i<element->getNBaseNodes(); ++i)
 		new_nodes[i] = nodes[element->getNode(i)->getID()];
-	return new E(new_nodes, element->getValue());
+	return new E(new_nodes);
 }
 
 } // namespace MeshLib
diff --git a/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp b/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp
index b0778b23891a4985c2021a59392e231481ad4e62..9fa0f45b899a885e69fd6b1d7729983d9d66ccc5 100644
--- a/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp
+++ b/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp
@@ -14,6 +14,7 @@
 
 #include <vector>
 #include <fstream>
+#include <boost/optional.hpp>
 
 #include "Mesh2MeshPropertyInterpolation.h"
 
@@ -73,9 +74,20 @@ void Mesh2MeshPropertyInterpolation::interpolatePropertiesForMesh(Mesh *dest_mes
 
 	GeoLib::Grid<MeshLib::Node> src_grid(src_nodes.begin(), src_nodes.end(), 64);
 
+	auto materialIds = dest_mesh->getProperties().getPropertyVector<int>("MaterialIDs");
+	if (!materialIds)
+		materialIds = dest_mesh->getProperties().createNewPropertyVector<int>(
+			"MaterialIDs", MeshLib::MeshItemType::Cell, 1);
+	if (!materialIds)
+	{
+		ERR("Could not create PropertyVector for MaterialIDs in Mesh.");
+		return;
+	}
+
 	std::vector<MeshLib::Element*> const& dest_elements(dest_mesh->getElements());
 	const std::size_t n_dest_elements(dest_elements.size());
-	for (std::size_t k(0); k<n_dest_elements; k++) {
+	for (std::size_t k(0); k<n_dest_elements; k++)
+	{
 		// compute axis aligned bounding box around the current element
 		const GeoLib::AABB<MeshLib::Node> elem_aabb(dest_elements[k]->getNodes(), dest_elements[k]->getNodes()+dest_elements[k]->getNBaseNodes());
 
@@ -101,7 +113,7 @@ void Mesh2MeshPropertyInterpolation::interpolatePropertiesForMesh(Mesh *dest_mes
 		}
 
 		dest_properties[k] /= cnt;
-		dest_elements[k]->setValue(k);
+		materialIds->push_back(k);
 
 		if (cnt == 0) {
 			std::string base_name("DebugData/Element-");
@@ -138,14 +150,17 @@ void Mesh2MeshPropertyInterpolation::interpolatePropertiesForMesh(Mesh *dest_mes
 
 void Mesh2MeshPropertyInterpolation::interpolateElementPropertiesToNodeProperties(std::vector<double> &interpolated_node_properties) const
 {
+	auto materialIds = _src_mesh->getProperties().getPropertyVector<int>("MaterialIDs");
+	if (!materialIds)
+		return;
+
 	std::vector<MeshLib::Node*> const& src_nodes(_src_mesh->getNodes());
 	const std::size_t n_src_nodes(src_nodes.size());
-
 	for (std::size_t k(0); k<n_src_nodes; k++) {
 		const std::size_t n_con_elems (src_nodes[k]->getNElements());
-		interpolated_node_properties[k] = (*_src_properties)[(src_nodes[k]->getElement(0))->getValue()];
+		interpolated_node_properties[k] = (*_src_properties)[(*materialIds)[src_nodes[k]->getElement(0)->getID()]];
 		for (std::size_t j(1); j<n_con_elems; j++) {
-			interpolated_node_properties[k] += (*_src_properties)[(src_nodes[k]->getElement(j))->getValue()];
+			interpolated_node_properties[k] += (*_src_properties)[(*materialIds)[src_nodes[k]->getElement(j)->getID()]];
 		}
 		interpolated_node_properties[k] /= n_con_elems;
 	}
diff --git a/MeshLib/MeshGenerators/ConvertRasterToMesh.cpp b/MeshLib/MeshGenerators/ConvertRasterToMesh.cpp
index 15dba953f235289a1c730b0228a57f4d2976d18f..6e75a32411506ca9a10f38b83763dd114de74e32 100644
--- a/MeshLib/MeshGenerators/ConvertRasterToMesh.cpp
+++ b/MeshLib/MeshGenerators/ConvertRasterToMesh.cpp
@@ -112,6 +112,7 @@ MeshLib::Mesh* ConvertRasterToMesh::constructMesh(const double* pix_vals, const
 		}
 	}
 
+	std::vector<int> mat_ids;
 	// set mesh elements
 	for (std::size_t i = 0; i < _raster.getNRows(); i++) {
 		for (std::size_t j = 0; j < _raster.getNCols(); j++) {
@@ -133,8 +134,12 @@ MeshLib::Mesh* ConvertRasterToMesh::constructMesh(const double* pix_vals, const
 					tri2_nodes[1] = nodes[node_idx_map[index + width + 1]];
 					tri2_nodes[2] = nodes[node_idx_map[index + width]];
 
-					elements.push_back(new MeshLib::Tri(tri1_nodes, mat)); // upper left triangle
-					elements.push_back(new MeshLib::Tri(tri2_nodes, mat)); // lower right triangle
+					// upper left triangle
+					elements.push_back(new MeshLib::Tri(tri1_nodes));
+					mat_ids.push_back(mat);
+					// lower right triangle
+					elements.push_back(new MeshLib::Tri(tri2_nodes));
+					mat_ids.push_back(mat);
 				}
 				if (_elem_type == MeshElemType::QUAD) {
 					MeshLib::Node** quad_nodes = new MeshLib::Node*[4];
@@ -142,14 +147,23 @@ MeshLib::Mesh* ConvertRasterToMesh::constructMesh(const double* pix_vals, const
 					quad_nodes[1] = nodes[node_idx_map[index + 1]];
 					quad_nodes[2] = nodes[node_idx_map[index + width + 1]];
 					quad_nodes[3] = nodes[node_idx_map[index + width]];
-					elements.push_back(new MeshLib::Quad(quad_nodes, mat));
+					elements.push_back(new MeshLib::Quad(quad_nodes));
+					mat_ids.push_back(mat);
 				}
 			}
 		}
 	}
 	delete [] node_idx_map;
 
-	return new MeshLib::Mesh("RasterDataMesh", nodes, elements); // the name is only a temp-name, the name given in the dialog is set later
+	// the name is only a temp-name, the name given in the dialog is set later
+	MeshLib::Properties properties;
+	boost::optional<MeshLib::PropertyVector<int>&> materials =
+	    properties.createNewPropertyVector<int>("MaterialIDs",
+	                                            MeshLib::MeshItemType::Cell);
+	assert(materials != boost::none);
+	materials->reserve(mat_ids.size());
+	std::copy(mat_ids.cbegin(), mat_ids.cend(), std::back_inserter(*materials));
+	return new MeshLib::Mesh("RasterDataMesh", nodes, elements, properties);
 }
 
 
diff --git a/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp b/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp
index 00f74050169e58962c83577e57bad51c3544dc9b..dbecc531732b0ff6b6834fc89cf967e55305c9ba 100644
--- a/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp
+++ b/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp
@@ -17,6 +17,8 @@
 #include <vector>
 #include <fstream>
 
+#include <logog/include/logog.hpp>
+
 #include "FileIO/AsciiRasterInterface.h"
 
 #include "GeoLib/Raster.h"
@@ -24,6 +26,8 @@
 #include "MeshLib/Mesh.h"
 #include "MeshLib/Node.h"
 #include "MeshLib/Elements/Element.h"
+#include "MeshLib/PropertyVector.h"
+#include "MeshLib/Properties.h"
 #include "MeshLib/MeshQuality/MeshValidation.h"
 #include "MeshLib/MeshSearch/NodeSearch.h"
 #include "MeshLib/MeshEditing/RemoveMeshComponents.h"
@@ -51,19 +55,35 @@ bool LayeredMeshGenerator::createLayers(MeshLib::Mesh const& mesh,
     return result;
 }
 
-MeshLib::Mesh* LayeredMeshGenerator::getMesh(std::string const& mesh_name) const
+std::unique_ptr<MeshLib::Mesh>
+LayeredMeshGenerator::getMesh(std::string const& mesh_name) const
 {
-    if (_nodes.empty() || _elements.empty())
-        return nullptr;
-
-    MeshLib::Mesh* result (new MeshLib::Mesh(mesh_name, _nodes, _elements));
-    MeshLib::NodeSearch ns(*result);
-    if (ns.searchUnused() > 0) {
-        auto new_mesh = MeshLib::removeNodes(*result, ns.getSearchedNodeIDs(), mesh_name);
-        delete result;
-        return new_mesh;
-    }
-    return result;
+	if (_nodes.empty() || _elements.empty())
+		return nullptr;
+
+	MeshLib::Properties properties;
+	if (_materials.size() == _elements.size())
+	{
+		boost::optional<MeshLib::PropertyVector<int>&> materials =
+		    properties.createNewPropertyVector<int>(
+		        "MaterialIDs", MeshLib::MeshItemType::Cell);
+		assert(materials != boost::none);
+		materials->reserve(_materials.size());
+		std::copy(_materials.cbegin(),
+		          _materials.cend(),
+		          std::back_inserter(*materials));
+	}
+	else
+		WARN ("Skipping MaterialID information, number of entries does not match element number");
+
+	std::unique_ptr<MeshLib::Mesh> result(new MeshLib::Mesh(mesh_name, _nodes, _elements, properties));
+	MeshLib::NodeSearch ns(*result.get());
+	if (ns.searchUnused() > 0) {
+		std::unique_ptr<MeshLib::Mesh> new_mesh(MeshLib::removeNodes(
+			*result.get(), ns.getSearchedNodeIDs(), mesh_name));
+		return new_mesh;
+	}
+	return result;
 }
 
 double LayeredMeshGenerator::calcEpsilon(GeoLib::Raster const& low, GeoLib::Raster const& high)
diff --git a/MeshLib/MeshGenerators/LayeredMeshGenerator.h b/MeshLib/MeshGenerators/LayeredMeshGenerator.h
index ec832bc47462029d2bc1185a5c9c84bb7a8e7ecf..1b5a4692ecf431e04825fa68acf21a27c4edc164 100644
--- a/MeshLib/MeshGenerators/LayeredMeshGenerator.h
+++ b/MeshLib/MeshGenerators/LayeredMeshGenerator.h
@@ -15,9 +15,10 @@
 #ifndef LAYEREDMESHGENERATOR_H
 #define LAYEREDMESHGENERATOR_H
 
+#include <memory>
+#include <limits>
 #include <string>
 #include <vector>
-#include <limits>
 
 namespace GeoLib {
     class Raster;
@@ -62,7 +63,7 @@ public:
 	                                double noDataReplacementValue) = 0;
 
 	/// Returns a mesh of the subsurface representation
-	MeshLib::Mesh* getMesh(std::string const& mesh_name) const;
+	std::unique_ptr<MeshLib::Mesh> getMesh(std::string const& mesh_name) const;
 
 protected:
 	LayeredMeshGenerator();
@@ -97,6 +98,7 @@ protected:
 
 	double _elevation_epsilon;
 	double _minimum_thickness;
+	std::vector<int> _materials;
 	std::vector<MeshLib::Node*> _nodes;
 	std::vector<MeshLib::Element*> _elements;
 };
diff --git a/MeshLib/MeshGenerators/LayeredVolume.cpp b/MeshLib/MeshGenerators/LayeredVolume.cpp
index 5bb71e3eee7a5f6f6f000557015ccd798a9e39cb..9e57d441717dc965afa6002aef852972530b5370 100644
--- a/MeshLib/MeshGenerators/LayeredVolume.cpp
+++ b/MeshLib/MeshGenerators/LayeredVolume.cpp
@@ -20,6 +20,7 @@
 
 #include "MeshLib/Elements/Tri.h"
 #include "MeshLib/Elements/Quad.h"
+#include "MeshLib/Properties.h"
 #include "MeshLib/MeshEditing/DuplicateMeshComponents.h"
 #include "MeshLib/MeshEditing/RemoveMeshComponents.h"
 #include "MeshLib/MeshGenerators/MeshLayerMapper.h"
@@ -41,24 +42,28 @@ bool LayeredVolume::createRasterLayers(const MeshLib::Mesh &mesh,
 	// remove line elements, only tri + quad remain
 	MeshLib::ElementSearch ex(mesh);
 	ex.searchByElementType(MeshLib::MeshElemType::LINE);
-	MeshLib::Mesh* top (removeElements(mesh, ex.getSearchedElementIDs(), "MeshLayer"));
+	std::unique_ptr<MeshLib::Mesh> top(
+	    removeElements(mesh, ex.getSearchedElementIDs(), "MeshLayer"));
 	if (top==nullptr)
-		top = new MeshLib::Mesh(mesh);
+		top.reset(new MeshLib::Mesh(mesh));
 
-	if (!MeshLib::MeshLayerMapper::layerMapping(*top, *rasters.back(), noDataReplacementValue))
+	if (!MeshLib::MeshLayerMapper::layerMapping(
+	        *top, *rasters.back(), noDataReplacementValue))
 		return false;
 
-	MeshLib::Mesh* bottom (new MeshLib::Mesh(*top));
+	std::unique_ptr<MeshLib::Mesh> bottom(new MeshLib::Mesh(*top));
 	if (!MeshLib::MeshLayerMapper::layerMapping(*bottom, *rasters[0], 0))
-	{
-		delete top;
 		return false;
-	}
 
 	this->_minimum_thickness = minimum_thickness;
 	_nodes = MeshLib::copyNodeVector(bottom->getNodes());
 	_elements = MeshLib::copyElementVector(bottom->getElements(), _nodes);
-	delete bottom;
+	if (!_materials.empty())
+	{
+		ERR("The materials vector is not empty.");
+		return false;
+	}
+	_materials.resize(_elements.size(), 0);
 
 	// map each layer and attach to subsurface mesh
 	const std::size_t nRasters (rasters.size());
@@ -68,7 +73,7 @@ bool LayeredVolume::createRasterLayers(const MeshLib::Mesh &mesh,
 	// close boundaries between layers
 	this->addLayerBoundaries(*top, nRasters);
 	this->removeCongruentElements(nRasters, top->getNElements());
-	delete top;
+
 	return true;
 }
 
@@ -90,7 +95,8 @@ void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer
 			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));
+			_elements.push_back(new MeshLib::Tri(tri_nodes));
+			_materials.push_back(layer_id);
 		}
 		else if (elem->getGeomType() == MeshLib::MeshElemType::QUAD)
 		{
@@ -98,7 +104,8 @@ void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer
 			                                             _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));
+			_elements.push_back(new MeshLib::Quad(quad_nodes));
+			_materials.push_back(layer_id);
 		}
 	}
 }
@@ -124,12 +131,14 @@ void LayeredVolume::addLayerBoundaries(const MeshLib::Mesh &layer, std::size_t n
 					if (MathLib::Vector3(*n1, *n2).getLength() > std::numeric_limits<double>::epsilon())
 					{
 						const std::array<MeshLib::Node*,3> tri_nodes = {{ n0, n2, n1 }};
-						_elements.push_back(new MeshLib::Tri(tri_nodes, nLayers+1+j));
+						_elements.push_back(new MeshLib::Tri(tri_nodes));
+						_materials.push_back(nLayers+j);
 					}
 					if (MathLib::Vector3(*n0, *n3).getLength() > std::numeric_limits<double>::epsilon())
 					{
 						const std::array<MeshLib::Node*,3> tri_nodes = {{ n0, n3, n2 }};
-						_elements.push_back(new MeshLib::Tri(tri_nodes, nLayers+1+j));
+						_elements.push_back(new MeshLib::Tri(tri_nodes));
+						_materials.push_back(nLayers+j);
 					}
 				}
 	}
@@ -158,15 +167,20 @@ void LayeredVolume::removeCongruentElements(std::size_t nLayers, std::size_t nEl
 			if (count == nElemNodes)
 			{
 				delete _elements[upper_offset+j];
+				// mark element and material entries for deletion
 				_elements[upper_offset+j] = nullptr;
+				_materials[upper_offset+j] = -1;
 			}
 			else
 			{
 				MeshLib::Node attr = high->getCenterOfGravity();
-				_attribute_points.push_back(MeshLib::Node(attr[0], attr[1], (attr[2] + low->getCenterOfGravity()[2])/2.0, low->getValue()));
+				_attribute_points.push_back(MeshLib::Node(attr[0], attr[1], (attr[2] + low->getCenterOfGravity()[2])/2.0, _materials[lower_offset+j]));
 			}
 		}
 	}
+	// delete marked entries
 	auto elem_vec_end = std::remove(_elements.begin(), _elements.end(), nullptr);
 	_elements.erase(elem_vec_end, _elements.end());
+	auto mat_vec_end = std::remove(_materials.begin(), _materials.end(), -1);
+	_materials.erase(mat_vec_end, _materials.end());
 }
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.cpp b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
index 3c3c5a5debd3fb49d8f574c77e346d5462fb8f76..c28a9769154ecf9c29e0f3cf9bd42b9bb5d3fc2d 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.cpp
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
@@ -29,6 +29,7 @@
 #include "MeshLib/Elements/Pyramid.h"
 #include "MeshLib/Elements/Prism.h"
 #include "MeshLib/MeshSurfaceExtraction.h"
+#include "MeshLib/Properties.h"
 
 namespace MeshLib
 {
@@ -60,6 +61,17 @@ MeshLib::Mesh* MeshLayerMapper::createStaticLayers(MeshLib::Mesh const& mesh, st
 	std::vector<MeshLib::Node*> new_nodes(nNodes + (nLayers * nNodes));
 	std::vector<MeshLib::Element*> new_elems;
 	new_elems.reserve(nElems * nLayers);
+	MeshLib::Properties properties;
+	boost::optional<PropertyVector<int>&> materials =
+	    properties.createNewPropertyVector<int>("MaterialIDs",
+	                                            MeshLib::MeshItemType::Cell);
+	if (!materials)
+	{
+		ERR("Could not create PropertyVector object \"MaterialIDs\".");
+		return nullptr;
+	}
+
+	materials->reserve(nElems * nLayers);
 	double z_offset (0.0);
 
 	for (unsigned layer_id = 0; layer_id <= nLayers; ++layer_id)
@@ -95,13 +107,14 @@ MeshLib::Mesh* MeshLayerMapper::createStaticLayers(MeshLib::Mesh const& mesh, st
 			}
 			// extrude triangles to prism
 			if (sfc_elem->getGeomType() == MeshLib::MeshElemType::TRIANGLE)
-				new_elems.push_back (new MeshLib::Prism(e_nodes, mat_id));
+				new_elems.push_back (new MeshLib::Prism(e_nodes));
 			// extrude quads to hexes
 			else if (sfc_elem->getGeomType() == MeshLib::MeshElemType::QUAD)
-				new_elems.push_back (new MeshLib::Hex(e_nodes, mat_id));
+				new_elems.push_back (new MeshLib::Hex(e_nodes));
+			materials->push_back(mat_id);
 		}
 	}
-	return new MeshLib::Mesh(mesh_name, new_nodes, new_elems);
+	return new MeshLib::Mesh(mesh_name, new_nodes, new_elems, properties);
 }
 
 bool MeshLayerMapper::createRasterLayers(
@@ -137,6 +150,7 @@ bool MeshLayerMapper::createRasterLayers(
 		[](MeshLib::Element const* elem)
 			{ return (elem->getGeomType() == MeshLib::MeshElemType::TRIANGLE);}));
 	_elements.reserve(nElems * (nLayers-1));
+	_materials.reserve(nElems *  (nLayers-1));
 
 	// add bottom layer
 	std::vector<MeshLib::Node*> const& nodes = bottom->getNodes();
@@ -145,8 +159,8 @@ bool MeshLayerMapper::createRasterLayers(
 	delete bottom;
 
 	// add the other layers
-	for (std::size_t i=1; i<nLayers; ++i)
-		addLayerToMesh(*top, i, *rasters[i]);
+	for (std::size_t i=0; i<nLayers-1; ++i)
+		addLayerToMesh(*top, i, *rasters[i+1]);
 
 	delete top;
 	return true;
@@ -163,7 +177,7 @@ 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;
+    int const last_layer_node_offset = layer_id * nNodes;
 
     // add nodes for new layer
     for (std::size_t i=0; i<nNodes; ++i)
@@ -192,7 +206,8 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
         switch (node_counter)
         {
         case 6:
-            _elements.push_back(new MeshLib::Prism(new_elem_nodes, layer_id));
+            _elements.push_back(new MeshLib::Prism(new_elem_nodes));
+            _materials.push_back(layer_id);
             break;
         case 5:
             std::array<MeshLib::Node*, 5> pyramid_nodes;
@@ -201,12 +216,14 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
             pyramid_nodes[2] = new_elem_nodes[pyramid_base[missing_idx][2]];
             pyramid_nodes[3] = new_elem_nodes[pyramid_base[missing_idx][3]];
             pyramid_nodes[4] = new_elem_nodes[missing_idx];
-            _elements.push_back(new MeshLib::Pyramid(pyramid_nodes, layer_id));
+            _elements.push_back(new MeshLib::Pyramid(pyramid_nodes));
+            _materials.push_back(layer_id);
             break;
         case 4:
             std::array<MeshLib::Node*, 4> tet_nodes;
             std::copy(new_elem_nodes.begin(), new_elem_nodes.begin() + node_counter, tet_nodes.begin());
-            _elements.push_back(new MeshLib::Tet(tet_nodes, layer_id));
+            _elements.push_back(new MeshLib::Tet(tet_nodes));
+            _materials.push_back(layer_id);
             break;
         default:
             continue;
diff --git a/MeshLib/MeshGenerators/VtkMeshConverter.cpp b/MeshLib/MeshGenerators/VtkMeshConverter.cpp
index 484f161efe7e72d9c2aa7ff2b05f9101edb39a7e..a6608337e7498a155cab8ea9e9a5bb71ac8421ff 100644
--- a/MeshLib/MeshGenerators/VtkMeshConverter.cpp
+++ b/MeshLib/MeshGenerators/VtkMeshConverter.cpp
@@ -41,12 +41,12 @@ namespace detail
 {
 template <class T_ELEMENT>
 MeshLib::Element* createElementWithSameNodeOrder(const std::vector<MeshLib::Node*> &nodes,
-		vtkIdList* const node_ids, unsigned material)
+		vtkIdList* const node_ids)
 {
 	MeshLib::Node** ele_nodes = new MeshLib::Node*[T_ELEMENT::n_all_nodes];
 	for (unsigned k(0); k<T_ELEMENT::n_all_nodes; k++)
 		ele_nodes[k] = nodes[node_ids->GetId(k)];
-	return new T_ELEMENT(ele_nodes, material);
+	return new T_ELEMENT(ele_nodes);
 }
 }
 
@@ -267,9 +267,11 @@ MeshLib::Mesh* VtkMeshConverter::constructMesh(const double* pixVal,
 
 	boost::optional< MeshLib::PropertyVector<int>& > materials =
 		properties.createNewPropertyVector<int>("MaterialIDs", MeshLib::MeshItemType::Cell, 1);
-	materials->insert(materials->end(), elements.size(), 0);
+	assert(materials != boost::none);
+	materials->resize(elements.size(), 0);
 
-	return new MeshLib::Mesh("RasterDataMesh", nodes, elements, properties); // the name is only a temp-name, the name given in the dialog is set later
+	// the name is only a temp-name, the name given in the dialog is set later
+	return new MeshLib::Mesh("RasterDataMesh", nodes, elements, properties);
 }
 
 MeshLib::Mesh* VtkMeshConverter::convertUnstructuredGrid(vtkUnstructuredGrid* grid, std::string const& mesh_name)
@@ -290,28 +292,25 @@ MeshLib::Mesh* VtkMeshConverter::convertUnstructuredGrid(vtkUnstructuredGrid* gr
 	// set mesh elements
 	const std::size_t nElems = grid->GetNumberOfCells();
 	std::vector<MeshLib::Element*> elements(nElems);
-	vtkDataArray* scalars = grid->GetCellData()->GetScalars("MaterialIDs");
 	auto node_ids = vtkSmartPointer<vtkIdList>::New();
 	for (std::size_t i = 0; i < nElems; i++)
 	{
 		MeshLib::Element* elem;
 		grid->GetCellPoints(i, node_ids);
 
-		const unsigned material = (scalars) ? static_cast<int>(scalars->GetComponent(i,0)) : 0;
-
 		int cell_type = grid->GetCellType(i);
 		switch (cell_type)
 		{
 		case VTK_LINE: {
-			elem = detail::createElementWithSameNodeOrder<MeshLib::Line>(nodes, node_ids, material);
+			elem = detail::createElementWithSameNodeOrder<MeshLib::Line>(nodes, node_ids);
 			break;
 		}
 		case VTK_TRIANGLE: {
-			elem = detail::createElementWithSameNodeOrder<MeshLib::Tri>(nodes, node_ids, material);
+			elem = detail::createElementWithSameNodeOrder<MeshLib::Tri>(nodes, node_ids);
 			break;
 		}
 		case VTK_QUAD: {
-			elem = detail::createElementWithSameNodeOrder<MeshLib::Quad>(nodes, node_ids, material);
+			elem = detail::createElementWithSameNodeOrder<MeshLib::Quad>(nodes, node_ids);
 			break;
 		}
 		case VTK_PIXEL: {
@@ -320,15 +319,15 @@ MeshLib::Mesh* VtkMeshConverter::convertUnstructuredGrid(vtkUnstructuredGrid* gr
 			quad_nodes[1] = nodes[node_ids->GetId(1)];
 			quad_nodes[2] = nodes[node_ids->GetId(3)];
 			quad_nodes[3] = nodes[node_ids->GetId(2)];
-			elem = new MeshLib::Quad(quad_nodes, material);
+			elem = new MeshLib::Quad(quad_nodes);
 			break;
 		}
 		case VTK_TETRA: {
-			elem = detail::createElementWithSameNodeOrder<MeshLib::Tet>(nodes, node_ids, material);
+			elem = detail::createElementWithSameNodeOrder<MeshLib::Tet>(nodes, node_ids);
 			break;
 		}
 		case VTK_HEXAHEDRON: {
-			elem = detail::createElementWithSameNodeOrder<MeshLib::Hex>(nodes, node_ids, material);
+			elem = detail::createElementWithSameNodeOrder<MeshLib::Hex>(nodes, node_ids);
 			break;
 		}
 		case VTK_VOXEL: {
@@ -341,11 +340,11 @@ MeshLib::Mesh* VtkMeshConverter::convertUnstructuredGrid(vtkUnstructuredGrid* gr
 			voxel_nodes[5] = nodes[node_ids->GetId(5)];
 			voxel_nodes[6] = nodes[node_ids->GetId(7)];
 			voxel_nodes[7] = nodes[node_ids->GetId(6)];
-			elem = new MeshLib::Hex(voxel_nodes, material);
+			elem = new MeshLib::Hex(voxel_nodes);
 			break;
 		}
 		case VTK_PYRAMID: {
-			elem = detail::createElementWithSameNodeOrder<MeshLib::Pyramid>(nodes, node_ids, material);
+			elem = detail::createElementWithSameNodeOrder<MeshLib::Pyramid>(nodes, node_ids);
 			break;
 		}
 		case VTK_WEDGE: {
@@ -355,35 +354,35 @@ MeshLib::Mesh* VtkMeshConverter::convertUnstructuredGrid(vtkUnstructuredGrid* gr
 				prism_nodes[i] = nodes[node_ids->GetId(i+3)];
 				prism_nodes[i+3] = nodes[node_ids->GetId(i)];
 			}
-			elem = new MeshLib::Prism(prism_nodes, material);
+			elem = new MeshLib::Prism(prism_nodes);
 			break;
 		}
 		case VTK_QUADRATIC_EDGE: {
-			elem = detail::createElementWithSameNodeOrder<MeshLib::Line3>(nodes, node_ids, material);
+			elem = detail::createElementWithSameNodeOrder<MeshLib::Line3>(nodes, node_ids);
 			break;
 		}
 		case VTK_QUADRATIC_TRIANGLE: {
-			elem = detail::createElementWithSameNodeOrder<MeshLib::Tri6>(nodes, node_ids, material);
+			elem = detail::createElementWithSameNodeOrder<MeshLib::Tri6>(nodes, node_ids);
 			break;
 		}
 		case VTK_QUADRATIC_QUAD: {
-			elem = detail::createElementWithSameNodeOrder<MeshLib::Quad8>(nodes, node_ids, material);
+			elem = detail::createElementWithSameNodeOrder<MeshLib::Quad8>(nodes, node_ids);
 			break;
 		}
 		case VTK_BIQUADRATIC_QUAD: {
-			elem = detail::createElementWithSameNodeOrder<MeshLib::Quad9>(nodes, node_ids, material);
+			elem = detail::createElementWithSameNodeOrder<MeshLib::Quad9>(nodes, node_ids);
 			break;
 		}
 		case VTK_QUADRATIC_TETRA: {
-			elem = detail::createElementWithSameNodeOrder<MeshLib::Tet10>(nodes, node_ids, material);
+			elem = detail::createElementWithSameNodeOrder<MeshLib::Tet10>(nodes, node_ids);
 			break;
 		}
 		case VTK_QUADRATIC_HEXAHEDRON: {
-			elem = detail::createElementWithSameNodeOrder<MeshLib::Hex20>(nodes, node_ids, material);
+			elem = detail::createElementWithSameNodeOrder<MeshLib::Hex20>(nodes, node_ids);
 			break;
 		}
 		case VTK_QUADRATIC_PYRAMID: {
-			elem = detail::createElementWithSameNodeOrder<MeshLib::Pyramid13>(nodes, node_ids, material);
+			elem = detail::createElementWithSameNodeOrder<MeshLib::Pyramid13>(nodes, node_ids);
 			break;
 		}
 		case VTK_QUADRATIC_WEDGE: {
@@ -400,7 +399,7 @@ MeshLib::Mesh* VtkMeshConverter::convertUnstructuredGrid(vtkUnstructuredGrid* gr
 			prism_nodes[11] = nodes[node_ids->GetId(13)];
 			for (unsigned i=0; i<3; ++i)
 				prism_nodes[12+i] = nodes[node_ids->GetId(11-i)];
-			elem = new MeshLib::Prism15(prism_nodes, material);
+			elem = new MeshLib::Prism15(prism_nodes);
 			break;
 		}
 		default:
diff --git a/MeshLib/MeshInformation.cpp b/MeshLib/MeshInformation.cpp
index cfbf7ae8039a55b4d3ed025ef6817942cd2880d9..cca73112f4162af060d1f1eaab2bc1913cd8f08d 100644
--- a/MeshLib/MeshInformation.cpp
+++ b/MeshLib/MeshInformation.cpp
@@ -19,15 +19,22 @@
 namespace MeshLib
 {
 
-const std::pair<unsigned, unsigned> MeshInformation::getValueBounds(const MeshLib::Mesh &mesh)
+const std::pair<int, int> MeshInformation::getValueBounds(const MeshLib::Mesh &mesh)
 {
-	const std::vector<MeshLib::Element*> &elements (mesh.getElements());
-	const auto minmax = std::minmax_element(elements.cbegin(), elements.cend(),
-        [](MeshLib::Element const*const a, MeshLib::Element const*const b)
-            {
-                return a->getValue() < b->getValue();
-        });
-	return std::make_pair<unsigned, unsigned>((*minmax.first)->getValue(), (*minmax.second)->getValue());
+	boost::optional<MeshLib::PropertyVector<int> const&> materialIds
+		= mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+	if (!materialIds) {
+		INFO("Mesh does not contain a property \"MaterialIDs\".");
+		return {std::numeric_limits<int>::max(),
+		        std::numeric_limits<int>::max()};
+	}
+	if (materialIds->empty()) {
+		INFO("Mesh does not contain values for the property \"MaterialIDs\".");
+		return {std::numeric_limits<int>::max(),
+		        std::numeric_limits<int>::max()};
+	}
+	auto mat_bounds = std::minmax_element(materialIds->cbegin(), materialIds->cend());
+	return {*(mat_bounds.first), *(mat_bounds.second)};
 }
 
 const GeoLib::AABB<MeshLib::Node> MeshInformation::getBoundingBox(const MeshLib::Mesh &mesh)
diff --git a/MeshLib/MeshInformation.h b/MeshLib/MeshInformation.h
index adec7ad3daa8c91eec0064e3bd71cea603641d40..fb542e4b6b97c9fbab8609958d4f88fd35c5c3fa 100644
--- a/MeshLib/MeshInformation.h
+++ b/MeshLib/MeshInformation.h
@@ -30,7 +30,7 @@ class MeshInformation
 {
 public:
 	/// Returns the smallest and largest MaterialID of the mesh.
-	static const std::pair<unsigned, unsigned> getValueBounds(const MeshLib::Mesh &mesh);
+	static const std::pair<int, int> getValueBounds(const MeshLib::Mesh &mesh);
 
 	/// Returns the bounding box of the mesh.
 	static const GeoLib::AABB<MeshLib::Node> getBoundingBox(const MeshLib::Mesh &mesh);
diff --git a/MeshLib/convertMeshToGeo.cpp b/MeshLib/convertMeshToGeo.cpp
index 1f0b91948058f7dff52f0ad5774f9c9316b3e229..a5dd7c16aaaa272e7d3c845779c8ec6a86e0b24e 100644
--- a/MeshLib/convertMeshToGeo.cpp
+++ b/MeshLib/convertMeshToGeo.cpp
@@ -51,7 +51,7 @@ bool convertMeshToGeo(const MeshLib::Mesh &mesh, GeoLib::GEOObjects &geo_objects
 	const std::vector<std::size_t> id_map (geo_objects.getPointVecObj(mesh_name)->getIDMap());
 
 	// elements to surface triangles conversion
-	const std::pair<unsigned, unsigned> bounds (MeshInformation::getValueBounds(mesh));
+	const auto bounds (MeshInformation::getValueBounds(mesh));
 	const unsigned nMatGroups(bounds.second-bounds.first+1);
 	auto sfcs = std::unique_ptr<std::vector<GeoLib::Surface*>>(new std::vector<GeoLib::Surface*>);
 	sfcs->reserve(nMatGroups);
@@ -62,15 +62,19 @@ bool convertMeshToGeo(const MeshLib::Mesh &mesh, GeoLib::GEOObjects &geo_objects
 	const std::vector<MeshLib::Element*> &elements = mesh.getElements();
 	const std::size_t nElems (mesh.getNElements());
 
+	auto materialIds = mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+	auto const materialIdExist = bool(materialIds);
+
 	for (unsigned i=0; i<nElems; ++i)
 	{
+		auto surfaceId = materialIdExist ? (*materialIds)[i]-bounds.first : 0;
 		MeshLib::Element* e (elements[i]);
 		if (e->getGeomType() == MeshElemType::TRIANGLE)
-			(*sfcs)[e->getValue()-bounds.first]->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(1)], id_map[e->getNodeIndex(2)]);
+			(*sfcs)[surfaceId]->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(1)], id_map[e->getNodeIndex(2)]);
 		if (e->getGeomType() == MeshElemType::QUAD)
 		{
-			(*sfcs)[e->getValue()-bounds.first]->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(1)], id_map[e->getNodeIndex(2)]);
-			(*sfcs)[e->getValue()-bounds.first]->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(2)], id_map[e->getNodeIndex(3)]);
+			(*sfcs)[surfaceId]->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(1)], id_map[e->getNodeIndex(2)]);
+			(*sfcs)[surfaceId]->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(2)], id_map[e->getNodeIndex(3)]);
 		}
 		// all other element types are ignored (i.e. lines)
 	}
@@ -95,7 +99,7 @@ MeshLib::Mesh* convertSurfaceToMesh(const GeoLib::Surface &sfc, const std::strin
 		MeshLib::Node** tri_nodes = new MeshLib::Node*[3];
 		for (unsigned j=0; j<3; j++)
 			tri_nodes[j] = new MeshLib::Node(tri->getPoint(j)->getCoords(), nodeId++);
-		elements.push_back(new MeshLib::Tri(tri_nodes, 0, i));
+		elements.push_back(new MeshLib::Tri(tri_nodes, i));
 		for (unsigned j=0; j<3; j++)
 			nodes.push_back(tri_nodes[j]);
 	}
diff --git a/Tests/MeshLib/TestElementStatus.cpp b/Tests/MeshLib/TestElementStatus.cpp
index c121267ac35b865f78a3d34a72f4489b915f1f06..06dbfc0cdc29d2cbdccd835ca42d6af0d5dd90f2 100644
--- a/Tests/MeshLib/TestElementStatus.cpp
+++ b/Tests/MeshLib/TestElementStatus.cpp
@@ -27,18 +27,25 @@ TEST(MeshLib, ElementStatus)
 	const unsigned width (100);
 	const unsigned elements_per_side (20);
 	auto const mesh = std::unique_ptr<MeshLib::Mesh>{
-        MeshLib::MeshGenerator::generateRegularQuadMesh(width, elements_per_side)};
+		MeshLib::MeshGenerator::generateRegularQuadMesh(width, elements_per_side)};
+
+	boost::optional<MeshLib::PropertyVector<int> &> material_id_properties(
+		mesh->getProperties().createNewPropertyVector<int>("MaterialIDs",
+		                                                   MeshLib::MeshItemType::Cell)
+	);
+	ASSERT_TRUE(bool(material_id_properties));
+	(*material_id_properties).resize(mesh->getNElements());
 
 	const std::vector<MeshLib::Element*> elements (mesh->getElements());
 
 	for (unsigned i=0; i<elements_per_side; ++i)
 	{
 		for (unsigned j=0; j<elements_per_side; ++j)
-			elements[i*elements_per_side + j]->setValue(i);
+			(*material_id_properties)[elements[i*elements_per_side + j]->getID()] = i;
 	}
 
 	{
-		// all elements and ndoes active
+		// all elements and nodes active
 		MeshLib::ElementStatus status(mesh.get());
 		ASSERT_EQ (elements.size(), status.getNActiveElements());
 		ASSERT_EQ (mesh->getNNodes(), status.getNActiveNodes());
diff --git a/Tests/MeshLib/TestVtkMeshConverter.cpp b/Tests/MeshLib/TestVtkMeshConverter.cpp
index 90c843c5bdc3f8ea64a4b288210994d0b8a28e49..72127e59b1066ab6c8c9ec616295f0974849f866 100644
--- a/Tests/MeshLib/TestVtkMeshConverter.cpp
+++ b/Tests/MeshLib/TestVtkMeshConverter.cpp
@@ -19,6 +19,7 @@
 #include <vtkUnsignedIntArray.h>
 #include <vtkUnstructuredGrid.h>
 
+#include "MeshLib/Elements/Element.h"
 #include "MeshLib/Mesh.h"
 #include "MeshLib/MeshGenerators/VtkMeshConverter.h"
 
@@ -96,6 +97,7 @@ TEST_F(TestVtkMeshConverter, Conversion)
 		MeshLib::VtkMeshConverter::convertUnstructuredGrid(vtu));
 	ASSERT_EQ(mesh->getNNodes(), vtu->GetNumberOfPoints());
 	ASSERT_EQ(mesh->getNElements(), vtu->GetNumberOfCells());
+	ASSERT_EQ(mesh->getElement(0)->getCellType(), MeshLib::CellType::HEX8);
 
 	auto meshProperties = mesh->getProperties();
 
diff --git a/Tests/NumLib/TestFunctionInterpolation.cpp b/Tests/NumLib/TestFunctionInterpolation.cpp
index dd6b3a353acb7e57df75e39be0e366ab88bdfcac..1792946cde20fd5cf9bec3f779d077cc1946b923 100644
--- a/Tests/NumLib/TestFunctionInterpolation.cpp
+++ b/Tests/NumLib/TestFunctionInterpolation.cpp
@@ -58,9 +58,7 @@ TEST(NumLibFunctionInterpolationTest, Linear1DElement)
     auto pt_a = MeshLib::Node({ 0.0, 0.0, 0.0 }, 0);
     auto pt_b = MeshLib::Node({ 1.0, 0.0, 0.0 }, 0);
 
-    auto const element = MeshLib::Line(std::array<MeshLib::Node*, 2>{
-                                           &pt_a, &pt_b
-                                       }, 0u, 0u);
+    auto const element = MeshLib::Line(std::array<MeshLib::Node*, 2>{&pt_a, &pt_b});
 
     // set up shape function
     FemType finite_element(*static_cast<const ShapeFunction::MeshElement*>(&element));