Merge branch 'master' of https://github.com/Loop3D/LoopStructural

Author: lachlangrose

LoopStructural/datasets/__init__.py

@@ -16,3 +16,4 @@
 from ._base import value_headers
 from ._base import load_unconformity
 from ._base import load_duplex
+from ._base import load_tabular_intrusion

LoopStructural/datasets/_base.py

@@ -142,3 +142,18 @@ def strike_dip_headers():
 
 def normal_vector_headers():
     return ["X", "Y", "Z", "nx", "ny", "nz"]
+
+
+def load_tabular_intrusion():
+    """Model dataset sampled for a model of a tabular intrusion
+
+
+    Returns
+    -------
+    tuple
+        pandas data frame with loopstructural dataset and numpy array for bounding box
+    """
+    module_path = dirname(__file__)
+    data = pd.read_csv(join(module_path, Path("data/tabular_intrusion.csv")))
+    bb = np.array([[0, 0, 0], [5, 5, 5]])
+    return data, bb

LoopStructural/datasets/data/tabular_intrusion.csv

@@ -0,0 +1,23 @@
+feature_name,X,Y,Z,coord,val,gx,gy,gz,intrusion_contact_type,intrusion_anisotropy,intrusion_side
+stratigraphy,0.00,0.00,2.00,,0,,,,,,
+stratigraphy,0.00,5.00,2.00,,0,,,,,,
+stratigraphy,5.00,0.00,2.00,,0,,,,,,
+stratigraphy,5.00,5.00,2.00,,0,,,,,,
+stratigraphy,2.50,2.50,2.00,,,0,0,1,,,
+tabular_intrusion,3.04,2.12,2.18,,,,,,roof,stratigraphy,
+tabular_intrusion,4.02,3.85,2.14,,,,,,roof,stratigraphy,TRUE
+tabular_intrusion,2.02,3.16,2.02,,,,,,roof,stratigraphy,
+tabular_intrusion,2.12,2.16,2.02,,,,,,roof,stratigraphy,
+tabular_intrusion,1.14,3.50,2.04,,,,,,roof,,TRUE
+tabular_intrusion,4.08,3.02,1.18,,,,,,floor,,TRUE
+tabular_intrusion,1.14,1.06,1.16,,,,,,floor,,TRUE
+tabular_intrusion,4.20,2.18,1.12,,,,,,floor,,TRUE
+tabular_intrusion,1.02,3.02,1.12,,,,,,floor,,TRUE
+tabular_intrusion_frame,2.00,2.00,2.00,0,,0,0,-1,,,
+tabular_intrusion_frame,3.00,1.00,2.00,0,,0,0,-1,,,
+tabular_intrusion_frame,1.00,3.00,2.00,0,,0,0,-1,,,
+tabular_intrusion_frame,3.00,2.00,1.00,1,0,,,,,,
+tabular_intrusion_frame,3.00,2.00,1.00,1,,0,1,0,,,
+tabular_intrusion_frame,2.50,1.00,1.00,2,0,,,,,,
+tabular_intrusion_frame,2.50,2.00,1.00,2,0,,,,,,
+tabular_intrusion_frame,2.50,2.00,1.00,2,,1,0,0,,,

LoopStructural/interpolators/__init__.py

@@ -50,11 +50,12 @@ class InterpolatorType(IntEnum):
 from LoopStructural.interpolators._discrete_fold_interpolator import (
     DiscreteFoldInterpolator,
 )
-from LoopStructural.interpolators._surfe_wrapper import SurfeRBFInterpolator
 
+try:
+    from LoopStructural.interpolators._surfe_wrapper import SurfeRBFInterpolator
+except ImportError:
+    logger.warning('Can\'t import surfepy - to install "pip install surfe"')
 
-logger.warning(
-    "Using experimental interpolators: P1Interpolator and P2Interpolator"
-)
+logger.warning("Using experimental interpolators: P1Interpolator and P2Interpolator")
 from ._p1interpolator import P1Interpolator
 from ._p2interpolator import P2Interpolator

LoopStructural/interpolators/_discrete_interpolator.py

@@ -705,7 +705,6 @@ def _solve(self, solver="cg", **kwargs):
             self.c[self.region] = self._solve_chol(A, B)
         if solver == "lu":
             logger.info("Solving using scipy LU")
-            print(self.region)
             self.c[self.region] = self._solve_lu(A, B)
         if solver == "pyamg":
             try:

LoopStructural/interpolators/supports/_3d_unstructured_tetra.py

@@ -380,12 +380,12 @@ def get_element_for_location(self, points):
         npts = 0
         npts_step = int(1e4)
         # break into blocks of 10k points
-        while npts<points.shape[0]:
-            
-            cell_index = np.array(self.aabb_grid.position_to_cell_index(points[:npts+npts_step,:])).swapaxes(
-                0, 1
-            )
-            inside = self.aabb_grid.inside(points[:npts+npts_step,:])
+        while npts < points.shape[0]:
+
+            cell_index = np.array(
+                self.aabb_grid.position_to_cell_index(points[: npts + npts_step, :])
+            ).swapaxes(0, 1)
+            inside = self.aabb_grid.inside(points[: npts + npts_step, :])
             global_index = (
                 cell_index[:, 0]
                 + self.aabb_grid.nsteps_cells[None, 0] * cell_index[:, 1]
@@ -423,13 +423,12 @@ def get_element_for_location(self, points):
             c[:, 3] = vd / v
             # inside = np.ones(c.shape[0],dtype=bool)
             mask = np.all(c >= 0, axis=1)
-        
 
-            verts[:npts+npts_step,:,:][row[mask], :, :] = vertices[mask, :, :]
-            bc[:npts+npts_step,:][row[mask], :] = c[mask, :]
-            tetras[:npts+npts_step][row[mask]] = col[mask]
-            inside[:npts+npts_step][row[mask]] = True
-            npts+=npts_step
+            verts[: npts + npts_step, :, :][row[mask], :, :] = vertices[mask, :, :]
+            bc[: npts + npts_step, :][row[mask], :] = c[mask, :]
+            tetras[: npts + npts_step][row[mask]] = col[mask]
+            inside[: npts + npts_step][row[mask]] = True
+            npts += npts_step
         return verts, bc, tetras, inside
 
     def get_element_gradients(self, elements=None):

LoopStructural/modelling/core/geological_model.py

@@ -13,6 +13,7 @@
 
 # if LoopStructural.experimental:
 from LoopStructural.interpolators import P2Interpolator
+
 try:
     from LoopStructural.interpolators import SurfeRBFInterpolator as Surfe
 
@@ -53,6 +54,16 @@
     get_vectors,
 )
 
+intrusions = True
+try:
+    # from LoopStructural.modelling.intrusions import IntrusionNetwork
+    from LoopStructural.modelling.intrusions import IntrusionBuilder
+
+    # from LoopStructural.modelling.intrusions import IntrusionBody
+    from LoopStructural.modelling.intrusions import IntrusionFrameBuilder
+except ImportError as e:
+    print(e)
+    intrusions = False
 from LoopStructural.utils import getLogger, log_to_file
 
 logger = getLogger(__name__)
@@ -1071,7 +1082,143 @@ def create_and_add_folded_fold_frame(
 
         self._add_feature(folded_fold_frame)
 
-        return folded_fold_frame
+        return fold_frame
+
+    def create_and_add_intrusion(
+        self,
+        intrusion_name,
+        intrusion_frame_name,
+        intrusion_lateral_extent_model=None,
+        intrusion_vertical_extent_model=None,
+        intrusion_network_parameters={},
+        lateral_extent_sgs_parameters={},
+        vertical_extent_sgs_parameters={},
+        **kwargs,
+    ):
+        """
+        An intrusion in built in two main steps:
+        (1) Intrusion builder: intrusion builder creates the intrusion structural frame.
+            This object is curvilinear coordinate system of the intrusion constrained with intrusion network points,
+            and flow and inflation measurements (provided by the user).
+            The intrusion network is a representation of the approximated location of roof or floor contact of the intrusion.
+            This object might be constrained using the anisotropies of the host rock if the roof (or floor) contact is not well constrained.
+
+        (2) Intrusion feature: simulation of lateral and vertical extent of intrusion within the model volume.
+            The simulations outcome consist in thresholds distances along the structural frame coordinates
+            that are used to constrained the extent of the intrusion.
+
+        Parameters
+        ----------
+        intrusion_name :  string,
+            name of intrusion feature in model data
+        intrusion_frame_name :  string,
+            name of intrusion frame in model data
+        intrusion_lateral_extent_model = function,
+            geometrical conceptual model for simulation of lateral extent
+        intrusion_vertical_extent_model = function,
+            geometrical conceptual model for simulation of vertical extent
+        intrusion_network_parameters : dictionary, optional
+            contact : string, contact of the intrusion to be used to create the network (roof or floor)
+            type : string, type of algorithm to create the intrusion network (interpolated or shortest path).
+                    Shortest path is recommended when intrusion contact is not well constrained
+            contacts_anisotropies : list of series-type features involved in intrusion emplacement
+            structures_anisotropies : list of fault-type features involved in intrusion emplacement
+            sequence_anisotropies : list of anisotropies to look for the shortest path. It could be only starting and end point.
+        lateral_extent_sgs_parameters = dictionary, optional
+            parameters for sequential gaussian simulation of lateral extent
+        vertical_extent_sgs_parameters = dictionary, optional
+            parameters for sequential gaussian simulation of vertical extent
+
+        kwargs
+
+        Returns
+        -------
+        intrusion feature
+
+        """
+        if intrusions == False:
+            logger.error("Libraries not installed")
+            raise Exception("Libraries not installed")
+
+        intrusion_data = self.data[self.data["feature_name"] == intrusion_name].copy()
+        intrusion_frame_data = self.data[
+            self.data["feature_name"] == intrusion_frame_name
+        ].copy()
+
+        # -- get variables for intrusion frame interpolation
+        gxxgz = kwargs.get("gxxgz", 0)
+        gxxgy = kwargs.get("gxxgy", 0)
+        gyxgz = kwargs.get("gyxgz", 0)
+
+        interpolatortype = kwargs.get("interpolatortype", "FDI")
+        nelements = kwargs.get("nelements", 1e2)
+
+        weights = [gxxgz, gxxgy, gyxgz]
+        interpolator = self.get_interpolator(interpolatortype=interpolatortype)
+
+        intrusion_frame_builder = IntrusionFrameBuilder(
+            interpolator, name=intrusion_frame_name, model=self, **kwargs
+        )
+
+        # -- create intrusion network
+        intrusion_frame_builder.set_intrusion_network_parameters(
+            intrusion_data, intrusion_network_parameters
+        )
+        intrusion_network_geometry = intrusion_frame_builder.create_intrusion_network()
+
+        # -- create intrusion frame using intrusion network points and flow/inflation measurements
+        intrusion_frame_builder.set_intrusion_frame_data(
+            intrusion_frame_data, intrusion_network_geometry
+        )
+
+        ## -- create intrusion frame
+        intrusion_frame_builder.setup(
+            nelements=nelements,
+            w2=weights[0],
+            w1=weights[1],
+            gxygz=weights[2],
+        )
+
+        intrusion_frame = intrusion_frame_builder.frame
+        intrusion_builder = IntrusionBuilder(
+            intrusion_frame, model=self, name=f"{intrusion_frame_name}_feature"
+        )
+        intrusion_builder.lateral_extent_model = intrusion_lateral_extent_model
+        intrusion_builder.vertical_extent_model = intrusion_vertical_extent_model
+
+        # # -- Create intrusion feature
+        # intrusion_feature = IntrusionFeature(
+        #     intrusion_name, model=self
+        # )
+
+        # if intrusion_lateral_extent_model == None:
+        #     logger.error(
+        #         "Specify conceptual model function for intrusion lateral extent"
+        #     )
+        # else:
+        #     intrusion_feature.lateral_extent_model = intrusion_lateral_extent_model
+
+        # if intrusion_vertical_extent_model == None:
+        #     logger.error(
+        #         "Specify conceptual model function for intrusion vertical extent"
+        #     )
+        # else:
+        #     intrusion_feature.vertical_extent_model = intrusion_vertical_extent_model
+
+        # logger.info("setting data for thresholds simulation")
+        intrusion_builder.set_data_for_extent_simulation(intrusion_data)
+        intrusion_builder.build_arguments = {
+            "lateral_extent_sgs_parameters": lateral_extent_sgs_parameters,
+            "vertical_extent_sgs_parameters": vertical_extent_sgs_parameters,
+        }
+        # intrusion_builder.vertical_extent_sgs_parameters = vertical_extent_sgs_parameters
+        # intrusion_builder.lateral_extent_sgs_parameters = lateral_extent_sgs_parameters
+        # intrusion_builder.set_l_sgs_GSLIBparameters(lateral_extent_sgs_parameters)
+        # intrusion_builder.set_g_sgs_GSLIBparameters(vertical_extent_sgs_parameters)
+        intrusion_feature = intrusion_builder.feature
+        self._add_feature(intrusion_feature)
+
+        return intrusion_feature
 
     def _add_faults(self, feature_builder, features=None):
         """Adds all existing faults to a geological feature builder
@@ -1360,12 +1507,12 @@ def create_and_add_fault(
         fault : FaultSegment
             created fault
         """
-        if 'fault_extent' in kwargs and major_axis == None:
-            major_axis = kwargs['fault_extent']
-        if 'fault_influence' in kwargs and minor_axis == None:
-            minor_axis = kwargs['fault_influence']
-        if 'fault_vectical_radius' in kwargs and intermediate_axis == None:
-            intermediate_axis = kwargs['fault_vectical_radius']
+        if "fault_extent" in kwargs and major_axis == None:
+            major_axis = kwargs["fault_extent"]
+        if "fault_influence" in kwargs and minor_axis == None:
+            minor_axis = kwargs["fault_influence"]
+        if "fault_vectical_radius" in kwargs and intermediate_axis == None:
+            intermediate_axis = kwargs["fault_vectical_radius"]
 
         logger.info(f'Creating fault "{fault_surface_data}"')
         logger.info(f"Displacement: {displacement}")

LoopStructural/modelling/fault/fault_function.py

@@ -21,11 +21,11 @@ def __init__(self):
         self.w = None
 
     def add_cstr(self, x, y):
-        self.A.append([x ** 3, x ** 2, x, 1.0])
+        self.A.append([x**3, x**2, x, 1.0])
         self.B.append(y)
 
     def add_grad(self, x, g):
-        self.A.append([3 * x ** 2, 2 * x, 1.0, 0.0])
+        self.A.append([3 * x**2, 2 * x, 1.0, 0.0])
         self.B.append(g)
 
     def add_max(self, max_v):
@@ -44,16 +44,16 @@ def __call__(self, v):
             ATA = A.T @ A
             ATB = A.T @ B
             self.w = np.linalg.lstsq(ATA, ATB, rcond=None)[0]
-        eva = self.w[0] * v ** 3 + self.w[1] * v ** 2 + self.w[2] * v + self.w[3]
+        eva = self.w[0] * v**3 + self.w[1] * v**2 + self.w[2] * v + self.w[3]
         eva[v > self.max_v] = (
-            self.w[0] * self.max_v ** 3
-            + self.w[1] * self.max_v ** 2
+            self.w[0] * self.max_v**3
+            + self.w[1] * self.max_v**2
             + self.w[2] * self.max_v
             + self.w[3]
         )
         eva[v < self.min_v] = (
-            self.w[0] * self.min_v ** 3
-            + self.w[1] * self.min_v ** 2
+            self.w[0] * self.min_v**3
+            + self.w[1] * self.min_v**2
             + self.w[2] * self.min_v
             + self.w[3]
         )

LoopStructural/modelling/features/structural_frame_builder.py

@@ -90,6 +90,7 @@ def __init__(
             ],
             fold=kwargs.get("fold", None),
         )
+        self._frame.builder = self
 
     @property
     def frame(self):

LoopStructural/modelling/input/process_data.py

@@ -141,9 +141,9 @@ def stratigraphic_column(self):
                         "id": unit_id,
                         "colour": self.colours[g],
                     }
-                    if i==0:
+                    if i == 0:
                         stratigraphic_column[name][g]["min"] = -np.inf
-                    if i==len(sg)-1:
+                    if i == len(sg) - 1:
                         stratigraphic_column[name][g]["max"] = np.inf
 
                 unit_id += 1