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