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