Home
Publications
Presentations
Support
Download
Examples
Developers
Class Docs
Namespaces
Classes
Files
Mailing Lists
Gallery
Wiki

mesh_communication.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2012 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



// C++ Includes   -----------------------------------
#include <numeric>

// Local Includes -----------------------------------
#include "libmesh/boundary_info.h"
#include "libmesh/elem.h"
#include "libmesh/libmesh_config.h"
#include "libmesh/libmesh_common.h"
#include "libmesh/libmesh_logging.h"
#include "libmesh/location_maps.h"
#include "libmesh/mesh_base.h"
#include "libmesh/mesh_communication.h"
#include "libmesh/mesh_inserter_iterator.h"
#include "libmesh/mesh_tools.h"
#include "libmesh/parallel.h"
#include "libmesh/parallel_mesh.h"
#include "libmesh/parallel_ghost_sync.h"
#include "libmesh/utility.h"
#include "libmesh/remote_elem.h"



//-----------------------------------------------
// anonymous namespace for implementation details
namespace {

  using libMesh::Elem;

  struct CompareElemIdsByLevel
  {
    bool operator()(const Elem *a,
                    const Elem *b) const
    {
      libmesh_assert (a);
      libmesh_assert (b);
      const unsigned int
        al = a->level(), bl = b->level();
      const dof_id_type
        aid = a->id(),   bid = b->id();

      return (al == bl) ? aid < bid : al < bl;
    }
  };
}


namespace libMesh
{


// ------------------------------------------------------------
// MeshCommunication class members
void MeshCommunication::clear ()
{
  //  _neighboring_processors.clear();
}


// #ifdef LIBMESH_HAVE_MPI
// void MeshCommunication::find_neighboring_processors (const MeshBase& mesh)
// {
//   // Don't need to do anything if there is
//   // only one processor.
//   if (libMesh::n_processors() == 1)
//     return;

//   _neighboring_processors.clear();

//   // Get the bounding sphere for the local processor
//   Sphere bounding_sphere =
//     MeshTools::processor_bounding_sphere (mesh, libMesh::processor_id());

//   // Just to be sure, increase its radius by 10%.  Sure would suck to
//   // miss a neighboring processor!
//   bounding_sphere.radius() *= 1.1;

//   // Collect the bounding spheres from all processors, test for intersection
//   {
//     std::vector<float>
//       send (4,                         0),
//       recv (4*libMesh::n_processors(), 0);

//     send[0] = bounding_sphere.center()(0);
//     send[1] = bounding_sphere.center()(1);
//     send[2] = bounding_sphere.center()(2);
//     send[3] = bounding_sphere.radius();

//     MPI_Allgather (&send[0], send.size(), MPI_FLOAT,
//                 &recv[0], send.size(), MPI_FLOAT,
//                 libMesh::COMM_WORLD);


//     for (processor_id_type proc=0; proc<libMesh::n_processors(); proc++)
//       {
//      const Point center (recv[4*proc+0],
//                          recv[4*proc+1],
//                          recv[4*proc+2]);

//      const Real radius = recv[4*proc+3];

//      const Sphere proc_sphere (center, radius);

//      if (bounding_sphere.intersects(proc_sphere))
//        _neighboring_processors.push_back(proc);
//       }

//     // Print out the _neighboring_processors list
//     libMesh::out << "Processor " << libMesh::processor_id()
//            << " intersects:" << std::endl;
//     for (unsigned int p=0; p<_neighboring_processors.size(); p++)
//       libMesh::out << " " << _neighboring_processors[p] << std::endl;
//   }
// }
// #else
// void MeshCommunication::find_neighboring_processors (const MeshBase&)
// {
// }
// #endif



#ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
// ------------------------------------------------------------
void MeshCommunication::redistribute (ParallelMesh &) const
{
  // no MPI == one processor, no redistribution
  return;
}
#else
// ------------------------------------------------------------
void MeshCommunication::redistribute (ParallelMesh &mesh) const
{
  // This method will be called after a new partitioning has been
  // assigned to the elements.  This partitioning was defined in
  // terms of the active elements, and "trickled down" to the
  // parents and nodes as to be consistent.
  //
  // The point is that the entire concept of local elements is
  // kinda shaky in this method.  Elements which were previously
  // local may now be assigned to other processors, so we need to
  // send those off.  Similarly, we need to accept elements from
  // other processors.
  //
  // The approach is as follows:
  // (1) send all elements we have stored to their proper homes
  // (2) receive elements from all processors, watching for duplicates
  // (3) deleting all nonlocal elements elements
  // (4) obtaining required ghost elements from neighboring processors
  parallel_only();
  libmesh_assert (!mesh.is_serial());
  libmesh_assert (MeshTools::n_elem(mesh.unpartitioned_elements_begin(),
                                    mesh.unpartitioned_elements_end()) == 0);

  START_LOG("redistribute()","MeshCommunication");

  // Get a few unique message tags to use in communications; we'll
  // default to some numbers around pi*1000
  Parallel::MessageTag
    nodestag   = Parallel::Communicator_World.get_unique_tag(3141),
    elemstag   = Parallel::Communicator_World.get_unique_tag(3142);

  // Figure out how many nodes and elements we have which are assigned to each
  // processor.  send_n_nodes_and_elem_per_proc contains the number of nodes/elements
  // we will be sending to each processor, recv_n_nodes_and_elem_per_proc contains
  // the number of nodes/elements we will be receiving from each processor.
  // Format:
  //  send_n_nodes_and_elem_per_proc[2*pid+0] = number of nodes to send to pid
  //  send_n_nodes_and_elem_per_proc[2*pid+1] = number of elements to send to pid
  std::vector<dof_id_type> send_n_nodes_and_elem_per_proc(2*libMesh::n_processors(), 0);

  std::vector<Parallel::Request>
    node_send_requests, element_send_requests;

  for (processor_id_type pid=0; pid<libMesh::n_processors(); pid++)
    if (pid != libMesh::processor_id()) // don't send to ourselves!!
      {
        // Build up a list of nodes and elements to send to processor pid.
        // We will certainly send all the elements assigned to this processor,
        // but we will also ship off any other elements which touch
        // their nodes.
        std::set<const Node*> connected_nodes;
        {
          MeshBase::const_element_iterator       elem_it  = mesh.pid_elements_begin(pid);
          const MeshBase::const_element_iterator elem_end = mesh.pid_elements_end(pid);

          for (; elem_it!=elem_end; ++elem_it)
            {
              const Elem* elem = *elem_it;

              for (unsigned int n=0; n<elem->n_nodes(); n++)
                connected_nodes.insert (elem->get_node(n));
            }
        }

        std::set<const Elem*, CompareElemIdsByLevel> elements_to_send;
        {
          MeshBase::const_element_iterator       elem_it  = mesh.elements_begin();
          const MeshBase::const_element_iterator elem_end = mesh.elements_end();

          for (; elem_it!=elem_end; ++elem_it)
            {
              const Elem* elem = *elem_it;

              for (unsigned int n=0; n<elem->n_nodes(); n++)
                if (connected_nodes.count(elem->get_node(n)))
                  elements_to_send.insert (elem);
            }
        }

        connected_nodes.clear();
        {
          std::set<const Elem*, CompareElemIdsByLevel>::iterator
            elem_it  = elements_to_send.begin(),
            elem_end = elements_to_send.end();

          for (; elem_it!=elem_end; ++elem_it)
            {
              const Elem* elem = *elem_it;

              for (unsigned int n=0; n<elem->n_nodes(); n++)
                connected_nodes.insert(elem->get_node(n));
            }
        }

        // the number of nodes we will ship to pid
        send_n_nodes_and_elem_per_proc[2*pid+0] = connected_nodes.size();

        // send any nodes off to the destination processor
        if (!connected_nodes.empty())
          {
            node_send_requests.push_back(Parallel::request());

            CommWorld.send_packed_range (pid,
                                         &mesh,
                                         connected_nodes.begin(),
                                         connected_nodes.end(),
                                         node_send_requests.back(),
                                         nodestag);
          }

        // the number of elements we will send to this processor
        send_n_nodes_and_elem_per_proc[2*pid+1] = elements_to_send.size();

        if (!elements_to_send.empty())
          {
            // send the elements off to the destination processor
            element_send_requests.push_back(Parallel::request());

            CommWorld.send_packed_range (pid,
                                         &mesh,
                                         elements_to_send.begin(),
                                         elements_to_send.end(),
                                         element_send_requests.back(),
                                         elemstag);
          }
      }

  std::vector<dof_id_type> recv_n_nodes_and_elem_per_proc(send_n_nodes_and_elem_per_proc);

  CommWorld.alltoall (recv_n_nodes_and_elem_per_proc);

  // In general we will only need to communicate with a subset of the other processors.
  // I can't immediately think of a case where we will send elements but not nodes, but
  // these are only bools and we have the information anyway...
  std::vector<bool>
    send_node_pair(libMesh::n_processors(),false), send_elem_pair(libMesh::n_processors(),false),
    recv_node_pair(libMesh::n_processors(),false), recv_elem_pair(libMesh::n_processors(),false);

  unsigned int
    n_send_node_pairs=0,      n_send_elem_pairs=0,
    n_recv_node_pairs=0,      n_recv_elem_pairs=0;

  for (processor_id_type pid=0; pid<libMesh::n_processors(); pid++)
    {
      if (send_n_nodes_and_elem_per_proc[2*pid+0]) // we have nodes to send
        {
          send_node_pair[pid] = true;
          n_send_node_pairs++;
        }

      if (send_n_nodes_and_elem_per_proc[2*pid+1]) // we have elements to send
        {
          send_elem_pair[pid] = true;
          n_send_elem_pairs++;
        }

      if (recv_n_nodes_and_elem_per_proc[2*pid+0]) // we have nodes to receive
        {
          recv_node_pair[pid] = true;
          n_recv_node_pairs++;
        }

      if (recv_n_nodes_and_elem_per_proc[2*pid+1]) // we have elements to receive
        {
          recv_elem_pair[pid] = true;
          n_recv_elem_pairs++;
        }
    }
  libmesh_assert_equal_to (n_send_node_pairs, node_send_requests.size());
  libmesh_assert_equal_to (n_send_elem_pairs, element_send_requests.size());

  // Receive nodes.
  // We now know how many processors will be sending us nodes.
  for (unsigned int node_comm_step=0; node_comm_step<n_recv_node_pairs; node_comm_step++)
    // but we don't necessarily want to impose an ordering, so
    // just process whatever message is available next.
    CommWorld.receive_packed_range (Parallel::any_source,
                                    &mesh,
                                    mesh_inserter_iterator<Node>(mesh),
                                    nodestag);

  // Receive elements.
  // Similarly we know how many processors are sending us elements,
  // but we don't really care in what order we receive them.
  for (unsigned int elem_comm_step=0; elem_comm_step<n_recv_elem_pairs; elem_comm_step++)
    CommWorld.receive_packed_range (Parallel::any_source,
                                    &mesh,
                                    mesh_inserter_iterator<Elem>(mesh),
                                    elemstag);

  // Wait for all sends to complete
  Parallel::wait (node_send_requests);
  Parallel::wait (element_send_requests);

  // Check on the redistribution consistency
#ifdef DEBUG
  MeshTools::libmesh_assert_equal_n_systems(mesh);

  MeshTools::libmesh_assert_valid_refinement_tree(mesh);
#endif

  STOP_LOG("redistribute()","MeshCommunication");
}
#endif // LIBMESH_HAVE_MPI



#ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
// ------------------------------------------------------------
void MeshCommunication::gather_neighboring_elements (ParallelMesh &) const
{
  // no MPI == one processor, no need for this method...
  return;
}
#else
// ------------------------------------------------------------
void MeshCommunication::gather_neighboring_elements (ParallelMesh &mesh) const
{
  // Don't need to do anything if there is
  // only one processor.
  if (libMesh::n_processors() == 1)
    return;

  // This function must be run on all processors at once
  parallel_only();

  START_LOG("gather_neighboring_elements()","MeshCommunication");

  //------------------------------------------------------------------
  // The purpose of this function is to provide neighbor data structure
  // consistency for a parallel, distributed mesh.  In libMesh we require
  // that each local element have access to a full set of valid face
  // neighbors.  In some cases this requires us to store "ghost elements" -
  // elements that belong to other processors but we store to provide
  // data structure consistency.  Also, it is assumed that any element
  // with a NULL neighbor resides on a physical domain boundary.  So,
  // even our "ghost elements" must have non-NULL neighbors.  To handle
  // this the concept of "RemoteElem" is used - a special construct which
  // is used to denote that an element has a face neighbor, but we do
  // not actually store detailed information about that neighbor.  This
  // is required to prevent data structure explosion.
  //
  // So when this method is called we should have only local elements.
  // These local elements will then find neighbors among the local
  // element set.  After this is completed, any element with a NULL
  // neighbor has either (i) a face on the physical boundary of the mesh,
  // or (ii) a neighboring element which lives on a remote processor.
  // To handle case (ii), we communicate the global node indices connected
  // to all such faces to our neighboring processors.  They then send us
  // all their elements with a NULL neighbor that are connected to any
  // of the nodes in our list.
  //------------------------------------------------------------------

  // Let's begin with finding consistent neighbor data information
  // for all the elements we currently have.  We'll use a clean
  // slate here - clear any existing information, including RemoteElem's.
//  mesh.find_neighbors (/* reset_remote_elements = */ true,
//                     /* reset_current_list    = */ true);

  // Get a unique message tag to use in communications; we'll default
  // to some numbers around pi*10000
  Parallel::MessageTag
    element_neighbors_tag = Parallel::Communicator_World.get_unique_tag(31416);

  // Now any element with a NULL neighbor either
  // (i) lives on the physical domain boundary, or
  // (ii) lives on an inter-processor boundary.
  // We will now gather all the elements from adjacent processors
  // which are of the same state, which should address all the type (ii)
  // elements.

  // A list of all the processors which *may* contain neighboring elements.
  // (for development simplicity, just make this the identity map)
  std::vector<processor_id_type> adjacent_processors;
  for (processor_id_type pid=0; pid<libMesh::n_processors(); pid++)
    if (pid != libMesh::processor_id())
      adjacent_processors.push_back (pid);


  const processor_id_type n_adjacent_processors = 
    libmesh_cast_int<processor_id_type>(adjacent_processors.size());

  //-------------------------------------------------------------------------
  // Let's build a list of all nodes which live on NULL-neighbor sides.
  // For simplicity, we will use a set to build the list, then transfer
  // it to a vector for communication.
  std::vector<dof_id_type> my_interface_node_list;
  std::vector<const Elem*>  my_interface_elements;
  {
    std::set<dof_id_type> my_interface_node_set;

    // since parent nodes are a subset of children nodes, this should be sufficient
    MeshBase::const_element_iterator       it     = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator it_end = mesh.active_local_elements_end();

    for (; it != it_end; ++it)
      {
        const Elem * const elem = *it;
        libmesh_assert(elem);

        if (elem->on_boundary()) // denotes *any* side has a NULL neighbor
          {
            my_interface_elements.push_back(elem); // add the element, but only once, even
                                                   // if there are multiple NULL neighbors
            for (unsigned int s=0; s<elem->n_sides(); s++)
              if (elem->neighbor(s) == NULL)
                {
                  AutoPtr<Elem> side(elem->build_side(s));

                  for (unsigned int n=0; n<side->n_vertices(); n++)
                    my_interface_node_set.insert (side->node(n));
                }
          }
      }

    my_interface_node_list.reserve (my_interface_node_set.size());
    my_interface_node_list.insert  (my_interface_node_list.end(),
                                    my_interface_node_set.begin(),
                                    my_interface_node_set.end());
  }

  // we will now send my_interface_node_list to all of the adjacent processors.
  // note that for the time being we will copy the list to a unique buffer for
  // each processor so that we can use a nonblocking send and not access the
  // buffer again until the send completes.  it is my understanding that the
  // MPI 2.1 standard seeks to remove this restriction as unnecessary, so in
  // the future we should change this to send the same buffer to each of the
  // adjacent processors. - BSK 11/17/2008
  std::vector<std::vector<dof_id_type> >
    my_interface_node_xfer_buffers (n_adjacent_processors, my_interface_node_list);
  std::map<processor_id_type, unsigned char> n_comm_steps;

  std::vector<Parallel::Request> send_requests (3*n_adjacent_processors);
  unsigned int current_request = 0;

  for (unsigned int comm_step=0; comm_step<n_adjacent_processors; comm_step++)
    {
      n_comm_steps[adjacent_processors[comm_step]]=1;
      CommWorld.send (adjacent_processors[comm_step],
                      my_interface_node_xfer_buffers[comm_step],
                      send_requests[current_request++],
                      element_neighbors_tag);
    }

  //-------------------------------------------------------------------------
  // processor pairings are symmetric - I expect to receive an interface node
  // list from each processor in adjacent_processors as well!
  // now we will catch an incoming node list for each of our adjacent processors.
  //
  // we are done with the adjacent_processors list - note that it is in general
  // a superset of the processors we truly share elements with.  so let's
  // clear the superset list, and we will fill it with the true list.
  adjacent_processors.clear();

  std::vector<dof_id_type> common_interface_node_list;

  // we expect two classess of messages -
  // (1) incoming interface node lists, to which we will reply with our elements
  //     touching nodes in the list, and
  // (2) replies from the requests we sent off previously.
  //  (2.a) - nodes
  //  (2.b) - elements
  // so we expect 3 communications from each adjacent processor.
  // by structuring the communication in this way we hopefully impose no
  // order on the handling of the arriving messages.  in particular, we
  // should be able to handle the case where we receive a request and
  // all replies from processor A before even receiving a request from
  // processor B.

  for (unsigned int comm_step=0; comm_step<3*n_adjacent_processors; comm_step++)
    {
      //------------------------------------------------------------------
      // catch incoming node list
      Parallel::Status
        status(CommWorld.probe (Parallel::any_source,
                                element_neighbors_tag));
      const processor_id_type
        source_pid_idx = status.source(),
        dest_pid_idx   = source_pid_idx;

      //------------------------------------------------------------------
      // first time - incoming request
      if (n_comm_steps[source_pid_idx] == 1)
        {
          n_comm_steps[source_pid_idx]++;

          CommWorld.receive (source_pid_idx,
                             common_interface_node_list,
                             element_neighbors_tag);
          const std::size_t
            their_interface_node_list_size = common_interface_node_list.size();

          // we now have the interface node list from processor source_pid_idx.
          // now we can find all of our elements which touch any of these nodes
          // and send copies back to this processor.  however, we can make our
          // search more efficient by first excluding all the nodes in
          // their list which are not also contained in
          // my_interface_node_list.  we can do this in place as a set
          // intersection.
          common_interface_node_list.erase
            (std::set_intersection (my_interface_node_list.begin(),
                                    my_interface_node_list.end(),
                                    common_interface_node_list.begin(),
                                    common_interface_node_list.end(),
                                    common_interface_node_list.begin()),
                                    common_interface_node_list.end());

          if (false)
            libMesh::out << "[" << libMesh::processor_id() << "] "
                          << "my_interface_node_list.size()="       << my_interface_node_list.size()
                          << ", [" << source_pid_idx << "] "
                          << "their_interface_node_list.size()="    << their_interface_node_list_size
                          << ", common_interface_node_list.size()=" << common_interface_node_list.size()
                          << std::endl;

          // Now we need to see which of our elements touch the nodes in the list.
          // We will certainly send all the active elements which intersect source_pid_idx,
          // but we will also ship off the other elements in the same family tree
          // as the active ones for data structure consistency.  
          //
          // FIXME - shipping full family trees is unnecessary and inefficient.
          //
          // We also ship any nodes connected to these elements.  Note
          // some of these nodes and elements may be replicated from
          // other processors, but that is OK.
          std::set<const Elem*, CompareElemIdsByLevel> elements_to_send;
          std::set<const Node*> connected_nodes;

          // Check for quick return?
          if (common_interface_node_list.empty())
            {
              // let's try to be smart here - if we have no nodes in common,
              // we cannot share elements.  so post the messages expected
              // from us here and go on about our business.
              // note that even though these are nonblocking sends
              // they should complete essentially instantly, because
              // in all cases the send buffers are empty
              CommWorld.send_packed_range (dest_pid_idx,
                                           &mesh,
                                           connected_nodes.begin(),
                                           connected_nodes.end(),
                                           send_requests[current_request++],
                                           element_neighbors_tag);

              CommWorld.send_packed_range (dest_pid_idx,
                                           &mesh,
                                           elements_to_send.begin(),
                                           elements_to_send.end(),
                                           send_requests[current_request++],
                                           element_neighbors_tag);

              continue;
            }
          // otherwise, this really *is* an adjacent processor.
          adjacent_processors.push_back(source_pid_idx);

          std::vector<const Elem*> family_tree;

          for (dof_id_type e=0, n_shared_nodes=0; e<my_interface_elements.size(); e++, n_shared_nodes=0)
            {
              const Elem * elem = my_interface_elements[e];

              for (unsigned int n=0; n<elem->n_vertices(); n++)
                if (std::binary_search (common_interface_node_list.begin(),
                                        common_interface_node_list.end(),
                                        elem->node(n)))
                  {
                    n_shared_nodes++;

                    // TBD - how many nodes do we need to share
                    // before we care?  certainly 2, but 1?  not
                    // sure, so let's play it safe...
                    if (n_shared_nodes > 0) break;
                  }

              if (n_shared_nodes) // share at least one node?
                {
                  elem = elem->top_parent();

                  // avoid a lot of duplicated effort -- if we already have elem
                  // in the set its entire family tree is already in the set.
                  if (!elements_to_send.count(elem))
                    {
#ifdef LIBMESH_ENABLE_AMR
                      elem->family_tree(family_tree);
#else
                      family_tree.clear();
                      family_tree.push_back(elem);
#endif
                      for (unsigned int leaf=0; leaf<family_tree.size(); leaf++)
                        {
                          elem = family_tree[leaf];
                          elements_to_send.insert (elem);

                          for (unsigned int n=0; n<elem->n_nodes(); n++)
                            connected_nodes.insert (elem->get_node(n));
                        }
                    }
                }
            }

          // The elements_to_send and connected_nodes sets now contain all
          // the elements and nodes we need to send to this processor.
          // All that remains is to pack up the objects (along with
          // any boundary conditions) and send the messages off.
          {
            libmesh_assert (connected_nodes.empty() || !elements_to_send.empty());
            libmesh_assert (!connected_nodes.empty() || elements_to_send.empty());

            // send the nodes off to the destination processor
            CommWorld.send_packed_range (dest_pid_idx,
                                         &mesh,
                                         connected_nodes.begin(),
                                         connected_nodes.end(),
                                         send_requests[current_request++],
                                         element_neighbors_tag);

            // send the elements off to the destination processor
            CommWorld.send_packed_range (dest_pid_idx,
                                         &mesh,
                                         elements_to_send.begin(),
                                         elements_to_send.end(),
                                         send_requests[current_request++],
                                         element_neighbors_tag);
          }
        }
      //------------------------------------------------------------------
      // second time - reply of nodes
      else if (n_comm_steps[source_pid_idx] == 2)
        {
          n_comm_steps[source_pid_idx]++;

          CommWorld.receive_packed_range (source_pid_idx,
                                          &mesh,
                                          mesh_inserter_iterator<Node>(mesh),
                                          element_neighbors_tag);
        }
      //------------------------------------------------------------------
      // third time - reply of elements
      else if (n_comm_steps[source_pid_idx] == 3)
        {
          n_comm_steps[source_pid_idx]++;

          CommWorld.receive_packed_range (source_pid_idx,
                                          &mesh,
                                          mesh_inserter_iterator<Elem>(mesh),
                                          element_neighbors_tag);
        }
      //------------------------------------------------------------------
      // fourth time - shouldn't happen
      else
        {
          libMesh::err << "ERROR:  unexpected number of replies: "
                        << n_comm_steps[source_pid_idx]
                    << std::endl;
        }
    } // done catching & processing replies associated with tag ~ 100,000pi

  // allow any pending requests to complete
  Parallel::wait (send_requests);

  STOP_LOG("gather_neighboring_elements()","MeshCommunication");
}
#endif // LIBMESH_HAVE_MPI


#ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
// ------------------------------------------------------------
void MeshCommunication::broadcast (MeshBase&) const
{
  // no MPI == one processor, no need for this method...
  return;
}
#else
// ------------------------------------------------------------
void MeshCommunication::broadcast (MeshBase& mesh) const
{
  // Don't need to do anything if there is
  // only one processor.
  if (libMesh::n_processors() == 1)
    return;

  // This function must be run on all processors at once
  parallel_only();

  START_LOG("broadcast()","MeshCommunication");

  // Explicitly clear the mesh on all but processor 0.
  if (libMesh::processor_id() != 0)
    mesh.clear();
  
  // Broadcast nodes
  CommWorld.broadcast_packed_range(&mesh,
                                   mesh.nodes_begin(),
                                   mesh.nodes_end(),
                                   &mesh,
                                   mesh_inserter_iterator<Node>(mesh));

  // Broadcast elements from coarsest to finest, so that child
  // elements will see their parents already in place.
  unsigned int n_levels = MeshTools::n_levels(mesh);
  CommWorld.broadcast(n_levels);

  for (unsigned int l=0; l != n_levels; ++l)
    CommWorld.broadcast_packed_range(&mesh,
                                     mesh.level_elements_begin(l),
                                     mesh.level_elements_end(l),
                                     &mesh,
                                     mesh_inserter_iterator<Elem>(mesh));

  // Make sure mesh dimension is consistent
  unsigned int mesh_dimension = mesh.mesh_dimension();
  CommWorld.broadcast(mesh_dimension);
  mesh.set_mesh_dimension(mesh_dimension);

  libmesh_assert (CommWorld.verify(mesh.n_elem()));
  libmesh_assert (CommWorld.verify(mesh.n_nodes()));

#ifdef DEBUG
  MeshTools::libmesh_assert_valid_procids<Elem>(mesh);
  MeshTools::libmesh_assert_valid_procids<Node>(mesh);
#endif

  STOP_LOG("broadcast()","MeshCommunication");
}
#endif // LIBMESH_HAVE_MPI



#ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
// ------------------------------------------------------------
void MeshCommunication::allgather (ParallelMesh&) const
{
  // no MPI == one processor, no need for this method...
  return;
}
#else
// ------------------------------------------------------------
void MeshCommunication::allgather (ParallelMesh& mesh) const
{
  // The mesh should know it's about to be serialized
  libmesh_assert (mesh.is_serial());

  // Check for quick return
  if (libMesh::n_processors() == 1)
    return;

  // This function must be run on all processors at once
  parallel_only();

  START_LOG("allgather()","MeshCommunication");

  CommWorld.allgather_packed_range (&mesh,
                                    mesh.nodes_begin(),
                                    mesh.nodes_end(),
                                    mesh_inserter_iterator<Node>(mesh));

  // Gather elements from coarsest to finest, so that child
  // elements will see their parents already in place.
  unsigned int n_levels = MeshTools::n_levels(mesh);
  CommWorld.broadcast(n_levels);

  for (unsigned int l=0; l != n_levels; ++l)
    CommWorld.allgather_packed_range (&mesh,
                                      mesh.level_elements_begin(l),
                                      mesh.level_elements_end(l),
                                      mesh_inserter_iterator<Elem>(mesh));

  libmesh_assert (CommWorld.verify(mesh.n_elem()));
  libmesh_assert (CommWorld.verify(mesh.n_nodes()));

  // Inform new elements of their neighbors,
  // while resetting all remote_elem links
  mesh.find_neighbors(true);

  // All done!
  STOP_LOG("allgather()","MeshCommunication");
}
#endif // LIBMESH_HAVE_MPI



// Functor for make_elems_parallel_consistent and
// make_node_ids_parallel_consistent
namespace {

struct SyncIds
{
typedef dof_id_type datum;
typedef void (MeshBase::*renumber_obj)(dof_id_type, dof_id_type);

SyncIds(MeshBase &_mesh, renumber_obj _renumberer) :
  mesh(_mesh),
  renumber(_renumberer) {}

MeshBase &mesh;
renumber_obj renumber;
// renumber_obj &renumber;

// Find the id of each requested DofObject -
// sync_dofobject_data_by_xyz already did the work for us
void gather_data (const std::vector<dof_id_type>& ids,
                  std::vector<datum>& ids_out)
{
  ids_out = ids;
}

void act_on_data (const std::vector<dof_id_type>& old_ids,
                  std::vector<datum>& new_ids)
{
  for (unsigned int i=0; i != old_ids.size(); ++i)
    if (old_ids[i] != new_ids[i])
      (mesh.*renumber)(old_ids[i], new_ids[i]);
}
};
}



// ------------------------------------------------------------
void MeshCommunication::make_node_ids_parallel_consistent
  (MeshBase &mesh,
   LocationMap<Node> &loc_map)
{
  // This function must be run on all processors at once
  parallel_only();

  START_LOG ("make_node_ids_parallel_consistent()", "MeshCommunication");

  SyncIds syncids(mesh, &MeshBase::renumber_node);
  Parallel::sync_dofobject_data_by_xyz
    (mesh.nodes_begin(), mesh.nodes_end(),
     loc_map, syncids);

  STOP_LOG ("make_node_ids_parallel_consistent()", "MeshCommunication");
}



// ------------------------------------------------------------
void MeshCommunication::make_elems_parallel_consistent(MeshBase &mesh)
{
  // This function must be run on all processors at once
  parallel_only();

  START_LOG ("make_elems_parallel_consistent()", "MeshCommunication");

  SyncIds syncids(mesh, &MeshBase::renumber_elem);
  Parallel::sync_element_data_by_parent_id
    (mesh, mesh.active_elements_begin(),
     mesh.active_elements_end(), syncids);

  STOP_LOG ("make_elems_parallel_consistent()", "MeshCommunication");
}



// Functors for make_node_proc_ids_parallel_consistent
namespace {

struct SyncProcIds
{
typedef processor_id_type datum;

SyncProcIds(MeshBase &_mesh) : mesh(_mesh) {}

MeshBase &mesh;

// ------------------------------------------------------------
void gather_data (const std::vector<dof_id_type>& ids,
                  std::vector<datum>& data)
{
  // Find the processor id of each requested node
  data.resize(ids.size());

  for (std::size_t i=0; i != ids.size(); ++i)
    {
      // Look for this point in the mesh
      // We'd better find every node we're asked for
      Node *node = mesh.node_ptr(ids[i]);

      // Return the node's correct processor id,
      data[i] = node->processor_id();
    }
}

// ------------------------------------------------------------
void act_on_data (const std::vector<dof_id_type>& ids,
                  std::vector<datum> proc_ids)
{
  // Set the ghost node processor ids we've now been informed of
  for (std::size_t i=0; i != ids.size(); ++i)
    {
      Node *node = mesh.node_ptr(ids[i]);
      node->processor_id() = proc_ids[i];
    }
}
};
}



// ------------------------------------------------------------
void MeshCommunication::make_node_proc_ids_parallel_consistent
  (MeshBase& mesh,
   LocationMap<Node>& loc_map)
{
  START_LOG ("make_node_proc_ids_parallel_consistent()", "MeshCommunication");

  // This function must be run on all processors at once
  parallel_only();

  // When this function is called, each section of a parallelized mesh
  // should be in the following state:
  //
  // All nodes should have the exact same physical location on every
  // processor where they exist.
  //
  // Local nodes should have unique authoritative ids,
  // and processor ids consistent with all processors which own
  // an element touching them.
  //
  // Ghost nodes touching local elements should have processor ids
  // consistent with all processors which own an element touching
  // them.

  SyncProcIds sync(mesh);
  Parallel::sync_dofobject_data_by_xyz
    (mesh.nodes_begin(), mesh.nodes_end(), loc_map, sync);

  STOP_LOG ("make_node_proc_ids_parallel_consistent()", "MeshCommunication");
}



// ------------------------------------------------------------
void MeshCommunication::make_nodes_parallel_consistent
  (MeshBase &mesh,
   LocationMap<Node> &loc_map)
{
  // This function must be run on all processors at once
  parallel_only();

  // Create the loc_map if it hasn't been done already
  bool need_map_update = (mesh.nodes_begin() != mesh.nodes_end() &&
                          loc_map.empty());
  CommWorld.max(need_map_update);

  if (need_map_update)
    loc_map.init(mesh);

  // When this function is called, each section of a parallelized mesh
  // should be in the following state:
  //
  // All nodes should have the exact same physical location on every
  // processor where they exist.
  //
  // Local nodes should have unique authoritative ids,
  // and processor ids consistent with all processors which own
  // an element touching them.
  //
  // Ghost nodes touching local elements should have processor ids
  // consistent with all processors which own an element touching
  // them.
  //
  // Ghost nodes should have ids which are either already correct
  // or which are in the "unpartitioned" id space.

  // First, let's sync up processor ids.  Some of these processor ids
  // may be "wrong" from coarsening, but they're right in the sense
  // that they'll tell us who has the authoritative dofobject ids for
  // each node.
  this->make_node_proc_ids_parallel_consistent(mesh, loc_map);

  // Second, sync up dofobject ids.
  this->make_node_ids_parallel_consistent(mesh, loc_map);

  // Finally, correct the processor ids to make DofMap happy
  MeshTools::correct_node_proc_ids(mesh, loc_map);
}



// ------------------------------------------------------------
void MeshCommunication::delete_remote_elements(ParallelMesh& mesh, const std::set<Elem *> & extra_ghost_elem_ids) const
{
  // The mesh should know it's about to be parallelized
  libmesh_assert (!mesh.is_serial());

  START_LOG("delete_remote_elements()", "MeshCommunication");

#ifdef DEBUG
  // We expect maximum ids to be in sync so we can use them to size
  // vectors
  CommWorld.verify(mesh.max_node_id());
  CommWorld.verify(mesh.max_elem_id());
  const dof_id_type par_max_node_id = mesh.parallel_max_node_id();
  const dof_id_type par_max_elem_id = mesh.parallel_max_elem_id();
  libmesh_assert_equal_to (par_max_node_id, mesh.max_node_id());
  libmesh_assert_equal_to (par_max_elem_id, mesh.max_elem_id());
#endif

  // FIXME - should these be "unsorted_set"s?  O(N) is O(N)...
  std::vector<bool> local_nodes(mesh.max_node_id(), false);
  std::vector<bool> semilocal_nodes(mesh.max_node_id(), false);
  std::vector<bool> semilocal_elems(mesh.max_elem_id(), false);

  // We don't want to delete any element that shares a node
  // with or is an ancestor of a local element.
  MeshBase::const_element_iterator l_elem_it = mesh.local_elements_begin(),
                                   l_end     = mesh.local_elements_end();
  for (; l_elem_it != l_end; ++l_elem_it)
    {
      const Elem *elem = *l_elem_it;
      for (unsigned int n=0; n != elem->n_nodes(); ++n)
        {
          dof_id_type nodeid = elem->node(n);
          libmesh_assert_less (nodeid, local_nodes.size());
          local_nodes[nodeid] = true;
        }
      while (elem)
        {
          dof_id_type elemid = elem->id();
          libmesh_assert_less (elemid, semilocal_elems.size());
          semilocal_elems[elemid] = true;

          for (unsigned int n=0; n != elem->n_nodes(); ++n)
            semilocal_nodes[elem->node(n)] = true;

          const Elem *parent = elem->parent();
          // Don't proceed from a boundary mesh to an interior mesh
          if (parent && parent->dim() != elem->dim())
            break;

          elem = parent;
        }
    }

  // We don't want to delete any element that shares a node
  // with or is an ancestor of an unpartitioned element either.
  MeshBase::const_element_iterator
    u_elem_it = mesh.unpartitioned_elements_begin(),
    u_end     = mesh.unpartitioned_elements_end();

  for (; u_elem_it != u_end; ++u_elem_it)
    {
      const Elem *elem = *u_elem_it;
      for (unsigned int n=0; n != elem->n_nodes(); ++n)
        local_nodes[elem->node(n)] = true;
      while (elem)
        {
          semilocal_elems[elem->id()] = true;

          for (unsigned int n=0; n != elem->n_nodes(); ++n)
            semilocal_nodes[elem->node(n)] = true;

          const Elem *parent = elem->parent();
          // Don't proceed from a boundary mesh to an interior mesh
          if (parent && parent->dim() != elem->dim())
            break;

          elem = parent;
        }
    }

  // Flag all the elements that share nodes with
  // local and unpartitioned elements, along with their ancestors
  MeshBase::element_iterator nl_elem_it = mesh.not_local_elements_begin(),
                             nl_end     = mesh.not_local_elements_end();
  for (; nl_elem_it != nl_end; ++nl_elem_it)
    {
      const Elem *elem = *nl_elem_it;
      for (unsigned int n=0; n != elem->n_nodes(); ++n)
        if (local_nodes[elem->node(n)])
          {
            while (elem)
              {
                semilocal_elems[elem->id()] = true;

                for (unsigned int nn=0; nn != elem->n_nodes(); ++nn)
                  semilocal_nodes[elem->node(nn)] = true;

                const Elem *parent = elem->parent();
                // Don't proceed from a boundary mesh to an interior mesh
                if (parent && parent->dim() != elem->dim())
                  break;

                elem = parent;
              }
            break;
          }
    }

  // Don't delete elements that we were explicitly told not to
  for(std::set<Elem *>::iterator it = extra_ghost_elem_ids.begin();
      it != extra_ghost_elem_ids.end();
      ++it)
    {
      const Elem *elem = *it;
      semilocal_elems[elem->id()] = true;
      for (unsigned int n=0; n != elem->n_nodes(); ++n)
        semilocal_nodes[elem->node(n)] = true;
    }

  // Delete all the elements we have no reason to save,
  // starting with the most refined so that the mesh
  // is valid at all intermediate steps
  unsigned int n_levels = MeshTools::n_levels(mesh);

  for (int l = n_levels - 1; l >= 0; --l)
    {
      MeshBase::element_iterator lev_elem_it = mesh.level_elements_begin(l),
                                 lev_end     = mesh.level_elements_end(l);
      for (; lev_elem_it != lev_end; ++lev_elem_it)
        {
          Elem *elem = *lev_elem_it;
          libmesh_assert (elem);
          // Make sure we don't leave any invalid pointers
          if (!semilocal_elems[elem->id()])
            elem->make_links_to_me_remote();

          // Subactive neighbor pointers aren't preservable here
          if (elem->subactive())
            for (unsigned int s=0; s != elem->n_sides(); ++s)
              elem->set_neighbor(s, NULL);

          // delete_elem doesn't currently invalidate element
          // iterators... that had better not change
          if (!semilocal_elems[elem->id()])
            mesh.delete_elem(elem);
        }
    }

  // Delete all the nodes we have no reason to save
  MeshBase::node_iterator node_it  = mesh.nodes_begin(),
                          node_end = mesh.nodes_end();
  for (node_it = mesh.nodes_begin(); node_it != node_end; ++node_it)
    {
      Node *node = *node_it;
      libmesh_assert(node);
      if (!semilocal_nodes[node->id()])
        mesh.delete_node(node);
    }

#ifdef DEBUG
  MeshTools::libmesh_assert_valid_refinement_tree(mesh);
#endif

  STOP_LOG("delete_remote_elements()", "MeshCommunication");
}

} // namespace libMesh

---

Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:47 UTC

Hosted By: