From c55fa4a1a7cc4d753f2402ca3e02021e64385965 Mon Sep 17 00:00:00 2001
From: Miguel de la Varga <miguel@terranigma-solutions.com>
Date: Mon, 21 Sep 2026 14:13:25 +0200
Subject: [PATCH 1/9] [ENH] Add optional octree corner evaluation deduplication

---
 docs/octree_refinement.md                     |  39 ++-
 .../API/interp_single/_interp_scalar_field.py |  71 +++++-
 .../interp_single/_interp_single_feature.py   |   2 +-
 gempy_engine/API/interp_single/_stack_ops.py  |   6 +
 .../core/data/options/evaluation_options.py   |   1 +
 .../modules/evaluator/generic_evaluator.py    |   2 +-
 .../test_modules/test_octree_optimizations.py | 235 ++++++++++++++++++
 7 files changed, 352 insertions(+), 4 deletions(-)
 create mode 100644 tests/test_common/test_modules/test_octree_optimizations.py

diff --git a/docs/octree_refinement.md b/docs/octree_refinement.md
index 4c965159..b4d65163 100644
--- a/docs/octree_refinement.md
+++ b/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,40 @@ 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 Corner Evaluation
+
+Corner deduplication is independent of the refinement mode and defaults to `False`:
+
+```python
+options.evaluation_options.deduplicate_octree_corners = True
+```
+
+Set the selector back to `False` to use its legacy path. The selector 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. With deduplication enabled for any
+stack in a flat chunk, that chunk uses per-stack evaluation instead of the fused
+PyKeOps evaluator; each stack still selects its usual dense or symbolic backend.
+External interpolation callbacks keep their existing path and full grid layout.
diff --git a/gempy_engine/API/interp_single/_interp_scalar_field.py b/gempy_engine/API/interp_single/_interp_scalar_field.py
index fc6d8961..1e350a5f 100644
--- a/gempy_engine/API/interp_single/_interp_scalar_field.py
+++ b/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,77 @@ 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:
+        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])
     return exported_fields
+
+
+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
diff --git a/gempy_engine/API/interp_single/_interp_single_feature.py b/gempy_engine/API/interp_single/_interp_single_feature.py
index 294cd6cb..ce7d32c2 100644
--- a/gempy_engine/API/interp_single/_interp_single_feature.py
+++ b/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,
diff --git a/gempy_engine/API/interp_single/_stack_ops.py b/gempy_engine/API/interp_single/_stack_ops.py
index d2b3969b..5b6a6c6c 100644
--- a/gempy_engine/API/interp_single/_stack_ops.py
+++ b/gempy_engine/API/interp_single/_stack_ops.py
@@ -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)
@@ -241,6 +242,11 @@ def _evaluate_optimized(interpolation_inputs: list[InterpolationInput], options:
                          options_per_stack: list[InterpolationOptions] | None = None) -> tuple[list[EvaluatorInput], list[ExportedFields]]:
     from gempy_engine.modules.evaluator.symbolic_evaluator import symbolic_evaluator_optimized_stacked
 
+    # Use the per-stack evaluator so expansion precedes metadata and segmentation.
+    if any(o.evaluation_options.deduplicate_octree_corners for o in (options_per_stack or [options])):
+        return _evaluate(interpolation_inputs, options, solver_inputs, stack_structure,
+                         tensor_structs, stack_indices, options_per_stack)
+
     eval_inputs: list[EvaluatorInput] = []
     for idx, global_i in enumerate(stack_indices):
         stack_structure.stack_number = global_i
diff --git a/gempy_engine/core/data/options/evaluation_options.py b/gempy_engine/core/data/options/evaluation_options.py
index b302004c..a41406eb 100644
--- a/gempy_engine/core/data/options/evaluation_options.py
+++ b/gempy_engine/core/data/options/evaluation_options.py
@@ -30,6 +30,7 @@ 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.
     
     mesh_extraction: bool = True
     mesh_extraction_masking_options: MeshExtractionMaskingOptions = MeshExtractionMaskingOptions.INTERSECT
diff --git a/gempy_engine/modules/evaluator/generic_evaluator.py b/gempy_engine/modules/evaluator/generic_evaluator.py
index c945a24f..1fb8402d 100644
--- a/gempy_engine/modules/evaluator/generic_evaluator.py
+++ b/gempy_engine/modules/evaluator/generic_evaluator.py
@@ -61,7 +61,7 @@ def generic_evaluator(
                 gz_field[slice_array] = gz_chunk  # type: ignore
 
     # Force garbage collection every few chunks to prevent memory buildup
-    if (i + 1) % 5 == 0 or i == n_chunks - 1:
+    if n_chunks:
         gc.collect()
         
     if n_chunks > 5:
diff --git a/tests/test_common/test_modules/test_octree_optimizations.py b/tests/test_common/test_modules/test_octree_optimizations.py
new file mode 100644
index 00000000..ef53db71
--- /dev/null
+++ b/tests/test_common/test_modules/test_octree_optimizations.py
@@ -0,0 +1,235 @@
+from unittest.mock import patch
+
+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 import InterpolationOptions
+from gempy_engine.core.data.engine_grid import EngineGrid
+from gempy_engine.core.data.generic_grid import GenericGrid
+from gempy_engine.core.data.internal_structs import SolverInput, EvaluatorInput
+from gempy_engine.core.data.kernel_classes.faults import FaultsData
+from gempy_engine.core.data.regular_grid import RegularGrid
+from gempy_engine.API.interp_single._octree_generation import _generate_corners
+from gempy_engine.API.interp_single._interp_scalar_field import _deduplicate_corners, _evaluate_sys_eq, _solve_interpolation
+from gempy_engine.API.interp_single._interp_single_feature import input_preprocess
+from gempy_engine.modules.octrees_topology._octree_common import _generate_next_level_centers
+from tests.fixtures.simple_models import simple_model_interpolation_input_factory
+
+
+@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.COMPUTE_GRADS)
+    BackendTensor._change_backend(engine_backend=AvailableBackends[request.param], use_gpu=False,
+                                 dtype='float64', use_pykeops=False, grads=True)
+    yield request.param
+    BackendTensor._change_backend(engine_backend=old[0], use_gpu=old[1], dtype=old[2], use_pykeops=old[3], grads=old[4])
+
+
+def corner_grid(sparse=False):
+    t = BackendTensor.t
+    root = RegularGrid([10, 13, -4, -2, 20, 24], [3, 2, 2])
+    if sparse:
+        active = t.array([True, True] + [False] * 10, dtype=bool)
+        xyz, bits = _generate_next_level_centers(root.values[active], root.dxdydz)
+        root = RegularGrid.from_octree_level(xyz, root, active, bits)
+    corners = _generate_corners(root)
+    return EngineGrid(octree_grid=root, corners_grid=GenericGrid(corners),
+                      custom_grid=GenericGrid(corners[:2]))
+
+
+@pytest.mark.parametrize('sparse', [False, True])
+@pytest.mark.parametrize('faulted', [False, True])
+def test_corner_mapping(backend, sparse, faulted):
+    t = BackendTensor.t
+    grid = corner_grid(sparse)
+    xyz = t.concatenate((grid.values, grid.corners_grid.values[:3]))
+    faults = FaultsData(xyz[:, 0][None, :], xyz[-3:, 0][None, :]) if faulted else None
+    original = SolverInput(None, None, xyz, faults)
+    reduced, inverse = _deduplicate_corners(original, grid)
+    assert inverse is not None
+    assert len(reduced.xyz_to_interpolate) < len(xyz)
+    np.testing.assert_allclose(t.to_numpy(reduced.xyz_to_interpolate[inverse]), t.to_numpy(xyz), atol=1e-13, rtol=0)
+    start = grid.corners_grid_slice.start
+    np.testing.assert_array_equal(t.to_numpy(reduced.xyz_to_interpolate[:start]), t.to_numpy(xyz[:start]))
+    np.testing.assert_array_equal(t.to_numpy(reduced.xyz_to_interpolate[-3:]), t.to_numpy(xyz[-3:]))
+    assert original.xyz_to_interpolate is xyz
+    if faulted:
+        np.testing.assert_array_equal(t.to_numpy(reduced.fault_internal.fault_values_everywhere[:, inverse]),
+                                      t.to_numpy(faults.fault_values_everywhere))
+        assert original.fault_internal is faults
+        # Different fault values at shared corners must not be merged.
+        faults.fault_values_everywhere = BackendTensor.arange(len(xyz), dtype='float64')[None, :]
+        assert _deduplicate_corners(original, grid)[1] is None
+
+
+def test_corner_fallbacks(backend):
+    t = BackendTensor.t
+    grid = corner_grid()
+    original = SolverInput(None, None, grid.values)
+    assert _deduplicate_corners(original, None)[1] is None
+    grid.corners_grid.values = grid.corners_grid.values[:0]
+    assert _deduplicate_corners(original, grid)[1] is None
+    grid = corner_grid()
+    grid.corners_grid.values[0] += 0.1
+    original.xyz_to_interpolate = grid.values
+    # The first corner is unique; alter a corner that is shared instead.
+    grid.corners_grid.values[7] += 0.1
+    original.xyz_to_interpolate = grid.values
+    assert _deduplicate_corners(original, grid)[1] is None
+    if backend == 'PYTORCH':
+        grid = corner_grid()
+        grid.corners_grid.values.requires_grad_()
+        assert _deduplicate_corners(original, grid)[1] is None
+        grid.corners_grid.values = grid.corners_grid.values.detach()
+        original.fault_internal = FaultsData(t.ones((1, len(original.xyz_to_interpolate))).requires_grad_(), t.ones((1, 3)))
+        assert _deduplicate_corners(original, grid)[1] is None
+
+
+@pytest.mark.parametrize('flat_input', [False, True])
+@pytest.mark.parametrize('symbolic', [False, True])
+def test_evaluation_and_gradient_parity(backend, flat_input, symbolic, monkeypatch):
+    if symbolic:
+        pytest.importorskip('pykeops')
+        monkeypatch.setattr(BackendTensor, 'use_pykeops', True)
+    t = BackendTensor.t
+    interp, options, descriptor = simple_model_interpolation_input_factory()
+    grid = interp.grid
+    grid.corners_grid = GenericGrid(_generate_corners(grid.octree_grid))
+    options.evaluation_options.compute_scalar_gradient = True
+    options.evaluation_options.evaluation_chunk_size = 300
+    if backend == 'PYTORCH':
+        interp.surface_points.sp_coords.requires_grad_()
+    solver = input_preprocess(descriptor.tensors_structure, interp)
+    weights = _solve_interpolation(solver, options.kernel_options)
+    eval_input = EvaluatorInput(solver, interp, descriptor.tensors_structure) if flat_input else solver
+    original_xyz = eval_input.xyz_to_interpolate
+    legacy = _evaluate_sys_eq(eval_input, weights, options, grid)
+    options.evaluation_options.deduplicate_octree_corners = True
+    reduced, inverse = _deduplicate_corners(eval_input, grid)
+    assert inverse is not None
+    optimized = _evaluate_sys_eq(eval_input, weights, options, grid)
+    assert eval_input.xyz_to_interpolate is original_xyz
+    for name in ('_scalar_field', '_gx_field', '_gy_field', '_gz_field'):
+        a, b = getattr(legacy, name), getattr(optimized, name)
+        np.testing.assert_allclose(t.to_numpy(a), t.to_numpy(b), atol=1e-10, rtol=1e-10)
+    for result in (legacy, optimized):
+        result.set_structure_values(descriptor.tensors_structure.reference_sp_position, interp.slice_feature, grid.len_all_grids)
+    np.testing.assert_allclose(t.to_numpy(legacy.scalar_field_at_surface_points), t.to_numpy(optimized.scalar_field_at_surface_points))
+    if backend == 'PYTORCH':
+        import torch
+        coefficients = torch.arange(len(original_xyz), dtype=weights.dtype) + 1
+        def loss(result):
+            return sum((getattr(result, name) * coefficients).sum()
+                       for name in ('_scalar_field', '_gx_field', '_gy_field', '_gz_field'))
+        targets = (weights, interp.surface_points.sp_coords)
+        a = torch.autograd.grad(loss(legacy), targets, retain_graph=True)
+        b = torch.autograd.grad(loss(optimized), targets)
+        for old, new in zip(a, b):
+            torch.testing.assert_close(old, new, atol=1e-7, rtol=1e-8)
+
+
+@pytest.mark.parametrize('flat', [False, True])
+def test_model_parity(backend, flat, monkeypatch):
+    from gempy_engine.API.model.model_api import compute_model
+    from gempy_engine.API.interp_single import _multi_scalar_field_manager as manager
+
+    monkeypatch.setenv('GEMPY_FLAT_STACKS', 'False')
+    interp, options, descriptor = simple_model_interpolation_input_factory()
+    options.evaluation_options.number_octree_levels = 2
+    legacy = compute_model(interp, options, descriptor)
+    interp, options, descriptor = simple_model_interpolation_input_factory()
+    options.evaluation_options.number_octree_levels = 2
+    options.evaluation_options.deduplicate_octree_corners = True
+    monkeypatch.setenv('GEMPY_FLAT_STACKS', str(flat))
+    if flat:
+        # Public flat dispatch requires PyKeOps; exercise the same stack manager
+        # with dense per-stack evaluation here without requiring the JIT compiler.
+        monkeypatch.setattr(manager, '_interpolate_stack', manager._interpolate_stack_flat)
+    optimized = compute_model(interp, options, descriptor)
+    t = BackendTensor.t
+    for old_level, new_level in zip(legacy.octrees_output, optimized.octrees_output):
+        np.testing.assert_allclose(t.to_numpy(old_level.grid.values), t.to_numpy(new_level.grid.values))
+        for old, new in zip(old_level.outputs, new_level.outputs):
+            np.testing.assert_allclose(t.to_numpy(old.exported_fields.scalar_field_everywhere),
+                                       t.to_numpy(new.exported_fields.scalar_field_everywhere), atol=1e-9)
+    assert len(legacy.dc_meshes) == len(optimized.dc_meshes)
+    for old, new in zip(legacy.dc_meshes, optimized.dc_meshes):
+        np.testing.assert_allclose(old.vertices, new.vertices, atol=1e-8)
+        np.testing.assert_array_equal(old.edges, new.edges)
+
+
+@pytest.mark.parametrize('flat', [False, True])
+def test_fault_model_parity(backend, flat, monkeypatch):
+    from gempy_engine.API.interp_single import _multi_scalar_field_manager as manager
+    from gempy_engine.API.model.model_api import compute_model
+    from tests.fixtures.complex_geometries import graben_fault_model
+
+    monkeypatch.setenv('GEMPY_FLAT_STACKS', 'False')
+    results = []
+    for optimized in (False, True):
+        interp, descriptor, options = graben_fault_model.__wrapped__()
+        options.evaluation_options.number_octree_levels = 2
+        options.evaluation_options.mesh_extraction = False
+        options.evaluation_options.deduplicate_octree_corners = optimized
+        if flat and optimized:
+            monkeypatch.setattr(manager, '_interpolate_stack', manager._interpolate_stack_flat)
+        results.append(compute_model(interp, options, descriptor))
+    t = BackendTensor.t
+    for old_level, new_level in zip(results[0].octrees_output, results[1].octrees_output):
+        for old, new in zip(old_level.outputs, new_level.outputs):
+            np.testing.assert_allclose(t.to_numpy(old.exported_fields.scalar_field_everywhere),
+                                       t.to_numpy(new.exported_fields.scalar_field_everywhere), atol=1e-8)
+            np.testing.assert_allclose(t.to_numpy(old.scalar_fields.values_block),
+                                       t.to_numpy(new.scalar_fields.values_block), atol=1e-8)
+
+
+def test_empty_evaluation_and_isolated_cell(backend):
+    t = BackendTensor.t
+    options = InterpolationOptions.from_args(range=1., c_o=1.)
+    options.evaluation_options.compute_scalar_gradient = True
+    options.evaluation_options.deduplicate_octree_corners = True
+    empty = SolverInput(None, None, t.zeros((0, 3)))
+    result = _evaluate_sys_eq(empty, t.ones(1), options)
+    assert result.scalar_field_everywhere.shape == (0,)
+    assert result.gx_field_everywhere.shape == (0,)
+    grid = corner_grid()
+    grid.octree_grid._integer_coordinates = grid.octree_grid.integer_coordinates[:1]
+    grid.octree_grid.values = grid.octree_grid.values[:1]
+    grid.corners_grid.values = grid.corners_grid.values[:8]
+    original = SolverInput(None, None, grid.values)
+    assert _deduplicate_corners(original, grid)[1] is None
+
+
+def test_public_pykeops_flat_parity(backend, monkeypatch):
+    pytest.importorskip('pykeops')
+    from gempy_engine.API.model.model_api import compute_model
+    from gempy_engine.API.interp_single import _stack_ops
+
+    monkeypatch.setattr(BackendTensor, 'use_pykeops', True)
+    monkeypatch.setenv('GEMPY_FLAT_STACKS', 'True')
+    results = []
+    with patch.object(_stack_ops, '_evaluate_optimized', wraps=_stack_ops._evaluate_optimized) as evaluate:
+        for optimized in (False, True):
+            interp, options, descriptor = simple_model_interpolation_input_factory()
+            options.evaluation_options.number_octree_levels = 2
+            options.evaluation_options.mesh_extraction = False
+            options.evaluation_options.deduplicate_octree_corners = optimized
+            results.append(compute_model(interp, options, descriptor))
+        assert evaluate.call_count == 4
+    t = BackendTensor.t
+    for old, new in zip(results[0].octrees_output, results[1].octrees_output):
+        np.testing.assert_allclose(t.to_numpy(old.outputs[0].exported_fields.scalar_field_everywhere),
+                                   t.to_numpy(new.outputs[0].exported_fields.scalar_field_everywhere), atol=1e-9)
+
+
+def test_selectors_serialization():
+    options = InterpolationOptions.from_args(range=1., c_o=1.)
+    assert not options.evaluation_options.deduplicate_octree_corners
+    options.evaluation_options.deduplicate_octree_corners = True
+    restored = InterpolationOptions.model_validate_json(options.model_dump_json())
+    assert restored.evaluation_options.deduplicate_octree_corners

From ac8d5eb0c05adb09c9adfaf73de9a7fad29dc09d Mon Sep 17 00:00:00 2001
From: Miguel de la Varga <miguel@terranigma-solutions.com>
Date: Mon, 21 Sep 2026 14:14:53 +0200
Subject: [PATCH 2/9] [ENH] Add optional sort-once triangulation lookup

---
 docs/octree_refinement.md                     | 13 ++++++--
 .../multi_scalar_dual_contouring.py           |  3 +-
 .../core/data/dual_contouring_data.py         |  2 +-
 .../core/data/options/evaluation_options.py   |  1 +
 .../dual_contouring/dual_contouring_v2.py     |  3 +-
 .../dual_contouring/fancy_triangulation.py    | 22 ++++++++-----
 .../test_modules/test_octree_optimizations.py | 31 +++++++++++++++++++
 7 files changed, 62 insertions(+), 13 deletions(-)

diff --git a/docs/octree_refinement.md b/docs/octree_refinement.md
index b4d65163..1bbb3b68 100644
--- a/docs/octree_refinement.md
+++ b/docs/octree_refinement.md
@@ -78,15 +78,16 @@ 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 Corner Evaluation
+## Opt-in Evaluation and Triangulation
 
-Corner deduplication is independent of the refinement mode and defaults to `False`:
+Both optimizations are independent of the refinement mode and default to `False`:
 
 ```python
 options.evaluation_options.deduplicate_octree_corners = True
+options.evaluation_options.triangulation_sort_once = True
 ```
 
-Set the selector back to `False` to use its legacy path. The selector is
+Set either selector back to `False` to use its legacy path. Both selectors are
 included in `InterpolationOptions` JSON serialization.
 
 Corner deduplication uses signed integer lattice coordinates and vectorized
@@ -114,3 +115,9 @@ Normal and flat stacks support the selector. With deduplication enabled for any
 stack in a flat chunk, that chunk uses per-stack evaluation instead of the fused
 PyKeOps evaluator; each stack still selects its usual dense or symbolic backend.
 External interpolation callbacks keep their existing path and full grid layout.
+
+The triangulation selector sorts the active voxel codes once per surface call and
+reuses that lookup across all six edge cases. It preserves triangle order,
+neighbor filtering, and normal correction; it does not rewrite edge topology or
+change vertex generation. Its lookup is local to the surface call and safe for
+parallel surface processing.
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 654a0cf5..ccd1a9a2 100644
--- a/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py
+++ b/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py
@@ -152,7 +152,8 @@ 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_sort_once=options.evaluation_options.triangulation_sort_once
             )
 
             dc_data_per_surface_all.append(dc_data_per_surface)
diff --git a/gempy_engine/core/data/dual_contouring_data.py b/gempy_engine/core/data/dual_contouring_data.py
index c879cf75..672b57c9 100644
--- a/gempy_engine/core/data/dual_contouring_data.py
+++ b/gempy_engine/core/data/dual_contouring_data.py
@@ -27,6 +27,7 @@ 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_sort_once: bool = False
 
     @property
     def valid_voxels(self):
@@ -39,4 +40,3 @@ def n_valid_edges(self):
     @property
     def n_evaluations_on_edges(self):
         return self.xyz_on_edge.shape[0]
-    
\ No newline at end of file
diff --git a/gempy_engine/core/data/options/evaluation_options.py b/gempy_engine/core/data/options/evaluation_options.py
index a41406eb..f5912bc4 100644
--- a/gempy_engine/core/data/options/evaluation_options.py
+++ b/gempy_engine/core/data/options/evaluation_options.py
@@ -31,6 +31,7 @@ class EvaluationOptions:
     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_sort_once: bool = False  #: Reuse one voxel-code sort across the six edge cases per surface.
     
     mesh_extraction: bool = True
     mesh_extraction_masking_options: MeshExtractionMaskingOptions = MeshExtractionMaskingOptions.INTERSECT
diff --git a/gempy_engine/modules/dual_contouring/dual_contouring_v2.py b/gempy_engine/modules/dual_contouring/dual_contouring_v2.py
index ae0c6c7e..cc28d1cd 100644
--- a/gempy_engine/modules/dual_contouring/dual_contouring_v2.py
+++ b/gempy_engine/modules/dual_contouring/dual_contouring_v2.py
@@ -88,7 +88,8 @@ def _compute_triangulation(dc_data_per_surface: DualContouringData,
         tree_depth=tree_depth_per_surface,
         voxel_normals=voxels_normals,
         vertex=vertex,
-        base_number=dc_data_per_surface.base_number
+        base_number=dc_data_per_surface.base_number,
+        sort_once=dc_data_per_surface.triangulation_sort_once
     )
 
     # @on
diff --git a/gempy_engine/modules/dual_contouring/fancy_triangulation.py b/gempy_engine/modules/dual_contouring/fancy_triangulation.py
index 638af625..be4b7f92 100644
--- a/gempy_engine/modules/dual_contouring/fancy_triangulation.py
+++ b/gempy_engine/modules/dual_contouring/fancy_triangulation.py
@@ -20,7 +20,7 @@ def _get_pack_factors(base_x, base_y, base_z):
     return BackendTensor.tfnp.stack([by * bz, bz, BackendTensor.t.array(1)], axis=0)
 
 
-def triangulate(left_right_array, valid_edges, tree_depth: int, voxel_normals, vertex, base_number: tuple[int, int, int] | list[int]):
+def triangulate(left_right_array, valid_edges, tree_depth: int, voxel_normals, vertex, base_number: tuple[int, int, int] | list[int], sort_once: bool = False):
     # * Variables
     # Determine base_number dynamically from the data to support arbitrary grid shapes
     base_x, base_y, base_z = base_number
@@ -32,6 +32,10 @@ def triangulate(left_right_array, valid_edges, tree_depth: int, voxel_normals, v
 
     # * Consts
     voxel_code = (left_right_array * pack_factors).sum(1).reshape(-1, 1)
+    sorted_lookup = None
+    if sort_once:
+        order = BackendTensor.tfnp.argsort(voxel_code.reshape(-1))
+        sorted_lookup = (voxel_code.reshape(-1)[order], order)
     # ----------
 
     indices = []
@@ -57,7 +61,8 @@ def triangulate(left_right_array, valid_edges, tree_depth: int, voxel_normals, v
             voxel_normals=voxel_normals,
             n=n,
             base_number=base_number,
-            pack_factors=pack_factors
+            pack_factors=pack_factors,
+            sorted_lookup=sorted_lookup
         )
 
         indices.append(indices_patch)
@@ -73,7 +78,7 @@ def triangulate(left_right_array, valid_edges, tree_depth: int, voxel_normals, v
 
 def compute_triangles_for_edge(edge_vector_a, edge_vector_b, edge_vector_c,
                                left_right_array_active_edge, voxel_code, voxel_normals, n,
-                               base_number, pack_factors):
+                               base_number, pack_factors, sorted_lookup=None):
     """
     Important concepts to understand this triangulation:
     - left_right_array (n_voxels, 3-directions) contains a unique number per direction describing if it is left (even) or right (odd) and the voxel level
@@ -114,7 +119,7 @@ def compute_triangles_for_edge(edge_vector_a, edge_vector_b, edge_vector_c,
 
     x, y, z = _get_indices_via_searchsorted(
         voxel_code,
-        compressed_idx_0, compressed_idx_1, compressed_idx_2
+        compressed_idx_0, compressed_idx_1, compressed_idx_2, sorted_lookup=sorted_lookup
     )
 
     # Step 5: Calculate normals and order triangles
@@ -257,7 +262,7 @@ def _calculate_normals_and_order_triangles(x, y, z, voxel_normals, n):
     return normal
 
 
-def _get_indices_via_searchsorted(voxel_code, compressed_0, compressed_1, compressed_2):
+def _get_indices_via_searchsorted(voxel_code, compressed_0, compressed_1, compressed_2, sorted_lookup=None):
     """
     Memory-efficient replacement for broadcasting, preserving exact original indices.
     """
@@ -267,8 +272,11 @@ def _get_indices_via_searchsorted(voxel_code, compressed_0, compressed_1, compre
         return empty, empty, empty
 
     # 1. Get sorting indices to map back to original positions later
-    sort_indices = BackendTensor.tfnp.argsort(vc_1d)
-    sorted_vc = vc_1d[sort_indices]
+    if sorted_lookup is None:
+        sort_indices = BackendTensor.tfnp.argsort(vc_1d)
+        sorted_vc = vc_1d[sort_indices]
+    else:
+        sorted_vc, sort_indices = sorted_lookup
 
     # 2. Search in the correctly sorted array (O(M log N) time)
     idx_0_sorted = BackendTensor.tfnp.searchsorted(sorted_vc, compressed_0)
diff --git a/tests/test_common/test_modules/test_octree_optimizations.py b/tests/test_common/test_modules/test_octree_optimizations.py
index ef53db71..842b3e0e 100644
--- a/tests/test_common/test_modules/test_octree_optimizations.py
+++ b/tests/test_common/test_modules/test_octree_optimizations.py
@@ -1,3 +1,4 @@
+from itertools import product
 from unittest.mock import patch
 
 import numpy as np
@@ -14,6 +15,7 @@
 from gempy_engine.API.interp_single._octree_generation import _generate_corners
 from gempy_engine.API.interp_single._interp_scalar_field import _deduplicate_corners, _evaluate_sys_eq, _solve_interpolation
 from gempy_engine.API.interp_single._interp_single_feature import input_preprocess
+from gempy_engine.modules.dual_contouring.fancy_triangulation import triangulate
 from gempy_engine.modules.octrees_topology._octree_common import _generate_next_level_centers
 from tests.fixtures.simple_models import simple_model_interpolation_input_factory
 
@@ -133,6 +135,31 @@ def loss(result):
             torch.testing.assert_close(old, new, atol=1e-7, rtol=1e-8)
 
 
+@pytest.mark.parametrize('case', ['dense', 'sparse', 'empty', 'no_edges'])
+def test_triangulation_parity_and_sort_count(backend, case):
+    t = BackendTensor.t
+    coords = np.array(list(product(range(3), range(4), range(5))), dtype=np.int64)
+    rng = np.random.default_rng(17)
+    rng.shuffle(coords)
+    if case == 'sparse':
+        coords = coords[::2]
+    if case == 'empty':
+        coords = coords[:0]
+    coords = t.array(coords, dtype='int64')
+    valid = t.array(rng.random((len(coords), 12)) > 0.3, dtype=bool)
+    if case == 'no_edges':
+        valid[:] = False
+    normals = t.array(rng.normal(size=(len(coords), 12, 3)))
+    vertices = t.array(coords, dtype='float64') + 0.5
+    with patch.object(BackendTensor.tfnp, 'argsort', wraps=BackendTensor.tfnp.argsort) as sort:
+        old = triangulate(coords, valid, 1, normals, vertices, (3, 4, 5))
+        assert sort.call_count == (0 if case == 'empty' else 6)
+        sort.reset_mock()
+        new = triangulate(coords, valid, 1, normals, vertices, (3, 4, 5), sort_once=True)
+        assert sort.call_count == 1
+    np.testing.assert_array_equal(t.to_numpy(old), t.to_numpy(new))
+
+
 @pytest.mark.parametrize('flat', [False, True])
 def test_model_parity(backend, flat, monkeypatch):
     from gempy_engine.API.model.model_api import compute_model
@@ -145,6 +172,7 @@ def test_model_parity(backend, flat, monkeypatch):
     interp, options, descriptor = simple_model_interpolation_input_factory()
     options.evaluation_options.number_octree_levels = 2
     options.evaluation_options.deduplicate_octree_corners = True
+    options.evaluation_options.triangulation_sort_once = True
     monkeypatch.setenv('GEMPY_FLAT_STACKS', str(flat))
     if flat:
         # Public flat dispatch requires PyKeOps; exercise the same stack manager
@@ -230,6 +258,9 @@ def test_public_pykeops_flat_parity(backend, monkeypatch):
 def test_selectors_serialization():
     options = InterpolationOptions.from_args(range=1., c_o=1.)
     assert not options.evaluation_options.deduplicate_octree_corners
+    assert not options.evaluation_options.triangulation_sort_once
     options.evaluation_options.deduplicate_octree_corners = True
+    options.evaluation_options.triangulation_sort_once = True
     restored = InterpolationOptions.model_validate_json(options.model_dump_json())
     assert restored.evaluation_options.deduplicate_octree_corners
+    assert restored.evaluation_options.triangulation_sort_once

From 879eb7c1a1691457c0795c79db3a0d7ee00457da Mon Sep 17 00:00:00 2001
From: Miguel de la Varga <miguel@terranigma-solutions.com>
Date: Mon, 21 Sep 2026 15:05:22 +0200
Subject: [PATCH 3/9] [ENH] Preserve fused evaluation with corner deduplication

---
 docs/octree_refinement.md                     |   9 +-
 .../API/interp_single/_interp_scalar_field.py |   5 +
 gempy_engine/API/interp_single/_stack_ops.py  |  15 +-
 .../modules/evaluator/symbolic_evaluator.py   |  64 +++----
 .../test_fused_corner_evaluation.py           | 165 ++++++++++++++++++
 5 files changed, 221 insertions(+), 37 deletions(-)
 create mode 100644 tests/test_common/test_modules/test_fused_corner_evaluation.py

diff --git a/docs/octree_refinement.md b/docs/octree_refinement.md
index 1bbb3b68..0cde850e 100644
--- a/docs/octree_refinement.md
+++ b/docs/octree_refinement.md
@@ -111,9 +111,12 @@ 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. With deduplication enabled for any
-stack in a flat chunk, that chunk uses per-stack evaluation instead of the fused
-PyKeOps evaluator; each stack still selects its usual dense or symbolic backend.
+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. Without
+PyKeOps, flat stacks use per-stack evaluation. Existing finite-fault dispatch
+restrictions remain unchanged.
 External interpolation callbacks keep their existing path and full grid layout.
 
 The triangulation selector sorts the active voxel codes once per surface call and
diff --git a/gempy_engine/API/interp_single/_interp_scalar_field.py b/gempy_engine/API/interp_single/_interp_scalar_field.py
index 1e350a5f..fe76e4b8 100644
--- a/gempy_engine/API/interp_single/_interp_scalar_field.py
+++ b/gempy_engine/API/interp_single/_interp_scalar_field.py
@@ -111,6 +111,11 @@ def _evaluate_sys_eq(eval_input: Union[SolverInput, EvaluatorInput], weights: np
     else:
         exported_fields = generic_evaluator(eval_input, weights, options)
 
+    return _restore_corner_fields(exported_fields, inverse)
+
+
+def _restore_corner_fields(exported_fields: ExportedFields, inverse) -> ExportedFields:
+    """Expand a reduced evaluation before attaching original grid metadata."""
     if inverse is not None:
         for name in ('_scalar_field', '_gx_field', '_gy_field', '_gz_field'):
             values = getattr(exported_fields, name)
diff --git a/gempy_engine/API/interp_single/_stack_ops.py b/gempy_engine/API/interp_single/_stack_ops.py
index 5b6a6c6c..6c0cc45c 100644
--- a/gempy_engine/API/interp_single/_stack_ops.py
+++ b/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
@@ -242,8 +242,7 @@ def _evaluate_optimized(interpolation_inputs: list[InterpolationInput], options:
                          options_per_stack: list[InterpolationOptions] | None = None) -> tuple[list[EvaluatorInput], list[ExportedFields]]:
     from gempy_engine.modules.evaluator.symbolic_evaluator import symbolic_evaluator_optimized_stacked
 
-    # Use the per-stack evaluator so expansion precedes metadata and segmentation.
-    if any(o.evaluation_options.deduplicate_octree_corners for o in (options_per_stack or [options])):
+    if not BackendTensor.use_pykeops:
         return _evaluate(interpolation_inputs, options, solver_inputs, stack_structure,
                          tensor_structs, stack_indices, options_per_stack)
 
@@ -263,15 +262,23 @@ 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):
+        _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
 
diff --git a/gempy_engine/modules/evaluator/symbolic_evaluator.py b/gempy_engine/modules/evaluator/symbolic_evaluator.py
index b0af7557..acde2c9a 100644
--- a/gempy_engine/modules/evaluator/symbolic_evaluator.py
+++ b/gempy_engine/modules/evaluator/symbolic_evaluator.py
@@ -178,8 +178,6 @@ def symbolic_evaluator_optimized_stacked(
     M_sizes = [ei.xyz_to_interpolate.shape[0] for ei in eval_inputs]
     N_sizes = [w.shape[0] for w in weights_list]
 
-    kernel_data_list = []
-
     # 2. Define a small wrapper function for the executor to map over
     def _run_prep(args):
         ei, axis, opt = args
@@ -195,33 +193,40 @@ def _run_prep(args):
     # For now, we take from the first one
     base_options = options_list[0]
 
-    if base_options.compute_scalar is True:
-        prep_tasks = [(ei, None, options_list[idx]) for idx, ei in enumerate(eval_inputs)]
-        with concurrent.futures.ThreadPoolExecutor() as executor:
-            kernel_data_list = list(executor.map(_run_prep, prep_tasks))
+    axes = ([None] if base_options.compute_scalar else []) + ([0, 1, 2] if base_options.compute_scalar_gradient else [])
+    if not axes:
+        raise ValueError("At least one of scalar or scalar gradient must be enabled")
+    tile_factor = len(axes)
+    prep_tasks = [(ei, axis, options_list[idx]) for axis in axes for idx, ei in enumerate(eval_inputs)]
+    with concurrent.futures.ThreadPoolExecutor() as executor:
+        kernel_data_list = list(executor.map(_run_prep, prep_tasks))
+    concat_kernel_data: KernelInput = _build_stacked_kernel_data(kernel_data_list)
 
-        concat_kernel_data: KernelInput = _build_stacked_kernel_data(kernel_data_list)
-        
+    if BackendTensor.engine_backend == gempy_engine.config.AvailableBackends.numpy:
+        from pykeops.numpy import LazyTensor
+    else:
+        from pykeops.torch import LazyTensor
+
+    if base_options.compute_scalar:
+        # The scalar constructor's fault selector assumes one stack. Select
+        # local fault rows here before the shared block-sparse reduction.
+        faults = concat_kernel_data.ref_fault
+        concat_kernel_data.ref_fault = None
         eval_kernel_scalar = create_scalar_kernel(
             concat_kernel_data,
             base_options.kernel_options,
             execution_mode=KernelExecutionMode.SYMBOLIC,
         )
-
-    if base_options.compute_scalar_gradient is True:
-        prep_tasks = []
-        for idx, ei in enumerate(eval_inputs):
-            prep_tasks.append((ei, 0, options_list[idx]))  # X gradient
-        for idx, ei in enumerate(eval_inputs):
-            prep_tasks.append((ei, 1, options_list[idx]))  # Y gradient
-        for idx, ei in enumerate(eval_inputs):
-            prep_tasks.append((ei, 2, options_list[idx]))  # Z gradient
-
-        with concurrent.futures.ThreadPoolExecutor() as executor:
-            kernel_data_list = list(executor.map(_run_prep, prep_tasks))
-
-        concat_kernel_data: KernelInput = _build_stacked_kernel_data(kernel_data_list)
-        
+        if faults is not None:
+            fault_rows = BackendTensor.t.concatenate([
+                BackendTensor.t.concatenate((BackendTensor.t.zeros(n - ei.fault_internal.n_faults, dtype=BackendTensor.dtype_obj),
+                                             BackendTensor.t.ones(ei.fault_internal.n_faults, dtype=BackendTensor.dtype_obj)))
+                for _ in axes for ei, n in zip(eval_inputs, N_sizes)
+            ])
+            fault_selector = LazyTensor(fault_rows.reshape(-1, 1), axis=0)
+            eval_kernel_scalar = eval_kernel_scalar + fault_selector * (faults.faults_i * faults.faults_j).sum(-1)
+
+    if base_options.compute_scalar_gradient:
         eval_kernel_grad = create_grad_kernel(
             concat_kernel_data,
             base_options.kernel_options,
@@ -231,17 +236,16 @@ def _run_prep(args):
     # region kernels
     match (base_options.compute_scalar, base_options.compute_scalar_gradient):
         case (True, True):
-            # Concatenate eval kernel
-            eval_kernel = BackendTensor.t.concatenate([eval_kernel_scalar, eval_kernel_grad], axis=1)
-            tile_factor = 4
+            # LazyTensor index axes cannot be concatenated. Both formulas use
+            # the same stacked data; select scalar for the first output block.
+            scalar_rows = BackendTensor.t.concatenate((BackendTensor.t.ones(sum(M_sizes), dtype=BackendTensor.dtype_obj),
+                                                       -BackendTensor.t.ones(3 * sum(M_sizes), dtype=BackendTensor.dtype_obj)))
+            scalar_selector = LazyTensor(scalar_rows.reshape(-1, 1), axis=1)
+            eval_kernel = scalar_selector.ifelse(eval_kernel_scalar, eval_kernel_grad)
         case (True, False):
             eval_kernel = eval_kernel_scalar
-            tile_factor = 1
         case (False, True):
             eval_kernel = eval_kernel_grad
-            tile_factor = 3
-        case (False, False):
-            raise ValueError("Cannot compute scalar and scalar gradient simultaneously")
     # endregion
 
     M_sizes = M_sizes * tile_factor
diff --git a/tests/test_common/test_modules/test_fused_corner_evaluation.py b/tests/test_common/test_modules/test_fused_corner_evaluation.py
new file mode 100644
index 00000000..56b335ba
--- /dev/null
+++ b/tests/test_common/test_modules/test_fused_corner_evaluation.py
@@ -0,0 +1,165 @@
+from unittest.mock import patch
+
+import numpy as np
+import pytest
+
+from gempy_engine.API.interp_single import _stack_ops
+from gempy_engine.API.interp_single._octree_generation import _generate_corners
+from gempy_engine.core.backend_tensor import BackendTensor
+from gempy_engine.core.data.engine_grid import EngineGrid
+from gempy_engine.core.data.generic_grid import GenericGrid
+from gempy_engine.core.data.kernel_classes.faults import FaultsData
+from gempy_engine.core.data.regular_grid import RegularGrid
+from gempy_engine.modules.evaluator import symbolic_evaluator as symbolic
+from tests.fixtures.simple_models import simple_model_interpolation_input_factory
+from tests.test_common.test_modules.test_octree_optimizations import backend
+
+
+FIELDS = ('_scalar_field', '_gx_field', '_gy_field', '_gz_field')
+
+
+@pytest.mark.parametrize('case', ['plain', 'faults', 'unequal_faults', 'corner_grad', 'fault_grad'])
+@pytest.mark.parametrize('mode', ['scalar', 'gradient', 'combined'])
+def test_fused_corner_evaluation(backend, case, mode, monkeypatch):
+    pytest.importorskip('pykeops')
+    if case in ('corner_grad', 'fault_grad') and backend != 'PYTORCH':
+        pytest.skip('Autograd safeguards require torch')
+    monkeypatch.setattr(BackendTensor, 'use_pykeops', True)
+    t = BackendTensor.t
+    inputs, solvers, structs, options = [], [], [], []
+    targets = []
+    for i in range(3):
+        interp, opt, descriptor = simple_model_interpolation_input_factory()
+        root = RegularGrid([0, 1, 0, 1, 0, 1], [i + 2, 2, 2])
+        corners = _generate_corners(root)
+        if case == 'corner_grad' and i == 2:
+            corners.requires_grad_()
+            targets.append(corners)
+        grid = EngineGrid(octree_grid=root, corners_grid=GenericGrid(corners),
+                          custom_grid=GenericGrid(t.array([[0.2, 0.3, 0.4]])))
+        interp.set_temp_grid(grid)
+        if backend == 'PYTORCH':
+            interp.surface_points.sp_coords.requires_grad_()
+            targets.append(interp.surface_points.sp_coords)
+        if case in ('faults', 'unequal_faults', 'fault_grad'):
+            xyz = t.concatenate((grid.values, interp.all_surface_points.sp_coords))
+            values = xyz[:, 0][None, :] ** 2
+            if backend == 'PYTORCH':
+                values = values.detach()
+            if case == 'unequal_faults' and i == 2:
+                values = BackendTensor.arange(len(xyz), dtype='float64')[None, :] / len(xyz)
+            if case == 'fault_grad' and i == 2:
+                values.requires_grad_()
+                # Spatial-gradient kernels do not depend on the fault values.
+                if mode != 'gradient':
+                    targets.append(values)
+            interp.fault_values = FaultsData(values, values[:, -interp.surface_points.n_points:])
+        # Third stack is ineligible in the plain case, independently of its option.
+        if case == 'plain' and i == 2:
+            grid.corners_grid.values[7] += 0.01
+        opt.evaluation_options.compute_scalar_gradient = mode != 'scalar'
+        opt.evaluation_options.compute_scalar = mode != 'gradient'
+        solver = _stack_ops.input_preprocess_v2(descriptor.tensors_structure, interp)
+        n = (solver.ori_internal.n_orientations_tiled + solver.sp_internal.n_points
+             + opt.kernel_options.n_uni_eq + solver.fault_internal.n_faults)
+        solver.weights_x0 = t.array(np.linspace(0.1 + i, 1.1 + i, n))
+        if backend == 'PYTORCH':
+            solver.weights_x0.requires_grad_()
+            targets.append(solver.weights_x0)
+        inputs.append(interp)
+        solvers.append(solver)
+        structs.append(descriptor.tensors_structure)
+        options.append(opt)
+
+    args = dict(interpolation_inputs=inputs, options=options[0], solver_inputs=solvers,
+                stack_structure=descriptor.stack_structure, tensor_structs=structs,
+                stack_indices=[0, 1, 2], options_per_stack=options)
+    _, expected = _stack_ops._evaluate_optimized(**args)
+    _, serial = _stack_ops._evaluate(**args)
+    fields = FIELDS if mode == 'combined' else FIELDS[1:] if mode == 'gradient' else FIELDS[:1]
+    for fused_result, serial_result in zip(expected, serial):
+        for name in fields:
+            np.testing.assert_allclose(t.to_numpy(getattr(fused_result, name)),
+                                       t.to_numpy(getattr(serial_result, name)), atol=1e-9, rtol=1e-9)
+    for i, opt in enumerate(options):
+        opt.evaluation_options.deduplicate_octree_corners = i != 1
+
+    if backend == 'numpy':
+        from pykeops.numpy import LazyTensor
+    else:
+        from pykeops.torch import LazyTensor
+    reductions = []
+    lazy_sum = LazyTensor.sum
+
+    def record_sum(self, *args, **kwargs):
+        if kwargs.get('axis') == 0:
+            reductions.append(kwargs)
+        return lazy_sum(self, *args, **kwargs)
+
+    with patch.object(symbolic, 'symbolic_evaluator_optimized_stacked',
+                      wraps=symbolic.symbolic_evaluator_optimized_stacked) as fused, \
+            patch.object(LazyTensor, 'sum', new=record_sum), \
+            patch.object(_stack_ops, '_evaluate', side_effect=AssertionError('Must stay fused')):
+        originals, actual = _stack_ops._evaluate_optimized(**args)
+    assert fused.call_count == 1
+    assert len(reductions) == 1
+    assert reductions[0]['backend'] == 'CPU'
+    assert reductions[0]['ranges'] is not None
+    assert fused.call_args.kwargs['options_list'] is options
+    reduced = fused.call_args.kwargs['eval_inputs']
+    assert len({len(e.xyz_to_interpolate) for e in reduced}) == 3
+    for i, (original, view, result, reference) in enumerate(zip(originals, reduced, actual, expected)):
+        eligible = i == 0 or (i == 2 and case == 'faults')
+        assert (len(view.xyz_to_interpolate) < len(original.xyz_to_interpolate)) == eligible
+        assert (view is not original) == eligible
+        assert original.solver_input is solvers[i]
+        assert fused.call_args.kwargs['weights_list'][i] is solvers[i].weights_x0
+        assert view.sp_internal is original.sp_internal
+        assert view.ori_internal is original.ori_internal
+        if original.fault_internal.n_faults:
+            assert original.fault_internal is inputs[i].fault_values
+            assert original.fault_internal.fault_values_everywhere.shape[1] == len(original.xyz_to_interpolate)
+            assert view.fault_internal.fault_values_everywhere.shape[1] == len(view.xyz_to_interpolate)
+            assert (view.fault_internal is not original.fault_internal) == eligible
+        np.testing.assert_array_equal(t.to_numpy(original.xyz_to_interpolate),
+                                      t.to_numpy(t.concatenate((inputs[i].grid.values,
+                                                               inputs[i].all_surface_points.sp_coords))))
+        assert result.grid_size == inputs[i].grid.len_all_grids
+        assert result.n_points_per_surface is original._n_points_per_surface
+        assert result.slice_feature is original._slice_feature
+        assert result.debug is solvers[i].debug
+        for name in FIELDS:
+            a, b = getattr(result, name), getattr(reference, name)
+            if b is None:
+                assert a is None
+            else:
+                np.testing.assert_allclose(t.to_numpy(a), t.to_numpy(b), atol=1e-9, rtol=1e-9)
+        np.testing.assert_allclose(t.to_numpy(result.scalar_field_at_surface_points),
+                                   t.to_numpy(reference.scalar_field_at_surface_points), atol=1e-9, rtol=1e-9)
+
+    if backend == 'PYTORCH':
+        import torch
+
+        def loss(results):
+            return sum((getattr(result, name) * (torch.arange(len(result._scalar_field),
+                                                             dtype=torch.float64) + 1)).sum()
+                       for result in results for name in fields)
+
+        old_grads = torch.autograd.grad(loss(expected), targets, retain_graph=True)
+        serial_grads = torch.autograd.grad(loss(serial), targets, retain_graph=True)
+        new_grads = torch.autograd.grad(loss(actual), targets)
+        for old, new, reference in zip(old_grads, new_grads, serial_grads):
+            torch.testing.assert_close(old, new, atol=1e-7, rtol=1e-8)
+            torch.testing.assert_close(new, reference, atol=1e-7, rtol=1e-8)
+
+
+@pytest.mark.parametrize('deduplicate', [False, True])
+def test_non_pykeops_dispatch(backend, deduplicate):
+    interp, options, descriptor = simple_model_interpolation_input_factory()
+    options.evaluation_options.deduplicate_octree_corners = deduplicate
+    with patch.object(_stack_ops, '_evaluate', return_value=([], [])) as evaluate, \
+            patch.object(symbolic, 'symbolic_evaluator_optimized_stacked',
+                         side_effect=AssertionError('PyKeOps is disabled')):
+        assert _stack_ops._evaluate_optimized([interp], options, [], descriptor.stack_structure,
+                                              [], [0]) == ([], [])
+    assert evaluate.call_count == 1

From 55d09948730bed288d0b17cc5473c950d53c5e0a Mon Sep 17 00:00:00 2001
From: Miguel de la Varga <miguel@terranigma-solutions.com>
Date: Mon, 21 Sep 2026 15:06:00 +0200
Subject: [PATCH 4/9] [ENH] Add optional unique-edge quad triangulation

---
 docs/octree_refinement.md                     |  33 +++
 .../multi_scalar_dual_contouring.py           |   4 +-
 .../core/data/dual_contouring_data.py         |   5 +
 .../core/data/options/evaluation_options.py   |   6 +
 .../dual_contouring/dual_contouring_v2.py     |  11 +
 .../dual_contouring/quad_triangulation.py     |  99 +++++++++
 .../test_modules/test_quad_triangulation.py   | 209 ++++++++++++++++++
 7 files changed, 366 insertions(+), 1 deletion(-)
 create mode 100644 gempy_engine/modules/dual_contouring/quad_triangulation.py
 create mode 100644 tests/test_common/test_modules/test_quad_triangulation.py

diff --git a/docs/octree_refinement.md b/docs/octree_refinement.md
index 0cde850e..76be72e9 100644
--- a/docs/octree_refinement.md
+++ b/docs/octree_refinement.md
@@ -124,3 +124,36 @@ reuses that lookup across all six edge cases. It preserves triangle order,
 neighbor filtering, and normal correction; it does not rewrite edge topology or
 change vertex generation. Its lookup is local to the surface call and safe for
 parallel surface processing.
+
+### 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. `triangulation_sort_once` applies only to the legacy
+method; quad mode already 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.
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 ccd1a9a2..4c86da70 100644
--- a/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py
+++ b/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py
@@ -153,7 +153,9 @@ def dual_contouring_multi_scalar(
                 n_surfaces_to_export=n_scalar_field,
                 tree_depth=options.number_octree_levels,
                 base_number=base_number,
-                triangulation_sort_once=options.evaluation_options.triangulation_sort_once
+                triangulation_sort_once=options.evaluation_options.triangulation_sort_once,
+                triangulation_method=options.evaluation_options.triangulation_method,
+                generated_cell_coordinates=left_right_codes
             )
 
             dc_data_per_surface_all.append(dc_data_per_surface)
diff --git a/gempy_engine/core/data/dual_contouring_data.py b/gempy_engine/core/data/dual_contouring_data.py
index 672b57c9..557f3001 100644
--- a/gempy_engine/core/data/dual_contouring_data.py
+++ b/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:
@@ -28,6 +30,9 @@ class DualContouringData:
     extra_edge_normals: Optional[np.ndarray] = None   # (n_valid_voxels, K, 3)
     extra_weights: Optional[np.ndarray] = None        # (n_valid_voxels, K)
     triangulation_sort_once: bool = False
+    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):
diff --git a/gempy_engine/core/data/options/evaluation_options.py b/gempy_engine/core/data/options/evaluation_options.py
index f5912bc4..73dc625d 100644
--- a/gempy_engine/core/data/options/evaluation_options.py
+++ b/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
@@ -32,6 +37,7 @@ class EvaluationOptions:
     octree_refinement_mode: OctreeRefinementMode = OctreeRefinementMode.FAST
     deduplicate_octree_corners: bool = False  #: Evaluate unique corners, then restore the original row layout.
     triangulation_sort_once: bool = False  #: Reuse one voxel-code sort across the six edge cases per surface.
+    triangulation_method: TriangulationMethod = TriangulationMethod.LEGACY
     
     mesh_extraction: bool = True
     mesh_extraction_masking_options: MeshExtractionMaskingOptions = MeshExtractionMaskingOptions.INTERSECT
diff --git a/gempy_engine/modules/dual_contouring/dual_contouring_v2.py b/gempy_engine/modules/dual_contouring/dual_contouring_v2.py
index cc28d1cd..42cc49ca 100644
--- a/gempy_engine/modules/dual_contouring/dual_contouring_v2.py
+++ b/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
 
diff --git a/gempy_engine/modules/dual_contouring/quad_triangulation.py b/gempy_engine/modules/dual_contouring/quad_triangulation.py
new file mode 100644
index 00000000..4fd6c34a
--- /dev/null
+++ b/gempy_engine/modules/dual_contouring/quad_triangulation.py
@@ -0,0 +1,99 @@
+"""Same-level dual quads keyed by (lower integer edge endpoint, direction)."""
+
+from ...config import AvailableBackends
+from ...core.backend_tensor import BackendTensor
+from ._aux import _correct_normals
+from .fancy_triangulation import _get_pack_factors
+
+
+def triangulate_quads(coordinates, valid_edges, voxel_normals, vertices, domain_shape,
+                      generated_coordinates=None, report=None):
+    """Emit two triangles per complete quad, never partial triangles.
+
+    Coordinates include retained non-crossing cells; vertices include only cells
+    with any valid edge, in the same order. Crossing flags come from the existing
+    tolerant intersection rule, not a new sign test. Conflicting flags on shared
+    edges are rejected. Missing-cell counts are incidences; edge categories can
+    overlap. Without pre-mask coordinates, interior omissions are 'unknown'.
+    This does not stitch different octree levels or resolve ambiguous topology.
+    """
+    t = BackendTensor.t
+    report = {} if report is None else report
+    report.clear()
+    report.update(crossing_edge_count=0, quad_count=0, missing_incident_cell_count=0,
+                  physical_boundary_edge_count=0, mask_boundary_edge_count=0,
+                  internal_refinement_boundary_edge_count=0, unknown_boundary_edge_count=0)
+    empty = t.zeros((0, 3), dtype='int64')
+    factors = _get_pack_factors(*domain_shape)
+    bounds = t.array(domain_shape, dtype='int64')
+    if coordinates.shape != (len(valid_edges), 3) or valid_edges.shape != (len(coordinates), 12):
+        raise ValueError("Expected cell coordinates (N, 3) and crossing flags (N, 12)")
+    if t.any((coordinates < 0) | (coordinates >= bounds)):
+        raise ValueError("Cell coordinates must lie inside the theoretical domain")
+    codes = (coordinates * factors).sum(axis=1)
+    order = t.argsort(codes)
+    sorted_codes = codes[order]
+    if t.any(sorted_codes[1:] == sorted_codes[:-1]):
+        raise ValueError("Duplicate cell coordinates are not supported")
+    active = t.any(valid_edges, axis=1)
+    if len(vertices) != int(active.sum()):
+        raise ValueError("Expected one vertex per crossing cell")
+    if len(coordinates) == 0:
+        return empty
+
+    # Edge numbering matches find_intersection_on_edge: x, then y, then z.
+    offsets = t.array([[0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1],
+                       [0, 0, 0], [0, 0, 1], [1, 0, 0], [1, 0, 1],
+                       [0, 0, 0], [0, 1, 0], [1, 0, 0], [1, 1, 0]], dtype='int64')
+    directions = t.array([0] * 4 + [1] * 4 + [2] * 4, dtype='int64')
+    origins = coordinates[:, None, :] + offsets[None, :, :]
+    ids = t.concatenate((origins, t.zeros(origins.shape[:2] + (1,), dtype='int64') +
+                         directions[None, :, None]), axis=2).reshape(-1, 4)
+    if BackendTensor.engine_backend == AvailableBackends.PYTORCH:
+        unique, inverse, counts = t.unique(ids, dim=0, sorted=True, return_inverse=True, return_counts=True)
+    else:
+        unique, inverse, counts = t.unique(ids, axis=0, return_inverse=True, return_counts=True)
+    crossings = t.bincount(inverse[valid_edges.reshape(-1)], minlength=len(unique))
+    if t.any((crossings != 0) & (crossings != counts)):
+        raise ValueError("Inconsistent crossing flags for duplicate canonical edges")
+    edges = unique[crossings > 0]
+    report['crossing_edge_count'] = len(edges)
+    if len(edges) == 0:
+        return empty
+
+    # Cyclic incident cells; the fixed 1--3 diagonal matches legacy complete quads.
+    incident_offsets = t.array([
+        [[0, -1, -1], [0, -1, 0], [0, 0, 0], [0, 0, -1]],
+        [[-1, 0, -1], [-1, 0, 0], [0, 0, 0], [0, 0, -1]],
+        [[-1, -1, 0], [-1, 0, 0], [0, 0, 0], [0, -1, 0]]
+    ], dtype='int64')
+    incident = edges[:, None, :3] + incident_offsets[edges[:, 3]]
+    inside = t.all((incident >= 0) & (incident < bounds), axis=2)
+    keys = (incident * factors).sum(axis=2)
+    positions = t.clip(t.searchsorted(sorted_codes, keys.reshape(-1)), 0, len(codes) - 1).reshape(-1, 4)
+    found = inside & (sorted_codes[positions] == keys)
+    cells = order[positions]
+    complete = t.all(found, axis=1)
+    report['quad_count'] = int(complete.sum())
+    report['missing_incident_cell_count'] = int((~found).sum())
+    report['physical_boundary_edge_count'] = int(t.any(~inside, axis=1).sum())
+    missing_inside = inside & ~found
+    if generated_coordinates is None:
+        report['unknown_boundary_edge_count'] = int(t.any(missing_inside, axis=1).sum())
+    else:
+        generated_codes = (generated_coordinates * factors).sum(axis=1)
+        generated = t.isin(keys, generated_codes)
+        report['mask_boundary_edge_count'] = int(t.any(missing_inside & generated, axis=1).sum())
+        report['internal_refinement_boundary_edge_count'] = int(t.any(missing_inside & ~generated, axis=1).sum())
+    if not t.any(complete):
+        return empty
+
+    # All four cells of a consistent crossing edge have a QEF vertex.
+    vertex_ids = t.cumsum(active, axis=0) - 1
+    quads = vertex_ids[cells[complete]]
+    split = t.array([[0, 1, 3], [2, 3, 1]], dtype='int64')
+    triangles = quads[:, split].reshape(-1, 3)
+    local_edges = t.array([[3, 2, 0, 1], [7, 6, 4, 5], [11, 10, 8, 9]], dtype='int64')
+    reference = voxel_normals[cells[complete], local_edges[edges[complete, 3]]].sum(axis=1)
+    reference = t.repeat(reference, 2, axis=0)
+    return _correct_normals(vertices, triangles, reference)[0]
diff --git a/tests/test_common/test_modules/test_quad_triangulation.py b/tests/test_common/test_modules/test_quad_triangulation.py
new file mode 100644
index 00000000..c7ed79ac
--- /dev/null
+++ b/tests/test_common/test_modules/test_quad_triangulation.py
@@ -0,0 +1,209 @@
+from itertools import product
+
+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 import InterpolationOptions
+from gempy_engine.core.data.options.evaluation_options import TriangulationMethod
+from gempy_engine.modules.dual_contouring.dual_contouring_interface import find_intersection_on_edge
+from gempy_engine.modules.dual_contouring.fancy_triangulation import triangulate
+from gempy_engine.modules.dual_contouring.quad_triangulation import triangulate_quads
+
+
+@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.COMPUTE_GRADS, BackendTensor.pykeops_enabled)
+    BackendTensor._change_backend(AvailableBackends[request.param], use_gpu=False,
+                                 dtype='float64', use_pykeops=False, grads=True)
+    BackendTensor.pykeops_enabled = False
+    yield request.param
+    BackendTensor._change_backend(old[0], use_gpu=old[1], dtype=old[2], use_pykeops=old[3], grads=old[4])
+    BackendTensor.pykeops_enabled = old[5]
+
+
+def plane_data(coords, axis=0, sign=1, iso=1.4):
+    t = BackendTensor.t
+    coords = t.array(coords, dtype='int64')
+    corners = coords[:, None, :] + t.array(list(product((0, 1), repeat=3)), dtype='float64')
+    _, valid = find_intersection_on_edge(corners.reshape(-1, 3),
+                                         sign * corners[:, :, axis].reshape(-1),
+                                         t.array([sign * iso], dtype='float64'))
+    valid = valid.reshape(-1, 12)
+    normals = t.zeros((len(coords), 12, 3), dtype='float64')
+    normals[:, :, axis] = sign
+    normals[~valid] = 0
+    active = t.any(valid, axis=1)
+    vertices = t.array(coords[active], dtype='float64') + 0.5
+    vertices[:, axis] = iso
+    return coords, valid, normals, vertices
+
+
+def canonical_faces(faces):
+    faces = np.sort(faces, axis=1)
+    return faces[np.lexsort(faces.T[::-1])]
+
+
+def oriented_faces(faces):
+    rotated = np.take_along_axis(faces, (np.argmin(faces, axis=1)[:, None] + np.arange(3)) % 3, axis=1)
+    return rotated[np.lexsort(rotated.T[::-1])]
+
+
+@pytest.mark.parametrize('axis', [0, 1, 2])
+@pytest.mark.parametrize('sign', [-1, 1])
+def test_complete_quads_legacy_parity_and_winding(backend, axis, sign):
+    t = BackendTensor.t
+    coords = np.array(list(product(range(3), range(4), range(5))), dtype=np.int64)
+    np.random.default_rng(42).shuffle(coords)
+    coords, valid, normals, vertices = plane_data(coords, axis, sign)
+    report = {}
+    faces = triangulate_quads(coords, valid, normals, vertices, (3, 4, 5), report=report)
+    active = t.any(valid, axis=1)
+    legacy = triangulate(coords[active], valid[active], 1, normals[active], vertices, (3, 4, 5))
+    faces, legacy, vertices = map(t.to_numpy, (faces, legacy, vertices))
+    np.testing.assert_array_equal(canonical_faces(faces), canonical_faces(legacy))
+    np.testing.assert_array_equal(oriented_faces(faces), oriented_faces(legacy))
+    assert len(faces) == 2 * report['quad_count'] > 0
+    assert len(np.unique(np.sort(faces, axis=1), axis=0)) == len(faces)
+    normal = np.cross(vertices[faces[:, 1]] - vertices[faces[:, 0]],
+                      vertices[faces[:, 2]] - vertices[faces[:, 0]])
+    assert np.all(normal[:, axis] * sign > 0)
+    assert faces.dtype == np.int64
+    assert report['physical_boundary_edge_count'] > 0
+    assert report['unknown_boundary_edge_count'] == 0
+
+
+@pytest.mark.parametrize('kind', ['masked', 'refinement', 'unknown'])
+def test_sparse_missing_quad_classification(backend, kind):
+    t = BackendTensor.t
+    full = np.array(list(product([1], [1, 2], [1, 2])), dtype=np.int64)
+    coords, valid, normals, vertices = plane_data(full[:-1])
+    generated = None if kind == 'unknown' else t.array(full if kind == 'masked' else full[:-1], dtype='int64')
+    report = {}
+    faces = triangulate_quads(coords, valid, normals, vertices, (4, 4, 4), generated, report)
+    assert faces.shape == (0, 3)  # Three vertices must not produce a partial quad.
+    assert report['quad_count'] == 0
+    assert report['physical_boundary_edge_count'] == 0
+    key = {'masked': 'mask_boundary_edge_count', 'refinement': 'internal_refinement_boundary_edge_count',
+           'unknown': 'unknown_boundary_edge_count'}[kind]
+    assert report[key] > 0
+    if kind == 'masked':
+        # Some other neighbors were never generated, rather than masked out.
+        assert report['internal_refinement_boundary_edge_count'] > 0
+
+
+def test_sparse_large_domain_and_row_order(backend):
+    t = BackendTensor.t
+    # Large theoretical extent, only four generated cells and no dense domain allocation.
+    coords = np.array(list(product([2 ** 32], [1, 2], [1, 2])), dtype=np.int64)
+    results = []
+    for order in ([0, 1, 2, 3], [3, 0, 2, 1]):
+        cells, valid, normals, vertices = plane_data(coords[order], iso=2 ** 32 + 0.4)
+        report = {}
+        faces = triangulate_quads(cells, valid, normals, vertices, (2 ** 32 + 2, 4, 4), cells, report)
+        assert faces.shape == (2, 3)
+        assert report == dict(crossing_edge_count=9, quad_count=1, missing_incident_cell_count=20,
+                              physical_boundary_edge_count=0, mask_boundary_edge_count=0,
+                              internal_refinement_boundary_edge_count=8, unknown_boundary_edge_count=0)
+        results.append(t.to_numpy(vertices[faces]))
+    np.testing.assert_array_equal(*results)
+
+
+def test_winding_uses_shared_edge_gradient(backend):
+    t = BackendTensor.t
+    coords = np.array(list(product([1], [1, 2], [1, 2])), dtype=np.int64)
+    coords, valid, normals, vertices = plane_data(coords)
+    normals[:, :, 0] = -100
+    # The four copies of the central crossing point have positive x gradients.
+    for cell, edge in enumerate([3, 2, 1, 0]):
+        normals[cell, edge, 0] = 1
+    faces = triangulate_quads(coords, valid, normals, vertices, (4, 4, 4))
+    xyz = t.to_numpy(vertices[faces])
+    assert np.all(np.cross(xyz[:, 1] - xyz[:, 0], xyz[:, 2] - xyz[:, 0])[:, 0] > 0)
+
+
+@pytest.mark.parametrize('case', ['empty', 'no_crossings', 'isolated', 'boundary'])
+def test_empty_and_incomplete(backend, case):
+    t = BackendTensor.t
+    coords = [] if case == 'empty' else [[0, 0, 0] if case == 'boundary' else [1, 1, 1]]
+    coords = np.array(coords, dtype=np.int64).reshape(-1, 3)
+    coords, valid, normals, vertices = plane_data(coords, iso=10 if case == 'no_crossings' else 0.4 if case == 'boundary' else 1.4)
+    report = {'stale': True}
+    faces = triangulate_quads(coords, valid, normals, vertices, (4, 4, 4), report=report)
+    assert faces.shape == (0, 3)
+    assert t.to_numpy(faces).dtype == np.int64
+    assert 'stale' not in report
+    if case == 'boundary':
+        assert report['physical_boundary_edge_count'] > 0
+    if backend == 'PYTORCH':
+        assert faces.device == coords.device
+
+
+def test_duplicate_consistency(backend):
+    coords = np.array(list(product([1], [1, 2], [1, 2])), dtype=np.int64)
+    coords, valid, normals, vertices = plane_data(coords)
+    valid[0, 3] = False  # Shared central x edge disagrees with three other cells.
+    with pytest.raises(ValueError, match='Inconsistent crossing'):
+        triangulate_quads(coords, valid, normals, vertices, (4, 4, 4))
+    # Also check duplicates from retained cells that have no QEF vertex at all.
+    valid[0, :] = False
+    with pytest.raises(ValueError, match='Inconsistent crossing'):
+        triangulate_quads(coords, valid, normals, vertices[1:], (4, 4, 4))
+    coords, valid, normals, vertices = plane_data([[1, 1, 1], [1, 1, 1]])
+    with pytest.raises(ValueError, match='Duplicate cell'):
+        triangulate_quads(coords, valid, normals, vertices, (4, 4, 4))
+
+
+def test_tolerant_crossing_rule(backend):
+    coords = np.array(list(product([1], [1, 2], [1, 2])), dtype=np.int64)
+    # Both endpoint values are below the isovalue, but inside the shipped tolerance.
+    coords, valid, normals, vertices = plane_data(coords, iso=2.005)
+    assert BackendTensor.t.all(valid[:, :4])
+    faces = triangulate_quads(coords, valid, normals, vertices, (4, 4, 4))
+    assert faces.shape == (2, 3)
+
+
+def test_bounds_and_overflow(backend):
+    coords, valid, normals, vertices = plane_data([[-1, 1, 1]])
+    with pytest.raises(ValueError, match='theoretical domain'):
+        triangulate_quads(coords, valid, normals, vertices, (4, 4, 4))
+    with pytest.raises(OverflowError):
+        triangulate_quads(coords, valid, normals, vertices, (2 ** 32, 2 ** 32, 2))
+
+
+def test_selector_serialization():
+    options = InterpolationOptions.from_args(range=1., c_o=1.)
+    assert options.evaluation_options.triangulation_method is TriangulationMethod.LEGACY
+    options.evaluation_options.triangulation_method = TriangulationMethod.QUADS
+    options.evaluation_options.triangulation_sort_once = True
+    restored = InterpolationOptions.model_validate_json(options.model_dump_json())
+    assert restored.evaluation_options.triangulation_method is TriangulationMethod.QUADS
+    assert restored.evaluation_options.triangulation_sort_once
+
+
+def test_public_model_mesh_parity(backend, monkeypatch):
+    from gempy_engine.API.model.model_api import compute_model
+    from tests.fixtures.simple_models import simple_model_interpolation_input_factory
+
+    monkeypatch.setenv('GEMPY_FLAT_STACKS', 'False')
+    monkeypatch.setenv('GEMPY_SKIP_TRIANGULATION', '0')
+    results = []
+    for method in (TriangulationMethod.LEGACY, TriangulationMethod.QUADS):
+        interp, options, descriptor = simple_model_interpolation_input_factory()
+        options.evaluation_options.number_octree_levels = 2
+        options.evaluation_options.triangulation_method = method
+        results.append(compute_model(interp, options, descriptor))
+    assert len(results[0].dc_meshes) == len(results[1].dc_meshes) > 0
+    for old, new in zip(results[0].dc_meshes, results[1].dc_meshes):
+        np.testing.assert_allclose(old.vertices, new.vertices, atol=1e-8)
+        np.testing.assert_array_equal(canonical_faces(old.edges), canonical_faces(new.edges))
+        np.testing.assert_array_equal(oriented_faces(old.edges), oriented_faces(new.edges))
+        assert len(new.edges) > 0
+        assert new.dc_data.triangulation_method is TriangulationMethod.QUADS
+        assert new.dc_data.generated_cell_coordinates is not None
+        assert new.dc_data.triangulation_report['unknown_boundary_edge_count'] == 0
+        assert new.dc_data.triangulation_report['quad_count'] * 2 == len(new.edges)

From 7801af78502ac8b219f27f4834a696d92b8e00e4 Mon Sep 17 00:00:00 2001
From: Miguel de la Varga <miguel@terranigma-solutions.com>
Date: Mon, 21 Sep 2026 15:41:46 +0200
Subject: [PATCH 5/9] [FIX] Keep fused corner deduplication scoped to
 evaluation views

---
 .../modules/evaluator/symbolic_evaluator.py   |  62 +++++-----
 .../test_fused_corner_evaluation.py           | 109 +++---------------
 2 files changed, 42 insertions(+), 129 deletions(-)

diff --git a/gempy_engine/modules/evaluator/symbolic_evaluator.py b/gempy_engine/modules/evaluator/symbolic_evaluator.py
index acde2c9a..c3588281 100644
--- a/gempy_engine/modules/evaluator/symbolic_evaluator.py
+++ b/gempy_engine/modules/evaluator/symbolic_evaluator.py
@@ -178,6 +178,8 @@ def symbolic_evaluator_optimized_stacked(
     M_sizes = [ei.xyz_to_interpolate.shape[0] for ei in eval_inputs]
     N_sizes = [w.shape[0] for w in weights_list]
 
+    kernel_data_list = []
+
     # 2. Define a small wrapper function for the executor to map over
     def _run_prep(args):
         ei, axis, opt = args
@@ -193,40 +195,33 @@ def _run_prep(args):
     # For now, we take from the first one
     base_options = options_list[0]
 
-    axes = ([None] if base_options.compute_scalar else []) + ([0, 1, 2] if base_options.compute_scalar_gradient else [])
-    if not axes:
-        raise ValueError("At least one of scalar or scalar gradient must be enabled")
-    tile_factor = len(axes)
-    prep_tasks = [(ei, axis, options_list[idx]) for axis in axes for idx, ei in enumerate(eval_inputs)]
-    with concurrent.futures.ThreadPoolExecutor() as executor:
-        kernel_data_list = list(executor.map(_run_prep, prep_tasks))
-    concat_kernel_data: KernelInput = _build_stacked_kernel_data(kernel_data_list)
+    if base_options.compute_scalar is True:
+        prep_tasks = [(ei, None, options_list[idx]) for idx, ei in enumerate(eval_inputs)]
+        with concurrent.futures.ThreadPoolExecutor() as executor:
+            kernel_data_list = list(executor.map(_run_prep, prep_tasks))
 
-    if BackendTensor.engine_backend == gempy_engine.config.AvailableBackends.numpy:
-        from pykeops.numpy import LazyTensor
-    else:
-        from pykeops.torch import LazyTensor
+        concat_kernel_data: KernelInput = _build_stacked_kernel_data(kernel_data_list)
 
-    if base_options.compute_scalar:
-        # The scalar constructor's fault selector assumes one stack. Select
-        # local fault rows here before the shared block-sparse reduction.
-        faults = concat_kernel_data.ref_fault
-        concat_kernel_data.ref_fault = None
         eval_kernel_scalar = create_scalar_kernel(
             concat_kernel_data,
             base_options.kernel_options,
             execution_mode=KernelExecutionMode.SYMBOLIC,
         )
-        if faults is not None:
-            fault_rows = BackendTensor.t.concatenate([
-                BackendTensor.t.concatenate((BackendTensor.t.zeros(n - ei.fault_internal.n_faults, dtype=BackendTensor.dtype_obj),
-                                             BackendTensor.t.ones(ei.fault_internal.n_faults, dtype=BackendTensor.dtype_obj)))
-                for _ in axes for ei, n in zip(eval_inputs, N_sizes)
-            ])
-            fault_selector = LazyTensor(fault_rows.reshape(-1, 1), axis=0)
-            eval_kernel_scalar = eval_kernel_scalar + fault_selector * (faults.faults_i * faults.faults_j).sum(-1)
-
-    if base_options.compute_scalar_gradient:
+
+    if base_options.compute_scalar_gradient is True:
+        prep_tasks = []
+        for idx, ei in enumerate(eval_inputs):
+            prep_tasks.append((ei, 0, options_list[idx]))  # X gradient
+        for idx, ei in enumerate(eval_inputs):
+            prep_tasks.append((ei, 1, options_list[idx]))  # Y gradient
+        for idx, ei in enumerate(eval_inputs):
+            prep_tasks.append((ei, 2, options_list[idx]))  # Z gradient
+
+        with concurrent.futures.ThreadPoolExecutor() as executor:
+            kernel_data_list = list(executor.map(_run_prep, prep_tasks))
+
+        concat_kernel_data: KernelInput = _build_stacked_kernel_data(kernel_data_list)
+
         eval_kernel_grad = create_grad_kernel(
             concat_kernel_data,
             base_options.kernel_options,
@@ -236,16 +231,17 @@ def _run_prep(args):
     # region kernels
     match (base_options.compute_scalar, base_options.compute_scalar_gradient):
         case (True, True):
-            # LazyTensor index axes cannot be concatenated. Both formulas use
-            # the same stacked data; select scalar for the first output block.
-            scalar_rows = BackendTensor.t.concatenate((BackendTensor.t.ones(sum(M_sizes), dtype=BackendTensor.dtype_obj),
-                                                       -BackendTensor.t.ones(3 * sum(M_sizes), dtype=BackendTensor.dtype_obj)))
-            scalar_selector = LazyTensor(scalar_rows.reshape(-1, 1), axis=1)
-            eval_kernel = scalar_selector.ifelse(eval_kernel_scalar, eval_kernel_grad)
+            # Concatenate eval kernel
+            eval_kernel = BackendTensor.t.concatenate([eval_kernel_scalar, eval_kernel_grad], axis=1)
+            tile_factor = 4
         case (True, False):
             eval_kernel = eval_kernel_scalar
+            tile_factor = 1
         case (False, True):
             eval_kernel = eval_kernel_grad
+            tile_factor = 3
+        case (False, False):
+            raise ValueError("Cannot compute scalar and scalar gradient simultaneously")
     # endregion
 
     M_sizes = M_sizes * tile_factor
diff --git a/tests/test_common/test_modules/test_fused_corner_evaluation.py b/tests/test_common/test_modules/test_fused_corner_evaluation.py
index 56b335ba..f095d48f 100644
--- a/tests/test_common/test_modules/test_fused_corner_evaluation.py
+++ b/tests/test_common/test_modules/test_fused_corner_evaluation.py
@@ -8,64 +8,31 @@
 from gempy_engine.core.backend_tensor import BackendTensor
 from gempy_engine.core.data.engine_grid import EngineGrid
 from gempy_engine.core.data.generic_grid import GenericGrid
-from gempy_engine.core.data.kernel_classes.faults import FaultsData
 from gempy_engine.core.data.regular_grid import RegularGrid
 from gempy_engine.modules.evaluator import symbolic_evaluator as symbolic
 from tests.fixtures.simple_models import simple_model_interpolation_input_factory
 from tests.test_common.test_modules.test_octree_optimizations import backend
 
 
-FIELDS = ('_scalar_field', '_gx_field', '_gy_field', '_gz_field')
-
-
-@pytest.mark.parametrize('case', ['plain', 'faults', 'unequal_faults', 'corner_grad', 'fault_grad'])
-@pytest.mark.parametrize('mode', ['scalar', 'gradient', 'combined'])
-def test_fused_corner_evaluation(backend, case, mode, monkeypatch):
+def test_fused_corner_evaluation_uses_reduced_views(backend, monkeypatch):
     pytest.importorskip('pykeops')
-    if case in ('corner_grad', 'fault_grad') and backend != 'PYTORCH':
-        pytest.skip('Autograd safeguards require torch')
     monkeypatch.setattr(BackendTensor, 'use_pykeops', True)
     t = BackendTensor.t
     inputs, solvers, structs, options = [], [], [], []
-    targets = []
     for i in range(3):
         interp, opt, descriptor = simple_model_interpolation_input_factory()
-        root = RegularGrid([0, 1, 0, 1, 0, 1], [i + 2, 2, 2])
+        root = RegularGrid([10, 13, -4, -2, 20, 24], [i + 2, 2, 2])
         corners = _generate_corners(root)
-        if case == 'corner_grad' and i == 2:
-            corners.requires_grad_()
-            targets.append(corners)
         grid = EngineGrid(octree_grid=root, corners_grid=GenericGrid(corners),
-                          custom_grid=GenericGrid(t.array([[0.2, 0.3, 0.4]])))
+                          custom_grid=GenericGrid(corners[:2]))
         interp.set_temp_grid(grid)
-        if backend == 'PYTORCH':
-            interp.surface_points.sp_coords.requires_grad_()
-            targets.append(interp.surface_points.sp_coords)
-        if case in ('faults', 'unequal_faults', 'fault_grad'):
-            xyz = t.concatenate((grid.values, interp.all_surface_points.sp_coords))
-            values = xyz[:, 0][None, :] ** 2
-            if backend == 'PYTORCH':
-                values = values.detach()
-            if case == 'unequal_faults' and i == 2:
-                values = BackendTensor.arange(len(xyz), dtype='float64')[None, :] / len(xyz)
-            if case == 'fault_grad' and i == 2:
-                values.requires_grad_()
-                # Spatial-gradient kernels do not depend on the fault values.
-                if mode != 'gradient':
-                    targets.append(values)
-            interp.fault_values = FaultsData(values, values[:, -interp.surface_points.n_points:])
-        # Third stack is ineligible in the plain case, independently of its option.
-        if case == 'plain' and i == 2:
-            grid.corners_grid.values[7] += 0.01
-        opt.evaluation_options.compute_scalar_gradient = mode != 'scalar'
-        opt.evaluation_options.compute_scalar = mode != 'gradient'
+        opt.evaluation_options.compute_scalar = True
+        opt.evaluation_options.compute_scalar_gradient = False
+        opt.evaluation_options.deduplicate_octree_corners = True
         solver = _stack_ops.input_preprocess_v2(descriptor.tensors_structure, interp)
         n = (solver.ori_internal.n_orientations_tiled + solver.sp_internal.n_points
              + opt.kernel_options.n_uni_eq + solver.fault_internal.n_faults)
         solver.weights_x0 = t.array(np.linspace(0.1 + i, 1.1 + i, n))
-        if backend == 'PYTORCH':
-            solver.weights_x0.requires_grad_()
-            targets.append(solver.weights_x0)
         inputs.append(interp)
         solvers.append(solver)
         structs.append(descriptor.tensors_structure)
@@ -74,53 +41,22 @@ def test_fused_corner_evaluation(backend, case, mode, monkeypatch):
     args = dict(interpolation_inputs=inputs, options=options[0], solver_inputs=solvers,
                 stack_structure=descriptor.stack_structure, tensor_structs=structs,
                 stack_indices=[0, 1, 2], options_per_stack=options)
-    _, expected = _stack_ops._evaluate_optimized(**args)
-    _, serial = _stack_ops._evaluate(**args)
-    fields = FIELDS if mode == 'combined' else FIELDS[1:] if mode == 'gradient' else FIELDS[:1]
-    for fused_result, serial_result in zip(expected, serial):
-        for name in fields:
-            np.testing.assert_allclose(t.to_numpy(getattr(fused_result, name)),
-                                       t.to_numpy(getattr(serial_result, name)), atol=1e-9, rtol=1e-9)
-    for i, opt in enumerate(options):
-        opt.evaluation_options.deduplicate_octree_corners = i != 1
-
-    if backend == 'numpy':
-        from pykeops.numpy import LazyTensor
-    else:
-        from pykeops.torch import LazyTensor
-    reductions = []
-    lazy_sum = LazyTensor.sum
-
-    def record_sum(self, *args, **kwargs):
-        if kwargs.get('axis') == 0:
-            reductions.append(kwargs)
-        return lazy_sum(self, *args, **kwargs)
 
     with patch.object(symbolic, 'symbolic_evaluator_optimized_stacked',
                       wraps=symbolic.symbolic_evaluator_optimized_stacked) as fused, \
-            patch.object(LazyTensor, 'sum', new=record_sum), \
             patch.object(_stack_ops, '_evaluate', side_effect=AssertionError('Must stay fused')):
         originals, actual = _stack_ops._evaluate_optimized(**args)
     assert fused.call_count == 1
-    assert len(reductions) == 1
-    assert reductions[0]['backend'] == 'CPU'
-    assert reductions[0]['ranges'] is not None
     assert fused.call_args.kwargs['options_list'] is options
     reduced = fused.call_args.kwargs['eval_inputs']
     assert len({len(e.xyz_to_interpolate) for e in reduced}) == 3
-    for i, (original, view, result, reference) in enumerate(zip(originals, reduced, actual, expected)):
-        eligible = i == 0 or (i == 2 and case == 'faults')
-        assert (len(view.xyz_to_interpolate) < len(original.xyz_to_interpolate)) == eligible
-        assert (view is not original) == eligible
+    for i, (original, view, result) in enumerate(zip(originals, reduced, actual)):
+        assert len(view.xyz_to_interpolate) < len(original.xyz_to_interpolate)
+        assert view is not original
         assert original.solver_input is solvers[i]
         assert fused.call_args.kwargs['weights_list'][i] is solvers[i].weights_x0
         assert view.sp_internal is original.sp_internal
         assert view.ori_internal is original.ori_internal
-        if original.fault_internal.n_faults:
-            assert original.fault_internal is inputs[i].fault_values
-            assert original.fault_internal.fault_values_everywhere.shape[1] == len(original.xyz_to_interpolate)
-            assert view.fault_internal.fault_values_everywhere.shape[1] == len(view.xyz_to_interpolate)
-            assert (view.fault_internal is not original.fault_internal) == eligible
         np.testing.assert_array_equal(t.to_numpy(original.xyz_to_interpolate),
                                       t.to_numpy(t.concatenate((inputs[i].grid.values,
                                                                inputs[i].all_surface_points.sp_coords))))
@@ -128,29 +64,10 @@ def record_sum(self, *args, **kwargs):
         assert result.n_points_per_surface is original._n_points_per_surface
         assert result.slice_feature is original._slice_feature
         assert result.debug is solvers[i].debug
-        for name in FIELDS:
-            a, b = getattr(result, name), getattr(reference, name)
-            if b is None:
-                assert a is None
-            else:
-                np.testing.assert_allclose(t.to_numpy(a), t.to_numpy(b), atol=1e-9, rtol=1e-9)
-        np.testing.assert_allclose(t.to_numpy(result.scalar_field_at_surface_points),
-                                   t.to_numpy(reference.scalar_field_at_surface_points), atol=1e-9, rtol=1e-9)
-
-    if backend == 'PYTORCH':
-        import torch
-
-        def loss(results):
-            return sum((getattr(result, name) * (torch.arange(len(result._scalar_field),
-                                                             dtype=torch.float64) + 1)).sum()
-                       for result in results for name in fields)
-
-        old_grads = torch.autograd.grad(loss(expected), targets, retain_graph=True)
-        serial_grads = torch.autograd.grad(loss(serial), targets, retain_graph=True)
-        new_grads = torch.autograd.grad(loss(actual), targets)
-        for old, new, reference in zip(old_grads, new_grads, serial_grads):
-            torch.testing.assert_close(old, new, atol=1e-7, rtol=1e-8)
-            torch.testing.assert_close(new, reference, atol=1e-7, rtol=1e-8)
+        assert len(result._scalar_field) == len(original.xyz_to_interpolate)
+        assert result._gx_field is None
+        assert result._gy_field is None
+        assert result._gz_field is None
 
 
 @pytest.mark.parametrize('deduplicate', [False, True])

From 1adde17cb30009738bebbcec603a1efc5017670f Mon Sep 17 00:00:00 2001
From: Miguel de la Varga <miguel@terranigma-solutions.com>
Date: Mon, 21 Sep 2026 15:51:52 +0200
Subject: [PATCH 6/9] [ENH] Remove non-PyKeOps fallback logic and update octree
 backend documentation

---
 docs/octree_refinement.md                            |  5 ++---
 gempy_engine/API/interp_single/_stack_ops.py         |  4 ----
 .../test_modules/test_fused_corner_evaluation.py     | 12 ------------
 3 files changed, 2 insertions(+), 19 deletions(-)

diff --git a/docs/octree_refinement.md b/docs/octree_refinement.md
index 76be72e9..3554b804 100644
--- a/docs/octree_refinement.md
+++ b/docs/octree_refinement.md
@@ -114,9 +114,8 @@ 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. Without
-PyKeOps, flat stacks use per-stack evaluation. Existing finite-fault dispatch
-restrictions remain unchanged.
+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.
 
 The triangulation selector sorts the active voxel codes once per surface call and
diff --git a/gempy_engine/API/interp_single/_stack_ops.py b/gempy_engine/API/interp_single/_stack_ops.py
index 6c0cc45c..0e605c5b 100644
--- a/gempy_engine/API/interp_single/_stack_ops.py
+++ b/gempy_engine/API/interp_single/_stack_ops.py
@@ -242,10 +242,6 @@ def _evaluate_optimized(interpolation_inputs: list[InterpolationInput], options:
                          options_per_stack: list[InterpolationOptions] | None = None) -> tuple[list[EvaluatorInput], list[ExportedFields]]:
     from gempy_engine.modules.evaluator.symbolic_evaluator import symbolic_evaluator_optimized_stacked
 
-    if not BackendTensor.use_pykeops:
-        return _evaluate(interpolation_inputs, options, solver_inputs, stack_structure,
-                         tensor_structs, stack_indices, options_per_stack)
-
     eval_inputs: list[EvaluatorInput] = []
     for idx, global_i in enumerate(stack_indices):
         stack_structure.stack_number = global_i
diff --git a/tests/test_common/test_modules/test_fused_corner_evaluation.py b/tests/test_common/test_modules/test_fused_corner_evaluation.py
index f095d48f..7aa8a405 100644
--- a/tests/test_common/test_modules/test_fused_corner_evaluation.py
+++ b/tests/test_common/test_modules/test_fused_corner_evaluation.py
@@ -68,15 +68,3 @@ def test_fused_corner_evaluation_uses_reduced_views(backend, monkeypatch):
         assert result._gx_field is None
         assert result._gy_field is None
         assert result._gz_field is None
-
-
-@pytest.mark.parametrize('deduplicate', [False, True])
-def test_non_pykeops_dispatch(backend, deduplicate):
-    interp, options, descriptor = simple_model_interpolation_input_factory()
-    options.evaluation_options.deduplicate_octree_corners = deduplicate
-    with patch.object(_stack_ops, '_evaluate', return_value=([], [])) as evaluate, \
-            patch.object(symbolic, 'symbolic_evaluator_optimized_stacked',
-                         side_effect=AssertionError('PyKeOps is disabled')):
-        assert _stack_ops._evaluate_optimized([interp], options, [], descriptor.stack_structure,
-                                              [], [0]) == ([], [])
-    assert evaluate.call_count == 1

From c66485367e594a3e1e47b513fc4a342b92294722 Mon Sep 17 00:00:00 2001
From: Miguel de la Varga <miguel@terranigma-solutions.com>
Date: Mon, 21 Sep 2026 16:00:47 +0200
Subject: [PATCH 7/9] [FIX] Add null check for inverse before restoring corner
 fields

---
 .../API/interp_single/_interp_scalar_field.py | 19 ++++++++++---------
 gempy_engine/API/interp_single/_stack_ops.py  |  3 ++-
 2 files changed, 12 insertions(+), 10 deletions(-)

diff --git a/gempy_engine/API/interp_single/_interp_scalar_field.py b/gempy_engine/API/interp_single/_interp_scalar_field.py
index fe76e4b8..012cc04c 100644
--- a/gempy_engine/API/interp_single/_interp_scalar_field.py
+++ b/gempy_engine/API/interp_single/_interp_scalar_field.py
@@ -111,18 +111,19 @@ def _evaluate_sys_eq(eval_input: Union[SolverInput, EvaluatorInput], weights: np
     else:
         exported_fields = generic_evaluator(eval_input, weights, options)
 
-    return _restore_corner_fields(exported_fields, inverse)
+    if inverse is not None:
+        _restore_corner_fields(exported_fields, inverse)
+
+    return exported_fields
 
 
-def _restore_corner_fields(exported_fields: ExportedFields, inverse) -> ExportedFields:
+def _restore_corner_fields(exported_fields: ExportedFields, inverse) -> None:
     """Expand a reduced evaluation before attaching original grid metadata."""
-    if inverse is not None:
-        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])
-    return exported_fields
+    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):
diff --git a/gempy_engine/API/interp_single/_stack_ops.py b/gempy_engine/API/interp_single/_stack_ops.py
index 0e605c5b..4a5c5b01 100644
--- a/gempy_engine/API/interp_single/_stack_ops.py
+++ b/gempy_engine/API/interp_single/_stack_ops.py
@@ -274,7 +274,8 @@ def _evaluate_optimized(interpolation_inputs: list[InterpolationInput], options:
     )
 
     for idx, exported_fields in enumerate(exported_fields_list):
-        _restore_corner_fields(exported_fields, inverses[idx])
+        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
 

From 4d1d687ec3e7d19d7dddc13ef1ff32bbc95fe7d9 Mon Sep 17 00:00:00 2001
From: Miguel de la Varga <miguel@terranigma-solutions.com>
Date: Mon, 21 Sep 2026 16:08:48 +0200
Subject: [PATCH 8/9] [ENH] Remove optional sort-once triangulation and update
 related tests/documentation

Removed `triangulation_sort_once` evaluation option and related logic. Updated the documentation and tests to reflect this removal, clarifying the independence of quad-based triangulation from sort-once behavior. Simplified triangulation code for consistency and maintainability.
---
 docs/octree_refinement.md                     | 16 +++-------
 .../multi_scalar_dual_contouring.py           |  1 -
 .../core/data/dual_contouring_data.py         |  1 -
 .../core/data/options/evaluation_options.py   |  1 -
 .../dual_contouring/dual_contouring_v2.py     |  3 +-
 .../dual_contouring/fancy_triangulation.py    | 22 +++++--------
 .../test_modules/test_octree_optimizations.py | 31 -------------------
 .../test_modules/test_quad_triangulation.py   |  2 --
 8 files changed, 13 insertions(+), 64 deletions(-)

diff --git a/docs/octree_refinement.md b/docs/octree_refinement.md
index 3554b804..23e9a348 100644
--- a/docs/octree_refinement.md
+++ b/docs/octree_refinement.md
@@ -80,14 +80,13 @@ time/memory benchmarks are still needed before recommending a different default.
 
 ## Opt-in Evaluation and Triangulation
 
-Both optimizations are independent of the refinement mode and default to `False`:
+Corner deduplication is independent of the refinement mode and defaults to `False`:
 
 ```python
 options.evaluation_options.deduplicate_octree_corners = True
-options.evaluation_options.triangulation_sort_once = True
 ```
 
-Set either selector back to `False` to use its legacy path. Both selectors are
+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
@@ -118,12 +117,6 @@ 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.
 
-The triangulation selector sorts the active voxel codes once per surface call and
-reuses that lookup across all six edge cases. It preserves triangle order,
-neighbor filtering, and normal correction; it does not rewrite edge topology or
-change vertex generation. Its lookup is local to the surface call and safe for
-parallel surface processing.
-
 ### Unique-Edge Quads
 
 To select quad-based connectivity instead of the legacy triangle construction:
@@ -135,8 +128,9 @@ options.evaluation_options.triangulation_method = TriangulationMethod.QUADS
 ```
 
 The default is `TriangulationMethod.LEGACY`. This selector is serialized with the
-other evaluation options. `triangulation_sort_once` applies only to the legacy
-method; quad mode already uses one sorted cell lookup.
+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
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 4c86da70..3a357871 100644
--- a/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py
+++ b/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py
@@ -153,7 +153,6 @@ def dual_contouring_multi_scalar(
                 n_surfaces_to_export=n_scalar_field,
                 tree_depth=options.number_octree_levels,
                 base_number=base_number,
-                triangulation_sort_once=options.evaluation_options.triangulation_sort_once,
                 triangulation_method=options.evaluation_options.triangulation_method,
                 generated_cell_coordinates=left_right_codes
             )
diff --git a/gempy_engine/core/data/dual_contouring_data.py b/gempy_engine/core/data/dual_contouring_data.py
index 557f3001..b43f22cd 100644
--- a/gempy_engine/core/data/dual_contouring_data.py
+++ b/gempy_engine/core/data/dual_contouring_data.py
@@ -29,7 +29,6 @@ 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_sort_once: bool = False
     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.
diff --git a/gempy_engine/core/data/options/evaluation_options.py b/gempy_engine/core/data/options/evaluation_options.py
index 73dc625d..c9aad945 100644
--- a/gempy_engine/core/data/options/evaluation_options.py
+++ b/gempy_engine/core/data/options/evaluation_options.py
@@ -36,7 +36,6 @@ class EvaluationOptions:
     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_sort_once: bool = False  #: Reuse one voxel-code sort across the six edge cases per surface.
     triangulation_method: TriangulationMethod = TriangulationMethod.LEGACY
     
     mesh_extraction: bool = True
diff --git a/gempy_engine/modules/dual_contouring/dual_contouring_v2.py b/gempy_engine/modules/dual_contouring/dual_contouring_v2.py
index 42cc49ca..82131f52 100644
--- a/gempy_engine/modules/dual_contouring/dual_contouring_v2.py
+++ b/gempy_engine/modules/dual_contouring/dual_contouring_v2.py
@@ -99,8 +99,7 @@ def _compute_triangulation(dc_data_per_surface: DualContouringData,
         tree_depth=tree_depth_per_surface,
         voxel_normals=voxels_normals,
         vertex=vertex,
-        base_number=dc_data_per_surface.base_number,
-        sort_once=dc_data_per_surface.triangulation_sort_once
+        base_number=dc_data_per_surface.base_number
     )
 
     # @on
diff --git a/gempy_engine/modules/dual_contouring/fancy_triangulation.py b/gempy_engine/modules/dual_contouring/fancy_triangulation.py
index be4b7f92..638af625 100644
--- a/gempy_engine/modules/dual_contouring/fancy_triangulation.py
+++ b/gempy_engine/modules/dual_contouring/fancy_triangulation.py
@@ -20,7 +20,7 @@ def _get_pack_factors(base_x, base_y, base_z):
     return BackendTensor.tfnp.stack([by * bz, bz, BackendTensor.t.array(1)], axis=0)
 
 
-def triangulate(left_right_array, valid_edges, tree_depth: int, voxel_normals, vertex, base_number: tuple[int, int, int] | list[int], sort_once: bool = False):
+def triangulate(left_right_array, valid_edges, tree_depth: int, voxel_normals, vertex, base_number: tuple[int, int, int] | list[int]):
     # * Variables
     # Determine base_number dynamically from the data to support arbitrary grid shapes
     base_x, base_y, base_z = base_number
@@ -32,10 +32,6 @@ def triangulate(left_right_array, valid_edges, tree_depth: int, voxel_normals, v
 
     # * Consts
     voxel_code = (left_right_array * pack_factors).sum(1).reshape(-1, 1)
-    sorted_lookup = None
-    if sort_once:
-        order = BackendTensor.tfnp.argsort(voxel_code.reshape(-1))
-        sorted_lookup = (voxel_code.reshape(-1)[order], order)
     # ----------
 
     indices = []
@@ -61,8 +57,7 @@ def triangulate(left_right_array, valid_edges, tree_depth: int, voxel_normals, v
             voxel_normals=voxel_normals,
             n=n,
             base_number=base_number,
-            pack_factors=pack_factors,
-            sorted_lookup=sorted_lookup
+            pack_factors=pack_factors
         )
 
         indices.append(indices_patch)
@@ -78,7 +73,7 @@ def triangulate(left_right_array, valid_edges, tree_depth: int, voxel_normals, v
 
 def compute_triangles_for_edge(edge_vector_a, edge_vector_b, edge_vector_c,
                                left_right_array_active_edge, voxel_code, voxel_normals, n,
-                               base_number, pack_factors, sorted_lookup=None):
+                               base_number, pack_factors):
     """
     Important concepts to understand this triangulation:
     - left_right_array (n_voxels, 3-directions) contains a unique number per direction describing if it is left (even) or right (odd) and the voxel level
@@ -119,7 +114,7 @@ def compute_triangles_for_edge(edge_vector_a, edge_vector_b, edge_vector_c,
 
     x, y, z = _get_indices_via_searchsorted(
         voxel_code,
-        compressed_idx_0, compressed_idx_1, compressed_idx_2, sorted_lookup=sorted_lookup
+        compressed_idx_0, compressed_idx_1, compressed_idx_2
     )
 
     # Step 5: Calculate normals and order triangles
@@ -262,7 +257,7 @@ def _calculate_normals_and_order_triangles(x, y, z, voxel_normals, n):
     return normal
 
 
-def _get_indices_via_searchsorted(voxel_code, compressed_0, compressed_1, compressed_2, sorted_lookup=None):
+def _get_indices_via_searchsorted(voxel_code, compressed_0, compressed_1, compressed_2):
     """
     Memory-efficient replacement for broadcasting, preserving exact original indices.
     """
@@ -272,11 +267,8 @@ def _get_indices_via_searchsorted(voxel_code, compressed_0, compressed_1, compre
         return empty, empty, empty
 
     # 1. Get sorting indices to map back to original positions later
-    if sorted_lookup is None:
-        sort_indices = BackendTensor.tfnp.argsort(vc_1d)
-        sorted_vc = vc_1d[sort_indices]
-    else:
-        sorted_vc, sort_indices = sorted_lookup
+    sort_indices = BackendTensor.tfnp.argsort(vc_1d)
+    sorted_vc = vc_1d[sort_indices]
 
     # 2. Search in the correctly sorted array (O(M log N) time)
     idx_0_sorted = BackendTensor.tfnp.searchsorted(sorted_vc, compressed_0)
diff --git a/tests/test_common/test_modules/test_octree_optimizations.py b/tests/test_common/test_modules/test_octree_optimizations.py
index 842b3e0e..ef53db71 100644
--- a/tests/test_common/test_modules/test_octree_optimizations.py
+++ b/tests/test_common/test_modules/test_octree_optimizations.py
@@ -1,4 +1,3 @@
-from itertools import product
 from unittest.mock import patch
 
 import numpy as np
@@ -15,7 +14,6 @@
 from gempy_engine.API.interp_single._octree_generation import _generate_corners
 from gempy_engine.API.interp_single._interp_scalar_field import _deduplicate_corners, _evaluate_sys_eq, _solve_interpolation
 from gempy_engine.API.interp_single._interp_single_feature import input_preprocess
-from gempy_engine.modules.dual_contouring.fancy_triangulation import triangulate
 from gempy_engine.modules.octrees_topology._octree_common import _generate_next_level_centers
 from tests.fixtures.simple_models import simple_model_interpolation_input_factory
 
@@ -135,31 +133,6 @@ def loss(result):
             torch.testing.assert_close(old, new, atol=1e-7, rtol=1e-8)
 
 
-@pytest.mark.parametrize('case', ['dense', 'sparse', 'empty', 'no_edges'])
-def test_triangulation_parity_and_sort_count(backend, case):
-    t = BackendTensor.t
-    coords = np.array(list(product(range(3), range(4), range(5))), dtype=np.int64)
-    rng = np.random.default_rng(17)
-    rng.shuffle(coords)
-    if case == 'sparse':
-        coords = coords[::2]
-    if case == 'empty':
-        coords = coords[:0]
-    coords = t.array(coords, dtype='int64')
-    valid = t.array(rng.random((len(coords), 12)) > 0.3, dtype=bool)
-    if case == 'no_edges':
-        valid[:] = False
-    normals = t.array(rng.normal(size=(len(coords), 12, 3)))
-    vertices = t.array(coords, dtype='float64') + 0.5
-    with patch.object(BackendTensor.tfnp, 'argsort', wraps=BackendTensor.tfnp.argsort) as sort:
-        old = triangulate(coords, valid, 1, normals, vertices, (3, 4, 5))
-        assert sort.call_count == (0 if case == 'empty' else 6)
-        sort.reset_mock()
-        new = triangulate(coords, valid, 1, normals, vertices, (3, 4, 5), sort_once=True)
-        assert sort.call_count == 1
-    np.testing.assert_array_equal(t.to_numpy(old), t.to_numpy(new))
-
-
 @pytest.mark.parametrize('flat', [False, True])
 def test_model_parity(backend, flat, monkeypatch):
     from gempy_engine.API.model.model_api import compute_model
@@ -172,7 +145,6 @@ def test_model_parity(backend, flat, monkeypatch):
     interp, options, descriptor = simple_model_interpolation_input_factory()
     options.evaluation_options.number_octree_levels = 2
     options.evaluation_options.deduplicate_octree_corners = True
-    options.evaluation_options.triangulation_sort_once = True
     monkeypatch.setenv('GEMPY_FLAT_STACKS', str(flat))
     if flat:
         # Public flat dispatch requires PyKeOps; exercise the same stack manager
@@ -258,9 +230,6 @@ def test_public_pykeops_flat_parity(backend, monkeypatch):
 def test_selectors_serialization():
     options = InterpolationOptions.from_args(range=1., c_o=1.)
     assert not options.evaluation_options.deduplicate_octree_corners
-    assert not options.evaluation_options.triangulation_sort_once
     options.evaluation_options.deduplicate_octree_corners = True
-    options.evaluation_options.triangulation_sort_once = True
     restored = InterpolationOptions.model_validate_json(options.model_dump_json())
     assert restored.evaluation_options.deduplicate_octree_corners
-    assert restored.evaluation_options.triangulation_sort_once
diff --git a/tests/test_common/test_modules/test_quad_triangulation.py b/tests/test_common/test_modules/test_quad_triangulation.py
index c7ed79ac..be5c3ef9 100644
--- a/tests/test_common/test_modules/test_quad_triangulation.py
+++ b/tests/test_common/test_modules/test_quad_triangulation.py
@@ -179,10 +179,8 @@ def test_selector_serialization():
     options = InterpolationOptions.from_args(range=1., c_o=1.)
     assert options.evaluation_options.triangulation_method is TriangulationMethod.LEGACY
     options.evaluation_options.triangulation_method = TriangulationMethod.QUADS
-    options.evaluation_options.triangulation_sort_once = True
     restored = InterpolationOptions.model_validate_json(options.model_dump_json())
     assert restored.evaluation_options.triangulation_method is TriangulationMethod.QUADS
-    assert restored.evaluation_options.triangulation_sort_once
 
 
 def test_public_model_mesh_parity(backend, monkeypatch):

From 1f05781ac334bc63c1711747f2780d36b950b02f Mon Sep 17 00:00:00 2001
From: Miguel de la Varga <miguel@terranigma-solutions.com>
Date: Mon, 21 Sep 2026 16:17:41 +0200
Subject: [PATCH 9/9] [FIX] Optimize garbage collection logic during chunk
 evaluation

---
 gempy_engine/modules/evaluator/generic_evaluator.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/gempy_engine/modules/evaluator/generic_evaluator.py b/gempy_engine/modules/evaluator/generic_evaluator.py
index 1fb8402d..13f8768f 100644
--- a/gempy_engine/modules/evaluator/generic_evaluator.py
+++ b/gempy_engine/modules/evaluator/generic_evaluator.py
@@ -60,9 +60,9 @@ def generic_evaluator(
             if gz_field is not None:
                 gz_field[slice_array] = gz_chunk  # type: ignore
 
-    # Force garbage collection every few chunks to prevent memory buildup
-    if n_chunks:
-        gc.collect()
+        # Collect every five chunks and after the final chunk.
+        if (i + 1) % 5 == 0 or i == n_chunks - 1:
+            gc.collect()
         
     if n_chunks > 5:
         print(f"Chunking done: {n_chunks} chunks")