diff --git a/MeshLib/Mesh.cpp b/MeshLib/Mesh.cpp
index 8b043cf8ae04d64679d49430c4da3b9a789f8028..b12854308095f2016a5c0f8acd713ab20d04e9ff 100644
--- a/MeshLib/Mesh.cpp
+++ b/MeshLib/Mesh.cpp
@@ -42,6 +42,8 @@ Mesh::Mesh(const std::string &name,
     assert(n_base_nodes <= nodes.size());
     this->resetNodeIDs();
     this->resetElementIDs();
+    if (isNonlinear())
+        this->checkNonlinearNodeIDs();
     this->setDimension();
     this->setElementsConnectedToNodes();
     //this->setNodesConnectedByEdges();
@@ -253,6 +255,29 @@ void Mesh::setNodesConnectedByElements()
     }
 }
 
+void Mesh::checkNonlinearNodeIDs() const
+{
+    for (MeshLib::Element const* e : _elements)
+    {
+        for (unsigned i=0; i<e->getNumberOfBaseNodes(); i++)
+        {
+            if (!(e->getNodeIndex(i) < getNumberOfBaseNodes()))
+                OGS_FATAL(
+                    "Node %d is a base/linear node, but the ID is not smaller "
+                    "than the number of base nodes %d. Please renumber node IDs in the mesh.",
+                    e->getNodeIndex(i), getNumberOfBaseNodes());
+        }
+        for (unsigned i=e->getNumberOfBaseNodes(); i<e->getNumberOfNodes(); i++)
+        {
+            if (!(e->getNodeIndex(i) >= getNumberOfBaseNodes()))
+                OGS_FATAL(
+                    "Node %d is a non-linear node, but the ID is smaller "
+                    "than the number of base nodes %d. Please renumber node IDs in the mesh.",
+                    e->getNodeIndex(i), getNumberOfBaseNodes());
+        }
+    }
+}
+
 void scaleMeshPropertyVector(MeshLib::Mesh & mesh,
                              std::string const& property_name,
                              double factor)
diff --git a/MeshLib/Mesh.h b/MeshLib/Mesh.h
index 0dbbc93cf0f2ded867f8a4c98e8bc8f433be9fcb..a3e1907049604ae289d4f4f84900a71ebbcd3eeb 100644
--- a/MeshLib/Mesh.h
+++ b/MeshLib/Mesh.h
@@ -162,6 +162,9 @@ protected:
     /// connected if they are shared by an element.
     void setNodesConnectedByElements();
 
+    /// Check if all the nonlinear nodes are stored at the end of the node vector
+    void checkNonlinearNodeIDs() const;
+
     std::size_t const _id;
     unsigned _mesh_dimension;
     /// The minimal and maximal edge length over all elements in the mesh