|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#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.
Defines | |
| #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 > °s2fit, 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[]) |
| #define nsamples 10 |
Definition at line 36 of file uref_hirec_test_parallel.cpp.
Referenced by 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.
References moab::Range::begin(), geomObject::compute_projecterror(), 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().
{
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;
}
| 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.
References moab::Interface::create_meshset(), 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.
{
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;
}
| int main | ( | int | argc, |
| char * | argv[] | ||
| ) |
Definition at line 420 of file uref_hirec_test_parallel.cpp.
References closedsurface_uref_hirec_convergence_study(), dim, moab::error(), ErrorCode, geom, MB_CHK_ERR, MPI_COMM_WORLD, rank, and usage.
{
#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
}
| 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;
}
1.7.6.1