fix: abutting fault polarity was being calculated inconsistently.

Have changed to use the estimated down_throw direction, if different <90 then use positive, otherwise use negative

Author: Lachlan Grose

LoopStructural/modelling/core/geological_model.py

@@ -235,7 +235,9 @@ def from_processor(cls, processor):
                     model[edge[1]].splay[model[edge[0]].name] = region
                     splay = True
             if splay == False:
-                model[edge[1]].add_abutting_fault(model[edge[0]])
+                model[edge[1]].add_abutting_fault(model[edge[0]],np.abs(processor.stratigraphic_column['faults'][edge[0]]['downthrow_dir']-
+                                                            processor.stratigraphic_column['faults'][edge[1]]['downthrow_dir']) < 90)
+                # model[edge[1]].add_abutting_fault(model[edge[0]])
         for s in processor.stratigraphic_column.keys():
             if s != 'faults':
                 f = model.create_and_add_foliation(s,**processor.foliation_properties[s])

LoopStructural/modelling/fault/fault_segment.py

@@ -2,7 +2,7 @@
 
 from LoopStructural.modelling.fault.fault_function_feature import FaultDisplacementFeature
 from LoopStructural.modelling.fault.fault_function import BaseFault
-from LoopStructural.utils import getLogger
+from LoopStructural.utils import getLogger, NegativeRegion, PositiveRegion
 logger = getLogger(__name__)
 from concurrent.futures import ThreadPoolExecutor
 import numpy as np
@@ -347,19 +347,24 @@ def apply_to_points(self, points):
             newp[mask, :] += g
         return newp
 
-    def add_abutting_fault(self,abutting_fault_feature):
-        pts = self.faultframe[0].builder.data[['X','Y','Z']].to_numpy()#get_value_constraints()
+    def add_abutting_fault(self,abutting_fault_feature,positive=None):
         # check whether the fault is on the hanging wall or footwall of abutting fault
-        
-        def abutting_region(pos):
+        abutting_region = None
+        if positive is None:
+            pts = self.faultframe[0].builder.data[['X','Y','Z']].to_numpy()#get_value_constraints()
             abut_value = np.nanmedian(abutting_fault_feature.evaluate_value(pts))
-            if abut_value > 0:
-                 ## adding the nan check avoids truncating the fault at the edge of the abutting fault bounding box.
-                    ## it makes the assumption that the abutted fault is not drawn across the abutting fault... but this should be ok
-                return np.logical_or(abutting_fault_feature.evaluate_value(pos) > 0, 
-                                        np.isnan(abutting_fault_feature.evaluate_value(pos)))
-            if abut_value < 0:
-                return np.logical_or(abutting_fault_feature.evaluate_value(pos) < 0, 
-                                        np.isnan(abutting_fault_feature.evaluate_value(pos)))
+            positive = abut_value > 0
+        if positive:
+            abutting_region = PositiveRegion(abutting_fault_feature)
+        if positive ==  False:
+            abutting_region = NegativeRegion(abutting_fault_feature)
+            # if positive == True:
+            #      ## adding the nan check avoids truncating the fault at the edge of the abutting fault bounding box.
+            #         ## it makes the assumption that the abutted fault is not drawn across the abutting fault... but this should be ok
+            #     return np.logical_or(abutting_fault_feature.evaluate_value(pos) > 0, 
+            #                             np.isnan(abutting_fault_feature.evaluate_value(pos)))
+            # if positive == False:
+            #     return np.logical_or(abutting_fault_feature.evaluate_value(pos) < 0, 
+            #                             np.isnan(abutting_fault_feature.evaluate_value(pos)))
         self.abut[abutting_fault_feature.name] = abutting_region
         self.faultframe[0].add_region(abutting_region)

LoopStructural/utils/__init__.py

@@ -6,5 +6,5 @@
 from .map2loop import process_map2loop, build_model
 from .helper import get_data_axis_aligned_bounding_box, get_data_bounding_box, get_data_bounding_box_map
 from .helper import get_dip_vector,get_strike_vector, get_vectors, strike_dip_vector
-from .regions import RegionEverywhere
+from .regions import RegionEverywhere, RegionFunction, NegativeRegion, PositiveRegion
 from .exceptions import LoopException, LoopImportError, InterpolatorError, LoopTypeError, LoopValueError

LoopStructural/utils/regions.py

@@ -1,11 +1,33 @@
 import numpy as np
+class BaseRegion:
+    def __init__(self):
+        self.type="BaseRegion"
 
-class RegionEverywhere:
+class RegionEverywhere(BaseRegion):
+    def __init__(self):
+        super().__init__()
+        self.type = "RegionEverywhere"
     def __call__(self,xyz):
         return np.ones(xyz.shape[0],dtype='bool')
 
-class RegionFunction:
+class RegionFunction(BaseRegion):
     def __init__(self,function):
+        super().__init__()
         self.function = function
     def __call__(self,xyz):
-        return self.function(xyz)
+        return self.function(xyz)
+    
+class PositiveRegion:
+    def __init__(self,feature):
+        self.feature = feature
+    def __call__(self,xyz):
+        return np.logical_or(self.feature.evaluate_value(xyz) > 0, 
+                                        np.isnan(self.feature.evaluate_value(xyz)))
+
+class NegativeRegion:
+    def __init__(self,feature):
+        self.feature = feature
+    def __call__(self,xyz):
+        return np.logical_or(self.feature.evaluate_value(xyz) < 0, 
+                                        np.isnan(self.feature.evaluate_value(xyz)))
+