diff --git a/Applications/DataExplorer/VtkVis/VtkPolylinesSource.cpp b/Applications/DataExplorer/VtkVis/VtkPolylinesSource.cpp
index 0719fd8c3a50536221a71a95b453fddbf6770c88..b27cd80316a696acdef8dc7b27b6d8db8ea3ed05 100644
--- a/Applications/DataExplorer/VtkVis/VtkPolylinesSource.cpp
+++ b/Applications/DataExplorer/VtkVis/VtkPolylinesSource.cpp
@@ -118,7 +118,7 @@ int VtkPolylinesSource::RequestData(vtkInformation* request,
         for (int i = 0; i < numVerts; i++)
         {
             const GeoLib::Point* point = (*_polylines)[j]->getPoint(i);
-            const double* coords = point->getCoords();
+            const double* coords = point->data();
             newPoints->InsertNextPoint(coords);
         }
 
diff --git a/Applications/DataExplorer/VtkVis/VtkStationSource.cpp b/Applications/DataExplorer/VtkVis/VtkStationSource.cpp
index c6c82543a69bc9b70e4cf0679ddf8e45506e4dc3..445a6a9aea19dfe4f3569105ec3e38aae914a0b7 100644
--- a/Applications/DataExplorer/VtkVis/VtkStationSource.cpp
+++ b/Applications/DataExplorer/VtkVis/VtkStationSource.cpp
@@ -141,8 +141,7 @@ int VtkStationSource::RequestData(vtkInformation* request,
     // Generate graphic objects
     for (auto station : *_stations)
     {
-        double coords[3] = {(*station)[0], (*station)[1], (*station)[2]};
-        vtkIdType sid = newStations->InsertNextPoint(coords);
+        vtkIdType sid = newStations->InsertNextPoint(station->data());
         station_ids->InsertNextValue(site_count);
         if (useStationValues)
         {
@@ -164,9 +163,7 @@ int VtkStationSource::RequestData(vtkInformation* request,
 
             for (std::size_t i = 1; i < nLayers; i++)
             {
-                auto* pCoords = const_cast<double*>(profile[i]->getCoords());
-                double loc[3] = {pCoords[0], pCoords[1], pCoords[2]};
-                newStations->InsertNextPoint(loc);
+                newStations->InsertNextPoint(profile[i]->data());
                 station_ids->InsertNextValue(site_count);
                 newLines->InsertNextCell(2);
                 newLines->InsertCellPoint(
diff --git a/Applications/DataExplorer/VtkVis/VtkSurfacesSource.cpp b/Applications/DataExplorer/VtkVis/VtkSurfacesSource.cpp
index 0d85259cf5f469b908e2bf569acb8dc8c9ffcd85..2b6ebedac5d7a8b31d37e6b4e99dc794a687b857 100644
--- a/Applications/DataExplorer/VtkVis/VtkSurfacesSource.cpp
+++ b/Applications/DataExplorer/VtkVis/VtkSurfacesSource.cpp
@@ -65,15 +65,14 @@ int VtkSurfacesSource::RequestData(vtkInformation* request,
     (void)request;
     (void)inputVector;
 
-    const int nSurfaces = _surfaces->size();
-    if (nSurfaces == 0)
+    if (_surfaces->empty())
     {
         return 0;
     }
 
     const std::vector<GeoLib::Point*>* surfacePoints =
         (*_surfaces)[0]->getPointVec();
-    std::size_t nPoints = surfacePoints->size();
+    std::size_t const nPoints = surfacePoints->size();
 
     vtkSmartPointer<vtkInformation> outInfo =
         outputVector->GetInformationObject(0);
@@ -98,9 +97,7 @@ int VtkSurfacesSource::RequestData(vtkInformation* request,
 
     for (std::size_t i = 0; i < nPoints; ++i)
     {
-        const double* coords =
-            const_cast<double*>((*surfacePoints)[i]->getCoords());
-        newPoints->SetPoint(i, coords);
+        newPoints->SetPoint(i, (*surfacePoints)[i]->data());
     }
 
     int count(0);
diff --git a/Applications/FileIO/GocadIO/GocadAsciiReader.cpp b/Applications/FileIO/GocadIO/GocadAsciiReader.cpp
index 49e5c1b7aa7f1cdb704680c23540775ad4059f82..f2c0cc1f0362c4c51355b20c8aed201b730f6c2a 100644
--- a/Applications/FileIO/GocadIO/GocadAsciiReader.cpp
+++ b/Applications/FileIO/GocadIO/GocadAsciiReader.cpp
@@ -344,8 +344,7 @@ bool parseNodes(std::ifstream& in,
             std::size_t ref_id;
             std::string keyword;
             sstr >> keyword >> new_id >> ref_id;
-            nodes.push_back(
-                new MeshLib::Node(nodes[ref_id]->getCoords(), new_id));
+            nodes.push_back(new MeshLib::Node(nodes[ref_id]->data(), new_id));
         }
         node_id_map[nodes.back()->getID()] = nodes.size() - 1;
         pos = in.tellg();
diff --git a/Applications/FileIO/GocadIO/GocadSGridReader.cpp b/Applications/FileIO/GocadIO/GocadSGridReader.cpp
index 3d7d96e0e922de391b3509c129c3cb38eb67f7df..7442ed3b60b2a286a48445c4a5d5a6e4066a21a7 100644
--- a/Applications/FileIO/GocadIO/GocadSGridReader.cpp
+++ b/Applications/FileIO/GocadIO/GocadSGridReader.cpp
@@ -639,8 +639,7 @@ void GocadSGridReader::applySplitInformation(
     for (auto split_node : _split_nodes)
     {
         std::size_t const new_node_pos(nodes.size());
-        nodes.push_back(
-            new MeshLib::Node(split_node->getCoords(), new_node_pos));
+        nodes.push_back(new MeshLib::Node(split_node->data(), new_node_pos));
 
         // get grid coordinates
         std::array<std::size_t, 3> const& gc(split_node->getGridCoords());
diff --git a/Applications/FileIO/TetGenInterface.cpp b/Applications/FileIO/TetGenInterface.cpp
index 0d7f0d54ebeb12b5bba46860e645a82582fb0abd..aef1f947b3a197f7265cd48f346b1c509a4324c9 100644
--- a/Applications/FileIO/TetGenInterface.cpp
+++ b/Applications/FileIO/TetGenInterface.cpp
@@ -826,7 +826,7 @@ void TetGenInterface::write3dElements(
         if (materialIds)
         {
             attribute_points.emplace_back(
-                MeshLib::getCenterOfGravity(*elements[i]).getCoords(),
+                MeshLib::getCenterOfGravity(*elements[i]).data(),
                 (*materialIds)[i]);
         }
     }
diff --git a/Applications/Utils/MeshEdit/AddFaultToVoxelGrid.cpp b/Applications/Utils/MeshEdit/AddFaultToVoxelGrid.cpp
index 1472a200dbca1eb0e4edf4abd1af8e3ea1ed1cf2..ed87ee390e8b5c7e09a454f8c55f165553e65c9e 100644
--- a/Applications/Utils/MeshEdit/AddFaultToVoxelGrid.cpp
+++ b/Applications/Utils/MeshEdit/AddFaultToVoxelGrid.cpp
@@ -48,9 +48,8 @@ bool testTriangleIntersectingAABB(MeshLib::Node const& n0,
 {
     // Translate triangle as conceptually moving AABB to origin
     Eigen::Matrix3d v;
-    v << Eigen::Vector3d(n0.getCoords()) - c,
-        Eigen::Vector3d(n1.getCoords()) - c,
-        Eigen::Vector3d(n2.getCoords()) - c;
+    v << n0.asEigenVector3d() - c, n1.asEigenVector3d() - c,
+        n2.asEigenVector3d() - c;
 
     // Test the three axes corresponding to the face normals of AABB b
     if (((v.rowwise().minCoeff() - e).array() > 0).any() ||
@@ -127,7 +126,7 @@ void markFaults(MeshLib::Mesh& mesh, MeshLib::Mesh const& fault,
         }
 
         // test if voxel is intersecting triangle
-        Eigen::Vector3d const c(centre_pnt.getCoords());
+        auto const& c(centre_pnt.asEigenVector3d());
         for (auto const* const fault_elem : felems)
         {
             if (fault_elem->getDimension() != 2)
diff --git a/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.cpp b/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.cpp
index 42f6a41ad657308a24f56db7a1ea529d13197568..04c0d50283ba27d0b91c9f29d36a389a2753a20d 100644
--- a/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.cpp
+++ b/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.cpp
@@ -65,7 +65,7 @@ std::ostream& Partition::writeNodes(
 
     for (const auto* node : nodes)
     {
-        double const* coords = node->getCoords();
+        double const* coords = node->data();
         nodes_buffer.emplace_back(global_node_ids[node->getID()], coords[0],
                                   coords[1], coords[2]);
     }
diff --git a/GeoLib/AnalyticalGeometry-impl.h b/GeoLib/AnalyticalGeometry-impl.h
index 69fbceaf0c751a50642cf0d72115619627b0da4a..6838ed5a81f31c35ca7440ce0b50ac23e414fb77 100644
--- a/GeoLib/AnalyticalGeometry-impl.h
+++ b/GeoLib/AnalyticalGeometry-impl.h
@@ -26,7 +26,7 @@ std::pair<Eigen::Vector3d, double> getNewellPlane(InputIterator pnts_begin,
         plane_normal[2] += (pt_i[0] - pt_j[0])
                            * (pt_i[1] + pt_j[1]); // projection on xy
 
-        centroid += Eigen::Map<Eigen::Vector3d const>(pt_j.getCoords());
+        centroid += pt_j.asEigenVector3d();
     }
 
     plane_normal.normalize();
@@ -62,7 +62,7 @@ std::pair<Eigen::Vector3d, double> getNewellPlane(
         plane_normal[2] += (pnts[i][0] - pnts[j][0])
                            * (pnts[i][1] + pnts[j][1]); // projection on xy
 
-        centroid += Eigen::Map<Eigen::Vector3d const>(pnts[j].getCoords());
+        centroid += pnts[j].asEigenVector3d();
     }
 
     plane_normal.normalize();
@@ -77,8 +77,7 @@ void rotatePoints(Eigen::Matrix3d const& rot_mat, InputIterator pnts_begin,
 {
     for (auto it = pnts_begin; it != pnts_end; ++it)
     {
-        Eigen::Map<Eigen::Vector3d>((*it)->getCoords()) =
-            rot_mat * Eigen::Map<Eigen::Vector3d const>((*it)->getCoords());
+        (*it)->asEigenVector3d() = rot_mat * (*it)->asEigenVector3d();
     }
 }
 
diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index 0fc29e46ddf10bc9be9a68fe7aa3d6161a540bfe..024e15ba8229b11560ce5bda2268f0a49e7021c0 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -32,18 +32,18 @@ namespace ExactPredicates
 double getOrientation2d(MathLib::Point3d const& a, MathLib::Point3d const& b,
                         MathLib::Point3d const& c)
 {
-    return orient2d(const_cast<double*>(a.getCoords()),
-                    const_cast<double*>(b.getCoords()),
-                    const_cast<double*>(c.getCoords()));
+    return orient2d(const_cast<double*>(a.data()),
+                    const_cast<double*>(b.data()),
+                    const_cast<double*>(c.data()));
 }
 
 double getOrientation2dFast(MathLib::Point3d const& a,
                             MathLib::Point3d const& b,
                             MathLib::Point3d const& c)
 {
-    return orient2dfast(const_cast<double*>(a.getCoords()),
-                        const_cast<double*>(b.getCoords()),
-                        const_cast<double*>(c.getCoords()));
+    return orient2dfast(const_cast<double*>(a.data()),
+                        const_cast<double*>(b.data()),
+                        const_cast<double*>(c.data()));
 }
 }  // namespace ExactPredicates
 
@@ -152,14 +152,10 @@ bool lineSegmentIntersect(GeoLib::LineSegment const& s0,
         return false;
     }
 
-    auto const a =
-        Eigen::Map<Eigen::Vector3d const>(s0.getBeginPoint().getCoords());
-    auto const b =
-        Eigen::Map<Eigen::Vector3d const>(s0.getEndPoint().getCoords());
-    auto const c =
-        Eigen::Map<Eigen::Vector3d const>(s1.getBeginPoint().getCoords());
-    auto const d =
-        Eigen::Map<Eigen::Vector3d const>(s1.getEndPoint().getCoords());
+    auto const& a = s0.getBeginPoint().asEigenVector3d();
+    auto const& b = s0.getEndPoint().asEigenVector3d();
+    auto const& c = s1.getBeginPoint().asEigenVector3d();
+    auto const& d = s1.getEndPoint().asEigenVector3d();
 
     Eigen::Vector3d const v = b - a;
     Eigen::Vector3d const w = d - c;
@@ -373,16 +369,10 @@ std::unique_ptr<GeoLib::Point> triangleLineIntersection(
     MathLib::Point3d const& c, MathLib::Point3d const& p,
     MathLib::Point3d const& q)
 {
-    auto const va = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
-    auto const vb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
-    auto const vc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
-    auto const vp = Eigen::Map<Eigen::Vector3d const>(p.getCoords());
-    auto const vq = Eigen::Map<Eigen::Vector3d const>(q.getCoords());
-
-    Eigen::Vector3d const pq = vq - vp;
-    Eigen::Vector3d const pa = va - vp;
-    Eigen::Vector3d const pb = vb - vp;
-    Eigen::Vector3d const pc = vc - vp;
+    Eigen::Vector3d const pq = q.asEigenVector3d() - p.asEigenVector3d();
+    Eigen::Vector3d const pa = a.asEigenVector3d() - p.asEigenVector3d();
+    Eigen::Vector3d const pb = b.asEigenVector3d() - p.asEigenVector3d();
+    Eigen::Vector3d const pc = c.asEigenVector3d() - p.asEigenVector3d();
 
     double u = pq.cross(pc).dot(pb);
     if (u < 0)
diff --git a/GeoLib/Grid.h b/GeoLib/Grid.h
index c8827c2fb7b267b185f32a042af1fb88243c079e..d748837d174334a185b8f88afc7fd1cea9bb4ff6 100644
--- a/GeoLib/Grid.h
+++ b/GeoLib/Grid.h
@@ -532,10 +532,8 @@ POINT* Grid<POINT>::getNearestPoint(P const& pnt) const
         }  // end while
     }      // end else
 
-    auto to_eigen = [](auto const& point)
-    { return Eigen::Map<Eigen::Vector3d const>(point.getCoords()); };
-
-    double len((to_eigen(pnt) - to_eigen(*nearest_pnt)).norm());
+    double const len(
+        (pnt.asEigenVector3d() - nearest_pnt->asEigenVector3d()).norm());
     // search all other grid cells within the cube with the edge nodes
     std::vector<std::vector<POINT*> const*> vecs_of_pnts(
         getPntVecsOfGridCellsIntersectingCube(pnt, len));
@@ -548,7 +546,8 @@ POINT* Grid<POINT>::getNearestPoint(P const& pnt) const
         for (std::size_t k(0); k < n_pnts; k++)
         {
             const double sqr_dist(
-                (to_eigen(pnt) - to_eigen(*pnts[k])).squaredNorm());
+                (pnt.asEigenVector3d() - pnts[k]->asEigenVector3d())
+                    .squaredNorm());
             if (sqr_dist < sqr_min_dist)
             {
                 sqr_min_dist = sqr_dist;
@@ -652,16 +651,14 @@ bool Grid<POINT>::calcNearestPointInGridCell(
     if (pnts.empty())
         return false;
 
-    auto to_eigen = [](auto const& point)
-    { return Eigen::Map<Eigen::Vector3d const>(point.getCoords()); };
-
     const std::size_t n_pnts(pnts.size());
-    sqr_min_dist = (to_eigen(*pnts[0]) - to_eigen(pnt)).squaredNorm();
+    sqr_min_dist =
+        (pnts[0]->asEigenVector3d() - pnt.asEigenVector3d()).squaredNorm();
     nearest_pnt = pnts[0];
     for (std::size_t i(1); i < n_pnts; i++)
     {
         const double sqr_dist(
-            (to_eigen(*pnts[i]) - to_eigen(pnt)).squaredNorm());
+            (pnts[i]->asEigenVector3d() - pnt.asEigenVector3d()).squaredNorm());
         if (sqr_dist < sqr_min_dist)
         {
             sqr_min_dist = sqr_dist;
@@ -681,15 +678,13 @@ std::vector<std::size_t> Grid<POINT>::getPointsInEpsilonEnvironment(
 
     double const sqr_eps(eps * eps);
 
-    auto to_eigen = [](auto const& point)
-    { return Eigen::Map<Eigen::Vector3d const>(point.getCoords()); };
-
     std::vector<std::size_t> pnts;
     for (auto vec : vec_pnts)
     {
         for (auto const p : *vec)
         {
-            if ((to_eigen(*p) - to_eigen(pnt)).squaredNorm() <= sqr_eps)
+            if ((p->asEigenVector3d() - pnt.asEigenVector3d()).squaredNorm() <=
+                sqr_eps)
             {
                 pnts.push_back(p->getID());
             }
diff --git a/GeoLib/MinimalBoundingSphere.cpp b/GeoLib/MinimalBoundingSphere.cpp
index 875c189428170be8cf047e80367dae5ff515322d..3195f3c24235f1493276df010ed8e2a00cdc3e2c 100644
--- a/GeoLib/MinimalBoundingSphere.cpp
+++ b/GeoLib/MinimalBoundingSphere.cpp
@@ -34,13 +34,11 @@ MinimalBoundingSphere::MinimalBoundingSphere(MathLib::Point3d const& p,
                                              MathLib::Point3d const& q)
     : _radius(std::numeric_limits<double>::epsilon()), _center(p)
 {
-    auto const vp = Eigen::Map<Eigen::Vector3d const>(p.getCoords());
-    auto const vq = Eigen::Map<Eigen::Vector3d const>(q.getCoords());
-    Eigen::Vector3d const a = vq - vp;
+    Eigen::Vector3d const a = q.asEigenVector3d() - p.asEigenVector3d();
 
     Eigen::Vector3d o = a / 2;
     _radius = o.norm() + std::numeric_limits<double>::epsilon();
-    o += vp;
+    o += p.asEigenVector3d();
     _center = MathLib::Point3d{{o[0], o[1], o[2]}};
 }
 
@@ -48,11 +46,9 @@ MinimalBoundingSphere::MinimalBoundingSphere(MathLib::Point3d const& p,
                                              MathLib::Point3d const& q,
                                              MathLib::Point3d const& r)
 {
-    auto const vp = Eigen::Map<Eigen::Vector3d const>(p.getCoords());
-    auto const vq = Eigen::Map<Eigen::Vector3d const>(q.getCoords());
-    auto const vr = Eigen::Map<Eigen::Vector3d const>(r.getCoords());
-    Eigen::Vector3d const a = vr - vp;
-    Eigen::Vector3d const b = vq - vp;
+    auto const& vp = p.asEigenVector3d();
+    Eigen::Vector3d const a = r.asEigenVector3d() - vp;
+    Eigen::Vector3d const b = q.asEigenVector3d() - vp;
     Eigen::Vector3d const axb = a.cross(b);
 
     if (axb.squaredNorm() > 0)
@@ -85,14 +81,9 @@ MinimalBoundingSphere::MinimalBoundingSphere(MathLib::Point3d const& p,
                                              MathLib::Point3d const& r,
                                              MathLib::Point3d const& s)
 {
-    auto const vp = Eigen::Map<Eigen::Vector3d const>(p.getCoords());
-    auto const vq = Eigen::Map<Eigen::Vector3d const>(q.getCoords());
-    auto const vr = Eigen::Map<Eigen::Vector3d const>(r.getCoords());
-    auto const vs = Eigen::Map<Eigen::Vector3d const>(s.getCoords());
-
-    Eigen::Vector3d const va = vq - vp;
-    Eigen::Vector3d const vb = vr - vp;
-    Eigen::Vector3d const vc = vs - vp;
+    Eigen::Vector3d const va = q.asEigenVector3d() - p.asEigenVector3d();
+    Eigen::Vector3d const vb = r.asEigenVector3d() - p.asEigenVector3d();
+    Eigen::Vector3d const vc = s.asEigenVector3d() - p.asEigenVector3d();
 
     if (!MathLib::isCoplanar(p, q, r, s))
     {
@@ -103,7 +94,7 @@ MinimalBoundingSphere::MinimalBoundingSphere(MathLib::Point3d const& p,
             denom;
 
         _radius = o.norm() + std::numeric_limits<double>::epsilon();
-        o += vp;
+        o += p.asEigenVector3d();
         _center = MathLib::Point3d({o[0], o[1], o[2]});
     }
     else
diff --git a/GeoLib/OctTree-impl.h b/GeoLib/OctTree-impl.h
index d431313beb0f04c2a4ea795828272ba1a5e221bd..15932f62140a1a037a87c65e96daf01ae934ccea 100644
--- a/GeoLib/OctTree-impl.h
+++ b/GeoLib/OctTree-impl.h
@@ -78,19 +78,16 @@ bool OctTree<POINT, MAX_POINTS>::addPoint(POINT* pnt, POINT*& ret_pnt)
 {
     // first do a range query using a epsilon box around the point pnt
     std::vector<POINT*> query_pnts;
-    Eigen::Vector3d const min =
-        Eigen::Map<Eigen::Vector3d>(pnt->getCoords()).array() - _eps;
-    Eigen::Vector3d const max =
-        Eigen::Map<Eigen::Vector3d>(pnt->getCoords()).array() + _eps;
+    Eigen::Vector3d const min = pnt->asEigenVector3d().array() - _eps;
+    Eigen::Vector3d const max = pnt->asEigenVector3d().array() + _eps;
     getPointsInRange(min, max, query_pnts);
-    auto const it = std::find_if(
-        query_pnts.begin(), query_pnts.end(),
-        [pnt, this](auto const* p)
-        {
-            return (Eigen::Map<Eigen::Vector3d const>(p->getCoords()) -
-                    Eigen::Map<Eigen::Vector3d const>(pnt->getCoords()))
-                       .squaredNorm() < _eps * _eps;
-        });
+    auto const it =
+        std::find_if(query_pnts.begin(), query_pnts.end(),
+                     [pnt, this](auto const* p)
+                     {
+                         return (p->asEigenVector3d() - pnt->asEigenVector3d())
+                                    .squaredNorm() < _eps * _eps;
+                     });
     if (it != query_pnts.end())
     {
         ret_pnt = *it;
diff --git a/GeoLib/Polyline.cpp b/GeoLib/Polyline.cpp
index 2388b4d30cbeac9a8b84ece9b02a1874df40164f..75ebbf69225243cb8447fc81150c3772ce0d533f 100644
--- a/GeoLib/Polyline.cpp
+++ b/GeoLib/Polyline.cpp
@@ -370,10 +370,8 @@ double Polyline::getDistanceAlongPolyline(const MathLib::Point3d& pnt,
     // loop over all line segments of the polyline
     for (std::size_t k = 0; k < getNumberOfSegments(); k++)
     {
-        auto const a =
-            Eigen::Map<Eigen::Vector3d const>(getPoint(k)->getCoords());
-        auto const b =
-            Eigen::Map<Eigen::Vector3d const>(getPoint(k + 1)->getCoords());
+        auto const& a = getPoint(k)->asEigenVector3d();
+        auto const& b = getPoint(k + 1)->asEigenVector3d();
         double const seg_length((b - a).norm());
         act_length_of_ply += seg_length;
         // is the orthogonal projection of the j-th node to the
diff --git a/GeoLib/Utils.cpp b/GeoLib/Utils.cpp
index 91dde5322653d88791b2dd2635a9810563d8fb8f..4772aefdc351f5cddb00162a9caf55437ed580d2 100644
--- a/GeoLib/Utils.cpp
+++ b/GeoLib/Utils.cpp
@@ -26,8 +26,8 @@ std::vector<GeoLib::Point*> generateEquidistantPoints(
             "be non-negative.");
     }
 
-    auto const start = Eigen::Map<Eigen::Vector3d const>(begin.getCoords());
-    auto const stop = Eigen::Map<Eigen::Vector3d const>(end.getCoords());
+    auto const& start = begin.asEigenVector3d();
+    auto const& stop = end.asEigenVector3d();
     auto const delta = (stop - start) / (number_of_subdivisions + 1);
 
     std::vector<GeoLib::Point*> points;
diff --git a/MathLib/GeometricBasics.cpp b/MathLib/GeometricBasics.cpp
index 22c39d5b74239a14662e8b24480fb1bb8d4125a4..d6cc807d2cbe2771e430ec1a6b23efa7c335f45c 100644
--- a/MathLib/GeometricBasics.cpp
+++ b/MathLib/GeometricBasics.cpp
@@ -21,14 +21,9 @@ double orientation3d(MathLib::Point3d const& p,
                      MathLib::Point3d const& b,
                      MathLib::Point3d const& c)
 {
-    auto const pp = Eigen::Map<Eigen::Vector3d const>(p.getCoords());
-    auto const pa = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
-    auto const pb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
-    auto const pc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
-
-    Eigen::Vector3d const u = pp - pa;
-    Eigen::Vector3d const v = pp - pb;
-    Eigen::Vector3d const w = pp - pc;
+    Eigen::Vector3d const u = p.asEigenVector3d() - a.asEigenVector3d();
+    Eigen::Vector3d const v = p.asEigenVector3d() - b.asEigenVector3d();
+    Eigen::Vector3d const w = p.asEigenVector3d() - c.asEigenVector3d();
     return u.cross(v).dot(w);
 }
 
@@ -37,24 +32,17 @@ double calcTetrahedronVolume(MathLib::Point3d const& a,
                              MathLib::Point3d const& c,
                              MathLib::Point3d const& d)
 {
-    auto const va = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
-    auto const vb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
-    auto const vc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
-    auto const vd = Eigen::Map<Eigen::Vector3d const>(d.getCoords());
-    Eigen::Vector3d const w = vb - va;
-    Eigen::Vector3d const u = vc - va;
-    Eigen::Vector3d const v = vd - va;
+    Eigen::Vector3d const w = b.asEigenVector3d() - a.asEigenVector3d();
+    Eigen::Vector3d const u = c.asEigenVector3d() - a.asEigenVector3d();
+    Eigen::Vector3d const v = d.asEigenVector3d() - a.asEigenVector3d();
     return std::abs(u.cross(v).dot(w)) / 6.0;
 }
 
 double calcTriangleArea(MathLib::Point3d const& a, MathLib::Point3d const& b,
                         MathLib::Point3d const& c)
 {
-    auto const va = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
-    auto const vb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
-    auto const vc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
-    Eigen::Vector3d const u = vc - va;
-    Eigen::Vector3d const v = vb - va;
+    Eigen::Vector3d const u = c.asEigenVector3d() - a.asEigenVector3d();
+    Eigen::Vector3d const v = b.asEigenVector3d() - a.asEigenVector3d();
     Eigen::Vector3d const w = u.cross(v);
     return 0.5 * w.norm();
 }
@@ -124,11 +112,9 @@ bool gaussPointInTriangle(MathLib::Point3d const& q,
                           double eps_pnt_out_of_plane,
                           double eps_pnt_out_of_tri)
 {
-    auto const pa = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
-    auto const pb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
-    auto const pc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
-    Eigen::Vector3d const v = pb - pa;
-    Eigen::Vector3d const w = pc - pa;
+    auto const& pa = a.asEigenVector3d();
+    Eigen::Vector3d const v = b.asEigenVector3d() - pa;
+    Eigen::Vector3d const w = c.asEigenVector3d() - pa;
 
     Eigen::Matrix2d mat;
     mat(0, 0) = v.squaredNorm();
@@ -171,13 +157,10 @@ bool barycentricPointInTriangle(MathLib::Point3d const& p,
         return false;
     }
 
-    auto const vp = Eigen::Map<Eigen::Vector3d const>(p.getCoords());
-    auto const va = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
-    auto const vb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
-    auto const vc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
-    Eigen::Vector3d const pa = va - vp;
-    Eigen::Vector3d const pb = vb - vp;
-    Eigen::Vector3d const pc = vc - vp;
+    auto const& vp = p.asEigenVector3d();
+    Eigen::Vector3d const& pa = a.asEigenVector3d() - vp;
+    Eigen::Vector3d const& pb = b.asEigenVector3d() - vp;
+    Eigen::Vector3d const& pc = c.asEigenVector3d() - vp;
     double const area_x_2(calcTriangleArea(a, b, c) * 2);
 
     double const alpha((pb.cross(pc).norm()) / area_x_2);
@@ -236,14 +219,9 @@ bool dividedByPlane(const MathLib::Point3d& a, const MathLib::Point3d& b,
 bool isCoplanar(const MathLib::Point3d& a, const MathLib::Point3d& b,
                 const MathLib::Point3d& c, const MathLib::Point3d& d)
 {
-    auto const pa = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
-    auto const pb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
-    auto const pc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
-    auto const pd = Eigen::Map<Eigen::Vector3d const>(d.getCoords());
-
-    Eigen::Vector3d const ab = pb - pa;
-    Eigen::Vector3d const ac = pc - pa;
-    Eigen::Vector3d const ad = pd - pa;
+    Eigen::Vector3d const ab = b.asEigenVector3d() - a.asEigenVector3d();
+    Eigen::Vector3d const ac = c.asEigenVector3d() - a.asEigenVector3d();
+    Eigen::Vector3d const ad = d.asEigenVector3d() - a.asEigenVector3d();
 
     auto const eps_squared =
         std::pow(std::numeric_limits<double>::epsilon(), 2);
diff --git a/MathLib/MathTools.cpp b/MathLib/MathTools.cpp
index cc8168011930c8781d8050bbd6833b0963f90f68..63b78854b209dd45206eed4982052cff23c1829b 100644
--- a/MathLib/MathTools.cpp
+++ b/MathLib/MathTools.cpp
@@ -19,9 +19,9 @@ namespace MathLib
 double calcProjPntToLineAndDists(Point3d const& pp, Point3d const& pa,
                                  Point3d const& pb, double& lambda, double& d0)
 {
-    auto const a = Eigen::Map<Eigen::Vector3d const>(pa.getCoords());
-    auto const b = Eigen::Map<Eigen::Vector3d const>(pb.getCoords());
-    auto const p = Eigen::Map<Eigen::Vector3d const>(pp.getCoords());
+    auto const& a = pa.asEigenVector3d();
+    auto const& b = pb.asEigenVector3d();
+    auto const& p = pp.asEigenVector3d();
 
     // g(lambda) = a + lambda v, v = b-a
     Eigen::Vector3d const v = b - a;
@@ -41,11 +41,9 @@ double calcProjPntToLineAndDists(Point3d const& pp, Point3d const& pa,
 
 double getAngle(Point3d const& p0, Point3d const& p1, Point3d const& p2)
 {
-    auto const a = Eigen::Map<Eigen::Vector3d const>(p0.getCoords());
-    auto const b = Eigen::Map<Eigen::Vector3d const>(p1.getCoords());
-    auto const c = Eigen::Map<Eigen::Vector3d const>(p2.getCoords());
-    Eigen::Vector3d const v0 = a - b;
-    Eigen::Vector3d const v1 = c - b;
+    auto const& b = p1.asEigenVector3d();
+    Eigen::Vector3d const v0 = p0.asEigenVector3d() - b;
+    Eigen::Vector3d const v1 = p2.asEigenVector3d() - b;
 
     // apply Cauchy Schwarz inequality
     return std::acos(v0.dot(v1) / (v0.norm() * v1.norm()));
diff --git a/MathLib/Point3d.cpp b/MathLib/Point3d.cpp
index 28b6f2394ddcc613eaeb960506595b70809adbba..c6a66613c4ccc24a9e02536d1d94fd116dcbc5ee 100644
--- a/MathLib/Point3d.cpp
+++ b/MathLib/Point3d.cpp
@@ -12,9 +12,35 @@
 
 namespace MathLib
 {
-Point3d::Point3d() : x_({{0}}) {}
+Point3d::Point3d() : x_({0, 0, 0}) {}
 
-Point3d::Point3d(std::array<double, 3> x) : x_(std::move(x)) {}
+Point3d::Point3d(std::array<double, 3> x) : x_(x[0], x[1], x[2]) {}
+
+double sqrDist(MathLib::Point3d const& p0, MathLib::Point3d const& p1)
+{
+    return (p0.asEigenVector3d() - p1.asEigenVector3d()).squaredNorm();
+}
+
+bool lessEq(Point3d const& a, Point3d const& b, double eps)
+{
+    auto absAndRelDiffLargerThanEps = [eps](double const u,
+                                            double const v) -> bool
+    {
+        return std::abs(u - v) > eps * std::min(std::abs(v), std::abs(u)) &&
+               std::abs(u - v) > eps;
+    };
+
+    return std::lexicographical_compare(
+        a.data(), a.data() + 3, b.data(), b.data() + 3,
+        [&absAndRelDiffLargerThanEps](auto const u, auto const v)
+        {
+            if (absAndRelDiffLargerThanEps(u, v))
+            {
+                return u <= v;
+            }
+            return true;
+        });
+}
 
 extern const Point3d ORIGIN{{{0.0, 0.0, 0.0}}};
 }  // namespace MathLib
diff --git a/MathLib/Point3d.h b/MathLib/Point3d.h
index 1375f8c82639337b840897b2eb1b85058679dd13..2fb478b34cf1bc002a4600c810c11c03edc211f9 100644
--- a/MathLib/Point3d.h
+++ b/MathLib/Point3d.h
@@ -60,32 +60,21 @@ public:
     }
 
     /** returns an array containing the coordinates of the point */
-    const double* getCoords() const { return x_.data(); }
+    const double* data() const { return x_.data(); }
 
-    double* getCoords() { return x_.data(); }
+    double* data() { return x_.data(); }
+
+    Eigen::Vector3d const& asEigenVector3d() const { return x_; }
+    Eigen::Vector3d& asEigenVector3d() { return x_; }
 
 private:
-    std::array<double, 3> x_;
+    Eigen::Vector3d x_;
 };
 
 inline bool operator<(Point3d const& a, Point3d const& b)
 {
-    for (std::size_t i = 0; i < 3; ++i)
-    {
-        if (a[i] > b[i])
-        {
-            return false;
-        }
-        if (a[i] < b[i])
-        {
-            return true;
-        }
-
-        // continue with next dimension, because a[0] == b[0]
-    }
-
-    // The values in all dimensions are equal.
-    return false;
+    return std::lexicographical_compare(a.data(), a.data() + 3, b.data(),
+                                        b.data() + 3);
 }
 
 /**
@@ -99,28 +88,8 @@ inline bool operator<(Point3d const& a, Point3d const& b)
  * test \f$ |a_i - b_i| > \epsilon \cdot \min (|a_i|, |b_i|) \f$ \b and
  * \f$  |a_i - b_i| > \epsilon \f$ for all coordinates \f$ 0 \le i < 3 \f$.
  */
-bool inline lessEq(Point3d const& a, Point3d const& b,
-                   double eps = std::numeric_limits<double>::epsilon())
-{
-    auto coordinateIsLargerEps = [&eps](double const u, double const v) -> bool
-    {
-        return std::abs(u - v) > eps * std::min(std::abs(v), std::abs(u)) &&
-               std::abs(u - v) > eps;
-    };
-
-    for (std::size_t i = 0; i < 3; ++i)
-    {
-        // test a relative and an absolute criterion
-        if (coordinateIsLargerEps(a[i], b[i]))
-        {
-            return static_cast<bool>(a[i] <= b[i]);
-        }
-        // a[i] ~= b[i] up to an epsilon. Compare next dimension.
-    }
-
-    // all coordinates are equal up to an epsilon.
-    return true;
-}
+bool lessEq(Point3d const& a, Point3d const& b,
+            double eps = std::numeric_limits<double>::epsilon());
 
 /** overload the output operator for class Point */
 inline std::ostream& operator<<(std::ostream& os, const Point3d& p)
@@ -152,12 +121,7 @@ inline MathLib::Point3d operator*(MATRIX const& mat, MathLib::Point3d const& p)
 
 /** Computes the squared dist between the two points p0 and p1.
  */
-inline double sqrDist(MathLib::Point3d const& p0, MathLib::Point3d const& p1)
-{
-    return (p0[0] - p1[0]) * (p0[0] - p1[0]) +
-           (p0[1] - p1[1]) * (p0[1] - p1[1]) +
-           (p0[2] - p1[2]) * (p0[2] - p1[2]);
-}
+double sqrDist(MathLib::Point3d const& p0, MathLib::Point3d const& p1);
 
 /** Equality of Point3d's up to an epsilon.
  */
diff --git a/MathLib/WeightedPoint.h b/MathLib/WeightedPoint.h
index 959820dd15a42fb6878d7e531508eefdfa577f27..ca1a9fdee466cca9a334d0f0504b479ac016baed 100644
--- a/MathLib/WeightedPoint.h
+++ b/MathLib/WeightedPoint.h
@@ -75,7 +75,7 @@ public:
         return !(*this == other);
     }
 
-    const double* getCoords() const { return coords_.data(); }
+    const double* data() const { return coords_.data(); }
 
     double getWeight() const { return weight_; }
 
diff --git a/MeshGeoToolsLib/GeoMapper.cpp b/MeshGeoToolsLib/GeoMapper.cpp
index a4fb484a2a20759671b274d1f76a592768def640..a98e64859c79bca98a94df1f38a69381c2bb6eab 100644
--- a/MeshGeoToolsLib/GeoMapper.cpp
+++ b/MeshGeoToolsLib/GeoMapper.cpp
@@ -450,8 +450,7 @@ static void mapPointOnSurfaceElement(MeshLib::Element const& elem,
                                      MathLib::Point3d& q)
 {
     // create plane equation: n*p = d
-    auto const p =
-        Eigen::Map<Eigen::Vector3d const>(elem.getNode(0)->getCoords());
+    auto const& p = elem.getNode(0)->asEigenVector3d();
     Eigen::Vector3d const n(MeshLib::FaceRule::getSurfaceNormal(elem));
     if (n[2] == 0.0)
     {  // vertical plane, z coordinate is arbitrary
diff --git a/MeshLib/ElementCoordinatesMappingLocal.cpp b/MeshLib/ElementCoordinatesMappingLocal.cpp
index c771fc5ffd152a5a828c6d7bbcd05bf6c81fa026..bac3131e117d0ae07c59cc14ba6656271f245851 100644
--- a/MeshLib/ElementCoordinatesMappingLocal.cpp
+++ b/MeshLib/ElementCoordinatesMappingLocal.cpp
@@ -42,8 +42,7 @@ MeshLib::RotationMatrix getRotationMatrixToGlobal(
     if (element_dimension == 1)
     {
         Eigen::Vector3d const xx =
-            (Eigen::Map<Eigen::Vector3d const>(points[1].getCoords()) -
-             Eigen::Map<Eigen::Vector3d const>(points[0].getCoords()))
+            (points[1].asEigenVector3d() - points[0].asEigenVector3d())
                 .normalized();
         if (global_dim == 2)
         {
diff --git a/MeshLib/Elements/CellRule.cpp b/MeshLib/Elements/CellRule.cpp
index 5fc840f2bd3f74ab40f6578857bd707e53374478..b5a88a9ef13b269b870e66edb690488f724b1b8c 100644
--- a/MeshLib/Elements/CellRule.cpp
+++ b/MeshLib/Elements/CellRule.cpp
@@ -18,8 +18,7 @@ namespace MeshLib
 {
 bool CellRule::testElementNodeOrder(Element const& e)
 {
-    Eigen::Vector3d const cc =
-        Eigen::Map<Eigen::Vector3d const>(getCenterOfGravity(e).getCoords());
+    auto const cc = getCenterOfGravity(e).asEigenVector3d();
     const unsigned nFaces(e.getNumberOfFaces());
     for (unsigned j = 0; j < nFaces; ++j)
     {
@@ -28,8 +27,7 @@ bool CellRule::testElementNodeOrder(Element const& e)
         // test at some point, while for node 0 at least one node in every
         // element type would be used for checking twice and one wouldn't be
         // checked at all. (based on the definition of the _face_nodes variable)
-        auto const x =
-            Eigen::Map<Eigen::Vector3d const>(face->getNode(1)->getCoords());
+        auto const& x = face->getNode(1)->asEigenVector3d();
         Eigen::Vector3d const cx = x - cc;
         const double s = FaceRule::getSurfaceNormal(*face).dot(cx);
         delete face;
diff --git a/MeshLib/Elements/FaceRule.cpp b/MeshLib/Elements/FaceRule.cpp
index d0955bdd564e25bae46d3557402e97944a2ea542..ad66b233b70a38317596ca572514d9ae2384a4b0 100644
--- a/MeshLib/Elements/FaceRule.cpp
+++ b/MeshLib/Elements/FaceRule.cpp
@@ -23,17 +23,15 @@ bool FaceRule::testElementNodeOrder(Element const& e)
 
 Eigen::Vector3d FaceRule::getFirstSurfaceVector(Element const& e)
 {
-    auto const a = Eigen::Map<Eigen::Vector3d const>(e.getNode(0)->getCoords());
-    auto const b = Eigen::Map<Eigen::Vector3d const>(e.getNode(1)->getCoords());
-    Eigen::Vector3d const v = a - b;
+    Eigen::Vector3d const v =
+        e.getNode(0)->asEigenVector3d() - e.getNode(1)->asEigenVector3d();
     return v;
 }
 
 Eigen::Vector3d FaceRule::getSecondSurfaceVector(Element const& e)
 {
-    auto const a = Eigen::Map<Eigen::Vector3d const>(e.getNode(1)->getCoords());
-    auto const b = Eigen::Map<Eigen::Vector3d const>(e.getNode(2)->getCoords());
-    Eigen::Vector3d const v = b - a;
+    Eigen::Vector3d const v =
+        e.getNode(2)->asEigenVector3d() - e.getNode(1)->asEigenVector3d();
     return v;
 }
 
diff --git a/MeshLib/Elements/Utils.h b/MeshLib/Elements/Utils.h
index d3916ad4072408e5517696ef5eafdb8a99eaac49..6fa00f580cbbaa0d0e3188a9b7844899ce7afe17 100644
--- a/MeshLib/Elements/Utils.h
+++ b/MeshLib/Elements/Utils.h
@@ -50,10 +50,8 @@ inline Eigen::Vector3d calculateNormalizedSurfaceNormal(
     {
         auto const bulk_element_normal =
             MeshLib::FaceRule::getSurfaceNormal(bulk_element);
-        auto const v0 = Eigen::Map<Eigen::Vector3d const>(
-            surface_element.getNode(0)->getCoords());
-        auto const v1 = Eigen::Map<Eigen::Vector3d const>(
-            surface_element.getNode(1)->getCoords());
+        auto const& v0 = surface_element.getNode(0)->asEigenVector3d();
+        auto const& v1 = surface_element.getNode(1)->asEigenVector3d();
         Eigen::Vector3d const edge_vector = v1 - v0;
         surface_element_normal = bulk_element_normal.cross(edge_vector);
     }
diff --git a/MeshLib/IO/XDMF/transformData.cpp b/MeshLib/IO/XDMF/transformData.cpp
index f4a52f033556c47410035e80b446703641ca343d..09f2e328db72c078ce37c6d2d6252425a23d8ded 100644
--- a/MeshLib/IO/XDMF/transformData.cpp
+++ b/MeshLib/IO/XDMF/transformData.cpp
@@ -238,7 +238,7 @@ std::vector<double> transformToXDMFGeometry(MeshLib::Mesh const& mesh)
     values.reserve(nodes.size() * point_size);
     for (auto const& n : nodes)
     {
-        const double* x = n->getCoords();
+        const double* x = n->data();
         values.insert(values.cend(), x, x + point_size);
     }
 
diff --git a/MeshLib/MeshEditing/DuplicateMeshComponents.cpp b/MeshLib/MeshEditing/DuplicateMeshComponents.cpp
index 9f90bc02bd5c924b29844d53481ee306ec2a7813..96c3b7488f3451f4f29ea1184afb0755a401e4e0 100644
--- a/MeshLib/MeshEditing/DuplicateMeshComponents.cpp
+++ b/MeshLib/MeshEditing/DuplicateMeshComponents.cpp
@@ -26,7 +26,7 @@ std::vector<Node*> copyNodeVector(const std::vector<Node*>& nodes)
     new_nodes.reserve(nNodes);
     for (std::size_t k = 0; k < nNodes; ++k)
     {
-        new_nodes.push_back(new Node(nodes[k]->getCoords(), new_nodes.size()));
+        new_nodes.push_back(new Node(nodes[k]->data(), new_nodes.size()));
     }
     return new_nodes;
 }
diff --git a/MeshLib/MeshEditing/ProjectPointOnMesh.cpp b/MeshLib/MeshEditing/ProjectPointOnMesh.cpp
index 3827355b5d9e2facf186d6d5d6fddda0acd4f023..095bb9d63b48d10510af9b716cd5f9a2fadc1469 100644
--- a/MeshLib/MeshEditing/ProjectPointOnMesh.cpp
+++ b/MeshLib/MeshEditing/ProjectPointOnMesh.cpp
@@ -58,8 +58,7 @@ Element const* getProjectedElement(std::vector<const Element*> const& elements,
 double getElevation(Element const& element, Node const& node)
 {
     Eigen::Vector3d const v =
-        Eigen::Map<Eigen::Vector3d const>(node.getCoords()) -
-        Eigen::Map<Eigen::Vector3d const>(element.getNode(0)->getCoords());
+        node.asEigenVector3d() - element.getNode(0)->asEigenVector3d();
     auto const n = FaceRule::getSurfaceNormal(element).normalized();
     return node[2] - n.dot(v) * n[2];
 }
diff --git a/MeshLib/MeshGenerators/LayeredVolume.cpp b/MeshLib/MeshGenerators/LayeredVolume.cpp
index af10eb32a391fbc349ae92ad380c52c73dc53542..f9c2c71cb3bc71091d698f5b13c26de430a90c5c 100644
--- a/MeshLib/MeshGenerators/LayeredVolume.cpp
+++ b/MeshLib/MeshGenerators/LayeredVolume.cpp
@@ -155,14 +155,10 @@ void LayeredVolume::addLayerBoundaries(const MeshLib::Mesh& layer,
                     MeshLib::Node* n3 =
                         _nodes[offset + nNodes + getNodeIndex(*elem, i)];
 
-                    auto const v0 =
-                        Eigen::Map<Eigen::Vector3d const>(n0->getCoords());
-                    auto const v1 =
-                        Eigen::Map<Eigen::Vector3d const>(n1->getCoords());
-                    auto const v2 =
-                        Eigen::Map<Eigen::Vector3d const>(n2->getCoords());
-                    auto const v3 =
-                        Eigen::Map<Eigen::Vector3d const>(n3->getCoords());
+                    auto const& v0 = n0->asEigenVector3d();
+                    auto const& v1 = n1->asEigenVector3d();
+                    auto const& v2 = n2->asEigenVector3d();
+                    auto const& v3 = n3->asEigenVector3d();
                     double const eps = std::numeric_limits<double>::epsilon();
 
                     if ((v2 - v1).norm() > eps)
diff --git a/MeshLib/Utils/GetSpaceDimension.cpp b/MeshLib/Utils/GetSpaceDimension.cpp
index b1ff1b29882a55b49ee51969e636b8e097179777..a0c9f3d673365e6568afb656a922fa441c3882af 100644
--- a/MeshLib/Utils/GetSpaceDimension.cpp
+++ b/MeshLib/Utils/GetSpaceDimension.cpp
@@ -24,10 +24,10 @@ int getSpaceDimension(std::vector<Node*> const& nodes)
 {
     std::array x_magnitude = {0.0, 0.0, 0.0};
 
-    double const* const x_ref = nodes[0]->getCoords();
+    double const* const x_ref = nodes[0]->data();
     for (auto const& node : nodes)
     {
-        auto const x = node->getCoords();
+        auto const x = node->data();
         for (int i = 0; i < 3; i++)
         {
             x_magnitude[i] += std::fabs(x[i] - x_ref[i]);
diff --git a/MeshLib/Utils/Is2DMeshOnRotatedVerticalPlane.cpp b/MeshLib/Utils/Is2DMeshOnRotatedVerticalPlane.cpp
index 98952108874416ee03f1cdb1ec06da1f2a61885a..2180dfeaa0f7af5f96bdf3c3b95b1df08c8e5156 100644
--- a/MeshLib/Utils/Is2DMeshOnRotatedVerticalPlane.cpp
+++ b/MeshLib/Utils/Is2DMeshOnRotatedVerticalPlane.cpp
@@ -49,9 +49,9 @@ bool is2DMeshOnRotatedVerticalPlane(Mesh const& mesh)
         [](auto const& element)
         {
             // 3 nodes are enough to make up a plane.
-            auto const x1 = element->getNode(0)->getCoords();
-            auto const x2 = element->getNode(1)->getCoords();
-            auto const x3 = element->getNode(2)->getCoords();
+            auto const x1 = element->getNode(0)->data();
+            auto const x2 = element->getNode(1)->data();
+            auto const x3 = element->getNode(2)->data();
 
             double const a0 = x2[0] - x1[0];
             double const a2 = x2[2] - x1[2];
@@ -68,9 +68,9 @@ bool is2DMeshOnRotatedVerticalPlane(Mesh const& mesh)
         [](auto const& element)
         {
             // 3 nodes are enough to make up a plane.
-            auto const x1 = element->getNode(0)->getCoords();
-            auto const x2 = element->getNode(1)->getCoords();
-            auto const x3 = element->getNode(2)->getCoords();
+            auto const x1 = element->getNode(0)->data();
+            auto const x2 = element->getNode(1)->data();
+            auto const x3 = element->getNode(2)->data();
 
             double const a0 = x2[0] - x1[0];
             double const a1 = x2[1] - x1[1];
diff --git a/MeshLib/convertMeshToGeo.cpp b/MeshLib/convertMeshToGeo.cpp
index f1d98e9e95fa7f16edc936c1bf9f49a5214472ee..c75b8fdfa9b728ea3965a4c274883d3ca5116bc9 100644
--- a/MeshLib/convertMeshToGeo.cpp
+++ b/MeshLib/convertMeshToGeo.cpp
@@ -161,7 +161,7 @@ MeshLib::Mesh* convertSurfaceToMesh(const GeoLib::Surface& sfc,
         for (unsigned j = 0; j < 3; j++)
         {
             tri_nodes[j] =
-                new MeshLib::Node(tri->getPoint(j)->getCoords(), nodeId++);
+                new MeshLib::Node(tri->getPoint(j)->data(), nodeId++);
         }
         elements.push_back(new MeshLib::Tri(tri_nodes, i));
         for (unsigned j = 0; j < 3; j++)
diff --git a/NumLib/Fem/InitShapeMatrices.h b/NumLib/Fem/InitShapeMatrices.h
index 772e492b97aadbf339feb8b08e7f768a7ea2e4dd..7e076ee4e82fae8c02150bbd717a624a29e0f1c2 100644
--- a/NumLib/Fem/InitShapeMatrices.h
+++ b/NumLib/Fem/InitShapeMatrices.h
@@ -39,8 +39,7 @@ computeShapeMatrices(MeshLib::Element const& e, bool const is_axially_symmetric,
         shape_matrices.emplace_back(ShapeFunction::DIM, GlobalDim,
                                     ShapeFunction::NPOINTS);
         fe.template computeShapeFunctions<SelectedShapeMatrixType>(
-            p.getCoords(), shape_matrices.back(), GlobalDim,
-            is_axially_symmetric);
+            p.data(), shape_matrices.back(), GlobalDim, is_axially_symmetric);
     }
 
     return shape_matrices;
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/NormalTractionBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryConditionAndSourceTerm/NormalTractionBoundaryConditionLocalAssembler.h
index 606f8775484de7eb0f3371d6c9a73849184b8582..5c5149004d31f358fd1378b9f6af092be00ad6d6 100644
--- a/ProcessLib/BoundaryConditionAndSourceTerm/NormalTractionBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/NormalTractionBoundaryConditionLocalAssembler.h
@@ -85,9 +85,8 @@ public:
         // TODO Extend to rotated 2d meshes and line elements.
         if (e.getGeomType() == MeshLib::MeshElemType::LINE)
         {
-            Eigen::Vector3d const v1 =
-                Eigen::Map<Eigen::Vector3d const>(e.getNode(1)->getCoords()) -
-                Eigen::Map<Eigen::Vector3d const>(e.getNode(0)->getCoords());
+            Eigen::Vector3d const v1 = e.getNode(1)->asEigenVector3d() -
+                                       e.getNode(0)->asEigenVector3d();
             element_normal[0] = -v1[1];
             element_normal[1] = v1[0];
             element_normal.normalize();
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/Python/PythonBoundaryCondition.cpp b/ProcessLib/BoundaryConditionAndSourceTerm/Python/PythonBoundaryCondition.cpp
index 278064fff9964ea3df250477fb7a225f1b5506da..a8621f77c746605ddbf131917571a9a2eb8d35d7 100644
--- a/ProcessLib/BoundaryConditionAndSourceTerm/Python/PythonBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/Python/PythonBoundaryCondition.cpp
@@ -129,10 +129,9 @@ void PythonBoundaryCondition::getEssentialBCValues(
 
         collectPrimaryVariables(primary_variables, *boundary_node, x);
 
-        auto* coords = boundary_node->getCoords();
-        auto const [apply_bc, bc_value] =
-            bc_object->getDirichletBCValue(t, {coords[0], coords[1], coords[2]},
-                                           boundary_node_id, primary_variables);
+        auto const [apply_bc, bc_value] = bc_object->getDirichletBCValue(
+            t, {(*boundary_node)[0], (*boundary_node)[1], (*boundary_node)[2]},
+            boundary_node_id, primary_variables);
 
         if (!bc_object->isOverriddenEssential())
         {
diff --git a/ProcessLib/DeactivatedSubdomain.cpp b/ProcessLib/DeactivatedSubdomain.cpp
index 3060846612a1402d586fa44377446d6b484186e9..0df6040cbdf038313ddaadeee13e18cab0b2f682 100644
--- a/ProcessLib/DeactivatedSubdomain.cpp
+++ b/ProcessLib/DeactivatedSubdomain.cpp
@@ -66,9 +66,8 @@ bool DeactivatedSubdomain::isDeactivated(MathLib::Point3d const& point,
 
     // Position r on the line at given time.
     Eigen::Vector3d const r = a + t * time_interval.getValue(time);
-    Eigen::Map<Eigen::Vector3d const> const p{point.getCoords(), 3};
 
     // Return true if p is "behind" the plane through r.
-    return (p - r).dot(t) <= 0;
+    return (point.asEigenVector3d() - r).dot(t) <= 0;
 }
 }  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
index 3db01406e44fa954bbe3e6685babe6573ec9405d..9c27859b3de308568f8e48d269bc300c78b0ff88 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
@@ -58,10 +58,8 @@ HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHEType>::
     }
 
     // calculate the element direction vector
-    auto const p0 =
-        Eigen::Map<Eigen::Vector3d const>(e.getNode(0)->getCoords(), 3);
-    auto const p1 =
-        Eigen::Map<Eigen::Vector3d const>(e.getNode(1)->getCoords(), 3);
+    auto const& p0 = e.getNode(0)->asEigenVector3d();
+    auto const& p1 = e.getNode(1)->asEigenVector3d();
 
     _element_direction = (p1 - p0).normalized();
 
diff --git a/ProcessLib/LIE/Common/BranchProperty.h b/ProcessLib/LIE/Common/BranchProperty.h
index 1b4fd08e3e11db15246e01b0230be1ef0eb4eaae..7e64e4daded3f315b7edfe2a587c7d5d7537fc7d 100644
--- a/ProcessLib/LIE/Common/BranchProperty.h
+++ b/ProcessLib/LIE/Common/BranchProperty.h
@@ -23,7 +23,7 @@ struct BranchProperty final
     BranchProperty(MeshLib::Node const& branchNode,
                    int const master_fracture_id_,
                    int const slave_fracture_id_)
-        : coords{branchNode.getCoords()},
+        : coords{branchNode.data()},
           node_id{branchNode.getID()},
           master_fracture_id{master_fracture_id_},
           slave_fracture_id{slave_fracture_id_}
diff --git a/ProcessLib/LIE/Common/FractureProperty.h b/ProcessLib/LIE/Common/FractureProperty.h
index a1fe3eb09bfb8092a7ea5750cf25587df283cdf3..1426ebd63364a96343e8bd78bb87480e7af9f02b 100644
--- a/ProcessLib/LIE/Common/FractureProperty.h
+++ b/ProcessLib/LIE/Common/FractureProperty.h
@@ -81,7 +81,7 @@ inline void setFractureProperty(int const dim, MeshLib::Element const& e,
     // a fracture is not curving
     for (int j = 0; j < 3; j++)
     {
-        frac_prop.point_on_fracture[j] = getCenterOfGravity(e).getCoords()[j];
+        frac_prop.point_on_fracture[j] = getCenterOfGravity(e).data()[j];
     }
 
     const MeshLib::ElementCoordinatesMappingLocal ele_local_coord(e, dim);
diff --git a/ProcessLib/LIE/Common/JunctionProperty.h b/ProcessLib/LIE/Common/JunctionProperty.h
index b89219722aa0a2e8af77a618d4a47084903e51b4..68d4bf5ed1e7933c2b058fbd4bd5a408713183f7 100644
--- a/ProcessLib/LIE/Common/JunctionProperty.h
+++ b/ProcessLib/LIE/Common/JunctionProperty.h
@@ -24,7 +24,7 @@ struct JunctionProperty final
     JunctionProperty(int const junction_id_,
                      MeshLib::Node const& junctionNode,
                      std::array<int, 2> const fracture_ids_)
-        : coords{junctionNode.getCoords()},
+        : coords{junctionNode.data()},
           node_id{junctionNode.getID()},
           fracture_ids{fracture_ids_},
           junction_id{junction_id_}
diff --git a/ProcessLib/LIE/Common/PostUtils.cpp b/ProcessLib/LIE/Common/PostUtils.cpp
index 84cfcd23b021f5575dd0001d3d8952367276f599..38cc4de63e38e7bd0be6078a452eadeac5811cae 100644
--- a/ProcessLib/LIE/Common/PostUtils.cpp
+++ b/ProcessLib/LIE/Common/PostUtils.cpp
@@ -72,7 +72,7 @@ PostProcessTool::PostProcessTool(
         for (auto const* org_node : vec_fracture_nodes)
         {
             auto duplicated_node =
-                new MeshLib::Node(org_node->getCoords(), new_nodes.size());
+                new MeshLib::Node(org_node->data(), new_nodes.size());
             new_nodes.push_back(duplicated_node);
             _map_dup_newNodeIDs[org_node->getID()].push_back(
                 duplicated_node->getID());
@@ -83,7 +83,7 @@ PostProcessTool::PostProcessTool(
     {
         auto* org_node = org_mesh.getNode(entry.first);
         auto duplicated_node =
-            new MeshLib::Node(org_node->getCoords(), new_nodes.size());
+            new MeshLib::Node(org_node->data(), new_nodes.size());
         new_nodes.push_back(duplicated_node);
         _map_dup_newNodeIDs[org_node->getID()].push_back(
             duplicated_node->getID());
diff --git a/ProcessLib/LIE/Common/Utils.h b/ProcessLib/LIE/Common/Utils.h
index bbf99f9fcba9ee579009dee439d97f50ea04a41c..45eb015a4b6e7ef8619321952f72fdae77a0d0fe 100644
--- a/ProcessLib/LIE/Common/Utils.h
+++ b/ProcessLib/LIE/Common/Utils.h
@@ -21,10 +21,10 @@ namespace LIE
 /// compute physical coordinates from the given shape vector, i.e. from the
 /// natural coordinates
 template <typename Derived>
-MathLib::Point3d computePhysicalCoordinates(
+Eigen::Vector3d computePhysicalCoordinates(
     MeshLib::Element const& e, Eigen::MatrixBase<Derived> const& shape)
 {
-    MathLib::Point3d pt;
+    Eigen::Vector3d pt = Eigen::Vector3d::Zero();
     for (unsigned i = 0; i < e.getNumberOfNodes(); i++)
     {
         MeshLib::Node const& node = *e.getNode(i);
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
index 67bd5710537b6be00dd42f6d0b66c04b1f193b6d..92fcd43173aee7d0f1f8e4e2f276ec5fcae1ee71 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
@@ -313,7 +313,7 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
             std::unordered_map<int, int> fracID_to_local({{0, 0}});
             std::vector<double> levelsets = uGlobalEnrichments(
                 fracture_props, junction_props, fracID_to_local,
-                Eigen::Vector3d(MeshLib::getCenterOfGravity(*e).getCoords()));
+                Eigen::Vector3d(MeshLib::getCenterOfGravity(*e).data()));
             (*mesh_prop_levelset)[e->getID()] = levelsets[0];
         }
 
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrixNearFracture-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrixNearFracture-impl.h
index e13e76ea88e25b04917fa09a71c04e9f7e042d42..451f06ede42fbac5d11576dc06ba67a2af8fbbd4 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrixNearFracture-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrixNearFracture-impl.h
@@ -36,7 +36,7 @@ HydroMechanicsLocalAssemblerMatrixNearFracture<ShapeFunctionDisplacement,
                                          IntegrationMethod, GlobalDim>(
           e, n_variables, local_matrix_size, dofIndex_to_localIndex,
           is_axially_symmetric, integration_order, process_data),
-      _e_center_coords(getCenterOfGravity(e).getCoords())
+      _e_center_coords(getCenterOfGravity(e).data())
 {
     // currently not supporting multiple fractures
     _fracture_props.push_back(process_data.fracture_property.get());
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h
index 6f5e79d90c30977d3d6305ed8b0d3a6a0df66399..65fd2a21862da15f54344c2542bd2115a5014065 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h
@@ -191,8 +191,7 @@ void SmallDeformationLocalAssemblerFracture<
         auto& state = *ip_data.material_state_variables;
         auto& N = _secondary_data.N[ip];
 
-        Eigen::Vector3d const ip_physical_coords(
-            computePhysicalCoordinates(_element, N).getCoords());
+        auto const ip_physical_coords(computePhysicalCoordinates(_element, N));
         std::vector<double> const levelsets(duGlobalEnrichments(
             _fracture_property->fracture_id, _fracture_props, _junction_props,
             _fracID_to_local, ip_physical_coords));
@@ -291,8 +290,7 @@ void SmallDeformationLocalAssemblerFracture<ShapeFunction, IntegrationMethod,
         auto& b_m = ip_data.aperture;
         auto& N = _secondary_data.N[ip];
 
-        Eigen::Vector3d const ip_physical_coords(
-            computePhysicalCoordinates(_element, N).getCoords());
+        auto const ip_physical_coords(computePhysicalCoordinates(_element, N));
         std::vector<double> const levelsets(duGlobalEnrichments(
             _fracture_property->fracture_id, _fracture_props, _junction_props,
             _fracID_to_local, ip_physical_coords));
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
index a4494f5cd435a4ca33936c5af88afb6c0c760b0b..93e568f2ea8f255646e558a36a7d4e60f798e350 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
@@ -211,7 +211,7 @@ void SmallDeformationLocalAssemblerMatrixNearFracture<
 
         // levelset functions
         Eigen::Vector3d const ip_physical_coords(
-            computePhysicalCoordinates(_element, N).getCoords());
+            computePhysicalCoordinates(_element, N).data());
         std::vector<double> const levelsets(
             uGlobalEnrichments(_fracture_props, _junction_props,
                                _fracID_to_local, ip_physical_coords));
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
index f185369cfe5e8051cbdf380776aa45a54ef21468..38ccd6f7cebf7f19b8f4e8dc4a8cfe63ff48e834 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
@@ -423,7 +423,7 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
             continue;
         }
 
-        Eigen::Vector3d const pt(getCenterOfGravity(*e).getCoords());
+        Eigen::Vector3d const pt(getCenterOfGravity(*e).asEigenVector3d());
         std::vector<FractureProperty*> e_fracture_props;
         std::unordered_map<int, int> e_fracID_to_local;
         unsigned tmpi = 0;
diff --git a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h
index 697318e7a2902ea12886f59d2a9294551f0899f4..a1677ddc9eaa9aed641d210bb37349f08ac4d697 100644
--- a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h
+++ b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h
@@ -284,8 +284,7 @@ public:
         auto* nodes = _element.getNodes();
         for (int i = 0; i < N.size(); ++i)
         {
-            Eigen::Map<Eigen::Vector3d const> node_coordinates{
-                nodes[i]->getCoords(), 3};
+            auto const& node_coordinates{nodes[i]->asEigenVector3d()};
             xyz += node_coordinates * N[i];
         }
         return xyz;
diff --git a/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h b/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
index 241242f25b269db58a3220300f3dcc867151b6a8..794dcb5d738ee0c173436e46cb8ac0963976d1cd 100644
--- a/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
+++ b/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
@@ -158,10 +158,8 @@ private:
         {
             auto const bulk_normal =
                 MeshLib::FaceRule::getSurfaceNormal(bulk_element);
-            auto const l0 = Eigen::Map<Eigen::Vector3d const>(
-                surface_element.getNode(0)->getCoords());
-            auto const l1 = Eigen::Map<Eigen::Vector3d const>(
-                surface_element.getNode(1)->getCoords());
+            auto const& l0 = surface_element.getNode(0)->asEigenVector3d();
+            auto const& l1 = surface_element.getNode(1)->asEigenVector3d();
             Eigen::Vector3d const line = l1 - l0;
             surface_element_normal = line.cross(bulk_normal);
         }
diff --git a/Tests/MeshLib/ConvertToLinearMesh.cpp b/Tests/MeshLib/ConvertToLinearMesh.cpp
index 1138c0778b47ebd506e82e6806cd2e06187520e5..f0065a12d07e32ebe770b8cb9a1be6d86e57dbb0 100644
--- a/Tests/MeshLib/ConvertToLinearMesh.cpp
+++ b/Tests/MeshLib/ConvertToLinearMesh.cpp
@@ -165,7 +165,7 @@ std::vector<Node*> permuteNodes(std::vector<std::size_t> const& permutation,
     for (std::size_t i = 0; i < permutation.size(); ++i)
     {
         auto const node_p_i = nodes[permutation[i]];
-        permuted_nodes.push_back(new Node{node_p_i->getCoords(), i});
+        permuted_nodes.push_back(new Node{node_p_i->data(), i});
     }
     return permuted_nodes;
 }
diff --git a/Tests/MeshLib/TestCoordinatesMappingLocal.cpp b/Tests/MeshLib/TestCoordinatesMappingLocal.cpp
index 5def7fd602dbb201d4d2c1f80f5caccd8eff399e..e94a48ee7a4d43c33c13984f7ccf4e5211aa7d45 100644
--- a/Tests/MeshLib/TestCoordinatesMappingLocal.cpp
+++ b/Tests/MeshLib/TestCoordinatesMappingLocal.cpp
@@ -136,13 +136,12 @@ void debugOutput(MeshLib::Element *ele, MeshLib::ElementCoordinatesMappingLocal
 #endif
 
 // check if using the rotation matrix results in the original coordinates
-#define CHECK_COORDS(ele, mapping)                                           \
-    for (unsigned ii = 0; ii < (ele)->getNumberOfNodes(); ii++)              \
-    {                                                                        \
-        MathLib::Point3d global(matR*(mapping).getMappedCoordinates(ii));    \
-        const double eps(std::numeric_limits<double>::epsilon());            \
-        ASSERT_ARRAY_NEAR(&(*(ele)->getNode(ii))[0], global.getCoords(), 3u, \
-                          eps);                                              \
+#define CHECK_COORDS(ele, mapping)                                            \
+    for (unsigned ii = 0; ii < (ele)->getNumberOfNodes(); ii++)               \
+    {                                                                         \
+        MathLib::Point3d global(matR*(mapping).getMappedCoordinates(ii));     \
+        const double eps(std::numeric_limits<double>::epsilon());             \
+        ASSERT_ARRAY_NEAR(&(*(ele)->getNode(ii))[0], global.data(), 3u, eps); \
     }
 
 }  // namespace
diff --git a/Tests/MeshLib/TestGetElementRotationMatrices.cpp b/Tests/MeshLib/TestGetElementRotationMatrices.cpp
index bf3c6434188c50aff4a793484bc36bc4f42a632b..508dc8b55bc5787085aa36749acc9376e0744a83 100644
--- a/Tests/MeshLib/TestGetElementRotationMatrices.cpp
+++ b/Tests/MeshLib/TestGetElementRotationMatrices.cpp
@@ -101,8 +101,8 @@ TEST(MeshLib, GetElementRotationMatrices3DMesh)
         element_rotation_matrices[3].transpose() * b;
     EXPECT_EQ(elements[3]->getDimension(), b_local_1D.size());
 
-    double const* const x_6 = nodes[6]->getCoords();
-    double const* const x_8 = nodes[8]->getCoords();
+    double const* const x_6 = nodes[6]->data();
+    double const* const x_8 = nodes[8]->data();
     // b_local_1D = |b| (x_8-x_6) * b/(|(x_8-x_6)| |b|)
     double dx[3];
     for (int i = 0; i < 3; i++)
@@ -181,8 +181,8 @@ TEST(MeshLib, GetElementRotationMatrices2DMesh)
         element_rotation_matrices[2].transpose() * b;
     EXPECT_EQ(elements[2]->getDimension(), b_local_1D.size());
 
-    double const* const x_0 = nodes[0]->getCoords();
-    double const* const x_2 = nodes[2]->getCoords();
+    double const* const x_0 = nodes[0]->data();
+    double const* const x_2 = nodes[2]->data();
     // b_local_1D = |b| (x_2-x_0) * b/(|(x_2-x_0)| |b|)
     double dx[2];
     for (int i = 0; i < 2; i++)
diff --git a/Tests/MeshLib/TestMapBulkElementPoint.cpp b/Tests/MeshLib/TestMapBulkElementPoint.cpp
index f5ad1fa1600b8ac5a28626fe2aab71571d74656e..f00e72052a93773ebdbf75d1c545a1f2f739d0b2 100644
--- a/Tests/MeshLib/TestMapBulkElementPoint.cpp
+++ b/Tests/MeshLib/TestMapBulkElementPoint.cpp
@@ -79,9 +79,8 @@ static Eigen::MatrixXd getSelectedNodes(Coords const& natural_coordss,
 static Eigen::MatrixXd appendNode(Eigen::MatrixXd const& mat,
                                   MathLib::Point3d const& pt)
 {
-    Eigen::Map<const Eigen::Vector3d> col(pt.getCoords());
     Eigen::MatrixXd nodes(mat.rows(), mat.cols() + 1);
-    nodes << mat, col;
+    nodes << mat, pt.asEigenVector3d();
     return nodes;
 }
 
diff --git a/Tests/NumLib/TestFunctionInterpolation.cpp b/Tests/NumLib/TestFunctionInterpolation.cpp
index 10f246e05192dcae0773024f4afb6907db59bc67..c449a61a4f94f6429f74b6ce23aa3cf894402bd0 100644
--- a/Tests/NumLib/TestFunctionInterpolation.cpp
+++ b/Tests/NumLib/TestFunctionInterpolation.cpp
@@ -65,7 +65,7 @@ TEST(NumLibFunctionInterpolationTest, Linear1DElement)
         ShapeFunction::DIM, ShapeFunction::DIM, ShapeFunction::NPOINTS);
 
     finite_element.computeShapeFunctions(
-        integration_method.getWeightedPoint(0).getCoords(), shape_matrix,
+        integration_method.getWeightedPoint(0).data(), shape_matrix,
         ShapeFunction::DIM, false);
     ASSERT_EQ(2, shape_matrix.N.size());