[ENH] Refactor mesh extraction and add scalar field at interface to structural elements | | GEN-12031 (#1019)

# Description
Refactored the mesh extraction process using marching cubes algorithm to improve accuracy and efficiency. The main changes include:

1. Renamed `scalar_field` to `scalar_field_at_interface` in the `StructuralElement` class for better clarity
2. Added proper assignment of scalar field values to structural elements
3. Simplified the mesh extraction logic by using the output group's scalar field matrix and mask directly
4. Improved the marching cubes implementation with proper parameters (allow_degenerate=False, method="lewiner")
5. Added comprehensive tests to verify the mesh extraction results

Relates to #mesh-extraction-improvements

# Checklist
- [x] My code uses type hinting for function and method arguments and return values.
- [x] I have created tests which cover my code.
- [x] The test code either 1. demonstrates at least one valuable use case (e.g. integration tests) 
or 2. verifies that outputs are as expected for given inputs (e.g. unit tests).
- [x] New tests pass locally with my changes.

Author: Leguark

gempy/core/data/geo_model.py

@@ -146,7 +146,7 @@ def solutions(self) -> Solutions:
         return self._solutions
 
     @solutions.setter
-    def solutions(self, value):
+    def solutions(self, value: Solutions):
         # * This is set  from the gempy engine
 
         self._solutions = value
@@ -161,6 +161,8 @@ def solutions(self, value):
 
         # * Set solutions per element
         for e, element in enumerate(self.structural_frame.structural_elements[:-1]):  # * Ignore basement
+            element.scalar_field_at_interface = value.scalar_field_at_surface_points[e]
+            
             if self._solutions.dc_meshes is None:
                 continue
             dc_mesh = self._solutions.dc_meshes[e]

gempy/core/data/structural_element.py

@@ -31,7 +31,7 @@ class StructuralElement:
     # ? Should we extract this to a separate class?
     vertices: Optional[np.ndarray] = None  #: The vertices of the element in 3D space.
     edges: Optional[np.ndarray] = None  #: The edges of the element in 3D space.
-    scalar_field: Optional[float] = None  #: The scalar field value for the element.
+    scalar_field_at_interface: Optional[float] = None  #: The scalar field value for the element.
 
     _id: int = -1
 

gempy/modules/mesh_extranction/marching_cubes.py

@@ -34,91 +34,27 @@ def set_meshes_with_marching_cubes(model: GeoModel) -> None:
 
     output_lvl0: list[InterpOutput] = model.solutions.octrees_output[0].outputs_centers
 
-    # TODO: How to get this properly in gempy
-    # get a list of indices of the lithological groups
-    lith_group_indices = []
-    fault_group_indices = []
-    index = 0
-    for i in model.structural_frame.structural_groups:
-        if i.is_fault:
-            fault_group_indices.append(index)
-        else:
-            lith_group_indices.append(index)
-        index += 1
-
-    # extract scalar field values at surface points
-    scalar_values = model.solutions.raw_arrays.scalar_field_at_surface_points
-
-    # TODO: Here I just get my own masks, cause the gempy masks dont work as expected
-    masks = _get_masking_arrays(lith_group_indices, model, scalar_values)
-
-    # TODO: Attribute of element.scalar_field was None, changed it to scalar field value of that element
-    #  This should probably be done somewhere else and maybe renamed to scalar_field_value?
-    #  This is just the most basic solution to be clear what I did
-    _set_scalar_field_to_element(model, output_lvl0, structural_groups)
-
-    # Trying to use the exiting gempy masks
-    # masks = []
-    # masks.append(
-    #     np.ones_like(model.solutions.raw_arrays.scalar_field_matrix[0].reshape(model.grid.regular_grid.resolution),
-    #                  dtype=bool))
-    # for idx in lith_group_indices:
-    #     output_group: InterpOutput = output_lvl0[idx]
-    #     masks.append(output_group.mask_components[8:].reshape(model.grid.regular_grid.resolution))
-
-    non_fault_counter = 0
     for e, structural_group in enumerate(structural_groups):
         if e >= len(output_lvl0):
             continue
 
-        # Outdated?
-        # output_group: InterpOutput = output_lvl0[e]
-        # scalar_field_matrix = output_group.exported_fields_dense_grid.scalar_field
-
-        # Specify the correct scalar field, can be removed in the future
-        scalar_field = model.solutions.raw_arrays.scalar_field_matrix[e].reshape(model.grid.regular_grid.resolution)
-
-        # pick mask depending on whether the structural group is a fault or not
-        if structural_group.is_fault:
-            mask = np.ones_like(scalar_field, dtype=bool)
+        output_group: InterpOutput = output_lvl0[e]
+        scalar_field_matrix = output_group.exported_fields_dense_grid.scalar_field
+        if structural_group.is_fault is False:
+            slice_: slice = output_group.grid.dense_grid_slice
+            mask = output_group.combined_scalar_field.squeezed_mask_array[slice_]
         else:
-            mask = masks[non_fault_counter]  # TODO: I need the entry without faults here
-            non_fault_counter += 1
+            mask = np.ones_like(scalar_field_matrix, dtype=bool)
 
         for element in structural_group.elements:
             extract_mesh_for_element(
                 structural_element=element,
                 regular_grid=regular_grid,
-                scalar_field=scalar_field,
+                scalar_field=scalar_field_matrix,
                 mask=mask
             )
 
 
-# TODO: This should be set somewhere else
-def _set_scalar_field_to_element(model, output_lvl0, structural_groups):
-    counter = 0
-    for e, structural_group in enumerate(structural_groups):
-        if e >= len(output_lvl0):
-            continue
-
-        for element in structural_group.elements:
-            element.scalar_field = model.solutions.scalar_field_at_surface_points[counter]
-            counter += 1
-
-
-# TODO: This should be set somewhere else
-def _get_masking_arrays(lith_group_indices, model, scalar_values):
-    masks = []
-    masks.append(np.ones_like(model.solutions.raw_arrays.scalar_field_matrix[0].reshape(model.grid.regular_grid.resolution),
-                              dtype=bool))
-    for idx in lith_group_indices:
-        mask = model.solutions.raw_arrays.scalar_field_matrix[idx].reshape(model.grid.regular_grid.resolution) <= \
-               scalar_values[idx][-1]
-
-        masks.append(mask)
-    return masks
-
-
 def extract_mesh_for_element(structural_element: StructuralElement,
                              regular_grid: RegularGrid,
                              scalar_field: np.ndarray,
@@ -138,10 +74,12 @@ def extract_mesh_for_element(structural_element: StructuralElement,
     """
     # Extract mesh using marching cubes
     verts, faces, _, _ = measure.marching_cubes(
-        volume=scalar_field,
-        level=structural_element.scalar_field,
+        volume=scalar_field.reshape(regular_grid.resolution),
+        level=structural_element.scalar_field_at_interface,
         spacing=(regular_grid.dx, regular_grid.dy, regular_grid.dz),
-        mask=mask
+        mask=mask.reshape(regular_grid.resolution) if mask is not None else None,
+        allow_degenerate=False,
+        method="lewiner"
     )
 
     # Adjust vertices to correct coordinates in the model's extent

test/test_modules/test_marching_cubes.py

@@ -38,12 +38,20 @@ def test_marching_cubes_implementation():
     # assert arrays.scalar_field_matrix.shape == (3, 8_000)  # * 3 surfaces, 8000 points
 
     marching_cubes.set_meshes_with_marching_cubes(model)
+    
+    # Assert
+    assert model.solutions.block_solution_type == RawArraysSolution.BlockSolutionType.DENSE_GRID
+    assert model.solutions.dc_meshes is  None
+    assert model.structural_frame.structural_groups[0].elements[0].vertices.shape == (600, 3)
+    assert model.structural_frame.structural_groups[1].elements[0].vertices.shape == (860, 3)
+    assert model.structural_frame.structural_groups[2].elements[0].vertices.shape == (1_256, 3)
+    assert model.structural_frame.structural_groups[2].elements[1].vertices.shape == (1_680, 3)
 
     if PLOT:
         gpv = require_gempy_viewer()
         gtv: gpv.GemPyToVista = gpv.plot_3d(
             model=model,
             show_data=True,
-            image=False,
+            image=True,
             show=True
         )