int generate_graph(Problem_Description* problem,
                   Mesh_Description<INT>* mesh,
                   Graph_Description<INT>* graph,
                   Weight_Description<INT>* weight,
                   Sphere_Info* sphere)
{
  double time1 = get_time();
  /* Find the elements surrounding a node */
  if(!find_surnd_elems(mesh, graph)) {
    Gen_Error(0, "fatal: could not find surrounding elements");
    return 0;
  }
  double time2 = get_time();
  printf("Time to find surrounding elements: %fs\n", time2-time1);

  /* Find the adjacency, if required */
  if(problem->alloc_graph == ELB_TRUE) {
    if(!find_adjacency(problem, mesh, graph, weight, sphere)) {
      Gen_Error(0, "fatal: could not find adjacency");
      return 0;
    }
    time1 = get_time();
    printf("Time to find the adjacency: %fs\n", time1-time2);
  }
  return 1;
}

static void create_Vtxdist()
{
/* Function that creates Vtxdist for the INITIAL_LINEAR distribution.
 * Vtxdist is an array of size Num_Proc+1 indicating range of
 * vertices assigned to a processor.  It is assumed that vertices are numbered 
 * consecutively.  The values stored in Vtxdist assume vertex IDs are 0-based 
 * (e.g., the lowest-numbered vertex has ID 0).  Proc p is assigned vertices 
 * Vtxdist[p] through (Vtxdist[p+1]-1), inclusive.         
 */
int rest, i, n;

  /* Set up Vtxdist data */
  if (Vtxdist == NULL){
    Vtxdist = (ZOLTAN_GNO_TYPE *) malloc((Num_Proc+1) * sizeof(ZOLTAN_GNO_TYPE));
    if (Vtxdist == NULL) {
      Gen_Error(0, "fatal: insufficient memory");
      return;
    }
  }
  /* Calculate uniform vertex distribution */
  Vtxdist[0] = 0;
  rest = Gnvtxs;
  for (i=0; i<Num_Proc_Dist; i++){
    n = rest/(Num_Proc_Dist-i);
    Vtxdist[i+1] = Vtxdist[i] + n;
    rest -= n;
  }
  /* Procs >= Num_Proc_Dist get no vertices */
  for (i=Num_Proc_Dist; i<Num_Proc; i++){
    Vtxdist[i+1] = Vtxdist[i];
  }

}

void ch_dist_vtx_list(
  int *vtx_list,
  int *nvtx,
  int target_proc,
  short *assignments          
)
{
/* Function that returns a list of vertices assigned to proc target_proc.
 * The list is returned in vtx_list.  The function assumes vertex ID
 * numbering is zero-based; e.g., the lowest-numbered vertex has ID 0.
 * This convention allows the list entries to be used as array indices
 * in the arrays containing chaco input file data.
 * The number of entries in the list is returned in nvtx.
 */
int i, proc;

  *nvtx = 0;
  proc = (ZOLTAN_GNO_TYPE)target_proc;

  switch(Initial_Method) {
  case INITIAL_FILE:
    for (i = 0; i < Gnvtxs; i++)
      if (assignments[i] == target_proc)
        vtx_list[(*nvtx)++] = i;
    break;
  case INITIAL_LINEAR:
  case INITIAL_OWNER:
    for (i = Vtxdist[target_proc]; i < Vtxdist[target_proc+1]; i++)
      vtx_list[(*nvtx)++] = i;
    break;
  case INITIAL_CYCLIC:
    if (target_proc < Num_Proc_Dist){
      for (i = proc; i < Gnvtxs; i+=Num_Proc_Dist) 
        vtx_list[(*nvtx)++] = i;
    }
    break;
  default:
    Gen_Error(0, "Invalid Initial Distribution Type in ch_dist_vtx_list");
    return;
  }
}

const char *get_elem_name(int itype) {
/* Function to return the name of an element given its type.
 * Inverse of get_elem_type().
 */
E_Type etype = (E_Type) itype;

  switch (etype) {
  case SPHERE:
    return("SPHERE");
  case BAR1:
  case BAR2:
    return("BAR");
  case QUAD1:
  case S_QUAD2:
  case QUAD2:
    return("QUAD");
  case HEX1:
  case HEXSHELL:
  case S_HEX2:
  case HEX2:
    return("HEX");
  case TRI1:
  case TSHELL1:
  case TRI2:
  case TSHELL2:
    return("TRI");
  case TET1:
  case TET2:
    return("TET");
  case SHELL1:
  case SHELL2:
    return("SHELL");
  case WEDGE1:
  case WEDGE2:
    return("WEDGE");
  default:
    Gen_Error(0, "fatal: unknown element type read");
    return("NULL");
  }
}

int ch_dist_num_vtx(
  int target_proc,
  short *assignments
)
{
/* Function that returns the number of vertices assigned to processor
 * target_proc.
 */
int num;

  switch(Initial_Method) {
  case INITIAL_FILE: {
    int i;
    num = 0;
    for (i = 0; i < Gnvtxs; i++)
      if (assignments[i] == target_proc) num++;
    break;
  }
  case INITIAL_LINEAR:
  case INITIAL_OWNER:
    num = Vtxdist[target_proc+1] - Vtxdist[target_proc];
    break;
  case INITIAL_CYCLIC:
    if (target_proc < Num_Proc_Dist){
      num = Gnvtxs / Num_Proc_Dist;
      if ((Gnvtxs % Num_Proc_Dist) > target_proc)
        num++;
    }
    else
      num = 0;
    break;
  default:
    Gen_Error(0, "Invalid Initial Distribution Type in ch_dist_num_vtx");
    return(-1);
  }
  return(num);
}

int ch_dist_proc(int v, short *assignments, int base)
{
/* Function that returns the processor to which a vertex v is assigned.
 * The function assumes the vertex numbering is "base"-based (i.e., lowest 
 * numbered vertex is vertex base); this convention is used since the function
 * is called primarily to find adjacent vertices' processor assignments and
 * the read-in adjacency information is "base"-based.
 * E.g., for Chaco input files, base == 1; vertex comparisons are one-based.
 *       for HG input files, base may be 0 or 1.
 */
int p;
ZOLTAN_GNO_TYPE b = (ZOLTAN_GNO_TYPE)base;

  switch(Initial_Method) {
  case INITIAL_FILE:
    /* return the appropriate entry from the assignments array. */
    p = assignments[v-b];
    break;
  case INITIAL_LINEAR:
  case INITIAL_OWNER:
    for (p = 0; p < Num_Proc_Dist; p++)
      /* Since v is "base"-based and Vtxdist is 0-based, add base to 
       * Vtxdist[p+1]. */
      if (v < Vtxdist[p+1]+b) break;
    break;
  case INITIAL_CYCLIC:
    /* test for (v-base) as v is "base"-based and 
     * INITIAL_CYCLIC equations are 0-based */
    p = (v-b) % Num_Proc_Dist;
    break;
  default:
    Gen_Error(0, "Invalid Initial Distribution Type in ch_dist_proc");
    return -1;
  }
  return p;
}

void migrate_unpack_elem(void *data, int num_gid_entries, ZOLTAN_ID_PTR elem_gid, 
                         int elem_data_size, char *buf, int *ierr)
{
  int gid = num_gid_entries-1;

  if (data == NULL) {
    *ierr = ZOLTAN_FATAL;
    return;
  }
  MESH_INFO_PTR mesh = (MESH_INFO_PTR) data;
  ELEM_INFO *elem = mesh->elements;
  ELEM_INFO *elem_mig = (ELEM_INFO *) buf;

  int proc = 0;
  MPI_Comm_rank(MPI_COMM_WORLD, &proc);

  int idx = 0;
  ZOLTAN_ID_TYPE egid = elem_gid[gid];

  if ((idx = in_list(egid, New_Elem_Index_Size, New_Elem_Index)) == -1) {
    Gen_Error(0, "fatal: Unable to locate position for element");
    *ierr = ZOLTAN_FATAL;
    return;
  }

  ELEM_INFO *current_elem = &(elem[idx]);
  /* now put the migrated information into the array */
  *current_elem = *elem_mig;
  int num_nodes = mesh->eb_nnodes[current_elem->elem_blk];

  int size = sizeof(ELEM_INFO);

  /*
   * copy the allocated integer fields for this element.
   */

  /* Pad the buffer so the following casts will work.  */
  size = Zoltan_Align(size);
  ZOLTAN_ID_TYPE *buf_id_type = (ZOLTAN_ID_TYPE *) (buf + size);

  /* copy the connect table */
  if (mesh->num_dims > 0) {

    // Don't use C++ new/delete here becuase these items are
    // malloc'd/free'd elsewhere in C code

    current_elem->connect = (ZOLTAN_ID_TYPE *) malloc(num_nodes * sizeof(ZOLTAN_ID_TYPE));
    if (current_elem->connect == NULL) {
      Gen_Error(0, "fatal: insufficient memory");
      *ierr = ZOLTAN_MEMERR;
      return;
    }
    for (int i = 0; i < num_nodes; i++) {
      current_elem->connect[i] = *buf_id_type++;
    }
    size += num_nodes * sizeof(ZOLTAN_ID_TYPE);
  }

  /* copy the adjacency info */
  /* globalIDs are received; convert to local IDs when adj elem is local */

  float *buf_float = NULL;
  int *buf_int = NULL;

  if (current_elem->adj_len > 0) {
    size_t adjgids_len = current_elem->adj_len * sizeof(ZOLTAN_ID_TYPE);
    size_t adjprocs_len = current_elem->adj_len * sizeof(int);

    current_elem->adj      = (ZOLTAN_ID_TYPE *)malloc(adjgids_len);
    current_elem->adj_proc = (int *)malloc(adjprocs_len);

    if (current_elem->adj == NULL || current_elem->adj_proc == NULL) {
      Gen_Error(0, "fatal: insufficient memory");
      *ierr = ZOLTAN_MEMERR;
      return;
    }

    buf_id_type = (ZOLTAN_ID_TYPE *) (buf + size);
    buf_int    = (int *) (buf + size + adjgids_len);

    for (int i =  0; i < current_elem->adj_len; i++) {
      current_elem->adj[i] =      *buf_id_type++;
      current_elem->adj_proc[i] = *buf_int++;

      if (current_elem->adj[i] != ZOLTAN_ID_INVALID && current_elem->adj_proc[i] == proc) {
        idx = in_list(current_elem->adj[i], New_Elem_Index_Size, New_Elem_Index);
        if (idx < 0)
          current_elem->adj[i] = ZOLTAN_ID_INVALID;
        else
          current_elem->adj[i] = (ZOLTAN_ID_TYPE)idx;
      }
    }
    size += (adjgids_len + adjprocs_len);

    /* copy the edge_wgt data */
    if (Use_Edge_Wgts) {

      /* Pad the buffer so the following casts will work.  */
      size = Zoltan_Align(size);
      buf_float = (float *) (buf + size);
//.........这里部分代码省略.........

int read_mesh(const std::string &exo_file,
              Problem_Description* problem,
              Mesh_Description<INT>* mesh,
              Weight_Description<INT>* weight
	      )
{
  float  version, *xptr, *yptr, *zptr;
  char   elem_type[MAX_STR_LENGTH+1];
  E_Type blk_elem_type;
  
  /*---------------------------Execution Begins--------------------------------*/

  /* Open the ExodusII file */
  int exoid, cpu_ws=0, io_ws=0;
  int mode = EX_READ | problem->int64api;
  if((exoid=ex_open(exo_file.c_str(), mode, &cpu_ws, &io_ws, &version)) < 0)
    {
      Gen_Error(0, "fatal: unable to open ExodusII mesh file");
      return 0;
    }

  /* Read the coordinates, if desired */
  xptr = yptr = zptr = NULL;
  if(problem->read_coords == ELB_TRUE)
    {
      switch(mesh->num_dims)
	{
	case 3:
	  zptr = (mesh->coords)+2*(mesh->num_nodes);
	  /* FALLTHRU */
	case 2:
	  yptr = (mesh->coords)+(mesh->num_nodes);
	  /* FALLTHRU */
	case 1:
	  xptr = mesh->coords;
	}

      if(ex_get_coord(exoid, xptr, yptr, zptr) < 0)
	{
	  Gen_Error(0, "fatal: unable to read coordinate values for mesh");
	  return 0;
	}

    } /* End "if(problem->read_coords == ELB_TRUE)" */

  /* Read the element block IDs */
  std::vector<INT> el_blk_ids(mesh->num_el_blks);
  std::vector<INT> el_blk_cnts(mesh->num_el_blks);

  if(ex_get_elem_blk_ids(exoid, &el_blk_ids[0]) < 0)
    {
      Gen_Error(0, "fatal: unable to read element block IDs");
      return 0;
    }

  /* Read the element connectivity */
  size_t gelem_cnt=0;
  for(size_t cnt=0; cnt < mesh->num_el_blks; cnt++) {
    INT nodes_per_elem, num_attr;
    if(ex_get_elem_block(exoid, el_blk_ids[cnt], elem_type,
                         &(el_blk_cnts[cnt]), &nodes_per_elem,
                         &num_attr) < 0)
      {
	Gen_Error(0, "fatal: unable to read element block");
	return 0;
      }

    blk_elem_type = get_elem_type(elem_type, nodes_per_elem, mesh->num_dims);

    INT *blk_connect = (INT*)malloc(sizeof(INT)*el_blk_cnts[cnt]*nodes_per_elem);
    if(!blk_connect)
      {
	Gen_Error(0, "fatal: insufficient memory");
	return 0;
      }

    /* Get the connectivity for this element block */
    if(ex_get_elem_conn(exoid, el_blk_ids[cnt], blk_connect) < 0)
      {
	Gen_Error(0, "fatal: failed to get element connectivity");
	return 0;
      }

    /* find out if this element block is weighted */
    int wgt = -1;
    if (weight->type & EL_BLK)
      wgt = in_list(el_blk_ids[cnt], weight->elemblk);
    
    /* Fill the 2D global connectivity array */
    if (((problem->type == ELEMENTAL) && (weight->type & EL_BLK)) ||
        ((problem->type == NODAL) && (weight->type & EL_BLK))) {
      
      for(int64_t cnt2=0; cnt2 < el_blk_cnts[cnt]; cnt2++) {
	mesh->elem_type[gelem_cnt] = blk_elem_type;
      
	/* while going through the blocks, take care of the weighting */
	if ((problem->type == ELEMENTAL) && (weight->type & EL_BLK)) {
	  /* is this block weighted */
	  if (wgt >= 0) {
	    /* check if there is a read value */
//.........这里部分代码省略.........

int write_vis(std::string &nemI_out_file,
	      std::string &exoII_inp_file,
	      Machine_Description* machine,
	      Problem_Description* prob,
	      Mesh_Description<INT>* mesh,
	      LB_Description<INT>* lb)
{
  int    exid_vis, exid_inp;

  char  title[MAX_LINE_LENGTH+1];
  const char   *coord_names[] = {"X", "Y", "Z"};

  /*-----------------------------Execution Begins------------------------------*/

  /* Generate the file name for the visualization file */
  std::string vis_file_name = remove_extension(nemI_out_file);
  vis_file_name += "-vis.exoII";

  /* Generate the title for the file */
  strcpy(title, UTIL_NAME);
  strcat(title, " ");
  strcat(title, ELB_VERSION);
  strcat(title, " load balance visualization file");

  /*
   * If the vis technique is to be by element block then calculate the
   * number of element blocks.
   */
  int    vis_nelem_blks;
  if(prob->type == ELEMENTAL)
    vis_nelem_blks = machine->num_procs;
  else
    vis_nelem_blks = machine->num_procs + 1;

  /* Create the ExodusII file */
  std::cout << "Outputting load balance visualization file " << vis_file_name.c_str() << "\n";
  int cpu_ws = 0;
  int io_ws = 0;
  int mode = EX_CLOBBER;
  if (prob->int64db|prob->int64api) {
    mode |= EX_NETCDF4|EX_NOCLASSIC|prob->int64db|prob->int64api;
  }
  if((exid_vis=ex_create(vis_file_name.c_str(), mode, &cpu_ws, &io_ws)) < 0) {
    Gen_Error(0, "fatal: unable to create visualization output file");
    return 0;
  }
  ON_BLOCK_EXIT(ex_close, exid_vis);

  /*
   * Open the original input ExodusII file, read the values for the
   * element blocks and output them to the visualization file.
   */
  int icpu_ws=0;
  int iio_ws=0;
  float vers=0.0;
  mode = EX_READ | prob->int64api;
  if((exid_inp=ex_open(exoII_inp_file.c_str(), mode, &icpu_ws, &iio_ws, &vers)) < 0) {
    Gen_Error(0, "fatal: unable to open input ExodusII file");
    return 0;
  }
  ON_BLOCK_EXIT(ex_close, exid_inp);
  
  char **elem_type  = (char**)array_alloc(2, mesh->num_el_blks, MAX_STR_LENGTH+1,
					  sizeof(char));
  if(!elem_type) {
    Gen_Error(0, "fatal: insufficient memory");
    return 0;
  }
  ON_BLOCK_EXIT(free, elem_type);

  std::vector<INT> el_blk_ids(mesh->num_el_blks);
  std::vector<INT> el_cnt_blk(mesh->num_el_blks);
  std::vector<INT> node_pel_blk(mesh->num_el_blks);
  std::vector<INT> nattr_el_blk(mesh->num_el_blks);

  if(ex_get_elem_blk_ids(exid_inp, TOPTR(el_blk_ids)) < 0) {
    Gen_Error(0, "fatal: unable to get element block IDs");
    return 0;
  }

  int acc_vis = ELB_TRUE; // Output a different element block per processor
  if (prob->vis_out == 2)
    acc_vis = ELB_FALSE; // Output a nodal/element variable showing processor

  size_t nsize = 0;

  /*
   * Find out if the mesh consists of mixed elements. If not then
   * element blocks will be used to visualize the partitioning. Otherwise
   * nodal/element results will be used.
   */
  for(size_t ecnt=0; ecnt < mesh->num_el_blks; ecnt++) {
    if(ex_get_elem_block(exid_inp, el_blk_ids[ecnt], elem_type[ecnt],
			 &el_cnt_blk[ecnt], &node_pel_blk[ecnt],
			 &nattr_el_blk[ecnt]) < 0) {
      Gen_Error(0, "fatal: unable to get element block parameters");
      return 0;
    }

    nsize += el_cnt_blk[ecnt]*node_pel_blk[ecnt];
//.........这里部分代码省略.........

void migrate_pre_process(void *data, int num_gid_entries, int num_lid_entries, 
                         int num_import, 
                         ZOLTAN_ID_PTR import_global_ids,
                         ZOLTAN_ID_PTR import_local_ids, int *import_procs,
                         int *import_to_part,
                         int num_export, ZOLTAN_ID_PTR export_global_ids,
                         ZOLTAN_ID_PTR export_local_ids, int *export_procs,
                         int *export_to_part,
                         int *ierr)
{
int lid = num_lid_entries-1;
int gid = num_gid_entries-1;
char msg[256];

  *ierr = ZOLTAN_OK;

  if (data == NULL) {
    *ierr = ZOLTAN_FATAL;
    return;
  }
  MESH_INFO_PTR mesh = (MESH_INFO_PTR) data;
  ELEM_INFO_PTR elements = mesh->elements;

  /*
   *  Set some flags. Assume if true for one element, true for all elements.
   *  Note that some procs may have no elements. 
   */

  int k = 0;

  if (elements[0].edge_wgt != NULL)
    k = 1;

  /* Make sure all procs have the same value */

  MPI_Allreduce(&k, &Use_Edge_Wgts, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);

  /*
   *  For all elements, update adjacent elements' processor information.
   *  That way, when perform migration, will be migrating updated adjacency
   *  information.  
   */
  
  int proc = 0;
  MPI_Comm_rank(MPI_COMM_WORLD, &proc);

  /*
   *  Build New_Elem_Index array and list of processor assignments.
   */

  New_Elem_Index_Size = mesh->num_elems + num_import - num_export;
  if (mesh->elem_array_len > New_Elem_Index_Size) 
    New_Elem_Index_Size = mesh->elem_array_len;

  New_Elem_Index = new ZOLTAN_ID_TYPE [New_Elem_Index_Size];

  int *proc_ids = NULL;
  char *change = NULL;

  if (mesh->num_elems > 0) {

    proc_ids = new int [mesh->num_elems];
    change   = new char [mesh->num_elems];

    if (New_Elem_Index == NULL || proc_ids == NULL || change == NULL) {
      Gen_Error(0, "fatal: insufficient memory");
      *ierr = ZOLTAN_MEMERR;
      if (proc_ids) delete [] proc_ids;
      if (change) delete [] change;
      if (New_Elem_Index)
        {
        delete [] New_Elem_Index;
        New_Elem_Index = NULL;
        }
      return;
    }

    for (int i = 0; i < mesh->num_elems; i++) {
      New_Elem_Index[i] = elements[i].globalID;
      proc_ids[i] = proc;
      change[i] = 0;
    }
  }

  for (int i = mesh->num_elems; i < New_Elem_Index_Size; i++) {
    New_Elem_Index[i] = ZOLTAN_ID_INVALID;
  }

  for (int i = 0; i < num_export; i++) {

    int exp_elem = 0;

    if (num_lid_entries)
      exp_elem = export_local_ids[lid+i*num_lid_entries];
    else  /* testing num_lid_entries == 0 */
      search_by_global_id(mesh, export_global_ids[gid+i*num_gid_entries], 
                          &exp_elem);

    if (export_procs[i] != proc) {
      /* Export is moving to a new processor */
//.........这里部分代码省略.........

int read_exo_weights(Problem_Description* prob, Weight_Description<INT>* weight)
{
  int    exoid, cpu_ws=0, io_ws=0;
  int    neblks;
  float  version, minval = 1.0f;
  char elem_type[MAX_STR_LENGTH+1];
  char ctemp[1024];
/*---------------------------Execution Begins--------------------------------*/

  /* Open the ExodusII file containing the weights */
  int mode = EX_READ | prob->int64api;
  if((exoid=ex_open(weight->exo_filename.c_str(), mode, &cpu_ws, &io_ws,
                    &version)) < 0)
  {
    sprintf(ctemp, "fatal: could not open ExodusII file %s",
            weight->exo_filename.c_str());
    Gen_Error(0, ctemp);
    return 0;
  }

  std::vector<float> values(weight->nvals);
  if(prob->type == NODAL)
  {
    size_t tmp_nodes = ex_inquire_int(exoid, EX_INQ_NODES);
    /* check to make sure the sizes agree */
    if ((size_t)weight->nvals != tmp_nodes) {
      Gen_Error(0, "fatal: different number of nodes in mesh and weight files");
      ex_close(exoid);
      return 0;
    }

    weight->ow.resize(weight->nvals);
    /* Read in the nodal values */
    if(ex_get_nodal_var(exoid, weight->exo_tindx, weight->exo_vindx,
                        weight->nvals, TOPTR(values)) < 0)
    {
      Gen_Error(0, "fatal: unable to read nodal values");
      ex_close(exoid);
      return 0;
    }
  }
  else
  {
    size_t tmp_elem = ex_inquire_int(exoid, EX_INQ_ELEM);
    /* check to make sure the sizes agree */
    if ((size_t)weight->nvals != tmp_elem) {
      Gen_Error(0, "fatal: different number of elems in mesh and weight files");
      ex_close(exoid);
      return 0;
    }

    /* Get the number of element blocks */
    neblks = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
    std::vector<INT> eblk_ids(neblks);
    std::vector<INT> eblk_ecnts(neblks);

    if(ex_get_ids(exoid, EX_ELEM_BLOCK, &eblk_ids[0]) < 0)
    {
      Gen_Error(0, "fatal: unable to get element block IDs");
      ex_close(exoid);
      return 0;
    }

    /* Get the count of elements in each element block */
    for(int cnt=0; cnt < neblks; cnt++) {
      INT dum1, dum2;
      if(ex_get_elem_block(exoid, eblk_ids[cnt], elem_type,
                           &(eblk_ecnts[cnt]), &dum1, &dum2) < 0)
      {
        Gen_Error(0, "fatal: unable to get element block");
        ex_close(exoid);
        return 0;
      }
    }

    /* Get the element variables */
    size_t offset = 0;
    for(int cnt=0; cnt < neblks; cnt++)
    {
      if(ex_get_elem_var(exoid, weight->exo_tindx, weight->exo_vindx,
                         eblk_ids[cnt], eblk_ecnts[cnt], &(values[offset])) < 0)
      {
        Gen_Error(0, "fatal: unable to get element variable");
        ex_close(exoid);
        return 0;
      }
      offset += eblk_ecnts[cnt];
    }
  }

  /* Close the ExodusII weighting file */
  if(ex_close(exoid) < 0)
  {
    sprintf(ctemp, "warning: failed to close ExodusII file %s",
            weight->exo_filename.c_str());
    Gen_Error(0, ctemp);
  }

  /* now I need to translate the values to positive integers */

//.........这里部分代码省略.........

int write_nemesis(std::string &nemI_out_file,
                  Machine_Description* machine,
                  Problem_Description* problem,
                  Mesh_Description<INT>* mesh,
                  LB_Description<INT>* lb,
                  Sphere_Info* sphere)
{
  int     exoid;
  char    title[MAX_LINE_LENGTH+1], method1[MAX_LINE_LENGTH+1];
  char    method2[MAX_LINE_LENGTH+1];

  int cpu_ws = sizeof(float);
  int io_ws  = sizeof(float);

  printf("Outputting load balance to file %s\n", nemI_out_file.c_str());

  /* Create the load balance file */
  /* Attempt to create a netcdf4-format file; if it fails, then assume
     that the netcdf library does not support that mode and fall back
     to classic netcdf3 format.  If that fails, issue an error and
     return failure.
  */
  int mode3 = EX_CLOBBER;
  int mode4 = mode3|EX_NETCDF4|EX_NOCLASSIC|problem->int64db|problem->int64api;

  ex_opts(EX_DEFAULT); // Eliminate misleading error if the first ex_create fails, but the second succeeds.
  if((exoid=ex_create(nemI_out_file.c_str(), mode4, &cpu_ws, &io_ws)) < 0) {
    /* If int64api or int64db non-zero, then netcdf-4 format is required, so
       fail now...
    */
    if (problem->int64db|problem->int64api) {
      Gen_Error(0, "fatal: failed to create Nemesis netcdf-4 file");
      return 0;
    }
    if((exoid=ex_create(nemI_out_file.c_str(), mode3, &cpu_ws, &io_ws)) < 0) {
      Gen_Error(0, "fatal: failed to create Nemesis file");
      return 0;
    }
  }
  ON_BLOCK_EXIT(ex_close, exoid);
  
  /* Set the error reporting value */
  if (error_lev > 1)
    ex_opts(EX_VERBOSE | EX_DEBUG);
  else
    ex_opts(EX_VERBOSE);

  /* Enable compression (if netcdf-4) */
  ex_set_option(exoid, EX_OPT_COMPRESSION_LEVEL, 1);
  ex_set_option(exoid, EX_OPT_COMPRESSION_SHUFFLE, 1);

  /* Create the title */
  if(problem->type == NODAL)
    strcpy(method1, "nodal");
  else
    strcpy(method1, "elemental");

  sprintf(title, "nem_slice %s load balance file", method1);

  strcpy(method1, "method1: ");
  strcpy(method2, "method2: ");

  switch(lb->type)
    {
    case MULTIKL:
      strcat(method1, "Multilevel-KL decomposition");
      strcat(method2, "With Kernighan-Lin refinement");
      break;
    case SPECTRAL:
      strcat(method1, "Spectral decomposition");
      break;
    case INERTIAL:
      strcat(method1, "Inertial decomposition");
      break;
    case ZPINCH:
      strcat(method1, "ZPINCH decomposition");
      break;
    case BRICK:
      strcat(method1, "BRICK decomposition");
      break;
    case ZOLTAN_RCB:
      strcat(method1, "RCB decomposition");
      break;
    case ZOLTAN_RIB:
      strcat(method1, "RIB decomposition");
      break;
    case ZOLTAN_HSFC:
      strcat(method1, "HSFC decomposition");
      break;
    case LINEAR:
      strcat(method1, "Linear decomposition");
      break;
    case RANDOM:
      strcat(method1, "Random decomposition");
      break;
    case SCATTERED:
      strcat(method1, "Scattered decomposition");
      break;
    }

//.........这里部分代码省略.........

int get_group_info(Machine_Description* machine,
                   Problem_Description* prob,
                   Mesh_Description<INT>* mesh,
                   Graph_Description<INT>* graph,
                   int elem2grp[],
                   int nprocg[],
                   int nelemg[],
                   size_t *max_vtx,
                   size_t *max_adj
  )
{
  int nproc=0;
  std::vector<int> nadj_per_grp;

  /* allocate array to hold adjacency counts, if necessary */
  if (prob->alloc_graph == ELB_TRUE) {
    nadj_per_grp.resize(prob->num_groups);
  }

  /* initialize the group counts arrays */
  for (int i = 0; i < prob->num_groups; i++) {
    nelemg[i] = 0;
  }

  /*
   * fill the vertex2proc array with the group number of the individual
   * elements, calculate how many elements in each group and determine
   * how many adjacencies are in each group (if necessary).
   */
  INT sum = 0;
  int iblk = 0;
  for (size_t i = 0; i < prob->num_vertices; i++) {

    /* figure out which element block this is */
    if (sum == mesh->eb_cnts[iblk]) {
      sum = 0;
      iblk++;
    }
    sum++;

    /*
     * use negative numbers to specify the groups, in order to
     * avoid having a problem with group 0, add 1 to the group
     * number
     */
    elem2grp[i] =  -(prob->group_no[iblk] + 1);

    nelemg[prob->group_no[iblk]]++;

    if (prob->alloc_graph == ELB_TRUE)
      nadj_per_grp[prob->group_no[iblk]] += graph->start[i+1] - graph->start[i];
  }

  /*
   * calculate how many processors to use for each group
   *   using method from the materials group, haven't really checked it
   */
  if (machine->type == MESH)
    nproc = machine->procs_per_box;
  else if (machine->type == HCUBE)
    nproc = ilog2i(machine->procs_per_box);
  for (int i = 0; i < prob->num_groups; i++) {
    nprocg[i] = int((nproc * (nelemg[i] + .5)) / (float) prob->num_vertices);
    if (nelemg[i] && !nprocg[i]) nprocg[i] = 1;
  }

  /*
   * check to see if correct number of processors have been allocated
   * and get the maximum number of vertices
   */
  sum = 0;
  size_t j = 0;
  *max_vtx = 0;
  *max_adj = 0;
  for (int i = 0; i < prob->num_groups; i++) {
    sum += nprocg[i];
    if (nprocg[i] > nprocg[j]) {
      j = i;
      *max_vtx = nelemg[j];  /* most processors implies most elements */
    }

    /* determine how large to make temporary arrays */
    if ((size_t)nelemg[i] > *max_vtx) *max_vtx = nelemg[i];
    if(prob->alloc_graph == ELB_TRUE)
      if ((size_t)nadj_per_grp[i] > *max_adj) *max_adj = nadj_per_grp[i];
  }

  if (sum != nproc) {
    /* correct group with most processors (j determined above) */
    nprocg[j] -= (sum - nproc);
    if (nprocg[j] <= 0) {
      Gen_Error(0,"Unable to balance # processors in get_group_info().");
      return 0;
    }
  }

  printf("Load balance information\n");
  for (int i = 0; i < prob->num_groups; i++)
    printf("group[%d]  #elements=%-10d  #proc=%d\n",i,nelemg[i],nprocg[i]);

//.........这里部分代码省略.........

void migrate_post_process(void *data, int num_gid_entries, int num_lid_entries,
                          int num_import, 
                          ZOLTAN_ID_PTR import_global_ids,
                          ZOLTAN_ID_PTR import_local_ids, int *import_procs,
                          int *import_to_part,
                          int num_export, ZOLTAN_ID_PTR export_global_ids,
                          ZOLTAN_ID_PTR export_local_ids, int *export_procs,
                          int *export_to_part,
                          int *ierr)
{
  if (data == NULL) {
    *ierr = ZOLTAN_FATAL;
    return;
  }
  MESH_INFO_PTR mesh = (MESH_INFO_PTR) data;
  ELEM_INFO *elements = mesh->elements;

  int proc = 0, num_proc = 0;
  MPI_Comm_rank(MPI_COMM_WORLD, &proc);
  MPI_Comm_size(MPI_COMM_WORLD, &num_proc);

  /* compact elements array, as the application expects the array to be dense */
  for (int i = 0; i < New_Elem_Index_Size; i++) {
    if (New_Elem_Index[i] != ZOLTAN_ID_INVALID) continue;

    /* Don't want to shift all elements down one position to fill the  */
    /* blank spot -- too much work to adjust adjacencies!  So find the */
    /* last element in the array and move it to the blank spot.        */

    int last = 0;

    for (last = New_Elem_Index_Size-1; last >= 0; last--)
      if (New_Elem_Index[last] != ZOLTAN_ID_INVALID) break;

    /* If (last < i), array is already dense; i is just in some blank spots  */
    /* at the end of the array.  Quit the compacting.                     */
    if (last < i) break;

    /* Copy elements[last] to elements[i]. */
    elements[i] = elements[last];

    /* Adjust adjacencies for local elements.  Off-processor adjacencies */
    /* don't matter here.                                                */

    for (int j = 0; j < elements[i].adj_len; j++) {

      /* Skip NULL adjacencies (sides that are not adjacent to another elem). */
      if (elements[i].adj[j] == ZOLTAN_ID_INVALID) continue;

      ZOLTAN_ID_TYPE adj_elem = elements[i].adj[j];

      /* See whether adjacent element is local; if so, adjust its entry */
      /* for local element i.                                           */
      if (elements[i].adj_proc[j] == proc) {
        for (int k = 0; k < elements[adj_elem].adj_len; k++) {
          if (elements[adj_elem].adj[k] == (ZOLTAN_ID_TYPE)last &&
              elements[adj_elem].adj_proc[k] == proc) {
            /* found adjacency entry for element last; change it to i */
            elements[adj_elem].adj[k] = (ZOLTAN_ID_TYPE)i;
            break;
          }
        }
      }
    }

    /* Update New_Elem_Index */
    New_Elem_Index[i] = New_Elem_Index[last];
    New_Elem_Index[last] = ZOLTAN_ID_INVALID;

    /* clear elements[last] */
    elements[last].globalID = ZOLTAN_ID_INVALID;
    elements[last].border = 0;
    elements[last].my_part = -1;
    elements[last].perm_value = -1;
    elements[last].invperm_value = -1;
    elements[last].nadj = 0;
    elements[last].adj_len = 0;
    elements[last].elem_blk = -1;
    for (int k=0; k<MAX_CPU_WGTS; k++)
      elements[last].cpu_wgt[k] = 0;
    elements[last].mem_wgt = 0;
    elements[last].avg_coord[0] = elements[last].avg_coord[1] 
                                = elements[last].avg_coord[2] = 0.;
    elements[last].coord = NULL;
    elements[last].connect = NULL;
    elements[last].adj = NULL;
    elements[last].adj_proc = NULL;
    elements[last].edge_wgt = NULL;
  }

  if (New_Elem_Index != NULL) {
    delete [] New_Elem_Index;
    New_Elem_Index = NULL;
   } 

  New_Elem_Index_Size = 0;

  if (!build_elem_comm_maps(proc, mesh)) {
    Gen_Error(0, "Fatal: error rebuilding elem comm maps");
  }
//.........这里部分代码省略.........

int parse_groups(INT *el_blk_ids,
                 INT *el_blk_cnts,
                 Mesh_Description<INT>* mesh,
                 Problem_Description* prob
  )
{
  char *id;
  int   last, found;

/*---------------------------Execution Begins--------------------------------*/

  /* allocate memory for the groups */
  prob->group_no = (int *) malloc (mesh->num_el_blks * sizeof(int));
  mesh->eb_cnts = (INT *) malloc (mesh->num_el_blks * sizeof(INT));
  if (!(prob->group_no) || !(mesh->eb_cnts))
  {
    Gen_Error(0, "fatal: insufficient memory");
    return 0;
  }

  /* prepare the group number array, and copy the element block counts */
  for (size_t i = 0; i < mesh->num_el_blks; i++) {
    prob->group_no[i] = -1;
    mesh->eb_cnts[i] = el_blk_cnts[i];
  }

  /* convert any comma's to blank spaces in the designator string */
  for (size_t i = 0; i < strlen(prob->groups); i++)
    if (prob->groups[i] == ',') prob->groups[i] = ' ';

  /* fill in the group identifier for each block */
  id = prob->groups;
  size_t i = 0;
  do {
    if (*id == '/') id++;
    scandescriptor(id, el_blk_ids, i, mesh->num_el_blks, prob);
    id = strchr(id, '/');
    i++;
  } while (id != NULL);
  last = i;

  /* set any remaining blocks to new group */
  found = 0;
  for (i = 0; i < mesh->num_el_blks; i++)
    if (prob->group_no[i] < 0) {
      prob->group_no[i] = last;
      found = 1;
     }

  if (found) last++;

  prob->num_groups = last;

  {
    size_t first_el = 0;
    printf("\nNumber of blocks: %lu\n", mesh->num_el_blks);
    printf("Block ID and associated groups:\n");
    printf("   block   #elems  group   type\n");
    for (i = 0; i < mesh->num_el_blks; i++) {
      printf("%8lu%8lu%8d%8s\n", (size_t)el_blk_ids[i], (size_t)mesh->eb_cnts[i], prob->group_no[i],
	     elem_name_from_enum(mesh->elem_type[first_el]));
      first_el += mesh->eb_cnts[i];
    }
    printf("There are %d groups of blocks\n", prob->num_groups);
  }

  /* finnished with the group designator string */
  free (prob->groups);

  return 1;

}

int migrate_elements(
  int Proc,
  MESH_INFO_PTR mesh,
  Zoltan &zz,
  int num_gid_entries, 
  int num_lid_entries,
  int num_imp,
  ZOLTAN_ID_PTR imp_gids,
  ZOLTAN_ID_PTR imp_lids,
  int *imp_procs,
  int *imp_to_part,
  int num_exp,
  ZOLTAN_ID_PTR exp_gids,
  ZOLTAN_ID_PTR exp_lids,
  int *exp_procs,
  int *exp_to_part)
{
/* Local declarations. */
const char *yo = "migrate_elements";

/***************************** BEGIN EXECUTION ******************************/
  DEBUG_TRACE_START(Proc, yo);

  /*
   * register migration functions
   */
  if (!Test.Null_Lists) {
    /* If not passing NULL lists, let Help_Migrate call the
     * pre-processing and post-processing routines.
     */
    if (zz.Set_Pre_Migrate_PP_Fn(migrate_pre_process,
                      (void *) mesh) == ZOLTAN_FATAL) {
      Gen_Error(0, "fatal:  error returned from Set_Pre_Migrate_PP_Fn()\n");
      return 0;
    }

    if (zz.Set_Post_Migrate_PP_Fn(migrate_post_process,
                      (void *) mesh) == ZOLTAN_FATAL) {
      Gen_Error(0, "fatal:  error returned from Set_Post_Migrate_PP_Fn()\n");
      return 0;
    }
  }

  if (Test.Multi_Callbacks) {
    if (zz.Set_Obj_Size_Multi_Fn(migrate_elem_size_multi,
                      (void *) mesh) == ZOLTAN_FATAL) {
      Gen_Error(0, "fatal:  error returned from Set_Obj_Size_Multi_Fn()\n");
      return 0;
    }

    if (zz.Set_Pack_Obj_Multi_Fn(migrate_pack_elem_multi,
                      (void *) mesh) == ZOLTAN_FATAL) {
      Gen_Error(0, "fatal:  error returned from Set_Pack_Obj_Multi_Fn()\n");
      return 0;
    }
  
    if (zz.Set_Unpack_Obj_Multi_Fn(migrate_unpack_elem_multi,
                      (void *) mesh) == ZOLTAN_FATAL) {
      Gen_Error(0, "fatal:  error returned from Set_Unpack_Obj_Multi_Fn()\n");
      return 0;
    }
  }
  else {
    if (zz.Set_Obj_Size_Fn(migrate_elem_size,
                      (void *) mesh) == ZOLTAN_FATAL) {
      Gen_Error(0, "fatal:  error returned from Set_Obj_Size_Fn()\n");
      return 0;
    }

    if (zz.Set_Pack_Obj_Fn(migrate_pack_elem,
                      (void *) mesh) == ZOLTAN_FATAL) {
      Gen_Error(0, "fatal:  error returned from Set_Pack_Obj_Fn()\n");
      return 0;
    }

    if (zz.Set_Unpack_Obj_Fn(migrate_unpack_elem,
                      (void *) mesh) == ZOLTAN_FATAL) {
      Gen_Error(0, "fatal:  error returned from Set_Unpack_Obj_Fn()\n");
      return 0;
    }
  }


  if (Test.Null_Lists == NONE) {
    if (zz.Migrate(num_imp, imp_gids, imp_lids, imp_procs, imp_to_part,
                   num_exp, exp_gids, exp_lids, exp_procs, exp_to_part)
                   == ZOLTAN_FATAL) {
      Gen_Error(0, "fatal:  error returned from Migrate()\n");
      return 0;
    }
  }
  else {
    /* Call zz.Help_Migrate with empty import lists. */
    /* Have to "manually" call migrate_pre_process and migrate_post_process. */
    int ierr = 0;
    migrate_pre_process((void *) mesh, 1, 1,
                        num_imp, imp_gids, imp_lids, imp_procs, imp_to_part,
                        num_exp, exp_gids, exp_lids, exp_procs, exp_to_part,
                        &ierr);
    if (Test.Null_Lists == IMPORT_LISTS) {
//.........这里部分代码省略.........

int read_mesh_params(const std::string &exo_file,
                     Problem_Description* problem,
                     Mesh_Description<INT>* mesh,
                     Sphere_Info* sphere
  )
{
  int    exoid, cpu_ws=0, io_ws=0;
  float  version;
  char   elem_type[MAX_STR_LENGTH+1];
/*---------------------------Execution Begins--------------------------------*/

  /* Open the ExodusII geometry file */
  int mode = EX_READ | problem->int64api;
  if((exoid=ex_open(exo_file.c_str(), mode, &cpu_ws, &io_ws, &version)) < 0)
  {
    Gen_Error(0, "fatal: unable to open ExodusII file for mesh params");
    return 0;
  }

  /* Get the init info */
  ex_init_params exo;
  if(ex_get_init_ext(exoid, &exo))
  {
    Gen_Error(0, "fatal: unable to get init info from ExodusII file");
    ex_close(exoid);
    return 0;
  }
  strcpy(mesh->title, exo.title);
  mesh->num_dims      = exo.num_dim;
  mesh->num_nodes     = exo.num_nodes;
  mesh->num_elems     = exo.num_elem;
  mesh->num_el_blks   = exo.num_elem_blk;
  mesh->num_node_sets = exo.num_node_sets;
  mesh->num_side_sets = exo.num_side_sets;

  /* Get the length of the concatenated node set node list */
  if(mesh->num_node_sets > 0)
  {
    mesh->ns_list_len = ex_inquire_int(exoid, EX_INQ_NS_NODE_LEN);
  }
  else
    mesh->ns_list_len = 0;

  /* Allocate and initialize memory for the sphere adjustment */
  sphere->adjust = (int*)malloc(sizeof(int)*3*(mesh->num_el_blks));
  if(!(sphere->adjust)) {
    Gen_Error(0, "fatal: insufficient memory");
    ex_close(exoid);
    return 0;
  }
  else {
    sphere->begin = sphere->adjust + mesh->num_el_blks;
    sphere->end   = sphere->begin  + mesh->num_el_blks;
    for(size_t cnt=0; cnt < mesh->num_el_blks; cnt++) {
      sphere->adjust[cnt] = 0;
      sphere->begin[cnt]  = 0;
      sphere->end[cnt]    = 0;
    }
  }

  std::vector<INT> el_blk_ids(mesh->num_el_blks);

  /* Read the element block IDs */
  if(ex_get_elem_blk_ids(exoid, &el_blk_ids[0]) < 0) {
    Gen_Error(0, "fatal: unable to get element block IDs");
    ex_close(exoid);
    return 0;
  }

  /* Determine the maximum number of nodes per element */
  mesh->max_np_elem = 0;
  for(size_t cnt=0; cnt < mesh->num_el_blks; cnt++) {
    INT num_elems, idum;
    INT nodes_in_elem;
    
    if(ex_get_elem_block(exoid, el_blk_ids[cnt], elem_type,
                         &num_elems, &nodes_in_elem, &idum)