Fix/fault displacements

LoopStructural/datatypes/_bounding_box.py

@@ -387,7 +387,7 @@ def regular_grid(
 
         if not local:
             coordinates = [
-                np.linspace(self.global_origin[i], self.global_maximum[i], nsteps[i])
+                np.linspace(self.global_origin[i]+self.origin[i], self.global_maximum[i], nsteps[i])
                 for i in range(self.dimensions)
             ]
         coordinate_grid = np.meshgrid(*coordinates, indexing="ij")

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__)
 
@@ -18,8 +19,8 @@ def __init__(
         name: str = "unnamed_lambda",
         gradient_function: Optional[Callable[[np.ndarray], np.ndarray]] = None,
         model=None,
-        regions: list = [],
-        faults: list = [],
+        regions: Optional[list] = None,
+        faults: Optional[list] = None,
         builder=None,
     ):
         """A lambda geological feature is a wrapper for a geological
@@ -43,10 +44,11 @@ def __init__(
         builder : _type_, optional
             _description_, by default None
         """
-        BaseFeature.__init__(self, name, model, faults, regions, builder)
+        BaseFeature.__init__(self, name, model, faults if faults is not None else [], regions if regions is not None else [], builder)
         self.type = FeatureType.LAMBDA
         self.function = function
         self.gradient_function = gradient_function
+        self.regions = regions if regions is not None else []
 
     def evaluate_value(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray:
         """_summary_
@@ -62,13 +64,16 @@ def evaluate_value(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray:
             _description_
         """
         v = np.zeros((pos.shape[0]))
-        if self.function is None:
-            v[:] = np.nan
-        else:
-            v[:] = self.function(pos)
+        v[:] = np.nan
+
+        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
@@ -81,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: