Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
mbslavepart.cpp
Go to the documentation of this file.
1 //
2 // Usage:
3 // tools/mbslavepart -d 2 -m mpas/x1.2562.grid.h5m -s mpas/x1.10242.grid.h5m -o mpas_slave.h5m -e
4 // 1e-8 -b 1e-6 -O
5 //
6 #include <iostream>
7 #include <exception>
8 #include <cmath>
9 #include <vector>
10 #include <string>
11 
12 #include "moab/ProgOptions.hpp"
13 #include "moab/Core.hpp"
14 #include "moab/AdaptiveKDTree.hpp"
15 #include "moab/BVHTree.hpp"
16 
17 #include "moab/IntxMesh/IntxUtils.hpp"
18 
19 #ifdef MOAB_HAVE_MPI
20 #include "moab_mpi.h"
21 #include "moab/ParallelComm.hpp"
22 #include "MBParallelConventions.h"
23 #endif
24 
25 using namespace moab;
26 
27 // A function to get the non-default value from a std::map
28 template < typename K, typename V >
29 static V get_map_value( const std::map< K, V >& m, const K& key, const V& defval )
30 {
31  typename std::map< K, V >::const_iterator it = m.find( key );
32  if( it == m.end() )
33  {
34  return defval;
35  }
36  else
37  {
38  return it->second;
39  }
40 }
41 
42 int main( int argc, char* argv[] )
43 {
44  int nprocs = 1, dimension = 3;
45 #ifdef MOAB_HAVE_MPI
46  int proc_id = 0;
47  MPI_Init( &argc, &argv );
48  MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
49  MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
50 #endif
51 
52  int defaultpart = 0;
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;
57  ProgOptions opts;
58 
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",
64  &defaultpart );
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)",
69  &keepsparts );
70  opts.addOpt< void >( "spherical", "Hint that the meshes are defined on a spherical surface (Climate problems)",
71  &use_spherical );
72  opts.parseCommandLine( argc, argv );
73 
74  if( masterfile.empty() || slavefile.empty() )
75  {
76  opts.printHelp();
77 #ifdef MOAB_HAVE_MPI
78  MPI_Finalize();
79 #endif
80  exit( 1 );
81  }
82 
83  ErrorCode error;
84  Core* mbCore = new Core();
85 
86  // Set the read options for parallel file loading
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( "" );
91 
92  if( nprocs > 1 )
93  {
94  read_options = "PARALLEL=READ_PART;PARTITION=" + partition_set_name + ";PARALLEL_RESOLVE_SHARED_ENTS";
95  write_options = "PARALLEL=WRITE_PART";
96  }
97 
98  EntityHandle masterfileset, slavefileset;
99  error = mbCore->create_meshset( moab::MESHSET_TRACK_OWNER | moab::MESHSET_SET, masterfileset );
100  MB_CHK_ERR( error );
101  error = mbCore->create_meshset( moab::MESHSET_TRACK_OWNER | moab::MESHSET_SET, slavefileset );
102  MB_CHK_ERR( error );
103 
104  // Load file
105  error = mbCore->load_file( masterfile.c_str(), &masterfileset, read_options.c_str() );
106  MB_CHK_ERR( error );
107  error = mbCore->load_file( slavefile.c_str(), &slavefileset, read_options.c_str() );
108  MB_CHK_ERR( error );
109  // if (error != MB_SUCCESS && size > 1)
110  // {
111  // std::string newread_options = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS";
112  // error = mbCore->load_file(slavefile.c_str(), &slavefileset, newread_options.c_str());
113  // }
114  // else MB_CHK_ERR(error);
115 
116  Tag gidtag = 0, parttag = 0, sparttag = 0;
117  int dum_id = -1;
118  error = mbCore->tag_get_handle( partition_set_name.c_str(), 1, MB_TYPE_INTEGER, parttag,
119  MB_TAG_SPARSE | MB_TAG_CREAT, &dum_id );
120  MB_CHK_ERR( error );
121  gidtag = mbCore->globalId_tag();
122  if( keepsparts )
123  {
124  error = mbCore->tag_get_handle( std::string( partition_set_name + "_SLAVE" ).c_str(), 1, MB_TYPE_INTEGER,
125  sparttag, MB_TAG_CREAT | MB_TAG_SPARSE, &dum_id );
126  MB_CHK_ERR( error );
127  }
128 
129  Range melems, msets, selems, ssets;
130 
131  // Get the partition sets on the master mesh
132  std::map< int, int > mpartvals;
133  error = mbCore->get_entities_by_type_and_tag( masterfileset, MBENTITYSET, &parttag, NULL, 1, msets,
134  moab::Interface::UNION, true );
135  MB_CHK_ERR( error );
136  if( msets.size() == 0 )
137  {
138  std::cout << "No partition sets found in the master mesh. Quitting..." << std::endl;
139  exit( 1 );
140  }
141 
142  for( unsigned i = 0; i < msets.size(); ++i )
143  {
144  EntityHandle mset = msets[i];
145 
146  moab::Range msetelems;
147  error = mbCore->get_entities_by_dimension( mset, dimension, msetelems );
148  MB_CHK_ERR( error );
149  melems.merge( msetelems );
150 
151  int partID;
152  error = mbCore->tag_get_data( parttag, &mset, 1, &partID );
153  MB_CHK_ERR( error );
154 
155  // Get the global ID and use that as the indicator
156  std::vector< int > gidMelems( msetelems.size() );
157  error = mbCore->tag_get_data( gidtag, msetelems, gidMelems.data() );
158  MB_CHK_ERR( error );
159 
160  for( unsigned j = 0; j < msetelems.size(); ++j )
161  mpartvals[gidMelems[j]] = partID;
162  // mpartvals[msetelems[j]]=partID;
163 #ifdef VERBOSE
164  std::cout << "Part " << partID << " has " << msetelems.size() << " elements." << std::endl;
165 #endif
166  }
167 
168  // Get information about the slave file set
169  error = mbCore->get_entities_by_type_and_tag( slavefileset, MBENTITYSET, &parttag, NULL, 1, ssets,
170  moab::Interface::UNION );
171  MB_CHK_ERR( error );
172  // TODO: expand and add other dimensional elements
173  error = mbCore->get_entities_by_dimension( slavefileset, dimension, selems );
174  MB_CHK_ERR( error );
175 
176  std::cout << "Master (elements, parts) : (" << melems.size() << ", " << msets.size()
177  << "), Slave (elements, parts) : (" << selems.size() << ", " << ssets.size() << ")" << std::endl;
178 
179  // double master_radius = 1.0;
180  double slave_radius = 1.0;
181  std::vector< double > mastercoords;
182  Range masterverts, slaveverts;
183  {
184  error = mbCore->get_entities_by_dimension( masterfileset, 0, masterverts );
185  MB_CHK_ERR( error );
186  error = mbCore->get_entities_by_dimension( slavefileset, 0, slaveverts );
187  MB_CHK_ERR( error );
188  }
189  if( use_spherical )
190  {
191  double points[6];
192  EntityHandle mfrontback[2] = { masterverts[0], masterverts[masterverts.size() - 1] };
193  error = mbCore->get_coords( &mfrontback[0], 2, points );
194  MB_CHK_ERR( error );
195  // master_radius = 0.5 * ( std::sqrt( points[0] * points[0] + points[1] * points[1] + points[2] * points[2] ) +
196  // std::sqrt( points[3] * points[3] + points[4] * points[4] + points[5] * points[5] ) );
197  EntityHandle sfrontback[2] = { slaveverts[0], slaveverts[slaveverts.size() - 1] };
198  error = mbCore->get_coords( &sfrontback[0], 2, points );
199  MB_CHK_ERR( error );
200  slave_radius = 0.5 * ( std::sqrt( points[0] * points[0] + points[1] * points[1] + points[2] * points[2] ) +
201  std::sqrt( points[3] * points[3] + points[4] * points[4] + points[5] * points[5] ) );
202  // Let us rescale both master and slave meshes to a unit sphere
203  error = moab::IntxUtils::ScaleToRadius( mbCore, masterfileset, 1.0 );
204  MB_CHK_ERR( error );
205  error = moab::IntxUtils::ScaleToRadius( mbCore, slavefileset, 1.0 );
206  MB_CHK_ERR( error );
207  }
208 
209  try
210  {
211  std::map< int, moab::Range > spartvals;
212  int npoints_notfound = 0;
213  {
214  FileOptions fopts( ( use_spherical ? "SPHERICAL" : "" ) );
215  moab::EntityHandle tree_root = masterfileset;
216 
217  moab::AdaptiveKDTree tree( mbCore );
218  error = tree.build_tree( melems, &tree_root, &fopts );
219 
220  // moab::BVHTree tree(mbCore);
221  // error = tree.build_tree(melems, &tree_root);MB_CHK_ERR(error);
222 
223  for( size_t ie = 0; ie < selems.size(); ie++ )
224  {
225  moab::EntityHandle selem, leaf;
226  double point[3];
227  selem = selems[ie];
228 
229  // Get the element centroid to be queried
230  error = mbCore->get_coords( &selem, 1, point );
231  MB_CHK_ERR( error );
232 
233  std::vector< moab::EntityHandle > leaf_elems;
234 
235  // Search for the closest source element in the master mesh corresponding
236  // to the target element centroid in the slave mesh
237  error = tree.point_search( point, leaf, treetolerance, btolerance );
238  MB_CHK_ERR( error );
239 
240  // We only care about the dimension that the user specified.
241  // MOAB partitions are ordered by elements anyway.
242  error = mbCore->get_entities_by_dimension( leaf, dimension, leaf_elems, true );
243  MB_CHK_ERR( error );
244 
245  if( leaf != 0 && leaf_elems.size() )
246  {
247 
248  // Now get the master element centroids so that we can compute
249  // the minimum distance to the target point
250  std::vector< double > centroids( leaf_elems.size() * 3 );
251  error = mbCore->get_coords( &leaf_elems[0], leaf_elems.size(), &centroids[0] );
252  MB_CHK_ERR( error );
253 
254  if( !leaf_elems.size() )
255  std::cout << ie << ": "
256  << " No leaf elements found." << std::endl;
257  // else std::cout << ie << " found " << leaf_elems.size() << " leaves for
258  // current point " << std::endl;
259 
260  double dist = 1e10;
261  int pinelem = -1;
262  for( size_t il = 0; il < leaf_elems.size(); ++il )
263  {
264  const double* centroid = &centroids[il * 3];
265  const double locdist = std::pow( point[0] - centroid[0], 2 ) +
266  std::pow( point[1] - centroid[1], 2 ) +
267  std::pow( point[2] - centroid[2], 2 );
268  if( locdist < dist && locdist < 1.0E-2 )
269  {
270  dist = locdist;
271  pinelem = il;
272 
273 #ifdef VERBOSE
274  int gidMelem;
275  error = mbCore->tag_get_data( gidtag, &leaf_elems[il], 1, &gidMelem );
276  MB_CHK_ERR( error );
277  std::cout << "\t Trial leaf " << il << " set " << gidMelem
278  << " and part = " << get_map_value( mpartvals, gidMelem, -1 )
279  << " with distance = " << locdist << std::endl;
280 #endif
281  }
282  }
283 
284  if( pinelem < 0 )
285  {
286 #ifdef VERBOSE
287  std::cout << ie
288  << ": [Error] - Could not find a minimum distance within the "
289  "leaf nodes."
290  << std::endl;
291 #endif
292  npoints_notfound++;
293  spartvals[defaultpart].insert( selems[ie] );
294  }
295  else
296  {
297  int gidMelem;
298  error = mbCore->tag_get_data( gidtag, &leaf_elems[pinelem], 1, &gidMelem );
299  MB_CHK_ERR( error );
300 
301  int mpartid = get_map_value( mpartvals, gidMelem, -1 );
302  if( mpartid < 0 )
303  std::cout << "[WARNING]: Part number for element " << leaf_elems[pinelem]
304  << " with global ID = " << gidMelem << " not found.\n";
305 
306 #ifdef VERBOSE
307  std::cout << "Adding element " << ie << " set " << mpartid << " with distance = " << dist
308  << std::endl;
309 #endif
310  spartvals[mpartid].insert( selems[ie] );
311  }
312  }
313  else
314  {
315 #ifdef VERBOSE
316  std::cout << "[WARNING]: Adding element " << ie << " to default (part) set " << defaultpart
317  << std::endl;
318 #endif
319  spartvals[defaultpart].insert( selems[ie] );
320  }
321  }
322  error = tree.reset_tree();
323  MB_CHK_ERR( error );
324  }
325  if( npoints_notfound )
326  std::cout << "Could not find " << npoints_notfound
327  << " points in the master mesh. Added to defaultpart set = " << defaultpart << std::endl;
328 
329  if( use_spherical )
330  {
331  error = moab::IntxUtils::ScaleToRadius( mbCore, slavefileset, slave_radius );
332  MB_CHK_ERR( error );
333  }
334 
335  error = mbCore->delete_entities( &masterfileset, 1 );
336  MB_CHK_ERR( error );
337  // Find parallel partition sets in the slave mesh - and delete it since we are going to
338  // overwrite the sets
339  if( !keepsparts )
340  {
341  std::cout << "Deleting " << ssets.size() << " sets in the slave mesh" << std::endl;
342  error = mbCore->remove_entities( slavefileset, ssets );
343  MB_CHK_ERR( error );
344  ssets.clear();
345  }
346 
347  size_t ntotslave_elems = 0, ntotslave_parts = 0;
348  for( std::map< int, moab::Range >::iterator it = spartvals.begin(); it != spartvals.end(); ++it )
349  {
350  int partID = it->first;
351  moab::EntityHandle pset;
352  error = mbCore->create_meshset( moab::MESHSET_SET, pset );
353  MB_CHK_ERR( error );
354  error = mbCore->add_entities( pset, it->second );
355  MB_CHK_ERR( error );
356  error = mbCore->add_parent_child( slavefileset, pset );
357  MB_CHK_ERR( error );
358 
359 #ifdef VERBOSE
360  std::cout << "Slave Part " << partID << " has " << it->second.size() << " elements." << std::endl;
361 #endif
362  ntotslave_elems += it->second.size();
363  ntotslave_parts++;
364 
365  if( keepsparts )
366  {
367  error = mbCore->tag_set_data( sparttag, &pset, 1, &partID );
368  MB_CHK_ERR( error );
369  }
370  else
371  {
372  error = mbCore->tag_set_data( parttag, &pset, 1, &partID );
373  MB_CHK_ERR( error );
374  }
375  }
376  std::cout << "Slave mesh: given " << selems.size() << " elements, and assigned " << ntotslave_elems
377  << " elements to " << ntotslave_parts << " parts.\n";
378  assert( ntotslave_elems == selems.size() );
379 
380  // mbCore->print_database();
381 
382  // Write the re-partitioned slave mesh to disk
383  if( nprocs == 1 )
384  {
385  error = mbCore->write_file( "slavemesh.vtk", "VTK", NULL, &slavefileset, 1 );
386  MB_CHK_ERR( error );
387  }
388  error = mbCore->write_file( outfile.c_str(), NULL, write_options.c_str(), &slavefileset, 1 );
389  MB_CHK_ERR( error );
390  // error = mbCore->write_file(outfile.c_str(), NULL,
391  // write_options.c_str());MB_CHK_ERR(error);
392  }
393  catch( std::exception& e )
394  {
395  std::cout << " exception caught during tree initialization " << e.what() << std::endl;
396  }
397  delete mbCore;
398 
399 #ifdef MOAB_HAVE_MPI
400  MPI_Finalize();
401 #endif
402  exit( 0 );
403 }