Fix/constant norm

LoopStructural/datatypes/_point.py

@@ -157,7 +157,6 @@ def vtk(
             try:
                 locations = bb.project(locations)
                 _projected = True
-                scale = bb.scale_by_projection_factor(scale)
             except Exception as e:
                 logger.error(f'Failed to project points to bounding box: {e}')
                 logger.error('Using unprojected points, this may cause issues with the glyphing')

LoopStructural/interpolators/__init__.py

@@ -3,6 +3,7 @@
 
 """
 
+
 __all__ = [
     "InterpolatorType",
     "GeologicalInterpolator",
@@ -21,39 +22,12 @@
     "StructuredGrid2D",
     "P2UnstructuredTetMesh",
 ]
-from enum import IntEnum
+from ._interpolatortype import InterpolatorType
 
 from ..utils import getLogger
 
 logger = getLogger(__name__)
 
-
-class InterpolatorType(IntEnum):
-    """
-    Enum for the different interpolator types
-
-    1-9 should cover interpolators with supports
-    9+ are data supported
-    """
-
-    BASE = 0
-    BASE_DISCRETE = 1
-    FINITE_DIFFERENCE = 2
-    DISCRETE_FOLD = 3
-    PIECEWISE_LINEAR = 4
-    PIECEWISE_QUADRATIC = 5
-    BASE_DATA_SUPPORTED = 10
-    SURFE = 11
-
-
-interpolator_string_map = {
-    "FDI": InterpolatorType.FINITE_DIFFERENCE,
-    "PLI": InterpolatorType.PIECEWISE_LINEAR,
-    "P2": InterpolatorType.PIECEWISE_QUADRATIC,
-    "P1": InterpolatorType.PIECEWISE_LINEAR,
-    "DFI": InterpolatorType.DISCRETE_FOLD,
-    'surfe': InterpolatorType.SURFE,
-}
 from ..interpolators._geological_interpolator import GeologicalInterpolator
 from ..interpolators._discrete_interpolator import DiscreteInterpolator
 from ..interpolators.supports import (
@@ -79,7 +53,7 @@ class InterpolatorType(IntEnum):
 )
 from ..interpolators._p2interpolator import P2Interpolator
 from ..interpolators._p1interpolator import P1Interpolator
-
+from ..interpolators._constant_norm import ConstantNormP1Interpolator, ConstantNormFDIInterpolator
 try:
     from ..interpolators._surfe_wrapper import SurfeRBFInterpolator
 except ImportError:
@@ -93,6 +67,24 @@ def __init__(self, *args, **kwargs):
             raise ImportError(
                 "Surfe cannot be imported. Please install Surfe. pip install surfe/ conda install -c loop3d surfe"
             ) 
+
+# Ensure compatibility between the fallback and imported class
+SurfeRBFInterpolator = SurfeRBFInterpolator
+
+
+interpolator_string_map = {
+    "FDI": InterpolatorType.FINITE_DIFFERENCE,
+    "PLI": InterpolatorType.PIECEWISE_LINEAR,
+    "P2": InterpolatorType.PIECEWISE_QUADRATIC,
+    "P1": InterpolatorType.PIECEWISE_LINEAR,
+    "DFI": InterpolatorType.DISCRETE_FOLD,
+    'surfe': InterpolatorType.SURFE,
+    "FDI_CN": InterpolatorType.FINITE_DIFFERENCE_CONSTANT_NORM,
+    "P1_CN": InterpolatorType.PIECEWISE_LINEAR_CONSTANT_NORM,
+
+}
+
+# Define the mapping after all imports
 interpolator_map = {
     InterpolatorType.BASE: GeologicalInterpolator,
     InterpolatorType.BASE_DISCRETE: DiscreteInterpolator,
@@ -102,6 +94,8 @@ def __init__(self, *args, **kwargs):
     InterpolatorType.PIECEWISE_QUADRATIC: P2Interpolator,
     InterpolatorType.BASE_DATA_SUPPORTED: GeologicalInterpolator,
     InterpolatorType.SURFE: SurfeRBFInterpolator,
+    InterpolatorType.PIECEWISE_LINEAR_CONSTANT_NORM: ConstantNormP1Interpolator,
+    InterpolatorType.FINITE_DIFFERENCE_CONSTANT_NORM: ConstantNormFDIInterpolator,
 }
 
 support_interpolator_map = {
@@ -119,9 +113,18 @@ def __init__(self, *args, **kwargs):
         3: SupportType.DataSupported,
         2: SupportType.DataSupported,
     },
+    InterpolatorType.PIECEWISE_LINEAR_CONSTANT_NORM:{
+        3: SupportType.TetMesh,
+        2: SupportType.P1Unstructured2d,
+    },
+    InterpolatorType.FINITE_DIFFERENCE_CONSTANT_NORM: {
+        3: SupportType.StructuredGrid,
+        2: SupportType.StructuredGrid2D,
+    }
 }
 
 from ._interpolator_factory import InterpolatorFactory
 from ._interpolator_builder import InterpolatorBuilder
 
-# from ._api import LoopInterpolator
+
+

LoopStructural/interpolators/_constant_norm.py

@@ -0,0 +1,195 @@
+import numpy as np
+
+from LoopStructural.interpolators._discrete_interpolator import DiscreteInterpolator
+from LoopStructural.interpolators._finite_difference_interpolator import FiniteDifferenceInterpolator
+from ._p1interpolator import P1Interpolator
+from typing import Optional, Union, Callable
+from scipy import sparse
+from LoopStructural.utils import rng
+
+class ConstantNormInterpolator:
+    """Adds a non linear constraint to an interpolator to constrain
+    the norm of the gradient to be a set value.
+
+    Returns
+    -------
+    _type_
+        _description_
+    """
+    def __init__(self, interpolator: DiscreteInterpolator,basetype):
+        """Initialise the constant norm inteprolator
+        with a discrete interpolator.
+
+        Parameters
+        ----------
+        interpolator : DiscreteInterpolator
+            The discrete interpolator to add constant norm to.
+        """
+        self.basetype = basetype
+        self.interpolator = interpolator
+        self.support = interpolator.support
+        self.random_subset = False
+        self.norm_length = 1.0
+        self.n_iterations = 20
+    def add_constant_norm(self, w:float):
+        """Add a constraint to the interpolator to constrain the norm of the gradient
+        to be a set value
+
+        Parameters
+        ----------
+        w : float
+            weighting of the constraint
+        """
+        if "constant norm" in self.interpolator.constraints:
+            _ = self.interpolator.constraints.pop("constant norm")
+
+        element_indices = np.arange(self.support.elements.shape[0])
+        if self.random_subset:
+            rng.shuffle(element_indices)
+            element_indices = element_indices[: int(0.1 * self.support.elements.shape[0])]
+        vertices, gradient, elements, inside = self.support.get_element_gradient_for_location(
+            self.support.barycentre[element_indices]
+        )
+
+        t_g = gradient[:, :, :]
+        # t_n = gradient[self.support.shared_element_relationships[:, 1], :, :]
+        v_t = np.einsum(
+            "ijk,ik->ij",
+            t_g,
+            self.interpolator.c[self.support.elements[elements]],
+        )
+
+        v_t = v_t / np.linalg.norm(v_t, axis=1)[:, np.newaxis]
+        A1 = np.einsum("ij,ijk->ik", v_t, t_g)
+
+        b = np.zeros(A1.shape[0]) + self.norm_length
+        idc = np.hstack(
+            [
+                self.support.elements[elements],
+            ]
+        )
+        self.interpolator.add_constraints_to_least_squares(A1, b, idc, w=w, name="constant norm")
+
+    def solve_system(
+        self,
+        solver: Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]] = None,
+        tol: Optional[float] = None,
+        solver_kwargs: dict = {},
+    ) -> bool:
+        """Solve the system of equations iteratively for the constant norm interpolator.
+
+        Parameters
+        ----------
+        solver : Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]], optional
+            Solver function or name, by default None
+        tol : Optional[float], optional
+            Tolerance for the solver, by default None
+        solver_kwargs : dict, optional
+            Additional arguments for the solver, by default {}
+
+        Returns
+        -------
+        bool
+            Success status of the solver
+        """
+        success = True
+        for i in range(self.n_iterations):
+            if i > 0:
+                self.add_constant_norm(w=(0.1 * i) ** 2 + 0.01)
+            # Ensure the interpolator is cast to P1Interpolator before calling solve_system
+            if isinstance(self.interpolator, self.basetype):
+                success = self.basetype.solve_system(self.interpolator, solver=solver, tol=tol, solver_kwargs=solver_kwargs)
+            else:
+                raise TypeError("self.interpolator is not an instance of P1Interpolator")
+            if not success:
+                break
+        return success
+
+class ConstantNormP1Interpolator(P1Interpolator, ConstantNormInterpolator):
+    """Constant norm interpolator using P1 base interpolator
+
+    Parameters
+    ----------
+    P1Interpolator : class
+        The P1Interpolator class.
+    ConstantNormInterpolator : class
+        The ConstantNormInterpolator class.
+    """
+    def __init__(self, support):
+        """Initialise the constant norm P1 interpolator.
+
+        Parameters
+        ----------
+        support : _type_
+            _description_
+        """
+        P1Interpolator.__init__(self, support)
+        ConstantNormInterpolator.__init__(self, self, P1Interpolator)
+
+    def solve_system(
+        self,
+        solver: Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]] = None,
+        tol: Optional[float] = None,
+        solver_kwargs: dict = {},
+    ) -> bool:
+        """Solve the system of equations for the constant norm P1 interpolator.
+
+        Parameters
+        ----------
+        solver : Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]], optional
+            Solver function or name, by default None
+        tol : Optional[float], optional
+            Tolerance for the solver, by default None
+        solver_kwargs : dict, optional
+            Additional arguments for the solver, by default {}
+
+        Returns
+        -------
+        bool
+            Success status of the solver
+        """
+        return ConstantNormInterpolator.solve_system(self, solver=solver, tol=tol, solver_kwargs=solver_kwargs)
+
+class ConstantNormFDIInterpolator(FiniteDifferenceInterpolator, ConstantNormInterpolator):
+    """Constant norm interpolator using finite difference base interpolator
+
+    Parameters
+    ----------
+    FiniteDifferenceInterpolator : class
+        The FiniteDifferenceInterpolator class.
+    ConstantNormInterpolator : class
+        The ConstantNormInterpolator class.
+    """
+    def __init__(self, support):
+        """Initialise the constant norm finite difference interpolator.
+
+        Parameters
+        ----------
+        support : _type_
+            _description_
+        """
+        FiniteDifferenceInterpolator.__init__(self, support)
+        ConstantNormInterpolator.__init__(self, self, FiniteDifferenceInterpolator)
+    def solve_system(
+        self,
+        solver: Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]] = None,
+        tol: Optional[float] = None,
+        solver_kwargs: dict = {},
+    ) -> bool:
+        """Solve the system of equations for the constant norm finite difference interpolator.
+
+        Parameters
+        ----------
+        solver : Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]], optional
+            Solver function or name, by default None
+        tol : Optional[float], optional
+            Tolerance for the solver, by default None
+        solver_kwargs : dict, optional
+            Additional arguments for the solver, by default {}
+
+        Returns
+        -------
+        bool
+            Success status of the solver
+        """
+        return ConstantNormInterpolator.solve_system(self, solver=solver, tol=tol, solver_kwargs=solver_kwargs)

LoopStructural/interpolators/_interpolatortype.py

@@ -0,0 +1,22 @@
+from enum import Enum
+
+class InterpolatorType(Enum):
+    """
+    Enum for the different interpolator types
+
+    Each value is a unique identifier.
+    """
+
+    BASE = "BASE"
+    BASE_DISCRETE = "BASE_DISCRETE"
+    FINITE_DIFFERENCE = "FINITE_DIFFERENCE"
+    DISCRETE_FOLD = "DISCRETE_FOLD"
+    PIECEWISE_LINEAR = "PIECEWISE_LINEAR"
+    PIECEWISE_QUADRATIC = "PIECEWISE_QUADRATIC"
+    BASE_DATA_SUPPORTED = "BASE_DATA_SUPPORTED"
+    SURFE = "SURFE"
+    PIECEWISE_LINEAR_CONSTANT_NORM = "PIECEWISE_LINEAR_CONSTANT_NORM"
+    FINITE_DIFFERENCE_CONSTANT_NORM = "FINITE_DIFFERENCE_CONSTANT_NORM"
+
+
+

LoopStructural/modelling/features/_base_geological_feature.py

@@ -347,6 +347,42 @@ def scalar_field(self, bounding_box=None):
         grid.cell_properties[self.name] = value
         return grid
 
+    def gradient_norm_scalar_field(self, bounding_box=None):
+        """Create a scalar field for the gradient norm of the feature
+
+        Parameters
+        ----------
+        bounding_box : Optional[BoundingBox], optional
+            bounding box to evaluate the scalar field in, by default None
+
+        Returns
+        -------
+        np.ndarray
+            scalar field of the gradient norm
+        """
+        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.structured_grid(name=self.name)
+        value = np.linalg.norm(
+            self.evaluate_gradient(bounding_box.regular_grid(local=False, order='F')),
+            axis=1,
+        )
+        if self.model is not None:
+            value = np.linalg.norm(
+                self.evaluate_gradient(
+                    self.model.scale(bounding_box.regular_grid(local=False, order='F'))
+                ),
+                axis=1,
+            )
+        grid.properties[self.name] = value
+
+        value = np.linalg.norm(
+            self.evaluate_gradient(bounding_box.cell_centres(order='F')), axis=1
+        )
+        grid.cell_properties[self.name] = value
+        return grid
     def vector_field(self, bounding_box=None, tolerance=0.05, scale=1.0):
         """Create a vector field for the feature