libMesh
mesh_communication_global_indices.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/libmesh_common.h"
23 #include "libmesh/libmesh_logging.h"
24 #include "libmesh/mesh_base.h"
25 #include "libmesh/mesh_tools.h"
26 #include "libmesh/mesh_communication.h"
27 #include "libmesh/parallel.h"
28 #include "libmesh/parallel_hilbert.h"
29 #include "libmesh/parallel_sort.h"
30 #include "libmesh/elem.h"
31 #include "libmesh/elem_range.h"
32 #include "libmesh/node_range.h"
33 #ifdef LIBMESH_HAVE_LIBHILBERT
34 # include "hilbert.h"
35 #endif
36 
37 #ifdef LIBMESH_HAVE_LIBHILBERT
38 namespace { // anonymous namespace for helper functions
39 
40 using namespace libMesh;
41 
42 // Utility function to map (x,y,z) in [bbox.min, bbox.max]^3 into
43 // [0,max_inttype]^3 for computing Hilbert keys
44 void get_hilbert_coords (const Point & p,
45  const libMesh::BoundingBox & bbox,
46  CFixBitVec icoords[3])
47 {
48  static const Hilbert::inttype max_inttype = static_cast<Hilbert::inttype>(-1);
49 
50  const long double // put (x,y,z) in [0,1]^3 (don't divide by 0)
51  x = ((bbox.first(0) == bbox.second(0)) ? 0. :
52  (p(0)-bbox.first(0))/(bbox.second(0)-bbox.first(0))),
53 
54 #if LIBMESH_DIM > 1
55  y = ((bbox.first(1) == bbox.second(1)) ? 0. :
56  (p(1)-bbox.first(1))/(bbox.second(1)-bbox.first(1))),
57 #else
58  y = 0.,
59 #endif
60 
61 #if LIBMESH_DIM > 2
62  z = ((bbox.first(2) == bbox.second(2)) ? 0. :
63  (p(2)-bbox.first(2))/(bbox.second(2)-bbox.first(2)));
64 #else
65  z = 0.;
66 #endif
67 
68  // (icoords) in [0,max_inttype]^3
69  icoords[0] = static_cast<Hilbert::inttype>(x*max_inttype);
70  icoords[1] = static_cast<Hilbert::inttype>(y*max_inttype);
71  icoords[2] = static_cast<Hilbert::inttype>(z*max_inttype);
72 }
73 
74 
75 
77 get_hilbert_index (const Elem * e,
78  const libMesh::BoundingBox & bbox)
79 {
80  static const unsigned int sizeof_inttype = sizeof(Hilbert::inttype);
81 
82  Hilbert::HilbertIndices index;
83  CFixBitVec icoords[3];
84  Hilbert::BitVecType bv;
85  get_hilbert_coords (e->centroid(), bbox, icoords);
86  Hilbert::coordsToIndex (icoords, 8*sizeof_inttype, 3, bv);
87  index = bv;
88 
89 #ifdef LIBMESH_ENABLE_UNIQUE_ID
90  return std::make_pair(index, e->unique_id());
91 #else
92  return index;
93 #endif
94 }
95 
96 
97 
98 // Compute the hilbert index
100 get_hilbert_index (const Node * n,
101  const libMesh::BoundingBox & bbox)
102 {
103  static const unsigned int sizeof_inttype = sizeof(Hilbert::inttype);
104 
105  Hilbert::HilbertIndices index;
106  CFixBitVec icoords[3];
107  Hilbert::BitVecType bv;
108  get_hilbert_coords (*n, bbox, icoords);
109  Hilbert::coordsToIndex (icoords, 8*sizeof_inttype, 3, bv);
110  index = bv;
111 
112 #ifdef LIBMESH_ENABLE_UNIQUE_ID
113  return std::make_pair(index, n->unique_id());
114 #else
115  return index;
116 #endif
117 }
118 
119 // Helper class for threaded Hilbert key computation
120 class ComputeHilbertKeys
121 {
122 public:
123  ComputeHilbertKeys (const libMesh::BoundingBox & bbox,
124  std::vector<Parallel::DofObjectKey> & keys) :
125  _bbox(bbox),
126  _keys(keys)
127  {}
128 
129  // computes the hilbert index for a node
130  void operator() (const ConstNodeRange & range) const
131  {
132  std::size_t pos = range.first_idx();
133  for (ConstNodeRange::const_iterator it = range.begin(); it!=range.end(); ++it)
134  {
135  const Node * node = (*it);
136  libmesh_assert(node);
137  libmesh_assert_less (pos, _keys.size());
138  _keys[pos++] = get_hilbert_index (node, _bbox);
139  }
140  }
141 
142  // computes the hilbert index for an element
143  void operator() (const ConstElemRange & range) const
144  {
145  std::size_t pos = range.first_idx();
146  for (ConstElemRange::const_iterator it = range.begin(); it!=range.end(); ++it)
147  {
148  const Elem * elem = (*it);
149  libmesh_assert(elem);
150  libmesh_assert_less (pos, _keys.size());
151  _keys[pos++] = get_hilbert_index (elem, _bbox);
152  }
153  }
154 
155 private:
156  const libMesh::BoundingBox & _bbox;
157  std::vector<Parallel::DofObjectKey> & _keys;
158 };
159 }
160 #endif
161 
162 
163 namespace libMesh
164 {
165 
166 // ------------------------------------------------------------
167 // MeshCommunication class members
168 #if defined(LIBMESH_HAVE_LIBHILBERT) && defined(LIBMESH_HAVE_MPI)
170 {
171  LOG_SCOPE ("assign_global_indices()", "MeshCommunication");
172 
173  // This method determines partition-agnostic global indices
174  // for nodes and elements.
175 
176  // Algorithm:
177  // (1) compute the Hilbert key for each local node/element
178  // (2) perform a parallel sort of the Hilbert key
179  // (3) get the min/max value on each processor
180  // (4) determine the position in the global ranking for
181  // each local object
182 
183  const Parallel::Communicator & communicator (mesh.comm());
184 
185  // Global bounding box. We choose the nodal bounding box for
186  // backwards compatibility; the element bounding box may be looser
187  // on curved elements.
188  BoundingBox bbox =
190 
191  //-------------------------------------------------------------
192  // (1) compute Hilbert keys
193  std::vector<Parallel::DofObjectKey>
194  node_keys, elem_keys;
195 
196  {
197  // Nodes first
198  {
200  mesh.local_nodes_end());
201  node_keys.resize (nr.size());
202  Threads::parallel_for (nr, ComputeHilbertKeys (bbox, node_keys));
203 
204  // // It's O(N^2) to check that these keys don't duplicate before the
205  // // sort...
206  // MeshBase::const_node_iterator nodei = mesh.local_nodes_begin();
207  // for (std::size_t i = 0; i != node_keys.size(); ++i, ++nodei)
208  // {
209  // MeshBase::const_node_iterator nodej = mesh.local_nodes_begin();
210  // for (std::size_t j = 0; j != i; ++j, ++nodej)
211  // {
212  // if (node_keys[i] == node_keys[j])
213  // {
214  // CFixBitVec icoords[3], jcoords[3];
215  // get_hilbert_coords(**nodej, bbox, jcoords);
216  // libMesh::err <<
217  // "node " << (*nodej)->id() << ", " <<
218  // *(Point *)(*nodej) << " has HilbertIndices " <<
219  // node_keys[j] << std::endl;
220  // get_hilbert_coords(**nodei, bbox, icoords);
221  // libMesh::err <<
222  // "node " << (*nodei)->id() << ", " <<
223  // *(Point *)(*nodei) << " has HilbertIndices " <<
224  // node_keys[i] << std::endl;
225  // libmesh_error_msg("Error: nodes with duplicate Hilbert keys!");
226  // }
227  // }
228  // }
229  }
230 
231  // Elements next
232  {
234  mesh.local_elements_end());
235  elem_keys.resize (er.size());
236  Threads::parallel_for (er, ComputeHilbertKeys (bbox, elem_keys));
237 
238  // // For elements, the keys can be (and in the case of TRI, are
239  // // expected to be) duplicates, but only if the elements are at
240  // // different levels
241  // MeshBase::const_element_iterator elemi = mesh.local_elements_begin();
242  // for (std::size_t i = 0; i != elem_keys.size(); ++i, ++elemi)
243  // {
244  // MeshBase::const_element_iterator elemj = mesh.local_elements_begin();
245  // for (std::size_t j = 0; j != i; ++j, ++elemj)
246  // {
247  // if ((elem_keys[i] == elem_keys[j]) &&
248  // ((*elemi)->level() == (*elemj)->level()))
249  // {
250  // libMesh::err <<
251  // "level " << (*elemj)->level() << " elem\n" <<
252  // (**elemj) << " centroid " <<
253  // (*elemj)->centroid() << " has HilbertIndices " <<
254  // elem_keys[j] << " or " <<
255  // get_hilbert_index((*elemj), bbox) <<
256  // std::endl;
257  // libMesh::err <<
258  // "level " << (*elemi)->level() << " elem\n" <<
259  // (**elemi) << " centroid " <<
260  // (*elemi)->centroid() << " has HilbertIndices " <<
261  // elem_keys[i] << " or " <<
262  // get_hilbert_index((*elemi), bbox) <<
263  // std::endl;
264  // libmesh_error_msg("Error: level " << (*elemi)->level() << " elements with duplicate Hilbert keys!");
265  // }
266  // }
267  // }
268  }
269  } // done computing Hilbert keys
270 
271 
272 
273  //-------------------------------------------------------------
274  // (2) parallel sort the Hilbert keys
276  node_keys);
277  node_sorter.sort(); /* done with node_keys */ //node_keys.clear();
278 
279  const std::vector<Parallel::DofObjectKey> & my_node_bin =
280  node_sorter.bin();
281 
283  elem_keys);
284  elem_sorter.sort(); /* done with elem_keys */ //elem_keys.clear();
285 
286  const std::vector<Parallel::DofObjectKey> & my_elem_bin =
287  elem_sorter.bin();
288 
289 
290 
291  //-------------------------------------------------------------
292  // (3) get the max value on each processor
293  std::vector<Parallel::DofObjectKey>
294  node_upper_bounds(communicator.size()),
295  elem_upper_bounds(communicator.size());
296 
297  { // limit scope of temporaries
298  std::vector<Parallel::DofObjectKey> recvbuf(2*communicator.size());
299  std::vector<unsigned short int> /* do not use a vector of bools here since it is not always so! */
300  empty_nodes (communicator.size()),
301  empty_elem (communicator.size());
302  std::vector<Parallel::DofObjectKey> my_max(2);
303 
304  communicator.allgather (static_cast<unsigned short int>(my_node_bin.empty()), empty_nodes);
305  communicator.allgather (static_cast<unsigned short int>(my_elem_bin.empty()), empty_elem);
306 
307  if (!my_node_bin.empty()) my_max[0] = my_node_bin.back();
308  if (!my_elem_bin.empty()) my_max[1] = my_elem_bin.back();
309 
310  communicator.allgather (my_max, /* identical_buffer_sizes = */ true);
311 
312  // Be careful here. The *_upper_bounds will be used to find the processor
313  // a given object belongs to. So, if a processor contains no objects (possible!)
314  // then copy the bound from the lower processor id.
315  for (processor_id_type p=0; p<communicator.size(); p++)
316  {
317  node_upper_bounds[p] = my_max[2*p+0];
318  elem_upper_bounds[p] = my_max[2*p+1];
319 
320  if (p > 0) // default hilbert index value is the OK upper bound for processor 0.
321  {
322  if (empty_nodes[p]) node_upper_bounds[p] = node_upper_bounds[p-1];
323  if (empty_elem[p]) elem_upper_bounds[p] = elem_upper_bounds[p-1];
324  }
325  }
326  }
327 
328 
329 
330  //-------------------------------------------------------------
331  // (4) determine the position in the global ranking for
332  // each local object
333  {
334  //----------------------------------------------
335  // Nodes first -- all nodes, not just local ones
336  {
337  // Request sets to send to each processor
338  std::vector<std::vector<Parallel::DofObjectKey>>
339  requested_ids (communicator.size());
340  // Results to gather from each processor
341  std::vector<std::vector<dof_id_type>>
342  filled_request (communicator.size());
343 
344  // build up list of requests
345  for (const auto & node : mesh.node_ptr_range())
346  {
347  libmesh_assert(node);
348  const Parallel::DofObjectKey hi =
349  get_hilbert_index (node, bbox);
350  const processor_id_type pid =
351  cast_int<processor_id_type>
352  (std::distance (node_upper_bounds.begin(),
353  std::lower_bound(node_upper_bounds.begin(),
354  node_upper_bounds.end(),
355  hi)));
356 
357  libmesh_assert_less (pid, communicator.size());
358 
359  requested_ids[pid].push_back(hi);
360  }
361 
362  // The number of objects in my_node_bin on each processor
363  std::vector<dof_id_type> node_bin_sizes(communicator.size());
364  communicator.allgather (static_cast<dof_id_type>(my_node_bin.size()), node_bin_sizes);
365 
366  // The offset of my first global index
367  dof_id_type my_offset = 0;
368  for (processor_id_type pid=0; pid<communicator.rank(); pid++)
369  my_offset += node_bin_sizes[pid];
370 
371  // start with pid=0, so that we will trade with ourself
372  for (processor_id_type pid=0; pid<communicator.size(); pid++)
373  {
374  // Trade my requests with processor procup and procdown
375  const processor_id_type procup = cast_int<processor_id_type>
376  ((communicator.rank() + pid) % communicator.size());
377  const processor_id_type procdown = cast_int<processor_id_type>
378  ((communicator.size() + communicator.rank() - pid) %
379  communicator.size());
380 
381  std::vector<Parallel::DofObjectKey> request_to_fill;
382  communicator.send_receive(procup, requested_ids[procup],
383  procdown, request_to_fill);
384 
385  // Fill the requests
386  std::vector<dof_id_type> global_ids; global_ids.reserve(request_to_fill.size());
387  for (std::size_t idx=0; idx<request_to_fill.size(); idx++)
388  {
389  const Parallel::DofObjectKey & hi = request_to_fill[idx];
390  libmesh_assert_less_equal (hi, node_upper_bounds[communicator.rank()]);
391 
392  // find the requested index in my node bin
393  std::vector<Parallel::DofObjectKey>::const_iterator pos =
394  std::lower_bound (my_node_bin.begin(), my_node_bin.end(), hi);
395  libmesh_assert (pos != my_node_bin.end());
396  libmesh_assert_equal_to (*pos, hi);
397 
398  // Finally, assign the global index based off the position of the index
399  // in my array, properly offset.
400  global_ids.push_back(cast_int<dof_id_type>(std::distance(my_node_bin.begin(), pos) + my_offset));
401  }
402 
403  // and trade back
404  communicator.send_receive (procdown, global_ids,
405  procup, filled_request[procup]);
406  }
407 
408  // We now have all the filled requests, so we can loop through our
409  // nodes once and assign the global index to each one.
410  {
411  std::vector<std::vector<dof_id_type>::const_iterator>
412  next_obj_on_proc; next_obj_on_proc.reserve(communicator.size());
413  for (processor_id_type pid=0; pid<communicator.size(); pid++)
414  next_obj_on_proc.push_back(filled_request[pid].begin());
415 
416  for (auto & node : mesh.node_ptr_range())
417  {
418  libmesh_assert(node);
419  const Parallel::DofObjectKey hi =
420  get_hilbert_index (node, bbox);
421  const processor_id_type pid =
422  cast_int<processor_id_type>
423  (std::distance (node_upper_bounds.begin(),
424  std::lower_bound(node_upper_bounds.begin(),
425  node_upper_bounds.end(),
426  hi)));
427 
428  libmesh_assert_less (pid, communicator.size());
429  libmesh_assert (next_obj_on_proc[pid] != filled_request[pid].end());
430 
431  const dof_id_type global_index = *next_obj_on_proc[pid];
432  libmesh_assert_less (global_index, mesh.n_nodes());
433  node->set_id() = global_index;
434 
435  ++next_obj_on_proc[pid];
436  }
437  }
438  }
439 
440  //---------------------------------------------------
441  // elements next -- all elements, not just local ones
442  {
443  // Request sets to send to each processor
444  std::vector<std::vector<Parallel::DofObjectKey>>
445  requested_ids (communicator.size());
446  // Results to gather from each processor
447  std::vector<std::vector<dof_id_type>>
448  filled_request (communicator.size());
449 
450  for (const auto & elem : mesh.element_ptr_range())
451  {
452  libmesh_assert(elem);
453  const Parallel::DofObjectKey hi =
454  get_hilbert_index (elem, bbox);
455  const processor_id_type pid =
456  cast_int<processor_id_type>
457  (std::distance (elem_upper_bounds.begin(),
458  std::lower_bound(elem_upper_bounds.begin(),
459  elem_upper_bounds.end(),
460  hi)));
461 
462  libmesh_assert_less (pid, communicator.size());
463 
464  requested_ids[pid].push_back(hi);
465  }
466 
467  // The number of objects in my_elem_bin on each processor
468  std::vector<dof_id_type> elem_bin_sizes(communicator.size());
469  communicator.allgather (static_cast<dof_id_type>(my_elem_bin.size()), elem_bin_sizes);
470 
471  // The offset of my first global index
472  dof_id_type my_offset = 0;
473  for (processor_id_type pid=0; pid<communicator.rank(); pid++)
474  my_offset += elem_bin_sizes[pid];
475 
476  // start with pid=0, so that we will trade with ourself
477  for (processor_id_type pid=0; pid<communicator.size(); pid++)
478  {
479  // Trade my requests with processor procup and procdown
480  const processor_id_type procup = cast_int<processor_id_type>
481  ((communicator.rank() + pid) % communicator.size());
482  const processor_id_type procdown = cast_int<processor_id_type>
483  ((communicator.size() + communicator.rank() - pid) %
484  communicator.size());
485 
486  std::vector<Parallel::DofObjectKey> request_to_fill;
487  communicator.send_receive(procup, requested_ids[procup],
488  procdown, request_to_fill);
489 
490  // Fill the requests
491  std::vector<dof_id_type> global_ids; global_ids.reserve(request_to_fill.size());
492  for (std::size_t idx=0; idx<request_to_fill.size(); idx++)
493  {
494  const Parallel::DofObjectKey & hi = request_to_fill[idx];
495  libmesh_assert_less_equal (hi, elem_upper_bounds[communicator.rank()]);
496 
497  // find the requested index in my elem bin
498  std::vector<Parallel::DofObjectKey>::const_iterator pos =
499  std::lower_bound (my_elem_bin.begin(), my_elem_bin.end(), hi);
500  libmesh_assert (pos != my_elem_bin.end());
501  libmesh_assert_equal_to (*pos, hi);
502 
503  // Finally, assign the global index based off the position of the index
504  // in my array, properly offset.
505  global_ids.push_back (cast_int<dof_id_type>(std::distance(my_elem_bin.begin(), pos) + my_offset));
506  }
507 
508  // and trade back
509  communicator.send_receive (procdown, global_ids,
510  procup, filled_request[procup]);
511  }
512 
513  // We now have all the filled requests, so we can loop through our
514  // elements once and assign the global index to each one.
515  {
516  std::vector<std::vector<dof_id_type>::const_iterator>
517  next_obj_on_proc; next_obj_on_proc.reserve(communicator.size());
518  for (processor_id_type pid=0; pid<communicator.size(); pid++)
519  next_obj_on_proc.push_back(filled_request[pid].begin());
520 
521  for (auto & elem : mesh.element_ptr_range())
522  {
523  libmesh_assert(elem);
524  const Parallel::DofObjectKey hi =
525  get_hilbert_index (elem, bbox);
526  const processor_id_type pid =
527  cast_int<processor_id_type>
528  (std::distance (elem_upper_bounds.begin(),
529  std::lower_bound(elem_upper_bounds.begin(),
530  elem_upper_bounds.end(),
531  hi)));
532 
533  libmesh_assert_less (pid, communicator.size());
534  libmesh_assert (next_obj_on_proc[pid] != filled_request[pid].end());
535 
536  const dof_id_type global_index = *next_obj_on_proc[pid];
537  libmesh_assert_less (global_index, mesh.n_elem());
538  elem->set_id() = global_index;
539 
540  ++next_obj_on_proc[pid];
541  }
542  }
543  }
544  }
545 }
546 #else // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
548 {
549 }
550 #endif // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
551 
552 
553 #if defined(LIBMESH_HAVE_LIBHILBERT) && defined(LIBMESH_HAVE_MPI)
555 {
556  LOG_SCOPE ("check_for_duplicate_global_indices()", "MeshCommunication");
557 
558  // Global bounding box. We choose the nodal bounding box for
559  // backwards compatibility; the element bounding box may be looser
560  // on curved elements.
561  BoundingBox bbox =
563 
564 
565  std::vector<Parallel::DofObjectKey>
566  node_keys, elem_keys;
567 
568  {
569  // Nodes first
570  {
572  mesh.local_nodes_end());
573  node_keys.resize (nr.size());
574  Threads::parallel_for (nr, ComputeHilbertKeys (bbox, node_keys));
575 
576  // It's O(N^2) to check that these keys don't duplicate before the
577  // sort...
579  for (std::size_t i = 0; i != node_keys.size(); ++i, ++nodei)
580  {
582  for (std::size_t j = 0; j != i; ++j, ++nodej)
583  {
584  if (node_keys[i] == node_keys[j])
585  {
586  CFixBitVec icoords[3], jcoords[3];
587  get_hilbert_coords(**nodej, bbox, jcoords);
588  libMesh::err <<
589  "node " << (*nodej)->id() << ", " <<
590  *(Point *)(*nodej) << " has HilbertIndices " <<
591  node_keys[j] << std::endl;
592  get_hilbert_coords(**nodei, bbox, icoords);
593  libMesh::err <<
594  "node " << (*nodei)->id() << ", " <<
595  *(Point *)(*nodei) << " has HilbertIndices " <<
596  node_keys[i] << std::endl;
597  libmesh_error_msg("Error: nodes with duplicate Hilbert keys!");
598  }
599  }
600  }
601  }
602 
603  // Elements next
604  {
606  mesh.local_elements_end());
607  elem_keys.resize (er.size());
608  Threads::parallel_for (er, ComputeHilbertKeys (bbox, elem_keys));
609 
610  // For elements, the keys can be (and in the case of TRI, are
611  // expected to be) duplicates, but only if the elements are at
612  // different levels
614  for (std::size_t i = 0; i != elem_keys.size(); ++i, ++elemi)
615  {
617  for (std::size_t j = 0; j != i; ++j, ++elemj)
618  {
619  if ((elem_keys[i] == elem_keys[j]) &&
620  ((*elemi)->level() == (*elemj)->level()))
621  {
622  libMesh::err <<
623  "level " << (*elemj)->level() << " elem\n" <<
624  (**elemj) << " centroid " <<
625  (*elemj)->centroid() << " has HilbertIndices " <<
626  elem_keys[j] << " or " <<
627  get_hilbert_index((*elemj), bbox) <<
628  std::endl;
629  libMesh::err <<
630  "level " << (*elemi)->level() << " elem\n" <<
631  (**elemi) << " centroid " <<
632  (*elemi)->centroid() << " has HilbertIndices " <<
633  elem_keys[i] << " or " <<
634  get_hilbert_index((*elemi), bbox) <<
635  std::endl;
636  libmesh_error_msg("Error: level " << (*elemi)->level() << " elements with duplicate Hilbert keys!");
637  }
638  }
639  }
640  }
641  } // done checking Hilbert keys
642 }
643 #else // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
645 {
646 }
647 #endif // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
648 
649 #if defined(LIBMESH_HAVE_LIBHILBERT) && defined(LIBMESH_HAVE_MPI)
650 template <typename ForwardIterator>
652  const libMesh::BoundingBox & bbox,
653  const ForwardIterator & begin,
654  const ForwardIterator & end,
655  std::vector<dof_id_type> & index_map) const
656 {
657  LOG_SCOPE ("find_global_indices()", "MeshCommunication");
658 
659  // This method determines partition-agnostic global indices
660  // for nodes and elements.
661 
662  // Algorithm:
663  // (1) compute the Hilbert key for each local node/element
664  // (2) perform a parallel sort of the Hilbert key
665  // (3) get the min/max value on each processor
666  // (4) determine the position in the global ranking for
667  // each local object
668  index_map.clear();
669  std::size_t n_objects = std::distance (begin, end);
670  index_map.reserve(n_objects);
671 
672  //-------------------------------------------------------------
673  // (1) compute Hilbert keys
674  // These aren't trivial to compute, and we will need them again.
675  // But the binsort will sort the input vector, trashing the order
676  // that we'd like to rely on. So, two vectors...
677  std::vector<Parallel::DofObjectKey>
678  sorted_hilbert_keys,
679  hilbert_keys;
680  sorted_hilbert_keys.reserve(n_objects);
681  hilbert_keys.reserve(n_objects);
682  {
683  LOG_SCOPE("compute_hilbert_indices()", "MeshCommunication");
684  for (ForwardIterator it=begin; it!=end; ++it)
685  {
686  const Parallel::DofObjectKey hi(get_hilbert_index (*it, bbox));
687  hilbert_keys.push_back(hi);
688 
689  if ((*it)->processor_id() == communicator.rank())
690  sorted_hilbert_keys.push_back(hi);
691 
692  // someone needs to take care of unpartitioned objects!
693  if ((communicator.rank() == 0) &&
694  ((*it)->processor_id() == DofObject::invalid_processor_id))
695  sorted_hilbert_keys.push_back(hi);
696  }
697  }
698 
699  //-------------------------------------------------------------
700  // (2) parallel sort the Hilbert keys
701  START_LOG ("parallel_sort()", "MeshCommunication");
702  Parallel::Sort<Parallel::DofObjectKey> sorter (communicator,
703  sorted_hilbert_keys);
704  sorter.sort();
705  STOP_LOG ("parallel_sort()", "MeshCommunication");
706  const std::vector<Parallel::DofObjectKey> & my_bin = sorter.bin();
707 
708  // The number of objects in my_bin on each processor
709  std::vector<unsigned int> bin_sizes(communicator.size());
710  communicator.allgather (static_cast<unsigned int>(my_bin.size()), bin_sizes);
711 
712  // The offset of my first global index
713  unsigned int my_offset = 0;
714  for (unsigned int pid=0; pid<communicator.rank(); pid++)
715  my_offset += bin_sizes[pid];
716 
717  //-------------------------------------------------------------
718  // (3) get the max value on each processor
719  std::vector<Parallel::DofObjectKey>
720  upper_bounds(1);
721 
722  if (!my_bin.empty())
723  upper_bounds[0] = my_bin.back();
724 
725  communicator.allgather (upper_bounds, /* identical_buffer_sizes = */ true);
726 
727  // Be careful here. The *_upper_bounds will be used to find the processor
728  // a given object belongs to. So, if a processor contains no objects (possible!)
729  // then copy the bound from the lower processor id.
730  for (unsigned int p=1; p<communicator.size(); p++)
731  if (!bin_sizes[p]) upper_bounds[p] = upper_bounds[p-1];
732 
733 
734  //-------------------------------------------------------------
735  // (4) determine the position in the global ranking for
736  // each local object
737  {
738  //----------------------------------------------
739  // all objects, not just local ones
740 
741  // Request sets to send to each processor
742  std::vector<std::vector<Parallel::DofObjectKey>>
743  requested_ids (communicator.size());
744  // Results to gather from each processor
745  std::vector<std::vector<dof_id_type>>
746  filled_request (communicator.size());
747 
748  // build up list of requests
749  std::vector<Parallel::DofObjectKey>::const_iterator hi =
750  hilbert_keys.begin();
751 
752  for (ForwardIterator it = begin; it != end; ++it)
753  {
754  libmesh_assert (hi != hilbert_keys.end());
755 
756  std::vector<Parallel::DofObjectKey>::iterator lb =
757  std::lower_bound(upper_bounds.begin(), upper_bounds.end(),
758  *hi);
759 
760  const processor_id_type pid =
761  cast_int<processor_id_type>
762  (std::distance (upper_bounds.begin(), lb));
763 
764  libmesh_assert_less (pid, communicator.size());
765 
766  requested_ids[pid].push_back(*hi);
767 
768  ++hi;
769  // go ahead and put pid in index_map, that way we
770  // don't have to repeat the std::lower_bound()
771  index_map.push_back(pid);
772  }
773 
774  // start with pid=0, so that we will trade with ourself
775  std::vector<Parallel::DofObjectKey> request_to_fill;
776  std::vector<dof_id_type> global_ids;
777  for (processor_id_type pid=0; pid<communicator.size(); pid++)
778  {
779  // Trade my requests with processor procup and procdown
780  const processor_id_type procup = cast_int<processor_id_type>
781  ((communicator.rank() + pid) % communicator.size());
782  const processor_id_type procdown = cast_int<processor_id_type>
783  ((communicator.size() + communicator.rank() - pid) %
784  communicator.size());
785 
786  communicator.send_receive(procup, requested_ids[procup],
787  procdown, request_to_fill);
788 
789  // Fill the requests
790  global_ids.clear(); global_ids.reserve(request_to_fill.size());
791  for (std::size_t idx=0; idx<request_to_fill.size(); idx++)
792  {
793  const Parallel::DofObjectKey & hilbert_indices = request_to_fill[idx];
794  libmesh_assert_less_equal (hilbert_indices, upper_bounds[communicator.rank()]);
795 
796  // find the requested index in my node bin
797  std::vector<Parallel::DofObjectKey>::const_iterator pos =
798  std::lower_bound (my_bin.begin(), my_bin.end(), hilbert_indices);
799  libmesh_assert (pos != my_bin.end());
800 #ifdef DEBUG
801  // If we could not find the requested Hilbert index in
802  // my_bin, something went terribly wrong, possibly the
803  // Mesh was displaced differently on different processors,
804  // and therefore the Hilbert indices don't agree.
805  if (*pos != hilbert_indices)
806  {
807  // The input will be hilbert_indices. We convert it
808  // to BitVecType using the operator= provided by the
809  // BitVecType class. BitVecType is a CBigBitVec!
810  Hilbert::BitVecType input;
811 #ifdef LIBMESH_ENABLE_UNIQUE_ID
812  input = hilbert_indices.first;
813 #else
814  input = hilbert_indices;
815 #endif
816 
817  // Get output in a vector of CBigBitVec
818  std::vector<CBigBitVec> output(3);
819 
820  // Call the indexToCoords function
821  Hilbert::indexToCoords(&output[0], 8*sizeof(Hilbert::inttype), 3, input);
822 
823  // The entries in the output racks are integers in the
824  // range [0, Hilbert::inttype::max] which can be
825  // converted to floating point values in [0,1] and
826  // finally to actual values using the bounding box.
827  Real max_int_as_real = static_cast<Real>(std::numeric_limits<Hilbert::inttype>::max());
828 
829  // Get the points in [0,1]^3. The zeroth rack of each entry in
830  // 'output' maps to the normalized x, y, and z locations,
831  // respectively.
832  Point p_hat(static_cast<Real>(output[0].racks()[0]) / max_int_as_real,
833  static_cast<Real>(output[1].racks()[0]) / max_int_as_real,
834  static_cast<Real>(output[2].racks()[0]) / max_int_as_real);
835 
836  // Convert the points from [0,1]^3 to their actual (x,y,z) locations
837  Real
838  xmin = bbox.first(0),
839  xmax = bbox.second(0),
840  ymin = bbox.first(1),
841  ymax = bbox.second(1),
842  zmin = bbox.first(2),
843  zmax = bbox.second(2);
844 
845  // Convert the points from [0,1]^3 to their actual (x,y,z) locations
846  Point p(xmin + (xmax-xmin)*p_hat(0),
847  ymin + (ymax-ymin)*p_hat(1),
848  zmin + (zmax-zmin)*p_hat(2));
849 
850  libmesh_error_msg("Could not find hilbert indices: "
851  << hilbert_indices
852  << " corresponding to point " << p);
853  }
854 #endif
855 
856  // Finally, assign the global index based off the position of the index
857  // in my array, properly offset.
858  global_ids.push_back (cast_int<dof_id_type>(std::distance(my_bin.begin(), pos) + my_offset));
859  }
860 
861  // and trade back
862  communicator.send_receive (procdown, global_ids,
863  procup, filled_request[procup]);
864  }
865 
866  // We now have all the filled requests, so we can loop through our
867  // nodes once and assign the global index to each one.
868  {
869  std::vector<std::vector<dof_id_type>::const_iterator>
870  next_obj_on_proc; next_obj_on_proc.reserve(communicator.size());
871  for (unsigned int pid=0; pid<communicator.size(); pid++)
872  next_obj_on_proc.push_back(filled_request[pid].begin());
873 
874  unsigned int cnt=0;
875  for (ForwardIterator it = begin; it != end; ++it, cnt++)
876  {
877  const processor_id_type pid = cast_int<processor_id_type>
878  (index_map[cnt]);
879 
880  libmesh_assert_less (pid, communicator.size());
881  libmesh_assert (next_obj_on_proc[pid] != filled_request[pid].end());
882 
883  const dof_id_type global_index = *next_obj_on_proc[pid];
884  index_map[cnt] = global_index;
885 
886  ++next_obj_on_proc[pid];
887  }
888  }
889  }
890 
891  libmesh_assert_equal_to(index_map.size(), n_objects);
892 }
893 #else // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
894 template <typename ForwardIterator>
896  const libMesh::BoundingBox &,
897  const ForwardIterator & begin,
898  const ForwardIterator & end,
899  std::vector<dof_id_type> & index_map) const
900 {
901  index_map.clear();
902  index_map.reserve(std::distance (begin, end));
903  unsigned int index = 0;
904  for (ForwardIterator it=begin; it!=end; ++it)
905  index_map.push_back(index++);
906 }
907 #endif // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
908 
909 
910 
911 //------------------------------------------------------------------
912 template void MeshCommunication::find_global_indices<MeshBase::const_node_iterator> (const Parallel::Communicator &,
913  const libMesh::BoundingBox &,
915  const MeshBase::const_node_iterator &,
916  std::vector<dof_id_type> &) const;
917 
918 template void MeshCommunication::find_global_indices<MeshBase::const_element_iterator> (const Parallel::Communicator &,
919  const libMesh::BoundingBox &,
921  const MeshBase::const_element_iterator &,
922  std::vector<dof_id_type> &) const;
923 template void MeshCommunication::find_global_indices<MeshBase::node_iterator> (const Parallel::Communicator &,
924  const libMesh::BoundingBox &,
925  const MeshBase::node_iterator &,
926  const MeshBase::node_iterator &,
927  std::vector<dof_id_type> &) const;
928 
929 template void MeshCommunication::find_global_indices<MeshBase::element_iterator> (const Parallel::Communicator &,
930  const libMesh::BoundingBox &,
932  const MeshBase::element_iterator &,
933  std::vector<dof_id_type> &) const;
934 
935 } // namespace libMesh
The definition of the element_iterator struct.
Definition: mesh_base.h:1476
OStreamProxy err
void assign_global_indices(MeshBase &) const
This method assigns globally unique, partition-agnostic indices to the nodes and elements in the mesh...
A Node is like a Point, but with more information.
Definition: node.h:52
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
const std::vector< KeyType > & bin()
Return a constant reference to _my_bin.
unsigned int size() const
Definition: parallel.h:726
void parallel_for(const Range &range, const Body &body)
Execute the provided function object in parallel on the specified range.
Definition: threads_none.h:73
The definition of the const_element_iterator struct.
Definition: mesh_base.h:1494
std::pair< Hilbert::HilbertIndices, unique_id_type > DofObjectKey
MPI_Comm communicator
Communicator object for talking with subsets of processors.
Definition: parallel.h:181
virtual element_iterator local_elements_begin()=0
void sort()
This is the only method which needs to be called by the user.
Definition: parallel_sort.C:58
std::size_t first_idx() const
Index in the stored vector of the first object.
Definition: stored_range.h:271
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
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
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 StoredRange class defines a contiguous, divisible set of objects.
Definition: stored_range.h:52
The libMesh namespace provides an interface to certain functionality in the library.
const_iterator begin() const
Beginning of the range.
Definition: stored_range.h:261
long double max(long double a, double b)
Real distance(const Point &p)
This is the MeshBase class.
Definition: mesh_base.h:68
libmesh_assert(j)
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...
PetscErrorCode Vec x
virtual SimpleRange< element_iterator > element_ptr_range()=0
virtual node_iterator local_nodes_end()=0
static const processor_id_type invalid_processor_id
An invalid processor_id to distinguish DoFs that have not been assigned to a processor.
Definition: dof_object.h:335
virtual SimpleRange< node_iterator > node_ptr_range()=0
virtual element_iterator local_elements_end()=0
void check_for_duplicate_global_indices(MeshBase &) const
Throw an error if we have any index clashes in the numbering used by assign_global_indices.
const_iterator end() const
End of the range.
Definition: stored_range.h:266
std::vector< object_type >::const_iterator const_iterator
Allows an StoredRange to behave like an STL container.
Definition: stored_range.h:58
Defines a Cartesian bounding box by the two corner extremum.
Definition: bounding_box.h:40
The definition of the node_iterator struct.
Definition: mesh_base.h:1528
virtual node_iterator local_nodes_begin()=0
Iterate over local nodes (nodes whose processor_id() matches the current processor).
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
The parallel sorting method is templated on the type of data which is to be sorted.
Definition: parallel_sort.h:53
const Parallel::Communicator & comm() const
The definition of the const_node_iterator struct.
Definition: mesh_base.h:1548
virtual Point centroid() const
Definition: elem.C:446
unsigned int rank() const
Definition: parallel.h:724
libMesh::BoundingBox create_nodal_bounding_box(const MeshBase &mesh)
Definition: mesh_tools.C:354
virtual dof_id_type n_nodes() const =0
virtual dof_id_type n_elem() const =0
unique_id_type unique_id() const
Definition: dof_object.h:649
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
A useful inline function which replaces the macros used previously.
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...