"""
NXOpen script — Extract sandbox face geometry for Adaptive Isogrid.

Runs from the .sim file context. Navigates:
    SIM → FEM → Idealized Part → find bodies with ISOGRID_SANDBOX attribute

For each sandbox body, exports `geometry_<sandbox_id>.json` containing:
    - outer_boundary: list of typed segments (line or arc) preserving exact geometry
    - inner_boundaries: same format for cutouts
    - transform: 3D <-> 2D mapping for reimporting geometry
    - thickness: from NX midsurface (if available)

Schema v2.0: segments are typed objects, not flat polylines.
  Line: {"type": "line", "start": [x,y], "end": [x,y]}
  Arc:  {"type": "arc", "start": [x,y], "end": [x,y], "center": [x,y], "radius": R, "clockwise": bool}

Inner loops are treated as boundary constraints (edges), NOT as holes to rib around,
because hole reservations are handled by separate solid cylinders in the fixed geometry.

Usage (NX Journal — just run it, no args needed):
    File > Execute > NX Journal > extract_sandbox.py

Author: Atomizer / Adaptive Isogrid
Created: 2026-02-16
"""

from __future__ import annotations

import json
import math
import os
import sys
from dataclasses import dataclass
from pathlib import Path
from typing import Any, Dict, List, Sequence, Tuple

# ---------------------------------------------------------------------------
# Geometry helpers (pure math, no NX dependency)
# ---------------------------------------------------------------------------

Point3D = Tuple[float, float, float]
Point2D = Tuple[float, float]


@dataclass
class LocalFrame:
    origin: Point3D
    x_axis: Point3D
    y_axis: Point3D
    normal: Point3D


def _norm(v: Sequence[float]) -> float:
    return math.sqrt(sum(c * c for c in v))


def _normalize(v: Sequence[float]) -> Tuple[float, float, float]:
    n = _norm(v)
    if n < 1e-12:
        return (0.0, 0.0, 1.0)
    return (v[0] / n, v[1] / n, v[2] / n)


def _dot(a: Sequence[float], b: Sequence[float]) -> float:
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]


def _cross(a: Sequence[float], b: Sequence[float]) -> Tuple[float, float, float]:
    return (
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0],
    )


def _sub(a: Sequence[float], b: Sequence[float]) -> Tuple[float, float, float]:
    return (a[0] - b[0], a[1] - b[1], a[2] - b[2])


def project_to_2d(points3d: Sequence[Point3D], frame: LocalFrame) -> List[Point2D]:
    out: List[Point2D] = []
    for p in points3d:
        v = _sub(p, frame.origin)
        out.append((_dot(v, frame.x_axis), _dot(v, frame.y_axis)))
    return out


def unproject_to_3d(points2d: Sequence[Point2D], frame: LocalFrame) -> List[Point3D]:
    """Inverse of project_to_2d — reconstruct 3D from local 2D coords."""
    out: List[Point3D] = []
    for x, y in points2d:
        px = frame.origin[0] + x * frame.x_axis[0] + y * frame.y_axis[0]
        py = frame.origin[1] + x * frame.x_axis[1] + y * frame.y_axis[1]
        pz = frame.origin[2] + x * frame.x_axis[2] + y * frame.y_axis[2]
        out.append((px, py, pz))
    return out


# ---------------------------------------------------------------------------
# Edge segment types
# ---------------------------------------------------------------------------

@dataclass
class EdgeSegment:
    """A typed geometry segment — either a line or an arc."""
    seg_type: str  # "line" or "arc"
    start_3d: Point3D
    end_3d: Point3D
    # Arc-specific (None for lines)
    center_3d: Point3D | None = None
    radius: float | None = None


# ---------------------------------------------------------------------------
# NX edge analysis — extract type + arc parameters
# ---------------------------------------------------------------------------

def _analyze_edge(edge: Any, lister: Any = None) -> EdgeSegment:
    """
    Analyze an NX edge and return a typed EdgeSegment.
    For arcs: extracts center + radius from UF.
    For lines: just start/end vertices.
    For unknown curves: falls back to line (vertices only) with a warning.
    """
    def _log(msg: str) -> None:
        if lister:
            lister.WriteLine(msg)

    # Get vertices
    try:
        v1, v2 = edge.GetVertices()
        p1 = (float(v1.X), float(v1.Y), float(v1.Z))
        p2 = (float(v2.X), float(v2.Y), float(v2.Z))
    except Exception as exc:
        raise RuntimeError(f"Edge.GetVertices() failed: {exc}")

    # Classify edge type
    edge_type_str = "?"
    try:
        edge_type_str = str(edge.SolidEdgeType)
    except Exception:
        pass

    is_linear = "Linear" in edge_type_str
    is_circular = "Circular" in edge_type_str

    # Linear edges — simple
    if is_linear:
        return EdgeSegment(seg_type="line", start_3d=p1, end_3d=p2)

    # Try UF Eval to detect arc
    try:
        import NXOpen
        import NXOpen.UF
        uf = NXOpen.UF.UFSession.GetUFSession()
        eval_obj = uf.Eval
        evaluator = None
        try:
            evaluator = eval_obj.Initialize2(edge.Tag)

            if eval_obj.IsArc(evaluator):
                arc_data = eval_obj.AskArc(evaluator)

                # Extract center and radius from arc_data
                center = None
                radius = None

                # Try named attributes (NXOpen struct)
                for attr in ('center', 'Center', 'arc_center'):
                    if hasattr(arc_data, attr):
                        c = getattr(arc_data, attr)
                        if hasattr(c, 'X'):
                            center = (float(c.X), float(c.Y), float(c.Z))
                        elif isinstance(c, (list, tuple)) and len(c) >= 3:
                            center = (float(c[0]), float(c[1]), float(c[2]))
                        break
                for attr in ('radius', 'Radius'):
                    if hasattr(arc_data, attr):
                        radius = float(getattr(arc_data, attr))
                        break

                # If named attrs didn't work, try UF Curve API
                if center is None or radius is None:
                    try:
                        curve_data = uf.Curve.AskArcData(edge.Tag)
                        if hasattr(curve_data, 'arc_center') and hasattr(curve_data, 'radius'):
                            c = curve_data.arc_center
                            center = (float(c[0]), float(c[1]), float(c[2]))
                            radius = float(curve_data.radius)
                    except Exception:
                        pass

                if center is not None and radius is not None and radius > 0.0:
                    _log(f"[edge] ARC: center=({center[0]:.3f},{center[1]:.3f},{center[2]:.3f}) "
                         f"r={radius:.3f}")
                    return EdgeSegment(
                        seg_type="arc", start_3d=p1, end_3d=p2,
                        center_3d=center, radius=radius,
                    )
                else:
                    _log(f"[edge] IsArc=True but could not extract center/radius. "
                         f"arc_data attrs: {[a for a in dir(arc_data) if not a.startswith('_')]}")
        finally:
            if evaluator is not None:
                try:
                    eval_obj.Free(evaluator)
                except Exception:
                    pass
    except Exception as exc:
        _log(f"[edge] UF arc detection failed: {exc}")

    # Fallback: try UF Curve.AskArcData directly (for circular edges not caught above)
    if is_circular:
        try:
            import NXOpen
            import NXOpen.UF
            uf = NXOpen.UF.UFSession.GetUFSession()
            curve_data = uf.Curve.AskArcData(edge.Tag)
            center = None
            radius = None
            if hasattr(curve_data, 'arc_center') and hasattr(curve_data, 'radius'):
                c = curve_data.arc_center
                center = (float(c[0]), float(c[1]), float(c[2]))
                radius = float(curve_data.radius)
            if center is not None and radius is not None and radius > 0.0:
                _log(f"[edge] ARC (UF Curve fallback): r={radius:.3f}")
                return EdgeSegment(
                    seg_type="arc", start_3d=p1, end_3d=p2,
                    center_3d=center, radius=radius,
                )
        except Exception as exc:
            _log(f"[edge] UF Curve.AskArcData fallback failed: {exc}")

    # Unknown curve type — warn and treat as line
    _log(f"[edge] WARNING: Non-line/arc edge type={edge_type_str}, treating as line (vertices only)")
    return EdgeSegment(seg_type="line", start_3d=p1, end_3d=p2)


# ---------------------------------------------------------------------------
# Face local frame
# ---------------------------------------------------------------------------

def _chain_edges_into_loops(
    edges: List[Any],
    lister: Any = None,
    tol: float = 0.01,
) -> List[Tuple[bool, List[EdgeSegment]]]:
    """
    Chain edges into closed loops by matching vertex endpoints.

    Returns list of (is_outer, segments) tuples where segments are EdgeSegment objects.
    The largest loop (by perimeter) is assumed to be the outer loop.
    """
    def _log(msg):
        if lister:
            lister.WriteLine(msg)

    if not edges:
        return []

    # Analyze each edge into a typed segment
    analyzed: List[Tuple[EdgeSegment, Any]] = []  # (segment, original_edge)
    for edge in edges:
        try:
            seg = _analyze_edge(edge, lister)
            analyzed.append((seg, edge))
        except Exception as exc:
            _log(f"[chain] Edge analysis failed: {exc}")
            continue

    _log(f"[chain] {len(analyzed)} edges analyzed ({sum(1 for s,_ in analyzed if s.seg_type == 'arc')} arcs, "
         f"{sum(1 for s,_ in analyzed if s.seg_type == 'line')} lines)")

    # Chain into loops
    used = [False] * len(analyzed)
    loops: List[List[EdgeSegment]] = []

    def pts_match(a: Point3D, b: Point3D) -> bool:
        return _norm(_sub(a, b)) < tol

    while True:
        start_idx = None
        for i, u in enumerate(used):
            if not u:
                start_idx = i
                break
        if start_idx is None:
            break

        chain: List[EdgeSegment] = []
        used[start_idx] = True
        seg, _ = analyzed[start_idx]
        chain.append(seg)

        current_end = seg.end_3d
        loop_start = seg.start_3d

        max_iters = len(analyzed) + 1
        for _ in range(max_iters):
            if pts_match(current_end, loop_start) and len(chain) > 1:
                break

            found = False
            for i, (s, _e) in enumerate(analyzed):
                if used[i]:
                    continue
                if pts_match(current_end, s.start_3d):
                    used[i] = True
                    chain.append(s)
                    current_end = s.end_3d
                    found = True
                    break
                elif pts_match(current_end, s.end_3d):
                    # Reversed edge — swap start/end
                    used[i] = True
                    reversed_seg = EdgeSegment(
                        seg_type=s.seg_type,
                        start_3d=s.end_3d,
                        end_3d=s.start_3d,
                        center_3d=s.center_3d,
                        radius=s.radius,
                    )
                    chain.append(reversed_seg)
                    current_end = s.start_3d
                    found = True
                    break

            if not found:
                _log(f"[chain] Warning: could not continue chain at "
                     f"({current_end[0]:.3f}, {current_end[1]:.3f}, {current_end[2]:.3f})")
                break

        loops.append(chain)

    _log(f"[chain] Built {len(loops)} loop(s)")

    if not loops:
        return []

    # Determine outer loop by perimeter
    def _loop_perimeter(segs: List[EdgeSegment]) -> float:
        total = 0.0
        for s in segs:
            if s.seg_type == "arc" and s.center_3d is not None and s.radius is not None:
                # Arc length = radius * angle
                r1 = _sub(s.start_3d, s.center_3d)
                r2 = _sub(s.end_3d, s.center_3d)
                cos_a = max(-1.0, min(1.0, _dot(_normalize(r1), _normalize(r2))))
                angle = math.acos(cos_a)
                total += s.radius * angle
            else:
                total += _norm(_sub(s.end_3d, s.start_3d))
        return total

    perimeters = [_loop_perimeter(segs) for segs in loops]
    outer_idx = perimeters.index(max(perimeters))

    result: List[Tuple[bool, List[EdgeSegment]]] = []
    for i, segs in enumerate(loops):
        is_outer = (i == outer_idx)
        n_arcs = sum(1 for s in segs if s.seg_type == "arc")
        n_lines = sum(1 for s in segs if s.seg_type == "line")
        result.append((is_outer, segs))
        _log(f"[chain]   loop {i}: {len(segs)} segments ({n_lines} lines, {n_arcs} arcs), "
             f"perimeter={perimeters[i]:.1f} mm {'(OUTER)' if is_outer else '(inner)'}")

    return result


def _face_local_frame(face: Any, lister: Any = None) -> LocalFrame:
    """
    Build a stable local frame on a planar face.
    """
    # Get a sample point from the first edge vertex
    edges = face.GetEdges()
    first_edge = edges[0]
    v1, v2 = first_edge.GetVertices()
    sample = (float(v1.X), float(v1.Y), float(v1.Z))

    # Get face normal
    normal = (0.0, 0.0, 1.0)
    try:
        import NXOpen
        pt = NXOpen.Point3d(sample[0], sample[1], sample[2])
        n = face.GetFaceNormal(pt)
        normal = _normalize((float(n.X), float(n.Y), float(n.Z)))
    except Exception:
        try:
            n = face.GetFaceNormal(sample[0], sample[1], sample[2])
            normal = _normalize((float(n.X), float(n.Y), float(n.Z)))
        except Exception:
            pass

    # Build orthonormal basis
    ref = (1.0, 0.0, 0.0) if abs(normal[0]) < 0.95 else (0.0, 1.0, 0.0)
    x_axis = _normalize(_cross(ref, normal))
    y_axis = _normalize(_cross(normal, x_axis))
    return LocalFrame(origin=sample, x_axis=x_axis, y_axis=y_axis, normal=normal)


# ---------------------------------------------------------------------------
# Attribute reading
# ---------------------------------------------------------------------------

def _get_string_attribute(obj: Any, title: str) -> str | None:
    """Try multiple NX API patterns to read a string attribute."""
    for method_name in ("GetStringUserAttribute", "GetUserAttributeAsString"):
        try:
            method = getattr(obj, method_name)
            val = method(title, -1)
            if val:
                return str(val)
        except Exception:
            continue
    return None


# ---------------------------------------------------------------------------
# SIM -> Idealized Part navigation
# ---------------------------------------------------------------------------

def _navigate_sim_to_idealized(session: Any) -> Any:
    """
    From the active .sim work part, navigate to the idealized part (_i.prt).
    Sets idealized part as work part and returns it.
    """
    work_part = session.Parts.Work
    part_name = work_part.Name if hasattr(work_part, "Name") else ""

    lister = session.ListingWindow
    lister.Open()
    lister.WriteLine(f"[extract_sandbox] Starting from: {part_name}")

    # Check if already in idealized part
    if part_name.endswith("_i"):
        lister.WriteLine("[extract_sandbox] Already in idealized part.")
        return work_part

    # Search loaded parts for the idealized part
    idealized_part = None
    for part in session.Parts:
        pname = part.Name if hasattr(part, "Name") else ""
        if pname.endswith("_i"):
            idealized_part = part
            lister.WriteLine(f"[extract_sandbox] Found idealized part: {pname}")
            break

    if idealized_part is None:
        raise RuntimeError(
            "Could not find idealized part (*_i.prt). "
            "Ensure the SIM is open with FEM + idealized part loaded."
        )

    # Set as work part
    try:
        session.Parts.SetWork(idealized_part)
        lister.WriteLine(f"[extract_sandbox] Set work part to: {idealized_part.Name}")
    except Exception as exc:
        lister.WriteLine(f"[extract_sandbox] Warning: SetWork failed: {exc}")

    return idealized_part


# ---------------------------------------------------------------------------
# Sandbox discovery
# ---------------------------------------------------------------------------

def find_sandbox_bodies(
    part: Any,
    lister: Any,
    attr_name: str = "ISOGRID_SANDBOX",
) -> List[Tuple[str, Any, Any]]:
    """
    Find bodies tagged with ISOGRID_SANDBOX attribute.

    Search order:
      1. Body-level attributes (part.Bodies)
      2. Face-level attributes
      3. Feature-level attributes (part history — Promote Body features)
      4. Feature name matching (e.g. 'Sandbox_1' in feature name)
      5. Body name matching

    Returns list of (sandbox_id, body, face) tuples.
    """
    tagged: List[Tuple[str, Any, Any]] = []
    found_ids: set = set()

    bodies = []
    try:
        bodies = list(part.Bodies.ToArray()) if hasattr(part.Bodies, "ToArray") else list(part.Bodies)
    except Exception:
        bodies = list(part.Bodies)

    lister.WriteLine(f"[extract_sandbox] Scanning {len(bodies)} bodies...")

    # --- Pass 1: body-level and face-level attributes ---
    for body in bodies:
        body_name = ""
        try:
            body_name = body.Name if hasattr(body, "Name") else str(body)
        except Exception:
            pass

        sandbox_id = _get_string_attribute(body, attr_name)
        if sandbox_id and sandbox_id not in found_ids:
            faces = body.GetFaces()
            if faces:
                tagged.append((sandbox_id, body, faces[0]))
                found_ids.add(sandbox_id)
                lister.WriteLine(f"[extract_sandbox] Found: {sandbox_id} (body attr on '{body_name}')")
                continue

        for face in body.GetFaces():
            sandbox_id = _get_string_attribute(face, attr_name)
            if sandbox_id and sandbox_id not in found_ids:
                tagged.append((sandbox_id, body, face))
                found_ids.add(sandbox_id)
                lister.WriteLine(f"[extract_sandbox] Found: {sandbox_id} (face attr on '{body_name}')")

    if tagged:
        return tagged

    # --- Pass 2: feature-level attributes (Promote Body features) ---
    lister.WriteLine("[extract_sandbox] No body/face attrs found, scanning features...")
    try:
        features = part.Features.ToArray() if hasattr(part.Features, "ToArray") else list(part.Features)
        lister.WriteLine(f"[extract_sandbox] Found {len(features)} features")

        for feat in features:
            feat_name = ""
            try:
                feat_name = feat.Name if hasattr(feat, "Name") else str(feat)
            except Exception:
                pass

            # Check feature attribute
            sandbox_id = _get_string_attribute(feat, attr_name)
            if sandbox_id and sandbox_id not in found_ids:
                # Get the body produced by this feature
                try:
                    feat_bodies = feat.GetBodies()
                    if feat_bodies:
                        body = feat_bodies[0]
                        faces = body.GetFaces()
                        if faces:
                            tagged.append((sandbox_id, body, faces[0]))
                            found_ids.add(sandbox_id)
                            lister.WriteLine(f"[extract_sandbox] Found: {sandbox_id} (feature attr on '{feat_name}')")
                except Exception as exc:
                    lister.WriteLine(f"[extract_sandbox] Feature '{feat_name}' has attr but GetBodies failed: {exc}")
    except Exception as exc:
        lister.WriteLine(f"[extract_sandbox] Feature scan error: {exc}")

    if tagged:
        return tagged

    # --- Pass 3: feature name matching (e.g. "Sandbox_1" in name) ---
    lister.WriteLine("[extract_sandbox] No feature attrs found, trying feature name matching...")
    try:
        features = part.Features.ToArray() if hasattr(part.Features, "ToArray") else list(part.Features)
        for feat in features:
            feat_name = ""
            try:
                feat_name = feat.Name if hasattr(feat, "Name") else str(feat)
            except Exception:
                continue

            if "sandbox" in feat_name.lower():
                try:
                    feat_bodies = feat.GetBodies()
                    if feat_bodies:
                        body = feat_bodies[0]
                        faces = body.GetFaces()
                        if faces:
                            sid = feat_name.lower().replace(" ", "_")
                            if sid not in found_ids:
                                tagged.append((sid, body, faces[0]))
                                found_ids.add(sid)
                                lister.WriteLine(f"[extract_sandbox] Found by feature name: {sid} ('{feat_name}')")
                except Exception as exc:
                    lister.WriteLine(f"[extract_sandbox] Feature '{feat_name}' name match but GetBodies failed: {exc}")
    except Exception:
        pass

    if tagged:
        return tagged

    # --- Pass 4: body name matching ---
    lister.WriteLine("[extract_sandbox] No features matched, trying body name matching...")
    for body in bodies:
        bname = ""
        try:
            bname = body.Name if hasattr(body, "Name") else str(body)
        except Exception:
            continue
        if "sandbox" in bname.lower():
            faces = body.GetFaces()
            if faces:
                sid = bname.lower().replace(" ", "_")
                if sid not in found_ids:
                    tagged.append((sid, body, faces[0]))
                    found_ids.add(sid)
                    lister.WriteLine(f"[extract_sandbox] Found by body name: {sid}")

    return tagged


# ---------------------------------------------------------------------------
# Core extraction
# ---------------------------------------------------------------------------

def _project_point_2d(pt3d: Point3D, frame: LocalFrame) -> Point2D:
    """Project a single 3D point to local 2D."""
    v = _sub(pt3d, frame.origin)
    return (_dot(v, frame.x_axis), _dot(v, frame.y_axis))


def _segments_to_json(segments: List[EdgeSegment], frame: LocalFrame) -> List[Dict[str, Any]]:
    """Convert a list of EdgeSegments to JSON-serializable dicts in 2D."""
    result = []
    for seg in segments:
        start_2d = _project_point_2d(seg.start_3d, frame)
        end_2d = _project_point_2d(seg.end_3d, frame)
        entry: Dict[str, Any] = {
            "type": seg.seg_type,
            "start": [round(start_2d[0], 6), round(start_2d[1], 6)],
            "end": [round(end_2d[0], 6), round(end_2d[1], 6)],
        }
        if seg.seg_type == "arc" and seg.center_3d is not None:
            center_2d = _project_point_2d(seg.center_3d, frame)
            entry["center"] = [round(center_2d[0], 6), round(center_2d[1], 6)]
            entry["radius"] = round(seg.radius, 6)
            # Determine clockwise/ccw: cross product of (start-center) × (end-center)
            # projected onto the face normal
            r1 = _sub(seg.start_3d, seg.center_3d)
            r2 = _sub(seg.end_3d, seg.center_3d)
            cross = _cross(r1, r2)
            dot_normal = _dot(cross, frame.normal)
            entry["clockwise"] = (dot_normal < 0)
        result.append(entry)
    return result


def extract_sandbox_geometry(
    face: Any,
    body: Any,
    sandbox_id: str,
    lister: Any,
) -> Dict[str, Any]:
    """
    Extract a sandbox face into a JSON-serializable dict.
    Schema v2.0: typed segments (line/arc) instead of polylines.
    Inner loops are boundary constraints (reserved geometry edges), not holes.
    """
    frame = _face_local_frame(face, lister)

    outer_segments: List[Dict[str, Any]] = []
    inner_boundaries: List[Dict[str, Any]] = []

    # Get all edges on the face and chain them into loops
    all_edges = list(face.GetEdges())
    lister.WriteLine(f"[extract_sandbox] {sandbox_id}: {len(all_edges)} edges on face")

    loops = _chain_edges_into_loops(all_edges, lister)
    lister.WriteLine(f"[extract_sandbox] {sandbox_id}: {len(loops)} loop(s) built")

    for loop_index, (is_outer, loop_segs) in enumerate(loops):
        seg_json = _segments_to_json(loop_segs, frame)
        n_arcs = sum(1 for s in seg_json if s["type"] == "arc")
        n_lines = sum(1 for s in seg_json if s["type"] == "line")

        if is_outer:
            outer_segments = seg_json
            lister.WriteLine(f"[extract_sandbox]   outer loop: {len(seg_json)} segments "
                           f"({n_lines} lines, {n_arcs} arcs)")
        else:
            inner_boundaries.append({
                "index": len(inner_boundaries),
                "segments": seg_json,
                "num_segments": len(seg_json),
            })
            lister.WriteLine(f"[extract_sandbox]   inner loop {len(inner_boundaries)}: "
                           f"{len(seg_json)} segments ({n_lines} lines, {n_arcs} arcs)")

    # Try thickness
    thickness = None
    try:
        thickness = float(body.GetThickness())
    except Exception:
        pass

    return {
        "schema_version": "2.0",
        "units": "mm",
        "sandbox_id": sandbox_id,
        "outer_boundary": outer_segments,
        "inner_boundaries": inner_boundaries,
        "num_inner_boundaries": len(inner_boundaries),
        "thickness": thickness,
        "transform": {
            "origin": [round(c, 6) for c in frame.origin],
            "x_axis": [round(c, 6) for c in frame.x_axis],
            "y_axis": [round(c, 6) for c in frame.y_axis],
            "normal": [round(c, 6) for c in frame.normal],
        },
    }


# ---------------------------------------------------------------------------
# Main — NX Journal entry point
# ---------------------------------------------------------------------------

def main():
    import NXOpen

    session = NXOpen.Session.GetSession()
    lister = session.ListingWindow
    lister.Open()

    lister.WriteLine("=" * 60)
    lister.WriteLine("  Adaptive Isogrid — Sandbox Geometry Extraction")
    lister.WriteLine("=" * 60)

    # Navigate to idealized part
    idealized_part = _navigate_sim_to_idealized(session)

    # Find sandboxes
    sandbox_entries = find_sandbox_bodies(idealized_part, lister)
    if not sandbox_entries:
        lister.WriteLine("[extract_sandbox] ERROR: No sandbox bodies found!")
        lister.WriteLine("Ensure bodies have ISOGRID_SANDBOX attribute set.")
        return

    lister.WriteLine(f"[extract_sandbox] Found {len(sandbox_entries)} sandbox(es)")

    # Output directory: next to the .sim file (or idealized part)
    try:
        part_dir = os.path.dirname(idealized_part.FullPath)
    except Exception:
        part_dir = os.getcwd()

    output_dir = os.path.join(part_dir, "adaptive_isogrid_data")
    os.makedirs(output_dir, exist_ok=True)
    lister.WriteLine(f"[extract_sandbox] Output dir: {output_dir}")

    # Extract each sandbox
    for sandbox_id, body, face in sandbox_entries:
        lister.WriteLine(f"\n--- Extracting {sandbox_id} ---")
        try:
            # Debug: print face info
            lister.WriteLine(f"[extract_sandbox] Face type: {type(face).__name__}")
            try:
                all_edges = face.GetEdges()
                lister.WriteLine(f"[extract_sandbox] Total edges on face: {len(all_edges)}")
            except Exception as exc:
                lister.WriteLine(f"[extract_sandbox] GetEdges failed: {exc}")

            geom = extract_sandbox_geometry(
                face=face,
                body=body,
                sandbox_id=sandbox_id,
                lister=lister,
            )

            out_path = os.path.join(output_dir, f"geometry_{sandbox_id}.json")
            with open(out_path, "w") as f:
                json.dump(geom, f, indent=2)
            lister.WriteLine(f"[extract_sandbox] Wrote: {out_path}")

            # Summary
            lister.WriteLine(f"  outer_boundary: {len(geom['outer_boundary'])} points")
            lister.WriteLine(f"  inner_boundaries: {geom['num_inner_boundaries']}")
            lister.WriteLine(f"  thickness: {geom['thickness']}")
        except Exception as exc:
            import traceback
            lister.WriteLine(f"[extract_sandbox] ERROR extracting {sandbox_id}: {exc}")
            lister.WriteLine(traceback.format_exc())

    lister.WriteLine("\n" + "=" * 60)
    lister.WriteLine(f"  Done — {len(sandbox_entries)} sandbox(es) exported")
    lister.WriteLine(f"  Output: {output_dir}")
    lister.WriteLine("=" * 60)


main()