diff --git a/setup.py b/setup.py index 95b5f1f55..19ce603f9 100644 --- a/setup.py +++ b/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", diff --git a/src/flownet/network_model/_network_model.py b/src/flownet/network_model/_network_model.py index c288f8d19..8a445d043 100644 --- a/src/flownet/network_model/_network_model.py +++ b/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.") diff --git a/src/flownet/utils/raytracing.py b/src/flownet/utils/raytracing.py index 0761cf2b2..ed816d6fb 100644 --- a/src/flownet/utils/raytracing.py +++ b/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