|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
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(), ¢roids[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 = ¢roids[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 }
1.7.6.1 
