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 0x7f857a0a1b40>, derivative_div_r=<function cubic_function_p_div_r at 0x7f857a0a1bd0>, second_derivative=<function cubic_function_a at 0x7f857a0a1c60>, consume_sq_distance=False)>, 'kernel_solver': <Solvers.DEFAULT: 1>, 'compute_condition_number': False, 'optimizing_condition_number': False, 'condition_number': None}, _number_octree_levels=6, current_octree_level=0, block_solutions_type=BlockSolutionType.OCTREE, compute_scalar_gradient=False, mesh_extraction=True, mesh_extraction_masking_options=MeshExtractionMaskingOptions.INTERSECT, mesh_extraction_fancy=True, debug=True, debug_water_tight=False, sigmoid_slope=50000, evaluation_chunk_size=50000, _number_octree_levels_surface=4, _model_name=None),
 'grid': <gempy.core.data.grid.Grid object at 0x7f851715ef50>,
 '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]],
,
 'transform': {'_is_default_transform': False,
 'position': array([-500., -500.,  650.]),
 'rotation': array([0., 0., 0.]),
 'scale': array([0.000625, 0.000625, 0.000625])}}

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.transform.apply_anisotropy(gp.data.GlobalAnisotropy.NONE)
if False:
    gp.compute_model(geo_model)
    # %%
    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
    # %%
    gpv.plot_2d(geo_model, show_data=True)
    gpv.plot_3d(geo_model, show_data=True, kwargs_plot_structured_grid={'opacity': .2})

# TODO: Add example of offsetting just one fault

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

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)

from gempy_engine.core.data.kernel_classes.solvers import Solvers

geo_model.interpolation_options.kernel_options.kernel_solver = Solvers.SCIPY_CG
geo_model.interpolation_options.kernel_options.compute_condition_number = True
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=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,  True],
 [False, False,  True],
 [False, False, False]],

/home/leguark/gempy/gempy/core/data/geo_model.py:118: UserWarning: Both octree levels and resolution are set. The default grid for the `raw_array_solution`and plots will be the dense regular grid. To use octrees instead, set resolution to None in the regular grid.
  warnings.warn("Both octree levels and resolution are set. The default grid for the `raw_array_solution`"
Setting Backend To: AvailableBackends.numpy
/home/leguark/gempy/gempy/core/data/geo_model.py:118: UserWarning: Both octree levels and resolution are set. The default grid for the `raw_array_solution`and plots will be the dense regular grid. To use octrees instead, set resolution to None in the regular grid.
  warnings.warn("Both octree levels and resolution are set. The default grid for the `raw_array_solution`"
A size: (14, 14)
CG iterations: 13
A size: (15, 15)
CG iterations: 10
A size: (85, 85)
CG iterations: 69

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

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

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/gempy/gempy/core/data/geo_model.py:118: UserWarning: Both octree levels and resolution are set. The default grid for the `raw_array_solution`and plots will be the dense regular grid. To use octrees instead, set resolution to None in the regular grid.
  warnings.warn("Both octree levels and resolution are set. The default grid for the `raw_array_solution`"
A size: (50, 50)
CG iterations: 100
A size: (38, 38)
CG iterations: 35

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