MOAB: Mesh Oriented datABase  (version 5.4.1)
high_aspect_ratio_test.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2007 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     retain 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     (2007) [email protected]
00024 
00025   ***************************************************************** */
00026 
00027 /** \file main.cpp
00028  *  \brief Test high aspect ratio case
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "MsqError.hpp"
00034 #include "MeshImpl.hpp"
00035 #include "XYRectangle.hpp"
00036 
00037 #include "ReferenceMesh.hpp"
00038 #include "RefMeshTargetCalculator.hpp"
00039 #include "TShapeB1.hpp"
00040 #include "TQualityMetric.hpp"
00041 
00042 #include "TerminationCriterion.hpp"
00043 #include "QualityAssessor.hpp"
00044 #include "PMeanPTemplate.hpp"
00045 #include "FeasibleNewton.hpp"
00046 #include "ConjugateGradient.hpp"
00047 #include "InstructionQueue.hpp"
00048 
00049 #include <iostream>
00050 #include <fstream>
00051 #include <cstdlib>
00052 #include <cstdio>
00053 
00054 using namespace MBMesquite;
00055 
00056 //! Struct in which to store mesh description
00057 struct MeshParams
00058 {
00059     double x, y, w, h;
00060 };
00061 std::ostream& operator<<( std::ostream& str, const MeshParams& p )
00062 {
00063     return str << p.x << ',' << p.y << ',' << p.w << ',' << p.h;
00064 }
00065 
00066 //! Default values for parameters.
00067 const char default_out_file[] = "high_aspect.vtk";
00068 const MeshParams default_mesh = { 0.6, 0.3, 2.0, 1.0 };
00069 const MeshParams default_ref  = { 0.5, 0.5, 1.0, 1.0 };
00070 
00071 const int INNER_ITERATES = 1;
00072 const int OUTER_ITERATES = 10;
00073 
00074 enum ExitCodes
00075 {
00076     NO_ERROR           = 0,
00077     USAGE_ERROR        = 1,
00078     NOT_AT_TARGET      = 2,
00079     FAR_FROM_TARGET    = 3,
00080     WRONG_DIRECTION    = 4,
00081     DEGENERATE_ELEMENT = 5,
00082     INVERTED_ELEMENT   = 6,
00083     LAST_EXIT_CODE     = INVERTED_ELEMENT
00084 };
00085 
00086 void usage( const char* argv0, bool brief = true )
00087 {
00088     std::ostream& str = brief ? std::cerr : std::cout;
00089 
00090     str << "Usage: " << argv0 << " [-o <output_file>]"
00091         << " [-f|-F] [-t|-T] [-n|-c] [-i <n>]"
00092         << " [-m <x>,<y>[,<w>,<h>]]"
00093         << " [-r <x>,<y>[,<w>,<h>]]" << std::endl;
00094     if( brief )
00095     {
00096         str << "       " << argv0 << " -h" << std::endl;
00097         std::exit( USAGE_ERROR );
00098     }
00099 
00100     str << "  -o  Specify output file (default is \"" << default_out_file << "\")" << std::endl
00101         << "  -f  Fixed boundary vertices" << std::endl
00102         << "  -F  Free boundary vertices (default)" << std::endl
00103         << "  -t  Write VTK timesteps" << std::endl
00104         << "  -T  Write GNUPlot timesteps" << std::endl
00105         << "  -m  Specify input mesh parameters (default " << default_mesh << ")" << std::endl
00106         << "  -r  Specify reference mesh parameters (default " << default_ref << ")" << std::endl
00107         << "  -n  Use FeasibleNewton solver" << std::endl
00108         << "  -c  Use ConjugateGradient solver (default)" << std::endl
00109         << "  -i  Specify number of iterations (default:" << OUTER_ITERATES << ")" << std::endl
00110         << std::endl;
00111 
00112     std::exit( NO_ERROR );
00113 }
00114 
00115 #define CHECKERR                                                                                        \
00116     if( err )                                                                                           \
00117     {                                                                                                   \
00118         std::cerr << "Internal error at line " << __LINE__ << ":" << std::endl << ( err ) << std::endl; \
00119         std::exit( LAST_EXIT_CODE + ( err ).error_code() );                                             \
00120     }
00121 
00122 /*    |<----- x ----->|
00123  *   (6)-------------(7)-------------(8)--
00124  *    |               |               | ^
00125  *    |               |               | |
00126  *    |      [2]      |      [3]      | |
00127  *    |               |               | |
00128  *    |               |               | |
00129  * --(3)-------------(4)-------------(5)h
00130  *  ^ |               |               | |
00131  *  | |               |               | |
00132  *  y |      [0]      |      [1]      | |
00133  *  | |               |               | |
00134  *  v |               |               | v
00135  * --(0)-------------(1)-------------(2)--
00136  *    |<------------- w ------------->|
00137  *
00138  * z = 0
00139  */
00140 void create_input_mesh( const MeshParams& params, bool all_fixed, MeshImpl& mesh, MsqError& err );
00141 
00142 void parse_options( char* argv[],
00143                     int argc,
00144                     MeshParams& mesh,
00145                     MeshParams& ref,
00146                     std::string& output_file,
00147                     bool& fixed_boundary,
00148                     TerminationCriterion::TimeStepFileType& write_timesteps,
00149                     bool& use_feas_newt,
00150                     int& num_iterations );
00151 
00152 std::string base_name( std::string filename );
00153 
00154 int main( int argc, char* argv[] )
00155 {
00156     MeshParams input_params, reference_params;
00157     bool fixed_boundary_vertices, feas_newt_solver;
00158     TerminationCriterion::TimeStepFileType write_timestep_files;
00159     std::string output_file_name;
00160     int num_iterations;
00161 
00162     parse_options( argv, argc, input_params, reference_params, output_file_name, fixed_boundary_vertices,
00163                    write_timestep_files, feas_newt_solver, num_iterations );
00164 
00165     MsqError err;
00166     MeshImpl mesh, refmesh;
00167     XYRectangle domain( input_params.w, input_params.h );
00168     create_input_mesh( input_params, fixed_boundary_vertices, mesh, err );
00169     CHECKERR
00170     create_input_mesh( reference_params, fixed_boundary_vertices, refmesh, err );
00171     CHECKERR
00172     domain.setup( &mesh, err );
00173     CHECKERR
00174 
00175     ReferenceMesh rmesh( &refmesh );
00176     RefMeshTargetCalculator tc( &rmesh );
00177     TShapeB1 tm;
00178     TQualityMetric qm( &tc, &tm );
00179 
00180     PMeanPTemplate of( 1.0, &qm );
00181     ConjugateGradient cg( &of );
00182     cg.use_element_on_vertex_patch();
00183     FeasibleNewton fn( &of );
00184     fn.use_element_on_vertex_patch();
00185     VertexMover* solver = feas_newt_solver ? (VertexMover*)&fn : (VertexMover*)&cg;
00186 
00187     TerminationCriterion inner, outer;
00188     inner.add_iteration_limit( INNER_ITERATES );
00189     outer.add_iteration_limit( num_iterations );
00190     if( write_timestep_files != TerminationCriterion::NOTYPE )
00191         outer.write_mesh_steps( base_name( output_file_name ).c_str(), write_timestep_files );
00192     solver->set_inner_termination_criterion( &inner );
00193     solver->set_outer_termination_criterion( &outer );
00194 
00195     QualityAssessor qa( &qm );
00196     InstructionQueue q;
00197     q.add_quality_assessor( &qa, err );
00198     q.set_master_quality_improver( solver, err );
00199     q.add_quality_assessor( &qa, err );
00200 
00201     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &domain );
00202     q.run_instructions( &mesh_and_domain, err );
00203     CHECKERR
00204 
00205     mesh.write_vtk( output_file_name.c_str(), err );
00206     CHECKERR
00207 
00208     // check for inverted elements
00209     int inv, unk;
00210     qa.get_inverted_element_count( inv, unk, err );
00211     if( inv )
00212     {
00213         std::cerr << inv << " inverted elements in final mesh" << std::endl;
00214         return INVERTED_ELEMENT;
00215     }
00216     else if( unk )
00217     {
00218         std::cerr << unk << " degenerate elements in final mesh" << std::endl;
00219         return DEGENERATE_ELEMENT;
00220     }
00221 
00222     // find the free vertex
00223     std::vector< Mesh::VertexHandle > vertices;
00224     mesh.get_all_vertices( vertices, err );
00225     if( vertices.empty() )
00226     {
00227         std::cerr << "Mesh contains no vertices" << std::endl;
00228         return USAGE_ERROR;
00229     }
00230     std::vector< unsigned short > dof( vertices.size(), 0 );
00231     domain.domain_DoF( arrptr( vertices ), arrptr( dof ), vertices.size(), err );
00232     CHECKERR
00233     int idx                              = std::find( dof.begin(), dof.end(), 2 ) - dof.begin();
00234     const Mesh::VertexHandle free_vertex = vertices[idx];
00235     MsqVertex coords;
00236     mesh.vertices_get_coordinates( &free_vertex, &coords, 1, err );
00237     CHECKERR
00238 
00239     // calculate optimal position for vertex
00240     const double xf = reference_params.x / reference_params.w;
00241     const double yf = reference_params.y / reference_params.h;
00242     Vector3D expect( xf * input_params.w, yf * input_params.h, 0 );
00243 
00244     // Test that we aren't further from the expected location
00245     // than when we started.
00246     const Vector3D init( input_params.x, input_params.y, 0 );
00247     if( ( coords - expect ).length() > ( init - expect ).length() )
00248     {
00249         std::cerr << "Vertex moved away from expected optimal position: "
00250                   << "(" << coords[0] << ", " << coords[1] << std::endl;
00251         return WRONG_DIRECTION;
00252     }
00253 
00254     // check if vertex moved MIN_FRACT of the way from the original position
00255     // to the desired one in the allowed iterations
00256     const double MIN_FRACT = 0.2;  // 20% of the way in 10 iterations
00257     const double fract     = ( coords - init ).length() / ( expect - init ).length();
00258     if( fract < MIN_FRACT )
00259     {
00260         std::cerr << "Vertex far from optimimal location" << std::endl
00261                   << "  Expected: (" << expect[0] << ", " << expect[1] << ", " << expect[2] << ")" << std::endl
00262                   << "  Actual:   (" << coords[0] << ", " << coords[1] << ", " << coords[2] << ")" << std::endl;
00263         return FAR_FROM_TARGET;
00264     }
00265 
00266     // check if vertex is at destired location
00267     const double EPS = 5e-2;  // allow a little leway
00268     if( fabs( coords[0] - expect[0] ) > EPS * input_params.w || fabs( coords[1] - expect[1] ) > EPS * input_params.h ||
00269         fabs( expect[2] ) > EPS )
00270     {
00271         std::cerr << "Vertex not at optimimal location" << std::endl
00272                   << "  Expected: (" << expect[0] << ", " << expect[1] << ", " << expect[2] << ")" << std::endl
00273                   << "  Actual:   (" << coords[0] << ", " << coords[1] << ", " << coords[2] << ")" << std::endl;
00274         return NOT_AT_TARGET;
00275     }
00276 
00277     return 0;
00278 }
00279 
00280 void parse_mesh_params( const char* argv, const char* arg, MeshParams& result )
00281 {
00282     int c = std::sscanf( arg, "%lf,%lf,%lf,%lf", &result.x, &result.y, &result.w, &result.h );
00283     if( c != 2 && c != 4 )
00284     {
00285         std::cerr << "Error parsing mesh dimensions: \"" << arg << '"' << std::endl;
00286         usage( argv );
00287     }
00288 }
00289 
00290 enum ParseState
00291 {
00292     OPEN,
00293     EXPECTING_M,
00294     EXPECTING_R,
00295     EXPECTING_O,
00296     EXPECTING_I
00297 };
00298 void parse_options( char* argv[],
00299                     int argc,
00300                     MeshParams& mesh,
00301                     MeshParams& ref,
00302                     std::string& output_file,
00303                     bool& fixed_boundary,
00304                     TerminationCriterion::TimeStepFileType& write_timesteps,
00305                     bool& feas_newt_solver,
00306                     int& num_iterations )
00307 {
00308     // begin with defaults
00309     mesh             = default_mesh;
00310     ref              = default_ref;
00311     output_file      = default_out_file;
00312     fixed_boundary   = false;
00313     write_timesteps  = TerminationCriterion::NOTYPE;
00314     feas_newt_solver = false;
00315     num_iterations   = OUTER_ITERATES;
00316 
00317     // parse CLI args
00318     ParseState state = OPEN;
00319     for( int i = 1; i < argc; ++i )
00320     {
00321         switch( state )
00322         {
00323             case EXPECTING_M:
00324                 parse_mesh_params( argv[0], argv[i], mesh );
00325                 state = OPEN;
00326                 break;
00327             case EXPECTING_R:
00328                 parse_mesh_params( argv[0], argv[i], ref );
00329                 state = OPEN;
00330                 break;
00331             case EXPECTING_O:
00332                 output_file = argv[i];
00333                 state       = OPEN;
00334                 break;
00335             case EXPECTING_I:
00336                 num_iterations = atoi( argv[i] );
00337                 state          = OPEN;
00338                 break;
00339             case OPEN:
00340                 if( argv[i][0] != '-' || argv[i][1] == '\0' || argv[i][2] != '\0' )
00341                 {
00342                     std::cerr << "Unexpected argument: \"" << argv[i] << '"' << std::endl;
00343                     usage( argv[0] );
00344                 }
00345 
00346                 switch( argv[i][1] )
00347                 {
00348                     default:
00349                         usage( argv[0], true );
00350                         break;
00351                     case 'h':
00352                         usage( argv[0], false );
00353                         break;
00354                     case 'o':
00355                         state = EXPECTING_O;
00356                         break;
00357                     case 'f':
00358                         fixed_boundary = true;
00359                         break;
00360                     case 'F':
00361                         fixed_boundary = false;
00362                         break;
00363                     case 't':
00364                         write_timesteps = TerminationCriterion::VTK;
00365                         break;
00366                     case 'T':
00367                         write_timesteps = TerminationCriterion::GNUPLOT;
00368                         break;
00369                     case 'm':
00370                         state = EXPECTING_M;
00371                         break;
00372                     case 'r':
00373                         state = EXPECTING_R;
00374                         break;
00375                     case 'n':
00376                         feas_newt_solver = true;
00377                         break;
00378                     case 'c':
00379                         feas_newt_solver = false;
00380                         break;
00381                     case 'i':
00382                         state = EXPECTING_I;
00383                         break;
00384                 }
00385                 break;
00386         }
00387     }
00388 }
00389 
00390 const char* temp_file = "high_aspect_input.vtk";
00391 void create_input_mesh( const MeshParams& p, bool all_fixed, MeshImpl& mesh, MsqError& err )
00392 {
00393     const double z = 0;
00394     const int F    = all_fixed;
00395     std::ofstream vtkfile( temp_file );
00396     double bx = all_fixed ? 0.5 * p.w : p.x;
00397     double by = all_fixed ? 0.5 * p.h : p.y;
00398     vtkfile << "# vtk DataFile Version 3.0" << std::endl
00399             << "Mesquite High Aspect Ratio test" << std::endl
00400             << "ASCII" << std::endl
00401             << "DATASET UNSTRUCTURED_GRID" << std::endl
00402             << "POINTS 9 float" << std::endl
00403             << 0.0 << ' ' << 0.0 << ' ' << z << std::endl
00404             << bx << ' ' << 0.0 << ' ' << z << std::endl
00405             << p.w << ' ' << 0.0 << ' ' << z << std::endl
00406             << 0.0 << ' ' << by << ' ' << z << std::endl
00407             << p.x << ' ' << p.y << ' ' << z << std::endl
00408             << p.w << ' ' << by << ' ' << z << std::endl
00409             << 0.0 << ' ' << p.h << ' ' << z << std::endl
00410             << bx << ' ' << p.h << ' ' << z << std::endl
00411             << p.w << ' ' << p.h << ' ' << z << std::endl
00412             << "CELLS 4 20" << std::endl
00413             << "4 0 1 4 3" << std::endl
00414             << "4 1 2 5 4" << std::endl
00415             << "4 4 5 8 7" << std::endl
00416             << "4 3 4 7 6" << std::endl
00417             << "CELL_TYPES 4" << std::endl
00418             << "9 9 9 9" << std::endl
00419             << "POINT_DATA 9" << std::endl
00420             << "SCALARS fixed int" << std::endl
00421             << "LOOKUP_TABLE default" << std::endl
00422             << "1 " << F << " 1" << std::endl
00423             << F << " 0 " << F << std::endl
00424             << "1 " << F << " 1" << std::endl;
00425     vtkfile.close();
00426     mesh.read_vtk( temp_file, err );
00427     std::remove( temp_file );MSQ_CHKERR( err );
00428 }
00429 
00430 std::string base_name( std::string filename )
00431 {
00432     std::string::size_type i = filename.rfind( "." );
00433     if( !i || i == std::string::npos )
00434         return filename;
00435     else
00436         return filename.substr( 0, i );
00437 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines