Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
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 )
00312     {
00313         queue1.run_instructions( &mesh_and_domain, err );
00314     }
00315     else
00316     {
00317         queue1.run_instructions( mesh, err );
00318     }
00319     if( err ) return 1;
00320 
00321     // Construct a MeshTSTT from the file
00322     int ierr = write_vtk_mesh( mesh, "feasible-newton-result.vtk" );
00323     if( ierr ) return 1;
00324     // MeshWriter::write_vtk(mesh, "feasible-newton-result.vtk", err);
00325     // if (err) return 1;
00326     cout << "Wrote \"feasible-newton-result.vtk\"" << endl;
00327 
00328     // print_timing_diagnostics( cout );
00329     return 0;
00330 }
00331 
00332 int write_vtk_mesh( Mesh* mesh, const char* filename )
00333 {
00334     moab::Interface* mbi =
00335         reinterpret_cast< MBiMesh* >( dynamic_cast< MsqIMesh* >( mesh )->get_imesh_instance() )->mbImpl;
00336 
00337     mbi->write_file( filename );
00338 
00339     return 0;
00340 }
00341 
00342 int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value );
00343 int run_local_smoother( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value )
00344 {
00345     Mesh* mesh = mesh_and_domain.get_mesh();
00346     moab::Interface* mbi =
00347         reinterpret_cast< MBiMesh* >( dynamic_cast< MsqIMesh* >( mesh )->get_imesh_instance() )->mbImpl;
00348 
00349     moab::Tag fixed;
00350     moab::ErrorCode rval = mbi->tag_get_handle( "fixed", 1, moab::MB_TYPE_INTEGER, fixed );MB_CHK_SET_ERR( rval, "Getting tag handle failed" );
00351     moab::Range cells;
00352     rval = mbi->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "Querying elements failed" );
00353 
00354     moab::LloydSmoother lloyd( mbi, 0, cells, 0, 0 /*fixed*/ );
00355 
00356     lloyd.perform_smooth();
00357 
00358     // run_local_smoother2(mesh_and_domain, err, OF_value);
00359 
00360     // Construct a MeshTSTT from the file
00361     int ierr = write_vtk_mesh( mesh, "smart-laplacian-result.vtk" );
00362     if( ierr ) return 1;
00363     // MeshWriter::write_vtk(mesh, "smart-laplacian-result.vtk", err);
00364     // if (err) return 1;
00365     cout << "Wrote \"smart-laplacian-result.vtk\"" << endl;
00366     return 0;
00367 }
00368 
00369 int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value )
00370 {
00371     // double OF_value = 1e-5;
00372 
00373     Mesh* mesh         = mesh_and_domain.get_mesh();
00374     MeshDomain* domain = mesh_and_domain.get_domain();
00375 
00376     // creates an intruction queue
00377     InstructionQueue queue1;
00378 
00379     // creates a mean ratio quality metric ...
00380     // IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio(err);
00381     ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric();
00382     // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric();
00383     if( err ) return 1;
00384     //  mean_ratio->set_gradient_type(QualityMetric::NUMERICAL_GRADIENT);
00385     //   mean_ratio->set_hessian_type(QualityMetric::NUMERICAL_HESSIAN);
00386     mean_ratio->set_averaging_method( QualityMetric::RMS, err );
00387     if( err ) return 1;
00388 
00389     // ... and builds an objective function with it
00390     LPtoPTemplate* obj_func = new LPtoPTemplate( mean_ratio, 1, err );
00391     if( err ) return 1;
00392 
00393     if( false )
00394     {
00395         InstructionQueue qtmp;
00396         Randomize rand( -0.005 );
00397         TerminationCriterion sc_rand;
00398         sc_rand.add_iteration_limit( 2 );
00399         rand.set_outer_termination_criterion( &sc_rand );
00400         qtmp.set_master_quality_improver( &rand, err );
00401         if( err ) return 1;
00402         if( domain )
00403         {
00404             qtmp.run_instructions( &mesh_and_domain, err );
00405         }
00406         else
00407         {
00408             qtmp.run_instructions( mesh, err );
00409         }
00410         if( err ) return 1;
00411     }
00412 
00413     // creates the smart laplacian optimization procedures
00414     SmartLaplacianSmoother* pass1 = new SmartLaplacianSmoother( obj_func );
00415     // SteepestDescent* pass1 = new SteepestDescent( obj_func );
00416 
00417     QualityAssessor stop_qa( mean_ratio );
00418 
00419     // **************Set stopping criterion****************
00420     TerminationCriterion tc_inner;
00421     tc_inner.add_absolute_vertex_movement( OF_value );
00422     TerminationCriterion tc_outer;
00423     tc_outer.add_iteration_limit( 10 );
00424     pass1->set_inner_termination_criterion( &tc_inner );
00425     pass1->set_outer_termination_criterion( &tc_outer );
00426 
00427     queue1.add_quality_assessor( &stop_qa, err );
00428     if( err ) return 1;
00429 
00430     // adds 1 pass of pass1 to mesh_set
00431     queue1.set_master_quality_improver( pass1, err );
00432     if( err ) return 1;
00433 
00434     queue1.add_quality_assessor( &stop_qa, err );
00435     if( err ) return 1;
00436 
00437     // launches optimization on mesh_set
00438     if( domain )
00439     {
00440         queue1.run_instructions( &mesh_and_domain, err );
00441     }
00442     else
00443     {
00444         queue1.run_instructions( mesh, err );
00445     }
00446     if( err ) return 1;
00447 
00448     // Construct a MeshTSTT from the file
00449     int ierr = write_vtk_mesh( mesh, "smart-laplacian-result.vtk" );
00450     if( ierr ) return 1;
00451     // MeshWriter::write_vtk(mesh, "smart-laplacian-result.vtk", err);
00452     // if (err) return 1;
00453     cout << "Wrote \"smart-laplacian-result.vtk\"" << endl;
00454 
00455     // print_timing_diagnostics( cout );
00456     return 0;
00457 }
00458 
00459 int run_quality_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err )
00460 {
00461     Mesh* mesh         = mesh_and_domain.get_mesh();
00462     MeshDomain* domain = mesh_and_domain.get_domain();
00463 
00464     // creates an intruction queue
00465     InstructionQueue q;
00466 
00467     // do optimization
00468     const double eps = 0.01;
00469     IdealShapeTarget w;
00470     TShapeB1 mu;
00471     TQualityMetric metric( &w, &mu );
00472     PMeanPTemplate func( 1.0, &metric );
00473 
00474     SteepestDescent solver( &func );
00475     solver.use_element_on_vertex_patch();
00476     solver.do_jacobi_optimization();
00477 
00478     TerminationCriterion inner, outer;
00479     inner.add_absolute_vertex_movement( 0.5 * eps );
00480     outer.add_absolute_vertex_movement( 0.5 * eps );
00481 
00482     QualityAssessor qa( &metric );
00483 
00484     q.add_quality_assessor( &qa, err );
00485     if( err ) return 1;
00486     q.set_master_quality_improver( &solver, err );
00487     if( err ) return 1;
00488     q.add_quality_assessor( &qa, err );
00489     if( err ) return 1;
00490 
00491     // launches optimization on mesh_set
00492     if( domain )
00493     {
00494         q.run_instructions( &mesh_and_domain, err );
00495     }
00496     else
00497     {
00498         q.run_instructions( mesh, err );
00499     }
00500     if( err ) return 1;
00501 
00502     // Construct a MeshTSTT from the file
00503     int ierr = write_vtk_mesh( mesh, "ideal-shape-result.vtk" );
00504     if( ierr ) return 1;
00505     // MeshWriter::write_vtk(mesh, "ideal-shape-result.vtk", err);
00506     // if (err) return 1;
00507     cout << "Wrote \"ideal-shape-result.vtk\"" << endl;
00508 
00509     print_timing_diagnostics( cout );
00510     return 0;
00511 }
00512 
00513 int run_solution_mesh_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err )
00514 {
00515     double OF_value = 0.01;
00516 
00517     Mesh* mesh         = mesh_and_domain.get_mesh();
00518     MeshDomain* domain = mesh_and_domain.get_domain();
00519 
00520     // creates an intruction queue
00521     InstructionQueue queue1;
00522 
00523     // creates a mean ratio quality metric ...
00524     // IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio(err);
00525     ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric();
00526     // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric();
00527     // AspectRatioGammaQualityMetric* mean_ratio = new AspectRatioGammaQualityMetric();
00528 
00529     // ElementSolIndQM* solindqm = new ElementSolIndQM(solution_indicator);
00530     // MultiplyQualityMetric* mean_ratio = new MultiplyQualityMetric(new
00531     // VertexConditionNumberQualityMetric(), solindqm, err);
00532     // ElementSolIndQM* mean_ratio = solindqm;
00533 
00534     // LocalSizeQualityMetric* mean_ratio = new LocalSizeQualityMetric();
00535 
00536     mean_ratio->set_averaging_method( QualityMetric::SUM_OF_RATIOS_SQUARED, err );
00537     if( err ) return 1;
00538 
00539     // ... and builds an objective function with it
00540     LPtoPTemplate* obj_func = new LPtoPTemplate( mean_ratio, 1, err );
00541     if( err ) return 1;
00542 
00543     // // creates the feas newt optimization procedures
00544     ConjugateGradient* pass1 = new ConjugateGradient( obj_func, err );
00545     // QuasiNewton* pass1 = new QuasiNewton( obj_func );
00546     // FeasibleNewton* pass1 = new FeasibleNewton( obj_func );
00547     pass1->use_global_patch();
00548 
00549     QualityAssessor stop_qa( mean_ratio );
00550 
00551     // **************Set stopping criterion****************
00552     TerminationCriterion tc_inner;
00553     tc_inner.add_absolute_vertex_movement( OF_value );
00554     if( err ) return 1;
00555     TerminationCriterion tc_outer;
00556     tc_inner.add_iteration_limit( 20 );
00557     tc_outer.add_iteration_limit( 5 );
00558     pass1->set_inner_termination_criterion( &tc_inner );
00559     pass1->set_outer_termination_criterion( &tc_outer );
00560     pass1->set_debugging_level( 3 );
00561 
00562     queue1.add_quality_assessor( &stop_qa, err );
00563     if( err ) return 1;
00564 
00565     // adds 1 pass of pass1 to mesh_set1
00566     queue1.set_master_quality_improver( pass1, err );
00567     if( err ) return 1;
00568 
00569     queue1.add_quality_assessor( &stop_qa, err );
00570     if( err ) return 1;
00571 
00572     // launches optimization on mesh_set
00573     if( domain )
00574     {
00575         queue1.run_instructions( &mesh_and_domain, err );
00576     }
00577     else
00578     {
00579         queue1.run_instructions( mesh, err );
00580     }
00581     if( err ) return 1;
00582 
00583     // Construct a MeshTSTT from the file
00584     int ierr = write_vtk_mesh( mesh, "solution-mesh-result.vtk" );
00585     if( ierr ) return 1;
00586     // MeshWriter::write_vtk(mesh, "solution-mesh-result.vtk", err);
00587     // if (err) return 1;
00588     cout << "Wrote \"solution-mesh-result.vtk\"" << endl;
00589 
00590     print_timing_diagnostics( cout );
00591     return 0;
00592 }
00593 
00594 int get_imesh_mesh( MBMesquite::Mesh** mesh, const char* file_name, int dimension )
00595 {
00596     int err;
00597     iMesh_Instance instance = 0;
00598     iMesh_newMesh( NULL, &instance, &err, 0 );
00599     CHECK_IMESH( "Creation of mesh instance failed" );
00600 
00601     iBase_EntitySetHandle root_set;
00602     iMesh_getRootSet( instance, &root_set, &err );
00603     CHECK_IMESH( "Could not get root set" );
00604 
00605     iMesh_load( instance, root_set, file_name, 0, &err, strlen( file_name ), 0 );
00606     CHECK_IMESH( "Could not load mesh from file" );
00607 
00608     iBase_TagHandle fixed_tag;
00609     iMesh_getTagHandle( instance, "fixed", &fixed_tag, &err, strlen( "fixed" ) );
00610     // if (iBase_SUCCESS != err)
00611     {
00612         // get the skin vertices of those cells and mark them as fixed; we don't want to fix the
00613         // vertices on a part boundary, but since we exchanged a layer of ghost cells, those
00614         // vertices aren't on the skin locally ok to mark non-owned skin vertices too, I won't move
00615         // those anyway use MOAB's skinner class to find the skin
00616 
00617         // get all vertices and cells
00618         // make tag to specify fixed vertices, since it's input to the algorithm; use a default
00619         // value of non-fixed so we only need to set the fixed tag for skin vertices
00620         moab::Interface* mbi       = reinterpret_cast< MBiMesh* >( instance )->mbImpl;
00621         moab::EntityHandle currset = 0;
00622         moab::Tag fixed;
00623         int def_val          = 0;
00624         err                  = 0;
00625         moab::ErrorCode rval = mbi->tag_get_handle( "fixed", 1, moab::MB_TYPE_INTEGER, fixed,
00626                                                     moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &def_val );MB_CHK_SET_ERR( rval, "Getting tag handle failed" );
00627         moab::Range verts, cells, skin_verts;
00628         rval = mbi->get_entities_by_type( currset, moab::MBVERTEX, verts );MB_CHK_SET_ERR( rval, "Querying vertices failed" );
00629         rval = mbi->get_entities_by_dimension( currset, dimension, cells );MB_CHK_SET_ERR( rval, "Querying elements failed" );
00630         std::cout << "Found " << verts.size() << " vertices and " << cells.size() << " elements" << std::endl;
00631 
00632         moab::Skinner skinner( mbi );
00633         rval = skinner.find_skin( currset, cells, true, skin_verts );MB_CHK_SET_ERR( rval,
00634                         "Finding the skin of the mesh failed" );  // 'true' param indicates we want
00635                                                                   // vertices back, not cells
00636 
00637         std::vector< int > fix_tag( skin_verts.size(), 1 );  // initialized to 1 to indicate fixed
00638         rval = mbi->tag_set_data( fixed, skin_verts, &fix_tag[0] );MB_CHK_SET_ERR( rval, "Setting tag data failed" );
00639         std::cout << "Found " << skin_verts.size() << " vertices on the skin of the domain." << std::endl;
00640 
00641         // fix_tag.resize(verts.size(),0);
00642         // rval = mbi->tag_get_data(fixed, verts, &fix_tag[0]); MB_CHK_SET_ERR(rval, "Getting tag
00643         // data failed");
00644 
00645         iMesh_getTagHandle( instance, "fixed", &fixed_tag, &err, strlen( "fixed" ) );
00646         CHECK_IMESH( "Getting tag handle (fixed) failed" );
00647 
00648         // Set some arbitrary solution indicator
00649         moab::Tag solindTag;
00650         double def_val_dbl = 0.0;
00651         rval               = mbi->tag_get_handle( "solution_indicator", 1, moab::MB_TYPE_DOUBLE, solindTag,
00652                                                   moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &def_val_dbl );MB_CHK_SET_ERR( rval, "Getting tag handle failed" );
00653         solution_indicator.resize( cells.size(), 0.01 );
00654         for( unsigned i = 0; i < cells.size() / 4; i++ )
00655             solution_indicator[i] = 0.1;
00656         for( unsigned i = cells.size() / 4; i < 2 * cells.size() / 4; i++ )
00657             solution_indicator[i] = 0.5;
00658         for( unsigned i = 2 * cells.size() / 4; i < 3 * cells.size() / 4; i++ )
00659             solution_indicator[i] = 0.5;
00660         for( unsigned i = 3 * cells.size() / 4; i < cells.size(); i++ )
00661             solution_indicator[i] = 0.5;
00662 
00663         rval = mbi->tag_set_data( solindTag, cells, &solution_indicator[0] );MB_CHK_SET_ERR( rval, "Setting tag data failed" );
00664     }
00665 
00666     MsqError ierr;
00667     MBMesquite::MsqIMesh* imesh =
00668         new MBMesquite::MsqIMesh( instance, root_set, ( dimension == 3 ? iBase_REGION : iBase_FACE ), ierr,
00669                                   &fixed_tag );
00670     if( MSQ_CHKERR( ierr ) )
00671     {
00672         delete imesh;
00673         cerr << err << endl;
00674         err = iBase_FAILURE;
00675         CHECK_IMESH( "Creation of MsqIMesh instance failed" );
00676         return 0;
00677     }
00678 
00679     *mesh = imesh;
00680     return iBase_SUCCESS;
00681 }
00682 
00683 int get_native_mesh( MBMesquite::Mesh** mesh, const char* file_name, int )
00684 {
00685     MsqError err;
00686     MBMesquite::MeshImpl* imesh = new MBMesquite::MeshImpl();
00687     imesh->read_vtk( file_name, err );
00688     if( err )
00689     {
00690         cerr << err << endl;
00691         exit( 3 );
00692     }
00693     *mesh = imesh;
00694 
00695     return iBase_SUCCESS;
00696 }
00697 
00698 int get_itaps_domain( MeshDomain** odomain, const char* filename )
00699 {
00700 
00701     if( filename == 0 || strlen( filename ) == 0 )
00702     {
00703         *odomain = new PlanarDomain( PlanarDomain::XY );
00704         return 0;
00705     }
00706 
00707     int err;
00708     iGeom_Instance geom;
00709     iGeom_newGeom( "", &geom, &err, 0 );
00710     CHECK_IGEOM( "ERROR: iGeom creation failed" );
00711 
00712 #ifdef MOAB_HAVE_CGM_FACET
00713     FacetModifyEngine::set_modify_enabled( CUBIT_TRUE );
00714 #endif
00715 
00716     iGeom_load( geom, filename, 0, &err, strlen( filename ), 0 );
00717     CHECK_IGEOM( "Cannot load the geometry" );
00718 
00719     iBase_EntitySetHandle root_set;
00720     iGeom_getRootSet( geom, &root_set, &err );
00721     CHECK_IGEOM( "getRootSet failed!" );
00722 
00723     // print out the number of entities
00724     std::cout << "Model contents: " << std::endl;
00725     const char* gtype[] = { "vertices: ", "edges: ", "faces: ", "regions: " };
00726     int nents[4];
00727     for( int i = 0; i <= 3; ++i )
00728     {
00729         iGeom_getNumOfType( geom, root_set, i, &nents[i], &err );
00730         CHECK_IGEOM( "Error: problem getting entities after gLoad." );
00731         std::cout << gtype[i] << nents[i] << std::endl;
00732     }
00733 
00734     iBase_EntityHandle* hd_geom_ents;
00735     int csize = 0, sizealloc = 0;
00736     if( nents[3] > 0 )
00737     {
00738         hd_geom_ents = (iBase_EntityHandle*)malloc( sizeof( iBase_EntityHandle ) * nents[2] );
00739         csize        = nents[2];
00740         iGeom_getEntities( geom, root_set, 2, &hd_geom_ents, &csize, &sizealloc, &err );
00741     }
00742     else
00743     {
00744         hd_geom_ents = (iBase_EntityHandle*)malloc( sizeof( iBase_EntityHandle ) * nents[1] );
00745         csize        = nents[1];
00746         iGeom_getEntities( geom, root_set, 1, &hd_geom_ents, &csize, &sizealloc, &err );
00747     }
00748     CHECK_IGEOM( "ERROR: Could not get entities" );
00749 
00750     *odomain = new MsqIGeom( geom, hd_geom_ents[0] );
00751     return iBase_SUCCESS;
00752 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines