From eebd65fe1a667b9a422d850939a78a16f78d6284 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@gmail.com>
Date: Mon, 28 Apr 2025 14:11:12 +1000
Subject: [PATCH 01/12] docs: update example 2

---
 examples/1_basic/plot_2_surface_modelling.py | 129 +++++++------------
 1 file changed, 45 insertions(+), 84 deletions(-)

diff --git a/examples/1_basic/plot_2_surface_modelling.py b/examples/1_basic/plot_2_surface_modelling.py
index 2e47a743..e01c61c5 100644
--- a/examples/1_basic/plot_2_surface_modelling.py
+++ b/examples/1_basic/plot_2_surface_modelling.py
@@ -41,53 +41,41 @@
 
 import numpy as np
 
-
 ######################################################################
+# Load Example Data
+# ~~~~~~~~~~~~~~~~~
 # The data for this example can be imported from the example datasets
-# module in loopstructural.
-#
+# module in LoopStructural. The `load_claudius` function provides both
+# the data and the bounding box (bb) for the model.
 
 data, bb = load_claudius()
 
-data["val"].unique()
-
+# Display unique scalar field values in the dataset
+print(data["val"].unique())
 
 ######################################################################
-# GeologicalModel
-# ~~~~~~~~~~~~~~~
-#
-# To create a ``GeologicalModel`` the origin (lower left) and max extent
-# (upper right) of the model area need to be specified. In this example
-# these are provided in the bb array.
-#
+# Create Geological Model
+# ~~~~~~~~~~~~~~~~~~~~~~~~
+# Define the model area using the bounding box (bb) and create a
+# GeologicalModel instance.
 
 model = GeologicalModel(bb[0, :], bb[1, :])
 
-
 ######################################################################
-# A pandas dataframe with appropriate columns can be used to link the data
-# to the geological model.
-#
-# * ``X`` is the x coordinate
-# * ``Y`` is the y # coordinate
-# * ``Z`` is the z coordinate
-# * ``feature_name`` is a name to link the data to a model object
-# * ``val`` is the value of the scalar field which represents the
-# distance from a reference horizon. It is comparable
-# to the relative thickness
-#
-# * ``nx`` is the x component of the normal vector to the surface gradient
-# * ``ny`` is the y component of the normal vector to the surface gradient
-# * ``nz`` is the z component of the normal vector to the surface gradeint
-# * ``strike`` is the strike angle
-# * ``dip`` is the dip angle
-#
-# Having a look at the data for this example by looking at the top of the
-# dataframe and then using a 3D plot
-#
-
-data["feature_name"].unique()
-
+# Link Data to Geological Model
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# A pandas dataframe with appropriate columns is used to link the data
+# to the geological model. Key columns include:
+# * `X`, `Y`, `Z`: Coordinates of the observation
+# * `feature_name`: Name linking the data to a model object
+# * `val`: Scalar field value representing distance from a reference horizon
+# * `nx`, `ny`, `nz`: Components of the normal vector to the surface gradient
+# * `strike`, `dip`: Strike and dip angles
+
+# Display unique feature names in the dataset
+print(data["feature_name"].unique())
+
+# Visualize the data points and orientation vectors in 3D
 viewer = Loop3DView(background="white")
 viewer.add_points(
     data[~np.isnan(data["val"])][["X", "Y", "Z"]].values,
@@ -99,50 +87,18 @@
 )
 viewer.display()
 
-
-######################################################################
-# The pandas dataframe can be linked to the ``GeologicalModel`` using
-# ``.set_model_data(dataframe``
-#
-
+# Link the data to the geological model
 model.set_model_data(data)
 
-
 ######################################################################
-# The ``GeologicalModel`` contains an ordered list of the different
-# geological features within a model and how these features interact. This
-# controls the topology of the different geological features in the model.
-# Different geological features can be added to the geological model such
-# as:
-# * Foliations
-# * Faults
-# * Unconformities
-# * Folded foliations
-# *  Structural Frames
-#
-# In this example we will only add a foliation using the function
-#
-# .. code:: python
-#
-#    model.create_and_add_foliation(name)
-#
-# where name is the name in the ``feature_name`` field, other parameters we
-# specified are the:
-# * ``interpolatortype`` - we can either use a
-# PiecewiseLinearInterpolator ``PLI``, a FiniteDifferenceInterpolator
-# ``FDI`` or a radial basis interpolator ``surfe``
-# * ``nelements - int`` is the how many elements are used to discretize the resulting solution
-# * ``buffer - float`` buffer percentage around the model area
-# * ``solver`` - the algorithm to solve the least squares problem e.g.
-# ``lu`` for lower upper decomposition, ``cg`` for conjugate gradient,
-# ``pyamg`` for an algorithmic multigrid solver
-# * ``damp - bool`` - whether to add a small number to the diagonal of the interpolation
-# matrix for discrete interpolators - this can help speed up the solver
-# and makes the solution more stable for some interpolators
-#
+# Add Geological Features
+# ~~~~~~~~~~~~~~~~~~~~~~~~
+# The GeologicalModel can include various geological features such as
+# foliations, faults, unconformities, and folds. In this example, we
+# add a foliation using the `create_and_add_foliation` method.
 
+# Define stratigraphic column with scalar field ranges for each unit
 vals = [0, 60, 250, 330, 600]
-
 strat_column = {"strati": {}}
 for i in range(len(vals) - 1):
     strat_column["strati"]["unit_{}".format(i)] = {
@@ -151,27 +107,32 @@
         "id": i,
     }
 
-
+# Set the stratigraphic column in the model
 model.set_stratigraphic_column(strat_column)
 
+# Add a foliation to the model
 strati = model.create_and_add_foliation(
     "strati",
-    interpolatortype="FDI",  # try changing this to 'PLI'
-    nelements=1e4,  # try changing between 1e3 and 5e4
-    buffer=0.3,
-    damp=True,
+    interpolatortype="FDI",  # Finite Difference Interpolator
+    nelements=int(1e4),  # Number of elements for discretization
+    buffer=0.3,  # Buffer percentage around the model area
+    damp=True,  # Add damping for stability
 )
+
 ######################################################################
-# Plot the surfaces
-# ------------------------------------
+# Visualize Model Surfaces
+# ~~~~~~~~~~~~~~~~~~~~~~~~
+# Plot the surfaces of the geological model using a 3D viewer.
 
 viewer = Loop3DView(model)
 viewer.plot_model_surfaces(cmap="tab20")
 viewer.display()
 
 ######################################################################
-# Plot block diagram
-# -------------------
+# Visualize Block Diagram
+# ~~~~~~~~~~~~~~~~~~~~~~~
+# Plot a block diagram of the geological model to visualize the
+# stratigraphic units in 3D.
 
 viewer = Loop3DView(model)
 viewer.plot_block_model(cmap="tab20")

From fbd8430a960b54a74979378d1da5a0a02085d21d Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@gmail.com>
Date: Mon, 28 Apr 2025 14:12:47 +1000
Subject: [PATCH 02/12] update formatting, docstrings and ensure that fault
 without normals works.

---
 .../features/builders/_fault_builder.py       | 255 ++++++++----------
 1 file changed, 110 insertions(+), 145 deletions(-)

diff --git a/LoopStructural/modelling/features/builders/_fault_builder.py b/LoopStructural/modelling/features/builders/_fault_builder.py
index a22fc689..563f8f5a 100644
--- a/LoopStructural/modelling/features/builders/_fault_builder.py
+++ b/LoopStructural/modelling/features/builders/_fault_builder.py
@@ -3,7 +3,6 @@
 from LoopStructural.utils.maths import rotation
 from ._structural_frame_builder import StructuralFrameBuilder
 from .. import AnalyticalGeologicalFeature
-from LoopStructural.utils import get_vectors
 import numpy as np
 import pandas as pd
 from ....utils import getLogger
@@ -52,8 +51,11 @@ def __init__(
         self.frame.model = model
         self.origin = np.array([np.nan, np.nan, np.nan])
         self.maximum = np.array([np.nan, np.nan, np.nan])  # self.model.bounding_box[1, :]
+
+        if bounding_box is None:
+            raise ValueError("BoundingBox cannot be None")
+
         # define a maximum area to mesh adding buffer to model
-        # buffer = .2
         self.minimum_origin = bounding_box.with_buffer(fault_bounding_box_buffer).origin
         self.maximum_maximum = bounding_box.with_buffer(fault_bounding_box_buffer).maximum
 
@@ -66,8 +68,23 @@ def __init__(
         self.fault_centre = None
 
     def update_geometry(self, points):
+        """
+        Update the geometry of the fault by adjusting the origin and maximum bounds
+        based on the provided points.
+
+        Parameters
+        ----------
+        points : numpy.ndarray
+            Array of points used to update the fault geometry.
+        """
+        if self.origin is None or self.maximum is None:
+            raise ValueError("Origin and maximum must be initialized before updating geometry.")
+
         self.origin = np.nanmin(np.array([np.min(points, axis=0), self.origin]), axis=0)
         self.maximum = np.nanmax(np.array([np.max(points, axis=0), self.maximum]), axis=0)
+        # add a small buffer 10% of current length to the origin and maximum
+        self.origin = self.origin - 0.1 * (self.maximum - self.origin)
+        self.maximum = self.maximum + 0.1 * (self.maximum - self.origin)
         self.origin[self.origin < self.minimum_origin] = self.minimum_origin[
             self.origin < self.minimum_origin
         ]
@@ -93,144 +110,87 @@ def create_data_from_geometry(
         fault_dip_anisotropy=1.0,
         fault_pitch=None,
     ):
-        """Generate the required data for building a fault frame for a fault with the
-        specified parameters
+        """
+        Generate the required data for building a fault frame with the specified parameters.
 
         Parameters
         ----------
-        data : DataFrame,
-            model data
-        fault_center : np.array(3)
-            x,y,z coordinates of the fault center
-        normal_vector : np.array(3)
-            x,y,z components of normal vector to fault, single observation usually
-            average direction
-        slip_vector : np.array(3)
-            x,y,z components of slip vector for the fault, single observation usually
-            average direction
-        minor_axis : double
-            distance away from fault for the fault volume
-        major_axis : double
-            fault extent
-        intermediate_axis : double
-            fault volume radius in the slip direction
+        fault_frame_data : pandas.DataFrame
+            DataFrame containing fault frame data.
+        fault_center : array-like, optional
+            Coordinates of the fault center.
+        fault_normal_vector : array-like, optional
+            Normal vector of the fault.
+        fault_slip_vector : array-like, optional
+            Slip vector of the fault.
+        minor_axis : float, optional
+            Minor axis length of the fault.
+        major_axis : float, optional
+            Major axis length of the fault.
+        intermediate_axis : float, optional
+            Intermediate axis length of the fault.
+        w : float, default=1.0
+            Weighting factor for the fault data.
+        points : bool, default=False
+            Whether to include points in the fault data.
+        force_mesh_geometry : bool, default=False
+            Whether to force the use of mesh geometry.
+        fault_buffer : float, default=0.2
+            Buffer size around the fault.
+        fault_trace_anisotropy : float, default=1.0
+            Anisotropy factor for the fault trace.
+        fault_dip : float, default=90
+            Dip angle of the fault in degrees.
+        fault_dip_anisotropy : float, default=1.0
+            Anisotropy factor for the fault dip.
+        fault_pitch : float, optional
+            Pitch angle of the fault.
         """
-        trace_mask = np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0)
-        logger.info(f"There are {np.sum(trace_mask)} points on the fault trace")
-        if np.sum(trace_mask) == 0 and fault_center is None:
-            logger.error("You cannot model a fault without defining the location of the fault")
-            raise ValueError("There are no points on the fault trace")
-        # find the middle point on the fault trace if center is not provided
-        if fault_center is None:
-            trace_mask = np.logical_and(
-                fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0
-            )
-            fault_center = fault_frame_data.loc[trace_mask, ["X", "Y", "Z"]].mean(axis=0).to_numpy()
-            dist = np.linalg.norm(
-                fault_center - fault_frame_data.loc[trace_mask, ["X", "Y", "Z"]].to_numpy(), axis=1
-            )
-            # make the nan points greater than the max dist 10 is arbitrary and doesn't matter
-            dist[np.isnan(dist)] = np.nanmax(dist) + 10
-            fault_center = fault_frame_data.loc[trace_mask, ["X", "Y", "Z"]].to_numpy()[
-                np.argmin(dist), :
-            ]
-        # get all of the gradient data associated with the fault trace
+        fault_trace = fault_frame_data.loc[
+            np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0),
+            ["X", "Y"],
+        ].to_numpy()
         if fault_normal_vector is None:
-            gradient_mask = np.logical_and(
-                fault_frame_data["coord"] == 0, ~np.isnan(fault_frame_data["gz"])
-            )
-            vector_data = fault_frame_data.loc[gradient_mask, ["gx", "gy", "gz"]].to_numpy()
-            normal_mask = np.logical_and(
-                fault_frame_data["coord"] == 0, ~np.isnan(fault_frame_data["nz"])
-            )
-            vector_data = np.vstack(
-                [
-                    vector_data,
-                    fault_frame_data.loc[normal_mask, ["nx", "ny", "nz"]].to_numpy(),
-                ]
-            )
-
-            if len(vector_data) == 0:
-                logger.warning(
-                    f"No orientation data for fault\n\
-                    Defaulting to a dip of {fault_dip}vertical fault"
-                )
-                # if the line is long enough, estimate the normal vector
-                # by finding the centre point of the line and calculating the tangnent
-                # of the two points
-                if fault_frame_data.loc[trace_mask, :].shape[0] > 3:
-
-                    pts = fault_frame_data.loc[trace_mask, ["X", "Y", "Z"]].to_numpy()
-                    dist = np.abs(np.linalg.norm(fault_center - pts, axis=1))
-                    # any nans just make them max distance + a bit
-                    dist[np.isnan(dist)] = np.nanmax(dist) + 10
-                    # idx = np.argsort(dist)
-                    # direction_vector = pts[idx[-1]] - pts[idx[-2]]
-                    # coefficients = np.polyfit(
-                    #     fault_frame_data.loc[trace_mask, "X"],
-                    #     fault_frame_data.loc[trace_mask, "Y"],
-                    #     1,
-                    # )
-                    # slope, intercept = coefficients
-                    slope, intercept = np.polyfit(
-                        pts[dist < 0.25 * np.nanmax(dist), 0],
-                        pts[dist < 0.25 * np.nanmax(dist), 1],
-                        1,
-                    )
-
-                    # # Create a direction vector using the slope
-                    direction_vector = np.array([1, slope, 0])
-                    direction_vector /= np.linalg.norm(direction_vector)
-                    rotation_matrix = rotation(direction_vector[None, :], [90 - fault_dip])
-                    vector_data = np.array(
-                        [
-                            [
-                                direction_vector[1],
-                                -direction_vector[0],
-                                0,
-                            ]
-                        ]
-                    )
-                    vector_data /= np.linalg.norm(vector_data, axis=1)
-                    vector_data = np.einsum("ijk,ik->ij", rotation_matrix, vector_data)
+            # Calculate fault strike using least squares fit
+            X = fault_trace[:, 0].reshape(-1, 1)
+            Y = fault_trace[:, 1]
+            # Fit line Y = mX + b
+            A = np.hstack([X, np.ones_like(X)])
+            m, _ = np.linalg.lstsq(A, Y, rcond=None)[0]
+            # Convert slope to strike vector
+            strike_vector = np.array([1, m, 0])
+            strike_vector /= np.linalg.norm(strike_vector)
+            fault_normal_vector = np.cross(strike_vector, [0, 0, 1])
 
-                    vector_data /= np.linalg.norm(vector_data, axis=1)
-            fault_normal_vector = np.mean(vector_data, axis=0)
+        if not isinstance(fault_normal_vector, np.ndarray):
+            fault_normal_vector = np.array(fault_normal_vector)
 
-        logger.info(f"Fault normal vector: {fault_normal_vector}")
+        if fault_pitch is not None:
+            rotation_matrix = rotation(fault_normal_vector[None, :], np.array([fault_pitch]))
+            fault_slip_vector = np.einsum("ijk,ik->ij", rotation_matrix, fault_normal_vector[None, :])[0]
 
-        # estimate the fault slip vector
         if fault_slip_vector is None:
-            slip_mask = np.logical_and(
-                fault_frame_data["coord"] == 1, ~np.isnan(fault_frame_data["gz"])
-            )
-            fault_slip_data = fault_frame_data.loc[slip_mask, ["gx", "gy", "gz"]]
-
-            if len(fault_slip_data) == 0:
-                logger.warning(
-                    "There is no slip vector data for the fault, using vertical slip vector\n\
-                          projected onto fault surface estimating from fault normal"
-                )
-                strike_vector, dip_vector = get_vectors(fault_normal_vector[None, :])
-                fault_slip_vector = dip_vector[:, 0]
-                if fault_pitch is not None:
-                    print('using pitch')
-                    rotm = rotation(fault_normal_vector[None,:],[fault_pitch])
-                    print(rotm.shape,fault_slip_vector.shape)
-                    fault_slip_vector = np.einsum("ijk,k->ij", rotm, fault_slip_vector)[0,:]
-                logger.info(f"Estimated fault slip vector: {fault_slip_vector}")
-            else:
-                fault_slip_vector = fault_slip_data.mean(axis=0).to_numpy()
-
-        self.fault_normal_vector = fault_normal_vector
-        self.fault_slip_vector = fault_slip_vector
-
-        self.fault_centre = fault_center
-        if major_axis is None:
+            fault_slip_vector = np.cross(fault_normal_vector, [1., 0., 0.])
+            if np.linalg.norm(fault_slip_vector) == 0:
+                fault_slip_vector = np.cross(fault_normal_vector, [0., 1., 0.])
+            fault_slip_vector /= np.linalg.norm(fault_slip_vector)
+        if fault_center is None:
             fault_trace = fault_frame_data.loc[
                 np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0),
                 ["X", "Y"],
             ].to_numpy()
+            fault_center = fault_trace.mean(axis=0)
+            fault_center = np.array([fault_center[0], fault_center[1], 0.0])
+        if not isinstance(fault_center, np.ndarray):
+            fault_center = np.array(fault_center)
+        if fault_center.shape[0] != 3:
+            raise ValueError("fault_center must be a 3 element array")
+        self.fault_normal_vector = fault_normal_vector / np.linalg.norm(fault_normal_vector)
+        self.fault_slip_vector = fault_slip_vector / np.linalg.norm(fault_slip_vector)
+
+        self.fault_centre = fault_center
+        if major_axis is None:
+
             distance = np.linalg.norm(fault_trace[:, None, :] - fault_trace[None, :, :], axis=2)
             if len(distance) == 0 or np.sum(distance) == 0:
                 logger.warning("There is no fault trace for {}".format(self.name))
@@ -382,6 +342,7 @@ def create_data_from_geometry(
                             0,
                             w,
                         ]
+
             if major_axis is not None:
                 fault_tips[0, :] = fault_center[:3] + strike_vector * 0.5 * major_axis
                 fault_tips[1, :] = fault_center[:3] - strike_vector * 0.5 * major_axis
@@ -503,16 +464,14 @@ def set_mesh_geometry(self, buffer, rotation):
             self.origin - length * buffer, self.maximum + length * buffer, rotation
         )
 
-    def add_splay(self, splay, splayregion=None):
-        if splayregion is None:
+    def add_splay(self, splay, splay_region=None):
+        if splay_region is None:
 
-            def splayregion(xyz):
+            def default_splay_region(xyz):
                 pts = (
-                    self.builders[0].data[["X", "Y", "Z", "val"]].to_numpy()
+                    self.builders[0].data["X", "Y", "Z", "val"].to_numpy()
                 )  # get_value_constraints()
                 pts = pts[pts[:, 3] == 0, :]
-                # check whether the fault is on the hanging wall or footwall of splay fault
-
                 ext_field = splay[2].evaluate_value(pts[:, :3])
                 surf_field = splay[0].evaluate_value(pts[:, :3])
                 intersection_value = ext_field[np.nanargmin(np.abs(surf_field))]
@@ -531,19 +490,24 @@ def splayregion(xyz):
                     )
                     return mask
 
+            splay_region = default_splay_region
+
         scalefactor = splay.fault_major_axis / self.fault_major_axis
-        self.builders[0].add_equality_constraints(splay, splayregion, scalefactor)
-        return splayregion
+        self.builders[0].add_equality_constraints(splay, splay_region, scalefactor)
+        return splay_region
 
     def add_fault_trace_anisotropy(self, w: float = 1.0):
-        """_summary_
+        """
+        Add fault trace anisotropy to the model.
 
         Parameters
         ----------
         w : float, optional
-            _description_, by default 1.0
+            Weighting factor for the anisotropy, by default 1.0
         """
         if w > 0:
+            if self.fault_normal_vector is None:
+                raise ValueError("fault_normal_vector must be initialized before adding anisotropy.")
 
             plane = np.array([0, 0, 1])
             strike_vector = (
@@ -553,24 +517,25 @@ def add_fault_trace_anisotropy(self, w: float = 1.0):
             strike_vector = np.array([strike_vector[1], -strike_vector[0], 0])
 
             anisotropy_feature = AnalyticalGeologicalFeature(
-                vector=strike_vector, origin=[0, 0, 0], name="fault_trace_anisotropy"
+                vector=strike_vector, origin=np.array([0, 0, 0]), name="fault_trace_anisotropy"
             )
-            # print('adding fault trace anisotropy')
             self.builders[0].add_orthogonal_feature(
                 anisotropy_feature, w=w, region=None, step=1, B=0
             )
 
     def add_fault_dip_anisotropy(self, w: float = 1.0):
-        """_summary_
+        """
+        Add fault dip anisotropy to the model.
 
         Parameters
         ----------
-        dip : np.ndarray
-            _description_
         w : float, optional
-            _description_, by default 1.0
+            Weighting factor for the anisotropy, by default 1.0
         """
         if w > 0:
+            if self.fault_normal_vector is None:
+                raise ValueError("fault_normal_vector must be initialized before adding anisotropy.")
+
             plane = np.array([0, 0, 1])
             strike_vector = (
                 self.fault_normal_vector - np.dot(self.fault_normal_vector, plane) * plane
@@ -579,11 +544,11 @@ def add_fault_dip_anisotropy(self, w: float = 1.0):
             strike_vector = np.array([strike_vector[1], -strike_vector[0], 0])
 
             dip_vector = np.cross(strike_vector, self.fault_normal_vector)
+            dip_vector /= np.linalg.norm(dip_vector)
 
             anisotropy_feature = AnalyticalGeologicalFeature(
-                vector=dip_vector, origin=[0, 0, 0], name="fault_dip_anisotropy"
+                vector=dip_vector, origin=np.array([0, 0, 0]), name="fault_dip_anisotropy"
             )
-            # print(f'adding fault dip anisotropy {anisotropy_feature.name}')
             self.builders[0].add_orthogonal_feature(
                 anisotropy_feature, w=w, region=None, step=1, B=0
             )

From 4d6b925615916f484f423395d109a819a1af8606 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@gmail.com>
Date: Mon, 28 Apr 2025 14:13:22 +1000
Subject: [PATCH 03/12] add local origin when creating regular grid

---
 LoopStructural/datatypes/_bounding_box.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/LoopStructural/datatypes/_bounding_box.py b/LoopStructural/datatypes/_bounding_box.py
index dc1d01e8..9e9901d0 100644
--- a/LoopStructural/datatypes/_bounding_box.py
+++ b/LoopStructural/datatypes/_bounding_box.py
@@ -387,7 +387,7 @@ def regular_grid(
 
         if not local:
             coordinates = [
-                np.linspace(self.global_origin[i], self.global_maximum[i], nsteps[i])
+                np.linspace(self.global_origin[i]+self.origin[i], self.global_maximum[i], nsteps[i])
                 for i in range(self.dimensions)
             ]
         coordinate_grid = np.meshgrid(*coordinates, indexing="ij")

From 75532dd4763674d4e70bc3edbe7d030b3fedb4e3 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@gmail.com>
Date: Mon, 28 Apr 2025 14:14:03 +1000
Subject: [PATCH 04/12] clip smooth peak at -1,1 to ensure fault displacement
 stops at edge of ellipsoid

---
 LoopStructural/modelling/features/fault/_fault_function.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/LoopStructural/modelling/features/fault/_fault_function.py b/LoopStructural/modelling/features/fault/_fault_function.py
index b5d58043..d1a5e7fd 100644
--- a/LoopStructural/modelling/features/fault/_fault_function.py
+++ b/LoopStructural/modelling/features/fault/_fault_function.py
@@ -10,8 +10,10 @@
 
 
 def smooth_peak(x):
-    return 0.25 * x**6 + 0.5 * x**4 - 1.75 * x**2 + 1
-
+    v = np.zeros(x.shape)
+    mask = np.logical_and(x >= -1, x <= 1)
+    v[mask] = 0.25 * x[mask] ** 2 + 0.5 * x[mask] ** 4 - 1.75 * x[mask] ** 2 + 1
+    return v
 
 class FaultProfileFunction(metaclass=ABCMeta):
     def __init__(self):

From eb1d3e12b777ac8e5ab3c3324e3fc63b9e6435b0 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@gmail.com>
Date: Mon, 28 Apr 2025 14:14:13 +1000
Subject: [PATCH 05/12] apply faults to lambda features

---
 .../modelling/features/_lambda_geological_feature.py  | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/LoopStructural/modelling/features/_lambda_geological_feature.py b/LoopStructural/modelling/features/_lambda_geological_feature.py
index 3894c15d..487a2bf4 100644
--- a/LoopStructural/modelling/features/_lambda_geological_feature.py
+++ b/LoopStructural/modelling/features/_lambda_geological_feature.py
@@ -62,10 +62,13 @@ def evaluate_value(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray:
             _description_
         """
         v = np.zeros((pos.shape[0]))
-        if self.function is None:
-            v[:] = np.nan
-        else:
-            v[:] = self.function(pos)
+        v[:] = np.nan
+
+        # mask = self._calculate_mask(pos, ignore_regions=ignore_regions)
+        pos = self._apply_faults(pos)
+        if self.function is not None:
+            
+            v = self.function(pos)
         return v
 
     def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray:

From fc272199503fe4c5cce01a6353bb6900052064ec Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@gmail.com>
Date: Mon, 28 Apr 2025 14:14:57 +1000
Subject: [PATCH 06/12] ensure regions follow a structure, and postive/negative
 are compatible with isosurfacing

---
 LoopStructural/utils/regions.py | 37 +++++++++++++++++++++------------
 1 file changed, 24 insertions(+), 13 deletions(-)

diff --git a/LoopStructural/utils/regions.py b/LoopStructural/utils/regions.py
index 890d30d4..792f71da 100644
--- a/LoopStructural/utils/regions.py
+++ b/LoopStructural/utils/regions.py
@@ -1,9 +1,20 @@
 import numpy as np
+from abc import ABC, abstractmethod
 
 
-class BaseRegion:
-    def __init__(self):
-        self.type = "BaseRegion"
+class BaseRegion(ABC):
+    @abstractmethod
+    def __init__(self, feature, vector=None, point=None):
+        self.feature = feature
+        self.vector = vector
+        self.point = point
+        self.name = None
+        self.parent = None
+
+    @abstractmethod
+    def __call__(self, xyz) -> np.ndarray:
+        """Evaluate the region based on the input coordinates."""
+        pass
 
 
 class RegionEverywhere(BaseRegion):
@@ -24,18 +35,18 @@ def __call__(self, xyz):
         return self.function(xyz)
 
 
-class PositiveRegion:
+class PositiveRegion(BaseRegion):
     """Helper class for evaluating whether you are in the positive  region of a scalar field.
     If its outside of the support it will interpolate the average gradient at a point on the 0 isovalue
     and calculate the distance from this. Alternatively, a point and vector can be used to save computational time
     """
 
     def __init__(self, feature, vector=None, point=None):
-        self.feature = feature
-        self.vector = vector
-        self.point = point
+        super().__init__(feature, vector, point)
+        self.name = 'PositiveRegion'
+        self.parent = feature
 
-    def __call__(self, xyz):
+    def __call__(self, xyz) -> np.ndarray:
         val = self.feature.evaluate_value(xyz)
         # find a point on/near 0 isosurface
         if self.point is None:
@@ -62,18 +73,18 @@ def __call__(self, xyz):
         )
 
 
-class NegativeRegion:
+class NegativeRegion(BaseRegion):
     """Helper class for evaluating whether you are in the positive  region of a scalar field.
     If its outside of the support it will interpolate the average gradient at a point on the 0 isovalue
     and calculate the distance from this. Alternatively, a point and vector can be used to save computational time
     """
 
     def __init__(self, feature, vector=None, point=None):
-        self.feature = feature
-        self.vector = vector
-        self.point = point
+        super().__init__(feature, vector, point)
+        self.name = 'NegativeRegion'
+        self.parent = feature
 
-    def __call__(self, xyz):
+    def __call__(self, xyz) -> np.ndarray:
         val = self.feature.evaluate_value(xyz)
         # find a point on/near 0 isosurface
         if self.point is None:

From 8db324ca3abb237c53431d73dad302b95112a160 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@gmail.com>
Date: Mon, 28 Apr 2025 14:15:28 +1000
Subject: [PATCH 07/12] add bb origin to surfaces. only relevant when using non
 0,0,0 local origin.

---
 LoopStructural/utils/_surface.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/LoopStructural/utils/_surface.py b/LoopStructural/utils/_surface.py
index a6c43111..efbf2103 100644
--- a/LoopStructural/utils/_surface.py
+++ b/LoopStructural/utils/_surface.py
@@ -145,7 +145,7 @@ def fit(
             values = np.zeros(verts.shape[0]) + isovalue
             # need to add both global and local origin. If the bb is a buffer the local
             # origin may not be 0
-            verts += self.bounding_box.global_origin
+            verts += self.bounding_box.global_origin+self.bounding_box.origin
             surfaces.append(
                 Surface(
                     vertices=verts,

From f8a1bb2a045cabff12e311e75bd99e7117e844b6 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@gmail.com>
Date: Mon, 28 Apr 2025 14:15:50 +1000
Subject: [PATCH 08/12] return abutting region when adding abutting region

---
 LoopStructural/modelling/features/fault/_fault_segment.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/LoopStructural/modelling/features/fault/_fault_segment.py b/LoopStructural/modelling/features/fault/_fault_segment.py
index 737b8f8a..e98f945a 100644
--- a/LoopStructural/modelling/features/fault/_fault_segment.py
+++ b/LoopStructural/modelling/features/fault/_fault_segment.py
@@ -467,7 +467,7 @@ def add_abutting_fault(self, abutting_fault_feature, positive=None):
             )
         self.abut[abutting_fault_feature.name] = abutting_region
         self.__getitem__(0).add_region(abutting_region)
-
+        return abutting_region
     def save(self, filename, scalar_field=True, slip_vector=True, surface=True):
         """
         Save the fault to a file

From 921ec09e42fa7f26636000aa7365407301857d19 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@gmail.com>
Date: Mon, 28 Apr 2025 15:17:23 +1000
Subject: [PATCH 09/12] fix: remove default initialisation with mutable.

Python being python makes these default arguments a shared
attribute which is undesirable.
---
 .../modelling/features/_lambda_geological_feature.py  | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/LoopStructural/modelling/features/_lambda_geological_feature.py b/LoopStructural/modelling/features/_lambda_geological_feature.py
index 487a2bf4..d8f0fb3a 100644
--- a/LoopStructural/modelling/features/_lambda_geological_feature.py
+++ b/LoopStructural/modelling/features/_lambda_geological_feature.py
@@ -18,8 +18,8 @@ def __init__(
         name: str = "unnamed_lambda",
         gradient_function: Optional[Callable[[np.ndarray], np.ndarray]] = None,
         model=None,
-        regions: list = [],
-        faults: list = [],
+        regions: Optional[list] = None,
+        faults: Optional[list] = None,
         builder=None,
     ):
         """A lambda geological feature is a wrapper for a geological
@@ -43,10 +43,11 @@ def __init__(
         builder : _type_, optional
             _description_, by default None
         """
-        BaseFeature.__init__(self, name, model, faults, regions, builder)
+        BaseFeature.__init__(self, name, model, faults if faults is not None else [], regions if regions is not None else [], builder)
         self.type = FeatureType.LAMBDA
         self.function = function
         self.gradient_function = gradient_function
+        self.regions = regions if regions is not None else []
 
     def evaluate_value(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray:
         """_summary_
@@ -64,11 +65,11 @@ def evaluate_value(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray:
         v = np.zeros((pos.shape[0]))
         v[:] = np.nan
 
-        # mask = self._calculate_mask(pos, ignore_regions=ignore_regions)
+        mask = self._calculate_mask(pos, ignore_regions=ignore_regions)
         pos = self._apply_faults(pos)
         if self.function is not None:
             
-            v = self.function(pos)
+            v[mask] = self.function(pos[mask,:])
         return v
 
     def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray:

From a866916643896f829cda1275144d25f0c452d658 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@gmail.com>
Date: Mon, 28 Apr 2025 15:22:59 +1000
Subject: [PATCH 10/12] abstract regions so that positive/negative regions
 aren't duplicating logic

---
 LoopStructural/utils/regions.py | 57 +++++++++++++++------------------
 1 file changed, 25 insertions(+), 32 deletions(-)

diff --git a/LoopStructural/utils/regions.py b/LoopStructural/utils/regions.py
index 792f71da..0c6ee543 100644
--- a/LoopStructural/utils/regions.py
+++ b/LoopStructural/utils/regions.py
@@ -1,6 +1,6 @@
 import numpy as np
 from abc import ABC, abstractmethod
-
+from typing import Tuple
 
 class BaseRegion(ABC):
     @abstractmethod
@@ -34,8 +34,7 @@ def __init__(self, function):
     def __call__(self, xyz):
         return self.function(xyz)
 
-
-class PositiveRegion(BaseRegion):
+class BaseSignRegion(BaseRegion):
     """Helper class for evaluating whether you are in the positive  region of a scalar field.
     If its outside of the support it will interpolate the average gradient at a point on the 0 isovalue
     and calculate the distance from this. Alternatively, a point and vector can be used to save computational time
@@ -46,7 +45,7 @@ def __init__(self, feature, vector=None, point=None):
         self.name = 'PositiveRegion'
         self.parent = feature
 
-    def __call__(self, xyz) -> np.ndarray:
+    def _calculate_value_and_distance(self, xyz)-> Tuple[np.ndarray, np.ndarray]:
         val = self.feature.evaluate_value(xyz)
         # find a point on/near 0 isosurface
         if self.point is None:
@@ -63,17 +62,33 @@ def __call__(self, xyz) -> np.ndarray:
             average_gradient /= np.linalg.norm(average_gradient)
         else:
             average_gradient = self.vector
-        # distance = ((xyz[:,0] - centre[None,0])*average_gradient[0] +
-        #             (xyz[:,1] - centre[None,1])*average_gradient[1] +
-        #             ( xyz[:,2] - centre[None,2])*average_gradient[2])
-        distance = np.einsum("ij,j->i", centre[None, :] - xyz, average_gradient)
+
+        distance = np.einsum(
+            "ij,j->i", centre[None, :] - xyz, average_gradient.reshape(-1, 3)[0, :]
+        )
+        return val, distance
+
+
+class PositiveRegion(BaseSignRegion):
+    """Helper class for evaluating whether you are in the positive  region of a scalar field.
+    If its outside of the support it will interpolate the average gradient at a point on the 0 isovalue
+    and calculate the distance from this. Alternatively, a point and vector can be used to save computational time
+    """
+
+    def __init__(self, feature, vector=None, point=None):
+        super().__init__(feature, vector, point)
+        self.name = 'PositiveRegion'
+        self.parent = feature
+
+    def __call__(self, xyz) -> np.ndarray:
+        val, distance = self._calculate_value_and_distance(xyz)
         return np.logical_or(
             np.logical_and(~np.isnan(val), val > 0),
             np.logical_and(np.isnan(val), distance > 0),
         )
 
 
-class NegativeRegion(BaseRegion):
+class NegativeRegion(BaseSignRegion):
     """Helper class for evaluating whether you are in the positive  region of a scalar field.
     If its outside of the support it will interpolate the average gradient at a point on the 0 isovalue
     and calculate the distance from this. Alternatively, a point and vector can be used to save computational time
@@ -85,29 +100,7 @@ def __init__(self, feature, vector=None, point=None):
         self.parent = feature
 
     def __call__(self, xyz) -> np.ndarray:
-        val = self.feature.evaluate_value(xyz)
-        # find a point on/near 0 isosurface
-        if self.point is None:
-            mask = np.zeros(xyz.shape[0], dtype="bool")
-            mask[:] = val < 0
-            if np.sum(mask) == 0:
-                raise ValueError("Cannot find point on surface")
-            centre = xyz[mask, :][0, :]
-        else:
-            centre = self.point
-        if self.vector is None:
-            average_gradient = self.feature.evaluate_gradient(np.array([centre]))[0]
-            average_gradient[2] = 0
-            average_gradient /= np.linalg.norm(average_gradient)
-        else:
-            average_gradient = self.vector
-        distance = np.einsum("ij,j->i", xyz - centre[None, :], average_gradient)
-        # distance = ((xyz[:,0] - centre[None,0])*average_gradient[0] +
-        #             (xyz[:,1] - centre[None,1])*average_gradient[1] +
-        #             ( xyz[:,2] - centre[None,2])*average_gradient[2])
-        # return np.logical_or(np.logical_and(~np.isnan(val),val
-        #      < 0),
-        #      np.logical_and(np.isnan(val),distance>0))
+        val, distance = self._calculate_value_and_distance(xyz)
         return np.logical_or(
             np.logical_and(~np.isnan(val), val < 0),
             np.logical_and(np.isnan(val), distance < 0),

From f2c51049b0a5489b3905dd4a4f0e360170415bb0 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@gmail.com>
Date: Mon, 28 Apr 2025 15:52:10 +1000
Subject: [PATCH 11/12] fix: use fault normals/slip vectors from data if
 available.

Otherwise calculate from trace
---
 .../features/builders/_fault_builder.py       | 50 ++++++++++----
 tests/unit/modelling/test__fault_builder.py   | 67 +++++++++++++++++--
 2 files changed, 97 insertions(+), 20 deletions(-)

diff --git a/LoopStructural/modelling/features/builders/_fault_builder.py b/LoopStructural/modelling/features/builders/_fault_builder.py
index 563f8f5a..b9974586 100644
--- a/LoopStructural/modelling/features/builders/_fault_builder.py
+++ b/LoopStructural/modelling/features/builders/_fault_builder.py
@@ -151,16 +151,30 @@ def create_data_from_geometry(
             ["X", "Y"],
         ].to_numpy()
         if fault_normal_vector is None:
-            # Calculate fault strike using least squares fit
-            X = fault_trace[:, 0].reshape(-1, 1)
-            Y = fault_trace[:, 1]
-            # Fit line Y = mX + b
-            A = np.hstack([X, np.ones_like(X)])
-            m, _ = np.linalg.lstsq(A, Y, rcond=None)[0]
-            # Convert slope to strike vector
-            strike_vector = np.array([1, m, 0])
-            strike_vector /= np.linalg.norm(strike_vector)
-            fault_normal_vector = np.cross(strike_vector, [0, 0, 1])
+            if fault_frame_data.loc[
+            np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["nx"].notna())].shape[0]>0:
+                fault_normal_vector = fault_frame_data.loc[
+                    np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["nx"].notna()),
+                    ["nx", "ny", "nz"],
+                ].to_numpy().mean(axis=0)
+
+            else:
+
+                # Calculate fault strike using eigenvectors
+                pts = fault_trace - fault_trace.mean(axis=0)
+                # Calculate covariance matrix
+                cov_matrix = pts.T @ pts
+                # Get eigenvectors and eigenvalues
+                eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
+                # Use eigenvector with largest eigenvalue as strike direction
+                strike_vector = eigenvectors[:, np.argmax(eigenvalues)]
+                strike_vector = np.append(strike_vector, 0)  # Add z component
+                strike_vector /= np.linalg.norm(strike_vector)
+
+                fault_normal_vector = np.cross(strike_vector, [0, 0, 1])
+                # Rotate the fault normal vector according to the fault dip
+                rotation_matrix = rotation(strike_vector[None, :], np.array([90 - fault_dip]))
+                fault_normal_vector = np.einsum("ijk,ik->ij", rotation_matrix, fault_normal_vector[None, :])[0]
 
         if not isinstance(fault_normal_vector, np.ndarray):
             fault_normal_vector = np.array(fault_normal_vector)
@@ -170,10 +184,18 @@ def create_data_from_geometry(
             fault_slip_vector = np.einsum("ijk,ik->ij", rotation_matrix, fault_normal_vector[None, :])[0]
 
         if fault_slip_vector is None:
-            fault_slip_vector = np.cross(fault_normal_vector, [1., 0., 0.])
-            if np.linalg.norm(fault_slip_vector) == 0:
-                fault_slip_vector = np.cross(fault_normal_vector, [0., 1., 0.])
-            fault_slip_vector /= np.linalg.norm(fault_slip_vector)
+            if fault_frame_data.loc[
+            np.logical_and(fault_frame_data["coord"] == 1, fault_frame_data["nx"].notna())].shape[0]>0:
+                fault_slip_vector = fault_frame_data.loc[
+                    np.logical_and(fault_frame_data["coord"] == 1, fault_frame_data["nx"].notna()),
+                    ["nx", "ny", "nz"],
+                ].to_numpy().mean(axis=0)
+
+            else:
+                fault_slip_vector = np.cross(fault_normal_vector, [1., 0., 0.])
+                if np.linalg.norm(fault_slip_vector) == 0:
+                    fault_slip_vector = np.cross(fault_normal_vector, [0., 1., 0.])
+                fault_slip_vector /= np.linalg.norm(fault_slip_vector)
         if fault_center is None:
             fault_trace = fault_frame_data.loc[
                 np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0),
diff --git a/tests/unit/modelling/test__fault_builder.py b/tests/unit/modelling/test__fault_builder.py
index 9006b61f..f5987a8e 100644
--- a/tests/unit/modelling/test__fault_builder.py
+++ b/tests/unit/modelling/test__fault_builder.py
@@ -21,13 +21,16 @@ def test_fault_builder_update_geometry(interpolatortype):
     fault_builder.update_geometry(points)
 
     # Check if the origin and maximum are updated correctly
-    expected_origin = np.array([0.5, 0.5, 0.5])
-    expected_maximum = np.array([0.7, 0.7, 0.7])
-    assert np.allclose(fault_builder.origin, expected_origin)
-    assert np.allclose(fault_builder.maximum, expected_maximum)
+    # Calculate buffer as 10% of the range
+    buffer = 0.1 * (np.array([0.7, 0.7, 0.7]) - np.array([0.5, 0.5, 0.5]))
+    expected_origin = np.array([0.5, 0.5, 0.5]) - buffer
+    expected_maximum = np.array([0.7, 0.7, 0.7]) + buffer
+    # rtol 1e-2 is 1% relative tolerance 
+    assert np.allclose(fault_builder.origin, expected_origin,rtol=1e-2)
+    assert np.allclose(fault_builder.maximum, expected_maximum, rtol=1e-2)
 
 
-def test_fault_builder_create_data_from_geometry(interpolatortype):
+def test_fault_builder_create_data_from_geometry_with_normals(interpolatortype):
     # Create a FaultBuilder instance
     bounding_box = BoundingBox([0, 0, 0], [1, 1, 1])
     nelements = 1000
@@ -60,7 +63,59 @@ def test_fault_builder_create_data_from_geometry(interpolatortype):
     # Check if the fault attributes are updated correctly
     expected_fault_normal_vector = np.array([0, 1.0, 0])
     expected_fault_slip_vector = np.array([0.0, 0.0, 1.0])
-    # expected_fault_centre = np.array([0.15, 0.6, 1.05])
+    expected_fault_centre = np.array([0.15, 0.6, 1.05])
+    expected_fault_normal_vector /= np.linalg.norm(expected_fault_normal_vector)
+    expected_fault_slip_vector /= np.linalg.norm(expected_fault_slip_vector)
+    # expected_fault_minor_axis = 0.5
+    # expected_fault_major_axis = 1.0
+    # expected_fault_intermediate_axis = 1.0
+    calculated_fault_normal_vector = fault_builder.fault_normal_vector
+    calculated_fault_slip_vector = fault_builder.fault_slip_vector
+    calculated_fault_normal_vector /= np.linalg.norm(calculated_fault_normal_vector)
+    calculated_fault_slip_vector /= np.linalg.norm(calculated_fault_slip_vector)
+    assert np.allclose(calculated_fault_normal_vector, expected_fault_normal_vector)
+    assert np.allclose(calculated_fault_slip_vector, expected_fault_slip_vector)
+    # assert np.allclose(fault_builder.fault_centre, expected_fault_centre)
+    # assert np.isclose(fault_builder.fault_minor_axis, expected_fault_minor_axis)
+    # assert np.isclose(fault_builder.fault_major_axis, expected_fault_major_axis)
+    # assert np.isclose(fault_builder.fault_intermediate_axis, expected_fault_intermediate_axis)
+
+
+def test_fault_builder_create_data_from_geometry_without_normals(interpolatortype):
+    # Create a FaultBuilder instance
+    bounding_box = BoundingBox([0, 0, 0], [1, 1, 1])
+    nelements = 1000
+    model = GeologicalModel([0, 0, 0], [1, 1, 1])
+    fault_bounding_box_buffer = 0.2
+    fault_builder = FaultBuilder(
+        interpolatortype, bounding_box, nelements, model, fault_bounding_box_buffer
+    )
+
+    # Create some test data
+    fault_frame_data = pd.DataFrame(
+        {
+            "coord": [0.0, 0.0, 1.0, 2.0],
+            "val": [0.0, 0.0, 1.0, 1.0],
+            "gx": [np.nan, np.nan, np.nan, np.nan],
+            "gy": [np.nan, np.nan, np.nan, np.nan],
+            "gz": [np.nan, np.nan, np.nan, np.nan],
+            "X": [0.1, 0.2, 0.3, 0.4],
+            "Y": [0.5, 0.6, 0.7, 0.8],
+            "Z": [0.9, 1.0, 1.1, 1.2],
+            "nx": [np.nan, np.nan, np.nan, np.nan],
+            "ny": [np.nan, np.nan, np.nan, np.nan],
+            "nz": [np.nan, np.nan, np.nan, np.nan],
+        }
+    )
+
+    # Call the create_data_from_geometry method
+    fault_builder.create_data_from_geometry(fault_frame_data)
+    strike_vector = np.array([-.1,-.1,0])
+    strike_vector /= np.linalg.norm(strike_vector)
+    expected_fault_normal_vector = np.cross(strike_vector, np.array([0, 0, 1]))
+    expected_fault_normal_vector /= np.linalg.norm(expected_fault_normal_vector)
+    # Check if the fault attributes are updated correctly
+    expected_fault_slip_vector = np.cross(expected_fault_normal_vector, [1.0, 0.0, 0.0])
     expected_fault_normal_vector /= np.linalg.norm(expected_fault_normal_vector)
     expected_fault_slip_vector /= np.linalg.norm(expected_fault_slip_vector)
     # expected_fault_minor_axis = 0.5

From 07ccd5411198018b39291c2dfefbe44225abaa48 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@gmail.com>
Date: Tue, 29 Apr 2025 15:07:31 +1000
Subject: [PATCH 12/12] fix: add faulted vector calc to lambda feature

---
 .../features/_base_geological_feature.py      |  2 +-
 .../features/_lambda_geological_feature.py    | 60 ++++++++++++++++++-
 2 files changed, 58 insertions(+), 4 deletions(-)

diff --git a/LoopStructural/modelling/features/_base_geological_feature.py b/LoopStructural/modelling/features/_base_geological_feature.py
index 39ea44bf..f1217f0b 100644
--- a/LoopStructural/modelling/features/_base_geological_feature.py
+++ b/LoopStructural/modelling/features/_base_geological_feature.py
@@ -229,7 +229,7 @@ def evaluate_gradient(self, pos, ignore_regions=False):
         """
         Evaluate the gradient of the feature at a given position.
         """
-
+    
         raise NotImplementedError
 
     def min(self):
diff --git a/LoopStructural/modelling/features/_lambda_geological_feature.py b/LoopStructural/modelling/features/_lambda_geological_feature.py
index d8f0fb3a..8e5ff331 100644
--- a/LoopStructural/modelling/features/_lambda_geological_feature.py
+++ b/LoopStructural/modelling/features/_lambda_geological_feature.py
@@ -1,12 +1,13 @@
 """
 Geological features
 """
-
+from LoopStructural.utils.maths import regular_tetraherdron_for_points, gradient_from_tetrahedron
 from ...modelling.features import BaseFeature
 from ...utils import getLogger
 from ...modelling.features import FeatureType
 import numpy as np
 from typing import Callable, Optional
+from ...utils import LoopValueError
 
 logger = getLogger(__name__)
 
@@ -68,11 +69,11 @@ def evaluate_value(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray:
         mask = self._calculate_mask(pos, ignore_regions=ignore_regions)
         pos = self._apply_faults(pos)
         if self.function is not None:
-            
+
             v[mask] = self.function(pos[mask,:])
         return v
 
-    def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray:
+    def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False,element_scale_parameter=None) -> np.ndarray:
         """_summary_
 
         Parameters
@@ -85,7 +86,60 @@ def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray
         np.ndarray
             _description_
         """
+        if pos.shape[1] != 3:
+            raise LoopValueError("Need Nx3 array of xyz points to evaluate gradient")
+        logger.info(f'Calculating gradient for {self.name}')
+        if element_scale_parameter is None:
+            if self.model is not None:
+                element_scale_parameter = np.min(self.model.bounding_box.step_vector) / 10
+            else:
+                element_scale_parameter = 1
+        else:
+            try:
+                element_scale_parameter = float(element_scale_parameter)
+            except ValueError:
+                logger.error("element_scale_parameter must be a float")
+                element_scale_parameter = 1
         v = np.zeros((pos.shape[0], 3))
+        v = np.zeros(pos.shape)
+        v[:] = np.nan
+        mask = self._calculate_mask(pos, ignore_regions=ignore_regions)
+        # evaluate the faults on the nodes of the faulted feature support
+        # then evaluate the gradient at these points
+        if len(self.faults) > 0:
+            # generate a regular tetrahedron for each point
+            # we will then move these points by the fault and then recalculate the gradient.
+            # this should work...
+            resolved = False
+            tetrahedron = regular_tetraherdron_for_points(pos, element_scale_parameter)
+
+            while resolved:
+                for f in self.faults:
+                    v = (
+                        f[0]
+                        .evaluate_value(tetrahedron.reshape(-1, 3), fillnan='nearest')
+                        .reshape(tetrahedron.shape[0], 4)
+                    )
+                    flag = np.logical_or(np.all(v > 0, axis=1), np.all(v < 0, axis=1))
+                    if np.any(~flag):
+                        logger.warning(
+                            f"Points are too close to fault {f[0].name}. Refining the tetrahedron"
+                        )
+                        element_scale_parameter *= 0.5
+                        tetrahedron = regular_tetraherdron_for_points(pos, element_scale_parameter)
+
+                resolved = True
+
+            tetrahedron_faulted = self._apply_faults(np.array(tetrahedron.reshape(-1, 3))).reshape(
+                tetrahedron.shape
+            )
+
+            values = self.function(tetrahedron_faulted.reshape(-1, 3)).reshape(
+                (-1, 4)
+            )
+            v[mask, :] = gradient_from_tetrahedron(tetrahedron[mask, :, :], values[mask])
+
+            return v
         if self.gradient_function is None:
             v[:, :] = np.nan
         else: