#include "vde/brep/modeling.h"
#include <cmath>
#include <map>
#include <algorithm>
#include <cassert>

namespace vde::brep {

namespace {

curves::NurbsSurface make_plane_surface(const Point3D& p0, const Point3D& p1,
                                         const Point3D& p2, const Point3D& p3) {
    std::vector<std::vector<Point3D>> grid = {{p0, p3}, {p1, p2}};
    return curves::NurbsSurface(grid, {0,0,1,1}, {0,0,1,1}, {}, 1, 1);
}

// Offset a vertex along its normal (average of adjacent face normals)
// 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;
    // 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);
            auto du = surf.derivative_u(0.5, 0.5);
            auto dv = surf.derivative_v(0.5, 0.5);
            Vector3D n = du.cross(dv);
            if (n.norm() > 1e-12) {
                avg_normal += n.normalized();
                count++;
            }
        }
    }
    if (count > 0) avg_normal /= count;
    return p + avg_normal * dist;
}

// Compute a representative face normal from its first 3 vertices
Vector3D compute_face_normal(const BrepModel& body, int face_id) {
    auto es = body.face_edges(face_id);
    if (es.size() < 3) return Vector3D::UnitZ();
    const auto& e0 = body.edge(es[0]);
    const auto& e1 = body.edge(es[1]);
    // 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();
}

// Direction in face plane, perpendicular to an edge
Vector3D perpendicular_in_face(const Vector3D& edge_dir, const Vector3D& face_normal) {
    Vector3D d = face_normal.cross(edge_dir);
    double len = d.norm();
    if (len > 1e-12) return d / len;
    return Vector3D::UnitY();
}

// Create a fillet arc surface (ruled surface along fillet blend)
curves::NurbsSurface make_fillet_surface(
    const std::vector<Point3D>& t1_pts,   // tangent points on face 1
    const std::vector<Point3D>& t2_pts,   // tangent points on face 2
    const std::vector<Point3D>& mid_pts)  // arc midpoint (for curvature)
{
    int n = static_cast<int>(t1_pts.size());
    if (n < 2) n = 2;
    // Build control grid: (n, 3) — n samples along edge, 3 across fillet
    std::vector<std::vector<Point3D>> grid(n);
    for (int i = 0; i < n; ++i) {
        grid[i] = {t1_pts[i], mid_pts[i], t2_pts[i]};
    }
    std::vector<double> knots_u(n + 3, 0.0);
    for (int i = 0; i < 3; ++i) knots_u[i] = 0.0;
    for (int i = n; i < n + 3; ++i) knots_u[i] = 1.0;
    for (int i = 3; i < n; ++i) knots_u[i] = static_cast<double>(i - 2) / (n - 2);
    // degree 2 in v for arc, degree 2 in u for smoothness along edge
    return curves::NurbsSurface(grid, knots_u, {0,0,0,1,1,1}, {}, 2, 2);
}

// Build vertices from a list of points into a model, return vertex IDs
std::vector<int> add_vertices(BrepModel& m, const std::vector<Point3D>& pts) {
    std::vector<int> ids;
    for (const auto& p : pts) ids.push_back(m.add_vertex(p));
    return ids;
}

// Build quad faces from vertex grid (n_rows × n_cols vertices), return face IDs
// Grid is stored row-major: grid[r][c]
std::vector<int> build_quad_faces(BrepModel& m,
    const std::vector<std::vector<int>>& grid)
{
    std::vector<int> face_ids;
    int rows = static_cast<int>(grid.size());
    int cols = static_cast<int>(grid[0].size());
    for (int r = 0; r + 1 < rows; ++r) {
        for (int c = 0; c + 1 < cols; ++c) {
            int v00 = grid[r][c], v10 = grid[r+1][c],
                v11 = grid[r+1][c+1], v01 = grid[r][c+1];
            int e1 = m.add_edge(v00, v10);
            int e2 = m.add_edge(v10, v11);
            int e3 = m.add_edge(v11, v01);
            int e4 = m.add_edge(v01, v00);
            Point3D p0 = m.vertex(v00).point, p1 = m.vertex(v10).point,
                    p2 = m.vertex(v11).point, p3 = m.vertex(v01).point;
            int s = m.add_surface(make_plane_surface(p0, p1, p2, p3));
            face_ids.push_back(m.add_face(s, {
                m.add_loop({e1, e2, e3, e4}, true)}));
        }
    }
    return face_ids;
}

}  // namespace

// ── Box/Cylinder/Sphere (same as before, reused) ──
BrepModel make_box(double w, double h, double d) {
    double hw=w/2, hh=h/2, hd=d/2;
    BrepModel model;
    Point3D corners[6][4] = {
        {{-hw,-hh,-hd},{hw,-hh,-hd},{hw,hh,-hd},{-hw,hh,-hd}},
        {{-hw,-hh,hd},{hw,-hh,hd},{hw,hh,hd},{-hw,hh,hd}},
        {{-hw,-hh,-hd},{hw,-hh,-hd},{hw,-hh,hd},{-hw,-hh,hd}},
        {{-hw,hh,-hd},{hw,hh,-hd},{hw,hh,hd},{-hw,hh,hd}},
        {{hw,-hh,-hd},{hw,hh,-hd},{hw,hh,hd},{hw,-hh,hd}},
        {{-hw,-hh,-hd},{-hw,hh,-hd},{-hw,hh,hd},{-hw,-hh,hd}},
    };
    // 8 unique corners (shared topology)
    int v000 = model.add_vertex({-hw,-hh,-hd}), v100 = model.add_vertex({ hw,-hh,-hd});
    int v110 = model.add_vertex({ hw, hh,-hd}), v010 = model.add_vertex({-hw, hh,-hd});
    int v001 = model.add_vertex({-hw,-hh, hd}), v101 = model.add_vertex({ hw,-hh, hd});
    int v111 = model.add_vertex({ hw, hh, hd}), v011 = model.add_vertex({-hw, hh, hd});
    int verts[6][4] = {
        {v000,v100,v110,v010},{v001,v101,v111,v011},{v000,v100,v101,v001},
        {v010,v110,v111,v011},{v100,v110,v111,v101},{v000,v010,v011,v001},
    };
    // Edge dedup helper
    auto find_edge = [&](int a, int b) -> int {
        for (size_t ei = 0; ei < model.num_edges(); ++ei) {
            auto& e = model.edge(static_cast<int>(ei));
            if ((e.v_start==a && e.v_end==b) || (e.v_start==b && e.v_end==a))
                return static_cast<int>(ei);
        }
        return model.add_edge(a, b);
    };
    std::vector<int> face_ids;
    for (int i = 0; i < 6; ++i) {
        int v0=verts[i][0],v1=verts[i][1],v2=verts[i][2],v3=verts[i][3];
        auto p0=model.vertex(v0).point,p1=model.vertex(v1).point;
        auto p2=model.vertex(v2).point,p3=model.vertex(v3).point;
        int s=model.add_surface(make_plane_surface(p0,p1,p2,p3));
        int e1=find_edge(v0,v1),e2=find_edge(v1,v2);
        int e3=find_edge(v2,v3),e4=find_edge(v3,v0);
        int lp=model.add_loop({e1,e2,e3,e4},true);
        face_ids.push_back(model.add_face(s,{lp}));
    }
    int sh=model.add_shell(face_ids,true);
    model.add_body({sh},"box");
    return model;
}

BrepModel make_cylinder(double r, double h, int segs) {
    BrepModel model;
    int vb=model.add_vertex({0,0,-h/2}), vt=model.add_vertex({0,0,h/2});
    std::vector<int> br(segs), tr(segs);
    for(int i=0;i<segs;++i){
        double a=2*M_PI*i/segs, x=r*cos(a), y=r*sin(a);
        br[i]=model.add_vertex({x,y,-h/2}); tr[i]=model.add_vertex({x,y,h/2});
    }
    std::vector<int> all_faces;
    // Bottom
    { std::vector<int> es; for(int i=0;i<segs;++i) es.push_back(model.add_edge(br[i],br[(i+1)%segs]));
      int s=model.add_surface(curves::NurbsSurface(
          {{{-r,-r,-h/2},{r,-r,-h/2}},{{-r,r,-h/2},{r,r,-h/2}}},{0,0,1,1},{0,0,1,1},{},1,1));
      all_faces.push_back(model.add_face(s,{model.add_loop(es,true)})); }
    // Top
    { std::vector<int> es; for(int i=0;i<segs;++i) es.push_back(model.add_edge(tr[(i+1)%segs],tr[i]));
      int s=model.add_surface(curves::NurbsSurface(
          {{{-r,-r,h/2},{r,-r,h/2}},{{-r,r,h/2},{r,r,h/2}}},{0,0,1,1},{0,0,1,1},{},1,1));
      all_faces.push_back(model.add_face(s,{model.add_loop(es,true)})); }
    // Sides
    for(int i=0;i<segs;++i){ int j=(i+1)%segs;
      int e1=model.add_edge(br[i],tr[i]), e2=model.add_edge(tr[i],tr[j]);
      int e3=model.add_edge(tr[j],br[j]), e4=model.add_edge(br[j],br[i]);
      double a0=2*M_PI*i/segs, a1=2*M_PI*j/segs;
      Point3D p0(r*cos(a0),r*sin(a0),-h/2), p1(r*cos(a1),r*sin(a1),-h/2);
      Point3D p2(r*cos(a1),r*sin(a1),h/2), p3(r*cos(a0),r*sin(a0),h/2);
      int s=model.add_surface(make_plane_surface(p0,p1,p2,p3));
      all_faces.push_back(model.add_face(s,{model.add_loop({e1,e2,e3,e4},true)})); }
    int sh=model.add_shell(all_faces,true);
    model.add_body({sh},"cylinder");
    return model;
}

BrepModel make_sphere(double r, int su, int sv) {
    BrepModel model;
    int vt=model.add_vertex({0,0,r}), vb=model.add_vertex({0,0,-r});
    std::vector<std::vector<int>> rings(sv+1);
    for(int j=1;j<sv;++j){
        double phi=M_PI*j/sv, z=r*cos(phi), rr=r*sin(phi);
        for(int i=0;i<su;++i){
            double t=2*M_PI*i/su;
            rings[j].push_back(model.add_vertex({rr*cos(t),rr*sin(t),z}));
        }
    }
    std::vector<int> all_faces;
    for(int i=0;i<su;++i){ int j=(i+1)%su;
      int e1=model.add_edge(vt,rings[1][i]), e2=model.add_edge(rings[1][i],rings[1][j]), e3=model.add_edge(rings[1][j],vt);
      auto g=std::vector<std::vector<Point3D>>{{{0,0,r},model.vertex(rings[1][j]).point},
          {{0,0,r},model.vertex(rings[1][i]).point}};
      int s=model.add_surface(curves::NurbsSurface(g,{0,0,1,1},{0,0,1,1},{},1,1));
      all_faces.push_back(model.add_face(s,{model.add_loop({e1,e2,e3},true)})); }
    for(int row=1;row<sv-1;++row) for(int i=0;i<su;++i){ int j=(i+1)%su;
      int v0=rings[row][i], v1=rings[row][j], v2=rings[row+1][j], v3=rings[row+1][i];
      int e1=model.add_edge(v0,v1), e2=model.add_edge(v1,v2), e3=model.add_edge(v2,v3), e4=model.add_edge(v3,v0);
      Point3D p0=model.vertex(v0).point, p1=model.vertex(v1).point, p2=model.vertex(v2).point, p3=model.vertex(v3).point;
      int s=model.add_surface(make_plane_surface(p0,p1,p2,p3));
      all_faces.push_back(model.add_face(s,{model.add_loop({e1,e2,e3,e4},true)})); }
    for(int i=0;i<su;++i){ int j=(i+1)%su;
      int e1=model.add_edge(rings[sv-1][i],vb), e2=model.add_edge(vb,rings[sv-1][j]), e3=model.add_edge(rings[sv-1][j],rings[sv-1][i]);
      auto g=std::vector<std::vector<Point3D>>{{{0,0,-r},model.vertex(rings[sv-1][j]).point},
          {{0,0,-r},model.vertex(rings[sv-1][i]).point}};
      int s=model.add_surface(curves::NurbsSurface(g,{0,0,1,1},{0,0,1,1},{},1,1));
      all_faces.push_back(model.add_face(s,{model.add_loop({e1,e2,e3},true)})); }
    int sh=model.add_shell(all_faces,true);
    model.add_body({sh},"sphere");
    return model;
}

// ── Extrude/Revolve/Sweep/Loft ──
BrepModel extrude(const curves::NurbsCurve& profile, const Vector3D& dir) {
    BrepModel result; const auto& cp=profile.control_points(); int n=cp.size();
    std::vector<int> bot,top;
    for(int i=0;i<n;++i){ bot.push_back(result.add_vertex(cp[i])); top.push_back(result.add_vertex(cp[i]+dir)); }
    std::vector<int> bot_es,top_es,all_faces;
    for(int i=0;i<n;++i){ int j=(i+1)%n; bot_es.push_back(result.add_edge(bot[i],bot[j])); }
    top_es.push_back(result.add_edge(top[n-1],top[0]));
    for(int i=1;i<n;++i) top_es.push_back(result.add_edge(top[i-1],top[i]));
    // Bottom face
    { Point3D bmin=cp[0],bmax=cp[0]; for(auto&p:cp){bmin=bmin.cwiseMin(p);bmax=bmax.cwiseMax(p);}
      curves::NurbsSurface bs({{{bmin.x(),bmin.y(),bmin.z()},{bmax.x(),bmin.y(),bmin.z()}},
                               {{bmin.x(),bmax.y(),bmin.z()},{bmax.x(),bmax.y(),bmin.z()}}},
                             {0,0,1,1},{0,0,1,1},{},1,1);
      int si=result.add_surface(bs); all_faces.push_back(result.add_face(si,{result.add_loop(bot_es,true)})); }
    // Top face
    { Point3D bmin=cp[0]+dir,bmax=cp[0]+dir; for(auto&p:cp){bmin=bmin.cwiseMin(p+dir);bmax=bmax.cwiseMax(p+dir);}
      curves::NurbsSurface ts({{{bmin.x(),bmin.y(),bmin.z()},{bmax.x(),bmin.y(),bmin.z()}},
                               {{bmin.x(),bmax.y(),bmin.z()},{bmax.x(),bmax.y(),bmin.z()}}},
                             {0,0,1,1},{0,0,1,1},{},1,1);
      int si=result.add_surface(ts); all_faces.push_back(result.add_face(si,{result.add_loop(top_es,true)})); }
    // Side faces
    for(int i=0;i<n;++i){ int j=(i+1)%n;
      int e1=result.add_edge(bot[i],top[i]), e2=result.add_edge(top[i],top[j]);
      int e3=result.add_edge(top[j],bot[j]), e4=result.add_edge(bot[j],bot[i]);
      int si=result.add_surface(make_plane_surface(cp[i],cp[j],cp[j]+dir,cp[i]+dir));
      all_faces.push_back(result.add_face(si,{result.add_loop({e1,e2,e3,e4},true)})); }
    int sh=result.add_shell(all_faces,true); result.add_body({sh},"extrude"); return result;
}

BrepModel revolve(const curves::NurbsCurve& profile, const Point3D& o, const Vector3D& ax, double ang) {
    BrepModel result; int n=profile.control_points().size(), samples=32;
    Vector3D a=ax.normalized();
    std::vector<std::vector<int>> rings(samples+1);
    for(int i=0;i<=samples;++i){
        double t=ang*i/samples; Eigen::AngleAxis<double> rot(t,a);
        for(int j=0;j<n;++j){ Vector3D rel=profile.control_points()[j]-o; rings[i].push_back(result.add_vertex(o+rot*rel)); }
    }
    std::vector<int> all_faces;
    for(int i=0;i<samples;++i) for(int j=0;j<n-1;++j){
        int v00=rings[i][j],v10=rings[i+1][j],v11=rings[i+1][j+1],v01=rings[i][j+1];
        int e1=result.add_edge(v00,v10),e2=result.add_edge(v10,v11),e3=result.add_edge(v11,v01),e4=result.add_edge(v01,v00);
        Point3D p0=result.vertex(v00).point,p1=result.vertex(v10).point,p2=result.vertex(v11).point,p3=result.vertex(v01).point;
        int s=result.add_surface(make_plane_surface(p0,p1,p2,p3));
        all_faces.push_back(result.add_face(s,{result.add_loop({e1,e2,e3,e4},true)})); }
    int sh=result.add_shell(all_faces,true); result.add_body({sh},"revolve"); return result;
}

BrepModel sweep(const curves::NurbsCurve& profile, const curves::NurbsCurve& path) {
    BrepModel result; int samples=32, pn=profile.control_points().size();
    auto [t0,t1]=path.domain();
    std::vector<Point3D> pos; std::vector<Vector3D> tans;
    for(int i=0;i<=samples;++i){ double t=t0+(t1-t0)*i/samples; pos.push_back(path.evaluate(t));
        auto dt=path.derivative(t,1); tans.push_back(dt.norm()>1e-10?dt.normalized():Vector3D::UnitZ()); }
    std::vector<std::vector<int>> rings(samples+1);
    for(int i=0;i<=samples;++i){
        Vector3D T=tans[i], N=std::abs(T.x())<0.9?Vector3D::UnitX().cross(T).normalized():Vector3D::UnitY().cross(T).normalized(), B=T.cross(N);
        for(int j=0;j<pn;++j){ const auto& cp=profile.control_points()[j]; rings[i].push_back(result.add_vertex(pos[i]+cp.x()*N+cp.y()*B+cp.z()*T)); }
    }
    std::vector<int> all_faces;
    for(int i=0;i<samples;++i) for(int j=0;j<pn-1;++j){
        int v00=rings[i][j],v10=rings[i+1][j],v11=rings[i+1][j+1],v01=rings[i][j+1];
        int e1=result.add_edge(v00,v10),e2=result.add_edge(v10,v11),e3=result.add_edge(v11,v01),e4=result.add_edge(v01,v00);
        Point3D p0=result.vertex(v00).point,p1=result.vertex(v10).point,p2=result.vertex(v11).point,p3=result.vertex(v01).point;
        int s=result.add_surface(make_plane_surface(p0,p1,p2,p3));
        all_faces.push_back(result.add_face(s,{result.add_loop({e1,e2,e3,e4},true)})); }
    int sh=result.add_shell(all_faces,true); result.add_body({sh},"sweep"); return result;
}

BrepModel loft(const std::vector<curves::NurbsCurve>& profiles) {
    if(profiles.size()<2) return BrepModel();
    BrepModel result; int pn=profiles[0].control_points().size();
    std::vector<std::vector<int>> rings;
    for(auto&p:profiles){ std::vector<int> ring; for(auto&cp:p.control_points()) ring.push_back(result.add_vertex(cp)); rings.push_back(ring); }
    std::vector<int> all_faces;
    for(size_t i=0;i+1<rings.size();++i) for(int j=0;j<pn;++j){ int k=(j+1)%pn;
        int v00=rings[i][j],v10=rings[i+1][j],v11=rings[i+1][k],v01=rings[i][k];
        int e1=result.add_edge(v00,v10),e2=result.add_edge(v10,v11),e3=result.add_edge(v11,v01),e4=result.add_edge(v01,v00);
        Point3D p0=result.vertex(v00).point,p1=result.vertex(v10).point,p2=result.vertex(v11).point,p3=result.vertex(v01).point;
        int s=result.add_surface(make_plane_surface(p0,p1,p2,p3));
        all_faces.push_back(result.add_face(s,{result.add_loop({e1,e2,e3,e4},true)})); }
    int sh=result.add_shell(all_faces,true); result.add_body({sh},"loft"); return result;
}

// ── Fillet ──
BrepModel fillet(const BrepModel& body, int edge_id, double radius) {
    // Validation
    int num_edges = static_cast<int>(body.num_edges());
    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.empty()) return body;
    int face_a = efs[0];

    // 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

    const auto& curve = *edge.curve;
    // Compute face normals
    Vector3D n_a = compute_face_normal(body, face_a);
    Vector3D n_b = compute_face_normal(body, face_b);

    // Check if faces are nearly coplanar (no meaningful fillet)
    double cos_angle = n_a.dot(n_b);
    if (std::abs(cos_angle) > 0.9999) return body;

    // Bisector direction
    Vector3D bisector = (n_a + n_b).normalized();
    double cos_half = bisector.dot(n_a);
    if (std::abs(cos_half) < 1e-8) return body;

    // Sample the edge
    int N = 8;  // Number of samples along the edge
    auto [t_min, t_max] = curve.domain();

    std::vector<Point3D> edge_pts(N + 1);
    std::vector<Point3D> t1_pts(N + 1);   // Tangent points on face_a
    std::vector<Point3D> t2_pts(N + 1);   // Tangent points on face_b
    std::vector<Point3D> mid_pts(N + 1);  // Fillet arc midpoints

    for (int i = 0; i <= N; ++i) {
        double t = t_min + (t_max - t_min) * static_cast<double>(i) / N;
        edge_pts[i] = curve.evaluate(t);

        // Fillet ball center
        double offset = radius / cos_half;
        Point3D C = edge_pts[i] + offset * bisector;

        // Tangent points
        t1_pts[i] = C - radius * n_a;
        t2_pts[i] = C - radius * n_b;

        // Arc midpoint (for better surface approximation)
        // Midpoint of the fillet arc: along the bisector direction from center
        Vector3D to_mid = (t1_pts[i] + t2_pts[i]) * 0.5 - C;
        double mid_len = to_mid.norm();
        if (mid_len > 1e-10) {
            mid_pts[i] = C + to_mid.normalized() * radius;
        } else {
            mid_pts[i] = (t1_pts[i] + t2_pts[i]) * 0.5;
        }
    }

    BrepModel result;

    // ── 1. Copy non-affected faces verbatim ──
    // For each face NOT adjacent to the filleted edge, copy all its geometry
    for (size_t fi = 0; fi < body.num_faces(); ++fi) {
        if (static_cast<int>(fi) == face_a || static_cast<int>(fi) == face_b)
            continue;

        // Copy this face into result
        const auto& f_src = body.face(static_cast<int>(fi));
        std::vector<int> new_loop_ids;
        for (int li : f_src.loops) {
            // Find the loop
            const auto& lop = body.loop_by_id(li);
            std::vector<int> new_es;
            for (int ei : lop.edges) {
                const auto& e_src = body.edge(ei);
                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));
            }
            new_loop_ids.push_back(result.add_loop(new_es, lop.is_outer));
        }
        // Copy surface
        int new_surf = result.add_surface(body.surface(f_src.surface_id));
        result.add_face(new_surf, new_loop_ids);
    }

    // ── 2. Rebuild adjacent faces with tangent offset ──
    auto rebuild_face = [&](int face_id, const Vector3D& n,
                            const std::vector<Point3D>& tangent_pts) {
        const auto& f_src = body.face(face_id);
        const auto& curve_copy = body.surface(f_src.surface_id);
        const auto& edge_src = body.edge(edge_id);

        // Get the face's edges in order
        auto f_edges = body.face_edges(face_id);

        // Find the index of the filleted edge in this face's edge list
        int fillet_edge_idx = -1;
        for (size_t ei = 0; ei < f_edges.size(); ++ei) {
            if (f_edges[ei] == edge_id) {
                fillet_edge_idx = static_cast<int>(ei);
                break;
            }
        }
        if (fillet_edge_idx < 0) return;  // Edge not found in this face — shouldn't happen

        // Rebuild the face: replace the filleted edge with the tangent offset edge
        std::vector<int> new_es;

        // We need to know the direction of the tangent edge relative
        // to the original edge orientation in this face's loop
        for (const auto& lop : body.all_loops()) {
            bool found = false;
            for (int li : f_src.loops) {
                if (lop.id == li) { found = true; break; }
            }
            if (!found) continue;

            // Determine the orientation of the edge in this loop
            for (size_t ei = 0; ei < lop.edges.size(); ++ei) {
                int eidx = lop.edges[ei];
                if (eidx != edge_id) continue;

                // Found the edge — now we can determine its orientation
                bool edge_reversed = (edge_src.v_start != body.edge(eidx).v_start);

                // Collect all edges in loop order
                std::vector<int> sorted_edges(lop.edges.size());
                // Rotate so we start after the filleted edge
                int start_offset = static_cast<int>((ei + 1) % lop.edges.size());
                for (size_t j = 0; j < lop.edges.size(); ++j) {
                    sorted_edges[j] = lop.edges[(start_offset + static_cast<int>(j)) % lop.edges.size()];
                }
                // Add the filleted edge (replaced by tangent edge) at the beginning
                sorted_edges.insert(sorted_edges.begin(), edge_id);

                // Edge direction for the tangent edge
                Point3D t0 = edge_reversed ? tangent_pts.back() : tangent_pts.front();
                Point3D t1 = edge_reversed ? tangent_pts.front() : tangent_pts.back();
                int nvt0 = result.add_vertex(t0);
                int nvt1 = result.add_vertex(t1);
                new_es.push_back(result.add_edge(nvt0, nvt1));

                // Add the other edges (these are the edges NOT being filleted)
                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_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 ||
                                       e_src.v_start == edge_src.v_end);
                    bool v1_on_edge = (e_src.v_end == edge_src.v_start ||
                                       e_src.v_end == edge_src.v_end);

                    // Determine which tangent point corresponds to which vertex
                    Point3D q0 = p0, q1 = p1;
                    if (v0_on_edge) {
                        bool is_start = (e_src.v_start == edge_src.v_start);
                        double dist0 = (is_start ? tangent_pts.front() : tangent_pts.back())
                                           .dot(n) - p0.dot(n);
                        q0 = p0 + n * dist0;
                        // Actually, use the tangent point directly if the vertex
                        // exactly matches the edge endpoint
                        q0 = (is_start ? tangent_pts.front() : tangent_pts.back());
                    }
                    if (v1_on_edge) {
                        bool is_start = (e_src.v_end == edge_src.v_start);
                        q1 = (is_start ? tangent_pts.front() : tangent_pts.back());
                    }

                    int nv0 = result.add_vertex(q0);
                    int nv1 = result.add_vertex(q1);
                    new_es.push_back(result.add_edge(nv0, nv1));
                }
                break;
            }
            break;
        }

        // Build the rebuilt face
        if (!new_es.empty()) {
            int new_loop = result.add_loop(new_es, true);
            int new_surf = result.add_surface(curve_copy);
            result.add_face(new_surf, {new_loop});
        }
    };

    rebuild_face(face_a, n_a, t1_pts);
    rebuild_face(face_b, n_b, t2_pts);

    // ── 3. Build fillet surface faces ──
    for (int i = 0; i < N; ++i) {
        // Create a quad face for the fillet strip between samples i and i+1
        Point3D q00 = t1_pts[i],     q01 = t2_pts[i];
        Point3D q10 = t1_pts[i + 1], q11 = t2_pts[i + 1];

        // Build small fillet surface patch (degree 2 × 2 for curvature)
        std::vector<std::vector<Point3D>> grid = {
            {q00, mid_pts[i],     q01},
            {q10, mid_pts[i + 1], q11}
        };
        curves::NurbsSurface fs(grid,
            {0, 0, 1, 1},
            {0, 0, 0, 1, 1, 1},
            {}, 1, 2);

        int v0 = result.add_vertex(q00), v1 = result.add_vertex(q10);
        int v2 = result.add_vertex(q11), v3 = result.add_vertex(q01);
        int e1 = result.add_edge(v0, v1);
        int e2 = result.add_edge(v1, v2);
        int e3 = result.add_edge(v2, v3);
        int e4 = result.add_edge(v3, v0);
        int surf_id = result.add_surface(fs);
        int loop_id = result.add_loop({e1, e2, e3, e4}, true);
        result.add_face(surf_id, {loop_id});
    }

    // ── 4. Build shell and body ──
    std::vector<int> all_face_ids;
    for (size_t fi = 0; fi < result.num_faces(); ++fi) {
        all_face_ids.push_back(result.face(static_cast<int>(fi)).id);
    }
    if (!all_face_ids.empty()) {
        int sh = result.add_shell(all_face_ids, true);
        result.add_body({sh}, "fillet");
    }
    return result;
}

// ── Variable-radius fillet ──
BrepModel fillet_variable(const BrepModel& body, int edge_id,
                          double r_start, double r_end, int samples) {
    int num_edges = static_cast<int>(body.num_edges());
    if (edge_id < 0 || edge_id >= num_edges) return body;
    if (r_start < 0.0 || r_end < 0.0) return body;
    if (r_start == 0.0 && r_end == 0.0) return body;

    auto efs = body.edge_faces(edge_id);
    if (efs.empty()) return body;
    int face_a = efs[0];

    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;

    const auto& curve = *edge.curve;
    Vector3D n_a = compute_face_normal(body, face_a);
    Vector3D n_b = compute_face_normal(body, face_b);

    double cos_angle = n_a.dot(n_b);
    if (std::abs(cos_angle) > 0.9999) return body;

    Vector3D bisector = (n_a + n_b).normalized();
    double cos_half = bisector.dot(n_a);
    if (std::abs(cos_half) < 1e-8) return body;

    auto [t_min, t_max] = curve.domain();

    int N = samples;
    std::vector<double> radii(N + 1);
    std::vector<Point3D> edge_pts(N + 1);
    std::vector<Point3D> t1_pts(N + 1);
    std::vector<Point3D> t2_pts(N + 1);
    std::vector<Point3D> mid_pts(N + 1);

    for (int i = 0; i <= N; ++i) {
        double t_val = static_cast<double>(i) / N;
        double r = r_start + (r_end - r_start) * t_val;
        double t = t_min + (t_max - t_min) * t_val;
        radii[i] = r;
        edge_pts[i] = curve.evaluate(t);

        double offset = r / cos_half;
        Point3D C = edge_pts[i] + offset * bisector;

        t1_pts[i] = C - r * n_a;
        t2_pts[i] = C - r * n_b;

        Vector3D to_mid = (t1_pts[i] + t2_pts[i]) * 0.5 - C;
        double mid_len = to_mid.norm();
        if (mid_len > 1e-10) {
            mid_pts[i] = C + to_mid.normalized() * r;
        } else {
            mid_pts[i] = (t1_pts[i] + t2_pts[i]) * 0.5;
        }
    }

    BrepModel result;

    // ── 1. Copy non-affected faces verbatim ──
    for (size_t fi = 0; fi < body.num_faces(); ++fi) {
        if (static_cast<int>(fi) == face_a || static_cast<int>(fi) == face_b)
            continue;
        const auto& f_src = body.face(static_cast<int>(fi));
        std::vector<int> new_loop_ids;
        for (int li : f_src.loops) {
            const auto& lop = body.loop_by_id(li);
            std::vector<int> new_es;
            for (int ei : lop.edges) {
                const auto& e_src = body.edge(ei);
                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));
            }
            new_loop_ids.push_back(result.add_loop(new_es, lop.is_outer));
        }
        int new_surf = result.add_surface(body.surface(f_src.surface_id));
        result.add_face(new_surf, new_loop_ids);
    }

    // ── 2. Rebuild adjacent faces with variable tangent offset ──
    auto rebuild_face = [&](int face_id, const Vector3D& n,
                            const std::vector<Point3D>& tangent_pts) {
        const auto& f_src = body.face(face_id);
        const auto& curve_copy = body.surface(f_src.surface_id);
        const auto& edge_src = body.edge(edge_id);

        auto f_edges = body.face_edges(face_id);
        int fillet_edge_idx = -1;
        for (size_t ei = 0; ei < f_edges.size(); ++ei) {
            if (f_edges[ei] == edge_id) {
                fillet_edge_idx = static_cast<int>(ei);
                break;
            }
        }
        if (fillet_edge_idx < 0) return;

        std::vector<int> new_es;
        for (const auto& lop : body.all_loops()) {
            bool found = false;
            for (int li : f_src.loops) {
                if (lop.id == li) { found = true; break; }
            }
            if (!found) continue;

            for (size_t ei = 0; ei < lop.edges.size(); ++ei) {
                int eidx = lop.edges[ei];
                if (eidx != edge_id) continue;

                bool edge_reversed = (edge_src.v_start != body.edge(eidx).v_start);

                std::vector<int> sorted_edges(lop.edges.size());
                int start_offset = static_cast<int>((ei + 1) % lop.edges.size());
                for (size_t j = 0; j < lop.edges.size(); ++j) {
                    sorted_edges[j] = lop.edges[(start_offset + static_cast<int>(j)) % lop.edges.size()];
                }
                sorted_edges.insert(sorted_edges.begin(), edge_id);

                Point3D t0 = edge_reversed ? tangent_pts.back() : tangent_pts.front();
                Point3D t1 = edge_reversed ? tangent_pts.front() : tangent_pts.back();
                new_es.push_back(result.add_edge(result.add_vertex(t0),
                                                  result.add_vertex(t1)));

                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_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);
                    bool v1_on = (e_src.v_end == edge_src.v_start ||
                                  e_src.v_end == edge_src.v_end);

                    Point3D q0 = p0, q1 = p1;
                    if (v0_on) q0 = (e_src.v_start == edge_src.v_start) ?
                                   tangent_pts.front() : tangent_pts.back();
                    if (v1_on) q1 = (e_src.v_end == edge_src.v_start) ?
                                   tangent_pts.front() : tangent_pts.back();

                    new_es.push_back(result.add_edge(
                        result.add_vertex(q0), result.add_vertex(q1)));
                }
                break;
            }
            break;
        }

        if (!new_es.empty()) {
            int new_loop = result.add_loop(new_es, true);
            int new_surf = result.add_surface(curve_copy);
            result.add_face(new_surf, {new_loop});
        }
    };

    rebuild_face(face_a, n_a, t1_pts);
    rebuild_face(face_b, n_b, t2_pts);

    // ── 3. Build variable-radius fillet surface faces ──
    for (int i = 0; i < N; ++i) {
        Point3D q00 = t1_pts[i],     q01 = t2_pts[i];
        Point3D q10 = t1_pts[i + 1], q11 = t2_pts[i + 1];

        std::vector<std::vector<Point3D>> grid = {
            {q00, mid_pts[i],     q01},
            {q10, mid_pts[i + 1], q11}
        };
        curves::NurbsSurface fs(grid,
            {0, 0, 1, 1},
            {0, 0, 0, 1, 1, 1},
            {}, 1, 2);

        int v0 = result.add_vertex(q00), v1 = result.add_vertex(q10);
        int v2 = result.add_vertex(q11), v3 = result.add_vertex(q01);
        int e1 = result.add_edge(v0, v1);
        int e2 = result.add_edge(v1, v2);
        int e3 = result.add_edge(v2, v3);
        int e4 = result.add_edge(v3, v0);
        int surf_id = result.add_surface(fs);
        result.add_face(surf_id,
                        {result.add_loop({e1, e2, e3, e4}, true)});
    }

    // ── 4. Build shell and body ──
    std::vector<int> all_face_ids;
    for (size_t fi = 0; fi < result.num_faces(); ++fi) {
        all_face_ids.push_back(result.face(static_cast<int>(fi)).id);
    }
    if (!all_face_ids.empty()) {
        int sh = result.add_shell(all_face_ids, true);
        result.add_body({sh}, "fillet_variable");
    }
    return result;
}

// ── Chamfer ──
BrepModel chamfer(const BrepModel& body, int edge_id, double dist) {
    // Validation
    int num_edges = static_cast<int>(body.num_edges());
    if (edge_id < 0 || edge_id >= num_edges) return body;
    if (dist <= 0.0) return body;

    auto efs = body.edge_faces(edge_id);
    if (efs.empty()) return body;
    int face_a = efs[0];

    // Find face_b geometrically (edges are not shared between faces in this impl)
    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;

    const auto& curve = *edge.curve;

    Vector3D n_a = compute_face_normal(body, face_a);
    Vector3D n_b = compute_face_normal(body, face_b);

    // Edge direction
    auto [t_min, t_max] = curve.domain();
    Point3D p_start = curve.evaluate(t_min);
    Point3D p_end   = curve.evaluate(t_max);
    Vector3D edge_dir = (p_end - p_start);
    if (edge_dir.norm() < 1e-12) return body;
    edge_dir.normalize();

    // Offset direction on each face: perpendicular to edge, in face plane
    // The offset should go "away" from the edge into the face
    Vector3D dir_a = perpendicular_in_face(edge_dir, n_a);
    Vector3D dir_b = perpendicular_in_face(edge_dir, n_b);

    // Ensure both offset directions point away from the sharp corner
    // The fillet center should be between the two offset points
    // Check: offset points should be on opposite sides of the edge
    Point3D mid_a = p_start + dir_a * dist;
    Point3D mid_b = p_start + dir_b * dist;
    // If dir_a and dir_b point toward each other (concave), the chamfer doesn't make sense
    Vector3D gap = mid_b - mid_a;
    if (gap.norm() < dist * 0.1) {
        // Try flipping one direction
        dir_b = -dir_b;
        mid_b = p_start + dir_b * dist;
    }

    int N = 4;  // Samples along edge (fewer needed for flat chamfer)
    std::vector<Point3D> a_pts(N + 1), b_pts(N + 1);

    for (int i = 0; i <= N; ++i) {
        double t = t_min + (t_max - t_min) * static_cast<double>(i) / N;
        Point3D pt = curve.evaluate(t);
        a_pts[i] = pt + dir_a * dist;
        b_pts[i] = pt + dir_b * dist;
    }

    BrepModel result;

    // ── 1. Copy non-affected faces ──
    for (size_t fi = 0; fi < body.num_faces(); ++fi) {
        if (static_cast<int>(fi) == face_a || static_cast<int>(fi) == face_b)
            continue;
        const auto& f_src = body.face(static_cast<int>(fi));
        std::vector<int> new_loop_ids;
        for (int li : f_src.loops) {
            const auto& lop = body.loop_by_id(li);
            std::vector<int> new_es;
            for (int ei : lop.edges) {
                const auto& e_src = body.edge(ei);
                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)));
            }
            new_loop_ids.push_back(result.add_loop(new_es, lop.is_outer));
        }
        result.add_face(result.add_surface(body.surface(f_src.surface_id)),
                        new_loop_ids);
    }

    // ── 2. Rebuild adjacent faces ──
    auto rebuild_chamfer_face = [&](int face_id,
                                     const Vector3D& dir,
                                     const std::vector<Point3D>& offset_pts) {
        const auto& f_src = body.face(face_id);
        const auto& edge_src = body.edge(edge_id);

        for (const auto& lop : body.all_loops()) {
            bool found = false;
            for (int li : f_src.loops) if (lop.id == li) { found = true; break; }
            if (!found) continue;

            for (size_t ei = 0; ei < lop.edges.size(); ++ei) {
                if (lop.edges[ei] != edge_id) continue;

                bool edge_reversed = (edge_src.v_start != body.edge(lop.edges[ei]).v_start);
                Point3D t0 = edge_reversed ? offset_pts.back() : offset_pts.front();
                Point3D t1 = edge_reversed ? offset_pts.front() : offset_pts.back();

                std::vector<int> new_es;
                new_es.push_back(result.add_edge(
                    result.add_vertex(t0), result.add_vertex(t1)));

                // Other edges
                for (size_t j = 1; j < lop.edges.size(); ++j) {
                    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_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);
                    bool v1_on = (e_src.v_end == edge_src.v_start ||
                                  e_src.v_end == edge_src.v_end);

                    Point3D q0 = p0, q1 = p1;
                    if (v0_on) q0 = (e_src.v_start == edge_src.v_start) ?
                                   offset_pts.front() : offset_pts.back();
                    if (v1_on) q1 = (e_src.v_end == edge_src.v_start) ?
                                   offset_pts.front() : offset_pts.back();

                    new_es.push_back(result.add_edge(
                        result.add_vertex(q0), result.add_vertex(q1)));
                }

                if (!new_es.empty()) {
                    int new_loop = result.add_loop(new_es, true);
                    int new_surf = result.add_surface(
                        body.surface(f_src.surface_id));
                    result.add_face(new_surf, {new_loop});
                }
                break;
            }
            break;
        }
    };

    rebuild_chamfer_face(face_a, dir_a, a_pts);
    rebuild_chamfer_face(face_b, dir_b, b_pts);

    // ── 3. Build chamfer faces (flat faces connecting offsets) ──
    for (int i = 0; i < N; ++i) {
        Point3D q00 = a_pts[i],     q01 = b_pts[i];
        Point3D q10 = a_pts[i + 1], q11 = b_pts[i + 1];

        // Chamfer face: flat quad (make_plane_surface)
        int v0 = result.add_vertex(q00), v1 = result.add_vertex(q10);
        int v2 = result.add_vertex(q11), v3 = result.add_vertex(q01);
        int e1 = result.add_edge(v0, v1);
        int e2 = result.add_edge(v1, v2);
        int e3 = result.add_edge(v2, v3);
        int e4 = result.add_edge(v3, v0);
        int surf_id = result.add_surface(make_plane_surface(q00, q10, q11, q01));
        result.add_face(surf_id,
                        {result.add_loop({e1, e2, e3, e4}, true)});
    }

    // ── 4. Shell and body ──
    std::vector<int> all_face_ids;
    for (size_t fi = 0; fi < result.num_faces(); ++fi) {
        all_face_ids.push_back(result.face(static_cast<int>(fi)).id);
    }
    if (!all_face_ids.empty()) {
        int sh = result.add_shell(all_face_ids, true);
        result.add_body({sh}, "chamfer");
    }
    return result;
}

// ── Shell ──
BrepModel shell(const BrepModel& body, int open_face, double thickness) {
    BrepModel result;

    // ── 1. Offset all vertices inward ──
    std::map<int, int> outer_vid;  // vertex ID → new outer vertex ID
    std::map<int, int> inner_vid;  // vertex ID → new inner vertex ID

    for (size_t i = 0; i < body.num_vertices(); ++i) {
        const auto& v = body.vertex(static_cast<int>(i));
        int vid = v.id;
        Point3D p_outer = v.point;
        Point3D p_inner = offset_vertex(body, static_cast<int>(i), -thickness);
        outer_vid[vid] = result.add_vertex(p_outer);
        inner_vid[vid] = result.add_vertex(p_inner);
    }

    std::vector<int> all_faces;

    // ── 2. Create inner offset faces ──
    for (size_t fi = 0; fi < body.num_faces(); ++fi) {
        if (static_cast<int>(fi) == open_face) continue;  // Skip open face

        const auto& face = body.face(static_cast<int>(fi));
        const auto& surf = body.surface(face.surface_id);

        auto es = body.face_edges(static_cast<int>(fi));

        // Create inner face (offset inward)
        std::vector<int> inner_es;
        for (int ei : es) {
            const auto& e = body.edge(ei);
            int vs = inner_vid[e.v_start];
            int ve = inner_vid[e.v_end];
            inner_es.push_back(result.add_edge(vs, ve));
        }
        int inner_surf = result.add_surface(surf);
        int inner_loop = result.add_loop(inner_es, face.loops.empty() ? true : true);
        all_faces.push_back(result.add_face(inner_surf, {inner_loop}));

        // Create outer face (original position)
        std::vector<int> outer_es;
        for (int ei : es) {
            const auto& e = body.edge(ei);
            int vs = outer_vid[e.v_start];
            int ve = outer_vid[e.v_end];
            outer_es.push_back(result.add_edge(vs, ve));
        }
        int outer_surf = result.add_surface(surf);
        int outer_loop = result.add_loop(outer_es,
            face.loops.empty() ? true : true);
        all_faces.push_back(result.add_face(outer_surf, {outer_loop}));
    }

    // ── 3. Create side walls connecting outer to inner ──
    // For each edge of the open face, create a quad wall:
    //   outer_start → outer_end → inner_end → inner_start
    if (open_face >= 0 && open_face < static_cast<int>(body.num_faces())) {
        auto open_edges = body.face_edges(open_face);

        for (int ei : open_edges) {
            const auto& e = body.edge(ei);
            int ovs = outer_vid[e.v_start];
            int ove = outer_vid[e.v_end];
            int ivs = inner_vid[e.v_start];
            int ive = inner_vid[e.v_end];

            Point3D p_os = result.vertex(ovs).point;
            Point3D p_oe = result.vertex(ove).point;
            Point3D p_is = result.vertex(ivs).point;
            Point3D p_ie = result.vertex(ive).point;

            // Side wall quad: outer_start → outer_end → inner_end → inner_start
            int wall_e1 = result.add_edge(ovs, ove);
            int wall_e2 = result.add_edge(ove, ive);
            int wall_e3 = result.add_edge(ive, ivs);
            int wall_e4 = result.add_edge(ivs, ovs);

            int wall_surf = result.add_surface(
                make_plane_surface(p_os, p_oe, p_ie, p_is));
            int wall_loop = result.add_loop(
                {wall_e1, wall_e2, wall_e3, wall_e4}, true);
            all_faces.push_back(
                result.add_face(wall_surf, {wall_loop}));
        }
    }

    // ── 4. Build shell and body ──
    if (!all_faces.empty()) {
        int sh = result.add_shell(all_faces, true);
        result.add_body({sh}, "shell");
    }
    return result;
}

}  // namespace vde::brep