Fix failing examples

examples/Epitaxy/Epitaxy.cpp

@@ -47,6 +47,7 @@ void writeSurface(SmartPointer<Domain<T, D>> domain,
 
 int main(int argc, char *argv[]) {
 
+  omp_set_num_threads(8);
   // Create hole geometry
   double bounds[2 * D] = {-1.0, 1.0, -1.0, 1.0, -1.0, 1.0};
   BoundaryConditionEnum boundaryConditions[2 * D] = {

examples/VolumeToLevelSets/VolumeToLevelSets.cpp

@@ -2,6 +2,7 @@
 
 #include <lsDomain.hpp>
 #include <lsFromVolumeMesh.hpp>
+#include <lsToMesh.hpp>
 #include <lsToSurfaceMesh.hpp>
 #include <lsVTKReader.hpp>
 #include <lsVTKWriter.hpp>
@@ -12,13 +13,27 @@ int main(int argc, char *argv[]) {
 
   using NumericType = double;
   constexpr int D = 3;
-  double gridDelta = 0.00023;
+  double gridDelta = 0.1;
 
   std::string fileName;
   if (argc > 1) {
     fileName = std::string(argv[1]);
   } else {
     fileName = "volumeInitial.vtu";
+    // Generate a simple volume mesh if no file is provided
+    auto mesh = ls::SmartPointer<ls::Mesh<NumericType>>::New();
+    // Create a simple cube [-1, 1]^3 consisting of 5 tetrahedra
+    mesh->nodes = {{-1., -1., -1.}, {1., -1., -1.}, {1., 1., -1.},
+                   {-1., 1., -1.},  {-1., -1., 1.}, {1., -1., 1.},
+                   {1., 1., 1.},    {-1., 1., 1.}};
+    // 5-tetra decomposition of a cube
+    mesh->tetras = {
+        {0, 1, 3, 4}, {1, 2, 3, 6}, {1, 4, 5, 6}, {3, 4, 6, 7}, {1, 3, 4, 6}};
+    // Assign materials
+    std::vector<double> materials = {0, 0, 1, 1, 1};
+    mesh->cellData.insertNextScalarData(materials, "Material");
+
+    ls::VTKWriter<NumericType>(mesh, ls::FileFormatEnum::VTU, fileName).apply();
   }
 
   auto mesh = ls::SmartPointer<ls::Mesh<NumericType>>::New();
@@ -32,7 +47,10 @@ int main(int argc, char *argv[]) {
     std::vector<int> translator = {3, 2, 4, 7, 7, 6, 5, 7, 1, 0};
     if (materialData != nullptr) {
       for (auto &cell : *materialData) {
-        cell = translator[std::round(cell)];
+        int id = std::round(cell);
+        if (id >= 0 && static_cast<std::size_t>(id) < translator.size()) {
+          cell = translator[id];
+        }
       }
     }
   }
@@ -41,12 +59,12 @@ int main(int argc, char *argv[]) {
 
   ls::VTKWriter(mesh, ls::FileFormatEnum::VTU, "ReadVolumeMesh.vtu").apply();
 
-  double bounds[2 * D] = {-6, 6, 1e-10, 0.078, -0.034, 0.034};
+  double bounds[2 * D] = {-2, 2, -2, 2, -2, 2};
   ls::BoundaryConditionEnum boundaryCons[D];
   for (unsigned i = 0; i < D; ++i) {
     boundaryCons[i] = ls::BoundaryConditionEnum::REFLECTIVE_BOUNDARY;
   }
-  boundaryCons[0] = ls::BoundaryConditionEnum::INFINITE_BOUNDARY;
+  // boundaryCons[0] = ls::BoundaryConditionEnum::INFINITE_BOUNDARY;
 
   auto domain = ls::SmartPointer<ls::Domain<NumericType, D>>::New(
       bounds, boundaryCons, gridDelta);
@@ -61,6 +79,15 @@ int main(int argc, char *argv[]) {
     ls::ToSurfaceMesh<double, D>(levelSets[i], mesh).apply();
     ls::VTKWriter<double>(mesh, "LSsurface-" + std::to_string(i) + ".vtp")
         .apply();
+
+    std::cout << "---------------------------------" << std::endl;
+    std::cout << "LevelSet " << i << std::endl;
+    levelSets[i]->print();
+
+    auto meshVol = ls::SmartPointer<ls::Mesh<>>::New();
+    ls::ToMesh<double, D>(levelSets[i], meshVol).apply();
+    ls::VTKWriter<double>(meshVol, "LSvolume-" + std::to_string(i) + ".vtu")
+        .apply();
   }
 
   return 0;

include/viennals/lsCompareSparseField.hpp

@@ -133,20 +133,12 @@ template <class T, int D = 2> class CompareSparseField {
   }
 
 public:
-  CompareSparseField() {
-    assert(
-        D == 2 &&
-        "CompareSparseField is currently only implemented for 2D level sets.");
-  }
+  CompareSparseField() {}
 
   CompareSparseField(SmartPointer<Domain<T, D>> passedLevelSetExpanded,
                      SmartPointer<Domain<T, D>> passedLevelSetIterated)
       : levelSetExpanded(passedLevelSetExpanded),
-        levelSetIterated(passedLevelSetIterated) {
-    assert(
-        D == 2 &&
-        "CompareSparseField is currently only implemented for 2D level sets.");
-  }
+        levelSetIterated(passedLevelSetIterated) {}
 
   void setLevelSetExpanded(SmartPointer<Domain<T, D>> passedLevelSet) {
     levelSetExpanded = passedLevelSet;
@@ -278,7 +270,8 @@ template <class T, int D = 2> class CompareSparseField {
       // Calculate coordinates
       T xCoord = indices[0] * gridDelta;
       T yCoord = indices[1] * gridDelta;
-      T zCoord = 0.0; // Always use 0 for z-coordinate in 2D
+      T zCoord = (D == 3) ? indices[2] * gridDelta
+                          : 0.0; // Always use 0 for z-coordinate in 2D
 
       // Skip if outside the specified x-range
       if (useXRange && (xCoord < xRangeMin || xCoord > xRangeMax)) {

include/viennals/lsConvexHull.hpp

@@ -42,6 +42,26 @@ template <class T, int D> class ConvexHull {
       ++nextIndex;
     }
 
+    // Robust initialization: find a point that is not collinear
+    if constexpr (D == 3) {
+      unsigned searchIndex = nextIndex;
+      for (; searchIndex < points.size(); ++searchIndex) {
+        if (searchIndex == currentEdge[0] || searchIndex == currentEdge[1])
+          continue;
+
+        auto v1 = points[currentEdge[1]] - points[currentEdge[0]];
+        auto v2 = points[searchIndex] - points[currentEdge[0]];
+        auto cross = CrossProduct(v1, v2);
+        if (DotProduct(cross, cross) > 1e-10) {
+          nextIndex = searchIndex;
+          break;
+        }
+      }
+    }
+
+    // Keep track of candidates we've already selected to prevent cycles
+    std::vector<unsigned> previousCandidates;
+
     for (unsigned i = 0; i < points.size(); ++i) {
       if (i == currentEdge[0] || i == nextIndex)
         continue;
@@ -67,6 +87,8 @@ template <class T, int D> class ConvexHull {
       }
 
       auto product = DotProduct(distance, normal);
+      bool update = false;
+
       // if dot product is very small, point is very close to plane
       // we need to check if we already have the correct point
       // or if it is the next correct one
@@ -88,12 +110,37 @@ template <class T, int D> class ConvexHull {
         edges[1][1] = i;
 
         if (!wasEdgeVisited(edges[0]) && !wasEdgeVisited(edges[1])) {
-          nextIndex = i;
+          // Tie-breaker: prefer closer points to avoid long chords
+          auto distI = DotProduct(distance, distance);
+          auto distNextVec = points[nextIndex] - points[currentEdge[0]];
+          auto distNext = DotProduct(distNextVec, distNextVec);
+          if (distI < distNext) {
+            update = true;
+          }
         }
       }
       // check if point is to the right of current element
       else if (product > 0) {
+        update = true;
+      }
+
+      if (update) {
+        // Check for cycles
+        bool cycleDetected = false;
+        for (unsigned prev : previousCandidates) {
+          if (prev == i) {
+            cycleDetected = true;
+            break;
+          }
+        }
+
+        if (cycleDetected) {
+          continue;
+        }
+
+        previousCandidates.push_back(nextIndex);
         nextIndex = i;
+        i = -1;
       }
     }
     return nextIndex;

tests/CompareSparseField/CompareSparseField.cpp

@@ -22,24 +22,27 @@
 
 namespace ls = viennals;
 
-int main() {
-  constexpr int D = 2;
-
-  omp_set_num_threads(4);
-
+template <int D> void runTest() {
   double extent = 15;
   double gridDelta = 0.5;
 
-  double bounds[2 * D] = {-extent, extent, -extent, extent};
-  ls::Domain<double, D>::BoundaryType boundaryCons[D];
+  double bounds[2 * D];
+  for (unsigned i = 0; i < D; ++i) {
+    bounds[2 * i] = -extent;
+    bounds[2 * i + 1] = extent;
+  }
+
+  typename ls::Domain<double, D>::BoundaryType boundaryCons[D];
   for (unsigned i = 0; i < D; ++i)
     boundaryCons[i] = ls::Domain<double, D>::BoundaryType::REFLECTIVE_BOUNDARY;
 
   // Create first sphere (target)
   auto sphere1 = ls::SmartPointer<ls::Domain<double, D>>::New(
       bounds, boundaryCons, gridDelta);
 
-  double origin1[D] = {0., 0.};
+  double origin1[D];
+  for (int i = 0; i < D; ++i)
+    origin1[i] = 0.;
   double radius1 = 5.0;
 
   ls::MakeGeometry<double, D>(
@@ -50,7 +53,12 @@ int main() {
   auto sphere2 = ls::SmartPointer<ls::Domain<double, D>>::New(
       bounds, boundaryCons, gridDelta);
 
-  double origin2[D] = {2., 1.}; // Shifted center
+  double origin2[D];
+  for (int i = 0; i < D; ++i)
+    origin2[i] = 0.;
+  origin2[0] = 2.;
+  if (D > 1)
+    origin2[1] = 1.;
   // double origin2[D] = {0., 0.}; // Same center
   double radius2 = 5.0; // Same radius
 
@@ -65,23 +73,29 @@ int main() {
   // Reduce the sample level set to a sparse field
   ls::Reduce<double, D>(sphere2, 1).apply();
 
+  std::string dimString = std::to_string(D) + "D";
+
   // Export both spheres as VTK files for visualization
   {
     auto mesh = ls::SmartPointer<ls::Mesh<>>::New();
     ls::ToMesh<double, D>(sphere1, mesh).apply();
-    ls::VTKWriter<double>(mesh, "sphere1_expanded.vtp").apply();
+    ls::VTKWriter<double>(mesh, "sphere1_expanded_" + dimString + ".vtp")
+        .apply();
     auto meshSurface = ls::SmartPointer<ls::Mesh<>>::New();
     ls::ToSurfaceMesh<double, D>(sphere1, meshSurface).apply();
-    ls::VTKWriter<double>(meshSurface, "sphere1_surface.vtp").apply();
+    ls::VTKWriter<double>(meshSurface, "sphere1_surface_" + dimString + ".vtp")
+        .apply();
   }
 
   {
     auto mesh = ls::SmartPointer<ls::Mesh<>>::New();
     ls::ToMesh<double, D>(sphere2, mesh).apply();
-    ls::VTKWriter<double>(mesh, "sphere2_sparse_iterated.vtp").apply();
+    ls::VTKWriter<double>(mesh, "sphere2_sparse_iterated_" + dimString + ".vtp")
+        .apply();
     auto meshSurface = ls::SmartPointer<ls::Mesh<>>::New();
     ls::ToSurfaceMesh<double, D>(sphere2, meshSurface).apply();
-    ls::VTKWriter<double>(meshSurface, "sphere2_surface.vtp").apply();
+    ls::VTKWriter<double>(meshSurface, "sphere2_surface_" + dimString + ".vtp")
+        .apply();
   }
 
   // Compare using sparse field comparison
@@ -102,11 +116,12 @@ int main() {
   auto meshWithPointData =
       ls::SmartPointer<ls::Mesh<>>::New(); // Mesh with point data
   ls::ToMesh<double, D>(sphere2, meshWithPointData).apply();
-  ls::VTKWriter<double>(meshWithPointData, "sphere2_LS_with_point_data.vtp")
+  ls::VTKWriter<double>(meshWithPointData,
+                        "sphere2_LS_with_point_data_" + dimString + ".vtp")
       .apply();
 
   // Save mesh to file
-  ls::VTKWriter<double>(mesh, "sparsefield.vtp").apply();
+  ls::VTKWriter<double>(mesh, "sparsefield_" + dimString + ".vtp").apply();
 
   // Get the calculated difference metrics
   double sumSquaredDifferences = compareSparseField.getSumSquaredDifferences();
@@ -117,11 +132,15 @@ int main() {
   unsigned numSkippedPoints =
       compareSparseField.getNumSkippedPoints(); // Number of skipped points
 
-  std::cout << "\nComparison Results:" << std::endl;
-  std::cout << "Sphere 1 center: (" << origin1[0] << ", " << origin1[1] << ")"
-            << std::endl;
-  std::cout << "Sphere 2 center: (" << origin2[0] << ", " << origin2[1] << ")"
-            << std::endl;
+  std::cout << "\nComparison Results (" << dimString << "):" << std::endl;
+  std::cout << "Sphere 1 center: (";
+  for (int i = 0; i < D; ++i)
+    std::cout << origin1[i] << ((i == D - 1) ? "" : ", ");
+  std::cout << ")" << std::endl;
+  std::cout << "Sphere 2 center: (";
+  for (int i = 0; i < D; ++i)
+    std::cout << origin2[i] << ((i == D - 1) ? "" : ", ");
+  std::cout << ")" << std::endl;
   std::cout << "Sphere 1 level set width after expansion: "
             << sphere1->getLevelSetWidth() << std::endl;
   std::cout << "Sum of squared differences: " << sumSquaredDifferences
@@ -177,7 +196,8 @@ int main() {
   // Create a mesh output with squared differences
   compareSparseField.setOutputMesh(mesh);
   compareSparseField.apply();
-  ls::VTKWriter<double>(mesh, "sparsefield_restricted.vtp").apply();
+  ls::VTKWriter<double>(mesh, "sparsefield_restricted_" + dimString + ".vtp")
+      .apply();
 
   // Test with different expansion widths
   std::cout << "\nTesting with different expansion widths:" << std::endl;
@@ -248,6 +268,11 @@ int main() {
   //           << std::endl;
   // std::cout << "Performance ratio: "
   //           << narrowband_ms.count() / sparse_ms.count() << "x" << std::endl;
+}
 
+int main() {
+  omp_set_num_threads(8);
+  runTest<2>();
+  runTest<3>();
   return 0;
 }