MOAB: Mesh Oriented datABase  (version 5.3.0)
hireconst_test.cpp
Go to the documentation of this file.
00001 /*This unit test is for high order reconstruction capability, which could be used for mesh
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 <cmath>
00022 
00023 #ifdef MOAB_HAVE_MPI
00024 #include "moab/ParallelComm.hpp"
00025 #include "MBParallelConventions.h"
00026 #include "ReadParallel.hpp"
00027 #include "moab/FileOptions.hpp"
00028 #include "MBTagConventions.hpp"
00029 #include "moab_mpi.h"
00030 #endif
00031 
00032 using namespace moab;
00033 
00034 #ifdef MOAB_HAVE_MPI
00035 std::string read_options;
00036 #endif
00037 
00038 ErrorCode load_meshset_hirec( const char* infile, Interface* mbimpl, EntityHandle& meshset, ParallelComm*& pc,
00039                               const int degree = 0, const int dim = 2 );
00040 ErrorCode test_mesh( const char* infile, const int degree, const bool interp, const int dim );
00041 
00042 void compute_linear_coords( const int nvpe, double* elemcoords, double* naturals, double* linearcoords );
00043 
00044 ErrorCode create_unitsq_tris( Interface* mbImpl, size_t n, std::vector< EntityHandle >& tris );
00045 ErrorCode create_unitsq_quads( Interface* mbImpl, size_t n, std::vector< EntityHandle >& quads );
00046 ErrorCode test_unitsq_tris();
00047 ErrorCode test_unitsq_quads();
00048 ErrorCode test_unitsphere();
00049 ErrorCode test_unitcircle();
00050 
00051 ErrorCode exact_error_torus( const double R, const double r, const double center[3], int npnts, double* pnt,
00052                              double& error_l1, double& error_l2, double& error_li );
00053 ErrorCode project_exact_torus( Interface* mbImpl, EntityHandle meshset, int dim, const double R, const double r,
00054                                const double center[3] );
00055 
00056 int main( int argc, char* argv[] )
00057 {
00058 #ifdef MOAB_HAVE_MPI
00059     MPI_Init( &argc, &argv );
00060     int nprocs, rank;
00061     MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
00062     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00063 #endif
00064     std::string infile;
00065     int degree = 2, dim = 0;
00066     bool interp = false;
00067     ErrorCode error;
00068 
00069     if( argc == 1 )
00070     {
00071         error = test_unitsq_tris();MB_CHK_ERR( error );
00072         error = test_unitsq_quads();MB_CHK_ERR( error );
00073         error = test_unitsphere();MB_CHK_ERR( error );
00074         error = test_unitcircle();MB_CHK_ERR( error );
00075 #ifdef MOAB_HAVE_MPI
00076         MPI_Finalize();
00077 #endif
00078         return 0;
00079     }
00080     else
00081     {
00082         infile      = std::string( argv[1] );
00083         bool hasdim = false;
00084 
00085         for( int i = 2; i < argc; ++i )
00086         {
00087             if( i + 1 != argc )
00088             {
00089                 if( std::string( argv[i] ) == "-degree" ) { degree = atoi( argv[++i] ); }
00090                 else if( std::string( argv[i] ) == "-interp" )
00091                 {
00092                     interp = atoi( argv[++i] );
00093                 }
00094                 else if( std::string( argv[i] ) == "-dim" )
00095                 {
00096                     dim    = atoi( argv[++i] );
00097                     hasdim = true;
00098                 }
00099                 else
00100                 {
00101                     std::cout << argv[i] << std::endl;
00102                     std::cout << "usage: " << argv[0]
00103                               << " <mesh file> -degree <degree> -interp <0=least square, "
00104                                  "1=interpolation> -dim <mesh dimension>"
00105                               << std::endl;
00106                     return 0;
00107                 }
00108             }
00109         }
00110 
00111         if( !hasdim )
00112         {
00113             std::cout << "Dimension of input mesh should be provided, positive and less than 3" << std::endl;
00114             return 0;
00115         }
00116 
00117         if( degree <= 0 || dim > 2 || dim <= 0 )
00118         {
00119             std::cout << "Input degree should be positive number;\n";
00120             std::cout << "Input dimesion should be positive and less than 3;" << std::endl;
00121             return -1;
00122         }
00123 
00124         std::cout << "Testing on " << infile << " with dimension " << dim << "\n";
00125         std::string opts = interp ? "interpolation" : "least square fitting";
00126         std::cout << "High order reconstruction with degree " << degree << " " << opts << std::endl;
00127     }
00128 
00129     error = test_mesh( infile.c_str(), degree, interp, dim );MB_CHK_ERR( error );
00130 #ifdef MOAB_HAVE_MPI
00131     MPI_Finalize();
00132 #endif
00133     return 0;
00134 }
00135 
00136 ErrorCode load_meshset_hirec( const char* infile, Interface* mbimpl, EntityHandle& meshset, ParallelComm*& pc,
00137                               const int degree, const int dim )
00138 {
00139     ErrorCode error;
00140     error = mbimpl->create_meshset( moab::MESHSET_SET, meshset );MB_CHK_ERR( error );
00141 #ifdef MOAB_HAVE_MPI
00142     int nprocs, rank;
00143     MPI_Comm comm = MPI_COMM_WORLD;
00144     MPI_Comm_size( comm, &nprocs );
00145     MPI_Comm_rank( comm, &rank );
00146     EntityHandle partnset;
00147     error = mbimpl->create_meshset( moab::MESHSET_SET, partnset );MB_CHK_ERR( error );
00148 
00149     if( nprocs > 1 ) { pc = moab::ParallelComm::get_pcomm( mbimpl, partnset, &comm ); }
00150 
00151     if( nprocs > 1 )
00152     {
00153         int nghlayers = degree > 0 ? HiReconstruction::estimate_num_ghost_layers( degree, true ) : 0;
00154 
00155         if( nghlayers )
00156         {
00157             // get ghost layers
00158             if( dim == 2 )
00159             {
00160                 read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
00161                                "SHARED_ENTS;PARALLEL_GHOST=2.0.";
00162             }
00163             else if( dim == 1 )
00164             {
00165                 read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
00166                                "SHARED_ENTS;PARALLEL_GHOST=1.0.";
00167             }
00168             else
00169             {
00170                 MB_SET_ERR( MB_FAILURE, "unsupported dimension" );
00171             }
00172 
00173             read_options += ( '0' + nghlayers );
00174         }
00175         else
00176         {
00177             read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
00178         }
00179 
00180         error = mbimpl->load_file( infile, &meshset, read_options.c_str() );MB_CHK_ERR( error );
00181     }
00182     else
00183     {
00184         error = mbimpl->load_file( infile, &meshset );MB_CHK_ERR( error );
00185     }
00186 
00187 #else
00188     assert( !pc && degree && dim );
00189     error = mbimpl->load_file( infile, &meshset );MB_CHK_ERR( error );
00190 #endif
00191     return error;
00192 }
00193 
00194 ErrorCode test_mesh( const char* infile, const int degree, const bool interp, const int dim )
00195 {
00196     Core moab;
00197     Interface* mbimpl = &moab;
00198     ParallelComm* pc  = NULL;
00199     EntityHandle meshset;
00200     // load mesh file
00201     ErrorCode error;
00202     error = load_meshset_hirec( infile, mbimpl, meshset, pc, degree, dim );MB_CHK_ERR( error );
00203     // project to exact surface: torus
00204     double center[3] = { 0, 0, 0 };
00205     double R = 1, r = 0.3;
00206     error = project_exact_torus( mbimpl, meshset, dim, R, r, center );CHECK_ERR( error );
00207     // initialize
00208     HiReconstruction hirec( dynamic_cast< Core* >( mbimpl ), pc, meshset );
00209     Range elems;
00210     error = mbimpl->get_entities_by_dimension( meshset, dim, elems );MB_CHK_ERR( error );
00211 
00212     // reconstruction
00213     if( dim == 2 )
00214     {
00215         // error = hirec.reconstruct3D_surf_geom(degree, interp, false); MB_CHK_ERR(error);
00216         error = hirec.reconstruct3D_surf_geom( degree, interp, true );MB_CHK_ERR( error );
00217     }
00218     else if( dim == 1 )
00219     {
00220         error = hirec.reconstruct3D_curve_geom( degree, interp, true );MB_CHK_ERR( error );
00221     }
00222 
00223     // fitting
00224     double mxdist = 0, errl1 = 0, errl2 = 0, errli = 0;
00225     double* pnts_proj = new double[elems.size() * 3];
00226 
00227     for( Range::iterator ielem = elems.begin(); ielem != elems.end(); ++ielem )
00228     {
00229         int nvpe;
00230         const EntityHandle* conn;
00231         error = mbimpl->get_connectivity( *ielem, conn, nvpe );MB_CHK_ERR( error );
00232         double w = 1.0 / (double)nvpe;
00233         std::vector< double > naturalcoords2fit( nvpe, w );
00234         double newcoords[3], linearcoords[3];
00235         error = hirec.hiproj_walf_in_element( *ielem, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
00236         pnts_proj[3 * ( *ielem - *elems.begin() )]     = newcoords[0];
00237         pnts_proj[3 * ( *ielem - *elems.begin() ) + 1] = newcoords[1];
00238         pnts_proj[3 * ( *ielem - *elems.begin() ) + 2] = newcoords[2];
00239         std::vector< double > coords( 3 * nvpe );
00240         error = mbimpl->get_coords( conn, nvpe, &( coords[0] ) );MB_CHK_ERR( error );
00241         compute_linear_coords( nvpe, &( coords[0] ), &( naturalcoords2fit[0] ), linearcoords );
00242         mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
00243     }
00244 
00245     delete[] pnts_proj;
00246     // compute error for torus
00247     error = exact_error_torus( R, r, center, (int)elems.size(), pnts_proj, errl1, errl2, errli );MB_CHK_ERR( error );
00248     std::cout << "Errors using exact torus for degree " << degree << " fit : L1 = " << errl1 << ", L2 = " << errl2
00249               << ", Linf = " << errli << std::endl;
00250     std::cout << "Maximum projection lift is " << mxdist << std::endl;
00251     return error;
00252 }
00253 
00254 ErrorCode create_unitsq_tris( Interface* mbImpl, size_t n, std::vector< EntityHandle >& tris )
00255 {
00256     if( n < 2 ) { MB_SET_ERR( MB_FAILURE, "n must be at least 2" ); }
00257 
00258     ErrorCode error;
00259     std::vector< EntityHandle > verts( n * n );
00260     size_t istr = tris.size();
00261     tris.resize( istr + 2 * ( n - 1 ) * ( n - 1 ) );
00262     double istep = 1.0 / (double)( n - 1 );
00263 
00264     for( size_t j = 0; j < n; ++j )
00265     {
00266         for( size_t i = 0; i < n; ++i )
00267         {
00268             double coord[3] = { i * istep, j * istep, 0 };
00269             EntityHandle temp;
00270             error = mbImpl->create_vertex( coord, temp );MB_CHK_ERR( error );
00271             verts[j * n + i] = temp;
00272         }
00273     }
00274 
00275     for( size_t jj = 0; jj < n - 1; ++jj )
00276     {
00277         for( size_t ii = 0; ii < n - 1; ++ii )
00278         {
00279             EntityHandle conn[3] = { verts[jj * n + ii], verts[( jj + 1 ) * n + ii + 1], verts[( jj + 1 ) * n + ii] };
00280             error                = mbImpl->create_element( MBTRI, conn, 3, tris[istr + 2 * jj * ( n - 1 ) + 2 * ii] );MB_CHK_ERR( error );
00281             conn[0] = verts[jj * n + ii];
00282             conn[1] = verts[jj * n + ii + 1];
00283             conn[2] = verts[( jj + 1 ) * n + ii + 1];
00284             error   = mbImpl->create_element( MBTRI, conn, 3, tris[istr + 2 * jj * ( n - 1 ) + 2 * ii + 1] );MB_CHK_ERR( error );
00285         }
00286     }
00287 
00288     return error;
00289 }
00290 
00291 ErrorCode create_unitsq_quads( Interface* mbImpl, size_t n, std::vector< EntityHandle >& quads )
00292 {
00293     if( n < 2 ) { MB_SET_ERR( MB_FAILURE, "n must be at least 2" ); }
00294 
00295     ErrorCode error;
00296     std::vector< EntityHandle > verts( n * n );
00297     size_t istr = quads.size();
00298     quads.resize( istr + ( n - 1 ) * ( n - 1 ) );
00299     double istep = 1.0 / (double)( n - 1 );
00300 
00301     for( size_t j = 0; j < n; ++j )
00302     {
00303         for( size_t i = 0; i < n; ++i )
00304         {
00305             double coord[3] = { i * istep, j * istep, 0 };
00306             error           = mbImpl->create_vertex( coord, verts[j * n + i] );MB_CHK_ERR( error );
00307         }
00308     }
00309 
00310     for( size_t jj = 0; jj < n - 1; ++jj )
00311     {
00312         for( size_t ii = 0; ii < n - 1; ++ii )
00313         {
00314             EntityHandle conn[4] = { verts[jj * n + ii], verts[jj * n + ii + 1], verts[( jj + 1 ) * n + ii + 1],
00315                                      verts[( jj + 1 ) * n + ii] };
00316             error                = mbImpl->create_element( MBQUAD, conn, 4, quads[istr + jj * ( n - 1 ) + ii] );MB_CHK_ERR( error );
00317         }
00318     }
00319 
00320     return error;
00321 }
00322 
00323 ErrorCode test_unitsq_tris()
00324 {
00325     ErrorCode error;
00326 
00327     for( size_t n = 2; n <= 2; ++n )
00328     {
00329         Core moab;
00330         Interface* mbImpl = &moab;
00331         std::vector< EntityHandle > tris;
00332         error = create_unitsq_tris( mbImpl, n, tris );MB_CHK_ERR( error );
00333         EntityHandle meshIn = 0;
00334         HiReconstruction hirec( dynamic_cast< Core* >( mbImpl ), 0, meshIn );
00335 
00336         for( int degree = 2; degree <= 2; ++degree )
00337         {
00338             // reconstruct geometry, interpolation
00339             hirec.reconstruct3D_surf_geom( degree, true, true, true );
00340             // test fitting result
00341             double mxdist = 0;
00342 
00343             for( size_t itri = 0; itri < tris.size(); ++itri )
00344             {
00345                 const int nvpe                 = 3;
00346                 double naturalcoords2fit[nvpe] = { 1.0 / (double)nvpe, 1.0 / (double)nvpe, 1.0 / (double)nvpe },
00347                        newcoords[3];
00348                 error = hirec.hiproj_walf_in_element( tris[itri], nvpe, 1, naturalcoords2fit, newcoords );MB_CHK_ERR( error );
00349                 std::vector< EntityHandle > conn;
00350                 error = mbImpl->get_connectivity( &( tris[itri] ), 1, conn );MB_CHK_ERR( error );
00351                 double coords[3 * nvpe], linearcoords[3];
00352                 error = mbImpl->get_coords( &( conn[0] ), nvpe, coords );MB_CHK_ERR( error );
00353                 compute_linear_coords( nvpe, coords, naturalcoords2fit, linearcoords );
00354                 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
00355             }
00356 
00357             std::cout << "triangulated unit square n= " << n << " degree= " << degree << " interpolation:\n";
00358             std::cout << "maximum projection lift is " << mxdist << std::endl;
00359             // for debug
00360             // return error;
00361             mxdist = 0;
00362             // reconstruct geometry, least square fitting
00363             hirec.reconstruct3D_surf_geom( degree, false, false, true );
00364 
00365             // test fitting result
00366             for( size_t itri = 0; itri < tris.size(); ++itri )
00367             {
00368                 const int nvpe                 = 3;
00369                 double naturalcoords2fit[nvpe] = { 1.0 / (double)nvpe, 1.0 / (double)nvpe, 1.0 / (double)nvpe },
00370                        newcoords[3];
00371                 error = hirec.hiproj_walf_in_element( tris[itri], nvpe, 1, naturalcoords2fit, newcoords );MB_CHK_ERR( error );
00372                 std::vector< EntityHandle > conn;
00373                 error = mbImpl->get_connectivity( &( tris[itri] ), 1, conn );MB_CHK_ERR( error );
00374                 double coords[3 * nvpe], linearcoords[3];
00375                 error = mbImpl->get_coords( &( conn[0] ), nvpe, coords );MB_CHK_ERR( error );
00376                 compute_linear_coords( nvpe, coords, naturalcoords2fit, linearcoords );
00377                 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
00378             }
00379 
00380             std::cout << "unit square n= " << n << " degree= " << degree << " least square:\n";
00381             std::cout << "maximum projection lift is " << mxdist << std::endl;
00382         }
00383     }
00384 
00385     return error;
00386 }
00387 
00388 void compute_linear_coords( const int nvpe, double* elemcoords, double* naturals, double* linearcoords )
00389 {
00390     assert( elemcoords && linearcoords );
00391 
00392     for( int i = 0; i < 3; ++i )
00393     {
00394         linearcoords[i] = 0;
00395 
00396         for( int j = 0; j < nvpe; ++j )
00397         {
00398             linearcoords[i] += naturals[j] * elemcoords[3 * j + i];
00399         }
00400     }
00401 }
00402 
00403 ErrorCode test_unitsq_quads()
00404 {
00405     ErrorCode error;
00406 
00407     for( size_t n = 2; n <= 8; ++n )
00408     {
00409         Core moab;
00410         Interface* mbImpl = &moab;
00411         std::vector< EntityHandle > quads;
00412         error = create_unitsq_quads( mbImpl, n, quads );MB_CHK_ERR( error );
00413         EntityHandle meshIn = 0;
00414         HiReconstruction hirec( dynamic_cast< Core* >( mbImpl ), 0, meshIn );
00415 
00416         for( int degree = 1; degree <= 6; ++degree )
00417         {
00418             // reconstruct geometry, interpolation
00419             hirec.reconstruct3D_surf_geom( degree, true, false, true );
00420             // test fitting result
00421             double mxdist = 0;
00422 
00423             for( size_t iquad = 0; iquad < quads.size(); ++iquad )
00424             {
00425                 const int nvpe                 = 4;
00426                 double w                       = 1.0 / (double)nvpe;
00427                 double naturalcoords2fit[nvpe] = { w, w, w, w }, newcoords[3];
00428                 error = hirec.hiproj_walf_in_element( quads[iquad], nvpe, 1, naturalcoords2fit, newcoords );MB_CHK_ERR( error );
00429                 std::vector< EntityHandle > conn;
00430                 error = mbImpl->get_connectivity( &( quads[iquad] ), 1, conn );MB_CHK_ERR( error );
00431                 double coords[3 * nvpe], linearcoords[3];
00432                 error = mbImpl->get_coords( &( conn[0] ), nvpe, coords );MB_CHK_ERR( error );
00433                 compute_linear_coords( nvpe, coords, naturalcoords2fit, linearcoords );
00434                 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
00435             }
00436 
00437             std::cout << "quadrilateral unit square n= " << n << " degree= " << degree << " interpolation:\n";
00438             std::cout << "maximum projection lift is " << mxdist << std::endl;
00439             mxdist = 0;
00440             // reconstruct geometry, least square fitting
00441             hirec.reconstruct3D_surf_geom( degree, false, false, true );
00442 
00443             // test fitting result
00444             for( size_t iquad = 0; iquad < quads.size(); ++iquad )
00445             {
00446                 const int nvpe                 = 4;
00447                 double w                       = 1.0 / (double)nvpe;
00448                 double naturalcoords2fit[nvpe] = { w, w, w, w }, newcoords[3];
00449                 error = hirec.hiproj_walf_in_element( quads[iquad], nvpe, 1, naturalcoords2fit, newcoords );MB_CHK_ERR( error );
00450                 std::vector< EntityHandle > conn;
00451                 error = mbImpl->get_connectivity( &( quads[iquad] ), 1, conn );MB_CHK_ERR( error );
00452                 double coords[3 * nvpe], linearcoords[3];
00453                 error = mbImpl->get_coords( &( conn[0] ), nvpe, coords );MB_CHK_ERR( error );
00454                 compute_linear_coords( nvpe, coords, naturalcoords2fit, linearcoords );
00455                 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
00456             }
00457 
00458             std::cout << "quadrilateral unit square n= " << n << " degree= " << degree << " least square:\n";
00459             std::cout << "maximum projection lift is " << mxdist << std::endl;
00460         }
00461     }
00462 
00463     return error;
00464 }
00465 
00466 ErrorCode test_unitsphere()
00467 {
00468     // path to test files
00469     int nfiles              = 4;
00470     std::string filenames[] = { TestDir + "/sphere_tris_5.vtk", TestDir + "/sphere_tris_20.vtk",
00471                                 TestDir + "/sphere_quads_5.vtk", TestDir + "/sphere_quads_20.vtk" };
00472     ErrorCode error;
00473     int maxdeg = 6;
00474 
00475     for( int ifile = 0; ifile < nfiles; ++ifile )
00476     {
00477         Core moab;
00478         Interface* mbimpl = &moab;
00479         ParallelComm* pc  = NULL;
00480         EntityHandle meshset;
00481         // load file
00482         error = load_meshset_hirec( filenames[ifile].c_str(), mbimpl, meshset, pc, maxdeg );MB_CHK_ERR( error );
00483         // initialize
00484         HiReconstruction hirec( &moab, pc, meshset );
00485         Range elems;
00486         error = mbimpl->get_entities_by_dimension( meshset, 2, elems );MB_CHK_ERR( error );
00487 
00488         // reconstruction
00489         for( int degree = 1; degree <= maxdeg; ++degree )
00490         {
00491             hirec.reconstruct3D_surf_geom( degree, true, false, true );
00492             // fitting
00493             double mxdist = 0, mxerr = 0;
00494 
00495             for( Range::iterator ielem = elems.begin(); ielem != elems.end(); ++ielem )
00496             {
00497                 int nvpe;
00498                 const EntityHandle* conn;
00499                 error = mbimpl->get_connectivity( *ielem, conn, nvpe );MB_CHK_ERR( error );
00500                 double w = 1.0 / (double)nvpe;
00501                 std::vector< double > naturalcoords2fit( nvpe, w );
00502                 double newcoords[3], linearcoords[3];
00503                 error = hirec.hiproj_walf_in_element( *ielem, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
00504                 std::vector< double > coords( 3 * nvpe );
00505                 error = mbimpl->get_coords( conn, nvpe, &( coords[0] ) );MB_CHK_ERR( error );
00506                 compute_linear_coords( nvpe, &( coords[0] ), &( naturalcoords2fit[0] ), linearcoords );
00507                 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
00508                 mxerr  = std::max( mxerr, fabs( DGMSolver::vec_2norm( 3, newcoords ) - 1 ) );
00509             }
00510 
00511             std::cout << filenames[ifile] << ": unit sphere"
00512                       << " degree= " << degree << " interpolation:\n";
00513             std::cout << "maximum projection lift is " << mxdist << ", maximum error is " << mxerr << std::endl;
00514             mxdist = 0;
00515             mxerr  = 0;
00516             hirec.reconstruct3D_surf_geom( degree, false, false, true );
00517 
00518             // fitting
00519             for( Range::iterator ielem = elems.begin(); ielem != elems.end(); ++ielem )
00520             {
00521                 int nvpe;
00522                 const EntityHandle* conn;
00523                 error = mbimpl->get_connectivity( *ielem, conn, nvpe );MB_CHK_ERR( error );
00524                 double w = 1.0 / (double)nvpe;
00525                 std::vector< double > naturalcoords2fit( nvpe, w );
00526                 double newcoords[3], linearcoords[3];
00527                 error = hirec.hiproj_walf_in_element( *ielem, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
00528                 std::vector< double > coords( 3 * nvpe );
00529                 error = mbimpl->get_coords( conn, nvpe, &( coords[0] ) );MB_CHK_ERR( error );
00530                 compute_linear_coords( nvpe, &( coords[0] ), &( naturalcoords2fit[0] ), linearcoords );
00531                 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
00532                 mxerr  = std::max( mxerr, fabs( DGMSolver::vec_2norm( 3, newcoords ) - 1 ) );
00533             }
00534 
00535             std::cout << filenames[ifile] << ": unit sphere"
00536                       << " degree= " << degree << " least square:\n";
00537             std::cout << "maximum projection lift is " << mxdist << ", maximum error is " << mxerr << std::endl;
00538         }
00539     }
00540 
00541     return error;
00542 }
00543 
00544 ErrorCode test_unitcircle()
00545 {
00546     // path to test files
00547     int nfiles              = 4;
00548     std::string filenames[] = { TestDir + "/circle_3.vtk", TestDir + "/circle_4.vtk", TestDir + "/circle_10.vtk",
00549                                 TestDir + "/circle_20.vtk" };
00550     ErrorCode error;
00551     int maxdeg = 6;
00552 
00553     for( int ifile = 0; ifile < nfiles; ++ifile )
00554     {
00555         Core moab;
00556         Interface* mbimpl = &moab;
00557         ParallelComm* pc  = 0;
00558         EntityHandle meshset;
00559         int dim = 1;
00560         // load file
00561         error = load_meshset_hirec( filenames[ifile].c_str(), mbimpl, meshset, pc, maxdeg, dim );MB_CHK_ERR( error );
00562         // initialize
00563         HiReconstruction hirec( &moab, pc, meshset );
00564         Range edges;
00565         error = mbimpl->get_entities_by_dimension( meshset, dim, edges );MB_CHK_ERR( error );
00566 
00567         // reconstruction
00568         for( int degree = 1; degree <= maxdeg; ++degree )
00569         {
00570             hirec.reconstruct3D_curve_geom( degree, true, false, true );
00571             // fitting
00572             double mxdist = 0, mxerr = 0;
00573 
00574             for( Range::iterator iedge = edges.begin(); iedge != edges.end(); ++iedge )
00575             {
00576                 int nvpe;
00577                 const EntityHandle* conn;
00578                 error = mbimpl->get_connectivity( *iedge, conn, nvpe );MB_CHK_ERR( error );
00579                 double w = 1.0 / (double)nvpe;
00580                 std::vector< double > naturalcoords2fit( nvpe, w );
00581                 double newcoords[3], linearcoords[3];
00582                 error = hirec.hiproj_walf_in_element( *iedge, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
00583                 std::vector< double > coords( 3 * nvpe );
00584                 error = mbimpl->get_coords( conn, nvpe, &( coords[0] ) );MB_CHK_ERR( error );
00585                 compute_linear_coords( nvpe, &( coords[0] ), &( naturalcoords2fit[0] ), linearcoords );
00586                 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
00587                 mxerr  = std::max( mxerr, fabs( DGMSolver::vec_2norm( 3, newcoords ) - 1 ) );
00588             }
00589 
00590             std::cout << filenames[ifile] << ": unit circle"
00591                       << " degree= " << degree << " interpolation:\n";
00592             std::cout << "maximum projection lift is " << mxdist << ", maximum error is " << mxerr << std::endl;
00593             mxdist = 0;
00594             mxerr  = 0;
00595             hirec.reconstruct3D_curve_geom( degree, false, false, true );
00596 
00597             // fitting
00598             for( Range::iterator iedge = edges.begin(); iedge != edges.end(); ++iedge )
00599             {
00600                 int nvpe;
00601                 const EntityHandle* conn;
00602                 error = mbimpl->get_connectivity( *iedge, conn, nvpe );MB_CHK_ERR( error );
00603                 double w = 1.0 / (double)nvpe;
00604                 std::vector< double > naturalcoords2fit( nvpe, w );
00605                 double newcoords[3], linearcoords[3];
00606                 error = hirec.hiproj_walf_in_element( *iedge, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
00607                 std::vector< double > coords( 3 * nvpe );
00608                 error = mbimpl->get_coords( conn, nvpe, &( coords[0] ) );MB_CHK_ERR( error );
00609                 compute_linear_coords( nvpe, &( coords[0] ), &( naturalcoords2fit[0] ), linearcoords );
00610                 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
00611                 mxerr  = std::max( mxerr, fabs( DGMSolver::vec_2norm( 3, newcoords ) - 1 ) );
00612             }
00613 
00614             std::cout << filenames[ifile] << ": unit circle"
00615                       << " degree= " << degree << " least square:\n";
00616             std::cout << "maximum projection lift is " << mxdist << ", maximum error is " << mxerr << std::endl;
00617         }
00618     }
00619 
00620     return error;
00621 }
00622 
00623 ErrorCode project_exact_torus( Interface* mbImpl, EntityHandle meshset, int dim, const double R, const double r,
00624                                const double center[] )
00625 {
00626     ErrorCode error;
00627     Range elems, verts;
00628     error = mbImpl->get_entities_by_dimension( meshset, dim, elems );CHECK_ERR( error );
00629     error = mbImpl->get_connectivity( elems, verts );CHECK_ERR( error );
00630     double pnts[3] = { 0, 0, 0 }, cnt[3] = { 0, 0, 0 }, nrms[3] = { 0, 0, 0 };
00631     double x = 0, y = 0, z = 0, d1 = 0, d2 = 0;
00632 
00633     for( int i = 0; i < (int)verts.size(); i++ )
00634     {
00635         EntityHandle v = verts[i];
00636         error          = mbImpl->get_coords( &v, 1, &pnts[0] );CHECK_ERR( error );
00637         x       = pnts[0] - center[0];
00638         y       = pnts[1] - center[0];
00639         z       = pnts[2] - center[2];
00640         d1      = sqrt( x * x + y * y );
00641         cnt[0]  = R * x / d1;
00642         cnt[1]  = R * y / d1;
00643         cnt[2]  = 0;
00644         d2      = sqrt( ( x - cnt[0] ) * ( x - cnt[0] ) + ( y - cnt[1] ) * ( y - cnt[1] ) + z * z );
00645         nrms[0] = ( x - cnt[0] ) / d2;
00646         nrms[1] = ( y - cnt[1] ) / d2;
00647         nrms[2] = z / d2;
00648         pnts[0] = cnt[0] + r * nrms[0];
00649         pnts[1] = cnt[1] + r * nrms[1];
00650         pnts[2] = cnt[2] + r * nrms[2];
00651         error   = mbImpl->set_coords( &v, 1, &pnts[0] );CHECK_ERR( error );
00652     }
00653 
00654     return MB_SUCCESS;
00655 }
00656 
00657 ErrorCode exact_error_torus( const double R, const double r, const double center[], int npnts, double* pnts,
00658                              double& error_l1, double& error_l2, double& error_li )
00659 {
00660     error_l1 = 0;
00661     error_l2 = 0;
00662     error_li = 0;
00663     double x = 0, y = 0, z = 0;
00664     double error_single = 0;
00665 
00666     for( int i = 0; i < npnts; i++ )
00667     {
00668         x            = pnts[3 * i] - center[0];
00669         y            = pnts[3 * i + 1] - center[1];
00670         z            = pnts[3 * i + 2] - center[2];
00671         error_single = fabs( sqrt( pow( sqrt( pow( x, 2 ) + pow( y, 2 ) ) - R, 2 ) + pow( z, 2 ) ) - r );
00672         error_l1     = error_l1 + error_single;
00673         error_l2     = error_l2 + error_single * error_single;
00674 
00675         if( error_li < error_single ) { error_li = error_single; }
00676     }
00677 
00678     error_l1 = error_l1 / (double)npnts;
00679     error_l2 = sqrt( error_l2 / (double)npnts );
00680     return MB_SUCCESS;
00681 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines