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)

{'grid': <gempy.core.data.grid.Grid object at 0x7f08ac52d1b0>,
 '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 0x7f090b7b1cf0>, derivative_div_r=<function cubic_function_p_div_r at 0x7f090b7b1d80>, second_derivative=<function cubic_function_a at 0x7f090b7b1e10>, consume_sq_distance=False)>, 'compute_condition_number': False, 'kernel_solver': <Solvers.DEFAULT: 1>}, number_octree_levels=6, current_octree_level=0, 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, _number_octree_levels_surface=4),
 '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]],

Setting Backend To: AvailableBackends.numpy
/home/leguark/gempy/gempy/core/data/geo_model.py:164: UserWarning: You are using refinement and passing a regular grid. The resolution of the regular grid will be overwritten
  warnings.warn(
A size: (14, 14)
CG iterations: 13
A size: (15, 15)
CG iterations: 10
A size: (85, 85)
CG iterations: 68

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

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

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.FAULT,
        elements=[
Element(
        name=rock1,
        color=#ffbe00,
        is_active=True
),
Element(
        name=rock2,
        color=#728f02,
        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
A size: (14, 14)
CG iterations: 13
A size: (15, 15)
CG iterations: 10
A size: (50, 50)
CG iterations: 100
A size: (38, 38)
CG iterations: 38

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