Speed up fault processing drastically

src/CornerPointGrid/faults.jl

@@ -1,6 +1,12 @@
 function mesh_add_fault_tags!(G::UnstructuredMesh, faults)
     ijk = map(x -> cell_ijk(G, x), 1:number_of_cells(G))
     N = G.faces.neighbors
+
+    sorted_face_pair_I = get_sorted_face_pairs(N, ijk, 1)
+    sorted_face_pair_J = get_sorted_face_pairs(N, ijk, 2)
+    sorted_face_pair_K = get_sorted_face_pairs(N, ijk, 3)
+
+    sorted_faces = (sorted_face_pair_I, sorted_face_pair_J, sorted_face_pair_K)
     for (fault, specs) in faults
         fault_faces = Int[]
         for (I, J, K, dir) in specs
@@ -33,50 +39,47 @@ function mesh_add_fault_tags!(G::UnstructuredMesh, faults)
             range_2 = IJK[ix_2]
             @assert length(range_self) == 1
             self = range_self[1]
-            match_fault_to_faces!(fault_faces, N, ijk, range_1, ix_1, range_2, ix_2, self, ix_self, inc)
+            match_fault_to_faces!(fault_faces, range_1, sorted_faces[ix_1], range_2, sorted_faces[ix_2], self, sorted_faces[ix_self], inc)
         end
         @debug "Fault $fault: Added $(length(fault_faces)) faces"
         set_mesh_entity_tag!(G, Faces(), :faults, Symbol(fault), fault_faces)
     end
 end
 
-function match_fault_to_faces!(fault_faces, N, ijk_cells, range_1, ix_1, range_2, ix_2, self, ix_self, inc)
-    function sorted_tuple(a, b)
-        if a < b
-            pair = (a, b)
-        else
-            pair = (b, a)
-        end
+function get_sorted_face_pairs(N, IJK, ijk_ix)
+    nf = length(N)
+    N_val = similar(N)
+    for i in 1:nf
+        l, r = N[i]
+        N_val[i] = sorted_neighbor_tuple(IJK[l][ijk_ix], IJK[r][ijk_ix])
     end
-    pair = sorted_tuple(self, self+inc)
-    face = 0
-    for (l, r) in N
-        face += 1
-        ijk_l = ijk_cells[l]
-        ijk_r = ijk_cells[r]
-        self_l = ijk_l[ix_self]
-        self_r = ijk_r[ix_self]
-        fpair = sorted_tuple(self_l, self_r)
-        if fpair != pair
-            continue
-        end
-        # Keep going unless fixed indices match
-        cr_1 = ijk_r[ix_1]
-        if !(cr_1 in range_1)
-            continue
-        end
-        cr_2 = ijk_r[ix_2]
-        if !(cr_2 in range_2)
-            continue
-        end
-        cl_1 = ijk_l[ix_1]
-        if !(cl_1 in range_1)
-            continue
-        end
-        cl_2 = ijk_l[ix_2]
-        if !(cl_2 in range_2)
-            continue
+    faces = sortperm(N_val)
+    return (
+        faces_sorted = faces, # Sorted face order
+        N_sorted = N_val[faces], # Sorted tuple pairs
+        N = N_val # Tuple pairs in original order
+    )
+end
+
+function sorted_neighbor_tuple(a, b)
+    if a < b
+        pair = (a, b)
+    else
+        pair = (b, a)
+    end
+end
+
+function match_fault_to_faces!(fault_faces, range_1, sorted_faces_1, range_2, sorted_faces_2, self, sorted_faces_self, inc)
+    pair = sorted_neighbor_tuple(self, self+inc)
+    rng = searchsorted(sorted_faces_self.N_sorted, pair)
+    faces = view(sorted_faces_self.faces_sorted, rng)
+    for face in faces
+        l_1, r_1 = sorted_faces_1.N[face]
+        if l_1 in range_1 && r_1 in range_1
+            l_2, r_2 = sorted_faces_2.N[face]
+            if l_2 in range_2 && r_2 in range_2
+                push!(fault_faces, face)
+            end
         end
-        push!(fault_faces, face)
     end
 end