diff --git a/firedrake/evaluate.h b/firedrake/evaluate.h index 4253c2174e..738a573867 100644 --- a/firedrake/evaluate.h +++ b/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, diff --git a/firedrake/locate.c b/firedrake/locate.c index 2588fa3ef3..a243e2119f 100644 --- a/firedrake/locate.c +++ b/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; diff --git a/firedrake/mesh.py b/firedrake/mesh.py index acc9031542..a844889274 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2494,7 +2494,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 @@ -2502,12 +2502,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. @@ -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 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] diff --git a/firedrake/pointeval_utils.py b/firedrake/pointeval_utils.py index ab2c3b0193..bf858f9c72 100644 --- a/firedrake/pointeval_utils.py +++ b/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; } diff --git a/tests/regression/test_interpolate_cross_mesh.py b/tests/regression/test_interpolate_cross_mesh.py index 1bd6ed6226..3f53287885 100644 --- a/tests/regression/test_interpolate_cross_mesh.py +++ b/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): diff --git a/tests/regression/test_locate_cell.py b/tests/regression/test_locate_cell.py index 21b35a7c66..fbbd39e709 100644 --- a/tests/regression/test_locate_cell.py +++ b/tests/regression/test_locate_cell.py @@ -30,14 +30,49 @@ 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): 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]])) + 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]])) + 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]])) + 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]])) + 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] diff --git a/tests/vertexonly/test_vertex_only_mesh_generation.py b/tests/vertexonly/test_vertex_only_mesh_generation.py index 3ac5bb43ba..559d46fef6 100644 --- a/tests/vertexonly/test_vertex_only_mesh_generation.py +++ b/tests/vertexonly/test_vertex_only_mesh_generation.py @@ -202,8 +202,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"