Merge branch 'master' of https://github.com/Loop3D/LoopStructural

Author: lachlangrose

LoopStructural/modelling/core/geological_model.py

@@ -253,7 +253,10 @@ def _add_feature(self, feature):
         self._add_domain_fault_above(feature)
         self._add_unconformity_above(feature)
         feature.set_model(self)
-
+    
+    def data_for_feature(self,feature_name):
+        return self.data.loc[self.data['feature_name'] == series_surface_data,:]
+        
     def set_model_data(self, data):
         """
         Set the data array for the model

LoopStructural/modelling/fault/fault_segment.py

@@ -138,6 +138,10 @@ def evaluate(self, locations):
 
         return self.faultframe.features[0].evaluate_value(locations) > 0
 
+    def inside_volume(self,locations):
+        v = self.faultframe.evaluate_value(locations)
+        return np.all(np.logical_and(v > -1,v<1),axis=1)
+
     def evaluate_value(self, locations):
         """
         Return the value of the fault surface scalar field
@@ -192,6 +196,27 @@ def evaluate_gradient(self, locations):
         # need to scale with fault displacement
         return self.faultframe[1].evaluate_gradient(locations[mask, :])
 
+    def evaluate_displacement(self, points):
+        newp = np.copy(points).astype(float)
+        # evaluate fault function for all points then define mask for only points affected by fault
+        with ThreadPoolExecutor(max_workers=8) as executor:
+            # all of these operations should be independent so just run as different threads
+            gx_future = executor.submit(self.faultframe.features[0].evaluate_value, newp)
+            gy_future = executor.submit(self.faultframe.features[1].evaluate_value, newp)
+            gz_future = executor.submit(self.faultframe.features[2].evaluate_value, newp)
+            gx = gx_future.result()
+            gy = gy_future.result()
+            gz = gz_future.result()
+        d = np.zeros(gx.shape)
+        mask = np.logical_and(~np.isnan(gx),~np.isnan(gy))
+        mask = np.logical_and(mask,~np.isnan(gz))
+        d[~mask] = 0
+        gx_mask = np.zeros_like(mask,dtype=bool)
+        gx_mask[mask] = gx[mask] > 0
+        d[gx_mask] = 1.
+        if self.faultfunction is not None:
+            d[mask] = self.faultfunction(gx[mask], gy[mask], gz[mask])
+        return d
     def apply_to_points(self, points):
         """
         Unfault the array of points

LoopStructural/modelling/features/geological_feature_builder.py

@@ -225,7 +225,7 @@ def get_gradient_constraints(self):
                 mask, xyz_names() + gradient_vec_names() + weight_name(
                 )].to_numpy()
         else:
-            return np.zeros(0, 7)
+            return np.zeros((0, 7))
 
     def get_tangent_constraints(self):
         """
@@ -242,7 +242,7 @@ def get_tangent_constraints(self):
                 mask, xyz_names() + tangent_vec_names() + weight_name(
                 )].to_numpy()
         else:
-            return np.zeros(0, 7)
+            return np.zeros((0, 7))
 
     def get_norm_constraints(self):
         """
@@ -260,15 +260,15 @@ def get_norm_constraints(self):
 
                 )].to_numpy()
         else:
-            return np.zeros(0,7)
+            return np.zeros((0,7))
 
     def get_interface_constraints(self):
         mask = np.all(~np.isnan(self.data.loc[:, interface_name()].to_numpy()), axis=1)
         if mask.shape[0] > 0:
             return self.data.loc[
                 mask, xyz_names() + interface_name() + weight_name()].to_numpy()
         else:
-            return np.zeros(0,5)
+            return np.zeros((0,5))
 
     def get_data_locations(self):
         """

LoopStructural/modelling/features/structural_frame.py

@@ -2,7 +2,7 @@
 Structural frames
 """
 import logging
-
+import numpy as np
 logger = logging.getLogger(__name__)
 
 
@@ -86,7 +86,7 @@ def set_data(self, data):
         """
         self.data = data
 
-    def evaluate_value(self, evaluation_points, i=None):
+    def evaluate_value(self, evaluation_points):
         """
         Evaluate the value of the structural frame for the points.
         Can optionally only evaluate one coordinate
@@ -100,11 +100,15 @@ def evaluate_value(self, evaluation_points, i=None):
         -------
 
         """
-        if i is not None:
-            self.features[i].support.evaluate_value(evaluation_points)
-        return (self.features[0].support.evaluate_value(evaluation_points),
-                self.features[1].support.evaluate_value(evaluation_points),
-                self.features[2].support.evaluate_value(evaluation_points))
+        v = np.zeros(evaluation_points.shape) #create new 3d array of correct length
+        v[:] = np.nan
+        v[:,0] = self.features[0].evaluate_value(evaluation_points)
+        v[:,1] = self.features[1].evaluate_value(evaluation_points)
+        v[:,2] = self.features[2].evaluate_value(evaluation_points)
+        return v
+        # return (self.features[0].evaluate_value(evaluation_points),
+        #         self.features[1].evaluate_value(evaluation_points),
+        #         self.features[2].evaluate_value(evaluation_points))
 
     def evaluate_gradient(self, evaluation_points, i=None):
         """

LoopStructural/utils/map2loop.py

@@ -92,7 +92,10 @@ def process_map2loop(m2l_directory, flags={}):
         contact_orientations['gy'] = np.nan
         contact_orientations['gz'] = np.nan
         contact_orientations[['gx', 'gy', 'gz']] = strike_dip_vector(contact_orientations['strike'],
-                                                                    contact_orientations['dip'])*max_thickness 
+                                                                    contact_orientations['dip'])*max_thickness   
+    if np.sum(contact_orientations['polarity']==0) >0 and np.sum(contact_orientations['polarity']==-1)==0:
+        # contact_orientations['polarity']+=1
+        contact_orientations.loc[contact_orientations['polarity']==0]=-1
     if not gradient:
         from LoopStructural.utils.helper import strike_dip_vector
         contact_orientations['strike'] = contact_orientations['azimuth'] - 90
@@ -202,7 +205,9 @@ def process_map2loop(m2l_directory, flags={}):
         fault_centers[3] = np.mean(fault_orientations.loc[fault_orientations['formation']==f,['DipDirection']])
         fault_centers[4] = fault_dimensions.loc[fault_dimensions['Fault']==f,'InfluenceDistance']
         fault_centers[5] = fault_dimensions.loc[fault_dimensions['Fault']==f,'HorizontalRadius']
-        stratigraphic_column['faults'][f] = {'InfluenceDistance':fault_dimensions.loc[fault_dimensions['Fault']==f,'InfluenceDistance'].to_numpy(),
+        stratigraphic_column['faults'][f] = {'FaultCenter':fault_centers[:3],
+                                            'FaultDipDirection':fault_centers[3],
+                                            'InfluenceDistance':fault_dimensions.loc[fault_dimensions['Fault']==f,'InfluenceDistance'].to_numpy(),
                                             'HorizontalRadius':fault_dimensions.loc[fault_dimensions['Fault']==f,'HorizontalRadius'].to_numpy(),
                                             'VerticalRadius':fault_dimensions.loc[fault_dimensions['Fault']==f,'VerticalRadius'].to_numpy()}
         if 'colour' in fault_dimensions.columns: