diff --git a/MeshLib/ElementStatus.cpp b/MeshLib/ElementStatus.cpp
index 1b46559feb9b251130cd83aa747c6b75b1384250..556e9c35b160906b9cfc86e919b63f6322a66304 100644
--- a/MeshLib/ElementStatus.cpp
+++ b/MeshLib/ElementStatus.cpp
@@ -1,5 +1,5 @@
 /**
- * \file
+ * \file   ElementStatus.cpp
  * \author Karsten Rink
  * \date   2012-12-18
  * \brief  Implementation of the ElementStatus class.
@@ -14,13 +14,108 @@
 
 #include "ElementStatus.h"
 
+#include "Mesh.h"
+#include "Node.h"
+#include "Elements/Element.h"
+
 namespace MeshLib {
 
 ElementStatus::ElementStatus(Mesh const*const mesh)
-: _mesh(mesh), _status(mesh->getNElements(), true)
+: _mesh(mesh), _element_status(mesh->getNElements(), true)
+{
+	const std::vector<MeshLib::Node*> &nodes (_mesh->getNodes());
+	for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
+		_active_nodes.push_back((*node)->getNElements());
+}
+
+
+std::vector<std::size_t> ElementStatus::getActiveElements() const
+{
+	std::vector<std::size_t> active_elements;
+	active_elements.reserve(this->getNActiveElements());
+	const std::size_t nElems (_mesh->getNElements());
+	for (std::size_t i=0; i<nElems; ++i)
+		if (_element_status[i])
+			active_elements.push_back(i);
+	return active_elements;
+}
+
+std::vector<std::size_t> ElementStatus::getActiveNodes() const 
+{
+	std::vector<std::size_t> active_nodes;
+	active_nodes.reserve(this->getNActiveNodes());
+	const std::size_t nNodes (_mesh->getNNodes());
+	for (std::size_t i=0; i<nNodes; ++i)
+		if (_active_nodes[i]>0)
+			active_nodes.push_back(i);
+	return active_nodes;
+}
+
+std::vector<std::size_t> ElementStatus::getActiveElementsAtNode(std::size_t node_id) const
+{
+	const std::size_t nActiveElements (_active_nodes[node_id]);
+	const std::vector<Element*> &elements (_mesh->getNode(node_id)->getElements());
+	std::vector<std::size_t> active_elements;
+	active_elements.reserve(nActiveElements);
+	for (auto elem = elements.cbegin(); elem != elements.cend(); ++elem)
+	{
+		if (active_elements.size() == nActiveElements)
+			return active_elements;
+		const std::size_t id ((*elem)->getID());
+		if (_element_status[id])
+			active_elements.push_back(id);
+	}
+	return active_elements;
+}
+
+std::size_t ElementStatus::getNActiveNodes() const 
+{
+	return _active_nodes.size() - std::count(_active_nodes.cbegin(), _active_nodes.cend(), 0);
+}
+
+std::size_t ElementStatus::getNActiveElements() const 
 {
+	return static_cast<std::size_t>(std::count(_element_status.cbegin(), _element_status.cend(), true));
 }
 
+void ElementStatus::setAll(bool status)
+{
+	std::fill(_element_status.begin(), _element_status.end(), status);
+
+	if (status)
+	{
+		const std::vector<MeshLib::Node*> &nodes (_mesh->getNodes());
+		const std::size_t nNodes (_mesh->getNNodes());
+		for (std::size_t i=0; i<nNodes; ++i)
+			_active_nodes[i] = nodes[i]->getNElements();
+	}
+	else
+		std::fill(_active_nodes.begin(), _active_nodes.end(), 0);
+}
+
+void ElementStatus::setElementStatus(std::size_t i, bool status)
+{
+	if (_element_status[i] != status)
+	{
+		const int change = (status) ? 1 : -1;
+		_element_status[i] = status;
+		const unsigned nElemNodes (_mesh->getElement(i)->getNNodes());
+		MeshLib::Node const*const*const nodes = _mesh->getElement(i)->getNodes();
+		for (unsigned i=0; i<nElemNodes; ++i)
+		{
+			assert(_active_nodes[i]<255); // if one node has >255 connected elements the data type is too small
+			_active_nodes[nodes[i]->getID()] += change;
+		}
+	}
+}
+
+void ElementStatus::setMaterialStatus(unsigned material_id, bool status)
+{
+	const std::size_t nElems (_mesh->getNElements());
+	for (std::size_t i=0; i<nElems; ++i)
+		if (_mesh->getElement(i)->getValue() == material_id)
+			this->setElementStatus(i, status);
+}
 
 }
 
diff --git a/MeshLib/ElementStatus.h b/MeshLib/ElementStatus.h
index 9fff2ef52a8393a0c48c3305ba2ab10fcb32c3b0..a838daa039e5c0fed99530aa1a7554a59b199848 100644
--- a/MeshLib/ElementStatus.h
+++ b/MeshLib/ElementStatus.h
@@ -1,5 +1,5 @@
 /**
- * \file
+ * \file   ElementStatus.h
  * \author Karsten Rink
  * \date   2012-12-18
  * \brief  Definition of the ElementStatus class.
@@ -15,35 +15,58 @@
 #ifndef ELEMENTSTATUS_H_
 #define ELEMENTSTATUS_H_
 
-#include "Mesh.h"
+#include <vector>
 
 namespace MeshLib {
+	class Mesh;
+	class Element;
 
 /**
- * An object to store which elements of the mesh are active and which are inactive
+ * Manages the active/inactive flags for mesh elements and their nodes
  */
 class ElementStatus
 {
 
 public:
-	/// Constructor using mesh
-	ElementStatus(Mesh const*const mesh);
+	/// Constructor
+	explicit ElementStatus(MeshLib::Mesh const*const mesh);
 
-	bool getElementStatus(unsigned i) { return _status[i]; };
+	/// Returns a vector of active element IDs
+	std::vector<std::size_t> getActiveElements() const;
 
-	bool isActive(unsigned i) { return _status[i]; };
+	/// Returns a vector of active node IDs
+	std::vector<std::size_t> getActiveNodes() const;
+	
+	/// Returns the status of element i
+	bool getElementStatus(std::size_t i) const { return _element_status[i]; }
 
-	void setActive(unsigned i) { _status[i] = true; };
+	/// Returns a vector of active elements connected to a node
+	std::vector<std::size_t> getActiveElementsAtNode(std::size_t node_id) const;
+	
+	/// Returns the total number of active nodes
+	std::size_t getNActiveNodes() const;
 
-	void setInactive(unsigned i) { _status[i] = false; };
+	/// Returns the total number of active elements
+	std::size_t getNActiveElements() const;
 
-	void setElementStatus(unsigned i, bool status) { _status[i] = status; };
+	/// Activates/Deactives all mesh elements
+	void setAll(bool status);
+
+	/// Sets the status of element i
+	void setElementStatus(std::size_t i, bool status);
+
+	/// Sets the status of material group i
+	void setMaterialStatus(unsigned material_id, bool status);
 
 	~ElementStatus() {};
 
 protected:
-	Mesh const*const _mesh;
-	std::vector<bool> _status;
+	/// The mesh for which the element status is administrated
+	MeshLib::Mesh const*const _mesh;
+	/// Element status for each mesh element (active/inactive = true/false)
+	std::vector<bool> _element_status;
+	/// Node status for each mesh node (value = number of active elements connected to node, 0 means inactive)
+	std::vector<unsigned char> _active_nodes;
 
 }; /* class */
 
diff --git a/MeshLib/Elements/Cell.cpp b/MeshLib/Elements/Cell.cpp
index 89314a2fa861095f65bc16f4bc19d2df718ef87d..075a5210740ae6b0e52c62efa8cbb8541d1b6dd4 100644
--- a/MeshLib/Elements/Cell.cpp
+++ b/MeshLib/Elements/Cell.cpp
@@ -21,8 +21,8 @@ Cell::Cell(Node** nodes, MeshElemType type, unsigned value)
 {
 }
 */
-Cell::Cell(unsigned value)
-	: Element(value), _volume(-1.0) // init with invalid value to detect errors
+Cell::Cell(unsigned value, std::size_t id)
+	: Element(value, id), _volume(-1.0) // init with invalid value to detect errors
 {
 }
 
diff --git a/MeshLib/Elements/Cell.h b/MeshLib/Elements/Cell.h
index e3983bfe6c4dd4a69d1a39636fafd2b566d93f9e..897558573f8b9739418604345d650aa1d6f13f44 100644
--- a/MeshLib/Elements/Cell.h
+++ b/MeshLib/Elements/Cell.h
@@ -54,7 +54,7 @@ protected:
 	Cell(Node** nodes, MeshElemType type, unsigned value = 0);
 */
 	/// Constructor for a generic mesh element without an array of mesh nodes.
-	Cell(unsigned value = 0);
+	Cell(unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	double _volume;
 
diff --git a/MeshLib/Elements/Edge.cpp b/MeshLib/Elements/Edge.cpp
index 1d1e2f31bd206986fc075a694bb4fe7a70bb494d..f818edcaa8cb911b1837d4576384c386171a7b49 100644
--- a/MeshLib/Elements/Edge.cpp
+++ b/MeshLib/Elements/Edge.cpp
@@ -20,8 +20,8 @@
 
 namespace MeshLib
 {
-Edge::Edge(unsigned value)
-    : Element(value), _length(-1.0) // init with invalid value to detect errors
+Edge::Edge(unsigned value, std::size_t id)
+    : Element(value, id), _length(-1.0) // init with invalid value to detect errors
 {
 }
 
diff --git a/MeshLib/Elements/Edge.h b/MeshLib/Elements/Edge.h
index c13d8d74637fbaf9b6ec30e2deb272569bdbb086..47e560fb45bd9f43dd05c4f40127149dce038dc6 100644
--- a/MeshLib/Elements/Edge.h
+++ b/MeshLib/Elements/Edge.h
@@ -85,7 +85,7 @@ protected:
     unsigned identifyFace(Node* [3]/*nodes[3]*/) const { return std::numeric_limits<unsigned>::max(); };
 
     /// Constructor for a generic mesh element without an array of mesh nodes.
-    Edge(unsigned value = 0);
+    Edge(unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
     double _length;
 
diff --git a/MeshLib/Elements/Element.cpp b/MeshLib/Elements/Element.cpp
index da7796bbefbfb300c275ecea44c85a696ed3beda..a70892710f8c18d78ffb026e057e1fc47378bb91 100644
--- a/MeshLib/Elements/Element.cpp
+++ b/MeshLib/Elements/Element.cpp
@@ -22,8 +22,8 @@
 
 namespace MeshLib {
 
-Element::Element(unsigned value)
-	: _nodes(nullptr), _value(value), _neighbors(nullptr)
+Element::Element(unsigned value, std::size_t id)
+	: _nodes(nullptr), _id(id), _value(value), _neighbors(nullptr)
 {
 }
 
diff --git a/MeshLib/Elements/Element.h b/MeshLib/Elements/Element.h
index 50017fcb044c7282e9e7b1940070df90778f8560..00cb03f8a3bcfa8e3d372f393bfec562f4c4f399 100644
--- a/MeshLib/Elements/Element.h
+++ b/MeshLib/Elements/Element.h
@@ -77,6 +77,9 @@ public:
 	/// Returns the i-th face of the element.
 	virtual const Element* getFace(unsigned i) const = 0;
 
+	/// Returns the ID of the element.
+	virtual std::size_t getID() const { return this->_id; }
+
 	/// Get the number of edges for this element.
 	virtual unsigned getNEdges() const = 0;
 
@@ -176,7 +179,7 @@ public:
 
 protected:
 	/// Constructor for a generic mesh element without an array of mesh nodes.
-	Element(unsigned value = 0);
+	Element(unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Return a specific edge node.
 	virtual Node* getEdgeNode(unsigned edge_id, unsigned node_id) const = 0;
@@ -184,8 +187,12 @@ protected:
 	/// Returns the ID of a face given an array of nodes.
 	virtual unsigned identifyFace(Node* nodes[3]) const = 0;
 
+	/// Sets the element ID.
+	virtual void setID(std::size_t id) { this->_id = id; }
+
 
 	Node** _nodes;
+	std::size_t _id;
 	/**
 	 * this is an index for external additional information like materials
 	 */
diff --git a/MeshLib/Elements/Face.cpp b/MeshLib/Elements/Face.cpp
index 3134989bc727bf227fe624addb5bc51631801e03..7915fdf68ac7c0cf6765d9410f82b6b35aac024a 100644
--- a/MeshLib/Elements/Face.cpp
+++ b/MeshLib/Elements/Face.cpp
@@ -29,8 +29,8 @@ Face::Face(Node** nodes, MeshElemType type, unsigned value)
 {
 }
 */
-Face::Face(unsigned value)
-	: Element(value), _area(-1.0) // init with invalid value to detect errors
+Face::Face(unsigned value, std::size_t id)
+	: Element(value, id), _area(-1.0) // init with invalid value to detect errors
 {
 }
 
diff --git a/MeshLib/Elements/Face.h b/MeshLib/Elements/Face.h
index 69c2b9126f286a0210ca481869f355bc444059a7..94729334f4d05622d2ab9f9959fcff05851a7618 100644
--- a/MeshLib/Elements/Face.h
+++ b/MeshLib/Elements/Face.h
@@ -78,7 +78,7 @@ protected:
 	Face(Node** nodes, MeshElemType type, unsigned value = 0);
 */
 	/// Constructor for a generic mesh element without an array of mesh nodes.
-	Face(unsigned value = 0);
+	Face(unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	double _area;
 
diff --git a/MeshLib/Elements/TemplateHex-impl.h b/MeshLib/Elements/TemplateHex-impl.h
index c31298d9f5290b98f2018e952cb6241ae75c9173..64bbb3c482df62c45d4884d66330f23945c77c66 100644
--- a/MeshLib/Elements/TemplateHex-impl.h
+++ b/MeshLib/Elements/TemplateHex-impl.h
@@ -51,8 +51,8 @@ const unsigned TemplateHex<NNODES,CELLHEXTYPE>::_edge_nodes[12][2] =
 };
 
 template <unsigned NNODES, CellType CELLHEXTYPE>
-TemplateHex<NNODES,CELLHEXTYPE>::TemplateHex(Node* nodes[NNODES], unsigned value)
-	: Cell(value)
+TemplateHex<NNODES,CELLHEXTYPE>::TemplateHex(Node* nodes[NNODES], unsigned value, std::size_t id)
+	: Cell(value, id)
 {
 	_nodes = nodes;
 
@@ -64,8 +64,8 @@ TemplateHex<NNODES,CELLHEXTYPE>::TemplateHex(Node* nodes[NNODES], unsigned value
 
 template<unsigned NNODES, CellType CELLHEXTYPE>
 TemplateHex<NNODES,CELLHEXTYPE>::TemplateHex(std::array<Node*, NNODES> const& nodes,
-                                             unsigned value)
-	: Cell(value)
+                                             unsigned value, std::size_t id)
+	: Cell(value, id)
 {
 	_nodes = new Node*[NNODES];
 	std::copy(nodes.begin(), nodes.end(), _nodes);
@@ -78,7 +78,7 @@ TemplateHex<NNODES,CELLHEXTYPE>::TemplateHex(std::array<Node*, NNODES> const& no
 
 template <unsigned NNODES, CellType CELLHEXTYPE>
 TemplateHex<NNODES,CELLHEXTYPE>::TemplateHex(const TemplateHex<NNODES,CELLHEXTYPE> &hex)
-	: Cell(hex.getValue())
+	: Cell(hex.getValue(), hex.getID())
 {
 	_nodes = new Node*[NNODES];
 	for (unsigned i=0; i<NNODES; i++)
diff --git a/MeshLib/Elements/TemplateHex.h b/MeshLib/Elements/TemplateHex.h
index 1adc68b9238220bbca79b51dedc581cce75f6cfb..ba3567a3b777a5ff89b9799d9e73998f49d6fe6b 100644
--- a/MeshLib/Elements/TemplateHex.h
+++ b/MeshLib/Elements/TemplateHex.h
@@ -51,10 +51,10 @@ class TemplateHex : public Cell
 {
 public:
 	/// Constructor with an array of mesh nodes.
-	TemplateHex(Node* nodes[NNODES], unsigned value = 0);
+	TemplateHex(Node* nodes[NNODES], unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Constructs a hex from array of Node pointers.
-	TemplateHex(std::array<Node*, NNODES> const& nodes, unsigned value = 0);
+	TemplateHex(std::array<Node*, NNODES> const& nodes, unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Copy constructor
 	TemplateHex(const TemplateHex &hex);
diff --git a/MeshLib/Elements/TemplateLine-impl.h b/MeshLib/Elements/TemplateLine-impl.h
index 328facc17d7a3ab92381eada463dd18e38f2b401..b3d098930c7102f8c766ee16f847e7606ff2042f 100644
--- a/MeshLib/Elements/TemplateLine-impl.h
+++ b/MeshLib/Elements/TemplateLine-impl.h
@@ -16,8 +16,8 @@ namespace MeshLib
 {
 template<unsigned NNODES, CellType CELLLINETYPE>
 TemplateLine<NNODES,CELLLINETYPE>::TemplateLine(std::array<Node*, NNODES> const& nodes,
-                                                unsigned value)
-	: Edge(value)
+                                                unsigned value, std::size_t id)
+	: Edge(value, id)
 {
 	_nodes = new Node*[NNODES];
 	std::copy(nodes.begin(), nodes.end(), _nodes);
@@ -26,8 +26,8 @@ TemplateLine<NNODES,CELLLINETYPE>::TemplateLine(std::array<Node*, NNODES> const&
 }
 
 template<unsigned NNODES, CellType CELLLINETYPE>
-TemplateLine<NNODES,CELLLINETYPE>::TemplateLine(Node* nodes[NNODES], unsigned value) :
-	Edge(value)
+TemplateLine<NNODES,CELLLINETYPE>::TemplateLine(Node* nodes[NNODES], unsigned value, std::size_t id) :
+	Edge(value, id)
 {
 	_nodes = nodes;
 	this->_length = this->computeVolume();
@@ -35,7 +35,7 @@ TemplateLine<NNODES,CELLLINETYPE>::TemplateLine(Node* nodes[NNODES], unsigned va
 
 template <unsigned NNODES, CellType CELLLINETYPE>
 TemplateLine<NNODES,CELLLINETYPE>::TemplateLine(const TemplateLine<NNODES,CELLLINETYPE> &line) :
-	Edge(line.getValue())
+	Edge(line.getValue(), line.getID())
 {
 	_nodes = new Node*[NNODES];
 	for (unsigned k(0); k<NNODES; k++)
diff --git a/MeshLib/Elements/TemplateLine.h b/MeshLib/Elements/TemplateLine.h
index ad07b20c2b668aef1cc15ee68571e6b12ba117df..dcdebc513edc46e40747717e9eb8eeca124601f3 100644
--- a/MeshLib/Elements/TemplateLine.h
+++ b/MeshLib/Elements/TemplateLine.h
@@ -40,10 +40,10 @@ class TemplateLine : public Edge
 {
 public:
 	/// Constructor with an array of mesh nodes.
-	TemplateLine(Node* nodes[NNODES], unsigned value = 0);
+	TemplateLine(Node* nodes[NNODES], unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Constructs a line from array of Node pointers.
-	TemplateLine(std::array<Node*, NNODES> const& nodes, unsigned value = 0);
+	TemplateLine(std::array<Node*, NNODES> const& nodes, unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Copy constructor
 	TemplateLine(const TemplateLine &line);
diff --git a/MeshLib/Elements/TemplatePrism-impl.h b/MeshLib/Elements/TemplatePrism-impl.h
index a81ce1ca08e5fb4ec8099e4a8d832a6f70e8598d..2decf3bcb5792b270c28e425ce8b47071a2d4a04 100644
--- a/MeshLib/Elements/TemplatePrism-impl.h
+++ b/MeshLib/Elements/TemplatePrism-impl.h
@@ -52,8 +52,8 @@ template <unsigned NNODES, CellType CELLPRISMTYPE>
 const unsigned TemplatePrism<NNODES,CELLPRISMTYPE>::_n_face_nodes[5] = { 3, 4, 4, 4, 3 };
 
 template <unsigned NNODES, CellType CELLPRISMTYPE>
-TemplatePrism<NNODES,CELLPRISMTYPE>::TemplatePrism(Node* nodes[NNODES], unsigned value)
-	: Cell(value)
+TemplatePrism<NNODES,CELLPRISMTYPE>::TemplatePrism(Node* nodes[NNODES], unsigned value, std::size_t id)
+	: Cell(value, id)
 {
 	_nodes = nodes;
 	_neighbors = new Element*[5];
@@ -63,8 +63,8 @@ TemplatePrism<NNODES,CELLPRISMTYPE>::TemplatePrism(Node* nodes[NNODES], unsigned
 
 template<unsigned NNODES, CellType CELLPRISMTYPE>
 TemplatePrism<NNODES,CELLPRISMTYPE>::TemplatePrism(std::array<Node*, NNODES> const& nodes,
-                                                   unsigned value)
-	: Cell(value)
+                                                   unsigned value, std::size_t id)
+	: Cell(value, id)
 {
 	_nodes = new Node*[NNODES];
 	std::copy(nodes.begin(), nodes.end(), _nodes);
@@ -77,7 +77,7 @@ TemplatePrism<NNODES,CELLPRISMTYPE>::TemplatePrism(std::array<Node*, NNODES> con
 
 template <unsigned NNODES, CellType CELLPRISMTYPE>
 TemplatePrism<NNODES,CELLPRISMTYPE>::TemplatePrism(const TemplatePrism<NNODES,CELLPRISMTYPE> &prism)
-	: Cell(prism.getValue())
+	: Cell(prism.getValue(), prism.getID())
 {
 	_nodes = new Node*[NNODES];
 	for (unsigned i=0; i<NNODES; i++)
diff --git a/MeshLib/Elements/TemplatePrism.h b/MeshLib/Elements/TemplatePrism.h
index 2faa97f38ea32ac23613650d7f2ab4676a2def47..98afdec004f1a087eeebe05eb8f42aa34426dd7c 100644
--- a/MeshLib/Elements/TemplatePrism.h
+++ b/MeshLib/Elements/TemplatePrism.h
@@ -49,10 +49,10 @@ class TemplatePrism : public Cell
 {
 public:
 	/// Constructor with an array of mesh nodes.
-	TemplatePrism(Node* nodes[NNODES], unsigned value = 0);
+	TemplatePrism(Node* nodes[NNODES], unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Constructs a prism from array of Node pointers.
-	TemplatePrism(std::array<Node*, NNODES> const& nodes, unsigned value = 0);
+	TemplatePrism(std::array<Node*, NNODES> const& nodes, unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Copy constructor
 	TemplatePrism(const TemplatePrism &prism);
diff --git a/MeshLib/Elements/TemplatePyramid-impl.h b/MeshLib/Elements/TemplatePyramid-impl.h
index 6dc3ed276910d04fec9680a3ab7fbe52aaf4d196..4df80819c0e3626950c924f46d78cda306274767 100644
--- a/MeshLib/Elements/TemplatePyramid-impl.h
+++ b/MeshLib/Elements/TemplatePyramid-impl.h
@@ -51,8 +51,8 @@ template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
 const unsigned TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::_n_face_nodes[5] = { 3, 3, 3, 3, 4 };
 
 template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::TemplatePyramid(Node* nodes[NNODES], unsigned value)
-	: Cell(value)
+TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::TemplatePyramid(Node* nodes[NNODES], unsigned value, std::size_t id)
+	: Cell(value, id)
 {
 	_nodes = nodes;
 
@@ -64,8 +64,8 @@ TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::TemplatePyramid(Node* nodes[NNODES], un
 
 template<unsigned NNODES, CellType CELLPYRAMIDTYPE>
 TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::TemplatePyramid(std::array<Node*, NNODES> const& nodes,
-                                                         unsigned value)
-	: Cell(value)
+                                                         unsigned value, std::size_t id)
+	: Cell(value, id)
 {
 	_nodes = new Node*[NNODES];
 	std::copy(nodes.begin(), nodes.end(), _nodes);
@@ -78,7 +78,7 @@ TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::TemplatePyramid(std::array<Node*, NNODE
 
 template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
 TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::TemplatePyramid(const TemplatePyramid<NNODES,CELLPYRAMIDTYPE> &pyramid)
-	: Cell(pyramid.getValue())
+	: Cell(pyramid.getValue(), pyramid.getID())
 {
 	_nodes = new Node*[NNODES];
 	for (unsigned i=0; i<NNODES; i++) {
diff --git a/MeshLib/Elements/TemplatePyramid.h b/MeshLib/Elements/TemplatePyramid.h
index 64f7a50c7518cb5e4f5ea1372a6b3a4ff82e1a7d..2630e3c65af5721b34eb271c94208eeb14d614fe 100644
--- a/MeshLib/Elements/TemplatePyramid.h
+++ b/MeshLib/Elements/TemplatePyramid.h
@@ -47,10 +47,10 @@ class TemplatePyramid : public Cell
 {
 public:
 	/// Constructor with an array of mesh nodes.
-	TemplatePyramid(Node* nodes[NNODES], unsigned value = 0);
+	TemplatePyramid(Node* nodes[NNODES], unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Constructs a pyramid from array of Node pointers.
-	TemplatePyramid(std::array<Node*, NNODES> const& nodes, unsigned value = 0);
+	TemplatePyramid(std::array<Node*, NNODES> const& nodes, unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Copy constructor
 	TemplatePyramid(const TemplatePyramid &pyramid);
diff --git a/MeshLib/Elements/TemplateQuad-impl.h b/MeshLib/Elements/TemplateQuad-impl.h
index 781ba81f9c12aaacae13b36f7eaf1e13c84d391d..ab017977b14584b388630e46023c70f1922e5ad7 100644
--- a/MeshLib/Elements/TemplateQuad-impl.h
+++ b/MeshLib/Elements/TemplateQuad-impl.h
@@ -25,8 +25,8 @@ namespace MeshLib
 {
 
 template <unsigned NNODES, CellType CELLQUADTYPE>
-TemplateQuad<NNODES,CELLQUADTYPE>::TemplateQuad(Node* nodes[NNODES], unsigned value)
-	: Face(value)
+TemplateQuad<NNODES,CELLQUADTYPE>::TemplateQuad(Node* nodes[NNODES], unsigned value, std::size_t id)
+	: Face(value, id)
 {
 	_nodes = nodes;
 
@@ -38,8 +38,8 @@ TemplateQuad<NNODES,CELLQUADTYPE>::TemplateQuad(Node* nodes[NNODES], unsigned va
 
 template<unsigned NNODES, CellType CELLQUADTYPE>
 TemplateQuad<NNODES,CELLQUADTYPE>::TemplateQuad(std::array<Node*, NNODES> const& nodes,
-                                                unsigned value)
-	: Face(value)
+                                                unsigned value, std::size_t id)
+	: Face(value, id)
 {
 	_nodes = new Node*[NNODES];
 	std::copy(nodes.begin(), nodes.end(), _nodes);
@@ -52,7 +52,7 @@ TemplateQuad<NNODES,CELLQUADTYPE>::TemplateQuad(std::array<Node*, NNODES> const&
 
 template <unsigned NNODES, CellType CELLQUADTYPE>
 TemplateQuad<NNODES,CELLQUADTYPE>::TemplateQuad(const TemplateQuad<NNODES,CELLQUADTYPE> &quad)
-	: Face(quad.getValue())
+	: Face(quad.getValue(), quad.getID())
 {
 	_nodes = new Node*[NNODES];
 	for (unsigned i=0; i<NNODES; i++) {
diff --git a/MeshLib/Elements/TemplateQuad.h b/MeshLib/Elements/TemplateQuad.h
index 869539a465ec82efa4208345d0218dd75d278bbd..909ecdce73598be27d47ffb42bfadf26b64a666d 100644
--- a/MeshLib/Elements/TemplateQuad.h
+++ b/MeshLib/Elements/TemplateQuad.h
@@ -47,10 +47,10 @@ public:
 	static const unsigned n_base_nodes;
 
 	/// Constructor with an array of mesh nodes.
-	TemplateQuad(Node* nodes[NNODES], unsigned value = 0);
+	TemplateQuad(Node* nodes[NNODES], unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Constructs an edge from array of Node pointers.
-	TemplateQuad(std::array<Node*, NNODES> const& nodes, unsigned value = 0);
+	TemplateQuad(std::array<Node*, NNODES> const& nodes, unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Constructs a quad from NNODES of Nodes initializing Face with
 	//  value = 0.
diff --git a/MeshLib/Elements/TemplateTet-impl.h b/MeshLib/Elements/TemplateTet-impl.h
index 5b738fee1996c09b17edd72c812cf1096e734587..49370381650325fddf5f445d4a1ce9d532a02098 100644
--- a/MeshLib/Elements/TemplateTet-impl.h
+++ b/MeshLib/Elements/TemplateTet-impl.h
@@ -42,8 +42,8 @@ const unsigned TemplateTet<NNODES,CELLTETTYPE>::_edge_nodes[6][2] =
 };
 
 template <unsigned NNODES, CellType CELLTETTYPE>
-TemplateTet<NNODES,CELLTETTYPE>::TemplateTet(Node* nodes[NNODES], unsigned value)
-	: Cell(value)
+TemplateTet<NNODES,CELLTETTYPE>::TemplateTet(Node* nodes[NNODES], unsigned value, std::size_t id)
+	: Cell(value, id)
 {
 	_nodes = nodes;
 
@@ -55,8 +55,8 @@ TemplateTet<NNODES,CELLTETTYPE>::TemplateTet(Node* nodes[NNODES], unsigned value
 
 template<unsigned NNODES, CellType CELLTETTYPE>
 TemplateTet<NNODES,CELLTETTYPE>::TemplateTet(std::array<Node*, NNODES> const& nodes,
-                                             unsigned value)
-	: Cell(value)
+                                             unsigned value, std::size_t id)
+	: Cell(value, id)
 {
 	_nodes = new Node*[NNODES];
 	std::copy(nodes.begin(), nodes.end(), _nodes);
@@ -69,7 +69,7 @@ TemplateTet<NNODES,CELLTETTYPE>::TemplateTet(std::array<Node*, NNODES> const& no
 
 template <unsigned NNODES, CellType CELLTETTYPE>
 TemplateTet<NNODES,CELLTETTYPE>::TemplateTet(const TemplateTet<NNODES,CELLTETTYPE> &tet)
-	: Cell(tet.getValue())
+	: Cell(tet.getValue(), tet.getID())
 {
 	_nodes = new Node*[NNODES];
 	for (unsigned i=0; i<NNODES; i++) {
diff --git a/MeshLib/Elements/TemplateTet.h b/MeshLib/Elements/TemplateTet.h
index b21a693b2a281f073f79254a434c9f4b6922d866..a1d3c32c5079ea379256f3aba5c5aa3c8df24291 100644
--- a/MeshLib/Elements/TemplateTet.h
+++ b/MeshLib/Elements/TemplateTet.h
@@ -46,10 +46,10 @@ class TemplateTet : public Cell
 {
 public:
 	/// Constructor with an array of mesh nodes.
-	TemplateTet(Node* nodes[NNODES], unsigned value = 0);
+	TemplateTet(Node* nodes[NNODES], unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Constructs a tetrahedron from array of Node pointers.
-	TemplateTet(std::array<Node*, NNODES> const& nodes, unsigned value = 0);
+	TemplateTet(std::array<Node*, NNODES> const& nodes, unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Copy constructor
 	TemplateTet(const TemplateTet &tet);
diff --git a/MeshLib/Elements/TemplateTri-impl.h b/MeshLib/Elements/TemplateTri-impl.h
index 870ef659aa3d62248ee7024fa5d65612dcadf32f..5bae658eba7b09ad57a20a11afff286abc1b0285 100644
--- a/MeshLib/Elements/TemplateTri-impl.h
+++ b/MeshLib/Elements/TemplateTri-impl.h
@@ -18,8 +18,8 @@
 namespace MeshLib {
 
 template <unsigned NNODES, CellType CELLTRITYPE>
-TemplateTri<NNODES,CELLTRITYPE>::TemplateTri(Node* nodes[NNODES], unsigned value) :
-	Face(value)
+TemplateTri<NNODES,CELLTRITYPE>::TemplateTri(Node* nodes[NNODES], unsigned value, std::size_t id) :
+	Face(value, id)
 {
 	_nodes = nodes;
 	_neighbors = new Element*[3];
@@ -29,8 +29,8 @@ TemplateTri<NNODES,CELLTRITYPE>::TemplateTri(Node* nodes[NNODES], unsigned value
 
 template<unsigned NNODES, CellType CELLTRITYPE>
 TemplateTri<NNODES,CELLTRITYPE>::TemplateTri(std::array<Node*, NNODES> const& nodes,
-                                             unsigned value)
-	: Face(value)
+                                             unsigned value, std::size_t id)
+	: Face(value, id)
 {
 	_nodes = new Node*[NNODES];
 	std::copy(nodes.begin(), nodes.end(), _nodes);
@@ -43,7 +43,7 @@ TemplateTri<NNODES,CELLTRITYPE>::TemplateTri(std::array<Node*, NNODES> const& no
 
 template <unsigned NNODES, CellType CELLTRITYPE>
 TemplateTri<NNODES,CELLTRITYPE>::TemplateTri(const TemplateTri<NNODES,CELLTRITYPE> &tri) :
-	Face(tri.getValue())
+	Face(tri.getValue(), tri.getID())
 {
 	_nodes = new Node*[NNODES];
 	for (unsigned i=0; i<NNODES; i++)
diff --git a/MeshLib/Elements/TemplateTri.h b/MeshLib/Elements/TemplateTri.h
index 499c551cd1c9e79e3118d2322cdfbe723b0aaa92..e0eee786f5bc1da979c5fb8ae3bc1a1261182710 100644
--- a/MeshLib/Elements/TemplateTri.h
+++ b/MeshLib/Elements/TemplateTri.h
@@ -48,10 +48,10 @@ class TemplateTri : public Face
 {
 public:
 	/// Constructor with an array of mesh nodes.
-	TemplateTri(Node* nodes[NNODES], unsigned value = 0);
+	TemplateTri(Node* nodes[NNODES], unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Constructs a triangle from array of Node pointers.
-	TemplateTri(std::array<Node*, NNODES> const& nodes, unsigned value = 0);
+	TemplateTri(std::array<Node*, NNODES> const& nodes, unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
 	/// Copy constructor
 	TemplateTri(const TemplateTri &tri);
diff --git a/MeshLib/Mesh.cpp b/MeshLib/Mesh.cpp
index 5b520e6ca1710272d8fd1358f03a97b511e64e56..7629cdc933246ec9d7979125421229673ccd492b 100644
--- a/MeshLib/Mesh.cpp
+++ b/MeshLib/Mesh.cpp
@@ -35,7 +35,8 @@ Mesh::Mesh(const std::string &name,
            const std::vector<Element*> &elements)
 	: _id(_counter_value), _mesh_dimension(0), _name(name), _nodes(nodes), _elements(elements)
 {
-	this->resetNodeIDs(); // reset node ids so they match the node position in the vector
+	this->resetNodeIDs();
+	this->resetElementIDs();
 	this->setDimension();
 	this->setElementsConnectedToNodes();
 	//this->setNodesConnectedByEdges();
@@ -108,6 +109,13 @@ void Mesh::resetNodeIDs()
 		_nodes[i]->setID(i);
 }
 
+void Mesh::resetElementIDs()
+{
+	const size_t nElements (this->_elements.size());
+	for (unsigned i=0; i<nElements; ++i)
+		_elements[i]->setID(i);
+}
+
 void Mesh::setDimension()
 {
 	const size_t nElements (_elements.size());
diff --git a/MeshLib/Mesh.h b/MeshLib/Mesh.h
index 3523a3deddc682c24e71dd0bac34ecd181bf7d33..03092d39306d639cf94b8953706533157fd152cc 100644
--- a/MeshLib/Mesh.h
+++ b/MeshLib/Mesh.h
@@ -84,6 +84,9 @@ public:
 	/// Get the element-vector for the mesh.
 	std::vector<Element*> const& getElements() const { return _elements; }
 
+	/// Resets the IDs of all mesh-elements to their position in the element vector
+	void resetElementIDs();
+	
 	/// Resets the IDs of all mesh-nodes to their position in the node vector
 	void resetNodeIDs();
 
diff --git a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
index 3dd163733569767075b6c6211745bf8fb5590e9f..5cdbb3ccc7d58d40b2a452e40609d9e6ccf772ba 100644
--- a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
+++ b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
@@ -11,6 +11,7 @@
 
 #include "FixedTimeStepping.h"
 
+#include <algorithm>
 #include <cmath>
 #include <numeric>
 #include <limits>
diff --git a/Tests/MeshLib/TestElementStatus.cpp b/Tests/MeshLib/TestElementStatus.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6b75ec19d1f54ecf860f9cda0cc5144103cdcccc
--- /dev/null
+++ b/Tests/MeshLib/TestElementStatus.cpp
@@ -0,0 +1,90 @@
+/**
+ * @file TestElementStatus.cpp
+ * @author Karsten Rink
+ * @date 2013-03-14
+ * @brief Tests for ElementStatus class
+ *
+ * @copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include "gtest/gtest.h"
+
+#include "Mesh.h"
+#include "Node.h"
+#include "Elements/Element.h"
+#include "ElementStatus.h"
+#include "MeshGenerators/MeshGenerator.h"
+
+
+TEST(MeshLib, ElementStatus)
+{
+	const unsigned width (100);
+	const unsigned elements_per_side (20);
+	MeshLib::Mesh* mesh (MeshLib::MeshGenerator::generateRegularQuadMesh(width, elements_per_side));
+	MeshLib::ElementStatus status(mesh);
+	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);
+	}
+
+	// all elements active
+	ASSERT_EQ (elements.size(), status.getNActiveElements());
+	// all nodes active
+	ASSERT_EQ (mesh->getNNodes(), status.getNActiveNodes());
+
+	// set material 1 to false
+	status.setMaterialStatus(1, false);
+	ASSERT_EQ (elements.size()-elements_per_side, status.getNActiveElements());
+
+	// set material 1 to false (again)
+	status.setMaterialStatus(1, false);
+	ASSERT_EQ (elements.size()-elements_per_side, status.getNActiveElements());
+
+	// set material 0 to false
+	status.setMaterialStatus(0, false);
+	ASSERT_EQ (elements.size()-(2*elements_per_side), status.getNActiveElements());
+
+	// active elements
+	std::vector<std::size_t> active_elements (status.getActiveElements());
+	ASSERT_EQ (active_elements.size(), status.getNActiveElements());
+
+	// active nodes
+	std::vector<std::size_t> active_nodes (status.getActiveNodes());
+	ASSERT_EQ (active_nodes.size(), status.getNActiveNodes());
+
+	// set element 1 to false (yet again)
+	status.setElementStatus(1, false);
+	status.getElementStatus(1);
+	ASSERT_EQ (elements.size()-(2*elements_per_side), status.getNActiveElements());
+	ASSERT_EQ (mesh->getNNodes()-(2*(elements_per_side+1)), status.getNActiveNodes());
+
+	// set element 1 to true
+	status.setElementStatus(1, true);
+	ASSERT_EQ (elements.size()-(2*elements_per_side)+1, status.getNActiveElements());
+	ASSERT_EQ (mesh->getNNodes()-(2*(elements_per_side+1))+4, status.getNActiveNodes());
+	ASSERT_EQ(status.getElementStatus(1), true);
+
+	std::vector<std::size_t> active_elements_at_node (status.getActiveElementsAtNode(2));
+	ASSERT_EQ(1, active_elements_at_node.size());
+	active_elements_at_node = status.getActiveElementsAtNode(22);
+	ASSERT_EQ(1, active_elements_at_node.size());
+	active_elements_at_node = status.getActiveElementsAtNode(44);
+	ASSERT_EQ(2, active_elements_at_node.size());
+	active_elements_at_node = status.getActiveElementsAtNode(102);
+	ASSERT_EQ(4, active_elements_at_node.size());
+
+	status.setAll(true);
+	ASSERT_EQ(elements.size(), status.getNActiveElements());
+	ASSERT_EQ(mesh->getNNodes(), status.getNActiveNodes());
+
+	status.setAll(false);
+	ASSERT_EQ(0, status.getNActiveElements());
+	ASSERT_EQ(0, status.getNActiveNodes());
+}