diff --git a/LoopStructural/__init__.py b/LoopStructural/__init__.py index 71b090963..91b4f2de0 100644 --- a/LoopStructural/__init__.py +++ b/LoopStructural/__init__.py @@ -21,7 +21,7 @@ from .modelling.core.geological_model import GeologicalModel from .interpolators._api import LoopInterpolator from .datatypes import BoundingBox -from .utils import log_to_console, log_to_file, getLogger +from .utils import log_to_console, log_to_file, getLogger, rng logger = getLogger(__name__) logger.info("Imported LoopStructural") diff --git a/LoopStructural/modelling/features/fault/_fault_function.py b/LoopStructural/modelling/features/fault/_fault_function.py index a5650649e..7774700dd 100644 --- a/LoopStructural/modelling/features/fault/_fault_function.py +++ b/LoopStructural/modelling/features/fault/_fault_function.py @@ -124,13 +124,27 @@ def solve(self): def __call__(self, v): self.solve() eva = self.w[0] * v**3 + self.w[1] * v**2 + self.w[2] * v + self.w[3] + eva[v > self.lim[1]] = ( + self.w[0] * self.lim[1] ** 3 + + self.w[1] * self.lim[1] ** 2 + + self.w[2] * self.lim[1] + + self.w[3] + ) + eva[v < self.lim[0]] = ( + self.w[0] * self.lim[0] ** 3 + + self.w[1] * self.lim[0] ** 2 + + self.w[2] * self.lim[0] + + self.w[3] + ) eva[eva > self.max_v] = self.max_v eva[eva < self.min_v] = self.min_v + return eva def to_dict(self) -> dict: """Export the function to a dictionary + Returns ------- dict diff --git a/examples/3_fault/plot_define_fault_displacement.py b/examples/3_fault/plot_define_fault_displacement.py new file mode 100644 index 000000000..a331e37ee --- /dev/null +++ b/examples/3_fault/plot_define_fault_displacement.py @@ -0,0 +1,129 @@ +""" +3c. Defining the fault displacement function +============================================ + +""" + +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", + 1, + 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"]] +) + + +import LoopStructural.visualisation as vis + +view = vis.LavaVuModelViewer(model) +view.add_isosurface(model.features[0], slices=[0]) +view.add_isosurface(model.features[1], nslices=5, paint_with=model.features[1]) +# view.add_vector_field(model["fault"][1], locations=model.regular_grid()[::100]) +view.camera = { + 'translate': [0.0, 0.0, -17.321], + 'rotate': [-0.703, -0.055, -0.043, 0.708], + 'xyzrotate': [-89.604, -8.007, 0.933], + 'fov': 45.0, +} +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.LavaVuModelViewer(model) +view.nelements = 1e5 +view.add_isosurface(model.features[0], slices=[0]) +view.add_isosurface(model['strati'], nslices=5) +# view.add_scalar_field(model['strati'], cmap='tab20') +# view.add_vector_field(model["fault"][1], locations=model.regular_grid()[::100]) +view.camera = { + 'translate': [0.0, 0.0, -17.321], + 'rotate': [-0.703, -0.055, -0.043, 0.708], + 'xyzrotate': [-89.604, -8.007, 0.933], + 'fov': 45.0, +} + +view.display()