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