diff --git a/docs/octree_refinement.md b/docs/octree_refinement.md
index 4c965159..23e9a348 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,76 @@ thresholds or measurements of interpolation/reporting overhead.
 
 The default remains fast. Representative curved, multi-stack, faulted, and GPU
 time/memory benchmarks are still needed before recommending a different default.
+
+## Opt-in Evaluation and Triangulation
+
+Corner deduplication is independent of the refinement mode and defaults to `False`:
+
+```python
+options.evaluation_options.deduplicate_octree_corners = True
+```
+
+Set this selector back to `False` to evaluate the full corner layout. It is
+included in `InterpolationOptions` JSON serialization.
+
+Corner deduplication uses signed integer lattice coordinates and vectorized
+unique/inverse operations on the active NumPy or Torch device. Each unique corner
+is evaluated at its first existing physical row, not reconstructed from the extent
+origin (which can shift across refinement levels). Scalar and gradient fields are
+gathered back to the full original layout before surface-point metadata, fault
+processing, segmentation, refinement, or mesh extraction. Centers, dense/custom
+grids, sections, topography, geophysics points, and appended surface points are
+never merged with corners or with each other. No lookup or evaluated field is
+cached across calls, stacks, or levels.
+
+Fault evaluation columns are gathered with the same indices only when duplicate
+corners have identical fault values. Otherwise that evaluation uses the legacy
+path. Differentiable corner coordinates or differentiable fault-value rows also
+use the legacy path: merging independent row derivatives would change autograd.
+Gradients with respect to weights, model inputs, and appended surface points are
+preserved. Empty, non-corner, and incompatible/custom corner layouts fall back
+safely. Physical duplicates can differ by roundoff; checks allow 32 dtype epsilons
+of relative/absolute error, so output parity is numerical rather than bitwise.
+Torch requires `scatter_reduce_` support. Small grids may not benefit from the
+unique operation, gathers, and equality checks (which can synchronize a GPU).
+
+Normal and flat stacks support the selector. Fused PyKeOps evaluation compresses
+each eligible stack independently, performs one block-sparse reduction with the
+different reduced lengths, and restores each result before attaching metadata.
+Ineligible stacks retain their full rows within the same fused call. Backend
+selection and existing finite-fault dispatch restrictions remain unchanged.
+External interpolation callbacks keep their existing path and full grid layout.
+
+### Unique-Edge Quads
+
+To select quad-based connectivity instead of the legacy triangle construction:
+
+```python
+from gempy_engine.core.data.options.evaluation_options import TriangulationMethod
+
+options.evaluation_options.triangulation_method = TriangulationMethod.QUADS
+```
+
+The default is `TriangulationMethod.LEGACY`. This selector is serialized with the
+other evaluation options and is independent of corner deduplication and refinement
+mode. Legacy triangulation sorts voxel codes locally for each edge case; quad mode
+uses one sorted cell lookup.
+
+Quad mode identifies each primal edge by its lower integer endpoint and direction,
+deduplicates these identities, and finds the four incident cells. Each complete
+crossing edge produces one quad, split deterministically into two triangles for
+the existing mesh output format. Winding follows the crossing-edge gradients.
+The existing tolerant crossing rule is preserved; inconsistent crossing flags
+on shared edges raise an error rather than silently creating inconsistent faces.
+
+Incomplete quads are skipped, never emitted as partial triangles.
+`mesh.dc_data.triangulation_report` records crossing edges, complete quads, and
+missing support at physical boundaries, geological masks, and internal refinement
+boundaries. Missing-cell counts are incidences and boundary categories can overlap.
+Direct callers without pre-mask cell coordinates receive an unknown interior
+boundary classification. Counts precede subsequent overlap/fault triangle removal.
+
+This is same-level connectivity, not coarse/fine transition stitching or extent
+capping, and it does not resolve ambiguous topology or guarantee watertightness.
+NumPy and CPU Torch parity are tested; GPU behavior and end-to-end performance
+still require validation.
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..3a357871 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,9 @@ def dual_contouring_multi_scalar(
                 gradients=output_on_edges[n_scalar_field][slice_object],
                 n_surfaces_to_export=n_scalar_field,
                 tree_depth=options.number_octree_levels,
-                base_number=base_number
+                base_number=base_number,
+                triangulation_method=options.evaluation_options.triangulation_method,
+                generated_cell_coordinates=left_right_codes
             )
 
             dc_data_per_surface_all.append(dc_data_per_surface)
diff --git a/gempy_engine/API/interp_single/_interp_scalar_field.py b/gempy_engine/API/interp_single/_interp_scalar_field.py
index fc6d8961..012cc04c 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,83 @@ def _solve_interpolation_result(
     return result
 
 
-def _evaluate_sys_eq(eval_input: Union[SolverInput, EvaluatorInput], weights: np.ndarray, options: InterpolationOptions) -> ExportedFields:
+def _evaluate_sys_eq(eval_input: Union[SolverInput, EvaluatorInput], weights: np.ndarray, options: InterpolationOptions,
+                     grid: EngineGrid | None = None) -> ExportedFields:
+    inverse = None
+    if options.evaluation_options.deduplicate_octree_corners:
+        eval_input, inverse = _deduplicate_corners(eval_input, grid)
     if BackendTensor.use_pykeops:
         exported_fields = symbolic_evaluator(eval_input, weights, options)
     else:
         exported_fields = generic_evaluator(eval_input, weights, options)
 
+    if inverse is not None:
+        _restore_corner_fields(exported_fields, inverse)
+
     return exported_fields
+
+
+def _restore_corner_fields(exported_fields: ExportedFields, inverse) -> None:
+    """Expand a reduced evaluation before attaching original grid metadata."""
+    for name in ('_scalar_field', '_gx_field', '_gy_field', '_gz_field'):
+        values = getattr(exported_fields, name)
+        if values is not None:
+            index = BackendTensor.t.to_numpy(inverse) if isinstance(values, np.ndarray) else inverse
+            setattr(exported_fields, name, values[index])
+
+
+def _deduplicate_corners(eval_input: SolverInput | EvaluatorInput, grid: EngineGrid | None):
+    """Call-local evaluation view; never change grid layout or surface-point metadata."""
+    if grid is None or grid.octree_grid is None or grid.corners_grid is None:
+        return eval_input, None
+    corners = grid.corners_grid.values
+    # Separate coordinate/fault row derivatives must not be redirected to a representative.
+    if len(corners) == 0 or getattr(corners, 'requires_grad', False):
+        return eval_input, None
+    faults = eval_input.fault_internal
+    fault_values = faults.fault_values_everywhere if faults.n_faults else None
+    if getattr(fault_values, 'requires_grad', False):
+        return eval_input, None
+
+    t = BackendTensor.t
+    offsets = t.array([[x, y, z] for x in (0, 1) for y in (0, 1) for z in (0, 1)], dtype='int64')
+    coordinates = (grid.octree_grid.integer_coordinates[:, None, :] + offsets).reshape(-1, 3)
+    if len(coordinates) != len(corners):
+        return eval_input, None
+    if BackendTensor.engine_backend == gempy_engine.config.AvailableBackends.PYTORCH:
+        import torch
+        unique, inverse = torch.unique(coordinates, dim=0, return_inverse=True)
+        first = torch.full((len(unique),), len(corners), dtype=torch.int64, device=coordinates.device)
+        first.scatter_reduce_(0, inverse, torch.arange(len(corners), device=coordinates.device), reduce='amin')
+    else:
+        _, first, inverse = np.unique(coordinates, axis=0, return_index=True, return_inverse=True)
+    if len(first) == len(corners):
+        return eval_input, None
+
+    start, stop = grid.corners_grid_slice.start, grid.corners_grid_slice.stop
+    xyz = eval_input.xyz_to_interpolate
+    # Use original physical rows, not extent + lattice * spacing: refined extents
+    # may carry a different origin shift. Reject non-lattice/custom corner layouts.
+    tolerance = 32 * np.finfo(BackendTensor.dtype).eps
+    if not t.allclose(xyz[start:stop], xyz[start + first][inverse], rtol=tolerance, atol=tolerance):
+        return eval_input, None
+    if fault_values is not None:
+        if not t.all(fault_values[:, start:stop] == fault_values[:, start + first][:, inverse]):
+            return eval_input, None
+
+    before = BackendTensor.arange(start, dtype='int64')
+    after = stop + BackendTensor.arange(len(xyz) - stop, dtype='int64')
+    keep = t.concatenate((before, start + first, after))
+    restore = t.concatenate((before, start + inverse,
+                             start + len(first) + BackendTensor.arange(len(xyz) - stop, dtype='int64')))
+    reduced = copy(eval_input)
+    reduced.xyz_to_interpolate = xyz[keep]
+    if fault_values is not None:
+        reduced_faults = copy(faults)
+        reduced_faults.fault_values_everywhere = fault_values[:, keep]
+        if isinstance(reduced, EvaluatorInput):
+            reduced.solver_input = copy(reduced.solver_input)
+            reduced.solver_input.fault_internal = reduced_faults
+        else:
+            reduced.fault_internal = reduced_faults
+    return reduced, restore
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..4a5c5b01 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
@@ -225,6 +225,7 @@ def _evaluate(interpolation_inputs: list[InterpolationInput], options: Interpola
             eval_input=eval_input,
             weights=eval_input.solver_input.weights_x0,
             options=options_per_stack[idx] if options_per_stack is not None else options,
+            grid=interpolation_inputs[idx].grid,
         )
 
         exported_fields.set_structure_values_from_eval_input(eval_input)
@@ -257,15 +258,24 @@ def _evaluate_optimized(interpolation_inputs: list[InterpolationInput], options:
     weights_list = [ei.solver_input.weights_x0 for ei in eval_inputs]
 
     options_list = options_per_stack or [options] * len(stack_indices)
+    reduced_inputs = []
+    inverses = []
+    for eval_input, interpolation_input, options_i in zip(eval_inputs, interpolation_inputs, options_list):
+        reduced, inverse = (_deduplicate_corners(eval_input, interpolation_input.grid)
+                            if options_i.evaluation_options.deduplicate_octree_corners else (eval_input, None))
+        reduced_inputs.append(reduced)
+        inverses.append(inverse)
 
     # Call the stacked evaluator (single PyKeOps call with block-sparse ranges)
     exported_fields_list: list[ExportedFields] = symbolic_evaluator_optimized_stacked(
-        eval_inputs=eval_inputs,
+        eval_inputs=reduced_inputs,
         weights_list=weights_list,
         options_list=options_list
     )
 
     for idx, exported_fields in enumerate(exported_fields_list):
+        if inverses[idx] is not None:
+            _restore_corner_fields(exported_fields, inverses[idx])
         exported_fields.set_structure_values_from_eval_input(eval_inputs[idx])
         exported_fields.debug = eval_inputs[idx].solver_input.debug
 
diff --git a/gempy_engine/core/data/dual_contouring_data.py b/gempy_engine/core/data/dual_contouring_data.py
index c879cf75..b43f22cd 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:
@@ -27,6 +29,9 @@ class DualContouringData:
     extra_edge_xyz: Optional[np.ndarray] = None      # (n_valid_voxels, K, 3)
     extra_edge_normals: Optional[np.ndarray] = None   # (n_valid_voxels, K, 3)
     extra_weights: Optional[np.ndarray] = None        # (n_valid_voxels, K)
+    triangulation_method: TriangulationMethod = TriangulationMethod.LEGACY
+    generated_cell_coordinates: Optional[np.ndarray] = None  # Before geological masking.
+    triangulation_report: dict = field(default_factory=dict)  # Quad support before overlap/fault triangle removal.
 
     @property
     def valid_voxels(self):
@@ -39,4 +44,3 @@ def n_valid_edges(self):
     @property
     def n_evaluations_on_edges(self):
         return self.xyz_on_edge.shape[0]
-    
\ 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 b302004c..c9aad945 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
@@ -30,6 +35,8 @@ class EvaluationOptions:
     octree_error_threshold: float = 1.  #: Number of standard deviations to consider a voxel as candidate to refine
     octree_min_level: int = 2
     octree_refinement_mode: OctreeRefinementMode = OctreeRefinementMode.FAST
+    deduplicate_octree_corners: bool = False  #: Evaluate unique corners, then restore the original row layout.
+    triangulation_method: TriangulationMethod = TriangulationMethod.LEGACY
     
     mesh_extraction: bool = True
     mesh_extraction_masking_options: MeshExtractionMaskingOptions = MeshExtractionMaskingOptions.INTERSECT
diff --git a/gempy_engine/modules/dual_contouring/dual_contouring_v2.py b/gempy_engine/modules/dual_contouring/dual_contouring_v2.py
index ae0c6c7e..82131f52 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/gempy_engine/modules/evaluator/generic_evaluator.py b/gempy_engine/modules/evaluator/generic_evaluator.py
index c945a24f..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 (i + 1) % 5 == 0 or i == n_chunks - 1:
-        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")
diff --git a/gempy_engine/modules/evaluator/symbolic_evaluator.py b/gempy_engine/modules/evaluator/symbolic_evaluator.py
index b0af7557..c3588281 100644
--- a/gempy_engine/modules/evaluator/symbolic_evaluator.py
+++ b/gempy_engine/modules/evaluator/symbolic_evaluator.py
@@ -201,7 +201,7 @@ def _run_prep(args):
             kernel_data_list = list(executor.map(_run_prep, prep_tasks))
 
         concat_kernel_data: KernelInput = _build_stacked_kernel_data(kernel_data_list)
-        
+
         eval_kernel_scalar = create_scalar_kernel(
             concat_kernel_data,
             base_options.kernel_options,
@@ -221,7 +221,7 @@ def _run_prep(args):
             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,
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..7aa8a405
--- /dev/null
+++ b/tests/test_common/test_modules/test_fused_corner_evaluation.py
@@ -0,0 +1,70 @@
+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.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
+
+
+def test_fused_corner_evaluation_uses_reduced_views(backend, monkeypatch):
+    pytest.importorskip('pykeops')
+    monkeypatch.setattr(BackendTensor, 'use_pykeops', True)
+    t = BackendTensor.t
+    inputs, solvers, structs, options = [], [], [], []
+    for i in range(3):
+        interp, opt, descriptor = simple_model_interpolation_input_factory()
+        root = RegularGrid([10, 13, -4, -2, 20, 24], [i + 2, 2, 2])
+        corners = _generate_corners(root)
+        grid = EngineGrid(octree_grid=root, corners_grid=GenericGrid(corners),
+                          custom_grid=GenericGrid(corners[:2]))
+        interp.set_temp_grid(grid)
+        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))
+        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)
+
+    with patch.object(symbolic, 'symbolic_evaluator_optimized_stacked',
+                      wraps=symbolic.symbolic_evaluator_optimized_stacked) as fused, \
+            patch.object(_stack_ops, '_evaluate', side_effect=AssertionError('Must stay fused')):
+        originals, actual = _stack_ops._evaluate_optimized(**args)
+    assert fused.call_count == 1
+    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) 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
+        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
+        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
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
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..be5c3ef9
--- /dev/null
+++ b/tests/test_common/test_modules/test_quad_triangulation.py
@@ -0,0 +1,207 @@
+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
+    restored = InterpolationOptions.model_validate_json(options.model_dump_json())
+    assert restored.evaluation_options.triangulation_method is TriangulationMethod.QUADS
+
+
+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)