From dab073b7ce1f8f8bc08e8cb5af62fef02e9fecf9 Mon Sep 17 00:00:00 2001
From: Norihiro Watanabe <norihiro.watanabe@ufz.de>
Date: Mon, 26 Jan 2015 12:16:15 +0100
Subject: [PATCH] replace mesh element class with new implementation using
 TemplateElement class

---
 MeshLib/Elements/Cell.cpp               |   8 +-
 MeshLib/Elements/Cell.h                 |  22 +--
 MeshLib/Elements/CellRule.cpp           |  41 +++++
 MeshLib/Elements/CellRule.h             |  36 +++++
 MeshLib/Elements/Edge.h                 |  52 ++++++
 MeshLib/Elements/EdgeReturn.cpp         |  49 ++++++
 MeshLib/Elements/EdgeReturn.h           |  51 ++++++
 MeshLib/Elements/Element.cpp            |   2 +-
 MeshLib/Elements/Element.h              |  11 +-
 MeshLib/Elements/Face.cpp               |   7 +-
 MeshLib/Elements/Face.h                 |  29 +---
 MeshLib/Elements/FaceRule.cpp           |  17 ++
 MeshLib/Elements/FaceRule.h             |  39 +++++
 MeshLib/Elements/Hex.h                  |   6 +-
 MeshLib/Elements/HexRule8.cpp           | 119 ++++++++++++++
 MeshLib/Elements/HexRule8.h             | 110 +++++++++++++
 MeshLib/Elements/Line.h                 |   6 +-
 MeshLib/Elements/LineRule2.cpp          |  60 +++++++
 MeshLib/Elements/LineRule2.h            |  90 +++++++++++
 MeshLib/Elements/Prism.h                |   6 +-
 MeshLib/Elements/PrismRule6.cpp         | 123 +++++++++++++++
 MeshLib/Elements/PrismRule6.h           | 111 +++++++++++++
 MeshLib/Elements/Pyramid.h              |   6 +-
 MeshLib/Elements/PyramidRule5.cpp       | 120 ++++++++++++++
 MeshLib/Elements/PyramidRule5.h         | 110 +++++++++++++
 MeshLib/Elements/Quad.h                 |  15 +-
 MeshLib/Elements/QuadRule4.cpp          |  73 +++++++++
 MeshLib/Elements/QuadRule4.h            |  92 +++++++++++
 MeshLib/Elements/QuadRule8.cpp          |  24 +++
 MeshLib/Elements/QuadRule8.h            |  56 +++++++
 MeshLib/Elements/QuadRule9.cpp          |  16 ++
 MeshLib/Elements/QuadRule9.h            |  48 ++++++
 MeshLib/Elements/TemplateElement-impl.h |  71 +++++++++
 MeshLib/Elements/TemplateElement.h      | 155 ++++++++++++++++++
 MeshLib/Elements/TemplateHex-impl.h     | 195 -----------------------
 MeshLib/Elements/TemplateHex.h          | 156 ------------------
 MeshLib/Elements/TemplateLine-impl.h    |  85 ----------
 MeshLib/Elements/TemplateLine.h         | 188 ----------------------
 MeshLib/Elements/TemplatePrism-impl.h   | 200 ------------------------
 MeshLib/Elements/TemplatePrism.h        | 156 ------------------
 MeshLib/Elements/TemplatePyramid-impl.h | 200 ------------------------
 MeshLib/Elements/TemplatePyramid.h      | 154 ------------------
 MeshLib/Elements/TemplateQuad-impl.h    | 150 ------------------
 MeshLib/Elements/TemplateQuad.h         | 140 -----------------
 MeshLib/Elements/TemplateTet-impl.h     | 168 --------------------
 MeshLib/Elements/TemplateTet.h          | 155 ------------------
 MeshLib/Elements/TemplateTri-impl.h     | 119 --------------
 MeshLib/Elements/TemplateTri.h          | 154 ------------------
 MeshLib/Elements/Tet.h                  |   6 +-
 MeshLib/Elements/TetRule4.cpp           |  91 +++++++++++
 MeshLib/Elements/TetRule4.h             | 105 +++++++++++++
 MeshLib/Elements/Tri.h                  |   6 +-
 MeshLib/Elements/TriRule3.cpp           |  64 ++++++++
 MeshLib/Elements/TriRule3.h             |  93 +++++++++++
 54 files changed, 2062 insertions(+), 2304 deletions(-)
 create mode 100644 MeshLib/Elements/CellRule.cpp
 create mode 100644 MeshLib/Elements/CellRule.h
 create mode 100644 MeshLib/Elements/Edge.h
 create mode 100644 MeshLib/Elements/EdgeReturn.cpp
 create mode 100644 MeshLib/Elements/EdgeReturn.h
 create mode 100644 MeshLib/Elements/FaceRule.cpp
 create mode 100644 MeshLib/Elements/FaceRule.h
 create mode 100644 MeshLib/Elements/HexRule8.cpp
 create mode 100644 MeshLib/Elements/HexRule8.h
 create mode 100644 MeshLib/Elements/LineRule2.cpp
 create mode 100644 MeshLib/Elements/LineRule2.h
 create mode 100644 MeshLib/Elements/PrismRule6.cpp
 create mode 100644 MeshLib/Elements/PrismRule6.h
 create mode 100644 MeshLib/Elements/PyramidRule5.cpp
 create mode 100644 MeshLib/Elements/PyramidRule5.h
 create mode 100644 MeshLib/Elements/QuadRule4.cpp
 create mode 100644 MeshLib/Elements/QuadRule4.h
 create mode 100644 MeshLib/Elements/QuadRule8.cpp
 create mode 100644 MeshLib/Elements/QuadRule8.h
 create mode 100644 MeshLib/Elements/QuadRule9.cpp
 create mode 100644 MeshLib/Elements/QuadRule9.h
 create mode 100644 MeshLib/Elements/TemplateElement-impl.h
 create mode 100644 MeshLib/Elements/TemplateElement.h
 delete mode 100644 MeshLib/Elements/TemplateHex-impl.h
 delete mode 100644 MeshLib/Elements/TemplateHex.h
 delete mode 100644 MeshLib/Elements/TemplateLine-impl.h
 delete mode 100644 MeshLib/Elements/TemplateLine.h
 delete mode 100644 MeshLib/Elements/TemplatePrism-impl.h
 delete mode 100644 MeshLib/Elements/TemplatePrism.h
 delete mode 100644 MeshLib/Elements/TemplatePyramid-impl.h
 delete mode 100644 MeshLib/Elements/TemplatePyramid.h
 delete mode 100644 MeshLib/Elements/TemplateQuad-impl.h
 delete mode 100644 MeshLib/Elements/TemplateQuad.h
 delete mode 100644 MeshLib/Elements/TemplateTet-impl.h
 delete mode 100644 MeshLib/Elements/TemplateTet.h
 delete mode 100644 MeshLib/Elements/TemplateTri-impl.h
 delete mode 100644 MeshLib/Elements/TemplateTri.h
 create mode 100644 MeshLib/Elements/TetRule4.cpp
 create mode 100644 MeshLib/Elements/TetRule4.h
 create mode 100644 MeshLib/Elements/TriRule3.cpp
 create mode 100644 MeshLib/Elements/TriRule3.h

diff --git a/MeshLib/Elements/Cell.cpp b/MeshLib/Elements/Cell.cpp
index 7b7ae70b979..01253939f01 100644
--- a/MeshLib/Elements/Cell.cpp
+++ b/MeshLib/Elements/Cell.cpp
@@ -16,14 +16,10 @@
 
 #include "MathLib/Vector3.h"
 #include "MeshLib/Node.h"
+#include "MeshLib/Elements/Face.h"
 
 namespace MeshLib {
 
-#ifndef WIN32
-/// \todo Windows compiler does not accept this definition and issues a linking error.
-const unsigned Cell::dimension;
-#endif	// WIN32
-
 /*
 Cell::Cell(Node** nodes, MeshElemType type, unsigned value)
 	: Element(nodes, type, value)
@@ -31,7 +27,7 @@ Cell::Cell(Node** nodes, MeshElemType type, unsigned value)
 }
 */
 Cell::Cell(unsigned value, std::size_t id)
-	: Element(value, id), _volume(-1.0) // init with invalid value to detect errors
+	: Element(value, id)
 {
 }
 
diff --git a/MeshLib/Elements/Cell.h b/MeshLib/Elements/Cell.h
index dcc04f76bb4..8c06822cf83 100644
--- a/MeshLib/Elements/Cell.h
+++ b/MeshLib/Elements/Cell.h
@@ -16,7 +16,6 @@
 #define CELL_H_
 
 #include "Element.h"
-#include "Face.h"
 
 namespace MeshLib {
 
@@ -26,28 +25,12 @@ namespace MeshLib {
 class Cell : public Element
 {
 public:
-	/// Constant: Dimension of this mesh element
-	static const unsigned dimension = 3;
-
-	/// Returns the length, area or volume of a 1D, 2D or 3D element
-	double getContent() const { return _volume; }
-
-	/// Get dimension of the mesh element.
-	unsigned getDimension() const { return dimension; }
-
 	/// Get the volume of this 3d element.
-	virtual double getVolume() const { return _volume; }
+	virtual double getVolume() const { return getContent(); }
 
 	/// Destructor
 	virtual ~Cell();
 
-	/**
-	 * This method is pure virtual and is inherited from class @sa Element.
-	 * It has to be implemented in the derived classes of class Cell!
-	 * @return a copy of the object
-	 */
-	virtual Element* clone() const = 0;
-
 	/**
 	 * Checks if the node order of an element is correct by testing surface normals.
 	 * For 3D elements true is returned if the normals of all faces points away from the centre of 
@@ -65,9 +48,6 @@ protected:
 */
 	/// Constructor for a generic mesh element without an array of mesh nodes.
 	Cell(unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
-
-	double _volume;
-
 }; /* class */
 
 }
diff --git a/MeshLib/Elements/CellRule.cpp b/MeshLib/Elements/CellRule.cpp
new file mode 100644
index 00000000000..9fa8ea6ba90
--- /dev/null
+++ b/MeshLib/Elements/CellRule.cpp
@@ -0,0 +1,41 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CellRule.h"
+
+#include "logog/include/logog.hpp"
+
+#include "GeoLib/AnalyticalGeometry.h"
+
+#include "MeshLib/Node.h"
+#include "Face.h"
+
+namespace MeshLib
+{
+
+const unsigned CellRule::dimension;
+
+bool CellRule::testElementNodeOrder(const Element* e)
+{
+	const MathLib::Vector3 c (e->getCenterOfGravity());
+	const unsigned nFaces (e->getNFaces());
+	for (unsigned j=0; j<nFaces; ++j)
+	{
+		MeshLib::Face const*const face (dynamic_cast<const MeshLib::Face*>(e->getFace(j)));
+		const MeshLib::Node x (*(face->getNode(1)));
+		const MathLib::Vector3 cx (c, x);
+		const double s = MathLib::scalarProduct(face->getSurfaceNormal(), cx);
+		delete face;
+		if (s >= 0)
+			return false;
+	}
+	return true;
+}
+
+} // end namespace MeshLib
diff --git a/MeshLib/Elements/CellRule.h b/MeshLib/Elements/CellRule.h
new file mode 100644
index 00000000000..f8ae75fbbc8
--- /dev/null
+++ b/MeshLib/Elements/CellRule.h
@@ -0,0 +1,36 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef CELLRULE_H_
+#define CELLRULE_H_
+
+#include "MeshLib/MeshEnums.h"
+#include "Element.h"
+
+namespace MeshLib {
+
+/**
+ */
+class CellRule
+{
+public:
+	/// Constant: Dimension of this mesh element
+	static const unsigned dimension = 3u;
+
+	/**
+	 * Checks if the node order of an element is correct by testing surface normals.
+	 */
+	static bool testElementNodeOrder(const Element* e);
+
+}; /* class */
+
+} /* namespace */
+
+#endif /* HEXRULE_H_ */
+
diff --git a/MeshLib/Elements/Edge.h b/MeshLib/Elements/Edge.h
new file mode 100644
index 00000000000..d7cc1ff77d8
--- /dev/null
+++ b/MeshLib/Elements/Edge.h
@@ -0,0 +1,52 @@
+/**
+ * \file
+ * \author Karsten Rink
+ * \date   2012-05-02
+ * \brief  Definition of the Face class.
+ *
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef EDGE_H_
+#define EDGE_H_
+
+#include <limits>
+
+#include "Element.h"
+
+namespace MeshLib {
+
+/**
+ * Virtual base class for 1d mesh elements.
+ */
+class Edge : public Element
+{
+public:
+	/// Get the length of this 1d element.
+	virtual double getLength() const { return _content; }
+
+	/// Destructor
+	virtual ~Edge() {}
+
+	/**
+		* Checks if the node order of an element is correct by testing surface normals.
+		* For 1D elements this always returns true.
+		*/
+	virtual bool testElementNodeOrder() const { return true; }
+
+protected:
+	/// Constructor for a generic mesh element without an array of mesh nodes.
+	Edge(unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max()) : Element(value, id) {}
+
+
+}; /* class */
+
+} /* namespace */
+
+#endif /* EDGE_H_ */
+
diff --git a/MeshLib/Elements/EdgeReturn.cpp b/MeshLib/Elements/EdgeReturn.cpp
new file mode 100644
index 00000000000..f0b8bb71a97
--- /dev/null
+++ b/MeshLib/Elements/EdgeReturn.cpp
@@ -0,0 +1,49 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "EdgeReturn.h"
+
+#include "logog/include/logog.hpp"
+
+#include "MeshLib/Node.h"
+#include "Element.h"
+#include "Line.h"
+
+namespace MeshLib
+{
+
+const Element* LinearEdgeReturn::getEdge(const Element* e, unsigned i)
+{
+	if (i < e->getNEdges())
+	{
+		Node** nodes = new Node*[2];
+		nodes[0] = const_cast<Node*>(e->getEdgeNode(i,0));
+		nodes[1] = const_cast<Node*>(e->getEdgeNode(i,1));
+		return new Line(nodes);
+	}
+	ERR("Error in MeshLib::Element::getEdge() - Index does not exist.");
+	return nullptr;
+}
+
+const Element* QuadraticEdgeReturn::getEdge(const Element* e, unsigned i)
+{
+	if (i < e->getNEdges())
+	{
+		Node** nodes = new Node*[3];
+		nodes[0] = const_cast<Node*>(e->getEdgeNode(i,0));
+		nodes[1] = const_cast<Node*>(e->getEdgeNode(i,1));
+		nodes[2] = const_cast<Node*>(e->getEdgeNode(i,2));
+		return new Line3(nodes);
+	}
+	ERR("Error in MeshLib::Element::getEdge() - Index does not exist.");
+	return nullptr;
+}
+
+} // end MeshLib
+
diff --git a/MeshLib/Elements/EdgeReturn.h b/MeshLib/Elements/EdgeReturn.h
new file mode 100644
index 00000000000..6becd0fb45a
--- /dev/null
+++ b/MeshLib/Elements/EdgeReturn.h
@@ -0,0 +1,51 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef EDGERETURN_H_
+#define EDGERETURN_H_
+
+#include "logog/include/logog.hpp"
+
+namespace MeshLib
+{
+
+class Element;
+
+/// Returns always null pointer
+class DummyEdgeReturn
+{
+public:
+	/// Returns i-th edge of the given element
+	static const Element* getEdge(const Element* /*e*/, unsigned /*i*/)
+	{
+		return nullptr;
+	}
+};
+
+/// Returns linear order edge
+class LinearEdgeReturn
+{
+public:
+	/// Returns i-th edge of the given element
+	static const Element* getEdge(const Element* e, unsigned i);
+};
+
+/// Returns quadratic order edge
+class QuadraticEdgeReturn
+{
+public:
+	/// Returns i-th edge of the given element
+	static const Element* getEdge(const Element* e, unsigned i);
+};
+
+} // end MeshLib
+
+
+#endif /* EDGERETURN_H_ */
+
diff --git a/MeshLib/Elements/Element.cpp b/MeshLib/Elements/Element.cpp
index 14a3d088db2..06b4a328bca 100644
--- a/MeshLib/Elements/Element.cpp
+++ b/MeshLib/Elements/Element.cpp
@@ -24,7 +24,7 @@
 namespace MeshLib {
 
 Element::Element(unsigned value, std::size_t id)
-	: _nodes(nullptr), _id(id), _value(value), _neighbors(nullptr)
+	: _nodes(nullptr), _id(id), _content(-1.0), _value(value), _neighbors(nullptr)
 {
 }
 
diff --git a/MeshLib/Elements/Element.h b/MeshLib/Elements/Element.h
index 39425430159..bcffc2cd401 100644
--- a/MeshLib/Elements/Element.h
+++ b/MeshLib/Elements/Element.h
@@ -57,7 +57,7 @@ public:
 	MeshLib::Node getCenterOfGravity() const;
 
 	/// Returns the length, area or volume of a 1D, 2D or 3D element
-	virtual double getContent() const = 0;
+	double getContent() const { return _content; }
 
 	/**
 	 * Get node with local index i where i should be at most the number
@@ -84,7 +84,7 @@ public:
 	virtual unsigned getDimension() const = 0;
 
 	/// Returns the i-th edge of the element.
-	const Element* getEdge(unsigned i) const;
+	virtual const Element* getEdge(unsigned i) const;
 
 	/// Returns the i-th face of the element.
 	virtual const Element* getFace(unsigned i) const = 0;
@@ -207,6 +207,9 @@ public:
 	 */
 	virtual bool testElementNodeOrder() const = 0;
 
+	/// Return a specific edge node.
+	virtual Node* getEdgeNode(unsigned edge_id, unsigned node_id) const = 0;
+
 #ifndef NDEBUG
 	friend std::ostream& operator<<(std::ostream& os, Element const& e);
 #endif  // NDEBUG
@@ -215,14 +218,12 @@ protected:
 	/// Constructor for a generic mesh element without an array of mesh nodes.
 	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;
-
 	/// Sets the element ID.
 	virtual void setID(std::size_t id) { this->_id = id; }
 
 	Node** _nodes;
 	std::size_t _id;
+	double _content;
 	/**
 	 * this is an index for external additional information like materials
 	 */
diff --git a/MeshLib/Elements/Face.cpp b/MeshLib/Elements/Face.cpp
index 53886a08cdd..600468da609 100644
--- a/MeshLib/Elements/Face.cpp
+++ b/MeshLib/Elements/Face.cpp
@@ -21,11 +21,6 @@
 
 namespace MeshLib {
 
-#ifndef WIN32
-	/// \todo Windows compiler does not accept this definition and issues a linking error.
-	const unsigned Face::dimension;
-#endif	// WIN32
-
 /*
 Face::Face(Node** nodes, MeshElemType type, unsigned value)
 	: Element(nodes, type, value)
@@ -33,7 +28,7 @@ Face::Face(Node** nodes, MeshElemType type, unsigned value)
 }
 */
 Face::Face(unsigned value, std::size_t id)
-	: Element(value, id), _area(-1.0) // init with invalid value to detect errors
+	: Element(value, id)
 {
 }
 
diff --git a/MeshLib/Elements/Face.h b/MeshLib/Elements/Face.h
index 9050f626f82..14209d906d4 100644
--- a/MeshLib/Elements/Face.h
+++ b/MeshLib/Elements/Face.h
@@ -31,26 +31,8 @@ namespace MeshLib {
 class Face : public Element
 {
 public:
-	/// Constant: Dimension of this mesh element
-	static const unsigned dimension = 2;
-
 	/// Get the area of this 2d element.
-	virtual double getArea() const { return _area; }
-
-	/// Returns the length, area or volume of a 1D, 2D or 3D element
-	double getContent() const { return _area; }
-
-	/// Get dimension of the mesh element.
-	unsigned getDimension() const { return dimension; }
-
-	/// Returns the face i of the element.
-	const Element* getFace(unsigned i) const { return this->getEdge(i); }
-
-	/// Get the number of nodes for face i.
-	unsigned getNFaceNodes(unsigned i) const { (void)i; return 2; }
-
-	/// 2D elements have no faces.
-	unsigned getNFaces() const { return 0; }
+	virtual double getArea() const { return _content; }
 
 	/// Returns the surface normal of a 2D element.
 	MathLib::Vector3 getSurfaceNormal() const;
@@ -58,13 +40,6 @@ public:
 	/// Destructor
 	virtual ~Face();
 
-	/**
-	 * This method is pure virtual and is inherited from class @sa Element.
-	 * It has to be implemented in the derived classes of class Face!
-	 * @return a copy of the object
-	 */
-	virtual Element* clone() const = 0;
-
 	/**
 	 * Checks if the node order of an element is correct by testing surface normals.
 	 * For 2D elements true is returned if the normal points (roughly) upwards.
@@ -79,8 +54,6 @@ protected:
 	/// Constructor for a generic mesh element without an array of mesh nodes.
 	Face(unsigned value = 0, std::size_t id = std::numeric_limits<std::size_t>::max());
 
-	double _area;
-
 private:
 
 
diff --git a/MeshLib/Elements/FaceRule.cpp b/MeshLib/Elements/FaceRule.cpp
new file mode 100644
index 00000000000..fb12a1fded1
--- /dev/null
+++ b/MeshLib/Elements/FaceRule.cpp
@@ -0,0 +1,17 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "FaceRule.h"
+
+namespace MeshLib {
+
+const unsigned FaceRule::dimension;
+
+} /* namespace */
+
diff --git a/MeshLib/Elements/FaceRule.h b/MeshLib/Elements/FaceRule.h
new file mode 100644
index 00000000000..5dc7da7971c
--- /dev/null
+++ b/MeshLib/Elements/FaceRule.h
@@ -0,0 +1,39 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef FACERULE_H_
+#define FACERULE_H_
+
+#include "MeshLib/MeshEnums.h"
+#include "Element.h"
+
+namespace MeshLib {
+
+/**
+ */
+class FaceRule
+{
+public:
+	/// Constant: Dimension of this mesh element
+	static const unsigned dimension = 2;
+
+	/// Returns the face i of the element.
+	static const Element* getFace(const Element* e, unsigned i) { return e->getEdge(i); }
+
+	/// Get the number of nodes for face i.
+	static unsigned getNFaceNodes(unsigned /*i*/) { return 2; }
+
+	/// 2D elements have no faces.
+	static unsigned getNFaces() { return 0; }
+}; /* class */
+
+} /* namespace */
+
+#endif /* FACERULE_H_ */
+
diff --git a/MeshLib/Elements/Hex.h b/MeshLib/Elements/Hex.h
index 8031d29073d..08be8f713d4 100644
--- a/MeshLib/Elements/Hex.h
+++ b/MeshLib/Elements/Hex.h
@@ -15,10 +15,12 @@
 #ifndef HEX_H_
 #define HEX_H_
 
-#include "TemplateHex.h"
+#include "TemplateElement.h"
+#include "Cell.h"
+#include "HexRule8.h"
 
 namespace MeshLib {
-typedef TemplateHex<8, CellType::HEX8> Hex;
+typedef TemplateElement<Cell, HexRule8> Hex;
 }
 
 #endif /* HEX_H_ */
diff --git a/MeshLib/Elements/HexRule8.cpp b/MeshLib/Elements/HexRule8.cpp
new file mode 100644
index 00000000000..43adf1dae3d
--- /dev/null
+++ b/MeshLib/Elements/HexRule8.cpp
@@ -0,0 +1,119 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "HexRule8.h"
+
+#include "logog/include/logog.hpp"
+
+#include "GeoLib/AnalyticalGeometry.h"
+
+#include "MeshLib/Node.h"
+#include "Quad.h"
+#include "Line.h"
+
+namespace MeshLib {
+
+const unsigned HexRule8::n_all_nodes;
+
+const unsigned HexRule8::n_base_nodes;
+
+const unsigned HexRule8::face_nodes[6][4] =
+{
+	{0, 3, 2, 1}, // Face 0
+	{0, 1, 5, 4}, // Face 1
+	{1, 2, 6, 5}, // Face 2
+	{2, 3, 7, 6}, // Face 3
+	{3, 0, 4, 7}, // Face 4
+	{4, 5, 6, 7}  // Face 5
+};
+
+const unsigned HexRule8::edge_nodes[12][2] =
+{
+	{0, 1}, // Edge 0
+	{1, 2}, // Edge 1
+	{2, 3}, // Edge 2
+	{0, 3}, // Edge 3
+	{4, 5}, // Edge 4
+	{5, 6}, // Edge 5
+	{6, 7}, // Edge 6
+	{4, 7}, // Edge 7
+	{0, 4}, // Edge 8
+	{1, 5}, // Edge 9
+	{2, 6}, // Edge 10
+	{3, 7}  // Edge 11
+};
+
+const Element* HexRule8::getFace(const Element* e, unsigned i)
+{
+	if (i < n_faces)
+	{
+		unsigned nFaceNodes (getNFaceNodes(i));
+		Node** nodes = new Node*[nFaceNodes];
+		for (unsigned j=0; j<nFaceNodes; j++)
+			nodes[j] = const_cast<Node*>(e->getNode(face_nodes[i][j]));
+		return new Quad(nodes);
+	}
+	ERR("Error in MeshLib::Element::getFace() - Index %d does not exist.", i);
+	return nullptr;
+}
+
+double HexRule8::computeVolume(Node const* const* _nodes)
+{
+	return GeoLib::calcTetrahedronVolume(_nodes[4]->getCoords(), _nodes[7]->getCoords(), _nodes[5]->getCoords(), _nodes[0]->getCoords())
+		 + GeoLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[3]->getCoords(), _nodes[1]->getCoords(), _nodes[0]->getCoords())
+		 + GeoLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[7]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords())
+		 + GeoLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[7]->getCoords(), _nodes[6]->getCoords(), _nodes[2]->getCoords())
+		 + GeoLib::calcTetrahedronVolume(_nodes[1]->getCoords(), _nodes[3]->getCoords(), _nodes[5]->getCoords(), _nodes[2]->getCoords())
+		 + GeoLib::calcTetrahedronVolume(_nodes[3]->getCoords(), _nodes[7]->getCoords(), _nodes[5]->getCoords(), _nodes[2]->getCoords());
+}
+
+bool HexRule8::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps)
+{
+	return (GeoLib::isPointInTetrahedron(pnt, *_nodes[4], *_nodes[7], *_nodes[5], *_nodes[0], eps) ||
+			GeoLib::isPointInTetrahedron(pnt, *_nodes[5], *_nodes[3], *_nodes[1], *_nodes[0], eps) ||
+			GeoLib::isPointInTetrahedron(pnt, *_nodes[5], *_nodes[7], *_nodes[3], *_nodes[0], eps) ||
+			GeoLib::isPointInTetrahedron(pnt, *_nodes[5], *_nodes[7], *_nodes[6], *_nodes[2], eps) ||
+			GeoLib::isPointInTetrahedron(pnt, *_nodes[1], *_nodes[3], *_nodes[5], *_nodes[2], eps) ||
+			GeoLib::isPointInTetrahedron(pnt, *_nodes[3], *_nodes[7], *_nodes[5], *_nodes[2], eps));
+}
+
+unsigned HexRule8::identifyFace(Node const* const* _nodes, Node* nodes[3])
+{
+	for (unsigned i=0; i<6; i++)
+	{
+		unsigned flag(0);
+		for (unsigned j=0; j<4; j++)
+			for (unsigned k=0; k<3; k++)
+				if (_nodes[face_nodes[i][j]] == nodes[k])
+					flag++;
+		if (flag==3)
+			return i;
+	}
+	return std::numeric_limits<unsigned>::max();
+}
+
+ElementErrorCode HexRule8::validate(const Element* e)
+{
+	ElementErrorCode error_code;
+	error_code[ElementErrorFlag::ZeroVolume] = e->hasZeroVolume();
+		
+	for (unsigned i=0; i<6; ++i)
+	{
+		if (error_code.all())
+			break;
+
+		const MeshLib::Element* quad (e->getFace(i));
+		error_code |= quad->validate();
+		delete quad;
+	}
+	error_code[ElementErrorFlag::NodeOrder]  = !e->testElementNodeOrder();
+	return error_code;
+}
+
+} // end namespace MeshLib
diff --git a/MeshLib/Elements/HexRule8.h b/MeshLib/Elements/HexRule8.h
new file mode 100644
index 00000000000..80169dc8953
--- /dev/null
+++ b/MeshLib/Elements/HexRule8.h
@@ -0,0 +1,110 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef HEXRULE8_H_
+#define HEXRULE8_H_
+
+#include "MeshLib/MeshEnums.h"
+#include "Element.h"
+#include "CellRule.h"
+#include "EdgeReturn.h"
+
+namespace MeshLib
+{
+
+/**
+ * A 8-nodes Hexahedron Element.
+ * @code
+ *
+ *  Hex:
+ *                6
+ *          7-----------6
+ *         /:          /|
+ *        / :         / |
+ *      7/  :        /5 |
+ *      / 11:       /   | 10
+ *     /    : 4    /    |
+ *    4-----------5     |
+ *    |     :     | 2   |
+ *    |     3.....|.....2
+ *    |    .      |    /
+ *  8 |   .       |9  /
+ *    | 3.        |  / 1
+ *    | .         | /
+ *    |.          |/
+ *    0-----------1
+ *          0
+ *
+ * @endcode
+ */
+class HexRule8 : public CellRule
+{
+public:
+	/// Constant: The number of base nodes for this element
+	static const unsigned n_base_nodes = 8u;
+
+	/// Constant: The number of all nodes for this element
+	static const unsigned n_all_nodes = 8u;
+
+	/// Constant: The geometric type of the element
+	static const MeshElemType mesh_elem_type = MeshElemType::HEXAHEDRON;
+
+	/// Constant: The FEM type of the element
+	static const CellType cell_type = CellType::HEX8;
+
+	/// Constant: The number of faces
+	static const unsigned n_faces = 6;
+
+	/// Constant: The number of edges
+	static const unsigned n_edges = 12;
+
+	/// Constant: The number of neighbors
+	static const unsigned n_neighbors = 6;
+
+	/// Constant: Local node index table for faces
+	static const unsigned face_nodes[6][4];
+
+	/// Constant: Local node index table for edge
+	static const unsigned edge_nodes[12][2];
+
+	/// Returns the i-th edge of the element.
+	typedef LinearEdgeReturn EdgeReturn;
+
+	/// Get the number of nodes for face i.
+	static unsigned getNFaceNodes(unsigned /*i*/) { return 4; }
+
+	/// Returns the i-th face of the element.
+	static const Element* getFace(const Element* e, unsigned i);
+
+	/**
+	 * Checks if a point is inside the element.
+	 * @param pnt a 3D GeoLib::Point object
+	 * @param eps tolerance for numerical algorithm used or computing the property
+	 * @return true if the point is not outside the element, false otherwise
+	 */
+	static bool isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps);
+
+	/**
+	 * Tests if the element is geometrically valid.
+	 * @param check_zero_volume indicates if volume == 0 should be checked
+	 */
+	static ElementErrorCode validate(const Element* e);
+
+	/// Returns the ID of a face given an array of nodes.
+	static unsigned identifyFace(Node const* const*, Node* nodes[3]);
+
+	/// Calculates the volume of a convex hexahedron by partitioning it into six tetrahedra.
+	static double computeVolume(Node const* const* _nodes);
+
+}; /* class */
+
+} /* namespace */
+
+#endif /* HEXRULE_H_ */
+
diff --git a/MeshLib/Elements/Line.h b/MeshLib/Elements/Line.h
index f236ac82ea0..bb5b9be5d0d 100644
--- a/MeshLib/Elements/Line.h
+++ b/MeshLib/Elements/Line.h
@@ -15,11 +15,13 @@
 #ifndef LINE_H_
 #define LINE_H_
 
-#include "TemplateLine.h"
+#include "TemplateElement.h"
+#include "Edge.h"
+#include "LineRule2.h"
 
 namespace MeshLib {
 
-typedef TemplateLine<2,CellType::LINE2> Line;
+typedef TemplateElement<Edge, LineRule2> Line;
 
 }
 
diff --git a/MeshLib/Elements/LineRule2.cpp b/MeshLib/Elements/LineRule2.cpp
new file mode 100644
index 00000000000..40022b29925
--- /dev/null
+++ b/MeshLib/Elements/LineRule2.cpp
@@ -0,0 +1,60 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "LineRule2.h"
+
+#include "logog/include/logog.hpp"
+
+#include "GeoLib/AnalyticalGeometry.h"
+
+#include "MeshLib/Node.h"
+
+namespace MeshLib {
+
+const unsigned LineRule2::n_all_nodes;
+
+const unsigned LineRule2::n_base_nodes;
+
+const unsigned LineRule2::dimension;
+
+const unsigned LineRule2::edge_nodes[1][2] =
+{
+	{0, 1} // Edge 0
+};
+
+double LineRule2::computeVolume(Node const* const* _nodes)
+{
+	return sqrt(MathLib::sqrDist(_nodes[0]->getCoords(), _nodes[1]->getCoords()));
+}
+
+bool LineRule2::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps)
+{
+	double tmp;
+	double tmp_dst(0);
+	double const dist =MathLib::calcProjPntToLineAndDists(pnt.getCoords(), _nodes[0]->getCoords(), _nodes[1]->getCoords(), tmp, tmp_dst);
+	return (dist < eps);
+}
+
+unsigned LineRule2::identifyFace(Node const* const* _nodes, Node* nodes[1])
+{
+	if (nodes[0] == _nodes[0])
+		return 0;
+	if (nodes[0] == _nodes[1])
+		return 1;
+	return std::numeric_limits<unsigned>::max();
+}
+
+ElementErrorCode LineRule2::validate(const Element* e)
+{
+	ElementErrorCode error_code;
+	error_code[ElementErrorFlag::ZeroVolume] = e->hasZeroVolume();
+	return error_code;
+}
+
+} // end namespace MeshLib
diff --git a/MeshLib/Elements/LineRule2.h b/MeshLib/Elements/LineRule2.h
new file mode 100644
index 00000000000..6dc8ec0241c
--- /dev/null
+++ b/MeshLib/Elements/LineRule2.h
@@ -0,0 +1,90 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef LINERULE2_H_
+#define LINERULE2_H_
+
+#include "MeshLib/MeshEnums.h"
+#include "Element.h"
+#include "EdgeReturn.h"
+
+namespace MeshLib
+{
+
+/**
+ * A 1d Edge or Line Element.
+ * @code
+ *  0--------1
+ * @endcode
+ */
+class LineRule2
+{
+public:
+	/// Constant: The number of base nodes for this element
+	static const unsigned n_base_nodes = 2u;
+
+	/// Constant: The number of all nodes for this element
+	static const unsigned n_all_nodes = 2u;
+
+	/// Constant: Dimension of this mesh element
+	static const unsigned dimension = 1;
+
+	/// Constant: The geometric type of the element
+	static const MeshElemType mesh_elem_type = MeshElemType::LINE;
+
+	/// Constant: The FEM type of the element
+	static const CellType cell_type = CellType::LINE2;
+
+	/// Constant: The number of faces
+	static const unsigned n_faces = 0;
+
+	/// Constant: The number of edges
+	static const unsigned n_edges = 0;
+
+	/// Constant: The number of neighbors
+	static const unsigned n_neighbors = 2;
+
+	/// Constant: Local node index table for edge
+	static const unsigned edge_nodes[1][2];
+
+	/// Edge rule
+	typedef DummyEdgeReturn EdgeReturn;
+
+	/// Get the number of nodes for face i.
+	static unsigned getNFaceNodes(unsigned /*i*/) { return 0; }
+
+	/// Returns the i-th face of the element.
+	static const Element* getFace(const Element* /*e*/, unsigned /*i*/) { return nullptr; }
+
+	/**
+	 * Checks if a point is inside the element.
+	 * @param pnt a 3D GeoLib::Point object
+	 * @param eps tolerance for numerical algorithm used or computing the property
+	 * @return true if the point is not outside the element, false otherwise
+	 */
+	static bool isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps);
+
+	/**
+	 * Tests if the element is geometrically valid.
+	 * @param check_zero_volume indicates if volume == 0 should be checked
+	 */
+	static ElementErrorCode validate(const Element* e);
+
+	/// Returns the ID of a face given an array of nodes.
+	static unsigned identifyFace(Node const* const*, Node* nodes[1]);
+
+	/// Calculates the length of a line
+	static double computeVolume(Node const* const* _nodes);
+
+}; /* class */
+
+} /* namespace */
+
+#endif /* LINERULE2_H_ */
+
diff --git a/MeshLib/Elements/Prism.h b/MeshLib/Elements/Prism.h
index 0dbfa0262d6..5fc96f59d37 100644
--- a/MeshLib/Elements/Prism.h
+++ b/MeshLib/Elements/Prism.h
@@ -15,11 +15,13 @@
 #ifndef PRISM_H_
 #define PRISM_H_
 
-#include "TemplatePrism.h"
+#include "TemplateElement.h"
+#include "Cell.h"
+#include "PrismRule6.h"
 
 namespace MeshLib {
 
-typedef TemplatePrism<6, CellType::PRISM6> Prism;
+typedef TemplateElement<Cell, PrismRule6> Prism;
 
 }
 
diff --git a/MeshLib/Elements/PrismRule6.cpp b/MeshLib/Elements/PrismRule6.cpp
new file mode 100644
index 00000000000..3c38deac334
--- /dev/null
+++ b/MeshLib/Elements/PrismRule6.cpp
@@ -0,0 +1,123 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "PrismRule6.h"
+
+#include "logog/include/logog.hpp"
+
+#include "GeoLib/AnalyticalGeometry.h"
+
+#include "MeshLib/Node.h"
+#include "Quad.h"
+#include "Tri.h"
+
+namespace MeshLib {
+
+const unsigned PrismRule6::n_all_nodes;
+
+const unsigned PrismRule6::n_base_nodes;
+
+const unsigned PrismRule6::face_nodes[5][4] =
+{
+	{0, 2, 1, 99}, // Face 0
+	{0, 1, 4,  3}, // Face 1
+	{1, 2, 5,  4}, // Face 2
+	{2, 0, 3,  5}, // Face 3
+	{3, 4, 5, 99}  // Face 4
+};
+
+const unsigned PrismRule6::edge_nodes[9][2] =
+{
+	{0, 1}, // Edge 0
+	{1, 2}, // Edge 1
+	{0, 2}, // Edge 2
+	{0, 3}, // Edge 3
+	{1, 4}, // Edge 4
+	{2, 5}, // Edge 5
+	{3, 4}, // Edge 6
+	{4, 5}, // Edge 7
+	{3, 5}  // Edge 8
+};
+
+const unsigned PrismRule6::n_face_nodes[5] = { 3, 4, 4, 4, 3 };
+
+unsigned PrismRule6::getNFaceNodes(unsigned i)
+{
+	if (i<5)
+		return n_face_nodes[i];
+	ERR("Error in MeshLib::Element::getNFaceNodes() - Index %d does not exist.", i);
+	return 0;
+}
+
+const Element* PrismRule6::getFace(const Element* e, unsigned i)
+{
+	if (i < n_faces)
+	{
+		unsigned nFaceNodes (e->getNFaceNodes(i));
+		Node** nodes = new Node*[nFaceNodes];
+		for (unsigned j=0; j<nFaceNodes; j++)
+			nodes[j] = const_cast<Node*>(e->getNode(face_nodes[i][j]));
+
+		if (i==0 || i==4)
+			return new Tri(nodes);
+		else
+			return new Quad(nodes);
+	}
+	ERR("Error in MeshLib::Element::getFace() - Index %d does not exist.", i);
+	return nullptr;
+}
+
+double PrismRule6::computeVolume(Node const* const* _nodes)
+{
+	return GeoLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords())
+		 + GeoLib::calcTetrahedronVolume(_nodes[1]->getCoords(), _nodes[4]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords())
+		 + GeoLib::calcTetrahedronVolume(_nodes[2]->getCoords(), _nodes[4]->getCoords(), _nodes[5]->getCoords(), _nodes[3]->getCoords());
+}
+
+bool PrismRule6::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps)
+{
+	return (GeoLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3], eps) ||
+			GeoLib::isPointInTetrahedron(pnt, *_nodes[1], *_nodes[4], *_nodes[2], *_nodes[3], eps) ||
+			GeoLib::isPointInTetrahedron(pnt, *_nodes[2], *_nodes[4], *_nodes[5], *_nodes[3], eps));
+}
+
+unsigned PrismRule6::identifyFace(Node const* const* _nodes, Node* nodes[3])
+{
+	for (unsigned i=0; i<5; i++)
+	{
+		unsigned flag(0);
+		for (unsigned j=0; j<4; j++)
+			for (unsigned k=0; k<3; k++)
+				if (face_nodes[i][j] != 99 && _nodes[face_nodes[i][j]] == nodes[k])
+					flag++;
+		if (flag==3)
+			return i;
+	}
+	return std::numeric_limits<unsigned>::max();
+}
+
+ElementErrorCode PrismRule6::validate(const Element* e)
+{
+	ElementErrorCode error_code;
+	error_code[ElementErrorFlag::ZeroVolume] = e->hasZeroVolume();
+
+	for (unsigned i=1; i<4; ++i)
+	{
+		const MeshLib::Quad* quad (dynamic_cast<const MeshLib::Quad*>(e->getFace(i)));
+		if (quad)
+			error_code |= quad->validate();
+		else
+			error_code.set(ElementErrorFlag::NodeOrder);
+		delete quad;
+	}
+	error_code[ElementErrorFlag::NodeOrder] = !e->testElementNodeOrder();
+	return error_code;
+}
+
+} // end namespace MeshLib
diff --git a/MeshLib/Elements/PrismRule6.h b/MeshLib/Elements/PrismRule6.h
new file mode 100644
index 00000000000..f533939681f
--- /dev/null
+++ b/MeshLib/Elements/PrismRule6.h
@@ -0,0 +1,111 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PRISMRULE6_H_
+#define PRISMRULE6_H_
+
+#include "MeshLib/MeshEnums.h"
+#include "Element.h"
+#include "CellRule.h"
+#include "EdgeReturn.h"
+
+namespace MeshLib
+{
+
+/**
+ * This class represents a 3d prism element. The following sketch shows the node and edge numbering.
+ * @anchor PrismNodeAndEdgeNumbering
+ * @code
+ *            5
+ *           / \
+ *          / : \
+ *        8/  :  \7
+ *        /   :5  \
+ *       /    :  6 \
+ *      3-----------4
+ *      |     :     |
+ *      |     2     |
+ *      |    . .    |
+ *     3|   .   .   |4
+ *      | 2.     .1 |
+ *      | .       . |
+ *      |.         .|
+ *      0-----------1
+ *            0
+ *
+ * @endcode
+ */
+class PrismRule6 : public CellRule
+{
+public:
+	/// Constant: The number of base nodes for this element
+	static const unsigned n_base_nodes = 6u;
+
+	/// Constant: The number of all nodes for this element
+	static const unsigned n_all_nodes = 6u;
+
+	/// Constant: The geometric type of the element
+	static const MeshElemType mesh_elem_type = MeshElemType::PRISM;
+
+	/// Constant: The FEM type of the element
+	static const CellType cell_type = CellType::PRISM6;
+
+	/// Constant: The number of faces
+	static const unsigned n_faces = 5;
+
+	/// Constant: The number of edges
+	static const unsigned n_edges = 9;
+
+	/// Constant: The number of neighbors
+	static const unsigned n_neighbors = 5;
+
+	/// Constant: Local node index table for faces
+	static const unsigned face_nodes[5][4];
+
+	/// Constant: Local node index table for edge
+	static const unsigned edge_nodes[9][2];
+
+	/// Constant: Table for the number of nodes for each face
+	static const unsigned n_face_nodes[5];
+
+	/// Returns the i-th edge of the element.
+	typedef LinearEdgeReturn EdgeReturn;
+
+	/// Get the number of nodes for face i.
+	static unsigned getNFaceNodes(unsigned i);
+
+	/// Returns the i-th face of the element.
+	static const Element* getFace(const Element* e, unsigned i);
+
+	/**
+	 * Checks if a point is inside the element.
+	 * @param pnt a 3D GeoLib::Point object
+	 * @param eps tolerance for numerical algorithm used or computing the property
+	 * @return true if the point is not outside the element, false otherwise
+	 */
+	static bool isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps);
+
+	/**
+	 * Tests if the element is geometrically valid.
+	 * @param check_zero_volume indicates if volume == 0 should be checked
+	 */
+	static ElementErrorCode validate(const Element* e);
+
+	/// Returns the ID of a face given an array of nodes.
+	static unsigned identifyFace(Node const* const*, Node* nodes[3]);
+
+	/// Calculates the volume of a convex hexahedron by partitioning it into six tetrahedra.
+	static double computeVolume(Node const* const* _nodes);
+
+}; /* class */
+
+} /* namespace */
+
+#endif /* PRISMRULE6_H_ */
+
diff --git a/MeshLib/Elements/Pyramid.h b/MeshLib/Elements/Pyramid.h
index de095cf13d5..42cf46a813a 100644
--- a/MeshLib/Elements/Pyramid.h
+++ b/MeshLib/Elements/Pyramid.h
@@ -15,11 +15,13 @@
 #ifndef PYRAMID_H_
 #define PYRAMID_H_
 
-#include "TemplatePyramid.h"
+#include "TemplateElement.h"
+#include "Cell.h"
+#include "PyramidRule5.h"
 
 namespace MeshLib {
 
-typedef TemplatePyramid<5,CellType::PYRAMID5> Pyramid;
+typedef TemplateElement<Cell, PyramidRule5> Pyramid;
 
 }
 
diff --git a/MeshLib/Elements/PyramidRule5.cpp b/MeshLib/Elements/PyramidRule5.cpp
new file mode 100644
index 00000000000..fee75fa986f
--- /dev/null
+++ b/MeshLib/Elements/PyramidRule5.cpp
@@ -0,0 +1,120 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "PyramidRule5.h"
+
+#include "logog/include/logog.hpp"
+
+#include "GeoLib/AnalyticalGeometry.h"
+
+#include "MeshLib/Node.h"
+#include "Quad.h"
+#include "Tri.h"
+
+namespace MeshLib {
+
+const unsigned PyramidRule5::n_all_nodes;
+
+const unsigned PyramidRule5::n_base_nodes;
+
+const unsigned PyramidRule5::face_nodes[5][4] =
+{
+	{0, 1, 4, 99}, // Face 0
+	{1, 2, 4, 99}, // Face 1
+	{2, 3, 4, 99}, // Face 2
+	{3, 0, 4, 99}, // Face 3
+	{0, 3, 2,  1}  // Face 4
+};
+
+const unsigned PyramidRule5::edge_nodes[8][2] =
+{
+	{0, 1}, // Edge 0
+	{1, 2}, // Edge 1
+	{2, 3}, // Edge 2
+	{0, 3}, // Edge 3
+	{0, 4}, // Edge 4
+	{1, 4}, // Edge 5
+	{2, 4}, // Edge 6
+	{3, 4}  // Edge 7
+};
+
+const unsigned PyramidRule5::n_face_nodes[5] = { 3, 3, 3, 3, 4 };
+
+unsigned PyramidRule5::getNFaceNodes(unsigned i)
+{
+	if (i<5)
+		return n_face_nodes[i];
+	ERR("Error in MeshLib::Element::getNFaceNodes() - Index %d does not exist.", i);
+	return 0;
+}
+
+const Element* PyramidRule5::getFace(const Element* e, unsigned i)
+{
+	if (i<e->getNFaces())
+	{
+		unsigned nFaceNodes (e->getNFaceNodes(i));
+		Node** nodes = new Node*[nFaceNodes];
+		for (unsigned j=0; j<nFaceNodes; j++)
+			nodes[j] = const_cast<Node*>(e->getNode(face_nodes[i][j]));
+
+		if (i<4)
+			return new Tri(nodes);
+		else
+			return new Quad(nodes);
+	}
+	ERR("Error in MeshLib::Element::getFace() - Index %d does not exist.", i);
+	return nullptr;
+}
+
+double PyramidRule5::computeVolume(Node const* const* _nodes)
+{
+	return GeoLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[4]->getCoords())
+		 + GeoLib::calcTetrahedronVolume(_nodes[2]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords(), _nodes[4]->getCoords());
+}
+
+bool PyramidRule5::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps)
+{
+	return (GeoLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[4], eps) ||
+			GeoLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[2], *_nodes[3], *_nodes[4], eps));
+}
+
+unsigned PyramidRule5::identifyFace(Node const* const* _nodes, Node* nodes[3])
+{
+	for (unsigned i=0; i<5; i++)
+	{
+		unsigned flag(0);
+		for (unsigned j=0; j<4; j++)
+			for (unsigned k=0; k<3; k++)
+				if (face_nodes[i][j] != 99 && _nodes[face_nodes[i][j]] == nodes[k])
+					flag++;
+		if (flag==3)
+			return i;
+	}
+	return std::numeric_limits<unsigned>::max();
+}
+
+ElementErrorCode PyramidRule5::validate(const Element* e)
+{
+	ElementErrorCode error_code;
+	error_code[ElementErrorFlag::ZeroVolume] = e->hasZeroVolume();
+
+	const MeshLib::Quad* base (dynamic_cast<const MeshLib::Quad*>(e->getFace(4)));
+	if (base)
+	{
+		error_code |= base->validate();
+		error_code[ElementErrorFlag::NodeOrder] = !e->testElementNodeOrder();
+	}
+	else
+		error_code.set(ElementErrorFlag::NodeOrder);
+	delete base;
+
+	return error_code;
+}
+
+} // end namespace MeshLib
diff --git a/MeshLib/Elements/PyramidRule5.h b/MeshLib/Elements/PyramidRule5.h
new file mode 100644
index 00000000000..069bdec6105
--- /dev/null
+++ b/MeshLib/Elements/PyramidRule5.h
@@ -0,0 +1,110 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PYRAMIDRULE5_H_
+#define PYRAMIDRULE5_H_
+
+#include "MeshLib/MeshEnums.h"
+#include "Element.h"
+#include "CellRule.h"
+#include "EdgeReturn.h"
+
+namespace MeshLib
+{
+
+/**
+ * This class represents a 3d pyramid element. The following sketch shows the node and edge numbering.
+ * @anchor PyramidNodeAndEdgeNumbering
+ * @code
+ *
+ *               4
+ *             //|\
+ *            // | \
+ *          7//  |  \6
+ *          //   |5  \
+ *         //    |    \
+ *        3/.... |.....2
+ *       ./      |  2 /
+ *      ./4      |   /
+ *    3./        |  /1
+ *    ./         | /
+ *   ./          |/
+ *  0------------1
+ *        0
+ * @endcode
+
+ */
+class PyramidRule5 : public CellRule
+{
+public:
+	/// Constant: The number of base nodes for this element
+	static const unsigned n_base_nodes = 5u;
+
+	/// Constant: The number of all nodes for this element
+	static const unsigned n_all_nodes = 5u;
+
+	/// Constant: The geometric type of the element
+	static const MeshElemType mesh_elem_type = MeshElemType::PYRAMID;
+
+	/// Constant: The FEM type of the element
+	static const CellType cell_type = CellType::PYRAMID5;
+
+	/// Constant: The number of faces
+	static const unsigned n_faces = 5;
+
+	/// Constant: The number of edges
+	static const unsigned n_edges = 8;
+
+	/// Constant: The number of neighbors
+	static const unsigned n_neighbors = 5;
+
+	/// Constant: Local node index table for faces
+	static const unsigned face_nodes[5][4];
+
+	/// Constant: Local node index table for edge
+	static const unsigned edge_nodes[8][2];
+
+	/// Constant: Table for the number of nodes for each face
+	static const unsigned n_face_nodes[5];
+
+	/// Returns the i-th edge of the element.
+	typedef LinearEdgeReturn EdgeReturn;
+
+	/// Get the number of nodes for face i.
+	static unsigned getNFaceNodes(unsigned i);
+
+	/// Returns the i-th face of the element.
+	static const Element* getFace(const Element* e, unsigned i);
+
+	/**
+	 * Checks if a point is inside the element.
+	 * @param pnt a 3D GeoLib::Point object
+	 * @param eps tolerance for numerical algorithm used or computing the property
+	 * @return true if the point is not outside the element, false otherwise
+	 */
+	static bool isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps);
+
+	/**
+	 * Tests if the element is geometrically valid.
+	 * @param check_zero_volume indicates if volume == 0 should be checked
+	 */
+	static ElementErrorCode validate(const Element* e);
+
+	/// Returns the ID of a face given an array of nodes.
+	static unsigned identifyFace(Node const* const*, Node* nodes[3]);
+
+	/// Calculates the volume of a convex hexahedron by partitioning it into six tetrahedra.
+	static double computeVolume(Node const* const* _nodes);
+
+}; /* class */
+
+} /* namespace */
+
+#endif /* PYRAMIDRULE5_H_ */
+
diff --git a/MeshLib/Elements/Quad.h b/MeshLib/Elements/Quad.h
index f7e4a6112f9..482a49f973f 100644
--- a/MeshLib/Elements/Quad.h
+++ b/MeshLib/Elements/Quad.h
@@ -15,13 +15,18 @@
 #ifndef QUAD_H_
 #define QUAD_H_
 
-#include "TemplateQuad.h"
+#include "TemplateElement.h"
+#include "Face.h"
+#include "QuadRule4.h"
+#include "QuadRule8.h"
+#include "QuadRule9.h"
 
-namespace MeshLib {
+namespace MeshLib
+{
 
-typedef TemplateQuad<4,CellType::QUAD4> Quad;
-typedef TemplateQuad<8, CellType::QUAD8> Quad8;
-typedef TemplateQuad<9, CellType::QUAD9> Quad9;
+typedef TemplateElement<Face, QuadRule4> Quad;
+typedef TemplateElement<Face, QuadRule8> Quad8;
+typedef TemplateElement<Face, QuadRule9> Quad9;
 
 }
 
diff --git a/MeshLib/Elements/QuadRule4.cpp b/MeshLib/Elements/QuadRule4.cpp
new file mode 100644
index 00000000000..c2c2325de99
--- /dev/null
+++ b/MeshLib/Elements/QuadRule4.cpp
@@ -0,0 +1,73 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "QuadRule4.h"
+
+#include "logog/include/logog.hpp"
+
+#include "GeoLib/AnalyticalGeometry.h"
+
+#include "MeshLib/Node.h"
+
+namespace MeshLib {
+
+const unsigned QuadRule4::n_all_nodes;
+
+const unsigned QuadRule4::n_base_nodes;
+
+const unsigned QuadRule4::edge_nodes[4][2] =
+{
+		{0, 1}, // Edge 0
+		{1, 2}, // Edge 1
+		{2, 3}, // Edge 2
+		{0, 3}  // Edge 3
+};
+
+double QuadRule4::computeVolume(Node const* const* _nodes)
+{
+	return GeoLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2])
+         + GeoLib::calcTriangleArea(*_nodes[2], *_nodes[3], *_nodes[0]);
+}
+
+bool QuadRule4::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps)
+{
+	return (GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[1], *_nodes[2], eps) ||
+	        GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[2], *_nodes[3], eps));
+}
+
+unsigned QuadRule4::identifyFace(Node const* const* _nodes, Node* nodes[3])
+{
+	for (unsigned i=0; i<4; i++)
+	{
+		unsigned flag(0);
+		for (unsigned j=0; j<2; j++)
+			for (unsigned k=0; k<2; k++)
+				if (_nodes[edge_nodes[i][j]] == nodes[k])
+					flag++;
+		if (flag==2)
+			return i;
+	}
+	return std::numeric_limits<unsigned>::max();
+}
+
+ElementErrorCode QuadRule4::validate(const Element* e)
+{
+	ElementErrorCode error_code;
+	error_code[ElementErrorFlag::ZeroVolume]  = e->hasZeroVolume();
+	Node const* const* _nodes = e->getNodes();
+	error_code[ElementErrorFlag::NonCoplanar] = (!GeoLib::isCoplanar(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3]));
+	// for collapsed quads (i.e. reduced to a line) this test might result "false" as all four points are actually located on a line.
+	if (!error_code[ElementErrorFlag::ZeroVolume])
+		error_code[ElementErrorFlag::NonConvex]   = (!(GeoLib::dividedByPlane(*_nodes[0], *_nodes[2], *_nodes[1], *_nodes[3]) &&
+			                                           GeoLib::dividedByPlane(*_nodes[1], *_nodes[3], *_nodes[0], *_nodes[2])));
+	error_code[ElementErrorFlag::NodeOrder]  = !e->testElementNodeOrder();
+	return error_code;
+}
+
+} // end namespace MeshLib
diff --git a/MeshLib/Elements/QuadRule4.h b/MeshLib/Elements/QuadRule4.h
new file mode 100644
index 00000000000..57386880994
--- /dev/null
+++ b/MeshLib/Elements/QuadRule4.h
@@ -0,0 +1,92 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef QUADRULE4_H_
+#define QUADRULE4_H_
+
+#include "MeshLib/MeshEnums.h"
+#include "Element.h"
+#include "FaceRule.h"
+#include "EdgeReturn.h"
+#include "Line.h"
+
+namespace MeshLib
+{
+
+/**
+ * This class represents a 2d quadliteral element. The following sketch shows the node and edge numbering.
+ * @anchor QuadNodeAndEdgeNumbering
+ * @code
+ *              2
+ *        3-----------2
+ *        |           |
+ *        |           |
+ *       3|           |1
+ *        |           |
+ *        |           |
+ *        0-----------1
+ *              0
+ * @endcode
+ */
+class QuadRule4 : public FaceRule
+{
+public:
+	/// Constant: The number of base nodes for this element
+	static const unsigned n_base_nodes = 4u;
+
+	/// Constant: The number of all nodes for this element
+	static const unsigned n_all_nodes = 4u;
+
+	/// Constant: The geometric type of the element
+	static const MeshElemType mesh_elem_type = MeshElemType::QUAD;
+
+	/// Constant: The FEM type of the element
+	static const CellType cell_type = CellType::QUAD4;
+
+	/// Constant: The number of faces
+	static const unsigned n_faces = 0;
+
+	/// Constant: The number of edges
+	static const unsigned n_edges = 4;
+
+	/// Constant: The number of neighbors
+	static const unsigned n_neighbors = 4;
+
+	/// Constant: Local node index table for edge
+	static const unsigned edge_nodes[4][2];
+
+	/// Returns the i-th edge of the element.
+	typedef LinearEdgeReturn EdgeReturn;
+
+	/**
+	 * Checks if a point is inside the element.
+	 * @param pnt a 3D GeoLib::Point object
+	 * @param eps tolerance for numerical algorithm used or computing the property
+	 * @return true if the point is not outside the element, false otherwise
+	 */
+	static bool isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps);
+
+	/**
+	 * Tests if the element is geometrically valid.
+	 * @param check_zero_volume indicates if volume == 0 should be checked
+	 */
+	static ElementErrorCode validate(const Element* e);
+
+	/// Returns the ID of a face given an array of nodes.
+	static unsigned identifyFace(Node const* const*, Node* nodes[3]);
+
+	/// Calculates the volume of a convex hexahedron by partitioning it into six tetrahedra.
+	static double computeVolume(Node const* const* _nodes);
+
+}; /* class */
+
+} /* namespace */
+
+#endif /* QUADRULE4_H_ */
+
diff --git a/MeshLib/Elements/QuadRule8.cpp b/MeshLib/Elements/QuadRule8.cpp
new file mode 100644
index 00000000000..7f4e8319d4c
--- /dev/null
+++ b/MeshLib/Elements/QuadRule8.cpp
@@ -0,0 +1,24 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "QuadRule8.h"
+
+namespace MeshLib {
+
+const unsigned QuadRule8::n_all_nodes;
+
+const unsigned QuadRule8::edge_nodes[4][3] =
+{
+		{0, 1, 4}, // Edge 0
+		{1, 2, 5}, // Edge 1
+		{2, 3, 6}, // Edge 2
+		{0, 3, 7}  // Edge 3
+};
+
+} // end namespace MeshLib
diff --git a/MeshLib/Elements/QuadRule8.h b/MeshLib/Elements/QuadRule8.h
new file mode 100644
index 00000000000..f45a44ebd2d
--- /dev/null
+++ b/MeshLib/Elements/QuadRule8.h
@@ -0,0 +1,56 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef QUADRULE8_H_
+#define QUADRULE8_H_
+
+#include "MeshLib/MeshEnums.h"
+#include "Element.h"
+#include "QuadRule4.h"
+#include "EdgeReturn.h"
+
+namespace MeshLib
+{
+
+/**
+ * This class represents a 2d quadliteral element. The following sketch shows the node and edge numbering.
+ * @anchor QuadNodeAndEdgeNumbering
+ * @code
+ *              2
+ *        3-----------2
+ *        |           |
+ *        |           |
+ *       3|           |1
+ *        |           |
+ *        |           |
+ *        0-----------1
+ *              0
+ * @endcode
+ */
+class QuadRule8 : public QuadRule4
+{
+public:
+	/// Constant: The number of all nodes for this element
+	static const unsigned n_all_nodes = 8u;
+
+	/// Constant: The FEM type of the element
+	static const CellType cell_type = CellType::QUAD8;
+
+	/// Constant: Local node index table for edge
+	static const unsigned edge_nodes[4][3];
+
+	/// Returns the i-th edge of the element.
+	typedef QuadraticEdgeReturn EdgeReturn;
+
+}; /* class */
+
+} /* namespace */
+
+#endif /* QUADRULE8_H_ */
+
diff --git a/MeshLib/Elements/QuadRule9.cpp b/MeshLib/Elements/QuadRule9.cpp
new file mode 100644
index 00000000000..7b999ed72f8
--- /dev/null
+++ b/MeshLib/Elements/QuadRule9.cpp
@@ -0,0 +1,16 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "QuadRule9.h"
+
+namespace MeshLib {
+
+const unsigned QuadRule9::n_all_nodes;
+
+} // end namespace MeshLib
diff --git a/MeshLib/Elements/QuadRule9.h b/MeshLib/Elements/QuadRule9.h
new file mode 100644
index 00000000000..5cac57cb571
--- /dev/null
+++ b/MeshLib/Elements/QuadRule9.h
@@ -0,0 +1,48 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef QUADRULE9_H_
+#define QUADRULE9_H_
+
+#include "MeshLib/MeshEnums.h"
+#include "Element.h"
+#include "QuadRule8.h"
+
+namespace MeshLib
+{
+
+/**
+ * This class represents a 2d quadliteral element. The following sketch shows the node and edge numbering.
+ * @anchor QuadNodeAndEdgeNumbering
+ * @code
+ *              2
+ *        3-----------2
+ *        |           |
+ *        |           |
+ *       3|           |1
+ *        |           |
+ *        |           |
+ *        0-----------1
+ *              0
+ * @endcode
+ */
+class QuadRule9 : public QuadRule8
+{
+public:
+	/// Constant: The number of all nodes for this element
+	static const unsigned n_all_nodes = 9u;
+
+	/// Constant: The FEM type of the element
+	static const CellType cell_type = CellType::QUAD9;
+}; /* class */
+
+} /* namespace */
+
+#endif /* QUADRULE9_H_ */
+
diff --git a/MeshLib/Elements/TemplateElement-impl.h b/MeshLib/Elements/TemplateElement-impl.h
new file mode 100644
index 00000000000..58007932e54
--- /dev/null
+++ b/MeshLib/Elements/TemplateElement-impl.h
@@ -0,0 +1,71 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <algorithm>
+
+namespace MeshLib
+{
+
+template <class T_BASE, class ELEMENT_RULE>
+const unsigned TemplateElement<T_BASE, ELEMENT_RULE>::n_all_nodes;
+
+template <class T_BASE, class ELEMENT_RULE>
+const unsigned TemplateElement<T_BASE, ELEMENT_RULE>::n_base_nodes;
+
+template <class T_BASE, class ELEMENT_RULE>
+const unsigned TemplateElement<T_BASE, ELEMENT_RULE>::dimension;
+
+template <class T_BASE, class ELEMENT_RULE>
+TemplateElement<T_BASE, ELEMENT_RULE>::TemplateElement(Node* nodes[n_all_nodes], unsigned value, std::size_t id)
+: T_BASE(value, 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 T_BASE, class ELEMENT_RULE>
+TemplateElement<T_BASE, ELEMENT_RULE>::TemplateElement(std::array<Node*, n_all_nodes> const& nodes, unsigned value, std::size_t id)
+: T_BASE(value, 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 T_BASE, class ELEMENT_RULE>
+TemplateElement<T_BASE, ELEMENT_RULE>::TemplateElement(const TemplateElement &e)
+: T_BASE(e.getValue(), e.getID())
+{
+	this->_nodes = new Node*[n_all_nodes];
+	for (unsigned i=0; i<n_all_nodes; i++)
+		this->_nodes[i] = e._nodes[i];
+	this->_neighbors = new Element*[getNNeighbors()];
+	for (unsigned i=0; i<getNNeighbors(); i++)
+		this->_neighbors[i] = e._neighbors[i];
+	this->_content = e.getContent();
+}
+
+template <class T_BASE, class ELEMENT_RULE>
+bool TemplateElement<T_BASE, ELEMENT_RULE>::isEdge(unsigned idx1, unsigned idx2) const
+{
+	for (unsigned i(0); i<getNEdges(); i++)
+	{
+		if (ELEMENT_RULE::edge_nodes[i][0]==idx1 && ELEMENT_RULE::edge_nodes[i][1]==idx2) return true;
+		if (ELEMENT_RULE::edge_nodes[i][1]==idx1 && ELEMENT_RULE::edge_nodes[i][0]==idx2) return true;
+	}
+	return false;
+}
+
+
+} // MeshLib
+
diff --git a/MeshLib/Elements/TemplateElement.h b/MeshLib/Elements/TemplateElement.h
new file mode 100644
index 00000000000..ed16bf2ab98
--- /dev/null
+++ b/MeshLib/Elements/TemplateElement.h
@@ -0,0 +1,155 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef TEMPLATEELEMENT_H_
+#define TEMPLATEELEMENT_H_
+
+#include <array>
+#include <limits>
+
+#include "MathLib/Point3d.h"
+
+#include "MeshLib/Node.h"
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/MeshQuality/ElementErrorCode.h"
+
+namespace MeshLib
+{
+
+/**
+ * Template for implementing mesh element classes
+ *
+ * \tparam T_BASE         Base element class, e.g. Face, Cell
+ * \tparam ELEMENT_RULE   Geometrical and topological rules of the element
+ */
+template <class T_BASE, class ELEMENT_RULE>
+class TemplateElement : public T_BASE
+{
+public:
+	/// Constant: The number of all nodes for this element
+	static const unsigned n_all_nodes = ELEMENT_RULE::n_all_nodes;
+
+	/// Constant: The number of base nodes for this element
+	static const unsigned n_base_nodes = ELEMENT_RULE::n_base_nodes;
+
+	/// Constant: The dimension of this element
+	static const unsigned dimension = ELEMENT_RULE::dimension;
+
+	/**
+	 * Constructor with an array of mesh nodes.
+	 *
+	 * @param nodes  an array of pointers of mesh nodes which form this element
+	 * @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());
+
+	/**
+	 * Constructor with an array of mesh nodes
+	 *
+	 * @param nodes  an array of pointers of mesh nodes which form this element
+	 * @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());
+
+	/// Copy constructor
+	TemplateElement(const TemplateElement &e);
+
+	/// Destructor
+	virtual ~TemplateElement() {}
+
+	/// Returns a copy of this object.
+	virtual Element* clone() const
+	{
+		return new TemplateElement(*this);
+	}
+
+	/// Get dimension of the mesh element.
+	unsigned getDimension() const { return dimension; }
+
+	/// Returns the edge i of the element.
+	const Element* getEdge(unsigned i) const { return ELEMENT_RULE::EdgeReturn::getEdge(this, i); }
+
+	/// Returns the face i of the element.
+	const Element* getFace(unsigned i) const { return ELEMENT_RULE::getFace(this, i); }
+
+	/// Get the number of edges for this element.
+	unsigned getNEdges() const { return ELEMENT_RULE::n_edges; }
+
+	/// Get the number of nodes for face i.
+	unsigned getNFaceNodes(unsigned i) const { return ELEMENT_RULE::getNFaceNodes(i); }
+
+	/// Get the number of faces for this element.
+	unsigned getNFaces() const { return ELEMENT_RULE::n_faces; }
+
+	/// Get the number of neighbors for this element.
+	unsigned getNNeighbors() const { return ELEMENT_RULE::n_neighbors; }
+
+	/// Get the number of linear nodes for this element.
+	virtual unsigned getNBaseNodes() const { return n_base_nodes; }
+
+	/// Get the number of all nodes for this element.
+	virtual unsigned getNNodes() const { return n_all_nodes; }
+
+	/// Get the type of this element.
+	virtual MeshElemType getGeomType() const { return ELEMENT_RULE::mesh_elem_type; }
+
+	/// Get the FEM type of this element.
+	virtual CellType getCellType() const { return ELEMENT_RULE::cell_type; }
+
+	/// Returns true if these two indices form an edge and false otherwise
+	bool isEdge(unsigned idx1, unsigned idx2) const;
+
+	/**
+	 * Checks if a point is inside the element.
+	 * @param pnt a 3D GeoLib::Point object
+	 * @param eps tolerance for numerical algorithm used or computing the property
+	 * @return true if the point is not outside the element, false otherwise
+	 */
+	bool isPntInElement(MathLib::Point3d const& pnt, double eps = std::numeric_limits<double>::epsilon()) const
+	{
+		return ELEMENT_RULE::isPntInElement(this->_nodes, pnt, eps);
+	}
+
+	/**
+	 * Tests if the element is geometrically valid.
+	 * @param check_zero_volume indicates if volume == 0 should be checked
+	 */
+	virtual ElementErrorCode validate() const
+	{
+		return ELEMENT_RULE::validate(this);
+	}
+
+	/// Returns the ID of a face given an array of nodes.
+	unsigned identifyFace(Node* nodes[3]) const
+	{
+		return ELEMENT_RULE::identifyFace(this->_nodes, nodes);
+	};
+
+	/// Calculates the volume of a convex hexahedron by partitioning it into six tetrahedra.
+	virtual double computeVolume() {return ELEMENT_RULE::computeVolume(this->_nodes);}
+
+	/// Return a specific edge node.
+	virtual inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const
+	{
+		if (getNEdges()>0)
+			return const_cast<Node*>(this->_nodes[ELEMENT_RULE::edge_nodes[edge_id][node_id]]);
+		else
+			return nullptr;
+	}
+
+};
+
+} // MeshLib
+
+#include "TemplateElement-impl.h"
+
+#endif /* TEMPLATEELEMENT_H_ */
+
diff --git a/MeshLib/Elements/TemplateHex-impl.h b/MeshLib/Elements/TemplateHex-impl.h
deleted file mode 100644
index cb1be254d82..00000000000
--- a/MeshLib/Elements/TemplateHex-impl.h
+++ /dev/null
@@ -1,195 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date   2012-05-02
- * \brief  Implementation of the TemplateHex class.
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include "logog/include/logog.hpp"
-
-#include "GeoLib/AnalyticalGeometry.h"
-#include "MeshLib/Node.h"
-
-#include "Quad.h"
-#include "Prism.h"
-
-
-namespace MeshLib {
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-const unsigned TemplateHex<NNODES, CELLHEXTYPE>::n_all_nodes;
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-const unsigned TemplateHex<NNODES, CELLHEXTYPE>::n_base_nodes;
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-const unsigned TemplateHex<NNODES,CELLHEXTYPE>::_face_nodes[6][4] =
-{
-	{0, 3, 2, 1}, // Face 0
-	{0, 1, 5, 4}, // Face 1
-	{1, 2, 6, 5}, // Face 2
-	{2, 3, 7, 6}, // Face 3
-	{3, 0, 4, 7}, // Face 4
-	{4, 5, 6, 7}  // Face 5
-};
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-const unsigned TemplateHex<NNODES,CELLHEXTYPE>::_edge_nodes[12][2] =
-{
-	{0, 1}, // Edge 0
-	{1, 2}, // Edge 1
-	{2, 3}, // Edge 2
-	{0, 3}, // Edge 3
-	{4, 5}, // Edge 4
-	{5, 6}, // Edge 5
-	{6, 7}, // Edge 6
-	{4, 7}, // Edge 7
-	{0, 4}, // Edge 8
-	{1, 5}, // Edge 9
-	{2, 6}, // Edge 10
-	{3, 7}  // Edge 11
-};
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-TemplateHex<NNODES,CELLHEXTYPE>::TemplateHex(Node* nodes[NNODES], unsigned value, std::size_t id)
-	: Cell(value, id)
-{
-	_nodes = nodes;
-
-	_neighbors = new Element*[6];
-	std::fill(_neighbors, _neighbors + 6, nullptr);
-
-	this->_volume = this->computeVolume();
-}
-
-template<unsigned NNODES, CellType CELLHEXTYPE>
-TemplateHex<NNODES,CELLHEXTYPE>::TemplateHex(std::array<Node*, NNODES> const& nodes,
-                                             unsigned value, std::size_t id)
-	: Cell(value, id)
-{
-	_nodes = new Node*[NNODES];
-	std::copy(nodes.begin(), nodes.end(), _nodes);
-
-	_neighbors = new Element*[6];
-	std::fill(_neighbors, _neighbors + 6, nullptr);
-
-	this->_volume = this->computeVolume();
-}
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-TemplateHex<NNODES,CELLHEXTYPE>::TemplateHex(const TemplateHex<NNODES,CELLHEXTYPE> &hex)
-	: Cell(hex.getValue(), hex.getID())
-{
-	_nodes = new Node*[NNODES];
-	for (unsigned i=0; i<NNODES; i++)
-		_nodes[i] = hex._nodes[i];
-
-	_neighbors = new Element*[6];
-	for (unsigned i=0; i<6; i++)
-		_neighbors[i] = hex._neighbors[i];
-
-	_volume = hex.getVolume();
-}
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-TemplateHex<NNODES,CELLHEXTYPE>::~TemplateHex()
-{
-}
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-double TemplateHex<NNODES,CELLHEXTYPE>::computeVolume()
-{
-	return GeoLib::calcTetrahedronVolume(_nodes[4]->getCoords(), _nodes[7]->getCoords(), _nodes[5]->getCoords(), _nodes[0]->getCoords())
-		 + GeoLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[3]->getCoords(), _nodes[1]->getCoords(), _nodes[0]->getCoords())
-		 + GeoLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[7]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords())
-		 + GeoLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[7]->getCoords(), _nodes[6]->getCoords(), _nodes[2]->getCoords())
-		 + GeoLib::calcTetrahedronVolume(_nodes[1]->getCoords(), _nodes[3]->getCoords(), _nodes[5]->getCoords(), _nodes[2]->getCoords())
-		 + GeoLib::calcTetrahedronVolume(_nodes[3]->getCoords(), _nodes[7]->getCoords(), _nodes[5]->getCoords(), _nodes[2]->getCoords());
-}
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-const Element* TemplateHex<NNODES,CELLHEXTYPE>::getFace(unsigned i) const
-{
-	if (i<this->getNFaces())
-	{
-		unsigned nFaceNodes (this->getNFaceNodes(i));
-		Node** nodes = new Node*[nFaceNodes];
-		for (unsigned j=0; j<nFaceNodes; j++)
-			nodes[j] = _nodes[_face_nodes[i][j]];
-		return new Quad(nodes);
-	}
-	ERR("Error in MeshLib::Element::getFace() - Index %d does not exist.", i);
-	return NULL;
-}
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-bool TemplateHex<NNODES,CELLHEXTYPE>::isEdge(unsigned idx1, unsigned idx2) const
-{
-	for (unsigned i(0); i<12; i++)
-	{
-		if (_edge_nodes[i][0]==idx1 && _edge_nodes[i][1]==idx2) return true;
-		if (_edge_nodes[i][1]==idx1 && _edge_nodes[i][0]==idx2) return true;
-	}
-	return false;
-}
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-bool TemplateHex<NNODES,CELLHEXTYPE>::isPntInElement(MathLib::Point3d const& pnt, double eps) const
-{
-    return (GeoLib::isPointInTetrahedron(pnt, *_nodes[4], *_nodes[7], *_nodes[5], *_nodes[0], eps) ||
-            GeoLib::isPointInTetrahedron(pnt, *_nodes[5], *_nodes[3], *_nodes[1], *_nodes[0], eps) ||
-            GeoLib::isPointInTetrahedron(pnt, *_nodes[5], *_nodes[7], *_nodes[3], *_nodes[0], eps) ||
-            GeoLib::isPointInTetrahedron(pnt, *_nodes[5], *_nodes[7], *_nodes[6], *_nodes[2], eps) ||
-            GeoLib::isPointInTetrahedron(pnt, *_nodes[1], *_nodes[3], *_nodes[5], *_nodes[2], eps) ||
-            GeoLib::isPointInTetrahedron(pnt, *_nodes[3], *_nodes[7], *_nodes[5], *_nodes[2], eps));
-}
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-Element* TemplateHex<NNODES,CELLHEXTYPE>::clone() const
-{
-	return new TemplateHex<NNODES,CELLHEXTYPE>(*this);
-}
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-unsigned TemplateHex<NNODES,CELLHEXTYPE>::identifyFace(Node* nodes[3]) const
-{
-	for (unsigned i=0; i<6; i++)
-	{
-		unsigned flag(0);
-		for (unsigned j=0; j<4; j++)
-			for (unsigned k=0; k<3; k++)
-				if (_nodes[_face_nodes[i][j]] == nodes[k])
-					flag++;
-		if (flag==3)
-			return i;
-	}
-	return std::numeric_limits<unsigned>::max();
-}
-
-template <unsigned NNODES, CellType CELLHEXTYPE>
-ElementErrorCode TemplateHex<NNODES,CELLHEXTYPE>::validate() const
-{
-	ElementErrorCode error_code;
-	error_code[ElementErrorFlag::ZeroVolume] = this->hasZeroVolume();
-
-	for (unsigned i=0; i<6; ++i)
-	{
-		if (error_code.all())
-			break;
-
-		const MeshLib::Element* quad (this->getFace(i));
-		error_code |= quad->validate();
-		delete quad;
-	}
-	error_code[ElementErrorFlag::NodeOrder]  = !this->testElementNodeOrder();
-	return error_code;
-}
-
-} // end namespace MeshLib
diff --git a/MeshLib/Elements/TemplateHex.h b/MeshLib/Elements/TemplateHex.h
deleted file mode 100644
index e781bea7154..00000000000
--- a/MeshLib/Elements/TemplateHex.h
+++ /dev/null
@@ -1,156 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date   2012-05-02
- * \brief  Definition of the TemplateHex class.
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef TEMPLATEHEX_H_
-#define TEMPLATEHEX_H_
-
-#include <array>
-#include "MeshLib/MeshEnums.h"
-#include "Cell.h"
-
-namespace MeshLib {
-
-/**
- * A 3d Hexahedron Element.
- * @code
- *
- *  Hex:
- *                6
- *          7-----------6
- *         /:          /|
- *        / :         / |
- *      7/  :        /5 |
- *      / 11:       /   | 10
- *     /    : 4    /    |
- *    4-----------5     |
- *    |     :     | 2   |
- *    |     3.....|.....2
- *    |    .      |    /
- *  8 |   .       |9  /
- *    | 3.        |  / 1
- *    | .         | /
- *    |.          |/
- *    0-----------1
- *          0
- *
- * @endcode
- */
-template <unsigned NNODES, CellType CELLHEXTYPE>
-class TemplateHex : public Cell
-{
-public:
-	/// Constant: The number of all nodes for this element
-	static const unsigned n_all_nodes = NNODES;
-
-	/// Constant: The number of base nodes for this element
-	static const unsigned n_base_nodes = 8u;
-
-	/// Constructor with an array of mesh nodes.
-	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, std::size_t id = std::numeric_limits<std::size_t>::max());
-
-	/// Copy constructor
-	TemplateHex(const TemplateHex &hex);
-
-	/// Destructor
-	virtual ~TemplateHex();
-
-	/// Returns the face i of the element.
-	const Element* getFace(unsigned i) const;
-
-	/// Get the number of edges for this element.
-	unsigned getNEdges() const { return 12; }
-
-	/// Get the number of nodes for face i.
-	unsigned getNFaceNodes(unsigned i) const { (void)i; return 4; }
-
-	/// Get the number of faces for this element.
-	unsigned getNFaces() const { return 6; }
-
-	/// Get the number of neighbors for this element.
-	unsigned getNNeighbors() const { return 6; }
-
-	/// Get the number of linear nodes for this element.
-	virtual unsigned getNBaseNodes() const
-	{
-		return n_base_nodes;
-	}
-
-	/// Get the number of all nodes for this element.
-	virtual unsigned getNNodes() const
-	{
-		return n_all_nodes;
-	}
-
-	/**
-	 * Method returns the type of the element. In this case HEXAHEDRON will be returned.
-	 * @return MeshElemType::HEXAHEDRON
-	 */
-	virtual MeshElemType getGeomType() const { return MeshElemType::HEXAHEDRON; }
-
-	/**
-	 * Method returns the FEM type of the element.
-	 * @return
-	 */
-	virtual CellType getCellType() const { return CELLHEXTYPE; }
-
-	/// Returns true if these two indices form an edge and false otherwise
-	bool isEdge(unsigned i, unsigned j) const;
-
-    /**
-	 * Checks if a point is inside the element.
-	 * @param pnt a 3D MathLib::Point3d object
-	 * @param eps tolerance for numerical algorithm used or computing the property
-	 * @return true if the point is not outside the element, false otherwise
-	 */
-    bool isPntInElement(MathLib::Point3d const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
-
-	/**
-	 * Tests if the element is geometrically valid.
-	 * @param check_zero_volume indicates if volume == 0 should be checked
-	 */
-	virtual ElementErrorCode validate() const;
-
-	/**
-	 * Method clone is inherited from class Element. It makes a deep copy of the Hex instance.
-	 * @return an exact copy of the object
-	 */
-	virtual Element* clone() const;
-
-	/// Returns the ID of a face given an array of nodes.
-	unsigned identifyFace(Node* nodes[3]) const;
-
-protected:
-	/// Calculates the volume of a convex hexahedron by partitioning it into six tetrahedra.
-	double computeVolume();
-
-	/// Return a specific edge node.
-	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const
-	{
-		return _nodes[_edge_nodes[edge_id][node_id]];
-	}
-
-	static const unsigned _face_nodes[6][4];
-	static const unsigned _edge_nodes[12][2];
-
-}; /* class */
-
-} /* namespace */
-
-#include "TemplateHex-impl.h"
-
-#endif /* TEMPLATEHEX_H_ */
-
diff --git a/MeshLib/Elements/TemplateLine-impl.h b/MeshLib/Elements/TemplateLine-impl.h
deleted file mode 100644
index 01ae92f2848..00000000000
--- a/MeshLib/Elements/TemplateLine-impl.h
+++ /dev/null
@@ -1,85 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date   Sep 27, 2012
- * \brief  Implementation of the TemplateLine class.
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-namespace MeshLib
-{
-template<unsigned NNODES, CellType CELLLINETYPE>
-TemplateLine<NNODES,CELLLINETYPE>::TemplateLine(std::array<Node*, NNODES> const& nodes,
-                                                unsigned value, std::size_t id)
-	: Element(value, id)
-{
-	_nodes = new Node*[NNODES];
-	std::copy(nodes.begin(), nodes.end(), _nodes);
-
-	_neighbors = new Element*[2];
-	std::fill_n(_neighbors, 2, nullptr);
-	this->_length = this->computeVolume();
-}
-
-template<unsigned NNODES, CellType CELLLINETYPE>
-TemplateLine<NNODES,CELLLINETYPE>::TemplateLine(Node* nodes[NNODES], unsigned value, std::size_t id)
-	: Element(value, id)
-{
-	_nodes = nodes;
-	_neighbors = new Element*[2];
-	std::fill_n(_neighbors, 2, nullptr);
-	this->_length = this->computeVolume();
-}
-
-template <unsigned NNODES, CellType CELLLINETYPE>
-TemplateLine<NNODES,CELLLINETYPE>::TemplateLine(const TemplateLine<NNODES,CELLLINETYPE> &line)
-	: Element(line.getValue(), line.getID())
-{
-	_nodes = new Node*[NNODES];
-	for (unsigned k(0); k<NNODES; k++)
-		_nodes[k] = line._nodes[k];
-
-	_neighbors = new Element*[2];
-	_neighbors[0] = line._neighbors[0];
-	_neighbors[1] = line._neighbors[1];
-	_length = line.getLength();
-}
-
-template <unsigned NNODES, CellType CELLLINETYPE>
-TemplateLine<NNODES,CELLLINETYPE>::~TemplateLine()
-{}
-
-template <unsigned NNODES, CellType CELLLINETYPE>
-bool TemplateLine<NNODES,CELLLINETYPE>::isPntInElement(MathLib::Point3d const& pnt, double eps) const
-{
-    double tmp;
-    double tmp_dst(0);
-    double const dist =MathLib::calcProjPntToLineAndDists(pnt.getCoords(), _nodes[0]->getCoords(), _nodes[1]->getCoords(), tmp, tmp_dst);
-    return (dist < eps);
-}
-
-template <unsigned NNODES, CellType CELLLINETYPE>
-ElementErrorCode TemplateLine<NNODES,CELLLINETYPE>::validate() const
-{ 
-	ElementErrorCode error_code;
-	error_code[ElementErrorFlag::ZeroVolume] = this->hasZeroVolume();
-	return error_code;
-}
-
-template <unsigned NNODES, CellType CELLLINETYPE>
-const unsigned TemplateLine<NNODES, CELLLINETYPE>::dimension;
-
-template <unsigned NNODES, CellType CELLLINETYPE>
-const unsigned TemplateLine<NNODES, CELLLINETYPE>::n_all_nodes;
-
-template <unsigned NNODES, CellType CELLLINETYPE>
-const unsigned TemplateLine<NNODES, CELLLINETYPE>::n_base_nodes;
-
-} // namespace MeshLib
-
diff --git a/MeshLib/Elements/TemplateLine.h b/MeshLib/Elements/TemplateLine.h
deleted file mode 100644
index 514e8b3c99a..00000000000
--- a/MeshLib/Elements/TemplateLine.h
+++ /dev/null
@@ -1,188 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date   2012-05-02
- * \brief  Definition of the Line class.
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef TEMPLATELINE_H_
-#define TEMPLATELINE_H_
-
-#include <array>
-#include <limits>
-#include <cmath>
-
-#include "MathLib/MathTools.h"
-
-#include "MeshLib/MeshEnums.h"
-#include "MeshLib/Node.h"
-
-#include "Element.h"
-
-namespace MeshLib {
-
-/**
- * A 1d Edge or Line Element.
- * @code
- *  0--------1
- * @endcode
- */
-template<unsigned NNODES, CellType CELLLINETYPE>
-class TemplateLine : public Element
-{
-public:
-	/// Constant: Dimension of this mesh element
-	static const unsigned dimension = 1u;
-
-	/// Constant: The number of all nodes for this element
-	static const unsigned n_all_nodes = NNODES;
-
-	/// Constant: The number of base nodes for this element
-	static const unsigned n_base_nodes = 2u;
-
-	/// Constructor with an array of mesh nodes.
-	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, std::size_t id = std::numeric_limits<std::size_t>::max());
-
-	/// Copy constructor
-	TemplateLine(const TemplateLine &line);
-
-	/// Destructor
-	virtual ~TemplateLine();
-
-	/// Compute the minimum and maximum squared edge length for this element
-	void computeSqrEdgeLengthRange(double &min, double &max) const { min = _length; max = _length; }
-
-	/// Returns the length, area or volume of a 1D, 2D or 3D element
-	double getContent() const { return _length; }
-
-	/// Get dimension of the mesh element.
-	unsigned getDimension() const { return dimension; }
-
-	/// Returns the edge i of the element.
-	const Element* getEdge(unsigned /*i*/) const { return nullptr; }
-
-	/// Returns the face i of the element.
-	const Element* getFace(unsigned /*i*/) const { return nullptr; }
-
-	/// Get the length of this 1d element.
-	double getLength() const { return _length; }
-
-	/// 1D elements have no edges
-	unsigned getNEdges() const { return 1; }
-
-	/// Get the number of nodes for face i.
-	unsigned getNFaceNodes(unsigned /*i*/) const { return 0; }
-
-	/// Get the number of faces for this element.
-	unsigned getNFaces() const { return 0; }
-
-	/// Get the number of neighbors for this element.
-	unsigned getNNeighbors() const { return 2; }
-
-	/// Get the number of linear nodes for this element.
-	virtual unsigned getNBaseNodes() const
-	{
-		return n_base_nodes;
-	}
-
-	/// Get the number of all nodes for this element.
-	virtual unsigned getNNodes() const
-	{
-		return n_all_nodes;
-	}
-
-	/**
-		* Method returns the type of the element. In this case LINE will be returned.
-		* @return MeshElemType::LINE
-		*/
-	virtual MeshElemType getGeomType() const { return MeshElemType::LINE; }
-
-	/**
-		* Get the type of the element in context of the finite element method.
-		* @return a value of the enum FEMElemType::type
-		*/
-	virtual CellType getCellType() const { return CELLLINETYPE; }
-
-	/// Returns true if these two indices form an edge and false otherwise
-	bool isEdge(unsigned idx1, unsigned idx2) const
-	{
-		if (0==idx1 && 1==idx2) return true;
-		if (1==idx1 && 0==idx2) return true;
-		return false;
-	}
-
-	/**
-	 * Checks if a point is located on the line
-	 * @param pnt a 3D MathLib::Point3d object
-	 * @param eps tolerance for numerical algorithm used or computing the property
-	 * @return true if the point is located on the line, false otherwise
-	 */
-	bool isPntInElement(MathLib::Point3d const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
-
-	/**
-		* Tests if the element is geometrically valid.
-		* @return error code (0 = okay, 1 = zero volume)
-		*/
-	virtual ElementErrorCode validate() const;
-
-	/**
-		* Checks if the node order of an element is correct by testing surface normals.
-		* For 1D elements this always returns true.
-		*/
-	virtual bool testElementNodeOrder() const { return true; }
-
-	/**
-		* Method clone is inherited from class Element. It makes a deep copy of the TemplateLine instance.
-		* @return an exact copy of the object
-		*/
-	virtual Element* clone() const
-	{
-		return new TemplateLine<NNODES,CELLLINETYPE>(*this);
-	}
-
-protected:
-	double computeVolume()
-	{
-		return sqrt(MathLib::sqrDist(_nodes[0]->getCoords(), _nodes[1]->getCoords()));
-	}
-
-	/// Returns the specified node.
-	Node* getEdgeNode(unsigned edge_id, unsigned node_id) const
-	{
-		if (edge_id==0 && node_id<2)
-			return _nodes[node_id];
-		return nullptr;
-	}
-
-	/// Returns the ID of a face given an array of nodes.
-	/// As element faces are always element->getDimensionality() - 1, the "face" of a line is just a node
-	/// and the method returns if another line is connected to node[0] or node[1] of the line.
-	unsigned identifyFace(Node* nodes[1]) const
-	{
-		if (nodes[0] == _nodes[0])
-			return 0;
-		if (nodes[0] == _nodes[1])
-			return 1;
-		return std::numeric_limits<unsigned>::max();
-	}
-
-	double _length;
-
-}; /* class */
-
-} /* namespace */
-
-#include "TemplateLine-impl.h"
-
-#endif /* TEMPLATELINE_H_ */
-
diff --git a/MeshLib/Elements/TemplatePrism-impl.h b/MeshLib/Elements/TemplatePrism-impl.h
deleted file mode 100644
index 1393be91522..00000000000
--- a/MeshLib/Elements/TemplatePrism-impl.h
+++ /dev/null
@@ -1,200 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date   2012-05-02
- * \brief  Implementation of the TemplatePrism class.
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-// Thirdparty
-#include "logog/include/logog.hpp"
-
-#include "GeoLib/AnalyticalGeometry.h"
-#include "MeshLib/Node.h"
-
-#include "Tri.h"
-#include "Quad.h"
-
-namespace MeshLib {
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-const unsigned TemplatePrism<NNODES, CELLPRISMTYPE>::n_all_nodes;
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-const unsigned TemplatePrism<NNODES, CELLPRISMTYPE>::n_base_nodes;
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-const unsigned TemplatePrism<NNODES,CELLPRISMTYPE>::_face_nodes[5][4] =
-{
-	{0, 2, 1, 99}, // Face 0
-	{0, 1, 4,  3}, // Face 1
-	{1, 2, 5,  4}, // Face 2
-	{2, 0, 3,  5}, // Face 3
-	{3, 4, 5, 99}  // Face 4
-};
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-const unsigned TemplatePrism<NNODES,CELLPRISMTYPE>::_edge_nodes[9][2] =
-{
-	{0, 1}, // Edge 0
-	{1, 2}, // Edge 1
-	{0, 2}, // Edge 2
-	{0, 3}, // Edge 3
-	{1, 4}, // Edge 4
-	{2, 5}, // Edge 5
-	{3, 4}, // Edge 6
-	{4, 5}, // Edge 7
-	{3, 5}  // Edge 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, std::size_t id)
-	: Cell(value, id)
-{
-	_nodes = nodes;
-	_neighbors = new Element*[5];
-	std::fill(_neighbors, _neighbors + 5, nullptr);
-	this->_volume = this->computeVolume();
-}
-
-template<unsigned NNODES, CellType CELLPRISMTYPE>
-TemplatePrism<NNODES,CELLPRISMTYPE>::TemplatePrism(std::array<Node*, NNODES> const& nodes,
-                                                   unsigned value, std::size_t id)
-	: Cell(value, id)
-{
-	_nodes = new Node*[NNODES];
-	std::copy(nodes.begin(), nodes.end(), _nodes);
-
-	_neighbors = new Element*[5];
-	std::fill(_neighbors, _neighbors + 5, nullptr);
-
-	this->_volume = this->computeVolume();
-}
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-TemplatePrism<NNODES,CELLPRISMTYPE>::TemplatePrism(const TemplatePrism<NNODES,CELLPRISMTYPE> &prism)
-	: Cell(prism.getValue(), prism.getID())
-{
-	_nodes = new Node*[NNODES];
-	for (unsigned i=0; i<NNODES; i++)
-		_nodes[i] = prism._nodes[i];
-
-	_neighbors = new Element*[5];
-	for (unsigned i=0; i<5; i++)
-		_neighbors[i] = prism._neighbors[i];
-
-	_volume = prism.getVolume();
-}
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-TemplatePrism<NNODES,CELLPRISMTYPE>::~TemplatePrism()
-{
-}
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-double TemplatePrism<NNODES,CELLPRISMTYPE>::computeVolume()
-{
-	return GeoLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords())
-		 + GeoLib::calcTetrahedronVolume(_nodes[1]->getCoords(), _nodes[4]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords())
-		 + GeoLib::calcTetrahedronVolume(_nodes[2]->getCoords(), _nodes[4]->getCoords(), _nodes[5]->getCoords(), _nodes[3]->getCoords());
-}
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-const Element* TemplatePrism<NNODES,CELLPRISMTYPE>::getFace(unsigned i) const
-{
-	if (i<this->getNFaces())
-	{
-		unsigned nFaceNodes (this->getNFaceNodes(i));
-		Node** nodes = new Node*[nFaceNodes];
-		for (unsigned j=0; j<nFaceNodes; j++)
-			nodes[j] = _nodes[_face_nodes[i][j]];
-
-		if (i==0 || i==4)
-			return new Tri(nodes);
-		else
-			return new Quad(nodes);
-	}
-	ERR("Error in MeshLib::Element::getFace() - Index %d does not exist.", i);
-	return NULL;
-}
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-unsigned TemplatePrism<NNODES,CELLPRISMTYPE>::getNFaceNodes(unsigned i) const
-{
-	if (i<5)
-		return _n_face_nodes[i];
-	ERR("Error in MeshLib::Element::getNFaceNodes() - Index %d does not exist.", i);
-	return 0;
-}
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-bool TemplatePrism<NNODES,CELLPRISMTYPE>::isEdge(unsigned idx1, unsigned idx2) const
-{
-	for (unsigned i(0); i<9; i++)
-	{
-		if (_edge_nodes[i][0]==idx1 && _edge_nodes[i][1]==idx2) return true;
-		if (_edge_nodes[i][1]==idx1 && _edge_nodes[i][0]==idx2) return true;
-	}
-	return false;
-}
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-bool TemplatePrism<NNODES,CELLPRISMTYPE>::isPntInElement(MathLib::Point3d const& pnt, double eps) const
-{
-    return (GeoLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3], eps) ||
-            GeoLib::isPointInTetrahedron(pnt, *_nodes[1], *_nodes[4], *_nodes[2], *_nodes[3], eps) ||
-            GeoLib::isPointInTetrahedron(pnt, *_nodes[2], *_nodes[4], *_nodes[5], *_nodes[3], eps));
-}
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-Element* TemplatePrism<NNODES,CELLPRISMTYPE>::clone() const
-{
-	return new TemplatePrism<NNODES,CELLPRISMTYPE>(*this);
-}
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-unsigned TemplatePrism<NNODES,CELLPRISMTYPE>::identifyFace(Node* nodes[3]) const
-{
-	for (unsigned i=0; i<5; i++)
-	{
-		unsigned flag(0);
-		for (unsigned j=0; j<4; j++)
-			for (unsigned k=0; k<3; k++)
-				if (_face_nodes[i][j] != 99 && _nodes[_face_nodes[i][j]] == nodes[k])
-					flag++;
-		if (flag==3)
-			return i;
-	}
-	return std::numeric_limits<unsigned>::max();
-}
-
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-ElementErrorCode TemplatePrism<NNODES,CELLPRISMTYPE>::validate() const
-{
-	ElementErrorCode error_code;
-	error_code[ElementErrorFlag::ZeroVolume] = this->hasZeroVolume();
-
-	for (unsigned i=1; i<4; ++i)
-	{
-		const MeshLib::Quad* quad (dynamic_cast<const MeshLib::Quad*>(this->getFace(i)));
-		if (quad)
-			error_code |= quad->validate();
-		else
-			error_code.set(ElementErrorFlag::NodeOrder);
-		delete quad;
-	}
-	error_code[ElementErrorFlag::NodeOrder] = !this->testElementNodeOrder();
-	return error_code;
-}
-
-} // end namespace MeshLib
-
diff --git a/MeshLib/Elements/TemplatePrism.h b/MeshLib/Elements/TemplatePrism.h
deleted file mode 100644
index 5f7aa36066c..00000000000
--- a/MeshLib/Elements/TemplatePrism.h
+++ /dev/null
@@ -1,156 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date   2012-05-02
- * \brief  Definition of the TemplatePrism class.
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef TEMPLATEPRISM_H_
-#define TEMPLATEPRISM_H_
-
-#include <array>
-#include "MeshLib/MeshEnums.h"
-#include "Cell.h"
-
-namespace MeshLib {
-
-/**
- * This class represents a 3d prism element. The following sketch shows the node and edge numbering.
- * @anchor PrismNodeAndEdgeNumbering
- * @code
- *            5
- *           / \
- *          / : \
- *        8/  :  \7
- *        /   :5  \
- *       /    :  6 \
- *      3-----------4
- *      |     :     |
- *      |     2     |
- *      |    . .    |
- *     3|   .   .   |4
- *      | 2.     .1 |
- *      | .       . |
- *      |.         .|
- *      0-----------1
- *            0
- *
- * @endcode
- */
-template <unsigned NNODES, CellType CELLPRISMTYPE>
-class TemplatePrism : public Cell
-{
-public:
-	/// Constant: The number of all nodes for this element
-	static const unsigned n_all_nodes = NNODES;
-
-	/// Constant: The number of base nodes for this element
-	static const unsigned n_base_nodes = 6u;
-
-	/// Constructor with an array of mesh nodes.
-	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, std::size_t id = std::numeric_limits<std::size_t>::max());
-
-	/// Copy constructor
-	TemplatePrism(const TemplatePrism &prism);
-
-	/// Destructor
-	virtual ~TemplatePrism();
-
-	/// Returns the face i of the element.
-	const Element* getFace(unsigned i) const;
-
-	/// Get the number of edges for this element.
-	unsigned getNEdges() const { return 9; }
-
-	/// Get the number of nodes for face i.
-	unsigned getNFaceNodes(unsigned i) const;
-
-	/// Get the number of faces for this element.
-	unsigned getNFaces() const { return 5; }
-
-	/// Get the number of neighbors for this element.
-	unsigned getNNeighbors() const { return 5; }
-
-	/// Get the number of linear nodes for this element.
-	virtual unsigned getNBaseNodes() const
-	{
-		return n_base_nodes;
-	}
-
-	/// Get the number of all nodes for this element.
-	virtual unsigned getNNodes() const
-	{
-		return n_all_nodes;
-	}
-
-	/**
-	 * Method returns the type of the element. In this case PRISM will be returned.
-	 * @return MeshElemType::PRISM
-	 */
-	virtual MeshElemType getGeomType() const { return MeshElemType::PRISM; }
-
-	/**
-	 * Get the type of the element in context of the finite element method.
-	 * @return a value of the enum CellType
-	 */
-	virtual CellType getCellType() const { return CELLPRISMTYPE; }
-
-	/// Returns true if these two indeces form an edge and false otherwise
-	bool isEdge(unsigned i, unsigned j) const;
-
-	/**
-	 * Checks if a point is inside the element.
-	 * @param pnt a 3D MathLib::Point3d object
-	 * @param eps tolerance for numerical algorithm used or computing the property
-	 * @return true if the point is not outside the element, false otherwise
-	 */
-	bool isPntInElement(MathLib::Point3d const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
-
-	/**
-	 * Tests if the element is geometrically valid.
-	 * @param check_zero_volume indicates if volume == 0 should be checked
-	 */
-	virtual ElementErrorCode validate() const;
-
-	/**
-	 * Method clone is inherited from class Element. It makes a deep copy of the
-	 * Hex instance employing the copy constructor of class TemplatePrism.
-	 * @return an exact copy of the object
-	 */
-	virtual Element* clone() const;
-
-	/// Returns the ID of a face given an array of nodes.
-	unsigned identifyFace(Node* nodes[3]) const;
-
-protected:
-	/// Calculates the volume of a prism by subdividing it into three tetrahedra.
-	double computeVolume();
-
-	/// Return a specific edge node.
-	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const
-	{
-		return _nodes[_edge_nodes[edge_id][node_id]];
-	}
-
-	static const unsigned _face_nodes[5][4];
-	static const unsigned _edge_nodes[9][2];
-	static const unsigned _n_face_nodes[5];
-
-}; /* class */
-
-} /* namespace */
-
-#include "TemplatePrism-impl.h"
-
-#endif /* TEMPLATEPRISM_H_ */
-
diff --git a/MeshLib/Elements/TemplatePyramid-impl.h b/MeshLib/Elements/TemplatePyramid-impl.h
deleted file mode 100644
index 9328c2109ac..00000000000
--- a/MeshLib/Elements/TemplatePyramid-impl.h
+++ /dev/null
@@ -1,200 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date   2012-05-02
- * \brief  Implementation of the TemplatePyramid class.
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-// Thirdparty
-#include "logog/include/logog.hpp"
-
-#include "GeoLib/AnalyticalGeometry.h"
-#include "MeshLib/Node.h"
-
-#include "Tri.h"
-#include "Quad.h"
-
-namespace MeshLib {
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-const unsigned TemplatePyramid<NNODES, CELLPYRAMIDTYPE>::n_all_nodes;
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-const unsigned TemplatePyramid<NNODES, CELLPYRAMIDTYPE>::n_base_nodes;
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-const unsigned TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::_face_nodes[5][4] =
-{
-	{0, 1, 4, 99}, // Face 0
-	{1, 2, 4, 99}, // Face 1
-	{2, 3, 4, 99}, // Face 2
-	{3, 0, 4, 99}, // Face 3
-	{0, 3, 2,  1}  // Face 4
-};
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-const unsigned TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::_edge_nodes[8][2] =
-{
-	{0, 1}, // Edge 0
-	{1, 2}, // Edge 1
-	{2, 3}, // Edge 2
-	{0, 3}, // Edge 3
-	{0, 4}, // Edge 4
-	{1, 4}, // Edge 5
-	{2, 4}, // Edge 6
-	{3, 4}  // Edge 7
-};
-
-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, std::size_t id)
-	: Cell(value, id)
-{
-	_nodes = nodes;
-
-	_neighbors = new Element*[5];
-	std::fill(_neighbors, _neighbors + 5, nullptr);
-
-	this->_volume = this->computeVolume();
-}
-
-template<unsigned NNODES, CellType CELLPYRAMIDTYPE>
-TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::TemplatePyramid(std::array<Node*, NNODES> const& nodes,
-                                                         unsigned value, std::size_t id)
-	: Cell(value, id)
-{
-	_nodes = new Node*[NNODES];
-	std::copy(nodes.begin(), nodes.end(), _nodes);
-
-	_neighbors = new Element*[5];
-	std::fill(_neighbors, _neighbors + 5, nullptr);
-
-	this->_volume = this->computeVolume();
-}
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::TemplatePyramid(const TemplatePyramid<NNODES,CELLPYRAMIDTYPE> &pyramid)
-	: Cell(pyramid.getValue(), pyramid.getID())
-{
-	_nodes = new Node*[NNODES];
-	for (unsigned i=0; i<NNODES; i++) {
-		_nodes[i] = pyramid._nodes[i];
-	}
-
-	_neighbors = new Element*[5];
-	for (unsigned i=0; i<5; i++) {
-		_neighbors[i] = pyramid._neighbors[i];
-	}
-
-	_volume = pyramid.getVolume();
-}
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::~TemplatePyramid()
-{
-}
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-double TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::computeVolume()
-{
-	return GeoLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[4]->getCoords())
-		 + GeoLib::calcTetrahedronVolume(_nodes[2]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords(), _nodes[4]->getCoords());
-}
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-const Element* TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::getFace(unsigned i) const
-{
-	if (i<this->getNFaces())
-	{
-		unsigned nFaceNodes (this->getNFaceNodes(i));
-		Node** nodes = new Node*[nFaceNodes];
-		for (unsigned j=0; j<nFaceNodes; j++)
-			nodes[j] = _nodes[_face_nodes[i][j]];
-
-		if (i<4)
-			return new Tri(nodes);
-		else
-			return new Quad(nodes);
-	}
-	ERR("Error in MeshLib::Element::getFace() - Index %d does not exist.", i);
-	return NULL;
-}
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-unsigned TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::getNFaceNodes(unsigned i) const
-{
-	if (i<5)
-		return _n_face_nodes[i];
-	ERR("Error in MeshLib::Element::getNFaceNodes() - Index %d does not exist.", i);
-	return 0;
-}
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-bool TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::isEdge(unsigned idx1, unsigned idx2) const
-{
-	for (unsigned i(0); i<8; i++)
-	{
-		if (_edge_nodes[i][0]==idx1 && _edge_nodes[i][1]==idx2) return true;
-		if (_edge_nodes[i][1]==idx1 && _edge_nodes[i][0]==idx2) return true;
-	}
-	return false;
-}
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-bool TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::isPntInElement(MathLib::Point3d const& pnt, double eps) const
-{
-    return (GeoLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[4], eps) ||
-            GeoLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[2], *_nodes[3], *_nodes[4], eps));
-}
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-Element* TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::clone() const
-{
-	return new TemplatePyramid<NNODES,CELLPYRAMIDTYPE>(*this);
-}
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-unsigned TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::identifyFace(Node* nodes[3]) const
-{
-	for (unsigned i=0; i<5; i++)
-	{
-		unsigned flag(0);
-		for (unsigned j=0; j<4; j++)
-			for (unsigned k=0; k<3; k++)
-				if (_face_nodes[i][j] != 99 && _nodes[_face_nodes[i][j]] == nodes[k])
-					flag++;
-		if (flag==3)
-			return i;
-	}
-	return std::numeric_limits<unsigned>::max();
-}
-
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-ElementErrorCode TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::validate() const
-{
-	ElementErrorCode error_code;
-	error_code[ElementErrorFlag::ZeroVolume] = this->hasZeroVolume();
-
-	const MeshLib::Quad* base (dynamic_cast<const MeshLib::Quad*>(this->getFace(4)));
-	if (base)
-	{
-		error_code |= base->validate();
-		error_code[ElementErrorFlag::NodeOrder] = !this->testElementNodeOrder();
-	}
-	else
-		error_code.set(ElementErrorFlag::NodeOrder);
-	delete base;
-
-	return error_code;
-}
-
-} // end namespace MeshLib
diff --git a/MeshLib/Elements/TemplatePyramid.h b/MeshLib/Elements/TemplatePyramid.h
deleted file mode 100644
index 8c874e7daae..00000000000
--- a/MeshLib/Elements/TemplatePyramid.h
+++ /dev/null
@@ -1,154 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date   2012-05-02
- * \brief  Definition of the TemplatePyramid class.
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef TEMPLATEPYRAMID_H_
-#define TEMPLATEPYRAMID_H_
-
-#include <array>
-#include "MeshLib/MeshEnums.h"
-#include "Cell.h"
-
-namespace MeshLib {
-
-/**
- * This class represents a 3d pyramid element. The following sketch shows the node and edge numbering.
- * @anchor PyramidNodeAndEdgeNumbering
- * @code
- *
- *               4
- *             //|\
- *            // | \
- *          7//  |  \6
- *          //   |5  \
- *         //    |    \
- *        3/.... |.....2
- *       ./      |  2 /
- *      ./4      |   /
- *    3./        |  /1
- *    ./         | /
- *   ./          |/
- *  0------------1
- *        0
- * @endcode
- */
-template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
-class TemplatePyramid : public Cell
-{
-public:
-	/// Constant: The number of all nodes for this element
-	static const unsigned n_all_nodes = NNODES;
-
-	/// Constant: The number of base nodes for this element
-	static const unsigned n_base_nodes = 5u;
-
-	/// Constructor with an array of mesh nodes.
-	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, std::size_t id = std::numeric_limits<std::size_t>::max());
-
-	/// Copy constructor
-	TemplatePyramid(const TemplatePyramid &pyramid);
-
-	/// Destructor
-	virtual ~TemplatePyramid();
-
-	/// Returns the face i of the element.
-	const Element* getFace(unsigned i) const;
-
-	/// Get the number of edges for this element.
-	unsigned getNEdges() const { return 8; }
-
-	/// Get the number of nodes for face i.
-	unsigned getNFaceNodes(unsigned i) const;
-
-	/// Get the number of faces for this element.
-	unsigned getNFaces() const { return 5; }
-
-	/// Get the number of neighbors for this element.
-	unsigned getNNeighbors() const { return 5; }
-
-	/// Get the number of linear nodes for this element.
-	virtual unsigned getNBaseNodes() const
-	{
-		return n_base_nodes;
-	}
-
-	/// Get the number of all nodes for this element.
-	virtual unsigned getNNodes() const
-	{
-		return n_all_nodes;
-	}
-
-	/**
-	 * Method returns the type of the element. In this case PYRAMID will be returned.
-	 * @return MeshElemType::PYRAMID
-	 */
-	virtual MeshElemType getGeomType() const { return MeshElemType::PYRAMID; }
-
-	/**
-	 * Get the type of the element in context of the finite element method.
-	 * @return a value of the enum CellType
-	 */
-	virtual CellType getCellType() const { return CELLPYRAMIDTYPE; }
-
-	/// Returns true if these two indeces form an edge and false otherwise
-	bool isEdge(unsigned i, unsigned j) const;
-
-	/**
-	 * Checks if a point is inside the element.
-	 * @param pnt a 3D MathLib::Point3d object
-	 * @param eps tolerance for numerical algorithm used or computing the property
-	 * @return true if the point is not outside the element, false otherwise
-	 */
-	bool isPntInElement(MathLib::Point3d const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
-
-	/**
-	 * Tests if the element is geometrically valid.
-	 * @param check_zero_volume indicates if volume == 0 should be checked
-	 */
-	virtual ElementErrorCode validate() const;
-
-	/**
-	 * Method clone is inherited from class Element. It makes a deep copy of the
-	 * TemplatePyramid instance employing the copy constructor of class TemplatePyramid.
-	 * @return an exact copy of the object
-	 */
-	virtual Element* clone() const;
-
-	/// Returns the ID of a face given an array of nodes.
-	unsigned identifyFace(Node* nodes[3]) const;
-
-protected:
-	/// Calculates the volume of a prism by subdividing it into two tetrahedra.
-	double computeVolume();
-
-	/// Return a specific edge node.
-	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const
-	{
-		return _nodes[_edge_nodes[edge_id][node_id]];
-	}
-
-	static const unsigned _face_nodes[5][4];
-	static const unsigned _edge_nodes[8][2];
-	static const unsigned _n_face_nodes[5];
-
-}; /* class */
-
-} /* namespace */
-
-#include "TemplatePyramid-impl.h"
-
-#endif /* TEMPLATEPYRAMID_H_ */
-
diff --git a/MeshLib/Elements/TemplateQuad-impl.h b/MeshLib/Elements/TemplateQuad-impl.h
deleted file mode 100644
index 864ab3f5385..00000000000
--- a/MeshLib/Elements/TemplateQuad-impl.h
+++ /dev/null
@@ -1,150 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date   2012-05-02
- * \brief  Implementation of the Quad class-
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include <array>
-
-#include "MeshLib/Node.h"
-
-#include "MathLib/MathTools.h"
-#include "GeoLib/AnalyticalGeometry.h"
-
-namespace MeshLib
-{
-
-template <unsigned NNODES, CellType CELLQUADTYPE>
-const unsigned TemplateQuad<NNODES, CELLQUADTYPE>::n_all_nodes;
-
-template <unsigned NNODES, CellType CELLQUADTYPE>
-const unsigned TemplateQuad<NNODES, CELLQUADTYPE>::n_base_nodes;
-
-template <unsigned NNODES, CellType CELLQUADTYPE>
-const unsigned TemplateQuad<NNODES, CELLQUADTYPE>::_edge_nodes[4][2] =
-{
-	{0, 1}, // Edge 0
-	{1, 2}, // Edge 1
-	{2, 3}, // Edge 2
-	{0, 3}  // Edge 3
-};
-
-template <unsigned NNODES, CellType CELLQUADTYPE>
-TemplateQuad<NNODES,CELLQUADTYPE>::TemplateQuad(Node* nodes[NNODES], unsigned value, std::size_t id)
-	: Face(value, id)
-{
-	_nodes = nodes;
-
-	_neighbors = new Element*[4];
-	std::fill(_neighbors, _neighbors + 4, nullptr);
-
-	this->_area = this->computeVolume();
-}
-
-template<unsigned NNODES, CellType CELLQUADTYPE>
-TemplateQuad<NNODES,CELLQUADTYPE>::TemplateQuad(std::array<Node*, NNODES> const& nodes,
-                                                unsigned value, std::size_t id)
-	: Face(value, id)
-{
-	_nodes = new Node*[NNODES];
-	std::copy(nodes.begin(), nodes.end(), _nodes);
-
-	_neighbors = new Element*[4];
-	std::fill(_neighbors, _neighbors + 4, nullptr);
-
-	this->_area = this->computeVolume();
-}
-
-template <unsigned NNODES, CellType CELLQUADTYPE>
-TemplateQuad<NNODES,CELLQUADTYPE>::TemplateQuad(const TemplateQuad<NNODES,CELLQUADTYPE> &quad)
-	: Face(quad.getValue(), quad.getID())
-{
-	_nodes = new Node*[NNODES];
-	for (unsigned i=0; i<NNODES; i++) {
-		_nodes[i] = quad._nodes[i];
-	}
-
-	_neighbors = new Element*[4];
-	for (unsigned i=0; i<4; i++) {
-		_neighbors[i] = quad._neighbors[i];
-	}
-
-	_area = quad.getArea();
-}
-
-template <unsigned NNODES, CellType CELLQUADTYPE>
-TemplateQuad<NNODES,CELLQUADTYPE>::~TemplateQuad()
-{
-}
-
-template <unsigned NNODES, CellType CELLQUADTYPE>
-double TemplateQuad<NNODES,CELLQUADTYPE>::computeVolume()
-{
-	return GeoLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2])
-         + GeoLib::calcTriangleArea(*_nodes[2], *_nodes[3], *_nodes[0]);
-}
-
-template <unsigned NNODES, CellType CELLQUADTYPE>
-bool TemplateQuad<NNODES,CELLQUADTYPE>::isEdge(unsigned idx1, unsigned idx2) const
-{
-	for (unsigned i(0); i<4; i++)
-	{
-		if (_edge_nodes[i][0]==idx1 && _edge_nodes[i][1]==idx2) return true;
-		if (_edge_nodes[i][1]==idx1 && _edge_nodes[i][0]==idx2) return true;
-	}
-	return false;
-}
-
-template <unsigned NNODES, CellType CELLQUADTYPE>
-bool TemplateQuad<NNODES,CELLQUADTYPE>::isPntInElement(MathLib::Point3d const& pnt, double eps) const
-{
-	return (GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[1], *_nodes[2], eps) ||
-	        GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[2], *_nodes[3], eps));
-}
-
-template <unsigned NNODES, CellType CELLQUADTYPE>
-Element* TemplateQuad<NNODES,CELLQUADTYPE>::clone() const
-{
-	return new TemplateQuad(*this);
-}
-
-template <unsigned NNODES, CellType CELLQUADTYPE>
-unsigned TemplateQuad<NNODES,CELLQUADTYPE>::identifyFace(Node* nodes[3]) const
-{
-	for (unsigned i=0; i<4; i++)
-	{
-		unsigned flag(0);
-		for (unsigned j=0; j<2; j++)
-			for (unsigned k=0; k<2; k++)
-				if (_nodes[_edge_nodes[i][j]] == nodes[k])
-					flag++;
-		if (flag==2)
-			return i;
-	}
-	return std::numeric_limits<unsigned>::max();
-}
-
-template <unsigned NNODES, CellType CELLQUADTYPE>
-ElementErrorCode TemplateQuad<NNODES,CELLQUADTYPE>::validate() const
-{
-	ElementErrorCode error_code;
-	error_code[ElementErrorFlag::ZeroVolume]  = this->hasZeroVolume();
-	error_code[ElementErrorFlag::NonCoplanar] = (!GeoLib::isCoplanar(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3]));
-	// for collapsed quads (i.e. reduced to a line) this test might result "false" as all four points are actually located on a line.
-	if (!error_code[ElementErrorFlag::ZeroVolume])
-		error_code[ElementErrorFlag::NonConvex]   = (!(GeoLib::dividedByPlane(*_nodes[0], *_nodes[2], *_nodes[1], *_nodes[3]) &&
-			                                           GeoLib::dividedByPlane(*_nodes[1], *_nodes[3], *_nodes[0], *_nodes[2])));
-	error_code[ElementErrorFlag::NodeOrder]  = !this->testElementNodeOrder();
-	return error_code;
-}
-
-}
-
diff --git a/MeshLib/Elements/TemplateQuad.h b/MeshLib/Elements/TemplateQuad.h
deleted file mode 100644
index 573e5474ece..00000000000
--- a/MeshLib/Elements/TemplateQuad.h
+++ /dev/null
@@ -1,140 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date   2012-05-02
- * \brief  Definition of the TemplateQuad class.
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef TEMPLATEQUAD_H_
-#define TEMPLATEQUAD_H_
-
-#include <array>
-#include "MeshLib/MeshEnums.h"
-#include "Face.h"
-
-namespace MeshLib {
-
-/**
- * This class represents a 2d quadliteral element. The following sketch shows the node and edge numbering.
- * @anchor QuadNodeAndEdgeNumbering
- * @code
- *              2
- *        3-----------2
- *        |           |
- *        |           |
- *       3|           |1
- *        |           |
- *        |           |
- *        0-----------1
- *              0
- * @endcode
- */
-template <unsigned NNODES, CellType CELLQUADTYPE>
-class TemplateQuad : public Face
-{
-public:
-	/// Constant: The number of all nodes for this element
-	static const unsigned n_all_nodes = NNODES;
-
-	/// Constant: The number of base nodes for this element
-	static const unsigned n_base_nodes = 4u;
-
-	/// Constructor with an array of mesh nodes.
-	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, std::size_t id = std::numeric_limits<std::size_t>::max());
-
-	/// Constructs a quad from NNODES of Nodes initializing Face with
-	//  value = 0.
-	TemplateQuad(Node* n0, Node* n1, Node* n2, Node* n3, ...);
-
-	/// Copy constructor
-	TemplateQuad(const TemplateQuad &quad);
-
-	/// Destructor
-	virtual ~TemplateQuad();
-
-	/// Get the number of edges for this element.
-	unsigned getNEdges() const { return 4; }
-
-	/// Get the number of neighbors for this element.
-	unsigned getNNeighbors() const { return 4; }
-
-	/// Get the number of linear nodes for this element.
-	virtual unsigned getNBaseNodes() const
-	{
-		return n_base_nodes;
-	}
-
-	/// Get the number of all nodes for this element.
-	virtual unsigned getNNodes() const
-	{
-		return n_all_nodes;
-	}
-
-	/**
-	 * Method returns the type of the element. In this case QUAD will be returned.
-	 * @return MeshElemType::QUAD
-	 */
-	virtual MeshElemType getGeomType() const { return MeshElemType::QUAD; }
-
-	/**
-	 * Get the type of the element in context of the finite element method.
-	 * @return a value of the enum CellType
-	 */
-	virtual CellType getCellType() const { return CELLQUADTYPE; }
-
-	/// Returns true if these two indeces form an edge and false otherwise
-	bool isEdge(unsigned i, unsigned j) const;
-
-	/**
-	 * Checks if a point is inside the element.
-	 * @param pnt a 3D MathLib::Point3d object
-	 * @param eps tolerance for numerical algorithm used or computing the property
-	 * @return true if the point is not outside the element, false otherwise
-	 */
-	bool isPntInElement(MathLib::Point3d const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
-
-	/**
-	 * Tests if the element is geometrically valid.
-	 * @param check_zero_volume indicates if area == 0 should be checked
-	 */
-	virtual ElementErrorCode validate() const;
-
-	/**
-	 * Method clone is inherited from class Element. It makes a deep copy of the TemplateQuad instance.
-	 * @return an exact copy of the object
-	 */
-	virtual Element* clone() const;
-
-	/// Returns the ID of a face given an array of nodes.
-	unsigned identifyFace(Node* nodes[3]) const;
-
-protected:
-	/// Calculates the area of a convex quadliteral by dividing it into two triangles.
-	double computeVolume();
-
-protected:
-	/// Return a specific edge node.
-	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const
-	{
-		return _nodes[_edge_nodes[edge_id][node_id]];
-	}
-
-	static const unsigned _edge_nodes[4][2];
-}; /* class */
-
-} /* namespace */
-
-#include "TemplateQuad-impl.h"
-
-#endif /* TEMPLATEQUAD_H_ */
-
diff --git a/MeshLib/Elements/TemplateTet-impl.h b/MeshLib/Elements/TemplateTet-impl.h
deleted file mode 100644
index 5ddb6c5176e..00000000000
--- a/MeshLib/Elements/TemplateTet-impl.h
+++ /dev/null
@@ -1,168 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date   2012-05-02
- * \brief  Implementation of the TemplateTet class.
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include "logog/include/logog.hpp"
-
-#include "GeoLib/AnalyticalGeometry.h"
-#include "MeshLib/Node.h"
-#include "Tri.h"
-
-namespace MeshLib {
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-const unsigned TemplateTet<NNODES, CELLTETTYPE>::n_all_nodes;
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-const unsigned TemplateTet<NNODES, CELLTETTYPE>::n_base_nodes;
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-const unsigned TemplateTet<NNODES,CELLTETTYPE>::_face_nodes[4][3] =
-{
-	{0, 2, 1}, // Face 0
-	{0, 1, 3}, // Face 1
-	{1, 2, 3}, // Face 2
-	{2, 0, 3}  // Face 3
-};
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-const unsigned TemplateTet<NNODES,CELLTETTYPE>::_edge_nodes[6][2] =
-{
-	{0, 1}, // Edge 0
-	{1, 2}, // Edge 1
-	{0, 2}, // Edge 2
-	{0, 3}, // Edge 3
-	{1, 3}, // Edge 4
-	{2, 3}  // Edge 5
-};
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-TemplateTet<NNODES,CELLTETTYPE>::TemplateTet(Node* nodes[NNODES], unsigned value, std::size_t id)
-	: Cell(value, id)
-{
-	_nodes = nodes;
-
-	_neighbors = new Element*[4];
-	std::fill(_neighbors, _neighbors + 4, nullptr);
-
-	this->_volume = this->computeVolume();
-}
-
-template<unsigned NNODES, CellType CELLTETTYPE>
-TemplateTet<NNODES,CELLTETTYPE>::TemplateTet(std::array<Node*, NNODES> const& nodes,
-                                             unsigned value, std::size_t id)
-	: Cell(value, id)
-{
-	_nodes = new Node*[NNODES];
-	std::copy(nodes.begin(), nodes.end(), _nodes);
-
-	_neighbors = new Element*[4];
-	std::fill(_neighbors, _neighbors + 4, nullptr);
-
-	this->_volume = this->computeVolume();
-}
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-TemplateTet<NNODES,CELLTETTYPE>::TemplateTet(const TemplateTet<NNODES,CELLTETTYPE> &tet)
-	: Cell(tet.getValue(), tet.getID())
-{
-	_nodes = new Node*[NNODES];
-	for (unsigned i=0; i<NNODES; i++) {
-		_nodes[i] = tet._nodes[i];
-	}
-
-	_neighbors = new Element*[4];
-	for (unsigned i=0; i<4; i++)
-	{
-		_neighbors[i] = tet._neighbors[i];
-	}
-
-	_volume = tet.getVolume();
-}
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-TemplateTet<NNODES,CELLTETTYPE>::~TemplateTet()
-{
-}
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-double TemplateTet<NNODES,CELLTETTYPE>::computeVolume()
-{
-	return GeoLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords());
-}
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-const Element* TemplateTet<NNODES,CELLTETTYPE>::getFace(unsigned i) const
-{
-	if (i<this->getNFaces())
-	{
-		unsigned nFaceNodes (this->getNFaceNodes(i));
-		Node** nodes = new Node*[nFaceNodes];
-		for (unsigned j=0; j<nFaceNodes; j++)
-			nodes[j] = _nodes[_face_nodes[i][j]];
-		return new Tri(nodes);
-	}
-	ERR("Error in MeshLib::Element::getFace() - Index %d does not exist.", i);
-	return NULL;
-}
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-bool TemplateTet<NNODES,CELLTETTYPE>::isEdge(unsigned idx1, unsigned idx2) const
-{
-	for (unsigned i(0); i<6; i++)
-	{
-		if (_edge_nodes[i][0]==idx1 && _edge_nodes[i][1]==idx2) return true;
-		if (_edge_nodes[i][1]==idx1 && _edge_nodes[i][0]==idx2) return true;
-	}
-	return false;
-}
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-bool TemplateTet<NNODES,CELLTETTYPE>::isPntInElement(MathLib::Point3d const& pnt, double eps) const
-{
-    return GeoLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3], eps);
-}
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-Element* TemplateTet<NNODES,CELLTETTYPE>::clone() const
-{
-	return new TemplateTet<NNODES,CELLTETTYPE>(*this);
-}
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-unsigned TemplateTet<NNODES,CELLTETTYPE>::identifyFace(Node* nodes[3]) const
-{
-	for (unsigned i=0; i<4; i++)
-	{
-		unsigned flag(0);
-		for (unsigned j=0; j<3; j++)
-			for (unsigned k=0; k<3; k++)
-				if (_nodes[_face_nodes[i][j]] == nodes[k])
-					flag++;
-		if (flag==3)
-			return i;
-	}
-	return std::numeric_limits<unsigned>::max();
-}
-
-template <unsigned NNODES, CellType CELLTETTYPE>
-ElementErrorCode TemplateTet<NNODES,CELLTETTYPE>::validate() const
-{
-	ElementErrorCode error_code;
-	error_code[ElementErrorFlag::ZeroVolume] = this->hasZeroVolume();
-	error_code[ElementErrorFlag::NodeOrder]  = !this->testElementNodeOrder();
-	return error_code;
-}
-
-} // end namespace MeshLib
-
diff --git a/MeshLib/Elements/TemplateTet.h b/MeshLib/Elements/TemplateTet.h
deleted file mode 100644
index dce10086ad0..00000000000
--- a/MeshLib/Elements/TemplateTet.h
+++ /dev/null
@@ -1,155 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date   2012-05-02
- * \brief  Definition of the TemplateTet class.
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef TEMPLATETET_H_
-#define TEMPLATETET_H_
-
-#include <array>
-#include "MeshLib/MeshEnums.h"
-#include "Cell.h"
-
-namespace MeshLib {
-
-/**
- * This class represents a 3d tetrahedron element. The following sketch shows the node and edge numbering.
- * @anchor TetrahedronNodeAndEdgeNumbering
- * @code
- *          3
- *         /|\
- *        / | \
- *      3/  |  \5
- *      /   |4  \
- *     /    |    \
- *    0.....|.....2
- *     \    |  2 /
- *      \   |   /
- *      0\  |  /1
- *        \ | /
- *         \|/
- *          1
- *
- * @endcode
- */
-template <unsigned NNODES, CellType CELLTETTYPE>
-class TemplateTet : public Cell
-{
-public:
-	/// Constant: The number of all nodes for this element
-	static const unsigned n_all_nodes = NNODES;
-
-	/// Constant: The number of base nodes for this element
-	static const unsigned n_base_nodes = 4u;
-
-	/// Constructor with an array of mesh nodes.
-	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, std::size_t id = std::numeric_limits<std::size_t>::max());
-
-	/// Copy constructor
-	TemplateTet(const TemplateTet &tet);
-
-	/// Destructor
-	virtual ~TemplateTet();
-
-	/// Returns the face i of the element.
-	const Element* getFace(unsigned i) const;
-
-	/// Get the number of edges for this element.
-	unsigned getNEdges() const { return 6; }
-
-	/// Get the number of nodes for face i.
-	unsigned getNFaceNodes(unsigned i) const { (void)i; return 3; }
-
-	/// Get the number of faces for this element.
-	unsigned getNFaces() const { return 4; }
-
-	/// Get the number of neighbors for this element.
-	unsigned getNNeighbors() const { return 4; }
-
-	/// Get the number of linear nodes for this element.
-	virtual unsigned getNBaseNodes() const
-	{
-		return n_base_nodes;
-	}
-
-	/// Get the number of all nodes for this element.
-	virtual unsigned getNNodes() const
-	{
-		return n_all_nodes;
-	}
-
-	/**
-	 * Method returns the type of the element. In this case TETRAHEDRON will be returned.
-	 * @return MeshElemType::TETRAHEDRON
-	 */
-	virtual MeshElemType getGeomType() const { return MeshElemType::TETRAHEDRON; }
-
-	/**
-	 * Get the type of the element in context of the finite element method.
-	 * @return a value of the enum CellType
-	 */
-	virtual CellType getCellType() const { return CELLTETTYPE; }
-
-	/// Returns true if these two indeces form an edge and false otherwise
-	bool isEdge(unsigned i, unsigned j) const;
-
-	/**
-	 * Checks if a point is inside the element.
-	 * @param pnt a 3D MathLib::Point3d object
-	 * @param eps tolerance for numerical algorithm used or computing the property
-	 * @return true if the point is not outside the element, false otherwise
-	 */
-	bool isPntInElement(MathLib::Point3d const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
-
-	/**
-	 * Tests if the element is geometrically valid.
-	 * @param check_zero_volume indicates if volume == 0 should be checked
-	 */
-	virtual ElementErrorCode validate() const;
-
-	/**
-	 * Method clone is inherited from class Element. It makes a deep copy of the TemplateTet instance.
-	 * @return an exact copy of the object
-	 */
-	virtual Element* clone() const;
-
-	/// Returns the ID of a face given an array of nodes.
-	unsigned identifyFace(Node* nodes[3]) const;
-
-protected:
-	/// Calculates the volume of a tetrahedron via the determinant of the matrix given by its four points.
-	double computeVolume();
-
-	/**
-	 * Return a specific edge node, see @ref TetrahedronNodeAndEdgeNumbering for numbering
-	 * @param edge_id the id/number of the edge, have to be an integer value in the interval [0,6]
-	 * @param node_id the id of the node within the edge (either 0 or 1)
-	 * @return a pointer to the internal Node
-	 */
-	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const {
-		return _nodes[_edge_nodes[edge_id][node_id]];
-	}
-
-	static const unsigned _face_nodes[4][3];
-	static const unsigned _edge_nodes[6][2];
-
-}; /* class */
-
-} /* namespace */
-
-#include "TemplateTet-impl.h"
-
-#endif /* TEMPLATETET_H_ */
-
diff --git a/MeshLib/Elements/TemplateTri-impl.h b/MeshLib/Elements/TemplateTri-impl.h
deleted file mode 100644
index 1e204207cd1..00000000000
--- a/MeshLib/Elements/TemplateTri-impl.h
+++ /dev/null
@@ -1,119 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date   Sep 27, 2012
- * \brief  Implementation of the TemplateTri class.
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-namespace MeshLib {
-
-template <unsigned NNODES, CellType CELLTRITYPE>
-const unsigned TemplateTri<NNODES, CELLTRITYPE>::n_all_nodes;
-
-template <unsigned NNODES, CellType CELLTRITYPE>
-const unsigned TemplateTri<NNODES, CELLTRITYPE>::n_base_nodes;
-
-template <unsigned NNODES, CellType CELLTRITYPE>
-const unsigned TemplateTri<NNODES,CELLTRITYPE>::_edge_nodes[3][2] = {
-		{0, 1}, // Edge 0
-		{1, 2}, // Edge 1
-		{0, 2}  // Edge 2
-	};
-
-template <unsigned NNODES, CellType CELLTRITYPE>
-TemplateTri<NNODES,CELLTRITYPE>::TemplateTri(Node* nodes[NNODES], unsigned value, std::size_t id) :
-	Face(value, id)
-{
-	_nodes = nodes;
-	_neighbors = new Element*[3];
-	std::fill(_neighbors, _neighbors + 3, nullptr);
-	this->_area = this->computeVolume();
-}
-
-template<unsigned NNODES, CellType CELLTRITYPE>
-TemplateTri<NNODES,CELLTRITYPE>::TemplateTri(std::array<Node*, NNODES> const& nodes,
-                                             unsigned value, std::size_t id)
-	: Face(value, id)
-{
-	_nodes = new Node*[NNODES];
-	std::copy(nodes.begin(), nodes.end(), _nodes);
-
-	_neighbors = new Element*[3];
-	std::fill(_neighbors, _neighbors + 3, nullptr);
-
-	this->_area = this->computeVolume();
-}
-
-template <unsigned NNODES, CellType CELLTRITYPE>
-TemplateTri<NNODES,CELLTRITYPE>::TemplateTri(const TemplateTri<NNODES,CELLTRITYPE> &tri) :
-	Face(tri.getValue(), tri.getID())
-{
-	_nodes = new Node*[NNODES];
-	for (unsigned i=0; i<NNODES; i++)
-	{
-		_nodes[i] = tri._nodes[i];
-	}
-
-	_neighbors = new Element*[3];
-	for (unsigned i=0; i<3; i++) {
-		_neighbors[i] = tri._neighbors[i];
-	}
-
-	_area = tri.getArea();
-}
-
-template <unsigned NNODES, CellType CELLTRITYPE>
-TemplateTri<NNODES,CELLTRITYPE>::~TemplateTri()
-{}
-
-template <unsigned NNODES, CellType CELLTRITYPE>
-bool TemplateTri<NNODES,CELLTRITYPE>::isEdge(unsigned idx1, unsigned idx2) const
-{
-	for (unsigned i(0); i<3; i++)
-	{
-		if (_edge_nodes[i][0]==idx1 && _edge_nodes[i][1]==idx2) return true;
-		if (_edge_nodes[i][1]==idx1 && _edge_nodes[i][0]==idx2) return true;
-	}
-	return false;
-}
-
-template <unsigned NNODES, CellType CELLTRITYPE>
-bool TemplateTri<NNODES,CELLTRITYPE>::isPntInElement(MathLib::Point3d const& pnt, double eps) const
-{
-	return GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[1], *_nodes[2], eps);
-}
-
-template <unsigned NNODES, CellType CELLTRITYPE>
-ElementErrorCode TemplateTri<NNODES,CELLTRITYPE>::validate() const 
-{ 
-	ElementErrorCode error_code;
-	error_code[ElementErrorFlag::ZeroVolume] = this->hasZeroVolume();
-	error_code[ElementErrorFlag::NodeOrder]  = !this->testElementNodeOrder();
-	return error_code;
-}
-
-template <unsigned NNODES, CellType CELLTRITYPE>
-unsigned TemplateTri<NNODES,CELLTRITYPE>::identifyFace(Node* nodes[3]) const
-{
-	for (unsigned i=0; i<3; i++)
-	{
-		unsigned flag(0);
-		for (unsigned j=0; j<2; j++)
-			for (unsigned k=0; k<2; k++)
-				if (_nodes[_edge_nodes[i][j]] == nodes[k])
-					flag++;
-		if (flag==2)
-			return i;
-	}
-	return std::numeric_limits<unsigned>::max();
-}
-
-} // namespace MeshLib
-
diff --git a/MeshLib/Elements/TemplateTri.h b/MeshLib/Elements/TemplateTri.h
deleted file mode 100644
index f5e20108809..00000000000
--- a/MeshLib/Elements/TemplateTri.h
+++ /dev/null
@@ -1,154 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date   2012-05-02
- * \brief  Definition of the TemplateTri class.
- *
- * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef TEMPLATETRI_H_
-#define TEMPLATETRI_H_
-
-#include <array>
-
-#include "GeoLib/AnalyticalGeometry.h"
-
-#include "MathLib/MathTools.h"
-
-#include "MeshLib/MeshEnums.h"
-#include "MeshLib/Node.h"
-
-#include "Face.h"
-#include "Line.h"
-
-namespace MeshLib {
-
-/**
- * This class represents a 2d triangle element. The following sketch shows the node and edge numbering.
- * @anchor TriNodeAndEdgeNumbering
- * @code
- *
- *          2
- *          o
- *         / \
- *        /   \
- *      2/     \1
- *      /       \
- *     /         \
- *    0-----------1
- *          0
- *
- * @endcode
- */
-template <unsigned NNODES, CellType CELLTRITYPE>
-class TemplateTri : public Face
-{
-public:
-	/// Constant: The number of all nodes for this element
-	static const unsigned n_all_nodes = NNODES;
-
-	/// Constant: The number of base nodes for this element
-	static const unsigned n_base_nodes = 3u;
-
-	/// Constructor with an array of mesh nodes.
-	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, std::size_t id = std::numeric_limits<std::size_t>::max());
-
-	/// Copy constructor
-	TemplateTri(const TemplateTri &tri);
-
-	/// Destructor
-	virtual ~TemplateTri();
-
-	/// Get the number of edges for this element.
-	unsigned getNEdges() const { return 3; }
-
-	/// Get the number of neighbors for this element.
-	unsigned getNNeighbors() const { return 3; }
-
-	/// Get the number of linear nodes for this element.
-	virtual unsigned getNBaseNodes() const
-	{
-		return n_base_nodes;
-	}
-
-	/// Get the number of all nodes for this element.
-	virtual unsigned getNNodes() const
-	{
-		return n_all_nodes;
-	}
-
-	/**
-	 * Method returns the type of the element. In this case TRIANGLE will be returned.
-	 * @return MeshElemType::TRIANGLE
-	 */
-	virtual MeshElemType getGeomType() const { return MeshElemType::TRIANGLE; }
-
-	/**
-	 * Get the type of the element in context of the finite element method.
-	 * @return a value of the enum CellType
-	 */
-	virtual CellType getCellType() const { return CELLTRITYPE; }
-
-	/// Returns true if these two indices form an edge and false otherwise
-	bool isEdge(unsigned idx1, unsigned idx2) const;
-
-	/**
-	 * Checks if a point is inside the element.
-	 * @param pnt a 3D MathLib::Point3d object
-	 * @param eps tolerance for numerical algorithm used or computing the property
-	 * @return true if the point is not outside the element, false otherwise
-	 */
-	bool isPntInElement(MathLib::Point3d const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
-
-	/**
-	 * Tests if the element is geometrically valid
-	 * @param check_zero_volume indicates if area == 0 should be checked
-	 * @return error code (0 = okay, 1 = zero volume)
-	 */
-	virtual ElementErrorCode validate() const;
-
-
-	/**
-	 * Method clone is inherited from class Element. It makes a deep copy of the TemplateTri instance.
-	 * @return an exact copy of the object
-	 */
-	virtual Element* clone() const
-	{
-		return new TemplateTri<NNODES,CELLTRITYPE>(*this);
-	}
-
-	/// Returns the ID of a face given an array of nodes.
-	unsigned identifyFace(Node* nodes[3]) const;
-
-protected:
-	/// Calculates the area of the triangle by returning half of the area of the corresponding parallelogram.
-	double computeVolume()
-	{
-		return GeoLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2]);
-	}
-
-protected:
-	/// Return a specific edge node.
-	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const
-	{
-		return _nodes[_edge_nodes[edge_id][node_id]];
-	}
-
-	static const unsigned _edge_nodes[3][2];
-}; /* class */
-
-} /* namespace */
-
-#include "TemplateTri-impl.h"
-
-#endif /* TEMPLATETRI_H_ */
-
diff --git a/MeshLib/Elements/Tet.h b/MeshLib/Elements/Tet.h
index eaf0ac2a05e..3505c2d849e 100644
--- a/MeshLib/Elements/Tet.h
+++ b/MeshLib/Elements/Tet.h
@@ -15,11 +15,13 @@
 #ifndef TET_H_
 #define TET_H_
 
-#include "TemplateTet.h"
+#include "TemplateElement.h"
+#include "Cell.h"
+#include "TetRule4.h"
 
 namespace MeshLib {
 
-typedef TemplateTet<4,CellType::TET4> Tet;
+typedef TemplateElement<Cell,TetRule4> Tet;
 
 }
 
diff --git a/MeshLib/Elements/TetRule4.cpp b/MeshLib/Elements/TetRule4.cpp
new file mode 100644
index 00000000000..5958b475d74
--- /dev/null
+++ b/MeshLib/Elements/TetRule4.cpp
@@ -0,0 +1,91 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "TetRule4.h"
+
+#include "logog/include/logog.hpp"
+
+#include "GeoLib/AnalyticalGeometry.h"
+
+#include "MeshLib/Node.h"
+#include "Tri.h"
+#include "Line.h"
+
+namespace MeshLib {
+
+const unsigned TetRule4::n_all_nodes;
+
+const unsigned TetRule4::n_base_nodes;
+
+const unsigned TetRule4::face_nodes[4][3] =
+{
+	{0, 2, 1}, // Face 0
+	{0, 1, 3}, // Face 1
+	{1, 2, 3}, // Face 2
+	{2, 0, 3}  // Face 3
+};
+
+const unsigned TetRule4::edge_nodes[6][2] =
+{
+	{0, 1}, // Edge 0
+	{1, 2}, // Edge 1
+	{0, 2}, // Edge 2
+	{0, 3}, // Edge 3
+	{1, 3}, // Edge 4
+	{2, 3}  // Edge 5
+};
+
+const Element* TetRule4::getFace(const Element* e, unsigned i)
+{
+	if (i<e->getNFaces())
+	{
+		unsigned nFaceNodes (e->getNFaceNodes(i));
+		Node** nodes = new Node*[nFaceNodes];
+		for (unsigned j=0; j<nFaceNodes; j++)
+			nodes[j] = const_cast<Node*>(e->getNode(face_nodes[i][j]));
+		return new Tri(nodes);
+	}
+	ERR("Error in MeshLib::Element::getFace() - Index %d does not exist.", i);
+	return nullptr;
+}
+
+double TetRule4::computeVolume(Node const* const* _nodes)
+{
+	return GeoLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords());
+}
+
+bool TetRule4::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps)
+{
+	return GeoLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3], eps);
+}
+
+unsigned TetRule4::identifyFace(Node const* const* _nodes, Node* nodes[3])
+{
+	for (unsigned i=0; i<4; i++)
+	{
+		unsigned flag(0);
+		for (unsigned j=0; j<3; j++)
+			for (unsigned k=0; k<3; k++)
+				if (_nodes[face_nodes[i][j]] == nodes[k])
+					flag++;
+		if (flag==3)
+			return i;
+	}
+	return std::numeric_limits<unsigned>::max();
+}
+
+ElementErrorCode TetRule4::validate(const Element* e)
+{
+	ElementErrorCode error_code;
+	error_code[ElementErrorFlag::ZeroVolume] = e->hasZeroVolume();
+	error_code[ElementErrorFlag::NodeOrder]  = !e->testElementNodeOrder();
+	return error_code;
+}
+
+} // end namespace MeshLib
diff --git a/MeshLib/Elements/TetRule4.h b/MeshLib/Elements/TetRule4.h
new file mode 100644
index 00000000000..4e94b3a09da
--- /dev/null
+++ b/MeshLib/Elements/TetRule4.h
@@ -0,0 +1,105 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef TETRULE4_H_
+#define TETRULE4_H_
+
+#include "MeshLib/MeshEnums.h"
+#include "Element.h"
+#include "CellRule.h"
+#include "EdgeReturn.h"
+
+namespace MeshLib
+{
+
+/**
+ * This class represents a 3d tetrahedron element. The following sketch shows the node and edge numbering.
+ * @anchor TetrahedronNodeAndEdgeNumbering
+ * @code
+ *          3
+ *         /|\
+ *        / | \
+ *      3/  |  \5
+ *      /   |4  \
+ *     /    |    \
+ *    0.....|.....2
+ *     \    |  2 /
+ *      \   |   /
+ *      0\  |  /1
+ *        \ | /
+ *         \|/
+ *          1
+ *
+ * @endcode
+ */
+class TetRule4 : public CellRule
+{
+public:
+	/// Constant: The number of base nodes for this element
+	static const unsigned n_base_nodes = 4u;
+
+	/// Constant: The number of all nodes for this element
+	static const unsigned n_all_nodes = 4u;
+
+	/// Constant: The geometric type of the element
+	static const MeshElemType mesh_elem_type = MeshElemType::TETRAHEDRON;
+
+	/// Constant: The FEM type of the element
+	static const CellType cell_type = CellType::TET4;
+
+	/// Constant: The number of faces
+	static const unsigned n_faces = 4;
+
+	/// Constant: The number of edges
+	static const unsigned n_edges = 6;
+
+	/// Constant: The number of neighbors
+	static const unsigned n_neighbors = 4;
+
+	/// Constant: Local node index table for faces
+	static const unsigned face_nodes[4][3];
+
+	/// Constant: Local node index table for edge
+	static const unsigned edge_nodes[6][2];
+
+	/// Returns the i-th edge of the element.
+	typedef LinearEdgeReturn EdgeReturn;
+
+	/// Get the number of nodes for face i.
+	static unsigned getNFaceNodes(unsigned /*i*/) { return 3; }
+
+	/// Returns the i-th face of the element.
+	static const Element* getFace(const Element* e, unsigned i);
+
+	/**
+	 * Checks if a point is inside the element.
+	 * @param pnt a 3D GeoLib::Point object
+	 * @param eps tolerance for numerical algorithm used or computing the property
+	 * @return true if the point is not outside the element, false otherwise
+	 */
+	static bool isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps);
+
+	/**
+	 * Tests if the element is geometrically valid.
+	 * @param check_zero_volume indicates if volume == 0 should be checked
+	 */
+	static ElementErrorCode validate(const Element* e);
+
+	/// Returns the ID of a face given an array of nodes.
+	static unsigned identifyFace(Node const* const*, Node* nodes[3]);
+
+	/// Calculates the volume of the element
+	static double computeVolume(Node const* const* _nodes);
+
+}; /* class */
+
+} /* namespace */
+
+#endif /* TETRULE4_H_ */
+
diff --git a/MeshLib/Elements/Tri.h b/MeshLib/Elements/Tri.h
index 54266624253..78c6ed84d39 100644
--- a/MeshLib/Elements/Tri.h
+++ b/MeshLib/Elements/Tri.h
@@ -15,11 +15,13 @@
 #ifndef TRI_H_
 #define TRI_H_
 
-#include "TemplateTri.h"
+#include "TemplateElement.h"
+#include "Face.h"
+#include "TriRule3.h"
 
 namespace MeshLib {
 
-typedef TemplateTri<3,CellType::TRI3> Tri;
+typedef TemplateElement<Face,TriRule3> Tri;
 
 }
 
diff --git a/MeshLib/Elements/TriRule3.cpp b/MeshLib/Elements/TriRule3.cpp
new file mode 100644
index 00000000000..b6d02af485b
--- /dev/null
+++ b/MeshLib/Elements/TriRule3.cpp
@@ -0,0 +1,64 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "TriRule3.h"
+
+#include "logog/include/logog.hpp"
+
+#include "GeoLib/AnalyticalGeometry.h"
+
+#include "MeshLib/Node.h"
+
+namespace MeshLib {
+
+const unsigned TriRule3::n_all_nodes;
+
+const unsigned TriRule3::n_base_nodes;
+
+const unsigned TriRule3::edge_nodes[3][2] =
+{
+		{0, 1}, // Edge 0
+		{1, 2}, // Edge 1
+		{2, 0}, // Edge 2
+};
+
+double TriRule3::computeVolume(Node const* const* _nodes)
+{
+	return GeoLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2]);
+}
+
+bool TriRule3::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps)
+{
+	return GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[1], *_nodes[2], eps);
+}
+
+unsigned TriRule3::identifyFace(Node const* const* _nodes, Node* nodes[3])
+{
+	for (unsigned i=0; i<3; i++)
+	{
+		unsigned flag(0);
+		for (unsigned j=0; j<2; j++)
+			for (unsigned k=0; k<2; k++)
+				if (_nodes[edge_nodes[i][j]] == nodes[k])
+					flag++;
+		if (flag==2)
+			return i;
+	}
+	return std::numeric_limits<unsigned>::max();
+}
+
+ElementErrorCode TriRule3::validate(const Element* e)
+{
+	ElementErrorCode error_code;
+	error_code[ElementErrorFlag::ZeroVolume] = e->hasZeroVolume();
+	error_code[ElementErrorFlag::NodeOrder]  = !e->testElementNodeOrder();
+	return error_code;
+}
+
+} // end namespace MeshLib
diff --git a/MeshLib/Elements/TriRule3.h b/MeshLib/Elements/TriRule3.h
new file mode 100644
index 00000000000..1352c6d0aaf
--- /dev/null
+++ b/MeshLib/Elements/TriRule3.h
@@ -0,0 +1,93 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef TRIRULE3_H_
+#define TRIRULE3_H_
+
+#include "MeshLib/MeshEnums.h"
+#include "Element.h"
+#include "FaceRule.h"
+#include "EdgeReturn.h"
+
+namespace MeshLib
+{
+
+/**
+ * This class represents a 2d triangle element. The following sketch shows the node and edge numbering.
+ * @anchor TriNodeAndEdgeNumbering
+ * @code
+ *
+ *          2
+ *          o
+ *         / \
+ *        /   \
+ *      2/     \1
+ *      /       \
+ *     /         \
+ *    0-----------1
+ *          0
+ *
+ * @endcode
+ */
+class TriRule3 : public FaceRule
+{
+public:
+	/// Constant: The number of base nodes for this element
+	static const unsigned n_base_nodes = 3u;
+
+	/// Constant: The number of all nodes for this element
+	static const unsigned n_all_nodes = 3u;
+
+	/// Constant: The geometric type of the element
+	static const MeshElemType mesh_elem_type = MeshElemType::TRIANGLE;
+
+	/// Constant: The FEM type of the element
+	static const CellType cell_type = CellType::TRI3;
+
+	/// Constant: The number of faces
+	static const unsigned n_faces = 0;
+
+	/// Constant: The number of edges
+	static const unsigned n_edges = 3;
+
+	/// Constant: The number of neighbors
+	static const unsigned n_neighbors = 3;
+
+	/// Constant: Local node index table for edge
+	static const unsigned edge_nodes[3][2];
+
+	/// Returns the i-th edge of the element.
+	typedef LinearEdgeReturn EdgeReturn;
+
+	/**
+	 * Checks if a point is inside the element.
+	 * @param pnt a 3D GeoLib::Point object
+	 * @param eps tolerance for numerical algorithm used or computing the property
+	 * @return true if the point is not outside the element, false otherwise
+	 */
+	static bool isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps);
+
+	/**
+	 * Tests if the element is geometrically valid.
+	 * @param check_zero_volume indicates if volume == 0 should be checked
+	 */
+	static ElementErrorCode validate(const Element* e);
+
+	/// Returns the ID of a face given an array of nodes.
+	static unsigned identifyFace(Node const* const*, Node* nodes[3]);
+
+	/// Calculates the volume of a convex hexahedron by partitioning it into six tetrahedra.
+	static double computeVolume(Node const* const* _nodes);
+
+}; /* class */
+
+} /* namespace */
+
+#endif /* TRIRULE3_H_ */
+
-- 
GitLab