MOAB: Mesh Oriented datABase  (version 5.3.1)
uref_hirec_test.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 #if defined( __MINGW32__ )
00007 #include <sys/time.h>
00008 #else
00009 #include <ctime>
00010 #endif
00011 
00012 #include <vector>
00013 #include <algorithm>
00014 #include "moab/Core.hpp"
00015 #include "moab/Range.hpp"
00016 #include "moab/MeshTopoUtil.hpp"
00017 #include "moab/NestedRefine.hpp"
00018 #include "moab/DiscreteGeometry/DGMSolver.hpp"
00019 #include "moab/DiscreteGeometry/HiReconstruction.hpp"
00020 #include "TestUtil.hpp"
00021 #include "geomObject.cpp"
00022 #include <cmath>
00023 #include <cstdlib>
00024 
00025 #ifdef MOAB_HAVE_MPI
00026 #include "moab/ParallelComm.hpp"
00027 #include "MBParallelConventions.h"
00028 #include "ReadParallel.hpp"
00029 #include "moab/FileOptions.hpp"
00030 #include "MBTagConventions.hpp"
00031 #include "moab_mpi.h"
00032 #endif
00033 
00034 using namespace moab;
00035 
00036 #define nsamples 10
00037 
00038 #ifdef MOAB_HAVE_MPI
00039 std::string read_options;
00040 #endif
00041 
00042 ErrorCode test_closedsurface_mesh( const char* filename,
00043                                    int* level_degrees,
00044                                    int num_levels,
00045                                    int degree,
00046                                    bool interp,
00047                                    int dim,
00048                                    geomObject* obj )
00049 {
00050     Core moab;
00051     Interface* mbImpl = &moab;
00052     ParallelComm* pc  = NULL;
00053     EntityHandle fileset;
00054     ErrorCode error;
00055     error = mbImpl->create_meshset( moab::MESHSET_SET, fileset );MB_CHK_ERR( error );
00056     // load mesh from file
00057 #ifdef MOAB_HAVE_MPI
00058     MPI_Comm comm = MPI_COMM_WORLD;
00059     EntityHandle partnset;
00060     error = mbImpl->create_meshset( moab::MESHSET_SET, partnset );MB_CHK_ERR( error );
00061     pc        = moab::ParallelComm::get_pcomm( mbImpl, partnset, &comm );
00062     int procs = 1;
00063     MPI_Comm_size( comm, &procs );
00064 
00065     if( procs > 1 )
00066     {
00067         read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
00068         error        = mbImpl->load_file( filename, &fileset, read_options.c_str() );MB_CHK_ERR( error );
00069     }
00070     else
00071     {
00072         error = mbImpl->load_file( filename, &fileset );MB_CHK_ERR( error );
00073     }
00074 
00075 #else
00076     error = mbImpl->load_file( filename, &fileset );MB_CHK_ERR( error );
00077 #endif
00078 
00079     // Generate hierarchy
00080     // error = refine_entities(&moab, pc, fileset, level_degrees, num_levels, true);
00081     // CHECK_ERR(error);
00082     NestedRefine uref( &moab, pc, fileset );
00083     std::vector< EntityHandle > meshes;
00084     error = uref.generate_mesh_hierarchy( num_levels, level_degrees, meshes );MB_CHK_ERR( error );
00085 
00086     // Perform high order reconstruction on level 0 mesh
00087     HiReconstruction hirec( &moab, pc, fileset );
00088     assert( dim == 2 );
00089     error = hirec.reconstruct3D_surf_geom( degree, interp, false );MB_CHK_ERR( error );
00090 
00091     // High order projection
00092     Range verts, vorig;
00093     error = mbImpl->get_entities_by_dimension( meshes.back(), 0, verts );MB_CHK_ERR( error );
00094     error = mbImpl->get_entities_by_dimension( meshes.front(), 0, vorig );MB_CHK_ERR( error );
00095     int nvorig   = vorig.size();
00096     double l1err = 0, l2err = 0, linferr = 0;
00097 
00098     for( Range::iterator ivert = verts.begin() + nvorig; ivert != verts.end(); ++ivert )
00099     {
00100         // locate the element in level 0 mesh, on which *ivert is lying
00101         EntityHandle currvert = *ivert;
00102 
00103         std::vector< EntityHandle > parentEntities;
00104         error = uref.get_adjacencies( *ivert, 2, parentEntities );MB_CHK_ERR( error );
00105         assert( parentEntities.size() );
00106 
00107         EntityHandle rootelem;
00108         error = uref.child_to_parent( parentEntities[0], num_levels, 0, &rootelem );MB_CHK_ERR( error );
00109 
00110         // compute the natural coordinates of *ivert in this element
00111         assert( TYPE_FROM_HANDLE( rootelem ) == MBTRI );
00112 
00113         const EntityHandle* conn;
00114         int nvpe = 3;
00115         error    = mbImpl->get_connectivity( rootelem, conn, nvpe );MB_CHK_ERR( error );
00116 
00117         std::vector< double > cornercoords( 3 * nvpe ), currcoords( 3 );
00118         error = mbImpl->get_coords( conn, nvpe, &( cornercoords[0] ) );MB_CHK_ERR( error );
00119         error = mbImpl->get_coords( &currvert, 1, &( currcoords[0] ) );MB_CHK_ERR( error );
00120 
00121         std::vector< double > naturalcoords2fit( 3 );
00122         DGMSolver::get_tri_natural_coords( 3, &( cornercoords[0] ), 1, &( currcoords[0] ), &( naturalcoords2fit[0] ) );
00123 
00124         // project onto the estimated geometry
00125         double hicoords[3];
00126         error = hirec.hiproj_walf_in_element( rootelem, nvpe, 1, &( naturalcoords2fit[0] ), hicoords );MB_CHK_ERR( error );
00127 
00128         // estimate error
00129         double err = obj->compute_projecterror( 3, hicoords );
00130         l1err += err;
00131         l2err += err * err;
00132         linferr = std::max( linferr, err );
00133     }
00134 
00135     l1err /= verts.size() - nvorig;
00136     l2err = sqrt( l2err / ( verts.size() - nvorig ) );
00137     std::cout << "L1 error " << l1err << " L2 error " << l2err << " Linf error " << linferr << std::endl;
00138     return error;
00139 }
00140 
00141 ErrorCode closedsurface_uref_hirec_convergence_study( const char* filename,
00142                                                       int* level_degrees,
00143                                                       int num_levels,
00144                                                       std::vector< int >& degs2fit,
00145                                                       bool interp,
00146                                                       geomObject* obj )
00147 {
00148     Core moab;
00149     Interface* mbImpl = &moab;
00150     ParallelComm* pc  = NULL;
00151     EntityHandle fileset;
00152     ErrorCode error;
00153     error = mbImpl->create_meshset( moab::MESHSET_SET, fileset );MB_CHK_ERR( error );
00154     // load mesh from file
00155 #ifdef MOAB_HAVE_MPI
00156     MPI_Comm comm = MPI_COMM_WORLD;
00157     EntityHandle partnset;
00158     error = mbImpl->create_meshset( moab::MESHSET_SET, partnset );MB_CHK_ERR( error );
00159     pc        = moab::ParallelComm::get_pcomm( mbImpl, partnset, &comm );
00160     int procs = 1;
00161     MPI_Comm_size( comm, &procs );
00162 
00163     if( procs > 1 )
00164     {
00165         read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
00166         error        = mbImpl->load_file( filename, &fileset, read_options.c_str() );MB_CHK_ERR( error );
00167     }
00168     else
00169     {
00170         error = mbImpl->load_file( filename, &fileset );MB_CHK_ERR( error );
00171     }
00172 
00173 #else
00174     error = mbImpl->load_file( filename, &fileset );MB_CHK_ERR( error );
00175 #endif
00176     // Generate hierarchy
00177     NestedRefine uref( &moab, pc, fileset );
00178     std::vector< EntityHandle > meshes;
00179     std::vector< int > meshsizes;
00180     error = uref.generate_mesh_hierarchy( num_levels, level_degrees, meshes );MB_CHK_ERR( error );
00181     std::vector< std::vector< double > > geoml1errs( 1 + degs2fit.size() ), geoml2errs( 1 + degs2fit.size() ),
00182         geomlinferrs( 1 + degs2fit.size() );
00183 
00184     // Perform high order reconstruction on each level of mesh (projected onto exact geometry) and
00185     // estimate geometric error with various degrees of fitting
00186     for( size_t i = 0; i < meshes.size(); ++i )
00187     {
00188         EntityHandle& mesh = meshes[i];
00189         // project onto exact geometry since each level with uref has only linear coordinates
00190         Range verts;
00191         error = mbImpl->get_entities_by_dimension( mesh, 0, verts );MB_CHK_ERR( error );
00192 
00193         for( Range::iterator ivert = verts.begin(); ivert != verts.end(); ++ivert )
00194         {
00195             EntityHandle currvert = *ivert;
00196             double currcoords[3], exactcoords[3];
00197             error = mbImpl->get_coords( &currvert, 1, currcoords );
00198             obj->project_points2geom( 3, currcoords, exactcoords, NULL );
00199             error = mbImpl->set_coords( &currvert, 1, exactcoords );
00200         }
00201 
00202         // generate random points on each elements, assument 3D coordinates
00203         Range elems;
00204         error = mbImpl->get_entities_by_dimension( mesh, 2, elems );MB_CHK_ERR( error );
00205         meshsizes.push_back( elems.size() );
00206         int nvpe = TYPE_FROM_HANDLE( *elems.begin() ) == MBTRI ? 3 : 4;
00207         std::vector< double > testpnts, testnaturalcoords;
00208         testpnts.reserve( 3 * elems.size() * nsamples );
00209         testnaturalcoords.reserve( nvpe * elems.size() * nsamples );
00210 
00211         for( Range::iterator ielem = elems.begin(); ielem != elems.end(); ++ielem )
00212         {
00213             EntityHandle currelem = *ielem;
00214             std::vector< EntityHandle > conn;
00215             error = mbImpl->get_connectivity( &currelem, 1, conn );MB_CHK_ERR( error );
00216             nvpe = conn.size();
00217             std::vector< double > elemcoords( 3 * conn.size() );
00218             error = mbImpl->get_coords( &( conn[0] ), conn.size(), &( elemcoords[0] ) );MB_CHK_ERR( error );
00219             EntityType type = TYPE_FROM_HANDLE( currelem );
00220 
00221             for( int s = 0; s < nsamples; ++s )
00222             {
00223                 if( type == MBTRI )
00224                 {
00225                     double a = (double)rand() / RAND_MAX, b = (double)rand() / RAND_MAX, c = (double)rand() / RAND_MAX,
00226                            sum;
00227                     sum = a + b + c;
00228 
00229                     if( sum < 1e-12 )
00230                     {
00231                         --s;
00232                         continue;
00233                     }
00234                     else
00235                     {
00236                         a /= sum, b /= sum, c /= sum;
00237                     }
00238 
00239                     testpnts.push_back( a * elemcoords[0] + b * elemcoords[3] + c * elemcoords[6] );
00240                     testpnts.push_back( a * elemcoords[1] + b * elemcoords[4] + c * elemcoords[7] );
00241                     testpnts.push_back( a * elemcoords[2] + b * elemcoords[5] + c * elemcoords[8] );
00242                     testnaturalcoords.push_back( a );
00243                     testnaturalcoords.push_back( b );
00244                     testnaturalcoords.push_back( c );
00245                 }
00246                 else if( type == MBQUAD )
00247                 {
00248                     double xi = (double)rand() / RAND_MAX, eta = (double)rand() / RAND_MAX, N[4];
00249                     xi   = 2 * xi - 1;
00250                     eta  = 2 * eta - 1;
00251                     N[0] = ( 1 - xi ) * ( 1 - eta ) / 4, N[1] = ( 1 + xi ) * ( 1 - eta ) / 4,
00252                     N[2] = ( 1 + xi ) * ( 1 + eta ) / 4, N[3] = ( 1 - xi ) * ( 1 + eta ) / 4;
00253                     testpnts.push_back( N[0] * elemcoords[0] + N[1] * elemcoords[3] + N[2] * elemcoords[6] +
00254                                         N[3] * elemcoords[9] );
00255                     testpnts.push_back( N[0] * elemcoords[1] + N[1] * elemcoords[4] + N[2] * elemcoords[7] +
00256                                         N[3] * elemcoords[10] );
00257                     testpnts.push_back( N[0] * elemcoords[2] + N[1] * elemcoords[5] + N[2] * elemcoords[8] +
00258                                         N[3] * elemcoords[11] );
00259                     testnaturalcoords.push_back( N[0] );
00260                     testnaturalcoords.push_back( N[1] );
00261                     testnaturalcoords.push_back( N[2] );
00262                     testnaturalcoords.push_back( N[3] );
00263                 }
00264                 else
00265                 {
00266                     throw std::invalid_argument( "NOT SUPPORTED ELEMENT TYPE" );
00267                 }
00268             }
00269         }
00270 
00271         // Compute error of linear interpolation
00272         double l1err, l2err, linferr;
00273         obj->compute_projecterror( 3, elems.size() * nsamples, &( testpnts[0] ), l1err, l2err, linferr );
00274         geoml1errs[0].push_back( l1err );
00275         geoml2errs[0].push_back( l2err );
00276         geomlinferrs[0].push_back( linferr );
00277         // Perform high order projection and compute error
00278         HiReconstruction hirec( &moab, pc, mesh );
00279 
00280         for( size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
00281         {
00282             // High order reconstruction
00283             error = hirec.reconstruct3D_surf_geom( degs2fit[ideg], interp, false, true );MB_CHK_ERR( error );
00284             int index = 0;
00285 
00286             for( Range::iterator ielem = elems.begin(); ielem != elems.end(); ++ielem, ++index )
00287             {
00288                 // Projection
00289                 error = hirec.hiproj_walf_in_element( *ielem, nvpe, nsamples,
00290                                                       &( testnaturalcoords[nvpe * nsamples * index] ),
00291                                                       &( testpnts[3 * nsamples * index] ) );MB_CHK_ERR( error );
00292             }
00293 
00294             // Estimate error
00295             obj->compute_projecterror( 3, elems.size() * nsamples, &( testpnts[0] ), l1err, l2err, linferr );
00296             geoml1errs[ideg + 1].push_back( l1err );
00297             geoml2errs[ideg + 1].push_back( l2err );
00298             geomlinferrs[ideg + 1].push_back( linferr );
00299         }
00300     }
00301 
00302     std::cout << "Mesh Size: ";
00303 
00304     for( size_t i = 0; i < meshsizes.size(); ++i )
00305     {
00306         std::cout << meshsizes[i] << " ";
00307     }
00308 
00309     std::cout << std::endl;
00310     std::cout << "Degrees: 0 ";
00311 
00312     for( size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
00313     {
00314         std::cout << degs2fit[ideg] << " ";
00315     }
00316 
00317     std::cout << std::endl;
00318     std::cout << "L1-norm error: \n";
00319 
00320     for( size_t i = 0; i < geoml1errs.size(); ++i )
00321     {
00322         for( size_t j = 0; j < geoml1errs[i].size(); ++j )
00323         {
00324             std::cout << geoml1errs[i][j] << " ";
00325         }
00326 
00327         std::cout << std::endl;
00328     }
00329 
00330     std::cout << "L2-norm error: \n";
00331 
00332     for( size_t i = 0; i < geoml2errs.size(); ++i )
00333     {
00334         for( size_t j = 0; j < geoml2errs[i].size(); ++j )
00335         {
00336             std::cout << geoml2errs[i][j] << " ";
00337         }
00338 
00339         std::cout << std::endl;
00340     }
00341 
00342     std::cout << "Linf-norm error: \n";
00343 
00344     for( size_t i = 0; i < geomlinferrs.size(); ++i )
00345     {
00346         for( size_t j = 0; j < geomlinferrs[i].size(); ++j )
00347         {
00348             std::cout << geomlinferrs[i][j] << " ";
00349         }
00350 
00351         std::cout << std::endl;
00352     }
00353 
00354     return error;
00355 }
00356 
00357 void usage()
00358 {
00359     std::cout << "usage: ./uref_hirec_test <mesh file> -degree <degree> -interp <0=least square, "
00360                  "1=interpolation> -dim <mesh dimension> -geom <s=sphere,t=torus>"
00361               << std::endl;
00362     std::cout << "Example: ./uref_hirec_test $MOAB_SOURCE_DIR/MeshFiles/unittest/sphere_tris_5.vtk "
00363                  "-degree 2 -interp 0 -dim 2 -geom s"
00364               << std::endl;
00365 }
00366 
00367 int main( int argc, char* argv[] )
00368 {
00369 #ifdef MOAB_HAVE_MPI
00370     MPI_Init( &argc, &argv );
00371     int nprocs, rank;
00372     MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
00373     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00374 #endif
00375 
00376     std::string infile = TestDir + "unittest/sphere_tris_5.vtk";
00377 
00378     int degree = 2, dim = 2, geom = 0;
00379     bool interp = false;
00380     ErrorCode error;
00381 
00382     if( argc == 10 )
00383     {
00384         infile      = std::string( argv[1] );
00385         bool hasdim = false;
00386 
00387         for( int i = 2; i < argc; ++i )
00388         {
00389             if( i + 1 != argc )
00390             {
00391                 if( std::string( argv[i] ) == "-degree" )
00392                 {
00393                     degree = atoi( argv[++i] );
00394                 }
00395                 else if( std::string( argv[i] ) == "-interp" )
00396                 {
00397                     interp = atoi( argv[++i] );
00398                 }
00399                 else if( std::string( argv[i] ) == "-dim" )
00400                 {
00401                     dim    = atoi( argv[++i] );
00402                     hasdim = true;
00403                 }
00404                 else if( std::string( argv[i] ) == "-geom" )
00405                 {
00406                     geom = std::string( argv[++i] ) == "s" ? 0 : 1;
00407                 }
00408                 else
00409                 {
00410 #ifdef MOAB_HAVE_MPI
00411 
00412                     if( 0 == rank )
00413                     {
00414                         usage();
00415                     }
00416                     MPI_Finalize();
00417 
00418 #else
00419                     usage();
00420 #endif
00421                     return -1;
00422                 }
00423             }
00424         }
00425 
00426         if( !hasdim )
00427         {
00428 #ifdef MOAB_HAVE_MPI
00429 
00430             if( 0 == rank )
00431             {
00432                 std::cout << "Dimension of input mesh should be provided, positive and less than 3" << std::endl;
00433             }
00434 
00435 #else
00436             std::cout << "Dimension of input mesh should be provided, positive and less than 3" << std::endl;
00437 #endif
00438             return -1;
00439         }
00440 
00441         if( degree <= 0 || dim > 2 || dim <= 0 )
00442         {
00443 #ifdef MOAB_HAVE_MPI
00444 
00445             if( 0 == rank )
00446             {
00447                 std::cout << "Input degree should be positive number;\n";
00448                 std::cout << "Input dimesion should be positive and less than 3;" << std::endl;
00449             }
00450 
00451 #else
00452             std::cout << "Input degree should be positive number;\n";
00453             std::cout << "Input dimesion should be positive and less than 3;" << std::endl;
00454 #endif
00455             return -1;
00456         }
00457 
00458 #ifdef MOAB_HAVE_MPI
00459 
00460         if( 0 == rank )
00461         {
00462             std::cout << "Testing on " << infile << " with dimension " << dim << "\n";
00463             std::string opts = interp ? "interpolation" : "least square fitting";
00464             std::cout << "High order reconstruction with degree " << degree << " " << opts << std::endl;
00465         }
00466 
00467 #else
00468         std::cout << "Testing on " << infile << " with dimension " << dim << "\n";
00469         std::string opts = interp ? "interpolation" : "least square fitting";
00470         std::cout << "High order reconstruction with degree " << degree << " " << opts << std::endl;
00471 #endif
00472     }
00473     else
00474     {
00475         if( argc > 1 )
00476         {
00477             usage();
00478             return -1;
00479         }
00480     }
00481 
00482     {
00483         int level_degrees[3] = { 2, 2, 2 };
00484         int num_levels       = 3;
00485 
00486         // create the geometry object
00487         geomObject* obj;
00488         if( geom )
00489             obj = new torus();
00490         else
00491             obj = new sphere();
00492 
00493         error = test_closedsurface_mesh( infile.c_str(), level_degrees, num_levels, degree, interp, dim, obj );MB_CHK_ERR( error );
00494 
00495         std::vector< int > degs2fit( 6 );
00496         for( int d = 1; d <= 6; ++d )
00497         {
00498             degs2fit[d - 1] = d;
00499         }
00500 
00501         // Call the higher order reconstruction routines and compute error convergence
00502         error = closedsurface_uref_hirec_convergence_study( infile.c_str(), level_degrees, num_levels, degs2fit, interp,
00503                                                             obj );MB_CHK_ERR( error );
00504         // cleanup memory
00505         delete obj;
00506     }
00507 
00508 #ifdef MOAB_HAVE_MPI
00509     MPI_Finalize();
00510 #endif
00511     return 0;
00512 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines