DataCoupler.cpp

Go to the documentation of this file.

#include "DataCoupler.hpp"
#include "moab/Tree.hpp"
#include "moab/TupleList.hpp"
#include "moab/SpatialLocator.hpp"
#include "moab/ElemEvaluator.hpp"
#include "moab/Error.hpp"
 
#ifdef MOAB_HAVE_MPI
#include "moab/ParallelComm.hpp"
#endif
 
#include "iostream"
#include <algorithm>
#include <sstream>
 
#include <cstdio>
#include <cassert>
 
namespace moab
{
 
DataCoupler::DataCoupler( Interface* impl,
 Range& source_ents,
 int coupler_id,
 ParallelComm* pc,
 bool init_locator,
 int dim )
 : mbImpl( impl ), myPcomm( pc ), myId( coupler_id ), myDim( dim )
{
 assert( NULL != mbImpl && ( myPcomm || !source_ents.empty() ) );
 
 // Now initialize the tree
 if( init_locator )
 {
 myLocator = new SpatialLocator( mbImpl, source_ents );
 myLocator->elem_eval( new ElemEvaluator( mbImpl ) );
 
 // Initialize element evaluator with the default for the entity types in source_ents;
 // can be replaced later by application if desired
 if( !source_ents.empty() )
 {
 Range::pair_iterator pit = source_ents.pair_begin();
 EntityType last_type = MBMAXTYPE;
 for( ; pit != source_ents.pair_end(); ++pit )
 {
 EntityType this_type = mbImpl->type_from_handle( pit->first );
 if( last_type == this_type ) continue;
 ErrorCode rval = myLocator->elem_eval()->set_eval_set( pit->first );
 if( MB_SUCCESS != rval ) throw( rval );
 last_type = this_type;
 }
 }
 }
 
 if( -1 == dim && !source_ents.empty() ) myDim = mbImpl->dimension_from_handle( *source_ents.rbegin() );
}
 
/* Destructor
 */
DataCoupler::~DataCoupler()
{
 delete myLocator->elem_eval();
 delete myLocator;
}
 
ErrorCode DataCoupler::locate_points( Range& targ_ents,
 const double rel_iter_tol,
 const double abs_iter_tol,
 const double inside_tol )
{
 targetEnts = targ_ents;
 
#ifdef MOAB_HAVE_MPI
 if( myPcomm && myPcomm->size() > 1 )
 return myLocator->par_locate_points( myPcomm, targ_ents, rel_iter_tol, abs_iter_tol, inside_tol );
#endif
 
 return myLocator->locate_points( targ_ents, rel_iter_tol, abs_iter_tol, inside_tol );
}
 
ErrorCode DataCoupler::locate_points( double* xyz,
 int num_points,
 const double rel_iter_tol,
 const double abs_iter_tol,
 const double inside_tol )
{
#ifdef MOAB_HAVE_MPI
 if( myPcomm && myPcomm->size() > 1 )
 return myLocator->par_locate_points( myPcomm, xyz, num_points, rel_iter_tol, abs_iter_tol, inside_tol );
#endif
 
 return myLocator->locate_points( xyz, num_points, rel_iter_tol, abs_iter_tol, inside_tol );
}
 
ErrorCode DataCoupler::interpolate( /*DataCoupler::Method*/ int method,
 const std::string& interp_tag,
 double* interp_vals,
 std::vector< int >* point_indices,
 bool normalize )
{
 // Tag name input, translate to tag handle and pass down the chain
 
 Tag tag;
 ErrorCode result = mbImpl->tag_get_handle( interp_tag.c_str(), tag );MB_CHK_SET_ERR( result, "Failed to get handle for interpolation tag \"" << interp_tag << "\"" );
 
 return interpolate( method, tag, interp_vals, point_indices, normalize );
}
 
ErrorCode DataCoupler::interpolate( /*DataCoupler::Method*/ int* /* methods */,
 Tag* tags,
 int* points_per_method,
 int num_methods,
 double* interp_vals,
 std::vector< int >* point_indices,
 bool /*normalize*/ )
{
 // Lowest-level interpolate function, does actual interpolation using calls to ElemEvaluator
 
 ErrorCode result = MB_SUCCESS;
 
 unsigned int pts_total = 0;
 for( int i = 0; i < num_methods; i++ )
 pts_total += ( points_per_method ? points_per_method[i] : targetEnts.size() );
 
 unsigned int num_indices = ( point_indices ? point_indices->size() : targetEnts.size() );
 // # points and indices should be identical
 if( pts_total != num_indices ) return MB_FAILURE;
 
 // Since each tuple contains one interpolated tag, if we're interpolating multiple tags, every
 // tuple needs to be able to store up to max tag size
 int max_tsize = -1;
 for( int i = 0; i < num_methods; i++ )
 {
 int tmp_tsize;
 result = mbImpl->tag_get_length( tags[i], tmp_tsize );
 if( MB_SUCCESS != result ) return MB_FAILURE;
 max_tsize = std::max( max_tsize, tmp_tsize );
 }
 
 if( myPcomm && myPcomm->size() > 1 )
 {
 // Build up the tuple list to distribute from my target points; assume that
 // all procs use same method/tag input
 TupleList TLob; // TLob structure: (pto_i, ridx_i, lidx_i, meth_tagidx_i)
 TLob.initialize( 4, 0, 0, 0, num_indices );
 int tn = 0; // The tuple number we're currently on
 TLob.enableWriteAccess();
 for( int m = 0; m < num_methods; m++ )
 {
 int num_points = ( points_per_method ? points_per_method[m] : targetEnts.size() );
 for( int j = 0; j < num_points; j++ )
 {
 int idx = ( point_indices ? ( *point_indices )[j]
 : j ); // The index in my targetEnts for this interpolation point
 
 // Remote proc/idx from myLocator->parLocTable
 TLob.vi_wr[4 * tn] = myLocator->par_loc_table().vi_rd[2 * idx]; // proc
 TLob.vi_wr[4 * tn + 1] = myLocator->par_loc_table().vi_rd[2 * idx + 1]; // Remote idx
 
 // Local entity index, tag/method index from my data
 TLob.vi_wr[4 * tn + 2] = idx;
 TLob.vi_wr[4 * tn + 3] = m;
 TLob.inc_n();
 tn++;
 }
 }
 
 // Scatter/gather interpolation points
 myPcomm->proc_config().crystal_router()->gs_transfer( 1, TLob, 0 );
 
 // Perform interpolation on local source mesh; put results into TLinterp
 TupleList TLinterp; // TLinterp structure: (pto_i, ridx_i, vals[max_tsize]_d)
 TLinterp.initialize( 2, 0, 0, max_tsize, TLob.get_n() );
 TLinterp.set_n( TLob.get_n() ); // Set the size right away
 TLinterp.enableWriteAccess();
 
 for( unsigned int i = 0; i < TLob.get_n(); i++ )
 {
 int lidx = TLob.vi_rd[4 * i + 1]; // Index into myLocator's local table
 // Method and tag indexed with same index
 ///*Method*/ int method = methods[TLob.vi_rd[4*i + 3]];
 Tag tag = tags[TLob.vi_rd[4 * i + 3]];
 // Copy proc/remote index from TLob
 TLinterp.vi_wr[2 * i] = TLob.vi_rd[4 * i];
 TLinterp.vi_wr[2 * i + 1] = TLob.vi_rd[4 * i + 2];
 
 // Set up the evaluator for the tag and entity, then interpolate, putting result in
 // TLinterp
 myLocator->elem_eval()->set_tag_handle( tag );
 myLocator->elem_eval()->set_ent_handle( myLocator->loc_table().vul_rd[lidx] );
 result =
 myLocator->elem_eval()->eval( myLocator->loc_table().vr_rd + 3 * lidx, TLinterp.vr_rd + i * max_tsize );
 if( MB_SUCCESS != result ) return result;
 }
 
 // Scatter/gather interpolation data
 myPcomm->proc_config().crystal_router()->gs_transfer( 1, TLinterp, 0 );
 
 // Copy the interpolated field as a unit
 std::copy( TLinterp.vr_rd, TLinterp.vr_rd + TLinterp.get_n() * max_tsize, interp_vals );
 }
 else
 {
 // If called serially, interpolate directly from local locations/entities,
 // into either interp_vals or by setting data directly on tags on those entities
 std::vector< double > tmp_vals;
 std::vector< EntityHandle > tmp_ents;
 double* tmp_dbl = interp_vals;
 for( int i = 0; i < num_methods; i++ )
 {
 int num_points = ( points_per_method ? points_per_method[i] : targetEnts.size() );
 
 // Interpolated data is tsize long, which is either max size (if data passed back to
 // caller in interp_vals) or tag size (if data will be set on entities, in which case it
 // shouldn't have padding) set sizes here, inside loop over methods, to reduce memory
 // usage
 int tsize = max_tsize, tsize_bytes = 0;
 if( !interp_vals )
 {
 tmp_vals.resize( num_points * max_tsize );
 tmp_dbl = &tmp_vals[0];
 tmp_ents.resize( num_points );
 result = mbImpl->tag_get_length( tags[i], tsize );
 if( MB_SUCCESS != result ) return result;
 result = mbImpl->tag_get_bytes( tags[i], tsize_bytes );
 if( MB_SUCCESS != result ) return result;
 }
 
 for( int j = 0; j < num_points; j++ )
 {
 int lidx;
 if( point_indices )
 lidx = ( *point_indices )[j];
 else
 lidx = j;
 
 myLocator->elem_eval()->set_tag_handle( tags[i] );
 myLocator->elem_eval()->set_ent_handle( myLocator->loc_table().vul_rd[lidx] );
 if( !interp_vals ) tmp_ents[j] = targetEnts[lidx]; // Could be performance-sensitive, thus the if test
 result = myLocator->elem_eval()->eval( myLocator->loc_table().vr_rd + 3 * lidx, tmp_dbl );
 tmp_dbl += tsize;
 if( MB_SUCCESS != result ) return result;
 } // for j over points
 
 if( !interp_vals )
 {
 // Set tags on tmp_ents; data is already w/o padding, due to tsize setting above
 result = mbImpl->tag_set_data( tags[i], &tmp_ents[0], tmp_ents.size(), &tmp_vals[0] );
 if( MB_SUCCESS != result ) return result;
 }
 
 } // for i over methods
 } // if myPcomm
 
 // Done
 return MB_SUCCESS;
}
 
} // namespace moab