// MOAB includes
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
#include "moab/Core.hpp"
#include "Coupler.hpp"
#include "moab_mpi.h"
#include "ElemUtil.hpp"

// C++ includes
#include <iostream>
#include <iomanip>
#include <sstream>
#include <assert.h>

using namespace moab;

#define STRINGIFY_(A) #A
#define STRINGIFY(A) STRINGIFY_(A)
#ifdef MESHDIR
std::string TestDir(STRINGIFY(MESHDIR));
#else
std::string TestDir(".");
#endif

/*
  Sample usages:
    1) P_0 interpolation: ./mbcoupler_test -meshes <src_mesh> <target_mesh> -itag <interp_tag> -meth 0 -outfile <output>
    2) P_1 interpolation: ./mbcoupler_test -meshes <src_mesh> <target_mesh> -itag <interp_tag> -meth 1 -outfile <output>
    3) P_1 interpolation with epsilon control: ./mbcoupler_test -meshes <src_mesh> <target_mesh> -itag <interp_tag> -meth 1 -eps <tolerance for locating points; say 0.01> -outfile <output>
    3) P_0 interpolation with global normalization: ./mbcoupler_test -meshes <src_mesh> <target_mesh> -itag <interp_tag> -meth 0 -gnorm <gnorm_tag_name> -outfile <output>
    4) P_1 interpolation with subset normalization: ./mbcoupler_test -meshes <src_mesh> <target_mesh> -itag <interp_tag> -meth 1 -ssnorm <snorm_tag_name> <snorm_criteria: MATERIAL_SET> -outfile <output>
*/

// Print usage
void print_usage()
{
  std::cerr << "Usage: ";
  std::cerr << "mbcoupler_test -meshes <source_mesh> <target_mesh> -itag <interp_tag> [-gnorm <gnorm_tag>] [-ssnorm <ssnorm_tag> <ssnorm_selection>] [-ropts <roptions>] [-outfile <out_file> [-wopts <woptions>]] [-dbgout [<dbg_file>]]" << std::endl;
  std::cerr << "    -meshes" << std::endl;
  std::cerr << "        Read in mesh files <source_mesh> and <target_mesh>." << std::endl;
  std::cerr << "    -itag" << std::endl;
  std::cerr << "        Interpolate tag <interp_tag> from source mesh to target mesh." << std::endl;
  std::cerr << "    -gnorm" << std::endl;
  std::cerr << "        Normalize the value of tag <gnorm_tag> over then entire mesh and save to" << std::endl;
  std::cerr << "        tag \"<gnorm_tag>_normf\" on the mesh set.  Do this for all meshes." << std::endl;
  std::cerr << "    -ssnorm" << std::endl;
  std::cerr << "        Normalize the value of tag <ssnorm_tag> over subsets of a mesh and save to" << std::endl;
  std::cerr << "        tag \"<ssnorm_tag>_normf\" on the Entity Set for each subset.  Subsets are selected" << std::endl;
  std::cerr << "        using criteria in <ssnorm_selection>.  Do this for all meshes." << std::endl;
  std::cerr << "    -ropts" << std::endl;
  std::cerr << "        Read in the mesh files using options in <roptions>." << std::endl;
  std::cerr << "    -outfile" << std::endl;
  std::cerr << "        Write out target mesh to <out_file>." << std::endl;
  std::cerr << "    -wopts" << std::endl;
  std::cerr << "        Write out mesh files using options in <woptions>." << std::endl;
  std::cerr << "    -dbgout" << std::endl;
  std::cerr << "        Write stdout and stderr streams to the file \'<dbg_file>.txt\'." << std::endl;
  std::cerr << "    -eps" << std::endl;
  std::cerr << "        epsilon" << std::endl;
  std::cerr << "    -meth <method> (0=CONSTANT, 1=LINEAR_FE, 2=QUADRATIC_FE, 3=SPECTRAL)" << std::endl;
}

#ifdef MOAB_HAVE_HDF5

ErrorCode get_file_options(int argc, char **argv, int rank,
                           std::vector<std::string> &meshFiles,
                           Coupler::Method &method,
                           std::string &interpTag,
                           std::string &gNormTag,
                           std::string &ssNormTag,
                           std::vector<const char*> &ssTagNames,
                           std::vector<const char*> &ssTagValues,
                           std::string &readOpts,
                           std::string &outFile,
                           std::string &writeOpts,
                           std::string &dbgFile,
                           bool &help,
                           double &epsilon);

ErrorCode report_iface_ents(Interface *mbImpl,
                            std::vector<ParallelComm*> &pcs,
                            bool print_results);

ErrorCode test_interpolation(Interface *mbImpl,
                             Coupler::Method method,
                             std::string &interpTag,
                             std::string &gNormTag,
                             std::string &ssNormTag,
                             std::vector<const char*> &ssTagNames,
                             std::vector<const char*> &ssTagValues,
                             EntityHandle *roots,
                             std::vector<ParallelComm*> &pcs,
                             double &instant_time,
                             double &pointloc_time,
                             double &interp_time,
                             double &gnorm_time,
                             double &ssnorm_time,
                             double &toler);

void reduceMax(double &v)
{
  double buf;

  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Allreduce(&v, &buf, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);

  v = buf;
}

int main(int argc, char **argv)
{
  // Need to init MPI first, to tell how many procs and rank
  int ierr = MPI_Init(&argc, &argv);

  std::vector<const char*> ssTagNames, ssTagValues;
  std::vector<std::string> meshFiles;
  std::string interpTag, gNormTag, ssNormTag, readOpts, outFile, writeOpts, dbgFile;
  Coupler::Method method = Coupler::CONSTANT;

  ErrorCode result;
  bool help = false;
  double toler = 5.e-10;
  int nprocs, rank;
  ierr = MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
  ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);

  result = get_file_options(argc, argv, rank, meshFiles, method, interpTag,
                            gNormTag, ssNormTag, ssTagNames, ssTagValues,
                            readOpts, outFile, writeOpts, dbgFile, help, toler);

  if (result != MB_SUCCESS || help) {
    print_usage();
    ierr = MPI_Finalize();
    return ierr;
  }

  // Redirect stdout and stderr if dbgFile is not null
  if (!dbgFile.empty()) {
    std::stringstream dfname;
    dfname << dbgFile << rank << ".txt";
    if (!std::freopen(dfname.str().c_str(), "a", stdout))
      return false;
    if (!std::freopen(dfname.str().c_str(), "a", stderr))
      return false;
  }

  // Create MOAB instance based on that
  Interface *mbImpl = new (std::nothrow) Core();
  if (NULL == mbImpl)
    return 1;

  // Read in mesh(es)
  std::vector<ParallelComm*> pcs(meshFiles.size());

  // Create root sets for each mesh using the iMesh API. Then pass these
  // to the load_file functions to be populated.
  EntityHandle *roots = (EntityHandle *)malloc(sizeof(EntityHandle) * meshFiles.size());

  for (unsigned int i = 0; i < meshFiles.size(); i++) {
    pcs[i] = new ParallelComm(mbImpl, MPI_COMM_WORLD);
    int index = pcs[i]->get_id();
    std::string newReadopts;
    std::ostringstream extraOpt;
    extraOpt << ";PARALLEL_COMM=" << index;
    newReadopts = readOpts+extraOpt.str();

    result = mbImpl->create_meshset(MESHSET_SET, roots[i]);MB_CHK_ERR(result);

    result = mbImpl->load_file( meshFiles[i].c_str(), &roots[i], newReadopts.c_str() );MB_CHK_ERR(result);
    //result = rps[i]->load_file(meshFiles[i].c_str(), &roots[i], FileOptions(readOpts.c_str()));MB_CHK_ERR(result);
  }

  result = report_iface_ents(mbImpl, pcs, true);MB_CHK_ERR(result);

  double instant_time = 0.0, pointloc_time = 0.0, interp_time = 0.0, gnorm_time = 0.0, ssnorm_time = 0.0;
  // Test interpolation and global normalization and subset normalization

  result = test_interpolation(mbImpl, method, interpTag, gNormTag, ssNormTag,
                              ssTagNames, ssTagValues, roots, pcs,
                              instant_time, pointloc_time, interp_time,
                              gnorm_time, ssnorm_time, toler);MB_CHK_ERR(result);

  reduceMax(instant_time);
  reduceMax(pointloc_time);
  reduceMax(interp_time);

  if (0 == rank)
    printf("\nMax time : %g %g %g (inst loc interp -- %d procs)\n", instant_time,
           pointloc_time, interp_time, nprocs);

  // Output mesh
  if (!outFile.empty()) {
    Range partSets;
    // Only save the target mesh
    partSets.insert((EntityHandle)roots[1]);
    std::string newwriteOpts;
    std::ostringstream extraOpt;
    extraOpt << ";PARALLEL_COMM=" << 1;
    newwriteOpts = writeOpts + extraOpt.str();
    result = mbImpl->write_file(outFile.c_str(), NULL, newwriteOpts.c_str(), partSets);MB_CHK_ERR(result);
    if(0==rank)
    {
      std::cout << "Wrote " << outFile << std::endl;
      std::cout << "mbcoupler_test complete." << std::endl;
    }
  }

  for (unsigned int i = 0; i < meshFiles.size(); i++)
    delete pcs[i];

  delete mbImpl;
  
  ierr = MPI_Finalize();
  return 0;
}

ErrorCode report_iface_ents(Interface *mbImpl,
                            std::vector<ParallelComm*> &pcs,
                            const bool print_results)
{
  Range iface_ents[6];
  ErrorCode result = MB_SUCCESS, tmp_result;

  // Now figure out which vertices are shared
  for (unsigned int p = 0; p < pcs.size(); p++) {
    for (int i = 0; i < 4; i++) {
      tmp_result = pcs[p]->get_iface_entities(-1, i, iface_ents[i]);

      if (MB_SUCCESS != tmp_result) {
        std::cerr << "get_iface_entities returned error on proc "
                  << pcs[p]->proc_config().proc_rank() << "; message: " << std::endl;
        std::string last_error;
        result = mbImpl->get_last_error(last_error);
        if (last_error.empty())
          std::cerr << "(none)" << std::endl;
        else
          std::cerr << last_error << std::endl;
        result = tmp_result;
      }
      if (0 != i)
        iface_ents[4].merge(iface_ents[i]);
    }
  }

  // Report # iface entities
  result = mbImpl->get_adjacencies(iface_ents[4], 0, false, iface_ents[5],
                                   Interface::UNION);MB_CHK_ERR(result);

  int rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);

  if (print_results || iface_ents[0].size() != iface_ents[5].size()) {
    std::cerr << "Proc " << rank << " iface entities: " << std::endl;
    for (int i = 0; i < 4; i++)
      std::cerr << "    " << iface_ents[i].size() << " "
                << i << "d iface entities." << std::endl;
    std::cerr << "    (" << iface_ents[5].size()
              << " verts adj to other iface ents)" << std::endl;
  }

  return result;
}

// Check first character for a '-'.
// Return true if one is found. False otherwise.
bool check_for_flag(const char *str) {
  if ('-' == str[0])
    return true;
  else
    return false;
}

// get_file_options() function with added possibilities for mbcoupler_test.
ErrorCode get_file_options(int argc, char **argv, int rank,
                           std::vector<std::string> &meshFiles,
                           Coupler::Method &method,
                           std::string &interpTag,
                           std::string &gNormTag,
                           std::string &ssNormTag,
                           std::vector<const char*> &ssTagNames,
                           std::vector<const char*> &ssTagValues,
                           std::string &readOpts,
                           std::string &outFile,
                           std::string &writeOpts,
                           std::string &dbgFile,
                           bool &help,
                           double &epsilon)
{
  // Initialize some of the outputs to null values indicating not present
  // in the argument list.
  gNormTag = "";
  ssNormTag = "";
  readOpts = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=3.0.1;CPUTIME";
  outFile = "";
  writeOpts = "PARALLEL=WRITE_PART;CPUTIME";
  dbgFile = "";
  std::string defaultDbgFile = argv[0]; // The executable name will be the default debug output file.

  // These will indicate if we've gotten our required parameters at the end of parsing.
  bool haveMeshes = false;
  bool haveInterpTag = false;

  // Loop over the values in argv pulling out an parsing each one
  int npos = 1;

  if (argc > 1 && argv[1] == std::string("-h")) {
    help = true;
    return MB_SUCCESS;
  }

  while (npos < argc) {
    if (argv[npos] == std::string("-meshes")) {
      // Parse out the mesh filenames
      npos++;
      int numFiles = 2;
      meshFiles.resize(numFiles);
      for (int i = 0; i < numFiles; i++) {
        if ((npos < argc) && (!check_for_flag(argv[npos])))
          meshFiles[i] = argv[npos++];
        else {
          std::cerr << "    ERROR - missing correct number of mesh filenames" << std::endl;
          return MB_FAILURE;
        }
      }

      haveMeshes = true;
    }
    else if (argv[npos] == std::string("-itag")) {
      // Parse out the interpolation tag
      npos++;
      if ((npos < argc) && (!check_for_flag(argv[npos])))
        interpTag = argv[npos++];
      else {
        std::cerr << "    ERROR - missing <interp_tag>" << std::endl;
        return MB_FAILURE;
      }

      haveInterpTag = true;
    }
    else if (argv[npos] == std::string("-meth")) {
      // Parse out the interpolation tag
      npos++;
      if (argv[npos][0] == '0')
        method = Coupler::CONSTANT;
      else if (argv[npos][0] == '1')
        method = Coupler::LINEAR_FE;
      else if (argv[npos][0] == '2')
        method = Coupler::QUADRATIC_FE;
      else if (argv[npos][0] == '3')
        method = Coupler::SPECTRAL;
      else {
        std::cerr << "    ERROR - unrecognized method number " << method << std::endl;
        return MB_FAILURE;
      }
      npos++;
    }
    else if (argv[npos] == std::string("-eps")) {
      // Parse out the tolerance
      npos++;
      if ((npos < argc) && (!check_for_flag(argv[npos])))
        epsilon = atof(argv[npos++]);
      else {
        std::cerr << "    ERROR - missing <epsilon>" << std::endl;
        return MB_FAILURE;
      }
    }
    else if (argv[npos] == std::string("-gnorm")) {
      // Parse out the global normalization tag
      npos++;
      if ((npos < argc) && (!check_for_flag(argv[npos])))
        gNormTag = argv[npos++];
      else {
        std::cerr << "    ERROR - missing <gnorm_tag>" << std::endl;
        return MB_FAILURE;
      }
    }
    else if (argv[npos] == std::string("-ssnorm")) {
      // Parse out the subset normalization tag and selection criteria
      npos++;
      if ((npos < argc) && (!check_for_flag(argv[npos])))
        ssNormTag = argv[npos++];
      else {
        std::cerr << "    ERROR - missing <ssnorm_tag>" << std::endl;
        return MB_FAILURE;
      }

      if ((npos < argc) && (!check_for_flag(argv[npos]))) {
        char* opts = argv[npos++];
        char sep1[1] = {';'};
        char sep2[1] = {'='};
        bool end_vals_seen = false;
        std::vector<char*> tmpTagOpts;

        // First get the options
        for (char* i = strtok(opts, sep1); i; i = strtok(0, sep1))
          tmpTagOpts.push_back(i);

        // Parse out the name and val or just name.
        for (unsigned int j = 0; j < tmpTagOpts.size(); j++) {
          char* e = strtok(tmpTagOpts[j], sep2);
          ssTagNames.push_back(e);
          e = strtok(0, sep2);
          if (e != NULL) {
            // We have a value
            if (end_vals_seen) {
              // ERROR we should not have a value after none are seen
              std::cerr << "    ERROR - new value seen after end of values in <ssnorm_selection>" << std::endl;
              return MB_FAILURE;
            }
            // Otherwise get the value string from e and convert it to an int
            int *valp = new int;
            *valp = atoi(e);
            ssTagValues.push_back((const char*)valp);
          }
          else {
            // Otherwise there is no '=' so push a null on the list
            end_vals_seen = true;
            ssTagValues.push_back((const char*)0);
          }
        }
      }
      else {
        std::cerr << "    ERROR - missing <ssnorm_selection>" << std::endl;
        return MB_FAILURE;
      }
    }
    else if (argv[npos] == std::string("-ropts")) {
      // Parse out the mesh file read options
      npos++;
      if ((npos < argc) && (!check_for_flag(argv[npos])))
        readOpts = argv[npos++];
      else {
        std::cerr << "    ERROR - missing <roptions>" << std::endl;
        return MB_FAILURE;
      }
    }
    else if (argv[npos] == std::string("-outfile")) {
      // Parse out the output file name
      npos++;
      if ((npos < argc) && (!check_for_flag(argv[npos])))
        outFile = argv[npos++];
      else {
        std::cerr << "    ERROR - missing <out_file>" << std::endl;
        return MB_FAILURE;
      }
    }
    else if (argv[npos] == std::string("-wopts")) {
      // Parse out the output file write options
      npos++;
      if ((npos < argc) && (!check_for_flag(argv[npos])))
        writeOpts = argv[npos++];
      else {
        std::cerr << "    ERROR - missing <woptions>" << std::endl;
        return MB_FAILURE;
      }
    }
    else if (argv[npos] == std::string("-dbgout")) {
      // Parse out the debug output file name.
      // If no name then use the default.
      npos++;
      if ((npos < argc) && (!check_for_flag(argv[npos])))
        dbgFile = argv[npos++];
      else
        dbgFile = defaultDbgFile;
    }
    else {
      // Unrecognized parameter.  Skip it and move along.
      std::cerr << "    ERROR - Unrecognized parameter:" << argv[npos] << std::endl;
      std::cerr << "            Skipping..." << std::endl;
      npos++;
    }
  }

  if (!haveMeshes) {
    meshFiles.resize(2);
    meshFiles[0] = std::string(TestDir + "/64bricks_1khex.h5m");
    meshFiles[1] = std::string(TestDir + "/64bricks_12ktet.h5m");
    if (0==rank)
      std::cout << "Mesh files not entered; using default files "
              << meshFiles[0] << " and " << meshFiles[1] << std::endl;
  }

  if (!haveInterpTag) {
    interpTag = "vertex_field";
    std::cout << "Interpolation field name not given, using default of " << interpTag << std::endl;
  }

#ifdef MOAB_HAVE_HDF5
  if (1 == argc) {
    if (0==rank) std::cout << "No arguments given; using output file dum.h5m." << std::endl;
    outFile = "dum.h5m";
  }
#endif

  return MB_SUCCESS;
}

// End get_file_options()


ErrorCode test_interpolation(Interface *mbImpl,
                             Coupler::Method method,
                             std::string &interpTag,
                             std::string &gNormTag,
                             std::string &ssNormTag,
                             std::vector<const char*> &ssTagNames,
                             std::vector<const char*> &ssTagValues,
                             EntityHandle *roots,
                             std::vector<ParallelComm*> &pcs,
                             double &instant_time,
                             double &pointloc_time,
                             double &interp_time,
                             double &gnorm_time,
                             double &ssnorm_time,
                             double &toler)
{
  assert(method >= Coupler::CONSTANT && method <= Coupler::SPECTRAL);

  // Source is 1st mesh, target is 2nd
  Range src_elems, targ_elems, targ_verts;
  ErrorCode result = pcs[0]->get_part_entities(src_elems, 3);MB_CHK_ERR(result);

  double start_time = MPI_Wtime();

  // Instantiate a coupler, which also initializes the tree
  Coupler mbc(mbImpl, pcs[0], src_elems, 0);

  // Initialize spectral elements, if they exist
  bool specSou=false, specTar = false;
  //result =  mbc.initialize_spectral_elements(roots[0], roots[1], specSou, specTar);

  instant_time = MPI_Wtime();

  // Get points from the target mesh to interpolate
  // We have to treat differently the case when the target is a spectral mesh
  // In that case, the points of interest are the GL points, not the vertex nodes
  std::vector<double> vpos; // This will have the positions we are interested in
  int numPointsOfInterest = 0;
  if (!specTar) { // Usual case
    Range tmp_verts;

    // First get all vertices adj to partition entities in target mesh
    result = pcs[1]->get_part_entities(targ_elems, 3);MB_CHK_ERR(result);
    if (Coupler::CONSTANT == method)
      targ_verts = targ_elems;
    else
      result = mbImpl->get_adjacencies(targ_elems, 0, false, targ_verts,
                                       Interface::UNION);MB_CHK_ERR(result);

    // Then get non-owned verts and subtract
    result = pcs[1]->get_pstatus_entities(0, PSTATUS_NOT_OWNED, tmp_verts);MB_CHK_ERR(result);
    targ_verts = subtract(targ_verts, tmp_verts);
    // get position of these entities; these are the target points
    numPointsOfInterest = (int)targ_verts.size();
    vpos.resize(3*targ_verts.size());
    result = mbImpl->get_coords(targ_verts, &vpos[0]);MB_CHK_ERR(result);
  }
  else {
    // In this case, the target mesh is spectral, we want values
    // interpolated on the GL positions; for each element, get the GL points, and construct CartVect!!!
    result = pcs[1]->get_part_entities(targ_elems, 3);MB_CHK_ERR(result);
    result = mbc.get_gl_points_on_elements(targ_elems, vpos, numPointsOfInterest);MB_CHK_ERR(result);
  }

  // Locate those points in the source mesh
  std::cout<<"rank "<< pcs[0]->proc_config().proc_rank() << " points of interest: " << numPointsOfInterest << "\n";
  result = mbc.locate_points(&vpos[0], numPointsOfInterest, 0, toler);MB_CHK_ERR(result);

  pointloc_time = MPI_Wtime();

  // Now interpolate tag onto target points
  std::vector<double> field(numPointsOfInterest);

  result = mbc.interpolate(method, interpTag, &field[0]);MB_CHK_ERR(result);

  interp_time = MPI_Wtime();

  // Do global normalization if specified
  if (!gNormTag.empty()) {
    // Normalize the source mesh
    result = mbc.normalize_mesh(roots[0], gNormTag.c_str(), Coupler::VOLUME, 4);MB_CHK_ERR(result);

    // Normalize the target mesh
    result = mbc.normalize_mesh(roots[1], gNormTag.c_str(), Coupler::VOLUME, 4);MB_CHK_ERR(result);
  }

  gnorm_time = MPI_Wtime();

  // Do subset normalization if specified

  if (!ssNormTag.empty()) {

    result = mbc.normalize_subset(roots[0],
                               ssNormTag.c_str(), 
                               &ssTagNames[0], 
                               ssTagNames.size(), 
                               &ssTagValues[0], 
                               Coupler::VOLUME, 
                               4);MB_CHK_ERR(result);

    result = mbc.normalize_subset(roots[1],
                               ssNormTag.c_str(), 
                               &ssTagNames[0], 
                               ssTagNames.size(), 
                               &ssTagValues[0], 
                               Coupler::VOLUME, 
                               4);MB_CHK_ERR(result);
  }

  ssnorm_time = MPI_Wtime();

  ssnorm_time -= gnorm_time;
  gnorm_time -= interp_time;
  interp_time -= pointloc_time;
  pointloc_time -= instant_time;
  instant_time -= start_time;

  // Set field values as tag on target vertices
  if (specSou) {
    // Create a new tag for the values on the target
    Tag tag;
    std::string newtag = interpTag + "_TAR";
    result = mbImpl->tag_get_handle(newtag.c_str(), 1, MB_TYPE_DOUBLE, tag, MB_TAG_CREAT | MB_TAG_DENSE);MB_CHK_ERR(result);
    result = mbImpl->tag_set_data(tag, targ_verts, &field[0]);MB_CHK_ERR(result);
  }
  else {
    if (!specTar) {
      // Use original tag
      Tag tag;
      result = mbImpl->tag_get_handle(interpTag.c_str(), 1, MB_TYPE_DOUBLE, tag);MB_CHK_ERR(result);
      result = mbImpl->tag_set_data(tag, targ_verts, &field[0]);MB_CHK_ERR(result);
    }
    else {
      // We have the field values computed at each GL points, on each element
      // in the target mesh
      // We need to create a new tag, on elements, of size _ntot, to hold
      // all those values.
      // So it turns out we need ntot. Maybe we can compute it from the
      // number of values computed, divided by number of elements
      int ntot = numPointsOfInterest / targ_elems.size();
      Tag tag;
      std::string newtag = interpTag + "_TAR";
      result = mbImpl->tag_get_handle(newtag.c_str(), ntot, MB_TYPE_DOUBLE, tag, MB_TAG_CREAT | MB_TAG_DENSE);MB_CHK_ERR(result);
      result = mbImpl->tag_set_data(tag, targ_elems, &field[0]);MB_CHK_ERR(result);
    }
  }

  // Done
  return MB_SUCCESS;
}

#else

int main(int /*argc*/, char** /*argv*/)
{
  print_usage();
  return 0;
}

#endif
