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