MOAB: Mesh Oriented datABase  (version 5.2.1)
DataCoupler.cpp
Go to the documentation of this file.
00001 #include "DataCoupler.hpp"
00002 #include "moab/Tree.hpp"
00003 #include "moab/TupleList.hpp"
00004 #include "moab/SpatialLocator.hpp"
00005 #include "moab/ElemEvaluator.hpp"
00006 #include "moab/Error.hpp"
00007 
00008 #ifdef MOAB_HAVE_MPI
00009 #include "moab/ParallelComm.hpp"
00010 #endif
00011 
00012 #include "iostream"
00013 #include <algorithm>
00014 #include <sstream>
00015 
00016 #include <stdio.h>
00017 #include "assert.h"
00018 
00019 namespace moab
00020 {
00021 
00022 DataCoupler::DataCoupler( Interface* impl, Range& source_ents, int coupler_id, ParallelComm* pc, bool init_locator,
00023                           int dim )
00024     : mbImpl( impl ), myPcomm( pc ), myId( coupler_id ), myDim( dim )
00025 {
00026     assert( NULL != mbImpl && ( myPcomm || !source_ents.empty() ) );
00027 
00028     // Now initialize the tree
00029     if( init_locator )
00030     {
00031         myLocator = new SpatialLocator( mbImpl, source_ents );
00032         myLocator->elem_eval( new ElemEvaluator( mbImpl ) );
00033 
00034         // Initialize element evaluator with the default for the entity types in source_ents;
00035         // can be replaced later by application if desired
00036         if( !source_ents.empty() )
00037         {
00038             Range::pair_iterator pit = source_ents.pair_begin();
00039             EntityType last_type     = MBMAXTYPE;
00040             for( ; pit != source_ents.pair_end(); ++pit )
00041             {
00042                 EntityType this_type = mbImpl->type_from_handle( pit->first );
00043                 if( last_type == this_type ) continue;
00044                 ErrorCode rval = myLocator->elem_eval()->set_eval_set( pit->first );
00045                 if( MB_SUCCESS != rval ) throw( rval );
00046                 last_type = this_type;
00047             }
00048         }
00049     }
00050 
00051     if( -1 == dim && !source_ents.empty() ) myDim = mbImpl->dimension_from_handle( *source_ents.rbegin() );
00052 }
00053 
00054 /* Destructor
00055  */
00056 DataCoupler::~DataCoupler()
00057 {
00058     delete myLocator->elem_eval();
00059     delete myLocator;
00060 }
00061 
00062 ErrorCode DataCoupler::locate_points( Range& targ_ents, const double rel_iter_tol, const double abs_iter_tol,
00063                                       const double inside_tol )
00064 {
00065     targetEnts = targ_ents;
00066 
00067 #ifdef MOAB_HAVE_MPI
00068     if( myPcomm && myPcomm->size() > 1 )
00069         return myLocator->par_locate_points( myPcomm, targ_ents, rel_iter_tol, abs_iter_tol, inside_tol );
00070 #endif
00071 
00072     return myLocator->locate_points( targ_ents, rel_iter_tol, abs_iter_tol, inside_tol );
00073 }
00074 
00075 ErrorCode DataCoupler::locate_points( double* xyz, int num_points, const double rel_iter_tol, const double abs_iter_tol,
00076                                       const double inside_tol )
00077 {
00078 #ifdef MOAB_HAVE_MPI
00079     if( myPcomm && myPcomm->size() > 1 )
00080         return myLocator->par_locate_points( myPcomm, xyz, num_points, rel_iter_tol, abs_iter_tol, inside_tol );
00081 #endif
00082 
00083     return myLocator->locate_points( xyz, num_points, rel_iter_tol, abs_iter_tol, inside_tol );
00084 }
00085 
00086 ErrorCode DataCoupler::interpolate( /*DataCoupler::Method*/ int method, const std::string& interp_tag,
00087                                     double* interp_vals, std::vector< int >* point_indices, bool normalize )
00088 {
00089     // Tag name input, translate to tag handle and pass down the chain
00090 
00091     Tag tag;
00092     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 << "\"" );
00093 
00094     return interpolate( method, tag, interp_vals, point_indices, normalize );
00095 }
00096 
00097 ErrorCode DataCoupler::interpolate( /*DataCoupler::Method*/ int* /* methods */, Tag* tags, int* points_per_method,
00098                                     int num_methods, double* interp_vals, std::vector< int >* point_indices,
00099                                     bool /*normalize*/ )
00100 {
00101     // Lowest-level interpolate function, does actual interpolation using calls to ElemEvaluator
00102 
00103     ErrorCode result = MB_SUCCESS;
00104 
00105     unsigned int pts_total = 0;
00106     for( int i = 0; i < num_methods; i++ )
00107         pts_total += ( points_per_method ? points_per_method[i] : targetEnts.size() );
00108 
00109     unsigned int num_indices = ( point_indices ? point_indices->size() : targetEnts.size() );
00110     // # points and indices should be identical
00111     if( pts_total != num_indices ) return MB_FAILURE;
00112 
00113     // Since each tuple contains one interpolated tag, if we're interpolating multiple tags, every
00114     // tuple needs to be able to store up to max tag size
00115     int max_tsize = -1;
00116     for( int i = 0; i < num_methods; i++ )
00117     {
00118         int tmp_tsize;
00119         result = mbImpl->tag_get_length( tags[i], tmp_tsize );
00120         if( MB_SUCCESS != result ) return MB_FAILURE;
00121         max_tsize = std::max( max_tsize, tmp_tsize );
00122     }
00123 
00124     if( myPcomm && myPcomm->size() > 1 )
00125     {
00126         // Build up the tuple list to distribute from my target points; assume that
00127         // all procs use same method/tag input
00128         TupleList TLob;  // TLob structure: (pto_i, ridx_i, lidx_i, meth_tagidx_i)
00129         TLob.initialize( 4, 0, 0, 0, num_indices );
00130         int tn = 0;  // The tuple number we're currently on
00131         TLob.enableWriteAccess();
00132         for( int m = 0; m < num_methods; m++ )
00133         {
00134             int num_points = ( points_per_method ? points_per_method[m] : targetEnts.size() );
00135             for( int j = 0; j < num_points; j++ )
00136             {
00137                 int idx = ( point_indices ? ( *point_indices )[j]
00138                                           : j );  // The index in my targetEnts for this interpolation point
00139 
00140                 // Remote proc/idx from myLocator->parLocTable
00141                 TLob.vi_wr[4 * tn]     = myLocator->par_loc_table().vi_rd[2 * idx];      // proc
00142                 TLob.vi_wr[4 * tn + 1] = myLocator->par_loc_table().vi_rd[2 * idx + 1];  // Remote idx
00143 
00144                 // Local entity index, tag/method index from my data
00145                 TLob.vi_wr[4 * tn + 2] = idx;
00146                 TLob.vi_wr[4 * tn + 3] = m;
00147                 TLob.inc_n();
00148                 tn++;
00149             }
00150         }
00151 
00152         // Scatter/gather interpolation points
00153         myPcomm->proc_config().crystal_router()->gs_transfer( 1, TLob, 0 );
00154 
00155         // Perform interpolation on local source mesh; put results into TLinterp
00156         TupleList TLinterp;  // TLinterp structure: (pto_i, ridx_i, vals[max_tsize]_d)
00157         TLinterp.initialize( 2, 0, 0, max_tsize, TLob.get_n() );
00158         TLinterp.set_n( TLob.get_n() );  // Set the size right away
00159         TLinterp.enableWriteAccess();
00160 
00161         for( unsigned int i = 0; i < TLob.get_n(); i++ )
00162         {
00163             int lidx = TLob.vi_rd[4 * i + 1];  // Index into myLocator's local table
00164             // Method and tag indexed with same index
00165             ///*Method*/ int method = methods[TLob.vi_rd[4*i + 3]];
00166             Tag tag = tags[TLob.vi_rd[4 * i + 3]];
00167             // Copy proc/remote index from TLob
00168             TLinterp.vi_wr[2 * i]     = TLob.vi_rd[4 * i];
00169             TLinterp.vi_wr[2 * i + 1] = TLob.vi_rd[4 * i + 2];
00170 
00171             // Set up the evaluator for the tag and entity, then interpolate, putting result in
00172             // TLinterp
00173             myLocator->elem_eval()->set_tag_handle( tag );
00174             myLocator->elem_eval()->set_ent_handle( myLocator->loc_table().vul_rd[lidx] );
00175             result =
00176                 myLocator->elem_eval()->eval( myLocator->loc_table().vr_rd + 3 * lidx, TLinterp.vr_rd + i * max_tsize );
00177             if( MB_SUCCESS != result ) return result;
00178         }
00179 
00180         // Scatter/gather interpolation data
00181         myPcomm->proc_config().crystal_router()->gs_transfer( 1, TLinterp, 0 );
00182 
00183         // Copy the interpolated field as a unit
00184         std::copy( TLinterp.vr_rd, TLinterp.vr_rd + TLinterp.get_n() * max_tsize, interp_vals );
00185     }
00186     else
00187     {
00188         // If called serially, interpolate directly from local locations/entities,
00189         // into either interp_vals or by setting data directly on tags on those entities
00190         std::vector< double > tmp_vals;
00191         std::vector< EntityHandle > tmp_ents;
00192         double* tmp_dbl = interp_vals;
00193         for( int i = 0; i < num_methods; i++ )
00194         {
00195             int num_points = ( points_per_method ? points_per_method[i] : targetEnts.size() );
00196 
00197             // Interpolated data is tsize long, which is either max size (if data passed back to
00198             // caller in interp_vals) or tag size (if data will be set on entities, in which case it
00199             // shouldn't have padding) set sizes here, inside loop over methods, to reduce memory
00200             // usage
00201             int tsize = max_tsize, tsize_bytes = 0;
00202             if( !interp_vals )
00203             {
00204                 tmp_vals.resize( num_points * max_tsize );
00205                 tmp_dbl = &tmp_vals[0];
00206                 tmp_ents.resize( num_points );
00207                 result = mbImpl->tag_get_length( tags[i], tsize );
00208                 if( MB_SUCCESS != result ) return result;
00209                 result = mbImpl->tag_get_bytes( tags[i], tsize_bytes );
00210                 if( MB_SUCCESS != result ) return result;
00211             }
00212 
00213             for( int j = 0; j < num_points; j++ )
00214             {
00215                 int lidx;
00216                 if( point_indices )
00217                     lidx = ( *point_indices )[j];
00218                 else
00219                     lidx = j;
00220 
00221                 myLocator->elem_eval()->set_tag_handle( tags[i] );
00222                 myLocator->elem_eval()->set_ent_handle( myLocator->loc_table().vul_rd[lidx] );
00223                 if( !interp_vals ) tmp_ents[j] = targetEnts[lidx];  // Could be performance-sensitive, thus the if test
00224                 result = myLocator->elem_eval()->eval( myLocator->loc_table().vr_rd + 3 * lidx, tmp_dbl );
00225                 tmp_dbl += tsize;
00226                 if( MB_SUCCESS != result ) return result;
00227             }  // for j over points
00228 
00229             if( !interp_vals )
00230             {
00231                 // Set tags on tmp_ents; data is already w/o padding, due to tsize setting above
00232                 result = mbImpl->tag_set_data( tags[i], &tmp_ents[0], tmp_ents.size(), &tmp_vals[0] );
00233                 if( MB_SUCCESS != result ) return result;
00234             }
00235 
00236         }  // for i over methods
00237     }      // if myPcomm
00238 
00239     // Done
00240     return MB_SUCCESS;
00241 }
00242 
00243 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines