MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 #include "Mesquite.hpp" 00002 #include "MsqIBase.hpp" 00003 #include "MsqIGeom.hpp" 00004 #include "MsqIMesh.hpp" 00005 #include "MBiMesh.hpp" 00006 #include "MeshImpl.hpp" 00007 #include "moab/Core.hpp" 00008 #include "moab/Skinner.hpp" 00009 #include "moab/LloydSmoother.hpp" 00010 #include "FacetModifyEngine.hpp" 00011 #include "MsqError.hpp" 00012 #include "InstructionQueue.hpp" 00013 #include "PatchData.hpp" 00014 #include "TerminationCriterion.hpp" 00015 #include "QualityAssessor.hpp" 00016 #include "SphericalDomain.hpp" 00017 #include "PlanarDomain.hpp" 00018 #include "MeshWriter.hpp" 00019 00020 #include "moab/Core.hpp" 00021 #include "moab/ProgOptions.hpp" 00022 #ifdef MOAB_HAVE_MPI 00023 #include "moab/ParallelComm.hpp" 00024 #endif 00025 00026 // algorithms 00027 #include "IdealWeightInverseMeanRatio.hpp" 00028 #include "TMPQualityMetric.hpp" 00029 #include "AspectRatioGammaQualityMetric.hpp" 00030 #include "ConditionNumberQualityMetric.hpp" 00031 #include "VertexConditionNumberQualityMetric.hpp" 00032 #include "MultiplyQualityMetric.hpp" 00033 #include "LPtoPTemplate.hpp" 00034 #include "LInfTemplate.hpp" 00035 #include "PMeanPTemplate.hpp" 00036 #include "SteepestDescent.hpp" 00037 #include "FeasibleNewton.hpp" 00038 #include "QuasiNewton.hpp" 00039 #include "ConjugateGradient.hpp" 00040 #include "SmartLaplacianSmoother.hpp" 00041 #include "Randomize.hpp" 00042 00043 #include "ReferenceMesh.hpp" 00044 #include "RefMeshTargetCalculator.hpp" 00045 #include "TShapeB1.hpp" 00046 #include "TQualityMetric.hpp" 00047 #include "IdealShapeTarget.hpp" 00048 00049 #include <sys/stat.h> 00050 #include <iostream> 00051 using std::cerr; 00052 using std::cout; 00053 using std::endl; 00054 00055 #include "iBase.h" 00056 00057 using namespace MBMesquite; 00058 00059 static int print_iGeom_error( const char* desc, int err, iGeom_Instance geom, const char* file, int line ) 00060 { 00061 char buffer[1024]; 00062 iGeom_getDescription( geom, buffer, sizeof( buffer ) ); 00063 buffer[sizeof( buffer ) - 1] = '\0'; 00064 00065 std::cerr << "ERROR: " << desc << std::endl 00066 << " Error code: " << err << std::endl 00067 << " Error desc: " << buffer << std::endl 00068 << " At : " << file << ':' << line << std::endl; 00069 00070 return -1; // must always return false or CHECK macro will break 00071 } 00072 00073 static int print_iMesh_error( const char* desc, int err, iMesh_Instance mesh, const char* file, int line ) 00074 { 00075 char buffer[1024]; 00076 iMesh_getDescription( mesh, buffer, sizeof( buffer ) ); 00077 buffer[sizeof( buffer ) - 1] = '\0'; 00078 00079 std::cerr << "ERROR: " << desc << std::endl 00080 << " Error code: " << err << std::endl 00081 << " Error desc: " << buffer << std::endl 00082 << " At : " << file << ':' << line << std::endl; 00083 00084 return -1; // must always return false or CHECK macro will break 00085 } 00086 00087 #define CHECK_IGEOM( STR ) \ 00088 if( err != iBase_SUCCESS ) return print_iGeom_error( STR, err, geom, __FILE__, __LINE__ ) 00089 00090 #define CHECK_IMESH( STR ) \ 00091 if( err != iBase_SUCCESS ) return print_iMesh_error( STR, err, instance, __FILE__, __LINE__ ) 00092 00093 const std::string default_file_name = std::string( MESH_DIR ) + std::string( "/surfrandomtris-4part.h5m" ); 00094 const std::string default_geometry_file_name = std::string( MESH_DIR ) + std::string( "/surfrandom.facet" ); 00095 00096 std::vector< double > solution_indicator; 00097 00098 int write_vtk_mesh( Mesh* mesh, const char* filename ); 00099 00100 // Construct a MeshTSTT from the file 00101 int get_imesh_mesh( MBMesquite::Mesh**, const char* file_name, int dimension ); 00102 00103 // Construct a MeshImpl from the file 00104 int get_native_mesh( MBMesquite::Mesh**, const char* file_name, int dimension ); 00105 00106 int get_itaps_domain( MeshDomain**, const char* ); 00107 int get_mesquite_domain( MeshDomain**, const char* ); 00108 00109 // Run FeasibleNewton solver 00110 int run_global_smoother( MeshDomainAssoc& mesh, MsqError& err, double OF_value = 1e-4 ); 00111 00112 // Run SmoothLaplacian solver 00113 int run_local_smoother( MeshDomainAssoc& mesh, MsqError& err, double OF_value = 1e-3 ); 00114 int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value = 1e-3 ); 00115 00116 int run_quality_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err ); 00117 00118 int run_solution_mesh_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err ); 00119 00120 bool file_exists( const std::string& name ) 00121 { 00122 struct stat buffer; 00123 return ( stat( name.c_str(), &buffer ) == 0 ); 00124 } 00125 00126 int main( int argc, char* argv[] ) 00127 { 00128 MBMesquite::MsqPrintError err( cout ); 00129 // command line arguments 00130 std::string file_name, geometry_file_name; 00131 bool use_native = false; 00132 int dimension = 2; 00133 00134 #ifdef MOAB_HAVE_MPI 00135 // MPI_Init(&argc, &argv); 00136 #endif 00137 ProgOptions opts; 00138 00139 opts.addOpt< void >( std::string( "native,N" ), std::string( "Use native representation (default=false)" ), 00140 &use_native ); 00141 opts.addOpt< int >( std::string( "dim,d" ), std::string( "Topological dimension of the mesh (default=2)" ), 00142 &dimension ); 00143 opts.addOpt< std::string >( std::string( "input_geo,i" ), 00144 std::string( "Input file name for the geometry (pattern=*.stl, *.facet)" ), 00145 &geometry_file_name ); 00146 opts.addOpt< std::string >( std::string( "input_mesh,m" ), 00147 std::string( "Input file name for the mesh (pattern=*.vtk, *.h5m)" ), &file_name ); 00148 00149 opts.parseCommandLine( argc, argv ); 00150 00151 if( !geometry_file_name.length() ) 00152 { 00153 file_name = default_file_name; 00154 geometry_file_name = default_geometry_file_name; 00155 cout << "No file specified: Using defaults.\n"; 00156 } 00157 cout << "\t Mesh filename = " << file_name << endl; 00158 cout << "\t Geometry filename = " << geometry_file_name << endl; 00159 00160 // Vector3D pnt(0,0,0); 00161 // Vector3D s_norm(0,0,1); 00162 // PlanarDomain plane(s_norm, pnt); 00163 00164 // PlanarDomain plane( PlanarDomain::XY ); 00165 MeshDomain* plane; 00166 int ierr; 00167 if( !file_exists( geometry_file_name ) ) geometry_file_name = ""; 00168 ierr = get_itaps_domain( &plane, geometry_file_name.c_str() ); // MB_CHK_ERR(ierr); 00169 00170 // Try running a global smoother on the mesh 00171 Mesh* mesh = 0; 00172 if( use_native ) 00173 { 00174 ierr = get_native_mesh( &mesh, file_name.c_str(), dimension ); // MB_CHK_ERR(ierr); 00175 } 00176 else 00177 { 00178 ierr = get_imesh_mesh( &mesh, file_name.c_str(), dimension ); // MB_CHK_ERR(ierr); 00179 } 00180 00181 if( !mesh ) 00182 { 00183 std::cerr << "Failed to load input file. Aborting." << std::endl; 00184 return 1; 00185 } 00186 00187 MeshDomainAssoc mesh_and_domain( mesh, plane ); 00188 00189 // ierr = write_vtk_mesh( mesh, "original.vtk");MB_CHK_ERR(ierr); 00190 // cout << "Wrote \"original.vtk\"" << endl; 00191 00192 // run_global_smoother( mesh_and_domain, err ); 00193 // if (err) return 1; 00194 00195 // Try running a local smoother on the mesh 00196 // Mesh* meshl=0; 00197 // if(use_native) 00198 // ierr = get_native_mesh(&meshl, file_name.c_str(), dimension); 00199 // else 00200 // ierr = get_imesh_mesh(&meshl, file_name.c_str(), dimension); 00201 // if (!mesh || ierr) { 00202 // std::cerr << "Failed to load input file. Aborting." << std::endl; 00203 // return 1; 00204 // } 00205 00206 // MeshDomainAssoc mesh_and_domain_local(meshl, plane); 00207 00208 // run_solution_mesh_optimizer( mesh_and_domain, err ); 00209 // if (err) return 1; 00210 00211 run_local_smoother( mesh_and_domain, err, 1e-4 ); // MB_CHK_ERR(err); 00212 if( err ) return 1; 00213 00214 double reps = 5e-2; 00215 for( int iter = 0; iter < 5; iter++ ) 00216 { 00217 00218 if( !( iter % 2 ) ) 00219 { 00220 run_local_smoother2( mesh_and_domain, err, 00221 reps * 10 ); // CHECK_IMESH("local smoother2 failed"); 00222 if( err ) return 1; 00223 } 00224 00225 // run_global_smoother( mesh_and_domain, err, reps );MB_CHK_ERR(ierr); 00226 00227 run_solution_mesh_optimizer( mesh_and_domain, 00228 err ); // CHECK_IMESH("solution mesh optimizer failed"); 00229 if( err ) return 1; 00230 00231 reps *= 0.01; 00232 } 00233 00234 run_local_smoother2( mesh_and_domain, err, 1e-4 ); // CHECK_IMESH("local smoother2 failed"); 00235 if( err ) return 1; 00236 00237 // run_quality_optimizer( mesh_and_domain, err );MB_CHK_ERR(ierr); 00238 00239 // run_local_smoother( mesh_and_domain, err );MB_CHK_ERR(ierr); 00240 00241 // Delete MOAB instance 00242 delete mesh; 00243 delete plane; 00244 00245 #ifdef MOAB_HAVE_MPI 00246 // MPI_Finalize(); 00247 #endif 00248 00249 return 0; 00250 } 00251 00252 int run_global_smoother( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value ) 00253 { 00254 // double OF_value = 1e-6; 00255 00256 Mesh* mesh = mesh_and_domain.get_mesh(); 00257 MeshDomain* domain = mesh_and_domain.get_domain(); 00258 00259 // creates an intruction queue 00260 InstructionQueue queue1; 00261 00262 // creates a mean ratio quality metric ... 00263 IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio( err ); 00264 // ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric(); 00265 // TMPQualityMetric* mean_ratio = new TMPQualityMetric(); 00266 00267 // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric(); 00268 if( err ) return 1; 00269 mean_ratio->set_averaging_method( QualityMetric::RMS, err ); 00270 if( err ) return 1; 00271 00272 // ... and builds an objective function with it 00273 LPtoPTemplate* obj_func = new LPtoPTemplate( mean_ratio, 1, err ); 00274 // LInfTemplate* obj_func = new LInfTemplate(mean_ratio); 00275 if( err ) return 1; 00276 00277 // creates the feas newt optimization procedures 00278 // ConjugateGradient* pass1 = new ConjugateGradient( obj_func, err ); 00279 // FeasibleNewton* pass1 = new FeasibleNewton( obj_func ); 00280 SteepestDescent* pass1 = new SteepestDescent( obj_func ); 00281 pass1->use_element_on_vertex_patch(); 00282 pass1->do_block_coordinate_descent_optimization(); 00283 pass1->use_global_patch(); 00284 if( err ) return 1; 00285 00286 QualityAssessor stop_qa( mean_ratio ); 00287 00288 // **************Set stopping criterion**************** 00289 TerminationCriterion tc_inner; 00290 tc_inner.add_absolute_vertex_movement( OF_value ); 00291 if( err ) return 1; 00292 TerminationCriterion tc_outer; 00293 tc_inner.add_iteration_limit( 10 ); 00294 tc_outer.add_iteration_limit( 5 ); 00295 tc_outer.set_debug_output_level( 3 ); 00296 tc_inner.set_debug_output_level( 3 ); 00297 pass1->set_inner_termination_criterion( &tc_inner ); 00298 pass1->set_outer_termination_criterion( &tc_outer ); 00299 00300 queue1.add_quality_assessor( &stop_qa, err ); 00301 if( err ) return 1; 00302 00303 // adds 1 pass of pass1 to mesh_set1 00304 queue1.set_master_quality_improver( pass1, err ); 00305 if( err ) return 1; 00306 00307 queue1.add_quality_assessor( &stop_qa, err ); 00308 if( err ) return 1; 00309 00310 // launches optimization on mesh_set 00311 if( domain ) { queue1.run_instructions( &mesh_and_domain, err ); } 00312 else 00313 { 00314 queue1.run_instructions( mesh, err ); 00315 } 00316 if( err ) return 1; 00317 00318 // Construct a MeshTSTT from the file 00319 int ierr = write_vtk_mesh( mesh, "feasible-newton-result.vtk" ); 00320 if( ierr ) return 1; 00321 // MeshWriter::write_vtk(mesh, "feasible-newton-result.vtk", err); 00322 // if (err) return 1; 00323 cout << "Wrote \"feasible-newton-result.vtk\"" << endl; 00324 00325 // print_timing_diagnostics( cout ); 00326 return 0; 00327 } 00328 00329 int write_vtk_mesh( Mesh* mesh, const char* filename ) 00330 { 00331 moab::Interface* mbi = 00332 reinterpret_cast< MBiMesh* >( dynamic_cast< MsqIMesh* >( mesh )->get_imesh_instance() )->mbImpl; 00333 00334 mbi->write_file( filename ); 00335 00336 return 0; 00337 } 00338 00339 int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value ); 00340 int run_local_smoother( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value ) 00341 { 00342 Mesh* mesh = mesh_and_domain.get_mesh(); 00343 moab::Interface* mbi = 00344 reinterpret_cast< MBiMesh* >( dynamic_cast< MsqIMesh* >( mesh )->get_imesh_instance() )->mbImpl; 00345 00346 moab::Tag fixed; 00347 moab::ErrorCode rval = mbi->tag_get_handle( "fixed", 1, moab::MB_TYPE_INTEGER, fixed );MB_CHK_SET_ERR( rval, "Getting tag handle failed" ); 00348 moab::Range cells; 00349 rval = mbi->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "Querying elements failed" ); 00350 00351 moab::LloydSmoother lloyd( mbi, 0, cells, 0, 0 /*fixed*/ ); 00352 00353 lloyd.perform_smooth(); 00354 00355 // run_local_smoother2(mesh_and_domain, err, OF_value); 00356 00357 // Construct a MeshTSTT from the file 00358 int ierr = write_vtk_mesh( mesh, "smart-laplacian-result.vtk" ); 00359 if( ierr ) return 1; 00360 // MeshWriter::write_vtk(mesh, "smart-laplacian-result.vtk", err); 00361 // if (err) return 1; 00362 cout << "Wrote \"smart-laplacian-result.vtk\"" << endl; 00363 return 0; 00364 } 00365 00366 int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value ) 00367 { 00368 // double OF_value = 1e-5; 00369 00370 Mesh* mesh = mesh_and_domain.get_mesh(); 00371 MeshDomain* domain = mesh_and_domain.get_domain(); 00372 00373 // creates an intruction queue 00374 InstructionQueue queue1; 00375 00376 // creates a mean ratio quality metric ... 00377 // IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio(err); 00378 ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric(); 00379 // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric(); 00380 if( err ) return 1; 00381 // mean_ratio->set_gradient_type(QualityMetric::NUMERICAL_GRADIENT); 00382 // mean_ratio->set_hessian_type(QualityMetric::NUMERICAL_HESSIAN); 00383 mean_ratio->set_averaging_method( QualityMetric::RMS, err ); 00384 if( err ) return 1; 00385 00386 // ... and builds an objective function with it 00387 LPtoPTemplate* obj_func = new LPtoPTemplate( mean_ratio, 1, err ); 00388 if( err ) return 1; 00389 00390 if( false ) 00391 { 00392 InstructionQueue qtmp; 00393 Randomize rand( -0.005 ); 00394 TerminationCriterion sc_rand; 00395 sc_rand.add_iteration_limit( 2 ); 00396 rand.set_outer_termination_criterion( &sc_rand ); 00397 qtmp.set_master_quality_improver( &rand, err ); 00398 if( err ) return 1; 00399 if( domain ) { qtmp.run_instructions( &mesh_and_domain, err ); } 00400 else 00401 { 00402 qtmp.run_instructions( mesh, err ); 00403 } 00404 if( err ) return 1; 00405 } 00406 00407 // creates the smart laplacian optimization procedures 00408 SmartLaplacianSmoother* pass1 = new SmartLaplacianSmoother( obj_func ); 00409 // SteepestDescent* pass1 = new SteepestDescent( obj_func ); 00410 00411 QualityAssessor stop_qa( mean_ratio ); 00412 00413 // **************Set stopping criterion**************** 00414 TerminationCriterion tc_inner; 00415 tc_inner.add_absolute_vertex_movement( OF_value ); 00416 TerminationCriterion tc_outer; 00417 tc_outer.add_iteration_limit( 10 ); 00418 pass1->set_inner_termination_criterion( &tc_inner ); 00419 pass1->set_outer_termination_criterion( &tc_outer ); 00420 00421 queue1.add_quality_assessor( &stop_qa, err ); 00422 if( err ) return 1; 00423 00424 // adds 1 pass of pass1 to mesh_set 00425 queue1.set_master_quality_improver( pass1, err ); 00426 if( err ) return 1; 00427 00428 queue1.add_quality_assessor( &stop_qa, err ); 00429 if( err ) return 1; 00430 00431 // launches optimization on mesh_set 00432 if( domain ) { queue1.run_instructions( &mesh_and_domain, err ); } 00433 else 00434 { 00435 queue1.run_instructions( mesh, err ); 00436 } 00437 if( err ) return 1; 00438 00439 // Construct a MeshTSTT from the file 00440 int ierr = write_vtk_mesh( mesh, "smart-laplacian-result.vtk" ); 00441 if( ierr ) return 1; 00442 // MeshWriter::write_vtk(mesh, "smart-laplacian-result.vtk", err); 00443 // if (err) return 1; 00444 cout << "Wrote \"smart-laplacian-result.vtk\"" << endl; 00445 00446 // print_timing_diagnostics( cout ); 00447 return 0; 00448 } 00449 00450 int run_quality_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err ) 00451 { 00452 Mesh* mesh = mesh_and_domain.get_mesh(); 00453 MeshDomain* domain = mesh_and_domain.get_domain(); 00454 00455 // creates an intruction queue 00456 InstructionQueue q; 00457 00458 // do optimization 00459 const double eps = 0.01; 00460 IdealShapeTarget w; 00461 TShapeB1 mu; 00462 TQualityMetric metric( &w, &mu ); 00463 PMeanPTemplate func( 1.0, &metric ); 00464 00465 SteepestDescent solver( &func ); 00466 solver.use_element_on_vertex_patch(); 00467 solver.do_jacobi_optimization(); 00468 00469 TerminationCriterion inner, outer; 00470 inner.add_absolute_vertex_movement( 0.5 * eps ); 00471 outer.add_absolute_vertex_movement( 0.5 * eps ); 00472 00473 QualityAssessor qa( &metric ); 00474 00475 q.add_quality_assessor( &qa, err ); 00476 if( err ) return 1; 00477 q.set_master_quality_improver( &solver, err ); 00478 if( err ) return 1; 00479 q.add_quality_assessor( &qa, err ); 00480 if( err ) return 1; 00481 00482 // launches optimization on mesh_set 00483 if( domain ) { q.run_instructions( &mesh_and_domain, err ); } 00484 else 00485 { 00486 q.run_instructions( mesh, err ); 00487 } 00488 if( err ) return 1; 00489 00490 // Construct a MeshTSTT from the file 00491 int ierr = write_vtk_mesh( mesh, "ideal-shape-result.vtk" ); 00492 if( ierr ) return 1; 00493 // MeshWriter::write_vtk(mesh, "ideal-shape-result.vtk", err); 00494 // if (err) return 1; 00495 cout << "Wrote \"ideal-shape-result.vtk\"" << endl; 00496 00497 print_timing_diagnostics( cout ); 00498 return 0; 00499 } 00500 00501 int run_solution_mesh_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err ) 00502 { 00503 double OF_value = 0.01; 00504 00505 Mesh* mesh = mesh_and_domain.get_mesh(); 00506 MeshDomain* domain = mesh_and_domain.get_domain(); 00507 00508 // creates an intruction queue 00509 InstructionQueue queue1; 00510 00511 // creates a mean ratio quality metric ... 00512 // IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio(err); 00513 ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric(); 00514 // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric(); 00515 // AspectRatioGammaQualityMetric* mean_ratio = new AspectRatioGammaQualityMetric(); 00516 00517 // ElementSolIndQM* solindqm = new ElementSolIndQM(solution_indicator); 00518 // MultiplyQualityMetric* mean_ratio = new MultiplyQualityMetric(new 00519 // VertexConditionNumberQualityMetric(), solindqm, err); 00520 // ElementSolIndQM* mean_ratio = solindqm; 00521 00522 // LocalSizeQualityMetric* mean_ratio = new LocalSizeQualityMetric(); 00523 00524 mean_ratio->set_averaging_method( QualityMetric::SUM_OF_RATIOS_SQUARED, err ); 00525 if( err ) return 1; 00526 00527 // ... and builds an objective function with it 00528 LPtoPTemplate* obj_func = new LPtoPTemplate( mean_ratio, 1, err ); 00529 if( err ) return 1; 00530 00531 // // creates the feas newt optimization procedures 00532 ConjugateGradient* pass1 = new ConjugateGradient( obj_func, err ); 00533 // QuasiNewton* pass1 = new QuasiNewton( obj_func ); 00534 // FeasibleNewton* pass1 = new FeasibleNewton( obj_func ); 00535 pass1->use_global_patch(); 00536 00537 QualityAssessor stop_qa( mean_ratio ); 00538 00539 // **************Set stopping criterion**************** 00540 TerminationCriterion tc_inner; 00541 tc_inner.add_absolute_vertex_movement( OF_value ); 00542 if( err ) return 1; 00543 TerminationCriterion tc_outer; 00544 tc_inner.add_iteration_limit( 20 ); 00545 tc_outer.add_iteration_limit( 5 ); 00546 pass1->set_inner_termination_criterion( &tc_inner ); 00547 pass1->set_outer_termination_criterion( &tc_outer ); 00548 pass1->set_debugging_level( 3 ); 00549 00550 queue1.add_quality_assessor( &stop_qa, err ); 00551 if( err ) return 1; 00552 00553 // adds 1 pass of pass1 to mesh_set1 00554 queue1.set_master_quality_improver( pass1, err ); 00555 if( err ) return 1; 00556 00557 queue1.add_quality_assessor( &stop_qa, err ); 00558 if( err ) return 1; 00559 00560 // launches optimization on mesh_set 00561 if( domain ) { queue1.run_instructions( &mesh_and_domain, err ); } 00562 else 00563 { 00564 queue1.run_instructions( mesh, err ); 00565 } 00566 if( err ) return 1; 00567 00568 // Construct a MeshTSTT from the file 00569 int ierr = write_vtk_mesh( mesh, "solution-mesh-result.vtk" ); 00570 if( ierr ) return 1; 00571 // MeshWriter::write_vtk(mesh, "solution-mesh-result.vtk", err); 00572 // if (err) return 1; 00573 cout << "Wrote \"solution-mesh-result.vtk\"" << endl; 00574 00575 print_timing_diagnostics( cout ); 00576 return 0; 00577 } 00578 00579 int get_imesh_mesh( MBMesquite::Mesh** mesh, const char* file_name, int dimension ) 00580 { 00581 int err; 00582 iMesh_Instance instance = 0; 00583 iMesh_newMesh( NULL, &instance, &err, 0 ); 00584 CHECK_IMESH( "Creation of mesh instance failed" ); 00585 00586 iBase_EntitySetHandle root_set; 00587 iMesh_getRootSet( instance, &root_set, &err ); 00588 CHECK_IMESH( "Could not get root set" ); 00589 00590 iMesh_load( instance, root_set, file_name, 0, &err, strlen( file_name ), 0 ); 00591 CHECK_IMESH( "Could not load mesh from file" ); 00592 00593 iBase_TagHandle fixed_tag; 00594 iMesh_getTagHandle( instance, "fixed", &fixed_tag, &err, strlen( "fixed" ) ); 00595 // if (iBase_SUCCESS != err) 00596 { 00597 // get the skin vertices of those cells and mark them as fixed; we don't want to fix the 00598 // vertices on a part boundary, but since we exchanged a layer of ghost cells, those 00599 // vertices aren't on the skin locally ok to mark non-owned skin vertices too, I won't move 00600 // those anyway use MOAB's skinner class to find the skin 00601 00602 // get all vertices and cells 00603 // make tag to specify fixed vertices, since it's input to the algorithm; use a default 00604 // value of non-fixed so we only need to set the fixed tag for skin vertices 00605 moab::Interface* mbi = reinterpret_cast< MBiMesh* >( instance )->mbImpl; 00606 moab::EntityHandle currset = 0; 00607 moab::Tag fixed; 00608 int def_val = 0; 00609 err = 0; 00610 moab::ErrorCode rval = mbi->tag_get_handle( "fixed", 1, moab::MB_TYPE_INTEGER, fixed, 00611 moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &def_val );MB_CHK_SET_ERR( rval, "Getting tag handle failed" ); 00612 moab::Range verts, cells, skin_verts; 00613 rval = mbi->get_entities_by_type( currset, moab::MBVERTEX, verts );MB_CHK_SET_ERR( rval, "Querying vertices failed" ); 00614 rval = mbi->get_entities_by_dimension( currset, dimension, cells );MB_CHK_SET_ERR( rval, "Querying elements failed" ); 00615 std::cout << "Found " << verts.size() << " vertices and " << cells.size() << " elements" << std::endl; 00616 00617 moab::Skinner skinner( mbi ); 00618 rval = skinner.find_skin( currset, cells, true, skin_verts );MB_CHK_SET_ERR( rval, 00619 "Finding the skin of the mesh failed" ); // 'true' param indicates we want 00620 // vertices back, not cells 00621 00622 std::vector< int > fix_tag( skin_verts.size(), 1 ); // initialized to 1 to indicate fixed 00623 rval = mbi->tag_set_data( fixed, skin_verts, &fix_tag[0] );MB_CHK_SET_ERR( rval, "Setting tag data failed" ); 00624 std::cout << "Found " << skin_verts.size() << " vertices on the skin of the domain." << std::endl; 00625 00626 // fix_tag.resize(verts.size(),0); 00627 // rval = mbi->tag_get_data(fixed, verts, &fix_tag[0]); MB_CHK_SET_ERR(rval, "Getting tag 00628 // data failed"); 00629 00630 iMesh_getTagHandle( instance, "fixed", &fixed_tag, &err, strlen( "fixed" ) ); 00631 CHECK_IMESH( "Getting tag handle (fixed) failed" ); 00632 00633 // Set some arbitrary solution indicator 00634 moab::Tag solindTag; 00635 double def_val_dbl = 0.0; 00636 rval = mbi->tag_get_handle( "solution_indicator", 1, moab::MB_TYPE_DOUBLE, solindTag, 00637 moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &def_val_dbl );MB_CHK_SET_ERR( rval, "Getting tag handle failed" ); 00638 solution_indicator.resize( cells.size(), 0.01 ); 00639 for( unsigned i = 0; i < cells.size() / 4; i++ ) 00640 solution_indicator[i] = 0.1; 00641 for( unsigned i = cells.size() / 4; i < 2 * cells.size() / 4; i++ ) 00642 solution_indicator[i] = 0.5; 00643 for( unsigned i = 2 * cells.size() / 4; i < 3 * cells.size() / 4; i++ ) 00644 solution_indicator[i] = 0.5; 00645 for( unsigned i = 3 * cells.size() / 4; i < cells.size(); i++ ) 00646 solution_indicator[i] = 0.5; 00647 00648 rval = mbi->tag_set_data( solindTag, cells, &solution_indicator[0] );MB_CHK_SET_ERR( rval, "Setting tag data failed" ); 00649 } 00650 00651 MsqError ierr; 00652 MBMesquite::MsqIMesh* imesh = 00653 new MBMesquite::MsqIMesh( instance, root_set, ( dimension == 3 ? iBase_REGION : iBase_FACE ), ierr, 00654 &fixed_tag ); 00655 if( MSQ_CHKERR( ierr ) ) 00656 { 00657 delete imesh; 00658 cerr << err << endl; 00659 err = iBase_FAILURE; 00660 CHECK_IMESH( "Creation of MsqIMesh instance failed" ); 00661 return 0; 00662 } 00663 00664 *mesh = imesh; 00665 return iBase_SUCCESS; 00666 } 00667 00668 int get_native_mesh( MBMesquite::Mesh** mesh, const char* file_name, int ) 00669 { 00670 MsqError err; 00671 MBMesquite::MeshImpl* imesh = new MBMesquite::MeshImpl(); 00672 imesh->read_vtk( file_name, err ); 00673 if( err ) 00674 { 00675 cerr << err << endl; 00676 exit( 3 ); 00677 } 00678 *mesh = imesh; 00679 00680 return iBase_SUCCESS; 00681 } 00682 00683 int get_itaps_domain( MeshDomain** odomain, const char* filename ) 00684 { 00685 00686 if( filename == 0 || strlen( filename ) == 0 ) 00687 { 00688 *odomain = new PlanarDomain( PlanarDomain::XY ); 00689 return 0; 00690 } 00691 00692 int err; 00693 iGeom_Instance geom; 00694 iGeom_newGeom( "", &geom, &err, 0 ); 00695 CHECK_IGEOM( "ERROR: iGeom creation failed" ); 00696 00697 #ifdef MOAB_HAVE_CGM_FACET 00698 FacetModifyEngine::set_modify_enabled( CUBIT_TRUE ); 00699 #endif 00700 00701 iGeom_load( geom, filename, 0, &err, strlen( filename ), 0 ); 00702 CHECK_IGEOM( "Cannot load the geometry" ); 00703 00704 iBase_EntitySetHandle root_set; 00705 iGeom_getRootSet( geom, &root_set, &err ); 00706 CHECK_IGEOM( "getRootSet failed!" ); 00707 00708 // print out the number of entities 00709 std::cout << "Model contents: " << std::endl; 00710 const char* gtype[] = { "vertices: ", "edges: ", "faces: ", "regions: " }; 00711 int nents[4]; 00712 for( int i = 0; i <= 3; ++i ) 00713 { 00714 iGeom_getNumOfType( geom, root_set, i, &nents[i], &err ); 00715 CHECK_IGEOM( "Error: problem getting entities after gLoad." ); 00716 std::cout << gtype[i] << nents[i] << std::endl; 00717 } 00718 00719 iBase_EntityHandle* hd_geom_ents; 00720 int csize = 0, sizealloc = 0; 00721 if( nents[3] > 0 ) 00722 { 00723 hd_geom_ents = (iBase_EntityHandle*)malloc( sizeof( iBase_EntityHandle ) * nents[2] ); 00724 csize = nents[2]; 00725 iGeom_getEntities( geom, root_set, 2, &hd_geom_ents, &csize, &sizealloc, &err ); 00726 } 00727 else 00728 { 00729 hd_geom_ents = (iBase_EntityHandle*)malloc( sizeof( iBase_EntityHandle ) * nents[1] ); 00730 csize = nents[1]; 00731 iGeom_getEntities( geom, root_set, 1, &hd_geom_ents, &csize, &sizealloc, &err ); 00732 } 00733 CHECK_IGEOM( "ERROR: Could not get entities" ); 00734 00735 *odomain = new MsqIGeom( geom, hd_geom_ents[0] ); 00736 return iBase_SUCCESS; 00737 }