fix: add isosurfacing and scalar field method to base feature

adds methods to get the isosurface or scalar field of
a geological feature. These return lists and can provide a list of names
to name the surfaces e.g. get names of units from stratigraphic column
or name using the fault name.

Author: lachlangrose

LoopStructural/interpolators/supports/_3d_structured_grid.py

@@ -301,7 +301,7 @@ def evaluate_value(self, evaluation_points, property_array):
 
         return np.sum(v, axis=1)
 
-    def evaluate_gradient(self, evaluation_points, property_array):
+    def evaluate_gradient(self, evaluation_points, property_array) -> np.ndarray:
         """Evaluate the gradient at a location given node values
 
         Parameters

LoopStructural/modelling/features/_base_geological_feature.py

@@ -1,7 +1,11 @@
+from __future__ import annotations
+
 from abc import ABCMeta, abstractmethod
+from typing import Union, List, Optional
 from LoopStructural.modelling.features import FeatureType
 from LoopStructural.utils import getLogger
 from LoopStructural.utils.typing import NumericInput
+from LoopStructural.api import LoopIsosurfacer, surface_list
 
 
 import numpy as np
@@ -263,3 +267,54 @@ def __tojson__(self):
             "regions": regions,
             "faults": faults,
         }
+
+    def surfaces(
+        self,
+        value: Union[float, int, List[Union[float, int]]],
+        bounding_box=None,
+        name=Optional[Union[List[str], str]],
+    ) -> surface_list:
+        """Find the surfaces of the geological feature at a given value
+
+        Parameters
+        ----------
+        value : Union[float, int, List[float, int]]
+            value or list of values to find the surface of the feature
+
+        Returns
+        -------
+        list
+            list of surfaces
+        """
+        if isinstance(value, (float, int)):
+            value = [value]
+        if bounding_box is None:
+            if self.model is None:
+                raise ValueError("Must specify bounding box")
+            bounding_box = self.model.bounding_box
+
+        isosurfacer = LoopIsosurfacer(bounding_box, self)
+        return isosurfacer.fit(value, name)
+
+    def scalar_field(self, bounding_box=None):
+        """Create a scalar field for the feature
+
+        Parameters
+        ----------
+        bounding_box : Optional[BoundingBox], optional
+            bounding box to evaluate the scalar field in, by default None
+
+        Returns
+        -------
+        np.ndarray
+            scalar field
+        """
+        if bounding_box is None:
+            if self.model is None:
+                raise ValueError("Must specify bounding box")
+            bounding_box = self.model.bounding_box
+        grid = bounding_box.vtk
+        points = grid.points
+        value = self.evaluate_value(points)
+        grid[self.name] = value
+        return grid