MOAB: Mesh Oriented datABase  (version 5.3.0)
synchronous_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 syncronous boundary case
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "ConditionNumberQualityMetric.hpp"
00034 #include "IdealWeightInverseMeanRatio.hpp"
00035 #include "LPtoPTemplate.hpp"
00036 #include "InstructionQueue.hpp"
00037 #include "ConjugateGradient.hpp"
00038 #include "FeasibleNewton.hpp"
00039 #include "MeshImpl.hpp"
00040 #include "QualityAssessor.hpp"
00041 #include "XYRectangle.hpp"
00042 
00043 #include <iostream>
00044 #include <fstream>
00045 #include <cstdlib>
00046 
00047 using namespace MBMesquite;
00048 
00049 // constants
00050 const double min_x = 0.0, max_x = 2.0;
00051 const double min_y = 0.0, mid_y = 1.0, max_y = 2.0;
00052 const double z = 0.0;
00053 
00054 const char default_out_file[] = "synchronous.vtk";
00055 double default_x              = 0.25;
00056 
00057 void create_input_mesh( double mid_x, MBMesquite::MeshImpl& mesh, MBMesquite::MsqError& );
00058 
00059 void usage( const char* argv0, bool brief = true )
00060 {
00061     std::ostream& str = brief ? std::cerr : std::cout;
00062 
00063     str << "Usage: " << argv0 << " [-x <coord>]"
00064         << " [-j|-n|-d]"
00065         << " [-r|-c]"
00066         << " [<output_file>]" << std::endl;
00067     if( brief )
00068     {
00069         str << "       " << argv0 << " -h" << std::endl;
00070         std::exit( 1 );
00071     }
00072 
00073     str << "  -x  Specify X coordinate value for mesh (default is " << default_x << ")" << std::endl
00074         << "  -j  Use ConjugateGradient solver (default)" << std::endl
00075         << "  -n  Use FeasibleNewton solver" << std::endl
00076         << "  -r  Use IdealWeightInverseMeanRation metric" << std::endl
00077         << "  -c  Use ConditionNumber metric (default)" << std::endl
00078         << "Default output file is \"" << default_out_file << '"' << std::endl;
00079 
00080     std::exit( 0 );
00081 }
00082 
00083 char mSolver = '\0', mMetric = '\0';
00084 const char* outputFile = default_out_file;
00085 double input_x         = default_x;
00086 
00087 void parse_options( char* argv[], int argc )
00088 {
00089     bool next_arg_is_x = false;
00090     for( int i = 1; i < argc; ++i )
00091     {
00092         if( next_arg_is_x )
00093         {
00094             next_arg_is_x = false;
00095             char* endptr  = 0;
00096             input_x       = std::strtod( argv[i], &endptr );
00097             if( endptr && *endptr ) usage( argv[0] );
00098             continue;
00099         }
00100 
00101         if( argv[i][0] != '-' )
00102         {
00103             if( outputFile != default_out_file ) usage( argv[0] );
00104             outputFile = argv[i];
00105             continue;
00106         }
00107 
00108         for( const char* p = argv[i] + 1; *p; ++p )
00109         {
00110             switch( *p )
00111             {
00112                 case 'x':
00113                     next_arg_is_x = true;
00114                     break;
00115 
00116                 case 'j':
00117                 case 'n':
00118                     if( mSolver ) usage( argv[0] );
00119                     mSolver = *p;
00120                     break;
00121 
00122                 case 'r':
00123                 case 'c':
00124                     if( mMetric ) usage( argv[0] );
00125                     mMetric = *p;
00126                     break;
00127 
00128                 default:
00129                     usage( argv[0], *p != 'h' );
00130                     break;
00131             }
00132         }
00133     }
00134 
00135     if( next_arg_is_x ) usage( argv[0] );
00136 
00137     // default values
00138     if( !mMetric ) mMetric = 'c';
00139     if( !mSolver ) mSolver = 'j';
00140 }
00141 
00142 int main( int argc, char* argv[] )
00143 {
00144     parse_options( argv, argc );
00145 
00146     MeshImpl mesh;
00147     XYRectangle domain( max_x - min_x, max_y - min_y, min_x, min_y );
00148     MsqError err;
00149 
00150     create_input_mesh( input_x, mesh, err );
00151     if( MSQ_CHKERR( err ) )
00152     {
00153         std::cerr << err << std::endl;
00154         return 2;
00155     }
00156 
00157     domain.setup( &mesh, err );
00158     if( MSQ_CHKERR( err ) )
00159     {
00160         std::cerr << err << std::endl;
00161         return 2;
00162     }
00163 
00164     QualityMetric* metric = 0;
00165     if( mMetric == 'c' )
00166         metric = new ConditionNumberQualityMetric;
00167     else
00168         metric = new IdealWeightInverseMeanRatio;
00169 
00170     LPtoPTemplate function( 1, metric );
00171 
00172     VertexMover* solver = 0;
00173     if( mSolver == 'j' )
00174         solver = new ConjugateGradient( &function );
00175     else
00176         solver = new FeasibleNewton( &function );
00177 
00178     if( PatchSetUser* psu = dynamic_cast< PatchSetUser* >( solver ) ) psu->use_global_patch();
00179 
00180     TerminationCriterion inner;
00181     inner.add_absolute_vertex_movement( 1e-4 );
00182     inner.write_mesh_steps( "synchronous", TerminationCriterion::GNUPLOT );
00183     solver->set_inner_termination_criterion( &inner );
00184 
00185     InstructionQueue q;
00186     QualityAssessor qa( metric, 10 );
00187     q.add_quality_assessor( &qa, err );
00188     q.set_master_quality_improver( solver, err );
00189     q.add_quality_assessor( &qa, err );
00190     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &domain );
00191     q.run_instructions( &mesh_and_domain, err );
00192     delete solver;
00193     delete metric;
00194 
00195     if( MSQ_CHKERR( err ) )
00196     {
00197         std::cerr << err << std::endl;
00198         return 3;
00199     }
00200 
00201     mesh.write_vtk( outputFile, err );
00202     if( MSQ_CHKERR( err ) )
00203     {
00204         std::cerr << err << std::endl;
00205         return 2;
00206     }
00207 
00208     return 0;
00209 }
00210 
00211 const char* temp_file = "syncrononous_input.vtk";
00212 void create_input_mesh( double mid_x, MeshImpl& mesh, MsqError& err )
00213 {
00214     std::ofstream vtkfile( temp_file );
00215     vtkfile << "# vtk DataFile Version 3.0" << std::endl
00216             << "Mesquite Syncronous Boundary test" << std::endl
00217             << "ASCII" << std::endl
00218             << "DATASET UNSTRUCTURED_GRID" << std::endl
00219             << "POINTS 9 float" << std::endl
00220             << min_x << ' ' << max_y << ' ' << z << std::endl
00221             << mid_x << ' ' << max_y << ' ' << z << std::endl
00222             << max_x << ' ' << max_y << ' ' << z << std::endl
00223             << min_x << ' ' << mid_y << ' ' << z << std::endl
00224             << mid_x << ' ' << mid_y << ' ' << z << std::endl
00225             << max_x << ' ' << mid_y << ' ' << z << std::endl
00226             << min_x << ' ' << min_y << ' ' << z << std::endl
00227             << mid_x << ' ' << min_y << ' ' << z << std::endl
00228             << max_x << ' ' << min_y << ' ' << z << std::endl
00229             << "CELLS 4 20" << std::endl
00230             << "4 1 0 3 4" << std::endl
00231             << "4 2 1 4 5" << std::endl
00232             << "4 4 3 6 7" << std::endl
00233             << "4 5 4 7 8" << std::endl
00234             << "CELL_TYPES 4" << std::endl
00235             << "9 9 9 9" << std::endl
00236             << "POINT_DATA 9" << std::endl
00237             << "SCALARS fixed int" << std::endl
00238             << "LOOKUP_TABLE default" << std::endl
00239             << "1 0 1" << std::endl
00240             << "1 0 1" << std::endl
00241             << "1 0 1" << std::endl;
00242     vtkfile.close();
00243 
00244     mesh.read_vtk( temp_file, err );
00245     remove( temp_file );
00246 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines