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