1 #include "moab/Core.hpp"
2 #include "moab/CpuTimer.hpp"
3 #include "DataCoupler.hpp"
4 #include "ElemUtil.hpp"
5 #include "TestUtil.hpp"
6
7 #include <iostream>
8 #include <iomanip>
9 #include <sstream>
10 #include <cassert>
11
12 #ifdef MOAB_HAVE_MPI
13 #include "moab/ParallelComm.hpp"
14 #include "MBParallelConventions.h"
15 #endif
16
17 using namespace moab;
18
19 #define PRINT_LAST_ERROR \
20 if( MB_SUCCESS != result ) \
21 { \
22 std::string tmp_str; \
23 std::cout << "Failure; message:" << std::endl; \
24 mbImpl->get_last_error( tmp_str ); \
25 std::cout << tmp_str << std::endl; \
26 MPI_Abort( MPI_COMM_WORLD, result ); \
27 return result; \
28 }
29
30
31 void print_usage( char** argv )
32 {
33 std::cerr << "Usage: ";
34 std::cerr << argv[0]
35 << " -meshes <source_mesh> <target_mesh> -itag <interp_tag> [-gnorm <gnorm_tag>] "
36 "[-ssnorm <ssnorm_tag> <ssnorm_selection>] [-ropts <roptions>] [-outfile "
37 "<out_file> [-wopts <woptions>]] [-dbgout [<dbg_file>]]"
38 << std::endl;
39 std::cerr << " -meshes" << std::endl;
40 std::cerr << " Read in mesh files <source_mesh> and <target_mesh>." << std::endl;
41 std::cerr << " -itag" << std::endl;
42 std::cerr << " Interpolate tag <interp_tag> from source mesh to target mesh." << std::endl;
43 std::cerr << " -gnorm" << std::endl;
44 std::cerr << " Normalize the value of tag <gnorm_tag> over then entire mesh and save to" << std::endl;
45 std::cerr << " tag \"<gnorm_tag>_normf\" on the mesh set. Do this for all meshes." << std::endl;
46 std::cerr << " -ssnorm" << std::endl;
47 std::cerr << " Normalize the value of tag <ssnorm_tag> over subsets of a mesh and save to" << std::endl;
48 std::cerr << " tag \"<ssnorm_tag>_normf\" on the Entity Set for each subset. Subsets "
49 "are selected"
50 << std::endl;
51 std::cerr << " using criteria in <ssnorm_selection>. Do this for all meshes." << std::endl;
52 std::cerr << " -ropts" << std::endl;
53 std::cerr << " Read in the mesh files using options in <roptions>." << std::endl;
54 std::cerr << " -outfile" << std::endl;
55 std::cerr << " Write out target mesh to <out_file>." << std::endl;
56 std::cerr << " -wopts" << std::endl;
57 std::cerr << " Write out mesh files using options in <woptions>." << std::endl;
58 std::cerr << " -dbgout" << std::endl;
59 std::cerr << " Write stdout and stderr streams to the file \'<dbg_file>.txt\'." << std::endl;
60 std::cerr << " -eps" << std::endl;
61 std::cerr << " epsilon" << std::endl;
62 std::cerr << " -meth <method> (0=CONSTANT, 1=LINEAR_FE, 2=QUADRATIC_FE, 3=SPECTRAL)" << std::endl;
63 }
64
65 #ifdef MOAB_HAVE_HDF5
66
67 ErrorCode get_file_options( int argc,
68 char** argv,
69 std::vector< std::string >& meshFiles,
70 DataCoupler::Method& method,
71 std::string& interpTag,
72 std::string& gNormTag,
73 std::string& ssNormTag,
74 std::vector< const char* >& ssTagNames,
75 std::vector< const char* >& ssTagValues,
76 std::string& readOpts,
77 std::string& outFile,
78 std::string& writeOpts,
79 std::string& dbgFile,
80 bool& help,
81 double& epsilon );
82
83
84
85
86
87
88
89 #ifdef MOAB_HAVE_MPI
90 ErrorCode report_iface_ents( Interface* mbImpl, std::vector< ParallelComm* >& pcs, bool print_results );
91 #endif
92
93 ErrorCode test_interpolation( Interface* mbImpl,
94 DataCoupler::Method method,
95 std::string& interpTag,
96 std::string& gNormTag,
97 std::string& ssNormTag,
98 std::vector< const char* >& ssTagNames,
99 std::vector< const char* >& ssTagValues,
100 std::vector< ParallelComm* >& pcs,
101 double& instant_time,
102 double& pointloc_time,
103 double& interp_time,
104 double& gnorm_time,
105 double& ssnorm_time,
106 double& toler );
107
108 void reduceMax( double& v )
109 {
110 double buf;
111
112 MPI_Barrier( MPI_COMM_WORLD );
113 MPI_Allreduce( &v, &buf, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
114
115 v = buf;
116 }
117
118 int main( int argc, char** argv )
119 {
120 #ifdef MOAB_HAVE_MPI
121
122 int err = MPI_Init( &argc, &argv );
123 if( err != 0 )
124 {
125 std::cout << "MPI Initialization did not succeed.\n";
126 exit( 1 );
127 }
128 #endif
129
130 std::vector< const char* > ssTagNames, ssTagValues;
131 std::vector< std::string > meshFiles;
132 std::string interpTag, gNormTag, ssNormTag, readOpts, outFile, writeOpts, dbgFile;
133 DataCoupler::Method method = DataCoupler::CONSTANT;
134
135 ErrorCode result = MB_SUCCESS;
136 bool help = false;
137 double toler = 5.e-10;
138 result = get_file_options( argc, argv, meshFiles, method, interpTag, gNormTag, ssNormTag, ssTagNames, ssTagValues,
139 readOpts, outFile, writeOpts, dbgFile, help, toler );
140
141 if( result != MB_SUCCESS || help )
142 {
143 print_usage( argv );
144 #ifdef MOAB_HAVE_MPI
145 err = MPI_Finalize();
146 if( err != 0 )
147 {
148 std::cout << "MPI Initialization did not succeed.\n";
149 exit( 1 );
150 }
151 #endif
152 return 1;
153 }
154
155 int nprocs = 1, rank = 0;
156
157 #ifdef MOAB_HAVE_MPI
158 err = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
159 err = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
160 std::vector< ParallelComm* > pcs( meshFiles.size() );
161 #endif
162
163
164 if( !dbgFile.empty() )
165 {
166 std::stringstream dfname;
167 dfname << dbgFile << rank << ".txt";
168 if( !std::freopen( dfname.str().c_str(), "a", stdout ) ) return 2;
169 if( !std::freopen( dfname.str().c_str(), "a", stderr ) ) return 2;
170 }
171
172
173 Interface* mbImpl = new( std::nothrow ) Core();
174 if( NULL == mbImpl ) return 1;
175
176
177
178
179 std::vector< EntityHandle > roots( meshFiles.size() );
180
181 for( unsigned int i = 0; i < meshFiles.size(); i++ )
182 {
183 std::string newReadopts;
184 std::ostringstream extraOpt;
185 #ifdef MOAB_HAVE_MPI
186 pcs[i] = new ParallelComm( mbImpl, MPI_COMM_WORLD );
187 int index = pcs[i]->get_id();
188 extraOpt << ";PARALLEL_COMM=" << index;
189 newReadopts = readOpts + extraOpt.str();
190 #endif
191
192 result = mbImpl->create_meshset( MESHSET_SET, roots[i] );
193 PRINT_LAST_ERROR;
194 result = mbImpl->load_file( meshFiles[i].c_str(), &roots[i], newReadopts.c_str() );
195 PRINT_LAST_ERROR;
196 }
197
198 #ifdef MOAB_HAVE_MPI
199 result = report_iface_ents( mbImpl, pcs, true );
200 PRINT_LAST_ERROR;
201 #endif
202
203 double instant_time = 0.0, pointloc_time = 0.0, interp_time = 0.0, gnorm_time = 0.0, ssnorm_time = 0.0;
204
205
206 result = test_interpolation( mbImpl, method, interpTag, gNormTag, ssNormTag, ssTagNames, ssTagValues, pcs,
207 instant_time, pointloc_time, interp_time, gnorm_time, ssnorm_time, toler );
208 PRINT_LAST_ERROR;
209
210 reduceMax( instant_time );
211 reduceMax( pointloc_time );
212 reduceMax( interp_time );
213
214 if( 0 == rank )
215 printf( "\nMax time : %g %g %g (inst loc interp -- %d procs)\n", instant_time, pointloc_time, interp_time,
216 nprocs );
217
218
219 if( !outFile.empty() )
220 {
221 Range partSets;
222
223 partSets.insert( (EntityHandle)roots[1] );
224 std::string newwriteOpts = writeOpts;
225 std::ostringstream extraOpt;
226 #ifdef MOAB_HAVE_MPI
227 extraOpt << ";PARALLEL_COMM=" << 1;
228 newwriteOpts += extraOpt.str();
229 #endif
230 result = mbImpl->write_file( outFile.c_str(), NULL, newwriteOpts.c_str(), partSets );
231 PRINT_LAST_ERROR;
232 std::cout << "Wrote " << outFile << std::endl;
233 std::cout << "mbcoupler_test complete." << std::endl;
234 }
235
236 #ifdef MOAB_HAVE_MPI
237 for( unsigned int i = 0; i < meshFiles.size(); i++ )
238 delete pcs[i];
239 #endif
240
241 delete mbImpl;
242
243
244 #ifdef MOAB_HAVE_MPI
245 err = MPI_Finalize();
246 #endif
247
248 return 0;
249 }
250
251 #ifdef MOAB_HAVE_MPI
252 ErrorCode report_iface_ents( Interface* mbImpl, std::vector< ParallelComm* >& pcs, const bool print_results )
253 {
254 Range iface_ents[6];
255 ErrorCode result = MB_SUCCESS, tmp_result;
256
257
258 for( unsigned int p = 0; p < pcs.size(); p++ )
259 {
260 for( int i = 0; i < 4; i++ )
261 {
262 tmp_result = pcs[p]->get_iface_entities( -1, i, iface_ents[i] );
263
264 if( MB_SUCCESS != tmp_result )
265 {
266 std::cerr << "get_iface_entities returned error on proc " << pcs[p]->proc_config().proc_rank()
267 << "; message: " << std::endl;
268 std::string last_error;
269 result = mbImpl->get_last_error( last_error );
270 if( last_error.empty() )
271 std::cerr << "(none)" << std::endl;
272 else
273 std::cerr << last_error << std::endl;
274 result = tmp_result;
275 }
276 if( 0 != i ) iface_ents[4].merge( iface_ents[i] );
277 }
278 }
279
280
281 result = mbImpl->get_adjacencies( iface_ents[4], 0, false, iface_ents[5], Interface::UNION );
282
283 int rank;
284 MPI_Comm_rank( MPI_COMM_WORLD, &rank );
285
286 if( print_results || iface_ents[0].size() != iface_ents[5].size() )
287 {
288 std::cerr << "Proc " << rank << " iface entities: " << std::endl;
289 for( int i = 0; i < 4; i++ )
290 std::cerr << " " << iface_ents[i].size() << " " << i << "d iface entities." << std::endl;
291 std::cerr << " (" << iface_ents[5].size() << " verts adj to other iface ents)" << std::endl;
292 }
293
294 return result;
295 }
296 #endif
297
298
299
300 bool check_for_flag( const char* str )
301 {
302 if( '-' == str[0] )
303 return true;
304 else
305 return false;
306 }
307
308
309 ErrorCode get_file_options( int argc,
310 char** argv,
311 std::vector< std::string >& meshFiles,
312 DataCoupler::Method& method,
313 std::string& interpTag,
314 std::string& gNormTag,
315 std::string& ssNormTag,
316 std::vector< const char* >& ssTagNames,
317 std::vector< const char* >& ssTagValues,
318 std::string& readOpts,
319 std::string& outFile,
320 std::string& writeOpts,
321 std::string& dbgFile,
322 bool& help,
323 double& epsilon )
324 {
325
326
327 gNormTag = "";
328 ssNormTag = "";
329 readOpts = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;PARALLEL_"
330 "RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=3.0.1;CPUTIME";
331 outFile = "";
332 writeOpts = "PARALLEL=WRITE_PART;CPUTIME";
333 dbgFile = "";
334 std::string defaultDbgFile = argv[0];
335
336
337 bool haveMeshes = false;
338 bool haveInterpTag = false;
339
340
341 int npos = 1;
342
343 if( argc > 1 && argv[1] == std::string( "-h" ) )
344 {
345 help = true;
346 return MB_SUCCESS;
347 }
348
349 while( npos < argc )
350 {
351 if( argv[npos] == std::string( "-meshes" ) )
352 {
353
354 npos++;
355 int numFiles = 2;
356 meshFiles.resize( numFiles );
357 for( int i = 0; i < numFiles; i++ )
358 {
359 if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
360 meshFiles[i] = argv[npos++];
361 else
362 {
363 std::cerr << " ERROR - missing correct number of mesh filenames" << std::endl;
364 return MB_FAILURE;
365 }
366 }
367
368 haveMeshes = true;
369 }
370 else if( argv[npos] == std::string( "-itag" ) )
371 {
372
373 npos++;
374 if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
375 interpTag = argv[npos++];
376 else
377 {
378 std::cerr << " ERROR - missing <interp_tag>" << std::endl;
379 return MB_FAILURE;
380 }
381
382 haveInterpTag = true;
383 }
384 else if( argv[npos] == std::string( "-meth" ) )
385 {
386
387 npos++;
388 if( argv[npos][0] == '0' )
389 method = DataCoupler::CONSTANT;
390 else if( argv[npos][0] == '1' )
391 method = DataCoupler::LINEAR_FE;
392 else if( argv[npos][0] == '2' )
393 method = DataCoupler::QUADRATIC_FE;
394 else if( argv[npos][0] == '3' )
395 method = DataCoupler::SPECTRAL;
396 else
397 {
398 std::cerr << " ERROR - unrecognized method number " << method << std::endl;
399 return MB_FAILURE;
400 }
401 npos++;
402 }
403 else if( argv[npos] == std::string( "-eps" ) )
404 {
405
406 npos++;
407 if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
408 epsilon = atof( argv[npos++] );
409 else
410 {
411 std::cerr << " ERROR - missing <epsilon>" << std::endl;
412 return MB_FAILURE;
413 }
414 }
415 else if( argv[npos] == std::string( "-gnorm" ) )
416 {
417
418 npos++;
419 if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
420 gNormTag = argv[npos++];
421 else
422 {
423 std::cerr << " ERROR - missing <gnorm_tag>" << std::endl;
424 return MB_FAILURE;
425 }
426 }
427 else if( argv[npos] == std::string( "-ssnorm" ) )
428 {
429
430 npos++;
431 if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
432 ssNormTag = argv[npos++];
433 else
434 {
435 std::cerr << " ERROR - missing <ssnorm_tag>" << std::endl;
436 return MB_FAILURE;
437 }
438
439 if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
440 {
441 char* opts = argv[npos++];
442 char sep1[1] = { ';' };
443 char sep2[1] = { '=' };
444 bool end_vals_seen = false;
445 std::vector< char* > tmpTagOpts;
446
447
448 for( char* i = strtok( opts, sep1 ); i; i = strtok( 0, sep1 ) )
449 tmpTagOpts.push_back( i );
450
451
452 for( unsigned int j = 0; j < tmpTagOpts.size(); j++ )
453 {
454 char* e = strtok( tmpTagOpts[j], sep2 );
455 ssTagNames.push_back( e );
456 e = strtok( 0, sep2 );
457 if( e != NULL )
458 {
459
460 if( end_vals_seen )
461 {
462
463 std::cerr << " ERROR - new value seen after end of values in "
464 "<ssnorm_selection>"
465 << std::endl;
466 return MB_FAILURE;
467 }
468
469 int* valp = new int;
470 *valp = atoi( e );
471 ssTagValues.push_back( (const char*)valp );
472 }
473 else
474 {
475
476 end_vals_seen = true;
477 ssTagValues.push_back( (const char*)0 );
478 }
479 }
480 }
481 else
482 {
483 std::cerr << " ERROR - missing <ssnorm_selection>" << std::endl;
484 return MB_FAILURE;
485 }
486 }
487 else if( argv[npos] == std::string( "-ropts" ) )
488 {
489
490 npos++;
491 if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
492 readOpts = argv[npos++];
493 else
494 {
495 std::cerr << " ERROR - missing <roptions>" << std::endl;
496 return MB_FAILURE;
497 }
498 }
499 else if( argv[npos] == std::string( "-outfile" ) )
500 {
501
502 npos++;
503 if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
504 outFile = argv[npos++];
505 else
506 {
507 std::cerr << " ERROR - missing <out_file>" << std::endl;
508 return MB_FAILURE;
509 }
510 }
511 else if( argv[npos] == std::string( "-wopts" ) )
512 {
513
514 npos++;
515 if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
516 writeOpts = argv[npos++];
517 else
518 {
519 std::cerr << " ERROR - missing <woptions>" << std::endl;
520 return MB_FAILURE;
521 }
522 }
523 else if( argv[npos] == std::string( "-dbgout" ) )
524 {
525
526
527 npos++;
528 if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
529 dbgFile = argv[npos++];
530 else
531 dbgFile = defaultDbgFile;
532 }
533 else
534 {
535
536 std::cerr << " ERROR - Unrecognized parameter:" << argv[npos] << std::endl;
537 std::cerr << " Skipping..." << std::endl;
538 npos++;
539 }
540 }
541
542 if( !haveMeshes )
543 {
544 meshFiles.resize( 2 );
545 meshFiles[0] = std::string( TestDir + "unittest/64bricks_1khex.h5m" );
546 meshFiles[1] = std::string( TestDir + "unittest/64bricks_12ktet.h5m" );
547 std::cout << "Mesh files not entered; using default files " << meshFiles[0] << " and " << meshFiles[1]
548 << std::endl;
549 }
550
551 if( !haveInterpTag )
552 {
553 interpTag = "vertex_field";
554 std::cout << "Interpolation field name not given, using default of " << interpTag << std::endl;
555 }
556
557 #ifdef MOAB_HAVE_HDF5
558 if( 1 == argc )
559 {
560 std::cout << "No arguments given; using output file dum.h5m." << std::endl;
561 outFile = "dum.h5m";
562 }
563 #endif
564
565 return MB_SUCCESS;
566 }
567
568
569
570 ErrorCode test_interpolation( Interface* mbImpl,
571 DataCoupler::Method method,
572 std::string& interpTag,
573 std::string& ,
574 std::string& ,
575 std::vector< const char* >& ,
576 std::vector< const char* >& ,
577 std::vector< ParallelComm* >& pcs,
578 double& instant_time,
579 double& pointloc_time,
580 double& interp_time,
581 double& ,
582 double& ,
583 double& toler )
584 {
585 assert( method >= DataCoupler::CONSTANT && method <= DataCoupler::SPECTRAL );
586
587
588 Range src_elems, targ_elems, targ_verts;
589 ErrorCode result = pcs[0]->get_part_entities( src_elems, 3 );
590 PRINT_LAST_ERROR;
591
592 CpuTimer timer;
593
594
595 DataCoupler dc( mbImpl, src_elems, 0, pcs[0] );
596
597
598
599
600
601
602 instant_time = timer.time_since_birth();
603
604
605
606
607 std::vector< double > vpos;
608 int numPointsOfInterest = 0;
609 #ifdef MOAB_HAVE_MPI
610 result = pcs[1]->get_part_entities( targ_elems, 3 );
611 #endif
612 PRINT_LAST_ERROR;
613
614
615 if( DataCoupler::CONSTANT == method )
616 targ_verts = targ_elems;
617 else
618 result = mbImpl->get_adjacencies( targ_elems, 0, false, targ_verts, Interface::UNION );
619 PRINT_LAST_ERROR;
620
621 #ifdef MOAB_HAVE_MPI
622
623 Range tmp_verts;
624 result = pcs[1]->get_pstatus_entities( 0, PSTATUS_NOT_OWNED, tmp_verts );
625 PRINT_LAST_ERROR;
626 targ_verts = subtract( targ_verts, tmp_verts );
627 #endif
628
629 numPointsOfInterest = (int)targ_verts.size();
630 vpos.resize( 3 * targ_verts.size() );
631 result = mbImpl->get_coords( targ_verts, &vpos[0] );
632 PRINT_LAST_ERROR;
633
634
635 #ifdef MOAB_HAVE_MPI
636 std::cout << "rank " << pcs[0]->proc_config().proc_rank();
637 #endif
638 std::cout << " points of interest: " << numPointsOfInterest << "\n";
639 result = dc.locate_points( &vpos[0], numPointsOfInterest, toler );
640 PRINT_LAST_ERROR;
641
642 pointloc_time = timer.time_elapsed();
643
644
645 std::vector< double > field( numPointsOfInterest );
646
647 result = dc.interpolate( method, interpTag, &field[0] );
648 PRINT_LAST_ERROR;
649
650 interp_time = timer.time_elapsed();
651
652
653
654 Tag tag;
655 result = mbImpl->tag_get_handle( interpTag.c_str(), 1, MB_TYPE_DOUBLE, tag );
656 PRINT_LAST_ERROR;
657 result = mbImpl->tag_set_data( tag, targ_verts, &field[0] );
658 PRINT_LAST_ERROR;
659
660
661 return MB_SUCCESS;
662 }
663
664 #else
665
666 int main( int , char** argv )
667 {
668 print_usage( argv );
669 return 0;
670 }
671
672 #endif