Merge branch 'intrusions' of github.com:Loop3D/LoopStructural into intrusions

Author: Lachlan Grose

LoopStructural/__init__.py

@@ -14,7 +14,7 @@
 
 ch = logging.StreamHandler()
 formatter = logging.Formatter(
-    "%(asctime)s ~ %(name)-12s ~ %(levelname)-10s ~ %(message)s"
+    "%(levelname)s: %(asctime)s: %(filename)s:%(lineno)d -- %(message)s"
 )
 ch.setFormatter(formatter)
 ch.setLevel(logging.WARNING)

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,17 @@ 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/discrete_interpolator.py

@@ -361,13 +361,17 @@ def build_matrix(self, square=True, damp=0.0):
         # which is not ideal.
         max_weight = 0
         for c in self.constraints.values():
+            if len(c["w"]) == 0:
+                continue
             if c["w"].max() > max_weight:
                 max_weight = c["w"].max()
         a = []
         b = []
         rows = []
         cols = []
         for c in self.constraints.values():
+            if len(c["w"]) == 0:
+                continue
             aa = (c["A"] * c["w"][:, None] / max_weight).flatten()
             b.extend((c["B"] * c["w"] / max_weight).tolist())
             mask = aa == 0

LoopStructural/interpolators/finite_difference_interpolator.py

@@ -177,7 +177,8 @@ def add_vaue_constraints(self, w=1.0):
                 w=w * points[inside, 4],
                 name="value",
             )
-
+            if np.sum(inside)<=0:
+                logger.warning(f"{self.propertyname}: {np.sum(~inside)} value constraints not added: outside of model bounding box")
     def add_interface_constraints(
         self, w=1.0
     ):  # for now weight all value points the same
@@ -242,40 +243,7 @@ def add_interface_constraints(
                     name="interface_{}".format(unique_id),
                 )
 
-        # if points.shape[0] > 1:
-        #     node_idx, inside = self.support.position_to_cell_corners(
-        #                     points[:, :3])
-        #     # print(points[inside,:].shape)
-
-        #     gi = np.zeros(self.support.n_nodes)
-        #     gi[:] = -1
-        #     gi[self.region] = np.arange(0, self.nx)
-        #     idc = np.zeros(node_idx.shape)
-        #     idc[:] = -1
-
-        #     idc[inside, :] = gi[node_idx[inside, :]]
-        #     inside = np.logical_and(~np.any(idc == -1, axis=1), inside)
-        #     a = self.support.position_to_dof_coefs(points[inside, :3]).T
-        #     # create oversided array for storing constraints
-        #     A = np.zeros((a.shape[0]*a.shape[0],a.shape[1]*2))
-        #     interface_idc = np.zeros((a.shape[0]*a.shape[0],a.shape[1]*2),dtype=int)
-        #     interface_idc[:] = -1
-        #     c_i = 0
-
-        #     for i in np.unique(points[np.logical_and(~np.isnan(points[:,3]),inside),3]):
-        #         mask = points[inside,3] == i
-        #         for p1 in range(points[inside][mask].shape[0]):
-        #             for p2 in range(p1+1,points[inside][mask].shape[0]):
-        #                 A[c_i,:8] = a[mask][p1,:]
-        #                 A[c_i,8:] -= a[mask][p2,:]
-        #                 interface_idc[c_i,:8] = idc[inside,:][mask,:][p1,:]
-        #                 interface_idc[c_i,8:] = idc[inside,:][mask,:][p2,:]
-        #                 c_i+=1
-        #     outside = ~np.any(interface_idc == -1, axis=1)
-
-        #     self.add_constraints_to_least_squares(A[outside,:] * w,
-        #                                             np.zeros(A[outside,:].shape[0]),
-        #                                             interface_idc[outside, :], name='interface')
+       
 
     def add_gradient_constraints(self, w=1.0):
         """
@@ -310,7 +278,7 @@ def add_gradient_constraints(self, w=1.0):
                 vertices,
                 T,
                 elements,
-                inside,
+                inside_,
             ) = self.support.get_element_gradient_for_location(points[inside, :3])
             # normalise constraint vector and scale element matrix by this
             norm = np.linalg.norm(points[:, 3:6], axis=1)
@@ -327,7 +295,8 @@ def add_gradient_constraints(self, w=1.0):
             self.add_constraints_to_least_squares(
                 A, B, idc[inside, :], w=w * self.vol, name="gradient"
             )
-
+            if np.sum(inside)<=0:
+                logger.warning(f"{self.propertyname}: {np.sum(~inside)} norm constraints not added: outside of model bounding box")
     def add_norm_constraints(self, w=1.0):
         """
         Add constraints to control the norm of the gradient of the scalar field
@@ -361,7 +330,7 @@ def add_norm_constraints(self, w=1.0):
                 vertices,
                 T,
                 elements,
-                inside,
+                inside_,
             ) = self.support.get_element_gradient_for_location(points[inside, :3])
             # T*=np.product(self.support.step_vector)
             # T/=self.support.step_vector[0]
@@ -387,6 +356,8 @@ def add_norm_constraints(self, w=1.0):
                 w=w * self.vol,
                 name="norm",
             )
+            if np.sum(inside)<=0:
+                logger.warning(f"{self.propertyname}: {np.sum(~inside)} norm constraints not added: outside of model bounding box")
 
     def add_gradient_orthogonal_constraints(self, points, vector, w=1.0, B=0):
         """
@@ -427,7 +398,7 @@ def add_gradient_orthogonal_constraints(self, points, vector, w=1.0, B=0):
                 vertices,
                 T,
                 elements,
-                inside,
+                inside_,
             ) = self.support.get_element_gradient_for_location(points[inside, :3])
             T /= norm[:, None, None]
 
@@ -437,6 +408,8 @@ def add_gradient_orthogonal_constraints(self, points, vector, w=1.0, B=0):
             self.add_constraints_to_least_squares(
                 A, B, idc[inside, :], w=w * self.vol, name="gradient orthogonal"
             )
+            if np.sum(inside)<=0:
+                logger.warning(f"{self.propertyname}: {np.sum(~inside)} gradient constraints not added: outside of model bounding box")
 
     def add_regularisation(self, operator, w=0.1):
         """

LoopStructural/interpolators/structured_tetra.py

@@ -451,11 +451,15 @@ def get_element_gradient_for_location(self, pos):
                 ],
             ]
         )
+        # m[~inside,:,:] = np.nan
         I = np.array(
             [[-1.0, 1.0, 0.0, 0.0], [-1.0, 0.0, 1.0, 0.0], [-1.0, 0.0, 0.0, 1.0]]
         )
         m = np.swapaxes(m, 0, 2)
-        element_gradients = np.linalg.inv(m)
+        element_gradients = np.zeros_like(m)
+        element_gradients[:] = np.nan
+        element_gradients[inside,:,:] = np.linalg.inv(m[inside,:,:])
+        # element_gradients = np.linalg.inv(m)
 
         element_gradients = element_gradients.swapaxes(1, 2)
         element_gradients = element_gradients @ I

LoopStructural/modelling/core/geological_model.py

@@ -231,7 +231,8 @@ def from_map2loop_directory(
             the created geological model and a dictionary of the map2loop data
         """
         from LoopStructural.modelling.input.map2loop_processor import Map2LoopProcessor
-
+        log_to_file(f"{m2l_directory}/loopstructural_log.txt")
+        logger.info('Creating model from m2l directory')
         processor = Map2LoopProcessor(m2l_directory, use_thickness)
         processor._gradient = gradient
         processor.vector_scale = vector_scale
@@ -263,9 +264,11 @@ def from_map2loop_directory(
 
     @classmethod
     def from_processor(cls, processor):
+        logger.info('Creating model from processor')
         model = GeologicalModel(processor.origin, processor.maximum)
         model.data = processor.data
         for i in processor.fault_network.faults:
+            logger.info(f"Adding fault {i}")
             model.create_and_add_fault(
                 i,
                 **processor.fault_properties.to_dict("index")[i],
@@ -311,6 +314,7 @@ def from_processor(cls, processor):
                 faults = None
                 if processor.fault_stratigraphy is not None:
                     faults = processor.fault_stratigraphy[s]
+                logger.info(f'Adding foliation {s}')
                 f = model.create_and_add_foliation(
                     s, **processor.foliation_properties[s], faults=faults
                 )
@@ -1073,6 +1077,11 @@ def create_and_add_intrusion(
         **kwargs,
     ):
         """
+        An intrusion in built in two main steps:
+        (1) Intrusion network and intrusion frame: the algorithm first identify the intrusion network, which is a set of points
+        representing the roof or floor contact of the intrusion. Then this set of points is used to contraint the main structural 
+        direction of the structural frame.
+        (2) Intrusion body: simulation of lateral and vertical extent of intrusion, using parameterization provided by the structural frame
 
         Parameters
         ----------
@@ -1130,28 +1139,50 @@ def create_and_add_intrusion(
         INet.build(**kwargs)
 
         # Create intrusion frame, using intrusion network points, propagation and inflation direction
+        if "gxxgz" in kwargs: # weight for orthogonality constraint between coord 0 and coord 2
+                gxxgz = kwargs["gxxgz"]
+        else:
+                gxxgz = 0
 
-        gxxgz = 0  # weight for orthogonality constraint between coord 0 and coord 2
-        gxxgy = 0  # weight for orthogonality constraint between coord 0 and coord 1
-        gyxgz = 0  # weight for orthogonality constraint between coord 1 and coord 2
+        if "gxxgy" in kwargs: # weight for orthogonality constraint between coord 0 and coord 1
+                gxxgy = kwargs["gxxgy"]
+        else:
+                gxxgy = 0
+        
+        if "gyxgz" in kwargs: # weight for orthogonality constraint between coord 1 and coord 2
+                gyxgz = kwargs["gyxgz"]
+        else:
+                gyxgz = 0
+        
+        if "interpolatortype" in kwargs: # number of elements for interpolation
+                interpolatortype = kwargs["interpolatortype"]
+        else:
+                interpolatortype = 'FDI'
+        
+        if "nelements" in kwargs: # number of elements for interpolation
+                nelements = kwargs["nelements"]
+        else:
+                nelements = 1e2
 
         weights = [gxxgz, gxxgy, gyxgz]
-        # reg=np.array([1,0.5,1])
         logger.info("building intrusion frame")
-        interpolator = self.get_interpolator(interpolatortype="FDI")
+        interpolator = self.get_interpolator(interpolatortype=interpolatortype)
         frame_data = self.data[self.data["feature_name"] == intrusion_frame_name].copy()
         IFrame_builder = IntrusionBuilder(
-            interpolator, model=self, feature_name=intrusion_frame_name
+            interpolator, 
+            model=self, 
+            name=intrusion_frame_name
         )
         IFrame_builder.set_data(frame_data, INet.intrusion_network_outcome)
-        IFrame = IFrame_builder.setup(
-            nelements=1e2,
-            solver="lu",
-            gxxgz=weights[0],
-            gxxgy=weights[1],
+        IFrame_builder.setup(
+            nelements = nelements,
+            # solver = solver,
+            w2=weights[0],
+            w1=weights[1],
             gxygz=weights[2],
         )
 
+        IFrame = IFrame_builder.frame
         # Create intrusion feature
         intrusion_feature = IntrusionFeature(
             intrusion_name, structural_frame=IFrame, model=self
@@ -1171,14 +1202,17 @@ def create_and_add_intrusion(
         intrusion_feature.set_intrusion_frame(IFrame)
         intrusion_feature.set_intrusion_body(IBody)
 
+        #set data for simulations
+        logger.info("setting data for lateral thresholds simulation")
+        IBody.set_data_for_s_simulation()
+        logger.info("setting data for vertical thresholds simulation")
+        IBody.set_data_for_g_simulation()
+
         if intrusion_lateral_extent_model == None:
             logger.error(
                 "Specify conceptual model function for intrusion lateral extent"
             )
-
-        else:
-            logger.info("setting data for lateral thresholds simulation")
-            IBody.set_data_for_s_simulation()
+        else:            
             IBody.set_lateral_extent_conceptual_model(intrusion_lateral_extent_model)
             IBody.set_s_simulation_GSLIBparameters(lateral_extent_sgs_parameters)
             IBody.make_s_simulation_variogram(lateral_extent_sgs_parameters)
@@ -1192,10 +1226,7 @@ def create_and_add_intrusion(
             logger.error(
                 "Specify conceptual model function for intrusion vertical extent"
             )
-
-        else:
-            logger.info("setting data for vertical thresholds simulation")
-            IBody.set_data_for_g_simulation()
+        else:           
             IBody.set_vertical_extent_conceptual_model(intrusion_vertical_extent_model)
             IBody.set_g_simulation_GSLIBparameters(vertical_extent_sgs_parameters)
             IBody.make_g_simulation_variogram(vertical_extent_sgs_parameters)
@@ -1493,6 +1524,18 @@ def create_and_add_fault(
         fault : FaultSegment
             created fault
         """
+        logger.info(f'Creating fault "{fault_surface_data}"')
+        logger.info(f'Displacement: {displacement}')
+        logger.info(f'Tolerance: {tol}')
+        logger.info(f'Fault function: {faultfunction}')
+        logger.info(f'Fault slip vector: {fault_slip_vector}')
+        logger.info(f'Fault center: {fault_center}')
+        logger.info(f'Major axis: {major_axis}')
+        logger.info(f'Minor axis: {minor_axis}')
+        logger.info(f'Intermediate axis: {intermediate_axis}')
+        for k, v in kwargs.items():
+            logger.info(f'{k}: {v}')
+        
         if tol is None:
             tol = self.tol
         self.parameters["features"].append(
@@ -1523,6 +1566,12 @@ def create_and_add_fault(
         fault_frame_data = self.data[
             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"])
         )
@@ -1562,9 +1611,7 @@ def create_and_add_fault(
             # if we haven't defined a fault centre take the center of mass for lines assocaited with
             # the fault trace
             if (
-                "centreEasting" in kwargs
-                and "centreNorthing" in kwargs
-                and "centreAltitude" in kwargs
+                ~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(
@@ -1598,7 +1645,7 @@ def create_and_add_fault(
             minor_axis=minor_axis,
             major_axis=major_axis,
             intermediate_axis=intermediate_axis,
-            points=kwargs.get("points", None),
+            points=kwargs.get("points", True),
         )
 
         # if minor_axis == None or major_axis == None or intermediate_axis == None: