MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 }