Coupler.cpp

Go to the documentation of this file.

#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, false );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, false );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