fix: gradient of faulted feature will not be aliased by the interpola…

…tion grid

LoopStructural/modelling/features/_geological_feature.py

@@ -2,6 +2,7 @@
 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
@@ -131,7 +132,9 @@ def evaluate_value(self, pos: np.ndarray, ignore_regions=False, fillnan=None) ->
             v[nanmask] = v[~nanmask][i]
         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:
         """
 
         Parameters
@@ -147,6 +150,17 @@ def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray
         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
 
         self.builder.up_to_date()
 
@@ -156,16 +170,38 @@ def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray
         # 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.interpolator.evaluate_value(tetrahedron_faulted.reshape(-1, 3)).reshape(
+                (-1, 4)
+            )
+            v[mask, :] = gradient_from_tetrahedron(tetrahedron_faulted[mask, :, :], values[mask])
 
-            if issubclass(type(self.interpolator), DiscreteInterpolator):
-                points = self.interpolator.support.nodes
-            else:
-                raise NotImplementedError(
-                    "Faulted feature gradients are only supported by DiscreteInterpolator at the moment."
-                )
-            points_faulted = self._apply_faults(points)
-            values = self.interpolator.evaluate_value(points_faulted)
-            v[mask, :] = self.interpolator.support.evaluate_gradient(pos[mask, :], values)
             return v
         pos = self._apply_faults(pos)
         if mask.dtype not in [int, bool]:

LoopStructural/modelling/features/fault/_fault_segment.py

@@ -1,3 +1,4 @@
+from LoopStructural.utils.maths import regular_tetraherdron_for_points, gradient_from_tetrahedron
 from ....modelling.features.fault._fault_function_feature import (
     FaultDisplacementFeature,
 )
@@ -390,51 +391,12 @@ def apply_to_vectors(self, vector: np.ndarray, scale_parameter: float = 1.0) ->
         # the nodes of the tetrahedron are then restored by the fault and then gradient is
         # recalculated using the updated node positions but the original corner values
 
-        regular_tetrahedron = np.array(
-            [
-                [np.sqrt(8 / 9), 0, -1 / 3],
-                [-np.sqrt(2 / 9), np.sqrt(2 / 3), -1 / 3],
-                [-np.sqrt(2 / 9), -np.sqrt(2 / 3), -1 / 3],
-                [0, 0, 1],
-            ]
-        )
-        regular_tetrahedron *= scale_parameter
-        xyz = vector[:, :3]
-        tetrahedron = np.zeros((xyz.shape[0], 4, 3))
-        tetrahedron[:] = xyz[:, None, :]
-        tetrahedron[:, :, :] += regular_tetrahedron[None, :, :]
-
-        vectors = vector[:, 3:]
-        corners = np.einsum('ikj,ij->ik', tetrahedron - xyz[:, None, :], vectors)
-        tetrahedron = tetrahedron.reshape(-1, 3)
-        tetrahedron = self.apply_to_points(tetrahedron)
-        tetrahedron = tetrahedron.reshape(-1, 4, 3)
-        m = np.array(
-            [
-                [
-                    (tetrahedron[:, 1, 0] - tetrahedron[:, 0, 0]),
-                    (tetrahedron[:, 1, 1] - tetrahedron[:, 0, 1]),
-                    (tetrahedron[:, 1, 2] - tetrahedron[:, 0, 2]),
-                ],
-                [
-                    (tetrahedron[:, 2, 0] - tetrahedron[:, 0, 0]),
-                    (tetrahedron[:, 2, 1] - tetrahedron[:, 0, 1]),
-                    (tetrahedron[:, 2, 2] - tetrahedron[:, 0, 2]),
-                ],
-                [
-                    (tetrahedron[:, 3, 0] - tetrahedron[:, 0, 0]),
-                    (tetrahedron[:, 3, 1] - tetrahedron[:, 0, 1]),
-                    (tetrahedron[:, 3, 2] - tetrahedron[:, 0, 2]),
-                ],
-            ]
-        )
-        I = np.array([[-1.0, 1.0, 0.0, 0.0], [-1.0, 0.0, 1.0, 0.0], [-1.0, 0.0, 0.0, 1.0]])
-        m = np.swapaxes(m, 0, 2)
-        element_gradients = np.linalg.inv(m)
-
-        element_gradients = element_gradients.swapaxes(1, 2)
-        element_gradients = element_gradients @ I
-        v = np.sum(element_gradients * corners[:, None, :], axis=2)
+        tetrahedron = regular_tetraherdron_for_points(vector[:, :3], scale_parameter)
+        corners = np.einsum('ikj,ij->ik', tetrahedron - vector[:, None, :3], vector[:, 3:])
+
+        tetrahedron = self.apply_to_points(tetrahedron.reshape(-1, 3)).reshape(-1, 4, 3)
+        v = gradient_from_tetrahedron(tetrahedron, corners)
+
         return v
 
     def add_abutting_fault(self, abutting_fault_feature, positive=None):

LoopStructural/utils/maths.py

@@ -244,3 +244,54 @@ def get_dip_vector(strike, dip):
         ]
     )
     return v
+
+
+def regular_tetraherdron_for_points(xyz, scale_parameter):
+    regular_tetrahedron = np.array(
+        [
+            [np.sqrt(8 / 9), 0, -1 / 3],
+            [-np.sqrt(2 / 9), np.sqrt(2 / 3), -1 / 3],
+            [-np.sqrt(2 / 9), -np.sqrt(2 / 3), -1 / 3],
+            [0, 0, 1],
+        ]
+    )
+    regular_tetrahedron *= scale_parameter
+    tetrahedron = np.zeros((xyz.shape[0], 4, 3))
+    tetrahedron[:] = xyz[:, None, :]
+    tetrahedron[:, :, :] += regular_tetrahedron[None, :, :]
+
+    return tetrahedron
+
+
+def gradient_from_tetrahedron(tetrahedron, value):
+    """
+    Calculate the gradient from a tetrahedron
+    """
+    tetrahedron = tetrahedron.reshape(-1, 4, 3)
+    m = np.array(
+        [
+            [
+                (tetrahedron[:, 1, 0] - tetrahedron[:, 0, 0]),
+                (tetrahedron[:, 1, 1] - tetrahedron[:, 0, 1]),
+                (tetrahedron[:, 1, 2] - tetrahedron[:, 0, 2]),
+            ],
+            [
+                (tetrahedron[:, 2, 0] - tetrahedron[:, 0, 0]),
+                (tetrahedron[:, 2, 1] - tetrahedron[:, 0, 1]),
+                (tetrahedron[:, 2, 2] - tetrahedron[:, 0, 2]),
+            ],
+            [
+                (tetrahedron[:, 3, 0] - tetrahedron[:, 0, 0]),
+                (tetrahedron[:, 3, 1] - tetrahedron[:, 0, 1]),
+                (tetrahedron[:, 3, 2] - tetrahedron[:, 0, 2]),
+            ],
+        ]
+    )
+    I = np.array([[-1.0, 1.0, 0.0, 0.0], [-1.0, 0.0, 1.0, 0.0], [-1.0, 0.0, 0.0, 1.0]])
+    m = np.swapaxes(m, 0, 2)
+    element_gradients = np.linalg.inv(m)
+
+    element_gradients = element_gradients.swapaxes(1, 2)
+    element_gradients = element_gradients @ I
+    v = np.sum(element_gradients * value[:, None, :], axis=2)
+    return v