fix: faults are more generic

Only normal vector and trace required
BaseFault is now default fault function.
Fault slip is estimated to be dip vector if not provided
Ellipsoid estimated from fault length

LoopStructural/modelling/core/geological_model.py

@@ -31,7 +31,7 @@
 
 from LoopStructural.utils.exceptions import LoopException, InterpolatorError
 from LoopStructural.utils.helper import (all_heading, gradient_vec_names,
-                                         strike_dip_vector)
+                                         strike_dip_vector, get_vectors)
 
 from LoopStructural.utils import getLogger, log_to_file
 logger = getLogger(__name__)
@@ -1177,7 +1177,7 @@ def create_and_add_fault(self,
                             major_axis = None, 
                             minor_axis = None, 
                             intermediate_axis = None,
-                            faultfunction=None, 
+                            faultfunction='BaseFault', 
                             **kwargs):
         """
         Parameters
@@ -1226,6 +1226,10 @@ def create_and_add_fault(self,
                 fault_slip_vector = np.array([kwargs['avgSlipDirEasting'],kwargs['avgSlipDirNorthing'],kwargs['avgSlipDirAltitude']],dtype=float)
             else:
                 fault_slip_vector = fault_frame_data.loc[mask,['gx','gy','gz']].mean(axis=0).to_numpy()
+        if np.any(np.isnan(fault_slip_vector)):
+            logger.warning("Fault slip vector is nan, estimating from fault normal")
+            strike_vector, dip_vector = get_vectors(fault_normal_vector[None,:])
+            fault_slip_vector = dip_vector[:,0]
         if fault_center is not None:
             fault_center = self.scale(fault_center,inplace=False)
         if fault_center is None:

LoopStructural/modelling/fault/fault_builder.py

@@ -67,12 +67,22 @@ def create_data_from_geometry(self,
         intermediate_axis : double
             fault volume radius in the slip direction
         """
+        if major_axis is None:
+            fault_trace = data.loc[np.logical_and(data['coord']==0,data['val']==0),['X','Y']].to_numpy()
+            distance = np.linalg.norm(fault_trace[:,None,:]-fault_trace[None,:,:],axis=2)
+            major_axis = np.max(distance)
+            logger.warning('Fault major axis using map length: {}'.format(major_axis))
+
+        if minor_axis is None:
+            minor_axis = major_axis /2.
+        if intermediate_axis is None:
+            intermediate_axis = major_axis
         normal_vector/=np.linalg.norm(normal_vector)
         slip_vector/=np.linalg.norm(slip_vector)
         # check if slip vector is inside fault plane, if not project onto fault plane
-        if not np.isclose(normal_vector @ slip_vector, 0):
-            logger.info("{} : projecting slip vector onto fault plane".format(self.name))
-            slip_vector = np.cross(normal_vector, np.cross(slip_vector ,normal_vector))
+        # if not np.isclose(normal_vector @ slip_vector, 0):
+        #     logger.info("{} : projecting slip vector onto fault plane".format(self.name))
+        #     slip_vector = np.cross(normal_vector, np.cross(slip_vector ,normal_vector))
         strike_vector = np.cross(normal_vector,slip_vector)
         fault_edges = np.zeros((2,3))
         fault_tips = np.zeros((2,3))

LoopStructural/modelling/features/geological_feature_builder.py

@@ -309,8 +309,9 @@ def install_equality_constraints(self):
                 mask = ~np.isnan(val)
                 self.interpolator.add_equality_constraints(
                     idc[mask], val[mask])
-            except:
-                logger.error("Could not add equality")
+            except as e:
+                logger.error("Could not add equality for {}".format(self.name))
+                logger.error("Exception: {}".format(e))
 
 
     def get_value_constraints(self):

LoopStructural/modelling/fold/fold_builder.py

@@ -75,7 +75,7 @@ def set_fold_limb_rotation(self):
             fold_limb_rotation.fit_fourier_series(wl=l_wl,**kwargs)
         self.fold.fold_limb_rotation = fold_limb_rotation
 
-    def build(self, data_region=None, **kwargs):
+    def build(self, data_region=None, constrained=True,**kwargs):
         """[summary]
 
         Parameters
@@ -87,7 +87,7 @@ def build(self, data_region=None, **kwargs):
         # gradient not norm, to prevent issues with fold norm constraint
         # TODO folding norm constraint should be minimising the difference in norm
         # not setting the norm
-        self.add_data_to_interpolator(constrained=True)
+        self.add_data_to_interpolator(constrained=constrained)
         if self.fold.foldframe[0].is_valid() == False:
             raise InterpolatorError("Fold frame main coordinate is not valid")
         self.set_fold_axis()

LoopStructural/visualisation/model_plotter.py

@@ -614,7 +614,7 @@ def add_vector_field(self, geological_feature, **kwargs):
         # normalise
         mask = ~np.any(np.isnan(vector), axis=1)
         vector[mask, :] /= np.linalg.norm(vector[mask, :], axis=1)[:, None]
-        self._add_vector_marker(locations, vector, name, **kwargs)
+        self._add_vector_marker(self.model.rescale(locations,inplace=False), vector, name, **kwargs)
 
         return