diff --git a/LoopStructural/__init__.py b/LoopStructural/__init__.py
index 44819da4e..147a984fb 100644
--- a/LoopStructural/__init__.py
+++ b/LoopStructural/__init__.py
@@ -20,6 +20,7 @@
 loggers = {}
 from .modelling.core.geological_model import GeologicalModel
 from .interpolators._api import LoopInterpolator
+from .interpolators import InterpolatorBuilder
 from .datatypes import BoundingBox
 from .utils import log_to_console, log_to_file, getLogger, rng, get_levels
 
diff --git a/LoopStructural/datatypes/_bounding_box.py b/LoopStructural/datatypes/_bounding_box.py
index 8518e7c88..f7e0c6922 100644
--- a/LoopStructural/datatypes/_bounding_box.py
+++ b/LoopStructural/datatypes/_bounding_box.py
@@ -43,7 +43,7 @@ def __init__(
         if maximum is None and nsteps is not None and step_vector is not None:
             maximum = origin + nsteps * step_vector
         if origin is not None and global_origin is None:
-            global_origin = origin
+            global_origin = np.zeros(3)
         self._origin = np.array(origin)
         self._maximum = np.array(maximum)
         self.dimensions = dimensions
@@ -90,7 +90,7 @@ def global_origin(self, global_origin):
 
     @property
     def global_maximum(self):
-        return self.maximum - self.origin + self._global_origin
+        return self.maximum + self.global_origin
 
     @property
     def valid(self):
@@ -242,6 +242,8 @@ def fit(self, locations: np.ndarray, local_coordinate: bool = False) -> Bounding
             )
         origin = locations.min(axis=0)
         maximum = locations.max(axis=0)
+        origin = np.array(origin)
+        maximum = np.array(maximum)
         if local_coordinate:
             self.global_origin = origin
             self.origin = np.zeros(3)
@@ -273,15 +275,50 @@ def with_buffer(self, buffer: float = 0.2) -> BoundingBox:
         if self.origin is None or self.maximum is None:
             raise LoopValueError("Cannot create bounding box with buffer, no origin or maximum")
         # local coordinates, rescale into the original bounding boxes global coordinates
-        origin = self.origin - buffer * (self.maximum - self.origin)
-        maximum = self.maximum + buffer * (self.maximum - self.origin)
+        origin = self.origin - buffer * np.max(self.maximum - self.origin)
+        maximum = self.maximum + buffer * np.max(self.maximum - self.origin)
         return BoundingBox(
             origin=origin,
             maximum=maximum,
-            global_origin=self.global_origin + origin,
+            global_origin=self.global_origin,
             dimensions=self.dimensions,
         )
 
+    # def __call__(self, xyz):
+    #     xyz = np.array(xyz)
+    #     if len(xyz.shape) == 1:
+    #         xyz = xyz.reshape((1, -1))
+
+    #     distances = np.maximum(0,
+    #                         np.maximum(self.global_origin+self.origin - xyz,
+    #                                 xyz - self.global_maximum))
+    #     distance = np.linalg.norm(distances, axis=1)
+    #     distance[self.is_inside(xyz)] = -1
+    #     return distance
+
+    def __call__(self, xyz):
+        # Calculate center and half-extents of the box
+        center = (self.maximum + self.global_origin + self.origin) / 2
+        half_extents = (self.maximum - self.global_origin + self.origin) / 2
+
+        # Calculate the distance from point to center
+        offset = np.abs(xyz - center) - half_extents
+
+        # Inside distance: negative value based on the smallest penetration
+        inside_distance = np.min(half_extents - np.abs(xyz - center), axis=1)
+
+        # Outside distance: length of the positive components of offset
+        outside_distance = np.linalg.norm(np.maximum(offset, 0))
+
+        # If any component of offset is positive, we're outside
+        # Otherwise, we're inside and return the negative penetration distance
+        distance = np.zeros(xyz.shape[0])
+        mask = np.any(offset > 0, axis=1)
+        distance[mask] = outside_distance
+        distance[~mask] = -inside_distance[~mask]
+        return distance
+        # return outside_distance if np.any(offset > 0) else -inside_distance
+
     def get_value(self, name):
         ix, iy = self.name_map.get(name, (-1, -1))
         if ix == -1 and iy == -1:
@@ -319,7 +356,7 @@ def regular_grid(
         self,
         nsteps: Optional[Union[list, np.ndarray]] = None,
         shuffle: bool = False,
-        order: str = "C",
+        order: str = "F",
         local: bool = True,
     ) -> np.ndarray:
         """Get the grid of points from the bounding box
@@ -361,8 +398,8 @@ def regular_grid(
             rng.shuffle(locs)
         return locs
 
-    def cell_centers(self, order: str = "F") -> np.ndarray:
-        """Get the cell centers of a regular grid
+    def cell_centres(self, order: str = "F") -> np.ndarray:
+        """Get the cell centres of a regular grid
 
         Parameters
         ----------
@@ -372,7 +409,7 @@ def cell_centers(self, order: str = "F") -> np.ndarray:
         Returns
         -------
         np.ndarray
-            array of cell centers
+            array of cell centres
         """
         locs = self.regular_grid(order=order, nsteps=self.nsteps - 1)
 
@@ -434,7 +471,7 @@ def structured_grid(
         _cell_data = copy.deepcopy(cell_data)
         _vertex_data = copy.deepcopy(vertex_data)
         return StructuredGrid(
-            origin=self.global_origin,
+            origin=self.global_origin + self.origin,
             step_vector=self.step_vector,
             nsteps=self.nsteps,
             cell_properties=_cell_data,
@@ -460,6 +497,9 @@ def project(self, xyz):
             (self.global_maximum - self.global_origin)
         )  # np.clip(xyz, self.origin, self.maximum)
 
+    def scale_by_projection_factor(self, value):
+        return value / np.max((self.global_maximum - self.global_origin))
+
     def reproject(self, xyz):
         """Reproject a point from the bounding box to the global space
 
diff --git a/LoopStructural/datatypes/_point.py b/LoopStructural/datatypes/_point.py
index 61645746d..c7829b2ad 100644
--- a/LoopStructural/datatypes/_point.py
+++ b/LoopStructural/datatypes/_point.py
@@ -129,9 +129,9 @@ def from_dict(self, d):
     def vtk(
         self,
         geom='arrow',
-        scale=0.10,
+        scale=1.0,
         scale_function=None,
-        normalise=True,
+        normalise=False,
         tolerance=0.05,
         bb=None,
         scalars=None,
@@ -140,9 +140,15 @@ def vtk(
 
         _projected = False
         vectors = np.copy(self.vectors)
+
         if normalise:
             norm = np.linalg.norm(vectors, axis=1)
             vectors[norm > 0, :] /= norm[norm > 0][:, None]
+        else:
+            norm = np.linalg.norm(vectors, axis=1)
+            vectors[norm > 0, :] /= norm[norm > 0][:, None]
+            norm = norm[norm > 0] / norm[norm > 0].max()
+            vectors *= norm[:, None]
         if scale_function is not None:
             # vectors /= np.linalg.norm(vectors, axis=1)[:, None]
             vectors *= scale_function(self.locations)[:, None]
@@ -151,6 +157,7 @@ def vtk(
             try:
                 locations = bb.project(locations)
                 _projected = True
+                scale = bb.scale_by_projection_factor(scale)
             except Exception as e:
                 logger.error(f'Failed to project points to bounding box: {e}')
                 logger.error('Using unprojected points, this may cause issues with the glyphing')
@@ -161,10 +168,10 @@ def vtk(
         if geom == 'arrow':
             geom = pv.Arrow(scale=scale)
         elif geom == 'disc':
-            geom = pv.Disc(inner=0, outer=scale).rotate_y(90)
+            geom = pv.Disc(inner=0, outer=scale * 0.5, c_res=50).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)
         if _projected:
             glyphed.points = bb.reproject(glyphed.points)
         return glyphed
diff --git a/LoopStructural/datatypes/_structured_grid.py b/LoopStructural/datatypes/_structured_grid.py
index ee4e34790..dd3229681 100644
--- a/LoopStructural/datatypes/_structured_grid.py
+++ b/LoopStructural/datatypes/_structured_grid.py
@@ -44,9 +44,9 @@ def vtk(self):
             z,
         )
         for name, data in self.properties.items():
-            grid[name] = data.flatten(order="F")
+            grid[name] = data.reshape((grid.n_points, -1), order="F")
         for name, data in self.cell_properties.items():
-            grid.cell_data[name] = data.flatten(order="F")
+            grid.cell_data[name] = data.reshape((grid.n_cells, -1), order="F")
         return grid
 
     def plot(self, pyvista_kwargs={}):
@@ -63,6 +63,34 @@ def plot(self, pyvista_kwargs={}):
         except ImportError:
             logger.error("pyvista is required for vtk")
 
+    @property
+    def cell_centres(self):
+        x = np.linspace(
+            self.origin[0] + self.step_vector[0] * 0.5,
+            self.maximum[0] + self.step_vector[0] * 0.5,
+            self.nsteps[0] - 1,
+        )
+        y = np.linspace(
+            self.origin[1] + self.step_vector[1] * 0.5,
+            self.maximum[1] - self.step_vector[1] * 0.5,
+            self.nsteps[1] - 1,
+        )
+        z = np.linspace(
+            self.origin[2] + self.step_vector[2] * 0.5,
+            self.maximum[2] - self.step_vector[2] * 0.5,
+            self.nsteps[2] - 1,
+        )
+        x, y, z = np.meshgrid(x, y, z, indexing="ij")
+        return np.vstack([x.flatten(order='f'), y.flatten(order='f'), z.flatten(order='f')]).T
+
+    @property
+    def nodes(self):
+        x = np.linspace(self.origin[0], self.maximum[0], self.nsteps[0])
+        y = np.linspace(self.origin[1], self.maximum[1], self.nsteps[1])
+        z = np.linspace(self.origin[2], self.maximum[2], self.nsteps[2])
+        x, y, z = np.meshgrid(x, y, z, indexing="ij")
+        return np.vstack([x.flatten(order='f'), y.flatten(order='f'), z.flatten(order='f')]).T
+
     def merge(self, other):
         if not np.all(np.isclose(self.origin, other.origin)):
             raise ValueError("Origin of grids must be the same")
diff --git a/LoopStructural/interpolators/_discrete_interpolator.py b/LoopStructural/interpolators/_discrete_interpolator.py
index 11f739546..60f19932f 100644
--- a/LoopStructural/interpolators/_discrete_interpolator.py
+++ b/LoopStructural/interpolators/_discrete_interpolator.py
@@ -63,6 +63,20 @@ def __init__(self, support, data={}, c=None, up_to_date=False):
         logger.info("Creating discrete interpolator with {} degrees of freedom".format(self.nx))
         self.type = InterpolatorType.BASE_DISCRETE
 
+    def set_nelements(self, nelements: int) -> int:
+        return self.support.set_nelements(nelements)
+
+    @property
+    def n_elements(self) -> int:
+        """Number of elements in the interpolator
+
+        Returns
+        -------
+        int
+            number of elements, positive
+        """
+        return self.support.n_elements
+
     @property
     def nx(self) -> int:
         """Number of degrees of freedom for the interpolator
@@ -161,6 +175,7 @@ def reset(self):
         """
         self.constraints = {}
         self.c_ = 0
+        self.regularisation_scale = np.ones(self.nx)
         logger.info("Resetting interpolation constraints")
 
     def add_constraints_to_least_squares(self, A, B, idc, w=1.0, name="undefined"):
@@ -737,3 +752,6 @@ def to_dict(self):
             **super().to_dict(),
             # 'region_function':self.region_function,
         }
+
+    def vtk(self):
+        return self.support.vtk({'c': self.c})
diff --git a/LoopStructural/interpolators/_finite_difference_interpolator.py b/LoopStructural/interpolators/_finite_difference_interpolator.py
index bdb687a4e..a0b03e11f 100644
--- a/LoopStructural/interpolators/_finite_difference_interpolator.py
+++ b/LoopStructural/interpolators/_finite_difference_interpolator.py
@@ -7,12 +7,37 @@
 from ..utils import get_vectors
 from ._discrete_interpolator import DiscreteInterpolator
 from ..interpolators import InterpolatorType
-
+from scipy.spatial import KDTree
 from LoopStructural.utils import getLogger
 
 logger = getLogger(__name__)
 
 
+def compute_weighting(grid_points, gradient_constraint_points, alpha=10.0, sigma=1.0):
+    """
+    Compute weights for second derivative regularization based on proximity to gradient constraints.
+
+    Parameters:
+        grid_points (ndarray): (N, 3) array of 3D coordinates for grid cells.
+        gradient_constraint_points (ndarray): (M, 3) array of 3D coordinates for gradient constraints.
+        alpha (float): Strength of weighting increase.
+        sigma (float): Decay parameter for Gaussian-like influence.
+
+    Returns:
+        weights (ndarray): (N,) array of weights for each grid point.
+    """
+    # Build a KDTree with the gradient constraint locations
+    tree = KDTree(gradient_constraint_points)
+
+    # Find the distance from each grid point to the nearest gradient constraint
+    distances, _ = tree.query(grid_points, k=1)
+
+    # Compute weighting function (higher weight for nearby points)
+    weights = 1 + alpha * np.exp(-(distances**2) / (2 * sigma**2))
+
+    return weights
+
+
 class FiniteDifferenceInterpolator(DiscreteInterpolator):
     def __init__(self, grid, data={}):
         """
@@ -44,6 +69,7 @@ def __init__(self, grid, data={}):
         )
 
         self.type = InterpolatorType.FINITE_DIFFERENCE
+        self.use_regularisation_weight_scale = False
 
     def setup_interpolator(self, **kwargs):
         """
@@ -76,20 +102,19 @@ def setup_interpolator(self, **kwargs):
         for key in kwargs:
             self.up_to_date = False
             if "regularisation" in kwargs:
-                self.interpolation_weights["dxy"] = 0.1 * kwargs["regularisation"]
-                self.interpolation_weights["dyz"] = 0.1 * kwargs["regularisation"]
-                self.interpolation_weights["dxz"] = 0.1 * kwargs["regularisation"]
-                self.interpolation_weights["dxx"] = 0.1 * kwargs["regularisation"]
-                self.interpolation_weights["dyy"] = 0.1 * kwargs["regularisation"]
-                self.interpolation_weights["dzz"] = 0.1 * kwargs["regularisation"]
+                self.interpolation_weights["dxy"] = kwargs["regularisation"]
+                self.interpolation_weights["dyz"] = kwargs["regularisation"]
+                self.interpolation_weights["dxz"] = kwargs["regularisation"]
+                self.interpolation_weights["dxx"] = kwargs["regularisation"]
+                self.interpolation_weights["dyy"] = kwargs["regularisation"]
+                self.interpolation_weights["dzz"] = kwargs["regularisation"]
             self.interpolation_weights[key] = kwargs[key]
         # either use the default operators or the ones passed to the function
         operators = kwargs.get(
             "operators", self.support.get_operators(weights=self.interpolation_weights)
         )
-        for k, o in operators.items():
-            self.assemble_inner(o[0], o[1], name=k)
 
+        self.use_regularisation_weight_scale = kwargs.get('use_regularisation_weight_scale', False)
         self.add_norm_constraints(self.interpolation_weights["npw"])
         self.add_gradient_constraints(self.interpolation_weights["gpw"])
         self.add_value_constraints(self.interpolation_weights["cpw"])
@@ -101,6 +126,8 @@ def setup_interpolator(self, **kwargs):
             upper_bound=kwargs.get('inequality_pair_upper_bound', np.finfo(float).eps),
             lower_bound=kwargs.get('inequality_pair_lower_bound', -np.inf),
         )
+        for k, o in operators.items():
+            self.assemble_inner(o[0], o[1], name=k)
 
     def copy(self):
         """
@@ -271,6 +298,11 @@ def add_gradient_constraints(self, w=1.0):
             self.add_constraints_to_least_squares(A, B, idc[inside, :], w=w, name="gradient")
             A = np.einsum("ij,ijk->ik", dip_vector.T, T)
             self.add_constraints_to_least_squares(A, B, idc[inside, :], w=w, name="gradient")
+            # self.regularisation_scale += compute_weighting(
+            #     self.support.nodes,
+            #     points[inside, : self.support.dimension],
+            #     sigma=self.support.nsteps[0] * 10,
+            # )
             if np.sum(inside) <= 0:
                 logger.warning(
                     f" {np.sum(~inside)} \
@@ -318,7 +350,24 @@ def add_norm_constraints(self, w=1.0):
             )
             # T*=np.product(self.support.step_vector)
             # T/=self.support.step_vector[0]
-
+            # indexes, inside2 = self.support.position_to_nearby_cell_indexes(
+            # points[inside, : self.support.dimension]
+            # )
+            # indexes = indexes[inside2, :]
+
+            # corners = self.support.cell_corner_indexes(indexes)
+            # node_indexes = corners.reshape(-1, 3)
+            # indexes = self.support.global_node_indices(indexes)
+            # self.regularisation_scale[indexes]  =10
+
+            self.regularisation_scale += compute_weighting(
+                self.support.nodes,
+                points[inside, : self.support.dimension],
+                sigma=self.support.nsteps[0] * 10,
+            )
+            # global_indexes = self.support.neighbour_global_indexes().T.astype(int)
+            # close_indexes =
+            # self.regularisation_scale[global_indexes[idc[inside,:].astype(int),]]=10
             w /= 3
             for d in range(self.support.dimension):
 
@@ -454,7 +503,11 @@ def assemble_inner(self, operator, w, name='regularisation'):
             a[inside, :],
             B[inside],
             idc[inside, :],
-            w=w,
+            w=(
+                self.regularisation_scale[idc[inside, 13].astype(int)] * w
+                if self.use_regularisation_weight_scale
+                else w
+            ),
             name=name,
         )
         return
diff --git a/LoopStructural/interpolators/_geological_interpolator.py b/LoopStructural/interpolators/_geological_interpolator.py
index 5014894d5..840faf6a1 100644
--- a/LoopStructural/interpolators/_geological_interpolator.py
+++ b/LoopStructural/interpolators/_geological_interpolator.py
@@ -46,6 +46,15 @@ def __init__(self, data={}, up_to_date=False):
         self.dimensions = 3  # default to 3d
         self.support = None
 
+    @abstractmethod
+    def set_nelements(self, nelements: int) -> int:
+        pass
+
+    @property
+    @abstractmethod
+    def n_elements(self) -> int:
+        pass
+
     @property
     def data(self):
         return self._data
diff --git a/LoopStructural/interpolators/_interpolator_builder.py b/LoopStructural/interpolators/_interpolator_builder.py
index 5e87ba582..b67946627 100644
--- a/LoopStructural/interpolators/_interpolator_builder.py
+++ b/LoopStructural/interpolators/_interpolator_builder.py
@@ -1,55 +1,134 @@
 from LoopStructural.interpolators import (
-    GeologicalInterpolator,
     InterpolatorFactory,
     InterpolatorType,
 )
 from LoopStructural.datatypes import BoundingBox
-from typing import Optional, Union
+from typing import  Union, Optional
 import numpy as np
 
+from LoopStructural.interpolators._geological_interpolator import GeologicalInterpolator
+
 
 class InterpolatorBuilder:
     def __init__(
         self,
         interpolatortype: Union[str, InterpolatorType],
         bounding_box: BoundingBox,
-        nelements: int = 1000,
-        buffer: float = 0.2,
+        nelements: Optional[int] = None,
+        buffer: Optional[float] = None,
         **kwargs,
     ):
+        """This class helps initialise and setup a geological interpolator.
+
+        Parameters
+        ----------
+        interpolatortype : Union[str, InterpolatorType]
+            type of interpolator
+        bounding_box : BoundingBox
+            bounding box of the area to interpolate
+        nelements : int, optional
+            degrees of freedom of the interpolator, by default 1000
+        buffer : float, optional
+            how much of a buffer around the bounding box should be used, by default 0.2
+        """
         self.interpolatortype = interpolatortype
         self.bounding_box = bounding_box
         self.nelements = nelements
         self.buffer = buffer
         self.kwargs = kwargs
-        self.interpolator : Optional[GeologicalInterpolator]= None
-
-    def create_interpolator(self) -> 'InterpolatorBuilder':
         self.interpolator = InterpolatorFactory.create_interpolator(
-            interpolatortype=self.interpolatortype,
-            boundingbox=self.bounding_box,
-            nelements=self.nelements,
-            buffer=self.buffer,
-            **self.kwargs,
-        )
-        return self
+                    interpolatortype=self.interpolatortype,
+                    boundingbox=self.bounding_box,
+                    nelements=self.nelements,
+                    buffer=self.buffer,
+                    **self.kwargs,
+                )
 
-    def set_value_constraints(self, value_constraints: np.ndarray) -> 'InterpolatorBuilder':
+    def add_value_constraints(self, value_constraints: np.ndarray) -> 'InterpolatorBuilder':
+        """Add value constraints to the interpolator
+
+        Parameters
+        ----------
+        value_constraints : np.ndarray
+            x,y,z,value of the constraints
+
+        Returns
+        -------
+        InterpolatorBuilder
+            reference to the builder 
+        """
         if self.interpolator:
             self.interpolator.set_value_constraints(value_constraints)
         return self
 
-    def set_gradient_constraints(self, gradient_constraints: np.ndarray) -> 'InterpolatorBuilder':
+    def add_gradient_constraints(self, gradient_constraints: np.ndarray) -> 'InterpolatorBuilder':
+        """Add gradient constraints to the interpolator
+        Where g1 and g2 are two vectors that are orthogonal to the gradient
+        $'(X)\cdot g1 = 0$ and $'(X)\cdot g2 = 0$
+
+        Parameters
+        ----------
+        gradient_constraints : np.ndarray
+            x,y,z,gradient_x,gradient_y,gradient_z of the constraints
+
+        Returns
+        -------
+        InterpolatorBuilder
+            reference to the builder
+        """
+
         if self.interpolator:
             self.interpolator.set_gradient_constraints(gradient_constraints)
         return self
 
-    def set_normal_constraints(self, normal_constraints: np.ndarray) -> 'InterpolatorBuilder':
+    def add_normal_constraints(self, normal_constraints: np.ndarray) -> 'InterpolatorBuilder':
+        """Add normal constraints to the interpolator
+        Where n is the normal vector to the surface
+        $f'(X).dx = nx$
+        $f'(X).dy = ny$
+        $f'(X).dz = nz$
+        Parameters
+        ----------
+        normal_constraints : np.ndarray
+            x,y,z,nx,ny,nz of the constraints
+
+        Returns
+        -------
+        InterpolatorBuilder
+            reference to the builder
+        """
         if self.interpolator:
             self.interpolator.set_normal_constraints(normal_constraints)
         return self
+    #TODO implement/check inequalities
+    # def add_inequality_constraints(self, inequality_constraints: np.ndarray) -> 'InterpolatorBuilder':
+    #     if self.interpolator:
+    #         self.interpolator.set_value_inequality_constraints(inequality_constraints)
+    #     return self
+    # def add_inequality_pair_constraints(self, inequality_pair_constraints: np.ndarray) -> 'InterpolatorBuilder':
+    #     if self.interpolator:
+    #         self.interpolator.set_inequality_pairs_constraints(inequality_pair_constraints)
+    #     return self
+
+    def setup_interpolator(self, **kwargs) -> 'InterpolatorBuilder':
+        """This adds all of the constraints to the interpolator and
+        sets the regularisation constraints
 
-    def setup_interpolator(self, **kwargs) -> Optional[GeologicalInterpolator]:
+        Returns
+        -------
+        InterpolatorBuilder
+            reference to the builder
+        """
         if self.interpolator:
             self.interpolator.setup(**kwargs)
-            return self.interpolator
+        return self
+
+    def build(self)->GeologicalInterpolator:
+        """Builds the interpolator and returns it
+
+        Returns
+        -------
+        GeologicalInterpolator
+            The interpolator fitting all of the constraints provided
+        """
+        return self.interpolator
diff --git a/LoopStructural/interpolators/_interpolator_factory.py b/LoopStructural/interpolators/_interpolator_factory.py
index f40077315..83fa9472d 100644
--- a/LoopStructural/interpolators/_interpolator_factory.py
+++ b/LoopStructural/interpolators/_interpolator_factory.py
@@ -18,14 +18,13 @@ def create_interpolator(
         nelements: Optional[int] = None,
         element_volume: Optional[float] = None,
         support=None,
-        buffer: float = 0.2,
+        buffer: Optional[float] = 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 = interpolator_string_map[interpolatortype]
         if support is None:
diff --git a/LoopStructural/interpolators/_surfe_wrapper.py b/LoopStructural/interpolators/_surfe_wrapper.py
index 437ad0f2b..6ebc1bcf7 100644
--- a/LoopStructural/interpolators/_surfe_wrapper.py
+++ b/LoopStructural/interpolators/_surfe_wrapper.py
@@ -35,6 +35,9 @@ def __init__(self, method="single_surface"):
     def set_region(self, **kwargs):
         pass
 
+    def set_nelements(self, nelements) -> int:
+        return 0
+
     def add_gradient_constraints(self, w=1):
         points = self.get_gradient_constraints()
         if points.shape[0] > 0:
diff --git a/LoopStructural/interpolators/supports/_2d_base_unstructured.py b/LoopStructural/interpolators/supports/_2d_base_unstructured.py
index aeafe49fe..82868d0b1 100644
--- a/LoopStructural/interpolators/supports/_2d_base_unstructured.py
+++ b/LoopStructural/interpolators/supports/_2d_base_unstructured.py
@@ -68,6 +68,9 @@ def aabb_table(self):
             _initialise_aabb(self)
         return self._aabb_table
 
+    def set_nelements(self, nelements) -> int:
+        raise NotImplementedError
+
     @property
     def shared_elements(self):
         if np.sum(self._shared_elements) == 0:
diff --git a/LoopStructural/interpolators/supports/_2d_structured_grid.py b/LoopStructural/interpolators/supports/_2d_structured_grid.py
index fec585801..61e8d825d 100644
--- a/LoopStructural/interpolators/supports/_2d_structured_grid.py
+++ b/LoopStructural/interpolators/supports/_2d_structured_grid.py
@@ -66,6 +66,9 @@ def nodes(self):
     def n_nodes(self):
         return self.nsteps[0] * self.nsteps[1]
 
+    def set_nelements(self, nelements) -> int:
+        raise NotImplementedError("Cannot set number of elements for 2D structured grid")
+
     @property
     def n_elements(self):
         return self.nsteps_cells[0] * self.nsteps_cells[1]
diff --git a/LoopStructural/interpolators/supports/_3d_base_structured.py b/LoopStructural/interpolators/supports/_3d_base_structured.py
index 495ad399a..4c91322ca 100644
--- a/LoopStructural/interpolators/supports/_3d_base_structured.py
+++ b/LoopStructural/interpolators/supports/_3d_base_structured.py
@@ -3,6 +3,7 @@
 import numpy as np
 from LoopStructural.utils import getLogger
 from . import SupportType
+from typing import Tuple
 
 logger = getLogger(__name__)
 
@@ -34,6 +35,11 @@ def __init__(
         # we use property decorators to update these when different parts of
         # the geometry need to change
         # inisialise the private attributes
+        # cast to numpy array, to allow list like input
+        origin = np.array(origin)
+        nsteps = np.array(nsteps)
+        step_vector = np.array(step_vector)
+
         self.type = SupportType.BaseStructured
         if np.any(step_vector == 0):
             logger.warning(f"Step vector {step_vector} has zero values")
@@ -41,10 +47,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"
-            )
+        # 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)
@@ -56,6 +62,23 @@ def __init__(
         self.rotation_xy = rotation_xy
         self.interpolator = None
 
+    @property
+    def volume(self):
+        return np.prod(self.maximum - self.origin)
+
+    def set_nelements(self, nelements) -> int:
+        box_vol = self.volume
+        ele_vol = box_vol / nelements
+        # calculate the step vector of a regular cube
+        step_vector = np.zeros(3)
+
+        step_vector[:] = ele_vol ** (1.0 / 3.0)
+
+        # number of steps is the length of the box / step vector
+        nsteps = np.ceil((self.maximum - self.origin) / step_vector).astype(int)
+        self.nsteps = nsteps
+        return self.n_elements
+
     def to_dict(self):
         return {
             "origin": self.origin,
@@ -229,7 +252,7 @@ def rotate(self, pos):
         """ """
         return np.einsum("ijk,ik->ij", self.rotation_xy[None, :, :], pos)
 
-    def position_to_cell_index(self, pos: np.ndarray) -> np.ndarray:
+    def position_to_cell_index(self, pos: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
         """Get the indexes (i,j,k) of a cell
         that a point is inside
 
diff --git a/LoopStructural/interpolators/supports/_3d_structured_grid.py b/LoopStructural/interpolators/supports/_3d_structured_grid.py
index 6f43d11b6..19cc641bc 100644
--- a/LoopStructural/interpolators/supports/_3d_structured_grid.py
+++ b/LoopStructural/interpolators/supports/_3d_structured_grid.py
@@ -40,6 +40,8 @@ def __init__(
         self.regions["everywhere"] = np.ones(self.n_nodes).astype(bool)
 
     def onGeometryChange(self):
+        if self.interpolator is not None:
+            self.interpolator.reset()
         pass
 
     @property
diff --git a/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py b/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py
index bd527ae01..f55e8deb3 100644
--- a/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py
+++ b/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py
@@ -61,19 +61,24 @@ def __init__(
         length = self.maximum - self.minimum
         self.minimum -= length * 0.1
         self.maximum += length * 0.1
-        if aabb_nsteps is None:
-            box_vol = np.prod(self.maximum - self.minimum)
-            element_volume = box_vol / (len(self.elements) / 20)
-            # calculate the step vector of a regular cube
-            step_vector = np.zeros(3)
-            step_vector[:] = element_volume ** (1.0 / 3.0)
-            # number of steps is the length of the box / step vector
-            aabb_nsteps = np.ceil((self.maximum - self.minimum) / step_vector).astype(int)
-            # make sure there is at least one cell in every dimension
-            aabb_nsteps[aabb_nsteps < 2] = 2
-        aabb_nsteps = np.array(aabb_nsteps, dtype=int)
-        step_vector = (self.maximum - self.minimum) / (aabb_nsteps - 1)
-        self.aabb_grid = StructuredGrid(self.minimum, nsteps=aabb_nsteps, step_vector=step_vector)
+        if self.elements.shape[0] < 2000:
+            self.aabb_grid = StructuredGrid(self.minimum, nsteps=[2, 2, 2], step_vector=[1, 1, 1])
+        else:
+            if aabb_nsteps is None:
+                box_vol = np.prod(self.maximum - self.minimum)
+                element_volume = box_vol / (len(self.elements) / 20)
+                # calculate the step vector of a regular cube
+                step_vector = np.zeros(3)
+                step_vector[:] = element_volume ** (1.0 / 3.0)
+                # number of steps is the length of the box / step vector
+                aabb_nsteps = np.ceil((self.maximum - self.minimum) / step_vector).astype(int)
+                # make sure there is at least one cell in every dimension
+                aabb_nsteps[aabb_nsteps < 2] = 2
+            aabb_nsteps = np.array(aabb_nsteps, dtype=int)
+            step_vector = (self.maximum - self.minimum) / (aabb_nsteps - 1)
+            self.aabb_grid = StructuredGrid(
+                self.minimum, nsteps=aabb_nsteps, step_vector=step_vector
+            )
         # make a big table to store which tetra are in which element.
         # if this takes up too much memory it could be simplified by using sparse matrices or dict but
         # at the expense of speed
@@ -87,6 +92,9 @@ def __init__(
         self._init_face_table()
         self._initialise_aabb()
 
+    def set_nelements(self, nelements):
+        raise NotImplementedError("Cannot set number of elements for unstructured mesh")
+
     @property
     def nodes(self):
         return self._nodes
diff --git a/LoopStructural/interpolators/supports/_base_support.py b/LoopStructural/interpolators/supports/_base_support.py
index b4777a7fc..1e1a1d093 100644
--- a/LoopStructural/interpolators/supports/_base_support.py
+++ b/LoopStructural/interpolators/supports/_base_support.py
@@ -119,3 +119,7 @@ def vtk(self, node_properties={}, cell_properties={}):
         Return a vtk object
         """
         pass
+
+    @abstractmethod
+    def set_nelements(self, nelements) -> int:
+        pass
diff --git a/LoopStructural/interpolators/supports/_support_factory.py b/LoopStructural/interpolators/supports/_support_factory.py
index d31190d11..1dadc2746 100644
--- a/LoopStructural/interpolators/supports/_support_factory.py
+++ b/LoopStructural/interpolators/supports/_support_factory.py
@@ -1,4 +1,6 @@
 from LoopStructural.interpolators.supports import support_map, SupportType
+import numpy as np
+from typing import Optional
 
 
 class SupportFactory:
@@ -20,13 +22,19 @@ def from_dict(d):
 
     @staticmethod
     def create_support_from_bbox(
-        support_type, bounding_box, nelements, element_volume=None, buffer: float = 0.2
+        support_type, bounding_box, nelements, element_volume=None, buffer: Optional[float] = None
     ):
         if isinstance(support_type, str):
             support_type = SupportType._member_map_[support_type].numerator
-        bbox = bounding_box.with_buffer(buffer=buffer)
-        bbox.nelements = nelements
+        if buffer is not None:
+            bounding_box = bounding_box.with_buffer(buffer=buffer)
+        if element_volume is not None:
+            nelements = int(np.prod(bounding_box.length) / element_volume)
+        if nelements is not None:
+            bounding_box.nelements = nelements
 
         return support_map[support_type](
-            origin=bbox.origin, step_vector=bbox.step_vector, nsteps=bbox.nsteps
+            origin=bounding_box.origin,
+            step_vector=bounding_box.step_vector,
+            nsteps=bounding_box.nsteps,
         )
diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py
index 2804f33b1..ab9c3925d 100644
--- a/LoopStructural/modelling/core/geological_model.py
+++ b/LoopStructural/modelling/core/geological_model.py
@@ -711,11 +711,10 @@ def create_and_add_foliation(
         # build feature
         # series_feature = series_builder.build(**kwargs)
         series_feature = series_builder.feature
-        series_builder.build_arguments = kwargs
+        series_builder.update_build_arguments(kwargs | {"domain": True, 'tol': tol})
         # this support is built for the entire model domain? Possibly would
         # could just pass a regular grid of points - mask by any above unconformities??
-        series_builder.build_arguments['domain'] = True
-        series_builder.build_arguments["tol"] = tol
+
         series_feature.type = FeatureType.INTERPOLATED
         self._add_feature(series_feature)
         return series_feature
@@ -850,7 +849,7 @@ def create_and_add_folded_foliation(
 
         # series_feature = series_builder.build(**kwargs)
         series_feature = series_builder.feature
-        series_builder.build_arguments = kwargs
+        series_builder.update_build_arguments(kwargs)
         series_feature.type = FeatureType.INTERPOLATED
         series_feature.fold = fold
 
@@ -1264,7 +1263,7 @@ def create_and_add_domain_fault(
         # build feature
         # domain_fault = domain_fault_feature_builder.build(**kwargs)
         domain_fault = domain_fault_feature_builder.feature
-        domain_fault_feature_builder.build_arguments = kwargs
+        domain_fault_feature_builder.update_build_arguments(kwargs)
         domain_fault.type = FeatureType.DOMAINFAULT
         self._add_feature(domain_fault)
         self._add_domain_fault_below(domain_fault)
@@ -1809,7 +1808,7 @@ def get_block_model(self, name='block model'):
         grid = self.bounding_box.structured_grid(name=name)
 
         grid.cell_properties['stratigraphy'] = self.evaluate_model(
-            self.rescale(self.bounding_box.cell_centers())
+            self.rescale(self.bounding_box.cell_centres())
         )
         return grid, self.stratigraphic_ids()
 
diff --git a/LoopStructural/modelling/features/_base_geological_feature.py b/LoopStructural/modelling/features/_base_geological_feature.py
index 03599cef5..d697aa98c 100644
--- a/LoopStructural/modelling/features/_base_geological_feature.py
+++ b/LoopStructural/modelling/features/_base_geological_feature.py
@@ -338,7 +338,7 @@ def scalar_field(self, bounding_box=None):
         )
         grid.properties[self.name] = value
 
-        value = self.evaluate_value(bounding_box.cell_centers(order='F'))
+        value = self.evaluate_value(bounding_box.cell_centres(order='F'))
         grid.cell_properties[self.name] = value
         return grid
 
@@ -359,7 +359,7 @@ def vector_field(self, bounding_box=None, tolerance=0.05, scale=1.0):
             if self.model is None:
                 raise ValueError("Must specify bounding box")
             bounding_box = self.model.bounding_box
-        points = bounding_box.cell_centers()
+        points = bounding_box.cell_centres()
         value = self.evaluate_gradient(points)
         if self.model is not None:
             points = self.model.rescale(points)
diff --git a/LoopStructural/modelling/features/_geological_feature.py b/LoopStructural/modelling/features/_geological_feature.py
index dc790da1b..0575764d7 100644
--- a/LoopStructural/modelling/features/_geological_feature.py
+++ b/LoopStructural/modelling/features/_geological_feature.py
@@ -6,7 +6,6 @@
 from ...modelling.features import BaseFeature
 from ...utils import getLogger
 from ...modelling.features import FeatureType
-from ...interpolators import GeologicalInterpolator
 import numpy as np
 from typing import Optional, List, Union
 from ...datatypes import ValuePoints, VectorPoints
@@ -40,8 +39,7 @@ class GeologicalFeature(BaseFeature):
     def __init__(
         self,
         name: str,
-        interpolator: GeologicalInterpolator,
-        builder=None,
+        builder,
         regions: list = [],
         faults: list = [],
         model=None,
@@ -61,8 +59,8 @@ def __init__(
         """
         BaseFeature.__init__(self, name, model, faults, regions, builder)
         self.name = name
-        self.interpolator = interpolator
         self.builder = builder
+        self.interpolator = self.builder.interpolator if self.builder is not None else None
         self.type = FeatureType.INTERPOLATED
 
     def to_json(self):
@@ -262,7 +260,6 @@ def copy(self, name=None):
             regions=[],  # feature.regions.copy(),  # don't want to share regionsbetween unconformity and # feature.regions,
             builder=self.builder,
             model=self.model,
-            interpolator=self.interpolator,
         )
         return feature
 
diff --git a/LoopStructural/modelling/features/_lambda_geological_feature.py b/LoopStructural/modelling/features/_lambda_geological_feature.py
index 78043f026..6c2f1d64c 100644
--- a/LoopStructural/modelling/features/_lambda_geological_feature.py
+++ b/LoopStructural/modelling/features/_lambda_geological_feature.py
@@ -12,7 +12,6 @@
 
 
 class LambdaGeologicalFeature(BaseFeature):
-
     def __init__(
         self,
         function: Optional[Callable[[np.ndarray], np.ndarray]] = None,
diff --git a/LoopStructural/modelling/features/_unconformity_feature.py b/LoopStructural/modelling/features/_unconformity_feature.py
index 79d1627a0..c3a8eae1e 100644
--- a/LoopStructural/modelling/features/_unconformity_feature.py
+++ b/LoopStructural/modelling/features/_unconformity_feature.py
@@ -24,7 +24,6 @@ def __init__(self, feature: GeologicalFeature, value: float, sign=True, onlap=Fa
             regions=[],  # feature.regions.copy(),  # don't want to share regionsbetween unconformity and # feature.regions,
             builder=feature.builder,
             model=feature.model,
-            interpolator=feature.interpolator,
         )
         self.value = value
         self.type = FeatureType.UNCONFORMITY if onlap is False else FeatureType.ONLAPUNCONFORMITY
diff --git a/LoopStructural/modelling/features/builders/_base_builder.py b/LoopStructural/modelling/features/builders/_base_builder.py
index 5f05b453a..9fd2ab073 100644
--- a/LoopStructural/modelling/features/builders/_base_builder.py
+++ b/LoopStructural/modelling/features/builders/_base_builder.py
@@ -47,10 +47,12 @@ def feature(self):
     def build_arguments(self):
         return self._build_arguments
 
-    @build_arguments.setter
-    def build_arguments(self, build_arguments):
+    def update_build_arguments(self, build_arguments):
         # self._build_arguments = {}
+        logger.info(f"Setting build arguments for {self.name}")
         for k, i in build_arguments.items():
+            logger.info(f"Setting {k} to {i} for {self.name}")
+            logger.info(f"{k} is currrently {self._build_arguments.get(k, None)}")
             if i != self._build_arguments.get(k, None):
                 logger.info(f"Setting {k} to {i} for {self.name}")
                 self._build_arguments[k] = i
diff --git a/LoopStructural/modelling/features/builders/_geological_feature_builder.py b/LoopStructural/modelling/features/builders/_geological_feature_builder.py
index 62d979eb9..016ffb758 100644
--- a/LoopStructural/modelling/features/builders/_geological_feature_builder.py
+++ b/LoopStructural/modelling/features/builders/_geological_feature_builder.py
@@ -25,7 +25,6 @@
 from ....utils.helper import (
     get_data_bounding_box_map as get_data_bounding_box,
 )
-from ....utils import RegionEverywhere
 from ....interpolators import DiscreteInterpolator
 from ....interpolators import InterpolatorFactory
 
@@ -39,7 +38,6 @@ def __init__(
         bounding_box,
         nelements: int = 1000,
         name="Feature",
-        interpolation_region=None,
         model=None,
         **kwargs,
     ):
@@ -61,7 +59,7 @@ def __init__(
             nelements=nelements,
             buffer=kwargs.get("buffer", 0.2),
         )
-
+        
         if not issubclass(type(interpolator), GeologicalInterpolator):
             raise TypeError(
                 "interpolator is {} and must be a GeologicalInterpolator".format(type(interpolator))
@@ -78,21 +76,17 @@ def __init__(
         )
         self.data = pd.DataFrame(columns=header)
         self.data_added = False
-        self._interpolation_region = None
-        self.interpolation_region = interpolation_region
-        if self.interpolation_region is not None:
-            self._interpolator.set_region(region=self.interpolation_region)
 
         self._feature = GeologicalFeature(
             self._name,
-            self._interpolator,
             builder=self,
             regions=[],
             faults=self.faults,
         )
         self._orthogonal_features = {}
         self._equality_constraints = {}
-
+        # add default parameters
+        self.update_build_arguments({'cpw':1.0,'npw':1.0,'regularisation':1.0,'nelements':self.interpolator.n_elements})
     def set_not_up_to_date(self, caller):
         logger.info(
             f"Setting {self.name} to not up to date from an instance of {caller.__class__.__name__}"
@@ -104,20 +98,12 @@ def set_not_up_to_date(self, caller):
     def interpolator(self):
         return self._interpolator
 
-    @property
-    def interpolation_region(self):
-        return self._interpolation_region
-
-    @interpolation_region.setter
-    def interpolation_region(self, interpolation_region):
-        if interpolation_region is not None:
-            self._interpolation_region = interpolation_region
-            self._interpolator.set_region(region=self._interpolation_region)
-        else:
-            self._interpolation_region = RegionEverywhere()
-            self._interpolator.set_region(region=self._interpolation_region)
-        logger.info(f'Setting interpolation region {self.name}')
-        self._up_to_date = False
+    @interpolator.setter
+    def interpolator(self, interpolator):
+        if not issubclass(type(interpolator), GeologicalInterpolator):
+            raise TypeError(
+                "interpolator is {} and must be a GeologicalInterpolator".format(type(interpolator))
+            )
 
     def add_data_from_data_frame(self, data_frame, overwrite=False):
         """
@@ -183,6 +169,7 @@ def add_data_to_interpolator(self, constrained=False, force_constrained=False, *
 
         """
         if self.data_added:
+            logger.info("Data already added to interpolator")
             return
         # first move the data for the fault
         logger.info(f"Adding {len(self.faults)} faults to {self.name}")
@@ -475,7 +462,7 @@ def check_interpolation_geometry(self, data, buffer=0.3):
         is big enough to capture the faulted feature.
         """
         if self.interpolator.support is not None:
-
+            logger.info(f"Checking interpolation geometry for {self.name}")
             origin = self.interpolator.support.origin
             maximum = self.interpolator.support.maximum
             pts = self.model.bounding_box.with_buffer(buffer).regular_grid(local=True)
@@ -488,7 +475,7 @@ def check_interpolation_geometry(self, data, buffer=0.3):
             ]
             self.interpolator.support.origin = origin
             self.interpolator.support.maximum = maximum
-
+            self.update_build_arguments({'nelements':self.interpolator.n_elements})
     def build(self, data_region=None, **kwargs):
         """
         Runs the interpolation and builds the geological feature
@@ -511,6 +498,15 @@ def build(self, data_region=None, **kwargs):
         for f in self.faults:
             f.builder.update()
         domain = kwargs.get("domain", None)
+        if 'nelements' in kwargs:
+            # if the number of elements has changed then update the interpolator
+            logger.info(f'Interpolator has {self.interpolator.n_elements} elements')
+            logger.info(f'Kwargs has {kwargs["nelements"]} elements')
+            if self.interpolator.n_elements != kwargs['nelements']:
+                logger.info('Setting nelements to {} for {}'.format(kwargs['nelements'], self.name))
+                self.build_arguments['nelements'] = self.interpolator.set_nelements(kwargs['nelements'])
+                logger.info(f'Interpolator nelements {self.interpolator.n_elements}')
+                self._up_to_date = False
         if domain:
             self.check_interpolation_geometry(None)
         self.add_data_to_interpolator(**kwargs)
diff --git a/LoopStructural/modelling/features/builders/_structural_frame_builder.py b/LoopStructural/modelling/features/builders/_structural_frame_builder.py
index f6601759b..b73907a65 100644
--- a/LoopStructural/modelling/features/builders/_structural_frame_builder.py
+++ b/LoopStructural/modelling/features/builders/_structural_frame_builder.py
@@ -120,7 +120,12 @@ def __init__(
             model=self.model,
         )
         self._frame.builder = self
-
+    @property
+    def build_arguments(self):
+        return self.builders[0].build_arguments
+    def update_build_arguments(self, kwargs):
+        for i in range(3):
+            self.builders[i].update_build_arguments(kwargs)
     @property
     def frame(self):
         return self._frame
@@ -197,7 +202,7 @@ def setup(self, w1=1.0, w2=1.0, w3=1.0, **kwargs):
         if len(self.builders[0].data) > 0:
             logger.info(f"Building {self.name} coordinate 0")
             kwargs["regularisation"] = regularisation[0]
-            self.builders[0].build_arguments = kwargs
+            self.builders[0].update_build_arguments(kwargs)
             kwargs.pop("fold", None)
 
         # make sure that all of the coordinates are using the same region
@@ -206,7 +211,7 @@ def setup(self, w1=1.0, w2=1.0, w3=1.0, **kwargs):
             if w2 > 0:
                 self.builders[2].add_orthogonal_feature(self.builders[0].feature, w2, step=step)
             kwargs["regularisation"] = regularisation[2]
-            self.builders[2].build_arguments = kwargs
+            self.builders[2].update_build_arguments(kwargs)
 
         if len(self.builders[1].data) > 0:
             logger.info(f"Building {self.name} coordinate 1")
@@ -215,7 +220,7 @@ def setup(self, w1=1.0, w2=1.0, w3=1.0, **kwargs):
             if w3 > 0 and len(self.builders[2].data) > 0:
                 self.builders[1].add_orthogonal_feature(self.builders[2].feature, w2, step=step)
             kwargs["regularisation"] = regularisation[1]
-            self.builders[1].build_arguments = kwargs
+            self.builders[1].update_build_arguments(kwargs)
 
         if len(self.builders[2].data) == 0:
             from LoopStructural.modelling.features import (
diff --git a/LoopStructural/modelling/features/fault/_fault_segment.py b/LoopStructural/modelling/features/fault/_fault_segment.py
index 9abd39ed6..c36e75bf0 100644
--- a/LoopStructural/modelling/features/fault/_fault_segment.py
+++ b/LoopStructural/modelling/features/fault/_fault_segment.py
@@ -270,7 +270,7 @@ def evaluate_gradient(self, locations):
         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]
+        v[mask, :] *= scale[:, None] * self.displacement
         return v
 
     def evaluate_displacement(self, points):
diff --git a/LoopStructural/modelling/features/fold/fold_function/_base_fold_rotation_angle.py b/LoopStructural/modelling/features/fold/fold_function/_base_fold_rotation_angle.py
index 486f1b236..54a537cb9 100644
--- a/LoopStructural/modelling/features/fold/fold_function/_base_fold_rotation_angle.py
+++ b/LoopStructural/modelling/features/fold/fold_function/_base_fold_rotation_angle.py
@@ -12,7 +12,6 @@
 
 
 class BaseFoldRotationAngleProfile(metaclass=ABCMeta):
-
     def __init__(
         self,
         rotation_angle: Optional[npt.NDArray[np.float64]] = None,
diff --git a/LoopStructural/modelling/features/fold/fold_function/_lambda_fold_rotation_angle.py b/LoopStructural/modelling/features/fold/fold_function/_lambda_fold_rotation_angle.py
index 2c1c5e326..79c01cba9 100644
--- a/LoopStructural/modelling/features/fold/fold_function/_lambda_fold_rotation_angle.py
+++ b/LoopStructural/modelling/features/fold/fold_function/_lambda_fold_rotation_angle.py
@@ -8,7 +8,6 @@
 
 
 class LambdaFoldRotationAngleProfile(BaseFoldRotationAngleProfile):
-
     def __init__(
         self,
         fn: Callable[[np.ndarray], np.ndarray],
diff --git a/LoopStructural/modelling/features/fold/fold_function/_trigo_fold_rotation_angle.py b/LoopStructural/modelling/features/fold/fold_function/_trigo_fold_rotation_angle.py
index 001ed5afc..1908179ec 100644
--- a/LoopStructural/modelling/features/fold/fold_function/_trigo_fold_rotation_angle.py
+++ b/LoopStructural/modelling/features/fold/fold_function/_trigo_fold_rotation_angle.py
@@ -8,7 +8,6 @@
 
 
 class TrigoFoldRotationAngleProfile(BaseFoldRotationAngleProfile):
-
     def __init__(
         self,
         rotation_angle: Optional[npt.NDArray[np.float64]] = None,
diff --git a/LoopStructural/utils/_surface.py b/LoopStructural/utils/_surface.py
index 020e91b2a..d9a94775e 100644
--- a/LoopStructural/utils/_surface.py
+++ b/LoopStructural/utils/_surface.py
@@ -89,7 +89,7 @@ def fit(
             raise ValueError("No interpolator of callable function set")
 
         surfaces = []
-        all_values = self.callable(self.bounding_box.regular_grid(local=local))
+        all_values = self.callable(self.bounding_box.regular_grid(local=local, order='C'))
         ## set value to mean value if its not specified
         if values is None:
             values = [((np.nanmax(all_values) - np.nanmin(all_values)) / 2) + np.nanmin(all_values)]
diff --git a/LoopStructural/utils/_transformation.py b/LoopStructural/utils/_transformation.py
index 75f88f2a1..af7116fec 100644
--- a/LoopStructural/utils/_transformation.py
+++ b/LoopStructural/utils/_transformation.py
@@ -6,7 +6,11 @@
 
 class EuclideanTransformation:
     def __init__(
-        self, dimensions: int = 2, angle: float = 0, translation: np.ndarray = np.zeros(3)
+        self,
+        dimensions: int = 2,
+        angle: float = 0,
+        translation: np.ndarray = np.zeros(3),
+        fit_rotation: bool = True,
     ):
         """Transforms points into a new coordinate
         system where the main eigenvector is aligned with x
@@ -23,6 +27,7 @@ def __init__(
         self.translation = translation[:dimensions]
         self.dimensions = dimensions
         self.angle = angle
+        self.fit_rotation = fit_rotation
 
     def fit(self, points: np.ndarray):
         """Fit the transformation to a point cloud
@@ -45,13 +50,16 @@ def fit(self, points: np.ndarray):
             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)
+        self.translation = np.mean(points[:, : self.dimensions], axis=0)
         # find main eigenvector and and calculate the angle of this with x
-        pca = decomposition.PCA(n_components=self.dimensions).fit(
-            points[:, : self.dimensions] - self.translation[None, : self.dimensions]
-        )
-        coeffs = pca.components_
-        self.angle = -np.arccos(np.dot(coeffs[0, :], [1, 0]))
+        if self.fit_rotation:
+            pca = decomposition.PCA(n_components=self.dimensions).fit(
+                points[:, : self.dimensions] - self.translation[None, : self.dimensions]
+            )
+            coeffs = pca.components_
+            self.angle = -np.arccos(np.dot(coeffs[0, :], [1, 0]))
+        else:
+            self.angle = 0
         return self
 
     @property
@@ -93,13 +101,15 @@ def transform(self, points: np.ndarray) -> np.ndarray:
         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(
+        centred = points[:, : self.dimensions] - self.translation[None, :]
+        rotated = np.einsum(
             'ik,jk->ij',
             centred,
             self.rotation[: self.dimensions, : self.dimensions],
         )
+        transformed_points = np.copy(points)
+        transformed_points[:, : self.dimensions] = rotated
+        return transformed_points
 
     def inverse_transform(self, points: np.ndarray) -> np.ndarray:
         """
@@ -115,15 +125,20 @@ def inverse_transform(self, points: np.ndarray) -> np.ndarray:
         np.ndarray
             xyz points in the original coordinate system
         """
-
-        return (
+        inversed = (
             np.einsum(
                 'ik,jk->ij',
-                points,
+                points[: self.dimensions],
                 self.inverse_rotation[: self.dimensions, : self.dimensions],
             )
             + self.translation
         )
+        inversed = (
+            np.vstack([inversed, points[self.dimensions :]])
+            if points.shape[1] > self.dimensions
+            else inversed
+        )
+        return inversed
 
     def __call__(self, points: np.ndarray) -> np.ndarray:
         """
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 686b817d6..33d246572 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -59,6 +59,9 @@
     # citations
     "myst_parser",
     "sphinxcontrib.bibtex",
+    "pyvista.ext.plot_directive",
+    "pyvista.ext.viewer_directive",
+    "sphinx_design",
 ]
 bibtex_bibfiles = ["docs_references.bib"]
 bibtex_default_style = "plain"
diff --git a/docs/source/getting_started/background.rst b/docs/source/getting_started/background.rst
index 12c41384c..1fef48d67 100644
--- a/docs/source/getting_started/background.rst
+++ b/docs/source/getting_started/background.rst
@@ -1,6 +1,6 @@
 Background
 ==========
-LoopStructural is an opensource Python 3.6+ library for 3D geological modelling where geological objects are represented by an implicit function. 
+LoopStructural is an opensource Python 3.9+ library for 3D geological modelling where geological objects are represented by an implicit function. 
 Where the implicit function represents the distance or pseudodistance to a reference horizion. 
 There is no analytical function that can be used to represent the geometry of geological objects, so the implicit function has to be approximated. 
 The implicit function is approximated from the observations of the surface
diff --git a/docs/source/getting_started/installation.rst b/docs/source/getting_started/installation.rst
index 70680bf8f..9f1934ddb 100644
--- a/docs/source/getting_started/installation.rst
+++ b/docs/source/getting_started/installation.rst
@@ -1,6 +1,9 @@
 Installation
 ====================
-LoopStructural is supported and tested on Python 3.6+ and can be installed on Linux, Windows and Mac. 
+
+**Last updated:** 4/02/2025
+
+LoopStructural is supported and tested on Python 3.9+ and can be installed on Linux, Windows and Mac. 
 We recommend installing LoopStructural into clean python environment. Either using anaconda or python virtual environments. 
 There are three ways of installing LoopStructural onto your system:
 
@@ -10,24 +13,19 @@ Installing from pip or conda
 .. code-block::
 
     pip install LoopStructural
+    pip install LoopStructural[all] # to include all optional dependencies
 
 
 .. code-block::
 
     conda install -c conda-forge -c loop3d loopstructural
-
+    
     
 Compiling LoopStructural from source
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 You can install the most recent version of LoopStructural by cloning it from GitHub. 
-You will need to have a C/C++ development environment for compiling cython extensions.
-
-If you are using a linux system you may need to install some dependencies for LavaVu.
 
-.. code-block::
-
-    sudo apt-get update  && sudo apt-get install python3 python3-venv python3-dev make pybind11-dev mesa-common-dev mesa-utils libgl1-mesa-dev gcc g++
 
 
 
@@ -35,7 +33,7 @@ If you are using a linux system you may need to install some dependencies for La
 
     git clone https://github.com/Loop3D/LoopStructural.git
     cd LoopStructural
-    pip install .
+    pip install -e . # -e installs as an editable package so you can modify the source code and see the changes immediately
 
 Dependencies
 ~~~~~~~~~~~~
@@ -51,18 +49,13 @@ Required dependencies:
 Optional dependencies:
 
 * matplotlib, 2D/3D visualisation
-* LavaVu, 3D visualisation
+* pyvista, 3D visualisation
 * surfepy, radial basis interpolation
 * map2loop, generation of input datasets from regional Australian maps
+* geoh5py, export to gocad hdf5 format
+* pyevtk, export to vtk format
+* dill, serialisation of python objects
+* loopsolver, solving of inequalities
+* tqdm, progress bar
 
 
-Docker
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-LoopStructural can be used either by compiling the docker image or by pulling the compiled
-docker image from docker hub.
-
-.. code-block::
-
-    docker pull loop3d/loopstructural
-    docker run -i -t -p 8888:8888 -v LOCALDIRPATH:/home/jovyan/shared_volume loop3d/loopstructural`.
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 480542b05..083f0d1a3 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -3,29 +3,11 @@
    You can adapt this file completely to your liking, but it should at least
    contain the root `toctree` directive.
 
-Loop Structural
+LoopStructural
 ===============
 .. image:: ./images/image823.png
 
-.. important:: **Upgrading from 1.5.x to 1.6.x**
 
-   Version 1.6.x of LoopStructural has been released and there are a number of changes to how the library is used. The main changes are:
-   
-   - The `LavaVuModelViewer` class has been removed and the model visualisation has been ported to `pyvista` using a LoopStructural 
-     wrapped `loopstructuralvisualisation`. Many of the functions have been re-named and can be combined with the `pyvista.Plotter` methods for advanced
-     visualisation. pyvista does not currently allow for interactive plots while using colab notebooks, so `jupyter_backend='static'` is needed.
-   - Datastructures for representing model input/output have been introduced including the `LoopStructural.datatypes.StructuredGrid`,
-     :class:`LoopStructural.datatypes.ValuePoints`, :class:`LoopStructural.datatypes.VectorPoints` and :class:`LoopStructural.datatypes.Surface`. These Datastructures
-     are intended to improve the accessibility of the modelling and allow for easier export into other software packages.
-     All data objects have a `.save` method that allows the user to save the points into various formats including `json`, `csv`, `vtk` and `geoh5`.
-     In the future exporters for gocad and other formats will be added. All objects also have a `.vtk()` method to return a `pyvista.PolyData` object.
-   - The solver options for the interpolators have been removed and only :meth:`scipy.sparse.linalg.cg` and :meth:`scipy.sparse.linalg.lsmr` are supported. 
-     To use an external solver the `external_solver` argument can be provided to the `interpolator.solve` method, usually passed through
-     the `model.create_and_add*` methods. To specify specific arguments to the solver these can be passed through `solver_kwargs` dict.
-   - Improvements have been made to how unconformities are handled and their interaction with faults
-   - A framework for including inequalities into the interpolation has been added but there is no solver that can be used for solving this problem.
-   - A convenience class :class:`LoopStructural.LoopInterpolator` was added for building an implicit function without requiring the :class:`LoopStructural.GeologicalModel`
-   **Please report any bugs or issues on the github repository.**
  
 Overview
 ========
@@ -38,9 +20,113 @@ Loop is an open source 3D probabilistic geological and geophysical modelling pla
 initiated by Geoscience Australia and the OneGeology consortium. The project is funded by Australian territory,
 State and Federal Geological Surveys, the Australian Research Council and the MinEx Collaborative Research Centre.
 
+Examples
+==========
 
+LoopStructural can be used to build a 3D geological model using geological relationships between geological objects
+e.g. faults, folds, unconformities and stratigraphic contacts. The library also provides a high level API to access
+the fast interpolation schemes that are used by LoopStructural.
+
+Using GeologicalModel
+----------------------
+
+The following example shows how to use the geological model interface to create a geological model from a dataset and 
+evaluate the scalar field and gradient of the interpolator at some random locations.
+
+.. pyvista-plot::
+   :context:
+   :include-source: true
+   :force_static:
+
+    from LoopStructural import GeologicalModel
+    from LoopStructural.datatypes import BoundingBox
+    from LoopStructural.visualisation import Loop3DView
+    from LoopStructural.datasets import load_claudius
+    
+    import numpy as np
+    data, bb = load_claudius()
+    
+    #bb constaints origin and maximum of axis aligned bounding box
+    #data is a pandas dataframe with X,Y,Z,val,nx,ny,nz, feature_name
+    
+    model = GeologicalModel(bb[0,:],bb[1,:])
+    model.data = data
+    # nelements specifies the number of discrete interpolation elements
+    # 'stratí' is the feature name in the data dataframe
+    model.create_and_add_foliation('strati',nelements=1e5)
+    model.update()
+    # get the value of the interpolator at some random locations
+    locations = np.array(
+        [
+            np.random.uniform(bb[0, 0], bb[1, 0],5),
+            np.random.uniform(bb[0, 1], bb[1, 1],5),
+            np.random.uniform(bb[0, 2], bb[1, 2],5),
+        ]
+    ).T
+    val = model.evaluate_feature_value('strati', locations)
+    # get the gradient of the interpolator
+    gradient = model.evaluate_feature_gradient('strati',locations)
+    
+    #Plot the scalar field of the model
+    model['strati'].scalar_field().plot()
+
+
+Using InterpolatorBuilder
+-------------------------
+
+To access the interpolators directly the `InterpolatorBuilder` can be used
+to help assemble an interpolator from a combination of value, gradient and inequality 
+constraints.
+
+.. pyvista-plot::
+   :context:
+   :include-source: true
+   :force_static:
+
+    from LoopStructural import BoundingBox, InterpolatorBuilder
+    from LoopStructural.utils import EuclideanTransformation
+    from LoopStructural.datasets import load_claudius
+
+    # load in a dataframe with x,y,z,val,nx,ny,nz to constrain the value and
+    # gradient normal of the interpolator
+    data, bb = load_claudius()
+    # create a transformer to move the data to the centre of mass of the dataset
+    # this avoid any numerical issues caused by large coordinates. This dataset
+    # is already axis aligned so we don't need to rotate the data but we can rotate it
+    # so that the main anisotropy of the dataset is aligned with the x axis.
+    transformer = EuclideanTransformation(dimensions=3, fit_rotation=False).fit(
+        data[["X", "Y", "Z"]]
+    )
+    data[["X", "Y", "Z"]] = transformer.transform(data[["X", "Y", "Z"]])
+    # Find the bounding box of the data by finding the extent
+    bounding_box = BoundingBox().fit(data[["X", "Y", "Z"]])
+    # assemble an interpolator using the bounding box and the data
+    interpolator = (
+        InterpolatorBuilder(interpolatortype="FDI", bounding_box=bounding_box)
+        .add_value_constraints(data.loc[data["val"].notna(), ["X", "Y", "Z", "val"]].values)
+        .add_normal_constraints(
+            data.loc[data["nx"].notna(), ["X", "Y", "Z", "nx", "ny", "nz"]].values
+        ).build()
+    )
+    # Set the number of elements in the bounding box to 10000 and create a structured grid
+    bounding_box.nelements = 10000
+    mesh = bounding_box.structured_grid()
+    # add the interpolated values to the mesh at the nodes
+    # mesh.properties["v"] = interpolator.evaluate_value(bounding_box.regular_grid(order='F'))
+    mesh.properties["v"] = interpolator.evaluate_value(mesh.nodes)
+
+    # or the cell centres
+    # mesh.cell_properties["v"] = interpolator.evaluate_value(mesh.cell_centres)
+
+
+    # visualise the scalar value
+
+    mesh.plot()
+
+    # We can also add gradient properties and visualise these
+    mesh.properties["grad"] = interpolator.evaluate_gradient(mesh.nodes)
+    mesh.vtk().glyph(orient="grad").plot()
 
-   
 .. toctree::
    :hidden:
 
diff --git a/docs/source/user_guide/index.rst b/docs/source/user_guide/index.rst
index 1c209014f..cccc80a91 100644
--- a/docs/source/user_guide/index.rst
+++ b/docs/source/user_guide/index.rst
@@ -1,9 +1,15 @@
 User Guide
 ===========
+This section of the documentation provides a detailed guide on how to use LoopStructural for users who may or may not be familiar with
+implicit 3D geological modelling, but are looking to utilise LoopStructural for their geological modelling needs. See the pane on the left
+for a list of topics covered in this section.
+
+
 
 .. toctree::
    :caption: User Guide
-
+   
+   what_is_a_geological_model
    input_data
    geological_model
    interpolation_options
diff --git a/docs/source/user_guide/what_is_a_geological_model.rst b/docs/source/user_guide/what_is_a_geological_model.rst
new file mode 100644
index 000000000..7e07fae26
--- /dev/null
+++ b/docs/source/user_guide/what_is_a_geological_model.rst
@@ -0,0 +1,197 @@
+Implicit Geological Modelling
+-----------------------------
+
+
+Implicit geological modelling is an approach for representing objects in a geological model using mathematical function :math:`f(x,y,z)` where the function
+returns the distance to the geological feature. Implicit modelling involves approximating the geometry of the geological object by using observations 
+of the location and shape of the geological objects. 
+
+Building a basic implicit function
+=======================================
+
+For example, if we have a fault on a map. The geometry of this fault is defined by a trace on the map. If we want to approximate an implicit function that
+represents this geometry then we would need to find the function :math:`f(x,y,z)` where :math:`f(x,y,z) = 0`` at the observations of the fault trace. This function
+is under constrained and the best solution would be a function that is 0 everywhere. To further constrain this function we either need to add additional 
+value constraints that are not 0 or a constraint to the gradient norm of the function. In a geological context, the gradient norm would be representative
+of the strike/dip of a geological object with some sense of younging (or a consistent direction for all observations).
+
+.. pyvista-plot::
+   :context:
+   
+    import numpy as np
+    import pandas as pd
+    from LoopStructural.visualisation import Loop3DView
+    x = np.linspace(0, 1, 10)
+    y = np.zeros(10)
+    z = np.ones(10)
+    v = np.zeros(10)
+
+    data = pd.DataFrame({"X": x, "Y": y, "Z": z, "val": v})
+    vector = pd.DataFrame(
+        [[0.5, 0, 1.0, 0, 1, 0]], columns=["X", "Y", "Z", "nx", "ny", "nz"]
+    )
+
+    view = Loop3DView()
+    view.add_points(data[["X", "Y", "Z"]].values)
+    view.add_arrows(
+        vector[["X", "Y", "Z"]].values,
+        vector[["nx", "ny", "nz"]].values,
+    )
+    view.show()
+
+To fit an implicit function using LoopStructural we can use the `InterpolatorBuilder` to create an interpolator that
+fits the value and orientation data as gradient norm constraints.
+The interpolator can be evaluated on a grid to visualise the geometry of the geological object. 
+
+.. pyvista-plot::
+   :context:
+
+    import numpy as np
+    import pandas as pd
+    from LoopStructural import InterpolatorBuilder, BoundingBox
+    from LoopStructural.visualisation import Loop3DView
+    
+    bounding_box = BoundingBox(np.array([-2, -2, -2]), np.array([2, 2, 2]))
+
+    interpolator = (
+        InterpolatorBuilder(
+            interpolatortype="FDI", nelements=1e4, bounding_box=bounding_box
+        ).add_value_constraints(data.values)
+        .add_normal_constraints(vector.values).setup_interpolator()
+        .build()
+    )
+    mesh = bounding_box.structured_grid()
+    mesh.properties['val'] = interpolator.evaluate_value(mesh.nodes)
+    view = Loop3DView()
+    view.add_mesh(mesh.vtk())
+    view.show()
+
+The resulting scalar field measures the distance to a reference horizon from -2 to 2. 
+
+We can extract an isosurface of the scalar field at 0 to visualise the geometry of the fault.
+
+.. pyvista-plot::
+   :context:
+
+    view = Loop3DView()
+    view.add_mesh(mesh.vtk().contour([0]))
+    view.add_points(data[["X", "Y", "Z"]].values)
+    view.add_arrows(
+        vector[["X", "Y", "Z"]].values,
+        vector[["nx", "ny", "nz"]].values,
+    )
+    view.show()
+
+The resulting isosurface represents the geometry of the fault. The norm of the gradient of the scalar field 
+controls how quickly we increase or decrease the scalar field value. For a single geological feature this
+is not particularly important but when we have multiple features represented by a single implicit function
+this can be important. Below shows the isovalues of -1,0,1 of the scalar field. 
+
+.. pyvista-plot::
+   :context:
+
+    view = Loop3DView()
+    view.add_mesh(mesh.vtk().contour([-1, 0, 1]))
+    view.remove_scalar_bar()
+    view.add_points(data[["X", "Y", "Z"]].values)
+    view.add_arrows(
+        vector[["X", "Y", "Z"]].values,
+        vector[["nx", "ny", "nz"]].values,
+    )
+    view.remove_scalar_bar()
+    view.view_yz()
+    view.show()
+
+We can see that the length of the vector (unit vector) coincides with the distance between the two surfaces. 
+If we were to change the length of the vector we would change the distance between the two surfaces.
+
+.. pyvista-plot::
+   :context:
+
+   vector[["nx", "ny", "nz"]] = vector[["nx", "ny", "nz"]]*2
+    bounding_box = BoundingBox(np.array([-2, -2, -2]), np.array([2, 2, 2]))
+
+    interpolator = (
+        InterpolatorBuilder(
+            interpolatortype="FDI", nelements=1e4, bounding_box=bounding_box
+        ).add_value_constraints(data.values)
+        .add_normal_constraints(vector.values).setup_interpolator()
+        .build()
+    )
+    mesh2 = bounding_box.structured_grid()
+    mesh2.properties['val'] = interpolator.evaluate_value(mesh2.nodes)
+    view = Loop3DView()
+    view.add_mesh(mesh2.vtk().contour([-1, 0, 1]),colour='red')
+    view.add_mesh(mesh.vtk().contour([-1, 0, 1]),colour='blue')
+    view.add_points(data[["X", "Y", "Z"]].values)
+    view.add_arrows(
+        vector[["X", "Y", "Z"]].values,
+        vector[["nx", "ny", "nz"]].values,
+    )
+    view.remove_scalar_bar()
+    view.view_yz()
+    view.show()
+
+
+This time the vector is twice as long and the distance between the surfaces is half. This is because we have specified that the gradient of the 
+scalar field has a magnitude of 2 which means that it grows twice as quickly. 
+
+Modelling with only value constraints
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+An alternative approach for constraining the scalar field is to use only value constraints. This is useful when modelling a stratigraphic series
+where we may try to interpolate the distance to a contact and use the cummulative thickness between the stratigraphic units to constrain the value 
+of the implicit function. 
+
+Following the example above we will use two lines of points with a value of 0 and 0.5 to represent two contacts between stratigraphic units.
+
+.. pyvista-plot::
+   :context:
+
+    import numpy as np
+    import pandas as pd
+    from LoopStructural import InterpolatorBuilder, BoundingBox
+    from LoopStructural.visualisation import Loop3DView
+
+    x = np.linspace(0, 1, 10)
+    y = np.zeros(10)
+    z = np.ones(10)
+    v = np.zeros(10)
+
+    data = pd.concat(
+        [
+            pd.DataFrame({"X": x, "Y": y, "Z": z, "val": v}),
+            pd.DataFrame({"X": x, "Y": y + 0.5, "Z": z, "val": v + 0.5}),
+        ]
+    )
+
+    bounding_box = BoundingBox(np.array([-2, -2, -2]), np.array([2, 2, 2]))
+
+    interpolator = (
+        InterpolatorBuilder(
+            interpolatortype="FDI", nelements=1e4, bounding_box=bounding_box
+        )
+        .add_value_constraints(data.values)
+        .setup_interpolator()
+        .build()
+    )
+    mesh = bounding_box.structured_grid()
+    mesh.properties["val"] = interpolator.evaluate_value(mesh.nodes)
+    view = Loop3DView()
+    view.add_mesh(mesh.vtk().contour([0, 0.5]), opacity=0.4)
+    view.remove_scalar_bar()
+    view.add_points(data[["X", "Y", "Z"]].values, scalars=data["val"].values)
+    view.remove_scalar_bar()
+
+    view.view_yz()
+    view.show()
+
+We can see that the scalar field fits the value constraints and interpolates between the two surfaces.
+The scalar field value constraints can have a significant impact on the geometry of the implicit function, especially
+where the the geometry of the contacts outlies a structure (e.g. folded layers). Larger differences between the scalar field
+value effectively increase the magnitude of the implicit functions gradient norm. 
+
+
+
+
+
diff --git a/pyproject.toml b/pyproject.toml
index 6df5be13b..6f9b40951 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -37,28 +37,31 @@ dependencies = [
     "scipy",
     "scikit-image",
     "scikit-learn",
-    "tqdm",
 ]
 dynamic = ['version']
 
 [project.optional-dependencies]
-all = ['loopstructural[visualisation,inequalities,export]']
-visualisation = [
-    "matplotlib",
-    "pyvista",
-    "loopstructuralviusualisation>0.1.4"
-    ]
-export = [
-    "geoh5py", 
-    "pyevtk", 
-    "dill"]
-jupyter = [
+all = ['loopstructural[visualisation,inequalities,export]', 'tqdm']
+visualisation = ["matplotlib", "pyvista", "loopstructuralviusualisation>0.1.4"]
+export = ["geoh5py", "pyevtk", "dill"]
+jupyter = ["pyvista[all]"]
+inequalities = ["loopsolver"]
+docs = [
     "pyvista[all]",
-    "tqdm"
-    ]
-inequalities = [
-"loopsolver"]
-docs = ["pyvista[all]", "pydata-sphinx-theme","meshio","loopstructuralvisualisation[all]","networkx","scikit-learn","scikit-image","sphinx","sphinx-gallery","geoh5py","geopandas","sphinxcontrib-bibtex","myst-parser"]
+    "pydata-sphinx-theme",
+    "meshio",
+    "loopstructuralvisualisation[all]",
+    "networkx",
+    "scikit-learn",
+    "scikit-image",
+    "sphinx",
+    "sphinx-gallery",
+    "geoh5py",
+    "geopandas",
+    "sphinxcontrib-bibtex",
+    "myst-parser",
+    "sphinx-design",
+]
 [project.urls]
 Documentation = 'https://Loop3d.org/LoopStructural/'
 "Bug Tracker" = 'https://github.com/loop3d/loopstructural/issues'
diff --git a/tests/unit/interpolator/test_interpolator_builder.py b/tests/unit/interpolator/test_interpolator_builder.py
index ee12de1b2..bd0f46155 100644
--- a/tests/unit/interpolator/test_interpolator_builder.py
+++ b/tests/unit/interpolator/test_interpolator_builder.py
@@ -22,56 +22,54 @@ def setup_builder():
 
 def test_create_interpolator(setup_builder):
     builder = setup_builder
-    builder.create_interpolator()
+    builder.build()
     assert builder.interpolator is not None, "Interpolator should be created"
 
 
 def test_set_value_constraints(setup_builder):
     builder = setup_builder
-    builder.create_interpolator()
-    value_constraints = np.array([[0.5, 0.5, 0.5, 1.0,1.]])
-    builder.set_value_constraints(value_constraints)
+    builder.build()
+    value_constraints = np.array([[0.5, 0.5, 0.5, 1.0, 1.0]])
+    builder.add_value_constraints(value_constraints)
     assert np.array_equal(
-        builder.interpolator.data['value'], value_constraints
+        builder.interpolator.data["value"], value_constraints
     ), "Value constraints should be set correctly"
 
 
 def test_set_gradient_constraints(setup_builder):
     builder = setup_builder
-    builder.create_interpolator()
-    gradient_constraints = np.array([[0.5, 0.5, 0.5, 1.0, 0.0, 0.0,1.]])
-    builder.set_gradient_constraints(gradient_constraints)
+    gradient_constraints = np.array([[0.5, 0.5, 0.5, 1.0, 0.0, 0.0, 1.0]])
+    builder.add_gradient_constraints(gradient_constraints)
     assert np.array_equal(
-        builder.interpolator.data['gradient'], gradient_constraints
+        builder.interpolator.data["gradient"], gradient_constraints
     ), "Gradient constraints should be set correctly"
 
 
 def test_set_normal_constraints(setup_builder):
     builder = setup_builder
-    builder.create_interpolator()
-    normal_constraints = np.array([[0.5, 0.5, 0.5, 1.0, 0.0, 0.0,1.]])
-    builder.set_normal_constraints(normal_constraints)
+    normal_constraints = np.array([[0.5, 0.5, 0.5, 1.0, 0.0, 0.0, 1.0]])
+    builder.add_normal_constraints(normal_constraints)
     assert np.array_equal(
-        builder.interpolator.data['normal'], normal_constraints
+        builder.interpolator.data["normal"], normal_constraints
     ), "Normal constraints should be set correctly"
 
 
 def test_setup_interpolator(setup_builder):
     builder = setup_builder
-    builder.create_interpolator()
-    value_constraints = np.array([[0.5, 0.5, 0.5, 1.0,1.]])
-    interpolator = builder.set_value_constraints(value_constraints).setup_interpolator()
+    builder.build()
+    value_constraints = np.array([[0.5, 0.5, 0.5, 1.0, 1.0]])
+    interpolator = builder.add_value_constraints(value_constraints).setup_interpolator().build()
     assert interpolator is not None, "Interpolator should be set up"
     assert np.array_equal(
-        interpolator.data['value'], value_constraints
+        interpolator.data["value"], value_constraints
     ), "Value constraints should be set correctly after setup"
 
 
 def test_evaluate_scalar_value(setup_builder):
     builder = setup_builder
-    builder.create_interpolator()
+    builder.build()
     value_constraints = np.array([[0.5, 0.5, 0.5, 1.0]])
-    interpolator = builder.set_value_constraints(value_constraints).setup_interpolator()
+    interpolator = builder.add_value_constraints(value_constraints).setup_interpolator().build()
     locations = np.array([[0.5, 0.5, 0.5]])
     values = interpolator.evaluate_value(locations)
     assert values is not None, "Evaluation should return values"
diff --git a/tests/unit/modelling/intrusions/test_intrusions.py b/tests/unit/modelling/intrusions/test_intrusions.py
index 3f7e62385..377ccc2e5 100644
--- a/tests/unit/modelling/intrusions/test_intrusions.py
+++ b/tests/unit/modelling/intrusions/test_intrusions.py
@@ -115,9 +115,11 @@ def test_intrusion_builder():
 
     intrusion_builder.set_data_for_extent_calculation(intrusion_data)
 
-    intrusion_builder.build_arguments = {
-        "geometric_scaling_parameters": {},
-    }
+    intrusion_builder.update_build_arguments(
+        {
+            "geometric_scaling_parameters": {},
+        }
+    )
 
     intrusion_builder.update()
 
diff --git a/tests/unit/modelling/test_structural_frame.py b/tests/unit/modelling/test_structural_frame.py
index c9dea67fa..7e51ecde4 100644
--- a/tests/unit/modelling/test_structural_frame.py
+++ b/tests/unit/modelling/test_structural_frame.py
@@ -67,13 +67,45 @@ def test_create_structural_frame_pli():
     model.update()
 
     assert np.all(
-        np.isclose(fault[0].evaluate_gradient(np.array([[5, 5, 5]])), [0, 0, 1], atol=1e-4)
+        np.isclose(fault[0].evaluate_gradient(np.array([[5, 5, 5]])), [0, 0, 1], atol=1e-1)
     )
     assert np.all(
-        np.isclose(fault[1].evaluate_gradient(np.array([[5, 5, 5]])), [0, 1, 0], atol=1e-4)
+        np.isclose(fault[1].evaluate_gradient(np.array([[5, 5, 5]])), [0, 1, 0], atol=1e-1)
     )
     assert np.all(
-        np.isclose(fault[2].evaluate_gradient(np.array([[5, 5, 5]])), [1, 0, 0], atol=1e-4)
+        np.isclose(fault[2].evaluate_gradient(np.array([[5, 5, 5]])), [1, 0, 0], atol=1e-1)
+    )
+
+
+def test_create_structural_frame_fdi():
+    data = pd.DataFrame(
+        [
+            [5.1, 5.1, 5, 0, 0, 1, 0, 0],
+            [5, 5.1, 5, 0, 1, 0, 1, 0],
+            [5.1, 5, 5, 1, 0, 0, 2, 0],
+        ],
+        columns=["X", "Y", "Z", "nx", "ny", "nz", "coord", "val"],
+    )
+    data["feature_name"] = "fault"
+
+    bb = BoundingBox(origin=np.zeros(3), maximum=np.ones(3) * 10)
+    model = GeologicalModel(bb.origin, bb.maximum)
+
+    model.data = data
+
+    fault = model.create_and_add_fault(
+        "fault", 10, nelements=2000, steps=4, interpolatortype="FDI", buffer=2
+    )
+    model.update()
+
+    assert np.all(
+        np.isclose(fault[0].evaluate_gradient(np.array([[5, 5, 5]])), [0, 0, 1], atol=1e-1)
+    )
+    assert np.all(
+        np.isclose(fault[1].evaluate_gradient(np.array([[5, 5, 5]])), [0, 1, 0], atol=1e-1)
+    )
+    assert np.all(
+        np.isclose(fault[2].evaluate_gradient(np.array([[5, 5, 5]])), [1, 0, 0], atol=1e-1)
     )