Merge branch 'fault-anisotropy-fix' into inequality-refresh

Author: lachlangrose

LoopStructural/interpolators/_p1interpolator.py

@@ -175,7 +175,12 @@ def setup_interpolator(self, **kwargs):
         # self.add_interface_constraints(self.interpolation_weights["ipw"])
 
     def add_gradient_orthogonal_constraints(
-        self, points: np.ndarray, vector: np.ndarray, w: float = 1.0, B: float = 0
+        self,
+        points: np.ndarray,
+        vector: np.ndarray,
+        w: float = 1.0,
+        B: float = 0,
+        name='undefined gradient orthogonal constraint',
     ):
         """
         constraints scalar field to be orthogonal to a given vector

LoopStructural/modelling/core/geological_model.py

@@ -1470,7 +1470,6 @@ def regular_grid(self, nsteps=None, shuffle=True, rescale=False, order="C"):
             locations of points in regular grid
         """
         return self.bounding_box.regular_grid(nsteps=nsteps, shuffle=shuffle, order=order)
-        
 
     def evaluate_model(self, xyz: np.ndarray, scale: bool = True) -> np.ndarray:
         """Evaluate the stratigraphic id at each location

LoopStructural/modelling/features/builders/_fault_builder.py

@@ -1,4 +1,6 @@
 from typing import Union
+
+from LoopStructural.utils.maths import rotation
 from ._structural_frame_builder import StructuralFrameBuilder
 from .. import AnalyticalGeologicalFeature
 from LoopStructural.utils import get_vectors
@@ -135,9 +137,9 @@ def create_data_from_geometry(
             )
 
             if len(vector_data) == 0:
-                logger.error(
-                    "You cannot model a fault without defining the orientation of the fault\n\
-                    Defaulting to a vertical fault"
+                logger.warning(
+                    f"No orientation data for fault\n\
+                    Defaulting to a dip of {fault_dip}vertical fault"
                 )
                 if fault_frame_data.loc[trace_mask, :].shape[0] > 3:
 
@@ -149,20 +151,22 @@ def create_data_from_geometry(
                     slope, intercept = coefficients
 
                     # Create a direction vector using the slope
-                    direction_vector = np.array([1, slope])
+                    direction_vector = np.array([1, slope, 0])
                     direction_vector /= np.linalg.norm(direction_vector)
-                    print(f"Fault dip: {fault_dip}")
+                    rotation_matrix = rotation(direction_vector[None, :], [90 - fault_dip])
                     vector_data = np.array(
                         [
                             [
                                 direction_vector[1],
                                 -direction_vector[0],
-                                np.sin(np.deg2rad(fault_dip)),
+                                0,
                             ]
                         ]
                     )
-                    vector_data /= np.linalg.norm(vector_data, axis=1)[:, None]
+                    vector_data /= np.linalg.norm(vector_data, axis=1)
+                    vector_data = np.einsum("ijk,ik->ij", rotation_matrix, vector_data)
 
+                    vector_data /= np.linalg.norm(vector_data, axis=1)
             fault_normal_vector = np.mean(vector_data, axis=0)
 
         logger.info(f"Fault normal vector: {fault_normal_vector}")
@@ -221,11 +225,11 @@ def create_data_from_geometry(
         self.fault_minor_axis = minor_axis
         self.fault_major_axis = major_axis
         self.fault_intermediate_axis = intermediate_axis
-        fault_normal_vector /= np.linalg.norm(fault_normal_vector)
+        self.fault_normal_vector /= np.linalg.norm(self.fault_normal_vector)
         fault_slip_vector /= np.linalg.norm(fault_slip_vector)
         # check if slip vector is inside fault plane, if not project onto fault plane
         # if not np.isclose(normal_vector @ slip_vector, 0):
-        strike_vector = np.cross(fault_normal_vector, fault_slip_vector)
+        strike_vector = np.cross(self.fault_normal_vector, fault_slip_vector)
         self.fault_strike_vector = strike_vector
 
         fault_edges = np.zeros((2, 3))
@@ -252,8 +256,8 @@ def create_data_from_geometry(
             )
         if fault_center is not None:
             if minor_axis is not None:
-                fault_edges[0, :] = fault_center[:3] + fault_normal_vector * minor_axis
-                fault_edges[1, :] = fault_center[:3] - fault_normal_vector * minor_axis
+                fault_edges[0, :] = fault_center[:3] + self.fault_normal_vector * minor_axis
+                fault_edges[1, :] = fault_center[:3] - self.fault_normal_vector * minor_axis
                 self.update_geometry(fault_edges)
 
                 # choose whether to add points -1,1 to constrain fault frame or a scaled
@@ -342,9 +346,9 @@ def create_data_from_geometry(
                             fault_center[1],
                             fault_center[2],
                             self.name,
-                            fault_normal_vector[0] / minor_axis * 0.5,
-                            fault_normal_vector[1] / minor_axis * 0.5,
-                            fault_normal_vector[2] / minor_axis * 0.5,
+                            self.fault_normal_vector[0] / minor_axis * 0.5,
+                            self.fault_normal_vector[1] / minor_axis * 0.5,
+                            self.fault_normal_vector[2] / minor_axis * 0.5,
                             np.nan,
                             0,
                             w,
@@ -439,9 +443,9 @@ def create_data_from_geometry(
 
         self.add_data_from_data_frame(fault_frame_data)
         if fault_trace_anisotropy > 0:
-            self.add_fault_trace_anisotropy(fault_trace_anisotropy)
+            self.add_fault_trace_anisotropy(w=fault_trace_anisotropy)
         if fault_dip_anisotropy > 0:
-            self.add_fault_dip_anisotropy(fault_dip_anisotropy)
+            self.add_fault_dip_anisotropy(w=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())
@@ -501,29 +505,24 @@ def add_fault_trace_anisotropy(self, w: float = 1.0):
         w : float, optional
             _description_, by default 1.0
         """
-        trace_data = self.builders[0].data.loc[self.builders[0].data["val"] == 0, :]
-        if trace_data.shape[0] < 2:
-            logger.warning(
-                f"Only {trace_data.shape[0]} point on fault trace, cannot add fault trace anisotropy. Need at least 3 points"
+        if w > 0:
+
+            plane = np.array([0, 0, 1])
+            strike_vector = (
+                self.fault_normal_vector - np.dot(self.fault_normal_vector, plane) * plane
+            )
+            strike_vector /= np.linalg.norm(strike_vector)
+            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"
+            )
+            print('adding fault trace anisotropy')
+            self.builders[0].add_orthogonal_feature(
+                anisotropy_feature, w=w, region=None, step=1, B=0
             )
-            return
-        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: np.ndarray, w: float = 1.0):
+    def add_fault_dip_anisotropy(self, w: float = 1.0):
         """_summary_
 
         Parameters
@@ -533,32 +532,23 @@ def add_fault_dip_anisotropy(self, dip: np.ndarray, w: float = 1.0):
         w : float, optional
             _description_, by default 1.0
         """
-
-        trace_data = self.builders[0].data.loc[self.builders[0].data["val"] == 0, :]
-        if trace_data.shape[0] < 2:
-            logger.warning(
-                f'Only {trace_data.shape[0]} point on fault trace, cannot add fault trace anisotropy. Need at least 3 points'
+        if w > 0:
+            plane = np.array([0, 0, 1])
+            strike_vector = (
+                self.fault_normal_vector - np.dot(self.fault_normal_vector, plane) * plane
             )
-            return
-        coefficients = np.polyfit(
-            trace_data["X"],
-            trace_data["Y"],
-            1,
-        )
-        slope, intercept = coefficients
+            strike_vector /= np.linalg.norm(strike_vector)
+            strike_vector = np.array([strike_vector[1], -strike_vector[0], 0])
 
-        # 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]])
+            dip_vector = np.cross(strike_vector, self.fault_normal_vector)
 
-        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)
+            anisotropy_feature = AnalyticalGeologicalFeature(
+                vector=dip_vector, origin=[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
+            )
 
     def update(self):
         for i in range(3):

LoopStructural/modelling/features/builders/_geological_feature_builder.py

@@ -304,6 +304,7 @@ def install_gradient_constraint(self):
                     vector[element_idx[::step], :],
                     w=w,
                     B=B,
+                    name=f'{feature.name}_orthogonal',
                 )
 
     def add_equality_constraints(self, feature, region, scalefactor=1.0):
@@ -506,9 +507,9 @@ def build(self, fold=None, fold_weights={}, data_region=None, **kwargs):
                 self.interpolator.support.maximum = newmax
             self.interpolator.set_region(region=region)
 
+        self.interpolator.setup_interpolator(**kwargs)
         self.install_gradient_constraint()
         self.install_equality_constraints()
-        self.interpolator.setup_interpolator(**kwargs)
         self.interpolator.solve_system(**kwargs)
         self._up_to_date = True
         return self._feature

LoopStructural/modelling/features/fault/_fault_function_feature.py

@@ -1,5 +1,5 @@
 from ....modelling.features import BaseFeature, StructuralFrame
-
+from typing import Optional
 from ....utils import getLogger
 
 logger = getLogger(__name__)
@@ -77,3 +77,6 @@ def evaluate_on_surface(self, location):
         fault_strike = self.fault_frame.features[2].evaluate_value(location)
         d = self.displacement.evaluate(fault_displacement, fault_strike)
         return d
+
+    def get_data(self, value_map: Optional[dict] = None):
+        pass

LoopStructural/modelling/features/fault/_fault_segment.py

@@ -423,7 +423,7 @@ def apply_to_vectors(self, vector: np.ndarray, scale_parameter: float = 1.0) ->
 
         element_gradients = element_gradients.swapaxes(1, 2)
         element_gradients = element_gradients @ I
-        v = np.sum(element_gradients * corners, axis=2)
+        v = np.sum(element_gradients * corners[:, None, :], axis=2)
         return v
 
     def add_abutting_fault(self, abutting_fault_feature, positive=None):