GemPy Finite Fault

Finite-length fault modelling in GemPy using ellipsoidal implicit masks.

This repository contains reproducible Jupyter Notebook examples demonstrating how finite faults can be represented in GemPy. Instead of allowing the displacement effect of a fault to extend across the entire geological model, an additional three-dimensional implicit function is used to restrict the fault to a user-defined spatial region.

The repository includes:

A minimal single-fault example based on GemPy's built-in ONE_FAULT model.
A multi-fault relay-ramp example containing four spatially restricted faults.
Example surface-point and orientation datasets.
PyVista-based visualisation of geological surfaces and finite-fault ellipsoids.

Status: Research prototype. The implementation uses APIs from both gempy and gempy_engine and may require adaptation when used with other GemPy versions.

Background
Relationship to GemPy
Finite-Fault Concept
Repository Structure
Installation
Quick Start
Example 1: Single Finite Fault
Example 2: Relay Ramp with Four Finite Faults
Input Data Format
Modelling Workflow
Important Parameters
Fault Relations
Visualisation
Using Your Own Data
CPU and GPU Execution
Known Limitations
Troubleshooting
Contributing
License
Acknowledgements

Background

Faults in implicit geological models are commonly represented by interpolated scalar fields. Without an additional spatial constraint, the mathematical fault surface and its displacement effect may extend farther than the geologically interpreted fault.

A finite fault should only affect the geological model within a limited three-dimensional region.

This repository demonstrates how such a region can be defined using an ellipsoidal implicit function. The ellipsoid acts as a spatial mask controlling where the fault displacement is active.

Typical applications include:

Fault tips.
Isolated fault segments.
Relay ramps.
Overlapping fault systems.
Synthetic structural-geology experiments.
Sensitivity studies of fault length, height, orientation, and interaction.

Relationship to GemPy

This repository is not a fork, replacement, or independent reimplementation of GemPy.

GemPy remains responsible for:

Interpolating geological interfaces and orientations.
Constructing implicit scalar fields.
Organising formations into structural groups.
Representing fault relationships.
Computing lithological models.
Extracting geological surfaces.
Providing model results for visualisation.

This repository adds example configurations for spatially restricting faults through the following GemPy components:

gp.implicit_functions.ellipsoid_3d_factory
gp.data.Transform
gp.data.FiniteFaultData
gp.data.FaultsData

It also accesses structural relationship definitions from gempy_engine:

from gempy_engine.core.data.stack_relation_type import StackRelationType

The finite-fault configuration is attached to a fault's structural group before calling:

gp.compute_model(...)

Conceptually, the workflow is:

flowchart LR
    A[Surface points and orientations] --> B[GemPy GeoModel]
    B --> C[Structural groups]
    C --> D[Fault relations]
    D --> E[Ellipsoidal finite-fault masks]
    E --> F[GemPy model computation]
    F --> G[3D geological model]
    F --> H[PyVista visualisation]

Finite-Fault Concept

An ellipsoid is defined by:

A centre:

[ \mathbf{c} = \begin{bmatrix} c_x & c_y & c_z \end{bmatrix} ]

Three radii:

[ \mathbf{r} = \begin{bmatrix} r_x & r_y & r_z \end{bmatrix} ]

An optional transformation containing translation, rotation, and scale.

In its local coordinate system, an ideal ellipsoid can be described by:

[ \left(\frac{x-c_x}{r_x}\right)^2+ \left(\frac{y-c_y}{r_y}\right)^2+ \left(\frac{z-c_z}{r_z}\right)^2 = 1 ]

GemPy's ellipsoid implicit function provides a continuous scalar field rather than only a geometric boundary. This field is used by FiniteFaultData to control the spatial activation of the fault.

Because GemPy applies internal coordinate transformations, the centre and radius must first be transformed into the model's internal coordinate system:

scaled_center = geo_model.input_transform.apply(
    center.reshape(1, -1)
)[0]

scaled_radius = geo_model.input_transform.scale_points(
    radius.reshape(1, -1)
)[0]

The original model-coordinate values should therefore be defined first and transformed only when constructing the finite-fault data.

Repository Structure

GemPy_finite_fault/
│
├── gempy_finite_fault_test_one_fault.ipynb
├── finite_fault_Relay_ramp_gempy_function.ipynb
│
├── model5_surface_points.csv
├── model5_orientations.csv
│
├── Relay_ramp_obs_points.csv
└── Relay_ramp_obs_orientation.csv

Files

File	Description
`gempy_finite_fault_test_one_fault.ipynb`	Minimal example comparing a conventional fault configuration with a spatially restricted finite fault.
`finite_fault_Relay_ramp_gempy_function.ipynb`	Four-fault relay-ramp example with an individual finite-fault ellipsoid assigned to each fault group.
`model5_surface_points.csv`	Surface points for the single-fault synthetic model.
`model5_orientations.csv`	Orientation data for the single-fault synthetic model.
`Relay_ramp_obs_points.csv`	Surface points for the relay-ramp model.
`Relay_ramp_obs_orientation.csv`	Orientation data for the relay-ramp model.

The notebooks use relative paths. They should therefore be launched from the repository root unless the input paths are changed manually.

Installation

1. Clone the repository

git clone https://github.com/yangjiandendi/GemPy_finite_fault.git
cd GemPy_finite_fault

2. Create a dedicated environment

Python 3.10 is recommended. The existing notebooks were created using Python 3.10.

Using Conda:

conda create -n gempy-finite-fault python=3.10
conda activate gempy-finite-fault

Alternatively, using Python's built-in virtual environment:

python -m venv .venv

On Windows:

.venv\Scripts\activate

On Linux or macOS:

source .venv/bin/activate

3. Install the required packages

python -m pip install --upgrade pip
python -m pip install gempy gempy-viewer numpy pyvista jupyterlab pytest
python -m pip install torch

The main dependencies are:

Package	Purpose
`gempy`	Geological modelling API and data structures.
`gempy-engine`	Numerical engine used internally by GemPy.
`gempy-viewer`	GemPy plotting and interactive visualisation.
`numpy`	Array operations and parameter definitions.
`torch`	PyTorch computation backend and optional GPU execution.
`pyvista`	Three-dimensional mesh and ellipsoid visualisation.
`jupyterlab`	Execution of the example notebooks.
`pytest`	Imported by the single-fault test notebook.

gempy-engine is normally installed as a dependency of gempy. It generally does not need to be installed separately.

4. Verify the installation

python -c "import gempy; import gempy_viewer; import pyvista; import torch; print('Installation successful')"

5. Start JupyterLab

jupyter lab

Open one of the two notebooks from the JupyterLab file browser.

Quick Start

For a first test, open:

gempy_finite_fault_test_one_fault.ipynb

This is the smaller example and is useful for understanding the finite-fault configuration.

After that, open:

finite_fault_Relay_ramp_gempy_function.ipynb

This notebook demonstrates:

Multiple finite faults.
Separate structural groups for each fault.
A fault-relations matrix.
Individual ellipsoidal masks.
PyTorch-based model computation.
Visualisation of the finite-fault ellipsoids.

Run the notebook cells from top to bottom.

Example 1: Single Finite Fault

The single-fault notebook starts from GemPy's built-in example model:

from gempy.core.data.enumerators import ExampleModel
import gempy as gp

geo_model = gp.generate_example_model(
    example_model=ExampleModel.ONE_FAULT,
    compute_model=False,
)

Define the ellipsoid

import numpy as np

center = np.array([500, 500, 300])
radius = np.array([100, 100, 200]) / 0.8
max_slope = np.array([1.0, 1.0, 1.0]) * 0.1

The values are given in the coordinate system of the geological model.

Transform the parameters

scaled_center = geo_model.input_transform.apply(
    center.reshape(1, -1)
)[0]

scaled_radius = geo_model.input_transform.scale_points(
    radius.reshape(1, -1)
)[0]

Create the implicit ellipsoid

scalar_function = gp.implicit_functions.ellipsoid_3d_factory(
    center=scaled_center,
    radius=scaled_radius,
    max_slope=max_slope,
)

Define an additional transformation

transform = gp.data.Transform(
    position=np.array([0.0, 0.0, 0.0]),
    rotation=np.array([0.0, 60.0, 0.0]),
    scale=np.ones(3),
)

The example rotates the finite-fault mask by 60 degrees around the second coordinate axis.

Construct the finite-fault data

finite_fault_data = gp.data.FiniteFaultData(
    implicit_function=scalar_function,
    implicit_function_transform=transform,
    pivot=scaled_center,
)

Attach it to the fault group

faults_data = gp.data.FaultsData(
    fault_values_everywhere=np.zeros(0),
    fault_values_on_sp=np.zeros(0),
    thickness=None,
    fault_values_ref=np.zeros(0),
    fault_values_rest=np.zeros(0),
    finite_fault_data=finite_fault_data,
)

geo_model.structural_frame.structural_groups[
    0
].faults_input_data = faults_data

Compute the model

solution = gp.compute_model(geo_model)

The fault displacement is then spatially controlled by the ellipsoidal mask.

Example 2: Relay Ramp with Four Finite Faults

The relay-ramp notebook creates a model with four faults and four stratigraphic units.

Create the model

data = gp.create_geomodel(
    project_name="fault",
    extent=[0, 10000, 0, 8000, 0, 5000],
    resolution=[50, 20, 50],
    importer_helper=gp.data.ImporterHelper(
        path_to_orientations="Relay_ramp_obs_orientation.csv",
        path_to_surface_points="Relay_ramp_obs_points.csv",
    ),
)

Separate the faults into structural groups

Each fault is placed in its own structural group:

gp.map_stack_to_surfaces(
    gempy_model=data,
    mapping_object={
        "Fault_f1": ("f4",),
        "Fault_f2": ("f3",),
        "Fault_f3": ("f2",),
        "Fault_f4": ("f1",),
        "Strat_Series1": ("4", "3", "2", "1"),
    },
)

Using one structural group per fault is important because every group can then receive its own FiniteFaultData object and therefore its own ellipsoid.

Mark the groups as faults

from gempy_engine.core.data.stack_relation_type import StackRelationType

structural_frame = data.structural_frame

for i in range(4):
    structural_frame.structural_groups[
        i
    ].structural_relation = StackRelationType.FAULT

Define fault relationships

n_groups = len(structural_frame.structural_groups)
n_faults = 4

fault_relations = np.zeros(
    (n_groups, n_groups),
    dtype=int,
)

for i in range(n_faults):
    fault_relations[i, i + 1:] = 1

structural_frame.fault_relations = fault_relations

Define the finite-fault ellipsoids

ellipsoids = [
    (
        3,
        np.array([1000, 0, 5000]),
        np.array([5000, 4000, 15000]),
    ),
    (
        2,
        np.array([4250, 3500, 5000]),
        np.array([3000, 2000, 10000]),
    ),
    (
        1,
        np.array([6000, 5000, 5000]),
        np.array([3000, 2000, 10000]),
    ),
    (
        0,
        np.array([8500, 8000, 5000]),
        np.array([4000, 4000, 10000]),
    ),
]

Each tuple contains:

(structural_group_index, centre_xyz, radius_xyz)

Create and assign each mask

max_slope = np.array([1.0, 1.0, 1.0]) * 1.3

for group_index, center, radius in ellipsoids:
    scaled_center = data.input_transform.apply(
        center.reshape(1, -1)
    )[0]

    scaled_radius = data.input_transform.scale_points(
        radius.reshape(1, -1)
    )[0]

    scalar_function = gp.implicit_functions.ellipsoid_3d_factory(
        center=scaled_center,
        radius=scaled_radius,
        max_slope=max_slope,
    )

    transform = gp.data.Transform(
        position=np.array([0.0, 0.0, 0.0]),
        rotation=np.array([0.0, 45.0, 0.0]),
        scale=np.ones(3) * 0.5,
    )

    finite_fault_data = gp.data.FiniteFaultData(
        implicit_function=scalar_function,
        implicit_function_transform=transform,
        pivot=scaled_center,
    )

    faults_data = gp.data.FaultsData(
        fault_values_everywhere=np.zeros(0),
        fault_values_on_sp=np.zeros(0),
        thickness=None,
        fault_values_ref=np.zeros(0),
        fault_values_rest=np.zeros(0),
        finite_fault_data=finite_fault_data,
    )

    structural_frame.structural_groups[
        group_index
    ].faults_input_data = faults_data

Compute using the PyTorch backend

solution = gp.compute_model(
    data,
    engine_config=gp.data.GemPyEngineConfig(
        backend=gp.data.AvailableBackends.PYTORCH,
        dtype="float64",
        use_gpu=True,
    ),
    to_numpy=True,
)

Set use_gpu=False when no compatible GPU is available.

Input Data Format

Surface points

X,Y,Z,formation

Example:

X,Y,Z,formation
0,200,600,rock1
200,500,600,rock1
450,500,400,fault

Column	Meaning
`X`	X coordinate of the geological interface point.
`Y`	Y coordinate of the geological interface point.
`Z`	Z coordinate of the geological interface point.
`formation`	Name of the geological surface or fault.

Orientations

X,Y,Z,azimuth,dip,polarity,formation

Example:

X,Y,Z,azimuth,dip,polarity,formation
450,500,400,90,60,1,fault

Column	Meaning
`X`, `Y`, `Z`	Position of the orientation measurement.
`azimuth`	Azimuth of the geological orientation.
`dip`	Dip angle.
`polarity`	Orientation polarity.
`formation`	Geological surface or fault represented by the measurement.

Formation names in the CSV files must match the names used in the stack mapping.

Modelling Workflow

A general finite-fault workflow consists of the following steps:

Prepare geological surface-point and orientation data.
Define the model domain and regular-grid resolution.
Create the GemPy model.
Map geological surfaces to structural groups.
Put faults requiring independent finite geometries into separate groups.
Mark the relevant groups as StackRelationType.FAULT.
Define the fault-relations matrix.
Define one ellipsoidal mask for each finite fault.
Attach the corresponding FiniteFaultData to each fault group.
Compute and visualise the model.

Important Parameters

Model parameters

Parameter	Description
`extent`	Minimum and maximum model coordinates in X, Y, and Z.
`resolution`	Resolution of the regular grid in the three coordinate directions.
`project_name`	Name assigned to the GemPy model.
`path_to_orientations`	Path to the orientation CSV file.
`path_to_surface_points`	Path to the surface-point CSV file.

Ellipsoid parameters

Parameter	Description
`center`	Centre of the ellipsoid in model coordinates.
`radius`	Semi-axis lengths in X, Y, and Z before the additional transformation.
`max_slope`	Controls the transition behaviour of the ellipsoidal implicit function.
`position`	Additional translation applied by `gp.data.Transform`.
`rotation`	Rotation angles passed to the finite-fault transformation.
`scale`	Additional scale applied to the finite-fault mask.
`pivot`	Point about which the transformation is applied.

Computation parameters

Parameter	Description
`backend`	Numerical backend used by GemPy Engine.
`dtype`	Floating-point precision, such as `float32` or `float64`.
`use_gpu`	Enables compatible GPU execution when using the PyTorch backend.
`to_numpy`	Converts the computed result arrays to NumPy arrays.

Fault Relations

The fault-relations matrix controls which structural groups are affected by each fault.

fault_relations[source_group, affected_group]

A true or non-zero entry indicates that the group represented by the row affects the group represented by the column.

Example:

fault_relations = np.array([
    [0, 1, 1],
    [0, 0, 1],
    [0, 0, 0],
])

This means:

Group 0 affects groups 1 and 2.
Group 1 affects group 2.
Group 2 affects no other group.

The matrix represents geological chronology and should be adapted to the intended model.

Inspect the structural-group order using:

for index, group in enumerate(
    model.structural_frame.structural_groups
):
    print(index, group.name, group.structural_relation)

Visualisation

GemPy Viewer

viewer = gpv.plot_3d(
    data,
    show=False,
    show_data=False,
    show_lith=False,
    show_boundaries=True,
    show_result=True,
)

plotter = viewer.p
plotter.show()

Finite-fault ellipsoids

The relay-ramp notebook contains a helper function:

plot_finite_fault_ellipsoid(...)

It reconstructs the ellipsoid, applies the finite-fault transform, converts the points back to model coordinates, and adds the surface to a PyVista plotter.

for i in range(4):
    finite_fault_data = (
        structural_frame
        .structural_groups[i]
        .faults_input_data
        .finite_fault_data
    )

    plot_finite_fault_ellipsoid(
        plotter=plotter,
        model=data,
        finite_fault_data=finite_fault_data,
        opacity=0.2,
    )

These transparent ellipsoids are diagnostic geometry and help validate the spatial extent assigned to each fault.

Using Your Own Data

To apply the workflow to another geological model:

Replace the example surface-point and orientation CSV files.
Update the model extent and resolution.
Map your surfaces and faults to structural groups.
Put independently restricted faults into separate groups.
Define the geological fault chronology.
Estimate the centre, radii, rotation, and scale for each finite-fault ellipsoid.
Attach one FiniteFaultData object to each fault group.
Plot the ellipsoids together with the geological model.
Perform sensitivity tests for mask dimensions and transformations.

All geological data and ellipsoid parameters must use a consistent coordinate system.

CPU and GPU Execution

The relay-ramp notebook uses the PyTorch backend with GPU execution:

backend=gp.data.AvailableBackends.PYTORCH
use_gpu=True

Verify GPU availability using:

import torch

print("PyTorch version:", torch.__version__)
print("CUDA available:", torch.cuda.is_available())

if torch.cuda.is_available():
    print("GPU:", torch.cuda.get_device_name(0))

For CPU execution:

solution = gp.compute_model(
    data,
    engine_config=gp.data.GemPyEngineConfig(
        backend=gp.data.AvailableBackends.PYTORCH,
        dtype="float64",
        use_gpu=False,
    ),
    to_numpy=True,
)

Known Limitations

The repository currently consists of notebooks and example datasets rather than an installable Python package.
There is currently no command-line interface.
Dependency versions are not pinned.
Some imports access gempy_engine internals directly.
Internal GemPy APIs may change between versions.
Finite-fault ellipsoid parameters must currently be selected manually.
The ellipsoids are not automatically inferred from fault observations.
An ellipsoid is an approximation and may not reproduce complex fault-tip geometry.
Interactive PyVista rendering may require additional configuration in remote or headless environments.
High grid resolutions may require substantial memory.
The example fault-relations matrix represents one specific chronology and is not universally applicable.

Troubleshooting

CSV files cannot be found

Start JupyterLab from the repository root:

cd GemPy_finite_fault
jupyter lab

Alternatively, replace the relative file paths with absolute paths.

`gempy_engine` import errors

Check the installed versions:

python -m pip show gempy
python -m pip show gempy-engine
python -m pip show gempy-viewer

The installed GemPy version may not provide the same internal API used by the notebooks.

CUDA is not available

Change:

use_gpu=True

to:

use_gpu=False

Model computation requires too much memory

Reduce the regular-grid resolution:

resolution=[30, 15, 30]

Increase it gradually after verifying that the workflow runs correctly.

The ellipsoid appears in the wrong position

Check that:

The centre is defined in model coordinates.
input_transform.apply(...) is used for the centre.
input_transform.scale_points(...) is used for the radius.
The transform uses the intended pivot.
The group index matches the intended fault.
All data use a consistent coordinate system.

The wrong fault receives a mask

Print the structural-group order:

for i, group in enumerate(
    model.structural_frame.structural_groups
):
    print(i, group.name)

The first value in each ellipsoid tuple is the structural-group index.

Recommended Future Improvements

Possible extensions include:

A pinned requirements.txt or environment.yml.
A reusable Python function for assigning finite-fault masks.
An installable Python package.
Automated ellipsoid estimation from fault observations.
Support for non-ellipsoidal masks.
Multiple masks for segmented or branching faults.
Automated validation of structural-group indices.
Unit tests.
VTK or VTU export.
Reproducible parameter files.
Additional real-data examples.
Continuous integration.
Quantitative comparison between finite and unrestricted fault displacement.

Contributing

Contributions are welcome, particularly for:

Additional finite-fault examples.
Improved GemPy version compatibility.
Alternative implicit finite-fault functions.
Automated parameter estimation.
Documentation and tutorials.
Tests and reproducible environments.
Improved three-dimensional visualisation.

Suggested workflow:

git checkout -b feature/my-improvement
git add .
git commit -m "Describe the improvement"
git push origin feature/my-improvement

Then open a pull request describing the change, the GemPy and Python versions used, and how the result was validated.

Acknowledgements

This repository builds on the open-source GemPy ecosystem:

GemPy for implicit geological modelling.
GemPy Viewer for geological visualisation.
GemPy Engine for numerical model computation.
PyVista for three-dimensional visualisation.
PyTorch for optional tensor and GPU-based computation.

Please cite GemPy and its associated publications when using this repository in scientific work.

Contact

For questions, suggestions, or bug reports, please open an issue:

https://github.com/yangjiandendi/GemPy_finite_fault/issues

Name		Name	Last commit message	Last commit date
README.md		README.md
Relay_ramp_obs_orientation.csv		Relay_ramp_obs_orientation.csv
Relay_ramp_obs_points.csv		Relay_ramp_obs_points.csv
finite_fault_Relay_ramp_gempy_function.ipynb		finite_fault_Relay_ramp_gempy_function.ipynb
gempy_finite_fault_test_one_fault.ipynb		gempy_finite_fault_test_one_fault.ipynb
model5_orientations.csv		model5_orientations.csv
model5_surface_points.csv		model5_surface_points.csv

Latest commit

History

Folders and files

Repository files navigation

GemPy Finite Fault

Table of Contents

Background

Relationship to GemPy

Finite-Fault Concept

Repository Structure

Files

Installation

1. Clone the repository

2. Create a dedicated environment

3. Install the required packages

4. Verify the installation

5. Start JupyterLab

Quick Start

Example 1: Single Finite Fault

Define the ellipsoid

Transform the parameters

Create the implicit ellipsoid

Define an additional transformation

Construct the finite-fault data

Attach it to the fault group

Compute the model

Example 2: Relay Ramp with Four Finite Faults

Create the model

Separate the faults into structural groups

Mark the groups as faults

Define fault relationships

Define the finite-fault ellipsoids

Create and assign each mask

Compute using the PyTorch backend

Input Data Format

Surface points

Orientations

Modelling Workflow

Important Parameters

Model parameters

Ellipsoid parameters

Computation parameters

Fault Relations

Visualisation

GemPy Viewer

Finite-fault ellipsoids

Using Your Own Data

CPU and GPU Execution

Known Limitations

Troubleshooting

CSV files cannot be found

gempy_engine import errors

CUDA is not available

Model computation requires too much memory

The ellipsoid appears in the wrong position

The wrong fault receives a mask

Recommended Future Improvements

Contributing

Acknowledgements

Contact

About

Resources

Stars

Watchers

Forks

Releases

Packages

Contributors

Languages

`gempy_engine` import errors