Allow velocity calculations between steps in high-order temporal schemes

PythonAPI.md

@@ -18,6 +18,11 @@ They must be imported from either `viennals.d2` or `viennals.d3`.
 - `CalculateNormalVectors`
 - `CalculateVisibilities`
 - `Check`
+- `CompareChamfer`
+- `CompareCriticalDimensions`
+- `CompareNarrowBand`
+- `CompareSparseField`
+- `CompareVolume`
 - `PointCloud`
 - `ConvexHull`
 - `DetectFeatures`
@@ -49,8 +54,6 @@ They must be imported from either `viennals.d2` or `viennals.d3`.
 
 ### **2D-Only Functions**
 - `CompareArea`
-- `CompareNarrowBand`
-- `CompareSparseField`
 
 ---
 

examples/AirGapDeposition/AirGapDeposition.cpp

@@ -86,7 +86,7 @@ double runSimulation(AdvectKernelType &kernel,
 int main() {
 
   constexpr int D = 2;
-  omp_set_num_threads(8);
+  omp_set_num_threads(16);
 
   NumericType extent = 30;
   NumericType gridDelta = 0.5;

include/viennals/lsAdvect.hpp

@@ -2,6 +2,7 @@
 
 #include <lsPreCompileMacros.hpp>
 
+#include <functional>
 #include <limits>
 #include <vector>
 
@@ -78,6 +79,8 @@ template <class T, int D> class Advect {
   unsigned adaptiveTimeStepSubdivisions = 20;
   static constexpr double wrappingLayerEpsilon = 1e-4;
   SmartPointer<Domain<T, D>> originalLevelSet = nullptr;
+  std::function<bool(SmartPointer<Domain<T, D>>)> velocityUpdateCallback =
+      nullptr;
 
   // this vector will hold the maximum time step for each point and the
   // corresponding velocity
@@ -951,6 +954,13 @@ template <class T, int D> class Advect {
   /// be translated to the advected LS. Defaults to true.
   void setUpdatePointData(bool update) { updatePointData = update; }
 
+  /// Set a callback function that is called after the level set has been
+  /// updated during intermediate time integration steps (e.g. RK2, RK3).
+  void setVelocityUpdateCallback(
+      std::function<bool(SmartPointer<Domain<T, D>>)> callback) {
+    velocityUpdateCallback = callback;
+  }
+
   // Prepare the levelset for advection, based on the provided spatial
   // discretization scheme.
   void prepareLS() {

include/viennals/lsAdvectIntegrationSchemes.hpp

@@ -74,6 +74,9 @@ template <class T, int D> struct AdvectTimeIntegration {
     // Stage 1: u^(1) = u^n + dt * L(u^n)
     kernel.updateLevelSet(dt);
 
+    if (kernel.velocityUpdateCallback)
+      kernel.velocityUpdateCallback(kernel.levelSets.back());
+
     // Stage 2: u^(n+1) = 1/2 u^n + 1/2 (u^(1) + dt * L(u^(1)))
     // Current level set is u^(1). Compute L(u^(1)).
     kernel.computeRates(dt);
@@ -82,6 +85,11 @@ template <class T, int D> struct AdvectTimeIntegration {
     // Combine: u^(n+1) = 0.5 * u^n + 0.5 * u*
     kernel.combineLevelSets(0.5, 0.5);
 
+    // If we are in single step mode, the level set will be rebuilt immediately
+    // after this, invalidating the velocity field. Thus, we skip the update.
+    if (kernel.velocityUpdateCallback && !kernel.performOnlySingleStep)
+      kernel.velocityUpdateCallback(kernel.levelSets.back());
+
     // Finalize
     kernel.rebuildLS();
 
@@ -109,18 +117,29 @@ template <class T, int D> struct AdvectTimeIntegration {
     // Stage 1: u^(1) = u^n + dt * L(u^n)
     kernel.updateLevelSet(dt);
 
+    if (kernel.velocityUpdateCallback)
+      kernel.velocityUpdateCallback(kernel.levelSets.back());
+
     // Stage 2: u^(2) = 3/4 u^n + 1/4 (u^(1) + dt * L(u^(1)))
     kernel.computeRates(dt);
     kernel.updateLevelSet(dt);
     // Combine to get u^(2) = 0.75 * u^n + 0.25 * u*.
     kernel.combineLevelSets(0.75, 0.25);
 
+    if (kernel.velocityUpdateCallback)
+      kernel.velocityUpdateCallback(kernel.levelSets.back());
+
     // Stage 3: u^(n+1) = 1/3 u^n + 2/3 (u^(2) + dt * L(u^(2)))
     kernel.computeRates(dt);
     kernel.updateLevelSet(dt);
     // Combine to get u^(n+1) = 1/3 * u^n + 2/3 * u**.
     kernel.combineLevelSets(1.0 / 3.0, 2.0 / 3.0);
 
+    // If we are in single step mode, the level set will be rebuilt immediately
+    // after this, invalidating the velocity field. Thus, we skip the update.
+    if (kernel.velocityUpdateCallback && !kernel.performOnlySingleStep)
+      kernel.velocityUpdateCallback(kernel.levelSets.back());
+
     // Finalize: Re-segment and renormalize the final result.
     kernel.rebuildLS();
 

include/viennals/lsCalculateVisibilities.hpp

@@ -26,109 +26,120 @@ template <class T, int D> class CalculateVisibilities {
     auto &domain = levelSet->getDomain();
     auto &grid = levelSet->getGrid();
 
-    // *** Determine extents of domain ***
-    Vec3D<T> minDefinedPoint;
-    Vec3D<T> maxDefinedPoint;
-    // Initialize with extreme values
+    // Invert the vector
+    auto dir = Normalize(Inv(direction));
+    std::vector<T> visibilities(domain.getNumberOfPoints(), static_cast<T>(-1));
+
+    // Check if the direction vector has a component in the simulation
+    // dimensions
+    T dirNorm = 0;
     for (int i = 0; i < D; ++i) {
-      minDefinedPoint[i] = std::numeric_limits<T>::max();
-      maxDefinedPoint[i] = std::numeric_limits<T>::lowest();
+      dirNorm += dir[i] * dir[i];
     }
-    // Iterate through all defined points in the domain
-    for (viennahrle::SparseIterator<hrleDomainType> it(domain);
-         !it.isFinished(); it.next()) {
-      if (!it.isDefined())
-        continue; // Skip undefined points
-
-      // Get the coordinate of the current point
-      auto point = it.getStartIndices();
+
+    if (dirNorm < epsilon) {
+      std::fill(visibilities.begin(), visibilities.end(), static_cast<T>(1.0));
+    } else {
+      // *** Determine extents of domain ***
+      Vec3D<T> minDefinedPoint;
+      Vec3D<T> maxDefinedPoint;
+      // Initialize with extreme values
       for (int i = 0; i < D; ++i) {
-        // Compare to update min and max defined points
-        T coord = point[i]; // * grid.getGridDelta();
-        minDefinedPoint[i] = std::min(minDefinedPoint[i], coord);
-        maxDefinedPoint[i] = std::max(maxDefinedPoint[i], coord);
+        minDefinedPoint[i] = std::numeric_limits<T>::max();
+        maxDefinedPoint[i] = std::numeric_limits<T>::lowest();
       }
-    }
-    //****************************
+      // Iterate through all defined points in the domain
+      for (viennahrle::SparseIterator<hrleDomainType> it(domain);
+           !it.isFinished(); it.next()) {
+        if (!it.isDefined())
+          continue; // Skip undefined points
 
-    // Invert the vector
-    auto dir = Normalize(Inv(direction));
-    std::vector<T> visibilities(domain.getNumberOfPoints(), static_cast<T>(-1));
+        // Get the coordinate of the current point
+        auto point = it.getStartIndices();
+        for (int i = 0; i < D; ++i) {
+          // Compare to update min and max defined points
+          T coord = point[i]; // * grid.getGridDelta();
+          minDefinedPoint[i] = std::min(minDefinedPoint[i], coord);
+          maxDefinedPoint[i] = std::max(maxDefinedPoint[i], coord);
+        }
+      }
+      //****************************
 
 #pragma omp parallel num_threads(levelSet->getNumberOfSegments())
-    {
-      int p = 0;
+      {
+        int p = 0;
 #ifdef _OPENMP
-      p = omp_get_thread_num();
+        p = omp_get_thread_num();
 #endif
 
-      const viennahrle::Index<D> startVector =
-          (p == 0) ? grid.getMinGridPoint() : domain.getSegmentation()[p - 1];
+        const viennahrle::Index<D> startVector =
+            (p == 0) ? grid.getMinGridPoint() : domain.getSegmentation()[p - 1];
 
-      const viennahrle::Index<D> endVector =
-          (p != static_cast<int>(domain.getNumberOfSegments() - 1))
-              ? domain.getSegmentation()[p]
-              : grid.incrementIndices(grid.getMaxGridPoint());
+        const viennahrle::Index<D> endVector =
+            (p != static_cast<int>(domain.getNumberOfSegments() - 1))
+                ? domain.getSegmentation()[p]
+                : grid.incrementIndices(grid.getMaxGridPoint());
 
-      for (viennahrle::SparseIterator<hrleDomainType> it(domain, startVector);
-           it.getStartIndices() < endVector; ++it) {
+        for (viennahrle::SparseIterator<hrleDomainType> it(domain, startVector);
+             it.getStartIndices() < endVector; ++it) {
 
-        if (!it.isDefined())
-          continue;
+          if (!it.isDefined())
+            continue;
 
-        // Starting position of the point
-        Vec3D<T> currentPos;
-        for (int i = 0; i < D; ++i) {
-          currentPos[i] = it.getStartIndices(i);
-        }
-
-        // Start tracing the ray
-        T minLevelSetValue = it.getValue(); // Starting level set value
-        Vec3D<T> rayPos = currentPos;
-
-        // Step once to skip immediate neighbor
-        for (int i = 0; i < D; ++i)
-          rayPos[i] += dir[i];
+          // Starting position of the point
+          Vec3D<T> currentPos;
+          for (int i = 0; i < D; ++i) {
+            currentPos[i] = it.getStartIndices(i);
+          }
 
-        bool visibility = true;
+          // Start tracing the ray
+          T minLevelSetValue = it.getValue(); // Starting level set value
+          Vec3D<T> rayPos = currentPos;
 
-        while (true) {
-          // Update the ray position
+          // Step once to skip immediate neighbor
           for (int i = 0; i < D; ++i)
             rayPos[i] += dir[i];
 
-          // Determine the nearest grid cell (round to nearest index)
-          viennahrle::Index<D> nearestCell;
-          for (int i = 0; i < D; ++i)
-            nearestCell[i] = static_cast<viennahrle::IndexType>(rayPos[i]);
-
-          // Check if the nearest cell is within bounds
-          bool outOfBounds = false;
-          for (int i = 0; i < D; ++i) {
-            if (nearestCell[i] < minDefinedPoint[i] ||
-                nearestCell[i] > maxDefinedPoint[i]) {
-              outOfBounds = true;
-              break;
+          bool visibility = true;
+
+          while (true) {
+            // Update the ray position
+            for (int i = 0; i < D; ++i)
+              rayPos[i] += dir[i];
+
+            // Determine the nearest grid cell (round to nearest index)
+            viennahrle::Index<D> nearestCell;
+            for (int i = 0; i < D; ++i)
+              nearestCell[i] = static_cast<viennahrle::IndexType>(rayPos[i]);
+
+            // Check if the nearest cell is within bounds
+            bool outOfBounds = false;
+            for (int i = 0; i < D; ++i) {
+              if (nearestCell[i] < minDefinedPoint[i] ||
+                  nearestCell[i] > maxDefinedPoint[i]) {
+                outOfBounds = true;
+                break;
+              }
             }
-          }
 
-          if (outOfBounds)
-            break; // Ray is outside the grid
+            if (outOfBounds)
+              break; // Ray is outside the grid
 
-          // Access the level set value at the nearest cell
-          T value =
-              viennahrle::SparseIterator<hrleDomainType>(domain, nearestCell)
-                  .getValue();
+            // Access the level set value at the nearest cell
+            T value =
+                viennahrle::SparseIterator<hrleDomainType>(domain, nearestCell)
+                    .getValue();
 
-          // Update the minimum value encountered
-          if (value < minLevelSetValue) {
-            visibility = false;
-            break;
+            // Update the minimum value encountered
+            if (value < minLevelSetValue) {
+              visibility = false;
+              break;
+            }
           }
-        }
 
-        // Update visibility for this point
-        visibilities[it.getPointId()] = visibility ? 1.0 : 0.0;
+          // Update visibility for this point
+          visibilities[it.getPointId()] = visibility ? 1.0 : 0.0;
+        }
       }
     }