28 template <
typename K,
typename V >
29 static V
get_map_value(
const std::map< K, V >& m,
const K& key,
const V& defval )
31 typename std::map< K, V >::const_iterator it = m.find( key );
42 int main(
int argc,
char* argv[] )
44 int nprocs = 1, dimension = 3;
47 MPI_Init( &argc, &argv );
48 MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
49 MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
53 double tolerance = 1e-2, treetolerance = 1e-13, btolerance = 1e-4;
54 std::string masterfile, slavefile, outfile(
"slavemesh.h5m" );
55 bool keepsparts =
false;
56 bool use_spherical =
false;
59 opts.
addOpt< std::string >(
"master,m",
"Master mesh filename", &masterfile );
60 opts.
addOpt< std::string >(
"slave,s",
"Slave mesh filename", &slavefile );
61 opts.
addOpt< std::string >(
"output,o",
"Output partitioned mesh filename", &outfile );
62 opts.
addOpt<
int >(
"dim,d",
"Dimension of entities to use for partitioning", &dimension );
63 opts.
addOpt<
int >(
"defaultpart,p",
"Default partition number if target element is not found on source grid",
65 opts.
addOpt<
double >(
"eps,e",
"Tolerance for the point search", &
tolerance );
66 opts.
addOpt<
double >(
"beps,b",
"Tolerance for the bounding box search", &btolerance );
67 opts.
addOpt<
void >(
"keep,K",
68 "Keep the existing partitions in the slave mesh (use PARALLEL_PARTITION_SLAVE instead)",
70 opts.
addOpt<
void >(
"spherical",
"Hint that the meshes are defined on a spherical surface (Climate problems)",
74 if( masterfile.empty() || slavefile.empty() )
87 const std::string partition_set_name =
"PARALLEL_PARTITION";
88 const std::string global_id_name =
"GLOBAL_ID";
89 std::vector< std::string > read_opts, write_opts;
90 std::string read_options(
"" ), write_options(
"" );
94 read_options =
"PARALLEL=READ_PART;PARTITION=" + partition_set_name +
";PARALLEL_RESOLVE_SHARED_ENTS";
95 write_options =
"PARALLEL=WRITE_PART";
112 Tag gidtag = 0, parttag = 0, sparttag = 0;
123 Range melems, msets, selems, ssets;
126 std::map< int, int > mpartvals;
129 if( msets.
size() == 0 )
131 std::cout <<
"No partition sets found in the master mesh. Quitting..." << std::endl;
135 for(
unsigned i = 0; i < msets.
size(); ++i )
141 melems.
merge( msetelems );
147 std::vector< int > gidMelems( msetelems.
size() );
150 for(
unsigned j = 0; j < msetelems.
size(); ++j )
151 mpartvals[gidMelems[j]] = partID;
154 std::cout <<
"Part " << partID <<
" has " << msetelems.
size() <<
" elements." << std::endl;
164 std::cout <<
"Master (elements, parts) : (" << melems.
size() <<
", " << msets.
size()
165 <<
"), Slave (elements, parts) : (" << selems.
size() <<
", " << ssets.
size() <<
")" << std::endl;
167 double master_radius = 1.0, slave_radius = 1.0;
168 std::vector< double > mastercoords;
169 Range masterverts, slaveverts;
177 EntityHandle mfrontback[2] = { masterverts[0], masterverts[masterverts.
size() - 1] };
179 master_radius = 0.5 * ( std::sqrt( points[0] * points[0] + points[1] * points[1] + points[2] * points[2] ) +
180 std::sqrt( points[3] * points[3] + points[4] * points[4] + points[5] * points[5] ) );
181 EntityHandle sfrontback[2] = { slaveverts[0], slaveverts[slaveverts.
size() - 1] };
183 slave_radius = 0.5 * ( std::sqrt( points[0] * points[0] + points[1] * points[1] + points[2] * points[2] ) +
184 std::sqrt( points[3] * points[3] + points[4] * points[4] + points[5] * points[5] ) );
192 std::map< int, moab::Range > spartvals;
193 int npoints_notfound = 0;
195 FileOptions fopts( ( use_spherical ?
"SPHERICAL" :
"" ) );
204 for(
size_t ie = 0; ie < selems.
size(); ie++ )
213 std::vector< moab::EntityHandle > leaf_elems;
223 if( leaf != 0 && leaf_elems.size() )
228 std::vector< double > centroids( leaf_elems.size() * 3 );
231 if( !leaf_elems.size() )
232 std::cout << ie <<
": "
233 <<
" No leaf elements found." << std::endl;
239 for(
size_t il = 0; il < leaf_elems.size(); ++il )
241 const double* centroid = ¢roids[il * 3];
242 const double locdist = std::pow( point[0] - centroid[0], 2 ) +
243 std::pow( point[1] - centroid[1], 2 ) +
244 std::pow( point[2] - centroid[2], 2 );
245 if( locdist < dist && locdist < 1.0E-2 )
253 std::cout <<
"\t Trial leaf " << il <<
" set " << gidMelem
254 <<
" and part = " <<
get_map_value( mpartvals, gidMelem, -1 )
255 <<
" with distance = " << locdist << std::endl;
264 <<
": [Error] - Could not find a minimum distance within the "
269 spartvals[defaultpart].insert( selems[ie] );
278 std::cout <<
"[WARNING]: Part number for element " << leaf_elems[pinelem]
279 <<
" with global ID = " << gidMelem <<
" not found.\n";
282 std::cout <<
"Adding element " << ie <<
" set " << mpartid <<
" with distance = " << dist
285 spartvals[mpartid].insert( selems[ie] );
291 std::cout <<
"[WARNING]: Adding element " << ie <<
" to default (part) set " << defaultpart
294 spartvals[defaultpart].insert( selems[ie] );
299 if( npoints_notfound )
300 std::cout <<
"Could not find " << npoints_notfound
301 <<
" points in the master mesh. Added to defaultpart set = " << defaultpart << std::endl;
313 std::cout <<
"Deleting " << ssets.
size() <<
" sets in the slave mesh" << std::endl;
318 size_t ntotslave_elems = 0, ntotslave_parts = 0;
319 for( std::map< int, moab::Range >::iterator it = spartvals.begin(); it != spartvals.end(); ++it )
321 int partID = it->first;
328 std::cout <<
"Slave Part " << partID <<
" has " << it->second.size() <<
" elements." << std::endl;
330 ntotslave_elems += it->second.size();
342 std::cout <<
"Slave mesh: given " << selems.
size() <<
" elements, and assigned " << ntotslave_elems
343 <<
" elements to " << ntotslave_parts <<
" parts.\n";
344 assert( ntotslave_elems == selems.
size() );
357 catch( std::exception& e )
359 std::cout <<
" exception caught during tree initialization " << e.what() << std::endl;