$extrastylesheet
#include <mesh_communication.h>
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 &) |
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.
Definition at line 52 of file mesh_communication.h.
| libMesh::MeshCommunication::MeshCommunication | ( | ) | [inline] |
| libMesh::MeshCommunication::~MeshCommunication | ( | ) | [inline] |
| 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().
| void libMesh::MeshCommunication::assign_global_indices | ( | MeshBase & | mesh | ) | const |
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");
}
| void libMesh::MeshCommunication::broadcast | ( | MeshBase & | mesh | ) | const |
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;
}
| void libMesh::MeshCommunication::clear | ( | ) |
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");
}
| 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;
}
| void libMesh::MeshCommunication::gather_neighboring_elements | ( | ParallelMesh & | mesh | ) | const |
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);
}
| void libMesh::MeshCommunication::redistribute | ( | ParallelMesh & | mesh | ) | const |
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;
}