From 2a483a61e70d3b8642c66229630a9171d8816b79 Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Fri, 9 Jul 2021 17:15:16 +0200
Subject: [PATCH] [MeL] Mv neighbor information from Node to Mesh.

---
 Applications/Utils/MeshEdit/queryMesh.cpp     |  11 +-
 .../IntegrateBoreholesIntoMesh.cpp            |   9 +-
 .../PartitionMesh/NodeWiseMeshPartitioner.cpp |  20 +-
 .../ModelPreparation/createNeumannBc.cpp      |   3 +-
 MeshGeoToolsLib/GeoMapper.cpp                 |   4 +-
 MeshGeoToolsLib/IdentifySubdomainMesh.cpp     |   3 +-
 MeshGeoToolsLib/MeshNodesOnPoint.cpp          |   3 +-
 MeshLib/ElementStatus.cpp                     |   6 +-
 MeshLib/Mesh.cpp                              | 102 +++++++---
 MeshLib/Mesh.h                                |  14 +-
 .../Mesh2MeshPropertyInterpolation.cpp        |  11 +-
 MeshLib/MeshSearch/ElementSearch.cpp          |   2 +-
 MeshLib/MeshSearch/NodeSearch.cpp             |   7 +-
 MeshLib/MeshSurfaceExtraction.cpp             |   2 +-
 MeshLib/Node.cpp                              |  22 --
 MeshLib/Node.h                                |  30 +--
 MeshLib/NodePartitionedMesh.h                 |   4 +-
 NumLib/DOF/MeshComponentMap.cpp               |   4 +-
 .../DeactivatedSubdomainDirichlet.cpp         |   6 +-
 .../DeactivatedSubdomainDirichlet.h           |   3 -
 ProcessLib/CreateDeactivatedSubdomain.cpp     |   5 +-
 .../HeatTransportBHEProcess.cpp               |   4 +-
 ProcessLib/LIE/Common/MeshUtils.cpp           |  35 ++--
 .../SmallDeformationProcess.cpp               |   7 +-
 .../MeshLib/TestFindElementsWithinRadius.cpp  |  16 +-
 Tests/MeshLib/TestLineMesh.cpp                |  16 +-
 Tests/MeshLib/TestNodeAdjacencyTable.cpp      |   9 +-
 Tests/MeshLib/TestQuadMesh.cpp                |  62 +++---
 Tests/MeshLib/TestQuadraticMesh.cpp           | 192 +++++++++++-------
 Tests/MeshLib/TestTriLineMesh.cpp             |  18 +-
 30 files changed, 358 insertions(+), 272 deletions(-)

diff --git a/Applications/Utils/MeshEdit/queryMesh.cpp b/Applications/Utils/MeshEdit/queryMesh.cpp
index 1ebd40fde92..5d00953d2ce 100644
--- a/Applications/Utils/MeshEdit/queryMesh.cpp
+++ b/Applications/Utils/MeshEdit/queryMesh.cpp
@@ -63,8 +63,10 @@ int main(int argc, char* argv[])
     {
         auto itr = std::max_element(
             mesh->getNodes().begin(), mesh->getNodes().end(),
-            [](MeshLib::Node* i, MeshLib::Node* j) {
-                return i->getNumberOfElements() < j->getNumberOfElements();
+            [&mesh](MeshLib::Node* i, MeshLib::Node* j)
+            {
+                return mesh->getElementsConnectedToNode(*i).size() <
+                       mesh->getElementsConnectedToNode(*j).size();
             });
         if (itr != mesh->getNodes().end())
         {
@@ -125,8 +127,9 @@ int main(int argc, char* argv[])
         MeshLib::Node const* node = mesh->getNode(node_id);
         out << "# Node " << node->getID() << std::endl;
         out << "Coordinates: " << *node << std::endl;
-        out << "Connected elements (" << node->getNumberOfElements() << "): ";
-        for (auto ele : node->getElements())
+        out << "Connected elements ("
+            << mesh->getElementsConnectedToNode(*node).size() << "): ";
+        for (auto ele : mesh->getElementsConnectedToNode(*node))
         {
             out << ele->getID() << " ";
         }
diff --git a/Applications/Utils/MeshGeoTools/IntegrateBoreholesIntoMesh.cpp b/Applications/Utils/MeshGeoTools/IntegrateBoreholesIntoMesh.cpp
index 460ff05dfb6..83346837d00 100644
--- a/Applications/Utils/MeshGeoTools/IntegrateBoreholesIntoMesh.cpp
+++ b/Applications/Utils/MeshGeoTools/IntegrateBoreholesIntoMesh.cpp
@@ -31,7 +31,8 @@ std::vector<std::size_t> getNodes(
     GeoLib::Point const& pnt, std::vector<MeshLib::Node*> const& nodes,
     MeshLib::PropertyVector<int> const& mat_ids,
     std::pair<int, int> const& mat_limits,
-    std::pair<double, double> const& elevation_limits)
+    std::pair<double, double> const& elevation_limits,
+    MeshLib::Mesh const& mesh)
 {
     std::vector<std::size_t> pnt_nodes;
     for (auto node : nodes)
@@ -40,7 +41,7 @@ std::vector<std::size_t> getNodes(
         if (std::abs((*node)[0] - pnt[0]) < eps &&
             std::abs((*node)[1] - pnt[1]) < eps)
         {
-            auto const& elems = node->getElements();
+            auto const& elems = mesh.getElementsConnectedToNode(*node);
             for (auto e : elems)
             {
                 if (mat_ids[e->getID()] >= mat_limits.first &&
@@ -207,8 +208,8 @@ int main(int argc, char* argv[])
 
     for (std::size_t i = 0; i < n_points; ++i)
     {
-        std::vector<std::size_t> const& line_nodes =
-            getNodes(*points[i], nodes, *mat_ids, mat_limits, elevation_limits);
+        std::vector<std::size_t> const& line_nodes = getNodes(
+            *points[i], nodes, *mat_ids, mat_limits, elevation_limits, *mesh);
         std::size_t const n_line_nodes = line_nodes.size();
         if (n_line_nodes < 2)
         {
diff --git a/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.cpp b/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.cpp
index e5f7f913256..85ac5a67c24 100644
--- a/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.cpp
+++ b/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.cpp
@@ -139,6 +139,7 @@ findRegularNodesInPartition(
     const bool is_mixed_high_order_linear_elems,
     std::vector<MeshLib::Node*> const& nodes,
     std::vector<std::size_t> const& partition_ids,
+    MeshLib::Mesh const& mesh,
     std::vector<std::size_t> const* node_id_mapping = nullptr)
 {
     // Find nodes belonging to a given partition id.
@@ -165,7 +166,10 @@ findRegularNodesInPartition(
         begin(partition_nodes), end(partition_nodes),
         std::back_inserter(base_nodes), std::back_inserter(higher_order_nodes),
         [&](MeshLib::Node* const n)
-        { return !is_mixed_high_order_linear_elems || isBaseNode(*n); });
+        {
+            return !is_mixed_high_order_linear_elems ||
+                   isBaseNode(*n, mesh.getElementsConnectedToNode(*n));
+        });
 
     return {base_nodes, higher_order_nodes};
 }
@@ -233,6 +237,7 @@ findGhostNodesInPartition(
     std::vector<MeshLib::Node*> const& nodes,
     std::vector<MeshLib::Element const*> const& ghost_elements,
     std::vector<std::size_t> const& partition_ids,
+    MeshLib::Mesh const& mesh,
     std::vector<std::size_t> const* node_id_mapping = nullptr)
 {
     std::vector<MeshLib::Node*> base_nodes;
@@ -251,7 +256,8 @@ findGhostNodesInPartition(
 
             if (partitionLookup(*n, partition_ids, node_id_mapping) != part_id)
             {
-                if (!is_mixed_high_order_linear_elems || isBaseNode(*n))
+                if (!is_mixed_high_order_linear_elems ||
+                    isBaseNode(*n, mesh.getElementsConnectedToNode(*n)))
                 {
                     base_nodes.push_back(nodes[n->getID()]);
                 }
@@ -274,7 +280,8 @@ void NodeWiseMeshPartitioner::processPartition(
     std::vector<MeshLib::Node*> higher_order_regular_nodes;
     std::tie(partition.nodes, higher_order_regular_nodes) =
         findRegularNodesInPartition(part_id, is_mixed_high_order_linear_elems,
-                                    _mesh->getNodes(), _nodes_partition_ids);
+                                    _mesh->getNodes(), _nodes_partition_ids,
+                                    *_mesh);
 
     partition.number_of_regular_base_nodes = partition.nodes.size();
     partition.number_of_regular_nodes = partition.number_of_regular_base_nodes +
@@ -288,7 +295,7 @@ void NodeWiseMeshPartitioner::processPartition(
     std::tie(base_ghost_nodes, higher_order_ghost_nodes) =
         findGhostNodesInPartition(part_id, is_mixed_high_order_linear_elems,
                                   _mesh->getNodes(), partition.ghost_elements,
-                                  _nodes_partition_ids);
+                                  _nodes_partition_ids, *_mesh);
 
     std::copy(begin(base_ghost_nodes), end(base_ghost_nodes),
               std::back_inserter(partition.nodes));
@@ -671,7 +678,7 @@ std::vector<Partition> NodeWiseMeshPartitioner::partitionOtherMesh(
         std::tie(partition.nodes, higher_order_regular_nodes) =
             findRegularNodesInPartition(
                 part_id, is_mixed_high_order_linear_elems, mesh.getNodes(),
-                _nodes_partition_ids, bulk_node_ids);
+                _nodes_partition_ids, mesh, bulk_node_ids);
 
         partition.number_of_regular_base_nodes = partition.nodes.size();
         partition.number_of_regular_nodes =
@@ -687,7 +694,8 @@ std::vector<Partition> NodeWiseMeshPartitioner::partitionOtherMesh(
         std::tie(base_ghost_nodes, higher_order_ghost_nodes) =
             findGhostNodesInPartition(part_id, is_mixed_high_order_linear_elems,
                                       mesh.getNodes(), partition.ghost_elements,
-                                      _nodes_partition_ids, bulk_node_ids);
+                                      _nodes_partition_ids, mesh,
+                                      bulk_node_ids);
 
         std::copy(begin(base_ghost_nodes), end(base_ghost_nodes),
                   std::back_inserter(partition.nodes));
diff --git a/Applications/Utils/ModelPreparation/createNeumannBc.cpp b/Applications/Utils/ModelPreparation/createNeumannBc.cpp
index 9e5ab188c8b..d6d1f70c6f1 100644
--- a/Applications/Utils/ModelPreparation/createNeumannBc.cpp
+++ b/Applications/Utils/ModelPreparation/createNeumannBc.cpp
@@ -54,7 +54,8 @@ std::vector<double> getSurfaceIntegratedValuesForNodes(
     {
         double node_area(0);
         double integrated_node_area(0);
-        for (auto const& connected_elem : node->getElements())
+        for (auto const& connected_elem :
+             mesh.getElementsConnectedToNode(*node))
         {
             double const area = connected_elem->getContent() /
                                 connected_elem->getNumberOfBaseNodes();
diff --git a/MeshGeoToolsLib/GeoMapper.cpp b/MeshGeoToolsLib/GeoMapper.cpp
index a43f3adc4b0..4ebe82656f0 100644
--- a/MeshGeoToolsLib/GeoMapper.cpp
+++ b/MeshGeoToolsLib/GeoMapper.cpp
@@ -220,8 +220,8 @@ double GeoMapper::getMeshElevation(double x, double y, double min_val,
 {
     const MeshLib::Node* pnt =
         _grid->getNearestPoint(MathLib::Point3d{{{x, y, 0}}});
-    const std::vector<MeshLib::Element*> elements(
-        _surface_mesh->getNode(pnt->getID())->getElements());
+    auto const elements(
+        _surface_mesh->getElementsConnectedToNode(pnt->getID()));
     std::unique_ptr<GeoLib::Point> intersection;
 
     for (auto const& element : elements)
diff --git a/MeshGeoToolsLib/IdentifySubdomainMesh.cpp b/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
index 62daa8e04d5..0cd1b936a97 100644
--- a/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
+++ b/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
@@ -59,7 +59,8 @@ std::vector<std::size_t> findElementsInMesh(
 
     for (auto const node_id : node_ids)
     {
-        auto const& connected_elements = mesh.getNode(node_id)->getElements();
+        auto const& connected_elements =
+            mesh.getElementsConnectedToNode(node_id);
         std::transform(
             begin(connected_elements), end(connected_elements),
             back_inserter(common_element_ids),
diff --git a/MeshGeoToolsLib/MeshNodesOnPoint.cpp b/MeshGeoToolsLib/MeshNodesOnPoint.cpp
index edeabd3664f..e70c1fa3af0 100644
--- a/MeshGeoToolsLib/MeshNodesOnPoint.cpp
+++ b/MeshGeoToolsLib/MeshNodesOnPoint.cpp
@@ -30,7 +30,8 @@ MeshNodesOnPoint::MeshNodesOnPoint(MeshLib::Mesh const& mesh,
     {
         for (auto id : vec_ids)
         {
-            if (isBaseNode(*mesh.getNode(id)))
+            if (MeshLib::isBaseNode(*_mesh.getNode(id),
+                                    _mesh.getElementsConnectedToNode(id)))
             {
                 _msh_node_ids.push_back(id);
             }
diff --git a/MeshLib/ElementStatus.cpp b/MeshLib/ElementStatus.cpp
index b82094acd83..1f894644833 100644
--- a/MeshLib/ElementStatus.cpp
+++ b/MeshLib/ElementStatus.cpp
@@ -26,9 +26,9 @@ ElementStatus::ElementStatus(Mesh const* const mesh, bool hasAnyInactive)
       _hasAnyInactive(hasAnyInactive)
 {
     const std::vector<MeshLib::Node*>& nodes(_mesh->getNodes());
-    std::transform(
-        begin(nodes), end(nodes), back_inserter(_active_nodes),
-        [](Node const* const n) { return n->getNumberOfElements(); });
+    std::transform(begin(nodes), end(nodes), back_inserter(_active_nodes),
+                   [this](Node const* const n)
+                   { return _mesh->getElementsConnectedToNode(*n).size(); });
 }
 
 ElementStatus::ElementStatus(Mesh const* const mesh,
diff --git a/MeshLib/Mesh.cpp b/MeshLib/Mesh.cpp
index 6db16060a72..8bd7c10bda6 100644
--- a/MeshLib/Mesh.cpp
+++ b/MeshLib/Mesh.cpp
@@ -32,6 +32,25 @@ static std::size_t global_mesh_counter = 0;
 
 namespace MeshLib
 {
+std::vector<std::vector<Element const*>> findElementsConnectedToNodes(
+    Mesh const& mesh)
+{
+    std::vector<std::vector<Element const*>> elements_connected_to_nodes;
+    auto const& nodes = mesh.getNodes();
+    elements_connected_to_nodes.resize(nodes.size());
+
+    for (auto const* element : mesh.getElements())
+    {
+        unsigned const number_nodes(element->getNumberOfNodes());
+        for (unsigned j = 0; j < number_nodes; ++j)
+        {
+            auto const node_id = element->getNode(j)->getID();
+            elements_connected_to_nodes[node_id].push_back(element);
+        }
+    }
+    return elements_connected_to_nodes;
+}
+
 Mesh::Mesh(std::string name,
            std::vector<Node*>
                nodes,
@@ -50,9 +69,10 @@ Mesh::Mesh(std::string name,
     this->resetNodeIDs();
     this->resetElementIDs();
     this->setDimension();
-    this->setElementsConnectedToNodes();
-    this->setElementNeighbors();
 
+    _elements_connected_to_nodes = findElementsConnectedToNodes(*this);
+
+    this->setElementNeighbors();
     this->calcEdgeLengthRange();
 }
 
@@ -89,7 +109,7 @@ Mesh::Mesh(const Mesh& mesh)
     {
         this->setDimension();
     }
-    this->setElementsConnectedToNodes();
+    _elements_connected_to_nodes = findElementsConnectedToNodes(*this);
     this->setElementNeighbors();
 }
 
@@ -111,13 +131,6 @@ Mesh::~Mesh()
 void Mesh::addElement(Element* elem)
 {
     _elements.push_back(elem);
-
-    // add element information to nodes
-    unsigned nNodes(elem->getNumberOfNodes());
-    for (unsigned i = 0; i < nNodes; ++i)
-    {
-        elem->getNode(i)->addElement(elem);
-    }
 }
 
 void Mesh::resetNodeIDs()
@@ -150,18 +163,6 @@ void Mesh::setDimension()
     }
 }
 
-void Mesh::setElementsConnectedToNodes()
-{
-    for (auto& element : _elements)
-    {
-        const unsigned nNodes(element->getNumberOfNodes());
-        for (unsigned j = 0; j < nNodes; ++j)
-        {
-            element->getNode(j)->addElement(element);
-        }
-    }
-}
-
 void Mesh::calcEdgeLengthRange()
 {
     const std::size_t nElems(getNumberOfElements());
@@ -178,7 +179,7 @@ void Mesh::calcEdgeLengthRange()
 
 void Mesh::setElementNeighbors()
 {
-    std::vector<Element*> neighbors;
+    std::vector<Element const*> neighbors;
     for (auto element : _elements)
     {
         // create vector with all elements connected to current element
@@ -186,8 +187,8 @@ void Mesh::setElementNeighbors()
         const std::size_t nNodes(element->getNumberOfBaseNodes());
         for (unsigned n(0); n < nNodes; ++n)
         {
-            std::vector<Element*> const& conn_elems(
-                (element->getNode(n)->getElements()));
+            auto const& conn_elems(
+                _elements_connected_to_nodes[element->getNode(n)->getID()]);
             neighbors.insert(neighbors.end(), conn_elems.begin(),
                              conn_elems.end());
         }
@@ -199,10 +200,11 @@ void Mesh::setElementNeighbors()
              ++neighbor)
         {
             std::optional<unsigned> const opposite_face_id =
-                element->addNeighbor(*neighbor);
+                element->addNeighbor(const_cast<Element*>(*neighbor));
             if (opposite_face_id)
             {
-                (*neighbor)->setNeighbor(element, *opposite_face_id);
+                const_cast<Element*>(*neighbor)->setNeighbor(element,
+                                                             *opposite_face_id);
             }
         }
         neighbors.clear();
@@ -212,16 +214,31 @@ void Mesh::setElementNeighbors()
 std::size_t Mesh::getNumberOfBaseNodes() const
 {
     return std::count_if(begin(_nodes), end(_nodes),
-                         [](auto const* const node)
-                         { return isBaseNode(*node); });
+                         [this](auto const* const node) {
+                             return isBaseNode(
+                                 *node,
+                                 _elements_connected_to_nodes[node->getID()]);
+                         });
 }
 
 bool Mesh::hasNonlinearElement() const
 {
     return std::any_of(
-        std::begin(_elements), std::end(_elements), [](Element const* const e) {
-            return e->getNumberOfNodes() != e->getNumberOfBaseNodes();
-        });
+        std::begin(_elements), std::end(_elements),
+        [](Element const* const e)
+        { return e->getNumberOfNodes() != e->getNumberOfBaseNodes(); });
+}
+
+std::vector<MeshLib::Element const*> const& Mesh::getElementsConnectedToNode(
+    std::size_t const node_id) const
+{
+    return _elements_connected_to_nodes[node_id];
+}
+
+std::vector<MeshLib::Element const*> const& Mesh::getElementsConnectedToNode(
+    Node const& node) const
+{
+    return _elements_connected_to_nodes[node.getID()];
 }
 
 void scaleMeshPropertyVector(MeshLib::Mesh& mesh,
@@ -314,6 +331,8 @@ std::unique_ptr<MeshLib::Mesh> createMeshFromElementSelection(
 std::vector<std::vector<Node*>> calculateNodesConnectedByElements(
     Mesh const& mesh)
 {
+    auto const elements_connected_to_nodes = findElementsConnectedToNodes(mesh);
+
     std::vector<std::vector<Node*>> nodes_connected_by_elements;
     auto const& nodes = mesh.getNodes();
     nodes_connected_by_elements.resize(nodes.size());
@@ -323,7 +342,8 @@ std::vector<std::vector<Node*>> calculateNodesConnectedByElements(
         auto const* node = nodes[i];
 
         // Get all elements, to which this node is connected.
-        auto const& connected_elements = node->getElements();
+        auto const& connected_elements =
+            elements_connected_to_nodes[node->getID()];
 
         // And collect all elements' nodes.
         for (Element const* const element : connected_elements)
@@ -347,4 +367,20 @@ std::vector<std::vector<Node*>> calculateNodesConnectedByElements(
     return nodes_connected_by_elements;
 }
 
+bool isBaseNode(Node const& node,
+                std::vector<Element const*> const& elements_connected_to_node)
+{
+    // Check if node is connected.
+    if (elements_connected_to_node.empty())
+    {
+        return true;
+    }
+
+    // In a mesh a node always belongs to at least one element.
+    auto const e = elements_connected_to_node[0];
+
+    auto const n_base_nodes = e->getNumberOfBaseNodes();
+    auto const local_index = getNodeIDinElement(*e, &node);
+    return local_index < n_base_nodes;
+}
 }  // namespace MeshLib
diff --git a/MeshLib/Mesh.h b/MeshLib/Mesh.h
index 2d384ac2185..7ef88772a6f 100644
--- a/MeshLib/Mesh.h
+++ b/MeshLib/Mesh.h
@@ -115,6 +115,11 @@ public:
     /// Check if the mesh contains any nonlinear element.
     bool hasNonlinearElement() const;
 
+    std::vector<Element const*> const& getElementsConnectedToNode(
+        std::size_t node_id) const;
+    std::vector<Element const*> const& getElementsConnectedToNode(
+        Node const& node) const;
+
     Properties& getProperties() { return _properties; }
     Properties const& getProperties() const { return _properties; }
 
@@ -130,9 +135,6 @@ protected:
     /// Sets the dimension of the mesh.
     void setDimension();
 
-    /// Fills in the neighbor-information for nodes (i.e. which element each node belongs to).
-    void setElementsConnectedToNodes();
-
     /// Fills in the neighbor-information for elements.
     /// Note: Using this implementation, an element e can only have neighbors that have the same dimensionality as e.
     void setElementNeighbors();
@@ -150,6 +152,8 @@ protected:
     std::vector<Element*> _elements;
     Properties _properties;
 
+    std::vector<std::vector<Element const*>> _elements_connected_to_nodes;
+
     bool _is_axially_symmetric = false;
 }; /* class */
 
@@ -292,4 +296,8 @@ PropertyVector<int> const* materialIDs(Mesh const& mesh);
 std::unique_ptr<Mesh> createMeshFromElementSelection(
     std::string mesh_name, std::vector<Element*> const& elements);
 
+/// Returns true if the given node is a base node of a (first) element, or if it
+/// is not connected to any element i.e. an unconnected node.
+bool isBaseNode(Node const& node,
+                std::vector<Element const*> const& elements_connected_to_node);
 }  // namespace MeshLib
diff --git a/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp b/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp
index 3504a8efec1..b94fd7792b4 100644
--- a/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp
+++ b/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp
@@ -165,13 +165,16 @@ void Mesh2MeshPropertyInterpolation::
     const std::size_t n_src_nodes(src_nodes.size());
     for (std::size_t k(0); k < n_src_nodes; k++)
     {
-        const std::size_t n_con_elems(src_nodes[k]->getNumberOfElements());
-        interpolated_properties[k] =
-            (*elem_props)[(src_nodes[k]->getElement(0))->getID()];
+        const std::size_t n_con_elems(
+            _src_mesh.getElementsConnectedToNode(*src_nodes[k]).size());
+        interpolated_properties[k] = (*elem_props)
+            [_src_mesh.getElementsConnectedToNode(*src_nodes[k])[0]->getID()];
         for (std::size_t j(1); j < n_con_elems; j++)
         {
             interpolated_properties[k] +=
-                (*elem_props)[(src_nodes[k]->getElement(j))->getID()];
+                (*elem_props)[_src_mesh
+                                  .getElementsConnectedToNode(*src_nodes[k])[j]
+                                  ->getID()];
         }
         interpolated_properties[k] /= n_con_elems;
     }
diff --git a/MeshLib/MeshSearch/ElementSearch.cpp b/MeshLib/MeshSearch/ElementSearch.cpp
index 7e5c8874790..513a69f7eee 100644
--- a/MeshLib/MeshSearch/ElementSearch.cpp
+++ b/MeshLib/MeshSearch/ElementSearch.cpp
@@ -78,7 +78,7 @@ std::size_t ElementSearch::searchByNodeIDs(
     std::vector<std::size_t> connected_elements;
     for (std::size_t node_id : nodes)
     {
-        auto const& elements = _mesh.getNode(node_id)->getElements();
+        auto const& elements = _mesh.getElementsConnectedToNode(node_id);
         std::transform(begin(elements), end(elements),
                        back_inserter(connected_elements),
                        [](Element const* const e) { return e->getID(); });
diff --git a/MeshLib/MeshSearch/NodeSearch.cpp b/MeshLib/MeshSearch/NodeSearch.cpp
index 7f3a4646154..866071dad6f 100644
--- a/MeshLib/MeshSearch/NodeSearch.cpp
+++ b/MeshLib/MeshSearch/NodeSearch.cpp
@@ -44,7 +44,8 @@ std::size_t NodeSearch::searchNodesConnectedToOnlyGivenElements(
     std::vector<std::size_t> connected_nodes;
     for (std::size_t i = 0; i < node_marked_counts.size(); i++)
     {
-        if (node_marked_counts[i] == _mesh.getNode(i)->getElements().size())
+        if (node_marked_counts[i] ==
+            _mesh.getElementsConnectedToNode(*_mesh.getNode(i)).size())
         {
             connected_nodes.push_back(i);
         }
@@ -62,7 +63,7 @@ std::size_t NodeSearch::searchUnused()
 
     for (unsigned i = 0; i < nNodes; ++i)
     {
-        if (nodes[i]->getNumberOfElements() == 0)
+        if (_mesh.getElementsConnectedToNode(*nodes[i]).empty())
         {
             del_node_idx.push_back(i);
         }
@@ -79,7 +80,7 @@ std::size_t NodeSearch::searchBoundaryNodes()
     {
         for (MeshLib::Node const* n : _mesh.getNodes())
         {
-            if (n->getElements().size() == 1)
+            if (_mesh.getElementsConnectedToNode(*n).size() == 1)
             {
                 vec_boundary_nodes.push_back(n->getID());
             }
diff --git a/MeshLib/MeshSurfaceExtraction.cpp b/MeshLib/MeshSurfaceExtraction.cpp
index 72eaff38961..d2a2871ebe6 100644
--- a/MeshLib/MeshSurfaceExtraction.cpp
+++ b/MeshLib/MeshSurfaceExtraction.cpp
@@ -205,7 +205,7 @@ std::vector<double> MeshSurfaceExtraction::getSurfaceAreaForNodes(
     {
         double node_area(0);
 
-        std::vector<MeshLib::Element*> conn_elems = nodes[n]->getElements();
+        auto const conn_elems = mesh.getElementsConnectedToNode(*nodes[n]);
         const std::size_t nConnElems(conn_elems.size());
 
         for (std::size_t i = 0; i < nConnElems; ++i)
diff --git a/MeshLib/Node.cpp b/MeshLib/Node.cpp
index 54d11e09de0..fa00d6ff279 100644
--- a/MeshLib/Node.cpp
+++ b/MeshLib/Node.cpp
@@ -41,27 +41,5 @@ void Node::updateCoordinates(double x, double y, double z)
     (*this)[0] = x;
     (*this)[1] = y;
     (*this)[2] = z;
-
-    const std::size_t nElements(this->_elements.size());
-    for (std::size_t i = 0; i < nElements; i++)
-    {
-        _elements[i]->computeVolume();
-    }
-}
-
-bool isBaseNode(Node const& node)
-{
-    // Check if node is connected.
-    if (node.getNumberOfElements() == 0)
-    {
-        return true;
-    }
-
-    // In a mesh a node always belongs to at least one element.
-    auto const e = node.getElement(0);
-
-    auto const n_base_nodes = e->getNumberOfBaseNodes();
-    auto const local_index = getNodeIDinElement(*e, &node);
-    return local_index < n_base_nodes;
 }
 }  // namespace MeshLib
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index cc250407242..e3840d8e918 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -25,10 +25,8 @@ namespace ApplicationUtils
     class NodeWiseMeshPartitioner;
 }
 
-namespace MeshLib {
-
-class Element;
-
+namespace MeshLib
+{
 /**
  * A mesh node with coordinates in 3D space.
  */
@@ -56,34 +54,10 @@ public:
     /// Copy constructor
     Node(const Node &node);
 
-    /// Get an element the node is part of.
-    const Element* getElement(std::size_t idx) const { return _elements[idx]; }
-
-    /// Get all elements the node is part of.
-    const std::vector<Element*>& getElements() const { return _elements; }
-
-    /// Get number of elements the node is part of.
-    std::size_t getNumberOfElements() const { return _elements.size(); }
-
 protected:
     /// Update coordinates of a node.
     /// This method automatically also updates the areas/volumes of all connected elements.
     void updateCoordinates(double x, double y, double z);
 
-    /**
-     * Add an element the node is part of.
-     * This method is called by Mesh::addElement(Element*), see friend definition.
-     */
-    void addElement(Element* elem) { _elements.push_back(elem); }
-
-    /// clear stored elements connecting to this node
-    void clearElements() { _elements.clear(); }
-
-    std::vector<Element*> _elements;
 }; /* class */
-
-/// Returns true if the given node is a base node of a (first) element, or if it
-/// is not connected to any element i.e. an unconnected node.
-bool isBaseNode(Node const& node);
-
 }  // namespace MeshLib
diff --git a/MeshLib/NodePartitionedMesh.h b/MeshLib/NodePartitionedMesh.h
index 73fc94c3ca8..8080bd0a7b0 100644
--- a/MeshLib/NodePartitionedMesh.h
+++ b/MeshLib/NodePartitionedMesh.h
@@ -110,7 +110,9 @@ public:
         {
             return false;
         }
-        if (!isBaseNode(*_nodes[node_id]) && node_id < getLargestActiveNodeID())
+        if (!isBaseNode(*_nodes[node_id],
+                        getElementsConnectedToNode(node_id)) &&
+            node_id < getLargestActiveNodeID())
         {
             return false;
         }
diff --git a/NumLib/DOF/MeshComponentMap.cpp b/NumLib/DOF/MeshComponentMap.cpp
index 14cc5cb0312..880435ca1a4 100644
--- a/NumLib/DOF/MeshComponentMap.cpp
+++ b/NumLib/DOF/MeshComponentMap.cpp
@@ -96,7 +96,9 @@ MeshComponentMap MeshComponentMap::getSubset(
     for (auto* const node : new_mesh_subset.getNodes())
     {
         auto const node_id = node->getID();
-        bool const is_base_node = isBaseNode(*node);
+        bool const is_base_node = MeshLib::isBaseNode(
+            *node,
+            new_mesh_subset.getMesh().getElementsConnectedToNode(node_id));
 
         MeshLib::Location const new_location{
             new_mesh_id, MeshLib::MeshItemType::Node, node_id};
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/DeactivatedSubdomainDirichlet.cpp b/ProcessLib/BoundaryConditionAndSourceTerm/DeactivatedSubdomainDirichlet.cpp
index d75678fd330..f6d9a30cda7 100644
--- a/ProcessLib/BoundaryConditionAndSourceTerm/DeactivatedSubdomainDirichlet.cpp
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/DeactivatedSubdomainDirichlet.cpp
@@ -72,7 +72,8 @@ void DeactivatedSubdomainDirichlet::getEssentialBCValues(
                  back_inserter(inactive_nodes_in_bc_mesh),
                  [&](MeshLib::Node* const n)
                  {
-                     const auto& connected_elements = n->getElements();
+                     const auto& connected_elements =
+                         _subdomain.mesh->getElementsConnectedToNode(*n);
 
                      return std::all_of(begin(connected_elements),
                                         end(connected_elements), is_inactive);
@@ -84,7 +85,8 @@ void DeactivatedSubdomainDirichlet::getEssentialBCValues(
                      back_inserter(inactive_nodes_in_bc_mesh),
                      [&](MeshLib::Node* const n)
                      {
-                         const auto& connected_elements = n->getElements();
+                         const auto& connected_elements =
+                             _subdomain.mesh->getElementsConnectedToNode(*n);
 
                          return std::all_of(begin(connected_elements),
                                             end(connected_elements),
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/DeactivatedSubdomainDirichlet.h b/ProcessLib/BoundaryConditionAndSourceTerm/DeactivatedSubdomainDirichlet.h
index 6b06c6aac66..73c4f920d17 100644
--- a/ProcessLib/BoundaryConditionAndSourceTerm/DeactivatedSubdomainDirichlet.h
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/DeactivatedSubdomainDirichlet.h
@@ -33,10 +33,7 @@ struct Parameter;
 namespace ProcessLib
 {
 struct DeactivatedSubdomainMesh;
-}
 
-namespace ProcessLib
-{
 class DeactivatedSubdomainDirichlet final : public BoundaryCondition
 {
 public:
diff --git a/ProcessLib/CreateDeactivatedSubdomain.cpp b/ProcessLib/CreateDeactivatedSubdomain.cpp
index 602fbbb863f..711999fc702 100644
--- a/ProcessLib/CreateDeactivatedSubdomain.cpp
+++ b/ProcessLib/CreateDeactivatedSubdomain.cpp
@@ -50,9 +50,8 @@ extractInnerAndOuterNodes(MeshLib::Mesh const& mesh,
         back_inserter(inner_nodes), back_inserter(outer_nodes),
         [&](MeshLib::Node* const n) {
             auto const bulk_node = mesh.getNode((*bulk_node_ids)[n->getID()]);
-            const auto& connected_elements = bulk_node->getElements();
-
-            // Check whether this node is connected only to an active
+            const auto& connected_elements =
+                mesh.getElementsConnectedToNode(*bulk_node);
             // elements. Then it is an inner node, and outer otherwise.
             return std::all_of(begin(connected_elements),
                                end(connected_elements), is_active);
diff --git a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
index d9738892f9e..918c3e88b8c 100644
--- a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
+++ b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
@@ -327,8 +327,10 @@ void HeatTransportBHEProcess::createBHEBoundaryConditionTopBottom(
         for (auto const& bhe_node : bhe_nodes)
         {
             // Count number of 1d elements connected with every BHE node.
+            auto const& connected_elements =
+                _mesh.getElementsConnectedToNode(*bhe_node);
             const std::size_t n_line_elements = std::count_if(
-                bhe_node->getElements().begin(), bhe_node->getElements().end(),
+                connected_elements.begin(), connected_elements.end(),
                 [](MeshLib::Element const* elem)
                 { return (elem->getDimension() == 1); });
 
diff --git a/ProcessLib/LIE/Common/MeshUtils.cpp b/ProcessLib/LIE/Common/MeshUtils.cpp
index 8eae30a0dbe..328f57613e9 100644
--- a/ProcessLib/LIE/Common/MeshUtils.cpp
+++ b/ProcessLib/LIE/Common/MeshUtils.cpp
@@ -29,9 +29,9 @@ class IsCrackTip
 {
 public:
     explicit IsCrackTip(MeshLib::Mesh const& mesh)
-        : _fracture_element_dim(mesh.getDimension() - 1)
+        : _mesh(mesh), _fracture_element_dim(_mesh.getDimension() - 1)
     {
-        _is_internal_node.resize(mesh.getNumberOfNodes(), true);
+        _is_internal_node.resize(_mesh.getNumberOfNodes(), true);
 
         MeshLib::NodeSearch nodeSearch(mesh);
         nodeSearch.searchBoundaryNodes();
@@ -43,29 +43,32 @@ public:
 
     bool operator()(MeshLib::Node const& node) const
     {
-        if (!_is_internal_node[node.getID()] || !isBaseNode(node))
+        if (!_is_internal_node[node.getID()] ||
+            !isBaseNode(node, _mesh.getElementsConnectedToNode(node)))
         {
             return false;
         }
 
+        auto const elements_connected_to_node =
+            _mesh.getElementsConnectedToNode(node);
         auto const n_connected_fracture_elements =
-            count_if(cbegin(node.getElements()), cend(node.getElements()),
-                     [&](auto* const e) {
-                         return e->getDimension() == _fracture_element_dim;
-                     });
+            count_if(cbegin(elements_connected_to_node),
+                     cend(elements_connected_to_node),
+                     [&](auto* const e)
+                     { return e->getDimension() == _fracture_element_dim; });
         assert(n_connected_fracture_elements > 0);
 
         return (n_connected_fracture_elements == 1);
     }
 
 private:
+    MeshLib::Mesh const& _mesh;
     unsigned const _fracture_element_dim;
     std::vector<bool> _is_internal_node;
 };
 
 void findFracutreIntersections(
-    MeshLib::Mesh const& mesh,
-    std::vector<int> const& vec_fracture_mat_IDs,
+    MeshLib::Mesh const& mesh, std::vector<int> const& vec_fracture_mat_IDs,
     std::vector<std::vector<MeshLib::Node*>> const& vec_fracture_nodes,
     std::vector<std::vector<MeshLib::Element*>>& intersected_fracture_elements,
     std::vector<std::pair<std::size_t, std::vector<int>>>&
@@ -113,9 +116,11 @@ void findFracutreIntersections(
             continue;  // no intersection
         }
 
+        auto const elements_connected_to_node =
+            mesh.getElementsConnectedToNode(*node);
         std::vector<MeshLib::Element*> conn_fracture_elements;
         {
-            for (auto const* e : node->getElements())
+            for (auto const* e : elements_connected_to_node)
             {
                 if (e->getDimension() == (mesh.getDimension() - 1))
                 {
@@ -315,16 +320,18 @@ void getFractureMatrixDataInMesh(
                 {
                     continue;
                 }
-                for (unsigned j = 0; j < node->getNumberOfElements(); j++)
+                auto const& elements_connected_to_node =
+                    mesh.getElementsConnectedToNode(*node);
+                for (unsigned j = 0; j < elements_connected_to_node.size(); j++)
                 {
                     // only matrix elements
-                    if (node->getElement(j)->getDimension() <
+                    if (elements_connected_to_node[j]->getDimension() <
                         mesh.getDimension())
                     {
                         continue;
                     }
-                    vec_ele.push_back(
-                        const_cast<MeshLib::Element*>(node->getElement(j)));
+                    vec_ele.push_back(const_cast<MeshLib::Element*>(
+                        elements_connected_to_node[j]));
                 }
             }
         }
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
index a6186f4f643..cc99fb75355 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
@@ -144,7 +144,7 @@ SmallDeformationProcess<DisplacementDim>::SmallDeformationProcess(
     for (unsigned i = 0; i < vec_junction_nodeID_matIDs.size(); i++)
     {
         auto node = mesh.getNode(vec_junction_nodeID_matIDs[i].first);
-        for (auto e : node->getElements())
+        for (auto e : mesh.getElementsConnectedToNode(*node))
         {
             _process_data._vec_ele_connected_junctionIDs[e->getID()].push_back(
                 i);
@@ -157,9 +157,10 @@ SmallDeformationProcess<DisplacementDim>::SmallDeformationProcess(
     for (unsigned i = 0; i < vec_junction_nodeID_matIDs.size(); i++)
     {
         auto node = mesh.getNode(vec_junction_nodeID_matIDs[i].first);
-        for (auto e : node->getElements())
+        for (auto e : mesh.getElementsConnectedToNode(*node))
         {
-            _vec_junction_fracture_matrix_elements[i].push_back(e);
+            _vec_junction_fracture_matrix_elements[i].push_back(
+                const_cast<MeshLib::Element*>(e));
         }
     }
 
diff --git a/Tests/MeshLib/TestFindElementsWithinRadius.cpp b/Tests/MeshLib/TestFindElementsWithinRadius.cpp
index 25accad7efa..db8c6b838ed 100644
--- a/Tests/MeshLib/TestFindElementsWithinRadius.cpp
+++ b/Tests/MeshLib/TestFindElementsWithinRadius.cpp
@@ -44,10 +44,11 @@ TEST_F(MeshLibFindElementWithinRadius, ZeroRadius)
                            ac::make_arbitrary<std::size_t>(), gtest_reporter);
 }
 
-std::vector<std::size_t> findElementIdsConnectedToNode(Node const& node)
+std::vector<std::size_t> findElementIdsConnectedToNode(
+    Node const& node, MeshLib::Mesh const& mesh)
 {
     std::vector<std::size_t> connected_elements;
-    auto const& elements_of_node = node.getElements();
+    auto const& elements_of_node = mesh.getElementsConnectedToNode(node);
     std::transform(std::begin(elements_of_node),
                    std::end(elements_of_node),
                    std::back_inserter(connected_elements),
@@ -57,13 +58,13 @@ std::vector<std::size_t> findElementIdsConnectedToNode(Node const& node)
 
 // Find all elements connected through all nodes of the element
 std::vector<std::size_t> findNextNeighborsConnectedByNode(
-    Element const& element)
+    Element const& element, MeshLib::Mesh const& mesh)
 {
     std::vector<std::size_t> connected_elements;
     for (unsigned n = 0; n < element.getNumberOfBaseNodes(); ++n)
     {
         auto const& node = *element.getNode(n);
-        auto const& elements = findElementIdsConnectedToNode(node);
+        auto const& elements = findElementIdsConnectedToNode(node, mesh);
         std::copy(std::begin(elements), std::end(elements),
                   std::back_inserter(connected_elements));
     }
@@ -120,7 +121,7 @@ std::vector<std::size_t> bruteForceFindElementIdsWithinRadius(
     auto const nodes_in_radius = findNodesWithinRadius(mesh, element, radius);
     for (auto n : nodes_in_radius)
     {
-        auto const& elements = findElementIdsConnectedToNode(*n);
+        auto const& elements = findElementIdsConnectedToNode(*n, mesh);
         std::copy(std::begin(elements), std::end(elements),
                   std::back_inserter(connected_elements));
     }
@@ -137,13 +138,14 @@ TEST_F(MeshLibFindElementWithinRadius, VerySmallRadius)
         std::unique_ptr<Mesh>(MeshGenerator::generateRegularQuadMesh(10., 10));
 
     auto neighboring_elements_returned =
-        [&mesh](std::size_t& element_id) -> bool {
+        [&mesh](std::size_t& element_id) -> bool
+    {
         auto const& element = *mesh->getElement(element_id);
         auto result = findElementsWithinRadius(element, 1e-5);
 
         std::sort(std::begin(result), std::end(result));
         auto const expected_elements =
-            findNextNeighborsConnectedByNode(element);
+            findNextNeighborsConnectedByNode(element, *mesh);
 
         return result.size() == expected_elements.size() &&
                std::includes(std::begin(result), std::end(result),
diff --git a/Tests/MeshLib/TestLineMesh.cpp b/Tests/MeshLib/TestLineMesh.cpp
index a9f43ddcb03..eafb0b122d2 100644
--- a/Tests/MeshLib/TestLineMesh.cpp
+++ b/Tests/MeshLib/TestLineMesh.cpp
@@ -104,16 +104,18 @@ TEST_F(MeshLibLineMesh, NodeToElementConnectivity)
 
     // Each node n of the line mesh is connected to two elements n-1 and n,
     // except the first and last nodes.
-    ASSERT_EQ(1u, nodes.front()->getNumberOfElements());
-    ASSERT_EQ(elements[0], nodes.front()->getElement(0));
+    ASSERT_EQ(1u, mesh->getElementsConnectedToNode(*nodes.front()).size());
+    ASSERT_EQ(elements[0], mesh->getElementsConnectedToNode(*nodes.front())[0]);
 
-    ASSERT_EQ(1u, nodes.back()->getNumberOfElements());
-    ASSERT_EQ(elements.back(), nodes.back()->getElement(0));
+    ASSERT_EQ(1u, mesh->getElementsConnectedToNode(*nodes.back()).size());
+    ASSERT_EQ(elements.back(),
+              mesh->getElementsConnectedToNode(*nodes.back())[0]);
 
     for (std::size_t i = 1; i < nodes.size() - 1; ++i)
     {
-        ASSERT_EQ(2u, nodes[i]->getNumberOfElements());
-        ASSERT_EQ(elements[i - 1], nodes[i]->getElement(0));
-        ASSERT_EQ(elements[i], nodes[i]->getElement(1));
+        ASSERT_EQ(2u, mesh->getElementsConnectedToNode(*nodes[i]).size());
+        ASSERT_EQ(elements[i - 1],
+                  mesh->getElementsConnectedToNode(*nodes[i])[0]);
+        ASSERT_EQ(elements[i], mesh->getElementsConnectedToNode(*nodes[i])[1]);
     }
 }
diff --git a/Tests/MeshLib/TestNodeAdjacencyTable.cpp b/Tests/MeshLib/TestNodeAdjacencyTable.cpp
index 984428ee634..b5488462ff0 100644
--- a/Tests/MeshLib/TestNodeAdjacencyTable.cpp
+++ b/Tests/MeshLib/TestNodeAdjacencyTable.cpp
@@ -36,7 +36,8 @@ TEST(MeshLib, CreateNodeAdjacencyTable1D)
     {
         std::size_t const n_connections = table.getNodeDegree(i);
 
-        std::size_t const n_elements = mesh->getNode(i)->getNumberOfElements();
+        std::size_t const n_elements =
+            mesh->getElementsConnectedToNode(i).size();
         switch (n_elements)
         {
             case 1:  // a corner node has 2 adjacent nodes.
@@ -70,7 +71,8 @@ TEST(MeshLib, CreateNodeAdjacencyTable2D)
     {
         std::size_t const n_connections = table.getNodeDegree(i);
 
-        std::size_t const n_elements = mesh->getNode(i)->getNumberOfElements();
+        std::size_t const n_elements =
+            mesh->getElementsConnectedToNode(i).size();
         switch (n_elements)
         {
             case 1:  // a corner node has 4 adjacent nodes.
@@ -109,7 +111,8 @@ TEST(MeshLib, CreateNodeAdjacencyTable3D)
     {
         std::size_t const n_connections = table.getNodeDegree(i);
 
-        std::size_t const n_elements = mesh->getNode(i)->getNumberOfElements();
+        std::size_t const n_elements =
+            mesh->getElementsConnectedToNode(i).size();
         switch (n_elements)
         {
             case 1:  // a corner node has 8 adjacent nodes.
diff --git a/Tests/MeshLib/TestQuadMesh.cpp b/Tests/MeshLib/TestQuadMesh.cpp
index 6fd6861c6fc..64a98ada4ec 100644
--- a/Tests/MeshLib/TestQuadMesh.cpp
+++ b/Tests/MeshLib/TestQuadMesh.cpp
@@ -286,8 +286,9 @@ TEST_F(MeshLibQuadMesh, ElementToNodeConnectivity)
 TEST_F(MeshLibQuadMesh, NodeToElementConnectivity)
 {
     testCornerNodes(
-        [this](MeshLib::Node const* const node, std::size_t i, std::size_t j) {
-            EXPECT_EQ(1u, node->getNumberOfElements());
+        [this](MeshLib::Node const* const node, std::size_t i, std::size_t j)
+        {
+            EXPECT_EQ(1u, mesh->getElementsConnectedToNode(*node).size());
 
             if (i == nodes_stride)
             {
@@ -298,45 +299,60 @@ TEST_F(MeshLibQuadMesh, NodeToElementConnectivity)
                 j--;
             }
 
-            EXPECT_EQ(getElement(i, j), node->getElement(0));
+            EXPECT_EQ(getElement(i, j),
+                      mesh->getElementsConnectedToNode(*node)[0]);
         });
 
     testBoundaryNodes(
-        [this](MeshLib::Node const* const node, std::size_t i, std::size_t j) {
-            EXPECT_EQ(2u, node->getNumberOfElements());
+        [this](MeshLib::Node const* const node, std::size_t i, std::size_t j)
+        {
+            EXPECT_EQ(2u, mesh->getElementsConnectedToNode(*node).size());
 
             if (i == 0)
             {
-                EXPECT_EQ(getElement(i, j - 1), node->getElement(0));
-                EXPECT_EQ(getElement(i, j), node->getElement(1));
+                EXPECT_EQ(getElement(i, j - 1),
+                          mesh->getElementsConnectedToNode(*node)[0]);
+                EXPECT_EQ(getElement(i, j),
+                          mesh->getElementsConnectedToNode(*node)[1]);
             }
             if (i == nodes_stride)
             {
                 EXPECT_EQ(getElement(elements_stride, j - 1),
-                          node->getElement(0));
-                EXPECT_EQ(getElement(elements_stride, j), node->getElement(1));
+                          mesh->getElementsConnectedToNode(*node)[0]);
+                EXPECT_EQ(getElement(elements_stride, j),
+                          mesh->getElementsConnectedToNode(*node)[1]);
             }
             if (j == 0)
             {
-                EXPECT_EQ(getElement(i - 1, j), node->getElement(0));
-                EXPECT_EQ(getElement(i, j), node->getElement(1));
+                EXPECT_EQ(getElement(i - 1, j),
+                          mesh->getElementsConnectedToNode(*node)[0]);
+                EXPECT_EQ(getElement(i, j),
+                          mesh->getElementsConnectedToNode(*node)[1]);
             }
             if (j == nodes_stride)
             {
                 j--;
-                EXPECT_EQ(getElement(i - 1, j), node->getElement(0));
-                EXPECT_EQ(getElement(i, j), node->getElement(1));
+                EXPECT_EQ(getElement(i - 1, j),
+                          mesh->getElementsConnectedToNode(*node)[0]);
+                EXPECT_EQ(getElement(i, j),
+                          mesh->getElementsConnectedToNode(*node)[1]);
             }
         });
 
-    testInsideNodes([this](MeshLib::Node const* const node,
-                           std::size_t const i,
-                           std::size_t const j) {
-        EXPECT_EQ(4u, node->getNumberOfElements());
-
-        EXPECT_EQ(getElement(i - 1, j - 1), node->getElement(0));
-        EXPECT_EQ(getElement(i - 1, j), node->getElement(1));
-        EXPECT_EQ(getElement(i, j - 1), node->getElement(2));
-        EXPECT_EQ(getElement(i, j), node->getElement(3));
-    });
+    testInsideNodes(
+        [this](MeshLib::Node const* const node,
+               std::size_t const i,
+               std::size_t const j)
+        {
+            EXPECT_EQ(4u, mesh->getElementsConnectedToNode(*node).size());
+
+            EXPECT_EQ(getElement(i - 1, j - 1),
+                      mesh->getElementsConnectedToNode(*node)[0]);
+            EXPECT_EQ(getElement(i - 1, j),
+                      mesh->getElementsConnectedToNode(*node)[1]);
+            EXPECT_EQ(getElement(i, j - 1),
+                      mesh->getElementsConnectedToNode(*node)[2]);
+            EXPECT_EQ(getElement(i, j),
+                      mesh->getElementsConnectedToNode(*node)[3]);
+        });
 }
diff --git a/Tests/MeshLib/TestQuadraticMesh.cpp b/Tests/MeshLib/TestQuadraticMesh.cpp
index 68b525b4aec..7895e84d137 100644
--- a/Tests/MeshLib/TestQuadraticMesh.cpp
+++ b/Tests/MeshLib/TestQuadraticMesh.cpp
@@ -40,12 +40,16 @@ TEST(MeshLib, QuadraticOrderMesh_Line)
 
         for (unsigned i = 0; i < e->getNumberOfBaseNodes(); i++)
         {
-            ASSERT_TRUE(isBaseNode(*e->getNode(i)));
+            ASSERT_TRUE(
+                isBaseNode(*e->getNode(i),
+                           mesh->getElementsConnectedToNode(*e->getNode(i))));
         }
         for (unsigned i = e->getNumberOfBaseNodes(); i < e->getNumberOfNodes();
              i++)
         {
-            ASSERT_FALSE(isBaseNode(*e->getNode(i)));
+            ASSERT_FALSE(
+                isBaseNode(*e->getNode(i),
+                           mesh->getElementsConnectedToNode(*e->getNode(i))));
         }
     }
 
@@ -54,12 +58,12 @@ TEST(MeshLib, QuadraticOrderMesh_Line)
     {
         if (node->getID() == 1)
         {
-            ASSERT_EQ(2u, node->getElements().size());
+            ASSERT_EQ(2u, mesh->getElementsConnectedToNode(*node).size());
             ASSERT_EQ(5u, connections[node->getID()].size());
         }
         else
         {
-            ASSERT_EQ(1u, node->getElements().size());
+            ASSERT_EQ(1u, mesh->getElementsConnectedToNode(*node).size());
             ASSERT_EQ(3u, connections[node->getID()].size());
         }
     }
@@ -85,12 +89,16 @@ TEST(MeshLib, QuadraticOrderMesh_Quad8)
 
         for (unsigned i = 0; i < e->getNumberOfBaseNodes(); i++)
         {
-            ASSERT_TRUE(isBaseNode(*e->getNode(i)));
+            ASSERT_TRUE(
+                isBaseNode(*e->getNode(i),
+                           mesh->getElementsConnectedToNode(*e->getNode(i))));
         }
         for (unsigned i = e->getNumberOfBaseNodes(); i < e->getNumberOfNodes();
              i++)
         {
-            ASSERT_FALSE(isBaseNode(*e->getNode(i)));
+            ASSERT_FALSE(
+                isBaseNode(*e->getNode(i),
+                           mesh->getElementsConnectedToNode(*e->getNode(i))));
         }
     }
 
@@ -98,33 +106,39 @@ TEST(MeshLib, QuadraticOrderMesh_Quad8)
     auto const& connections = MeshLib::calculateNodesConnectedByElements(*mesh);
 
     // Count nodes shared by four elements and also connected to all 21 nodes.
-    ASSERT_EQ(1,
-              std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
-                            [&connections](Node* const n)
-                            {
-                                return (n->getElements().size() == 4) &&
-                                       (connections[n->getID()].size() == 21);
-                            }));
+    ASSERT_EQ(
+        1,
+        std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+                      [&connections, &mesh](Node* const n)
+                      {
+                          return (mesh->getElementsConnectedToNode(*n).size() ==
+                                  4) &&
+                                 (connections[n->getID()].size() == 21);
+                      }));
 
     // Count nodes belonging to one element and also connected to all 8 nodes of
     // that corner element.
-    ASSERT_EQ(12,
-              std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
-                            [&connections](Node* const n)
-                            {
-                                return (n->getElements().size() == 1) &&
-                                       (connections[n->getID()].size() == 8);
-                            }));
+    ASSERT_EQ(
+        12,
+        std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+                      [&connections, &mesh](Node* const n)
+                      {
+                          return (mesh->getElementsConnectedToNode(*n).size() ==
+                                  1) &&
+                                 (connections[n->getID()].size() == 8);
+                      }));
 
     // Count nodes shared by two elements and also connected to the 13 nodes of
     // the two elements.
-    ASSERT_EQ(8,
-              std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
-                            [&connections](Node* const n)
-                            {
-                                return (n->getElements().size() == 2) &&
-                                       (connections[n->getID()].size() == 13);
-                            }));
+    ASSERT_EQ(
+        8,
+        std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+                      [&connections, &mesh](Node* const n)
+                      {
+                          return (mesh->getElementsConnectedToNode(*n).size() ==
+                                  2) &&
+                                 (connections[n->getID()].size() == 13);
+                      }));
 }
 
 TEST(MeshLib, QuadraticOrderMesh_Quad9)
@@ -147,12 +161,16 @@ TEST(MeshLib, QuadraticOrderMesh_Quad9)
 
         for (unsigned i = 0; i < e->getNumberOfBaseNodes(); i++)
         {
-            ASSERT_TRUE(isBaseNode(*e->getNode(i)));
+            ASSERT_TRUE(
+                isBaseNode(*e->getNode(i),
+                           mesh->getElementsConnectedToNode(*e->getNode(i))));
         }
         for (unsigned i = e->getNumberOfBaseNodes(); i < e->getNumberOfNodes();
              i++)
         {
-            ASSERT_FALSE(isBaseNode(*e->getNode(i)));
+            ASSERT_FALSE(
+                isBaseNode(*e->getNode(i),
+                           mesh->getElementsConnectedToNode(*e->getNode(i))));
         }
     }
 
@@ -161,32 +179,38 @@ TEST(MeshLib, QuadraticOrderMesh_Quad9)
     auto const& connections = MeshLib::calculateNodesConnectedByElements(*mesh);
     // Count nodes shared by four elements and also connected to all 25 nodes.
     ASSERT_EQ(
-        1, std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
-                         [&connections](Node* const n)
-                         {
-                             return (n->getElements().size() == 4) &&
-                                    (connections[n->getID()].size() == 21 + 4);
-                         }));
+        1,
+        std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+                      [&connections, &mesh](Node* const n)
+                      {
+                          return (mesh->getElementsConnectedToNode(*n).size() ==
+                                  4) &&
+                                 (connections[n->getID()].size() == 21 + 4);
+                      }));
 
     // Count nodes belonging to one element and also connected to all 9 nodes of
     // that corner element---three corners and the centre node.
-    ASSERT_EQ(12 + 4,
-              std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
-                            [&connections](Node* const n)
-                            {
-                                return (n->getElements().size() == 1) &&
-                                       (connections[n->getID()].size() == 9);
-                            }));
+    ASSERT_EQ(
+        12 + 4,
+        std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+                      [&connections, &mesh](Node* const n)
+                      {
+                          return (mesh->getElementsConnectedToNode(*n).size() ==
+                                  1) &&
+                                 (connections[n->getID()].size() == 9);
+                      }));
 
     // Count nodes shared by two elements and also connected to the 15 nodes of
     // the two elements.
     ASSERT_EQ(
-        8, std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
-                         [&connections](Node* const n)
-                         {
-                             return (n->getElements().size() == 2) &&
-                                    (connections[n->getID()].size() == 13 + 2);
-                         }));
+        8,
+        std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+                      [&connections, &mesh](Node* const n)
+                      {
+                          return (mesh->getElementsConnectedToNode(*n).size() ==
+                                  2) &&
+                                 (connections[n->getID()].size() == 13 + 2);
+                      }));
 }
 
 TEST(MeshLib, QuadraticOrderMesh_LineQuad)
@@ -243,12 +267,16 @@ TEST(MeshLib, QuadraticOrderMesh_LineQuad)
 
         for (unsigned i = 0; i < e->getNumberOfBaseNodes(); i++)
         {
-            ASSERT_TRUE(isBaseNode(*e->getNode(i)));
+            ASSERT_TRUE(
+                isBaseNode(*e->getNode(i),
+                           mesh->getElementsConnectedToNode(*e->getNode(i))));
         }
         for (unsigned i = e->getNumberOfBaseNodes(); i < e->getNumberOfNodes();
              i++)
         {
-            ASSERT_FALSE(isBaseNode(*e->getNode(i)));
+            ASSERT_FALSE(
+                isBaseNode(*e->getNode(i),
+                           mesh->getElementsConnectedToNode(*e->getNode(i))));
         }
     }
 
@@ -256,41 +284,49 @@ TEST(MeshLib, QuadraticOrderMesh_LineQuad)
 
     auto const& connections = MeshLib::calculateNodesConnectedByElements(*mesh);
     // Count nodes shared by six elements and also connected to all 21 nodes.
-    ASSERT_EQ(1,
-              std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
-                            [&connections](Node* const n)
-                            {
-                                return (n->getElements().size() == 6) &&
-                                       (connections[n->getID()].size() == 21);
-                            }));
+    ASSERT_EQ(
+        1,
+        std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+                      [&connections, &mesh](Node* const n)
+                      {
+                          return (mesh->getElementsConnectedToNode(*n).size() ==
+                                  6) &&
+                                 (connections[n->getID()].size() == 21);
+                      }));
 
     // Count nodes belonging to one element and also connected to all 8 nodes of
     // that corner element.
-    ASSERT_EQ(12,
-              std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
-                            [&connections](Node* const n)
-                            {
-                                return (n->getElements().size() == 1) &&
-                                       (connections[n->getID()].size() == 8);
-                            }));
+    ASSERT_EQ(
+        12,
+        std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+                      [&connections, &mesh](Node* const n)
+                      {
+                          return (mesh->getElementsConnectedToNode(*n).size() ==
+                                  1) &&
+                                 (connections[n->getID()].size() == 8);
+                      }));
 
     // Count nodes shared by three elements (quads and the line) and also
     // connected to the 13 nodes of the two elements.
-    ASSERT_EQ(4,
-              std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
-                            [&connections](Node* const n)
-                            {
-                                return (n->getElements().size() == 3) &&
-                                       (connections[n->getID()].size() == 13);
-                            }));
+    ASSERT_EQ(
+        4,
+        std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+                      [&connections, &mesh](Node* const n)
+                      {
+                          return (mesh->getElementsConnectedToNode(*n).size() ==
+                                  3) &&
+                                 (connections[n->getID()].size() == 13);
+                      }));
 
     // Count nodes shared by two elements (quads) and also connected to the 13
     // nodes of the two elements.
-    ASSERT_EQ(4,
-              std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
-                            [&connections](Node* const n)
-                            {
-                                return (n->getElements().size() == 2) &&
-                                       (connections[n->getID()].size() == 13);
-                            }));
+    ASSERT_EQ(
+        4,
+        std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+                      [&connections, &mesh](Node* const n)
+                      {
+                          return (mesh->getElementsConnectedToNode(*n).size() ==
+                                  2) &&
+                                 (connections[n->getID()].size() == 13);
+                      }));
 }
diff --git a/Tests/MeshLib/TestTriLineMesh.cpp b/Tests/MeshLib/TestTriLineMesh.cpp
index 776ec355a4f..29002a1d791 100644
--- a/Tests/MeshLib/TestTriLineMesh.cpp
+++ b/Tests/MeshLib/TestTriLineMesh.cpp
@@ -64,9 +64,9 @@ public:
 
     bool isConnectedToNode(std::size_t const n, std::size_t const e) const
     {
-        return std::find(nodes[n]->getElements().cbegin(),
-                         nodes[n]->getElements().cend(),
-                         elements[e]) != nodes[n]->getElements().cend();
+        auto const& node_elements = mesh->getElementsConnectedToNode(n);
+        return std::find(node_elements.cbegin(), node_elements.cend(),
+                         elements[e]) != node_elements.cend();
     }
 };
 
@@ -83,18 +83,18 @@ TEST_F(MeshLibTriLineMesh, Construction)
 TEST_F(MeshLibTriLineMesh, NodeToElementConnectivity)
 {
     // Nodes 0 and 3 are connected only to triangles.
-    EXPECT_EQ(1u, nodes[0]->getNumberOfElements());
-    EXPECT_EQ(1u, nodes[3]->getNumberOfElements());
-    EXPECT_EQ(elements[0], nodes[0]->getElement(0));
-    EXPECT_EQ(elements[1], nodes[3]->getElement(0));
+    EXPECT_EQ(1u, mesh->getElementsConnectedToNode(*nodes[0]).size());
+    EXPECT_EQ(1u, mesh->getElementsConnectedToNode(*nodes[3]).size());
+    EXPECT_EQ(elements[0], mesh->getElementsConnectedToNode(*nodes[0])[0]);
+    EXPECT_EQ(elements[1], mesh->getElementsConnectedToNode(*nodes[3])[0]);
 
     // Nodes 1 and 2 are connected to all elements.
-    EXPECT_EQ(3u, nodes[1]->getNumberOfElements());
+    EXPECT_EQ(3u, mesh->getElementsConnectedToNode(*nodes[1]).size());
     EXPECT_TRUE(isConnectedToNode(1, 0));
     EXPECT_TRUE(isConnectedToNode(1, 1));
     EXPECT_TRUE(isConnectedToNode(1, 2));
 
-    EXPECT_EQ(3u, nodes[2]->getNumberOfElements());
+    EXPECT_EQ(3u, mesh->getElementsConnectedToNode(*nodes[2]).size());
     EXPECT_TRUE(isConnectedToNode(2, 0));
     EXPECT_TRUE(isConnectedToNode(2, 1));
     EXPECT_TRUE(isConnectedToNode(2, 2));
-- 
GitLab