Update dependencies, fix element to point data conversion radius

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 project(
   ViennaPS
   LANGUAGES CXX C
-  VERSION 4.1.2)
+  VERSION 4.1.3)
 
 # --------------------------------------------------------------------------------------------------------
 # Library switches
@@ -101,7 +101,7 @@ include("cmake/cpm.cmake")
 
 CPMAddPackage(
   NAME ViennaCore
-  VERSION 1.7.2
+  VERSION 1.7.3
   GIT_REPOSITORY "https://github.com/ViennaTools/ViennaCore"
   EXCLUDE_FROM_ALL ${VIENNAPS_BUILD_PYTHON}
   OPTIONS "VIENNACORE_USE_GPU ${VIENNAPS_USE_GPU}")
@@ -113,7 +113,7 @@ CPMAddPackage(
 
 CPMFindPackage(
   NAME ViennaRay
-  VERSION 3.8.3
+  VERSION 3.8.4
   GIT_REPOSITORY "https://github.com/ViennaTools/ViennaRay"
   EXCLUDE_FROM_ALL ${VIENNAPS_BUILD_PYTHON}
   OPTIONS "VIENNARAY_USE_GPU ${VIENNAPS_USE_GPU}")

README.md

@@ -123,7 +123,7 @@ We recommend using [CPM.cmake](https://github.com/cpm-cmake/CPM.cmake) to consum
 
 * Installation with CPM
   ```cmake
-  CPMAddPackage("gh:viennatools/viennaps@4.1.2")
+  CPMAddPackage("gh:viennatools/viennaps@4.1.3")
   ```
 
 * With a local installation

include/viennaps/models/psFaradayCageEtching.hpp

@@ -287,6 +287,8 @@ class FaradayCageEtching final : public ProcessModelGPU<NumericType, D> {
     this->hasGPU = true;
   }
 
+  ~FaradayCageEtching() override { this->processData.free(); }
+
 private:
   FaradayCageParameters<NumericType> params_;
   std::vector<Material> maskMaterials_;

include/viennaps/psCreateSurfaceMesh.hpp

@@ -184,7 +184,7 @@ template <class LsNT, class MeshNT = LsNT, int D = 3> class CreateSurfaceMesh {
           } else {
             // if node does not exist yet
             // calculate coordinate of new node
-            Vec3D<MeshNT> cc{}; // initialise with zeros
+            Vec3D<MeshNT> cc{0., 0., 0.}; // initialise with zeros
             for (int z = 0; z < D; z++) {
               if (z != dir) {
                 // TODO might not need BitMaskToVector here, just check if z

include/viennaps/psElementToPointData.hpp

@@ -109,7 +109,6 @@ class ElementToPointDataBase {
     const auto numElements = elementKdTree_->getNumberOfPoints();
     const auto normals = diskMesh_->cellData.getVectorData("Normals");
     const auto elementNormals = surfaceMesh_->cellData.getVectorData("Normals");
-    const auto sqrdDist = conversionRadius_ * conversionRadius_;
 
     // retrieve data from device
     std::vector<MeshNT> elementData;
@@ -133,20 +132,22 @@ class ElementToPointDataBase {
     for (unsigned i = 0; i < numPoints; i++) {
 
       // we have to use the squared distance here
-      auto closePoints =
-          elementKdTree_->findNearestWithinRadius(points[i], sqrdDist).value();
+      auto closeElements =
+          elementKdTree_->findNearestWithinRadius(points[i], conversionRadius_)
+              .value();
 
-      std::vector<NumericType> weights(closePoints.size(), 0.);
+      std::vector<double> weights(closeElements.size(), 0.);
 
+      // compute weights based on normal alignment
       unsigned numClosePoints = 0;
-      for (std::size_t n = 0; n < closePoints.size(); ++n) {
-        const auto &p = closePoints[n];
+      for (std::size_t n = 0; n < closeElements.size(); ++n) {
+        const auto &p = closeElements[n];
         assert(p.first < numElements);
 
-        const auto weight =
+        const double weight =
             DotProductNT(normals->at(i), elementNormals->at(p.first));
 
-        if (weight > NumericType(1e-6) && !std::isnan(weight)) {
+        if (weight > 1e-6 && !std::isnan(weight)) {
           weights[n] = weight;
           ++numClosePoints;
         }
@@ -162,22 +163,22 @@ class ElementToPointDataBase {
         NumericType value = NumericType(0);
 
         if (numClosePoints > 0) {
-          auto weightsCopy = weights;
 
-          // Discard outliers if enabled
-          discardOutliers(weightsCopy, weights, closePoints, elementData, j,
-                          numElements, numClosePoints);
+          // Discard outlier values if enabled
+          const auto weightsCopy =
+              discardOutliers(weights, closeElements, elementData, j,
+                              numElements, numClosePoints);
 
           // Compute weighted average
-          const NumericType sum = std::accumulate(
-              weightsCopy.begin(), weightsCopy.end(), NumericType(0));
-
-          if (sum > NumericType(0)) {
-            for (std::size_t k = 0; k < closePoints.size(); ++k) {
-              if (weightsCopy[k] > NumericType(0)) {
-                const unsigned elementIdx =
-                    closePoints[k].first + j * numElements;
-                value += weightsCopy[k] * elementData[elementIdx];
+          const double sum =
+              std::accumulate(weightsCopy.cbegin(), weightsCopy.cend(), 0.0);
+
+          if (sum > 1e-6) {
+            for (std::size_t k = 0; k < closeElements.size(); ++k) {
+              if (weightsCopy[k] > 0.0) {
+                const unsigned flattIdx =
+                    closeElements[k].first + j * numElements;
+                value += weightsCopy[k] * elementData[flattIdx];
               }
             }
             value /= sum;
@@ -188,6 +189,7 @@ class ElementToPointDataBase {
           }
         } else {
           // Fallback to nearest point
+          assert(numClosePoints == 0);
           value = elementData[nearestIdx + j * numElements];
         }
 
@@ -211,15 +213,15 @@ class ElementToPointDataBase {
           minIdx1(-1), minIdx2(-1), maxIdx1(-1), maxIdx2(-1) {}
   };
 
-  MinMaxInfo findMinMaxValues(
-      const std::vector<NumericType> &weights,
+  static MinMaxInfo findMinMaxValues(
+      const std::vector<double> &weights,
       const std::vector<std::pair<std::size_t, NumericType>> &closePoints,
       const std::vector<MeshNT> &elementData, unsigned j,
       unsigned numElements) {
     MinMaxInfo info;
 
     for (std::size_t k = 0; k < closePoints.size(); ++k) {
-      if (weights[k] != NumericType(0)) {
+      if (weights[k] > 0.0) {
         const unsigned elementIdx = closePoints[k].first + j * numElements;
         const auto value = elementData[elementIdx];
 
@@ -250,33 +252,37 @@ class ElementToPointDataBase {
     return info;
   }
 
-  void discardOutliers(
-      std::vector<NumericType> &weightsCopy,
-      const std::vector<NumericType> &weights,
+  static auto discardOutliers(
+      const std::vector<double> &weights,
       const std::vector<std::pair<std::size_t, NumericType>> &closePoints,
       const std::vector<MeshNT> &elementData, unsigned j, unsigned numElements,
       unsigned numClosePoints) {
 
+    // copy weights to modify
+    auto weightsCopy = weights;
+
     if (discard4 && numClosePoints > 4) {
       const auto info =
           findMinMaxValues(weights, closePoints, elementData, j, numElements);
 
       if (info.maxIdx1 != -1 && info.maxIdx2 != -1 && info.minIdx1 != -1 &&
           info.minIdx2 != -1) {
-        weightsCopy[info.minIdx1] = NumericType(0);
-        weightsCopy[info.minIdx2] = NumericType(0);
-        weightsCopy[info.maxIdx1] = NumericType(0);
-        weightsCopy[info.maxIdx2] = NumericType(0);
+        weightsCopy[info.minIdx1] = 0.0;
+        weightsCopy[info.minIdx2] = 0.0;
+        weightsCopy[info.maxIdx1] = 0.0;
+        weightsCopy[info.maxIdx2] = 0.0;
       }
     } else if (discard2 && numClosePoints > 2) {
       const auto info =
           findMinMaxValues(weights, closePoints, elementData, j, numElements);
 
       if (info.minIdx1 != -1 && info.maxIdx1 != -1) {
-        weightsCopy[info.minIdx1] = NumericType(0);
-        weightsCopy[info.maxIdx1] = NumericType(0);
+        weightsCopy[info.minIdx1] = 0.0;
+        weightsCopy[info.maxIdx1] = 0.0;
       }
     }
+
+    return weightsCopy;
   }
 
   template <class AT, class BT, std::size_t D>

include/viennaps/psVersion.hpp

@@ -3,10 +3,10 @@
 namespace viennaps {
 
 // Version information generated by CMake
-inline constexpr const char *version = "4.1.2";
+inline constexpr const char *version = "4.1.3";
 inline constexpr int versionMajor = 4;
 inline constexpr int versionMinor = 1;
-inline constexpr int versionPatch = 2;
+inline constexpr int versionPatch = 3;
 
 // Utility functions for version comparison
 inline constexpr uint32_t versionAsInteger() {

pyproject.toml

@@ -3,7 +3,7 @@ requires = ["scikit-build-core", "pybind11"]
 build-backend = "scikit_build_core.build"
 
 [project]
-version = "4.1.2"
+version = "4.1.3"
 name = "ViennaPS"
 readme = "README.md"
 license = { file = "LICENSE" }

python/viennaps/__init__.pyi

@@ -235,7 +235,7 @@ def setDimension(d: int):
     """
 
 PROXY_DIM: int = 2
-__version__: str = "4.1.2"
+__version__: str = "4.1.3"
 ptxPath: str = ""
-version: str = "4.1.2"
+version: str = "4.1.3"
 _C = _core

python/viennaps/_core/__init__.pyi

@@ -9,8 +9,8 @@ import typing
 import viennals._core
 from viennaps import d2
 import viennaps.d2
-from viennaps import d3
 import viennaps.d3
+from viennaps import d3
 from . import constants
 from . import gpu
 from . import util
@@ -1162,5 +1162,5 @@ def gpuAvailable() -> bool:
 
 def setNumThreads(arg0: typing.SupportsInt) -> None: ...
 
-__version__: str = "4.1.2"
-version: str = "4.1.2"
+__version__: str = "4.1.3"
+version: str = "4.1.3"

tests/removeStrayPoints/CMakeLists.txt

@@ -1,7 +1,10 @@
-# project(removeStrayPoints LANGUAGES CXX)
+project(removeStrayPoints LANGUAGES CXX)
 
-# add_executable(${PROJECT_NAME} "removeStrayPoints.cpp")
-# target_link_libraries(${PROJECT_NAME} PRIVATE ViennaPS)
+# TODO: This test gets stuck on MSVC for some reason...
+if(NOT MSVC)
+  add_executable(${PROJECT_NAME} "removeStrayPoints.cpp")
+  target_link_libraries(${PROJECT_NAME} PRIVATE ViennaPS)
 
-# add_dependencies(ViennaPS_Tests ${PROJECT_NAME})
-# add_test(NAME ${PROJECT_NAME} COMMAND $<TARGET_FILE:${PROJECT_NAME}>)
+  add_dependencies(ViennaPS_Tests ${PROJECT_NAME})
+  add_test(NAME ${PROJECT_NAME} COMMAND $<TARGET_FILE:${PROJECT_NAME}>)
+endif()