MOAB: Mesh Oriented datABase  (version 5.2.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) kraftche@cae.wisc.edu
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[], int argc, MeshParams& mesh, MeshParams& ref, std::string& output_file,
00143                     bool& fixed_boundary, TerminationCriterion::TimeStepFileType& write_timesteps, bool& use_feas_newt,
00144                     int& num_iterations );
00145 
00146 std::string base_name( std::string filename );
00147 
00148 int main( int argc, char* argv[] )
00149 {
00150     MeshParams input_params, reference_params;
00151     bool fixed_boundary_vertices, feas_newt_solver;
00152     TerminationCriterion::TimeStepFileType write_timestep_files;
00153     std::string output_file_name;
00154     int num_iterations;
00155 
00156     parse_options( argv, argc, input_params, reference_params, output_file_name, fixed_boundary_vertices,
00157                    write_timestep_files, feas_newt_solver, num_iterations );
00158 
00159     MsqError err;
00160     MeshImpl mesh, refmesh;
00161     XYRectangle domain( input_params.w, input_params.h );
00162     create_input_mesh( input_params, fixed_boundary_vertices, mesh, err );
00163     CHECKERR
00164     create_input_mesh( reference_params, fixed_boundary_vertices, refmesh, err );
00165     CHECKERR
00166     domain.setup( &mesh, err );
00167     CHECKERR
00168 
00169     ReferenceMesh rmesh( &refmesh );
00170     RefMeshTargetCalculator tc( &rmesh );
00171     TShapeB1 tm;
00172     TQualityMetric qm( &tc, &tm );
00173 
00174     PMeanPTemplate of( 1.0, &qm );
00175     ConjugateGradient cg( &of );
00176     cg.use_element_on_vertex_patch();
00177     FeasibleNewton fn( &of );
00178     fn.use_element_on_vertex_patch();
00179     VertexMover* solver = feas_newt_solver ? (VertexMover*)&fn : (VertexMover*)&cg;
00180 
00181     TerminationCriterion inner, outer;
00182     inner.add_iteration_limit( INNER_ITERATES );
00183     outer.add_iteration_limit( num_iterations );
00184     if( write_timestep_files != TerminationCriterion::NOTYPE )
00185         outer.write_mesh_steps( base_name( output_file_name ).c_str(), write_timestep_files );
00186     solver->set_inner_termination_criterion( &inner );
00187     solver->set_outer_termination_criterion( &outer );
00188 
00189     QualityAssessor qa( &qm );
00190     InstructionQueue q;
00191     q.add_quality_assessor( &qa, err );
00192     q.set_master_quality_improver( solver, err );
00193     q.add_quality_assessor( &qa, err );
00194 
00195     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &domain );
00196     q.run_instructions( &mesh_and_domain, err );
00197     CHECKERR
00198 
00199     mesh.write_vtk( output_file_name.c_str(), err );
00200     CHECKERR
00201 
00202     // check for inverted elements
00203     int inv, unk;
00204     qa.get_inverted_element_count( inv, unk, err );
00205     if( inv )
00206     {
00207         std::cerr << inv << " inverted elements in final mesh" << std::endl;
00208         return INVERTED_ELEMENT;
00209     }
00210     else if( unk )
00211     {
00212         std::cerr << unk << " degenerate elements in final mesh" << std::endl;
00213         return DEGENERATE_ELEMENT;
00214     }
00215 
00216     // find the free vertex
00217     std::vector< Mesh::VertexHandle > vertices;
00218     mesh.get_all_vertices( vertices, err );
00219     if( vertices.empty() )
00220     {
00221         std::cerr << "Mesh contains no vertices" << std::endl;
00222         return USAGE_ERROR;
00223     }
00224     std::vector< unsigned short > dof( vertices.size(), 0 );
00225     domain.domain_DoF( arrptr( vertices ), arrptr( dof ), vertices.size(), err );
00226     CHECKERR
00227     int idx                              = std::find( dof.begin(), dof.end(), 2 ) - dof.begin();
00228     const Mesh::VertexHandle free_vertex = vertices[idx];
00229     MsqVertex coords;
00230     mesh.vertices_get_coordinates( &free_vertex, &coords, 1, err );
00231     CHECKERR
00232 
00233     // calculate optimal position for vertex
00234     const double xf = reference_params.x / reference_params.w;
00235     const double yf = reference_params.y / reference_params.h;
00236     Vector3D expect( xf * input_params.w, yf * input_params.h, 0 );
00237 
00238     // Test that we aren't further from the expected location
00239     // than when we started.
00240     const Vector3D init( input_params.x, input_params.y, 0 );
00241     if( ( coords - expect ).length() > ( init - expect ).length() )
00242     {
00243         std::cerr << "Vertex moved away from expected optimal position: "
00244                   << "(" << coords[0] << ", " << coords[1] << std::endl;
00245         return WRONG_DIRECTION;
00246     }
00247 
00248     // check if vertex moved MIN_FRACT of the way from the original position
00249     // to the desired one in the allowed iterations
00250     const double MIN_FRACT = 0.2;  // 20% of the way in 10 iterations
00251     const double fract     = ( coords - init ).length() / ( expect - init ).length();
00252     if( fract < MIN_FRACT )
00253     {
00254         std::cerr << "Vertex far from optimimal location" << std::endl
00255                   << "  Expected: (" << expect[0] << ", " << expect[1] << ", " << expect[2] << ")" << std::endl
00256                   << "  Actual:   (" << coords[0] << ", " << coords[1] << ", " << coords[2] << ")" << std::endl;
00257         return FAR_FROM_TARGET;
00258     }
00259 
00260     // check if vertex is at destired location
00261     const double EPS = 5e-2;  // allow a little leway
00262     if( fabs( coords[0] - expect[0] ) > EPS * input_params.w || fabs( coords[1] - expect[1] ) > EPS * input_params.h ||
00263         fabs( expect[2] ) > EPS )
00264     {
00265         std::cerr << "Vertex not at optimimal location" << std::endl
00266                   << "  Expected: (" << expect[0] << ", " << expect[1] << ", " << expect[2] << ")" << std::endl
00267                   << "  Actual:   (" << coords[0] << ", " << coords[1] << ", " << coords[2] << ")" << std::endl;
00268         return NOT_AT_TARGET;
00269     }
00270 
00271     return 0;
00272 }
00273 
00274 void parse_mesh_params( const char* argv, const char* arg, MeshParams& result )
00275 {
00276     int c = std::sscanf( arg, "%lf,%lf,%lf,%lf", &result.x, &result.y, &result.w, &result.h );
00277     if( c != 2 && c != 4 )
00278     {
00279         std::cerr << "Error parsing mesh dimensions: \"" << arg << '"' << std::endl;
00280         usage( argv );
00281     }
00282 }
00283 
00284 enum ParseState
00285 {
00286     OPEN,
00287     EXPECTING_M,
00288     EXPECTING_R,
00289     EXPECTING_O,
00290     EXPECTING_I
00291 };
00292 void parse_options( char* argv[], int argc, MeshParams& mesh, MeshParams& ref, std::string& output_file,
00293                     bool& fixed_boundary, TerminationCriterion::TimeStepFileType& write_timesteps,
00294                     bool& feas_newt_solver, int& num_iterations )
00295 {
00296     // begin with defaults
00297     mesh             = default_mesh;
00298     ref              = default_ref;
00299     output_file      = default_out_file;
00300     fixed_boundary   = false;
00301     write_timesteps  = TerminationCriterion::NOTYPE;
00302     feas_newt_solver = false;
00303     num_iterations   = OUTER_ITERATES;
00304 
00305     // parse CLI args
00306     ParseState state = OPEN;
00307     for( int i = 1; i < argc; ++i )
00308     {
00309         switch( state )
00310         {
00311             case EXPECTING_M:
00312                 parse_mesh_params( argv[0], argv[i], mesh );
00313                 state = OPEN;
00314                 break;
00315             case EXPECTING_R:
00316                 parse_mesh_params( argv[0], argv[i], ref );
00317                 state = OPEN;
00318                 break;
00319             case EXPECTING_O:
00320                 output_file = argv[i];
00321                 state       = OPEN;
00322                 break;
00323             case EXPECTING_I:
00324                 num_iterations = atoi( argv[i] );
00325                 state          = OPEN;
00326                 break;
00327             case OPEN:
00328                 if( argv[i][0] != '-' || argv[i][1] == '\0' || argv[i][2] != '\0' )
00329                 {
00330                     std::cerr << "Unexpected argument: \"" << argv[i] << '"' << std::endl;
00331                     usage( argv[0] );
00332                 }
00333 
00334                 switch( argv[i][1] )
00335                 {
00336                     default:
00337                         usage( argv[0], true );
00338                         break;
00339                     case 'h':
00340                         usage( argv[0], false );
00341                         break;
00342                     case 'o':
00343                         state = EXPECTING_O;
00344                         break;
00345                     case 'f':
00346                         fixed_boundary = true;
00347                         break;
00348                     case 'F':
00349                         fixed_boundary = false;
00350                         break;
00351                     case 't':
00352                         write_timesteps = TerminationCriterion::VTK;
00353                         break;
00354                     case 'T':
00355                         write_timesteps = TerminationCriterion::GNUPLOT;
00356                         break;
00357                     case 'm':
00358                         state = EXPECTING_M;
00359                         break;
00360                     case 'r':
00361                         state = EXPECTING_R;
00362                         break;
00363                     case 'n':
00364                         feas_newt_solver = true;
00365                         break;
00366                     case 'c':
00367                         feas_newt_solver = false;
00368                         break;
00369                     case 'i':
00370                         state = EXPECTING_I;
00371                         break;
00372                 }
00373                 break;
00374         }
00375     }
00376 }
00377 
00378 const char* temp_file = "high_aspect_input.vtk";
00379 void create_input_mesh( const MeshParams& p, bool all_fixed, MeshImpl& mesh, MsqError& err )
00380 {
00381     const double z = 0;
00382     const int F    = all_fixed;
00383     std::ofstream vtkfile( temp_file );
00384     double bx = all_fixed ? 0.5 * p.w : p.x;
00385     double by = all_fixed ? 0.5 * p.h : p.y;
00386     vtkfile << "# vtk DataFile Version 3.0" << std::endl
00387             << "Mesquite High Aspect Ratio test" << std::endl
00388             << "ASCII" << std::endl
00389             << "DATASET UNSTRUCTURED_GRID" << std::endl
00390             << "POINTS 9 float" << std::endl
00391             << 0.0 << ' ' << 0.0 << ' ' << z << std::endl
00392             << bx << ' ' << 0.0 << ' ' << z << std::endl
00393             << p.w << ' ' << 0.0 << ' ' << z << std::endl
00394             << 0.0 << ' ' << by << ' ' << z << std::endl
00395             << p.x << ' ' << p.y << ' ' << z << std::endl
00396             << p.w << ' ' << by << ' ' << z << std::endl
00397             << 0.0 << ' ' << p.h << ' ' << z << std::endl
00398             << bx << ' ' << p.h << ' ' << z << std::endl
00399             << p.w << ' ' << p.h << ' ' << z << std::endl
00400             << "CELLS 4 20" << std::endl
00401             << "4 0 1 4 3" << std::endl
00402             << "4 1 2 5 4" << std::endl
00403             << "4 4 5 8 7" << std::endl
00404             << "4 3 4 7 6" << std::endl
00405             << "CELL_TYPES 4" << std::endl
00406             << "9 9 9 9" << std::endl
00407             << "POINT_DATA 9" << std::endl
00408             << "SCALARS fixed int" << std::endl
00409             << "LOOKUP_TABLE default" << std::endl
00410             << "1 " << F << " 1" << std::endl
00411             << F << " 0 " << F << std::endl
00412             << "1 " << F << " 1" << std::endl;
00413     vtkfile.close();
00414     mesh.read_vtk( temp_file, err );
00415     std::remove( temp_file );MSQ_CHKERR( err );
00416 }
00417 
00418 std::string base_name( std::string filename )
00419 {
00420     std::string::size_type i = filename.rfind( "." );
00421     if( !i || i == std::string::npos )
00422         return filename;
00423     else
00424         return filename.substr( 0, i );
00425 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines