Add NNC to mesh processing

sintefmath · Apr 18, 2024 · b763b99 · b763b99
1 parent 40afb3a
commit b763b99
Show file tree

Hide file tree

Showing 2 changed files with 29 additions and 2 deletions.
diff --git a/src/CornerPointGrid/interface.jl b/src/CornerPointGrid/interface.jl
@@ -27,16 +27,18 @@ function mesh_from_grid_section(f, actnum = missing)
         actnum = get_effective_actnum(grid)
     end
     cartdims = grid["cartDims"]
+    nnc = get(grid, "NNC", missing)
     if haskey(grid, "COORD")
         coord = grid["COORD"]
         zcorn = grid["ZCORN"]
         primitives = cpgrid_primitives(coord, zcorn, cartdims, actnum = actnum)
-        G = grid_from_primitives(primitives)
+        G = grid_from_primitives(primitives, nnc = nnc)
     else
         @assert haskey(grid, "DX")
         @assert haskey(grid, "DY")
         @assert haskey(grid, "DZ")
         @assert haskey(grid, "TOPS")
+        ismissing(nnc) || throw(ArgumentError("NNC is not supported together with DX/DY/DZ/TOPS mesh."))
         @warn "DX+DY+DZ+TOPS format is only supported if all cells are equally sized and at same TOPS depth. If you get an error, this is the cause."
         @assert all(actnum)
         dx = only(unique(grid["DX"]))

diff --git a/src/CornerPointGrid/processing.jl b/src/CornerPointGrid/processing.jl
@@ -312,7 +312,7 @@ function process_lines!(lines)
     return (nodes, active_lines)
 end
 
-function grid_from_primitives(primitives)
+function grid_from_primitives(primitives; nnc = missing)
     (;
         lines,
         lines_active,
@@ -515,6 +515,31 @@ function grid_from_primitives(primitives)
         end
     end
 
+    if !ismissing(nnc)
+        to_active_ix = zeros(Int, nx*ny*nz)
+        to_active_ix[active] = eachindex(active)
+        function cell_index(i, j, k)
+            ix = ijk_to_linear(i, j, k, cartdims)
+            aix = to_active_ix[ix]
+            return aix
+        end
+
+        for nnc_entry in nnc
+            c1 = cell_index(nnc_entry[1], nnc_entry[2], nnc_entry[3])
+            c2 = cell_index(nnc_entry[4], nnc_entry[5], nnc_entry[6])
+            if c1 > 0 && c2 > 0
+                @assert c1 != c2 "NNC cell pair must be distinct."
+                push!(face_pos, face_pos[end])
+                push!(face_neighbors, (c1, c2))
+                push!(cell_faces[c1], faceno)
+                push!(cell_faces[c2], faceno)
+                faceno += 1
+            else
+                @warn "NNC connects inactive cells, skipped: $(Tuple(nnc_entry[1:3])) -> $(Tuple(nnc_entry[4:6]))"
+            end
+        end
+    end
+
     function convert_to_flat(v)
         flat_vals = Int[]
         flat_pos = Int[1]