DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
utils.h
1// Copyright (C) 2019-2021 Garth N. Wells and Jørgen S. Dokken
2//
3// This file is part of DOLFINx (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#pragma once
8
9#include "BoundingBoxTree.h"
10#include "gjk.h"
11#include <algorithm>
12#include <array>
13#include <concepts>
14#include <cstdint>
15#include <deque>
16#include <dolfinx/graph/AdjacencyList.h>
17#include <dolfinx/mesh/Mesh.h>
18#include <map>
19#include <numeric>
20#include <span>
21#include <vector>
22
23namespace dolfinx::geometry
24{
28template <std::floating_point T>
30{
31 std::vector<int> src_owner;
33 std::vector<int>
35 std::vector<T> dest_points;
36 std::vector<std::int32_t>
39};
40
51template <std::floating_point T>
52std::vector<T> shortest_vector(const mesh::Mesh<T>& mesh, int dim,
53 std::span<const std::int32_t> entities,
54 std::span<const T> points)
55{
56 const int tdim = mesh.topology()->dim();
57 const mesh::Geometry<T>& geometry = mesh.geometry();
58
59 std::span<const T> geom_dofs = geometry.x();
60 auto x_dofmap = geometry.dofmap();
61 std::vector<T> shortest_vectors;
62 shortest_vectors.reserve(3 * entities.size());
63 if (dim == tdim)
64 {
65 for (std::size_t e = 0; e < entities.size(); e++)
66 {
67 // Check that we have sent in valid entities, i.e. that they exist in the
68 // local dofmap. One gets a cryptical memory segfault if entities is -1
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)
73 {
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];
77 }
78
79 std::array<T, 3> d
80 = compute_distance_gjk<T>(points.subspan(3 * e, 3), nodes);
81 shortest_vectors.insert(shortest_vectors.end(), d.begin(), d.end());
82 }
83 }
84 else
85 {
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);
89 assert(e_to_c);
90 auto c_to_e = mesh.topology_mutable()->connectivity(tdim, dim);
91 assert(c_to_e);
92 for (std::size_t e = 0; e < entities.size(); e++)
93 {
94 const std::int32_t index = entities[e];
95
96 // Find attached cell
97 assert(e_to_c->num_links(index) > 0);
98 const std::int32_t c = e_to_c->links(index)[0];
99
100 // Find local number of entity wrt cell
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);
105
106 // Tabulate geometry dofs for the entity
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++)
113 {
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];
117 }
118
119 std::array<T, 3> d
120 = compute_distance_gjk<T>(points.subspan(3 * e, 3), nodes);
121 shortest_vectors.insert(shortest_vectors.end(), d.begin(), d.end());
122 }
123 }
124
125 return shortest_vectors;
126}
127
134template <std::floating_point T>
135T compute_squared_distance_bbox(std::span<const T, 6> b,
136 std::span<const T, 3> x)
137{
138 auto b0 = b.template subspan<0, 3>();
139 auto b1 = b.template subspan<3, 3>();
140 return std::transform_reduce(x.begin(), x.end(), b0.begin(), 0.0,
141 std::plus<>{},
142 [](auto x, auto b)
143 {
144 auto dx = x - b;
145 return dx > 0 ? 0 : dx * dx;
146 })
147 + std::transform_reduce(x.begin(), x.end(), b1.begin(), 0.0,
148 std::plus<>{},
149 [](auto x, auto b)
150 {
151 auto dx = x - b;
152 return dx < 0 ? 0 : dx * dx;
153 });
154}
155
171template <std::floating_point T>
172std::vector<T> squared_distance(const mesh::Mesh<T>& mesh, int dim,
173 std::span<const std::int32_t> entities,
174 std::span<const T> points)
175{
176 std::vector<T> v = shortest_vector(mesh, dim, entities, points);
177 std::vector<T> d(v.size() / 3, 0);
178 for (std::size_t i = 0; i < d.size(); ++i)
179 for (std::size_t j = 0; j < 3; ++j)
180 d[i] += v[3 * i + j] * v[3 * i + j];
181 return d;
182}
183
184namespace impl
185{
187constexpr bool is_leaf(std::array<int, 2> bbox)
188{
189 // Leaf nodes are marked by setting child_0 equal to child_1
190 return bbox[0] == bbox[1];
191}
192
197template <std::floating_point T>
198constexpr bool point_in_bbox(std::span<const T, 6> b, std::span<const T, 3> x)
199{
200 constexpr T rtol = 1e-14;
201 bool in = true;
202 for (std::size_t i = 0; i < 3; i++)
203 {
204 T eps = rtol * (b[i + 3] - b[i]);
205 in &= (x[i] >= (b[i] - eps)) && (x[i] <= (b[i + 3] + eps));
206 if (!in)
207 break;
208 }
209
210 return in;
211}
212
216template <std::floating_point T>
217constexpr bool bbox_in_bbox(std::span<const T, 6> a, std::span<const T, 6> b)
218{
219 constexpr T rtol = 1e-14;
220 auto a0 = a.template subspan<0, 3>();
221 auto a1 = a.template subspan<3, 3>();
222 auto b0 = b.template subspan<0, 3>();
223 auto b1 = b.template subspan<3, 3>();
224
225 bool in = true;
226 for (std::size_t i = 0; i < 3; i++)
227 {
228 T eps = rtol * (b1[i] - b0[i]);
229 in &= a1[i] >= (b0[i] - eps);
230 in &= a0[i] <= (b1[i] + eps);
231 }
232
233 return in;
234}
235
237template <std::floating_point T>
238std::pair<std::int32_t, T>
239_compute_closest_entity(const geometry::BoundingBoxTree<T>& tree,
240 std::span<const T, 3> point, std::int32_t node,
241 const mesh::Mesh<T>& mesh, std::int32_t closest_entity,
242 T R2)
243{
244 // Get children of current bounding box node (child_1 denotes entity
245 // index for leaves)
246 const std::array<int, 2> bbox = tree.bbox(node);
247 T r2;
248 if (is_leaf(bbox))
249 {
250 // If point cloud tree the exact distance is easy to compute
251 if (tree.tdim() == 0)
252 {
253 std::array<T, 6> diff = tree.get_bbox(node);
254 for (std::size_t k = 0; k < 3; ++k)
255 diff[k] -= point[k];
256 r2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2];
257 }
258 else
259 {
260 r2 = compute_squared_distance_bbox<T>(tree.get_bbox(node), point);
261
262 // If bounding box closer than previous closest entity, use gjk to
263 // obtain exact distance to the convex hull of the entity
264 if (r2 <= R2)
265 {
266 r2 = squared_distance<T>(mesh, tree.tdim(),
267 std::span(std::next(bbox.begin(), 1), 1),
268 point)
269 .front();
270 }
271 }
272
273 // If entity is closer than best result so far, return it
274 if (r2 <= R2)
275 {
276 closest_entity = bbox.back();
277 R2 = r2;
278 }
279
280 return {closest_entity, R2};
281 }
282 else
283 {
284 // If bounding box is outside radius, then don't search further
285 r2 = compute_squared_distance_bbox<T>(tree.get_bbox(node), point);
286 if (r2 > R2)
287 return {closest_entity, R2};
288
289 // Check both children. We use R2 (as opposed to r2), as a bounding
290 // box can be closer than the actual entity.
291 std::pair<int, T> p0 = _compute_closest_entity(tree, point, bbox.front(),
292 mesh, closest_entity, R2);
293 std::pair<int, T> p1 = _compute_closest_entity(tree, point, bbox.back(),
294 mesh, p0.first, p0.second);
295 return p1;
296 }
297}
298
304template <std::floating_point T>
305void _compute_collisions_point(const geometry::BoundingBoxTree<T>& tree,
306 std::span<const T, 3> p,
307 std::vector<std::int32_t>& entities)
308{
309 std::deque<std::int32_t> stack;
310 std::int32_t next = tree.num_bboxes() - 1;
311 std::span<const T> coords = tree.bbox_coordinates();
312 auto view_bbox = [&coords](std::size_t node)
313 { return std::span<const T, 6>(coords.data() + 6 * node, 6); };
314 while (next != -1)
315 {
316 if (std::array bbox = tree.bbox(next);
317 is_leaf(bbox) and point_in_bbox(view_bbox(next), p))
318 {
319 // If box is a leaf node then add it to the list of colliding
320 // entities
321 entities.push_back(bbox[1]);
322 next = -1;
323 }
324 else
325 {
326 // Check whether the point collides with child nodes (left and
327 // right)
328 bool left = point_in_bbox(view_bbox(bbox[0]), p);
329 bool right = point_in_bbox(view_bbox(bbox[1]), p);
330 if (left and right)
331 {
332 // If the point collides with both child nodes, add the right
333 // node to the stack (for later visiting) and continue the tree
334 // traversal with the left subtree
335 stack.push_back(bbox[1]);
336 next = bbox[0];
337 }
338 else if (left)
339 {
340 // Traverse the current node's left subtree
341 next = bbox[0];
342 }
343 else if (right)
344 {
345 // Traverse the current node's right subtree
346 next = bbox[1];
347 }
348 else
349 next = -1;
350 }
351
352 // If tree traversal reaches a dead end (box is a leaf node or no
353 // collision detected), check the stack for deferred subtrees
354 if (next == -1 and !stack.empty())
355 {
356 next = stack.back();
357 stack.pop_back();
358 }
359 }
360}
361
362// Compute collisions with tree (recursive)
363template <std::floating_point T>
364void _compute_collisions_tree(const geometry::BoundingBoxTree<T>& A,
365 const geometry::BoundingBoxTree<T>& B,
366 std::int32_t node_A, std::int32_t node_B,
367 std::vector<std::int32_t>& entities)
368{
369 // If bounding boxes don't collide, then don't search further
370 if (!bbox_in_bbox<T>(A.get_bbox(node_A), B.get_bbox(node_B)))
371 return;
372
373 // Get bounding boxes for current nodes
374 const std::array<std::int32_t, 2> bbox_A = A.bbox(node_A);
375 const std::array<std::int32_t, 2> bbox_B = B.bbox(node_B);
376
377 // Check whether we've reached a leaf in A or B
378 const bool is_leaf_A = is_leaf(bbox_A);
379 const bool is_leaf_B = is_leaf(bbox_B);
380 if (is_leaf_A and is_leaf_B)
381 {
382 // If both boxes are leaves (which we know collide), then add them
383 // child_1 denotes entity for leaves
384 entities.push_back(bbox_A[1]);
385 entities.push_back(bbox_B[1]);
386 }
387 else if (is_leaf_A)
388 {
389 // If we reached the leaf in A, then descend B
390 _compute_collisions_tree(A, B, node_A, bbox_B[0], entities);
391 _compute_collisions_tree(A, B, node_A, bbox_B[1], entities);
392 }
393 else if (is_leaf_B)
394 {
395 // If we reached the leaf in B, then descend A
396 _compute_collisions_tree(A, B, bbox_A[0], node_B, entities);
397 _compute_collisions_tree(A, B, bbox_A[1], node_B, entities);
398 }
399 else if (node_A > node_B)
400 {
401 // At this point, we know neither is a leaf so descend the largest
402 // tree first. Note that nodes are added in reverse order with the
403 // top bounding box at the end so the largest tree (the one with the
404 // the most boxes left to traverse) has the largest node number.
405 _compute_collisions_tree(A, B, bbox_A[0], node_B, entities);
406 _compute_collisions_tree(A, B, bbox_A[1], node_B, entities);
407 }
408 else
409 {
410 _compute_collisions_tree(A, B, node_A, bbox_B[0], entities);
411 _compute_collisions_tree(A, B, node_A, bbox_B[1], entities);
412 }
413
414 // Note that cases above can be collected in fewer cases but this way
415 // the logic is easier to follow.
416}
417
418} // namespace impl
419
426template <std::floating_point T>
428 std::span<const std::int32_t> entities)
429{
430 spdlog::info("Building point search tree to accelerate distance queries for "
431 "a given topological dimension and subset of entities.");
432
433 const std::vector<T> midpoints
434 = mesh::compute_midpoints(mesh, tdim, entities);
435 std::vector<std::pair<std::array<T, 3>, std::int32_t>> points(
436 entities.size());
437 for (std::size_t i = 0; i < points.size(); ++i)
438 {
439 for (std::size_t j = 0; j < 3; ++j)
440 points[i].first[j] = midpoints[3 * i + j];
441 points[i].second = entities[i];
442 }
443
444 // Build tree
445 return BoundingBoxTree(points);
446}
447
453template <std::floating_point T>
454std::vector<std::int32_t> compute_collisions(const BoundingBoxTree<T>& tree0,
455 const BoundingBoxTree<T>& tree1)
456{
457 // Call recursive find function
458 std::vector<std::int32_t> entities;
459 if (tree0.num_bboxes() > 0 and tree1.num_bboxes() > 0)
460 {
461 impl::_compute_collisions_tree(tree0, tree1, tree0.num_bboxes() - 1,
462 tree1.num_bboxes() - 1, entities);
463 }
464
465 return entities;
466}
467
478template <std::floating_point T>
480compute_collisions(const BoundingBoxTree<T>& tree, std::span<const T> points)
481{
482 if (tree.num_bboxes() > 0)
483 {
484 std::vector<std::int32_t> entities, offsets(points.size() / 3 + 1, 0);
485 entities.reserve(points.size() / 3);
486 for (std::size_t p = 0; p < points.size() / 3; ++p)
487 {
488 impl::_compute_collisions_point(
489 tree, std::span<const T, 3>(points.data() + 3 * p, 3), entities);
490 offsets[p + 1] = entities.size();
491 }
492
493 return graph::AdjacencyList(std::move(entities), std::move(offsets));
494 }
495 else
496 {
498 std::vector<std::int32_t>(),
499 std::vector<std::int32_t>(points.size() / 3 + 1, 0));
500 }
501}
502
520template <std::floating_point T>
522 std::span<const std::int32_t> cells,
523 std::array<T, 3> point, T tol)
524{
525 if (cells.empty())
526 return -1;
527 else
528 {
529 const mesh::Geometry<T>& geometry = mesh.geometry();
530 std::span<const T> geom_dofs = geometry.x();
531 auto x_dofmap = geometry.dofmap();
532 const std::size_t num_nodes = x_dofmap.extent(1);
533 std::vector<T> coordinate_dofs(num_nodes * 3);
534 for (auto cell : cells)
535 {
536 auto dofs = md::submdspan(x_dofmap, cell, md::full_extent);
537 for (std::size_t i = 0; i < num_nodes; ++i)
538 {
539 std::copy_n(std::next(geom_dofs.begin(), 3 * dofs[i]), 3,
540 std::next(coordinate_dofs.begin(), 3 * i));
541 }
542
543 std::array<T, 3> shortest_vector
544 = compute_distance_gjk<T>(point, coordinate_dofs);
545 T d2 = std::reduce(shortest_vector.begin(), shortest_vector.end(), T(0),
546 [](auto d, auto e) { return d + e * e; });
547 if (d2 < tol)
548 return cell;
549 }
550
551 return -1;
552 }
553}
554
567template <std::floating_point T>
568std::vector<std::int32_t>
570 const BoundingBoxTree<T>& midpoint_tree,
571 const mesh::Mesh<T>& mesh, std::span<const T> points)
572{
573 if (tree.num_bboxes() == 0)
574 return std::vector<std::int32_t>(points.size() / 3, -1);
575
576 std::vector<std::int32_t> entities;
577 entities.reserve(points.size() / 3);
578 for (std::size_t i = 0; i < points.size() / 3; ++i)
579 {
580 // Use midpoint tree to find initial closest entity to the point.
581 // Start by using a leaf node as the initial guess for the input
582 // entity
583 std::array<int, 2> leaf0 = midpoint_tree.bbox(0);
584 assert(impl::is_leaf(leaf0));
585 std::array<T, 6> diff = midpoint_tree.get_bbox(0);
586 for (std::size_t k = 0; k < 3; ++k)
587 diff[k] -= points[3 * i + k];
588 T R2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2];
589
590 // Use a recursive search through the bounding box tree
591 // to find determine the entity with the closest midpoint.
592 // As the midpoint tree only consist of points, the distance
593 // queries are lightweight.
594 const auto [m_index, m_distance2] = impl::_compute_closest_entity(
595 midpoint_tree, std::span<const T, 3>(points.data() + 3 * i, 3),
596 midpoint_tree.num_bboxes() - 1, mesh, leaf0[0], R2);
597
598 // Use a recursives search through the bounding box tree to
599 // determine which entity is actually closest.
600 // Uses the entity with the closest midpoint as initial guess, and
601 // the distance from the midpoint to the point of interest as the
602 // initial search radius.
603 const auto [index, distance2] = impl::_compute_closest_entity(
604 tree, std::span<const T, 3>(points.data() + 3 * i, 3),
605 tree.num_bboxes() - 1, mesh, m_index, m_distance2);
606
607 entities.push_back(index);
608 }
609
610 return entities;
611}
612
628template <std::floating_point T>
630 const mesh::Mesh<T>& mesh,
631 const graph::AdjacencyList<std::int32_t>& candidate_cells,
632 std::span<const T> points)
633{
634 std::vector<std::int32_t> offsets = {0};
635 offsets.reserve(candidate_cells.num_nodes() + 1);
636 std::vector<std::int32_t> colliding_cells;
637 constexpr T eps2 = 1e-12;
638 const int tdim = mesh.topology()->dim();
639 for (std::int32_t i = 0; i < candidate_cells.num_nodes(); i++)
640 {
641 auto cells = candidate_cells.links(i);
642 std::vector<T> _point(3 * cells.size());
643 for (std::size_t j = 0; j < cells.size(); ++j)
644 for (std::size_t k = 0; k < 3; ++k)
645 _point[3 * j + k] = points[3 * i + k];
646
647 std::vector distances_sq = squared_distance<T>(mesh, tdim, cells, _point);
648 for (std::size_t j = 0; j < cells.size(); j++)
649 if (distances_sq[j] < eps2)
650 colliding_cells.push_back(cells[j]);
651
652 offsets.push_back(colliding_cells.size());
653 }
654
655 return graph::AdjacencyList(std::move(colliding_cells), std::move(offsets));
656}
657
685template <std::floating_point T>
687 std::span<const T> points,
688 T padding)
689{
690 MPI_Comm comm = mesh.comm();
691
692 // Create a global bounding-box tree to find candidate processes with
693 // cells that could collide with the points
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();
697 // NOTE: Should we send the cells in as input?
698 std::vector<std::int32_t> cells(num_cells, 0);
699 std::iota(cells.begin(), cells.end(), 0);
700 BoundingBoxTree bb(mesh, tdim, cells, padding);
701 BoundingBoxTree global_bbtree = bb.create_global_tree(comm);
702
703 // Compute collisions:
704 // For each point in `points` get the processes it should be sent to
705 graph::AdjacencyList collisions = compute_collisions(global_bbtree, points);
706
707 // Get unique list of outgoing ranks
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);
712
713 // Compute incoming edges (source processes)
714 std::vector in_ranks = dolfinx::MPI::compute_graph_edges_nbx(comm, out_ranks);
715 std::ranges::sort(in_ranks);
716
717 // Create neighborhood communicator in forward direction
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);
722
723 // Compute map from global mpi rank to neighbor rank, "collisions"
724 // uses global rank
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;
728
729 // Count the number of points to send per neighbor process
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;
734
735 // Compute receive sizes
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);
742
743 // Compute sending offsets
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));
747
748 // Pack data to send and store unpack map
749 std::vector<T> send_data(send_offsets.back());
750 std::vector<std::int32_t> counter(send_sizes.size(), 0);
751 // unpack map: [index in adj list][pos in x]
752 std::vector<std::int32_t> unpack_map(send_offsets.back() / 3);
753 for (std::size_t i = 0; i < points.size(); i += 3)
754 {
755 for (auto p : collisions.links(i / 3))
756 {
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;
763 }
764 }
765
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));
770
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(),
774 dolfinx::MPI::mpi_t<T>, received_points.data(), recv_sizes.data(),
775 recv_offsets.data(), dolfinx::MPI::mpi_t<T>, forward_comm);
776
777 // Get mesh geometry for closest entity
778 const mesh::Geometry<T>& geometry = mesh.geometry();
779 std::span<const T> geom_dofs = geometry.x();
780 auto x_dofmap = geometry.dofmap();
781
782 // Compute candidate cells for collisions (and extrapolation)
783 const graph::AdjacencyList<std::int32_t> candidate_collisions
784 = compute_collisions(bb, std::span<const T>(received_points.data(),
785 received_points.size()));
786
787 // Each process checks which points collide with a cell on the process
788 const int rank = dolfinx::MPI::rank(comm);
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)
792 {
793 std::array<T, 3> point;
794 std::copy_n(std::next(received_points.begin(), p), 3, point.begin());
795 // Find first colliding cell among the cells with colliding bounding boxes
796 const int colliding_cell = geometry::compute_first_colliding_cell(
797 mesh, candidate_collisions.links(p / 3), point,
798 10 * std::numeric_limits<T>::epsilon());
799 // If a collding cell is found, store the rank of the current process
800 // which will be sent back to the owner of the point
801 cell_indicator[p / 3] = (colliding_cell >= 0) ? rank : -1;
802 // Store the cell index for lookup once the owning processes has determined
803 // the ownership of the point
804 closest_cells[p / 3] = colliding_cell;
805 }
806
807 // Create neighborhood communicator in the reverse direction: send
808 // back col to requesting processes
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);
813
814 // Reuse sizes and offsets from first communication set
815 // but divide by three
816 {
817 auto rescale = [](auto& x)
818 { std::ranges::transform(x, x.begin(), [](auto e) { return (e / 3); }); };
819 rescale(recv_sizes);
820 rescale(recv_offsets);
821 rescale(send_sizes);
822 rescale(send_offsets);
823
824 // The communication is reversed, so swap recv to send offsets
825 std::swap(recv_sizes, send_sizes);
826 std::swap(recv_offsets, send_offsets);
827 }
828
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,
833 reverse_comm);
834
835 std::vector<int> point_owners(points.size() / 3, -1);
836 for (std::size_t i = 0; i < unpack_map.size(); i++)
837 {
838 const std::int32_t pos = unpack_map[i];
839 // Only insert new owner if no owner has previously been found
840 if (recv_ranks[i] >= 0 && point_owners[pos] == -1)
841 point_owners[pos] = recv_ranks[i];
842 }
843
844 // Create extrapolation marker for those points already sent to other
845 // process
846 std::vector<std::uint8_t> send_extrapolate(recv_offsets.back());
847 for (std::int32_t i = 0; i < recv_offsets.back(); i++)
848 {
849 const std::int32_t pos = unpack_map[i];
850 send_extrapolate[i] = point_owners[pos] == -1;
851 }
852
853 // Swap communication direction, to send extrapolation marker to other
854 // processes
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);
862
863 std::vector<T> squared_distances(received_points.size() / 3, -1);
864
865 for (std::size_t i = 0; i < dest_extrapolate.size(); i++)
866 {
867 if (dest_extrapolate[i] == 1)
868 {
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());
872
873 // Find shortest distance among cells with colldiing bounding box
874 T shortest_distance = std::numeric_limits<T>::max();
875 std::int32_t closest_cell = -1;
876 for (auto cell : candidate_collisions.links(i))
877 {
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)
881 {
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];
885 }
886 const std::array<T, 3> d = compute_distance_gjk<T>(
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)
890 {
891 shortest_distance = current_distance;
892 closest_cell = cell;
893 }
894 }
895 closest_cells[i] = closest_cell;
896 squared_distances[i] = shortest_distance;
897 }
898 }
899
900 std::swap(recv_sizes, send_sizes);
901 std::swap(recv_offsets, send_offsets);
902
903 // Get distances from closest entity of points that were on the other process
904 std::vector<T> recv_distances(recv_offsets.back());
905 MPI_Neighbor_alltoallv(
906 squared_distances.data(), send_sizes.data(), send_offsets.data(),
907 dolfinx::MPI::mpi_t<T>, recv_distances.data(), recv_sizes.data(),
908 recv_offsets.data(), dolfinx::MPI::mpi_t<T>, reverse_comm);
909
910 // Update point ownership with extrapolation information
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++)
914 {
915 for (std::int32_t j = recv_offsets[i]; j < recv_offsets[i + 1]; j++)
916 {
917 const std::int32_t pos = unpack_map[j];
918 auto current_dist = recv_distances[j];
919 // Update if closer than previous guess and was found
920 if (auto d = closest_distance[pos];
921 (current_dist > 0) and (current_dist < d))
922 {
923 point_owners[pos] = out_ranks[i];
924 closest_distance[pos] = current_dist;
925 }
926 }
927 }
928
929 // Communication is reversed again to send dest ranks to all processes
930 std::swap(send_sizes, recv_sizes);
931 std::swap(send_offsets, recv_offsets);
932
933 // Pack ownership data
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)
937 {
938 for (auto p : collisions.links(i))
939 {
940 int neighbor = rank_to_neighbor[p];
941 send_owners[send_offsets[neighbor] + counter[neighbor]++]
942 = point_owners[i];
943 }
944 }
945
946 // Send ownership info
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,
951 forward_comm);
952
953 // Unpack dest ranks if point owner is this rank
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++)
959 {
960 for (std::int32_t j = recv_offsets[i]; j < recv_offsets[i + 1]; j++)
961 {
962 if (rank == dest_ranks[j])
963 {
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]);
969 }
970 }
971 }
972
973 MPI_Comm_free(&forward_comm);
974 MPI_Comm_free(&reverse_comm);
975 return PointOwnershipData<T>{.src_owner = std::move(point_owners),
976 .dest_owners = std::move(owned_recv_ranks),
977 .dest_points = std::move(owned_recv_points),
978 .dest_cells = std::move(owned_recv_cells)};
979}
980
981} // namespace dolfinx::geometry
Definition BoundingBoxTree.h:206
BoundingBoxTree create_global_tree(MPI_Comm comm) const
Definition BoundingBoxTree.h:335
std::int32_t num_bboxes() const
Return number of bounding boxes.
Definition BoundingBoxTree.h:371
std::array< T, 6 > get_bbox(std::size_t node) const
Return bounding box coordinates for a given node in the tree,.
Definition BoundingBoxTree.h:323
std::array< std::int32_t, 2 > bbox(std::size_t node) const
Definition BoundingBoxTree.h:405
Definition AdjacencyList.h:27
std::span< T > links(std::size_t node)
Definition AdjacencyList.h:111
const std::vector< T > & array() const
Return contiguous array of links for all nodes (const version)
Definition AdjacencyList.h:128
std::int32_t num_nodes() const
Definition AdjacencyList.h:96
Geometry stores the geometry imposed on a mesh.
Definition Geometry.h:34
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
MPI_Datatype mpi_t
Retrieves the MPI data type associated to the provided type.
Definition MPI.h:281
std::vector< int > compute_graph_edges_nbx(MPI_Comm comm, std::span< const int > edges, int tag=static_cast< int >(tag::consensus_nbx))
Determine incoming graph edges using the NBX consensus algorithm.
Definition MPI.cpp:162
int rank(MPI_Comm comm)
Return process rank for the communicator.
Definition MPI.cpp:64
Geometry data structures and algorithms.
Definition BoundingBoxTree.h:22
PointOwnershipData< T > determine_point_ownership(const mesh::Mesh< T > &mesh, std::span< const T > points, T padding)
Given a set of points, determine which process is colliding, using the GJK algorithm on cells to dete...
Definition utils.h:686
std::vector< T > squared_distance(const mesh::Mesh< T > &mesh, int dim, std::span< const std::int32_t > entities, std::span< const T > points)
Compute the squared distance between a point and a mesh entity.
Definition utils.h:172
std::vector< std::int32_t > compute_collisions(const BoundingBoxTree< T > &tree0, const BoundingBoxTree< T > &tree1)
Compute all collisions between two bounding box trees.
Definition utils.h:454
std::int32_t compute_first_colliding_cell(const mesh::Mesh< T > &mesh, std::span< const std::int32_t > cells, std::array< T, 3 > point, T tol)
Given a set of cells, find the first one that collides with a point.
Definition utils.h:521
T compute_squared_distance_bbox(std::span< const T, 6 > b, std::span< const T, 3 > x)
Compute squared distance between point and bounding box.
Definition utils.h:135
graph::AdjacencyList< std::int32_t > compute_colliding_cells(const mesh::Mesh< T > &mesh, const graph::AdjacencyList< std::int32_t > &candidate_cells, std::span< const T > points)
Compute which cells collide with a point.
Definition utils.h:629
std::vector< T > shortest_vector(const mesh::Mesh< T > &mesh, int dim, std::span< const std::int32_t > entities, std::span< const T > points)
Compute the shortest vector from a mesh entity to a point.
Definition utils.h:52
std::array< T, 3 > compute_distance_gjk(std::span< const T > p, std::span< const T > q)
Compute the distance between two convex bodies p and q, each defined by a set of points.
Definition gjk.h:256
BoundingBoxTree< T > create_midpoint_tree(const mesh::Mesh< T > &mesh, int tdim, std::span< const std::int32_t > entities)
Create a bounding box tree for the midpoints of a subset of entities.
Definition utils.h:427
std::vector< std::int32_t > compute_closest_entity(const BoundingBoxTree< T > &tree, const BoundingBoxTree< T > &midpoint_tree, const mesh::Mesh< T > &mesh, std::span< const T > points)
Compute closest mesh entity to a point.
Definition utils.h:569
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32
std::vector< T > compute_midpoints(const Mesh< T > &mesh, int dim, std::span< const std::int32_t > entities)
Compute the midpoints for mesh entities of a given dimension.
Definition utils.h:408
Information on the ownership of points distributed across processes.
Definition utils.h:30
std::vector< T > dest_points
Points that are owned by current process.
Definition utils.h:35
std::vector< std::int32_t > dest_cells
Definition utils.h:37
std::vector< int > dest_owners
Ranks that sent dest_points to current process.
Definition utils.h:34
std::vector< int > src_owner
Definition utils.h:31