diff --git a/Applications/DataExplorer/DataView/DirectConditionGenerator.cpp b/Applications/DataExplorer/DataView/DirectConditionGenerator.cpp
index 1607922052d4c29860da79d867a7ebcace70e1e7..aea1f39dfa4e68616f7844884d6be6cfde9d1758 100644
--- a/Applications/DataExplorer/DataView/DirectConditionGenerator.cpp
+++ b/Applications/DataExplorer/DataView/DirectConditionGenerator.cpp
@@ -83,7 +83,7 @@ const std::vector< std::pair<std::size_t,double> >& DirectConditionGenerator::di
         return _direct_values;
     }
 
-    MathLib::Vector3 const dir(0.0, 0.0, -1.0);
+    Eigen::Vector3d const dir({0.0, 0.0, -1.0});
     double const angle(90);
     std::string const prop_name("bulk_node_ids");
     std::unique_ptr<MeshLib::Mesh> surface_mesh(
diff --git a/Applications/DataExplorer/DataView/MeshView.cpp b/Applications/DataExplorer/DataView/MeshView.cpp
index 20a1da06844a56e21931248854b666eb2806d487..11d10e2d0659f70b1d436d59a93ac69c854cf165 100644
--- a/Applications/DataExplorer/DataView/MeshView.cpp
+++ b/Applications/DataExplorer/DataView/MeshView.cpp
@@ -319,7 +319,7 @@ void MeshView::extractSurfaceMesh()
         return;
     }
 
-    MathLib::Vector3 const& dir (dlg.getNormal());
+    Eigen::Vector3d const& dir(dlg.getNormal());
     int const tolerance (dlg.getTolerance());
     std::unique_ptr<MeshLib::Mesh> sfc_mesh(
         MeshLib::MeshSurfaceExtraction::getMeshSurface(
diff --git a/Applications/DataExplorer/DataView/SurfaceExtractionDialog.cpp b/Applications/DataExplorer/DataView/SurfaceExtractionDialog.cpp
index 809840ab095aa45031e990d21c35c771fb7de3a3..5df9ef941582ade0bd19721f97ded7849668103c 100644
--- a/Applications/DataExplorer/DataView/SurfaceExtractionDialog.cpp
+++ b/Applications/DataExplorer/DataView/SurfaceExtractionDialog.cpp
@@ -27,9 +27,9 @@ SurfaceExtractionDialog::SurfaceExtractionDialog(QDialog* parent)
 
 void SurfaceExtractionDialog::accept()
 {
-    _dir = MathLib::Vector3(xNormalEdit->text().toDouble(),
+    _dir = Eigen::Vector3d({xNormalEdit->text().toDouble(),
                             yNormalEdit->text().toDouble(),
-                            zNormalEdit->text().toDouble());
+                            zNormalEdit->text().toDouble()});
     _tolerance = degreesSpinBox->text().toInt();
 
     this->done(QDialog::Accepted);
diff --git a/Applications/DataExplorer/DataView/SurfaceExtractionDialog.h b/Applications/DataExplorer/DataView/SurfaceExtractionDialog.h
index e53278525c1355de74c9c564b79140a27e4c86ea..eef26f2a6e1f007d0542c65bd20aa3dc1f1f5737 100644
--- a/Applications/DataExplorer/DataView/SurfaceExtractionDialog.h
+++ b/Applications/DataExplorer/DataView/SurfaceExtractionDialog.h
@@ -35,7 +35,7 @@ public:
     ~SurfaceExtractionDialog() override = default;
 
     int getTolerance() const { return _tolerance; }
-    MathLib::Vector3 const& getNormal() const { return _dir; }
+    Eigen::Vector3d const& getNormal() const { return _dir; }
 
 private slots:
     /// Instructions if the OK-Button has been pressed.
@@ -46,5 +46,5 @@ private slots:
 
 private:
     int _tolerance{90};
-    MathLib::Vector3 _dir{0, 0, -1};
+    Eigen::Vector3d _dir{0, 0, -1};
 };
diff --git a/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp b/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp
index 14c0c5b8d51b1f280b9469ccf48b1faee7747cd1..4cf726f5b33563cb3b4323903b2a835b3855c0cf 100644
--- a/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp
+++ b/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp
@@ -179,7 +179,7 @@ int main (int argc, char* argv[])
         MeshLib::IO::readMeshFromFile(mesh_arg.getValue()));
     INFO("done.");
     INFO("Extracting top surface of mesh '{:s}' ... ", mesh_arg.getValue());
-    const MathLib::Vector3 dir(0,0,-1);
+    Eigen::Vector3d const dir({0, 0, -1});
     double const angle(90);
     std::unique_ptr<MeshLib::Mesh> surface_mesh(
         MeshLib::MeshSurfaceExtraction::getMeshSurface(*subsurface_mesh, dir,
diff --git a/Applications/Utils/MeshEdit/ExtractSurface.cpp b/Applications/Utils/MeshEdit/ExtractSurface.cpp
index 79dfab15521c0f4de0f69af19c2304ab9a96b000..ce30f87a1ecf75a192c2f67ae0530b7d06a7626f 100644
--- a/Applications/Utils/MeshEdit/ExtractSurface.cpp
+++ b/Applications/Utils/MeshEdit/ExtractSurface.cpp
@@ -111,7 +111,7 @@ int main (int argc, char* argv[])
          mesh->getNumberOfElements());
 
     // extract surface
-    MathLib::Vector3 const dir(x.getValue(), y.getValue(), z.getValue());
+    Eigen::Vector3d const dir({x.getValue(), y.getValue(), z.getValue()});
     double const angle(angle_arg.getValue());
     std::unique_ptr<MeshLib::Mesh> surface_mesh(
         MeshLib::MeshSurfaceExtraction::getMeshSurface(
diff --git a/Applications/Utils/MeshGeoTools/computeSurfaceNodeIDsInPolygonalRegion.cpp b/Applications/Utils/MeshGeoTools/computeSurfaceNodeIDsInPolygonalRegion.cpp
index 671d896427c817e8dc0164dd4664d8018dee3b19..b3a5a3d449a5ef5c7ff7128ac87fc24bd566d84a 100644
--- a/Applications/Utils/MeshGeoTools/computeSurfaceNodeIDsInPolygonalRegion.cpp
+++ b/Applications/Utils/MeshGeoTools/computeSurfaceNodeIDsInPolygonalRegion.cpp
@@ -102,11 +102,11 @@ int main (int argc, char* argv[])
          geo_objs.getPointVec(geo_name)->size(),
          geo_objs.getPolylineVec(geo_name)->size());
 
-    MathLib::Vector3 const dir(0.0, 0.0, -1.0);
+    Eigen::Vector3d const dir({0.0, 0.0, -1.0});
     double angle(90);
 
     auto computeElementTopSurfaceAreas = [](MeshLib::Mesh const& mesh,
-        MathLib::Vector3 const& d, double angle)
+        Eigen::Vector3d const& d, double angle)
     {
         std::unique_ptr<MeshLib::Mesh> surface_mesh(
             MeshLib::MeshSurfaceExtraction::getMeshSurface(mesh, d, angle));
@@ -115,8 +115,9 @@ int main (int argc, char* argv[])
     };
 
     std::vector<double> areas(computeElementTopSurfaceAreas(*mesh, dir, angle));
+    MathLib::Vector3 const mal_dir(dir[0], dir[1], dir[2]);
     std::vector<MeshLib::Node*> all_sfc_nodes(
-        MeshLib::MeshSurfaceExtraction::getSurfaceNodes(*mesh, dir, angle)
+        MeshLib::MeshSurfaceExtraction::getSurfaceNodes(*mesh, mal_dir, angle)
     );
 
     std::for_each(all_sfc_nodes.begin(), all_sfc_nodes.end(),
diff --git a/MeshGeoToolsLib/GeoMapper.cpp b/MeshGeoToolsLib/GeoMapper.cpp
index 41890c56780b293aeb16f33d9d080ec9b8c93952..97dcd824304a631ac2d14a93c5e264f799c49b63 100644
--- a/MeshGeoToolsLib/GeoMapper.cpp
+++ b/MeshGeoToolsLib/GeoMapper.cpp
@@ -82,8 +82,9 @@ void GeoMapper::mapOnMesh(MeshLib::Mesh const*const mesh)
     }
     else
     {
-        const MathLib::Vector3 dir(0,0,-1);
-        _surface_mesh = MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh, dir, 90);
+        Eigen::Vector3d const dir(0,0,-1);
+        _surface_mesh =
+            MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh, dir, 90);
     }
 
     // init grid
@@ -577,7 +578,7 @@ void GeoMapper::advancedMapOnMesh(MeshLib::Mesh const& mesh)
     if (mesh.getDimension()<3) {
         _surface_mesh = new MeshLib::Mesh(mesh);
     } else {
-        const MathLib::Vector3 dir(0,0,-1);
+        Eigen::Vector3d const dir({0, 0, -1});
         _surface_mesh =
             MeshLib::MeshSurfaceExtraction::getMeshSurface(mesh, dir, 90+1e-6);
     }
diff --git a/MeshLib/MeshEditing/AddLayerToMesh.cpp b/MeshLib/MeshEditing/AddLayerToMesh.cpp
index 8a4b25cd1b4a70fdd34d7b0b131d37031ed63aa2..4b7d855a9d6f5ec36eeaede2e03a41315df3c440 100644
--- a/MeshLib/MeshEditing/AddLayerToMesh.cpp
+++ b/MeshLib/MeshEditing/AddLayerToMesh.cpp
@@ -88,8 +88,8 @@ MeshLib::Mesh* addLayerToMesh(MeshLib::Mesh const& mesh, double thickness,
                               bool copy_material_ids)
 {
     INFO("Extracting top surface of mesh '{:s}' ... ", mesh.getName());
-    int const flag = (on_top) ? -1 : 1;
-    const MathLib::Vector3 dir(0, 0, flag);
+    double const flag = on_top ? -1.0 : 1.0;
+    Eigen::Vector3d const dir({0, 0, flag});
     double const angle(90);
     std::unique_ptr<MeshLib::Mesh> sfc_mesh (nullptr);
 
diff --git a/MeshLib/MeshSurfaceExtraction.cpp b/MeshLib/MeshSurfaceExtraction.cpp
index 76f792587ac5a29d6684693349fb8cab99a78901..2c39902520d296f49ddd62a9973fb37354785d9c 100644
--- a/MeshLib/MeshSurfaceExtraction.cpp
+++ b/MeshLib/MeshSurfaceExtraction.cpp
@@ -227,7 +227,7 @@ std::vector<double> MeshSurfaceExtraction::getSurfaceAreaForNodes(
 }
 
 MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(
-    const MeshLib::Mesh& subsfc_mesh, const MathLib::Vector3& dir, double angle,
+    const MeshLib::Mesh& subsfc_mesh, Eigen::Vector3d const& dir, double angle,
     std::string const& subsfc_node_id_prop_name,
     std::string const& subsfc_element_id_prop_name,
     std::string const& face_id_prop_name)
@@ -287,20 +287,18 @@ void MeshSurfaceExtraction::get2DSurfaceElements(
     std::vector<MeshLib::Element*>& sfc_elements,
     std::vector<std::size_t>& element_to_bulk_element_id_map,
     std::vector<std::size_t>& element_to_bulk_face_id_map,
-    const MathLib::Vector3& dir, double angle, unsigned mesh_dimension)
+    Eigen::Vector3d const& dir, double angle, unsigned mesh_dimension)
 {
-    auto const d = Eigen::Map<Eigen::Vector3d const>(dir.getCoords());
-
     if (mesh_dimension < 2 || mesh_dimension > 3)
     {
         ERR("Cannot handle meshes of dimension {:i}", mesh_dimension);
     }
 
-    bool const complete_surface = (MathLib::scalarProduct(dir, dir) == 0);
+    bool const complete_surface = (dir.dot(dir) == 0);
 
     double const pi(boost::math::constants::pi<double>());
     double const cos_theta(std::cos(angle * pi / 180.0));
-    Eigen::Vector3d const norm_dir(d.normalized());
+    Eigen::Vector3d const norm_dir(dir.normalized());
 
     for (auto const* elem : all_elements)
     {
@@ -374,9 +372,11 @@ std::vector<MeshLib::Node*> MeshSurfaceExtraction::getSurfaceNodes(
     std::vector<MeshLib::Element*> sfc_elements;
     std::vector<std::size_t> element_to_bulk_element_id_map;
     std::vector<std::size_t> element_to_bulk_face_id_map;
+
+    auto const edir = Eigen::Map<Eigen::Vector3d const>(dir.getCoords());
     get2DSurfaceElements(
         mesh.getElements(), sfc_elements, element_to_bulk_element_id_map,
-        element_to_bulk_face_id_map, dir, angle, mesh.getDimension());
+        element_to_bulk_face_id_map, edir, angle, mesh.getDimension());
 
     std::vector<MeshLib::Node*> surface_nodes;
     std::tie(surface_nodes, std::ignore) =
diff --git a/MeshLib/MeshSurfaceExtraction.h b/MeshLib/MeshSurfaceExtraction.h
index ec01996ef0fb95b7f574b497ec271a207460ae8a..5735d30a4199615ffe8779324039b5b71cd67b69 100644
--- a/MeshLib/MeshSurfaceExtraction.h
+++ b/MeshLib/MeshSurfaceExtraction.h
@@ -64,7 +64,7 @@ public:
      */
     static MeshLib::Mesh* getMeshSurface(
         const MeshLib::Mesh& subsfc_mesh,
-        const MathLib::Vector3& dir,
+        Eigen::Vector3d const& dir,
         double angle,
         std::string const& subsfc_node_id_prop_name = "",
         std::string const& subsfc_element_id_prop_name = "",
@@ -77,7 +77,7 @@ private:
         std::vector<MeshLib::Element*>& sfc_elements,
         std::vector<std::size_t>& element_to_bulk_element_id_map,
         std::vector<std::size_t>& element_to_bulk_face_id_map,
-        const MathLib::Vector3& dir,
+        Eigen::Vector3d const& dir,
         double angle,
         unsigned mesh_dimension);
 };
diff --git a/Tests/MeshLib/TestAddLayerToMesh.cpp b/Tests/MeshLib/TestAddLayerToMesh.cpp
index 0f9c2d25a2f7691ca5397a53df21cf22919c4fa1..e23cea96a10a18d757d77ba4f4b3143c3cf97ced 100644
--- a/Tests/MeshLib/TestAddLayerToMesh.cpp
+++ b/Tests/MeshLib/TestAddLayerToMesh.cpp
@@ -219,7 +219,7 @@ TEST(MeshLib, AddTopLayerToHexMesh)
         MeshLib::MeshInformation::getNumberOfElementTypes(*result);
     ASSERT_EQ(150, n_elems.at(MeshLib::MeshElemType::HEXAHEDRON));
 
-    MathLib::Vector3 const dir(0, 0, -1);
+    Eigen::Vector3d const dir({0, 0, -1});
     std::unique_ptr<MeshLib::Mesh> const test_input (
         MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh, dir, 90));
     std::unique_ptr<MeshLib::Mesh> const test_output (
@@ -243,7 +243,7 @@ TEST(MeshLib, AddBottomLayerToHexMesh)
         MeshLib::MeshInformation::getNumberOfElementTypes(*result);
     ASSERT_EQ(150, n_elems.at(MeshLib::MeshElemType::HEXAHEDRON));
 
-    MathLib::Vector3 const dir(0, 0, 1);
+    Eigen::Vector3d const dir({0, 0, 1});
     std::unique_ptr<MeshLib::Mesh> const test_input (
         MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh, dir, 90));
     std::unique_ptr<MeshLib::Mesh> const test_output (
@@ -270,7 +270,7 @@ TEST(MeshLib, AddTopLayerToPrismMesh)
     ASSERT_EQ(50, n_elems.at(MeshLib::MeshElemType::TRIANGLE));
     ASSERT_EQ(100, n_elems.at(MeshLib::MeshElemType::PRISM));
 
-    MathLib::Vector3 const dir(0, 0, -1);
+    Eigen::Vector3d const dir({0, 0, -1});
     std::unique_ptr<MeshLib::Mesh> test_input (
         MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh2, dir, 90));
     std::unique_ptr<MeshLib::Mesh> test_output (
@@ -312,7 +312,7 @@ TEST(MeshLib, AddBottomLayerToPrismMesh)
     ASSERT_EQ(mesh2->getNumberOfElements(), std::count(new_mats->cbegin(), new_mats->cend(), 0));
     ASSERT_EQ(mesh->getNumberOfElements(), std::count(new_mats->cbegin(), new_mats->cend(), 1));
 
-    MathLib::Vector3 const dir(0, 0, 1);
+    Eigen::Vector3d const dir({0, 0, 1});
     std::unique_ptr<MeshLib::Mesh> test_input (
         MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh2, dir, 90));
     std::unique_ptr<MeshLib::Mesh> test_output (