Speed up ray tracing

setup.py

@@ -19,7 +19,6 @@
     "pyyaml>=5.2",
     "scikit-learn>=0.22",
     "scipy>=1.4",
-    "tqdm>=4.43",
     "webviz-config>=0.0.42",
     "webviz-config-equinor>=0.0.9",
     "webviz-subsurface>=0.0.24",

src/flownet/network_model/_network_model.py

@@ -1,11 +1,10 @@
 from typing import List, Tuple, Optional, Dict, Any
-from itertools import combinations
+from itertools import combinations, repeat, compress
 
 import numpy as np
 import pandas as pd
 import scipy.spatial
 import pyvista as pv
-from tqdm import tqdm
 
 from ._one_dimensional_model import OneDimensionalModel
 from ..utils.types import Coordinate
@@ -233,13 +232,13 @@ def _calculate_nncs(self) -> List[Tuple[int, int]]:
     # pylint: disable=too-many-locals
     def _calculate_faults(
         self, fault_tolerance: float = 1.0e-05
-    ) -> Optional[Dict[Any, List[int]]]:
+    ) -> Optional[Dict[str, List[int]]]:
         """
         Calculates fault definitions using the following approach:
 
             1) Loop through all faults
-            2) Perform a triangulation of all points belonging to a fault
-            3) For each triangle, perform ray tracing using the
+            2) Perform a triangulation of all points belonging to a fault plane and store the triangles
+            3) For each connection, find all triangles in its bounding box, perform ray tracing using the
                Möller-Trumbore intersection algorithm.
             4) If an intersection is found, identify the grid blocks that are
                associated with the intersection.
@@ -254,17 +253,19 @@ def _calculate_faults(
             Listing of tuples of FAULTS entries with zero-offset i-coordinates, or None if no faults are present.
 
         """
-
-        dict_fault_keyword = {}
+        dict_fault_keyword: Dict[str, List[int]] = {}
         if self._fault_planes is not None:
             fault_names = self._fault_planes["NAME"].unique().tolist()
         if not fault_names:
             return None
 
-        print("Performing fault ray tracing...")
+        print("Performing fault ray tracing...", end=" ", flush=True)
 
-        for fault_name in list(tqdm(fault_names, total=len(fault_names))):
+        # Gather all triangles for all faults and keep track of fault names
+        all_triangles_fault_names: List = []
+        all_triangles = np.empty(shape=[0, 9])
 
+        for fault_name in fault_names:
             data = self._fault_planes.loc[self._fault_planes["NAME"] == fault_name][
                 ["X", "Y", "Z"]
             ].values
@@ -275,30 +276,66 @@ def _calculate_faults(
             vertices = surf.points[surf.faces.reshape(-1, 4)[:, 1:4].ravel()]
 
             triangles = np.array(vertices).reshape(-1, 9)
-            connections = np.hstack(
-                (
-                    np.arange(len(self._df_entity_connections)).reshape(-1, 1),
-                    self._df_entity_connections[
-                        ["xstart", "ystart", "zstart", "xend", "yend", "zend"]
-                    ].values,
-                )
-            )
 
-            combinations_triangles = np.array(triangles).repeat(
-                len(connections), axis=0
+            all_triangles_fault_names.extend(repeat(fault_name, np.shape(triangles)[0]))
+            all_triangles = np.append(all_triangles, triangles, axis=0)
+            dict_fault_keyword[fault_name] = []
+
+        # Loop through all connections and select all triangles inside of the bounding box of the connection
+        # Perform ray tracing on all resulting triangles.
+        for index, row in self._df_entity_connections.iterrows():
+            bx1, bx2 = sorted([row["xstart"], row["xend"]])
+            by1, by2 = sorted([row["ystart"], row["yend"]])
+            bz1, bz2 = sorted([row["zstart"], row["zend"]])
+
+            corner1 = np.array([bx1, by1, bz1])
+            corner2 = np.array([bx2, by2, bz2])
+
+            vertex1_in_box = np.all(
+                np.logical_and(
+                    corner1 <= all_triangles[:, 0:3], all_triangles[:, 0:3] <= corner2
+                ),
+                axis=1,
+            )
+            vertex2_in_box = np.all(
+                np.logical_and(
+                    corner1 <= all_triangles[:, 3:6], all_triangles[:, 3:6] <= corner2
+                ),
+                axis=1,
+            )
+            vertex3_in_box = np.all(
+                np.logical_and(
+                    corner1 <= all_triangles[:, 6:9], all_triangles[:, 6:9] <= corner2
+                ),
+                axis=1,
+            )
+            triangle_in_box = np.any(
+                np.column_stack((vertex1_in_box, vertex2_in_box, vertex3_in_box)),
+                axis=1,
             )
-            connections_connections = np.tile(connections, (len(triangles), 1))
 
-            parameters = list(
-                map(tuple, np.hstack((connections_connections, combinations_triangles)))
+            triangles_in_bounding_box = all_triangles[triangle_in_box]
+            fault_names_in_bounding_box = list(
+                compress(all_triangles_fault_names, triangle_in_box)
             )
 
+            # Loop through all triangles inside of the bounding box and perform ray tracing
             cells_in_fault = []
+            for (triangle, fault_name) in list(
+                zip(triangles_in_bounding_box, fault_names_in_bounding_box)
+            ):
+
+                distance = moller_trumbore(
+                    row["xstart"],
+                    row["ystart"],
+                    row["zstart"],
+                    row["xend"],
+                    row["yend"],
+                    row["zend"],
+                    *triangle
+                )
 
-            for row in list(tqdm(parameters, total=len(parameters), desc=fault_name)):
-                distance, index = moller_trumbore(*row)
-
-                if distance is not False:
+                if distance:
                     indices = self._grid.index[self.active_mask(index)].tolist()
                     tube_cell_index = min(
                         max(0, int(distance * len(indices)) - 1), len(indices) - 2
@@ -308,7 +345,16 @@ def _calculate_faults(
                     cells_in_fault.append(cell_i_index)
 
             if len(cells_in_fault) > 0:
-                dict_fault_keyword[fault_name] = list(set(cells_in_fault))
+                dict_fault_keyword[fault_name].extend(cells_in_fault)
+
+        # Remove empty and double entries
+        for fault_name in fault_names:
+            if not dict_fault_keyword[fault_name]:
+                dict_fault_keyword.pop(fault_name)
+            else:
+                dict_fault_keyword[fault_name] = list(
+                    set(dict_fault_keyword[fault_name])
+                )
 
         print("done.")
 

src/flownet/utils/raytracing.py

@@ -1,11 +1,10 @@
-from typing import Tuple, Any
+from typing import Optional
 
 import numpy as np
 
 
 # pylint: disable=too-many-arguments,too-many-locals
 def moller_trumbore(
-    index: int,
     xstart: float,
     ystart: float,
     zstart: float,
@@ -21,7 +20,7 @@ def moller_trumbore(
     v31: float,
     v32: float,
     v33: float,
-) -> Tuple[Any, Any]:
+) -> Optional[float]:
     """
     The Möller–Trumbore ray-triangle intersection algorithm, named after its inventors Tomas Möller and Ben Trumbore,
     is a fast method for calculating the intersection of a ray and a triangle in three dimensions without
@@ -31,7 +30,6 @@ def moller_trumbore(
     (a collection of triangles).
 
     Args:
-        index: Tube index
         xstart: X-coordinate of the start point of the tube
         ystart: Y-coordinate of the start point of the tube
         zstart: Z-coordinate of the start point of the tube
@@ -49,8 +47,7 @@ def moller_trumbore(
         v33: Triangle corner 3 Z-coordinate
 
     Returns:
-        Tuple of the float representing the fraction of the tube where the intersection happened as well as
-        the tube index it happened in, otherwise False, False
+        Either a float representing the fraction of the tube where the intersection happened or None
 
     """
     eps = 0.000001
@@ -64,27 +61,27 @@ def moller_trumbore(
     det = edge1.dot(pvec)
 
     if abs(det) < eps:
-        return False, False
+        return None
 
     inv_det = 1.0 / det
     t_vector = ray_start - triangle[0]
     u_value = t_vector.dot(pvec) * inv_det
 
     if u_value < 0.0 or u_value > 1.0:
-        return False, False
+        return None
 
     q_vector = np.cross(t_vector, edge1)
     v_value = ray_direction.dot(q_vector) * inv_det
 
     if v_value < 0.0 or u_value + v_value > 1.0:
-        return False, False
+        return None
 
     ray_fraction = edge2.dot(q_vector) * inv_det
 
     if ray_fraction < eps:
-        return False, False
+        return None
 
     if ray_fraction > 1:
-        return False, False
+        return None
 
-    return ray_fraction, index
+    return ray_fraction