MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 }