parmetis_partitioner.c

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

    // create the unique mapping for all active elements independent of partitioning
    {
      MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::const_element_iterator end = mesh.active_elements_end(); 
	
      // Calling this on all processors a unique range in [0,n_active_elem) is constructed.  
      // Only the indices for the elements we pass in are returned in the array.
      MeshCommunication().find_global_indices (bbox, it, end, 
					       global_index);
      
      for (unsigned int cnt=0; it != end; ++it)
	{
	  const Elem *elem = *it;
	  libmesh_assert (!global_index_map.count(elem->id()));
	  libmesh_assert (cnt < global_index.size());
	  libmesh_assert (global_index[cnt] < n_active_elem);

	  global_index_map[elem->id()]  = global_index[cnt++];
	}
      }
    // really, shouldn't be close!
    libmesh_assert (global_index_map.size() <= n_active_elem);
    libmesh_assert (_global_index_by_pid_map.size() <= n_active_elem);

    // At this point the two maps should be the same size.  If they are not
    // then the number of active elements is not the same as the sum over all 
    // processors of the number of active elements per processor, which means
    // there must be some unpartitioned objects out there.
    if (global_index_map.size() != _global_index_by_pid_map.size())
      {
	std::cerr << "ERROR:  ParmetisPartitioner cannot handle unpartitioned objects!"
		  << std::endl;
	libmesh_error();
      }    
  }

  // Finally, we need to initialize the vertex (partition) weights and the initial subdomain
  // mapping.  The subdomain mapping will be independent of the processor mapping, and is 
  // defined by a simple mapping of the global indices we just found.
  {
    std::vector<unsigned int> subdomain_bounds(libMesh::n_processors());

    const unsigned int first_local_elem = _vtxdist[libMesh::processor_id()];

    for (unsigned int pid=0; pid<libMesh::n_processors(); pid++)
      {
	unsigned int tgt_subdomain_size = 0;
	
	// watch out for the case that n_subdomains < n_processors
	if (pid < static_cast<unsigned int>(_nparts))
	  {	
	    tgt_subdomain_size =
	      n_active_elem/std::min(libMesh::n_processors(),
				     static_cast<unsigned int>(_nparts));
	    
	    if (pid < n_active_elem%_nparts)
	      tgt_subdomain_size++;
	  }
	if (pid == 0)
	  subdomain_bounds[0] = tgt_subdomain_size;
	else
	  subdomain_bounds[pid] = subdomain_bounds[pid-1] + tgt_subdomain_size;
      }

    libmesh_assert (subdomain_bounds.back() == n_active_elem);

    MeshBase::const_element_iterator       elem_it  = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end();
    
    for (; elem_it != elem_end; ++elem_it)
      {
	const Elem *elem = *elem_it;
	
	libmesh_assert (_global_index_by_pid_map.count(elem->id()));
	const unsigned int global_index_by_pid = 
	  _global_index_by_pid_map[elem->id()];	
	libmesh_assert (global_index_by_pid < n_active_elem);
	
	const unsigned int local_index =
	  global_index_by_pid - first_local_elem;
	
	libmesh_assert (local_index < n_active_local_elem);
	libmesh_assert (local_index < _vwgt.size());
	
	// TODO:[BSK] maybe there is a better weight?
	_vwgt[local_index] = elem->n_nodes();
	
	// find the subdomain this element belongs in
	libmesh_assert (global_index_map.count(elem->id()));
	const unsigned int global_index =
	  global_index_map[elem->id()];
	
	libmesh_assert (global_index < subdomain_bounds.back());
	
	const unsigned int subdomain_id =
	  std::distance(subdomain_bounds.begin(),
			std::lower_bound(subdomain_bounds.begin(),
					 subdomain_bounds.end(),
					 global_index));
	libmesh_assert (subdomain_id < static_cast<unsigned int>(_nparts));
	libmesh_assert (local_index < _part.size());
	  
	_part[local_index] = subdomain_id;
      }    
  }
}



void ParmetisPartitioner::build_graph (const MeshBase& mesh)
{
  // build the graph in distributed CSR format.  Note that
  // the edges in the graph will correspond to
  // face neighbors
  const unsigned int n_active_local_elem  = mesh.n_active_local_elem();
    
  std::vector<const Elem*> neighbors_offspring;
  
  std::vector<std::vector<unsigned int> > graph(n_active_local_elem);
  unsigned int graph_size=0;

  const unsigned int first_local_elem = _vtxdist[libMesh::processor_id()];
  
  MeshBase::const_element_iterator       elem_it  = mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end(); 
  
  for (; elem_it != elem_end; ++elem_it)
    {
      const Elem* elem = *elem_it;
      	
      libmesh_assert (_global_index_by_pid_map.count(elem->id()));
      const unsigned int global_index_by_pid = 
	_global_index_by_pid_map[elem->id()];	
	
      const unsigned int local_index =
	global_index_by_pid - first_local_elem;
      libmesh_assert (local_index < n_active_local_elem);

      std::vector<unsigned int> &graph_row = graph[local_index];

      // Loop over the element's neighbors.  An element
      // adjacency corresponds to a face neighbor
      for (unsigned int ms=0; ms<elem->n_neighbors(); ms++)
	{
	  const Elem* neighbor = elem->neighbor(ms);
	  
	  if (neighbor != NULL)
	    {  
	      // If the neighbor is active treat it
	      // as a connection
	      if (neighbor->active())
		{
		  libmesh_assert(_global_index_by_pid_map.count(neighbor->id()));
		  const unsigned int neighbor_global_index_by_pid =
		    _global_index_by_pid_map[neighbor->id()];
		  
		  graph_row.push_back(neighbor_global_index_by_pid);
		  graph_size++;
		}
  
#ifdef ENABLE_AMR
	      
	      // Otherwise we need to find all of the
	      // neighbor's children that are connected to
	      // us and add them
	      else
		{
		  // The side of the neighbor to which
		  // we are connected
		  const unsigned int ns =
		    neighbor->which_neighbor_am_i (elem);
                  libmesh_assert (ns < neighbor->n_neighbors());
		  
		  // Get all the active children (& grandchildren, etc...)
		  // of the neighbor.
		  neighbor->active_family_tree (neighbors_offspring);
		  
		  // Get all the neighbor's children that
		  // live on that side and are thus connected
		  // to us
		  for (unsigned int nc=0; nc<neighbors_offspring.size(); nc++)
		    {
		      const Elem* child =
			neighbors_offspring[nc];
		      
		      // This does not assume a level-1 mesh.
		      // Note that since children have sides numbered
		      // coincident with the parent then this is a sufficient test.
		      if (child->neighbor(ns) == elem)
			{
			  libmesh_assert (child->active());
			  libmesh_assert (_global_index_by_pid_map.count(child->id()));
			  const unsigned int child_global_index_by_pid =
			    _global_index_by_pid_map[child->id()];

			  graph_row.push_back(child_global_index_by_pid);
			  graph_size++;
			}
		    }
		}

#endif /* ifdef ENABLE_AMR */


	    }
	}
    }

  // Reserve space in the adjacency array
  _xadj.clear();
  _xadj.reserve (n_active_local_elem + 1);
  _adjncy.clear();
  _adjncy.reserve (graph_size);

  for (unsigned int r=0; r<graph.size(); r++)
    {
      _xadj.push_back(_adjncy.size());
      std::vector<unsigned int> graph_row; // build this emtpy
      graph_row.swap(graph[r]); // this will deallocate at the end of scope
      _adjncy.insert(_adjncy.end(),
		     graph_row.begin(),
		     graph_row.end());
    }
  
  // The end of the adjacency array for the last elem
  _xadj.push_back(_adjncy.size());
  
  libmesh_assert (_xadj.size() == n_active_local_elem+1);
  libmesh_assert (_adjncy.size() == graph_size);
}



void ParmetisPartitioner::assign_partitioning (MeshBase& mesh)
{
  const unsigned int 
    first_local_elem = _vtxdist[libMesh::processor_id()];
    
  std::vector<std::vector<unsigned int> > 
    requested_ids(libMesh::n_processors()),
    requests_to_fill(libMesh::n_processors());
  
  MeshBase::element_iterator elem_it  = mesh.active_elements_begin();
  MeshBase::element_iterator elem_end = mesh.active_elements_end();
    
  for (; elem_it != elem_end; ++elem_it)
    {	
      Elem *elem = *elem_it;

      // we need to get the index from the owning processor
      // (note we cannot assign it now -- we are iterating
      // over elements again and this will be bad!)
      libmesh_assert (elem->processor_id() < requested_ids.size());
      requested_ids[elem->processor_id()].push_back(elem->id());
    }
  
  // Trade with all processors (including self) to get their indices
  for (unsigned int pid=0; pid<libMesh::n_processors(); pid++)
    {
      // Trade my requests with processor procup and procdown
      const unsigned int procup = (libMesh::processor_id() + pid) %
	                           libMesh::n_processors();
      const unsigned int procdown = (libMesh::n_processors() +
				     libMesh::processor_id() - pid) %
	                             libMesh::n_processors();

      Parallel::send_receive (procup,   requested_ids[procup],
			      procdown, requests_to_fill[procdown]); 

      // we can overwrite these requested ids in-place.
      for (unsigned int i=0; i<requests_to_fill[procdown].size(); i++)
	{
	  const unsigned int requested_elem_index = 
	    requests_to_fill[procdown][i];

	  libmesh_assert(_global_index_by_pid_map.count(requested_elem_index));

	  const unsigned int global_index_by_pid = 
	    _global_index_by_pid_map[requested_elem_index];

	  const unsigned int local_index =
	    global_index_by_pid - first_local_elem;
	    
	  libmesh_assert (local_index < _part.size());
	  libmesh_assert (local_index < mesh.n_active_local_elem());
	  
	  const unsigned int elem_procid =
	    static_cast<unsigned int>(_part[local_index]);

	  libmesh_assert (elem_procid < static_cast<unsigned int>(_nparts));

	  requests_to_fill[procdown][i] = elem_procid;
	}

      // Trade back
      Parallel::send_receive (procdown, requests_to_fill[procdown],
			      procup,   requested_ids[procup]);
    }

   // and finally assign the partitioning.
   // note we are iterating in exactly the same order
   // used to build up the request, so we can expect the
   // required entries to be in the proper sequence.
   elem_it  = mesh.active_elements_begin();
   elem_end = mesh.active_elements_end();      
   
   for (std::vector<unsigned int> counters(libMesh::n_processors(), 0);  
	elem_it != elem_end; ++elem_it)
     {	
       Elem *elem = *elem_it;
       
       const unsigned int current_pid = elem->processor_id();
       
       libmesh_assert (counters[current_pid] < requested_ids[current_pid].size());
       
       const unsigned int elem_procid = 
	 requested_ids[current_pid][counters[current_pid]++];
       
       libmesh_assert (elem_procid < static_cast<unsigned int>(_nparts));
       elem->processor_id() = elem_procid;
     }
}

#endif // #ifdef HAVE_PARMETIS