Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
mbcoupler_test.cpp
Go to the documentation of this file.
1 // MOAB includes
2 #include "moab/ParallelComm.hpp"
4 #include "moab/Core.hpp"
5 #include "Coupler.hpp"
6 #include "moab_mpi.h"
7 #include "ElemUtil.hpp"
8 #include "TestUtil.hpp"
9 
10 // C++ includes
11 #include <iostream>
12 #include <iomanip>
13 #include <sstream>
14 #include <cassert>
15 
16 using namespace moab;
17 
18 /*
19  Sample usages:
20  1) P_0 interpolation: ./mbcoupler_test -meshes <src_mesh> <target_mesh> -itag <interp_tag> -meth
21  0 -outfile <output> 2) P_1 interpolation: ./mbcoupler_test -meshes <src_mesh> <target_mesh> -itag
22  <interp_tag> -meth 1 -outfile <output> 3) P_1 interpolation with epsilon control: ./mbcoupler_test
23  -meshes <src_mesh> <target_mesh> -itag <interp_tag> -meth 1 -eps <tolerance for locating points;
24  say 0.01> -outfile <output> 3) P_0 interpolation with global normalization: ./mbcoupler_test
25  -meshes <src_mesh> <target_mesh> -itag <interp_tag> -meth 0 -gnorm <gnorm_tag_name> -outfile
26  <output> 4) P_1 interpolation with subset normalization: ./mbcoupler_test -meshes <src_mesh>
27  <target_mesh> -itag <interp_tag> -meth 1 -ssnorm <snorm_tag_name> <snorm_criteria: MATERIAL_SET>
28  -outfile <output> 5) P_1 interpolation for meshes on sphere ./mbcoupler_test -meshes <src_mesh>
29  <target_mesh> -itag <interp_tag> -meth 4 -outfile <output>
30 */
31 
32 // Print usage
34 {
35  std::cerr << "Usage: ";
36  std::cerr << "mbcoupler_test -meshes <source_mesh> <target_mesh> -itag <interp_tag> [-gnorm "
37  "<gnorm_tag>] [-ssnorm <ssnorm_tag> <ssnorm_selection>] [-ropts <roptions>] "
38  "[-outfile <out_file> [-wopts <woptions>]] [-dbgout [<dbg_file>]]"
39  << std::endl;
40  std::cerr << " -meshes" << std::endl;
41  std::cerr << " Read in mesh files <source_mesh> and <target_mesh>." << std::endl;
42  std::cerr << " -itag" << std::endl;
43  std::cerr << " Interpolate tag <interp_tag> from source mesh to target mesh." << std::endl;
44  std::cerr << " -gnorm" << std::endl;
45  std::cerr << " Normalize the value of tag <gnorm_tag> over then entire mesh and save to" << std::endl;
46  std::cerr << " tag \"<gnorm_tag>_normf\" on the mesh set. Do this for all meshes." << std::endl;
47  std::cerr << " -ssnorm" << std::endl;
48  std::cerr << " Normalize the value of tag <ssnorm_tag> over subsets of a mesh and save to" << std::endl;
49  std::cerr << " tag \"<ssnorm_tag>_normf\" on the Entity Set for each subset. Subsets "
50  "are selected"
51  << std::endl;
52  std::cerr << " using criteria in <ssnorm_selection>. Do this for all meshes." << std::endl;
53  std::cerr << " -ropts" << std::endl;
54  std::cerr << " Read in the mesh files using options in <roptions>." << std::endl;
55  std::cerr << " -outfile" << std::endl;
56  std::cerr << " Write out target mesh to <out_file>." << std::endl;
57  std::cerr << " -wopts" << std::endl;
58  std::cerr << " Write out mesh files using options in <woptions>." << std::endl;
59  std::cerr << " -dbgout" << std::endl;
60  std::cerr << " Write stdout and stderr streams to the file \'<dbg_file>.txt\'." << std::endl;
61  std::cerr << " -eps" << std::endl;
62  std::cerr << " epsilon" << std::endl;
63  std::cerr << " -meth <method> (0=CONSTANT, 1=LINEAR_FE, 2=QUADRATIC_FE, 3=SPECTRAL, 4=SPHERICAL)" << std::endl;
64 }
65 
66 #ifdef MOAB_HAVE_HDF5
67 
68 ErrorCode get_file_options( int argc,
69  char** argv,
70  int nprocs,
71  int rank,
72  std::vector< std::string >& meshFiles,
73  Coupler::Method& method,
74  std::string& interpTag,
75  std::string& gNormTag,
76  std::string& ssNormTag,
77  std::vector< const char* >& ssTagNames,
78  std::vector< const char* >& ssTagValues,
79  std::string& readOpts,
80  std::string& outFile,
81  std::string& writeOpts,
82  std::string& dbgFile,
83  bool& help,
84  double& epsilon );
85 
86 ErrorCode report_iface_ents( Interface* mbImpl, std::vector< ParallelComm* >& pcs, bool print_results );
87 
88 ErrorCode test_interpolation( Interface* mbImpl,
89  Coupler::Method method,
90  std::string& interpTag,
91  std::string& gNormTag,
92  std::string& ssNormTag,
93  std::vector< const char* >& ssTagNames,
94  std::vector< const char* >& ssTagValues,
95  EntityHandle* roots,
96  std::vector< ParallelComm* >& pcs,
97  double& instant_time,
98  double& pointloc_time,
99  double& interp_time,
100  double& gnorm_time,
101  double& ssnorm_time,
102  double& toler );
103 
104 void reduceMax( double& v )
105 {
106  double buf;
107 
108  MPI_Barrier( MPI_COMM_WORLD );
109  MPI_Allreduce( &v, &buf, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
110 
111  v = buf;
112 }
113 
114 int main( int argc, char** argv )
115 {
116  // Need to init MPI first, to tell how many procs and rank
117  int ierr = MPI_Init( &argc, &argv );
118  assert( MPI_SUCCESS == ierr );
119 
120  std::vector< const char* > ssTagNames, ssTagValues;
121  std::vector< std::string > meshFiles;
122  std::string interpTag, gNormTag, ssNormTag, readOpts, outFile, writeOpts, dbgFile;
124 
125  ErrorCode result;
126  bool help = false;
127  double toler = 5.e-10;
128  int nprocs, rank;
129  ierr = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
130  assert( MPI_SUCCESS == ierr );
131  ierr = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
132  assert( MPI_SUCCESS == ierr );
133 
134  result = get_file_options( argc, argv, nprocs, rank, meshFiles, method, interpTag, gNormTag, ssNormTag, ssTagNames,
135  ssTagValues, readOpts, outFile, writeOpts, dbgFile, help, toler );
136 
137  if( result != MB_SUCCESS || help )
138  {
139  print_usage();
140  ierr = MPI_Finalize();
141  return ierr;
142  }
143 
144  // Redirect stdout and stderr if dbgFile is not null
145  if( !dbgFile.empty() )
146  {
147  std::stringstream dfname;
148  dfname << dbgFile << rank << ".txt";
149  if( !std::freopen( dfname.str().c_str(), "a", stdout ) ) return -1;
150  if( !std::freopen( dfname.str().c_str(), "a", stderr ) ) return -1;
151  }
152 
153  // Create MOAB instance based on that
154  Interface* mbImpl = new( std::nothrow ) Core();
155  if( NULL == mbImpl ) return 1;
156 
157  // Read in mesh(es)
158  std::vector< ParallelComm* > pcs( meshFiles.size() );
159 
160  // Create root sets for each mesh using the iMesh API. Then pass these
161  // to the load_file functions to be populated.
162  EntityHandle* roots = (EntityHandle*)malloc( sizeof( EntityHandle ) * meshFiles.size() );
163 
164  for( unsigned int i = 0; i < meshFiles.size(); i++ )
165  {
166  pcs[i] = new ParallelComm( mbImpl, MPI_COMM_WORLD );
167  int index = pcs[i]->get_id();
168  std::string newReadopts;
169  std::ostringstream extraOpt;
170  extraOpt << ";PARALLEL_COMM=" << index;
171  newReadopts = readOpts + extraOpt.str();
172 
173  result = mbImpl->create_meshset( MESHSET_SET, roots[i] );MB_CHK_ERR( result );
174 
175  result = mbImpl->load_file( meshFiles[i].c_str(), &roots[i], newReadopts.c_str() );MB_CHK_ERR( result );
176  // result = rps[i]->load_file(meshFiles[i].c_str(), &roots[i],
177  // FileOptions(readOpts.c_str()));MB_CHK_ERR(result);
178  }
179 
180  result = report_iface_ents( mbImpl, pcs, true );MB_CHK_ERR( result );
181 
182  double instant_time = 0.0, pointloc_time = 0.0, interp_time = 0.0, gnorm_time = 0.0, ssnorm_time = 0.0;
183  // Test interpolation and global normalization and subset normalization
184 
185  result = test_interpolation( mbImpl, method, interpTag, gNormTag, ssNormTag, ssTagNames, ssTagValues, roots, pcs,
186  instant_time, pointloc_time, interp_time, gnorm_time, ssnorm_time, toler );MB_CHK_ERR( result );
187 
188  reduceMax( instant_time );
189  reduceMax( pointloc_time );
190  reduceMax( interp_time );
191 
192  if( 0 == rank )
193  printf( "\nMax time : %g %g %g (inst loc interp -- %d procs)\n", instant_time, pointloc_time, interp_time,
194  nprocs );
195 
196  // Output mesh
197  if( !outFile.empty() )
198  {
199  Range partSets;
200  // Only save the target mesh
201  partSets.insert( (EntityHandle)roots[1] );
202  std::string newwriteOpts;
203  std::ostringstream extraOpt;
204  if( nprocs > 1 ) extraOpt << ";PARALLEL_COMM=" << 1;
205  newwriteOpts = writeOpts + extraOpt.str();
206  result = mbImpl->write_file( outFile.c_str(), NULL, newwriteOpts.c_str(), partSets );MB_CHK_ERR( result );
207  if( 0 == rank )
208  {
209  std::cout << "Wrote " << outFile << std::endl;
210  std::cout << "mbcoupler_test complete." << std::endl;
211  }
212  }
213 
214  free( roots );
215 
216  for( unsigned int i = 0; i < meshFiles.size(); i++ )
217  delete pcs[i];
218 
219  delete mbImpl;
220 
221  ierr = MPI_Finalize();
222  assert( MPI_SUCCESS == ierr );
223  return 0;
224 }
225 
226 ErrorCode report_iface_ents( Interface* mbImpl, std::vector< ParallelComm* >& pcs, const bool print_results )
227 {
228  Range iface_ents[6];
229  ErrorCode result = MB_SUCCESS, tmp_result;
230 
231  // Now figure out which vertices are shared
232  for( unsigned int p = 0; p < pcs.size(); p++ )
233  {
234  for( int i = 0; i < 4; i++ )
235  {
236  tmp_result = pcs[p]->get_iface_entities( -1, i, iface_ents[i] );
237 
238  if( MB_SUCCESS != tmp_result )
239  {
240  std::cerr << "get_iface_entities returned error on proc " << pcs[p]->proc_config().proc_rank()
241  << "; message: " << std::endl;
242  std::string last_error;
243  result = mbImpl->get_last_error( last_error );
244  if( last_error.empty() )
245  std::cerr << "(none)" << std::endl;
246  else
247  std::cerr << last_error << std::endl;
248  result = tmp_result;
249  }
250  if( 0 != i ) iface_ents[4].merge( iface_ents[i] );
251  }
252  }
253 
254  // Report # iface entities
255  result = mbImpl->get_adjacencies( iface_ents[4], 0, false, iface_ents[5], Interface::UNION );MB_CHK_ERR( result );
256 
257  int rank;
258  MPI_Comm_rank( MPI_COMM_WORLD, &rank );
259 
260  if( print_results || iface_ents[0].size() != iface_ents[5].size() )
261  {
262  std::cerr << "Proc " << rank << " iface entities: " << std::endl;
263  for( int i = 0; i < 4; i++ )
264  std::cerr << " " << iface_ents[i].size() << " " << i << "d iface entities." << std::endl;
265  std::cerr << " (" << iface_ents[5].size() << " verts adj to other iface ents)" << std::endl;
266  }
267 
268  return result;
269 }
270 
271 // Check first character for a '-'.
272 // Return true if one is found. False otherwise.
273 bool check_for_flag( const char* str )
274 {
275  if( '-' == str[0] )
276  return true;
277  else
278  return false;
279 }
280 
281 // get_file_options() function with added possibilities for mbcoupler_test.
282 ErrorCode get_file_options( int argc,
283  char** argv,
284  int nprocs,
285  int rank,
286  std::vector< std::string >& meshFiles,
287  Coupler::Method& method,
288  std::string& interpTag,
289  std::string& gNormTag,
290  std::string& ssNormTag,
291  std::vector< const char* >& ssTagNames,
292  std::vector< const char* >& ssTagValues,
293  std::string& readOpts,
294  std::string& outFile,
295  std::string& writeOpts,
296  std::string& dbgFile,
297  bool& help,
298  double& epsilon )
299 {
300  // Initialize some of the outputs to null values indicating not present
301  // in the argument list.
302  gNormTag = "";
303  ssNormTag = "";
304  readOpts = ( nprocs > 1 ? "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;"
305  "PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=3.0.1;CPUTIME"
306  : "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;"
307  "PARALLEL_RESOLVE_SHARED_ENTS;CPUTIME" );
308  outFile = "";
309  writeOpts = ( nprocs > 1 ? "PARALLEL=WRITE_PART;CPUTIME" : "" );
310  dbgFile = "";
311  std::string defaultDbgFile = argv[0]; // The executable name will be the default debug output file.
312 
313  // These will indicate if we've gotten our required parameters at the end of parsing.
314  bool haveMeshes = false;
315  bool haveInterpTag = false;
316 
317  // Loop over the values in argv pulling out an parsing each one
318  int npos = 1;
319 
320  if( argc > 1 && argv[1] == std::string( "-h" ) )
321  {
322  help = true;
323  return MB_SUCCESS;
324  }
325 
326  while( npos < argc )
327  {
328  if( argv[npos] == std::string( "-meshes" ) )
329  {
330  // Parse out the mesh filenames
331  npos++;
332  int numFiles = 2;
333  meshFiles.resize( numFiles );
334  for( int i = 0; i < numFiles; i++ )
335  {
336  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
337  meshFiles[i] = argv[npos++];
338  else
339  {
340  std::cerr << " ERROR - missing correct number of mesh filenames" << std::endl;
341  return MB_FAILURE;
342  }
343  }
344 
345  haveMeshes = true;
346  }
347  else if( argv[npos] == std::string( "-itag" ) )
348  {
349  // Parse out the interpolation tag
350  npos++;
351  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
352  interpTag = argv[npos++];
353  else
354  {
355  std::cerr << " ERROR - missing <interp_tag>" << std::endl;
356  return MB_FAILURE;
357  }
358 
359  haveInterpTag = true;
360  }
361  else if( argv[npos] == std::string( "-meth" ) )
362  {
363  // Parse out the interpolation tag
364  npos++;
365  if( argv[npos][0] == '0' )
366  method = Coupler::CONSTANT;
367  else if( argv[npos][0] == '1' )
368  method = Coupler::LINEAR_FE;
369  else if( argv[npos][0] == '2' )
370  method = Coupler::QUADRATIC_FE;
371  else if( argv[npos][0] == '3' )
372  method = Coupler::SPECTRAL;
373  else if( argv[npos][0] == '4' )
374  method = Coupler::SPHERICAL;
375  else
376  {
377  std::cerr << " ERROR - unrecognized method number " << method << std::endl;
378  return MB_FAILURE;
379  }
380  npos++;
381  }
382  else if( argv[npos] == std::string( "-eps" ) )
383  {
384  // Parse out the tolerance
385  npos++;
386  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
387  epsilon = atof( argv[npos++] );
388  else
389  {
390  std::cerr << " ERROR - missing <epsilon>" << std::endl;
391  return MB_FAILURE;
392  }
393  }
394  else if( argv[npos] == std::string( "-gnorm" ) )
395  {
396  // Parse out the global normalization tag
397  npos++;
398  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
399  gNormTag = argv[npos++];
400  else
401  {
402  std::cerr << " ERROR - missing <gnorm_tag>" << std::endl;
403  return MB_FAILURE;
404  }
405  }
406  else if( argv[npos] == std::string( "-ssnorm" ) )
407  {
408  // Parse out the subset normalization tag and selection criteria
409  npos++;
410  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
411  ssNormTag = argv[npos++];
412  else
413  {
414  std::cerr << " ERROR - missing <ssnorm_tag>" << std::endl;
415  return MB_FAILURE;
416  }
417 
418  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
419  {
420  char* opts = argv[npos++];
421  char sep1[1] = { ';' };
422  char sep2[1] = { '=' };
423  bool end_vals_seen = false;
424  std::vector< char* > tmpTagOpts;
425 
426  // First get the options
427  for( char* i = strtok( opts, sep1 ); i; i = strtok( 0, sep1 ) )
428  tmpTagOpts.push_back( i );
429 
430  // Parse out the name and val or just name.
431  for( unsigned int j = 0; j < tmpTagOpts.size(); j++ )
432  {
433  char* e = strtok( tmpTagOpts[j], sep2 );
434  ssTagNames.push_back( e );
435  e = strtok( 0, sep2 );
436  if( e != NULL )
437  {
438  // We have a value
439  if( end_vals_seen )
440  {
441  // ERROR we should not have a value after none are seen
442  std::cerr << " ERROR - new value seen after end of values in "
443  "<ssnorm_selection>"
444  << std::endl;
445  return MB_FAILURE;
446  }
447  // Otherwise get the value string from e and convert it to an int
448  int* valp = new int;
449  *valp = atoi( e );
450  ssTagValues.push_back( (const char*)valp );
451  }
452  else
453  {
454  // Otherwise there is no '=' so push a null on the list
455  end_vals_seen = true;
456  ssTagValues.push_back( (const char*)0 );
457  }
458  }
459  }
460  else
461  {
462  std::cerr << " ERROR - missing <ssnorm_selection>" << std::endl;
463  return MB_FAILURE;
464  }
465  }
466  else if( argv[npos] == std::string( "-ropts" ) )
467  {
468  // Parse out the mesh file read options
469  npos++;
470  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
471  readOpts = argv[npos++];
472  else
473  {
474  std::cerr << " ERROR - missing <roptions>" << std::endl;
475  return MB_FAILURE;
476  }
477  }
478  else if( argv[npos] == std::string( "-outfile" ) )
479  {
480  // Parse out the output file name
481  npos++;
482  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
483  outFile = argv[npos++];
484  else
485  {
486  std::cerr << " ERROR - missing <out_file>" << std::endl;
487  return MB_FAILURE;
488  }
489  }
490  else if( argv[npos] == std::string( "-wopts" ) )
491  {
492  // Parse out the output file write options
493  npos++;
494  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
495  writeOpts = argv[npos++];
496  else
497  {
498  std::cerr << " ERROR - missing <woptions>" << std::endl;
499  return MB_FAILURE;
500  }
501  }
502  else if( argv[npos] == std::string( "-dbgout" ) )
503  {
504  // Parse out the debug output file name.
505  // If no name then use the default.
506  npos++;
507  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
508  dbgFile = argv[npos++];
509  else
510  dbgFile = defaultDbgFile;
511  }
512  else
513  {
514  // Unrecognized parameter. Skip it and move along.
515  std::cerr << " ERROR - Unrecognized parameter:" << argv[npos] << std::endl;
516  std::cerr << " Skipping..." << std::endl;
517  npos++;
518  }
519  }
520 
521  if( !haveMeshes )
522  {
523  meshFiles.resize( 2 );
524  meshFiles[0] = std::string( TestDir + "unittest/64bricks_1khex.h5m" );
525  meshFiles[1] = std::string( TestDir + "unittest/64bricks_12ktet.h5m" );
526  if( 0 == rank )
527  std::cout << "Mesh files not entered; using default files " << meshFiles[0] << " and " << meshFiles[1]
528  << std::endl;
529  }
530 
531  if( !haveInterpTag )
532  {
533  interpTag = "vertex_field";
534  std::cout << "Interpolation field name not given, using default of " << interpTag << std::endl;
535  }
536 
537 #ifdef MOAB_HAVE_HDF5
538  if( 1 == argc )
539  {
540  std::stringstream dfname;
541  dfname << "dum" << nprocs << ".h5m";
542  outFile = dfname.str();
543  if( 0 == rank ) std::cout << "No arguments given; using output file " << outFile << std::endl;
544  }
545 #endif
546 
547  return MB_SUCCESS;
548 }
549 
550 // End get_file_options()
551 
552 ErrorCode test_interpolation( Interface* mbImpl,
553  Coupler::Method method,
554  std::string& interpTag,
555  std::string& gNormTag,
556  std::string& ssNormTag,
557  std::vector< const char* >& ssTagNames,
558  std::vector< const char* >& ssTagValues,
559  EntityHandle* roots,
560  std::vector< ParallelComm* >& pcs,
561  double& instant_time,
562  double& pointloc_time,
563  double& interp_time,
564  double& gnorm_time,
565  double& ssnorm_time,
566  double& toler )
567 {
568  assert( method >= Coupler::CONSTANT && method <= Coupler::SPHERICAL );
569 
570  // Source is 1st mesh, target is 2nd
571  Range src_elems, targ_elems, targ_verts;
572  ErrorCode result = pcs[0]->get_part_entities( src_elems, 3 );MB_CHK_ERR( result );
573 
574  double start_time = MPI_Wtime();
575 
576  // Instantiate a coupler, which does not initialize the tree yet
577  Coupler mbc( mbImpl, pcs[0], src_elems, 0, false ); // do not initialize tree yet
578  if( Coupler::SPHERICAL == method ) mbc.set_spherical();
579  mbc.initialize_tree(); // it is expensive, but do something different for spherical
580 
581  // Initialize spectral elements, if they exist
582  bool specSou = false, specTar = false;
583  result = mbc.initialize_spectral_elements( roots[0], roots[1], specSou, specTar );
584 
585  instant_time = MPI_Wtime();
586 
587  // Get points from the target mesh to interpolate
588  // We have to treat differently the case when the target is a spectral mesh
589  // In that case, the points of interest are the GL points, not the vertex nodes
590  std::vector< double > vpos; // This will have the positions we are interested in
591  int numPointsOfInterest = 0;
592  if( !specTar )
593  { // Usual case
594  Range tmp_verts;
595 
596  // First get all vertices adj to partition entities in target mesh
597  result = pcs[1]->get_part_entities( targ_elems, 3 );MB_CHK_ERR( result );
598  if( Coupler::SPHERICAL == method ) result = pcs[1]->get_part_entities( targ_elems, 2 );MB_CHK_ERR( result ); // get the polygons/quads on a sphere.
599  if( Coupler::CONSTANT == method )
600  targ_verts = targ_elems;
601  else
602  result = mbImpl->get_adjacencies( targ_elems, 0, false, targ_verts, Interface::UNION );MB_CHK_ERR( result );
603 
604  // Then get non-owned verts and subtract
605  result = pcs[1]->get_pstatus_entities( 0, PSTATUS_NOT_OWNED, tmp_verts );MB_CHK_ERR( result );
606  targ_verts = subtract( targ_verts, tmp_verts );
607  // get position of these entities; these are the target points
608  numPointsOfInterest = (int)targ_verts.size();
609  vpos.resize( 3 * targ_verts.size() );
610  result = mbImpl->get_coords( targ_verts, &vpos[0] );MB_CHK_ERR( result );
611  // Locate those points in the source mesh
612  std::cout << "rank " << pcs[0]->proc_config().proc_rank() << " points of interest: " << numPointsOfInterest
613  << "\n";
614  result = mbc.locate_points( &vpos[0], numPointsOfInterest, 0, toler );MB_CHK_ERR( result );
615  }
616  else
617  {
618  // In this case, the target mesh is spectral, we want values
619  // interpolated on the GL positions; for each element, get the GL points, and construct
620  // CartVect!!!
621  result = pcs[1]->get_part_entities( targ_elems, 3 );MB_CHK_ERR( result );
622  result = mbc.get_gl_points_on_elements( targ_elems, vpos, numPointsOfInterest );MB_CHK_ERR( result );
623  std::cout << "rank " << pcs[0]->proc_config().proc_rank() << " points of interest: " << numPointsOfInterest
624  << "\n";
625  }
626 
627  pointloc_time = MPI_Wtime();
628 
629  // Now interpolate tag onto target points
630  std::vector< double > field( numPointsOfInterest );
631 
632  result = mbc.interpolate( method, interpTag, &field[0] );MB_CHK_ERR( result );
633 
634  interp_time = MPI_Wtime();
635 
636  // Do global normalization if specified
637  if( !gNormTag.empty() )
638  {
639  // Normalize the source mesh
640  result = mbc.normalize_mesh( roots[0], gNormTag.c_str(), Coupler::VOLUME, 4 );MB_CHK_ERR( result );
641 
642  // Normalize the target mesh
643  result = mbc.normalize_mesh( roots[1], gNormTag.c_str(), Coupler::VOLUME, 4 );MB_CHK_ERR( result );
644  }
645 
646  gnorm_time = MPI_Wtime();
647 
648  // Do subset normalization if specified
649 
650  if( !ssNormTag.empty() )
651  {
652 
653  result = mbc.normalize_subset( roots[0], ssNormTag.c_str(), &ssTagNames[0], ssTagNames.size(), &ssTagValues[0],
654  Coupler::VOLUME, 4 );MB_CHK_ERR( result );
655 
656  result = mbc.normalize_subset( roots[1], ssNormTag.c_str(), &ssTagNames[0], ssTagNames.size(), &ssTagValues[0],
657  Coupler::VOLUME, 4 );MB_CHK_ERR( result );
658  }
659 
660  ssnorm_time = MPI_Wtime();
661 
662  ssnorm_time -= gnorm_time;
663  gnorm_time -= interp_time;
664  interp_time -= pointloc_time;
665  pointloc_time -= instant_time;
666  instant_time -= start_time;
667 
668  // Set field values as tag on target vertices
669  if( specSou )
670  {
671  // Create a new tag for the values on the target
672  Tag tag;
673  std::string newtag = interpTag + "_TAR";
674  result = mbImpl->tag_get_handle( newtag.c_str(), 1, MB_TYPE_DOUBLE, tag, MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_ERR( result );
675  result = mbImpl->tag_set_data( tag, targ_verts, &field[0] );MB_CHK_ERR( result );
676  }
677  else
678  {
679  if( !specTar )
680  {
681  // Use original tag
682  Tag tag;
683  result = mbImpl->tag_get_handle( interpTag.c_str(), 1, MB_TYPE_DOUBLE, tag );MB_CHK_ERR( result );
684  result = mbImpl->tag_set_data( tag, targ_verts, &field[0] );MB_CHK_ERR( result );
685  }
686  else
687  {
688  // We have the field values computed at each GL points, on each element
689  // in the target mesh
690  // We need to create a new tag, on elements, of size _ntot, to hold
691  // all those values.
692  // So it turns out we need ntot. Maybe we can compute it from the
693  // number of values computed, divided by number of elements
694  int ntot = numPointsOfInterest / targ_elems.size();
695  Tag tag;
696  std::string newtag = interpTag + "_TAR";
697  result = mbImpl->tag_get_handle( newtag.c_str(), ntot, MB_TYPE_DOUBLE, tag, MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_ERR( result );
698  result = mbImpl->tag_set_data( tag, targ_elems, &field[0] );MB_CHK_ERR( result );
699  }
700  }
701 
702  // Done
703  return MB_SUCCESS;
704 }
705 
706 #else
707 
708 int main( int /*argc*/, char** /*argv*/ )
709 {
710  print_usage();
711  return 0;
712 }
713 
714 #endif