diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 5d24b799..898c79de 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -89,7 +89,7 @@ jobs:
       - name: 🏗️ Compile (other)
         if: ${{ matrix.os != 'windows-latest' }}
         run: |
-          cmake -DVIENNALS_BUILD_TESTS=ON -B build
+          cmake -DVIENNALS_BUILD_TESTS=ON -DCMAKE_BUILD_TYPE=${{ matrix.config }} -B build
           cmake --build build --config ${{ matrix.config }}
 
       - name: 🧪 Test
diff --git a/CMakeLists.txt b/CMakeLists.txt
index caaa0b7c..39d4883f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 project(
   ViennaLS
   LANGUAGES CXX
-  VERSION 5.8.5)
+  VERSION 5.9.0)
 
 # --------------------------------------------------------------------------------------------------------
 # Library options
@@ -521,6 +521,11 @@ endif()
 # Install Target
 # --------------------------------------------------------------------------------------------------------
 
+set(VIENNALS_DEPENDENCIES "ViennaHRLE;ViennaCore")
+if(VIENNALS_USE_VTK)
+  list(APPEND VIENNALS_DEPENDENCIES "VTK")
+endif()
+
 packageProject(
   NAME ${PROJECT_NAME} NAMESPACE ViennaTools
   VERSION ${PROJECT_VERSION}
@@ -528,4 +533,4 @@ packageProject(
   INCLUDE_DIR ${PROJECT_SOURCE_DIR}/include/viennals
   INCLUDE_DESTINATION include/viennals-${PROJECT_VERSION}
   COMPATIBILITY SameMajorVersion
-  DEPENDENCIES "VTK;ViennaHRLE;ViennaCore")
+  DEPENDENCIES ${VIENNALS_DEPENDENCIES})
diff --git a/README.md b/README.md
index 63197824..631fa8d6 100644
--- a/README.md
+++ b/README.md
@@ -90,9 +90,6 @@ This will install the necessary headers and CMake files to the specified path. I
 
 In order to install ViennaLS without VTK, run:
 ```bash
-git clone https://github.com/ViennaTools/ViennaLS.git
-cd ViennaLS
-
 cmake -B build -D CMAKE_INSTALL_PREFIX=/path/to/your/custom/install/ -D VIENNALS_USE_VTK=OFF
 cmake --install build
 ```
@@ -152,9 +149,6 @@ ViennaLS uses CTest to run its tests.
 In order to check whether ViennaLS runs without issues on your system, you can run:
 
 ```bash
-git clone https://github.com/ViennaTools/ViennaLS.git
-cd ViennaLS
-
 cmake -B build -DVIENNALS_BUILD_TESTS=ON
 cmake --build build
 ctest -E "Benchmark|Performance" --test-dir build
@@ -175,7 +169,9 @@ We recommend using [CPM.cmake](https://github.com/cpm-cmake/CPM.cmake) to consum
 
 * Installation with CPM
   ```cmake
-  CPMAddPackage("gh:viennatools/viennals@5.8.5")
+  CPMAddPackage("gh:viennatools/viennals@5.9.0")
+
+  target_link_libraries(${PROJECT_NAME} PUBLIC ViennaTools::ViennaLS)
   ```
 
 * With a local installation
@@ -216,10 +212,6 @@ cmake --build build --target format
 
 ## Authors
 
-Current contributors: Tobias Reiter, Roman Kostal, Lado Filipovic
-
-Founder and initial developer: Otmar Ertl
-
 Contact us via: viennatools@iue.tuwien.ac.at
 
 ViennaLS was developed under the aegis of the 'Institute for Microelectronics' at the 'TU Wien'.  
diff --git a/include/viennals/lsAdvect.hpp b/include/viennals/lsAdvect.hpp
index 0d9ad9b5..0cf97aab 100644
--- a/include/viennals/lsAdvect.hpp
+++ b/include/viennals/lsAdvect.hpp
@@ -98,13 +98,10 @@ template <class T, int D> class Advect {
 
     VectorType<T, D> finalAlphas{};
 
-#pragma omp parallel num_threads((levelSets.back())->getNumberOfSegments())
-    {
+#pragma omp parallel for
+    for (unsigned p = 0; p < levelSets.back()->getNumberOfSegments(); ++p) {
       VectorType<T, D> localAlphas{};
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+
       viennahrle::Index<D> startVector =
           (p == 0) ? grid.getMinGridPoint()
                    : topDomain.getSegmentation()[p - 1];
@@ -274,12 +271,8 @@ template <class T, int D> class Advect {
     }
 #endif
 
-#pragma omp parallel num_threads(newDomain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < newDomain.getNumberOfSegments(); ++p) {
       auto &domainSegment = newDomain.getDomainSegment(p);
 
       viennahrle::Index<D> startVector =
@@ -287,7 +280,7 @@ template <class T, int D> class Advect {
                    : newDomain.getSegmentation()[p - 1];
 
       viennahrle::Index<D> endVector =
-          (p != static_cast<int>(newDomain.getNumberOfSegments() - 1))
+          (p != newDomain.getNumberOfSegments() - 1)
               ? newDomain.getSegmentation()[p]
               : grid.incrementIndices(grid.getMaxGridPoint());
 
@@ -461,18 +454,15 @@ template <class T, int D> class Advect {
 
     storedRates.resize(topDomain.getNumberOfSegments());
 
-#pragma omp parallel num_threads(topDomain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < topDomain.getNumberOfSegments(); ++p) {
+
       viennahrle::Index<D> startVector =
           (p == 0) ? grid.getMinGridPoint()
                    : topDomain.getSegmentation()[p - 1];
 
       viennahrle::Index<D> endVector =
-          (p != static_cast<int>(topDomain.getNumberOfSegments() - 1))
+          (p != topDomain.getNumberOfSegments() - 1)
               ? topDomain.getSegmentation()[p]
               : grid.incrementIndices(grid.getMaxGridPoint());
 
@@ -707,12 +697,9 @@ template <class T, int D> class Advect {
 
     const bool checkDiss = checkDissipation;
 
-#pragma omp parallel num_threads(topDomain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < topDomain.getNumberOfSegments(); ++p) {
+
       auto itRS = storedRates[p].cbegin();
       auto &segment = topDomain.getDomainSegment(p);
       const unsigned maxId = segment.getNumberOfPoints();
@@ -738,8 +725,8 @@ template <class T, int D> class Advect {
         T velocity = gradient - dissipation;
         // check if dissipation is too high and would cause a change in
         // direction of the velocity
-        if (checkDiss && (gradient < 0 && velocity > 0) ||
-            (gradient > 0 && velocity < 0)) {
+        if (checkDiss && ((gradient < 0 && velocity > 0) ||
+                          (gradient > 0 && velocity < 0))) {
           velocity = 0;
         }
 
@@ -751,8 +738,8 @@ template <class T, int D> class Advect {
 
           // recalculate velocity and rate
           velocity = itRS->first.first - itRS->first.second;
-          if (checkDiss && (itRS->first.first < 0 && velocity > 0) ||
-              (itRS->first.first > 0 && velocity < 0)) {
+          if (checkDiss && ((itRS->first.first < 0 && velocity > 0) ||
+                            (itRS->first.first > 0 && velocity < 0))) {
             velocity = 0;
           }
           rate = time * velocity;
diff --git a/include/viennals/lsBooleanOperation.hpp b/include/viennals/lsBooleanOperation.hpp
index 2595a95d..7d4c6fe5 100644
--- a/include/viennals/lsBooleanOperation.hpp
+++ b/include/viennals/lsBooleanOperation.hpp
@@ -74,12 +74,8 @@ template <class T, int D> class BooleanOperation {
       newDataLS.resize(newDataSourceIds.size());
     }
 
-#pragma omp parallel num_threads(newDomain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < newDomain.getNumberOfSegments(); ++p) {
 
       auto &domainSegment = newDomain.getDomainSegment(p);
 
@@ -214,12 +210,8 @@ template <class T, int D> class BooleanOperation {
 
   void invert() {
     auto &hrleDomain = levelSetA->getDomain();
-#pragma omp parallel num_threads(hrleDomain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < hrleDomain.getNumberOfSegments(); ++p) {
       auto &domainSegment = hrleDomain.getDomainSegment(p);
 
       // change all defined values
diff --git a/include/viennals/lsCalculateCurvatures.hpp b/include/viennals/lsCalculateCurvatures.hpp
index 3d53499a..80176088 100644
--- a/include/viennals/lsCalculateCurvatures.hpp
+++ b/include/viennals/lsCalculateCurvatures.hpp
@@ -89,12 +89,8 @@ template <class T, int D> class CalculateCurvatures {
         (type == CurvatureEnum::MEAN_AND_GAUSSIAN_CURVATURE);
 
     //! Calculate Curvatures
-#pragma omp parallel num_threads(levelSet->getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < levelSet->getNumberOfSegments(); ++p) {
 
       auto &meanCurvatures = meanCurvaturesVector[p];
       auto &gaussCurvatures = gaussCurvaturesVector[p];
diff --git a/include/viennals/lsCalculateNormalVectors.hpp b/include/viennals/lsCalculateNormalVectors.hpp
index fd8f8473..e13f9266 100644
--- a/include/viennals/lsCalculateNormalVectors.hpp
+++ b/include/viennals/lsCalculateNormalVectors.hpp
@@ -131,12 +131,8 @@ template <class T, int D> class CalculateNormalVectors {
     auto grid = levelSet->getGrid();
 
     // Calculate Normalvectors
-#pragma omp parallel num_threads(levelSet->getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < levelSet->getNumberOfSegments(); ++p) {
 
       auto &normalVectors = normalVectorsVector[p];
       normalVectors.reserve(pointsPerSegment);
@@ -209,12 +205,8 @@ template <class T, int D> class CalculateNormalVectors {
     // points.
     std::vector<Vec3D<T>> normalVectors(levelSet->getNumberOfPoints());
 
-#pragma omp parallel num_threads(levelSet->getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < levelSet->getNumberOfSegments(); ++p) {
 
       viennahrle::Index<D> startVector =
           (p == 0) ? grid.getMinGridPoint()
diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp
index 45bee138..4f68c9c9 100644
--- a/include/viennals/lsCalculateVisibilities.hpp
+++ b/include/viennals/lsCalculateVisibilities.hpp
@@ -65,12 +65,8 @@ template <class T, int D> class CalculateVisibilities {
       }
       //****************************
 
-#pragma omp parallel num_threads(levelSet->getNumberOfSegments())
-      {
-        int p = 0;
-#ifdef _OPENMP
-        p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+      for (unsigned p = 0; p < domain.getNumberOfSegments(); ++p) {
 
         const viennahrle::Index<D> startVector =
             (p == 0) ? grid.getMinGridPoint() : domain.getSegmentation()[p - 1];
diff --git a/include/viennals/lsDetectFeatures.hpp b/include/viennals/lsDetectFeatures.hpp
index a69de7ea..e843ea5c 100644
--- a/include/viennals/lsDetectFeatures.hpp
+++ b/include/viennals/lsDetectFeatures.hpp
@@ -94,12 +94,8 @@ template <class T, int D> class DetectFeatures {
     typename Domain<T, D>::DomainType &domain = levelSet->getDomain();
     std::vector<std::vector<T>> flagsReserve(levelSet->getNumberOfSegments());
 
-#pragma omp parallel num_threads((levelSet)->getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < levelSet->getNumberOfSegments(); ++p) {
 
       auto &flagsSegment = flagsReserve[p];
       flagsSegment.reserve(
@@ -169,12 +165,8 @@ template <class T, int D> class DetectFeatures {
     std::vector<std::vector<T>> flagsReserve(levelSet->getNumberOfSegments());
 
     // Compare angles between normal vectors
-#pragma omp parallel num_threads((levelSet)->getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < levelSet->getNumberOfSegments(); ++p) {
 
       Vec3D<T> zeroVector{};
 
diff --git a/include/viennals/lsExpand.hpp b/include/viennals/lsExpand.hpp
index 3f9fa8cc..83e2a1c4 100644
--- a/include/viennals/lsExpand.hpp
+++ b/include/viennals/lsExpand.hpp
@@ -81,12 +81,8 @@ template <class T, int D> class Expand {
       if (updateData)
         newDataSourceIds.resize(newDomain.getNumberOfSegments());
 
-#pragma omp parallel num_threads(newDomain.getNumberOfSegments())
-      {
-        int p = 0;
-#ifdef _OPENMP
-        p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+      for (unsigned p = 0; p < newDomain.getNumberOfSegments(); ++p) {
 
         auto &domainSegment = newDomain.getDomainSegment(p);
 
diff --git a/include/viennals/lsGeometricAdvect.hpp b/include/viennals/lsGeometricAdvect.hpp
index 66f4d62b..33dcf4c9 100644
--- a/include/viennals/lsGeometricAdvect.hpp
+++ b/include/viennals/lsGeometricAdvect.hpp
@@ -297,12 +297,8 @@ template <class T, int D> class GeometricAdvect {
     }
 #endif
 // set up multithreading
-#pragma omp parallel num_threads(domain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < domain.getNumberOfSegments(); ++p) {
 
       viennahrle::Index<D> startVector;
       if (p == 0) {
diff --git a/include/viennals/lsInterior.hpp b/include/viennals/lsInterior.hpp
index 47df1cf4..d42faadb 100644
--- a/include/viennals/lsInterior.hpp
+++ b/include/viennals/lsInterior.hpp
@@ -73,14 +73,8 @@ template <class T, int D> class Interior {
       if (updateData)
         newDataSourceIds.resize(newDomain.getNumberOfSegments());
 
-#pragma omp parallel num_threads(newDomain.getNumberOfSegments())              \
-    reduction(+ : addedPoints)
-      {
-        int p = 0;
-#ifdef _OPENMP
-        p = omp_get_thread_num();
-#endif
-
+#pragma omp parallel for reduction(+ : addedPoints)
+      for (unsigned p = 0; p < newDomain.getNumberOfSegments(); ++p) {
         auto &domainSegment = newDomain.getDomainSegment(p);
 
         viennahrle::Index<D> const startVector =
diff --git a/include/viennals/lsMaterialMap.hpp b/include/viennals/lsMaterialMap.hpp
index 8adb907a..8e8e2625 100644
--- a/include/viennals/lsMaterialMap.hpp
+++ b/include/viennals/lsMaterialMap.hpp
@@ -45,6 +45,9 @@ class MaterialMap {
   void setMaterialId(const std::size_t index, const int materialId) {
     if (index >= materialMap.size()) {
       materialMap.resize(index + 1, -1); // Initialize new elements with -1
+    } else {
+      auto oldMaterialId = materialMap[index];
+      materials.erase(oldMaterialId); // Remove old material ID if it exists
     }
     materialMap[index] = materialId;
     materials.insert(materialId);
diff --git a/include/viennals/lsMesh.hpp b/include/viennals/lsMesh.hpp
index 6f95e5ff..15ec19bf 100644
--- a/include/viennals/lsMesh.hpp
+++ b/include/viennals/lsMesh.hpp
@@ -3,7 +3,11 @@
 #include <lsPreCompileMacros.hpp>
 
 #include <array>
+#include <cmath>
+#include <functional>
 #include <iostream>
+#include <stdexcept>
+#include <unordered_map>
 #include <vector>
 
 #include <vcPointData.hpp>
@@ -34,33 +38,6 @@ template <class T = double> class Mesh {
   constexpr static const char *materialIdsLabel = "MaterialIds";
   constexpr static const char *normalsLabel = "Normals";
 
-private:
-  // iterator typedef
-  using VectorIt = typename PointData<T>::VectorDataType::iterator;
-  // find function to avoid including the whole algorithm header
-  static VectorIt find(VectorIt first, VectorIt last, const Vec3D<T> &value) {
-    for (; first != last; ++first) {
-      if (*first == value) {
-        return first;
-      }
-    }
-    return last;
-  }
-
-  // helper function for duplicate removal
-  template <class ElementType>
-  static void replaceNode(ElementType &elements,
-                          std::pair<unsigned, unsigned> node) {
-    for (unsigned i = 0; i < elements.size(); ++i) {
-      for (unsigned j = 0; j < elements[i].size(); ++j) {
-        if (elements[i][j] == node.first) {
-          elements[i][j] = node.second;
-        }
-      }
-    }
-  };
-
-public:
   // Convenience function to create a new mesh smart pointer.
   static auto New() { return SmartPointer<Mesh>::New(); }
 
@@ -166,39 +143,88 @@ template <class T = double> class Mesh {
     return hexas.size() - 1;
   }
 
+  /// Remove exactly equal nodes, preserving first-occurrence order and the
+  /// first node's scalar/vector point data. Remap all element node IDs;
+  /// elements and cell data are retained, including degenerate elements.
+  /// Nodes containing NaNs remain distinct. Expected linear time in the number
+  /// of nodes, connectivity entries, and point-data values, with linear scratch
+  /// storage. When nodes are merged, point-data arrays must contain one value
+  /// per input node; otherwise throws std::invalid_argument without changing
+  /// the mesh.
   void removeDuplicateNodes() {
+    if (nodes.size() < 2)
+      return;
+
+    struct NodeHash {
+      std::size_t operator()(const Vec3D<T> &node) const {
+        std::size_t seed = 0;
+        for (const T coordinate : node) {
+          seed ^= std::hash<T>{}(coordinate) + std::size_t(0x9e3779b9) +
+                  (seed << 6) + (seed >> 2);
+        }
+        return seed;
+      }
+    };
+
+    // Keep full coordinates as keys so hash collisions cannot merge nodes.
+    // std::hash<T> also gives equal hashes for +0 and -0, which compare equal.
+    std::unordered_map<Vec3D<T>, unsigned, NodeHash> uniqueNodes;
+    uniqueNodes.reserve(nodes.size());
     std::vector<Vec3D<T>> newNodes;
-    // can just push first point since it cannot be duplicate
-    newNodes.push_back(nodes[0]);
-    // now check for duplicates
-    // pair of oldId <-> newId
-    std::vector<std::pair<unsigned, unsigned>> duplicates;
-    bool adjusted = false;
-    for (unsigned i = 1; i < nodes.size(); ++i) {
-      auto it = find(newNodes.begin(), newNodes.end(), nodes[i]);
-      if (it != newNodes.end()) {
-        adjusted = true;
-        // if duplicate point, save it to be replaced
-        unsigned nodeId =
-            static_cast<unsigned>(std::distance(newNodes.begin(), it));
-        duplicates.emplace_back(i, nodeId);
+    newNodes.reserve(nodes.size());
+    std::vector<unsigned> oldToNew(nodes.size());
+    std::vector<unsigned> retainedIndices;
+    retainedIndices.reserve(nodes.size());
+
+    for (std::size_t i = 0; i < nodes.size(); ++i) {
+      const auto &node = nodes[i];
+      const auto newId = static_cast<unsigned>(newNodes.size());
+      // NaNs are not equal to themselves and cannot serve as hash-table keys.
+      if (!std::isnan(node[0]) && !std::isnan(node[1]) &&
+          !std::isnan(node[2])) {
+        const auto result = uniqueNodes.try_emplace(node, newId);
+        oldToNew[i] = result.first->second;
+        if (!result.second)
+          continue;
       } else {
-        if (adjusted)
-          duplicates.emplace_back(i, static_cast<unsigned>(newNodes.size()));
-        newNodes.push_back(nodes[i]);
+        oldToNew[i] = newId;
       }
+      newNodes.push_back(node);
+      retainedIndices.push_back(static_cast<unsigned>(i));
     }
-    nodes = newNodes;
-
-    // now replace in vertices
-    // TODO also need to shift down all other nodes
-    for (auto &duplicate : duplicates) {
-      replaceNode(vertices, duplicate);
-      replaceNode(lines, duplicate);
-      replaceNode(triangles, duplicate);
-      replaceNode(tetras, duplicate);
-      replaceNode(hexas, duplicate);
-    }
+
+    if (newNodes.size() == nodes.size())
+      return;
+
+    // Validate before translating data or changing any mesh connectivity.
+    const auto validateData = [this](const auto &arrays) {
+      for (const auto &data : arrays) {
+        if (data.size() != nodes.size()) {
+          throw std::invalid_argument(
+              "Mesh::removeDuplicateNodes: point-data size must match the "
+              "number of nodes.");
+        }
+      }
+    };
+    validateData(pointData.getScalarData());
+    validateData(pointData.getVectorData());
+    PointData<T> newPointData;
+    newPointData.translateFromData(pointData, retainedIndices);
+
+    const auto remapElements = [&oldToNew](auto &elements) {
+      for (auto &element : elements) {
+        for (auto &nodeId : element)
+          nodeId = oldToNew[nodeId];
+      }
+    };
+    remapElements(vertices);
+    remapElements(lines);
+    remapElements(triangles);
+    remapElements(tetras);
+    remapElements(hexas);
+
+    nodes = std::move(newNodes);
+    pointData = std::move(newPointData);
   }
 
   void append(const Mesh<T> &passedMesh) {
diff --git a/include/viennals/lsOxidation.hpp b/include/viennals/lsOxidation.hpp
index ce29cb45..c4a96f45 100644
--- a/include/viennals/lsOxidation.hpp
+++ b/include/viennals/lsOxidation.hpp
@@ -150,6 +150,8 @@ LOCOSConservationDiagnostics<T> computeLOCOSOpenWindowConservation(
 /// @endcode
 template <class T, int D> class Oxidation {
   using IndexType = viennahrle::Index<D>;
+  using IndexCacheMap =
+      std::unordered_map<IndexType, T, typename IndexType::hash>;
 
   SmartPointer<Domain<T, D>> siInterface = nullptr;
   SmartPointer<Domain<T, D>> ambientInterface = nullptr;
@@ -184,13 +186,13 @@ template <class T, int D> class Oxidation {
 
   GpuMode gpuMode_ = GpuMode::Cpu;
   GpuPreconditioner gpuPreconditioner_ = GpuPreconditioner::Jacobi;
-  std::unordered_map<std::size_t, T> concentrationCache_;
+  IndexCacheMap concentrationCache_;
 
 public:
-  const std::unordered_map<std::size_t, T> &getConcentrationCache() const {
+  const IndexCacheMap &getConcentrationCache() const {
     return concentrationCache_;
   }
-  void setConcentrationCache(std::unordered_map<std::size_t, T> cache) {
+  void setConcentrationCache(IndexCacheMap cache) {
     concentrationCache_ = std::move(cache);
   }
 
diff --git a/include/viennals/lsOxidationDeformation.hpp b/include/viennals/lsOxidationDeformation.hpp
index f8e44f82..4736c2cc 100644
--- a/include/viennals/lsOxidationDeformation.hpp
+++ b/include/viennals/lsOxidationDeformation.hpp
@@ -130,7 +130,7 @@ class OxidationDeformation final : public VelocityField<T>,
   bool nodesDirty_ = true;
   std::array<T, D> maxVelocity_{};
   bool useRequestedBounds = false;
-  std::unordered_map<IndexType, std::array<T, 9>, detail::IndexTypeHasher<D>>
+  std::unordered_map<IndexType, std::array<T, 9>, typename IndexType::hash>
       deviatoricStressHistory;
 
   // Warm-start storage: solutions from previous time step used as initial guess
@@ -1133,25 +1133,24 @@ class OxidationDeformation final : public VelocityField<T>,
         neighbor[direction] += offset;
 
         if (!inBounds(neighbor)) {
-          detail::vecAddTo(sum, toT(v[nodeId])); // zero-flux: ghost = self
+          sum = sum + toT(v[nodeId]); // zero-flux: ghost = self
           continue;
         }
 
         const std::size_t neighborId = lookupNode(neighbor);
         if (neighborId != noNode) {
-          detail::vecAddTo(sum, toT(v[neighborId]));
+          sum = sum + toT(v[neighborId]);
           continue;
         }
 
         const unsigned fi = direction * 2u + (offset == 1 ? 1u : 0u);
         const Boundary boundary = faceBCTypes_[fi * nodes.size() + nodeId];
         if (boundary == Boundary::REACTION) {
-          detail::vecAddTo(sum, reactionBoundaryVelocity(node.index));
+          sum = sum + reactionBoundaryVelocity(node.index);
         } else if (boundary == Boundary::MASK) {
-          detail::vecAddTo(sum,
-                           maskVelocityBoundary(node.index, toT(v[nodeId])));
+          sum = sum + maskVelocityBoundary(node.index, toT(v[nodeId]));
         } else {
-          detail::vecAddTo(sum, toT(v[nodeId])); // AMBIENT/NONE: zero-flux
+          sum = sum + toT(v[nodeId]); // AMBIENT/NONE: zero-flux
         }
       }
     }
@@ -2166,8 +2165,7 @@ class OxidationDeformation final : public VelocityField<T>,
       const T dSum = plus.distance + minus.distance;
       const T plusCoeff = T(2) / (plus.distance * dSum);
       const T minusCoeff = T(2) / (minus.distance * dSum);
-      detail::vecAddTo(rhs, detail::vecScaled(plus.value, plusCoeff));
-      detail::vecAddTo(rhs, detail::vecScaled(minus.value, minusCoeff));
+      rhs = rhs + plusCoeff * plus.value + minusCoeff * minus.value;
       diag += plusCoeff + minusCoeff;
     }
   }
@@ -2803,8 +2801,7 @@ class OxidationDeformation final : public VelocityField<T>,
     for (unsigned i = 0; i < D; ++i)
       coordinate[i] = index[i] * gridDelta;
     const T expansionVelocity = localExpansionSpeed(coordinate);
-    return detail::vecScaled(reactionNormal(index),
-                             reactionSign * expansionVelocity);
+    return reactionNormal(index) * (reactionSign * expansionVelocity);
   }
 
   Vec3D<T> unresolvedAmbientVelocity(const Vec3D<T> &coordinate) const {
@@ -2817,7 +2814,7 @@ class OxidationDeformation final : public VelocityField<T>,
 
     ConstSparseIterator ambientIt(ambientInterface->getDomain());
     const auto normal = levelSetNormal(ambientIt, index);
-    return detail::vecScaled(normal, localExpansionSpeed(coordinate));
+    return normal * localExpansionSpeed(coordinate);
   }
 
   Vec3D<T> estimateMaxUnresolvedAmbientVelocity() const {
@@ -3005,7 +3002,7 @@ class OxidationDeformation final : public VelocityField<T>,
       historyEntries[i] = {node.index, deviatoricStress};
     }
 
-    std::unordered_map<IndexType, std::array<T, 9>, detail::IndexTypeHasher<D>>
+    std::unordered_map<IndexType, std::array<T, 9>, typename IndexType::hash>
         nextHistory;
     nextHistory.reserve(nodes.size());
     for (const auto &entry : historyEntries)
diff --git a/include/viennals/lsOxidationDiffusion.hpp b/include/viennals/lsOxidationDiffusion.hpp
index 293a36fd..a0d5b07d 100644
--- a/include/viennals/lsOxidationDiffusion.hpp
+++ b/include/viennals/lsOxidationDiffusion.hpp
@@ -99,6 +99,8 @@ class OxidationDiffusion final : public VelocityField<T>,
   using IndexType = viennahrle::Index<D>;
   using ConstSparseIterator =
       viennahrle::ConstSparseIterator<typename Domain<T, D>::DomainType>;
+  using IndexCacheMap =
+      std::unordered_map<IndexType, T, typename IndexType::hash>;
 
 private:
   static constexpr T boltzmannConstant = T(1.380649e-23);
@@ -159,8 +161,8 @@ class OxidationDiffusion final : public VelocityField<T>,
   bool useRequestedBounds = false;
   bool warmStartable_ =
       false; // true when nodes[i].concentration holds a prior solution
-  std::unordered_map<std::size_t, T> pressureLookup;
-  std::unordered_map<std::size_t, T> concentrationCache_;
+  IndexCacheMap pressureLookup;
+  IndexCacheMap concentrationCache_;
   std::vector<Node> nodes;
   // Face-major flat BC arrays: index = fi * n + nodeId, where fi in [0, 2*D).
   // Face-major layout gives coalesced GPU reads when all warp threads access
@@ -211,6 +213,10 @@ class OxidationDiffusion final : public VelocityField<T>,
 #endif
   }
 
+  // Disable copy and assignment to avoid double-free of GPU buffers
+  OxidationDiffusion(const OxidationDiffusion &) = delete;
+  OxidationDiffusion &operator=(const OxidationDiffusion &) = delete;
+
   template <class... Args> static auto New(Args &&...args) {
     return SmartPointer<OxidationDiffusion>::New(std::forward<Args>(args)...);
   }
@@ -268,11 +274,10 @@ class OxidationDiffusion final : public VelocityField<T>,
   }
 
   void setPressure(const IndexType &index, T pressure) {
-    const auto key = detail::gridIndexHash<D>(index);
     if (std::isfinite(pressure))
-      pressureLookup[key] = pressure;
+      pressureLookup[index] = pressure;
     else
-      pressureLookup.erase(key);
+      pressureLookup.erase(index);
     solved = false;
   }
 
@@ -333,8 +338,7 @@ class OxidationDiffusion final : public VelocityField<T>,
 
     concentrationCache_.clear();
     for (const auto &node : nodes)
-      concentrationCache_[detail::gridIndexHash<D>(node.index)] =
-          node.concentration;
+      concentrationCache_[node.index] = node.concentration;
 
     maxScalarVelocity_ = 0.;
     ConstSparseIterator reactionIt(reactionInterface->getDomain());
@@ -432,11 +436,11 @@ class OxidationDiffusion final : public VelocityField<T>,
     return true;
   }
 
-  const std::unordered_map<std::size_t, T> &getConcentrationCache() const {
+  const IndexCacheMap &getConcentrationCache() const {
     return concentrationCache_;
   }
 
-  void setConcentrationCache(std::unordered_map<std::size_t, T> cache) {
+  void setConcentrationCache(IndexCacheMap cache) {
     concentrationCache_ = std::move(cache);
   }
 
@@ -493,7 +497,7 @@ class OxidationDiffusion final : public VelocityField<T>,
       if (!it.isDefined())
         continue;
       const IndexType idx = it.getStartIndices();
-      const auto pIt = pressureLookup.find(detail::gridIndexHash<D>(idx));
+      const auto pIt = pressureLookup.find(idx);
       const T value = pIt != pressureLookup.end() ? pIt->second : T(0);
       pressures.push_back(std::isfinite(value) ? value : T(0));
     }
@@ -564,8 +568,7 @@ class OxidationDiffusion final : public VelocityField<T>,
           if (!it.isDefined())
             continue;
           const auto ptId = it.getPointId();
-          const std::size_t key =
-              detail::gridIndexHash<D>(it.getStartIndices());
+          const auto key = it.getStartIndices();
           if (cd != nullptr && ptId < static_cast<decltype(ptId)>(cd->size()) &&
               std::isfinite((*cd)[ptId]))
             concentrationCache_[key] = (*cd)[ptId];
@@ -590,8 +593,7 @@ class OxidationDiffusion final : public VelocityField<T>,
         const std::size_t id = nodes.size();
         nodeLookupFlat[linearIndex(index)] = id;
         T seedConc = parameters.equilibriumConcentration;
-        auto cacheIt =
-            concentrationCache_.find(detail::gridIndexHash<D>(index));
+        auto cacheIt = concentrationCache_.find(index);
         if (cacheIt != concentrationCache_.end())
           seedConc = cacheIt->second;
         if (!std::isfinite(seedConc))
@@ -1395,8 +1397,7 @@ class OxidationDiffusion final : public VelocityField<T>,
 
     if (parameters.reactionActivationVolume != T(0)) {
       T pressure = parameters.referencePressure;
-      const auto foundPressure =
-          pressureLookup.find(detail::gridIndexHash<D>(index));
+      const auto foundPressure = pressureLookup.find(index);
       if (foundPressure != pressureLookup.end())
         pressure = foundPressure->second;
       if (!std::isfinite(pressure))
@@ -1426,7 +1427,7 @@ class OxidationDiffusion final : public VelocityField<T>,
       return parameters.diffusionCoefficient;
 
     T pressure = parameters.referencePressure;
-    const auto found = pressureLookup.find(detail::gridIndexHash<D>(index));
+    const auto found = pressureLookup.find(index);
     if (found != pressureLookup.end())
       pressure = found->second;
     if (!std::isfinite(pressure))
diff --git a/include/viennals/lsOxidationMask.hpp b/include/viennals/lsOxidationMask.hpp
index 4f4fb1b3..c28d9225 100644
--- a/include/viennals/lsOxidationMask.hpp
+++ b/include/viennals/lsOxidationMask.hpp
@@ -174,7 +174,7 @@ class OxidationMaskBending final : public VelocityField<T>,
   std::vector<Vec3D<T>> contactFaceVelocity_; // Dirichlet velocity at contact
   std::vector<Vec3D<T>> contactFaceTraction_; // Oxide traction on mask face
   std::vector<T> contactFaceDistance_;        // Node-to-interface distance
-  std::unordered_map<std::size_t, T> ambientPhiCache_;
+  std::unordered_map<IndexType, T, typename IndexType::hash> ambientPhiCache_;
 
   // Cached multigrid hierarchy.  The stiffness matrix depends on the node
   // geometry and the contact face classification (active/inactive) but NOT on
@@ -408,7 +408,7 @@ class OxidationMaskBending final : public VelocityField<T>,
       const T dt = parameters.stressTimeStep;
       if (dt <= T(0))
         return {T(0), T(0), T(0)};
-      return detail::vecScaled(getVelocity(index), T(1) / dt);
+      return getVelocity(index) / dt;
     }
 
     return getVelocity(index);
@@ -802,7 +802,7 @@ class OxidationMaskBending final : public VelocityField<T>,
          !it.isFinished(); ++it) {
       if (!it.isDefined())
         continue;
-      ambientPhiCache_[detail::gridIndexHash<D>(it.getStartIndices())] =
+      ambientPhiCache_[it.getStartIndices()] =
           static_cast<T>(ambientSign) * it.getValue();
     }
   }
@@ -946,7 +946,7 @@ class OxidationMaskBending final : public VelocityField<T>,
               // Scale by dt: elastic solver uses dt_ref=1hr so the Dirichlet BC
               // must be the physical displacement (v_oxide × dt), not the
               // velocity.
-              oxVel = detail::vecScaled(oxVel, parameters.stressTimeStep);
+              oxVel = oxVel * parameters.stressTimeStep;
             }
             contactFaceVelocity_[faceIdx * n + id] = oxVel;
           }
@@ -958,7 +958,7 @@ class OxidationMaskBending final : public VelocityField<T>,
 
   std::size_t contactFaceKey(const IndexType &index, unsigned direction,
                              int offset) const {
-    std::size_t seed = detail::gridIndexHash<D>(index);
+    std::size_t seed = typename IndexType::hash{}(index);
     seed ^= std::hash<unsigned>{}(direction) +
             std::size_t(0x9e3779b97f4a7c15ULL) + (seed << 6) + (seed >> 2);
     seed ^= std::hash<int>{}(offset) + std::size_t(0x9e3779b97f4a7c15ULL) +
@@ -1056,12 +1056,10 @@ class OxidationMaskBending final : public VelocityField<T>,
     const std::size_t nn = nodes.size();
     if (contactFaceActive_[faceIdx * nn + nodeId]) {
       if (usesKinematicContactBoundary())
-        return detail::vecSubtract(
-            detail::vecScaled(contactFaceVelocity_[faceIdx * nn + nodeId],
-                              T(2)),
-            Vec3D<T>{static_cast<T>(velocity[nodeId][0]),
-                     static_cast<T>(velocity[nodeId][1]),
-                     static_cast<T>(velocity[nodeId][2])});
+        return contactFaceVelocity_[faceIdx * nn + nodeId] * T(2) -
+               Vec3D<T>{static_cast<T>(velocity[nodeId][0]),
+                        static_cast<T>(velocity[nodeId][1]),
+                        static_cast<T>(velocity[nodeId][2])};
 
       return stressBoundaryGhost(velocity, nodeId, direction, offset,
                                  contactFaceDistance_[faceIdx * nn + nodeId],
@@ -1115,16 +1113,16 @@ class OxidationMaskBending final : public VelocityField<T>,
     Vec3D<T> laplaceAverage{T(0), T(0), T(0)};
     for (unsigned direction = 0; direction < D; ++direction)
       for (int offset : {-1, 1})
-        detail::vecAddTo(laplaceAverage,
-                         neighborVelocity(v, nodeId, direction, offset));
+        laplaceAverage =
+            laplaceAverage + neighborVelocity(v, nodeId, direction, offset);
 
     const T count = static_cast<T>(2 * D);
-    const Vec3D<T> lapAvg = detail::vecScaled(laplaceAverage, T(1) / count);
-    const Vec3D<T> gradDivCorr = detail::vecScaled(
-        divergenceGradient(v, nodes[nodeId].index),
-        gradDivWeight * gridDelta * gridDelta / (T(2) * static_cast<T>(D)));
+    const Vec3D<T> lapAvg = laplaceAverage / count;
+    const Vec3D<T> gradDivCorr =
+        divergenceGradient(v, nodes[nodeId].index) *
+        (gradDivWeight * gridDelta * gridDelta / (T(2) * static_cast<T>(D)));
 
-    return detail::vecAdd(lapAvg, gradDivCorr);
+    return lapAvg + gradDivCorr;
   }
 
   // (Av)[i] = v[i] - F(v)[i] + b[i], stored as SolverT.
@@ -1178,18 +1176,18 @@ class OxidationMaskBending final : public VelocityField<T>,
       coarse.children.reserve(previous.indices.size() / 2 + 1);
       coarse.fineToCoarse.assign(previous.indices.size(), mgNoNode);
 
-      std::unordered_map<std::size_t, std::size_t> coarseLookup;
+      std::unordered_map<IndexType, std::size_t, typename IndexType::hash>
+          coarseLookup;
       coarseLookup.reserve(previous.indices.size());
       for (std::size_t fineId = 0; fineId < previous.indices.size(); ++fineId) {
         const IndexType coarseIndex = coarsenIndex(previous.indices[fineId]);
-        const std::size_t key = detail::gridIndexHash<D>(coarseIndex);
-        auto found = coarseLookup.find(key);
+        auto found = coarseLookup.find(coarseIndex);
         if (found == coarseLookup.end()) {
           const std::size_t coarseId = coarse.indices.size();
-          coarseLookup.emplace(key, coarseId);
+          coarseLookup.emplace(coarseIndex, coarseId);
           coarse.indices.push_back(coarseIndex);
           coarse.children.emplace_back();
-          found = coarseLookup.find(key);
+          found = coarseLookup.find(coarseIndex);
         }
 
         const std::size_t coarseId = found->second;
@@ -2075,7 +2073,7 @@ class OxidationMaskBending final : public VelocityField<T>,
   bool isInsideOxide(const IndexType &index) const {
     if (ambientInterface == nullptr)
       return false;
-    const auto it = ambientPhiCache_.find(detail::gridIndexHash<D>(index));
+    const auto it = ambientPhiCache_.find(index);
     return it != ambientPhiCache_.end() && it->second >= T(0);
   }
 
@@ -2230,7 +2228,7 @@ class OxidationConstrainedAmbient final : public VelocityField<T> {
   SmartPointer<Domain<T, D>> maskInterface = nullptr;
   SmartPointer<Domain<T, D>> ambientInterface = nullptr;
   int maskSign = 1;
-  std::unordered_map<std::size_t, T> maskPhiCache_;
+  std::unordered_map<IndexType, T, typename IndexType::hash> maskPhiCache_;
   T maskGridDelta_ = 1.;
   std::array<T, D> maxVelocity_{};
 
@@ -2344,7 +2342,7 @@ class OxidationConstrainedAmbient final : public VelocityField<T> {
     IndexType index;
     for (unsigned i = 0; i < D; ++i)
       index[i] = std::llround(coordinate[i] / maskGridDelta_);
-    const auto it = maskPhiCache_.find(detail::gridIndexHash<D>(index));
+    const auto it = maskPhiCache_.find(index);
     return (it != maskPhiCache_.end()) ? it->second
                                        : std::numeric_limits<T>::lowest();
   }
@@ -2358,7 +2356,7 @@ class OxidationConstrainedAmbient final : public VelocityField<T> {
          ++it) {
       if (!it.isDefined())
         continue;
-      const auto key = detail::gridIndexHash<D>(it.getStartIndices());
+      const auto key = it.getStartIndices();
       maskPhiCache_[key] = static_cast<T>(maskSign) * it.getValue();
     }
   }
diff --git a/include/viennals/lsOxidationSolverBase.hpp b/include/viennals/lsOxidationSolverBase.hpp
index 5554e5f5..22a7c99a 100644
--- a/include/viennals/lsOxidationSolverBase.hpp
+++ b/include/viennals/lsOxidationSolverBase.hpp
@@ -19,53 +19,6 @@ using namespace viennacore;
 
 namespace detail {
 
-template <int D>
-inline std::size_t gridIndexHash(const viennahrle::Index<D> &index) {
-  std::size_t seed = 0;
-  for (unsigned i = 0; i < static_cast<unsigned>(D); ++i) {
-    seed ^= std::hash<long long>{}(static_cast<long long>(index[i])) +
-            std::size_t(0x9e3779b97f4a7c15ULL) + (seed << 6) + (seed >> 2);
-  }
-  return seed;
-}
-
-// Hash functor for viennahrle::Index<D> — safe to use as unordered_map key
-// because collisions are resolved by operator== on the full index.
-template <int D> struct IndexTypeHasher {
-  std::size_t operator()(const viennahrle::Index<D> &idx) const {
-    return gridIndexHash<D>(idx);
-  }
-};
-
-template <class T> inline Vec3D<T> vecScaled(const Vec3D<T> &source, T factor) {
-  Vec3D<T> result{0., 0., 0.};
-  for (unsigned i = 0; i < 3; ++i)
-    result[i] = source[i] * factor;
-  return result;
-}
-
-template <class T>
-inline Vec3D<T> vecAdd(const Vec3D<T> &a, const Vec3D<T> &b) {
-  Vec3D<T> result{0., 0., 0.};
-  for (unsigned i = 0; i < 3; ++i)
-    result[i] = a[i] + b[i];
-  return result;
-}
-
-template <class T>
-inline Vec3D<T> vecSubtract(const Vec3D<T> &a, const Vec3D<T> &b) {
-  Vec3D<T> result{0., 0., 0.};
-  for (unsigned i = 0; i < 3; ++i)
-    result[i] = a[i] - b[i];
-  return result;
-}
-
-template <class T>
-inline void vecAddTo(Vec3D<T> &target, const Vec3D<T> &source) {
-  for (unsigned i = 0; i < 3; ++i)
-    target[i] += source[i];
-}
-
 /// Clamp HRLE far-field sentinels (±DBL_MAX) to ±1 before differencing to
 /// prevent DBL_MAX² overflow that silently returns the zero vector.
 template <class T> inline T clampLevelSetPhi(T v) {
diff --git a/include/viennals/lsPrune.hpp b/include/viennals/lsPrune.hpp
index 75486ac6..99af7c4d 100644
--- a/include/viennals/lsPrune.hpp
+++ b/include/viennals/lsPrune.hpp
@@ -97,12 +97,8 @@ template <class T, int D> class Prune {
     if (updateData)
       newDataSourceIds.resize(newDomain.getNumberOfSegments());
 
-#pragma omp parallel num_threads(newDomain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < newDomain.getNumberOfSegments(); ++p) {
 
       auto &domainSegment = newDomain.getDomainSegment(p);
 
diff --git a/include/viennals/lsReduce.hpp b/include/viennals/lsReduce.hpp
index 0c0917d0..976953a0 100644
--- a/include/viennals/lsReduce.hpp
+++ b/include/viennals/lsReduce.hpp
@@ -82,12 +82,8 @@ template <class T, int D> class Reduce {
     if (updateData)
       newDataSourceIds.resize(newDomain.getNumberOfSegments());
 
-#pragma omp parallel num_threads(newDomain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < newDomain.getNumberOfSegments(); ++p) {
 
       auto &domainSegment = newDomain.getDomainSegment(p);
 
diff --git a/include/viennals/lsRemoveStrayPoints.hpp b/include/viennals/lsRemoveStrayPoints.hpp
index a09f7f2f..b22d1cd4 100644
--- a/include/viennals/lsRemoveStrayPoints.hpp
+++ b/include/viennals/lsRemoveStrayPoints.hpp
@@ -73,12 +73,8 @@ template <class T, int D> class RemoveStrayPoints {
     std::vector<typename Domain<T, D>::PointValueVectorType> newPoints;
     newPoints.resize(newDomain.getNumberOfSegments());
 
-#pragma omp parallel num_threads(newDomain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < newDomain.getNumberOfSegments(); ++p) {
 
       auto &domainSegment = newDomain.getDomainSegment(p);
 
diff --git a/include/viennals/lsVersion.hpp b/include/viennals/lsVersion.hpp
index d19776a4..5e930de7 100644
--- a/include/viennals/lsVersion.hpp
+++ b/include/viennals/lsVersion.hpp
@@ -5,10 +5,10 @@
 namespace viennals {
 
 // Version information generated by CMake
-inline constexpr const char *version = "5.8.5";
+inline constexpr const char *version = "5.9.0";
 inline constexpr int versionMajor = 5;
-inline constexpr int versionMinor = 8;
-inline constexpr int versionPatch = 5;
+inline constexpr int versionMinor = 9;
+inline constexpr int versionPatch = 0;
 
 // Utility functions for version comparison
 inline constexpr uint32_t versionAsInteger() {
diff --git a/pyproject.toml b/pyproject.toml
index d3336c7d..1bda7252 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,12 +1,12 @@
 [build-system]
 requires = [
     "scikit-build-core>=0.12.2",
-    "pybind11"
+    "pybind11>=3.1.0"
 ]
 build-backend = "scikit_build_core.build"
 
 [project]
-version = "5.8.5"
+version = "5.9.0"
 name = "ViennaLS"
 readme = "README.md"
 license = {file = "LICENSE"}
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index c67869f4..387d2625 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -42,7 +42,7 @@ set(PYBIND11_FINDPYTHON ON)
 
 CPMFindPackage(
   NAME pybind11
-  VERSION 3.0.0
+  VERSION 3.1.0
   GIT_REPOSITORY "https://github.com/pybind/pybind11")
 
 # --------------------------------------------------------------------------------------------------------
diff --git a/python/pyWrap.cpp b/python/pyWrap.cpp
index 33b7f2cd..6ac60e3d 100644
--- a/python/pyWrap.cpp
+++ b/python/pyWrap.cpp
@@ -259,8 +259,9 @@ PYBIND11_MODULE(VIENNALS_MODULE_NAME, module) {
       .def("insertNextHexa", &Mesh<T>::insertNextHexa,
            "Insert a hexahedron in the mesh.")
       .def("removeDuplicateNodes", &Mesh<T>::removeDuplicateNodes,
-           "Remove nodes which occur twice in the mesh, and replace their IDs "
-           "in the mesh elements.")
+           "Remove exactly equal nodes and remap mesh elements, preserving "
+           "first-occurrence order and the first node's point data. Nodes "
+           "containing NaNs remain distinct. Cell data is unchanged.")
       .def("append", &Mesh<T>::append, "Append another mesh to this mesh.")
       .def("print", &Mesh<T>::print, "Print basic information about the mesh.")
       .def("clear", &Mesh<T>::clear, "Clear all data in the mesh.");
diff --git a/python/viennals/__init__.pyi b/python/viennals/__init__.pyi
index 6df3496d..ce457f3b 100644
--- a/python/viennals/__init__.pyi
+++ b/python/viennals/__init__.pyi
@@ -100,7 +100,7 @@ from viennals.d2 import hrleGrid
 from . import _core
 from . import d2
 from . import d3
-__all__: list[str] = ['Advect', 'BooleanOperation', 'BooleanOperationEnum', 'BoundaryConditionEnum', 'Box', 'BoxDistribution', 'CalculateCurvatures', 'CalculateNormalVectors', 'CalculateVisibilities', 'Check', 'CompareArea', 'CompareChamfer', 'CompareCriticalDimensions', 'CompareNarrowBand', 'CompareSparseField', 'CompareVolume', 'ConvexHull', 'Cpu', 'CurvatureEnum', 'CustomSphereDistribution', 'Cylinder', 'DetectFeatures', 'Domain', 'Expand', 'Extrude', 'FeatureDetectionEnum', 'FileFormatEnum', 'FinalizeStencilLocalLaxFriedrichs', 'FromMesh', 'FromSurfaceMesh', 'FromVolumeMesh', 'GeometricAdvect', 'GeometricAdvectDistribution', 'Gpu', 'GpuMode', 'GpuPreconditioner', 'ILU0', 'IntegrationSchemeEnum', 'Jacobi', 'LOCOSConservationDiagnostics', 'LogLevel', 'Logger', 'MakeGeometry', 'MarkVoidPoints', 'MaterialMap', 'Mesh', 'NormalCalculationMethodEnum', 'Oxidation', 'OxidationConstrainedAmbient', 'OxidationCouplingParameters', 'OxidationDeformation', 'OxidationDeformationParameters', 'OxidationDiffusion', 'OxidationMaskBending', 'OxidationMaskParameters', 'OxidationModel', 'OxidationParameters', 'OxidationPresets', 'PROXY_DIM', 'Plane', 'PointCloud', 'PointData', 'PrepareStencilLocalLaxFriedrichs', 'Prune', 'ReactionBoundarySample', 'Reader', 'Reduce', 'RemoveStrayPoints', 'Slice', 'SpatialSchemeEnum', 'Sphere', 'SphereDistribution', 'StencilLocalLaxFriedrichsScalar', 'TemporalSchemeEnum', 'ToDiskMesh', 'ToHullMesh', 'ToMesh', 'ToMultiSurfaceMesh', 'ToSurfaceMesh', 'ToVoxelMesh', 'TransformEnum', 'TransformMesh', 'VTKReader', 'VTKRenderWindow', 'VTKWriter', 'VelocityField', 'VoidTopSurfaceEnum', 'WriteVisualizationMesh', 'Writer', 'computeLOCOSOpenWindowConservation', 'd2', 'd3', 'getDimension', 'hrleGrid', 'setDimension', 'setNumThreads', 'version']
+__all__: list[str] = ['Advect', 'Auto', 'BooleanOperation', 'BooleanOperationEnum', 'BoundaryConditionEnum', 'Box', 'BoxDistribution', 'CalculateCurvatures', 'CalculateNormalVectors', 'CalculateVisibilities', 'Check', 'CompareArea', 'CompareChamfer', 'CompareCriticalDimensions', 'CompareNarrowBand', 'CompareSparseField', 'CompareVolume', 'ConvexHull', 'Cpu', 'CurvatureEnum', 'CustomSphereDistribution', 'Cylinder', 'DetectFeatures', 'Domain', 'Expand', 'Extrude', 'FeatureDetectionEnum', 'FileFormatEnum', 'FinalizeStencilLocalLaxFriedrichs', 'FromMesh', 'FromSurfaceMesh', 'FromVolumeMesh', 'GeometricAdvect', 'GeometricAdvectDistribution', 'Gpu', 'GpuMode', 'GpuPreconditioner', 'ILU0', 'IntegrationSchemeEnum', 'Jacobi', 'LOCOSConservationDiagnostics', 'LogLevel', 'Logger', 'MakeGeometry', 'MarkVoidPoints', 'MaterialMap', 'Mesh', 'NormalCalculationMethodEnum', 'Oxidation', 'OxidationConstrainedAmbient', 'OxidationCouplingParameters', 'OxidationDeformation', 'OxidationDeformationParameters', 'OxidationDiffusion', 'OxidationMaskBending', 'OxidationMaskParameters', 'OxidationModel', 'OxidationParameters', 'OxidationPresets', 'PROXY_DIM', 'Plane', 'PointCloud', 'PointData', 'PrepareStencilLocalLaxFriedrichs', 'Prune', 'ReactionBoundarySample', 'Reader', 'Reduce', 'RemoveStrayPoints', 'Slice', 'SpatialSchemeEnum', 'Sphere', 'SphereDistribution', 'StencilLocalLaxFriedrichsScalar', 'TemporalSchemeEnum', 'ToDiskMesh', 'ToHullMesh', 'ToMesh', 'ToMultiSurfaceMesh', 'ToSurfaceMesh', 'ToVoxelMesh', 'TransformEnum', 'TransformMesh', 'VTKReader', 'VTKRenderWindow', 'VTKWriter', 'VelocityField', 'VoidTopSurfaceEnum', 'WriteVisualizationMesh', 'Writer', 'computeLOCOSOpenWindowConservation', 'd2', 'd3', 'getDimension', 'hrleGrid', 'setDimension', 'setNumThreads', 'version']
 def __dir__():
     ...
 def __getattr__(name):
@@ -127,13 +127,14 @@ def setDimension(d: int):
             Dimension of the simulation (2 or 3).
         
     """
+Auto: _core.GpuMode  # value = <GpuMode.Auto: 2>
 Cpu: _core.GpuMode  # value = <GpuMode.Cpu: 0>
 Gpu: _core.GpuMode  # value = <GpuMode.Gpu: 1>
 ILU0: _core.GpuPreconditioner  # value = <GpuPreconditioner.ILU0: 1>
 Jacobi: _core.GpuPreconditioner  # value = <GpuPreconditioner.Jacobi: 0>
 PROXY_DIM: int = 2
 _SHARED_OXIDATION_TYPES: tuple = ('OxidationParameters', 'OxidationPresets', 'OxidationDeformationParameters', 'OxidationMaskParameters', 'OxidationCouplingParameters')
-__version__: str = '5.8.5'
+__version__: str = '5.9.0'
 _name: str = 'OxidationCouplingParameters'
-version: str = '5.8.5'
+version: str = '5.9.0'
 _C = _core
diff --git a/python/viennals/_core.pyi b/python/viennals/_core.pyi
index f86251b4..16695138 100644
--- a/python/viennals/_core.pyi
+++ b/python/viennals/_core.pyi
@@ -9,7 +9,7 @@ from viennals import d2
 import viennals.d2
 from viennals import d3
 import viennals.d3
-__all__: list[str] = ['BooleanOperationEnum', 'BoundaryConditionEnum', 'Cpu', 'CurvatureEnum', 'Extrude', 'FeatureDetectionEnum', 'FileFormatEnum', 'Gpu', 'GpuMode', 'GpuPreconditioner', 'ILU0', 'IntegrationSchemeEnum', 'Jacobi', 'LOCOSConservationDiagnostics', 'LogLevel', 'Logger', 'MaterialMap', 'Mesh', 'NormalCalculationMethodEnum', 'OxidationCouplingParameters', 'OxidationDeformationParameters', 'OxidationMaskParameters', 'OxidationParameters', 'OxidationPresets', 'PointData', 'Slice', 'SpatialSchemeEnum', 'TemporalSchemeEnum', 'TransformEnum', 'TransformMesh', 'VTKReader', 'VTKRenderWindow', 'VTKWriter', 'VelocityField', 'VoidTopSurfaceEnum', 'd2', 'd3', 'setNumThreads', 'version']
+__all__: list[str] = ['Auto', 'BooleanOperationEnum', 'BoundaryConditionEnum', 'Cpu', 'CurvatureEnum', 'Extrude', 'FeatureDetectionEnum', 'FileFormatEnum', 'Gpu', 'GpuMode', 'GpuPreconditioner', 'ILU0', 'IntegrationSchemeEnum', 'Jacobi', 'LOCOSConservationDiagnostics', 'LogLevel', 'Logger', 'MaterialMap', 'Mesh', 'NormalCalculationMethodEnum', 'OxidationCouplingParameters', 'OxidationDeformationParameters', 'OxidationMaskParameters', 'OxidationParameters', 'OxidationPresets', 'PointData', 'Slice', 'SpatialSchemeEnum', 'TemporalSchemeEnum', 'TransformEnum', 'TransformMesh', 'VTKReader', 'VTKRenderWindow', 'VTKWriter', 'VelocityField', 'VoidTopSurfaceEnum', 'd2', 'd3', 'setNumThreads', 'version']
 class BooleanOperationEnum(enum.IntEnum):
     INTERSECT: typing.ClassVar[BooleanOperationEnum]  # value = <BooleanOperationEnum.INTERSECT: 0>
     INVERT: typing.ClassVar[BooleanOperationEnum]  # value = <BooleanOperationEnum.INVERT: 3>
@@ -111,10 +111,17 @@ class GpuMode:
       Cpu
     
       Gpu
+    
+      Auto
     """
+    Auto: typing.ClassVar[GpuMode]  # value = <GpuMode.Auto: 2>
     Cpu: typing.ClassVar[GpuMode]  # value = <GpuMode.Cpu: 0>
     Gpu: typing.ClassVar[GpuMode]  # value = <GpuMode.Gpu: 1>
-    __members__: typing.ClassVar[dict[str, GpuMode]]  # value = {'Cpu': <GpuMode.Cpu: 0>, 'Gpu': <GpuMode.Gpu: 1>}
+    __members__: typing.ClassVar[dict[str, GpuMode]]  # value = {'Cpu': <GpuMode.Cpu: 0>, 'Gpu': <GpuMode.Gpu: 1>, 'Auto': <GpuMode.Auto: 2>}
+    @typing.overload
+    def __eq__(self, other: GpuMode) -> bool:
+        ...
+    @typing.overload
     def __eq__(self, other: typing.Any) -> bool:
         ...
     def __getstate__(self) -> int:
@@ -127,6 +134,10 @@ class GpuMode:
         ...
     def __int__(self) -> int:
         ...
+    @typing.overload
+    def __ne__(self, other: GpuMode) -> bool:
+        ...
+    @typing.overload
     def __ne__(self, other: typing.Any) -> bool:
         ...
     def __repr__(self) -> str:
@@ -152,6 +163,10 @@ class GpuPreconditioner:
     ILU0: typing.ClassVar[GpuPreconditioner]  # value = <GpuPreconditioner.ILU0: 1>
     Jacobi: typing.ClassVar[GpuPreconditioner]  # value = <GpuPreconditioner.Jacobi: 0>
     __members__: typing.ClassVar[dict[str, GpuPreconditioner]]  # value = {'Jacobi': <GpuPreconditioner.Jacobi: 0>, 'ILU0': <GpuPreconditioner.ILU0: 1>}
+    @typing.overload
+    def __eq__(self, other: GpuPreconditioner) -> bool:
+        ...
+    @typing.overload
     def __eq__(self, other: typing.Any) -> bool:
         ...
     def __getstate__(self) -> int:
@@ -164,6 +179,10 @@ class GpuPreconditioner:
         ...
     def __int__(self) -> int:
         ...
+    @typing.overload
+    def __ne__(self, other: GpuPreconditioner) -> bool:
+        ...
+    @typing.overload
     def __ne__(self, other: typing.Any) -> bool:
         ...
     def __repr__(self) -> str:
@@ -240,6 +259,10 @@ class LogLevel:
     TIMING: typing.ClassVar[LogLevel]  # value = <LogLevel.TIMING: 4>
     WARNING: typing.ClassVar[LogLevel]  # value = <LogLevel.WARNING: 1>
     __members__: typing.ClassVar[dict[str, LogLevel]]  # value = {'ERROR': <LogLevel.ERROR: 0>, 'WARNING': <LogLevel.WARNING: 1>, 'INFO': <LogLevel.INFO: 2>, 'INTERMEDIATE': <LogLevel.INTERMEDIATE: 3>, 'TIMING': <LogLevel.TIMING: 4>, 'DEBUG': <LogLevel.DEBUG: 5>}
+    @typing.overload
+    def __eq__(self, other: LogLevel) -> bool:
+        ...
+    @typing.overload
     def __eq__(self, other: typing.Any) -> bool:
         ...
     def __getstate__(self) -> int:
@@ -252,6 +275,10 @@ class LogLevel:
         ...
     def __int__(self) -> int:
         ...
+    @typing.overload
+    def __ne__(self, other: LogLevel) -> bool:
+        ...
+    @typing.overload
     def __ne__(self, other: typing.Any) -> bool:
         ...
     def __repr__(self) -> str:
@@ -407,7 +434,7 @@ class Mesh:
         """
     def removeDuplicateNodes(self) -> None:
         """
-        Remove nodes which occur twice in the mesh, and replace their IDs in the mesh elements.
+        Remove exactly equal nodes and remap mesh elements, preserving first-occurrence order and the first node's point data. Nodes containing NaNs remain distinct. Cell data is unchanged.
         """
 class NormalCalculationMethodEnum(enum.IntEnum):
     CENTRAL_DIFFERENCES: typing.ClassVar[NormalCalculationMethodEnum]  # value = <NormalCalculationMethodEnum.CENTRAL_DIFFERENCES: 0>
@@ -1059,10 +1086,11 @@ class VoidTopSurfaceEnum(enum.IntEnum):
         """
 def setNumThreads(arg0: typing.SupportsInt | typing.SupportsIndex) -> None:
     ...
+Auto: GpuMode  # value = <GpuMode.Auto: 2>
 Cpu: GpuMode  # value = <GpuMode.Cpu: 0>
 Gpu: GpuMode  # value = <GpuMode.Gpu: 1>
 ILU0: GpuPreconditioner  # value = <GpuPreconditioner.ILU0: 1>
 Jacobi: GpuPreconditioner  # value = <GpuPreconditioner.Jacobi: 0>
-__version__: str = '5.8.5'
-version: str = '5.8.5'
+__version__: str = '5.9.0'
+version: str = '5.9.0'
 IntegrationSchemeEnum = SpatialSchemeEnum
diff --git a/python/viennals/d2.pyi b/python/viennals/d2.pyi
index d7047e01..e942061f 100644
--- a/python/viennals/d2.pyi
+++ b/python/viennals/d2.pyi
@@ -1022,7 +1022,7 @@ class OxidationDiffusion(viennals._core.VelocityField):
         ...
     def getConcentration(self, arg0: typing.Annotated[collections.abc.Sequence[typing.SupportsFloat | typing.SupportsIndex], "FixedSize(3)"]) -> float:
         ...
-    def getConcentrationCache(self) -> dict[int, float]:
+    def getConcentrationCache(self) -> dict[..., float]:
         ...
     def getEffectiveReactionRate(self, arg0: typing.Annotated[collections.abc.Sequence[typing.SupportsFloat | typing.SupportsIndex], "FixedSize(3)"]) -> float:
         ...
@@ -1044,7 +1044,7 @@ class OxidationDiffusion(viennals._core.VelocityField):
         ...
     def setAmbientInterface(self, arg0: Domain) -> None:
         ...
-    def setConcentrationCache(self, arg0: collections.abc.Mapping[typing.SupportsInt | typing.SupportsIndex, typing.SupportsFloat | typing.SupportsIndex]) -> None:
+    def setConcentrationCache(self, arg0: collections.abc.Mapping[..., typing.SupportsFloat | typing.SupportsIndex]) -> None:
         ...
     def setMaskInterface(self, maskInterface: Domain, maskSign: typing.SupportsInt | typing.SupportsIndex = 1) -> None:
         ...
diff --git a/python/viennals/d3.pyi b/python/viennals/d3.pyi
index bcaf1c2c..8346e8c4 100644
--- a/python/viennals/d3.pyi
+++ b/python/viennals/d3.pyi
@@ -1022,7 +1022,7 @@ class OxidationDiffusion(viennals._core.VelocityField):
         ...
     def getConcentration(self, arg0: typing.Annotated[collections.abc.Sequence[typing.SupportsFloat | typing.SupportsIndex], "FixedSize(3)"]) -> float:
         ...
-    def getConcentrationCache(self) -> dict[int, float]:
+    def getConcentrationCache(self) -> dict[..., float]:
         ...
     def getEffectiveReactionRate(self, arg0: typing.Annotated[collections.abc.Sequence[typing.SupportsFloat | typing.SupportsIndex], "FixedSize(3)"]) -> float:
         ...
@@ -1044,7 +1044,7 @@ class OxidationDiffusion(viennals._core.VelocityField):
         ...
     def setAmbientInterface(self, arg0: Domain) -> None:
         ...
-    def setConcentrationCache(self, arg0: collections.abc.Mapping[typing.SupportsInt | typing.SupportsIndex, typing.SupportsFloat | typing.SupportsIndex]) -> None:
+    def setConcentrationCache(self, arg0: collections.abc.Mapping[..., typing.SupportsFloat | typing.SupportsIndex]) -> None:
         ...
     def setMaskInterface(self, maskInterface: Domain, maskSign: typing.SupportsInt | typing.SupportsIndex = 1) -> None:
         ...
diff --git a/tests/RemoveDuplicateNodes/CMakeLists.txt b/tests/RemoveDuplicateNodes/CMakeLists.txt
new file mode 100644
index 00000000..9a6bfd95
--- /dev/null
+++ b/tests/RemoveDuplicateNodes/CMakeLists.txt
@@ -0,0 +1,7 @@
+project(RemoveDuplicateNodes LANGUAGES CXX)
+
+add_executable(${PROJECT_NAME} "${PROJECT_NAME}.cpp")
+target_link_libraries(${PROJECT_NAME} PRIVATE ViennaLS)
+
+add_dependencies(ViennaLS_Tests ${PROJECT_NAME})
+add_test(NAME ${PROJECT_NAME} COMMAND $<TARGET_FILE:${PROJECT_NAME}>)
diff --git a/tests/RemoveDuplicateNodes/RemoveDuplicateNodes.cpp b/tests/RemoveDuplicateNodes/RemoveDuplicateNodes.cpp
new file mode 100644
index 00000000..d5071a34
--- /dev/null
+++ b/tests/RemoveDuplicateNodes/RemoveDuplicateNodes.cpp
@@ -0,0 +1,150 @@
+#include <lsMesh.hpp>
+#include <vcTestAsserts.hpp>
+
+#include <cmath>
+#include <limits>
+
+template <class T> void testEmptyAndUnique() {
+  viennals::Mesh<T> mesh;
+  mesh.removeDuplicateNodes();
+  VC_TEST_ASSERT(mesh.nodes.empty());
+
+  mesh.nodes = {{1, 2, 3}};
+  mesh.pointData.insertNextScalarData({T(7)}, "Values");
+  mesh.removeDuplicateNodes();
+  VC_TEST_ASSERT(mesh.nodes.size() == 1);
+  VC_TEST_ASSERT(mesh.pointData.getScalarData("Values")->at(0) == T(7));
+
+  mesh.nodes.push_back({4, 5, 6});
+  mesh.pointData.getScalarData("Values")->push_back(T(8));
+  mesh.lines = {{1, 0}};
+  const auto originalNodes = mesh.nodes;
+  const auto originalLines = mesh.lines;
+  const auto originalValues = *mesh.pointData.getScalarData("Values");
+  mesh.removeDuplicateNodes();
+  VC_TEST_ASSERT(mesh.nodes == originalNodes);
+  VC_TEST_ASSERT(mesh.lines == originalLines);
+  VC_TEST_ASSERT(*mesh.pointData.getScalarData("Values") == originalValues);
+}
+
+template <class T> void testConnectivityAndData() {
+  viennals::Mesh<T> mesh;
+  mesh.nodes = {{0, 0, 0}, {1, 0, 0}, {0, 0, 0}, {0, 1, 0}, {1, 0, 0},
+                {0, 0, 1}, {1, 1, 0}, {1, 0, 1}, {0, 1, 1}, {1, 1, 1}};
+  mesh.vertices = {{2}, {9}};
+  mesh.lines = {{4, 3}, {0, 2}};
+  mesh.triangles = {{2, 3, 9}};
+  mesh.tetras = {{0, 4, 5, 9}};
+  mesh.hexas = {{0, 1, 3, 5, 6, 7, 8, 9}};
+  mesh.minimumExtent = {0, 0, 0};
+  mesh.maximumExtent = {1, 1, 1};
+  mesh.pointData.insertNextScalarData({10, 20, 99, 30, 98, 40, 50, 60, 70, 80},
+                                      "Values");
+  std::vector<viennals::Vec3D<T>> vectors;
+  for (const T value : *mesh.pointData.getScalarData("Values"))
+    vectors.push_back({value, 0, -value});
+  mesh.pointData.insertNextVectorData(vectors, "Vectors");
+  const std::vector<T> cellValues = {1, 2, 3, 4, 5, 6, 7};
+  mesh.cellData.insertNextScalarData(cellValues, "Materials");
+
+  viennals::Mesh<T> expected;
+  expected.nodes = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1},
+                    {1, 1, 0}, {1, 0, 1}, {0, 1, 1}, {1, 1, 1}};
+  expected.vertices = {{0}, {7}};
+  expected.lines = {{1, 2}, {0, 0}};
+  expected.triangles = {{0, 2, 7}};
+  expected.tetras = {{0, 1, 3, 7}};
+  expected.hexas = {{0, 1, 2, 3, 4, 5, 6, 7}};
+  const std::vector<T> expectedValues = {10, 20, 30, 40, 50, 60, 70, 80};
+  std::vector<viennals::Vec3D<T>> expectedVectors;
+  for (const T value : expectedValues)
+    expectedVectors.push_back({value, 0, -value});
+
+  // Verify both the first application and idempotence.
+  for (int pass = 0; pass < 2; ++pass) {
+    mesh.removeDuplicateNodes();
+    VC_TEST_ASSERT(mesh.nodes == expected.nodes);
+    VC_TEST_ASSERT(mesh.vertices == expected.vertices);
+    VC_TEST_ASSERT(mesh.lines == expected.lines);
+    VC_TEST_ASSERT(mesh.triangles == expected.triangles);
+    VC_TEST_ASSERT(mesh.tetras == expected.tetras);
+    VC_TEST_ASSERT(mesh.hexas == expected.hexas);
+    VC_TEST_ASSERT(*mesh.pointData.getScalarData("Values") == expectedValues);
+    VC_TEST_ASSERT(*mesh.pointData.getVectorData("Vectors") == expectedVectors);
+    VC_TEST_ASSERT(*mesh.cellData.getScalarData("Materials") == cellValues);
+    VC_TEST_ASSERT(mesh.minimumExtent == expected.nodes.front());
+    VC_TEST_ASSERT(mesh.maximumExtent == expected.nodes.back());
+  }
+}
+
+template <class T> void testFloatingPointEquality() {
+  const T inf = std::numeric_limits<T>::infinity();
+  const T nan = std::numeric_limits<T>::quiet_NaN();
+  viennals::Mesh<T> mesh;
+  mesh.nodes = {{-T(0), 0, 0},
+                {T(0), 0, 0},
+                {inf, 0, 0},
+                {inf, 0, 0},
+                {-inf, 0, 0},
+                {nan, 0, 0},
+                {nan, 0, 0},
+                {0, nan, 0},
+                {0, nan, 0},
+                {0, 0, nan},
+                {0, 0, nan},
+                {1, 0, 0},
+                {std::nextafter(T(1), T(2)), 0, 0}};
+  for (unsigned i = 0; i < mesh.nodes.size(); ++i)
+    mesh.vertices.push_back({i});
+  const std::vector<std::array<unsigned, 1>> expectedVertices = {
+      {0}, {0}, {1}, {1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, {10}};
+  mesh.removeDuplicateNodes();
+  VC_TEST_ASSERT(mesh.nodes.size() == 11);
+  VC_TEST_ASSERT(mesh.vertices == expectedVertices);
+  VC_TEST_ASSERT(std::signbit(mesh.nodes[0][0]));
+  VC_TEST_ASSERT(mesh.nodes[1][0] == inf);
+  VC_TEST_ASSERT(mesh.nodes[2][0] == -inf);
+  VC_TEST_ASSERT(std::isnan(mesh.nodes[3][0]));
+  VC_TEST_ASSERT(std::isnan(mesh.nodes[5][1]));
+  VC_TEST_ASSERT(std::isnan(mesh.nodes[7][2]));
+  VC_TEST_ASSERT(mesh.nodes[9][0] != mesh.nodes[10][0]);
+}
+
+template <class T> void testAllDuplicatesAndInvalidData() {
+  viennals::Mesh<T> mesh;
+  mesh.nodes.assign(100, {1, 2, 3});
+  mesh.lines = {{99, 50}};
+  // Reject incomplete data without partially modifying the mesh.
+  mesh.pointData.insertNextScalarData({T(7)}, "Values");
+  const auto originalNodes = mesh.nodes;
+  const auto originalLines = mesh.lines;
+  bool rejected = false;
+  try {
+    mesh.removeDuplicateNodes();
+  } catch (const std::invalid_argument &) {
+    rejected = true;
+  }
+  VC_TEST_ASSERT(rejected);
+  VC_TEST_ASSERT(mesh.nodes == originalNodes);
+  VC_TEST_ASSERT(mesh.lines == originalLines);
+  VC_TEST_ASSERT(mesh.pointData.getScalarData("Values")->size() == 1);
+
+  mesh.pointData.getScalarData("Values")->resize(100, T(99));
+  mesh.removeDuplicateNodes();
+  VC_TEST_ASSERT(mesh.nodes.size() == 1);
+  VC_TEST_ASSERT(mesh.lines.front()[0] == 0 && mesh.lines.front()[1] == 0);
+  VC_TEST_ASSERT(mesh.pointData.getScalarData("Values")->size() == 1);
+  VC_TEST_ASSERT(mesh.pointData.getScalarData("Values")->front() == T(7));
+}
+
+template <class T> void runTests() {
+  testEmptyAndUnique<T>();
+  testConnectivityAndData<T>();
+  testFloatingPointEquality<T>();
+  testAllDuplicatesAndInvalidData<T>();
+}
+
+int main() {
+  runTests<float>();
+  runTests<double>();
+}