diff --git a/Wrapping/Python/pyWrap.cpp b/Wrapping/Python/pyWrap.cpp
index 91fa0461..9eb934fa 100644
--- a/Wrapping/Python/pyWrap.cpp
+++ b/Wrapping/Python/pyWrap.cpp
@@ -100,9 +100,10 @@ class PylsGeometricAdvectDistribution
     PYBIND11_OVERLOAD(bool, ClassType, isInside, initial, candidate, eps);
   }
 
-  T getSignedDistance(const vectorType &initial,
-                      const vectorType &candidate) const override {
-    PYBIND11_OVERLOAD_PURE(T, ClassType, getSignedDistance, initial, candidate);
+  T getSignedDistance(const vectorType &initial, const vectorType &candidate,
+                      unsigned long initialPointId) const override {
+    PYBIND11_OVERLOAD_PURE(T, ClassType, getSignedDistance, initial, candidate,
+                           initialPointId);
   }
 
   boundsType getBounds() const override {
@@ -407,7 +408,7 @@ PYBIND11_MODULE(VIENNALS_MODULE_NAME, module) {
           pybind11::detail::pythonbuf buf(fileHandle);
           std::ostream stream(&buf);
           d.print(stream);
-          }, pybind11::arg("stream") = pybind11::module_::import("sys").attr("stdout"));
+          }, pybind11::arg("stream") = pybind11::module::import("sys").attr("stdout"));
 
   // enums
   pybind11::enum_<lsBoundaryConditionEnum<D>>(module, "lsBoundaryConditionEnum")
@@ -605,7 +606,7 @@ PYBIND11_MODULE(VIENNALS_MODULE_NAME, module) {
                           const std::vector<std::vector<T>> &>))
       // methods
       .def("insertNextPoint",
-           (void (lsPointCloud<T, D>::*)(const std::vector<T> &)) &
+           (void(lsPointCloud<T, D>::*)(const std::vector<T> &)) &
                lsPointCloud<T, D>::insertNextPoint);
 
   // lsMakeGeometry
@@ -631,22 +632,22 @@ PYBIND11_MODULE(VIENNALS_MODULE_NAME, module) {
       .def("setLevelSet", &lsMakeGeometry<T, D>::setLevelSet,
            "Set the levelset in which to create the geometry.")
       .def("setGeometry",
-           (void (lsMakeGeometry<T, D>::*)(lsSmartPointer<lsSphere<T, D>>)) &
+           (void(lsMakeGeometry<T, D>::*)(lsSmartPointer<lsSphere<T, D>>)) &
                lsMakeGeometry<T, D>::setGeometry)
       .def("setGeometry",
-           (void (lsMakeGeometry<T, D>::*)(lsSmartPointer<lsPlane<T, D>>)) &
+           (void(lsMakeGeometry<T, D>::*)(lsSmartPointer<lsPlane<T, D>>)) &
                lsMakeGeometry<T, D>::setGeometry)
       .def("setGeometry",
-           (void (lsMakeGeometry<T, D>::*)(lsSmartPointer<lsBox<T, D>>)) &
+           (void(lsMakeGeometry<T, D>::*)(lsSmartPointer<lsBox<T, D>>)) &
+               lsMakeGeometry<T, D>::setGeometry)
+      .def("setGeometry",
+           (void(lsMakeGeometry<T, D>::*)(lsSmartPointer<lsPointCloud<T, D>>)) &
                lsMakeGeometry<T, D>::setGeometry)
-      .def("setGeometry", (void (lsMakeGeometry<T, D>::*)(
-                              lsSmartPointer<lsPointCloud<T, D>>)) &
-                              lsMakeGeometry<T, D>::setGeometry)
       .def("setIgnoreBoundaryConditions",
-           (void (lsMakeGeometry<T, D>::*)(bool)) &
+           (void(lsMakeGeometry<T, D>::*)(bool)) &
                lsMakeGeometry<T, D>::setIgnoreBoundaryConditions)
       .def("setIgnoreBoundaryConditions",
-           (void (lsMakeGeometry<T, D>::*)(std::array<bool, 3>)) &
+           (void(lsMakeGeometry<T, D>::*)(std::array<bool, 3>)) &
                lsMakeGeometry<T, D>::setIgnoreBoundaryConditions)
       .def("apply", &lsMakeGeometry<T, D>::apply, "Generate the geometry.");
 
@@ -692,13 +693,13 @@ PYBIND11_MODULE(VIENNALS_MODULE_NAME, module) {
       .def(pybind11::init(&lsSmartPointer<lsPointData<T>>::New<>))
       // methods
       .def("insertNextScalarData",
-           (void (lsPointData<T>::*)(const lsPointData<T>::ScalarDataType &,
-                                     std::string)) &
+           (void(lsPointData<T>::*)(const lsPointData<T>::ScalarDataType &,
+                                    std::string)) &
                lsPointData<T>::insertNextScalarData,
            pybind11::arg("scalars"), pybind11::arg("label") = "Scalars")
       .def("insertNextVectorData",
-           (void (lsPointData<T>::*)(const lsPointData<T>::VectorDataType &,
-                                     std::string)) &
+           (void(lsPointData<T>::*)(const lsPointData<T>::VectorDataType &,
+                                    std::string)) &
                lsPointData<T>::insertNextVectorData,
            pybind11::arg("vectors"), pybind11::arg("label") = "Vectors")
       .def("getScalarDataSize", &lsPointData<T>::getScalarDataSize)
@@ -728,7 +729,7 @@ PYBIND11_MODULE(VIENNALS_MODULE_NAME, module) {
                lsMesh<T>::getNodes,
            "Get all nodes of the mesh as a list.")
       .def("getNodes",
-           (const std::vector<std::array<double, 3>> &(lsMesh<T>::*)() const) &
+           (std::vector<std::array<double, 3>> & (lsMesh<T>::*)()) &
                lsMesh<T>::getNodes,
            "Get all nodes of the mesh as a list.")
       .def("getVerticies",
diff --git a/include/lsGeometricAdvect.hpp b/include/lsGeometricAdvect.hpp
index 5f0cd103..081af244 100644
--- a/include/lsGeometricAdvect.hpp
+++ b/include/lsGeometricAdvect.hpp
@@ -52,6 +52,15 @@ template <class T, int D> class lsGeometricAdvect {
     ++indices[dim];
   }
 
+  template <class K, class V, template <class...> class MapType, class... Ts>
+  MapType<V, K> inverseTranslator(MapType<K, V, Ts...> &map) {
+    MapType<V, K> inv;
+    std::for_each(map.begin(), map.end(), [&inv](const std::pair<K, V> &p) {
+      inv.insert(std::make_pair(p.second, p.first));
+    });
+    return inv;
+  }
+
 public:
   lsGeometricAdvect() {}
 
@@ -123,7 +132,11 @@ template <class T, int D> class lsGeometricAdvect {
     // Extract the original surface as a point cloud of grid
     // points shifted to the surface (disk mesh)
     auto surfaceMesh = lsSmartPointer<lsMesh<hrleCoordType>>::New();
-    lsToDiskMesh<T, D, hrleCoordType>(levelSet, surfaceMesh).apply();
+    auto pointIdTranslator =
+        lsSmartPointer<typename lsToDiskMesh<T, D>::TranslatorType>::New();
+    lsToDiskMesh<T, D, hrleCoordType>(levelSet, surfaceMesh, pointIdTranslator)
+        .apply();
+    *pointIdTranslator = inverseTranslator(*pointIdTranslator);
 
     // find bounds of distribution
     auto distBounds = dist->getBounds();
@@ -359,10 +372,11 @@ template <class T, int D> class lsGeometricAdvect {
 
         T distance = initialDistance;
 
+        unsigned long currentPointId = 0;
         // now check which surface points contribute to currentIndex
         for (typename SurfaceNodesType::const_iterator surfIt =
                  surfaceNodes.begin();
-             surfIt != surfaceNodes.end(); ++surfIt) {
+             surfIt != surfaceNodes.end(); ++surfIt, ++currentPointId) {
 
           auto &currentNode = *surfIt;
 
@@ -387,8 +401,10 @@ template <class T, int D> class lsGeometricAdvect {
           }
 
           // get filling fraction from distance to dist surface
-          T tmpDistance =
-              dist->getSignedDistance(currentNode, currentCoords) / gridDelta;
+          T tmpDistance = dist->getSignedDistance(
+                              currentNode, currentCoords,
+                              pointIdTranslator->find(currentPointId)->second) /
+                          gridDelta;
 
           // if cell is far within a distribution, set it filled
           if (distIsPositive) {
diff --git a/include/lsGeometricAdvectDistributions.hpp b/include/lsGeometricAdvectDistributions.hpp
index 936d010a..0cf19ba6 100644
--- a/include/lsGeometricAdvectDistributions.hpp
+++ b/include/lsGeometricAdvectDistributions.hpp
@@ -24,9 +24,9 @@ template <class T, int D> class lsGeometricAdvectDistribution {
   /// Returns the signed distance of a point relative to the distributions
   /// center. This is the signed manhatten distance to the nearest surface
   /// point.
-  virtual T
-  getSignedDistance(const std::array<hrleCoordType, 3> &initial,
-                    const std::array<hrleCoordType, 3> &candidate) const = 0;
+  virtual T getSignedDistance(const std::array<hrleCoordType, 3> &initial,
+                              const std::array<hrleCoordType, 3> &candidate,
+                              unsigned long initialPointId) const = 0;
 
   /// Sets bounds to the bounding box of the distribution.
   virtual std::array<hrleCoordType, 6> getBounds() const = 0;
@@ -48,7 +48,7 @@ class lsSphereDistribution : public lsGeometricAdvectDistribution<T, D> {
 
   bool isInside(const std::array<hrleCoordType, 3> &initial,
                 const std::array<hrleCoordType, 3> &candidate,
-                double eps = 0.) const {
+                double eps = 0.) const override {
     hrleCoordType dot = 0.;
     for (unsigned i = 0; i < D; ++i) {
       double tmp = candidate[i] - initial[i];
@@ -62,7 +62,8 @@ class lsSphereDistribution : public lsGeometricAdvectDistribution<T, D> {
   }
 
   T getSignedDistance(const std::array<hrleCoordType, 3> &initial,
-                      const std::array<hrleCoordType, 3> &candidate) const {
+                      const std::array<hrleCoordType, 3> &candidate,
+                      unsigned long /*initialPointId*/) const override {
     T distance = std::numeric_limits<T>::max();
     std::array<hrleCoordType, 3> v{};
     for (unsigned i = 0; i < D; ++i) {
@@ -95,7 +96,7 @@ class lsSphereDistribution : public lsGeometricAdvectDistribution<T, D> {
     }
   }
 
-  std::array<hrleCoordType, 6> getBounds() const {
+  std::array<hrleCoordType, 6> getBounds() const override {
     std::array<hrleCoordType, 6> bounds = {};
     for (unsigned i = 0; i < D; ++i) {
       bounds[2 * i] = -radius;
@@ -128,7 +129,7 @@ class lsBoxDistribution : public lsGeometricAdvectDistribution<T, D> {
 
   bool isInside(const std::array<hrleCoordType, 3> &initial,
                 const std::array<hrleCoordType, 3> &candidate,
-                double eps = 0.) const {
+                double eps = 0.) const override {
     for (unsigned i = 0; i < D; ++i) {
       if (std::abs(candidate[i] - initial[i]) >
           (std::abs(posExtent[i]) + eps)) {
@@ -139,7 +140,8 @@ class lsBoxDistribution : public lsGeometricAdvectDistribution<T, D> {
   }
 
   T getSignedDistance(const std::array<hrleCoordType, 3> &initial,
-                      const std::array<hrleCoordType, 3> &candidate) const {
+                      const std::array<hrleCoordType, 3> &candidate,
+                      unsigned long /*initialPointId*/) const override {
     T distance = std::numeric_limits<T>::lowest();
     for (unsigned i = 0; i < D; ++i) {
       T vector = std::abs(candidate[i] - initial[i]);
@@ -148,7 +150,7 @@ class lsBoxDistribution : public lsGeometricAdvectDistribution<T, D> {
     return (posExtent[0] < 0) ? -distance : distance;
   }
 
-  std::array<hrleCoordType, 6> getBounds() const {
+  std::array<hrleCoordType, 6> getBounds() const override {
     std::array<hrleCoordType, 6> bounds = {};
     for (unsigned i = 0; i < D; ++i) {
       bounds[2 * i] = -posExtent[i];
diff --git a/include/lsSmartPointer.hpp b/include/lsSmartPointer.hpp
index 89bf0af0..d8c56422 100644
--- a/include/lsSmartPointer.hpp
+++ b/include/lsSmartPointer.hpp
@@ -14,11 +14,11 @@ template <class T> class lsSmartPointer : public std::shared_ptr<T> {
   // Make visible all constructors of std::shared_ptr
   // including copy constructors
   template <typename... Args>
-  lsSmartPointer(Args &&...args)
+  lsSmartPointer(Args &&... args)
       : std::shared_ptr<T>(std::forward<Args>(args)...) {}
 
   /// Use this function to create new objects when using ViennaLS
-  template <typename... TArgs> static lsSmartPointer New(TArgs &&...targs) {
+  template <typename... TArgs> static lsSmartPointer New(TArgs &&... targs) {
     return lsSmartPointer(std::make_shared<T>(std::forward<TArgs>(targs)...));
   }
 };
diff --git a/include/lsToDiskMesh.hpp b/include/lsToDiskMesh.hpp
index 67dcb130..8e2f8810 100644
--- a/include/lsToDiskMesh.hpp
+++ b/include/lsToDiskMesh.hpp
@@ -19,11 +19,14 @@
 /// This allows for a simple setup of disks for ray tracing.
 template <class T, int D, class N = T> class lsToDiskMesh {
   typedef typename lsDomain<T, D>::DomainType hrleDomainType;
-  typedef std::unordered_map<unsigned long, unsigned long> translatorType;
 
+public:
+  using TranslatorType = std::unordered_map<unsigned long, unsigned long>;
+
+private:
   std::vector<lsSmartPointer<lsDomain<T, D>>> levelSets;
   lsSmartPointer<lsMesh<N>> mesh = nullptr;
-  lsSmartPointer<translatorType> translator = nullptr;
+  lsSmartPointer<TranslatorType> translator = nullptr;
   T maxValue = 0.5;
   bool buildTranslator = false;
   static constexpr double wrappingLayerEpsilon = 1e-4;
@@ -42,7 +45,7 @@ template <class T, int D, class N = T> class lsToDiskMesh {
 
   lsToDiskMesh(lsSmartPointer<lsDomain<T, D>> passedLevelSet,
                lsSmartPointer<lsMesh<N>> passedMesh,
-               lsSmartPointer<translatorType> passedTranslator,
+               lsSmartPointer<TranslatorType> passedTranslator,
                T passedMaxValue = 0.5)
       : mesh(passedMesh), translator(passedTranslator),
         maxValue(passedMaxValue) {
@@ -61,7 +64,7 @@ template <class T, int D, class N = T> class lsToDiskMesh {
 
   void setMesh(lsSmartPointer<lsMesh<N>> passedMesh) { mesh = passedMesh; }
 
-  void setTranslator(lsSmartPointer<translatorType> passedTranslator) {
+  void setTranslator(lsSmartPointer<TranslatorType> passedTranslator) {
     translator = passedTranslator;
     buildTranslator = true;
   }