From bb61595c1255888954699de079bd5a8857846d3c Mon Sep 17 00:00:00 2001
From: Miguel de la Varga <miguel@terranigma-solutions.com>
Date: Sat, 19 Sep 2026 11:41:38 +0200
Subject: [PATCH 1/2] [ENH] Add dual-contouring extent capping implementation
 with tests and documentation

Introduce a dual-contouring extent capping feature for mesh closure, including scalar-based boundary handling and outward-facing triangulation. Provides geometry, performance diagnostics, and failure reporting via `capping_report`. Add tests validating closure success, geometry correctness, and boundary lattice properties. Include detailed documentation for usage, best-effort constraints, and configuration options.
---
 docs/extent_capping.md                        | 109 +++++
 .../multi_scalar_dual_contouring.py           |  79 +++-
 gempy_engine/API/model/model_api.py           |  20 +
 gempy_engine/core/data/__init__.py            |   2 +-
 .../core/data/dual_contouring_data.py         |   1 +
 .../core/data/dual_contouring_mesh.py         |   6 +
 gempy_engine/core/data/engine_grid.py         |   4 +-
 .../core/data/options/evaluation_options.py   |   6 +
 gempy_engine/core/data/raw_arrays_solution.py |   3 +-
 gempy_engine/core/data/regular_grid.py        |   3 +
 .../dual_contouring/_extent_capping.py        | 368 ++++++++++++++++
 .../modules/dual_contouring/_gen_vertices.py  |  13 +-
 .../dual_contouring/_scalar_crossing.py       |  32 ++
 .../dual_contouring/_support_report.py        |   7 +-
 .../dual_contouring_interface.py              |  34 +-
 .../test_modules/test_extent_capping.py       | 342 +++++++++++++++
 .../test_extent_capping_integration.py        | 412 ++++++++++++++++++
 .../test_modules/test_scalar_crossing.py      | 111 +++++
 18 files changed, 1539 insertions(+), 13 deletions(-)
 create mode 100644 docs/extent_capping.md
 create mode 100644 gempy_engine/modules/dual_contouring/_extent_capping.py
 create mode 100644 gempy_engine/modules/dual_contouring/_scalar_crossing.py
 create mode 100644 tests/test_common/test_modules/test_extent_capping.py
 create mode 100644 tests/test_common/test_modules/test_extent_capping_integration.py
 create mode 100644 tests/test_common/test_modules/test_scalar_crossing.py

diff --git a/docs/extent_capping.md b/docs/extent_capping.md
new file mode 100644
index 00000000..4e948132
--- /dev/null
+++ b/docs/extent_capping.md
@@ -0,0 +1,109 @@
+# Best-effort extent capping
+
+Extent capping is opt-in and operates independently on each exported isosurface:
+
+```python
+from gempy_engine.core.data import MeshExtentCapping
+
+options.evaluation_options.mesh_extraction_extent_capping = (
+    MeshExtentCapping.SCALAR_LESS_EQUAL
+)
+```
+
+The equivalent string is `"scalar_less_equal"`; the default is `"none"`.
+The inside convention is `scalar <= isovalue`. Nested surfaces therefore produce
+overlapping solids, not a partition into geological unit volumes.
+
+## Geometry
+
+The capper constructs three connected pieces:
+
+1. The ordinary interior dual-contouring triangles.
+2. Transition triangles from boundary-cell QEF vertices to the extent contour,
+   including primal-edge triangles between adjacent boundary cells.
+3. Flat, outward-facing triangles covering inside regions of the six box faces.
+
+Boundary squares use a consistent diagonal and piecewise-linear triangle
+clipping. Lattice points and edge intersections have shared integer identities,
+including intersections on face diagonals and on box edges and corners.
+
+Only cap components connected to a contour with an available, contained QEF
+vertex are retained. An enclosed isosurface receives no disconnected box shell.
+If no isosurface intersects the box, no box shell is manufactured, even when
+the whole box satisfies the inside condition. This is closure of existing
+isosurfaces, not unconditional extraction of the complete clipped sublevel set.
+
+`RegularGrid.physical_extent` retains the requested extent through refinement.
+Enabled `compute_model` calls use a private root grid sampled against that exact
+box, rather than the legacy `1e-6` translated box. Caps are generated in engine
+coordinates, before downstream output transforms. Original vertex indices are
+preserved and cap vertices are appended.
+
+Enabled extraction and cap clipping share strict scalar crossing rules. An
+endpoint equal to the isovalue is inside; an entire iso-valued edge is not a
+crossing. Interpolation parameters are clamped and snapped within `1e-12` of an
+endpoint. Nonfinite scalar inputs are rejected with a diagnostic. Enabled QEF
+mass points use edge-validity masks, including genuine zero coordinates, and
+avoid PyTorch's additional origin-centered regularization.
+
+With capping disabled, legacy grid sampling, crossings, QEF solving, and mesh
+arrays retain their previous behavior. Enabled interior geometry can differ
+because its sampling and crossing conventions are deliberately stricter.
+
+## Best-Effort Contract
+
+Inspect `mesh.capping_report["closure_success"]` before treating a result as a
+closed solid. Unsupported geometry emits `RuntimeWarning` and returns the
+available mesh; it is not silently declared successful.
+
+- Fault stacks, fault-affected stacks, and stacks with removed extraction-mask
+  cells currently skip capping. Scalar samples alone do not encode enough cut
+  ownership to safely prevent bridging intentional fault or relation boundaries.
+  These meshes retain their extracted geometry and report `skipped_reason`.
+- Missing refinement support is never repaired by extending a cap into the
+  interior. Missing boundary QEFs are reported and transitions are omitted.
+- Escaped QEF vertices are reported and excluded from transition construction;
+  this implementation does not introduce a bounded QEF solver.
+- Exact ties and multiple contour components in one cell remain best effort.
+  Edge-incidence, directed-winding, and vertex-link audits detect many failures,
+  but there is no geometric self-intersection test.
+- New duplicate or zero-area triangles are removed. Invalid original geometry
+  is preserved and reported unsuccessful rather than silently repaired.
+- `GEMPY_SKIP_TRIANGULATION` also skips capping and boundary evaluation.
+
+The full issue's general fault/mask-aware capping and all ambiguous topology
+cases are not implemented by this first version.
+
+## Reports and Metadata
+
+Each mesh exposes `stack_index`, `surface_index`, `exported_surface_index`,
+`isovalue`, and, when enabled, `inside_convention`.
+
+The capping report includes before/after open-edge counts, inferred non-extent
+openings, added cap and transition triangle counts, removed triangle counts,
+missing/escaped QEF counts, manifoldness and winding checks, maximum cap-plane
+error, evaluated boundary-point count, and the original vertex count.
+
+`watertight` describes topology; `closure_success` additionally checks finite,
+nondegenerate geometry and unresolved construction diagnostics. Neither proves
+absence of self-intersections. Physical-opening classification for original
+QEF edges uses owning-cell membership, not the interior QEF positions; it is
+therefore a cell-level inference rather than proof of the opening's cause.
+Cap-plane checks use `64 * float64_epsilon * max(1, abs(extent))` tolerance.
+
+`vertices_tensor` remains the pre-overlap differentiable QEF snapshot. Final
+triangle indices reference `mesh.vertices`, **not** `vertices_tensor`. Appended
+cap geometry and topology are NumPy postprocessing and are not differentiable.
+
+## Cost
+
+One unique finest-level six-face lattice is shared across surfaces. Boundary
+scalar values are evaluated once per stack per batch and reused across that
+stack's isosurfaces. `evaluation_chunk_size` bounds the boundary point batch;
+the existing evaluator also applies its own kernel-workload chunking.
+
+Boundary storage scales as `O(nx*ny + nx*nz + ny*nz)`, but Python dictionaries,
+triangle connectivity, and topology audits have substantial additional memory
+cost. This is a correctness-first implementation: there is no adaptive cap
+simplification, cross-call cache, or reuse of shifted legacy corner samples.
+Large-depth GPU and production-scale memory benchmarks remain outstanding.
diff --git a/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py b/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py
index 3a357871..05ddc042 100644
--- a/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py
+++ b/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py
@@ -1,4 +1,5 @@
 import copy
+import os
 import warnings
 from typing import List, Any
 
@@ -24,7 +25,8 @@
                                                                   get_masked_codes, mask_generation)
 from ...modules.dual_contouring.overlapping import average_overlapping_vertices, remove_fault_overlap_triangles
 from ...modules.dual_contouring._support_report import mesh_support_report
-from ...core.data.options.evaluation_options import OctreeRefinementMode
+from ...core.data.options.evaluation_options import OctreeRefinementMode, MeshExtentCapping
+from ...modules.dual_contouring._extent_capping import boundary_lattice, cap_mesh
 
 
 @gempy_profiler_decorator
@@ -52,6 +54,7 @@ def dual_contouring_multi_scalar(
 
     octree_leaves = octree_list[-1]
     all_meshes: List[DualContouringMesh] = []
+    cap_enabled = MeshExtentCapping(options.evaluation_options.mesh_extraction_extent_capping) != MeshExtentCapping.NONE
 
     dual_contouring_options = copy.deepcopy(options)
     dual_contouring_options.evaluation_options.compute_scalar_gradient = True
@@ -86,7 +89,8 @@ def dual_contouring_multi_scalar(
             _xyz_corners=octree_leaves.grid.corners_grid.values,
             scalar_field_on_corners=output.exported_fields.scalar_field[output.grid.corners_grid_slice],
             scalar_at_sp=output.scalar_field_at_sp,
-            masking=mask
+            masking=mask,
+            strict_crossings=cap_enabled
         )
 
         all_surfaces_intersection.append(intersection_xyz)
@@ -110,6 +114,7 @@ def dual_contouring_multi_scalar(
     # Generate meshes for each scalar field
     dc_data_per_surface_all = []
     support_reports = []
+    surface_metadata = []
     stack_relations = data_descriptor.stack_structure.masking_descriptor
     for n_scalar_field in range(data_descriptor.stack_structure.n_stacks):
         if stack_relations[n_scalar_field] is StackRelationType.NULL_SPACE:
@@ -125,7 +130,8 @@ def dual_contouring_multi_scalar(
                     output.exported_fields.scalar_field[output.grid.corners_grid_slice],
                     output.scalar_field_at_sp[surface_i], base_number, mask,
                     surface_index=surface_i,
-                    ancestor_coordinates=[level.grid.octree_grid.integer_coordinates for level in octree_list[:-1]]
+                    ancestor_coordinates=[level.grid.octree_grid.integer_coordinates for level in octree_list[:-1]],
+                    strict_crossings=cap_enabled
                 )
                 report['stack_index'] = n_scalar_field
                 if report['internal_refinement_boundary_edge_count']:
@@ -154,11 +160,13 @@ def dual_contouring_multi_scalar(
                 tree_depth=options.number_octree_levels,
                 base_number=base_number,
                 triangulation_method=options.evaluation_options.triangulation_method,
-                generated_cell_coordinates=left_right_codes
+                generated_cell_coordinates=left_right_codes,
+                strict_crossings=cap_enabled
             )
 
             dc_data_per_surface_all.append(dc_data_per_surface)
             surface_to_stack.append(n_scalar_field)
+            surface_metadata.append((n_scalar_field, surface_i, float(output.scalar_field_at_sp[surface_i])))
             if compute_overlap:
                 left_right_per_mesh.append(all_left_right_codes[n_scalar_field][dc_data_per_surface.valid_voxels])
 
@@ -213,9 +221,72 @@ def dual_contouring_multi_scalar(
             mesh.vertices = BackendTensor.t.to_numpy(mesh.vertices)
             mesh.edges = BackendTensor.t.to_numpy(mesh.edges)
 
+    for index, (mesh, (stack_index, surface_index, isovalue)) in enumerate(zip(all_meshes, surface_metadata)):
+        mesh.stack_index = stack_index
+        mesh.surface_index = surface_index
+        mesh.exported_surface_index = index
+        mesh.isovalue = isovalue
+        mesh.inside_convention = "scalar <= isovalue" if cap_enabled else None
+
+    skip_triangles = os.getenv("GEMPY_SKIP_TRIANGULATION", "0").lower() in ("true", "1", "t", "y", "yes")
+    if cap_enabled and all_meshes and not skip_triangles:
+        extent = BackendTensor.t.to_numpy(octree_list[0].grid.octree_grid.physical_extent)
+        skip_reasons = []
+        fault_relations = data_descriptor.stack_structure.faults_relations
+        for mesh in all_meshes:
+            stack = mesh.stack_index
+            mask = all_mask_arrays[stack]
+            faulted = (stack_relations[stack] is StackRelationType.FAULT
+                       or (fault_relations is not None and np.any(fault_relations[:, stack])))
+            masked = mask is not None and not bool(mask.all())
+            skip_reasons.append(
+                "Capping skipped: fault or extraction-mask boundary ownership is not supported"
+                if faulted or masked else None
+            )
+        geometry = None
+        coordinates = points = np.empty((0, 3))
+        boundary_scalars = None
+        if any(reason is None for reason in skip_reasons):
+            geometry = boundary_lattice(base_number, extent)
+            coordinates, points, _ = geometry
+            boundary_scalars = _interp_on_boundary(points, interpolation_input, options, data_descriptor)
+        for mesh, dc_data, reason in zip(all_meshes, dc_data_per_surface_all, skip_reasons):
+            cell_coordinates = BackendTensor.t.to_numpy(dc_data.left_right_codes[dc_data.valid_voxels])
+            cap_mesh(mesh, cell_coordinates, base_number, extent, coordinates,
+                     boundary_scalars[mesh.stack_index] if reason is None else np.empty(0),
+                     mesh.isovalue, boundary_geometry=geometry, skip_reason=reason)
+            mesh.capping_report["boundary_scalar_points"] = len(points) if reason is None else 0
+            mesh.capping_report["original_vertex_count"] = len(cell_coordinates)
+
     return all_meshes
 
 
+def _interp_on_boundary(points, interpolation_input, options, data_descriptor):
+    """Evaluate each stack once per boundary batch, reusing it for its surfaces."""
+    boundary_options = copy.deepcopy(options)
+    boundary_options.evaluation_options.compute_scalar = True
+    boundary_options.evaluation_options.compute_scalar_gradient = False
+    saved_grid = interpolation_input.grid
+    saved_stack = data_descriptor.stack_structure.stack_number
+    scalars = [np.empty(len(points)) for _ in range(data_descriptor.stack_structure.n_stacks)]
+    batch_size = max(1, int(options.evaluation_options.evaluation_chunk_size))
+    try:
+        for start in range(0, len(points), batch_size):
+            stop = min(start + batch_size, len(points))
+            interpolation_input.set_temp_grid(EngineGrid(custom_grid=GenericGrid(
+                values=BackendTensor.t.array(points[start:stop], dtype=BackendTensor.dtype)
+            )))
+            outputs = interpolate_all_fields_no_octree(interpolation_input, boundary_options, data_descriptor)
+            for stack_index, output in enumerate(outputs):
+                scalars[stack_index][start:stop] = BackendTensor.t.to_numpy(
+                    output.exported_fields.scalar_field[output.grid.custom_grid_slice]
+                )
+    finally:
+        interpolation_input.set_temp_grid(saved_grid)
+        data_descriptor.stack_structure.stack_number = saved_stack
+    return scalars
+
+
 def _validate_stack_relations(data_descriptor: InputDataDescriptor, n_scalar_field: int) -> None:
     """
     Validate stack relations for the given scalar field.
diff --git a/gempy_engine/API/model/model_api.py b/gempy_engine/API/model/model_api.py
index 89cb670a..033bb4ab 100644
--- a/gempy_engine/API/model/model_api.py
+++ b/gempy_engine/API/model/model_api.py
@@ -20,6 +20,7 @@
 from ...core.utils import gempy_profiler_decorator
 from ...core.exceptions import GemPyEngineInputError
 from ...core.data.options.temp_interpolation_values import TempInterpolationValues
+from ...core.data.options.evaluation_options import MeshExtentCapping
 from ...modules.geophysics.fw_gravity import compute_gravity
 from ...modules.geophysics.fw_magnetic import compute_tmi
 from ...modules.weights_cache.weights_cache_interface import WeightCache
@@ -43,6 +44,25 @@ def compute_model(interpolation_input: InterpolationInput, options: Interpolatio
         # Check input is valid
         _check_input_validity(interpolation_input, options, data_descriptor)  # TODO
 
+        if (options.evaluation_options.mesh_extraction
+                and MeshExtentCapping(options.evaluation_options.mesh_extraction_extent_capping) != MeshExtentCapping.NONE
+                and interpolation_input.grid.octree_grid is not None):
+            # Keep the legacy offset and caller-owned grids untouched when capping
+            # is disabled. Enabled extraction samples the exact physical box.
+            interpolation_input = copy.copy(interpolation_input)
+            grid = copy.copy(interpolation_input.grid)
+            root = copy.copy(grid.octree_grid)
+            root.orthogonal_extent = BackendTensor.t.copy(root.physical_extent)
+            coordinates = BackendTensor.t.array(root.integer_coordinates, dtype=BackendTensor.dtype)
+            root.values = root.physical_extent[::2] + (coordinates + 0.5) * (
+                (root.physical_extent[1::2] - root.physical_extent[::2]) / root.regular_grid_shape
+            )
+            root.original_values = BackendTensor.t.copy(root.values)
+            grid.octree_grid = root
+            grid.corners_grid = None
+            interpolation_input._original_grid = grid
+            interpolation_input.set_temp_grid(grid)
+
         output: list[OctreeLevel] = interpolate_n_octree_levels(
             interpolation_input=interpolation_input,
             options=options,
diff --git a/gempy_engine/core/data/__init__.py b/gempy_engine/core/data/__init__.py
index 2715672e..acc87081 100644
--- a/gempy_engine/core/data/__init__.py
+++ b/gempy_engine/core/data/__init__.py
@@ -3,5 +3,5 @@
 from .kernel_classes.orientations import Orientations, OrientationsInternals
 from .kernel_classes.surface_points import SurfacePoints, SurfacePointsInternals
 from .options.interpolation_options import InterpolationOptions
-from .options.evaluation_options import OctreeRefinementMode
+from .options.evaluation_options import OctreeRefinementMode, MeshExtentCapping
 from .solutions import Solutions
diff --git a/gempy_engine/core/data/dual_contouring_data.py b/gempy_engine/core/data/dual_contouring_data.py
index b43f22cd..b794e422 100644
--- a/gempy_engine/core/data/dual_contouring_data.py
+++ b/gempy_engine/core/data/dual_contouring_data.py
@@ -32,6 +32,7 @@ class DualContouringData:
     triangulation_method: TriangulationMethod = TriangulationMethod.LEGACY
     generated_cell_coordinates: Optional[np.ndarray] = None  # Before geological masking.
     triangulation_report: dict = field(default_factory=dict)  # Quad support before overlap/fault triangle removal.
+    strict_crossings: bool = False
 
     @property
     def valid_voxels(self):
diff --git a/gempy_engine/core/data/dual_contouring_mesh.py b/gempy_engine/core/data/dual_contouring_mesh.py
index 692ad997..25a6108e 100644
--- a/gempy_engine/core/data/dual_contouring_mesh.py
+++ b/gempy_engine/core/data/dual_contouring_mesh.py
@@ -11,6 +11,12 @@ class DualContouringMesh:
     edges: np.ndarray
     dc_data: Optional[DualContouringData] = None  # * In principle we need this just for testing
     support_report: Optional[dict] = None
+    capping_report: Optional[dict] = None
+    stack_index: Optional[int] = None
+    surface_index: Optional[int] = None
+    exported_surface_index: Optional[int] = None
+    isovalue: Optional[float] = None
+    inside_convention: Optional[str] = None
 
     def __repr__(self):
         return f"DualContouringMesh({self.vertices.shape[0]} vertices, {self.edges.shape[0]} edges)"
diff --git a/gempy_engine/core/data/engine_grid.py b/gempy_engine/core/data/engine_grid.py
index 1806ccc5..8b42168b 100644
--- a/gempy_engine/core/data/engine_grid.py
+++ b/gempy_engine/core/data/engine_grid.py
@@ -59,10 +59,12 @@ def from_xyz_coords(cls, xyz_coords: ndarray) -> "EngineGrid":
 
     @classmethod
     def from_regular_grid(cls, regular_grid: RegularGrid) -> "EngineGrid":
-        return cls(
+        grid = cls(
             dense_grid=regular_grid,
             octree_grid=RegularGrid(regular_grid.orthogonal_extent, np.array([2, 2, 2]))
         )
+        grid.octree_grid.physical_extent = BackendTensor.t.copy(regular_grid.physical_extent)
+        return grid
 
     @property
     def values(self) -> np.ndarray:
diff --git a/gempy_engine/core/data/options/evaluation_options.py b/gempy_engine/core/data/options/evaluation_options.py
index c9aad945..6ed35f60 100644
--- a/gempy_engine/core/data/options/evaluation_options.py
+++ b/gempy_engine/core/data/options/evaluation_options.py
@@ -15,6 +15,11 @@ class OctreeRefinementMode(str, enum.Enum):
     PRECISE = "precise"
 
 
+class MeshExtentCapping(str, enum.Enum):
+    NONE = "none"
+    SCALAR_LESS_EQUAL = "scalar_less_equal"
+
+
 class MeshExtractionMaskingOptions(enum.Enum):
     NOTHING = enum.auto()  # * This is only for testing
     DISJOINT = enum.auto()
@@ -39,6 +44,7 @@ class EvaluationOptions:
     triangulation_method: TriangulationMethod = TriangulationMethod.LEGACY
     
     mesh_extraction: bool = True
+    mesh_extraction_extent_capping: MeshExtentCapping = MeshExtentCapping.NONE
     mesh_extraction_masking_options: MeshExtractionMaskingOptions = MeshExtractionMaskingOptions.INTERSECT
     mesh_extraction_fancy: Annotated[bool, deprecated("Old extraction method not in use anymore")] = True
 
diff --git a/gempy_engine/core/data/raw_arrays_solution.py b/gempy_engine/core/data/raw_arrays_solution.py
index b7f7706c..ef7f31ce 100644
--- a/gempy_engine/core/data/raw_arrays_solution.py
+++ b/gempy_engine/core/data/raw_arrays_solution.py
@@ -123,7 +123,7 @@ def meshes_to_subsurface(self, input_transform: Transform | None = None, element
         pd = require_pandas()
 
         vertex: list[np.ndarray] = self.vertices
-        simplex_list: list[np.ndarray] = self.edges
+        simplex_list: list[np.ndarray] = [triangles.copy() for triangles in self.edges]
         
         idx_max = 0
         for i, simplex_array in enumerate(simplex_list):
@@ -287,4 +287,3 @@ def _fill_finite_fault_scalar_fields_with_dense_grid(
     ]
     raw_arrays_solution.finite_fault_scalar_field_matrix = np.vstack(fields)
     raw_arrays_solution.finite_fault_stack_indices = stack_indices
-
diff --git a/gempy_engine/core/data/regular_grid.py b/gempy_engine/core/data/regular_grid.py
index 80cabd71..d2f7bb09 100644
--- a/gempy_engine/core/data/regular_grid.py
+++ b/gempy_engine/core/data/regular_grid.py
@@ -18,6 +18,7 @@ class RegularGrid:
     left_right: np.ndarray = field(default=None, repr=False, init=False)
     _integer_coordinates: np.ndarray = field(default=None, repr=False, init=False)
     refinement_debug: dict = field(default=None, repr=False, init=False)
+    physical_extent: np.ndarray = field(default=None, repr=False, init=False)
     
     values: np.ndarray = field(default=None, repr=False, init=False)
     original_values: np.ndarray = field(default=None, repr=False, init=False)  #: When the regular grid is representing a octree level, only active cells are stored in values. This is the original values of the regular grid.
@@ -28,6 +29,7 @@ def __len__(self):
 
     def __post_init__(self):
         self.regular_grid_shape = BackendTensor.t.array(self.regular_grid_shape)
+        self.physical_extent = BackendTensor.t.array(self.orthogonal_extent)
         self.orthogonal_extent = BackendTensor.t.array(self.orthogonal_extent) + 1e-6  # * This to avoid some errors evaluating in 0 (e.g. bias in dual contouring)
 
         self._create_regular_grid_3d()
@@ -68,6 +70,7 @@ def from_octree_level(cls, xyz_coords_octree: np.ndarray, previous_regular_grid:
         )
 
         regular_grid_for_octree_level.values = xyz_coords_octree  # ! Overwrite the common values
+        regular_grid_for_octree_level.physical_extent = BackendTensor.t.copy(previous_regular_grid.physical_extent)
         regular_grid_for_octree_level._active_cells = active_cells
         regular_grid_for_octree_level.left_right = left_right
         regular_grid_for_octree_level._integer_coordinates = (
diff --git a/gempy_engine/modules/dual_contouring/_extent_capping.py b/gempy_engine/modules/dual_contouring/_extent_capping.py
new file mode 100644
index 00000000..96e94211
--- /dev/null
+++ b/gempy_engine/modules/dual_contouring/_extent_capping.py
@@ -0,0 +1,368 @@
+"""Best-effort NumPy caps for a uniform dual-contouring lattice.
+
+This is not a mesh repair operation: only the six supplied extent faces are
+sampled. Internal mask/fault holes remain untouched. Ambiguous/tied contours,
+missing cells and unconstrained QEFs can leave holes; inspect ``capping_report``.
+"""
+
+from collections import defaultdict
+from itertools import product
+import warnings
+
+import numpy as np
+
+from ._scalar_crossing import scalar_crossing_parameters
+
+
+def boundary_lattice(shape, extent):
+    """Return integer corners, physical positions, and face metadata.
+
+    ``shape`` is the three *cell* counts. ``extent`` is
+    ``(xmin, xmax, ymin, ymax, zmin, zmax)``. Corners are unique globally,
+    including box edges/corners. Metadata contains ``triangles`` (indices into
+    these arrays) and ``cells`` (one owning in-domain cell per triangle).
+    Faces are triangulated with a consistent low-to-high diagonal and outward
+    winding. Positions use the exact supplied bounds, not QEF-derived bounds.
+
+    Storage scales with the full boundary area, not the active contour. The
+    Python lists/dictionaries used during construction have substantial peak
+    memory overhead on large grids. Reuse this tuple via ``boundary_geometry``
+    rather than constructing it twice; this is not a streaming implementation.
+    """
+    shape = np.asarray(shape)
+    bounds = np.asarray(extent, dtype=float).reshape(3, 2)
+    if shape.shape != (3,) or np.any(shape < 1) or np.any(shape != shape.astype(int)):
+        raise ValueError("shape must contain three positive integer cell counts")
+    if not np.isfinite(bounds).all() or np.any(bounds[:, 1] <= bounds[:, 0]):
+        raise ValueError("extent must contain finite increasing bounds")
+    shape = shape.astype(int)
+    coordinates, triangles, cells, lookup = [], [], [], {}
+    for axis in range(3):
+        u, v = (axis + 1) % 3, (axis + 2) % 3
+        for side in (0, 1):
+            for i, j in product(range(shape[u]), range(shape[v])):
+                corners = []
+                for du, dv in ((0, 0), (1, 0), (1, 1), (0, 1)):
+                    c = [0, 0, 0]
+                    c[axis], c[u], c[v] = side * shape[axis], i + du, j + dv
+                    key = tuple(c)
+                    if key not in lookup:
+                        lookup[key] = len(coordinates)
+                        coordinates.append(c)
+                    corners.append(lookup[key])
+                cell = [0, 0, 0]
+                cell[axis], cell[u], cell[v] = side * (shape[axis] - 1), i, j
+                for a, b, c in ((0, 1, 2), (0, 2, 3)):
+                    triangles.append([corners[a], corners[b if side else c], corners[c if side else b]])
+                    cells.append(cell.copy())
+    coordinates = np.asarray(coordinates, dtype=np.int64)
+    positions = bounds[:, 0] + coordinates / shape * (bounds[:, 1] - bounds[:, 0])
+    for axis in range(3):
+        positions[coordinates[:, axis] == shape[axis], axis] = bounds[axis, 1]
+    return coordinates, positions, {
+        "triangles": np.asarray(triangles, dtype=np.int64),
+        "cells": np.asarray(cells, dtype=np.int64),
+    }
+
+
+def cap_mesh(mesh, cell_coordinates, shape, extent, boundary_coordinates,
+             boundary_scalar, isovalue, *, boundary_geometry=None, skip_reason=None):
+    """Mutate NumPy ``mesh.vertices``/``mesh.edges`` and return ``mesh``.
+
+    ``edges`` means triangle indices (M, 3). ``cell_coordinates`` is (N, 3)
+    in the exact row order of the original QEF vertices, not sorted order.
+    Boundary samples may be permuted, but must cover ``boundary_lattice``.
+    The retained solid is ``scalar <= isovalue``. Original vertices/triangles
+    keep their indices and winding. ``mesh.capping_report`` records counts and
+    whole-mesh topology audits; a RuntimeWarning summarizes limitations.
+    Invalid API array shapes raise ValueError; unsupported geometry does not.
+    ``boundary_geometry`` may be the unmodified tuple from ``boundary_lattice``
+    for this same shape/extent, avoiding a second lattice construction.
+
+    Report ``open_edges_before/after`` counts incidence-one edges.
+    ``non_extent_open_edges_before/after`` excludes edges whose endpoints share
+    an extent face: original QEF membership comes from owning cells, whereas
+    new cap points use physical positions. This is a cell-level classification,
+    not proof that an opening in a boundary cell was caused by the extent.
+    ``watertight`` audits topology only; ``closure_success`` additionally
+    requires valid finite geometry and no missing/escaped/unsupported cases.
+    ``added_cap_triangles`` and ``added_transition_triangles`` count retained
+    triangles after cleanup. ``cap_max_plane_error`` is the maximum, over cap
+    triangles, of the distance to their nearest common extent plane (physical
+    units, zero if no cap). Invalid original geometry is preserved and reported
+    unsuccessful without attempting additions.
+    """
+    report = dict(added_vertices=0, added_triangles=0, contour_edges=0,
+                  missing_qef=0, escaped_qef=0, unsupported_topology=0,
+                  discarded_components=0, removed_triangles=0, warnings=[],
+                  added_cap_triangles=0, added_transition_triangles=0,
+                  cap_max_plane_error=0.0, closure_success=False)
+    mesh.capping_report = report
+
+    def audit():
+        result = {}
+        v = np.asarray(mesh.vertices)
+        triangles = np.asarray(mesh.edges)
+        valid = np.all(np.isfinite(triangles) & (triangles >= 0) & (triangles < len(v))
+                       & (triangles == np.floor(triangles)), axis=1)
+        result["invalid_triangles"] = int(np.count_nonzero(~valid))
+        result["nonfinite_vertices"] = int(np.count_nonzero(~np.isfinite(v).all(axis=1)))
+        triangles = triangles[valid].astype(np.int64)
+        p = v[triangles]
+        with np.errstate(over="ignore", invalid="ignore"):
+            cross = np.cross(p[:, 1] - p[:, 0], p[:, 2] - p[:, 0])
+        result["nonfinite_triangles"] = int(np.count_nonzero(~np.isfinite(cross).all(axis=1)))
+        result["degenerate_triangles"] = int(np.count_nonzero(np.all(cross == 0, axis=1)))
+        membership = np.zeros((len(v), 3, 2), dtype=bool)
+        membership[:len(cells), :, 0] = cells == 0
+        membership[:len(cells), :, 1] = cells == shape - 1
+        tolerance = 64 * np.finfo(float).eps * np.maximum(1, np.abs(bounds))
+        membership[len(cells):] = np.abs(v[len(cells):, :, None] - bounds) <= tolerance
+        incidence = defaultdict(list)
+        links = defaultdict(list)
+        for a, b, c in triangles:
+            for x, y in ((a, b), (b, c), (c, a)):
+                incidence[tuple(sorted((int(x), int(y))))].append((int(x), int(y)))
+            for x, y, z in ((a, b, c), (b, c, a), (c, a, b)):
+                links[int(x)].append((int(y), int(z)))
+        openings = [edge for edge, entries in incidence.items() if len(entries) == 1]
+        result["open_edges"] = len(openings)
+        result["non_extent_open_edges"] = sum(not np.any(membership[a] & membership[b]) for a, b in openings)
+        result["nonmanifold_edges"] = sum(len(e) > 2 for e in incidence.values())
+        result["orientation_conflicts"] = sum(len(e) == 2 and e[0] == e[1] for e in incidence.values())
+        bad_links = 0
+        for edges in links.values():
+            graph = defaultdict(list)
+            for a, b in edges:
+                graph[a].append(b)
+                graph[b].append(a)
+            seen, todo = set(), [next(iter(graph))]
+            while todo:
+                node = todo.pop()
+                if node not in seen:
+                    seen.add(node)
+                    todo.extend(graph[node])
+            degrees = [len(neighbors) for neighbors in graph.values()]
+            if len(seen) != len(graph) or max(degrees) > 2 or degrees.count(1) not in (0, 2):
+                bad_links += 1
+        result["nonmanifold_vertices"] = bad_links
+        result["watertight"] = bool(incidence) and not any(result[k] for k in (
+            "invalid_triangles", "open_edges", "nonmanifold_edges", "orientation_conflicts", "nonmanifold_vertices"))
+        return result
+
+    def finish():
+        after = audit()
+        report.update(after)
+        report["open_edges_after"] = after["open_edges"]
+        report["non_extent_open_edges_after"] = after["non_extent_open_edges"]
+        report["boundary_edges"] = after["open_edges"]
+        report["closure_success"] = report["watertight"] and not any(report[k] for k in (
+            "missing_qef", "escaped_qef", "unsupported_topology", "invalid_triangles",
+            "nonfinite_vertices", "nonfinite_triangles", "degenerate_triangles"))
+        plane_tolerance = 64 * np.finfo(float).eps * max(1, np.max(np.abs(bounds)))
+        if report["cap_max_plane_error"] > plane_tolerance:
+            report["closure_success"] = False
+            report["warnings"].append("Cap triangles are off the supplied extent planes")
+        for key in ("missing_qef", "escaped_qef", "unsupported_topology", "boundary_edges",
+                    "nonmanifold_edges", "orientation_conflicts", "nonmanifold_vertices",
+                    "invalid_triangles", "nonfinite_vertices", "nonfinite_triangles", "degenerate_triangles"):
+            if report[key]:
+                report["warnings"].append(f"{key}: {report[key]}")
+        if report["warnings"]:
+            warnings.warn("Extent capping best effort: " + "; ".join(report["warnings"]),
+                          RuntimeWarning, stacklevel=2)
+        return mesh
+
+    if not isinstance(mesh.vertices, np.ndarray) or not isinstance(mesh.edges, np.ndarray):
+        report["unsupported_topology"] += 1
+        report["warnings"].append("Only NumPy meshes are supported")
+        warnings.warn(report["warnings"][0], RuntimeWarning, stacklevel=2)
+        return mesh
+    vertices = np.asarray(mesh.vertices)
+    cells = np.asarray(cell_coordinates)
+    if vertices.ndim != 2 or vertices.shape[1] != 3 or cells.shape != vertices.shape:
+        raise ValueError("cell_coordinates must have one (x, y, z) row per QEF vertex")
+    if mesh.edges.ndim != 2 or mesh.edges.shape[1] != 3:
+        raise ValueError("mesh.edges must be triangle indices (M, 3)")
+    bounds = np.asarray(extent, dtype=float).reshape(3, 2)
+    shape = np.asarray(shape)
+    if shape.shape != (3,) or np.any(shape < 1) or not np.isfinite(shape).all() or np.any(shape != np.floor(shape)):
+        raise ValueError("shape must contain three positive integer cell counts")
+    if not np.isfinite(bounds).all() or np.any(bounds[:, 1] <= bounds[:, 0]):
+        raise ValueError("extent must contain finite increasing bounds")
+    shape = shape.astype(int)
+    before = audit()
+    report.update({key + "_before": value for key, value in before.items()})
+    if skip_reason is not None:
+        report["unsupported_topology"] += 1
+        report["warnings"].append(skip_reason)
+        report["skipped_reason"] = skip_reason
+        return finish()
+    if any(before[k] for k in ("invalid_triangles", "nonfinite_vertices", "nonfinite_triangles", "degenerate_triangles")):
+        return finish()
+    coords, positions, faces = (boundary_lattice(shape, extent) if boundary_geometry is None
+                                else boundary_geometry)
+    supplied = np.asarray(boundary_coordinates)
+    scalar = np.asarray(boundary_scalar, dtype=float).reshape(-1)
+    if supplied.shape != (len(scalar), 3):
+        raise ValueError("boundary_coordinates and boundary_scalar must have matching rows")
+    samples = {tuple(c): s for c, s in zip(supplied, scalar)}
+    if len(samples) != len(supplied) or len(samples) != len(coords) or any(tuple(c) not in samples for c in coords):
+        report["unsupported_topology"] += 1
+        report["warnings"].append("Incomplete or duplicate boundary lattice")
+        return finish()
+    values = np.asarray([samples[tuple(c)] for c in coords])
+    # Canonical endpoint order gives every shared edge the identical parameter.
+    source = faces["triangles"]
+    ends = np.roll(source, -1, axis=1)
+    starts, ends = np.minimum(source, ends), np.maximum(source, ends)
+    try:
+        with np.errstate(over="ignore", invalid="ignore"):
+            crossing, parameters = scalar_crossing_parameters(values[starts], values[ends], isovalue)
+    except ValueError as error:
+        report["unsupported_topology"] += 1
+        report["warnings"].append(f"Nonfinite boundary crossing: {error}")
+        return finish()
+    step = (bounds[:, 1] - bounds[:, 0]) / shape
+    qef = {}
+    seen_cells = set()
+    for i, cell in enumerate(cells):
+        key = tuple(cell)
+        if np.any(cell != np.floor(cell)) or np.any(cell < 0) or np.any(cell >= shape) or key in seen_cells:
+            report["unsupported_topology"] += 1
+            continue
+        seen_cells.add(key)
+        low = bounds[:, 0] + cell * step
+        tolerance = 64 * np.finfo(float).eps * np.maximum(1, np.maximum(abs(low), abs(low + step)))
+        if not np.isfinite(vertices[i]).all() or np.any(vertices[i] < low - tolerance) or np.any(vertices[i] > low + step + tolerance):
+            report["escaped_qef"] += 1
+            continue
+        qef[key] = i
+
+    points = list(vertices.astype(float))
+    point_ids = {}
+
+    def point(a, b=None, t=None):
+        # Endpoint ties must share the lattice-point key, not an edge key.
+        if b is None or values[a] == isovalue:
+            key, p = ("p", a), positions[a]
+        elif values[b] == isovalue:
+            key, p = ("p", b), positions[b]
+        else:
+            a, b = sorted((a, b))
+            if t == 0 or t == 1:
+                endpoint = a if t == 0 else b
+                key, p = ("p", endpoint), positions[endpoint]
+            else:
+                key = ("e", a, b)  # Shared by primal edges and face diagonals.
+                p = positions[a] + t * (positions[b] - positions[a])
+        if key not in point_ids:
+            point_ids[key] = len(points)
+            points.append(p)
+        return point_ids[key]
+
+    cap, owners = [], []
+    for i, (triangle, cell) in enumerate(zip(faces["triangles"], faces["cells"])):
+        polygon = []
+        for j, (a, b) in enumerate(zip(triangle, np.roll(triangle, -1))):
+            if values[a] <= isovalue:
+                polygon.append(point(a))
+            if crossing[i, j]:
+                polygon.append(point(a, b, parameters[i, j]))
+        polygon = list(dict.fromkeys(polygon))
+        for j in range(1, len(polygon) - 1):
+            cap.append((polygon[0], polygon[j], polygon[j + 1]))
+            owners.append(tuple(cell))
+
+    # Connected cap patches are seeded only by contour edges with a usable QEF.
+    incidence = defaultdict(list)
+    for i, triangle in enumerate(cap):
+        for a, b in zip(triangle, np.roll(triangle, -1)):
+            incidence[tuple(sorted((a, b)))].append((i, a, b))
+    adjacency = defaultdict(set)
+    contours, seeds = [], set()
+    for entries in incidence.values():
+        if len(entries) == 2:
+            i, j = entries[0][0], entries[1][0]
+            adjacency[i].add(j)
+            adjacency[j].add(i)
+        elif len(entries) == 1:
+            i, a, b = entries[0]
+            report["contour_edges"] += 1
+            q = qef.get(owners[i])
+            if q is None:
+                report["missing_qef"] += 1
+            else:
+                seeds.add(i)
+                contours.append((i, a, b, q))
+        else:
+            report["unsupported_topology"] += 1
+    retained, visited = set(), set()
+    for start in range(len(cap)):
+        if start in visited:
+            continue
+        component, todo = set(), [start]
+        while todo:
+            i = todo.pop()
+            if i not in component:
+                component.add(i)
+                todo.extend(adjacency[i] - component)
+        visited.update(component)
+        if component & seeds:
+            retained.update(component)
+        else:
+            report["discarded_components"] += 1
+    added = [cap[i] for i in sorted(retained)]
+    cap_count = len(added)
+    fans = [(b, a, q) for i, a, b, q in contours if i in retained]
+    added.extend(fans)
+
+    # The two in-domain cells around a boundary primal edge need one more
+    # triangle. Derive its winding from the existing directed face-fan edges.
+    spokes = defaultdict(list)
+    for b, a, q in fans:
+        spokes[(a, q)].append((a, q))
+        spokes[(b, q)].append((q, b))
+    residual = defaultdict(list)
+    for (p, q), edges in spokes.items():
+        forward = sum(a == p for a, b in edges)
+        reverse = len(edges) - forward
+        if forward != reverse:
+            residual[p].extend([(p, q) if forward > reverse else (q, p)] * abs(forward - reverse))
+    for p, edges in residual.items():
+        if len(edges) == 2:
+            incoming = [a for a, b in edges if b == p]
+            outgoing = [b for a, b in edges if a == p]
+            if (len(incoming) == len(outgoing) == 1
+                    and np.sum(np.abs(cells[incoming[0]] - cells[outgoing[0]])) == 1):
+                added.append((p, incoming[0], outgoing[0]))
+                continue
+        report["unsupported_topology"] += 1
+
+    # Never remove/reindex original geometry, even if it already has defects.
+    seen = {tuple(sorted(t)) for t in mesh.edges}
+    clean = []
+    for index, triangle in enumerate(added):
+        key = tuple(sorted(triangle))
+        a, b, c = (points[i] for i in triangle)
+        cross = np.cross(b - a, c - a)
+        if len(set(triangle)) < 3 or key in seen or not np.any(cross) or not np.isfinite(cross).all():
+            report["removed_triangles"] += 1
+            continue
+        seen.add(key)
+        clean.append(triangle)
+        if index < cap_count:
+            report["added_cap_triangles"] += 1
+            # A cap triangle must have all three vertices on the same plane.
+            distance = np.abs(np.asarray([a, b, c])[:, :, None] - bounds)
+            report["cap_max_plane_error"] = max(report["cap_max_plane_error"], float(distance.max(axis=0).min()))
+        else:
+            report["added_transition_triangles"] += 1
+    used = sorted({i for t in clean for i in t if i >= len(vertices)})
+    remap = {old: len(vertices) + i for i, old in enumerate(used)}
+    if clean:
+        mesh.vertices = np.concatenate((vertices, np.asarray([points[i] for i in used]).reshape(-1, 3)))
+        new_edges = np.asarray([[remap.get(i, i) for i in t] for t in clean], dtype=np.int64)
+        mesh.edges = np.concatenate((mesh.edges, new_edges))
+    report["added_vertices"], report["added_triangles"] = len(used), len(clean)
+    return finish()
diff --git a/gempy_engine/modules/dual_contouring/_gen_vertices.py b/gempy_engine/modules/dual_contouring/_gen_vertices.py
index a622e25b..dca24d13 100644
--- a/gempy_engine/modules/dual_contouring/_gen_vertices.py
+++ b/gempy_engine/modules/dual_contouring/_gen_vertices.py
@@ -32,7 +32,13 @@ def generate_dual_contouring_vertices(dc_data_per_stack: DualContouringData, sli
     # Use nanmean directly without intermediate copy
     bias_xyz_slice = edges_xyz[:, :12]
     
-    if BackendTensor.engine_backend == AvailableBackends.PYTORCH:
+    if dc_data_per_stack.strict_crossings:
+        # Zero coordinates are valid samples, not missing edge constraints.
+        mask = valid_edges_bool[:, :, None]
+        sum_valid = (bias_xyz_slice * mask).sum(axis=1)
+        count_valid = mask.sum(axis=1)
+        mass_points = sum_valid / count_valid
+    elif BackendTensor.engine_backend == AvailableBackends.PYTORCH:
         mask = bias_xyz_slice == 0
         bias_xyz_masked = BackendTensor.tfnp.where(mask, float('nan'), bias_xyz_slice)
         mass_points = BackendTensor.tfnp.nanmean(bias_xyz_masked, axis=1)
@@ -96,7 +102,10 @@ def generate_dual_contouring_vertices(dc_data_per_stack: DualContouringData, sli
         
         # Solve ATA @ x = ATb  (use solve instead of inv for numerical stability)
         import torch
-        reg = 1e-4 * torch.eye(3, device=ATA.device, dtype=ATA.dtype).unsqueeze(0)
+        # The mass-point constraints already make ATA positive definite. Extra
+        # origin-centered regularization biases the capped solid's volume.
+        strength = 0 if dc_data_per_stack.strict_crossings else 1e-4
+        reg = strength * torch.eye(3, device=ATA.device, dtype=ATA.dtype).unsqueeze(0)
         vertices = torch.linalg.solve(ATA + reg, ATb)
     else:
         # NumPy: use efficient einsum
diff --git a/gempy_engine/modules/dual_contouring/_scalar_crossing.py b/gempy_engine/modules/dual_contouring/_scalar_crossing.py
new file mode 100644
index 00000000..e17c1a71
--- /dev/null
+++ b/gempy_engine/modules/dual_contouring/_scalar_crossing.py
@@ -0,0 +1,32 @@
+"""Shared scalar-only edge classification for contouring and NumPy cappers."""
+
+import numpy as np
+
+
+def scalar_crossing_parameters(start, end, iso, *, xp=np):
+    """Return broadcast ``(crossing, t)`` for ``(1-t)*start + t*end == iso``.
+
+    Inputs are floating arrays in the supplied NumPy or PyTorch namespace.
+    ``scalar <= iso`` is inside: an entirely iso-valued edge never crosses,
+    and an equal endpoint crosses only if the other endpoint is outside.
+    Parameters are clamped to [0, 1], snapped within 1e-12 of an endpoint,
+    and zero on non-crossing edges.
+    Nonfinite inputs or interpolation arithmetic raise ValueError.
+    The default namespace is NumPy, independent of the engine backend.
+    """
+    for name, value in (("start", start), ("end", end), ("iso", iso)):
+        if not bool(xp.isfinite(value).all()):
+            raise ValueError(f"Scalar crossing requires finite {name} values")
+
+    crossing = (start <= iso) != (end <= iso)
+    # Neutralize unused edges before dividing, without perturbing crossings.
+    denominator = xp.where(crossing, end, 0) - xp.where(crossing, start, 0)
+    numerator = xp.where(crossing, iso, 0) - xp.where(crossing, start, 0)
+    if not bool((xp.isfinite(denominator) & xp.isfinite(numerator)).all()):
+        raise ValueError("Nonfinite scalar crossing interpolation arithmetic")
+    t = numerator / xp.where(crossing, denominator, 1)
+    if not bool(xp.isfinite(t).all()):
+        raise ValueError("Nonfinite scalar crossing interpolation parameter")
+    t = xp.clip(t, 0, 1)
+    t = xp.where(t <= 1e-12, 0, xp.where(t >= 1 - 1e-12, 1, t))
+    return crossing, t
diff --git a/gempy_engine/modules/dual_contouring/_support_report.py b/gempy_engine/modules/dual_contouring/_support_report.py
index 08a32bac..20f8aac7 100644
--- a/gempy_engine/modules/dual_contouring/_support_report.py
+++ b/gempy_engine/modules/dual_contouring/_support_report.py
@@ -7,7 +7,7 @@
 
 
 def mesh_support_report(coordinates, scalar_corners, isovalue, domain_shape,
-                        mask=None, surface_index=0, ancestor_coordinates=()):
+                        mask=None, surface_index=0, ancestor_coordinates=(), strict_crossings=False):
     """Classify missing incident cells for unique sampled sign-changing edges.
 
     This diagnoses sampled crossings, not unsampled components or later triangle
@@ -30,7 +30,10 @@ def mesh_support_report(coordinates, scalar_corners, isovalue, domain_shape,
         ((0, 1), (2, 3), (4, 5), (6, 7)),
     )):
         for a, b in pairs:
-            crossing = (scalar[:, a] >= iso) != (scalar[:, b] >= iso)
+            if strict_crossings:
+                crossing = (scalar[:, a] <= iso) != (scalar[:, b] <= iso)
+            else:
+                crossing = (scalar[:, a] >= iso) != (scalar[:, b] >= iso)
             edges.update((direction, *p) for p in coords[crossing] + corners[a])
     report = dict(surface_index=surface_index, crossing_edge_count=len(edges),
                   missing_incident_cell_count=0, physical_boundary_edge_count=0,
diff --git a/gempy_engine/modules/dual_contouring/dual_contouring_interface.py b/gempy_engine/modules/dual_contouring/dual_contouring_interface.py
index 6bcb9ce1..d1f75aef 100644
--- a/gempy_engine/modules/dual_contouring/dual_contouring_interface.py
+++ b/gempy_engine/modules/dual_contouring/dual_contouring_interface.py
@@ -5,6 +5,7 @@
 
 from ._find_vertex_overlap import find_repeated_voxels_across_stacks
 from ._apply_vertex_overlap_logic import apply_relations_to_overlaps
+from ._scalar_crossing import scalar_crossing_parameters
 from .fancy_triangulation import get_left_right_array
 from ...config import AvailableBackends
 from ...core.backend_tensor import BackendTensor
@@ -20,13 +21,44 @@
 # region edges
 
 def find_intersection_on_edge(_xyz_corners, scalar_field_on_corners,
-                              scalar_at_sp, masking=None) -> Tuple:
+                              scalar_at_sp, masking=None, *, strict_crossings=False) -> Tuple:
+    """Find edge intersections, optionally using exact scalar-side classification."""
+    if strict_crossings:
+        return _find_strict_intersections(_xyz_corners, scalar_field_on_corners, scalar_at_sp, masking)
     if BackendTensor.engine_backend == AvailableBackends.PYTORCH:
         return find_intersection_on_edge_torch(_xyz_corners, scalar_field_on_corners, scalar_at_sp, masking)
     else:
         return find_intersection_on_edge_numpy(_xyz_corners, scalar_field_on_corners, scalar_at_sp, masking)
 
 
+def _find_strict_intersections(xyz_corners, scalars, isovalues, masking):
+    xp = BackendTensor.t
+    xyz = xyz_corners.reshape(-1, 8, 3)
+    scalars = scalars.reshape(1, -1, 8)
+    if masking is not None:
+        xyz = xyz[masking]
+        scalars = scalars[:, masking]
+    if not bool(xp.isfinite(xyz).all()):
+        raise ValueError("Strict edge crossings require finite corner coordinates")
+
+    # Preserve the legacy x/y/z edge ordering and start-to-end direction.
+    start = [4, 5, 6, 7, 2, 3, 6, 7, 1, 3, 5, 7]
+    end = [0, 1, 2, 3, 0, 1, 4, 5, 0, 2, 4, 6]
+    valid, t = scalar_crossing_parameters(
+        scalars[:, :, start], scalars[:, :, end], isovalues.reshape(-1, 1, 1), xp=xp
+    )
+    shape = (*valid.shape, 3)
+    a = xp.broadcast_to(xyz[None, :, start, :], shape)[valid]
+    b = xp.broadcast_to(xyz[None, :, end, :], shape)[valid]
+    weights = t[valid][:, None]
+    intersections = (1 - weights) * a + weights * b
+    if not bool(xp.isfinite(intersections).all()):
+        raise ValueError("Nonfinite strict edge intersection coordinates")
+    if BackendTensor.engine_backend == AvailableBackends.PYTORCH:
+        return intersections, valid.reshape(-1)
+    return intersections, valid.reshape(-1, 12)
+
+
 def find_intersection_on_edge_numpy(_xyz_corners, scalar_field_on_corners,
                                     scalar_at_sp, masking=None) -> Tuple:
     """This function finds all the intersections for multiple layers per series
diff --git a/tests/test_common/test_modules/test_extent_capping.py b/tests/test_common/test_modules/test_extent_capping.py
new file mode 100644
index 00000000..df24c748
--- /dev/null
+++ b/tests/test_common/test_modules/test_extent_capping.py
@@ -0,0 +1,342 @@
+"""Analytic interior DC meshes, independent of the capping triangulator."""
+
+from itertools import product
+from types import SimpleNamespace
+import warnings
+
+import numpy as np
+import pytest
+
+from gempy_engine.modules.dual_contouring._extent_capping import boundary_lattice, cap_mesh
+from gempy_engine.modules.dual_contouring import _extent_capping
+
+
+def plane_mesh(shape, extent, normal, level):
+    bounds = np.asarray(extent).reshape(3, 2)
+    shape = np.asarray(shape)
+    normal = np.asarray(normal)
+    step = (bounds[:, 1] - bounds[:, 0]) / shape
+    positions = lambda c: bounds[:, 0] + np.asarray(c) * step
+    scalar = lambda c: positions(c) @ normal - level
+    vertices, cells, lookup = [], [], {}
+    for cell in product(*(range(n) for n in shape)):
+        intersections = []
+        for axis in range(3):
+            others = [i for i in range(3) if i != axis]
+            for offsets in product((0, 1), repeat=2):
+                a = np.array(cell)
+                a[others] += offsets
+                b = a.copy()
+                b[axis] += 1
+                sa, sb = scalar(a), scalar(b)
+                if (sa <= 0) != (sb <= 0):
+                    intersections.append(positions(a) + sa / (sa - sb) * (positions(b) - positions(a)))
+        if intersections:
+            lookup[cell] = len(vertices)
+            cells.append(cell)
+            vertices.append(np.mean(intersections, axis=0))
+    triangles = []
+    for axis in range(3):
+        u, v = (axis + 1) % 3, (axis + 2) % 3
+        ranges = [range(n + 1) for n in shape]
+        ranges[axis] = range(shape[axis])
+        for corner in product(*ranges):
+            a = np.array(corner)
+            b = a.copy()
+            b[axis] += 1
+            if (scalar(a) <= 0) == (scalar(b) <= 0):
+                continue
+            ring = []
+            for du, dv in ((-1, -1), (0, -1), (0, 0), (-1, 0)):
+                c = a.copy()
+                c[u] += du
+                c[v] += dv
+                if tuple(c) in lookup:
+                    ring.append(lookup[tuple(c)])
+            if len(ring) != 4:
+                continue
+            for triangle in ((ring[0], ring[1], ring[2]), (ring[0], ring[2], ring[3])):
+                p, q, r = [vertices[i] for i in triangle]
+                if np.dot(np.cross(q - p, r - p), normal) < 0:
+                    triangle = triangle[::-1]
+                triangles.append(triangle)
+    return SimpleNamespace(vertices=np.asarray(vertices).reshape(-1, 3),
+                           edges=np.asarray(triangles, dtype=int).reshape(-1, 3)), np.asarray(cells).reshape(-1, 3)
+
+
+def assert_closed(mesh):
+    report = mesh.capping_report
+    assert report["watertight"], report
+    assert report["closure_success"], report
+    assert report["added_triangles"] == report["added_cap_triangles"] + report["added_transition_triangles"]
+    assert report["cap_max_plane_error"] == 0
+    assert not report["warnings"], report
+    assert len({tuple(sorted(t)) for t in mesh.edges}) == len(mesh.edges)
+    triangles = mesh.vertices[mesh.edges]
+    assert np.all(np.linalg.norm(np.cross(triangles[:, 1] - triangles[:, 0],
+                                         triangles[:, 2] - triangles[:, 0]), axis=1) > 0)
+
+
+def test_boundary_lattice_shared_points_and_outward_faces():
+    shape, extent = (2, 3, 4), (-7, 2, 10, 16, -2, 6)
+    coordinates, positions, metadata = boundary_lattice(shape, extent)
+    assert coordinates.dtype.kind == "i"
+    assert len(coordinates) == 3 * 4 * 5 - 1 * 2 * 3
+    assert len(np.unique(coordinates, axis=0)) == len(coordinates)
+    assert metadata["triangles"].shape == (4 * (2 * 3 + 3 * 4 + 4 * 2), 3)
+    center = np.asarray(extent).reshape(3, 2).mean(axis=1)
+    for triangle in positions[metadata["triangles"]]:
+        a, b, c = triangle
+        assert np.dot(np.cross(b - a, c - a), triangle.mean(axis=0) - center) > 0
+    np.testing.assert_array_equal(positions.min(axis=0), np.asarray(extent)[::2])
+    np.testing.assert_array_equal(positions.max(axis=0), np.asarray(extent)[1::2])
+
+
+@pytest.mark.parametrize("normal,level", [((1, .37, -.21), .731), ((1, 0, 0), .43),
+                                        ((1, 1, 1), .71), ((-.3, 1, .7), .52)])
+def test_plane_closed_and_original_indices_preserved(normal, level):
+    shape, extent = (4, 5, 3), (0, 1, 0, 1, 0, 1)
+    mesh, cells = plane_mesh(shape, extent, normal, level)
+    original_vertices, original_triangles = mesh.vertices.copy(), mesh.edges.copy()
+    # Neither the QEF rows nor boundary sample order is assumed sorted.
+    permutation = np.arange(len(cells))[::-1]
+    inverse = np.argsort(permutation)
+    mesh.vertices = mesh.vertices[permutation]
+    mesh.edges = inverse[mesh.edges]
+    coordinates, positions, _ = boundary_lattice(shape, extent)
+    assert cap_mesh(mesh, cells[permutation], shape, extent, coordinates[::-1],
+                    (positions @ normal)[::-1], level) is mesh
+    assert_closed(mesh)
+    assert mesh.capping_report["open_edges_before"] > 0
+    assert mesh.capping_report["non_extent_open_edges_before"] == 0
+    assert mesh.capping_report["open_edges_after"] == 0
+    assert mesh.capping_report["non_extent_open_edges_after"] == 0
+    assert mesh.capping_report["added_cap_triangles"] > 0
+    assert mesh.capping_report["added_transition_triangles"] > 0
+    np.testing.assert_array_equal(mesh.vertices[:len(cells)], original_vertices[permutation])
+    np.testing.assert_array_equal(mesh.edges[:len(original_triangles)], inverse[original_triangles])
+    new_triangles = mesh.edges[len(original_triangles):]
+    assert np.count_nonzero(np.all(new_triangles >= len(cells), axis=1)) == mesh.capping_report["added_cap_triangles"]
+    assert np.count_nonzero(np.any(new_triangles < len(cells), axis=1)) == mesh.capping_report["added_transition_triangles"]
+    assert np.all(mesh.vertices >= 0) and np.all(mesh.vertices <= 1)
+    # Shared box edges have one vertex per position, not one per face.
+    new = mesh.vertices[len(cells):]
+    assert len(np.unique(new, axis=0)) == len(new)
+    t = mesh.vertices[mesh.edges]
+    assert np.sum(np.einsum("ij,ij->i", t[:, 0], np.cross(t[:, 1], t[:, 2]))) > 0
+
+
+def test_exact_extent_nonunit_plane():
+    shape, extent = (3, 4, 5), (-2, 4, 10, 14, -8, -3)
+    normal, level = (1, .31, .17), 3.19
+    mesh, cells = plane_mesh(shape, extent, normal, level)
+    coordinates, positions, _ = boundary_lattice(shape, extent)
+    cap_mesh(mesh, cells, shape, extent, coordinates, positions @ normal, level)
+    assert_closed(mesh)
+    np.testing.assert_array_equal(mesh.vertices.min(axis=0), [-2, 10, -8])
+    np.testing.assert_array_equal(mesh.vertices.max(axis=0)[1:], [14, -3])
+
+
+@pytest.mark.parametrize("level", [1., 1.5])
+def test_corner_ties_are_shared_and_limitations_reported(level):
+    shape, extent = (4, 4, 4), (0, 1, 0, 1, 0, 1)
+    mesh, cells = plane_mesh(shape, extent, (1, 1, 1), level)
+    coordinates, positions, _ = boundary_lattice(shape, extent)
+    with warnings.catch_warnings(record=True) as emitted:
+        warnings.simplefilter("always")
+        cap_mesh(mesh, cells, shape, extent, coordinates, positions.sum(axis=1), level)
+    new = mesh.vertices[len(cells):]
+    assert len(np.unique(new, axis=0)) == len(new)
+    assert mesh.capping_report["added_triangles"] > 0
+    assert mesh.capping_report["watertight"] or (emitted and mesh.capping_report["warnings"])
+
+
+@pytest.mark.parametrize("inside", [False, True])
+def test_no_boundary_contour_adds_no_shell(inside):
+    # An enclosed sphere represented by an outward-oriented octahedron.
+    vertices = np.vstack((np.eye(3), -np.eye(3))) * .2 + .5
+    triangles = []
+    for x, y, z in product((0, 3), (1, 4), (2, 5)):
+        a, b, c = vertices[[x, y, z]]
+        t = [x, y, z]
+        if np.dot(np.cross(b - a, c - a), a - .5) < 0:
+            t.reverse()
+        triangles.append(t)
+    mesh = SimpleNamespace(vertices=vertices.copy(), edges=np.asarray(triangles))
+    coordinates, positions, _ = boundary_lattice((10, 10, 10), (0, 1, 0, 1, 0, 1))
+    values = np.linalg.norm(positions - .5, axis=1) - .2
+    cap_mesh(mesh, np.floor(vertices * 10).astype(int), (10, 10, 10), (0, 1, 0, 1, 0, 1),
+             coordinates, -values if inside else values, 0)
+    np.testing.assert_array_equal(mesh.vertices, vertices)
+    np.testing.assert_array_equal(mesh.edges, triangles)
+    assert mesh.capping_report["added_triangles"] == 0
+    assert mesh.capping_report["contour_edges"] == 0
+    assert_closed(mesh)
+
+
+def test_missing_escaped_qef_and_internal_hole_are_not_repaired():
+    shape, extent = (5, 5, 5), (0, 1, 0, 1, 0, 1)
+    mesh, cells = plane_mesh(shape, extent, (1, .37, -.21), .73)
+    mesh.vertices[0] = [-2, -2, -2]
+    # Remove an interior triangle, not a boundary feature.
+    mesh.edges = mesh.edges[1:].copy()
+    original = mesh.edges.copy()
+    coordinates, positions, _ = boundary_lattice(shape, extent)
+    with pytest.warns(RuntimeWarning, match="best effort"):
+        cap_mesh(mesh, cells, shape, extent, coordinates, positions @ (1, .37, -.21), .73)
+    report = mesh.capping_report
+    assert report["escaped_qef"] == 1
+    assert report["missing_qef"] > 0
+    assert report["boundary_edges"] > 0
+    assert not report["watertight"]
+    np.testing.assert_array_equal(mesh.edges[:len(original)], original)
+    assert all(np.any(t >= len(cells)) for t in mesh.edges[len(original):])
+
+
+def test_no_available_qef_discards_all_cap_components():
+    mesh = SimpleNamespace(vertices=np.empty((0, 3)), edges=np.empty((0, 3), dtype=int))
+    coordinates, positions, _ = boundary_lattice((2, 2, 2), (0, 1, 0, 1, 0, 1))
+    with pytest.warns(RuntimeWarning, match="missing_qef"):
+        cap_mesh(mesh, np.empty((0, 3)), (2, 2, 2), (0, 1, 0, 1, 0, 1),
+                 coordinates, positions[:, 0], .3)
+    assert mesh.capping_report["discarded_components"] > 0
+    assert mesh.vertices.shape == (0, 3)
+    assert mesh.edges.shape == (0, 3)
+
+
+def test_nonfinite_boundary_is_diagnostic_not_geometry():
+    mesh = SimpleNamespace(vertices=np.empty((0, 3)), edges=np.empty((0, 3), dtype=int))
+    coordinates, positions, _ = boundary_lattice((1, 1, 1), (0, 1, 0, 1, 0, 1))
+    values = positions[:, 0].copy()
+    values[0] = np.nan
+    with pytest.warns(RuntimeWarning, match="Nonfinite"):
+        cap_mesh(mesh, np.empty((0, 3)), (1, 1, 1), (0, 1, 0, 1, 0, 1), coordinates, values, .3)
+    assert mesh.capping_report["unsupported_topology"] == 1
+    assert mesh.capping_report["added_triangles"] == 0
+
+
+def test_only_component_connected_to_available_qef_is_retained():
+    shape, extent = (2, 2, 2), (0, 1, 0, 1, 0, 1)
+    coordinates, positions, _ = boundary_lattice(shape, extent)
+    values = np.ones(len(coordinates))
+    values[np.all(coordinates == 0, axis=1) | np.all(coordinates == 2, axis=1)] = -1
+    mesh = SimpleNamespace(vertices=np.asarray([[.1, .1, .1]]),
+                           edges=np.empty((0, 3), dtype=int))
+    with pytest.warns(RuntimeWarning, match="missing_qef"):
+        cap_mesh(mesh, np.asarray([[0, 0, 0]]), shape, extent, coordinates, values, 0)
+    assert mesh.capping_report["discarded_components"] == 1
+    assert mesh.capping_report["watertight"]
+    assert not mesh.capping_report["closure_success"]
+    assert np.all(mesh.vertices <= .25)
+
+
+def test_internal_hole_without_missing_qef_is_preserved():
+    shape, extent = (6, 6, 6), (0, 1, 0, 1, 0, 1)
+    mesh, cells = plane_mesh(shape, extent, (1, .37, -.21), .731)
+    interior = np.flatnonzero(np.all((cells[mesh.edges] > 0) & (cells[mesh.edges] < 5), axis=(1, 2)))
+    assert len(interior)
+    removed = mesh.edges[interior[0]].copy()
+    mesh.edges = np.delete(mesh.edges, interior[0], axis=0)
+    coordinates, positions, _ = boundary_lattice(shape, extent)
+    with pytest.warns(RuntimeWarning, match="boundary_edges: 3"):
+        cap_mesh(mesh, cells, shape, extent, coordinates, positions @ (1, .37, -.21), .731)
+    assert mesh.capping_report["boundary_edges"] == 3
+    assert mesh.capping_report["missing_qef"] == 0
+    assert mesh.capping_report["escaped_qef"] == 0
+    assert mesh.capping_report["non_extent_open_edges_before"] == 3
+    assert mesh.capping_report["non_extent_open_edges_after"] == 3
+    assert not mesh.capping_report["closure_success"]
+    assert tuple(sorted(removed)) not in {tuple(sorted(t)) for t in mesh.edges}
+
+
+def test_prebuilt_geometry_and_single_batched_strict_crossing(monkeypatch):
+    shape, extent = (3, 4, 5), (0, 1, 0, 1, 0, 1)
+    geometry = boundary_lattice(shape, extent)
+    coordinates, positions, metadata = geometry
+    original_coordinates, original_positions = coordinates.copy(), positions.copy()
+    original_triangles, original_cells = metadata["triangles"].copy(), metadata["cells"].copy()
+    mesh, cells = plane_mesh(shape, extent, (1, 0, 0), .43)
+    # Every QEF on this plane is inset from the actual extent planes, including
+    # the rim. Its owning cell, not its position, identifies an extent opening.
+    assert np.all((mesh.vertices > 0) & (mesh.vertices < 1))
+    original_crossing = _extent_capping.scalar_crossing_parameters
+    calls = []
+
+    def crossing(start, end, iso, *, xp=np):
+        calls.append(start.shape)
+        return original_crossing(start, end, iso, xp=xp)
+
+    def unexpected_lattice(*args):
+        pytest.fail("prebuilt boundary geometry must not be regenerated")
+
+    monkeypatch.setattr(_extent_capping, "boundary_lattice", unexpected_lattice)
+    monkeypatch.setattr(_extent_capping, "scalar_crossing_parameters", crossing)
+    # Tiny but nonzero scalar differences must not get tolerance-snapped.
+    values = (positions[:, 0] - .43) * 1e-200
+    cap_mesh(mesh, cells, shape, extent, coordinates, values, 0, boundary_geometry=geometry)
+    assert calls == [metadata["triangles"].shape]
+    assert_closed(mesh)
+    assert mesh.capping_report["open_edges_before"] > 0
+    assert mesh.capping_report["non_extent_open_edges_before"] == 0
+    new = mesh.vertices[len(cells):]
+    assert np.any(np.isclose(new[:, 0], .43, atol=1e-15))
+    np.testing.assert_array_equal(coordinates, original_coordinates)
+    np.testing.assert_array_equal(positions, original_positions)
+    np.testing.assert_array_equal(metadata["triangles"], original_triangles)
+    np.testing.assert_array_equal(metadata["cells"], original_cells)
+
+
+def test_crossing_arithmetic_overflow_is_reported():
+    shape, extent = (1, 1, 1), (0, 1, 0, 1, 0, 1)
+    coordinates, positions, _ = boundary_lattice(shape, extent)
+    values = np.where(positions[:, 0] == 0, -1e308, 1e308)
+    mesh = SimpleNamespace(vertices=np.empty((0, 3)), edges=np.empty((0, 3), dtype=int))
+    with pytest.warns(RuntimeWarning, match="Nonfinite scalar crossing interpolation arithmetic"):
+        cap_mesh(mesh, np.empty((0, 3)), shape, extent, coordinates, values, 0)
+    assert mesh.capping_report["unsupported_topology"] == 1
+    assert not mesh.capping_report["closure_success"]
+    assert mesh.capping_report["added_triangles"] == 0
+
+
+@pytest.mark.parametrize("defect,field", [
+    ("nonfinite", "nonfinite_vertices"),
+    ("negative_index", "invalid_triangles"),
+    ("out_of_range", "invalid_triangles"),
+    ("fractional_index", "invalid_triangles"),
+    ("nonfinite_index", "invalid_triangles"),
+    ("degenerate", "degenerate_triangles"),
+    ("escaped", "escaped_qef"),
+    ("unsupported", "unsupported_topology"),
+])
+def test_original_geometry_defects_never_report_success(defect, field):
+    shape, extent = (3, 3, 3), (0, 1, 0, 1, 0, 1)
+    vertices = np.asarray([[.2, .2, .2], [.7, .2, .2], [.2, .7, .2], [.2, .2, .7]])
+    cells = np.floor(vertices * 3).astype(int)
+    triangles = np.asarray([[0, 2, 1], [0, 1, 3], [0, 3, 2], [1, 2, 3]])
+    if defect == "nonfinite":
+        vertices[0, 0] = np.nan
+    elif defect == "negative_index":
+        triangles[0, 0] = -1
+    elif defect == "out_of_range":
+        triangles[0, 0] = len(vertices)
+    elif defect in ("fractional_index", "nonfinite_index"):
+        triangles = triangles.astype(float)
+        triangles[0, 0] = .5 if defect == "fractional_index" else np.nan
+    elif defect == "degenerate":
+        vertices[0] = vertices[1]
+    elif defect == "escaped":
+        vertices[0, 0] = -.1
+    elif defect == "unsupported":
+        cells[0] = cells[1]
+    mesh = SimpleNamespace(vertices=vertices.copy(), edges=triangles.copy())
+    coordinates, positions, _ = boundary_lattice(shape, extent)
+    with pytest.warns(RuntimeWarning, match=field):
+        cap_mesh(mesh, cells, shape, extent, coordinates, np.ones(len(positions)), 0)
+    report = mesh.capping_report
+    assert report[field] > 0
+    assert not report["closure_success"]
+    assert report["added_triangles"] == 0
+    np.testing.assert_array_equal(mesh.vertices, vertices)
+    np.testing.assert_array_equal(mesh.edges, triangles)
+    if defect in ("nonfinite", "degenerate", "escaped", "unsupported"):
+        assert report["watertight"]  # Topological closure alone is insufficient.
diff --git a/tests/test_common/test_modules/test_extent_capping_integration.py b/tests/test_common/test_modules/test_extent_capping_integration.py
new file mode 100644
index 00000000..4dcc4c21
--- /dev/null
+++ b/tests/test_common/test_modules/test_extent_capping_integration.py
@@ -0,0 +1,412 @@
+"""Extent capping through the production model and boundary evaluation APIs."""
+
+import importlib
+import json
+import warnings
+from copy import deepcopy
+from types import SimpleNamespace
+from unittest.mock import Mock
+
+import numpy as np
+import pytest
+
+from gempy_engine.API.model.model_api import compute_model
+from gempy_engine.config import AvailableBackends
+from gempy_engine.core.backend_tensor import BackendTensor
+from gempy_engine.core.data import TensorsStructure
+from gempy_engine.core.data.engine_grid import EngineGrid
+from gempy_engine.core.data.input_data_descriptor import InputDataDescriptor
+from gempy_engine.core.data.interpolation_functions import CustomInterpolationFunctions
+from gempy_engine.core.data.interpolation_input import InterpolationInput
+from gempy_engine.core.data.kernel_classes.orientations import Orientations
+from gempy_engine.core.data.kernel_classes.surface_points import SurfacePoints
+from gempy_engine.core.data.options import InterpolationOptions
+from gempy_engine.core.data.options.evaluation_options import MeshExtentCapping
+from gempy_engine.core.data.regular_grid import RegularGrid
+from gempy_engine.core.data.stack_relation_type import StackRelationType
+from gempy_engine.core.data.stacks_structure import StacksStructure
+
+
+dc = importlib.import_module("gempy_engine.API.dual_contouring.multi_scalar_dual_contouring")
+EXTENT = np.array([-2., 4., 10., 14., -8., -3.])
+
+
+@pytest.fixture(params=[AvailableBackends.numpy, AvailableBackends.PYTORCH], ids=["numpy", "pytorch-cpu"])
+def backend(request, monkeypatch):
+    torch = pytest.importorskip("torch") if request.param is AvailableBackends.PYTORCH else None
+    grad_enabled = torch.is_grad_enabled() if torch is not None else None
+    saved = dict(engine_backend=BackendTensor.engine_backend, use_gpu=BackendTensor.use_gpu,
+                 use_pykeops=BackendTensor.use_pykeops, dtype=BackendTensor.dtype,
+                 grads=BackendTensor.COMPUTE_GRADS)
+    saved_pykeops = BackendTensor.pykeops_enabled
+    monkeypatch.setenv("GEMPY_SKIP_TRIANGULATION", "0")
+    try:
+        BackendTensor._change_backend(request.param, use_gpu=False, use_pykeops=False, dtype="float64")
+        BackendTensor.pykeops_enabled = False
+        yield request.param
+    finally:
+        BackendTensor._change_backend(**saved)
+        BackendTensor.COMPUTE_GRADS = saved["grads"]
+        BackendTensor.pykeops_enabled = saved_pykeops
+        if torch is not None:
+            torch.set_grad_enabled(grad_enabled)
+
+
+@pytest.fixture
+def plane_model(backend):
+    def make(levels=(-5.7,), normal=(0., 0., 1.), extent=EXTENT, resolution=(3, 3, 3)):
+        functions = CustomInterpolationFunctions(
+            scalar_field_at_surface_points=np.asarray(levels),
+            implicit_function=lambda xyz: xyz[:, 0] * normal[0] + xyz[:, 1] * normal[1] + xyz[:, 2] * normal[2],
+            gx_function=lambda xyz: xyz[:, 0] * 0 + normal[0],
+            gy_function=lambda xyz: xyz[:, 0] * 0 + normal[1],
+            gz_function=lambda xyz: xyz[:, 0] * 0 + normal[2],
+        )
+        structure = InputDataDescriptor(
+            TensorsStructure(number_of_points_per_surface=np.array([], dtype=int)),
+            StacksStructure(
+                number_of_points_per_stack=np.array([0]),
+                number_of_orientations_per_stack=np.array([0]),
+                number_of_surfaces_per_stack=np.array([len(levels)]),
+                masking_descriptor=[StackRelationType.BASEMENT],
+                interp_functions_per_stack=[functions],
+            ),
+        )
+        inputs = InterpolationInput(
+            SurfacePoints(np.empty((0, 3))),
+            Orientations(np.empty((0, 3)), np.empty((0, 3))),
+            EngineGrid(octree_grid=RegularGrid(np.array(extent, dtype=float), list(resolution))),
+            np.arange(len(levels) + 1),
+        )
+        options = InterpolationOptions.from_args(10., 1.)
+        options.evaluation_options.number_octree_levels = 2
+        options.evaluation_options.number_octree_levels_surface = 2
+        return inputs, options, structure
+
+    return make
+
+
+@pytest.mark.parametrize("mode", list(MeshExtentCapping))
+def test_capping_option_json_roundtrip(mode):
+    options = InterpolationOptions.from_args(10., 1.)
+    assert options.evaluation_options.mesh_extraction_extent_capping is MeshExtentCapping.NONE
+    options.evaluation_options.mesh_extraction_extent_capping = mode
+    serialized = options.model_dump_json()
+    assert json.loads(serialized)["evaluation_options"]["mesh_extraction_extent_capping"] == mode.value
+    restored = InterpolationOptions.model_validate_json(serialized)
+    assert restored.evaluation_options.mesh_extraction_extent_capping is mode
+
+
+def test_default_and_explicit_none_have_identical_arrays(plane_model, monkeypatch):
+    boundary = Mock(side_effect=AssertionError("Disabled capping must not evaluate the boundary"))
+    monkeypatch.setattr(dc, "_interp_on_boundary", boundary)
+    inputs, options, descriptor = plane_model()
+    default = compute_model(inputs, options, descriptor)
+    inputs, options, descriptor = plane_model()
+    options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.NONE
+    explicit = compute_model(inputs, options, descriptor)
+    assert len(default.dc_meshes) == len(explicit.dc_meshes) == 1
+    np.testing.assert_array_equal(default.dc_meshes[0].vertices, explicit.dc_meshes[0].vertices)
+    np.testing.assert_array_equal(default.dc_meshes[0].edges, explicit.dc_meshes[0].edges)
+    assert default.dc_meshes[0].capping_report is None
+    for first, second in zip(default.octrees_output, explicit.octrees_output):
+        np.testing.assert_array_equal(BackendTensor.t.to_numpy(first.grid.octree_grid.values),
+                                      BackendTensor.t.to_numpy(second.grid.octree_grid.values))
+        np.testing.assert_array_equal(BackendTensor.t.to_numpy(first.outputs[0].exported_fields.scalar_field),
+                                      BackendTensor.t.to_numpy(second.outputs[0].exported_fields.scalar_field))
+    boundary.assert_not_called()
+
+
+def test_enabled_preserves_physical_extent_and_caller_grid(plane_model, backend):
+    inputs, options, descriptor = plane_model()
+    grid = inputs.grid
+    root = grid.octree_grid
+    original_extent = BackendTensor.t.to_numpy(root.orthogonal_extent).copy()
+    original_values = BackendTensor.t.to_numpy(root.values).copy()
+    options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL
+    solution = compute_model(inputs, options, descriptor)
+    assert len(solution.octrees_output) == 2
+    root_values = BackendTensor.t.to_numpy(solution.octrees_output[0].grid.octree_grid.values)
+    assert root_values.dtype == np.float64
+    coordinate_tolerance = 2 * np.finfo(root_values.dtype).eps * np.max(np.abs(EXTENT))
+    for level in solution.octrees_output:
+        octree = level.grid.octree_grid
+        np.testing.assert_array_equal(BackendTensor.t.to_numpy(octree.physical_extent), EXTENT)
+        if backend is AvailableBackends.PYTORCH:
+            assert octree.physical_extent.device.type == "cpu"
+        coordinates = BackendTensor.t.to_numpy(octree.integer_coordinates)
+        shape = BackendTensor.t.to_numpy(octree.regular_grid_shape)
+        expected = EXTENT[::2] + (coordinates + .5) * (EXTENT[1::2] - EXTENT[::2]) / shape
+        np.testing.assert_allclose(BackendTensor.t.to_numpy(octree.values), expected,
+                                   rtol=0, atol=coordinate_tolerance)
+    np.testing.assert_array_equal(
+        BackendTensor.t.to_numpy(solution.octrees_output[0].grid.octree_grid.orthogonal_extent), EXTENT)
+    assert inputs.grid is grid
+    assert grid.octree_grid is root
+    np.testing.assert_array_equal(BackendTensor.t.to_numpy(root.physical_extent), EXTENT)
+    np.testing.assert_array_equal(BackendTensor.t.to_numpy(root.orthogonal_extent), original_extent)
+    np.testing.assert_array_equal(BackendTensor.t.to_numpy(root.values), original_values)
+
+
+def test_enabled_planes_are_closed_with_analytic_volume_and_reuse_boundary(plane_model, monkeypatch):
+    levels = (-4.7, -6.3)
+    inputs, options, descriptor = plane_model(levels)
+    options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL
+    boundary = Mock(wraps=dc._interp_on_boundary)
+    monkeypatch.setattr(dc, "_interp_on_boundary", boundary)
+    solution = compute_model(inputs, options, descriptor)
+    boundary.assert_called_once()
+    assert len(solution.dc_meshes) == len(levels)
+    volumes = []
+    for index, (mesh, height) in enumerate(zip(solution.dc_meshes, levels)):
+        report = mesh.capping_report
+        assert report["watertight"], report
+        assert report["closure_success"], report
+        assert not report["warnings"], report
+        assert report["open_edges_before"] > 0
+        assert report["open_edges_after"] == 0
+        assert report["cap_max_plane_error"] == 0
+        assert report["added_cap_triangles"] > 0
+        assert report["added_transition_triangles"] > 0
+        assert report["boundary_scalar_points"] == len(boundary.call_args.args[0])
+        assert (mesh.stack_index, mesh.surface_index, mesh.exported_surface_index) == (0, index, index)
+        assert mesh.isovalue == height
+        assert mesh.inside_convention == "scalar <= isovalue"
+        np.testing.assert_array_equal(mesh.vertices.min(axis=0), EXTENT[::2])
+        np.testing.assert_array_equal(mesh.vertices.max(axis=0)[:2], EXTENT[1:4:2])
+        assert np.all(mesh.vertices <= EXTENT[1::2])
+        triangles = mesh.vertices[mesh.edges]
+        edges = np.concatenate([mesh.edges[:, [0, 1]], mesh.edges[:, [1, 2]], mesh.edges[:, [2, 0]]])
+        _, counts = np.unique(np.sort(edges, axis=1), axis=0, return_counts=True)
+        np.testing.assert_array_equal(counts, 2)
+        # Shift the origin to avoid cancellation for a translated model box.
+        triangles = triangles - EXTENT[::2]
+        volume = np.einsum("ij,ij->i", triangles[:, 0], np.cross(triangles[:, 1], triangles[:, 2])).sum() / 6
+        volumes.append(volume)
+    np.testing.assert_allclose(volumes, 6 * 4 * (np.asarray(levels) + 8), rtol=1e-10, atol=0)
+
+
+@pytest.mark.parametrize("fail", [False, True], ids=["success", "failure"])
+def test_boundary_chunks_restore_grid_stack_and_options(plane_model, monkeypatch, fail):
+    inputs, options, _ = plane_model()
+    descriptor = SimpleNamespace(stack_structure=SimpleNamespace(n_stacks=2, stack_number=1))
+    saved_grid = inputs.grid
+    options.evaluation_options.evaluation_chunk_size = 3
+    options.evaluation_options.compute_scalar = False
+    options.evaluation_options.compute_scalar_gradient = True
+    saved_options = options.model_dump_json()
+    points = np.arange(21, dtype=float).reshape(7, 3)
+    seen = []
+
+    def interpolate(actual_input, boundary_options, actual_descriptor):
+        assert actual_input is inputs
+        assert actual_descriptor is descriptor
+        assert boundary_options is not options
+        assert boundary_options.evaluation_options.compute_scalar is True
+        assert boundary_options.evaluation_options.compute_scalar_gradient is False
+        xyz = actual_input.grid.custom_grid.values
+        seen.append(BackendTensor.t.to_numpy(xyz).copy())
+        descriptor.stack_structure.stack_number = 0
+        if fail and len(seen) == 2:
+            raise RuntimeError("boundary evaluation failed")
+        # Include sentinel values outside the custom-grid slice.
+        return [SimpleNamespace(
+            grid=SimpleNamespace(custom_grid_slice=slice(1, len(xyz) + 1)),
+            exported_fields=SimpleNamespace(scalar_field=BackendTensor.t.concatenate([
+                BackendTensor.t.array([-999.]), xyz[:, 2] + offset, BackendTensor.t.array([-999.])
+            ])),
+        ) for offset in (0, 100)]
+
+    interpolation = Mock(side_effect=interpolate)
+    monkeypatch.setattr(dc, "interpolate_all_fields_no_octree", interpolation)
+    if fail:
+        with pytest.raises(RuntimeError, match="boundary evaluation failed"):
+            dc._interp_on_boundary(points, inputs, options, descriptor)
+    else:
+        scalars = dc._interp_on_boundary(points, inputs, options, descriptor)
+        assert len(scalars) == 2
+        for actual, offset in zip(scalars, (0, 100)):
+            np.testing.assert_array_equal(actual, points[:, 2] + offset)
+    assert interpolation.call_count == (2 if fail else 3)
+    np.testing.assert_array_equal(np.concatenate(seen), points[:6] if fail else points)
+    assert inputs.grid is saved_grid
+    assert descriptor.stack_structure.stack_number == 1
+    assert options.model_dump_json() == saved_options
+
+
+def assert_report_matches_mesh(mesh, emitted):
+    """Best-effort closure must not conceal topology or geometry defects."""
+    report = mesh.capping_report
+    assert report is not None
+    assert np.isfinite(mesh.vertices).all()
+    assert mesh.edges.ndim == 2 and mesh.edges.shape[1] == 3
+    assert np.all((mesh.edges >= 0) & (mesh.edges < len(mesh.vertices)))
+    edges = np.concatenate([mesh.edges[:, [0, 1]], mesh.edges[:, [1, 2]], mesh.edges[:, [2, 0]]])
+    _, inverse, counts = np.unique(np.sort(edges, axis=1), axis=0, return_inverse=True, return_counts=True)
+    assert report["open_edges_after"] == np.count_nonzero(counts == 1), report
+    assert report["nonmanifold_edges"] == np.count_nonzero(counts > 2), report
+    directions = np.bincount(inverse, weights=np.where(edges[:, 0] < edges[:, 1], 1, -1))
+    assert report["orientation_conflicts"] == np.count_nonzero((counts == 2) & (np.abs(directions) == 2)), report
+    triangles = mesh.vertices[mesh.edges]
+    cross = np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0])
+    assert report["degenerate_triangles"] == np.count_nonzero(np.all(cross == 0, axis=1)), report
+    defects = ("missing_qef", "escaped_qef", "unsupported_topology", "boundary_edges",
+               "nonmanifold_edges", "orientation_conflicts", "nonmanifold_vertices",
+               "invalid_triangles", "nonfinite_vertices", "nonfinite_triangles", "degenerate_triangles")
+    for key in defects:
+        if report[key]:
+            assert not report["closure_success"], report
+            assert f"{key}: {report[key]}" in report["warnings"], report
+    if report["warnings"]:
+        assert any(issubclass(w.category, RuntimeWarning)
+                   and all(message in str(w.message) for message in report["warnings"])
+                   for w in emitted), report
+    if report["closure_success"]:
+        assert report["watertight"] and not report["warnings"], report
+    elif len(mesh.edges):
+        assert report["warnings"], report
+    if report["watertight"]:
+        assert len(counts) > 0
+        np.testing.assert_array_equal(counts, 2)
+    assert report["added_triangles"] == report["added_cap_triangles"] + report["added_transition_triangles"]
+    assert report["cap_max_plane_error"] == 0
+
+
+@pytest.mark.parametrize("normal,level,volume", [
+    ((-.2, -.1, 1.), .313, .463),
+    ((.2, .1, -1.), -.313, .537),
+])
+def test_oblique_plane_production_volume(plane_model, normal, level, volume):
+    inputs, options, descriptor = plane_model((level,), normal, (0, 1, 0, 1, 0, 1))
+    options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL
+    with warnings.catch_warnings(record=True) as emitted:
+        warnings.simplefilter("always")
+        solution = compute_model(inputs, options, descriptor)
+    assert len(solution.dc_meshes) == 1
+    mesh = solution.dc_meshes[0]
+    assert_report_matches_mesh(mesh, emitted)
+    assert mesh.capping_report["closure_success"], mesh.capping_report
+    original = mesh.vertices[:mesh.capping_report["original_vertex_count"]]
+    np.testing.assert_allclose(original @ normal, level, rtol=0, atol=2e-13)
+    assert np.all(mesh.vertices >= -2e-13) and np.all(mesh.vertices <= 1 + 2e-13)
+    triangles = mesh.vertices[mesh.edges]
+    actual_volume = np.einsum("ij,ij->i", triangles[:, 0], np.cross(triangles[:, 1], triangles[:, 2])).sum() / 6
+    assert actual_volume == pytest.approx(volume, rel=1e-11)
+
+
+@pytest.mark.parametrize("normal,level", [
+    ((0., 0., 1.), .5),
+    ((1., 1., 1.), 1.),
+    ((1., 1., 1.), 1.5),
+], ids=["lattice-plane", "box-corners", "lattice-points"])
+def test_exact_plane_ties_have_honest_diagnostics(plane_model, normal, level):
+    inputs, options, descriptor = plane_model((level,), normal, (0, 1, 0, 1, 0, 1), (2, 2, 2))
+    options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL
+    with warnings.catch_warnings(record=True) as emitted:
+        warnings.simplefilter("always")
+        solution = compute_model(inputs, options, descriptor)
+    assert len(solution.dc_meshes) == 1
+    mesh = solution.dc_meshes[0]
+    assert len(mesh.vertices) > 0 and len(mesh.edges) > 0
+    assert_report_matches_mesh(mesh, emitted)
+    report = mesh.capping_report
+    assert report["closure_success"] or report["warnings"], report
+    added = mesh.vertices[report["original_vertex_count"]:]
+    assert len(np.unique(added, axis=0)) == len(added)
+
+
+@pytest.mark.parametrize("height", [-9., -2.], ids=["fully-outside", "fully-inside"])
+def test_nonintersecting_plane_produces_no_artificial_box(plane_model, height):
+    inputs, options, descriptor = plane_model((height,))
+    options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL
+    with warnings.catch_warnings(record=True) as emitted:
+        warnings.simplefilter("always")
+        solution = compute_model(inputs, options, descriptor)
+    assert len(solution.dc_meshes) == 1
+    mesh = solution.dc_meshes[0]
+    assert mesh.vertices.shape == (0, 3)
+    assert mesh.edges.shape == (0, 3)
+    assert_report_matches_mesh(mesh, emitted)
+    assert mesh.capping_report["added_vertices"] == 0
+    assert mesh.capping_report["added_triangles"] == 0
+    assert not mesh.capping_report["closure_success"]
+
+
+@pytest.mark.parametrize("skip", ["1", "true"])
+def test_skip_triangulation_skips_capping(plane_model, monkeypatch, skip):
+    inputs, options, descriptor = plane_model()
+    options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL
+    monkeypatch.setenv("GEMPY_SKIP_TRIANGULATION", skip)
+    boundary = Mock(side_effect=AssertionError("Skipped triangulation must not evaluate caps"))
+    cap = Mock(side_effect=AssertionError("Skipped triangulation must not add triangles"))
+    monkeypatch.setattr(dc, "_interp_on_boundary", boundary)
+    monkeypatch.setattr(dc, "cap_mesh", cap)
+    solution = compute_model(inputs, options, descriptor)
+    assert len(solution.dc_meshes) == 1
+    mesh = solution.dc_meshes[0]
+    assert len(mesh.vertices) > 0
+    assert mesh.edges.size == 0
+    assert mesh.capping_report is None
+    boundary.assert_not_called()
+    cap.assert_not_called()
+
+
+@pytest.mark.parametrize("backend", [AvailableBackends.numpy], indirect=True)
+@pytest.mark.parametrize("normal,level", [((0., 0., 1.), .43), ((-.2, -.1, 1.), .313)])
+def test_numpy_pytorch_mesh_equivalence(plane_model, normal, level):
+    torch = pytest.importorskip("torch")
+    saved_grad = torch.is_grad_enabled()
+    meshes = []
+    try:
+        for engine in (AvailableBackends.numpy, AvailableBackends.PYTORCH):
+            BackendTensor._change_backend(engine, use_gpu=False, use_pykeops=False,
+                                          dtype="float64", grads=saved_grad)
+            inputs, options, descriptor = plane_model((level,), normal, (0, 1, 0, 1, 0, 1))
+            options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL
+            meshes.append(compute_model(inputs, options, descriptor).dc_meshes[0])
+        first, second = meshes
+        assert first.capping_report["closure_success"], first.capping_report
+        assert second.capping_report["closure_success"], second.capping_report
+        np.testing.assert_array_equal(first.edges, second.edges)
+        np.testing.assert_allclose(first.vertices, second.vertices, rtol=0, atol=2e-13)
+        assert first.capping_report == second.capping_report
+    finally:
+        torch.set_grad_enabled(saved_grad)
+
+
+def test_fault_model_capping_reports_all_limitations(graben_fault_model, monkeypatch):
+    inputs, descriptor, options = deepcopy(graben_fault_model)
+    options.evaluation_options.number_octree_levels = 3
+    options.evaluation_options.number_octree_levels_surface = 3
+    options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL
+    monkeypatch.setenv("GEMPY_SKIP_TRIANGULATION", "0")
+    with warnings.catch_warnings(record=True) as emitted:
+        warnings.simplefilter("always")
+        solution = compute_model(inputs, options, descriptor)
+    assert len(solution.dc_meshes) == 6
+    assert {mesh.stack_index for mesh in solution.dc_meshes} == {0, 1, 2}
+    for mesh in solution.dc_meshes:
+        assert_report_matches_mesh(mesh, emitted)
+        assert mesh.capping_report["boundary_scalar_points"] == 0
+        assert mesh.capping_report["added_triangles"] == 0
+        assert "boundary ownership" in mesh.capping_report["skipped_reason"]
+        assert mesh.inside_convention == "scalar <= isovalue"
+
+
+def test_capped_export_offsets_do_not_mutate_meshes(plane_model, monkeypatch):
+    raw_module = importlib.import_module("gempy_engine.core.data.raw_arrays_solution")
+    monkeypatch.setattr(raw_module, "require_subsurface", lambda: SimpleNamespace(
+        UnstructuredData=SimpleNamespace(from_array=lambda **kwargs: kwargs)))
+    monkeypatch.setattr(raw_module, "require_pandas", lambda: SimpleNamespace(DataFrame=lambda data: data))
+    inputs, options, descriptor = plane_model((-4.7, -6.3))
+    options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL
+    meshes = compute_model(inputs, options, descriptor).dc_meshes
+    raw = raw_module.RawArraysSolution(vertices=[mesh.vertices for mesh in meshes],
+                                       edges=[mesh.edges for mesh in meshes])
+    original = [mesh.edges.copy() for mesh in meshes]
+    expected = np.concatenate((original[0], original[1] + len(meshes[0].vertices)))
+    for _ in range(2):
+        exported = raw.meshes_to_subsurface()
+        np.testing.assert_array_equal(exported["cells"], expected)
+        for mesh, triangles in zip(meshes, original):
+            np.testing.assert_array_equal(mesh.edges, triangles)
diff --git a/tests/test_common/test_modules/test_scalar_crossing.py b/tests/test_common/test_modules/test_scalar_crossing.py
new file mode 100644
index 00000000..e892f70d
--- /dev/null
+++ b/tests/test_common/test_modules/test_scalar_crossing.py
@@ -0,0 +1,111 @@
+import itertools
+
+import numpy as np
+import pytest
+
+from gempy_engine.core.backend_tensor import BackendTensor, AvailableBackends
+from gempy_engine.core.data.dual_contouring_data import DualContouringData
+from gempy_engine.modules.dual_contouring._gen_vertices import generate_dual_contouring_vertices
+from gempy_engine.modules.dual_contouring._scalar_crossing import scalar_crossing_parameters
+
+
+def test_crossing_snaps_near_endpoints():
+    crossing, parameters = scalar_crossing_parameters(
+        np.array([-1e-14, -1.]), np.array([1., 1e-14]), 0.
+    )
+    np.testing.assert_array_equal(crossing, [True, True])
+    np.testing.assert_array_equal(parameters, [0., 1.])
+from gempy_engine.modules.dual_contouring.dual_contouring_interface import find_intersection_on_edge
+
+
+@pytest.fixture(params=[AvailableBackends.numpy, AvailableBackends.PYTORCH])
+def backend(request):
+    if request.param == AvailableBackends.PYTORCH:
+        pytest.importorskip("torch")
+    BackendTensor._change_backend(request.param, use_gpu=False)
+    return BackendTensor.t
+
+
+def as_numpy(value):
+    return value.detach().cpu().numpy() if hasattr(value, "detach") else value
+
+
+def test_parameters(backend):
+    start = backend.array([-1., 1., 0., 1., 0., -1., 0., 1., -1., -1e-14])
+    end = backend.array([1., -1., 1., 0., -1., 0., 0., 2., -1., 3e-14])
+    iso = backend.array(0.)
+    valid, t = scalar_crossing_parameters(start, end, iso, xp=backend)
+    np.testing.assert_array_equal(as_numpy(valid), [1, 1, 1, 1, 0, 0, 0, 0, 0, 1])
+    np.testing.assert_allclose(as_numpy(t), [.5, .5, 0, 1, 0, 0, 0, 0, 0, .25])
+    reverse_valid, reverse_t = scalar_crossing_parameters(end, start, iso, xp=backend)
+    np.testing.assert_array_equal(as_numpy(reverse_valid), as_numpy(valid))
+    np.testing.assert_allclose(as_numpy(reverse_t)[as_numpy(valid)], 1 - as_numpy(t)[as_numpy(valid)])
+
+
+@pytest.mark.parametrize("bad", [np.nan, np.inf, -np.inf])
+@pytest.mark.parametrize("argument", [0, 1, 2])
+def test_nonfinite(backend, bad, argument):
+    args = [backend.array([-1.]), backend.array([1.]), backend.array([0.])]
+    args[argument] = backend.array([bad])
+    with pytest.raises(ValueError, match="finite"):
+        scalar_crossing_parameters(*args, xp=backend)
+
+
+def test_numpy_utility_independent_of_backend(backend):
+    valid, t = scalar_crossing_parameters(np.array([-1.]), np.array([3.]), np.array([0.]))
+    assert isinstance(t, np.ndarray)
+    np.testing.assert_array_equal(valid, [True])
+    np.testing.assert_array_equal(t, [.25])
+
+
+def test_intersections_masking_and_multiple_surfaces(backend):
+    cube = np.array(list(itertools.product([0., 1.], repeat=3)))
+    xyz = backend.array(np.concatenate([cube, cube + 2]))
+    scalars = xyz[:, 0]
+    iso = backend.array([0., .25, 1.])
+    points, valid = find_intersection_on_edge(
+        xyz, scalars, iso, masking=backend.array([True, False]), strict_crossings=True
+    )
+    expected_mask = np.zeros((3, 12), dtype=bool)
+    expected_mask[:2, :4] = True
+    np.testing.assert_array_equal(as_numpy(valid).reshape(3, 12), expected_mask)
+    expected = np.concatenate([cube[:4], cube[:4] + [.25, 0, 0]])
+    np.testing.assert_allclose(as_numpy(points), expected)
+    assert valid.ndim == (1 if BackendTensor.engine_backend == AvailableBackends.PYTORCH else 2)
+
+
+def test_empty_and_iso_edges(backend):
+    xyz = backend.array(np.zeros((8, 3)))
+    for masking in (None, backend.array([False])):
+        points, valid = find_intersection_on_edge(
+            xyz, xyz[:, 0], backend.array([0.]), masking=masking, strict_crossings=True
+        )
+        assert points.shape == (0, 3)
+        assert not as_numpy(valid).any()
+
+
+def test_default_retains_tolerant_crossings(backend):
+    xyz = backend.array(np.array(list(itertools.product([0., 1.], repeat=3))))
+    iso = backend.array([1.005])
+    default = find_intersection_on_edge(xyz, xyz[:, 0], iso)
+    disabled = find_intersection_on_edge(xyz, xyz[:, 0], iso, strict_crossings=False)
+    for actual, expected in zip(default, disabled):
+        np.testing.assert_array_equal(as_numpy(actual), as_numpy(expected))
+    assert as_numpy(default[1]).any()
+    strict = find_intersection_on_edge(xyz, xyz[:, 0], iso, strict_crossings=True)
+    assert not as_numpy(strict[1]).any()
+
+
+def test_mass_point_includes_zero_coordinates(backend):
+    valid = np.zeros((2, 12), dtype=bool)
+    valid[0, :2] = True
+    data = DualContouringData(
+        xyz_on_edge=backend.array([[0., 0., 0.], [2., 4., 0.]], dtype=BackendTensor.dtype_obj),
+        valid_edges=backend.array(valid), xyz_on_centers=None, dxdydz=(1., 1., 1.),
+        n_surfaces_to_export=1, left_right_codes=None,
+        gradients=backend.array(np.zeros((2, 3)), dtype=BackendTensor.dtype_obj), strict_crossings=True
+    )
+    vertices = generate_dual_contouring_vertices(data, debug=True)
+    np.testing.assert_allclose(as_numpy(data.bias_center_mass), [[1., 2., 0.]] * 3)
+    assert np.isfinite(as_numpy(vertices)).all()
+    assert DualContouringData.__dataclass_fields__["strict_crossings"].default is False

From 524e4b8200ceee169dac0b72ed2280ecc3d6b6ea Mon Sep 17 00:00:00 2001
From: Miguel de la Varga <miguel@terranigma-solutions.com>
Date: Mon, 21 Sep 2026 17:05:39 +0200
Subject: [PATCH 2/2] [ENH] Add modular extent capping implementation with
 tests and documentation updates

Refactor extent capping logic into a reusable function for dual-contouring workflows. Add integration tests ensuring consistency across modes and improve orthogonal extent handling. Update documentation to reflect boundary evaluation mechanics and configuration options.
---
 docs/extent_capping.md                        | 31 ++++----
 .../API/dual_contouring/extent_capping.py     | 76 +++++++++++++++++++
 .../multi_scalar_dual_contouring.py           | 64 ++--------------
 gempy_engine/API/model/model_api.py           | 20 -----
 gempy_engine/core/data/engine_grid.py         |  4 +-
 gempy_engine/core/data/regular_grid.py        | 11 +--
 .../modules/dual_contouring/_gen_vertices.py  | 25 ++----
 .../test_extent_capping_integration.py        | 13 ++--
 .../test_modules/test_grid_zero_handling.py   | 67 ++++++++++++++++
 9 files changed, 185 insertions(+), 126 deletions(-)
 create mode 100644 gempy_engine/API/dual_contouring/extent_capping.py
 create mode 100644 tests/test_common/test_modules/test_grid_zero_handling.py

diff --git a/docs/extent_capping.md b/docs/extent_capping.md
index 4e948132..5a7d9860 100644
--- a/docs/extent_capping.md
+++ b/docs/extent_capping.md
@@ -33,22 +33,25 @@ If no isosurface intersects the box, no box shell is manufactured, even when
 the whole box satisfies the inside condition. This is closure of existing
 isosurfaces, not unconditional extraction of the complete clipped sublevel set.
 
-`RegularGrid.physical_extent` retains the requested extent through refinement.
-Enabled `compute_model` calls use a private root grid sampled against that exact
-box, rather than the legacy `1e-6` translated box. Caps are generated in engine
-coordinates, before downstream output transforms. Original vertex indices are
-preserved and cap vertices are appended.
+`RegularGrid.orthogonal_extent` is the single requested extent used for grid
+sampling, refinement, and capping. All octree levels retain that exact extent,
+whether capping is enabled or disabled. There is no conditional private root-grid
+resampling or extent translation. Caps are generated in engine coordinates,
+before downstream output transforms. Original vertex indices are preserved and
+cap vertices are appended.
 
 Enabled extraction and cap clipping share strict scalar crossing rules. An
 endpoint equal to the isovalue is inside; an entire iso-valued edge is not a
 crossing. Interpolation parameters are clamped and snapped within `1e-12` of an
-endpoint. Nonfinite scalar inputs are rejected with a diagnostic. Enabled QEF
-mass points use edge-validity masks, including genuine zero coordinates, and
-avoid PyTorch's additional origin-centered regularization.
+endpoint. Nonfinite scalar inputs are rejected with a diagnostic. QEF mass
+points use edge-validity masks globally, including genuine zero coordinates,
+regardless of capping mode. Enabled QEF solving avoids PyTorch's additional
+origin-centered regularization.
 
-With capping disabled, legacy grid sampling, crossings, QEF solving, and mesh
-arrays retain their previous behavior. Enabled interior geometry can differ
-because its sampling and crossing conventions are deliberately stricter.
+With capping disabled, no boundary evaluation or cap construction is performed.
+Exact grid sampling and zero-coordinate handling apply in both modes; uncapped
+mesh arrays are not guaranteed to match legacy output. Enabled interior geometry
+can differ because its crossing and QEF-solving conventions are stricter.
 
 ## Best-Effort Contract
 
@@ -99,11 +102,13 @@ cap geometry and topology are NumPy postprocessing and are not differentiable.
 
 One unique finest-level six-face lattice is shared across surfaces. Boundary
 scalar values are evaluated once per stack per batch and reused across that
-stack's isosurfaces. `evaluation_chunk_size` bounds the boundary point batch;
+stack's isosurfaces. The API module `API/dual_contouring/extent_capping.py`
+orchestrates capping and restores the grid and stack cursor after boundary
+evaluation, including failures. `evaluation_chunk_size` bounds the boundary point batch;
 the existing evaluator also applies its own kernel-workload chunking.
 
 Boundary storage scales as `O(nx*ny + nx*nz + ny*nz)`, but Python dictionaries,
 triangle connectivity, and topology audits have substantial additional memory
 cost. This is a correctness-first implementation: there is no adaptive cap
-simplification, cross-call cache, or reuse of shifted legacy corner samples.
+simplification, cross-call cache, or reuse of existing corner samples.
 Large-depth GPU and production-scale memory benchmarks remain outstanding.
diff --git a/gempy_engine/API/dual_contouring/extent_capping.py b/gempy_engine/API/dual_contouring/extent_capping.py
new file mode 100644
index 00000000..de205e65
--- /dev/null
+++ b/gempy_engine/API/dual_contouring/extent_capping.py
@@ -0,0 +1,76 @@
+"""Orchestrate extent capping and batched boundary scalar evaluation."""
+
+import copy
+import os
+
+import numpy as np
+
+from ..interp_single.interp_features import interpolate_all_fields_no_octree
+from ...core.backend_tensor import BackendTensor
+from ...core.data.engine_grid import EngineGrid
+from ...core.data.generic_grid import GenericGrid
+from ...core.data.stack_relation_type import StackRelationType
+from ...modules.dual_contouring._extent_capping import boundary_lattice, cap_mesh
+
+
+def cap_meshes_at_extent(all_meshes, dc_data_per_surface, all_mask_arrays, base_number,
+                         orthogonal_extent, interpolation_input, options, data_descriptor):
+    """Cap extracted meshes in place, sharing boundary samples across surfaces."""
+    skip_triangles = os.getenv("GEMPY_SKIP_TRIANGULATION", "0").lower() in ("true", "1", "t", "y", "yes")
+    if not all_meshes or skip_triangles:
+        return
+
+    extent = BackendTensor.t.to_numpy(orthogonal_extent)
+    skip_reasons = []
+    stack_relations = data_descriptor.stack_structure.masking_descriptor
+    fault_relations = data_descriptor.stack_structure.faults_relations
+    for mesh in all_meshes:
+        stack = mesh.stack_index
+        mask = all_mask_arrays[stack]
+        faulted = (stack_relations[stack] is StackRelationType.FAULT
+                   or (fault_relations is not None and np.any(fault_relations[:, stack])))
+        masked = mask is not None and not bool(mask.all())
+        skip_reasons.append(
+            "Capping skipped: fault or extraction-mask boundary ownership is not supported"
+            if faulted or masked else None
+        )
+    geometry = None
+    coordinates = points = np.empty((0, 3))
+    boundary_scalars = None
+    if any(reason is None for reason in skip_reasons):
+        geometry = boundary_lattice(base_number, extent)
+        coordinates, points, _ = geometry
+        boundary_scalars = _interp_on_boundary(points, interpolation_input, options, data_descriptor)
+    for mesh, dc_data, reason in zip(all_meshes, dc_data_per_surface, skip_reasons):
+        cell_coordinates = BackendTensor.t.to_numpy(dc_data.left_right_codes[dc_data.valid_voxels])
+        cap_mesh(mesh, cell_coordinates, base_number, extent, coordinates,
+                 boundary_scalars[mesh.stack_index] if reason is None else np.empty(0),
+                 mesh.isovalue, boundary_geometry=geometry, skip_reason=reason)
+        mesh.capping_report["boundary_scalar_points"] = len(points) if reason is None else 0
+        mesh.capping_report["original_vertex_count"] = len(cell_coordinates)
+
+
+def _interp_on_boundary(points, interpolation_input, options, data_descriptor):
+    """Evaluate each stack once per boundary batch, reusing it for its surfaces."""
+    boundary_options = copy.deepcopy(options)
+    boundary_options.evaluation_options.compute_scalar = True
+    boundary_options.evaluation_options.compute_scalar_gradient = False
+    saved_grid = interpolation_input.grid
+    saved_stack = data_descriptor.stack_structure.stack_number
+    scalars = [np.empty(len(points)) for _ in range(data_descriptor.stack_structure.n_stacks)]
+    batch_size = max(1, int(options.evaluation_options.evaluation_chunk_size))
+    try:
+        for start in range(0, len(points), batch_size):
+            stop = min(start + batch_size, len(points))
+            interpolation_input.set_temp_grid(EngineGrid(custom_grid=GenericGrid(
+                values=BackendTensor.t.array(points[start:stop], dtype=BackendTensor.dtype)
+            )))
+            outputs = interpolate_all_fields_no_octree(interpolation_input, boundary_options, data_descriptor)
+            for stack_index, output in enumerate(outputs):
+                scalars[stack_index][start:stop] = BackendTensor.t.to_numpy(
+                    output.exported_fields.scalar_field[output.grid.custom_grid_slice]
+                )
+    finally:
+        interpolation_input.set_temp_grid(saved_grid)
+        data_descriptor.stack_structure.stack_number = saved_stack
+    return scalars
diff --git a/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py b/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py
index 05ddc042..c5f0b12c 100644
--- a/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py
+++ b/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py
@@ -1,5 +1,4 @@
 import copy
-import os
 import warnings
 from typing import List, Any
 
@@ -26,7 +25,7 @@
 from ...modules.dual_contouring.overlapping import average_overlapping_vertices, remove_fault_overlap_triangles
 from ...modules.dual_contouring._support_report import mesh_support_report
 from ...core.data.options.evaluation_options import OctreeRefinementMode, MeshExtentCapping
-from ...modules.dual_contouring._extent_capping import boundary_lattice, cap_mesh
+from .extent_capping import cap_meshes_at_extent
 
 
 @gempy_profiler_decorator
@@ -228,65 +227,16 @@ def dual_contouring_multi_scalar(
         mesh.isovalue = isovalue
         mesh.inside_convention = "scalar <= isovalue" if cap_enabled else None
 
-    skip_triangles = os.getenv("GEMPY_SKIP_TRIANGULATION", "0").lower() in ("true", "1", "t", "y", "yes")
-    if cap_enabled and all_meshes and not skip_triangles:
-        extent = BackendTensor.t.to_numpy(octree_list[0].grid.octree_grid.physical_extent)
-        skip_reasons = []
-        fault_relations = data_descriptor.stack_structure.faults_relations
-        for mesh in all_meshes:
-            stack = mesh.stack_index
-            mask = all_mask_arrays[stack]
-            faulted = (stack_relations[stack] is StackRelationType.FAULT
-                       or (fault_relations is not None and np.any(fault_relations[:, stack])))
-            masked = mask is not None and not bool(mask.all())
-            skip_reasons.append(
-                "Capping skipped: fault or extraction-mask boundary ownership is not supported"
-                if faulted or masked else None
-            )
-        geometry = None
-        coordinates = points = np.empty((0, 3))
-        boundary_scalars = None
-        if any(reason is None for reason in skip_reasons):
-            geometry = boundary_lattice(base_number, extent)
-            coordinates, points, _ = geometry
-            boundary_scalars = _interp_on_boundary(points, interpolation_input, options, data_descriptor)
-        for mesh, dc_data, reason in zip(all_meshes, dc_data_per_surface_all, skip_reasons):
-            cell_coordinates = BackendTensor.t.to_numpy(dc_data.left_right_codes[dc_data.valid_voxels])
-            cap_mesh(mesh, cell_coordinates, base_number, extent, coordinates,
-                     boundary_scalars[mesh.stack_index] if reason is None else np.empty(0),
-                     mesh.isovalue, boundary_geometry=geometry, skip_reason=reason)
-            mesh.capping_report["boundary_scalar_points"] = len(points) if reason is None else 0
-            mesh.capping_report["original_vertex_count"] = len(cell_coordinates)
+    if cap_enabled:
+        cap_meshes_at_extent(
+            all_meshes, dc_data_per_surface_all, all_mask_arrays, base_number,
+            octree_list[0].grid.octree_grid.orthogonal_extent,
+            interpolation_input, options, data_descriptor
+        )
 
     return all_meshes
 
 
-def _interp_on_boundary(points, interpolation_input, options, data_descriptor):
-    """Evaluate each stack once per boundary batch, reusing it for its surfaces."""
-    boundary_options = copy.deepcopy(options)
-    boundary_options.evaluation_options.compute_scalar = True
-    boundary_options.evaluation_options.compute_scalar_gradient = False
-    saved_grid = interpolation_input.grid
-    saved_stack = data_descriptor.stack_structure.stack_number
-    scalars = [np.empty(len(points)) for _ in range(data_descriptor.stack_structure.n_stacks)]
-    batch_size = max(1, int(options.evaluation_options.evaluation_chunk_size))
-    try:
-        for start in range(0, len(points), batch_size):
-            stop = min(start + batch_size, len(points))
-            interpolation_input.set_temp_grid(EngineGrid(custom_grid=GenericGrid(
-                values=BackendTensor.t.array(points[start:stop], dtype=BackendTensor.dtype)
-            )))
-            outputs = interpolate_all_fields_no_octree(interpolation_input, boundary_options, data_descriptor)
-            for stack_index, output in enumerate(outputs):
-                scalars[stack_index][start:stop] = BackendTensor.t.to_numpy(
-                    output.exported_fields.scalar_field[output.grid.custom_grid_slice]
-                )
-    finally:
-        interpolation_input.set_temp_grid(saved_grid)
-        data_descriptor.stack_structure.stack_number = saved_stack
-    return scalars
-
-
 def _validate_stack_relations(data_descriptor: InputDataDescriptor, n_scalar_field: int) -> None:
     """
     Validate stack relations for the given scalar field.
diff --git a/gempy_engine/API/model/model_api.py b/gempy_engine/API/model/model_api.py
index 033bb4ab..89cb670a 100644
--- a/gempy_engine/API/model/model_api.py
+++ b/gempy_engine/API/model/model_api.py
@@ -20,7 +20,6 @@
 from ...core.utils import gempy_profiler_decorator
 from ...core.exceptions import GemPyEngineInputError
 from ...core.data.options.temp_interpolation_values import TempInterpolationValues
-from ...core.data.options.evaluation_options import MeshExtentCapping
 from ...modules.geophysics.fw_gravity import compute_gravity
 from ...modules.geophysics.fw_magnetic import compute_tmi
 from ...modules.weights_cache.weights_cache_interface import WeightCache
@@ -44,25 +43,6 @@ def compute_model(interpolation_input: InterpolationInput, options: Interpolatio
         # Check input is valid
         _check_input_validity(interpolation_input, options, data_descriptor)  # TODO
 
-        if (options.evaluation_options.mesh_extraction
-                and MeshExtentCapping(options.evaluation_options.mesh_extraction_extent_capping) != MeshExtentCapping.NONE
-                and interpolation_input.grid.octree_grid is not None):
-            # Keep the legacy offset and caller-owned grids untouched when capping
-            # is disabled. Enabled extraction samples the exact physical box.
-            interpolation_input = copy.copy(interpolation_input)
-            grid = copy.copy(interpolation_input.grid)
-            root = copy.copy(grid.octree_grid)
-            root.orthogonal_extent = BackendTensor.t.copy(root.physical_extent)
-            coordinates = BackendTensor.t.array(root.integer_coordinates, dtype=BackendTensor.dtype)
-            root.values = root.physical_extent[::2] + (coordinates + 0.5) * (
-                (root.physical_extent[1::2] - root.physical_extent[::2]) / root.regular_grid_shape
-            )
-            root.original_values = BackendTensor.t.copy(root.values)
-            grid.octree_grid = root
-            grid.corners_grid = None
-            interpolation_input._original_grid = grid
-            interpolation_input.set_temp_grid(grid)
-
         output: list[OctreeLevel] = interpolate_n_octree_levels(
             interpolation_input=interpolation_input,
             options=options,
diff --git a/gempy_engine/core/data/engine_grid.py b/gempy_engine/core/data/engine_grid.py
index 8b42168b..1806ccc5 100644
--- a/gempy_engine/core/data/engine_grid.py
+++ b/gempy_engine/core/data/engine_grid.py
@@ -59,12 +59,10 @@ def from_xyz_coords(cls, xyz_coords: ndarray) -> "EngineGrid":
 
     @classmethod
     def from_regular_grid(cls, regular_grid: RegularGrid) -> "EngineGrid":
-        grid = cls(
+        return cls(
             dense_grid=regular_grid,
             octree_grid=RegularGrid(regular_grid.orthogonal_extent, np.array([2, 2, 2]))
         )
-        grid.octree_grid.physical_extent = BackendTensor.t.copy(regular_grid.physical_extent)
-        return grid
 
     @property
     def values(self) -> np.ndarray:
diff --git a/gempy_engine/core/data/regular_grid.py b/gempy_engine/core/data/regular_grid.py
index d2f7bb09..eb510785 100644
--- a/gempy_engine/core/data/regular_grid.py
+++ b/gempy_engine/core/data/regular_grid.py
@@ -18,7 +18,6 @@ class RegularGrid:
     left_right: np.ndarray = field(default=None, repr=False, init=False)
     _integer_coordinates: np.ndarray = field(default=None, repr=False, init=False)
     refinement_debug: dict = field(default=None, repr=False, init=False)
-    physical_extent: np.ndarray = field(default=None, repr=False, init=False)
     
     values: np.ndarray = field(default=None, repr=False, init=False)
     original_values: np.ndarray = field(default=None, repr=False, init=False)  #: When the regular grid is representing a octree level, only active cells are stored in values. This is the original values of the regular grid.
@@ -29,8 +28,7 @@ def __len__(self):
 
     def __post_init__(self):
         self.regular_grid_shape = BackendTensor.t.array(self.regular_grid_shape)
-        self.physical_extent = BackendTensor.t.array(self.orthogonal_extent)
-        self.orthogonal_extent = BackendTensor.t.array(self.orthogonal_extent) + 1e-6  # * This to avoid some errors evaluating in 0 (e.g. bias in dual contouring)
+        self.orthogonal_extent = BackendTensor.t.array(self.orthogonal_extent)
 
         self._create_regular_grid_3d()
 
@@ -50,15 +48,15 @@ def dz(self):
     
     @property
     def x_coord(self):
-        return BackendTensor.t.linspace(self.orthogonal_extent[0] + self.dx / 2, self.orthogonal_extent[1] - self.dx / 2, self.resolution[0])
+        return BackendTensor.t.linspace(self.orthogonal_extent[0] + self.dx / 2, self.orthogonal_extent[1] - self.dx / 2, self.resolution[0], dtype=BackendTensor.dtype_obj)
 
     @property
     def y_coord(self):
-        return BackendTensor.t.linspace(self.orthogonal_extent[2] + self.dy / 2, self.orthogonal_extent[3] - self.dy / 2, self.resolution[1])
+        return BackendTensor.t.linspace(self.orthogonal_extent[2] + self.dy / 2, self.orthogonal_extent[3] - self.dy / 2, self.resolution[1], dtype=BackendTensor.dtype_obj)
 
     @property
     def z_coord(self):
-        return BackendTensor.t.linspace(self.orthogonal_extent[4] + self.dz / 2, self.orthogonal_extent[5] - self.dz / 2, self.resolution[2])
+        return BackendTensor.t.linspace(self.orthogonal_extent[4] + self.dz / 2, self.orthogonal_extent[5] - self.dz / 2, self.resolution[2], dtype=BackendTensor.dtype_obj)
 
     @classmethod
     def from_octree_level(cls, xyz_coords_octree: np.ndarray, previous_regular_grid: "RegularGrid",
@@ -70,7 +68,6 @@ def from_octree_level(cls, xyz_coords_octree: np.ndarray, previous_regular_grid:
         )
 
         regular_grid_for_octree_level.values = xyz_coords_octree  # ! Overwrite the common values
-        regular_grid_for_octree_level.physical_extent = BackendTensor.t.copy(previous_regular_grid.physical_extent)
         regular_grid_for_octree_level._active_cells = active_cells
         regular_grid_for_octree_level.left_right = left_right
         regular_grid_for_octree_level._integer_coordinates = (
diff --git a/gempy_engine/modules/dual_contouring/_gen_vertices.py b/gempy_engine/modules/dual_contouring/_gen_vertices.py
index dca24d13..a0378975 100644
--- a/gempy_engine/modules/dual_contouring/_gen_vertices.py
+++ b/gempy_engine/modules/dual_contouring/_gen_vertices.py
@@ -29,27 +29,12 @@ def generate_dual_contouring_vertices(dc_data_per_stack: DualContouringData, sli
     edges_xyz[:, :12][valid_edges_bool] = xyz_on_edge
     edges_normals[:, :12][valid_edges_bool] = gradients
 
-    # Use nanmean directly without intermediate copy
+    # Zero coordinates are valid samples, not missing edge constraints.
     bias_xyz_slice = edges_xyz[:, :12]
-    
-    if dc_data_per_stack.strict_crossings:
-        # Zero coordinates are valid samples, not missing edge constraints.
-        mask = valid_edges_bool[:, :, None]
-        sum_valid = (bias_xyz_slice * mask).sum(axis=1)
-        count_valid = mask.sum(axis=1)
-        mass_points = sum_valid / count_valid
-    elif BackendTensor.engine_backend == AvailableBackends.PYTORCH:
-        mask = bias_xyz_slice == 0
-        bias_xyz_masked = BackendTensor.tfnp.where(mask, float('nan'), bias_xyz_slice)
-        mass_points = BackendTensor.tfnp.nanmean(bias_xyz_masked, axis=1)
-    else:
-        # NumPy: more efficient approach using sum and count
-        mask = bias_xyz_slice != 0
-        sum_valid = (bias_xyz_slice * mask).sum(axis=1)
-        count_valid = mask.sum(axis=1)
-        # Avoid division by zero
-        count_valid = BackendTensor.tfnp.maximum(count_valid, 1)
-        mass_points = sum_valid / count_valid
+    mask = valid_edges_bool[:, :, None]
+    sum_valid = (bias_xyz_slice * mask).sum(axis=1)
+    count_valid = mask.sum(axis=1)
+    mass_points = sum_valid / count_valid
 
     # Assign mass points to bias positions
     edges_xyz[:, 12:15] = mass_points[:, None, :]
diff --git a/tests/test_common/test_modules/test_extent_capping_integration.py b/tests/test_common/test_modules/test_extent_capping_integration.py
index 4dcc4c21..9eaa5dca 100644
--- a/tests/test_common/test_modules/test_extent_capping_integration.py
+++ b/tests/test_common/test_modules/test_extent_capping_integration.py
@@ -27,7 +27,7 @@
 from gempy_engine.core.data.stacks_structure import StacksStructure
 
 
-dc = importlib.import_module("gempy_engine.API.dual_contouring.multi_scalar_dual_contouring")
+dc = importlib.import_module("gempy_engine.API.dual_contouring.extent_capping")
 EXTENT = np.array([-2., 4., 10., 14., -8., -3.])
 
 
@@ -117,13 +117,14 @@ def test_default_and_explicit_none_have_identical_arrays(plane_model, monkeypatc
     boundary.assert_not_called()
 
 
-def test_enabled_preserves_physical_extent_and_caller_grid(plane_model, backend):
+@pytest.mark.parametrize("mode", list(MeshExtentCapping))
+def test_all_modes_preserve_exact_extent_and_caller_grid(plane_model, backend, mode):
     inputs, options, descriptor = plane_model()
     grid = inputs.grid
     root = grid.octree_grid
     original_extent = BackendTensor.t.to_numpy(root.orthogonal_extent).copy()
     original_values = BackendTensor.t.to_numpy(root.values).copy()
-    options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL
+    options.evaluation_options.mesh_extraction_extent_capping = mode
     solution = compute_model(inputs, options, descriptor)
     assert len(solution.octrees_output) == 2
     root_values = BackendTensor.t.to_numpy(solution.octrees_output[0].grid.octree_grid.values)
@@ -131,9 +132,9 @@ def test_enabled_preserves_physical_extent_and_caller_grid(plane_model, backend)
     coordinate_tolerance = 2 * np.finfo(root_values.dtype).eps * np.max(np.abs(EXTENT))
     for level in solution.octrees_output:
         octree = level.grid.octree_grid
-        np.testing.assert_array_equal(BackendTensor.t.to_numpy(octree.physical_extent), EXTENT)
+        np.testing.assert_array_equal(BackendTensor.t.to_numpy(octree.orthogonal_extent), EXTENT)
         if backend is AvailableBackends.PYTORCH:
-            assert octree.physical_extent.device.type == "cpu"
+            assert octree.orthogonal_extent.device.type == "cpu"
         coordinates = BackendTensor.t.to_numpy(octree.integer_coordinates)
         shape = BackendTensor.t.to_numpy(octree.regular_grid_shape)
         expected = EXTENT[::2] + (coordinates + .5) * (EXTENT[1::2] - EXTENT[::2]) / shape
@@ -143,7 +144,7 @@ def test_enabled_preserves_physical_extent_and_caller_grid(plane_model, backend)
         BackendTensor.t.to_numpy(solution.octrees_output[0].grid.octree_grid.orthogonal_extent), EXTENT)
     assert inputs.grid is grid
     assert grid.octree_grid is root
-    np.testing.assert_array_equal(BackendTensor.t.to_numpy(root.physical_extent), EXTENT)
+    np.testing.assert_array_equal(BackendTensor.t.to_numpy(root.orthogonal_extent), EXTENT)
     np.testing.assert_array_equal(BackendTensor.t.to_numpy(root.orthogonal_extent), original_extent)
     np.testing.assert_array_equal(BackendTensor.t.to_numpy(root.values), original_values)
 
diff --git a/tests/test_common/test_modules/test_grid_zero_handling.py b/tests/test_common/test_modules/test_grid_zero_handling.py
new file mode 100644
index 00000000..3a852c9d
--- /dev/null
+++ b/tests/test_common/test_modules/test_grid_zero_handling.py
@@ -0,0 +1,67 @@
+import numpy as np
+import pytest
+
+from gempy_engine.config import AvailableBackends
+from gempy_engine.core.backend_tensor import BackendTensor
+from gempy_engine.core.data.dual_contouring_data import DualContouringData
+from gempy_engine.core.data.engine_grid import EngineGrid
+from gempy_engine.core.data.regular_grid import RegularGrid
+from gempy_engine.modules.dual_contouring._gen_vertices import generate_dual_contouring_vertices
+from gempy_engine.modules.octrees_topology._octree_common import _generate_next_level_centers
+
+
+@pytest.fixture(params=['numpy', 'PYTORCH'])
+def backend(request):
+    if request.param == 'PYTORCH':
+        pytest.importorskip('torch')
+    old = BackendTensor.engine_backend, BackendTensor.use_gpu, BackendTensor.dtype, BackendTensor.use_pykeops
+    BackendTensor._change_backend(engine_backend=AvailableBackends[request.param], use_gpu=False, dtype='float64')
+    yield BackendTensor.t
+    BackendTensor._change_backend(engine_backend=old[0], use_gpu=old[1], dtype=old[2], use_pykeops=old[3])
+
+
+def test_grid_extent_and_centers_are_not_translated(backend):
+    t = backend
+    extent = [-1., 1., -2., 2., -3., 3.]
+    dense = RegularGrid(extent, [3, 3, 3])
+    np.testing.assert_array_equal(t.to_numpy(dense.values[13]), [0., 0., 0.])
+    grid = EngineGrid.from_regular_grid(dense)
+    assert grid.dense_grid is dense
+    root = grid.octree_grid
+    np.testing.assert_array_equal(t.to_numpy(root.orthogonal_extent), extent)
+    np.testing.assert_array_equal(t.to_numpy(dense.orthogonal_extent), extent)
+    np.testing.assert_array_equal(t.to_numpy(root.values[0]), [-0.5, -1., -1.5])
+    for _ in range(3):
+        active = root.integer_coordinates[:, 0] == 0
+        xyz, bits = _generate_next_level_centers(root.values[active], root.dxdydz)
+        root = RegularGrid.from_octree_level(xyz, root, active, bits)
+        np.testing.assert_array_equal(t.to_numpy(root.orthogonal_extent), extent)
+        expected = root.orthogonal_extent[::2] + (root.integer_coordinates + 0.5) * (
+            (root.orthogonal_extent[1::2] - root.orthogonal_extent[::2]) / root.regular_grid_shape
+        )
+        np.testing.assert_array_equal(t.to_numpy(root.values), t.to_numpy(expected))
+
+
+@pytest.mark.parametrize('strict_crossings', [False, True])
+@pytest.mark.parametrize('points', [
+    [[0., 0., 0.], [2., 0., 4.]],
+    [[0., 0., 0.], [0., 0., 0.]],
+])
+def test_mass_points_include_zero_coordinates(backend, strict_crossings, points):
+    t = backend
+    valid = t.zeros((2, 12), dtype=bool)
+    valid[1, 0] = valid[1, 5] = True
+    data = DualContouringData(
+        xyz_on_edge=t.array(points, dtype=BackendTensor.dtype_obj), valid_edges=valid,
+        xyz_on_centers=t.zeros((2, 3)), dxdydz=(1., 1., 1.),
+        n_surfaces_to_export=1, left_right_codes=None,
+        gradients=t.array([[0., 1., 0.], [0., 1., 0.]], dtype=BackendTensor.dtype_obj),
+        strict_crossings=strict_crossings,
+    )
+    vertices = generate_dual_contouring_vertices(data, debug=True)
+    expected = np.mean(points, axis=0)
+    np.testing.assert_array_equal(t.to_numpy(data.bias_center_mass), np.tile(expected, (3, 1)))
+    # Preserve the existing non-strict PyTorch origin regularization.
+    if BackendTensor.engine_backend == AvailableBackends.PYTORCH and not strict_crossings:
+        expected = expected / 1.0001
+    np.testing.assert_allclose(t.to_numpy(vertices), expected[None, :], atol=1e-15)