tools/mbcoupler/Coupler.cpp

#include "Coupler.hpp"
#include "moab/ParallelComm.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "ElemUtil.hpp"
#include "moab/CN.hpp"
#include "moab/gs.hpp"
#include "moab/TupleList.hpp"
#include "moab/Error.hpp"

/* C++ includes */
#include <cstdio>
#include <cassert>
#include <iostream>
#include <algorithm>
#include <sstream>

#define ERROR(a) {if (MB_SUCCESS != err) std::cerr << a << std::endl;}
#define ERRORR(a,b) {if (MB_SUCCESS != b) {std::cerr << a << std::endl; return b;}}
#define ERRORMPI(a,b) {if (MPI_SUCCESS != b) {std::cerr << a << std::endl; return MB_FAILURE;}}

#define MASTER_PROC 0

namespace moab {

bool debug = false;
int pack_tuples(TupleList *tl, void **ptr);
void unpack_tuples(void *ptr, TupleList **tlp);

Coupler::Coupler(Interface *impl,
                 ParallelComm *pc,
                 Range &local_elems,
                 int coupler_id,
                 bool init_tree,
                 int max_ent_dim)
  : mbImpl(impl), myPc(pc), myId(coupler_id), numIts(3), max_dim(max_ent_dim), _ntot(0), spherical(false)
{
  assert(NULL != impl && (pc || !local_elems.empty()));

  // Keep track of the local points, at least for now
  myRange = local_elems;
  myTree = NULL;

  // Now initialize the tree
  if (init_tree)
    initialize_tree();

  // Initialize tuple lists to indicate not initialized
  mappedPts = NULL;
  targetPts = NULL;
  _spectralSource = _spectralTarget = NULL;
}

  /* Destructor
   */
Coupler::~Coupler()
{
  // This will clear the cache
  delete (moab::Element::SpectralHex*)_spectralSource;
  delete (moab::Element::SpectralHex*)_spectralTarget;
  delete myTree;
  delete targetPts;
  delete mappedPts;
}

ErrorCode Coupler::initialize_tree()
{
  Range local_ents;

  // Get entities on the local part
  ErrorCode result = MB_SUCCESS;
  if (myPc)
  {
    result = myPc->get_part_entities(local_ents, max_dim);
    if (local_ents.empty())
    {
      max_dim--;
      result = myPc->get_part_entities(local_ents, max_dim);// go one dimension lower
      // right now, this is used for spherical meshes only
    }
  }
  else local_ents = myRange;
  if (MB_SUCCESS != result || local_ents.empty()) {
    std::cout << "Problems getting source entities" << std::endl;
    return result;
  }

  // Build the tree for local processor
  int max_per_leaf = 6;
  for (int i = 0; i < numIts; i++) {
    std::ostringstream str;
    str << "PLANE_SET=0;"
        << "MAX_PER_LEAF=" << max_per_leaf << ";";
    if (spherical && !local_ents.empty())
    {
      // get coordinates of one vertex, and use for radius computation
      EntityHandle elem= local_ents[0];
      const EntityHandle * conn;
      int numn=0;
      mbImpl->get_connectivity(elem, conn, numn);
      CartVect pos0;
      mbImpl -> get_coords(conn, 1, &(pos0[0]) );
      double radius = pos0.length();
      str<<"SPHERICAL=true;RADIUS="<<radius<<";";
    }
    FileOptions opts(str.str().c_str());
    myTree = new AdaptiveKDTree(mbImpl);
    result = myTree->build_tree(local_ents, &localRoot, &opts);
    if (MB_SUCCESS != result) {
      std::cout << "Problems building tree";
      if (numIts != i) {
        delete myTree;
        max_per_leaf *= 2;
        std::cout << "; increasing elements/leaf to "
                  << max_per_leaf << std::endl;
      }
      else {
        std::cout << "; exiting" << std::endl;
        return result;
      }
    }
    else
      break; // Get out of tree building
  }

  // Get the bounding box for local tree
  if (myPc)
    allBoxes.resize(6*myPc->proc_config().proc_size());
  else
    allBoxes.resize(6);

  unsigned int my_rank = (myPc ? myPc->proc_config().proc_rank() : 0);
  BoundBox box;
  result = myTree->get_bounding_box(box, &localRoot);
  if (MB_SUCCESS != result)
    return result;
  box.bMin.get(&allBoxes[6*my_rank]);
  box.bMax.get(&allBoxes[6*my_rank + 3]);

  // Now communicate to get all boxes
  if (myPc) {
    int mpi_err;
#if (MPI_VERSION >= 2)
    // Use "in place" option
    mpi_err = MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL,
                            &allBoxes[0], 6, MPI_DOUBLE,
                            myPc->proc_config().proc_comm());
#else
    {
      std::vector<double> allBoxes_tmp(6*myPc->proc_config().proc_size());
      mpi_err = MPI_Allgather(&allBoxes[6*my_rank], 6, MPI_DOUBLE,
                              &allBoxes_tmp[0], 6, MPI_DOUBLE,
                              myPc->proc_config().proc_comm());
      allBoxes = allBoxes_tmp;
    }
#endif
    if (MPI_SUCCESS != mpi_err)
      return MB_FAILURE;
  }

  /*  std::ostringstream blah;
  for(int i = 0; i < allBoxes.size(); i++)
  blah << allBoxes[i] << " ";
  std::cout << blah.str() << "\n";*/

#ifdef VERBOSE
  double min[3] = {0, 0, 0}, max[3] = {0, 0, 0};
  unsigned int dep;
  myTree->get_info(localRoot, min, max, dep);
  std::cout << "Proc " << my_rank << ": box min/max, tree depth = ("
            << min[0] << "," << min[1] << "," << min[2] << "), ("
            << max[0] << "," << max[1] << "," << max[2] << "), "
            << dep << std::endl;
#endif

  return result;
}

ErrorCode Coupler::initialize_spectral_elements(EntityHandle rootSource, EntityHandle rootTarget,
                                                bool &specSou, bool &specTar)
{
  /*void * _spectralSource;
    void * _spectralTarget;*/

  moab::Range spectral_sets;
  moab::Tag  sem_tag;
  int sem_dims[3];
  ErrorCode rval = mbImpl->tag_get_handle("SEM_DIMS", 3, moab::MB_TYPE_INTEGER, sem_tag);
  if (moab::MB_SUCCESS != rval) {
    std::cout << "Can't find tag, no spectral set\n";
    return MB_SUCCESS; // Nothing to do, no spectral elements
  }
  rval = mbImpl->get_entities_by_type_and_tag(rootSource, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets);
  if (moab::MB_SUCCESS != rval || spectral_sets.empty())
    std::cout << "Can't get sem set on source\n";
  else {
    moab::EntityHandle firstSemSet=spectral_sets[0];
    rval = mbImpl->tag_get_data(sem_tag, &firstSemSet, 1, (void*)sem_dims);
    if (moab::MB_SUCCESS != rval)
      return MB_FAILURE;

    if (sem_dims[0]!=sem_dims[1] || sem_dims[0] != sem_dims[2]) {
      std::cout << " dimensions are different. bail out\n";
      return MB_FAILURE;
    }
    // Repeat for target sets
    spectral_sets.empty();
    // Now initialize a source spectral element !
    _spectralSource = new moab::Element::SpectralHex(sem_dims[0]);
    specSou = true;
  }
  // Repeat for target source
  rval = mbImpl->get_entities_by_type_and_tag(rootTarget, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets);
  if (moab::MB_SUCCESS != rval || spectral_sets.empty())
    std::cout << "Can't get sem set on target\n";
  else {
    moab::EntityHandle firstSemSet=spectral_sets[0];
    rval = mbImpl->tag_get_data(sem_tag, &firstSemSet, 1, (void*)sem_dims);
    if (moab::MB_SUCCESS != rval)
      return MB_FAILURE;

    if (sem_dims[0]!=sem_dims[1] || sem_dims[0] != sem_dims[2]) {
      std::cout << " dimensions are different. bail out\n";
      return MB_FAILURE;
    }
    // Repeat for target sets
    spectral_sets.empty();
    // Now initialize a target spectral element !
    _spectralTarget = new moab::Element::SpectralHex(sem_dims[0]);
    specTar = true;
  }

  _ntot = sem_dims[0]*sem_dims[1]*sem_dims[2];
  rval = mbImpl->tag_get_handle("SEM_X", _ntot, moab::MB_TYPE_DOUBLE, _xm1Tag);
  if (moab::MB_SUCCESS != rval) {
     std::cout << "Can't get xm1tag \n";
     return MB_FAILURE;
  }
  rval = mbImpl->tag_get_handle("SEM_Y", _ntot, moab::MB_TYPE_DOUBLE, _ym1Tag);
  if (moab::MB_SUCCESS != rval) {
     std::cout << "Can't get ym1tag \n";
     return MB_FAILURE;
  }
  rval = mbImpl->tag_get_handle("SEM_Z", _ntot, moab::MB_TYPE_DOUBLE, _zm1Tag);
  if (moab::MB_SUCCESS != rval) {
     std::cout << "Can't get zm1tag \n";
     return MB_FAILURE;
  }

  return MB_SUCCESS;
}

ErrorCode Coupler::locate_points(Range &targ_ents,
                                 double rel_eps,
                                 double abs_eps,
                                 TupleList *tl,
                                 bool store_local)
{
  // Get locations
  std::vector<double> locs(3*targ_ents.size());
  Range verts = targ_ents.subset_by_type(MBVERTEX);
  ErrorCode rval = mbImpl->get_coords(verts, &locs[0]);
  if (MB_SUCCESS != rval)
    return rval;
  // Now get other ents; reuse verts
  unsigned int num_verts = verts.size();
  verts = subtract(targ_ents, verts);
  // Compute centroids
  std::vector<EntityHandle> dum_conn(CN::MAX_NODES_PER_ELEMENT);
  std::vector<double> dum_pos(CN::MAX_NODES_PER_ELEMENT);
  const EntityHandle *conn;
  int num_conn;
  double *coords = &locs[num_verts];
  // Do this here instead of a function to allow reuse of dum_pos and dum_conn
  for (Range::const_iterator rit = verts.begin(); rit != verts.end(); ++rit) {
    rval = mbImpl->get_connectivity(*rit, conn, num_conn, false, &dum_conn);
    if (MB_SUCCESS != rval)
      return rval;
    rval = mbImpl->get_coords(conn, num_conn, &dum_pos[0]);
    if (MB_SUCCESS != rval)
      return rval;
    coords[0] = coords[1] = coords[2] = 0.0;
    for (int i = 0; i < num_conn; i++) {
      coords[0] += dum_pos[3*i];
      coords[1] += dum_pos[3*i + 1];
      coords[2] += dum_pos[3*i + 2];
    }
    coords[0] /= num_conn;
    coords[1] /= num_conn;
    coords[2] /= num_conn;
    coords += 3;
  }

  if (store_local)
    targetEnts = targ_ents;

  return locate_points(&locs[0], targ_ents.size(), rel_eps, abs_eps, tl, store_local);
}

ErrorCode Coupler::locate_points(double *xyz, unsigned int num_points,
                                 double rel_eps,
                                 double abs_eps,
                                 TupleList *tl,
                                 bool store_local)
{
  assert(tl || store_local);

  // target_pts: TL(to_proc, tgt_index, x, y, z): tuples sent to source mesh procs representing pts to be located
  // source_pts: TL(from_proc, tgt_index, src_index): results of source mesh proc point location, ready to send
  //             back to tgt procs; src_index of -1 indicates point not located (arguably not useful...)
  TupleList target_pts;
  target_pts.initialize(2, 0, 0, 3, num_points);
  target_pts.enableWriteAccess();

  TupleList source_pts;
  mappedPts = new TupleList(0, 0, 1, 3, target_pts.get_max());
  mappedPts->enableWriteAccess();

  source_pts.initialize(3, 0, 0, 0, target_pts.get_max());
  source_pts.enableWriteAccess();

  mappedPts->set_n(0);
  source_pts.set_n(0);
  ErrorCode result;

  unsigned int my_rank = (myPc ? myPc->proc_config().proc_rank() : 0);
  bool point_located;

  // For each point, find box(es) containing the point,
  // appending results to tuple_list;
  // keep local points separately, in local_pts, which has pairs
  // of <local_index, mapped_index>, where mapped_index is the index
  // into the mappedPts tuple list
  for (unsigned int i = 0; i < 3*num_points; i += 3) {

    std::vector<int>  procs_to_send_to;
    for (unsigned int j = 0; j < (myPc ? myPc->proc_config().proc_size() : 0); j++) {
      // Test if point is in proc's box
      if ((allBoxes[6*j] <= xyz[i] + abs_eps) && (xyz[i] <= allBoxes[6*j + 3] + abs_eps) &&
          (allBoxes[6*j + 1] <= xyz[i + 1] + abs_eps) && (xyz[i + 1] <= allBoxes[6*j + 4] + abs_eps) &&
          (allBoxes[6*j + 2] <= xyz[i + 2] + abs_eps) && (xyz[i + 2] <= allBoxes[6*j + 5] + abs_eps)) {
        // If in this proc's box, will send to proc to test further
        procs_to_send_to.push_back(j);
      }
    }
    if (procs_to_send_to.empty())
    {
#ifdef VERBOSE
      std::cout << " point index " << i/3 << ": " << xyz[i] << " " << xyz[i+1] << " " << xyz[i+2] << " not found in any box\n";
#endif
      // try to find the closest box, and put it in that box, anyway
      double min_dist = 1.e+20;
      int index = -1;
      for (unsigned int j = 0; j < (myPc ? myPc->proc_config().proc_size() : 0); j++)
      {
        BoundBox box( &allBoxes[6*j]); // form back the box
        double distance = box.distance(&xyz[i]); // will compute the distance in 3d, from the box
        if (distance < min_dist)
        {
          index = j;
          min_dist = distance;
        }
      }
      if (index ==-1)
      {
        // need to abort early, nothing we can do
        assert("cannot locate any box for some points");
        // need a better exit strategy
      }
#ifdef VERBOSE
      std::cout << " point index " << i/3 << " added to box for proc j:" << index << "\n";
#endif
      procs_to_send_to.push_back(index); // will send to just one proc, that has the closest box
    }
    // we finally decided to populate the tuple list for a list of processors
    for (size_t k=0; k<procs_to_send_to.size(); k++)
    {
      unsigned int j = procs_to_send_to[k];
      // Check size of tuple list, grow if we're at max
      if (target_pts.get_n() == target_pts.get_max())
             target_pts.resize(std::max(10.0, 1.5*target_pts.get_max()));

      target_pts.vi_wr[2*target_pts.get_n()] = j;
      target_pts.vi_wr[2*target_pts.get_n() + 1] = i / 3;

      target_pts.vr_wr[3*target_pts.get_n()] = xyz[i];
      target_pts.vr_wr[3*target_pts.get_n() + 1] = xyz[i + 1];
      target_pts.vr_wr[3*target_pts.get_n() + 2] = xyz[i + 2];
      target_pts.inc_n();
    }

  } // end for (unsigned int i = 0; ..

  int num_to_me = 0;
  for (unsigned int i = 0; i < target_pts.get_n(); i++)
    if (target_pts.vi_rd[2*i] == (int)my_rank)
      num_to_me++;
#ifdef VERBOSE
  printf("rank: %u local points: %u, nb sent target pts: %u mappedPts: %u num to me: %d \n",
         my_rank, num_points, target_pts.get_n(), mappedPts->get_n(), num_to_me);
#endif
  // Perform scatter/gather, to gather points to source mesh procs
  if (myPc) {
    (myPc->proc_config().crystal_router())->gs_transfer(1, target_pts, 0);

    num_to_me = 0;
    for (unsigned int i = 0; i < target_pts.get_n(); i++) {
      if (target_pts.vi_rd[2*i] == (int)my_rank)
        num_to_me++;
    }
#ifdef VERBOSE
    printf("rank: %u after first gs nb received_pts: %u; num_from_me = %d\n",
           my_rank, target_pts.get_n(), num_to_me);
#endif
    // After scatter/gather:
    // target_pts.set_n(# points local proc has to map);
    // target_pts.vi_wr[2*i] = proc sending point i
    // target_pts.vi_wr[2*i + 1] = index of point i on sending proc
    // target_pts.vr_wr[3*i..3*i + 2] = xyz of point i
    //
    // Mapping builds the tuple list:
    // source_pts.set_n(target_pts.get_n())
    // source_pts.vi_wr[3*i] = target_pts.vi_wr[2*i] = sending proc
    // source_pts.vi_wr[3*i + 1] = index of point i on sending proc
    // source_pts.vi_wr[3*i + 2] = index of mapped point (-1 if not mapped)
    //
    // Also, mapping builds local tuple_list mappedPts:
    // mappedPts->set_n( # mapped points );
    // mappedPts->vul_wr[i] = local handle of mapped entity
    // mappedPts->vr_wr[3*i..3*i + 2] = natural coordinates in mapped entity

    // Test target points against my elements
    for (unsigned i = 0; i < target_pts.get_n(); i++) {
      result = test_local_box(target_pts.vr_wr + 3*i,
                              target_pts.vi_rd[2*i], target_pts.vi_rd[2*i + 1], i,
                              point_located, rel_eps, abs_eps, &source_pts);
      if (MB_SUCCESS != result)
        return result;
    }

    // No longer need target_pts
    target_pts.reset();
#ifdef VERBOSE
    printf("rank: %u nb sent source pts: %u, mappedPts now: %u\n",
           my_rank, source_pts.get_n(),  mappedPts->get_n());
#endif
      // Send target points back to target procs
    (myPc->proc_config().crystal_router())->gs_transfer(1, source_pts, 0);

#ifdef VERBOSE
    printf("rank: %u nb received source pts: %u\n",
           my_rank, source_pts.get_n());
#endif
  }

  // Store proc/index tuples in targetPts, and/or pass back to application;
  // the tuple this gets stored to looks like:
  // tl.set_n(# mapped points);
  // tl.vi_wr[3*i] = remote proc mapping point
  // tl.vi_wr[3*i + 1] = local index of mapped point
  // tl.vi_wr[3*i + 2] = remote index of mapped point
  //
  // Local index is mapped into either myRange, holding the handles of
  // local mapped entities, or myXyz, holding locations of mapped pts

  // Store information about located points
  TupleList *tl_tmp;
  if (!store_local)
    tl_tmp = tl;
  else {
    targetPts = new TupleList();
    tl_tmp = targetPts;
  }

  tl_tmp->initialize(3, 0, 0, 0, num_points);
  tl_tmp->set_n(num_points); // Automatically sets tl to write_enabled
  // Initialize so we know afterwards how many pts weren't located
  std::fill(tl_tmp->vi_wr, tl_tmp->vi_wr + 3*num_points, -1);

  unsigned int local_pts = 0;
  for (unsigned int i = 0; i < source_pts.get_n(); i++) {
    if (-1 != source_pts.vi_rd[3*i + 2]) { // Why bother sending message saying "i don't have the point" if it gets discarded?
      int tgt_index = 3*source_pts.vi_rd[3*i + 1];
      // Prefer always entities that are local, from the source_pts
      // if a local entity was already found to contain the target point, skip
      // tl_tmp->vi_wr[tgt_index] is -1 initially, but it could already be set with
      // a remote processor
      if (tl_tmp->vi_wr[tgt_index] != (int)my_rank) {
        tl_tmp->vi_wr[tgt_index]     = source_pts.vi_rd[3*i];
        tl_tmp->vi_wr[tgt_index + 1] = source_pts.vi_rd[3*i + 1];
        tl_tmp->vi_wr[tgt_index + 2] = source_pts.vi_rd[3*i + 2];
      }
    }
  }

  // Count missing points
  unsigned int missing_pts = 0;
  for (unsigned int i = 0; i < num_points; i++) {
    if (tl_tmp->vi_rd[3*i + 1] == -1) {
      missing_pts++;
#ifdef VERBOSE
      printf("missing point at index i:  %d -> %15.10f %15.10f %15.10f\n", i, xyz[3*i], xyz[3*i + 1], xyz[3*i + 2]);
#endif
    }
    else if (tl_tmp->vi_rd[3*i] == (int)my_rank)
      local_pts++;
  }
#ifdef VERBOSE
  printf("rank: %u point location: wanted %u got %u locally, %u remote, missing %u\n",
         my_rank, num_points, local_pts, num_points - missing_pts - local_pts, missing_pts);
#endif
  assert(0 == missing_pts); // Will likely break on curved geometries

  // No longer need source_pts
  source_pts.reset();

  // Copy into tl if passed in and storing locally
  if (tl && store_local) {
    tl->initialize(3, 0, 0, 0, num_points);
    tl->enableWriteAccess();
    memcpy(tl->vi_wr, tl_tmp->vi_rd, 3 * tl_tmp->get_n() * sizeof(int));
    tl->set_n(tl_tmp->get_n());
    tl->disableWriteAccess();
  }

  tl_tmp->disableWriteAccess();

  // Done
  return MB_SUCCESS;
}

ErrorCode Coupler::test_local_box(double *xyz,
                                  int from_proc, int remote_index, int /*index*/,
                                  bool &point_located,
                                  double rel_eps, double abs_eps,
                                  TupleList *tl)
{
  std::vector<EntityHandle> entities;
  std::vector<CartVect> nat_coords;
  bool canWrite = false;
  if (tl) {
    canWrite = tl->get_writeEnabled();
    if (!canWrite) {
      tl->enableWriteAccess();
      canWrite = true;
    }
  }

  if (rel_eps && !abs_eps) {
    // Relative epsilon given, translate to absolute epsilon using box dimensions
    BoundBox box;
    myTree->get_bounding_box(box, &localRoot);
    abs_eps = rel_eps * box.diagonal_length();
  }

  ErrorCode result = nat_param(xyz, entities, nat_coords, abs_eps);
  if (MB_SUCCESS != result)
    return result;

  // If we didn't find any ents and we're looking locally, nothing more to do
  if (entities.empty()) {
    if (tl->get_n() == tl->get_max())
      tl->resize(std::max(10.0, 1.5 * tl->get_max()));

    tl->vi_wr[3 * tl->get_n()] = from_proc;
    tl->vi_wr[3 * tl->get_n() + 1] = remote_index;
    tl->vi_wr[3 * tl->get_n() + 2] = -1;
    tl->inc_n();

    point_located = false;
    return MB_SUCCESS;
  }

  // Grow if we know we'll exceed size
  if (mappedPts->get_n() + entities.size() >= mappedPts->get_max())
    mappedPts->resize(std::max(10.0, 1.5 * mappedPts->get_max()));

  std::vector<EntityHandle>::iterator eit = entities.begin();
  std::vector<CartVect>::iterator ncit = nat_coords.begin();

  mappedPts->enableWriteAccess();
  for (; eit != entities.end(); ++eit, ++ncit) {
    // Store in tuple mappedPts
    mappedPts->vr_wr[3*mappedPts->get_n()] = (*ncit)[0];
    mappedPts->vr_wr[3*mappedPts->get_n() + 1] = (*ncit)[1];
    mappedPts->vr_wr[3*mappedPts->get_n() + 2] = (*ncit)[2];
    mappedPts->vul_wr[mappedPts->get_n()] = *eit;
    mappedPts->inc_n();

    // Also store local point, mapped point indices
    if (tl->get_n() == tl->get_max())
      tl->resize(std::max(10.0, 1.5*tl->get_max()));

    // Store in tuple source_pts
    tl->vi_wr[3*tl->get_n()] = from_proc;
    tl->vi_wr[3*tl->get_n() + 1] = remote_index;
    tl->vi_wr[3*tl->get_n() + 2] = mappedPts->get_n()-1;
    tl->inc_n();
  }

  point_located = true;

  if (tl && !canWrite)
    tl->disableWriteAccess();

  return MB_SUCCESS;
}

ErrorCode Coupler::interpolate( Coupler::Method method,
                                const std::string &interp_tag,
                                double *interp_vals,
                                TupleList *tl,
                                bool normalize)
{
  Tag tag;
  ErrorCode result ;
  if (_spectralSource) {
    result = mbImpl->tag_get_handle(interp_tag.c_str(), _ntot, MB_TYPE_DOUBLE, tag);MB_CHK_SET_ERR(result, "Failed to get handle for interpolation tag \"" << interp_tag << "\"");
  }
  else {
    result = mbImpl->tag_get_handle(interp_tag.c_str(), 1, MB_TYPE_DOUBLE, tag);MB_CHK_SET_ERR(result, "Failed to get handle for interpolation tag \"" << interp_tag << "\"");
  }

  return interpolate(method, tag, interp_vals, tl, normalize);
}

ErrorCode Coupler::interpolate(Coupler::Method *methods,
                               Tag *tags,
                               int *points_per_method,
                               int num_methods,
                               double *interp_vals,
                               TupleList *tl,
                               bool /* normalize */)
{
  //if (!((LINEAR_FE == method) || (CONSTANT == method)))
  // return MB_FAILURE;

  // remote pts first
  TupleList *tl_tmp = (tl ? tl : targetPts);

  ErrorCode result = MB_SUCCESS;

  unsigned int pts_total = 0;
  for (int i = 0; i < num_methods; i++)
    pts_total += points_per_method[i];

  // If tl was passed in non-NULL, just have those points, otherwise have targetPts plus
  // locally mapped pts
  if (pts_total != tl_tmp->get_n())
    return MB_FAILURE;

  TupleList tinterp;
  tinterp.initialize(5, 0, 0, 1, tl_tmp->get_n());
  int t = 0;
  tinterp.enableWriteAccess();
  for (int i = 0; i < num_methods; i++) {
    for (int j = 0; j < points_per_method[i]; j++) {
      tinterp.vi_wr[5*t] = tl_tmp->vi_rd[3*t];
      tinterp.vi_wr[5*t + 1] = tl_tmp->vi_rd[3*t + 1];
      tinterp.vi_wr[5*t + 2] = tl_tmp->vi_rd[3*t + 2];
      tinterp.vi_wr[5*t + 3] = methods[i];
      tinterp.vi_wr[5*t + 4] = i;
      tinterp.vr_wr[t] = 0.0;
      tinterp.inc_n();
      t++;
    }
  }

  // Scatter/gather interpolation points
  if (myPc) {
    (myPc->proc_config().crystal_router())->gs_transfer(1, tinterp, 0);

    // Perform interpolation on local source mesh; put results into
    // tinterp.vr_wr[i]

    mappedPts->enableWriteAccess();
    for (unsigned int i = 0; i < tinterp.get_n(); i++) {
      int mindex = tinterp.vi_rd[5*i + 2];
      Method method = (Method)tinterp.vi_rd[5*i + 3];
      Tag tag = tags[tinterp.vi_rd[5*i + 4]];

      result = MB_FAILURE;
      if (LINEAR_FE == method || QUADRATIC_FE == method || SPHERICAL==method) {
        result = interp_field(mappedPts->vul_rd[mindex],
                              CartVect(mappedPts->vr_wr + 3*mindex),
                              tag, tinterp.vr_wr[i]);
      }
      else if (CONSTANT == method) {
        result = constant_interp(mappedPts->vul_rd[mindex],
                                 tag, tinterp.vr_wr[i]);
      }

      if (MB_SUCCESS != result)
        return result;
    }

    // Scatter/gather interpolation data
    myPc->proc_config().crystal_router()->gs_transfer(1, tinterp, 0);
  }

  // Copy the interpolated field as a unit
  for (unsigned int i = 0; i < tinterp.get_n(); i++)
    interp_vals[tinterp.vi_rd[5*i + 1]] = tinterp.vr_rd[i];

  // Done
  return MB_SUCCESS;
}

ErrorCode Coupler::nat_param(double xyz[3],
                             std::vector<EntityHandle> &entities,
                             std::vector<CartVect> &nat_coords,
                             double epsilon)
{
  if (!myTree)
    return MB_FAILURE;

  AdaptiveKDTreeIter treeiter;
  ErrorCode result = myTree->get_tree_iterator(localRoot, treeiter);
  if (MB_SUCCESS != result) {
    std::cout << "Problems getting iterator" << std::endl;
    return result;
  }

  EntityHandle closest_leaf;
  if (epsilon) {
    std::vector<double> dists;
    std::vector<EntityHandle> leaves;

    // Two tolerances
    result = myTree->distance_search(xyz, epsilon, leaves,
                                     /*iter_tol*/ epsilon,
                                     /*inside_tol*/ 10*epsilon,
                                     &dists, NULL, &localRoot);
    if (leaves.empty())
      // Not found returns success here, with empty list, just like case with no epsilon
      return MB_SUCCESS;
    // Get closest leaf
    double min_dist = *dists.begin();
    closest_leaf = *leaves.begin();
    std::vector<EntityHandle>::iterator vit = leaves.begin() + 1;
    std::vector<double>::iterator dit = dists.begin() + 1;
    for (; vit != leaves.end() && min_dist; ++vit, ++dit) {
      if (*dit < min_dist) {
        min_dist = *dit;
        closest_leaf = *vit;
      }
    }
  }
  else {
    result = myTree->point_search(xyz, treeiter, 1.0e-10, 1.0e-6, NULL, &localRoot);
    if (MB_ENTITY_NOT_FOUND == result) // Point is outside of myTree's bounding box
      return MB_SUCCESS;
    else if (MB_SUCCESS != result) {
      std::cout << "Problems getting leaf \n";
      return result;
    }
    closest_leaf = treeiter.handle();
  }

  // Find natural coordinates of point in element(s) in that leaf
  CartVect tmp_nat_coords;
  Range range_leaf;
  result = mbImpl->get_entities_by_dimension(closest_leaf, max_dim, range_leaf, false);
  if (MB_SUCCESS != result)
    std::cout << "Problem getting leaf in a range" << std::endl;

  // Loop over the range_leaf
  for (Range::iterator iter = range_leaf.begin(); iter != range_leaf.end(); ++iter) {
    // Test to find out in which entity the point is
    // Get the EntityType and create the appropriate Element::Map subtype
    // If spectral, do not need coordinates, just the GL points
    EntityType etype = mbImpl->type_from_handle(*iter);
    if (NULL != this->_spectralSource && MBHEX == etype) {
      EntityHandle eh = *iter;
      const double * xval;
      const double * yval;
      const double * zval;
      ErrorCode rval = mbImpl->tag_get_by_ptr(_xm1Tag, &eh, 1, (const void**)&xval);
      if (moab::MB_SUCCESS != rval) {
        std::cout << "Can't get xm1 values \n";
        return MB_FAILURE;
      }
      rval = mbImpl->tag_get_by_ptr(_ym1Tag, &eh, 1, (const void**)&yval);
      if (moab::MB_SUCCESS != rval) {
        std::cout << "Can't get ym1 values \n";
        return MB_FAILURE;
      }
      rval = mbImpl->tag_get_by_ptr(_zm1Tag, &eh, 1, (const void**)&zval);
      if (moab::MB_SUCCESS != rval) {
        std::cout << "Can't get zm1 values \n";
        return MB_FAILURE;
      }
      Element::SpectralHex* spcHex = (Element::SpectralHex*)_spectralSource;

      spcHex->set_gl_points((double*)xval, (double*)yval, (double*)zval);
      try {
        tmp_nat_coords = spcHex->ievaluate(CartVect(xyz), epsilon ); // introduce
        bool inside = spcHex->inside_nat_space(CartVect(tmp_nat_coords), epsilon);
        if (!inside) {
#ifdef VERBOSE
          std::cout << "point " << xyz[0] << " " << xyz[1] << " " << xyz[2] <<
              " is not converging inside hex " << mbImpl->id_from_handle(eh) << "\n";
#endif
          continue; // It is possible that the point is outside, so it will not converge
        }
      }
      catch (Element::Map::EvaluationError&) {
        continue;
      }


    }
    else {
      const EntityHandle *connect;
      int num_connect;

      // Get connectivity
      result = mbImpl->get_connectivity(*iter, connect, num_connect, true);
      if (MB_SUCCESS != result)
        return result;

      // Get coordinates of the vertices
      std::vector<CartVect> coords_vert(num_connect);
      result = mbImpl->get_coords(connect, num_connect, &(coords_vert[0][0]));
      if (MB_SUCCESS != result) {
        std::cout << "Problems getting coordinates of vertices\n";
        return result;
      }
      CartVect pos(xyz);
      if (MBHEX == etype) {
        if (8 == num_connect) {
          Element::LinearHex hexmap(coords_vert);
          if (!hexmap.inside_box(pos, epsilon))
            continue;
          try {
            tmp_nat_coords = hexmap.ievaluate(pos, epsilon);
            bool inside = hexmap.inside_nat_space(tmp_nat_coords, epsilon);
            if (!inside)
              continue;
          }
          catch (Element::Map::EvaluationError&) {
            continue;
          }
        }
        else if (27 == num_connect) {
          Element::QuadraticHex hexmap(coords_vert);
         if (!hexmap.inside_box(pos, epsilon))
           continue;
         try {
           tmp_nat_coords = hexmap.ievaluate(pos, epsilon);
           bool inside = hexmap.inside_nat_space(tmp_nat_coords, epsilon);
           if (!inside)
             continue;
         }
         catch (Element::Map::EvaluationError&) {
           continue;
         }
        }
        else // TODO this case not treated yet, no interpolation
          continue;
      }
      else if (MBTET == etype) {
        Element::LinearTet tetmap(coords_vert);
        // This is just a linear solve; unless degenerate, will not except
        tmp_nat_coords = tetmap.ievaluate(pos);
        bool inside = tetmap.inside_nat_space(tmp_nat_coords, epsilon);
        if (!inside)
          continue;
      }
      else if (MBQUAD == etype && spherical) {
        Element::SphericalQuad sphermap(coords_vert);
        /* skip box test, because it can filter out good elements with high curvature
         * if (!sphermap.inside_box(pos, epsilon))
          continue;*/
        try {
          tmp_nat_coords = sphermap.ievaluate(pos, epsilon);
          bool inside = sphermap.inside_nat_space(tmp_nat_coords, epsilon);
          if (!inside)
            continue;
        }
        catch (Element::Map::EvaluationError&) {
          continue;
        }

      }
      else if (MBTRI == etype && spherical) {
        Element::SphericalTri sphermap(coords_vert);
        /* skip box test, because it can filter out good elements with high curvature
         * if (!sphermap.inside_box(pos, epsilon))
            continue;*/
        try {
          tmp_nat_coords = sphermap.ievaluate(pos, epsilon);
          bool inside = sphermap.inside_nat_space(tmp_nat_coords, epsilon);
          if (!inside)
            continue;
        }
        catch (Element::Map::EvaluationError&) {
          continue;
        }
      }

      else if (MBQUAD == etype) {
        Element::LinearQuad quadmap(coords_vert);
        if (!quadmap.inside_box(pos, epsilon))
          continue;
        try {
          tmp_nat_coords = quadmap.ievaluate(pos, epsilon);
          bool inside = quadmap.inside_nat_space(tmp_nat_coords, epsilon);
          if (!inside)
            continue;
        }
        catch (Element::Map::EvaluationError&) {
          continue;
        }
        if (!quadmap.inside_nat_space(tmp_nat_coords, epsilon))
          continue;
      }
      /*
      else if (etype == MBTRI){
        Element::LinearTri trimap(coords_vert);
        if (!trimap.inside_box( pos, epsilon))
          continue;
        try {
          tmp_nat_coords = trimap.ievaluate(pos, epsilon);
          bool inside = trimap.inside_nat_space(tmp_nat_coords, epsilon);
          if (!inside) continue;
        }
        catch (Element::Map::EvaluationError) {
          continue;
        }
        if (!trimap.inside_nat_space(tmp_nat_coords, epsilon))
          continue;
      }
      */
      else if (etype == MBEDGE){
        Element::LinearEdge edgemap(coords_vert);
        try {
          tmp_nat_coords = edgemap.ievaluate(CartVect(xyz), epsilon);
        }
        catch (Element::Map::EvaluationError) {
          continue;
        }
        if (!edgemap.inside_nat_space(tmp_nat_coords, epsilon))
          continue;
      }
      else {
        std::cout << "Entity not Hex/Tet/Quad/Tri/Edge. Please verify." << std::endl;
        continue;
      }
    }
    // If we get here then we've found the coordinates.
    // Save them and the entity and return success.
    entities.push_back(*iter);
    nat_coords.push_back(tmp_nat_coords);
    return MB_SUCCESS;
  }

  // Didn't find any elements containing the point
  return MB_SUCCESS;
}

ErrorCode Coupler::interp_field(EntityHandle elem,
                                CartVect nat_coord,
                                Tag tag,
                                double &field)
{
  if (_spectralSource) {
    // Get tag values at the GL points for some field (Tag)
    const double * vx;
    ErrorCode rval = mbImpl-> tag_get_by_ptr(tag, &elem, 1, (const void**)&vx);
    if (moab::MB_SUCCESS != rval) {
      std::cout << "Can't get field values for the tag \n";
      return MB_FAILURE;
    }
    Element::SpectralHex * spcHex = (Element::SpectralHex*)_spectralSource;
    field = spcHex->evaluate_scalar_field(nat_coord, vx);
  }
  else {
    double vfields[27]; // Will work for linear hex, quadratic hex or Tets
    moab::Element::Map *elemMap = NULL;
    int num_verts = 0;
    // Get the EntityType
    // Get the tag values at the vertices
    const EntityHandle *connect;
    int num_connect;
    ErrorCode result = mbImpl->get_connectivity(elem, connect, num_connect);
    if (MB_SUCCESS != result)
      return result;
    EntityType etype = mbImpl->type_from_handle(elem);
    if (MBHEX == etype) {
      if (8 == num_connect) {
        elemMap = new moab::Element::LinearHex();
        num_verts = 8;
      }
      else { /* (MBHEX == etype && 27 == num_connect) */
        elemMap = new moab::Element::QuadraticHex();
        num_verts = 27;
      }
    }
    else if (MBTET == etype) {
      elemMap = new moab::Element::LinearTet();
      num_verts = 4;
    }
    else if (MBQUAD == etype) {
      elemMap = new moab::Element::LinearQuad();
      num_verts = 4;
    }
    else if (MBTRI == etype) {
      elemMap = new moab::Element::LinearTri();
      num_verts = 3;
    }
    else
      return MB_FAILURE;

    result = mbImpl->tag_get_data(tag, connect, std::min(num_verts, num_connect), vfields);
    if (MB_SUCCESS != result) {
      delete elemMap;
      return result;
    }

    // Function for the interpolation
    field = 0;

    // Check the number of vertices
    assert(num_connect >= num_verts);

    // Calculate the field
    try {
      field = elemMap->evaluate_scalar_field(nat_coord, vfields);
    }
    catch (moab::Element::Map::EvaluationError&) {
      delete elemMap;
      return MB_FAILURE;
    }

    delete elemMap;
  }

  return MB_SUCCESS;
}

// Simplest "interpolation" for element-based source fields. Set the value of the field
// at the target point to that of the field in the source element it lies in.
ErrorCode Coupler::constant_interp(EntityHandle elem,
                                   Tag tag,
                                   double &field)
{
  double tempField;

  // Get the tag values at the vertices
  ErrorCode result = mbImpl->tag_get_data(tag, &elem, 1, &tempField);
  if (MB_SUCCESS != result)
    return result;

  field = tempField;

  return MB_SUCCESS;
}

// Normalize a field over the entire mesh represented by the root_set.
ErrorCode Coupler::normalize_mesh(EntityHandle          root_set,
                                  const char            *norm_tag,
                                  Coupler::IntegType    integ_type,
                                  int                   num_integ_pts)
{
  ErrorCode err;

  // SLAVE START ****************************************************************
  // Search for entities based on tag_handles and tag_values
  std::vector< std::vector<EntityHandle> > entity_sets;
  std::vector< std::vector<EntityHandle> > entity_groups;

  // put the root_set into entity_sets
  std::vector<EntityHandle> ent_set;
  ent_set.push_back(root_set);
  entity_sets.push_back(ent_set);

  // get all entities from root_set and put into entity_groups
  std::vector<EntityHandle> entities;
  err = mbImpl->get_entities_by_handle(root_set, entities, true);
  ERRORR("Failed to get entities in root_set.", err);

  entity_groups.push_back(entities);

  // Call do_normalization() to continue common normalization processing
  err = do_normalization(norm_tag, entity_sets, entity_groups, integ_type, num_integ_pts);
  ERRORR("Failure in do_normalization().", err);
  // SLAVE END   ****************************************************************

  return err;
}

// Normalize a field over the subset of entities identified by the tags and values passed
ErrorCode Coupler::normalize_subset(EntityHandle        root_set,
                                    const char          *norm_tag,
                                    const char          **tag_names,
                                    int                 num_tags,
                                    const char          **tag_values,
                                    Coupler::IntegType  integ_type,
                                    int                 num_integ_pts)
{
  moab::ErrorCode err;
  std::vector<Tag> tag_handles;

  // Lookup tag handles from tag names
  for (int t = 0; t < num_tags; t++) {
    // get tag handle & size
    Tag th;
    err = mbImpl->tag_get_handle(tag_names[t], 1, moab::MB_TYPE_DOUBLE, th, moab::MB_TAG_ANY);
    ERRORR("Failed to get tag handle.", err);
    tag_handles.push_back(th);
  }

  return normalize_subset(root_set,
                          norm_tag,
                          &tag_handles[0],
                          num_tags,
                          tag_values,
                          integ_type,
                          num_integ_pts);
}

ErrorCode Coupler::normalize_subset(EntityHandle       root_set,
                                    const char         *norm_tag,
                                    Tag                *tag_handles,
                                    int                num_tags,
                                    const char         **tag_values,
                                    Coupler::IntegType integ_type,
                                    int                num_integ_pts)
{
  ErrorCode err;

  // SLAVE START ****************************************************************
  // Search for entities based on tag_handles and tag_values
  std::vector< std::vector<EntityHandle> > entity_sets;
  std::vector< std::vector<EntityHandle> > entity_groups;

  err = get_matching_entities(root_set, tag_handles, tag_values, num_tags,
                              &entity_sets, &entity_groups);
  ERRORR("Failed to get matching entities.", err);

  // Call do_normalization() to continue common normalization processing
  err = do_normalization(norm_tag, entity_sets, entity_groups, integ_type, num_integ_pts);
  ERRORR("Failure in do_normalization().", err);
  // SLAVE END   ****************************************************************

  return err;
}

ErrorCode Coupler::do_normalization(const char                               *norm_tag,
                                    std::vector<std::vector<EntityHandle> >  &entity_sets,
                                    std::vector<std::vector<EntityHandle> >  &entity_groups,
                                    Coupler::IntegType                       integ_type,
                                    int                                      num_integ_pts)
{
  // SLAVE START ****************************************************************
  ErrorCode err;
  int ierr = 0;

  // Setup data for parallel computing
  int nprocs, rank;
  ierr = MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
  ERRORMPI("Getting number of procs failed.", ierr);
  ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  ERRORMPI("Getting rank failed.", ierr);

  // Get the integrated field value for each group(vector) of entities.
  // If no entities are in a group then a zero will be put in the list
  // of return values.
  unsigned int num_ent_grps = entity_groups.size();
  std::vector<double> integ_vals(num_ent_grps);

  err = get_group_integ_vals(entity_groups, integ_vals, norm_tag, num_integ_pts, integ_type);
  ERRORR("Failed to get integrated field values for groups in mesh.", err);
  // SLAVE END   ****************************************************************

  // SLAVE/MASTER START #########################################################
  // Send list of integrated values back to master proc. The ordering of the
  // values will match the ordering of the entity groups (i.e. vector of vectors)
  // sent from master to slaves earlier.  The values for each entity group will
  // be summed during the transfer.
  std::vector<double> sum_integ_vals(num_ent_grps);

  if (nprocs > 1) {
    // If parallel then send the values back to the master.
    ierr = MPI_Reduce(&integ_vals[0], &sum_integ_vals[0], num_ent_grps, MPI_DOUBLE,
                     MPI_SUM, MASTER_PROC, myPc->proc_config().proc_comm());
    ERRORMPI("Transfer and reduction of integrated values failed.", ierr);
  }
  else {
    // Otherwise just copy the vector
    sum_integ_vals = integ_vals;
  }
  // SLAVE/MASTER END   #########################################################

  // MASTER START ***************************************************************
  // Calculate the normalization factor for each group by taking the
  // inverse of each integrated field value. Put the normalization factor
  // for each group back into the list in the same order.
  for (unsigned int i = 0; i < num_ent_grps; i++) {
    double val = sum_integ_vals[i];
    if (fabs(val) > 1e-8) sum_integ_vals[i] = 1.0/val;
    else {
      sum_integ_vals[i] = 0.0; /* VSM: not sure what we should do here ? */
      /* commenting out error below since if integral(value)=0.0, then normalization
         is probably unnecessary to start with ? */
      /* ERRORR("Integrating an invalid field -- integral("<<norm_tag<<") = "<<val<<".", err); */
    }
  }
  // MASTER END   ***************************************************************

  // MASTER/SLAVE START #########################################################
  if (nprocs > 1) {
    // If parallel then broadcast the normalization factors to the procs.
    ierr = MPI_Bcast(&sum_integ_vals[0], num_ent_grps, MPI_DOUBLE, MASTER_PROC,
                    myPc->proc_config().proc_comm());
    ERRORMPI("Broadcast of normalization factors failed.", ierr);
  }
  // MASTER/SLAVE END   #########################################################

  // SLAVE START ****************************************************************
  // Save the normalization factors to a new tag with name of norm_tag's value
  // and the string "_normF" appended.  This new tag will be created on the entity
  // set that contains all of the entities from a group.

  err = apply_group_norm_factor(entity_sets, sum_integ_vals, norm_tag, integ_type);
  ERRORR("Failed to set the normalization factor for groups in mesh.", err);
  // SLAVE END   ****************************************************************

  return err;
}

// Functions supporting the subset normalization function

// Retrieve groups of entities matching tags and values if present
ErrorCode Coupler::get_matching_entities(EntityHandle                            root_set,
                                         const char                              **tag_names,
                                         const char                              **tag_values,
                                         int                                     num_tags,
                                         std::vector<std::vector<EntityHandle> > *entity_sets,
                                         std::vector<std::vector<EntityHandle> > *entity_groups)
{
  ErrorCode err;
  std::vector<Tag> tag_handles;

  for (int t = 0; t < num_tags; t++) {
    // Get tag handle & size
    Tag th;
    err = mbImpl->tag_get_handle(tag_names[t], 1, moab::MB_TYPE_DOUBLE, th, moab::MB_TAG_ANY);
    ERRORR("Failed to get tag handle.", err);
    tag_handles.push_back(th);
  }

  return get_matching_entities(root_set, &tag_handles[0], tag_values, num_tags,
                               entity_sets, entity_groups);
}

// Retrieve groups of entities matching tags and values if present
ErrorCode Coupler::get_matching_entities(EntityHandle                            root_set,
                                         Tag                                     *tag_handles,
                                         const char                              **tag_values,
                                         int                                     num_tags,
                                         std::vector<std::vector<EntityHandle> > *entity_sets,
                                         std::vector<std::vector<EntityHandle> > *entity_groups)
{
  // SLAVE START ****************************************************************

  // Setup data for parallel computing
  ErrorCode err;
  int ierr = 0;
  int nprocs, rank;
  ierr = MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
  ERRORMPI("Getting number of procs failed.", ierr);
  ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  ERRORMPI("Getting rank failed.", ierr);

  Range ent_sets;
  err = mbImpl->get_entities_by_type_and_tag(root_set, moab::MBENTITYSET, tag_handles,
                                             (const void* const *)tag_values,
                                             num_tags, ent_sets, Interface::INTERSECT, 0);
  ERRORR("Core::get_entities_by_type_and_tag failed.", err);

  TupleList *tag_list = NULL;
  err = create_tuples(ent_sets, tag_handles, num_tags, &tag_list);
  ERRORR("Failed to create tuples from entity sets.", err);

  // Free up range
  ent_sets.clear();
  // SLAVE END   ****************************************************************

  // If we are running in a multi-proc session then send tuple list back to master
  // proc for consolidation. Otherwise just copy the pointer to the tuple_list.
  TupleList *cons_tuples;
  if (nprocs > 1) {
    // SLAVE/MASTER START #########################################################

    // Pack the tuple_list in a buffer.
    uint *tuple_buf;
    int tuple_buf_len;
    tuple_buf_len = pack_tuples(tag_list, (void**)&tuple_buf);

    // Free tag_list here as its not used again if nprocs > 1
    tag_list->reset();

    // Send back the buffer sizes to the master proc
    int *recv_cnts = (int*)malloc(nprocs * sizeof(int));
    int *offsets   = (int*)malloc(nprocs * sizeof(int));
    uint *all_tuples_buf = NULL;

    MPI_Gather(&tuple_buf_len, 1, MPI_INT, recv_cnts, 1, MPI_INT, MASTER_PROC,
               myPc->proc_config().proc_comm());
    ERRORMPI("Gathering buffer sizes failed.", err);

    // Allocate a buffer large enough for all the data
    if (rank == MASTER_PROC) {
      int all_tuples_len = recv_cnts[0];
      offsets[0] = 0;
      for (int i = 1; i < nprocs; i++) {
        offsets[i] = offsets[i - 1] + recv_cnts[i - 1];
        all_tuples_len += recv_cnts[i];
      }

      all_tuples_buf = (uint*)malloc(all_tuples_len * sizeof(uint));
    }

    // Send all buffers to the master proc for consolidation
    MPI_Gatherv((void*)tuple_buf, tuple_buf_len, MPI_INT,
                (void*)all_tuples_buf, recv_cnts, offsets, MPI_INT, MASTER_PROC,
                myPc->proc_config().proc_comm());
    ERRORMPI("Gathering tuple_lists failed.", err);
    free(tuple_buf); // malloc'd in pack_tuples

    if (rank == MASTER_PROC) {
      // Unpack the tuple_list from the buffer.
      TupleList **tl_array = (TupleList**)malloc(nprocs * sizeof(TupleList*));
      for (int i = 0; i < nprocs; i++)
        unpack_tuples((void*) &all_tuples_buf[offsets[i]], &tl_array[i]);

      // Free all_tuples_buf here as it is only allocated on the MASTER_PROC
      free(all_tuples_buf);
      // SLAVE/MASTER END   #########################################################

      // MASTER START ***************************************************************
      // Consolidate all tuple_lists into one tuple_list with no duplicates.
      err = consolidate_tuples(tl_array, nprocs, &cons_tuples);
      ERRORR("Failed to consolidate tuples.", err);

      for (int i = 0; i < nprocs; i++)
        tl_array[i]->reset();
      free(tl_array);
      // MASTER END   ***************************************************************
    }

    // Free offsets and recv_cnts as they are allocated on all procs
    free(offsets);
    free(recv_cnts);

    // MASTER/SLAVE START #########################################################
    // Broadcast condensed tuple list back to all procs.
    uint *ctl_buf;
    int ctl_buf_sz;
    if (rank == MASTER_PROC)
      ctl_buf_sz = pack_tuples(cons_tuples, (void**)&ctl_buf);

    // Send buffer size
    ierr = MPI_Bcast(&ctl_buf_sz, 1, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm());
    ERRORMPI("Broadcasting tuple_list size failed.", ierr);

    // Allocate a buffer in the other procs
    if (rank != MASTER_PROC)
      ctl_buf = (uint*)malloc(ctl_buf_sz * sizeof(uint));

    ierr = MPI_Bcast((void*)ctl_buf, ctl_buf_sz, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm());
    ERRORMPI("Broadcasting tuple_list failed.", ierr);

    if (rank != MASTER_PROC)
      unpack_tuples(ctl_buf, &cons_tuples);
    free(ctl_buf);
    // MASTER/SLAVE END   #########################################################
  }
  else
    cons_tuples = tag_list;

  // SLAVE START ****************************************************************
  // Loop over the tuple list getting the entities with the tags in the tuple_list entry
  uint mi, ml, mul, mr;
  cons_tuples->getTupleSize(mi, ml, mul, mr);

  for (unsigned int i = 0; i < cons_tuples->get_n(); i++) {
    // Get Entity Sets that match the tags and values.

    // Convert the data in the tuple_list to an array of pointers to the data
    // in the tuple_list as that is what the iMesh API call is expecting.
    int **vals = (int**)malloc(mi * sizeof(int*));
    for (unsigned int j = 0; j < mi; j++)
      vals[j] = (int *)&(cons_tuples->vi_rd[(i*mi) + j]);

    // Get entities recursively based on type and tag data
    err = mbImpl->get_entities_by_type_and_tag(root_set, moab::MBENTITYSET, tag_handles,
                                               (const void* const *)vals,
                                               mi, ent_sets, Interface::INTERSECT, 0);
    ERRORR("Core::get_entities_by_type_and_tag failed.", err);
    if (debug) std::cout << "ent_sets_size=" << ent_sets.size() << std::endl;

    // Free up the array of pointers
    free(vals);

    // Loop over the entity sets and then free the memory for ent_sets.
    std::vector<EntityHandle> ent_set_hdls;
    std::vector<EntityHandle> ent_hdls;
    for (unsigned int j = 0; j < ent_sets.size(); j++) {
      // Save the entity set
      ent_set_hdls.push_back(ent_sets[j]);

      // Get all entities for the entity set
      Range ents;

      /* VSM: do we need to filter out entity sets ? */
      err = mbImpl->get_entities_by_handle(ent_sets[j], ents, false);
      ERRORR("Core::get_entities_by_handle failed.", err);
      if (debug) std::cout << "ents_size=" << ents.size() << std::endl;

      // Save all of the entities from the entity set and free the memory for ents.
      for (unsigned int k = 0; k < ents.size(); k++) {
        ent_hdls.push_back(ents[k]);
      }
      ents.clear();
      if (debug) std::cout << "ent_hdls.size=" << ent_hdls.size() << std::endl;
    }

    // Free the entity set list for next tuple iteration.
    ent_sets.clear();

    // Push ent_set_hdls onto entity_sets, ent_hdls onto entity_groups
    // and clear both ent_set_hdls and ent_hdls.
    entity_sets->push_back(ent_set_hdls);
    ent_set_hdls.clear();
    entity_groups->push_back(ent_hdls);
    ent_hdls.clear();
    if (debug) std::cout << "entity_sets->size=" << entity_sets->size()
                         << ", entity_groups->size=" << entity_groups->size() << std::endl;
  }

  cons_tuples->reset();
  // SLAVE END   ****************************************************************

  return err;
}

// Return a tuple_list containing  tag values for each Entity Set
// The tuple_list will have a column for each tag and a row for each
// Entity Set. It is assumed all of the tags are integer tags.
ErrorCode Coupler::create_tuples(Range        &ent_sets,
                                 const char   **tag_names,
                                 unsigned int num_tags,
                                 TupleList    **tuple_list)
{
  ErrorCode err;
  std::vector<Tag> tag_handles;

  for (unsigned int t = 0; t < num_tags; t++) {
    // Get tag handle & size
    Tag th;
    err = mbImpl->tag_get_handle(tag_names[t], 1, moab::MB_TYPE_DOUBLE, th, moab::MB_TAG_ANY);
    ERRORR("Failed to get tag handle.", err);
    tag_handles.push_back(th);
  }

  return create_tuples(ent_sets, &tag_handles[0], num_tags, tuple_list);
}

// Return a tuple_list containing  tag values for each Entity Set
// The tuple_list will have a column for each tag and a row for each
// Entity Set.  It is assumed all of the tags are integer tags.
ErrorCode Coupler::create_tuples(Range        &ent_sets,
                                 Tag          *tag_handles,
                                 unsigned int num_tags,
                                 TupleList    **tuples)
{
  // ASSUMPTION: All tags are of type integer.  This may need to be expanded in future.
  ErrorCode err;

  // Allocate a tuple_list for the number of entity sets passed in
  TupleList *tag_tuples = new TupleList(num_tags, 0, 0, 0, (int)ent_sets.size());
  //tag_tuples->initialize(num_tags, 0, 0, 0, num_sets);
  uint mi, ml, mul, mr;
  tag_tuples->getTupleSize(mi, ml, mul, mr);
  tag_tuples->enableWriteAccess();

  if (mi == 0)
    ERRORR("Failed to initialize tuple_list.", MB_FAILURE);

  // Loop over the filtered entity sets retrieving each matching tag value one by one.
  int val;
  for (unsigned int i = 0; i < ent_sets.size(); i++) {
    for (unsigned int j = 0; j < num_tags; j++) {
      EntityHandle set_handle = ent_sets[i];
      err = mbImpl->tag_get_data(tag_handles[j], &set_handle, 1, &val);
      ERRORR("Failed to get integer tag data.", err);
      tag_tuples->vi_wr[i*mi + j] = val;
    }

    // If we get here there was no error so increment n in the tuple_list
    tag_tuples->inc_n();
  }
  tag_tuples->disableWriteAccess();
  *tuples = tag_tuples;

  return MB_SUCCESS;
}

// Consolidate tuple_lists into one list with no duplicates
ErrorCode Coupler::consolidate_tuples(TupleList     **all_tuples,
                                      unsigned int  num_tuples,
                                      TupleList     **unique_tuples)
{
  int total_rcv_tuples = 0;
  int offset = 0, copysz = 0;
  unsigned num_tags = 0;

  uint ml, mul, mr;
  uint *mi = (uint *)malloc(sizeof(uint) * num_tuples);

  for(unsigned int i = 0; i < num_tuples; i++) {
    all_tuples[i]->getTupleSize(mi[i], ml, mul, mr);
  }

  for (unsigned int i = 0; i < num_tuples; i++) {
    if (all_tuples[i] != NULL) {
      total_rcv_tuples += all_tuples[i]->get_n();
      num_tags = mi[i];
    }
  }
  const unsigned int_size = sizeof(sint);
  const unsigned int_width = num_tags * int_size;

  // Get the total size of all of the tuple_lists in all_tuples.
  for (unsigned int i = 0; i < num_tuples; i++) {
    if (all_tuples[i] != NULL)
      total_rcv_tuples += all_tuples[i]->get_n();
  }

  // Copy the tuple_lists into a single tuple_list.
  TupleList *all_tuples_list = new TupleList(num_tags, 0, 0, 0, total_rcv_tuples);
  all_tuples_list->enableWriteAccess();
  //all_tuples_list->initialize(num_tags, 0, 0, 0, total_rcv_tuples);
  for (unsigned int i = 0; i < num_tuples; i++) {
    if (all_tuples[i] != NULL) {
      copysz = all_tuples[i]->get_n() * int_width;
      memcpy(all_tuples_list->vi_wr+offset, all_tuples[i]->vi_rd, copysz);
      offset = offset + (all_tuples[i]->get_n() * mi[i]);
      all_tuples_list->set_n(all_tuples_list->get_n() + all_tuples[i]->get_n());
    }
  }

  // Sort the new tuple_list. Use a radix type sort, starting with the last (or least significant)
  // tag column in the vi array and working towards the first (or most significant) tag column.
  TupleList::buffer sort_buffer;
  sort_buffer.buffer_init(2 * total_rcv_tuples * int_width);
  for (int i = num_tags - 1; i >= 0; i--) {
    all_tuples_list->sort(i, &sort_buffer);
  }

  // Cycle through the sorted list eliminating duplicates.
  // Keep counters to the current end of the tuple_list (w/out dups) and the last tuple examined.
  unsigned int end_idx = 0, last_idx = 1;
  while (last_idx < all_tuples_list->get_n()) {
    if (memcmp(all_tuples_list->vi_rd + (end_idx*num_tags), all_tuples_list->vi_rd + (last_idx*num_tags), int_width) == 0) {
      // Values equal - skip
      last_idx += 1;
    }
    else {
      // Values different - copy
      // Move up the end index
      end_idx += 1;
      memcpy(all_tuples_list->vi_wr + (end_idx*num_tags), all_tuples_list->vi_rd + (last_idx*num_tags), int_width);
      last_idx += 1;
    }
  }
  // Update the count in all_tuples_list
  all_tuples_list->set_n(end_idx + 1);

  // Resize the tuple_list
  all_tuples_list->resize(all_tuples_list->get_n());

  // Set the output parameter
  *unique_tuples = all_tuples_list;

  return MB_SUCCESS;
}

// Calculate integrated field values for groups of entities
ErrorCode Coupler::get_group_integ_vals(std::vector<std::vector<EntityHandle> > &groups,
                                        std::vector<double> &integ_vals,
                                        const char *norm_tag,
                                        int /*num_integ_vals*/,
                                        Coupler::IntegType integ_type)
{
  ErrorCode err;

  std::vector<std::vector<EntityHandle> >::iterator iter_i;
  std::vector<EntityHandle>::iterator iter_j;
  double grp_intrgr_val, intgr_val;

  // Get the tag handle for norm_tag
  Tag norm_hdl;
  err = mbImpl->tag_get_handle(norm_tag, 1, moab::MB_TYPE_DOUBLE, norm_hdl, moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT);
  ERRORR("Failed to get norm_tag handle.", err);

  // Check size of integ_vals vector
  if (integ_vals.size() != groups.size())
    integ_vals.resize(groups.size());

  // Loop over the groups(vectors) of entities
  unsigned int i;
  for (i = 0, iter_i = groups.begin(); iter_i != groups.end(); i++, ++iter_i) {
    grp_intrgr_val = 0;

    // Loop over the all the entities in the group, integrating
    // the field_fn over the entity in iter_j
    for (iter_j = (*iter_i).begin(); iter_j != (*iter_i).end(); ++iter_j) {
      EntityHandle ehandle = (*iter_j);

      // Check that the entity in iter_j is of the same dimension as the
      // integ_type we are performing
      EntityType j_type;
      j_type = mbImpl->type_from_handle(ehandle);
      ERRORR("Failed to get entity type.", err);
      // Skip any entities in the group that are not of the type being considered
      if ((integ_type == VOLUME) && (j_type < MBTET || j_type >= MBENTITYSET))
        continue;

      intgr_val = 0;

      // Retrieve the vertices from the element
      const EntityHandle* verts = NULL;
      int connectivity_size = 0;

      err = mbImpl->get_connectivity(ehandle, verts, connectivity_size, false);
      ERRORR("Failed to get vertices from entity.", err);

      // Get the vertex coordinates and the field values at the vertices.
      double *coords = (double*) malloc(sizeof(double) * (3*connectivity_size));
      /* TODO: VSM: check if this works for lower dimensions also without problems */
      /* if (3 == geom_dim) */
      err = mbImpl->get_coords(verts, connectivity_size, coords);
      ERRORR("Failed to get vertex coordinates.", err);

      /* allocate the field data array */
      double *vfield = (double*) malloc(sizeof(double) * (connectivity_size));
      err = mbImpl->tag_get_data(norm_hdl, verts, connectivity_size, vfield);
      if (MB_SUCCESS != err) {
        free(coords);
      }
      ERRORR("Failed to get vertex coordinates.", err);

      // Get coordinates of all corner vertices (in normal order) and
      // put in array of CartVec.
      std::vector<CartVect> vertices(connectivity_size);

      // Put the vertices into a CartVect vector
      double *x = coords;
      for (int j = 0; j < connectivity_size; j++, x += 3) {
        vertices[j] = CartVect(x);
      }
      free(coords);

      moab::Element::Map *elemMap;
      if (j_type == MBHEX) {
        if (connectivity_size == 8)
          elemMap = new moab::Element::LinearHex(vertices);
        else
          elemMap = new moab::Element::QuadraticHex(vertices);
      }
      else if (j_type == MBTET) {
        elemMap = new moab::Element::LinearTet(vertices);
      }
      else if (j_type == MBQUAD) {
        elemMap = new moab::Element::LinearQuad(vertices);
      }
      /*
      else if (j_type == MBTRI) {
        elemMap = new moab::Element::LinearTri(vertices);
      }
      */
      else if (j_type == MBEDGE) {
        elemMap = new moab::Element::LinearEdge(vertices);
      }
      else
        ERRORR("Unknown topology type.", MB_UNSUPPORTED_OPERATION);

      // Set the vertices in the Map and perform the integration
      try {
        /* VSM: Do we need this call ?? */
        // elemMap->set_vertices(vertices);

        // Perform the actual integration over the element
        intgr_val = elemMap->integrate_scalar_field(vfield);

        // Combine the result with those of the group
        grp_intrgr_val += intgr_val;
      }
      catch (moab::Element::Map::ArgError&) {
        std::cerr << "Failed to set vertices on Element::Map." << std::endl;
      }
      catch (moab::Element::Map::EvaluationError&) {
        std::cerr << "Failed to get inverse evaluation of coordinate on Element::Map." << std::endl;
      }

      delete(elemMap);
      free(vfield);
    }

    // Set the group integrated value in the vector
    integ_vals[i] = grp_intrgr_val;
  }

  return err;
}

// Apply a normalization factor to group of entities
ErrorCode Coupler::apply_group_norm_factor(std::vector<std::vector<EntityHandle> > &entity_sets,
                                           std::vector<double>                     &norm_factors,
                                           const char                              *norm_tag,
                                           Coupler::IntegType                      /*integ_type*/)
{
  ErrorCode err;

  // Construct the new tag for the normalization factor from the norm_tag name
  // and "_normf".
  int norm_tag_len = strlen(norm_tag);
  const char* normf_appd = "_normf";
  int normf_appd_len = strlen(normf_appd);

  char* normf_tag = (char*)malloc(norm_tag_len + normf_appd_len + 1);
  char* tmp_ptr = normf_tag;

  memcpy(tmp_ptr, norm_tag, norm_tag_len);
  tmp_ptr += norm_tag_len;
  memcpy(tmp_ptr, normf_appd, normf_appd_len);
  tmp_ptr += normf_appd_len;
  *tmp_ptr = '\0';

  Tag normf_hdl;
  // Check to see if the tag exists.  If not then create it and get the handle.
  err = mbImpl->tag_get_handle(normf_tag, 1, moab::MB_TYPE_DOUBLE, normf_hdl, moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT);
  ERRORR("Failed to create normalization factor tag.", err);
  if (normf_hdl == NULL) {
    std::string msg("Failed to create normalization factor tag named '");
    msg += std::string(normf_tag) + std::string("'");
    ERRORR(msg.c_str(), MB_FAILURE);
  }
  free(normf_tag);

  std::vector< std::vector<EntityHandle> >::iterator iter_i;
  std::vector<EntityHandle>::iterator iter_j;
  std::vector<double>::iterator iter_f;
  double grp_norm_factor = 0.0;

  // Loop over the entity sets
  for (iter_i = entity_sets.begin(), iter_f = norm_factors.begin();
       (iter_i != entity_sets.end()) && (iter_f != norm_factors.end());
       ++iter_i, ++iter_f) {
    grp_norm_factor = *iter_f;

    // Loop over the all the entity sets in iter_i and set the
    // new normf_tag with the norm factor value on each
    for (iter_j = (*iter_i).begin(); iter_j != (*iter_i).end(); ++iter_j) {
      EntityHandle entset = *iter_j;

      std::cout << "Coupler: applying normalization for entity set="
                << entset << ",  normalization_factor=" << grp_norm_factor << std::endl;

      err = mbImpl->tag_set_data(normf_hdl, &entset, 1, &grp_norm_factor);
      ERRORR("Failed to set normalization factor on entity set.", err);
    }
  }

  return MB_SUCCESS;
}

#define UINT_PER_X(X) ((sizeof(X) + sizeof(uint) - 1) / sizeof(uint))
#define UINT_PER_REAL UINT_PER_X(realType)
#define UINT_PER_LONG UINT_PER_X(slong)
#define UINT_PER_UNSIGNED UINT_PER_X(unsigned)

// Function for packing tuple_list.  Returns number of uints copied into buffer.
int pack_tuples(TupleList* tl, void** ptr)
{
  uint mi, ml, mul, mr;
  tl->getTupleSize(mi, ml, mul, mr);

  uint n = tl->get_n();

  int sz_buf = 1 + 4*UINT_PER_UNSIGNED +
               tl->get_n() * (mi +
                              ml*UINT_PER_LONG +
                              mul*UINT_PER_LONG +
                              mr*UINT_PER_REAL);

  uint *buf = (uint*) malloc(sz_buf*sizeof(uint));
  *ptr = (void*) buf;

  // Copy n
  memcpy(buf, &n,   sizeof(uint)),     buf += 1;
  // Copy mi
  memcpy(buf, &mi,  sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
  // Copy ml
  memcpy(buf, &ml,  sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
  // Copy mul
  memcpy(buf, &mul, sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
  // Copy mr
  memcpy(buf, &mr,  sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
  // Copy vi_wr
  memcpy(buf, tl->vi_rd,  tl->get_n()*mi*sizeof(sint)),   buf += tl->get_n()*mi;
  // Copy vl_wr
  memcpy(buf, tl->vl_rd,  tl->get_n()*ml*sizeof(slong)),  buf += tl->get_n()*ml*UINT_PER_LONG;
  // Copy vul_wr
  memcpy(buf, tl->vul_rd, tl->get_n()*mul*sizeof(ulong)), buf += tl->get_n()*mul*UINT_PER_LONG;
  // Copy vr_wr
  memcpy(buf, tl->vr_rd,  tl->get_n()*mr*sizeof(realType)),   buf += tl->get_n()*mr*UINT_PER_REAL;

  return sz_buf;
}

// Function for packing tuple_list
void unpack_tuples(void* ptr, TupleList** tlp)
{
  TupleList *tl = new TupleList();
  *tlp = tl;

  uint nt;
  unsigned mit, mlt, mult, mrt;
  uint *buf = (uint*)ptr;

  // Get n
  memcpy(&nt,   buf, sizeof(uint)),     buf += 1;
  // Get mi
  memcpy(&mit,  buf, sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
  // Get ml
  memcpy(&mlt,  buf, sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
  // Get mul
  memcpy(&mult, buf, sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
  // Get mr
  memcpy(&mrt,  buf, sizeof(unsigned)), buf += UINT_PER_UNSIGNED;

  // Initialize tl
  tl->initialize(mit, mlt, mult, mrt, nt);
  tl->enableWriteAccess();
  tl->set_n(nt);

  uint mi, ml, mul, mr;
  tl->getTupleSize(mi, ml, mul, mr);

  // Get vi_rd
  memcpy(tl->vi_wr,  buf, tl->get_n()*mi*sizeof(sint)),   buf += tl->get_n()*mi;
  // Get vl_rd
  memcpy(tl->vl_wr,  buf, tl->get_n()*ml*sizeof(slong)),  buf += tl->get_n()*ml*UINT_PER_LONG;
  // Get vul_rd
  memcpy(tl->vul_wr, buf, tl->get_n()*mul*sizeof(ulong)), buf += tl->get_n()*mul*UINT_PER_LONG;
  // Get vr_rd
  memcpy(tl->vr_wr,  buf, tl->get_n()*mr*sizeof(realType)),   buf += tl->get_n()*mr*UINT_PER_REAL;

  tl->disableWriteAccess();
  return;
}

ErrorCode Coupler::get_gl_points_on_elements(Range &targ_elems, std::vector<double> &vpos, int &numPointsOfInterest)
{
  numPointsOfInterest = targ_elems.size() * _ntot;
  vpos.resize(3 * numPointsOfInterest);
  int ielem = 0;
  for (Range::iterator eit = targ_elems.begin(); eit != targ_elems.end(); ++eit, ielem += _ntot * 3) {
    EntityHandle eh = *eit;
    const double * xval;
    const double * yval;
    const double * zval;
    ErrorCode rval = mbImpl-> tag_get_by_ptr(_xm1Tag, &eh, 1, (const void**)&xval);
    if (moab::MB_SUCCESS != rval) {
      std::cout << "Can't get xm1 values \n";
      return MB_FAILURE;
    }
    rval = mbImpl-> tag_get_by_ptr(_ym1Tag, &eh, 1, (const void**)&yval);
    if (moab::MB_SUCCESS != rval) {
      std::cout << "Can't get ym1 values \n";
      return MB_FAILURE;
    }
    rval = mbImpl-> tag_get_by_ptr(_zm1Tag, &eh, 1, (const void**)&zval);
    if (moab::MB_SUCCESS != rval) {
      std::cout << "Can't get zm1 values \n";
      return MB_FAILURE;
    }
    // Now, in a stride, populate vpos
    for (int i = 0; i < _ntot; i++) {
      vpos[ielem + 3*i    ] = xval[i];
      vpos[ielem + 3*i + 1] = yval[i];
      vpos[ielem + 3*i + 2] = zval[i];
    }
  }

  return MB_SUCCESS;
}

} // namespace_moab