fix: 🐛 adding experimental flag for loopstructural

a boolean flag that is used to define whether experimental modules are used

LoopStructural/__init__.py

@@ -22,3 +22,4 @@
 
 logger = getLogger(__name__)
 logger.info("Imported LoopStructural")
+experimental = False

LoopStructural/analysis/__init__.py

@@ -5,10 +5,10 @@
 Various tools for analysing loopstructural models, including calculating fault intersections and fault toplogies
 """
 from LoopStructural.utils import getLogger
-
+import LoopStructural
 logger = getLogger(__name__)
-
-logger.warning("LoopStructural.analysis is experimental and may not perform as expected")
-from ._fault_displacement import displacement_missfit
-from ._fault_intersection import calculate_fault_intersections
-from ._topology import calculate_fault_topology_matrix
+if LoopStructural.experimental:
+    logger.warning("LoopStructural.analysis is experimental and may not perform as expected")
+    from ._fault_displacement import displacement_missfit
+    from ._fault_intersection import calculate_fault_intersections
+    from ._topology import calculate_fault_topology_matrix

LoopStructural/interpolators/__init__.py

@@ -4,7 +4,7 @@
 """
 from enum import IntEnum
 
-
+import LoopStructural
 class InterpolatorType(IntEnum):
     """
     Enum for the different interpolator types
@@ -37,5 +37,7 @@ class InterpolatorType(IntEnum):
 from LoopStructural.interpolators._discrete_fold_interpolator import (
     DiscreteFoldInterpolator,
 )
-from ._p1interpolator import P1Interpolator
-from ._p2interpolator import P2Interpolator
+if LoopStructural.experimental:
+    logger.warning("Using experimental interpolators: P1Interpolator and P2Interpolator")
+    from ._p1interpolator import P1Interpolator
+    from ._p2interpolator import P2Interpolator

LoopStructural/interpolators/_p2interpolator.py

@@ -86,6 +86,11 @@ def _setup_interpolator(self, **kwargs):
         self.add_value_constraints(self.interpolation_weights["cpw"])
         self.add_tangent_constraints(self.interpolation_weights["tpw"])
         # self.add_interface_constraints(self.interpolation_weights["ipw"])
+
+
+    def copy(self):
+        return P2Interpolator(self.support)
+
     def add_gradient_constraints(self, w=1.0):
         points = self.get_gradient_constraints()
         if points.shape[0] > 0:
@@ -100,7 +105,32 @@ def add_gradient_constraints(self, w=1.0):
             self.add_constraints_to_least_squares(
                 A*wt[:,None], B, elements, name="gradient")
 
-
+    def add_gradient_orthogonal_constraints(self, points, vector, w=1.0, B=0):
+        """
+        constraints scalar field to be orthogonal to a given vector
+
+        Parameters
+        ----------
+        position
+        normals
+        w
+        B
+
+        Returns
+        -------
+
+        """
+        if points.shape[0] > 0:
+            grad, elements = self.support.evaluate_shape_derivatives(points[:, :3])
+            inside = elements > -1
+            area = self.support.element_size[elements[inside]]
+            wt = np.ones(area.shape[0])
+            wt *= w * area
+            A = np.einsum("ijk,ij->ik", grad[inside, :], vector[inside, :])
+            B = np.zeros(A.shape[0])
+            elements = self.support.elements[elements[inside]]
+            self.add_constraints_to_least_squares(
+                A*wt[:,None], B, elements, name="gradient orthogonal")
 
     def add_norm_constraints(self, w=1.0):
         points = self.get_norm_constraints()