MOAB: Mesh Oriented datABase  (version 5.2.1)
2d_metrics_test.cpp
Go to the documentation of this file.
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     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     (2006) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file main.cpp
00028  *  \brief Implement all experiments from 2D Metrics Paper
00029  *  \author Jason Kraftcheck
00030  */
00031 #include "TestUtil.hpp"
00032 #include "Mesquite.hpp"
00033 #include "MeshImpl.hpp"
00034 #include "MsqError.hpp"
00035 #include "MsqVertex.hpp"
00036 #include "PMeanPTemplate.hpp"
00037 #include "ConjugateGradient.hpp"
00038 #include "FeasibleNewton.hpp"
00039 #include "PlanarDomain.hpp"
00040 #include "XYPlanarDomain.hpp"
00041 #include "ReferenceMesh.hpp"
00042 #include "RefMeshTargetCalculator.hpp"
00043 #include "TQualityMetric.hpp"
00044 #include "InstructionQueue.hpp"
00045 #include "IdealShapeTarget.hpp"
00046 #include "MeshWriter.hpp"
00047 
00048 #include "TSquared.hpp"
00049 #include "TShapeNB1.hpp"
00050 #include "TShapeB1.hpp"
00051 #include "TShapeOrientNB1.hpp"
00052 #include "TShapeOrientNB2.hpp"
00053 #include "TShapeOrientB1.hpp"
00054 #include "TShapeSize2DNB1.hpp"
00055 #include "TShapeSizeB1.hpp"
00056 #include "TShapeSize2DB2.hpp"
00057 #include "TShapeSizeB3.hpp"
00058 #include "TShapeSizeOrientNB1.hpp"
00059 #include "TShapeSizeOrientB1.hpp"
00060 #include "TShapeSizeOrientB2.hpp"
00061 #include "InvTransBarrier.hpp"
00062 
00063 #include <iostream>
00064 #include <sstream>
00065 #include <stdlib.h>
00066 #include <assert.h>
00067 using namespace std;
00068 using namespace MBMesquite;
00069 
00070 #define CHKERR( A )                 \
00071     if( MSQ_CHKERR( ( A ) ) )       \
00072     {                               \
00073         cerr << ( A ) << std::endl; \
00074         exit( 3 );                  \
00075     }
00076 
00077 const double P            = 1.0;  // objective function power
00078 const bool ONE_TRI_SAMPLE = true;
00079 const bool LOCAL_PATCHES  = false;
00080 const bool USE_FEAS_NEWT  = true;
00081 const bool WRITE_VTK      = true;
00082 const bool WRITE_GNUPLOT  = false;
00083 
00084 void usage()
00085 {
00086     cerr << "main [e[.m]]" << endl;
00087     cerr << "Run experiments from 2d metrics paper" << endl;
00088     cerr << "e : expiriment" << endl;
00089     cerr << "m : metric number within experiment" << endl;
00090     cerr << "default is all" << endl;
00091 }
00092 
00093 typedef void ( *mesh_reader_t )( MeshImpl* mesh );
00094 
00095 bool run_smoother( mesh_reader_t input_mesh, mesh_reader_t reference_mesh, int exp, int n, TMetric* metric );
00096 
00097 void write_mesh( mesh_reader_t mesh, const char* filename );
00098 void write_mesh( MeshImpl* mesh, const char* filename );
00099 
00100 void quarter_annulus( double inner, double outer, Mesh* mesh );
00101 void parabolic_squash( double height, Mesh* mesh );
00102 void horseshoe( double x_inner, double x_outer, double y_inner, double y_outer, Mesh* mesh );
00103 void scale( double s, Mesh* mesh );
00104 
00105 void reference( MeshImpl* mesh )
00106 {
00107     MsqError err;
00108     std::string filename = std::string( STRINGIFY( SRCDIR ) ) + "/2d_metrics_reference.vtk";
00109     mesh->read_vtk( filename.c_str(), err );CHKERR( err )
00110 }
00111 
00112 void exp_1_init( MeshImpl* mesh )
00113 {
00114     MsqError err;
00115     std::string filename = std::string( STRINGIFY( SRCDIR ) ) + "/2d_metrics_reference.vtk";
00116     mesh->read_vtk( filename.c_str(), err );CHKERR( err )
00117 
00118     vector< Mesh::VertexHandle > handles;
00119     mesh->get_all_vertices( handles, err );CHKERR( err )
00120 
00121     // vector<Mesh::VertexHandle>::iterator i;
00122     for( size_t i = 0; i < 8; ++i )
00123     {
00124         MsqVertex vtx;
00125         mesh->vertices_get_coordinates( &handles[i], &vtx, 1, err );CHKERR( err )
00126 
00127         vtx[1] = -0.4 * ( vtx[0] - 0.5 ) * ( vtx[0] - 0.5 ) + 0.1;
00128 
00129         mesh->vertex_set_coordinates( handles[i], vtx, err );CHKERR( err )
00130     }
00131 }
00132 
00133 void exp_2_init( MeshImpl* mesh )
00134 {
00135     exp_1_init( mesh );
00136     scale( 8.0, mesh );
00137 }
00138 
00139 void exp_3_ref( MeshImpl* mesh )
00140 {
00141     reference( mesh );
00142     quarter_annulus( 0.65, 1.30, mesh );
00143 }
00144 
00145 void exp_4_init( MeshImpl* mesh )
00146 {
00147     reference( mesh );
00148     parabolic_squash( 0.8, mesh );
00149 }
00150 
00151 void exp_5_init( MeshImpl* mesh )
00152 {
00153     reference( mesh );
00154     horseshoe( 0.26, 0.53, 0.26, 1.3, mesh );
00155 }
00156 
00157 bool run_exp_1( int n );
00158 bool run_exp_2( int n );
00159 bool run_exp_3( int n );
00160 bool run_exp_4( int n );
00161 bool run_exp_5( int n );
00162 
00163 typedef bool ( *exp_func )( int );
00164 const exp_func experiment[] = { &run_exp_1, &run_exp_2, &run_exp_3, &run_exp_4, &run_exp_5 };
00165 
00166 TShapeSizeOrientNB1 nb1;
00167 TShapeOrientNB1 nb2;
00168 TShapeSize2DNB1 nb3;
00169 TShapeNB1 nb4;
00170 TSquared nb5;
00171 TShapeOrientNB2 nb6;
00172 TMetric* nb[] = { &nb1, &nb2, &nb3, &nb4, &nb5, &nb6 };
00173 
00174 TShapeSizeOrientB1 b1;
00175 TShapeOrientB1 b2;
00176 TShapeSize2DB2 b3;
00177 TShapeB1 b4;
00178 TShapeSizeOrientB2 b5;
00179 TShapeSizeB1 b6;
00180 TShapeSizeB3 b7;
00181 TMetric* b[] = { &b1, &b2, &b3, &b4, &b5, &b6, &b7 };
00182 
00183 bool run_exp_1( int n )
00184 {
00185     if( n > 6 ) return false;
00186 
00187     write_mesh( exp_1_init, "Exp1-initial" );
00188 
00189     bool r  = true, tmpr;
00190     int beg = 0, end = 5;
00191     if( n ) beg = end = n - 1;
00192 
00193     for( int i = beg; i <= end; ++i )
00194     {
00195         if( i == 5 )
00196             tmpr = run_smoother( &exp_1_init, 0, 1, 6, nb[0] );
00197         else
00198             tmpr = run_smoother( &exp_1_init, &reference, 1, i + 1, nb[i] );
00199         r = r && tmpr;
00200     }
00201     return r;
00202 }
00203 
00204 bool run_exp_2( int n )
00205 {
00206     if( n > 6 ) return false;
00207 
00208     write_mesh( exp_2_init, "Exp2-initial" );
00209 
00210     bool r  = true, tmpr;
00211     int beg = 0, end = 5;
00212     if( n ) beg = end = n - 1;
00213 
00214     for( int i = beg; i <= end; ++i )
00215     {
00216         if( i == 5 )
00217             tmpr = run_smoother( &exp_2_init, 0, 2, 6, nb[0] );
00218         else
00219             tmpr = run_smoother( &exp_2_init, &reference, 2, i + 1, nb[i] );
00220         r = r && tmpr;
00221     }
00222     return r;
00223 }
00224 
00225 bool run_exp_3( int n )
00226 {
00227     if( n > 6 ) return false;
00228 
00229     write_mesh( exp_3_ref, "Exp3-reference" );
00230 
00231     bool r  = true, tmpr;
00232     int beg = 0, end = 5;
00233     if( n ) beg = end = n - 1;
00234 
00235     for( int i = beg; i <= end; ++i )
00236     {
00237         if( i == 5 )
00238             tmpr = run_smoother( &exp_1_init, 0, 3, 6, nb[0] );
00239         else
00240             tmpr = run_smoother( &exp_1_init, &exp_3_ref, 3, i + 1, nb[i] );
00241         r = r && tmpr;
00242     }
00243     return r;
00244 }
00245 
00246 bool run_exp_4( int n )
00247 {
00248     if( n > 3 ) return false;
00249 
00250     write_mesh( exp_4_init, "Exp4-initial" );
00251     TMetric* m[3] = { &nb1, &b1, &b5 };
00252 
00253     bool r  = true, tmpr;
00254     int beg = 0, end = 2;
00255     if( n ) beg = end = n - 1;
00256 
00257     for( int i = beg; i <= end; ++i )
00258     {
00259         tmpr = run_smoother( &exp_4_init, &reference, 4, i + 1, m[i] );
00260         r    = r && tmpr;
00261     }
00262     return r;
00263 }
00264 
00265 bool run_exp_5( int n )
00266 {
00267     if( n > 12 || n == 5 ) return false;
00268 
00269     write_mesh( exp_5_init, "Exp5-initial" );
00270 
00271     bool r = true, tmpr = false;
00272     int beg = 0, end = 11;
00273     if( n ) beg = end = n - 1;
00274 
00275     for( int i = beg; i <= end; ++i )
00276     {
00277         switch( i )
00278         {
00279             case 0:
00280             case 1:
00281             case 2:
00282             case 3:
00283             case 4:
00284                 tmpr = run_smoother( &exp_5_init, &reference, 5, i + 1, b[i] );
00285                 break;
00286             case 5:
00287                 break;
00288             case 6:
00289             case 7:
00290             case 8:
00291             case 9:
00292             case 10: {
00293                 InvTransBarrier metric( nb[i - 6] );
00294                 tmpr = run_smoother( &exp_5_init, &reference, 5, i + 1, &metric );
00295             }
00296             break;
00297             case 11: {
00298                 InvTransBarrier metric( nb[0] );
00299                 tmpr = run_smoother( &exp_5_init, 0, 5, 12, &metric );
00300             }
00301             break;
00302         }
00303         r = r && tmpr;
00304     }
00305     return r;
00306 }
00307 
00308 void scale( double s, Mesh* mesh )
00309 {
00310     MsqError err;
00311     vector< Mesh::VertexHandle > handles;
00312     mesh->get_all_vertices( handles, err );CHKERR( err )
00313 
00314     vector< Mesh::VertexHandle >::iterator i;
00315     for( i = handles.begin(); i != handles.end(); ++i )
00316     {
00317         MsqVertex vtx;
00318         mesh->vertices_get_coordinates( &*i, &vtx, 1, err );CHKERR( err )
00319 
00320         vtx *= s;
00321 
00322         mesh->vertex_set_coordinates( *i, vtx, err );CHKERR( err )
00323     }
00324 }
00325 
00326 void quarter_annulus( double inner, double outer, Mesh* mesh )
00327 {
00328     MsqError err;
00329     vector< Mesh::VertexHandle > handles;
00330     mesh->get_all_vertices( handles, err );CHKERR( err )
00331 
00332     vector< Mesh::VertexHandle >::iterator i;
00333     for( i = handles.begin(); i != handles.end(); ++i )
00334     {
00335         MsqVertex vtx;
00336         mesh->vertices_get_coordinates( &*i, &vtx, 1, err );CHKERR( err )
00337 
00338         double r = inner + ( outer - inner ) * vtx[1];
00339         double a = M_PI * 0.5 * ( 1.0 - vtx[0] );
00340         vtx[0]   = r * cos( a );
00341         vtx[1]   = r * sin( a );
00342 
00343         mesh->vertex_set_coordinates( *i, vtx, err );CHKERR( err )
00344     }
00345 }
00346 
00347 void parabolic_squash( double height, Mesh* mesh )
00348 {
00349     MsqError err;
00350     vector< Mesh::VertexHandle > handles;
00351     mesh->get_all_vertices( handles, err );CHKERR( err )
00352 
00353     vector< Mesh::VertexHandle >::iterator i;
00354     for( i = handles.begin(); i != handles.end(); ++i )
00355     {
00356         MsqVertex vtx;
00357         mesh->vertices_get_coordinates( &*i, &vtx, 1, err );CHKERR( err )
00358 
00359         const double br = ( 1.0 - vtx[1] ) * height;
00360         const double a  = -4 * br;
00361         vtx[1] += a * ( vtx[0] - 0.5 ) * ( vtx[0] - 0.5 ) + br;
00362 
00363         mesh->vertex_set_coordinates( *i, vtx, err );CHKERR( err )
00364     }
00365 }
00366 
00367 void horseshoe( double x_inner, double x_outer, double y_inner, double y_outer, Mesh* mesh )
00368 {
00369     MsqError err;
00370     vector< Mesh::VertexHandle > handles;
00371     mesh->get_all_vertices( handles, err );CHKERR( err )
00372 
00373     vector< Mesh::VertexHandle >::iterator i;
00374     for( i = handles.begin(); i != handles.end(); ++i )
00375     {
00376         MsqVertex vtx;
00377         mesh->vertices_get_coordinates( &*i, &vtx, 1, err );CHKERR( err )
00378 
00379         double a = M_PI * ( 1 - vtx[0] );
00380         vtx[0]   = ( x_inner + ( x_outer - x_inner ) * vtx[1] ) * cos( a );
00381         vtx[1]   = ( y_inner + ( y_outer - y_inner ) * vtx[1] ) * sin( a );
00382 
00383         mesh->vertex_set_coordinates( *i, vtx, err );CHKERR( err )
00384     }
00385 }
00386 
00387 void write_mesh( MeshImpl* mesh, const char* filename )
00388 {
00389     MsqError err;
00390     if( WRITE_VTK )
00391     {
00392         string vfile = string( filename ) + ".vtk";
00393         mesh->write_vtk( vfile.c_str(), err );CHKERR( err )
00394         cout << "Wrote: \"" << vfile << '"' << endl;
00395     }
00396     if( WRITE_GNUPLOT )
00397     {
00398         string vfile = string( filename ) + ".eps";
00399         MeshWriter::write_eps( mesh, vfile.c_str(), MeshWriter::Projection( MeshWriter::X, MeshWriter::Y ), err );CHKERR( err )
00400         cout << "Wrote: \"" << vfile << '"' << endl;
00401     }
00402 }
00403 
00404 void write_mesh( mesh_reader_t mesh_func, const char* filename )
00405 {
00406     MeshImpl mesh;
00407     mesh_func( &mesh );
00408     write_mesh( &mesh, filename );
00409 }
00410 
00411 bool run_smoother( mesh_reader_t input_mesh, mesh_reader_t reference_mesh, int exp, int n, TMetric* target_metric )
00412 {
00413     MsqError err;
00414     MeshImpl active, reference;
00415     input_mesh( &active );
00416 
00417     ReferenceMesh refmesh( &reference );
00418     RefMeshTargetCalculator ref_target( &refmesh );
00419     IdealShapeTarget ident_target;
00420     TargetCalculator* target;
00421     if( reference_mesh )
00422     {
00423         reference_mesh( &reference );
00424         target = &ref_target;
00425     }
00426     else
00427     {
00428         target = &ident_target;
00429     }
00430 
00431     TQualityMetric metric( target, target_metric );
00432 
00433     TerminationCriterion outer, inner;
00434     if( LOCAL_PATCHES )
00435     {
00436         outer.add_iteration_limit( 100 );
00437         outer.add_absolute_vertex_movement( 1e-3 );
00438         inner.add_iteration_limit( 3 );
00439     }
00440     else
00441     {
00442         outer.add_iteration_limit( 1 );
00443         inner.add_absolute_vertex_movement( 1e-4 );
00444         inner.add_iteration_limit( 100 );
00445     }
00446 
00447     PMeanPTemplate of( P, &metric );
00448     ConjugateGradient cg( &of );
00449     FeasibleNewton fn( &of );
00450     if( LOCAL_PATCHES )
00451     {
00452         cg.use_element_on_vertex_patch();
00453         fn.use_element_on_vertex_patch();
00454     }
00455     else
00456     {
00457         cg.use_global_patch();
00458         fn.use_global_patch();
00459     }
00460     VertexMover* solver = USE_FEAS_NEWT ? (VertexMover*)&fn : (VertexMover*)&cg;
00461     XYPlanarDomain plane2;
00462     solver->set_inner_termination_criterion( &inner );
00463     solver->set_outer_termination_criterion( &outer );
00464 
00465     InstructionQueue q;
00466     q.set_master_quality_improver( solver, err );CHKERR( err )
00467 
00468     cout << "Running " << exp << "." << n << " ...";
00469 
00470     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &active, &plane2 );
00471     q.run_instructions( &mesh_and_domain, err );
00472     if( MSQ_CHKERR( err ) )
00473     {
00474         cout << "######## EXPERIMENT " << exp << "." << n << " FAILED! ##########" << endl;
00475         return false;
00476     }
00477 
00478     ostringstream s;
00479     s << "Exp" << exp << "-" << n;
00480     write_mesh( &active, s.str().c_str() );
00481     return true;
00482 }
00483 
00484 int main( int argc, char* argv[] )
00485 {
00486     if( argc > 2 ) usage();
00487     int exp = 0, n = 0;
00488     if( argc == 2 && !sscanf( argv[1], "%d.%d", &exp, &n ) ) usage();
00489 
00490     if( exp > 5 ) usage();
00491 
00492     int lower, upper;
00493     if( exp > 0 ) { lower = upper = exp - 1; }
00494     else
00495     {
00496         lower = 0;
00497         upper = 4;
00498     }
00499 
00500     if( WRITE_GNUPLOT ) write_mesh( reference, "reference" );
00501 
00502     int fail = 0;
00503     for( int e = lower; e <= upper; ++e )
00504         fail += !( experiment[e]( n ) );
00505 
00506     return fail;
00507 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines