Fix voting algorithm

firedrake/evaluate.h

@@ -52,7 +52,9 @@ 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,
+			   size_t ncells_ignore,
+			   int* cells_ignore);
 
 extern int evaluate(struct Function *f,
 		    double *x,

firedrake/locate.c

@@ -11,10 +11,13 @@ 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,
+        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
@@ -42,6 +45,16 @@ 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);
+                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) {
                     /* Found cell! */
                     cell = ids[i];
@@ -67,6 +80,16 @@ 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);
+                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) {
                     /* Found cell! */
                     cell = ids[i];
@@ -91,6 +114,16 @@ 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);
+                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) {
                     /* Found cell! */
                     cell = c;
@@ -114,6 +147,16 @@ 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);
+                    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) {
                         /* Found cell! */
                         cell = l;

firedrake/mesh.py

@@ -2494,20 +2494,21 @@ 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
         :kwarg tolerance: Tolerance for checking if a point is in a cell.
             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.
 
@@ -2516,12 +2517,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.
 
@@ -2530,6 +2532,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).
@@ -2538,12 +2541,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.
 
@@ -2552,6 +2555,11 @@ 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. 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.
         :returns: tuple either
             (cell numbers array, reference coordinates array, ref_cell_dists_l1 array)
             of type
@@ -2572,6 +2580,13 @@ 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][0] is None:
+            cells_ignore = np.full((npoints, 1), -1, dtype=IntType, order="C")
+        else:
+            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 == (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()
@@ -2580,7 +2595,9 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None):
                                              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,
+                                             cells_ignore.shape[1],
+                                             cells_ignore)
         return cells, Xs, ref_cell_dists_l1
 
     def _c_locator(self, tolerance=None):
@@ -2596,7 +2613,7 @@ def _c_locator(self, tolerance=None):
             IntTypeC = as_cstr(IntType)
             src = pq_utils.src_locate_cell(self, tolerance=tolerance)
             src += dedent(f"""
-                int locator(struct Function *f, double *x, double *X, double *ref_cell_dists_l1, {IntTypeC} *cells, {IntTypeC} npoints)
+                int locator(struct Function *f, double *x, double *X, double *ref_cell_dists_l1, {IntTypeC} *cells, {IntTypeC} npoints, size_t ncells_ignore, int* cells_ignore)
                 {{
                     {IntTypeC} j = 0;  /* index into x and X */
                     for({IntTypeC} i=0; i<npoints; i++) {{
@@ -2609,7 +2626,9 @@ 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], {self.geometric_dimension()}, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i]);
+                        /* 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], {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];
@@ -2643,7 +2662,9 @@ 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_size_t,
+                                np.ctypeslib.ndpointer(ctypes.c_int, flags="C_CONTIGUOUS")]
             locator.restype = ctypes.c_int
             return cache.setdefault(tolerance, locator)
 
@@ -4100,12 +4121,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
@@ -4124,23 +4151,73 @@ 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."
-        )
-    locally_visible[extra_missing_points] = False
-    parent_cell_nums[extra_missing_points] = -1
+    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 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
+    # 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.
+    changed_ranks_tied = changed_ranks & ~changed_ref_cell_dists_l1
+    if any(changed_ranks_tied):
+        cells_ignore_T = np.asarray([np.copy(parent_cell_nums)])
+        while any(changed_ranks_tied):
+            (
+                parent_cell_nums[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=cells_ignore_T.T[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
+            )
+            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] &= (
+                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]
+            )
+            # remove these rank matches from changed_ranks_tied
+            changed_ranks_tied &= ~locally_visible
+            # add more cells to ignore
+            cells_ignore_T = np.vstack((
+                cells_ignore_T,
+                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]

firedrake/pointeval_utils.py

@@ -140,7 +140,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);
+    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;
     }

tests/regression/test_interpolate_cross_mesh.py

@@ -756,6 +756,41 @@ 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(2)
+
+
+@pytest.mark.parallel(nprocs=3)
+def test_voting_algorithm_edgecases_3_ranks():
+    voting_algorithm_edgecases(3)
+
+
+@pytest.mark.parallel(nprocs=4)
+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):