diff --git a/examples/Epitaxy/Epitaxy.cpp b/examples/Epitaxy/Epitaxy.cpp index d744c25c..9cab0203 100644 --- a/examples/Epitaxy/Epitaxy.cpp +++ b/examples/Epitaxy/Epitaxy.cpp @@ -47,6 +47,7 @@ void writeSurface(SmartPointer> 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] = { diff --git a/examples/VolumeToLevelSets/VolumeToLevelSets.cpp b/examples/VolumeToLevelSets/VolumeToLevelSets.cpp index fbabab4c..68ef3236 100644 --- a/examples/VolumeToLevelSets/VolumeToLevelSets.cpp +++ b/examples/VolumeToLevelSets/VolumeToLevelSets.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include #include @@ -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>::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 materials = {0, 0, 1, 1, 1}; + mesh->cellData.insertNextScalarData(materials, "Material"); + + ls::VTKWriter(mesh, ls::FileFormatEnum::VTU, fileName).apply(); } auto mesh = ls::SmartPointer>::New(); @@ -32,7 +47,10 @@ int main(int argc, char *argv[]) { std::vector 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(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>::New( bounds, boundaryCons, gridDelta); @@ -61,6 +79,15 @@ int main(int argc, char *argv[]) { ls::ToSurfaceMesh(levelSets[i], mesh).apply(); ls::VTKWriter(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>::New(); + ls::ToMesh(levelSets[i], meshVol).apply(); + ls::VTKWriter(meshVol, "LSvolume-" + std::to_string(i) + ".vtu") + .apply(); } return 0; diff --git a/include/viennals/lsCompareSparseField.hpp b/include/viennals/lsCompareSparseField.hpp index e4b70e9c..0b012476 100644 --- a/include/viennals/lsCompareSparseField.hpp +++ b/include/viennals/lsCompareSparseField.hpp @@ -133,20 +133,12 @@ template class CompareSparseField { } public: - CompareSparseField() { - assert( - D == 2 && - "CompareSparseField is currently only implemented for 2D level sets."); - } + CompareSparseField() {} CompareSparseField(SmartPointer> passedLevelSetExpanded, SmartPointer> passedLevelSetIterated) : levelSetExpanded(passedLevelSetExpanded), - levelSetIterated(passedLevelSetIterated) { - assert( - D == 2 && - "CompareSparseField is currently only implemented for 2D level sets."); - } + levelSetIterated(passedLevelSetIterated) {} void setLevelSetExpanded(SmartPointer> passedLevelSet) { levelSetExpanded = passedLevelSet; @@ -278,7 +270,8 @@ template 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)) { diff --git a/include/viennals/lsConvexHull.hpp b/include/viennals/lsConvexHull.hpp index 89450c60..fe583f12 100644 --- a/include/viennals/lsConvexHull.hpp +++ b/include/viennals/lsConvexHull.hpp @@ -42,6 +42,26 @@ template 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 previousCandidates; + for (unsigned i = 0; i < points.size(); ++i) { if (i == currentEdge[0] || i == nextIndex) continue; @@ -67,6 +87,8 @@ template 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 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; diff --git a/tests/CompareSparseField/CompareSparseField.cpp b/tests/CompareSparseField/CompareSparseField.cpp index 370d6c4b..a43ddb43 100644 --- a/tests/CompareSparseField/CompareSparseField.cpp +++ b/tests/CompareSparseField/CompareSparseField.cpp @@ -22,16 +22,17 @@ namespace ls = viennals; -int main() { - constexpr int D = 2; - - omp_set_num_threads(4); - +template void runTest() { double extent = 15; double gridDelta = 0.5; - double bounds[2 * D] = {-extent, extent, -extent, extent}; - ls::Domain::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::BoundaryType boundaryCons[D]; for (unsigned i = 0; i < D; ++i) boundaryCons[i] = ls::Domain::BoundaryType::REFLECTIVE_BOUNDARY; @@ -39,7 +40,9 @@ int main() { auto sphere1 = ls::SmartPointer>::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( @@ -50,7 +53,12 @@ int main() { auto sphere2 = ls::SmartPointer>::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(sphere2, 1).apply(); + std::string dimString = std::to_string(D) + "D"; + // Export both spheres as VTK files for visualization { auto mesh = ls::SmartPointer>::New(); ls::ToMesh(sphere1, mesh).apply(); - ls::VTKWriter(mesh, "sphere1_expanded.vtp").apply(); + ls::VTKWriter(mesh, "sphere1_expanded_" + dimString + ".vtp") + .apply(); auto meshSurface = ls::SmartPointer>::New(); ls::ToSurfaceMesh(sphere1, meshSurface).apply(); - ls::VTKWriter(meshSurface, "sphere1_surface.vtp").apply(); + ls::VTKWriter(meshSurface, "sphere1_surface_" + dimString + ".vtp") + .apply(); } { auto mesh = ls::SmartPointer>::New(); ls::ToMesh(sphere2, mesh).apply(); - ls::VTKWriter(mesh, "sphere2_sparse_iterated.vtp").apply(); + ls::VTKWriter(mesh, "sphere2_sparse_iterated_" + dimString + ".vtp") + .apply(); auto meshSurface = ls::SmartPointer>::New(); ls::ToSurfaceMesh(sphere2, meshSurface).apply(); - ls::VTKWriter(meshSurface, "sphere2_surface.vtp").apply(); + ls::VTKWriter(meshSurface, "sphere2_surface_" + dimString + ".vtp") + .apply(); } // Compare using sparse field comparison @@ -102,11 +116,12 @@ int main() { auto meshWithPointData = ls::SmartPointer>::New(); // Mesh with point data ls::ToMesh(sphere2, meshWithPointData).apply(); - ls::VTKWriter(meshWithPointData, "sphere2_LS_with_point_data.vtp") + ls::VTKWriter(meshWithPointData, + "sphere2_LS_with_point_data_" + dimString + ".vtp") .apply(); // Save mesh to file - ls::VTKWriter(mesh, "sparsefield.vtp").apply(); + ls::VTKWriter(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(mesh, "sparsefield_restricted.vtp").apply(); + ls::VTKWriter(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; }