From 8ab6bca43386b7a7f337bb12d734d6f19db66f96 Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Thu, 8 Jan 2026 11:27:00 +1100 Subject: [PATCH 1/4] fix: add matrix scaling to improve performance --- .../interpolators/_discrete_interpolator.py | 43 +++++++++++++++++-- 1 file changed, 39 insertions(+), 4 deletions(-) diff --git a/LoopStructural/interpolators/_discrete_interpolator.py b/LoopStructural/interpolators/_discrete_interpolator.py index 594ebec7..883e8a43 100644 --- a/LoopStructural/interpolators/_discrete_interpolator.py +++ b/LoopStructural/interpolators/_discrete_interpolator.py @@ -62,7 +62,9 @@ def __init__(self, support, data={}, c=None, up_to_date=False): self.interpolation_weights = {} logger.info("Creating discrete interpolator with {} degrees of freedom".format(self.dof)) self.type = InterpolatorType.BASE_DISCRETE - + self.apply_scaling_matrix = True + self.add_ridge_regulatisation = True + self.ridge_factor = 1e-8 def set_nelements(self, nelements: int) -> int: return self.support.set_nelements(nelements) @@ -511,6 +513,25 @@ def build_matrix(self): B = np.hstack(bs) return A, B + def compute_column_scaling_matrix(self, A: sparse.csr_matrix) -> sparse.dia_matrix: + """Compute column scaling matrix S for matrix A so that A @ S has columns with unit norm. + + Parameters + ---------- + A : sparse.csr_matrix + interpolation matrix + + Returns + ------- + scipy.sparse.dia_matrix + diagonal scaling matrix S + """ + col_norms = sparse.linalg.norm(A, axis=0) + scaling_factors = np.ones(A.shape[1]) + mask = col_norms > 0 + scaling_factors[mask] = 1.0 / col_norms[mask] + S = sparse.diags(scaling_factors) + return S def add_equality_block(self, A, B): if len(self.equal_constraints) > 0: @@ -591,6 +612,15 @@ def solve_system( raise ValueError("Pre solve failed") A, b = self.build_matrix() + if self.add_ridge_regulatisation: + ridge = sparse.eye(A.shape[1]) * self.ridge_factor + A = sparse.vstack([A, ridge]) + b = np.hstack([b, np.zeros(A.shape[1])]) + logger.info("Adding ridge regularisation to interpolation matrix") + if self.apply_scaling_matrix: + S = self.compute_column_scaling_matrix(A) + A = A @ S + Q, bounds = self.build_inequality_matrix() if callable(solver): logger.warning('Using custom solver') @@ -620,12 +650,14 @@ def solve_system( elif solver == 'lsmr': logger.info("Solving using lsmr") - if 'atol' not in solver_kwargs: - if tol is not None: - solver_kwargs['atol'] = tol + # if 'atol' not in solver_kwargs: + # if tol is not None: + # solver_kwargs['atol'] = tol if 'btol' not in solver_kwargs: if tol is not None: solver_kwargs['btol'] = tol + solver_kwargs['atol'] = 0. + logger.info(f"Setting lsmr btol to {tol}") logger.info(f"Solver kwargs: {solver_kwargs}") res = sparse.linalg.lsmr(A, b, **solver_kwargs) if res[1] == 1 or res[1] == 4 or res[1] == 2 or res[1] == 5: @@ -684,6 +716,9 @@ def solve_system( logger.error(f"Unknown solver {solver}") self.up_to_date = False # self._post_solve() + # apply scaling matrix to solution + if self.apply_scaling_matrix: + self.c = S @ self.c return self.up_to_date def update(self) -> bool: From c049fa77bb172152428473cdb54cc82c7ad5b445 Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Thu, 8 Jan 2026 11:28:02 +1100 Subject: [PATCH 2/4] fix: add framework to rebuild fault data if a geometry parameter changes e.g. fault_normal or ellipsoid --- .../features/builders/_fault_builder.py | 752 +++++++++++------- .../3_fault/plot_update_fault_geometry.py | 122 +++ 2 files changed, 571 insertions(+), 303 deletions(-) create mode 100644 examples/3_fault/plot_update_fault_geometry.py diff --git a/LoopStructural/modelling/features/builders/_fault_builder.py b/LoopStructural/modelling/features/builders/_fault_builder.py index 0ee391aa..265db8b0 100644 --- a/LoopStructural/modelling/features/builders/_fault_builder.py +++ b/LoopStructural/modelling/features/builders/_fault_builder.py @@ -51,7 +51,7 @@ def __init__( self.frame.model = model self.origin = np.array([np.nan, np.nan, np.nan]) self.maximum = np.array([np.nan, np.nan, np.nan]) # self.model.bounding_box[1, :] - + self.raw_data = pd.DataFrame() if bounding_box is None: raise ValueError("BoundingBox cannot be None") @@ -59,13 +59,17 @@ def __init__( self.minimum_origin = bounding_box.with_buffer(fault_bounding_box_buffer).origin self.maximum_maximum = bounding_box.with_buffer(fault_bounding_box_buffer).maximum - self.fault_normal_vector = None - self.fault_slip_vector = None - self.fault_strike_vector = None - self.fault_minor_axis = None - self.fault_major_axis = None - self.fault_intermediate_axis = None - self.fault_centre = None + # Use private attributes and property setters to auto-rebuild when geometrical + # attributes change. Suspend rebuild during initialization. + self._suspend_rebuild = True + self._fault_normal_vector = None + self._fault_slip_vector = None + self._fault_strike_vector = None + self._fault_minor_axis = None + self._fault_major_axis = None + self._fault_intermediate_axis = None + self._fault_centre = None + self._suspend_rebuild = False def update_geometry(self, points): """ @@ -92,6 +96,139 @@ def update_geometry(self, points): self.maximum > self.maximum_maximum ] + def _on_geometry_change( + self, + force_mesh_geometry: bool = False, + fault_buffer: float | None = None, + fault_trace_anisotropy: float = 1.0, + fault_dip_anisotropy: float = 1.0, + fault_dip: float | None = None, + fault_pitch: float | None = None, + ): + """Rebuild data from the current stored geometrical attributes. + + This calls create_data_from_geometry with a copy of the last raw_data so that + changing properties triggers a rebuild. + """ + if getattr(self, "raw_data", None) is None or self.raw_data.empty: + return + + # Use a copy of the DataFrame so callers don't observe mutations + self.create_data_from_geometry( + self.raw_data.copy(), + fault_center=self._fault_centre, + fault_normal_vector=self._fault_normal_vector, + fault_slip_vector=self._fault_slip_vector, + minor_axis=self._fault_minor_axis, + major_axis=self._fault_major_axis, + intermediate_axis=self._fault_intermediate_axis, + w=1.0, + points=False, + force_mesh_geometry=force_mesh_geometry, + fault_buffer=fault_buffer if fault_buffer is not None else 0.2, + fault_trace_anisotropy=fault_trace_anisotropy, + fault_dip=fault_dip if fault_dip is not None else getattr(self, "fault_dip", 90), + fault_dip_anisotropy=fault_dip_anisotropy, + fault_pitch=fault_pitch, + ) + + @property + def fault_normal_vector(self): + return self._fault_normal_vector + + @fault_normal_vector.setter + def fault_normal_vector(self, v): + if v is None: + self._fault_normal_vector = None + else: + if len(v.shape) != 1 or v.shape[0] != 3: + raise ValueError("fault_normal_vector must be a 3 element array") + if v.ndim != 1: + v = v[0,:] + if v.shape[0] != 3: + raise ValueError("fault_normal_vector must be a 3 element array") + arr = np.array(v, dtype=float) + norm = np.linalg.norm(arr) + if norm == 0: + raise ValueError("fault_normal_vector cannot be the zero vector") + self._fault_normal_vector = arr + # trigger rebuild unless suspended + if not getattr(self, "_suspend_rebuild", False): + self._on_geometry_change() + + @property + def fault_slip_vector(self): + return self._fault_slip_vector + + @fault_slip_vector.setter + def fault_slip_vector(self, v): + if v is None: + self._fault_slip_vector = None + else: + arr = np.array(v, dtype=float) + norm = np.linalg.norm(arr) + if norm == 0: + raise ValueError("fault_slip_vector cannot be the zero vector") + self._fault_slip_vector = arr / norm + if not getattr(self, "_suspend_rebuild", False): + self._on_geometry_change() + + @property + def fault_strike_vector(self): + return self._fault_strike_vector + + @fault_strike_vector.setter + def fault_strike_vector(self, v): + self._fault_strike_vector = None if v is None else np.array(v, dtype=float) + if not getattr(self, "_suspend_rebuild", False): + self._on_geometry_change() + + @property + def fault_minor_axis(self): + return self._fault_minor_axis + + @fault_minor_axis.setter + def fault_minor_axis(self, v): + self._fault_minor_axis = None if v is None else float(v) + if not getattr(self, "_suspend_rebuild", False): + self._on_geometry_change() + + @property + def fault_major_axis(self): + return self._fault_major_axis + + @fault_major_axis.setter + def fault_major_axis(self, v): + self._fault_major_axis = None if v is None else float(v) + if not getattr(self, "_suspend_rebuild", False): + self._on_geometry_change() + + @property + def fault_intermediate_axis(self): + return self._fault_intermediate_axis + + @fault_intermediate_axis.setter + def fault_intermediate_axis(self, v): + self._fault_intermediate_axis = None if v is None else float(v) + if not getattr(self, "_suspend_rebuild", False): + self._on_geometry_change() + + @property + def fault_centre(self): + return self._fault_centre + + @fault_centre.setter + def fault_centre(self, v): + if v is None: + self._fault_centre = None + else: + arr = np.array(v, dtype=float) + if arr.shape[0] != 3: + raise ValueError("fault_center must be a 3 element array") + self._fault_centre = arr + if not getattr(self, "_suspend_rebuild", False): + self._on_geometry_change() + def create_data_from_geometry( self, fault_frame_data: pd.DataFrame, @@ -146,325 +283,334 @@ def create_data_from_geometry( fault_pitch : float, optional Pitch angle of the fault. """ - fault_trace = fault_frame_data.loc[ - np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0), - ["X", "Y"], - ].to_numpy() - self.fault_dip = fault_dip - if fault_normal_vector is None: - if fault_frame_data.loc[ - np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["nx"].notna())].shape[0]>0: - fault_normal_vector = fault_frame_data.loc[ - np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["nx"].notna()), - ["nx", "ny", "nz"], - ].to_numpy().mean(axis=0) + # Suspend automatic rebuilds while we populate internal attributes to avoid recursion + prev_suspend = getattr(self, "_suspend_rebuild", False) + self._suspend_rebuild = True + try: + self.raw_data = fault_frame_data.copy() + fault_trace = fault_frame_data.loc[ + np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0), + ["X", "Y"], + ].to_numpy() + self.fault_dip = fault_dip + if fault_normal_vector is None: + if fault_frame_data.loc[ + np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["nx"].notna())].shape[0]>0: + fault_normal_vector = fault_frame_data.loc[ + np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["nx"].notna()), + ["nx", "ny", "nz"], + ].to_numpy().mean(axis=0) - else: + else: - # Calculate fault strike using eigenvectors - pts = fault_trace - fault_trace.mean(axis=0) - # Calculate covariance matrix - cov_matrix = pts.T @ pts - # Get eigenvectors and eigenvalues - eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix) - # Use eigenvector with largest eigenvalue as strike direction - strike_vector = eigenvectors[:, np.argmax(eigenvalues)] - strike_vector = np.append(strike_vector, 0) # Add z component - strike_vector /= np.linalg.norm(strike_vector) - - fault_normal_vector = np.cross(strike_vector, [0, 0, 1]) - # Rotate the fault normal vector according to the fault dip - rotation_matrix = rotation(strike_vector[None, :], np.array([90 - fault_dip])) - fault_normal_vector = np.einsum("ijk,ik->ij", rotation_matrix, fault_normal_vector[None, :])[0] - - if not isinstance(fault_normal_vector, np.ndarray): - fault_normal_vector = np.array(fault_normal_vector) - - if fault_pitch is not None: - rotation_matrix = rotation(fault_normal_vector[None, :], np.array([fault_pitch])) - fault_slip_vector = np.einsum("ijk,ik->ij", rotation_matrix, fault_normal_vector[None, :])[0] - - if fault_slip_vector is None: - if fault_frame_data.loc[ - np.logical_and(fault_frame_data["coord"] == 1, fault_frame_data["nx"].notna())].shape[0]>0: - fault_slip_vector = fault_frame_data.loc[ - np.logical_and(fault_frame_data["coord"] == 1, fault_frame_data["nx"].notna()), - ["nx", "ny", "nz"], - ].to_numpy().mean(axis=0) + # Calculate fault strike using eigenvectors + pts = fault_trace - fault_trace.mean(axis=0) + # Calculate covariance matrix + cov_matrix = pts.T @ pts + # Get eigenvectors and eigenvalues + eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix) + # Use eigenvector with largest eigenvalue as strike direction + strike_vector = eigenvectors[:, np.argmax(eigenvalues)] + strike_vector = np.append(strike_vector, 0) # Add z component + strike_vector /= np.linalg.norm(strike_vector) + + fault_normal_vector = np.cross(strike_vector, [0, 0, 1]) + # Rotate the fault normal vector according to the fault dip + rotation_matrix = rotation(strike_vector[None, :], np.array([90 - fault_dip])) + fault_normal_vector = np.einsum("ijk,ik->ij", rotation_matrix, fault_normal_vector[None, :])[0] + + if not isinstance(fault_normal_vector, np.ndarray): + fault_normal_vector = np.array(fault_normal_vector) + + if fault_pitch is not None: + rotation_matrix = rotation(fault_normal_vector[None, :], np.array([fault_pitch])) + fault_slip_vector = np.einsum("ijk,ik->ij", rotation_matrix, fault_normal_vector[None, :])[0] + + if fault_slip_vector is None: + if fault_frame_data.loc[ + np.logical_and(fault_frame_data["coord"] == 1, fault_frame_data["nx"].notna())].shape[0]>0: + fault_slip_vector = fault_frame_data.loc[ + np.logical_and(fault_frame_data["coord"] == 1, fault_frame_data["nx"].notna()), + ["nx", "ny", "nz"], + ].to_numpy().mean(axis=0) - else: - fault_slip_vector = np.cross(fault_normal_vector, [1., 0., 0.]) - if np.linalg.norm(fault_slip_vector) == 0: - fault_slip_vector = np.cross(fault_normal_vector, [0., 1., 0.]) - fault_slip_vector /= np.linalg.norm(fault_slip_vector) - if fault_center is None: - fault_trace = fault_frame_data.loc[ + else: + fault_slip_vector = np.cross(fault_normal_vector, [1., 0., 0.]) + if np.linalg.norm(fault_slip_vector) == 0: + fault_slip_vector = np.cross(fault_normal_vector, [0., 1., 0.]) + fault_slip_vector /= np.linalg.norm(fault_slip_vector) + if fault_center is None: + fault_trace = fault_frame_data.loc[ np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0), ["X", "Y"], ].to_numpy() - fault_center = fault_trace.mean(axis=0) - fault_center = np.array([fault_center[0], fault_center[1], 0.0]) - if not isinstance(fault_center, np.ndarray): - fault_center = np.array(fault_center) - if fault_center.shape[0] != 3: - raise ValueError("fault_center must be a 3 element array") - self.fault_normal_vector = fault_normal_vector / np.linalg.norm(fault_normal_vector) - self.fault_slip_vector = fault_slip_vector / np.linalg.norm(fault_slip_vector) - - self.fault_centre = fault_center - if major_axis is None: - - distance = np.linalg.norm(fault_trace[:, None, :] - fault_trace[None, :, :], axis=2) - if len(distance) == 0 or np.sum(distance) == 0: - logger.warning("There is no fault trace for {}".format(self.name)) - # this can mean there is only a single data point for - # the fault, its not critical - # but probably means the fault isn't well defined. - # add any data anyway - usually just orientation data - self.add_data_from_data_frame(fault_frame_data) - self.origin = self.model.bounding_box.origin - self.maximum = self.model.bounding_box.maximum - return - major_axis = np.max(distance) - logger.info(f"Fault major axis using map length: {major_axis}") - - if minor_axis is None: - logger.info(f"Fault minor axis not set, using half major axis: {major_axis/2}") - minor_axis = major_axis / 2.0 - if intermediate_axis is None: - intermediate_axis = major_axis - logger.info(f"Fault intermediate axis not set, using major axis: {intermediate_axis}") - self.fault_minor_axis = minor_axis - self.fault_major_axis = major_axis - self.fault_intermediate_axis = intermediate_axis - self.fault_normal_vector /= np.linalg.norm(self.fault_normal_vector) - fault_slip_vector /= np.linalg.norm(fault_slip_vector) - # check if slip vector is inside fault plane, if not project onto fault plane - # if not np.isclose(normal_vector @ slip_vector, 0): - strike_vector = np.cross(self.fault_normal_vector, fault_slip_vector) - self.fault_strike_vector = strike_vector - - fault_edges = np.zeros((2, 3)) - fault_tips = np.zeros((2, 3)) - fault_depth = np.zeros((2, 3)) - fault_frame_data.reset_index(inplace=True) - if not self.fault_major_axis: - logger.warning( - "Fault major axis is not set and cannot be determined from the fault trace. \ - This will result in a fault that is represented by a 1 unit major axis. \ - If this is not intended add major_axis to fault parameters." - ) - if not self.fault_intermediate_axis: - logger.warning( - "Fault intermediate axis is not set and cannot be determined from the fault trace. \ - This will result in a fault that is represented by a 1 unit intermediate axis. \ - If this is not intended add intermediate_axis to fault parameters." - ) - if not self.fault_minor_axis: - logger.warning( - "Fault minor axis is not set and cannot be determined from the fault trace. \ - This will result in a fault that is represented by a 1 unit minor axis. \ - If this is not intended add minor_axis to fault parameters." - ) - if fault_center is not None: - if minor_axis is not None: - fault_edges[0, :] = fault_center[:3] + self.fault_normal_vector * minor_axis - fault_edges[1, :] = fault_center[:3] - self.fault_normal_vector * minor_axis - self.update_geometry(fault_edges) - - # choose whether to add points -1,1 to constrain fault frame or a scaled - # vector - if points: + fault_center = fault_trace.mean(axis=0) + fault_center = np.array([fault_center[0], fault_center[1], 0.0]) + if not isinstance(fault_center, np.ndarray): + fault_center = np.array(fault_center) + if fault_center.shape[0] != 3: + raise ValueError("fault_center must be a 3 element array") + # Assign to properties (these won't trigger _on_geometry_change because suspend flag is True) + self.fault_normal_vector = fault_normal_vector / np.linalg.norm(fault_normal_vector) + self.fault_slip_vector = fault_slip_vector / np.linalg.norm(fault_slip_vector) + + self.fault_centre = fault_center + if major_axis is None: + + distance = np.linalg.norm(fault_trace[:, None, :] - fault_trace[None, :, :], axis=2) + if len(distance) == 0 or np.sum(distance) == 0: + logger.warning("There is no fault trace for {}".format(self.name)) + # this can mean there is only a single data point for + # the fault, its not critical + # but probably means the fault isn't well defined. + # add any data anyway - usually just orientation data + self.add_data_from_data_frame(fault_frame_data) + self.origin = self.model.bounding_box.origin + self.maximum = self.model.bounding_box.maximum + return + major_axis = np.max(distance) + logger.info(f"Fault major axis using map length: {major_axis}") + + if minor_axis is None: + logger.info(f"Fault minor axis not set, using half major axis: {major_axis/2}") + minor_axis = major_axis / 2.0 + if intermediate_axis is None: + intermediate_axis = major_axis + logger.info(f"Fault intermediate axis not set, using major axis: {intermediate_axis}") + self.fault_minor_axis = minor_axis + self.fault_major_axis = major_axis + self.fault_intermediate_axis = intermediate_axis + self.fault_normal_vector /= np.linalg.norm(self.fault_normal_vector) + fault_slip_vector /= np.linalg.norm(fault_slip_vector) + # check if slip vector is inside fault plane, if not project onto fault plane + # if not np.isclose(normal_vector @ slip_vector, 0): + strike_vector = np.cross(self.fault_normal_vector, fault_slip_vector) + self.fault_strike_vector = strike_vector + + fault_edges = np.zeros((2, 3)) + fault_tips = np.zeros((2, 3)) + fault_depth = np.zeros((2, 3)) + fault_frame_data.reset_index(inplace=True) + if not self.fault_major_axis: + logger.warning( + "Fault major axis is not set and cannot be determined from the fault trace. \ + This will result in a fault that is represented by a 1 unit major axis. \ + If this is not intended add major_axis to fault parameters." + ) + if not self.fault_intermediate_axis: + logger.warning( + "Fault intermediate axis is not set and cannot be determined from the fault trace. \ + This will result in a fault that is represented by a 1 unit intermediate axis. \ + If this is not intended add intermediate_axis to fault parameters." + ) + if not self.fault_minor_axis: + logger.warning( + "Fault minor axis is not set and cannot be determined from the fault trace. \ + This will result in a fault that is represented by a 1 unit minor axis. \ + If this is not intended add minor_axis to fault parameters." + ) + if fault_center is not None: + if minor_axis is not None: + fault_edges[0, :] = fault_center[:3] + self.fault_normal_vector * minor_axis + fault_edges[1, :] = fault_center[:3] - self.fault_normal_vector * minor_axis + self.update_geometry(fault_edges) + + # choose whether to add points -1,1 to constrain fault frame or a scaled + # vector + if points: + fault_frame_data.loc[ + len(fault_frame_data), + ["X", "Y", "Z", "feature_name", "val", "coord", "w"], + ] = [ + fault_edges[0, 0], + fault_edges[0, 1], + fault_edges[0, 2], + self.name, + 1, + 0, + w, + ] + fault_frame_data.loc[ + len(fault_frame_data), + ["X", "Y", "Z", "feature_name", "val", "coord", "w"], + ] = [ + fault_edges[1, 0], + fault_edges[1, 1], + fault_edges[1, 2], + self.name, + -1, + 0, + w, + ] + logger.info("Converting fault norm data to gradient data") + mask = np.logical_and( + fault_frame_data["coord"] == 0, + ~np.isnan(fault_frame_data["nx"]), + ) + fault_frame_data.loc[mask, ["gx", "gy", "gz"]] = fault_frame_data.loc[ + mask, ["nx", "ny", "nz"] + ] + + fault_frame_data.loc[mask, ["nx", "ny", "nz"]] = np.nan + mask = np.logical_and( + fault_frame_data["coord"] == 0, + ~np.isnan(fault_frame_data["gx"]), + ) + fault_frame_data.loc[mask, ["gx", "gy", "gz"]] /= minor_axis * 0.5 + if not points: + logger.info("Rescaling fault norm constraint length for fault frame") + mask = np.logical_and( + fault_frame_data["coord"] == 0, + ~np.isnan(fault_frame_data["gx"]), + ) + fault_frame_data.loc[mask, ["gx", "gy", "gz"]] /= np.linalg.norm( + fault_frame_data.loc[mask, ["gx", "gy", "gz"]], axis=1 + )[:, None] + # scale vector so that the distance between -1 + # and 1 is the minor axis length + fault_frame_data.loc[mask, ["gx", "gy", "gz"]] /= minor_axis * 0.5 + mask = np.logical_and( + fault_frame_data["coord"] == 0, + ~np.isnan(fault_frame_data["nx"]), + ) + fault_frame_data.loc[mask, ["nx", "ny", "nz"]] /= np.linalg.norm( + fault_frame_data.loc[mask, ["nx", "ny", "nz"]], axis=1 + )[:, None] + # scale vector so that the distance between -1 + # and 1 is the minor axis length + fault_frame_data.loc[mask, ["nx", "ny", "nz"]] /= minor_axis * 0.5 + # self.builders[0].add_orthogonal_feature(self, + # feature, w=1.0, region=None, step=1, B=0): + if np.sum(mask) == 0: + fault_frame_data.loc[ + len(fault_frame_data), + [ + "X", + "Y", + "Z", + "feature_name", + "nx", + "ny", + "nz", + "val", + "coord", + "w", + ], + ] = [ + fault_center[0], + fault_center[1], + fault_center[2], + self.name, + self.fault_normal_vector[0] / minor_axis * 0.5, + self.fault_normal_vector[1] / minor_axis * 0.5, + self.fault_normal_vector[2] / minor_axis * 0.5, + np.nan, + 0, + w, + ] + + if major_axis is not None: + fault_tips[0, :] = fault_center[:3] + strike_vector * 0.5 * major_axis + fault_tips[1, :] = fault_center[:3] - strike_vector * 0.5 * major_axis + self.update_geometry(fault_tips) + # we want the tips of the fault to be -1 and 1 fault_frame_data.loc[ len(fault_frame_data), ["X", "Y", "Z", "feature_name", "val", "coord", "w"], ] = [ - fault_edges[0, 0], - fault_edges[0, 1], - fault_edges[0, 2], + fault_center[0], + fault_center[1], + fault_center[2], self.name, - 1, 0, + 2, w, ] fault_frame_data.loc[ len(fault_frame_data), ["X", "Y", "Z", "feature_name", "val", "coord", "w"], ] = [ - fault_edges[1, 0], - fault_edges[1, 1], - fault_edges[1, 2], + fault_tips[1, 0], + fault_tips[1, 1], + fault_tips[1, 2], self.name, - -1, - 0, + -0.5, + 2, w, ] - logger.info("Converting fault norm data to gradient data") - mask = np.logical_and( - fault_frame_data["coord"] == 0, - ~np.isnan(fault_frame_data["nx"]), - ) - fault_frame_data.loc[mask, ["gx", "gy", "gz"]] = fault_frame_data.loc[ - mask, ["nx", "ny", "nz"] + fault_frame_data.loc[ + len(fault_frame_data), + ["X", "Y", "Z", "feature_name", "val", "coord", "w"], + ] = [ + fault_tips[0, 0], + fault_tips[0, 1], + fault_tips[0, 2], + self.name, + 0.5, + 2, + w, + ] + strike_vector /= major_axis + if intermediate_axis is not None: + fault_depth[0, :] = fault_center[:3] + fault_slip_vector * intermediate_axis + fault_depth[1, :] = fault_center[:3] - fault_slip_vector * intermediate_axis + fault_frame_data.loc[ + len(fault_frame_data), + ["X", "Y", "Z", "feature_name", "val", "coord", "w"], + ] = [ + fault_center[0], + fault_center[1], + fault_center[2], + self.name, + 0, + 1, + w, ] - fault_frame_data.loc[mask, ["nx", "ny", "nz"]] = np.nan - mask = np.logical_and( - fault_frame_data["coord"] == 0, - ~np.isnan(fault_frame_data["gx"]), - ) - fault_frame_data.loc[mask, ["gx", "gy", "gz"]] /= minor_axis * 0.5 - if not points: - logger.info("Rescaling fault norm constraint length for fault frame") - mask = np.logical_and( - fault_frame_data["coord"] == 0, - ~np.isnan(fault_frame_data["gx"]), - ) - fault_frame_data.loc[mask, ["gx", "gy", "gz"]] /= np.linalg.norm( - fault_frame_data.loc[mask, ["gx", "gy", "gz"]], axis=1 - )[:, None] - # scale vector so that the distance between -1 - # and 1 is the minor axis length - fault_frame_data.loc[mask, ["gx", "gy", "gz"]] /= minor_axis * 0.5 - mask = np.logical_and( - fault_frame_data["coord"] == 0, - ~np.isnan(fault_frame_data["nx"]), - ) - fault_frame_data.loc[mask, ["nx", "ny", "nz"]] /= np.linalg.norm( - fault_frame_data.loc[mask, ["nx", "ny", "nz"]], axis=1 - )[:, None] - # scale vector so that the distance between -1 - # and 1 is the minor axis length - fault_frame_data.loc[mask, ["nx", "ny", "nz"]] /= minor_axis * 0.5 - # self.builders[0].add_orthogonal_feature(self, - # feature, w=1.0, region=None, step=1, B=0): - if np.sum(mask) == 0: - fault_frame_data.loc[ - len(fault_frame_data), - [ - "X", - "Y", - "Z", - "feature_name", - "nx", - "ny", - "nz", - "val", - "coord", - "w", - ], - ] = [ - fault_center[0], - fault_center[1], - fault_center[2], - self.name, - self.fault_normal_vector[0] / minor_axis * 0.5, - self.fault_normal_vector[1] / minor_axis * 0.5, - self.fault_normal_vector[2] / minor_axis * 0.5, - np.nan, - 0, - w, - ] + self.update_geometry(fault_depth) + # TODO need to add data here + fault_slip_vector /= intermediate_axis + fault_frame_data.loc[ + len(fault_frame_data), + [ + "X", + "Y", + "Z", + "feature_name", + "nx", + "ny", + "nz", + "val", + "coord", + "w", + ], + ] = [ + fault_center[0], + fault_center[1], + fault_center[2], + self.name, + fault_slip_vector[0], + fault_slip_vector[1], + fault_slip_vector[2], + 0, + 1, + w, + ] - if major_axis is not None: - fault_tips[0, :] = fault_center[:3] + strike_vector * 0.5 * major_axis - fault_tips[1, :] = fault_center[:3] - strike_vector * 0.5 * major_axis - self.update_geometry(fault_tips) - # we want the tips of the fault to be -1 and 1 - fault_frame_data.loc[ - len(fault_frame_data), - ["X", "Y", "Z", "feature_name", "val", "coord", "w"], - ] = [ - fault_center[0], - fault_center[1], - fault_center[2], - self.name, - 0, - 2, - w, - ] - fault_frame_data.loc[ - len(fault_frame_data), - ["X", "Y", "Z", "feature_name", "val", "coord", "w"], - ] = [ - fault_tips[1, 0], - fault_tips[1, 1], - fault_tips[1, 2], - self.name, - -0.5, - 2, - w, - ] - fault_frame_data.loc[ - len(fault_frame_data), - ["X", "Y", "Z", "feature_name", "val", "coord", "w"], - ] = [ - fault_tips[0, 0], - fault_tips[0, 1], - fault_tips[0, 2], - self.name, - 0.5, - 2, - w, - ] - strike_vector /= major_axis - if intermediate_axis is not None: - fault_depth[0, :] = fault_center[:3] + fault_slip_vector * intermediate_axis - fault_depth[1, :] = fault_center[:3] - fault_slip_vector * intermediate_axis - fault_frame_data.loc[ - len(fault_frame_data), - ["X", "Y", "Z", "feature_name", "val", "coord", "w"], - ] = [ - fault_center[0], - fault_center[1], - fault_center[2], - self.name, - 0, - 1, - w, - ] - - self.update_geometry(fault_depth) - # TODO need to add data here - fault_slip_vector /= intermediate_axis - fault_frame_data.loc[ - len(fault_frame_data), - [ - "X", - "Y", - "Z", - "feature_name", - "nx", - "ny", - "nz", - "val", - "coord", - "w", - ], - ] = [ - fault_center[0], - fault_center[1], - fault_center[2], - self.name, - fault_slip_vector[0], - fault_slip_vector[1], - fault_slip_vector[2], - 0, - 1, - w, - ] - - self.add_data_from_data_frame(fault_frame_data) - if fault_trace_anisotropy > 0: - self.add_fault_trace_anisotropy(w=fault_trace_anisotropy) - if fault_dip_anisotropy > 0: - self.add_fault_dip_anisotropy(w=fault_dip_anisotropy) - if force_mesh_geometry: - self.origin = self.model.bounding_box.origin - self.maximum = self.model.bounding_box.maximum - else: - self.update_geometry(fault_frame_data[["X", "Y", "Z"]].to_numpy()) - self.set_mesh_geometry(fault_buffer, None) + self.add_data_from_data_frame(fault_frame_data) + if fault_trace_anisotropy > 0: + self.add_fault_trace_anisotropy(w=fault_trace_anisotropy) + if fault_dip_anisotropy > 0: + self.add_fault_dip_anisotropy(w=fault_dip_anisotropy) + if force_mesh_geometry: + self.origin = self.model.bounding_box.origin + self.maximum = self.model.bounding_box.maximum + else: + self.update_geometry(fault_frame_data[["X", "Y", "Z"]].to_numpy()) + self.set_mesh_geometry(fault_buffer, None) + finally: + # restore previous suspend state so property setters behave as before + self._suspend_rebuild = prev_suspend def set_mesh_geometry(self, buffer, rotation): """set the mesh geometry diff --git a/examples/3_fault/plot_update_fault_geometry.py b/examples/3_fault/plot_update_fault_geometry.py new file mode 100644 index 00000000..21c2dfe9 --- /dev/null +++ b/examples/3_fault/plot_update_fault_geometry.py @@ -0,0 +1,122 @@ +""" +3d. Updating fault geometry +============================================ + +""" + +import numpy as np +import pandas as pd +import LoopStructural as LS + +# Define a dataset for a fault + +origin = [0, 0, 0] +extent = [10, 10, 10] + +data = pd.DataFrame( + [ + [5, 5, 5, 0, 0.70710678, 0.0, 0.70710678, 0, "fault"], + [5, 5, 5, 0, -0.70710678, 0.0, 0.70710678, 1, "fault"], + [8, 5, 5, 0, 0, 0, 1, np.nan, "strati"], + ], + columns=["X", "Y", "Z", "val", "nx", "ny", "nz", "coord", "feature_name"], +) + +data + +###################################################################### +# Create model using the standard fault displacement model + +model = LS.GeologicalModel(origin, extent) +model.data = data +model.create_and_add_fault( + "fault", + 10, + nelements=1000, + interpolator_type="PLI", + buffer=0.5, + major_axis=10, + minor_axis=3, + intermediate_axis=10, +) +model.create_and_add_foliation( + "strati", nelements=1000, interpolator_type="PLI", faults=[model["fault"]] +) +model.update() +print(model['fault'].builder.fault_minor_axis, model['fault'].builder.up_to_date) + +model['fault'].builder.fault_minor_axis = 6.0 +print(model['fault'].builder.fault_minor_axis, model['fault'].builder.up_to_date) +model.update() +print(model['fault'].builder.fault_minor_axis, model['fault'].builder.up_to_date) + +# import LoopStructural.visualisation as vis + +# view = vis.Loop3DView(model) +# view.plot_surface(model.features[0], value=[0]) +# view.plot_surface(model.features[1], value=5, paint_with=model.features[1]) +# # view.add_vector_field(model["fault"][1], locations=model.regular_grid()[::100]) + +# view.display() + +# ###################################################################### +# # Define a fault displacement profile which +# # is a drag fault only on the footwall side. +# # In LoopStructural the displacement is defined by a function of the three +# # coordinates of the fault frame. +# # The fault profile in the fault surface field + +# model['fault'].faultfunction.gx.plot() + +# ###################################################################### +# # The fault profile in the fault extent +# model['fault'].faultfunction.gy.plot() + + +# ###################################################################### +# # The fault profile down dip is kept constant. +# # We will modify this profile so that the hanging wall is displaced by a constant value + +# from LoopStructural.modelling.features.fault._fault_function import ( +# FaultDisplacement, +# CubicFunction, +# Ones, +# ) + +# fw = CubicFunction() +# fw.add_cstr(0, -1) +# fw.add_grad(0, 0) +# fw.add_cstr(-1, 0) +# fw.add_grad(-1, 0) +# fw.add_min(-1) +# hw = Ones() +# drag_fault = FaultDisplacement(hw=hw, fw=fw) + +# drag_fault.gx.plot() +# drag_fault.gy.plot() +# drag_fault.gz.plot() + +# model = LS.GeologicalModel(origin, extent) +# model.data = data +# model.create_and_add_fault( +# "fault", +# -1, +# nelements=1000, +# interpolator_type="PLI", +# buffer=0.5, +# major_axis=10, +# minor_axis=6, +# intermediate_axis=10, +# faultfunction=drag_fault, +# ) +# model.create_and_add_foliation( +# "strati", nelements=1000, interpolator_type="PLI", faults=[model["fault"]] +# ) + + +# view = vis.Loop3DView(model) +# model.bounding_box.nelements = 1e5 +# view.plot_surface(model.features[0], value=[0]) +# view.plot_surface(model['strati'], value=5) + +# view.display() From 48823fe533665338758e229a05a409cf80279944 Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Thu, 8 Jan 2026 11:28:47 +1100 Subject: [PATCH 3/4] fix: strikedip2vector works for single floats not only list like --- LoopStructural/utils/maths.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/LoopStructural/utils/maths.py b/LoopStructural/utils/maths.py index 545e36db..6ec8d305 100644 --- a/LoopStructural/utils/maths.py +++ b/LoopStructural/utils/maths.py @@ -19,6 +19,15 @@ def strikedip2vector(strike: NumericInput, dip: NumericInput) -> np.ndarray: _type_ _description_ """ + if isinstance(strike, numbers.Number): + strike = np.array([strike]) + else: + strike = np.array(strike) + if isinstance(dip, numbers.Number): + dip = np.array([dip]) + else: + dip = np.array(dip) + vec = np.zeros((len(strike), 3)) s_r = np.deg2rad(strike) d_r = np.deg2rad((dip)) From 23e0b2104e5fcada715c53886b2bd006444b4e9e Mon Sep 17 00:00:00 2001 From: Copilot <198982749+Copilot@users.noreply.github.com> Date: Thu, 15 Jan 2026 10:54:42 +1100 Subject: [PATCH 4/4] Initial plan (#278) Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com>