diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py index c915722d8..c781626ca 100644 --- a/LoopStructural/modelling/core/geological_model.py +++ b/LoopStructural/modelling/core/geological_model.py @@ -1288,12 +1288,16 @@ def create_and_add_fault( interpolatortype="FDI", tol=None, fault_slip_vector=None, + fault_normal_vector=None, fault_center=None, major_axis=None, minor_axis=None, intermediate_axis=None, faultfunction="BaseFault", faults=[], + force_mesh_geometry: bool = False, + points: bool = False, + fault_buffer=0.2, **kwargs, ): """ @@ -1357,10 +1361,6 @@ def create_and_add_fault( kwargs.pop("data_region") logger.error("kwarg data_region currently not supported, disabling") displacement_scaled = displacement / self.scale_factor - # create fault frame - # interpolator = self.get_interpolator(**kwargs) - # faults arent supported for surfe - fault_frame_builder = FaultBuilder( interpolatortype, bounding_box=self.bounding_box, @@ -1369,91 +1369,17 @@ def create_and_add_fault( model=self, **kwargs, ) + fault_frame_data = self.data.loc[ + self.data["feature_name"] == fault_surface_data + ].copy() self._add_faults(fault_frame_builder, features=faults) # add data fault_frame_data = self.data.loc[ self.data["feature_name"] == fault_surface_data ].copy() - trace_mask = np.logical_and( - fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0 - ) - logger.info(f"There are {np.sum(trace_mask)} points on the fault trace") - if np.sum(trace_mask) == 0: - logger.error( - "You cannot model a fault without defining the location of the fault" - ) - raise ValueError(f"There are no points on the fault trace") - - mask = np.logical_and( - fault_frame_data["coord"] == 0, ~np.isnan(fault_frame_data["gz"]) - ) - vector_data = fault_frame_data.loc[mask, ["gx", "gy", "gz"]].to_numpy() - mask2 = np.logical_and( - fault_frame_data["coord"] == 0, ~np.isnan(fault_frame_data["nz"]) - ) - vector_data = np.vstack( - [vector_data, fault_frame_data.loc[mask2, ["nx", "ny", "nz"]].to_numpy()] - ) - fault_normal_vector = np.mean(vector_data, axis=0) - logger.info(f"Fault normal vector: {fault_normal_vector}") - - mask = np.logical_and( - fault_frame_data["coord"] == 1, ~np.isnan(fault_frame_data["gz"]) - ) - if fault_slip_vector is None: - if ( - "avgSlipDirEasting" in kwargs - and "avgSlipDirNorthing" in kwargs - and "avgSlipDirAltitude" in kwargs - ): - fault_slip_vector = np.array( - [ - kwargs["avgSlipDirEasting"], - kwargs["avgSlipDirNorthing"], - kwargs["avgSlipDirAltitude"], - ], - dtype=float, - ) - else: - fault_slip_vector = ( - fault_frame_data.loc[mask, ["gx", "gy", "gz"]] - .mean(axis=0) - .to_numpy() - ) - if np.any(np.isnan(fault_slip_vector)): - logger.info("Fault slip vector is nan, estimating from fault normal") - strike_vector, dip_vector = get_vectors(fault_normal_vector[None, :]) - fault_slip_vector = dip_vector[:, 0] - logger.info(f"Estimated fault slip vector: {fault_slip_vector}") if fault_center is not None and ~np.isnan(fault_center).any(): fault_center = self.scale(fault_center, inplace=False) - else: - # if we haven't defined a fault centre take the - # center of mass for lines assocaited with the fault trace - if ( - ~np.isnan(kwargs.get("centreEasting", np.nan)) - and ~np.isnan(kwargs.get("centreNorthing", np.nan)) - and ~np.isnan(kwargs.get("centreAltitude", np.nan)) - ): - fault_center = self.scale( - np.array( - [ - kwargs["centreEasting"], - kwargs["centreNorthing"], - kwargs["centreAltitude"], - ], - dtype=float, - ), - inplace=False, - ) - else: - mask = np.logical_and( - fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0 - ) - fault_center = ( - fault_frame_data.loc[mask, ["X", "Y", "Z"]].mean(axis=0).to_numpy() - ) if minor_axis: minor_axis = minor_axis / self.scale_factor if major_axis: @@ -1461,17 +1387,18 @@ def create_and_add_fault( if intermediate_axis: intermediate_axis = intermediate_axis / self.scale_factor fault_frame_builder.create_data_from_geometry( - fault_frame_data, - fault_center, - fault_normal_vector, - fault_slip_vector, + fault_frame_data=fault_frame_data, + fault_center=fault_center, + fault_normal_vector=fault_normal_vector, + fault_slip_vector=fault_slip_vector, minor_axis=minor_axis, major_axis=major_axis, intermediate_axis=intermediate_axis, - points=kwargs.get("points", False), + points=points, + force_mesh_geometry=force_mesh_geometry, + fault_buffer=fault_buffer, ) - if "force_mesh_geometry" not in kwargs: - fault_frame_builder.set_mesh_geometry(kwargs.get("fault_buffer", 0.2), 0) + if "splay" in kwargs and "splayregion" in kwargs: fault_frame_builder.add_splay(kwargs["splay"], kwargs["splayregion"]) diff --git a/LoopStructural/modelling/features/builders/_fault_builder.py b/LoopStructural/modelling/features/builders/_fault_builder.py index 119041d74..2715bf8a3 100644 --- a/LoopStructural/modelling/features/builders/_fault_builder.py +++ b/LoopStructural/modelling/features/builders/_fault_builder.py @@ -1,7 +1,8 @@ from typing import Union from ._structural_frame_builder import StructuralFrameBuilder - +from LoopStructural.utils import get_vectors import numpy as np +import pandas as pd from ....utils import getLogger, BoundingBox logger = getLogger(__name__) @@ -78,15 +79,17 @@ def update_geometry(self, points): def create_data_from_geometry( self, - data, - fault_center, - normal_vector, - slip_vector, + fault_frame_data: pd.DataFrame, + fault_center=None, + fault_normal_vector=None, + fault_slip_vector=None, minor_axis=None, major_axis=None, intermediate_axis=None, w=1.0, points=False, + force_mesh_geometry=False, + fault_buffer=0.2, ): """Generate the required data for building a fault frame for a fault with the specified parameters @@ -110,13 +113,82 @@ def create_data_from_geometry( intermediate_axis : double fault volume radius in the slip direction """ - self.fault_normal_vector = normal_vector - self.fault_slip_vector = slip_vector + trace_mask = np.logical_and( + fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0 + ) + logger.info(f"There are {np.sum(trace_mask)} points on the fault trace") + if np.sum(trace_mask) == 0: + logger.error( + "You cannot model a fault without defining the location of the fault" + ) + raise ValueError(f"There are no points on the fault trace") + + # get all of the gradient data associated with the fault trace + if fault_normal_vector is None: + gradient_mask = np.logical_and( + fault_frame_data["coord"] == 0, ~np.isnan(fault_frame_data["gz"]) + ) + vector_data = fault_frame_data.loc[ + gradient_mask, ["gx", "gy", "gz"] + ].to_numpy() + normal_mask = np.logical_and( + fault_frame_data["coord"] == 0, ~np.isnan(fault_frame_data["nz"]) + ) + vector_data = np.vstack( + [ + vector_data, + fault_frame_data.loc[normal_mask, ["nx", "ny", "nz"]].to_numpy(), + ] + ) + if len(vector_data) == 0: + logger.error( + "You cannot model a fault without defining the orientation of the fault\n\ + Either add orientation data or define the fault normal vector argument" + ) + raise ValueError(f"There are no points on the fault trace") + fault_normal_vector = np.mean(vector_data, axis=0) + + logger.info(f"Fault normal vector: {fault_normal_vector}") + + # estimate the fault slip vector + + if fault_slip_vector is None: + slip_mask = np.logical_and( + fault_frame_data["coord"] == 1, ~np.isnan(fault_frame_data["gz"]) + ) + fault_slip_data = fault_frame_data.loc[slip_mask, ["gx", "gy", "gz"]] + if len(fault_slip_data) == 0: + logger.warning( + "There is no slip vector data for the fault, using vertical slip vector\n\ + projected onto fault surface estimating from fault normal" + ) + strike_vector, dip_vector = get_vectors(fault_normal_vector[None, :]) + fault_slip_vector = dip_vector[:, 0] + logger.info(f"Estimated fault slip vector: {fault_slip_vector}") + + fault_slip_vector = fault_slip_data.mean(axis=0).to_numpy() + if fault_center is None: + trace_mask = np.logical_and( + fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0 + ) + fault_center = ( + fault_frame_data.loc[trace_mask, ["X", "Y", "Z"]] + .mean(axis=0) + .to_numpy() + ) + if np.any(np.isnan(fault_slip_vector)): + logger.info("Fault slip vector is nan, estimating from fault normal") + + self.fault_normal_vector = fault_normal_vector + self.fault_slip_vector = fault_slip_vector self.fault_centre = fault_center if major_axis is None: - fault_trace = data.loc[ - np.logical_and(data["coord"] == 0, data["val"] == 0), ["X", "Y"] + fault_trace = fault_frame_data.loc[ + np.logical_and( + fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0 + ), + ["X", "Y"], ].to_numpy() distance = np.linalg.norm( fault_trace[:, None, :] - fault_trace[None, :, :], axis=2 @@ -127,9 +199,9 @@ def create_data_from_geometry( # 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(data) - self.origin = self.model.bounding_box[0, :] - self.maximum = self.model.bounding_box[1, :] + 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.warning(f"Fault major axis using map length: {major_axis}") @@ -147,18 +219,18 @@ def create_data_from_geometry( self.fault_minor_axis = minor_axis self.fault_major_axis = major_axis self.fault_intermediate_axis = intermediate_axis - normal_vector /= np.linalg.norm(normal_vector) - slip_vector /= np.linalg.norm(slip_vector) + fault_normal_vector /= np.linalg.norm(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(normal_vector, slip_vector) + strike_vector = np.cross(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)) - data.reset_index(inplace=True) + fault_frame_data.reset_index(inplace=True) if not self.fault_major_axis: logger.warning( f"Fault major axis is not set and cannot be determined from the fault trace. \ @@ -186,8 +258,9 @@ def create_data_from_geometry( # choose whether to add points -1,1 to constrain fault frame or a scaled # vector if points: - data.loc[ - len(data), ["X", "Y", "Z", "feature_name", "val", "coord", "w"] + fault_frame_data.loc[ + len(fault_frame_data), + ["X", "Y", "Z", "feature_name", "val", "coord", "w"], ] = [ fault_edges[0, 0], fault_edges[0, 1], @@ -197,8 +270,9 @@ def create_data_from_geometry( 0, w, ] - data.loc[ - len(data), ["X", "Y", "Z", "feature_name", "val", "coord", "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], @@ -209,32 +283,44 @@ def create_data_from_geometry( w, ] logger.info("Converting fault norm data to gradient data") - mask = np.logical_and(data["coord"] == 0, ~np.isnan(data["nx"])) - data.loc[mask, ["gx", "gy", "gz"]] = data.loc[ - mask, ["nx", "ny", "nz"] - ] + 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"]] - data.loc[mask, ["nx", "ny", "nz"]] = np.nan - mask = np.logical_and(data["coord"] == 0, ~np.isnan(data["gx"])) - data.loc[mask, ["gx", "gy", "gz"]] /= minor_axis * 0.5 + 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 points == False: logger.info( "Rescaling fault norm constraint length for fault frame" ) - mask = np.logical_and(data["coord"] == 0, ~np.isnan(data["gx"])) - data.loc[mask, ["gx", "gy", "gz"]] /= np.linalg.norm( - data.loc[mask, ["gx", "gy", "gz"]], axis=1 + 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 - data.loc[mask, ["gx", "gy", "gz"]] /= minor_axis * 0.5 - mask = np.logical_and(data["coord"] == 0, ~np.isnan(data["nx"])) - data.loc[mask, ["nx", "ny", "nz"]] /= np.linalg.norm( - data.loc[mask, ["nx", "ny", "nz"]], axis=1 + 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 - data.loc[mask, ["nx", "ny", "nz"]] /= minor_axis * 0.5 + 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): @@ -243,8 +329,9 @@ def create_data_from_geometry( 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 - data.loc[ - len(data), ["X", "Y", "Z", "feature_name", "val", "coord", "w"] + fault_frame_data.loc[ + len(fault_frame_data), + ["X", "Y", "Z", "feature_name", "val", "coord", "w"], ] = [ fault_center[0], fault_center[1], @@ -254,8 +341,9 @@ def create_data_from_geometry( 2, w, ] - data.loc[ - len(data), ["X", "Y", "Z", "feature_name", "val", "coord", "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], @@ -265,8 +353,9 @@ def create_data_from_geometry( 2, w, ] - data.loc[ - len(data), ["X", "Y", "Z", "feature_name", "val", "coord", "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], @@ -280,8 +369,9 @@ def create_data_from_geometry( if intermediate_axis is not None: fault_depth[0, :] = fault_center[:3] + slip_vector * intermediate_axis fault_depth[1, :] = fault_center[:3] - slip_vector * intermediate_axis - data.loc[ - len(data), ["X", "Y", "Z", "feature_name", "val", "coord", "w"] + fault_frame_data.loc[ + len(fault_frame_data), + ["X", "Y", "Z", "feature_name", "val", "coord", "w"], ] = [ fault_center[0], fault_center[1], @@ -295,8 +385,8 @@ def create_data_from_geometry( self.update_geometry(fault_depth) # TODO need to add data here slip_vector /= intermediate_axis - data.loc[ - len(data), + fault_frame_data.loc[ + len(fault_frame_data), [ "X", "Y", @@ -321,8 +411,10 @@ def create_data_from_geometry( 1, w, ] - self.add_data_from_data_frame(data) - self.update_geometry(data[["X", "Y", "Z"]].to_numpy()) + self.add_data_from_data_frame(fault_frame_data) + if force_mesh_geometry: + self.set_mesh_geometry(fault_buffer, None) + self.update_geometry(fault_frame_data[["X", "Y", "Z"]].to_numpy()) def set_mesh_geometry(self, buffer, rotation): """set the mesh geometry diff --git a/tests/unit/modelling/test_structural_frame.py b/tests/unit/modelling/test_structural_frame.py index 28bd77f98..73ce63901 100644 --- a/tests/unit/modelling/test_structural_frame.py +++ b/tests/unit/modelling/test_structural_frame.py @@ -3,9 +3,11 @@ GeologicalFeature, ) from LoopStructural.modelling.features.builders import StructuralFrameBuilder +from LoopStructural.datatypes import BoundingBox from LoopStructural import GeologicalModel import numpy as np +import pandas as pd def test_structural_frame(): @@ -48,3 +50,38 @@ def get_item(): def test_structural_frame_cross_product(): pass + + +def test_create_structural_frame_pli(): + data = pd.DataFrame( + [ + [5.1, 5.1, 5, 0, 0, 1, 0, 0], + [5, 5.1, 5, 0, 1, 0, 1, 0], + [5.1, 5, 5, 1, 0, 0, 2, 0], + ], + columns=["X", "Y", "Z", "nx", "ny", "nz", "coord", "val"], + ) + data["feature_name"] = "fault" + + bb = BoundingBox(3, origin=np.zeros(3), maximum=np.ones(3) * 10) + model = GeologicalModel(bb.origin, bb.maximum) + + model.data = data + + fault = model.create_and_add_fault( + "fault", 10, nelements=2000, steps=4, interpolatortype="PLI", buffer=2 + ) + model.update() + assert np.all( + np.isclose(fault[0].evaluate_gradient(np.array([[5, 5, 5]])), [0, 0, 1]) + ) + assert np.all( + np.isclose(fault[1].evaluate_gradient(np.array([[5, 5, 5]])), [0, 1, 0]) + ) + assert np.all( + np.isclose(fault[2].evaluate_gradient(np.array([[5, 5, 5]])), [1, 0, 0]) + ) + + +if __name__ == "__main__": + test_create_structural_frame_pli()