MOAB: Mesh Oriented datABase  (version 5.2.1)
mbslavepart.cpp
Go to the documentation of this file.
00001 //
00002 // Usage:
00003 // tools/mbslavepart -d 2 -m mpas/x1.2562.grid.h5m -s mpas/x1.10242.grid.h5m -o mpas_slave.h5m -e
00004 // 1e-8 -b 1e-6 -O
00005 //
00006 #include <iostream>
00007 #include <exception>
00008 #include <cmath>
00009 #include <vector>
00010 #include <string>
00011 
00012 #include "moab/ProgOptions.hpp"
00013 #include "moab/Core.hpp"
00014 #include "moab/AdaptiveKDTree.hpp"
00015 #include "moab/BVHTree.hpp"
00016 
00017 #include "moab/IntxMesh/IntxUtils.hpp"
00018 
00019 #ifdef MOAB_HAVE_MPI
00020 #include "moab_mpi.h"
00021 #include "moab/ParallelComm.hpp"
00022 #include "MBParallelConventions.h"
00023 #endif
00024 
00025 using namespace moab;
00026 
00027 // A function to get the non-default value from a std::map
00028 template < typename K, typename V >
00029 static V get_map_value( const std::map< K, V >& m, const K& key, const V& defval )
00030 {
00031     typename std::map< K, V >::const_iterator it = m.find( key );
00032     if( it == m.end() ) { return defval; }
00033     else
00034     {
00035         return it->second;
00036     }
00037 }
00038 
00039 int main( int argc, char* argv[] )
00040 {
00041     int proc_id = 0, size = 1, dimension = 3;
00042 #ifdef MOAB_HAVE_MPI
00043     MPI_Init( &argc, &argv );
00044     MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
00045     MPI_Comm_size( MPI_COMM_WORLD, &size );
00046 #endif
00047 
00048     int defaultpart  = 0;
00049     double tolerance = 1e-2, treetolerance = 1e-13, btolerance = 1e-4;
00050     std::string masterfile, slavefile, outfile( "slavemesh.h5m" );
00051     bool keepsparts    = false;
00052     bool use_spherical = false;
00053     ProgOptions opts;
00054 
00055     opts.addOpt< std::string >( "master,m", "Master mesh filename", &masterfile );
00056     opts.addOpt< std::string >( "slave,s", "Slave mesh filename", &slavefile );
00057     opts.addOpt< std::string >( "output,o", "Output partitioned mesh filename", &outfile );
00058     opts.addOpt< int >( "dim,d", "Dimension of entities to use for partitioning", &dimension );
00059     opts.addOpt< int >( "defaultpart,p", "Default partition number if target element is not found on source grid",
00060                         &defaultpart );
00061     opts.addOpt< double >( "eps,e", "Tolerance for the point search", &tolerance );
00062     opts.addOpt< double >( "beps,b", "Tolerance for the bounding box search", &btolerance );
00063     opts.addOpt< void >( "keep,K",
00064                          "Keep the existing partitions in the slave mesh (use PARALLEL_PARTITION_SLAVE instead)",
00065                          &keepsparts );
00066     opts.addOpt< void >( "spherical", "Hint that the meshes are defined on a spherical surface (Climate problems)",
00067                          &use_spherical );
00068     opts.parseCommandLine( argc, argv );
00069 
00070     if( masterfile.empty() || slavefile.empty() )
00071     {
00072         opts.printHelp();
00073 #ifdef MOAB_HAVE_MPI
00074         MPI_Finalize();
00075 #endif
00076         exit( 1 );
00077     }
00078 
00079     ErrorCode error;
00080     Core* mbCore = new Core();
00081 
00082     // Set the read options for parallel file loading
00083     const std::string partition_set_name = "PARALLEL_PARTITION";
00084     const std::string global_id_name     = "GLOBAL_ID";
00085     std::vector< std::string > read_opts, write_opts;
00086     std::string read_options( "" ), write_options( "" );
00087 
00088     if( size > 1 )
00089     {
00090         read_options  = "PARALLEL=READ_PART;PARTITION=" + partition_set_name + ";PARALLEL_RESOLVE_SHARED_ENTS";
00091         write_options = "PARALLEL=WRITE_PART";
00092     }
00093 
00094     EntityHandle masterfileset, slavefileset;
00095     error = mbCore->create_meshset( moab::MESHSET_TRACK_OWNER | moab::MESHSET_SET, masterfileset );MB_CHK_ERR( error );
00096     error = mbCore->create_meshset( moab::MESHSET_TRACK_OWNER | moab::MESHSET_SET, slavefileset );MB_CHK_ERR( error );
00097 
00098     // Load file
00099     error = mbCore->load_file( masterfile.c_str(), &masterfileset, read_options.c_str() );MB_CHK_ERR( error );
00100     error = mbCore->load_file( slavefile.c_str(), &slavefileset, read_options.c_str() );MB_CHK_ERR( error );
00101     // if (error != MB_SUCCESS && size > 1)
00102     // {
00103     //   std::string newread_options = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS";
00104     //   error = mbCore->load_file(slavefile.c_str(), &slavefileset, newread_options.c_str());
00105     // }
00106     // else MB_CHK_ERR(error);
00107 
00108     Tag gidtag = 0, parttag = 0, sparttag = 0;
00109     int dum_id = -1;
00110     error      = mbCore->tag_get_handle( partition_set_name.c_str(), 1, MB_TYPE_INTEGER, parttag,
00111                                     MB_TAG_SPARSE | MB_TAG_CREAT, &dum_id );MB_CHK_ERR( error );
00112     gidtag = mbCore->globalId_tag();
00113     if( keepsparts )
00114     {
00115         error = mbCore->tag_get_handle( std::string( partition_set_name + "_SLAVE" ).c_str(), 1, MB_TYPE_INTEGER,
00116                                         sparttag, MB_TAG_CREAT | MB_TAG_SPARSE, &dum_id );MB_CHK_ERR( error );
00117     }
00118 
00119     Range melems, msets, selems, ssets;
00120 
00121     // Get the partition sets on the master mesh
00122     std::map< int, int > mpartvals;
00123     error = mbCore->get_entities_by_type_and_tag( masterfileset, MBENTITYSET, &parttag, NULL, 1, msets,
00124                                                   moab::Interface::UNION, true );MB_CHK_ERR( error );
00125     if( msets.size() == 0 )
00126     {
00127         std::cout << "No partition sets found in the master mesh. Quitting..." << std::endl;
00128         exit( 1 );
00129     }
00130 
00131     for( unsigned i = 0; i < msets.size(); ++i )
00132     {
00133         EntityHandle mset = msets[i];
00134 
00135         moab::Range msetelems;
00136         error = mbCore->get_entities_by_dimension( mset, dimension, msetelems );MB_CHK_ERR( error );
00137         melems.merge( msetelems );
00138 
00139         int partID;
00140         error = mbCore->tag_get_data( parttag, &mset, 1, &partID );MB_CHK_ERR( error );
00141 
00142         // Get the global ID and use that as the indicator
00143         std::vector< int > gidMelems( msetelems.size() );
00144         error = mbCore->tag_get_data( gidtag, msetelems, gidMelems.data() );MB_CHK_ERR( error );
00145 
00146         for( unsigned j = 0; j < msetelems.size(); ++j )
00147             mpartvals[gidMelems[j]] = partID;
00148             // mpartvals[msetelems[j]]=partID;
00149 #ifdef VERBOSE
00150         std::cout << "Part " << partID << " has " << msetelems.size() << " elements." << std::endl;
00151 #endif
00152     }
00153 
00154     // Get information about the slave file set
00155     error = mbCore->get_entities_by_type_and_tag( slavefileset, MBENTITYSET, &parttag, NULL, 1, ssets,
00156                                                   moab::Interface::UNION );MB_CHK_ERR( error );
00157     // TODO: expand and add other dimensional elements
00158     error = mbCore->get_entities_by_dimension( slavefileset, dimension, selems );MB_CHK_ERR( error );
00159 
00160     std::cout << "Master (elements, parts) : (" << melems.size() << ", " << msets.size()
00161               << "), Slave (elements, parts) : (" << selems.size() << ", " << ssets.size() << ")" << std::endl;
00162 
00163     double master_radius = 1.0, slave_radius = 1.0;
00164     std::vector< double > mastercoords;
00165     Range masterverts, slaveverts;
00166     {
00167         error = mbCore->get_entities_by_dimension( masterfileset, 0, masterverts );MB_CHK_ERR( error );
00168         error = mbCore->get_entities_by_dimension( slavefileset, 0, slaveverts );MB_CHK_ERR( error );
00169     }
00170     if( use_spherical )
00171     {
00172         double points[6];
00173         EntityHandle mfrontback[2] = { masterverts[0], masterverts[masterverts.size() - 1] };
00174         error                      = mbCore->get_coords( &mfrontback[0], 2, points );MB_CHK_ERR( error );
00175         master_radius = 0.5 * ( std::sqrt( points[0] * points[0] + points[1] * points[1] + points[2] * points[2] ) +
00176                                 std::sqrt( points[3] * points[3] + points[4] * points[4] + points[5] * points[5] ) );
00177         EntityHandle sfrontback[2] = { slaveverts[0], slaveverts[slaveverts.size() - 1] };
00178         error                      = mbCore->get_coords( &sfrontback[0], 2, points );MB_CHK_ERR( error );
00179         slave_radius = 0.5 * ( std::sqrt( points[0] * points[0] + points[1] * points[1] + points[2] * points[2] ) +
00180                                std::sqrt( points[3] * points[3] + points[4] * points[4] + points[5] * points[5] ) );
00181         // Let us rescale both master and slave meshes to a unit sphere
00182         error = moab::IntxUtils::ScaleToRadius( mbCore, masterfileset, 1.0 );MB_CHK_ERR( error );
00183         error = moab::IntxUtils::ScaleToRadius( mbCore, slavefileset, 1.0 );MB_CHK_ERR( error );
00184     }
00185 
00186     try
00187     {
00188         std::map< int, moab::Range > spartvals;
00189         int npoints_notfound = 0;
00190         {
00191             FileOptions fopts( ( use_spherical ? "SPHERICAL" : "" ) );
00192             moab::EntityHandle tree_root = masterfileset;
00193 
00194             moab::AdaptiveKDTree tree( mbCore );
00195             error = tree.build_tree( melems, &tree_root, &fopts );
00196 
00197             // moab::BVHTree tree(mbCore);
00198             // error = tree.build_tree(melems, &tree_root);MB_CHK_ERR(error);
00199 
00200             for( size_t ie = 0; ie < selems.size(); ie++ )
00201             {
00202                 moab::EntityHandle selem, leaf;
00203                 double point[3];
00204                 selem = selems[ie];
00205 
00206                 // Get the element centroid to be queried
00207                 error = mbCore->get_coords( &selem, 1, point );MB_CHK_ERR( error );
00208 
00209                 std::vector< moab::EntityHandle > leaf_elems;
00210 
00211                 // Search for the closest source element in the master mesh corresponding
00212                 // to the target element centroid in the slave mesh
00213                 error = tree.point_search( point, leaf, treetolerance, btolerance );MB_CHK_ERR( error );
00214 
00215                 // We only care about the dimension that the user specified.
00216                 // MOAB partitions are ordered by elements anyway.
00217                 error = mbCore->get_entities_by_dimension( leaf, dimension, leaf_elems, true );MB_CHK_ERR( error );
00218 
00219                 if( leaf != 0 && leaf_elems.size() )
00220                 {
00221 
00222                     // Now get the master element centroids so that we can compute
00223                     // the minimum distance to the target point
00224                     std::vector< double > centroids( leaf_elems.size() * 3 );
00225                     error = mbCore->get_coords( &leaf_elems[0], leaf_elems.size(), &centroids[0] );MB_CHK_ERR( error );
00226 
00227                     if( !leaf_elems.size() )
00228                         std::cout << ie << ": "
00229                                   << " No leaf elements found." << std::endl;
00230                     // else std::cout << ie << " found " << leaf_elems.size() << " leaves for
00231                     // current point " << std::endl;
00232 
00233                     double dist = 1e10;
00234                     int pinelem = -1;
00235                     for( size_t il = 0; il < leaf_elems.size(); ++il )
00236                     {
00237                         const double* centroid = &centroids[il * 3];
00238                         const double locdist   = std::pow( point[0] - centroid[0], 2 ) +
00239                                                std::pow( point[1] - centroid[1], 2 ) +
00240                                                std::pow( point[2] - centroid[2], 2 );
00241                         if( locdist < dist && locdist < 1.0E-2 )
00242                         {
00243                             dist    = locdist;
00244                             pinelem = il;
00245 
00246 #ifdef VERBOSE
00247                             int gidMelem;
00248                             error = mbCore->tag_get_data( gidtag, &leaf_elems[il], 1, &gidMelem );MB_CHK_ERR( error );
00249                             std::cout << "\t Trial leaf " << il << " set " << gidMelem
00250                                       << " and part = " << get_map_value( mpartvals, gidMelem, -1 )
00251                                       << " with distance = " << locdist << std::endl;
00252 #endif
00253                         }
00254                     }
00255 
00256                     if( pinelem < 0 )
00257                     {
00258 #ifdef VERBOSE
00259                         std::cout << ie
00260                                   << ": [Error] - Could not find a minimum distance within the "
00261                                      "leaf nodes."
00262                                   << std::endl;
00263 #endif
00264                         npoints_notfound++;
00265                         spartvals[defaultpart].insert( selems[ie] );
00266                     }
00267                     else
00268                     {
00269                         int gidMelem;
00270                         error = mbCore->tag_get_data( gidtag, &leaf_elems[pinelem], 1, &gidMelem );MB_CHK_ERR( error );
00271 
00272                         int mpartid = get_map_value( mpartvals, gidMelem, -1 );
00273                         if( mpartid < 0 )
00274                             std::cout << "[WARNING]: Part number for element " << leaf_elems[pinelem]
00275                                       << " with global ID = " << gidMelem << " not found.\n";
00276 
00277 #ifdef VERBOSE
00278                         std::cout << "Adding element " << ie << " set " << mpartid << " with distance = " << dist
00279                                   << std::endl;
00280 #endif
00281                         spartvals[mpartid].insert( selems[ie] );
00282                     }
00283                 }
00284                 else
00285                 {
00286 #ifdef VERBOSE
00287                     std::cout << "[WARNING]: Adding element " << ie << " to default (part) set " << defaultpart
00288                               << std::endl;
00289 #endif
00290                     spartvals[defaultpart].insert( selems[ie] );
00291                 }
00292             }
00293             error = tree.reset_tree();MB_CHK_ERR( error );
00294         }
00295         if( npoints_notfound )
00296             std::cout << "Could not find " << npoints_notfound
00297                       << " points in the master mesh. Added to defaultpart set = " << defaultpart << std::endl;
00298 
00299         if( use_spherical )
00300         {
00301             error = moab::IntxUtils::ScaleToRadius( mbCore, slavefileset, slave_radius );MB_CHK_ERR( error );
00302         }
00303 
00304         error = mbCore->delete_entities( &masterfileset, 1 );MB_CHK_ERR( error );
00305         // Find parallel partition sets in the slave mesh - and delete it since we are going to
00306         // overwrite the sets
00307         if( !keepsparts )
00308         {
00309             std::cout << "Deleting " << ssets.size() << " sets in the slave mesh" << std::endl;
00310             error = mbCore->remove_entities( slavefileset, ssets );MB_CHK_ERR( error );
00311             ssets.clear();
00312         }
00313 
00314         size_t ntotslave_elems = 0, ntotslave_parts = 0;
00315         for( std::map< int, moab::Range >::iterator it = spartvals.begin(); it != spartvals.end(); ++it )
00316         {
00317             int partID = it->first;
00318             moab::EntityHandle pset;
00319             error = mbCore->create_meshset( moab::MESHSET_SET, pset );MB_CHK_ERR( error );
00320             error = mbCore->add_entities( pset, it->second );MB_CHK_ERR( error );
00321             error = mbCore->add_parent_child( slavefileset, pset );MB_CHK_ERR( error );
00322 
00323 #ifdef VERBOSE
00324             std::cout << "Slave Part " << partID << " has " << it->second.size() << " elements." << std::endl;
00325 #endif
00326             ntotslave_elems += it->second.size();
00327             ntotslave_parts++;
00328 
00329             if( keepsparts )
00330             {
00331                 error = mbCore->tag_set_data( sparttag, &pset, 1, &partID );MB_CHK_ERR( error );
00332             }
00333             else
00334             {
00335                 error = mbCore->tag_set_data( parttag, &pset, 1, &partID );MB_CHK_ERR( error );
00336             }
00337         }
00338         std::cout << "Slave mesh: given " << selems.size() << " elements, and assigned " << ntotslave_elems
00339                   << " elements to " << ntotslave_parts << " parts.\n";
00340         assert( ntotslave_elems == selems.size() );
00341 
00342         // mbCore->print_database();
00343 
00344         // Write the re-partitioned slave mesh to disk
00345         if( size == 1 )
00346         {
00347             error = mbCore->write_file( "slavemesh.vtk", "VTK", NULL, &slavefileset, 1 );MB_CHK_ERR( error );
00348         }
00349         error = mbCore->write_file( outfile.c_str(), NULL, write_options.c_str(), &slavefileset, 1 );MB_CHK_ERR( error );
00350         // error = mbCore->write_file(outfile.c_str(), NULL,
00351         // write_options.c_str());MB_CHK_ERR(error);
00352     }
00353     catch( std::exception& e )
00354     {
00355         std::cout << " exception caught during tree initialization " << e.what() << std::endl;
00356     }
00357     delete mbCore;
00358 
00359 #ifdef MOAB_HAVE_MPI
00360     MPI_Finalize();
00361 #endif
00362     exit( 0 );
00363 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines