Zoltan2_PamgenMeshAdapter.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//   Zoltan2: A package of combinatorial algorithms for scientific computing
//                  Copyright 2012 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Karen Devine      (kddevin@sandia.gov)
//                    Erik Boman        (egboman@sandia.gov)
//                    Siva Rajamanickam (srajama@sandia.gov)
//
// ***********************************************************************
//
// @HEADER

#ifndef _ZOLTAN2_PAMGENMESHADAPTER_HPP_
#define _ZOLTAN2_PAMGENMESHADAPTER_HPP_

#include <Zoltan2_MeshAdapter.hpp>
#include <Zoltan2_StridedData.hpp>
#include <vector>
#include <string>

#include <pamgen_im_exodusII.h>
#include "pamgen_im_ne_nemesisI.h"

namespace Zoltan2 {

template <typename User>
  class PamgenMeshAdapter: public MeshAdapter<User> {

public:

  typedef typename InputTraits<User>::scalar_t    scalar_t;
  typedef typename InputTraits<User>::lno_t    lno_t;
  typedef typename InputTraits<User>::gno_t    gno_t;
  typedef typename InputTraits<User>::part_t   part_t;
  typedef typename InputTraits<User>::node_t   node_t;
  typedef MeshAdapter<User>       base_adapter_t;
  typedef User user_t;
  typedef std::map<int, int> MapType;

  PamgenMeshAdapter(const Comm<int> &comm, std::string typestr="region");

  void print(int);

  // The MeshAdapter interface.
  // This is the interface that would be called by a model or a problem .

  bool areEntityIDsUnique(MeshEntityType etype) const {
    return etype==MESH_REGION;
  }

  size_t getLocalNumOf(MeshEntityType etype) const
  {
    if ((MESH_REGION == etype && 3 == dimension_) ||
  (MESH_FACE == etype && 2 == dimension_)) {
      return num_elem_;
    }

    if (MESH_VERTEX == etype) {
      return num_nodes_;
    }

    return 0;
  }
   
  void getIDsViewOf(MeshEntityType etype, const gno_t *&Ids) const
  {
    if ((MESH_REGION == etype && 3 == dimension_) ||
  (MESH_FACE == etype && 2 == dimension_)) {
      Ids = element_num_map_;
    }

    else if (MESH_VERTEX == etype) {
      Ids = node_num_map_;
    }

    else Ids = NULL;
  }
  
  void getTopologyViewOf(MeshEntityType etype,
       enum EntityTopologyType const *&Types) const {
    if ((MESH_REGION == etype && 3 == dimension_) ||
  (MESH_FACE == etype && 2 == dimension_)) {
      Types = elemTopology;
    }

    else if (MESH_VERTEX == etype) {
      Types = nodeTopology;
    }

    else Types = NULL;
  }

  void getWeightsViewOf(MeshEntityType etype, const scalar_t *&weights,
      int &stride, int idx = 0) const
  {
    weights = NULL;
    stride = 0;
  }

  int getDimension() const { return dimension_; }

  void getCoordinatesViewOf(MeshEntityType etype, const scalar_t *&coords,
          int &stride, int dim) const {
    if ((MESH_REGION == etype && 3 == dimension_) ||
         (MESH_FACE == etype && 2 == dimension_)) {
      if (dim == 0) {
  coords = Acoords_;
      } else if (dim == 1) {
  coords = Acoords_ + num_elem_;
      } else if (dim == 2) {
  coords = Acoords_ + 2 * num_elem_;
      }
      stride = 1;
    } else if (MESH_REGION == etype && 2 == dimension_) {
      coords = NULL;
      stride = 0;
    } else if (MESH_VERTEX == etype) {
      if (dim == 0) {
  coords = coords_;
      } else if (dim == 1) {
  coords = coords_ + num_nodes_;
      } else if (dim == 2) {
  coords = coords_ + 2 * num_nodes_;
      }
      stride = 1;
    } else {
      coords = NULL;
      stride = 0;
      Z2_THROW_NOT_IMPLEMENTED_IN_ADAPTER
    }
  }

  bool availAdjs(MeshEntityType source, MeshEntityType target) const {
    if ((MESH_REGION == source && MESH_VERTEX == target && 3 == dimension_) ||
  (MESH_FACE == source && MESH_VERTEX == target && 2 == dimension_) ||
  (MESH_VERTEX == source && MESH_REGION == target && 3 == dimension_) ||
  (MESH_VERTEX == source && MESH_FACE == target && 2 == dimension_)) {
      return TRUE;
    }

    return false;
  }

  size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
  {
    if ((MESH_REGION == source && MESH_VERTEX == target && 3 == dimension_) ||
  (MESH_FACE == source && MESH_VERTEX == target && 2 == dimension_)) {
      return tnoct_;
    }

    if ((MESH_VERTEX == source && MESH_REGION == target && 3 == dimension_) ||
  (MESH_VERTEX == source && MESH_FACE == target && 2 == dimension_)) {
      return telct_;
    }
    
    return 0;
  }

  void getAdjsView(MeshEntityType source, MeshEntityType target,
       const lno_t *&offsets, const gno_t *& adjacencyIds) const
  {
    if ((MESH_REGION == source && MESH_VERTEX == target && 3 == dimension_) ||
  (MESH_FACE == source && MESH_VERTEX == target && 2 == dimension_)) {
      offsets = elemOffsets_;
      adjacencyIds = (gno_t *)elemToNode_;
    } else if ((MESH_REGION==target && MESH_VERTEX==source && 3==dimension_) ||
         (MESH_FACE==target && MESH_VERTEX==source && 2==dimension_)) {
      offsets = nodeOffsets_;
      adjacencyIds = (gno_t *)nodeToElem_;
    } else if (MESH_REGION == source && 2 == dimension_) {
      offsets = NULL;
      adjacencyIds = NULL;
    } else {
      offsets = NULL;
      adjacencyIds = NULL;
      Z2_THROW_NOT_IMPLEMENTED_IN_ADAPTER
    }
  }

  //#define USE_MESH_ADAPTER
#ifndef USE_MESH_ADAPTER
  bool avail2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
  {
    if (through == MESH_VERTEX) {
      if (sourcetarget == MESH_REGION && dimension_ == 3) return true;
      if (sourcetarget == MESH_FACE && dimension_ == 2) return true;
    }
    if (sourcetarget == MESH_VERTEX) {
      if (through == MESH_REGION && dimension_ == 3) return true;
      if (through == MESH_FACE && dimension_ == 2) return true;
    }
    return false;
  }

  size_t getLocalNum2ndAdjs(MeshEntityType sourcetarget, 
          MeshEntityType through) const
  {
    if (through == MESH_VERTEX &&
  ((sourcetarget == MESH_REGION && dimension_ == 3) ||
   (sourcetarget == MESH_FACE && dimension_ == 2))) {
      return nEadj_;
    }

    if (sourcetarget == MESH_VERTEX &&
  ((through == MESH_REGION && dimension_ == 3) ||
   (through == MESH_FACE && dimension_ == 2))) {
      return nNadj_;
    }

    return 0;

  }

  void get2ndAdjsView(MeshEntityType sourcetarget, MeshEntityType through, 
          const lno_t *&offsets, const gno_t *&adjacencyIds) const
  {
    if (through == MESH_VERTEX &&
  ((sourcetarget == MESH_REGION && dimension_ == 3) ||
   (sourcetarget == MESH_FACE && dimension_ == 2))) {
      offsets = eStart_;
      adjacencyIds = eAdj_;
    } else if (sourcetarget == MESH_VERTEX &&
         ((through == MESH_REGION && dimension_ == 3) ||
    (through == MESH_FACE && dimension_ == 2))) {
      offsets = nStart_;
      adjacencyIds = nAdj_;
    } else {
      offsets = NULL;
      adjacencyIds = NULL;
      Z2_THROW_NOT_IMPLEMENTED_IN_ADAPTER
    }
  }
#endif

private:
  int dimension_, num_nodes_global_, num_elems_global_, num_nodes_, num_elem_;
  gno_t *element_num_map_, *node_num_map_;
  int *elemToNode_, tnoct_, *elemOffsets_;
  int *nodeToElem_, telct_, *nodeOffsets_;
  double *coords_, *Acoords_;
  lno_t *eStart_, *nStart_;
  gno_t *eAdj_, *nAdj_;
  size_t nEadj_, nNadj_;
  EntityTopologyType* nodeTopology;
  EntityTopologyType* elemTopology;

};

// Definitions

template <typename User>
PamgenMeshAdapter<User>::PamgenMeshAdapter(const Comm<int> &comm,
             std::string typestr):
  dimension_(0)
{
  int error = 0;
  char title[100];
  int exoid = 0;
  int num_elem_blk, num_node_sets, num_side_sets;
  error += im_ex_get_init(exoid, title, &dimension_,
        &num_nodes_, &num_elem_, &num_elem_blk,
        &num_node_sets, &num_side_sets);

  if (typestr.compare("region") == 0) {
    if (dimension_ == 3)
      this->setEntityTypes(typestr, "vertex", "vertex");
    else
      // automatically downgrade primary entity if problem is only 2D
      this->setEntityTypes("face", "vertex", "vertex");
  }
  else if (typestr.compare("vertex") == 0) {
    if (dimension_ == 3)
      this->setEntityTypes(typestr, "region", "region");
    else 
      this->setEntityTypes(typestr, "face", "face");
  }
  else {
    Z2_THROW_NOT_IMPLEMENTED_IN_ADAPTER
  }

  coords_ = new double [num_nodes_ * dimension_];

  error += im_ex_get_coord(exoid, coords_, coords_ + num_nodes_,
         coords_ + 2 * num_nodes_);
  
  element_num_map_ = new gno_t [num_elem_];
  std::vector<int> tmp;
  tmp.resize(num_elem_);
  
  // BDD cast to int did not always work!
  // error += im_ex_get_elem_num_map(exoid, (int *)element_num_map_)
  // This may be a case of calling the wrong method
  error += im_ex_get_elem_num_map(exoid, &tmp[0]);
  for(size_t i = 0; i < tmp.size(); i++)
    element_num_map_[i] = static_cast<gno_t>(tmp[i]);
    
  tmp.clear();
  tmp.resize(num_nodes_);
  node_num_map_ = new gno_t [num_nodes_];
  
  // BDD cast to int did not always work!
  // error += im_ex_get_node_num_map(exoid, (int *)node_num_map_);
  // This may be a case of calling the wrong method
  error += im_ex_get_node_num_map(exoid, &tmp[0]);
  for(size_t i = 0; i < tmp.size(); i++)
    node_num_map_[i] = static_cast<gno_t>(tmp[i]);
  
  nodeTopology = new enum EntityTopologyType[num_nodes_];
  for (int i=0;i<num_nodes_;i++)
    nodeTopology[i] = POINT;
  elemTopology = new enum EntityTopologyType[num_elem_];
  for (int i=0;i<num_elem_;i++) {
    if (dimension_==2)
      elemTopology[i] = QUADRILATERAL;
    else
      elemTopology[i] = HEXAHEDRON;
  }
  
  int *elem_blk_ids       = new int [num_elem_blk];
  error += im_ex_get_elem_blk_ids(exoid, elem_blk_ids);

  int *num_nodes_per_elem = new int [num_elem_blk];
  int *num_attr           = new int [num_elem_blk];
  int *num_elem_this_blk  = new int [num_elem_blk];
  char **elem_type        = new char * [num_elem_blk];
  int **connect           = new int * [num_elem_blk];

  for(int i = 0; i < num_elem_blk; i++){
    elem_type[i] = new char [MAX_STR_LENGTH + 1];
    error += im_ex_get_elem_block(exoid, elem_blk_ids[i], elem_type[i],
          (int*)&(num_elem_this_blk[i]),
          (int*)&(num_nodes_per_elem[i]),
          (int*)&(num_attr[i]));
    delete[] elem_type[i];
  }

  delete[] elem_type;
  elem_type = NULL;
  delete[] num_attr;
  num_attr = NULL;
  Acoords_ = new double [num_elem_ * dimension_];
  int a = 0;
  std::vector<std::vector<int> > sur_elem;
  sur_elem.resize(num_nodes_);

  for(int b = 0; b < num_elem_blk; b++) {
    connect[b] = new int [num_nodes_per_elem[b]*num_elem_this_blk[b]];
    error += im_ex_get_elem_conn(exoid, elem_blk_ids[b], connect[b]);

    for(int i = 0; i < num_elem_this_blk[b]; i++) {
      Acoords_[a] = 0;
      Acoords_[num_elem_ + a] = 0;

      if (3 == dimension_) {
  Acoords_[2 * num_elem_ + a] = 0;
      }

      for(int j = 0; j < num_nodes_per_elem[b]; j++) {
  int node = connect[b][i * num_nodes_per_elem[b] + j] - 1;
  Acoords_[a] += coords_[node];
  Acoords_[num_elem_ + a] += coords_[num_nodes_ + node];

  if(3 == dimension_) {
    Acoords_[2 * num_elem_ + a] += coords_[2 * num_nodes_ + node];
  }

  /*
   * in the case of degenerate elements, where a node can be
   * entered into the connect table twice, need to check to
   * make sure that this element is not already listed as
   * surrounding this node
   */
  if (sur_elem[node].empty() ||
      element_num_map_[a] != sur_elem[node][sur_elem[node].size()-1]) {
    /* Add the element to the list */
    sur_elem[node].push_back(element_num_map_[a]);
  }
      }

      Acoords_[a] /= num_nodes_per_elem[b];
      Acoords_[num_elem_ + a] /= num_nodes_per_elem[b];

      if(3 == dimension_) {
  Acoords_[2 * num_elem_ + a] /= num_nodes_per_elem[b];
      }

      a++;
    }

  }

  delete[] elem_blk_ids;
  elem_blk_ids = NULL;
  int nnodes_per_elem = num_nodes_per_elem[0];
  elemToNode_ = new int [num_elem_ * nnodes_per_elem];
  int telct = 0;
  elemOffsets_ = new int [num_elem_+1];
  tnoct_ = 0;
  int **reconnect = new int * [num_elem_];
  size_t max_nsur = 0;

  for (int b = 0; b < num_elem_blk; b++) {
    for (int i = 0; i < num_elem_this_blk[b]; i++) {
      elemOffsets_[telct] = tnoct_;
      reconnect[telct] = new int [num_nodes_per_elem[b]];

      for (int j = 0; j < num_nodes_per_elem[b]; j++) {
  elemToNode_[tnoct_]=
    node_num_map_[connect[b][i*num_nodes_per_elem[b] + j]-1];
  reconnect[telct][j] = connect[b][i*num_nodes_per_elem[b] + j];
  ++tnoct_;
      }

      ++telct;
    }
  }
  elemOffsets_[telct] = tnoct_;

  delete[] num_nodes_per_elem;
  num_nodes_per_elem = NULL;
  delete[] num_elem_this_blk;
  num_elem_this_blk = NULL;

  for(int b = 0; b < num_elem_blk; b++) {
    delete[] connect[b];
  }

  delete[] connect;
  connect = NULL;

  int max_side_nodes = nnodes_per_elem;
  int *side_nodes = new int [max_side_nodes];
  int *mirror_nodes = new int [max_side_nodes];

  /* Allocate memory necessary for the adjacency */
  eStart_ = new lno_t [num_elem_+1];
  nStart_ = new lno_t [num_nodes_+1];
  std::vector<int> eAdj;
  std::vector<int> nAdj;

  for (int i=0; i < max_side_nodes; i++) {
    side_nodes[i]=-999;
    mirror_nodes[i]=-999;
  }

  /* Find the adjacency for a nodal based decomposition */
  nEadj_ = 0;
  nNadj_ = 0;
  for(int ncnt=0; ncnt < num_nodes_; ncnt++) {
    if(sur_elem[ncnt].empty()) {
      std::cout << "WARNING: Node = " << ncnt+1 << " has no elements"
                << std::endl;
    } else {
      size_t nsur = sur_elem[ncnt].size();
      if (nsur > max_nsur)
  max_nsur = nsur;
    }
  }

  nodeToElem_ = new int [num_nodes_ * max_nsur];
  nodeOffsets_ = new int [num_nodes_+1];
  telct_ = 0;

  for (int ncnt = 0; ncnt < num_nodes_; ncnt++) {
    nodeOffsets_[ncnt] = telct_;
    nStart_[ncnt] = nNadj_;
    MapType nAdjMap;

    for (size_t i = 0; i < sur_elem[ncnt].size(); i++) {
      nodeToElem_[telct_] = sur_elem[ncnt][i];
      ++telct_;

#ifndef USE_MESH_ADAPTER
      for(int ecnt = 0; ecnt < num_elem_; ecnt++) {
  if (element_num_map_[ecnt] == sur_elem[ncnt][i]) {
    for (int j = 0; j < nnodes_per_elem; j++) {
      MapType::iterator iter =
        nAdjMap.find(elemToNode_[elemOffsets_[ecnt]+j]);

      if (node_num_map_[ncnt] != elemToNode_[elemOffsets_[ecnt]+j] &&
    iter == nAdjMap.end() ) {
        nAdj.push_back(elemToNode_[elemOffsets_[ecnt]+j]);
        nNadj_++;
        nAdjMap.insert({elemToNode_[elemOffsets_[ecnt]+j],
        elemToNode_[elemOffsets_[ecnt]+j]});
      }
    }

    break;
  }
      }
#endif
    }

    nAdjMap.clear();
  }

  nodeOffsets_[num_nodes_] = telct_;
  nStart_[num_nodes_] = nNadj_;

  nAdj_ = new gno_t [nNadj_];

  for (size_t i=0; i < nNadj_; i++) {
    nAdj_[i] = nAdj[i];
  }

  int nprocs = comm.getSize();
  //if (nprocs > 1) {
    int neid=0,num_elem_blks_global,num_node_sets_global,num_side_sets_global;
    error += im_ne_get_init_global(neid,&num_nodes_global_,&num_elems_global_,
           &num_elem_blks_global,&num_node_sets_global,
           &num_side_sets_global);

    int num_internal_nodes, num_border_nodes, num_external_nodes;
    int num_internal_elems, num_border_elems, num_node_cmaps, num_elem_cmaps;
    int proc = 0;
    error += im_ne_get_loadbal_param(neid, &num_internal_nodes,
             &num_border_nodes, &num_external_nodes,
             &num_internal_elems, &num_border_elems,
             &num_node_cmaps, &num_elem_cmaps, proc);

    int *node_cmap_ids = new int [num_node_cmaps];
    int *node_cmap_node_cnts = new int [num_node_cmaps];
    int *elem_cmap_ids = new int [num_elem_cmaps];
    int *elem_cmap_elem_cnts = new int [num_elem_cmaps];
    error += im_ne_get_cmap_params(neid, node_cmap_ids, node_cmap_node_cnts,
           elem_cmap_ids, elem_cmap_elem_cnts, proc);
    delete[] elem_cmap_ids;
    elem_cmap_ids = NULL;
    delete[] elem_cmap_elem_cnts;
    elem_cmap_elem_cnts = NULL;

    int **node_ids = new int * [num_node_cmaps];
    int **node_proc_ids = new int * [num_node_cmaps];
    for(int j = 0; j < num_node_cmaps; j++) {
      node_ids[j] = new int [node_cmap_node_cnts[j]];
      node_proc_ids[j] = new int [node_cmap_node_cnts[j]];
      error += im_ne_get_node_cmap(neid, node_cmap_ids[j], node_ids[j],
           node_proc_ids[j], proc);
    }
    delete[] node_cmap_ids;
    node_cmap_ids = NULL;
    int *sendCount = new int [nprocs];
    int *recvCount = new int [nprocs];

    // Post receives
    RCP<CommRequest<int> >*requests=new RCP<CommRequest<int> >[num_node_cmaps];
    for (int cnt = 0, i = 0; i < num_node_cmaps; i++) {
      try {
  requests[cnt++] =
    Teuchos::ireceive<int,int>(comm,
             rcp(&(recvCount[node_proc_ids[i][0]]),
           false),
             node_proc_ids[i][0]);
      }
      Z2_FORWARD_EXCEPTIONS;
    }

    Teuchos::barrier<int>(comm);
    size_t totalsend = 0;

    for(int j = 0; j < num_node_cmaps; j++) {
      sendCount[node_proc_ids[j][0]] = 1;
      for(int i = 0; i < node_cmap_node_cnts[j]; i++) {
  sendCount[node_proc_ids[j][i]] += sur_elem[node_ids[j][i]-1].size()+2;
      }
      totalsend += sendCount[node_proc_ids[j][0]];
    }

    // Send data; can use readySend since receives are posted.
    for (int i = 0; i < num_node_cmaps; i++) {
      try {
  Teuchos::readySend<int,int>(comm, sendCount[node_proc_ids[i][0]],
            node_proc_ids[i][0]);
      }
      Z2_FORWARD_EXCEPTIONS;
    }

    // Wait for messages to return.
    try {
      Teuchos::waitAll<int>(comm, arrayView(requests, num_node_cmaps));
    }
    Z2_FORWARD_EXCEPTIONS;

    delete [] requests;

    // Allocate the receive buffer.
    size_t totalrecv = 0;
    int maxMsg = 0;
    int nrecvranks = 0;
    for(int i = 0; i < num_node_cmaps; i++) {
      if (recvCount[node_proc_ids[i][0]] > 0) {
  totalrecv += recvCount[node_proc_ids[i][0]];
  nrecvranks++;
  if (recvCount[node_proc_ids[i][0]] > maxMsg)
    maxMsg = recvCount[node_proc_ids[i][0]];
      }
    }

    int *rbuf = NULL;
    if (totalrecv) rbuf = new int[totalrecv];

    requests = new RCP<CommRequest<int> > [nrecvranks];

    // Error checking for memory and message size.
    int OK[2] = {1,1};
    // OK[0] -- true/false indicating whether each message size fits in an int
    //          (for MPI).
    // OK[1] -- true/false indicating whether memory allocs are OK
    int gOK[2]; // For global reduce of OK.

    if (size_t(maxMsg) * sizeof(int) > INT_MAX) OK[0] = false;
    if (totalrecv && !rbuf) OK[1] = 0;
    if (!requests) OK[1] = 0;

    // Post receives

    size_t offset = 0;

    if (OK[0] && OK[1]) {
      int rcnt = 0;
      for (int i = 0; i < num_node_cmaps; i++) {
  if (recvCount[node_proc_ids[i][0]]) {
    try {
      requests[rcnt++] =
        Teuchos::
        ireceive<int,int>(comm,
        Teuchos::arcp(&rbuf[offset], 0,
                recvCount[node_proc_ids[i][0]],
                false),
        node_proc_ids[i][0]);
    }
    Z2_FORWARD_EXCEPTIONS;
  }
  offset += recvCount[node_proc_ids[i][0]];
      }
    }

    delete[] recvCount;

    // Use barrier for error checking
    Teuchos::reduceAll<int>(comm, Teuchos::REDUCE_MIN, 2, OK, gOK);
    if (!gOK[0] || !gOK[1]) {
      delete [] rbuf;
      delete [] requests;
      if (!gOK[0])
  throw std::runtime_error("Max single message length exceeded");
      else
  throw std::bad_alloc();
    }

    int *sbuf = NULL;
    if (totalsend) sbuf = new int[totalsend];
    a = 0;

    for(int j = 0; j < num_node_cmaps; j++) {
      sbuf[a++] = node_cmap_node_cnts[j];
      for(int i = 0; i < node_cmap_node_cnts[j]; i++) {
  sbuf[a++] = node_num_map_[node_ids[j][i]-1];
  sbuf[a++] = sur_elem[node_ids[j][i]-1].size();
  for(size_t ecnt=0; ecnt < sur_elem[node_ids[j][i]-1].size(); ecnt++) {
    sbuf[a++] = sur_elem[node_ids[j][i]-1][ecnt];
  }
      }
    }

    delete[] node_cmap_node_cnts;
    node_cmap_node_cnts = NULL;

    for(int j = 0; j < num_node_cmaps; j++) {
      delete[] node_ids[j];
    }

    delete[] node_ids;
    node_ids = NULL;
    ArrayRCP<int> sendBuf;

    if (totalsend)
      sendBuf = ArrayRCP<int>(sbuf, 0, totalsend, true);
    else
      sendBuf = Teuchos::null;

    // Send data; can use readySend since receives are posted.
    offset = 0;
    for (int i = 0; i < num_node_cmaps; i++) {
      if (sendCount[node_proc_ids[i][0]]) {
  try{
    Teuchos::readySend<int,
      int>(comm,
     Teuchos::arrayView(&sendBuf[offset],
            sendCount[node_proc_ids[i][0]]),
     node_proc_ids[i][0]);
  }
  Z2_FORWARD_EXCEPTIONS;
      }
      offset += sendCount[node_proc_ids[i][0]];
    }

    for(int j = 0; j < num_node_cmaps; j++) {
      delete[] node_proc_ids[j];
    }

    delete[] node_proc_ids;
    node_proc_ids = NULL;
    delete[] sendCount;

    // Wait for messages to return.
    try{
      Teuchos::waitAll<int>(comm, Teuchos::arrayView(requests, nrecvranks));
    }
    Z2_FORWARD_EXCEPTIONS;

    delete[] requests;
    a = 0;

    for (int i = 0; i < num_node_cmaps; i++) {
      int num_nodes_this_processor = rbuf[a++];

      for (int j = 0; j < num_nodes_this_processor; j++) {
  int this_node = rbuf[a++];
  int num_elem_this_node = rbuf[a++];

  for (int ncnt = 0; ncnt < num_nodes_; ncnt++) {
    if (node_num_map_[ncnt] == this_node) {
      for (int ecnt = 0; ecnt < num_elem_this_node; ecnt++) {
        sur_elem[ncnt].push_back(rbuf[a++]);
      }

      break;
    }
  }
      }
    }

    delete[] rbuf;
    //}

#ifndef USE_MESH_ADAPTER
  for(int ecnt=0; ecnt < num_elem_; ecnt++) {
    eStart_[ecnt] = nEadj_;
    MapType eAdjMap;
    int nnodes = nnodes_per_elem;
    for(int ncnt=0; ncnt < nnodes; ncnt++) {
      int node = reconnect[ecnt][ncnt]-1;
      for(size_t i=0; i < sur_elem[node].size(); i++) {
  int entry = sur_elem[node][i];
  MapType::iterator iter = eAdjMap.find(entry);

  if(element_num_map_[ecnt] != entry &&
     iter == eAdjMap.end()) {
    eAdj.push_back(entry);
    nEadj_++;
    eAdjMap.insert({entry, entry});
  }
      }
    }

    eAdjMap.clear();
  }
#endif

  for(int b = 0; b < num_elem_; b++) {
    delete[] reconnect[b];
  }

  delete[] reconnect;
  reconnect = NULL;
  eStart_[num_elem_] = nEadj_;

  eAdj_ = new gno_t [nEadj_];

  for (size_t i=0; i < nEadj_; i++) {
    eAdj_[i] = eAdj[i];
  }

  delete[] side_nodes;
  side_nodes = NULL;
  delete[] mirror_nodes;
  mirror_nodes = NULL;
}

template <typename User>
void PamgenMeshAdapter<User>::print(int me)
{
  std::string fn(" PamgenMesh ");
  std::cout << me << fn
            << " dim = " << dimension_
            << " nodes = " << num_nodes_
            << " nelems = " << num_elem_
            << std::endl;

  for (int i = 0; i < num_elem_; i++) {
    std::cout << me << fn << i 
              << " Elem " << element_num_map_[i]
              << " Coords: ";
    for (int j = 0; j < dimension_; j++)
      std::cout << Acoords_[i + j * num_elem_] << " ";
    std::cout << std::endl;
  }

#ifndef USE_MESH_ADAPTER
  for (int i = 0; i < num_elem_; i++) {
    std::cout << me << fn << i 
              << " Elem " << element_num_map_[i]
              << " Graph: ";
    for (int j = eStart_[i]; j < eStart_[i+1]; j++)
      std::cout << eAdj_[j] << " ";
    std::cout << std::endl;
  }
#endif
}
  
}  //namespace Zoltan2
  
#endif