[ENH] Add octree refinement modes, diagnostics, and support stencils with tests

docs/octree_refinement.md

@@ -0,0 +1,79 @@
+# Octree refinement support modes
+
+```python
+from gempy_engine.core.data import OctreeRefinementMode
+
+options.evaluation_options.number_octree_levels_surface = 4
+options.evaluation_options.octree_refinement_mode = OctreeRefinementMode.PRECISE
+```
+
+The mode also accepts the strings `"fast"`, `"balanced"`, and `"precise"`.
+
+| Mode | Support stencil | Isolated interior parent count |
+| --- | --- | --- |
+| `fast` (default) | No additional support; existing selector | 1 |
+| `balanced` | Six face neighbors | 7 |
+| `precise` | 26 face, edge, and corner neighbors | 27 |
+
+This setting does not change octree depth, minimum-level refinement, curvature
+thresholds, or categorical surface selection. One dilation of their combined
+selection is applied per transition, only while generating surface extraction
+levels and only when mesh extraction is enabled. Support cells never seed another
+dilation in the same transition. Each selected parent generates eight children.
+
+`balanced` is a cost/coverage compromise, **not 2:1 octree balancing**. It does not
+provide diagonal support around primal edges. `precise` supplies the full touching
+neighborhood among existing sparse cells, clipped to the physical domain. Neither
+mode guarantees discovery of unsampled features or watertight surfaces. A surface
+newly detected at the outer edge of an earlier support band can still request an
+absent branch; this produces a warning rather than silently claiming closure.
+Physical extent capping and extraction-mask stitching are separate concerns.
+
+## Diagnostics
+
+With `options.debug = True` (or evaluation `verbose`), generated octree grids expose
+`refinement_debug`: primary surface, additional, and support-only parent masks;
+their counts; final count; closure multiplier; generated children; and unique
+missing in-domain support requests. These masks index the **previous** generation,
+like `active_cells`. Counts for surface and additional selection may overlap.
+
+For non-fast modes, or when `options.debug` is true, each extracted mesh exposes
+`support_report`. This CPU diagnostic enumerates unique sampled sign-changing
+primal edges before extraction masking. Missing incident cells are classified as
+physical, mask-removed, or never generated. Records include edge direction and
+integer coordinate, missing coordinates, and the last existing ancestor level.
+An edge can have more than one boundary classification. Internal refinement
+failures emit a warning; tests can assert the corresponding count is zero.
+
+The report checks true sign changes, not the extractor's slightly extrapolated
+edge-intersection tolerance. It does not audit later fault/overlap triangle
+removal or all final mesh edge incidences. Debug reports and coordinate sets add
+CPU memory/time overhead, especially for many surfaces.
+
+Integer extraction coordinates now use signed int64 and theoretical domain
+bounds rather than byte packing and maximum active coordinates. Normal fast-mode
+selection is unchanged; meshes affected by the old coordinate/bounds defects can
+change even in fast mode.
+
+## Performance
+
+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.
+
+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):
+
+| Selection | Primary | Balanced final | Precise final |
+| --- | ---: | ---: | ---: |
+| Plane `ix == 16` | 1,024 | 3,072 (3x) | 3,072 (3x) |
+| Sphere band `abs(norm(coord - 15.5) - 9) < 0.5` | 1,032 | 2,696 (2.61x) | 4,024 (3.90x) |
+| Isolated seeds `all(coord % 4 == 2)` | 512 | 3,584 (7x) | 13,824 (27x) |
+
+Single-call lookup times on the development machine were 0.36–0.63 ms for
+balanced and 1.05–1.68 ms for precise; these are indicative, not performance
+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.

gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py

@@ -1,4 +1,5 @@
 import copy
+import warnings
 from typing import List, Any
 
 import numpy as np
@@ -22,6 +23,8 @@
 from ...modules.dual_contouring.dual_contouring_interface import (find_intersection_on_edge, get_triangulation_codes,
                                                                   get_masked_codes, mask_generation)
 from ...modules.dual_contouring.overlapping import average_overlapping_vertices, remove_fault_overlap_triangles
+from ...modules.dual_contouring._support_report import mesh_support_report
+from ...core.data.options.evaluation_options import OctreeRefinementMode
 
 
 @gempy_profiler_decorator
@@ -106,6 +109,7 @@ def dual_contouring_multi_scalar(
     surface_to_stack = []  # track which stack each surface belongs to
     # Generate meshes for each scalar field
     dc_data_per_surface_all = []
+    support_reports = []
     stack_relations = data_descriptor.stack_structure.masking_descriptor
     for n_scalar_field in range(data_descriptor.stack_structure.n_stacks):
         if stack_relations[n_scalar_field] is StackRelationType.NULL_SPACE:
@@ -114,6 +118,24 @@ def dual_contouring_multi_scalar(
         mask = all_mask_arrays[n_scalar_field]
         n_surfaces_to_export = output.scalar_field_at_sp.shape[0]
         for surface_i in range(n_surfaces_to_export):
+            report = None
+            if options.debug or options.evaluation_options.octree_refinement_mode != OctreeRefinementMode.FAST:
+                report = mesh_support_report(
+                    left_right_codes,
+                    output.exported_fields.scalar_field[output.grid.corners_grid_slice],
+                    output.scalar_field_at_sp[surface_i], base_number, mask,
+                    surface_index=surface_i,
+                    ancestor_coordinates=[level.grid.octree_grid.integer_coordinates for level in octree_list[:-1]]
+                )
+                report['stack_index'] = n_scalar_field
+                if report['internal_refinement_boundary_edge_count']:
+                    warnings.warn(
+                        f"Stack {n_scalar_field}, surface {surface_i}: "
+                        f"{report['internal_refinement_boundary_edge_count']} crossing edges lack "
+                        "in-domain refinement support; see mesh.support_report.",
+                        RuntimeWarning, stacklevel=2
+                    )
+            support_reports.append(report)
             valid_edges = all_valid_edges[n_scalar_field]
             valid_edges_per_surface = valid_edges.reshape((n_surfaces_to_export, -1, 12))
             slice_object = _surface_slicer(surface_i, valid_edges_per_surface)
@@ -153,6 +175,8 @@ def dual_contouring_multi_scalar(
             dc_data_list=dc_data_per_surface_all,
             max_workers=None
         )
+        for mesh, report in zip(all_meshes, support_reports):
+            mesh.support_report = report
         # endregion
         # Save differentiable vertices before in-place overlap modifications,
         # then replace mesh.vertices with a detached copy so averaging doesn't

gempy_engine/API/interp_single/interp_features.py

@@ -28,7 +28,8 @@ def interpolate_n_octree_levels(interpolation_input: InterpolationInput, options
             grid_1_centers: EngineGrid = get_next_octree_grid(
                 prev_octree=next_octree,
                 evaluation_options=options.evaluation_options,
-                current_octree_level=i
+                current_octree_level=i,
+                debug=options.debug
             )
             interpolation_input.set_temp_grid(grid_1_centers)
         octree_list.append(next_octree)

gempy_engine/core/data/__init__.py

@@ -3,4 +3,5 @@
 from .kernel_classes.orientations import Orientations, OrientationsInternals
 from .kernel_classes.surface_points import SurfacePoints, SurfacePointsInternals
 from .options.interpolation_options import InterpolationOptions
+from .options.evaluation_options import OctreeRefinementMode
 from .solutions import Solutions

gempy_engine/core/data/dual_contouring_mesh.py

@@ -10,6 +10,7 @@ class DualContouringMesh:
     vertices: np.ndarray
     edges: np.ndarray
     dc_data: Optional[DualContouringData] = None  # * In principle we need this just for testing
+    support_report: Optional[dict] = None
 
     def __repr__(self):
         return f"DualContouringMesh({self.vertices.shape[0]} vertices, {self.edges.shape[0]} edges)"

gempy_engine/core/data/options/evaluation_options.py

@@ -8,6 +8,13 @@
 from .micro_anisotropic_options import MicroAnisotropicOptions
 
 
+class OctreeRefinementMode(str, enum.Enum):
+    """Support stencil, not 2:1 balancing or a guarantee of watertightness."""
+    FAST = "fast"
+    BALANCED = "balanced"
+    PRECISE = "precise"
+
+
 class MeshExtractionMaskingOptions(enum.Enum):
     NOTHING = enum.auto()  # * This is only for testing
     DISJOINT = enum.auto()
@@ -22,6 +29,7 @@ class EvaluationOptions:
     octree_curvature_threshold: float = -1.  #: Threshold to do octree refinement due to curvature to deal with angular geometries. This curvature assumes that 1 is the maximum curvature of any voxel
     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
 
     mesh_extraction: bool = True
     mesh_extraction_masking_options: MeshExtractionMaskingOptions = MeshExtractionMaskingOptions.INTERSECT

gempy_engine/core/data/regular_grid.py

@@ -16,6 +16,8 @@ class RegularGrid:
 
     _active_cells: np.ndarray = field(default=None, repr=False, init=False)
     left_right: np.ndarray = field(default=None, repr=False, init=False)
+    _integer_coordinates: np.ndarray = field(default=None, repr=False, init=False)
+    refinement_debug: dict = field(default=None, repr=False, init=False)
 
     values: np.ndarray = field(default=None, repr=False, init=False)
     original_values: np.ndarray = field(default=None, repr=False, init=False)  #: When the regular grid is representing a octree level, only active cells are stored in values. This is the original values of the regular grid.
@@ -68,6 +70,10 @@ def from_octree_level(cls, xyz_coords_octree: np.ndarray, previous_regular_grid:
         regular_grid_for_octree_level.values = xyz_coords_octree  # ! Overwrite the common values
         regular_grid_for_octree_level._active_cells = active_cells
         regular_grid_for_octree_level.left_right = left_right
+        regular_grid_for_octree_level._integer_coordinates = (
+            2 * BackendTensor.tfnp.repeat(previous_regular_grid.integer_coordinates[active_cells], 8, axis=0)
+            + BackendTensor.t.array(left_right, dtype='int64')
+        )
 
         return regular_grid_for_octree_level
 
@@ -86,6 +92,16 @@ def from_schema(cls, schema: GridSchema):
             left_right=None
         )
 
+    @property
+    def integer_coordinates(self):
+        """Signed lattice coordinates in sparse row order (Z fastest at root)."""
+        if self._integer_coordinates is None:
+            axes = [BackendTensor.arange(int(n), dtype='int64') for n in self.regular_grid_shape]
+            self._integer_coordinates = BackendTensor.t.stack(
+                BackendTensor.t.meshgrid(*axes, indexing='ij'), axis=-1
+            ).reshape(-1, 3)
+        return self._integer_coordinates
+
     @property
     def active_cells(self) -> np.ndarray:
         if self._active_cells is not None:

gempy_engine/modules/dual_contouring/_support_report.py

@@ -0,0 +1,69 @@
+"""CPU-only topology diagnostics, independent of triangle-candidate filtering."""
+from itertools import product
+
+import numpy as np
+
+from ...core.backend_tensor import BackendTensor
+
+
+def mesh_support_report(coordinates, scalar_corners, isovalue, domain_shape,
+                        mask=None, surface_index=0, ancestor_coordinates=()):
+    """Classify missing incident cells for unique sampled sign-changing edges.
+
+    This diagnoses sampled crossings, not unsampled components or later triangle
+    removal by overlap/fault processing. Counts of missing cells are incidences.
+    """
+    to_numpy = BackendTensor.t.to_numpy
+    coords = np.asarray(to_numpy(coordinates), dtype=np.int64)
+    scalar = np.asarray(to_numpy(scalar_corners)).reshape(-1, 8)
+    iso = float(isovalue)
+    retained = np.ones(len(coords), dtype=bool) if mask is None else np.asarray(to_numpy(mask), dtype=bool)
+    generated = set(map(tuple, coords))
+    kept = set(map(tuple, coords[retained]))
+    ancestors = [set(map(tuple, to_numpy(c))) for c in ancestor_coordinates]
+    bounds = tuple(int(n) for n in domain_shape)
+    corners = np.array(list(product((0, 1), repeat=3)), dtype=np.int64)
+    edges = set()
+    for direction, pairs in enumerate((
+        ((0, 4), (1, 5), (2, 6), (3, 7)),
+        ((0, 2), (1, 3), (4, 6), (5, 7)),
+        ((0, 1), (2, 3), (4, 5), (6, 7)),
+    )):
+        for a, b in pairs:
+            crossing = (scalar[:, a] >= iso) != (scalar[:, b] >= iso)
+            edges.update((direction, *p) for p in coords[crossing] + corners[a])
+    report = dict(surface_index=surface_index, crossing_edge_count=len(edges),
+                  missing_incident_cell_count=0, physical_boundary_edge_count=0,
+                  mask_boundary_edge_count=0, internal_refinement_boundary_edge_count=0,
+                  violations=[])
+    for edge in sorted(edges):
+        direction, *origin = edge
+        transverse = [i for i in range(3) if i != direction]
+        missing = []
+        kinds = set()
+        for offsets in product((-1, 0), repeat=2):
+            cell = list(origin)
+            for axis, offset in zip(transverse, offsets):
+                cell[axis] += offset
+            cell = tuple(cell)
+            if cell in kept:
+                continue
+            outside = any(c < 0 or c >= n for c, n in zip(cell, bounds))
+            kind = 'physical' if outside else 'mask' if cell in generated else 'refinement'
+            kinds.add(kind)
+            stopped = None
+            if kind == 'refinement':
+                for level, existing in enumerate(ancestors):
+                    shift = len(ancestors) - level
+                    if tuple(c >> shift for c in cell) in existing:
+                        stopped = level
+            missing.append(dict(coordinate=cell, kind=kind, outside_extent=outside,
+                                ancestor_stop_level=stopped))
+        if missing:
+            report['missing_incident_cell_count'] += len(missing)
+            for kind, key in (('physical', 'physical_boundary_edge_count'),
+                              ('mask', 'mask_boundary_edge_count'),
+                              ('refinement', 'internal_refinement_boundary_edge_count')):
+                report[key] += kind in kinds
+            report['violations'].append(dict(direction=direction, coordinate=tuple(origin), missing=missing))
+    return report