1.5: Fault relations

Importing gempy

import gempy as gp
import gempy_viewer as gpv

# Aux imports
import numpy as np
import os

np.random.seed(1515)

We import a model from an existing folder.

data_path = os.path.abspath('../../')

geo_model: gp.data.GeoModel = gp.create_geomodel(
    project_name='Faults_relations',
    extent=[0, 1000, 0, 1000, -1000, -400],
    resolution=[20, 20, 20],
    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=data_path + "/data/input_data/tut-ch1-5/tut_ch1-5_orientations.csv",
        path_to_surface_points=data_path + "/data/input_data/tut-ch1-5/tut_ch1-5_points.csv",
    )
)

print(geo_model)

meta=GeoModelMeta(name='Faults_relations', creation_date='2025-05-21T11:32:27.213985', last_modification_date=None, owner=None) structural_frame=StructuralFrame(
        structural_groups=[
StructuralGroup(
        name=default_formation,
        structural_relation=StackRelationType.ERODE,
        elements=[
Element(
        name=fault1,
        color=#015482,
        is_active=True
),
Element(
        name=fault2,
        color=#9f0052,
        is_active=True
),
Element(
        name=rock1,
        color=#ffbe00,
        is_active=True
),
Element(
        name=rock2,
        color=#728f02,
        is_active=True
),
Element(
        name=rock3,
        color=#443988,
        is_active=True
),
Element(
        name=rock4,
        color=#ff3f20,
        is_active=True
)
]
)
],
        fault_relations=
[[False]],
 grid=Grid(values=array([[  25.,   25., -985.],
       [  25.,   25., -955.],
       [  25.,   25., -925.],
       ...,
       [ 975.,  975., -475.],
       [ 975.,  975., -445.],
       [ 975.,  975., -415.]], shape=(8000, 3)), length=array([], dtype=float64), _octree_grid=None, _dense_grid=RegularGrid(resolution=array([20, 20, 20]), extent=array([    0.,  1000.,     0.,  1000., -1000.,  -400.]), values=array([[  25.,   25., -985.],
       [  25.,   25., -955.],
       [  25.,   25., -925.],
       ...,
       [ 975.,  975., -475.],
       [ 975.,  975., -445.],
       [ 975.,  975., -415.]], shape=(8000, 3)), mask_topo=array([], shape=(0, 3), dtype=bool), _transform=None), _custom_grid=None, _topography=None, _sections=None, _centered_grid=None, _transform=None, _octree_levels=-1) geophysics_input=None input_transform={'_cached_pivot': None,
 '_is_default_transform': False,
 'position': array([-500., -500.,  650.]),
 'rotation': array([0., 0., 0.]),
 'scale': array([0.000625, 0.000625, 0.000625])} interpolation_grid=None interpolation_options=InterpolationOptions(kernel_options=KernelOptions(range=1.7, c_o=10.0, uni_degree=1, i_res=4.0, gi_res=2.0, number_dimensions=3, kernel_function=AvailableKernelFunctions.cubic, kernel_solver=Solvers.DEFAULT, compute_condition_number=False, optimizing_condition_number=False, condition_number=None), evaluation_options=EvaluationOptions(_number_octree_levels=6, _number_octree_levels_surface=4, octree_curvature_threshold=-1.0, octree_error_threshold=1.0, octree_min_level=2, mesh_extraction=True, mesh_extraction_masking_options=<MeshExtractionMaskingOptions.INTERSECT: 3>, mesh_extraction_fancy=True, evaluation_chunk_size=500000, compute_scalar_gradient=False, verbose=False), debug=True, cache_mode=<CacheMode.IN_MEMORY_CACHE: 3>, cache_model_name='Faults_relations', block_solutions_type=<BlockSolutionType.DENSE_GRID: 2>, sigmoid_slope=5000000, debug_water_tight=False, temp_interpolation_values=TempInterpolationValues(current_octree_level=0))

One fault model

Setting the structural frame

fault1: gp.data.StructuralElement = geo_model.structural_frame.get_element_by_name("fault1")
fault2: gp.data.StructuralElement = geo_model.structural_frame.get_element_by_name("fault2")

# Remove the faults from the default group
default_group: gp.data.StructuralGroup = geo_model.structural_frame.get_group_by_name("default_formation")
default_group.elements.remove(fault1)
default_group.elements.remove(fault2)

# Add a new group for the fault
gp.add_structural_group(
    model=geo_model,
    group_index=0,
    structural_group_name="fault_series_1",
    elements=[fault1],
    structural_relation=gp.data.StackRelationType.FAULT,
    fault_relations=gp.data.FaultsRelationSpecialCase.OFFSET_ALL
)

print(geo_model.structural_frame)

StructuralFrame(
        structural_groups=[
StructuralGroup(
        name=fault_series_1,
        structural_relation=StackRelationType.FAULT,
        elements=[
Element(
        name=fault1,
        color=#015482,
        is_active=True
)
]
),
StructuralGroup(
        name=default_formation,
        structural_relation=StackRelationType.ERODE,
        elements=[
Element(
        name=rock1,
        color=#ffbe00,
        is_active=True
),
Element(
        name=rock2,
        color=#728f02,
        is_active=True
),
Element(
        name=rock3,
        color=#443988,
        is_active=True
),
Element(
        name=rock4,
        color=#ff3f20,
        is_active=True
)
]
)
],
        fault_relations=
[[False,  True],
 [False, False]],

geo_model.input_transform.apply_anisotropy(gp.data.GlobalAnisotropy.NONE)
gp.compute_model(geo_model)

Setting Backend To: AvailableBackends.numpy
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))
Chunking done: 6 chunks
Chunking done: 31 chunks
Chunking done: 6 chunks
Chunking done: 31 chunks
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))

Solutions: 6 Octree Levels, 5 DualContouringMeshes

print(geo_model.solutions.raw_arrays.block_matrix[0])  # This contains the block values for the fault1
print(geo_model.solutions.raw_arrays.block_matrix[1])  # This contains the block values for the formations

[2. 2. 2. ... 1. 1. 1.]
[6. 6. 6. ... 2. 2. 2.]

gpv.plot_2d(geo_model, show_data=True)
gpv.plot_3d(geo_model, show_data=True, kwargs_plot_structured_grid={'opacity': .2})


# %5
# Graben example
# --------------

<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7ff2af877fa0>

gp.add_structural_group(
    model=geo_model,
    group_index=1,
    structural_group_name="fault_series_2",
    elements=[fault2],
    structural_relation=gp.data.StackRelationType.FAULT,
    fault_relations=gp.data.FaultsRelationSpecialCase.OFFSET_ALL
)
print(geo_model.structural_frame)

gp.compute_model(geo_model)

StructuralFrame(
        structural_groups=[
StructuralGroup(
        name=fault_series_1,
        structural_relation=StackRelationType.FAULT,
        elements=[
Element(
        name=fault1,
        color=#015482,
        is_active=True
)
]
),
StructuralGroup(
        name=fault_series_2,
        structural_relation=StackRelationType.FAULT,
        elements=[
Element(
        name=fault2,
        color=#9f0052,
        is_active=True
)
]
),
StructuralGroup(
        name=default_formation,
        structural_relation=StackRelationType.ERODE,
        elements=[
Element(
        name=rock4,
        color=#ff3f20,
        is_active=True
),
Element(
        name=rock3,
        color=#443988,
        is_active=True
),
Element(
        name=rock2,
        color=#728f02,
        is_active=True
),
Element(
        name=rock1,
        color=#ffbe00,
        is_active=True
)
]
)
],
        fault_relations=
[[False,  True,  True],
 [False, False,  True],
 [False, False, False]],

Setting Backend To: AvailableBackends.numpy
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))
Chunking done: 7 chunks
Chunking done: 7 chunks
Chunking done: 38 chunks
Chunking done: 7 chunks
Chunking done: 7 chunks
Chunking done: 38 chunks
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))

Solutions: 6 Octree Levels, 6 DualContouringMeshes

gpv.plot_2d(geo_model, show_data=True)
gpv.plot_3d(geo_model, show_data=True, image=True, kwargs_plot_structured_grid={'opacity': .2})

<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7ff2a5440670>

gpv.plot_2d(geo_model, show_scalar=True, show_lith=False, series_n=0)
gpv.plot_2d(geo_model, show_scalar=True, show_lith=False, series_n=1)
gpv.plot_2d(geo_model, show_scalar=True, show_lith=False, series_n=2)

• 
• 
• 

<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7ff2a53bcbe0>

Finite Faults

Faults relations

Let’s split the formations in two groups

gp.add_structural_group(
    model=geo_model,
    group_index=2,
    structural_group_name="series_1",
    elements=[
        geo_model.structural_frame.get_element_by_name("rock4"),
        geo_model.structural_frame.get_element_by_name("rock3")
    ],
    structural_relation=gp.data.StackRelationType.ERODE
)

default_group.elements.remove(geo_model.structural_frame.get_element_by_name("rock4"))
default_group.elements.remove(geo_model.structural_frame.get_element_by_name("rock3"))

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]
    ]
    )
)
print(geo_model.structural_frame)

StructuralFrame(
        structural_groups=[
StructuralGroup(
        name=fault_series_1,
        structural_relation=StackRelationType.FAULT,
        elements=[
Element(
        name=fault1,
        color=#015482,
        is_active=True
)
]
),
StructuralGroup(
        name=fault_series_2,
        structural_relation=StackRelationType.FAULT,
        elements=[
Element(
        name=fault2,
        color=#9f0052,
        is_active=True
)
]
),
StructuralGroup(
        name=series_1,
        structural_relation=StackRelationType.ERODE,
        elements=[
Element(
        name=rock4,
        color=#ff3f20,
        is_active=True
),
Element(
        name=rock3,
        color=#443988,
        is_active=True
)
]
),
StructuralGroup(
        name=default_formation,
        structural_relation=StackRelationType.ERODE,
        elements=[
Element(
        name=rock2,
        color=#728f02,
        is_active=True
),
Element(
        name=rock1,
        color=#ffbe00,
        is_active=True
)
]
)
],
        fault_relations=
[[False,  True,  True,  True],
 [False, False, False, False],
 [False, False, False, False],
 [False, False, False, False]],

gp.compute_model(geo_model)

Setting Backend To: AvailableBackends.numpy
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))
Chunking done: 7 chunks
Chunking done: 7 chunks
Chunking done: 23 chunks
Chunking done: 17 chunks
Chunking done: 7 chunks
Chunking done: 7 chunks
Chunking done: 22 chunks
Chunking done: 17 chunks
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))

Solutions: 6 Octree Levels, 6 DualContouringMeshes

gpv.plot_2d(geo_model, show_data=True)
gpv.plot_3d(geo_model, show_data=True, image=False, kwargs_plot_structured_grid={'opacity': .2})

<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7ff2a55e83d0>