mesh_communication.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 // C++ Includes -----------------------------------
21 #include <numeric>
22 
23 // Local Includes -----------------------------------
24 #include "libmesh/boundary_info.h"
25 #include "libmesh/elem.h"
26 #include "libmesh/libmesh_config.h"
27 #include "libmesh/libmesh_common.h"
29 #include "libmesh/location_maps.h"
30 #include "libmesh/mesh_base.h"
33 #include "libmesh/mesh_tools.h"
34 #include "libmesh/parallel.h"
35 #include "libmesh/parallel_elem.h"
36 #include "libmesh/parallel_mesh.h"
37 #include "libmesh/parallel_node.h"
39 #include "libmesh/utility.h"
40 #include "libmesh/remote_elem.h"
41 
42 
43 
44 //-----------------------------------------------
45 // anonymous namespace for implementation details
46 namespace {
47 
48  using libMesh::Elem;
49 
56  struct CompareElemIdsByLevel
57  {
58  bool operator()(const Elem *a,
59  const Elem *b) const
60  {
61  libmesh_assert (a);
62  libmesh_assert (b);
63  const unsigned int
64  al = a->level(), bl = b->level();
65  const dof_id_type
66  aid = a->id(), bid = b->id();
67 
68  return (al == bl) ? aid < bid : al < bl;
69  }
70  };
71 }
72 
73 
74 namespace libMesh
75 {
76 
77 
78 // ------------------------------------------------------------
79 // MeshCommunication class members
81 {
82  // _neighboring_processors.clear();
83 }
84 
85 
86 
87 #ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
88 // ------------------------------------------------------------
90 {
91  // no MPI == one processor, no redistribution
92  return;
93 }
94 #else
95 // ------------------------------------------------------------
97 {
98  // This method will be called after a new partitioning has been
99  // assigned to the elements. This partitioning was defined in
100  // terms of the active elements, and "trickled down" to the
101  // parents and nodes as to be consistent.
102  //
103  // The point is that the entire concept of local elements is
104  // kinda shaky in this method. Elements which were previously
105  // local may now be assigned to other processors, so we need to
106  // send those off. Similarly, we need to accept elements from
107  // other processors.
108  //
109  // The approach is as follows:
110  // (1) send all elements we have stored to their proper homes
111  // (2) receive elements from all processors, watching for duplicates
112  // (3) deleting all nonlocal elements elements
113  // (4) obtaining required ghost elements from neighboring processors
114  libmesh_parallel_only(mesh.comm());
115  libmesh_assert (!mesh.is_serial());
117  mesh.unpartitioned_elements_end()) == 0);
118 
119  START_LOG("redistribute()","MeshCommunication");
120 
121  // Get a few unique message tags to use in communications; we'll
122  // default to some numbers around pi*1000
124  nodestag = mesh.comm().get_unique_tag(3141),
125  elemstag = mesh.comm().get_unique_tag(3142);
126 
127  // Figure out how many nodes and elements we have which are assigned to each
128  // processor. send_n_nodes_and_elem_per_proc contains the number of nodes/elements
129  // we will be sending to each processor, recv_n_nodes_and_elem_per_proc contains
130  // the number of nodes/elements we will be receiving from each processor.
131  // Format:
132  // send_n_nodes_and_elem_per_proc[2*pid+0] = number of nodes to send to pid
133  // send_n_nodes_and_elem_per_proc[2*pid+1] = number of elements to send to pid
134  std::vector<dof_id_type> send_n_nodes_and_elem_per_proc(2*mesh.n_processors(), 0);
135 
136  std::vector<Parallel::Request>
137  node_send_requests, element_send_requests;
138 
139  for (processor_id_type pid=0; pid<mesh.n_processors(); pid++)
140  if (pid != mesh.processor_id()) // don't send to ourselves!!
141  {
142  // Build up a list of nodes and elements to send to processor pid.
143  // We will certainly send all the elements assigned to this processor,
144  // but we will also ship off any other elements which touch
145  // their nodes.
146  std::set<const Node*> connected_nodes;
147  {
149  const MeshBase::const_element_iterator elem_end = mesh.pid_elements_end(pid);
150 
151  for (; elem_it!=elem_end; ++elem_it)
152  {
153  const Elem* elem = *elem_it;
154 
155  for (unsigned int n=0; n<elem->n_nodes(); n++)
156  connected_nodes.insert (elem->get_node(n));
157  }
158  }
159 
160  std::set<const Elem*, CompareElemIdsByLevel> elements_to_send;
161  {
162  MeshBase::const_element_iterator elem_it = mesh.elements_begin();
163  const MeshBase::const_element_iterator elem_end = mesh.elements_end();
164 
165  for (; elem_it!=elem_end; ++elem_it)
166  {
167  const Elem* elem = *elem_it;
168 
169  for (unsigned int n=0; n<elem->n_nodes(); n++)
170  if (connected_nodes.count(elem->get_node(n)))
171  elements_to_send.insert (elem);
172  }
173  }
174 
175  connected_nodes.clear();
176  {
177  std::set<const Elem*, CompareElemIdsByLevel>::iterator
178  elem_it = elements_to_send.begin(),
179  elem_end = elements_to_send.end();
180 
181  for (; elem_it!=elem_end; ++elem_it)
182  {
183  const Elem* elem = *elem_it;
184 
185  for (unsigned int n=0; n<elem->n_nodes(); n++)
186  connected_nodes.insert(elem->get_node(n));
187  }
188  }
189 
190  // the number of nodes we will ship to pid
191  send_n_nodes_and_elem_per_proc[2*pid+0] = connected_nodes.size();
192 
193  // send any nodes off to the destination processor
194  if (!connected_nodes.empty())
195  {
196  node_send_requests.push_back(Parallel::request());
197 
198  mesh.comm().send_packed_range (pid,
199  &mesh,
200  connected_nodes.begin(),
201  connected_nodes.end(),
202  node_send_requests.back(),
203  nodestag);
204  }
205 
206  // the number of elements we will send to this processor
207  send_n_nodes_and_elem_per_proc[2*pid+1] = elements_to_send.size();
208 
209  if (!elements_to_send.empty())
210  {
211  // send the elements off to the destination processor
212  element_send_requests.push_back(Parallel::request());
213 
214  mesh.comm().send_packed_range (pid,
215  &mesh,
216  elements_to_send.begin(),
217  elements_to_send.end(),
218  element_send_requests.back(),
219  elemstag);
220  }
221  }
222 
223  std::vector<dof_id_type> recv_n_nodes_and_elem_per_proc(send_n_nodes_and_elem_per_proc);
224 
225  mesh.comm().alltoall (recv_n_nodes_and_elem_per_proc);
226 
227  // In general we will only need to communicate with a subset of the other processors.
228  // I can't immediately think of a case where we will send elements but not nodes, but
229  // these are only bools and we have the information anyway...
230  std::vector<bool>
231  send_node_pair(mesh.n_processors(),false), send_elem_pair(mesh.n_processors(),false),
232  recv_node_pair(mesh.n_processors(),false), recv_elem_pair(mesh.n_processors(),false);
233 
234  unsigned int
235  n_send_node_pairs=0, n_send_elem_pairs=0,
236  n_recv_node_pairs=0, n_recv_elem_pairs=0;
237 
238  for (processor_id_type pid=0; pid<mesh.n_processors(); pid++)
239  {
240  if (send_n_nodes_and_elem_per_proc[2*pid+0]) // we have nodes to send
241  {
242  send_node_pair[pid] = true;
243  n_send_node_pairs++;
244  }
245 
246  if (send_n_nodes_and_elem_per_proc[2*pid+1]) // we have elements to send
247  {
248  send_elem_pair[pid] = true;
249  n_send_elem_pairs++;
250  }
251 
252  if (recv_n_nodes_and_elem_per_proc[2*pid+0]) // we have nodes to receive
253  {
254  recv_node_pair[pid] = true;
255  n_recv_node_pairs++;
256  }
257 
258  if (recv_n_nodes_and_elem_per_proc[2*pid+1]) // we have elements to receive
259  {
260  recv_elem_pair[pid] = true;
261  n_recv_elem_pairs++;
262  }
263  }
264  libmesh_assert_equal_to (n_send_node_pairs, node_send_requests.size());
265  libmesh_assert_equal_to (n_send_elem_pairs, element_send_requests.size());
266 
267  // Receive nodes.
268  // We now know how many processors will be sending us nodes.
269  for (unsigned int node_comm_step=0; node_comm_step<n_recv_node_pairs; node_comm_step++)
270  // but we don't necessarily want to impose an ordering, so
271  // just process whatever message is available next.
273  &mesh,
274  mesh_inserter_iterator<Node>(mesh),
275  nodestag);
276 
277  // Receive elements.
278  // Similarly we know how many processors are sending us elements,
279  // but we don't really care in what order we receive them.
280  for (unsigned int elem_comm_step=0; elem_comm_step<n_recv_elem_pairs; elem_comm_step++)
282  &mesh,
283  mesh_inserter_iterator<Elem>(mesh),
284  elemstag);
285 
286  // Wait for all sends to complete
287  Parallel::wait (node_send_requests);
288  Parallel::wait (element_send_requests);
289 
290  // Check on the redistribution consistency
291 #ifdef DEBUG
293 
295 #endif
296 
297  STOP_LOG("redistribute()","MeshCommunication");
298 }
299 #endif // LIBMESH_HAVE_MPI
300 
301 
302 
303 #ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
304 // ------------------------------------------------------------
306 {
307  // no MPI == one processor, no need for this method...
308  return;
309 }
310 #else
311 // ------------------------------------------------------------
313 {
314  // Don't need to do anything if there is
315  // only one processor.
316  if (mesh.n_processors() == 1)
317  return;
318 
319  // This function must be run on all processors at once
320  libmesh_parallel_only(mesh.comm());
321 
322  START_LOG("gather_neighboring_elements()","MeshCommunication");
323 
324  //------------------------------------------------------------------
325  // The purpose of this function is to provide neighbor data structure
326  // consistency for a parallel, distributed mesh. In libMesh we require
327  // that each local element have access to a full set of valid face
328  // neighbors. In some cases this requires us to store "ghost elements" -
329  // elements that belong to other processors but we store to provide
330  // data structure consistency. Also, it is assumed that any element
331  // with a NULL neighbor resides on a physical domain boundary. So,
332  // even our "ghost elements" must have non-NULL neighbors. To handle
333  // this the concept of "RemoteElem" is used - a special construct which
334  // is used to denote that an element has a face neighbor, but we do
335  // not actually store detailed information about that neighbor. This
336  // is required to prevent data structure explosion.
337  //
338  // So when this method is called we should have only local elements.
339  // These local elements will then find neighbors among the local
340  // element set. After this is completed, any element with a NULL
341  // neighbor has either (i) a face on the physical boundary of the mesh,
342  // or (ii) a neighboring element which lives on a remote processor.
343  // To handle case (ii), we communicate the global node indices connected
344  // to all such faces to our neighboring processors. They then send us
345  // all their elements with a NULL neighbor that are connected to any
346  // of the nodes in our list.
347  //------------------------------------------------------------------
348 
349  // Let's begin with finding consistent neighbor data information
350  // for all the elements we currently have. We'll use a clean
351  // slate here - clear any existing information, including RemoteElem's.
352 // mesh.find_neighbors (/* reset_remote_elements = */ true,
353 // /* reset_current_list = */ true);
354 
355  // Get a unique message tag to use in communications; we'll default
356  // to some numbers around pi*10000
358  element_neighbors_tag = mesh.comm().get_unique_tag(31416);
359 
360  // Now any element with a NULL neighbor either
361  // (i) lives on the physical domain boundary, or
362  // (ii) lives on an inter-processor boundary.
363  // We will now gather all the elements from adjacent processors
364  // which are of the same state, which should address all the type (ii)
365  // elements.
366 
367  // A list of all the processors which *may* contain neighboring elements.
368  // (for development simplicity, just make this the identity map)
369  std::vector<processor_id_type> adjacent_processors;
370  for (processor_id_type pid=0; pid<mesh.n_processors(); pid++)
371  if (pid != mesh.processor_id())
372  adjacent_processors.push_back (pid);
373 
374 
375  const processor_id_type n_adjacent_processors =
376  libmesh_cast_int<processor_id_type>(adjacent_processors.size());
377 
378  //-------------------------------------------------------------------------
379  // Let's build a list of all nodes which live on NULL-neighbor sides.
380  // For simplicity, we will use a set to build the list, then transfer
381  // it to a vector for communication.
382  std::vector<dof_id_type> my_interface_node_list;
383  std::vector<const Elem*> my_interface_elements;
384  {
385  std::set<dof_id_type> my_interface_node_set;
386 
387  // since parent nodes are a subset of children nodes, this should be sufficient
390 
391  for (; it != it_end; ++it)
392  {
393  const Elem * const elem = *it;
394  libmesh_assert(elem);
395 
396  if (elem->on_boundary()) // denotes *any* side has a NULL neighbor
397  {
398  my_interface_elements.push_back(elem); // add the element, but only once, even
399  // if there are multiple NULL neighbors
400  for (unsigned int s=0; s<elem->n_sides(); s++)
401  if (elem->neighbor(s) == NULL)
402  {
403  AutoPtr<Elem> side(elem->build_side(s));
404 
405  for (unsigned int n=0; n<side->n_vertices(); n++)
406  my_interface_node_set.insert (side->node(n));
407  }
408  }
409  }
410 
411  my_interface_node_list.reserve (my_interface_node_set.size());
412  my_interface_node_list.insert (my_interface_node_list.end(),
413  my_interface_node_set.begin(),
414  my_interface_node_set.end());
415  }
416 
417  // we will now send my_interface_node_list to all of the adjacent processors.
418  // note that for the time being we will copy the list to a unique buffer for
419  // each processor so that we can use a nonblocking send and not access the
420  // buffer again until the send completes. it is my understanding that the
421  // MPI 2.1 standard seeks to remove this restriction as unnecessary, so in
422  // the future we should change this to send the same buffer to each of the
423  // adjacent processors. - BSK 11/17/2008
424  std::vector<std::vector<dof_id_type> >
425  my_interface_node_xfer_buffers (n_adjacent_processors, my_interface_node_list);
426  std::map<processor_id_type, unsigned char> n_comm_steps;
427 
428  std::vector<Parallel::Request> send_requests (3*n_adjacent_processors);
429  unsigned int current_request = 0;
430 
431  for (unsigned int comm_step=0; comm_step<n_adjacent_processors; comm_step++)
432  {
433  n_comm_steps[adjacent_processors[comm_step]]=1;
434  mesh.comm().send (adjacent_processors[comm_step],
435  my_interface_node_xfer_buffers[comm_step],
436  send_requests[current_request++],
437  element_neighbors_tag);
438  }
439 
440  //-------------------------------------------------------------------------
441  // processor pairings are symmetric - I expect to receive an interface node
442  // list from each processor in adjacent_processors as well!
443  // now we will catch an incoming node list for each of our adjacent processors.
444  //
445  // we are done with the adjacent_processors list - note that it is in general
446  // a superset of the processors we truly share elements with. so let's
447  // clear the superset list, and we will fill it with the true list.
448  adjacent_processors.clear();
449 
450  std::vector<dof_id_type> common_interface_node_list;
451 
452  // we expect two classess of messages -
453  // (1) incoming interface node lists, to which we will reply with our elements
454  // touching nodes in the list, and
455  // (2) replies from the requests we sent off previously.
456  // (2.a) - nodes
457  // (2.b) - elements
458  // so we expect 3 communications from each adjacent processor.
459  // by structuring the communication in this way we hopefully impose no
460  // order on the handling of the arriving messages. in particular, we
461  // should be able to handle the case where we receive a request and
462  // all replies from processor A before even receiving a request from
463  // processor B.
464 
465  for (unsigned int comm_step=0; comm_step<3*n_adjacent_processors; comm_step++)
466  {
467  //------------------------------------------------------------------
468  // catch incoming node list
469  Parallel::Status
471  element_neighbors_tag));
472  const processor_id_type
473  source_pid_idx = status.source(),
474  dest_pid_idx = source_pid_idx;
475 
476  //------------------------------------------------------------------
477  // first time - incoming request
478  if (n_comm_steps[source_pid_idx] == 1)
479  {
480  n_comm_steps[source_pid_idx]++;
481 
482  mesh.comm().receive (source_pid_idx,
483  common_interface_node_list,
484  element_neighbors_tag);
485  const std::size_t
486  their_interface_node_list_size = common_interface_node_list.size();
487 
488  // we now have the interface node list from processor source_pid_idx.
489  // now we can find all of our elements which touch any of these nodes
490  // and send copies back to this processor. however, we can make our
491  // search more efficient by first excluding all the nodes in
492  // their list which are not also contained in
493  // my_interface_node_list. we can do this in place as a set
494  // intersection.
495  common_interface_node_list.erase
496  (std::set_intersection (my_interface_node_list.begin(),
497  my_interface_node_list.end(),
498  common_interface_node_list.begin(),
499  common_interface_node_list.end(),
500  common_interface_node_list.begin()),
501  common_interface_node_list.end());
502 
503  if (false)
504  libMesh::out << "[" << mesh.processor_id() << "] "
505  << "my_interface_node_list.size()=" << my_interface_node_list.size()
506  << ", [" << source_pid_idx << "] "
507  << "their_interface_node_list.size()=" << their_interface_node_list_size
508  << ", common_interface_node_list.size()=" << common_interface_node_list.size()
509  << std::endl;
510 
511  // Now we need to see which of our elements touch the nodes in the list.
512  // We will certainly send all the active elements which intersect source_pid_idx,
513  // but we will also ship off the other elements in the same family tree
514  // as the active ones for data structure consistency.
515  //
516  // FIXME - shipping full family trees is unnecessary and inefficient.
517  //
518  // We also ship any nodes connected to these elements. Note
519  // some of these nodes and elements may be replicated from
520  // other processors, but that is OK.
521  std::set<const Elem*, CompareElemIdsByLevel> elements_to_send;
522  std::set<const Node*> connected_nodes;
523 
524  // Check for quick return?
525  if (common_interface_node_list.empty())
526  {
527  // let's try to be smart here - if we have no nodes in common,
528  // we cannot share elements. so post the messages expected
529  // from us here and go on about our business.
530  // note that even though these are nonblocking sends
531  // they should complete essentially instantly, because
532  // in all cases the send buffers are empty
533  mesh.comm().send_packed_range (dest_pid_idx,
534  &mesh,
535  connected_nodes.begin(),
536  connected_nodes.end(),
537  send_requests[current_request++],
538  element_neighbors_tag);
539 
540  mesh.comm().send_packed_range (dest_pid_idx,
541  &mesh,
542  elements_to_send.begin(),
543  elements_to_send.end(),
544  send_requests[current_request++],
545  element_neighbors_tag);
546 
547  continue;
548  }
549  // otherwise, this really *is* an adjacent processor.
550  adjacent_processors.push_back(source_pid_idx);
551 
552  std::vector<const Elem*> family_tree;
553 
554  for (dof_id_type e=0, n_shared_nodes=0; e<my_interface_elements.size(); e++, n_shared_nodes=0)
555  {
556  const Elem * elem = my_interface_elements[e];
557 
558  for (unsigned int n=0; n<elem->n_vertices(); n++)
559  if (std::binary_search (common_interface_node_list.begin(),
560  common_interface_node_list.end(),
561  elem->node(n)))
562  {
563  n_shared_nodes++;
564 
565  // TBD - how many nodes do we need to share
566  // before we care? certainly 2, but 1? not
567  // sure, so let's play it safe...
568  if (n_shared_nodes > 0) break;
569  }
570 
571  if (n_shared_nodes) // share at least one node?
572  {
573  elem = elem->top_parent();
574 
575  // avoid a lot of duplicated effort -- if we already have elem
576  // in the set its entire family tree is already in the set.
577  if (!elements_to_send.count(elem))
578  {
579 #ifdef LIBMESH_ENABLE_AMR
580  elem->family_tree(family_tree);
581 #else
582  family_tree.clear();
583  family_tree.push_back(elem);
584 #endif
585  for (unsigned int leaf=0; leaf<family_tree.size(); leaf++)
586  {
587  elem = family_tree[leaf];
588  elements_to_send.insert (elem);
589 
590  for (unsigned int n=0; n<elem->n_nodes(); n++)
591  connected_nodes.insert (elem->get_node(n));
592  }
593  }
594  }
595  }
596 
597  // The elements_to_send and connected_nodes sets now contain all
598  // the elements and nodes we need to send to this processor.
599  // All that remains is to pack up the objects (along with
600  // any boundary conditions) and send the messages off.
601  {
602  libmesh_assert (connected_nodes.empty() || !elements_to_send.empty());
603  libmesh_assert (!connected_nodes.empty() || elements_to_send.empty());
604 
605  // send the nodes off to the destination processor
606  mesh.comm().send_packed_range (dest_pid_idx,
607  &mesh,
608  connected_nodes.begin(),
609  connected_nodes.end(),
610  send_requests[current_request++],
611  element_neighbors_tag);
612 
613  // send the elements off to the destination processor
614  mesh.comm().send_packed_range (dest_pid_idx,
615  &mesh,
616  elements_to_send.begin(),
617  elements_to_send.end(),
618  send_requests[current_request++],
619  element_neighbors_tag);
620  }
621  }
622  //------------------------------------------------------------------
623  // second time - reply of nodes
624  else if (n_comm_steps[source_pid_idx] == 2)
625  {
626  n_comm_steps[source_pid_idx]++;
627 
628  mesh.comm().receive_packed_range (source_pid_idx,
629  &mesh,
630  mesh_inserter_iterator<Node>(mesh),
631  element_neighbors_tag);
632  }
633  //------------------------------------------------------------------
634  // third time - reply of elements
635  else if (n_comm_steps[source_pid_idx] == 3)
636  {
637  n_comm_steps[source_pid_idx]++;
638 
639  mesh.comm().receive_packed_range (source_pid_idx,
640  &mesh,
641  mesh_inserter_iterator<Elem>(mesh),
642  element_neighbors_tag);
643  }
644  //------------------------------------------------------------------
645  // fourth time - shouldn't happen
646  else
647  {
648  libMesh::err << "ERROR: unexpected number of replies: "
649  << n_comm_steps[source_pid_idx]
650  << std::endl;
651  }
652  } // done catching & processing replies associated with tag ~ 100,000pi
653 
654  // allow any pending requests to complete
655  Parallel::wait (send_requests);
656 
657  STOP_LOG("gather_neighboring_elements()","MeshCommunication");
658 }
659 #endif // LIBMESH_HAVE_MPI
660 
661 
662 #ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
663 // ------------------------------------------------------------
665 {
666  // no MPI == one processor, no need for this method...
667  return;
668 }
669 #else
670 // ------------------------------------------------------------
671 void MeshCommunication::broadcast (MeshBase& mesh) const
672 {
673  // Don't need to do anything if there is
674  // only one processor.
675  if (mesh.n_processors() == 1)
676  return;
677 
678  // This function must be run on all processors at once
679  libmesh_parallel_only(mesh.comm());
680 
681  START_LOG("broadcast()","MeshCommunication");
682 
683  // Explicitly clear the mesh on all but processor 0.
684  if (mesh.processor_id() != 0)
685  mesh.clear();
686 
687  // Broadcast nodes
688  mesh.comm().broadcast_packed_range(&mesh,
689  mesh.nodes_begin(),
690  mesh.nodes_end(),
691  &mesh,
693 
694  // Broadcast elements from coarsest to finest, so that child
695  // elements will see their parents already in place.
696  unsigned int n_levels = MeshTools::n_levels(mesh);
697  mesh.comm().broadcast(n_levels);
698 
699  for (unsigned int l=0; l != n_levels; ++l)
700  mesh.comm().broadcast_packed_range(&mesh,
701  mesh.level_elements_begin(l),
702  mesh.level_elements_end(l),
703  &mesh,
705 
706  // Make sure mesh dimension is consistent
707  unsigned int mesh_dimension = mesh.mesh_dimension();
708  mesh.comm().broadcast(mesh_dimension);
709  mesh.set_mesh_dimension(mesh_dimension);
710 
711  // Broadcast all of the named entity information
712  mesh.comm().broadcast(mesh.set_subdomain_name_map());
713  mesh.comm().broadcast(mesh.boundary_info->set_sideset_name_map());
714  mesh.comm().broadcast(mesh.boundary_info->set_nodeset_name_map());
715 
716  libmesh_assert (mesh.comm().verify(mesh.n_elem()));
717  libmesh_assert (mesh.comm().verify(mesh.n_nodes()));
718 
719 #ifdef DEBUG
722 #endif
723 
724  STOP_LOG("broadcast()","MeshCommunication");
725 }
726 #endif // LIBMESH_HAVE_MPI
727 
728 
729 
730 #ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
731 // ------------------------------------------------------------
733 {
734  // no MPI == one processor, no need for this method...
735  return;
736 }
737 #else
738 // ------------------------------------------------------------
739 void MeshCommunication::gather (const processor_id_type root_id, ParallelMesh& mesh) const
740 {
741  // The mesh should know it's about to be serialized
742  libmesh_assert (mesh.is_serial());
743 
744  // Check for quick return
745  if (mesh.n_processors() == 1)
746  return;
747 
748  // This function must be run on all processors at once
749  libmesh_parallel_only(mesh.comm());
750 
751  START_LOG("(all)gather()","MeshCommunication");
752 
753  (root_id == DofObject::invalid_processor_id) ?
754 
755  mesh.comm().allgather_packed_range (&mesh,
756  mesh.nodes_begin(),
757  mesh.nodes_end(),
759 
760  mesh.comm().gather_packed_range (root_id,
761  &mesh,
762  mesh.nodes_begin(),
763  mesh.nodes_end(),
765 
766  // Gather elements from coarsest to finest, so that child
767  // elements will see their parents already in place.
768  // rank 0 should know n_levels regardless, so this is
769  // safe independent of root_id
770  unsigned int n_levels = MeshTools::n_levels(mesh);
771  mesh.comm().broadcast(n_levels);
772 
773  for (unsigned int l=0; l != n_levels; ++l)
774  (root_id == DofObject::invalid_processor_id) ?
775 
776  mesh.comm().allgather_packed_range (&mesh,
777  mesh.level_elements_begin(l),
778  mesh.level_elements_end(l),
780 
781  mesh.comm().gather_packed_range (root_id,
782  &mesh,
783  mesh.level_elements_begin(l),
784  mesh.level_elements_end(l),
786 
787 
788  // If we are doing an allgather(), perform sanity check on the result.
789  if (root_id == DofObject::invalid_processor_id)
790  {
791  libmesh_assert (mesh.comm().verify(mesh.n_elem()));
792  libmesh_assert (mesh.comm().verify(mesh.n_nodes()));
793  }
794 
795  // Inform new elements of their neighbors,
796  // while resetting all remote_elem links on
797  // the appropriate ranks.
798  if ((root_id == DofObject::invalid_processor_id) ||
799  (mesh.comm().rank() == root_id))
800  mesh.find_neighbors(true);
801 
802  // All done!
803  STOP_LOG("(all)gather()","MeshCommunication");
804 }
805 #endif // LIBMESH_HAVE_MPI
806 
807 
808 
809 // Functor for make_elems_parallel_consistent and
810 // make_node_ids_parallel_consistent
811 namespace {
812 
813 struct SyncIds
814 {
815 typedef dof_id_type datum;
816 typedef void (MeshBase::*renumber_obj)(dof_id_type, dof_id_type);
817 
818 SyncIds(MeshBase &_mesh, renumber_obj _renumberer) :
819  mesh(_mesh),
820  renumber(_renumberer) {}
821 
822 MeshBase &mesh;
823 renumber_obj renumber;
824 // renumber_obj &renumber;
825 
826 // Find the id of each requested DofObject -
827 // sync_dofobject_data_by_xyz already did the work for us
828 void gather_data (const std::vector<dof_id_type>& ids,
829  std::vector<datum>& ids_out)
830 {
831  ids_out = ids;
832 }
833 
834 void act_on_data (const std::vector<dof_id_type>& old_ids,
835  std::vector<datum>& new_ids)
836 {
837  for (unsigned int i=0; i != old_ids.size(); ++i)
838  if (old_ids[i] != new_ids[i])
839  (mesh.*renumber)(old_ids[i], new_ids[i]);
840 }
841 };
842 }
843 
844 
845 
846 // ------------------------------------------------------------
848  (MeshBase &mesh,
849  LocationMap<Node> &loc_map)
850 {
851  // This function must be run on all processors at once
852  libmesh_parallel_only(mesh.comm());
853 
854  START_LOG ("make_node_ids_parallel_consistent()", "MeshCommunication");
855 
856  SyncIds syncids(mesh, &MeshBase::renumber_node);
858  (mesh.comm(),
859  mesh.nodes_begin(), mesh.nodes_end(),
860  loc_map, syncids);
861 
862  STOP_LOG ("make_node_ids_parallel_consistent()", "MeshCommunication");
863 }
864 
865 
866 
867 // ------------------------------------------------------------
869 {
870  // This function must be run on all processors at once
871  libmesh_parallel_only(mesh.comm());
872 
873  START_LOG ("make_elems_parallel_consistent()", "MeshCommunication");
874 
875  SyncIds syncids(mesh, &MeshBase::renumber_elem);
877  (mesh, mesh.active_elements_begin(),
878  mesh.active_elements_end(), syncids);
879 
880  STOP_LOG ("make_elems_parallel_consistent()", "MeshCommunication");
881 }
882 
883 
884 
885 // Functors for make_node_proc_ids_parallel_consistent
886 namespace {
887 
888 struct SyncProcIds
889 {
890 typedef processor_id_type datum;
891 
892 SyncProcIds(MeshBase &_mesh) : mesh(_mesh) {}
893 
894 MeshBase &mesh;
895 
896 // ------------------------------------------------------------
897 void gather_data (const std::vector<dof_id_type>& ids,
898  std::vector<datum>& data)
899 {
900  // Find the processor id of each requested node
901  data.resize(ids.size());
902 
903  for (std::size_t i=0; i != ids.size(); ++i)
904  {
905  // Look for this point in the mesh
906  // We'd better find every node we're asked for
907  Node *node = mesh.node_ptr(ids[i]);
908 
909  // Return the node's correct processor id,
910  data[i] = node->processor_id();
911  }
912 }
913 
914 // ------------------------------------------------------------
915 void act_on_data (const std::vector<dof_id_type>& ids,
916  std::vector<datum> proc_ids)
917 {
918  // Set the ghost node processor ids we've now been informed of
919  for (std::size_t i=0; i != ids.size(); ++i)
920  {
921  Node *node = mesh.node_ptr(ids[i]);
922  node->processor_id() = proc_ids[i];
923  }
924 }
925 };
926 }
927 
928 
929 
930 // ------------------------------------------------------------
932  (MeshBase& mesh,
933  LocationMap<Node>& loc_map)
934 {
935  START_LOG ("make_node_proc_ids_parallel_consistent()", "MeshCommunication");
936 
937  // This function must be run on all processors at once
938  libmesh_parallel_only(mesh.comm());
939 
940  // When this function is called, each section of a parallelized mesh
941  // should be in the following state:
942  //
943  // All nodes should have the exact same physical location on every
944  // processor where they exist.
945  //
946  // Local nodes should have unique authoritative ids,
947  // and processor ids consistent with all processors which own
948  // an element touching them.
949  //
950  // Ghost nodes touching local elements should have processor ids
951  // consistent with all processors which own an element touching
952  // them.
953 
954  SyncProcIds sync(mesh);
956  (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), loc_map, sync);
957 
958  STOP_LOG ("make_node_proc_ids_parallel_consistent()", "MeshCommunication");
959 }
960 
961 
962 
963 // ------------------------------------------------------------
965  (MeshBase &mesh,
966  LocationMap<Node> &loc_map)
967 {
968  // This function must be run on all processors at once
969  libmesh_parallel_only(mesh.comm());
970 
971  // Create the loc_map if it hasn't been done already
972  bool need_map_update = (mesh.nodes_begin() != mesh.nodes_end() &&
973  loc_map.empty());
974  mesh.comm().max(need_map_update);
975 
976  if (need_map_update)
977  loc_map.init(mesh);
978 
979  // When this function is called, each section of a parallelized mesh
980  // should be in the following state:
981  //
982  // All nodes should have the exact same physical location on every
983  // processor where they exist.
984  //
985  // Local nodes should have unique authoritative ids,
986  // and processor ids consistent with all processors which own
987  // an element touching them.
988  //
989  // Ghost nodes touching local elements should have processor ids
990  // consistent with all processors which own an element touching
991  // them.
992  //
993  // Ghost nodes should have ids which are either already correct
994  // or which are in the "unpartitioned" id space.
995 
996  // First, let's sync up processor ids. Some of these processor ids
997  // may be "wrong" from coarsening, but they're right in the sense
998  // that they'll tell us who has the authoritative dofobject ids for
999  // each node.
1000  this->make_node_proc_ids_parallel_consistent(mesh, loc_map);
1001 
1002  // Second, sync up dofobject ids.
1003  this->make_node_ids_parallel_consistent(mesh, loc_map);
1004 
1005  // Finally, correct the processor ids to make DofMap happy
1006  MeshTools::correct_node_proc_ids(mesh, loc_map);
1007 }
1008 
1009 
1010 
1011 // ------------------------------------------------------------
1012 void MeshCommunication::delete_remote_elements(ParallelMesh& mesh, const std::set<Elem *> & extra_ghost_elem_ids) const
1013 {
1014  // The mesh should know it's about to be parallelized
1015  libmesh_assert (!mesh.is_serial());
1016 
1017  START_LOG("delete_remote_elements()", "MeshCommunication");
1018 
1019 #ifdef DEBUG
1020  // We expect maximum ids to be in sync so we can use them to size
1021  // vectors
1022  mesh.comm().verify(mesh.max_node_id());
1023  mesh.comm().verify(mesh.max_elem_id());
1024  const dof_id_type par_max_node_id = mesh.parallel_max_node_id();
1025  const dof_id_type par_max_elem_id = mesh.parallel_max_elem_id();
1026  libmesh_assert_equal_to (par_max_node_id, mesh.max_node_id());
1027  libmesh_assert_equal_to (par_max_elem_id, mesh.max_elem_id());
1028 #endif
1029 
1030  // FIXME - should these be "unsorted_set"s? O(N) is O(N)...
1031  std::vector<bool> local_nodes(mesh.max_node_id(), false);
1032  std::vector<bool> semilocal_nodes(mesh.max_node_id(), false);
1033  std::vector<bool> semilocal_elems(mesh.max_elem_id(), false);
1034 
1035  // We don't want to delete any element that shares a node
1036  // with or is an ancestor of a local element.
1038  l_end = mesh.local_elements_end();
1039  for (; l_elem_it != l_end; ++l_elem_it)
1040  {
1041  const Elem *elem = *l_elem_it;
1042  for (unsigned int n=0; n != elem->n_nodes(); ++n)
1043  {
1044  dof_id_type nodeid = elem->node(n);
1045  libmesh_assert_less (nodeid, local_nodes.size());
1046  local_nodes[nodeid] = true;
1047  }
1048  while (elem)
1049  {
1050  dof_id_type elemid = elem->id();
1051  libmesh_assert_less (elemid, semilocal_elems.size());
1052  semilocal_elems[elemid] = true;
1053 
1054  for (unsigned int n=0; n != elem->n_nodes(); ++n)
1055  semilocal_nodes[elem->node(n)] = true;
1056 
1057  const Elem *parent = elem->parent();
1058  // Don't proceed from a boundary mesh to an interior mesh
1059  if (parent && parent->dim() != elem->dim())
1060  break;
1061 
1062  elem = parent;
1063  }
1064  }
1065 
1066  // We don't want to delete any element that shares a node
1067  // with or is an ancestor of an unpartitioned element either.
1069  u_elem_it = mesh.unpartitioned_elements_begin(),
1070  u_end = mesh.unpartitioned_elements_end();
1071 
1072  for (; u_elem_it != u_end; ++u_elem_it)
1073  {
1074  const Elem *elem = *u_elem_it;
1075  for (unsigned int n=0; n != elem->n_nodes(); ++n)
1076  local_nodes[elem->node(n)] = true;
1077  while (elem)
1078  {
1079  semilocal_elems[elem->id()] = true;
1080 
1081  for (unsigned int n=0; n != elem->n_nodes(); ++n)
1082  semilocal_nodes[elem->node(n)] = true;
1083 
1084  const Elem *parent = elem->parent();
1085  // Don't proceed from a boundary mesh to an interior mesh
1086  if (parent && parent->dim() != elem->dim())
1087  break;
1088 
1089  elem = parent;
1090  }
1091  }
1092 
1093  // Flag all the elements that share nodes with
1094  // local and unpartitioned elements, along with their ancestors
1096  nl_end = mesh.not_local_elements_end();
1097  for (; nl_elem_it != nl_end; ++nl_elem_it)
1098  {
1099  const Elem *elem = *nl_elem_it;
1100  for (unsigned int n=0; n != elem->n_nodes(); ++n)
1101  if (local_nodes[elem->node(n)])
1102  {
1103  while (elem)
1104  {
1105  semilocal_elems[elem->id()] = true;
1106 
1107  for (unsigned int nn=0; nn != elem->n_nodes(); ++nn)
1108  semilocal_nodes[elem->node(nn)] = true;
1109 
1110  const Elem *parent = elem->parent();
1111  // Don't proceed from a boundary mesh to an interior mesh
1112  if (parent && parent->dim() != elem->dim())
1113  break;
1114 
1115  elem = parent;
1116  }
1117  break;
1118  }
1119  }
1120 
1121  // Don't delete elements that we were explicitly told not to
1122  for(std::set<Elem *>::iterator it = extra_ghost_elem_ids.begin();
1123  it != extra_ghost_elem_ids.end();
1124  ++it)
1125  {
1126  const Elem *elem = *it;
1127  semilocal_elems[elem->id()] = true;
1128  for (unsigned int n=0; n != elem->n_nodes(); ++n)
1129  semilocal_nodes[elem->node(n)] = true;
1130  }
1131 
1132  // Delete all the elements we have no reason to save,
1133  // starting with the most refined so that the mesh
1134  // is valid at all intermediate steps
1135  unsigned int n_levels = MeshTools::n_levels(mesh);
1136 
1137  for (int l = n_levels - 1; l >= 0; --l)
1138  {
1139  MeshBase::element_iterator lev_elem_it = mesh.level_elements_begin(l),
1140  lev_end = mesh.level_elements_end(l);
1141  for (; lev_elem_it != lev_end; ++lev_elem_it)
1142  {
1143  Elem *elem = *lev_elem_it;
1144  libmesh_assert (elem);
1145  // Make sure we don't leave any invalid pointers
1146  if (!semilocal_elems[elem->id()])
1147  elem->make_links_to_me_remote();
1148 
1149  // Subactive neighbor pointers aren't preservable here
1150  if (elem->subactive())
1151  for (unsigned int s=0; s != elem->n_sides(); ++s)
1152  elem->set_neighbor(s, NULL);
1153 
1154  // delete_elem doesn't currently invalidate element
1155  // iterators... that had better not change
1156  if (!semilocal_elems[elem->id()])
1157  mesh.delete_elem(elem);
1158  }
1159  }
1160 
1161  // Delete all the nodes we have no reason to save
1162  MeshBase::node_iterator node_it = mesh.nodes_begin(),
1163  node_end = mesh.nodes_end();
1164  for (node_it = mesh.nodes_begin(); node_it != node_end; ++node_it)
1165  {
1166  Node *node = *node_it;
1167  libmesh_assert(node);
1168  if (!semilocal_nodes[node->id()])
1169  mesh.delete_node(node);
1170  }
1171 
1172 #ifdef DEBUG
1174 #endif
1175 
1176  STOP_LOG("delete_remote_elements()", "MeshCommunication");
1177 }
1178 
1179 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo