fix: under constrained faults now work

automatically adds an anisotropy to the fault trace
if no dip data assumes 90 degrees and adds an anisotropy along dip vector
fault_dip parameter can be used to tune the dip value
This should make faults more reliable in loopstructural without having to add lots of points
fault_dip_anisotropy and fault_trace_anisotropy are weights for the least squares system
default to 1.0 but can be increased/decreased to change how much these constraints are used.

Author: lachlangrose

LoopStructural/modelling/core/geological_model.py

@@ -1,6 +1,7 @@
 """
 Main entry point for creating a geological model
 """
+
 from ...utils import getLogger, log_to_file
 
 import numpy as np
@@ -134,7 +135,6 @@ def __init__(
         )
         logger.info(maximumstr)
 
-        lengths = self.maximum - self.origin
         self.scale_factor = 1.0
 
         self.bounding_box = BoundingBox(
@@ -310,7 +310,6 @@ def from_processor(cls, processor):
         model.data = processor.data
         if processor.fault_properties is not None:
             for i in processor.fault_network.faults:
-                logger.info(f"Adding fault {i}")
                 model.create_and_add_fault(
                     i,
                     **processor.fault_properties.to_dict("index")[i],
@@ -642,27 +641,6 @@ def set_model_data(self, data):
         )
         self.data = data
 
-    def extend_model_data(self, newdata):
-        """
-        Extends the data frame
-
-        Parameters
-        ----------
-        newdata : pandas data frame
-            data to add to the existing dataframe
-        Returns
-        -------
-        """
-        logger.warning("Extend data is untested and may have unexpected consequences")
-        data_temp = newdata.copy()
-        data_temp["X"] -= self.origin[0]
-        data_temp["Y"] -= self.origin[1]
-        data_temp["Z"] -= self.origin[2]
-        data_temp["X"] /= self.scale_factor
-        data_temp["Y"] /= self.scale_factor
-        data_temp["Z"] /= self.scale_factor
-        self.data.concat([self.data, data_temp], sort=True)
-
     def set_stratigraphic_column(self, stratigraphic_column, cmap="tab20"):
         """
         Adds a stratigraphic column to the model
@@ -833,7 +811,7 @@ def create_and_add_folded_foliation(
         fold_frame=None,
         svario=True,
         tol=None,
-        invert_fold_norm = False,
+        invert_fold_norm=False,
         **kwargs,
     ):
         """
@@ -873,8 +851,10 @@ def create_and_add_folded_foliation(
             fold_frame = self.features[-1]
         assert type(fold_frame) == FoldFrame, "Please specify a FoldFrame"
 
-        fold = FoldEvent(fold_frame, name=f"Fold_{foliation_data}", invert_norm=invert_fold_norm)
-        
+        fold = FoldEvent(
+            fold_frame, name=f"Fold_{foliation_data}", invert_norm=invert_fold_norm
+        )
+
         if "fold_weights" not in kwargs:
             kwargs["fold_weights"] = {}
         if interpolatortype != "DFI":
@@ -1342,6 +1322,9 @@ def create_and_add_fault(
         force_mesh_geometry: bool = False,
         points: bool = False,
         fault_buffer=0.2,
+        fault_trace_anisotropy=1.0,
+        fault_dip=90,
+        fault_dip_anisotropy=1.0,
         **kwargs,
     ):
         """
@@ -1443,6 +1426,9 @@ def create_and_add_fault(
             points=points,
             force_mesh_geometry=force_mesh_geometry,
             fault_buffer=fault_buffer,
+            fault_trace_anisotropy=fault_trace_anisotropy,
+            fault_dip=fault_dip,
+            fault_dip_anisotropy=fault_dip_anisotropy,
         )
         if "force_mesh_geometry" not in kwargs:
             fault_frame_builder.set_mesh_geometry(kwargs.get("fault_buffer", 0.2), 0)

LoopStructural/modelling/features/builders/_fault_builder.py

@@ -1,5 +1,6 @@
 from typing import Union
 from ._structural_frame_builder import StructuralFrameBuilder
+from .. import AnalyticalGeologicalFeature
 from LoopStructural.utils import get_vectors
 import numpy as np
 import pandas as pd
@@ -90,6 +91,9 @@ def create_data_from_geometry(
         points=False,
         force_mesh_geometry=False,
         fault_buffer=0.2,
+        fault_trace_anisotropy=1.0,
+        fault_dip=90,
+        fault_dip_anisotropy=1.0,
     ):
         """Generate the required data for building a fault frame for a fault with the
         specified parameters
@@ -140,12 +144,34 @@ def create_data_from_geometry(
                     fault_frame_data.loc[normal_mask, ["nx", "ny", "nz"]].to_numpy(),
                 ]
             )
+
             if len(vector_data) == 0:
                 logger.error(
                     "You cannot model a fault without defining the orientation of the fault\n\
-                    Either add orientation data or define the fault normal vector argument"
+                    Defaulting to a vertical fault"
+                )
+                coefficients = np.polyfit(
+                    fault_frame_data.loc[trace_mask, "X"],
+                    fault_frame_data.loc[trace_mask, "Y"],
+                    1,
                 )
-                raise ValueError(f"There are no points on the fault trace")
+                slope, intercept = coefficients
+
+                # Create a direction vector using the slope
+                direction_vector = np.array([1, slope])
+                direction_vector /= np.linalg.norm(direction_vector)
+                print(f"Fault dip: {fault_dip}")
+                vector_data = np.array(
+                    [
+                        [
+                            direction_vector[1],
+                            -direction_vector[0],
+                            np.sin(np.deg2rad(fault_dip)),
+                        ]
+                    ]
+                )
+                vector_data /= np.linalg.norm(vector_data, axis=1)[:, None]
+
             fault_normal_vector = np.mean(vector_data, axis=0)
 
         logger.info(f"Fault normal vector: {fault_normal_vector}")
@@ -283,9 +309,9 @@ def create_data_from_geometry(
                         fault_frame_data["coord"] == 0,
                         ~np.isnan(fault_frame_data["nx"]),
                     )
-                    fault_frame_data.loc[
-                        mask, ["gx", "gy", "gz"]
-                    ] = fault_frame_data.loc[mask, ["nx", "ny", "nz"]]
+                    fault_frame_data.loc[mask, ["gx", "gy", "gz"]] = (
+                        fault_frame_data.loc[mask, ["nx", "ny", "nz"]]
+                    )
 
                     fault_frame_data.loc[mask, ["nx", "ny", "nz"]] = np.nan
                     mask = np.logical_and(
@@ -411,7 +437,37 @@ def create_data_from_geometry(
                     1,
                     w,
                 ]
+                fault_frame_data.loc[
+                    len(fault_frame_data),
+                    [
+                        "X",
+                        "Y",
+                        "Z",
+                        "feature_name",
+                        "nx",
+                        "ny",
+                        "nz",
+                        "val",
+                        "coord",
+                        "w",
+                    ],
+                ] = [
+                    fault_center[0],
+                    fault_center[1],
+                    fault_center[2],
+                    self.name,
+                    fault_normal_vector[0],
+                    fault_normal_vector[1],
+                    fault_normal_vector[2],
+                    np.nan,
+                    0,
+                    w,
+                ]
         self.add_data_from_data_frame(fault_frame_data)
+        if fault_trace_anisotropy > 0:
+            self.add_fault_trace_anisotropy(fault_trace_anisotropy)
+        if fault_dip_anisotropy > 0:
+            self.add_fault_dip_anisotropy(fault_dip_anisotropy)
         if force_mesh_geometry:
             self.set_mesh_geometry(fault_buffer, None)
         self.update_geometry(fault_frame_data[["X", "Y", "Z"]].to_numpy())
@@ -425,7 +481,6 @@ def set_mesh_geometry(self, buffer, rotation):
             percentage of length to add to edges
         """
         length = np.max(self.maximum - self.origin)
-        print(self.origin, length * buffer, self.maximum)
         # for builder in self.builders:
         # all three coordinates share the same support
         self.builders[0].set_interpolation_geometry(
@@ -464,6 +519,50 @@ def splayregion(xyz):
         self.builders[0].add_equality_constraints(splay, splayregion, scalefactor)
         return splayregion
 
+    def add_fault_trace_anisotropy(self, w=1.0):
+        trace_data = self.builders[0].data.loc[self.builders[0].data["val"] == 0, :]
+        coefficients = np.polyfit(
+            trace_data["X"],
+            trace_data["Y"],
+            1,
+        )
+        slope, intercept = coefficients
+
+        # Create a direction vector using the slope
+        direction_vector = np.array([1, slope])
+        direction_vector /= np.linalg.norm(direction_vector)
+        vector_data = np.array([[direction_vector[1], -direction_vector[0], 0]])
+        anisotropy_feature = AnalyticalGeologicalFeature(
+            vector=vector_data, origin=[0, 0, 0], name="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, dip, w=1.0):
+        trace_data = self.builders[0].data.loc[self.builders[0].data["val"] == 0, :]
+        coefficients = np.polyfit(
+            trace_data["X"],
+            trace_data["Y"],
+            1,
+        )
+        slope, intercept = coefficients
+
+        # Create a direction vector using the slope
+        direction_vector = np.array([1, slope])
+        direction_vector /= np.linalg.norm(direction_vector)
+        vector_data = np.array([[direction_vector[0], direction_vector[1], 0]])
+
+        vector_data = np.array(
+            [[direction_vector[1], -direction_vector[0], np.sin(np.deg2rad(dip))]]
+        )
+        anisotropy_feature = AnalyticalGeologicalFeature(
+            vector=vector_data, origin=[0, 0, 0], name="fault_dip_anisotropy"
+        )
+        self.builders[0].add_orthogonal_feature(
+            anisotropy_feature, w=w, region=None, step=1, B=0
+        )
+
     def update(self):
         for i in range(3):
             self.builders[i].update()