hecmw_dist_refine.c

Go to the documentation of this file.

/*****************************************************************************
 * Copyright (c) 2019 FrontISTR Commons
 * This software is released under the MIT License, see LICENSE.txt
 *****************************************************************************/
 
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>
#include "hecmw_struct.h"
#include "hecmw_util.h"
#include "hecmw_common_define.h"
#include "hecmw_etype.h"
#include "hecmw_bit_array.h"
#include "hecmw_set_int.h"
#include "hecmw_varray_int.h"
#include "hecmw_dist.h"
#include "hecmw_dist_alloc.h"
#include "hecmw_dist_free.h"
#include "hecmw_dist_refine.h"
 
#ifdef HECMW_WITH_REFINER
 
#include "rcapRefiner.h"
 
/*============================================================================*/
/*                                                                            */
/*  convert elem_node_item                                                    */
/*                                                                            */
/*============================================================================*/
static const int *get_enode_h2r(int etype, int *ierror) {
  /* tetra, prism, hexa: no need to convert */
  static const int enode_pyr_h2r[13] = {4,  0,  1, 2, 3, 9, 10,
                                        11, 12, 5, 6, 7, 8};
  static const int enode_rod2_h2r[3] = {0, 2, 1};
  const int *enode_h2r;
 
  *ierror = 0;
 
  if (HECMW_is_etype_link(etype)) return NULL;
  if (HECMW_is_etype_truss(etype)) return NULL;
 
  switch (etype) {
    case HECMW_ETYPE_ROD1:
    /* case HECMW_ETYPE_ROD31: */
    case HECMW_ETYPE_TRI1:
    case HECMW_ETYPE_TRI2:
    case HECMW_ETYPE_QUA1:
    case HECMW_ETYPE_QUA2:
    case HECMW_ETYPE_TET1:
    case HECMW_ETYPE_TET1_4:
    case HECMW_ETYPE_TET2:
    case HECMW_ETYPE_PRI1:
    case HECMW_ETYPE_PRI2:
    case HECMW_ETYPE_HEX1:
    case HECMW_ETYPE_HEX1_4:
    case HECMW_ETYPE_HEX2:
    case HECMW_ETYPE_BEM1:
    case HECMW_ETYPE_BEM3:
    case HECMW_ETYPE_SHT1:
    case HECMW_ETYPE_SHT2:
    case HECMW_ETYPE_SHT6:
    case HECMW_ETYPE_SHQ1:
    case HECMW_ETYPE_SHQ2:
    case HECMW_ETYPE_SHQ8:
    case HECMW_ETYPE_PTT1:
    case HECMW_ETYPE_PTT2:
    case HECMW_ETYPE_PTQ1:
    case HECMW_ETYPE_PTQ2:
      enode_h2r = NULL;
      break;
    case HECMW_ETYPE_PYR1:
    case HECMW_ETYPE_PYR2:
      enode_h2r = enode_pyr_h2r;
      break;
    case HECMW_ETYPE_ROD2:
    case HECMW_ETYPE_BEM2:
      enode_h2r = enode_rod2_h2r;
      break;
    default:
      HECMW_log(HECMW_LOG_ERROR,
                "Element type %d not supported for rerinement.\n", etype);
      *ierror   = 1;
      enode_h2r = NULL;
      break;
  }
  return enode_h2r;
}
 
static const int *get_enode_r2h(int etype, int *ierror) {
  /* tetra, prism, hexa: no need to convert */
  static const int enode_pyr_r2h[13] = {1,  2,  3, 4, 0, 9, 10,
                                        11, 12, 5, 6, 7, 8};
  static const int enode_rod2_r2h[3] = {0, 2, 1};
  const int *enode_r2h;
 
  *ierror = 0;
 
  if (HECMW_is_etype_link(etype)) return NULL;
  if (HECMW_is_etype_truss(etype)) return NULL;
 
  switch (etype) {
    case HECMW_ETYPE_ROD1:
    case HECMW_ETYPE_ROD31:
    case HECMW_ETYPE_TRI1:
    case HECMW_ETYPE_TRI2:
    case HECMW_ETYPE_QUA1:
    case HECMW_ETYPE_QUA2:
    case HECMW_ETYPE_TET1:
    case HECMW_ETYPE_TET1_4:
    case HECMW_ETYPE_TET2:
    case HECMW_ETYPE_PRI1:
    case HECMW_ETYPE_PRI2:
    case HECMW_ETYPE_HEX1:
    case HECMW_ETYPE_HEX1_4:
    case HECMW_ETYPE_HEX2:
    case HECMW_ETYPE_BEM1:
    case HECMW_ETYPE_BEM3:
    case HECMW_ETYPE_SHT1:
    case HECMW_ETYPE_SHT2:
    case HECMW_ETYPE_SHT6:
    case HECMW_ETYPE_SHQ1:
    case HECMW_ETYPE_SHQ2:
    case HECMW_ETYPE_SHQ8:
    case HECMW_ETYPE_PTT1:
    case HECMW_ETYPE_PTT2:
    case HECMW_ETYPE_PTQ1:
    case HECMW_ETYPE_PTQ2:
      enode_r2h = NULL;
      break;
    case HECMW_ETYPE_PYR1:
    case HECMW_ETYPE_PYR2:
      enode_r2h = enode_pyr_r2h;
      break;
    case HECMW_ETYPE_ROD2:
    case HECMW_ETYPE_BEM2:
      enode_r2h = enode_rod2_r2h;
      break;
    default:
      HECMW_log(HECMW_LOG_ERROR,
                "Element type %d not supported for refinement.\n", etype);
      *ierror   = 1;
      enode_r2h = NULL;
  }
  return enode_r2h;
}
 
static int elem_node_item_hecmw2rcap(struct hecmwST_local_mesh *mesh) {
  int i, j;
  int etype, ierror;
  const int *enode_h2r;
  int *enode;
  int nn, istart, iend;
  int temp_enode[HECMW_MAX_NODE_MAX];
 
  for (i = 0; i < mesh->n_elem_type; i++) {
    etype     = mesh->elem_type_item[i];
    enode_h2r = get_enode_h2r(etype, &ierror);
    if (ierror) return HECMW_ERROR;
    if (enode_h2r == NULL) continue;
    nn     = HECMW_get_max_node(etype);
    istart = mesh->elem_type_index[i];
    iend   = mesh->elem_type_index[i + 1];
    for (j = istart; j < iend; j++) {
      enode = mesh->elem_node_item + mesh->elem_node_index[i];
      for (j = 0; j < nn; j++) temp_enode[j] = enode[enode_h2r[j]];
      for (j = 0; j < nn; j++) enode[j] = temp_enode[j];
    }
  }
  return HECMW_SUCCESS;
}
 
static int elem_node_item_rcap2hecmw(struct hecmwST_local_mesh *mesh) {
  int i, j;
  int etype, ierror;
  const int *enode_r2h;
  int *enode;
  int nn, istart, iend;
  int temp_enode[HECMW_MAX_NODE_MAX];
 
  for (i = 0; i < mesh->n_elem_type; i++) {
    etype     = mesh->elem_type_item[i];
    enode_r2h = get_enode_r2h(etype, &ierror);
    if (ierror) return HECMW_ERROR;
    if (enode_r2h == NULL) continue;
    nn     = HECMW_get_max_node(etype);
    istart = mesh->elem_type_index[i];
    iend   = mesh->elem_type_index[i + 1];
    for (j = istart; j < iend; j++) {
      enode = mesh->elem_node_item + mesh->elem_node_index[i];
      for (j = 0; j < nn; j++) temp_enode[j] = enode[enode_r2h[j]];
      for (j = 0; j < nn; j++) enode[j] = temp_enode[j];
    }
  }
  return HECMW_SUCCESS;
}
 
/*============================================================================*/
/*                                                                            */
/*  convert surface IDs                                                       */
/*                                                                            */
/*============================================================================*/
static const int *get_sid_h2r(int etype, int *ierror) {
  /* HECMW to RCAP */
  /*
  static const int sid_tri_h2r[4] = {-1, 0, 1, 2};
  static const int sid_qua_h2r[5] = {-1, 3, 1, 0, 2};
  static const int sid_tet_h2r[5] = {-1, 0, 1, 2, 3};
  static const int sid_pri_h2r[6] = {-1, 3, 4, 2, 0, 1};
  static const int sid_hex_h2r[7] = {-1, 5, 3, 2, 4, 0, 1};
  static const int sid_pyr_h2r[6] = {-1, 3, 1, 0, 2, 4};
  */
  /* FSTR to RCAP*/
  static const int sid_tri_h2r[4] = {-1, 2, 0, 1};
  static const int sid_qua_h2r[5] = {-1, 0, 1, 2, 3};
  static const int sid_tet_h2r[5] = {-1, 3, 2, 0, 1};
  static const int sid_pri_h2r[6] = {-1, 0, 1, 2, 3, 4};
  static const int sid_hex_h2r[7] = {-1, 0, 1, 2, 3, 4, 5};
  static const int sid_pyr_h2r[6] = {-1, 3, 1,
                                     0,  2, 4}; /* TEMPORARY: same as HECMW */
 
  const int *sid_h2r;
  *ierror = 0;
  switch (etype) {
    case HECMW_ETYPE_TRI1:
    case HECMW_ETYPE_TRI2:
      sid_h2r = sid_tri_h2r;
      break;
    case HECMW_ETYPE_QUA1:
    case HECMW_ETYPE_QUA2:
      sid_h2r = sid_qua_h2r;
      break;
    case HECMW_ETYPE_TET1:
    case HECMW_ETYPE_TET1_4:
    case HECMW_ETYPE_TET2:
      sid_h2r = sid_tet_h2r;
      break;
    case HECMW_ETYPE_PRI1:
    case HECMW_ETYPE_PRI2:
      sid_h2r = sid_pri_h2r;
      break;
    case HECMW_ETYPE_HEX1:
    case HECMW_ETYPE_HEX1_4:
    case HECMW_ETYPE_HEX2:
      sid_h2r = sid_hex_h2r;
      break;
    case HECMW_ETYPE_PYR1:
    case HECMW_ETYPE_PYR2:
      sid_h2r = sid_pyr_h2r;
      break;
    case HECMW_ETYPE_SHT1:
    case HECMW_ETYPE_SHT2:
    case HECMW_ETYPE_SHT6:
    case HECMW_ETYPE_SHQ1:
    case HECMW_ETYPE_SHQ2:
    case HECMW_ETYPE_SHQ8:
    case HECMW_ETYPE_PTT1:
    case HECMW_ETYPE_PTT2:
    case HECMW_ETYPE_PTQ1:
    case HECMW_ETYPE_PTQ2:
      sid_h2r = NULL;
      break;
    default:
      HECMW_log(HECMW_LOG_ERROR,
                "Element type %d not supported for refinement.\n", etype);
      *ierror = 1;
      sid_h2r = NULL;
  }
  return sid_h2r;
}
 
static const int *get_sid_r2h(int etype, int *ierror) {
  /* RCAP to HECMW */
  /*
  static const int sid_tri_r2h[3] = {1, 2, 3};
  static const int sid_qua_r2h[4] = {3, 2, 4, 1};
  static const int sid_tet_r2h[4] = {1, 2, 3, 4};
  static const int sid_pri_r2h[5] = {4, 5, 3, 1, 2};
  static const int sid_hex_r2h[6] = {1, 2, 3, 4, 5, 6};
  static const int sid_pyr_r2h[5] = {3, 2, 4, 1, 5};
  */
  /* RCAP to FSTR */
  static const int sid_tri_r2h[3] = {2, 3, 1};
  static const int sid_qua_r2h[4] = {1, 2, 3, 4};
  static const int sid_tet_r2h[4] = {3, 4, 2, 1};
  static const int sid_pri_r2h[5] = {1, 2, 3, 4, 5};
  static const int sid_hex_r2h[6] = {1, 2, 3, 4, 5, 6};
  static const int sid_pyr_r2h[5] = {3, 2, 4, 1,
                                     5}; /* TEMPORARY: same as HECMW */
 
  const int *sid_r2h;
  *ierror = 0;
  switch (etype) {
    case HECMW_ETYPE_TRI1:
    case HECMW_ETYPE_TRI2:
      sid_r2h = sid_tri_r2h;
      break;
    case HECMW_ETYPE_QUA1:
    case HECMW_ETYPE_QUA2:
      sid_r2h = sid_qua_r2h;
      break;
    case HECMW_ETYPE_TET1:
    case HECMW_ETYPE_TET1_4:
    case HECMW_ETYPE_TET2:
      sid_r2h = sid_tet_r2h;
      break;
    case HECMW_ETYPE_PRI1:
    case HECMW_ETYPE_PRI2:
      sid_r2h = sid_pri_r2h;
      break;
    case HECMW_ETYPE_HEX1:
    case HECMW_ETYPE_HEX1_4:
    case HECMW_ETYPE_HEX2:
      sid_r2h = sid_hex_r2h;
      break;
    case HECMW_ETYPE_PYR1:
    case HECMW_ETYPE_PYR2:
      sid_r2h = sid_pyr_r2h;
      break;
    case HECMW_ETYPE_SHT1:
    case HECMW_ETYPE_SHT2:
    case HECMW_ETYPE_SHT6:
    case HECMW_ETYPE_SHQ1:
    case HECMW_ETYPE_SHQ2:
    case HECMW_ETYPE_SHQ8:
    case HECMW_ETYPE_PTT1:
    case HECMW_ETYPE_PTT2:
    case HECMW_ETYPE_PTQ1:
    case HECMW_ETYPE_PTQ2:
      sid_r2h = NULL;
      break;
    default:
      HECMW_log(HECMW_LOG_ERROR,
                "Element type %d not supported for refinement.\n", etype);
      *ierror = 1;
      sid_r2h = NULL;
  }
  return sid_r2h;
}
 
static int surf_ID_hecmw2rcap(struct hecmwST_local_mesh *mesh) {
  struct hecmwST_surf_grp *grp = mesh->surf_group;
  int i;
 
  if (grp->n_grp == 0) return HECMW_SUCCESS;
 
  for (i = 0; i < grp->grp_index[grp->n_grp]; i++) {
    int eid    = grp->grp_item[2 * i];
    int *sid_p = &(grp->grp_item[2 * i + 1]);
    const int *sid_h2r;
    int ierror;
 
    HECMW_assert(0 < eid && eid <= mesh->n_elem_gross);
 
    sid_h2r = get_sid_h2r(mesh->elem_type[eid - 1], &ierror);
    if (ierror) return HECMW_ERROR;
    if (sid_h2r) *sid_p = sid_h2r[*sid_p];
  }
  return HECMW_SUCCESS;
}
 
static int surf_ID_rcap2hecmw(struct hecmwST_local_mesh *mesh) {
  struct hecmwST_surf_grp *grp = mesh->surf_group;
  int i;
 
  if (grp->n_grp == 0) return HECMW_SUCCESS;
 
  for (i = 0; i < grp->grp_index[grp->n_grp]; i++) {
    int eid    = grp->grp_item[2 * i];
    int *sid_p = &(grp->grp_item[2 * i + 1]);
    const int *sid_r2h;
    int ierror;
 
    HECMW_assert(0 < eid && eid <= mesh->n_elem_gross);
 
    sid_r2h = get_sid_r2h(mesh->elem_type[eid - 1], &ierror);
    if (ierror) return HECMW_ERROR;
    if (sid_r2h) *sid_p = sid_r2h[*sid_p];
  }
  return HECMW_SUCCESS;
}
 
/*============================================================================*/
/*                                                                            */
/* for models with multiple element types                                     */
/*                                                                            */
/*============================================================================*/
#define NDIV_SEG 2
#define NDIV_SURF 4
#define NDIV_BODY 8
#define NDIV_OTHER 1
#define NDIV_MAX 8
 
static int get_elem_ndiv(int etype, int *ierror) {
  int ndiv;
  *ierror = 0;
  if (HECMW_is_etype_truss(etype)) {
    ndiv = NDIV_OTHER;
  } else if (HECMW_is_etype_rod(etype) || HECMW_is_etype_beam(etype)) {
    ndiv = NDIV_SEG;
  } else if (HECMW_is_etype_surface(etype) || HECMW_is_etype_shell(etype)
             || HECMW_is_etype_patch(etype)) {
    ndiv = NDIV_SURF;
  } else if (HECMW_is_etype_solid(etype)) {
    ndiv = NDIV_BODY;
  } else if (HECMW_is_etype_link(etype)) {
    ndiv = NDIV_OTHER;
  } else {
    HECMW_log(HECMW_LOG_ERROR,
              "Element type %d not supported for refinement.\n", etype);
    *ierror = 1;
    ndiv    = NDIV_OTHER;
  }
  return ndiv;
}
 
/*============================================================================*/
/*                                                                            */
/*  prepare, clear and terminate REVOCAP_Refiner                              */
/*                                                                            */
/*============================================================================*/
static void register_node(struct hecmwST_local_mesh *mesh) {
  HECMW_log(HECMW_LOG_DEBUG, "Original Node Count = %d\n", mesh->n_node_gross);
 
#ifdef DEBUG_REFINER
  HECMW_log(HECMW_LOG_DEBUG, "RCAP: calling rcapSetNode64(num=%d, coords, globalIds, NULL)\n", mesh->n_node_gross);
  {
    char fname[256];
    FILE *fp_dump;
    snprintf(fname, sizeof(fname), "dump_rcapSetNode64.%d", HECMW_comm_get_rank());
    fp_dump = fopen(fname, "w");
    if (fp_dump == NULL) {
      fprintf(stderr, "Error: cannot open file %s\n", fname);
    } else {
      HECMW_log(HECMW_LOG_DEBUG, "RCAP: writing num, coords[0..3*num-1], globalIds[0..num-1] to %s\n", fname);
      fwrite(&(mesh->n_node_gross), sizeof(int), 1, fp_dump);
      fwrite(mesh->node, sizeof(double), 3*mesh->n_node_gross, fp_dump);
      fwrite(mesh->global_node_ID, sizeof(int), mesh->n_node_gross, fp_dump);
      fclose(fp_dump);
    }
  }
#endif
  rcapSetNode64(mesh->n_node_gross, mesh->node, mesh->global_node_ID, NULL);
 
  HECMW_assert(rcapGetNodeCount() == mesh->n_node_gross);
}
 
static void register_node_groups(struct hecmwST_node_grp *grp) {
  int i;
  char rcap_name[80];
 
  HECMW_assert(strcmp(grp->grp_name[0], "ALL") == 0);
 
  for (i = 1; i < grp->n_grp; i++) {
    int start  = grp->grp_index[i];
    int num    = grp->grp_index[i + 1] - start;
    int *array = grp->grp_item + start;
    snprintf(rcap_name, sizeof(rcap_name), "NG_%s", grp->grp_name[i]);
 
#ifdef DEBUG_REFINER
    HECMW_log(HECMW_LOG_DEBUG, "RCAP: calling rcapAppendNodeGroup(dataname=%s, num=%d, nodeArray)\n", rcap_name, num);
    {
      char fname[256];
      FILE *fp_dump;
      snprintf(fname, sizeof(fname), "dump_rcapAppendNodeGroup_%s.%d", grp->grp_name[i], HECMW_comm_get_rank());
      fp_dump = fopen(fname, "w");
      if (fp_dump == NULL) {
        fprintf(stderr, "Error: cannot open file %s\n", fname);
      } else {
        int len = strlen(rcap_name);
        HECMW_log(HECMW_LOG_DEBUG, "RCAP: writing strlen(dataname), dataname[], num, nodeArray[0..num-1] to %s\n", fname);
        fwrite(&len, sizeof(int), 1, fp_dump);
        fwrite(rcap_name, sizeof(char), len, fp_dump);
        fwrite(&num, sizeof(int), 1, fp_dump);
        fwrite(array, sizeof(int), num, fp_dump);
        fclose(fp_dump);
      }
    }
#endif
    rcapAppendNodeGroup(rcap_name, num, array);
  }
}
 
static int register_surf_groups(struct hecmwST_local_mesh *mesh) {
  struct hecmwST_surf_grp *grp = mesh->surf_group;
  int i, j;
  char rcap_name[80];
 
  struct hecmw_varray_int other;
 
  for (i = 0; i < grp->n_grp; i++) {
    int start  = grp->grp_index[i];
    int num    = grp->grp_index[i + 1] - start;
    int *array = grp->grp_item + start * 2;
 
    if (HECMW_varray_int_init(&other) != HECMW_SUCCESS) return HECMW_ERROR;
 
    for (j = 0; j < num; j++) {
      int elem  = array[j * 2];
      int surf  = array[j * 2 + 1];
      int etype = mesh->elem_type[elem - 1];
 
      /* ignore shell/patch surface */
      if (HECMW_is_etype_shell(etype) || HECMW_is_etype_patch(etype)) continue;
 
      if (HECMW_varray_int_append(&other, elem) != HECMW_SUCCESS)
        return HECMW_ERROR;
      if (HECMW_varray_int_append(&other, surf) != HECMW_SUCCESS)
        return HECMW_ERROR;
    }
 
    snprintf(rcap_name, sizeof(rcap_name), "SG_%s", grp->grp_name[i]);
    num   = HECMW_varray_int_nval(&other) / 2;
    array = HECMW_varray_int_get_v(&other);
 
#ifdef DEBUG_REFINER
    HECMW_log(HECMW_LOG_DEBUG, "RCAP: calling rcapAppendFaceGroup(dataname=%s, num=%d, faceArray)\n", rcap_name, num);
    {
      char fname[256];
      FILE *fp_dump;
      snprintf(fname, sizeof(fname), "dump_rcapAppendFaceGroup_%s.%d", grp->grp_name[i], HECMW_comm_get_rank());
      fp_dump = fopen(fname, "w");
      if (fp_dump == NULL) {
        fprintf(stderr, "Error: cannot open file %s\n", fname);
      } else {
        int len = strlen(rcap_name);
        HECMW_log(HECMW_LOG_DEBUG, "RCAP: writing strlen(dataname), dataname[], num, faceArray[0..2*num-1] to %s\n", fname);
        fwrite(&len, sizeof(int), 1, fp_dump);
        fwrite(rcap_name, sizeof(char), len, fp_dump);
        fwrite(&num, sizeof(int), 1, fp_dump);
        fwrite(array, sizeof(int), num*2, fp_dump);
        fclose(fp_dump);
      }
    }
#endif
    rcapAppendFaceGroup(rcap_name, num, array);
 
    HECMW_varray_int_finalize(&other);
  }
  return HECMW_SUCCESS;
}
 
static int prepare_refiner(struct hecmwST_local_mesh *mesh,
                           const char *cad_filename,
                           const char *part_filename) {
  HECMW_log(HECMW_LOG_DEBUG, "Started preparing refiner...\n");
 
  if (elem_node_item_hecmw2rcap(mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (surf_ID_hecmw2rcap(mesh) != HECMW_SUCCESS) return HECMW_ERROR;
 
#ifdef DEBUG_REFINER
  HECMW_log(HECMW_LOG_DEBUG, "RCAP: calling rcapInitRefiner(nodeOffset=1, elementOffset=1)\n");
#endif
  rcapInitRefiner(1, 1);
  if (cad_filename != NULL) {
#ifdef DEBUG_REFINER
    HECMW_log(HECMW_LOG_DEBUG, "RCAP: calling rcapSetCADFilename(filename=%s)\n", cad_filename);
#endif
    rcapSetCADFilename(cad_filename);
  }
  if (part_filename != NULL) {
#ifdef DEBUG_REFINER
    HECMW_log(HECMW_LOG_DEBUG, "RCAP: calling rcapSetPartitionFilename(filename=%s)\n", part_filename);
#endif
    rcapSetPartitionFilename(part_filename);
  }
  register_node(mesh);
 
  /* node groups are refined by REVOCAP_Refiner */
  register_node_groups(mesh->node_group);
 
  /* element groups are refined by myself */
 
  /* surface groups are refined by REVOCAP_Refiner except for shell/patch surface */
  if (register_surf_groups(mesh) != HECMW_SUCCESS) return HECMW_ERROR;
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished preparing refiner.\n");
  return HECMW_SUCCESS;
}
 
static void clear_refiner(void) {
#ifdef DEBUG_REFINER
  HECMW_log(HECMW_LOG_DEBUG, "RCAP: calling rcapClearRefiner()\n");
#endif
  rcapClearRefiner();
}
 
static int terminate_refiner(struct hecmwST_local_mesh *mesh) {
#ifdef DEBUG_REFINER
  HECMW_log(HECMW_LOG_DEBUG, "RCAP: calling rcapTermRefiner()\n");
#endif
  rcapTermRefiner();
  if (surf_ID_rcap2hecmw(mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (elem_node_item_rcap2hecmw(mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  return HECMW_SUCCESS;
}
 
/*============================================================================*/
/*                                                                            */
/*  call REVOCAP_Refiner                                                      */
/*                                                                            */
/*============================================================================*/
static int elem_type_hecmw2rcap(int etype) {
  switch (etype) {
    case HECMW_ETYPE_ROD1:
      return RCAP_SEGMENT;
    case HECMW_ETYPE_ROD2:
      return RCAP_SEGMENT2;
    case HECMW_ETYPE_TRI1:
      return RCAP_TRIANGLE;
    case HECMW_ETYPE_TRI2:
      return RCAP_TRIANGLE2;
    case HECMW_ETYPE_QUA1:
      return RCAP_QUAD;
    case HECMW_ETYPE_QUA2:
      return RCAP_QUAD2;
    case HECMW_ETYPE_ROD31:
      return RCAP_SEGMENT;
    case HECMW_ETYPE_TET1:
      return RCAP_TETRAHEDRON;
    case HECMW_ETYPE_TET1_4:
      return RCAP_TETRAHEDRON;
    case HECMW_ETYPE_TET2:
      return RCAP_TETRAHEDRON2;
    case HECMW_ETYPE_PRI1:
      return RCAP_WEDGE;
    case HECMW_ETYPE_PRI2:
      return RCAP_WEDGE2;
    case HECMW_ETYPE_PYR1:
      return RCAP_PYRAMID;
    case HECMW_ETYPE_PYR2:
      return RCAP_PYRAMID2;
    case HECMW_ETYPE_HEX1:
      return RCAP_HEXAHEDRON;
    case HECMW_ETYPE_HEX1_4:
      return RCAP_HEXAHEDRON;
    case HECMW_ETYPE_HEX2:
      return RCAP_HEXAHEDRON2;
    case HECMW_ETYPE_BEM1:
      return RCAP_SEGMENT;
    case HECMW_ETYPE_BEM2:
      return RCAP_SEGMENT2;
    case HECMW_ETYPE_BEM3:
      return RCAP_SEGMENT;
    case HECMW_ETYPE_SHT1:
      return RCAP_TRIANGLE;
    case HECMW_ETYPE_SHT2:
      return RCAP_TRIANGLE2;
    case HECMW_ETYPE_SHT6:
      return RCAP_TRIANGLE;
    case HECMW_ETYPE_SHQ1:
      return RCAP_QUAD;
    case HECMW_ETYPE_SHQ2:
      return RCAP_QUAD2;
    case HECMW_ETYPE_SHQ8:
      return RCAP_QUAD;
    case HECMW_ETYPE_PTT1:
      return RCAP_TRIANGLE;
    case HECMW_ETYPE_PTT2:
      return RCAP_TRIANGLE2;
    case HECMW_ETYPE_PTQ1:
      return RCAP_QUAD;
    case HECMW_ETYPE_PTQ2:
      return RCAP_QUAD2;
  }
  return RCAP_UNKNOWNTYPE;
}
 
static int reorder_enode_33struct(int etype, int n_elem,
                                  int *elem_node_item_ref) {
  int nn, nn_2, ndiv, nn_tot, ierror;
  int j, k, l;
  int tmp[32];
  int *enode1, *enode2;
  nn   = HECMW_get_max_node(etype);
  nn_2 = nn / 2;
  ndiv = get_elem_ndiv(etype, &ierror);
  if (ierror) return HECMW_ERROR;
  nn_tot = ndiv * nn;
  enode1 = elem_node_item_ref;
  enode2 = enode1 + nn_2 * ndiv;
  for (j = 0; j < n_elem; j++) {
    int *tmp_p = tmp;
    int *en1   = enode1;
    for (k = 0; k < ndiv; k++) {
      for (l = 0; l < nn_2; l++) {
        tmp_p[l]        = enode1[l];
        tmp_p[nn_2 + l] = enode2[l];
      }
      tmp_p += nn;
      enode1 += nn_2;
      enode2 += nn_2;
    }
    for (k = 0; k < nn_tot; k++) {
      en1[k] = tmp[k];
    }
    enode1 += nn_2 * ndiv;
    enode2 += nn_2 * ndiv;
  }
  return HECMW_SUCCESS;
}
 
static int refine_element(struct hecmwST_local_mesh *mesh,
                          struct hecmwST_local_mesh *ref_mesh) {
  int i, j;
  int n_elem_ref_tot;
  size_t n_enode_ref_tot;
  long long *elem_node_index_ref;
  int *elem_node_item_ref;
 
  HECMW_log(HECMW_LOG_DEBUG, "Original Element Count = %d\n", mesh->n_elem_gross);
 
  /*
   * count num of elems after refinement
   */
  n_elem_ref_tot  = 0;
  n_enode_ref_tot = 0;
  for (i = 0; i < mesh->n_elem_type; i++) {
    int etype, istart, iend, n_elem, etype_rcap, nn;
    size_t n_elem_ref;
    int *elem_node_item;
    etype = mesh->elem_type_item[i];
    istart = mesh->elem_type_index[i];
    iend   = mesh->elem_type_index[i + 1];
    n_elem = iend - istart;
    if (HECMW_is_etype_link(etype) || HECMW_is_etype_truss(etype)) {
      n_elem_ref = n_elem;
    } else {
      int is33 = HECMW_is_etype_33struct(etype);
      if (is33) n_elem *= 2;
      etype_rcap = elem_type_hecmw2rcap(etype);
      HECMW_assert(etype_rcap != RCAP_UNKNOWNTYPE);
      elem_node_item = mesh->elem_node_item + mesh->elem_node_index[istart];
 
#ifdef DEBUG_REFINER
      HECMW_log(HECMW_LOG_DEBUG, "RCAP: calling rcapRefineElement(num=%d, etype=%d, nodeArray, NULL)\n", n_elem, etype_rcap);
      {
        char fname[256];
        FILE *fp_dump;
        snprintf(fname, sizeof(fname), "dump_rcapRefineElement_pre_%d.%d", etype, HECMW_comm_get_rank());
        fp_dump = fopen(fname, "w");
        if (fp_dump == NULL) {
          fprintf(stderr, "Error: cannot open file %s\n", fname);
        } else {
          int nn = HECMW_get_max_node(etype);
          HECMW_log(HECMW_LOG_DEBUG, "RCAP: writing num, etype, nodeArray[0..%d*num-1] to %s\n", nn, fname);
          fwrite(&n_elem, sizeof(int), 1, fp_dump);
          fwrite(&etype_rcap, sizeof(int), 1, fp_dump);
          fwrite(elem_node_item, sizeof(int), nn * n_elem, fp_dump);
          fclose(fp_dump);
        }
      }
#endif
      n_elem_ref = rcapRefineElement(n_elem, etype_rcap, elem_node_item, NULL);
#ifdef DEBUG_REFINER
      HECMW_log(HECMW_LOG_DEBUG, "RCAP: return value was %d\n", n_elem_ref);
#endif
      if (is33) {
        n_elem_ref /= 2;
        n_elem /= 2;
      }
    }
#ifdef DEBUG
    {
      int ndiv, ierror;
      ndiv  = get_elem_ndiv(etype, &ierror);
      if (ierror) return HECMW_ERROR;
      HECMW_assert(n_elem_ref == n_elem * ndiv);
    }
#endif
    n_elem_ref_tot += n_elem_ref;
    nn = HECMW_get_max_node(etype);
    n_enode_ref_tot += n_elem_ref * nn;
  }
  ref_mesh->n_elem_gross = n_elem_ref_tot;
 
  HECMW_log(HECMW_LOG_DEBUG, "Refined Element Count will be %d\n", ref_mesh->n_elem_gross);
 
  /*
   * allocate memory
   */
  ref_mesh->elem_node_index =
      (long long *)HECMW_malloc(sizeof(long long) * (n_elem_ref_tot + 1));
  if (ref_mesh->elem_node_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  ref_mesh->elem_node_item = (int *)HECMW_calloc(n_enode_ref_tot, sizeof(int));
  if (ref_mesh->elem_node_item == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  HECMW_assert(rcapGetNodeCount() == mesh->n_node_gross);
 
  /*
   * perform actual refinement
   */
  ref_mesh->elem_node_index[0] = 0;
  elem_node_index_ref          = ref_mesh->elem_node_index;
  elem_node_item_ref           = ref_mesh->elem_node_item;
  n_elem_ref_tot               = 0;
  for (i = 0; i < mesh->n_elem_type; i++) {
    int etype, istart, iend, n_elem, etype_rcap, n_elem_ref, nn;
    long long jstart;
    int *elem_node_item;
    etype  = mesh->elem_type_item[i];
    istart = mesh->elem_type_index[i];
    iend   = mesh->elem_type_index[i + 1];
    n_elem = iend - istart;
    /* rcapRefineElement reads the first items of both node arrays even when
     * there is no element, which can be past the end of the arrays */
    if (n_elem == 0) continue;
    nn     = HECMW_get_max_node(etype);
    if (HECMW_is_etype_link(etype) || HECMW_is_etype_truss(etype)) {
      n_elem_ref = n_elem;
      for (j = 0; j < n_elem * nn; j++) {
        /* elem_node_item_ref[j] = elem_node_item[j]; */
        jstart = mesh->elem_node_index[istart];
        /* printf("aa %d %d\n",jstart,j); */
        elem_node_item_ref[j] = mesh->elem_node_item[jstart + j];
      }
    } else {
      int is33 = HECMW_is_etype_33struct(etype);
      if (is33) n_elem *= 2;
      etype_rcap     = elem_type_hecmw2rcap(etype);
      elem_node_item = mesh->elem_node_item + mesh->elem_node_index[istart];
 
#ifdef DEBUG_REFINER
      HECMW_log(HECMW_LOG_DEBUG, "RCAP: calling rcapRefineElement(num=%d, etype=%d, nodeArray, resultNodeArray)\n", n_elem, etype_rcap);
#endif
      n_elem_ref     = rcapRefineElement(n_elem, etype_rcap, elem_node_item,
                                     elem_node_item_ref);
#ifdef DEBUG_REFINER
      {
        char fname[256];
        FILE *fp_dump;
        snprintf(fname, sizeof(fname), "dump_rcapRefineElement_%d.%d", etype, HECMW_comm_get_rank());
        fp_dump = fopen(fname, "w");
        if (fp_dump == NULL) {
          fprintf(stderr, "Error: cannot open file %s\n", fname);
        } else {
          int nn = HECMW_get_max_node(etype);
          HECMW_log(HECMW_LOG_DEBUG,
                    "RCAP: writing num, etype, nodeArray[0..%d*num-1], num_ref, resultNodeArray[0..%d*num_ref-1] to %s\n",
                    nn, nn, fname);
          fwrite(&n_elem, sizeof(int), 1, fp_dump);
          fwrite(&etype_rcap, sizeof(int), 1, fp_dump);
          fwrite(elem_node_item, sizeof(int), nn * n_elem, fp_dump);
          fwrite(&n_elem_ref, sizeof(int), 1, fp_dump);
          fwrite(elem_node_item_ref, sizeof(int), nn * n_elem_ref, fp_dump);
          fclose(fp_dump);
        }
      }
#endif
 
      if (is33) {
        n_elem /= 2;
        n_elem_ref /= 2;
        reorder_enode_33struct(etype, n_elem, elem_node_item_ref);
      }
    }
    for (j = 0; j < n_elem_ref; j++) {
      elem_node_index_ref[j + 1] = elem_node_index_ref[j] + nn;
      if (elem_node_index_ref[j + 1] < elem_node_index_ref[j]) {
        HECMW_log(HECMW_LOG_ERROR,
                  "Integer overflow detected while creating elem_node_index\n");
        return HECMW_ERROR;
      }
    }
    elem_node_index_ref += n_elem_ref;
    elem_node_item_ref += nn * n_elem_ref;
    n_elem_ref_tot += n_elem_ref;
  }
 
  HECMW_log(HECMW_LOG_DEBUG, "Refined Element Count = %d\n", n_elem_ref_tot);
 
  HECMW_assert(n_elem_ref_tot == ref_mesh->n_elem_gross);
 
#ifdef DEBUG_REFINER
  HECMW_log(HECMW_LOG_DEBUG, "RCAP: calling rcapCommit()\n");
#endif
  rcapCommit();
  return HECMW_SUCCESS;
}
 
static int refine_node(struct hecmwST_local_mesh *mesh,
                       struct hecmwST_local_mesh *ref_mesh) {
  /* n_node_gross */
#ifdef DEBUG_REFINER
  HECMW_log(HECMW_LOG_DEBUG, "RCAP: calling rcapGetNodeCount()\n");
#endif
  ref_mesh->n_node_gross = rcapGetNodeCount();
#ifdef DEBUG_REFINER
  HECMW_log(HECMW_LOG_DEBUG, "RCAP: return value was %d\n", ref_mesh->n_node_gross);
#endif
 
  HECMW_log(HECMW_LOG_DEBUG, "Refined Node Count = %d\n", ref_mesh->n_node_gross);
 
  HECMW_assert(ref_mesh->n_node_gross > mesh->n_node_gross);
 
  /* node */
  ref_mesh->node =
      (double *)HECMW_malloc(sizeof(double) * ref_mesh->n_node_gross * 3);
  if (ref_mesh->node == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
#ifdef DEBUG_REFINER
  HECMW_log(HECMW_LOG_DEBUG, "RCAP: calling rcapGetNodeSeq64(num=%d, initId=1, coords)\n", ref_mesh->n_node_gross);
#endif
  rcapGetNodeSeq64(ref_mesh->n_node_gross, 1, ref_mesh->node);
#ifdef DEBUG_REFINER
  {
    char fname[256];
    FILE *fp_dump;
    snprintf(fname, sizeof(fname), "dump_rcapGetNodeSeq64.%d", HECMW_comm_get_rank());
    fp_dump = fopen(fname, "w");
    if (fp_dump == NULL) {
      fprintf(stderr, "Error: cannot open file %s\n", fname);
    } else {
      int one = 1;
      HECMW_log(HECMW_LOG_DEBUG, "RCAP: writing num, initId, coords[0..3*num-1] to %s\n", fname);
      fwrite(&(ref_mesh->n_node_gross), sizeof(int), 1, fp_dump);
      fwrite(&one, sizeof(int), 1, fp_dump);
      fwrite(ref_mesh->node, sizeof(int), 3 * ref_mesh->n_node_gross, fp_dump);
      fclose(fp_dump);
    }
  }
#endif
 
  return HECMW_SUCCESS;
}
 
static int refine_node_group_info(struct hecmwST_node_grp *grp,
                                  struct hecmwST_node_grp *ref_grp,
                                  int ref_n_node_gross) {
  int i;
  char rcap_name[80];
 
  ref_grp->n_grp = grp->n_grp;
 
  /* grp_name (COPY) */
  ref_grp->grp_name = (char **)HECMW_malloc(sizeof(char *) * grp->n_grp);
  if (ref_grp->grp_name == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  /* ALL */
  ref_grp->grp_name[0] = HECMW_strdup("ALL");
  if (ref_grp->grp_name[0] == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  /* other groups */
  for (i = 1; i < grp->n_grp; i++) {
    ref_grp->grp_name[i] = HECMW_strdup(grp->grp_name[i]);
    if (ref_grp->grp_name[i] == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
  }
 
  /* grp_index */
  ref_grp->grp_index = (int *)HECMW_malloc(sizeof(int) * (grp->n_grp + 1));
  if (ref_grp->grp_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  ref_grp->grp_index[0] = 0;
  /* ALL */
  ref_grp->grp_index[1] = ref_n_node_gross;
  /* other groups */
  for (i = 1; i < grp->n_grp; i++) {
    snprintf(rcap_name, sizeof(rcap_name), "NG_%s", grp->grp_name[i]);
    ref_grp->grp_index[i + 1] =
        ref_grp->grp_index[i] + rcapGetNodeGroupCount(rcap_name);
  }
 
  /* grp_item */
  if (ref_grp->grp_index[grp->n_grp] > 0) {
    ref_grp->grp_item =
        (int *)HECMW_calloc(ref_grp->grp_index[grp->n_grp], sizeof(int));
    if (ref_grp->grp_item == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    /* ALL */
    for (i = 0; i < ref_n_node_gross; i++) {
      ref_grp->grp_item[i] = i + 1;
    }
    /* other groups */
    for (i = 1; i < grp->n_grp; i++) {
      int start = ref_grp->grp_index[i];
      int num   = ref_grp->grp_index[i + 1] - start;
      if (num > 0) {
        int *array = ref_grp->grp_item + start;
        snprintf(rcap_name, sizeof(rcap_name), "NG_%s", grp->grp_name[i]);
        rcapGetNodeGroup(rcap_name, num, array);
      }
    }
  } else {
    ref_grp->grp_item = NULL;
  }
  return HECMW_SUCCESS;
}
 
/*
 * static functions for refining element groups
 */
 
static int get_refined_element_count(const struct hecmwST_local_mesh *mesh,
                                     int n_elem, int *elem_list) {
  int i, eid, etype, ierror, ndiv;
  int n_elem_ref = 0;
  for (i = 0; i < n_elem; i++) {
    eid = elem_list[i];
    HECMW_assert(0 < eid && eid <= mesh->n_elem_gross);
    etype = mesh->elem_type[eid - 1];
    ndiv  = get_elem_ndiv(etype, &ierror);
    if (ierror) return -1;
    n_elem_ref += ndiv;
  }
  return n_elem_ref;
}
 
static int get_refined_elements(const struct hecmwST_local_mesh *mesh,
                                struct hecmwST_local_mesh *ref_mesh, int eid,
                                int *eid_ref) {
  int i, etype, ierror, ndiv;
  int etype_idx, eid_type;
  int offset;
 
  HECMW_assert(0 < eid);
  HECMW_assert(eid <= mesh->n_elem_gross);
  HECMW_assert(eid_ref);
 
  etype = mesh->elem_type[eid - 1];
  ndiv  = get_elem_ndiv(etype, &ierror);
  if (ierror) return -1;
  etype_idx = -1;
  eid_type  = -1;
  for (i = 0; i < mesh->n_elem_type; i++) {
    if (mesh->elem_type_item[i] == etype) {
      etype_idx = i;
      HECMW_assert(mesh->elem_type_index[i] < eid);
      HECMW_assert(eid <= mesh->elem_type_index[i + 1]);
      eid_type = eid - mesh->elem_type_index[i];
      break;
    }
  }
  HECMW_assert(etype_idx >= 0);
  HECMW_assert(eid_type > 0);
 
  offset = ref_mesh->elem_type_index[etype_idx] + (eid_type - 1) * ndiv + 1;
  for (i = 0; i < ndiv; i++) {
    eid_ref[i] = offset + i;
  }
  return ndiv;
}
 
static int get_refined_elem_list(const struct hecmwST_local_mesh *mesh,
                                 struct hecmwST_local_mesh *ref_mesh,
                                 int n_elem, int *elem_list,
                                 int *elem_list_ref) {
  int i, j;
  int eid_ref[NDIV_MAX];
  int *el_ref    = elem_list_ref;
  int n_elem_ref = 0;
  for (i = 0; i < n_elem; i++) {
    int ndiv = get_refined_elements(mesh, ref_mesh, elem_list[i], eid_ref);
    if (ndiv < 0) return -1;
    for (j = 0; j < ndiv; j++) {
      el_ref[j] = eid_ref[j];
    }
    el_ref += ndiv;
    n_elem_ref += ndiv;
  }
  return n_elem_ref;
}
 
/*
 * refinement of element groups is done by myself
 */
 
static int refine_elem_group_info(struct hecmwST_local_mesh *mesh,
                                  struct hecmwST_local_mesh *ref_mesh) {
  int i;
  struct hecmwST_elem_grp *grp     = mesh->elem_group;
  struct hecmwST_elem_grp *ref_grp = ref_mesh->elem_group;
  int ref_n_elem_gross             = ref_mesh->n_elem_gross;
 
  ref_grp->n_grp = grp->n_grp;
 
  /* grp_name (COPY) */
  ref_grp->grp_name = (char **)HECMW_malloc(sizeof(char *) * grp->n_grp);
  if (ref_grp->grp_name == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  /* ALL */
  HECMW_assert(strcmp(grp->grp_name[0], "ALL") == 0);
  ref_grp->grp_name[0] = HECMW_strdup("ALL");
  if (ref_grp->grp_name[0] == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  /* other groups */
  for (i = 1; i < grp->n_grp; i++) {
    ref_grp->grp_name[i] = HECMW_strdup(grp->grp_name[i]);
    if (ref_grp->grp_name[i] == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
  }
 
  /* grp_index */
  ref_grp->grp_index = (int *)HECMW_malloc(sizeof(int) * (grp->n_grp + 1));
  if (ref_grp->grp_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  ref_grp->grp_index[0] = 0;
  /* ALL */
  ref_grp->grp_index[1] = ref_n_elem_gross;
  /* other groups */
  for (i = 1; i < grp->n_grp; i++) {
    int start, n_elem, n_elem_ref;
    int *elem_list;
    start      = grp->grp_index[i];
    n_elem     = grp->grp_index[i + 1] - start;
    elem_list  = grp->grp_item + start;
    n_elem_ref = get_refined_element_count(mesh, n_elem, elem_list);
    if (n_elem_ref < 0) return HECMW_ERROR;
    ref_grp->grp_index[i + 1] = ref_grp->grp_index[i] + n_elem_ref;
  }
 
  /* grp_item */
  if (ref_grp->grp_index[grp->n_grp] > 0) {
    ref_grp->grp_item =
        (int *)HECMW_calloc(ref_grp->grp_index[grp->n_grp], sizeof(int));
    if (ref_grp->grp_item == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    /* ALL */
    for (i = 0; i < ref_n_elem_gross; i++) {
      ref_grp->grp_item[i] = i + 1;
    }
    /* other groups */
    for (i = 1; i < grp->n_grp; i++) {
      int start, start_ref, n_elem, n_elem_ref, ret;
      int *elem_list, *elem_list_ref;
      start_ref  = ref_grp->grp_index[i];
      n_elem_ref = ref_grp->grp_index[i + 1] - start_ref;
      HECMW_assert(n_elem_ref >= 0);
      if (n_elem_ref == 0) continue;
      start         = grp->grp_index[i];
      n_elem        = grp->grp_index[i + 1] - start;
      elem_list     = grp->grp_item + start;
      elem_list_ref = ref_grp->grp_item + start_ref;
      ret           = get_refined_elem_list(mesh, ref_mesh, n_elem, elem_list,
                                  elem_list_ref);
      HECMW_assert(ret == n_elem_ref);
    }
  } else {
    ref_grp->grp_item = NULL;
  }
  return HECMW_SUCCESS;
}
 
/*
 * refinement of surface groups except for shell/patch surface is done
 * by REVOCAP_Refiner
 * refinement of shell/patch surfaces is done by myself
 */
 
static int refine_surf_group_info(struct hecmwST_local_mesh *mesh,
                                  struct hecmwST_local_mesh *ref_mesh) {
  struct hecmwST_surf_grp *grp     = mesh->surf_group;
  struct hecmwST_surf_grp *ref_grp = ref_mesh->surf_group;
  int i, j;
  char rcap_name[80];
 
  ref_grp->n_grp = grp->n_grp;
 
  /* grp_name (COPY) */
  ref_grp->grp_name = (char **)HECMW_malloc(sizeof(char *) * grp->n_grp);
  if (ref_grp->grp_name == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < grp->n_grp; i++) {
    ref_grp->grp_name[i] = HECMW_strdup(grp->grp_name[i]);
    if (ref_grp->grp_name[i] == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
  }
 
  /* grp_index */
  ref_grp->grp_index = (int *)HECMW_malloc(sizeof(int) * (grp->n_grp + 1));
  if (ref_grp->grp_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  ref_grp->grp_index[0] = 0;
  for (i = 0; i < grp->n_grp; i++) {
    int start  = grp->grp_index[i];
    int num    = grp->grp_index[i + 1] - start;
    int *array = grp->grp_item + start * 2;
    int num_sh, num_other;
 
    /* count surfaces except for shell/patch */
    snprintf(rcap_name, sizeof(rcap_name), "SG_%s", grp->grp_name[i]);
    num_other = rcapGetFaceGroupCount(rcap_name);
 
    /* count shell/patch surfaces */
    num_sh = 0;
    for (j = 0; j < num; j++) {
      int elem  = array[j * 2];
      int etype = mesh->elem_type[elem - 1];
      if (!(HECMW_is_etype_shell(etype) || HECMW_is_etype_patch(etype))) continue;
      num_sh += 1;
    }
    num_sh *= NDIV_SURF;
 
    ref_grp->grp_index[i + 1] = ref_grp->grp_index[i] + num_other + num_sh;
  }
  HECMW_assert(ref_grp->grp_index[grp->n_grp] >= 0);
 
  /* grp_item */
  if (ref_grp->grp_index[grp->n_grp] == 0) {
    ref_grp->grp_item = NULL;
    return HECMW_SUCCESS;
  }
  ref_grp->grp_item =
      (int *)HECMW_calloc(ref_grp->grp_index[grp->n_grp] * 2, sizeof(int));
  if (ref_grp->grp_item == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < grp->n_grp; i++) {
    struct hecmw_varray_int sh1, sh2;
    int start, num_tot, n_elem, ret;
    int start_ref, num_ref, num_tot_ref, n_elem_ref;
    int *array, *array_ref, *elem_list, *elem_list_ref;
 
    start_ref   = ref_grp->grp_index[i];
    num_tot_ref = ref_grp->grp_index[i + 1] - start_ref;
    if (num_tot_ref == 0) continue;
    array_ref = ref_grp->grp_item + start_ref * 2;
 
    /* get surfaces from REVOCAP_Refiner */
    snprintf(rcap_name, sizeof(rcap_name), "SG_%s", grp->grp_name[i]);
    num_ref = rcapGetFaceGroupCount(rcap_name);
    rcapGetFaceGroup(rcap_name, num_ref, array_ref);
 
    num_tot_ref -= num_ref;
    if (num_tot_ref == 0) continue;
    array_ref += num_ref * 2;
 
    /*
     * make refined shell/patch surfaces by myself
     */
 
    if (HECMW_varray_int_init(&sh1) != HECMW_SUCCESS) return HECMW_ERROR;
    if (HECMW_varray_int_init(&sh2) != HECMW_SUCCESS) return HECMW_ERROR;
 
    /* collect shell/patch surface */
    start   = grp->grp_index[i];
    num_tot = grp->grp_index[i + 1] - start;
    array   = grp->grp_item + start * 2;
    for (j = 0; j < num_tot; j++) {
      int elem, etype, surf;
      elem  = array[j * 2];
      surf  = array[j * 2 + 1];
      etype = mesh->elem_type[elem - 1];
      if (!(HECMW_is_etype_shell(etype) || HECMW_is_etype_patch(etype))) continue;
      if (surf == 1) {
        if (HECMW_varray_int_append(&sh1, elem) != HECMW_SUCCESS)
          return HECMW_ERROR;
      } else {
        HECMW_assert(surf == 2);
        if (HECMW_varray_int_append(&sh2, elem) != HECMW_SUCCESS)
          return HECMW_ERROR;
      }
    }
 
    /* 1st surface of shells/patches */
    n_elem        = HECMW_varray_int_nval(&sh1);
    if (n_elem > 0) {
      elem_list     = HECMW_varray_int_get_v(&sh1);
      n_elem_ref    = n_elem * NDIV_SURF;
      elem_list_ref = (int *)HECMW_calloc(n_elem_ref, sizeof(int));
      if (elem_list_ref == NULL) {
        HECMW_set_error(errno, "");
        return HECMW_ERROR;
      }
      ret = get_refined_elem_list(mesh, ref_mesh, n_elem, elem_list, elem_list_ref);
      HECMW_assert(ret == n_elem_ref);
      for (j = 0; j < n_elem_ref; j++) {
        array_ref[j * 2]     = elem_list_ref[j];
        array_ref[j * 2 + 1] = 1;
      }
      HECMW_free(elem_list_ref);
 
      num_tot_ref -= n_elem_ref;
      array_ref += n_elem_ref * 2;
    }
 
    /* 2nd surface of shells */
    n_elem        = HECMW_varray_int_nval(&sh2);
    if (n_elem > 0) {
      elem_list     = HECMW_varray_int_get_v(&sh2);
      n_elem_ref    = n_elem * NDIV_SURF;
      elem_list_ref = (int *)HECMW_calloc(n_elem_ref, sizeof(int));
      if (elem_list_ref == NULL) {
        HECMW_set_error(errno, "");
        return HECMW_ERROR;
      }
      ret = get_refined_elem_list(mesh, ref_mesh, n_elem, elem_list, elem_list_ref);
      HECMW_assert(ret == n_elem_ref);
      for (j = 0; j < n_elem_ref; j++) {
        array_ref[j * 2]     = elem_list_ref[j];
        array_ref[j * 2 + 1] = 2;
      }
      HECMW_free(elem_list_ref);
 
      num_tot_ref -= n_elem_ref;
    }
 
    HECMW_varray_int_finalize(&sh1);
    HECMW_varray_int_finalize(&sh2);
 
    HECMW_assert(num_tot_ref == 0);
  }
  return HECMW_SUCCESS;
}
 
static int call_refiner(struct hecmwST_local_mesh *mesh,
                        struct hecmwST_local_mesh *ref_mesh) {
  HECMW_log(HECMW_LOG_DEBUG, "Started calling refiner...\n");
 
  if (refine_element(mesh, ref_mesh) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
  if (refine_node(mesh, ref_mesh) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
  if (refine_node_group_info(mesh->node_group, ref_mesh->node_group,
                             ref_mesh->n_node_gross) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
  if (refine_elem_group_info(mesh, ref_mesh) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
  if (refine_surf_group_info(mesh, ref_mesh) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
  HECMW_log(HECMW_LOG_DEBUG, "Finished calling refiner.\n");
  return HECMW_SUCCESS;
}
 
/*============================================================================*/
/*                                                                            */
/*  rebuild information                                                       */
/*                                                                            */
/*============================================================================*/
static int rebuild_elem_ID(const struct hecmwST_local_mesh *mesh,
                           struct hecmwST_local_mesh *ref_mesh) {
  int i, j, k;
  int count;
  int *elem_ID_org;
  int *elem_ID_ref;
 
  ref_mesh->elem_ID =
      (int *)HECMW_malloc(sizeof(int) * ref_mesh->n_elem_gross * 2);
  if (ref_mesh->elem_ID == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  count       = 0;
  elem_ID_org = mesh->elem_ID;
  elem_ID_ref = ref_mesh->elem_ID;
  for (i = 0; i < mesh->n_elem_type; i++) {
    int etype, istart, iend, n_elem, ierror, ndiv;
    etype  = mesh->elem_type_item[i];
    istart = mesh->elem_type_index[i];
    iend   = mesh->elem_type_index[i + 1];
    n_elem = iend - istart;
    ndiv   = get_elem_ndiv(etype, &ierror);
    if (ierror) return HECMW_ERROR;
 
    for (j = 0; j < n_elem; j++) {
      int rank = elem_ID_org[1];
      for (k = 0; k < ndiv; k++) {
        count++;
        elem_ID_ref[0] = -count; /* TEMPORARY */
        elem_ID_ref[1] = rank;   /* to be corrected later */
        elem_ID_ref += 2;
      }
      elem_ID_org += 2;
    }
  }
  return HECMW_SUCCESS;
}
 
/* new IDs are numbered downward from *min_ID; shift the start so that the IDs
 * of each rank do not overlap with those of the other ranks */
static int shift_new_global_ID(int *min_ID, int n_new) {
  HECMW_Comm comm = HECMW_comm_get_comm();
  int n_proc      = HECMW_comm_get_size();
  int my_rank     = HECMW_comm_get_rank();
  int min_ID_g, i;
  int *n_new_all;
 
  n_new_all = (int *)HECMW_malloc(sizeof(int) * n_proc);
  if (n_new_all == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  HECMW_Allreduce(min_ID, &min_ID_g, 1, HECMW_INT, HECMW_MIN, comm);
  HECMW_Allgather(&n_new, 1, HECMW_INT, n_new_all, 1, HECMW_INT, comm);
  *min_ID = min_ID_g;
  for (i = 0; i < my_rank; i++) *min_ID -= n_new_all[i];
  HECMW_free(n_new_all);
  return HECMW_SUCCESS;
}
 
static int rebuild_global_elem_ID(const struct hecmwST_local_mesh *mesh,
                                  struct hecmwST_local_mesh *ref_mesh) {
  int i, j, k;
  int min_ID;
  int *global_elem_ID_org;
  int *global_elem_ID_ref;
 
  ref_mesh->global_elem_ID =
      (int *)HECMW_malloc(sizeof(int) * ref_mesh->n_elem_gross);
  if (ref_mesh->global_elem_ID == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  min_ID = 0;
  global_elem_ID_org = mesh->global_elem_ID;
  for (i = 0; i < mesh->n_elem_gross; i++) {
    if (global_elem_ID_org[i] < min_ID)
      min_ID = global_elem_ID_org[i];
  }
  /* contact search tells master elements of different ranks apart by IDs */
  if (shift_new_global_ID(&min_ID, ref_mesh->n_elem_gross -
                                       mesh->n_elem_gross) != HECMW_SUCCESS)
    return HECMW_ERROR;
  global_elem_ID_ref = ref_mesh->global_elem_ID;
  for (i = 0; i < mesh->n_elem_type; i++) {
    int etype, istart, iend, n_elem, ierror, ndiv;
    etype  = mesh->elem_type_item[i];
    istart = mesh->elem_type_index[i];
    iend   = mesh->elem_type_index[i + 1];
    n_elem = iend - istart;
    ndiv   = get_elem_ndiv(etype, &ierror);
    if (ierror) return HECMW_ERROR;
 
    for (j = 0; j < n_elem; j++) {
      global_elem_ID_ref[0] = global_elem_ID_org[0];
      for (k = 1; k < ndiv; k++) {
        min_ID--;
        global_elem_ID_ref[k] = min_ID;
      }
      global_elem_ID_org += 1;
      global_elem_ID_ref += ndiv;
    }
  }
  return HECMW_SUCCESS;
}
 
static int rebuild_elem_type(const struct hecmwST_local_mesh *mesh,
                             struct hecmwST_local_mesh *ref_mesh) {
  int i, j, k;
  int *elem_type_ref;
 
  ref_mesh->elem_type =
      (int *)HECMW_malloc(sizeof(int) * ref_mesh->n_elem_gross);
  if (ref_mesh->elem_type == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  elem_type_ref = ref_mesh->elem_type;
  for (i = 0; i < mesh->n_elem_type; i++) {
    int etype, istart, iend, n_elem, ierror, ndiv;
    etype  = mesh->elem_type_item[i];
    istart = mesh->elem_type_index[i];
    iend   = mesh->elem_type_index[i + 1];
    n_elem = iend - istart;
    ndiv   = get_elem_ndiv(etype, &ierror);
    if (ierror) return HECMW_ERROR;
 
    for (j = 0; j < n_elem; j++) {
      for (k = 0; k < ndiv; k++) {
        elem_type_ref[k] = etype;
      }
      elem_type_ref += ndiv;
    }
  }
  return HECMW_SUCCESS;
}
 
static int rebuild_section_ID(const struct hecmwST_local_mesh *mesh,
                              struct hecmwST_local_mesh *ref_mesh) {
  int i, j, k;
  int *section_ID_org;
  int *section_ID_ref;
  ref_mesh->section_ID =
      (int *)HECMW_malloc(sizeof(int) * ref_mesh->n_elem_gross);
  if (ref_mesh->section_ID == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  section_ID_org = mesh->section_ID;
  section_ID_ref = ref_mesh->section_ID;
  for (i = 0; i < mesh->n_elem_type; i++) {
    int etype, istart, iend, n_elem, ierror, ndiv;
    etype  = mesh->elem_type_item[i];
    istart = mesh->elem_type_index[i];
    iend   = mesh->elem_type_index[i + 1];
    n_elem = iend - istart;
    ndiv   = get_elem_ndiv(etype, &ierror);
    if (ierror) return HECMW_ERROR;
 
    for (j = 0; j < n_elem; j++) {
      int sect_id = section_ID_org[0];
      for (k = 0; k < ndiv; k++) {
        section_ID_ref[k] = sect_id;
      }
      section_ID_org += 1;
      section_ID_ref += ndiv;
    }
  }
  return HECMW_SUCCESS;
}
 
static int rebuild_elem_mat_ID_index(const struct hecmwST_local_mesh *mesh,
                                     struct hecmwST_local_mesh *ref_mesh) {
  int i, j, k;
  int *elem_mat_ID_index_org;
  int *elem_mat_ID_index_ref;
 
  ref_mesh->elem_mat_ID_index =
      (int *)HECMW_malloc(sizeof(int) * (ref_mesh->n_elem_gross + 1));
  if (ref_mesh->elem_mat_ID_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  ref_mesh->elem_mat_ID_index[0] = 0;
  elem_mat_ID_index_org          = mesh->elem_mat_ID_index;
  elem_mat_ID_index_ref          = ref_mesh->elem_mat_ID_index;
  for (i = 0; i < mesh->n_elem_type; i++) {
    int etype, istart, iend, n_elem, ierror, ndiv;
    etype  = mesh->elem_type_item[i];
    istart = mesh->elem_type_index[i];
    iend   = mesh->elem_type_index[i + 1];
    n_elem = iend - istart;
    ndiv   = get_elem_ndiv(etype, &ierror);
    if (ierror) return HECMW_ERROR;
 
    for (j = 0; j < n_elem; j++) {
      int n = elem_mat_ID_index_org[1] - elem_mat_ID_index_org[0];
      for (k = 0; k < ndiv; k++) {
        elem_mat_ID_index_ref[k + 1] = elem_mat_ID_index_ref[k] + n;
      }
      elem_mat_ID_index_org += 1;
      elem_mat_ID_index_ref += ndiv;
    }
  }
  return HECMW_SUCCESS;
}
 
static int rebuild_elem_mat_ID_item(const struct hecmwST_local_mesh *mesh,
                                    struct hecmwST_local_mesh *ref_mesh) {
  int i, j, k, l;
  int *elem_mat_ID_index_org;
  int *elem_mat_ID_item_org;
  int *elem_mat_ID_item_ref;
 
  ref_mesh->elem_mat_ID_item = (int *)HECMW_malloc(
      sizeof(int) * ref_mesh->elem_mat_ID_index[ref_mesh->n_elem_gross]);
  if (ref_mesh->elem_mat_ID_item == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  elem_mat_ID_index_org = mesh->elem_mat_ID_index;
  elem_mat_ID_item_org  = mesh->elem_mat_ID_item;
  elem_mat_ID_item_ref  = ref_mesh->elem_mat_ID_item;
  for (i = 0; i < mesh->n_elem_type; i++) {
    int etype, istart, iend, n_elem, ierror, ndiv;
    etype  = mesh->elem_type_item[i];
    istart = mesh->elem_type_index[i];
    iend   = mesh->elem_type_index[i + 1];
    n_elem = iend - istart;
    ndiv   = get_elem_ndiv(etype, &ierror);
    if (ierror) return HECMW_ERROR;
 
    for (j = 0; j < n_elem; j++) {
      int n = elem_mat_ID_index_org[1] - elem_mat_ID_index_org[0];
      for (k = 0; k < ndiv; k++) {
        for (l = 0; l < n; l++) {
          elem_mat_ID_item_ref[l] = elem_mat_ID_item_org[l];
        }
        elem_mat_ID_item_ref += n;
      }
      elem_mat_ID_index_org += 1;
      elem_mat_ID_item_org += n;
    }
  }
  return HECMW_SUCCESS;
}
 
static int rebuild_elem_info(const struct hecmwST_local_mesh *mesh,
                             struct hecmwST_local_mesh *ref_mesh) {
  HECMW_log(HECMW_LOG_DEBUG, "Started rebuilding element info...\n");
 
  if (rebuild_elem_ID(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (rebuild_global_elem_ID(mesh, ref_mesh) != HECMW_SUCCESS)
    return HECMW_ERROR;
  if (rebuild_elem_type(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (rebuild_section_ID(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (rebuild_elem_mat_ID_index(mesh, ref_mesh) != HECMW_SUCCESS)
    return HECMW_ERROR;
  if (rebuild_elem_mat_ID_item(mesh, ref_mesh) != HECMW_SUCCESS)
    return HECMW_ERROR;
 
  ref_mesh->n_elem_mat_ID = ref_mesh->elem_mat_ID_index[ref_mesh->n_elem_gross];
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished rebuilding element info.\n");
  return HECMW_SUCCESS;
}
 
static int rebuild_node_info(const struct hecmwST_local_mesh *mesh,
                             struct hecmwST_local_mesh *ref_mesh) {
  int i;
  int count, min_ID;
 
  HECMW_log(HECMW_LOG_DEBUG, "Started rebuilding node info...\n");
 
  /* allocate node_ID */
  ref_mesh->node_ID =
      (int *)HECMW_malloc(sizeof(int) * ref_mesh->n_node_gross * 2);
  if (ref_mesh->node_ID == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  /* node_ID */
  for (i = 0; i < mesh->n_node_gross; i++) {
    ref_mesh->node_ID[2 * i]     = mesh->node_ID[2 * i];
    ref_mesh->node_ID[2 * i + 1] = mesh->node_ID[2 * i + 1];
  }
  count = 0;
  for (i = mesh->n_node_gross; i < ref_mesh->n_node_gross; i++) {
    count++;
    ref_mesh->node_ID[2 * i]     = -count;        /* TEMPORARY */
    ref_mesh->node_ID[2 * i + 1] = mesh->my_rank; /* to be corrected later */
  }
 
  /* global_node_ID */
  ref_mesh->global_node_ID =
      (int *)HECMW_malloc(sizeof(int) * ref_mesh->n_node_gross);
  if (ref_mesh->global_node_ID == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  min_ID = 0;
  for (i = 0; i < mesh->n_node_gross; i++) {
    ref_mesh->global_node_ID[i] = mesh->global_node_ID[i];
    if (mesh->global_node_ID[i] < min_ID)
      min_ID = mesh->global_node_ID[i];
  }
  for (i = mesh->n_node_gross; i < ref_mesh->n_node_gross; i++) {
    min_ID--;
    ref_mesh->global_node_ID[i] = min_ID;
  }
 
  if (mesh->hecmw_flag_initcon && mesh->n_node_gross > 0) {
    /* node_init_val_index */
    ref_mesh->node_init_val_index =
        (int *)HECMW_malloc(sizeof(int) * (ref_mesh->n_node_gross + 1));
    if (ref_mesh->node_init_val_index == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < mesh->n_node_gross + 1; i++) {
      ref_mesh->node_init_val_index[i] = mesh->node_init_val_index[i];
    }
    for (i = mesh->n_node_gross + 1; i < ref_mesh->n_node_gross + 1; i++) {
      ref_mesh->node_init_val_index[i] =
          mesh->node_init_val_index[mesh->n_node_gross];
    }
 
    /* node_init_val_item */
    if (ref_mesh->node_init_val_index[ref_mesh->n_node_gross] == 0) {
      ref_mesh->node_init_val_item = NULL;
    } else {
      ref_mesh->node_init_val_item = (double *)HECMW_malloc(
          sizeof(double) * mesh->node_init_val_index[mesh->n_node_gross]);
      if (ref_mesh->node_init_val_item == NULL) {
        HECMW_set_error(errno, "");
        return HECMW_ERROR;
      }
      for (i = 0; i < mesh->node_init_val_index[mesh->n_node_gross]; i++) {
        ref_mesh->node_init_val_item[i] = mesh->node_init_val_item[i];
      }
    }
  }
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished rebuilding node info.\n");
  return HECMW_SUCCESS;
}
 
/*
 * static functions called by rebuild_comm_tables
 */
 
static int get_refined_shared(const struct hecmwST_local_mesh *mesh,
                              struct hecmwST_local_mesh *ref_mesh,
                              struct hecmw_varray_int *shared) {
  int i;
 
  for (i = 0; i < mesh->n_neighbor_pe; i++) {
    int start, n_elem, n_elem_ref, ret;
    int *elem_list, *elem_list_ref;
    start      = mesh->shared_index[i];
    n_elem     = mesh->shared_index[i + 1] - start;
    elem_list  = mesh->shared_item + start;
    n_elem_ref = get_refined_element_count(mesh, n_elem, elem_list);
    if (n_elem_ref < 0) return HECMW_ERROR;
    HECMW_varray_int_resize(shared + i, n_elem_ref);
    elem_list_ref = HECMW_varray_int_get_v(shared + i);
    ret =
        get_refined_elem_list(mesh, ref_mesh, n_elem, elem_list, elem_list_ref);
    HECMW_assert(ret == n_elem_ref);
  }
  return HECMW_SUCCESS;
}
 
static int elem_type_rcap2hecmw(int etype) {
  switch (etype) {
    case RCAP_SEGMENT:
      return HECMW_ETYPE_ROD1;
    case RCAP_SEGMENT2:
      return HECMW_ETYPE_ROD2;
    case RCAP_TRIANGLE:
      return HECMW_ETYPE_TRI1;
    case RCAP_TRIANGLE2:
      return HECMW_ETYPE_TRI2;
    case RCAP_QUAD:
      return HECMW_ETYPE_QUA1;
    case RCAP_QUAD2:
      return HECMW_ETYPE_QUA2;
    case RCAP_TETRAHEDRON:
      return HECMW_ETYPE_TET1;
    case RCAP_TETRAHEDRON2:
      return HECMW_ETYPE_TET2;
    case RCAP_WEDGE:
      return HECMW_ETYPE_PRI1;
    case RCAP_WEDGE2:
      return HECMW_ETYPE_PRI2;
    case RCAP_PYRAMID:
      return HECMW_ETYPE_PYR1;
    case RCAP_PYRAMID2:
      return HECMW_ETYPE_PYR2;
    case RCAP_HEXAHEDRON:
      return HECMW_ETYPE_HEX1;
    case RCAP_HEXAHEDRON2:
      return HECMW_ETYPE_HEX2;
  }
  return HECMW_ERROR;
}
 
static int get_rcap_elem_max_node(int etype_rcap) {
  int etype_hecmw = elem_type_rcap2hecmw(etype_rcap);
  return HECMW_get_max_node(etype_hecmw);
}
 
static int determine_node_rank_of_purely_external_elems(
  const struct hecmwST_local_mesh *mesh,
  struct hecmwST_local_mesh *ref_mesh) {
  int eid, rank, k, nid;
  long long nn;
  int *elem_nodes;
 
  for (eid = 1; eid <= ref_mesh->n_elem_gross; eid++) {
    rank = ref_mesh->elem_ID[2 * (eid - 1) + 1];
 
    if (rank >= 0) continue; /* purely external elems have rank < 0 */
 
    nn = ref_mesh->elem_node_index[eid] - ref_mesh->elem_node_index[eid - 1];
    elem_nodes = ref_mesh->elem_node_item + ref_mesh->elem_node_index[eid - 1];
    for (k = 0; k < nn; k++) {
      nid = elem_nodes[k];
      if (nid <= mesh->n_node_gross) continue; /* skip old nodes */
 
      /* new nodes that have not been assigned new rank have default rank = my_rank */
      if (ref_mesh->node_ID[2 * (nid - 1) + 1] == ref_mesh->my_rank) {
        /* assign out-of-range rank = n_subdomain as purely external number = -rank-1 */
        ref_mesh->node_ID[2 * (nid - 1) + 1] = - ref_mesh->n_subdomain - 1;
      }
    }
  }
  return HECMW_SUCCESS;
}
 
static int determine_node_rank(const struct hecmwST_local_mesh *mesh,
                               const struct hecmw_varray_int *shared,
                               struct hecmwST_local_mesh *ref_mesh) {
  int i, j, k;
  struct hecmw_set_int boundary_nodes;
  int bnode;
 
  HECMW_log(HECMW_LOG_DEBUG, "Started determine_node_rank...\n");
 
  HECMW_set_int_init(&boundary_nodes);
 
  /* Collect nodes of the shared elements, except for the old nodes */
  for (i = 0; i < mesh->n_neighbor_pe; i++) {
    for (j = 0; j < HECMW_varray_int_nval(shared + i); j++) {
      int eid = HECMW_varray_int_get(shared + i, j);
      int nn =
          ref_mesh->elem_node_index[eid] - ref_mesh->elem_node_index[eid - 1];
      int *elem_nodes =
          ref_mesh->elem_node_item + ref_mesh->elem_node_index[eid - 1];
 
      for (k = 0; k < nn; k++) {
        if (elem_nodes[k] <= mesh->n_node_gross) continue;
 
        HECMW_set_int_add(&boundary_nodes, elem_nodes[k]);
      }
    }
  }
  HECMW_set_int_check_dup(&boundary_nodes);
 
  /* for each node: determine which domain it belongs to; save in node_ID */
  HECMW_set_int_iter_init(&boundary_nodes);
  while (HECMW_set_int_iter_next(&boundary_nodes, &bnode)) {
    int original[HECMW_MAX_NODE_MAX];
    int orig_type_rcap = rcapGetOriginal(bnode, original);
    int min_rank       = mesh->n_subdomain;
    int max_rank       = -1;
    int nn;
 
    HECMW_assert(orig_type_rcap != RCAP_UNKNOWNTYPE);
 
    nn = get_rcap_elem_max_node(orig_type_rcap);
 
    for (k = 0; k < nn; k++) {
      int rank;
 
      HECMW_assert(1 <= original[k] && original[k] < bnode);
 
      rank = ref_mesh->node_ID[2 * (original[k] - 1) + 1];
 
      /* for purely external nodes, rank (0, 1, ...) is stored as -rank-1 (-1,
       * -2, ...) */
      if (rank < 0) rank = -rank - 1;
 
      if (rank < min_rank) min_rank = rank;
      if (max_rank < rank) max_rank = rank;
    }
    HECMW_assert(0 <= min_rank && min_rank < mesh->n_subdomain);
    HECMW_assert(0 <= max_rank && max_rank < mesh->n_subdomain);
 
    /* min-rule for the first step, max-rule for the following steps */
    if (mesh->n_refine == 0) {
      ref_mesh->node_ID[2 * (bnode - 1) + 1] = min_rank;
    } else {
      ref_mesh->node_ID[2 * (bnode - 1) + 1] = max_rank;
    }
  }
 
  HECMW_set_int_finalize(&boundary_nodes);
 
  /* new nodes in purely external elems */
  if (determine_node_rank_of_purely_external_elems(mesh, ref_mesh) != HECMW_SUCCESS)
    return HECMW_ERROR;
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished determine_node_rank.\n");
  return HECMW_SUCCESS;
}
 
static int *conv_index_ucd2hec(int etype)
{
  static int conv_index_ucd2hec_rod1[] = {0, 1};
  static int conv_index_ucd2hec_rod2[] = {0, -1, 2};
  static int conv_index_ucd2hec_tri1[] = {0, 1, 2};
  static int conv_index_ucd2hec_tri2[] = {0, 1, 2, -1, -1, -1};
  static int conv_index_ucd2hec_qua1[] = {0, 1, 2, 3};
  static int conv_index_ucd2hec_qua2[] = {0, 1, 2, 3, -1, -1, -1, -1};
  static int conv_index_ucd2hec_tet1[] = {0, 3, 2, 1};
  static int conv_index_ucd2hec_tet2[] = {0, 3, 2, 1, -1, -1, -1, -1, -1, -1};
  static int conv_index_ucd2hec_pri1[] = {3, 4, 5, 0, 1, 2};
  static int conv_index_ucd2hec_pri2[] = {3,  4,  5,  0,  1,  2,  -1, -1,
                                          -1, -1, -1, -1, -1, -1, -1};
  static int conv_index_ucd2hec_hex1[] = {4, 5, 6, 7, 0, 1, 2, 3};
  static int conv_index_ucd2hec_hex2[] = {4,  5,  6,  7,  0,  1,  2,  3,  -1, -1,
                                          -1, -1, -1, -1, -1, -1, -1, -1, -1, -1};
  static int conv_index_ucd2hec_pyr1[] = {4, 0, 1, 2, 3};
  static int conv_index_ucd2hec_pyr2[] = {4,  0,  1,  2,  3,  -1, -1,
                                          -1, -1, -1, -1, -1, -1};
  static int conv_index_ucd2hec_mst1[] = {0, 1, 2, 3};
  static int conv_index_ucd2hec_mst2[] = {0, 1, 2, 3, -1, -1, -1, -1};
  static int conv_index_ucd2hec_msq1[] = {0, 1, 2, 3, 4};
  static int conv_index_ucd2hec_msq2[] = {0, 1, 2, 3, 4, -1, -1, -1, -1};
  static int conv_index_ucd2hec_jtt1[] = {3, 4, 5, 0, 1, 2};
  static int conv_index_ucd2hec_jtt2[] = {3,  4,  5,  0,  1,  2,
                                          -1, -1, -1, -1, -1, -1};
  static int conv_index_ucd2hec_jtq1[] = {4, 5, 6, 7, 0, 1, 2, 3};
  static int conv_index_ucd2hec_jtq2[] = {4,  5,  6,  7,  0,  1,  2,  3,
                                          -1, -1, -1, -1, -1, -1, -1, -1};
#if 0
  static int conv_index_ucd2hec_bem1[] = {0, 1};
  static int conv_index_ucd2hec_bem2[] = {0, -1, 1};
#endif
  static int conv_index_ucd2hec_sht1[] = {0, 1, 2};
  static int conv_index_ucd2hec_sht2[] = {0, 1, 2, -1, -1, -1};
  static int conv_index_ucd2hec_shq1[] = {0, 1, 2, 3};
  static int conv_index_ucd2hec_shq2[] = {0, 1, 2, 3, -1, -1, -1, -1};
  static int conv_index_ucd2hec_ln[] = {0, 1};
  switch (etype) {
  case HECMW_ETYPE_ROD1:
    return conv_index_ucd2hec_rod1;
  case HECMW_ETYPE_ROD2:
    return conv_index_ucd2hec_rod2;
  case HECMW_ETYPE_TRI1:
    return conv_index_ucd2hec_tri1;
  case HECMW_ETYPE_TRI2:
    return conv_index_ucd2hec_tri2;
  case HECMW_ETYPE_QUA1:
    return conv_index_ucd2hec_qua1;
  case HECMW_ETYPE_QUA2:
    return conv_index_ucd2hec_qua2;
  case HECMW_ETYPE_TET1:
    return conv_index_ucd2hec_tet1;
  case HECMW_ETYPE_TET2:
    return conv_index_ucd2hec_tet2;
  case HECMW_ETYPE_PRI1:
    return conv_index_ucd2hec_pri1;
  case HECMW_ETYPE_PRI2:
    return conv_index_ucd2hec_pri2;
  case HECMW_ETYPE_HEX1:
    return conv_index_ucd2hec_hex1;
  case HECMW_ETYPE_HEX2:
    return conv_index_ucd2hec_hex2;
  case HECMW_ETYPE_PYR1:
    return conv_index_ucd2hec_pyr1;
  case HECMW_ETYPE_PYR2:
    return conv_index_ucd2hec_pyr2;
  case HECMW_ETYPE_MST1:
    return conv_index_ucd2hec_mst1;
  case HECMW_ETYPE_MST2:
    return conv_index_ucd2hec_mst2;
  case HECMW_ETYPE_MSQ1:
    return conv_index_ucd2hec_msq1;
  case HECMW_ETYPE_MSQ2:
    return conv_index_ucd2hec_msq2;
  case HECMW_ETYPE_JTT1:
    return conv_index_ucd2hec_jtt1;
  case HECMW_ETYPE_JTT2:
    return conv_index_ucd2hec_jtt2;
  case HECMW_ETYPE_JTQ1:
    return conv_index_ucd2hec_jtq1;
  case HECMW_ETYPE_JTQ2:
    return conv_index_ucd2hec_jtq2;
  case HECMW_ETYPE_SHT1:
    return conv_index_ucd2hec_sht1;
  case HECMW_ETYPE_SHT2:
    return conv_index_ucd2hec_sht2;
  case HECMW_ETYPE_SHQ1:
    return conv_index_ucd2hec_shq1;
  case HECMW_ETYPE_SHQ2:
    return conv_index_ucd2hec_shq2;
  case HECMW_ETYPE_LN11:
  case HECMW_ETYPE_LN12:
  case HECMW_ETYPE_LN13:
  case HECMW_ETYPE_LN14:
  case HECMW_ETYPE_LN15:
  case HECMW_ETYPE_LN16:
  case HECMW_ETYPE_LN21:
  case HECMW_ETYPE_LN22:
  case HECMW_ETYPE_LN23:
  case HECMW_ETYPE_LN24:
  case HECMW_ETYPE_LN25:
  case HECMW_ETYPE_LN26:
  case HECMW_ETYPE_LN31:
  case HECMW_ETYPE_LN32:
  case HECMW_ETYPE_LN33:
  case HECMW_ETYPE_LN34:
  case HECMW_ETYPE_LN35:
  case HECMW_ETYPE_LN36:
  case HECMW_ETYPE_LN41:
  case HECMW_ETYPE_LN42:
  case HECMW_ETYPE_LN43:
  case HECMW_ETYPE_LN44:
  case HECMW_ETYPE_LN45:
  case HECMW_ETYPE_LN46:
  case HECMW_ETYPE_LN51:
  case HECMW_ETYPE_LN52:
  case HECMW_ETYPE_LN53:
  case HECMW_ETYPE_LN54:
  case HECMW_ETYPE_LN55:
  case HECMW_ETYPE_LN56:
  case HECMW_ETYPE_LN61:
  case HECMW_ETYPE_LN62:
  case HECMW_ETYPE_LN63:
  case HECMW_ETYPE_LN64:
  case HECMW_ETYPE_LN65:
  case HECMW_ETYPE_LN66:
    return conv_index_ucd2hec_ln;
  default:
    return NULL;
  }
}
 
static void write_elem_with_mismatching_refinement(const struct hecmwST_local_mesh *mesh,
                                                  const struct hecmwST_local_mesh *ref_mesh,
                                                  const struct hecmw_varray_int *shared,
                                                  const struct hecmw_varray_int *shared_nodes,
                                                  int i)
{
  int irank = mesh->neighbor_pe[i];
  int j, k, l, nid;
  int n_shared_elem = HECMW_varray_int_nval(shared + i);
  int n_shared_node = HECMW_varray_int_nval(shared_nodes + i);
  struct hecmw_varray_int shared_elem_node_index;
  struct hecmw_varray_int shared_elem_node_item;
  int n_shared_elem_node;
  int n_shared_elem_r;
  int n_shared_elem_node_r;
  int *shared_elem_node_index_l = NULL;
  int *shared_elem_node_item_l = NULL;
  int *shared_elem_node_index_r = NULL;
  int *shared_elem_node_item_r = NULL;
  { /* create elem_node_index/item for shared elements */
    if (HECMW_varray_int_init(&shared_elem_node_index) != HECMW_SUCCESS) goto exit;
    if (HECMW_varray_int_init(&shared_elem_node_item) != HECMW_SUCCESS) goto exit;
    if (HECMW_varray_int_append(&shared_elem_node_index, 0) != HECMW_SUCCESS) goto exit;
    n_shared_elem_node = 0;
    for (j = 0; j < n_shared_elem; j++) {
      int eid = HECMW_varray_int_get(shared + i, j);
      int *elem_nodes = ref_mesh->elem_node_item + ref_mesh->elem_node_index[eid-1];
      long long nn = ref_mesh->elem_node_index[eid] - ref_mesh->elem_node_index[eid-1];
      for (k = 0; k < nn; k++) {
        nid = elem_nodes[k];
        /* find local ID within shared_nodes */
        for (l = 0; l < n_shared_node; l++) {
          if (HECMW_varray_int_get(shared_nodes + i, l) == nid) break;
        }
        if (HECMW_varray_int_append(&shared_elem_node_item, l+1) != HECMW_SUCCESS) goto exit;
      }
      n_shared_elem_node += nn;
      if (HECMW_varray_int_append(&shared_elem_node_index, n_shared_elem_node) != HECMW_SUCCESS) goto exit;
    }
  }
  shared_elem_node_index_l = HECMW_varray_int_get_v(&shared_elem_node_index);
  shared_elem_node_item_l = HECMW_varray_int_get_v(&shared_elem_node_item);
  { /* send/receive shared_elem_node_index/item */
    HECMW_Comm comm;
    HECMW_Request requests[3];
    HECMW_Status statuses[3];
    int tag;
    { /* send elem_node_index, elem_node_item */
      int send_buf[2];
      comm = HECMW_comm_get_comm();
      send_buf[0] = n_shared_elem;
      send_buf[1] = n_shared_elem_node;
      tag = 9801;
      HECMW_Isend(send_buf, 2, HECMW_INT, irank, tag, comm, &requests[0]);
      tag = 9802;
      HECMW_Isend(shared_elem_node_index_l, n_shared_elem + 1, HECMW_INT, irank, tag, comm, &requests[1]);
      tag = 9803;
      HECMW_Isend(shared_elem_node_item_l, n_shared_elem_node, HECMW_INT, irank, tag, comm, &requests[2]);
    }
    { /* receive elem_node_index, elem_node_item */
      int recv_buf[2];
      /* receive n_shared_elem/n_shared_elem_node */
      tag = 9801;
      HECMW_Recv(recv_buf, 2, HECMW_INT, irank, tag, comm, &statuses[0]);
      n_shared_elem_r = recv_buf[0];
      n_shared_elem_node_r = recv_buf[1];
      /* allocate */
      shared_elem_node_index_r = (int *)HECMW_malloc(sizeof(int) * (n_shared_elem_r + 1));
      if (shared_elem_node_index_r == NULL) {
        HECMW_set_error(errno, "");
        goto exit;
      }
      shared_elem_node_item_r = (int *)HECMW_malloc(sizeof(int) * n_shared_elem_node_r);
      if (shared_elem_node_item_r == NULL) {
        HECMW_set_error(errno, "");
        goto exit;
      }
      /* receive index/item */
      tag = 9802;
      HECMW_Recv(shared_elem_node_index_r, n_shared_elem_r + 1, HECMW_INT, irank, tag, comm, &statuses[1]);
      tag = 9803;
      HECMW_Recv(shared_elem_node_item_r, n_shared_elem_node_r, HECMW_INT, irank, tag, comm, &statuses[2]);
    }
    HECMW_Waitall(3, requests, statuses);
  }
  { /* compare */
    if (n_shared_elem != n_shared_elem_r) {
      HECMW_log(HECMW_LOG_ERROR, "rank %d has different number of shared elements\n", irank);
      goto exit;
    }
    for (j = 0; j < n_shared_elem; j++) {
      int is_error = 0;
      if (shared_elem_node_index_l[j+1] != shared_elem_node_index_r[j+1]) {
        HECMW_log(HECMW_LOG_ERROR, "rank %d has different number of nodes for %dth shared element\n", irank, j+1);
        goto exit;
      }
      for (k = shared_elem_node_index_l[j]; k < shared_elem_node_index_l[j+1]; k++) {
        if (shared_elem_node_item_l[k] != shared_elem_node_item_r[k]) {
          HECMW_log(HECMW_LOG_ERROR, "rank %d has different node for %dth node of %dth shared element\n", irank,
                    k - shared_elem_node_index_l[j] + 1, j+1);
          is_error = 1;
          break;
        }
      }
      if (is_error) {
        int eid = HECMW_varray_int_get(shared + i, j);
        int *elem_nodes = ref_mesh->elem_node_item + ref_mesh->elem_node_index[eid-1];
        long long nn = ref_mesh->elem_node_index[eid] - ref_mesh->elem_node_index[eid-1];
        HECMW_log(HECMW_LOG_ERROR, "elem local ID: %d\n", eid);
        for (k = 0; k < nn; k++) {
          nid = elem_nodes[k];
          HECMW_log(HECMW_LOG_ERROR, "elem node (%2d) local ID: %10d, coordinates: %14.6e, %14.6e, %14.6e\n", k+1, nid,
                    ref_mesh->node[3*(nid-1)], ref_mesh->node[3*(nid-1)+1], ref_mesh->node[3*(nid-1)+2]);
        }
        /* TODO: original element */
        /* TODO: origin of the node */
        goto exit;
      }
    }
  }
 
exit:
  if (shared_elem_node_index_r != NULL) HECMW_free(shared_elem_node_index_r);
  if (shared_elem_node_item_r != NULL) HECMW_free(shared_elem_node_item_r);
  HECMW_varray_int_finalize(&shared_elem_node_index);
  HECMW_varray_int_finalize(&shared_elem_node_item);
}
 
static void write_refined_shared_in_ucd(const struct hecmwST_local_mesh *mesh,
                                       const struct hecmwST_local_mesh *ref_mesh,
                                       const struct hecmw_varray_int *shared,
                                       const struct hecmw_varray_int *shared_nodes,
                                       int i)
{
  char fname[64];
  FILE *fp;
  int irank = mesh->neighbor_pe[i];
  int j, k, l, nid;
  int n_shared_elem = HECMW_varray_int_nval(shared + i);
  int n_shared_node = HECMW_varray_int_nval(shared_nodes + i);
  snprintf(fname, sizeof(fname), "shared_%d_%d.inp", mesh->my_rank, irank);
  HECMW_log(HECMW_LOG_DEBUG, "writing refined shared elements to %s\n", fname);
  fp = fopen(fname, "w");
  if (fp == NULL) {
    HECMW_log(HECMW_LOG_ERROR, "failed to open file %s\n", fname);
  } else {
    fprintf(fp, "%d %d 1 1 0\n", n_shared_node, n_shared_elem);
    for (j = 0; j < n_shared_node; j++) {
      nid = HECMW_varray_int_get(shared_nodes + i, j);
      fprintf(fp, "%d %e %e %e\n", j+1, ref_mesh->node[3*(nid-1)], ref_mesh->node[3*(nid-1)+1], ref_mesh->node[3*(nid-1)+2]);
    }
    for (j = 0; j < n_shared_elem; j++) {
      int eid = HECMW_varray_int_get(shared + i, j);
      int etype = ref_mesh->elem_type[eid-1];
      int *elem_nodes = ref_mesh->elem_node_item + ref_mesh->elem_node_index[eid-1];
      long long nn = ref_mesh->elem_node_index[eid] - ref_mesh->elem_node_index[eid-1];
      int *index_ucd = conv_index_ucd2hec(etype);
      HECMW_assert(index_ucd);
      fprintf(fp, "%d 0 %s", j+1, HECMW_get_ucd_label(etype));
      for (k = 0; k < nn; k++) {
        if (index_ucd[k] < 0) continue;
        nid = elem_nodes[index_ucd[k]];
        for (l = 0; l < n_shared_node; l++) {
          if (HECMW_varray_int_get(shared_nodes + i, l) == nid) break;
        }
        fprintf(fp, " %d", l+1);
      }
      fprintf(fp, "\n");
    }
    fclose(fp);
  }
}
 
static int check_node_rank(const struct hecmwST_local_mesh *mesh,
                           const struct hecmw_varray_int *shared,
                           const struct hecmwST_local_mesh *ref_mesh)
{
  int *n_shared_nodes;
  struct hecmw_varray_int *shared_nodes, *shared_ranks;
  HECMW_Request *requests;
  HECMW_Status *statuses;
  HECMW_Comm comm;
  int *send_buf;
  int recv_buf[4];
  int n_error = 0, n_error_g;
  int i, j, k, irank, tag, nid, rank;
 
  HECMW_log(HECMW_LOG_DEBUG, "Started check_node_rank...\n");
 
  /* a rank without neighbors has nothing to check, but has to join the
   * agreements of the other ranks */
  if (mesh->n_neighbor_pe == 0) {
    comm = HECMW_comm_get_comm();
    HECMW_Allreduce(&n_error, &n_error_g, 1, HECMW_INT, HECMW_SUM, comm);
    if (n_error_g > 0) return HECMW_ERROR;
    HECMW_Allreduce(&n_error, &n_error_g, 1, HECMW_INT, HECMW_SUM, comm);
    if (n_error_g > 0) return HECMW_ERROR;
    return HECMW_SUCCESS;
  }
 
  /* alloc memory */
  n_shared_nodes = (int *)HECMW_malloc(sizeof(int) * mesh->n_neighbor_pe);
  if (n_shared_nodes == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  shared_nodes = (struct hecmw_varray_int *)HECMW_malloc(
    sizeof(struct hecmw_varray_int) * mesh->n_neighbor_pe);
  if (shared_nodes == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  shared_ranks = (struct hecmw_varray_int *)HECMW_malloc(
    sizeof(struct hecmw_varray_int) * mesh->n_neighbor_pe);
  if (shared_ranks == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  send_buf = (int *)HECMW_malloc(sizeof(int) * 4 * mesh->n_neighbor_pe);
  if (send_buf == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  comm = HECMW_comm_get_comm();
 
  requests = (HECMW_Request *)HECMW_malloc(sizeof(HECMW_Request) * mesh->n_neighbor_pe);
  if (requests == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  statuses = (HECMW_Status *)HECMW_malloc(sizeof(HECMW_Status) * mesh->n_neighbor_pe);
  if (statuses == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  /*
   * check number of refined shared nodes
   */
 
  /* Collect nodes of the shared elements, and send size */
  for (i = 0; i < mesh->n_neighbor_pe; i++) {
    /* number of shared elements in original mesh */
    send_buf[4*i] = mesh->shared_index[i+1] - mesh->shared_index[i];
 
    /* number of refined shared elements */
    send_buf[4*i+1] = HECMW_varray_int_nval(shared + i);
 
    HECMW_varray_int_init(shared_nodes + i);
    for (j = 0; j < HECMW_varray_int_nval(shared + i); j++) {
      int eid = HECMW_varray_int_get(shared + i, j);
      long long nn = ref_mesh->elem_node_index[eid] - ref_mesh->elem_node_index[eid - 1];
      int *elem_nodes = ref_mesh->elem_node_item + ref_mesh->elem_node_index[eid - 1];
 
      for (k = 0; k < nn; k++) {
#if 0  /* keep old nodes for plotting */
        if (elem_nodes[k] <= mesh->n_node_gross) continue; /* skip old nodes */
#endif
 
        HECMW_varray_int_append(shared_nodes + i, elem_nodes[k]);
      }
    }
    /* number of shared nodes before removing duplication */
    send_buf[4*i+2] = HECMW_varray_int_nval(shared_nodes + i);
 
#if 1
    HECMW_varray_int_rmdup(shared_nodes + i);
#endif
    n_shared_nodes[i] = HECMW_varray_int_nval(shared_nodes + i);
 
    /* number of shared nodes after removing duplication */
    send_buf[4*i+3] = n_shared_nodes[i];
 
    /* send sizes */
    irank = mesh->neighbor_pe[i];
    tag = 9901;
    HECMW_Isend(send_buf + 4*i, 4, HECMW_INT, irank, tag, comm, requests + i);
  }
 
  /* receive size and check */
  for (i = 0; i < mesh->n_neighbor_pe; i++) {
    irank = mesh->neighbor_pe[i];
    tag = 9901;
    HECMW_Recv(recv_buf, 4, HECMW_INT, irank, tag, comm, statuses + i);
 
    if (recv_buf[0] != send_buf[4*i]) {
      HECMW_log(HECMW_LOG_ERROR,
                "number of shared elements %d not match with %d in rank %d\n",
                send_buf[4*i], recv_buf[0], irank);
      n_error++;
    }
    if (recv_buf[1] != send_buf[4*i+1]) {
      HECMW_log(HECMW_LOG_ERROR,
                "number of refined shared elements %d not match with %d in rank %d\n",
                send_buf[4*i+1], recv_buf[1], irank);
    }
    if (recv_buf[2] != send_buf[4*i+2]) {
      HECMW_log(HECMW_LOG_ERROR,
                "number of new shared nodes before rmdup %d not match with %d in rank %d\n",
                send_buf[4*i+2], recv_buf[2], irank);
      n_error++;
    }
    if (recv_buf[3] != send_buf[4*i+3]) {
      HECMW_log(HECMW_LOG_ERROR,
                "number of new shared nodes %d not match with %d in rank %d\n",
                send_buf[4*i+3], recv_buf[3], irank);
      n_error++;
      write_elem_with_mismatching_refinement(mesh, ref_mesh, shared, shared_nodes, i);
      write_refined_shared_in_ucd(mesh, ref_mesh, shared, shared_nodes, i);
    }
  }
  HECMW_Waitall(mesh->n_neighbor_pe, requests, statuses);
 
  HECMW_Allreduce(&n_error, &n_error_g, 1, HECMW_INT, HECMW_SUM, comm);
  if (n_error_g > 0) return HECMW_ERROR;
 
  /*
   * check rank of new shared nodes
   */
 
  n_error = 0;
 
  /* send rank of new shared nodes */
  for (i = 0; i < mesh->n_neighbor_pe; i++) {
    HECMW_varray_int_init(shared_ranks + i);
    /* HECMW_varray_int_resize(shared_ranks + i, n_shared_nodes[i]); */
    for (j = 0; j < n_shared_nodes[i]; j++) {
      nid = HECMW_varray_int_get(shared_nodes + i, j);
      rank = ref_mesh->node_ID[2 * (nid - 1) + 1];
      /* HECMW_varray_int_set(shared_ranks + i, j, rank); */
      HECMW_varray_int_append(shared_ranks + i, rank);
    }
 
    irank = mesh->neighbor_pe[i];
    tag = 9902;
    int *send_ranks = HECMW_varray_int_get_v(shared_ranks + i);
    HECMW_Isend(send_ranks, n_shared_nodes[i], HECMW_INT, irank, tag, comm, requests + i);
  }
 
  /* receive rank of new shared nodes */
  for (i = 0; i < mesh->n_neighbor_pe; i++) {
    int *recv_ranks = (int *)HECMW_malloc(sizeof(int) * n_shared_nodes[i]);
    if (recv_ranks == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
 
    irank = mesh->neighbor_pe[i];
    tag = 9902;
    HECMW_Recv(recv_ranks, n_shared_nodes[i], HECMW_INT, irank, tag, comm, statuses + i);
 
    for (j = 0; j < n_shared_nodes[i]; j++) {
      if (recv_ranks[j] != HECMW_varray_int_get(shared_ranks + i, j)) {
        HECMW_log(HECMW_LOG_ERROR,
                  "rank of new shared node %d: %d not match with %d in rank %d\n",
                  HECMW_varray_int_get(shared_nodes + i, j),
                  HECMW_varray_int_get(shared_ranks + i, j), recv_ranks[j], irank);
        n_error++;
      }
    }
    HECMW_free(recv_ranks);
  }
  HECMW_Waitall(mesh->n_neighbor_pe, requests, statuses);
 
  HECMW_Allreduce(&n_error, &n_error_g, 1, HECMW_INT, HECMW_SUM, comm);
  if (n_error_g > 0) return HECMW_ERROR;
 
  /* free memory */
  for (i = 0; i < mesh->n_neighbor_pe; i++) {
    HECMW_varray_int_finalize(shared_nodes + i);
    HECMW_varray_int_finalize(shared_ranks + i);
  }
  HECMW_free(shared_nodes);
  HECMW_free(shared_ranks);
  HECMW_free(n_shared_nodes);
  HECMW_free(send_buf);
  HECMW_free(requests);
  HECMW_free(statuses);
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished check_node_rank.\n");
 
  return HECMW_SUCCESS;
}
 
static int count_internal(int my_rank, int n, int *ID) {
  int i, count = 0;
  for (i = 0; i < n; i++) {
    if (ID[2 * i + 1] == my_rank) {
      ID[2 * i] = ++count;
    }
  }
  return count;
}
 
static int count_purely_external(int n, int *ID) {
  int i, count = 0;
  for (i = 0; i < n; i++) {
    if (ID[2 * i + 1] < 0) {
      count++;
    }
  }
  return count;
}
 
static int *new_internal_list(int my_rank, int n, int n_internal, int *ID) {
  int i, k;
  int *internal_list;
 
  internal_list = (int *)HECMW_malloc(sizeof(int) * n_internal);
  if (internal_list == NULL) {
    HECMW_set_error(errno, "");
    return NULL;
  }
 
  for (i = 0, k = 0; i < n; i++) {
    if (ID[2 * i + 1] == my_rank) {
      internal_list[k] = i + 1;
      k++;
    }
  }
  return internal_list;
}
 
static int create_index_item(int n, const struct hecmw_varray_int *arrays,
                             int **index, int **item) {
  int i, j, k;
 
  /* index */
  *index = (int *)HECMW_malloc(sizeof(int) * (n + 1));
  if (*index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  (*index)[0] = 0;
  for (i = 0; i < n; i++) {
    (*index)[i + 1] = (*index)[i] + HECMW_varray_int_nval(arrays + i);
    if ((*index)[i + 1] < (*index)[i]) {
      HECMW_log(HECMW_LOG_ERROR,
                "Integer overflow detected while creating index array\n");
      return HECMW_ERROR;
    }
  }
 
  /* item */
  if ((*index)[n] == 0) {
    (*item) = NULL;
    return HECMW_SUCCESS;
  }
  *item = (int *)HECMW_malloc(sizeof(int) * (*index)[n]);
  if (*item == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  k = 0;
  for (i = 0; i < n; i++) {
    for (j = 0; j < HECMW_varray_int_nval(arrays + i); j++) {
      (*item)[k] = HECMW_varray_int_get(arrays + i, j);
      k++;
    }
  }
  return HECMW_SUCCESS;
}
 
static int new_shared_export_import(const struct hecmwST_local_mesh *mesh,
                                    const struct hecmw_varray_int *shared,
                                    struct hecmwST_local_mesh *ref_mesh) {
  struct hecmw_varray_int *new_shared, *new_export, *new_import;
  int i;
 
  /* prepare new shared, import and export lists */
  new_shared = (struct hecmw_varray_int *)HECMW_malloc(
      sizeof(struct hecmw_varray_int) * mesh->n_neighbor_pe);
  if (new_shared == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  new_export = (struct hecmw_varray_int *)HECMW_malloc(
      sizeof(struct hecmw_varray_int) * mesh->n_neighbor_pe);
  if (new_export == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  new_import = (struct hecmw_varray_int *)HECMW_malloc(
      sizeof(struct hecmw_varray_int) * mesh->n_neighbor_pe);
  if (new_import == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  /* for each refined shared list */
  for (i = 0; i < mesh->n_neighbor_pe; i++) {
    int j;
    int rank_i = mesh->neighbor_pe[i];
 
    HECMW_varray_int_init(new_shared + i);
    HECMW_varray_int_init(new_export + i);
    HECMW_varray_int_init(new_import + i);
 
    /* for each element in the list */
    for (j = 0; j < HECMW_varray_int_nval(shared + i); j++) {
      struct hecmw_varray_int exp_cand, imp_cand;
 
      int eid = HECMW_varray_int_get(shared + i, j);
      int nn =
          ref_mesh->elem_node_index[eid] - ref_mesh->elem_node_index[eid - 1];
      int *elem_nodes =
          ref_mesh->elem_node_item + ref_mesh->elem_node_index[eid - 1];
      int etype = ref_mesh->elem_type[eid - 1];
 
      int min_rank        = mesh->n_subdomain;
      int myrank_included = 0;
      int rank_i_included = 0;
      int k;
 
      HECMW_varray_int_init(&exp_cand);
      HECMW_varray_int_init(&imp_cand);
 
      /* for each node in the element */
      for (k = 0; k < nn; k++) {
        int rank = ref_mesh->node_ID[2 * (elem_nodes[k] - 1) + 1];
 
        /* for purely external nodes, rank (0, 1, ...) is stored as -rank-1 (-1,
         * -2, ...) */
        if (rank < 0) rank = -rank - 1;
 
        if (rank < min_rank) min_rank = rank;
 
        if (rank == mesh->my_rank) {
          myrank_included = 1;
          HECMW_varray_int_append(&exp_cand, elem_nodes[k]);
        } else if (rank == rank_i) {
          rank_i_included = 1;
          HECMW_varray_int_append(&imp_cand, elem_nodes[k]);
        }
      }
      HECMW_assert(0 <= min_rank && min_rank < mesh->n_subdomain);
 
      ref_mesh->elem_ID[2 * (eid - 1) + 1] = min_rank;
 
      if (HECMW_is_etype_patch(etype)) {
        if (myrank_included || rank_i_included) HECMW_varray_int_append(new_shared + i, eid);
        if (myrank_included) HECMW_varray_int_cat(new_export + i, &exp_cand);
        if (rank_i_included) HECMW_varray_int_cat(new_import + i, &imp_cand);
      } else {
        if (myrank_included && rank_i_included) { /* the element is shared */
          HECMW_varray_int_append(new_shared + i, eid);
          HECMW_varray_int_cat(new_export + i, &exp_cand);
          HECMW_varray_int_cat(new_import + i, &imp_cand);
        } else if (myrank_included == 0) {
          /* for purely external elements, rank (0, 1, ...) is stored as -rank-1
           * (-1, -2, ...) */
          ref_mesh->elem_ID[2 * (eid - 1) + 1] = -min_rank - 1;
        }
      }
 
      HECMW_varray_int_finalize(&exp_cand);
      HECMW_varray_int_finalize(&imp_cand);
    }
    /* remove duplication */
    HECMW_varray_int_rmdup(new_export + i);
    HECMW_varray_int_rmdup(new_import + i);
  }
 
  /* new shared */
  if (create_index_item(mesh->n_neighbor_pe, new_shared, &(ref_mesh->shared_index),
                        &(ref_mesh->shared_item)) != HECMW_SUCCESS) {
    HECMW_log(HECMW_LOG_ERROR, "Create shared_index and shared_item failed\n");
    return HECMW_ERROR;
  }
  HECMW_log(HECMW_LOG_DEBUG, "Total number of shared elements = %d\n",
            ref_mesh->shared_index[mesh->n_neighbor_pe]);
 
  /* new export */
  if (create_index_item(mesh->n_neighbor_pe, new_export, &(ref_mesh->export_index),
                        &(ref_mesh->export_item)) != HECMW_SUCCESS) {
    HECMW_log(HECMW_LOG_ERROR, "Create export_index and export_item failed\n");
    return HECMW_ERROR;
  }
  HECMW_log(HECMW_LOG_DEBUG, "Total number of export nodes = %d\n",
            ref_mesh->export_index[mesh->n_neighbor_pe]);
 
  /* new import */
  if (create_index_item(mesh->n_neighbor_pe, new_import, &(ref_mesh->import_index),
                        &(ref_mesh->import_item)) != HECMW_SUCCESS) {
    HECMW_log(HECMW_LOG_ERROR, "Create import_index and import_item failed\n");
    return HECMW_ERROR;
  }
  HECMW_log(HECMW_LOG_DEBUG, "Total number of import nodes = %d\n",
            ref_mesh->import_index[mesh->n_neighbor_pe]);
 
  /* deallocate new shared, import and export lists */
  for (i = 0; i < mesh->n_neighbor_pe; i++) {
    HECMW_varray_int_finalize(new_shared + i);
    HECMW_varray_int_finalize(new_export + i);
    HECMW_varray_int_finalize(new_import + i);
  }
  HECMW_free(new_shared);
  HECMW_free(new_export);
  HECMW_free(new_import);
 
  return HECMW_SUCCESS;
}
 
static int rebuild_comm_tables_serial(const struct hecmwST_local_mesh *mesh,
                                      struct hecmwST_local_mesh *ref_mesh) {
  int i;
 
  HECMW_assert(mesh->hecmw_flag_parttype == HECMW_FLAG_PARTTYPE_UNKNOWN ||
               mesh->hecmw_flag_parttype == HECMW_FLAG_PARTTYPE_NODEBASED);
 
  /* nn_internal, n_node */
  ref_mesh->nn_internal = count_internal(mesh->my_rank, ref_mesh->n_node_gross, ref_mesh->node_ID);
  ref_mesh->n_node = ref_mesh->n_node_gross;
  ref_mesh->nn_middle = ref_mesh->n_node;
 
  HECMW_log(HECMW_LOG_DEBUG, "nn_internal = %d, n_node = %d, nn_middle = %d\n",
            ref_mesh->nn_internal, ref_mesh->n_node, ref_mesh->nn_middle);
 
  /* node_internal_list */
  ref_mesh->node_internal_list =
      (int *)HECMW_malloc(sizeof(int) * ref_mesh->nn_internal);
  if (ref_mesh->node_internal_list == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < ref_mesh->nn_internal; i++) {
    ref_mesh->node_internal_list[i] = i + 1;
  }
 
  /* ne_internal, n_elem */
  ref_mesh->ne_internal = count_internal(mesh->my_rank, ref_mesh->n_elem_gross, ref_mesh->elem_ID);
  ref_mesh->n_elem = ref_mesh->n_elem_gross;
 
  HECMW_log(HECMW_LOG_DEBUG, "ne_internal = %d, n_elem = %d\n",
            ref_mesh->ne_internal, ref_mesh->n_elem);
 
  /* elem_internal_list */
  ref_mesh->elem_internal_list =
      (int *)HECMW_malloc(sizeof(int) * ref_mesh->ne_internal);
  if (ref_mesh->elem_internal_list == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < ref_mesh->ne_internal; i++) {
    ref_mesh->elem_internal_list[i] = i + 1;
  }
  return HECMW_SUCCESS;
}
 
static int mark_purely_external_nodes(int my_rank,
                                      struct hecmwST_local_mesh *ref_mesh) {
  struct hecmw_bit_array mark;
  int i;
 
  if (HECMW_bit_array_init(&mark, ref_mesh->n_node_gross) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
  /* mark all nodes in the new import lists */
  for (i = 0; i < ref_mesh->import_index[ref_mesh->n_neighbor_pe]; i++) {
    HECMW_bit_array_set(&mark, ref_mesh->import_item[i] - 1);
  }
  /* scan for all external, unmarked nodes and set their rank to -rank-1 */
  for (i = 0; i < ref_mesh->n_node_gross; i++) {
    int rank = ref_mesh->node_ID[2 * i + 1];
    if (rank == my_rank || rank < 0) continue;
    if (!HECMW_bit_array_get(&mark, i)) {
      /* for purely external nodes, rank (0, 1, ...) is stored as -rank-1 (-1,
       * -2, ...) */
      ref_mesh->node_ID[2 * i + 1] = -rank - 1;
    }
  }
  HECMW_bit_array_finalize(&mark);
  return HECMW_SUCCESS;
}
 
static int rebuild_comm_tables(const struct hecmwST_local_mesh *mesh,
                               struct hecmwST_local_mesh *ref_mesh) {
  int i;
  struct hecmw_varray_int *shared;
 
  if (mesh->n_neighbor_pe == 0) { /* SERIAL */
    if (check_node_rank(mesh, NULL, ref_mesh) != HECMW_SUCCESS) {
      return HECMW_ERROR;
    }
    return rebuild_comm_tables_serial(mesh, ref_mesh);
  }
 
  /* PARALLEL */
  HECMW_log(HECMW_LOG_DEBUG, "Started rebuilding communication tables...\n");
 
  /* check part-type */
  if (mesh->hecmw_flag_parttype != HECMW_FLAG_PARTTYPE_NODEBASED) {
    HECMW_log(HECMW_LOG_ERROR,
              "Partitioning-type must be node-based for refinement.\n");
    return HECMW_ERROR;
  }
 
  /* Get refined shared element list */
  shared = (struct hecmw_varray_int *)HECMW_malloc(
      sizeof(struct hecmw_varray_int) * mesh->n_neighbor_pe);
  if (shared == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < mesh->n_neighbor_pe; i++) {
    HECMW_varray_int_init(shared + i);
  }
  if (get_refined_shared(mesh, ref_mesh, shared) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
 
  /* Determine rank of new nodes */
  if (determine_node_rank(mesh, shared, ref_mesh) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
 
  /* Check rank of new nodes with neighbor processes */
  if (check_node_rank(mesh, shared, ref_mesh) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
 
  /* nn_internal */
  ref_mesh->nn_internal =
      count_internal(mesh->my_rank, ref_mesh->n_node_gross, ref_mesh->node_ID);
 
  HECMW_assert(mesh->node_internal_list ==
               NULL); /* because PARTTYPE is NODEBASED. */
 
  /* new shared, import and export lists */
  if (new_shared_export_import(mesh, shared, ref_mesh) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
 
  /* deallocate refined shared list */
  for (i = 0; i < mesh->n_neighbor_pe; i++) {
    HECMW_varray_int_finalize(shared + i);
  }
  HECMW_free(shared);
 
  /* ne_internal */
  ref_mesh->ne_internal =
      count_internal(mesh->my_rank, ref_mesh->n_elem_gross, ref_mesh->elem_ID);
 
  /* n_elem */
  ref_mesh->n_elem =
      ref_mesh->n_elem_gross -
      count_purely_external(ref_mesh->n_elem_gross, ref_mesh->elem_ID);
 
  /* elem_internal_list */
  if (ref_mesh->ne_internal > 0) {
    ref_mesh->elem_internal_list =
        new_internal_list(mesh->my_rank, ref_mesh->n_elem_gross,
                          ref_mesh->ne_internal, ref_mesh->elem_ID);
    if (ref_mesh->elem_internal_list == NULL) return HECMW_ERROR;
  }
 
  /* mark purely external nodes */
  if (mark_purely_external_nodes(mesh->my_rank, ref_mesh) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
 
  /* n_node */
  ref_mesh->n_node =
      ref_mesh->n_node_gross -
      count_purely_external(ref_mesh->n_node_gross, ref_mesh->node_ID);
  ref_mesh->nn_middle = ref_mesh->n_node;  /* not set properly (since it's never used) */
 
  HECMW_log(HECMW_LOG_DEBUG, "nn_internal = %d, n_node = %d nn_middle = %d, \n",
            ref_mesh->nn_internal, ref_mesh->n_node, ref_mesh->nn_middle);
  HECMW_log(HECMW_LOG_DEBUG, "ne_internal = %d, nelem = %d\n",
            ref_mesh->ne_internal, ref_mesh->n_elem);
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished rebuilding communication tables.\n");
  return HECMW_SUCCESS;
}
 
static int check_comm_table_len(const struct hecmwST_local_mesh *ref_mesh) {
  int i, irank, js, je, len, tag, nsend;
  HECMW_Request *requests;
  HECMW_Status *statuses;
  HECMW_Comm comm;
  int *send_buf;
  int n_error = 0, n_error_g;
 
  /* a rank without neighbors has nothing to check, but has to join the
   * agreement of the other ranks */
  if (ref_mesh->n_neighbor_pe == 0) {
    comm = HECMW_comm_get_comm();
    HECMW_Allreduce(&n_error, &n_error_g, 1, HECMW_INT, HECMW_SUM, comm);
    if (n_error_g > 0) return HECMW_ERROR;
    return HECMW_SUCCESS;
  }
 
  HECMW_log(HECMW_LOG_DEBUG, "Started checking communication tables...\n");
 
  comm = HECMW_comm_get_comm();
 
  send_buf = (int *)HECMW_malloc(sizeof(int) * ref_mesh->n_neighbor_pe);
  if (send_buf == 0) {
    HECMW_set_error(errno, "");
    HECMW_abort(comm);
  }
 
  requests = (HECMW_Request *)HECMW_malloc(sizeof(HECMW_Request) *
                                           ref_mesh->n_neighbor_pe);
  if (requests == NULL) {
    HECMW_set_error(errno, "");
    HECMW_abort(comm);
  }
  statuses = (HECMW_Status *)HECMW_malloc(sizeof(HECMW_Status) *
                                          ref_mesh->n_neighbor_pe);
  if (statuses == NULL) {
    HECMW_set_error(errno, "");
    HECMW_abort(comm);
  }
 
  /* export and import */
 
  for (i = 0; i < ref_mesh->n_neighbor_pe; i++) {
    irank     = ref_mesh->neighbor_pe[i];
    js        = ref_mesh->export_index[i];
    je        = ref_mesh->export_index[i + 1];
    send_buf[i] = je - js;
    /* isend number of export nodes */
    tag = 0;
    HECMW_Isend(send_buf + i, 1, HECMW_INT, irank, tag, comm, requests + i);
  }
  for (i = 0; i < ref_mesh->n_neighbor_pe; i++) {
    irank     = ref_mesh->neighbor_pe[i];
    js        = ref_mesh->import_index[i];
    je        = ref_mesh->import_index[i + 1];
    /* recv number of import nodes */
    tag = 0;
    HECMW_Recv(&len, 1, HECMW_INT, irank, tag, comm, statuses + i);
    if (len != je - js) {
      HECMW_log(HECMW_LOG_ERROR,
                "inconsistent length of import (%d) with export on rank %d (%d)\n",
                je-js, irank, len);
      n_error++;
    }
  }
  HECMW_Waitall(ref_mesh->n_neighbor_pe, requests, statuses);
 
  /* shared */
 
  nsend = 0;
  for (i = 0; i < ref_mesh->n_neighbor_pe; i++) {
    irank   = ref_mesh->neighbor_pe[i];
    if (irank < ref_mesh->my_rank) continue;
    js      = ref_mesh->shared_index[i];
    je      = ref_mesh->shared_index[i + 1];
    send_buf[i] = je - js;
    /* isend number of shared elems */
    tag = 1;
    HECMW_Isend(send_buf + i, 1, HECMW_INT, irank, tag, comm, requests + nsend);
    nsend++;
  }
  for (i = 0; i < ref_mesh->n_neighbor_pe; i++) {
    irank   = ref_mesh->neighbor_pe[i];
    if (ref_mesh->my_rank < irank) continue;
    js      = ref_mesh->shared_index[i];
    je      = ref_mesh->shared_index[i + 1];
    /* recv number of shared elems */
    tag = 1;
    HECMW_Recv(&len, 1, HECMW_INT, irank, tag, comm, statuses);
    if (len != je - js) {
      HECMW_log(HECMW_LOG_ERROR,
                "inconsistent length of shared (%d) with rank %d (%d)\n",
                je-js, irank, len);
      n_error++;
    }
  }
  HECMW_Waitall(nsend, requests, statuses);
 
  HECMW_free(send_buf);
  HECMW_free(requests);
  HECMW_free(statuses);
 
  HECMW_Allreduce(&n_error, &n_error_g, 1, HECMW_INT, HECMW_SUM, comm);
  if (n_error_g > 0) return HECMW_ERROR;
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished checking communication tables.\n");
  return HECMW_SUCCESS;
}
 
static int rebuild_ID_external(struct hecmwST_local_mesh *ref_mesh) {
  int len_tot;
  int *sendbuf, *recvbuf, *srbuf;
  int i, j, irank, js, je, len, nid, tag, nsend, rank, lid;
  int *item_p, *sendbuf_p, *recvbuf_p, *srbuf_p;
  HECMW_Request *requests;
  HECMW_Status *statuses;
  HECMW_Comm comm;
 
  if (ref_mesh->n_neighbor_pe == 0) return HECMW_SUCCESS;
 
  HECMW_log(HECMW_LOG_DEBUG, "Started rebuilding external IDs...\n");
 
  comm = HECMW_comm_get_comm();
 
  requests = (HECMW_Request *)HECMW_malloc(sizeof(HECMW_Request) *
                                           ref_mesh->n_neighbor_pe);
  if (requests == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  statuses = (HECMW_Status *)HECMW_malloc(sizeof(HECMW_Status) *
                                          ref_mesh->n_neighbor_pe);
  if (statuses == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  /* node_ID (for external nodes) */
 
  /* send local IDs of export nodes */
  len_tot = ref_mesh->export_index[ref_mesh->n_neighbor_pe];
  sendbuf = (int *)HECMW_malloc(sizeof(int) * len_tot);
  if (sendbuf == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < ref_mesh->n_neighbor_pe; i++) {
    irank     = ref_mesh->neighbor_pe[i];
    js        = ref_mesh->export_index[i];
    je        = ref_mesh->export_index[i + 1];
    len       = je - js;
    item_p    = ref_mesh->export_item + js;
    sendbuf_p = sendbuf + js;
    for (j = 0; j < len; j++) {
      nid = item_p[j] - 1;
      lid = ref_mesh->node_ID[2 * nid];
      HECMW_assert(0 < lid);
      HECMW_assert(lid <= ref_mesh->nn_internal);
      sendbuf_p[j] = lid;
    }
    /* isend local id of export nodes */
    tag = 0;
    HECMW_Isend(sendbuf_p, len, HECMW_INT, irank, tag, comm, requests + i);
  }
 
  /* recv local IDs of import nodes */
  len_tot = ref_mesh->import_index[ref_mesh->n_neighbor_pe];
  recvbuf = (int *)HECMW_malloc(sizeof(int) * len_tot);
  if (recvbuf == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < ref_mesh->n_neighbor_pe; i++) {
    irank     = ref_mesh->neighbor_pe[i];
    js        = ref_mesh->import_index[i];
    je        = ref_mesh->import_index[i + 1];
    len       = je - js;
    item_p    = ref_mesh->import_item + js;
    recvbuf_p = recvbuf + js;
    /* recv local id of import nodes */
    tag = 0;
    HECMW_Recv(recvbuf_p, len, HECMW_INT, irank, tag, comm, statuses + i);
    /* set node_ID[2*j] */
    for (j = 0; j < len; j++) {
      nid = item_p[j] - 1;
      lid = recvbuf_p[j];
      HECMW_assert(0 < lid);
      ref_mesh->node_ID[2 * nid] = lid;
    }
  }
 
  HECMW_Waitall(ref_mesh->n_neighbor_pe, requests, statuses);
  HECMW_free(sendbuf);
  HECMW_free(recvbuf);
 
  /* elem_ID (for external elements) */
 
  len_tot = ref_mesh->shared_index[ref_mesh->n_neighbor_pe];
  srbuf   = (int *)HECMW_malloc(sizeof(int) * len_tot);
  if (srbuf == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  nsend = 0;
  for (i = 0; i < ref_mesh->n_neighbor_pe; i++) {
    irank   = ref_mesh->neighbor_pe[i];
    if (irank < ref_mesh->my_rank) continue;
    js      = ref_mesh->shared_index[i];
    je      = ref_mesh->shared_index[i + 1];
    len     = je - js;
    srbuf_p = srbuf + js;
    item_p  = ref_mesh->shared_item + js;
    for (j = 0; j < len; j++) {
      nid  = item_p[j] - 1;
      rank = ref_mesh->elem_ID[2 * nid + 1];
      lid  = ref_mesh->elem_ID[2 * nid];
      if (rank == ref_mesh->my_rank && lid > 0) {
        HECMW_assert(lid <= ref_mesh->ne_internal);
        srbuf_p[j] = lid;
      } else {
        srbuf_p[j] = -1;
      }
    }
    /* isend list of local_ID of those elems */
    tag = 1;
    HECMW_Isend(srbuf_p, len, HECMW_INT, irank, tag, comm, requests + nsend);
    nsend++;
  }
  for (i = 0; i < ref_mesh->n_neighbor_pe; i++) {
    irank   = ref_mesh->neighbor_pe[i];
    if (ref_mesh->my_rank < irank) continue;
    js      = ref_mesh->shared_index[i];
    je      = ref_mesh->shared_index[i + 1];
    len     = je - js;
    srbuf_p = srbuf + js;
    item_p  = ref_mesh->shared_item + js;
    /* recv list of local_ID of those elems */
    tag = 1;
    HECMW_Recv(srbuf_p, len, HECMW_INT, irank, tag, comm, statuses);
    /* set elem_ID[2*j] */
    for (j = 0; j < len; j++) {
      lid = srbuf_p[j];
      if (lid < 0) continue;
      nid  = item_p[j] - 1;
      rank = ref_mesh->elem_ID[2 * nid + 1];
      HECMW_assert(rank == irank || rank == -irank - 1);
      if (rank < 0) continue;
      ref_mesh->elem_ID[2 * nid] = lid;
    }
  }
  HECMW_Waitall(nsend, requests, statuses);
  HECMW_free(srbuf);
 
  HECMW_free(requests);
  HECMW_free(statuses);
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished rebuilding external IDs.\n");
  return HECMW_SUCCESS;
}
 
static int rebuild_refine_origin(const struct hecmwST_local_mesh *mesh,
                                 struct hecmwST_local_mesh *ref_mesh) {
  const struct hecmwST_refine_origin *reforg;
  struct hecmwST_refine_origin *ref_reforg;
  struct hecmw_varray_int ro_item;
  const int *ro_item_v;
  int n_refine, offset, cnt, cnt2, i, j;
 
  reforg   = mesh->refine_origin;
  n_refine = mesh->n_refine + 1;
 
  ref_reforg = (struct hecmwST_refine_origin *)HECMW_malloc(
      sizeof(struct hecmwST_refine_origin));
  if (ref_reforg == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  ref_reforg->index = (int *)HECMW_malloc(sizeof(int) * (n_refine + 1));
  if (ref_reforg->index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  ref_reforg->index[0] = 0;
  for (i = 1; i < n_refine; i++) {
    ref_reforg->index[i] = reforg->index[i];
  }
  ref_reforg->index[n_refine] =
      ref_reforg->index[n_refine - 1] + ref_mesh->n_node;
 
  ref_reforg->item_index =
      (int *)HECMW_malloc(sizeof(int) * (ref_reforg->index[n_refine] + 1));
  if (ref_reforg->item_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  /* copy original-node information up to previous refinement */
  ref_reforg->item_index[0] = 0;
  for (i = 1; i <= ref_reforg->index[n_refine - 1]; i++) {
    ref_reforg->item_index[i] = reforg->item_index[i];
  }
  offset = ref_reforg->index[n_refine - 1];
  cnt    = ref_reforg->item_index[offset];
  /* fprintf(stderr, "offset=%d, cnt=%d\n", offset, cnt); */
 
  /* get original nodes of newly generated nodes at current refinement */
  HECMW_varray_int_init(&ro_item);
  for (i = 0; i < ref_mesh->n_node; i++) {
    int iold;
    iold = ref_mesh->node_new2old ? ref_mesh->node_new2old[i] : i;
    if (iold < mesh->n_node_gross) {
      int org;
      org = mesh->node_old2new ? mesh->node_old2new[iold] + 1 : iold + 1;
      HECMW_varray_int_append(&ro_item, org);
      cnt++;
    } else {
      int original[HECMW_MAX_NODE_MAX];
      int orig_type_rcap = rcapGetOriginal(iold + 1, original);
      int nn;
      HECMW_assert(orig_type_rcap != RCAP_UNKNOWNTYPE);
      nn = get_rcap_elem_max_node(orig_type_rcap);
      for (j = 0; j < nn; j++) {
        int org                     = original[j];
        if (mesh->node_old2new) org = mesh->node_old2new[org - 1] + 1;
        HECMW_varray_int_append(&ro_item, org);
      }
      cnt += nn;
    }
    ref_reforg->item_index[offset + i + 1] = cnt;
  }
 
  ref_reforg->item_item = (int *)HECMW_malloc(sizeof(int) * cnt);
  if (ref_reforg->item_item == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  /* copy original nodes generated until previous refinements */
  cnt = ref_reforg->item_index[ref_reforg->index[n_refine - 1]];
  for (i = 0; i < cnt; i++) {
    ref_reforg->item_item[i] = reforg->item_item[i];
  }
  /* copy original nodes generated at current refinement */
  cnt2      = HECMW_varray_int_nval(&ro_item);
  ro_item_v = HECMW_varray_int_get_cv(&ro_item);
  for (i = 0; i < cnt2; i++) {
    ref_reforg->item_item[cnt + i] = ro_item_v[i];
  }
  HECMW_varray_int_finalize(&ro_item);
#if 0
    fprintf(stderr, "refine_origin->index:\n");
    for( i=0; i <= n_refine; i++ ) {
        fprintf(stderr, " %d", ref_reforg->index[i]);
        if( i % 10 == 9 ) fprintf(stderr, "\n");
    }
    if( i % 10 != 0 ) fprintf(stderr, "\n");
    fprintf(stderr, "refine_origin->item_index:\n");
    for( i=0; i <= ref_reforg->index[n_refine]; i++) {
        fprintf(stderr, " %d", ref_reforg->item_index[i]);
        if( i % 10 == 9 ) fprintf(stderr, "\n");
    }
    if( i % 10 != 0 ) fprintf(stderr, "\n");
    fprintf(stderr, "refine_origin->item_item:\n");
    for( i=0; i < ref_reforg->item_index[ref_reforg->index[n_refine]]; i++ ) {
        fprintf(stderr, " %d", ref_reforg->item_item[i]);
        if( i % 10 == 9 ) fprintf(stderr, "\n");
    }
    if( i % 10 != 0 ) fprintf(stderr, "\n");
#endif
  ref_mesh->refine_origin = ref_reforg;
  return HECMW_SUCCESS;
}
 
static int rebuild_n_node_refine_hist(const struct hecmwST_local_mesh *mesh,
                                      struct hecmwST_local_mesh *ref_mesh) {
  int i;
 
  ref_mesh->n_node_refine_hist =
    (int *)HECMW_malloc((mesh->n_refine + 1) * sizeof(int));
  if (ref_mesh->n_node_refine_hist == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  for (i = 0; i < mesh->n_refine; i++) {
    ref_mesh->n_node_refine_hist[i] = mesh->n_node_refine_hist[i];
  }
  ref_mesh->n_node_refine_hist[mesh->n_refine] = ref_mesh->n_node_gross;
 
  return HECMW_SUCCESS;
}
 
static int rebuild_refine_info(const struct hecmwST_local_mesh *mesh,
                               struct hecmwST_local_mesh *ref_mesh) {
  ref_mesh->n_refine = mesh->n_refine + 1;
 
  if (rebuild_refine_origin(mesh, ref_mesh) != HECMW_SUCCESS)
    return HECMW_ERROR;
  if (rebuild_n_node_refine_hist(mesh, ref_mesh) != HECMW_SUCCESS)
    return HECMW_ERROR;
 
  return HECMW_SUCCESS;
}
 
static int renumber_nodes_generate_tables(struct hecmwST_local_mesh *mesh);
static int renumber_elements_generate_tables(struct hecmwST_local_mesh *mesh);
 
static int rebuild_info(const struct hecmwST_local_mesh *mesh,
                        struct hecmwST_local_mesh *ref_mesh) {
  HECMW_log(HECMW_LOG_DEBUG, "Started rebuilding info...\n");
 
  if (rebuild_elem_info(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (rebuild_node_info(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (rebuild_comm_tables(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (check_comm_table_len(ref_mesh) != HECMW_SUCCESS) {
    HECMW_log(HECMW_LOG_ERROR, "Check communication table failed. Checking original mesh...\n");
    if (check_comm_table_len(mesh) != HECMW_SUCCESS) {
      HECMW_log(HECMW_LOG_ERROR, "Original mesh has error\n");
    } else {
      HECMW_log(HECMW_LOG_ERROR, "Original mesh is OK\n");
    }
    return HECMW_ERROR;
  }
  if (rebuild_ID_external(ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
 
  if (renumber_nodes_generate_tables(ref_mesh) != HECMW_SUCCESS)
    return HECMW_ERROR;
  if (renumber_elements_generate_tables(ref_mesh) != HECMW_SUCCESS)
    return HECMW_ERROR;
 
  if (rebuild_refine_info(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished rebuilding info.\n");
  return HECMW_SUCCESS;
}
 
/*============================================================================*/
/*                                                                            */
/*  copy unchanging information                                               */
/*                                                                            */
/*============================================================================*/
static int copy_global_info(const struct hecmwST_local_mesh *mesh,
                            struct hecmwST_local_mesh *ref_mesh) {
  int i;
 
  ref_mesh->hecmw_flag_adapt       = mesh->hecmw_flag_adapt;
  ref_mesh->hecmw_flag_initcon     = mesh->hecmw_flag_initcon;
  ref_mesh->hecmw_flag_parttype    = mesh->hecmw_flag_parttype;
  ref_mesh->hecmw_flag_partdepth   = mesh->hecmw_flag_partdepth;
  ref_mesh->hecmw_flag_version     = mesh->hecmw_flag_version;
  ref_mesh->hecmw_flag_partcontact = mesh->hecmw_flag_partcontact;
  snprintf(ref_mesh->gridfile, sizeof(ref_mesh->gridfile), "%s", mesh->gridfile);
  ref_mesh->hecmw_n_file       = mesh->hecmw_n_file;
 
  /* files */
  if (mesh->hecmw_n_file > 0) {
    ref_mesh->files = (char **)HECMW_calloc(mesh->hecmw_n_file, sizeof(char *));
    if (ref_mesh->files == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < mesh->hecmw_n_file; i++) {
      ref_mesh->files[i] = HECMW_strdup(mesh->files[i]);
      if (ref_mesh->files[i] == NULL) {
        HECMW_set_error(errno, "");
        return HECMW_ERROR;
      }
    }
  } else {
    ref_mesh->files = NULL;
  }
 
  snprintf(ref_mesh->header, sizeof(ref_mesh->header), "%s", mesh->header);
  ref_mesh->zero_temp = mesh->zero_temp;
 
  return HECMW_SUCCESS;
}
 
static int copy_node_ndof(const struct hecmwST_local_mesh *mesh,
                          struct hecmwST_local_mesh *ref_mesh) {
  int i;
 
  ref_mesh->n_dof     = mesh->n_dof;
  ref_mesh->n_dof_grp = mesh->n_dof_grp;
  ref_mesh->n_dof_tot = mesh->n_dof_tot;
 
  if (mesh->n_dof_grp > 0) {
    /* node_dof_index */
    ref_mesh->node_dof_index =
        (int *)HECMW_malloc(sizeof(int) * (mesh->n_dof_grp + 1));
    if (ref_mesh->node_dof_index == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < mesh->n_dof_grp + 1; i++) {
      ref_mesh->node_dof_index[i] = mesh->node_dof_index[i];
    }
 
    /* node_dof_item */
    ref_mesh->node_dof_item =
        (int *)HECMW_malloc(sizeof(int) * mesh->n_dof_grp);
    if (ref_mesh->node_dof_item == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < mesh->n_dof_grp; i++) {
      ref_mesh->node_dof_item[i] = mesh->node_dof_item[i];
    }
  }
  return HECMW_SUCCESS;
}
 
static int copy_elem_etype(const struct hecmwST_local_mesh *mesh,
                           struct hecmwST_local_mesh *ref_mesh) {
  int i;
 
  ref_mesh->n_elem_type = mesh->n_elem_type;
 
  if (mesh->n_elem_type > 0) {
    /* elem_type_index */
    ref_mesh->elem_type_index =
        (int *)HECMW_malloc(sizeof(int) * (mesh->n_elem_type + 1));
    if (ref_mesh->elem_type_index == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
 
    /* elem_type_item */
    ref_mesh->elem_type_item =
        (int *)HECMW_malloc(sizeof(int) * mesh->n_elem_type);
    if (ref_mesh->elem_type_item == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
 
    ref_mesh->elem_type_index[0] = 0;
    for (i = 0; i < mesh->n_elem_type; i++) {
      int etype, ierror, ndiv, nelem;
      etype = mesh->elem_type_item[i];
      ndiv  = get_elem_ndiv(etype, &ierror);
      if (ierror) return HECMW_ERROR;
      nelem = mesh->elem_type_index[i + 1] - mesh->elem_type_index[i];
      ref_mesh->elem_type_index[i + 1] =
          ref_mesh->elem_type_index[i] + nelem * ndiv;
      ref_mesh->elem_type_item[i] = etype;
    }
  }
  return HECMW_SUCCESS;
}
 
static int copy_comm_info(const struct hecmwST_local_mesh *mesh,
                          struct hecmwST_local_mesh *ref_mesh) {
  int i;
 
  ref_mesh->zero          = mesh->zero;
  ref_mesh->HECMW_COMM    = mesh->HECMW_COMM;
  ref_mesh->PETOT         = mesh->PETOT;
  ref_mesh->PEsmpTOT      = mesh->PEsmpTOT;
  ref_mesh->my_rank       = mesh->my_rank;
  ref_mesh->errnof        = mesh->errnof;
  ref_mesh->n_subdomain   = mesh->n_subdomain;
  ref_mesh->n_neighbor_pe = mesh->n_neighbor_pe;
 
  if (mesh->n_neighbor_pe == 0) {
    ref_mesh->neighbor_pe = NULL;
    ref_mesh->import_item = NULL;
    ref_mesh->export_item = NULL;
    ref_mesh->shared_item = NULL;
 
    ref_mesh->import_index = (int *)HECMW_malloc(sizeof(int));
    if (ref_mesh->import_index == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    ref_mesh->import_index[0] = 0;
 
    ref_mesh->export_index = (int *)HECMW_malloc(sizeof(int));
    if (ref_mesh->export_index == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    ref_mesh->export_index[0] = 0;
 
    ref_mesh->shared_index = (int *)HECMW_malloc(sizeof(int));
    if (ref_mesh->shared_index == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    ref_mesh->shared_index[0] = 0;
 
    return HECMW_SUCCESS;
  }
 
  /* neighbor_pe */
  ref_mesh->neighbor_pe =
      (int *)HECMW_malloc(sizeof(int) * mesh->n_neighbor_pe);
  if (ref_mesh->neighbor_pe == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < mesh->n_neighbor_pe; i++) {
    ref_mesh->neighbor_pe[i] = mesh->neighbor_pe[i];
  }
 
  return HECMW_SUCCESS;
}
 
static int copy_adapt_info(const struct hecmwST_local_mesh *mesh,
                           struct hecmwST_local_mesh *ref_mesh) {
  if (mesh->hecmw_flag_adapt != 0) {
    HECMW_log(HECMW_LOG_ERROR, "Refinement of adaptation not supported\n");
    return HECMW_ERROR;
  }
  ref_mesh->coarse_grid_level       = 0;
  ref_mesh->n_adapt                 = 0;
  ref_mesh->when_i_was_refined_node = NULL;
  ref_mesh->when_i_was_refined_elem = NULL;
  ref_mesh->adapt_parent_type       = NULL;
  ref_mesh->adapt_type              = NULL;
  ref_mesh->adapt_level             = NULL;
  ref_mesh->adapt_parent            = NULL;
  ref_mesh->adapt_children_index    = NULL;
  ref_mesh->adapt_children_item     = NULL;
  return HECMW_SUCCESS;
}
 
static int copy_section_info(const struct hecmwST_section *sect,
                             struct hecmwST_section *ref_sect) {
  int i;
 
  ref_sect->n_sect = sect->n_sect;
 
  if (sect->n_sect == 0) {
    ref_sect->sect_type         = NULL;
    ref_sect->sect_opt          = NULL;
    ref_sect->sect_mat_ID_index = NULL;
    ref_sect->sect_mat_ID_item  = NULL;
    ref_sect->sect_I_index      = NULL;
    ref_sect->sect_I_item       = NULL;
    ref_sect->sect_R_index      = NULL;
    ref_sect->sect_R_item       = NULL;
    return HECMW_SUCCESS;
  }
 
  /* sect_type */
  ref_sect->sect_type = (int *)HECMW_malloc(sizeof(int) * sect->n_sect);
  if (ref_sect->sect_type == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < sect->n_sect; i++) {
    ref_sect->sect_type[i] = sect->sect_type[i];
  }
 
  /* sect_opt */
  ref_sect->sect_opt = (int *)HECMW_malloc(sizeof(int) * sect->n_sect);
  if (ref_sect->sect_opt == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < sect->n_sect; i++) {
    ref_sect->sect_opt[i] = sect->sect_opt[i];
  }
 
  /* sect_mat_ID_index */
  ref_sect->sect_mat_ID_index =
      (int *)HECMW_malloc(sizeof(int) * (sect->n_sect + 1));
  if (ref_sect->sect_mat_ID_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < sect->n_sect + 1; i++) {
    ref_sect->sect_mat_ID_index[i] = sect->sect_mat_ID_index[i];
  }
 
  /* sect_mat_ID_item */
  if (sect->sect_mat_ID_index[sect->n_sect] > 0) {
    ref_sect->sect_mat_ID_item = (int *)HECMW_malloc(
        sizeof(int) * sect->sect_mat_ID_index[sect->n_sect]);
    if (ref_sect->sect_mat_ID_item == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < sect->sect_mat_ID_index[sect->n_sect]; i++) {
      ref_sect->sect_mat_ID_item[i] = sect->sect_mat_ID_item[i];
    }
  }
 
  /* sect_I_index */
  ref_sect->sect_I_index =
      (int *)HECMW_malloc(sizeof(int) * (sect->n_sect + 1));
  if (ref_sect->sect_I_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < sect->n_sect + 1; i++) {
    ref_sect->sect_I_index[i] = sect->sect_I_index[i];
  }
 
  /* sect_I_item */
  if (sect->sect_I_index[sect->n_sect] > 0) {
    ref_sect->sect_I_item =
        (int *)HECMW_malloc(sizeof(int) * sect->sect_I_index[sect->n_sect]);
    if (ref_sect->sect_I_item == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < sect->sect_I_index[sect->n_sect]; i++) {
      ref_sect->sect_I_item[i] = sect->sect_I_item[i];
    }
  }
 
  /* sect_R_index */
  ref_sect->sect_R_index =
      (int *)HECMW_malloc(sizeof(int) * (sect->n_sect + 1));
  if (ref_sect->sect_R_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < sect->n_sect + 1; i++) {
    ref_sect->sect_R_index[i] = sect->sect_R_index[i];
  }
 
  /* sect_R_item */
  if (sect->sect_R_index[sect->n_sect] > 0) {
    ref_sect->sect_R_item = (double *)HECMW_malloc(
        sizeof(double) * sect->sect_R_index[sect->n_sect]);
    if (ref_sect->sect_R_item == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < sect->sect_R_index[sect->n_sect]; i++) {
      ref_sect->sect_R_item[i] = sect->sect_R_item[i];
    }
  }
 
  return HECMW_SUCCESS;
}
 
static int copy_material_info(const struct hecmwST_material *mat,
                              struct hecmwST_material *ref_mat) {
  int i;
 
  ref_mat->n_mat = mat->n_mat;
 
  if (mat->n_mat == 0) {
    ref_mat->n_mat_item        = 0;
    ref_mat->n_mat_subitem     = 0;
    ref_mat->n_mat_table       = 0;
    ref_mat->mat_name          = NULL;
    ref_mat->mat_item_index    = NULL;
    ref_mat->mat_subitem_index = NULL;
    ref_mat->mat_table_index   = NULL;
    ref_mat->mat_val           = NULL;
    ref_mat->mat_temp          = NULL;
    return HECMW_SUCCESS;
  }
 
  ref_mat->n_mat_item    = mat->n_mat_item;
  ref_mat->n_mat_subitem = mat->n_mat_subitem;
  ref_mat->n_mat_table   = mat->n_mat_table;
 
  /* mat_name */
  ref_mat->mat_name = (char **)HECMW_malloc(sizeof(char *) * mat->n_mat);
  if (ref_mat->mat_name == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < mat->n_mat; i++) {
    ref_mat->mat_name[i] = HECMW_strdup(mat->mat_name[i]);
    if (ref_mat->mat_name[i] == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
  }
 
  /* mat_item_index */
  ref_mat->mat_item_index = (int *)HECMW_malloc(sizeof(int) * (mat->n_mat + 1));
  if (ref_mat->mat_item_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < mat->n_mat + 1; i++) {
    ref_mat->mat_item_index[i] = mat->mat_item_index[i];
  }
 
  /* mat_subitem_index */
  ref_mat->mat_subitem_index =
      (int *)HECMW_malloc(sizeof(int) * (mat->n_mat_item + 1));
  if (ref_mat->mat_subitem_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < mat->n_mat_item + 1; i++) {
    ref_mat->mat_subitem_index[i] = mat->mat_subitem_index[i];
  }
 
  /* mat_table_index */
  ref_mat->mat_table_index =
      (int *)HECMW_malloc(sizeof(int) * (mat->n_mat_subitem + 1));
  if (ref_mat->mat_table_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < mat->n_mat_subitem + 1; i++) {
    ref_mat->mat_table_index[i] = mat->mat_table_index[i];
  }
 
  /* mat_val */
  ref_mat->mat_val = (double *)HECMW_malloc(sizeof(double) * mat->n_mat_table);
  if (ref_mat->mat_val == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < mat->n_mat_table; i++) {
    ref_mat->mat_val[i] = mat->mat_val[i];
  }
 
  /* mat_temp */
  ref_mat->mat_temp = (double *)HECMW_malloc(sizeof(double) * mat->n_mat_table);
  if (ref_mat->mat_temp == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < mat->n_mat_table; i++) {
    ref_mat->mat_temp[i] = mat->mat_temp[i];
  }
 
  return HECMW_SUCCESS;
}
 
static int copy_mpc_info(const struct hecmwST_mpc *mpc,
                         struct hecmwST_mpc *ref_mpc) {
  int i;
 
  ref_mpc->n_mpc = mpc->n_mpc;
 
  if (mpc->n_mpc == 0) {
    ref_mpc->mpc_item  = NULL;
    ref_mpc->mpc_dof   = NULL;
    ref_mpc->mpc_val   = NULL;
    ref_mpc->mpc_const = NULL;
    ref_mpc->mpc_index = HECMW_malloc(sizeof(int));
    if (ref_mpc->mpc_index == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    ref_mpc->mpc_index[0] = 0;
    return HECMW_SUCCESS;
  }
 
#if 0
  HECMW_log(HECMW_LOG_ERROR, "refinement of MPC information is not supported\n");
  return HECMW_ERROR;
 
#else
  HECMW_log(HECMW_LOG_WARN, "MPC information is not refined\n");
 
  /* mpc_index */
  ref_mpc->mpc_index = (int *) HECMW_malloc( sizeof(int) * (mpc->n_mpc+1) );
  if( ref_mpc->mpc_index == NULL ) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for ( i = 0; i < mpc->n_mpc + 1; i++ ) {
    ref_mpc->mpc_index[i] = mpc->mpc_index[i];
  }
 
  /* mpc_item */
  ref_mpc->mpc_item = (int *) HECMW_malloc( sizeof(int) * mpc->mpc_index[mpc->n_mpc] );
  if( ref_mpc->mpc_item == NULL ) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for ( i = 0; i < mpc->mpc_index[mpc->n_mpc]; i++ ) {
    ref_mpc->mpc_item[i] = mpc->mpc_item[i];
  }
 
  /* mpc_dof */
  ref_mpc->mpc_dof = (int *) HECMW_malloc( sizeof(int) * mpc->mpc_index[mpc->n_mpc] );
  if( ref_mpc->mpc_dof == NULL ) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for ( i = 0; i < mpc->mpc_index[mpc->n_mpc]; i++ ) {
    ref_mpc->mpc_dof[i] = mpc->mpc_dof[i];
  }
 
  /* mpc_val */
  ref_mpc->mpc_val = (double *) HECMW_malloc( sizeof(double) * mpc->mpc_index[mpc->n_mpc] );
  if( ref_mpc->mpc_val == NULL ) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for ( i = 0; i < mpc->mpc_index[mpc->n_mpc]; i++ ) {
    ref_mpc->mpc_val[i] = mpc->mpc_val[i];
  }
 
  /* mpc_const */
  ref_mpc->mpc_const = (double *) HECMW_malloc( sizeof(double) * mpc->n_mpc );
  if( ref_mpc->mpc_const == NULL ) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for ( i = 0; i < mpc->n_mpc; i++ ) {
    ref_mpc->mpc_const[i] = mpc->mpc_const[i];
  }
 
  return HECMW_SUCCESS;
#endif
}
 
static int copy_amp_info(const struct hecmwST_amplitude *amp,
                         struct hecmwST_amplitude *ref_amp) {
  int i;
 
  ref_amp->n_amp = amp->n_amp;
 
  if (amp->n_amp <= 0) {
    ref_amp->amp_name            = NULL;
    ref_amp->amp_type_definition = NULL;
    ref_amp->amp_type_time       = NULL;
    ref_amp->amp_type_value      = NULL;
    ref_amp->amp_val             = NULL;
    ref_amp->amp_table           = NULL;
    ref_amp->amp_index           = (int *)HECMW_malloc(sizeof(int));
    if (ref_amp->amp_index == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    ref_amp->amp_index[0] = 0;
    return HECMW_SUCCESS;
  }
 
  /* amp_name */
  ref_amp->amp_name = (char **)HECMW_malloc(sizeof(char *) * amp->n_amp);
  if (ref_amp->amp_name == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < amp->n_amp; i++) {
    ref_amp->amp_name[i] = HECMW_strdup(amp->amp_name[i]);
    if (ref_amp->amp_name[i] == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
  }
 
  /* amp_type_definition */
  ref_amp->amp_type_definition = (int *)HECMW_malloc(sizeof(int) * amp->n_amp);
  if (ref_amp->amp_type_definition == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < amp->n_amp; i++) {
    ref_amp->amp_type_definition[i] = amp->amp_type_definition[i];
  }
 
  /* amp_type_time */
  ref_amp->amp_type_time = (int *)HECMW_malloc(sizeof(int) * amp->n_amp);
  if (ref_amp->amp_type_time == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < amp->n_amp; i++) {
    ref_amp->amp_type_time[i] = amp->amp_type_time[i];
  }
 
  /* amp_type_value */
  ref_amp->amp_type_value = (int *)HECMW_malloc(sizeof(int) * amp->n_amp);
  if (ref_amp->amp_type_value == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < amp->n_amp; i++) {
    ref_amp->amp_type_value[i] = amp->amp_type_value[i];
  }
 
  /* amp_index */
  ref_amp->amp_index = (int *)HECMW_malloc(sizeof(int) * (amp->n_amp + 1));
  if (ref_amp->amp_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < amp->n_amp + 1; i++) {
    ref_amp->amp_index[i] = amp->amp_index[i];
  }
 
  /* amp_val */
  ref_amp->amp_val =
      (double *)HECMW_malloc(sizeof(double) * amp->amp_index[amp->n_amp]);
  if (ref_amp->amp_val == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < amp->amp_index[amp->n_amp]; i++) {
    ref_amp->amp_val[i] = amp->amp_val[i];
  }
 
  /* amp_table */
  ref_amp->amp_table =
      (double *)HECMW_malloc(sizeof(double) * amp->amp_index[amp->n_amp]);
  if (ref_amp->amp_table == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < amp->amp_index[amp->n_amp]; i++) {
    ref_amp->amp_table[i] = amp->amp_table[i];
  }
 
  return HECMW_SUCCESS;
}
 
static int copy_contact_pair(const struct hecmwST_contact_pair *cp,
                             struct hecmwST_contact_pair *ref_cp) {
  int i;
 
  ref_cp->n_pair = cp->n_pair;
 
  /* name */
  ref_cp->name = (char **)HECMW_malloc(sizeof(char *) * cp->n_pair);
  if (ref_cp->name == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < cp->n_pair; i++) {
    ref_cp->name[i] = HECMW_strdup(cp->name[i]);
    if (ref_cp->name[i] == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
  }
 
  /* type */
  ref_cp->type = (int *)HECMW_malloc(sizeof(int) * cp->n_pair);
  if (ref_cp->type == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < cp->n_pair; i++) {
    ref_cp->type[i] = cp->type[i];
  }
 
  /* slave_grp_id */
  ref_cp->slave_grp_id = (int *)HECMW_malloc(sizeof(int) * cp->n_pair);
  if (ref_cp->slave_grp_id == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < cp->n_pair; i++) {
    ref_cp->slave_grp_id[i] = cp->slave_grp_id[i];
  }
 
  /* slave_orisgrp_id */
  ref_cp->slave_orisgrp_id = (int *)HECMW_malloc(sizeof(int) * cp->n_pair);
  if (ref_cp->slave_orisgrp_id == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < cp->n_pair; i++) {
    ref_cp->slave_orisgrp_id[i] = cp->slave_orisgrp_id[i];
  }
 
  /* master_grp_id */
  ref_cp->master_grp_id = (int *)HECMW_malloc(sizeof(int) * cp->n_pair);
  if (ref_cp->master_grp_id == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i < cp->n_pair; i++) {
    ref_cp->master_grp_id[i] = cp->master_grp_id[i];
  }
 
  return HECMW_SUCCESS;
}
 
static int copy_unchanging_info(const struct hecmwST_local_mesh *mesh,
                                struct hecmwST_local_mesh *ref_mesh) {
  HECMW_log(HECMW_LOG_DEBUG, "Started copying unchanging info...\n");
 
  if (copy_global_info(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (copy_node_ndof(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (copy_elem_etype(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (copy_comm_info(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (copy_adapt_info(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (copy_section_info(mesh->section, ref_mesh->section) != HECMW_SUCCESS)
    return HECMW_ERROR;
  if (copy_material_info(mesh->material, ref_mesh->material) != HECMW_SUCCESS)
    return HECMW_ERROR;
  if (copy_mpc_info(mesh->mpc, ref_mesh->mpc) != HECMW_SUCCESS)
    return HECMW_ERROR;
  if (copy_amp_info(mesh->amp, ref_mesh->amp) != HECMW_SUCCESS)
    return HECMW_ERROR;
  if (copy_contact_pair(mesh->contact_pair, ref_mesh->contact_pair) !=
      HECMW_SUCCESS)
    return HECMW_ERROR;
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished copying unchanging info.\n");
  return HECMW_SUCCESS;
}
 
/*============================================================================*/
/*                                                                            */
/*  renumber nodes                                                            */
/*                                                                            */
/*============================================================================*/
static int renumber_nodes_generate_tables(struct hecmwST_local_mesh *mesh) {
  int i;
  int count_in, count_ex, count_pex;
  int *old2new, *new2old;
 
  if (mesh->n_node_gross == mesh->nn_internal) return HECMW_SUCCESS;
 
  HECMW_log(HECMW_LOG_DEBUG, "Started generating renumbering tables...\n");
 
  old2new = (int *)HECMW_malloc(sizeof(int) * mesh->n_node_gross);
  if (old2new == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  new2old = (int *)HECMW_malloc(sizeof(int) * mesh->n_node_gross);
  if (new2old == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  count_in  = 0;
  count_ex  = mesh->nn_internal;
  count_pex = mesh->n_node;
  for (i = 0; i < mesh->n_node_gross; i++) {
    int rank = mesh->node_ID[i * 2 + 1];
    if (rank == mesh->my_rank) {
      old2new[i]        = count_in;
      new2old[count_in] = i;
      count_in++;
    } else if (rank >= 0) {
      old2new[i]        = count_ex;
      new2old[count_ex] = i;
      count_ex++;
    } else {
      old2new[i]         = count_pex;
      new2old[count_pex] = i;
      count_pex++;
    }
  }
  mesh->node_old2new = old2new;
  mesh->node_new2old = new2old;
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished generating renumbering tables.\n");
  return HECMW_SUCCESS;
}
 
static int reorder_nodes(struct hecmwST_local_mesh *mesh, int *old2new,
                         int *new2old) {
  int i;
 
#ifndef NDEBUG
  for (i = 0; i < mesh->n_node_gross; i++) {
    HECMW_assert(new2old[old2new[i]] == i);
  }
#endif
 
  /*
   * Reorder using new2old
   */
  /* node_ID */
  {
    int *new_node_ID =
        (int *)HECMW_malloc(sizeof(int) * mesh->n_node_gross * 2);
    if (new_node_ID == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < mesh->n_node_gross; i++) {
      int old                = new2old[i];
      new_node_ID[2 * i]     = mesh->node_ID[2 * old];
      new_node_ID[2 * i + 1] = mesh->node_ID[2 * old + 1];
    }
    HECMW_free(mesh->node_ID);
    mesh->node_ID = new_node_ID;
  }
  /* global_node_ID */
  {
    int *new_global_node_ID =
        (int *)HECMW_malloc(sizeof(int) * mesh->n_node_gross);
    if (new_global_node_ID == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < mesh->n_node_gross; i++) {
      int old               = new2old[i];
      new_global_node_ID[i] = mesh->global_node_ID[old];
    }
    HECMW_free(mesh->global_node_ID);
    mesh->global_node_ID = new_global_node_ID;
  }
  /* node */
  {
    double *new_node =
        (double *)HECMW_malloc(sizeof(double) * mesh->n_node_gross * 3);
    if (new_node == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < mesh->n_node_gross; i++) {
      int old = new2old[i];
      int j;
      for (j = 0; j < 3; j++) new_node[3 * i + j] = mesh->node[3 * old + j];
    }
    HECMW_free(mesh->node);
    mesh->node = new_node;
  }
  /* node_init_val_index, node_init_val_item */
  if (mesh->hecmw_flag_initcon && mesh->n_node_gross > 0) {
    /* node_init_val_index */
    int *new_node_init_val_index =
        (int *)HECMW_malloc(sizeof(int) * (mesh->n_node_gross + 1));
    if (new_node_init_val_index == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    new_node_init_val_index[0] = 0;
    for (i = 0; i < mesh->n_node_gross; i++) {
      int old = new2old[i];
      int ninit =
          mesh->node_init_val_index[old + 1] - mesh->node_init_val_index[old];
      new_node_init_val_index[i + 1] = new_node_init_val_index[i] + ninit;
    }
 
    /* node_init_val_item */
    if (mesh->node_init_val_index[mesh->n_node_gross] != 0) {
      double *new_node_init_val_item = (double *)HECMW_malloc(
          sizeof(double) * mesh->node_init_val_index[mesh->n_node_gross]);
      if (new_node_init_val_item == NULL) {
        HECMW_set_error(errno, "");
        return HECMW_ERROR;
      }
      for (i = 0; i < mesh->n_node_gross; i++) {
        int old     = new2old[i];
        int idx_new = new_node_init_val_index[i];
        int idx_old = mesh->node_init_val_index[old];
        int ninit   = mesh->node_init_val_index[old + 1] - idx_old;
        int j;
        for (j = 0; j < ninit; j++) {
          new_node_init_val_item[idx_new + j] =
              mesh->node_init_val_item[idx_old + j];
        }
      }
      HECMW_free(mesh->node_init_val_item);
      mesh->node_init_val_item = new_node_init_val_item;
    }
    HECMW_free(mesh->node_init_val_index);
    mesh->node_init_val_index = new_node_init_val_index;
  }
 
  /*
   * Update using old2new
   */
  /* elem_node_item */
  {
    long long ii;
    for (ii = 0; ii < mesh->elem_node_index[mesh->n_elem_gross]; ii++) {
      int old                 = mesh->elem_node_item[ii];
      int new                 = old2new[old - 1] + 1;
      mesh->elem_node_item[ii] = new;
    }
  }
  /* import_item */
  for (i = 0; i < mesh->import_index[mesh->n_neighbor_pe]; i++) {
    int old              = mesh->import_item[i];
    int new              = old2new[old - 1] + 1;
    mesh->import_item[i] = new;
  }
  /* export_item */
  for (i = 0; i < mesh->export_index[mesh->n_neighbor_pe]; i++) {
    int old              = mesh->export_item[i];
    int new              = old2new[old - 1] + 1;
    mesh->export_item[i] = new;
  }
  /* mpc->mpc_item */
  for (i = 0; i < mesh->mpc->mpc_index[mesh->mpc->n_mpc]; i++) {
    int old                = mesh->mpc->mpc_item[i];
    int new                = old2new[old - 1] + 1;
    mesh->mpc->mpc_item[i] = new;
  }
  /* node_group->grp_item */
  for (i = 0; i < mesh->node_group->grp_index[mesh->node_group->n_grp]; i++) {
    int old                       = mesh->node_group->grp_item[i];
    int new                       = old2new[old - 1] + 1;
    mesh->node_group->grp_item[i] = new;
  }
  return HECMW_SUCCESS;
}
 
static int delete_external_items(int n_elem, int n_grp, int *index, int *item,
                                 int n_memb) {
  int i, j, k;
  int n_del = 0;
  int start = 0;
  for (i = 0; i < n_grp; i++) {
    int end = index[i + 1];
    for (j = start; j < end; j++) {
      if (item[j * n_memb] > n_elem) {
        n_del++;
      } else {
        for (k                           = 0; k < n_memb; k++)
          item[(j - n_del) * n_memb + k] = item[j * n_memb + k];
      }
    }
    /* start for next group */
    start        = end;
    index[i + 1] = end - n_del;
  }
  return HECMW_SUCCESS;
}
 
static int delete_external_equations(int n_node, struct hecmwST_mpc *mpc) {
  int iold, inew, jold, jnew;
  int *active_eqn;
  int n_mpc_new;
  int *mpc_index_new, *mpc_item_new, *mpc_dof_new;
  double *mpc_val_new, *mpc_const_new;
  active_eqn = (int *) HECMW_malloc(sizeof(int) * mpc->n_mpc);
  if (active_eqn == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  n_mpc_new = 0;
  for (iold = 0; iold < mpc->n_mpc; iold++) {
    int is_ext = 0;
    for (jold = mpc->mpc_index[iold]; jold < mpc->mpc_index[iold+1]; jold++) {
      if (mpc->mpc_item[jold] > n_node) {
        is_ext = 1;
        break;
      }
    }
    if (is_ext == 0) {
      active_eqn[n_mpc_new] = iold;
      n_mpc_new++;
    }
  }
  mpc_index_new = (int *) HECMW_malloc(sizeof(int) * (n_mpc_new + 1));
  if (mpc_index_new == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  mpc_index_new[0] = 0;
  if (n_mpc_new == 0) {
    mpc_item_new = NULL;
    mpc_dof_new = NULL;
    mpc_val_new = NULL;
    mpc_const_new = NULL;
  } else {
    for (inew = 0; inew < n_mpc_new; inew++) {
      iold = active_eqn[inew];
      mpc_index_new[inew+1] = mpc_index_new[inew] + (mpc->mpc_index[iold+1] - mpc->mpc_index[iold]);
    }
    mpc_item_new = (int *) HECMW_malloc(sizeof(int) * mpc_index_new[n_mpc_new]);
    if (mpc_item_new == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    mpc_dof_new = (int *) HECMW_malloc(sizeof(int) * mpc_index_new[n_mpc_new]);
    if (mpc_dof_new == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    mpc_val_new = (double *) HECMW_malloc(sizeof(double) * mpc_index_new[n_mpc_new]);
    if (mpc_val_new == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    mpc_const_new = (double *) HECMW_malloc(sizeof(double) * n_mpc_new);
    if (mpc_const_new == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (inew = 0; inew < n_mpc_new; inew++) {
      iold = active_eqn[inew];
      for (jold = mpc->mpc_index[iold], jnew = mpc_index_new[inew]; jold < mpc->mpc_index[iold+1]; jold++, jnew++) {
        mpc_item_new[jnew] = mpc->mpc_item[jold];
        mpc_dof_new[jnew]  = mpc->mpc_dof[jold];
        mpc_val_new[jnew]  = mpc->mpc_val[jold];
      }
      mpc_const_new[inew] = mpc->mpc_const[iold];
    }
  }
  HECMW_free(active_eqn);
  mpc->n_mpc = n_mpc_new;
  HECMW_free(mpc->mpc_index);
  HECMW_free(mpc->mpc_item);
  HECMW_free(mpc->mpc_dof);
  HECMW_free(mpc->mpc_val);
  HECMW_free(mpc->mpc_const);
  mpc->mpc_index = mpc_index_new;
  mpc->mpc_item  = mpc_item_new;
  mpc->mpc_dof   = mpc_dof_new;
  mpc->mpc_val   = mpc_val_new;
  mpc->mpc_const = mpc_const_new;
  return HECMW_SUCCESS;
}
 
static int renumber_nodes(struct hecmwST_local_mesh *mesh) {
  if (mesh->n_node_gross == mesh->nn_internal) return HECMW_SUCCESS;
 
  HECMW_log(HECMW_LOG_DEBUG, "Started renumbering nodes...\n");
 
  if (reorder_nodes(mesh, mesh->node_old2new, mesh->node_new2old) !=
      HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
 
  if (delete_external_items(mesh->n_node, mesh->node_group->n_grp,
                            mesh->node_group->grp_index, mesh->node_group->grp_item,
                            1) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
 
  if (delete_external_equations(mesh->n_node, mesh->mpc) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished renumbering nodes.\n");
  return HECMW_SUCCESS;
}
 
static int renumber_back_nodes(struct hecmwST_local_mesh *mesh) {
  if (mesh->n_node_gross == mesh->nn_internal) return HECMW_SUCCESS;
 
  HECMW_log(HECMW_LOG_DEBUG, "Started renumbering back nodes...\n");
 
  if (reorder_nodes(mesh, mesh->node_new2old, mesh->node_old2new) !=
      HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished renumbering back nodes.\n");
  return HECMW_SUCCESS;
}
 
/*============================================================================*/
/*                                                                            */
/*  renumber elements                                                         */
/*                                                                            */
/*============================================================================*/
static int renumber_elements_generate_tables(struct hecmwST_local_mesh *mesh) {
  int i;
  int count_rel, count_pex;
  int *old2new, *new2old;
 
  if (mesh->n_node_gross == mesh->nn_internal) return HECMW_SUCCESS;
 
  HECMW_log(HECMW_LOG_DEBUG, "Started generating renumbering tables...\n");
 
  old2new = (int *)HECMW_malloc(sizeof(int) * mesh->n_elem_gross);
  if (old2new == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  new2old = (int *)HECMW_malloc(sizeof(int) * mesh->n_elem_gross);
  if (new2old == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  count_rel = 0;
  count_pex = mesh->n_elem;
  for (i = 0; i < mesh->n_elem_gross; i++) {
    int rank = mesh->elem_ID[i * 2 + 1];
    if (rank >= 0) {
      old2new[i]         = count_rel;
      new2old[count_rel] = i;
      count_rel++;
    } else {
      old2new[i]         = count_pex;
      new2old[count_pex] = i;
      count_pex++;
    }
  }
  mesh->elem_old2new = old2new;
  mesh->elem_new2old = new2old;
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished generating renumbering tables.\n");
  return HECMW_SUCCESS;
}
 
static int reorder_elems(struct hecmwST_local_mesh *mesh, int *old2new,
                         int *new2old) {
  int i;
 
#ifndef NDEBUG
  for (i = 0; i < mesh->n_elem_gross; i++) {
    HECMW_assert(new2old[old2new[i]] == i);
  }
#endif
 
  /*
   * Reorder using new2old
   */
  /* elem_ID */
  {
    int *new_elem_ID =
        (int *)HECMW_malloc(sizeof(int) * mesh->n_elem_gross * 2);
    if (new_elem_ID == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < mesh->n_elem_gross; i++) {
      int old                = new2old[i];
      new_elem_ID[2 * i]     = mesh->elem_ID[2 * old];
      new_elem_ID[2 * i + 1] = mesh->elem_ID[2 * old + 1];
    }
    HECMW_free(mesh->elem_ID);
    mesh->elem_ID = new_elem_ID;
  }
  /* global_elem_ID */
  {
    int *new_global_elem_ID =
        (int *)HECMW_malloc(sizeof(int) * mesh->n_elem_gross);
    if (new_global_elem_ID == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < mesh->n_elem_gross; i++) {
      int old               = new2old[i];
      new_global_elem_ID[i] = mesh->global_elem_ID[old];
    }
    HECMW_free(mesh->global_elem_ID);
    mesh->global_elem_ID = new_global_elem_ID;
  }
  /* elem_type */
  {
    int *new_elem_type = (int *)HECMW_malloc(sizeof(int) * mesh->n_elem_gross);
    if (new_elem_type == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < mesh->n_elem_gross; i++) {
      int old          = new2old[i];
      new_elem_type[i] = mesh->elem_type[old];
    }
    HECMW_free(mesh->elem_type);
    mesh->elem_type = new_elem_type;
  }
  /* elem_type_index, elem_type_item */
  {
    for (i = 0; i < mesh->n_elem_type; i++) {
      int j;
      for (j = mesh->elem_type_index[i]; j < mesh->n_elem; j++) {
        if (mesh->elem_type[j] != mesh->elem_type_item[i]) {
          break;
        }
      }
      mesh->elem_type_index[i + 1] = j;
    }
  }
  /* elem_node_index, elem_node_item */
  {
    long long *new_elem_node_index =
        (long long *)HECMW_malloc(sizeof(long long) * (mesh->n_elem_gross + 1));
    int *new_elem_node_item = (int *)HECMW_malloc(
        sizeof(int) * mesh->elem_node_index[mesh->n_elem_gross]);
    if (new_elem_node_index == NULL || new_elem_node_item == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    new_elem_node_index[0] = 0;
    for (i = 0; i < mesh->n_elem_gross; i++) {
      int old = new2old[i];
      long long nn  = mesh->elem_node_index[old + 1] - mesh->elem_node_index[old];
      new_elem_node_index[i + 1] = new_elem_node_index[i] + nn;
    }
    for (i = 0; i < mesh->n_elem_gross; i++) {
      int old       = new2old[i];
      long long start_old = mesh->elem_node_index[old];
      long long start_new = new_elem_node_index[i];
      long long nn        = new_elem_node_index[i + 1] - start_new;
      long long j;
      for (j                              = 0; j < nn; j++)
        new_elem_node_item[start_new + j] = mesh->elem_node_item[start_old + j];
    }
    HECMW_free(mesh->elem_node_index);
    HECMW_free(mesh->elem_node_item);
    mesh->elem_node_index = new_elem_node_index;
    mesh->elem_node_item  = new_elem_node_item;
  }
  /* section_ID */
  {
    int *new_section_ID = (int *)HECMW_malloc(sizeof(int) * mesh->n_elem_gross);
    if (new_section_ID == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < mesh->n_elem_gross; i++) {
      int old           = new2old[i];
      new_section_ID[i] = mesh->section_ID[old];
    }
    HECMW_free(mesh->section_ID);
    mesh->section_ID = new_section_ID;
  }
  /* elem_mat_ID_index, elem_mat_ID_item */
  {
    int *new_emID_index, *new_emID_item;
    new_emID_index =
        (int *)HECMW_malloc(sizeof(int) * (mesh->n_elem_gross + 1));
    if (new_emID_index == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    new_emID_item = (int *)HECMW_malloc(
        sizeof(int) * mesh->elem_mat_ID_index[mesh->n_elem_gross]);
    if (new_emID_item == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    new_emID_index[0] = 0;
    for (i = 0; i < mesh->n_elem_gross; i++) {
      int old, js_old, je_old, len, js_new, j;
      old                   = new2old[i];
      js_old                = mesh->elem_mat_ID_index[old];
      je_old                = mesh->elem_mat_ID_index[old + 1];
      len                   = je_old - js_old;
      js_new                = new_emID_index[i];
      new_emID_index[i + 1] = js_new + len;
      for (j = 0; j < len; j++) {
        new_emID_item[js_new + j] = mesh->elem_mat_ID_item[js_old + j];
      }
    }
    HECMW_free(mesh->elem_mat_ID_index);
    HECMW_free(mesh->elem_mat_ID_item);
    mesh->elem_mat_ID_index = new_emID_index;
    mesh->elem_mat_ID_item  = new_emID_item;
  }
 
  /*
   * Update using old2new
   */
  /* elem_internal_list */
  for (i = 0; i < mesh->ne_internal; i++) {
    int old                     = mesh->elem_internal_list[i];
    int new                     = old2new[old - 1] + 1;
    mesh->elem_internal_list[i] = new;
  }
  /* shared_item */
  for (i = 0; i < mesh->shared_index[mesh->n_neighbor_pe]; i++) {
    int old              = mesh->shared_item[i];
    int new              = old2new[old - 1] + 1;
    mesh->shared_item[i] = new;
  }
  /* elem_groups->grp_item */
  for (i = 0; i < mesh->elem_group->grp_index[mesh->elem_group->n_grp]; i++) {
    int old                       = mesh->elem_group->grp_item[i];
    int new                       = old2new[old - 1] + 1;
    mesh->elem_group->grp_item[i] = new;
  }
  /* surf_groups->grp_item */
  for (i = 0; i < mesh->surf_group->grp_index[mesh->surf_group->n_grp]; i++) {
    int old                           = mesh->surf_group->grp_item[2 * i];
    int new                           = old2new[old - 1] + 1;
    mesh->surf_group->grp_item[2 * i] = new;
  }
  return HECMW_SUCCESS;
}
 
static int renumber_elements(struct hecmwST_local_mesh *mesh) {
  if (mesh->n_elem_gross == mesh->n_elem) return HECMW_SUCCESS;
 
  HECMW_log(HECMW_LOG_DEBUG, "Started renumbering elements...\n");
 
  if (reorder_elems(mesh, mesh->elem_old2new, mesh->elem_new2old) !=
      HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
 
  /* if (delete_external_items( mesh->n_elem, mesh->n_neighbor_pe, */
  /*                            mesh->shared_index, mesh->shared_item, 1 ) != */
  /*     HECMW_SUCCESS) { */
  /*   return HECMW_ERROR; */
  /* } */
  if (delete_external_items(mesh->n_elem, mesh->elem_group->n_grp,
                            mesh->elem_group->grp_index, mesh->elem_group->grp_item,
                            1) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
  if (delete_external_items(mesh->n_elem, mesh->surf_group->n_grp,
                            mesh->surf_group->grp_index, mesh->surf_group->grp_item,
                            2) != HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished renumbering elements.\n");
  return HECMW_SUCCESS;
}
 
static int renumber_back_elements(struct hecmwST_local_mesh *mesh) {
  if (mesh->n_elem_gross == mesh->n_elem) return HECMW_SUCCESS;
 
  HECMW_log(HECMW_LOG_DEBUG, "Started renumbering back elements...\n");
 
  if (reorder_elems(mesh, mesh->elem_new2old, mesh->elem_old2new) !=
      HECMW_SUCCESS) {
    return HECMW_ERROR;
  }
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished renumbering back elements.\n");
  return HECMW_SUCCESS;
}
 
/*============================================================================*/
/*                                                                            */
/*  replicate slave nodes of contact pairs                                    */
/*                                                                            */
/*============================================================================*/
/*
 * The partitioner gives each domain owning a master element of a contact pair
 * all the slave nodes as external nodes, which are not shared by any element.
 * The communication tables rebuilt from the refined shared elements drop them,
 * and the slave nodes generated by refinement exist only in the domains having
 * the slave elements.  This gives them again to the domains owning master
 * elements, as the partitioner does for the master-owner scheme.  SURF_SURF
 * pairs are not handled here; their slave elements are replicated instead.
 */
static int is_node_slave_pair(const struct hecmwST_contact_pair *cp, int i) {
  return cp->type[i] == HECMW_CONTACT_TYPE_NODE_SURF ||
         cp->type[i] == HECMW_CONTACT_TYPE_NODE_ELEM;
}
 
static int has_internal_master(const struct hecmwST_local_mesh *mesh,
                               int ipair) {
  const struct hecmwST_contact_pair *cp = mesh->contact_pair;
  int gid = cp->master_grp_id[ipair];
  int j, elem;
 
  if (cp->type[ipair] == HECMW_CONTACT_TYPE_NODE_ELEM) {
    const struct hecmwST_elem_grp *grp = mesh->elem_group;
    for (j = grp->grp_index[gid - 1]; j < grp->grp_index[gid]; j++) {
      elem = grp->grp_item[j];
      if (mesh->elem_ID[2 * (elem - 1) + 1] == mesh->my_rank) return 1;
    }
  } else {
    const struct hecmwST_surf_grp *grp = mesh->surf_group;
    for (j = grp->grp_index[gid - 1]; j < grp->grp_index[gid]; j++) {
      elem = grp->grp_item[2 * j];
      if (mesh->elem_ID[2 * (elem - 1) + 1] == mesh->my_rank) return 1;
    }
  }
  return 0;
}
 
/* group membership of internal nodes, except for group ALL */
static int get_internal_node_membership(const struct hecmwST_local_mesh *mesh,
                                        int **memb_index, int **memb_item) {
  const struct hecmwST_node_grp *grp = mesh->node_group;
  int nn_int = mesh->nn_internal;
  int i, j, node;
  int *cnt;
 
  *memb_index = (int *)HECMW_calloc(nn_int + 1, sizeof(int));
  if (*memb_index == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 1; i < grp->n_grp; i++) {
    for (j = grp->grp_index[i]; j < grp->grp_index[i + 1]; j++) {
      node = grp->grp_item[j];
      if (node <= nn_int) (*memb_index)[node]++;
    }
  }
  for (i = 0; i < nn_int; i++) (*memb_index)[i + 1] += (*memb_index)[i];
 
  *memb_item = (int *)HECMW_malloc(sizeof(int) * ((*memb_index)[nn_int] + 1));
  cnt        = (int *)HECMW_calloc(nn_int, sizeof(int));
  if (*memb_item == NULL || cnt == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 1; i < grp->n_grp; i++) {
    for (j = grp->grp_index[i]; j < grp->grp_index[i + 1]; j++) {
      node = grp->grp_item[j];
      if (node > nn_int) continue;
      (*memb_item)[(*memb_index)[node - 1] + cnt[node - 1]] = i;
      cnt[node - 1]++;
    }
  }
  HECMW_free(cnt);
  return HECMW_SUCCESS;
}
 
/*
 * pack the internal slave nodes of the pairs whose master the neighbor owns:
 *   ibuf: n, then for each node: local ID, global ID, number of groups, groups,
 *         number of initial values
 *   dbuf: for each node: coordinates, initial values
 */
static int pack_contact_slave_nodes(const struct hecmwST_local_mesh *mesh,
                                    const int *nb_flag, const int *memb_index,
                                    const int *memb_item,
                                    struct hecmw_varray_int *ibuf,
                                    double **dbuf, int *dlen) {
  const struct hecmwST_contact_pair *cp = mesh->contact_pair;
  const struct hecmwST_node_grp *grp    = mesh->node_group;
  struct hecmw_bit_array mark;
  struct hecmw_varray_int nodes;
  int nn_int = mesh->nn_internal;
  int i, j, k, node, n, ninit, d;
 
  if (HECMW_bit_array_init(&mark, nn_int) != HECMW_SUCCESS) return HECMW_ERROR;
  for (i = 0; i < cp->n_pair; i++) {
    int gid;
    if (!nb_flag[i] || !is_node_slave_pair(cp, i)) continue;
    gid = cp->slave_grp_id[i];
    for (j = grp->grp_index[gid - 1]; j < grp->grp_index[gid]; j++) {
      node = grp->grp_item[j];
      if (node <= nn_int) HECMW_bit_array_set(&mark, node - 1);
    }
  }
  HECMW_varray_int_init(&nodes);
  for (i = 0; i < nn_int; i++) {
    if (HECMW_bit_array_get(&mark, i)) HECMW_varray_int_append(&nodes, i + 1);
  }
  HECMW_bit_array_finalize(&mark);
 
  n     = HECMW_varray_int_nval(&nodes);
  *dlen = 0;
  HECMW_varray_int_append(ibuf, n);
  for (k = 0; k < n; k++) {
    node  = HECMW_varray_int_get(&nodes, k);
    ninit = 0;
    if (mesh->hecmw_flag_initcon) {
      ninit = mesh->node_init_val_index[node] -
              mesh->node_init_val_index[node - 1];
    }
    HECMW_varray_int_append(ibuf, node);
    HECMW_varray_int_append(ibuf, mesh->global_node_ID[node - 1]);
    HECMW_varray_int_append(ibuf, memb_index[node] - memb_index[node - 1]);
    for (j = memb_index[node - 1]; j < memb_index[node]; j++) {
      HECMW_varray_int_append(ibuf, memb_item[j]);
    }
    HECMW_varray_int_append(ibuf, ninit);
    *dlen += 3 + ninit;
  }
 
  *dbuf = (double *)HECMW_malloc(sizeof(double) * (*dlen + 1));
  if (*dbuf == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  d = 0;
  for (k = 0; k < n; k++) {
    node = HECMW_varray_int_get(&nodes, k);
    for (j = 0; j < 3; j++) (*dbuf)[d++] = mesh->node[3 * (node - 1) + j];
    if (mesh->hecmw_flag_initcon) {
      for (j = mesh->node_init_val_index[node - 1];
           j < mesh->node_init_val_index[node]; j++) {
        (*dbuf)[d++] = mesh->node_init_val_item[j];
      }
    }
  }
  HECMW_varray_int_finalize(&nodes);
  return HECMW_SUCCESS;
}
 
/* nodes to be added as external nodes */
struct added_nodes {
  struct hecmw_varray_int node_ID;  /* local ID in the owner, owner rank */
  struct hecmw_varray_int global_ID;
  struct hecmw_varray_int memb_index;
  struct hecmw_varray_int memb_item;
  struct hecmw_varray_int init_index;
  int n_coord, n_init;
  double *coord;
  double *init;
};
 
/*
 * select the received nodes not present yet, and append them to added; the
 * local IDs of the selected nodes in the owner are appended to lid_added
 */
static int unpack_contact_slave_nodes(const struct hecmwST_local_mesh *mesh,
                                      int rank, const int *ibuf,
                                      const double *dbuf, int dlen,
                                      struct added_nodes *added,
                                      struct hecmw_varray_int *lid_added) {
  int n = ibuf[0];
  int i, j, k, ip, dp, lid, nmemb, ninit, max_lid;
  int *lid2node;
 
  added->coord = (double *)HECMW_realloc(
      added->coord, sizeof(double) * (added->n_coord + 3 * n + 1));
  added->init = (double *)HECMW_realloc(
      added->init, sizeof(double) * (added->n_init + dlen - 3 * n + 1));
  if (added->coord == NULL || added->init == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
 
  /* external nodes from the owner, looked up by the local ID in the owner */
  max_lid = 0;
  for (ip = 1, k = 0; k < n; k++) {
    lid = ibuf[ip];
    if (lid > max_lid) max_lid = lid;
    ip += 3 + ibuf[ip + 2] + 1;
  }
  lid2node = (int *)HECMW_malloc(sizeof(int) * (max_lid + 1));
  if (lid2node == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  for (i = 0; i <= max_lid; i++) lid2node[i] = -1;
  for (i = mesh->nn_internal; i < mesh->n_node; i++) {
    if (mesh->node_ID[2 * i + 1] != rank) continue;
    lid = mesh->node_ID[2 * i];
    if (lid <= max_lid) lid2node[lid] = i;
  }
 
  for (ip = 1, dp = 0, k = 0; k < n; k++) {
    lid   = ibuf[ip];
    nmemb = ibuf[ip + 2];
    ninit = ibuf[ip + 3 + nmemb];
    if (lid2node[lid] < 0) {
      HECMW_varray_int_append(lid_added, lid);
      HECMW_varray_int_append(&added->node_ID, lid);
      HECMW_varray_int_append(&added->node_ID, rank);
      HECMW_varray_int_append(&added->global_ID, ibuf[ip + 1]);
      for (j = 0; j < nmemb; j++) {
        HECMW_varray_int_append(&added->memb_item, ibuf[ip + 3 + j]);
      }
      HECMW_varray_int_append(&added->memb_index,
                              HECMW_varray_int_nval(&added->memb_item));
      for (j = 0; j < 3; j++) added->coord[added->n_coord++] = dbuf[dp + j];
      for (j = 0; j < ninit; j++) {
        added->init[added->n_init++] = dbuf[dp + 3 + j];
      }
      HECMW_varray_int_append(&added->init_index, added->n_init);
    }
    ip += 3 + nmemb + 1;
    dp += 3 + ninit;
  }
  HECMW_free(lid2node);
  return HECMW_SUCCESS;
}
 
/* insert the added nodes after the external nodes */
static int insert_external_nodes(struct hecmwST_local_mesh *mesh,
                                 const struct added_nodes *added,
                                 struct hecmw_varray_int *imp_added,
                                 struct hecmw_varray_int *exp_added) {
  int n_add  = HECMW_varray_int_nval(&added->global_ID);
  int nn_old = mesh->n_node;
  int ng_old = mesh->n_node_gross;
  int ng_new = ng_old + n_add;
  int n_exp_add = 0;
  int i, j, k;
 
  for (i = 0; i < mesh->n_neighbor_pe; i++) {
    n_exp_add += HECMW_varray_int_nval(exp_added + i);
  }
  if (n_add == 0 && n_exp_add == 0) return HECMW_SUCCESS;
 
#define SHIFTED(node) ((node) > nn_old ? (node) + n_add : (node))
  /* node, node_ID, global_node_ID */
  {
    double *new_node    = (double *)HECMW_malloc(sizeof(double) * ng_new * 3);
    int *new_node_ID    = (int *)HECMW_malloc(sizeof(int) * ng_new * 2);
    int *new_global_ID  = (int *)HECMW_malloc(sizeof(int) * ng_new);
    if (new_node == NULL || new_node_ID == NULL || new_global_ID == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < ng_old; i++) {
      int inew = SHIFTED(i + 1) - 1;
      for (j = 0; j < 3; j++) new_node[3 * inew + j] = mesh->node[3 * i + j];
      new_node_ID[2 * inew]     = mesh->node_ID[2 * i];
      new_node_ID[2 * inew + 1] = mesh->node_ID[2 * i + 1];
      new_global_ID[inew]       = mesh->global_node_ID[i];
    }
    for (k = 0; k < n_add; k++) {
      int inew = nn_old + k;
      for (j = 0; j < 3; j++) new_node[3 * inew + j] = added->coord[3 * k + j];
      new_node_ID[2 * inew] = HECMW_varray_int_get(&added->node_ID, 2 * k);
      new_node_ID[2 * inew + 1] =
          HECMW_varray_int_get(&added->node_ID, 2 * k + 1);
      new_global_ID[inew]       = HECMW_varray_int_get(&added->global_ID, k);
    }
    HECMW_free(mesh->node);
    HECMW_free(mesh->node_ID);
    HECMW_free(mesh->global_node_ID);
    mesh->node           = new_node;
    mesh->node_ID        = new_node_ID;
    mesh->global_node_ID = new_global_ID;
  }
  /* node_init_val_index, node_init_val_item */
  if (mesh->hecmw_flag_initcon && ng_old > 0) {
    int n_val       = mesh->node_init_val_index[ng_old] + added->n_init;
    int *new_index  = (int *)HECMW_malloc(sizeof(int) * (ng_new + 1));
    double *new_val = (double *)HECMW_malloc(sizeof(double) * (n_val + 1));
    int pos         = 0;
    if (new_index == NULL || new_val == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    new_index[0] = 0;
    for (i = 0; i < ng_new; i++) {
      if (i < nn_old || i >= nn_old + n_add) {
        int iold = i < nn_old ? i : i - n_add;
        for (j = mesh->node_init_val_index[iold];
             j < mesh->node_init_val_index[iold + 1]; j++) {
          new_val[pos++] = mesh->node_init_val_item[j];
        }
      } else {
        k = i - nn_old;
        for (j = (k == 0 ? 0 : HECMW_varray_int_get(&added->init_index, k - 1));
             j < HECMW_varray_int_get(&added->init_index, k); j++) {
          new_val[pos++] = added->init[j];
        }
      }
      new_index[i + 1] = pos;
    }
    HECMW_free(mesh->node_init_val_index);
    HECMW_free(mesh->node_init_val_item);
    mesh->node_init_val_index = new_index;
    mesh->node_init_val_item  = new_val;
  }
  /* elem_node_item */
  {
    long long ii;
    for (ii = 0; ii < mesh->elem_node_index[mesh->n_elem_gross]; ii++) {
      mesh->elem_node_item[ii] = SHIFTED(mesh->elem_node_item[ii]);
    }
  }
  /* mpc->mpc_item */
  for (i = 0; i < mesh->mpc->mpc_index[mesh->mpc->n_mpc]; i++) {
    mesh->mpc->mpc_item[i] = SHIFTED(mesh->mpc->mpc_item[i]);
  }
  /* node_group: ALL takes every added node, the others as in the owner */
  {
    struct hecmwST_node_grp *grp = mesh->node_group;
    int n_item = grp->grp_index[grp->n_grp] + n_add +
                 (int)HECMW_varray_int_nval(&added->memb_item);
    int *new_index = (int *)HECMW_malloc(sizeof(int) * (grp->n_grp + 1));
    int *new_item  = (int *)HECMW_malloc(sizeof(int) * (n_item + 1));
    int pos        = 0;
    if (new_index == NULL || new_item == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    new_index[0] = 0;
    for (i = 0; i < grp->n_grp; i++) {
      for (j = grp->grp_index[i]; j < grp->grp_index[i + 1]; j++) {
        new_item[pos++] = SHIFTED(grp->grp_item[j]);
      }
      for (k = 0; k < n_add; k++) {
        int js = k == 0 ? 0 : HECMW_varray_int_get(&added->memb_index, k - 1);
        int je = HECMW_varray_int_get(&added->memb_index, k);
        int in_grp = (i == 0);
        for (j = js; j < je && !in_grp; j++) {
          if (HECMW_varray_int_get(&added->memb_item, j) == i) in_grp = 1;
        }
        if (in_grp) new_item[pos++] = nn_old + k + 1;
      }
      new_index[i + 1] = pos;
    }
    HECMW_free(grp->grp_index);
    HECMW_free(grp->grp_item);
    grp->grp_index = new_index;
    grp->grp_item  = new_item;
  }
  /* import_index, import_item, export_index, export_item */
  {
    int n_nb = mesh->n_neighbor_pe;
    int n_imp = mesh->import_index[n_nb] + n_add;
    int n_exp = mesh->export_index[n_nb];
    int *new_imp_index, *new_imp_item, *new_exp_index, *new_exp_item;
    int pos;
    for (i = 0; i < n_nb; i++) n_exp += HECMW_varray_int_nval(exp_added + i);
    new_imp_index = (int *)HECMW_malloc(sizeof(int) * (n_nb + 1));
    new_imp_item  = (int *)HECMW_malloc(sizeof(int) * (n_imp + 1));
    new_exp_index = (int *)HECMW_malloc(sizeof(int) * (n_nb + 1));
    new_exp_item  = (int *)HECMW_malloc(sizeof(int) * (n_exp + 1));
    if (new_imp_index == NULL || new_imp_item == NULL ||
        new_exp_index == NULL || new_exp_item == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    new_imp_index[0] = 0;
    for (pos = 0, i = 0; i < n_nb; i++) {
      for (j = mesh->import_index[i]; j < mesh->import_index[i + 1]; j++) {
        new_imp_item[pos++] = mesh->import_item[j];
      }
      for (j = 0; j < (int)HECMW_varray_int_nval(imp_added + i); j++) {
        new_imp_item[pos++] = HECMW_varray_int_get(imp_added + i, j);
      }
      new_imp_index[i + 1] = pos;
    }
    new_exp_index[0] = 0;
    for (pos = 0, i = 0; i < n_nb; i++) {
      for (j = mesh->export_index[i]; j < mesh->export_index[i + 1]; j++) {
        new_exp_item[pos++] = mesh->export_item[j];
      }
      for (j = 0; j < (int)HECMW_varray_int_nval(exp_added + i); j++) {
        new_exp_item[pos++] = HECMW_varray_int_get(exp_added + i, j);
      }
      new_exp_index[i + 1] = pos;
    }
    HECMW_free(mesh->import_index);
    HECMW_free(mesh->import_item);
    HECMW_free(mesh->export_index);
    HECMW_free(mesh->export_item);
    mesh->import_index = new_imp_index;
    mesh->import_item  = new_imp_item;
    mesh->export_index = new_exp_index;
    mesh->export_item  = new_exp_item;
  }
  /* node_old2new, node_new2old */
  if (mesh->node_old2new != NULL) {
    int *old2new = (int *)HECMW_malloc(sizeof(int) * ng_new);
    int *new2old = (int *)HECMW_malloc(sizeof(int) * ng_new);
    if (old2new == NULL || new2old == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (i = 0; i < ng_old; i++) {
      old2new[i] = SHIFTED(mesh->node_old2new[i] + 1) - 1;
    }
    for (k = 0; k < n_add; k++) old2new[ng_old + k] = nn_old + k;
    for (i = 0; i < ng_new; i++) new2old[old2new[i]] = i;
    HECMW_free(mesh->node_old2new);
    HECMW_free(mesh->node_new2old);
    mesh->node_old2new = old2new;
    mesh->node_new2old = new2old;
  }
#undef SHIFTED
  /* refine_origin: the added nodes have no origin */
  if (mesh->n_refine > 0) {
    struct hecmwST_refine_origin *reforg = mesh->refine_origin;
    int last  = reforg->index[mesh->n_refine];
    int count = reforg->item_index[last];
    reforg->item_index = (int *)HECMW_realloc(
        reforg->item_index, sizeof(int) * (last + n_add + 1));
    if (reforg->item_index == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    for (k = 1; k <= n_add; k++) reforg->item_index[last + k] = count;
    reforg->index[mesh->n_refine] = last + n_add;
    mesh->n_node_refine_hist[mesh->n_refine - 1] = ng_new;
  }
 
  mesh->n_node       = nn_old + n_add;
  mesh->nn_middle    = mesh->n_node;
  mesh->n_node_gross = ng_new;
  return HECMW_SUCCESS;
}
 
static int replicate_contact_slave_nodes(struct hecmwST_local_mesh *mesh) {
  const struct hecmwST_contact_pair *cp = mesh->contact_pair;
  int n_nb   = mesh->n_neighbor_pe;
  int n_pair = cp->n_pair;
  HECMW_Comm comm;
  HECMW_Request *requests;
  HECMW_Status *statuses;
  int *my_flag, *nb_flag, *len_send, *len_recv;
  int *memb_index, *memb_item;
  struct hecmw_varray_int *ibuf, *lid_added, *imp_added, *exp_added;
  double **dbuf;
  struct added_nodes added;
  int i, k, n_node_pair, tag, n_add, n_add_prev;
 
  if (n_nb == 0) return HECMW_SUCCESS;
  if (HECMW_partcontact_get_owner(mesh->hecmw_flag_partcontact) !=
      HECMW_FLAG_PARTCONTACT_OWNER_MASTER)
    return HECMW_SUCCESS;
  for (n_node_pair = 0, i = 0; i < n_pair; i++) {
    if (is_node_slave_pair(cp, i)) n_node_pair++;
  }
  if (n_node_pair == 0) return HECMW_SUCCESS;
 
  HECMW_log(HECMW_LOG_DEBUG, "Started replicating contact slave nodes...\n");
 
  comm     = HECMW_comm_get_comm();
  requests = (HECMW_Request *)HECMW_malloc(sizeof(HECMW_Request) * n_nb * 2);
  statuses = (HECMW_Status *)HECMW_malloc(sizeof(HECMW_Status) * n_nb * 2);
  my_flag  = (int *)HECMW_malloc(sizeof(int) * n_pair);
  nb_flag  = (int *)HECMW_malloc(sizeof(int) * n_pair * n_nb);
  len_send = (int *)HECMW_malloc(sizeof(int) * n_nb * 2);
  len_recv = (int *)HECMW_malloc(sizeof(int) * n_nb * 2);
  ibuf     = (struct hecmw_varray_int *)HECMW_malloc(
      sizeof(struct hecmw_varray_int) * n_nb * 2);
  dbuf = (double **)HECMW_malloc(sizeof(double *) * n_nb * 2);
  imp_added = (struct hecmw_varray_int *)HECMW_malloc(
      sizeof(struct hecmw_varray_int) * n_nb * 3);
  if (requests == NULL || statuses == NULL || my_flag == NULL ||
      nb_flag == NULL || len_send == NULL || len_recv == NULL ||
      ibuf == NULL || dbuf == NULL || imp_added == NULL) {
    HECMW_set_error(errno, "");
    return HECMW_ERROR;
  }
  lid_added = imp_added + n_nb;
  exp_added = imp_added + n_nb * 2;
  for (i = 0; i < n_nb * 3; i++) HECMW_varray_int_init(imp_added + i);
  for (i = 0; i < n_nb * 2; i++) HECMW_varray_int_init(ibuf + i);
 
  /* which pairs have master elements in the neighbors */
  for (i = 0; i < n_pair; i++) {
    my_flag[i] = is_node_slave_pair(cp, i) && has_internal_master(mesh, i);
  }
  tag = 0;
  for (i = 0; i < n_nb; i++) {
    HECMW_Isend(my_flag, n_pair, HECMW_INT, mesh->neighbor_pe[i], tag, comm,
                requests + i);
  }
  for (i = 0; i < n_nb; i++) {
    HECMW_Recv(nb_flag + n_pair * i, n_pair, HECMW_INT, mesh->neighbor_pe[i],
               tag, comm, statuses + i);
  }
  HECMW_Waitall(n_nb, requests, statuses);
 
  /* send the slave nodes to the neighbors owning master elements */
  if (get_internal_node_membership(mesh, &memb_index, &memb_item) !=
      HECMW_SUCCESS)
    return HECMW_ERROR;
  for (i = 0; i < n_nb; i++) {
    if (pack_contact_slave_nodes(mesh, nb_flag + n_pair * i, memb_index,
                                 memb_item, ibuf + i, dbuf + i,
                                 len_send + 2 * i + 1) != HECMW_SUCCESS)
      return HECMW_ERROR;
    len_send[2 * i] = HECMW_varray_int_nval(ibuf + i);
  }
  HECMW_free(memb_index);
  HECMW_free(memb_item);
  for (i = 0; i < n_nb; i++) {
    HECMW_Isend(len_send + 2 * i, 2, HECMW_INT, mesh->neighbor_pe[i], tag,
                comm, requests + i);
  }
  for (i = 0; i < n_nb; i++) {
    HECMW_Recv(len_recv + 2 * i, 2, HECMW_INT, mesh->neighbor_pe[i], tag, comm,
               statuses + i);
  }
  HECMW_Waitall(n_nb, requests, statuses);
  for (i = 0; i < n_nb; i++) {
    HECMW_Isend(HECMW_varray_int_get_v(ibuf + i), len_send[2 * i], HECMW_INT,
                mesh->neighbor_pe[i], tag, comm, requests + 2 * i);
    HECMW_Isend(dbuf[i], len_send[2 * i + 1], HECMW_DOUBLE,
                mesh->neighbor_pe[i], tag, comm, requests + 2 * i + 1);
  }
 
  /* receive the slave nodes and select those to add */
  HECMW_varray_int_init(&added.node_ID);
  HECMW_varray_int_init(&added.global_ID);
  HECMW_varray_int_init(&added.memb_index);
  HECMW_varray_int_init(&added.memb_item);
  HECMW_varray_int_init(&added.init_index);
  added.n_coord = 0;
  added.n_init  = 0;
  added.coord   = NULL;
  added.init    = NULL;
  for (i = 0; i < n_nb; i++) {
    HECMW_varray_int_resize(ibuf + n_nb + i, len_recv[2 * i]);
    dbuf[n_nb + i] =
        (double *)HECMW_malloc(sizeof(double) * (len_recv[2 * i + 1] + 1));
    if (dbuf[n_nb + i] == NULL) {
      HECMW_set_error(errno, "");
      return HECMW_ERROR;
    }
    HECMW_Recv(HECMW_varray_int_get_v(ibuf + n_nb + i), len_recv[2 * i],
               HECMW_INT, mesh->neighbor_pe[i], tag, comm, statuses);
    HECMW_Recv(dbuf[n_nb + i], len_recv[2 * i + 1], HECMW_DOUBLE,
               mesh->neighbor_pe[i], tag, comm, statuses);
    n_add_prev = HECMW_varray_int_nval(&added.global_ID);
    if (unpack_contact_slave_nodes(mesh, mesh->neighbor_pe[i],
                                   HECMW_varray_int_get_cv(ibuf + n_nb + i),
                                   dbuf[n_nb + i], len_recv[2 * i + 1],
                                   &added, lid_added + i) != HECMW_SUCCESS)
      return HECMW_ERROR;
    n_add = HECMW_varray_int_nval(&added.global_ID);
    for (k = n_add_prev; k < n_add; k++) {
      HECMW_varray_int_append(imp_added + i, mesh->n_node + k + 1);
    }
  }
  HECMW_Waitall(n_nb * 2, requests, statuses);
 
  /* tell the owners which nodes have been added, to export them */
  for (i = 0; i < n_nb; i++) {
    len_send[i] = HECMW_varray_int_nval(lid_added + i);
    HECMW_Isend(len_send + i, 1, HECMW_INT, mesh->neighbor_pe[i], tag, comm,
                requests + i);
  }
  for (i = 0; i < n_nb; i++) {
    HECMW_Recv(len_recv + i, 1, HECMW_INT, mesh->neighbor_pe[i], tag, comm,
               statuses + i);
  }
  HECMW_Waitall(n_nb, requests, statuses);
  for (i = 0; i < n_nb; i++) {
    HECMW_Isend(HECMW_varray_int_get_v(lid_added + i), len_send[i], HECMW_INT,
                mesh->neighbor_pe[i], tag, comm, requests + i);
  }
  for (i = 0; i < n_nb; i++) {
    HECMW_varray_int_resize(exp_added + i, len_recv[i]);
    HECMW_Recv(HECMW_varray_int_get_v(exp_added + i), len_recv[i], HECMW_INT,
               mesh->neighbor_pe[i], tag, comm, statuses + i);
  }
  HECMW_Waitall(n_nb, requests, statuses);
 
  if (insert_external_nodes(mesh, &added, imp_added, exp_added) !=
      HECMW_SUCCESS)
    return HECMW_ERROR;
 
  for (i = 0; i < n_nb * 2; i++) {
    HECMW_varray_int_finalize(ibuf + i);
    HECMW_free(dbuf[i]);
  }
  for (i = 0; i < n_nb * 3; i++) HECMW_varray_int_finalize(imp_added + i);
  HECMW_varray_int_finalize(&added.node_ID);
  HECMW_varray_int_finalize(&added.global_ID);
  HECMW_varray_int_finalize(&added.memb_index);
  HECMW_varray_int_finalize(&added.memb_item);
  HECMW_varray_int_finalize(&added.init_index);
  HECMW_free(added.coord);
  HECMW_free(added.init);
  HECMW_free(ibuf);
  HECMW_free(dbuf);
  HECMW_free(imp_added);
  HECMW_free(my_flag);
  HECMW_free(nb_flag);
  HECMW_free(len_send);
  HECMW_free(len_recv);
  HECMW_free(requests);
  HECMW_free(statuses);
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished replicating contact slave nodes.\n");
  return HECMW_SUCCESS;
}
 
/*============================================================================*/
/*                                                                            */
/*  refine distributed mesh ONCE                                              */
/*                                                                            */
/*============================================================================*/
static int refine_dist_mesh(struct hecmwST_local_mesh *mesh,
                            struct hecmwST_local_mesh *ref_mesh) {
  HECMW_log(HECMW_LOG_DEBUG, "Started refining mesh...\n");
 
  if (copy_unchanging_info(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (call_refiner(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (rebuild_info(mesh, ref_mesh) != HECMW_SUCCESS) return HECMW_ERROR;
 
  HECMW_log(HECMW_LOG_DEBUG, "Finished refining mesh.\n");
  return HECMW_SUCCESS;
}
 
/*============================================================================*/
/*                                                                            */
/*  refine HEC-MW distributed mesh data                                       */
/*                                                                            */
/*============================================================================*/
int HECMW_dist_refine(struct hecmwST_local_mesh **mesh, int refine,
                      const char *cad_filename, const char *part_filename) {
  int error_flag = HECMW_SUCCESS;
  int i;
 
  if (refine <= 0) return HECMW_SUCCESS;
 
  HECMW_log(HECMW_LOG_DEBUG, "Refinement requested; starting...\n");
 
  if ((*mesh)->n_refine > 0) {
    if (renumber_back_elements(*mesh) != HECMW_SUCCESS) return HECMW_ERROR;
    if (renumber_back_nodes(*mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  }
 
  if (prepare_refiner(*mesh, cad_filename, part_filename) != HECMW_SUCCESS)
    return HECMW_ERROR;
 
  for (i = 0; i < refine; i++) {
    struct hecmwST_local_mesh *ref_mesh = HECMW_dist_alloc();
    if (ref_mesh == NULL) {
      error_flag = HECMW_ERROR;
      break;
    }
 
    if (i > 0) {
      clear_refiner();
    }
 
    HECMW_log(HECMW_LOG_DEBUG, "Refining(%d)...\n", i + 1);
 
    if (refine_dist_mesh(*mesh, ref_mesh) != HECMW_SUCCESS) {
      HECMW_log(HECMW_LOG_ERROR, "Refinement failed\n");
      error_flag = HECMW_ERROR;
      break;
    }
 
    HECMW_dist_free(*mesh);
    *mesh = ref_mesh;
  }
 
  if (terminate_refiner(*mesh) != HECMW_SUCCESS) return HECMW_ERROR;
 
  if (error_flag != HECMW_SUCCESS) return HECMW_ERROR;
 
  if (renumber_nodes(*mesh) != HECMW_SUCCESS) return HECMW_ERROR;
  if (renumber_elements(*mesh) != HECMW_SUCCESS) return HECMW_ERROR;
 
  if (replicate_contact_slave_nodes(*mesh) != HECMW_SUCCESS) return HECMW_ERROR;
 
  HECMW_log(HECMW_LOG_DEBUG, "Refinement finished.\n");
 
  return HECMW_SUCCESS;
}
 
#else  /* HECMW_WITH_REFINER */
 
int HECMW_dist_refine(struct hecmwST_local_mesh **mesh, int refine,
                      const char *cad_filename, const char *part_filename) {
  if (refine > 0) {
    HECMW_log(HECMW_LOG_WARN, "Refiner not enabled; ignoring...\n");
  }
  return HECMW_SUCCESS;
}
#endif /* HECMW_WITH_REFINER */