From e0c0f1fef7554b1d06aafcc943210b480dc4ba7f Mon Sep 17 00:00:00 2001
From: El Hadi Moussi <moussi@phimeca.com>
Date: Tue, 6 Aug 2024 17:55:26 +0200
Subject: [PATCH] Check if the neighbors can be empty

If there is double nodes some nodes don't have neighbors
---
 src/ShapeRecogn/NodeCurvatureCalculator.cxx | 96 +++++++++++----------
 1 file changed, 50 insertions(+), 46 deletions(-)

diff --git a/src/ShapeRecogn/NodeCurvatureCalculator.cxx b/src/ShapeRecogn/NodeCurvatureCalculator.cxx
index 46b545415..22fda87ae 100644
--- a/src/ShapeRecogn/NodeCurvatureCalculator.cxx
+++ b/src/ShapeRecogn/NodeCurvatureCalculator.cxx
@@ -48,8 +48,6 @@ void NodeCurvatureCalculator::computeNormals()
     mesh->getReverseNodalConnectivity(revNodal, revNodalIdx);
     for (mcIdType nodeId = 0; nodeId < nbNodes; nodeId++)
     {
-        // nodeIds[0] = nodeId;
-        // auto cellIds = mesh->getCellIdsLyingOnNodes(nodeIds.data(), nodeIds.data() + 1, false);
         int nbCells = revNodalIdx->getIJ(nodeId + 1, 0) - revNodalIdx->getIJ(nodeId, 0);
         std::vector<mcIdType> cellIds(nbCells, 0);
         int start = revNodalIdx->getIJ(nodeId, 0);
@@ -92,55 +90,61 @@ void NodeCurvatureCalculator::computeCurvatures(mcIdType nodeId, double tol)
 {
     std::array<double, 3> normal = nodes->getNormal(nodeId);
     std::vector<mcIdType> neighborIds = nodes->getNeighbors(nodeId);
-    std::vector<double> discreteCurvatures = computeDiscreteCurvatures(nodeId, neighborIds);
-    std::vector<double> tangents = computeTangentDirections(nodeId, neighborIds);
-    mcIdType maxCurvatureId = std::distance(
-        discreteCurvatures.begin(),
-        std::max_element(discreteCurvatures.begin(), discreteCurvatures.end()));
-    std::array<double, 3> e1{
-        tangents[3 * maxCurvatureId],
-        tangents[3 * maxCurvatureId + 1],
-        tangents[3 * maxCurvatureId + 2]};
-    std::array<double, 3> e2 = MathOps::normalize(MathOps::cross(e1, normal));
-    std::vector<double> coeffs = computeNormalCurvatureCoefficients(
-        discreteCurvatures, tangents, normal, e1);
-    double a = coeffs[0], b = coeffs[1], c = coeffs[2];
-    double h = (a + c) / 2.0;
-    double kg = a * c - pow(b, 2) / 4.0;
-    double k1 = h + sqrt(fabs(pow(h, 2) - kg));
-    double k2 = h - sqrt(fabs(pow(h, 2) - kg));
     double theta0 = 0.0;
-    if (fabs(k1 - k2) >= tol)
-        theta0 = 0.5 * asin(b / (k1 - k2));
-    std::array<double, 3> direction1{0.0, 0.0, 0.0};
-    std::array<double, 3> direction2{0.0, 0.0, 0.0};
-    for (size_t i = 0; i < 3; ++i)
-    {
-        direction1[i] = cos(theta0) * e1[i] + sin(theta0) * e2[i];
-        direction2[i] = -sin(theta0) * e1[i] + cos(theta0) * e2[i];
-    }
-    bool isConvex;
-    double kdiff0;
-    double kis0;
+    bool isConvex = false;
+    double k1 = 0.0;
+    double k2 = 0.0;
+    double kdiff0 = 0.0;
+    double kis0 = 0.0;
+    PrimitiveType primitive = PrimitiveType::Unknown;
     std::array<double, 3> mainDirection{0.0, 0.0, 0.0};
     std::array<double, 3> weakDirection{0.0, 0.0, 0.0};
-    if (fabs(k1) > fabs(k2))
-    {
-        isConvex = true;
-        kdiff0 = k1;
-        kis0 = k2;
-        mainDirection = direction1;
-        weakDirection = direction2;
-    }
-    else
+    if (neighborIds.size() > 0)
     {
-        isConvex = false;
-        kdiff0 = k2;
-        kis0 = k1;
-        mainDirection = direction2;
-        weakDirection = direction1;
+        std::vector<double> discreteCurvatures = computeDiscreteCurvatures(nodeId, neighborIds);
+        std::vector<double> tangents = computeTangentDirections(nodeId, neighborIds);
+        mcIdType maxCurvatureId = std::distance(
+            discreteCurvatures.begin(),
+            std::max_element(discreteCurvatures.begin(), discreteCurvatures.end()));
+        std::array<double, 3> e1{
+            tangents[3 * maxCurvatureId],
+            tangents[3 * maxCurvatureId + 1],
+            tangents[3 * maxCurvatureId + 2]};
+        std::array<double, 3> e2 = MathOps::normalize(MathOps::cross(e1, normal));
+        std::vector<double> coeffs = computeNormalCurvatureCoefficients(
+            discreteCurvatures, tangents, normal, e1);
+        double a = coeffs[0], b = coeffs[1], c = coeffs[2];
+        double h = (a + c) / 2.0;
+        double kg = a * c - pow(b, 2) / 4.0;
+        k1 = h + sqrt(fabs(pow(h, 2) - kg));
+        k2 = h - sqrt(fabs(pow(h, 2) - kg));
+        if (fabs(k1 - k2) >= tol)
+            theta0 = 0.5 * asin(b / (k1 - k2));
+        std::array<double, 3> direction1{0.0, 0.0, 0.0};
+        std::array<double, 3> direction2{0.0, 0.0, 0.0};
+        for (size_t i = 0; i < 3; ++i)
+        {
+            direction1[i] = cos(theta0) * e1[i] + sin(theta0) * e2[i];
+            direction2[i] = -sin(theta0) * e1[i] + cos(theta0) * e2[i];
+        }
+        if (fabs(k1) > fabs(k2))
+        {
+            isConvex = true;
+            kdiff0 = k1;
+            kis0 = k2;
+            mainDirection = direction1;
+            weakDirection = direction2;
+        }
+        else
+        {
+            isConvex = false;
+            kdiff0 = k2;
+            kis0 = k1;
+            mainDirection = direction2;
+            weakDirection = direction1;
+        }
+        primitive = findPrimitiveType(k1, k2, kdiff0, kis0);
     }
-    PrimitiveType primitive = findPrimitiveType(k1, k2, kdiff0, kis0);
     // Populate nodes
     nodes->isConvex[nodeId] = isConvex;
     nodes->theta0[nodeId] = theta0;
-- 
2.39.2