53 std::span<const std::int32_t> entities,
54 std::span<const T> points)
56 const int tdim =
mesh.topology()->dim();
59 std::span<const T> geom_dofs =
geometry.x();
61 std::vector<T> shortest_vectors;
62 shortest_vectors.reserve(3 * entities.size());
65 for (std::size_t e = 0; e < entities.size(); e++)
69 assert(entities[e] >= 0);
70 auto dofs = md::submdspan(x_dofmap, entities[e], md::full_extent);
71 std::vector<T> nodes(3 * dofs.size());
72 for (std::size_t i = 0; i < dofs.size(); ++i)
74 const std::int32_t pos = 3 * dofs[i];
75 for (std::size_t j = 0; j < 3; ++j)
76 nodes[3 * i + j] = geom_dofs[pos + j];
81 shortest_vectors.insert(shortest_vectors.end(), d.begin(), d.end());
86 mesh.topology_mutable()->create_connectivity(dim, tdim);
87 mesh.topology_mutable()->create_connectivity(tdim, dim);
88 auto e_to_c =
mesh.topology()->connectivity(dim, tdim);
90 auto c_to_e =
mesh.topology_mutable()->connectivity(tdim, dim);
92 for (std::size_t e = 0; e < entities.size(); e++)
94 const std::int32_t index = entities[e];
97 assert(e_to_c->num_links(index) > 0);
98 const std::int32_t c = e_to_c->links(index)[0];
101 auto cell_entities = c_to_e->links(c);
102 auto it0 = std::find(cell_entities.begin(), cell_entities.end(), index);
103 assert(it0 != cell_entities.end());
104 const int local_cell_entity = std::distance(cell_entities.begin(), it0);
107 auto dofs = md::submdspan(x_dofmap, c, md::full_extent);
108 const std::vector<int> entity_dofs
109 =
geometry.cmap().create_dof_layout().entity_closure_dofs(
110 dim, local_cell_entity);
111 std::vector<T> nodes(3 * entity_dofs.size());
112 for (std::size_t i = 0; i < entity_dofs.size(); i++)
114 const std::int32_t pos = 3 * dofs[entity_dofs[i]];
115 for (std::size_t j = 0; j < 3; ++j)
116 nodes[3 * i + j] = geom_dofs[pos + j];
121 shortest_vectors.insert(shortest_vectors.end(), d.begin(), d.end());
125 return shortest_vectors;
687 std::span<const T> points,
690 MPI_Comm comm =
mesh.comm();
694 const int tdim =
mesh.topology()->dim();
695 auto cell_map =
mesh.topology()->index_map(tdim);
696 const std::int32_t num_cells = cell_map->size_local();
698 std::vector<std::int32_t> cells(num_cells, 0);
699 std::iota(cells.begin(), cells.end(), 0);
708 std::vector<std::int32_t> out_ranks = collisions.
array();
709 std::ranges::sort(out_ranks);
710 auto [unique_end, range_end] = std::ranges::unique(out_ranks);
711 out_ranks.erase(unique_end, range_end);
715 std::ranges::sort(in_ranks);
718 MPI_Comm forward_comm;
719 MPI_Dist_graph_create_adjacent(
720 comm, in_ranks.size(), in_ranks.data(), MPI_UNWEIGHTED, out_ranks.size(),
721 out_ranks.data(), MPI_UNWEIGHTED, MPI_INFO_NULL,
false, &forward_comm);
725 std::map<std::int32_t, std::int32_t> rank_to_neighbor;
726 for (std::size_t i = 0; i < out_ranks.size(); i++)
727 rank_to_neighbor[out_ranks[i]] = i;
730 std::vector<std::int32_t> send_sizes(out_ranks.size());
731 for (std::size_t i = 0; i < points.size() / 3; ++i)
732 for (
auto p : collisions.
links(i))
733 send_sizes[rank_to_neighbor[p]] += 3;
736 std::vector<std::int32_t> recv_sizes(in_ranks.size());
737 send_sizes.reserve(1);
738 recv_sizes.reserve(1);
739 MPI_Request sizes_request;
740 MPI_Ineighbor_alltoall(send_sizes.data(), 1, MPI_INT, recv_sizes.data(), 1,
741 MPI_INT, forward_comm, &sizes_request);
744 std::vector<std::int32_t> send_offsets(send_sizes.size() + 1, 0);
745 std::partial_sum(send_sizes.begin(), send_sizes.end(),
746 std::next(send_offsets.begin(), 1));
749 std::vector<T> send_data(send_offsets.back());
750 std::vector<std::int32_t> counter(send_sizes.size(), 0);
752 std::vector<std::int32_t> unpack_map(send_offsets.back() / 3);
753 for (std::size_t i = 0; i < points.size(); i += 3)
755 for (
auto p : collisions.
links(i / 3))
757 int neighbor = rank_to_neighbor[p];
758 int pos = send_offsets[neighbor] + counter[neighbor];
759 auto it = std::next(send_data.begin(), pos);
760 std::copy_n(std::next(points.begin(), i), 3, it);
761 unpack_map[pos / 3] = i / 3;
762 counter[neighbor] += 3;
766 MPI_Wait(&sizes_request, MPI_STATUS_IGNORE);
767 std::vector<std::int32_t> recv_offsets(in_ranks.size() + 1, 0);
768 std::partial_sum(recv_sizes.begin(), recv_sizes.end(),
769 std::next(recv_offsets.begin(), 1));
771 std::vector<T> received_points((std::size_t)recv_offsets.back());
772 MPI_Neighbor_alltoallv(
773 send_data.data(), send_sizes.data(), send_offsets.data(),
779 std::span<const T> geom_dofs =
geometry.x();
785 received_points.size()));
789 std::vector<std::int32_t> cell_indicator(received_points.size() / 3);
790 std::vector<std::int32_t> closest_cells(received_points.size() / 3);
791 for (std::size_t p = 0; p < received_points.size(); p += 3)
793 std::array<T, 3> point;
794 std::copy_n(std::next(received_points.begin(), p), 3, point.begin());
797 mesh, candidate_collisions.
links(p / 3), point,
798 10 * std::numeric_limits<T>::epsilon());
801 cell_indicator[p / 3] = (colliding_cell >= 0) ? rank : -1;
804 closest_cells[p / 3] = colliding_cell;
809 MPI_Comm reverse_comm;
810 MPI_Dist_graph_create_adjacent(
811 comm, out_ranks.size(), out_ranks.data(), MPI_UNWEIGHTED, in_ranks.size(),
812 in_ranks.data(), MPI_UNWEIGHTED, MPI_INFO_NULL,
false, &reverse_comm);
817 auto rescale = [](
auto& x)
818 { std::ranges::transform(x, x.begin(), [](
auto e) { return (e / 3); }); };
820 rescale(recv_offsets);
822 rescale(send_offsets);
825 std::swap(recv_sizes, send_sizes);
826 std::swap(recv_offsets, send_offsets);
829 std::vector<std::int32_t> recv_ranks(recv_offsets.back());
830 MPI_Neighbor_alltoallv(cell_indicator.data(), send_sizes.data(),
831 send_offsets.data(), MPI_INT32_T, recv_ranks.data(),
832 recv_sizes.data(), recv_offsets.data(), MPI_INT32_T,
835 std::vector<int> point_owners(points.size() / 3, -1);
836 for (std::size_t i = 0; i < unpack_map.size(); i++)
838 const std::int32_t pos = unpack_map[i];
840 if (recv_ranks[i] >= 0 && point_owners[pos] == -1)
841 point_owners[pos] = recv_ranks[i];
846 std::vector<std::uint8_t> send_extrapolate(recv_offsets.back());
847 for (std::int32_t i = 0; i < recv_offsets.back(); i++)
849 const std::int32_t pos = unpack_map[i];
850 send_extrapolate[i] = point_owners[pos] == -1;
855 std::swap(send_sizes, recv_sizes);
856 std::swap(send_offsets, recv_offsets);
857 std::vector<std::uint8_t> dest_extrapolate(recv_offsets.back());
858 MPI_Neighbor_alltoallv(send_extrapolate.data(), send_sizes.data(),
859 send_offsets.data(), MPI_UINT8_T,
860 dest_extrapolate.data(), recv_sizes.data(),
861 recv_offsets.data(), MPI_UINT8_T, forward_comm);
863 std::vector<T> squared_distances(received_points.size() / 3, -1);
865 for (std::size_t i = 0; i < dest_extrapolate.size(); i++)
867 if (dest_extrapolate[i] == 1)
869 assert(closest_cells[i] == -1);
870 std::array<T, 3> point;
871 std::copy_n(std::next(received_points.begin(), 3 * i), 3, point.begin());
874 T shortest_distance = std::numeric_limits<T>::max();
875 std::int32_t closest_cell = -1;
876 for (
auto cell : candidate_collisions.
links(i))
878 auto dofs = md::submdspan(x_dofmap, cell, md::full_extent);
879 std::vector<T> nodes(3 * dofs.size());
880 for (std::size_t j = 0; j < dofs.size(); ++j)
882 const int pos = 3 * dofs[j];
883 for (std::size_t k = 0; k < 3; ++k)
884 nodes[3 * j + k] = geom_dofs[pos + k];
887 std::span<const T>(point.data(), point.size()), nodes);
888 if (T current_distance = d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
889 current_distance < shortest_distance)
891 shortest_distance = current_distance;
895 closest_cells[i] = closest_cell;
896 squared_distances[i] = shortest_distance;
900 std::swap(recv_sizes, send_sizes);
901 std::swap(recv_offsets, send_offsets);
904 std::vector<T> recv_distances(recv_offsets.back());
905 MPI_Neighbor_alltoallv(
906 squared_distances.data(), send_sizes.data(), send_offsets.data(),
911 std::vector<T> closest_distance(point_owners.size(),
912 std::numeric_limits<T>::max());
913 for (std::size_t i = 0; i < out_ranks.size(); i++)
915 for (std::int32_t j = recv_offsets[i]; j < recv_offsets[i + 1]; j++)
917 const std::int32_t pos = unpack_map[j];
918 auto current_dist = recv_distances[j];
920 if (
auto d = closest_distance[pos];
921 (current_dist > 0) and (current_dist < d))
923 point_owners[pos] = out_ranks[i];
924 closest_distance[pos] = current_dist;
930 std::swap(send_sizes, recv_sizes);
931 std::swap(send_offsets, recv_offsets);
934 std::vector<std::int32_t> send_owners(send_offsets.back());
935 std::ranges::fill(counter, 0);
936 for (std::size_t i = 0; i < points.size() / 3; ++i)
938 for (
auto p : collisions.
links(i))
940 int neighbor = rank_to_neighbor[p];
941 send_owners[send_offsets[neighbor] + counter[neighbor]++]
947 std::vector<std::int32_t> dest_ranks(recv_offsets.back());
948 MPI_Neighbor_alltoallv(send_owners.data(), send_sizes.data(),
949 send_offsets.data(), MPI_INT32_T, dest_ranks.data(),
950 recv_sizes.data(), recv_offsets.data(), MPI_INT32_T,
954 std::vector<int> owned_recv_ranks;
955 owned_recv_ranks.reserve(recv_offsets.back());
956 std::vector<T> owned_recv_points;
957 std::vector<std::int32_t> owned_recv_cells;
958 for (std::size_t i = 0; i < in_ranks.size(); i++)
960 for (std::int32_t j = recv_offsets[i]; j < recv_offsets[i + 1]; j++)
962 if (rank == dest_ranks[j])
964 owned_recv_ranks.push_back(in_ranks[i]);
965 owned_recv_points.insert(
966 owned_recv_points.end(), std::next(received_points.cbegin(), 3 * j),
967 std::next(received_points.cbegin(), 3 * (j + 1)));
968 owned_recv_cells.push_back(closest_cells[j]);
973 MPI_Comm_free(&forward_comm);
974 MPI_Comm_free(&reverse_comm);
976 .dest_owners = std::move(owned_recv_ranks),
977 .dest_points = std::move(owned_recv_points),
978 .dest_cells = std::move(owned_recv_cells)};