MOAB: Mesh Oriented datABase  (version 5.2.1)
uref_hirec_test_parallel.cpp
Go to the documentation of this file.
00001 /*This unit test is for the convergence study of high order reconstruction under uniform
00002  * refinement*/
00003 #include <iostream>
00004 #include <string>
00005 #include <sstream>
00006 #include <sys/time.h>
00007 #include <vector>
00008 #include <algorithm>
00009 #include "moab/Core.hpp"
00010 #include "moab/Range.hpp"
00011 #include "moab/MeshTopoUtil.hpp"
00012 #include "moab/NestedRefine.hpp"
00013 #include "moab/DiscreteGeometry/HiReconstruction.hpp"
00014 #include "TestUtil.hpp"
00015 #include "geomObject.cpp"
00016 #include <math.h>
00017 #include <stdlib.h>
00018 #include <assert.h>
00019 
00020 #ifdef MOAB_HAVE_MPI
00021 #include "moab/ParallelComm.hpp"
00022 #include "MBParallelConventions.h"
00023 #include "ReadParallel.hpp"
00024 #include "moab/FileOptions.hpp"
00025 #include "MBTagConventions.hpp"
00026 #include "moab_mpi.h"
00027 #endif
00028 
00029 using namespace moab;
00030 
00031 #define nsamples 10
00032 
00033 #ifdef MOAB_HAVE_MPI
00034 std::string read_options;
00035 #endif
00036 
00037 ErrorCode load_meshset_hirec( const char* infile, Interface* mbimpl, EntityHandle& meshset, ParallelComm*& pc,
00038                               const int degree, const int dim )
00039 {
00040     ErrorCode error;
00041     error = mbimpl->create_meshset( moab::MESHSET_SET, meshset );MB_CHK_ERR( error );
00042 #ifdef MOAB_HAVE_MPI
00043     int nprocs, rank;
00044     MPI_Comm comm = MPI_COMM_WORLD;
00045     MPI_Comm_size( comm, &nprocs );
00046     MPI_Comm_rank( comm, &rank );
00047     EntityHandle partnset;
00048     error = mbimpl->create_meshset( moab::MESHSET_SET, partnset );MB_CHK_ERR( error );
00049 
00050     if( nprocs > 1 ) { pc = moab::ParallelComm::get_pcomm( mbimpl, partnset, &comm ); }
00051 
00052     if( nprocs > 1 )
00053     {
00054         int nghlayers = degree > 0 ? HiReconstruction::estimate_num_ghost_layers( degree, true ) : 0;
00055         assert( nghlayers < 10 );
00056 
00057         if( nghlayers )
00058         {
00059             // get ghost layers
00060             if( dim == 2 )
00061             {
00062                 read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
00063                                "SHARED_ENTS;PARALLEL_GHOSTS=2.0.";
00064             }
00065             else if( dim == 1 )
00066             {
00067                 read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
00068                                "SHARED_ENTS;PARALLEL_GHOSTS=1.0.";
00069             }
00070             else
00071             {
00072                 MB_SET_ERR( MB_FAILURE, "unsupported dimension" );
00073             }
00074 
00075             read_options += (char)( '0' + nghlayers );
00076             // std::cout << "On processor " << rank << " Degree=" << degree << " "<< read_options <<
00077             // std::endl;
00078         }
00079         else
00080         {
00081             read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
00082         }
00083 
00084         /*for debug*/
00085         /*std::string outfile = std::string(infile); std::size_t dotpos = outfile.find_last_of(".");
00086         std::ostringstream convert; convert << rank;
00087         std::string localfile = outfile.substr(0,dotpos)+convert.str()+".vtk";
00088         //write local mesh
00089         EntityHandle local;
00090         error = mbimpl->create_meshset(moab::MESHSET_SET,local);CHECK_ERR(error);
00091         std::string local_options =
00092         "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"; error =
00093         mbimpl->load_file(infile,&local,local_options.c_str()); MB_CHK_ERR(error); error =
00094         mbimpl->write_file(localfile.c_str(),0,0,&local,1);*/
00095 
00096         error = mbimpl->load_file( infile, &meshset, read_options.c_str() );MB_CHK_ERR( error );
00097         /*for debug
00098         //write local mesh with ghost layers
00099         localfile = outfile.substr(0,dotpos)+convert.str()+"_ghost.vtk";
00100         error = mbimpl->write_file(localfile.c_str(),0,0,&meshset,1);*/
00101     }
00102     else
00103     {
00104         error = mbimpl->load_file( infile, &meshset );MB_CHK_ERR( error );
00105     }
00106 
00107 #else
00108     assert( !pc && degree && dim );
00109     error = mbimpl->load_file( infile, &meshset );MB_CHK_ERR( error );
00110 #endif
00111     return error;
00112 }
00113 
00114 ErrorCode closedsurface_uref_hirec_convergence_study( const char* infile, std::vector< int >& degs2fit, bool interp,
00115                                                       int dim, geomObject* obj, int& ntestverts,
00116                                                       std::vector< double >& geoml1errs,
00117                                                       std::vector< double >& geoml2errs,
00118                                                       std::vector< double >& geomlinferrs )
00119 {
00120     Core moab;
00121     Interface* mbImpl = &moab;
00122     ParallelComm* pc  = NULL;
00123     EntityHandle meshset;
00124 
00125 #ifdef MOAB_HAVE_MPI
00126     int nprocs, rank;
00127     MPI_Comm comm = MPI_COMM_WORLD;
00128     MPI_Comm_size( comm, &nprocs );
00129     MPI_Comm_rank( comm, &rank );
00130 #endif
00131 
00132     ErrorCode error;
00133     // mesh will be loaded and communicator pc will be updated
00134     int mxdeg = 1;
00135     for( size_t i = 0; i < degs2fit.size(); ++i )
00136     {
00137         mxdeg = std::max( degs2fit[i], mxdeg );
00138     }
00139     error = load_meshset_hirec( infile, mbImpl, meshset, pc, mxdeg, dim );MB_CHK_ERR( error );
00140 
00141     Range elems, elems_owned;
00142     error = mbImpl->get_entities_by_dimension( meshset, dim, elems );MB_CHK_ERR( error );
00143 
00144 #ifdef MOAB_HAVE_MPI
00145 
00146     if( pc )
00147     {
00148         error = pc->filter_pstatus( elems, PSTATUS_GHOST, PSTATUS_NOT, -1, &elems_owned );MB_CHK_ERR( error );
00149     }
00150     else
00151     {
00152         elems_owned = elems;
00153     }
00154 
00155 #endif
00156 
00157 #ifdef MOAB_HAVE_MPI
00158     // std::cout << "Mesh has " << nelems << " elements on Processor " << rank << " in total;";
00159     // std::cout << elems_owned.size() << " of which are locally owned elements" << std::endl;
00160 #else
00161     // std::cout << "Mesh has " << nelems << " elements" << std::endl;
00162 #endif
00163 
00164     /************************
00165      *   convergence study   *
00166      *************************/
00167     // project onto exact geometry since each level with uref has only linear coordinates
00168     Range verts;
00169     error = mbImpl->get_entities_by_dimension( meshset, 0, verts );MB_CHK_ERR( error );
00170 
00171     for( Range::iterator ivert = verts.begin(); ivert != verts.end(); ++ivert )
00172     {
00173         EntityHandle currvert = *ivert;
00174         double currcoords[3], exactcoords[3];
00175         error = mbImpl->get_coords( &currvert, 1, currcoords );MB_CHK_ERR( error );
00176         obj->project_points2geom( 3, currcoords, exactcoords, NULL );
00177 
00178         error = mbImpl->set_coords( &currvert, 1, exactcoords );MB_CHK_ERR( error );
00179         // for debug
00180         /*error = mbImpl->get_coords(&currvert,1,currcoords); MB_CHK_ERR(error);
00181         assert(currcoords[0]==exactcoords[0]&&currcoords[1]==exactcoords[1]&&currcoords[2]==exactcoords[2]);*/
00182     }
00183 
00184     // test ahf on mesh with ghost layers
00185     /*Range verts_owned;
00186     #ifdef MOAB_HAVE_MPI
00187         if(pc){
00188             error =
00189     pc->filter_pstatus(verts,PSTATUS_GHOST,PSTATUS_NOT,-1,&verts_owned);MB_CHK_ERR(error); }else{
00190     verts_owned = verts;
00191         }
00192 
00193         HalfFacetRep *ahf = new HalfFacetRep(&moab,pc,meshset);
00194         error = ahf->initialize(); MB_CHK_ERR(error);
00195 
00196         for(int i=0;i<nprocs;++i){
00197             if(rank==i){
00198                 std::cout << "Processor " << rank << " Local elements: \n";
00199                 for(Range::iterator ielem=elems.begin();ielem!=elems.end();++ielem){
00200                     EntityHandle currelem = *ielem;
00201                     std::vector<EntityHandle> conn;
00202                     error = mbImpl->get_connectivity(&currelem,1,conn); MB_CHK_ERR(error);
00203                     std::cout << *ielem << ": ";
00204                     for(size_t k=0;k<conn.size();++k) std::cout << conn[k] << " ";
00205                     std::cout << std::endl;
00206                 }
00207 
00208                 std::cout << "Processor " << rank << " Local owned elements: \n";
00209                 for(Range::iterator ielem=elems_owned.begin();ielem!=elems_owned.end();++ielem){
00210                     EntityHandle currelem = *ielem;
00211                     std::vector<EntityHandle> conn;
00212                     error = mbImpl->get_connectivity(&currelem,1,conn); MB_CHK_ERR(error);
00213                     std::cout << *ielem << ": ";
00214                     for(size_t k=0;k<conn.size();++k) std::cout << conn[k] << " ";
00215                     std::cout << std::endl;
00216                 }
00217 
00218 
00219                 std::cout << "Processor " << rank << " Local vertices: ";
00220                 for(Range::iterator ivert=verts.begin();ivert!=verts.end();++ivert){
00221                     std::cout << *ivert << " ";
00222                 }
00223                 std::cout << std::endl;
00224                 std::cout << "Processor " << rank << " Local owned vertices: ";
00225                 for(Range::iterator ivert=verts_owned.begin();ivert!=verts_owned.end();++ivert){
00226                     std::cout << *ivert << " ";
00227                 }
00228                 std::cout << std::endl;
00229 
00230                 for(Range::iterator ivert=verts_owned.begin();ivert!=verts_owned.end();++ivert){
00231                     EntityHandle currvid = *ivert;
00232                     std::cout << "Processor " << rank << " local verts: " << *ivert << " has
00233     adjfaces: "; std::vector<EntityHandle> adjfaces;
00234                     //error = ahf->get_up_adjacencies(currvid,2,adjfaces); MB_CHK_ERR(error);
00235                     error = mbImpl->get_adjacencies(&currvid,1,2,false,adjfaces); MB_CHK_ERR(error);
00236                     for(size_t k=0;k<adjfaces.size();++k){
00237                         std::cout << adjfaces[k] << " ";
00238                     }
00239                     std::cout << std::endl;
00240                 }
00241 
00242                 for(Range::iterator ivert=verts.begin();ivert!=verts.end();++ivert){
00243                     EntityHandle currvid = *ivert;
00244                     std::cout << "Processor " << rank << " all verts: " << *ivert << " has adjfaces:
00245     "; std::vector<EntityHandle> adjfaces;
00246                     //error = ahf->get_up_adjacencies(*ivert,2,adjfaces); MB_CHK_ERR(error);
00247                     error = mbImpl->get_adjacencies(&currvid,1,2,false,adjfaces); MB_CHK_ERR(error);
00248                     for(size_t k=0;k<adjfaces.size();++k){
00249                         std::cout << adjfaces[k] << " ";
00250                     }
00251                     std::cout << std::endl;
00252                 }
00253 
00254             }else{
00255                 MPI_Barrier(MPI_COMM_WORLD);
00256             }
00257         }
00258         exit(0);
00259     #endif*/
00260 
00261     // generate random points on each elements, assument 3D coordinates
00262     int nvpe = TYPE_FROM_HANDLE( *elems.begin() ) == MBTRI ? 3 : 4;
00263     std::vector< double > testpnts, testnaturalcoords;
00264     ntestverts = elems_owned.size() * nsamples;
00265     testpnts.reserve( 3 * elems_owned.size() * nsamples );
00266     testnaturalcoords.reserve( nvpe * elems_owned.size() * nsamples );
00267 
00268     for( Range::iterator ielem = elems_owned.begin(); ielem != elems_owned.end(); ++ielem )
00269     {
00270         EntityHandle currelem = *ielem;
00271         std::vector< EntityHandle > conn;
00272         error = mbImpl->get_connectivity( &currelem, 1, conn );MB_CHK_ERR( error );
00273         std::vector< double > elemcoords( 3 * conn.size() );
00274         error = mbImpl->get_coords( &( conn[0] ), conn.size(), &( elemcoords[0] ) );MB_CHK_ERR( error );
00275         EntityType type = TYPE_FROM_HANDLE( currelem );
00276 
00277         for( int s = 0; s < nsamples; ++s )
00278         {
00279             if( type == MBTRI )
00280             {
00281                 double a = (double)rand() / RAND_MAX, b = (double)rand() / RAND_MAX, c = (double)rand() / RAND_MAX, sum;
00282                 sum = a + b + c;
00283 
00284                 if( sum < 1e-12 )
00285                 {
00286                     --s;
00287                     continue;
00288                 }
00289                 else
00290                 {
00291                     a /= sum, b /= sum, c /= sum;
00292                 }
00293 
00294                 assert( a > 0 && b > 0 && c > 0 && fabs( a + b + c - 1 ) < 1e-12 );
00295                 testpnts.push_back( a * elemcoords[0] + b * elemcoords[3] + c * elemcoords[6] );
00296                 testpnts.push_back( a * elemcoords[1] + b * elemcoords[4] + c * elemcoords[7] );
00297                 testpnts.push_back( a * elemcoords[2] + b * elemcoords[5] + c * elemcoords[8] );
00298                 testnaturalcoords.push_back( a ), testnaturalcoords.push_back( b ), testnaturalcoords.push_back( c );
00299             }
00300             else if( type == MBQUAD )
00301             {
00302                 double xi = (double)rand() / RAND_MAX, eta = (double)rand() / RAND_MAX, N[4];
00303                 xi   = 2 * xi - 1;
00304                 eta  = 2 * eta - 1;
00305                 N[0] = ( 1 - xi ) * ( 1 - eta ) / 4, N[1] = ( 1 + xi ) * ( 1 - eta ) / 4,
00306                 N[2] = ( 1 + xi ) * ( 1 + eta ) / 4, N[3] = ( 1 - xi ) * ( 1 + eta ) / 4;
00307                 testpnts.push_back( N[0] * elemcoords[0] + N[1] * elemcoords[3] + N[2] * elemcoords[6] +
00308                                     N[3] * elemcoords[9] );
00309                 testpnts.push_back( N[0] * elemcoords[1] + N[1] * elemcoords[4] + N[2] * elemcoords[7] +
00310                                     N[3] * elemcoords[10] );
00311                 testpnts.push_back( N[0] * elemcoords[2] + N[1] * elemcoords[5] + N[2] * elemcoords[8] +
00312                                     N[3] * elemcoords[11] );
00313                 testnaturalcoords.push_back( N[0] ), testnaturalcoords.push_back( N[1] ),
00314                     testnaturalcoords.push_back( N[2] ), testnaturalcoords.push_back( N[3] );
00315             }
00316             else
00317             {
00318                 throw std::invalid_argument( "NOT SUPPORTED ELEMENT TYPE" );
00319             }
00320         }
00321     }
00322 
00323     // Compute error of linear interpolation in PARALLEL
00324     double l1err, l2err, linferr;
00325     geoml1errs.clear();
00326     geoml2errs.clear();
00327     geomlinferrs.clear();
00328     obj->compute_projecterror( 3, elems_owned.size() * nsamples, &( testpnts[0] ), l1err, l2err, linferr );
00329     geoml1errs.push_back( l1err );
00330     geoml2errs.push_back( l2err );
00331     geomlinferrs.push_back( linferr );
00332     /*Perform high order projection and compute error*/
00333     // initialize
00334     HiReconstruction hirec( &moab, pc, meshset );
00335 
00336     for( size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
00337     {
00338         // High order reconstruction
00339         error = hirec.reconstruct3D_surf_geom( degs2fit[ideg], interp, true, true );MB_CHK_ERR( error );
00340 
00341         int index = 0;
00342         // for debug
00343         /*double maxlinferr=0,elemlinferr=0,eleml1err=0,eleml2err;
00344         EntityHandle elem_maxerr;*/
00345 
00346         for( Range::iterator ielem = elems_owned.begin(); ielem != elems_owned.end(); ++ielem, ++index )
00347         {
00348             // Projection
00349             error =
00350                 hirec.hiproj_walf_in_element( *ielem, nvpe, nsamples, &( testnaturalcoords[nvpe * nsamples * index] ),
00351                                               &( testpnts[3 * nsamples * index] ) );MB_CHK_ERR( error );
00352             // for debug
00353             /*obj->compute_projecterror(3,nsamples,&(testpnts[3*nsamples*index]),eleml1err,eleml2err,elemlinferr);
00354             if(elemlinferr>maxlinferr){
00355                 maxlinferr = elemlinferr; elem_maxerr = *ielem;
00356             }*/
00357         }
00358 
00359         assert( (size_t)index == elems_owned.size() );
00360         // Estimate error
00361         obj->compute_projecterror( 3, elems_owned.size() * nsamples, &( testpnts[0] ), l1err, l2err, linferr );
00362         geoml1errs.push_back( l1err );
00363         geoml2errs.push_back( l2err );
00364         geomlinferrs.push_back( linferr );
00365 
00366         // for debug
00367         /*std::cout << "Degree " << degs2fit[ideg] << " has L-inf error " << maxlinferr << " at
00368         element " << elems_owned.index(elem_maxerr) << ":"; std::vector<EntityHandle> conn; error =
00369         mbImpl->get_connectivity(&elem_maxerr,1,conn); MB_CHK_ERR(error); for(size_t
00370         ii=0;ii<conn.size();++ii) std::cout << verts.index(conn[ii]) << " "; std::cout << std::endl;
00371 
00372         for(size_t ii=0;ii<conn.size();++ii){
00373             //assume triangle
00374             std::vector<double> vertexcoords(3,0);
00375             error = mbImpl->get_coords(&(conn[ii]),1,&(vertexcoords[0])); MB_CHK_ERR(error);
00376             std::cout << verts.index(conn[ii]) << ":" << vertexcoords[0] << " "<< vertexcoords[1] <<
00377         " "<< vertexcoords[2] << "\n"; GEOMTYPE geomtype; std::vector<double> local_coords_system,
00378         local_coeffs; int deg_out; bool local_interp; bool hasfit =
00379         hirec.get_fittings_data(conn[ii],geomtype,local_coords_system,deg_out,local_coeffs,local_interp);
00380         assert(hasfit); std::cout << "Local fitting result: " << std::endl; std::cout << "Geom Type:
00381         " << geomtype << std::endl; std::cout << "Local coordinates system: " << std::endl;
00382             std::cout << local_coords_system[0] << "\t" << local_coords_system[1] << "\t" <<
00383         local_coords_system[2] << std::endl; std::cout << local_coords_system[3] << "\t" <<
00384         local_coords_system[4] << "\t" << local_coords_system[5] << std::endl; std::cout <<
00385         local_coords_system[6] << "\t" << local_coords_system[7] << "\t" << local_coords_system[8]
00386         << std::endl; std::cout << "Fitting degree is " << deg_out << " interp is " << local_interp
00387         << std::endl; std::cout << "Fitting coefficients are :" << std::endl; for(size_t
00388         kk=0;kk<local_coeffs.size();++kk) std::cout << local_coeffs[kk] << "\t"; std::cout <<
00389         std::endl;
00390         }*/
00391     }
00392 
00393     return error;
00394 }
00395 
00396 void usage()
00397 {
00398     std::cout << "usage: mpirun -np <number of processors> ./uref_hirec_test_parallel <prefix of "
00399                  "files> -str <first number> -end <lastnumber> -suffix <suffix> -interp <0=least "
00400                  "square, 1=interpolation> -dim <mesh dimension> -geom <s=sphere,t=torus>"
00401               << std::endl;
00402 }
00403 
00404 int main( int argc, char* argv[] )
00405 {
00406 #ifdef MOAB_HAVE_MPI
00407     MPI_Init( &argc, &argv );
00408     int nprocs, rank;
00409     MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
00410     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00411 #endif
00412     std::string prefix, suffix, istr, iend;
00413     int dim = 2, geom = 0;
00414     bool interp = false;
00415     ErrorCode error;
00416 
00417     if( argc == 1 )
00418     {
00419         usage();
00420         prefix = TestDir + "/sphere_quads_5_ML_";
00421         suffix = ".h5m";
00422         istr   = "0";
00423         iend   = "1";
00424         std::cout << "Using defaults: MeshFiles/unittest/sphere_quads_5_ML_*.h5m -str 0 -end 1 "
00425                      "-suffix \".h5m\" -interp 0 -dim 2 -geom s"
00426                   << std::endl;
00427     }
00428     else
00429     {
00430         prefix      = std::string( argv[1] );
00431         bool hasdim = false;
00432 
00433         for( int i = 2; i < argc; ++i )
00434         {
00435             if( i + 1 != argc )
00436             {
00437                 if( std::string( argv[i] ) == "-suffix" ) { suffix = std::string( argv[++i] ); }
00438                 else if( std::string( argv[i] ) == "-str" )
00439                 {
00440                     istr = std::string( argv[++i] );
00441                 }
00442                 else if( std::string( argv[i] ) == "-end" )
00443                 {
00444                     iend = std::string( argv[++i] );
00445                 }
00446                 else if( std::string( argv[i] ) == "-interp" )
00447                 {
00448                     interp = atoi( argv[++i] );
00449                 }
00450                 else if( std::string( argv[i] ) == "-dim" )
00451                 {
00452                     dim    = atoi( argv[++i] );
00453                     hasdim = true;
00454                 }
00455                 else if( std::string( argv[i] ) == "-geom" )
00456                 {
00457                     geom = std::string( argv[++i] ) == "s" ? 0 : 1;
00458                 }
00459                 else
00460                 {
00461 #ifdef MOAB_HAVE_MPI
00462 
00463                     if( 0 == rank ) { usage(); }
00464 
00465                     MPI_Finalize();
00466 #else
00467                     usage();
00468 #endif
00469                     return 0;
00470                 }
00471             }
00472         }
00473 
00474         if( !hasdim )
00475         {
00476 #ifdef MOAB_HAVE_MPI
00477 
00478             if( 0 == rank )
00479             { std::cout << "Dimension of input mesh should be provided, positive and less than 3" << std::endl; }
00480 
00481 #else
00482             std::cout << "Dimension of input mesh should be provided, positive and less than 3" << std::endl;
00483 #endif
00484             return 0;
00485         }
00486 
00487         if( dim > 2 || dim <= 0 )
00488         {
00489 #ifdef MOAB_HAVE_MPI
00490 
00491             if( 0 == rank ) { std::cout << "Input dimesion should be positive and less than 3;" << std::endl; }
00492 
00493 #else
00494             std::cout << "Input dimesion should be positive and less than 3;" << std::endl;
00495 #endif
00496             return 0;
00497         }
00498 
00499 #ifdef MOAB_HAVE_MPI
00500 
00501         if( 0 == rank )
00502         {
00503             std::cout << "Testing on " << prefix + istr + "-" + iend + suffix << " with dimension " << dim << "\n";
00504             std::string opts = interp ? "interpolation" : "least square fitting";
00505             std::cout << "High order reconstruction in " << opts << std::endl;
00506         }
00507 
00508 #else
00509         std::cout << "Testing on " << prefix + istr + "-" + iend + suffix << " with dimension " << dim << "\n";
00510         std::string opts = interp ? "interpolation" : "least square fitting";
00511         std::cout << "High order reconstruction in " << opts << std::endl;
00512 #endif
00513     }
00514 
00515     geomObject* obj;
00516 
00517     if( geom == 0 ) { obj = new sphere(); }
00518     else
00519     {
00520         obj = new torus();
00521     }
00522 
00523     std::vector< int > degs2fit;
00524     for( int d = 1; d <= 6; ++d )
00525     {
00526         degs2fit.push_back( d );
00527     }
00528 
00529     int begin = atoi( istr.c_str() ), end = atoi( iend.c_str() ) + 1;
00530 #ifdef MOAB_HAVE_MPI
00531     std::vector< std::vector< double > > geoml1errs_global, geoml2errs_global, geomlinferrs_global;
00532 
00533     if( rank == 0 )
00534     {
00535         geoml1errs_global =
00536             std::vector< std::vector< double > >( 1 + degs2fit.size(), std::vector< double >( end - begin, 0 ) );
00537         geoml2errs_global =
00538             std::vector< std::vector< double > >( 1 + degs2fit.size(), std::vector< double >( end - begin, 0 ) );
00539         geomlinferrs_global =
00540             std::vector< std::vector< double > >( 1 + degs2fit.size(), std::vector< double >( end - begin, 0 ) );
00541     }
00542 
00543 #else
00544     std::vector< std::vector< double > > geoml1errs_global( 1 + degs2fit.size(),
00545                                                             std::vector< double >( end - begin, 0 ) );
00546     std::vector< std::vector< double > > geoml2errs_global( 1 + degs2fit.size(),
00547                                                             std::vector< double >( end - begin, 0 ) );
00548     std::vector< std::vector< double > > geomlinferrs_global( 1 + degs2fit.size(),
00549                                                               std::vector< double >( end - begin, 0 ) );
00550 #endif
00551 
00552     for( int i = begin; i < end; ++i )
00553     {
00554         std::ostringstream convert;
00555         convert << i;
00556         std::string infile = prefix + convert.str() + suffix;
00557         int ntestverts;
00558         std::vector< double > geoml1errs, geoml2errs, geomlinferrs;
00559 
00560 #ifdef MOAB_HAVE_MPI
00561         std::cout << "Processor " << rank << " is working on file " << infile << std::endl;
00562 #endif
00563         error = closedsurface_uref_hirec_convergence_study( infile.c_str(), degs2fit, interp, dim, obj, ntestverts,
00564                                                             geoml1errs, geoml2errs, geomlinferrs );MB_CHK_ERR( error );
00565         assert( geoml1errs.size() == 1 + degs2fit.size() && geoml2errs.size() == 1 + degs2fit.size() &&
00566                 geomlinferrs.size() == 1 + degs2fit.size() );
00567 #ifdef MOAB_HAVE_MPI
00568 
00569         if( nprocs > 1 )
00570         {
00571             int ntestverts_global = 0;
00572             MPI_Reduce( &ntestverts, &ntestverts_global, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD );
00573 
00574             for( size_t d = 0; d < degs2fit.size() + 1; ++d )
00575             {
00576                 double local_l1err = ntestverts * geoml1errs[d],
00577                        local_l2err = geoml2errs[d] * ( geoml2errs[d] * ntestverts ), local_linferr = geomlinferrs[d];
00578                 std::cout << "On Processor " << rank << " with mesh " << i
00579                           << " Degree = " << ( d == 0 ? 0 : degs2fit[d - 1] ) << " L1:" << geoml1errs[d]
00580                           << " L2:" << geoml2errs[d] << " Li:" << geomlinferrs[d] << std::endl;
00581                 double global_l1err = 0, global_l2err = 0, global_linferr = 0;
00582                 MPI_Reduce( &local_l1err, &global_l1err, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
00583                 MPI_Reduce( &local_l2err, &global_l2err, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
00584                 MPI_Reduce( &local_linferr, &global_linferr, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD );
00585 
00586                 if( rank == 0 )
00587                 {
00588                     geoml1errs_global[d][i - begin]   = global_l1err / ntestverts_global;
00589                     geoml2errs_global[d][i - begin]   = sqrt( global_l2err / ntestverts_global );
00590                     geomlinferrs_global[d][i - begin] = global_linferr;
00591                 }
00592             }
00593         }
00594         else
00595         {
00596             for( size_t d = 0; d < degs2fit.size() + 1; ++d )
00597             {
00598                 geoml1errs_global[d][i - begin]   = geoml1errs[d];
00599                 geoml2errs_global[d][i - begin]   = geoml2errs[d];
00600                 geomlinferrs_global[d][i - begin] = geomlinferrs[d];
00601             }
00602         }
00603 
00604 #else
00605 
00606         for( size_t d = 0; d < degs2fit.size() + 1; ++d )
00607         {
00608             geoml1errs_global[d][i - begin]   = geoml1errs[d];
00609             geoml2errs_global[d][i - begin]   = geoml2errs[d];
00610             geomlinferrs_global[d][i - begin] = geomlinferrs[d];
00611         }
00612 
00613 #endif
00614     }
00615 
00616 #ifdef MOAB_HAVE_MPI
00617 
00618     if( rank == 0 )
00619     {
00620         std::cout << "Degrees: 0 ";
00621 
00622         for( size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
00623         {
00624             std::cout << degs2fit[ideg] << " ";
00625         }
00626 
00627         std::cout << std::endl;
00628         std::cout << "L1-norm error: \n";
00629 
00630         for( size_t i = 0; i < geoml1errs_global.size(); ++i )
00631         {
00632             for( size_t j = 0; j < geoml1errs_global[i].size(); ++j )
00633             {
00634                 std::cout << geoml1errs_global[i][j] << " ";
00635             }
00636 
00637             std::cout << std::endl;
00638         }
00639 
00640         std::cout << "L2-norm error: \n";
00641 
00642         for( size_t i = 0; i < geoml2errs_global.size(); ++i )
00643         {
00644             for( size_t j = 0; j < geoml2errs_global[i].size(); ++j )
00645             {
00646                 std::cout << geoml2errs_global[i][j] << " ";
00647             }
00648 
00649             std::cout << std::endl;
00650         }
00651 
00652         std::cout << "Linf-norm error: \n";
00653 
00654         for( size_t i = 0; i < geomlinferrs_global.size(); ++i )
00655         {
00656             for( size_t j = 0; j < geomlinferrs_global[i].size(); ++j )
00657             {
00658                 std::cout << geomlinferrs_global[i][j] << " ";
00659             }
00660 
00661             std::cout << std::endl;
00662         }
00663     }
00664 
00665 #else
00666     std::cout << "Degrees: 0 ";
00667 
00668     for( size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
00669     {
00670         std::cout << degs2fit[ideg] << " ";
00671     }
00672 
00673     std::cout << std::endl;
00674     std::cout << "L1-norm error: \n";
00675 
00676     for( size_t i = 0; i < geoml1errs_global.size(); ++i )
00677     {
00678         for( size_t j = 0; j < geoml1errs_global[i].size(); ++j )
00679         {
00680             std::cout << geoml1errs_global[i][j] << " ";
00681         }
00682 
00683         std::cout << std::endl;
00684     }
00685 
00686     std::cout << "L2-norm error: \n";
00687 
00688     for( size_t i = 0; i < geoml2errs_global.size(); ++i )
00689     {
00690         for( size_t j = 0; j < geoml2errs_global[i].size(); ++j )
00691         {
00692             std::cout << geoml2errs_global[i][j] << " ";
00693         }
00694 
00695         std::cout << std::endl;
00696     }
00697 
00698     std::cout << "Linf-norm error: \n";
00699 
00700     for( size_t i = 0; i < geomlinferrs_global.size(); ++i )
00701     {
00702         for( size_t j = 0; j < geomlinferrs_global[i].size(); ++j )
00703         {
00704             std::cout << geomlinferrs_global[i][j] << " ";
00705         }
00706 
00707         std::cout << std::endl;
00708     }
00709 
00710 #endif
00711 #ifdef MOAB_HAVE_MPI
00712     MPI_Finalize();
00713 #endif
00714 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines