Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/mesh_communication.h"
22 :
23 : // libMesh includes
24 : #include "libmesh/libmesh_config.h"
25 : #include "libmesh/libmesh_common.h"
26 : #include "libmesh/libmesh_logging.h"
27 : #include "libmesh/mesh_base.h"
28 : #include "libmesh/mesh_tools.h"
29 : #include "libmesh/parallel_hilbert.h"
30 : #include "libmesh/parallel_sort.h"
31 : #include "libmesh/elem.h"
32 : #include "libmesh/elem_range.h"
33 : #include "libmesh/node_range.h"
34 :
35 : // TIMPI includes
36 : #include "timpi/parallel_implementation.h"
37 : #include "timpi/parallel_sync.h"
38 :
39 : // C/C++ includes
40 : #ifdef LIBMESH_HAVE_LIBHILBERT
41 : # include "hilbert.h"
42 : #endif
43 :
44 :
45 : #if defined(LIBMESH_HAVE_LIBHILBERT) && defined(LIBMESH_HAVE_MPI)
46 :
47 : namespace { // anonymous namespace for helper functions
48 :
49 : using namespace libMesh;
50 :
51 :
52 : // Trying to divide per mesh, not per point
53 26100 : Point invert_bbox (const BoundingBox & bbox)
54 : {
55 13137 : Point p;
56 :
57 52374 : p(0) = (bbox.first(0) == bbox.second(0)) ? 0. :
58 39237 : 1/(bbox.second(0)-bbox.first(0));
59 :
60 : #if LIBMESH_DIM > 1
61 51530 : p(1) = (bbox.first(1) == bbox.second(1)) ? 0. :
62 36705 : 1/(bbox.second(1)-bbox.first(1));
63 : #endif
64 :
65 : #if LIBMESH_DIM > 2
66 52245 : p(2) = (bbox.first(2) == bbox.second(2)) ? 0. :
67 19450 : 1/(bbox.second(2)-bbox.first(2));
68 : #endif
69 :
70 26100 : return p;
71 : }
72 :
73 :
74 :
75 : // Utility function to map (x,y,z) in [bbox.min, bbox.max]^3 into
76 : // [0,max_inttype]^3 for computing Hilbert keys
77 18976192 : void get_hilbert_coords (const Point & p,
78 : const BoundingBox & bbox,
79 : const Point & bboxinv,
80 : CFixBitVec icoords[3])
81 : {
82 : static const Hilbert::inttype max_inttype = static_cast<Hilbert::inttype>(-1);
83 :
84 : const Real // put (x,y,z) in [0,1]^3 (don't divide by 0)
85 18976192 : x = (p(0)-bbox.first(0)) * bboxinv(0),
86 :
87 : #if LIBMESH_DIM > 1
88 18976192 : y = (p(1)-bbox.first(1)) * bboxinv(1),
89 : #else
90 : y = 0.,
91 : #endif
92 :
93 : #if LIBMESH_DIM > 2
94 18976192 : z = (p(2)-bbox.first(2)) * bboxinv(2);
95 : #else
96 : z = 0.;
97 : #endif
98 :
99 : // (icoords) in [0,max_inttype]^3
100 18976192 : icoords[0] = static_cast<Hilbert::inttype>(x*max_inttype);
101 18976192 : icoords[1] = static_cast<Hilbert::inttype>(y*max_inttype);
102 18976192 : icoords[2] = static_cast<Hilbert::inttype>(z*max_inttype);
103 18976192 : }
104 :
105 :
106 :
107 : Parallel::DofObjectKey
108 12007594 : get_dofobject_key (const Elem & e,
109 : const BoundingBox & bbox,
110 : const Point & bboxinv)
111 : {
112 : static const unsigned int sizeof_inttype = sizeof(Hilbert::inttype);
113 :
114 4026716 : Hilbert::HilbertIndices index;
115 48030376 : CFixBitVec icoords[3];
116 4026716 : Hilbert::BitVecType bv;
117 12007594 : get_hilbert_coords (e.vertex_average(), bbox, bboxinv, icoords);
118 4026716 : Hilbert::coordsToIndex (icoords, 8*sizeof_inttype, 3, bv);
119 12007594 : index = bv;
120 :
121 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
122 16034310 : return std::make_pair(index, e.unique_id());
123 : #else
124 : // A parent and a child can share a vertex average; their levels differ
125 : return std::make_pair(index, unique_id_type(e.level()));
126 : #endif
127 : }
128 :
129 :
130 :
131 : // Compute the hilbert index
132 : Parallel::DofObjectKey
133 6968598 : get_dofobject_key (const Node & n,
134 : const BoundingBox & bbox,
135 : const Point & bboxinv)
136 : {
137 : static const unsigned int sizeof_inttype = sizeof(Hilbert::inttype);
138 :
139 2303645 : Hilbert::HilbertIndices index;
140 27874392 : CFixBitVec icoords[3];
141 2303645 : Hilbert::BitVecType bv;
142 6968598 : get_hilbert_coords (n, bbox, bboxinv, icoords);
143 2303645 : Hilbert::coordsToIndex (icoords, 8*sizeof_inttype, 3, bv);
144 6968598 : index = bv;
145 :
146 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
147 9272243 : return std::make_pair(index, n.unique_id());
148 : #else
149 : return std::make_pair(index, unique_id_type(0));
150 : #endif
151 : }
152 :
153 : // Helper class for threaded Hilbert key computation
154 : class ComputeHilbertKeys
155 : {
156 : public:
157 1140 : ComputeHilbertKeys (const BoundingBox & bbox,
158 : const Point & bboxinv,
159 3444 : std::vector<Parallel::DofObjectKey> & keys) :
160 1164 : _bbox(bbox),
161 1164 : _bboxinv(bboxinv),
162 3444 : _keys(keys)
163 1140 : {}
164 :
165 : // computes the hilbert index for a node
166 2334 : void operator() (const ConstNodeRange & range) const
167 : {
168 1548 : std::size_t pos = range.first_idx();
169 1405468 : for (const auto & node : range)
170 : {
171 460729 : libmesh_assert(node);
172 460729 : libmesh_assert_less (pos, _keys.size());
173 1403134 : _keys[pos++] = get_dofobject_key (*node, _bbox, _bboxinv);
174 : }
175 2334 : }
176 :
177 : // computes the hilbert index for an element
178 2334 : void operator() (const ConstElemRange & range) const
179 : {
180 1548 : std::size_t pos = range.first_idx();
181 1725708 : for (const auto & elem : range)
182 : {
183 571452 : libmesh_assert(elem);
184 571452 : libmesh_assert_less (pos, _keys.size());
185 1723374 : _keys[pos++] = get_dofobject_key (*elem, _bbox, _bboxinv);
186 : }
187 2334 : }
188 :
189 : private:
190 : const BoundingBox _bbox;
191 : const Point _bboxinv;
192 : std::vector<Parallel::DofObjectKey> & _keys;
193 : };
194 :
195 : }
196 : #endif // defined(LIBMESH_HAVE_LIBHILBERT) && defined(LIBMESH_HAVE_MPI)
197 :
198 :
199 : namespace libMesh
200 : {
201 :
202 : // ------------------------------------------------------------
203 : // MeshCommunication class members
204 : #if defined(LIBMESH_HAVE_LIBHILBERT) && defined(LIBMESH_HAVE_MPI)
205 1722 : void MeshCommunication::assign_global_indices (MeshBase & mesh) const
206 : {
207 1710 : LOG_SCOPE ("assign_global_indices()", "MeshCommunication");
208 :
209 : // This method determines partition-agnostic global indices
210 : // for nodes and elements.
211 :
212 : // Algorithm:
213 : // (1) compute the Hilbert key for each local node/element
214 : // (2) perform a parallel sort of the Hilbert key
215 : // (3) get the min/max value on each processor
216 : // (4) determine the position in the global ranking for
217 : // each local object
218 :
219 1140 : const Parallel::Communicator & communicator (mesh.comm());
220 :
221 : #ifdef DEBUG
222 : // This is going to be a mess if geometry is out of sync
223 570 : MeshTools::libmesh_assert_equal_points(mesh);
224 570 : MeshTools::libmesh_assert_equal_connectivity(mesh);
225 : #endif
226 :
227 : // Global bounding box. We choose the nodal bounding box for
228 : // backwards compatibility; the element bounding box may be looser
229 : // on curved elements.
230 : const BoundingBox bbox =
231 1722 : MeshTools::create_nodal_bounding_box (mesh);
232 :
233 1152 : const Point bboxinv = invert_bbox(bbox);
234 :
235 : //-------------------------------------------------------------
236 : // (1) compute Hilbert keys
237 : std::vector<Parallel::DofObjectKey>
238 1140 : node_keys, elem_keys;
239 :
240 : {
241 : // Nodes first
242 : {
243 2862 : ConstNodeRange nr (mesh.local_nodes_begin(),
244 6888 : mesh.local_nodes_end());
245 1722 : node_keys.resize (nr.size());
246 1722 : Threads::parallel_for (nr, ComputeHilbertKeys (bbox, bboxinv, node_keys));
247 :
248 : // // It's O(N^2) to check that these keys don't duplicate before the
249 : // // sort...
250 : //
251 : // MeshBase::const_node_iterator nodei = mesh.local_nodes_begin();
252 : // for (std::size_t i = 0; i != node_keys.size(); ++i, ++nodei)
253 : // {
254 : // MeshBase::const_node_iterator nodej = mesh.local_nodes_begin();
255 : // for (std::size_t j = 0; j != i; ++j, ++nodej)
256 : // {
257 : // if (node_keys[i] == node_keys[j])
258 : // {
259 : // CFixBitVec icoords[3], jcoords[3];
260 : // get_hilbert_coords(**nodej, bbox, bboxinv, jcoords);
261 : // libMesh::err << "node " << (*nodej)->id() << ", " << static_cast<Point &>(**nodej) << " has HilbertIndices " << node_keys[j] << std::endl;
262 : // get_hilbert_coords(**nodei, bbox, bboxinv, icoords);
263 : // libMesh::err << "node " << (*nodei)->id() << ", " << static_cast<Point &>(**nodei) << " has HilbertIndices " << node_keys[i] << std::endl;
264 : // libmesh_error_msg("Error: nodes with duplicate Hilbert keys!");
265 : // }
266 : // }
267 : // }
268 : }
269 :
270 : // Elements next
271 : {
272 2862 : ConstElemRange er (mesh.local_elements_begin(),
273 6888 : mesh.local_elements_end());
274 1722 : elem_keys.resize (er.size());
275 1722 : Threads::parallel_for (er, ComputeHilbertKeys (bbox, bboxinv, elem_keys));
276 :
277 : // // For elements, the keys can be (and in the case of TRI, are
278 : // // expected to be) duplicates, but only if the elements are at
279 : // // different levels
280 : // MeshBase::const_element_iterator elemi = mesh.local_elements_begin();
281 : // for (std::size_t i = 0; i != elem_keys.size(); ++i, ++elemi)
282 : // {
283 : // MeshBase::const_element_iterator elemj = mesh.local_elements_begin();
284 : // for (std::size_t j = 0; j != i; ++j, ++elemj)
285 : // {
286 : // if ((elem_keys[i] == elem_keys[j]) &&
287 : // ((*elemi)->level() == (*elemj)->level()))
288 : // {
289 : // libMesh::err << "level " << (*elemj)->level()
290 : // << " elem\n" << (**elemj)
291 : // << " vertex average " << (*elemj)->vertex_average()
292 : // << " has HilbertIndices " << elem_keys[j]
293 : // << " or " << get_dofobject_key((**elemj), bbox, bboxinv)
294 : // << std::endl;
295 : // libMesh::err << "level " << (*elemi)->level()
296 : // << " elem\n" << (**elemi)
297 : // << " vertex average " << (*elemi)->vertex_average()
298 : // << " has HilbertIndices " << elem_keys[i]
299 : // << " or " << get_dofobject_key((**elemi), bbox, bboxinv)
300 : // << std::endl;
301 : // libmesh_error_msg("Error: level " << (*elemi)->level() << " elements with duplicate Hilbert keys!");
302 : // }
303 : // }
304 : // }
305 : }
306 : } // done computing Hilbert keys
307 :
308 :
309 :
310 : //-------------------------------------------------------------
311 : // (2) parallel sort the Hilbert keys
312 : Parallel::Sort<Parallel::DofObjectKey> node_sorter (communicator,
313 2862 : node_keys);
314 1722 : node_sorter.sort(); /* done with node_keys */ //node_keys.clear();
315 :
316 : const std::vector<Parallel::DofObjectKey> & my_node_bin =
317 1722 : node_sorter.bin();
318 :
319 : Parallel::Sort<Parallel::DofObjectKey> elem_sorter (communicator,
320 2862 : elem_keys);
321 1722 : elem_sorter.sort(); /* done with elem_keys */ //elem_keys.clear();
322 :
323 : const std::vector<Parallel::DofObjectKey> & my_elem_bin =
324 1722 : elem_sorter.bin();
325 :
326 :
327 :
328 : //-------------------------------------------------------------
329 : // (3) get the max value on each processor
330 : std::vector<Parallel::DofObjectKey>
331 2292 : node_upper_bounds(communicator.size()),
332 2292 : elem_upper_bounds(communicator.size());
333 :
334 : { // limit scope of temporaries
335 2292 : std::vector<Parallel::DofObjectKey> recvbuf(2*communicator.size());
336 : std::vector<unsigned short int> /* do not use a vector of bools here since it is not always so! */
337 2292 : empty_nodes (communicator.size()),
338 2292 : empty_elem (communicator.size());
339 2292 : std::vector<Parallel::DofObjectKey> my_max(2);
340 :
341 2292 : communicator.allgather (static_cast<unsigned short int>(my_node_bin.empty()), empty_nodes);
342 1722 : communicator.allgather (static_cast<unsigned short int>(my_elem_bin.empty()), empty_elem);
343 :
344 1722 : if (!my_node_bin.empty()) my_max[0] = my_node_bin.back();
345 1722 : if (!my_elem_bin.empty()) my_max[1] = my_elem_bin.back();
346 :
347 1722 : communicator.allgather (my_max, /* identical_buffer_sizes = */ true);
348 :
349 : // Be careful here. The *_upper_bounds will be used to find the processor
350 : // a given object belongs to. So, if a processor contains no objects (possible!)
351 : // then copy the bound from the lower processor id.
352 5148 : for (auto p : make_range(communicator.size()))
353 : {
354 4566 : node_upper_bounds[p] = my_max[2*p+0];
355 4566 : elem_upper_bounds[p] = my_max[2*p+1];
356 :
357 3426 : if (p > 0) // default hilbert index value is the OK upper bound for processor 0.
358 : {
359 2274 : if (empty_nodes[p]) node_upper_bounds[p] = node_upper_bounds[p-1];
360 2274 : if (empty_elem[p]) elem_upper_bounds[p] = elem_upper_bounds[p-1];
361 : }
362 : }
363 : }
364 :
365 :
366 :
367 : //-------------------------------------------------------------
368 : // (4) determine the position in the global ranking for
369 : // each local object
370 : {
371 : //----------------------------------------------
372 : // Nodes first -- all nodes, not just local ones
373 : {
374 : // Request sets to send to each processor
375 : std::map<processor_id_type, std::vector<Parallel::DofObjectKey>>
376 1140 : requested_ids;
377 : // Results to gather from each processor - kept in a map so we
378 : // do only one loop over nodes after all receives are done.
379 : std::map<dof_id_type, std::vector<dof_id_type>>
380 1140 : filled_request;
381 :
382 : // build up list of requests
383 3725422 : for (const auto & node : mesh.node_ptr_range())
384 : {
385 921458 : libmesh_assert(node);
386 : const Parallel::DofObjectKey hi =
387 2782732 : get_dofobject_key (*node, bbox, bboxinv);
388 : const processor_id_type pid =
389 : cast_int<processor_id_type>
390 2782732 : (std::distance (node_upper_bounds.begin(),
391 : std::lower_bound(node_upper_bounds.begin(),
392 : node_upper_bounds.end(),
393 2782732 : hi)));
394 :
395 921458 : libmesh_assert_less (pid, communicator.size());
396 :
397 2782732 : requested_ids[pid].push_back(hi);
398 582 : }
399 :
400 : // The number of objects in my_node_bin on each processor
401 2292 : std::vector<dof_id_type> node_bin_sizes(communicator.size());
402 2292 : communicator.allgather (static_cast<dof_id_type>(my_node_bin.size()), node_bin_sizes);
403 :
404 : // The offset of my first global index
405 570 : dof_id_type my_offset = 0;
406 2574 : for (auto pid : make_range(communicator.rank()))
407 1137 : my_offset += node_bin_sizes[pid];
408 :
409 : auto gather_functor =
410 1146 : [
411 : #ifndef NDEBUG
412 : & node_upper_bounds,
413 : & communicator,
414 : #endif
415 : & my_node_bin,
416 : my_offset
417 : ]
418 : (processor_id_type,
419 : const std::vector<Parallel::DofObjectKey> & keys,
420 9235218 : std::vector<dof_id_type> & global_ids)
421 : {
422 : // Fill the requests
423 2280 : const std::size_t keys_size = keys.size();
424 3426 : global_ids.reserve(keys_size);
425 2786158 : for (std::size_t idx=0; idx != keys_size; idx++)
426 : {
427 1842916 : const Parallel::DofObjectKey & hi = keys[idx];
428 921458 : libmesh_assert_less_equal (hi, node_upper_bounds[communicator.rank()]);
429 :
430 : // find the requested index in my node bin
431 : std::vector<Parallel::DofObjectKey>::const_iterator pos =
432 2782732 : std::lower_bound (my_node_bin.begin(), my_node_bin.end(), hi);
433 921458 : libmesh_assert (pos != my_node_bin.end());
434 921458 : libmesh_assert_equal_to (*pos, hi);
435 :
436 : // Finally, assign the global index based off the position of the index
437 : // in my array, properly offset.
438 4625648 : global_ids.push_back(cast_int<dof_id_type>(std::distance(my_node_bin.begin(), pos) + my_offset));
439 : }
440 4578 : };
441 :
442 : auto action_functor =
443 1146 : [&filled_request]
444 : (processor_id_type pid,
445 : const std::vector<Parallel::DofObjectKey> &,
446 2286 : const std::vector<dof_id_type> & global_ids)
447 : {
448 3426 : filled_request[pid] = global_ids;
449 4008 : };
450 :
451 : // Trade requests with other processors
452 570 : const dof_id_type * ex = nullptr;
453 : Parallel::pull_parallel_vector_data
454 1722 : (communicator, requested_ids, gather_functor, action_functor, ex);
455 :
456 : // We now have all the filled requests, so we can loop through our
457 : // nodes once and assign the global index to each one.
458 : {
459 : std::map<dof_id_type, std::vector<dof_id_type>::const_iterator>
460 1140 : next_obj_on_proc;
461 5148 : for (auto & p : filled_request)
462 3426 : next_obj_on_proc[p.first] = p.second.begin();
463 :
464 2809142 : for (auto & node : mesh.node_ptr_range())
465 : {
466 921458 : libmesh_assert(node);
467 : const Parallel::DofObjectKey hi =
468 2782732 : get_dofobject_key (*node, bbox, bboxinv);
469 : const processor_id_type pid =
470 : cast_int<processor_id_type>
471 2782732 : (std::distance (node_upper_bounds.begin(),
472 : std::lower_bound(node_upper_bounds.begin(),
473 : node_upper_bounds.end(),
474 1866452 : hi)));
475 :
476 921458 : libmesh_assert_less (pid, communicator.size());
477 921458 : libmesh_assert (next_obj_on_proc[pid] != filled_request[pid].end());
478 :
479 2782732 : const dof_id_type global_index = *next_obj_on_proc[pid];
480 921458 : libmesh_assert_less (global_index, mesh.n_nodes());
481 2782732 : node->set_id() = global_index;
482 :
483 2782732 : ++next_obj_on_proc[pid];
484 582 : }
485 : }
486 : }
487 :
488 : //---------------------------------------------------
489 : // elements next -- all elements, not just local ones
490 : {
491 : // Request sets to send to each processor
492 : std::map<processor_id_type, std::vector<Parallel::DofObjectKey>>
493 1140 : requested_ids;
494 : // Results to gather from each processor - kept in a map so we
495 : // do only one loop over elements after all receives are done.
496 : std::map<dof_id_type, std::vector<dof_id_type>>
497 1140 : filled_request;
498 :
499 4586878 : for (const auto & elem : mesh.element_ptr_range())
500 : {
501 1142904 : libmesh_assert(elem);
502 : const Parallel::DofObjectKey hi =
503 3434906 : get_dofobject_key (*elem, bbox, bboxinv);
504 : const processor_id_type pid =
505 : cast_int<processor_id_type>
506 3434906 : (std::distance (elem_upper_bounds.begin(),
507 : std::lower_bound(elem_upper_bounds.begin(),
508 : elem_upper_bounds.end(),
509 3434906 : hi)));
510 :
511 1142904 : libmesh_assert_less (pid, communicator.size());
512 :
513 3434906 : requested_ids[pid].push_back(hi);
514 582 : }
515 :
516 : // The number of objects in my_elem_bin on each processor
517 2292 : std::vector<dof_id_type> elem_bin_sizes(communicator.size());
518 2292 : communicator.allgather (static_cast<dof_id_type>(my_elem_bin.size()), elem_bin_sizes);
519 :
520 : // The offset of my first global index
521 570 : dof_id_type my_offset = 0;
522 2574 : for (auto pid : make_range(communicator.rank()))
523 1137 : my_offset += elem_bin_sizes[pid];
524 :
525 : auto gather_functor =
526 1146 : [
527 : #ifndef NDEBUG
528 : & elem_upper_bounds,
529 : & communicator,
530 : #endif
531 : & my_elem_bin,
532 : my_offset
533 : ]
534 : (processor_id_type,
535 : const std::vector<Parallel::DofObjectKey> & keys,
536 11437514 : std::vector<dof_id_type> & global_ids)
537 : {
538 : // Fill the requests
539 2280 : const std::size_t keys_size = keys.size();
540 3426 : global_ids.reserve(keys_size);
541 3438332 : for (std::size_t idx=0; idx != keys_size; idx++)
542 : {
543 2285808 : const Parallel::DofObjectKey & hi = keys[idx];
544 1142904 : libmesh_assert_less_equal (hi, elem_upper_bounds[communicator.rank()]);
545 :
546 : // find the requested index in my elem bin
547 : std::vector<Parallel::DofObjectKey>::const_iterator pos =
548 3434906 : std::lower_bound (my_elem_bin.begin(), my_elem_bin.end(), hi);
549 1142904 : libmesh_assert (pos != my_elem_bin.end());
550 1142904 : libmesh_assert_equal_to (*pos, hi);
551 :
552 : // Finally, assign the global index based off the position of the index
553 : // in my array, properly offset.
554 5720714 : global_ids.push_back (cast_int<dof_id_type>(std::distance(my_elem_bin.begin(), pos) + my_offset));
555 : }
556 4578 : };
557 :
558 : auto action_functor =
559 1146 : [&filled_request]
560 : (processor_id_type pid,
561 : const std::vector<Parallel::DofObjectKey> &,
562 2286 : const std::vector<dof_id_type> & global_ids)
563 : {
564 3426 : filled_request[pid] = global_ids;
565 4008 : };
566 :
567 : // Trade requests with other processors
568 570 : const dof_id_type * ex = nullptr;
569 : Parallel::pull_parallel_vector_data
570 1722 : (communicator, requested_ids, gather_functor, action_functor, ex);
571 :
572 : // We now have all the filled requests, so we can loop through our
573 : // elements once and assign the global index to each one.
574 : {
575 : std::vector<std::vector<dof_id_type>::const_iterator>
576 2292 : next_obj_on_proc; next_obj_on_proc.reserve(communicator.size());
577 5148 : for (auto pid : make_range(communicator.size()))
578 5712 : next_obj_on_proc.push_back(filled_request[pid].begin());
579 :
580 4586878 : for (auto & elem : mesh.element_ptr_range())
581 : {
582 1142904 : libmesh_assert(elem);
583 : const Parallel::DofObjectKey hi =
584 3434906 : get_dofobject_key (*elem, bbox, bboxinv);
585 : const processor_id_type pid =
586 : cast_int<processor_id_type>
587 3434906 : (std::distance (elem_upper_bounds.begin(),
588 : std::lower_bound(elem_upper_bounds.begin(),
589 : elem_upper_bounds.end(),
590 1142904 : hi)));
591 :
592 1142904 : libmesh_assert_less (pid, communicator.size());
593 1142904 : libmesh_assert (next_obj_on_proc[pid] != filled_request[pid].end());
594 :
595 3434906 : const dof_id_type global_index = *next_obj_on_proc[pid];
596 1142904 : libmesh_assert_less (global_index, mesh.n_elem());
597 3434906 : elem->set_id() = global_index;
598 :
599 1142904 : ++next_obj_on_proc[pid];
600 582 : }
601 : }
602 : }
603 : }
604 2304 : }
605 : #else // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
606 : void MeshCommunication::assign_global_indices (MeshBase &) const
607 : {
608 : }
609 : #endif // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
610 :
611 :
612 : #if defined(LIBMESH_HAVE_LIBHILBERT) && defined(LIBMESH_HAVE_MPI)
613 0 : void MeshCommunication::check_for_duplicate_global_indices (MeshBase & mesh) const
614 : {
615 0 : LOG_SCOPE ("check_for_duplicate_global_indices()", "MeshCommunication");
616 :
617 : // Global bounding box. We choose the nodal bounding box for
618 : // backwards compatibility; the element bounding box may be looser
619 : // on curved elements.
620 : const BoundingBox bbox =
621 0 : MeshTools::create_nodal_bounding_box (mesh);
622 :
623 0 : const Point bboxinv = invert_bbox(bbox);
624 :
625 : std::vector<Parallel::DofObjectKey>
626 0 : node_keys, elem_keys;
627 :
628 : {
629 : // Nodes first
630 : {
631 0 : ConstNodeRange nr (mesh.local_nodes_begin(),
632 0 : mesh.local_nodes_end());
633 0 : node_keys.resize (nr.size());
634 0 : Threads::parallel_for (nr, ComputeHilbertKeys (bbox, bboxinv, node_keys));
635 :
636 : // It's O(N^2) to check that these keys don't duplicate before the
637 : // sort...
638 0 : MeshBase::const_node_iterator nodei = mesh.local_nodes_begin();
639 0 : for (std::size_t i = 0; i != node_keys.size(); ++i, ++nodei)
640 : {
641 0 : MeshBase::const_node_iterator nodej = mesh.local_nodes_begin();
642 0 : for (std::size_t j = 0; j != i; ++j, ++nodej)
643 : {
644 0 : if (node_keys[i] == node_keys[j])
645 : {
646 0 : CFixBitVec icoords[3], jcoords[3];
647 0 : get_hilbert_coords(**nodej, bbox, bboxinv, jcoords);
648 : libMesh::err <<
649 0 : "node " << (*nodej)->id() << ", " <<
650 0 : *(const Point *)(*nodej) << " has HilbertIndices " <<
651 0 : node_keys[j] << std::endl;
652 0 : get_hilbert_coords(**nodei, bbox, bboxinv, icoords);
653 : libMesh::err <<
654 0 : "node " << (*nodei)->id() << ", " <<
655 0 : *(const Point *)(*nodei) << " has HilbertIndices " <<
656 0 : node_keys[i] << std::endl;
657 0 : libmesh_error_msg("Error: nodes with duplicate Hilbert keys!");
658 : }
659 : }
660 : }
661 : }
662 :
663 : // Elements next
664 : {
665 0 : ConstElemRange er (mesh.local_elements_begin(),
666 0 : mesh.local_elements_end());
667 0 : elem_keys.resize (er.size());
668 0 : Threads::parallel_for (er, ComputeHilbertKeys (bbox, bboxinv, elem_keys));
669 :
670 : // For elements, the keys can be (and in the case of TRI, are
671 : // expected to be) duplicates, but only if the elements are at
672 : // different levels
673 0 : MeshBase::const_element_iterator elemi = mesh.local_elements_begin();
674 0 : for (std::size_t i = 0; i != elem_keys.size(); ++i, ++elemi)
675 : {
676 0 : MeshBase::const_element_iterator elemj = mesh.local_elements_begin();
677 0 : for (std::size_t j = 0; j != i; ++j, ++elemj)
678 : {
679 0 : if ((elem_keys[i] == elem_keys[j]) &&
680 0 : ((*elemi)->level() == (*elemj)->level()))
681 : {
682 : libMesh::err <<
683 0 : "level " << (*elemj)->level() << " elem\n" <<
684 0 : (**elemj) << " vertex average " <<
685 0 : (*elemj)->vertex_average() << " has HilbertIndices " <<
686 0 : elem_keys[j] << " or " <<
687 0 : get_dofobject_key((**elemj), bbox, bboxinv) <<
688 0 : std::endl;
689 : libMesh::err <<
690 0 : "level " << (*elemi)->level() << " elem\n" <<
691 0 : (**elemi) << " vertex average " <<
692 0 : (*elemi)->vertex_average() << " has HilbertIndices " <<
693 0 : elem_keys[i] << " or " <<
694 0 : get_dofobject_key((**elemi), bbox, bboxinv) <<
695 0 : std::endl;
696 0 : libmesh_error_msg("Error: level " << (*elemi)->level() << " elements with duplicate Hilbert keys!");
697 : }
698 : }
699 : }
700 : }
701 : } // done checking Hilbert keys
702 0 : }
703 : #else // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
704 : void MeshCommunication::check_for_duplicate_global_indices (MeshBase &) const
705 : {
706 : }
707 : #endif // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
708 :
709 : #if defined(LIBMESH_HAVE_LIBHILBERT) && defined(LIBMESH_HAVE_MPI)
710 : template <typename ForwardIterator>
711 2340 : void MeshCommunication::find_local_indices (const BoundingBox & bbox,
712 : const ForwardIterator & begin,
713 : const ForwardIterator & end,
714 : std::unordered_map<dof_id_type, dof_id_type> & index_map) const
715 : {
716 2340 : LOG_SCOPE ("find_local_indices()", "MeshCommunication");
717 :
718 : // This method determines id-agnostic local indices
719 : // for nodes and elements by sorting Hilbert keys.
720 :
721 780 : index_map.clear();
722 :
723 1560 : const Point bboxinv = invert_bbox(bbox);
724 :
725 : //-------------------------------------------------------------
726 : // (1) compute Hilbert keys
727 : // These aren't trivial to compute, and we will need them again.
728 : // But the binsort will sort the input vector, trashing the order
729 : // that we'd like to rely on. So, two vectors...
730 1560 : std::map<Parallel::DofObjectKey, dof_id_type> hilbert_keys;
731 : {
732 2340 : LOG_SCOPE("local_hilbert_indices", "MeshCommunication");
733 73344 : for (ForwardIterator it=begin; it!=end; ++it)
734 : {
735 71004 : const Parallel::DofObjectKey hi(get_dofobject_key ((**it), bbox, bboxinv));
736 71004 : hilbert_keys.emplace(hi, (*it)->id());
737 : }
738 : }
739 :
740 : {
741 780 : dof_id_type cnt = 0;
742 73344 : for (auto key_val : hilbert_keys)
743 71004 : index_map[key_val.second] = cnt++;
744 : }
745 2340 : }
746 :
747 :
748 : template <typename ForwardIterator>
749 35175 : void MeshCommunication::find_global_indices (const Parallel::Communicator & communicator,
750 : const BoundingBox & bbox,
751 : const ForwardIterator & begin,
752 : const ForwardIterator & end,
753 : std::vector<dof_id_type> & index_map) const
754 : {
755 35361 : LOG_SCOPE ("find_global_indices()", "MeshCommunication");
756 :
757 : // This method determines partition-agnostic global indices
758 : // for nodes and elements.
759 :
760 : // Algorithm:
761 : // (1) compute the Hilbert key for each local node/element
762 : // (2) perform a parallel sort of the Hilbert key
763 : // (3) get the min/max value on each processor
764 : // (4) determine the position in the global ranking for
765 : // each local object
766 11787 : index_map.clear();
767 35175 : std::size_t n_objects = std::distance (begin, end);
768 35175 : index_map.reserve(n_objects);
769 :
770 23388 : const Point bboxinv = invert_bbox(bbox);
771 :
772 : //-------------------------------------------------------------
773 : // (1) compute Hilbert keys
774 : // These aren't trivial to compute, and we will need them again.
775 : // But the binsort will sort the input vector, trashing the order
776 : // that we'd like to rely on. So, two vectors...
777 : std::vector<Parallel::DofObjectKey>
778 23574 : sorted_hilbert_keys,
779 23574 : hilbert_keys;
780 35175 : sorted_hilbert_keys.reserve(n_objects);
781 35175 : hilbert_keys.reserve(n_objects);
782 : {
783 35361 : LOG_SCOPE("compute_hilbert_indices()", "MeshCommunication");
784 23574 : const processor_id_type my_pid = communicator.rank();
785 4430407 : for (ForwardIterator it=begin; it!=end; ++it)
786 : {
787 3343404 : const Parallel::DofObjectKey hi(get_dofobject_key (**it, bbox, bboxinv));
788 3343404 : hilbert_keys.push_back(hi);
789 :
790 3343404 : const processor_id_type pid = (*it)->processor_id();
791 :
792 3343404 : if (pid == my_pid)
793 1667816 : sorted_hilbert_keys.push_back(hi);
794 :
795 : // someone needs to take care of unpartitioned objects!
796 4489192 : if ((my_pid == 0) &&
797 2197616 : (pid == DofObject::invalid_processor_id))
798 9762 : sorted_hilbert_keys.push_back(hi);
799 : }
800 : }
801 :
802 : //-------------------------------------------------------------
803 : // (2) parallel sort the Hilbert keys
804 23574 : START_LOG ("parallel_sort()", "MeshCommunication");
805 58749 : Parallel::Sort<Parallel::DofObjectKey> sorter (communicator,
806 : sorted_hilbert_keys);
807 35175 : sorter.sort();
808 23574 : STOP_LOG ("parallel_sort()", "MeshCommunication");
809 35175 : const std::vector<Parallel::DofObjectKey> & my_bin = sorter.bin();
810 :
811 : // The number of objects in my_bin on each processor
812 46962 : std::vector<unsigned int> bin_sizes(communicator.size());
813 46962 : communicator.allgather (static_cast<unsigned int>(my_bin.size()), bin_sizes);
814 :
815 : // The offset of my first global index
816 11787 : unsigned int my_offset = 0;
817 52758 : for (auto pid : make_range(communicator.rank()))
818 23476 : my_offset += bin_sizes[pid];
819 :
820 : //-------------------------------------------------------------
821 : // (3) get the max value on each processor
822 : std::vector<Parallel::DofObjectKey>
823 46962 : upper_bounds(1);
824 :
825 35175 : if (!my_bin.empty())
826 23574 : upper_bounds[0] = my_bin.back();
827 :
828 35175 : communicator.allgather (upper_bounds, /* identical_buffer_sizes = */ true);
829 :
830 : // Be careful here. The *_upper_bounds will be used to find the processor
831 : // a given object belongs to. So, if a processor contains no objects (possible!)
832 : // then copy the bound from the lower processor id.
833 70341 : for (auto p : IntRange<processor_id_type>(1, communicator.size()))
834 46952 : if (!bin_sizes[p]) upper_bounds[p] = upper_bounds[p-1];
835 :
836 :
837 : //-------------------------------------------------------------
838 : // (4) determine the position in the global ranking for
839 : // each local object
840 : {
841 : //----------------------------------------------
842 : // all objects, not just local ones
843 :
844 : // Request sets to send to each processor
845 : std::map<processor_id_type, std::vector<Parallel::DofObjectKey>>
846 23574 : requested_ids;
847 : // Results to gather from each processor
848 : std::map<processor_id_type, std::vector<dof_id_type>>
849 23574 : filled_request;
850 :
851 : // build up list of requests
852 11787 : std::vector<Parallel::DofObjectKey>::const_iterator hi =
853 11787 : hilbert_keys.begin();
854 :
855 3378696 : for (ForwardIterator it = begin; it != end; ++it)
856 : {
857 1145788 : libmesh_assert (hi != hilbert_keys.end());
858 :
859 : std::vector<Parallel::DofObjectKey>::iterator lb =
860 3343404 : std::lower_bound(upper_bounds.begin(), upper_bounds.end(),
861 1145788 : *hi);
862 :
863 3343404 : const processor_id_type pid =
864 : cast_int<processor_id_type>
865 1145788 : (std::distance (upper_bounds.begin(), lb));
866 :
867 1145788 : libmesh_assert_less (pid, communicator.size());
868 :
869 3343404 : requested_ids[pid].push_back(*hi);
870 :
871 1145788 : ++hi;
872 : // go ahead and put pid in index_map, that way we
873 : // don't have to repeat the std::lower_bound()
874 3343404 : index_map.push_back(pid);
875 : }
876 :
877 3519683 : auto gather_functor =
878 1075022 : [
879 : #ifndef NDEBUG
880 : & upper_bounds,
881 : & communicator,
882 : #endif
883 : & bbox,
884 : & my_bin,
885 : my_offset
886 : ]
887 : (processor_id_type, const std::vector<Parallel::DofObjectKey> & keys,
888 47144 : std::vector<dof_id_type> & global_ids)
889 : {
890 : // Ignore unused lambda capture warnings in devel mode
891 23572 : libmesh_ignore(bbox);
892 :
893 : // Fill the requests
894 47144 : const std::size_t keys_size = keys.size();
895 23572 : global_ids.clear();
896 70338 : global_ids.reserve(keys_size);
897 3413742 : for (std::size_t idx=0; idx != keys_size; idx++)
898 : {
899 2291576 : const Parallel::DofObjectKey & hilbert_indices = keys[idx];
900 1145788 : libmesh_assert_less_equal (hilbert_indices, upper_bounds[communicator.rank()]);
901 :
902 : // find the requested index in my node bin
903 : std::vector<Parallel::DofObjectKey>::const_iterator pos =
904 3343404 : std::lower_bound (my_bin.begin(), my_bin.end(), hilbert_indices);
905 1145788 : libmesh_assert (pos != my_bin.end());
906 : #ifdef DEBUG
907 : // If we could not find the requested Hilbert index in
908 : // my_bin, something went terribly wrong, possibly the
909 : // Mesh was displaced differently on different processors,
910 : // and therefore the Hilbert indices don't agree.
911 1145788 : if (*pos != hilbert_indices)
912 : {
913 : // The input will be hilbert_indices. We convert it
914 : // to BitVecType using the operator= provided by the
915 : // BitVecType class. BitVecType is a CBigBitVec!
916 0 : Hilbert::BitVecType input;
917 0 : input = hilbert_indices.first;
918 :
919 : // Get output in a vector of CBigBitVec
920 0 : std::vector<CBigBitVec> output(3);
921 :
922 : // Call the indexToCoords function
923 0 : Hilbert::indexToCoords(output.data(), 8*sizeof(Hilbert::inttype), 3, input);
924 :
925 : // The entries in the output racks are integers in the
926 : // range [0, Hilbert::inttype::max] which can be
927 : // converted to floating point values in [0,1] and
928 : // finally to actual values using the bounding box.
929 0 : const Real max_int_as_real =
930 : static_cast<Real>(std::numeric_limits<Hilbert::inttype>::max());
931 :
932 : // Get the points in [0,1]^3. The zeroth rack of each entry in
933 : // 'output' maps to the normalized x, y, and z locations,
934 : // respectively.
935 0 : Point p_hat(static_cast<Real>(output[0].racks()[0]) / max_int_as_real,
936 0 : static_cast<Real>(output[1].racks()[0]) / max_int_as_real,
937 0 : static_cast<Real>(output[2].racks()[0]) / max_int_as_real);
938 :
939 : // Convert the points from [0,1]^3 to their actual (x,y,z) locations
940 : Real
941 0 : xmin = bbox.first(0),
942 0 : xmax = bbox.second(0),
943 0 : ymin = bbox.first(1),
944 0 : ymax = bbox.second(1),
945 0 : zmin = bbox.first(2),
946 0 : zmax = bbox.second(2);
947 :
948 : // Convert the points from [0,1]^3 to their actual (x,y,z) locations
949 0 : Point p(xmin + (xmax-xmin)*p_hat(0),
950 0 : ymin + (ymax-ymin)*p_hat(1),
951 0 : zmin + (zmax-zmin)*p_hat(2));
952 :
953 0 : libmesh_error_msg("Could not find hilbert indices: "
954 : << hilbert_indices
955 : << " corresponding to point " << p);
956 : }
957 : #endif
958 :
959 : // Finally, assign the global index based off the position of the index
960 : // in my array, properly offset.
961 5634980 : global_ids.push_back (cast_int<dof_id_type>(std::distance(my_bin.begin(), pos) + my_offset));
962 : }
963 : };
964 :
965 82319 : auto action_functor =
966 46388 : [&filled_request]
967 : (processor_id_type pid,
968 : const std::vector<Parallel::DofObjectKey> &,
969 23572 : const std::vector<dof_id_type> & global_ids)
970 : {
971 70338 : filled_request[pid] = global_ids;
972 : };
973 :
974 11787 : const dof_id_type * ex = nullptr;
975 : Parallel::pull_parallel_vector_data
976 35175 : (communicator, requested_ids, gather_functor, action_functor, ex);
977 :
978 : // We now have all the filled requests, so we can loop through our
979 : // nodes once and assign the global index to each one.
980 : {
981 : std::vector<std::vector<dof_id_type>::const_iterator>
982 46962 : next_obj_on_proc; next_obj_on_proc.reserve(communicator.size());
983 105516 : for (auto pid : make_range(communicator.size()))
984 117109 : next_obj_on_proc.push_back(filled_request[pid].begin());
985 :
986 11787 : unsigned int cnt=0;
987 4524367 : for (ForwardIterator it = begin; it != end; ++it, cnt++)
988 : {
989 1145788 : const processor_id_type pid = cast_int<processor_id_type>
990 3343404 : (index_map[cnt]);
991 :
992 1145788 : libmesh_assert_less (pid, communicator.size());
993 1145788 : libmesh_assert (next_obj_on_proc[pid] != filled_request[pid].end());
994 :
995 3343404 : const dof_id_type global_index = *next_obj_on_proc[pid];
996 3343404 : index_map[cnt] = global_index;
997 :
998 1145788 : ++next_obj_on_proc[pid];
999 : }
1000 : }
1001 : }
1002 :
1003 11787 : libmesh_assert_equal_to(index_map.size(), n_objects);
1004 46776 : }
1005 : #else // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
1006 : template <typename ForwardIterator>
1007 : void MeshCommunication::find_local_indices (const libMesh::BoundingBox &,
1008 : const ForwardIterator & begin,
1009 : const ForwardIterator & end,
1010 : std::unordered_map<dof_id_type, dof_id_type> & index_map) const
1011 : {
1012 : index_map.clear();
1013 : dof_id_type index = 0;
1014 : for (ForwardIterator it=begin; it!=end; ++it)
1015 : index_map[(*it)->id()] = (index++);
1016 : }
1017 :
1018 : template <typename ForwardIterator>
1019 : void MeshCommunication::find_global_indices (const Parallel::Communicator &,
1020 : const libMesh::BoundingBox &,
1021 : const ForwardIterator & begin,
1022 : const ForwardIterator & end,
1023 : std::vector<dof_id_type> & index_map) const
1024 : {
1025 : index_map.clear();
1026 : index_map.reserve(std::distance (begin, end));
1027 : unsigned int index = 0;
1028 : for (ForwardIterator it=begin; it!=end; ++it)
1029 : index_map.push_back(index++);
1030 : }
1031 : #endif // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
1032 :
1033 :
1034 :
1035 : //------------------------------------------------------------------
1036 : template LIBMESH_EXPORT void MeshCommunication::find_global_indices<MeshBase::const_node_iterator> (const Parallel::Communicator &,
1037 : const libMesh::BoundingBox &,
1038 : const MeshBase::const_node_iterator &,
1039 : const MeshBase::const_node_iterator &,
1040 : std::vector<dof_id_type> &) const;
1041 :
1042 : template LIBMESH_EXPORT void MeshCommunication::find_global_indices<MeshBase::const_element_iterator> (const Parallel::Communicator &,
1043 : const libMesh::BoundingBox &,
1044 : const MeshBase::const_element_iterator &,
1045 : const MeshBase::const_element_iterator &,
1046 : std::vector<dof_id_type> &) const;
1047 : template LIBMESH_EXPORT void MeshCommunication::find_global_indices<MeshBase::node_iterator> (const Parallel::Communicator &,
1048 : const libMesh::BoundingBox &,
1049 : const MeshBase::node_iterator &,
1050 : const MeshBase::node_iterator &,
1051 : std::vector<dof_id_type> &) const;
1052 :
1053 : template LIBMESH_EXPORT void MeshCommunication::find_global_indices<MeshBase::element_iterator> (const Parallel::Communicator &,
1054 : const libMesh::BoundingBox &,
1055 : const MeshBase::element_iterator &,
1056 : const MeshBase::element_iterator &,
1057 : std::vector<dof_id_type> &) const;
1058 : template LIBMESH_EXPORT void MeshCommunication::find_local_indices<MeshBase::const_element_iterator> (const libMesh::BoundingBox &,
1059 : const MeshBase::const_element_iterator &,
1060 : const MeshBase::const_element_iterator &,
1061 : std::unordered_map<dof_id_type, dof_id_type> &) const;
1062 :
1063 : } // namespace libMesh
|
