fix: return axis and angle for movement of fault points

Author: lachlangrose

LoopStructural/modelling/features/_geological_feature.py

@@ -1,13 +1,14 @@
 """
 Geological features
 """
+
 from ...modelling.features import BaseFeature
 from ...utils import getLogger
 from ...modelling.features import FeatureType
 from ...interpolators import GeologicalInterpolator
 import numpy as np
 
-from ...utils import getLogger, LoopValueError
+from ...utils import getLogger, LoopValueError, rotate
 
 logger = getLogger(__name__)
 
@@ -114,14 +115,14 @@ def evaluate_value(self, evaluation_points: np.ndarray) -> np.ndarray:
         v = np.zeros(evaluation_points.shape[0])
         v[:] = np.nan
         mask = self._calculate_mask(evaluation_points)
-        evaluation_points = self._apply_faults(evaluation_points)
+        evaluation_points, axis, angle = self._apply_faults(evaluation_points)
         if mask.dtype not in [int, bool]:
             logger.error(f"Unable to evaluate value for {self.name}")
         else:
             v[mask] = self.interpolator.evaluate_value(evaluation_points[mask, :])
         return v
 
-    def evaluate_gradient(self, evaluation_points: np.ndarray) -> np.ndarray:
+    def evaluate_gradient(self, pos: np.ndarray) -> np.ndarray:
         """
 
         Parameters
@@ -134,18 +135,24 @@ def evaluate_gradient(self, evaluation_points: np.ndarray) -> np.ndarray:
 
 
         """
-        if evaluation_points.shape[1] != 3:
+        if pos.shape[1] != 3:
             raise LoopValueError("Need Nx3 array of xyz points to evaluate gradient")
         self.builder.up_to_date()
-        v = np.zeros(evaluation_points.shape)
+        v = np.zeros(pos.shape)
         v[:] = np.nan
-        mask = self._calculate_mask(evaluation_points)
-        evaluation_points = self._apply_faults(evaluation_points)
+        mask = self._calculate_mask(pos)
+        original_pos = pos.copy()
+        pos, axis, angle = self._apply_faults(pos)
         if mask.dtype not in [int, bool]:
             logger.error(f"Unable to evaluate gradient for {self.name}")
         else:
-            v[mask, :] = self.interpolator.evaluate_gradient(evaluation_points[mask, :])
-
+            v[mask, :] = self.interpolator.evaluate_gradient(pos[mask, :])
+            v_norm = np.linalg.norm(v[mask, :], axis=1)
+            v[mask, :] /= v_norm[:, None]
+            if axis is not None:
+                for i in reversed(range(axis.shape[1])):
+                    v[mask, :] = rotate(v[mask, :], axis[:, i, :], angle[:, i])
+            v[mask, :] *= v_norm[:, None]
         return v
 
     def evaluate_gradient_misfit(self):

LoopStructural/modelling/features/fault/_fault_segment.py

@@ -19,9 +19,7 @@ class FaultSegment(StructuralFrame):
     Class for representing a slip event of a fault
     """
 
-    def __init__(
-        self, features, name, faultfunction=None, steps=3, displacement=1.0, fold=None
-    ):
+    def __init__(self, features, name, faultfunction=None, steps=10, displacement=1.0, fold=None):
         """
         A slip event of a fault
 
@@ -263,7 +261,7 @@ def evaluate_displacement(self, points):
             d[mask] = self.faultfunction(gx[mask], gy[mask], gz[mask])
         return d * self.displacement
 
-    def apply_to_points(self, points):
+    def apply_to_points(self, points, reverse=False):
         """
         Unfault the array of points
 
@@ -308,6 +306,13 @@ def apply_to_points(self, points):
         mask = np.abs(d) > 0.0
 
         d *= self.displacement
+        if reverse:
+            d *= -1.0
+        # calculate the anglee and vector to rotate the points
+        axis = np.zeros((points.shape[0], steps, 3))
+        axis[:] = np.nan
+        angle = np.zeros((points.shape[0], steps))
+        angle[:] = np.nan
         # calculate the fault frame for the evaluation points
         for i in range(steps):
             gx = None
@@ -318,18 +323,10 @@ def apply_to_points(self, points):
                 with ThreadPoolExecutor(max_workers=8) as executor:
                     # all of these operations should be
                     # independent so just run as different threads
-                    gx_future = executor.submit(
-                        self.__getitem__(0).evaluate_value, newp[mask, :]
-                    )
-                    g_future = executor.submit(
-                        self.__getitem__(1).evaluate_gradient, newp[mask, :]
-                    )
-                    gy_future = executor.submit(
-                        self.__getitem__(1).evaluate_value, newp[mask, :]
-                    )
-                    gz_future = executor.submit(
-                        self.__getitem__(2).evaluate_value, newp[mask, :]
-                    )
+                    gx_future = executor.submit(self.__getitem__(0).evaluate_value, newp[mask, :])
+                    g_future = executor.submit(self.__getitem__(1).evaluate_gradient, newp[mask, :])
+                    gy_future = executor.submit(self.__getitem__(1).evaluate_value, newp[mask, :])
+                    gz_future = executor.submit(self.__getitem__(2).evaluate_value, newp[mask, :])
                     gx = gx_future.result()
                     g = g_future.result()
                     gy = gy_future.result()
@@ -363,10 +360,18 @@ def apply_to_points(self, points):
             # multiply displacement vector by the displacement magnitude for
             # step
             g *= (1.0 / steps) * d[:, None]
-
+            prev_p = newp[mask, :].copy()
             # apply displacement
             newp[mask, :] += g
-        return newp
+            # axis[mask, i, :] = np.cross(prev_p, newp[mask, :], axisa=1, axisb=1)
+            # angle[mask, i] = 2 * np.arctan(
+            #     0.5 * (np.linalg.norm(prev_p - newp[mask, :], axis=1))
+            # )
+            # calculate the angle between the previous and new points
+            # and the axis of rotation
+            # g /= np.linalg.norm(g, axis=1)[:, None]
+            axis[mask, i, :] /= np.linalg.norm(axis[mask, i, :], axis=1)[:, None]
+        return newp, axis, angle
 
     def add_abutting_fault(self, abutting_fault_feature, positive=None):