fix: updating scaling for plotting

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
 

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
 

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

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

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