uref_hirec_test_parallel.cpp File Reference

#include <iostream>
#include <string>
#include <sstream>
#include <ctime>
#include <vector>
#include <algorithm>
#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/MeshTopoUtil.hpp"
#include "moab/NestedRefine.hpp"
#include "moab/DiscreteGeometry/HiReconstruction.hpp"
#include "TestUtil.hpp"
#include "geomObject.cpp"
#include <cmath>
#include <cstdlib>
#include <cassert>

Include dependency graph for uref_hirec_test_parallel.cpp:

Go to the source code of this file.

Macros
#define	nsamples   10

Functions
ErrorCode	load_meshset_hirec (const char *infile, Interface *mbimpl, EntityHandle &meshset, ParallelComm *&pc, const int degree, const int dim)
ErrorCode	closedsurface_uref_hirec_convergence_study (const char *infile, std::vector< int > &degs2fit, bool interp, int dim, geomObject *obj, int &ntestverts, std::vector< double > &geoml1errs, std::vector< double > &geoml2errs, std::vector< double > &geomlinferrs)
void	usage ()
int	main (int argc, char *argv[])

Macro Definition Documentation

nsamples

#define nsamples   10

Definition at line 36 of file uref_hirec_test_parallel.cpp.

Function Documentation

closedsurface_uref_hirec_convergence_study()

ErrorCode closedsurface_uref_hirec_convergence_study	(	const char *	infile,
		std::vector< int > &	degs2fit,
		bool	interp,
		int	dim,
		geomObject *	obj,
		int &	ntestverts,
		std::vector< double > &	geoml1errs,
		std::vector< double > &	geoml2errs,
		std::vector< double > &	geomlinferrs
	)

Definition at line 126 of file uref_hirec_test_parallel.cpp.

{
 Core moab;
 Interface* mbImpl = &moab;
 ParallelComm* pc = NULL;
 EntityHandle meshset;
 
#ifdef MOAB_HAVE_MPI
 int nprocs, rank;
 MPI_Comm comm = MPI_COMM_WORLD;
 MPI_Comm_size( comm, &nprocs );
 MPI_Comm_rank( comm, &rank );
#endif
 
 ErrorCode error;
 // mesh will be loaded and communicator pc will be updated
 int mxdeg = 1;
 for( size_t i = 0; i < degs2fit.size(); ++i )
 {
 mxdeg = std::max( degs2fit[i], mxdeg );
 }
 error = load_meshset_hirec( infile, mbImpl, meshset, pc, mxdeg, dim );MB_CHK_ERR( error );
 
 Range elems, elems_owned;
 error = mbImpl->get_entities_by_dimension( meshset, dim, elems );MB_CHK_ERR( error );
 
#ifdef MOAB_HAVE_MPI
 
 if( pc )
 {
 error = pc->filter_pstatus( elems, PSTATUS_GHOST, PSTATUS_NOT, -1, &elems_owned );MB_CHK_ERR( error );
 }
 else
 {
 elems_owned = elems;
 }
 
#endif
 
#ifdef MOAB_HAVE_MPI
 // std::cout << "Mesh has " << nelems << " elements on Processor " << rank << " in total;";
 // std::cout << elems_owned.size() << " of which are locally owned elements" << std::endl;
#else
 // std::cout << "Mesh has " << nelems << " elements" << std::endl;
#endif
 
 /************************
 * convergence study *
 *************************/
 // project onto exact geometry since each level with uref has only linear coordinates
 Range verts;
 error = mbImpl->get_entities_by_dimension( meshset, 0, verts );MB_CHK_ERR( error );
 
 for( Range::iterator ivert = verts.begin(); ivert != verts.end(); ++ivert )
 {
 EntityHandle currvert = *ivert;
 double currcoords[3], exactcoords[3];
 error = mbImpl->get_coords( &currvert, 1, currcoords );MB_CHK_ERR( error );
 obj->project_points2geom( 3, currcoords, exactcoords, NULL );
 
 error = mbImpl->set_coords( &currvert, 1, exactcoords );MB_CHK_ERR( error );
 // for debug
 /*error = mbImpl->get_coords(&currvert,1,currcoords); MB_CHK_ERR(error);
 assert(currcoords[0]==exactcoords[0]&&currcoords[1]==exactcoords[1]&&currcoords[2]==exactcoords[2]);*/
 }
 
 // test ahf on mesh with ghost layers
 /*Range verts_owned;
 #ifdef MOAB_HAVE_MPI
 if(pc){
 error =
 pc->filter_pstatus(verts,PSTATUS_GHOST,PSTATUS_NOT,-1,&verts_owned);MB_CHK_ERR(error); }else{
 verts_owned = verts;
 }
 
 HalfFacetRep *ahf = new HalfFacetRep(&moab,pc,meshset);
 error = ahf->initialize(); MB_CHK_ERR(error);
 
 for(int i=0;i<nprocs;++i){
 if(rank==i){
 std::cout << "Processor " << rank << " Local elements: \n";
 for(Range::iterator ielem=elems.begin();ielem!=elems.end();++ielem){
 EntityHandle currelem = *ielem;
 std::vector<EntityHandle> conn;
 error = mbImpl->get_connectivity(&currelem,1,conn); MB_CHK_ERR(error);
 std::cout << *ielem << ": ";
 for(size_t k=0;k<conn.size();++k) std::cout << conn[k] << " ";
 std::cout << std::endl;
 }
 
 std::cout << "Processor " << rank << " Local owned elements: \n";
 for(Range::iterator ielem=elems_owned.begin();ielem!=elems_owned.end();++ielem){
 EntityHandle currelem = *ielem;
 std::vector<EntityHandle> conn;
 error = mbImpl->get_connectivity(&currelem,1,conn); MB_CHK_ERR(error);
 std::cout << *ielem << ": ";
 for(size_t k=0;k<conn.size();++k) std::cout << conn[k] << " ";
 std::cout << std::endl;
 }
 
 
 std::cout << "Processor " << rank << " Local vertices: ";
 for(Range::iterator ivert=verts.begin();ivert!=verts.end();++ivert){
 std::cout << *ivert << " ";
 }
 std::cout << std::endl;
 std::cout << "Processor " << rank << " Local owned vertices: ";
 for(Range::iterator ivert=verts_owned.begin();ivert!=verts_owned.end();++ivert){
 std::cout << *ivert << " ";
 }
 std::cout << std::endl;
 
 for(Range::iterator ivert=verts_owned.begin();ivert!=verts_owned.end();++ivert){
 EntityHandle currvid = *ivert;
 std::cout << "Processor " << rank << " local verts: " << *ivert << " has
 adjfaces: "; std::vector<EntityHandle> adjfaces;
 //error = ahf->get_up_adjacencies(currvid,2,adjfaces); MB_CHK_ERR(error);
 error = mbImpl->get_adjacencies(&currvid,1,2,false,adjfaces); MB_CHK_ERR(error);
 for(size_t k=0;k<adjfaces.size();++k){
 std::cout << adjfaces[k] << " ";
 }
 std::cout << std::endl;
 }
 
 for(Range::iterator ivert=verts.begin();ivert!=verts.end();++ivert){
 EntityHandle currvid = *ivert;
 std::cout << "Processor " << rank << " all verts: " << *ivert << " has adjfaces:
 "; std::vector<EntityHandle> adjfaces;
 //error = ahf->get_up_adjacencies(*ivert,2,adjfaces); MB_CHK_ERR(error);
 error = mbImpl->get_adjacencies(&currvid,1,2,false,adjfaces); MB_CHK_ERR(error);
 for(size_t k=0;k<adjfaces.size();++k){
 std::cout << adjfaces[k] << " ";
 }
 std::cout << std::endl;
 }
 
 }else{
 MPI_Barrier(MPI_COMM_WORLD);
 }
 }
 exit(0);
 #endif*/
 
 // generate random points on each elements, assument 3D coordinates
 int nvpe = TYPE_FROM_HANDLE( *elems.begin() ) == MBTRI ? 3 : 4;
 std::vector< double > testpnts, testnaturalcoords;
 ntestverts = elems_owned.size() * nsamples;
 testpnts.reserve( 3 * elems_owned.size() * nsamples );
 testnaturalcoords.reserve( nvpe * elems_owned.size() * nsamples );
 
 for( Range::iterator ielem = elems_owned.begin(); ielem != elems_owned.end(); ++ielem )
 {
 EntityHandle currelem = *ielem;
 std::vector< EntityHandle > conn;
 error = mbImpl->get_connectivity( &currelem, 1, conn );MB_CHK_ERR( error );
 std::vector< double > elemcoords( 3 * conn.size() );
 error = mbImpl->get_coords( &( conn[0] ), conn.size(), &( elemcoords[0] ) );MB_CHK_ERR( error );
 EntityType type = TYPE_FROM_HANDLE( currelem );
 
 for( int s = 0; s < nsamples; ++s )
 {
 if( type == MBTRI )
 {
 double a = (double)rand() / RAND_MAX, b = (double)rand() / RAND_MAX, c = (double)rand() / RAND_MAX, sum;
 sum = a + b + c;
 
 if( sum < 1e-12 )
 {
 --s;
 continue;
 }
 else
 {
 a /= sum, b /= sum, c /= sum;
 }
 
 assert( a > 0 && b > 0 && c > 0 && fabs( a + b + c - 1 ) < 1e-12 );
 testpnts.push_back( a * elemcoords[0] + b * elemcoords[3] + c * elemcoords[6] );
 testpnts.push_back( a * elemcoords[1] + b * elemcoords[4] + c * elemcoords[7] );
 testpnts.push_back( a * elemcoords[2] + b * elemcoords[5] + c * elemcoords[8] );
 testnaturalcoords.push_back( a ), testnaturalcoords.push_back( b ), testnaturalcoords.push_back( c );
 }
 else if( type == MBQUAD )
 {
 double xi = (double)rand() / RAND_MAX, eta = (double)rand() / RAND_MAX, N[4];
 xi = 2 * xi - 1;
 eta = 2 * eta - 1;
 N[0] = ( 1 - xi ) * ( 1 - eta ) / 4, N[1] = ( 1 + xi ) * ( 1 - eta ) / 4,
 N[2] = ( 1 + xi ) * ( 1 + eta ) / 4, N[3] = ( 1 - xi ) * ( 1 + eta ) / 4;
 testpnts.push_back( N[0] * elemcoords[0] + N[1] * elemcoords[3] + N[2] * elemcoords[6] +
 N[3] * elemcoords[9] );
 testpnts.push_back( N[0] * elemcoords[1] + N[1] * elemcoords[4] + N[2] * elemcoords[7] +
 N[3] * elemcoords[10] );
 testpnts.push_back( N[0] * elemcoords[2] + N[1] * elemcoords[5] + N[2] * elemcoords[8] +
 N[3] * elemcoords[11] );
 testnaturalcoords.push_back( N[0] ), testnaturalcoords.push_back( N[1] ),
 testnaturalcoords.push_back( N[2] ), testnaturalcoords.push_back( N[3] );
 }
 else
 {
 throw std::invalid_argument( "NOT SUPPORTED ELEMENT TYPE" );
 }
 }
 }
 
 // Compute error of linear interpolation in PARALLEL
 double l1err, l2err, linferr;
 geoml1errs.clear();
 geoml2errs.clear();
 geomlinferrs.clear();
 obj->compute_projecterror( 3, elems_owned.size() * nsamples, &( testpnts[0] ), l1err, l2err, linferr );
 geoml1errs.push_back( l1err );
 geoml2errs.push_back( l2err );
 geomlinferrs.push_back( linferr );
 /*Perform high order projection and compute error*/
 // initialize
 HiReconstruction hirec( &moab, pc, meshset );
 
 for( size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
 {
 // High order reconstruction
 error = hirec.reconstruct3D_surf_geom( degs2fit[ideg], interp, true, true );MB_CHK_ERR( error );
 
 int index = 0;
 // for debug
 /*double maxlinferr=0,elemlinferr=0,eleml1err=0,eleml2err;
 EntityHandle elem_maxerr;*/
 
 for( Range::iterator ielem = elems_owned.begin(); ielem != elems_owned.end(); ++ielem, ++index )
 {
 // Projection
 error =
 hirec.hiproj_walf_in_element( *ielem, nvpe, nsamples, &( testnaturalcoords[nvpe * nsamples * index] ),
 &( testpnts[3 * nsamples * index] ) );MB_CHK_ERR( error );
 // for debug
 /*obj->compute_projecterror(3,nsamples,&(testpnts[3*nsamples*index]),eleml1err,eleml2err,elemlinferr);
 if(elemlinferr>maxlinferr){
 maxlinferr = elemlinferr; elem_maxerr = *ielem;
 }*/
 }
 
 assert( (size_t)index == elems_owned.size() );
 // Estimate error
 obj->compute_projecterror( 3, elems_owned.size() * nsamples, &( testpnts[0] ), l1err, l2err, linferr );
 geoml1errs.push_back( l1err );
 geoml2errs.push_back( l2err );
 geomlinferrs.push_back( linferr );
 
 // for debug
 /*std::cout << "Degree " << degs2fit[ideg] << " has L-inf error " << maxlinferr << " at
 element " << elems_owned.index(elem_maxerr) << ":"; std::vector<EntityHandle> conn; error =
 mbImpl->get_connectivity(&elem_maxerr,1,conn); MB_CHK_ERR(error); for(size_t
 ii=0;ii<conn.size();++ii) std::cout << verts.index(conn[ii]) << " "; std::cout << std::endl;
 
 for(size_t ii=0;ii<conn.size();++ii){
 //assume triangle
 std::vector<double> vertexcoords(3,0);
 error = mbImpl->get_coords(&(conn[ii]),1,&(vertexcoords[0])); MB_CHK_ERR(error);
 std::cout << verts.index(conn[ii]) << ":" << vertexcoords[0] << " "<< vertexcoords[1] <<
 " "<< vertexcoords[2] << "\n"; GEOMTYPE geomtype; std::vector<double> local_coords_system,
 local_coeffs; int deg_out; bool local_interp; bool hasfit =
 hirec.get_fittings_data(conn[ii],geomtype,local_coords_system,deg_out,local_coeffs,local_interp);
 assert(hasfit); std::cout << "Local fitting result: " << std::endl; std::cout << "Geom Type:
 " << geomtype << std::endl; std::cout << "Local coordinates system: " << std::endl;
 std::cout << local_coords_system[0] << "\t" << local_coords_system[1] << "\t" <<
 local_coords_system[2] << std::endl; std::cout << local_coords_system[3] << "\t" <<
 local_coords_system[4] << "\t" << local_coords_system[5] << std::endl; std::cout <<
 local_coords_system[6] << "\t" << local_coords_system[7] << "\t" << local_coords_system[8]
 << std::endl; std::cout << "Fitting degree is " << deg_out << " interp is " << local_interp
 << std::endl; std::cout << "Fitting coefficients are :" << std::endl; for(size_t
 kk=0;kk<local_coeffs.size();++kk) std::cout << local_coeffs[kk] << "\t"; std::cout <<
 std::endl;
 }*/
 }
 
 return error;
}

References moab::Range::begin(), geomObject::compute_projecterror(), dim, moab::Range::end(), moab::error(), ErrorCode, moab::ParallelComm::filter_pstatus(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::HiReconstruction::hiproj_walf_in_element(), load_meshset_hirec(), MB_CHK_ERR, MBQUAD, MBTRI, MPI_COMM_WORLD, nsamples, geomObject::project_points2geom(), PSTATUS_GHOST, PSTATUS_NOT, rank, moab::HiReconstruction::reconstruct3D_surf_geom(), moab::Interface::set_coords(), moab::Range::size(), moab::sum(), and moab::TYPE_FROM_HANDLE().

Referenced by main().

load_meshset_hirec()

ErrorCode load_meshset_hirec	(	const char *	infile,
		Interface *	mbimpl,
		EntityHandle &	meshset,
		ParallelComm *&	pc,
		const int	degree,
		const int	dim
	)

Definition at line 42 of file uref_hirec_test_parallel.cpp.

{
 ErrorCode error;
 error = mbimpl->create_meshset( moab::MESHSET_SET, meshset );MB_CHK_ERR( error );
#ifdef MOAB_HAVE_MPI
 int nprocs, rank;
 MPI_Comm comm = MPI_COMM_WORLD;
 MPI_Comm_size( comm, &nprocs );
 MPI_Comm_rank( comm, &rank );
 EntityHandle partnset;
 error = mbimpl->create_meshset( moab::MESHSET_SET, partnset );MB_CHK_ERR( error );
 
 if( nprocs > 1 )
 {
 pc = moab::ParallelComm::get_pcomm( mbimpl, partnset, &comm );
 }
 
 if( nprocs > 1 )
 {
 int nghlayers = degree > 0 ? HiReconstruction::estimate_num_ghost_layers( degree, true ) : 0;
 assert( nghlayers < 10 );
 
 if( nghlayers )
 {
 // get ghost layers
 if( dim == 2 )
 {
 read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
 "SHARED_ENTS;PARALLEL_GHOSTS=2.0.";
 }
 else if( dim == 1 )
 {
 read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
 "SHARED_ENTS;PARALLEL_GHOSTS=1.0.";
 }
 else
 {
 MB_SET_ERR( MB_FAILURE, "unsupported dimension" );
 }
 
 read_options += (char)( '0' + nghlayers );
 // std::cout << "On processor " << rank << " Degree=" << degree << " "<< read_options <<
 // std::endl;
 }
 else
 {
 read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
 }
 
 /*for debug*/
 /*std::string outfile = std::string(infile); std::size_t dotpos = outfile.find_last_of(".");
 std::ostringstream convert; convert << rank;
 std::string localfile = outfile.substr(0,dotpos)+convert.str()+".vtk";
 //write local mesh
 EntityHandle local;
 error = mbimpl->create_meshset(moab::MESHSET_SET,local);CHECK_ERR(error);
 std::string local_options =
 "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"; error =
 mbimpl->load_file(infile,&local,local_options.c_str()); MB_CHK_ERR(error); error =
 mbimpl->write_file(localfile.c_str(),0,0,&local,1);*/
 
 error = mbimpl->load_file( infile, &meshset, read_options.c_str() );MB_CHK_ERR( error );
 /*for debug
 //write local mesh with ghost layers
 localfile = outfile.substr(0,dotpos)+convert.str()+"_ghost.vtk";
 error = mbimpl->write_file(localfile.c_str(),0,0,&meshset,1);*/
 }
 else
 {
 error = mbimpl->load_file( infile, &meshset );MB_CHK_ERR( error );
 }
 
#else
 assert( !pc && degree && dim );
 error = mbimpl->load_file( infile, &meshset );MB_CHK_ERR( error );
#endif
 return error;
}

References moab::Interface::create_meshset(), dim, moab::error(), ErrorCode, moab::HiReconstruction::estimate_num_ghost_layers(), moab::ParallelComm::get_pcomm(), moab::Interface::load_file(), MB_CHK_ERR, MB_SET_ERR, MESHSET_SET, MPI_COMM_WORLD, rank, and read_options.

Referenced by closedsurface_uref_hirec_convergence_study().

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 420 of file uref_hirec_test_parallel.cpp.

{
#ifdef MOAB_HAVE_MPI
 MPI_Init( &argc, &argv );
 int nprocs, rank;
 MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
 MPI_Comm_rank( MPI_COMM_WORLD, &rank );
#endif
 std::string prefix, suffix, istr, iend;
 int dim = 2, geom = 0;
 bool interp = false;
 ErrorCode error;
 
 if( argc == 1 )
 {
 usage();
 prefix = TestDir + "unittest/sphere_quads_5_ML_";
 suffix = ".h5m";
 istr = "0";
 iend = "1";
 std::cout << "Using defaults: MeshFiles/unittest/sphere_quads_5_ML_*.h5m -str 0 -end 1 "
 "-suffix \".h5m\" -interp 0 -dim 2 -geom s"
 << std::endl;
 }
 else
 {
 prefix = std::string( argv[1] );
 bool hasdim = false;
 
 for( int i = 2; i < argc; ++i )
 {
 if( i + 1 != argc )
 {
 if( std::string( argv[i] ) == "-suffix" )
 {
 suffix = std::string( argv[++i] );
 }
 else if( std::string( argv[i] ) == "-str" )
 {
 istr = std::string( argv[++i] );
 }
 else if( std::string( argv[i] ) == "-end" )
 {
 iend = std::string( argv[++i] );
 }
 else if( std::string( argv[i] ) == "-interp" )
 {
 interp = atoi( argv[++i] );
 }
 else if( std::string( argv[i] ) == "-dim" )
 {
 dim = atoi( argv[++i] );
 hasdim = true;
 }
 else if( std::string( argv[i] ) == "-geom" )
 {
 geom = std::string( argv[++i] ) == "s" ? 0 : 1;
 }
 else
 {
#ifdef MOAB_HAVE_MPI
 
 if( 0 == rank )
 {
 usage();
 }
 
 MPI_Finalize();
#else
 usage();
#endif
 return 0;
 }
 }
 }
 
 if( !hasdim )
 {
#ifdef MOAB_HAVE_MPI
 
 if( 0 == rank )
 {
 std::cout << "Dimension of input mesh should be provided, positive and less than 3" << std::endl;
 }
 
#else
 std::cout << "Dimension of input mesh should be provided, positive and less than 3" << std::endl;
#endif
 return 0;
 }
 
 if( dim > 2 || dim <= 0 )
 {
#ifdef MOAB_HAVE_MPI
 
 if( 0 == rank )
 {
 std::cout << "Input dimesion should be positive and less than 3;" << std::endl;
 }
 
#else
 std::cout << "Input dimesion should be positive and less than 3;" << std::endl;
#endif
 return 0;
 }
 
#ifdef MOAB_HAVE_MPI
 
 if( 0 == rank )
 {
 std::cout << "Testing on " << prefix + istr + "-" + iend + suffix << " with dimension " << dim << "\n";
 std::string opts = interp ? "interpolation" : "least square fitting";
 std::cout << "High order reconstruction in " << opts << std::endl;
 }
 
#else
 std::cout << "Testing on " << prefix + istr + "-" + iend + suffix << " with dimension " << dim << "\n";
 std::string opts = interp ? "interpolation" : "least square fitting";
 std::cout << "High order reconstruction in " << opts << std::endl;
#endif
 }
 
 geomObject* obj;
 
 if( geom == 0 )
 {
 obj = new sphere();
 }
 else
 {
 obj = new torus();
 }
 
 std::vector< int > degs2fit;
 for( int d = 1; d <= 6; ++d )
 {
 degs2fit.push_back( d );
 }
 
 int begin = atoi( istr.c_str() ), end = atoi( iend.c_str() ) + 1;
#ifdef MOAB_HAVE_MPI
 std::vector< std::vector< double > > geoml1errs_global, geoml2errs_global, geomlinferrs_global;
 
 if( rank == 0 )
 {
 geoml1errs_global =
 std::vector< std::vector< double > >( 1 + degs2fit.size(), std::vector< double >( end - begin, 0 ) );
 geoml2errs_global =
 std::vector< std::vector< double > >( 1 + degs2fit.size(), std::vector< double >( end - begin, 0 ) );
 geomlinferrs_global =
 std::vector< std::vector< double > >( 1 + degs2fit.size(), std::vector< double >( end - begin, 0 ) );
 }
 
#else
 std::vector< std::vector< double > > geoml1errs_global( 1 + degs2fit.size(),
 std::vector< double >( end - begin, 0 ) );
 std::vector< std::vector< double > > geoml2errs_global( 1 + degs2fit.size(),
 std::vector< double >( end - begin, 0 ) );
 std::vector< std::vector< double > > geomlinferrs_global( 1 + degs2fit.size(),
 std::vector< double >( end - begin, 0 ) );
#endif
 
 for( int i = begin; i < end; ++i )
 {
 std::ostringstream convert;
 convert << i;
 std::string infile = prefix + convert.str() + suffix;
 int ntestverts;
 std::vector< double > geoml1errs, geoml2errs, geomlinferrs;
 
#ifdef MOAB_HAVE_MPI
 std::cout << "Processor " << rank << " is working on file " << infile << std::endl;
#endif
 error = closedsurface_uref_hirec_convergence_study( infile.c_str(), degs2fit, interp, dim, obj, ntestverts,
 geoml1errs, geoml2errs, geomlinferrs );MB_CHK_ERR( error );
 assert( geoml1errs.size() == 1 + degs2fit.size() && geoml2errs.size() == 1 + degs2fit.size() &&
 geomlinferrs.size() == 1 + degs2fit.size() );
#ifdef MOAB_HAVE_MPI
 
 if( nprocs > 1 )
 {
 int ntestverts_global = 0;
 MPI_Reduce( &ntestverts, &ntestverts_global, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD );
 
 for( size_t d = 0; d < degs2fit.size() + 1; ++d )
 {
 double local_l1err = ntestverts * geoml1errs[d],
 local_l2err = geoml2errs[d] * ( geoml2errs[d] * ntestverts ), local_linferr = geomlinferrs[d];
 std::cout << "On Processor " << rank << " with mesh " << i
 << " Degree = " << ( d == 0 ? 0 : degs2fit[d - 1] ) << " L1:" << geoml1errs[d]
 << " L2:" << geoml2errs[d] << " Li:" << geomlinferrs[d] << std::endl;
 double global_l1err = 0, global_l2err = 0, global_linferr = 0;
 MPI_Reduce( &local_l1err, &global_l1err, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
 MPI_Reduce( &local_l2err, &global_l2err, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
 MPI_Reduce( &local_linferr, &global_linferr, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD );
 
 if( rank == 0 )
 {
 geoml1errs_global[d][i - begin] = global_l1err / ntestverts_global;
 geoml2errs_global[d][i - begin] = sqrt( global_l2err / ntestverts_global );
 geomlinferrs_global[d][i - begin] = global_linferr;
 }
 }
 }
 else
 {
 for( size_t d = 0; d < degs2fit.size() + 1; ++d )
 {
 geoml1errs_global[d][i - begin] = geoml1errs[d];
 geoml2errs_global[d][i - begin] = geoml2errs[d];
 geomlinferrs_global[d][i - begin] = geomlinferrs[d];
 }
 }
 
#else
 
 for( size_t d = 0; d < degs2fit.size() + 1; ++d )
 {
 geoml1errs_global[d][i - begin] = geoml1errs[d];
 geoml2errs_global[d][i - begin] = geoml2errs[d];
 geomlinferrs_global[d][i - begin] = geomlinferrs[d];
 }
 
#endif
 }
 
#ifdef MOAB_HAVE_MPI
 
 if( rank == 0 )
 {
 std::cout << "Degrees: 0 ";
 
 for( size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
 {
 std::cout << degs2fit[ideg] << " ";
 }
 
 std::cout << std::endl;
 std::cout << "L1-norm error: \n";
 
 for( size_t i = 0; i < geoml1errs_global.size(); ++i )
 {
 for( size_t j = 0; j < geoml1errs_global[i].size(); ++j )
 {
 std::cout << geoml1errs_global[i][j] << " ";
 }
 
 std::cout << std::endl;
 }
 
 std::cout << "L2-norm error: \n";
 
 for( size_t i = 0; i < geoml2errs_global.size(); ++i )
 {
 for( size_t j = 0; j < geoml2errs_global[i].size(); ++j )
 {
 std::cout << geoml2errs_global[i][j] << " ";
 }
 
 std::cout << std::endl;
 }
 
 std::cout << "Linf-norm error: \n";
 
 for( size_t i = 0; i < geomlinferrs_global.size(); ++i )
 {
 for( size_t j = 0; j < geomlinferrs_global[i].size(); ++j )
 {
 std::cout << geomlinferrs_global[i][j] << " ";
 }
 
 std::cout << std::endl;
 }
 }
 
#else
 std::cout << "Degrees: 0 ";
 
 for( size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
 {
 std::cout << degs2fit[ideg] << " ";
 }
 
 std::cout << std::endl;
 std::cout << "L1-norm error: \n";
 
 for( size_t i = 0; i < geoml1errs_global.size(); ++i )
 {
 for( size_t j = 0; j < geoml1errs_global[i].size(); ++j )
 {
 std::cout << geoml1errs_global[i][j] << " ";
 }
 
 std::cout << std::endl;
 }
 
 std::cout << "L2-norm error: \n";
 
 for( size_t i = 0; i < geoml2errs_global.size(); ++i )
 {
 for( size_t j = 0; j < geoml2errs_global[i].size(); ++j )
 {
 std::cout << geoml2errs_global[i][j] << " ";
 }
 
 std::cout << std::endl;
 }
 
 std::cout << "Linf-norm error: \n";
 
 for( size_t i = 0; i < geomlinferrs_global.size(); ++i )
 {
 for( size_t j = 0; j < geomlinferrs_global[i].size(); ++j )
 {
 std::cout << geomlinferrs_global[i][j] << " ";
 }
 
 std::cout << std::endl;
 }
 
#endif
#ifdef MOAB_HAVE_MPI
 MPI_Finalize();
#endif
}

References closedsurface_uref_hirec_convergence_study(), dim, moab::error(), ErrorCode, geom, MB_CHK_ERR, MPI_COMM_WORLD, rank, and usage().

usage()

void usage

(

)

Definition at line 412 of file uref_hirec_test_parallel.cpp.

{
 std::cout << "usage: mpirun -np <number of processors> ./uref_hirec_test_parallel <prefix of "
 "files> -str <first number> -end <lastnumber> -suffix <suffix> -interp <0=least "
 "square, 1=interpolation> -dim <mesh dimension> -geom <s=sphere,t=torus>"
 << std::endl;
}

Referenced by main().