From 433b004229a88bd37adc96c39e4016a7216ce6d7 Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Thu, 14 Dec 2023 10:55:02 +0000
Subject: [PATCH 01/19] Voting algorithm: prioritise owned cells

Not yet implemented second cell search
---
 firedrake/mesh.py | 58 +++++++++++++++++++++++++++++++++++------------
 1 file changed, 43 insertions(+), 15 deletions(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 1bbbf32d91..c74beff8d2 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -3689,12 +3689,18 @@ def _parent_mesh_embedding(
     if parent_mesh.geometric_dimension() > parent_mesh.topological_dimension():
         # The reference coordinates contain an extra unnecessary dimension
         # which we can safely delete
-        reference_coords = reference_coords[:, :parent_mesh.topological_dimension()]
+        reference_coords = reference_coords[:, : parent_mesh.topological_dimension()]
 
     locally_visible[:] = parent_cell_nums != -1
     ranks[locally_visible] = visible_ranks[parent_cell_nums[locally_visible]]
     # see below for why np.inf is used here.
     ref_cell_dists_l1[~locally_visible] = np.inf
+
+    # ensure that points which a rank thinks it owns are always chosen in a tie
+    # break by setting the rank to be negative. If multiple ranks think they
+    # own a point then the one with the highest rank will be chosen.
+    on_this_rank = ranks == parent_mesh.comm.rank
+    ranks[on_this_rank] = -parent_mesh.comm.rank
     ref_cell_dists_l1_and_ranks = np.stack((ref_cell_dists_l1, ranks), axis=1)
 
     # In parallel there will regularly be disagreements about which cell owns a
@@ -3713,23 +3719,45 @@ def _parent_mesh_embedding(
         ref_cell_dists_l1_and_ranks, op=array_lexicographic_mpi_op
     )
 
+    # switch ranks back to positive
+    owned_ref_cell_dists_l1_and_ranks[:, 1] = np.abs(
+        owned_ref_cell_dists_l1_and_ranks[:, 1]
+    )
+    ref_cell_dists_l1_and_ranks[:, 1] = np.abs(ref_cell_dists_l1_and_ranks[:, 1])
+    ranks = np.abs(ranks)
+
+    owned_ref_cell_dists_l1 = owned_ref_cell_dists_l1_and_ranks[:, 0]
     owned_ranks = owned_ref_cell_dists_l1_and_ranks[:, 1]
 
-    # Any rows where owned_ref_cell_dists_l1_and_ranks and
-    # ref_cell_dists_l1_and_ranks differ in distance or rank correspond to
-    # points which are claimed by a cell that we cannot see. We should now
-    # update our information accordingly. This should only happen for points
-    # which we've already marked as being owned by a different rank.
-    extra_missing_points = ~np.all(
-        owned_ref_cell_dists_l1_and_ranks == ref_cell_dists_l1_and_ranks, axis=1
-    )
-    if any(owned_ranks[extra_missing_points] == parent_mesh.comm.rank):
-        raise RuntimeError(
-            "Some points have been claimed by a cell that we cannot see, "
-            "but which we think we own. This should not happen."
+    changed_ref_cell_dists_l1 = owned_ref_cell_dists_l1 != ref_cell_dists_l1
+    changed_ranks = owned_ranks != ranks
+
+    # If distance has changed the the point is not in local mesh partition
+    # since some other cell is closer.
+    locally_visible[changed_ref_cell_dists_l1] = False
+    parent_cell_nums[changed_ref_cell_dists_l1] = -1
+    # If the rank has changed but the distance hasn't then there was a tie
+    # break and we need to search for the point again, this time disallowing
+    # the previously identified cell: if we match the identified owned_rank AND
+    # the distance is the same then we have found the correct cell. If we
+    # cannot make a match to owned_rank and distance then we can't see the
+    # point. This should only happen for halo points.
+    changed_ranks_tied = changed_ranks & ~changed_ref_cell_dists_l1
+    if any(changed_ranks_tied):
+        # check that only halos are affected
+        if any(owned_ranks[changed_ranks_tied] == parent_mesh.comm.rank):
+            raise RuntimeError(
+                "A tie break has occurred in the voting algorithm which has "
+                "resulted in a non-halo point changing MPI rank. This should "
+                "not happen."
+            )
+        # do a second search for the points which have changed rank but not
+        # distance
+        raise NotImplementedError(
+            "Cell re-identification for points which have changed rank but "
+            "not distance is not yet implemented. This only affects halo "
+            "points."
         )
-    locally_visible[extra_missing_points] = False
-    parent_cell_nums[extra_missing_points] = -1
 
     # Any ranks which are still np.inf are not in the mesh
     missing_global_idxs = np.where(owned_ranks == np.inf)[0]

From c0b007aec60dceaf58d71a69a8085f04c055a13f Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Thu, 14 Dec 2023 11:12:50 +0000
Subject: [PATCH 02/19] update locate_cell.c to allow point ignoring

---
 firedrake/evaluate.h         |  3 ++-
 firedrake/locate.c           | 15 ++++++++++++++-
 firedrake/mesh.py            | 11 +++++++----
 firedrake/pointeval_utils.py |  2 +-
 4 files changed, 24 insertions(+), 7 deletions(-)

diff --git a/firedrake/evaluate.h b/firedrake/evaluate.h
index 4253c2174e..e10f38ae1b 100644
--- a/firedrake/evaluate.h
+++ b/firedrake/evaluate.h
@@ -52,7 +52,8 @@ extern int locate_cell(struct Function *f,
 		       ref_cell_l1_dist_xtr try_candidate_xtr,
 		       void *temp_ref_coords,
 		       void *found_ref_coords,
-		       double *found_ref_cell_dist_l1);
+		       double *found_ref_cell_dist_l1,
+			   int cell_ignore);
 
 extern int evaluate(struct Function *f,
 		    double *x,
diff --git a/firedrake/locate.c b/firedrake/locate.c
index 2588fa3ef3..db3c8a04a3 100644
--- a/firedrake/locate.c
+++ b/firedrake/locate.c
@@ -11,7 +11,8 @@ int locate_cell(struct Function *f,
         ref_cell_l1_dist_xtr try_candidate_xtr,
         void *temp_ref_coords,
         void *found_ref_coords,
-        double *found_ref_cell_dist_l1)
+        double *found_ref_cell_dist_l1,
+        int cell_ignore)
 {
     RTError err;
     int cell = -1;
@@ -42,6 +43,9 @@ int locate_cell(struct Function *f,
         if (f->extruded == 0) {
             for (uint64_t i = 0; i < nids; i++) {
                 current_ref_cell_dist_l1 = (*try_candidate)(temp_ref_coords, f, ids[i], x);
+                if (ids[i] == cell_ignore) {
+                    continue;
+                }
                 if (current_ref_cell_dist_l1 <= 0.0) {
                     /* Found cell! */
                     cell = ids[i];
@@ -67,6 +71,9 @@ int locate_cell(struct Function *f,
                 int c = ids[i] / nlayers;
                 int l = ids[i] % nlayers;
                 current_ref_cell_dist_l1 = (*try_candidate_xtr)(temp_ref_coords, f, c, l, x);
+                if (ids[i] == cell_ignore) {
+                    continue;
+                }
                 if (current_ref_cell_dist_l1 <= 0.0) {
                     /* Found cell! */
                     cell = ids[i];
@@ -91,6 +98,9 @@ int locate_cell(struct Function *f,
         if (f->extruded == 0) {
             for (int c = 0; c < f->n_cols; c++) {
                 current_ref_cell_dist_l1 = (*try_candidate)(temp_ref_coords, f, c, x);
+                if (c == cell_ignore) {
+                    continue;
+                }
                 if (current_ref_cell_dist_l1 <= 0.0) {
                     /* Found cell! */
                     cell = c;
@@ -114,6 +124,9 @@ int locate_cell(struct Function *f,
             for (int c = 0; c < f->n_cols; c++) {
                 for (int l = 0; l < f->n_layers; l++) {
                     current_ref_cell_dist_l1 = (*try_candidate_xtr)(temp_ref_coords, f, c, l, x);
+                    if (l == cell_ignore) {
+                        continue;
+                    }
                     if (current_ref_cell_dist_l1 <= 0.0) {
                         /* Found cell! */
                         cell = l;
diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index c74beff8d2..8a5983001c 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2225,12 +2225,14 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None):
         ref_cell_dists_l1 = np.empty(npoints, dtype=utils.RealType)
         cells = np.empty(npoints, dtype=IntType)
         assert xs.size == npoints * self.geometric_dimension()
+        cell_ignore = -1
         self._c_locator(tolerance=tolerance)(self.coordinates._ctypes,
                                              xs.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                                              Xs.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                                              ref_cell_dists_l1.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                                              cells.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
-                                             npoints)
+                                             npoints,
+                                             cell_ignore)
         return cells, Xs, ref_cell_dists_l1
 
     def _c_locator(self, tolerance=None):
@@ -2245,7 +2247,7 @@ def _c_locator(self, tolerance=None):
         except KeyError:
             src = pq_utils.src_locate_cell(self, tolerance=tolerance)
             src += """
-    int locator(struct Function *f, double *x, double *X, double *ref_cell_dists_l1, int *cells, size_t npoints)
+    int locator(struct Function *f, double *x, double *X, double *ref_cell_dists_l1, int *cells, size_t npoints, int cell_ignore)
     {
         size_t j = 0;  /* index into x and X */
         for(size_t i=0; i<npoints; i++) {
@@ -2258,7 +2260,7 @@ def _c_locator(self, tolerance=None):
             /* to_reference_coords and to_reference_coords_xtr are defined in
             pointquery_utils.py. If they contain python calls, this loop will
             not run at c-loop speed. */
-            cells[i] = locate_cell(f, &x[j], %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i]);
+            cells[i] = locate_cell(f, &x[j], %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i], cell_ignore);
 
             for (int k = 0; k < %(geometric_dimension)d; k++) {
                 X[j] = found_reference_coords.X[k];
@@ -2282,7 +2284,8 @@ def _c_locator(self, tolerance=None):
                                 ctypes.POINTER(ctypes.c_double),
                                 ctypes.POINTER(ctypes.c_double),
                                 ctypes.POINTER(ctypes.c_int),
-                                ctypes.c_size_t]
+                                ctypes.c_size_t,
+                                ctypes.c_int]
             locator.restype = ctypes.c_int
             return cache.setdefault(tolerance, locator)
 
diff --git a/firedrake/pointeval_utils.py b/firedrake/pointeval_utils.py
index ab6df3ddb4..61947887e5 100644
--- a/firedrake/pointeval_utils.py
+++ b/firedrake/pointeval_utils.py
@@ -139,7 +139,7 @@ def predicate(index):
     /* The type definitions and arguments used here are defined as statics in pointquery_utils.py */
     double found_ref_cell_dist_l1 = DBL_MAX;
     struct ReferenceCoords temp_reference_coords, found_reference_coords;
-    %(IntType)s cell = locate_cell(f, x, %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &found_ref_cell_dist_l1);
+    %(IntType)s cell = locate_cell(f, x, %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &found_ref_cell_dist_l1, -1);
     if (cell == -1) {
         return -1;
     }

From 280d0098629b1f20739b7bc3fdac2c307cdf978e Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Thu, 14 Dec 2023 12:01:26 +0000
Subject: [PATCH 03/19] update to use multiple points and test

---
 firedrake/mesh.py                    | 36 ++++++++++++++++++----------
 tests/regression/test_locate_cell.py | 18 ++++++++++++--
 2 files changed, 40 insertions(+), 14 deletions(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 8a5983001c..c11c15cfcf 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2144,7 +2144,7 @@ def spatial_index(self):
         return spatialindex.from_regions(coords_min, coords_max)
 
     @PETSc.Log.EventDecorator()
-    def locate_cell(self, x, tolerance=None):
+    def locate_cell(self, x, tolerance=None, cell_ignore=None):
         """Locate cell containing a given point.
 
         :arg x: point coordinates
@@ -2152,12 +2152,13 @@ def locate_cell(self, x, tolerance=None):
             Default is this mesh's :attr:`tolerance` property. Changing
             this from default will cause the spatial index to be rebuilt which
             can take some time.
+        :kwarg cell_ignore: Cell number to ignore in the search.
         :returns: cell number (int), or None (if the point is not
             in the domain)
         """
-        return self.locate_cell_and_reference_coordinate(x, tolerance=tolerance)[0]
+        return self.locate_cell_and_reference_coordinate(x, tolerance=tolerance, cell_ignore=cell_ignore)[0]
 
-    def locate_reference_coordinate(self, x, tolerance=None):
+    def locate_reference_coordinate(self, x, tolerance=None, cell_ignore=None):
         """Get reference coordinates of a given point in its cell. Which
         cell the point is in can be queried with the locate_cell method.
 
@@ -2166,12 +2167,13 @@ def locate_reference_coordinate(self, x, tolerance=None):
             Default is this mesh's :attr:`tolerance` property. Changing
             this from default will cause the spatial index to be rebuilt which
             can take some time.
+        :kwarg cell_ignore: Cell number to ignore in the search.
         :returns: reference coordinates within cell (numpy array) or
             None (if the point is not in the domain)
         """
-        return self.locate_cell_and_reference_coordinate(x, tolerance=tolerance)[1]
+        return self.locate_cell_and_reference_coordinate(x, tolerance=tolerance, cell_ignore=cell_ignore)[1]
 
-    def locate_cell_and_reference_coordinate(self, x, tolerance=None):
+    def locate_cell_and_reference_coordinate(self, x, tolerance=None, cell_ignore=None):
         """Locate cell containing a given point and the reference
         coordinates of the point within the cell.
 
@@ -2180,6 +2182,7 @@ def locate_cell_and_reference_coordinate(self, x, tolerance=None):
             Default is this mesh's :attr:`tolerance` property. Changing
             this from default will cause the spatial index to be rebuilt which
             can take some time.
+        :kwarg cell_ignore: Cell number to ignore in the search.
         :returns: tuple either
             (cell number, reference coordinates) of type (int, numpy array),
             or, when point is not in the domain, (None, None).
@@ -2188,12 +2191,12 @@ def locate_cell_and_reference_coordinate(self, x, tolerance=None):
         if x.size != self.geometric_dimension():
             raise ValueError("Point must have the same geometric dimension as the mesh")
         x = x.reshape((1, self.geometric_dimension()))
-        cells, ref_coords, _ = self.locate_cells_ref_coords_and_dists(x, tolerance=tolerance)
+        cells, ref_coords, _ = self.locate_cells_ref_coords_and_dists(x, tolerance=tolerance, cells_ignore=[cell_ignore])
         if cells[0] == -1:
             return None, None
         return cells[0], ref_coords[0]
 
-    def locate_cells_ref_coords_and_dists(self, xs, tolerance=None):
+    def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=None):
         """Locate cell containing a given point and the reference
         coordinates of the point within the cell.
 
@@ -2202,6 +2205,10 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None):
             Default is this mesh's :attr:`tolerance` property. Changing
             this from default will cause the spatial index to be rebuilt which
             can take some time.
+        :kwarg cells_ignore: Cell numbers to ignore in the search for each
+            point in xs. Each entry corresponds to a single coordinate in xs.
+            To not ignore any cells, pass None. To ensure a full cell search
+            for any given point, set the corresponding list entry to -1.
         :returns: tuple either
             (cell numbers array, reference coordinates array, ref_cell_dists_l1 array)
             of type
@@ -2222,17 +2229,22 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None):
             raise ValueError("Point coordinate dimension does not match mesh geometric dimension")
         Xs = np.empty_like(xs)
         npoints = len(xs)
+        if cells_ignore is None or cells_ignore[0] is None:
+            cells_ignore = np.full(npoints, -1, dtype=IntType)
+        else:
+            cells_ignore = np.asarray(cells_ignore, dtype=IntType)
+        if len(cells_ignore) != npoints:
+            raise ValueError("Number of cells to ignore does not match number of points")
         ref_cell_dists_l1 = np.empty(npoints, dtype=utils.RealType)
         cells = np.empty(npoints, dtype=IntType)
         assert xs.size == npoints * self.geometric_dimension()
-        cell_ignore = -1
         self._c_locator(tolerance=tolerance)(self.coordinates._ctypes,
                                              xs.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                                              Xs.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                                              ref_cell_dists_l1.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                                              cells.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
                                              npoints,
-                                             cell_ignore)
+                                             cells_ignore.ctypes.data_as(ctypes.POINTER(ctypes.c_int)))
         return cells, Xs, ref_cell_dists_l1
 
     def _c_locator(self, tolerance=None):
@@ -2247,7 +2259,7 @@ def _c_locator(self, tolerance=None):
         except KeyError:
             src = pq_utils.src_locate_cell(self, tolerance=tolerance)
             src += """
-    int locator(struct Function *f, double *x, double *X, double *ref_cell_dists_l1, int *cells, size_t npoints, int cell_ignore)
+    int locator(struct Function *f, double *x, double *X, double *ref_cell_dists_l1, int *cells, size_t npoints, int* cells_ignore)
     {
         size_t j = 0;  /* index into x and X */
         for(size_t i=0; i<npoints; i++) {
@@ -2260,7 +2272,7 @@ def _c_locator(self, tolerance=None):
             /* to_reference_coords and to_reference_coords_xtr are defined in
             pointquery_utils.py. If they contain python calls, this loop will
             not run at c-loop speed. */
-            cells[i] = locate_cell(f, &x[j], %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i], cell_ignore);
+            cells[i] = locate_cell(f, &x[j], %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i], cells_ignore[i]);
 
             for (int k = 0; k < %(geometric_dimension)d; k++) {
                 X[j] = found_reference_coords.X[k];
@@ -2285,7 +2297,7 @@ def _c_locator(self, tolerance=None):
                                 ctypes.POINTER(ctypes.c_double),
                                 ctypes.POINTER(ctypes.c_int),
                                 ctypes.c_size_t,
-                                ctypes.c_int]
+                                ctypes.POINTER(ctypes.c_int)]
             locator.restype = ctypes.c_int
             return cache.setdefault(tolerance, locator)
 
diff --git a/tests/regression/test_locate_cell.py b/tests/regression/test_locate_cell.py
index 21b35a7c66..f84bf42927 100644
--- a/tests/regression/test_locate_cell.py
+++ b/tests/regression/test_locate_cell.py
@@ -30,11 +30,25 @@ def meshdata(request):
 def test_locate_cell(meshdata, point, value):
     m, f = meshdata
 
-    def value_at(p):
-        cell = m.locate_cell(p)
+    def value_at(p, cell_ignore=None):
+        cell = m.locate_cell(p, cell_ignore=cell_ignore)
         return f.dat.data[cell]
 
+    def value_at_and_dist(p, cell_ignore=None):
+        if cell_ignore is not None:
+            cell_ignore = [cell_ignore]
+        cells, _, l1_dists = m.locate_cells_ref_coords_and_dists([p], cells_ignore=cell_ignore)
+        return f.dat.data[cells[0]], l1_dists[0]
+
     assert np.allclose(value, value_at(point))
+    cell = m.locate_cell(point)
+    assert ~np.allclose(value, value_at(point, cell_ignore=cell))
+    value_at, l1_dist = value_at_and_dist(point)
+    assert np.allclose(value, value_at)
+    assert np.isclose(l1_dist, 0.0)
+    value_at, l1_dist = value_at_and_dist(point, cell_ignore=cell)
+    assert ~np.allclose(value, value_at)
+    assert l1_dist > 0.0
 
 
 def test_locate_cell_not_found(meshdata):

From aa278e33619cb9e6d4f9bf4de70a93e72457c3d3 Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Thu, 14 Dec 2023 12:57:29 +0000
Subject: [PATCH 04/19] use new search ability

---
 firedrake/mesh.py | 35 ++++++++++++++++++++++++++++-------
 1 file changed, 28 insertions(+), 7 deletions(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index c11c15cfcf..c355ae4fef 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -3766,13 +3766,34 @@ def _parent_mesh_embedding(
                 "resulted in a non-halo point changing MPI rank. This should "
                 "not happen."
             )
-        # do a second search for the points which have changed rank but not
-        # distance
-        raise NotImplementedError(
-            "Cell re-identification for points which have changed rank but "
-            "not distance is not yet implemented. This only affects halo "
-            "points."
-        )
+        while any(changed_ranks_tied):
+            (
+                parent_cell_nums[changed_ranks_tied],
+                reference_coords[changed_ranks_tied, :],
+                ref_cell_dists_l1[changed_ranks_tied],
+            ) = parent_mesh.locate_cells_ref_coords_and_dists(
+                coords_global[changed_ranks_tied],
+                tolerance,
+                cells_ignore=parent_cell_nums[changed_ranks_tied],
+            )
+            locally_visible[changed_ranks_tied] = (
+                parent_cell_nums[changed_ranks_tied] != -1
+            )
+            changed_ranks_tied = changed_ranks_tied & locally_visible
+            # if new ref_cell_dists_l1 > owned_ref_cell_dists_l1 then we should
+            # disregard the point
+            locally_visible[changed_ranks_tied] = locally_visible[
+                changed_ranks_tied
+            ] & (
+                ref_cell_dists_l1[changed_ranks_tied]
+                <= owned_ref_cell_dists_l1[changed_ranks_tied]
+            )
+            changed_ranks_tied = changed_ranks_tied & locally_visible
+            # if the rank now matches then we have found the correct cell
+            locally_visible[changed_ranks_tied] = locally_visible[
+                changed_ranks_tied
+            ] & (owned_ranks[changed_ranks_tied] == ranks[changed_ranks_tied])
+            changed_ranks_tied = changed_ranks_tied & locally_visible
 
     # Any ranks which are still np.inf are not in the mesh
     missing_global_idxs = np.where(owned_ranks == np.inf)[0]

From df8eb23f48da16952418af94380daa64de190736 Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Thu, 14 Dec 2023 13:01:20 +0000
Subject: [PATCH 05/19] tidy

---
 firedrake/mesh.py | 16 +++++++---------
 1 file changed, 7 insertions(+), 9 deletions(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index c355ae4fef..60d9460584 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -3779,21 +3779,19 @@ def _parent_mesh_embedding(
             locally_visible[changed_ranks_tied] = (
                 parent_cell_nums[changed_ranks_tied] != -1
             )
-            changed_ranks_tied = changed_ranks_tied & locally_visible
+            changed_ranks_tied &= locally_visible
             # if new ref_cell_dists_l1 > owned_ref_cell_dists_l1 then we should
             # disregard the point
-            locally_visible[changed_ranks_tied] = locally_visible[
-                changed_ranks_tied
-            ] & (
+            locally_visible[changed_ranks_tied] &= (
                 ref_cell_dists_l1[changed_ranks_tied]
                 <= owned_ref_cell_dists_l1[changed_ranks_tied]
             )
-            changed_ranks_tied = changed_ranks_tied & locally_visible
+            changed_ranks_tied &= locally_visible
             # if the rank now matches then we have found the correct cell
-            locally_visible[changed_ranks_tied] = locally_visible[
-                changed_ranks_tied
-            ] & (owned_ranks[changed_ranks_tied] == ranks[changed_ranks_tied])
-            changed_ranks_tied = changed_ranks_tied & locally_visible
+            locally_visible[changed_ranks_tied] &= (
+                owned_ranks[changed_ranks_tied] == ranks[changed_ranks_tied]
+            )
+            changed_ranks_tied &= locally_visible
 
     # Any ranks which are still np.inf are not in the mesh
     missing_global_idxs = np.where(owned_ranks == np.inf)[0]

From aa5aaca7bd69b41685f66a9beccc97e530ffd672 Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Thu, 14 Dec 2023 17:05:44 +0000
Subject: [PATCH 06/19] Apparently working voting algorithm

---
 firedrake/evaluate.h                 |  3 +-
 firedrake/locate.c                   | 40 ++++++++++++++--
 firedrake/mesh.py                    | 68 +++++++++++++++++++---------
 firedrake/pointeval_utils.py         |  3 +-
 tests/regression/test_locate_cell.py |  2 +-
 5 files changed, 87 insertions(+), 29 deletions(-)

diff --git a/firedrake/evaluate.h b/firedrake/evaluate.h
index e10f38ae1b..738a573867 100644
--- a/firedrake/evaluate.h
+++ b/firedrake/evaluate.h
@@ -53,7 +53,8 @@ extern int locate_cell(struct Function *f,
 		       void *temp_ref_coords,
 		       void *found_ref_coords,
 		       double *found_ref_cell_dist_l1,
-			   int cell_ignore);
+			   size_t ncells_ignore,
+			   int* cells_ignore);
 
 extern int evaluate(struct Function *f,
 		    double *x,
diff --git a/firedrake/locate.c b/firedrake/locate.c
index db3c8a04a3..a243e2119f 100644
--- a/firedrake/locate.c
+++ b/firedrake/locate.c
@@ -12,10 +12,12 @@ int locate_cell(struct Function *f,
         void *temp_ref_coords,
         void *found_ref_coords,
         double *found_ref_cell_dist_l1,
-        int cell_ignore)
+        size_t ncells_ignore,
+        int* cells_ignore)
 {
     RTError err;
     int cell = -1;
+    int cell_ignore_found = 0;
     /* NOTE: temp_ref_coords and found_ref_coords are actually of type
     struct ReferenceCoords but can't be declared as such in the function
     signature because the dimensions of the reference coordinates in the
@@ -43,7 +45,14 @@ int locate_cell(struct Function *f,
         if (f->extruded == 0) {
             for (uint64_t i = 0; i < nids; i++) {
                 current_ref_cell_dist_l1 = (*try_candidate)(temp_ref_coords, f, ids[i], x);
-                if (ids[i] == cell_ignore) {
+                for (uint64_t j = 0; j < ncells_ignore; j++) {
+                    if (ids[i] == cells_ignore[j]) {
+                        cell_ignore_found = 1;
+                        break;
+                    }
+                }
+                if (cell_ignore_found) {
+                    cell_ignore_found = 0;
                     continue;
                 }
                 if (current_ref_cell_dist_l1 <= 0.0) {
@@ -71,7 +80,14 @@ int locate_cell(struct Function *f,
                 int c = ids[i] / nlayers;
                 int l = ids[i] % nlayers;
                 current_ref_cell_dist_l1 = (*try_candidate_xtr)(temp_ref_coords, f, c, l, x);
-                if (ids[i] == cell_ignore) {
+                for (uint64_t j = 0; j < ncells_ignore; j++) {
+                    if (ids[i] == cells_ignore[j]) {
+                        cell_ignore_found = 1;
+                        break;
+                    }
+                }
+                if (cell_ignore_found) {
+                    cell_ignore_found = 0;
                     continue;
                 }
                 if (current_ref_cell_dist_l1 <= 0.0) {
@@ -98,7 +114,14 @@ int locate_cell(struct Function *f,
         if (f->extruded == 0) {
             for (int c = 0; c < f->n_cols; c++) {
                 current_ref_cell_dist_l1 = (*try_candidate)(temp_ref_coords, f, c, x);
-                if (c == cell_ignore) {
+                for (uint64_t j = 0; j < ncells_ignore; j++) {
+                    if (c == cells_ignore[j]) {
+                        cell_ignore_found = 1;
+                        break;
+                    }
+                }
+                if (cell_ignore_found) {
+                    cell_ignore_found = 0;
                     continue;
                 }
                 if (current_ref_cell_dist_l1 <= 0.0) {
@@ -124,7 +147,14 @@ int locate_cell(struct Function *f,
             for (int c = 0; c < f->n_cols; c++) {
                 for (int l = 0; l < f->n_layers; l++) {
                     current_ref_cell_dist_l1 = (*try_candidate_xtr)(temp_ref_coords, f, c, l, x);
-                    if (l == cell_ignore) {
+                    for (uint64_t j = 0; j < ncells_ignore; j++) {
+                        if (l == cells_ignore[j]) {
+                            cell_ignore_found = 1;
+                            break;
+                        }
+                    }
+                    if (cell_ignore_found) {
+                        cell_ignore_found = 0;
                         continue;
                     }
                     if (current_ref_cell_dist_l1 <= 0.0) {
diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 60d9460584..1ddf21a14e 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2191,7 +2191,7 @@ def locate_cell_and_reference_coordinate(self, x, tolerance=None, cell_ignore=No
         if x.size != self.geometric_dimension():
             raise ValueError("Point must have the same geometric dimension as the mesh")
         x = x.reshape((1, self.geometric_dimension()))
-        cells, ref_coords, _ = self.locate_cells_ref_coords_and_dists(x, tolerance=tolerance, cells_ignore=[cell_ignore])
+        cells, ref_coords, _ = self.locate_cells_ref_coords_and_dists(x, tolerance=tolerance, cells_ignore=[[cell_ignore]])
         if cells[0] == -1:
             return None, None
         return cells[0], ref_coords[0]
@@ -2206,9 +2206,10 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non
             this from default will cause the spatial index to be rebuilt which
             can take some time.
         :kwarg cells_ignore: Cell numbers to ignore in the search for each
-            point in xs. Each entry corresponds to a single coordinate in xs.
-            To not ignore any cells, pass None. To ensure a full cell search
-            for any given point, set the corresponding list entry to -1.
+            point in xs. Shape should be (n_ignore_pts, npoints). Each column
+            corresponds to a single coordinate in xs. To not ignore any cells,
+            pass None. To ensure a full cell search for any given point, set
+            the corresponding entries to -1.
         :returns: tuple either
             (cell numbers array, reference coordinates array, ref_cell_dists_l1 array)
             of type
@@ -2229,12 +2230,13 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non
             raise ValueError("Point coordinate dimension does not match mesh geometric dimension")
         Xs = np.empty_like(xs)
         npoints = len(xs)
-        if cells_ignore is None or cells_ignore[0] is None:
-            cells_ignore = np.full(npoints, -1, dtype=IntType)
+        if cells_ignore is None or cells_ignore[0][0] is None:
+            cells_ignore = np.full((1, npoints), -1, dtype=IntType)
         else:
             cells_ignore = np.asarray(cells_ignore, dtype=IntType)
-        if len(cells_ignore) != npoints:
+        if cells_ignore.shape[1] != npoints:
             raise ValueError("Number of cells to ignore does not match number of points")
+        assert cells_ignore.shape == (cells_ignore.shape[0], npoints)
         ref_cell_dists_l1 = np.empty(npoints, dtype=utils.RealType)
         cells = np.empty(npoints, dtype=IntType)
         assert xs.size == npoints * self.geometric_dimension()
@@ -2244,6 +2246,7 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non
                                              ref_cell_dists_l1.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                                              cells.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
                                              npoints,
+                                             cells_ignore.shape[0],
                                              cells_ignore.ctypes.data_as(ctypes.POINTER(ctypes.c_int)))
         return cells, Xs, ref_cell_dists_l1
 
@@ -2259,7 +2262,7 @@ def _c_locator(self, tolerance=None):
         except KeyError:
             src = pq_utils.src_locate_cell(self, tolerance=tolerance)
             src += """
-    int locator(struct Function *f, double *x, double *X, double *ref_cell_dists_l1, int *cells, size_t npoints, int* cells_ignore)
+    int locator(struct Function *f, double *x, double *X, double *ref_cell_dists_l1, int *cells, size_t npoints, size_t ncells_ignore, int* cells_ignore)
     {
         size_t j = 0;  /* index into x and X */
         for(size_t i=0; i<npoints; i++) {
@@ -2272,7 +2275,10 @@ def _c_locator(self, tolerance=None):
             /* to_reference_coords and to_reference_coords_xtr are defined in
             pointquery_utils.py. If they contain python calls, this loop will
             not run at c-loop speed. */
-            cells[i] = locate_cell(f, &x[j], %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i], cells_ignore[i]);
+
+            /* cells_ignore has shape (ncells_ignore, npoints) - find the ith column */
+            int *cells_ignore_i = cells_ignore + i*ncells_ignore;
+            cells[i] = locate_cell(f, &x[j], %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i], ncells_ignore, cells_ignore_i);
 
             for (int k = 0; k < %(geometric_dimension)d; k++) {
                 X[j] = found_reference_coords.X[k];
@@ -2297,6 +2303,7 @@ def _c_locator(self, tolerance=None):
                                 ctypes.POINTER(ctypes.c_double),
                                 ctypes.POINTER(ctypes.c_int),
                                 ctypes.c_size_t,
+                                ctypes.c_size_t,
                                 ctypes.POINTER(ctypes.c_int)]
             locator.restype = ctypes.c_int
             return cache.setdefault(tolerance, locator)
@@ -3748,7 +3755,7 @@ def _parent_mesh_embedding(
     changed_ranks = owned_ranks != ranks
 
     # If distance has changed the the point is not in local mesh partition
-    # since some other cell is closer.
+    # since some other cell on another rank is closer.
     locally_visible[changed_ref_cell_dists_l1] = False
     parent_cell_nums[changed_ref_cell_dists_l1] = -1
     # If the rank has changed but the distance hasn't then there was a tie
@@ -3756,16 +3763,24 @@ def _parent_mesh_embedding(
     # the previously identified cell: if we match the identified owned_rank AND
     # the distance is the same then we have found the correct cell. If we
     # cannot make a match to owned_rank and distance then we can't see the
-    # point. This should only happen for halo points.
+    # point.
     changed_ranks_tied = changed_ranks & ~changed_ref_cell_dists_l1
+
+    # # parent_mesh.comm.rank
+    # # parent_mesh.comm.rank
+    # # locally_visible[9846]
+    # # ranks[9846]
+    # # owned_ranks[9846]
+    # # ref_cell_dists_l1[9846]
+    # # parent_cell_nums[9846]
+    # # for debugging: find index into owned ranks where changed ranks ties equals
+    # idxs = []
+    # for i in range(len(owned_ranks)):
+    #     if (owned_ranks[i] == parent_mesh.comm.rank) & changed_ranks_tied[i]:
+    #         idxs.append(i)
+
     if any(changed_ranks_tied):
-        # check that only halos are affected
-        if any(owned_ranks[changed_ranks_tied] == parent_mesh.comm.rank):
-            raise RuntimeError(
-                "A tie break has occurred in the voting algorithm which has "
-                "resulted in a non-halo point changing MPI rank. This should "
-                "not happen."
-            )
+        old_parent_cell_nums = np.asarray([np.copy(parent_cell_nums)])
         while any(changed_ranks_tied):
             (
                 parent_cell_nums[changed_ranks_tied],
@@ -3774,24 +3789,35 @@ def _parent_mesh_embedding(
             ) = parent_mesh.locate_cells_ref_coords_and_dists(
                 coords_global[changed_ranks_tied],
                 tolerance,
-                cells_ignore=parent_cell_nums[changed_ranks_tied],
+                cells_ignore=old_parent_cell_nums[:, changed_ranks_tied],
             )
             locally_visible[changed_ranks_tied] = (
                 parent_cell_nums[changed_ranks_tied] != -1
             )
             changed_ranks_tied &= locally_visible
             # if new ref_cell_dists_l1 > owned_ref_cell_dists_l1 then we should
-            # disregard the point
+            # disregard the point.
             locally_visible[changed_ranks_tied] &= (
                 ref_cell_dists_l1[changed_ranks_tied]
                 <= owned_ref_cell_dists_l1[changed_ranks_tied]
             )
             changed_ranks_tied &= locally_visible
+            # update the identified rank
+            ranks[changed_ranks_tied] = visible_ranks[
+                parent_cell_nums[changed_ranks_tied]
+            ]
             # if the rank now matches then we have found the correct cell
             locally_visible[changed_ranks_tied] &= (
                 owned_ranks[changed_ranks_tied] == ranks[changed_ranks_tied]
             )
-            changed_ranks_tied &= locally_visible
+            # remove these rank matches from changed_ranks_tied
+            changed_ranks_tied &= ~locally_visible
+
+            # add more cells to ignore
+            old_parent_cell_nums = np.vstack((
+                old_parent_cell_nums,
+                parent_cell_nums)
+            )
 
     # Any ranks which are still np.inf are not in the mesh
     missing_global_idxs = np.where(owned_ranks == np.inf)[0]
diff --git a/firedrake/pointeval_utils.py b/firedrake/pointeval_utils.py
index 61947887e5..954c695a65 100644
--- a/firedrake/pointeval_utils.py
+++ b/firedrake/pointeval_utils.py
@@ -139,7 +139,8 @@ def predicate(index):
     /* The type definitions and arguments used here are defined as statics in pointquery_utils.py */
     double found_ref_cell_dist_l1 = DBL_MAX;
     struct ReferenceCoords temp_reference_coords, found_reference_coords;
-    %(IntType)s cell = locate_cell(f, x, %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &found_ref_cell_dist_l1, -1);
+    int cells_ignore[1] = {-1};
+    %(IntType)s cell = locate_cell(f, x, %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &found_ref_cell_dist_l1, 1, cells_ignore);
     if (cell == -1) {
         return -1;
     }
diff --git a/tests/regression/test_locate_cell.py b/tests/regression/test_locate_cell.py
index f84bf42927..d238d1727f 100644
--- a/tests/regression/test_locate_cell.py
+++ b/tests/regression/test_locate_cell.py
@@ -36,7 +36,7 @@ def value_at(p, cell_ignore=None):
 
     def value_at_and_dist(p, cell_ignore=None):
         if cell_ignore is not None:
-            cell_ignore = [cell_ignore]
+            cell_ignore = [[cell_ignore]]
         cells, _, l1_dists = m.locate_cells_ref_coords_and_dists([p], cells_ignore=cell_ignore)
         return f.dat.data[cells[0]], l1_dists[0]
 

From 5c4c835cbddb7121bedbf5c3bb5f915ebae3e50e Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Fri, 15 Dec 2023 10:12:03 +0000
Subject: [PATCH 07/19] add test

---
 tests/regression/test_locate_cell.py | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/tests/regression/test_locate_cell.py b/tests/regression/test_locate_cell.py
index d238d1727f..ece8dfe013 100644
--- a/tests/regression/test_locate_cell.py
+++ b/tests/regression/test_locate_cell.py
@@ -55,3 +55,24 @@ def test_locate_cell_not_found(meshdata):
     m, f = meshdata
 
     assert m.locate_cell((0.2, -0.4)) is None
+
+
+def test_locate_cells_ref_coords_and_dists(meshdata):
+    m, f = meshdata
+
+    points = [(0.2, 0.1), (0.5, 0.2), (0.7, 0.1), (0.2, 0.4), (0.4, 0.4), (0.8, 0.5), (0.1, 0.7), (0.5, 0.9), (0.9, 0.8)]
+    cells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points)
+    assert np.allclose(f.dat.data[cells], [1, 2, 3, 4, 5, 6, 7, 8, 9])
+    assert np.allclose(l1_dists, 0.0)
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:1], cells[1:2]]).T)
+    assert fcells[0] == -1 or fcells[0] in cells[1:]
+    assert fcells[1] == -1 or fcells[1] in cells[2:] or fcells[1] in cells[:1]
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:2], cells[1:3]]).T)
+    assert fcells[0] == -1 or fcells[0] in cells[2:]
+    assert fcells[1] == -1 or fcells[1] in cells[3:] or fcells[1] in cells[:1]
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:3], cells[1:4]]).T)
+    assert fcells[0] == -1 or fcells[0] in cells[3:]
+    assert fcells[1] == -1 or fcells[1] in cells[4:] or fcells[1] in cells[:1]
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:4], cells[1:5]]).T)
+    assert fcells[0] == -1 or fcells[0] in cells[4:]
+    assert fcells[1] == -1 or fcells[1] in cells[5:] or fcells[1] in cells[:1]

From c2bb720be7eef7950e6d35cf17ab5152d17ce4df Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Fri, 15 Dec 2023 10:59:25 +0000
Subject: [PATCH 08/19] add explicit test

---
 .../regression/test_interpolate_cross_mesh.py | 35 +++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/tests/regression/test_interpolate_cross_mesh.py b/tests/regression/test_interpolate_cross_mesh.py
index ef8343ec61..3ee5c42a2b 100644
--- a/tests/regression/test_interpolate_cross_mesh.py
+++ b/tests/regression/test_interpolate_cross_mesh.py
@@ -750,3 +750,38 @@ def test_missing_dofs_parallel():
 @pytest.mark.parallel
 def test_exact_refinement_parallel():
     test_exact_refinement()
+
+
+def voting_algorithm_edgecases(nprocs):
+    # this triggers lots of cases where the VOM voting algorithm has to deal
+    # with points being claimed by multiple ranks: there are cases where each
+    # rank will claim another one owns a point, for example, and yet also all
+    # claim zero distance to the reference cell!
+    s = nprocs
+    nx = 2 * s
+    mx = 3 * nx
+    mh = [UnitCubeMesh(nx, nx, nx),
+          UnitCubeMesh(mx, mx, mx)]
+    family = "Lagrange"
+    degree = 1
+    Vc = FunctionSpace(mh[0], family, degree=degree)
+    Vf = FunctionSpace(mh[1], family, degree=degree)
+    uc = Function(Vc).interpolate(SpatialCoordinate(mh[0])[0])
+    uf = Function(Vf).interpolate(uc)
+    uf2 = Function(Vf).interpolate(SpatialCoordinate(mh[1])[0])
+    assert np.isclose(errornorm(uf, uf2), 0.0)
+
+
+@pytest.mark.parallel(nprocs=2)
+def test_voting_algorithm_edgecases_2_ranks():
+    voting_algorithm_edgecases(4)
+
+
+@pytest.mark.parallel(nprocs=3)
+def test_voting_algorithm_edgecases_3_ranks():
+    voting_algorithm_edgecases(4)
+
+
+@pytest.mark.parallel(nprocs=4)
+def test_voting_algorithm_edgecases_4_ranks():
+    voting_algorithm_edgecases(4)

From 4969d37c7f2937067471b913d7409669cc6cb69f Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Fri, 15 Dec 2023 11:06:35 +0000
Subject: [PATCH 09/19] wip change cells_ignore ordering

---
 firedrake/mesh.py                    | 18 +++++++++---------
 tests/regression/test_locate_cell.py |  8 ++++----
 2 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 1ddf21a14e..830c24c2d8 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2206,7 +2206,7 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non
             this from default will cause the spatial index to be rebuilt which
             can take some time.
         :kwarg cells_ignore: Cell numbers to ignore in the search for each
-            point in xs. Shape should be (n_ignore_pts, npoints). Each column
+            point in xs. Shape should be (npoints, n_ignore_pts). Each column
             corresponds to a single coordinate in xs. To not ignore any cells,
             pass None. To ensure a full cell search for any given point, set
             the corresponding entries to -1.
@@ -2231,12 +2231,12 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non
         Xs = np.empty_like(xs)
         npoints = len(xs)
         if cells_ignore is None or cells_ignore[0][0] is None:
-            cells_ignore = np.full((1, npoints), -1, dtype=IntType)
+            cells_ignore = np.full((npoints, 1), -1, dtype=IntType, order="C")
         else:
-            cells_ignore = np.asarray(cells_ignore, dtype=IntType)
-        if cells_ignore.shape[1] != npoints:
+            cells_ignore = np.asarray(cells_ignore, dtype=IntType, order="C")
+        if cells_ignore.shape[0] != npoints:
             raise ValueError("Number of cells to ignore does not match number of points")
-        assert cells_ignore.shape == (cells_ignore.shape[0], npoints)
+        assert cells_ignore.shape == (npoints, cells_ignore.shape[1])
         ref_cell_dists_l1 = np.empty(npoints, dtype=utils.RealType)
         cells = np.empty(npoints, dtype=IntType)
         assert xs.size == npoints * self.geometric_dimension()
@@ -2246,8 +2246,8 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non
                                              ref_cell_dists_l1.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                                              cells.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
                                              npoints,
-                                             cells_ignore.shape[0],
-                                             cells_ignore.ctypes.data_as(ctypes.POINTER(ctypes.c_int)))
+                                             cells_ignore.shape[1],
+                                             cells_ignore)
         return cells, Xs, ref_cell_dists_l1
 
     def _c_locator(self, tolerance=None):
@@ -2276,7 +2276,7 @@ def _c_locator(self, tolerance=None):
             pointquery_utils.py. If they contain python calls, this loop will
             not run at c-loop speed. */
 
-            /* cells_ignore has shape (ncells_ignore, npoints) - find the ith column */
+            /* cells_ignore has shape (npoints, ncells_ignore) - find the ith row */
             int *cells_ignore_i = cells_ignore + i*ncells_ignore;
             cells[i] = locate_cell(f, &x[j], %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i], ncells_ignore, cells_ignore_i);
 
@@ -2304,7 +2304,7 @@ def _c_locator(self, tolerance=None):
                                 ctypes.POINTER(ctypes.c_int),
                                 ctypes.c_size_t,
                                 ctypes.c_size_t,
-                                ctypes.POINTER(ctypes.c_int)]
+                                np.ctypeslib.ndpointer(ctypes.c_int, flags="C_CONTIGUOUS")]
             locator.restype = ctypes.c_int
             return cache.setdefault(tolerance, locator)
 
diff --git a/tests/regression/test_locate_cell.py b/tests/regression/test_locate_cell.py
index ece8dfe013..fbbd39e709 100644
--- a/tests/regression/test_locate_cell.py
+++ b/tests/regression/test_locate_cell.py
@@ -64,15 +64,15 @@ def test_locate_cells_ref_coords_and_dists(meshdata):
     cells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points)
     assert np.allclose(f.dat.data[cells], [1, 2, 3, 4, 5, 6, 7, 8, 9])
     assert np.allclose(l1_dists, 0.0)
-    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:1], cells[1:2]]).T)
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:1], cells[1:2]]))
     assert fcells[0] == -1 or fcells[0] in cells[1:]
     assert fcells[1] == -1 or fcells[1] in cells[2:] or fcells[1] in cells[:1]
-    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:2], cells[1:3]]).T)
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:2], cells[1:3]]))
     assert fcells[0] == -1 or fcells[0] in cells[2:]
     assert fcells[1] == -1 or fcells[1] in cells[3:] or fcells[1] in cells[:1]
-    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:3], cells[1:4]]).T)
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:3], cells[1:4]]))
     assert fcells[0] == -1 or fcells[0] in cells[3:]
     assert fcells[1] == -1 or fcells[1] in cells[4:] or fcells[1] in cells[:1]
-    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:4], cells[1:5]]).T)
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:4], cells[1:5]]))
     assert fcells[0] == -1 or fcells[0] in cells[4:]
     assert fcells[1] == -1 or fcells[1] in cells[5:] or fcells[1] in cells[:1]

From c46ed0061404712f76823bb947993086585df855 Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Fri, 15 Dec 2023 13:13:49 +0000
Subject: [PATCH 10/19] Revert "wip change cells_ignore ordering"

This reverts commit 4969d37c7f2937067471b913d7409669cc6cb69f.
---
 firedrake/mesh.py                    | 18 +++++++++---------
 tests/regression/test_locate_cell.py |  8 ++++----
 2 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 830c24c2d8..1ddf21a14e 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2206,7 +2206,7 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non
             this from default will cause the spatial index to be rebuilt which
             can take some time.
         :kwarg cells_ignore: Cell numbers to ignore in the search for each
-            point in xs. Shape should be (npoints, n_ignore_pts). Each column
+            point in xs. Shape should be (n_ignore_pts, npoints). Each column
             corresponds to a single coordinate in xs. To not ignore any cells,
             pass None. To ensure a full cell search for any given point, set
             the corresponding entries to -1.
@@ -2231,12 +2231,12 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non
         Xs = np.empty_like(xs)
         npoints = len(xs)
         if cells_ignore is None or cells_ignore[0][0] is None:
-            cells_ignore = np.full((npoints, 1), -1, dtype=IntType, order="C")
+            cells_ignore = np.full((1, npoints), -1, dtype=IntType)
         else:
-            cells_ignore = np.asarray(cells_ignore, dtype=IntType, order="C")
-        if cells_ignore.shape[0] != npoints:
+            cells_ignore = np.asarray(cells_ignore, dtype=IntType)
+        if cells_ignore.shape[1] != npoints:
             raise ValueError("Number of cells to ignore does not match number of points")
-        assert cells_ignore.shape == (npoints, cells_ignore.shape[1])
+        assert cells_ignore.shape == (cells_ignore.shape[0], npoints)
         ref_cell_dists_l1 = np.empty(npoints, dtype=utils.RealType)
         cells = np.empty(npoints, dtype=IntType)
         assert xs.size == npoints * self.geometric_dimension()
@@ -2246,8 +2246,8 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non
                                              ref_cell_dists_l1.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                                              cells.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
                                              npoints,
-                                             cells_ignore.shape[1],
-                                             cells_ignore)
+                                             cells_ignore.shape[0],
+                                             cells_ignore.ctypes.data_as(ctypes.POINTER(ctypes.c_int)))
         return cells, Xs, ref_cell_dists_l1
 
     def _c_locator(self, tolerance=None):
@@ -2276,7 +2276,7 @@ def _c_locator(self, tolerance=None):
             pointquery_utils.py. If they contain python calls, this loop will
             not run at c-loop speed. */
 
-            /* cells_ignore has shape (npoints, ncells_ignore) - find the ith row */
+            /* cells_ignore has shape (ncells_ignore, npoints) - find the ith column */
             int *cells_ignore_i = cells_ignore + i*ncells_ignore;
             cells[i] = locate_cell(f, &x[j], %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i], ncells_ignore, cells_ignore_i);
 
@@ -2304,7 +2304,7 @@ def _c_locator(self, tolerance=None):
                                 ctypes.POINTER(ctypes.c_int),
                                 ctypes.c_size_t,
                                 ctypes.c_size_t,
-                                np.ctypeslib.ndpointer(ctypes.c_int, flags="C_CONTIGUOUS")]
+                                ctypes.POINTER(ctypes.c_int)]
             locator.restype = ctypes.c_int
             return cache.setdefault(tolerance, locator)
 
diff --git a/tests/regression/test_locate_cell.py b/tests/regression/test_locate_cell.py
index fbbd39e709..ece8dfe013 100644
--- a/tests/regression/test_locate_cell.py
+++ b/tests/regression/test_locate_cell.py
@@ -64,15 +64,15 @@ def test_locate_cells_ref_coords_and_dists(meshdata):
     cells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points)
     assert np.allclose(f.dat.data[cells], [1, 2, 3, 4, 5, 6, 7, 8, 9])
     assert np.allclose(l1_dists, 0.0)
-    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:1], cells[1:2]]))
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:1], cells[1:2]]).T)
     assert fcells[0] == -1 or fcells[0] in cells[1:]
     assert fcells[1] == -1 or fcells[1] in cells[2:] or fcells[1] in cells[:1]
-    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:2], cells[1:3]]))
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:2], cells[1:3]]).T)
     assert fcells[0] == -1 or fcells[0] in cells[2:]
     assert fcells[1] == -1 or fcells[1] in cells[3:] or fcells[1] in cells[:1]
-    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:3], cells[1:4]]))
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:3], cells[1:4]]).T)
     assert fcells[0] == -1 or fcells[0] in cells[3:]
     assert fcells[1] == -1 or fcells[1] in cells[4:] or fcells[1] in cells[:1]
-    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:4], cells[1:5]]))
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:4], cells[1:5]]).T)
     assert fcells[0] == -1 or fcells[0] in cells[4:]
     assert fcells[1] == -1 or fcells[1] in cells[5:] or fcells[1] in cells[:1]

From 5a08ac22251d3ea974fb7bb1958b1faab379e8ac Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Fri, 15 Dec 2023 13:25:53 +0000
Subject: [PATCH 11/19] remember to delete extra dim if needed

---
 firedrake/mesh.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 1ddf21a14e..87c4307a96 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -3784,13 +3784,18 @@ def _parent_mesh_embedding(
         while any(changed_ranks_tied):
             (
                 parent_cell_nums[changed_ranks_tied],
-                reference_coords[changed_ranks_tied, :],
+                new_reference_coords,
                 ref_cell_dists_l1[changed_ranks_tied],
             ) = parent_mesh.locate_cells_ref_coords_and_dists(
                 coords_global[changed_ranks_tied],
                 tolerance,
                 cells_ignore=old_parent_cell_nums[:, changed_ranks_tied],
             )
+            # delete extra dimension if necessary
+            if parent_mesh.geometric_dimension() > parent_mesh.topological_dimension():
+                new_reference_coords = new_reference_coords[:, : parent_mesh.topological_dimension()]
+            reference_coords[changed_ranks_tied, :] = new_reference_coords
+            # remove newly lost points
             locally_visible[changed_ranks_tied] = (
                 parent_cell_nums[changed_ranks_tied] != -1
             )
@@ -3812,7 +3817,6 @@ def _parent_mesh_embedding(
             )
             # remove these rank matches from changed_ranks_tied
             changed_ranks_tied &= ~locally_visible
-
             # add more cells to ignore
             old_parent_cell_nums = np.vstack((
                 old_parent_cell_nums,

From 4e577dd0d14e55382604e7e09a93612ea8d20374 Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Fri, 15 Dec 2023 13:34:20 +0000
Subject: [PATCH 12/19] skip test no longer guaranteed to be true

---
 tests/vertexonly/test_vertex_only_mesh_generation.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/tests/vertexonly/test_vertex_only_mesh_generation.py b/tests/vertexonly/test_vertex_only_mesh_generation.py
index 5ea52d60dc..db87183339 100644
--- a/tests/vertexonly/test_vertex_only_mesh_generation.py
+++ b/tests/vertexonly/test_vertex_only_mesh_generation.py
@@ -180,8 +180,11 @@ def verify_vertexonly_mesh(m, vm, inputvertexcoords, name):
     stored_parent_cell_nums = np.copy(vm.topology_dm.getField("parentcellnum"))
     vm.topology_dm.restoreField("parentcellnum")
     assert len(stored_vertex_coords) == len(stored_parent_cell_nums)
-    for i in range(len(stored_vertex_coords)):
-        assert m.locate_cell(stored_vertex_coords[i]) == stored_parent_cell_nums[i]
+    if MPI.COMM_WORLD.size == 1:
+        for i in range(len(stored_vertex_coords)):
+            # this will only be true if no extra point searches were done,
+            # which is only guaranteed to be true in serial.
+            assert m.locate_cell(stored_vertex_coords[i]) == stored_parent_cell_nums[i]
     # Input is correct (and includes points that were out of bounds)
     vm_input = vm.input_ordering
     assert vm_input.name == name + "_input_ordering"

From 15d3cda981d343389701eac4efd814ad26b445b6 Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Fri, 15 Dec 2023 13:34:51 +0000
Subject: [PATCH 13/19] Revert "Revert "wip change cells_ignore ordering""

This reverts commit c46ed0061404712f76823bb947993086585df855.
---
 firedrake/mesh.py                    | 18 +++++++++---------
 tests/regression/test_locate_cell.py |  8 ++++----
 2 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 87c4307a96..2440f09629 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2206,7 +2206,7 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non
             this from default will cause the spatial index to be rebuilt which
             can take some time.
         :kwarg cells_ignore: Cell numbers to ignore in the search for each
-            point in xs. Shape should be (n_ignore_pts, npoints). Each column
+            point in xs. Shape should be (npoints, n_ignore_pts). Each column
             corresponds to a single coordinate in xs. To not ignore any cells,
             pass None. To ensure a full cell search for any given point, set
             the corresponding entries to -1.
@@ -2231,12 +2231,12 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non
         Xs = np.empty_like(xs)
         npoints = len(xs)
         if cells_ignore is None or cells_ignore[0][0] is None:
-            cells_ignore = np.full((1, npoints), -1, dtype=IntType)
+            cells_ignore = np.full((npoints, 1), -1, dtype=IntType, order="C")
         else:
-            cells_ignore = np.asarray(cells_ignore, dtype=IntType)
-        if cells_ignore.shape[1] != npoints:
+            cells_ignore = np.asarray(cells_ignore, dtype=IntType, order="C")
+        if cells_ignore.shape[0] != npoints:
             raise ValueError("Number of cells to ignore does not match number of points")
-        assert cells_ignore.shape == (cells_ignore.shape[0], npoints)
+        assert cells_ignore.shape == (npoints, cells_ignore.shape[1])
         ref_cell_dists_l1 = np.empty(npoints, dtype=utils.RealType)
         cells = np.empty(npoints, dtype=IntType)
         assert xs.size == npoints * self.geometric_dimension()
@@ -2246,8 +2246,8 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non
                                              ref_cell_dists_l1.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                                              cells.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
                                              npoints,
-                                             cells_ignore.shape[0],
-                                             cells_ignore.ctypes.data_as(ctypes.POINTER(ctypes.c_int)))
+                                             cells_ignore.shape[1],
+                                             cells_ignore)
         return cells, Xs, ref_cell_dists_l1
 
     def _c_locator(self, tolerance=None):
@@ -2276,7 +2276,7 @@ def _c_locator(self, tolerance=None):
             pointquery_utils.py. If they contain python calls, this loop will
             not run at c-loop speed. */
 
-            /* cells_ignore has shape (ncells_ignore, npoints) - find the ith column */
+            /* cells_ignore has shape (npoints, ncells_ignore) - find the ith row */
             int *cells_ignore_i = cells_ignore + i*ncells_ignore;
             cells[i] = locate_cell(f, &x[j], %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i], ncells_ignore, cells_ignore_i);
 
@@ -2304,7 +2304,7 @@ def _c_locator(self, tolerance=None):
                                 ctypes.POINTER(ctypes.c_int),
                                 ctypes.c_size_t,
                                 ctypes.c_size_t,
-                                ctypes.POINTER(ctypes.c_int)]
+                                np.ctypeslib.ndpointer(ctypes.c_int, flags="C_CONTIGUOUS")]
             locator.restype = ctypes.c_int
             return cache.setdefault(tolerance, locator)
 
diff --git a/tests/regression/test_locate_cell.py b/tests/regression/test_locate_cell.py
index ece8dfe013..fbbd39e709 100644
--- a/tests/regression/test_locate_cell.py
+++ b/tests/regression/test_locate_cell.py
@@ -64,15 +64,15 @@ def test_locate_cells_ref_coords_and_dists(meshdata):
     cells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points)
     assert np.allclose(f.dat.data[cells], [1, 2, 3, 4, 5, 6, 7, 8, 9])
     assert np.allclose(l1_dists, 0.0)
-    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:1], cells[1:2]]).T)
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:1], cells[1:2]]))
     assert fcells[0] == -1 or fcells[0] in cells[1:]
     assert fcells[1] == -1 or fcells[1] in cells[2:] or fcells[1] in cells[:1]
-    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:2], cells[1:3]]).T)
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:2], cells[1:3]]))
     assert fcells[0] == -1 or fcells[0] in cells[2:]
     assert fcells[1] == -1 or fcells[1] in cells[3:] or fcells[1] in cells[:1]
-    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:3], cells[1:4]]).T)
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:3], cells[1:4]]))
     assert fcells[0] == -1 or fcells[0] in cells[3:]
     assert fcells[1] == -1 or fcells[1] in cells[4:] or fcells[1] in cells[:1]
-    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:4], cells[1:5]]).T)
+    fcells, ref_coords, l1_dists = m.locate_cells_ref_coords_and_dists(points[:2], cells_ignore=np.array([cells[:4], cells[1:5]]))
     assert fcells[0] == -1 or fcells[0] in cells[4:]
     assert fcells[1] == -1 or fcells[1] in cells[5:] or fcells[1] in cells[:1]

From b39005f2826b5003521530c1c748f6efef4f10f2 Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Fri, 15 Dec 2023 13:48:22 +0000
Subject: [PATCH 14/19] update voting algorithm to use new cell_ignore argument
 order

---
 firedrake/mesh.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 2440f09629..29fbed32e1 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -3789,7 +3789,7 @@ def _parent_mesh_embedding(
             ) = parent_mesh.locate_cells_ref_coords_and_dists(
                 coords_global[changed_ranks_tied],
                 tolerance,
-                cells_ignore=old_parent_cell_nums[:, changed_ranks_tied],
+                cells_ignore=old_parent_cell_nums.T[changed_ranks_tied, :],
             )
             # delete extra dimension if necessary
             if parent_mesh.geometric_dimension() > parent_mesh.topological_dimension():

From e97478f598587f69044d74f57bd6ad90265c3bda Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Fri, 15 Dec 2023 13:48:33 +0000
Subject: [PATCH 15/19] fix test

---
 tests/regression/test_interpolate_cross_mesh.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/regression/test_interpolate_cross_mesh.py b/tests/regression/test_interpolate_cross_mesh.py
index 3ee5c42a2b..4fe57c4d2f 100644
--- a/tests/regression/test_interpolate_cross_mesh.py
+++ b/tests/regression/test_interpolate_cross_mesh.py
@@ -774,12 +774,12 @@ def voting_algorithm_edgecases(nprocs):
 
 @pytest.mark.parallel(nprocs=2)
 def test_voting_algorithm_edgecases_2_ranks():
-    voting_algorithm_edgecases(4)
+    voting_algorithm_edgecases(2)
 
 
 @pytest.mark.parallel(nprocs=3)
 def test_voting_algorithm_edgecases_3_ranks():
-    voting_algorithm_edgecases(4)
+    voting_algorithm_edgecases(3)
 
 
 @pytest.mark.parallel(nprocs=4)

From d3ae063d132e9a14f4a82260af5310a7a6031da3 Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Fri, 15 Dec 2023 13:55:59 +0000
Subject: [PATCH 16/19] tidy

---
 firedrake/mesh.py | 14 --------------
 1 file changed, 14 deletions(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 29fbed32e1..e2ca3c7965 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -3765,20 +3765,6 @@ def _parent_mesh_embedding(
     # cannot make a match to owned_rank and distance then we can't see the
     # point.
     changed_ranks_tied = changed_ranks & ~changed_ref_cell_dists_l1
-
-    # # parent_mesh.comm.rank
-    # # parent_mesh.comm.rank
-    # # locally_visible[9846]
-    # # ranks[9846]
-    # # owned_ranks[9846]
-    # # ref_cell_dists_l1[9846]
-    # # parent_cell_nums[9846]
-    # # for debugging: find index into owned ranks where changed ranks ties equals
-    # idxs = []
-    # for i in range(len(owned_ranks)):
-    #     if (owned_ranks[i] == parent_mesh.comm.rank) & changed_ranks_tied[i]:
-    #         idxs.append(i)
-
     if any(changed_ranks_tied):
         old_parent_cell_nums = np.asarray([np.copy(parent_cell_nums)])
         while any(changed_ranks_tied):

From ed4e262271cafd102b0c5874dd66529bfd488089 Mon Sep 17 00:00:00 2001
From: Reuben Nixon-Hill <ReubenHill@users.noreply.github.com>
Date: Fri, 15 Dec 2023 14:00:46 +0000
Subject: [PATCH 17/19] rename parent_cell_nums cells_ignore_T

---
 firedrake/mesh.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index e2ca3c7965..b8bfd2384a 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -3766,7 +3766,7 @@ def _parent_mesh_embedding(
     # point.
     changed_ranks_tied = changed_ranks & ~changed_ref_cell_dists_l1
     if any(changed_ranks_tied):
-        old_parent_cell_nums = np.asarray([np.copy(parent_cell_nums)])
+        cells_ignore_T = np.asarray([np.copy(parent_cell_nums)])
         while any(changed_ranks_tied):
             (
                 parent_cell_nums[changed_ranks_tied],
@@ -3775,7 +3775,7 @@ def _parent_mesh_embedding(
             ) = parent_mesh.locate_cells_ref_coords_and_dists(
                 coords_global[changed_ranks_tied],
                 tolerance,
-                cells_ignore=old_parent_cell_nums.T[changed_ranks_tied, :],
+                cells_ignore=cells_ignore_T.T[changed_ranks_tied, :],
             )
             # delete extra dimension if necessary
             if parent_mesh.geometric_dimension() > parent_mesh.topological_dimension():
@@ -3804,8 +3804,8 @@ def _parent_mesh_embedding(
             # remove these rank matches from changed_ranks_tied
             changed_ranks_tied &= ~locally_visible
             # add more cells to ignore
-            old_parent_cell_nums = np.vstack((
-                old_parent_cell_nums,
+            cells_ignore_T = np.vstack((
+                cells_ignore_T,
                 parent_cell_nums)
             )
 

From cb7939a0c008fedfe11b911e44d83da937560649 Mon Sep 17 00:00:00 2001
From: Connor <c.ward20@imperial.ac.uk>
Date: Wed, 23 Oct 2024 15:55:01 +0100
Subject: [PATCH 18/19] linting

---
 tests/regression/test_interpolate_cross_mesh.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/regression/test_interpolate_cross_mesh.py b/tests/regression/test_interpolate_cross_mesh.py
index a6e98dee55..3f53287885 100644
--- a/tests/regression/test_interpolate_cross_mesh.py
+++ b/tests/regression/test_interpolate_cross_mesh.py
@@ -790,6 +790,7 @@ def test_voting_algorithm_edgecases_3_ranks():
 def test_voting_algorithm_edgecases_4_ranks():
     voting_algorithm_edgecases(4)
 
+
 @pytest.mark.parallel
 @pytest.mark.parametrize('periodic', [False, True])
 def test_interpolate_cross_mesh_interval(periodic):

From 68a42e10bb1130937455f9d054c151cb1b53267a Mon Sep 17 00:00:00 2001
From: Connor <c.ward20@imperial.ac.uk>
Date: Wed, 23 Oct 2024 16:27:44 +0100
Subject: [PATCH 19/19] fstring fixup

---
 firedrake/mesh.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 53a47e2314..a844889274 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2628,7 +2628,7 @@ def _c_locator(self, tolerance=None):
                         not run at c-loop speed. */
                         /* cells_ignore has shape (npoints, ncells_ignore) - find the ith row */
                         int *cells_ignore_i = cells_ignore + i*ncells_ignore;
-                        cells[i] = locate_cell(f, &x[j], %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i], ncells_ignore, cells_ignore_i);
+                        cells[i] = locate_cell(f, &x[j], {self.geometric_dimension()}, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i], ncells_ignore, cells_ignore_i);
 
                         for (int k = 0; k < {self.geometric_dimension()}; k++) {{
                             X[j] = found_reference_coords.X[k];