 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 // MOAB includes
2 #include "moab/ParallelComm.hpp"
3 #include "MBParallelConventions.h"
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
33 void 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;
123 Coupler::Method method = Coupler::CONSTANT;
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
1.9.1.