metis_partitioner.c

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

// $Id: metis_partitioner.C 2789 2008-04-13 02:24:40Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
  
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



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

// Local Includes -----------------------------------
#include "libmesh_config.h"
#include "mesh_base.h"
#include "metis_partitioner.h"
#include "libmesh_logging.h"
#include "elem.h"
#include "mesh_communication.h"

#ifdef HAVE_METIS
// MIPSPro 7.4.2 gets confused about these nested namespaces
# ifdef __sgi
#  include <cstdarg>
# endif
  namespace Metis {
    extern "C" {
#     include "metis.h"
    }
  }
#else
#  include "sfc_partitioner.h"
#endif



// ------------------------------------------------------------
// MetisPartitioner implementation
void MetisPartitioner::_do_partition (MeshBase& mesh,
				      const unsigned int n_pieces)
{
  libmesh_assert (n_pieces > 0);
  libmesh_assert (mesh.is_serial());

  // Check for an easy return
  if (n_pieces == 1)
    {
      this->single_partition (mesh);
      return;
    }
  
// What to do if the Metis library IS NOT present
#ifndef HAVE_METIS

  here();
  std::cerr << "ERROR: The library has been built without"    << std::endl
	    << "Metis support.  Using a space-filling curve"  << std::endl
	    << "partitioner instead!"                         << std::endl;

  SFCPartitioner sfcp;

  sfcp.partition (mesh, n_pieces);
  
// What to do if the Metis library IS present
#else
  
  START_LOG("partition()", "MetisPartitioner");

  const unsigned int n_active_elem = mesh.n_active_elem();
  
  // build the graph  
  std::vector<int> options(5);
  std::vector<int> vwgt(n_active_elem);
  std::vector<int> part(n_active_elem);
  
  int
    n = static_cast<int>(n_active_elem),  // number of "nodes" (elements)
                                          //   in the graph
    wgtflag = 2,                          // weights on vertices only,
                                          //   none on edges
    numflag = 0,                          // C-style 0-based numbering
    nparts  = static_cast<int>(n_pieces), // number of subdomains to create
    edgecut = 0;                          // the numbers of edges cut by the
                                          //   resulting partition
  
  // Set the options
  options[0] = 0; // use default options

  // Metis will only consider the active elements.
  // We need to map the active element ids into a
  // contiguous range.  Further, we want the unique range indexing to be
  // independednt of the element ordering, otherwise a circular dependency
  // can result in which the partitioning depends on the ordering which 
  // depends on the partitioning...
  std::map<const Elem*, unsigned int> global_index_map;
  {
    std::vector<unsigned int> global_index;
    
    MeshBase::element_iterator       it  = mesh.active_elements_begin();
    const MeshBase::element_iterator end = mesh.active_elements_end(); 
    
    MeshCommunication().find_global_indices (MeshTools::bounding_box(mesh),
					     it, end, global_index);
    
    libmesh_assert (global_index.size() == n_active_elem);    
    
    for (unsigned int cnt=0; it != end; ++it)
      {
	const Elem *elem = *it;
	libmesh_assert (!global_index_map.count(elem));
	
	global_index_map[elem]  = global_index[cnt++];
      }
    libmesh_assert (global_index_map.size() == n_active_elem);    
  }
  
  
  // build the graph in CSR format.  Note that
  // the edges in the graph will correspond to
  // face neighbors
  std::vector<int> xadj, adjncy;
  {
    std::vector<const Elem*> neighbors_offspring;    
    
    MeshBase::element_iterator       elem_it  = mesh.active_elements_begin();
    const MeshBase::element_iterator elem_end = mesh.active_elements_end(); 
    
    // This will be exact when there is no refinement and all the
    // elements are of the same type.
    unsigned int graph_size=0;
    std::vector<std::vector<unsigned int> > graph(n_active_elem);
    
    for (; elem_it != elem_end; ++elem_it)
      {
	const Elem* elem = *elem_it;
	
	libmesh_assert (global_index_map.count(elem));
	
	const unsigned int elem_global_index = 
	  global_index_map[elem];
	
	libmesh_assert (elem_global_index < vwgt.size());
	libmesh_assert (elem_global_index < graph.size());

	// maybe there is a better weight?
	// The weight is used to define what a balanced graph is
	vwgt[elem_global_index] = elem->n_nodes(); 

	// 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_map.count(neighbor));
		    
		    const unsigned int neighbor_global_index = 
		      global_index_map[neighbor];
		    
		    graph[elem_global_index].push_back(neighbor_global_index);
		    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_map.count(child));

			    const unsigned int child_global_index =
			      global_index_map[child];
			    
			    graph[elem_global_index].push_back(child_global_index);
			    graph_size++;
			  }
		      }
		  }

#endif /* ifdef ENABLE_AMR */

	      }
	  }
      }
    
    // Convert the graph into the format Metis wants
    xadj.reserve(n_active_elem+1);
    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 (adjncy.size() == graph_size);
    libmesh_assert (xadj.size() == n_active_elem+1);    
  } // done building the graph
  
  
  if (adjncy.empty())
    adjncy.push_back(0);
  
  // Select which type of partitioning to create
  
  // Use recursive if the number of partitions is less than or equal to 8
  if (n_pieces <= 8)
    Metis::METIS_PartGraphRecursive(&n, &xadj[0], &adjncy[0], &vwgt[0], NULL,
				    &wgtflag, &numflag, &nparts, &options[0],
				    &edgecut, &part[0]);
  
  // Otherwise  use kway
  else
    Metis::METIS_PartGraphKway(&n, &xadj[0], &adjncy[0], &vwgt[0], NULL,
			       &wgtflag, &numflag, &nparts, &options[0],
			       &edgecut, &part[0]);
  
  
  // Assign the returned processor ids.  The part array contains
  // the processor id for each active element, but in terms of
  // the contiguous indexing we defined above
  {
    MeshBase::element_iterator       it  = mesh.active_elements_begin();
    const MeshBase::element_iterator end = mesh.active_elements_end(); 
    
    for (; it!=end; ++it)
      {
	Elem* elem = *it;
	
	libmesh_assert (global_index_map.count(elem));
	
	const unsigned int elem_global_index = 
	  global_index_map[elem];
	
	libmesh_assert (elem_global_index < part.size());
	const unsigned int elem_procid =  static_cast<short int>(part[elem_global_index]);
    
        elem->processor_id() = elem_procid;
      }
  }

  STOP_LOG("partition()", "MetisPartitioner");
#endif
}