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)

{'_interpolation_options': InterpolationOptions(kernel_options={'range': 5, 'c_o': 10, 'uni_degree': 1, 'i_res': 4, 'gi_res': 2, 'number_dimensions': 3, 'kernel_function': <AvailableKernelFunctions.cubic: KernelFunction(base_function=<function cubic_function at 0x7f8833a9d3f0>, derivative_div_r=<function cubic_function_p_div_r at 0x7f8833a9d480>, second_derivative=<function cubic_function_a at 0x7f8833a9d510>, consume_sq_distance=False)>, 'kernel_solver': <Solvers.DEFAULT: 1>, 'compute_condition_number': False, 'optimizing_condition_number': False, 'condition_number': None}, evaluation_options={'_number_octree_levels': 6, '_number_octree_levels_surface': 4, 'curvature_threshold': -1, 'error_threshold': 1.0, 'min_octree_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}, temp_interpolation_values=<gempy_engine.core.data.options.temp_interpolation_values.TempInterpolationValues object at 0x7f87d3e93730>, debug=True, cache_mode=CacheMode.IN_MEMORY_CACHE, cache_model_name=, block_solutions_type=BlockSolutionType.OCTREE, sigmoid_slope=50000, debug_water_tight=False),
 'grid': Grid(_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.]]),
                                      mask_topo=array([], shape=(0, 3), dtype=bool),
                                      _transform=None),
              _custom_grid=None,
              _topography=None,
              _sections=None,
              _centered_grid=None,
              values=array([[  25.,   25., -985.],
       [  25.,   25., -955.],
       [  25.,   25., -925.],
       ...,
       [ 975.,  975., -475.],
       [ 975.,  975., -445.],
       [ 975.,  975., -415.]]),
              length=array([], dtype=float64),
              _transform=None,
              _octree_levels=-1),
 '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])},
 'meta': GeoModelMeta(name='Faults_relations',
                      creation_date=None,
                      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]],
}

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
Chunking done: 6 chunks
Chunking done: 31 chunks
Chunking done: 6 chunks
Chunking done: 31 chunks

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 0x7f87daf85210>

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
Chunking done: 7 chunks
Chunking done: 7 chunks
Chunking done: 38 chunks
Chunking done: 7 chunks
Chunking done: 7 chunks
Chunking done: 38 chunks

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 0x7f87d3ec5690>

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 0x7f87dbd8e140>

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
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

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 0x7f87dbca7430>