mesh_communication.c

一个用来实现偏微分方程中网格的计算库

  {    
    std::vector<unsigned int>       el_id;
    std::vector<unsigned short int> side_id;
    std::vector<short int>          bc_id;
    
    boundary_info.build_side_list (el_id, side_id, bc_id);

    libmesh_assert (el_id.size() == side_id.size());
    libmesh_assert (el_id.size() == bc_id.size());
    
    const unsigned int n_bcs = el_id.size();

    // reserve an upper bound for the number of BCs
    xfer_elem_bcs.reserve(3*n_bcs);

    // populate the xfer_elem_bcs list with *local* elements only.
    for (unsigned int bc=0; bc<n_bcs; bc++)
      {
	const Elem* elem = mesh.elem(el_id[bc]);
	
	// sanity: be sure that the element returned by mesh.elem() really has id()==el_id[e]
	libmesh_assert(elem != NULL);
	libmesh_assert(elem->id() == el_id[bc]);
	libmesh_assert(elem->level() == 0);
	libmesh_assert(side_id[bc] < elem->n_sides());

	if (elem->processor_id() == libMesh::processor_id())
	  {
	    xfer_elem_bcs.push_back(el_id[bc]);
	    xfer_elem_bcs.push_back(side_id[bc]);
	    xfer_elem_bcs.push_back(bc_id[bc]);
	  }
      }
  
  } // done with element boundary conditions


  // Get the nodal boundary conditions
  {
    std::vector<unsigned int> node_id;
    std::vector<short int>    bc_id;
    
    boundary_info.build_node_list (node_id, bc_id);

    libmesh_assert (node_id.size() == bc_id.size());

    const unsigned int n_bcs = node_id.size();

    // reserve an upper bound for the number of BCs
    xfer_node_bcs.reserve(2*n_bcs);

    // populate the xfer_node_bcs witl *local* nodes only
    for (unsigned int bc=0; bc<n_bcs; bc++)
      {
	const Node* node = mesh.node_ptr(node_id[bc]);
	
	libmesh_assert(node != NULL);
	libmesh_assert(node->id() == node_id[bc]);

	if (node->processor_id() == libMesh::processor_id())
	  {
	    xfer_node_bcs.push_back(node_id[bc]);
	    xfer_node_bcs.push_back(bc_id[bc]);
	  }
      }
  } // done with nodal boundary conditions


  // The xfer arrays now contain all the information for our
  // local bcs, and we are about to get all the information for
  // remote bcs.  Go ahead and clear the current boundary_info
  // information and rebuild it after we get the remote data.
  boundary_info.clear();
  
  // Get the boundary condition information from adacent processors
  Parallel::allgather (xfer_elem_bcs);
  Parallel::allgather (xfer_node_bcs);


  // Insert the elements
  {
    const unsigned int n_bcs = xfer_elem_bcs.size()/3;

    for (unsigned int bc=0, cnt=0; bc<n_bcs; bc++)
      {
	const Elem* elem = mesh.elem(xfer_elem_bcs[cnt++]);
	const unsigned short int side_id = xfer_elem_bcs[cnt++];
	const short int bc_id = xfer_elem_bcs[cnt++];

	boundary_info.add_side (elem, side_id, bc_id);	
      }

    // no need for this any more
    xfer_elem_bcs.resize(0);
  }

  
  // Insert the nodes
  {
    const unsigned int n_bcs = xfer_node_bcs.size()/2;

    for (unsigned int bc=0, cnt=0; bc<n_bcs; bc++)
      {
	const Node* node = mesh.node_ptr (xfer_node_bcs[cnt++]);
	const short int bc_id = xfer_node_bcs[cnt++];

	boundary_info.add_node (node, bc_id);
      }
  }

 
#ifndef NDEBUG
  
  // Make sure all processors agree on the number of boundary ids.
  const unsigned int n_bc_ids = boundary_info.n_boundary_ids();
  unsigned int global_n_bc_ids = n_bc_ids;
  
  Parallel::max (global_n_bc_ids);
  libmesh_assert (n_bc_ids == global_n_bc_ids);
  
#endif
  

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


#endif // HAVE_MPI



// Functor for make_elems_parallel_consistent and
// make_node_ids_parallel_consistent
namespace {

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

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<unsigned int>& ids,
                  std::vector<datum>& ids_out)
{
  ids_out = ids;
}

void act_on_data (const std::vector<unsigned int>& 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 unsigned int datum;

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

MeshBase &mesh;

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

  for (unsigned int i=0; i != ids.size(); ++i)
    {
      // Look for this point in the mesh
      Node *node = mesh.node_ptr(ids[i]);

      // We'd better find every node we're asked for
      libmesh_assert (node);

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

void act_on_data (const std::vector<unsigned int>& ids,
                  std::vector<datum> proc_ids)
{
  // Set the ghost node processor ids we've now been informed of
  for (unsigned int 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());
  Parallel::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
{
  // The mesh should know it's about to be parallelized
  libmesh_assert (!mesh.is_serial());

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

  // 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_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)
        local_nodes[elem->node(n)] = true;
      while (elem)
        {
          semilocal_elems[elem->id()] = true;
          elem = 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;
          elem = 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;
                elem = elem->parent();
              }
            break;
          }
    }

  // 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);
          if (!semilocal_elems[elem->id()])
            {
              // Make sure we don't leave any invalid pointers
              elem->make_links_to_me_remote();

              // delete_elem doesn't currently invalidate element
              // iterators... that had better not change
              mesh.delete_elem(elem);
            }
        }
    }

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

  // Now make sure the containers actually shrink - strip
  // any newly-created NULL voids out of the element array
  mesh.renumber_nodes_and_elements();
}




// // Pack all this information into one communication to avoid two latency hits
// // For each element it is of the form
// // [ level p_level r_flag p_flag etype subdomain_id 
// //   self_ID parent_ID which_child node_0 node_1 ... node_n]
// // We cannot use unsigned int because parent_ID can be negative
// void MeshCommunication::pack_element (std::vector<int> &conn, const Elem* elem) const
// {
//   libmesh_assert (elem != NULL);

// #ifdef ENABLE_AMR
//   conn.push_back (static_cast<int>(elem->level()));
//   conn.push_back (static_cast<int>(elem->p_level()));
//   conn.push_back (static_cast<int>(elem->refinement_flag()));
//   conn.push_back (static_cast<int>(elem->p_refinement_flag()));
// #endif
//   conn.push_back (static_cast<int>(elem->type()));
//   conn.push_back (static_cast<int>(elem->processor_id()));
//   conn.push_back (static_cast<int>(elem->subdomain_id()));
//   conn.push_back (elem->id());
		
// #ifdef ENABLE_AMR
//   // use parent_ID of -1 to indicate a level 0 element
//   if (elem->level() == 0)
//     {
//       conn.push_back(-1);
//       conn.push_back(-1);
//     }
//   else
//     {
//       conn.push_back(elem->parent()->id());
//       conn.push_back(elem->parent()->which_child_am_i(elem));
//     }
// #endif
  
//   for (unsigned int n=0; n<elem->n_nodes(); n++)
//     conn.push_back (elem->node(n));		
// }