MOAB: Mesh Oriented datABase  (version 5.2.1)
2d_target_test.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retian certain rights to this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     (2006) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 #include "Mesquite.hpp"
00028 #include <iostream>
00029 using std::cerr;
00030 using std::cout;
00031 using std::endl;
00032 using std::ostream;
00033 #include <memory>
00034 using std::auto_ptr;
00035 #include <ctype.h>
00036 
00037 #include "MeshImpl.hpp"
00038 #include "MsqError.hpp"
00039 #include "InstructionQueue.hpp"
00040 #include "TerminationCriterion.hpp"
00041 #include "QualityAssessor.hpp"
00042 #include "PMeanPTemplate.hpp"
00043 #include "PatchPowerMeanP.hpp"
00044 #include "ConjugateGradient.hpp"
00045 #include "PlanarDomain.hpp"
00046 #include "IdealShapeTarget.hpp"
00047 #include "ConditionNumberQualityMetric.hpp"
00048 #include "ReferenceMesh.hpp"
00049 #include "RefMeshTargetCalculator.hpp"
00050 #include "TestUtil.hpp"
00051 
00052 #include "TQualityMetric.hpp"
00053 #include "ElementPMeanP.hpp"
00054 #include "VertexPMeanP.hpp"
00055 
00056 #include "TSquared.hpp"
00057 #include "TShapeNB1.hpp"
00058 #include "TShapeB1.hpp"
00059 #include "TShapeOrientNB1.hpp"
00060 #include "TShapeOrientNB2.hpp"
00061 #include "TShapeOrientB1.hpp"
00062 #include "TShapeSize2DNB1.hpp"
00063 #include "TShapeSizeB1.hpp"
00064 #include "TShapeSize2DB2.hpp"
00065 #include "TShapeSizeB3.hpp"
00066 #include "TShapeSizeOrientNB1.hpp"
00067 #include "TShapeSizeOrientB1.hpp"
00068 #include "TShapeSizeOrientB2.hpp"
00069 
00070 using namespace MBMesquite;
00071 
00072 static const struct
00073 {
00074     TMetric* u;
00075     const char* n;
00076 } metrics[] = { { new TSquared, "TSquared" },
00077                 { new TShapeNB1, "Shape" },
00078                 { new TShapeB1, "ShapeBarrier" },
00079                 { new TShapeOrientNB1, "ShapeOrient1" },
00080                 { new TShapeOrientNB2, "ShapeOrient2" },
00081                 { new TShapeOrientB1, "ShapeOrientBarrier" },
00082                 { new TShapeSize2DNB1, "ShapeSize" },
00083                 { new TShapeSizeB1, "ShapeSizeBarrier1" },
00084                 { new TShapeSize2DB2, "ShapeSizeBarrier2" },
00085                 { new TShapeSizeB3, "ShapeSizeBarrier3" },
00086                 { new TShapeSizeOrientB1, "ShapeSizeOrient1" },
00087                 { new TShapeSizeOrientB1, "ShapeSizeOrientBarrier1" },
00088                 { new TShapeSizeOrientB2, "ShapeSizeOrientBarrier2" },
00089                 { 0, 0 } };
00090 
00091 enum AveragingScheme
00092 {
00093     NONE = 0,
00094     ELEMENT,
00095     VERTEX,
00096     PATCH
00097 };
00098 const char* const averaging_names[] = { "none", "element", "vertex", "patch" };
00099 
00100 // default values
00101 const double DEFAULT_OF_POWER            = 1.0;
00102 const unsigned DEFAULT_METRIC_IDX        = 0;
00103 const AveragingScheme DEFAULT_AVG_SCHEME = NONE;
00104 std::string DEFAULT_INPUT_FILE           = TestDir + "/2D/vtk/quads/untangled/quads_4by2_bad.vtk";
00105 std::string DEFAULT_OUTPUT_FILE          = "./out.vtk";
00106 
00107 static PlanarDomain make_domain( Mesh* mesh, MsqError& );
00108 
00109 static int do_smoother( const char* input_file, const char* output_file, const char* ref_mesh_file, double of_power,
00110                         unsigned metric_idx, AveragingScheme avg_scheme )
00111 {
00112     MsqPrintError err( cerr );
00113 
00114     TMetric* const target_metric = metrics[metric_idx].u;
00115     cout << "Input file:  " << input_file << endl;
00116     cout << "Metric:      ";
00117     if( avg_scheme != NONE ) cout << averaging_names[avg_scheme] << " average of ";
00118     cout << metrics[metric_idx].n << endl;
00119     cout << "Of Power:    " << of_power << endl;
00120 
00121     auto_ptr< TargetCalculator > tc;
00122     auto_ptr< MeshImpl > ref_mesh_impl;
00123     auto_ptr< ReferenceMesh > ref_mesh;
00124     if( ref_mesh_file )
00125     {
00126         ref_mesh_impl.reset( new MeshImpl );
00127         ref_mesh_impl->read_vtk( ref_mesh_file, err );
00128         if( MSQ_CHKERR( err ) ) return 2;
00129         ref_mesh.reset( new ReferenceMesh( ref_mesh_impl.get() ) );
00130         tc.reset( new RefMeshTargetCalculator( ref_mesh.get() ) );
00131     }
00132     else
00133     {
00134         tc.reset( new IdealShapeTarget() );
00135     }
00136 
00137     TQualityMetric jacobian_metric( tc.get(), target_metric );
00138     ElementPMeanP elem_avg( of_power, &jacobian_metric );
00139     VertexPMeanP vtx_avg( of_power, &jacobian_metric );
00140     QualityMetric* mmetrics[] = { &jacobian_metric, &elem_avg, &vtx_avg, &jacobian_metric };
00141     QualityMetric* metric     = mmetrics[avg_scheme];
00142 
00143     TerminationCriterion outer, inner;
00144     outer.add_iteration_limit( 1 );
00145     inner.add_absolute_vertex_movement( 1e-4 );
00146     inner.add_iteration_limit( 100 );
00147 
00148     PMeanPTemplate obj1( of_power, metric );
00149     PatchPowerMeanP obj2( of_power, metric );
00150     ObjectiveFunction& objective = *( ( avg_scheme == PATCH ) ? (ObjectiveFunction*)&obj2 : (ObjectiveFunction*)&obj1 );
00151 
00152     ConjugateGradient solver( &objective, err );
00153     if( MSQ_CHKERR( err ) ) return 1;
00154     solver.set_inner_termination_criterion( &inner );
00155     solver.set_outer_termination_criterion( &outer );
00156     solver.use_global_patch();
00157 
00158     ConditionNumberQualityMetric qm_metric;
00159     QualityAssessor before_assessor;
00160     QualityAssessor after_assessor;
00161     before_assessor.add_quality_assessment( metric, 10 );
00162     before_assessor.add_quality_assessment( &qm_metric );
00163     after_assessor.add_quality_assessment( metric, 10 );
00164 
00165     InstructionQueue q;
00166     q.add_quality_assessor( &before_assessor, err );
00167     q.set_master_quality_improver( &solver, err );
00168     q.add_quality_assessor( &after_assessor, err );
00169 
00170     MeshImpl mesh;
00171     mesh.read_vtk( input_file, err );
00172     if( MSQ_CHKERR( err ) ) return 2;
00173     PlanarDomain geom = make_domain( &mesh, err );
00174     if( MSQ_CHKERR( err ) ) return 1;
00175     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &geom );
00176     q.run_instructions( &mesh_and_domain, err );
00177     if( MSQ_CHKERR( err ) ) return 3;
00178     mesh.write_vtk( output_file, err );
00179     if( MSQ_CHKERR( err ) ) return 2;
00180     cout << "Wrote: " << output_file << endl;
00181 
00182     before_assessor.scale_histograms( &after_assessor );
00183 
00184     return 0;
00185 }
00186 
00187 static PlanarDomain make_domain( Mesh* mesh, MsqError& err )
00188 {
00189     // calculate bounding box of mesh vertices
00190     Vector3D minimum( HUGE_VAL, HUGE_VAL, HUGE_VAL );
00191     Vector3D maximum( -HUGE_VAL, -HUGE_VAL, -HUGE_VAL );
00192     std::vector< Mesh::VertexHandle > vertices;
00193     mesh->get_all_vertices( vertices, err );
00194     if( MSQ_CHKERR( err ) ) { return PlanarDomain( minimum, maximum ); }
00195     if( vertices.empty() )
00196     {
00197         std::cerr << "Mesh contains no vertices" << std::endl;
00198         exit( 1 );
00199     }
00200     std::vector< MsqVertex > coords( vertices.size() );
00201     mesh->vertices_get_coordinates( arrptr( vertices ), arrptr( coords ), vertices.size(), err );
00202     if( MSQ_CHKERR( err ) ) { return PlanarDomain( minimum, maximum ); }
00203     std::vector< MsqVertex >::const_iterator i;
00204     for( i = coords.begin(); i != coords.end(); ++i )
00205     {
00206         const MsqVertex& v = *i;
00207         for( unsigned j = 0; j < 3; ++j )
00208         {
00209             if( v[j] < minimum[j] ) minimum[j] = v[j];
00210             if( v[j] > maximum[j] ) maximum[j] = v[j];
00211         }
00212     }
00213 
00214     // Look for a "zero" plane
00215     int k;
00216     maximum -= minimum;
00217     for( k = 2; k >= 0 && maximum[k] > 1e-6; --k )
00218         ;
00219     if( k < 0 )
00220     {
00221         MSQ_SETERR( err )( "Cannot determine plane of mesh.", MsqError::INVALID_STATE );
00222         return PlanarDomain( minimum, maximum );
00223     }
00224 
00225     Vector3D point( 0.0, 0.0, 0.0 ), normal( 0.0, 0.0, 0.0 );
00226     normal[k] = 1.0;
00227     point[k]  = minimum[k];
00228     return PlanarDomain( normal, point );
00229 }
00230 
00231 static void usage( const char* argv0, bool help = false )
00232 {
00233     ostream& str = help ? cout : cerr;
00234     str << argv0 << " [-p <power>] [-m metric_name] [-a averaging] [-e]"
00235         << " -r [ref_mesh] [input_file] [output_file]" << endl
00236         << argv0 << " <-l|-h>" << std::endl;
00237     if( help )
00238     {
00239         str << "     -p  : specify exponent value for p-mean^p OF template (default: " << DEFAULT_OF_POWER << ")"
00240             << endl
00241             << "     -m  : specify 2D target metric to use (default: " << metrics[DEFAULT_METRIC_IDX].n << ")" << endl
00242             << "     -a  : specify metric averaging scheme (default: " << averaging_names[DEFAULT_AVG_SCHEME] << ")"
00243             << endl
00244             << "              (none,vertex,element,patch)" << endl
00245             << "     -e  : sample at mid-edge points (default is corners only)" << endl
00246             << "     -r  : use reference mesh instead of ideal elements for targets" << endl
00247             << "     -l  : list available metrics" << endl
00248             << "     -h  : this help output" << endl
00249             << " default input file:  " << DEFAULT_INPUT_FILE << endl
00250             << " default output file: " << DEFAULT_OUTPUT_FILE << endl
00251             << endl;
00252     }
00253     exit( !help );
00254 }
00255 
00256 static void check_next_arg( int argc, char* argv[], int& i )
00257 {
00258     if( i == argc )
00259     {
00260         cerr << "Expected argument following \"" << argv[i] << '"' << endl;
00261         usage( argv[0] );
00262     }
00263     ++i;
00264 }
00265 
00266 static double parse_double( int argc, char* argv[], int& i )
00267 {
00268     check_next_arg( argc, argv, i );
00269     char* endptr;
00270     double result = strtod( argv[i], &endptr );
00271     if( endptr == argv[i] || *endptr )
00272     {
00273         cerr << "Expected real number following \"" << argv[i - 1] << '"' << endl;
00274         usage( argv[0] );
00275     }
00276     return result;
00277 }
00278 
00279 static int comp_string_start( const char* p, const char* s )
00280 {
00281     int i;
00282     for( i = 0; p[i]; ++i )
00283         if( tolower( p[i] ) != tolower( s[i] ) ) return 0;
00284     return s[i] ? -1 : 1;
00285 }
00286 
00287 static AveragingScheme parse_averaging( int argc, char* argv[], int& i )
00288 {
00289     check_next_arg( argc, argv, i );
00290     for( unsigned j = 0; j < 4; ++j )
00291         if( comp_string_start( argv[i], averaging_names[j] ) ) return (AveragingScheme)j;
00292     cerr << "Expected one of { ";
00293     for( unsigned j = 0; j < 4; ++j )
00294         cerr << '"' << averaging_names[j] << "\" ";
00295     cerr << "} following \"" << argv[i - 1] << '"' << endl;
00296     usage( argv[0] );
00297     return NONE;
00298 }
00299 
00300 static unsigned parse_metric( int argc, char* argv[], int& i )
00301 {
00302     check_next_arg( argc, argv, i );
00303     unsigned part = 0, all = 0, count = 0, have_all = 0;
00304     for( unsigned j = 0; metrics[j].u; ++j )
00305     {
00306         if( unsigned k = comp_string_start( argv[i], metrics[j].n ) )
00307         {
00308             if( k > 0 )
00309             {
00310                 all      = j;
00311                 have_all = 1;
00312             }
00313             else
00314             {
00315                 part = j;
00316                 ++count;
00317             }
00318         }
00319     }
00320 
00321     if( have_all ) return all;
00322 
00323     if( count )
00324     {
00325         if( count == 1 ) return part;
00326         cerr << "Ambiguous metric name: \"" << argv[i] << '"' << endl;
00327         usage( argv[0] );
00328     }
00329 
00330     cerr << "Invalid metric name following \"" << argv[i - 1] << "\" option" << endl;
00331     usage( argv[0] );
00332     return (unsigned)-1;
00333 }
00334 
00335 static void list_metrics()
00336 {
00337     cout << "Available metrics:" << endl;
00338     for( unsigned j = 0; metrics[j].u; ++j )
00339         cout << "\t" << metrics[j].n << endl;
00340     exit( 0 );
00341 }
00342 
00343 int main( int argc, char* argv[] )
00344 {
00345     MsqPrintError err( cout );
00346 
00347     // CL settings
00348     double of_power            = DEFAULT_OF_POWER;
00349     unsigned metric_idx        = DEFAULT_METRIC_IDX;
00350     AveragingScheme avg_scheme = DEFAULT_AVG_SCHEME;
00351     const char* input_file     = 0;
00352     const char* output_file    = 0;
00353     const char* ref_mesh_file  = 0;
00354 
00355     bool proc_opts = true;
00356     for( int i = 1; i < argc; ++i )
00357     {
00358         if( !proc_opts || argv[i][0] != '-' )
00359         {
00360             if( output_file )
00361             {
00362                 cerr << "Unexpected file name: \"" << argv[i] << '"' << endl;
00363                 usage( argv[0] );
00364             }
00365             else if( input_file )
00366                 output_file = argv[i];
00367             else
00368                 input_file = argv[i];
00369             continue;
00370         }
00371 
00372         if( !argv[i][1] || argv[i][2] )
00373         {
00374             cerr << "Invalid option: \"" << argv[i] << '"' << endl;
00375             usage( argv[0] );
00376         }
00377 
00378         switch( argv[i][1] )
00379         {
00380             case 'p':
00381                 of_power = parse_double( argc, argv, i );
00382                 break;
00383             case 'm':
00384                 metric_idx = parse_metric( argc, argv, i );
00385                 break;
00386             case 'a':
00387                 avg_scheme = parse_averaging( argc, argv, i );
00388                 break;
00389             case 'r':
00390                 check_next_arg( argc, argv, i );
00391                 ref_mesh_file = argv[i];
00392                 break;
00393             case '-':
00394                 proc_opts = false;
00395                 break;
00396             case 'h':
00397                 usage( argv[0], true );
00398                 break;
00399             case 'l':
00400                 list_metrics();
00401                 break;
00402             default:
00403                 cerr << "Invalid option: \"" << argv[i] << '"' << endl;
00404                 usage( argv[0] );
00405         }
00406     }
00407 
00408     if( !input_file ) input_file = DEFAULT_INPUT_FILE.c_str();
00409     if( !output_file ) output_file = DEFAULT_OUTPUT_FILE.c_str();
00410 
00411     return do_smoother( input_file, output_file, ref_mesh_file, of_power, metric_idx, avg_scheme );
00412 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines