[ENH] Improve fault handling in dual contouring with arbitrary grid shapes and conservative triangle removal

gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py

@@ -62,7 +62,7 @@ def dual_contouring_multi_scalar(
         return all_meshes
 
     # * 1) Triangulation code
-    left_right_codes = get_triangulation_codes(octree_list)
+    left_right_codes, base_number = get_triangulation_codes(octree_list)
 
     # * 2) Dual contouring mask
     # ? I guess this mask is different that erosion mask
@@ -110,7 +110,7 @@ def dual_contouring_multi_scalar(
     # Generate meshes for each scalar field
     if LEGACY:=False:
         for n_scalar_field in range(data_descriptor.stack_structure.n_stacks):
-            _compute_meshes_legacy(all_left_right_codes, all_mask_arrays, all_meshes, all_surfaces_intersection, all_valid_edges, n_scalar_field, octree_leaves, options, output_on_edges)
+            _compute_meshes_legacy(all_left_right_codes, all_mask_arrays, all_meshes, all_surfaces_intersection, all_valid_edges, n_scalar_field, octree_leaves, options, output_on_edges, base_number)
     else:
         dc_data_per_surface_all = []
         for n_scalar_field in range(data_descriptor.stack_structure.n_stacks):
@@ -133,7 +133,8 @@ def dual_contouring_multi_scalar(
                     left_right_codes=all_left_right_codes[n_scalar_field],
                     gradients=output_on_edges[n_scalar_field][slice_object],
                     n_surfaces_to_export=n_scalar_field,
-                    tree_depth=options.number_octree_levels
+                    tree_depth=options.number_octree_levels,
+                    base_number=base_number
                 )
 
                 dc_data_per_surface_all.append(dc_data_per_surface)
@@ -153,7 +154,8 @@ def dual_contouring_multi_scalar(
 def _compute_meshes_legacy(all_left_right_codes: list[Any], all_mask_arrays: np.ndarray,
                            all_meshes: list[DualContouringMesh], all_stack_intersection: list[Any],
                            all_valid_edges: list[Any], n_scalar_field: int,
-                           octree_leaves: OctreeLevel, options: InterpolationOptions, output_on_edges: list[np.ndarray]):
+                           octree_leaves: OctreeLevel, options: InterpolationOptions, output_on_edges: list[np.ndarray],
+                           base_number: tuple[int, int, int]):
     output: InterpOutput = octree_leaves.outputs_centers[n_scalar_field]
     mask = all_mask_arrays[n_scalar_field]
 
@@ -168,6 +170,8 @@ def _compute_meshes_legacy(all_left_right_codes: list[Any], all_mask_arrays: np.
         gradients=output_on_edges[n_scalar_field],
         n_surfaces_to_export=output.scalar_field_at_sp.shape[0],
         tree_depth=options.number_octree_levels,
+        base_number=base_number,
+        left_right_codes=all_left_right_codes[n_scalar_field]
     )
 
     meshes: List[DualContouringMesh] = compute_dual_contouring(

gempy_engine/core/data/dual_contouring_data.py

@@ -16,6 +16,7 @@ class DualContouringData:
     gradients: np.ndarray = None
 
     tree_depth: int = -1
+    base_number: tuple[int, int, int] = None
     # Water tight 
 
     bias_center_mass: np.ndarray = None  # * Only for testing

gempy_engine/modules/dual_contouring/_dual_contouring.py

@@ -90,6 +90,8 @@ def _parallel_process_surfaces(dc_data_per_stack, left_right_codes, debug, num_w
             'gradients': dc_data_per_stack.gradients,
             'n_surfaces_to_export'    : dc_data_per_stack.n_surfaces_to_export,
             'tree_depth'              : dc_data_per_stack.tree_depth,
+            'base_number'             : dc_data_per_stack.base_number,
+            'left_right_codes'        : dc_data_per_stack.left_right_codes
     }
 
     # Create surface index chunks

gempy_engine/modules/dual_contouring/_dual_contouring_v2.py

@@ -81,7 +81,8 @@ def _parallel_process_surfaces_v2(dc_data_list: list[DualContouringData], num_wo
                 'gradients'           : dc_data.gradients,
                 'left_right_codes'    : dc_data.left_right_codes,
                 'n_surfaces_to_export': dc_data.n_surfaces_to_export,
-                'tree_depth'          : dc_data.tree_depth
+                'tree_depth'          : dc_data.tree_depth,
+                'base_number'         : dc_data.base_number
         }
         dc_data_dicts.append(dc_data_dict)
 
@@ -132,7 +133,8 @@ def _process_surface_batch_v2(surface_indices, dc_data_dicts, left_right_codes):
             gradients=dc_data_dict['gradients'],
             left_right_codes=dc_data_dict['left_right_codes'],
             n_surfaces_to_export=dc_data_dict['n_surfaces_to_export'],
-            tree_depth=dc_data_dict['tree_depth']
+            tree_depth=dc_data_dict['tree_depth'],
+            base_number=dc_data_dict['base_number']
         )
         # Process the surface
         mesh = _process_one_surface(dc_data, dc_data.left_right_codes)
@@ -155,22 +157,5 @@ def _process_one_surface(dc_data: DualContouringData, left_right_codes) -> DualC
     )
 
     vertices_numpy = BackendTensor.t.to_numpy(vertices)
-    if TRIMESH_LAST_PASS := True:
-        vertices_numpy, indices_numpy = _last_pass(vertices_numpy, indices_numpy)
-
     mesh = DualContouringMesh(vertices_numpy, indices_numpy, dc_data)
     return mesh
-
-
-def _last_pass(vertices, indices):
-    """Apply trimesh post-processing if available."""
-    if not TRIMESH_AVAILABLE:
-        return vertices, indices
-
-    try:
-        mesh = trimesh.Trimesh(vertices=vertices, faces=indices)
-        mesh.fill_holes()
-        return mesh.vertices, mesh.faces
-    except Exception as e:
-        print(f"Warning: Trimesh post-processing failed: {e}")
-        return vertices, indices

gempy_engine/modules/dual_contouring/_gen_vertices.py

@@ -24,7 +24,9 @@ def _compute_vertices(dc_data_per_stack: DualContouringData,
         dxdydz=dc_data_per_stack.dxdydz,
         gradients=dc_data_per_stack.gradients,
         n_surfaces_to_export=dc_data_per_stack.n_surfaces_to_export,
-        tree_depth=dc_data_per_stack.tree_depth
+        tree_depth=dc_data_per_stack.tree_depth,
+        base_number=dc_data_per_stack.base_number,
+        left_right_codes=dc_data_per_stack.left_right_codes
     )
 
     vertices_numpy = generate_dual_contouring_vertices(dc_data_per_surface, slice_object, debug)

gempy_engine/modules/dual_contouring/_parallel_triangulation.py

@@ -98,10 +98,13 @@ def _process_single_surface(i, dc_data_per_stack, valid_edges_per_surface, left_
             dxdydz=dc_data_per_stack.dxdydz,
             gradients=dc_data_per_stack.gradients,
             n_surfaces_to_export=dc_data_per_stack.n_surfaces_to_export,
-            tree_depth=dc_data_per_stack.tree_depth
+            tree_depth=dc_data_per_stack.tree_depth,
+            base_number=dc_data_per_stack.base_number,
+            left_right_codes=dc_data_per_stack.left_right_codes
         )
 
-        if left_right_codes is None:
+        from ._dual_contouring import triangulate_dual_contouring
+        if dc_data_per_surface.left_right_codes is None:
             # Legacy triangulation
             indices = triangulate_dual_contouring(dc_data_per_surface)
         else:
@@ -147,7 +150,8 @@ def _process_single_surface(i, dc_data_per_stack, valid_edges_per_surface, left_
                 valid_edges=valid_voxels_per_surface,
                 tree_depth=tree_depth_per_surface,
                 voxel_normals=voxel_normal,
-                vertex=vertex
+                vertex=None,
+                base_number=dc_data_per_surface.base_number
             )
             indices = BackendTensor.t.concatenate(indices, axis=0)
 

gempy_engine/modules/dual_contouring/_sequential_triangulation.py

@@ -58,7 +58,8 @@ def _compute_triangulation(dc_data_per_surface: DualContouringData,
         valid_edges=valid_voxels_per_surface,
         tree_depth=tree_depth_per_surface,
         voxel_normals=voxels_normals,
-        vertex=vertex
+        vertex=vertex,
+        base_number=dc_data_per_surface.base_number
     )
 
     # @on