Fix vector point visualisation

LoopStructural/datatypes/_bounding_box.py

@@ -326,9 +326,7 @@ def get_value(self, name):
         if iy == -1:
             return self.origin[ix]
 
-        return self.bb[
-            ix,
-        ]
+        return self.bb[ix,]
 
     def __getitem__(self, name):
         if isinstance(name, str):

LoopStructural/datatypes/_point.py

@@ -147,8 +147,8 @@ def vtk(
         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]
+            norm[norm > 0] = norm[norm > 0] / norm[norm > 0].max()
+            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]

LoopStructural/interpolators/_finite_difference_interpolator.py

@@ -272,7 +272,12 @@ def add_gradient_constraints(self, w=1.0):
             idc[inside, :] = gi[node_idx[inside, :]]
             inside = np.logical_and(~np.any(idc == -1, axis=1), inside)
 
-            (vertices, T, elements, inside_,) = self.support.get_element_gradient_for_location(
+            (
+                vertices,
+                T,
+                elements,
+                inside_,
+            ) = self.support.get_element_gradient_for_location(
                 points[inside, : self.support.dimension]
             )
             # normalise constraint vector and scale element matrix by this
@@ -335,7 +340,12 @@ def add_norm_constraints(self, w=1.0):
             # calculate unit vector for node gradients
             # this means we are only constraining direction of grad not the
             # magnitude
-            (vertices, T, elements, inside_,) = self.support.get_element_gradient_for_location(
+            (
+                vertices,
+                T,
+                elements,
+                inside_,
+            ) = self.support.get_element_gradient_for_location(
                 points[inside, : self.support.dimension]
             )
             # T*=np.product(self.support.step_vector)
@@ -422,7 +432,12 @@ def add_gradient_orthogonal_constraints(
             vectors[norm > 0, :] /= norm[norm > 0, None]
 
             # normalise element vector to unit vector for dot product
-            (vertices, T, elements, inside_,) = self.support.get_element_gradient_for_location(
+            (
+                vertices,
+                T,
+                elements,
+                inside_,
+            ) = self.support.get_element_gradient_for_location(
                 points[inside, : self.support.dimension]
             )
             T[norm > 0, :, :] /= norm[norm > 0, None, None]

LoopStructural/interpolators/supports/_3d_base_structured.py

@@ -162,9 +162,9 @@ def origin(self, origin):
         length = self.maximum - origin
         length /= self.step_vector
         self._nsteps = np.ceil(length).astype(np.int64)
-        self._nsteps[
-            self._nsteps == 0
-        ] = 3  # need to have a minimum of 3 elements to apply the finite difference mask
+        self._nsteps[self._nsteps == 0] = (
+            3  # need to have a minimum of 3 elements to apply the finite difference mask
+        )
         if np.any(~(self._nsteps > 0)):
             logger.error(
                 f"Cannot resize the interpolation support. The proposed number of steps is {self._nsteps}, these must be all > 0"

LoopStructural/interpolators/supports/_3d_structured_tetra.py

@@ -166,9 +166,9 @@ def _init_face_table(self):
         shared_face_index[:] = -1
         shared_face_index[row.reshape(-1, 3)[:, 0], :] = col.reshape(-1, 3)
 
-        self.shared_elements[
-            np.arange(self.shared_element_relationships.shape[0]), :
-        ] = shared_face_index
+        self.shared_elements[np.arange(self.shared_element_relationships.shape[0]), :] = (
+            shared_face_index
+        )
         # resize
         self.shared_elements = self.shared_elements[: len(self.shared_element_relationships), :]
 

LoopStructural/interpolators/supports/_3d_unstructured_tetra.py

@@ -173,9 +173,9 @@ def _init_face_table(self):
         shared_face_index[:] = -1
         shared_face_index[row.reshape(-1, 3)[:, 0], :] = col.reshape(-1, 3)
 
-        self.shared_elements[
-            np.arange(self.shared_element_relationships.shape[0]), :
-        ] = shared_face_index
+        self.shared_elements[np.arange(self.shared_element_relationships.shape[0]), :] = (
+            shared_face_index
+        )
         # resize
         self.shared_elements = self.shared_elements[: len(self.shared_element_relationships), :]
         # flag = np.zeros(self.elements.shape[0])

LoopStructural/interpolators/supports/_face_table.py

@@ -63,8 +63,8 @@ def _init_face_table(grid):
     shared_face_index = np.zeros((shared_faces.shape[0], grid.dimension), dtype=int)
     shared_face_index[:] = -1
     shared_face_index[row.reshape(-1, grid.dimension)[:, 0], :] = col.reshape(-1, grid.dimension)
-    grid._shared_elements[
-        np.arange(grid.shared_element_relationships.shape[0]), :
-    ] = shared_face_index
+    grid._shared_elements[np.arange(grid.shared_element_relationships.shape[0]), :] = (
+        shared_face_index
+    )
     # resize
     grid._shared_elements = grid.shared_elements[: len(grid.shared_element_relationships), :]

LoopStructural/modelling/core/geological_model.py

@@ -598,12 +598,10 @@ def data(self, data: pd.DataFrame):
                 * self._data.loc[mask, "polarity"].to_numpy()[:, None]
             )
             self._data.drop(["strike", "dip"], axis=1, inplace=True)
-        self._data[
-            ['X', 'Y', 'Z', 'val', 'nx', 'ny', 'nz', 'gx', 'gy', 'gz', 'tx', 'ty', 'tz']
-        ] = self._data[
-            ['X', 'Y', 'Z', 'val', 'nx', 'ny', 'nz', 'gx', 'gy', 'gz', 'tx', 'ty', 'tz']
-        ].astype(
-            float
+        self._data[['X', 'Y', 'Z', 'val', 'nx', 'ny', 'nz', 'gx', 'gy', 'gz', 'tx', 'ty', 'tz']] = (
+            self._data[
+                ['X', 'Y', 'Z', 'val', 'nx', 'ny', 'nz', 'gx', 'gy', 'gz', 'tx', 'ty', 'tz']
+            ].astype(float)
         )
 
     def set_model_data(self, data):
@@ -1805,7 +1803,9 @@ def get_stratigraphic_surfaces(self, units: List[str] = [], bottoms: bool = True
                     values.to_list(),
                     self.bounding_box,
                     name=names.loc[values.index].to_list(),
-                    colours=unit_table.loc[unit_table['feature_name'] == u, 'colour'].tolist()[1:], #we don't isosurface basement, no value
+                    colours=unit_table.loc[unit_table['feature_name'] == u, 'colour'].tolist()[
+                        1:
+                    ],  # we don't isosurface basement, no value
                 )
             )
 

LoopStructural/modelling/features/fault/_fault_function.py

@@ -9,6 +9,10 @@
 logger = getLogger(__name__)
 
 
+def smooth_peak(x):
+    return 0.25 * x**6 + 0.5 * x**4 - 1.75 * x**2 + 1
+
+
 class FaultProfileFunction(metaclass=ABCMeta):
     def __init__(self):
         self.lim = [-1, 1]
@@ -412,14 +416,7 @@ class BaseFault(object):
     # gyf.add_min(-1)
     # gyf.add_max(1)
     gyf = Ones()
-    gzf = CubicFunction()
-    gzf.add_cstr(-1, 0)
-    gzf.add_cstr(1, 0)
-    gzf.add_cstr(-0.2, 1)
-    gzf.add_cstr(0.2, 1)
-    gzf.add_grad(0, 0)
-    gzf.add_min(-1)
-    gzf.add_max(1)
+    gzf = smooth_peak
     gxf = Composite(hw, fw)
     fault_displacement = FaultDisplacement(gx=gxf, gy=gyf, gz=gzf)
 
@@ -441,22 +438,17 @@ class BaseFault3D(object):
     fw.add_cstr(-1, 0)
     fw.add_grad(-1, 0)
     fw.add_min(-1)
-    gyf = CubicFunction()
-    gyf.add_cstr(-1, 0)
-    gyf.add_cstr(1, 0)
-    gyf.add_cstr(-0.2, 1)
-    gyf.add_cstr(0.2, 1)
-    gyf.add_grad(0, 0)
-    gyf.add_min(-1)
-    gyf.add_max(1)
+
+    gyf = smooth_peak
+    # CubicFunction()
+    # gyf.add_cstr(-1, 0)
+    # gyf.add_cstr(1, 0)
+    # gyf.add_cstr(-0.2, 1)
+    # gyf.add_cstr(0.2, 1)
+    # gyf.add_grad(0, 0)
+    # gyf.add_min(-1)
+    # gyf.add_max(1)
     # gyf = Ones()
-    gzf = CubicFunction()
-    gzf.add_cstr(-1, 0)
-    gzf.add_cstr(1, 0)
-    gzf.add_cstr(-0.2, 1)
-    gzf.add_cstr(0.2, 1)
-    gzf.add_grad(0, 0)
-    gzf.add_min(-1)
-    gzf.add_max(1)
+    gzf = smooth_peak
     gxf = Composite(hw, fw)
     fault_displacement = FaultDisplacement(gx=gxf, gy=gyf, gz=gzf)

LoopStructural/modelling/features/fault/_fault_segment.py

@@ -125,11 +125,17 @@ def fault_ellipsoid(self, **kwargs):
         try:
             import pyvista as pv
 
-            fault_ellipsoid = pv.PolyData(
-                self.model.rescale(self.fault_centre[None, :], inplace=False)
-            )
-            fault_ellipsoid["norm"] = self.builder.fault_normal_vector[None, :]
-
+            if self.model is None:
+                pts = self.fault_centre[None, :]
+            else:
+                pts = self.model.rescale(self.fault_centre[None, :], inplace=False)
+            # pts = self.fault_centre[None, :]
+            fault_ellipsoid = pv.PolyData(pts)
+            # fault_ellipsoid = pv.PolyData(
+            #     self.model.rescale(self.fault_centre[None, :], inplace=False)
+            # )
+            fault_ellipsoid["norm"] = self.fault_normal_vector[None, :]
+            fault_ellipsoid['norm'] /= np.linalg.norm(fault_ellipsoid['norm'], axis=1)[:, None]
             geom = pv.ParametricEllipsoid(
                 self.fault_minor_axis,
                 self.fault_major_axis,

LoopStructural/modelling/input/process_data.py

@@ -299,9 +299,9 @@ def fault_properties(self, fault_properties):
             pts = self.fault_locations.loc[
                 self.fault_locations["feature_name"] == fname, ["X", "Y", "Z"]
             ]
-            fault_properties.loc[
-                fname, ["centreEasting", "centreNorthing", "centreAltitude"]
-            ] = np.nanmean(pts, axis=0)
+            fault_properties.loc[fname, ["centreEasting", "centreNorthing", "centreAltitude"]] = (
+                np.nanmean(pts, axis=0)
+            )
         if (
             "avgNormalEasting" not in fault_properties.columns
             or "avgNormalNorthing" not in fault_properties.columns
@@ -449,9 +449,9 @@ def _stratigraphic_value(self):
         for _name, sg in self.stratigraphic_order:
             value = 0.0  # reset for each supergroup
             if sg[0] not in self.thicknesses or self.thicknesses[sg[0]] <= 0:
-                self.thicknesses[
-                    sg[0]
-                ] = np.inf  # make the top unit infinite as it should extend to the top of the model
+                self.thicknesses[sg[0]] = (
+                    np.inf
+                )  # make the top unit infinite as it should extend to the top of the model
             for g in reversed(
                 sg[:-1]
             ):  # don't add the last unit as we never see the base of this unit.

tests/unit/datatypes/test__surface.py

@@ -10,7 +10,9 @@ def test_surface_creation():
     name = "surface1"
     values = np.array([1, 2, 3])
 
-    surface = Surface(vertices=vertices, triangles=triangles, normals=normals, name=name, values=values)
+    surface = Surface(
+        vertices=vertices, triangles=triangles, normals=normals, name=name, values=values
+    )
 
     assert np.array_equal(surface.vertices, vertices)
     assert np.array_equal(surface.triangles, triangles)