diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index b8006bd6c43f744aa8af341580431739bbe5b913..51e46fefd1777995f3d435a90c5a3cf76ce8c094 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -321,15 +321,6 @@ bool barycentricPointInTriangle(MathLib::Point3d const& p,
     return true;
 }
 
-double calcTriangleArea(MathLib::Point3d const& a,
-    MathLib::Point3d const& b, MathLib::Point3d const& c)
-{
-    MathLib::Vector3 const u(a,c);
-    MathLib::Vector3 const v(a,b);
-    MathLib::Vector3 const w(MathLib::crossProduct(u, v));
-    return 0.5 * w.getLength();
-}
-
 void computeRotationMatrixToXZ(MathLib::Vector3 const& plane_normal, MathLib::DenseMatrix<double> & rot_mat)
 {
     // *** some frequently used terms ***
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index cf5f1d83c5a974578792c1b3302a16b85853275a..a02c7cdd8e9f9f319a6a0c15d4075ce37931ec9e 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -176,14 +176,6 @@ MathLib::DenseMatrix<double> rotatePointsToXY(InputIterator1 p_pnts_begin,
  */
 void rotatePointsToXZ(std::vector<GeoLib::Point*> &pnts);
 
-/**
- * Calculates the area of the triangle defined by its edge nodes a, b and c.
- * The formula is \f$A= \frac{1}{2} \cdot |u \times v|\f$, i.e. half of the area of the
- * parallelogram specified by the vectors\f$u=b-a\f$ and \f$v=c-a\f$.
- */
-double calcTriangleArea(MathLib::Point3d const& a, MathLib::Point3d const& b,
-    MathLib::Point3d const& c);
-
 /**
  * Tests if the given point p is within the triangle, defined by its edge nodes a, b and c.
  * Using two eps-values it is possible to test an 'epsilon' neighbourhood around the triangle
diff --git a/GeoLib/IO/TINInterface.cpp b/GeoLib/IO/TINInterface.cpp
index 566e2ea6ec279fd9106fdfbfa4ee4a2251440025..ab9a892d431e8489fbf5b0811e89dbeb0db44f93 100644
--- a/GeoLib/IO/TINInterface.cpp
+++ b/GeoLib/IO/TINInterface.cpp
@@ -20,6 +20,8 @@
 #include "GeoLib/Surface.h"
 #include "GeoLib/Triangle.h"
 
+#include "MathLib/GeometricBasics.h"
+
 namespace GeoLib
 {
 namespace IO
@@ -92,7 +94,7 @@ GeoLib::Surface* TINInterface::readTIN(std::string const& fname,
 
         // check area of triangle
         double const d_eps(std::numeric_limits<double>::epsilon());
-        if (GeoLib::calcTriangleArea(p0, p1, p2) < d_eps) {
+        if (MathLib::calcTriangleArea(p0, p1, p2) < d_eps) {
             ERR("readTIN: Triangle %d has zero area.", id);
             if (errors)
                 errors->push_back (std::string("readTIN: Triangle ")
diff --git a/MathLib/GeometricBasics.cpp b/MathLib/GeometricBasics.cpp
index caab906d2d31b68da916f1d376d3763a2fd05d53..3d26f8e877f08863d8d200d1562be7bf475514d6 100644
--- a/MathLib/GeometricBasics.cpp
+++ b/MathLib/GeometricBasics.cpp
@@ -36,6 +36,15 @@ double calcTetrahedronVolume(MathLib::Point3d const& a,
     return std::abs(MathLib::scalarTriple(ac, ad, ab)) / 6.0;
 }
 
+double calcTriangleArea(MathLib::Point3d const& a, MathLib::Point3d const& b,
+                        MathLib::Point3d const& c)
+{
+    MathLib::Vector3 const u(a, c);
+    MathLib::Vector3 const v(a, b);
+    MathLib::Vector3 const w(MathLib::crossProduct(u, v));
+    return 0.5 * w.getLength();
+}
+
 bool isPointInTetrahedron(MathLib::Point3d const& p, MathLib::Point3d const& a,
                           MathLib::Point3d const& b, MathLib::Point3d const& c,
                           MathLib::Point3d const& d, double eps)
diff --git a/MathLib/GeometricBasics.h b/MathLib/GeometricBasics.h
index e97a9b5e407b2ba461c298cbb9c9084994676bb5..ba569dddbf31abce433629e5c7399c61d2db1828 100644
--- a/MathLib/GeometricBasics.h
+++ b/MathLib/GeometricBasics.h
@@ -45,6 +45,15 @@ double calcTetrahedronVolume(MathLib::Point3d const& x1,
                              MathLib::Point3d const& x2,
                              MathLib::Point3d const& x3,
                              MathLib::Point3d const& x4);
+
+/**
+ * Calculates the area of the triangle defined by its edge nodes a, b and c.
+ * The formula is \f$A= \frac{1}{2} \cdot |u \times v|\f$, i.e. half of the area of the
+ * parallelogram specified by the vectors\f$u=b-a\f$ and \f$v=c-a\f$.
+ */
+double calcTriangleArea(MathLib::Point3d const& a, MathLib::Point3d const& b,
+                        MathLib::Point3d const& c);
+
 /**
  * Tests if the given point p is located within a tetrahedron spanned by points
  * a, b, c, d.
diff --git a/MeshLib/Elements/QuadRule4.cpp b/MeshLib/Elements/QuadRule4.cpp
index b1ec6e43cb65342b5c70cc99690af15b17955a07..fec6900320583cb0297e6e37b226278b0d8e1450 100644
--- a/MeshLib/Elements/QuadRule4.cpp
+++ b/MeshLib/Elements/QuadRule4.cpp
@@ -12,6 +12,7 @@
 #include "logog/include/logog.hpp"
 
 #include "GeoLib/AnalyticalGeometry.h"
+#include "MathLib/GeometricBasics.h"
 
 #include "MeshLib/Node.h"
 
@@ -27,8 +28,8 @@ const unsigned QuadRule4::edge_nodes[4][2] =
 
 double QuadRule4::computeVolume(Node const* const* _nodes)
 {
-    return GeoLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2])
-         + GeoLib::calcTriangleArea(*_nodes[2], *_nodes[3], *_nodes[0]);
+    return MathLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2])
+         + MathLib::calcTriangleArea(*_nodes[2], *_nodes[3], *_nodes[0]);
 }
 
 bool QuadRule4::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps)
diff --git a/MeshLib/Elements/TriRule3.cpp b/MeshLib/Elements/TriRule3.cpp
index 2da9a0877c058acf73aaafd62de099de4b80c21c..0927588134cac1a6860d07766bff4918e184bfab 100644
--- a/MeshLib/Elements/TriRule3.cpp
+++ b/MeshLib/Elements/TriRule3.cpp
@@ -12,6 +12,7 @@
 #include "logog/include/logog.hpp"
 
 #include "GeoLib/AnalyticalGeometry.h"
+#include "MathLib/GeometricBasics.h"
 
 #include "MeshLib/Node.h"
 
@@ -26,7 +27,7 @@ const unsigned TriRule3::edge_nodes[3][2] =
 
 double TriRule3::computeVolume(Node const* const* _nodes)
 {
-    return GeoLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2]);
+    return MathLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2]);
 }
 
 bool TriRule3::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps)
diff --git a/NumLib/Function/LinearInterpolationOnSurface.cpp b/NumLib/Function/LinearInterpolationOnSurface.cpp
index 27fd3cc3455699e4e312a700127babe410f697e8..a9a0da32ffcbad0649e9c49ed98043a50bcb7893 100644
--- a/NumLib/Function/LinearInterpolationOnSurface.cpp
+++ b/NumLib/Function/LinearInterpolationOnSurface.cpp
@@ -17,10 +17,11 @@
 #include <cassert>
 #include <numeric>
 
-#include "MathLib/Vector3.h"
 #include "GeoLib/Surface.h"
 #include "GeoLib/Triangle.h"
 #include "GeoLib/AnalyticalGeometry.h"
+#include "MathLib/Vector3.h"
+#include "MathLib/GeometricBasics.h"
 #include "MeshLib/Mesh.h"
 #include "MeshLib/Node.h"
 #include "MeshGeoToolsLib/MeshNodesAlongSurface.h"
@@ -77,7 +78,7 @@ double LinearInterpolationOnSurface::interpolateInTri(
     GeoLib::Point const& v2(pnts[1]);
     GeoLib::Point const& v3(pnts[2]);
     GeoLib::Point const& v_pnt(pnts[3]);
-    const double area = GeoLib::calcTriangleArea(v1, v2, v3);
+    const double area = MathLib::calcTriangleArea(v1, v2, v3);
 
     if (area==.0) {
         // take average if all points have the same coordinates