fix: add faulted vector calc to lambda feature

Author: Lachlan Grose

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):

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: