MOAB: Mesh Oriented datABase  (version 5.2.1)
Coupler.cpp
Go to the documentation of this file.
00001 #include "Coupler.hpp"
00002 #include "moab/ParallelComm.hpp"
00003 #include "moab/AdaptiveKDTree.hpp"
00004 #include "ElemUtil.hpp"
00005 #include "moab/CN.hpp"
00006 #include "moab/gs.hpp"
00007 #include "moab/TupleList.hpp"
00008 #include "moab/Error.hpp"
00009 
00010 /* C++ includes */
00011 #include <cstdio>
00012 #include <cassert>
00013 #include <iostream>
00014 #include <algorithm>
00015 #include <sstream>
00016 
00017 #define ERROR( a )                                           \
00018     {                                                        \
00019         if( MB_SUCCESS != err ) std::cerr << a << std::endl; \
00020     }
00021 #define ERRORR( a, b )                   \
00022     {                                    \
00023         if( MB_SUCCESS != b )            \
00024         {                                \
00025             std::cerr << a << std::endl; \
00026             return b;                    \
00027         }                                \
00028     }
00029 #define ERRORMPI( a, b )                 \
00030     {                                    \
00031         if( MPI_SUCCESS != b )           \
00032         {                                \
00033             std::cerr << a << std::endl; \
00034             return MB_FAILURE;           \
00035         }                                \
00036     }
00037 
00038 #define MASTER_PROC 0
00039 
00040 namespace moab
00041 {
00042 
00043 bool debug = false;
00044 int pack_tuples( TupleList* tl, void** ptr );
00045 void unpack_tuples( void* ptr, TupleList** tlp );
00046 
00047 Coupler::Coupler( Interface* impl, ParallelComm* pc, Range& local_elems, int coupler_id, bool init_tree,
00048                   int max_ent_dim )
00049     : mbImpl( impl ), myPc( pc ), myId( coupler_id ), numIts( 3 ), max_dim( max_ent_dim ), _ntot( 0 ),
00050       spherical( false )
00051 {
00052     assert( NULL != impl && ( pc || !local_elems.empty() ) );
00053 
00054     // Keep track of the local points, at least for now
00055     myRange = local_elems;
00056     myTree  = NULL;
00057 
00058     // Now initialize the tree
00059     if( init_tree ) initialize_tree();
00060 
00061     // Initialize tuple lists to indicate not initialized
00062     mappedPts       = NULL;
00063     targetPts       = NULL;
00064     _spectralSource = _spectralTarget = NULL;
00065 }
00066 
00067 /* Destructor
00068  */
00069 Coupler::~Coupler()
00070 {
00071     // This will clear the cache
00072     delete(moab::Element::SpectralHex*)_spectralSource;
00073     delete(moab::Element::SpectralHex*)_spectralTarget;
00074     delete myTree;
00075     delete targetPts;
00076     delete mappedPts;
00077 }
00078 
00079 ErrorCode Coupler::initialize_tree()
00080 {
00081     Range local_ents;
00082 
00083     // Get entities on the local part
00084     ErrorCode result = MB_SUCCESS;
00085     if( myPc )
00086     {
00087         result = myPc->get_part_entities( local_ents, max_dim );
00088         if( local_ents.empty() )
00089         {
00090             max_dim--;
00091             result = myPc->get_part_entities( local_ents, max_dim );  // go one dimension lower
00092             // right now, this is used for spherical meshes only
00093         }
00094     }
00095     else
00096         local_ents = myRange;
00097     if( MB_SUCCESS != result || local_ents.empty() )
00098     {
00099         std::cout << "Problems getting source entities" << std::endl;
00100         return result;
00101     }
00102 
00103     // Build the tree for local processor
00104     int max_per_leaf = 6;
00105     for( int i = 0; i < numIts; i++ )
00106     {
00107         std::ostringstream str;
00108         str << "PLANE_SET=0;"
00109             << "MAX_PER_LEAF=" << max_per_leaf << ";";
00110         if( spherical && !local_ents.empty() )
00111         {
00112             // get coordinates of one vertex, and use for radius computation
00113             EntityHandle elem = local_ents[0];
00114             const EntityHandle* conn;
00115             int numn = 0;
00116             mbImpl->get_connectivity( elem, conn, numn );
00117             CartVect pos0;
00118             mbImpl->get_coords( conn, 1, &( pos0[0] ) );
00119             double radius = pos0.length();
00120             str << "SPHERICAL=true;RADIUS=" << radius << ";";
00121         }
00122         FileOptions opts( str.str().c_str() );
00123         myTree = new AdaptiveKDTree( mbImpl );
00124         result = myTree->build_tree( local_ents, &localRoot, &opts );
00125         if( MB_SUCCESS != result )
00126         {
00127             std::cout << "Problems building tree";
00128             if( numIts != i )
00129             {
00130                 delete myTree;
00131                 max_per_leaf *= 2;
00132                 std::cout << "; increasing elements/leaf to " << max_per_leaf << std::endl;
00133             }
00134             else
00135             {
00136                 std::cout << "; exiting" << std::endl;
00137                 return result;
00138             }
00139         }
00140         else
00141             break;  // Get out of tree building
00142     }
00143 
00144     // Get the bounding box for local tree
00145     if( myPc )
00146         allBoxes.resize( 6 * myPc->proc_config().proc_size() );
00147     else
00148         allBoxes.resize( 6 );
00149 
00150     unsigned int my_rank = ( myPc ? myPc->proc_config().proc_rank() : 0 );
00151     BoundBox box;
00152     result = myTree->get_bounding_box( box, &localRoot );
00153     if( MB_SUCCESS != result ) return result;
00154     box.bMin.get( &allBoxes[6 * my_rank] );
00155     box.bMax.get( &allBoxes[6 * my_rank + 3] );
00156 
00157     // Now communicate to get all boxes
00158     if( myPc )
00159     {
00160         int mpi_err;
00161 #if( MPI_VERSION >= 2 )
00162         // Use "in place" option
00163         mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 6, MPI_DOUBLE,
00164                                  myPc->proc_config().proc_comm() );
00165 #else
00166         {
00167             std::vector< double > allBoxes_tmp( 6 * myPc->proc_config().proc_size() );
00168             mpi_err  = MPI_Allgather( &allBoxes[6 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 6, MPI_DOUBLE,
00169                                      myPc->proc_config().proc_comm() );
00170             allBoxes = allBoxes_tmp;
00171         }
00172 #endif
00173         if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
00174     }
00175 
00176     /*  std::ostringstream blah;
00177     for(int i = 0; i < allBoxes.size(); i++)
00178     blah << allBoxes[i] << " ";
00179     std::cout << blah.str() << "\n";*/
00180 
00181 #ifdef VERBOSE
00182     double min[3] = { 0, 0, 0 }, max[3] = { 0, 0, 0 };
00183     unsigned int dep;
00184     myTree->get_info( localRoot, min, max, dep );
00185     std::cout << "Proc " << my_rank << ": box min/max, tree depth = (" << min[0] << "," << min[1] << "," << min[2]
00186               << "), (" << max[0] << "," << max[1] << "," << max[2] << "), " << dep << std::endl;
00187 #endif
00188 
00189     return result;
00190 }
00191 
00192 ErrorCode Coupler::initialize_spectral_elements( EntityHandle rootSource, EntityHandle rootTarget, bool& specSou,
00193                                                  bool& specTar )
00194 {
00195     /*void * _spectralSource;
00196       void * _spectralTarget;*/
00197 
00198     moab::Range spectral_sets;
00199     moab::Tag sem_tag;
00200     int sem_dims[3];
00201     ErrorCode rval = mbImpl->tag_get_handle( "SEM_DIMS", 3, moab::MB_TYPE_INTEGER, sem_tag );
00202     if( moab::MB_SUCCESS != rval )
00203     {
00204         std::cout << "Can't find tag, no spectral set\n";
00205         return MB_SUCCESS;  // Nothing to do, no spectral elements
00206     }
00207     rval = mbImpl->get_entities_by_type_and_tag( rootSource, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets );
00208     if( moab::MB_SUCCESS != rval || spectral_sets.empty() )
00209         std::cout << "Can't get sem set on source\n";
00210     else
00211     {
00212         moab::EntityHandle firstSemSet = spectral_sets[0];
00213         rval                           = mbImpl->tag_get_data( sem_tag, &firstSemSet, 1, (void*)sem_dims );
00214         if( moab::MB_SUCCESS != rval ) return MB_FAILURE;
00215 
00216         if( sem_dims[0] != sem_dims[1] || sem_dims[0] != sem_dims[2] )
00217         {
00218             std::cout << " dimensions are different. bail out\n";
00219             return MB_FAILURE;
00220         }
00221         // Repeat for target sets
00222         spectral_sets.empty();
00223         // Now initialize a source spectral element !
00224         _spectralSource = new moab::Element::SpectralHex( sem_dims[0] );
00225         specSou         = true;
00226     }
00227     // Repeat for target source
00228     rval = mbImpl->get_entities_by_type_and_tag( rootTarget, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets );
00229     if( moab::MB_SUCCESS != rval || spectral_sets.empty() )
00230         std::cout << "Can't get sem set on target\n";
00231     else
00232     {
00233         moab::EntityHandle firstSemSet = spectral_sets[0];
00234         rval                           = mbImpl->tag_get_data( sem_tag, &firstSemSet, 1, (void*)sem_dims );
00235         if( moab::MB_SUCCESS != rval ) return MB_FAILURE;
00236 
00237         if( sem_dims[0] != sem_dims[1] || sem_dims[0] != sem_dims[2] )
00238         {
00239             std::cout << " dimensions are different. bail out\n";
00240             return MB_FAILURE;
00241         }
00242         // Repeat for target sets
00243         spectral_sets.empty();
00244         // Now initialize a target spectral element !
00245         _spectralTarget = new moab::Element::SpectralHex( sem_dims[0] );
00246         specTar         = true;
00247     }
00248 
00249     _ntot = sem_dims[0] * sem_dims[1] * sem_dims[2];
00250     rval  = mbImpl->tag_get_handle( "SEM_X", _ntot, moab::MB_TYPE_DOUBLE, _xm1Tag );
00251     if( moab::MB_SUCCESS != rval )
00252     {
00253         std::cout << "Can't get xm1tag \n";
00254         return MB_FAILURE;
00255     }
00256     rval = mbImpl->tag_get_handle( "SEM_Y", _ntot, moab::MB_TYPE_DOUBLE, _ym1Tag );
00257     if( moab::MB_SUCCESS != rval )
00258     {
00259         std::cout << "Can't get ym1tag \n";
00260         return MB_FAILURE;
00261     }
00262     rval = mbImpl->tag_get_handle( "SEM_Z", _ntot, moab::MB_TYPE_DOUBLE, _zm1Tag );
00263     if( moab::MB_SUCCESS != rval )
00264     {
00265         std::cout << "Can't get zm1tag \n";
00266         return MB_FAILURE;
00267     }
00268 
00269     return MB_SUCCESS;
00270 }
00271 
00272 ErrorCode Coupler::locate_points( Range& targ_ents, double rel_eps, double abs_eps, TupleList* tl, bool store_local )
00273 {
00274     // Get locations
00275     std::vector< double > locs( 3 * targ_ents.size() );
00276     Range verts    = targ_ents.subset_by_type( MBVERTEX );
00277     ErrorCode rval = mbImpl->get_coords( verts, &locs[0] );
00278     if( MB_SUCCESS != rval ) return rval;
00279     // Now get other ents; reuse verts
00280     unsigned int num_verts = verts.size();
00281     verts                  = subtract( targ_ents, verts );
00282     // Compute centroids
00283     std::vector< EntityHandle > dum_conn( CN::MAX_NODES_PER_ELEMENT );
00284     std::vector< double > dum_pos( CN::MAX_NODES_PER_ELEMENT );
00285     const EntityHandle* conn;
00286     int num_conn;
00287     double* coords = &locs[num_verts];
00288     // Do this here instead of a function to allow reuse of dum_pos and dum_conn
00289     for( Range::const_iterator rit = verts.begin(); rit != verts.end(); ++rit )
00290     {
00291         rval = mbImpl->get_connectivity( *rit, conn, num_conn, false, &dum_conn );
00292         if( MB_SUCCESS != rval ) return rval;
00293         rval = mbImpl->get_coords( conn, num_conn, &dum_pos[0] );
00294         if( MB_SUCCESS != rval ) return rval;
00295         coords[0] = coords[1] = coords[2] = 0.0;
00296         for( int i = 0; i < num_conn; i++ )
00297         {
00298             coords[0] += dum_pos[3 * i];
00299             coords[1] += dum_pos[3 * i + 1];
00300             coords[2] += dum_pos[3 * i + 2];
00301         }
00302         coords[0] /= num_conn;
00303         coords[1] /= num_conn;
00304         coords[2] /= num_conn;
00305         coords += 3;
00306     }
00307 
00308     if( store_local ) targetEnts = targ_ents;
00309 
00310     return locate_points( &locs[0], targ_ents.size(), rel_eps, abs_eps, tl, store_local );
00311 }
00312 
00313 ErrorCode Coupler::locate_points( double* xyz, unsigned int num_points, double rel_eps, double abs_eps, TupleList* tl,
00314                                   bool store_local )
00315 {
00316     assert( tl || store_local );
00317 
00318     // target_pts: TL(to_proc, tgt_index, x, y, z): tuples sent to source mesh procs representing
00319     // pts to be located source_pts: TL(from_proc, tgt_index, src_index): results of source mesh
00320     // proc point location, ready to send
00321     //             back to tgt procs; src_index of -1 indicates point not located (arguably not
00322     //             useful...)
00323     TupleList target_pts;
00324     target_pts.initialize( 2, 0, 0, 3, num_points );
00325     target_pts.enableWriteAccess();
00326 
00327     TupleList source_pts;
00328     mappedPts = new TupleList( 0, 0, 1, 3, target_pts.get_max() );
00329     mappedPts->enableWriteAccess();
00330 
00331     source_pts.initialize( 3, 0, 0, 0, target_pts.get_max() );
00332     source_pts.enableWriteAccess();
00333 
00334     mappedPts->set_n( 0 );
00335     source_pts.set_n( 0 );
00336     ErrorCode result;
00337 
00338     unsigned int my_rank = ( myPc ? myPc->proc_config().proc_rank() : 0 );
00339     bool point_located;
00340 
00341     // For each point, find box(es) containing the point,
00342     // appending results to tuple_list;
00343     // keep local points separately, in local_pts, which has pairs
00344     // of <local_index, mapped_index>, where mapped_index is the index
00345     // into the mappedPts tuple list
00346     for( unsigned int i = 0; i < 3 * num_points; i += 3 )
00347     {
00348 
00349         std::vector< int > procs_to_send_to;
00350         for( unsigned int j = 0; j < ( myPc ? myPc->proc_config().proc_size() : 0 ); j++ )
00351         {
00352             // Test if point is in proc's box
00353             if( ( allBoxes[6 * j] <= xyz[i] + abs_eps ) && ( xyz[i] <= allBoxes[6 * j + 3] + abs_eps ) &&
00354                 ( allBoxes[6 * j + 1] <= xyz[i + 1] + abs_eps ) && ( xyz[i + 1] <= allBoxes[6 * j + 4] + abs_eps ) &&
00355                 ( allBoxes[6 * j + 2] <= xyz[i + 2] + abs_eps ) && ( xyz[i + 2] <= allBoxes[6 * j + 5] + abs_eps ) )
00356             {
00357                 // If in this proc's box, will send to proc to test further
00358                 procs_to_send_to.push_back( j );
00359             }
00360         }
00361         if( procs_to_send_to.empty() )
00362         {
00363 #ifdef VERBOSE
00364             std::cout << " point index " << i / 3 << ": " << xyz[i] << " " << xyz[i + 1] << " " << xyz[i + 2]
00365                       << " not found in any box\n";
00366 #endif
00367             // try to find the closest box, and put it in that box, anyway
00368             double min_dist = 1.e+20;
00369             int index       = -1;
00370             for( unsigned int j = 0; j < ( myPc ? myPc->proc_config().proc_size() : 0 ); j++ )
00371             {
00372                 BoundBox box( &allBoxes[6 * j] );           // form back the box
00373                 double distance = box.distance( &xyz[i] );  // will compute the distance in 3d, from the box
00374                 if( distance < min_dist )
00375                 {
00376                     index    = j;
00377                     min_dist = distance;
00378                 }
00379             }
00380             if( index == -1 )
00381             {
00382                 // need to abort early, nothing we can do
00383                 assert( "cannot locate any box for some points" );
00384                 // need a better exit strategy
00385             }
00386 #ifdef VERBOSE
00387             std::cout << " point index " << i / 3 << " added to box for proc j:" << index << "\n";
00388 #endif
00389             procs_to_send_to.push_back( index );  // will send to just one proc, that has the closest box
00390         }
00391         // we finally decided to populate the tuple list for a list of processors
00392         for( size_t k = 0; k < procs_to_send_to.size(); k++ )
00393         {
00394             unsigned int j = procs_to_send_to[k];
00395             // Check size of tuple list, grow if we're at max
00396             if( target_pts.get_n() == target_pts.get_max() )
00397                 target_pts.resize( std::max( 10.0, 1.5 * target_pts.get_max() ) );
00398 
00399             target_pts.vi_wr[2 * target_pts.get_n()]     = j;
00400             target_pts.vi_wr[2 * target_pts.get_n() + 1] = i / 3;
00401 
00402             target_pts.vr_wr[3 * target_pts.get_n()]     = xyz[i];
00403             target_pts.vr_wr[3 * target_pts.get_n() + 1] = xyz[i + 1];
00404             target_pts.vr_wr[3 * target_pts.get_n() + 2] = xyz[i + 2];
00405             target_pts.inc_n();
00406         }
00407 
00408     }  // end for (unsigned int i = 0; ..
00409 
00410     int num_to_me = 0;
00411     for( unsigned int i = 0; i < target_pts.get_n(); i++ )
00412         if( target_pts.vi_rd[2 * i] == (int)my_rank ) num_to_me++;
00413 #ifdef VERBOSE
00414     printf( "rank: %u local points: %u, nb sent target pts: %u mappedPts: %u num to me: %d \n", my_rank, num_points,
00415             target_pts.get_n(), mappedPts->get_n(), num_to_me );
00416 #endif
00417     // Perform scatter/gather, to gather points to source mesh procs
00418     if( myPc )
00419     {
00420         ( myPc->proc_config().crystal_router() )->gs_transfer( 1, target_pts, 0 );
00421 
00422         num_to_me = 0;
00423         for( unsigned int i = 0; i < target_pts.get_n(); i++ )
00424         {
00425             if( target_pts.vi_rd[2 * i] == (int)my_rank ) num_to_me++;
00426         }
00427 #ifdef VERBOSE
00428         printf( "rank: %u after first gs nb received_pts: %u; num_from_me = %d\n", my_rank, target_pts.get_n(),
00429                 num_to_me );
00430 #endif
00431         // After scatter/gather:
00432         // target_pts.set_n(# points local proc has to map);
00433         // target_pts.vi_wr[2*i] = proc sending point i
00434         // target_pts.vi_wr[2*i + 1] = index of point i on sending proc
00435         // target_pts.vr_wr[3*i..3*i + 2] = xyz of point i
00436         //
00437         // Mapping builds the tuple list:
00438         // source_pts.set_n(target_pts.get_n())
00439         // source_pts.vi_wr[3*i] = target_pts.vi_wr[2*i] = sending proc
00440         // source_pts.vi_wr[3*i + 1] = index of point i on sending proc
00441         // source_pts.vi_wr[3*i + 2] = index of mapped point (-1 if not mapped)
00442         //
00443         // Also, mapping builds local tuple_list mappedPts:
00444         // mappedPts->set_n( # mapped points );
00445         // mappedPts->vul_wr[i] = local handle of mapped entity
00446         // mappedPts->vr_wr[3*i..3*i + 2] = natural coordinates in mapped entity
00447 
00448         // Test target points against my elements
00449         for( unsigned i = 0; i < target_pts.get_n(); i++ )
00450         {
00451             result = test_local_box( target_pts.vr_wr + 3 * i, target_pts.vi_rd[2 * i], target_pts.vi_rd[2 * i + 1], i,
00452                                      point_located, rel_eps, abs_eps, &source_pts );
00453             if( MB_SUCCESS != result ) return result;
00454         }
00455 
00456         // No longer need target_pts
00457         target_pts.reset();
00458 #ifdef VERBOSE
00459         printf( "rank: %u nb sent source pts: %u, mappedPts now: %u\n", my_rank, source_pts.get_n(),
00460                 mappedPts->get_n() );
00461 #endif
00462         // Send target points back to target procs
00463         ( myPc->proc_config().crystal_router() )->gs_transfer( 1, source_pts, 0 );
00464 
00465 #ifdef VERBOSE
00466         printf( "rank: %u nb received source pts: %u\n", my_rank, source_pts.get_n() );
00467 #endif
00468     }
00469 
00470     // Store proc/index tuples in targetPts, and/or pass back to application;
00471     // the tuple this gets stored to looks like:
00472     // tl.set_n(# mapped points);
00473     // tl.vi_wr[3*i] = remote proc mapping point
00474     // tl.vi_wr[3*i + 1] = local index of mapped point
00475     // tl.vi_wr[3*i + 2] = remote index of mapped point
00476     //
00477     // Local index is mapped into either myRange, holding the handles of
00478     // local mapped entities, or myXyz, holding locations of mapped pts
00479 
00480     // Store information about located points
00481     TupleList* tl_tmp;
00482     if( !store_local )
00483         tl_tmp = tl;
00484     else
00485     {
00486         targetPts = new TupleList();
00487         tl_tmp    = targetPts;
00488     }
00489 
00490     tl_tmp->initialize( 3, 0, 0, 0, num_points );
00491     tl_tmp->set_n( num_points );  // Automatically sets tl to write_enabled
00492     // Initialize so we know afterwards how many pts weren't located
00493     std::fill( tl_tmp->vi_wr, tl_tmp->vi_wr + 3 * num_points, -1 );
00494 
00495     unsigned int local_pts = 0;
00496     for( unsigned int i = 0; i < source_pts.get_n(); i++ )
00497     {
00498         if( -1 != source_pts.vi_rd[3 * i + 2] )
00499         {  // Why bother sending message saying "i don't have the point" if it gets discarded?
00500             int tgt_index = 3 * source_pts.vi_rd[3 * i + 1];
00501             // Prefer always entities that are local, from the source_pts
00502             // if a local entity was already found to contain the target point, skip
00503             // tl_tmp->vi_wr[tgt_index] is -1 initially, but it could already be set with
00504             // a remote processor
00505             if( tl_tmp->vi_wr[tgt_index] != (int)my_rank )
00506             {
00507                 tl_tmp->vi_wr[tgt_index]     = source_pts.vi_rd[3 * i];
00508                 tl_tmp->vi_wr[tgt_index + 1] = source_pts.vi_rd[3 * i + 1];
00509                 tl_tmp->vi_wr[tgt_index + 2] = source_pts.vi_rd[3 * i + 2];
00510             }
00511         }
00512     }
00513 
00514     // Count missing points
00515     unsigned int missing_pts = 0;
00516     for( unsigned int i = 0; i < num_points; i++ )
00517     {
00518         if( tl_tmp->vi_rd[3 * i + 1] == -1 )
00519         {
00520             missing_pts++;
00521 #ifdef VERBOSE
00522             printf( "missing point at index i:  %d -> %15.10f %15.10f %15.10f\n", i, xyz[3 * i], xyz[3 * i + 1],
00523                     xyz[3 * i + 2] );
00524 #endif
00525         }
00526         else if( tl_tmp->vi_rd[3 * i] == (int)my_rank )
00527             local_pts++;
00528     }
00529 #ifdef VERBOSE
00530     printf( "rank: %u point location: wanted %u got %u locally, %u remote, missing %u\n", my_rank, num_points,
00531             local_pts, num_points - missing_pts - local_pts, missing_pts );
00532 #endif
00533     assert( 0 == missing_pts );  // Will likely break on curved geometries
00534 
00535     // No longer need source_pts
00536     source_pts.reset();
00537 
00538     // Copy into tl if passed in and storing locally
00539     if( tl && store_local )
00540     {
00541         tl->initialize( 3, 0, 0, 0, num_points );
00542         tl->enableWriteAccess();
00543         memcpy( tl->vi_wr, tl_tmp->vi_rd, 3 * tl_tmp->get_n() * sizeof( int ) );
00544         tl->set_n( tl_tmp->get_n() );
00545         tl->disableWriteAccess();
00546     }
00547 
00548     tl_tmp->disableWriteAccess();
00549 
00550     // Done
00551     return MB_SUCCESS;
00552 }
00553 
00554 ErrorCode Coupler::test_local_box( double* xyz, int from_proc, int remote_index, int /*index*/, bool& point_located,
00555                                    double rel_eps, double abs_eps, TupleList* tl )
00556 {
00557     std::vector< EntityHandle > entities;
00558     std::vector< CartVect > nat_coords;
00559     bool canWrite = false;
00560     if( tl )
00561     {
00562         canWrite = tl->get_writeEnabled();
00563         if( !canWrite )
00564         {
00565             tl->enableWriteAccess();
00566             canWrite = true;
00567         }
00568     }
00569 
00570     if( rel_eps && !abs_eps )
00571     {
00572         // Relative epsilon given, translate to absolute epsilon using box dimensions
00573         BoundBox box;
00574         myTree->get_bounding_box( box, &localRoot );
00575         abs_eps = rel_eps * box.diagonal_length();
00576     }
00577 
00578     ErrorCode result = nat_param( xyz, entities, nat_coords, abs_eps );
00579     if( MB_SUCCESS != result ) return result;
00580 
00581     // If we didn't find any ents and we're looking locally, nothing more to do
00582     if( entities.empty() )
00583     {
00584         if( tl->get_n() == tl->get_max() ) tl->resize( std::max( 10.0, 1.5 * tl->get_max() ) );
00585 
00586         tl->vi_wr[3 * tl->get_n()]     = from_proc;
00587         tl->vi_wr[3 * tl->get_n() + 1] = remote_index;
00588         tl->vi_wr[3 * tl->get_n() + 2] = -1;
00589         tl->inc_n();
00590 
00591         point_located = false;
00592         return MB_SUCCESS;
00593     }
00594 
00595     // Grow if we know we'll exceed size
00596     if( mappedPts->get_n() + entities.size() >= mappedPts->get_max() )
00597         mappedPts->resize( std::max( 10.0, 1.5 * mappedPts->get_max() ) );
00598 
00599     std::vector< EntityHandle >::iterator eit = entities.begin();
00600     std::vector< CartVect >::iterator ncit    = nat_coords.begin();
00601 
00602     mappedPts->enableWriteAccess();
00603     for( ; eit != entities.end(); ++eit, ++ncit )
00604     {
00605         // Store in tuple mappedPts
00606         mappedPts->vr_wr[3 * mappedPts->get_n()]     = ( *ncit )[0];
00607         mappedPts->vr_wr[3 * mappedPts->get_n() + 1] = ( *ncit )[1];
00608         mappedPts->vr_wr[3 * mappedPts->get_n() + 2] = ( *ncit )[2];
00609         mappedPts->vul_wr[mappedPts->get_n()]        = *eit;
00610         mappedPts->inc_n();
00611 
00612         // Also store local point, mapped point indices
00613         if( tl->get_n() == tl->get_max() ) tl->resize( std::max( 10.0, 1.5 * tl->get_max() ) );
00614 
00615         // Store in tuple source_pts
00616         tl->vi_wr[3 * tl->get_n()]     = from_proc;
00617         tl->vi_wr[3 * tl->get_n() + 1] = remote_index;
00618         tl->vi_wr[3 * tl->get_n() + 2] = mappedPts->get_n() - 1;
00619         tl->inc_n();
00620     }
00621 
00622     point_located = true;
00623 
00624     if( tl && !canWrite ) tl->disableWriteAccess();
00625 
00626     return MB_SUCCESS;
00627 }
00628 
00629 ErrorCode Coupler::interpolate( Coupler::Method method, const std::string& interp_tag, double* interp_vals,
00630                                 TupleList* tl, bool normalize )
00631 {
00632     Tag tag;
00633     ErrorCode result;
00634     if( _spectralSource )
00635     {
00636         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 << "\"" );
00637     }
00638     else
00639     {
00640         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 << "\"" );
00641     }
00642 
00643     return interpolate( method, tag, interp_vals, tl, normalize );
00644 }
00645 
00646 ErrorCode Coupler::interpolate( Coupler::Method* methods, Tag* tags, int* points_per_method, int num_methods,
00647                                 double* interp_vals, TupleList* tl, bool /* normalize */ )
00648 {
00649     // if (!((LINEAR_FE == method) || (CONSTANT == method)))
00650     // return MB_FAILURE;
00651 
00652     // remote pts first
00653     TupleList* tl_tmp = ( tl ? tl : targetPts );
00654 
00655     ErrorCode result = MB_SUCCESS;
00656 
00657     unsigned int pts_total = 0;
00658     for( int i = 0; i < num_methods; i++ )
00659         pts_total += points_per_method[i];
00660 
00661     // If tl was passed in non-NULL, just have those points, otherwise have targetPts plus
00662     // locally mapped pts
00663     if( pts_total != tl_tmp->get_n() ) return MB_FAILURE;
00664 
00665     TupleList tinterp;
00666     tinterp.initialize( 5, 0, 0, 1, tl_tmp->get_n() );
00667     int t = 0;
00668     tinterp.enableWriteAccess();
00669     for( int i = 0; i < num_methods; i++ )
00670     {
00671         for( int j = 0; j < points_per_method[i]; j++ )
00672         {
00673             tinterp.vi_wr[5 * t]     = tl_tmp->vi_rd[3 * t];
00674             tinterp.vi_wr[5 * t + 1] = tl_tmp->vi_rd[3 * t + 1];
00675             tinterp.vi_wr[5 * t + 2] = tl_tmp->vi_rd[3 * t + 2];
00676             tinterp.vi_wr[5 * t + 3] = methods[i];
00677             tinterp.vi_wr[5 * t + 4] = i;
00678             tinterp.vr_wr[t]         = 0.0;
00679             tinterp.inc_n();
00680             t++;
00681         }
00682     }
00683 
00684     // Scatter/gather interpolation points
00685     if( myPc )
00686     {
00687         ( myPc->proc_config().crystal_router() )->gs_transfer( 1, tinterp, 0 );
00688 
00689         // Perform interpolation on local source mesh; put results into
00690         // tinterp.vr_wr[i]
00691 
00692         mappedPts->enableWriteAccess();
00693         for( unsigned int i = 0; i < tinterp.get_n(); i++ )
00694         {
00695             int mindex    = tinterp.vi_rd[5 * i + 2];
00696             Method method = (Method)tinterp.vi_rd[5 * i + 3];
00697             Tag tag       = tags[tinterp.vi_rd[5 * i + 4]];
00698 
00699             result = MB_FAILURE;
00700             if( LINEAR_FE == method || QUADRATIC_FE == method || SPHERICAL == method )
00701             {
00702                 result = interp_field( mappedPts->vul_rd[mindex], CartVect( mappedPts->vr_wr + 3 * mindex ), tag,
00703                                        tinterp.vr_wr[i] );
00704             }
00705             else if( CONSTANT == method )
00706             {
00707                 result = constant_interp( mappedPts->vul_rd[mindex], tag, tinterp.vr_wr[i] );
00708             }
00709 
00710             if( MB_SUCCESS != result ) return result;
00711         }
00712 
00713         // Scatter/gather interpolation data
00714         myPc->proc_config().crystal_router()->gs_transfer( 1, tinterp, 0 );
00715     }
00716 
00717     // Copy the interpolated field as a unit
00718     for( unsigned int i = 0; i < tinterp.get_n(); i++ )
00719         interp_vals[tinterp.vi_rd[5 * i + 1]] = tinterp.vr_rd[i];
00720 
00721     // Done
00722     return MB_SUCCESS;
00723 }
00724 
00725 ErrorCode Coupler::nat_param( double xyz[3], std::vector< EntityHandle >& entities, std::vector< CartVect >& nat_coords,
00726                               double epsilon )
00727 {
00728     if( !myTree ) return MB_FAILURE;
00729 
00730     AdaptiveKDTreeIter treeiter;
00731     ErrorCode result = myTree->get_tree_iterator( localRoot, treeiter );
00732     if( MB_SUCCESS != result )
00733     {
00734         std::cout << "Problems getting iterator" << std::endl;
00735         return result;
00736     }
00737 
00738     EntityHandle closest_leaf;
00739     if( epsilon )
00740     {
00741         std::vector< double > dists;
00742         std::vector< EntityHandle > leaves;
00743 
00744         // Two tolerances
00745         result = myTree->distance_search( xyz, epsilon, leaves,
00746                                           /*iter_tol*/ epsilon,
00747                                           /*inside_tol*/ 10 * epsilon, &dists, NULL, &localRoot );
00748         if( leaves.empty() )
00749             // Not found returns success here, with empty list, just like case with no epsilon
00750             return MB_SUCCESS;
00751         // Get closest leaf
00752         double min_dist                           = *dists.begin();
00753         closest_leaf                              = *leaves.begin();
00754         std::vector< EntityHandle >::iterator vit = leaves.begin() + 1;
00755         std::vector< double >::iterator dit       = dists.begin() + 1;
00756         for( ; vit != leaves.end() && min_dist; ++vit, ++dit )
00757         {
00758             if( *dit < min_dist )
00759             {
00760                 min_dist     = *dit;
00761                 closest_leaf = *vit;
00762             }
00763         }
00764     }
00765     else
00766     {
00767         result = myTree->point_search( xyz, treeiter, 1.0e-10, 1.0e-6, NULL, &localRoot );
00768         if( MB_ENTITY_NOT_FOUND == result )  // Point is outside of myTree's bounding box
00769             return MB_SUCCESS;
00770         else if( MB_SUCCESS != result )
00771         {
00772             std::cout << "Problems getting leaf \n";
00773             return result;
00774         }
00775         closest_leaf = treeiter.handle();
00776     }
00777 
00778     // Find natural coordinates of point in element(s) in that leaf
00779     CartVect tmp_nat_coords;
00780     Range range_leaf;
00781     result = mbImpl->get_entities_by_dimension( closest_leaf, max_dim, range_leaf, false );
00782     if( MB_SUCCESS != result ) std::cout << "Problem getting leaf in a range" << std::endl;
00783 
00784     // Loop over the range_leaf
00785     for( Range::iterator iter = range_leaf.begin(); iter != range_leaf.end(); ++iter )
00786     {
00787         // Test to find out in which entity the point is
00788         // Get the EntityType and create the appropriate Element::Map subtype
00789         // If spectral, do not need coordinates, just the GL points
00790         EntityType etype = mbImpl->type_from_handle( *iter );
00791         if( NULL != this->_spectralSource && MBHEX == etype )
00792         {
00793             EntityHandle eh = *iter;
00794             const double* xval;
00795             const double* yval;
00796             const double* zval;
00797             ErrorCode rval = mbImpl->tag_get_by_ptr( _xm1Tag, &eh, 1, (const void**)&xval );
00798             if( moab::MB_SUCCESS != rval )
00799             {
00800                 std::cout << "Can't get xm1 values \n";
00801                 return MB_FAILURE;
00802             }
00803             rval = mbImpl->tag_get_by_ptr( _ym1Tag, &eh, 1, (const void**)&yval );
00804             if( moab::MB_SUCCESS != rval )
00805             {
00806                 std::cout << "Can't get ym1 values \n";
00807                 return MB_FAILURE;
00808             }
00809             rval = mbImpl->tag_get_by_ptr( _zm1Tag, &eh, 1, (const void**)&zval );
00810             if( moab::MB_SUCCESS != rval )
00811             {
00812                 std::cout << "Can't get zm1 values \n";
00813                 return MB_FAILURE;
00814             }
00815             Element::SpectralHex* spcHex = (Element::SpectralHex*)_spectralSource;
00816 
00817             spcHex->set_gl_points( (double*)xval, (double*)yval, (double*)zval );
00818             try
00819             {
00820                 tmp_nat_coords = spcHex->ievaluate( CartVect( xyz ), epsilon );  // introduce
00821                 bool inside    = spcHex->inside_nat_space( CartVect( tmp_nat_coords ), epsilon );
00822                 if( !inside )
00823                 {
00824 #ifdef VERBOSE
00825                     std::cout << "point " << xyz[0] << " " << xyz[1] << " " << xyz[2]
00826                               << " is not converging inside hex " << mbImpl->id_from_handle( eh ) << "\n";
00827 #endif
00828                     continue;  // It is possible that the point is outside, so it will not converge
00829                 }
00830             }
00831             catch( Element::Map::EvaluationError& )
00832             {
00833                 continue;
00834             }
00835         }
00836         else
00837         {
00838             const EntityHandle* connect;
00839             int num_connect;
00840 
00841             // Get connectivity
00842             result = mbImpl->get_connectivity( *iter, connect, num_connect, true );
00843             if( MB_SUCCESS != result ) return result;
00844 
00845             // Get coordinates of the vertices
00846             std::vector< CartVect > coords_vert( num_connect );
00847             result = mbImpl->get_coords( connect, num_connect, &( coords_vert[0][0] ) );
00848             if( MB_SUCCESS != result )
00849             {
00850                 std::cout << "Problems getting coordinates of vertices\n";
00851                 return result;
00852             }
00853             CartVect pos( xyz );
00854             if( MBHEX == etype )
00855             {
00856                 if( 8 == num_connect )
00857                 {
00858                     Element::LinearHex hexmap( coords_vert );
00859                     if( !hexmap.inside_box( pos, epsilon ) ) continue;
00860                     try
00861                     {
00862                         tmp_nat_coords = hexmap.ievaluate( pos, epsilon );
00863                         bool inside    = hexmap.inside_nat_space( tmp_nat_coords, epsilon );
00864                         if( !inside ) continue;
00865                     }
00866                     catch( Element::Map::EvaluationError& )
00867                     {
00868                         continue;
00869                     }
00870                 }
00871                 else if( 27 == num_connect )
00872                 {
00873                     Element::QuadraticHex hexmap( coords_vert );
00874                     if( !hexmap.inside_box( pos, epsilon ) ) continue;
00875                     try
00876                     {
00877                         tmp_nat_coords = hexmap.ievaluate( pos, epsilon );
00878                         bool inside    = hexmap.inside_nat_space( tmp_nat_coords, epsilon );
00879                         if( !inside ) continue;
00880                     }
00881                     catch( Element::Map::EvaluationError& )
00882                     {
00883                         continue;
00884                     }
00885                 }
00886                 else  // TODO this case not treated yet, no interpolation
00887                     continue;
00888             }
00889             else if( MBTET == etype )
00890             {
00891                 Element::LinearTet tetmap( coords_vert );
00892                 // This is just a linear solve; unless degenerate, will not except
00893                 tmp_nat_coords = tetmap.ievaluate( pos );
00894                 bool inside    = tetmap.inside_nat_space( tmp_nat_coords, epsilon );
00895                 if( !inside ) continue;
00896             }
00897             else if( MBQUAD == etype && spherical )
00898             {
00899                 Element::SphericalQuad sphermap( coords_vert );
00900                 /* skip box test, because it can filter out good elements with high curvature
00901                  * if (!sphermap.inside_box(pos, epsilon))
00902                   continue;*/
00903                 try
00904                 {
00905                     tmp_nat_coords = sphermap.ievaluate( pos, epsilon );
00906                     bool inside    = sphermap.inside_nat_space( tmp_nat_coords, epsilon );
00907                     if( !inside ) continue;
00908                 }
00909                 catch( Element::Map::EvaluationError& )
00910                 {
00911                     continue;
00912                 }
00913             }
00914             else if( MBTRI == etype && spherical )
00915             {
00916                 Element::SphericalTri sphermap( coords_vert );
00917                 /* skip box test, because it can filter out good elements with high curvature
00918                  * if (!sphermap.inside_box(pos, epsilon))
00919                     continue;*/
00920                 try
00921                 {
00922                     tmp_nat_coords = sphermap.ievaluate( pos, epsilon );
00923                     bool inside    = sphermap.inside_nat_space( tmp_nat_coords, epsilon );
00924                     if( !inside ) continue;
00925                 }
00926                 catch( Element::Map::EvaluationError& )
00927                 {
00928                     continue;
00929                 }
00930             }
00931 
00932             else if( MBQUAD == etype )
00933             {
00934                 Element::LinearQuad quadmap( coords_vert );
00935                 if( !quadmap.inside_box( pos, epsilon ) ) continue;
00936                 try
00937                 {
00938                     tmp_nat_coords = quadmap.ievaluate( pos, epsilon );
00939                     bool inside    = quadmap.inside_nat_space( tmp_nat_coords, epsilon );
00940                     if( !inside ) continue;
00941                 }
00942                 catch( Element::Map::EvaluationError& )
00943                 {
00944                     continue;
00945                 }
00946                 if( !quadmap.inside_nat_space( tmp_nat_coords, epsilon ) ) continue;
00947             }
00948             /*
00949             else if (etype == MBTRI){
00950               Element::LinearTri trimap(coords_vert);
00951               if (!trimap.inside_box( pos, epsilon))
00952                 continue;
00953               try {
00954                 tmp_nat_coords = trimap.ievaluate(pos, epsilon);
00955                 bool inside = trimap.inside_nat_space(tmp_nat_coords, epsilon);
00956                 if (!inside) continue;
00957               }
00958               catch (Element::Map::EvaluationError) {
00959                 continue;
00960               }
00961               if (!trimap.inside_nat_space(tmp_nat_coords, epsilon))
00962                 continue;
00963             }
00964             */
00965             else if( etype == MBEDGE )
00966             {
00967                 Element::LinearEdge edgemap( coords_vert );
00968                 try
00969                 {
00970                     tmp_nat_coords = edgemap.ievaluate( CartVect( xyz ), epsilon );
00971                 }
00972                 catch( Element::Map::EvaluationError )
00973                 {
00974                     continue;
00975                 }
00976                 if( !edgemap.inside_nat_space( tmp_nat_coords, epsilon ) ) continue;
00977             }
00978             else
00979             {
00980                 std::cout << "Entity not Hex/Tet/Quad/Tri/Edge. Please verify." << std::endl;
00981                 continue;
00982             }
00983         }
00984         // If we get here then we've found the coordinates.
00985         // Save them and the entity and return success.
00986         entities.push_back( *iter );
00987         nat_coords.push_back( tmp_nat_coords );
00988         return MB_SUCCESS;
00989     }
00990 
00991     // Didn't find any elements containing the point
00992     return MB_SUCCESS;
00993 }
00994 
00995 ErrorCode Coupler::interp_field( EntityHandle elem, CartVect nat_coord, Tag tag, double& field )
00996 {
00997     if( _spectralSource )
00998     {
00999         // Get tag values at the GL points for some field (Tag)
01000         const double* vx;
01001         ErrorCode rval = mbImpl->tag_get_by_ptr( tag, &elem, 1, (const void**)&vx );
01002         if( moab::MB_SUCCESS != rval )
01003         {
01004             std::cout << "Can't get field values for the tag \n";
01005             return MB_FAILURE;
01006         }
01007         Element::SpectralHex* spcHex = (Element::SpectralHex*)_spectralSource;
01008         field                        = spcHex->evaluate_scalar_field( nat_coord, vx );
01009     }
01010     else
01011     {
01012         double vfields[27];  // Will work for linear hex, quadratic hex or Tets
01013         moab::Element::Map* elemMap = NULL;
01014         int num_verts               = 0;
01015         // Get the EntityType
01016         // Get the tag values at the vertices
01017         const EntityHandle* connect;
01018         int num_connect;
01019         ErrorCode result = mbImpl->get_connectivity( elem, connect, num_connect );
01020         if( MB_SUCCESS != result ) return result;
01021         EntityType etype = mbImpl->type_from_handle( elem );
01022         if( MBHEX == etype )
01023         {
01024             if( 8 == num_connect )
01025             {
01026                 elemMap   = new moab::Element::LinearHex();
01027                 num_verts = 8;
01028             }
01029             else
01030             { /* (MBHEX == etype && 27 == num_connect) */
01031                 elemMap   = new moab::Element::QuadraticHex();
01032                 num_verts = 27;
01033             }
01034         }
01035         else if( MBTET == etype )
01036         {
01037             elemMap   = new moab::Element::LinearTet();
01038             num_verts = 4;
01039         }
01040         else if( MBQUAD == etype )
01041         {
01042             elemMap   = new moab::Element::LinearQuad();
01043             num_verts = 4;
01044         }
01045         else if( MBTRI == etype )
01046         {
01047             elemMap   = new moab::Element::LinearTri();
01048             num_verts = 3;
01049         }
01050         else
01051             return MB_FAILURE;
01052 
01053         result = mbImpl->tag_get_data( tag, connect, std::min( num_verts, num_connect ), vfields );
01054         if( MB_SUCCESS != result )
01055         {
01056             delete elemMap;
01057             return result;
01058         }
01059 
01060         // Function for the interpolation
01061         field = 0;
01062 
01063         // Check the number of vertices
01064         assert( num_connect >= num_verts );
01065 
01066         // Calculate the field
01067         try
01068         {
01069             field = elemMap->evaluate_scalar_field( nat_coord, vfields );
01070         }
01071         catch( moab::Element::Map::EvaluationError& )
01072         {
01073             delete elemMap;
01074             return MB_FAILURE;
01075         }
01076 
01077         delete elemMap;
01078     }
01079 
01080     return MB_SUCCESS;
01081 }
01082 
01083 // Simplest "interpolation" for element-based source fields. Set the value of the field
01084 // at the target point to that of the field in the source element it lies in.
01085 ErrorCode Coupler::constant_interp( EntityHandle elem, Tag tag, double& field )
01086 {
01087     double tempField;
01088 
01089     // Get the tag values at the vertices
01090     ErrorCode result = mbImpl->tag_get_data( tag, &elem, 1, &tempField );
01091     if( MB_SUCCESS != result ) return result;
01092 
01093     field = tempField;
01094 
01095     return MB_SUCCESS;
01096 }
01097 
01098 // Normalize a field over the entire mesh represented by the root_set.
01099 ErrorCode Coupler::normalize_mesh( EntityHandle root_set, const char* norm_tag, Coupler::IntegType integ_type,
01100                                    int num_integ_pts )
01101 {
01102     ErrorCode err;
01103 
01104     // SLAVE START ****************************************************************
01105     // Search for entities based on tag_handles and tag_values
01106     std::vector< std::vector< EntityHandle > > entity_sets;
01107     std::vector< std::vector< EntityHandle > > entity_groups;
01108 
01109     // put the root_set into entity_sets
01110     std::vector< EntityHandle > ent_set;
01111     ent_set.push_back( root_set );
01112     entity_sets.push_back( ent_set );
01113 
01114     // get all entities from root_set and put into entity_groups
01115     std::vector< EntityHandle > entities;
01116     err = mbImpl->get_entities_by_handle( root_set, entities, true );ERRORR( "Failed to get entities in root_set.", err );
01117 
01118     entity_groups.push_back( entities );
01119 
01120     // Call do_normalization() to continue common normalization processing
01121     err = do_normalization( norm_tag, entity_sets, entity_groups, integ_type, num_integ_pts );ERRORR( "Failure in do_normalization().", err );
01122     // SLAVE END   ****************************************************************
01123 
01124     return err;
01125 }
01126 
01127 // Normalize a field over the subset of entities identified by the tags and values passed
01128 ErrorCode Coupler::normalize_subset( EntityHandle root_set, const char* norm_tag, const char** tag_names, int num_tags,
01129                                      const char** tag_values, Coupler::IntegType integ_type, int num_integ_pts )
01130 {
01131     moab::ErrorCode err;
01132     std::vector< Tag > tag_handles;
01133 
01134     // Lookup tag handles from tag names
01135     for( int t = 0; t < num_tags; t++ )
01136     {
01137         // get tag handle & size
01138         Tag th;
01139         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 );
01140         tag_handles.push_back( th );
01141     }
01142 
01143     return normalize_subset( root_set, norm_tag, &tag_handles[0], num_tags, tag_values, integ_type, num_integ_pts );
01144 }
01145 
01146 ErrorCode Coupler::normalize_subset( EntityHandle root_set, const char* norm_tag, Tag* tag_handles, int num_tags,
01147                                      const char** tag_values, Coupler::IntegType integ_type, int num_integ_pts )
01148 {
01149     ErrorCode err;
01150 
01151     // SLAVE START ****************************************************************
01152     // Search for entities based on tag_handles and tag_values
01153     std::vector< std::vector< EntityHandle > > entity_sets;
01154     std::vector< std::vector< EntityHandle > > entity_groups;
01155 
01156     err = get_matching_entities( root_set, tag_handles, tag_values, num_tags, &entity_sets, &entity_groups );ERRORR( "Failed to get matching entities.", err );
01157 
01158     // Call do_normalization() to continue common normalization processing
01159     err = do_normalization( norm_tag, entity_sets, entity_groups, integ_type, num_integ_pts );ERRORR( "Failure in do_normalization().", err );
01160     // SLAVE END   ****************************************************************
01161 
01162     return err;
01163 }
01164 
01165 ErrorCode Coupler::do_normalization( const char* norm_tag, std::vector< std::vector< EntityHandle > >& entity_sets,
01166                                      std::vector< std::vector< EntityHandle > >& entity_groups,
01167                                      Coupler::IntegType integ_type, int num_integ_pts )
01168 {
01169     // SLAVE START ****************************************************************
01170     ErrorCode err;
01171     int ierr = 0;
01172 
01173     // Setup data for parallel computing
01174     int nprocs, rank;
01175     ierr = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
01176     ERRORMPI( "Getting number of procs failed.", ierr );
01177     ierr = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
01178     ERRORMPI( "Getting rank failed.", ierr );
01179 
01180     // Get the integrated field value for each group(vector) of entities.
01181     // If no entities are in a group then a zero will be put in the list
01182     // of return values.
01183     unsigned int num_ent_grps = entity_groups.size();
01184     std::vector< double > integ_vals( num_ent_grps );
01185 
01186     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 );
01187     // SLAVE END   ****************************************************************
01188 
01189     // SLAVE/MASTER START #########################################################
01190     // Send list of integrated values back to master proc. The ordering of the
01191     // values will match the ordering of the entity groups (i.e. vector of vectors)
01192     // sent from master to slaves earlier.  The values for each entity group will
01193     // be summed during the transfer.
01194     std::vector< double > sum_integ_vals( num_ent_grps );
01195 
01196     if( nprocs > 1 )
01197     {
01198         // If parallel then send the values back to the master.
01199         ierr = MPI_Reduce( &integ_vals[0], &sum_integ_vals[0], num_ent_grps, MPI_DOUBLE, MPI_SUM, MASTER_PROC,
01200                            myPc->proc_config().proc_comm() );
01201         ERRORMPI( "Transfer and reduction of integrated values failed.", ierr );
01202     }
01203     else
01204     {
01205         // Otherwise just copy the vector
01206         sum_integ_vals = integ_vals;
01207     }
01208     // SLAVE/MASTER END   #########################################################
01209 
01210     // MASTER START ***************************************************************
01211     // Calculate the normalization factor for each group by taking the
01212     // inverse of each integrated field value. Put the normalization factor
01213     // for each group back into the list in the same order.
01214     for( unsigned int i = 0; i < num_ent_grps; i++ )
01215     {
01216         double val = sum_integ_vals[i];
01217         if( fabs( val ) > 1e-8 )
01218             sum_integ_vals[i] = 1.0 / val;
01219         else
01220         {
01221             sum_integ_vals[i] = 0.0; /* VSM: not sure what we should do here ? */
01222             /* commenting out error below since if integral(value)=0.0, then normalization
01223                is probably unnecessary to start with ? */
01224             /* ERRORR("Integrating an invalid field -- integral("<<norm_tag<<") = "<<val<<".", err);
01225              */
01226         }
01227     }
01228     // MASTER END   ***************************************************************
01229 
01230     // MASTER/SLAVE START #########################################################
01231     if( nprocs > 1 )
01232     {
01233         // If parallel then broadcast the normalization factors to the procs.
01234         ierr = MPI_Bcast( &sum_integ_vals[0], num_ent_grps, MPI_DOUBLE, MASTER_PROC, myPc->proc_config().proc_comm() );
01235         ERRORMPI( "Broadcast of normalization factors failed.", ierr );
01236     }
01237     // MASTER/SLAVE END   #########################################################
01238 
01239     // SLAVE START ****************************************************************
01240     // Save the normalization factors to a new tag with name of norm_tag's value
01241     // and the string "_normF" appended.  This new tag will be created on the entity
01242     // set that contains all of the entities from a group.
01243 
01244     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 );
01245     // SLAVE END   ****************************************************************
01246 
01247     return err;
01248 }
01249 
01250 // Functions supporting the subset normalization function
01251 
01252 // Retrieve groups of entities matching tags and values if present
01253 ErrorCode Coupler::get_matching_entities( EntityHandle root_set, const char** tag_names, const char** tag_values,
01254                                           int num_tags, std::vector< std::vector< EntityHandle > >* entity_sets,
01255                                           std::vector< std::vector< EntityHandle > >* entity_groups )
01256 {
01257     ErrorCode err;
01258     std::vector< Tag > tag_handles;
01259 
01260     for( int t = 0; t < num_tags; t++ )
01261     {
01262         // Get tag handle & size
01263         Tag th;
01264         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 );
01265         tag_handles.push_back( th );
01266     }
01267 
01268     return get_matching_entities( root_set, &tag_handles[0], tag_values, num_tags, entity_sets, entity_groups );
01269 }
01270 
01271 // Retrieve groups of entities matching tags and values if present
01272 ErrorCode Coupler::get_matching_entities( EntityHandle root_set, Tag* tag_handles, const char** tag_values,
01273                                           int num_tags, std::vector< std::vector< EntityHandle > >* entity_sets,
01274                                           std::vector< std::vector< EntityHandle > >* entity_groups )
01275 {
01276     // SLAVE START ****************************************************************
01277 
01278     // Setup data for parallel computing
01279     ErrorCode err;
01280     int ierr = 0;
01281     int nprocs, rank;
01282     ierr = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
01283     ERRORMPI( "Getting number of procs failed.", ierr );
01284     ierr = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
01285     ERRORMPI( "Getting rank failed.", ierr );
01286 
01287     Range ent_sets;
01288     err =
01289         mbImpl->get_entities_by_type_and_tag( root_set, moab::MBENTITYSET, tag_handles, (const void* const*)tag_values,
01290                                               num_tags, ent_sets, Interface::INTERSECT, 0 );ERRORR( "Core::get_entities_by_type_and_tag failed.", err );
01291 
01292     TupleList* tag_list = NULL;
01293     err                 = create_tuples( ent_sets, tag_handles, num_tags, &tag_list );ERRORR( "Failed to create tuples from entity sets.", err );
01294 
01295     // Free up range
01296     ent_sets.clear();
01297     // SLAVE END   ****************************************************************
01298 
01299     // If we are running in a multi-proc session then send tuple list back to master
01300     // proc for consolidation. Otherwise just copy the pointer to the tuple_list.
01301     TupleList* cons_tuples;
01302     if( nprocs > 1 )
01303     {
01304         // SLAVE/MASTER START #########################################################
01305 
01306         // Pack the tuple_list in a buffer.
01307         uint* tuple_buf;
01308         int tuple_buf_len;
01309         tuple_buf_len = pack_tuples( tag_list, (void**)&tuple_buf );
01310 
01311         // Free tag_list here as its not used again if nprocs > 1
01312         tag_list->reset();
01313 
01314         // Send back the buffer sizes to the master proc
01315         int* recv_cnts       = (int*)malloc( nprocs * sizeof( int ) );
01316         int* offsets         = (int*)malloc( nprocs * sizeof( int ) );
01317         uint* all_tuples_buf = NULL;
01318 
01319         MPI_Gather( &tuple_buf_len, 1, MPI_INT, recv_cnts, 1, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm() );
01320         ERRORMPI( "Gathering buffer sizes failed.", err );
01321 
01322         // Allocate a buffer large enough for all the data
01323         if( rank == MASTER_PROC )
01324         {
01325             int all_tuples_len = recv_cnts[0];
01326             offsets[0]         = 0;
01327             for( int i = 1; i < nprocs; i++ )
01328             {
01329                 offsets[i] = offsets[i - 1] + recv_cnts[i - 1];
01330                 all_tuples_len += recv_cnts[i];
01331             }
01332 
01333             all_tuples_buf = (uint*)malloc( all_tuples_len * sizeof( uint ) );
01334         }
01335 
01336         // Send all buffers to the master proc for consolidation
01337         MPI_Gatherv( (void*)tuple_buf, tuple_buf_len, MPI_INT, (void*)all_tuples_buf, recv_cnts, offsets, MPI_INT,
01338                      MASTER_PROC, myPc->proc_config().proc_comm() );
01339         ERRORMPI( "Gathering tuple_lists failed.", err );
01340         free( tuple_buf );  // malloc'd in pack_tuples
01341 
01342         if( rank == MASTER_PROC )
01343         {
01344             // Unpack the tuple_list from the buffer.
01345             TupleList** tl_array = (TupleList**)malloc( nprocs * sizeof( TupleList* ) );
01346             for( int i = 0; i < nprocs; i++ )
01347                 unpack_tuples( (void*)&all_tuples_buf[offsets[i]], &tl_array[i] );
01348 
01349             // Free all_tuples_buf here as it is only allocated on the MASTER_PROC
01350             free( all_tuples_buf );
01351             // SLAVE/MASTER END   #########################################################
01352 
01353             // MASTER START ***************************************************************
01354             // Consolidate all tuple_lists into one tuple_list with no duplicates.
01355             err = consolidate_tuples( tl_array, nprocs, &cons_tuples );ERRORR( "Failed to consolidate tuples.", err );
01356 
01357             for( int i = 0; i < nprocs; i++ )
01358                 tl_array[i]->reset();
01359             free( tl_array );
01360             // MASTER END   ***************************************************************
01361         }
01362 
01363         // Free offsets and recv_cnts as they are allocated on all procs
01364         free( offsets );
01365         free( recv_cnts );
01366 
01367         // MASTER/SLAVE START #########################################################
01368         // Broadcast condensed tuple list back to all procs.
01369         uint* ctl_buf;
01370         int ctl_buf_sz;
01371         if( rank == MASTER_PROC ) ctl_buf_sz = pack_tuples( cons_tuples, (void**)&ctl_buf );
01372 
01373         // Send buffer size
01374         ierr = MPI_Bcast( &ctl_buf_sz, 1, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm() );
01375         ERRORMPI( "Broadcasting tuple_list size failed.", ierr );
01376 
01377         // Allocate a buffer in the other procs
01378         if( rank != MASTER_PROC ) ctl_buf = (uint*)malloc( ctl_buf_sz * sizeof( uint ) );
01379 
01380         ierr = MPI_Bcast( (void*)ctl_buf, ctl_buf_sz, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm() );
01381         ERRORMPI( "Broadcasting tuple_list failed.", ierr );
01382 
01383         if( rank != MASTER_PROC ) unpack_tuples( ctl_buf, &cons_tuples );
01384         free( ctl_buf );
01385         // MASTER/SLAVE END   #########################################################
01386     }
01387     else
01388         cons_tuples = tag_list;
01389 
01390     // SLAVE START ****************************************************************
01391     // Loop over the tuple list getting the entities with the tags in the tuple_list entry
01392     uint mi, ml, mul, mr;
01393     cons_tuples->getTupleSize( mi, ml, mul, mr );
01394 
01395     for( unsigned int i = 0; i < cons_tuples->get_n(); i++ )
01396     {
01397         // Get Entity Sets that match the tags and values.
01398 
01399         // Convert the data in the tuple_list to an array of pointers to the data
01400         // in the tuple_list as that is what the iMesh API call is expecting.
01401         int** vals = (int**)malloc( mi * sizeof( int* ) );
01402         for( unsigned int j = 0; j < mi; j++ )
01403             vals[j] = (int*)&( cons_tuples->vi_rd[( i * mi ) + j] );
01404 
01405         // Get entities recursively based on type and tag data
01406         err = mbImpl->get_entities_by_type_and_tag( root_set, moab::MBENTITYSET, tag_handles, (const void* const*)vals,
01407                                                     mi, ent_sets, Interface::INTERSECT, 0 );ERRORR( "Core::get_entities_by_type_and_tag failed.", err );
01408         if( debug ) std::cout << "ent_sets_size=" << ent_sets.size() << std::endl;
01409 
01410         // Free up the array of pointers
01411         free( vals );
01412 
01413         // Loop over the entity sets and then free the memory for ent_sets.
01414         std::vector< EntityHandle > ent_set_hdls;
01415         std::vector< EntityHandle > ent_hdls;
01416         for( unsigned int j = 0; j < ent_sets.size(); j++ )
01417         {
01418             // Save the entity set
01419             ent_set_hdls.push_back( ent_sets[j] );
01420 
01421             // Get all entities for the entity set
01422             Range ents;
01423 
01424             /* VSM: do we need to filter out entity sets ? */
01425             err = mbImpl->get_entities_by_handle( ent_sets[j], ents, false );ERRORR( "Core::get_entities_by_handle failed.", err );
01426             if( debug ) std::cout << "ents_size=" << ents.size() << std::endl;
01427 
01428             // Save all of the entities from the entity set and free the memory for ents.
01429             for( unsigned int k = 0; k < ents.size(); k++ )
01430             {
01431                 ent_hdls.push_back( ents[k] );
01432             }
01433             ents.clear();
01434             if( debug ) std::cout << "ent_hdls.size=" << ent_hdls.size() << std::endl;
01435         }
01436 
01437         // Free the entity set list for next tuple iteration.
01438         ent_sets.clear();
01439 
01440         // Push ent_set_hdls onto entity_sets, ent_hdls onto entity_groups
01441         // and clear both ent_set_hdls and ent_hdls.
01442         entity_sets->push_back( ent_set_hdls );
01443         ent_set_hdls.clear();
01444         entity_groups->push_back( ent_hdls );
01445         ent_hdls.clear();
01446         if( debug )
01447             std::cout << "entity_sets->size=" << entity_sets->size()
01448                       << ", entity_groups->size=" << entity_groups->size() << std::endl;
01449     }
01450 
01451     cons_tuples->reset();
01452     // SLAVE END   ****************************************************************
01453 
01454     return err;
01455 }
01456 
01457 // Return a tuple_list containing  tag values for each Entity Set
01458 // The tuple_list will have a column for each tag and a row for each
01459 // Entity Set. It is assumed all of the tags are integer tags.
01460 ErrorCode Coupler::create_tuples( Range& ent_sets, const char** tag_names, unsigned int num_tags,
01461                                   TupleList** tuple_list )
01462 {
01463     ErrorCode err;
01464     std::vector< Tag > tag_handles;
01465 
01466     for( unsigned int t = 0; t < num_tags; t++ )
01467     {
01468         // Get tag handle & size
01469         Tag th;
01470         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 );
01471         tag_handles.push_back( th );
01472     }
01473 
01474     return create_tuples( ent_sets, &tag_handles[0], num_tags, tuple_list );
01475 }
01476 
01477 // Return a tuple_list containing  tag values for each Entity Set
01478 // The tuple_list will have a column for each tag and a row for each
01479 // Entity Set.  It is assumed all of the tags are integer tags.
01480 ErrorCode Coupler::create_tuples( Range& ent_sets, Tag* tag_handles, unsigned int num_tags, TupleList** tuples )
01481 {
01482     // ASSUMPTION: All tags are of type integer.  This may need to be expanded in future.
01483     ErrorCode err;
01484 
01485     // Allocate a tuple_list for the number of entity sets passed in
01486     TupleList* tag_tuples = new TupleList( num_tags, 0, 0, 0, (int)ent_sets.size() );
01487     // tag_tuples->initialize(num_tags, 0, 0, 0, num_sets);
01488     uint mi, ml, mul, mr;
01489     tag_tuples->getTupleSize( mi, ml, mul, mr );
01490     tag_tuples->enableWriteAccess();
01491 
01492     if( mi == 0 ) ERRORR( "Failed to initialize tuple_list.", MB_FAILURE );
01493 
01494     // Loop over the filtered entity sets retrieving each matching tag value one by one.
01495     int val;
01496     for( unsigned int i = 0; i < ent_sets.size(); i++ )
01497     {
01498         for( unsigned int j = 0; j < num_tags; j++ )
01499         {
01500             EntityHandle set_handle = ent_sets[i];
01501             err                     = mbImpl->tag_get_data( tag_handles[j], &set_handle, 1, &val );ERRORR( "Failed to get integer tag data.", err );
01502             tag_tuples->vi_wr[i * mi + j] = val;
01503         }
01504 
01505         // If we get here there was no error so increment n in the tuple_list
01506         tag_tuples->inc_n();
01507     }
01508     tag_tuples->disableWriteAccess();
01509     *tuples = tag_tuples;
01510 
01511     return MB_SUCCESS;
01512 }
01513 
01514 // Consolidate tuple_lists into one list with no duplicates
01515 ErrorCode Coupler::consolidate_tuples( TupleList** all_tuples, unsigned int num_tuples, TupleList** unique_tuples )
01516 {
01517     int total_rcv_tuples = 0;
01518     int offset = 0, copysz = 0;
01519     unsigned num_tags = 0;
01520 
01521     uint ml, mul, mr;
01522     uint* mi = (uint*)malloc( sizeof( uint ) * num_tuples );
01523 
01524     for( unsigned int i = 0; i < num_tuples; i++ )
01525     {
01526         all_tuples[i]->getTupleSize( mi[i], ml, mul, mr );
01527     }
01528 
01529     for( unsigned int i = 0; i < num_tuples; i++ )
01530     {
01531         if( all_tuples[i] != NULL )
01532         {
01533             total_rcv_tuples += all_tuples[i]->get_n();
01534             num_tags = mi[i];
01535         }
01536     }
01537     const unsigned int_size  = sizeof( sint );
01538     const unsigned int_width = num_tags * int_size;
01539 
01540     // Get the total size of all of the tuple_lists in all_tuples.
01541     for( unsigned int i = 0; i < num_tuples; i++ )
01542     {
01543         if( all_tuples[i] != NULL ) total_rcv_tuples += all_tuples[i]->get_n();
01544     }
01545 
01546     // Copy the tuple_lists into a single tuple_list.
01547     TupleList* all_tuples_list = new TupleList( num_tags, 0, 0, 0, total_rcv_tuples );
01548     all_tuples_list->enableWriteAccess();
01549     // all_tuples_list->initialize(num_tags, 0, 0, 0, total_rcv_tuples);
01550     for( unsigned int i = 0; i < num_tuples; i++ )
01551     {
01552         if( all_tuples[i] != NULL )
01553         {
01554             copysz = all_tuples[i]->get_n() * int_width;
01555             memcpy( all_tuples_list->vi_wr + offset, all_tuples[i]->vi_rd, copysz );
01556             offset = offset + ( all_tuples[i]->get_n() * mi[i] );
01557             all_tuples_list->set_n( all_tuples_list->get_n() + all_tuples[i]->get_n() );
01558         }
01559     }
01560 
01561     // Sort the new tuple_list. Use a radix type sort, starting with the last (or least significant)
01562     // tag column in the vi array and working towards the first (or most significant) tag column.
01563     TupleList::buffer sort_buffer;
01564     sort_buffer.buffer_init( 2 * total_rcv_tuples * int_width );
01565     for( int i = num_tags - 1; i >= 0; i-- )
01566     {
01567         all_tuples_list->sort( i, &sort_buffer );
01568     }
01569 
01570     // Cycle through the sorted list eliminating duplicates.
01571     // Keep counters to the current end of the tuple_list (w/out dups) and the last tuple examined.
01572     unsigned int end_idx = 0, last_idx = 1;
01573     while( last_idx < all_tuples_list->get_n() )
01574     {
01575         if( memcmp( all_tuples_list->vi_rd + ( end_idx * num_tags ), all_tuples_list->vi_rd + ( last_idx * num_tags ),
01576                     int_width ) == 0 )
01577         {
01578             // Values equal - skip
01579             last_idx += 1;
01580         }
01581         else
01582         {
01583             // Values different - copy
01584             // Move up the end index
01585             end_idx += 1;
01586             memcpy( all_tuples_list->vi_wr + ( end_idx * num_tags ), all_tuples_list->vi_rd + ( last_idx * num_tags ),
01587                     int_width );
01588             last_idx += 1;
01589         }
01590     }
01591     // Update the count in all_tuples_list
01592     all_tuples_list->set_n( end_idx + 1 );
01593 
01594     // Resize the tuple_list
01595     all_tuples_list->resize( all_tuples_list->get_n() );
01596 
01597     // Set the output parameter
01598     *unique_tuples = all_tuples_list;
01599 
01600     return MB_SUCCESS;
01601 }
01602 
01603 // Calculate integrated field values for groups of entities
01604 ErrorCode Coupler::get_group_integ_vals( std::vector< std::vector< EntityHandle > >& groups,
01605                                          std::vector< double >& integ_vals, const char* norm_tag,
01606                                          int /*num_integ_vals*/, Coupler::IntegType integ_type )
01607 {
01608     ErrorCode err;
01609 
01610     std::vector< std::vector< EntityHandle > >::iterator iter_i;
01611     std::vector< EntityHandle >::iterator iter_j;
01612     double grp_intrgr_val, intgr_val;
01613 
01614     // Get the tag handle for norm_tag
01615     Tag norm_hdl;
01616     err =
01617         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 );
01618 
01619     // Check size of integ_vals vector
01620     if( integ_vals.size() != groups.size() ) integ_vals.resize( groups.size() );
01621 
01622     // Loop over the groups(vectors) of entities
01623     unsigned int i;
01624     for( i = 0, iter_i = groups.begin(); iter_i != groups.end(); i++, ++iter_i )
01625     {
01626         grp_intrgr_val = 0;
01627 
01628         // Loop over the all the entities in the group, integrating
01629         // the field_fn over the entity in iter_j
01630         for( iter_j = ( *iter_i ).begin(); iter_j != ( *iter_i ).end(); ++iter_j )
01631         {
01632             EntityHandle ehandle = ( *iter_j );
01633 
01634             // Check that the entity in iter_j is of the same dimension as the
01635             // integ_type we are performing
01636             EntityType j_type;
01637             j_type = mbImpl->type_from_handle( ehandle );ERRORR( "Failed to get entity type.", err );
01638             // Skip any entities in the group that are not of the type being considered
01639             if( ( integ_type == VOLUME ) && ( j_type < MBTET || j_type >= MBENTITYSET ) ) continue;
01640 
01641             intgr_val = 0;
01642 
01643             // Retrieve the vertices from the element
01644             const EntityHandle* verts = NULL;
01645             int connectivity_size     = 0;
01646 
01647             err = mbImpl->get_connectivity( ehandle, verts, connectivity_size, false );ERRORR( "Failed to get vertices from entity.", err );
01648 
01649             // Get the vertex coordinates and the field values at the vertices.
01650             double* coords = (double*)malloc( sizeof( double ) * ( 3 * connectivity_size ) );
01651             /* TODO: VSM: check if this works for lower dimensions also without problems */
01652             /* if (3 == geom_dim) */
01653             err = mbImpl->get_coords( verts, connectivity_size, coords );ERRORR( "Failed to get vertex coordinates.", err );
01654 
01655             /* allocate the field data array */
01656             double* vfield = (double*)malloc( sizeof( double ) * ( connectivity_size ) );
01657             err            = mbImpl->tag_get_data( norm_hdl, verts, connectivity_size, vfield );
01658             if( MB_SUCCESS != err ) { free( coords ); }
01659             ERRORR( "Failed to get vertex coordinates.", err );
01660 
01661             // Get coordinates of all corner vertices (in normal order) and
01662             // put in array of CartVec.
01663             std::vector< CartVect > vertices( connectivity_size );
01664 
01665             // Put the vertices into a CartVect vector
01666             double* x = coords;
01667             for( int j = 0; j < connectivity_size; j++, x += 3 )
01668             {
01669                 vertices[j] = CartVect( x );
01670             }
01671             free( coords );
01672 
01673             moab::Element::Map* elemMap;
01674             if( j_type == MBHEX )
01675             {
01676                 if( connectivity_size == 8 )
01677                     elemMap = new moab::Element::LinearHex( vertices );
01678                 else
01679                     elemMap = new moab::Element::QuadraticHex( vertices );
01680             }
01681             else if( j_type == MBTET )
01682             {
01683                 elemMap = new moab::Element::LinearTet( vertices );
01684             }
01685             else if( j_type == MBQUAD )
01686             {
01687                 elemMap = new moab::Element::LinearQuad( vertices );
01688             }
01689             /*
01690             else if (j_type == MBTRI) {
01691               elemMap = new moab::Element::LinearTri(vertices);
01692             }
01693             */
01694             else if( j_type == MBEDGE )
01695             {
01696                 elemMap = new moab::Element::LinearEdge( vertices );
01697             }
01698             else
01699                 ERRORR( "Unknown topology type.", MB_UNSUPPORTED_OPERATION );
01700 
01701             // Set the vertices in the Map and perform the integration
01702             try
01703             {
01704                 /* VSM: Do we need this call ?? */
01705                 // elemMap->set_vertices(vertices);
01706 
01707                 // Perform the actual integration over the element
01708                 intgr_val = elemMap->integrate_scalar_field( vfield );
01709 
01710                 // Combine the result with those of the group
01711                 grp_intrgr_val += intgr_val;
01712             }
01713             catch( moab::Element::Map::ArgError& )
01714             {
01715                 std::cerr << "Failed to set vertices on Element::Map." << std::endl;
01716             }
01717             catch( moab::Element::Map::EvaluationError& )
01718             {
01719                 std::cerr << "Failed to get inverse evaluation of coordinate on Element::Map." << std::endl;
01720             }
01721 
01722             delete( elemMap );
01723             free( vfield );
01724         }
01725 
01726         // Set the group integrated value in the vector
01727         integ_vals[i] = grp_intrgr_val;
01728     }
01729 
01730     return err;
01731 }
01732 
01733 // Apply a normalization factor to group of entities
01734 ErrorCode Coupler::apply_group_norm_factor( std::vector< std::vector< EntityHandle > >& entity_sets,
01735                                             std::vector< double >& norm_factors, const char* norm_tag,
01736                                             Coupler::IntegType /*integ_type*/ )
01737 {
01738     ErrorCode err;
01739 
01740     // Construct the new tag for the normalization factor from the norm_tag name
01741     // and "_normf".
01742     int norm_tag_len       = strlen( norm_tag );
01743     const char* normf_appd = "_normf";
01744     int normf_appd_len     = strlen( normf_appd );
01745 
01746     char* normf_tag = (char*)malloc( norm_tag_len + normf_appd_len + 1 );
01747     char* tmp_ptr   = normf_tag;
01748 
01749     memcpy( tmp_ptr, norm_tag, norm_tag_len );
01750     tmp_ptr += norm_tag_len;
01751     memcpy( tmp_ptr, normf_appd, normf_appd_len );
01752     tmp_ptr += normf_appd_len;
01753     *tmp_ptr = '\0';
01754 
01755     Tag normf_hdl;
01756     // Check to see if the tag exists.  If not then create it and get the handle.
01757     err = mbImpl->tag_get_handle( normf_tag, 1, moab::MB_TYPE_DOUBLE, normf_hdl,
01758                                   moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT );ERRORR( "Failed to create normalization factor tag.", err );
01759     if( normf_hdl == NULL )
01760     {
01761         std::string msg( "Failed to create normalization factor tag named '" );
01762         msg += std::string( normf_tag ) + std::string( "'" );ERRORR( msg.c_str(), MB_FAILURE );
01763     }
01764     free( normf_tag );
01765 
01766     std::vector< std::vector< EntityHandle > >::iterator iter_i;
01767     std::vector< EntityHandle >::iterator iter_j;
01768     std::vector< double >::iterator iter_f;
01769     double grp_norm_factor = 0.0;
01770 
01771     // Loop over the entity sets
01772     for( iter_i = entity_sets.begin(), iter_f = norm_factors.begin();
01773          ( iter_i != entity_sets.end() ) && ( iter_f != norm_factors.end() ); ++iter_i, ++iter_f )
01774     {
01775         grp_norm_factor = *iter_f;
01776 
01777         // Loop over the all the entity sets in iter_i and set the
01778         // new normf_tag with the norm factor value on each
01779         for( iter_j = ( *iter_i ).begin(); iter_j != ( *iter_i ).end(); ++iter_j )
01780         {
01781             EntityHandle entset = *iter_j;
01782 
01783             std::cout << "Coupler: applying normalization for entity set=" << entset
01784                       << ",  normalization_factor=" << grp_norm_factor << std::endl;
01785 
01786             err = mbImpl->tag_set_data( normf_hdl, &entset, 1, &grp_norm_factor );ERRORR( "Failed to set normalization factor on entity set.", err );
01787         }
01788     }
01789 
01790     return MB_SUCCESS;
01791 }
01792 
01793 #define UINT_PER_X( X )   ( ( sizeof( X ) + sizeof( uint ) - 1 ) / sizeof( uint ) )
01794 #define UINT_PER_REAL     UINT_PER_X( realType )
01795 #define UINT_PER_LONG     UINT_PER_X( slong )
01796 #define UINT_PER_UNSIGNED UINT_PER_X( unsigned )
01797 
01798 // Function for packing tuple_list.  Returns number of uints copied into buffer.
01799 int pack_tuples( TupleList* tl, void** ptr )
01800 {
01801     uint mi, ml, mul, mr;
01802     tl->getTupleSize( mi, ml, mul, mr );
01803 
01804     uint n = tl->get_n();
01805 
01806     int sz_buf = 1 + 4 * UINT_PER_UNSIGNED +
01807                  tl->get_n() * ( mi + ml * UINT_PER_LONG + mul * UINT_PER_LONG + mr * UINT_PER_REAL );
01808 
01809     uint* buf = (uint*)malloc( sz_buf * sizeof( uint ) );
01810     *ptr      = (void*)buf;
01811 
01812     // Copy n
01813     memcpy( buf, &n, sizeof( uint ) ), buf += 1;
01814     // Copy mi
01815     memcpy( buf, &mi, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
01816     // Copy ml
01817     memcpy( buf, &ml, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
01818     // Copy mul
01819     memcpy( buf, &mul, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
01820     // Copy mr
01821     memcpy( buf, &mr, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
01822     // Copy vi_wr
01823     memcpy( buf, tl->vi_rd, tl->get_n() * mi * sizeof( sint ) ), buf += tl->get_n() * mi;
01824     // Copy vl_wr
01825     memcpy( buf, tl->vl_rd, tl->get_n() * ml * sizeof( slong ) ), buf += tl->get_n() * ml * UINT_PER_LONG;
01826     // Copy vul_wr
01827     memcpy( buf, tl->vul_rd, tl->get_n() * mul * sizeof( ulong ) ), buf += tl->get_n() * mul * UINT_PER_LONG;
01828     // Copy vr_wr
01829     memcpy( buf, tl->vr_rd, tl->get_n() * mr * sizeof( realType ) ), buf += tl->get_n() * mr * UINT_PER_REAL;
01830 
01831     return sz_buf;
01832 }
01833 
01834 // Function for packing tuple_list
01835 void unpack_tuples( void* ptr, TupleList** tlp )
01836 {
01837     TupleList* tl = new TupleList();
01838     *tlp          = tl;
01839 
01840     uint nt;
01841     unsigned mit, mlt, mult, mrt;
01842     uint* buf = (uint*)ptr;
01843 
01844     // Get n
01845     memcpy( &nt, buf, sizeof( uint ) ), buf += 1;
01846     // Get mi
01847     memcpy( &mit, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
01848     // Get ml
01849     memcpy( &mlt, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
01850     // Get mul
01851     memcpy( &mult, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
01852     // Get mr
01853     memcpy( &mrt, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
01854 
01855     // Initialize tl
01856     tl->initialize( mit, mlt, mult, mrt, nt );
01857     tl->enableWriteAccess();
01858     tl->set_n( nt );
01859 
01860     uint mi, ml, mul, mr;
01861     tl->getTupleSize( mi, ml, mul, mr );
01862 
01863     // Get vi_rd
01864     memcpy( tl->vi_wr, buf, tl->get_n() * mi * sizeof( sint ) ), buf += tl->get_n() * mi;
01865     // Get vl_rd
01866     memcpy( tl->vl_wr, buf, tl->get_n() * ml * sizeof( slong ) ), buf += tl->get_n() * ml * UINT_PER_LONG;
01867     // Get vul_rd
01868     memcpy( tl->vul_wr, buf, tl->get_n() * mul * sizeof( ulong ) ), buf += tl->get_n() * mul * UINT_PER_LONG;
01869     // Get vr_rd
01870     memcpy( tl->vr_wr, buf, tl->get_n() * mr * sizeof( realType ) ), buf += tl->get_n() * mr * UINT_PER_REAL;
01871 
01872     tl->disableWriteAccess();
01873     return;
01874 }
01875 
01876 ErrorCode Coupler::get_gl_points_on_elements( Range& targ_elems, std::vector< double >& vpos, int& numPointsOfInterest )
01877 {
01878     numPointsOfInterest = targ_elems.size() * _ntot;
01879     vpos.resize( 3 * numPointsOfInterest );
01880     int ielem = 0;
01881     for( Range::iterator eit = targ_elems.begin(); eit != targ_elems.end(); ++eit, ielem += _ntot * 3 )
01882     {
01883         EntityHandle eh = *eit;
01884         const double* xval;
01885         const double* yval;
01886         const double* zval;
01887         ErrorCode rval = mbImpl->tag_get_by_ptr( _xm1Tag, &eh, 1, (const void**)&xval );
01888         if( moab::MB_SUCCESS != rval )
01889         {
01890             std::cout << "Can't get xm1 values \n";
01891             return MB_FAILURE;
01892         }
01893         rval = mbImpl->tag_get_by_ptr( _ym1Tag, &eh, 1, (const void**)&yval );
01894         if( moab::MB_SUCCESS != rval )
01895         {
01896             std::cout << "Can't get ym1 values \n";
01897             return MB_FAILURE;
01898         }
01899         rval = mbImpl->tag_get_by_ptr( _zm1Tag, &eh, 1, (const void**)&zval );
01900         if( moab::MB_SUCCESS != rval )
01901         {
01902             std::cout << "Can't get zm1 values \n";
01903             return MB_FAILURE;
01904         }
01905         // Now, in a stride, populate vpos
01906         for( int i = 0; i < _ntot; i++ )
01907         {
01908             vpos[ielem + 3 * i]     = xval[i];
01909             vpos[ielem + 3 * i + 1] = yval[i];
01910             vpos[ielem + 3 * i + 2] = zval[i];
01911         }
01912     }
01913 
01914     return MB_SUCCESS;
01915 }
01916 
01917 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines