From 5f7c44268bb2a1856ef9bdde83297b6bdcd81a9e Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Mon, 6 Jan 2025 10:31:38 +1100
Subject: [PATCH 01/24] fix: fault segment structural frame was intiialised
 incorrectly

---
 LoopStructural/modelling/features/fault/_fault_segment.py | 3 ++-
 docs/Dockerfile                                           | 5 ++++-
 2 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/LoopStructural/modelling/features/fault/_fault_segment.py b/LoopStructural/modelling/features/fault/_fault_segment.py
index 58757e58..057ef73c 100644
--- a/LoopStructural/modelling/features/fault/_fault_segment.py
+++ b/LoopStructural/modelling/features/fault/_fault_segment.py
@@ -37,7 +37,7 @@ def __init__(
             how many integration steps for faults
         kwargs
         """
-        StructuralFrame.__init__(self, features, name, fold, model)
+        StructuralFrame.__init__(self, name=name, features=features, fold=fold,model= model)
         self.type = FeatureType.FAULT
         self.displacement = displacement
         self._faultfunction = BaseFault().fault_displacement
@@ -261,6 +261,7 @@ def evaluate_gradient(self, locations):
                 logger.error("nan slicing ")
         # need to scale with fault displacement
         v[mask, :] = self.__getitem__(1).evaluate_gradient(locations[mask, :])
+        v[mask,:]/=np.linalg.norm(v[mask,:],axis=1)[:,None]
         scale = self.displacementfeature.evaluate_value(locations[mask, :])
         v[mask, :] *= scale[:, None]
         return v
diff --git a/docs/Dockerfile b/docs/Dockerfile
index 1837a56e..a062bd12 100644
--- a/docs/Dockerfile
+++ b/docs/Dockerfile
@@ -15,6 +15,7 @@ RUN apt-get update -qq && \
     libgl1-mesa-glx\
     xvfb
 RUN conda install -c conda-forge\
+    -c loop3d\
     # python"<=3.8"\
     cython\
     numpy\ 
@@ -32,9 +33,11 @@ RUN conda install -c conda-forge\
     meshio\
     python=3.10\
     pydata-sphinx-theme\
+    pyvista\
+    loopstructuralvisualisation\
     -y
 RUN pip install git+https://github.com/geopandas/geopandas.git@v0.10.2
-RUN pip install loopstructuralvisualisation[all] geoh5py
+RUN pip install geoh5py
 RUN pip install sphinxcontrib-bibtex
 ENV TINI_VERSION v0.19.0
 ADD https://github.com/krallin/tini/releases/download/${TINI_VERSION}/tini /tini

From 9b2880e22e5c9a2356cc352b94aa90dd98e799ec Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Mon, 6 Jan 2025 10:32:14 +1100
Subject: [PATCH 02/24] fix: adding copy method as a required method for all
 geological features

---
 .../modelling/features/_analytical_feature.py    | 16 ++++++++++++++--
 .../features/_base_geological_feature.py         | 13 ++++++++++++-
 .../_cross_product_geological_feature.py         |  4 ++++
 .../modelling/features/_structural_frame.py      |  5 +++++
 .../features/fault/_fault_function_feature.py    |  3 +++
 5 files changed, 38 insertions(+), 3 deletions(-)

diff --git a/LoopStructural/modelling/features/_analytical_feature.py b/LoopStructural/modelling/features/_analytical_feature.py
index db967d02..55b9d98e 100644
--- a/LoopStructural/modelling/features/_analytical_feature.py
+++ b/LoopStructural/modelling/features/_analytical_feature.py
@@ -39,8 +39,16 @@ def __init__(
         builder=None,
     ):
         BaseFeature.__init__(self, name, model, faults, regions, builder)
-        self.vector = np.array(vector, dtype=float)
-        self.origin = np.array(origin, dtype=float)
+        try:
+            self.vector = np.array(vector, dtype=float).reshape(3)
+        except ValueError:
+            logger.error("AnalyticalGeologicalFeature: vector must be a 3 element array")
+            raise ValueError("vector must be a 3 element array")
+        try:
+            self.origin = np.array(origin, dtype=float).reshape(3)
+        except ValueError:
+            logger.error("AnalyticalGeologicalFeature: origin must be a 3 element array")
+            raise ValueError("origin must be a 3 element array")
         self.type = FeatureType.ANALYTICAL
 
     def to_json(self):
@@ -86,3 +94,7 @@ def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False):
 
     def get_data(self, value_map: Optional[dict] = None):
         return
+    def copy(self,name:Optional[str]=None):
+        if name is None:
+            name = self.name
+        return AnalyticalGeologicalFeature(name,self.vector.copy(),self.origin.copy(),list(self.regions),list(self.faults),self.model,self.builder)
\ No newline at end of file
diff --git a/LoopStructural/modelling/features/_base_geological_feature.py b/LoopStructural/modelling/features/_base_geological_feature.py
index ee770151..dcfe3819 100644
--- a/LoopStructural/modelling/features/_base_geological_feature.py
+++ b/LoopStructural/modelling/features/_base_geological_feature.py
@@ -296,7 +296,8 @@ def surfaces(
             self.regions = [
                 r for r in self.regions if r.name != self.name and r.parent.name != self.name
             ]
-            callable = lambda xyz: self.evaluate_value(self.model.scale(xyz))
+            
+            callable = lambda xyz: self.evaluate_value(self.model.scale(xyz)) if self.model is not None else self.evaluate_value(xyz)
             isosurfacer = LoopIsosurfacer(bounding_box, callable=callable)
             if name is None and self.name is not None:
                 name = self.name
@@ -375,3 +376,13 @@ def get_data(self, value_map: Optional[dict] = None):
             dictionary of data
         """
         raise NotImplementedError
+    @abstractmethod
+    def copy(self,name:Optional[str]=None):
+        """Copy the feature
+
+        Returns
+        -------
+        BaseFeature
+            copied feature
+        """
+        raise NotImplementedError
\ No newline at end of file
diff --git a/LoopStructural/modelling/features/_cross_product_geological_feature.py b/LoopStructural/modelling/features/_cross_product_geological_feature.py
index a8fa5ad4..c1fab54e 100644
--- a/LoopStructural/modelling/features/_cross_product_geological_feature.py
+++ b/LoopStructural/modelling/features/_cross_product_geological_feature.py
@@ -98,3 +98,7 @@ def max(self):
 
     def get_data(self, value_map: Optional[dict] = None):
         return
+    def copy(self,name:Optional[str]=None):
+        if name is None:
+            name = f'{self.name}_copy'
+        return CrossProductGeologicalFeature(name,self.geological_feature_a,self.geological_feature_b)
\ No newline at end of file
diff --git a/LoopStructural/modelling/features/_structural_frame.py b/LoopStructural/modelling/features/_structural_frame.py
index 466632ee..19d47c24 100644
--- a/LoopStructural/modelling/features/_structural_frame.py
+++ b/LoopStructural/modelling/features/_structural_frame.py
@@ -176,3 +176,8 @@ def get_data(self, value_map: Optional[dict] = None) -> List[Union[ValuePoints,
         for f in self.features:
             data.extend(f.get_data(value_map))
         return data
+    def copy(self,name:Optional[str]=None):
+        if name is None:
+            name = f'{self.name}_copy'
+        # !TODO check if this needs to be a deep copy
+        return StructuralFrame(name,self.features,self.fold,self.model)
\ No newline at end of file
diff --git a/LoopStructural/modelling/features/fault/_fault_function_feature.py b/LoopStructural/modelling/features/fault/_fault_function_feature.py
index f4f8bcc6..6be1343b 100644
--- a/LoopStructural/modelling/features/fault/_fault_function_feature.py
+++ b/LoopStructural/modelling/features/fault/_fault_function_feature.py
@@ -80,3 +80,6 @@ def evaluate_on_surface(self, location):
 
     def get_data(self, value_map: Optional[dict] = None):
         pass
+
+    def copy(self, name: Optional[str] = None):
+        raise NotImplementedError("Not implemented yet")
\ No newline at end of file

From d5bf487834084a68b81a242342c59488052d5763 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Mon, 6 Jan 2025 10:33:01 +1100
Subject: [PATCH 03/24] fix: allow feature to be initialised without a builder

this means the feature won't check if the interpolator is up to date. In order for this to occur a builder needs to have an up_to_date method
---
 LoopStructural/modelling/features/_geological_feature.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/LoopStructural/modelling/features/_geological_feature.py b/LoopStructural/modelling/features/_geological_feature.py
index a26f9e3a..1a1543b7 100644
--- a/LoopStructural/modelling/features/_geological_feature.py
+++ b/LoopStructural/modelling/features/_geological_feature.py
@@ -113,7 +113,9 @@ def evaluate_value(self, pos: np.ndarray, ignore_regions=False, fillnan=None) ->
         # if evaluation_points is not a numpy array try and convert
         # otherwise error
         evaluation_points = np.asarray(pos)
-        self.builder.up_to_date()
+        # if there is a builder lets make sure that the feature is up to date 
+        if self.builder is not None:
+            self.builder.up_to_date()
         # check if the points are within the display region
         v = np.zeros(evaluation_points.shape[0])
         v[:] = np.nan

From b554522fd16be060b074b5d8e1c147c6e85a66c9 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Mon, 6 Jan 2025 10:35:11 +1100
Subject: [PATCH 04/24] fix: bug with pyvista glyphs when working in a real
 world coordinate system. Fixed by reprojecting into local coordinates for the
 glyphing and then reprojecting out of the local coords after glyping.

requires a bounding box to be passed to vtk method for any vector calls.
Added a project/reproject method to the bounding box
---
 LoopStructural/datatypes/_bounding_box.py | 43 +++++++++++++++++++++++
 LoopStructural/datatypes/_point.py        | 30 ++++++++++++----
 2 files changed, 66 insertions(+), 7 deletions(-)

diff --git a/LoopStructural/datatypes/_bounding_box.py b/LoopStructural/datatypes/_bounding_box.py
index 0d91a571..b4ab6bdf 100644
--- a/LoopStructural/datatypes/_bounding_box.py
+++ b/LoopStructural/datatypes/_bounding_box.py
@@ -435,3 +435,46 @@ def structured_grid(
             properties=_vertex_data,
             name=name,
         )
+
+    def project(self, xyz):
+        """Project a point into the bounding box
+
+        Parameters
+        ----------
+        xyz : np.ndarray
+            point to project
+
+        Returns
+        -------
+        np.ndarray
+            projected point
+        """
+    
+        return (xyz - self.global_origin) / np.max((self.global_maximum-self.global_origin))#np.clip(xyz, self.origin, self.maximum) 
+
+    def reproject(self, xyz):
+        """Reproject a point from the bounding box to the global space
+
+        Parameters
+        ----------
+        xyz : np.ndarray
+            point to reproject
+
+        Returns
+        -------
+        np.ndarray
+            reprojected point
+        """ 
+        
+        return xyz * np.max((self.global_maximum - self.global_origin)) + self.global_origin
+
+    def __repr__(self):
+        return f"BoundingBox({self.origin}, {self.maximum}, {self.nsteps})"
+
+    def __str__(self):
+        return f"BoundingBox({self.origin}, {self.maximum}, {self.nsteps})"
+
+    def __eq__(self, other):
+        if not isinstance(other, BoundingBox):
+            return False
+        return np.allclose(self.origin, other.origin) and np.allclose(self.maximum, other.maximum) and np.allclose(self.nsteps, other.nsteps)
diff --git a/LoopStructural/datatypes/_point.py b/LoopStructural/datatypes/_point.py
index b48a7e79..5b52452a 100644
--- a/LoopStructural/datatypes/_point.py
+++ b/LoopStructural/datatypes/_point.py
@@ -25,11 +25,14 @@ def to_dict(self):
             ),
         }
 
-    def vtk(self):
+    def vtk(self, scalars=None):
         import pyvista as pv
 
         points = pv.PolyData(self.locations)
-        points["values"] = self.values
+        if scalars is not None and len(scalars) == len(self.locations):
+            points.point_data['scalars'] = scalars
+        else:
+            points["values"] = self.values
         return points
 
     def plot(self, pyvista_kwargs={}):
@@ -123,9 +126,9 @@ def to_dict(self):
     def from_dict(self, d):
         return VectorPoints(d['locations'], d['vectors'], d['name'], d.get('properties', None))
 
-    def vtk(self, geom='arrow', scale=1.0, scale_function=None, normalise=True, tolerance=0.05):
+    def vtk(self, geom='arrow', scale=.10, scale_function=None, normalise=True, tolerance=0.05, bb=None, scalars=None):
         import pyvista as pv
-
+        _projected=False
         vectors = np.copy(self.vectors)
         if normalise:
             norm = np.linalg.norm(vectors, axis=1)
@@ -133,15 +136,28 @@ def vtk(self, geom='arrow', scale=1.0, scale_function=None, normalise=True, tole
         if scale_function is not None:
             # vectors /= np.linalg.norm(vectors, axis=1)[:, None]
             vectors *= scale_function(self.locations)[:, None]
-        points = pv.PolyData(self.locations)
+        locations = self.locations
+        if bb is not None:
+            try:
+                locations = bb.project(locations)
+                _projected=True
+            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')
+        points = pv.PolyData(locations)
+        if scalars is not None and len(scalars) == len(self.locations):
+            points['scalars'] = scalars
         points.point_data.set_vectors(vectors, 'vectors')
         if geom == 'arrow':
             geom = pv.Arrow(scale=scale)
         elif geom == 'disc':
-            geom = pv.Disc(inner=0, outer=scale)
+            geom = pv.Disc(inner=0, outer=scale).rotate_y(90)
 
         # Perform the glyph
-        return points.glyph(orient="vectors", geom=geom, tolerance=tolerance)
+        glyphed = points.glyph(orient="vectors", geom=geom, tolerance=tolerance,scale=False)
+        if _projected:
+            glyphed.points = bb.reproject(glyphed.points)
+        return glyphed
 
     def plot(self, pyvista_kwargs={}):
         """Calls pyvista plot on the vtk object

From d4f3e12f3c9cf857ea7b057f190c247120759e1c Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Mon, 6 Jan 2025 10:35:39 +1100
Subject: [PATCH 05/24] fix: rescale/scale no longer default to in place and
 will copy the points given to them

---
 LoopStructural/modelling/core/geological_model.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py
index a73c9cec..fa27234d 100644
--- a/LoopStructural/modelling/core/geological_model.py
+++ b/LoopStructural/modelling/core/geological_model.py
@@ -1417,7 +1417,7 @@ def create_and_add_fault(
         return fault
 
     # TODO move rescale to bounding box/transformer
-    def rescale(self, points: np.ndarray, inplace: bool = True) -> np.ndarray:
+    def rescale(self, points: np.ndarray, inplace: bool = False) -> np.ndarray:
         """
         Convert from model scale to real world scale - in the future this
         should also do transformations?
@@ -1440,7 +1440,7 @@ def rescale(self, points: np.ndarray, inplace: bool = True) -> np.ndarray:
         return points
 
     # TODO move scale to bounding box/transformer
-    def scale(self, points: np.ndarray, inplace: bool = True) -> np.ndarray:
+    def scale(self, points: np.ndarray, inplace: bool = False) -> np.ndarray:
         """Take points in UTM coordinates and reproject
         into scaled model space
 

From bed530fa34a836773218a9eca1d94ec607a6ab0a Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Mon, 6 Jan 2025 10:36:06 +1100
Subject: [PATCH 06/24] fix: interpolator factory with data works

---
 .../interpolators/_interpolator_factory.py    | 25 ++++++-------------
 1 file changed, 7 insertions(+), 18 deletions(-)

diff --git a/LoopStructural/interpolators/_interpolator_factory.py b/LoopStructural/interpolators/_interpolator_factory.py
index 7882d141..f4007731 100644
--- a/LoopStructural/interpolators/_interpolator_factory.py
+++ b/LoopStructural/interpolators/_interpolator_factory.py
@@ -65,25 +65,14 @@ def create_interpolator_with_data(
         gradient_norm_constraints: Optional[np.ndarray] = None,
         gradient_constraints: Optional[np.ndarray] = None,
     ):
-        if interpolatortype is None:
-            raise ValueError("No interpolator type specified")
-        if boundingbox is None:
-            raise ValueError("No bounding box specified")
-        if nelements is None:
-            raise ValueError("No number of elements specified")
-        if isinstance(interpolatortype, str):
-            interpolatortype = InterpolatorType._member_map_[interpolatortype].numerator
-        if support is None:
-            raise Exception("Support must be specified")
-            # supporttype = support_interpolator_map[interpolatortype]
-            # support = SupportFactory.create_support(
-            #     supporttype, boundingbox, nelements, element_volume
-            # )
-        interpolator = interpolator_map[interpolatortype](support)
+        interpolator = InterpolatorFactory.create_interpolator(
+            interpolatortype, boundingbox, nelements, element_volume, support
+        )
         if value_constraints is not None:
-            interpolator.add_value_constraints(value_constraints)
+            interpolator.set_value_constraints(value_constraints)
         if gradient_norm_constraints is not None:
-            interpolator.add_gradient_constraints(gradient_norm_constraints)
+            interpolator.set_normal_constraints(gradient_norm_constraints)
         if gradient_constraints is not None:
-            interpolator.add_gradient_constraints(gradient_constraints)
+            interpolator.set_gradient_constraints(gradient_constraints)
+        interpolator.setup()
         return interpolator

From 4d000755bbd4ee5f276cb0dd6c46f65baee1dfa4 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Mon, 6 Jan 2025 10:36:26 +1100
Subject: [PATCH 07/24] ci: adding icon to docs

---
 .github/workflows/documentation.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
index 7a2a5253..f39f370f 100644
--- a/.github/workflows/documentation.yml
+++ b/.github/workflows/documentation.yml
@@ -1,4 +1,4 @@
-name: Build documentation and deploy
+name: "📚 Build documentation and deploy "
 
 on:
   workflow_dispatch:

From 8764278cec1a4cbd384787e000971edef528c76f Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Mon, 6 Jan 2025 10:36:41 +1100
Subject: [PATCH 08/24] fix: removing config details for release please

---
 release-please-config.json | 8 +-------
 1 file changed, 1 insertion(+), 7 deletions(-)

diff --git a/release-please-config.json b/release-please-config.json
index 4f4b4ee0..7d88c468 100644
--- a/release-please-config.json
+++ b/release-please-config.json
@@ -1,13 +1,7 @@
 {
   "packages": {
     "LoopStructural": {
-      "release-type": "python",
-      "types": ["feat", "fix"],
-      "bump-minor-pre-major": true,
-      "bump-minor-pre-major-pattern": "feat",
-      "bump-patch-for-minor-pre-major": true,
-      "bump-patch-for-minor-pre-major-pattern": "fix",
-      "include-v-in-tag": true
+      "release-type": "python"
     }
   }
 }

From b95fec0b84c480f6b148aa48a113a0db2e4b2da0 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Tue, 7 Jan 2025 08:59:38 +1100
Subject: [PATCH 09/24] fix: interface constraints used wrong reference to
 interpolation matrix

---
 .../_finite_difference_interpolator.py        | 29 +++++++++++++------
 1 file changed, 20 insertions(+), 9 deletions(-)

diff --git a/LoopStructural/interpolators/_finite_difference_interpolator.py b/LoopStructural/interpolators/_finite_difference_interpolator.py
index c50e332d..01a28f64 100644
--- a/LoopStructural/interpolators/_finite_difference_interpolator.py
+++ b/LoopStructural/interpolators/_finite_difference_interpolator.py
@@ -172,16 +172,28 @@ def add_interface_constraints(self, w=1.0):
         # get elements for points
         points = self.get_interface_constraints()
         if points.shape[0] > 1:
-            vertices, c, tetras, inside = self.support.get_element_for_location(
+            # vertices, c, tetras, inside = self.support.get_element_for_location(
+            #     points[:, : self.support.dimension]
+            # )
+            # # calculate volume of tetras
+            # # vecs = vertices[inside, 1:, :] - vertices[inside, 0, None, :]
+            # # vol = np.abs(np.linalg.det(vecs)) / 6
+            # A = c[inside, :]
+            # # A *= vol[:,None]
+            # idc = tetras[inside, :]
+            node_idx, inside = self.support.position_to_cell_corners(
                 points[:, : self.support.dimension]
             )
-            # calculate volume of tetras
-            # vecs = vertices[inside, 1:, :] - vertices[inside, 0, None, :]
-            # vol = np.abs(np.linalg.det(vecs)) / 6
-            A = c[inside, :]
-            # A *= vol[:,None]
-            idc = tetras[inside, :]
-
+            # print(points[inside,:].shape)
+            gi = np.zeros(self.support.n_nodes, dtype=int)
+            gi[:] = -1
+            gi[self.region] = np.arange(0, self.nx, dtype=int)
+            idc = np.zeros(node_idx.shape).astype(int)
+            idc[:] = -1
+            idc[inside, :] = gi[node_idx[inside, :]]
+            inside = np.logical_and(~np.any(idc == -1, axis=1), inside)
+            idc = idc[inside, :]
+            A = self.support.position_to_dof_coefs(points[inside, : self.support.dimension])
             for unique_id in np.unique(
                 points[
                     np.logical_and(~np.isnan(points[:, self.support.dimension]), inside),
@@ -197,7 +209,6 @@ def add_interface_constraints(self, w=1.0):
                 ).T.reshape(-1, 2)
                 interface_A = np.hstack([A[mask, :][ij[:, 0], :], -A[mask, :][ij[:, 1], :]])
                 interface_idc = np.hstack([idc[mask, :][ij[:, 0], :], idc[mask, :][ij[:, 1], :]])
-
                 # now map the index from global to region create array size of mesh
                 # initialise as np.nan, then map points inside region to 0->nx
                 gi = np.zeros(self.support.n_nodes).astype(int)

From d0a9c5b392183f52d83c4078f5c1f0eae94223e7 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Tue, 7 Jan 2025 09:00:05 +1100
Subject: [PATCH 10/24] fix: throw error if interpolation support doesn't have
 3 steps

---
 LoopStructural/interpolators/supports/_3d_base_structured.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/LoopStructural/interpolators/supports/_3d_base_structured.py b/LoopStructural/interpolators/supports/_3d_base_structured.py
index 92785b34..3dd06c1d 100644
--- a/LoopStructural/interpolators/supports/_3d_base_structured.py
+++ b/LoopStructural/interpolators/supports/_3d_base_structured.py
@@ -41,6 +41,8 @@ def __init__(
             raise LoopException("nsteps cannot be zero")
         if np.any(nsteps < 0):
             raise LoopException("nsteps cannot be negative")
+        if np.any(step_vector < 3):
+            raise LoopException("step vector cannot be less than 3. Try increasing the resolution of the interpolator")
         self._nsteps = np.array(nsteps, dtype=int) + 1
         self._step_vector = np.array(step_vector)
         self._origin = np.array(origin)

From 6e378cb7e08dd391b4ca13da9f7582cf26611acb Mon Sep 17 00:00:00 2001
From: lachlangrose <lachlangrose@users.noreply.github.com>
Date: Mon, 6 Jan 2025 22:01:27 +0000
Subject: [PATCH 11/24] style: style fixes by ruff and autoformatting by black

---
 LoopStructural/datatypes/_bounding_box.py      | 16 +++++++++++-----
 LoopStructural/datatypes/_point.py             | 18 ++++++++++++++----
 .../supports/_3d_base_structured.py            |  4 +++-
 .../modelling/features/_analytical_feature.py  | 13 +++++++++++--
 .../features/_base_geological_feature.py       | 13 +++++++++----
 .../_cross_product_geological_feature.py       |  7 +++++--
 .../modelling/features/_geological_feature.py  |  2 +-
 .../modelling/features/_structural_frame.py    |  5 +++--
 .../features/fault/_fault_function_feature.py  |  2 +-
 .../modelling/features/fault/_fault_segment.py |  4 ++--
 10 files changed, 60 insertions(+), 24 deletions(-)

diff --git a/LoopStructural/datatypes/_bounding_box.py b/LoopStructural/datatypes/_bounding_box.py
index b4ab6bdf..0db03d27 100644
--- a/LoopStructural/datatypes/_bounding_box.py
+++ b/LoopStructural/datatypes/_bounding_box.py
@@ -449,8 +449,10 @@ def project(self, xyz):
         np.ndarray
             projected point
         """
-    
-        return (xyz - self.global_origin) / np.max((self.global_maximum-self.global_origin))#np.clip(xyz, self.origin, self.maximum) 
+
+        return (xyz - self.global_origin) / np.max(
+            (self.global_maximum - self.global_origin)
+        )  # np.clip(xyz, self.origin, self.maximum)
 
     def reproject(self, xyz):
         """Reproject a point from the bounding box to the global space
@@ -464,8 +466,8 @@ def reproject(self, xyz):
         -------
         np.ndarray
             reprojected point
-        """ 
-        
+        """
+
         return xyz * np.max((self.global_maximum - self.global_origin)) + self.global_origin
 
     def __repr__(self):
@@ -477,4 +479,8 @@ def __str__(self):
     def __eq__(self, other):
         if not isinstance(other, BoundingBox):
             return False
-        return np.allclose(self.origin, other.origin) and np.allclose(self.maximum, other.maximum) and np.allclose(self.nsteps, other.nsteps)
+        return (
+            np.allclose(self.origin, other.origin)
+            and np.allclose(self.maximum, other.maximum)
+            and np.allclose(self.nsteps, other.nsteps)
+        )
diff --git a/LoopStructural/datatypes/_point.py b/LoopStructural/datatypes/_point.py
index 5b52452a..c8c3f668 100644
--- a/LoopStructural/datatypes/_point.py
+++ b/LoopStructural/datatypes/_point.py
@@ -126,9 +126,19 @@ def to_dict(self):
     def from_dict(self, d):
         return VectorPoints(d['locations'], d['vectors'], d['name'], d.get('properties', None))
 
-    def vtk(self, geom='arrow', scale=.10, scale_function=None, normalise=True, tolerance=0.05, bb=None, scalars=None):
+    def vtk(
+        self,
+        geom='arrow',
+        scale=0.10,
+        scale_function=None,
+        normalise=True,
+        tolerance=0.05,
+        bb=None,
+        scalars=None,
+    ):
         import pyvista as pv
-        _projected=False
+
+        _projected = False
         vectors = np.copy(self.vectors)
         if normalise:
             norm = np.linalg.norm(vectors, axis=1)
@@ -140,7 +150,7 @@ def vtk(self, geom='arrow', scale=.10, scale_function=None, normalise=True, tole
         if bb is not None:
             try:
                 locations = bb.project(locations)
-                _projected=True
+                _projected = True
             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')
@@ -154,7 +164,7 @@ def vtk(self, geom='arrow', scale=.10, scale_function=None, normalise=True, tole
             geom = pv.Disc(inner=0, outer=scale).rotate_y(90)
 
         # Perform the glyph
-        glyphed = points.glyph(orient="vectors", geom=geom, tolerance=tolerance,scale=False)
+        glyphed = points.glyph(orient="vectors", geom=geom, tolerance=tolerance, scale=False)
         if _projected:
             glyphed.points = bb.reproject(glyphed.points)
         return glyphed
diff --git a/LoopStructural/interpolators/supports/_3d_base_structured.py b/LoopStructural/interpolators/supports/_3d_base_structured.py
index 3dd06c1d..dba04bbb 100644
--- a/LoopStructural/interpolators/supports/_3d_base_structured.py
+++ b/LoopStructural/interpolators/supports/_3d_base_structured.py
@@ -42,7 +42,9 @@ def __init__(
         if np.any(nsteps < 0):
             raise LoopException("nsteps cannot be negative")
         if np.any(step_vector < 3):
-            raise LoopException("step vector cannot be less than 3. Try increasing the resolution of the interpolator")
+            raise LoopException(
+                "step vector cannot be less than 3. Try increasing the resolution of the interpolator"
+            )
         self._nsteps = np.array(nsteps, dtype=int) + 1
         self._step_vector = np.array(step_vector)
         self._origin = np.array(origin)
diff --git a/LoopStructural/modelling/features/_analytical_feature.py b/LoopStructural/modelling/features/_analytical_feature.py
index 55b9d98e..27bc71d5 100644
--- a/LoopStructural/modelling/features/_analytical_feature.py
+++ b/LoopStructural/modelling/features/_analytical_feature.py
@@ -94,7 +94,16 @@ def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False):
 
     def get_data(self, value_map: Optional[dict] = None):
         return
-    def copy(self,name:Optional[str]=None):
+
+    def copy(self, name: Optional[str] = None):
         if name is None:
             name = self.name
-        return AnalyticalGeologicalFeature(name,self.vector.copy(),self.origin.copy(),list(self.regions),list(self.faults),self.model,self.builder)
\ No newline at end of file
+        return AnalyticalGeologicalFeature(
+            name,
+            self.vector.copy(),
+            self.origin.copy(),
+            list(self.regions),
+            list(self.faults),
+            self.model,
+            self.builder,
+        )
diff --git a/LoopStructural/modelling/features/_base_geological_feature.py b/LoopStructural/modelling/features/_base_geological_feature.py
index dcfe3819..03599cef 100644
--- a/LoopStructural/modelling/features/_base_geological_feature.py
+++ b/LoopStructural/modelling/features/_base_geological_feature.py
@@ -296,8 +296,12 @@ def surfaces(
             self.regions = [
                 r for r in self.regions if r.name != self.name and r.parent.name != self.name
             ]
-            
-            callable = lambda xyz: self.evaluate_value(self.model.scale(xyz)) if self.model is not None else self.evaluate_value(xyz)
+
+            callable = lambda xyz: (
+                self.evaluate_value(self.model.scale(xyz))
+                if self.model is not None
+                else self.evaluate_value(xyz)
+            )
             isosurfacer = LoopIsosurfacer(bounding_box, callable=callable)
             if name is None and self.name is not None:
                 name = self.name
@@ -376,8 +380,9 @@ def get_data(self, value_map: Optional[dict] = None):
             dictionary of data
         """
         raise NotImplementedError
+
     @abstractmethod
-    def copy(self,name:Optional[str]=None):
+    def copy(self, name: Optional[str] = None):
         """Copy the feature
 
         Returns
@@ -385,4 +390,4 @@ def copy(self,name:Optional[str]=None):
         BaseFeature
             copied feature
         """
-        raise NotImplementedError
\ No newline at end of file
+        raise NotImplementedError
diff --git a/LoopStructural/modelling/features/_cross_product_geological_feature.py b/LoopStructural/modelling/features/_cross_product_geological_feature.py
index c1fab54e..12ad3321 100644
--- a/LoopStructural/modelling/features/_cross_product_geological_feature.py
+++ b/LoopStructural/modelling/features/_cross_product_geological_feature.py
@@ -98,7 +98,10 @@ def max(self):
 
     def get_data(self, value_map: Optional[dict] = None):
         return
-    def copy(self,name:Optional[str]=None):
+
+    def copy(self, name: Optional[str] = None):
         if name is None:
             name = f'{self.name}_copy'
-        return CrossProductGeologicalFeature(name,self.geological_feature_a,self.geological_feature_b)
\ No newline at end of file
+        return CrossProductGeologicalFeature(
+            name, self.geological_feature_a, self.geological_feature_b
+        )
diff --git a/LoopStructural/modelling/features/_geological_feature.py b/LoopStructural/modelling/features/_geological_feature.py
index 1a1543b7..dc790da1 100644
--- a/LoopStructural/modelling/features/_geological_feature.py
+++ b/LoopStructural/modelling/features/_geological_feature.py
@@ -113,7 +113,7 @@ def evaluate_value(self, pos: np.ndarray, ignore_regions=False, fillnan=None) ->
         # if evaluation_points is not a numpy array try and convert
         # otherwise error
         evaluation_points = np.asarray(pos)
-        # if there is a builder lets make sure that the feature is up to date 
+        # if there is a builder lets make sure that the feature is up to date
         if self.builder is not None:
             self.builder.up_to_date()
         # check if the points are within the display region
diff --git a/LoopStructural/modelling/features/_structural_frame.py b/LoopStructural/modelling/features/_structural_frame.py
index 19d47c24..9a884536 100644
--- a/LoopStructural/modelling/features/_structural_frame.py
+++ b/LoopStructural/modelling/features/_structural_frame.py
@@ -176,8 +176,9 @@ def get_data(self, value_map: Optional[dict] = None) -> List[Union[ValuePoints,
         for f in self.features:
             data.extend(f.get_data(value_map))
         return data
-    def copy(self,name:Optional[str]=None):
+
+    def copy(self, name: Optional[str] = None):
         if name is None:
             name = f'{self.name}_copy'
         # !TODO check if this needs to be a deep copy
-        return StructuralFrame(name,self.features,self.fold,self.model)
\ No newline at end of file
+        return StructuralFrame(name, self.features, self.fold, self.model)
diff --git a/LoopStructural/modelling/features/fault/_fault_function_feature.py b/LoopStructural/modelling/features/fault/_fault_function_feature.py
index 6be1343b..e70dab1f 100644
--- a/LoopStructural/modelling/features/fault/_fault_function_feature.py
+++ b/LoopStructural/modelling/features/fault/_fault_function_feature.py
@@ -82,4 +82,4 @@ def get_data(self, value_map: Optional[dict] = None):
         pass
 
     def copy(self, name: Optional[str] = None):
-        raise NotImplementedError("Not implemented yet")
\ No newline at end of file
+        raise NotImplementedError("Not implemented yet")
diff --git a/LoopStructural/modelling/features/fault/_fault_segment.py b/LoopStructural/modelling/features/fault/_fault_segment.py
index 057ef73c..32003e71 100644
--- a/LoopStructural/modelling/features/fault/_fault_segment.py
+++ b/LoopStructural/modelling/features/fault/_fault_segment.py
@@ -37,7 +37,7 @@ def __init__(
             how many integration steps for faults
         kwargs
         """
-        StructuralFrame.__init__(self, name=name, features=features, fold=fold,model= model)
+        StructuralFrame.__init__(self, name=name, features=features, fold=fold, model=model)
         self.type = FeatureType.FAULT
         self.displacement = displacement
         self._faultfunction = BaseFault().fault_displacement
@@ -261,7 +261,7 @@ def evaluate_gradient(self, locations):
                 logger.error("nan slicing ")
         # need to scale with fault displacement
         v[mask, :] = self.__getitem__(1).evaluate_gradient(locations[mask, :])
-        v[mask,:]/=np.linalg.norm(v[mask,:],axis=1)[:,None]
+        v[mask, :] /= np.linalg.norm(v[mask, :], axis=1)[:, None]
         scale = self.displacementfeature.evaluate_value(locations[mask, :])
         v[mask, :] *= scale[:, None]
         return v

From ee7f267d8e80bce304629e51fd282a480ff6c967 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Tue, 7 Jan 2025 09:29:20 +1100
Subject: [PATCH 12/24] typo, steps_vector instead of nsteps

---
 LoopStructural/interpolators/supports/_3d_base_structured.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/LoopStructural/interpolators/supports/_3d_base_structured.py b/LoopStructural/interpolators/supports/_3d_base_structured.py
index 3dd06c1d..1d4fdcec 100644
--- a/LoopStructural/interpolators/supports/_3d_base_structured.py
+++ b/LoopStructural/interpolators/supports/_3d_base_structured.py
@@ -41,7 +41,7 @@ def __init__(
             raise LoopException("nsteps cannot be zero")
         if np.any(nsteps < 0):
             raise LoopException("nsteps cannot be negative")
-        if np.any(step_vector < 3):
+        if np.any(nsteps < 3):
             raise LoopException("step vector cannot be less than 3. Try increasing the resolution of the interpolator")
         self._nsteps = np.array(nsteps, dtype=int) + 1
         self._step_vector = np.array(step_vector)

From f427ad005c5f5c30905e7026b102efbe3b8bb2ab Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Tue, 7 Jan 2025 09:29:51 +1100
Subject: [PATCH 13/24] fix: order of arguments incorrect for fault segment

---
 .../features/builders/_structural_frame_builder.py       | 4 ++--
 .../modelling/features/fault/_fault_segment.py           | 9 ++++++++-
 2 files changed, 10 insertions(+), 3 deletions(-)

diff --git a/LoopStructural/modelling/features/builders/_structural_frame_builder.py b/LoopStructural/modelling/features/builders/_structural_frame_builder.py
index 27372c17..f6601759 100644
--- a/LoopStructural/modelling/features/builders/_structural_frame_builder.py
+++ b/LoopStructural/modelling/features/builders/_structural_frame_builder.py
@@ -110,8 +110,8 @@ def __init__(
         )  # ,region=self.region))
 
         self._frame = frame(
-            self.name,
-            [
+            name=self.name,
+            features=[
                 self.builders[0].feature,
                 self.builders[1].feature,
                 self.builders[2].feature,
diff --git a/LoopStructural/modelling/features/fault/_fault_segment.py b/LoopStructural/modelling/features/fault/_fault_segment.py
index 057ef73c..24dc7d11 100644
--- a/LoopStructural/modelling/features/fault/_fault_segment.py
+++ b/LoopStructural/modelling/features/fault/_fault_segment.py
@@ -22,7 +22,14 @@ class FaultSegment(StructuralFrame):
     """
 
     def __init__(
-        self, features, name, faultfunction=None, steps=10, displacement=1.0, fold=None, model=None
+        self,
+        name: str,
+        features: list,
+        faultfunction=None,
+        steps=10,
+        displacement=1.0,
+        fold=None,
+        model=None,
     ):
         """
         A slip event of a fault

From 08e595782290a45efd164f545a969fdf75376ed6 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Tue, 7 Jan 2025 09:30:03 +1100
Subject: [PATCH 14/24] fix: add copy to intrusion

---
 LoopStructural/modelling/intrusions/intrusion_feature.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/LoopStructural/modelling/intrusions/intrusion_feature.py b/LoopStructural/modelling/intrusions/intrusion_feature.py
index a2ebd5cd..7046d757 100644
--- a/LoopStructural/modelling/intrusions/intrusion_feature.py
+++ b/LoopStructural/modelling/intrusions/intrusion_feature.py
@@ -408,3 +408,6 @@ def evaluate_value_test(self, points):
 
     def get_data(self, value_map: Optional[dict] = None):
         pass
+
+    def copy(self):
+        pass
\ No newline at end of file

From a1613b075611925789b7f8c3fe8b356497d27c4e Mon Sep 17 00:00:00 2001
From: lachlangrose <lachlangrose@users.noreply.github.com>
Date: Mon, 6 Jan 2025 22:30:53 +0000
Subject: [PATCH 15/24] style: style fixes by ruff and autoformatting by black

---
 LoopStructural/interpolators/supports/_3d_base_structured.py | 4 +++-
 LoopStructural/modelling/intrusions/intrusion_feature.py     | 2 +-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/LoopStructural/interpolators/supports/_3d_base_structured.py b/LoopStructural/interpolators/supports/_3d_base_structured.py
index 1d4fdcec..495ad399 100644
--- a/LoopStructural/interpolators/supports/_3d_base_structured.py
+++ b/LoopStructural/interpolators/supports/_3d_base_structured.py
@@ -42,7 +42,9 @@ def __init__(
         if np.any(nsteps < 0):
             raise LoopException("nsteps cannot be negative")
         if np.any(nsteps < 3):
-            raise LoopException("step vector cannot be less than 3. Try increasing the resolution of the interpolator")
+            raise LoopException(
+                "step vector cannot be less than 3. Try increasing the resolution of the interpolator"
+            )
         self._nsteps = np.array(nsteps, dtype=int) + 1
         self._step_vector = np.array(step_vector)
         self._origin = np.array(origin)
diff --git a/LoopStructural/modelling/intrusions/intrusion_feature.py b/LoopStructural/modelling/intrusions/intrusion_feature.py
index 7046d757..60a2a13c 100644
--- a/LoopStructural/modelling/intrusions/intrusion_feature.py
+++ b/LoopStructural/modelling/intrusions/intrusion_feature.py
@@ -410,4 +410,4 @@ def get_data(self, value_map: Optional[dict] = None):
         pass
 
     def copy(self):
-        pass
\ No newline at end of file
+        pass

From a4c98ed7dab98e20d74215f5d5b20f6de9f3b612 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Wed, 8 Jan 2025 14:54:48 +1100
Subject: [PATCH 16/24] fix: bb plotting in correct place when not using global
 origin

---
 LoopStructural/datatypes/_bounding_box.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/LoopStructural/datatypes/_bounding_box.py b/LoopStructural/datatypes/_bounding_box.py
index 0db03d27..e36059fd 100644
--- a/LoopStructural/datatypes/_bounding_box.py
+++ b/LoopStructural/datatypes/_bounding_box.py
@@ -410,9 +410,9 @@ def vtk(self):
             import pyvista as pv
         except ImportError:
             raise ImportError("pyvista is required for vtk support")
-        x = np.linspace(self.global_origin[0], self.global_maximum[0], self.nsteps[0])
-        y = np.linspace(self.global_origin[1], self.global_maximum[1], self.nsteps[1])
-        z = np.linspace(self.global_origin[2], self.global_maximum[2], self.nsteps[2])
+        x = np.linspace(self.global_origin[0]+self.origin[0], self.global_maximum[0], self.nsteps[0])
+        y = np.linspace(self.global_origin[1]+self.origin[1], self.global_maximum[1], self.nsteps[1])
+        z = np.linspace(self.global_origin[2]+self.origin[2], self.global_maximum[2], self.nsteps[2])
         return pv.RectilinearGrid(
             x,
             y,

From 8fed1f1c060a33fabd1aa4b0d44350f769539985 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Wed, 8 Jan 2025 14:55:04 +1100
Subject: [PATCH 17/24] fix: inequality pairs being used by FD interpolator

---
 LoopStructural/interpolators/_finite_difference_interpolator.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/LoopStructural/interpolators/_finite_difference_interpolator.py b/LoopStructural/interpolators/_finite_difference_interpolator.py
index 01a28f64..0661dbc2 100644
--- a/LoopStructural/interpolators/_finite_difference_interpolator.py
+++ b/LoopStructural/interpolators/_finite_difference_interpolator.py
@@ -97,7 +97,7 @@ def setup_interpolator(self, **kwargs):
         self.add_interface_constraints(self.interpolation_weights["ipw"])
         self.add_value_inequality_constraints()
         self.add_inequality_pairs_constraints(
-            kwargs.get('inequality_pairs', None),
+            pairs=kwargs.get('inequality_pairs', None),
             upper_bound=kwargs.get('inequality_pair_upper_bound', np.finfo(float).eps),
             lower_bound=kwargs.get('inequality_pair_lower_bound', -np.inf),
         )

From ded85d6488e407eab64d59ab205c6ca89e073ac2 Mon Sep 17 00:00:00 2001
From: lachlangrose <lachlangrose@users.noreply.github.com>
Date: Wed, 8 Jan 2025 03:55:40 +0000
Subject: [PATCH 18/24] style: style fixes by ruff and autoformatting by black

---
 LoopStructural/datatypes/_bounding_box.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/LoopStructural/datatypes/_bounding_box.py b/LoopStructural/datatypes/_bounding_box.py
index e36059fd..8518e7c8 100644
--- a/LoopStructural/datatypes/_bounding_box.py
+++ b/LoopStructural/datatypes/_bounding_box.py
@@ -410,9 +410,15 @@ def vtk(self):
             import pyvista as pv
         except ImportError:
             raise ImportError("pyvista is required for vtk support")
-        x = np.linspace(self.global_origin[0]+self.origin[0], self.global_maximum[0], self.nsteps[0])
-        y = np.linspace(self.global_origin[1]+self.origin[1], self.global_maximum[1], self.nsteps[1])
-        z = np.linspace(self.global_origin[2]+self.origin[2], self.global_maximum[2], self.nsteps[2])
+        x = np.linspace(
+            self.global_origin[0] + self.origin[0], self.global_maximum[0], self.nsteps[0]
+        )
+        y = np.linspace(
+            self.global_origin[1] + self.origin[1], self.global_maximum[1], self.nsteps[1]
+        )
+        z = np.linspace(
+            self.global_origin[2] + self.origin[2], self.global_maximum[2], self.nsteps[2]
+        )
         return pv.RectilinearGrid(
             x,
             y,

From 04d6f4a3bec62daa098def3b353b7ef5bc485ce3 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Wed, 8 Jan 2025 15:30:09 +1100
Subject: [PATCH 19/24] style: remove comments/old code

---
 .../interpolators/_finite_difference_interpolator.py  | 11 -----------
 1 file changed, 11 deletions(-)

diff --git a/LoopStructural/interpolators/_finite_difference_interpolator.py b/LoopStructural/interpolators/_finite_difference_interpolator.py
index 0661dbc2..bdb687a4 100644
--- a/LoopStructural/interpolators/_finite_difference_interpolator.py
+++ b/LoopStructural/interpolators/_finite_difference_interpolator.py
@@ -172,19 +172,9 @@ def add_interface_constraints(self, w=1.0):
         # get elements for points
         points = self.get_interface_constraints()
         if points.shape[0] > 1:
-            # vertices, c, tetras, inside = self.support.get_element_for_location(
-            #     points[:, : self.support.dimension]
-            # )
-            # # calculate volume of tetras
-            # # vecs = vertices[inside, 1:, :] - vertices[inside, 0, None, :]
-            # # vol = np.abs(np.linalg.det(vecs)) / 6
-            # A = c[inside, :]
-            # # A *= vol[:,None]
-            # idc = tetras[inside, :]
             node_idx, inside = self.support.position_to_cell_corners(
                 points[:, : self.support.dimension]
             )
-            # print(points[inside,:].shape)
             gi = np.zeros(self.support.n_nodes, dtype=int)
             gi[:] = -1
             gi[self.region] = np.arange(0, self.nx, dtype=int)
@@ -240,7 +230,6 @@ def add_gradient_constraints(self, w=1.0):
         points = self.get_gradient_constraints()
         if points.shape[0] > 0:
             # calculate unit vector for orientation data
-            # points[:,3:]/=np.linalg.norm(points[:,3:],axis=1)[:,None]
 
             node_idx, inside = self.support.position_to_cell_corners(
                 points[:, : self.support.dimension]

From 3fa75a4fbf69838a0ec9eb1ec484090922124f44 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Mon, 13 Jan 2025 09:33:28 +1100
Subject: [PATCH 20/24] fix: add a feature to project a vector onto a plane

---
 LoopStructural/modelling/features/__init__.py |   1 +
 .../features/_projected_vector_feature.py     | 112 ++++++++++++++++++
 2 files changed, 113 insertions(+)
 create mode 100644 LoopStructural/modelling/features/_projected_vector_feature.py

diff --git a/LoopStructural/modelling/features/__init__.py b/LoopStructural/modelling/features/__init__.py
index cb6f2ad3..af93f1b0 100644
--- a/LoopStructural/modelling/features/__init__.py
+++ b/LoopStructural/modelling/features/__init__.py
@@ -30,3 +30,4 @@ class FeatureType(IntEnum):
 
 from ._unconformity_feature import UnconformityFeature
 from ._analytical_feature import AnalyticalGeologicalFeature
+from ._projected_vector_feature import ProjectedVectorFeature
diff --git a/LoopStructural/modelling/features/_projected_vector_feature.py b/LoopStructural/modelling/features/_projected_vector_feature.py
new file mode 100644
index 00000000..42d7ccc3
--- /dev/null
+++ b/LoopStructural/modelling/features/_projected_vector_feature.py
@@ -0,0 +1,112 @@
+"""
+"""
+
+import numpy as np
+from typing import Optional
+
+from ...modelling.features import BaseFeature
+
+from ...utils import getLogger
+
+logger = getLogger(__name__)
+
+
+class ProjectedVectorFeature(BaseFeature):
+    def __init__(
+        self,
+        name: str,
+        vector: np.ndarray,
+        plane_feature: BaseFeature,
+    ):
+        """
+
+        Create a geological feature by projecting a vector onto a feature representing a plane
+        E.g. project a thickness vector onto an axial surface
+        
+        Parameters
+        ----------
+        name: feature name
+        vector: the vector to project
+        plane_feature: the plane
+
+
+        Parameters
+        ----------
+        name : str
+            name of the feature
+        geological_feature_a : BaseFeature
+            Left hand side of cross product
+        geological_feature_b : BaseFeature
+            Right hand side of cross product
+        """
+        super().__init__(name)
+        self.plane_feature = plane_feature
+        self.vector = vector
+        
+        self.value_feature = None
+
+    def evaluate_gradient(self, locations: np.ndarray, ignore_regions=False) -> np.ndarray:
+        """
+        Calculate the gradient of the geological feature by using numpy to
+        calculate the cross
+        product between the two existing feature gradients.
+        This means both features have to be evaluated for the locations
+
+        Parameters
+        ----------
+        locations
+
+        Returns
+        -------
+
+        """
+        
+        # project s0 onto axis plane B X A X B
+        plane_normal = self.plane_feature.evaluate_gradient(locations, ignore_regions)
+        vector = np.tile(self.vector,(locations.shape[0],1))
+
+        projected_vector = np.cross(
+            plane_normal, np.cross(vector, plane_normal, axisa=1, axisb=1), axisa=1, axisb=1
+        )
+        return projected_vector
+
+    def evaluate_value(self, evaluation_points: np.ndarray, ignore_regions=False) -> np.ndarray:
+        """
+        Return 0 because there is no value for this feature
+        Parameters
+        ----------
+        evaluation_points
+
+        Returns
+        -------
+
+        """
+        values = np.zeros(evaluation_points.shape[0])
+        if self.value_feature is not None:
+            values[:] = self.value_feature.evaluate_value(evaluation_points, ignore_regions)
+        return values
+
+    def mean(self):
+        if self.value_feature:
+            return self.value_feature.mean()
+        return 0.0
+
+    def min(self):
+        if self.value_feature:
+            return self.value_feature.min()
+        return 0.0
+
+    def max(self):
+        if self.value_feature:
+            return self.value_feature.max()
+        return 0.0
+
+    def get_data(self, value_map: Optional[dict] = None):
+        return
+
+    def copy(self, name: Optional[str] = None):
+        if name is None:
+            name = f'{self.name}_copy'
+        return ProjectedVectorFeature(
+            name=name, vector = self.vector, plane_feature=self.plane_feature
+        )

From 628e6d89c3cca3dd94f8a6a01dd072e982663993 Mon Sep 17 00:00:00 2001
From: lachlangrose <7371904+lachlangrose@users.noreply.github.com>
Date: Sun, 12 Jan 2025 22:33:56 +0000
Subject: [PATCH 21/24] style: style fixes by ruff and autoformatting by black

---
 .../modelling/features/_projected_vector_feature.py    | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/LoopStructural/modelling/features/_projected_vector_feature.py b/LoopStructural/modelling/features/_projected_vector_feature.py
index 42d7ccc3..a8636725 100644
--- a/LoopStructural/modelling/features/_projected_vector_feature.py
+++ b/LoopStructural/modelling/features/_projected_vector_feature.py
@@ -22,7 +22,7 @@ def __init__(
 
         Create a geological feature by projecting a vector onto a feature representing a plane
         E.g. project a thickness vector onto an axial surface
-        
+
         Parameters
         ----------
         name: feature name
@@ -42,7 +42,7 @@ def __init__(
         super().__init__(name)
         self.plane_feature = plane_feature
         self.vector = vector
-        
+
         self.value_feature = None
 
     def evaluate_gradient(self, locations: np.ndarray, ignore_regions=False) -> np.ndarray:
@@ -60,10 +60,10 @@ def evaluate_gradient(self, locations: np.ndarray, ignore_regions=False) -> np.n
         -------
 
         """
-        
+
         # project s0 onto axis plane B X A X B
         plane_normal = self.plane_feature.evaluate_gradient(locations, ignore_regions)
-        vector = np.tile(self.vector,(locations.shape[0],1))
+        vector = np.tile(self.vector, (locations.shape[0], 1))
 
         projected_vector = np.cross(
             plane_normal, np.cross(vector, plane_normal, axisa=1, axisb=1), axisa=1, axisb=1
@@ -108,5 +108,5 @@ def copy(self, name: Optional[str] = None):
         if name is None:
             name = f'{self.name}_copy'
         return ProjectedVectorFeature(
-            name=name, vector = self.vector, plane_feature=self.plane_feature
+            name=name, vector=self.vector, plane_feature=self.plane_feature
         )

From 47cb0a5424438399c2477a9c6846d8457639e692 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Mon, 13 Jan 2025 10:39:48 +1100
Subject: [PATCH 22/24] fix: updating transformer to take angle and translation
 as input

---
 LoopStructural/utils/_transformation.py | 106 ++++++++++++++++++++----
 1 file changed, 90 insertions(+), 16 deletions(-)

diff --git a/LoopStructural/utils/_transformation.py b/LoopStructural/utils/_transformation.py
index 7c3b7de3..8a1c6ddf 100644
--- a/LoopStructural/utils/_transformation.py
+++ b/LoopStructural/utils/_transformation.py
@@ -1,9 +1,9 @@
 import numpy as np
-from sklearn import decomposition
-
+from . import getLogger
 
+logger = getLogger(__name__)
 class EuclideanTransformation:
-    def __init__(self, dimensions=2):
+    def __init__(self, dimensions: int=2, angle: float=0, translation: np.ndarray=np.zeros(3)):
         """Transforms points into a new coordinate
         system where the main eigenvector is aligned with x
 
@@ -11,12 +11,15 @@ def __init__(self, dimensions=2):
         ----------
         dimensions : int, optional
             Do transformation in map view or on 3d volume, by default 2
+        angle : float, optional
+            Angle to rotate the points by, by default 0
+        translation : np.ndarray, default zeros
+            Translation to apply to the points, by default
         """
-        self.rotation = None
-        self.translation = None
+        self.translation = translation[:dimensions]
         self.dimensions = dimensions
-        self.angle = 0
-
+        self.angle = angle
+    
     def fit(self, points: np.ndarray):
         """Fit the transformation to a point cloud
         This function will find the main eigenvector of the point cloud
@@ -28,10 +31,16 @@ def fit(self, points: np.ndarray):
         points : np.ndarray
             xyz points as as numpy array
         """
+        try:
+            from sklearn import decomposition
+        except ImportError:
+            logger.error('scikit-learn is required for this function')
+            return
         points = np.array(points)
         if points.shape[1] < self.dimensions:
             raise ValueError("Points must have at least {} dimensions".format(self.dimensions))
         # standardise the points so that centre is 0
+        # self.translation = np.zeros(3)
         self.translation = np.mean(points, axis=0)
         # find main eigenvector and and calculate the angle of this with x
         pca = decomposition.PCA(n_components=self.dimensions).fit(
@@ -39,38 +48,103 @@ def fit(self, points: np.ndarray):
         )
         coeffs = pca.components_
         self.angle = -np.arccos(np.dot(coeffs[0, :], [1, 0]))
-        self.rotation = self._rotation(self.angle)
+        return self
+
+    @property
+    def rotation(self):
+        return self._rotation(self.angle)
+    @property
+    def inverse_rotation(self):
+        return self._rotation(-self.angle)
 
     def _rotation(self, angle):
         return np.array(
             [
                 [np.cos(angle), -np.sin(angle), 0],
                 [np.sin(angle), np.cos(angle), 0],
-                [0, 0, 1],
+                [0, 0, -1],
             ]
         )
 
     def fit_transform(self, points: np.ndarray) -> np.ndarray:
+        """Fit the transformation and transform the points"""
+
         self.fit(points)
         return self.transform(points)
 
     def transform(self, points: np.ndarray) -> np.ndarray:
-        """_summary_
+        """Transform points using the transformation and rotation
 
         Parameters
         ----------
-        points : _type_
-            _description_
+        points : np.ndarray
+            xyz points as as numpy array
 
         Returns
         -------
-        _type_
-            _description_
+        np.ndarray
+            xyz points in the transformed coordinate system
         """
-        return np.dot(points - self.translation, self.rotation)
+        points = np.array(points)
+        if points.shape[1] < self.dimensions:
+            raise ValueError("Points must have at least {} dimensions".format(self.dimensions))
+        centred = points - self.translation
 
+        return np.einsum(
+            'ik,jk->ij',
+            centred,
+            self.rotation[: self.dimensions, : self.dimensions],
+        )
     def inverse_transform(self, points: np.ndarray) -> np.ndarray:
-        return np.dot(points, self._rotation(-self.angle)) + self.translation
+        """
+        Transform points back to the original coordinate system
+
+        Parameters
+        ----------
+        points : np.ndarray
+            xyz points as as numpy array
+        
+        Returns
+        -------
+        np.ndarray
+            xyz points in the original coordinate system
+        """
+
+        return np.einsum(
+            'ik,jk->ij',
+            points,
+            self.inverse_rotation[: self.dimensions, : self.dimensions],
+        )+self.translation
 
     def __call__(self, points: np.ndarray) -> np.ndarray:
+        """
+        Transform points into the transformed space
+
+        Parameters
+        ----------
+        points : np.ndarray
+            xyz points as as numpy array
+
+        Returns
+        -------
+        np.ndarray
+            xyz points in the transformed coordinate system
+        """
+
         return self.transform(points)
+
+    def _repr_html_(self):
+
+        """
+        Provides an HTML representation of the TransRotator.
+        """
+        html_str = """
+        <div class="collapsible">
+          <button class="collapsible-button">{self.__class__.__name__}</button>
+          <div class="content">
+            <p>Translation: {self.translation}</p>
+            <p>Rotation Angle: {self.angle} degrees</p>
+          </div>
+        </div>
+        """.format(self=self)
+        return html_str

From 3dd2dafc8983e0388a9716f655f751342c3cf87b Mon Sep 17 00:00:00 2001
From: lachlangrose <7371904+lachlangrose@users.noreply.github.com>
Date: Sun, 12 Jan 2025 23:40:22 +0000
Subject: [PATCH 23/24] style: style fixes by ruff and autoformatting by black

---
 LoopStructural/utils/_transformation.py | 30 ++++++++++++++++---------
 1 file changed, 20 insertions(+), 10 deletions(-)

diff --git a/LoopStructural/utils/_transformation.py b/LoopStructural/utils/_transformation.py
index 8a1c6ddf..75f88f2a 100644
--- a/LoopStructural/utils/_transformation.py
+++ b/LoopStructural/utils/_transformation.py
@@ -2,8 +2,12 @@
 from . import getLogger
 
 logger = getLogger(__name__)
+
+
 class EuclideanTransformation:
-    def __init__(self, dimensions: int=2, angle: float=0, translation: np.ndarray=np.zeros(3)):
+    def __init__(
+        self, dimensions: int = 2, angle: float = 0, translation: np.ndarray = np.zeros(3)
+    ):
         """Transforms points into a new coordinate
         system where the main eigenvector is aligned with x
 
@@ -19,7 +23,7 @@ def __init__(self, dimensions: int=2, angle: float=0, translation: np.ndarray=np
         self.translation = translation[:dimensions]
         self.dimensions = dimensions
         self.angle = angle
-    
+
     def fit(self, points: np.ndarray):
         """Fit the transformation to a point cloud
         This function will find the main eigenvector of the point cloud
@@ -53,6 +57,7 @@ def fit(self, points: np.ndarray):
     @property
     def rotation(self):
         return self._rotation(self.angle)
+
     @property
     def inverse_rotation(self):
         return self._rotation(-self.angle)
@@ -95,6 +100,7 @@ def transform(self, points: np.ndarray) -> np.ndarray:
             centred,
             self.rotation[: self.dimensions, : self.dimensions],
         )
+
     def inverse_transform(self, points: np.ndarray) -> np.ndarray:
         """
         Transform points back to the original coordinate system
@@ -103,18 +109,21 @@ def inverse_transform(self, points: np.ndarray) -> np.ndarray:
         ----------
         points : np.ndarray
             xyz points as as numpy array
-        
+
         Returns
         -------
         np.ndarray
             xyz points in the original coordinate system
         """
 
-        return np.einsum(
-            'ik,jk->ij',
-            points,
-            self.inverse_rotation[: self.dimensions, : self.dimensions],
-        )+self.translation
+        return (
+            np.einsum(
+                'ik,jk->ij',
+                points,
+                self.inverse_rotation[: self.dimensions, : self.dimensions],
+            )
+            + self.translation
+        )
 
     def __call__(self, points: np.ndarray) -> np.ndarray:
         """
@@ -134,7 +143,6 @@ def __call__(self, points: np.ndarray) -> np.ndarray:
         return self.transform(points)
 
     def _repr_html_(self):
-
         """
         Provides an HTML representation of the TransRotator.
         """
@@ -146,5 +154,7 @@ def _repr_html_(self):
             <p>Rotation Angle: {self.angle} degrees</p>
           </div>
         </div>
-        """.format(self=self)
+        """.format(
+            self=self
+        )
         return html_str

From b22c0e91d781aa61a3c1716d8a97a70b7542b9ac Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Mon, 13 Jan 2025 12:12:42 +1100
Subject: [PATCH 24/24] fix: disable axis_function and limb_function for folds

---
 .../modelling/features/builders/_folded_feature_builder.py    | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/LoopStructural/modelling/features/builders/_folded_feature_builder.py b/LoopStructural/modelling/features/builders/_folded_feature_builder.py
index d38b9525..f33367b1 100644
--- a/LoopStructural/modelling/features/builders/_folded_feature_builder.py
+++ b/LoopStructural/modelling/features/builders/_folded_feature_builder.py
@@ -89,7 +89,7 @@ def set_fold_axis(self):
             fold_axis_rotation = get_fold_rotation_profile(self.axis_profile_type, far, fad)
             if "axis_function" in kwargs:
                 # allow predefined function to be used
-                fold_axis_rotation.set_function(kwargs["axis_function"])
+                logger.error("axis_function is deprecated, use a specific fold rotation angle profile type")
             else:
                 fold_axis_rotation.fit(params={'wavelength': kwargs.get("axis_wl", None)})
             self.fold.fold_axis_rotation = fold_axis_rotation
@@ -107,7 +107,7 @@ def set_fold_limb_rotation(self):
         fold_limb_rotation = get_fold_rotation_profile(self.limb_profile_type, flr, fld)
         if "limb_function" in kwargs:
             # allow for predefined functions to be used
-            fold_limb_rotation.set_function(kwargs["limb_function"])
+            logger.error("limb_function is deprecated, use a specific fold rotation angle profile type")
         else:
             fold_limb_rotation.fit(params={'wavelength': kwargs.get("limb_wl", None)})