libMesh
parmetis_partitioner.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // Local Includes
21 #include "libmesh/libmesh_config.h"
22 #include "libmesh/mesh_base.h"
23 #include "libmesh/parallel.h" // also includes mpi.h
24 #include "libmesh/mesh_serializer.h"
25 #include "libmesh/mesh_tools.h"
26 #include "libmesh/mesh_communication.h"
27 #include "libmesh/parmetis_partitioner.h"
28 #include "libmesh/metis_partitioner.h"
29 #include "libmesh/libmesh_logging.h"
30 #include "libmesh/elem.h"
31 #include "libmesh/parmetis_helper.h"
32 
33 // Include the ParMETIS header file.
34 #ifdef LIBMESH_HAVE_PARMETIS
35 
36 namespace Parmetis {
37 extern "C" {
38 # include "libmesh/ignore_warnings.h"
39 # include "parmetis.h"
40 # include "libmesh/restore_warnings.h"
41 }
42 }
43 
44 #endif
45 
46 
47 // C++ includes
48 #include <unordered_map>
49 
50 
51 namespace libMesh
52 {
53 
54 // Minimum elements on each processor required for us to choose
55 // Parmetis over Metis.
56 #ifdef LIBMESH_HAVE_PARMETIS
57 const unsigned int MIN_ELEM_PER_PROC = 4;
58 #endif
59 
60 // ------------------------------------------------------------
61 // ParmetisPartitioner implementation
62 ParmetisPartitioner::ParmetisPartitioner()
63 #ifdef LIBMESH_HAVE_PARMETIS
64  : _pmetis(new ParmetisHelper)
65 #endif
66 {}
67 
68 
69 
70 ParmetisPartitioner::~ParmetisPartitioner()
71 {
72 }
73 
74 
75 
76 void ParmetisPartitioner::_do_partition (MeshBase & mesh,
77  const unsigned int n_sbdmns)
78 {
79  this->_do_repartition (mesh, n_sbdmns);
80 }
81 
82 
83 
84 void ParmetisPartitioner::_do_repartition (MeshBase & mesh,
85  const unsigned int n_sbdmns)
86 {
87  libmesh_assert_greater (n_sbdmns, 0);
88 
89  // Check for an easy return
90  if (n_sbdmns == 1)
91  {
92  this->single_partition(mesh);
93  return;
94  }
95 
96  // This function must be run on all processors at once
97  libmesh_parallel_only(mesh.comm());
98 
99  // What to do if the Parmetis library IS NOT present
100 #ifndef LIBMESH_HAVE_PARMETIS
101 
102  libmesh_here();
103  libMesh::err << "ERROR: The library has been built without" << std::endl
104  << "Parmetis support. Using a Metis" << std::endl
105  << "partitioner instead!" << std::endl;
106 
107  MetisPartitioner mp;
108 
109  mp.partition (mesh, n_sbdmns);
110 
111  // What to do if the Parmetis library IS present
112 #else
113 
114  // Revert to METIS on one processor.
115  if (mesh.n_processors() == 1)
116  {
117  // Make sure the mesh knows it's serial
118  mesh.allgather();
119 
120  MetisPartitioner mp;
121  mp.partition (mesh, n_sbdmns);
122  return;
123  }
124 
125  LOG_SCOPE("repartition()", "ParmetisPartitioner");
126 
127  // Initialize the data structures required by ParMETIS
128  this->initialize (mesh, n_sbdmns);
129 
130  // Make sure all processors have enough active local elements.
131  // Parmetis tends to crash when it's given only a couple elements
132  // per partition.
133  {
134  bool all_have_enough_elements = true;
135  for (std::size_t pid=0; pid<_n_active_elem_on_proc.size(); pid++)
136  if (_n_active_elem_on_proc[pid] < MIN_ELEM_PER_PROC)
137  all_have_enough_elements = false;
138 
139  // Parmetis will not work unless each processor has some
140  // elements. Specifically, it will abort when passed a NULL
141  // partition array on *any* of the processors.
142  if (!all_have_enough_elements)
143  {
144  // FIXME: revert to METIS, although this requires a serial mesh
145  MeshSerializer serialize(mesh);
146  MetisPartitioner mp;
147  mp.partition (mesh, n_sbdmns);
148  return;
149  }
150  }
151 
152  // build the graph corresponding to the mesh
153  this->build_graph (mesh);
154 
155 
156  // Partition the graph
157  std::vector<Parmetis::idx_t> vsize(_pmetis->vwgt.size(), 1);
158  Parmetis::real_t itr = 1000000.0;
159  MPI_Comm mpi_comm = mesh.comm().get();
160 
161  // Call the ParMETIS adaptive repartitioning method. This respects the
162  // original partitioning when computing the new partitioning so as to
163  // minimize the required data redistribution.
164  Parmetis::ParMETIS_V3_AdaptiveRepart(_pmetis->vtxdist.empty() ? libmesh_nullptr : &_pmetis->vtxdist[0],
165  _pmetis->xadj.empty() ? libmesh_nullptr : &_pmetis->xadj[0],
166  _pmetis->adjncy.empty() ? libmesh_nullptr : &_pmetis->adjncy[0],
167  _pmetis->vwgt.empty() ? libmesh_nullptr : &_pmetis->vwgt[0],
168  vsize.empty() ? libmesh_nullptr : &vsize[0],
170  &_pmetis->wgtflag,
171  &_pmetis->numflag,
172  &_pmetis->ncon,
173  &_pmetis->nparts,
174  _pmetis->tpwgts.empty() ? libmesh_nullptr : &_pmetis->tpwgts[0],
175  _pmetis->ubvec.empty() ? libmesh_nullptr : &_pmetis->ubvec[0],
176  &itr,
177  &_pmetis->options[0],
178  &_pmetis->edgecut,
179  _pmetis->part.empty() ? libmesh_nullptr : &_pmetis->part[0],
180  &mpi_comm);
181 
182  // Assign the returned processor ids
183  this->assign_partitioning (mesh);
184 
185 #endif // #ifndef LIBMESH_HAVE_PARMETIS ... else ...
186 
187 }
188 
189 
190 
191 // Only need to compile these methods if ParMETIS is present
192 #ifdef LIBMESH_HAVE_PARMETIS
193 
195  const unsigned int n_sbdmns)
196 {
197  LOG_SCOPE("initialize()", "ParmetisPartitioner");
198 
199  const dof_id_type n_active_local_elem = mesh.n_active_local_elem();
200 
201  // Set parameters.
202  _pmetis->wgtflag = 2; // weights on vertices only
203  _pmetis->ncon = 1; // one weight per vertex
204  _pmetis->numflag = 0; // C-style 0-based numbering
205  _pmetis->nparts = static_cast<Parmetis::idx_t>(n_sbdmns); // number of subdomains to create
206  _pmetis->edgecut = 0; // the numbers of edges cut by the
207  // partition
208 
209  // Initialize data structures for ParMETIS
210  _pmetis->vtxdist.assign (mesh.n_processors()+1, 0);
211  _pmetis->tpwgts.assign (_pmetis->nparts, 1./_pmetis->nparts);
212  _pmetis->ubvec.assign (_pmetis->ncon, 1.05);
213  _pmetis->part.assign (n_active_local_elem, 0);
214  _pmetis->options.resize (5);
215  _pmetis->vwgt.resize (n_active_local_elem);
216 
217  // Set the options
218  _pmetis->options[0] = 1; // don't use default options
219  _pmetis->options[1] = 0; // default (level of timing)
220  _pmetis->options[2] = 15; // random seed (default)
221  _pmetis->options[3] = 2; // processor distribution and subdomain distribution are decoupled
222 
223  // Find the number of active elements on each processor. We cannot use
224  // mesh.n_active_elem_on_proc(pid) since that only returns the number of
225  // elements assigned to pid which are currently stored on the calling
226  // processor. This will not in general be correct for parallel meshes
227  // when (pid!=mesh.processor_id()).
228  _n_active_elem_on_proc.resize(mesh.n_processors());
229  mesh.comm().allgather(n_active_local_elem, _n_active_elem_on_proc);
230 
231  // count the total number of active elements in the mesh. Note we cannot
232  // use mesh.n_active_elem() in general since this only returns the number
233  // of active elements which are stored on the calling processor.
234  // We should not use n_active_elem for any allocation because that will
235  // be inherently unscalable, but it can be useful for libmesh_assertions.
236  dof_id_type n_active_elem=0;
237 
238  // Set up the vtxdist array. This will be the same on each processor.
239  // ***** Consult the Parmetis documentation. *****
240  libmesh_assert_equal_to (_pmetis->vtxdist.size(),
241  cast_int<std::size_t>(mesh.n_processors()+1));
242  libmesh_assert_equal_to (_pmetis->vtxdist[0], 0);
243 
244  for (processor_id_type pid=0; pid<mesh.n_processors(); pid++)
245  {
246  _pmetis->vtxdist[pid+1] = _pmetis->vtxdist[pid] + _n_active_elem_on_proc[pid];
247  n_active_elem += _n_active_elem_on_proc[pid];
248  }
249  libmesh_assert_equal_to (_pmetis->vtxdist.back(), static_cast<Parmetis::idx_t>(n_active_elem));
250 
251  // ParMetis expects the elements to be numbered in contiguous blocks
252  // by processor, i.e. [0, ne0), [ne0, ne0+ne1), ...
253  // Since we only partition active elements we should have no expectation
254  // that we currently have such a distribution. So we need to create it.
255  // Also, at the same time we are going to map all the active elements into a globally
256  // unique range [0,n_active_elem) which is *independent* of the current partitioning.
257  // This can be fed to ParMetis as the initial partitioning of the subdomains (decoupled
258  // from the partitioning of the objects themselves). This allows us to get the same
259  // resultant partitioning independent of the input partitioning.
260  libMesh::BoundingBox bbox =
262 
263  _global_index_by_pid_map.clear();
264 
265  // Maps active element ids into a contiguous range independent of partitioning.
266  // (only needs local scope)
267  vectormap<dof_id_type, dof_id_type> global_index_map;
268 
269  {
270  std::vector<dof_id_type> global_index;
271 
272  // create the mapping which is contiguous by processor
273  dof_id_type pid_offset=0;
274  for (processor_id_type pid=0; pid<mesh.n_processors(); pid++)
275  {
278 
279  // note that we may not have all (or any!) the active elements which belong on this processor,
280  // but by calling this on all processors a unique range in [0,_n_active_elem_on_proc[pid])
281  // is constructed. Only the indices for the elements we pass in are returned in the array.
283  bbox, it, end,
284  global_index);
285 
286  for (dof_id_type cnt=0; it != end; ++it)
287  {
288  const Elem * elem = *it;
289  // vectormap::count forces a sort, which is too expensive
290  // in a loop
291  // libmesh_assert (!_global_index_by_pid_map.count(elem->id()));
292  libmesh_assert_less (cnt, global_index.size());
293  libmesh_assert_less (global_index[cnt], _n_active_elem_on_proc[pid]);
294 
295  _global_index_by_pid_map.insert(std::make_pair(elem->id(), global_index[cnt++] + pid_offset));
296  }
297 
298  pid_offset += _n_active_elem_on_proc[pid];
299  }
300 
301  // create the unique mapping for all active elements independent of partitioning
302  {
305 
306  // Calling this on all processors a unique range in [0,n_active_elem) is constructed.
307  // Only the indices for the elements we pass in are returned in the array.
309  bbox, it, end,
310  global_index);
311 
312  for (dof_id_type cnt=0; it != end; ++it)
313  {
314  const Elem * elem = *it;
315  // vectormap::count forces a sort, which is too expensive
316  // in a loop
317  // libmesh_assert (!global_index_map.count(elem->id()));
318  libmesh_assert_less (cnt, global_index.size());
319  libmesh_assert_less (global_index[cnt], n_active_elem);
320 
321  global_index_map.insert(std::make_pair(elem->id(), global_index[cnt++]));
322  }
323  }
324  // really, shouldn't be close!
325  libmesh_assert_less_equal (global_index_map.size(), n_active_elem);
326  libmesh_assert_less_equal (_global_index_by_pid_map.size(), n_active_elem);
327 
328  // At this point the two maps should be the same size. If they are not
329  // then the number of active elements is not the same as the sum over all
330  // processors of the number of active elements per processor, which means
331  // there must be some unpartitioned objects out there.
332  if (global_index_map.size() != _global_index_by_pid_map.size())
333  libmesh_error_msg("ERROR: ParmetisPartitioner cannot handle unpartitioned objects!");
334  }
335 
336  // Finally, we need to initialize the vertex (partition) weights and the initial subdomain
337  // mapping. The subdomain mapping will be independent of the processor mapping, and is
338  // defined by a simple mapping of the global indices we just found.
339  {
340  std::vector<dof_id_type> subdomain_bounds(mesh.n_processors());
341 
342  const dof_id_type first_local_elem = _pmetis->vtxdist[mesh.processor_id()];
343 
344  for (processor_id_type pid=0; pid<mesh.n_processors(); pid++)
345  {
346  dof_id_type tgt_subdomain_size = 0;
347 
348  // watch out for the case that n_subdomains < n_processors
349  if (pid < static_cast<unsigned int>(_pmetis->nparts))
350  {
351  tgt_subdomain_size = n_active_elem/std::min
352  (cast_int<Parmetis::idx_t>(mesh.n_processors()), _pmetis->nparts);
353 
354  if (pid < n_active_elem%_pmetis->nparts)
355  tgt_subdomain_size++;
356  }
357  if (pid == 0)
358  subdomain_bounds[0] = tgt_subdomain_size;
359  else
360  subdomain_bounds[pid] = subdomain_bounds[pid-1] + tgt_subdomain_size;
361  }
362 
363  libmesh_assert_equal_to (subdomain_bounds.back(), n_active_elem);
364 
365  for (const auto & elem : mesh.active_local_element_ptr_range())
366  {
367  libmesh_assert (_global_index_by_pid_map.count(elem->id()));
368  const dof_id_type global_index_by_pid =
369  _global_index_by_pid_map[elem->id()];
370  libmesh_assert_less (global_index_by_pid, n_active_elem);
371 
372  const dof_id_type local_index =
373  global_index_by_pid - first_local_elem;
374 
375  libmesh_assert_less (local_index, n_active_local_elem);
376  libmesh_assert_less (local_index, _pmetis->vwgt.size());
377 
378  // TODO:[BSK] maybe there is a better weight?
379  _pmetis->vwgt[local_index] = elem->n_nodes();
380 
381  // find the subdomain this element belongs in
382  libmesh_assert (global_index_map.count(elem->id()));
383  const dof_id_type global_index =
384  global_index_map[elem->id()];
385 
386  libmesh_assert_less (global_index, subdomain_bounds.back());
387 
388  const unsigned int subdomain_id =
389  std::distance(subdomain_bounds.begin(),
390  std::lower_bound(subdomain_bounds.begin(),
391  subdomain_bounds.end(),
392  global_index));
393  libmesh_assert_less (subdomain_id, static_cast<unsigned int>(_pmetis->nparts));
394  libmesh_assert_less (local_index, _pmetis->part.size());
395 
396  _pmetis->part[local_index] = subdomain_id;
397  }
398  }
399 }
400 
401 
402 
403 void ParmetisPartitioner::build_graph (const MeshBase & mesh)
404 {
405  LOG_SCOPE("build_graph()", "ParmetisPartitioner");
406 
407  // build the graph in distributed CSR format. Note that
408  // the edges in the graph will correspond to
409  // face neighbors
410  const dof_id_type n_active_local_elem = mesh.n_active_local_elem();
411 
412  // If we have boundary elements in this mesh, we want to account for
413  // the connectivity between them and interior elements. We can find
414  // interior elements from boundary elements, but we need to build up
415  // a lookup map to do the reverse.
416 
417  typedef std::unordered_multimap<const Elem *, const Elem *> map_type;
418  map_type interior_to_boundary_map;
419 
420  for (const auto & elem : mesh.active_element_ptr_range())
421  {
422  // If we don't have an interior_parent then there's nothing to look us
423  // up.
424  if ((elem->dim() >= LIBMESH_DIM) ||
425  !elem->interior_parent())
426  continue;
427 
428  // get all relevant interior elements
429  std::set<const Elem *> neighbor_set;
430  elem->find_interior_neighbors(neighbor_set);
431 
432  std::set<const Elem *>::iterator n_it = neighbor_set.begin();
433  for (; n_it != neighbor_set.end(); ++n_it)
434  interior_to_boundary_map.insert(std::make_pair(*n_it, elem));
435  }
436 
437 #ifdef LIBMESH_ENABLE_AMR
438  std::vector<const Elem *> neighbors_offspring;
439 #endif
440 
441  std::vector<std::vector<dof_id_type>> graph(n_active_local_elem);
442  dof_id_type graph_size=0;
443 
444  const dof_id_type first_local_elem = _pmetis->vtxdist[mesh.processor_id()];
445 
446  for (const auto & elem : mesh.active_local_element_ptr_range())
447  {
448  libmesh_assert (_global_index_by_pid_map.count(elem->id()));
449  const dof_id_type global_index_by_pid =
450  _global_index_by_pid_map[elem->id()];
451 
452  const dof_id_type local_index =
453  global_index_by_pid - first_local_elem;
454  libmesh_assert_less (local_index, n_active_local_elem);
455 
456  std::vector<dof_id_type> & graph_row = graph[local_index];
457 
458  // Loop over the element's neighbors. An element
459  // adjacency corresponds to a face neighbor
460  for (auto neighbor : elem->neighbor_ptr_range())
461  {
462  if (neighbor != libmesh_nullptr)
463  {
464  // If the neighbor is active treat it
465  // as a connection
466  if (neighbor->active())
467  {
468  libmesh_assert(_global_index_by_pid_map.count(neighbor->id()));
469  const dof_id_type neighbor_global_index_by_pid =
470  _global_index_by_pid_map[neighbor->id()];
471 
472  graph_row.push_back(neighbor_global_index_by_pid);
473  graph_size++;
474  }
475 
476 #ifdef LIBMESH_ENABLE_AMR
477 
478  // Otherwise we need to find all of the
479  // neighbor's children that are connected to
480  // us and add them
481  else
482  {
483  // The side of the neighbor to which
484  // we are connected
485  const unsigned int ns =
486  neighbor->which_neighbor_am_i (elem);
487  libmesh_assert_less (ns, neighbor->n_neighbors());
488 
489  // Get all the active children (& grandchildren, etc...)
490  // of the neighbor
491 
492  // FIXME - this is the wrong thing, since we
493  // should be getting the active family tree on
494  // our side only. But adding too many graph
495  // links may cause hanging nodes to tend to be
496  // on partition interiors, which would reduce
497  // communication overhead for constraint
498  // equations, so we'll leave it.
499 
500  neighbor->active_family_tree (neighbors_offspring);
501 
502  // Get all the neighbor's children that
503  // live on that side and are thus connected
504  // to us
505  for (std::size_t nc=0; nc<neighbors_offspring.size(); nc++)
506  {
507  const Elem * child =
508  neighbors_offspring[nc];
509 
510  // This does not assume a level-1 mesh.
511  // Note that since children have sides numbered
512  // coincident with the parent then this is a sufficient test.
513  if (child->neighbor_ptr(ns) == elem)
514  {
515  libmesh_assert (child->active());
516  libmesh_assert (_global_index_by_pid_map.count(child->id()));
517  const dof_id_type child_global_index_by_pid =
518  _global_index_by_pid_map[child->id()];
519 
520  graph_row.push_back(child_global_index_by_pid);
521  graph_size++;
522  }
523  }
524  }
525 
526 #endif /* ifdef LIBMESH_ENABLE_AMR */
527 
528 
529  }
530  }
531 
532  if ((elem->dim() < LIBMESH_DIM) &&
533  elem->interior_parent())
534  {
535  // get all relevant interior elements
536  std::set<const Elem *> neighbor_set;
537  elem->find_interior_neighbors(neighbor_set);
538 
539  std::set<const Elem *>::iterator n_it = neighbor_set.begin();
540  for (; n_it != neighbor_set.end(); ++n_it)
541  {
542  // FIXME - non-const versions of the Elem set methods
543  // would be nice
544  Elem * neighbor = const_cast<Elem *>(*n_it);
545 
546  const dof_id_type neighbor_global_index_by_pid =
547  _global_index_by_pid_map[neighbor->id()];
548 
549  graph_row.push_back(neighbor_global_index_by_pid);
550  graph_size++;
551  }
552  }
553 
554  // Check for any boundary neighbors
555  for (const auto & pr : as_range(interior_to_boundary_map.equal_range(elem)))
556  {
557  const Elem * neighbor = pr.second;
558 
559  const dof_id_type neighbor_global_index_by_pid =
560  _global_index_by_pid_map[neighbor->id()];
561 
562  graph_row.push_back(neighbor_global_index_by_pid);
563  graph_size++;
564  }
565  }
566 
567  // Reserve space in the adjacency array
568  _pmetis->xadj.clear();
569  _pmetis->xadj.reserve (n_active_local_elem + 1);
570  _pmetis->adjncy.clear();
571  _pmetis->adjncy.reserve (graph_size);
572 
573  for (std::size_t r=0; r<graph.size(); r++)
574  {
575  _pmetis->xadj.push_back(_pmetis->adjncy.size());
576  std::vector<dof_id_type> graph_row; // build this empty
577  graph_row.swap(graph[r]); // this will deallocate at the end of scope
578  _pmetis->adjncy.insert(_pmetis->adjncy.end(),
579  graph_row.begin(),
580  graph_row.end());
581  }
582 
583  // The end of the adjacency array for the last elem
584  _pmetis->xadj.push_back(_pmetis->adjncy.size());
585 
586  libmesh_assert_equal_to (_pmetis->xadj.size(), n_active_local_elem+1);
587  libmesh_assert_equal_to (_pmetis->adjncy.size(), graph_size);
588 }
589 
590 
591 
592 void ParmetisPartitioner::assign_partitioning (MeshBase & mesh)
593 {
594  LOG_SCOPE("assign_partitioning()", "ParmetisPartitioner");
595 
596  // This function must be run on all processors at once
597  libmesh_parallel_only(mesh.comm());
598 
599  const dof_id_type
600  first_local_elem = _pmetis->vtxdist[mesh.processor_id()];
601 
602 #ifndef NDEBUG
603  const dof_id_type n_active_local_elem = mesh.n_active_local_elem();
604 #endif
605 
606  std::vector<std::vector<dof_id_type>>
607  requested_ids(mesh.n_processors()),
608  requests_to_fill(mesh.n_processors());
609 
610  for (auto & elem : mesh.active_element_ptr_range())
611  {
612  // we need to get the index from the owning processor
613  // (note we cannot assign it now -- we are iterating
614  // over elements again and this will be bad!)
615  libmesh_assert_less (elem->processor_id(), requested_ids.size());
616  requested_ids[elem->processor_id()].push_back(elem->id());
617  }
618 
619  // Trade with all processors (including self) to get their indices
620  for (processor_id_type pid=0; pid<mesh.n_processors(); pid++)
621  {
622  // Trade my requests with processor procup and procdown
623  const processor_id_type procup = (mesh.processor_id() + pid) % mesh.n_processors();
624  const processor_id_type procdown = (mesh.n_processors() +
625  mesh.processor_id() - pid) % mesh.n_processors();
626 
627  mesh.comm().send_receive (procup, requested_ids[procup],
628  procdown, requests_to_fill[procdown]);
629 
630  // we can overwrite these requested ids in-place.
631  for (std::size_t i=0; i<requests_to_fill[procdown].size(); i++)
632  {
633  const dof_id_type requested_elem_index =
634  requests_to_fill[procdown][i];
635 
636  libmesh_assert(_global_index_by_pid_map.count(requested_elem_index));
637 
638  const dof_id_type global_index_by_pid =
639  _global_index_by_pid_map[requested_elem_index];
640 
641  const dof_id_type local_index =
642  global_index_by_pid - first_local_elem;
643 
644  libmesh_assert_less (local_index, _pmetis->part.size());
645  libmesh_assert_less (local_index, n_active_local_elem);
646 
647  const unsigned int elem_procid =
648  static_cast<unsigned int>(_pmetis->part[local_index]);
649 
650  libmesh_assert_less (elem_procid, static_cast<unsigned int>(_pmetis->nparts));
651 
652  requests_to_fill[procdown][i] = elem_procid;
653  }
654 
655  // Trade back
656  mesh.comm().send_receive (procdown, requests_to_fill[procdown],
657  procup, requested_ids[procup]);
658  }
659 
660  // and finally assign the partitioning.
661  // note we are iterating in exactly the same order
662  // used to build up the request, so we can expect the
663  // required entries to be in the proper sequence.
664  std::vector<unsigned int> counters(mesh.n_processors(), 0);
665  for (auto & elem : mesh.active_element_ptr_range())
666  {
667  const processor_id_type current_pid = elem->processor_id();
668 
669  libmesh_assert_less (counters[current_pid], requested_ids[current_pid].size());
670 
671  const processor_id_type elem_procid =
672  requested_ids[current_pid][counters[current_pid]++];
673 
674  libmesh_assert_less (elem_procid, static_cast<unsigned int>(_pmetis->nparts));
675  elem->processor_id() = elem_procid;
676  }
677 }
678 
679 #endif // #ifdef LIBMESH_HAVE_PARMETIS
680 
681 } // namespace libMesh
OStreamProxy err
difference_type count(const key_type &key) const
Definition: vectormap.h:210
bool active() const
Definition: elem.h:2257
libMesh::BoundingBox create_bounding_box(const MeshBase &mesh)
The same functionality as the deprecated MeshTools::bounding_box().
Definition: mesh_tools.C:329
The definition of the const_element_iterator struct.
Definition: mesh_base.h:1494
virtual void allgather()
Gathers all elements and nodes of the mesh onto every processor.
Definition: mesh_base.h:169
processor_id_type n_processors() const
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
uint8_t processor_id_type
Definition: id_types.h:99
const class libmesh_nullptr_t libmesh_nullptr
The ParmetisHelper class allows us to use a &#39;pimpl&#39; strategy in the ParmetisPartitioner class...
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
void find_global_indices(const Parallel::Communicator &communicator, const libMesh::BoundingBox &, const ForwardIterator &, const ForwardIterator &, std::vector< dof_id_type > &) const
This method determines a globally unique, partition-agnostic index for each object in the input range...
The MetisPartitioner uses the Metis graph partitioner to partition the elements.
The libMesh namespace provides an interface to certain functionality in the library.
void initialize(EquationSystems &es, const std::string &system_name)
Definition: main.C:52
Real distance(const Point &p)
This is the MeshBase class.
Definition: mesh_base.h:68
libmesh_assert(j)
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1967
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
void send_receive(const unsigned int dest_processor_id, const T1 &send, const unsigned int source_processor_id, T2 &recv, const MessageTag &send_tag=no_tag, const MessageTag &recv_tag=any_tag) const
Send data send to one processor while simultaneously receiving other data recv from a (potentially di...
This is the MeshCommunication class.
virtual element_iterator active_pid_elements_begin(processor_id_type proc_id)=0
void insert(const value_type &x)
Inserts x into the vectormap.
Definition: vectormap.h:116
SimpleRange< I > as_range(const std::pair< I, I > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
virtual element_iterator active_elements_begin()=0
Active, local, and negation forms of the element iterators described above.
dof_id_type n_active_local_elem() const
Definition: mesh_base.h:389
Defines a Cartesian bounding box by the two corner extremum.
Definition: bounding_box.h:40
virtual element_iterator active_elements_end()=0
const unsigned int MIN_ELEM_PER_PROC
virtual element_iterator active_pid_elements_end(processor_id_type proc_id)=0
Temporarily serialize a DistributedMesh for output; a distributed mesh is allgathered by the MeshSeri...
virtual void partition(MeshBase &mesh, const unsigned int n)
Partitions the MeshBase into n parts by setting processor_id() on Nodes and Elems.
Definition: partitioner.C:49
const Parallel::Communicator & comm() const
dof_id_type id() const
Definition: dof_object.h:632
long double min(long double a, double b)
processor_id_type processor_id() const
uint8_t dof_id_type
Definition: id_types.h:64
void allgather(const T &send, std::vector< T > &recv) const
Take a vector of length this->size(), and fill in recv[processor_id] = the value of send on that proc...