Various small fixes

.github/workflows/documentation.yml

@@ -1,4 +1,4 @@
-name: Build documentation and deploy
+name: "📚 Build documentation and deploy "
 
 on:
   workflow_dispatch:

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)
+        )

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,25 +126,48 @@ 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)
             vectors[norm > 0, :] /= norm[norm > 0][:, None]
         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

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]

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

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)

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
 

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

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,
+        )

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

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
+        )

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