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