[ENH] Add optional octree corner evaluation deduplication

docs/octree_refinement.md

@@ -60,7 +60,7 @@ change even in fast mode.
 For an interior planar one-cell selection, the full halo approaches 3x the parent
 count; isolated selections can reach 27x (7x in balanced mode). Each selected
 parent contributes eight centers and 64 stored corner rows at the next level.
-No corner deduplication or interpolation caching is introduced here.
+Stored corner rows retain that layout even when evaluation deduplication is enabled.
 
 A NumPy float64 lookup-only smoke benchmark on a dense `32 x 32 x 32` lattice
 gave the following counts (not an end-to-end interpolation benchmark):
@@ -77,3 +77,76 @@ thresholds or measurements of interpolation/reporting overhead.
 
 The default remains fast. Representative curved, multi-stack, faulted, and GPU
 time/memory benchmarks are still needed before recommending a different default.
+
+## Opt-in Evaluation and Triangulation
+
+Corner deduplication is independent of the refinement mode and defaults to `False`:
+
+```python
+options.evaluation_options.deduplicate_octree_corners = True
+```
+
+Set this selector back to `False` to evaluate the full corner layout. It is
+included in `InterpolationOptions` JSON serialization.
+
+Corner deduplication uses signed integer lattice coordinates and vectorized
+unique/inverse operations on the active NumPy or Torch device. Each unique corner
+is evaluated at its first existing physical row, not reconstructed from the extent
+origin (which can shift across refinement levels). Scalar and gradient fields are
+gathered back to the full original layout before surface-point metadata, fault
+processing, segmentation, refinement, or mesh extraction. Centers, dense/custom
+grids, sections, topography, geophysics points, and appended surface points are
+never merged with corners or with each other. No lookup or evaluated field is
+cached across calls, stacks, or levels.
+
+Fault evaluation columns are gathered with the same indices only when duplicate
+corners have identical fault values. Otherwise that evaluation uses the legacy
+path. Differentiable corner coordinates or differentiable fault-value rows also
+use the legacy path: merging independent row derivatives would change autograd.
+Gradients with respect to weights, model inputs, and appended surface points are
+preserved. Empty, non-corner, and incompatible/custom corner layouts fall back
+safely. Physical duplicates can differ by roundoff; checks allow 32 dtype epsilons
+of relative/absolute error, so output parity is numerical rather than bitwise.
+Torch requires `scatter_reduce_` support. Small grids may not benefit from the
+unique operation, gathers, and equality checks (which can synchronize a GPU).
+
+Normal and flat stacks support the selector. Fused PyKeOps evaluation compresses
+each eligible stack independently, performs one block-sparse reduction with the
+different reduced lengths, and restores each result before attaching metadata.
+Ineligible stacks retain their full rows within the same fused call. Backend
+selection and existing finite-fault dispatch restrictions remain unchanged.
+External interpolation callbacks keep their existing path and full grid layout.
+
+### Unique-Edge Quads
+
+To select quad-based connectivity instead of the legacy triangle construction:
+
+```python
+from gempy_engine.core.data.options.evaluation_options import TriangulationMethod
+
+options.evaluation_options.triangulation_method = TriangulationMethod.QUADS
+```
+
+The default is `TriangulationMethod.LEGACY`. This selector is serialized with the
+other evaluation options and is independent of corner deduplication and refinement
+mode. Legacy triangulation sorts voxel codes locally for each edge case; quad mode
+uses one sorted cell lookup.
+
+Quad mode identifies each primal edge by its lower integer endpoint and direction,
+deduplicates these identities, and finds the four incident cells. Each complete
+crossing edge produces one quad, split deterministically into two triangles for
+the existing mesh output format. Winding follows the crossing-edge gradients.
+The existing tolerant crossing rule is preserved; inconsistent crossing flags
+on shared edges raise an error rather than silently creating inconsistent faces.
+
+Incomplete quads are skipped, never emitted as partial triangles.
+`mesh.dc_data.triangulation_report` records crossing edges, complete quads, and
+missing support at physical boundaries, geological masks, and internal refinement
+boundaries. Missing-cell counts are incidences and boundary categories can overlap.
+Direct callers without pre-mask cell coordinates receive an unknown interior
+boundary classification. Counts precede subsequent overlap/fault triangle removal.
+
+This is same-level connectivity, not coarse/fine transition stitching or extent
+capping, and it does not resolve ambiguous topology or guarantee watertightness.
+NumPy and CPU Torch parity are tested; GPU behavior and end-to-end performance
+still require validation.

gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py

@@ -152,7 +152,9 @@ def dual_contouring_multi_scalar(
                 gradients=output_on_edges[n_scalar_field][slice_object],
                 n_surfaces_to_export=n_scalar_field,
                 tree_depth=options.number_octree_levels,
-                base_number=base_number
+                base_number=base_number,
+                triangulation_method=options.evaluation_options.triangulation_method,
+                generated_cell_coordinates=left_right_codes
             )
 
             dc_data_per_surface_all.append(dc_data_per_surface)

gempy_engine/API/interp_single/_interp_scalar_field.py

@@ -1,10 +1,12 @@
 from typing import Any, Union
+from copy import copy
 
 import numpy as np
 from numpy import dtype, ndarray
 
 import gempy_engine.config
 from ...core.backend_tensor import BackendTensor
+from ...core.data.engine_grid import EngineGrid
 from ...core.data.exported_fields import ExportedFields
 from ...core.data.internal_structs import SolverInput, SolverInput_v2, EvaluatorInput
 from ...core.data.options import KernelOptions, InterpolationOptions
@@ -99,10 +101,83 @@ def _solve_interpolation_result(
     return result
 
 
-def _evaluate_sys_eq(eval_input: Union[SolverInput, EvaluatorInput], weights: np.ndarray, options: InterpolationOptions) -> ExportedFields:
+def _evaluate_sys_eq(eval_input: Union[SolverInput, EvaluatorInput], weights: np.ndarray, options: InterpolationOptions,
+                     grid: EngineGrid | None = None) -> ExportedFields:
+    inverse = None
+    if options.evaluation_options.deduplicate_octree_corners:
+        eval_input, inverse = _deduplicate_corners(eval_input, grid)
     if BackendTensor.use_pykeops:
         exported_fields = symbolic_evaluator(eval_input, weights, options)
     else:
         exported_fields = generic_evaluator(eval_input, weights, options)
 
+    if inverse is not None:
+        _restore_corner_fields(exported_fields, inverse)
+
     return exported_fields
+
+
+def _restore_corner_fields(exported_fields: ExportedFields, inverse) -> None:
+    """Expand a reduced evaluation before attaching original grid metadata."""
+    for name in ('_scalar_field', '_gx_field', '_gy_field', '_gz_field'):
+        values = getattr(exported_fields, name)
+        if values is not None:
+            index = BackendTensor.t.to_numpy(inverse) if isinstance(values, np.ndarray) else inverse
+            setattr(exported_fields, name, values[index])
+
+
+def _deduplicate_corners(eval_input: SolverInput | EvaluatorInput, grid: EngineGrid | None):
+    """Call-local evaluation view; never change grid layout or surface-point metadata."""
+    if grid is None or grid.octree_grid is None or grid.corners_grid is None:
+        return eval_input, None
+    corners = grid.corners_grid.values
+    # Separate coordinate/fault row derivatives must not be redirected to a representative.
+    if len(corners) == 0 or getattr(corners, 'requires_grad', False):
+        return eval_input, None
+    faults = eval_input.fault_internal
+    fault_values = faults.fault_values_everywhere if faults.n_faults else None
+    if getattr(fault_values, 'requires_grad', False):
+        return eval_input, None
+
+    t = BackendTensor.t
+    offsets = t.array([[x, y, z] for x in (0, 1) for y in (0, 1) for z in (0, 1)], dtype='int64')
+    coordinates = (grid.octree_grid.integer_coordinates[:, None, :] + offsets).reshape(-1, 3)
+    if len(coordinates) != len(corners):
+        return eval_input, None
+    if BackendTensor.engine_backend == gempy_engine.config.AvailableBackends.PYTORCH:
+        import torch
+        unique, inverse = torch.unique(coordinates, dim=0, return_inverse=True)
+        first = torch.full((len(unique),), len(corners), dtype=torch.int64, device=coordinates.device)
+        first.scatter_reduce_(0, inverse, torch.arange(len(corners), device=coordinates.device), reduce='amin')
+    else:
+        _, first, inverse = np.unique(coordinates, axis=0, return_index=True, return_inverse=True)
+    if len(first) == len(corners):
+        return eval_input, None
+
+    start, stop = grid.corners_grid_slice.start, grid.corners_grid_slice.stop
+    xyz = eval_input.xyz_to_interpolate
+    # Use original physical rows, not extent + lattice * spacing: refined extents
+    # may carry a different origin shift. Reject non-lattice/custom corner layouts.
+    tolerance = 32 * np.finfo(BackendTensor.dtype).eps
+    if not t.allclose(xyz[start:stop], xyz[start + first][inverse], rtol=tolerance, atol=tolerance):
+        return eval_input, None
+    if fault_values is not None:
+        if not t.all(fault_values[:, start:stop] == fault_values[:, start + first][:, inverse]):
+            return eval_input, None
+
+    before = BackendTensor.arange(start, dtype='int64')
+    after = stop + BackendTensor.arange(len(xyz) - stop, dtype='int64')
+    keep = t.concatenate((before, start + first, after))
+    restore = t.concatenate((before, start + inverse,
+                             start + len(first) + BackendTensor.arange(len(xyz) - stop, dtype='int64')))
+    reduced = copy(eval_input)
+    reduced.xyz_to_interpolate = xyz[keep]
+    if fault_values is not None:
+        reduced_faults = copy(faults)
+        reduced_faults.fault_values_everywhere = fault_values[:, keep]
+        if isinstance(reduced, EvaluatorInput):
+            reduced.solver_input = copy(reduced.solver_input)
+            reduced.solver_input.fault_internal = reduced_faults
+        else:
+            reduced.fault_internal = reduced_faults
+    return reduced, restore

gempy_engine/API/interp_single/_interp_single_feature.py

@@ -35,7 +35,7 @@ def interpolate_feature_with_cokrig(interpolation_input: InterpolationInput,
     xyz = solver_input.xyz_to_interpolate
 
     weights = compute_weights(solver_input, stack_number, options)
-    exported_fields: ExportedFields = _evaluate_sys_eq(solver_input, weights, options)
+    exported_fields: ExportedFields = _evaluate_sys_eq(solver_input, weights, options, grid=grid)
 
     exported_fields.set_structure_values(
         reference_sp_position=data_shape.reference_sp_position,

gempy_engine/API/interp_single/_stack_ops.py

@@ -1,7 +1,7 @@
 import concurrent.futures
 
 from ._aux_faults_ops import _grab_stack_fault_data, _modify_faults_values_output, _options_with_finite_fault_gradients
-from ._interp_scalar_field import _evaluate_sys_eq, compute_weights
+from ._interp_scalar_field import _deduplicate_corners, _evaluate_sys_eq, _restore_corner_fields, compute_weights
 from ._interp_single_feature import interpolate_feature_with_external_function
 from ...config import AvailableBackends
 from ...core.backend_tensor import BackendTensor
@@ -225,6 +225,7 @@ def _evaluate(interpolation_inputs: list[InterpolationInput], options: Interpola
             eval_input=eval_input,
             weights=eval_input.solver_input.weights_x0,
             options=options_per_stack[idx] if options_per_stack is not None else options,
+            grid=interpolation_inputs[idx].grid,
         )
 
         exported_fields.set_structure_values_from_eval_input(eval_input)
@@ -257,15 +258,24 @@ def _evaluate_optimized(interpolation_inputs: list[InterpolationInput], options:
     weights_list = [ei.solver_input.weights_x0 for ei in eval_inputs]
 
     options_list = options_per_stack or [options] * len(stack_indices)
+    reduced_inputs = []
+    inverses = []
+    for eval_input, interpolation_input, options_i in zip(eval_inputs, interpolation_inputs, options_list):
+        reduced, inverse = (_deduplicate_corners(eval_input, interpolation_input.grid)
+                            if options_i.evaluation_options.deduplicate_octree_corners else (eval_input, None))
+        reduced_inputs.append(reduced)
+        inverses.append(inverse)
 
     # Call the stacked evaluator (single PyKeOps call with block-sparse ranges)
     exported_fields_list: list[ExportedFields] = symbolic_evaluator_optimized_stacked(
-        eval_inputs=eval_inputs,
+        eval_inputs=reduced_inputs,
         weights_list=weights_list,
         options_list=options_list
     )
 
     for idx, exported_fields in enumerate(exported_fields_list):
+        if inverses[idx] is not None:
+            _restore_corner_fields(exported_fields, inverses[idx])
         exported_fields.set_structure_values_from_eval_input(eval_inputs[idx])
         exported_fields.debug = eval_inputs[idx].solver_input.debug
 

gempy_engine/core/data/dual_contouring_data.py

@@ -3,6 +3,8 @@
 
 import numpy as np
 
+from .options.evaluation_options import TriangulationMethod
+
 
 @dataclass(init=True)
 class DualContouringData:
@@ -27,6 +29,9 @@ class DualContouringData:
     extra_edge_xyz: Optional[np.ndarray] = None      # (n_valid_voxels, K, 3)
     extra_edge_normals: Optional[np.ndarray] = None   # (n_valid_voxels, K, 3)
     extra_weights: Optional[np.ndarray] = None        # (n_valid_voxels, K)
+    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.
 
     @property
     def valid_voxels(self):
@@ -39,4 +44,3 @@ def n_valid_edges(self):
     @property
     def n_evaluations_on_edges(self):
         return self.xyz_on_edge.shape[0]
-

gempy_engine/core/data/options/evaluation_options.py

@@ -22,6 +22,11 @@ class MeshExtractionMaskingOptions(enum.Enum):
     RAW = enum.auto()
 
 
+class TriangulationMethod(str, enum.Enum):
+    LEGACY = "legacy"
+    QUADS = "quads"
+
+
 @dataclass
 class EvaluationOptions:
     _number_octree_levels: int = 1
@@ -30,6 +35,8 @@ class EvaluationOptions:
     octree_error_threshold: float = 1.  #: Number of standard deviations to consider a voxel as candidate to refine
     octree_min_level: int = 2
     octree_refinement_mode: OctreeRefinementMode = OctreeRefinementMode.FAST
+    deduplicate_octree_corners: bool = False  #: Evaluate unique corners, then restore the original row layout.
+    triangulation_method: TriangulationMethod = TriangulationMethod.LEGACY
 
     mesh_extraction: bool = True
     mesh_extraction_masking_options: MeshExtractionMaskingOptions = MeshExtractionMaskingOptions.INTERSECT

gempy_engine/modules/dual_contouring/dual_contouring_v2.py

@@ -7,8 +7,10 @@
 from ...core.backend_tensor import BackendTensor
 from ...core.data.dual_contouring_data import DualContouringData
 from ...core.data.dual_contouring_mesh import DualContouringMesh
+from ...core.data.options.evaluation_options import TriangulationMethod
 from ...core.utils import gempy_profiler_decorator
 from ...modules.dual_contouring.fancy_triangulation import triangulate
+from ...modules.dual_contouring.quad_triangulation import triangulate_quads
 
 
 @gempy_profiler_decorator
@@ -73,6 +75,15 @@ def _compute_triangulation(dc_data_per_surface: DualContouringData,
                            left_right_codes, edges_normals, vertex):
     """Compute triangulation indices for a specific surface."""
 
+    method = TriangulationMethod(dc_data_per_surface.triangulation_method)
+    if method == TriangulationMethod.QUADS:
+        return triangulate_quads(
+            left_right_codes, dc_data_per_surface.valid_edges, edges_normals, vertex,
+            dc_data_per_surface.base_number,
+            generated_coordinates=dc_data_per_surface.generated_cell_coordinates,
+            report=dc_data_per_surface.triangulation_report
+        )
+
     # * Fancy triangulation 👗
     valid_voxels = dc_data_per_surface.valid_voxels