diff --git a/tools/adaptive-isogrid/docs/technical-spec.md b/tools/adaptive-isogrid/docs/technical-spec.md
index 7a6599c5..0e36a24e 100644
--- a/tools/adaptive-isogrid/docs/technical-spec.md
+++ b/tools/adaptive-isogrid/docs/technical-spec.md
@@ -21,11 +21,11 @@ After extensive brainstorming, the following decisions are locked:
 | Decision | Choice | Rationale |
 |---|---|---|
 | Geometry generation | External Python (Constrained Delaunay) | Full access to scipy/triangle/gmsh, debuggable, fast |
-| FEA strategy | Assembly FEM with superposed models | Decouples fixed interfaces (loads/BCs) from variable rib pattern |
+| FEA strategy | Reserved-region monolithic remesh | Keep load/BC topology stable while allowing local rib updates |
 | FEA solver | NX Simcenter + Nastran (2D shell) | Existing expertise, handles complex BCs, extensible to modal/buckling |
-| NX role | Remesh plate model + merge nodes + solve | Geometry is agnostic to rib pattern; loads/BCs never re-associated |
+| NX role | Extract sandbox faces, reimport profile, remesh + solve | Reserved regions preserve associations; no assembly merge pipeline needed |
 | Optimization | Atomizer (Optuna TPE), pure parametric v1 | One FEA per trial, ~2 min/iteration, stress feedback deferred to v2 |
-| Manufacturing | Through-cuts only (waterjet + CNC finish) | Simplifies geometry to 2D profile problem |
+| Geometry transfer | JSON-only round trip | Deterministic, scriptable, no DXF/STEP conversion drift |
 | Plate type | Flat, 200–600 mm scale, 6–15 mm thick, 16–30 holes | Shell elements appropriate, fast solves |
 
 ### System Diagram
@@ -34,27 +34,11 @@ After extensive brainstorming, the following decisions are locked:
 ONE-TIME SETUP
 ══════════════
 User in NX:
-  ├── Select plate face → tool extracts boundary + holes → geometry.json
+  ├── Partition plate mid-surface into regions:
+  │   ├── sandbox face(s): ISOGRID_SANDBOX = sandbox_1, sandbox_2, ...
+  │   └── reserved regions: edges/functional zones to stay untouched
   ├── Assign hole weights (interactive UI or table)
-  │
-  ├── Build Assembly FEM:
-  │   ├── Model A — "Interface Model" (PERMANENT, never changes)
-  │   │   ├── Spider elements (RBE2/RBE3) at each hole center → circumference
-  │   │   ├── All loads applied to spider center nodes
-  │   │   ├── All BCs applied to spider center nodes or edge nodes
-  │   │   └── Edge BC nodes along plate boundary
-  │   │
-  │   ├── Model B — "Plate Model" (VARIABLE, rebuilt each iteration)
-  │   │   ├── 2D shell mesh of the ribbed plate profile
-  │   │   ├── Mesh seeds at hole circumference nodes (match Model A spiders)
-  │   │   └── Mesh seeds at edge BC nodes (match Model A edge nodes)
-  │   │
-  │   └── Assembly FEM
-  │       ├── Superpose Model A + Model B
-  │       ├── Merge nodes at hole circumferences + edges
-  │       └── Solver settings (SOL 101, SOL 103, etc.)
-  │
-  └── Verify: solve dummy Model B (solid plate), confirm results
+  └── Verify baseline solve and saved solution setup
 
 OPTIMIZATION LOOP (~2 min/iteration)
 ════════════════════════════════════
@@ -64,58 +48,49 @@ Atomizer (Optuna TPE)
   │  (η₀, α, β, p, R₀, κ, s_min, s_max, t_min, t₀, γ, w_frame, r_f, ...)
   │
   ▼
+NXOpen Extractor (~seconds)
+  ├── Find all sandbox faces by attribute
+  ├── Extract each sandbox loop set to local 2D
+  └── Write geometry_sandbox_n.json
+  │
+  ▼
 External Python — "The Brain" (~1-3 sec)
-  ├── Load plate geometry + hole table (fixed)
+  ├── Load geometry_sandbox_n.json
   ├── Compute density field η(x)
   ├── Generate constrained Delaunay triangulation
-  ├── Compute rib thicknesses from density
   ├── Apply manufacturing constraints
-  ├── Validate geometry (min web, min pocket, no degenerates)
-  ├── Output: rib_profile.json (curve coordinates)
-  └── Output: validity_flag + mass_estimate
+  └── Write rib_profile_sandbox_n.json
   │
   ▼
-NXOpen Journal — "The Hands" (~60-90 sec)
-  ├── Delete old Model B geometry + mesh
-  ├── Import new 2D ribbed profile (from rib_profile.json)
-  ├── Create sheet body from profile curves
-  ├── Mesh Model B with seeds at fixed interface nodes
-  ├── Merge nodes in Assembly FEM (holes + edges)
+NXOpen Import + Solve (~60-90 sec)
+  ├── Replace sandbox geometry only
+  ├── Keep reserved regions unchanged (loads/BCs persist)
+  ├── Remesh full plate as monolithic model
   ├── Solve (Nastran)
-  └── Extract results (stress, displacement, strain, mass)
-  │
-  ▼
-Result Extraction (~5-10 sec)
-  ├── Von Mises stress field (nodal)
-  ├── Displacement magnitude field (nodal)
-  ├── Strain field (elemental)
-  ├── Total mass
-  ├── (Future: modal frequencies, buckling load factors)
-  └── Report metrics to Atomizer
+  └── Extract results.json (stress/disp/strain/mass)
   │
   ▼
 Atomizer updates surrogate, samples next trial
 Repeat until convergence (~500-2000 trials)
 ```
 
-### Key Architectural Insight: Why Assembly FEM
+### Key Architectural Insight: Why Reserved Regions
 
 The plate lightweighting problem has a natural separation:
 
-**What doesn't change:** Hole positions, hole diameters, plate boundary, loads, boundary conditions. These are the structural interfaces — where forces enter and leave the plate.
+**What should stay fixed:** load/BC attachment topology around edges, holes, and functional interfaces.
 
-**What changes every iteration:** The rib pattern between holes. This is purely the internal load-path topology.
+**What should evolve:** the rib network inside the designable interior.
 
-The Assembly FEM with superposed models exploits this separation directly. Model A captures everything fixed (interfaces, loads, BCs) using spider elements at holes and edge nodes at boundaries. Model B captures everything variable (the ribbed plate mesh). They connect through node merging at the fixed interface locations (hole circumferences and plate edges).
+The reserved-region workflow maps directly to this separation. By only replacing
+sandbox faces, all references in reserved zones remain valid while the interior
+rib topology changes every iteration.
 
 This means:
-- Loads and BCs **never need re-association** — they live permanently on Model A
-- Only Model B's mesh is rebuilt each iteration
-- Node merging at fixed geometric locations is a reliable, automated operation in NX
-- The solver sees one connected model with proper load paths through the spiders into the ribs
-Atomizer updates surrogate, samples next trial
-Repeat until convergence (~500-2000 trials)
-```
+- Loads and BCs persist naturally because reserved geometry is untouched
+- No multi-model coupling layer or interface node reconciliation workflow
+- Full-plate remesh remains straightforward and robust
+- Geometry extraction runs each iteration, enabling hole-position optimization
 
 ---
 
@@ -577,341 +552,100 @@ These constraints are enforced during geometry generation, not as FEA post-check
 
 ---
 
-## 4. The NX Hands: Assembly FEM Architecture
+## 4. The NX Hands: Reserved-Region FEM Architecture
 
-### 4.1 The Two-Model Structure
+### 4.1 Core Concept: Sandbox + Reserved Regions
 
-The FEA uses an Assembly FEM (AFEM) with two superposed finite element models. This is the standard aerospace approach for decoupling load introduction from structural detail, and it perfectly fits our problem where interfaces are fixed but internal topology varies.
+The FEA workflow uses direct reserved-region remeshing rather than interface-coupling constructs.
 
-**Model A — Interface Model (built once, never modified):**
+Instead, each plate mid-surface is explicitly partitioned in NX into:
 
-Model A contains only the structural interface elements. It has no plate geometry — just the load introduction and boundary condition hardware.
+- **Sandbox region(s)** — where the adaptive isogrid is allowed to change each iteration
+- **Reserved regions** — geometry that must remain untouched (edges, functional features, local areas around critical holes/BC interfaces)
 
-For each hole:
-- One node at the hole center (the "master" node)
-- N nodes equally spaced on the hole circumference (N = 12–24 depending on hole diameter, typically 1 node per ~2 mm of circumference)
-- RBE2 or RBE3 spider element connecting center node to circumference nodes
-  - RBE2 (rigid): for bolted/pinned connections where the hole is rigidly constrained
-  - RBE3 (distributing): for bearing loads or distributed force introduction
+Each sandbox face is tagged with an NX user attribute:
 
-For plate boundary (edge BCs):
-- Nodes distributed along the plate outer boundary at regular spacing (~2–5 mm)
-- These serve as merge targets for Model B's edge mesh nodes
+- **Title:** `ISOGRID_SANDBOX`
+- **Value:** `sandbox_1`, `sandbox_2`, ...
 
-All loads and boundary conditions are applied to Model A's nodes:
-- Bolt forces → hole center nodes
-- Fixed constraints → hole center nodes or edge nodes
-- Bearing loads → hole center nodes (via RBE3)
-- Enforced displacements → relevant center/edge nodes
-- Pressures → applied directly to Model B elements (only exception)
+Multiple sandbox faces are supported on the same plate.
 
-**Model B — Plate Model (rebuilt each iteration):**
+### 4.2 Iterative Geometry Round-Trip (JSON-only)
 
-Model B is the 2D shell mesh of the current ribbed plate profile. It is the only model that changes during optimization. Key requirements:
+Per optimization iteration, NXOpen performs a full geometry round-trip for each sandbox:
 
-- Shell elements (CQUAD4/CTRIA3) with PSHELL property at plate thickness
-- Mesh must have nodes at exact locations matching Model A's spider circumference nodes and edge BC nodes
-- These "interface nodes" are enforced via mesh seed points (hard nodes) in NX's mesher
-- The rest of the mesh fills in freely, adapting to the rib pattern geometry
+1. **Extract current sandbox geometry** from NX into `geometry_sandbox_n.json`
+   - outer loop boundary
+   - inner loops (holes/cutouts)
+   - local 2D transform (face-local XY frame)
+2. **Python Brain generates isogrid profile** inside sandbox boundary
+   - outputs `rib_profile_sandbox_n.json`
+3. **NXOpen imports profile** and replaces only sandbox geometry
+   - reserved faces remain unchanged
+4. **Full plate remesh** (single monolithic mesh)
+5. **Solve existing Simcenter solution setup**
+6. **Extract field + scalar results** to `results.json`
 
-**Assembly FEM:**
-- Superimposes Model A and Model B
-- Merges coincident nodes at hole circumferences and plate edges (tolerance ~0.01 mm)
-- After merge, loads flow from Model A's spiders through the merged nodes into Model B's shell mesh
-- Solver settings, solution sequences, and output requests live at the assembly level
+This JSON-only flow replaces DXF/STEP transfer for v1 and keeps CAD/CAE handoff deterministic.
 
-### 4.2 One-Time Setup Procedure
+### 4.3 Why Reserved-Region Is Robust
 
-This is done once per plate project, before any optimization runs.
+This architecture addresses the load/BC persistence problem by preserving topology where constraints live:
 
-**Step 1 — Extract geometry and build Model A:**
+- Loads and boundary conditions are assigned to entities in reserved regions that do not change between iterations
+- Only sandbox geometry is replaced, so reserved-region references remain valid
+- The solver always runs on one connected, remeshed plate model (no multi-model merge operations)
+- Hole positions can be optimized because geometry extraction happens **every iteration** from the latest NX model state
 
-```python
-# NXOpen setup script — build the interface model
-import NXOpen
-import json
-import math
+### 4.4 NXOpen Phase-2 Script Responsibilities
 
-def build_interface_model(geometry_json_path, fem_part):
-    """
-    Build Model A: spider elements at each hole + edge BC nodes.
-    """
-    with open(geometry_json_path) as f:
-        geometry = json.load(f)
+#### A) `extract_sandbox.py`
 
-    interface_nodes = {
-        'hole_centers': [],     # (hole_idx, node_id, x, y)
-        'hole_circumferences': [],  # (hole_idx, node_id, x, y)
-        'edge_nodes': []        # (node_id, x, y)
-    }
+- Discover all faces with `ISOGRID_SANDBOX` attribute
+- For each sandbox face:
+  - read outer + inner loops
+  - sample edges to polyline (chord tolerance 0.1 mm)
+  - fit circles on inner loops when circular
+  - project 3D points to face-local 2D
+  - write `geometry_sandbox_n.json` using Python Brain input schema
 
-    # --- Create spider elements for each hole ---
-    for hole in geometry['holes']:
-        cx, cy = hole['center']
-        radius = hole['diameter'] / 2.0
+#### B) `import_profile.py`
 
-        # Create center node
-        center_node = create_node(fem_part, cx, cy, 0.0)
-        interface_nodes['hole_centers'].append({
-            'hole_index': hole['index'],
-            'node_id': center_node.Label,
-            'x': cx, 'y': cy
-        })
+- Read `rib_profile_sandbox_n.json`
+- Rebuild NX curves from coordinate arrays
+- Create/update sheet body for the sandbox zone
+- Replace sandbox face geometry while preserving surrounding reserved faces
+- Sew/unite resulting geometry into a watertight plate body
 
-        # Create circumference nodes
-        n_circ = max(12, int(math.pi * hole['diameter'] / 2.0))  # ~1 node per 2mm
-        circ_nodes = []
-        for j in range(n_circ):
-            angle = 2 * math.pi * j / n_circ
-            nx_ = cx + radius * math.cos(angle)
-            ny_ = cy + radius * math.sin(angle)
-            node = create_node(fem_part, nx_, ny_, 0.0)
-            circ_nodes.append(node)
-            interface_nodes['hole_circumferences'].append({
-                'hole_index': hole['index'],
-                'node_id': node.Label,
-                'x': nx_, 'y': ny_
-            })
+#### C) `solve_and_extract.py`
 
-        # Create RBE2 spider (rigid) — default for mounting holes
-        # Use RBE3 (distributing) for bearing load holes
-        spider_type = 'RBE2'  # could be per-hole from hole table
-        create_spider(fem_part, center_node, circ_nodes, spider_type)
+- Regenerate FE mesh for full plate
+- Trigger solve using existing solution object(s)
+- Extract from sandbox-associated nodes/elements:
+  - nodal Von Mises stress
+  - nodal displacement magnitude
+  - elemental strain
+  - total mass
+- Write `results.json`:
 
-    # --- Create edge BC nodes ---
-    boundary = geometry['outer_boundary']
-    edge_spacing = 3.0  # mm between edge nodes
-    for i in range(len(boundary)):
-        p1 = boundary[i]
-        p2 = boundary[(i + 1) % len(boundary)]
-        edge_length = distance(p1, p2)
-        n_nodes = max(2, int(edge_length / edge_spacing))
-        for j in range(n_nodes):
-            t = j / n_nodes
-            ex = p1[0] + t * (p2[0] - p1[0])
-            ey = p1[1] + t * (p2[1] - p1[1])
-            node = create_node(fem_part, ex, ey, 0.0)
-            interface_nodes['edge_nodes'].append({
-                'node_id': node.Label,
-                'x': ex, 'y': ey
-            })
-
-    # Save interface node map for the iteration script
-    with open('interface_nodes.json', 'w') as f:
-        json.dump(interface_nodes, f, indent=2)
-
-    return interface_nodes
+```json
+{
+  "nodes_xy": [[...], [...]],
+  "stress_values": [...],
+  "disp_values": [...],
+  "strain_values": [...],
+  "mass": 0.0
+}
 ```
 
-**Step 2 — Apply loads and BCs to Model A:**
+### 4.5 Results Extraction Strategy
 
-This is done manually in NX Simcenter (or scripted for standard load cases). The user applies all structural loads and boundary conditions to Model A's center nodes and edge nodes. These never change.
+Two output backends are acceptable:
 
-Examples:
-- Bolt M8 at hole 3: axial force 5000 N on center node of hole 3
-- Fixed constraint at holes 0, 1: fix all DOFs on center nodes of holes 0, 1
-- Simply supported edge: constrain Z-displacement on all edge nodes along one plate edge
-- Bearing load at hole 7: 2000 N in X-direction via RBE3 at hole 7 center
+1. **Direct NXOpen post-processing API** (preferred when available)
+2. **Simcenter CSV export + parser** (robust fallback)
 
-**Step 3 — Build dummy Model B and verify:**
-
-Create a simple Model B (e.g., the un-lightweighted plate meshed as shells), set up the assembly FEM with merging, and solve to verify the setup is correct. Compare results against a monolithic (non-assembly) FEA to confirm the spider elements and merging behave as expected.
-
-**Step 4 — Save the template:**
-
-Save Model A, the assembly FEM structure, and the interface node map. This is the reusable template for all optimization iterations.
-
-### 4.3 Per-Iteration NXOpen Journal Script
-
-This is the script that runs inside the optimization loop. It receives the ribbed plate profile from the Python brain and handles Model B rebuild + solve.
-
-```python
-# NXOpen iteration script — rebuild Model B, merge, solve, extract
-import NXOpen
-import json
-import numpy as np
-
-def iteration_solve(profile_path, interface_nodes_path, afem_part):
-    """
-    Single optimization iteration:
-    1. Delete old Model B geometry + mesh
-    2. Import new profile
-    3. Mesh with interface node seeds
-    4. Merge nodes in AFEM
-    5. Solve
-    6. Extract results
-    """
-    session = NXOpen.Session.GetSession()
-
-    # Load inputs
-    with open(profile_path) as f:
-        profile = json.load(f)
-    with open(interface_nodes_path) as f:
-        interface_nodes = json.load(f)
-
-    if not profile['valid']:
-        return {'status': 'invalid_geometry', 'mass': float('inf')}
-
-    # --- Step 1: Clean old Model B ---
-    model_b = get_model_b_fem(afem_part)
-    delete_all_mesh(model_b)
-    delete_all_geometry(model_b)
-
-    # --- Step 2: Import new 2D profile ---
-    # Create curves from profile coordinate arrays
-    outer_coords = profile['outer_boundary']
-    create_closed_polyline(model_b, outer_coords)
-
-    for pocket_coords in profile['pockets']:
-        create_closed_polyline(model_b, pocket_coords)
-
-    # Create sheet body from bounded regions
-    sheet_body = create_sheet_from_curves(model_b)
-
-    # --- Step 3: Mesh with interface seeds ---
-    mesh_control = get_mesh_collector(model_b)
-
-    # Add hard-point mesh seeds at all interface locations
-    # These force the mesher to place nodes exactly where
-    # Model A's spiders attach
-    for node_info in interface_nodes['hole_circumferences']:
-        add_mesh_seed_point(mesh_control, node_info['x'], node_info['y'])
-
-    for node_info in interface_nodes['edge_nodes']:
-        add_mesh_seed_point(mesh_control, node_info['x'], node_info['y'])
-
-    # Set mesh parameters
-    set_element_size(mesh_control, target=2.0, min_size=0.5)  # mm
-    set_element_type(mesh_control, 'CQUAD4')  # prefer quads, allow tris
-
-    # Generate mesh
-    mesh_control.GenerateMesh()
-
-    # --- Step 4: Merge nodes in Assembly FEM ---
-    afem = get_assembly_fem(afem_part)
-    merge_coincident_nodes(afem, tolerance=0.05)  # mm
-
-    # Verify merge count matches expected
-    expected_merges = (
-        len(interface_nodes['hole_circumferences']) +
-        len(interface_nodes['edge_nodes'])
-    )
-    actual_merges = get_merge_count(afem)
-    if actual_merges < expected_merges * 0.95:
-        # Merge failed for some nodes — flag as warning
-        print(f"WARNING: Expected {expected_merges} merges, got {actual_merges}")
-
-    # --- Step 5: Solve ---
-    solution = get_solution(afem, 'static_analysis')
-    solve_result = solution.Solve()
-
-    if not solve_result.success:
-        return {'status': 'solve_failed', 'mass': float('inf')}
-
-    # --- Step 6: Extract results ---
-    results = extract_results(afem, solution)
-
-    return results
-
-
-def extract_results(afem, solution):
-    """
-    Extract field results from the solved assembly FEM.
-    Only extracts from Model B elements (the plate mesh),
-    ignoring Model A's spider elements.
-    """
-    post = get_post_processor(afem)
-
-    # Get results only from Model B's element group
-    model_b_elements = get_model_b_element_group(afem)
-
-    # Von Mises stress (nodal, averaged at nodes)
-    stress_data = post.GetNodalResults(
-        solution,
-        result_type='Stress',
-        component='Von Mises',
-        element_group=model_b_elements
-    )
-
-    # Displacement magnitude (nodal)
-    disp_data = post.GetNodalResults(
-        solution,
-        result_type='Displacement',
-        component='Magnitude',
-        element_group=model_b_elements
-    )
-
-    # Strain (elemental)
-    strain_data = post.GetElementalResults(
-        solution,
-        result_type='Strain',
-        component='Von Mises',
-        element_group=model_b_elements
-    )
-
-    # Mass from Model B mesh (plate material only, not spiders)
-    mass = compute_shell_mass(model_b_elements)
-
-    results = {
-        'status': 'solved',
-        'mass': mass,
-        'max_von_mises': float(np.max(stress_data['values'])),
-        'max_displacement': float(np.max(disp_data['values'])),
-        'mean_von_mises': float(np.mean(stress_data['values'])),
-        'stress_field': {
-            'nodes_xy': stress_data['coordinates'].tolist(),
-            'values': stress_data['values'].tolist()
-        },
-        'displacement_field': {
-            'nodes_xy': disp_data['coordinates'].tolist(),
-            'values': disp_data['values'].tolist()
-        },
-        'strain_field': {
-            'elements_xy': strain_data['centroids'].tolist(),
-            'values': strain_data['values'].tolist()
-        }
-    }
-
-    return results
-```
-
-### 4.4 Why This Approach Is Robust
-
-The AFEM node-merge strategy solves the hardest problem in automated FEA iteration: **load and BC persistence across geometry changes.** Here's why it works reliably:
-
-**Fixed merge locations:** Hole centers, hole circumferences, and plate edges don't move between iterations. The merge is always at the same physical coordinates. NX's node merge by tolerance is a simple, reliable geometric operation.
-
-**Mesh seed enforcement:** By placing hard-point seeds at all interface locations, the mesher is forced to create nodes at exactly those coordinates. This guarantees that every merge finds its partner node.
-
-**Rib pattern agnostic:** Model A doesn't know or care what the rib pattern looks like. It only knows where the holes and edges are. Whether there are 50 or 200 pockets, the spiders connect the same way.
-
-**Easy validation:** After merging, a simple check (did we get the expected number of merged node pairs?) catches any meshing or geometry issues before wasting time on a solve.
-
-**Extensible:** Adding new load cases, new holes, or new BC types only requires modifying Model A (once). The optimization loop and Model B generation are unaffected.
-
-### 4.5 Simcenter Configuration Details
-
-**Shell property (PSHELL):**
-- Thickness: from `geometry.json` (e.g., 10.0 mm)
-- Material: MAT1 referencing the plate material (e.g., AL6061-T6)
-- Applied to all Model B elements
-
-**Spider elements:**
-- RBE2 (rigid): 6 DOF coupling from center to circumference nodes. Use for fixed/bolted connections where the hole acts as a rigid interface.
-- RBE3 (weighted average): Distributes loads from center to circumference. Use for bearing loads where the hole deforms under load and you want realistic load distribution.
-- Choice is per-hole, set during one-time setup based on connection type.
-
-**Mesh controls for Model B:**
-- Target element size: 1.5–3.0 mm (captures rib geometry adequately)
-- Minimum element size: 0.5 mm (allows refinement at narrow rib junctions)
-- Element type: CQUAD4 dominant with CTRIA3 fill
-- Mesh seeds: hard points at all interface node locations
-- Edge mesh control on hole circumferences: N elements matching spider node count
-
-**Solution sequences:**
-- SOL 101: Static analysis (v1)
-- SOL 103: Normal modes / modal analysis (v2, for natural frequency constraints)
-- SOL 105: Buckling (v2, for thin-rib stability checks)
+Both must produce consistent arrays for downstream optimization and optional stress-feedback loops.
 
 ---
 
@@ -1054,33 +788,32 @@ Build and test the geometry generator independently of NX:
 
 **Deliverable:** A Python module that takes `geometry.json` + parameters → outputs `rib_profile.json` + visualization plots.
 
-### Phase 2 — NX Geometry Extraction + AFEM Setup Scripts (1-2 weeks)
+### Phase 2 — NX Sandbox Extraction + Profile Import Scripts (1-2 weeks)
 
-Build the NXOpen scripts for one-time project setup:
+Build the NXOpen scripts for reserved-region geometry handling:
 
-- Face selection and hole detection script → exports `geometry.json`
-- Hole weight assignment UI (or table import)
-- Model A builder: spider elements at all holes + edge BC nodes → exports `interface_nodes.json`
-- Assembly FEM creation with Model A + dummy Model B
-- Load/BC application to Model A (manual or scripted for standard cases)
-- Verification solve on dummy Model B
-- Save as reusable template
+- Sandbox discovery via NX attribute (`ISOGRID_SANDBOX = sandbox_n`)
+- Per-sandbox extraction script → `geometry_sandbox_n.json`
+- Local 2D projection + loop sampling (outer + inner loops)
+- Rib profile import script (`rib_profile_sandbox_n.json`)
+- Sandbox-only geometry replacement + sew/unite with reserved regions
+- End-to-end JSON round-trip validation on multi-sandbox plates
 
-**Deliverable:** Complete one-time setup pipeline. Given any plate model, produces the AFEM template ready for optimization.
+**Deliverable:** Complete NX geometry round-trip pipeline (extract/import) with reserved regions preserved.
 
-### Phase 3 — NX Iteration Script (1-2 weeks)
+### Phase 3 — NX Iteration Solve + Results Extraction (1-2 weeks)
 
-Build the NXOpen journal script for the per-iteration loop:
+Build the NXOpen per-iteration analysis script:
 
-- Model B geometry cleanup (delete old mesh + geometry)
-- Profile import from `rib_profile.json` → NX curves → sheet body
-- Mesh with hard-point seeds at interface node locations
-- Assembly node merge with verification
-- Nastran solve trigger
-- Result extraction (stress/displacement/strain fields + scalar metrics) → `results.json`
-- End-to-end test: Python Brain → NX journal → results → validate against manual FEA
+- Extract sandbox geometry every iteration (supports moving/optimized holes)
+- Import regenerated rib profile into sandbox region(s)
+- Remesh the **full** plate as one monolithic model
+- Trigger existing Simcenter solution setup
+- Extract nodal stress/displacement + elemental strain + mass
+- Serialize standardized `results.json` for Atomizer
+- End-to-end test: Python Brain → NX scripts → `results.json` vs manual benchmark
 
-**Deliverable:** Complete single-iteration pipeline, verified against known-good manual analysis.
+**Deliverable:** Production-ready per-iteration NX pipeline with stable result export for optimization.
 
 ### Phase 4 — Atomizer Integration (1 week)
 
diff --git a/tools/adaptive-isogrid/src/nx/__init__.py b/tools/adaptive-isogrid/src/nx/__init__.py
index b1e66c31..8a6572ab 100644
--- a/tools/adaptive-isogrid/src/nx/__init__.py
+++ b/tools/adaptive-isogrid/src/nx/__init__.py
@@ -1,6 +1,9 @@
 """
 Adaptive Isogrid — NX Hands
 
-NXOpen journal scripts for geometry extraction, AFEM setup, and per-iteration solve.
-These scripts run inside NX Simcenter via the NXOpen Python API.
+Reserved-region NXOpen scripts:
+- extract_sandbox.py: sandbox loop extraction to geometry JSON
+- import_profile.py: profile reimport and sandbox replacement
+- solve_and_extract.py: remesh, solve, and result export
+- run_iteration.py: one-iteration orchestrator
 """
diff --git a/tools/adaptive-isogrid/src/nx/extract_sandbox.py b/tools/adaptive-isogrid/src/nx/extract_sandbox.py
new file mode 100644
index 00000000..eec6af39
--- /dev/null
+++ b/tools/adaptive-isogrid/src/nx/extract_sandbox.py
@@ -0,0 +1,318 @@
+"""
+NXOpen script — extract sandbox face geometry for Adaptive Isogrid.
+
+Finds faces tagged with user attribute:
+    ISOGRID_SANDBOX = sandbox_1, sandbox_2, ...
+
+For each sandbox face, exports `geometry_<sandbox_id>.json` in the same schema
+expected by the Python Brain (`outer_boundary`, `holes`, transform metadata, etc.).
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import math
+from dataclasses import dataclass
+from pathlib import Path
+from typing import Any, Dict, Iterable, List, Sequence, Tuple
+
+
+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 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:
+        v = _sub(p, frame.origin)
+        out.append((_dot(v, frame.x_axis), _dot(v, frame.y_axis)))
+    return out
+
+
+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
+
+
+def _sample_edge_polyline(edge: Any, chord_tol_mm: float) -> List[Point3D]:
+    """
+    Sample an NX edge as a polyline.
+    NOTE: NX APIs vary by curve type; this helper intentionally keeps fallback logic.
+    """
+    # Preferred path: use evaluator where available.
+    try:
+        evaluator = edge.CreateEvaluator()
+        t0, t1 = evaluator.GetLimits()
+        length = edge.GetLength()
+        n = max(2, int(length / max(chord_tol_mm, 1e-3)))
+        pts: List[Point3D] = []
+        for i in range(n + 1):
+            t = t0 + (t1 - t0) * (i / n)
+            p, _ = evaluator.Evaluate(t)
+            pts.append((float(p.X), float(p.Y), float(p.Z)))
+        return pts
+    except Exception:
+        pass
+
+    # Fallback: edge vertices only (less accurate, but safe fallback).
+    try:
+        verts = edge.GetVertices()
+        pts = []
+        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 _face_local_frame(face: Any) -> LocalFrame:
+    """
+    Build a stable local frame on a face:
+    - origin: first loop first point
+    - normal: face normal near origin
+    - x/y axes: orthonormal basis on tangent plane
+    """
+    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.
+    normal = (0.0, 0.0, 1.0)
+    try:
+        n = face.GetFaceNormal(sample[0], sample[1], sample[2])
+        normal = _normalize((float(n.X), float(n.Y), float(n.Z)))
+    except Exception:
+        pass
+
+    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)
+
+
+def _get_string_attribute(obj: Any, title: str) -> str | None:
+    try:
+        return obj.GetStringUserAttribute(title, -1)
+    except Exception:
+        pass
+    try:
+        return obj.GetUserAttributeAsString(title, -1)
+    except Exception:
+        return None
+
+
+def find_sandbox_faces(work_part: Any, attr_name: str = "ISOGRID_SANDBOX") -> List[Tuple[str, Any]]:
+    tagged: List[Tuple[str, Any]] = []
+    for body in getattr(work_part.Bodies, "ToArray", lambda: work_part.Bodies)():
+        for face in body.GetFaces():
+            sandbox_id = _get_string_attribute(face, attr_name)
+            if sandbox_id and sandbox_id.startswith("sandbox_"):
+                tagged.append((sandbox_id, face))
+    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] = []
+    holes: List[Dict[str, Any]] = []
+
+    loops = face.GetLoops()
+    for loop_index, loop in enumerate(loops):
+        loop_pts3d: List[Point3D] = []
+        for edge in loop.GetEdges():
+            pts = _sample_edge_polyline(edge, chord_tol_mm)
+            if loop_pts3d and pts:
+                pts = pts[1:]  # avoid duplicate joining point
+            loop_pts3d.extend(pts)
+
+        loop_pts3d = _close_polyline(loop_pts3d)
+        loop_pts2d = project_to_2d(loop_pts3d, frame)
+
+        is_outer = False
+        try:
+            is_outer = loop.IsOuter()
+        except Exception:
+            is_outer = (loop_index == 0)
+
+        if is_outer:
+            outer_2d = loop_pts2d
+            continue
+
+        is_circ, center, diameter = is_loop_circular(loop_pts2d)
+        holes.append(
+            {
+                "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())
+    except Exception:
+        pass
+
+    return geom
+
+
+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)
+    return written
+
+
+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)
+    if not sandbox_faces:
+        raise RuntimeError("No faces found with ISOGRID_SANDBOX attribute.")
+
+    payloads: Dict[str, Dict[str, Any]] = {}
+    for sandbox_id, face in sandbox_faces:
+        payloads[sandbox_id] = extract_sandbox_geometry(face, sandbox_id, chord_tol_mm=chord_tol_mm)
+
+    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.add_argument("--output-dir", default=".", help="Directory for geometry_sandbox_*.json")
+    parser.add_argument("--chord-tol", type=float, default=0.1, help="Edge sampling chord tolerance (mm)")
+    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)
+    for p in written:
+        print(f"[extract_sandbox] wrote: {p}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/tools/adaptive-isogrid/src/nx/import_profile.py b/tools/adaptive-isogrid/src/nx/import_profile.py
new file mode 100644
index 00000000..65cc5c18
--- /dev/null
+++ b/tools/adaptive-isogrid/src/nx/import_profile.py
@@ -0,0 +1,163 @@
+"""
+NXOpen script — import rib profile JSON and replace sandbox geometry.
+
+Input:
+  rib_profile_<sandbox_id>.json (or rib_profile.json)
+
+Responsibilities:
+- Recreate closed polylines from profile coordinate arrays
+- Build sheet region for sandbox
+- Replace sandbox face geometry only
+- Sew/unite with neighboring reserved faces
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+from pathlib import Path
+from typing import Any, Dict, Iterable, List, Sequence, Tuple
+
+
+Point2D = Tuple[float, float]
+Point3D = Tuple[float, float, float]
+
+
+def _add(a: Sequence[float], b: Sequence[float]) -> Point3D:
+    return (a[0] + b[0], a[1] + b[1], a[2] + b[2])
+
+
+def _mul(v: Sequence[float], s: float) -> Point3D:
+    return (v[0] * s, v[1] * s, v[2] * s)
+
+
+def load_json(path: Path) -> Dict[str, Any]:
+    return json.loads(path.read_text())
+
+
+def map_2d_to_3d(p: Point2D, transform: Dict[str, List[float]]) -> Point3D:
+    origin = transform["origin"]
+    x_axis = transform["x_axis"]
+    y_axis = transform["y_axis"]
+    return _add(_add(origin, _mul(x_axis, p[0])), _mul(y_axis, p[1]))
+
+
+def _ensure_closed(coords: List[Point2D]) -> List[Point2D]:
+    if not coords:
+        return coords
+    if coords[0] != coords[-1]:
+        coords.append(coords[0])
+    return coords
+
+
+def _create_polyline_curve(work_part: Any, pts3d: List[Point3D]) -> Any:
+    """
+    Create a closed polyline curve in NX.
+    API notes: this can be implemented with StudioSplineBuilderEx, PolygonBuilder,
+    or line segments + composite curve depending on NX version/license.
+    """
+    # Line-segment fallback (works in all NX versions)
+    curves = []
+    for i in range(len(pts3d) - 1):
+        p1 = work_part.Points.CreatePoint(pts3d[i])
+        p2 = work_part.Points.CreatePoint(pts3d[i + 1])
+        curves.append(work_part.Curves.CreateLine(p1, p2))
+    return curves
+
+
+def build_profile_curves(work_part: Any, profile: Dict[str, Any], transform: Dict[str, List[float]]) -> Dict[str, List[Any]]:
+    created: Dict[str, List[Any]] = {"outer": [], "pockets": [], "holes": []}
+
+    outer = _ensure_closed([(float(x), float(y)) for x, y in profile["outer_boundary"]])
+    outer_3d = [map_2d_to_3d(p, transform) for p in outer]
+    created["outer"] = _create_polyline_curve(work_part, outer_3d)
+
+    for pocket in profile.get("pockets", []):
+        coords = _ensure_closed([(float(x), float(y)) for x, y in pocket])
+        pts3d = [map_2d_to_3d(p, transform) for p in coords]
+        created["pockets"].extend(_create_polyline_curve(work_part, pts3d))
+
+    for hole in profile.get("hole_boundaries", []):
+        coords = _ensure_closed([(float(x), float(y)) for x, y in hole])
+        pts3d = [map_2d_to_3d(p, transform) for p in coords]
+        created["holes"].extend(_create_polyline_curve(work_part, pts3d))
+
+    return created
+
+
+def _find_sandbox_face(work_part: Any, sandbox_id: str) -> Any:
+    for body in getattr(work_part.Bodies, "ToArray", lambda: work_part.Bodies)():
+        for face in body.GetFaces():
+            try:
+                tag = face.GetStringUserAttribute("ISOGRID_SANDBOX", -1)
+            except Exception:
+                tag = None
+            if tag == sandbox_id:
+                return face
+    raise RuntimeError(f"Sandbox face not found for id={sandbox_id}")
+
+
+def replace_sandbox_face_geometry(work_part: Any, sandbox_face: Any, created_curves: Dict[str, List[Any]]) -> None:
+    """
+    Replace sandbox surface region from generated profile curves.
+
+    This operation depends on the model topology and NX license package.
+    Typical implementation:
+      1) Build bounded plane/sheet from outer and inner loops
+      2) Trim/split host face by new boundaries
+      3) Delete old sandbox patch
+      4) Sew new patch with reserved neighboring faces
+      5) Unite if multiple sheet bodies are produced
+    """
+    # Recommended implementation hook points.
+    # - Through Curve Mesh / Bounded Plane builders in NXOpen.Features
+    # - SewBuilder in NXOpen.Features
+    # - DeleteFace + ReplaceFace in synchronous modeling toolkit
+    raise NotImplementedError(
+        "Sandbox face replacement is model-specific. Implement with NXOpen feature builders "
+        "(bounded sheet + replace face + sew/unite) in target NX environment."
+    )
+
+
+def run_in_nx(
+    profile_path: Path,
+    geometry_path: Path,
+    sandbox_id: str,
+) -> None:
+    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.")
+
+    profile = load_json(profile_path)
+    geometry = load_json(geometry_path)
+    transform = geometry.get("transform")
+    if not transform:
+        raise ValueError(f"Missing transform in {geometry_path}")
+
+    sandbox_face = _find_sandbox_face(work_part, sandbox_id)
+    created_curves = build_profile_curves(work_part, profile, transform)
+    replace_sandbox_face_geometry(work_part, sandbox_face, created_curves)
+
+    print(f"[import_profile] Imported profile for {sandbox_id}: {profile_path}")
+
+
+def main(argv: Sequence[str] | None = None) -> int:
+    parser = argparse.ArgumentParser(description="Import rib profile JSON into NX sandbox face")
+    parser.add_argument("--profile", required=True, help="Path to rib_profile json")
+    parser.add_argument("--geometry", required=True, help="Path to geometry_sandbox json")
+    parser.add_argument("--sandbox-id", required=True, help="Sandbox id (e.g. sandbox_1)")
+    args = parser.parse_args(argv)
+
+    run_in_nx(
+        profile_path=Path(args.profile),
+        geometry_path=Path(args.geometry),
+        sandbox_id=args.sandbox_id,
+    )
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/tools/adaptive-isogrid/src/nx/run_iteration.py b/tools/adaptive-isogrid/src/nx/run_iteration.py
new file mode 100644
index 00000000..f6aa8278
--- /dev/null
+++ b/tools/adaptive-isogrid/src/nx/run_iteration.py
@@ -0,0 +1,120 @@
+"""
+NXOpen orchestrator — single Adaptive Isogrid iteration.
+
+Pipeline (inside NX batch or GUI session):
+  1) extract_sandbox.py
+  2) external Python Brain call
+  3) import_profile.py
+  4) solve_and_extract.py
+
+This script is designed to be launched from external Python with subprocess,
+for example via run_journal.exe / NX batch mode.
+"""
+
+from __future__ import annotations
+
+import argparse
+import subprocess
+from pathlib import Path
+from typing import List, Sequence
+
+try:
+    from .extract_sandbox import run_in_nx as run_extract_sandbox
+    from .import_profile import run_in_nx as run_import_profile
+    from .solve_and_extract import run_in_nx as run_solve_and_extract
+except ImportError:
+    # NX run_journal often executes files as scripts (no package context).
+    from extract_sandbox import run_in_nx as run_extract_sandbox  # type: ignore
+    from import_profile import run_in_nx as run_import_profile  # type: ignore
+    from solve_and_extract import run_in_nx as run_solve_and_extract  # type: ignore
+
+
+def _detect_sandbox_ids_from_files(work_dir: Path) -> List[str]:
+    ids = []
+    for p in sorted(work_dir.glob("geometry_sandbox_*.json")):
+        # geometry_sandbox_1.json -> sandbox_1
+        suffix = p.stem.replace("geometry_", "", 1)
+        ids.append(suffix)
+    return ids
+
+
+def call_python_brain(brain_cmd: str, work_dir: Path, sandbox_ids: List[str]) -> None:
+    """
+    Executes external Python Brain command.
+
+    The command should read geometry_sandbox_*.json from work_dir and write
+    rib_profile_sandbox_*.json files in the same folder.
+    """
+    env = None
+    cmd = [brain_cmd, "--work-dir", str(work_dir)]
+    for sid in sandbox_ids:
+        cmd += ["--sandbox-id", sid]
+
+    result = subprocess.run(cmd, cwd=str(work_dir), capture_output=True, text=True)
+    if result.returncode != 0:
+        raise RuntimeError(
+            "Python Brain call failed\n"
+            f"cmd: {' '.join(cmd)}\n"
+            f"stdout:\n{result.stdout}\n"
+            f"stderr:\n{result.stderr}"
+        )
+
+
+def run_iteration(
+    work_dir: Path,
+    brain_cmd: str,
+    solution_name: str | None = None,
+    csv_fallback: bool = False,
+) -> Path:
+    work_dir.mkdir(parents=True, exist_ok=True)
+
+    # 1) Extract all sandbox geometries from active NX model.
+    run_extract_sandbox(output_dir=work_dir, chord_tol_mm=0.1)
+    sandbox_ids = _detect_sandbox_ids_from_files(work_dir)
+    if not sandbox_ids:
+        raise RuntimeError("No geometry_sandbox_*.json files generated.")
+
+    # 2) Call external Brain.
+    call_python_brain(brain_cmd=brain_cmd, work_dir=work_dir, sandbox_ids=sandbox_ids)
+
+    # 3) Import each rib profile back into NX.
+    for sid in sandbox_ids:
+        geom_path = work_dir / f"geometry_{sid}.json"
+        profile_path = work_dir / f"rib_profile_{sid}.json"
+        if not profile_path.exists():
+            # fallback name used by some scripts
+            profile_path = work_dir / "rib_profile.json"
+        run_import_profile(profile_path=profile_path, geometry_path=geom_path, sandbox_id=sid)
+
+    # 4) Remesh + solve + extract results.
+    results_path = work_dir / "results.json"
+    run_solve_and_extract(
+        work_dir=work_dir,
+        result_path=results_path,
+        sandbox_ids=sandbox_ids,
+        solution_name=solution_name,
+        use_csv_fallback=csv_fallback,
+    )
+    return results_path
+
+
+def main(argv: Sequence[str] | None = None) -> int:
+    parser = argparse.ArgumentParser(description="Run one adaptive-isogrid NX iteration")
+    parser.add_argument("--work-dir", required=True, help="Iteration working directory")
+    parser.add_argument("--brain-cmd", required=True, help="Executable/command for external Python Brain")
+    parser.add_argument("--solution", default=None, help="Optional NX solution name")
+    parser.add_argument("--csv-fallback", action="store_true", help="Use CSV fallback for results extraction")
+    args = parser.parse_args(argv)
+
+    results_path = run_iteration(
+        work_dir=Path(args.work_dir),
+        brain_cmd=args.brain_cmd,
+        solution_name=args.solution,
+        csv_fallback=args.csv_fallback,
+    )
+    print(f"[run_iteration] completed -> {results_path}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/tools/adaptive-isogrid/src/nx/solve_and_extract.py b/tools/adaptive-isogrid/src/nx/solve_and_extract.py
new file mode 100644
index 00000000..46fd36c8
--- /dev/null
+++ b/tools/adaptive-isogrid/src/nx/solve_and_extract.py
@@ -0,0 +1,166 @@
+"""
+NXOpen script — remesh full plate, solve, and export results.json.
+
+Outputs:
+{
+  "nodes_xy": [[x, y], ...],
+  "stress_values": [...],
+  "disp_values": [...],
+  "strain_values": [...],
+  "mass": 0.0
+}
+"""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import json
+from pathlib import Path
+from typing import Any, Dict, List, Sequence, Tuple
+
+
+def remesh_full_plate(sim_part: Any) -> None:
+    """
+    Trigger mesh regeneration for all FEM collectors associated with the plate.
+    """
+    # This hook should iterate mesh managers/collectors in your SIM/FEM template.
+    # API often used: FEModel.UpdateFemodel(), mesh manager GenerateMesh/Update.
+    fe_model = sim_part.FindObject("FEModel")
+    fe_model.UpdateFemodel()
+
+
+def solve_active_solution(sim_part: Any, solution_name: str | None = None) -> Any:
+    """
+    Solve the requested solution (or first available).
+    """
+    import NXOpen  # type: ignore
+    import NXOpen.CAE  # type: ignore
+
+    simulation = sim_part.FindObject("Simulation")
+
+    target_solution = None
+    if solution_name:
+        try:
+            target_solution = simulation.FindObject(f"Solution[{solution_name}]")
+        except Exception:
+            target_solution = None
+
+    if target_solution is None:
+        target_solution = simulation.FindObject("Solution[Solution 1]")
+
+    solve_mgr = NXOpen.CAE.SimSolveManager.GetSimSolveManager(NXOpen.Session.GetSession())
+    solve_mgr.SubmitSolves([target_solution])
+    return target_solution
+
+
+def _parse_csv_results(path: Path) -> Tuple[List[List[float]], List[float]]:
+    """
+    Parse generic CSV with columns: x,y,value (header-insensitive).
+    """
+    coords: List[List[float]] = []
+    values: List[float] = []
+    with path.open(newline="") as f:
+        reader = csv.DictReader(f)
+        for row in reader:
+            keys = {k.lower(): k for k in row.keys()}
+            x = float(row[keys.get("x", keys.get("x_mm", "x"))])
+            y = float(row[keys.get("y", keys.get("y_mm", "y"))])
+            v = float(row[keys.get("value", keys.get("von_mises", "value"))])
+            coords.append([x, y])
+            values.append(v)
+    return coords, values
+
+
+def extract_results_nxopen(sim_part: Any, sandbox_ids: List[str]) -> Dict[str, Any]:
+    """
+    Preferred extractor: NXOpen post-processing API.
+
+    NOTE: Result object names/components vary by template. Keep this function as the
+    primary integration point for project-specific API wiring.
+    """
+    raise NotImplementedError(
+        "Wire NXOpen post-processing calls here (nodal stress/displacement, elemental strain, mass)."
+    )
+
+
+def extract_results_csv_fallback(work_dir: Path) -> Dict[str, Any]:
+    """
+    Fallback extractor: parse Simcenter-exported CSV files in work_dir.
+
+    Expected files:
+      - nodal_stress.csv
+      - nodal_disp.csv
+      - elemental_strain.csv
+      - mass.json (optional: {"mass": ...})
+    """
+    stress_coords, stress_vals = _parse_csv_results(work_dir / "nodal_stress.csv")
+    disp_coords, disp_vals = _parse_csv_results(work_dir / "nodal_disp.csv")
+    _, strain_vals = _parse_csv_results(work_dir / "elemental_strain.csv")
+
+    # Use stress nodal coordinates as canonical nodes_xy
+    nodes_xy = stress_coords if stress_coords else disp_coords
+
+    mass = 0.0
+    mass_file = work_dir / "mass.json"
+    if mass_file.exists():
+        mass = float(json.loads(mass_file.read_text()).get("mass", 0.0))
+
+    return {
+        "nodes_xy": nodes_xy,
+        "stress_values": stress_vals,
+        "disp_values": disp_vals,
+        "strain_values": strain_vals,
+        "mass": mass,
+    }
+
+
+def run_in_nx(
+    work_dir: Path,
+    result_path: Path,
+    sandbox_ids: List[str],
+    solution_name: str | None = None,
+    use_csv_fallback: bool = False,
+) -> Dict[str, Any]:
+    import NXOpen  # type: ignore
+
+    session = NXOpen.Session.GetSession()
+    sim_part = session.Parts.BaseWork
+    if sim_part is None:
+        raise RuntimeError("No active NX SIM/FEM work part.")
+
+    remesh_full_plate(sim_part)
+    solve_active_solution(sim_part, solution_name=solution_name)
+
+    if use_csv_fallback:
+        results = extract_results_csv_fallback(work_dir)
+    else:
+        results = extract_results_nxopen(sim_part, sandbox_ids)
+
+    result_path.parent.mkdir(parents=True, exist_ok=True)
+    result_path.write_text(json.dumps(results, indent=2))
+    print(f"[solve_and_extract] wrote {result_path}")
+    return results
+
+
+def main(argv: Sequence[str] | None = None) -> int:
+    parser = argparse.ArgumentParser(description="Remesh, solve, and export results.json")
+    parser.add_argument("--work-dir", default=".", help="Working directory for CSV fallback artifacts")
+    parser.add_argument("--results", default="results.json", help="Output JSON path")
+    parser.add_argument("--sandbox-id", action="append", default=[], help="Sandbox id filter (repeatable)")
+    parser.add_argument("--solution", default=None, help="NX solution name (defaults to Solution 1)")
+    parser.add_argument("--csv-fallback", action="store_true", help="Parse CSV files instead of NXOpen post API")
+    args = parser.parse_args(argv)
+
+    run_in_nx(
+        work_dir=Path(args.work_dir),
+        result_path=Path(args.results),
+        sandbox_ids=args.sandbox_id,
+        solution_name=args.solution,
+        use_csv_fallback=args.csv_fallback,
+    )
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())