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 : #include "libmesh/partitioner.h"
21 :
22 : // libMesh includes
23 : #include "libmesh/compare_elems_by_level.h"
24 : #include "libmesh/elem.h"
25 : #include "libmesh/enum_to_string.h"
26 : #include "libmesh/int_range.h"
27 : #include "libmesh/libmesh_logging.h"
28 : #include "libmesh/mesh_base.h"
29 : #include "libmesh/mesh_tools.h"
30 : #include "libmesh/mesh_communication.h"
31 : #include "libmesh/node_range.h"
32 : #include "libmesh/parallel_ghost_sync.h"
33 : #include "libmesh/wrapped_petsc.h"
34 : #include "libmesh/boundary_info.h"
35 :
36 : // Subclasses to build()
37 : #include "libmesh/enum_partitioner_type.h"
38 : #include "libmesh/centroid_partitioner.h"
39 : #include "libmesh/hilbert_sfc_partitioner.h"
40 : #include "libmesh/linear_partitioner.h"
41 : #include "libmesh/mapped_subdomain_partitioner.h"
42 : #include "libmesh/metis_partitioner.h"
43 : #include "libmesh/morton_sfc_partitioner.h"
44 : #include "libmesh/parmetis_partitioner.h"
45 : #include "libmesh/sfc_partitioner.h"
46 : #include "libmesh/subdomain_partitioner.h"
47 :
48 : // TIMPI includes
49 : #include "timpi/parallel_implementation.h"
50 : #include "timpi/parallel_sync.h"
51 :
52 : // C/C++ includes
53 : #ifdef LIBMESH_HAVE_PETSC
54 : #include "libmesh/ignore_warnings.h"
55 : #include "petscmat.h"
56 : #include "libmesh/restore_warnings.h"
57 : #include "libmesh/petsc_solver_exception.h"
58 : #endif
59 :
60 :
61 : namespace {
62 :
63 : using namespace libMesh;
64 :
65 : struct CorrectProcIds
66 : {
67 : typedef processor_id_type datum;
68 :
69 146 : CorrectProcIds(MeshBase & _mesh,
70 445 : std::unordered_set<dof_id_type> & _bad_pids) :
71 445 : mesh(_mesh), bad_pids(_bad_pids) {}
72 :
73 : MeshBase & mesh;
74 : std::unordered_set<dof_id_type> & bad_pids;
75 :
76 : // ------------------------------------------------------------
77 936 : void gather_data (const std::vector<dof_id_type> & ids,
78 : std::vector<datum> & data)
79 : {
80 : // Find the processor id of each requested node
81 936 : data.resize(ids.size());
82 :
83 172794 : for (auto i : index_range(ids))
84 : {
85 : // Look for this point in the mesh and make sure we have a
86 : // good pid for it before sending that on
87 257787 : if (ids[i] != DofObject::invalid_id &&
88 171858 : !bad_pids.count(ids[i]))
89 : {
90 171858 : Node & node = mesh.node_ref(ids[i]);
91 :
92 : // Return the node's correct processor id,
93 171858 : data[i] = node.processor_id();
94 : }
95 : else
96 0 : data[i] = DofObject::invalid_processor_id;
97 : }
98 936 : }
99 :
100 : // ------------------------------------------------------------
101 936 : bool act_on_data (const std::vector<dof_id_type> & ids,
102 : const std::vector<datum> proc_ids)
103 : {
104 468 : bool data_changed = false;
105 :
106 : // Set the ghost node processor ids we've now been informed of
107 172794 : for (auto i : index_range(ids))
108 : {
109 171858 : Node & node = mesh.node_ref(ids[i]);
110 :
111 85929 : processor_id_type & proc_id = node.processor_id();
112 171858 : const processor_id_type new_proc_id = proc_ids[i];
113 :
114 171858 : if (const auto it = bad_pids.find(ids[i]);
115 171858 : (new_proc_id != proc_id || it != bad_pids.end()) &&
116 : new_proc_id != DofObject::invalid_processor_id)
117 : {
118 212 : if (it != bad_pids.end())
119 : {
120 212 : proc_id = new_proc_id;
121 106 : bad_pids.erase(it);
122 : }
123 : else
124 0 : proc_id = node.choose_processor_id(proc_id, new_proc_id);
125 :
126 212 : if (proc_id == new_proc_id)
127 106 : data_changed = true;
128 : }
129 : }
130 :
131 936 : return data_changed;
132 : }
133 : };
134 :
135 : }
136 :
137 : namespace libMesh
138 : {
139 :
140 :
141 :
142 : // ------------------------------------------------------------
143 : // Partitioner static data
144 : const dof_id_type Partitioner::communication_blocksize =
145 : dof_id_type(1000000);
146 :
147 :
148 :
149 : // ------------------------------------------------------------
150 : // Partitioner implementation
151 :
152 :
153 0 : PartitionerType Partitioner::type() const
154 : {
155 0 : return INVALID_PARTITIONER;
156 : }
157 :
158 :
159 : std::unique_ptr<Partitioner>
160 48045 : Partitioner::build (const PartitionerType partitioner_type)
161 : {
162 48045 : switch (partitioner_type)
163 : {
164 7 : case CENTROID_PARTITIONER:
165 7 : return std::make_unique<CentroidPartitioner>();
166 21 : case LINEAR_PARTITIONER:
167 21 : return std::make_unique<LinearPartitioner>();
168 0 : case MAPPED_SUBDOMAIN_PARTITIONER:
169 0 : return std::make_unique<MappedSubdomainPartitioner>();
170 47071 : case METIS_PARTITIONER:
171 47071 : return std::make_unique<MetisPartitioner>();
172 925 : case PARMETIS_PARTITIONER:
173 925 : return std::make_unique<ParmetisPartitioner>();
174 7 : case HILBERT_SFC_PARTITIONER:
175 7 : return std::make_unique<HilbertSFCPartitioner>();
176 7 : case MORTON_SFC_PARTITIONER:
177 7 : return std::make_unique<MortonSFCPartitioner>();
178 7 : case SFC_PARTITIONER:
179 7 : return std::make_unique<SFCPartitioner>();
180 0 : case SUBDOMAIN_PARTITIONER:
181 0 : return std::make_unique<SubdomainPartitioner>();
182 0 : default:
183 0 : libmesh_error_msg("Invalid partitioner type: " <<
184 : Utility::enum_to_string(partitioner_type));
185 : }
186 : }
187 :
188 :
189 :
190 0 : void Partitioner::partition (MeshBase & mesh)
191 : {
192 0 : this->partition(mesh,mesh.n_processors());
193 0 : }
194 :
195 :
196 :
197 60163 : void Partitioner::partition (MeshBase & mesh,
198 : const unsigned int n)
199 : {
200 17565 : libmesh_parallel_only(mesh.comm());
201 :
202 : // BSK - temporary fix while redistribution is integrated 6/26/2008
203 : // Uncomment this to not repartition in parallel
204 : // if (!mesh.is_serial())
205 : // return;
206 :
207 : // we cannot partition into more pieces than we have
208 : // active elements!
209 : const unsigned int n_parts =
210 : static_cast<unsigned int>
211 60163 : (std::min(mesh.n_active_elem(), static_cast<dof_id_type>(n)));
212 :
213 : // Set the number of partitions in the mesh
214 60163 : mesh.set_n_partitions()=n_parts;
215 :
216 60163 : if (n_parts == 1)
217 : {
218 24957 : this->single_partition (mesh);
219 24957 : return;
220 : }
221 :
222 : // First assign a temporary partitioning to any unpartitioned elements. A spatially
223 : // coherent assignment is only worth its parallel sort for a partitioner that refines what
224 : // it is handed; the other partitioners just need one that is consistent across processors,
225 : // which the cheap traversal-order assignment already is on a replicated mesh.
226 : Partitioner::partition_unpartitioned_elements
227 35206 : (mesh, n_parts,
228 35206 : /*spatial_indexing=*/!(mesh.is_replicated() && this->partitions_from_scratch()));
229 :
230 : // Call the partitioning function
231 35206 : this->_do_partition(mesh,n_parts);
232 :
233 : // Set the parent's processor ids
234 35206 : Partitioner::set_parent_processor_ids(mesh);
235 :
236 : // Ghosting functors may iterate over local elements while we
237 : // redistribute, so cached local ranges must reflect the new element
238 : // processor ids
239 35206 : mesh.update_post_partitioning();
240 :
241 : // Redistribute elements if necessary, before setting node processor
242 : // ids, to make sure those will be set consistently
243 35206 : mesh.redistribute();
244 :
245 : #ifdef DEBUG
246 11823 : MeshTools::libmesh_assert_valid_remote_elems(mesh);
247 :
248 : // Messed up elem processor_id()s can leave us without the child
249 : // elements we need to restrict vectors on a distributed mesh
250 11823 : MeshTools::libmesh_assert_valid_procids<Elem>(mesh);
251 : #endif
252 :
253 : // Set the node's processor ids
254 35206 : Partitioner::set_node_processor_ids(mesh);
255 :
256 : #ifdef DEBUG
257 11823 : MeshTools::libmesh_assert_valid_procids<Elem>(mesh);
258 : #endif
259 :
260 : // Give derived Mesh classes a chance to update any cached data to
261 : // reflect the new partitioning
262 35206 : mesh.update_post_partitioning();
263 : }
264 :
265 :
266 :
267 0 : void Partitioner::repartition (MeshBase & mesh)
268 : {
269 0 : this->repartition(mesh,mesh.n_processors());
270 0 : }
271 :
272 :
273 :
274 0 : void Partitioner::repartition (MeshBase & mesh,
275 : const unsigned int n)
276 : {
277 : // we cannot partition into more pieces than we have
278 : // active elements!
279 : const unsigned int n_parts =
280 : static_cast<unsigned int>
281 0 : (std::min(mesh.n_active_elem(), static_cast<dof_id_type>(n)));
282 :
283 : // Set the number of partitions in the mesh
284 0 : mesh.set_n_partitions()=n_parts;
285 :
286 0 : if (n_parts == 1)
287 : {
288 0 : this->single_partition (mesh);
289 0 : return;
290 : }
291 :
292 : // First assign a temporary partitioning to any unpartitioned elements; see partition()
293 : // for why the indexing is chosen this way.
294 : Partitioner::partition_unpartitioned_elements
295 0 : (mesh, n_parts,
296 0 : /*spatial_indexing=*/!(mesh.is_replicated() && this->partitions_from_scratch()));
297 :
298 : // Call the partitioning function
299 0 : this->_do_repartition(mesh,n_parts);
300 :
301 : // Set the parent's processor ids
302 0 : Partitioner::set_parent_processor_ids(mesh);
303 :
304 : // Set the node's processor ids
305 0 : Partitioner::set_node_processor_ids(mesh);
306 :
307 : // Give derived Mesh classes a chance to update any cached data to
308 : // reflect the new partitioning
309 0 : mesh.update_post_partitioning();
310 : }
311 :
312 :
313 :
314 :
315 :
316 24957 : bool Partitioner::single_partition (MeshBase & mesh)
317 : {
318 : bool changed_pid =
319 44172 : this->single_partition_range(mesh.elements_begin(),
320 63387 : mesh.elements_end());
321 :
322 : // If we have a distributed mesh with an empty rank (or where rank
323 : // 0 has only its own component of a disconnected mesh, I guess),
324 : // that rank might need to be informed of a change.
325 24957 : mesh.comm().max(changed_pid);
326 :
327 24957 : if (changed_pid)
328 : {
329 : // Ghosting functors may iterate over local elements while we
330 : // redistribute, so cached local ranges must reflect the new
331 : // element processor ids
332 22834 : mesh.update_post_partitioning();
333 :
334 : // We may need to redistribute, in case someone (like our unit
335 : // tests) is doing something silly (like moving a whole
336 : // already-distributed mesh back onto rank 0).
337 22834 : mesh.redistribute();
338 :
339 : // Even single_partition can be a change, if we're starting with
340 : // any unpartitioned elements or nodes
341 22834 : mesh.update_post_partitioning();
342 : }
343 :
344 24957 : return changed_pid;
345 : }
346 :
347 :
348 :
349 24957 : bool Partitioner::single_partition_range (MeshBase::element_iterator it,
350 : MeshBase::element_iterator end)
351 : {
352 11484 : LOG_SCOPE("single_partition_range()", "Partitioner");
353 :
354 5742 : bool changed_pid = false;
355 :
356 1200886 : for (auto & elem : as_range(it, end))
357 : {
358 590018 : if (elem->processor_id())
359 10342 : changed_pid = true;
360 590018 : elem->processor_id() = 0;
361 :
362 : // Assign all this element's nodes to processor 0 as well.
363 5135063 : for (Node & node : elem->node_ref_range())
364 : {
365 4545045 : if (node.processor_id())
366 52439 : changed_pid = true;
367 4545045 : node.processor_id() = 0;
368 : }
369 13473 : }
370 :
371 30699 : return changed_pid;
372 : }
373 :
374 0 : void Partitioner::partition_unpartitioned_elements (MeshBase & mesh)
375 : {
376 0 : Partitioner::partition_unpartitioned_elements(mesh, mesh.n_processors());
377 0 : }
378 :
379 :
380 :
381 0 : void Partitioner::partition_unpartitioned_elements (MeshBase & mesh,
382 : const unsigned int n_subdomains)
383 : {
384 0 : Partitioner::partition_unpartitioned_elements(mesh, n_subdomains,
385 : /*spatial_indexing=*/true);
386 0 : }
387 :
388 :
389 :
390 35206 : void Partitioner::partition_unpartitioned_elements (MeshBase & mesh,
391 : const unsigned int n_subdomains,
392 : const bool spatial_indexing)
393 : {
394 35206 : MeshBase::element_iterator it = mesh.unpartitioned_elements_begin();
395 35206 : const MeshBase::element_iterator end = mesh.unpartitioned_elements_end();
396 :
397 58589 : const dof_id_type n_unpartitioned_elements = MeshTools::n_elem (it, end);
398 :
399 : // the unpartitioned elements must exist on all processors. If the range is empty on one
400 : // it is empty on all, and we can quit right here.
401 35206 : if (!n_unpartitioned_elements)
402 3787 : return;
403 :
404 : // find the target subdomain sizes
405 40144 : std::vector<dof_id_type> subdomain_bounds(mesh.n_processors());
406 :
407 72216 : for (auto pid : make_range(mesh.n_processors()))
408 : {
409 16072 : dof_id_type tgt_subdomain_size = 0;
410 :
411 : // watch out for the case that n_subdomains < n_processors
412 48144 : if (pid < n_subdomains)
413 : {
414 48144 : tgt_subdomain_size = n_unpartitioned_elements/n_subdomains;
415 :
416 48144 : if (pid < n_unpartitioned_elements%n_subdomains)
417 1426 : tgt_subdomain_size++;
418 :
419 : }
420 :
421 : //libMesh::out << "pid, #= " << pid << ", " << tgt_subdomain_size << std::endl;
422 48144 : if (pid == 0)
423 24072 : subdomain_bounds[0] = tgt_subdomain_size;
424 : else
425 40144 : subdomain_bounds[pid] = subdomain_bounds[pid-1] + tgt_subdomain_size;
426 : }
427 :
428 8036 : libmesh_assert_equal_to (subdomain_bounds.back(), n_unpartitioned_elements);
429 :
430 : // create the unique mapping for all unpartitioned elements independent of partitioning
431 : // determine the global indexing for all the unpartitioned elements
432 16072 : std::vector<dof_id_type> global_indices;
433 :
434 : // Calling this on all processors a unique range in [0,n_unpartitioned_elements) is constructed.
435 : // Only the indices for the elements we pass in are returned in the array.
436 : // Without the spatial indexing an element's position in the traversal of the range is used
437 : // instead, which a replicated mesh reaches in the same order on every processor and which
438 : // therefore indexes [0, n_unpartitioned_elements) just as uniquely, at no cost.
439 24072 : if (spatial_indexing)
440 666 : MeshCommunication().find_global_indices (mesh.comm(),
441 1332 : MeshTools::create_bounding_box(mesh), it, end,
442 : global_indices);
443 : else
444 7814 : libmesh_assert (mesh.is_replicated());
445 :
446 8036 : dof_id_type cnt=0;
447 1193074 : for (auto & elem : as_range(it, end))
448 : {
449 377882 : libmesh_assert (!spatial_indexing || cnt < global_indices.size());
450 : const dof_id_type global_index =
451 1152966 : spatial_indexing ? global_indices[cnt++] : cnt++;
452 :
453 377882 : libmesh_assert_less (global_index, subdomain_bounds.back());
454 377882 : libmesh_assert_less (global_index, n_unpartitioned_elements);
455 :
456 : const processor_id_type subdomain_id =
457 : cast_int<processor_id_type>
458 775084 : (std::distance(subdomain_bounds.begin(),
459 : std::upper_bound(subdomain_bounds.begin(),
460 : subdomain_bounds.end(),
461 377882 : global_index)));
462 377882 : libmesh_assert_less (subdomain_id, n_subdomains);
463 :
464 1152966 : elem->processor_id() = subdomain_id;
465 : //libMesh::out << "assigning " << global_index << " to " << subdomain_id << std::endl;
466 8000 : }
467 :
468 : // Even a temporary partitioning may invalidate some mesh caches
469 24072 : mesh.update_post_partitioning();
470 : }
471 :
472 :
473 :
474 35206 : void Partitioner::set_parent_processor_ids(MeshBase & mesh)
475 : {
476 : // Ignore the parameter when !LIBMESH_ENABLE_AMR
477 11823 : libmesh_ignore(mesh);
478 :
479 35469 : LOG_SCOPE("set_parent_processor_ids()", "Partitioner");
480 :
481 : #ifdef LIBMESH_ENABLE_AMR
482 :
483 : // If the mesh is serial we have access to all the elements,
484 : // in particular all the active ones. We can therefore set
485 : // the parent processor ids indirectly through their children, and
486 : // set the subactive processor ids while examining their active
487 : // ancestors.
488 : // By convention a parent is assigned to the minimum processor
489 : // of all its children, and a subactive is assigned to the processor
490 : // of its active ancestor.
491 35206 : if (mesh.is_serial())
492 : {
493 4418277 : for (auto & elem : mesh.active_element_ptr_range())
494 : {
495 : // First set descendents
496 2279056 : std::vector<Elem *> subactive_family;
497 3319702 : elem->total_family_tree(subactive_family);
498 6771932 : for (const auto & f : subactive_family)
499 3452230 : f->processor_id() = elem->processor_id();
500 :
501 : // Then set ancestors
502 3319702 : Elem * parent = elem->parent();
503 :
504 8752966 : while (parent)
505 : {
506 : // invalidate the parent id, otherwise the min below
507 : // will not work if the current parent id is less
508 : // than all the children!
509 1988386 : parent->invalidate_processor_id();
510 :
511 31257528 : for (auto & child : parent->child_ref_range())
512 : {
513 9497036 : libmesh_assert(!child.is_remote());
514 9497036 : libmesh_assert_not_equal_to (child.processor_id(), DofObject::invalid_processor_id);
515 25824264 : parent->processor_id() = std::min(parent->processor_id(),
516 9497036 : child.processor_id());
517 : }
518 3976772 : parent = parent->parent();
519 : }
520 11428 : }
521 : }
522 :
523 : // When the mesh is parallel we cannot guarantee that parents have access to
524 : // all their children.
525 : else
526 : {
527 : // Setting subactive processor ids is easy: we can guarantee
528 : // that children have access to all their parents.
529 :
530 : // Loop over all the active elements in the mesh
531 46234 : for (auto & child : mesh.active_element_ptr_range())
532 : {
533 22788 : std::vector<Elem *> subactive_family;
534 34181 : child->total_family_tree(subactive_family);
535 68362 : for (const auto & f : subactive_family)
536 34181 : f->processor_id() = child->processor_id();
537 132 : }
538 :
539 : // When the mesh is parallel we cannot guarantee that parents have access to
540 : // all their children.
541 :
542 : // We will use a brute-force approach here. Each processor finds its parent
543 : // elements and sets the parent pid to the minimum of its
544 : // semilocal descendants.
545 : // A global reduction is then performed to make sure the true minimum is found.
546 : // As noted, this is required because we cannot guarantee that a parent has
547 : // access to all its children on any single processor.
548 132 : libmesh_parallel_only(mesh.comm());
549 132 : libmesh_assert(MeshTools::n_elem(mesh.unpartitioned_elements_begin(),
550 : mesh.unpartitioned_elements_end()) == 0);
551 :
552 396 : const dof_id_type max_elem_id = mesh.max_elem_id();
553 :
554 : std::vector<processor_id_type>
555 264 : parent_processor_ids (std::min(communication_blocksize,
556 528 : max_elem_id));
557 :
558 792 : for (dof_id_type blk=0, last_elem_id=0; last_elem_id<max_elem_id; blk++)
559 : {
560 396 : last_elem_id =
561 528 : std::min(static_cast<dof_id_type>((blk+1)*communication_blocksize),
562 132 : max_elem_id);
563 396 : const dof_id_type first_elem_id = blk*communication_blocksize;
564 :
565 132 : std::fill (parent_processor_ids.begin(),
566 : parent_processor_ids.end(),
567 : DofObject::invalid_processor_id);
568 :
569 : // first build up local contributions to parent_processor_ids
570 132 : bool have_parent_in_block = false;
571 :
572 1652 : for (auto & parent : as_range(mesh.ancestor_elements_begin(),
573 6280 : mesh.ancestor_elements_end()))
574 : {
575 2976 : const dof_id_type parent_idx = parent->id();
576 992 : libmesh_assert_less (parent_idx, max_elem_id);
577 :
578 3968 : if ((parent_idx >= first_elem_id) &&
579 2976 : (parent_idx < last_elem_id))
580 : {
581 992 : have_parent_in_block = true;
582 2976 : processor_id_type parent_pid = DofObject::invalid_processor_id;
583 :
584 992 : std::vector<const Elem *> active_family;
585 2976 : parent->active_family_tree(active_family);
586 29112 : for (const auto & f : active_family)
587 28376 : parent_pid = std::min (parent_pid, f->processor_id());
588 :
589 2976 : const dof_id_type packed_idx = parent_idx - first_elem_id;
590 992 : libmesh_assert_less (packed_idx, parent_processor_ids.size());
591 :
592 3968 : parent_processor_ids[packed_idx] = parent_pid;
593 : }
594 132 : }
595 :
596 : // then find the global minimum
597 396 : mesh.comm().min (parent_processor_ids);
598 :
599 : // and assign the ids, if we have a parent in this block.
600 396 : if (have_parent_in_block)
601 1252 : for (auto & parent : as_range(mesh.ancestor_elements_begin(),
602 5480 : mesh.ancestor_elements_end()))
603 : {
604 2976 : const dof_id_type parent_idx = parent->id();
605 :
606 3968 : if ((parent_idx >= first_elem_id) &&
607 1984 : (parent_idx < last_elem_id))
608 : {
609 2976 : const dof_id_type packed_idx = parent_idx - first_elem_id;
610 992 : libmesh_assert_less (packed_idx, parent_processor_ids.size());
611 :
612 : const processor_id_type parent_pid =
613 2976 : parent_processor_ids[packed_idx];
614 :
615 992 : libmesh_assert_not_equal_to (parent_pid, DofObject::invalid_processor_id);
616 :
617 2976 : parent->processor_id() = parent_pid;
618 : }
619 52 : }
620 : }
621 : }
622 :
623 : #endif // LIBMESH_ENABLE_AMR
624 35206 : }
625 :
626 : void
627 0 : Partitioner::processor_pairs_to_interface_nodes(MeshBase & mesh,
628 : std::map<std::pair<processor_id_type, processor_id_type>, std::set<dof_id_type>> & processor_pair_to_nodes)
629 : {
630 : // This function must be run on all processors at once
631 0 : libmesh_parallel_only(mesh.comm());
632 :
633 0 : processor_pair_to_nodes.clear();
634 :
635 0 : std::set<dof_id_type> mynodes;
636 0 : std::set<dof_id_type> neighbor_nodes;
637 0 : std::vector<dof_id_type> common_nodes;
638 :
639 : // Loop over all the active elements
640 0 : for (auto & elem : mesh.active_element_ptr_range())
641 : {
642 0 : libmesh_assert(elem);
643 :
644 0 : libmesh_assert_not_equal_to (elem->processor_id(), DofObject::invalid_processor_id);
645 :
646 0 : auto n_nodes = elem->n_nodes();
647 :
648 : // prepare data for this element
649 0 : mynodes.clear();
650 0 : neighbor_nodes.clear();
651 0 : common_nodes.clear();
652 :
653 0 : for (unsigned int inode = 0; inode < n_nodes; inode++)
654 0 : mynodes.insert(elem->node_id(inode));
655 :
656 0 : for (auto i : elem->side_index_range())
657 : {
658 0 : auto neigh = elem->neighbor_ptr(i);
659 0 : if (neigh && !neigh->is_remote() && neigh->processor_id() != elem->processor_id())
660 : {
661 0 : neighbor_nodes.clear();
662 0 : common_nodes.clear();
663 0 : auto neigh_n_nodes = neigh->n_nodes();
664 0 : for (unsigned int inode = 0; inode < neigh_n_nodes; inode++)
665 0 : neighbor_nodes.insert(neigh->node_id(inode));
666 :
667 0 : std::set_intersection(mynodes.begin(), mynodes.end(),
668 : neighbor_nodes.begin(), neighbor_nodes.end(),
669 0 : std::back_inserter(common_nodes));
670 :
671 0 : auto & map_set = processor_pair_to_nodes[std::make_pair(std::min(elem->processor_id(), neigh->processor_id()),
672 0 : std::max(elem->processor_id(), neigh->processor_id()))];
673 0 : for (auto global_node_id : common_nodes)
674 0 : map_set.insert(global_node_id);
675 : }
676 : }
677 0 : }
678 0 : }
679 :
680 0 : void Partitioner::set_interface_node_processor_ids_linear(MeshBase & mesh)
681 : {
682 : // This function must be run on all processors at once
683 0 : libmesh_parallel_only(mesh.comm());
684 :
685 0 : std::map<std::pair<processor_id_type, processor_id_type>, std::set<dof_id_type>> processor_pair_to_nodes;
686 :
687 0 : processor_pairs_to_interface_nodes(mesh, processor_pair_to_nodes);
688 :
689 0 : for (auto & pmap : processor_pair_to_nodes)
690 : {
691 0 : std::size_t n_own_nodes = pmap.second.size()/2, i = 0;
692 :
693 0 : for (dof_id_type id : pmap.second)
694 : {
695 0 : auto & node = mesh.node_ref(id);
696 0 : if (i <= n_own_nodes)
697 0 : node.processor_id() = pmap.first.first;
698 : else
699 0 : node.processor_id() = pmap.first.second;
700 0 : i++;
701 : }
702 : }
703 0 : }
704 :
705 0 : void Partitioner::set_interface_node_processor_ids_BFS(MeshBase & mesh)
706 : {
707 : // This function must be run on all processors at once
708 0 : libmesh_parallel_only(mesh.comm());
709 :
710 : // I see occasional consistency failures when using this on a
711 : // distributed mesh
712 : libmesh_experimental();
713 :
714 0 : std::map<std::pair<processor_id_type, processor_id_type>, std::set<dof_id_type>> processor_pair_to_nodes;
715 :
716 0 : processor_pairs_to_interface_nodes(mesh, processor_pair_to_nodes);
717 :
718 0 : std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
719 :
720 0 : MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
721 :
722 0 : std::vector<const Node *> neighbors;
723 0 : std::set<dof_id_type> neighbors_order;
724 0 : std::vector<dof_id_type> common_nodes;
725 0 : std::queue<dof_id_type> nodes_queue;
726 0 : std::set<dof_id_type> visted_nodes;
727 :
728 0 : for (auto & pmap : processor_pair_to_nodes)
729 : {
730 0 : std::size_t n_own_nodes = pmap.second.size()/2;
731 :
732 : // Initialize node assignment
733 0 : for (dof_id_type id : pmap.second)
734 0 : mesh.node_ref(id).processor_id() = pmap.first.second;
735 :
736 0 : visted_nodes.clear();
737 0 : for (dof_id_type id : pmap.second)
738 : {
739 0 : mesh.node_ref(id).processor_id() = pmap.first.second;
740 :
741 0 : if (visted_nodes.count(id))
742 0 : continue;
743 : else
744 : {
745 0 : nodes_queue.push(id);
746 0 : visted_nodes.insert(id);
747 0 : if (visted_nodes.size() >= n_own_nodes)
748 0 : break;
749 : }
750 :
751 0 : while (!nodes_queue.empty())
752 : {
753 0 : auto & node = mesh.node_ref(nodes_queue.front());
754 0 : nodes_queue.pop();
755 :
756 0 : neighbors.clear();
757 0 : MeshTools::find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);
758 0 : neighbors_order.clear();
759 0 : for (auto & neighbor : neighbors)
760 0 : neighbors_order.insert(neighbor->id());
761 :
762 0 : common_nodes.clear();
763 0 : std::set_intersection(pmap.second.begin(), pmap.second.end(),
764 : neighbors_order.begin(), neighbors_order.end(),
765 0 : std::back_inserter(common_nodes));
766 :
767 0 : for (auto c_node : common_nodes)
768 0 : if (!visted_nodes.count(c_node))
769 : {
770 0 : nodes_queue.push(c_node);
771 0 : visted_nodes.insert(c_node);
772 0 : if (visted_nodes.size() >= n_own_nodes)
773 0 : goto queue_done;
774 : }
775 :
776 0 : if (visted_nodes.size() >= n_own_nodes)
777 0 : goto queue_done;
778 : }
779 : }
780 0 : queue_done:
781 0 : for (auto node : visted_nodes)
782 0 : mesh.node_ref(node).processor_id() = pmap.first.first;
783 : }
784 0 : }
785 :
786 0 : void Partitioner::set_interface_node_processor_ids_petscpartitioner(MeshBase & mesh)
787 : {
788 0 : libmesh_ignore(mesh); // Only used if LIBMESH_HAVE_PETSC
789 :
790 : // This function must be run on all processors at once
791 0 : libmesh_parallel_only(mesh.comm());
792 :
793 : #ifdef LIBMESH_HAVE_PETSC
794 0 : std::map<std::pair<processor_id_type, processor_id_type>, std::set<dof_id_type>> processor_pair_to_nodes;
795 :
796 0 : processor_pairs_to_interface_nodes(mesh, processor_pair_to_nodes);
797 :
798 0 : std::vector<std::vector<const Elem *>> nodes_to_elem_map;
799 :
800 0 : MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
801 :
802 0 : std::vector<const Node *> neighbors;
803 0 : std::set<dof_id_type> neighbors_order;
804 0 : std::vector<dof_id_type> common_nodes;
805 :
806 0 : std::vector<dof_id_type> rows;
807 0 : std::vector<dof_id_type> cols;
808 :
809 0 : std::map<dof_id_type, dof_id_type> global_to_local;
810 :
811 0 : for (auto & pmap : processor_pair_to_nodes)
812 : {
813 0 : unsigned int i = 0;
814 :
815 0 : rows.clear();
816 0 : rows.resize(pmap.second.size()+1);
817 0 : cols.clear();
818 0 : for (dof_id_type id : pmap.second)
819 0 : global_to_local[id] = i++;
820 :
821 0 : i = 0;
822 0 : for (auto id : pmap.second)
823 : {
824 0 : auto & node = mesh.node_ref(id);
825 0 : neighbors.clear();
826 0 : MeshTools::find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);
827 0 : neighbors_order.clear();
828 0 : for (auto & neighbor : neighbors)
829 0 : neighbors_order.insert(neighbor->id());
830 :
831 0 : common_nodes.clear();
832 0 : std::set_intersection(pmap.second.begin(), pmap.second.end(),
833 : neighbors_order.begin(), neighbors_order.end(),
834 0 : std::back_inserter(common_nodes));
835 :
836 0 : rows[i+1] = rows[i] + cast_int<dof_id_type>(common_nodes.size());
837 :
838 0 : for (auto c_node : common_nodes)
839 0 : cols.push_back(global_to_local[c_node]);
840 :
841 0 : i++;
842 : }
843 :
844 : // Next we construct an IS from a MatPartitioning
845 0 : WrappedPetsc<IS> is;
846 : {
847 : PetscInt *adj_i, *adj_j;
848 0 : LibmeshPetscCall2(mesh.comm(), PetscCalloc1(rows.size(), &adj_i));
849 0 : LibmeshPetscCall2(mesh.comm(), PetscCalloc1(cols.size(), &adj_j));
850 0 : PetscInt rows_size = cast_int<PetscInt>(rows.size());
851 0 : for (PetscInt ii=0; ii<rows_size; ii++)
852 0 : adj_i[ii] = rows[ii];
853 :
854 0 : PetscInt cols_size = cast_int<PetscInt>(cols.size());
855 0 : for (PetscInt ii=0; ii<cols_size; ii++)
856 0 : adj_j[ii] = cols[ii];
857 :
858 0 : const PetscInt sz = cast_int<PetscInt>(pmap.second.size());
859 :
860 : // Create sparse matrix representing an adjacency list
861 0 : WrappedPetsc<Mat> adj;
862 0 : LibmeshPetscCall2(mesh.comm(), MatCreateMPIAdj(PETSC_COMM_SELF, sz, sz, adj_i, adj_j, nullptr, adj.get()));
863 :
864 : // Create MatPartitioning object
865 0 : WrappedPetsc<MatPartitioning> part;
866 0 : LibmeshPetscCall2(mesh.comm(), MatPartitioningCreate(PETSC_COMM_SELF, part.get()));
867 :
868 : // Apply MatPartitioning, storing results in "is"
869 0 : LibmeshPetscCall2(mesh.comm(), MatPartitioningSetAdjacency(part, adj));
870 0 : LibmeshPetscCall2(mesh.comm(), MatPartitioningSetNParts(part, 2));
871 0 : LibmeshPetscCall2(mesh.comm(), PetscObjectSetOptionsPrefix((PetscObject)(*part), "balance_"));
872 0 : LibmeshPetscCall2(mesh.comm(), MatPartitioningSetFromOptions(part));
873 0 : LibmeshPetscCall2(mesh.comm(), MatPartitioningApply(part, is.get()));
874 : }
875 :
876 : PetscInt local_size;
877 : const PetscInt *indices;
878 0 : LibmeshPetscCall2(mesh.comm(), ISGetLocalSize(is, &local_size));
879 0 : LibmeshPetscCall2(mesh.comm(), ISGetIndices(is, &indices));
880 :
881 0 : i = 0;
882 0 : for (auto id : pmap.second)
883 : {
884 0 : auto & node = mesh.node_ref(id);
885 0 : if (indices[i])
886 0 : node.processor_id() = pmap.first.second;
887 : else
888 0 : node.processor_id() = pmap.first.first;
889 :
890 0 : i++;
891 : }
892 0 : LibmeshPetscCall2(mesh.comm(), ISRestoreIndices(is, &indices));
893 : }
894 : #else
895 0 : libmesh_error_msg("PETSc is required");
896 : #endif
897 0 : }
898 :
899 :
900 41697 : void Partitioner::set_node_processor_ids(MeshBase & mesh)
901 : {
902 41421 : LOG_SCOPE("set_node_processor_ids()","Partitioner");
903 :
904 : // This function must be run on all processors at once
905 13807 : libmesh_parallel_only(mesh.comm());
906 :
907 : // If we have any unpartitioned elements at this
908 : // stage there is a problem
909 13807 : libmesh_assert (MeshTools::n_elem(mesh.unpartitioned_elements_begin(),
910 : mesh.unpartitioned_elements_end()) == 0);
911 :
912 : // Start from scratch here: nodes we used to own may not be
913 : // eligible for us to own any more.
914 : Threads::parallel_for_each
915 41697 : (mesh.node_stored_range(), [](Node * node)
916 8263281 : { node->processor_id() = DofObject::invalid_processor_id; });
917 :
918 : // Loop over all the active elements
919 7232059 : for (auto & elem : mesh.active_element_ptr_range())
920 : {
921 1785462 : libmesh_assert(elem);
922 :
923 1785462 : libmesh_assert_not_equal_to (elem->processor_id(), DofObject::invalid_processor_id);
924 :
925 : // Consider updating the processor id on this element's nodes
926 46067732 : for (Node & node : elem->node_ref_range())
927 : {
928 12721665 : processor_id_type & pid = node.processor_id();
929 38915572 : pid = node.choose_processor_id(pid, elem->processor_id());
930 : }
931 14083 : }
932 :
933 : // How we finish off the node partitioning depends on our command
934 : // line options.
935 :
936 : const bool load_balanced_nodes_linear =
937 41697 : libMesh::on_command_line ("--load-balanced-nodes-linear");
938 :
939 : const bool load_balanced_nodes_bfs =
940 41697 : libMesh::on_command_line ("--load-balanced-nodes-bfs");
941 :
942 : const bool load_balanced_nodes_petscpartition =
943 41697 : libMesh::on_command_line ("--load-balanced-nodes-petscpartitioner");
944 :
945 41697 : unsigned int n_load_balance_options = load_balanced_nodes_linear;
946 41697 : n_load_balance_options += load_balanced_nodes_bfs;
947 41697 : n_load_balance_options += load_balanced_nodes_petscpartition;
948 41697 : libmesh_error_msg_if(n_load_balance_options > 1,
949 : "Cannot perform more than one load balancing type at a time");
950 :
951 41697 : if (load_balanced_nodes_linear)
952 0 : set_interface_node_processor_ids_linear(mesh);
953 41697 : else if (load_balanced_nodes_bfs)
954 0 : set_interface_node_processor_ids_BFS(mesh);
955 41697 : else if (load_balanced_nodes_petscpartition)
956 0 : set_interface_node_processor_ids_petscpartitioner(mesh);
957 :
958 : // Node balancing algorithm will response to assign owned nodes.
959 : // We still need to sync PIDs
960 : {
961 : // For inactive elements, we will have already gotten most of
962 : // these nodes, *except* for the case of a parent with a subset
963 : // of active descendants which are remote elements. In that
964 : // case some of the parent nodes will not have been properly
965 : // handled yet on our processor.
966 : //
967 : // We don't want to inadvertently give one of them an incorrect
968 : // processor id, but if we're not in serial then we have to
969 : // assign them temporary pids to make querying work, so we'll
970 : // save our *valid* pids before assigning temporaries.
971 : //
972 : // Even in serial we'll want to check and make sure we're not
973 : // overwriting valid active node pids with pids from subactive
974 : // elements.
975 27614 : std::unordered_set<dof_id_type> bad_pids;
976 :
977 14865397 : for (auto & node : mesh.node_ptr_range())
978 10960591 : if (node->processor_id() == DofObject::invalid_processor_id)
979 187913 : bad_pids.insert(node->id());
980 :
981 : // If we assign our temporary ids by looping from finer elements
982 : // to coarser elements, we'll always get an id from the finest
983 : // ghost element we can see, which will usually be "closer" to
984 : // the true processor we want to query and so will reduce query
985 : // cycles that don't reach that processor.
986 :
987 : // But we can still end up with a query cycle that dead-ends, so
988 : // we need to prepare a "push" communication step here.
989 :
990 41697 : const bool is_serial = mesh.is_serial();
991 : std::unordered_map
992 : <processor_id_type,
993 : std::unordered_map<dof_id_type, processor_id_type>>
994 27614 : potential_pids;
995 :
996 41697 : const unsigned int n_levels = MeshTools::n_levels(mesh);
997 109180 : for (unsigned int level = n_levels; level > 0; --level)
998 : {
999 2287467 : for (auto & elem : as_range(mesh.level_elements_begin(level-1),
1000 11089890 : mesh.level_elements_end(level-1)))
1001 : {
1002 2174918 : libmesh_assert_not_equal_to (elem->processor_id(),
1003 : DofObject::invalid_processor_id);
1004 :
1005 6519855 : const processor_id_type elem_pid = elem->processor_id();
1006 :
1007 : // Consider updating the processor id on this element's nodes
1008 51263988 : for (Node & node : elem->node_ref_range())
1009 : {
1010 14666607 : processor_id_type & pid = node.processor_id();
1011 44744133 : if (bad_pids.count(node.id()))
1012 522398 : pid = node.choose_processor_id(pid, elem_pid);
1013 44221735 : else if (!is_serial)
1014 369605 : potential_pids[elem_pid][node.id()] = pid;
1015 : }
1016 22649 : }
1017 : }
1018 :
1019 41697 : if (!is_serial)
1020 : {
1021 : std::unordered_map
1022 : <processor_id_type,
1023 : std::vector<std::pair<dof_id_type, processor_id_type>>>
1024 146 : potential_pids_vecs;
1025 :
1026 1268 : for (auto & pair : potential_pids)
1027 823 : potential_pids_vecs[pair.first].assign(pair.second.begin(), pair.second.end());
1028 :
1029 : auto pids_action_functor =
1030 279 : [& mesh, & bad_pids]
1031 : (processor_id_type /* src_pid */,
1032 211622 : const std::vector<std::pair<dof_id_type, processor_id_type>> & data)
1033 : {
1034 106520 : for (auto pair : data)
1035 : {
1036 105697 : Node & node = mesh.node_ref(pair.first);
1037 35209 : processor_id_type & pid = node.processor_id();
1038 105697 : if (const auto it = bad_pids.find(pair.first);
1039 35209 : it != bad_pids.end())
1040 : {
1041 78 : pid = pair.second;
1042 26 : bad_pids.erase(it);
1043 : }
1044 : else
1045 105619 : pid = node.choose_processor_id(pid, pair.second);
1046 : }
1047 843 : };
1048 :
1049 : Parallel::push_parallel_vector_data
1050 445 : (mesh.comm(), potential_pids_vecs, pids_action_functor);
1051 :
1052 : // Using default libMesh options, we'll just need to sync
1053 : // between processors now. The catch here is that we can't
1054 : // initially trust Node::choose_processor_id() because some
1055 : // of those node processor ids are the temporary ones.
1056 146 : CorrectProcIds correct_pids(mesh, bad_pids);
1057 : Parallel::sync_node_data_by_element_id
1058 1043 : (mesh, mesh.elements_begin(), mesh.elements_end(),
1059 299 : Parallel::SyncEverything(), Parallel::SyncEverything(),
1060 : correct_pids);
1061 :
1062 : // But once we've got all the non-temporary pids synced, we
1063 : // may need to sync again to get any pids on nodes only
1064 : // connected to subactive elements, for which *only*
1065 : // "temporary" pids are possible.
1066 146 : bad_pids.clear();
1067 : Parallel::sync_node_data_by_element_id
1068 890 : (mesh,
1069 1481 : mesh.elements_begin(), mesh.elements_end(),
1070 445 : Parallel::SyncEverything(), Parallel::SyncEverything(),
1071 : correct_pids);
1072 : }
1073 : }
1074 :
1075 : // We can't assert that all nodes are connected to elements, because
1076 : // a DistributedMesh with NodeConstraints might have pulled in some
1077 : // remote nodes solely for evaluating those constraints.
1078 : // MeshTools::libmesh_assert_connected_nodes(mesh);
1079 :
1080 : #ifdef DEBUG
1081 13807 : MeshTools::libmesh_assert_valid_procids<Node>(mesh);
1082 : //MeshTools::libmesh_assert_canonical_node_procids(mesh);
1083 : #endif
1084 41697 : }
1085 :
1086 :
1087 :
1088 : struct SyncLocalIDs
1089 : {
1090 : typedef dof_id_type datum;
1091 :
1092 : typedef std::unordered_map<dof_id_type, dof_id_type> map_type;
1093 :
1094 2340 : SyncLocalIDs(map_type & _id_map) : id_map(_id_map) {}
1095 :
1096 : map_type & id_map;
1097 :
1098 2298 : void gather_data (const std::vector<dof_id_type> & ids,
1099 : std::vector<datum> & local_ids) const
1100 : {
1101 3064 : local_ids.resize(ids.size());
1102 :
1103 50056 : for (auto i : index_range(ids))
1104 47758 : local_ids[i] = id_map[ids[i]];
1105 2298 : }
1106 :
1107 2298 : void act_on_data (const std::vector<dof_id_type> & ids,
1108 : const std::vector<datum> & local_ids)
1109 : {
1110 50056 : for (auto i : index_range(local_ids))
1111 63678 : id_map[ids[i]] = local_ids[i];
1112 2298 : }
1113 : };
1114 :
1115 2340 : void Partitioner::_find_global_index_by_pid_map(const MeshBase & mesh)
1116 : {
1117 2340 : const dof_id_type n_active_local_elem = mesh.n_active_local_elem();
1118 :
1119 : // Find the number of active elements on each processor. We cannot use
1120 : // mesh.n_active_elem_on_proc(pid) since that only returns the number of
1121 : // elements assigned to pid which are currently stored on the calling
1122 : // processor. This will not in general be correct for parallel meshes
1123 : // when (pid!=mesh.processor_id()).
1124 1560 : auto n_proc = mesh.n_processors();
1125 2340 : _n_active_elem_on_proc.resize(n_proc);
1126 2340 : mesh.comm().allgather(n_active_local_elem, _n_active_elem_on_proc);
1127 :
1128 3900 : std::vector<dof_id_type> n_active_elem_before_proc(mesh.n_processors());
1129 :
1130 4680 : for (auto i : make_range(n_proc-1))
1131 2340 : n_active_elem_before_proc[i+1] =
1132 3900 : n_active_elem_before_proc[i] + _n_active_elem_on_proc[i];
1133 :
1134 : libMesh::BoundingBox bbox =
1135 2340 : MeshTools::create_bounding_box(mesh);
1136 :
1137 780 : _global_index_by_pid_map.clear();
1138 :
1139 : // create the mapping which is contiguous by processor
1140 2340 : MeshCommunication().find_local_indices (bbox,
1141 4680 : mesh.active_local_elements_begin(),
1142 3900 : mesh.active_local_elements_end(),
1143 2340 : _global_index_by_pid_map);
1144 :
1145 780 : SyncLocalIDs sync(_global_index_by_pid_map);
1146 :
1147 : Parallel::sync_dofobject_data_by_id
1148 3900 : (mesh.comm(), mesh.active_elements_begin(), mesh.active_elements_end(), sync);
1149 :
1150 122662 : for (const auto & elem : mesh.active_element_ptr_range())
1151 : {
1152 118762 : const processor_id_type pid = elem->processor_id();
1153 39588 : libmesh_assert_less (_global_index_by_pid_map[elem->id()], _n_active_elem_on_proc[pid]);
1154 :
1155 158350 : _global_index_by_pid_map[elem->id()] += n_active_elem_before_proc[pid];
1156 780 : }
1157 2340 : }
1158 :
1159 1170 : void Partitioner::build_graph (const MeshBase & mesh)
1160 : {
1161 1170 : LOG_SCOPE("build_graph()", "Partitioner");
1162 :
1163 390 : const dof_id_type n_active_local_elem = mesh.n_active_local_elem();
1164 :
1165 : // If we have boundary elements in this mesh, we want to account for
1166 : // the connectivity between them and interior elements. We can find
1167 : // interior elements from boundary elements, but we need to build up
1168 : // a lookup map to do the reverse.
1169 : typedef std::unordered_multimap<const Elem *, const Elem *> map_type;
1170 780 : map_type interior_to_boundary_map;
1171 :
1172 : // If we have spline nodes in this mesh, we want to account for the
1173 : // connectivity between them and integration elements. We can find
1174 : // spline nodes from integration elements, but need a reverse map
1175 : // {integration_elements} = elems_constrained_by[spline_nodeelem]
1176 780 : map_type elems_constrained_by;
1177 :
1178 390 : const auto & mesh_constrained_nodes = mesh.get_constraint_rows();
1179 :
1180 120711 : for (const Elem * elem : mesh.active_element_ptr_range())
1181 : {
1182 59381 : if (!mesh_constrained_nodes.empty()) // quick test for non-IGA cases
1183 : {
1184 0 : const auto end_it = mesh_constrained_nodes.end();
1185 :
1186 : // Use a set to avoid duplicates. Use a well-defined method
1187 : // of ordering that set to make debugging easier.
1188 0 : std::set<const Elem *, CompareElemIdsByLevel> constraining_elems;
1189 0 : for (const Node & node : elem->node_ref_range())
1190 : {
1191 0 : if (const auto row_it = mesh_constrained_nodes.find(&node);
1192 0 : row_it != end_it)
1193 0 : for (const auto & [pr, coef] : row_it->second)
1194 : {
1195 0 : libmesh_ignore(coef); // avoid gcc 7 warning
1196 0 : constraining_elems.insert(pr.first);
1197 : }
1198 : }
1199 0 : for (const Elem * constraining_elem : constraining_elems)
1200 0 : elems_constrained_by.emplace(constraining_elem, elem);
1201 : }
1202 :
1203 : // If we don't have an interior_parent and we don't have any
1204 : // constrained nodes then there's nothing else to look up.
1205 59381 : if (elem->interior_parent())
1206 : {
1207 : // get all relevant interior elements
1208 0 : std::set<const Elem *> neighbor_set;
1209 0 : elem->find_interior_neighbors(neighbor_set);
1210 :
1211 0 : for (const auto & neighbor : neighbor_set)
1212 0 : interior_to_boundary_map.emplace(neighbor, elem);
1213 : }
1214 390 : }
1215 :
1216 : #ifdef LIBMESH_ENABLE_AMR
1217 780 : std::vector<const Elem *> neighbors_offspring;
1218 : #endif
1219 :
1220 : // This is costly, and we only need to do it if the mesh has
1221 : // changed since we last partitioned... but the mesh probably has
1222 : // changed since we last partitioned, and if it hasn't we don't
1223 : // have a reliable way to be sure of that.
1224 1170 : _find_global_index_by_pid_map(mesh);
1225 :
1226 390 : dof_id_type first_local_elem = 0;
1227 1755 : for (auto pid : make_range(mesh.processor_id()))
1228 780 : first_local_elem += _n_active_elem_on_proc[pid];
1229 :
1230 780 : _dual_graph.clear();
1231 1170 : _dual_graph.resize(n_active_local_elem);
1232 1170 : _local_id_to_elem.resize(n_active_local_elem);
1233 :
1234 : // We may need to communicate constraint-row-based connections
1235 : // between processors
1236 : std::unordered_map<processor_id_type,
1237 780 : std::vector<std::pair<dof_id_type, dof_id_type>>> connections_to_push;
1238 :
1239 49286 : for (const auto & elem : mesh.active_local_element_ptr_range())
1240 : {
1241 11834 : libmesh_assert (_global_index_by_pid_map.count(elem->id()));
1242 : const dof_id_type global_index_by_pid =
1243 35502 : _global_index_by_pid_map[elem->id()];
1244 :
1245 35502 : const dof_id_type local_index =
1246 11834 : global_index_by_pid - first_local_elem;
1247 11834 : libmesh_assert_less (local_index, n_active_local_elem);
1248 :
1249 23668 : std::vector<dof_id_type> & graph_row = _dual_graph[local_index];
1250 :
1251 : // Save this off to make it easy to index later
1252 35502 : _local_id_to_elem[local_index] = const_cast<Elem*>(elem);
1253 :
1254 : // Loop over the element's neighbors. An element
1255 : // adjacency corresponds to a face neighbor
1256 216932 : for (auto neighbor : elem->neighbor_ptr_range())
1257 : {
1258 169596 : if (neighbor != nullptr)
1259 : {
1260 : // If the neighbor is active treat it
1261 : // as a connection
1262 50897 : if (neighbor->active())
1263 : {
1264 50770 : libmesh_assert(_global_index_by_pid_map.count(neighbor->id()));
1265 : const dof_id_type neighbor_global_index_by_pid =
1266 152310 : _global_index_by_pid_map[neighbor->id()];
1267 :
1268 152310 : graph_row.push_back(neighbor_global_index_by_pid);
1269 : }
1270 :
1271 : #ifdef LIBMESH_ENABLE_AMR
1272 :
1273 : // Otherwise we need to find all of the
1274 : // neighbor's children that are connected to
1275 : // us and add them
1276 : else
1277 : {
1278 : // The side of the neighbor to which
1279 : // we are connected
1280 : const unsigned int ns =
1281 381 : neighbor->which_neighbor_am_i (elem);
1282 127 : libmesh_assert_less (ns, neighbor->n_neighbors());
1283 :
1284 : // Get all the active children (& grandchildren, etc...)
1285 : // of the neighbor
1286 :
1287 : // FIXME - this is the wrong thing, since we
1288 : // should be getting the active family tree on
1289 : // our side only. But adding too many graph
1290 : // links may cause hanging nodes to tend to be
1291 : // on partition interiors, which would reduce
1292 : // communication overhead for constraint
1293 : // equations, so we'll leave it.
1294 :
1295 381 : neighbor->active_family_tree (neighbors_offspring);
1296 :
1297 : // Get all the neighbor's children that
1298 : // live on that side and are thus connected
1299 : // to us
1300 2148 : for (const auto & child : neighbors_offspring)
1301 : {
1302 : // This does not assume a level-1 mesh.
1303 : // Note that since children have sides numbered
1304 : // coincident with the parent then this is a sufficient test.
1305 2356 : if (child->neighbor_ptr(ns) == elem)
1306 : {
1307 254 : libmesh_assert (child->active());
1308 254 : libmesh_assert (_global_index_by_pid_map.count(child->id()));
1309 : const dof_id_type child_global_index_by_pid =
1310 762 : _global_index_by_pid_map[child->id()];
1311 :
1312 762 : graph_row.push_back(child_global_index_by_pid);
1313 : }
1314 : }
1315 : }
1316 :
1317 : #endif /* ifdef LIBMESH_ENABLE_AMR */
1318 :
1319 :
1320 : }
1321 : }
1322 :
1323 45885 : if ((elem->dim() < LIBMESH_DIM) &&
1324 10383 : elem->interior_parent())
1325 : {
1326 : // get all relevant interior elements
1327 0 : std::set<const Elem *> neighbor_set;
1328 0 : elem->find_interior_neighbors(neighbor_set);
1329 :
1330 0 : for (const auto & neighbor : neighbor_set)
1331 : {
1332 : const dof_id_type neighbor_global_index_by_pid =
1333 0 : _global_index_by_pid_map[neighbor->id()];
1334 :
1335 0 : graph_row.push_back(neighbor_global_index_by_pid);
1336 : }
1337 : }
1338 :
1339 : // Check for any boundary neighbors
1340 35502 : for (const auto & pr : as_range(interior_to_boundary_map.equal_range(elem)))
1341 : {
1342 0 : const Elem * neighbor = pr.second;
1343 :
1344 : const dof_id_type neighbor_global_index_by_pid =
1345 0 : _global_index_by_pid_map[neighbor->id()];
1346 :
1347 0 : graph_row.push_back(neighbor_global_index_by_pid);
1348 : }
1349 :
1350 : // Check for any constraining elements
1351 35502 : if (!mesh_constrained_nodes.empty()) // quick test for non-IGA cases
1352 : {
1353 0 : const auto end_it = mesh_constrained_nodes.end();
1354 :
1355 : // Use a set to avoid duplicates. Use a well-defined method
1356 : // of ordering that set to make debugging easier.
1357 0 : std::set<const Elem *, CompareElemIdsByLevel> constraining_elems;
1358 0 : for (const Node & node : elem->node_ref_range())
1359 : {
1360 0 : if (const auto row_it = mesh_constrained_nodes.find(&node);
1361 0 : row_it != end_it)
1362 0 : for (const auto & [pr, coef] : row_it->second)
1363 : {
1364 0 : libmesh_ignore(coef); // avoid gcc 7 warning
1365 0 : constraining_elems.insert(pr.first);
1366 : }
1367 : }
1368 0 : for (const Elem * constraining_elem : constraining_elems)
1369 : {
1370 : const dof_id_type constraining_global_index_by_pid =
1371 0 : _global_index_by_pid_map[constraining_elem->id()];
1372 :
1373 0 : graph_row.push_back(constraining_global_index_by_pid);
1374 :
1375 : // We can't be sure if the constraining element's owner sees
1376 : // the assembly element, so to get a symmetric connectivity
1377 : // graph we'll need to tell them about us to be safe.
1378 0 : if (constraining_elem->processor_id() != mesh.processor_id())
1379 0 : connections_to_push[constraining_elem->processor_id()].emplace_back
1380 0 : (global_index_by_pid, constraining_global_index_by_pid);
1381 : }
1382 : }
1383 :
1384 : // Check for any constrained elements
1385 35502 : for (const auto & pr : as_range(elems_constrained_by.equal_range(elem)))
1386 : {
1387 0 : const Elem * constrained = pr.second;
1388 : const dof_id_type constrained_global_index_by_pid =
1389 0 : _global_index_by_pid_map[constrained->id()];
1390 :
1391 0 : graph_row.push_back(constrained_global_index_by_pid);
1392 :
1393 : // We can't be sure if the constrained element's owner sees
1394 : // the assembly element, so to get a symmetric connectivity
1395 : // graph we'll need to tell them about us to be safe.
1396 0 : if (constrained->processor_id() != mesh.processor_id())
1397 0 : connections_to_push[constrained->processor_id()].emplace_back
1398 0 : (global_index_by_pid, constrained_global_index_by_pid);
1399 : }
1400 390 : }
1401 :
1402 : // Partitioners like Parmetis require a symmetric adjacency matrix,
1403 : // but if we have an assembly element constrained by a spline
1404 : // NodeElem owned by another processor, it's possible that that
1405 : // processor doesn't see our assembly element. Let's push those
1406 : // entries, to ensure they're counted on both sides.
1407 : auto symmetrize_entries =
1408 0 : [this, first_local_elem]
1409 : (processor_id_type /*src_pid*/,
1410 0 : const std::vector<std::pair<dof_id_type, dof_id_type>> & incoming_entries)
1411 : {
1412 0 : for (auto [i, j] : incoming_entries)
1413 : {
1414 0 : libmesh_assert_greater_equal(j, first_local_elem);
1415 0 : const std::size_t jl = j - first_local_elem;
1416 0 : libmesh_assert_less(jl, _dual_graph.size());
1417 0 : std::vector<dof_id_type> & graph_row = _dual_graph[jl];
1418 0 : if (std::find(graph_row.begin(), graph_row.end(), i) == graph_row.end())
1419 : {
1420 : // std::cerr << "Pushing back (" << j << ", " << i << ") from " << src_pid << std::endl;
1421 0 : graph_row.push_back(i);
1422 : }
1423 : }
1424 780 : };
1425 :
1426 1170 : Parallel::push_parallel_vector_data(mesh.comm(),
1427 : connections_to_push,
1428 : symmetrize_entries);
1429 :
1430 : // *Now* we should have a symmetric adjacency matrix. Let's check
1431 : // that, so any failures get caught before they e.g. confuse
1432 : // Parmetis in hard-to-debug ways. That's a global communication so
1433 : // we'll only do it in debug mode.
1434 : #ifdef DEBUG
1435 390 : auto n_proc = mesh.n_processors();
1436 780 : std::vector<dof_id_type> first_local_index_on_proc(n_proc, 0);
1437 780 : for (auto pid : make_range(1u,n_proc))
1438 390 : first_local_index_on_proc[pid] = first_local_index_on_proc[pid-1] + _n_active_elem_on_proc[pid-1];
1439 :
1440 390 : libmesh_assert_equal_to(first_local_index_on_proc[mesh.processor_id()],
1441 : first_local_elem);
1442 :
1443 780 : std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, dof_id_type>>> entries_to_send;
1444 12224 : for (auto il : index_range(_dual_graph))
1445 : {
1446 11834 : const std::vector<dof_id_type> & graph_row = _dual_graph[il];
1447 :
1448 11834 : const auto i = il + first_local_elem;
1449 :
1450 62858 : for (auto j : graph_row)
1451 : {
1452 : // Stupid graph rows aren't sorted yet...
1453 51024 : processor_id_type target_pid = 0;
1454 127340 : while (target_pid+1 < n_proc &&
1455 51024 : j >= first_local_index_on_proc[target_pid+1])
1456 25292 : ++target_pid;
1457 :
1458 51024 : entries_to_send[target_pid].emplace_back(i,j);
1459 : }
1460 : }
1461 :
1462 780 : std::vector<std::tuple<processor_id_type, dof_id_type, dof_id_type>> bad_entries;
1463 :
1464 : auto check_incoming_entries =
1465 : [this, first_local_elem, &bad_entries]
1466 : (processor_id_type src_pid,
1467 204824 : const std::vector<std::pair<dof_id_type, dof_id_type>> & incoming_entries)
1468 : {
1469 51752 : for (auto [i, j] : incoming_entries)
1470 : {
1471 51024 : if (j < first_local_elem)
1472 : {
1473 0 : bad_entries.emplace_back(src_pid,i,j);
1474 0 : continue;
1475 : }
1476 102048 : const std::size_t jl = j - first_local_elem;
1477 51024 : if (jl >= _dual_graph.size())
1478 : {
1479 0 : bad_entries.emplace_back(src_pid,i,j);
1480 0 : continue;
1481 : }
1482 51024 : const std::vector<dof_id_type> & graph_row = _dual_graph[jl];
1483 51024 : if (std::find(graph_row.begin(), graph_row.end(), i) ==
1484 102048 : graph_row.end())
1485 : {
1486 0 : bad_entries.emplace_back(src_pid,i,j);
1487 0 : continue;
1488 : }
1489 : }
1490 728 : };
1491 :
1492 : // Keep any failures in sync in parallel, for easier debugging in
1493 : // unit tests where the failure exception will be caught.
1494 : Parallel::push_parallel_vector_data
1495 390 : (mesh.comm(), entries_to_send, check_incoming_entries);
1496 390 : bool bad_entries_exist = !bad_entries.empty();
1497 390 : mesh.comm().max(bad_entries_exist);
1498 390 : if (bad_entries_exist)
1499 : {
1500 : #if 0 // Optional verbosity for if this breaks again...
1501 : if (!bad_entries.empty())
1502 : {
1503 : std::cerr << "Bad entries on processor " << mesh.processor_id() << ": ";
1504 : for (auto [p, i, j] : bad_entries)
1505 : std::cerr << '(' << p << ", " << i << ", " << j << "), ";
1506 : std::cerr << std::endl;
1507 : }
1508 : #endif
1509 0 : libmesh_error_msg("Asymmetric partitioner graph detected");
1510 : }
1511 : #endif
1512 1170 : }
1513 :
1514 882 : void Partitioner::assign_partitioning (MeshBase & mesh, const std::vector<dof_id_type> & parts)
1515 : {
1516 882 : LOG_SCOPE("assign_partitioning()", "Partitioner");
1517 :
1518 : // This function must be run on all processors at once
1519 294 : libmesh_parallel_only(mesh.comm());
1520 :
1521 294 : dof_id_type first_local_elem = 0;
1522 1323 : for (auto pid : make_range(mesh.processor_id()))
1523 588 : first_local_elem += _n_active_elem_on_proc[pid];
1524 :
1525 : #ifndef NDEBUG
1526 294 : const dof_id_type n_active_local_elem = mesh.n_active_local_elem();
1527 : #endif
1528 :
1529 : std::map<processor_id_type, std::vector<dof_id_type>>
1530 588 : requested_ids;
1531 :
1532 : // Results to gather from each processor - kept in a map so we
1533 : // do only one loop over elements after all receives are done.
1534 : std::map<processor_id_type, std::vector<processor_id_type>>
1535 588 : filled_request;
1536 :
1537 59045 : for (auto & elem : mesh.active_element_ptr_range())
1538 : {
1539 : // we need to get the index from the owning processor
1540 : // (note we cannot assign it now -- we are iterating
1541 : // over elements again and this will be bad!)
1542 76767 : requested_ids[elem->processor_id()].push_back(elem->id());
1543 294 : }
1544 :
1545 : auto gather_functor =
1546 588 : [this,
1547 : & parts,
1548 : #ifndef NDEBUG
1549 : & mesh,
1550 : n_active_local_elem,
1551 : #endif
1552 : first_local_elem]
1553 : (processor_id_type, const std::vector<dof_id_type> & ids,
1554 212286 : std::vector<processor_id_type> & data)
1555 : {
1556 1176 : const std::size_t ids_size = ids.size();
1557 1764 : data.resize(ids.size());
1558 :
1559 59339 : for (std::size_t i=0; i != ids_size; i++)
1560 : {
1561 76767 : const dof_id_type requested_elem_index = ids[i];
1562 :
1563 19192 : libmesh_assert(_global_index_by_pid_map.count(requested_elem_index));
1564 :
1565 : const dof_id_type global_index_by_pid =
1566 57575 : _global_index_by_pid_map[requested_elem_index];
1567 :
1568 57575 : const dof_id_type local_index =
1569 19192 : global_index_by_pid - first_local_elem;
1570 :
1571 19192 : libmesh_assert_less (local_index, parts.size());
1572 19192 : libmesh_assert_less (local_index, n_active_local_elem);
1573 :
1574 : const processor_id_type elem_procid =
1575 57575 : cast_int<processor_id_type>(parts[local_index]);
1576 :
1577 19192 : libmesh_assert_less (elem_procid, mesh.n_partitions());
1578 :
1579 76767 : data[i] = elem_procid;
1580 : }
1581 2352 : };
1582 :
1583 : auto action_functor =
1584 588 : [&filled_request]
1585 : (processor_id_type pid,
1586 : const std::vector<dof_id_type> &,
1587 1176 : const std::vector<processor_id_type> & new_procids)
1588 : {
1589 1764 : filled_request[pid] = new_procids;
1590 1470 : };
1591 :
1592 : // Trade requests with other processors
1593 294 : const processor_id_type * ex = nullptr;
1594 : Parallel::pull_parallel_vector_data
1595 882 : (mesh.comm(), requested_ids, gather_functor, action_functor, ex);
1596 :
1597 : // and finally assign the partitioning.
1598 : // note we are iterating in exactly the same order
1599 : // used to build up the request, so we can expect the
1600 : // required entries to be in the proper sequence.
1601 1470 : std::vector<unsigned int> counters(mesh.n_processors(), 0);
1602 59045 : for (auto & elem : mesh.active_element_ptr_range())
1603 : {
1604 57575 : const processor_id_type current_pid = elem->processor_id();
1605 :
1606 19192 : libmesh_assert_less (counters[current_pid], requested_ids[current_pid].size());
1607 :
1608 : const processor_id_type elem_procid =
1609 57575 : filled_request[current_pid][counters[current_pid]++];
1610 :
1611 19192 : libmesh_assert_less (elem_procid, mesh.n_partitions());
1612 57575 : elem->processor_id() = elem_procid;
1613 294 : }
1614 882 : }
1615 :
1616 :
1617 : } // namespace libMesh
|
