[ENH] Add fault relation example model and fix fault relation implementation

examples/data/input_data/tests/fault_relations_test_ori.csv

@@ -0,0 +1,13 @@
+X,Y,Z,G_x,G_y,G_z,id,nugget,formation
+350.0,500.0,-700.0,0.8949343616020251,5.4798925069145e-17,0.4461978131098087,26244822,0.01,fault1
+550.0,500.0,-600.0,-0.8949343616015264,-1.6439677520734337e-16,0.44619781311080886,124053343,0.01,fault2
+400.0,500.0,-625.0,0.0,0.0,1.0,563652887,0.01,rock4
+150.0,500.0,-550.0,0.0,0.0,1.0,563652887,0.01,rock4
+850.0,500.0,-525.0,0.0,0.0,1.0,563652887,0.01,rock4
+400.0,500.0,-675.0,0.0,0.0,1.0,438217026,0.01,rock3
+150.0,500.0,-600.0,0.0,0.0,1.0,438217026,0.01,rock3
+850.0,500.0,-575.0,0.0,0.0,1.0,438217026,0.01,rock3
+850.0,500.0,-625.0,0.0,0.0,1.0,317776925,0.01,rock2
+150.0,500.0,-650.0,0.0,0.0,1.0,317776925,0.01,rock2
+150.0,500.0,-700.0,0.0,0.0,1.0,267239155,0.01,rock1
+850.0,500.0,-675.0,0.0,0.0,1.0,267239155,0.01,rock1

examples/data/input_data/tests/fault_relations_test_surf.csv

@@ -0,0 +1,73 @@
+X,Y,Z,id,nugget,formation
+300.0,200.0,-600.0,26244822,2e-05,fault1
+300.0,500.0,-600.0,26244822,2e-05,fault1
+300.0,800.0,-600.0,26244822,2e-05,fault1
+350.0,200.0,-700.0,26244822,2e-05,fault1
+350.0,500.0,-700.0,26244822,2e-05,fault1
+350.0,800.0,-700.0,26244822,2e-05,fault1
+400.0,800.0,-800.0,26244822,2e-05,fault1
+400.0,500.0,-800.0,26244822,2e-05,fault1
+400.0,200.0,-800.0,26244822,2e-05,fault1
+500.0,500.0,-700.0,124053343,2e-05,fault2
+500.0,200.0,-700.0,124053343,2e-05,fault2
+500.0,800.0,-700.0,124053343,2e-05,fault2
+550.0,200.0,-600.0,124053343,2e-05,fault2
+550.0,500.0,-600.0,124053343,2e-05,fault2
+550.0,800.0,-600.0,124053343,2e-05,fault2
+600.0,200.0,-500.0,124053343,2e-05,fault2
+600.0,500.0,-500.0,124053343,2e-05,fault2
+600.0,800.0,-500.0,124053343,2e-05,fault2
+400.0,200.0,-625.0,563652887,2e-05,rock4
+100.0,500.0,-550.0,563652887,2e-05,rock4
+100.0,800.0,-550.0,563652887,2e-05,rock4
+100.0,200.0,-550.0,563652887,2e-05,rock4
+200.0,200.0,-550.0,563652887,2e-05,rock4
+200.0,500.0,-550.0,563652887,2e-05,rock4
+200.0,800.0,-550.0,563652887,2e-05,rock4
+800.0,200.0,-525.0,563652887,2e-05,rock4
+800.0,500.0,-525.0,563652887,2e-05,rock4
+800.0,800.0,-525.0,563652887,2e-05,rock4
+900.0,200.0,-525.0,563652887,2e-05,rock4
+900.0,500.0,-525.0,563652887,2e-05,rock4
+900.0,800.0,-525.0,563652887,2e-05,rock4
+400.0,500.0,-625.0,563652887,2e-05,rock4
+400.0,800.0,-625.0,563652887,2e-05,rock4
+400.0,800.0,-675.0,438217026,2e-05,rock3
+100.0,200.0,-600.0,438217026,2e-05,rock3
+100.0,500.0,-600.0,438217026,2e-05,rock3
+100.0,800.0,-600.0,438217026,2e-05,rock3
+200.0,200.0,-600.0,438217026,2e-05,rock3
+200.0,500.0,-600.0,438217026,2e-05,rock3
+200.0,800.0,-600.0,438217026,2e-05,rock3
+800.0,200.0,-575.0,438217026,2e-05,rock3
+800.0,500.0,-575.0,438217026,2e-05,rock3
+800.0,800.0,-575.0,438217026,2e-05,rock3
+900.0,500.0,-575.0,438217026,2e-05,rock3
+900.0,200.0,-575.0,438217026,2e-05,rock3
+900.0,800.0,-575.0,438217026,2e-05,rock3
+400.0,200.0,-675.0,438217026,2e-05,rock3
+400.0,500.0,-675.0,438217026,2e-05,rock3
+200.0,800.0,-650.0,317776925,2e-05,rock2
+200.0,500.0,-650.0,317776925,2e-05,rock2
+200.0,200.0,-650.0,317776925,2e-05,rock2
+100.0,500.0,-650.0,317776925,2e-05,rock2
+100.0,800.0,-650.0,317776925,2e-05,rock2
+100.0,200.0,-650.0,317776925,2e-05,rock2
+900.0,200.0,-625.0,317776925,2e-05,rock2
+900.0,500.0,-625.0,317776925,2e-05,rock2
+900.0,800.0,-625.0,317776925,2e-05,rock2
+800.0,200.0,-625.0,317776925,2e-05,rock2
+800.0,500.0,-625.0,317776925,2e-05,rock2
+800.0,800.0,-625.0,317776925,2e-05,rock2
+100.0,200.0,-700.0,267239155,2e-05,rock1
+100.0,500.0,-700.0,267239155,2e-05,rock1
+100.0,800.0,-700.0,267239155,2e-05,rock1
+200.0,200.0,-700.0,267239155,2e-05,rock1
+200.0,500.0,-700.0,267239155,2e-05,rock1
+200.0,800.0,-700.0,267239155,2e-05,rock1
+900.0,200.0,-675.0,267239155,2e-05,rock1
+900.0,500.0,-675.0,267239155,2e-05,rock1
+900.0,800.0,-675.0,267239155,2e-05,rock1
+800.0,200.0,-675.0,267239155,2e-05,rock1
+800.0,500.0,-675.0,267239155,2e-05,rock1
+800.0,800.0,-675.0,267239155,2e-05,rock1

gempy/API/examples_generator.py

@@ -7,8 +7,6 @@
 from gempy.core.data.enumerators import ExampleModel
 
 
-
-
 def generate_example_model(example_model: ExampleModel, compute_model: bool = True) -> gp.data.GeoModel:
     match example_model:
         case ExampleModel.HORIZONTAL_STRAT:
@@ -27,6 +25,8 @@ def generate_example_model(example_model: ExampleModel, compute_model: bool = Tr
             return _generate_graben_model(compute_model)
         case ExampleModel.GREENSTONE:
             return _generate_greenstone_model(compute_model)
+        case ExampleModel.FAULT_RELATION:
+            return _generate_fault_relation_model(compute_model)
         case _:
             raise NotImplementedError(f"Example model {example_model} not implemented.")
 
@@ -296,7 +296,7 @@ def _generate_one_fault_model_gravity(compute_model):
     from gempy_engine.core.backend_tensor import BackendTensor
     from gempy.optional_dependencies import require_pandas
     pd = require_pandas()
-    
+
     resolution = [150, 10, 150]
     extent = [0, 200, -100, 100, -100, 0]
 
@@ -394,7 +394,7 @@ def _generate_one_fault_model_gravity(compute_model):
     ## Compute model
     # %%
     geo_model.update_transform(gp.data.GlobalAnisotropy.NONE)
-    
+
     interesting_columns = pd.DataFrame()
     # x_vals = np.arange(20, 191, 10)
 
@@ -411,7 +411,7 @@ def _generate_one_fault_model_gravity(compute_model):
     gp.set_centered_grid(
         grid=geo_model.grid,
         centers=device_location,
-        resolution=np.array([75/3, 5, 150/3]),
+        resolution=np.array([75 / 3, 5, 150 / 3]),
         radius=np.array([150, 10, 300])
     )
 
@@ -443,7 +443,7 @@ def _generate_one_fault_model_gravity(compute_model):
         )
         grav = - sol.gravity
         grav[0].backward()
-    
+
     return geo_model
 
 
@@ -458,15 +458,15 @@ def _generate_graben_model(compute_model: bool) -> gp.data.GeoModel:
 
     data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/'
     path_to_data = data_path + "/data/input_data/lisa_models/"
-    
+
     geo_data: gp.data.GeoModel = gp.create_geomodel(
         project_name="Graben",
         extent=[0, 2000, 0, 2000, 0, 1600],
         resolution=[50, 50, 50],
         refinement=6,  # * For this model is better not to use octrees because we want to see what is happening in the scalar fields
         importer_helper=gp.data.ImporterHelper(
-            path_to_orientations= path_to_data + "foliations7.csv",
-            path_to_surface_points= path_to_data + "interfaces7.csv"
+            path_to_orientations=path_to_data + "foliations7.csv",
+            path_to_surface_points=path_to_data + "interfaces7.csv"
         )
     )
 
@@ -495,15 +495,66 @@ def _generate_greenstone_model(compute_model: bool) -> gp.data.GeoModel:
 
     from gempy.modules.serialization.save_load import _load_model_from_bytes
     geo_model: gp.data.GeoModel = _load_model_from_bytes(binary_file)
-    
+
     if compute_model:
         sol = gp.compute_model(
             gempy_model=geo_model,
             engine_config=gp.data.GemPyEngineConfig(
                 backend=gp.data.AvailableBackends.numpy,
                 dtype='float32'
             )
-    )
+        )
     return geo_model
-
 
+
+def _generate_fault_relation_model(compute_model: bool) -> gp.data.GeoModel:
+    # Path to input data
+    test_dir = os.path.dirname(os.path.abspath(__file__))
+    data_path = os.path.join(test_dir, '..', '..', 'examples', 'data', 'input_data', 'tests')
+    # data_path = os.path.join(test_dir, "..\\..\\examples\\data\\input_data\\tests\\")
+
+    # Create instance of new geomodel
+    geo_model = gp.create_geomodel(
+        project_name='fault_reations_test',
+        extent=[0, 1000, 0, 1000, -1000, -400],
+        resolution=[20, 20, 20],
+        refinement=4,
+        importer_helper=gp.data.ImporterHelper(
+            path_to_orientations=data_path + "/fault_relations_test_ori.csv",
+            path_to_surface_points=data_path + "/fault_relations_test_surf.csv"
+        )
+    )
+    # Map geological series to surfaces
+    gp.map_stack_to_surfaces(
+        gempy_model=geo_model,
+        mapping_object={
+                "fault_series_1": ('fault1'),
+                "fault_series_2": ('fault2'),
+                "series_1"      : ('rock4', 'rock3'),
+                "default_series": ('rock2', 'rock1')
+        }
+    )
+
+    gp.set_is_fault(geo_model, ["fault_series_1", "fault_series_2"])
+
+    gp.set_fault_relation(
+        frame=geo_model.structural_frame,
+        rel_matrix=np.array([
+                [0, 1, 1, 1],
+                [0, 0, 0, 1],
+                [0, 0, 0, 0],
+                [0, 0, 0, 0]
+        ]
+        )
+    )
+    geo_model.input_transform.apply_anisotropy(gp.data.GlobalAnisotropy.NONE)
+    if compute_model:
+        sol = gp.compute_model(
+            gempy_model=geo_model,
+            engine_config=gp.data.GemPyEngineConfig(
+                backend=gp.data.AvailableBackends.numpy,
+                dtype='float32'
+            ),
+            validate_serialization=False
+        )
+    return geo_model

gempy/core/data/_data_points_helpers.py

@@ -1,4 +1,5 @@
 import hashlib
+import sys
 
 from typing import Sequence
 
@@ -17,6 +18,8 @@ def generate_ids_from_names(name_id_map, names, x):
         ids = np.array([name_id_map[names]] * len(x))
     elif isinstance(names, Sequence) or isinstance(names, np.ndarray):
         ids = np.array([name_id_map[name] for name in names])
+    elif type(names).__module__.startswith('pandas') and type(names).__name__ == 'StringArray':
+        ids = np.array([name_id_map[name] for name in names])
     else:
         raise TypeError(f"Names should be a string or a NumPy array, not {type(names)}")
     return ids, name_id_map

gempy/core/data/enumerators.py

@@ -10,3 +10,4 @@ class ExampleModel(Enum):
     ONE_FAULT_GRAVITY = auto()
     GRABEN = auto()
     GREENSTONE = auto()
+    FAULT_RELATION = auto()

gempy/core/data/structural_frame.py

@@ -265,6 +265,7 @@ def fault_relations(self) -> np.ndarray:
                         j = self.structural_groups.index(fault_group)
                         if j <= i:  # Only consider groups that are 
                             raise ValueError(f"Fault {group.name} cannot affect older fault {fault_group.name}")
+                        fault_relations[i, j] = True
                 case (StackRelationType.FAULT, _):
                     raise ValueError(f"Fault {group.name} has an invalid fault relation")
                 case _:

gempy/modules/serialization/save_load.py

@@ -133,6 +133,7 @@ def load_model(path: str) -> GeoModel:
 
     return _load_model_from_bytes(binary_file)
 
+
 def model_to_bytes(model: GeoModel) -> bytes:
     # 1) Make a fully deterministic JSON header
     # header_dict = model.model_dump(by_alias=True)

requirements/base-requirements.txt

@@ -2,4 +2,4 @@
 
 # This install matplotlib
 gempy_viewer~=2025.1.4
-pandas
+pandas>=2.2.0,<3.0.0

test/test_modules/test_faults/test_fault_relations.py

@@ -0,0 +1,30 @@
+import numpy as np
+import gempy as gp
+from gempy.core.data.enumerators import ExampleModel
+from gempy.optional_dependencies import require_gempy_viewer
+
+PLOT = True
+
+
+def test_fault_relations_implementation():
+    # TODO! (Miguel Dec25) These fault description are not serializing!
+    model = gp.generate_example_model(ExampleModel.FAULT_RELATION, compute_model=True)
+
+    correct_relations = np.array([
+        [0, 1, 1, 1],
+        [0, 0, 0, 1],
+        [0, 0, 0, 0],
+        [0, 0, 0, 0]
+    ], dtype=bool)
+
+    # Assert
+    assert np.array_equal(model.structural_frame.fault_relations, correct_relations) == True
+
+    if PLOT:
+        gpv = require_gempy_viewer()
+        gtv: gpv.GemPyToVista = gpv.plot_3d(
+            model=model,
+            show_data=True,
+            image=True,
+            show=True
+        )

test/test_modules/_geophysics_TO_UPDATE/test_gravity.test_gravity.verify/2-layers.approved.txt → test/test_modules/test_geophysics/test_gravity.test_gravity.verify/2-layers.approved.txt

@@ -112,12 +112,15 @@
         "active_grids": 1058
     },
     "geophysics_input": {
-        "tz": [],
-        "densities": [
-            2.6,
-            2.4,
-            3.2
-        ]
+        "gravity_input": {
+            "tz": [],
+            "densities": [
+                2.6,
+                2.4,
+                3.2
+            ]
+        },
+        "magnetics_input": null
     },
     "input_transform": {
         "position": [

test/test_modules/test_pile/test_stratigraphic_pile.py

@@ -1,4 +1,3 @@
-import numpy as np
 import os
 import pandas as pd
 

test/verify_helper.py

@@ -1,3 +1,4 @@
+import os
 from typing import Literal
 
 import subprocess
@@ -13,6 +14,7 @@
 from gempy.core.data import GeoModel
 from gempy.modules.serialization.save_load import _load_model_from_bytes, model_to_bytes
 
+
 class WSLWindowsDiffReporter(GenericDiffReporter):
     def get_command(self, received, approved):
         # Convert WSL paths to Windows paths
@@ -22,20 +24,27 @@ def get_command(self, received, approved):
         cmd = [self.path] + self.extra_args + [win_received, win_approved]
         return cmd
 
-def verify_json(item, name: str):
 
+class LinuxDiffReporter(GenericDiffReporter):
+    pass
+
+
+def verify_json(item, name: str):
     config = GenericDiffReporterConfig(
         name="custom",
         path=r"pycharm",
-        extra_args= ["diff"]
+        extra_args=["diff"]
     )
 
+    if os.environ.get("WSL_DISTRO_NAME"):
+        reporter = WSLWindowsDiffReporter(config)
+    else:
+        reporter = LinuxDiffReporter(config)
+
     parameters: Options = NamerFactory \
         .with_parameters(name) \
-        .with_reporter(
-        reporter=(WSLWindowsDiffReporter(config))
-    )
-
+        .with_reporter(reporter)
+
     verify(item, options=parameters)
 
 
@@ -91,8 +100,10 @@ def compare(self, received_path: str, approved_path: str) -> bool:
         except BaseException:
             return False
 
+
 class JsonSerializer:
     """Serializer that writes JSON with an indent and declares its own extension."""
+
     def get_default_extension(self) -> str:
         return "json"