From dab7339622145c51762ae31e68a0edf98b37bfda Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 13 Feb 2026 13:06:41 +0000
Subject: [PATCH 01/14] Initial plan


From 3c9814aee95b326ffb5d7f6e09f229593db05599 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 13 Feb 2026 14:13:25 +0000
Subject: [PATCH 02/14] Add partial_cut mode to cut_mesh and layered_mesh
 functionality

- Added partial_cut keyword to cut_mesh (:positive/:negative/:none)
  for keeping only one side of the cut plane
- Added depth_grid_to_surface to convert depth grids to PolygonalSurface
- Added layered_mesh to build layered reservoir models from surfaces
- Added _point_below_surface helper for classifying cell layers
- Updated module exports and imports

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>
---
 src/Jutul.jl                              |   2 +-
 src/meshes/CutCellMeshes/CutCellMeshes.jl |   2 +
 src/meshes/CutCellMeshes/cutting.jl       | 233 +++++++++++++++++-----
 src/meshes/CutCellMeshes/layered.jl       | 184 +++++++++++++++++
 4 files changed, 371 insertions(+), 50 deletions(-)
 create mode 100644 src/meshes/CutCellMeshes/layered.jl

diff --git a/src/Jutul.jl b/src/Jutul.jl
index cde6a4c02..5d0dfc0f4 100644
--- a/src/Jutul.jl
+++ b/src/Jutul.jl
@@ -156,7 +156,7 @@ module Jutul
 
     # Cut-cell meshes
     include("meshes/CutCellMeshes/CutCellMeshes.jl")
-    import Jutul.CutCellMeshes: cut_mesh, PlaneCut, PolygonalSurface, glue_mesh, cut_and_displace_mesh
+    import Jutul.CutCellMeshes: cut_mesh, PlaneCut, PolygonalSurface, glue_mesh, cut_and_displace_mesh, layered_mesh, depth_grid_to_surface
 
     # This is to make Jutul simulators work nicely with nested ForwardDiff.
     JutulSimulateTag = ForwardDiff.Tag{typeof(simulate), <:JutulEntity}
diff --git a/src/meshes/CutCellMeshes/CutCellMeshes.jl b/src/meshes/CutCellMeshes/CutCellMeshes.jl
index d2be86468..7f6e2ba12 100644
--- a/src/meshes/CutCellMeshes/CutCellMeshes.jl
+++ b/src/meshes/CutCellMeshes/CutCellMeshes.jl
@@ -7,6 +7,8 @@ module CutCellMeshes
     include("geometry.jl")
     include("cutting.jl")
     include("gluing.jl")
+    include("layered.jl")
 
     export cut_mesh, glue_mesh, cut_and_displace_mesh
+    export layered_mesh, depth_grid_to_surface
 end
diff --git a/src/meshes/CutCellMeshes/cutting.jl b/src/meshes/CutCellMeshes/cutting.jl
index 93ca9d3e3..83697113e 100644
--- a/src/meshes/CutCellMeshes/cutting.jl
+++ b/src/meshes/CutCellMeshes/cutting.jl
@@ -76,12 +76,23 @@ function cut_mesh(mesh::UnstructuredMesh{3}, surface::PolygonalSurface; extra_ou
     end
 end
 
+"""
+    cut_mesh(mesh, plane; partial_cut=:none, ...)
+
+When `partial_cut` is `:positive` or `:negative`, only cells on the specified
+side of the cutting plane are kept after cutting.  Cut cells lose their
+sub-cell on the discarded side and the cut face becomes a new boundary face.
+Uncut cells that lie entirely on the discarded side are removed.  When
+`partial_cut=:none` (default), both sub-cells are kept (standard behaviour).
+"""
 function cut_mesh(mesh::UnstructuredMesh{3}, plane::PlaneCut{T};
         min_cut_fraction::Real = 0.05,
         bounding_polygon = nothing,
         clip_to_polygon::Bool = false,
-        extra_out::Bool = false
+        extra_out::Bool = false,
+        partial_cut::Symbol = :none
     ) where T
+    @assert partial_cut in (:none, :positive, :negative) "partial_cut must be :none, :positive, or :negative"
     nc = number_of_cells(mesh)
     nn = length(mesh.node_points)
 
@@ -305,7 +316,7 @@ function cut_mesh(mesh::UnstructuredMesh{3}, plane::PlaneCut{T};
         end
     end
 
-    return build_cut_mesh(mesh, plane, new_node_points, is_cut, cut_cell_infos, node_class, get_or_create_intersection, extra_out)
+    return build_cut_mesh(mesh, plane, new_node_points, is_cut, cut_cell_infos, node_class, get_or_create_intersection, extra_out, partial_cut)
 end
 
 """
@@ -334,9 +345,14 @@ function dominant_side(fnodes, node_class)
 end
 
 """
-    build_cut_mesh(mesh, plane, node_points, is_cut, cut_infos, node_class, get_intersection, extra_out)
+    build_cut_mesh(mesh, plane, node_points, is_cut, cut_infos, node_class, get_intersection, extra_out, partial_cut)
 
 Build the final mesh after cutting.
+
+When `partial_cut` is `:positive` or `:negative`, only cells on the specified
+side are kept.  Cut cells lose the sub-cell on the discarded side and the cut
+face becomes a boundary face.  Uncut cells entirely on the discarded side are
+removed.
 """
 function build_cut_mesh(
     mesh::UnstructuredMesh{3},
@@ -346,12 +362,37 @@ function build_cut_mesh(
     cut_infos::Dict{Int, CutCellInfo},
     node_class::Dict{Int, Int},
     get_intersection::Function,
-    extra_out::Bool
+    extra_out::Bool,
+    partial_cut::Symbol = :none
 ) where T
     nc_old = number_of_cells(mesh)
     nf_old = number_of_faces(mesh)
     nb_old = number_of_boundary_faces(mesh)
 
+    # ------------------------------------------------------------------
+    # For partial_cut, classify uncut cells and decide which to discard
+    # ------------------------------------------------------------------
+    # cell_side[c]: :positive, :negative, or :cut
+    cell_side = Vector{Symbol}(undef, nc_old)
+    for c in 1:nc_old
+        if is_cut[c]
+            cell_side[c] = :cut
+        else
+            cell_side[c] = classify_cell(mesh, c, plane)
+        end
+    end
+
+    # keep_cell[c] - whether this old cell produces at least one new cell
+    keep_cell = trues(nc_old)
+    if partial_cut != :none
+        discard = partial_cut == :positive ? :negative : :positive
+        for c in 1:nc_old
+            if !is_cut[c] && cell_side[c] == discard
+                keep_cell[c] = false
+            end
+        end
+    end
+
     # Compute new cell numbering
     new_cell_count = 0
     old_to_new = Dict{Int, Vector{Int}}()
@@ -360,13 +401,27 @@ function build_cut_mesh(
     cell_index = Int[]
 
     for c in 1:nc_old
+        if !keep_cell[c]
+            old_to_new[c] = Int[]
+            continue
+        end
         if is_cut[c]
-            pos_cell = new_cell_count + 1
-            neg_cell = new_cell_count + 2
-            old_to_new[c] = [pos_cell, neg_cell]
-            new_cell_count += 2
-            push!(cell_index, c)  # pos sub-cell maps to original cell c
-            push!(cell_index, c)  # neg sub-cell maps to original cell c
+            if partial_cut == :none
+                pos_cell = new_cell_count + 1
+                neg_cell = new_cell_count + 2
+                old_to_new[c] = [pos_cell, neg_cell]
+                new_cell_count += 2
+                push!(cell_index, c)  # pos sub-cell
+                push!(cell_index, c)  # neg sub-cell
+            elseif partial_cut == :positive
+                new_cell_count += 1
+                old_to_new[c] = [new_cell_count, 0]  # [pos, no neg]
+                push!(cell_index, c)
+            else  # :negative
+                new_cell_count += 1
+                old_to_new[c] = [0, new_cell_count]  # [no pos, neg]
+                push!(cell_index, c)
+            end
         else
             new_cell_count += 1
             old_to_new[c] = [new_cell_count]
@@ -452,6 +507,22 @@ function build_cut_mesh(
         return has_pos && has_neg
     end
 
+    # Helper: get the new cell index for a given old cell and side
+    # For cut cells in partial_cut mode, old_to_new has [pos, neg] where one is 0
+    # For cut cells in normal mode, old_to_new has [pos, neg]
+    # For uncut cells, old_to_new has [cell]
+    function get_new_cell(old_cell::Int, side::Symbol)
+        mapping = old_to_new[old_cell]
+        if length(mapping) == 0
+            return 0  # cell was removed
+        elseif length(mapping) == 1
+            return mapping[1]  # uncut cell - single new cell
+        else
+            # Cut cell: mapping = [pos, neg]
+            return side == :positive ? mapping[1] : mapping[2]
+        end
+    end
+
     # Process interior faces
     for face in 1:nf_old
         l_old, r_old = mesh.faces.neighbors[face]
@@ -459,75 +530,102 @@ function build_cut_mesh(
         l_cut = is_cut[l_old]
         r_cut = is_cut[r_old]
 
+        if !keep_cell[l_old] && !keep_cell[r_old]
+            continue  # both cells removed
+        end
+
         if !l_cut && !r_cut
-            l_new = old_to_new[l_old][1]
-            r_new = old_to_new[r_old][1]
-            add_interior_face!(fnodes, l_new, r_new; old_face = face)
+            l_new = get_new_cell(l_old, cell_side[l_old])
+            r_new = get_new_cell(r_old, cell_side[r_old])
+            if l_new == 0 && r_new == 0
+                continue
+            elseif l_new == 0
+                # Left cell removed → face becomes boundary on right cell
+                add_boundary_face!(fnodes, r_new; old_bf = 0)
+            elseif r_new == 0
+                # Right cell removed → face becomes boundary on left cell
+                add_boundary_face!(fnodes, l_new; old_bf = 0)
+            else
+                add_interior_face!(fnodes, l_new, r_new; old_face = face)
+            end
         elseif l_cut && !r_cut
-            r_new = old_to_new[r_old][1]
-            l_pos = old_to_new[l_old][1]
-            l_neg = old_to_new[l_old][2]
+            r_new = get_new_cell(r_old, cell_side[r_old])
+            l_pos = get_new_cell(l_old, :positive)
+            l_neg = get_new_cell(l_old, :negative)
 
             if face_needs_split(fnodes)
                 pos_fn, neg_fn = split_face_cached(fnodes)
                 if length(pos_fn) >= 3
-                    add_interior_face!(pos_fn, l_pos, r_new; old_face = face)
+                    _add_face_or_bnd!(pos_fn, l_pos, r_new, face,
+                        add_interior_face!, add_boundary_face!)
                 end
                 if length(neg_fn) >= 3
-                    add_interior_face!(neg_fn, l_neg, r_new; old_face = face)
+                    _add_face_or_bnd!(neg_fn, l_neg, r_new, face,
+                        add_interior_face!, add_boundary_face!)
                 end
             else
                 side = dominant_side(fnodes, node_class)
                 if side >= 0
-                    add_interior_face!(fnodes, l_pos, r_new; old_face = face)
+                    _add_face_or_bnd!(fnodes, l_pos, r_new, face,
+                        add_interior_face!, add_boundary_face!)
                 else
-                    add_interior_face!(fnodes, l_neg, r_new; old_face = face)
+                    _add_face_or_bnd!(fnodes, l_neg, r_new, face,
+                        add_interior_face!, add_boundary_face!)
                 end
             end
         elseif !l_cut && r_cut
-            l_new = old_to_new[l_old][1]
-            r_pos = old_to_new[r_old][1]
-            r_neg = old_to_new[r_old][2]
+            l_new = get_new_cell(l_old, cell_side[l_old])
+            r_pos = get_new_cell(r_old, :positive)
+            r_neg = get_new_cell(r_old, :negative)
 
             if face_needs_split(fnodes)
                 pos_fn, neg_fn = split_face_cached(fnodes)
                 if length(pos_fn) >= 3
-                    add_interior_face!(pos_fn, l_new, r_pos; old_face = face)
+                    _add_face_or_bnd!(pos_fn, l_new, r_pos, face,
+                        add_interior_face!, add_boundary_face!)
                 end
                 if length(neg_fn) >= 3
-                    add_interior_face!(neg_fn, l_new, r_neg; old_face = face)
+                    _add_face_or_bnd!(neg_fn, l_new, r_neg, face,
+                        add_interior_face!, add_boundary_face!)
                 end
             else
                 side = dominant_side(fnodes, node_class)
                 if side >= 0
-                    add_interior_face!(fnodes, l_new, r_pos; old_face = face)
+                    _add_face_or_bnd!(fnodes, l_new, r_pos, face,
+                        add_interior_face!, add_boundary_face!)
                 else
-                    add_interior_face!(fnodes, l_new, r_neg; old_face = face)
+                    _add_face_or_bnd!(fnodes, l_new, r_neg, face,
+                        add_interior_face!, add_boundary_face!)
                 end
             end
         else
             # Both sides cut
-            l_pos = old_to_new[l_old][1]
-            l_neg = old_to_new[l_old][2]
-            r_pos = old_to_new[r_old][1]
-            r_neg = old_to_new[r_old][2]
+            l_pos = get_new_cell(l_old, :positive)
+            l_neg = get_new_cell(l_old, :negative)
+            r_pos = get_new_cell(r_old, :positive)
+            r_neg = get_new_cell(r_old, :negative)
 
             if face_needs_split(fnodes)
                 pos_fn, neg_fn = split_face_cached(fnodes)
                 if length(pos_fn) >= 3
-                    add_interior_face!(pos_fn, l_pos, r_pos; old_face = face)
+                    _add_face_or_bnd!(pos_fn, l_pos, r_pos, face,
+                        add_interior_face!, add_boundary_face!)
                 end
                 if length(neg_fn) >= 3
-                    add_interior_face!(neg_fn, l_neg, r_neg; old_face = face)
+                    _add_face_or_bnd!(neg_fn, l_neg, r_neg, face,
+                        add_interior_face!, add_boundary_face!)
                 end
             else
                 side = dominant_side(fnodes, node_class)
                 if side >= 0
-                    add_interior_face!(fnodes, l_pos, r_pos; old_face = face)
+                    _add_face_or_bnd!(fnodes, l_pos, r_pos, face,
+                        add_interior_face!, add_boundary_face!)
                 elseif side < 0
-                    add_interior_face!(fnodes, l_neg, r_neg; old_face = face)
+                    _add_face_or_bnd!(fnodes, l_neg, r_neg, face,
+                        add_interior_face!, add_boundary_face!)
                 else
-                    add_interior_face!(fnodes, l_pos, r_pos; old_face = face)
+                    _add_face_or_bnd!(fnodes, l_pos, r_pos, face,
+                        add_interior_face!, add_boundary_face!)
                 end
             end
         end
@@ -538,43 +636,63 @@ function build_cut_mesh(
         cell_old = mesh.boundary_faces.neighbors[bf]
         fnodes = collect(mesh.boundary_faces.faces_to_nodes[bf])
 
+        if !keep_cell[cell_old]
+            continue  # cell removed
+        end
+
         if !is_cut[cell_old]
-            cell_new = old_to_new[cell_old][1]
-            add_boundary_face!(fnodes, cell_new; old_bf = bf)
+            cell_new = get_new_cell(cell_old, cell_side[cell_old])
+            if cell_new != 0
+                add_boundary_face!(fnodes, cell_new; old_bf = bf)
+            end
         else
-            pos_cell = old_to_new[cell_old][1]
-            neg_cell = old_to_new[cell_old][2]
+            pos_cell = get_new_cell(cell_old, :positive)
+            neg_cell = get_new_cell(cell_old, :negative)
 
             if face_needs_split(fnodes)
                 pos_fn, neg_fn = split_face_cached(fnodes)
-                if length(pos_fn) >= 3
+                if length(pos_fn) >= 3 && pos_cell != 0
                     add_boundary_face!(pos_fn, pos_cell; old_bf = bf)
                 end
-                if length(neg_fn) >= 3
+                if length(neg_fn) >= 3 && neg_cell != 0
                     add_boundary_face!(neg_fn, neg_cell; old_bf = bf)
                 end
             else
                 side = dominant_side(fnodes, node_class)
-                if side >= 0
+                if side >= 0 && pos_cell != 0
                     add_boundary_face!(fnodes, pos_cell; old_bf = bf)
-                else
+                elseif side < 0 && neg_cell != 0
                     add_boundary_face!(fnodes, neg_cell; old_bf = bf)
+                elseif side >= 0 && neg_cell != 0
+                    # Face is on positive side but only negative cell kept - skip
+                elseif side < 0 && pos_cell != 0
+                    # Face is on negative side but only positive cell kept - skip
                 end
             end
         end
     end
 
-    # Add cut faces (internal faces between pos and neg sub-cells)
+    # Add cut faces
     for c in 1:nc_old
         if !is_cut[c]
             continue
         end
         info = cut_infos[c]
         if length(info.cut_face_nodes) >= 3
-            pos_cell = old_to_new[c][1]
-            neg_cell = old_to_new[c][2]
-            fi = add_interior_face!(info.cut_face_nodes, neg_cell, pos_cell; old_face = 0)
-            push!(new_faces_list, fi)
+            pos_cell = get_new_cell(c, :positive)
+            neg_cell = get_new_cell(c, :negative)
+            if partial_cut == :none
+                # Both sub-cells exist: cut face is interior between them
+                fi = add_interior_face!(info.cut_face_nodes, neg_cell, pos_cell; old_face = 0)
+                push!(new_faces_list, fi)
+            else
+                # Only one sub-cell: cut face becomes boundary
+                kept = partial_cut == :positive ? pos_cell : neg_cell
+                if kept != 0
+                    bi = add_boundary_face!(info.cut_face_nodes, kept; old_bf = 0)
+                    push!(new_faces_list, bi)
+                end
+            end
         end
     end
 
@@ -632,3 +750,20 @@ function build_cut_mesh(
     end
     return new_mesh
 end
+
+"""
+    _add_face_or_bnd!(nodes, left, right, old_face, add_int!, add_bnd!)
+
+Helper: if both `left` and `right` are valid (non-zero) cell indices, add an
+interior face; if only one is valid, add a boundary face on that cell; if
+neither is valid, skip.
+"""
+function _add_face_or_bnd!(nodes, left, right, old_face, add_int!, add_bnd!)
+    if left != 0 && right != 0
+        add_int!(nodes, left, right; old_face = old_face)
+    elseif left != 0
+        add_bnd!(nodes, left; old_bf = 0)
+    elseif right != 0
+        add_bnd!(nodes, right; old_bf = 0)
+    end
+end
diff --git a/src/meshes/CutCellMeshes/layered.jl b/src/meshes/CutCellMeshes/layered.jl
new file mode 100644
index 000000000..07c858ce2
--- /dev/null
+++ b/src/meshes/CutCellMeshes/layered.jl
@@ -0,0 +1,184 @@
+"""
+    depth_grid_to_surface(xs, ys, depths)
+
+Convert a regular depth grid to a triangulated [`PolygonalSurface`](@ref).
+
+`xs` and `ys` are vectors of coordinates defining the grid axes.  `depths` is a
+matrix of size `(length(xs), length(ys))` giving the depth (z-coordinate) at
+each grid node.  The grid is triangulated into pairs of triangles per
+rectangular cell.
+
+Entries that are `NaN` or `missing` are skipped, so incomplete horizons are
+supported.
+
+The resulting surface normal points upward (positive z) for each triangle.
+"""
+function depth_grid_to_surface(
+    xs::AbstractVector{<:Real},
+    ys::AbstractVector{<:Real},
+    depths::AbstractMatrix{<:Real}
+)
+    T = Float64
+    nx = length(xs)
+    ny = length(ys)
+    @assert size(depths) == (nx, ny) "depths must be (length(xs), length(ys))"
+
+    polygons = Vector{SVector{3, T}}[]
+
+    for i in 1:(nx - 1)
+        for j in 1:(ny - 1)
+            z00 = depths[i,     j    ]
+            z10 = depths[i + 1, j    ]
+            z01 = depths[i,     j + 1]
+            z11 = depths[i + 1, j + 1]
+
+            p00 = SVector{3, T}(xs[i],     ys[j],     z00)
+            p10 = SVector{3, T}(xs[i + 1], ys[j],     z10)
+            p01 = SVector{3, T}(xs[i],     ys[j + 1], z01)
+            p11 = SVector{3, T}(xs[i + 1], ys[j + 1], z11)
+
+            # Lower-left triangle: p00-p10-p01
+            if _depth_valid(z00) && _depth_valid(z10) && _depth_valid(z01)
+                push!(polygons, [p00, p10, p01])
+            end
+            # Upper-right triangle: p10-p11-p01
+            if _depth_valid(z10) && _depth_valid(z11) && _depth_valid(z01)
+                push!(polygons, [p10, p11, p01])
+            end
+        end
+    end
+
+    return PolygonalSurface(polygons)
+end
+
+_depth_valid(z::Real) = isfinite(z) && !isnan(z)
+_depth_valid(::Missing) = false
+
+"""
+    layered_mesh(mesh, surfaces; min_cut_fraction=0.01)
+
+Build a layered reservoir mesh from a base mesh and a set of depth surfaces.
+
+`mesh` is a 3D `UnstructuredMesh`.  `surfaces` is a vector of
+[`PolygonalSurface`](@ref) objects ordered by increasing depth (shallowest
+first).  Each surface defines a horizon that separates two layers.
+
+The function returns `(result_mesh, info)` where:
+- `result_mesh` is an `UnstructuredMesh` with the same total volume as the
+  original mesh (all cells are preserved, just split at surfaces).
+- `info` is a `Dict{String, Any}` containing:
+  - `"layer_indices"`: A `Vector{Int}` of length `number_of_cells(result_mesh)`.
+    Each entry gives the layer number of the corresponding cell.  Layer `0` is
+    above all surfaces, layer `k` (1 ≤ k ≤ N-1 where N = length(surfaces))
+    is between surface `k` and surface `k+1`, and layer `N` is below all
+    surfaces.
+  - `"cell_index"`: A `Vector{Int}` mapping each new cell to the original
+    cell index in the input mesh.
+
+# Algorithm
+
+Each surface is applied sequentially from shallowest to deepest using the
+standard `cut_mesh` (keeping both halves of every cut cell).  After all cuts,
+each cell's centroid is compared against every surface to determine its layer.
+
+If a shallower surface dips below a deeper one at some location, the shallower
+surface takes priority (the cell is classified by the topmost surface it is
+below).
+"""
+function layered_mesh(
+    mesh::UnstructuredMesh{3},
+    surfaces::Vector{<:PolygonalSurface};
+    min_cut_fraction::Real = 0.01
+)
+    T = Float64
+    n_surfaces = length(surfaces)
+    nc_orig = number_of_cells(mesh)
+
+    # ----------------------------------------------------------------
+    # 1. Apply all surface cuts sequentially
+    # ----------------------------------------------------------------
+    current_mesh = mesh
+    cell_map = collect(1:nc_orig)  # maps new cell → original cell
+
+    for (si, surface) in enumerate(surfaces)
+        for (pi, poly) in enumerate(surface.polygons)
+            n = surface.normals[pi]
+            c = sum(poly) / length(poly)
+            plane = PlaneCut(c, n)
+
+            result, step_info = cut_mesh(current_mesh, plane;
+                extra_out = true,
+                min_cut_fraction = min_cut_fraction
+            )
+
+            # Compose cell maps
+            cell_map = [cell_map[j] for j in step_info["cell_index"]]
+            current_mesh = result
+        end
+    end
+
+    # ----------------------------------------------------------------
+    # 2. Classify each cell by layer using centroid position
+    # ----------------------------------------------------------------
+    nc_final = number_of_cells(current_mesh)
+    geo = tpfv_geometry(current_mesh)
+    layer_indices = Vector{Int}(undef, nc_final)
+
+    for c in 1:nc_final
+        cx = geo.cell_centroids[1, c]
+        cy = geo.cell_centroids[2, c]
+        cz = geo.cell_centroids[3, c]
+        centroid = SVector{3, T}(cx, cy, cz)
+
+        # Find the first surface this cell is above (positive side)
+        # layer = number of surfaces that the centroid is below
+        layer = 0
+        for (si, surface) in enumerate(surfaces)
+            if _point_below_surface(centroid, surface)
+                layer = si
+            else
+                break
+            end
+        end
+        layer_indices[c] = layer
+    end
+
+    info = Dict{String, Any}(
+        "layer_indices" => layer_indices,
+        "cell_index" => cell_map
+    )
+
+    return (current_mesh, info)
+end
+
+"""
+    _point_below_surface(pt, surface)
+
+Check whether a 3D point is deeper than (below) a polygonal surface.
+
+For each polygon the signed distance from the polygon's plane is computed.  The
+polygon whose centroid is closest in the horizontal (xy) plane is used.  A
+positive signed distance from that plane (in the direction of the polygon
+normal) corresponds to *greater depth* when the normal is oriented in the
+depth-increasing direction — which is the default produced by
+[`depth_grid_to_surface`](@ref).
+"""
+function _point_below_surface(pt::SVector{3, T}, surface::PolygonalSurface) where T
+    best_dist = T(Inf)
+    best_sd = zero(T)
+
+    for (i, poly) in enumerate(surface.polygons)
+        c = sum(poly) / length(poly)
+        dx = pt[1] - c[1]
+        dy = pt[2] - c[2]
+        xy_dist = dx^2 + dy^2
+
+        if xy_dist < best_dist
+            best_dist = xy_dist
+            plane = PlaneCut(c, surface.normals[i])
+            best_sd = signed_distance(plane, pt)
+        end
+    end
+
+    return best_sd > 0
+end

From c093e6a15c330b4d4b5f83eca6a698458021113e Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 13 Feb 2026 14:20:54 +0000
Subject: [PATCH 03/14] Fix partial_cut boundary normals and add comprehensive
 tests

- Reverse cut face node ordering when partial_cut=:positive to ensure
  correct outward boundary normal orientation
- Handle edge case where plane is outside mesh with partial_cut
- Add tests for partial_cut, depth_grid_to_surface, and layered_mesh

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>
---
 src/meshes/CutCellMeshes/cutting.jl |  47 +++++-
 test/cut_cell_meshes.jl             | 212 +++++++++++++++++++++++++++-
 2 files changed, 253 insertions(+), 6 deletions(-)

diff --git a/src/meshes/CutCellMeshes/cutting.jl b/src/meshes/CutCellMeshes/cutting.jl
index 83697113e..17b2e6dcc 100644
--- a/src/meshes/CutCellMeshes/cutting.jl
+++ b/src/meshes/CutCellMeshes/cutting.jl
@@ -140,7 +140,7 @@ function cut_mesh(mesh::UnstructuredMesh{3}, plane::PlaneCut{T};
         end
     end
 
-    if isempty(cut_cells)
+    if isempty(cut_cells) && partial_cut == :none
         if extra_out
             nf = number_of_faces(mesh)
             nb = number_of_boundary_faces(mesh)
@@ -155,6 +155,36 @@ function cut_mesh(mesh::UnstructuredMesh{3}, plane::PlaneCut{T};
         return mesh
     end
 
+    if isempty(cut_cells) && partial_cut != :none
+        # No cells are cut.  Check whether any cells would be discarded.
+        discard = partial_cut == :positive ? :negative : :positive
+        any_discarded = false
+        for c in 1:nc
+            side = classify_cell(mesh, c, plane)
+            if side == discard
+                any_discarded = true
+                break
+            end
+        end
+        if !any_discarded
+            # All cells are on the kept side → return unchanged mesh
+            if extra_out
+                nf = number_of_faces(mesh)
+                nb = number_of_boundary_faces(mesh)
+                info = Dict{String, Any}(
+                    "cell_index" => collect(1:nc),
+                    "face_index" => collect(1:nf),
+                    "boundary_face_index" => collect(1:nb),
+                    "new_faces" => Int[]
+                )
+                return (mesh, info)
+            end
+            return mesh
+        end
+        # Otherwise fall through to build_cut_mesh which will remove
+        # cells on the discarded side.
+    end
+
     # Create mutable node list and edge intersection cache
     new_node_points = copy(mesh.node_points)
     edge_cache = Dict{Tuple{Int, Int}, Int}()
@@ -685,11 +715,18 @@ function build_cut_mesh(
                 # Both sub-cells exist: cut face is interior between them
                 fi = add_interior_face!(info.cut_face_nodes, neg_cell, pos_cell; old_face = 0)
                 push!(new_faces_list, fi)
+            elseif partial_cut == :negative
+                # Keeping neg sub-cell: outward normal should point neg→pos
+                # (the original ordering already gives that)
+                if neg_cell != 0
+                    bi = add_boundary_face!(info.cut_face_nodes, neg_cell; old_bf = 0)
+                    push!(new_faces_list, bi)
+                end
             else
-                # Only one sub-cell: cut face becomes boundary
-                kept = partial_cut == :positive ? pos_cell : neg_cell
-                if kept != 0
-                    bi = add_boundary_face!(info.cut_face_nodes, kept; old_bf = 0)
+                # Keeping pos sub-cell: outward normal should point pos→neg,
+                # i.e. the reverse of the original neg→pos ordering
+                if pos_cell != 0
+                    bi = add_boundary_face!(reverse(info.cut_face_nodes), pos_cell; old_bf = 0)
                     push!(new_faces_list, bi)
                 end
             end
diff --git a/test/cut_cell_meshes.jl b/test/cut_cell_meshes.jl
index 9ac8d15cc..eed790be3 100644
--- a/test/cut_cell_meshes.jl
+++ b/test/cut_cell_meshes.jl
@@ -3,7 +3,7 @@ using Test
 using LinearAlgebra
 using StaticArrays
 
-import Jutul.CutCellMeshes: PlaneCut, PolygonalSurface, cut_mesh
+import Jutul.CutCellMeshes: PlaneCut, PolygonalSurface, cut_mesh, layered_mesh, depth_grid_to_surface
 
 @testset "CutCellMeshes" begin
     @testset "PlaneCut construction" begin
@@ -454,4 +454,214 @@ import Jutul.CutCellMeshes: PlaneCut, PolygonalSurface, cut_mesh
         @test length(info["boundary_face_index"]) == number_of_boundary_faces(cut)
         @test length(info["new_faces"]) > 0
     end
+
+    @testset "partial_cut negative" begin
+        g = CartesianMesh((3, 3, 3))
+        mesh = UnstructuredMesh(g)
+        nc_orig = number_of_cells(mesh)
+
+        plane = PlaneCut([0.0, 0.0, 0.5], [0.0, 0.0, 1.0])
+        cut = cut_mesh(mesh, plane; min_cut_fraction = 0.01, partial_cut = :negative)
+
+        # Keeping negative side (below z=0.5): 9 bottom cells + 9 cut half-cells = 18
+        @test number_of_cells(cut) == 18
+        geo = tpfv_geometry(cut)
+        @test all(geo.volumes .> 0)
+        vol_orig = sum(tpfv_geometry(mesh).volumes)
+        vol_cut = sum(geo.volumes)
+        @test vol_cut < vol_orig
+        @test vol_cut ≈ 0.5 rtol=1e-8
+    end
+
+    @testset "partial_cut positive" begin
+        g = CartesianMesh((3, 3, 3))
+        mesh = UnstructuredMesh(g)
+
+        plane = PlaneCut([0.0, 0.0, 0.5], [0.0, 0.0, 1.0])
+        cut_pos = cut_mesh(mesh, plane; min_cut_fraction = 0.01, partial_cut = :positive)
+        cut_neg = cut_mesh(mesh, plane; min_cut_fraction = 0.01, partial_cut = :negative)
+
+        geo_pos = tpfv_geometry(cut_pos)
+        geo_neg = tpfv_geometry(cut_neg)
+        vol_orig = sum(tpfv_geometry(mesh).volumes)
+
+        @test all(geo_pos.volumes .> 0)
+        @test all(geo_neg.volumes .> 0)
+        # Volumes of positive and negative sides should sum to original
+        @test sum(geo_pos.volumes) + sum(geo_neg.volumes) ≈ vol_orig rtol=1e-8
+    end
+
+    @testset "partial_cut geometry consistency" begin
+        g = CartesianMesh((4, 4, 4))
+        mesh = UnstructuredMesh(g)
+        plane = PlaneCut([0.5, 0.5, 0.55], [0.0, 0.0, 1.0])
+
+        for side in [:negative, :positive]
+            cut = cut_mesh(mesh, plane; min_cut_fraction = 0.01, partial_cut = side)
+            geo = tpfv_geometry(cut)
+            @test all(geo.volumes .> 0)
+
+            # Interior normals point from left to right cell
+            for f in 1:number_of_faces(cut)
+                l, r = cut.faces.neighbors[f]
+                cl = geo.cell_centroids[:, l]
+                cr = geo.cell_centroids[:, r]
+                N = geo.normals[:, f]
+                @test dot(N, cr - cl) > 0
+            end
+
+            # Boundary normals point outward
+            for f in 1:number_of_boundary_faces(cut)
+                c = cut.boundary_faces.neighbors[f]
+                cc = geo.cell_centroids[:, c]
+                fc = geo.boundary_centroids[:, f]
+                N = geo.boundary_normals[:, f]
+                @test dot(N, fc - cc) > 0
+            end
+        end
+    end
+
+    @testset "partial_cut extra_out" begin
+        g = CartesianMesh((3, 3, 3))
+        mesh = UnstructuredMesh(g)
+        plane = PlaneCut([0.0, 0.0, 0.5], [0.0, 0.0, 1.0])
+
+        cut, info = cut_mesh(mesh, plane; min_cut_fraction = 0.01,
+            partial_cut = :negative, extra_out = true)
+
+        @test length(info["cell_index"]) == number_of_cells(cut)
+        @test length(info["face_index"]) == number_of_faces(cut)
+        @test length(info["boundary_face_index"]) == number_of_boundary_faces(cut)
+    end
+
+    @testset "partial_cut no-cut case" begin
+        g = CartesianMesh((2, 2, 2))
+        mesh = UnstructuredMesh(g)
+
+        # Plane outside mesh on negative side — all cells are positive
+        plane = PlaneCut([0.0, 0.0, -1.0], [0.0, 0.0, 1.0])
+
+        # partial_cut = :positive should keep all cells (they're all positive)
+        cut_pos = cut_mesh(mesh, plane; partial_cut = :positive)
+        @test number_of_cells(cut_pos) == number_of_cells(mesh)
+
+        # Plane outside mesh on positive side — all cells are negative
+        plane2 = PlaneCut([0.0, 0.0, 5.0], [0.0, 0.0, 1.0])
+
+        # partial_cut = :negative should keep all cells (they're all negative)
+        cut_neg = cut_mesh(mesh, plane2; partial_cut = :negative)
+        @test number_of_cells(cut_neg) == number_of_cells(mesh)
+    end
+
+    @testset "depth_grid_to_surface basic" begin
+        xs = collect(range(0.0, 1.0, length=4))
+        ys = collect(range(0.0, 1.0, length=4))
+        depths = fill(0.5, 4, 4)
+        surface = depth_grid_to_surface(xs, ys, depths)
+        @test length(surface.polygons) == 18  # 3x3 grid = 9 quads = 18 triangles
+        @test length(surface.normals) == 18
+    end
+
+    @testset "depth_grid_to_surface with NaN" begin
+        xs = collect(range(0.0, 1.0, length=3))
+        ys = collect(range(0.0, 1.0, length=3))
+        depths = fill(0.5, 3, 3)
+        depths[2, 2] = NaN
+        surface = depth_grid_to_surface(xs, ys, depths)
+        @test length(surface.polygons) < 8  # Less than full grid
+        @test length(surface.polygons) > 0
+    end
+
+    @testset "layered_mesh flat surfaces" begin
+        g = CartesianMesh((3, 3, 4), (1.0, 1.0, 1.0))
+        mesh = UnstructuredMesh(g)
+        vol_orig = sum(tpfv_geometry(mesh).volumes)
+
+        xs = collect(range(0.0, 1.0, length=4))
+        ys = collect(range(0.0, 1.0, length=4))
+        d1 = fill(0.3, 4, 4)
+        d2 = fill(0.7, 4, 4)
+        s1 = depth_grid_to_surface(xs, ys, d1)
+        s2 = depth_grid_to_surface(xs, ys, d2)
+
+        result, info = layered_mesh(mesh, [s1, s2])
+
+        geo = tpfv_geometry(result)
+        @test all(geo.volumes .> 0)
+        @test sum(geo.volumes) ≈ vol_orig rtol=1e-8
+
+        layer_vals = sort(unique(info["layer_indices"]))
+        @test length(layer_vals) == 3  # layers 0, 1, 2
+        @test layer_vals == [0, 1, 2]
+
+        @test length(info["layer_indices"]) == number_of_cells(result)
+        @test length(info["cell_index"]) == number_of_cells(result)
+    end
+
+    @testset "layered_mesh single surface" begin
+        g = CartesianMesh((3, 3, 4), (1.0, 1.0, 1.0))
+        mesh = UnstructuredMesh(g)
+
+        xs = collect(range(0.0, 1.0, length=4))
+        ys = collect(range(0.0, 1.0, length=4))
+        d = fill(0.5, 4, 4)
+        s = depth_grid_to_surface(xs, ys, d)
+
+        result, info = layered_mesh(mesh, [s])
+        layer_vals = sort(unique(info["layer_indices"]))
+        @test layer_vals == [0, 1]
+
+        vol_orig = sum(tpfv_geometry(mesh).volumes)
+        vol_result = sum(tpfv_geometry(result).volumes)
+        @test vol_result ≈ vol_orig rtol=1e-8
+    end
+
+    @testset "layered_mesh tilted surface" begin
+        g = CartesianMesh((4, 4, 6), (100.0, 100.0, 60.0))
+        mesh = UnstructuredMesh(g)
+        vol_orig = sum(tpfv_geometry(mesh).volumes)
+
+        xs = collect(range(0.0, 100.0, length=5))
+        ys = collect(range(0.0, 100.0, length=5))
+        d1 = [15.0 + (i - 1) * 2.5 for i in 1:5, j in 1:5]
+        d2 = fill(40.0, 5, 5)
+        s1 = depth_grid_to_surface(xs, ys, d1)
+        s2 = depth_grid_to_surface(xs, ys, d2)
+
+        result, info = layered_mesh(mesh, [s1, s2])
+
+        geo = tpfv_geometry(result)
+        @test all(geo.volumes .> 0)
+        @test sum(geo.volumes) ≈ vol_orig rtol=1e-6
+
+        layer_vals = sort(unique(info["layer_indices"]))
+        @test length(layer_vals) == 3
+    end
+
+    @testset "layered_mesh volume per layer" begin
+        g = CartesianMesh((3, 3, 4), (1.0, 1.0, 1.0))
+        mesh = UnstructuredMesh(g)
+
+        xs = collect(range(0.0, 1.0, length=4))
+        ys = collect(range(0.0, 1.0, length=4))
+        d1 = fill(0.3, 4, 4)
+        d2 = fill(0.7, 4, 4)
+        s1 = depth_grid_to_surface(xs, ys, d1)
+        s2 = depth_grid_to_surface(xs, ys, d2)
+
+        result, info = layered_mesh(mesh, [s1, s2])
+        geo = tpfv_geometry(result)
+        nc = number_of_cells(result)
+
+        # Layer 0: above z=0.3 → volume ≈ 0.3
+        # Layer 1: between z=0.3 and z=0.7 → volume ≈ 0.4
+        # Layer 2: below z=0.7 → volume ≈ 0.3
+        vol0 = sum(geo.volumes[i] for i in 1:nc if info["layer_indices"][i] == 0)
+        vol1 = sum(geo.volumes[i] for i in 1:nc if info["layer_indices"][i] == 1)
+        vol2 = sum(geo.volumes[i] for i in 1:nc if info["layer_indices"][i] == 2)
+
+        @test vol0 ≈ 0.3 rtol=1e-6
+        @test vol1 ≈ 0.4 rtol=1e-6
+        @test vol2 ≈ 0.3 rtol=1e-6
+    end
 end

From 912e061300212ceb166da35d721d0f42aa41b6e7 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 13 Feb 2026 14:23:13 +0000
Subject: [PATCH 04/14] Address code review: simplify _depth_valid, remove
 empty branches, add comment

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>
---
 src/meshes/CutCellMeshes/cutting.jl | 4 ----
 src/meshes/CutCellMeshes/layered.jl | 8 +++++---
 2 files changed, 5 insertions(+), 7 deletions(-)

diff --git a/src/meshes/CutCellMeshes/cutting.jl b/src/meshes/CutCellMeshes/cutting.jl
index 17b2e6dcc..c5a069325 100644
--- a/src/meshes/CutCellMeshes/cutting.jl
+++ b/src/meshes/CutCellMeshes/cutting.jl
@@ -693,10 +693,6 @@ function build_cut_mesh(
                     add_boundary_face!(fnodes, pos_cell; old_bf = bf)
                 elseif side < 0 && neg_cell != 0
                     add_boundary_face!(fnodes, neg_cell; old_bf = bf)
-                elseif side >= 0 && neg_cell != 0
-                    # Face is on positive side but only negative cell kept - skip
-                elseif side < 0 && pos_cell != 0
-                    # Face is on negative side but only positive cell kept - skip
                 end
             end
         end
diff --git a/src/meshes/CutCellMeshes/layered.jl b/src/meshes/CutCellMeshes/layered.jl
index 07c858ce2..c2babb0ce 100644
--- a/src/meshes/CutCellMeshes/layered.jl
+++ b/src/meshes/CutCellMeshes/layered.jl
@@ -51,7 +51,7 @@ function depth_grid_to_surface(
     return PolygonalSurface(polygons)
 end
 
-_depth_valid(z::Real) = isfinite(z) && !isnan(z)
+_depth_valid(z::Real) = isfinite(z)
 _depth_valid(::Missing) = false
 
 """
@@ -130,8 +130,10 @@ function layered_mesh(
         cz = geo.cell_centroids[3, c]
         centroid = SVector{3, T}(cx, cy, cz)
 
-        # Find the first surface this cell is above (positive side)
-        # layer = number of surfaces that the centroid is below
+        # Determine layer by counting how many surfaces the centroid is below.
+        # Surfaces are ordered by increasing depth, so once the centroid is
+        # above a surface we can stop: it cannot be below any deeper surface
+        # without contradicting the ordering requirement.
         layer = 0
         for (si, surface) in enumerate(surfaces)
             if _point_below_surface(centroid, surface)

From f047d468a04a3d3c6eabf1bcfd55e6c6ccd2869c Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 13 Feb 2026 19:17:25 +0000
Subject: [PATCH 05/14] Fix layered_mesh hanging: use bounding polygons and fix
 cut face construction

Two fixes for the layered_mesh hanging issue:

1. Use each polygon as a bounding_polygon when cutting by a PolygonalSurface.
   Previously each polygon's cutting plane extended to infinity, causing
   O(n_polygons) redundant cuts per cell and exponential mesh growth.
   Now each polygon only cuts cells within its spatial footprint.

2. Fix cut face construction for cells adjacent to already-cut neighbors.
   When a non-split face has nodes on the cutting plane (from a previous
   cut), those nodes must still be collected into the cut face polygon.
   Without this fix, the cut face was incomplete, causing volume loss.

Also simplified layered_mesh to delegate to cut_mesh(mesh, surface)
instead of manually iterating over individual polygons.

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>
---
 src/meshes/CutCellMeshes/cutting.jl | 49 +++++++++++++++++++++++++++--
 src/meshes/CutCellMeshes/layered.jl | 22 +++++--------
 test/cut_cell_meshes.jl             | 27 ++++++++++++++++
 3 files changed, 81 insertions(+), 17 deletions(-)

diff --git a/src/meshes/CutCellMeshes/cutting.jl b/src/meshes/CutCellMeshes/cutting.jl
index c5a069325..ffe4cc176 100644
--- a/src/meshes/CutCellMeshes/cutting.jl
+++ b/src/meshes/CutCellMeshes/cutting.jl
@@ -33,7 +33,10 @@ of the plane.
 
 Returns a new `UnstructuredMesh`, or `(UnstructuredMesh, Dict)` if `extra_out=true`.
 """
-function cut_mesh(mesh::UnstructuredMesh{3}, surface::PolygonalSurface; extra_out::Bool = false, kwargs...)
+function cut_mesh(mesh::UnstructuredMesh{3}, surface::PolygonalSurface;
+        extra_out::Bool = false,
+        min_cut_fraction::Real = 0.05,
+        partial_cut::Symbol = :none)
     result = mesh
     if extra_out
         nc = number_of_cells(mesh)
@@ -48,7 +51,13 @@ function cut_mesh(mesh::UnstructuredMesh{3}, surface::PolygonalSurface; extra_ou
             n = surface.normals[i]
             c = sum(poly) / length(poly)
             plane = PlaneCut(c, n)
-            result, step_info = cut_mesh(result, plane; extra_out = true, kwargs...)
+            bpoly = _expand_polygon(poly)
+            result, step_info = cut_mesh(result, plane;
+                extra_out = true,
+                min_cut_fraction = min_cut_fraction,
+                partial_cut = partial_cut,
+                bounding_polygon = bpoly,
+                clip_to_polygon = true)
             # Compose mappings: new cell → intermediate cell → original cell
             cell_idx = [cell_idx[j] for j in step_info["cell_index"]]
             face_idx = [j == 0 ? 0 : face_idx[j] for j in step_info["face_index"]]
@@ -70,12 +79,31 @@ function cut_mesh(mesh::UnstructuredMesh{3}, surface::PolygonalSurface; extra_ou
             n = surface.normals[i]
             c = sum(poly) / length(poly)
             plane = PlaneCut(c, n)
-            result = cut_mesh(result, plane; extra_out = false, kwargs...)
+            bpoly = _expand_polygon(poly)
+            result = cut_mesh(result, plane;
+                extra_out = false,
+                min_cut_fraction = min_cut_fraction,
+                partial_cut = partial_cut,
+                bounding_polygon = bpoly,
+                clip_to_polygon = true)
         end
         return result
     end
 end
 
+"""
+    _expand_polygon(poly, frac=0.02)
+
+Expand a polygon slightly outward from its centroid.  Each vertex is pushed
+away from the centroid by `frac` times the centroid-to-vertex distance.  This
+ensures that the point-in-polygon test reliably includes cells at shared edges
+and vertices of adjacent polygons.
+"""
+function _expand_polygon(poly::Vector{SVector{3, T}}, frac = 0.02) where T
+    c = sum(poly) / length(poly)
+    return [p + T(frac) * (p - c) for p in poly]
+end
+
 """
     cut_mesh(mesh, plane; partial_cut=:none, ...)
 
@@ -278,6 +306,16 @@ function cut_mesh(mesh::UnstructuredMesh{3}, plane::PlaneCut{T};
                 if side <= 0
                     push!(neg_faces, fnodes)
                 end
+                # Collect on-plane nodes even from non-split faces so that
+                # the cut face polygon is complete.  This matters when a
+                # neighbor was already split by a previous cut: the shared
+                # face may have nodes exactly on the plane that belong to
+                # the cut face boundary.
+                for n in fnodes
+                    if get(node_class, n, 0) == 0
+                        push!(all_plane_nodes, n)
+                    end
+                end
             end
         end
 
@@ -301,6 +339,11 @@ function cut_mesh(mesh::UnstructuredMesh{3}, plane::PlaneCut{T};
                 if side <= 0
                     push!(neg_faces, fnodes)
                 end
+                for n in fnodes
+                    if get(node_class, n, 0) == 0
+                        push!(all_plane_nodes, n)
+                    end
+                end
             end
         end
 
diff --git a/src/meshes/CutCellMeshes/layered.jl b/src/meshes/CutCellMeshes/layered.jl
index c2babb0ce..d6f850759 100644
--- a/src/meshes/CutCellMeshes/layered.jl
+++ b/src/meshes/CutCellMeshes/layered.jl
@@ -101,20 +101,14 @@ function layered_mesh(
     cell_map = collect(1:nc_orig)  # maps new cell → original cell
 
     for (si, surface) in enumerate(surfaces)
-        for (pi, poly) in enumerate(surface.polygons)
-            n = surface.normals[pi]
-            c = sum(poly) / length(poly)
-            plane = PlaneCut(c, n)
-
-            result, step_info = cut_mesh(current_mesh, plane;
-                extra_out = true,
-                min_cut_fraction = min_cut_fraction
-            )
-
-            # Compose cell maps
-            cell_map = [cell_map[j] for j in step_info["cell_index"]]
-            current_mesh = result
-        end
+        result, step_info = cut_mesh(current_mesh, surface;
+            extra_out = true,
+            min_cut_fraction = min_cut_fraction
+        )
+
+        # Compose cell maps
+        cell_map = [cell_map[j] for j in step_info["cell_index"]]
+        current_mesh = result
     end
 
     # ----------------------------------------------------------------
diff --git a/test/cut_cell_meshes.jl b/test/cut_cell_meshes.jl
index eed790be3..20ada629f 100644
--- a/test/cut_cell_meshes.jl
+++ b/test/cut_cell_meshes.jl
@@ -664,4 +664,31 @@ import Jutul.CutCellMeshes: PlaneCut, PolygonalSurface, cut_mesh, layered_mesh,
         @test vol1 ≈ 0.4 rtol=1e-6
         @test vol2 ≈ 0.3 rtol=1e-6
     end
+
+    @testset "layered_mesh perturbed surface performance" begin
+        # Regression test: this scenario used to hang or produce far too many
+        # cells because each polygon's cut plane extended to infinity.
+        g = CartesianMesh((10, 10, 10), (1000.0, 1000.0, 100.0))
+        mesh = UnstructuredMesh(g)
+
+        xs = collect(range(0.0, 1000.0, length=11))
+        ys = collect(range(0.0, 1000.0, length=11))
+        d1 = [33.0 + 5.0 * sin(0.1 * i + 0.2 * j) for i in 1:11, j in 1:11]
+        d2 = fill(61.5, 11, 11)
+        s1 = depth_grid_to_surface(xs, ys, d1)
+        s2 = depth_grid_to_surface(xs, ys, d2)
+
+        t = @elapsed begin
+            result, info = layered_mesh(mesh, [s1, s2])
+        end
+
+        geo = tpfv_geometry(result)
+        vol_orig = sum(tpfv_geometry(mesh).volumes)
+        vol_result = sum(geo.volumes)
+
+        @test all(geo.volumes .> 0)
+        @test vol_result ≈ vol_orig rtol=1e-6
+        @test sort(unique(info["layer_indices"])) == [0, 1, 2]
+        @test t < 60.0  # Must complete in reasonable time
+    end
 end

From e4450ceccfd2171d1d59fd3811673305e82537a7 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 13 Feb 2026 19:19:06 +0000
Subject: [PATCH 06/14] Address code review: tighten time limit, document
 return type

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>
---
 src/meshes/CutCellMeshes/cutting.jl | 2 +-
 test/cut_cell_meshes.jl             | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/meshes/CutCellMeshes/cutting.jl b/src/meshes/CutCellMeshes/cutting.jl
index ffe4cc176..b6650697a 100644
--- a/src/meshes/CutCellMeshes/cutting.jl
+++ b/src/meshes/CutCellMeshes/cutting.jl
@@ -92,7 +92,7 @@ function cut_mesh(mesh::UnstructuredMesh{3}, surface::PolygonalSurface;
 end
 
 """
-    _expand_polygon(poly, frac=0.02)
+    _expand_polygon(poly, frac=0.02) -> Vector{SVector{3, T}}
 
 Expand a polygon slightly outward from its centroid.  Each vertex is pushed
 away from the centroid by `frac` times the centroid-to-vertex distance.  This
diff --git a/test/cut_cell_meshes.jl b/test/cut_cell_meshes.jl
index 20ada629f..7e7deea10 100644
--- a/test/cut_cell_meshes.jl
+++ b/test/cut_cell_meshes.jl
@@ -689,6 +689,6 @@ import Jutul.CutCellMeshes: PlaneCut, PolygonalSurface, cut_mesh, layered_mesh,
         @test all(geo.volumes .> 0)
         @test vol_result ≈ vol_orig rtol=1e-6
         @test sort(unique(info["layer_indices"])) == [0, 1, 2]
-        @test t < 60.0  # Must complete in reasonable time
+        @test t < 30.0  # Must complete in reasonable time (was hanging before fix)
     end
 end

From b85e6e3be8b4ffd50dbfadc8456980bdd1fdd05e Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 13 Feb 2026 20:19:37 +0000
Subject: [PATCH 07/14] Add merge_coplanar_faces post-processing to cut cell
 meshes
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

- New merge_faces.jl: merge_coplanar_faces function that merges interior
  and boundary faces sharing ≥2 nodes, same cell pair/cell, and coplanar
  normals. Re-orders merged polygon nodes for correct orientation.
- New merge_faces keyword (default true) on cut_mesh for both PlaneCut
  and PolygonalSurface methods, passed through to build_cut_mesh.
- Tests for merge on/off, orientation, geometry preservation, extra_out.

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>
---
 src/meshes/CutCellMeshes/CutCellMeshes.jl |   1 +
 src/meshes/CutCellMeshes/cutting.jl       |  27 +-
 src/meshes/CutCellMeshes/merge_faces.jl   | 437 ++++++++++++++++++++++
 test/cut_cell_meshes.jl                   | 138 +++++++
 4 files changed, 596 insertions(+), 7 deletions(-)
 create mode 100644 src/meshes/CutCellMeshes/merge_faces.jl

diff --git a/src/meshes/CutCellMeshes/CutCellMeshes.jl b/src/meshes/CutCellMeshes/CutCellMeshes.jl
index 7f6e2ba12..b01abbf80 100644
--- a/src/meshes/CutCellMeshes/CutCellMeshes.jl
+++ b/src/meshes/CutCellMeshes/CutCellMeshes.jl
@@ -6,6 +6,7 @@ module CutCellMeshes
     include("types.jl")
     include("geometry.jl")
     include("cutting.jl")
+    include("merge_faces.jl")
     include("gluing.jl")
     include("layered.jl")
 
diff --git a/src/meshes/CutCellMeshes/cutting.jl b/src/meshes/CutCellMeshes/cutting.jl
index b6650697a..145e3c934 100644
--- a/src/meshes/CutCellMeshes/cutting.jl
+++ b/src/meshes/CutCellMeshes/cutting.jl
@@ -36,7 +36,8 @@ Returns a new `UnstructuredMesh`, or `(UnstructuredMesh, Dict)` if `extra_out=tr
 function cut_mesh(mesh::UnstructuredMesh{3}, surface::PolygonalSurface;
         extra_out::Bool = false,
         min_cut_fraction::Real = 0.05,
-        partial_cut::Symbol = :none)
+        partial_cut::Symbol = :none,
+        merge_faces::Bool = true)
     result = mesh
     if extra_out
         nc = number_of_cells(mesh)
@@ -57,7 +58,8 @@ function cut_mesh(mesh::UnstructuredMesh{3}, surface::PolygonalSurface;
                 min_cut_fraction = min_cut_fraction,
                 partial_cut = partial_cut,
                 bounding_polygon = bpoly,
-                clip_to_polygon = true)
+                clip_to_polygon = true,
+                merge_faces = merge_faces)
             # Compose mappings: new cell → intermediate cell → original cell
             cell_idx = [cell_idx[j] for j in step_info["cell_index"]]
             face_idx = [j == 0 ? 0 : face_idx[j] for j in step_info["face_index"]]
@@ -85,7 +87,8 @@ function cut_mesh(mesh::UnstructuredMesh{3}, surface::PolygonalSurface;
                 min_cut_fraction = min_cut_fraction,
                 partial_cut = partial_cut,
                 bounding_polygon = bpoly,
-                clip_to_polygon = true)
+                clip_to_polygon = true,
+                merge_faces = merge_faces)
         end
         return result
     end
@@ -118,7 +121,8 @@ function cut_mesh(mesh::UnstructuredMesh{3}, plane::PlaneCut{T};
         bounding_polygon = nothing,
         clip_to_polygon::Bool = false,
         extra_out::Bool = false,
-        partial_cut::Symbol = :none
+        partial_cut::Symbol = :none,
+        merge_faces::Bool = true
     ) where T
     @assert partial_cut in (:none, :positive, :negative) "partial_cut must be :none, :positive, or :negative"
     nc = number_of_cells(mesh)
@@ -389,7 +393,7 @@ function cut_mesh(mesh::UnstructuredMesh{3}, plane::PlaneCut{T};
         end
     end
 
-    return build_cut_mesh(mesh, plane, new_node_points, is_cut, cut_cell_infos, node_class, get_or_create_intersection, extra_out, partial_cut)
+    return build_cut_mesh(mesh, plane, new_node_points, is_cut, cut_cell_infos, node_class, get_or_create_intersection, extra_out, partial_cut, merge_faces)
 end
 
 """
@@ -418,7 +422,7 @@ function dominant_side(fnodes, node_class)
 end
 
 """
-    build_cut_mesh(mesh, plane, node_points, is_cut, cut_infos, node_class, get_intersection, extra_out, partial_cut)
+    build_cut_mesh(mesh, plane, node_points, is_cut, cut_infos, node_class, get_intersection, extra_out, partial_cut, do_merge_faces)
 
 Build the final mesh after cutting.
 
@@ -426,6 +430,10 @@ When `partial_cut` is `:positive` or `:negative`, only cells on the specified
 side are kept.  Cut cells lose the sub-cell on the discarded side and the cut
 face becomes a boundary face.  Uncut cells entirely on the discarded side are
 removed.
+
+When `do_merge_faces` is `true`, coplanar faces that share the same cell pair
+(or same cell for boundary faces) and share at least two nodes are merged into
+a single face.
 """
 function build_cut_mesh(
     mesh::UnstructuredMesh{3},
@@ -436,7 +444,8 @@ function build_cut_mesh(
     node_class::Dict{Int, Int},
     get_intersection::Function,
     extra_out::Bool,
-    partial_cut::Symbol = :none
+    partial_cut::Symbol = :none,
+    do_merge_faces::Bool = true
 ) where T
     nc_old = number_of_cells(mesh)
     nf_old = number_of_faces(mesh)
@@ -815,6 +824,10 @@ function build_cut_mesh(
         all_bnd_cells
     )
 
+    if do_merge_faces
+        new_mesh = merge_coplanar_faces(new_mesh)
+    end
+
     if extra_out
         info_dict = Dict{String, Any}(
             "cell_index" => cell_index,
diff --git a/src/meshes/CutCellMeshes/merge_faces.jl b/src/meshes/CutCellMeshes/merge_faces.jl
new file mode 100644
index 000000000..f5617a439
--- /dev/null
+++ b/src/meshes/CutCellMeshes/merge_faces.jl
@@ -0,0 +1,437 @@
+"""
+    merge_coplanar_faces(mesh; coplanar_tol=1e-8)
+
+Post-process a mesh by merging coplanar faces that share the same cell pair
+(for interior faces) or the same cell (for boundary faces) and share at least
+two nodes.  The merged polygon is re-ordered so that its normal points from
+left to right cell (interior) or outward from the cell (boundary).
+
+Returns a new `UnstructuredMesh`.
+"""
+function merge_coplanar_faces(
+    mesh::UnstructuredMesh{3};
+    coplanar_tol::Real = 1e-8
+)
+    nc = number_of_cells(mesh)
+    nf = number_of_faces(mesh)
+    nb = number_of_boundary_faces(mesh)
+    node_points = mesh.node_points
+
+    # ------------------------------------------------------------------
+    # 1. Merge interior faces
+    # ------------------------------------------------------------------
+    # Group faces by their (sorted) neighbor pair
+    pair_to_faces = Dict{Tuple{Int,Int}, Vector{Int}}()
+    for f in 1:nf
+        l, r = mesh.faces.neighbors[f]
+        key = l < r ? (l, r) : (r, l)
+        push!(get!(pair_to_faces, key, Int[]), f)
+    end
+
+    merged_face_nodes = Vector{Int}[]          # nodes per merged face
+    merged_face_neighbors = Tuple{Int,Int}[]   # (left, right) per merged face
+
+    face_done = falses(nf)
+    for (pair, fgroup) in pair_to_faces
+        _merge_face_group!(
+            merged_face_nodes, merged_face_neighbors,
+            fgroup, face_done, mesh, node_points, coplanar_tol, false
+        )
+    end
+    # Any un-merged face is added as-is
+    for f in 1:nf
+        if !face_done[f]
+            push!(merged_face_nodes, collect(mesh.faces.faces_to_nodes[f]))
+            push!(merged_face_neighbors, mesh.faces.neighbors[f])
+        end
+    end
+
+    # ------------------------------------------------------------------
+    # 2. Merge boundary faces
+    # ------------------------------------------------------------------
+    cell_to_bfaces = Dict{Int, Vector{Int}}()
+    for bf in 1:nb
+        c = mesh.boundary_faces.neighbors[bf]
+        push!(get!(cell_to_bfaces, c, Int[]), bf)
+    end
+
+    merged_bnd_nodes = Vector{Int}[]
+    merged_bnd_cells = Int[]
+
+    bnd_done = falses(nb)
+    for (c, bgroup) in cell_to_bfaces
+        _merge_face_group!(
+            merged_bnd_nodes, merged_bnd_cells,
+            bgroup, bnd_done, mesh, node_points, coplanar_tol, true
+        )
+    end
+    for bf in 1:nb
+        if !bnd_done[bf]
+            push!(merged_bnd_nodes, collect(mesh.boundary_faces.faces_to_nodes[bf]))
+            push!(merged_bnd_cells, mesh.boundary_faces.neighbors[bf])
+        end
+    end
+
+    # ------------------------------------------------------------------
+    # 3. Build the new mesh
+    # ------------------------------------------------------------------
+    return _rebuild_mesh(mesh, node_points, merged_face_nodes, merged_face_neighbors,
+                         merged_bnd_nodes, merged_bnd_cells)
+end
+
+# ==========================================================================
+#  Internal helpers
+# ==========================================================================
+
+"""
+    _merge_face_group!(out_nodes, out_meta, fgroup, done, mesh, node_points, tol, is_boundary)
+
+Given a group of faces that all belong to the same cell pair (interior) or
+same cell (boundary), find clusters of coplanar faces that share nodes and
+merge each cluster into a single face.
+"""
+function _merge_face_group!(
+    out_nodes::Vector{Vector{Int}},
+    out_meta,
+    fgroup::Vector{Int},
+    done::BitVector,
+    mesh::UnstructuredMesh{3},
+    node_points::Vector{SVector{3, T}},
+    tol::Real,
+    is_boundary::Bool
+) where T
+    nf = length(fgroup)
+    if nf <= 1
+        return  # nothing to merge
+    end
+
+    # Compute normal for each face
+    face_normals = Vector{SVector{3, T}}(undef, nf)
+    face_node_sets = Vector{Set{Int}}(undef, nf)
+    face_node_lists = Vector{Vector{Int}}(undef, nf)
+    for (i, f) in enumerate(fgroup)
+        if is_boundary
+            fnodes = collect(mesh.boundary_faces.faces_to_nodes[f])
+        else
+            fnodes = collect(mesh.faces.faces_to_nodes[f])
+        end
+        face_node_lists[i] = fnodes
+        face_node_sets[i] = Set(fnodes)
+        pts = [node_points[n] for n in fnodes]
+        face_normals[i] = polygon_normal(pts)
+    end
+
+    # Union-Find clustering: merge faces that are coplanar AND share ≥ 2 nodes
+    parent = collect(1:nf)
+    function find(x)
+        while parent[x] != x
+            parent[x] = parent[parent[x]]
+            x = parent[x]
+        end
+        return x
+    end
+    function unite(a, b)
+        ra, rb = find(a), find(b)
+        if ra != rb
+            parent[ra] = rb
+        end
+    end
+
+    for i in 1:nf
+        for j in (i+1):nf
+            shared = length(intersect(face_node_sets[i], face_node_sets[j]))
+            if shared >= 2
+                # Check coplanarity
+                n1 = face_normals[i]
+                n2 = face_normals[j]
+                if abs(abs(dot(n1, n2)) - 1) < tol
+                    unite(i, j)
+                end
+            end
+        end
+    end
+
+    # Group faces by cluster root
+    clusters = Dict{Int, Vector{Int}}()
+    for i in 1:nf
+        r = find(i)
+        push!(get!(clusters, r, Int[]), i)
+    end
+
+    for (_, cluster) in clusters
+        if length(cluster) == 1
+            continue  # nothing to merge; leave for the fallback loop
+        end
+        # Collect all nodes from the cluster
+        all_nodes = Set{Int}()
+        for idx in cluster
+            union!(all_nodes, face_node_sets[idx])
+        end
+
+        # Mark source faces as done
+        for idx in cluster
+            done[fgroup[idx]] = true
+        end
+
+        # Find interior nodes (shared by ≥ 2 faces in this cluster) that are
+        # not on the boundary of the merged polygon.
+        # The boundary of the merged polygon consists of edges that appear
+        # exactly once across all faces in the cluster.
+        edge_count = Dict{Tuple{Int,Int}, Int}()
+        for idx in cluster
+            fnodes = face_node_lists[idx]
+            nn = length(fnodes)
+            for k in 1:nn
+                a = fnodes[k]
+                b = fnodes[mod1(k + 1, nn)]
+                edge = a < b ? (a, b) : (b, a)
+                edge_count[edge] = get(edge_count, edge, 0) + 1
+            end
+        end
+
+        # Boundary edges appear exactly once
+        boundary_edges = Tuple{Int,Int}[]
+        for (edge, cnt) in edge_count
+            if cnt == 1
+                push!(boundary_edges, edge)
+            end
+        end
+
+        # Extract the boundary polygon by chaining boundary edges
+        merged_nodes = _chain_boundary_edges(boundary_edges)
+        if isempty(merged_nodes)
+            # Fallback: use all nodes sorted by angle
+            merged_nodes = collect(all_nodes)
+        end
+
+        # Order the polygon nodes correctly
+        pts = [node_points[n] for n in merged_nodes]
+        avg_normal = normalize(sum(face_normals[idx] for idx in cluster))
+        ordered_pts = order_polygon_points(pts, avg_normal)
+
+        # Map ordered points back to node indices
+        ordered_nodes = Int[]
+        for pt in ordered_pts
+            for n in merged_nodes
+                if node_points[n] ≈ pt
+                    push!(ordered_nodes, n)
+                    break
+                end
+            end
+        end
+
+        # Fix orientation for interior faces: normal should point from left to right
+        if !is_boundary
+            ref_face = fgroup[cluster[1]]
+            l, r = mesh.faces.neighbors[ref_face]
+            _fix_interior_orientation!(ordered_nodes, node_points, l, r, mesh)
+            push!(out_nodes, ordered_nodes)
+            push!(out_meta, (l, r))
+        else
+            ref_face = fgroup[cluster[1]]
+            c = mesh.boundary_faces.neighbors[ref_face]
+            _fix_boundary_orientation!(ordered_nodes, node_points, c, mesh)
+            push!(out_nodes, ordered_nodes)
+            push!(out_meta, c)
+        end
+    end
+end
+
+"""
+    _chain_boundary_edges(edges) -> Vector{Int}
+
+Given a set of directed/undirected edges, chain them into a single polygon
+boundary.  Returns an ordered list of node indices.
+"""
+function _chain_boundary_edges(edges::Vector{Tuple{Int,Int}})
+    if isempty(edges)
+        return Int[]
+    end
+    # Build adjacency: each node connects to its neighbors
+    adj = Dict{Int, Vector{Int}}()
+    for (a, b) in edges
+        push!(get!(adj, a, Int[]), b)
+        push!(get!(adj, b, Int[]), a)
+    end
+
+    # Chain starting from the first edge's first node
+    start = edges[1][1]
+    chain = [start]
+    visited = Set{Int}([start])
+
+    current = start
+    while true
+        neighbors = get(adj, current, Int[])
+        found = false
+        for n in neighbors
+            if !(n in visited)
+                push!(chain, n)
+                push!(visited, n)
+                current = n
+                found = true
+                break
+            end
+        end
+        if !found
+            break
+        end
+    end
+
+    return chain
+end
+
+"""
+    _fix_interior_orientation!(nodes, node_points, left, right, mesh)
+
+Ensure interior face normal (computed from node ordering) points from cell
+`left` toward cell `right`.  If not, reverse the node order.
+"""
+function _fix_interior_orientation!(
+    nodes::Vector{Int},
+    node_points::Vector{SVector{3, T}},
+    left::Int, right::Int,
+    mesh::UnstructuredMesh{3}
+) where T
+    pts = [node_points[n] for n in nodes]
+    face_normal = polygon_normal(pts)
+    face_centroid = sum(pts) / length(pts)
+
+    # Approximate cell centroids from their face centroids
+    left_centroid = _approx_cell_centroid(mesh, left)
+    right_centroid = _approx_cell_centroid(mesh, right)
+
+    # Normal should point from left to right
+    dir = right_centroid - left_centroid
+    if dot(face_normal, dir) < 0
+        reverse!(nodes)
+    end
+end
+
+"""
+    _fix_boundary_orientation!(nodes, node_points, cell, mesh)
+
+Ensure boundary face normal points outward (away from the cell).
+"""
+function _fix_boundary_orientation!(
+    nodes::Vector{Int},
+    node_points::Vector{SVector{3, T}},
+    cell::Int,
+    mesh::UnstructuredMesh{3}
+) where T
+    pts = [node_points[n] for n in nodes]
+    face_normal = polygon_normal(pts)
+    face_centroid = sum(pts) / length(pts)
+
+    cell_centroid = _approx_cell_centroid(mesh, cell)
+
+    # Normal should point away from cell (outward)
+    dir = face_centroid - cell_centroid
+    if dot(face_normal, dir) < 0
+        reverse!(nodes)
+    end
+end
+
+"""
+    _approx_cell_centroid(mesh, cell)
+
+Compute an approximate cell centroid by averaging the centroids of all
+its faces (interior + boundary).
+"""
+function _approx_cell_centroid(
+    mesh::UnstructuredMesh{3},
+    cell::Int
+)
+    T = eltype(eltype(mesh.node_points))
+    centroid = zero(SVector{3, T})
+    count = 0
+
+    for f in mesh.faces.cells_to_faces[cell]
+        fnodes = collect(mesh.faces.faces_to_nodes[f])
+        pts = [mesh.node_points[n] for n in fnodes]
+        centroid += sum(pts) / length(pts)
+        count += 1
+    end
+    for bf in mesh.boundary_faces.cells_to_faces[cell]
+        fnodes = collect(mesh.boundary_faces.faces_to_nodes[bf])
+        pts = [mesh.node_points[n] for n in fnodes]
+        centroid += sum(pts) / length(pts)
+        count += 1
+    end
+    if count > 0
+        centroid /= count
+    end
+    return centroid
+end
+
+"""
+    _rebuild_mesh(mesh, node_points, face_nodes, face_neighbors, bnd_nodes, bnd_cells)
+
+Build a new UnstructuredMesh from merged face data.
+"""
+function _rebuild_mesh(
+    mesh::UnstructuredMesh{3},
+    node_points::Vector{SVector{3, T}},
+    face_nodes_list::Vector{Vector{Int}},
+    face_neighbors::Vector{Tuple{Int,Int}},
+    bnd_nodes_list::Vector{Vector{Int}},
+    bnd_cells::Vector{Int}
+) where T
+    nc = number_of_cells(mesh)
+    nf_new = length(face_nodes_list)
+    nb_new = length(bnd_nodes_list)
+
+    # Build cell-to-face mappings
+    cell_int_faces = [Int[] for _ in 1:nc]
+    for (fi, (l, r)) in enumerate(face_neighbors)
+        push!(cell_int_faces[l], fi)
+        push!(cell_int_faces[r], fi)
+    end
+
+    cell_bnd_faces = [Int[] for _ in 1:nc]
+    for (bi, c) in enumerate(bnd_cells)
+        push!(cell_bnd_faces[c], bi)
+    end
+
+    # Flatten into indirection-map format
+    faces_nodes = Int[]
+    faces_nodespos = Int[1]
+    for fnodes in face_nodes_list
+        append!(faces_nodes, fnodes)
+        push!(faces_nodespos, faces_nodespos[end] + length(fnodes))
+    end
+
+    bnd_faces_nodes = Int[]
+    bnd_faces_nodespos = Int[1]
+    for fnodes in bnd_nodes_list
+        append!(bnd_faces_nodes, fnodes)
+        push!(bnd_faces_nodespos, bnd_faces_nodespos[end] + length(fnodes))
+    end
+
+    cells_faces = Int[]
+    cells_facepos = Int[1]
+    for c in 1:nc
+        append!(cells_faces, cell_int_faces[c])
+        push!(cells_facepos, cells_facepos[end] + length(cell_int_faces[c]))
+    end
+
+    bnd_cells_faces = Int[]
+    bnd_cells_facepos = Int[1]
+    for c in 1:nc
+        append!(bnd_cells_faces, cell_bnd_faces[c])
+        push!(bnd_cells_facepos, bnd_cells_facepos[end] + length(cell_bnd_faces[c]))
+    end
+
+    return UnstructuredMesh(
+        cells_faces,
+        cells_facepos,
+        bnd_cells_faces,
+        bnd_cells_facepos,
+        faces_nodes,
+        faces_nodespos,
+        bnd_faces_nodes,
+        bnd_faces_nodespos,
+        node_points,
+        face_neighbors,
+        bnd_cells
+    )
+end
diff --git a/test/cut_cell_meshes.jl b/test/cut_cell_meshes.jl
index 7e7deea10..d9681f5c6 100644
--- a/test/cut_cell_meshes.jl
+++ b/test/cut_cell_meshes.jl
@@ -691,4 +691,142 @@ import Jutul.CutCellMeshes: PlaneCut, PolygonalSurface, cut_mesh, layered_mesh,
         @test sort(unique(info["layer_indices"])) == [0, 1, 2]
         @test t < 30.0  # Must complete in reasonable time (was hanging before fix)
     end
+
+    @testset "merge_faces default on" begin
+        # Verify merge_faces=true is the default and produces valid meshes
+        g = CartesianMesh((3, 3, 3))
+        mesh = UnstructuredMesh(g)
+
+        plane = PlaneCut([0.0, 0.0, 0.5], [0.0, 0.0, 1.0])
+        cut = cut_mesh(mesh, plane; min_cut_fraction = 0.01)  # default merge_faces=true
+
+        geo = tpfv_geometry(cut)
+        vol_orig = sum(tpfv_geometry(mesh).volumes)
+        vol_cut = sum(geo.volumes)
+
+        @test all(geo.volumes .> 0)
+        @test vol_cut ≈ vol_orig rtol=1e-10
+
+        # Interior normals consistency
+        for f in 1:number_of_faces(cut)
+            l, r = cut.faces.neighbors[f]
+            cl = geo.cell_centroids[:, l]
+            cr = geo.cell_centroids[:, r]
+            N = geo.normals[:, f]
+            @test dot(N, cr - cl) > 0
+        end
+
+        # Boundary normals consistency
+        for f in 1:number_of_boundary_faces(cut)
+            c = cut.boundary_faces.neighbors[f]
+            cc = geo.cell_centroids[:, c]
+            fc = geo.boundary_centroids[:, f]
+            N = geo.boundary_normals[:, f]
+            @test dot(N, fc - cc) > 0
+        end
+    end
+
+    @testset "merge_faces=false produces valid mesh" begin
+        g = CartesianMesh((3, 3, 3))
+        mesh = UnstructuredMesh(g)
+
+        plane = PlaneCut([0.0, 0.0, 0.5], [0.0, 0.0, 1.0])
+        cut = cut_mesh(mesh, plane; min_cut_fraction = 0.01, merge_faces = false)
+
+        geo = tpfv_geometry(cut)
+        vol_orig = sum(tpfv_geometry(mesh).volumes)
+        vol_cut = sum(geo.volumes)
+
+        @test all(geo.volumes .> 0)
+        @test vol_cut ≈ vol_orig rtol=1e-10
+    end
+
+    @testset "merge_faces with PolygonalSurface" begin
+        g = CartesianMesh((3, 3, 3))
+        mesh = UnstructuredMesh(g)
+
+        poly = [
+            SVector{3, Float64}(-1.0, -1.0, 0.5),
+            SVector{3, Float64}(2.0, -1.0, 0.5),
+            SVector{3, Float64}(2.0, 2.0, 0.5),
+            SVector{3, Float64}(-1.0, 2.0, 0.5)
+        ]
+        surface = PolygonalSurface([poly])
+
+        cut_merged = cut_mesh(mesh, surface; min_cut_fraction = 0.01, merge_faces = true)
+        cut_nomerge = cut_mesh(mesh, surface; min_cut_fraction = 0.01, merge_faces = false)
+
+        # Both should have valid geometry
+        geo_m = tpfv_geometry(cut_merged)
+        geo_n = tpfv_geometry(cut_nomerge)
+        @test all(geo_m.volumes .> 0)
+        @test all(geo_n.volumes .> 0)
+        # Same cell count
+        @test number_of_cells(cut_merged) == number_of_cells(cut_nomerge)
+        # Merged should have <= faces (equal if no merges possible)
+        @test number_of_faces(cut_merged) <= number_of_faces(cut_nomerge)
+        @test number_of_boundary_faces(cut_merged) <= number_of_boundary_faces(cut_nomerge)
+    end
+
+    @testset "merge_coplanar_faces on synthetic mesh" begin
+        import Jutul.CutCellMeshes: merge_coplanar_faces
+
+        # Create a simple mesh, then manually verify merge_coplanar_faces works
+        # on a Cartesian mesh (which has no mergeable faces — it's a no-op test)
+        g = CartesianMesh((2, 2, 2))
+        mesh = UnstructuredMesh(g)
+        merged = merge_coplanar_faces(mesh)
+
+        @test number_of_cells(merged) == number_of_cells(mesh)
+        @test number_of_faces(merged) == number_of_faces(mesh)
+        @test number_of_boundary_faces(merged) == number_of_boundary_faces(mesh)
+
+        # Verify geometry is preserved
+        geo_orig = tpfv_geometry(mesh)
+        geo_merged = tpfv_geometry(merged)
+        @test sum(geo_orig.volumes) ≈ sum(geo_merged.volumes) rtol=1e-10
+    end
+
+    @testset "merge_coplanar_faces orientation" begin
+        import Jutul.CutCellMeshes: merge_coplanar_faces
+
+        # Use a cut mesh where merging occurs, then verify normals
+        g = CartesianMesh((4, 4, 4))
+        mesh = UnstructuredMesh(g)
+        plane = PlaneCut([0.0, 0.0, 0.5], [0.0, 0.0, 1.0])
+        cut = cut_mesh(mesh, plane; min_cut_fraction = 0.01, merge_faces = false)
+        merged = merge_coplanar_faces(cut)
+
+        geo = tpfv_geometry(merged)
+
+        # Interior normals should point from left to right cell
+        for f in 1:number_of_faces(merged)
+            l, r = merged.faces.neighbors[f]
+            cl = geo.cell_centroids[:, l]
+            cr = geo.cell_centroids[:, r]
+            N = geo.normals[:, f]
+            @test dot(N, cr - cl) > 0
+        end
+
+        # Boundary normals should point outward
+        for f in 1:number_of_boundary_faces(merged)
+            c = merged.boundary_faces.neighbors[f]
+            cc = geo.cell_centroids[:, c]
+            fc = geo.boundary_centroids[:, f]
+            N = geo.boundary_normals[:, f]
+            @test dot(N, fc - cc) > 0
+        end
+    end
+
+    @testset "merge_faces with extra_out" begin
+        g = CartesianMesh((3, 3, 3))
+        mesh = UnstructuredMesh(g)
+
+        plane = PlaneCut([0.0, 0.0, 0.5], [0.0, 0.0, 1.0])
+        cut, info = cut_mesh(mesh, plane; min_cut_fraction = 0.01, extra_out = true, merge_faces = true)
+
+        # Info dict should still work (cell_index refers to pre-merge cell counts)
+        @test length(info["cell_index"]) > 0
+        @test all(1 .<= info["cell_index"] .<= number_of_cells(mesh))
+    end
 end

From 22f4cad372428b1574e75735f3a2ec28bef91d30 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 13 Feb 2026 20:22:57 +0000
Subject: [PATCH 08/14] Address code review: rename union-find, optimize point
 lookup, add boundary merge test

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>
---
 src/meshes/CutCellMeshes/merge_faces.jl | 30 ++++++++++++++++---------
 test/cut_cell_meshes.jl                 | 30 +++++++++++++++++++++++++
 2 files changed, 50 insertions(+), 10 deletions(-)

diff --git a/src/meshes/CutCellMeshes/merge_faces.jl b/src/meshes/CutCellMeshes/merge_faces.jl
index f5617a439..7535f9bc9 100644
--- a/src/meshes/CutCellMeshes/merge_faces.jl
+++ b/src/meshes/CutCellMeshes/merge_faces.jl
@@ -123,15 +123,15 @@ function _merge_face_group!(
 
     # Union-Find clustering: merge faces that are coplanar AND share ≥ 2 nodes
     parent = collect(1:nf)
-    function find(x)
+    function uf_find(x)
         while parent[x] != x
             parent[x] = parent[parent[x]]
             x = parent[x]
         end
         return x
     end
-    function unite(a, b)
-        ra, rb = find(a), find(b)
+    function uf_unite(a, b)
+        ra, rb = uf_find(a), uf_find(b)
         if ra != rb
             parent[ra] = rb
         end
@@ -145,7 +145,7 @@ function _merge_face_group!(
                 n1 = face_normals[i]
                 n2 = face_normals[j]
                 if abs(abs(dot(n1, n2)) - 1) < tol
-                    unite(i, j)
+                    uf_unite(i, j)
                 end
             end
         end
@@ -154,7 +154,7 @@ function _merge_face_group!(
     # Group faces by cluster root
     clusters = Dict{Int, Vector{Int}}()
     for i in 1:nf
-        r = find(i)
+        r = uf_find(i)
         push!(get!(clusters, r, Int[]), i)
     end
 
@@ -209,13 +209,23 @@ function _merge_face_group!(
         avg_normal = normalize(sum(face_normals[idx] for idx in cluster))
         ordered_pts = order_polygon_points(pts, avg_normal)
 
-        # Map ordered points back to node indices
+        # Map ordered points back to node indices using a point→node dictionary
+        pt_to_node = Dict{SVector{3, T}, Int}()
+        for n in merged_nodes
+            pt_to_node[node_points[n]] = n
+        end
         ordered_nodes = Int[]
         for pt in ordered_pts
-            for n in merged_nodes
-                if node_points[n] ≈ pt
-                    push!(ordered_nodes, n)
-                    break
+            n = get(pt_to_node, pt, 0)
+            if n != 0
+                push!(ordered_nodes, n)
+            else
+                # Fallback: linear scan with approximate match
+                for mn in merged_nodes
+                    if node_points[mn] ≈ pt
+                        push!(ordered_nodes, mn)
+                        break
+                    end
                 end
             end
         end
diff --git a/test/cut_cell_meshes.jl b/test/cut_cell_meshes.jl
index d9681f5c6..e2d370a11 100644
--- a/test/cut_cell_meshes.jl
+++ b/test/cut_cell_meshes.jl
@@ -787,6 +787,36 @@ import Jutul.CutCellMeshes: PlaneCut, PolygonalSurface, cut_mesh, layered_mesh,
         @test sum(geo_orig.volumes) ≈ sum(geo_merged.volumes) rtol=1e-10
     end
 
+    @testset "merge_coplanar_faces reduces boundary face count" begin
+        import Jutul.CutCellMeshes: merge_coplanar_faces
+
+        # A single z-cut on a 2x1x1 mesh creates cut cells whose boundary
+        # faces on the ±y and ±x sides are split into sub-faces that share
+        # the same cell AND are coplanar.  merge_coplanar_faces should
+        # recombine them.
+        g = CartesianMesh((1, 1, 1))
+        mesh = UnstructuredMesh(g)
+        plane = PlaneCut([0.0, 0.0, 0.5], [0.0, 0.0, 1.0])
+        # cut without merging
+        cut_no = cut_mesh(mesh, plane; min_cut_fraction = 0.01, merge_faces = false)
+        nb_no = number_of_boundary_faces(cut_no)
+
+        # cut with merging (default)
+        cut_yes = cut_mesh(mesh, plane; min_cut_fraction = 0.01)
+        nb_yes = number_of_boundary_faces(cut_yes)
+
+        # Each cut sub-cell has 5 boundary faces (4 sides + 1 z-face) = 10
+        # but the 4 side-faces per sub-cell are unique (different cells).
+        # If any coplanar pairs exist they will be merged; if not, counts equal.
+        @test nb_yes <= nb_no
+        # Both meshes should have valid geometry
+        geo_no = tpfv_geometry(cut_no)
+        geo_yes = tpfv_geometry(cut_yes)
+        @test all(geo_no.volumes .> 0)
+        @test all(geo_yes.volumes .> 0)
+        @test sum(geo_no.volumes) ≈ sum(geo_yes.volumes) rtol=1e-10
+    end
+
     @testset "merge_coplanar_faces orientation" begin
         import Jutul.CutCellMeshes: merge_coplanar_faces
 

From 166890c1c9ef666f0486da3bf786662921ae7475 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Sat, 14 Feb 2026 17:30:54 +0000
Subject: [PATCH 09/14] Add embed_mesh function for embedding one mesh inside
 another

Algorithm:
1. Collect boundary face polygons/normals from mesh B
2. Cut mesh A with each B boundary face (keeping both sides)
3. Remove cells of cut-A whose centroids are inside B
4. Glue trimmed A with B using glue_mesh
5. Track cell origins (mesh_a/mesh_b) and original indices

Includes _remove_cells helper and _point_inside_mesh ray-casting test.

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>
---
 src/Jutul.jl                              |   2 +-
 src/meshes/CutCellMeshes/CutCellMeshes.jl |   2 +
 src/meshes/CutCellMeshes/embedding.jl     | 362 ++++++++++++++++++++++
 3 files changed, 365 insertions(+), 1 deletion(-)
 create mode 100644 src/meshes/CutCellMeshes/embedding.jl

diff --git a/src/Jutul.jl b/src/Jutul.jl
index 5d0dfc0f4..31d832b8f 100644
--- a/src/Jutul.jl
+++ b/src/Jutul.jl
@@ -156,7 +156,7 @@ module Jutul
 
     # Cut-cell meshes
     include("meshes/CutCellMeshes/CutCellMeshes.jl")
-    import Jutul.CutCellMeshes: cut_mesh, PlaneCut, PolygonalSurface, glue_mesh, cut_and_displace_mesh, layered_mesh, depth_grid_to_surface
+    import Jutul.CutCellMeshes: cut_mesh, PlaneCut, PolygonalSurface, glue_mesh, cut_and_displace_mesh, layered_mesh, depth_grid_to_surface, embed_mesh
 
     # This is to make Jutul simulators work nicely with nested ForwardDiff.
     JutulSimulateTag = ForwardDiff.Tag{typeof(simulate), <:JutulEntity}
diff --git a/src/meshes/CutCellMeshes/CutCellMeshes.jl b/src/meshes/CutCellMeshes/CutCellMeshes.jl
index b01abbf80..d97521347 100644
--- a/src/meshes/CutCellMeshes/CutCellMeshes.jl
+++ b/src/meshes/CutCellMeshes/CutCellMeshes.jl
@@ -9,7 +9,9 @@ module CutCellMeshes
     include("merge_faces.jl")
     include("gluing.jl")
     include("layered.jl")
+    include("embedding.jl")
 
     export cut_mesh, glue_mesh, cut_and_displace_mesh
     export layered_mesh, depth_grid_to_surface
+    export embed_mesh
 end
diff --git a/src/meshes/CutCellMeshes/embedding.jl b/src/meshes/CutCellMeshes/embedding.jl
new file mode 100644
index 000000000..36da77d09
--- /dev/null
+++ b/src/meshes/CutCellMeshes/embedding.jl
@@ -0,0 +1,362 @@
+"""
+    embed_mesh(mesh_a::UnstructuredMesh{3}, mesh_b::UnstructuredMesh{3}; kwargs...)
+
+Embed mesh B inside mesh A by removing overlapping cells from A, cutting
+partially overlapping cells, and gluing the two meshes together.
+
+The structure of mesh B is preserved exactly.  Cells in A that are entirely
+inside B are removed.  Cells in A that are partially cut by the boundary of B
+are split so that only the portion outside B remains.  The trimmed A and B are
+then stitched together with `glue_mesh`.
+
+# Keyword arguments
+- `extra_out::Bool = false`: If `true`, return `(mesh, info)` where `info` is a
+  `Dict{String, Any}` containing:
+  - `"cell_origin"`:  `Vector{Symbol}` — `:mesh_a` or `:mesh_b` for each cell.
+  - `"cell_index_a"`: `Vector{Int}` — maps each cell to its index in the original
+    mesh A (0 for cells from B).
+  - `"cell_index_b"`: `Vector{Int}` — maps each cell to its index in the original
+    mesh B (0 for cells from A).
+- `min_cut_fraction::Real = 0.01`: Passed to `cut_mesh` for volume filtering.
+- `tol::Real = 1e-6`: Node merging tolerance for gluing.
+- `face_tol::Real = 1e-4`: Face centroid proximity tolerance for gluing.
+- `coplanar_tol::Real = 1e-3`: Coplanarity tolerance for gluing.
+- `merge_faces::Bool = true`: Whether to merge coplanar faces after cutting.
+
+Returns a new `UnstructuredMesh{3}`, or `(UnstructuredMesh{3}, Dict)` when
+`extra_out=true`.
+"""
+function embed_mesh(
+    mesh_a::UnstructuredMesh{3},
+    mesh_b::UnstructuredMesh{3};
+    extra_out::Bool = false,
+    min_cut_fraction::Real = 0.01,
+    tol::Real = 1e-6,
+    face_tol::Real = 1e-4,
+    coplanar_tol::Real = 1e-3,
+    merge_faces::Bool = true
+)
+    T = Float64
+
+    # ------------------------------------------------------------------
+    # 1. Collect boundary face polygons and normals from mesh B
+    # ------------------------------------------------------------------
+    nb_b = number_of_boundary_faces(mesh_b)
+    bnd_polys = Vector{Vector{SVector{3, T}}}(undef, nb_b)
+    bnd_normals = Vector{SVector{3, T}}(undef, nb_b)
+    bnd_centroids = Vector{SVector{3, T}}(undef, nb_b)
+
+    for bf in 1:nb_b
+        fnodes = collect(mesh_b.boundary_faces.faces_to_nodes[bf])
+        pts = SVector{3, T}[mesh_b.node_points[n] for n in fnodes]
+        bnd_polys[bf] = pts
+        bnd_normals[bf] = polygon_normal(pts)
+        bnd_centroids[bf] = sum(pts) / length(pts)
+    end
+
+    # Compute bounding box of mesh B for quick rejection
+    b_lo = SVector{3, T}(T(Inf), T(Inf), T(Inf))
+    b_hi = SVector{3, T}(T(-Inf), T(-Inf), T(-Inf))
+    for pt in mesh_b.node_points
+        b_lo = min.(b_lo, pt)
+        b_hi = max.(b_hi, pt)
+    end
+
+    # ------------------------------------------------------------------
+    # 2. Cut A with each boundary face of B, keeping both sides.
+    #    This splits cells of A that straddle B's boundary without
+    #    removing anything yet.
+    # ------------------------------------------------------------------
+    nc_a = number_of_cells(mesh_a)
+    current_mesh = mesh_a
+    cell_map = collect(1:nc_a)
+
+    for bf in 1:nb_b
+        poly = bnd_polys[bf]
+        normal = bnd_normals[bf]
+        c = bnd_centroids[bf]
+
+        plane = PlaneCut(c, normal)
+        bpoly = _expand_polygon(poly)
+
+        current_mesh, step_info = cut_mesh(current_mesh, plane;
+            extra_out = true,
+            min_cut_fraction = min_cut_fraction,
+            bounding_polygon = bpoly,
+            clip_to_polygon = true,
+            partial_cut = :none,
+            merge_faces = merge_faces
+        )
+        cell_map = [cell_map[j] for j in step_info["cell_index"]]
+    end
+
+    # ------------------------------------------------------------------
+    # 4. Remove all cells whose centroids are inside B
+    # ------------------------------------------------------------------
+    nc_cut = number_of_cells(current_mesh)
+    geo = tpfv_geometry(current_mesh)
+    keep_cell = trues(nc_cut)
+
+    for c in 1:nc_cut
+        cx = geo.cell_centroids[1, c]
+        cy = geo.cell_centroids[2, c]
+        cz = geo.cell_centroids[3, c]
+        centroid = SVector{3, T}(cx, cy, cz)
+
+        # Quick bounding box test
+        if any(centroid .< b_lo) || any(centroid .> b_hi)
+            continue
+        end
+
+        if _point_inside_mesh(centroid, bnd_polys, bnd_normals, bnd_centroids)
+            keep_cell[c] = false
+        end
+    end
+
+    # If there are cells to remove, build trimmed mesh
+    trimmed_a = current_mesh
+    trimmed_cell_map = cell_map
+    if !all(keep_cell)
+        trimmed_a, trim_info = _remove_cells(current_mesh, keep_cell)
+        trimmed_cell_map = [cell_map[j] for j in trim_info]
+    end
+
+    # ------------------------------------------------------------------
+    # 5. Glue trimmed A with B
+    # ------------------------------------------------------------------
+    result, glue_info = glue_mesh(trimmed_a, mesh_b;
+        tol = tol,
+        face_tol = face_tol,
+        coplanar_tol = coplanar_tol,
+        extra_out = true
+    )
+
+    if extra_out
+        nc_result = number_of_cells(result)
+        nc_b = number_of_cells(mesh_b)
+
+        cell_origin = Vector{Symbol}(undef, nc_result)
+        cell_index_a = Vector{Int}(undef, nc_result)
+        cell_index_b = Vector{Int}(undef, nc_result)
+
+        for c in 1:nc_result
+            idx_a = glue_info["cell_index_a"][c]
+            idx_b = glue_info["cell_index_b"][c]
+            if idx_a > 0
+                cell_origin[c] = :mesh_a
+                cell_index_a[c] = trimmed_cell_map[idx_a]
+                cell_index_b[c] = 0
+            else
+                cell_origin[c] = :mesh_b
+                cell_index_a[c] = 0
+                cell_index_b[c] = idx_b
+            end
+        end
+
+        info = Dict{String, Any}(
+            "cell_origin" => cell_origin,
+            "cell_index_a" => cell_index_a,
+            "cell_index_b" => cell_index_b,
+        )
+        return (result, info)
+    end
+
+    return result
+end
+
+# ==========================================================================
+#  Helpers
+# ==========================================================================
+
+"""
+    _point_inside_mesh(pt, polys, normals, centroids)
+
+Test whether a 3D point is inside a closed mesh defined by boundary face
+polygons.  Uses ray-casting along the +x direction and counts crossings with
+boundary face triangles.
+"""
+function _point_inside_mesh(
+    pt::SVector{3, T},
+    polys::Vector{Vector{SVector{3, T}}},
+    normals::Vector{SVector{3, T}},
+    centroids::Vector{SVector{3, T}}
+) where T
+    crossings = 0
+    ray_dir = SVector{3, T}(1, 0, 0)
+
+    for (i, poly) in enumerate(polys)
+        n = normals[i]
+
+        # Ray-plane intersection: t = dot(n, centroid - pt) / dot(n, ray_dir)
+        denom = dot(n, ray_dir)
+        if abs(denom) < eps(T) * 100
+            continue  # Ray is parallel to face
+        end
+
+        t = dot(n, centroids[i] - pt) / denom
+        if t <= 0
+            continue  # Behind the ray origin
+        end
+
+        # Intersection point
+        hit = pt + t * ray_dir
+
+        # Check if hit point is inside the face polygon using
+        # projection to the face plane
+        if _point_in_face_polygon(hit, poly, n)
+            crossings += 1
+        end
+    end
+
+    return isodd(crossings)
+end
+
+"""
+    _point_in_face_polygon(pt, poly, normal)
+
+Test whether a 3D point (assumed to be on the polygon's plane) lies inside
+the polygon.  Projects to the polygon's local 2D frame and uses ray-casting.
+"""
+function _point_in_face_polygon(
+    pt::SVector{3, T},
+    poly::Vector{SVector{3, T}},
+    normal::SVector{3, T}
+) where T
+    # Build local 2D frame
+    ref = abs(normal[1]) < 0.9 ? SVector{3, T}(1, 0, 0) : SVector{3, T}(0, 1, 0)
+    u = normalize(cross(normal, ref))
+    v = cross(normal, u)
+
+    # Project polygon and point to 2D
+    c = sum(poly) / length(poly)
+    poly_2d = [SVector{2, T}(dot(p - c, u), dot(p - c, v)) for p in poly]
+    pt_2d = SVector{2, T}(dot(pt - c, u), dot(pt - c, v))
+
+    return point_in_polygon_2d(pt_2d, poly_2d)
+end
+
+"""
+    _remove_cells(mesh, keep) -> (new_mesh, cell_mapping)
+
+Remove cells from a mesh.  `keep[c]` is `true` for cells to retain.
+Returns the new mesh and a vector mapping new cell indices to old cell indices.
+"""
+function _remove_cells(
+    mesh::UnstructuredMesh{3},
+    keep::BitVector
+)
+    T = eltype(eltype(mesh.node_points))
+    nc_old = number_of_cells(mesh)
+    nf_old = number_of_faces(mesh)
+    nb_old = number_of_boundary_faces(mesh)
+
+    # Build old→new cell mapping
+    old_to_new = zeros(Int, nc_old)
+    new_count = 0
+    cell_mapping = Int[]  # new cell → old cell
+    for c in 1:nc_old
+        if keep[c]
+            new_count += 1
+            old_to_new[c] = new_count
+            push!(cell_mapping, c)
+        end
+    end
+
+    # Rebuild face data
+    all_face_nodes = Vector{Int}[]
+    all_face_neighbors = Tuple{Int, Int}[]
+    all_bnd_nodes = Vector{Int}[]
+    all_bnd_cells = Int[]
+
+    cell_int_faces = [Int[] for _ in 1:new_count]
+    cell_bnd_faces = [Int[] for _ in 1:new_count]
+
+    # Process interior faces
+    for f in 1:nf_old
+        l_old, r_old = mesh.faces.neighbors[f]
+        l_new = old_to_new[l_old]
+        r_new = old_to_new[r_old]
+
+        fnodes = collect(mesh.faces.faces_to_nodes[f])
+
+        if l_new > 0 && r_new > 0
+            # Both cells kept → interior face
+            push!(all_face_nodes, fnodes)
+            push!(all_face_neighbors, (l_new, r_new))
+            fi = length(all_face_nodes)
+            push!(cell_int_faces[l_new], fi)
+            push!(cell_int_faces[r_new], fi)
+        elseif l_new > 0
+            # Only left kept → boundary face
+            push!(all_bnd_nodes, fnodes)
+            push!(all_bnd_cells, l_new)
+            bi = length(all_bnd_nodes)
+            push!(cell_bnd_faces[l_new], bi)
+        elseif r_new > 0
+            # Only right kept → boundary face (reverse node order for outward normal)
+            push!(all_bnd_nodes, reverse(fnodes))
+            push!(all_bnd_cells, r_new)
+            bi = length(all_bnd_nodes)
+            push!(cell_bnd_faces[r_new], bi)
+        end
+        # else: both removed, skip
+    end
+
+    # Process boundary faces
+    for bf in 1:nb_old
+        c_old = mesh.boundary_faces.neighbors[bf]
+        c_new = old_to_new[c_old]
+        if c_new > 0
+            fnodes = collect(mesh.boundary_faces.faces_to_nodes[bf])
+            push!(all_bnd_nodes, fnodes)
+            push!(all_bnd_cells, c_new)
+            bi = length(all_bnd_nodes)
+            push!(cell_bnd_faces[c_new], bi)
+        end
+    end
+
+    # Flatten into CSR format
+    faces_nodes = Int[]
+    faces_nodespos = Int[1]
+    for fnodes in all_face_nodes
+        append!(faces_nodes, fnodes)
+        push!(faces_nodespos, faces_nodespos[end] + length(fnodes))
+    end
+
+    bnd_faces_nodes = Int[]
+    bnd_faces_nodespos = Int[1]
+    for fnodes in all_bnd_nodes
+        append!(bnd_faces_nodes, fnodes)
+        push!(bnd_faces_nodespos, bnd_faces_nodespos[end] + length(fnodes))
+    end
+
+    cells_faces = Int[]
+    cells_facepos = Int[1]
+    for c in 1:new_count
+        append!(cells_faces, cell_int_faces[c])
+        push!(cells_facepos, cells_facepos[end] + length(cell_int_faces[c]))
+    end
+
+    bnd_cells_faces = Int[]
+    bnd_cells_facepos = Int[1]
+    for c in 1:new_count
+        append!(bnd_cells_faces, cell_bnd_faces[c])
+        push!(bnd_cells_facepos, bnd_cells_facepos[end] + length(cell_bnd_faces[c]))
+    end
+
+    new_mesh = UnstructuredMesh(
+        cells_faces,
+        cells_facepos,
+        bnd_cells_faces,
+        bnd_cells_facepos,
+        faces_nodes,
+        faces_nodespos,
+        bnd_faces_nodes,
+        bnd_faces_nodespos,
+        mesh.node_points,
+        all_face_neighbors,
+        all_bnd_cells
+    )
+
+    return (new_mesh, cell_mapping)
+end

From 6f5ce8f697e27b11fbacfd5256cc7535d583d970 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Sat, 14 Feb 2026 17:33:01 +0000
Subject: [PATCH 10/14] Add tests for embed_mesh: aligned, misaligned,
 protruding, outside, no extra_out

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>
---
 test/cut_cell_meshes.jl | 168 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 168 insertions(+)

diff --git a/test/cut_cell_meshes.jl b/test/cut_cell_meshes.jl
index e2d370a11..b934966de 100644
--- a/test/cut_cell_meshes.jl
+++ b/test/cut_cell_meshes.jl
@@ -859,4 +859,172 @@ import Jutul.CutCellMeshes: PlaneCut, PolygonalSurface, cut_mesh, layered_mesh,
         @test length(info["cell_index"]) > 0
         @test all(1 .<= info["cell_index"] .<= number_of_cells(mesh))
     end
+
+    # ==================================================================
+    # embed_mesh tests
+    # ==================================================================
+    @testset "embed_mesh grid-aligned" begin
+        import Jutul.CutCellMeshes: embed_mesh
+
+        # B is entirely inside A and aligns with A's grid lines
+        g_a = CartesianMesh((4, 4, 4), (4.0, 4.0, 4.0))
+        mesh_a = UnstructuredMesh(g_a)
+
+        g_b = CartesianMesh((2, 2, 2), (2.0, 2.0, 2.0))
+        mesh_b_raw = UnstructuredMesh(g_b)
+        offset = SVector{3, Float64}(1.0, 1.0, 1.0)
+        new_points = [p + offset for p in mesh_b_raw.node_points]
+        mesh_b = Jutul.CutCellMeshes._rebuild_mesh_with_nodes(mesh_b_raw, new_points)
+
+        result, info = embed_mesh(mesh_a, mesh_b; extra_out = true)
+        geo = tpfv_geometry(result)
+
+        vol_a = sum(tpfv_geometry(mesh_a).volumes)
+        vol_result = sum(geo.volumes)
+
+        # Volume conservation
+        @test vol_result ≈ vol_a rtol=1e-8
+
+        # All volumes positive
+        @test all(geo.volumes .> 0)
+
+        # B cells preserved
+        n_b = count(x -> x == :mesh_b, info["cell_origin"])
+        @test n_b == number_of_cells(mesh_b)
+
+        # Total cells: 56 from A + 8 from B = 64
+        n_a = count(x -> x == :mesh_a, info["cell_origin"])
+        @test n_a + n_b == number_of_cells(result)
+
+        # Cell indices valid
+        for c in 1:number_of_cells(result)
+            if info["cell_origin"][c] == :mesh_a
+                @test 1 <= info["cell_index_a"][c] <= number_of_cells(mesh_a)
+                @test info["cell_index_b"][c] == 0
+            else
+                @test info["cell_index_a"][c] == 0
+                @test 1 <= info["cell_index_b"][c] <= number_of_cells(mesh_b)
+            end
+        end
+    end
+
+    @testset "embed_mesh misaligned" begin
+        import Jutul.CutCellMeshes: embed_mesh
+
+        # B does NOT align with A's grid lines, causing cells to be split
+        g_a = CartesianMesh((4, 4, 4), (4.0, 4.0, 4.0))
+        mesh_a = UnstructuredMesh(g_a)
+
+        g_b = CartesianMesh((2, 2, 2), (1.0, 1.0, 1.0))
+        mesh_b_raw = UnstructuredMesh(g_b)
+        offset = SVector{3, Float64}(0.5, 0.5, 0.5)
+        new_points = [p + offset for p in mesh_b_raw.node_points]
+        mesh_b = Jutul.CutCellMeshes._rebuild_mesh_with_nodes(mesh_b_raw, new_points)
+
+        result, info = embed_mesh(mesh_a, mesh_b; extra_out = true)
+        geo = tpfv_geometry(result)
+
+        vol_a = sum(tpfv_geometry(mesh_a).volumes)
+        vol_result = sum(geo.volumes)
+
+        # Volume conservation
+        @test vol_result ≈ vol_a rtol=1e-6
+
+        # All volumes positive
+        @test all(geo.volumes .> 0)
+
+        # B cells preserved exactly
+        n_b = count(x -> x == :mesh_b, info["cell_origin"])
+        @test n_b == number_of_cells(mesh_b)
+
+        # More A cells than original (splitting occurred)
+        n_a = count(x -> x == :mesh_a, info["cell_origin"])
+        @test n_a > number_of_cells(mesh_a) - number_of_cells(mesh_b)
+
+        # Interior normals consistency
+        for f in 1:number_of_faces(result)
+            l, r = result.faces.neighbors[f]
+            cl = geo.cell_centroids[:, l]
+            cr = geo.cell_centroids[:, r]
+            N = geo.normals[:, f]
+            @test dot(N, cr - cl) > 0
+        end
+    end
+
+    @testset "embed_mesh B protruding" begin
+        import Jutul.CutCellMeshes: embed_mesh
+
+        # B extends outside A on one side
+        g_a = CartesianMesh((4, 4, 4), (4.0, 4.0, 4.0))
+        mesh_a = UnstructuredMesh(g_a)
+
+        g_b = CartesianMesh((2, 2, 2), (2.0, 2.0, 2.0))
+        mesh_b_raw = UnstructuredMesh(g_b)
+        offset = SVector{3, Float64}(3.0, 1.0, 1.0)
+        new_points = [p + offset for p in mesh_b_raw.node_points]
+        mesh_b = Jutul.CutCellMeshes._rebuild_mesh_with_nodes(mesh_b_raw, new_points)
+
+        result, info = embed_mesh(mesh_a, mesh_b; extra_out = true)
+        geo = tpfv_geometry(result)
+
+        vol_a = sum(tpfv_geometry(mesh_a).volumes)
+        vol_b = sum(tpfv_geometry(mesh_b).volumes)
+        vol_result = sum(geo.volumes)
+
+        # Volume should be A + part of B outside A
+        expected_vol = vol_a + vol_b / 2
+        @test vol_result ≈ expected_vol rtol=1e-6
+
+        # All volumes positive
+        @test all(geo.volumes .> 0)
+
+        # B cells preserved exactly
+        n_b = count(x -> x == :mesh_b, info["cell_origin"])
+        @test n_b == number_of_cells(mesh_b)
+    end
+
+    @testset "embed_mesh without extra_out" begin
+        import Jutul.CutCellMeshes: embed_mesh
+
+        g_a = CartesianMesh((3, 3, 3), (3.0, 3.0, 3.0))
+        mesh_a = UnstructuredMesh(g_a)
+
+        g_b = CartesianMesh((1, 1, 1), (1.0, 1.0, 1.0))
+        mesh_b_raw = UnstructuredMesh(g_b)
+        offset = SVector{3, Float64}(1.0, 1.0, 1.0)
+        new_points = [p + offset for p in mesh_b_raw.node_points]
+        mesh_b = Jutul.CutCellMeshes._rebuild_mesh_with_nodes(mesh_b_raw, new_points)
+
+        # Without extra_out should just return mesh
+        result = embed_mesh(mesh_a, mesh_b)
+        @test isa(result, UnstructuredMesh)
+
+        vol_a = sum(tpfv_geometry(mesh_a).volumes)
+        vol_result = sum(tpfv_geometry(result).volumes)
+        @test vol_result ≈ vol_a rtol=1e-6
+    end
+
+    @testset "embed_mesh B fully outside A" begin
+        import Jutul.CutCellMeshes: embed_mesh
+
+        # B is completely outside A -- should just add B's cells
+        g_a = CartesianMesh((2, 2, 2), (2.0, 2.0, 2.0))
+        mesh_a = UnstructuredMesh(g_a)
+
+        g_b = CartesianMesh((2, 2, 2), (1.0, 1.0, 1.0))
+        mesh_b_raw = UnstructuredMesh(g_b)
+        offset = SVector{3, Float64}(5.0, 5.0, 5.0)
+        new_points = [p + offset for p in mesh_b_raw.node_points]
+        mesh_b = Jutul.CutCellMeshes._rebuild_mesh_with_nodes(mesh_b_raw, new_points)
+
+        result, info = embed_mesh(mesh_a, mesh_b; extra_out = true)
+
+        # Total cells = A cells + B cells (no overlap)
+        @test number_of_cells(result) == number_of_cells(mesh_a) + number_of_cells(mesh_b)
+
+        vol_a = sum(tpfv_geometry(mesh_a).volumes)
+        vol_b = sum(tpfv_geometry(mesh_b).volumes)
+        vol_result = sum(tpfv_geometry(result).volumes)
+        @test vol_result ≈ vol_a + vol_b rtol=1e-8
+    end
 end

From b80cba2852e891718dbca9b2de127b13c225e4ce Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Sat, 14 Feb 2026 17:35:37 +0000
Subject: [PATCH 11/14] Address code review: fix step numbering and add comment
 for epsilon scaling

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>
---
 src/meshes/CutCellMeshes/embedding.jl | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/meshes/CutCellMeshes/embedding.jl b/src/meshes/CutCellMeshes/embedding.jl
index 36da77d09..5bddce68f 100644
--- a/src/meshes/CutCellMeshes/embedding.jl
+++ b/src/meshes/CutCellMeshes/embedding.jl
@@ -91,7 +91,7 @@ function embed_mesh(
     end
 
     # ------------------------------------------------------------------
-    # 4. Remove all cells whose centroids are inside B
+    # 3. Remove all cells whose centroids are inside B
     # ------------------------------------------------------------------
     nc_cut = number_of_cells(current_mesh)
     geo = tpfv_geometry(current_mesh)
@@ -122,7 +122,7 @@ function embed_mesh(
     end
 
     # ------------------------------------------------------------------
-    # 5. Glue trimmed A with B
+    # 4. Glue trimmed A with B
     # ------------------------------------------------------------------
     result, glue_info = glue_mesh(trimmed_a, mesh_b;
         tol = tol,
@@ -189,8 +189,9 @@ function _point_inside_mesh(
 
         # Ray-plane intersection: t = dot(n, centroid - pt) / dot(n, ray_dir)
         denom = dot(n, ray_dir)
-        if abs(denom) < eps(T) * 100
-            continue  # Ray is parallel to face
+        # Skip faces nearly parallel to the ray (scaled epsilon for robustness)
+        if abs(denom) < 100 * eps(T)
+            continue
         end
 
         t = dot(n, centroids[i] - pt) / denom

From 89cdfc5172d188f09d14c7ee1a6ab43ae7c741d6 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Sat, 14 Feb 2026 21:17:16 +0000
Subject: [PATCH 12/14] Add merge_split_cells, merge_small_cells, and fix
 merge_coplanar_faces

- merge_split_cells: merges cells from same original cell in same layer
- merge_small_cells: merges cells below volume threshold with best neighbor
- merge_coplanar_faces: add same-plane check, convexity check, fix moved
  done-marking after convexity check
- layered_mesh gains merge_cells keyword (default true)
- _dedup_nodes helper removes consecutive duplicate nodes in face polygons

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>
---
 src/Jutul.jl                              |   2 +-
 src/meshes/CutCellMeshes/CutCellMeshes.jl |   1 +
 src/meshes/CutCellMeshes/layered.jl       |  13 +-
 src/meshes/CutCellMeshes/merge_cells.jl   | 333 ++++++++++++++++++++++
 src/meshes/CutCellMeshes/merge_faces.jl   |  62 +++-
 5 files changed, 399 insertions(+), 12 deletions(-)
 create mode 100644 src/meshes/CutCellMeshes/merge_cells.jl

diff --git a/src/Jutul.jl b/src/Jutul.jl
index 31d832b8f..168c61ca5 100644
--- a/src/Jutul.jl
+++ b/src/Jutul.jl
@@ -156,7 +156,7 @@ module Jutul
 
     # Cut-cell meshes
     include("meshes/CutCellMeshes/CutCellMeshes.jl")
-    import Jutul.CutCellMeshes: cut_mesh, PlaneCut, PolygonalSurface, glue_mesh, cut_and_displace_mesh, layered_mesh, depth_grid_to_surface, embed_mesh
+    import Jutul.CutCellMeshes: cut_mesh, PlaneCut, PolygonalSurface, glue_mesh, cut_and_displace_mesh, layered_mesh, depth_grid_to_surface, embed_mesh, merge_split_cells, merge_small_cells
 
     # This is to make Jutul simulators work nicely with nested ForwardDiff.
     JutulSimulateTag = ForwardDiff.Tag{typeof(simulate), <:JutulEntity}
diff --git a/src/meshes/CutCellMeshes/CutCellMeshes.jl b/src/meshes/CutCellMeshes/CutCellMeshes.jl
index d97521347..94e5d09cf 100644
--- a/src/meshes/CutCellMeshes/CutCellMeshes.jl
+++ b/src/meshes/CutCellMeshes/CutCellMeshes.jl
@@ -7,6 +7,7 @@ module CutCellMeshes
     include("geometry.jl")
     include("cutting.jl")
     include("merge_faces.jl")
+    include("merge_cells.jl")
     include("gluing.jl")
     include("layered.jl")
     include("embedding.jl")
diff --git a/src/meshes/CutCellMeshes/layered.jl b/src/meshes/CutCellMeshes/layered.jl
index d6f850759..7e5a96897 100644
--- a/src/meshes/CutCellMeshes/layered.jl
+++ b/src/meshes/CutCellMeshes/layered.jl
@@ -88,7 +88,8 @@ below).
 function layered_mesh(
     mesh::UnstructuredMesh{3},
     surfaces::Vector{<:PolygonalSurface};
-    min_cut_fraction::Real = 0.01
+    min_cut_fraction::Real = 0.01,
+    merge_cells::Bool = true
 )
     T = Float64
     n_surfaces = length(surfaces)
@@ -144,6 +145,16 @@ function layered_mesh(
         "cell_index" => cell_map
     )
 
+    # ----------------------------------------------------------------
+    # 3. Optionally merge cells that came from the same original cell
+    #    and ended up in the same layer (unnecessarily split).
+    # ----------------------------------------------------------------
+    if merge_cells
+        current_mesh, info = merge_split_cells(current_mesh, info;
+            category_key = "layer_indices"
+        )
+    end
+
     return (current_mesh, info)
 end
 
diff --git a/src/meshes/CutCellMeshes/merge_cells.jl b/src/meshes/CutCellMeshes/merge_cells.jl
new file mode 100644
index 000000000..00259c99b
--- /dev/null
+++ b/src/meshes/CutCellMeshes/merge_cells.jl
@@ -0,0 +1,333 @@
+"""
+    merge_split_cells(mesh, info; category_key="layer_indices")
+
+Merge cells that were unnecessarily split during cutting.  Two cells are
+merged when they
+
+  1. originated from the same cell in the input mesh (`info["cell_index"]`),
+  2. belong to the same category (`info[category_key]`), and
+  3. share at least one interior face.
+
+The merge is transitive (union-find): if cells A–B and B–C satisfy the
+conditions they are all merged into one cell.
+
+Returns `(merged_mesh, merged_info)` where `merged_info` contains updated
+`"cell_index"` and `category_key` vectors.
+"""
+function merge_split_cells(
+    mesh::UnstructuredMesh{3},
+    info::Dict{String, Any};
+    category_key::String = "layer_indices"
+)
+    cell_index = info["cell_index"]::Vector{Int}
+    categories = info[category_key]::Vector{Int}
+    nc = number_of_cells(mesh)
+    @assert length(cell_index) == nc
+    @assert length(categories) == nc
+
+    # Union-Find
+    parent = collect(1:nc)
+    function uf_find(x)
+        while parent[x] != x
+            parent[x] = parent[parent[x]]
+            x = parent[x]
+        end
+        return x
+    end
+    function uf_unite(a, b)
+        ra, rb = uf_find(a), uf_find(b)
+        if ra != rb
+            parent[ra] = rb
+        end
+    end
+
+    # Unite cells that share a face, same original cell, same category
+    nf = number_of_faces(mesh)
+    for f in 1:nf
+        l, r = mesh.faces.neighbors[f]
+        if cell_index[l] == cell_index[r] && categories[l] == categories[r]
+            uf_unite(l, r)
+        end
+    end
+
+    # Build groups
+    groups = Dict{Int, Vector{Int}}()
+    for c in 1:nc
+        r = uf_find(c)
+        push!(get!(groups, r, Int[]), c)
+    end
+
+    return _do_merge_cells(mesh, groups, cell_index, categories, category_key)
+end
+
+"""
+    merge_small_cells(mesh, original_mesh, info;
+                      threshold=0.1, category_key="layer_indices")
+
+Merge cells whose volume is less than `threshold` times their original cell
+volume with the neighbouring cell of the same category that shares the
+largest face area.
+
+`original_mesh` is the mesh before cutting.  `info` must contain
+`"cell_index"` (mapping to original cells) and a category vector under
+`category_key`.
+
+Returns `(merged_mesh, merged_info)`.
+"""
+function merge_small_cells(
+    mesh::UnstructuredMesh{3},
+    original_mesh::UnstructuredMesh{3},
+    info::Dict{String, Any};
+    threshold::Real = 0.1,
+    category_key::String = "layer_indices"
+)
+    cell_index = info["cell_index"]::Vector{Int}
+    categories = info[category_key]::Vector{Int}
+    nc = number_of_cells(mesh)
+
+    # Compute volumes
+    geo = tpfv_geometry(mesh)
+    geo_orig = tpfv_geometry(original_mesh)
+
+    # Original cell volumes
+    orig_vol = geo_orig.volumes
+
+    # Union-Find
+    parent = collect(1:nc)
+    function uf_find(x)
+        while parent[x] != x
+            parent[x] = parent[parent[x]]
+            x = parent[x]
+        end
+        return x
+    end
+    function uf_unite(a, b)
+        ra, rb = uf_find(a), uf_find(b)
+        if ra != rb
+            parent[ra] = rb
+        end
+    end
+
+    # Identify small cells and merge with best neighbour
+    nf = number_of_faces(mesh)
+
+    # Build face areas
+    face_areas = zeros(nf)
+    for f in 1:nf
+        face_areas[f] = geo.areas[f]
+    end
+
+    # Process cells in order of increasing volume so that small cells get
+    # absorbed first.
+    order = sortperm(geo.volumes)
+
+    for c in order
+        rc = uf_find(c)
+        # Check if this cell (or its merged group representative) is small
+        v = geo.volumes[c]
+        ov = orig_vol[cell_index[c]]
+        if ov > 0 && v / ov < threshold
+            # Find the neighbor with same category and largest shared face
+            best_neighbor = 0
+            best_area = 0.0
+            for f in mesh.faces.cells_to_faces[c]
+                l, r = mesh.faces.neighbors[f]
+                other = l == c ? r : l
+                ro = uf_find(other)
+                if ro != rc && categories[other] == categories[c]
+                    a = face_areas[f]
+                    if a > best_area
+                        best_area = a
+                        best_neighbor = other
+                    end
+                end
+            end
+            if best_neighbor != 0
+                uf_unite(c, best_neighbor)
+            end
+        end
+    end
+
+    # Build groups
+    groups = Dict{Int, Vector{Int}}()
+    for c in 1:nc
+        r = uf_find(c)
+        push!(get!(groups, r, Int[]), c)
+    end
+
+    return _do_merge_cells(mesh, groups, cell_index, categories, category_key)
+end
+
+# =========================================================================
+#  Shared implementation: build a new mesh from cell groups
+# =========================================================================
+
+"""
+    _do_merge_cells(mesh, groups, cell_index, categories, category_key)
+
+Build a new mesh by merging cells according to `groups`.  `groups` maps a
+representative cell index to the vector of all cells in that group.
+Interior faces between cells in the same group are removed; boundary faces
+and inter-group interior faces are preserved.
+"""
+function _do_merge_cells(
+    mesh::UnstructuredMesh{3},
+    groups::Dict{Int, Vector{Int}},
+    cell_index::Vector{Int},
+    categories::Vector{Int},
+    category_key::String
+)
+    T = eltype(eltype(mesh.node_points))
+    nc_old = number_of_cells(mesh)
+    nf_old = number_of_faces(mesh)
+    nb_old = number_of_boundary_faces(mesh)
+
+    # Old cell → new cell mapping
+    old_to_new = zeros(Int, nc_old)
+    new_cell_index = Int[]
+    new_categories = Int[]
+    new_cell_id = 0
+    for (rep, members) in groups
+        new_cell_id += 1
+        for c in members
+            old_to_new[c] = new_cell_id
+        end
+        # Use representative cell's metadata
+        push!(new_cell_index, cell_index[rep])
+        push!(new_categories, categories[rep])
+    end
+    nc_new = new_cell_id
+
+    # Collect faces
+    all_face_nodes = Vector{Int}[]
+    all_face_neighbors = Tuple{Int, Int}[]
+    all_bnd_nodes = Vector{Int}[]
+    all_bnd_cells = Int[]
+
+    cell_int_faces = [Int[] for _ in 1:nc_new]
+    cell_bnd_faces = [Int[] for _ in 1:nc_new]
+
+    # Interior faces: keep only if the two neighbors map to different new cells.
+    for f in 1:nf_old
+        l_old, r_old = mesh.faces.neighbors[f]
+        l_new = old_to_new[l_old]
+        r_new = old_to_new[r_old]
+        if l_new != r_new
+            fnodes = _dedup_nodes(collect(mesh.faces.faces_to_nodes[f]))
+            if length(fnodes) >= 3
+                push!(all_face_nodes, fnodes)
+                push!(all_face_neighbors, (l_new, r_new))
+                fi = length(all_face_nodes)
+                push!(cell_int_faces[l_new], fi)
+                push!(cell_int_faces[r_new], fi)
+            end
+        end
+    end
+
+    # Boundary faces: keep all, remap cell
+    for bf in 1:nb_old
+        c_old = mesh.boundary_faces.neighbors[bf]
+        c_new = old_to_new[c_old]
+        fnodes = _dedup_nodes(collect(mesh.boundary_faces.faces_to_nodes[bf]))
+        if length(fnodes) >= 3
+            push!(all_bnd_nodes, fnodes)
+            push!(all_bnd_cells, c_new)
+            bi = length(all_bnd_nodes)
+            push!(cell_bnd_faces[c_new], bi)
+        end
+    end
+
+    # Build mesh arrays
+    new_mesh = _rebuild_mesh_from_data(
+        mesh.node_points, nc_new,
+        all_face_nodes, all_face_neighbors,
+        all_bnd_nodes, all_bnd_cells,
+        cell_int_faces, cell_bnd_faces
+    )
+
+    new_info = Dict{String, Any}(
+        "cell_index" => new_cell_index,
+        category_key => new_categories
+    )
+    return (new_mesh, new_info)
+end
+
+"""
+    _rebuild_mesh_from_data(node_points, nc, face_nodes, face_neighbors,
+                            bnd_nodes, bnd_cells, cell_int_faces, cell_bnd_faces)
+
+Low-level helper: assemble an `UnstructuredMesh` from pre-computed arrays.
+"""
+function _rebuild_mesh_from_data(
+    node_points::Vector{SVector{3, T}},
+    nc::Int,
+    face_nodes_list::Vector{Vector{Int}},
+    face_neighbors::Vector{Tuple{Int, Int}},
+    bnd_nodes_list::Vector{Vector{Int}},
+    bnd_cells::Vector{Int},
+    cell_int_faces::Vector{Vector{Int}},
+    cell_bnd_faces::Vector{Vector{Int}}
+) where T
+    faces_nodes = Int[]
+    faces_nodespos = Int[1]
+    for fnodes in face_nodes_list
+        append!(faces_nodes, fnodes)
+        push!(faces_nodespos, faces_nodespos[end] + length(fnodes))
+    end
+
+    bnd_faces_nodes = Int[]
+    bnd_faces_nodespos = Int[1]
+    for fnodes in bnd_nodes_list
+        append!(bnd_faces_nodes, fnodes)
+        push!(bnd_faces_nodespos, bnd_faces_nodespos[end] + length(fnodes))
+    end
+
+    cells_faces = Int[]
+    cells_facepos = Int[1]
+    for c in 1:nc
+        append!(cells_faces, cell_int_faces[c])
+        push!(cells_facepos, cells_facepos[end] + length(cell_int_faces[c]))
+    end
+
+    bnd_cells_faces = Int[]
+    bnd_cells_facepos = Int[1]
+    for c in 1:nc
+        append!(bnd_cells_faces, cell_bnd_faces[c])
+        push!(bnd_cells_facepos, bnd_cells_facepos[end] + length(cell_bnd_faces[c]))
+    end
+
+    return UnstructuredMesh(
+        cells_faces,
+        cells_facepos,
+        bnd_cells_faces,
+        bnd_cells_facepos,
+        faces_nodes,
+        faces_nodespos,
+        bnd_faces_nodes,
+        bnd_faces_nodespos,
+        node_points,
+        face_neighbors,
+        bnd_cells
+    )
+end
+
+"""
+    _dedup_nodes(nodes) -> Vector{Int}
+
+Remove consecutive duplicate node indices from a polygon node list
+(treating the list as cyclic).
+"""
+function _dedup_nodes(nodes::Vector{Int})
+    n = length(nodes)
+    if n <= 1
+        return nodes
+    end
+    result = Int[]
+    for i in 1:n
+        j = mod1(i + 1, n)
+        if nodes[i] != nodes[j]
+            push!(result, nodes[i])
+        end
+    end
+    return result
+end
diff --git a/src/meshes/CutCellMeshes/merge_faces.jl b/src/meshes/CutCellMeshes/merge_faces.jl
index 7535f9bc9..6aa471a77 100644
--- a/src/meshes/CutCellMeshes/merge_faces.jl
+++ b/src/meshes/CutCellMeshes/merge_faces.jl
@@ -141,11 +141,20 @@ function _merge_face_group!(
         for j in (i+1):nf
             shared = length(intersect(face_node_sets[i], face_node_sets[j]))
             if shared >= 2
-                # Check coplanarity
+                # Check coplanarity: normals must be parallel AND faces
+                # must lie on the same plane (a shared node's position
+                # relative to the other face's plane must be near-zero).
                 n1 = face_normals[i]
                 n2 = face_normals[j]
                 if abs(abs(dot(n1, n2)) - 1) < tol
-                    uf_unite(i, j)
+                    # Verify same-plane: pick a node from face j and
+                    # check its distance from face i's plane.
+                    ref_pt = node_points[face_node_lists[i][1]]
+                    test_pt = node_points[face_node_lists[j][1]]
+                    plane_dist = abs(dot(n1, test_pt - ref_pt))
+                    if plane_dist < tol * max(1.0, norm(test_pt))
+                        uf_unite(i, j)
+                    end
                 end
             end
         end
@@ -168,13 +177,6 @@ function _merge_face_group!(
             union!(all_nodes, face_node_sets[idx])
         end
 
-        # Mark source faces as done
-        for idx in cluster
-            done[fgroup[idx]] = true
-        end
-
-        # Find interior nodes (shared by ≥ 2 faces in this cluster) that are
-        # not on the boundary of the merged polygon.
         # The boundary of the merged polygon consists of edges that appear
         # exactly once across all faces in the cluster.
         edge_count = Dict{Tuple{Int,Int}, Int}()
@@ -220,7 +222,6 @@ function _merge_face_group!(
             if n != 0
                 push!(ordered_nodes, n)
             else
-                # Fallback: linear scan with approximate match
                 for mn in merged_nodes
                     if node_points[mn] ≈ pt
                         push!(ordered_nodes, mn)
@@ -230,6 +231,16 @@ function _merge_face_group!(
             end
         end
 
+        # Only merge if the resulting polygon is convex on its plane
+        if !_is_convex_polygon([node_points[n] for n in ordered_nodes], avg_normal)
+            continue  # leave faces un-merged for the fallback loop
+        end
+
+        # Mark source faces as done
+        for idx in cluster
+            done[fgroup[idx]] = true
+        end
+
         # Fix orientation for interior faces: normal should point from left to right
         if !is_boundary
             ref_face = fgroup[cluster[1]]
@@ -445,3 +456,34 @@ function _rebuild_mesh(
         bnd_cells
     )
 end
+
+"""
+    _is_convex_polygon(pts, normal; tol=1e-10)
+
+Check whether a 3D planar polygon (given as ordered vertices) is convex when
+projected onto its plane.  `normal` is the face normal used to define the
+winding direction.  All cross products of consecutive edge pairs must point
+in the same direction as `normal` (or be zero for collinear edges).
+"""
+function _is_convex_polygon(
+    pts::Vector{SVector{3, T}},
+    normal::SVector{3, T};
+    tol::Real = 1e-10
+) where T
+    n = length(pts)
+    if n <= 3
+        return true  # triangles are always convex
+    end
+    for i in 1:n
+        j = mod1(i + 1, n)
+        k = mod1(i + 2, n)
+        e1 = pts[j] - pts[i]
+        e2 = pts[k] - pts[j]
+        c = cross(e1, e2)
+        d = dot(c, normal)
+        if d < -tol
+            return false
+        end
+    end
+    return true
+end

From cb07f934bcdd2e6257b27aa288b1a53c39479f7c Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Sat, 14 Feb 2026 21:19:31 +0000
Subject: [PATCH 13/14] Add tests for merge_split_cells, merge_small_cells, and
 convexity check

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>
---
 test/cut_cell_meshes.jl | 175 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 175 insertions(+)

diff --git a/test/cut_cell_meshes.jl b/test/cut_cell_meshes.jl
index b934966de..a04f44277 100644
--- a/test/cut_cell_meshes.jl
+++ b/test/cut_cell_meshes.jl
@@ -1027,4 +1027,179 @@ import Jutul.CutCellMeshes: PlaneCut, PolygonalSurface, cut_mesh, layered_mesh,
         vol_result = sum(tpfv_geometry(result).volumes)
         @test vol_result ≈ vol_a + vol_b rtol=1e-8
     end
+
+    # ==================================================================
+    # merge_split_cells tests
+    # ==================================================================
+    @testset "merge_split_cells basic" begin
+        import Jutul.CutCellMeshes: merge_split_cells
+
+        g = CartesianMesh((1, 1, 1), (100.0, 100.0, 100.0))
+        mesh = UnstructuredMesh(g)
+
+        N = 2
+        xs = range(0.0, 100.0, length=N)
+        ys = range(0.0, 100.0, length=N)
+        # Non-flat surface that causes over-splitting
+        s1 = depth_grid_to_surface(xs, ys, [50.0 50.0; 55.0 50.0])
+
+        result_no, info_no = layered_mesh(mesh, [s1]; merge_cells=false)
+        @test number_of_cells(result_no) > 2  # over-split
+
+        result_yes, info_yes = layered_mesh(mesh, [s1]; merge_cells=true)
+        @test number_of_cells(result_yes) == 2  # correctly merged
+
+        # Volume preserved
+        geo_orig = tpfv_geometry(mesh)
+        geo_yes = tpfv_geometry(result_yes)
+        @test sum(geo_yes.volumes) ≈ sum(geo_orig.volumes) rtol=1e-8
+        @test all(geo_yes.volumes .> 0)
+
+        # Info has correct keys
+        @test haskey(info_yes, "cell_index")
+        @test haskey(info_yes, "layer_indices")
+        @test length(info_yes["cell_index"]) == 2
+        @test length(info_yes["layer_indices"]) == 2
+    end
+
+    @testset "merge_split_cells reduces cells on perturbed surface" begin
+        import Jutul.CutCellMeshes: merge_split_cells
+
+        g = CartesianMesh((3, 3, 3), (100.0, 100.0, 100.0))
+        mesh = UnstructuredMesh(g)
+        N = 4
+        xs = range(0.0, 100.0, length=N)
+        ys = range(0.0, 100.0, length=N)
+        depths = [50.0 50.0 50.0 50.0;
+                  53.0 48.0 55.0 50.0;
+                  50.0 52.0 47.0 54.0;
+                  50.0 50.0 50.0 50.0]
+        s1 = depth_grid_to_surface(xs, ys, depths)
+
+        result_no, _ = layered_mesh(mesh, [s1]; merge_cells=false)
+        result_yes, info = layered_mesh(mesh, [s1]; merge_cells=true)
+
+        # Merge should significantly reduce cell count
+        @test number_of_cells(result_yes) < number_of_cells(result_no)
+
+        # Geometry is valid
+        geo = tpfv_geometry(result_yes)
+        geo_orig = tpfv_geometry(mesh)
+        @test sum(geo.volumes) ≈ sum(geo_orig.volumes) rtol=1e-6
+        @test all(geo.volumes .> 0)
+
+        # Interior normals point from L to R
+        for f in 1:number_of_faces(result_yes)
+            l, r = result_yes.faces.neighbors[f]
+            cl = geo.cell_centroids[:, l]
+            cr = geo.cell_centroids[:, r]
+            N_f = geo.normals[:, f]
+            @test dot(N_f, cr - cl) > 0
+        end
+    end
+
+    @testset "merge_split_cells no-op on flat surface" begin
+        import Jutul.CutCellMeshes: merge_split_cells
+
+        # Flat surface: all triangles are coplanar, so no over-splitting
+        g = CartesianMesh((3, 3, 3), (100.0, 100.0, 100.0))
+        mesh = UnstructuredMesh(g)
+        N = 4
+        xs = range(0.0, 100.0, length=N)
+        ys = range(0.0, 100.0, length=N)
+        s1 = depth_grid_to_surface(xs, ys, fill(50.0, N, N))
+
+        result_no, _ = layered_mesh(mesh, [s1]; merge_cells=false)
+        result_yes, _ = layered_mesh(mesh, [s1]; merge_cells=true)
+
+        # No change because no over-splitting
+        @test number_of_cells(result_yes) == number_of_cells(result_no)
+    end
+
+    # ==================================================================
+    # merge_small_cells tests
+    # ==================================================================
+    @testset "merge_small_cells basic" begin
+        import Jutul.CutCellMeshes: merge_small_cells
+
+        g = CartesianMesh((5, 5, 5), (100.0, 100.0, 100.0))
+        mesh = UnstructuredMesh(g)
+        N = 6
+        xs = range(0.0, 100.0, length=N)
+        ys = range(0.0, 100.0, length=N)
+        depths = fill(50.0, N, N)
+        for i in 1:N, j in 1:N
+            depths[i,j] += sin(i*0.5) * 5.0 + cos(j*0.3) * 3.0
+        end
+        s1 = depth_grid_to_surface(xs, ys, depths)
+
+        result, info = layered_mesh(mesh, [s1]; merge_cells=true)
+        merged, minfo = merge_small_cells(result, mesh, info; threshold=0.5)
+
+        # Should reduce cell count
+        @test number_of_cells(merged) <= number_of_cells(result)
+
+        # Volume preserved
+        geo_orig = tpfv_geometry(mesh)
+        geo_merged = tpfv_geometry(merged)
+        @test sum(geo_merged.volumes) ≈ sum(geo_orig.volumes) rtol=1e-6
+        @test all(geo_merged.volumes .> 0)
+
+        # Info has correct keys
+        @test haskey(minfo, "cell_index")
+        @test haskey(minfo, "layer_indices")
+    end
+
+    @testset "merge_small_cells does nothing with low threshold" begin
+        import Jutul.CutCellMeshes: merge_small_cells
+
+        g = CartesianMesh((3, 3, 3), (100.0, 100.0, 100.0))
+        mesh = UnstructuredMesh(g)
+        N = 4
+        xs = range(0.0, 100.0, length=N)
+        ys = range(0.0, 100.0, length=N)
+        s1 = depth_grid_to_surface(xs, ys, fill(50.0, N, N))
+
+        result, info = layered_mesh(mesh, [s1])
+        merged, _ = merge_small_cells(result, mesh, info; threshold=0.001)
+
+        # Very low threshold: nothing to merge
+        @test number_of_cells(merged) == number_of_cells(result)
+    end
+
+    # ==================================================================
+    # merge_coplanar_faces convexity check
+    # ==================================================================
+    @testset "merge_coplanar_faces convexity check" begin
+        import Jutul.CutCellMeshes: merge_coplanar_faces, _is_convex_polygon
+
+        # Test the convexity checker directly
+        # Convex square
+        pts_convex = [
+            SVector{3, Float64}(0.0, 0.0, 0.0),
+            SVector{3, Float64}(1.0, 0.0, 0.0),
+            SVector{3, Float64}(1.0, 1.0, 0.0),
+            SVector{3, Float64}(0.0, 1.0, 0.0)
+        ]
+        @test _is_convex_polygon(pts_convex, SVector{3, Float64}(0.0, 0.0, 1.0))
+
+        # Non-convex (L-shape approximation)
+        pts_nonconvex = [
+            SVector{3, Float64}(0.0, 0.0, 0.0),
+            SVector{3, Float64}(2.0, 0.0, 0.0),
+            SVector{3, Float64}(2.0, 1.0, 0.0),
+            SVector{3, Float64}(1.0, 1.0, 0.0),
+            SVector{3, Float64}(1.0, 2.0, 0.0),
+            SVector{3, Float64}(0.0, 2.0, 0.0)
+        ]
+        @test !_is_convex_polygon(pts_nonconvex, SVector{3, Float64}(0.0, 0.0, 1.0))
+
+        # Triangle is always convex
+        pts_tri = [
+            SVector{3, Float64}(0.0, 0.0, 0.0),
+            SVector{3, Float64}(1.0, 0.0, 0.0),
+            SVector{3, Float64}(0.5, 1.0, 0.0)
+        ]
+        @test _is_convex_polygon(pts_tri, SVector{3, Float64}(0.0, 0.0, 1.0))
+    end
 end

From 73ac486eb2163aab512f515cd12a3f09df8ad71c Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Sat, 14 Feb 2026 21:21:43 +0000
Subject: [PATCH 14/14] Address code review: use face diagonal for coplanarity
 tolerance scaling

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>
---
 src/meshes/CutCellMeshes/merge_faces.jl | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/meshes/CutCellMeshes/merge_faces.jl b/src/meshes/CutCellMeshes/merge_faces.jl
index 6aa471a77..22a21e6e9 100644
--- a/src/meshes/CutCellMeshes/merge_faces.jl
+++ b/src/meshes/CutCellMeshes/merge_faces.jl
@@ -148,11 +148,15 @@ function _merge_face_group!(
                 n2 = face_normals[j]
                 if abs(abs(dot(n1, n2)) - 1) < tol
                     # Verify same-plane: pick a node from face j and
-                    # check its distance from face i's plane.
+                    # check its distance from face i's plane.  Scale
+                    # tolerance by the face diagonal to be independent
+                    # of the coordinate system origin.
                     ref_pt = node_points[face_node_lists[i][1]]
                     test_pt = node_points[face_node_lists[j][1]]
                     plane_dist = abs(dot(n1, test_pt - ref_pt))
-                    if plane_dist < tol * max(1.0, norm(test_pt))
+                    diag = norm(test_pt - ref_pt)
+                    char_len = max(one(T), diag)
+                    if plane_dist < tol * char_len
                         uf_unite(i, j)
                     end
                 end