[ENH] Add dual-contouring extent capping implementation with tests and documentation

docs/extent_capping.md

@@ -0,0 +1,114 @@
+# 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.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. 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, 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
+
+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. 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 existing corner samples.
+Large-depth GPU and production-scale memory benchmarks remain outstanding.

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

gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py

@@ -24,7 +24,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 .extent_capping import cap_meshes_at_extent
 
 
 @gempy_profiler_decorator
@@ -52,6 +53,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 +88,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 +113,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 +129,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 +159,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,6 +220,20 @@ 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
+
+    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
 
 

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

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):

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)"

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
 

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
-

gempy_engine/core/data/regular_grid.py

@@ -28,7 +28,7 @@ def __len__(self):
 
     def __post_init__(self):
         self.regular_grid_shape = BackendTensor.t.array(self.regular_grid_shape)
-        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()
 
@@ -48,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",