diff --git a/LoopStructural/datatypes/_bounding_box.py b/LoopStructural/datatypes/_bounding_box.py index 5b6f662d..dc1d01e8 100644 --- a/LoopStructural/datatypes/_bounding_box.py +++ b/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): diff --git a/LoopStructural/datatypes/_point.py b/LoopStructural/datatypes/_point.py index c7829b2a..4542262d 100644 --- a/LoopStructural/datatypes/_point.py +++ b/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] diff --git a/LoopStructural/interpolators/_finite_difference_interpolator.py b/LoopStructural/interpolators/_finite_difference_interpolator.py index 1a169747..a0b03e11 100644 --- a/LoopStructural/interpolators/_finite_difference_interpolator.py +++ b/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] diff --git a/LoopStructural/interpolators/supports/_3d_base_structured.py b/LoopStructural/interpolators/supports/_3d_base_structured.py index 74a48e5a..4c91322c 100644 --- a/LoopStructural/interpolators/supports/_3d_base_structured.py +++ b/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" diff --git a/LoopStructural/interpolators/supports/_3d_structured_tetra.py b/LoopStructural/interpolators/supports/_3d_structured_tetra.py index 052d6397..70ecba37 100644 --- a/LoopStructural/interpolators/supports/_3d_structured_tetra.py +++ b/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), :] diff --git a/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py b/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py index 7ecae726..f55e8deb 100644 --- a/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py +++ b/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]) diff --git a/LoopStructural/interpolators/supports/_face_table.py b/LoopStructural/interpolators/supports/_face_table.py index fa818a53..40746719 100644 --- a/LoopStructural/interpolators/supports/_face_table.py +++ b/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), :] diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py index 7f9acb7f..72233b09 100644 --- a/LoopStructural/modelling/core/geological_model.py +++ b/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 ) ) diff --git a/LoopStructural/modelling/features/fault/_fault_function.py b/LoopStructural/modelling/features/fault/_fault_function.py index 890ff9d0..b5d58043 100644 --- a/LoopStructural/modelling/features/fault/_fault_function.py +++ b/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) diff --git a/LoopStructural/modelling/features/fault/_fault_segment.py b/LoopStructural/modelling/features/fault/_fault_segment.py index c36e75bf..737b8f8a 100644 --- a/LoopStructural/modelling/features/fault/_fault_segment.py +++ b/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, diff --git a/LoopStructural/modelling/input/process_data.py b/LoopStructural/modelling/input/process_data.py index 2c2a030e..bc985b23 100644 --- a/LoopStructural/modelling/input/process_data.py +++ b/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. diff --git a/tests/unit/datatypes/test__surface.py b/tests/unit/datatypes/test__surface.py index 99a8add0..db0396b4 100644 --- a/tests/unit/datatypes/test__surface.py +++ b/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)