$extrastylesheet
libMesh::MeshCommunication Class Reference

#include <mesh_communication.h>

List of all members.

Public Member Functions

 MeshCommunication ()
 ~MeshCommunication ()
void clear ()
void broadcast (MeshBase &) const
void redistribute (ParallelMesh &) const
void gather_neighboring_elements (ParallelMesh &) const
void gather (const processor_id_type root_id, ParallelMesh &) const
void allgather (ParallelMesh &mesh) const
void delete_remote_elements (ParallelMesh &, const std::set< Elem * > &) const
void assign_global_indices (MeshBase &) const
template<typename ForwardIterator >
void find_global_indices (const Parallel::Communicator &communicator, const MeshTools::BoundingBox &, const ForwardIterator &, const ForwardIterator &, std::vector< dof_id_type > &) const
void make_elems_parallel_consistent (MeshBase &)
void make_node_ids_parallel_consistent (MeshBase &)
void make_node_proc_ids_parallel_consistent (MeshBase &)
void make_nodes_parallel_consistent (MeshBase &)

Detailed Description

This is the MeshCommunication class. It handles all the details of communicating mesh information from one processor to another. All parallelization of the Mesh data structures is done via this class.

Author:
Benjamin S. Kirk, 2003

Definition at line 52 of file mesh_communication.h.


Constructor & Destructor Documentation

Constructor.

Definition at line 59 of file mesh_communication.h.

{}

Destructor.

Definition at line 64 of file mesh_communication.h.

{}

Member Function Documentation

void libMesh::MeshCommunication::allgather ( ParallelMesh mesh) const [inline]

This method takes an input ParallelMesh which may be distributed among all the processors. Each processor then sends its local nodes and elements to the other processors. The end result is that a previously distributed ParallelMesh will be serialized on each processor. Since this method is collective it must be called by all processors.

Definition at line 127 of file mesh_communication.h.

References gather(), and libMesh::DofObject::invalid_processor_id.

Referenced by libMesh::ParallelMesh::allgather().

This method assigns globally unique, partition-agnostic indices to the nodes and elements in the mesh. The approach is to compute the Hilbert space-filling curve key and use its value to assign an index in [0,N_global). Since the Hilbert key is unique for each spatial location, two objects occupying the same location will be assigned the same global id. Thus, this method can also be useful for identifying duplicate nodes which may occur during parallel refinement.

Definition at line 183 of file mesh_communication_global_indices.C.

References libMesh::Parallel::Sort< KeyType, IdxType >::bin(), libMesh::MeshTools::bounding_box(), libMesh::Elem::centroid(), libMesh::ParallelObject::comm(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), end, libMesh::err, libMesh::MeshTools::Generation::Private::idx(), libMesh::libmesh_assert(), libMesh::MeshBase::local_elements_begin(), libMesh::MeshBase::local_elements_end(), libMesh::MeshBase::local_nodes_begin(), libMesh::MeshBase::local_nodes_end(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::Threads::parallel_for(), libMesh::DofObject::set_id(), libMesh::Parallel::Sort< KeyType, IdxType >::sort(), and libMesh::START_LOG().

Referenced by libMesh::MeshTools::Private::globally_renumber_nodes_and_elements().

{
  START_LOG ("assign_global_indices()", "MeshCommunication");

  // This method determines partition-agnostic global indices
  // for nodes and elements.

  // Algorithm:
  // (1) compute the Hilbert key for each local node/element
  // (2) perform a parallel sort of the Hilbert key
  // (3) get the min/max value on each processor
  // (4) determine the position in the global ranking for
  //     each local object

  const Parallel::Communicator &communicator (mesh.comm());

  // Global bounding box
  MeshTools::BoundingBox bbox =
    MeshTools::bounding_box (mesh);

  //-------------------------------------------------------------
  // (1) compute Hilbert keys
  std::vector<Hilbert::HilbertIndices>
    node_keys, elem_keys;

  {
    // Nodes first
    {
      ConstNodeRange nr (mesh.local_nodes_begin(),
                         mesh.local_nodes_end());
      node_keys.resize (nr.size());
      Threads::parallel_for (nr, ComputeHilbertKeys (bbox, node_keys));

#if 0
      // It's O(N^2) to check that these keys don't duplicate before the
      // sort...
      MeshBase::const_node_iterator nodei = mesh.local_nodes_begin();
      for (std::size_t i = 0; i != node_keys.size(); ++i, ++nodei)
        {
          MeshBase::const_node_iterator nodej = mesh.local_nodes_begin();
          for (std::size_t j = 0; j != i; ++j, ++nodej)
            {
              if (node_keys[i] == node_keys[j])
                {
                  CFixBitVec icoords[3], jcoords[3];
                  get_hilbert_coords(**nodej, bbox, jcoords);
                  libMesh::err <<
                    "node " << (*nodej)->id() << ", " <<
                    *(Point*)(*nodej) << " has HilbertIndices " <<
                    node_keys[j] << std::endl;
                  get_hilbert_coords(**nodei, bbox, icoords);
                  libMesh::err <<
                    "node " << (*nodei)->id() << ", " <<
                    *(Point*)(*nodei) << " has HilbertIndices " <<
                    node_keys[i] << std::endl;
                  libmesh_error_msg("Error: nodes with duplicate Hilbert keys!");
                }
            }
        }
#endif
    }

    // Elements next
    {
      ConstElemRange er (mesh.local_elements_begin(),
                         mesh.local_elements_end());
      elem_keys.resize (er.size());
      Threads::parallel_for (er, ComputeHilbertKeys (bbox, elem_keys));

#if 0
      // For elements, the keys can be (and in the case of TRI, are
      // expected to be) duplicates, but only if the elements are at
      // different levels
      MeshBase::const_element_iterator elemi = mesh.local_elements_begin();
      for (std::size_t i = 0; i != elem_keys.size(); ++i, ++elemi)
        {
          MeshBase::const_element_iterator elemj = mesh.local_elements_begin();
          for (std::size_t j = 0; j != i; ++j, ++elemj)
            {
              if ((elem_keys[i] == elem_keys[j]) &&
                  ((*elemi)->level() == (*elemj)->level()))
                {
                  libMesh::err <<
                    "level " << (*elemj)->level() << " elem\n" <<
                    (**elemj) << " centroid " <<
                    (*elemj)->centroid() << " has HilbertIndices " <<
                    elem_keys[j] << " or " <<
                    get_hilbert_index((*elemj)->centroid(), bbox) <<
                    std::endl;
                  libMesh::err <<
                    "level " << (*elemi)->level() << " elem\n" <<
                    (**elemi) << " centroid " <<
                    (*elemi)->centroid() << " has HilbertIndices " <<
                    elem_keys[i] << " or " <<
                    get_hilbert_index((*elemi)->centroid(), bbox) <<
                    std::endl;
                  libmesh_error_msg("Error: level " << (*elemi)->level() << " elements with duplicate Hilbert keys!");
                }
            }
        }
#endif
    }
  } // done computing Hilbert keys



  //-------------------------------------------------------------
  // (2) parallel sort the Hilbert keys
  Parallel::Sort<Hilbert::HilbertIndices> node_sorter (communicator,
                                                       node_keys);
  node_sorter.sort(); /* done with node_keys */ //node_keys.clear();

  const std::vector<Hilbert::HilbertIndices> &my_node_bin =
    node_sorter.bin();

  Parallel::Sort<Hilbert::HilbertIndices> elem_sorter (communicator,
                                                       elem_keys);
  elem_sorter.sort(); /* done with elem_keys */ //elem_keys.clear();

  const std::vector<Hilbert::HilbertIndices> &my_elem_bin =
    elem_sorter.bin();



  //-------------------------------------------------------------
  // (3) get the max value on each processor
  std::vector<Hilbert::HilbertIndices>
    node_upper_bounds(communicator.size()),
    elem_upper_bounds(communicator.size());

  { // limit scope of temporaries
    std::vector<Hilbert::HilbertIndices> recvbuf(2*communicator.size());
    std::vector<unsigned short int> /* do not use a vector of bools here since it is not always so! */
      empty_nodes (communicator.size()),
      empty_elem  (communicator.size());
    std::vector<Hilbert::HilbertIndices> my_max(2);

    communicator.allgather (static_cast<unsigned short int>(my_node_bin.empty()), empty_nodes);
    communicator.allgather (static_cast<unsigned short int>(my_elem_bin.empty()),  empty_elem);

    if (!my_node_bin.empty()) my_max[0] = my_node_bin.back();
    if (!my_elem_bin.empty()) my_max[1] = my_elem_bin.back();

    communicator.allgather (my_max, /* identical_buffer_sizes = */ true);

    // Be cereful here.  The *_upper_bounds will be used to find the processor
    // a given object belongs to.  So, if a processor contains no objects (possible!)
    // then copy the bound from the lower processor id.
    for (processor_id_type p=0; p<communicator.size(); p++)
      {
        node_upper_bounds[p] = my_max[2*p+0];
        elem_upper_bounds[p] = my_max[2*p+1];

        if (p > 0) // default hilbert index value is the OK upper bound for processor 0.
          {
            if (empty_nodes[p]) node_upper_bounds[p] = node_upper_bounds[p-1];
            if (empty_elem[p])  elem_upper_bounds[p] = elem_upper_bounds[p-1];
          }
      }
  }



  //-------------------------------------------------------------
  // (4) determine the position in the global ranking for
  //     each local object
  {
    //----------------------------------------------
    // Nodes first -- all nodes, not just local ones
    {
      // Request sets to send to each processor
      std::vector<std::vector<Hilbert::HilbertIndices> >
        requested_ids (communicator.size());
      // Results to gather from each processor
      std::vector<std::vector<dof_id_type> >
        filled_request (communicator.size());

      {
        MeshBase::const_node_iterator       it  = mesh.nodes_begin();
        const MeshBase::const_node_iterator end = mesh.nodes_end();

        // build up list of requests
        for (; it != end; ++it)
          {
            const Node* node = (*it);
            libmesh_assert(node);
            const Hilbert::HilbertIndices hi =
              get_hilbert_index (*node, bbox);
            const processor_id_type pid =
              cast_int<processor_id_type>
              (std::distance (node_upper_bounds.begin(),
                              std::lower_bound(node_upper_bounds.begin(),
                                               node_upper_bounds.end(),
                                               hi)));

            libmesh_assert_less (pid, communicator.size());

            requested_ids[pid].push_back(hi);
          }
      }

      // The number of objects in my_node_bin on each processor
      std::vector<dof_id_type> node_bin_sizes(communicator.size());
      communicator.allgather (static_cast<dof_id_type>(my_node_bin.size()), node_bin_sizes);

      // The offset of my first global index
      dof_id_type my_offset = 0;
      for (processor_id_type pid=0; pid<communicator.rank(); pid++)
        my_offset += node_bin_sizes[pid];

      // start with pid=0, so that we will trade with ourself
      for (processor_id_type pid=0; pid<communicator.size(); pid++)
        {
          // Trade my requests with processor procup and procdown
          const processor_id_type procup = cast_int<processor_id_type>
            ((communicator.rank() + pid) % communicator.size());
          const processor_id_type procdown = cast_int<processor_id_type>
            ((communicator.size() + communicator.rank() - pid) %
             communicator.size());

          std::vector<Hilbert::HilbertIndices> request_to_fill;
          communicator.send_receive(procup, requested_ids[procup],
                                    procdown, request_to_fill);

          // Fill the requests
          std::vector<dof_id_type> global_ids;  global_ids.reserve(request_to_fill.size());
          for (std::size_t idx=0; idx<request_to_fill.size(); idx++)
            {
              const Hilbert::HilbertIndices &hi = request_to_fill[idx];
              libmesh_assert_less_equal (hi, node_upper_bounds[communicator.rank()]);

              // find the requested index in my node bin
              std::vector<Hilbert::HilbertIndices>::const_iterator pos =
                std::lower_bound (my_node_bin.begin(), my_node_bin.end(), hi);
              libmesh_assert (pos != my_node_bin.end());
              libmesh_assert_equal_to (*pos, hi);

              // Finally, assign the global index based off the position of the index
              // in my array, properly offset.
              global_ids.push_back(cast_int<dof_id_type>(std::distance(my_node_bin.begin(), pos) + my_offset));
            }

          // and trade back
          communicator.send_receive (procdown, global_ids,
                                     procup,   filled_request[procup]);
        }

      // We now have all the filled requests, so we can loop through our
      // nodes once and assign the global index to each one.
      {
        std::vector<std::vector<dof_id_type>::const_iterator>
          next_obj_on_proc; next_obj_on_proc.reserve(communicator.size());
        for (processor_id_type pid=0; pid<communicator.size(); pid++)
          next_obj_on_proc.push_back(filled_request[pid].begin());

        {
          MeshBase::node_iterator       it  = mesh.nodes_begin();
          const MeshBase::node_iterator end = mesh.nodes_end();

          for (; it != end; ++it)
            {
              Node* node = (*it);
              libmesh_assert(node);
              const Hilbert::HilbertIndices hi =
                get_hilbert_index (*node, bbox);
              const processor_id_type pid =
                cast_int<processor_id_type>
                (std::distance (node_upper_bounds.begin(),
                                std::lower_bound(node_upper_bounds.begin(),
                                                 node_upper_bounds.end(),
                                                 hi)));

              libmesh_assert_less (pid, communicator.size());
              libmesh_assert (next_obj_on_proc[pid] != filled_request[pid].end());

              const dof_id_type global_index = *next_obj_on_proc[pid];
              libmesh_assert_less (global_index, mesh.n_nodes());
              node->set_id() = global_index;

              ++next_obj_on_proc[pid];
            }
        }
      }
    }

    //---------------------------------------------------
    // elements next -- all elements, not just local ones
    {
      // Request sets to send to each processor
      std::vector<std::vector<Hilbert::HilbertIndices> >
        requested_ids (communicator.size());
      // Results to gather from each processor
      std::vector<std::vector<dof_id_type> >
        filled_request (communicator.size());

      {
        MeshBase::const_element_iterator       it  = mesh.elements_begin();
        const MeshBase::const_element_iterator end = mesh.elements_end();

        for (; it != end; ++it)
          {
            const Elem* elem = (*it);
            libmesh_assert(elem);
            const Hilbert::HilbertIndices hi =
              get_hilbert_index (elem->centroid(), bbox);
            const processor_id_type pid =
              cast_int<processor_id_type>
              (std::distance (elem_upper_bounds.begin(),
                              std::lower_bound(elem_upper_bounds.begin(),
                                               elem_upper_bounds.end(),
                                               hi)));

            libmesh_assert_less (pid, communicator.size());

            requested_ids[pid].push_back(hi);
          }
      }

      // The number of objects in my_elem_bin on each processor
      std::vector<dof_id_type> elem_bin_sizes(communicator.size());
      communicator.allgather (static_cast<dof_id_type>(my_elem_bin.size()), elem_bin_sizes);

      // The offset of my first global index
      dof_id_type my_offset = 0;
      for (processor_id_type pid=0; pid<communicator.rank(); pid++)
        my_offset += elem_bin_sizes[pid];

      // start with pid=0, so that we will trade with ourself
      for (processor_id_type pid=0; pid<communicator.size(); pid++)
        {
          // Trade my requests with processor procup and procdown
          const processor_id_type procup = cast_int<processor_id_type>
            ((communicator.rank() + pid) % communicator.size());
          const processor_id_type procdown = cast_int<processor_id_type>
            ((communicator.size() + communicator.rank() - pid) %
             communicator.size());

          std::vector<Hilbert::HilbertIndices> request_to_fill;
          communicator.send_receive(procup, requested_ids[procup],
                                    procdown, request_to_fill);

          // Fill the requests
          std::vector<dof_id_type> global_ids;  global_ids.reserve(request_to_fill.size());
          for (std::size_t idx=0; idx<request_to_fill.size(); idx++)
            {
              const Hilbert::HilbertIndices &hi = request_to_fill[idx];
              libmesh_assert_less_equal (hi, elem_upper_bounds[communicator.rank()]);

              // find the requested index in my elem bin
              std::vector<Hilbert::HilbertIndices>::const_iterator pos =
                std::lower_bound (my_elem_bin.begin(), my_elem_bin.end(), hi);
              libmesh_assert (pos != my_elem_bin.end());
              libmesh_assert_equal_to (*pos, hi);

              // Finally, assign the global index based off the position of the index
              // in my array, properly offset.
              global_ids.push_back (cast_int<dof_id_type>(std::distance(my_elem_bin.begin(), pos) + my_offset));
            }

          // and trade back
          communicator.send_receive (procdown, global_ids,
                                     procup,   filled_request[procup]);
        }

      // We now have all the filled requests, so we can loop through our
      // elements once and assign the global index to each one.
      {
        std::vector<std::vector<dof_id_type>::const_iterator>
          next_obj_on_proc; next_obj_on_proc.reserve(communicator.size());
        for (processor_id_type pid=0; pid<communicator.size(); pid++)
          next_obj_on_proc.push_back(filled_request[pid].begin());

        {
          MeshBase::element_iterator       it  = mesh.elements_begin();
          const MeshBase::element_iterator end = mesh.elements_end();

          for (; it != end; ++it)
            {
              Elem* elem = (*it);
              libmesh_assert(elem);
              const Hilbert::HilbertIndices hi =
                get_hilbert_index (elem->centroid(), bbox);
              const processor_id_type pid =
                cast_int<processor_id_type>
                (std::distance (elem_upper_bounds.begin(),
                                std::lower_bound(elem_upper_bounds.begin(),
                                                 elem_upper_bounds.end(),
                                                 hi)));

              libmesh_assert_less (pid, communicator.size());
              libmesh_assert (next_obj_on_proc[pid] != filled_request[pid].end());

              const dof_id_type global_index = *next_obj_on_proc[pid];
              libmesh_assert_less (global_index, mesh.n_elem());
              elem->set_id() = global_index;

              ++next_obj_on_proc[pid];
            }
        }
      }
    }
  }

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

Finds all the processors that may contain elements that neighbor my elements. This list is guaranteed to include all processors that border any of my elements, but may include additional ones as well. This method computes bounding boxes for the elements on each processor and checks for overlaps. This method takes a mesh (which is assumed to reside on processor 0) and broadcasts it to all the other processors. It also broadcasts any boundary information the mesh has associated with it.

Definition at line 665 of file mesh_communication.C.

Referenced by libMesh::NameBasedIO::read(), and libMesh::CheckpointIO::read().

{
  // no MPI == one processor, no need for this method...
  return;
}

Clears all data structures and returns to a pristine state.

Definition at line 79 of file mesh_communication.C.

{
  //  _neighboring_processors.clear();
}
void libMesh::MeshCommunication::delete_remote_elements ( ParallelMesh mesh,
const std::set< Elem * > &  extra_ghost_elem_ids 
) const

This method takes an input ParallelMesh which may be distributed among all the processors. Each processor deletes all elements which are neither local elements nor "ghost" elements which touch local elements, and deletes all nodes which are not contained in local or ghost elements. The end result is that a previously serial ParallelMesh will be distributed between processors. Since this method is collective it must be called by all processors.

The std::set is a list of extra elements that you _don't_ want to delete. These will be left on the current processor along with local elements and ghosted neighbors.

Definition at line 998 of file mesh_communication.C.

References libMesh::ParallelObject::comm(), libMesh::ParallelMesh::delete_elem(), libMesh::ParallelMesh::delete_node(), libMesh::Elem::dim(), libMesh::DofObject::id(), libMesh::ParallelMesh::is_serial(), libMesh::ParallelMesh::level_elements_begin(), libMesh::ParallelMesh::level_elements_end(), libMesh::libmesh_assert(), libMesh::MeshTools::libmesh_assert_valid_refinement_tree(), libMesh::ParallelMesh::local_elements_begin(), libMesh::ParallelMesh::local_elements_end(), libMesh::Elem::make_links_to_me_remote(), libMesh::ParallelMesh::max_elem_id(), libMesh::ParallelMesh::max_node_id(), libMesh::MeshTools::n_levels(), libMesh::Elem::n_nodes(), libMesh::Elem::n_sides(), libMesh::Elem::node(), libMesh::ParallelMesh::nodes_begin(), libMesh::ParallelMesh::nodes_end(), libMesh::ParallelMesh::not_local_elements_begin(), libMesh::ParallelMesh::not_local_elements_end(), libMesh::ParallelMesh::parallel_max_elem_id(), libMesh::ParallelMesh::parallel_max_node_id(), libMesh::Elem::parent(), libMesh::Elem::set_neighbor(), libMesh::START_LOG(), libMesh::Elem::subactive(), libMesh::ParallelMesh::unpartitioned_elements_begin(), libMesh::ParallelMesh::unpartitioned_elements_end(), and libMesh::Parallel::Communicator::verify().

Referenced by libMesh::ParallelMesh::delete_remote_elements().

{
  // 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
  mesh.comm().verify(mesh.max_node_id());
  mesh.comm().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");
}
template<typename ForwardIterator >
void libMesh::MeshCommunication::find_global_indices ( const Parallel::Communicator communicator,
const MeshTools::BoundingBox bbox,
const ForwardIterator &  begin,
const ForwardIterator &  end,
std::vector< dof_id_type > &  index_map 
) const

This method determines a globally unique, partition-agnostic index for each object in the input range.

Definition at line 598 of file mesh_communication_global_indices.C.

References libMesh::Parallel::Communicator::allgather(), libMesh::Parallel::Sort< KeyType, IdxType >::bin(), end, libMesh::MeshTools::Generation::Private::idx(), libMesh::DofObject::invalid_processor_id, libMesh::libmesh_assert(), std::max(), libMesh::Parallel::Communicator::rank(), libMesh::Real, libMesh::Parallel::Communicator::send_receive(), libMesh::Parallel::Communicator::size(), libMesh::Parallel::Sort< KeyType, IdxType >::sort(), and libMesh::START_LOG().

Referenced by libMesh::MetisPartitioner::_do_partition(), libMesh::ParmetisPartitioner::initialize(), and libMesh::Partitioner::partition_unpartitioned_elements().

{
  START_LOG ("find_global_indices()", "MeshCommunication");

  // This method determines partition-agnostic global indices
  // for nodes and elements.

  // Algorithm:
  // (1) compute the Hilbert key for each local node/element
  // (2) perform a parallel sort of the Hilbert key
  // (3) get the min/max value on each processor
  // (4) determine the position in the global ranking for
  //     each local object
  index_map.clear();
  index_map.reserve(std::distance (begin, end));

  //-------------------------------------------------------------
  // (1) compute Hilbert keys
  // These aren't trivial to compute, and we will need them again.
  // But the binsort will sort the input vector, trashing the order
  // that we'd like to rely on.  So, two vectors...
  std::vector<Hilbert::HilbertIndices>
    sorted_hilbert_keys,
    hilbert_keys;
  sorted_hilbert_keys.reserve(index_map.capacity());
  hilbert_keys.reserve(index_map.capacity());
  {
    START_LOG("compute_hilbert_indices()", "MeshCommunication");
    for (ForwardIterator it=begin; it!=end; ++it)
      {
        const Hilbert::HilbertIndices hi(get_hilbert_index (*it, bbox));
        hilbert_keys.push_back(hi);

        if ((*it)->processor_id() == communicator.rank())
          sorted_hilbert_keys.push_back(hi);

        // someone needs to take care of unpartitioned objects!
        if ((communicator.rank() == 0) &&
            ((*it)->processor_id() == DofObject::invalid_processor_id))
          sorted_hilbert_keys.push_back(hi);
      }
    STOP_LOG("compute_hilbert_indices()", "MeshCommunication");
  }

  //-------------------------------------------------------------
  // (2) parallel sort the Hilbert keys
  START_LOG ("parallel_sort()", "MeshCommunication");
  Parallel::Sort<Hilbert::HilbertIndices> sorter (communicator,
                                                  sorted_hilbert_keys);
  sorter.sort();
  STOP_LOG ("parallel_sort()", "MeshCommunication");
  const std::vector<Hilbert::HilbertIndices> &my_bin = sorter.bin();

  // The number of objects in my_bin on each processor
  std::vector<unsigned int> bin_sizes(communicator.size());
  communicator.allgather (static_cast<unsigned int>(my_bin.size()), bin_sizes);

  // The offset of my first global index
  unsigned int my_offset = 0;
  for (unsigned int pid=0; pid<communicator.rank(); pid++)
    my_offset += bin_sizes[pid];

  //-------------------------------------------------------------
  // (3) get the max value on each processor
  std::vector<Hilbert::HilbertIndices>
    upper_bounds(1);

  if (!my_bin.empty())
    upper_bounds[0] = my_bin.back();

  communicator.allgather (upper_bounds, /* identical_buffer_sizes = */ true);

  // Be cereful here.  The *_upper_bounds will be used to find the processor
  // a given object belongs to.  So, if a processor contains no objects (possible!)
  // then copy the bound from the lower processor id.
  for (unsigned int p=1; p<communicator.size(); p++)
    if (!bin_sizes[p]) upper_bounds[p] = upper_bounds[p-1];


  //-------------------------------------------------------------
  // (4) determine the position in the global ranking for
  //     each local object
  {
    //----------------------------------------------
    // all objects, not just local ones

    // Request sets to send to each processor
    std::vector<std::vector<Hilbert::HilbertIndices> >
      requested_ids (communicator.size());
    // Results to gather from each processor
    std::vector<std::vector<dof_id_type> >
      filled_request (communicator.size());

    // build up list of requests
    std::vector<Hilbert::HilbertIndices>::const_iterator hi =
      hilbert_keys.begin();

    for (ForwardIterator it = begin; it != end; ++it)
      {
        libmesh_assert (hi != hilbert_keys.end());
        const processor_id_type pid =
          cast_int<processor_id_type>
          (std::distance (upper_bounds.begin(),
                          std::lower_bound(upper_bounds.begin(),
                                           upper_bounds.end(),
                                           *hi)));

        libmesh_assert_less (pid, communicator.size());

        requested_ids[pid].push_back(*hi);

        ++hi;
        // go ahead and put pid in index_map, that way we
        // don't have to repeat the std::lower_bound()
        index_map.push_back(pid);
      }

    // start with pid=0, so that we will trade with ourself
    std::vector<Hilbert::HilbertIndices> request_to_fill;
    std::vector<dof_id_type> global_ids;
    for (processor_id_type pid=0; pid<communicator.size(); pid++)
      {
        // Trade my requests with processor procup and procdown
        const processor_id_type procup = cast_int<processor_id_type>
          ((communicator.rank() + pid) % communicator.size());
        const processor_id_type procdown = cast_int<processor_id_type>
          ((communicator.size() + communicator.rank() - pid) %
           communicator.size());

        communicator.send_receive(procup, requested_ids[procup],
                                  procdown, request_to_fill);

        // Fill the requests
        global_ids.clear();  global_ids.reserve(request_to_fill.size());
        for (unsigned int idx=0; idx<request_to_fill.size(); idx++)
          {
            const Hilbert::HilbertIndices &hilbert_indices = request_to_fill[idx];
            libmesh_assert_less_equal (hilbert_indices, upper_bounds[communicator.rank()]);

            // find the requested index in my node bin
            std::vector<Hilbert::HilbertIndices>::const_iterator pos =
              std::lower_bound (my_bin.begin(), my_bin.end(), hilbert_indices);
            libmesh_assert (pos != my_bin.end());
#ifdef DEBUG
            // If we could not find the requested Hilbert index in
            // my_bin, something went terribly wrong, possibly the
            // Mesh was displaced differently on different processors,
            // and therefore the Hilbert indices don't agree.
            if (*pos != hilbert_indices)
              {
                // The input will be hilbert_indices.  We convert it
                // to BitVecType using the operator= provided by the
                // BitVecType class. BitVecType is a CBigBitVec!
                Hilbert::BitVecType input;
                input = hilbert_indices;

                // Get output in a vector of CBigBitVec
                std::vector<CBigBitVec> output(3);

                // Call the indexToCoords function
                Hilbert::indexToCoords(&output[0], 8*sizeof(Hilbert::inttype), 3, input);

                // The entries in the output racks are integers in the
                // range [0, Hilbert::inttype::max] which can be
                // converted to floating point values in [0,1] and
                // finally to actual values using the bounding box.
                Real max_int_as_real = static_cast<Real>(std::numeric_limits<Hilbert::inttype>::max());

                // Get the points in [0,1]^3.  The zeroth rack of each entry in
                // 'output' maps to the normalized x, y, and z locations,
                // respectively.
                Point p_hat(static_cast<Real>(output[0].racks()[0]) / max_int_as_real,
                            static_cast<Real>(output[1].racks()[0]) / max_int_as_real,
                            static_cast<Real>(output[2].racks()[0]) / max_int_as_real);

                // Convert the points from [0,1]^3 to their actual (x,y,z) locations
                Real
                  xmin = bbox.first(0),
                  xmax = bbox.second(0),
                  ymin = bbox.first(1),
                  ymax = bbox.second(1),
                  zmin = bbox.first(2),
                  zmax = bbox.second(2);

                // Convert the points from [0,1]^3 to their actual (x,y,z) locations
                Point p(xmin + (xmax-xmin)*p_hat(0),
                        ymin + (ymax-ymin)*p_hat(1),
                        zmin + (zmax-zmin)*p_hat(2));

                libmesh_error_msg("Could not find hilbert indices: "
                                  << hilbert_indices
                                  << " corresponding to point " << p);
              }
#endif

            // Finally, assign the global index based off the position of the index
            // in my array, properly offset.
            global_ids.push_back (cast_int<dof_id_type>(std::distance(my_bin.begin(), pos) + my_offset));
          }

        // and trade back
        communicator.send_receive (procdown, global_ids,
                                   procup,   filled_request[procup]);
      }

    // We now have all the filled requests, so we can loop through our
    // nodes once and assign the global index to each one.
    {
      std::vector<std::vector<dof_id_type>::const_iterator>
        next_obj_on_proc; next_obj_on_proc.reserve(communicator.size());
      for (unsigned int pid=0; pid<communicator.size(); pid++)
        next_obj_on_proc.push_back(filled_request[pid].begin());

      unsigned int cnt=0;
      for (ForwardIterator it = begin; it != end; ++it, cnt++)
        {
          const processor_id_type pid = cast_int<processor_id_type>
            (index_map[cnt]);

          libmesh_assert_less (pid, communicator.size());
          libmesh_assert (next_obj_on_proc[pid] != filled_request[pid].end());

          const dof_id_type global_index = *next_obj_on_proc[pid];
          index_map[cnt] = global_index;

          ++next_obj_on_proc[pid];
        }
    }
  }

  STOP_LOG ("find_global_indices()", "MeshCommunication");
}
void libMesh::MeshCommunication::gather ( const processor_id_type  root_id,
ParallelMesh mesh 
) const

This method takes an input ParallelMesh which may be distributed among all the processors. Each processor then sends its local nodes and elements to processor root_id. The end result is that a previously distributed ParallelMesh will be serialized on processor root_id. Since this method is collective it must be called by all processors. For the special case of root_id equal to DofObject::invalid_processor_id this function performs an allgather.

Definition at line 731 of file mesh_communication.C.

Referenced by allgather().

{
  // no MPI == one processor, no need for this method...
  return;
}

Definition at line 306 of file mesh_communication.C.

Referenced by libMesh::Nemesis_IO::read().

{
  // no MPI == one processor, no need for this method...
  return;
}

Copy ids of ghost elements from their local processors.

Definition at line 864 of file mesh_communication.C.

References libMesh::MeshBase::active_elements_begin(), libMesh::MeshBase::active_elements_end(), libMesh::ParallelObject::comm(), libMesh::START_LOG(), and libMesh::Parallel::sync_element_data_by_parent_id().

Referenced by libMesh::MeshRefinement::_refine_elements().

{
  // This function must be run on all processors at once
  libmesh_parallel_only(mesh.comm());

  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");
}

Assuming all ids on local nodes are globally unique, and assuming all processor ids are parallel consistent, this function makes all other ids parallel consistent.

Definition at line 847 of file mesh_communication.C.

References libMesh::ParallelObject::comm(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), libMesh::START_LOG(), and libMesh::Parallel::sync_node_data_by_element_id().

{
  // This function must be run on all processors at once
  libmesh_parallel_only(mesh.comm());

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

  SyncIds syncids(mesh, &MeshBase::renumber_node);
  Parallel::sync_node_data_by_element_id
    (mesh, mesh.elements_begin(), mesh.elements_end(), syncids);

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

Assuming all processor ids on nodes touching local elements are parallel consistent, this function makes all other processor ids parallel consistent as well.

Definition at line 928 of file mesh_communication.C.

References libMesh::ParallelObject::comm(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), libMesh::START_LOG(), and libMesh::Parallel::sync_node_data_by_element_id().

Referenced by libMesh::MeshTools::correct_node_proc_ids().

{
  START_LOG ("make_node_proc_ids_parallel_consistent()", "MeshCommunication");

  // This function must be run on all processors at once
  libmesh_parallel_only(mesh.comm());

  // 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_node_data_by_element_id
    (mesh, mesh.elements_begin(), mesh.elements_end(), sync);

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

Copy processor_ids and ids on ghost nodes from their local processors. This is an internal function of MeshRefinement which turns out to be useful for other code which wants to add nodes to a distributed mesh.

Definition at line 960 of file mesh_communication.C.

References libMesh::ParallelObject::comm(), and libMesh::MeshTools::correct_node_proc_ids().

Referenced by libMesh::MeshRefinement::_coarsen_elements(), libMesh::MeshRefinement::_refine_elements(), and libMesh::UnstructuredMesh::all_second_order().

{
  // This function must be run on all processors at once
  libmesh_parallel_only(mesh.comm());

  // 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);

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

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

This method takes a parallel distributed mesh and redistributes the elements. Specifically, any elements stored on a given processor are sent to the processor which "owns" them. Similarly, any elements assigned to the current processor but stored on another are received. Once this step is completed any required ghost elements are updated. The final result is that each processor stores only the elements it actually owns and any ghost elements required to satisfy data dependencies. This method can be invoked after a partitioning step to affect the new partitioning.

Definition at line 88 of file mesh_communication.C.

Referenced by libMesh::ParallelMesh::redistribute().

{
  // no MPI == one processor, no redistribution
  return;
}

The documentation for this class was generated from the following files: