MOAB: Mesh Oriented datABase  (version 5.2.1)
OptimizeMeshMesquite.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines