diff --git a/MeshLib/Elements/Edge.h b/MeshLib/Elements/Edge.h
index 4f23c6b762ccc574732e7cfb04046c1b4f56ec66..d0f9e5ae5d3caf827b5660480e66277f44a330a3 100644
--- a/MeshLib/Elements/Edge.h
+++ b/MeshLib/Elements/Edge.h
@@ -16,7 +16,7 @@
 
 namespace MeshLib {
 
-typedef TemplateEdge<2> Edge;
+typedef TemplateEdge<2,FEMElemType::EDGE2> Edge;
 
 }
 
diff --git a/MeshLib/Elements/Element.h b/MeshLib/Elements/Element.h
index 6fef4bf7fe436eb18e8a906c24796d924c48d991..56fc9d62868e301f4a7dacbb24a77feb97034de6 100644
--- a/MeshLib/Elements/Element.h
+++ b/MeshLib/Elements/Element.h
@@ -113,9 +113,17 @@ public:
 	 */
 	unsigned getNodeIndex(unsigned i) const;
 
-	/// Get the type of the mesh element (as a MshElemType-enum).
+	/**
+	 * Get the type of the mesh element in geometric context (as a MshElemType-enum).
+	 */
 	virtual MshElemType::type getGeoType() const = 0;
 
+	/**
+	 * Get the type of the element in context of the finite element method.
+	 * @return a value of the enum FEMElemType::type
+	 */
+	virtual FEMElemType::type getFEMType() const = 0;
+
 	/// Get the value for this element.
 	unsigned getValue() const { return _value; };
 
diff --git a/MeshLib/Elements/Hex.h b/MeshLib/Elements/Hex.h
index 98c082ab04e3c38d3855783bf77182438f8525cf..5e42c3230b8f41f37d1f9016438dca50714b566f 100644
--- a/MeshLib/Elements/Hex.h
+++ b/MeshLib/Elements/Hex.h
@@ -15,7 +15,7 @@
 #include "TemplateHex.h"
 
 namespace MeshLib {
-typedef TemplateHex<8> Hex;
+typedef TemplateHex<8, FEMElemType::HEX8> Hex;
 }
 
 #endif /* HEX_H_ */
diff --git a/MeshLib/Elements/Prism.h b/MeshLib/Elements/Prism.h
index b17fd500b8071071a0f92be7efab58aef98dc522..23bebf882e606f9dc4ba2c075265d674890f5984 100644
--- a/MeshLib/Elements/Prism.h
+++ b/MeshLib/Elements/Prism.h
@@ -16,7 +16,7 @@
 
 namespace MeshLib {
 
-typedef TemplatePrism<6> Prism;
+typedef TemplatePrism<6, FEMElemType::PRISM6> Prism;
 
 }
 
diff --git a/MeshLib/Elements/Pyramid.h b/MeshLib/Elements/Pyramid.h
index 96343c38730721b9e0f4fc818727bf13428f585e..eb194e8e989c9defbdff9dbc778ed8f9c170a043 100644
--- a/MeshLib/Elements/Pyramid.h
+++ b/MeshLib/Elements/Pyramid.h
@@ -16,7 +16,7 @@
 
 namespace MeshLib {
 
-typedef TemplatePyramid<5> Pyramid;
+typedef TemplatePyramid<5,FEMElemType::PYRAMID5> Pyramid;
 
 }
 
diff --git a/MeshLib/Elements/Quad.h b/MeshLib/Elements/Quad.h
index 8b1e3af10e283c12414983799d264fac94f4e315..d8339208c75d73bd90ef68b37c05ceffe5146d69 100644
--- a/MeshLib/Elements/Quad.h
+++ b/MeshLib/Elements/Quad.h
@@ -16,7 +16,7 @@
 
 namespace MeshLib {
 
-typedef TemplateQuad<4> Quad;
+typedef TemplateQuad<4,FEMElemType::QUAD4> Quad;
 
 }
 
diff --git a/MeshLib/Elements/TemplateEdge.h b/MeshLib/Elements/TemplateEdge.h
index 47386c34916a688930e20a82d6c2f2589ebf7c96..e04e88932971f81c94bbb3808757da5ef12c3ea2 100644
--- a/MeshLib/Elements/TemplateEdge.h
+++ b/MeshLib/Elements/TemplateEdge.h
@@ -15,6 +15,7 @@
 
 #include <limits>
 
+#include "MshEnums.h"
 #include "Element.h"
 #include "Node.h"
 
@@ -29,7 +30,7 @@ namespace MeshLib {
  *  0--------1
  * @endcode
  */
-template <unsigned NNODES>
+template<unsigned NNODES, FEMElemType::type FEMEDGETYPE>
 class TemplateEdge : public Element
 {
 public:
@@ -37,7 +38,7 @@ public:
 	TemplateEdge(Node* nodes[NNODES], unsigned value = 0);
 
 	/// Copy constructor
-	TemplateEdge(const TemplateEdge<NNODES> &edge);
+	TemplateEdge(const TemplateEdge &edge);
 
 	/// Destructor
 	virtual ~TemplateEdge();
@@ -84,6 +85,12 @@ public:
 	 */
 	virtual MshElemType::type getGeoType() const { return MshElemType::EDGE; }
 
+	/**
+	 * Get the type of the element in context of the finite element method.
+	 * @return a value of the enum FEMElemType::type
+	 */
+	virtual FEMElemType::type getFEMType() const { return FEMEDGETYPE; }
+
 	/// Returns true if these two indices form an edge and false otherwise
 	bool isEdge(unsigned idx1, unsigned idx2) const
 	{
@@ -94,7 +101,7 @@ public:
 
 	virtual Element* clone() const
 	{
-		return new TemplateEdge<NNODES>(*this);
+		return new TemplateEdge<NNODES,FEMEDGETYPE>(*this);
 	}
 
 	/**
diff --git a/MeshLib/Elements/TemplateEdge.tpp b/MeshLib/Elements/TemplateEdge.tpp
index 19a44c28e8d94a48a7baccf125512cb10a83b4a9..d33393bc2cf316c01fbcc71b473f2128b1e3e454 100644
--- a/MeshLib/Elements/TemplateEdge.tpp
+++ b/MeshLib/Elements/TemplateEdge.tpp
@@ -11,16 +11,16 @@
 
 namespace MeshLib {
 
-template <unsigned NNODES>
-TemplateEdge<NNODES>::TemplateEdge(Node* nodes[NNODES], unsigned value) :
+template <unsigned NNODES, FEMElemType::type FEMEDGETYPE>
+TemplateEdge<NNODES,FEMEDGETYPE>::TemplateEdge(Node* nodes[NNODES], unsigned value) :
 	Element(value)
 {
 	_nodes = nodes;
 	this->_length = this->computeVolume();
 }
 
-template <unsigned NNODES>
-TemplateEdge<NNODES>::TemplateEdge(const TemplateEdge<NNODES> &edge) :
+template <unsigned NNODES, FEMElemType::type FEMEDGETYPE>
+TemplateEdge<NNODES,FEMEDGETYPE>::TemplateEdge(const TemplateEdge<NNODES,FEMEDGETYPE> &edge) :
 	Element(edge.getValue())
 {
 	_nodes = new Node*[NNODES];
@@ -29,8 +29,8 @@ TemplateEdge<NNODES>::TemplateEdge(const TemplateEdge<NNODES> &edge) :
 	_length = edge.getLength();
 }
 
-template <unsigned NNODES>
-TemplateEdge<NNODES>::~TemplateEdge()
+template <unsigned NNODES, FEMElemType::type FEMEDGETYPE>
+TemplateEdge<NNODES,FEMEDGETYPE>::~TemplateEdge()
 {}
 
 } // namespace MeshLib
diff --git a/MeshLib/Elements/TemplateHex.h b/MeshLib/Elements/TemplateHex.h
index 0a67ccf42620ed8beaac1d02305598672b454cfd..3aab5c49738fc3335a977534184a7b6cb480ebad 100644
--- a/MeshLib/Elements/TemplateHex.h
+++ b/MeshLib/Elements/TemplateHex.h
@@ -13,6 +13,7 @@
 #ifndef TEMPLATEHEX_H_
 #define TEMPLATEHEX_H_
 
+#include "MshEnums.h"
 #include "Cell.h"
 
 namespace MeshLib {
@@ -42,15 +43,15 @@ namespace MeshLib {
  *
  * @endcode
  */
-template <unsigned NNODES>
+template <unsigned NNODES, FEMElemType::type FEMHEXTYPE>
 class TemplateHex : public Cell
 {
 public:
 	/// Constructor with an array of mesh nodes.
-	TemplateHex(Node* nodes[8], unsigned value = 0);
+	TemplateHex(Node* nodes[NNODES], unsigned value = 0);
 
 	/// Copy constructor
-	TemplateHex(const TemplateHex<NNODES> &hex);
+	TemplateHex(const TemplateHex &hex);
 
 	/// Destructor
 	virtual ~TemplateHex();
@@ -82,6 +83,12 @@ public:
 	 */
 	virtual MshElemType::type getGeoType() const { return MshElemType::HEXAHEDRON; }
 
+	/**
+	 * Method returns the FEM type of the element.
+	 * @return
+	 */
+	virtual FEMElemType::type getFEMType() const { return FEMHEXTYPE; }
+
 	/// Returns true if these two indices form an edge and false otherwise
 	bool isEdge(unsigned i, unsigned j) const;
 
diff --git a/MeshLib/Elements/TemplateHex.tpp b/MeshLib/Elements/TemplateHex.tpp
index cd05fdebcb073c8f36e5b0305cc51bf844474fd6..43d4d0b4475571e0db398babf4b1301321af6f96 100644
--- a/MeshLib/Elements/TemplateHex.tpp
+++ b/MeshLib/Elements/TemplateHex.tpp
@@ -18,8 +18,8 @@
 
 namespace MeshLib {
 
-template <unsigned NNODES>
-const unsigned TemplateHex<NNODES>::_face_nodes[6][4] =
+template <unsigned NNODES, FEMElemType::type FEMHEXTYPE>
+const unsigned TemplateHex<NNODES,FEMHEXTYPE>::_face_nodes[6][4] =
 {
 	{0, 3, 2, 1}, // Face 0
 	{0, 1, 5, 4}, // Face 1
@@ -29,8 +29,8 @@ const unsigned TemplateHex<NNODES>::_face_nodes[6][4] =
 	{4, 5, 6, 7}  // Face 5
 };
 
-template <unsigned NNODES>
-const unsigned TemplateHex<NNODES>::_edge_nodes[12][2] =
+template <unsigned NNODES, FEMElemType::type FEMHEXTYPE>
+const unsigned TemplateHex<NNODES,FEMHEXTYPE>::_edge_nodes[12][2] =
 {
 	{0, 1}, // Edge 0
 	{1, 2}, // Edge 1
@@ -46,8 +46,8 @@ const unsigned TemplateHex<NNODES>::_edge_nodes[12][2] =
 	{4, 7}  // Edge 11
 };
 
-template <unsigned NNODES>
-TemplateHex<NNODES>::TemplateHex(Node* nodes[8], unsigned value)
+template <unsigned NNODES, FEMElemType::type FEMHEXTYPE>
+TemplateHex<NNODES,FEMHEXTYPE>::TemplateHex(Node* nodes[8], unsigned value)
 	: Cell(value)
 {
 	_nodes = nodes;
@@ -59,8 +59,8 @@ TemplateHex<NNODES>::TemplateHex(Node* nodes[8], unsigned value)
 	this->_volume = this->computeVolume();
 }
 
-template <unsigned NNODES>
-TemplateHex<NNODES>::TemplateHex(const TemplateHex<NNODES> &hex)
+template <unsigned NNODES, FEMElemType::type FEMHEXTYPE>
+TemplateHex<NNODES,FEMHEXTYPE>::TemplateHex(const TemplateHex<NNODES,FEMHEXTYPE> &hex)
 	: Cell(hex.getValue())
 {
 	_nodes = new Node*[NNODES];
@@ -74,13 +74,13 @@ TemplateHex<NNODES>::TemplateHex(const TemplateHex<NNODES> &hex)
 	_volume = hex.getVolume();
 }
 
-template <unsigned NNODES>
-TemplateHex<NNODES>::~TemplateHex()
+template <unsigned NNODES, FEMElemType::type FEMHEXTYPE>
+TemplateHex<NNODES,FEMHEXTYPE>::~TemplateHex()
 {
 }
 
-template <unsigned NNODES>
-double TemplateHex<NNODES>::computeVolume()
+template <unsigned NNODES, FEMElemType::type FEMHEXTYPE>
+double TemplateHex<NNODES,FEMHEXTYPE>::computeVolume()
 {
 	return MathLib::calcTetrahedronVolume(_nodes[4]->getCoords(), _nodes[7]->getCoords(), _nodes[5]->getCoords(), _nodes[0]->getCoords())
 		 + MathLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[3]->getCoords(), _nodes[1]->getCoords(), _nodes[0]->getCoords())
@@ -90,8 +90,8 @@ double TemplateHex<NNODES>::computeVolume()
 		 + MathLib::calcTetrahedronVolume(_nodes[3]->getCoords(), _nodes[7]->getCoords(), _nodes[5]->getCoords(), _nodes[2]->getCoords());
 }
 
-template <unsigned NNODES>
-const Element* TemplateHex<NNODES>::getFace(unsigned i) const
+template <unsigned NNODES, FEMElemType::type FEMHEXTYPE>
+const Element* TemplateHex<NNODES,FEMHEXTYPE>::getFace(unsigned i) const
 {
 	if (i<this->getNFaces())
 	{
@@ -105,8 +105,8 @@ const Element* TemplateHex<NNODES>::getFace(unsigned i) const
 	return NULL;
 }
 
-template <unsigned NNODES>
-bool TemplateHex<NNODES>::isEdge(unsigned idx1, unsigned idx2) const
+template <unsigned NNODES, FEMElemType::type FEMHEXTYPE>
+bool TemplateHex<NNODES,FEMHEXTYPE>::isEdge(unsigned idx1, unsigned idx2) const
 {
 	for (unsigned i(0); i<12; i++)
 	{
@@ -116,14 +116,14 @@ bool TemplateHex<NNODES>::isEdge(unsigned idx1, unsigned idx2) const
 	return false;
 }
 
-template <unsigned NNODES>
-Element* TemplateHex<NNODES>::clone() const
+template <unsigned NNODES, FEMElemType::type FEMHEXTYPE>
+Element* TemplateHex<NNODES,FEMHEXTYPE>::clone() const
 {
-	return new TemplateHex<NNODES>(*this);
+	return new TemplateHex<NNODES,FEMHEXTYPE>(*this);
 }
 
-template <unsigned NNODES>
-unsigned TemplateHex<NNODES>::identifyFace(Node* nodes[3]) const
+template <unsigned NNODES, FEMElemType::type FEMHEXTYPE>
+unsigned TemplateHex<NNODES,FEMHEXTYPE>::identifyFace(Node* nodes[3]) const
 {
 	for (unsigned i=0; i<6; i++)
 	{
@@ -138,8 +138,8 @@ unsigned TemplateHex<NNODES>::identifyFace(Node* nodes[3]) const
 	return std::numeric_limits<unsigned>::max();
 }
 
-template <unsigned NNODES>
-Element* TemplateHex<NNODES>::reviseElement() const
+template <unsigned NNODES, FEMElemType::type FEMHEXTYPE>
+Element* TemplateHex<NNODES,FEMHEXTYPE>::reviseElement() const
 {
 	std::vector<size_t> collapsed_edges;
 	for (size_t edge(0); edge<getNEdges(); edge++) {
@@ -149,7 +149,7 @@ Element* TemplateHex<NNODES>::reviseElement() const
 	}
 
 	if (collapsed_edges.size() == 1) {
-		std::cerr << "[TemplateHex<NNODES>::reviseElement()] collapsing of one edge in hexahedron not handled" << std::endl;
+		std::cerr << "[TemplateHex<NNODES,FEMHEXTYPE>::reviseElement()] collapsing of one edge in hexahedron not handled" << std::endl;
 		return NULL;
 	}
 
diff --git a/MeshLib/Elements/TemplatePrism.h b/MeshLib/Elements/TemplatePrism.h
index f9522d505e8024028da499c8fb70db1b2118b794..ccc5614f9a7b1ff63a02f91e804cb8fe45fa4157 100644
--- a/MeshLib/Elements/TemplatePrism.h
+++ b/MeshLib/Elements/TemplatePrism.h
@@ -13,6 +13,7 @@
 #ifndef TEMPLATEPRISM_H_
 #define TEMPLATEPRISM_H_
 
+#include "MshEnums.h"
 #include "Cell.h"
 
 namespace MeshLib {
@@ -40,7 +41,7 @@ namespace MeshLib {
  *
  * @endcode
  */
-template <unsigned NNODES>
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
 class TemplatePrism : public Cell
 {
 public:
@@ -48,7 +49,7 @@ public:
 	TemplatePrism(Node* nodes[6], unsigned value = 0);
 
 	/// Copy constructor
-	TemplatePrism(const TemplatePrism<NNODES> &prism);
+	TemplatePrism(const TemplatePrism &prism);
 
 	/// Destructor
 	virtual ~TemplatePrism();
@@ -80,6 +81,12 @@ public:
 	 */
 	virtual MshElemType::type getGeoType() const { return MshElemType::PRISM; }
 
+	/**
+	 * Get the type of the element in context of the finite element method.
+	 * @return a value of the enum FEMElemType::type
+	 */
+	virtual FEMElemType::type getFEMType() const { return FEMPRISMTYPE; };
+
 	/// Returns true if these two indeces form an edge and false otherwise
 	bool isEdge(unsigned i, unsigned j) const;
 
diff --git a/MeshLib/Elements/TemplatePrism.tpp b/MeshLib/Elements/TemplatePrism.tpp
index cf1534cb586cc46dc36c39531d66a9c59ece3691..c640fb81d83f3dbdbe81f3b124487b28c290c3f6 100644
--- a/MeshLib/Elements/TemplatePrism.tpp
+++ b/MeshLib/Elements/TemplatePrism.tpp
@@ -19,8 +19,8 @@
 
 namespace MeshLib {
 
-template <unsigned NNODES>
-const unsigned TemplatePrism<NNODES>::_face_nodes[5][4] =
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
+const unsigned TemplatePrism<NNODES,FEMPRISMTYPE>::_face_nodes[5][4] =
 {
 	{0, 2, 1, 99}, // Face 0
 	{0, 1, 4,  3}, // Face 1
@@ -29,8 +29,8 @@ const unsigned TemplatePrism<NNODES>::_face_nodes[5][4] =
 	{3, 4, 5, 99}  // Face 4
 };
 
-template <unsigned NNODES>
-const unsigned TemplatePrism<NNODES>::_edge_nodes[9][2] =
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
+const unsigned TemplatePrism<NNODES,FEMPRISMTYPE>::_edge_nodes[9][2] =
 {
 	{0, 1}, // Edge 0
 	{1, 2}, // Edge 1
@@ -43,11 +43,11 @@ const unsigned TemplatePrism<NNODES>::_edge_nodes[9][2] =
 	{3, 5}  // Edge 8
 };
 
-template <unsigned NNODES>
-const unsigned TemplatePrism<NNODES>::_n_face_nodes[5] = { 3, 4, 4, 4, 3 };
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
+const unsigned TemplatePrism<NNODES,FEMPRISMTYPE>::_n_face_nodes[5] = { 3, 4, 4, 4, 3 };
 
-template <unsigned NNODES>
-TemplatePrism<NNODES>::TemplatePrism(Node* nodes[6], unsigned value)
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
+TemplatePrism<NNODES,FEMPRISMTYPE>::TemplatePrism(Node* nodes[6], unsigned value)
 	: Cell(value)
 {
 	_nodes = nodes;
@@ -57,8 +57,8 @@ TemplatePrism<NNODES>::TemplatePrism(Node* nodes[6], unsigned value)
 	this->_volume = this->computeVolume();
 }
 
-template <unsigned NNODES>
-TemplatePrism<NNODES>::TemplatePrism(const TemplatePrism<NNODES> &prism)
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
+TemplatePrism<NNODES,FEMPRISMTYPE>::TemplatePrism(const TemplatePrism<NNODES,FEMPRISMTYPE> &prism)
 	: Cell(prism.getValue())
 {
 	_nodes = new Node*[NNODES];
@@ -72,21 +72,21 @@ TemplatePrism<NNODES>::TemplatePrism(const TemplatePrism<NNODES> &prism)
 	_volume = prism.getVolume();
 }
 
-template <unsigned NNODES>
-TemplatePrism<NNODES>::~TemplatePrism()
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
+TemplatePrism<NNODES,FEMPRISMTYPE>::~TemplatePrism()
 {
 }
 
-template <unsigned NNODES>
-double TemplatePrism<NNODES>::computeVolume()
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
+double TemplatePrism<NNODES,FEMPRISMTYPE>::computeVolume()
 {
 	return MathLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords())
 		 + MathLib::calcTetrahedronVolume(_nodes[1]->getCoords(), _nodes[4]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords())
 		 + MathLib::calcTetrahedronVolume(_nodes[2]->getCoords(), _nodes[4]->getCoords(), _nodes[5]->getCoords(), _nodes[3]->getCoords());
 }
 
-template <unsigned NNODES>
-const Element* TemplatePrism<NNODES>::getFace(unsigned i) const
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
+const Element* TemplatePrism<NNODES,FEMPRISMTYPE>::getFace(unsigned i) const
 {
 	if (i<this->getNFaces())
 	{
@@ -104,8 +104,8 @@ const Element* TemplatePrism<NNODES>::getFace(unsigned i) const
 	return NULL;
 }
 
-template <unsigned NNODES>
-unsigned TemplatePrism<NNODES>::getNFaceNodes(unsigned i) const
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
+unsigned TemplatePrism<NNODES,FEMPRISMTYPE>::getNFaceNodes(unsigned i) const
 {
 	if (i<5)
 		return _n_face_nodes[i];
@@ -113,8 +113,8 @@ unsigned TemplatePrism<NNODES>::getNFaceNodes(unsigned i) const
 	return 0;
 }
 
-template <unsigned NNODES>
-bool TemplatePrism<NNODES>::isEdge(unsigned idx1, unsigned idx2) const
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
+bool TemplatePrism<NNODES,FEMPRISMTYPE>::isEdge(unsigned idx1, unsigned idx2) const
 {
 	for (unsigned i(0); i<9; i++)
 	{
@@ -124,14 +124,14 @@ bool TemplatePrism<NNODES>::isEdge(unsigned idx1, unsigned idx2) const
 	return false;
 }
 
-template <unsigned NNODES>
-Element* TemplatePrism<NNODES>::clone() const
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
+Element* TemplatePrism<NNODES,FEMPRISMTYPE>::clone() const
 {
-	return new TemplatePrism<NNODES>(*this);
+	return new TemplatePrism<NNODES,FEMPRISMTYPE>(*this);
 }
 
-template <unsigned NNODES>
-unsigned TemplatePrism<NNODES>::identifyFace(Node* nodes[3]) const
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
+unsigned TemplatePrism<NNODES,FEMPRISMTYPE>::identifyFace(Node* nodes[3]) const
 {
 	for (unsigned i=0; i<5; i++)
 	{
@@ -146,8 +146,8 @@ unsigned TemplatePrism<NNODES>::identifyFace(Node* nodes[3]) const
 	return std::numeric_limits<unsigned>::max();
 }
 
-template <unsigned NNODES>
-Element* TemplatePrism<NNODES>::reviseElement() const
+template <unsigned NNODES, FEMElemType::type FEMPRISMTYPE>
+Element* TemplatePrism<NNODES,FEMPRISMTYPE>::reviseElement() const
 {
 	// try to create Pyramid
 	if (_nodes[_edge_nodes[3][0]] == _nodes[_edge_nodes[3][1]]) {
diff --git a/MeshLib/Elements/TemplatePyramid.h b/MeshLib/Elements/TemplatePyramid.h
index 58503f9d64b5123daf82b50984a87d7649015aee..03abf05b734e726e35f64ffa54251ddcaddc6221 100644
--- a/MeshLib/Elements/TemplatePyramid.h
+++ b/MeshLib/Elements/TemplatePyramid.h
@@ -13,6 +13,7 @@
 #ifndef TEMPLATEPYRAMID_H_
 #define TEMPLATEPYRAMID_H_
 
+#include "MshEnums.h"
 #include "Cell.h"
 
 namespace MeshLib {
@@ -38,7 +39,7 @@ namespace MeshLib {
  *        0
  * @endcode
  */
-template <unsigned NNODES>
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
 class TemplatePyramid : public Cell
 {
 public:
@@ -46,7 +47,7 @@ public:
 	TemplatePyramid(Node* nodes[NNODES], unsigned value = 0);
 
 	/// Copy constructor
-	TemplatePyramid(const TemplatePyramid<NNODES> &pyramid);
+	TemplatePyramid(const TemplatePyramid &pyramid);
 
 	/// Destructor
 	virtual ~TemplatePyramid();
@@ -78,6 +79,12 @@ public:
 	 */
 	virtual MshElemType::type getGeoType() const { return MshElemType::PYRAMID; }
 
+	/**
+	 * Get the type of the element in context of the finite element method.
+	 * @return a value of the enum FEMElemType::type
+	 */
+	virtual FEMElemType::type getFEMType() const { return FEMPYRAMIDTYPE; }
+
 	/// Returns true if these two indeces form an edge and false otherwise
 	bool isEdge(unsigned i, unsigned j) const;
 
diff --git a/MeshLib/Elements/TemplatePyramid.tpp b/MeshLib/Elements/TemplatePyramid.tpp
index 1ddcb4c5b289be7fb2be76a8cb93c9e15130877c..cd3646e874131fe6e709ae5ebf458f494773062e 100644
--- a/MeshLib/Elements/TemplatePyramid.tpp
+++ b/MeshLib/Elements/TemplatePyramid.tpp
@@ -19,8 +19,8 @@
 
 namespace MeshLib {
 
-template <unsigned NNODES>
-const unsigned TemplatePyramid<NNODES>::_face_nodes[5][4] =
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
+const unsigned TemplatePyramid<NNODES,FEMPYRAMIDTYPE>::_face_nodes[5][4] =
 {
 	{0, 1, 4, 99}, // Face 0
 	{1, 2, 4, 99}, // Face 1
@@ -29,8 +29,8 @@ const unsigned TemplatePyramid<NNODES>::_face_nodes[5][4] =
 	{0, 3, 2,  1}  // Face 4
 };
 
-template <unsigned NNODES>
-const unsigned TemplatePyramid<NNODES>::_edge_nodes[8][2] =
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
+const unsigned TemplatePyramid<NNODES,FEMPYRAMIDTYPE>::_edge_nodes[8][2] =
 {
 	{0, 1}, // Edge 0
 	{1, 2}, // Edge 1
@@ -42,11 +42,11 @@ const unsigned TemplatePyramid<NNODES>::_edge_nodes[8][2] =
 	{3, 4}  // Edge 7
 };
 
-template <unsigned NNODES>
-const unsigned TemplatePyramid<NNODES>::_n_face_nodes[5] = { 3, 3, 3, 3, 4 };
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
+const unsigned TemplatePyramid<NNODES,FEMPYRAMIDTYPE>::_n_face_nodes[5] = { 3, 3, 3, 3, 4 };
 
-template <unsigned NNODES>
-TemplatePyramid<NNODES>::TemplatePyramid(Node* nodes[NNODES], unsigned value)
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
+TemplatePyramid<NNODES,FEMPYRAMIDTYPE>::TemplatePyramid(Node* nodes[NNODES], unsigned value)
 	: Cell(value)
 {
 	_nodes = nodes;
@@ -58,8 +58,8 @@ TemplatePyramid<NNODES>::TemplatePyramid(Node* nodes[NNODES], unsigned value)
 	this->_volume = this->computeVolume();
 }
 
-template <unsigned NNODES>
-TemplatePyramid<NNODES>::TemplatePyramid(const TemplatePyramid<NNODES> &pyramid)
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
+TemplatePyramid<NNODES,FEMPYRAMIDTYPE>::TemplatePyramid(const TemplatePyramid<NNODES,FEMPYRAMIDTYPE> &pyramid)
 	: Cell(pyramid.getValue())
 {
 	_nodes = new Node*[NNODES];
@@ -75,20 +75,20 @@ TemplatePyramid<NNODES>::TemplatePyramid(const TemplatePyramid<NNODES> &pyramid)
 	_volume = pyramid.getVolume();
 }
 
-template <unsigned NNODES>
-TemplatePyramid<NNODES>::~TemplatePyramid()
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
+TemplatePyramid<NNODES,FEMPYRAMIDTYPE>::~TemplatePyramid()
 {
 }
 
-template <unsigned NNODES>
-double TemplatePyramid<NNODES>::computeVolume()
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
+double TemplatePyramid<NNODES,FEMPYRAMIDTYPE>::computeVolume()
 {
 	return MathLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[4]->getCoords())
 		 + MathLib::calcTetrahedronVolume(_nodes[2]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords(), _nodes[4]->getCoords());
 }
 
-template <unsigned NNODES>
-const Element* TemplatePyramid<NNODES>::getFace(unsigned i) const
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
+const Element* TemplatePyramid<NNODES,FEMPYRAMIDTYPE>::getFace(unsigned i) const
 {
 	if (i<this->getNFaces())
 	{
@@ -106,8 +106,8 @@ const Element* TemplatePyramid<NNODES>::getFace(unsigned i) const
 	return NULL;
 }
 
-template <unsigned NNODES>
-unsigned TemplatePyramid<NNODES>::getNFaceNodes(unsigned i) const
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
+unsigned TemplatePyramid<NNODES,FEMPYRAMIDTYPE>::getNFaceNodes(unsigned i) const
 {
 	if (i<5)
 		return _n_face_nodes[i];
@@ -115,8 +115,8 @@ unsigned TemplatePyramid<NNODES>::getNFaceNodes(unsigned i) const
 	return 0;
 }
 
-template <unsigned NNODES>
-bool TemplatePyramid<NNODES>::isEdge(unsigned idx1, unsigned idx2) const
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
+bool TemplatePyramid<NNODES,FEMPYRAMIDTYPE>::isEdge(unsigned idx1, unsigned idx2) const
 {
 	for (unsigned i(0); i<8; i++)
 	{
@@ -126,14 +126,14 @@ bool TemplatePyramid<NNODES>::isEdge(unsigned idx1, unsigned idx2) const
 	return false;
 }
 
-template <unsigned NNODES>
-Element* TemplatePyramid<NNODES>::clone() const
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
+Element* TemplatePyramid<NNODES,FEMPYRAMIDTYPE>::clone() const
 {
-	return new TemplatePyramid<NNODES>(*this);
+	return new TemplatePyramid<NNODES,FEMPYRAMIDTYPE>(*this);
 }
 
-template <unsigned NNODES>
-unsigned TemplatePyramid<NNODES>::identifyFace(Node* nodes[3]) const
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
+unsigned TemplatePyramid<NNODES,FEMPYRAMIDTYPE>::identifyFace(Node* nodes[3]) const
 {
 	for (unsigned i=0; i<5; i++)
 	{
@@ -148,8 +148,8 @@ unsigned TemplatePyramid<NNODES>::identifyFace(Node* nodes[3]) const
 	return std::numeric_limits<unsigned>::max();
 }
 
-template <unsigned NNODES>
-Element* TemplatePyramid<NNODES>::reviseElement() const
+template <unsigned NNODES, FEMElemType::type FEMPYRAMIDTYPE>
+Element* TemplatePyramid<NNODES,FEMPYRAMIDTYPE>::reviseElement() const
 {
 	// try to create tetrahedron
 	if (_nodes[_edge_nodes[0][0]] == _nodes[_edge_nodes[0][1]]
diff --git a/MeshLib/Elements/TemplateQuad.h b/MeshLib/Elements/TemplateQuad.h
index f438c86ebb6946454551457f84a3898e2a92b33f..5acef9471e1000714888885dbfd3db423df849e3 100644
--- a/MeshLib/Elements/TemplateQuad.h
+++ b/MeshLib/Elements/TemplateQuad.h
@@ -5,7 +5,7 @@
  *              http://www.opengeosys.org/project/license
  *
  *
- * \file Quad.h
+ * \file TemplateQuad.h
  *
  * Created on 2012-05-02 by Karsten Rink
  */
@@ -13,6 +13,7 @@
 #ifndef TEMPLATEQUAD_H_
 #define TEMPLATEQUAD_H_
 
+#include "MshEnums.h"
 #include "Face.h"
 
 namespace MeshLib {
@@ -32,7 +33,7 @@ namespace MeshLib {
  *              0
  * @endcode
  */
-template <unsigned NNODES>
+template <unsigned NNODES, FEMElemType::type FEMQUADTYPE>
 class TemplateQuad : public Face
 {
 public:
@@ -63,6 +64,12 @@ public:
 	 */
 	virtual MshElemType::type getGeoType() const { return MshElemType::QUAD; }
 
+	/**
+	 * Get the type of the element in context of the finite element method.
+	 * @return a value of the enum FEMElemType::type
+	 */
+	virtual FEMElemType::type getFEMType() const { return FEMQUADTYPE; }
+
 	/// Returns true if these two indeces form an edge and false otherwise
 	bool isEdge(unsigned i, unsigned j) const;
 
@@ -98,8 +105,8 @@ protected:
 
 }; /* class */
 
-template <unsigned NNODES>
-const unsigned TemplateQuad<NNODES>::_edge_nodes[4][2] =
+template <unsigned NNODES, FEMElemType::type FEMQUADTYPE>
+const unsigned TemplateQuad<NNODES, FEMQUADTYPE>::_edge_nodes[4][2] =
 {
 	{0, 1}, // Edge 0
 	{1, 2}, // Edge 1
diff --git a/MeshLib/Elements/TemplateQuad.tpp b/MeshLib/Elements/TemplateQuad.tpp
index f9afea0632b6d44e846b58d862e696edc1454099..eeb46b579d35565e159009d36f1c3fa9cfbef83b 100644
--- a/MeshLib/Elements/TemplateQuad.tpp
+++ b/MeshLib/Elements/TemplateQuad.tpp
@@ -17,8 +17,8 @@
 
 namespace MeshLib {
 
-template <unsigned NNODES>
-TemplateQuad<NNODES>::TemplateQuad(Node* nodes[NNODES], unsigned value)
+template <unsigned NNODES, FEMElemType::type FEMQUADTYPE>
+TemplateQuad<NNODES,FEMQUADTYPE>::TemplateQuad(Node* nodes[NNODES], unsigned value)
 	: Face(value)
 {
 	_nodes = nodes;
@@ -30,8 +30,8 @@ TemplateQuad<NNODES>::TemplateQuad(Node* nodes[NNODES], unsigned value)
 	this->_area = this->computeVolume();
 }
 
-template <unsigned NNODES>
-TemplateQuad<NNODES>::TemplateQuad(const TemplateQuad<NNODES> &quad)
+template <unsigned NNODES, FEMElemType::type FEMQUADTYPE>
+TemplateQuad<NNODES,FEMQUADTYPE>::TemplateQuad(const TemplateQuad<NNODES,FEMQUADTYPE> &quad)
 	: Face(quad.getValue())
 {
 	_nodes = new Node*[NNODES];
@@ -47,20 +47,20 @@ TemplateQuad<NNODES>::TemplateQuad(const TemplateQuad<NNODES> &quad)
 	_area = quad.getArea();
 }
 
-template <unsigned NNODES>
-TemplateQuad<NNODES>::~TemplateQuad()
+template <unsigned NNODES, FEMElemType::type FEMQUADTYPE>
+TemplateQuad<NNODES,FEMQUADTYPE>::~TemplateQuad()
 {
 }
 
-template <unsigned NNODES>
-double TemplateQuad<NNODES>::computeVolume()
+template <unsigned NNODES, FEMElemType::type FEMQUADTYPE>
+double TemplateQuad<NNODES,FEMQUADTYPE>::computeVolume()
 {
 	return MathLib::calcTriangleArea(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords())
          + MathLib::calcTriangleArea(_nodes[2]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords());
 }
 
-template <unsigned NNODES>
-bool TemplateQuad<NNODES>::isEdge(unsigned idx1, unsigned idx2) const
+template <unsigned NNODES, FEMElemType::type FEMQUADTYPE>
+bool TemplateQuad<NNODES,FEMQUADTYPE>::isEdge(unsigned idx1, unsigned idx2) const
 {
 	for (unsigned i(0); i<4; i++)
 	{
@@ -70,14 +70,14 @@ bool TemplateQuad<NNODES>::isEdge(unsigned idx1, unsigned idx2) const
 	return false;
 }
 
-template <unsigned NNODES>
-Element* TemplateQuad<NNODES>::clone() const
+template <unsigned NNODES, FEMElemType::type FEMQUADTYPE>
+Element* TemplateQuad<NNODES,FEMQUADTYPE>::clone() const
 {
 	return new TemplateQuad(*this);
 }
 
-template <unsigned NNODES>
-unsigned TemplateQuad<NNODES>::identifyFace(Node* nodes[3]) const
+template <unsigned NNODES, FEMElemType::type FEMQUADTYPE>
+unsigned TemplateQuad<NNODES,FEMQUADTYPE>::identifyFace(Node* nodes[3]) const
 {
 	for (unsigned i=0; i<4; i++)
 	{
@@ -92,8 +92,8 @@ unsigned TemplateQuad<NNODES>::identifyFace(Node* nodes[3]) const
 	return std::numeric_limits<unsigned>::max();
 }
 
-template <unsigned NNODES>
-Element* TemplateQuad<NNODES>::reviseElement() const
+template <unsigned NNODES, FEMElemType::type FEMQUADTYPE>
+Element* TemplateQuad<NNODES,FEMQUADTYPE>::reviseElement() const
 {
 	if (_nodes[0] == _nodes[1] || _nodes[1] == _nodes[2]) {
 		MeshLib::Node** tri_nodes = new MeshLib::Node*[3];
diff --git a/MeshLib/Elements/TemplateTet.h b/MeshLib/Elements/TemplateTet.h
index 802568e56efbfb3803648d300ad62addb8769e15..bc03f5012d8d66c35a4ff8501c18e4ccb10e47bc 100644
--- a/MeshLib/Elements/TemplateTet.h
+++ b/MeshLib/Elements/TemplateTet.h
@@ -13,6 +13,7 @@
 #ifndef TEMPLATETET_H_
 #define TEMPLATETET_H_
 
+#include "MshEnums.h"
 #include "Cell.h"
 
 namespace MeshLib {
@@ -37,7 +38,7 @@ namespace MeshLib {
  *
  * @endcode
  */
-template <unsigned NNODES>
+template <unsigned NNODES, FEMElemType::type FEMTETTYPE>
 class TemplateTet : public Cell
 {
 public:
@@ -45,7 +46,7 @@ public:
 	TemplateTet(Node* nodes[NNODES], unsigned value = 0);
 
 	/// Copy constructor
-	TemplateTet(const TemplateTet<NNODES> &tet);
+	TemplateTet(const TemplateTet &tet);
 
 	/// Destructor
 	virtual ~TemplateTet();
@@ -77,6 +78,12 @@ public:
 	 */
 	virtual MshElemType::type getGeoType() const { return MshElemType::TETRAHEDRON; }
 
+	/**
+	 * Get the type of the element in context of the finite element method.
+	 * @return a value of the enum FEMElemType::type
+	 */
+	virtual FEMElemType::type getFEMType() const { return FEMTETTYPE; }
+
 	/// Returns true if these two indeces form an edge and false otherwise
 	bool isEdge(unsigned i, unsigned j) const;
 
diff --git a/MeshLib/Elements/TemplateTet.tpp b/MeshLib/Elements/TemplateTet.tpp
index 4e26e01b53e029a5d80ba2d780d4bada2dd2dfce..e5f1245091bc63685c5df03b6ed014b48d8c1ea7 100644
--- a/MeshLib/Elements/TemplateTet.tpp
+++ b/MeshLib/Elements/TemplateTet.tpp
@@ -17,8 +17,8 @@
 
 namespace MeshLib {
 
-template <unsigned NNODES>
-const unsigned TemplateTet<NNODES>::_face_nodes[4][3] =
+template <unsigned NNODES, FEMElemType::type FEMTETTYPE>
+const unsigned TemplateTet<NNODES,FEMTETTYPE>::_face_nodes[4][3] =
 {
 	{0, 2, 1}, // Face 0
 	{0, 1, 3}, // Face 1
@@ -26,8 +26,8 @@ const unsigned TemplateTet<NNODES>::_face_nodes[4][3] =
 	{2, 0, 3}  // Face 3
 };
 
-template <unsigned NNODES>
-const unsigned TemplateTet<NNODES>::_edge_nodes[6][2] =
+template <unsigned NNODES, FEMElemType::type FEMTETTYPE>
+const unsigned TemplateTet<NNODES,FEMTETTYPE>::_edge_nodes[6][2] =
 {
 	{0, 1}, // Edge 0
 	{1, 2}, // Edge 1
@@ -37,8 +37,8 @@ const unsigned TemplateTet<NNODES>::_edge_nodes[6][2] =
 	{2, 3}  // Edge 5
 };
 
-template <unsigned NNODES>
-TemplateTet<NNODES>::TemplateTet(Node* nodes[NNODES], unsigned value)
+template <unsigned NNODES, FEMElemType::type FEMTETTYPE>
+TemplateTet<NNODES,FEMTETTYPE>::TemplateTet(Node* nodes[NNODES], unsigned value)
 	: Cell(value)
 {
 	_nodes = nodes;
@@ -50,8 +50,8 @@ TemplateTet<NNODES>::TemplateTet(Node* nodes[NNODES], unsigned value)
 	this->_volume = this->computeVolume();
 }
 
-template <unsigned NNODES>
-TemplateTet<NNODES>::TemplateTet(const TemplateTet<NNODES> &tet)
+template <unsigned NNODES, FEMElemType::type FEMTETTYPE>
+TemplateTet<NNODES,FEMTETTYPE>::TemplateTet(const TemplateTet<NNODES,FEMTETTYPE> &tet)
 	: Cell(tet.getValue())
 {
 	_nodes = new Node*[NNODES];
@@ -68,19 +68,19 @@ TemplateTet<NNODES>::TemplateTet(const TemplateTet<NNODES> &tet)
 	_volume = tet.getVolume();
 }
 
-template <unsigned NNODES>
-TemplateTet<NNODES>::~TemplateTet()
+template <unsigned NNODES, FEMElemType::type FEMTETTYPE>
+TemplateTet<NNODES,FEMTETTYPE>::~TemplateTet()
 {
 }
 
-template <unsigned NNODES>
-double TemplateTet<NNODES>::computeVolume()
+template <unsigned NNODES, FEMElemType::type FEMTETTYPE>
+double TemplateTet<NNODES,FEMTETTYPE>::computeVolume()
 {
 	return MathLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords());
 }
 
-template <unsigned NNODES>
-const Element* TemplateTet<NNODES>::getFace(unsigned i) const
+template <unsigned NNODES, FEMElemType::type FEMTETTYPE>
+const Element* TemplateTet<NNODES,FEMTETTYPE>::getFace(unsigned i) const
 {
 	if (i<this->getNFaces())
 	{
@@ -94,8 +94,8 @@ const Element* TemplateTet<NNODES>::getFace(unsigned i) const
 	return NULL;
 }
 
-template <unsigned NNODES>
-bool TemplateTet<NNODES>::isEdge(unsigned idx1, unsigned idx2) const
+template <unsigned NNODES, FEMElemType::type FEMTETTYPE>
+bool TemplateTet<NNODES,FEMTETTYPE>::isEdge(unsigned idx1, unsigned idx2) const
 {
 	for (unsigned i(0); i<6; i++)
 	{
@@ -105,14 +105,14 @@ bool TemplateTet<NNODES>::isEdge(unsigned idx1, unsigned idx2) const
 	return false;
 }
 
-template <unsigned NNODES>
-Element* TemplateTet<NNODES>::clone() const
+template <unsigned NNODES, FEMElemType::type FEMTETTYPE>
+Element* TemplateTet<NNODES,FEMTETTYPE>::clone() const
 {
-	return new TemplateTet<NNODES>(*this);
+	return new TemplateTet<NNODES,FEMTETTYPE>(*this);
 }
 
-template <unsigned NNODES>
-unsigned TemplateTet<NNODES>::identifyFace(Node* nodes[3]) const
+template <unsigned NNODES, FEMElemType::type FEMTETTYPE>
+unsigned TemplateTet<NNODES,FEMTETTYPE>::identifyFace(Node* nodes[3]) const
 {
 	for (unsigned i=0; i<4; i++)
 	{
@@ -127,8 +127,8 @@ unsigned TemplateTet<NNODES>::identifyFace(Node* nodes[3]) const
 	return std::numeric_limits<unsigned>::max();
 }
 
-template <unsigned NNODES>
-Element* TemplateTet<NNODES>::reviseElement() const
+template <unsigned NNODES, FEMElemType::type FEMTETTYPE>
+Element* TemplateTet<NNODES,FEMTETTYPE>::reviseElement() const
 {
 	if (_nodes[0] == _nodes[1] || _nodes[1] == _nodes[2]) {
 		MeshLib::Node** tri_nodes = new MeshLib::Node*[3];
diff --git a/MeshLib/Elements/TemplateTri.h b/MeshLib/Elements/TemplateTri.h
index 50ba490d0fe1918b4b4d65e86c781f00de4b3f0e..12ccd1ea1cdcad05ba0db15eff3300d9922da0ef 100644
--- a/MeshLib/Elements/TemplateTri.h
+++ b/MeshLib/Elements/TemplateTri.h
@@ -16,6 +16,7 @@
 #include "Edge.h"
 #include "Node.h"
 #include "Face.h"
+#include "MshEnums.h"
 
 #include "MathTools.h"
 
@@ -39,7 +40,7 @@ namespace MeshLib {
  *
  * @endcode
  */
-template <unsigned NNODES>
+template <unsigned NNODES, FEMElemType::type FEMTRITYPE>
 class TemplateTri : public Face
 {
 public:
@@ -47,7 +48,7 @@ public:
 	TemplateTri(Node* nodes[NNODES], unsigned value = 0);
 
 	/// Copy constructor
-	TemplateTri(const TemplateTri<NNODES> &tri);
+	TemplateTri(const TemplateTri &tri);
 
 	/// Destructor
 	virtual ~TemplateTri();
@@ -70,6 +71,12 @@ public:
 	 */
 	virtual MshElemType::type getGeoType() const { return MshElemType::TRIANGLE; }
 
+	/**
+	 * Get the type of the element in context of the finite element method.
+	 * @return a value of the enum FEMElemType::type
+	 */
+	virtual FEMElemType::type getFEMType() const { return FEMTRITYPE; }
+
 	/// Returns true if these two indices form an edge and false otherwise
 	bool isEdge(unsigned idx1, unsigned idx2) const;
 
@@ -79,7 +86,7 @@ public:
 	 */
 	virtual Element* clone() const
 	{
-		return new TemplateTri<NNODES>(*this);
+		return new TemplateTri<NNODES,FEMTRITYPE>(*this);
 	}
 
 
@@ -113,8 +120,8 @@ protected:
 	static const unsigned _edge_nodes[3][2];
 }; /* class */
 
-template <unsigned NNODES>
-const unsigned TemplateTri<NNODES>::_edge_nodes[3][2] = {
+template <unsigned NNODES, FEMElemType::type FEMTRITYPE>
+const unsigned TemplateTri<NNODES,FEMTRITYPE>::_edge_nodes[3][2] = {
 		{0, 1}, // Edge 0
 		{1, 2}, // Edge 1
 		{0, 2}  // Edge 2
diff --git a/MeshLib/Elements/TemplateTri.tpp b/MeshLib/Elements/TemplateTri.tpp
index 0138e0bd9ef53ff2b0375e905567ef0a614a8e0d..ab1d48e8baac97c53a1ad757d6dc5bef956f731e 100644
--- a/MeshLib/Elements/TemplateTri.tpp
+++ b/MeshLib/Elements/TemplateTri.tpp
@@ -11,8 +11,8 @@
 
 namespace MeshLib {
 
-template <unsigned NNODES>
-TemplateTri<NNODES>::TemplateTri(Node* nodes[NNODES], unsigned value) :
+template <unsigned NNODES, FEMElemType::type FEMTRITYPE>
+TemplateTri<NNODES,FEMTRITYPE>::TemplateTri(Node* nodes[NNODES], unsigned value) :
 	Face(value)
 {
 	_nodes = nodes;
@@ -22,8 +22,8 @@ TemplateTri<NNODES>::TemplateTri(Node* nodes[NNODES], unsigned value) :
 	this->_area = this->computeVolume();
 }
 
-template <unsigned NNODES>
-TemplateTri<NNODES>::TemplateTri(const TemplateTri<NNODES> &tri) :
+template <unsigned NNODES, FEMElemType::type FEMTRITYPE>
+TemplateTri<NNODES,FEMTRITYPE>::TemplateTri(const TemplateTri<NNODES,FEMTRITYPE> &tri) :
 	Face(tri.getValue())
 {
 	_nodes = new Node*[NNODES];
@@ -40,12 +40,12 @@ TemplateTri<NNODES>::TemplateTri(const TemplateTri<NNODES> &tri) :
 	_area = tri.getArea();
 }
 
-template <unsigned NNODES>
-TemplateTri<NNODES>::~TemplateTri()
+template <unsigned NNODES, FEMElemType::type FEMTRITYPE>
+TemplateTri<NNODES,FEMTRITYPE>::~TemplateTri()
 {}
 
-template <unsigned NNODES>
-bool TemplateTri<NNODES>::isEdge(unsigned idx1, unsigned idx2) const
+template <unsigned NNODES, FEMElemType::type FEMTRITYPE>
+bool TemplateTri<NNODES,FEMTRITYPE>::isEdge(unsigned idx1, unsigned idx2) const
 {
 	for (unsigned i(0); i<3; i++)
 	{
@@ -55,8 +55,8 @@ bool TemplateTri<NNODES>::isEdge(unsigned idx1, unsigned idx2) const
 	return false;
 }
 
-template <unsigned NNODES>
-Element* TemplateTri<NNODES>::reviseElement() const
+template <unsigned NNODES, FEMElemType::type FEMTRITYPE>
+Element* TemplateTri<NNODES,FEMTRITYPE>::reviseElement() const
 {
 	// try to create an edge
 	if (_nodes[0] == _nodes[1] || _nodes[1] == _nodes[2]) {
@@ -76,8 +76,8 @@ Element* TemplateTri<NNODES>::reviseElement() const
 	return NULL;
 }
 
-template <unsigned NNODES>
-unsigned TemplateTri<NNODES>::identifyFace(Node* nodes[3]) const
+template <unsigned NNODES, FEMElemType::type FEMTRITYPE>
+unsigned TemplateTri<NNODES,FEMTRITYPE>::identifyFace(Node* nodes[3]) const
 {
 	for (unsigned i=0; i<3; i++)
 	{
diff --git a/MeshLib/Elements/Tet.h b/MeshLib/Elements/Tet.h
index 1aaa64e70975eebb48ca6648221ff1b146669ab8..e146722708c6247cc092848f140dfc5caa3fa1df 100644
--- a/MeshLib/Elements/Tet.h
+++ b/MeshLib/Elements/Tet.h
@@ -16,7 +16,7 @@
 
 namespace MeshLib {
 
-typedef TemplateTet<4> Tet;
+typedef TemplateTet<4,FEMElemType::TET4> Tet;
 
 }
 
diff --git a/MeshLib/Elements/Tri.h b/MeshLib/Elements/Tri.h
index 1cdcbf0ea27b8718513f3e6058ff1156028017dd..4ea28d0f60af3dd830b281bf0d0ade5ae8b4e6cb 100644
--- a/MeshLib/Elements/Tri.h
+++ b/MeshLib/Elements/Tri.h
@@ -16,7 +16,7 @@
 
 namespace MeshLib {
 
-typedef TemplateTri<3> Tri;
+typedef TemplateTri<3,FEMElemType::TRI3> Tri;
 
 }
 
diff --git a/MeshLib/MshEnums.h b/MeshLib/MshEnums.h
index 437fddffed41de42741e634e0c200282df3be5d6..efcef912af5ad60f8a3902818233a05c4999a2da 100644
--- a/MeshLib/MshEnums.h
+++ b/MeshLib/MshEnums.h
@@ -32,6 +32,32 @@ struct MshElemType
 	};
 };
 
+/**
+ * \brief Types of mesh elements supported by OpenGeoSys.
+ */
+struct FEMElemType
+{
+	enum type {
+		INVALID,
+		EDGE2,
+		EDGE3,
+		TRI3,
+		TRI6,
+		QUAD4,
+		QUAD8,
+		QUAD9,
+		TET4,
+		TET10,
+		HEX8,
+		HEX20,
+		HEX27,
+		PRISM6,
+		PRISM15,
+		PRISM18,
+		PYRAMID5
+	};
+};
+
 /**
  * \brief Describes a mesh quality criteria.
  */