diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index fec6cbbdcfef3d399004354642637a2f2572671d..16f9877c11d0ffac4c7d6c31099bca8759709266 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -355,4 +355,18 @@ bool dividedByPlane(const GeoLib::Point& a, const GeoLib::Point& b, const GeoLib
 	return false;
 }
 
+bool pointsOnAPlane(const GeoLib::Point& a, const GeoLib::Point& b, const GeoLib::Point& c, const GeoLib::Point& d)
+{
+	const GeoLib::Point AB(b[0]-a[0], b[1]-a[1], b[2]-a[2]);
+	const GeoLib::Point AC(c[0]-a[0], c[1]-a[1], c[2]-a[2]);
+	const GeoLib::Point AD(d[0]-a[0], d[1]-a[1], d[2]-a[2]);
+
+	double squared_scalar_triple = pow(GeoLib::scalarTriple(AC, AD, AB), 2);
+	double normalisation_factor  = (AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]) * 
+			                        (AC[0]*AC[0]+AC[1]*AC[1]+AC[2]*AC[2]) * 
+									(AD[0]*AD[0]+AD[1]*AD[1]+AD[2]*AD[2]);
+
+	return (squared_scalar_triple/normalisation_factor < std::numeric_limits<double>::epsilon());
+}
+
 } // end namespace GeoLib
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index 9623083b768a644a16d7ee871ba84a0336eea0e4..7ff3f05b271c0930be5760fd93321d926b8b4471 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -150,6 +150,10 @@ double scalarTriple(GeoLib::Point const& u, GeoLib::Point const& v, GeoLib::Poin
  bool dividedByPlane(const GeoLib::Point& a, const GeoLib::Point& b, 
 	 const GeoLib::Point& c, const GeoLib::Point& d);
 
+ /// Checks if the four given points are located on a plane.
+ bool pointsOnAPlane(const GeoLib::Point& a, const GeoLib::Point& b, 
+	 const GeoLib::Point& c, const GeoLib::Point& d);
+
 
 } // end namespace GeoLib
 
diff --git a/MeshLib/Elements/TemplateQuad.tpp b/MeshLib/Elements/TemplateQuad.tpp
index 8f8fe18782cad454db41105b811b163ab246cda8..70cd0aa073bd7da650a8e54e89e0beae21406f40 100644
--- a/MeshLib/Elements/TemplateQuad.tpp
+++ b/MeshLib/Elements/TemplateQuad.tpp
@@ -125,17 +125,7 @@ bool TemplateQuad<NNODES,CELLQUADTYPE>::isValid() const
 	if (this->_area <= std::numeric_limits<double>::epsilon())
 		return false;
 
-	const GeoLib::Point AB((*_nodes[1])[0]-(*_nodes[0])[0], (*_nodes[1])[1]-(*_nodes[0])[1], (*_nodes[1])[2]-(*_nodes[0])[2]);
-	const GeoLib::Point AC((*_nodes[2])[0]-(*_nodes[0])[0], (*_nodes[2])[1]-(*_nodes[0])[1], (*_nodes[2])[2]-(*_nodes[0])[2]);
-	const GeoLib::Point AD((*_nodes[3])[0]-(*_nodes[0])[0], (*_nodes[3])[1]-(*_nodes[0])[1], (*_nodes[3])[2]-(*_nodes[0])[2]);
-
-	// check if all points lie on the same plane
-	double squared_scalar_triple = pow(GeoLib::scalarTriple(AC, AD, AB), 2);
-	double normalisation_factor  = (AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]) * 
-			                        (AC[0]*AC[0]+AC[1]*AC[1]+AC[2]*AC[2]) * 
-									(AD[0]*AD[0]+AD[1]*AD[1]+AD[2]*AD[2]);
-
-	if (squared_scalar_triple/normalisation_factor < std::numeric_limits<double>::epsilon())
+	if (GeoLib::pointsOnAPlane(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3]))
 	{
 		// check if quad is convex
 		if (GeoLib::dividedByPlane(*_nodes[0], *_nodes[2], *_nodes[1], *_nodes[3]) &&