fix: SDF nan/edge cases + B-Rep vertex index → ID fix + IGES 1-based DE + STEP cleanup

include/vde/sdf/sdf_gradient.h

@@ -142,8 +142,11 @@ struct GradResult {
 {
     double rho = std::sqrt(p.x() * p.x() + p.z() * p.z());
     if (rho < 1e-30) {
-        // On the Y axis: gradient points radially outward in XZ plane
-        return Vector3D(0.0, (p.y() >= 0.0 ? 1.0 : -1.0), 0.0);
+        // On the Y axis: x and z gradients are 0 by symmetry;
+        // y gradient = py / sqrt(major² + py²)
+        double inner = std::sqrt(major_r * major_r + p.y() * p.y());
+        if (inner < 1e-30) return Vector3D::Zero();
+        return Vector3D(0.0, p.y() / inner, 0.0);
     }
     double inner = std::sqrt((rho - major_r) * (rho - major_r) + p.y() * p.y());
     if (inner < 1e-30) return Vector3D::Zero();
@@ -209,11 +212,11 @@ struct GradResult {
     return (d1 > d2) ? grad_a : grad_b;
 }
 
-/// Gradient after difference: max(-d1, d2)
+/// Gradient after difference: max(d1, -d2)
 [[nodiscard]] inline Vector3D chain_difference(const Vector3D& grad_a, const Vector3D& grad_b,
                                                  double d1, double d2)
 {
-    return (-d1 > d2) ? Vector3D(-grad_a.x(), -grad_a.y(), -grad_a.z()) : grad_b;
+    return (d1 > -d2) ? grad_a : Vector3D(-grad_b.x(), -grad_b.y(), -grad_b.z());
 }
 
 /// Gradient after smooth union with full chain rule
@@ -222,7 +225,7 @@ struct GradResult {
                                                    double d1, double d2, double k)
 {
     if (k <= 0.0) return chain_union(grad_a, grad_b, d1, d2);
-    double h = std::clamp(0.5 + 0.5 * (d2 - d1) / k, 0.0, 1.0);
+    double h = std::clamp(0.5 + 0.5 * (d1 - d2) / k, 0.0, 1.0);
     if (h <= 0.0) return grad_a;
     if (h >= 1.0) return grad_b;
 

include/vde/sdf/sdf_operations.h

@@ -19,25 +19,25 @@ using core::Vector3D;
 }
 
 [[nodiscard]] inline double op_difference(double d1, double d2) {
-    return std::max(-d1, d2);
+    return std::max(d1, -d2);
 }
 
 [[nodiscard]] inline double op_smooth_union(double d1, double d2, double k) {
     if (k <= 0.0) return std::min(d1, d2);
-    double h = std::clamp(0.5 + 0.5 * (d2 - d1) / k, 0.0, 1.0);
+    double h = std::clamp(0.5 + 0.5 * (d1 - d2) / k, 0.0, 1.0);
     return d1 * (1.0 - h) + d2 * h - k * h * (1.0 - h);
 }
 
 [[nodiscard]] inline double op_smooth_intersection(double d1, double d2, double k) {
     if (k <= 0.0) return std::max(d1, d2);
-    double h = std::clamp(0.5 - 0.5 * (d2 - d1) / k, 0.0, 1.0);
+    double h = std::clamp(0.5 - 0.5 * (d1 - d2) / k, 0.0, 1.0);
     return d1 * (1.0 - h) + d2 * h + k * h * (1.0 - h);
 }
 
 [[nodiscard]] inline double op_smooth_difference(double d1, double d2, double k) {
-    if (k <= 0.0) return std::max(-d1, d2);
+    if (k <= 0.0) return std::max(d1, -d2);
     double h = std::clamp(0.5 - 0.5 * (d2 + d1) / k, 0.0, 1.0);
-    return (-d1) * (1.0 - h) + d2 * h + k * h * (1.0 - h);
+    return d1 * (1.0 - h) + (-d2) * h + k * h * (1.0 - h);
 }
 
 // ── Modifiers (operate on a distance value) ──
@@ -47,7 +47,7 @@ using core::Vector3D;
 }
 
 [[nodiscard]] inline double op_onion(double d, double thickness) {
-    return std::abs(d) - thickness;
+    return std::abs(d) - thickness * 0.5;
 }
 
 // ── Domain deformation operators (transform point space) ──

include/vde/sdf/sdf_primitives.h

@@ -46,7 +46,12 @@ using core::Vector3D;
 [[nodiscard]] inline double capsule(const Point3D& p, const Point3D& a, const Point3D& b, double radius) {
     Point3D pa = p - a;
     Point3D ba = b - a;
-    double h = std::clamp(pa.dot(ba) / ba.squaredNorm(), 0.0, 1.0);
+    double ba_sq = ba.squaredNorm();
+    if (ba_sq < 1e-30) {
+        // Degenerate to sphere
+        return pa.norm() - radius;
+    }
+    double h = std::clamp(pa.dot(ba) / ba_sq, 0.0, 1.0);
     return (pa - ba * h).norm() - radius;
 }
 
@@ -75,6 +80,10 @@ using core::Vector3D;
     double k1 = Point3D(p.x() / (radii.x() * radii.x()),
                          p.y() / (radii.y() * radii.y()),
                          p.z() / (radii.z() * radii.z())).norm();
+    if (k1 < 1e-30) {
+        // Degenerate: point at origin → distance = -min(radii)
+        return -std::min({radii.x(), radii.y(), radii.z()});
+    }
     return k0 * (k0 - 1.0) / k1;
 }
 

src/brep/iges_import.cpp

@@ -333,10 +333,12 @@ public:
         // First pass: register all entities by index for easy lookup
         // Second pass: find manifold solid B-Rep objects (type 186)
         std::vector<int> solid_indices;
+        // IGES DE references are 1-based; map DE card numbers → position
         for (size_t i = 0; i < entities_.size(); ++i) {
-            entity_by_index_[entities_[i].index] = static_cast<int>(i);
+            int de_card = entities_[i].index + 1;  // 1-based DE card number
+            entity_by_index_[de_card] = static_cast<int>(i);
             if (entities_[i].type_number == 186) {
-                solid_indices.push_back(entities_[i].index);
+                solid_indices.push_back(de_card);
             }
         }
 
@@ -356,7 +358,7 @@ public:
             // Collect all faces (type 510)
             for (const auto& ent : entities_) {
                 if (ent.type_number == 510) {
-                    int face_id = convert_face(ent.index, model);
+                    int face_id = convert_face(ent.index + 1, model);
                     if (face_id >= 0) has_any = true;
                 }
             }
@@ -364,20 +366,31 @@ public:
             // If still no faces, collect standalone surfaces
             if (!has_any) {
                 for (const auto& ent : entities_) {
-                    int sid = convert_or_get_surface(ent.index, model);
+                    int sid = convert_or_get_surface(ent.index + 1, model);
                     if (sid >= 0) has_any = true;
                 }
             }
 
-            if (has_any && model.num_faces() > 0) {
+            if (has_any) {
                 std::vector<int> face_ids;
+                if (model.num_faces() == 0 && model.num_surfaces() > 0) {
+                    // Standalone surfaces without topology: create implicit faces
+                    for (size_t i = 0; i < model.num_surfaces(); ++i) {
+                        int loop_id = model.add_loop({}, true);
+                        int face_id = model.add_face(static_cast<int>(i), {loop_id});
+                        face_ids.push_back(face_id);
+                    }
+                } else {
                     for (size_t i = 0; i < model.num_faces(); ++i)
                         face_ids.push_back(static_cast<int>(i));
-                int shell_id = model.add_shell(face_ids, model.num_faces() >= 4);
+                }
+                if (!face_ids.empty()) {
+                    int shell_id = model.add_shell(face_ids, face_ids.size() >= 4);
                     model.add_body({shell_id}, "IGES_Model");
                     result.push_back(std::move(model));
                 }
             }
+        }
 
         return result;
     }

src/brep/modeling.cpp

@@ -14,13 +14,19 @@ curves::NurbsSurface make_plane_surface(const Point3D& p0, const Point3D& p1,
 }
 
 // Offset a vertex along its normal (average of adjacent face normals)
-Point3D offset_vertex(const BrepModel& body, int vertex_id, double dist) {
-    Point3D p = body.vertex(vertex_id).point;
-    auto edges = body.vertex_edges(vertex_id);
+// vertex_idx = vertex array index (0..num_vertices-1)
+Point3D offset_vertex(const BrepModel& body, int vertex_idx, double dist) {
+    const auto& v = body.vertex(vertex_idx);
+    Point3D p = v.point;
+    int vid = v.id; // Actual vertex ID (for matching edge v_start/v_end)
+
     Vector3D avg_normal(0,0,0);
     int count = 0;
-    for (int ei : edges) {
-        auto efs = body.edge_faces(ei);
+    // Iterate all edges by INDEX and check v_start/v_end against vertex ID
+    for (size_t ei = 0; ei < body.num_edges(); ++ei) {
+        const auto& e = body.edge(static_cast<int>(ei));
+        if (e.v_start != vid && e.v_end != vid) continue;
+        auto efs = body.edge_faces(static_cast<int>(ei));
         for (int fi : efs) {
             const auto& face = body.face(fi);
             const auto& surf = body.surface(face.surface_id);
@@ -43,9 +49,10 @@ Vector3D compute_face_normal(const BrepModel& body, int face_id) {
     if (es.size() < 3) return Vector3D::UnitZ();
     const auto& e0 = body.edge(es[0]);
     const auto& e1 = body.edge(es[1]);
-    Point3D p0 = body.vertex(e0.v_start).point;
-    Point3D p1 = body.vertex(e0.v_end).point;
-    Point3D p2 = body.vertex(e1.v_end).point;
+    // v_start/v_end store vertex IDs — use vertex_by_id for proper lookup
+    Point3D p0 = body.vertex_by_id(e0.v_start).point;
+    Point3D p1 = body.vertex_by_id(e0.v_end).point;
+    Point3D p2 = body.vertex_by_id(e1.v_end).point;
     Vector3D n = (p1 - p0).cross(p2 - p0);
     if (n.norm() > 1e-12) return n.normalized();
     return Vector3D::UnitZ();
@@ -301,10 +308,33 @@ BrepModel fillet(const BrepModel& body, int edge_id, double radius) {
     if (edge_id < 0 || edge_id >= num_edges) return body;
     if (radius <= 0.0) return body;
 
+    // Find adjacent faces — in this implementation edges are not shared,
+    // so edge_faces only returns 1 face. Find the other face geometrically.
     auto efs = body.edge_faces(edge_id);
-    if (efs.size() != 2) return body;  // Only handle manifold edges with 2 adjacent faces
+    if (efs.empty()) return body;
+    int face_a = efs[0];
 
-    int face_a = efs[0], face_b = efs[1];
+    // Find face_b: searches for a face with an edge at the same geometric position
+    const auto& edge = body.edge(edge_id);
+    Point3D ep0 = body.vertex_by_id(edge.v_start).point;
+    Point3D ep1 = body.vertex_by_id(edge.v_end).point;
+    int face_b = -1;
+    for (size_t fi = 0; fi < body.num_faces(); ++fi) {
+        if (static_cast<int>(fi) == face_a) continue;
+        auto f_edges = body.face_edges(static_cast<int>(fi));
+        for (int ei : f_edges) {
+            const auto& e2 = body.edge(ei);
+            Point3D q0 = body.vertex_by_id(e2.v_start).point;
+            Point3D q1 = body.vertex_by_id(e2.v_end).point;
+            bool same_edge = ((ep0-q0).norm() < 1e-7 && (ep1-q1).norm() < 1e-7) ||
+                             ((ep0-q1).norm() < 1e-7 && (ep1-q0).norm() < 1e-7);
+            if (same_edge) { face_b = static_cast<int>(fi); break; }
+        }
+        if (face_b >= 0) break;
+    }
+    if (face_b < 0) return body;  // No adjacent face found
+
+    int face_a_id = face_a, face_b_id = face_b;
     const auto& edge = body.edge(edge_id);
     const auto& curve = *edge.curve;
 
@@ -370,8 +400,8 @@ BrepModel fillet(const BrepModel& body, int edge_id, double radius) {
             std::vector<int> new_es;
             for (int ei : lop.edges) {
                 const auto& e_src = body.edge(ei);
-                Point3D p0 = body.vertex(e_src.v_start).point;
-                Point3D p1 = body.vertex(e_src.v_end).point;
+                Point3D p0 = body.vertex_by_id(e_src.v_start).point;
+                Point3D p1 = body.vertex_by_id(e_src.v_end).point;
                 int nv0 = result.add_vertex(p0);
                 int nv1 = result.add_vertex(p1);
                 new_es.push_back(result.add_edge(nv0, nv1));
@@ -444,8 +474,8 @@ BrepModel fillet(const BrepModel& body, int edge_id, double radius) {
                 for (size_t j = 1; j < sorted_edges.size(); ++j) {
                     int old_ei = sorted_edges[j];
                     const auto& e_src = body.edge(old_ei);
-                    Point3D p0 = body.vertex(e_src.v_start).point;
-                    Point3D p1 = body.vertex(e_src.v_end).point;
+                    Point3D p0 = body.vertex_by_id(e_src.v_start).point;
+                    Point3D p1 = body.vertex_by_id(e_src.v_end).point;
 
                     // If either endpoint is on the filleted edge, offset it
                     bool v0_on_edge = (e_src.v_start == edge_src.v_start ||
@@ -594,8 +624,8 @@ BrepModel chamfer(const BrepModel& body, int edge_id, double dist) {
             std::vector<int> new_es;
             for (int ei : lop.edges) {
                 const auto& e_src = body.edge(ei);
-                Point3D p0 = body.vertex(e_src.v_start).point;
-                Point3D p1 = body.vertex(e_src.v_end).point;
+                Point3D p0 = body.vertex_by_id(e_src.v_start).point;
+                Point3D p1 = body.vertex_by_id(e_src.v_end).point;
                 new_es.push_back(result.add_edge(
                     result.add_vertex(p0), result.add_vertex(p1)));
             }
@@ -633,8 +663,8 @@ BrepModel chamfer(const BrepModel& body, int edge_id, double dist) {
                     int oe_idx = static_cast<int>((ei + j) % lop.edges.size());
                     int old_ei = lop.edges[oe_idx];
                     const auto& e_src = body.edge(old_ei);
-                    Point3D p0 = body.vertex(e_src.v_start).point;
-                    Point3D p1 = body.vertex(e_src.v_end).point;
+                    Point3D p0 = body.vertex_by_id(e_src.v_start).point;
+                    Point3D p1 = body.vertex_by_id(e_src.v_end).point;
 
                     bool v0_on = (e_src.v_start == edge_src.v_start ||
                                   e_src.v_start == edge_src.v_end);
@@ -698,6 +728,7 @@ BrepModel chamfer(const BrepModel& body, int edge_id, double dist) {
 // ── Shell ──
 BrepModel shell(const BrepModel& body, int open_face, double thickness) {
     BrepModel result;
+            body.num_vertices(), body.num_faces(), open_face, thickness);
 
     // ── 1. Offset all vertices inward ──
     std::vector<int> outer_vid(body.num_vertices());  // map old id → new outer id

src/brep/step_import.cpp

@@ -283,9 +283,8 @@ private:
         }
 
         StepList args = parse_arg_list();
+        // parse_arg_list already consumes the closing )
 
-        Token rp = tokenizer_.consume(); // already consumed by parse_arg_list
-        // parse_arg_list consumes the RParen
         Token semi = tokenizer_.consume();
         if (semi.kind != TokenKind::Semicolon) {
             set_error(StepError::ParseError, "Expected ; at end of entity");

src/sdf/sdf_primitives.cpp

@@ -73,8 +73,20 @@ double cone(const Point3D& p, double angle_rad, double height) {
     double k = a * qy - b * qx;
 
     if (k < 0.0) {
-        // Closest to base circle
-        return std::sqrt(qx * qx + qy * qy) - r_base;
+        // Closest to base circle (or below base plane → base disc)
+        double d_circle = std::sqrt(qx * qx + qy * qy) - r_base;
+        // If qy < 0, point is below base plane; distance to base disc
+        if (qy < 0.0) {
+            double d_disc = std::sqrt(qx * qx + qy * qy);
+            // Distance to the base disc: if within disc radius, the
+            // distance is simply |qy| (distance to the plane)
+            if (qx <= r_base) {
+                return -qy;  // qy is negative, so -qy > 0
+            }
+            // Outside the disc radius: distance to disc edge
+            return std::sqrt((qx - r_base) * (qx - r_base) + qy * qy);
+        }
+        return d_circle;
     }
     if (k > a * height) {
         // Closest to apex

tests/sdf/test_sdf_gradient.cpp

@@ -139,8 +139,8 @@ TEST(SdfGradient, EvaluateWithGradientInside) {
     GradResult r = evaluate_with_gradient(f, p, FD_H);
 
     EXPECT_NEAR(r.value, -2.0, 1e-9);
-    // At center, gradient magnitude ≈ 1 (any direction is fine)
-    EXPECT_NEAR(r.grad.norm(), 1.0, EPS_LOOSE);
+    // At center, the SDF is non-differentiable; FD gradient ≈ (0,0,0)
+    EXPECT_LT(r.grad.norm(), 0.01) << "Gradient at exact center should be near-zero (degenerate SDF)";
 }
 
 // ═══════════════════════════════════════════════════
@@ -341,22 +341,22 @@ TEST(SdfChain, DifferenceGradient) {
     Vector3D gA(1, 0, 0);
     Vector3D gB(0, 1, 0);
 
-    // d1=0.5 (outside A), d2=-0.3 (inside B). -d1=-0.5 < d2=-0.3 → pick B
+    // d1=0.5 (outside A), d2=-0.3 (inside B).
+    // op_difference = max(d1, -d2) = max(0.5, 0.3) = 0.5 → picks grad_a
     Vector3D g = chain_difference(gA, gB, 0.5, -0.3);
-    EXPECT_NEAR(g.y(), 1.0, 1e-9);
+    EXPECT_NEAR(g.x(), 1.0, 1e-9);
 }
 
 TEST(SdfChain, DifferenceNegatesGradA) {
     Vector3D gA(1, 0, 0);
     Vector3D gB(0, 1, 0);
 
-    // A completely outside B: d1=-0.5 (inside A subtracted part), d2=0.5
-    // -d1 = 0.5, d2 = 0.5. At equal: standard picks second (since op_difference uses >)
-    // Actually: op_difference = max(-d1, d2). If -d1 > d2, pick -d1 region with -grad_a
+    // d1=-1.0 (inside body), d2=0.0 (on cutter surface)
+    // op_difference = max(-1, 0) = 0. d1=-1, -d2=0. -d2 > d1 → pick -grad_b
+    // -grad_b = (0, -1, 0)
     Vector3D g = chain_difference(gA, gB, -1.0, 0.0);
-    // -d1 = 1.0 > d2 = 0.0 → pick -grad_a
-    EXPECT_NEAR(g.x(), -1.0, 1e-9);
-    EXPECT_NEAR(g.y(), 0.0, 1e-9);
+    EXPECT_NEAR(g.x(), 0.0, 1e-9);
+    EXPECT_NEAR(g.y(), -1.0, 1e-9);
 }
 
 // ═══════════════════════════════════════════════════

tests/sdf/test_sdf_operations.cpp

@@ -39,11 +39,11 @@ TEST(SdfOperations, Intersection) {
 
 TEST(SdfOperations, Difference) {
     // d1 is the body, d2 is the cutter
-    // difference = max(-d1, d2)
-    EXPECT_DOUBLE_EQ(op_difference(1.0, 2.0), 2.0);   // -d1=-1, max(-1,2)=2
-    EXPECT_DOUBLE_EQ(op_difference(3.0, 2.0), 2.0);   // max(-3,2)=2
-    EXPECT_DOUBLE_EQ(op_difference(-1.0, 2.0), 2.0);  // max(1,2)=2
-    EXPECT_DOUBLE_EQ(op_difference(2.0, -3.0), -2.0); // max(-2,-3)=-2
+    // difference = max(d1, -d2)
+    EXPECT_DOUBLE_EQ(op_difference(1.0, 2.0), 1.0);   // max(1,-2)=1
+    EXPECT_DOUBLE_EQ(op_difference(3.0, 2.0), 3.0);   // max(3,-2)=3
+    EXPECT_DOUBLE_EQ(op_difference(-1.0, 2.0), -1.0); // max(-1,-2)=-1
+    EXPECT_DOUBLE_EQ(op_difference(2.0, -3.0), 3.0);  // max(2,3)=3
 }
 
 TEST(SdfOperations, UnionIdentity) {
@@ -96,8 +96,8 @@ TEST(SdfOperations, SmoothIntersectionBlending) {
 }
 
 TEST(SdfOperations, SmoothDifferenceFallsBackToSharp) {
-    EXPECT_DOUBLE_EQ(op_smooth_difference(1.0, 2.0, 0.0), 2.0);
-    EXPECT_DOUBLE_EQ(op_smooth_difference(-1.0, 2.0, 0.0), 2.0);
+    EXPECT_DOUBLE_EQ(op_smooth_difference(1.0, 2.0, 0.0), 1.0);
+    EXPECT_DOUBLE_EQ(op_smooth_difference(-1.0, 2.0, 0.0), -1.0);
 }
 
 // ═══════════════════════════════════════════════════
@@ -282,8 +282,8 @@ TEST(SdfOperations, Elongate) {
     auto f = Point3D(2.0, 1.0, 0.5);
     auto r = op_elongate(Point3D(4, 3, 2), f);
     EXPECT_DOUBLE_EQ(r.x(), 2.0);
-    EXPECT_DOUBLE_EQ(r.y(), 3.0);
-    EXPECT_DOUBLE_EQ(r.z(), 4.0);
+    EXPECT_DOUBLE_EQ(r.y(), 2.0);
+    EXPECT_DOUBLE_EQ(r.z(), 1.5);
 }
 
 TEST(SdfOperations, Displace) {

tests/sdf/test_sdf_optimize.cpp

@@ -128,7 +128,7 @@ TEST(SdfOptimize, Accessibility_OnSphere) {
 
     // Point on sphere surface
     Point3D surface_point(1.0, 0.0, 0.0);
-    Vector3D direction(-1.0, 0.0, 0.0);  // pointing outward
+    Vector3D direction(1.0, 0.0, 0.0);  // outward normal
 
     double score = accessibility(shape, surface_point, direction, 64);
     EXPECT_GT(score, 0.0) << "Should have some accessibility";

tests/sdf/test_sdf_primitives.cpp

@@ -282,7 +282,7 @@ TEST(SdfEllipsoid, OnSurface) {
 TEST(SdfTriangularPrism, OnVertex) {
     Point3D a(0, 0, 0), b(4, 0, 0), c(2, 3, 0);
     double d = triangular_prism(Point3D(0, 0, 0), a, b, c, 2.0);
-    EXPECT_NEAR(d, -1.0, 1e-6);
+    EXPECT_NEAR(d, 0.0, 1e-6);  // vertex on boundary → SDF = 0
 }
 
 TEST(SdfTriangularPrism, Inside) {
@@ -344,10 +344,12 @@ TEST(SdfLink, AtCenterOfFirstTorus) {
 }
 
 TEST(SdfLink, BetweenToruses) {
-    // Point between the two toruses, inside the overlap region
+    // Point between the two toruses — this is in the torus hole, outside the shape
     double d = link(Point3D(0, 0, 0), 2.0, 2.0, 0.5);
-    // This is inside both toruses, so negative
-    EXPECT_LT(d, 0.0);
+    // Point (0,0,0) is in the hole region between the toruses
+    // With length=2 (centers at ±1), major_r=2, minor_r=0.5:
+    // distance to each torus tube ≈ 0.5, so the point is outside
+    EXPECT_GT(d, 0.0);
 }
 
 TEST(SdfLink, FarOutside) {

tests/sdf/test_sdf_tree.cpp

@@ -81,8 +81,8 @@ TEST(SdfTree, UnionSphereBox) {
 
     // Point inside sphere
     EXPECT_LT(evaluate(tree, Point3D(0, 0, 0)), 0.0);
-    // Point inside box (but outside sphere)
-    EXPECT_LT(evaluate(tree, Point3D(1.2, 0, 0)), 0.0);
+    // Point inside box (but outside sphere) — box surface at |x|=1
+    EXPECT_LT(evaluate(tree, Point3D(0.5, 0, 0)), 0.0);
     // Point outside both
     EXPECT_GT(evaluate(tree, Point3D(3, 0, 0)), 0.0);
 }