|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
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 }
1.7.6.1 
