fix: use fault normals/slip vectors from data if available.

Otherwise calculate from trace

Author: Lachlan Grose

LoopStructural/modelling/features/builders/_fault_builder.py

@@ -151,16 +151,30 @@ def create_data_from_geometry(
             ["X", "Y"],
         ].to_numpy()
         if fault_normal_vector is None:
-            # Calculate fault strike using least squares fit
-            X = fault_trace[:, 0].reshape(-1, 1)
-            Y = fault_trace[:, 1]
-            # Fit line Y = mX + b
-            A = np.hstack([X, np.ones_like(X)])
-            m, _ = np.linalg.lstsq(A, Y, rcond=None)[0]
-            # Convert slope to strike vector
-            strike_vector = np.array([1, m, 0])
-            strike_vector /= np.linalg.norm(strike_vector)
-            fault_normal_vector = np.cross(strike_vector, [0, 0, 1])
+            if fault_frame_data.loc[
+            np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["nx"].notna())].shape[0]>0:
+                fault_normal_vector = fault_frame_data.loc[
+                    np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["nx"].notna()),
+                    ["nx", "ny", "nz"],
+                ].to_numpy().mean(axis=0)
+
+            else:
+
+                # Calculate fault strike using eigenvectors
+                pts = fault_trace - fault_trace.mean(axis=0)
+                # Calculate covariance matrix
+                cov_matrix = pts.T @ pts
+                # Get eigenvectors and eigenvalues
+                eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
+                # Use eigenvector with largest eigenvalue as strike direction
+                strike_vector = eigenvectors[:, np.argmax(eigenvalues)]
+                strike_vector = np.append(strike_vector, 0)  # Add z component
+                strike_vector /= np.linalg.norm(strike_vector)
+
+                fault_normal_vector = np.cross(strike_vector, [0, 0, 1])
+                # Rotate the fault normal vector according to the fault dip
+                rotation_matrix = rotation(strike_vector[None, :], np.array([90 - fault_dip]))
+                fault_normal_vector = np.einsum("ijk,ik->ij", rotation_matrix, fault_normal_vector[None, :])[0]
 
         if not isinstance(fault_normal_vector, np.ndarray):
             fault_normal_vector = np.array(fault_normal_vector)
@@ -170,10 +184,18 @@ def create_data_from_geometry(
             fault_slip_vector = np.einsum("ijk,ik->ij", rotation_matrix, fault_normal_vector[None, :])[0]
 
         if fault_slip_vector is None:
-            fault_slip_vector = np.cross(fault_normal_vector, [1., 0., 0.])
-            if np.linalg.norm(fault_slip_vector) == 0:
-                fault_slip_vector = np.cross(fault_normal_vector, [0., 1., 0.])
-            fault_slip_vector /= np.linalg.norm(fault_slip_vector)
+            if fault_frame_data.loc[
+            np.logical_and(fault_frame_data["coord"] == 1, fault_frame_data["nx"].notna())].shape[0]>0:
+                fault_slip_vector = fault_frame_data.loc[
+                    np.logical_and(fault_frame_data["coord"] == 1, fault_frame_data["nx"].notna()),
+                    ["nx", "ny", "nz"],
+                ].to_numpy().mean(axis=0)
+
+            else:
+                fault_slip_vector = np.cross(fault_normal_vector, [1., 0., 0.])
+                if np.linalg.norm(fault_slip_vector) == 0:
+                    fault_slip_vector = np.cross(fault_normal_vector, [0., 1., 0.])
+                fault_slip_vector /= np.linalg.norm(fault_slip_vector)
         if fault_center is None:
             fault_trace = fault_frame_data.loc[
                 np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0),

tests/unit/modelling/test__fault_builder.py

@@ -21,13 +21,16 @@ def test_fault_builder_update_geometry(interpolatortype):
     fault_builder.update_geometry(points)
 
     # Check if the origin and maximum are updated correctly
-    expected_origin = np.array([0.5, 0.5, 0.5])
-    expected_maximum = np.array([0.7, 0.7, 0.7])
-    assert np.allclose(fault_builder.origin, expected_origin)
-    assert np.allclose(fault_builder.maximum, expected_maximum)
+    # Calculate buffer as 10% of the range
+    buffer = 0.1 * (np.array([0.7, 0.7, 0.7]) - np.array([0.5, 0.5, 0.5]))
+    expected_origin = np.array([0.5, 0.5, 0.5]) - buffer
+    expected_maximum = np.array([0.7, 0.7, 0.7]) + buffer
+    # rtol 1e-2 is 1% relative tolerance 
+    assert np.allclose(fault_builder.origin, expected_origin,rtol=1e-2)
+    assert np.allclose(fault_builder.maximum, expected_maximum, rtol=1e-2)
 
 
-def test_fault_builder_create_data_from_geometry(interpolatortype):
+def test_fault_builder_create_data_from_geometry_with_normals(interpolatortype):
     # Create a FaultBuilder instance
     bounding_box = BoundingBox([0, 0, 0], [1, 1, 1])
     nelements = 1000
@@ -60,7 +63,59 @@ def test_fault_builder_create_data_from_geometry(interpolatortype):
     # Check if the fault attributes are updated correctly
     expected_fault_normal_vector = np.array([0, 1.0, 0])
     expected_fault_slip_vector = np.array([0.0, 0.0, 1.0])
-    # expected_fault_centre = np.array([0.15, 0.6, 1.05])
+    expected_fault_centre = np.array([0.15, 0.6, 1.05])
+    expected_fault_normal_vector /= np.linalg.norm(expected_fault_normal_vector)
+    expected_fault_slip_vector /= np.linalg.norm(expected_fault_slip_vector)
+    # expected_fault_minor_axis = 0.5
+    # expected_fault_major_axis = 1.0
+    # expected_fault_intermediate_axis = 1.0
+    calculated_fault_normal_vector = fault_builder.fault_normal_vector
+    calculated_fault_slip_vector = fault_builder.fault_slip_vector
+    calculated_fault_normal_vector /= np.linalg.norm(calculated_fault_normal_vector)
+    calculated_fault_slip_vector /= np.linalg.norm(calculated_fault_slip_vector)
+    assert np.allclose(calculated_fault_normal_vector, expected_fault_normal_vector)
+    assert np.allclose(calculated_fault_slip_vector, expected_fault_slip_vector)
+    # assert np.allclose(fault_builder.fault_centre, expected_fault_centre)
+    # assert np.isclose(fault_builder.fault_minor_axis, expected_fault_minor_axis)
+    # assert np.isclose(fault_builder.fault_major_axis, expected_fault_major_axis)
+    # assert np.isclose(fault_builder.fault_intermediate_axis, expected_fault_intermediate_axis)
+
+
+def test_fault_builder_create_data_from_geometry_without_normals(interpolatortype):
+    # Create a FaultBuilder instance
+    bounding_box = BoundingBox([0, 0, 0], [1, 1, 1])
+    nelements = 1000
+    model = GeologicalModel([0, 0, 0], [1, 1, 1])
+    fault_bounding_box_buffer = 0.2
+    fault_builder = FaultBuilder(
+        interpolatortype, bounding_box, nelements, model, fault_bounding_box_buffer
+    )
+
+    # Create some test data
+    fault_frame_data = pd.DataFrame(
+        {
+            "coord": [0.0, 0.0, 1.0, 2.0],
+            "val": [0.0, 0.0, 1.0, 1.0],
+            "gx": [np.nan, np.nan, np.nan, np.nan],
+            "gy": [np.nan, np.nan, np.nan, np.nan],
+            "gz": [np.nan, np.nan, np.nan, np.nan],
+            "X": [0.1, 0.2, 0.3, 0.4],
+            "Y": [0.5, 0.6, 0.7, 0.8],
+            "Z": [0.9, 1.0, 1.1, 1.2],
+            "nx": [np.nan, np.nan, np.nan, np.nan],
+            "ny": [np.nan, np.nan, np.nan, np.nan],
+            "nz": [np.nan, np.nan, np.nan, np.nan],
+        }
+    )
+
+    # Call the create_data_from_geometry method
+    fault_builder.create_data_from_geometry(fault_frame_data)
+    strike_vector = np.array([-.1,-.1,0])
+    strike_vector /= np.linalg.norm(strike_vector)
+    expected_fault_normal_vector = np.cross(strike_vector, np.array([0, 0, 1]))
+    expected_fault_normal_vector /= np.linalg.norm(expected_fault_normal_vector)
+    # Check if the fault attributes are updated correctly
+    expected_fault_slip_vector = np.cross(expected_fault_normal_vector, [1.0, 0.0, 0.0])
     expected_fault_normal_vector /= np.linalg.norm(expected_fault_normal_vector)
     expected_fault_slip_vector /= np.linalg.norm(expected_fault_slip_vector)
     # expected_fault_minor_axis = 0.5