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