feat(adaptive-isogrid): rewrite extract_sandbox.py - start from .sim, navigate to idealized part, find sandbox solid bodies by ISOGRID_SANDBOX attribute, inner loops as boundary constraints

tools/adaptive-isogrid/src/nx/extract_sandbox.py

@@ -1,11 +1,23 @@
 """
-NXOpen script — extract sandbox face geometry for Adaptive Isogrid.
+NXOpen script — Extract sandbox sheet geometry for Adaptive Isogrid.
 
-Finds faces tagged with user attribute:
+Runs from the .sim file context. Navigates:
-    ISOGRID_SANDBOX = sandbox_1, sandbox_2, ...
+    SIM → FEM → Idealized Part → find sheet bodies with ISOGRID_SANDBOX attribute
 
-For each sandbox face, exports `geometry_<sandbox_id>.json` in the same schema
+For each sandbox sheet body, exports `geometry_<sandbox_id>.json` containing:
-expected by the Python Brain (`outer_boundary`, `holes`, transform metadata, etc.).
+    - outer_boundary: 2D polyline of the sandbox outline
+    - inner_boundaries: 2D polylines of cutouts (reserved cylinder intersections, etc.)
+    - transform: 3D ↔ 2D mapping for reimporting geometry
+    - thickness: from NX midsurface (if available)
+
+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):
+    run_journal.exe extract_sandbox.py [--output-dir <path>] [--chord-tol <mm>]
+
+Author: Atomizer / Adaptive Isogrid
+Created: 2026-02-16
 """
 
 from __future__ import annotations
@@ -15,8 +27,11 @@ import json
 import math
 from dataclasses import dataclass
 from pathlib import Path
-from typing import Any, Dict, Iterable, List, Sequence, Tuple
+from typing import Any, Dict, List, Sequence, Tuple
 
+# ---------------------------------------------------------------------------
+# Geometry helpers (pure math, no NX dependency)
+# ---------------------------------------------------------------------------
 
 Point3D = Tuple[float, float, float]
 Point2D = Tuple[float, float]
@@ -57,60 +72,6 @@ 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 fit_circle(points: Sequence[Point2D]) -> Tuple[Point2D, float, float]:
-    """
-    Least-squares circle fit.
-    Returns (center, diameter, rms_error).
-    """
-    if len(points) < 3:
-        return ((0.0, 0.0), 0.0, float("inf"))
-
-    sx = sy = sxx = syy = sxy = 0.0
-    sxxx = syyy = sxxy = sxyy = 0.0
-
-    for x, y in points:
-        xx = x * x
-        yy = y * y
-        sx += x
-        sy += y
-        sxx += xx
-        syy += yy
-        sxy += x * y
-        sxxx += xx * x
-        syyy += yy * y
-        sxxy += xx * y
-        sxyy += x * yy
-
-    n = float(len(points))
-    c = n * sxx - sx * sx
-    d = n * sxy - sx * sy
-    e = n * (sxxx + sxyy) - (sxx + syy) * sx
-    g = n * syy - sy * sy
-    h = n * (sxxy + syyy) - (sxx + syy) * sy
-    denom = (c * g - d * d)
-
-    if abs(denom) < 1e-12:
-        return ((0.0, 0.0), 0.0, float("inf"))
-
-    a = (h * d - e * g) / denom
-    b = (h * c - e * d) / (d * d - g * c)
-    cx = -a / 2.0
-    cy = -b / 2.0
-    r = math.sqrt(max((a * a + b * b) / 4.0 - (sx * sx + sy * sy - n * (sxx + syy)) / n, 0.0))
-
-    errs = []
-    for x, y in points:
-        errs.append(abs(math.hypot(x - cx, y - cy) - r))
-    rms = math.sqrt(sum(e * e for e in errs) / len(errs)) if errs else float("inf")
-
-    return ((cx, cy), 2.0 * r, rms)
-
-
-def is_loop_circular(points2d: Sequence[Point2D], tol_mm: float = 0.5) -> Tuple[bool, Point2D, float]:
-    center, dia, rms = fit_circle(points2d)
-    return (rms <= tol_mm, center, dia)
-
-
 def project_to_2d(points3d: Sequence[Point3D], frame: LocalFrame) -> List[Point2D]:
     out: List[Point2D] = []
     for p in points3d:
@@ -119,20 +80,27 @@ def project_to_2d(points3d: Sequence[Point3D], frame: LocalFrame) -> List[Point2
     return out
 
 
-def _close_polyline(points: List[Point3D]) -> List[Point3D]:
+def unproject_to_3d(points2d: Sequence[Point2D], frame: LocalFrame) -> List[Point3D]:
-    if not points:
+    """Inverse of project_to_2d — reconstruct 3D from local 2D coords."""
-        return points
+    out: List[Point3D] = []
-    if _norm(_sub(points[0], points[-1])) > 1e-6:
+    for x, y in points2d:
-        points.append(points[0])
+        px = frame.origin[0] + x * frame.x_axis[0] + y * frame.y_axis[0]
-    return points
+        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
 
 
+# ---------------------------------------------------------------------------
+# NX edge sampling
+# ---------------------------------------------------------------------------
+
 def _sample_edge_polyline(edge: Any, chord_tol_mm: float) -> List[Point3D]:
     """
-    Sample an NX edge as a polyline.
+    Sample an NX edge as a polyline with adaptive point density.
-    NOTE: NX APIs vary by curve type; this helper intentionally keeps fallback logic.
+    Falls back to vertex extraction if evaluator is unavailable.
     """
-    # Preferred path: use evaluator where available.
+    # Preferred: parametric evaluator
     try:
         evaluator = edge.CreateEvaluator()
         t0, t1 = evaluator.GetLimits()
@@ -147,80 +115,253 @@ def _sample_edge_polyline(edge: Any, chord_tol_mm: float) -> List[Point3D]:
     except Exception:
         pass
 
-    # Fallback: edge vertices only (less accurate, but safe fallback).
+    # Fallback: edge vertices only
     try:
         verts = edge.GetVertices()
-        pts = []
+        return [(float(v.Coordinates.X), float(v.Coordinates.Y), float(v.Coordinates.Z)) for v in verts]
-        for v in verts:
-            p = v.Coordinates
-            pts.append((float(p.X), float(p.Y), float(p.Z)))
-        return pts
     except Exception as exc:
         raise RuntimeError(f"Could not sample edge polyline: {exc}")
 
 
+def _close_polyline(points: List[Point3D]) -> List[Point3D]:
+    if not points:
+        return points
+    if _norm(_sub(points[0], points[-1])) > 1e-6:
+        points.append(points[0])
+    return points
+
+
+# ---------------------------------------------------------------------------
+# Face local frame
+# ---------------------------------------------------------------------------
+
 def _face_local_frame(face: Any) -> LocalFrame:
     """
-    Build a stable local frame on a face:
+    Build a stable local frame on a planar face:
     - origin: first loop first point
-    - normal: face normal near origin
+    - normal: face normal at origin
-    - x/y axes: orthonormal basis on tangent plane
+    - x/y: orthonormal tangent basis
     """
     loops = face.GetLoops()
     first_edge = loops[0].GetEdges()[0]
     sample = _sample_edge_polyline(first_edge, chord_tol_mm=1.0)[0]
 
-    # Try direct normal from face API.
+    # Get face normal
     normal = (0.0, 0.0, 1.0)
     try:
-        n = face.GetFaceNormal(sample[0], sample[1], sample[2])
+        import NXOpen  # noqa
+        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:
-        pass
+        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:
+    """Try multiple NX API patterns to read a string attribute."""
-        return obj.GetStringUserAttribute(title, -1)
+    for method_name in ("GetStringUserAttribute", "GetUserAttributeAsString"):
-    except Exception:
+        try:
-        pass
+            method = getattr(obj, method_name)
-    try:
+            val = method(title, -1)
-        return obj.GetUserAttributeAsString(title, -1)
+            if val:
-    except Exception:
+                return str(val)
-        return None
+        except Exception:
+            continue
+    return None
 
 
-def find_sandbox_faces(work_part: Any, attr_name: str = "ISOGRID_SANDBOX") -> List[Tuple[str, Any]]:
+# ---------------------------------------------------------------------------
-    tagged: List[Tuple[str, Any]] = []
+# SIM → FEM → Idealized Part navigation
-    for body in getattr(work_part.Bodies, "ToArray", lambda: work_part.Bodies)():
+# ---------------------------------------------------------------------------
+
+def _navigate_sim_to_idealized(session: Any) -> Any:
+    """
+    From the active .sim work part, navigate to the idealized part.
+
+    Path: SIM → FEM component → idealized part (_i.prt)
+
+    Returns the idealized part (NXOpen.BasePart), and sets it as work part.
+    """
+    import NXOpen.CAE  # noqa
+
+    work_part = session.Parts.Work
+    part_name = work_part.Name if hasattr(work_part, "Name") else ""
+    part_path = work_part.FullPath if hasattr(work_part, "FullPath") else ""
+
+    print(f"[extract_sandbox] Starting from: {part_name}")
+    print(f"[extract_sandbox] Full path: {part_path}")
+
+    # Check if we're already in an idealized part
+    if "_i" in part_name.lower() or part_name.endswith("_i"):
+        print("[extract_sandbox] Already in idealized part.")
+        return work_part
+
+    # If in .sim, navigate to FEM then to idealized part
+    # Strategy: iterate through loaded parts to find the _i.prt
+    idealized_part = None
+    fem_part = None
+
+    for part in session.Parts:
+        pname = part.Name if hasattr(part, "Name") else ""
+        ppath = ""
+        try:
+            ppath = part.FullPath
+        except Exception:
+            pass
+
+        if "_i" in pname and pname.endswith("_i"):
+            # Found idealized part candidate
+            idealized_part = part
+            print(f"[extract_sandbox] Found idealized part: {pname}")
+        elif "_fem" in pname.lower() and not pname.endswith("_i"):
+            fem_part = part
+            print(f"[extract_sandbox] Found FEM part: {pname}")
+
+    if idealized_part is None:
+        # Try alternative: look for parts whose name contains "fem" + "_i"
+        for part in session.Parts:
+            pname = part.Name if hasattr(part, "Name") else ""
+            if "fem" in pname.lower() and "_i" in pname:
+                idealized_part = part
+                print(f"[extract_sandbox] Found idealized part (alt): {pname}")
+                break
+
+    if idealized_part is None:
+        raise RuntimeError(
+            "Could not find idealized part (_i.prt). "
+            "Ensure the SIM file is open and its FEM + idealized part are loaded."
+        )
+
+    # Set idealized part as work part
+    try:
+        session.Parts.SetWork(idealized_part)
+        print(f"[extract_sandbox] Set work part to: {idealized_part.Name}")
+    except Exception as exc:
+        print(f"[extract_sandbox] Warning: Could not set work part: {exc}")
+        print("[extract_sandbox] Proceeding with part reference anyway...")
+
+    return idealized_part
+
+
+# ---------------------------------------------------------------------------
+# Sandbox discovery — find sheet bodies with ISOGRID_SANDBOX attribute
+# ---------------------------------------------------------------------------
+
+def find_sandbox_bodies(
+    part: Any,
+    attr_name: str = "ISOGRID_SANDBOX",
+) -> List[Tuple[str, Any, Any]]:
+    """
+    Find sheet bodies tagged with ISOGRID_SANDBOX attribute.
+
+    Searches:
+      1. Body-level attributes
+      2. Face-level attributes (fallback)
+
+    Returns list of (sandbox_id, body, face) tuples.
+    The face is the primary face of the sheet body.
+    """
+    tagged: List[Tuple[str, Any, Any]] = []
+
+    bodies = []
+    try:
+        bodies = list(part.Bodies.ToArray()) if hasattr(part.Bodies, "ToArray") else list(part.Bodies)
+    except Exception:
+        bodies = list(part.Bodies)
+
+    for body in bodies:
+        # Check body-level attribute first
+        sandbox_id = _get_string_attribute(body, attr_name)
+        if sandbox_id:
+            # Sheet body — get its face(s)
+            faces = body.GetFaces()
+            if faces:
+                face = faces[0]  # Sheet body typically has one face
+                tagged.append((sandbox_id, body, face))
+                print(f"[extract_sandbox] Found sandbox body: {sandbox_id} "
+                      f"(body attr, {len(faces)} face(s))")
+                continue
+
+        # Fallback: check face-level attributes
         for face in body.GetFaces():
             sandbox_id = _get_string_attribute(face, attr_name)
-            if sandbox_id and sandbox_id.startswith("sandbox_"):
+            if sandbox_id:
-                tagged.append((sandbox_id, face))
+                tagged.append((sandbox_id, body, face))
+                print(f"[extract_sandbox] Found sandbox face: {sandbox_id} (face attr)")
+
+    if not tagged:
+        # Last resort: try matching body names
+        print(f"[extract_sandbox] No attributes found, 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:
+                    sandbox_id = bname.lower().replace(" ", "_")
+                    tagged.append((sandbox_id, body, faces[0]))
+                    print(f"[extract_sandbox] Found sandbox by name: {sandbox_id}")
+
     return tagged
 
 
-def _extract_face_loops(face: Any, chord_tol_mm: float, frame: LocalFrame) -> Tuple[List[Point2D], List[Dict[str, Any]]]:
+# ---------------------------------------------------------------------------
-    outer_2d: List[Point2D] = []
+# Core extraction
-    holes: List[Dict[str, Any]] = []
+# ---------------------------------------------------------------------------
+
+def extract_sandbox_geometry(
+    face: Any,
+    body: Any,
+    sandbox_id: str,
+    chord_tol_mm: float = 0.1,
+) -> Dict[str, Any]:
+    """
+    Extract a sandbox sheet face into a JSON-serializable dict.
+
+    - outer_boundary: the sandbox outline (2D polyline)
+    - inner_boundaries: cutouts from reserved cylinders etc. (2D polylines)
+      These are boundary constraints, NOT optimization holes.
+    - transform: 3D ↔ 2D coordinate mapping
+    """
+    frame = _face_local_frame(face)
+
+    outer_2d: List[List[float]] = []
+    inner_boundaries: List[Dict[str, Any]] = []
 
     loops = face.GetLoops()
+    print(f"[extract_sandbox] {sandbox_id}: {len(loops)} loop(s) on face")
+
     for loop_index, loop in enumerate(loops):
+        # Collect all edge points for this loop
         loop_pts3d: List[Point3D] = []
-        for edge in loop.GetEdges():
+        edges = loop.GetEdges()
+        for edge in edges:
             pts = _sample_edge_polyline(edge, chord_tol_mm)
             if loop_pts3d and pts:
-                pts = pts[1:]  # avoid duplicate joining point
+                pts = pts[1:]  # avoid duplicate at junction
             loop_pts3d.extend(pts)
 
         loop_pts3d = _close_polyline(loop_pts3d)
         loop_pts2d = project_to_2d(loop_pts3d, frame)
 
+        # Determine outer vs inner
         is_outer = False
         try:
             is_outer = loop.IsOuter()
@@ -228,89 +369,119 @@ def _extract_face_loops(face: Any, chord_tol_mm: float, frame: LocalFrame) -> Tu
             is_outer = (loop_index == 0)
 
         if is_outer:
-            outer_2d = loop_pts2d
+            outer_2d = [[round(x, 6), round(y, 6)] for x, y in loop_pts2d]
-            continue
+            print(f"[extract_sandbox]   outer loop: {len(outer_2d)} points")
+        else:
+            boundary = [[round(x, 6), round(y, 6)] for x, y in loop_pts2d]
+            inner_boundaries.append({
+                "index": len(inner_boundaries),
+                "boundary": boundary,
+                "num_points": len(boundary),
+            })
+            print(f"[extract_sandbox]   inner loop {len(inner_boundaries)}: "
+                  f"{len(boundary)} points")
 
-        is_circ, center, diameter = is_loop_circular(loop_pts2d)
+    # Try to get thickness from body or expression
-        holes.append(
+    thickness = None
-            {
-                "index": len(holes),
-                "boundary": [[x, y] for x, y in loop_pts2d],
-                "center": [center[0], center[1]] if is_circ else None,
-                "diameter": diameter if is_circ else None,
-                "is_circular": bool(is_circ),
-                "weight": 0.0,
-            }
-        )
-
-    return outer_2d, holes
-
-
-def extract_sandbox_geometry(face: Any, sandbox_id: str, chord_tol_mm: float = 0.1) -> Dict[str, Any]:
-    frame = _face_local_frame(face)
-    outer, holes = _extract_face_loops(face, chord_tol_mm=chord_tol_mm, frame=frame)
-
-    geom = {
-        "units": "mm",
-        "sandbox_id": sandbox_id,
-        "outer_boundary": [[x, y] for x, y in outer],
-        "holes": holes,
-        "transform": {
-            "origin": list(frame.origin),
-            "x_axis": list(frame.x_axis),
-            "y_axis": list(frame.y_axis),
-            "normal": list(frame.normal),
-        },
-    }
-
-    # Optional thickness hint if available.
     try:
-        geom["thickness"] = float(face.GetBody().GetThickness())
+        # NX sheet bodies may expose thickness
+        thickness = float(body.GetThickness())
     except Exception:
         pass
 
+    geom: Dict[str, Any] = {
+        "schema_version": "1.0",
+        "units": "mm",
+        "sandbox_id": sandbox_id,
+        "outer_boundary": outer_2d,
+        "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],
+        },
+    }
+
     return geom
 
 
-def export_sandbox_geometries(output_dir: Path, geometries: Dict[str, Dict[str, Any]]) -> List[Path]:
+# ---------------------------------------------------------------------------
+# Export
+# ---------------------------------------------------------------------------
+
+def export_sandbox_geometries(
+    output_dir: Path,
+    geometries: Dict[str, Dict[str, Any]],
+) -> List[Path]:
     output_dir.mkdir(parents=True, exist_ok=True)
     written: List[Path] = []
     for sandbox_id, payload in geometries.items():
         out = output_dir / f"geometry_{sandbox_id}.json"
         out.write_text(json.dumps(payload, indent=2))
         written.append(out)
+        print(f"[extract_sandbox] Wrote: {out}")
     return written
 
 
+# ---------------------------------------------------------------------------
+# Main entry point — runs inside NX
+# ---------------------------------------------------------------------------
+
 def run_in_nx(output_dir: Path, chord_tol_mm: float = 0.1) -> List[Path]:
     import NXOpen  # type: ignore
 
     session = NXOpen.Session.GetSession()
-    work_part = session.Parts.Work
-    if work_part is None:
-        raise RuntimeError("No active NX work part.")
 
-    sandbox_faces = find_sandbox_faces(work_part)
+    # Navigate from .sim to idealized part
-    if not sandbox_faces:
+    idealized_part = _navigate_sim_to_idealized(session)
-        raise RuntimeError("No faces found with ISOGRID_SANDBOX attribute.")
 
+    # Find sandbox bodies
+    sandbox_entries = find_sandbox_bodies(idealized_part)
+    if not sandbox_entries:
+        raise RuntimeError(
+            "No sandbox bodies found. Ensure sheet bodies have "
+            "the ISOGRID_SANDBOX attribute set (body or face level)."
+        )
+
+    print(f"[extract_sandbox] Found {len(sandbox_entries)} sandbox(es)")
+
+    # Extract geometry for each sandbox
     payloads: Dict[str, Dict[str, Any]] = {}
-    for sandbox_id, face in sandbox_faces:
+    for sandbox_id, body, face in sandbox_entries:
-        payloads[sandbox_id] = extract_sandbox_geometry(face, sandbox_id, chord_tol_mm=chord_tol_mm)
+        payloads[sandbox_id] = extract_sandbox_geometry(
+            face=face,
+            body=body,
+            sandbox_id=sandbox_id,
+            chord_tol_mm=chord_tol_mm,
+        )
 
+    # Export
     return export_sandbox_geometries(output_dir=output_dir, geometries=payloads)
 
 
 def main(argv: Sequence[str] | None = None) -> int:
-    parser = argparse.ArgumentParser(description="Extract NX sandbox face geometry to JSON")
+    parser = argparse.ArgumentParser(
-    parser.add_argument("--output-dir", default=".", help="Directory for geometry_sandbox_*.json")
+        description="Extract NX sandbox sheet geometry to JSON (run from .sim)"
-    parser.add_argument("--chord-tol", type=float, default=0.1, help="Edge sampling chord tolerance (mm)")
+    )
+    parser.add_argument(
+        "--output-dir", default=".",
+        help="Directory for geometry_sandbox_*.json output files",
+    )
+    parser.add_argument(
+        "--chord-tol", type=float, default=0.1,
+        help="Edge sampling chord tolerance in mm (default: 0.1)",
+    )
     args = parser.parse_args(argv)
 
     out_dir = Path(args.output_dir)
     written = run_in_nx(output_dir=out_dir, chord_tol_mm=args.chord_tol)
+
+    print(f"\n[extract_sandbox] Done — {len(written)} file(s) exported.")
     for p in written:
-        print(f"[extract_sandbox] wrote: {p}")
+        print(f"  → {p}")
     return 0