6 #if defined( __MINGW32__ )
54 MPI_Comm_size( comm, &nprocs );
55 MPI_Comm_rank( comm, &
rank );
67 assert( nghlayers < 10 );
74 read_options =
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
75 "SHARED_ENTS;PARALLEL_GHOSTS=2.0.";
79 read_options =
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
80 "SHARED_ENTS;PARALLEL_GHOSTS=1.0.";
84 MB_SET_ERR( MB_FAILURE,
"unsupported dimension" );
93 read_options =
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
120 assert( !pc && degree &&
dim );
127 std::vector< int >& degs2fit,
132 std::vector< double >& geoml1errs,
133 std::vector< double >& geoml2errs,
134 std::vector< double >& geomlinferrs )
144 MPI_Comm_size( comm, &nprocs );
145 MPI_Comm_rank( comm, &
rank );
151 for(
size_t i = 0; i < degs2fit.size(); ++i )
153 mxdeg = std::max( degs2fit[i], mxdeg );
157 Range elems, elems_owned;
190 double currcoords[3], exactcoords[3];
279 std::vector< double > testpnts, testnaturalcoords;
282 testnaturalcoords.reserve( nvpe * elems_owned.
size() *
nsamples );
287 std::vector< EntityHandle > conn;
289 std::vector< double > elemcoords( 3 * conn.size() );
297 double a = (double)rand() / RAND_MAX, b = (double)rand() / RAND_MAX, c = (double)rand() / RAND_MAX,
sum;
310 assert( a > 0 && b > 0 && c > 0 && fabs( a + b + c - 1 ) < 1e-12 );
311 testpnts.push_back( a * elemcoords[0] + b * elemcoords[3] + c * elemcoords[6] );
312 testpnts.push_back( a * elemcoords[1] + b * elemcoords[4] + c * elemcoords[7] );
313 testpnts.push_back( a * elemcoords[2] + b * elemcoords[5] + c * elemcoords[8] );
314 testnaturalcoords.push_back( a ), testnaturalcoords.push_back( b ), testnaturalcoords.push_back( c );
318 double xi = (double)rand() / RAND_MAX, eta = (double)rand() / RAND_MAX, N[4];
321 N[0] = ( 1 - xi ) * ( 1 - eta ) / 4, N[1] = ( 1 + xi ) * ( 1 - eta ) / 4,
322 N[2] = ( 1 + xi ) * ( 1 + eta ) / 4, N[3] = ( 1 - xi ) * ( 1 + eta ) / 4;
323 testpnts.push_back( N[0] * elemcoords[0] + N[1] * elemcoords[3] + N[2] * elemcoords[6] +
324 N[3] * elemcoords[9] );
325 testpnts.push_back( N[0] * elemcoords[1] + N[1] * elemcoords[4] + N[2] * elemcoords[7] +
326 N[3] * elemcoords[10] );
327 testpnts.push_back( N[0] * elemcoords[2] + N[1] * elemcoords[5] + N[2] * elemcoords[8] +
328 N[3] * elemcoords[11] );
329 testnaturalcoords.push_back( N[0] ), testnaturalcoords.push_back( N[1] ),
330 testnaturalcoords.push_back( N[2] ), testnaturalcoords.push_back( N[3] );
334 throw std::invalid_argument(
"NOT SUPPORTED ELEMENT TYPE" );
340 double l1err, l2err, linferr;
343 geomlinferrs.clear();
345 geoml1errs.push_back( l1err );
346 geoml2errs.push_back( l2err );
347 geomlinferrs.push_back( linferr );
352 for(
size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
375 assert( (
size_t)index == elems_owned.
size() );
378 geoml1errs.push_back( l1err );
379 geoml2errs.push_back( l2err );
380 geomlinferrs.push_back( linferr );
414 std::cout <<
"usage: mpirun -np <number of processors> ./uref_hirec_test_parallel <prefix of "
415 "files> -str <first number> -end <lastnumber> -suffix <suffix> -interp <0=least "
416 "square, 1=interpolation> -dim <mesh dimension> -geom <s=sphere,t=torus>"
420 int main(
int argc,
char* argv[] )
423 MPI_Init( &argc, &argv );
428 std::string prefix, suffix, istr, iend;
436 prefix = TestDir +
"unittest/sphere_quads_5_ML_";
440 std::cout <<
"Using defaults: MeshFiles/unittest/sphere_quads_5_ML_*.h5m -str 0 -end 1 "
441 "-suffix \".h5m\" -interp 0 -dim 2 -geom s"
446 prefix = std::string( argv[1] );
449 for(
int i = 2; i < argc; ++i )
453 if( std::string( argv[i] ) ==
"-suffix" )
455 suffix = std::string( argv[++i] );
457 else if( std::string( argv[i] ) ==
"-str" )
459 istr = std::string( argv[++i] );
461 else if( std::string( argv[i] ) ==
"-end" )
463 iend = std::string( argv[++i] );
465 else if( std::string( argv[i] ) ==
"-interp" )
467 interp = atoi( argv[++i] );
469 else if( std::string( argv[i] ) ==
"-dim" )
471 dim = atoi( argv[++i] );
474 else if( std::string( argv[i] ) ==
"-geom" )
476 geom = std::string( argv[++i] ) ==
"s" ? 0 : 1;
502 std::cout <<
"Dimension of input mesh should be provided, positive and less than 3" << std::endl;
506 std::cout <<
"Dimension of input mesh should be provided, positive and less than 3" << std::endl;
517 std::cout <<
"Input dimesion should be positive and less than 3;" << std::endl;
521 std::cout <<
"Input dimesion should be positive and less than 3;" << std::endl;
530 std::cout <<
"Testing on " << prefix + istr +
"-" + iend + suffix <<
" with dimension " <<
dim <<
"\n";
531 std::string opts = interp ?
"interpolation" :
"least square fitting";
532 std::cout <<
"High order reconstruction in " << opts << std::endl;
536 std::cout <<
"Testing on " << prefix + istr +
"-" + iend + suffix <<
" with dimension " <<
dim <<
"\n";
537 std::string opts = interp ?
"interpolation" :
"least square fitting";
538 std::cout <<
"High order reconstruction in " << opts << std::endl;
553 std::vector< int > degs2fit;
554 for(
int d = 1; d <= 6; ++d )
556 degs2fit.push_back( d );
559 int begin = atoi( istr.c_str() ), end = atoi( iend.c_str() ) + 1;
561 std::vector< std::vector< double > > geoml1errs_global, geoml2errs_global, geomlinferrs_global;
566 std::vector< std::vector< double > >( 1 + degs2fit.size(), std::vector< double >( end - begin, 0 ) );
568 std::vector< std::vector< double > >( 1 + degs2fit.size(), std::vector< double >( end - begin, 0 ) );
569 geomlinferrs_global =
570 std::vector< std::vector< double > >( 1 + degs2fit.size(), std::vector< double >( end - begin, 0 ) );
574 std::vector< std::vector< double > > geoml1errs_global( 1 + degs2fit.size(),
575 std::vector< double >( end - begin, 0 ) );
576 std::vector< std::vector< double > > geoml2errs_global( 1 + degs2fit.size(),
577 std::vector< double >( end - begin, 0 ) );
578 std::vector< std::vector< double > > geomlinferrs_global( 1 + degs2fit.size(),
579 std::vector< double >( end - begin, 0 ) );
582 for(
int i = begin; i < end; ++i )
584 std::ostringstream convert;
586 std::string infile = prefix + convert.str() + suffix;
588 std::vector< double > geoml1errs, geoml2errs, geomlinferrs;
591 std::cout <<
"Processor " <<
rank <<
" is working on file " << infile << std::endl;
595 assert( geoml1errs.size() == 1 + degs2fit.size() && geoml2errs.size() == 1 + degs2fit.size() &&
596 geomlinferrs.size() == 1 + degs2fit.size() );
601 int ntestverts_global = 0;
602 MPI_Reduce( &ntestverts, &ntestverts_global, 1, MPI_INT, MPI_SUM, 0,
MPI_COMM_WORLD );
604 for(
size_t d = 0; d < degs2fit.size() + 1; ++d )
606 double local_l1err = ntestverts * geoml1errs[d],
607 local_l2err = geoml2errs[d] * ( geoml2errs[d] * ntestverts ), local_linferr = geomlinferrs[d];
608 std::cout <<
"On Processor " <<
rank <<
" with mesh " << i
609 <<
" Degree = " << ( d == 0 ? 0 : degs2fit[d - 1] ) <<
" L1:" << geoml1errs[d]
610 <<
" L2:" << geoml2errs[d] <<
" Li:" << geomlinferrs[d] << std::endl;
611 double global_l1err = 0, global_l2err = 0, global_linferr = 0;
612 MPI_Reduce( &local_l1err, &global_l1err, 1, MPI_DOUBLE, MPI_SUM, 0,
MPI_COMM_WORLD );
613 MPI_Reduce( &local_l2err, &global_l2err, 1, MPI_DOUBLE, MPI_SUM, 0,
MPI_COMM_WORLD );
614 MPI_Reduce( &local_linferr, &global_linferr, 1, MPI_DOUBLE, MPI_MAX, 0,
MPI_COMM_WORLD );
618 geoml1errs_global[d][i - begin] = global_l1err / ntestverts_global;
619 geoml2errs_global[d][i - begin] = sqrt( global_l2err / ntestverts_global );
620 geomlinferrs_global[d][i - begin] = global_linferr;
626 for(
size_t d = 0; d < degs2fit.size() + 1; ++d )
628 geoml1errs_global[d][i - begin] = geoml1errs[d];
629 geoml2errs_global[d][i - begin] = geoml2errs[d];
630 geomlinferrs_global[d][i - begin] = geomlinferrs[d];
636 for(
size_t d = 0; d < degs2fit.size() + 1; ++d )
638 geoml1errs_global[d][i - begin] = geoml1errs[d];
639 geoml2errs_global[d][i - begin] = geoml2errs[d];
640 geomlinferrs_global[d][i - begin] = geomlinferrs[d];
650 std::cout <<
"Degrees: 0 ";
652 for(
size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
654 std::cout << degs2fit[ideg] <<
" ";
657 std::cout << std::endl;
658 std::cout <<
"L1-norm error: \n";
660 for(
size_t i = 0; i < geoml1errs_global.size(); ++i )
662 for(
size_t j = 0; j < geoml1errs_global[i].size(); ++j )
664 std::cout << geoml1errs_global[i][j] <<
" ";
667 std::cout << std::endl;
670 std::cout <<
"L2-norm error: \n";
672 for(
size_t i = 0; i < geoml2errs_global.size(); ++i )
674 for(
size_t j = 0; j < geoml2errs_global[i].size(); ++j )
676 std::cout << geoml2errs_global[i][j] <<
" ";
679 std::cout << std::endl;
682 std::cout <<
"Linf-norm error: \n";
684 for(
size_t i = 0; i < geomlinferrs_global.size(); ++i )
686 for(
size_t j = 0; j < geomlinferrs_global[i].size(); ++j )
688 std::cout << geomlinferrs_global[i][j] <<
" ";
691 std::cout << std::endl;
696 std::cout <<
"Degrees: 0 ";
698 for(
size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
700 std::cout << degs2fit[ideg] <<
" ";
703 std::cout << std::endl;
704 std::cout <<
"L1-norm error: \n";
706 for(
size_t i = 0; i < geoml1errs_global.size(); ++i )
708 for(
size_t j = 0; j < geoml1errs_global[i].size(); ++j )
710 std::cout << geoml1errs_global[i][j] <<
" ";
713 std::cout << std::endl;
716 std::cout <<
"L2-norm error: \n";
718 for(
size_t i = 0; i < geoml2errs_global.size(); ++i )
720 for(
size_t j = 0; j < geoml2errs_global[i].size(); ++j )
722 std::cout << geoml2errs_global[i][j] <<
" ";
725 std::cout << std::endl;
728 std::cout <<
"Linf-norm error: \n";
730 for(
size_t i = 0; i < geomlinferrs_global.size(); ++i )
732 for(
size_t j = 0; j < geomlinferrs_global[i].size(); ++j )
734 std::cout << geomlinferrs_global[i][j] <<
" ";
737 std::cout << std::endl;