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:
diff --git a/LoopStructural/datatypes/_bounding_box.py b/LoopStructural/datatypes/_bounding_box.py
index 0d91a571..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.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,
@@ -435,3 +441,52 @@ 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..c8c3f668 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,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=1.0, scale_function=None, normalise=True, tolerance=0.05):
+    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
         vectors = np.copy(self.vectors)
         if normalise:
             norm = np.linalg.norm(vectors, axis=1)
@@ -133,15 +146,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
diff --git a/LoopStructural/interpolators/_finite_difference_interpolator.py b/LoopStructural/interpolators/_finite_difference_interpolator.py
index c50e332d..bdb687a4 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),
         )
@@ -172,16 +172,18 @@ 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(
+            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, :]
-
+            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 +199,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)
@@ -229,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]
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
diff --git a/LoopStructural/interpolators/supports/_3d_base_structured.py b/LoopStructural/interpolators/supports/_3d_base_structured.py
index 92785b34..495ad399 100644
--- a/LoopStructural/interpolators/supports/_3d_base_structured.py
+++ b/LoopStructural/interpolators/supports/_3d_base_structured.py
@@ -41,6 +41,10 @@ def __init__(
             raise LoopException("nsteps cannot be zero")
         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"
+            )
         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/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
 
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/_analytical_feature.py b/LoopStructural/modelling/features/_analytical_feature.py
index db967d02..27bc71d5 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,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):
+        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,
+        )
diff --git a/LoopStructural/modelling/features/_base_geological_feature.py b/LoopStructural/modelling/features/_base_geological_feature.py
index ee770151..03599cef 100644
--- a/LoopStructural/modelling/features/_base_geological_feature.py
+++ b/LoopStructural/modelling/features/_base_geological_feature.py
@@ -296,7 +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))
+
+            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 +380,14 @@ 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
diff --git a/LoopStructural/modelling/features/_cross_product_geological_feature.py b/LoopStructural/modelling/features/_cross_product_geological_feature.py
index a8fa5ad4..12ad3321 100644
--- a/LoopStructural/modelling/features/_cross_product_geological_feature.py
+++ b/LoopStructural/modelling/features/_cross_product_geological_feature.py
@@ -98,3 +98,10 @@ 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
+        )
diff --git a/LoopStructural/modelling/features/_geological_feature.py b/LoopStructural/modelling/features/_geological_feature.py
index a26f9e3a..dc790da1 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
diff --git a/LoopStructural/modelling/features/_projected_vector_feature.py b/LoopStructural/modelling/features/_projected_vector_feature.py
new file mode 100644
index 00000000..a8636725
--- /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
+        )
diff --git a/LoopStructural/modelling/features/_structural_frame.py b/LoopStructural/modelling/features/_structural_frame.py
index 466632ee..9a884536 100644
--- a/LoopStructural/modelling/features/_structural_frame.py
+++ b/LoopStructural/modelling/features/_structural_frame.py
@@ -176,3 +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):
+        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)
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)})
 
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_function_feature.py b/LoopStructural/modelling/features/fault/_fault_function_feature.py
index f4f8bcc6..e70dab1f 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")
diff --git a/LoopStructural/modelling/features/fault/_fault_segment.py b/LoopStructural/modelling/features/fault/_fault_segment.py
index 58757e58..9abd39ed 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
@@ -37,7 +44,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 +268,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/LoopStructural/modelling/intrusions/intrusion_feature.py b/LoopStructural/modelling/intrusions/intrusion_feature.py
index a2ebd5cd..60a2a13c 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
diff --git a/LoopStructural/utils/_transformation.py b/LoopStructural/utils/_transformation.py
index 7c3b7de3..75f88f2a 100644
--- a/LoopStructural/utils/_transformation.py
+++ b/LoopStructural/utils/_transformation.py
@@ -1,9 +1,13 @@
 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,11 +15,14 @@ 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
@@ -28,10 +35,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 +52,109 @@ 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
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
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"
     }
   }
 }