diff --git a/PythonAPI.md b/PythonAPI.md index 002965d9..2c335a26 100644 --- a/PythonAPI.md +++ b/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` --- diff --git a/examples/AirGapDeposition/AirGapDeposition.cpp b/examples/AirGapDeposition/AirGapDeposition.cpp index aa49d222..59448c29 100644 --- a/examples/AirGapDeposition/AirGapDeposition.cpp +++ b/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; diff --git a/include/viennals/lsAdvect.hpp b/include/viennals/lsAdvect.hpp index f6ecf3ad..33309638 100644 --- a/include/viennals/lsAdvect.hpp +++ b/include/viennals/lsAdvect.hpp @@ -2,6 +2,7 @@ #include +#include #include #include @@ -78,6 +79,8 @@ template class Advect { unsigned adaptiveTimeStepSubdivisions = 20; static constexpr double wrappingLayerEpsilon = 1e-4; SmartPointer> originalLevelSet = nullptr; + std::function>)> velocityUpdateCallback = + nullptr; // this vector will hold the maximum time step for each point and the // corresponding velocity @@ -951,6 +954,13 @@ template 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>)> callback) { + velocityUpdateCallback = callback; + } + // Prepare the levelset for advection, based on the provided spatial // discretization scheme. void prepareLS() { diff --git a/include/viennals/lsAdvectIntegrationSchemes.hpp b/include/viennals/lsAdvectIntegrationSchemes.hpp index beb80d41..691c3ba5 100644 --- a/include/viennals/lsAdvectIntegrationSchemes.hpp +++ b/include/viennals/lsAdvectIntegrationSchemes.hpp @@ -74,6 +74,9 @@ template 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 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 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(); diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp index d068cee0..33f10fb9 100644 --- a/include/viennals/lsCalculateVisibilities.hpp +++ b/include/viennals/lsCalculateVisibilities.hpp @@ -26,109 +26,120 @@ template class CalculateVisibilities { auto &domain = levelSet->getDomain(); auto &grid = levelSet->getGrid(); - // *** Determine extents of domain *** - Vec3D minDefinedPoint; - Vec3D maxDefinedPoint; - // Initialize with extreme values + // Invert the vector + auto dir = Normalize(Inv(direction)); + std::vector visibilities(domain.getNumberOfPoints(), static_cast(-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::max(); - maxDefinedPoint[i] = std::numeric_limits::lowest(); + dirNorm += dir[i] * dir[i]; } - // Iterate through all defined points in the domain - for (viennahrle::SparseIterator 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(1.0)); + } else { + // *** Determine extents of domain *** + Vec3D minDefinedPoint; + Vec3D 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::max(); + maxDefinedPoint[i] = std::numeric_limits::lowest(); } - } - //**************************** + // Iterate through all defined points in the domain + for (viennahrle::SparseIterator it(domain); + !it.isFinished(); it.next()) { + if (!it.isDefined()) + continue; // Skip undefined points - // Invert the vector - auto dir = Normalize(Inv(direction)); - std::vector visibilities(domain.getNumberOfPoints(), static_cast(-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 startVector = - (p == 0) ? grid.getMinGridPoint() : domain.getSegmentation()[p - 1]; + const viennahrle::Index startVector = + (p == 0) ? grid.getMinGridPoint() : domain.getSegmentation()[p - 1]; - const viennahrle::Index endVector = - (p != static_cast(domain.getNumberOfSegments() - 1)) - ? domain.getSegmentation()[p] - : grid.incrementIndices(grid.getMaxGridPoint()); + const viennahrle::Index endVector = + (p != static_cast(domain.getNumberOfSegments() - 1)) + ? domain.getSegmentation()[p] + : grid.incrementIndices(grid.getMaxGridPoint()); - for (viennahrle::SparseIterator it(domain, startVector); - it.getStartIndices() < endVector; ++it) { + for (viennahrle::SparseIterator it(domain, startVector); + it.getStartIndices() < endVector; ++it) { - if (!it.isDefined()) - continue; + if (!it.isDefined()) + continue; - // Starting position of the point - Vec3D 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 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 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 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 nearestCell; - for (int i = 0; i < D; ++i) - nearestCell[i] = static_cast(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 nearestCell; + for (int i = 0; i < D; ++i) + nearestCell[i] = static_cast(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(domain, nearestCell) - .getValue(); + // Access the level set value at the nearest cell + T value = + viennahrle::SparseIterator(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; + } } } diff --git a/include/viennals/lsCompareChamfer.hpp b/include/viennals/lsCompareChamfer.hpp index 4c1d4b89..dc34dbe6 100644 --- a/include/viennals/lsCompareChamfer.hpp +++ b/include/viennals/lsCompareChamfer.hpp @@ -27,8 +27,8 @@ using namespace viennacore; /// - RMS Chamfer distance: root mean square of nearest-neighbor distances /// - Maximum distance: maximum nearest-neighbor distance across both directions /// -/// The code is currently intended for 2D level sets only, where surfaces are -/// represented as line segments. +/// The code works for 2D and 3D level sets. Surfaces are represented as line +/// segments in 2D and triangles in 3D. /// /// Both level sets must have a width of at least 2 to extract surfaces. If /// not, they will be automatically expanded. @@ -73,20 +73,12 @@ template class CompareChamfer { } public: - CompareChamfer() { - static_assert( - D == 2, - "CompareChamfer is currently only implemented for 2D level sets."); - } + CompareChamfer() {} CompareChamfer(SmartPointer> passedLevelSetTarget, SmartPointer> passedLevelSetSample) : levelSetTarget(passedLevelSetTarget), - levelSetSample(passedLevelSetSample) { - static_assert( - D == 2, - "CompareChamfer is currently only implemented for 2D level sets."); - } + levelSetSample(passedLevelSetSample) {} /// Set the target level set void setLevelSetTarget(SmartPointer> passedLevelSet) { @@ -171,8 +163,7 @@ template class CompareChamfer { return; } - // Convert nodes to format suitable for KDTree (only using x, y - // coordinates) + // Convert nodes to format suitable for KDTree std::vector> targetPoints(numTargetPoints); std::vector> samplePoints(numSamplePoints); @@ -255,10 +246,11 @@ template class CompareChamfer { outputMeshTarget->clear(); outputMeshTarget->nodes = targetNodes; - // Create vertices for visualization - outputMeshTarget->vertices.reserve(numTargetPoints); - for (unsigned i = 0; i < numTargetPoints; ++i) { - outputMeshTarget->vertices.push_back({i}); + // Copy topology for visualization + if constexpr (D == 2) { + outputMeshTarget->lines = targetSurfaceMesh->lines; + } else { + outputMeshTarget->triangles = targetSurfaceMesh->triangles; } // Add distance data @@ -266,9 +258,13 @@ template class CompareChamfer { std::move(targetDistances), "DistanceToSample"); // Set mesh extent + for (unsigned d = 0; d < 3; ++d) { + outputMeshTarget->minimumExtent[d] = + (d < D) ? std::numeric_limits::max() : 0.0; + outputMeshTarget->maximumExtent[d] = + (d < D) ? std::numeric_limits::lowest() : 0.0; + } for (unsigned d = 0; d < D; ++d) { - outputMeshTarget->minimumExtent[d] = std::numeric_limits::max(); - outputMeshTarget->maximumExtent[d] = std::numeric_limits::lowest(); for (const auto &node : targetNodes) { outputMeshTarget->minimumExtent[d] = std::min(outputMeshTarget->minimumExtent[d], node[d]); @@ -282,10 +278,11 @@ template class CompareChamfer { outputMeshSample->clear(); outputMeshSample->nodes = sampleNodes; - // Create vertices for visualization - outputMeshSample->vertices.reserve(numSamplePoints); - for (unsigned i = 0; i < numSamplePoints; ++i) { - outputMeshSample->vertices.push_back({i}); + // Copy topology for visualization + if constexpr (D == 2) { + outputMeshSample->lines = sampleSurfaceMesh->lines; + } else { + outputMeshSample->triangles = sampleSurfaceMesh->triangles; } // Add distance data @@ -293,9 +290,13 @@ template class CompareChamfer { std::move(sampleDistances), "DistanceToTarget"); // Set mesh extent + for (unsigned d = 0; d < 3; ++d) { + outputMeshSample->minimumExtent[d] = + (d < D) ? std::numeric_limits::max() : 0.0; + outputMeshSample->maximumExtent[d] = + (d < D) ? std::numeric_limits::lowest() : 0.0; + } for (unsigned d = 0; d < D; ++d) { - outputMeshSample->minimumExtent[d] = std::numeric_limits::max(); - outputMeshSample->maximumExtent[d] = std::numeric_limits::lowest(); for (const auto &node : sampleNodes) { outputMeshSample->minimumExtent[d] = std::min(outputMeshSample->minimumExtent[d], node[d]); diff --git a/include/viennals/lsCompareCriticalDimensions.hpp b/include/viennals/lsCompareCriticalDimensions.hpp index fd851b7a..98e8270d 100644 --- a/include/viennals/lsCompareCriticalDimensions.hpp +++ b/include/viennals/lsCompareCriticalDimensions.hpp @@ -27,8 +27,6 @@ using namespace viennacore; /// zero. Multiple ranges can be specified to compare different critical /// dimensions. /// -/// The code is currently intended for 2D level sets only. -/// /// Note for the future: lsToDiskMesh could be used instead of lsToSurfaceMesh, /// which is probably more efficient but slightly less accurate. @@ -40,9 +38,9 @@ template class CompareCriticalDimensions { // Structure to hold a range specification struct RangeSpec { - bool isXRange; // true if X range, false if Y range - T rangeMin; - T rangeMax; + int measureDimension; + std::array minBounds; + std::array maxBounds; bool findMaximum; // true for maximum, false for minimum }; @@ -50,9 +48,9 @@ template class CompareCriticalDimensions { // Structure to hold critical dimension results struct CriticalDimensionResult { - bool isXRange; - T rangeMin; - T rangeMax; + int measureDimension; + std::array minBounds; + std::array maxBounds; bool findMaximum; T positionTarget; // Critical dimension position in target LS T positionSample; // Critical dimension position in sample LS @@ -94,32 +92,27 @@ template class CompareCriticalDimensions { // Extract surface positions from mesh nodes within the specified range // Returns the position coordinates (Y if isXRange=true, X if isXRange=false) std::vector findSurfaceCrossings(SmartPointer> surfaceMesh, - bool isXRange, T scanMin, T scanMax, - T perpMin, T perpMax) { + const RangeSpec &spec) { std::vector crossings; // Iterate through all surface mesh nodes for (const auto &node : surfaceMesh->nodes) { - T xCoord = node[0]; - T yCoord = node[1]; - // Check if point is in our scan range - bool inRange = false; - if (isXRange) { - // X range specified - check if X is in scan range - inRange = (xCoord >= scanMin && xCoord <= scanMax); - } else { - // Y range specified - check if Y is in scan range - inRange = (yCoord >= scanMin && yCoord <= scanMax); + bool inRange = true; + for (int i = 0; i < D; ++i) { + if (i == spec.measureDimension) + continue; + if (node[i] < spec.minBounds[i] || node[i] > spec.maxBounds[i]) { + inRange = false; + break; + } } if (inRange) { - // Extract perpendicular coordinate - T perpCoord = isXRange ? yCoord : xCoord; - - // Check if perpendicular coordinate is in range - if (perpCoord >= perpMin && perpCoord <= perpMax) { - crossings.push_back(perpCoord); + T val = node[spec.measureDimension]; + if (val >= spec.minBounds[spec.measureDimension] && + val <= spec.maxBounds[spec.measureDimension]) { + crossings.push_back(val); } } } @@ -142,18 +135,12 @@ template class CompareCriticalDimensions { } public: - CompareCriticalDimensions() { - static_assert(D == 2 && "CompareCriticalDimensions is currently only " - "implemented for 2D level sets."); - } + CompareCriticalDimensions() {} CompareCriticalDimensions(SmartPointer> passedLevelSetTarget, SmartPointer> passedLevelSetSample) : levelSetTarget(passedLevelSetTarget), - levelSetSample(passedLevelSetSample) { - static_assert(D == 2 && "CompareCriticalDimensions is currently only " - "implemented for 2D level sets."); - } + levelSetSample(passedLevelSetSample) {} void setLevelSetTarget(SmartPointer> passedLevelSet) { levelSetTarget = passedLevelSet; @@ -163,24 +150,39 @@ template class CompareCriticalDimensions { levelSetSample = passedLevelSet; } - /// Add an X range to find maximum or minimum Y position - void addXRange(T minX, T maxX, bool findMaximum = true) { + /// Add a generic range specification + void addRange(int measureDimension, const std::array &minBounds, + const std::array &maxBounds, bool findMaximum = true) { RangeSpec spec; - spec.isXRange = true; - spec.rangeMin = minX; - spec.rangeMax = maxX; + spec.measureDimension = measureDimension; + spec.minBounds = minBounds; + spec.maxBounds = maxBounds; spec.findMaximum = findMaximum; rangeSpecs.push_back(spec); } + /// Add an X range to find maximum or minimum Y position + void addXRange(T minX, T maxX, bool findMaximum = true) { + if constexpr (D == 2) { + std::array minBounds = {minX, std::numeric_limits::lowest()}; + std::array maxBounds = {maxX, std::numeric_limits::max()}; + addRange(1, minBounds, maxBounds, findMaximum); + } else { + VIENNACORE_LOG_WARNING( + "addXRange is only supported for 2D. Use addRange for 3D."); + } + } + /// Add a Y range to find maximum or minimum X position void addYRange(T minY, T maxY, bool findMaximum = true) { - RangeSpec spec; - spec.isXRange = false; - spec.rangeMin = minY; - spec.rangeMax = maxY; - spec.findMaximum = findMaximum; - rangeSpecs.push_back(spec); + if constexpr (D == 2) { + std::array minBounds = {std::numeric_limits::lowest(), minY}; + std::array maxBounds = {std::numeric_limits::max(), maxY}; + addRange(0, minBounds, maxBounds, findMaximum); + } else { + VIENNACORE_LOG_WARNING( + "addYRange is only supported for 2D. Use addRange for 3D."); + } } /// Clear all range specifications @@ -209,27 +211,6 @@ template class CompareCriticalDimensions { ToSurfaceMesh(levelSetTarget, surfaceMeshRef).apply(); ToSurfaceMesh(levelSetSample, surfaceMeshCmp).apply(); - // Get actual mesh extents instead of grid bounds - // This ensures we don't filter out surface points that extend beyond grid - // bounds - T xMin = std::numeric_limits::max(); - T xMax = std::numeric_limits::lowest(); - T yMin = std::numeric_limits::max(); - T yMax = std::numeric_limits::lowest(); - - for (const auto &node : surfaceMeshRef->nodes) { - xMin = std::min(xMin, node[0]); - xMax = std::max(xMax, node[0]); - yMin = std::min(yMin, node[1]); - yMax = std::max(yMax, node[1]); - } - for (const auto &node : surfaceMeshCmp->nodes) { - xMin = std::min(xMin, node[0]); - xMax = std::max(xMax, node[0]); - yMin = std::min(yMin, node[1]); - yMax = std::max(yMax, node[1]); - } - // Pre-allocate results vector to avoid race conditions results.resize(rangeSpecs.size()); @@ -238,33 +219,15 @@ template class CompareCriticalDimensions { for (size_t specIdx = 0; specIdx < rangeSpecs.size(); ++specIdx) { const auto &spec = rangeSpecs[specIdx]; CriticalDimensionResult result; - result.isXRange = spec.isXRange; - result.rangeMin = spec.rangeMin; - result.rangeMax = spec.rangeMax; + result.measureDimension = spec.measureDimension; + result.minBounds = spec.minBounds; + result.maxBounds = spec.maxBounds; result.findMaximum = spec.findMaximum; result.valid = false; - // Determine scan and perpendicular ranges - T scanMin, scanMax, perpMin, perpMax; - if (spec.isXRange) { - // X range specified - scan along X, find Y positions - scanMin = spec.rangeMin; - scanMax = spec.rangeMax; - perpMin = yMin; - perpMax = yMax; - } else { - // Y range specified - scan along Y, find X positions - scanMin = spec.rangeMin; - scanMax = spec.rangeMax; - perpMin = xMin; - perpMax = xMax; - } - // Find all surface crossings from the mesh nodes - auto crossingsRef = findSurfaceCrossings( - surfaceMeshRef, spec.isXRange, scanMin, scanMax, perpMin, perpMax); - auto crossingsCmp = findSurfaceCrossings( - surfaceMeshCmp, spec.isXRange, scanMin, scanMax, perpMin, perpMax); + auto crossingsRef = findSurfaceCrossings(surfaceMeshRef, spec); + auto crossingsCmp = findSurfaceCrossings(surfaceMeshCmp, spec); // Find critical dimensions auto [validRef, cdRef] = @@ -363,9 +326,11 @@ template class CompareCriticalDimensions { std::vector targetValues; std::vector sampleValues; - for (unsigned i = 0; i < D; ++i) { - outputMesh->minimumExtent[i] = std::numeric_limits::max(); - outputMesh->maximumExtent[i] = std::numeric_limits::lowest(); + for (unsigned i = 0; i < 3; ++i) { + outputMesh->minimumExtent[i] = + (i < D) ? std::numeric_limits::max() : 0.0; + outputMesh->maximumExtent[i] = + (i < D) ? std::numeric_limits::lowest() : 0.0; } unsigned pointId = 0; @@ -374,18 +339,21 @@ template class CompareCriticalDimensions { continue; // Create points for target and sample positions - Vec3D coordTarget, coordSample; - - if (result.isXRange) { - // Critical dimension is in Y, position is along X range - T xMid = (result.rangeMin + result.rangeMax) / 2.0; - coordTarget = {xMid, result.positionTarget, 0.0}; - coordSample = {xMid, result.positionSample, 0.0}; - } else { - // Critical dimension is in X, position is along Y range - T yMid = (result.rangeMin + result.rangeMax) / 2.0; - coordTarget = {result.positionTarget, yMid, 0.0}; - coordSample = {result.positionSample, yMid, 0.0}; + Vec3D coordTarget = {0.0, 0.0, 0.0}, coordSample = {0.0, 0.0, 0.0}; + + for (int i = 0; i < D; ++i) { + if (i == result.measureDimension) { + coordTarget[i] = result.positionTarget; + coordSample[i] = result.positionSample; + } else { + // Use midpoint of range for other dimensions + T mid = (result.minBounds[i] + result.maxBounds[i]) / 2.0; + // If bounds are infinite, use 0 + if (std::abs(mid) > std::numeric_limits::max() / 2.0) + mid = 0.0; + coordTarget[i] = mid; + coordSample[i] = mid; + } } // Add target point diff --git a/include/viennals/lsCompareNarrowBand.hpp b/include/viennals/lsCompareNarrowBand.hpp index c75f343f..6b952793 100644 --- a/include/viennals/lsCompareNarrowBand.hpp +++ b/include/viennals/lsCompareNarrowBand.hpp @@ -15,7 +15,7 @@ using namespace viennacore; /// Calculate distance measure between two level sets by comparing their SDF /// values on a narrow band. Returns the sum of squared differences between /// corresponding grid points. -/// The code is currently tended for 2D level sets only. +/// The code is intended for 2D and 3D level sets. template class CompareNarrowBand { using hrleIndexType = viennahrle::IndexType; SmartPointer> levelSetTarget = nullptr; @@ -27,8 +27,11 @@ template class CompareNarrowBand { T xRangeMax = std::numeric_limits::max(); T yRangeMin = std::numeric_limits::lowest(); T yRangeMax = std::numeric_limits::max(); + T zRangeMin = std::numeric_limits::lowest(); + T zRangeMax = std::numeric_limits::max(); bool useXRange = false; bool useYRange = false; + bool useZRange = false; // Fields to store the calculation results T sumSquaredDifferences = 0.0; @@ -135,20 +138,12 @@ template class CompareNarrowBand { } public: - CompareNarrowBand() { - assert( - D == 2 && - "CompareNarrowBand is currently only implemented for 2D level sets."); - } + CompareNarrowBand() {} CompareNarrowBand(SmartPointer> passedLevelSetTarget, SmartPointer> passedlevelSetSample) : levelSetTarget(passedLevelSetTarget), - levelSetSample(passedlevelSetSample) { - assert( - D == 2 && - "CompareNarrowBand is currently only implemented for 2D level sets."); - } + levelSetSample(passedlevelSetSample) {} void setLevelSetTarget(SmartPointer> passedLevelSet) { levelSetTarget = passedLevelSet; @@ -186,6 +181,20 @@ template class CompareNarrowBand { yRangeMax = std::numeric_limits::max(); } + /// Set the z-coordinate range to restrict the comparison area + void setZRange(T minZRange, T maxZRange) { + zRangeMin = minZRange; + zRangeMax = maxZRange; + useZRange = true; + } + + /// Clear the z-range restriction + void clearZRange() { + useZRange = false; + zRangeMin = std::numeric_limits::lowest(); + zRangeMax = std::numeric_limits::max(); + } + /// Set the output mesh where difference values will be stored void setOutputMesh(SmartPointer> passedMesh, bool outputMeshSquaredDiffs = true) { @@ -233,9 +242,11 @@ template class CompareNarrowBand { outputMesh->clear(); // Initialize mesh extent - for (unsigned i = 0; i < D; ++i) { - outputMesh->minimumExtent[i] = std::numeric_limits::max(); - outputMesh->maximumExtent[i] = std::numeric_limits::lowest(); + for (unsigned i = 0; i < 3; ++i) { + outputMesh->minimumExtent[i] = + (i < D) ? std::numeric_limits::max() : 0.0; + outputMesh->maximumExtent[i] = + (i < D) ? std::numeric_limits::lowest() : 0.0; } } @@ -244,6 +255,7 @@ template class CompareNarrowBand { // Check if current point is within specified x and y ranges T xCoord = itSample.getIndices()[0] * gridDelta; T yCoord = (D > 1) ? itSample.getIndices()[1] * gridDelta : 0; + T zCoord = (D > 2) ? itSample.getIndices()[2] * gridDelta : 0; // Skip if outside the specified x-range if (useXRange && (xCoord < xRangeMin || xCoord > xRangeMax)) { @@ -255,6 +267,11 @@ template class CompareNarrowBand { continue; } + // Skip if outside the specified z-range (only check in 3D) + if (D > 2 && useZRange && (zCoord < zRangeMin || zCoord > zRangeMax)) { + continue; + } + // Move the second iterator to the same position itTarget.goToIndicesSequential(itSample.getIndices()); diff --git a/include/viennals/lsCompareSparseField.hpp b/include/viennals/lsCompareSparseField.hpp index 0b012476..1041400b 100644 --- a/include/viennals/lsCompareSparseField.hpp +++ b/include/viennals/lsCompareSparseField.hpp @@ -29,7 +29,7 @@ using namespace viennacore; /// The iterated level set is expected to be sparse. The reduction is performed /// automatically if this is not the case. /// -/// The code is currently intended for 2D level sets only. +/// The code is intended for 2D and 3D level sets. template class CompareSparseField { using hrleIndexType = viennahrle::IndexType; @@ -42,8 +42,11 @@ template class CompareSparseField { T xRangeMax = std::numeric_limits::max(); T yRangeMin = std::numeric_limits::lowest(); T yRangeMax = std::numeric_limits::max(); + T zRangeMin = std::numeric_limits::lowest(); + T zRangeMax = std::numeric_limits::max(); bool useXRange = false; bool useYRange = false; + bool useZRange = false; // Fields to store the calculation results T sumSquaredDifferences = 0.0; @@ -176,6 +179,20 @@ template class CompareSparseField { yRangeMax = std::numeric_limits::max(); } + /// Set the z-coordinate range to restrict the comparison area + void setZRange(T minZRange, T maxZRange) { + zRangeMin = minZRange; + zRangeMax = maxZRange; + useZRange = true; + } + + /// Clear the z-range restriction + void clearZRange() { + useZRange = false; + zRangeMin = std::numeric_limits::lowest(); + zRangeMax = std::numeric_limits::max(); + } + /// Set the output mesh where difference values will be stored void setOutputMesh(SmartPointer> passedMesh) { outputMesh = passedMesh; @@ -231,9 +248,11 @@ template class CompareSparseField { outputMesh->clear(); // Initialize mesh extent - for (unsigned i = 0; i < D; ++i) { - outputMesh->minimumExtent[i] = std::numeric_limits::max(); - outputMesh->maximumExtent[i] = std::numeric_limits::lowest(); + for (unsigned i = 0; i < 3; ++i) { + outputMesh->minimumExtent[i] = + (i < D) ? std::numeric_limits::max() : 0.0; + outputMesh->maximumExtent[i] = + (i < D) ? std::numeric_limits::lowest() : 0.0; } // Reserve space for mesh data @@ -285,6 +304,12 @@ template class CompareSparseField { continue; } + // Skip if outside the specified z-range + if (D == 3 && useZRange && (zCoord < zRangeMin || zCoord > zRangeMax)) { + itIterated.next(); + continue; + } + // Get iterated value T valueIterated = itIterated.getValue(); diff --git a/include/viennals/lsCompareArea.hpp b/include/viennals/lsCompareVolume.hpp similarity index 76% rename from include/viennals/lsCompareArea.hpp rename to include/viennals/lsCompareVolume.hpp index 7115c223..f568ef0b 100644 --- a/include/viennals/lsCompareArea.hpp +++ b/include/viennals/lsCompareVolume.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include @@ -13,15 +14,15 @@ namespace viennals { using namespace viennacore; -/// Computes an estimate of the area where two level sets differ. -/// The area is calculated by iterating through the bounding box of the two +/// Computes an estimate of the volume/area where two level sets differ. +/// The volume is calculated by iterating through the bounding box of the two /// level sets and comparing the cell values. The grid delta is used as the unit -/// of area. Custom increment values can be set for specific x and y ranges, -/// allowing to count certain areas multiple times or skip them. Optionally, a -/// passed mesh can be filled with the area information, allowing for -/// visualization of the differences. -/// The code is currently itended for 2D level sets only. -template class CompareArea { +/// of volume. Custom increment values can be set for specific x, y and z +/// ranges, allowing to count certain areas multiple times or skip them. +/// Optionally, a passed mesh can be filled with the volume information, +/// allowing for visualization of the differences. The code is intended for 2D +/// and 3D level sets. +template class CompareVolume { using hrleDomainType = typename Domain::DomainType; using hrleIndexType = viennahrle::IndexType; @@ -36,11 +37,15 @@ template class CompareArea { hrleIndexType xRangeMax = std::numeric_limits::max(); hrleIndexType yRangeMin = std::numeric_limits::lowest(); hrleIndexType yRangeMax = std::numeric_limits::max(); + hrleIndexType zRangeMin = std::numeric_limits::lowest(); + hrleIndexType zRangeMax = std::numeric_limits::max(); bool useCustomXIncrement = false; bool useCustomYIncrement = false; + bool useCustomZIncrement = false; unsigned short int customXIncrement = 0; unsigned short int customYIncrement = 0; + unsigned short int customZIncrement = 0; unsigned short int defaultIncrement = 1; double gridDelta = 0.0; @@ -48,10 +53,17 @@ template class CompareArea { // Mesh output related members SmartPointer> outputMesh = nullptr; + // Helper to get the class name for logging + static const char *getClassName() { + if constexpr (D == 2) + return "CompareArea"; + return "CompareVolume"; + } + bool checkAndCalculateBounds() { if (levelSetTarget == nullptr || levelSetSample == nullptr) { Logger::getInstance() - .addError("Missing level set in CompareArea.") + .addError(std::string("Missing level set in ") + getClassName() + ".") .print(); return false; } @@ -62,8 +74,8 @@ template class CompareArea { if (gridTarget.getGridDelta() != gridSample.getGridDelta()) { Logger::getInstance() - .addError("Grid delta mismatch in CompareArea. The grid deltas of " - "the two level sets must be equal.") + .addError(std::string("Grid delta mismatch in ") + getClassName() + + ". The grid deltas of the two level sets must be equal.") .print(); return false; } else { @@ -103,18 +115,12 @@ template class CompareArea { } public: - CompareArea() { - assert(D == 2 && - "CompareArea is currently only implemented for 2D level sets."); - } + CompareVolume() {} - CompareArea(SmartPointer> passedLevelSetTarget, - SmartPointer> passedLevelSetSample) + CompareVolume(SmartPointer> passedLevelSetTarget, + SmartPointer> passedLevelSetSample) : levelSetTarget(passedLevelSetTarget), - levelSetSample(passedLevelSetSample) { - assert(D == 2 && - "CompareArea is currently only implemented for 2D level sets."); - } + levelSetSample(passedLevelSetSample) {} /// Sets the target level set. void setLevelSetTarget(SmartPointer> passedLevelSet) { @@ -149,6 +155,15 @@ template class CompareArea { useCustomYIncrement = true; } + /// Sets the z-range and custom increment value + void setZRangeAndIncrement(hrleIndexType minZRange, hrleIndexType maxZRange, + unsigned short int Zincrement) { + zRangeMin = minZRange; + zRangeMax = maxZRange; + customZIncrement = Zincrement; + useCustomZIncrement = true; + } + /// Set the output mesh where difference areas will be stored for /// visualization. Each cell in the mesh will have a cell data: /// 0: Areas where both level sets are inside @@ -157,17 +172,23 @@ template class CompareArea { outputMesh = passedMesh; } - /// Returns the computed area mismatch. - double getAreaMismatch() const { - return static_cast(differentCellsCount) * gridDelta * gridDelta; + /// Returns the computed volume/area mismatch. + double getVolumeMismatch() const { + return static_cast(differentCellsCount) * std::pow(gridDelta, D); } - /// Returns the computed area mismatch, with custom increments applied. - double getCustomAreaMismatch() const { - return static_cast(customDifferentCellCount) * gridDelta * - gridDelta; + /// Alias for getVolumeMismatch for 2D compatibility + double getAreaMismatch() const { return getVolumeMismatch(); } + + /// Returns the computed volume/area mismatch, with custom increments applied. + double getCustomVolumeMismatch() const { + return static_cast(customDifferentCellCount) * + std::pow(gridDelta, D); } + /// Alias for getCustomVolumeMismatch for 2D compatibility + double getCustomAreaMismatch() const { return getCustomVolumeMismatch(); } + /// Returns the number of cells where the level sets differ. unsigned long int getCellCount() const { return differentCellsCount; } @@ -177,7 +198,7 @@ template class CompareArea { return customDifferentCellCount; } - /// Computes the area difference between the two level sets. + /// Computes the volume/area difference between the two level sets. void apply() { // Calculate the bounds for iteration if (!checkAndCalculateBounds()) { @@ -200,7 +221,8 @@ template class CompareArea { if (levelSetTarget->getLevelSetWidth() < minimumWidth) { workingTarget = SmartPointer>::New(levelSetTarget); Expand(workingTarget, minimumWidth).apply(); - VIENNACORE_LOG_INFO("CompareArea: Expanded target level set to width " + + VIENNACORE_LOG_INFO(std::string(getClassName()) + + ": Expanded target level set to width " + std::to_string(minimumWidth) + " to avoid undefined values."); } @@ -208,7 +230,8 @@ template class CompareArea { if (levelSetSample->getLevelSetWidth() < minimumWidth) { workingSample = SmartPointer>::New(levelSetSample); Expand(workingSample, minimumWidth).apply(); - VIENNACORE_LOG_INFO("CompareArea: Expanded sample level set to width " + + VIENNACORE_LOG_INFO(std::string(getClassName()) + + ": Expanded sample level set to width " + std::to_string(minimumWidth) + " to avoid undefined values."); } @@ -227,9 +250,11 @@ template class CompareArea { if (generateMesh) { // Save the extent of the resulting mesh outputMesh->clear(); - for (unsigned i = 0; i < D; ++i) { - outputMesh->minimumExtent[i] = std::numeric_limits::max(); - outputMesh->maximumExtent[i] = std::numeric_limits::lowest(); + for (unsigned i = 0; i < 3; ++i) { + outputMesh->minimumExtent[i] = + (i < D) ? std::numeric_limits::max() : 0.0; + outputMesh->maximumExtent[i] = + (i < D) ? std::numeric_limits::lowest() : 0.0; } } @@ -270,16 +295,24 @@ template class CompareArea { bool inYRange = useCustomYIncrement && (itTarget.getIndices()[1] * gridDelta >= yRangeMin && itTarget.getIndices()[1] * gridDelta <= yRangeMax); + bool inZRange = + (D < 3) || (useCustomZIncrement && + (itTarget.getIndices()[2] * gridDelta >= zRangeMin && + itTarget.getIndices()[2] * gridDelta <= zRangeMax)); // Calculate increment to add based on ranges unsigned short int incrementToAdd = defaultIncrement; - if (inXRange && inYRange) { - // Apply both increments - incrementToAdd = customXIncrement + customYIncrement; - } else if (inXRange) { - incrementToAdd = customXIncrement; - } else if (inYRange) { - incrementToAdd = customYIncrement; + + // If any custom range is active, start from 0 and add the custom + // increments + if (inXRange || inYRange || (D == 3 && inZRange)) { + incrementToAdd = 0; + if (inXRange) + incrementToAdd += customXIncrement; + if (inYRange) + incrementToAdd += customYIncrement; + if (D == 3 && inZRange) + incrementToAdd += customZIncrement; } // If cells differ, update the counters @@ -350,7 +383,7 @@ template class CompareArea { // Insert points into the mesh outputMesh->nodes.resize(pointIdMapping.size()); for (auto it = pointIdMapping.begin(); it != pointIdMapping.end(); ++it) { - Vec3D coords; + Vec3D coords = {0.0, 0.0, 0.0}; for (unsigned i = 0; i < D; ++i) { coords[i] = gridDelta * it->first[i]; @@ -374,7 +407,10 @@ template class CompareArea { } }; +// Aliases for backward compatibility and clarity +template using CompareArea = CompareVolume; + // Precompile for common precision and dimensions -PRECOMPILE_PRECISION_DIMENSION(CompareArea) +PRECOMPILE_PRECISION_DIMENSION(CompareVolume) } // namespace viennals diff --git a/include/viennals/lsDetectFeatures.hpp b/include/viennals/lsDetectFeatures.hpp index 537fa616..a69de7ea 100644 --- a/include/viennals/lsDetectFeatures.hpp +++ b/include/viennals/lsDetectFeatures.hpp @@ -204,9 +204,12 @@ template class DetectFeatures { bool flag = false; - for (unsigned dir = 0; dir < (D * D * D); dir++) { - Vec3D currentNormal = - normals[neighborIt.getNeighbor(dir).getPointId()]; + constexpr unsigned numNeighbors = (D == 3) ? 27 : 9; + for (unsigned dir = 0; dir < numNeighbors; dir++) { + auto neighbor = neighborIt.getNeighbor(dir); + if (!neighbor.isDefined()) + continue; + Vec3D currentNormal = normals[neighbor.getPointId()]; if (currentNormal != zeroVector) { T skp = 0.; diff --git a/include/viennals/lsFromMesh.hpp b/include/viennals/lsFromMesh.hpp index 2d57cdac..fcbfeddf 100644 --- a/include/viennals/lsFromMesh.hpp +++ b/include/viennals/lsFromMesh.hpp @@ -53,7 +53,10 @@ template class FromMesh { auto &domain = levelSet->getDomain(); auto &nodes = mesh->getNodes(); - auto values = mesh->cellData.getScalarData("LSValues", true); + auto values = mesh->getPointData().getScalarData("LSValues", true); + if (values == nullptr) { + values = mesh->getCellData().getScalarData("LSValues", true); + } if (values == nullptr) { Logger::getInstance() diff --git a/include/viennals/lsMakeGeometry.hpp b/include/viennals/lsMakeGeometry.hpp index 50daf16a..2a45bbd3 100644 --- a/include/viennals/lsMakeGeometry.hpp +++ b/include/viennals/lsMakeGeometry.hpp @@ -419,115 +419,171 @@ template class MakeGeometry { } void makeCylinder(SmartPointer> cylinder) { - if (D != 3) { - Logger::getInstance() - .addError("MakeGeometry: Cylinder can only be created in 3D!") - .print(); - return; - } - // generate the points on the edges of the cylinders and mesh - // them manually - // cylinder axis will be (0,0,1) - auto gridDelta = levelSet->getGrid().getGridDelta(); - - auto points = PointCloud::New(); - const unsigned numPoints = - std::ceil(2 * M_PI * cylinder->radius / gridDelta); - const double smallAngle = 2.0 * M_PI / static_cast(numPoints); - - auto mesh = Mesh::New(); - // insert midpoint at base - mesh->insertNextNode(Vec3D{0.0, 0.0, 0.0}); - { - constexpr double limit = 2 * M_PI - 1e-6; - std::vector> points; - if (cylinder->topRadius) - std::vector> pointsTop; - - // create and insert points at base - for (double angle = 0.; angle < limit; angle += smallAngle) { - Vec3D point; - point[0] = cylinder->radius * std::cos(angle); - point[1] = cylinder->radius * std::sin(angle); - point[2] = 0.0; - points.push_back(point); - mesh->insertNextNode(point); - } - - // insert midpoint at top - mesh->insertNextNode(Vec3D{0.0, 0.0, cylinder->height}); - - double angle = 0; - for (unsigned i = 0; i < numPoints; ++i) { - // create triangles at base - std::array triangle{}; - triangle[0] = (i + 1) % numPoints + 1; - triangle[1] = i + 1; - triangle[2] = 0; - mesh->insertNextTriangle(triangle); + if constexpr (D == 3) { + // generate the points on the edges of the cylinders and mesh + // them manually + // cylinder axis will be (0,0,1) + auto gridDelta = levelSet->getGrid().getGridDelta(); + + auto points = PointCloud::New(); + const unsigned numPoints = + std::ceil(2 * M_PI * cylinder->radius / gridDelta); + const double smallAngle = 2.0 * M_PI / static_cast(numPoints); + + auto mesh = Mesh::New(); + // insert midpoint at base + mesh->insertNextNode(Vec3D{0.0, 0.0, 0.0}); + { + constexpr double limit = 2 * M_PI - 1e-6; + std::vector> points; + if (cylinder->topRadius) + std::vector> pointsTop; + + // create and insert points at base + for (double angle = 0.; angle < limit; angle += smallAngle) { + Vec3D point; + point[0] = cylinder->radius * std::cos(angle); + point[1] = cylinder->radius * std::sin(angle); + point[2] = 0.0; + points.push_back(point); + mesh->insertNextNode(point); + } - // insert points at top - // If topRadius is specified, update the first two coordinates of the - // points - if (cylinder->topRadius) { - points[i][0] = cylinder->topRadius * std::cos(angle); - points[i][1] = cylinder->topRadius * std::sin(angle); - angle += smallAngle; + // insert midpoint at top + mesh->insertNextNode(Vec3D{0.0, 0.0, cylinder->height}); + + double angle = 0; + for (unsigned i = 0; i < numPoints; ++i) { + // create triangles at base + std::array triangle{}; + triangle[0] = (i + 1) % numPoints + 1; + triangle[1] = i + 1; + triangle[2] = 0; + mesh->insertNextTriangle(triangle); + + // insert points at top + // If topRadius is specified, update the first two coordinates of the + // points + if (cylinder->topRadius) { + points[i][0] = cylinder->topRadius * std::cos(angle); + points[i][1] = cylinder->topRadius * std::sin(angle); + angle += smallAngle; + } + points[i][2] = cylinder->height; + mesh->insertNextNode(points[i]); + + // insert triangles at top + triangle[0] = numPoints + 1; + triangle[1] = numPoints + i + 2; + triangle[2] = (i + 1) % numPoints + 2 + numPoints; + mesh->insertNextTriangle(triangle); } - points[i][2] = cylinder->height; - mesh->insertNextNode(points[i]); - // insert triangles at top - triangle[0] = numPoints + 1; - triangle[1] = numPoints + i + 2; - triangle[2] = (i + 1) % numPoints + 2 + numPoints; - mesh->insertNextTriangle(triangle); + // insert sidewall triangles + for (unsigned i = 0; i < numPoints; ++i) { + std::array triangle{}; + triangle[0] = i + 1; + triangle[1] = (i + 1) % numPoints + 1; + triangle[2] = i + numPoints + 2; + mesh->insertNextTriangle(triangle); + + triangle[0] = (i + 1) % numPoints + 1; + triangle[1] = (i + 1) % numPoints + 2 + numPoints; + triangle[2] = i + numPoints + 2; + mesh->insertNextTriangle(triangle); + } } - // insert sidewall triangles - for (unsigned i = 0; i < numPoints; ++i) { - std::array triangle{}; - triangle[0] = i + 1; - triangle[1] = (i + 1) % numPoints + 1; - triangle[2] = i + numPoints + 2; - mesh->insertNextTriangle(triangle); - - triangle[0] = (i + 1) % numPoints + 1; - triangle[1] = (i + 1) % numPoints + 2 + numPoints; - triangle[2] = i + numPoints + 2; - mesh->insertNextTriangle(triangle); + // rotate mesh + // normalise axis vector + T unit = std::sqrt( + DotProduct(cylinder->axisDirection, cylinder->axisDirection)); + Vec3D cylinderAxis; + for (int i = 0; i < 3; ++i) { + cylinderAxis[i] = cylinder->axisDirection[i] / unit; + } + // get rotation axis via cross product of (0,0,1) and axis of cylinder + Vec3D rotAxis = {-cylinderAxis[1], cylinderAxis[0], 0.0}; + // angle is acos of dot product + T rotationAngle = std::acos(cylinderAxis[2]); + + // rotate mesh + TransformMesh(mesh, TransformEnum::ROTATION, rotAxis, rotationAngle) + .apply(); + + // translate mesh + Vec3D translationVector; + for (int i = 0; i < 3; ++i) { + translationVector[i] = cylinder->origin[i]; + } + TransformMesh(mesh, TransformEnum::TRANSLATION, translationVector) + .apply(); + + // read mesh from surface + FromSurfaceMesh mesher(levelSet, mesh); + mesher.setRemoveBoundaryTriangles(ignoreBoundaryConditions); + mesher.apply(); + } else if constexpr (D == 2) { + VIENNACORE_LOG_WARNING( + "MakeGeometry: Cylinder in 2D creates a box, not a cylinder."); + + auto mesh = Mesh::New(); + + // Calculate axis length for normalization + T axisLen = 0; + for (unsigned i = 0; i < D; ++i) + axisLen += cylinder->axisDirection[i] * cylinder->axisDirection[i]; + axisLen = std::sqrt(axisLen); + + // Normalized axis and perpendicular vector + // In 2D, if axis is (u, v), perpendicular is (v, -u) for CCW winding + T axis[2] = {cylinder->axisDirection[0] / axisLen, + cylinder->axisDirection[1] / axisLen}; + T perp[2] = {axis[1], -axis[0]}; + + T topRadius = + (cylinder->topRadius != 0.) ? cylinder->topRadius : cylinder->radius; + + // Define the 4 corners of the rectangle/trapezoid + std::array, 4> corners; + + // Base Right + corners[0][0] = cylinder->origin[0] + cylinder->radius * perp[0]; + corners[0][1] = cylinder->origin[1] + cylinder->radius * perp[1]; + corners[0][2] = 0.; + + // Top Right + corners[1][0] = cylinder->origin[0] + cylinder->height * axis[0] + + topRadius * perp[0]; + corners[1][1] = cylinder->origin[1] + cylinder->height * axis[1] + + topRadius * perp[1]; + corners[1][2] = 0.; + + // Top Left + corners[2][0] = cylinder->origin[0] + cylinder->height * axis[0] - + topRadius * perp[0]; + corners[2][1] = cylinder->origin[1] + cylinder->height * axis[1] - + topRadius * perp[1]; + corners[2][2] = 0.; + + // Base Left + corners[3][0] = cylinder->origin[0] - cylinder->radius * perp[0]; + corners[3][1] = cylinder->origin[1] - cylinder->radius * perp[1]; + corners[3][2] = 0.; + + for (const auto &p : corners) { + mesh->insertNextNode(p); } - } - // rotate mesh - // normalise axis vector - T unit = - std::sqrt(DotProduct(cylinder->axisDirection, cylinder->axisDirection)); - Vec3D cylinderAxis; - for (int i = 0; i < 3; ++i) { - cylinderAxis[i] = cylinder->axisDirection[i] / unit; - } - // get rotation axis via cross product of (0,0,1) and axis of cylinder - Vec3D rotAxis = {-cylinderAxis[1], cylinderAxis[0], 0.0}; - // angle is acos of dot product - T rotationAngle = std::acos(cylinderAxis[2]); + mesh->insertNextLine({0, 1}); + mesh->insertNextLine({1, 2}); + mesh->insertNextLine({2, 3}); + mesh->insertNextLine({3, 0}); - // rotate mesh - TransformMesh(mesh, TransformEnum::ROTATION, rotAxis, rotationAngle) - .apply(); - - // translate mesh - Vec3D translationVector; - for (int i = 0; i < 3; ++i) { - translationVector[i] = cylinder->origin[i]; + FromSurfaceMesh mesher(levelSet, mesh); + mesher.setRemoveBoundaryTriangles(ignoreBoundaryConditions); + mesher.apply(); } - TransformMesh(mesh, TransformEnum::TRANSLATION, translationVector) - .apply(); - - // read mesh from surface - FromSurfaceMesh mesher(levelSet, mesh); - mesher.setRemoveBoundaryTriangles(ignoreBoundaryConditions); - mesher.apply(); } void makeCustom(SmartPointer> pointCloud) { diff --git a/include/viennals/lsMarkVoidPoints.hpp b/include/viennals/lsMarkVoidPoints.hpp index d2a2484d..92e1753e 100644 --- a/include/viennals/lsMarkVoidPoints.hpp +++ b/include/viennals/lsMarkVoidPoints.hpp @@ -317,7 +317,7 @@ template class MarkVoidPoints { auto maxComponent = *std::max_element(componentMarkers.begin(), componentMarkers.end()); assert(maxComponent >= 0); - numComponents = maxComponent; + numComponents = maxComponent + 1; auto componentMarkersPointer = pointData.getScalarData("ConnectedComponentId", true); diff --git a/lib/specialisations.cpp b/lib/specialisations.cpp index 32c4eb4e..3e82cdcb 100644 --- a/lib/specialisations.cpp +++ b/lib/specialisations.cpp @@ -9,9 +9,10 @@ #include #include #include -#include +#include #include #include +#include #include #include #include @@ -46,7 +47,8 @@ PRECOMPILE_SPECIALIZE(CalculateCurvatures) PRECOMPILE_SPECIALIZE(CalculateNormalVectors) PRECOMPILE_SPECIALIZE(Check) PRECOMPILE_SPECIALIZE(ConvexHull) -PRECOMPILE_SPECIALIZE(CompareArea) +PRECOMPILE_SPECIALIZE(CompareChamfer) +PRECOMPILE_SPECIALIZE(CompareVolume) PRECOMPILE_SPECIALIZE(CompareNarrowBand) PRECOMPILE_SPECIALIZE(CompareSparseField) PRECOMPILE_SPECIALIZE(Domain) diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 192e67c1..dd14b5d0 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -77,7 +77,6 @@ add_dependencies(${PROJECT_NAME} ${MOD}) set(VIENNALS_LIB_FOLDER ${CMAKE_BINARY_DIR}/${PKG}.libs) -# Not required for both targets, one will suffice viennacore_setup_vtk_env(${MOD} ${VIENNALS_LIB_FOLDER}) install( diff --git a/python/pyWrap.cpp b/python/pyWrap.cpp index 844c7fbe..f231b5c7 100644 --- a/python/pyWrap.cpp +++ b/python/pyWrap.cpp @@ -231,9 +231,11 @@ PYBIND11_MODULE(VIENNALS_MODULE_NAME, module) { "Get a list of hexahedrons of the mesh.") .def("getPointData", (PointData & (Mesh::*)()) & Mesh::getPointData, + py::return_value_policy::reference_internal, "Return a reference to the point data of the mesh.") .def("getCellData", (PointData & (Mesh::*)()) & Mesh::getCellData, + py::return_value_policy::reference_internal, "Return a reference to the cell data of the mesh.") .def("insertNextNode", &Mesh::insertNextNode, "Insert a node in the mesh.") diff --git a/python/pyWrap.hpp b/python/pyWrap.hpp index bd56a936..f6582d83 100644 --- a/python/pyWrap.hpp +++ b/python/pyWrap.hpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -9,11 +10,11 @@ #include #include #include -#include #include #include #include #include +#include #include #include #include @@ -232,6 +233,10 @@ template void bindApi(py::module &module) { .def("setUpdatePointData", &Advect::setUpdatePointData, "Set whether the point data in the old LS should be translated to " "the advected LS. Defaults to true.") + .def( + "setVelocityUpdateCallback", &Advect::setVelocityUpdateCallback, + "Set a callback function that is called after the level set has been " + "updated during intermediate time integration steps (e.g. RK2, RK3).") .def("prepareLS", &Advect::prepareLS, "Prepare the level-set.") // need scoped release since we are calling a python method from // parallelised C++ code here @@ -648,7 +653,7 @@ template void bindApi(py::module &module) { .def("setVoidTopSurface", &MarkVoidPoints::setVoidTopSurface, "Set the logic by which to choose the surface which is non-void. " "All other connected surfaces will then be marked as void points.") - .def("setSaveComponentsId", &MarkVoidPoints::setSaveComponentIds, + .def("setSaveComponentIds", &MarkVoidPoints::setSaveComponentIds, "Save the connectivity information of all LS points in the " "pointData of the level set.") .def("getNumberOfComponents", @@ -843,6 +848,58 @@ template void bindApi(py::module &module) { "Set the filename for the output file.") .def("apply", &Writer::apply, "Write to file."); + // CompareSparseField + py::class_, SmartPointer>>( + module, "CompareSparseField") + // constructors + .def(py::init(&SmartPointer>::template New<>)) + .def( + py::init(&SmartPointer>::template New< + SmartPointer> &, SmartPointer> &>)) + // methods + .def("setLevelSetExpanded", + &CompareSparseField::setLevelSetExpanded, + "Sets the expanded level set for comparison.") + .def("setLevelSetIterated", + &CompareSparseField::setLevelSetIterated, + "Sets the iterated level set to compare against the expanded one.") + .def("setXRange", &CompareSparseField::setXRange, + "Set the x-coordinate range to restrict the comparison area") + .def("setYRange", &CompareSparseField::setYRange, + "Set the y-coordinate range to restrict the comparison area") + .def("clearXRange", &CompareSparseField::clearXRange, + "Clear the x-range restriction") + .def("clearYRange", &CompareSparseField::clearYRange, + "Clear the y-range restriction") + .def("setZRange", &CompareSparseField::setZRange, + "Set the z-coordinate range to restrict the comparison area") + .def("clearZRange", &CompareSparseField::clearZRange, + "Clear the z-range restriction") + .def("setOutputMesh", &CompareSparseField::setOutputMesh, + "Set the output mesh where difference values will be stored") + .def("setFillIteratedWithDistances", + &CompareSparseField::setFillIteratedWithDistances, + "Set whether to fill the iterated level set with distance values") + .def("setExpandedLevelSetWidth", + &CompareSparseField::setExpandedLevelSetWidth, + "Set the expansion width for the expanded level set") + .def("apply", &CompareSparseField::apply, + "Apply the comparison and calculate the sum of squared " + "differences.") + .def("getSumSquaredDifferences", + &CompareSparseField::getSumSquaredDifferences, + "Return the sum of squared differences calculated by apply().") + .def("getSumDifferences", &CompareSparseField::getSumDifferences, + "Return the sum of absolute differences calculated by apply().") + .def("getNumPoints", &CompareSparseField::getNumPoints, + "Return the number of points used in the comparison.") + .def("getNumSkippedPoints", + &CompareSparseField::getNumSkippedPoints, + "Return the number of points skipped during comparison.") + .def("getRMSE", &CompareSparseField::getRMSE, + "Calculate the root mean square error from previously computed " + "values."); + // WriteVisualizationMesh #ifdef VIENNALS_USE_VTK py::class_, @@ -885,227 +942,196 @@ template void bindApi(py::module &module) { "Make and write mesh."); #endif - if constexpr (D == 2) { - // CompareArea - py::class_, SmartPointer>>( - module, "CompareArea") - // constructors - .def(py::init(&SmartPointer>::template New<>)) - .def(py::init( - &SmartPointer>::template New< - SmartPointer> &, SmartPointer> &>)) - // methods - .def("setLevelSetTarget", &CompareArea::setLevelSetTarget, - "Sets the target level set.") - .def("setLevelSetSample", &CompareArea::setLevelSetSample, - "Sets the sample level set.") - .def("setDefaultIncrement", &CompareArea::setDefaultIncrement, - "Set default increment value") - .def("setXRangeAndIncrement", &CompareArea::setXRangeAndIncrement, - "Sets the x-range and custom increment value") - .def("setYRangeAndIncrement", &CompareArea::setYRangeAndIncrement, - "Sets the y-range and custom increment value") - .def("setOutputMesh", &CompareArea::setOutputMesh, - "Set the output mesh where difference areas will be stored") - .def("getAreaMismatch", &CompareArea::getAreaMismatch, - "Returns the computed area mismatch.") - .def("getCustomAreaMismatch", &CompareArea::getCustomAreaMismatch, - "Returns the computed area mismatch, with custom increments " - "applied.") - .def("getCellCount", &CompareArea::getCellCount, - "Returns the number of cells where the level sets differ.") - .def("getCustomCellCount", &CompareArea::getCustomCellCount, - "Returns the number of cells where the level sets differ, with " - "custom increments applied.") - .def("apply", &CompareArea::apply, - "Computes the area difference between the two level sets."); + // CompareVolume + py::class_, SmartPointer>>( + module, "CompareVolume") + // constructors + .def(py::init(&SmartPointer>::template New<>)) + .def( + py::init(&SmartPointer>::template New< + SmartPointer> &, SmartPointer> &>)) + // methods + .def("setLevelSetTarget", &CompareVolume::setLevelSetTarget, + "Sets the target level set.") + .def("setLevelSetSample", &CompareVolume::setLevelSetSample, + "Sets the sample level set.") + .def("setDefaultIncrement", &CompareVolume::setDefaultIncrement, + "Set default increment value") + .def("setXRangeAndIncrement", &CompareVolume::setXRangeAndIncrement, + "Sets the x-range and custom increment value") + .def("setYRangeAndIncrement", &CompareVolume::setYRangeAndIncrement, + "Sets the y-range and custom increment value") + .def("setZRangeAndIncrement", &CompareVolume::setZRangeAndIncrement, + "Sets the z-range and custom increment value") + .def("setOutputMesh", &CompareVolume::setOutputMesh, + "Set the output mesh where difference areas will be stored") + .def("getVolumeMismatch", &CompareVolume::getVolumeMismatch, + "Returns the computed volume mismatch.") + .def("getAreaMismatch", &CompareVolume::getAreaMismatch, + "Returns the computed area mismatch.") + .def("getCustomVolumeMismatch", + &CompareVolume::getCustomVolumeMismatch, + "Returns the computed volume mismatch, with custom increments " + "applied.") + .def("getCustomAreaMismatch", &CompareVolume::getCustomAreaMismatch, + "Returns the computed area mismatch, with custom increments " + "applied.") + .def("getCellCount", &CompareVolume::getCellCount, + "Returns the number of cells where the level sets differ.") + .def("getCustomCellCount", &CompareVolume::getCustomCellCount, + "Returns the number of cells where the level sets differ, with " + "custom increments applied.") + .def("apply", &CompareVolume::apply, + "Computes the volume difference between the two level sets."); - // CompareChamfer - py::class_, SmartPointer>>( - module, "CompareChamfer") - // constructors - .def(py::init(&SmartPointer>::template New<>)) - .def(py::init( - &SmartPointer>::template New< - SmartPointer> &, SmartPointer> &>)) - // methods - .def("setLevelSetTarget", &CompareChamfer::setLevelSetTarget, - "Set the target level set.") - .def("setLevelSetSample", &CompareChamfer::setLevelSetSample, - "Set the sample level set.") - .def("setOutputMeshTarget", &CompareChamfer::setOutputMeshTarget, - "Set output mesh for target surface points with distance data.") - .def("setOutputMeshSample", &CompareChamfer::setOutputMeshSample, - "Set output mesh for sample surface points with distance data.") - .def("apply", &CompareChamfer::apply, - "Apply the Chamfer distance calculation.") - .def("getForwardDistance", &CompareChamfer::getForwardDistance, - "Get the forward distance (average distance from target to " - "sample).") - .def("getBackwardDistance", &CompareChamfer::getBackwardDistance, - "Get the backward distance (average distance from sample to " - "target).") - .def("getChamferDistance", &CompareChamfer::getChamferDistance, - "Get the Chamfer distance (average of forward and backward).") - .def("getRMSChamferDistance", - &CompareChamfer::getRMSChamferDistance, - "Get the RMS Chamfer distance.") - .def("getMaxDistance", &CompareChamfer::getMaxDistance, - "Get the maximum nearest-neighbor distance.") - .def("getNumTargetPoints", &CompareChamfer::getNumTargetPoints, - "Get the number of target surface points.") - .def("getNumSamplePoints", &CompareChamfer::getNumSamplePoints, - "Get the number of sample surface points."); + if constexpr (D == 2) { + module.attr("CompareArea") = module.attr("CompareVolume"); + } - // CompareCriticalDimensions - py::class_, - SmartPointer>>( - module, "CompareCriticalDimensions") - // constructors - .def(py::init( - &SmartPointer>::template New<>)) - .def(py::init( - &SmartPointer>::template New< - SmartPointer> &, SmartPointer> &>)) - // methods - .def("setLevelSetTarget", - &CompareCriticalDimensions::setLevelSetTarget, - "Sets the target level set.") - .def("setLevelSetSample", - &CompareCriticalDimensions::setLevelSetSample, - "Sets the sample level set.") - .def("addXRange", &CompareCriticalDimensions::addXRange, - py::arg("minX"), py::arg("maxX"), py::arg("findMaximum") = true, - "Add an X range to find maximum or minimum Y position.") - .def("addYRange", &CompareCriticalDimensions::addYRange, - py::arg("minY"), py::arg("maxY"), py::arg("findMaximum") = true, - "Add a Y range to find maximum or minimum X position.") - .def("clearRanges", &CompareCriticalDimensions::clearRanges, - "Clear all range specifications.") - .def("setOutputMesh", &CompareCriticalDimensions::setOutputMesh, - "Set the output mesh where critical dimension locations will be " - "stored.") - .def("apply", &CompareCriticalDimensions::apply, - "Apply the comparison.") - .def("getNumCriticalDimensions", - &CompareCriticalDimensions::getNumCriticalDimensions, - "Get the number of critical dimensions compared.") - .def( - "getCriticalDimensionResult", - [](CompareCriticalDimensions &self, size_t index) { - T posRef, posCmp, diff; - bool valid = - self.getCriticalDimensionResult(index, posRef, posCmp, diff); - if (valid) { - return py::make_tuple(true, posRef, posCmp, diff); - } else { - return py::make_tuple(false, 0.0, 0.0, 0.0); - } - }, - py::arg("index"), - "Get a specific critical dimension result. Returns (valid, " - "positionTarget, positionSample, difference).") - .def("getMeanDifference", - &CompareCriticalDimensions::getMeanDifference, - "Get mean absolute difference across all valid critical " - "dimensions.") - .def("getMaxDifference", - &CompareCriticalDimensions::getMaxDifference, - "Get maximum difference across all valid critical dimensions.") - .def("getRMSE", &CompareCriticalDimensions::getRMSE, - "Get RMSE across all valid critical dimensions.") - .def("getAllDifferences", - &CompareCriticalDimensions::getAllDifferences, - "Get all valid differences as a list."); + // CompareChamfer + py::class_, SmartPointer>>( + module, "CompareChamfer") + // constructors + .def(py::init(&SmartPointer>::template New<>)) + .def( + py::init(&SmartPointer>::template New< + SmartPointer> &, SmartPointer> &>)) + // methods + .def("setLevelSetTarget", &CompareChamfer::setLevelSetTarget, + "Set the target level set.") + .def("setLevelSetSample", &CompareChamfer::setLevelSetSample, + "Set the sample level set.") + .def("setOutputMeshTarget", &CompareChamfer::setOutputMeshTarget, + "Set output mesh for target surface points with distance data.") + .def("setOutputMeshSample", &CompareChamfer::setOutputMeshSample, + "Set output mesh for sample surface points with distance data.") + .def("apply", &CompareChamfer::apply, + "Apply the Chamfer distance calculation.") + .def("getForwardDistance", &CompareChamfer::getForwardDistance, + "Get the forward distance (average distance from target to " + "sample).") + .def("getBackwardDistance", &CompareChamfer::getBackwardDistance, + "Get the backward distance (average distance from sample to " + "target).") + .def("getChamferDistance", &CompareChamfer::getChamferDistance, + "Get the Chamfer distance (average of forward and backward).") + .def("getRMSChamferDistance", + &CompareChamfer::getRMSChamferDistance, + "Get the RMS Chamfer distance.") + .def("getMaxDistance", &CompareChamfer::getMaxDistance, + "Get the maximum nearest-neighbor distance.") + .def("getNumTargetPoints", &CompareChamfer::getNumTargetPoints, + "Get the number of target surface points.") + .def("getNumSamplePoints", &CompareChamfer::getNumSamplePoints, + "Get the number of sample surface points."); - // CompareNarrowBand - py::class_, SmartPointer>>( - module, "CompareNarrowBand") - // constructors - .def(py::init(&SmartPointer>::template New<>)) - .def(py::init( - &SmartPointer>::template New< - SmartPointer> &, SmartPointer> &>)) - // methods - .def("setLevelSetTarget", &CompareNarrowBand::setLevelSetTarget, - "Sets the target level set.") - .def("setLevelSetSample", &CompareNarrowBand::setLevelSetSample, - "Sets the sample level set.") - .def("setXRange", &CompareNarrowBand::setXRange, - "Set the x-coordinate range to restrict the comparison area") - .def("setYRange", &CompareNarrowBand::setYRange, - "Set the y-coordinate range to restrict the comparison area") - .def("clearXRange", &CompareNarrowBand::clearXRange, - "Clear the x-range restriction") - .def("clearYRange", &CompareNarrowBand::clearYRange, - "Clear the y-range restriction") - .def("setOutputMesh", &CompareNarrowBand::setOutputMesh, - "Set the output mesh where difference values will be stored") - .def("setOutputMeshSquaredDifferences", - &CompareNarrowBand::setOutputMeshSquaredDifferences, - "Set whether to output squared differences (true) or absolute " - "differences (false)") - .def("apply", &CompareNarrowBand::apply, - "Apply the comparison and calculate the sum of squared " - "differences.") - .def("getSumSquaredDifferences", - &CompareNarrowBand::getSumSquaredDifferences, - "Return the sum of squared differences calculated by apply().") - .def("getSumDifferences", &CompareNarrowBand::getSumDifferences, - "Return the sum of absolute differences calculated by apply().") - .def("getNumPoints", &CompareNarrowBand::getNumPoints, - "Return the number of points used in the comparison.") - .def("getRMSE", &CompareNarrowBand::getRMSE, - "Calculate the root mean square error from previously computed " - "values."); + // CompareCriticalDimensions + py::class_, + SmartPointer>>( + module, "CompareCriticalDimensions") + // constructors + .def(py::init( + &SmartPointer>::template New<>)) + .def( + py::init(&SmartPointer>::template New< + SmartPointer> &, SmartPointer> &>)) + // methods + .def("setLevelSetTarget", + &CompareCriticalDimensions::setLevelSetTarget, + "Sets the target level set.") + .def("setLevelSetSample", + &CompareCriticalDimensions::setLevelSetSample, + "Sets the sample level set.") + .def("addRange", &CompareCriticalDimensions::addRange, + py::arg("measureDimension"), py::arg("minBounds"), + py::arg("maxBounds"), py::arg("findMaximum") = true, + "Add a generic range specification.") + .def("addXRange", &CompareCriticalDimensions::addXRange, + py::arg("minX"), py::arg("maxX"), py::arg("findMaximum") = true, + "Add an X range to find maximum or minimum Y position.") + .def("addYRange", &CompareCriticalDimensions::addYRange, + py::arg("minY"), py::arg("maxY"), py::arg("findMaximum") = true, + "Add a Y range to find maximum or minimum X position.") + .def("clearRanges", &CompareCriticalDimensions::clearRanges, + "Clear all range specifications.") + .def("setOutputMesh", &CompareCriticalDimensions::setOutputMesh, + "Set the output mesh where critical dimension locations will be " + "stored.") + .def("apply", &CompareCriticalDimensions::apply, + "Apply the comparison.") + .def("getNumCriticalDimensions", + &CompareCriticalDimensions::getNumCriticalDimensions, + "Get the number of critical dimensions compared.") + .def( + "getCriticalDimensionResult", + [](CompareCriticalDimensions &self, size_t index) { + T posRef, posCmp, diff; + bool valid = + self.getCriticalDimensionResult(index, posRef, posCmp, diff); + if (valid) { + return py::make_tuple(true, posRef, posCmp, diff); + } else { + return py::make_tuple(false, 0.0, 0.0, 0.0); + } + }, + py::arg("index"), + "Get a specific critical dimension result. Returns (valid, " + "positionTarget, positionSample, difference).") + .def("getMeanDifference", + &CompareCriticalDimensions::getMeanDifference, + "Get mean absolute difference across all valid critical " + "dimensions.") + .def("getMaxDifference", + &CompareCriticalDimensions::getMaxDifference, + "Get maximum difference across all valid critical dimensions.") + .def("getRMSE", &CompareCriticalDimensions::getRMSE, + "Get RMSE across all valid critical dimensions.") + .def("getAllDifferences", + &CompareCriticalDimensions::getAllDifferences, + "Get all valid differences as a list."); - // CompareSparseField - py::class_, - SmartPointer>>(module, - "CompareSparseField") - // constructors - .def(py::init(&SmartPointer>::template New<>)) - .def(py::init( - &SmartPointer>::template New< - SmartPointer> &, SmartPointer> &>)) - // methods - .def("setLevelSetExpanded", - &CompareSparseField::setLevelSetExpanded, - "Sets the expanded level set for comparison.") - .def("setLevelSetIterated", - &CompareSparseField::setLevelSetIterated, - "Sets the iterated level set to compare against the expanded one.") - .def("setXRange", &CompareSparseField::setXRange, - "Set the x-coordinate range to restrict the comparison area") - .def("setYRange", &CompareSparseField::setYRange, - "Set the y-coordinate range to restrict the comparison area") - .def("clearXRange", &CompareSparseField::clearXRange, - "Clear the x-range restriction") - .def("clearYRange", &CompareSparseField::clearYRange, - "Clear the y-range restriction") - .def("setOutputMesh", &CompareSparseField::setOutputMesh, - "Set the output mesh where difference values will be stored") - .def("setFillIteratedWithDistances", - &CompareSparseField::setFillIteratedWithDistances, - "Set whether to fill the iterated level set with distance values") - .def("setExpandedLevelSetWidth", - &CompareSparseField::setExpandedLevelSetWidth, - "Set the expansion width for the expanded level set") - .def("apply", &CompareSparseField::apply, - "Apply the comparison and calculate the sum of squared " - "differences.") - .def("getSumSquaredDifferences", - &CompareSparseField::getSumSquaredDifferences, - "Return the sum of squared differences calculated by apply().") - .def("getSumDifferences", &CompareSparseField::getSumDifferences, - "Return the sum of absolute differences calculated by apply().") - .def("getNumPoints", &CompareSparseField::getNumPoints, - "Return the number of points used in the comparison.") - .def("getNumSkippedPoints", - &CompareSparseField::getNumSkippedPoints, - "Return the number of points skipped during comparison.") - .def("getRMSE", &CompareSparseField::getRMSE, - "Calculate the root mean square error from previously computed " - "values."); - } + // CompareNarrowBand + py::class_, SmartPointer>>( + module, "CompareNarrowBand") + // constructors + .def(py::init(&SmartPointer>::template New<>)) + .def( + py::init(&SmartPointer>::template New< + SmartPointer> &, SmartPointer> &>)) + // methods + .def("setLevelSetTarget", &CompareNarrowBand::setLevelSetTarget, + "Sets the target level set.") + .def("setLevelSetSample", &CompareNarrowBand::setLevelSetSample, + "Sets the sample level set.") + .def("setXRange", &CompareNarrowBand::setXRange, + "Set the x-coordinate range to restrict the comparison area") + .def("setYRange", &CompareNarrowBand::setYRange, + "Set the y-coordinate range to restrict the comparison area") + .def("clearXRange", &CompareNarrowBand::clearXRange, + "Clear the x-range restriction") + .def("clearYRange", &CompareNarrowBand::clearYRange, + "Clear the y-range restriction") + .def("setZRange", &CompareNarrowBand::setZRange, + "Set the z-coordinate range to restrict the comparison area") + .def("clearZRange", &CompareNarrowBand::clearZRange, + "Clear the z-range restriction") + .def("setOutputMesh", &CompareNarrowBand::setOutputMesh, + "Set the output mesh where difference values will be stored") + .def("setOutputMeshSquaredDifferences", + &CompareNarrowBand::setOutputMeshSquaredDifferences, + "Set whether to output squared differences (true) or absolute " + "differences (false)") + .def("apply", &CompareNarrowBand::apply, + "Apply the comparison and calculate the sum of squared " + "differences.") + .def("getSumSquaredDifferences", + &CompareNarrowBand::getSumSquaredDifferences, + "Return the sum of squared differences calculated by apply().") + .def("getSumDifferences", &CompareNarrowBand::getSumDifferences, + "Return the sum of absolute differences calculated by apply().") + .def("getNumPoints", &CompareNarrowBand::getNumPoints, + "Return the number of points used in the comparison.") + .def("getRMSE", &CompareNarrowBand::getRMSE, + "Calculate the root mean square error from previously computed " + "values."); } diff --git a/python/tests/run_all_tests.py b/python/tests/run_all_tests.py new file mode 100644 index 00000000..be29aec2 --- /dev/null +++ b/python/tests/run_all_tests.py @@ -0,0 +1,810 @@ +import unittest +import os +import viennals as vls +import viennals.d2 as d2 +import viennals.d3 as d3 + +class TestVelocityField(vls.VelocityField): + def __init__(self): + super().__init__() + + def getScalarVelocity(self, coord, material, normal, pointId): + return 1.0 + + def getVectorVelocity(self, coord, material, normal, pointId): + return [0.0, 0.0, 0.0] + +class WrappingTest(unittest.TestCase): + def test_enums_and_constants(self): + print("\n[TEST] test_enums_and_constants") + # Verify Enums exist + self.assertTrue(hasattr(vls, "SpatialSchemeEnum")) + self.assertTrue(hasattr(vls, "TemporalSchemeEnum")) + self.assertTrue(hasattr(vls, "BooleanOperationEnum")) + self.assertTrue(hasattr(vls, "CurvatureEnum")) + self.assertTrue(hasattr(vls, "FeatureDetectionEnum")) + self.assertTrue(hasattr(vls, "BoundaryConditionEnum")) + self.assertTrue(hasattr(vls, "FileFormatEnum")) + self.assertTrue(hasattr(vls, "VoidTopSurfaceEnum")) + self.assertTrue(hasattr(vls, "TransformEnum")) + + # Verify Logger + logger = vls.Logger.getInstance() + logger.setLogLevel(vls.LogLevel.INFO) + logger.addInfo("Testing Logger") + print(" > Done test_enums_and_constants") + + def test_mesh_and_data(self): + print("\n[TEST] test_mesh_and_data") + # PointData + print(" > Testing PointData scalar insertion") + pd = vls.PointData() + pd.insertNextScalarData([1.0, 2.0], "Scalars") + self.assertEqual(pd.getScalarDataSize(), 1) + + # Mesh + print(" > Testing Mesh insertion") + mesh = vls.Mesh() + mesh.insertNextNode([0.0, 0.0, 0.0]) + mesh.insertNextNode([1.0, 0.0, 0.0]) + mesh.insertNextNode([0.0, 1.0, 0.0]) + mesh.insertNextTriangle([0, 1, 2]) + self.assertEqual(len(mesh.getNodes()), 3) + self.assertEqual(len(mesh.getTriangles()), 1) + + # MaterialMap + print(" > Testing MaterialMap") + mm = vls.MaterialMap() + mm.insertNextMaterial(0) + self.assertEqual(mm.getNumberOfMaterials(), 1) + + # Vector Data + print(" > Testing PointData vector insertion") + pd.insertNextVectorData([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], "Vectors") + self.assertEqual(pd.getVectorDataSize(), 1) + + print(" > Done test_mesh_and_data") + + def test_algorithms(self): + print("\n[TEST] test_algorithms (2D & 3D)") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + gridDelta = 0.2 + bounds = [-10.0, 10.0] * dim + dom = module.Domain(bounds, [vls.BoundaryConditionEnum.REFLECTIVE_BOUNDARY] * dim, gridDelta) + + # Geometries + origin = [0.0] * dim + radius = 5.0 + print(" >> MakeGeometry") + sphere_geom = module.Sphere(origin, radius) + module.MakeGeometry(dom, sphere_geom).apply() + self.assertGreater(dom.getNumberOfPoints(), 0) + + print(" >> Advect") + # Advect + vel = TestVelocityField() + adv = module.Advect() + adv.insertNextLevelSet(dom) + adv.setVelocityField(vel) + adv.setAdvectionTime(0.1) + adv.setSpatialScheme(vls.SpatialSchemeEnum.ENGQUIST_OSHER_1ST_ORDER) + adv.setTemporalScheme(vls.TemporalSchemeEnum.FORWARD_EULER) + adv.apply() + + print(" >> GeometricAdvect") + # Geometric Advect & Distributions + dist = module.SphereDistribution(radius) + gadv = module.GeometricAdvect(dom, dist) + gadv.apply() + + # Check other distributions instantiation + print(" >> Checking Distribution instantiation") + _ = module.BoxDistribution([1.0] * 3) # BoxDistribution always takes 3D bounds in constructor? No, usually dim. + # Checking pyWrap.hpp: BoxDistribution takes std::array in constructor? + # py::class_ ... .def(py::init(... std::array)) + # It seems BoxDistribution constructor in Python bindings expects 3 args regardless of D? + # Let's check pyWrap.hpp again. + # .def(py::init(&SmartPointer>::template New>)) + # Yes, it seems hardcoded to 3 in the binding signature for BoxDistribution. + _ = module.BoxDistribution([1.0, 1.0, 1.0]) + _ = module.CustomSphereDistribution([1.0, 2.0]) + + print(" >> Boolean") + # Boolean + dom2 = module.Domain(dom) + bool_op = module.BooleanOperation(dom, dom2, vls.BooleanOperationEnum.UNION) + bool_op.apply() + + print(" >> CompareVolume") + # CompareVolume + dom_vol1 = module.Domain(gridDelta) + module.MakeGeometry(dom_vol1, module.Sphere([0.0] * dim, 5.0)).apply() + dom_vol2 = module.Domain(gridDelta) + module.MakeGeometry(dom_vol2, module.Sphere([0.0] * dim, 6.0)).apply() + cv = module.CompareVolume(dom_vol1, dom_vol2) + cv.apply() + self.assertGreater(cv.getVolumeMismatch(), 0.0) + + print(" >> Calc Curvatures") + # Calc Curvatures + curv = module.CalculateCurvatures(dom) + curv.apply() + + print(" >> Calc Normals") + # Calc Normals + norms = module.CalculateNormalVectors(dom) + norms.apply() + + print(" >> Calc Visibilities") + # Calc Visibilities + # CalculateVisibilities expects Vec3D (3 components) even in 2D + module.CalculateVisibilities(dom, [0.0, 0.0, 1.0], "Visibility").apply() + + print(" >> Check") + # Check + module.Check(dom).apply() + + print(" >> Prune") + # Prune + module.Prune(dom).apply() + + print(" >> Expand") + # Expand + module.Expand(dom, 1).apply() + + print(" >> Reduce") + # Reduce + module.Reduce(dom, 1).apply() + + print(" >> MarkVoidPoints") + # Mark Void Points + module.MarkVoidPoints(dom).apply() + + # Expand to ensure sufficient width for feature detection (curvature) + module.Prune(dom).apply() + module.Expand(dom, 5).apply() + + print(" >> DetectFeatures") + # Detect Features + module.DetectFeatures(dom).apply() + + print(" >> Mesh Conversions") + # Mesh Conversions + # Expand back to a safe width for meshing + module.Expand(dom, 3).apply() + + mesh = vls.Mesh() + print(" >>> ToMesh") + module.ToMesh(dom, mesh).apply() + print(" >>> ToSurfaceMesh") + module.ToSurfaceMesh(dom, mesh).apply() + print(" >>> ToDiskMesh") + module.ToDiskMesh(dom, mesh).apply() + + print(" >> IO") + # IO + fname = f"test_{dim}d.vtp" + vls.VTKWriter(mesh, fname).apply() + if os.path.exists(fname): + os.remove(fname) + + print(" >> RemoveStrayPoints") + # Remove Stray Points + module.RemoveStrayPoints(dom).apply() + + print(" > Done test_algorithms") + + def test_geometries(self): + print("\n[TEST] test_geometries") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + # Initialize domain with bounds to ensure MakeGeometry(Plane) works efficiently + bounds = [-10.0, 10.0] * dim + boundaryCons = [vls.BoundaryConditionEnum.REFLECTIVE_BOUNDARY] * dim + # Make one boundary infinite for variety + boundaryCons[-1] = vls.BoundaryConditionEnum.INFINITE_BOUNDARY + dom = module.Domain(bounds, boundaryCons, 0.2) + + # Plane + print(" >> Plane") + normal = [0.0] * dim + normal[-1] = 1.0 + plane = module.Plane([0.0] * dim, normal) + module.MakeGeometry(dom, plane).apply() + self.assertGreater(dom.getNumberOfPoints(), 0) + + # Box + print(" >> Box") + box = module.Box([-1.0] * dim, [1.0] * dim) + module.MakeGeometry(dom, box).apply() + + # Cylinder + print(" >> Cylinder") + # Cylinder constructor takes 3D vectors for origin and axis even in 2D? + # pyWrap.hpp: .def(py::init(&SmartPointer>::template New & ...>)) + # It uses std::vector, so it should adapt to D. + axis = [0.0] * dim + axis[-1] = 1.0 + cyl = module.Cylinder([0.0] * dim, axis, 2.0, 1.0, 1.0) + module.MakeGeometry(dom, cyl).apply() + + print(" > Done test_geometries") + + def test_io_and_transforms(self): + print("\n[TEST] test_io_and_transforms") + # Create a mesh + mesh = vls.Mesh() + mesh.insertNextNode([0.0, 0.0, 0.0]) + mesh.insertNextNode([1.0, 0.0, 0.0]) + mesh.insertNextNode([0.0, 1.0, 0.0]) + mesh.insertNextTriangle([0, 1, 2]) + + # Transform + print(" > TransformMesh") + trans = vls.TransformMesh(mesh, vls.TransformEnum.TRANSLATION, [1.0, 0.0, 0.0], 0.0) + trans.apply() + nodes = mesh.getNodes() + self.assertEqual(nodes[0][0], 1.0) + + # VTK Writer + print(" > VTKWriter") + vls.VTKWriter(mesh, "test_io.vtp").apply() + + # VTK Reader + print(" > VTKReader") + mesh2 = vls.Mesh() + vls.VTKReader(mesh2, vls.FileFormatEnum.VTP, "test_io.vtp").apply() + self.assertEqual(len(mesh2.getNodes()), 3) + + if os.path.exists("test_io.vtp"): + os.remove("test_io.vtp") + print(" > Done test_io_and_transforms") + + def test_mesh_conversion_roundtrip(self): + print("\n[TEST] test_mesh_conversion_roundtrip") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + dom = module.Domain(0.2) + module.MakeGeometry(dom, module.Sphere([0.0]*dim, 5.0)).apply() + + mesh = vls.Mesh() + module.ToSurfaceMesh(dom, mesh).apply() + + dom2 = module.Domain(0.2) + module.FromSurfaceMesh(dom2, mesh).apply() + self.assertGreater(dom2.getNumberOfPoints(), 0) + print(" > Done test_mesh_conversion_roundtrip") + + def test_native_io(self): + print("\n[TEST] test_native_io") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + dom = module.Domain(0.2) + module.MakeGeometry(dom, module.Sphere([0.0]*dim, 5.0)).apply() + + filename = f"test_native_{dim}d.lvst" + print(" >> Writer") + module.Writer(dom, filename).apply() + + dom_in = module.Domain(0.2) + print(" >> Reader") + module.Reader(dom_in, filename).apply() + self.assertGreater(dom_in.getNumberOfPoints(), 0) + + if os.path.exists(filename): + os.remove(filename) + print(" > Done test_native_io") + + def test_advanced_meshing(self): + print("\n[TEST] test_advanced_meshing") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + dom1 = module.Domain(0.2) + module.MakeGeometry(dom1, module.Sphere([0.0]*dim, 5.0)).apply() + + dom2 = module.Domain(0.2) + center = [0.0]*dim + center[0] = 2.0 + module.MakeGeometry(dom2, module.Sphere(center, 5.0)).apply() + + # ToMultiSurfaceMesh + print(" >> ToMultiSurfaceMesh") + mesh_multi = vls.Mesh() + module.ToMultiSurfaceMesh([dom1, dom2], mesh_multi).apply() + self.assertGreater(len(mesh_multi.getNodes()), 0) + + # ToVoxelMesh + print(" >> ToVoxelMesh") + mesh_voxel = vls.Mesh() + module.ToVoxelMesh([dom1, dom2], mesh_voxel).apply() + self.assertGreater(len(mesh_voxel.getNodes()), 0) + print(" > Done test_advanced_meshing") + + def test_stencil_functions(self): + print("\n[TEST] test_stencil_functions") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + dom = module.Domain(0.2) + module.MakeGeometry(dom, module.Sphere([0.0]*dim, 5.0)).apply() + module.Expand(dom, 5).apply() + + # Create a second domain to form a valid multi-layer system + dom2 = module.Domain(0.2) + module.MakeGeometry(dom2, module.Sphere([0.0]*dim, 10.0)).apply() + module.Expand(dom2, 5).apply() + + # PrepareStencilLocalLaxFriedrichs + print(" >> PrepareStencilLocalLaxFriedrichs") + module.PrepareStencilLocalLaxFriedrichs([dom, dom2], [False, True]) + + # Finalize + print(" >> FinalizeStencilLocalLaxFriedrichs") + module.FinalizeStencilLocalLaxFriedrichs([dom, dom2]) + print(" > Done test_stencil_functions") + + def test_geometric_advect_box(self): + print("\n[TEST] test_geometric_advect_box") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + dom = module.Domain(0.2) + module.MakeGeometry(dom, module.Sphere([0.0]*dim, 5.0)).apply() + + print(" >> BoxDistribution") + # BoxDistribution constructor in Python bindings expects 3 args regardless of D + dist = module.BoxDistribution([1.0, 1.0, 1.0]) + + print(" >> GeometricAdvect") + module.GeometricAdvect(dom, dist).apply() + self.assertGreater(dom.getNumberOfPoints(), 0) + print(" > Done test_geometric_advect_box") + + def test_from_mesh_volume(self): + print("\n[TEST] test_from_mesh_volume") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + mesh = vls.Mesh() + if dim == 3: + # Create a simple tetrahedral mesh + mesh.insertNextNode([0.0, 0.0, 0.0]) + mesh.insertNextNode([5.0, 0.0, 0.0]) + mesh.insertNextNode([0.0, 5.0, 0.0]) + mesh.insertNextNode([0.0, 0.0, 5.0]) + mesh.insertNextTetra([0, 1, 2, 3]) + ls_values = [-1.0, 1.0, 1.0, 1.0] + else: + # Create a simple triangle mesh for 2D volume + mesh.insertNextNode([0.0, 0.0, 0.0]) + mesh.insertNextNode([5.0, 0.0, 0.0]) + mesh.insertNextNode([0.0, 5.0, 0.0]) + mesh.insertNextTriangle([0, 1, 2]) + ls_values = [-1.0, 1.0, 1.0] + + # FromMesh requires "LSValues" in cellData corresponding to nodes + mesh.getPointData().insertNextScalarData(ls_values, "LSValues") + + dom = module.Domain(0.2) + # FromMesh should handle volume meshes if implemented in the wrapper logic + print(" >> FromMesh (Volume)") + module.FromMesh(dom, mesh).apply() + self.assertGreater(dom.getNumberOfPoints(), 0) + print(" > Done test_from_mesh_volume") + + def test_boolean_variants(self): + print("\n[TEST] test_boolean_variants") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + dom1 = module.Domain(0.2) + module.MakeGeometry(dom1, module.Sphere([0.0]*dim, 5.0)).apply() + + dom2 = module.Domain(0.2) + center = [0.0]*dim + center[0] = 2.0 + module.MakeGeometry(dom2, module.Sphere(center, 5.0)).apply() + + # Intersect + print(" >> Intersect") + dom_int = module.Domain(dom1) + module.BooleanOperation(dom_int, dom2, vls.BooleanOperationEnum.INTERSECT).apply() + self.assertGreater(dom_int.getNumberOfPoints(), 0) + + # Relative Complement (Difference) + print(" >> Relative Complement") + dom_diff = module.Domain(dom1) + module.BooleanOperation(dom_diff, dom2, vls.BooleanOperationEnum.RELATIVE_COMPLEMENT).apply() + self.assertGreater(dom_diff.getNumberOfPoints(), 0) + + # Invert + print(" >> Invert") + dom_inv = module.Domain(dom1) + module.BooleanOperation(dom_inv, vls.BooleanOperationEnum.INVERT).apply() + # Inverted domain is infinite, but points should exist + self.assertGreater(dom_inv.getNumberOfPoints(), 0) + print(" > Done test_boolean_variants") + + def test_domain_manipulation(self): + print("\n[TEST] test_domain_manipulation") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + dom = module.Domain(0.2) + module.MakeGeometry(dom, module.Sphere([0.0]*dim, 5.0)).apply() + + # Deep Copy + print(" >> deepCopy") + dom_copy = module.Domain(0.2) + dom_copy.deepCopy(dom) + self.assertEqual(dom.getNumberOfPoints(), dom_copy.getNumberOfPoints()) + + # Level Set Width + print(" >> setLevelSetWidth") + dom.setLevelSetWidth(3) + self.assertEqual(dom.getLevelSetWidth(), 3) + print(" > Done test_domain_manipulation") + + def test_vtk_metadata(self): + print("\n[TEST] test_vtk_metadata") + mesh = vls.Mesh() + mesh.insertNextNode([0.0, 0.0, 0.0]) + + writer = vls.VTKWriter(mesh, "test_meta.vtp") + writer.addMetaData("Value", 42.0) + writer.apply() + + if os.path.exists("test_meta.vtp"): + os.remove("test_meta.vtp") + print(" > Done test_vtk_metadata") + + def test_write_visualization_mesh(self): + print("\n[TEST] test_write_visualization_mesh") + if hasattr(d3, "WriteVisualizationMesh"): + dom = d3.Domain(0.2) + d3.MakeGeometry(dom, d3.Sphere([0.0, 0.0, 0.0], 5.0)).apply() + + wvm = d3.WriteVisualizationMesh(dom) + wvm.setFileName("test_vis") + wvm.apply() + else: + print(" > WriteVisualizationMesh not available (VIENNALS_USE_VTK off?)") + print(" > Done test_write_visualization_mesh") + + def test_comparisons(self): + print("\n[TEST] test_comparisons") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + gridDelta = 0.5 + dom = module.Domain(gridDelta) + + # Make Geometry + module.MakeGeometry(dom, module.Sphere([0.0]*dim, 5.0)).apply() + + dom2 = module.Domain(dom) + # Shift dom2 slightly + center = [0.0]*dim + center[0] = 0.1 + module.MakeGeometry(dom2, module.Sphere(center, 5.0)).apply() + + print(" >> Comparisons") + # Comparisons + module.CompareVolume(dom, dom2).apply() + module.CompareChamfer(dom, dom2).apply() + + comp_crit = module.CompareCriticalDimensions(dom, dom2) + comp_crit.addXRange(-1.0, 1.0, True) + comp_crit.apply() + + module.CompareNarrowBand(dom, dom2).apply() + + print(" >> Sparse Field") + # Sparse Field (needs expanded/reduced) + dom_exp = module.Domain(dom) + module.Expand(dom_exp, 55).apply() + module.CompareSparseField(dom_exp, dom).apply() + print(" > Done test_comparisons") + + def test_cross_dimension(self): + print("\n[TEST] test_cross_dimension") + # Slice 3D -> 2D + print(" > Creating 3D domain") + dom3 = d3.Domain(0.2) + d3.MakeGeometry(dom3, d3.Sphere([0.0, 0.0, 0.0], 5.0)).apply() + + dom2 = d2.Domain(0.2) + + # Slice at z=0 + print(" > Slice") + s = vls.Slice(dom3, dom2, 2, 0.0) + s.apply() + self.assertGreater(dom2.getNumberOfPoints(), 0) + + # Extrude 2D -> 3D + print(" > Extrude") + dom3_out = d3.Domain(0.2) + extent = [-1.0, 1.0] + bcs = [vls.BoundaryConditionEnum.INFINITE_BOUNDARY] * 3 + + e = vls.Extrude(dom2, dom3_out, extent, 2, bcs) + e.apply() + self.assertGreater(dom3_out.getNumberOfPoints(), 0) + print(" > Done test_cross_dimension") + + def test_point_cloud_convex_hull(self): + print("\n[TEST] test_point_cloud_convex_hull") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + # PointCloud + print(" >> Creating PointCloud") + # Create points on axes + points = [[0.0]*dim] + [[1.0 if i == j else 0.0 for i in range(dim)] for j in range(dim)] + pc = module.PointCloud(points) + + # ConvexHull + print(" >> ConvexHull") + mesh = vls.Mesh() + hull = module.ConvexHull(mesh, pc) + hull.apply() + self.assertGreater(len(mesh.getNodes()), 0) + print(" > Done test_point_cloud_convex_hull") + + def test_advection_callback(self): + print("\n[TEST] test_advection_callback") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + print(" >> Setting up advection") + dom = module.Domain(0.2) + module.MakeGeometry(dom, module.Sphere([0.0]*dim, 5.0)).apply() + + vel = TestVelocityField() + adv = module.Advect() + adv.insertNextLevelSet(dom) + adv.setVelocityField(vel) + adv.setAdvectionTime(0.1) + adv.setTemporalScheme(vls.TemporalSchemeEnum.RUNGE_KUTTA_2ND_ORDER) + + # Use a mutable container to capture callback execution + callback_data = {'called': False} + def callback(d): + callback_data['called'] = True + return True + + print(" >> Running advection with callback") + adv.setVelocityUpdateCallback(callback) + adv.apply() + + self.assertTrue(callback_data['called']) + print(" > Done test_advection_callback") + + def test_make_geometry_boundary(self): + print("\n[TEST] test_make_geometry_boundary") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + print(" >> Creating domain with bounds") + bounds = [-5.0, 5.0] * dim + boundaryCons = [vls.BoundaryConditionEnum.REFLECTIVE_BOUNDARY] * dim + dom = module.Domain(bounds, boundaryCons, 1.0) + + print(" >> MakeGeometry with IgnoreBoundaryConditions") + maker = module.MakeGeometry(dom, module.Sphere([0.0]*dim, 4.0)) + maker.setIgnoreBoundaryConditions(True) + maker.apply() + + self.assertGreater(dom.getNumberOfPoints(), 0) + print(" > Done test_make_geometry_boundary") + + def test_mesh_operations(self): + print("\n[TEST] test_mesh_operations") + print(" > Inserting duplicate nodes") + mesh = vls.Mesh() + # Insert duplicate nodes + mesh.insertNextNode([0.0, 0.0, 0.0]) + mesh.insertNextNode([0.0, 0.0, 0.0]) + self.assertEqual(len(mesh.getNodes()), 2) + + print(" > removeDuplicateNodes") + mesh.removeDuplicateNodes() + self.assertEqual(len(mesh.getNodes()), 1) + + # Append + print(" > append") + mesh2 = vls.Mesh() + mesh2.insertNextNode([1.0, 1.0, 1.0]) + + mesh.append(mesh2) + self.assertEqual(len(mesh.getNodes()), 2) + print(" > Done test_mesh_operations") + + def test_advect_configuration(self): + print("\n[TEST] test_advect_configuration") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + print(" >> Creating domain") + dom = module.Domain(0.2) + module.MakeGeometry(dom, module.Sphere([0.0]*dim, 5.0)).apply() + + vel = TestVelocityField() + adv = module.Advect() + adv.insertNextLevelSet(dom) + adv.setVelocityField(vel) + + # Test Single Step + print(" >> Testing Single Step") + adv.setAdvectionTime(10.0) # Large time that would normally require multiple steps + adv.setSingleStep(True) + adv.apply() + self.assertEqual(adv.getNumberOfTimeSteps(), 1) + print(" > Done test_advect_configuration") + + def test_mark_void_points_advanced(self): + print("\n[TEST] test_mark_void_points_advanced") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + print(" >> Creating disjoint spheres") + dom = module.Domain(0.2) + # Create two disjoint spheres + module.MakeGeometry(dom, module.Sphere([0.0]*dim, 3.0)).apply() + + dom2 = module.Domain(0.2) + center = [0.0]*dim + center[0] = 10.0 + module.MakeGeometry(dom2, module.Sphere(center, 3.0)).apply() + + # Union them to create a domain with two separate components + print(" >> Union") + module.BooleanOperation(dom, dom2, vls.BooleanOperationEnum.UNION).apply() + + print(" >> MarkVoidPoints") + marker = module.MarkVoidPoints(dom) + marker.setSaveComponentIds(True) + marker.apply() + + # Should find 3 connected components (2 spheres + 1 background) + self.assertEqual(marker.getNumberOfComponents(), 3) + print(" > Done test_mark_void_points_advanced") + + def test_adaptive_time_stepping(self): + print("\n[TEST] test_adaptive_time_stepping") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + print(" >> Setting up advection") + dom = module.Domain(0.2) + module.MakeGeometry(dom, module.Sphere([0.0]*dim, 5.0)).apply() + + vel = TestVelocityField() + adv = module.Advect() + adv.insertNextLevelSet(dom) + adv.setVelocityField(vel) + adv.setAdvectionTime(0.1) + + # Enable adaptive time stepping + print(" >> Enabling Adaptive Time Stepping") + adv.setAdaptiveTimeStepping(True, 10) + adv.apply() + + self.assertGreater(dom.getNumberOfPoints(), 0) + print(" > Done test_adaptive_time_stepping") + + def test_geometric_distribution_functions(self): + print("\n[TEST] test_geometric_distribution_functions") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + print(" >> SphereDistribution") + # Sphere Distribution + dist = module.SphereDistribution(5.0) + # Point inside (0,0,0) relative to center (0,0,0) + print(" >> Checking isInside") + # isInside takes 3D points in bindings? + # pyWrap.hpp: .def("isInside", &GeometricAdvectDistribution::isInside) + # GeometricAdvectDistribution uses std::array in pyWrap.hpp for PylsGeometricAdvectDistribution trampoline + # but the base class is templated on D. + # However, the trampoline class PylsGeometricAdvectDistribution hardcodes vectorType to std::array. + # This suggests that even for 2D, the distribution functions might expect 3D coordinates in the Python binding layer if they go through the trampoline. + # But SphereDistribution is a direct binding. + # Let's assume it takes 3D points because of the trampoline definition in pyWrap.hpp which seems to force 3D for the virtual methods. + # Wait, SphereDistribution binding doesn't use the trampoline unless it inherits from it in Python, which it doesn't here. + # But the C++ class SphereDistribution uses VectorType. + # Let's try passing 3D points to be safe, as extra dimensions are usually ignored in 2D logic if passed as std::array to a function expecting std::array? No, pybind11 would complain. + # Actually, looking at pyWrap.hpp, PylsGeometricAdvectDistribution is defined with `typedef std::array vectorType;` + # This looks like a bug/limitation in the bindings for 2D if it enforces 3D. + # However, let's try using `[0.0]*3` which is safe for 3D and might work for 2D if the binding expects 3 args. + # If the binding expects 2 args for 2D, `[0.0]*3` will fail. + # Let's check `GeometricAdvectDistribution` binding in `pyWrap.hpp`. + # It uses `PylsGeometricAdvectDistribution` as trampoline. + # `PylsGeometricAdvectDistribution` has `typedef std::array vectorType;` regardless of D. + # This strongly suggests that the Python side expects 3-element lists/tuples even for 2D distributions. + + pt_zero = [0.0] * 3 + pt_out = [0.0] * 3 + pt_out[0] = 6.0 + + self.assertTrue(dist.isInside(pt_zero, pt_zero, 1e-6)) + self.assertFalse(dist.isInside(pt_zero, pt_out, 1e-6)) + + # Signed Distance + print(" >> Checking getSignedDistance") + dist_val = dist.getSignedDistance(pt_zero, pt_zero, 0) + self.assertAlmostEqual(dist_val, -5.0, delta=1e-5) + print(" > Done test_geometric_distribution_functions") + + def test_advect_additional_options(self): + print("\n[TEST] test_advect_additional_options") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + print(" >> Setting up advection") + dom = module.Domain(0.2) + module.MakeGeometry(dom, module.Sphere([0.0]*dim, 5.0)).apply() + + vel = TestVelocityField() + adv = module.Advect() + adv.insertNextLevelSet(dom) + adv.setVelocityField(vel) + adv.setAdvectionTime(0.1) + + print(" >> Setting IgnoreVoids") + # Test Ignore Voids + adv.setIgnoreVoids(True) + + print(" >> Setting DissipationAlpha") + # Test Dissipation Alpha + adv.setDissipationAlpha(0.5) + + print(" >> Setting CheckDissipation") + # Test Check Dissipation + adv.setCheckDissipation(False) + + adv.apply() + + self.assertGreater(dom.getNumberOfPoints(), 0) + print(" > Done test_advect_additional_options") + + def test_additional_options(self): + print("\n[TEST] test_additional_options") + for module, dim in [(d2, 2), (d3, 3)]: + print(f" > Testing {dim}D") + print(" >> Creating domain with bounds") + bounds = [-10.0, 10.0] * dim + boundaryCons = [vls.BoundaryConditionEnum.REFLECTIVE_BOUNDARY] * dim + dom = module.Domain(bounds, boundaryCons, 0.2) + module.MakeGeometry(dom, module.Sphere([0.0]*dim, 5.0)).apply() + + # ToMesh options + print(" >> ToMesh options (OnlyDefined, OnlyActive)") + mesh = vls.Mesh() + tomesh = module.ToMesh(dom, mesh) + tomesh.setOnlyDefined(True) + tomesh.setOnlyActive(True) + tomesh.apply() + self.assertGreater(len(mesh.getNodes()), 0) + + # DetectFeatures options + print(" >> DetectFeatures options (Threshold, Method)") + df = module.DetectFeatures(dom) + df.setDetectionThreshold(10.0) + df.setDetectionMethod(vls.FeatureDetectionEnum.NORMALS_ANGLE) + df.apply() + + # MarkVoidPoints options + print(" >> MarkVoidPoints options (Reverse, LargestSurface)") + mvp = module.MarkVoidPoints(dom) + mvp.setReverseVoidDetection(True) + mvp.setDetectLargestSurface(True) + mvp.apply() + + # RemoveStrayPoints options + print(" >> RemoveStrayPoints options (VoidTopSurface)") + rsp = module.RemoveStrayPoints(dom) + rsp.setVoidTopSurface(vls.VoidTopSurfaceEnum.LARGEST) + rsp.apply() + + print(" > Done test_additional_options") + +def run_test_method(method_name): + """Run a single test method in a separate process.""" + suite = unittest.TestSuite() + suite.addTest(WrappingTest(method_name)) + # buffer=True captures stdout/stderr to prevent interleaving + runner = unittest.TextTestRunner(verbosity=2, buffer=True) + result = runner.run(suite) + return result.wasSuccessful() + +if __name__ == "__main__": + vls.setNumThreads(8) + unittest.main() diff --git a/python/viennals/__init__.pyi b/python/viennals/__init__.pyi index 79b08550..550f2a77 100644 --- a/python/viennals/__init__.pyi +++ b/python/viennals/__init__.pyi @@ -22,8 +22,8 @@ from viennals._core import MaterialMap from viennals._core import Mesh from viennals._core import PointData from viennals._core import Slice -from viennals._core import SpatialSchemeEnum as IntegrationSchemeEnum from viennals._core import SpatialSchemeEnum +from viennals._core import SpatialSchemeEnum as IntegrationSchemeEnum from viennals._core import TemporalSchemeEnum from viennals._core import TransformEnum from viennals._core import TransformMesh @@ -41,11 +41,12 @@ from viennals.d2 import CalculateCurvatures from viennals.d2 import CalculateNormalVectors from viennals.d2 import CalculateVisibilities from viennals.d2 import Check -from viennals.d2 import CompareArea from viennals.d2 import CompareChamfer from viennals.d2 import CompareCriticalDimensions from viennals.d2 import CompareNarrowBand from viennals.d2 import CompareSparseField +from viennals.d2 import CompareVolume as CompareArea +from viennals.d2 import CompareVolume from viennals.d2 import ConvexHull from viennals.d2 import CustomSphereDistribution from viennals.d2 import Cylinder @@ -81,7 +82,7 @@ from viennals.d2 import hrleGrid from . import _core from . import d2 from . import d3 -__all__: list[str] = ['Advect', 'BooleanOperation', 'BooleanOperationEnum', 'BoundaryConditionEnum', 'Box', 'BoxDistribution', 'CalculateCurvatures', 'CalculateNormalVectors', 'CalculateVisibilities', 'Check', 'CompareArea', 'CompareChamfer', 'CompareCriticalDimensions', 'CompareNarrowBand', 'CompareSparseField', 'ConvexHull', 'CurvatureEnum', 'CustomSphereDistribution', 'Cylinder', 'DetectFeatures', 'Domain', 'Expand', 'Extrude', 'FeatureDetectionEnum', 'FileFormatEnum', 'FinalizeStencilLocalLaxFriedrichs', 'FromMesh', 'FromSurfaceMesh', 'FromVolumeMesh', 'GeometricAdvect', 'GeometricAdvectDistribution', 'IntegrationSchemeEnum', 'LogLevel', 'Logger', 'MakeGeometry', 'MarkVoidPoints', 'MaterialMap', 'Mesh', 'PROXY_DIM', 'Plane', 'PointCloud', 'PointData', 'PrepareStencilLocalLaxFriedrichs', 'Prune', 'Reader', 'Reduce', 'RemoveStrayPoints', 'Slice', 'SpatialSchemeEnum', 'Sphere', 'SphereDistribution', 'StencilLocalLaxFriedrichsScalar', 'TemporalSchemeEnum', 'ToDiskMesh', 'ToMesh', 'ToMultiSurfaceMesh', 'ToSurfaceMesh', 'ToVoxelMesh', 'TransformEnum', 'TransformMesh', 'VTKReader', 'VTKRenderWindow', 'VTKWriter', 'VelocityField', 'VoidTopSurfaceEnum', 'WriteVisualizationMesh', 'Writer', 'd2', 'd3', 'getDimension', 'hrleGrid', 'setDimension', 'setNumThreads', 'version'] +__all__: list[str] = ['Advect', 'BooleanOperation', 'BooleanOperationEnum', 'BoundaryConditionEnum', 'Box', 'BoxDistribution', 'CalculateCurvatures', 'CalculateNormalVectors', 'CalculateVisibilities', 'Check', 'CompareArea', 'CompareChamfer', 'CompareCriticalDimensions', 'CompareNarrowBand', 'CompareSparseField', 'CompareVolume', 'ConvexHull', 'CurvatureEnum', 'CustomSphereDistribution', 'Cylinder', 'DetectFeatures', 'Domain', 'Expand', 'Extrude', 'FeatureDetectionEnum', 'FileFormatEnum', 'FinalizeStencilLocalLaxFriedrichs', 'FromMesh', 'FromSurfaceMesh', 'FromVolumeMesh', 'GeometricAdvect', 'GeometricAdvectDistribution', 'IntegrationSchemeEnum', 'LogLevel', 'Logger', 'MakeGeometry', 'MarkVoidPoints', 'MaterialMap', 'Mesh', 'PROXY_DIM', 'Plane', 'PointCloud', 'PointData', 'PrepareStencilLocalLaxFriedrichs', 'Prune', 'Reader', 'Reduce', 'RemoveStrayPoints', 'Slice', 'SpatialSchemeEnum', 'Sphere', 'SphereDistribution', 'StencilLocalLaxFriedrichsScalar', 'TemporalSchemeEnum', 'ToDiskMesh', 'ToMesh', 'ToMultiSurfaceMesh', 'ToSurfaceMesh', 'ToVoxelMesh', 'TransformEnum', 'TransformMesh', 'VTKReader', 'VTKRenderWindow', 'VTKWriter', 'VelocityField', 'VoidTopSurfaceEnum', 'WriteVisualizationMesh', 'Writer', 'd2', 'd3', 'getDimension', 'hrleGrid', 'setDimension', 'setNumThreads', 'version'] def __dir__(): ... def __getattr__(name): diff --git a/python/viennals/_core.pyi b/python/viennals/_core.pyi index adb22c02..3adc3a1c 100644 --- a/python/viennals/_core.pyi +++ b/python/viennals/_core.pyi @@ -7,8 +7,8 @@ import enum import typing from viennals import d2 import viennals.d2 -import viennals.d3 from viennals import d3 +import viennals.d3 __all__: list[str] = ['BooleanOperationEnum', 'BoundaryConditionEnum', 'CurvatureEnum', 'Extrude', 'FeatureDetectionEnum', 'FileFormatEnum', 'IntegrationSchemeEnum', 'LogLevel', 'Logger', 'MaterialMap', 'Mesh', 'PointData', 'Slice', 'SpatialSchemeEnum', 'TemporalSchemeEnum', 'TransformEnum', 'TransformMesh', 'VTKReader', 'VTKRenderWindow', 'VTKWriter', 'VelocityField', 'VoidTopSurfaceEnum', 'd2', 'd3', 'setNumThreads', 'version'] class BooleanOperationEnum(enum.IntEnum): INTERSECT: typing.ClassVar[BooleanOperationEnum] # value = diff --git a/python/viennals/d2.pyi b/python/viennals/d2.pyi index d79668d9..fcfd1336 100644 --- a/python/viennals/d2.pyi +++ b/python/viennals/d2.pyi @@ -5,7 +5,7 @@ from __future__ import annotations import collections.abc import typing import viennals._core -__all__: list[str] = ['Advect', 'BooleanOperation', 'Box', 'BoxDistribution', 'CalculateCurvatures', 'CalculateNormalVectors', 'CalculateVisibilities', 'Check', 'CompareArea', 'CompareChamfer', 'CompareCriticalDimensions', 'CompareNarrowBand', 'CompareSparseField', 'ConvexHull', 'CustomSphereDistribution', 'Cylinder', 'DetectFeatures', 'Domain', 'Expand', 'FinalizeStencilLocalLaxFriedrichs', 'FromMesh', 'FromSurfaceMesh', 'FromVolumeMesh', 'GeometricAdvect', 'GeometricAdvectDistribution', 'MakeGeometry', 'MarkVoidPoints', 'Plane', 'PointCloud', 'PrepareStencilLocalLaxFriedrichs', 'Prune', 'Reader', 'Reduce', 'RemoveStrayPoints', 'Sphere', 'SphereDistribution', 'StencilLocalLaxFriedrichsScalar', 'ToDiskMesh', 'ToMesh', 'ToMultiSurfaceMesh', 'ToSurfaceMesh', 'ToVoxelMesh', 'WriteVisualizationMesh', 'Writer', 'hrleGrid'] +__all__: list[str] = ['Advect', 'BooleanOperation', 'Box', 'BoxDistribution', 'CalculateCurvatures', 'CalculateNormalVectors', 'CalculateVisibilities', 'Check', 'CompareArea', 'CompareChamfer', 'CompareCriticalDimensions', 'CompareNarrowBand', 'CompareSparseField', 'CompareVolume', 'ConvexHull', 'CustomSphereDistribution', 'Cylinder', 'DetectFeatures', 'Domain', 'Expand', 'FinalizeStencilLocalLaxFriedrichs', 'FromMesh', 'FromSurfaceMesh', 'FromVolumeMesh', 'GeometricAdvect', 'GeometricAdvectDistribution', 'MakeGeometry', 'MarkVoidPoints', 'Plane', 'PointCloud', 'PrepareStencilLocalLaxFriedrichs', 'Prune', 'Reader', 'Reduce', 'RemoveStrayPoints', 'Sphere', 'SphereDistribution', 'StencilLocalLaxFriedrichsScalar', 'ToDiskMesh', 'ToMesh', 'ToMultiSurfaceMesh', 'ToSurfaceMesh', 'ToVoxelMesh', 'WriteVisualizationMesh', 'Writer', 'hrleGrid'] class Advect: @typing.overload def __init__(self) -> None: @@ -114,6 +114,10 @@ class Advect: """ Set the velocity to use for advection. """ + def setVelocityUpdateCallback(self, arg0: collections.abc.Callable[[Domain], bool]) -> None: + """ + Set a callback function that is called after the level set has been updated during intermediate time integration steps (e.g. RK2, RK3). + """ class BooleanOperation: @typing.overload def __init__(self) -> None: @@ -225,57 +229,6 @@ class Check: """ Set levelset for which to calculate normal vectors. """ -class CompareArea: - @typing.overload - def __init__(self) -> None: - ... - @typing.overload - def __init__(self, arg0: Domain, arg1: Domain) -> None: - ... - def apply(self) -> None: - """ - Computes the area difference between the two level sets. - """ - def getAreaMismatch(self) -> float: - """ - Returns the computed area mismatch. - """ - def getCellCount(self) -> int: - """ - Returns the number of cells where the level sets differ. - """ - def getCustomAreaMismatch(self) -> float: - """ - Returns the computed area mismatch, with custom increments applied. - """ - def getCustomCellCount(self) -> int: - """ - Returns the number of cells where the level sets differ, with custom increments applied. - """ - def setDefaultIncrement(self, arg0: typing.SupportsInt) -> None: - """ - Set default increment value - """ - def setLevelSetSample(self, arg0: Domain) -> None: - """ - Sets the sample level set. - """ - def setLevelSetTarget(self, arg0: Domain) -> None: - """ - Sets the target level set. - """ - def setOutputMesh(self, arg0: viennals._core.Mesh) -> None: - """ - Set the output mesh where difference areas will be stored - """ - def setXRangeAndIncrement(self, arg0: typing.SupportsInt, arg1: typing.SupportsInt, arg2: typing.SupportsInt) -> None: - """ - Sets the x-range and custom increment value - """ - def setYRangeAndIncrement(self, arg0: typing.SupportsInt, arg1: typing.SupportsInt, arg2: typing.SupportsInt) -> None: - """ - Sets the y-range and custom increment value - """ class CompareChamfer: @typing.overload def __init__(self) -> None: @@ -338,6 +291,10 @@ class CompareCriticalDimensions: @typing.overload def __init__(self, arg0: Domain, arg1: Domain) -> None: ... + def addRange(self, measureDimension: typing.SupportsInt, minBounds: typing.Annotated[collections.abc.Sequence[typing.SupportsFloat], "FixedSize(2)"], maxBounds: typing.Annotated[collections.abc.Sequence[typing.SupportsFloat], "FixedSize(2)"], findMaximum: bool = True) -> None: + """ + Add a generic range specification. + """ def addXRange(self, minX: typing.SupportsFloat, maxX: typing.SupportsFloat, findMaximum: bool = True) -> None: """ Add an X range to find maximum or minimum Y position. @@ -409,6 +366,10 @@ class CompareNarrowBand: """ Clear the y-range restriction """ + def clearZRange(self) -> None: + """ + Clear the z-range restriction + """ def getNumPoints(self) -> int: """ Return the number of points used in the comparison. @@ -449,6 +410,10 @@ class CompareNarrowBand: """ Set the y-coordinate range to restrict the comparison area """ + def setZRange(self, arg0: typing.SupportsFloat, arg1: typing.SupportsFloat) -> None: + """ + Set the z-coordinate range to restrict the comparison area + """ class CompareSparseField: @typing.overload def __init__(self) -> None: @@ -468,6 +433,10 @@ class CompareSparseField: """ Clear the y-range restriction """ + def clearZRange(self) -> None: + """ + Clear the z-range restriction + """ def getNumPoints(self) -> int: """ Return the number of points used in the comparison. @@ -516,6 +485,73 @@ class CompareSparseField: """ Set the y-coordinate range to restrict the comparison area """ + def setZRange(self, arg0: typing.SupportsFloat, arg1: typing.SupportsFloat) -> None: + """ + Set the z-coordinate range to restrict the comparison area + """ +class CompareVolume: + @typing.overload + def __init__(self) -> None: + ... + @typing.overload + def __init__(self, arg0: Domain, arg1: Domain) -> None: + ... + def apply(self) -> None: + """ + Computes the volume difference between the two level sets. + """ + def getAreaMismatch(self) -> float: + """ + Returns the computed area mismatch. + """ + def getCellCount(self) -> int: + """ + Returns the number of cells where the level sets differ. + """ + def getCustomAreaMismatch(self) -> float: + """ + Returns the computed area mismatch, with custom increments applied. + """ + def getCustomCellCount(self) -> int: + """ + Returns the number of cells where the level sets differ, with custom increments applied. + """ + def getCustomVolumeMismatch(self) -> float: + """ + Returns the computed volume mismatch, with custom increments applied. + """ + def getVolumeMismatch(self) -> float: + """ + Returns the computed volume mismatch. + """ + def setDefaultIncrement(self, arg0: typing.SupportsInt) -> None: + """ + Set default increment value + """ + def setLevelSetSample(self, arg0: Domain) -> None: + """ + Sets the sample level set. + """ + def setLevelSetTarget(self, arg0: Domain) -> None: + """ + Sets the target level set. + """ + def setOutputMesh(self, arg0: viennals._core.Mesh) -> None: + """ + Set the output mesh where difference areas will be stored + """ + def setXRangeAndIncrement(self, arg0: typing.SupportsInt, arg1: typing.SupportsInt, arg2: typing.SupportsInt) -> None: + """ + Sets the x-range and custom increment value + """ + def setYRangeAndIncrement(self, arg0: typing.SupportsInt, arg1: typing.SupportsInt, arg2: typing.SupportsInt) -> None: + """ + Sets the y-range and custom increment value + """ + def setZRangeAndIncrement(self, arg0: typing.SupportsInt, arg1: typing.SupportsInt, arg2: typing.SupportsInt) -> None: + """ + Sets the z-range and custom increment value + """ class ConvexHull: @typing.overload def __init__(self) -> None: @@ -843,7 +879,7 @@ class MarkVoidPoints: """ Reverse the logic of detecting the top surface. """ - def setSaveComponentsId(self, arg0: bool) -> None: + def setSaveComponentIds(self, arg0: bool) -> None: """ Save the connectivity information of all LS points in the pointData of the level set. """ @@ -1189,3 +1225,4 @@ def FinalizeStencilLocalLaxFriedrichs(levelSets: collections.abc.Sequence[Domain ... def PrepareStencilLocalLaxFriedrichs(levelSets: collections.abc.Sequence[Domain], isDepo: collections.abc.Sequence[bool]) -> None: ... +CompareArea = CompareVolume diff --git a/python/viennals/d3.pyi b/python/viennals/d3.pyi index 495756b0..3ad263f3 100644 --- a/python/viennals/d3.pyi +++ b/python/viennals/d3.pyi @@ -5,7 +5,7 @@ from __future__ import annotations import collections.abc import typing import viennals._core -__all__: list[str] = ['Advect', 'BooleanOperation', 'Box', 'BoxDistribution', 'CalculateCurvatures', 'CalculateNormalVectors', 'CalculateVisibilities', 'Check', 'ConvexHull', 'CustomSphereDistribution', 'Cylinder', 'DetectFeatures', 'Domain', 'Expand', 'FinalizeStencilLocalLaxFriedrichs', 'FromMesh', 'FromSurfaceMesh', 'FromVolumeMesh', 'GeometricAdvect', 'GeometricAdvectDistribution', 'MakeGeometry', 'MarkVoidPoints', 'Plane', 'PointCloud', 'PrepareStencilLocalLaxFriedrichs', 'Prune', 'Reader', 'Reduce', 'RemoveStrayPoints', 'Sphere', 'SphereDistribution', 'StencilLocalLaxFriedrichsScalar', 'ToDiskMesh', 'ToMesh', 'ToMultiSurfaceMesh', 'ToSurfaceMesh', 'ToVoxelMesh', 'WriteVisualizationMesh', 'Writer', 'hrleGrid'] +__all__: list[str] = ['Advect', 'BooleanOperation', 'Box', 'BoxDistribution', 'CalculateCurvatures', 'CalculateNormalVectors', 'CalculateVisibilities', 'Check', 'CompareChamfer', 'CompareCriticalDimensions', 'CompareNarrowBand', 'CompareSparseField', 'CompareVolume', 'ConvexHull', 'CustomSphereDistribution', 'Cylinder', 'DetectFeatures', 'Domain', 'Expand', 'FinalizeStencilLocalLaxFriedrichs', 'FromMesh', 'FromSurfaceMesh', 'FromVolumeMesh', 'GeometricAdvect', 'GeometricAdvectDistribution', 'MakeGeometry', 'MarkVoidPoints', 'Plane', 'PointCloud', 'PrepareStencilLocalLaxFriedrichs', 'Prune', 'Reader', 'Reduce', 'RemoveStrayPoints', 'Sphere', 'SphereDistribution', 'StencilLocalLaxFriedrichsScalar', 'ToDiskMesh', 'ToMesh', 'ToMultiSurfaceMesh', 'ToSurfaceMesh', 'ToVoxelMesh', 'WriteVisualizationMesh', 'Writer', 'hrleGrid'] class Advect: @typing.overload def __init__(self) -> None: @@ -114,6 +114,10 @@ class Advect: """ Set the velocity to use for advection. """ + def setVelocityUpdateCallback(self, arg0: collections.abc.Callable[[Domain], bool]) -> None: + """ + Set a callback function that is called after the level set has been updated during intermediate time integration steps (e.g. RK2, RK3). + """ class BooleanOperation: @typing.overload def __init__(self) -> None: @@ -225,6 +229,329 @@ class Check: """ Set levelset for which to calculate normal vectors. """ +class CompareChamfer: + @typing.overload + def __init__(self) -> None: + ... + @typing.overload + def __init__(self, arg0: Domain, arg1: Domain) -> None: + ... + def apply(self) -> None: + """ + Apply the Chamfer distance calculation. + """ + def getBackwardDistance(self) -> float: + """ + Get the backward distance (average distance from sample to target). + """ + def getChamferDistance(self) -> float: + """ + Get the Chamfer distance (average of forward and backward). + """ + def getForwardDistance(self) -> float: + """ + Get the forward distance (average distance from target to sample). + """ + def getMaxDistance(self) -> float: + """ + Get the maximum nearest-neighbor distance. + """ + def getNumSamplePoints(self) -> int: + """ + Get the number of sample surface points. + """ + def getNumTargetPoints(self) -> int: + """ + Get the number of target surface points. + """ + def getRMSChamferDistance(self) -> float: + """ + Get the RMS Chamfer distance. + """ + def setLevelSetSample(self, arg0: Domain) -> None: + """ + Set the sample level set. + """ + def setLevelSetTarget(self, arg0: Domain) -> None: + """ + Set the target level set. + """ + def setOutputMeshSample(self, arg0: viennals._core.Mesh) -> None: + """ + Set output mesh for sample surface points with distance data. + """ + def setOutputMeshTarget(self, arg0: viennals._core.Mesh) -> None: + """ + Set output mesh for target surface points with distance data. + """ +class CompareCriticalDimensions: + @typing.overload + def __init__(self) -> None: + ... + @typing.overload + def __init__(self, arg0: Domain, arg1: Domain) -> None: + ... + def addRange(self, measureDimension: typing.SupportsInt, minBounds: typing.Annotated[collections.abc.Sequence[typing.SupportsFloat], "FixedSize(3)"], maxBounds: typing.Annotated[collections.abc.Sequence[typing.SupportsFloat], "FixedSize(3)"], findMaximum: bool = True) -> None: + """ + Add a generic range specification. + """ + def addXRange(self, minX: typing.SupportsFloat, maxX: typing.SupportsFloat, findMaximum: bool = True) -> None: + """ + Add an X range to find maximum or minimum Y position. + """ + def addYRange(self, minY: typing.SupportsFloat, maxY: typing.SupportsFloat, findMaximum: bool = True) -> None: + """ + Add a Y range to find maximum or minimum X position. + """ + def apply(self) -> None: + """ + Apply the comparison. + """ + def clearRanges(self) -> None: + """ + Clear all range specifications. + """ + def getAllDifferences(self) -> list[float]: + """ + Get all valid differences as a list. + """ + def getCriticalDimensionResult(self, index: typing.SupportsInt) -> tuple: + """ + Get a specific critical dimension result. Returns (valid, positionTarget, positionSample, difference). + """ + def getMaxDifference(self) -> float: + """ + Get maximum difference across all valid critical dimensions. + """ + def getMeanDifference(self) -> float: + """ + Get mean absolute difference across all valid critical dimensions. + """ + def getNumCriticalDimensions(self) -> int: + """ + Get the number of critical dimensions compared. + """ + def getRMSE(self) -> float: + """ + Get RMSE across all valid critical dimensions. + """ + def setLevelSetSample(self, arg0: Domain) -> None: + """ + Sets the sample level set. + """ + def setLevelSetTarget(self, arg0: Domain) -> None: + """ + Sets the target level set. + """ + def setOutputMesh(self, arg0: viennals._core.Mesh) -> None: + """ + Set the output mesh where critical dimension locations will be stored. + """ +class CompareNarrowBand: + @typing.overload + def __init__(self) -> None: + ... + @typing.overload + def __init__(self, arg0: Domain, arg1: Domain) -> None: + ... + def apply(self) -> None: + """ + Apply the comparison and calculate the sum of squared differences. + """ + def clearXRange(self) -> None: + """ + Clear the x-range restriction + """ + def clearYRange(self) -> None: + """ + Clear the y-range restriction + """ + def clearZRange(self) -> None: + """ + Clear the z-range restriction + """ + def getNumPoints(self) -> int: + """ + Return the number of points used in the comparison. + """ + def getRMSE(self) -> float: + """ + Calculate the root mean square error from previously computed values. + """ + def getSumDifferences(self) -> float: + """ + Return the sum of absolute differences calculated by apply(). + """ + def getSumSquaredDifferences(self) -> float: + """ + Return the sum of squared differences calculated by apply(). + """ + def setLevelSetSample(self, arg0: Domain) -> None: + """ + Sets the sample level set. + """ + def setLevelSetTarget(self, arg0: Domain) -> None: + """ + Sets the target level set. + """ + def setOutputMesh(self, arg0: viennals._core.Mesh, arg1: bool) -> None: + """ + Set the output mesh where difference values will be stored + """ + def setOutputMeshSquaredDifferences(self, arg0: bool) -> None: + """ + Set whether to output squared differences (true) or absolute differences (false) + """ + def setXRange(self, arg0: typing.SupportsFloat, arg1: typing.SupportsFloat) -> None: + """ + Set the x-coordinate range to restrict the comparison area + """ + def setYRange(self, arg0: typing.SupportsFloat, arg1: typing.SupportsFloat) -> None: + """ + Set the y-coordinate range to restrict the comparison area + """ + def setZRange(self, arg0: typing.SupportsFloat, arg1: typing.SupportsFloat) -> None: + """ + Set the z-coordinate range to restrict the comparison area + """ +class CompareSparseField: + @typing.overload + def __init__(self) -> None: + ... + @typing.overload + def __init__(self, arg0: Domain, arg1: Domain) -> None: + ... + def apply(self) -> None: + """ + Apply the comparison and calculate the sum of squared differences. + """ + def clearXRange(self) -> None: + """ + Clear the x-range restriction + """ + def clearYRange(self) -> None: + """ + Clear the y-range restriction + """ + def clearZRange(self) -> None: + """ + Clear the z-range restriction + """ + def getNumPoints(self) -> int: + """ + Return the number of points used in the comparison. + """ + def getNumSkippedPoints(self) -> int: + """ + Return the number of points skipped during comparison. + """ + def getRMSE(self) -> float: + """ + Calculate the root mean square error from previously computed values. + """ + def getSumDifferences(self) -> float: + """ + Return the sum of absolute differences calculated by apply(). + """ + def getSumSquaredDifferences(self) -> float: + """ + Return the sum of squared differences calculated by apply(). + """ + def setExpandedLevelSetWidth(self, arg0: typing.SupportsInt) -> None: + """ + Set the expansion width for the expanded level set + """ + def setFillIteratedWithDistances(self, arg0: bool) -> None: + """ + Set whether to fill the iterated level set with distance values + """ + def setLevelSetExpanded(self, arg0: Domain) -> None: + """ + Sets the expanded level set for comparison. + """ + def setLevelSetIterated(self, arg0: Domain) -> None: + """ + Sets the iterated level set to compare against the expanded one. + """ + def setOutputMesh(self, arg0: viennals._core.Mesh) -> None: + """ + Set the output mesh where difference values will be stored + """ + def setXRange(self, arg0: typing.SupportsFloat, arg1: typing.SupportsFloat) -> None: + """ + Set the x-coordinate range to restrict the comparison area + """ + def setYRange(self, arg0: typing.SupportsFloat, arg1: typing.SupportsFloat) -> None: + """ + Set the y-coordinate range to restrict the comparison area + """ + def setZRange(self, arg0: typing.SupportsFloat, arg1: typing.SupportsFloat) -> None: + """ + Set the z-coordinate range to restrict the comparison area + """ +class CompareVolume: + @typing.overload + def __init__(self) -> None: + ... + @typing.overload + def __init__(self, arg0: Domain, arg1: Domain) -> None: + ... + def apply(self) -> None: + """ + Computes the volume difference between the two level sets. + """ + def getAreaMismatch(self) -> float: + """ + Returns the computed area mismatch. + """ + def getCellCount(self) -> int: + """ + Returns the number of cells where the level sets differ. + """ + def getCustomAreaMismatch(self) -> float: + """ + Returns the computed area mismatch, with custom increments applied. + """ + def getCustomCellCount(self) -> int: + """ + Returns the number of cells where the level sets differ, with custom increments applied. + """ + def getCustomVolumeMismatch(self) -> float: + """ + Returns the computed volume mismatch, with custom increments applied. + """ + def getVolumeMismatch(self) -> float: + """ + Returns the computed volume mismatch. + """ + def setDefaultIncrement(self, arg0: typing.SupportsInt) -> None: + """ + Set default increment value + """ + def setLevelSetSample(self, arg0: Domain) -> None: + """ + Sets the sample level set. + """ + def setLevelSetTarget(self, arg0: Domain) -> None: + """ + Sets the target level set. + """ + def setOutputMesh(self, arg0: viennals._core.Mesh) -> None: + """ + Set the output mesh where difference areas will be stored + """ + def setXRangeAndIncrement(self, arg0: typing.SupportsInt, arg1: typing.SupportsInt, arg2: typing.SupportsInt) -> None: + """ + Sets the x-range and custom increment value + """ + def setYRangeAndIncrement(self, arg0: typing.SupportsInt, arg1: typing.SupportsInt, arg2: typing.SupportsInt) -> None: + """ + Sets the y-range and custom increment value + """ + def setZRangeAndIncrement(self, arg0: typing.SupportsInt, arg1: typing.SupportsInt, arg2: typing.SupportsInt) -> None: + """ + Sets the z-range and custom increment value + """ class ConvexHull: @typing.overload def __init__(self) -> None: @@ -552,7 +879,7 @@ class MarkVoidPoints: """ Reverse the logic of detecting the top surface. """ - def setSaveComponentsId(self, arg0: bool) -> None: + def setSaveComponentIds(self, arg0: bool) -> None: """ Save the connectivity information of all LS points in the pointData of the level set. """ diff --git a/tests/CompareArea/CompareArea.cpp b/tests/CompareArea/CompareArea.cpp deleted file mode 100644 index 9bb6d8ae..00000000 --- a/tests/CompareArea/CompareArea.cpp +++ /dev/null @@ -1,125 +0,0 @@ -#include - -#include -#include -#include -#include -#include -#include - -/** - Test for lsCompareArea that compares the area difference between two level - sets. This test creates two different spheres and measures their area - difference. \example CompareArea.cpp -*/ - -namespace ls = viennals; - -int main() { - constexpr int D = 2; - - omp_set_num_threads(4); - - double extent = 15; - double gridDelta = 0.5; - - double bounds[2 * D] = {-extent, extent, -extent, extent}; - ls::Domain::BoundaryType boundaryCons[D]; - for (unsigned i = 0; i < D; ++i) - boundaryCons[i] = ls::Domain::BoundaryType::REFLECTIVE_BOUNDARY; - - // Create first sphere (target) - auto sphere1 = ls::SmartPointer>::New( - bounds, boundaryCons, gridDelta); - - double origin1[D] = {0., 0.}; - double radius1 = 5.0; - - ls::MakeGeometry( - sphere1, ls::SmartPointer>::New(origin1, radius1)) - .apply(); - - // Create second sphere (sample) with different radius - auto sphere2 = ls::SmartPointer>::New( - bounds, boundaryCons, gridDelta); - - double origin2[D] = {0., 0.}; - double radius2 = 8.0; - - ls::MakeGeometry( - sphere2, ls::SmartPointer>::New(origin2, radius2)) - .apply(); - - // Export both spheres as VTK files for visualization - { - auto mesh = ls::SmartPointer>::New(); - ls::ToMesh(sphere1, mesh).apply(); - ls::VTKWriter(mesh, "sphere1.vtp").apply(); - } - - { - auto mesh = ls::SmartPointer>::New(); - ls::ToMesh(sphere2, mesh).apply(); - ls::VTKWriter(mesh, "sphere2.vtp").apply(); - } - - // Compare the areas with mesh Output - ls::CompareArea compareArea(sphere1, sphere2); - auto mesh = ls::SmartPointer>::New(); // Create mesh for output - compareArea.setOutputMesh(mesh); - compareArea.apply(); - // save mesh to file - ls::VTKWriter(mesh, "areaDifference.vtu").apply(); - - // Calculate theoretical area difference (for circles in 2D) - // Area of circle = π * r² - double theoreticalArea1 = M_PI * radius1 * radius1; - double theoreticalArea2 = M_PI * radius2 * radius2; - double theoreticalDifference = std::abs(theoreticalArea2 - theoreticalArea1); - - // Get the calculated area difference - double calculatedDifference = compareArea.getAreaMismatch(); - unsigned long int cellCount = compareArea.getCellCount(); - - std::cout << "Sphere 1 radius: " << radius1 << std::endl; - std::cout << "Sphere 2 radius: " << radius2 << std::endl; - std::cout << "Theoretical area difference: " << theoreticalDifference - << std::endl; - std::cout << "Calculated area difference: " << calculatedDifference - << std::endl; - std::cout << "Number of differing cells: " << cellCount << std::endl; - std::cout << "Error: " - << std::abs(calculatedDifference - theoreticalDifference) - << std::endl; - - // Test custom increment and range functionality - std::cout << "\nTesting custom increments and ranges:" << std::endl; - - // Set custom increment for whole domain - compareArea.setDefaultIncrement(2); - compareArea.apply(); - std::cout << "Area difference with default increment of 2: " - << compareArea.getCustomAreaMismatch() << std::endl; - std::cout << "Cell count with default increment of 2: " - << compareArea.getCustomCellCount() << std::endl; - - // Set range-specific increment for x-range - compareArea.setDefaultIncrement(1); - compareArea.setXRangeAndIncrement(-5, 5, 3); - compareArea.apply(); - std::cout << "Area difference with x-range increment of 3: " - << compareArea.getCustomAreaMismatch() << std::endl; - std::cout << "Cell count with x-range increment of 3: " - << compareArea.getCustomCellCount() << std::endl; - - // Set range-specific increment for y-range - compareArea.setDefaultIncrement(1); - compareArea.setYRangeAndIncrement(-5, 5, 4); - compareArea.apply(); - std::cout << "Area difference with y-range increment of 4: " - << compareArea.getCustomAreaMismatch() << std::endl; - std::cout << "Cell count with y-range increment of 4: " - << compareArea.getCustomCellCount() << std::endl; - - return 0; -} diff --git a/tests/CompareChamfer/CompareChamfer.cpp b/tests/CompareChamfer/CompareChamfer.cpp index 82fc6ab9..a8886c9f 100644 --- a/tests/CompareChamfer/CompareChamfer.cpp +++ b/tests/CompareChamfer/CompareChamfer.cpp @@ -1,10 +1,11 @@ #include #include #include +#include -#include #include #include +#include #include #include #include @@ -24,72 +25,74 @@ namespace ls = viennals; -int main() { - constexpr int D = 2; - - omp_set_num_threads(4); - +template void runTest() { + std::cout << "Running " << D << "D Test..." << std::endl; 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 (int i = 0; i < 2 * D; ++i) + bounds[i] = (i % 2 == 0) ? -extent : extent; + + typename ls::Domain::BoundaryType boundaryCons[D]; for (unsigned i = 0; i < D; ++i) boundaryCons[i] = ls::Domain::BoundaryType::REFLECTIVE_BOUNDARY; - // Create first circle (target) - auto circle1 = ls::SmartPointer>::New( + // Create first sphere (target) + auto sphere1 = ls::SmartPointer>::New( bounds, boundaryCons, gridDelta); - double origin1[D] = {0., 0.}; + std::vector origin1(D, 0.0); double radius1 = 5.0; ls::MakeGeometry( - circle1, ls::SmartPointer>::New(origin1, radius1)) + sphere1, ls::SmartPointer>::New(origin1, radius1)) .apply(); - // Create second circle (sample) with different center - auto circle2 = ls::SmartPointer>::New( + // Create second sphere (sample) with different center + auto sphere2 = ls::SmartPointer>::New( bounds, boundaryCons, gridDelta); - double origin2[D] = {2., 1.}; // Shifted center - // double origin2[D] = {0., 0.}; // Same center (for testing) + std::vector origin2(D, 0.0); + origin2[0] = 2.; + origin2[1] = 1.; double radius2 = 5.0; // Same radius ls::MakeGeometry( - circle2, ls::SmartPointer>::New(origin2, radius2)) + sphere2, ls::SmartPointer>::New(origin2, radius2)) .apply(); - // Export both circles as VTK files for visualization - std::cout << "Exporting surface meshes..." << std::endl; + // Export both spheres as VTK files for visualization + std::string suffix = "_" + std::to_string(D) + "D.vtp"; + std::cout << "Exporting surface meshes to *" << suffix << "..." << std::endl; { auto meshSurface = ls::SmartPointer>::New(); - ls::ToSurfaceMesh(circle1, meshSurface).apply(); - ls::VTKWriter(meshSurface, "circle1_surface.vtp").apply(); - std::cout << " Circle 1 surface points: " << meshSurface->nodes.size() + ls::ToSurfaceMesh(sphere1, meshSurface).apply(); + ls::VTKWriter(meshSurface, "sphere1_surface" + suffix).apply(); + std::cout << " Sphere 1 surface points: " << meshSurface->nodes.size() << std::endl; } { auto meshSurface = ls::SmartPointer>::New(); - ls::ToSurfaceMesh(circle2, meshSurface).apply(); - ls::VTKWriter(meshSurface, "circle2_surface.vtp").apply(); - std::cout << " Circle 2 surface points: " << meshSurface->nodes.size() + ls::ToSurfaceMesh(sphere2, meshSurface).apply(); + ls::VTKWriter(meshSurface, "sphere2_surface" + suffix).apply(); + std::cout << " Sphere 2 surface points: " << meshSurface->nodes.size() << std::endl; } // Test 1: Basic Chamfer distance calculation std::cout << "\n=== Test 1: Basic Chamfer Distance ===" << std::endl; - std::cout << "Circle 1 center: (" << origin1[0] << ", " << origin1[1] << ")" + std::cout << "Sphere 1 center: (" << origin1[0] << ", " << origin1[1] << ")" << std::endl; - std::cout << "Circle 2 center: (" << origin2[0] << ", " << origin2[1] << ")" - << std::endl; - std::cout << "Expected geometric shift: " - << std::sqrt((origin2[0] - origin1[0]) * (origin2[0] - origin1[0]) + - (origin2[1] - origin1[1]) * (origin2[1] - origin1[1])) + std::cout << "Sphere 2 center: (" << origin2[0] << ", " << origin2[1] << ")" << std::endl; + double distSq = 0.0; + for (int i = 0; i < D; ++i) + distSq += (origin2[i] - origin1[i]) * (origin2[i] - origin1[i]); + std::cout << "Expected geometric shift: " << std::sqrt(distSq) << std::endl; - ls::CompareChamfer compareChamfer(circle1, circle2); + ls::CompareChamfer compareChamfer(sphere1, sphere2); // Create output meshes with distance information auto targetMesh = ls::SmartPointer>::New(); @@ -121,20 +124,22 @@ int main() { std::cout << " Execution time: " << chamfer_ms.count() << " ms" << std::endl; // Save meshes with distance data - ls::VTKWriter(targetMesh, "chamfer_target_distances.vtp").apply(); - ls::VTKWriter(sampleMesh, "chamfer_sample_distances.vtp").apply(); + ls::VTKWriter(targetMesh, "chamfer_target_distances" + suffix) + .apply(); + ls::VTKWriter(sampleMesh, "chamfer_sample_distances" + suffix) + .apply(); // Test 2: Compare with other metrics std::cout << "\n=== Test 2: Comparison with Other Metrics ===" << std::endl; // Sparse Field comparison - auto circle1_expanded = ls::SmartPointer>::New(circle1); - ls::Expand(circle1_expanded, 50).apply(); - auto circle2_reduced = ls::SmartPointer>::New(circle2); - ls::Reduce(circle2_reduced, 1).apply(); + auto sphere1_expanded = ls::SmartPointer>::New(sphere1); + ls::Expand(sphere1_expanded, 50).apply(); + auto sphere2_reduced = ls::SmartPointer>::New(sphere2); + ls::Reduce(sphere2_reduced, 1).apply(); - ls::CompareSparseField compareSparseField(circle1_expanded, - circle2_reduced); + ls::CompareSparseField compareSparseField(sphere1_expanded, + sphere2_reduced); auto t3 = std::chrono::high_resolution_clock::now(); compareSparseField.apply(); auto t4 = std::chrono::high_resolution_clock::now(); @@ -147,47 +152,48 @@ int main() { << std::endl; std::cout << " Execution time: " << sparse_ms.count() << " ms" << std::endl; - // Area comparison - ls::CompareArea compareArea(circle1, circle2); + // Area/Volume comparison + ls::CompareVolume compareVolume(sphere1, sphere2); auto t5 = std::chrono::high_resolution_clock::now(); - compareArea.apply(); + compareVolume.apply(); auto t6 = std::chrono::high_resolution_clock::now(); std::chrono::duration area_ms = t6 - t5; - std::cout << "\nArea Comparison Results:" << std::endl; - std::cout << " Area mismatch: " << compareArea.getAreaMismatch() + std::cout << "\nArea/Volume Comparison Results:" << std::endl; + std::cout << " Area/Volume mismatch: " << compareVolume.getVolumeMismatch() + << std::endl; + std::cout << " Different cells: " << compareVolume.getCellCount() << std::endl; - std::cout << " Different cells: " << compareArea.getCellCount() << std::endl; std::cout << " Execution time: " << area_ms.count() << " ms" << std::endl; // Test 3: Different geometric configurations std::cout << "\n=== Test 3: Different Geometric Configurations ===" << std::endl; - // Test 3a: Identical circles (should give near-zero Chamfer distance) - auto circle3 = ls::SmartPointer>::New( + // Test 3a: Identical spheres (should give near-zero Chamfer distance) + auto sphere3 = ls::SmartPointer>::New( bounds, boundaryCons, gridDelta); ls::MakeGeometry( - circle3, ls::SmartPointer>::New(origin1, radius1)) + sphere3, ls::SmartPointer>::New(origin1, radius1)) .apply(); - ls::CompareChamfer compareIdentical(circle1, circle3); + ls::CompareChamfer compareIdentical(sphere1, sphere3); compareIdentical.apply(); - std::cout << "Identical circles:" << std::endl; + std::cout << "Identical spheres:" << std::endl; std::cout << " Chamfer distance: " << compareIdentical.getChamferDistance() << " (expected ~0)" << std::endl; // Test 3b: Different radii - auto circle4 = ls::SmartPointer>::New( + auto sphere4 = ls::SmartPointer>::New( bounds, boundaryCons, gridDelta); double radius4 = 7.0; // Larger radius ls::MakeGeometry( - circle4, ls::SmartPointer>::New(origin1, radius4)) + sphere4, ls::SmartPointer>::New(origin1, radius4)) .apply(); - ls::CompareChamfer compareDifferentSize(circle1, circle4); + ls::CompareChamfer compareDifferentSize(sphere1, sphere4); compareDifferentSize.apply(); std::cout << "\nDifferent radii (r1=" << radius1 << ", r2=" << radius4 @@ -201,15 +207,29 @@ int main() { std::cout << " Backward distance: " << compareDifferentSize.getBackwardDistance() << std::endl; + // Validation: Check if Chamfer distance is close to expected radius + // difference + double expectedDiff = std::abs(radius4 - radius1); + double actualChamfer = compareDifferentSize.getChamferDistance(); + double tolerance = 0.01 * expectedDiff; // 1% tolerance + if (std::abs(actualChamfer - expectedDiff) > tolerance) { + std::cerr << "ERROR: Chamfer distance " << actualChamfer + << " differs significantly from expected " << expectedDiff + << " (tolerance: " << tolerance << ")" << std::endl; + return; + } + std::cout << " Validation: PASSED (within 1% tolerance)" << std::endl; + // Test 3c: Large shift - auto circle5 = ls::SmartPointer>::New( + auto sphere5 = ls::SmartPointer>::New( bounds, boundaryCons, gridDelta); - double origin5[D] = {5., 0.}; // Larger shift + std::vector origin5(D, 0.0); + origin5[0] = 5.; // Larger shift ls::MakeGeometry( - circle5, ls::SmartPointer>::New(origin5, radius1)) + sphere5, ls::SmartPointer>::New(origin5, radius1)) .apply(); - ls::CompareChamfer compareLargeShift(circle1, circle5); + ls::CompareChamfer compareLargeShift(sphere1, sphere5); compareLargeShift.apply(); std::cout << "\nLarge shift (5 units in x-direction):" << std::endl; @@ -219,11 +239,20 @@ int main() { // Test 4: Performance summary std::cout << "\n=== Performance Summary ===" << std::endl; - std::cout << "Chamfer distance: " << chamfer_ms.count() << " ms" << std::endl; - std::cout << "Sparse field: " << sparse_ms.count() << " ms" << std::endl; - std::cout << "Area comparison: " << area_ms.count() << " ms" << std::endl; + std::cout << "Chamfer distance: " << chamfer_ms.count() << " ms" + << std::endl; + std::cout << "Sparse Field: " << sparse_ms.count() << " ms" + << std::endl; + std::cout << "Area/Volume: " << area_ms.count() << " ms" + << std::endl; + std::cout << "==========================" << std::endl; +} + +int main() { + omp_set_num_threads(8); - std::cout << "\n=== All Tests Completed Successfully ===" << std::endl; + runTest<2>(); + runTest<3>(); return 0; } diff --git a/tests/CompareCriticalDimensions/CompareCriticalDimensions.cpp b/tests/CompareCriticalDimensions/CompareCriticalDimensions.cpp index b3d0c0f0..be2488e3 100644 --- a/tests/CompareCriticalDimensions/CompareCriticalDimensions.cpp +++ b/tests/CompareCriticalDimensions/CompareCriticalDimensions.cpp @@ -1,4 +1,6 @@ #include +#include +#include #include #include @@ -16,16 +18,16 @@ namespace ls = viennals; -int main() { - constexpr int D = 2; - - omp_set_num_threads(4); - +template void runTest() { + std::cout << "Running " << D << "D Test..." << std::endl; double extent = 15; double gridDelta = 0.1; - double bounds[2 * D] = {-extent, extent, -extent, extent}; - ls::Domain::BoundaryType boundaryCons[D]; + double bounds[2 * D]; + for (int i = 0; i < 2 * D; ++i) + bounds[i] = (i % 2 == 0) ? -extent : extent; + + typename ls::Domain::BoundaryType boundaryCons[D]; for (unsigned i = 0; i < D; ++i) boundaryCons[i] = ls::Domain::BoundaryType::REFLECTIVE_BOUNDARY; @@ -33,7 +35,7 @@ int main() { auto circle1 = ls::SmartPointer>::New( bounds, boundaryCons, gridDelta); - double origin1[D] = {0., 0.}; + std::vector origin1(D, 0.0); double radius1 = 5.0; ls::MakeGeometry( @@ -44,39 +46,63 @@ int main() { auto circle2 = ls::SmartPointer>::New( bounds, boundaryCons, gridDelta); - double origin2[D] = {1.5, 0.5}; // Shifted center - double radius2 = 5.0; // Same radius + std::vector origin2(D, 0.0); + origin2[0] = 1.5; + origin2[1] = 0.5; + if (D > 2) + origin2[2] = 0.3; + double radius2 = 5.0; // Same radius ls::MakeGeometry( circle2, ls::SmartPointer>::New(origin2, radius2)) .apply(); + std::string suffix = "_" + std::to_string(D) + "D.vtp"; // Export both circles as VTK files for visualization { auto mesh = ls::SmartPointer>::New(); ls::ToSurfaceMesh(circle1, mesh).apply(); - ls::VTKWriter(mesh, "circle1_target.vtp").apply(); + ls::VTKWriter(mesh, "circle1_target" + suffix).apply(); } { auto mesh = ls::SmartPointer>::New(); ls::ToSurfaceMesh(circle2, mesh).apply(); - ls::VTKWriter(mesh, "circle2_sample.vtp").apply(); + ls::VTKWriter(mesh, "circle2_sample" + suffix).apply(); } // Compare critical dimensions ls::CompareCriticalDimensions compareCriticalDims(circle1, circle2); - // Add X ranges to find maximum and minimum Y positions (top and bottom) - // Search in the central X range where both circles overlap - compareCriticalDims.addXRange(-0, 0, true); // Find maximum Y (top) - compareCriticalDims.addXRange(-0, 0, false); // Find minimum Y (bottom) - - // Add Y ranges to find maximum and minimum X positions (right and left) - // Search in the central Y range where both circles overlap - compareCriticalDims.addYRange(-0, 0, true); // Find maximum X (right) - compareCriticalDims.addYRange(-0, 0, false); // Find minimum X (left) + if constexpr (D == 2) { + // Add X ranges to find maximum and minimum Y positions (top and bottom) + // Search in the central X range where both circles overlap + compareCriticalDims.addXRange(-0.1, 0.1, true); // Find maximum Y (top) + compareCriticalDims.addXRange(-0.1, 0.1, false); // Find minimum Y (bottom) + + // Add Y ranges to find maximum and minimum X positions (right and left) + // Search in the central Y range where both circles overlap + compareCriticalDims.addYRange(-0.1, 0.1, true); // Find maximum X (right) + compareCriticalDims.addYRange(-0.1, 0.1, false); // Find minimum X (left) + } else { + // For 3D, measure Z extent (top/bottom) at center (X=0, Y=0) + // And X extent (right/left) at center (Y=0, Z=0) + double inf = std::numeric_limits::max(); + double lowest = std::numeric_limits::lowest(); + + // Measure Z (dim 2) at X~0, Y~0 + std::array minB = {-0.1, -0.1, lowest}; + std::array maxB = {0.1, 0.1, inf}; + compareCriticalDims.addRange(2, minB, maxB, true); // Max Z + compareCriticalDims.addRange(2, minB, maxB, false); // Min Z + + // Measure X (dim 0) at Y~0, Z~0 + minB = {lowest, -0.1, -0.1}; + maxB = {inf, 0.1, 0.1}; + compareCriticalDims.addRange(0, minB, maxB, true); // Max X + compareCriticalDims.addRange(0, minB, maxB, false); // Min X + } // Create mesh for output auto mesh = ls::SmartPointer>::New(); @@ -86,7 +112,7 @@ int main() { compareCriticalDims.apply(); // Save mesh to file - ls::VTKWriter(mesh, "criticalDimensions.vtp").apply(); + ls::VTKWriter(mesh, "criticalDimensions" + suffix).apply(); // Debug: Print some surface mesh nodes to see actual positions std::cout << "\nDebug - Sample surface nodes from circle1:" << std::endl; @@ -95,18 +121,26 @@ int main() { std::cout << "Total nodes in circle1 surface: " << debugMesh1->nodes.size() << std::endl; for (size_t i = 0; i < std::min(size_t(10), debugMesh1->nodes.size()); ++i) { - std::cout << " Node " << i << ": (" << debugMesh1->nodes[i][0] << ", " - << debugMesh1->nodes[i][1] << ")" << std::endl; + std::cout << " Node " << i << ": ("; + for (int j = 0; j < D; ++j) + std::cout << debugMesh1->nodes[i][j] << (j == D - 1 ? "" : ", "); + std::cout << ")" << std::endl; } // Print results - std::cout << "Circle 1 center: (" << origin1[0] << ", " << origin1[1] << ")" - << std::endl; - std::cout << "Circle 2 center: (" << origin2[0] << ", " << origin2[1] << ")" - << std::endl; + std::cout << "Circle 1 center: ("; + for (int i = 0; i < D; ++i) + std::cout << origin1[i] << (i == D - 1 ? "" : ", "); + std::cout << ")" << std::endl; + std::cout << "Circle 2 center: ("; + for (int i = 0; i < D; ++i) + std::cout << origin2[i] << (i == D - 1 ? "" : ", "); + std::cout << ")" << std::endl; std::cout << "Radius: " << radius1 << std::endl; - std::cout << "Center shift: (" << (origin2[0] - origin1[0]) << ", " - << (origin2[1] - origin1[1]) << ")" << std::endl; + std::cout << "Center shift: ("; + for (int i = 0; i < D; ++i) + std::cout << (origin2[i] - origin1[i]) << (i == D - 1 ? "" : ", "); + std::cout << ")" << std::endl; std::cout << std::endl; // Get statistics @@ -114,6 +148,14 @@ int main() { std::cout << "Number of critical dimensions compared: " << numCriticalDims << std::endl; + // Compute shifts for expected value calculations + double xShift = origin2[0] - origin1[0]; + double yShift = origin2[1] - origin1[1]; + double zShift = (D > 2) ? origin2[2] - origin1[2] : 0.0; + + // Analytical differences based on geometry: + // For max: diff = shift + sqrt(r² - perp²) - r + // For min: diff = |shift - sqrt(r² - perp²) + r| std::cout << "\nIndividual critical dimension results:" << std::endl; for (size_t i = 0; i < numCriticalDims; ++i) { double posRef, posCmp, diff; @@ -123,6 +165,35 @@ int main() { std::cout << " Target position: " << posRef << std::endl; std::cout << " Sample position: " << posCmp << std::endl; std::cout << " Difference: " << diff << std::endl; + + // Compute expected difference based on dimension index + double expected = 0.0; + if constexpr (D == 2) { + double rEffectiveY = std::sqrt(radius1 * radius1 - xShift * xShift); + double rEffectiveX = std::sqrt(radius1 * radius1 - yShift * yShift); + if (i == 0) // max Y + expected = std::abs(yShift + rEffectiveY - radius1); + else if (i == 1) // min Y + expected = std::abs(yShift - rEffectiveY + radius1); + else if (i == 2) // max X + expected = std::abs(xShift + rEffectiveX - radius1); + else if (i == 3) // min X + expected = std::abs(xShift - rEffectiveX + radius1); + } else if constexpr (D == 3) { + double rEffectiveZ = + std::sqrt(radius1 * radius1 - xShift * xShift - yShift * yShift); + double rEffectiveX = + std::sqrt(radius1 * radius1 - yShift * yShift - zShift * zShift); + if (i == 0) // max Z + expected = std::abs(zShift + rEffectiveZ - radius1); + else if (i == 1) // min Z + expected = std::abs(zShift - rEffectiveZ + radius1); + else if (i == 2) // max X + expected = std::abs(xShift + rEffectiveX - radius1); + else if (i == 3) // min X + expected = std::abs(xShift - rEffectiveX + radius1); + } + std::cout << " Analytical: " << expected << std::endl; } else { std::cout << " Dimension " << i << ": Invalid (not found)" << std::endl; } @@ -135,24 +206,22 @@ int main() { << std::endl; std::cout << "RMSE: " << compareCriticalDims.getRMSE() << std::endl; - // Theoretical validation - // For circles with the same radius but shifted centers, - // the critical dimensions should differ by approximately the shift amount - // in the corresponding direction - std::cout << "\nExpected differences based on geometry:" << std::endl; - std::cout << "Top/Bottom (max/min Y): should be close to Y-shift = " - << std::abs(origin2[1] - origin1[1]) << std::endl; - std::cout << "Left/Right (max/min X): should be close to X-shift = " - << std::abs(origin2[0] - origin1[0]) << std::endl; - - // Additional test: Test with wider ranges - std::cout << "\n--- Testing with wider X range ---" << std::endl; - compareCriticalDims.clearRanges(); - compareCriticalDims.addXRange(-10, 10, true); // Find maximum Y - compareCriticalDims.addXRange(-10, 10, false); // Find minimum Y - compareCriticalDims.setOutputMesh(nullptr); // Don't create mesh - compareCriticalDims.apply(); + if constexpr (D == 2) { + // Additional test: Test with wider ranges + std::cout << "\n--- Testing with wider X range ---" << std::endl; + compareCriticalDims.clearRanges(); + compareCriticalDims.addXRange(-10, 10, true); // Find maximum Y + compareCriticalDims.addXRange(-10, 10, false); // Find minimum Y + compareCriticalDims.setOutputMesh(nullptr); // Don't create mesh + compareCriticalDims.apply(); + } else { + // Skip additional tests for 3D to keep it simple + return; + } + // With wide X range covering entire sphere, expected diff = yShift + std::cout << "Analytical difference (wide range): " << std::abs(yShift) + << std::endl; std::cout << "Number of critical dimensions: " << compareCriticalDims.getNumCriticalDimensions() << std::endl; for (size_t i = 0; i < compareCriticalDims.getNumCriticalDimensions(); ++i) { @@ -171,6 +240,9 @@ int main() { compareCriticalDims.addYRange(-10, 10, false); // Find minimum X compareCriticalDims.apply(); + // With wide Y range covering entire sphere, expected diff = xShift + std::cout << "Analytical difference (wide range): " << std::abs(xShift) + << std::endl; std::cout << "Number of critical dimensions: " << compareCriticalDims.getNumCriticalDimensions() << std::endl; for (size_t i = 0; i < compareCriticalDims.getNumCriticalDimensions(); ++i) { @@ -181,6 +253,12 @@ int main() { << std::endl; } } +} + +int main() { + omp_set_num_threads(4); + runTest<2>(); + runTest<3>(); return 0; } diff --git a/tests/CompareNarrowBand/CompareNarrowBand.cpp b/tests/CompareNarrowBand/CompareNarrowBand.cpp index 522f4b57..8efa10af 100644 --- a/tests/CompareNarrowBand/CompareNarrowBand.cpp +++ b/tests/CompareNarrowBand/CompareNarrowBand.cpp @@ -1,5 +1,7 @@ #include #include +#include +#include #include #include @@ -17,16 +19,16 @@ namespace ls = viennals; -int main() { - constexpr int D = 2; - - omp_set_num_threads(4); - +template void runTest() { + std::cout << "Running " << D << "D Test..." << std::endl; 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 (int i = 0; i < 2 * D; ++i) + bounds[i] = (i % 2 == 0) ? -extent : extent; + + typename ls::Domain::BoundaryType boundaryCons[D]; for (unsigned i = 0; i < D; ++i) boundaryCons[i] = ls::Domain::BoundaryType::REFLECTIVE_BOUNDARY; @@ -34,7 +36,7 @@ int main() { auto sphere1 = ls::SmartPointer>::New( bounds, boundaryCons, gridDelta); - double origin1[D] = {0., 0.}; + std::vector origin1(D, 0.0); double radius1 = 5.0; ls::MakeGeometry( @@ -45,24 +47,27 @@ int main() { auto sphere2 = ls::SmartPointer>::New( bounds, boundaryCons, gridDelta); - double origin2[D] = {2., 1.}; // Shifted center - double radius2 = 5.0; // Same radius + std::vector origin2(D, 0.0); + origin2[0] = 2.; + origin2[1] = 1.; + double radius2 = 5.0; // Same radius ls::MakeGeometry( sphere2, ls::SmartPointer>::New(origin2, radius2)) .apply(); + std::string suffix = "_" + std::to_string(D) + "D.vtp"; // Export both spheres as VTK files for visualization { auto mesh = ls::SmartPointer>::New(); ls::ToMesh(sphere1, mesh).apply(); - ls::VTKWriter(mesh, "sphere1_narrowband.vtp").apply(); + ls::VTKWriter(mesh, "sphere1_narrowband" + suffix).apply(); } { auto mesh = ls::SmartPointer>::New(); ls::ToMesh(sphere2, mesh).apply(); - ls::VTKWriter(mesh, "sphere2_narrowband.vtp").apply(); + ls::VTKWriter(mesh, "sphere2_narrowband" + suffix).apply(); } // Compare the narrow bands @@ -75,7 +80,9 @@ int main() { compareNarrowBand.apply(); // Save mesh to file - ls::VTKWriter(mesh, "narrowband_absolute_differences.vtu").apply(); + ls::VTKWriter(mesh, + "narrowband_absolute_differences" + suffix + ".vtu") + .apply(); // Get the calculated difference metrics double sumSquaredDifferences = compareNarrowBand.getSumSquaredDifferences(); @@ -123,13 +130,33 @@ int main() { std::cout << "Number of points in both ranges: " << compareNarrowBand.getNumPoints() << std::endl; + if constexpr (D == 3) { + // Test with restricted Z range + compareNarrowBand.clearXRange(); + compareNarrowBand.clearYRange(); + compareNarrowBand.setZRange(-5, 5); + compareNarrowBand.apply(); + std::cout << "RMSE with Z range [-5, 5]: " << compareNarrowBand.getRMSE() + << std::endl; + std::cout << "Number of points in Z range: " + << compareNarrowBand.getNumPoints() << std::endl; + compareNarrowBand.clearZRange(); + } + // Create a mesh output with squared differences compareNarrowBand.setOutputMesh(mesh); compareNarrowBand.setOutputMeshSquaredDifferences(true); compareNarrowBand.apply(); - ls::VTKWriter(mesh, - "narrowband_resctricted-range_squared_differences.vtu") + ls::VTKWriter(mesh, "narrowband_resctricted-range_squared_" + "differences" + + suffix + ".vtu") .apply(); +} + +int main() { + omp_set_num_threads(4); + runTest<2>(); + runTest<3>(); return 0; } diff --git a/tests/CompareSparseField/CompareSparseField.cpp b/tests/CompareSparseField/CompareSparseField.cpp index a43ddb43..f5c673be 100644 --- a/tests/CompareSparseField/CompareSparseField.cpp +++ b/tests/CompareSparseField/CompareSparseField.cpp @@ -175,14 +175,26 @@ template void runTest() { std::cout << "Number of points in X range: " << compareSparseField.getNumPoints() << std::endl; - // // Test with restricted Y range - // compareSparseField.clearXRange(); - // compareSparseField.setYRange(-5, 5); - // compareSparseField.apply(); - // std::cout << "RMSE with Y range [-5, 5]: " << compareSparseField.getRMSE() - // << std::endl; - // std::cout << "Number of points in Y range: " - // << compareSparseField.getNumPoints() << std::endl; + // Test with restricted Y range + compareSparseField.clearXRange(); + compareSparseField.setYRange(-5, 5); + compareSparseField.apply(); + std::cout << "RMSE with Y range [-5, 5]: " << compareSparseField.getRMSE() + << std::endl; + std::cout << "Number of points in Y range: " + << compareSparseField.getNumPoints() << std::endl; + + if constexpr (D == 3) { + // Test with restricted Z range + compareSparseField.clearYRange(); + compareSparseField.setZRange(-5, 5); + compareSparseField.apply(); + std::cout << "RMSE with Z range [-5, 5]: " << compareSparseField.getRMSE() + << std::endl; + std::cout << "Number of points in Z range: " + << compareSparseField.getNumPoints() << std::endl; + compareSparseField.clearZRange(); + } // Test with both X and Y range restrictions compareSparseField.setXRange(-3, 3); diff --git a/tests/CompareArea/CMakeLists.txt b/tests/CompareVolume/CMakeLists.txt similarity index 86% rename from tests/CompareArea/CMakeLists.txt rename to tests/CompareVolume/CMakeLists.txt index fdf14a4e..3311cd62 100644 --- a/tests/CompareArea/CMakeLists.txt +++ b/tests/CompareVolume/CMakeLists.txt @@ -1,4 +1,4 @@ -project(CompareArea LANGUAGES CXX) +project(CompareVolume LANGUAGES CXX) add_executable(${PROJECT_NAME} "${PROJECT_NAME}.cpp") target_link_libraries(${PROJECT_NAME} PRIVATE ViennaLS) diff --git a/tests/CompareVolume/CompareVolume.cpp b/tests/CompareVolume/CompareVolume.cpp new file mode 100644 index 00000000..04f803ff --- /dev/null +++ b/tests/CompareVolume/CompareVolume.cpp @@ -0,0 +1,179 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +/** + Test for lsCompareVolume that compares the area/volume difference between two + level sets. This test creates two different spheres and measures their + area/volume difference. \example CompareVolume.cpp +*/ + +namespace ls = viennals; + +template void runTest() { + std::cout << "Running " << D << "D Test..." << std::endl; + double extent = 15; + double gridDelta = 0.5; + + double bounds[2 * D]; + for (int i = 0; i < 2 * D; ++i) + bounds[i] = (i % 2 == 0) ? -extent : extent; + + typename ls::Domain::BoundaryType boundaryCons[D]; + for (unsigned i = 0; i < D; ++i) + boundaryCons[i] = ls::Domain::BoundaryType::REFLECTIVE_BOUNDARY; + + // Create first sphere (target) + auto sphere1 = ls::SmartPointer>::New( + bounds, boundaryCons, gridDelta); + + std::vector origin1(D, 0.0); + double radius1 = 5.0; + + ls::MakeGeometry( + sphere1, ls::SmartPointer>::New(origin1, radius1)) + .apply(); + + // Create second sphere (sample) with different radius + auto sphere2 = ls::SmartPointer>::New( + bounds, boundaryCons, gridDelta); + + std::vector origin2(D, 0.0); + double radius2 = 8.0; + + ls::MakeGeometry( + sphere2, ls::SmartPointer>::New(origin2, radius2)) + .apply(); + + std::string suffix = "_" + 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" + suffix + ".vtp").apply(); + } + + { + auto mesh = ls::SmartPointer>::New(); + ls::ToMesh(sphere2, mesh).apply(); + ls::VTKWriter(mesh, "sphere2" + suffix + ".vtp").apply(); + } + + // Compare the volumes/areas with mesh Output + auto mesh = ls::SmartPointer>::New(); // Create mesh for output + + if constexpr (D == 2) { + ls::CompareArea compareArea(sphere1, sphere2); + compareArea.setOutputMesh(mesh); + compareArea.apply(); + ls::VTKWriter(mesh, "volumeDifference" + suffix + ".vtu").apply(); + + // Area of circle = π * r² + double area1 = M_PI * radius1 * radius1; + double area2 = M_PI * radius2 * radius2; + double theoreticalDiff = std::abs(area2 - area1); + + double calculatedDifference = compareArea.getAreaMismatch(); + unsigned long int cellCount = compareArea.getCellCount(); + + std::cout << "Sphere 1 radius: " << radius1 << std::endl; + std::cout << "Sphere 2 radius: " << radius2 << std::endl; + std::cout << "Theoretical difference: " << theoreticalDiff << std::endl; + std::cout << "Calculated difference: " << calculatedDifference << std::endl; + std::cout << "Number of differing cells: " << cellCount << std::endl; + std::cout << "Error: " << std::abs(calculatedDifference - theoreticalDiff) + << std::endl; + + std::cout << "\nTesting custom increments and ranges:" << std::endl; + compareArea.setDefaultIncrement(2); + compareArea.apply(); + std::cout << "Difference with default increment of 2: " + << compareArea.getCustomAreaMismatch() << std::endl; + std::cout << "Cell count with default increment of 2: " + << compareArea.getCustomCellCount() << std::endl; + + compareArea.setDefaultIncrement(1); + compareArea.setXRangeAndIncrement(-5, 5, 3); + compareArea.apply(); + std::cout << "Difference with x-range increment of 3: " + << compareArea.getCustomAreaMismatch() << std::endl; + std::cout << "Cell count with x-range increment of 3: " + << compareArea.getCustomCellCount() << std::endl; + + compareArea.setDefaultIncrement(1); + compareArea.setYRangeAndIncrement(-5, 5, 4); + compareArea.apply(); + std::cout << "Difference with y-range increment of 4: " + << compareArea.getCustomAreaMismatch() << std::endl; + std::cout << "Cell count with y-range increment of 4: " + << compareArea.getCustomCellCount() << std::endl; + + } else { + ls::CompareVolume compareVolume(sphere1, sphere2); + compareVolume.setOutputMesh(mesh); + compareVolume.apply(); + ls::VTKWriter(mesh, "volumeDifference" + suffix + ".vtu").apply(); + + // Volume of sphere = 4/3 * π * r³ + double vol1 = (4.0 / 3.0) * M_PI * std::pow(radius1, 3); + double vol2 = (4.0 / 3.0) * M_PI * std::pow(radius2, 3); + double theoreticalDiff = std::abs(vol2 - vol1); + + double calculatedDifference = compareVolume.getVolumeMismatch(); + unsigned long int cellCount = compareVolume.getCellCount(); + + std::cout << "Sphere 1 radius: " << radius1 << std::endl; + std::cout << "Sphere 2 radius: " << radius2 << std::endl; + std::cout << "Theoretical difference: " << theoreticalDiff << std::endl; + std::cout << "Calculated difference: " << calculatedDifference << std::endl; + std::cout << "Number of differing cells: " << cellCount << std::endl; + std::cout << "Error: " << std::abs(calculatedDifference - theoreticalDiff) + << std::endl; + + std::cout << "\nTesting custom increments and ranges:" << std::endl; + compareVolume.setDefaultIncrement(2); + compareVolume.apply(); + std::cout << "Difference with default increment of 2: " + << compareVolume.getCustomVolumeMismatch() << std::endl; + std::cout << "Cell count with default increment of 2: " + << compareVolume.getCustomCellCount() << std::endl; + + compareVolume.setDefaultIncrement(1); + compareVolume.setXRangeAndIncrement(-5, 5, 3); + compareVolume.apply(); + std::cout << "Difference with x-range increment of 3: " + << compareVolume.getCustomVolumeMismatch() << std::endl; + std::cout << "Cell count with x-range increment of 3: " + << compareVolume.getCustomCellCount() << std::endl; + + compareVolume.setDefaultIncrement(1); + compareVolume.setYRangeAndIncrement(-5, 5, 4); + compareVolume.apply(); + std::cout << "Difference with y-range increment of 4: " + << compareVolume.getCustomVolumeMismatch() << std::endl; + std::cout << "Cell count with y-range increment of 4: " + << compareVolume.getCustomCellCount() << std::endl; + + compareVolume.setDefaultIncrement(1); + compareVolume.setZRangeAndIncrement(-5, 5, 5); + compareVolume.apply(); + std::cout << "Difference with z-range increment of 5: " + << compareVolume.getCustomVolumeMismatch() << std::endl; + std::cout << "Cell count with z-range increment of 5: " + << compareVolume.getCustomCellCount() << std::endl; + } +} + +int main() { + omp_set_num_threads(4); + runTest<2>(); + runTest<3>(); + return 0; +}