diff --git a/MeshLib/Elements/CellRule.cpp b/MeshLib/Elements/CellRule.cpp
index 136bebd583dbf81cbe047decf937db3d80d90dc1..7381f66debada2f56b20a2f059fc1d87acfbf469 100644
--- a/MeshLib/Elements/CellRule.cpp
+++ b/MeshLib/Elements/CellRule.cpp
@@ -12,7 +12,6 @@
 
 #include "Element.h"
 #include "FaceRule.h"
-#include "MeshLib/Node.h"
 
 namespace MeshLib
 {
diff --git a/MeshLib/Elements/CellRule.h b/MeshLib/Elements/CellRule.h
index 99c6412ddec3620d0f54f94e18452f3910f6a4e2..cd79a0311f124787e0f07fafa33b79c468ec292a 100644
--- a/MeshLib/Elements/CellRule.h
+++ b/MeshLib/Elements/CellRule.h
@@ -10,6 +10,8 @@
 
 #pragma once
 
+#include "MeshLib/Node.h"
+
 namespace MeshLib
 {
 
@@ -30,6 +32,36 @@ public:
      * center of gravity to lie outside of the actual element
      */
     static bool testElementNodeOrder(Element const& e);
+
+protected:
+    /// Returns the ID of a face given an array of nodes.
+    template <typename ElementRule>
+    static unsigned identifyFace(Node const* const* _nodes,
+                                 Node const* nodes[3])
+    {
+        for (unsigned i = 0; i < ElementRule::n_faces; i++)
+        {
+            unsigned flag(0);
+            constexpr std::size_t n = sizeof(ElementRule::face_nodes[0]) /
+                                      sizeof(ElementRule::face_nodes[0][0]);
+            for (unsigned j = 0; j < n; j++)
+            {
+                for (unsigned k = 0; k < 3; k++)
+                {
+                    if (ElementRule::face_nodes[i][j] != 99 &&
+                        _nodes[ElementRule::face_nodes[i][j]] == nodes[k])
+                    {
+                        flag++;
+                    }
+                }
+            }
+            if (flag == 3)
+            {
+                return i;
+            }
+        }
+        return std::numeric_limits<unsigned>::max();
+    }
 }; /* class */
 
 }  // namespace MeshLib
diff --git a/MeshLib/Elements/HexRule.h b/MeshLib/Elements/HexRule.h
index 203c3b9c3a62c4a04e094ea858c6aed348a76ceb..30ba9c469cbea31105fe63dcf6e46e4cdd442005 100644
--- a/MeshLib/Elements/HexRule.h
+++ b/MeshLib/Elements/HexRule.h
@@ -49,32 +49,5 @@ public:
     /// Calculates the volume of a convex hexahedron by partitioning it into six
     /// tetrahedra.
     static double computeVolume(Node const* const* _nodes);
-
-protected:
-    template <std::size_t N>
-    static unsigned identifyFace(Node const* const* _nodes,
-                                 Node const* nodes[3],
-                                 unsigned const face_nodes[6][N])
-    {
-        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();
-    }
 };
 }  // namespace MeshLib
diff --git a/MeshLib/Elements/HexRule20.h b/MeshLib/Elements/HexRule20.h
index b5024f5aaa976050e7cd88d903b26e64ea6fb014..cfa67d6d6c53edbdddbcf38786506da139d3abf7 100644
--- a/MeshLib/Elements/HexRule20.h
+++ b/MeshLib/Elements/HexRule20.h
@@ -66,7 +66,7 @@ public:
     static unsigned identifyFace(Node const* const* _nodes,
                                  Node const* nodes[3])
     {
-        return HexRule::identifyFace(_nodes, nodes, face_nodes);
+        return CellRule::identifyFace<HexRule20>(_nodes, nodes);
     }
 }; /* class */
 
diff --git a/MeshLib/Elements/HexRule8.h b/MeshLib/Elements/HexRule8.h
index d1bb1ef2c79d9f43e1a4a4a540678c94f01a8d10..c4e5fdb9fbd46c1e1539df05c53317239ab0834d 100644
--- a/MeshLib/Elements/HexRule8.h
+++ b/MeshLib/Elements/HexRule8.h
@@ -65,7 +65,7 @@ public:
     static unsigned identifyFace(Node const* const* _nodes,
                                  Node const* nodes[3])
     {
-        return HexRule::identifyFace(_nodes, nodes, face_nodes);
+        return CellRule::identifyFace<HexRule8>(_nodes, nodes);
     }
 }; /* class */
 
diff --git a/MeshLib/Elements/PrismRule.h b/MeshLib/Elements/PrismRule.h
index d4ef88a753e907bb185e4f692c19006b4e3ea8af..9da0786033a7e1ac699e5e49ba409b6e770e981a 100644
--- a/MeshLib/Elements/PrismRule.h
+++ b/MeshLib/Elements/PrismRule.h
@@ -49,33 +49,5 @@ public:
     /// Calculates the volume of a convex hexahedron by partitioning it into six
     /// tetrahedra.
     static double computeVolume(Node const* const* _nodes);
-
-protected:
-    template <std::size_t N>
-    static unsigned identifyFace(Node const* const* _nodes,
-                                 Node const* nodes[3],
-                                 unsigned const face_nodes[5][N])
-    {
-        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();
-    }
 };
 }  // namespace MeshLib
diff --git a/MeshLib/Elements/PrismRule15.h b/MeshLib/Elements/PrismRule15.h
index de747689bfefbf1940e52bb55edd7759c5ec6b16..b277b84ee0828319b075c87162b65eece601da58 100644
--- a/MeshLib/Elements/PrismRule15.h
+++ b/MeshLib/Elements/PrismRule15.h
@@ -69,7 +69,7 @@ public:
     static unsigned identifyFace(Node const* const* _nodes,
                                  Node const* nodes[3])
     {
-        return PrismRule::identifyFace(_nodes, nodes, face_nodes);
+        return CellRule::identifyFace<PrismRule15>(_nodes, nodes);
     }
 }; /* class */
 
diff --git a/MeshLib/Elements/PrismRule6.h b/MeshLib/Elements/PrismRule6.h
index 37393b4eec99caa0d7f8a52e5eeaf1de5ce2b80a..c115ef586aed917859ff28c349256a70b1e75608 100644
--- a/MeshLib/Elements/PrismRule6.h
+++ b/MeshLib/Elements/PrismRule6.h
@@ -69,7 +69,7 @@ public:
     static unsigned identifyFace(Node const* const* _nodes,
                                  Node const* nodes[3])
     {
-        return PrismRule::identifyFace(_nodes, nodes, face_nodes);
+        return CellRule::identifyFace<PrismRule6>(_nodes, nodes);
     }
 }; /* class */
 
diff --git a/MeshLib/Elements/PyramidRule.h b/MeshLib/Elements/PyramidRule.h
index 6496bc513fe262277172b0a855059d72b199536e..58b2224277077cbab73e7bf7a2915b92477d846c 100644
--- a/MeshLib/Elements/PyramidRule.h
+++ b/MeshLib/Elements/PyramidRule.h
@@ -50,33 +50,5 @@ public:
     /// Calculates the volume of a convex hexahedron by partitioning it into six
     /// tetrahedra.
     static double computeVolume(Node const* const* _nodes);
-
-protected:
-    template <std::size_t N>
-    static unsigned identifyFace(Node const* const* _nodes,
-                                 Node const* nodes[3],
-                                 unsigned const face_nodes[5][N])
-    {
-        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();
-    }
 };
 }  // namespace MeshLib
diff --git a/MeshLib/Elements/PyramidRule13.h b/MeshLib/Elements/PyramidRule13.h
index 3a7028d6013dde1e1e52fe4bb0d35993abe88121..f8da93345786f342a83046d329648c625d9bd940 100644
--- a/MeshLib/Elements/PyramidRule13.h
+++ b/MeshLib/Elements/PyramidRule13.h
@@ -68,7 +68,7 @@ public:
     static unsigned identifyFace(Node const* const* _nodes,
                                  Node const* nodes[3])
     {
-        return PyramidRule::identifyFace(_nodes, nodes, face_nodes);
+        return CellRule::identifyFace<PyramidRule13>(_nodes, nodes);
     }
 }; /* class */
 
diff --git a/MeshLib/Elements/PyramidRule5.h b/MeshLib/Elements/PyramidRule5.h
index a80da0a5a2ed2a93e2a7b191c5d279515a569dbc..9cd1f8c433400743787f3e2a23f2da555effb4a0 100644
--- a/MeshLib/Elements/PyramidRule5.h
+++ b/MeshLib/Elements/PyramidRule5.h
@@ -68,7 +68,7 @@ public:
     static unsigned identifyFace(Node const* const* _nodes,
                                  Node const* nodes[3])
     {
-        return PyramidRule::identifyFace(_nodes, nodes, face_nodes);
+        return CellRule::identifyFace<PyramidRule5>(_nodes, nodes);
     }
 }; /* class */
 
diff --git a/MeshLib/Elements/TetRule.h b/MeshLib/Elements/TetRule.h
index 159b4a71a16283ac0d7b5f52a476bf0ab5a33463..117184b6562ce2309d30ee9590eb3f6ef1c4fbdf 100644
--- a/MeshLib/Elements/TetRule.h
+++ b/MeshLib/Elements/TetRule.h
@@ -49,32 +49,5 @@ public:
 
     /// Calculates the volume of the element
     static double computeVolume(Node const* const* _nodes);
-
-protected:
-    template <std::size_t N>
-    static unsigned identifyFace(Node const* const* _nodes,
-                                 Node const* nodes[3],
-                                 unsigned const face_nodes[4][N])
-    {
-        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();
-    }
 };
 }  // namespace MeshLib
diff --git a/MeshLib/Elements/TetRule10.h b/MeshLib/Elements/TetRule10.h
index 263540b16a0ae8f04c7724ba0301f6f5dcbcd6cd..9eeea81b5d307af8305c57aaca3266bafb5afaeb 100644
--- a/MeshLib/Elements/TetRule10.h
+++ b/MeshLib/Elements/TetRule10.h
@@ -63,7 +63,7 @@ public:
     static unsigned identifyFace(Node const* const* _nodes,
                                  Node const* nodes[3])
     {
-        return TetRule::identifyFace(_nodes, nodes, face_nodes);
+        return CellRule::identifyFace<TetRule10>(_nodes, nodes);
     }
 }; /* class */
 
diff --git a/MeshLib/Elements/TetRule4.h b/MeshLib/Elements/TetRule4.h
index de806f568bb5a29a2e7273ed03991522ba9eeb25..707d792c0ce47fc1d9ece774434026261100855b 100644
--- a/MeshLib/Elements/TetRule4.h
+++ b/MeshLib/Elements/TetRule4.h
@@ -63,7 +63,7 @@ public:
     static unsigned identifyFace(Node const* const* _nodes,
                                  Node const* nodes[3])
     {
-        return TetRule::identifyFace(_nodes, nodes, face_nodes);
+        return CellRule::identifyFace<TetRule4>(_nodes, nodes);
     }
 }; /* class */