fix: 🐛 adding scaling to splay faults to reduce bubbles

LoopStructural/modelling/core/geological_model.py

@@ -292,20 +292,20 @@ def from_processor(cls, processor):
                 )
                 continue
             splay = False
-            # if "angle" in properties:
-            #     if (
-            #         float(properties["angle"]) < 30
-            #         and np.abs(
-            #             processor.stratigraphic_column["faults"][edge[0]]["dip_dir"]
-            #             - processor.stratigraphic_column["faults"][edge[1]]["dip_dir"]
-            #         )
-            #         < 90
-            #     ):
-            #         # splay
-            #         region = model[edge[1]].builder.add_splay(model[edge[0]])
-
-            #         model[edge[1]].splay[model[edge[0]].name] = region
-            #         splay = True
+            if "angle" in properties:
+                if float(properties["angle"]) < 30:
+                    # if 'dip_dir' in processor.stratigraphic_column["faults"][edge[0]] and 'dip_dir' in processor.stratigraphic_column["faults"][edge[1]]:
+                    #     if np.abs(
+                    #     processor.stratigraphic_column["faults"][edge[0]]["dip_dir"]
+                    #     - processor.stratigraphic_column["faults"][edge[1]]["dip_dir"] > 90
+                    #     ):
+                    #         pass
+                    #     else:
+                    #     # splay
+                    region = model[edge[1]].builder.add_splay(model[edge[0]])
+
+                    model[edge[1]].splay[model[edge[0]].name] = region
+                    splay = True
             if splay == False:
                 model[edge[1]].add_abutting_fault(
                     model[edge[0]],
@@ -319,17 +319,18 @@ def from_processor(cls, processor):
                     )
                     < 90,
                 )
-        for s in processor.stratigraphic_column.keys():
-            if s != "faults":
-                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
-                )
-                model.add_unconformity(f, 0)
-        model.stratigraphic_column = processor.stratigraphic_column
+        if processor.stratigraphy:
+            for s in processor.stratigraphic_column.keys():
+                if s != "faults":
+                    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
+                    )
+                    model.add_unconformity(f, 0)
+            model.stratigraphic_column = processor.stratigraphic_column
         return model
 
     @classmethod
@@ -387,7 +388,7 @@ def dtm(self, dtm):
     def faults(self):
         faults = []
         for f in self.features:
-            if f.type == "fault":
+            if type(f) == FaultSegment:
                 faults.append(f)
 
         return faults
@@ -1608,7 +1609,7 @@ def create_and_add_fault(
             minor_axis=minor_axis,
             major_axis=major_axis,
             intermediate_axis=intermediate_axis,
-            points=kwargs.get("points", True),
+            points=kwargs.get("points", False),
         )
         if "force_mesh_geometry" not in kwargs:
 
@@ -1621,13 +1622,6 @@ def create_and_add_fault(
         fault = fault_frame_builder.frame
         fault.displacement = displacement_scaled
         fault.faultfunction = faultfunction
-        # fault = FaultSegment(
-        #     fault_frame,
-        #     displacement=displacement_scaled,
-        #     faultfunction=faultfunction,
-        #     **kwargs,
-        # )
-        # fault.builder = fault_frame_builder
 
         for f in reversed(self.features):
             if f.type == "unconformity":

LoopStructural/modelling/core/stratigraphic_column.py

@@ -1,4 +1,4 @@
-class StratigraphicColumn:
+class StratigraphicColumn(dict):
     def __init__(self):
         pass
 
@@ -7,3 +7,11 @@ def add_series(self, feature, name, maximum, minimum, colour):
 
     def get_cmap(self):
         pass
+
+    @property
+    def faults(self):
+        return
+
+    @property
+    def stratigraphy(self):
+        return

LoopStructural/modelling/features/geological_feature_builder.py

@@ -330,9 +330,9 @@ def install_gradient_constraint(self):
                 B=B,
             )
 
-    def add_equality_constraints(self, feature, region):
+    def add_equality_constraints(self, feature, region, scalefactor=1.0):
 
-        self._equality_constraints[feature.name] = [feature, region]
+        self._equality_constraints[feature.name] = [feature, region, scalefactor]
         self._up_to_date = False
 
     def install_equality_constraints(self):
@@ -346,7 +346,7 @@ def install_equality_constraints(self):
                 idc = np.arange(0, support.n_nodes)[e[1](support.nodes)]
                 val = e[0].evaluate_value(support.nodes[e[1](support.nodes), :])
                 mask = ~np.isnan(val)
-                self.interpolator.add_equality_constraints(idc[mask], val[mask])
+                self.interpolator.add_equality_constraints(idc[mask], val[mask] * e[2])
             except BaseException as e:
                 logger.error(f"Could not add equality for {self.name}")
                 logger.error(f"Exception: {e}")
@@ -416,6 +416,18 @@ def get_norm_constraints(self):
         else:
             return np.zeros((0, 7))
 
+    def get_orientation_constraints(self):
+        """
+        Get the orientation constraints
+
+        Returns
+        -------
+        numpy array
+        """
+        gradient_constraints = self.get_gradient_constraints()
+        normal_constraints = self.get_norm_constraints()
+        return np.vstack([gradient_constraints, normal_constraints])
+
     def get_interface_constraints(self):
         mask = np.all(
             ~np.isnan(self.data.loc[:, interface_name()].to_numpy(float)), axis=1

LoopStructural/modelling/input/process_data.py

@@ -65,6 +65,7 @@ def __init__(
         self.stratigraphic_order = stratigraphic_order
         self._thicknesses = thicknesses
         self._use_thickness = use_thickness
+        self.stratigraphy = False
         if self.thicknesses is None:
             self._use_thickness = False
         if self._use_thickness is None:
@@ -85,7 +86,13 @@ def __init__(
         # self._fault_dimensions = None
         self._fault_network = None
         self._fault_properties = None
-        self.fault_properties = fault_properties
+        if fault_properties is not None:
+            self.fault_properties = fault_properties
+        elif fault_locations is not None:
+            self.fault_properties = pd.DataFrame(
+                fault_locations["fault_name"].unique(), columns=["name"]
+            ).set_index("name")
+
         if fault_edges is not None and fault_edge_properties is not None:
             self.set_fault_network(fault_edges, fault_edge_properties)  # = fault_graph
         self._fault_stratigraphy = fault_stratigraphy
@@ -124,6 +131,9 @@ def stratigraphic_order(self):
 
     @stratigraphic_order.setter
     def stratigraphic_order(self, stratigraphic_order):
+        if stratigraphic_order is None:
+            logger.warning("No stratigraphic order provided")
+            return
         if isinstance(stratigraphic_order[0][1], list) == False:
             raise TypeError(
                 "Stratigraphic_order must of the format [[('group_name',['unit1','unit2']),('group_name2',['unit3','unit4'])]]"
@@ -132,6 +142,7 @@ def stratigraphic_order(self, stratigraphic_order):
             raise TypeError("Stratigraphic_order must be a list")
         if isinstance(stratigraphic_order[0][1][0], str) == False:
             raise TypeError("Stratigraphic_order elements must be strings")
+        self.stratigraphy = True
         self._stratigraphic_order = stratigraphic_order
 
     @property
@@ -140,6 +151,7 @@ def colours(self):
 
     @colours.setter
     def colours(self, colours):
+
         if colours is None:
             self._colours = {}
             for s in self.stratigraphic_name:
@@ -151,25 +163,26 @@ def colours(self, colours):
     def stratigraphic_column(self):
         stratigraphic_column = {}
         # add stratigraphy into the column
-        unit_id = 0
-        val = self._stratigraphic_value()
-        for name, sg in self._stratigraphic_order:
-            stratigraphic_column[name] = {}
-            for i, g in enumerate(reversed(sg)):
-                if g in self.thicknesses:
-                    stratigraphic_column[name][g] = {
-                        "max": val[g] + self.thicknesses[g],
-                        "min": val[g],
-                        "id": unit_id,
-                        "colour": self.colours[g],
-                    }
-                    if i == 0:
-                        stratigraphic_column[name][g]["min"] = -np.inf
-                    if i == len(sg) - 1:
-                        stratigraphic_column[name][g]["max"] = np.inf
-
-                unit_id += 1
-        # add faults into the column
+        if self.stratigraphy:
+            unit_id = 0
+            val = self._stratigraphic_value()
+            for name, sg in self._stratigraphic_order:
+                stratigraphic_column[name] = {}
+                for i, g in enumerate(reversed(sg)):
+                    if g in self.thicknesses:
+                        stratigraphic_column[name][g] = {
+                            "max": val[g] + self.thicknesses[g],
+                            "min": val[g],
+                            "id": unit_id,
+                            "colour": self.colours[g],
+                        }
+                        if i == 0:
+                            stratigraphic_column[name][g]["min"] = -np.inf
+                        if i == len(sg) - 1:
+                            stratigraphic_column[name][g]["max"] = np.inf
+
+                    unit_id += 1
+            # add faults into the column
         if self.fault_properties is not None:
             stratigraphic_column["faults"] = self.fault_properties.to_dict("index")
         return stratigraphic_column
@@ -249,6 +262,8 @@ def foliation_properties(self):
 
     @foliation_properties.setter
     def foliation_properties(self, foliation_properties):
+        if self.stratigraphic_order == None:
+            return
         if foliation_properties is None:
             for k in self.stratigraphic_column.keys():
                 if k != "faults":
@@ -334,10 +349,16 @@ def fault_properties(self, fault_properties):
                             "avgSlipDirAltitude",
                         ],
                     ] = np.mean(pts, axis=0)
+        if "displacment" not in fault_properties.columns:
+            fault_properties["displacement"] = 0
+            logger.warning(
+                "Fault displacement not provided, setting to 0. \n\
+                This will result in only a fault surface, no displacement on older features"
+            )
+            # faults need a displacement, if we don't know it set as 0 so that the surface is built
+
         self._fault_properties = fault_properties
-        print(list(self.fault_properties.index))
         self._fault_network = FaultNetwork(list(self.fault_properties.index))
-        print(self.fault_network)
 
     @property
     def fault_network(self):
@@ -401,7 +422,9 @@ def vector_scale(self, vector_scale):
     @property
     def stratigraphic_name(self):
         names = []
-        for name, sg in self._stratigraphic_order:
+        if self.stratigraphic_order is None:
+            return names
+        for name, sg in self.stratigraphic_order:
             for g in sg:
                 names.append(g)
         return names

LoopStructural/visualisation/lavavu.py

@@ -91,14 +91,22 @@ def _add_surface(
             # calculate the mode value, just to get the most common value
             surfaceval = np.zeros(vertices.shape[0])
             if isinstance(paint_with, GeologicalFeature):
+                # paint with a geological feature
+                # TODO make sure everything that could be
+                # a feature derives from the same base
+                # class, currently structuralframes and
+                # faults don't work here..
+                # or just change to using __call__
                 surfaceval[:] = paint_with.evaluate_value(
                     self.model.scale(vertices, inplace=False)
                 )
                 surf.values(surfaceval, "paint_with")
             if callable(paint_with):
+                # paint with a callable function e.g. (xyz)->value
                 surfaceval[:] = paint_with(self.model.scale(vertices))
                 surf.values(surfaceval, "paint_with")
             if isinstance(paint_with, (float, int)):
+                # paint with array a constant value
                 surfaceval[:] = paint_with
                 surf.values(surfaceval, "paint_with")
             surf["colourby"] = "paint_with"