Merge pull request #130 from Loop3D/intrusions

Merging intrusions, although they may not work

Author: lachlangrose

.gitignore

@@ -125,6 +125,7 @@ docs/source/LoopStructural.utils.rst
 docs/source/LoopStructural.visualisation.rst
 docs/source/_autosummary/*
 docs/source/_auto_examples/*
+docs/source/auto_examples/*
 examples/*/*.png
 docs/source/test/*
 examples/*.png

LoopStructural/__init__.py

@@ -6,7 +6,8 @@
 
 import logging
 from logging.config import dictConfig
-__all__ = ['GeologicalModel']
+
+__all__ = ["GeologicalModel"]
 import tempfile
 from pathlib import Path
 from .version import __version__

LoopStructural/interpolators/__init__.py

@@ -2,13 +2,24 @@
 Interpolators and interpolation supports
 
 """
-__all__ = ['InterpolatorType','GeologicalInterpolator','DiscreteInterpolator','FiniteDifferenceInterpolator','PiecewiseLinearInterpolator','DiscreteFoldInterpolator','SurfeRBFInterpolator','P1Interpolator','P2Interpolator','TetMesh',
-    'StructuredGrid',
-    'UnStructuredTetMesh',
-    'P1Unstructured2d',
-    'P2Unstructured2d',
-    'StructuredGrid2D',
-    'P2UnstructuredTetMesh',]
+__all__ = [
+    "InterpolatorType",
+    "GeologicalInterpolator",
+    "DiscreteInterpolator",
+    "FiniteDifferenceInterpolator",
+    "PiecewiseLinearInterpolator",
+    "DiscreteFoldInterpolator",
+    "SurfeRBFInterpolator",
+    "P1Interpolator",
+    "P2Interpolator",
+    "TetMesh",
+    "StructuredGrid",
+    "UnStructuredTetMesh",
+    "P1Unstructured2d",
+    "P2Unstructured2d",
+    "StructuredGrid2D",
+    "P2UnstructuredTetMesh",
+]
 from enum import IntEnum
 
 from LoopStructural.utils import getLogger

LoopStructural/interpolators/_discrete_interpolator.py

@@ -85,7 +85,7 @@ def region(self) -> np.ndarray:
         Returns
         -------
         np.ndarray
-            
+
         """
         return self.region_function(self.support.nodes).astype(bool)
 
@@ -352,7 +352,7 @@ def add_inequality_feature(self, feature, lower=True, mask=None):
         Parameters
         ----------
         feature : BaseFeature
-            the feature that will be used to constraint the interpolator 
+            the feature that will be used to constraint the interpolator
         lower : bool, optional
             lower or upper constraint, by default True
         mask : np.ndarray, optional
@@ -819,7 +819,7 @@ def update(self):
             self.setup_interpolator()
             return self._solve(self.solver)
 
-    def evaluate_value(self, evaluation_points:np.ndarray) ->np.ndarray:
+    def evaluate_value(self, evaluation_points: np.ndarray) -> np.ndarray:
         """Evaluate the value of the interpolator at location
 
         Parameters
@@ -842,7 +842,7 @@ def evaluate_value(self, evaluation_points:np.ndarray) ->np.ndarray:
             )
         return evaluated
 
-    def evaluate_gradient(self, evaluation_points:np.ndarray)->np.ndarray:
+    def evaluate_gradient(self, evaluation_points: np.ndarray) -> np.ndarray:
         """
         Evaluate the gradient of the scalar field at the evaluation points
         Parameters

LoopStructural/interpolators/_finite_difference_interpolator.py

@@ -17,7 +17,6 @@
 
 
 class FiniteDifferenceInterpolator(DiscreteInterpolator):
-
     def __init__(self, grid):
         """
         Finite difference interpolation on a regular cartesian grid
@@ -224,9 +223,7 @@ def add_inequality_constraints(self, w=1.0):
                         value constraints not added: outside of model bounding box"
                 )
 
-    def add_interface_constraints(
-        self, w=1.0
-    ):  
+    def add_interface_constraints(self, w=1.0):
         """
         Adds a constraint that defines all points
         with the same 'id' to be the same value

LoopStructural/interpolators/_geological_interpolator.py

@@ -71,7 +71,7 @@ def set_property_name(self, name):
         """
         self.propertyname = name
 
-    def set_value_constraints(self, points : np.ndarray):
+    def set_value_constraints(self, points: np.ndarray):
         """
 
         Parameters
@@ -85,12 +85,12 @@ def set_value_constraints(self, points : np.ndarray):
 
         """
         if points.shape[1] < 4:
-            raise ValueError('Value points must at least have X,Y,Z,val')
+            raise ValueError("Value points must at least have X,Y,Z,val")
         self.data["value"] = points
         self.n_i = points.shape[0]
         self.up_to_date = False
 
-    def set_gradient_constraints(self, points : np.ndarray):
+    def set_gradient_constraints(self, points: np.ndarray):
         """
 
         Parameters
@@ -104,7 +104,7 @@ def set_gradient_constraints(self, points : np.ndarray):
 
         """
         if points.shape[1] < 7:
-            raise ValueError('Gradient constraints must at least have X,Y,Z,gx,gy,gz')
+            raise ValueError("Gradient constraints must at least have X,Y,Z,gx,gy,gz")
         self.n_g = points.shape[0]
         self.data["gradient"] = points
         self.up_to_date = False
@@ -123,7 +123,7 @@ def set_normal_constraints(self, points: np.ndarray):
 
         """
         if points.shape[1] < 7:
-            raise ValueError('Nonrmal constraints must at least have X,Y,Z,nx,ny,nz')
+            raise ValueError("Nonrmal constraints must at least have X,Y,Z,nx,ny,nz")
         self.n_n = points.shape[0]
         self.data["normal"] = points
         self.up_to_date = False
@@ -142,7 +142,7 @@ def set_tangent_constraints(self, points: np.ndarray):
 
         """
         if points.shape[1] < 7:
-            raise ValueError('Tangent constraints must at least have X,Y,Z,tx,ty,tz')
+            raise ValueError("Tangent constraints must at least have X,Y,Z,tx,ty,tz")
         self.data["tangent"] = points
         self.up_to_date = False
 
@@ -200,7 +200,7 @@ def get_data_locations(self):
         numpy array
             Nx3 - X,Y,Z location of all data points
         """
-        return np.vstack([d for d in self.data.values()[:,:3]])
+        return np.vstack([d for d in self.data.values()[:, :3]])
 
     def get_interface_constraints(self):
         """Get the location of interface constraints

LoopStructural/interpolators/supports/_3d_structured_grid.py

@@ -314,7 +314,7 @@ def evaluate_value(self, evaluation_points, property_array):
         # print(idc[inside,:], self.n_nodes,inside)
 
         if idc.shape[0] != inside.shape[0]:
-            raise ValueError('index does not match number of nodes')
+            raise ValueError("index does not match number of nodes")
         v = np.zeros(idc.shape)
         v[:, :] = np.nan
 
@@ -364,19 +364,32 @@ def evaluate_gradient(self, evaluation_points, property_array):
         # indices = np.array([self.position_to_cell_index(evaluation_points)])
         # idc = self.global_indicies(indices.swapaxes(0,1))
         # print(idc)
-        if np.max(idc[inside,:]) > property_array.shape[0]:
+        if np.max(idc[inside, :]) > property_array.shape[0]:
             cix, ciy, ciz = self.position_to_cell_index(evaluation_points)
             if np.all(cix[inside] < self.nsteps_cells[0]) == False:
-                print(evaluation_points[inside,:][cix[inside] < self.nsteps_cells[0],0],self.origin[0],self.maximum[0])
+                print(
+                    evaluation_points[inside, :][cix[inside] < self.nsteps_cells[0], 0],
+                    self.origin[0],
+                    self.maximum[0],
+                )
             if np.all(ciy[inside] < self.nsteps_cells[1]) == False:
-                print(evaluation_points[inside,:][ciy[inside] < self.nsteps_cells[1],1],self.origin[1],self.maximum[1])
+                print(
+                    evaluation_points[inside, :][ciy[inside] < self.nsteps_cells[1], 1],
+                    self.origin[1],
+                    self.maximum[1],
+                )
             if np.all(ciz[inside] < self.nsteps_cells[2]) == False:
-                print(ciz[inside],self.nsteps_cells[2])
-                print(self.step_vector, self.nsteps_cells,self.nsteps)
-                print(evaluation_points[inside,:][~(ciz[inside] < self.nsteps_cells[2]),2],self.origin[2],self.maximum[2])
-            
+                print(ciz[inside], self.nsteps_cells[2])
+                print(self.step_vector, self.nsteps_cells, self.nsteps)
+                print(
+                    evaluation_points[inside, :][
+                        ~(ciz[inside] < self.nsteps_cells[2]), 2
+                    ],
+                    self.origin[2],
+                    self.maximum[2],
+                )
 
-            raise ValueError('index does not match number of nodes')
+            raise ValueError("index does not match number of nodes")
         T[inside, 0, :] *= property_array[idc[inside, :]]
         T[inside, 1, :] *= property_array[idc[inside, :]]
         T[inside, 2, :] *= property_array[idc[inside, :]]

LoopStructural/interpolators/supports/_3d_structured_tetra.py

@@ -14,7 +14,9 @@
 class TetMesh(BaseStructuredSupport):
     """ """
 
-    def __init__(self, origin=np.zeros(3), nsteps=np.ones(3)*10, step_vector=np.ones(3)):
+    def __init__(
+        self, origin=np.zeros(3), nsteps=np.ones(3) * 10, step_vector=np.ones(3)
+    ):
         BaseStructuredSupport.__init__(self, origin, nsteps, step_vector)
 
         self.tetra_mask_even = np.array(

LoopStructural/modelling/__init__.py

@@ -2,7 +2,13 @@
 Geological modelling classes and functions
 
 """
-__all__ = ['GeologicalModel', 'ProcessInputData', 'LavaVuModelViewer', 'Map2LoopProcessor','LoopProjectfileProcessor']
+__all__ = [
+    "GeologicalModel",
+    "ProcessInputData",
+    "LavaVuModelViewer",
+    "Map2LoopProcessor",
+    "LoopProjectfileProcessor",
+]
 from LoopStructural.utils import getLogger
 from LoopStructural.utils import LoopImportError
 

LoopStructural/modelling/core/geological_model.py

@@ -398,13 +398,13 @@ def dtm(self):
     def dtm(self, dtm):
         """Set a dtm to the model.
         The dtm is a function that can be called for dtm(xy) where xy is
-        a numpy array of xy locations. The function will return an array of 
+        a numpy array of xy locations. The function will return an array of
         z values corresponding to the elevation at xy.
 
         Parameters
         ----------
         dtm : callable
-            
+
         """
         if not callable(dtm):
             raise BaseException(
@@ -431,6 +431,14 @@ def series(self):
                 series.append(f)
         return series
 
+    @property
+    def intrusions(self):
+        intrusions = []
+        for f in self.features:
+            if f.type == "intrusion":
+                intrusions.append(f)
+        return intrusions
+
     @property
     def faults_displacement_magnitude(self):
         displacements = []
@@ -500,7 +508,7 @@ def _add_feature(self, feature):
             self._add_unconformity_above(feature)
         feature.model = self
 
-    def data_for_feature(self, feature_name:str)-> pd.DataFrame:
+    def data_for_feature(self, feature_name: str) -> pd.DataFrame:
         return self.data.loc[self.data["feature_name"] == feature_name, :]
 
     @property
@@ -1008,7 +1016,7 @@ def create_and_add_folded_fold_frame(
         fold_frame : StructuralFrame, optional
             the fold frame for the fold if not specified uses last feature added
 
-        kwargs : dict  
+        kwargs : dict
             parameters passed to child functions
 
         Returns
@@ -1080,9 +1088,14 @@ def create_and_add_intrusion(
         intrusion_network_parameters={},
         lateral_extent_sgs_parameters={},
         vertical_extent_sgs_parameters={},
+        geometric_scaling_parameters={},
+        faults=None,  # LG seems unused?
         **kwargs,
     ):
         """
+
+        Note
+        -----
         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,
@@ -1126,7 +1139,7 @@ def create_and_add_intrusion(
         intrusion feature
 
         """
-        if intrusions == False:
+        if intrusions is False:
             logger.error("Libraries not installed")
             raise Exception("Libraries not installed")
 
@@ -1149,6 +1162,7 @@ def create_and_add_intrusion(
         intrusion_frame_builder = IntrusionFrameBuilder(
             interpolator, name=intrusion_frame_name, model=self, **kwargs
         )
+        intrusion_frame_builder.post_intrusion_faults = faults  # LG unused?
 
         # -- create intrusion network
         intrusion_frame_builder.set_intrusion_network_parameters(
@@ -1172,8 +1186,9 @@ def create_and_add_intrusion(
 
         intrusion_frame = intrusion_frame_builder.frame
 
-        # -- create intrusion builder to simulate distance thresholds
-        # along frame coordinates
+        self._add_faults(intrusion_frame_builder, features=faults)
+
+        # -- create intrusion builder to simulate distance thresholds along frame coordinates
         intrusion_builder = IntrusionBuilder(
             intrusion_frame, model=self, name=f"{intrusion_name}_feature"
         )
@@ -1185,8 +1200,10 @@ def create_and_add_intrusion(
         intrusion_builder.build_arguments = {
             "lateral_extent_sgs_parameters": lateral_extent_sgs_parameters,
             "vertical_extent_sgs_parameters": vertical_extent_sgs_parameters,
+            "geometric_scaling_parameters": geometric_scaling_parameters,
         }
         intrusion_feature = intrusion_builder.feature
+        # self._add_faults(intrusion_feature, features = faults)
         self._add_feature(intrusion_feature)
 
         return intrusion_feature
@@ -1647,7 +1664,7 @@ def create_and_add_fault(
 
         return fault
 
-    def rescale(self, points:np.ndarray, inplace:bool=True) ->np.ndarray:
+    def rescale(self, points: np.ndarray, inplace: bool = True) -> np.ndarray:
         """
         Convert from model scale to real world scale - in the future this
         should also do transformations?
@@ -1669,7 +1686,7 @@ def rescale(self, points:np.ndarray, inplace:bool=True) ->np.ndarray:
         points += self.origin
         return points
 
-    def scale(self, points:np.ndarray, inplace:bool=True)->np.ndarray:
+    def scale(self, points: np.ndarray, inplace: bool = True) -> np.ndarray:
         """Take points in UTM coordinates and reproject
         into scaled model space