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