MOAB: Mesh Oriented datABase  (version 5.3.1)
idft_time_test.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in 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     kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 #include <iostream>
00027 #include <sstream>
00028 using std::cerr;
00029 using std::cout;
00030 using std::endl;
00031 
00032 #include <cstdlib>
00033 
00034 #include "Mesquite.hpp"
00035 #include "MeshImpl.hpp"
00036 #include "MsqError.hpp"
00037 #include "InstructionQueue.hpp"
00038 #include "TerminationCriterion.hpp"
00039 #include "QualityAssessor.hpp"
00040 #include "PlanarDomain.hpp"
00041 
00042 #include "IdealWeightInverseMeanRatio.hpp"
00043 #include "ElementPMeanP.hpp"
00044 #include "TInverseMeanRatio.hpp"
00045 #include "TQualityMetric.hpp"
00046 #include "TOffset.hpp"
00047 #include "TOffset.hpp"
00048 
00049 #include "IdealShapeTarget.hpp"
00050 #include "PMeanPTemplate.hpp"
00051 #include "SteepestDescent.hpp"
00052 #include "ConjugateGradient.hpp"
00053 #include "QuasiNewton.hpp"
00054 
00055 #include "MeshWriter.hpp"
00056 #include "MsqTimer.hpp"
00057 #include "TestUtil.hpp"
00058 
00059 using namespace MBMesquite;
00060 
00061 std::string DEFAULT_INPUT = TestDir + "unittest/mesquite/2D/vtk/tris/untangled/equil_tri2.vtk";
00062 
00063 /* Print usage or help information: exits if err == true */
00064 void usage( const char* argv0 = 0, bool err = true )
00065 {
00066     std::ostream& s     = err ? cerr : cout;
00067     const char* defname = "main";
00068     if( !argv0 ) argv0 = defname;
00069 
00070     s << "Usage: " << defname << " [-n|-c|-q] [-e] [-t] [-a] [-N] [{-v|-g|-p} <output_file>] [<input_file>]" << endl
00071       << "       " << defname << " -h" << std::endl;
00072 
00073     if( err ) exit( 1 );
00074 
00075     s << "  -n : Use SteepestDescent solver (default)" << endl
00076       << "  -c : Use ConjugateGradient solver" << endl
00077       << "  -q : Use QuasiNewton solver" << endl
00078       << "  -e : Test IdealWeightInverseMeanRatio metric" << endl
00079       << "  -t : Test InverseMeanRatio3D target metric" << endl
00080       << "  -a : Test ElementPMeanP(InverseMeanRatio3D)" << endl
00081       << "  -N : Test InverseMeanRatio3D with finite difference derivatives" << endl
00082       << "  -C : Compare InverseMeanRatio3D and IdealWeightInverseMeanRatio" << endl
00083       << "  -v : Write output mesh to specified VTK file" << endl
00084       << "  -g : Write output mesh to specified GnuPlot file" << endl
00085       << "  -p : Write output mesh to specified EPS file" << endl
00086       << "  -P : Write solver plot data to specified file." << endl
00087       << "Default is all metrics if non specified." << endl
00088       << "Default input file: " << DEFAULT_INPUT << endl;
00089 }
00090 
00091 const char* vtk_file  = 0; /* vtk output file name */
00092 const char* eps_file  = 0; /* eps output file name */
00093 const char* gpt_file  = 0; /* GNUPlot output file name */
00094 const char* plot_file = 0; /* Time-dependent plot of solver data */
00095 
00096 enum Solver
00097 {
00098     STEEP_DESCENT,
00099     CONJ_GRAD,
00100     QUASI_NEWT
00101 };
00102 
00103 /* Run an optimization: returns average quality of final mesh */
00104 double run( QualityMetric* metric, Solver solver, const char* input_file, double& seconds_out, int& iterations_out );
00105 
00106 /* Force finite difference approximation of derivatives for
00107  * an existing quality metric */
00108 class NumericQM : public QualityMetric
00109 {
00110   private:
00111     QualityMetric* realMetric;
00112 
00113   public:
00114     NumericQM( QualityMetric* real_metric ) : realMetric( real_metric ) {}
00115     virtual MetricType get_metric_type() const;
00116     virtual std::string get_name() const;
00117     virtual int get_negate_flag() const;
00118     virtual void get_evaluations( PatchData& pd,
00119                                   std::vector< size_t >& handles,
00120                                   bool free_vertices_only,
00121                                   MsqError& err );
00122     virtual bool evaluate( PatchData& pd, size_t handle, double& value, MsqError& err );
00123     virtual bool evaluate_with_indices( PatchData& pd,
00124                                         size_t handle,
00125                                         double& value,
00126                                         std::vector< size_t >& indices,
00127                                         MsqError& err );
00128 };
00129 QualityMetric::MetricType NumericQM::get_metric_type() const
00130 {
00131     return realMetric->get_metric_type();
00132 }
00133 std::string NumericQM::get_name() const
00134 {
00135     std::string r = realMetric->get_name();
00136     r += " (FD)";
00137     return r;
00138 }
00139 
00140 /* At each evaluation of this metric, compare the values resulting
00141  * from the evaluation of two other metrics: flag an error if they
00142  * differ or return the common result if the same */
00143 class CompareMetric : public QualityMetric
00144 {
00145   private:
00146     QualityMetric *metric1, *metric2;
00147     bool maskPlane;
00148     int maskAxis;
00149     std::vector< size_t > m2Handles;
00150     std::vector< Vector3D > m2Grad;
00151     std::vector< SymMatrix3D > m2Diag;
00152     std::vector< Matrix3D > m2Hess;
00153     static const double epsilon;
00154 
00155   public:
00156     CompareMetric( QualityMetric* qm1, QualityMetric* qm2, bool mask_qm2_coord = false )
00157         : metric1( qm1 ), metric2( qm2 ), maskPlane( mask_qm2_coord ), maskAxis( -1 )
00158     {
00159     }
00160 
00161     MetricType get_metric_type() const;
00162 
00163     std::string get_name() const;
00164 
00165     int get_negate_flag() const;
00166 
00167     void get_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free_vertices_only, MsqError& err );
00168 
00169     bool evaluate( PatchData& pd, size_t handle, double& value, MsqError& err );
00170 
00171     bool evaluate_with_indices( PatchData& pd,
00172                                 size_t handle,
00173                                 double& value,
00174                                 std::vector< size_t >& indices,
00175                                 MsqError& err );
00176 
00177     bool evaluate_with_gradient( PatchData& pd,
00178                                  size_t handle,
00179                                  double& value,
00180                                  std::vector< size_t >& indices,
00181                                  std::vector< Vector3D >& gradient,
00182                                  MsqError& err );
00183 
00184     bool evaluate_with_Hessian_diagonal( PatchData& pd,
00185                                          size_t handle,
00186                                          double& value,
00187                                          std::vector< size_t >& indices,
00188                                          std::vector< Vector3D >& gradient,
00189                                          std::vector< SymMatrix3D >& Hessian_diagonal,
00190                                          MsqError& err );
00191 
00192     bool evaluate_with_Hessian( PatchData& pd,
00193                                 size_t handle,
00194                                 double& value,
00195                                 std::vector< size_t >& indices,
00196                                 std::vector< Vector3D >& gradient,
00197                                 std::vector< Matrix3D >& Hessian,
00198                                 MsqError& err );
00199 
00200     void get_mask_axis( PatchData& pd );
00201     bool equal( Vector3D grad1, const Vector3D& grad2 ) const;
00202     bool equal( SymMatrix3D hess1, const SymMatrix3D& hess2 ) const;
00203     bool equal( Matrix3D hess1, const Matrix3D& hess2 ) const;
00204 };
00205 const double CompareMetric::epsilon = 5e-2;
00206 
00207 /* Parse command line options and call 'run' */
00208 int main( int argc, char* argv[] )
00209 {
00210     Solver solver              = STEEP_DESCENT;
00211     bool do_non_target_metric  = false;
00212     bool do_new_target_metric  = false;
00213     bool do_new_target_average = false;
00214     bool do_new_target_numeric = false;
00215     bool do_compare_metric     = false;
00216     const char* input_file     = DEFAULT_INPUT.c_str();
00217     bool no_more_flags         = false;
00218 
00219     std::list< const char** > exp_list;
00220     for( int i = 1; i < argc; ++i )
00221     {
00222         if( argv[i][0] == '-' && !no_more_flags )
00223         {
00224             for( int k = 1; argv[i][k]; ++k )
00225             {
00226                 switch( argv[i][k] )
00227                 {
00228                     case 'n':
00229                         solver = STEEP_DESCENT;
00230                         break;
00231                     case 'c':
00232                         solver = CONJ_GRAD;
00233                         break;
00234                     case 'q':
00235                         solver = QUASI_NEWT;
00236                         break;
00237                     case 'e':
00238                         do_non_target_metric = true;
00239                         break;
00240                     case 't':
00241                         do_new_target_metric = true;
00242                         break;
00243                     case 'a':
00244                         do_new_target_average = true;
00245                         break;
00246                     case 'N':
00247                         do_new_target_numeric = true;
00248                         break;
00249                     case 'C':
00250                         do_compare_metric = true;
00251                         break;
00252                     case 'v':
00253                         exp_list.push_back( &vtk_file );
00254                         break;
00255                     case 'g':
00256                         exp_list.push_back( &gpt_file );
00257                         break;
00258                     case 'p':
00259                         exp_list.push_back( &eps_file );
00260                         break;
00261                     case 'P':
00262                         exp_list.push_back( &plot_file );
00263                         break;
00264                     case 'h':
00265                         usage( argv[0], false );
00266                         return 0;
00267                     case '-':
00268                         no_more_flags = true;
00269                         break;
00270                     default:
00271                         cerr << "Invalid flag: '" << argv[i][k] << "'" << endl;
00272                         usage( argv[0] );
00273                 }
00274             }
00275         }
00276         else if( !exp_list.empty() )
00277         {
00278             const char** ptr = exp_list.front();
00279             exp_list.pop_front();
00280             *ptr = argv[i];
00281         }
00282         else
00283         {
00284             if( !DEFAULT_INPUT.compare( input_file ) )
00285             {
00286                 cerr << "Unexpected argument: \"" << argv[i] << '"' << endl;
00287                 usage( argv[0] );
00288             }
00289             input_file = argv[i];
00290         }
00291     }
00292 
00293     int count = 0;
00294     if( do_non_target_metric ) ++count;
00295     if( do_new_target_metric ) ++count;
00296     if( do_new_target_average ) ++count;
00297     if( do_new_target_numeric ) ++count;
00298     if( do_compare_metric ) ++count;
00299 
00300     if( !count )
00301     {
00302         do_compare_metric = true;
00303         count             = 1;
00304     }
00305 
00306     if( ( vtk_file || gpt_file || eps_file ) && count != 1 )
00307     {
00308         cerr << "Error: Cannot write output file if running multiple tests" << endl;
00309         return 2;
00310     }
00311 
00312     IdealWeightInverseMeanRatio non_target_metric;
00313     IdealShapeTarget new_target;
00314     TInverseMeanRatio tmp;
00315     TOffset tmp_off( 1.0, &tmp );  // target metrics are zero-based
00316     TQualityMetric new_target_metric( &new_target, &tmp_off );
00317     ElementPMeanP new_target_average( 1.0, &new_target_metric );
00318     NumericQM new_target_numeric( &new_target_metric );
00319     CompareMetric comp_metric( &non_target_metric, &new_target_average, true );
00320 
00321     std::ostringstream os;
00322     double secs, qual;
00323     if( do_non_target_metric )
00324     {
00325         qual = run( &non_target_metric, solver, input_file, secs, count );
00326         os << "IdealWeightInverseMeanRatio: " << qual << " after " << count << " iterations in " << secs << " seconds"
00327            << endl;
00328     }
00329     if( do_new_target_metric )
00330     {
00331         qual = run( &new_target_metric, solver, input_file, secs, count );
00332         os << "TQualityMetric              : " << qual << " after " << count << " iterations in " << secs << " seconds"
00333            << endl;
00334     }
00335     if( do_new_target_average )
00336     {
00337         qual = run( &new_target_average, solver, input_file, secs, count );
00338         os << "ElementPMeanP              : " << qual << " after " << count << " iterations in " << secs << " seconds"
00339            << endl;
00340     }
00341     if( do_new_target_numeric )
00342     {
00343         qual = run( &new_target_numeric, solver, input_file, secs, count );
00344         os << "TQualityMetric (FD)         : " << qual << " after " << count << " iterations in " << secs << " seconds"
00345            << endl;
00346     }
00347     if( do_compare_metric )
00348     {
00349         qual = run( &comp_metric, solver, input_file, secs, count );
00350         os << "Metric comparison      : " << qual << " after " << count << " iterations in " << secs << " seconds"
00351            << endl;
00352     }
00353 
00354     cout << endl << os.str() << endl;
00355     return 0;
00356 }
00357 
00358 double run( QualityMetric* metric,
00359             Solver solver_type,
00360             const char* input_file,
00361             double& seconds_out,
00362             int& iterations_out )
00363 {
00364     MsqPrintError err( cerr );
00365     IdealWeightInverseMeanRatio qa_metric;
00366     TerminationCriterion inner, outer;
00367     outer.add_iteration_limit( 1 );
00368     inner.add_absolute_vertex_movement( 1e-4 );
00369     inner.add_iteration_limit( 100 );
00370     PMeanPTemplate of( 1.0, metric );
00371     QualityAssessor qa( &qa_metric );
00372     qa.add_quality_assessment( metric );
00373     InstructionQueue q;
00374     SteepestDescent steep( &of );
00375     QuasiNewton quasi( &of );
00376     ConjugateGradient conj( &of );
00377     VertexMover* solver = 0;
00378     switch( solver_type )
00379     {
00380         case STEEP_DESCENT:
00381             solver = &steep;
00382             break;
00383         case QUASI_NEWT:
00384             solver = &quasi;
00385             break;
00386         case CONJ_GRAD:
00387             solver = &conj;
00388             break;
00389     }
00390     q.set_master_quality_improver( solver, err );
00391     q.add_quality_assessor( &qa, err );
00392     solver->set_inner_termination_criterion( &inner );
00393     solver->set_outer_termination_criterion( &outer );
00394 
00395     if( plot_file ) inner.write_iterations( plot_file, err );
00396 
00397     MeshImpl mesh;
00398     mesh.read_vtk( input_file, err );
00399     if( err )
00400     {
00401         cerr << "Failed to read input file: \"" << input_file << '"' << endl;
00402         exit( 1 );
00403     }
00404 
00405     std::vector< Mesh::VertexHandle > handles;
00406     mesh.get_all_vertices( handles, err );
00407     if( handles.empty() )
00408     {
00409         cerr << "no veritces in mesh" << endl;
00410         exit( 1 );
00411     }
00412     std::vector< MsqVertex > coords( handles.size() );
00413     mesh.vertices_get_coordinates( arrptr( handles ), arrptr( coords ), handles.size(), err );
00414     Vector3D min( HUGE_VAL ), max( -HUGE_VAL );
00415     for( size_t i = 0; i < coords.size(); ++i )
00416     {
00417         for( int j = 0; j < 3; ++j )
00418         {
00419             if( coords[i][j] < min[j] ) min[j] = coords[i][j];
00420             if( coords[i][j] > max[j] ) max[j] = coords[i][j];
00421         }
00422     }
00423 
00424     Vector3D size        = max - min;
00425     PlanarDomain* domain = 0;
00426     if( size[0] < 1e-4 )
00427         domain = new PlanarDomain( PlanarDomain::YZ, min[0] );
00428     else if( size[1] < 1e-4 )
00429         domain = new PlanarDomain( PlanarDomain::XZ, min[1] );
00430     else if( size[2] < 1e-4 )
00431         domain = new PlanarDomain( PlanarDomain::XY, min[2] );
00432 
00433     Timer timer;
00434     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, domain );
00435     q.run_instructions( &mesh_and_domain, err );
00436     seconds_out = timer.since_birth();
00437     if( err )
00438     {
00439         cerr << "Optimization failed." << endl << err << endl;
00440         abort();
00441     }
00442 
00443     if( vtk_file )
00444     {
00445         MeshWriter::write_vtk( &mesh, vtk_file, err );
00446         if( err ) cerr << vtk_file << ": failed to write file." << endl;
00447     }
00448     if( gpt_file )
00449     {
00450         MeshWriter::write_gnuplot( &mesh, gpt_file, err );
00451         if( err ) cerr << gpt_file << ": failed to write file." << endl;
00452     }
00453     if( eps_file )
00454     {
00455         PlanarDomain xy( PlanarDomain::XY );
00456         MeshWriter::Projection proj( domain ? domain : &xy );
00457         MeshWriter::write_eps( &mesh, eps_file, proj, err );
00458         if( err ) cerr << eps_file << ": failed to write file." << endl;
00459     }
00460     delete domain;
00461 
00462     iterations_out = inner.get_iteration_count();
00463 
00464     const QualityAssessor::Assessor* a = qa.get_results( &qa_metric );
00465     return a->get_average();
00466 }
00467 
00468 int NumericQM::get_negate_flag() const
00469 {
00470     return realMetric->get_negate_flag();
00471 }
00472 
00473 void NumericQM::get_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free_vertices_only, MsqError& err )
00474 {
00475     realMetric->get_evaluations( pd, handles, free_vertices_only, err );
00476 }
00477 
00478 bool NumericQM::evaluate( PatchData& pd, size_t handle, double& value, MsqError& err )
00479 {
00480     return realMetric->evaluate( pd, handle, value, err );
00481 }
00482 
00483 bool NumericQM::evaluate_with_indices( PatchData& pd,
00484                                        size_t handle,
00485                                        double& value,
00486                                        std::vector< size_t >& indices,
00487                                        MsqError& err )
00488 {
00489     return realMetric->evaluate_with_indices( pd, handle, value, indices, err );
00490 }
00491 
00492 QualityMetric::MetricType CompareMetric::get_metric_type() const
00493 {
00494     MetricType t1 = metric1->get_metric_type();
00495     assert( metric2->get_metric_type() == t1 );
00496     return t1;
00497 }
00498 
00499 std::string CompareMetric::get_name() const
00500 {
00501     std::string n = metric1->get_name();
00502     n += " =? ";
00503     n += metric2->get_name();
00504     return n;
00505 }
00506 
00507 int CompareMetric::get_negate_flag() const
00508 {
00509     assert( metric1->get_negate_flag() == metric2->get_negate_flag() );
00510     return metric1->get_negate_flag();
00511 }
00512 
00513 void CompareMetric::get_evaluations( PatchData& pd,
00514                                      std::vector< size_t >& handles,
00515                                      bool free_vertices_only,
00516                                      MsqError& err )
00517 {
00518     if( maskPlane ) get_mask_axis( pd );
00519 
00520     m2Handles.clear();
00521     metric1->get_evaluations( pd, handles, free_vertices_only, err );MSQ_ERRRTN( err );
00522     metric2->get_evaluations( pd, m2Handles, free_vertices_only, err );MSQ_ERRRTN( err );
00523     bool same = ( handles.size() == m2Handles.size() );
00524     std::sort( m2Handles.begin(), m2Handles.end() );
00525     for( std::vector< size_t >::iterator i = handles.begin(); i != handles.end(); ++i )
00526         if( !std::binary_search( m2Handles.begin(), m2Handles.end(), *i ) ) same = false;
00527     if( !same )
00528     {
00529         MSQ_SETERR( err )
00530         ( "Metrics have incompatible lists of evaluation handles.\n", MsqError::INVALID_STATE );
00531     }
00532 }
00533 
00534 bool CompareMetric::evaluate( PatchData& pd, size_t handle, double& value, MsqError& err )
00535 {
00536     double m2val;
00537     bool r1, r2;
00538     r1 = metric1->evaluate( pd, handle, value, err );
00539     MSQ_ERRZERO( err );
00540     r2 = metric2->evaluate( pd, handle, m2val, err );
00541     MSQ_ERRZERO( err );
00542     if( r1 != r2 || ( r1 && fabs( value - m2val ) > epsilon ) )
00543     {
00544         MSQ_SETERR( err )
00545         ( MsqError::INVALID_STATE,
00546           "Metrics returned different values for handle %lu "
00547           "in evaluate:\n"
00548           "\t%s %f vs. %s %f\n",
00549           (unsigned long)handle, r1 ? "true" : "false", value, r2 ? "true" : "false", m2val );
00550     }
00551 
00552     return r1 && !err;
00553 }
00554 
00555 bool CompareMetric::evaluate_with_indices( PatchData& pd,
00556                                            size_t handle,
00557                                            double& value,
00558                                            std::vector< size_t >& indices,
00559                                            MsqError& err )
00560 {
00561     double m2val;
00562     bool r1, r2;
00563     m2Handles.clear();
00564     r1 = metric1->evaluate_with_indices( pd, handle, value, indices, err );
00565     MSQ_ERRZERO( err );
00566     r2 = metric2->evaluate_with_indices( pd, handle, m2val, m2Handles, err );
00567     MSQ_ERRZERO( err );
00568     if( r1 != r2 || ( r1 && fabs( value - m2val ) > epsilon ) )
00569     {
00570         MSQ_SETERR( err )
00571         ( MsqError::INVALID_STATE,
00572           "Metrics returned different values for handle %lu "
00573           "in evaluate_with_indices:\n"
00574           "\t%s %f vs. %s %f\n",
00575           (unsigned long)handle, r1 ? "true" : "false", value, r2 ? "true" : "false", m2val );
00576     }
00577     else
00578     {
00579         bool same = ( indices.size() == m2Handles.size() );
00580         std::sort( m2Handles.begin(), m2Handles.end() );
00581         for( std::vector< size_t >::iterator i = indices.begin(); i != indices.end(); ++i )
00582             if( !std::binary_search( m2Handles.begin(), m2Handles.end(), *i ) ) same = false;
00583         if( !same )
00584         {
00585             MSQ_SETERR( err )
00586             ( MsqError::INVALID_STATE,
00587               "Metrics returned incompatible lists of vertex indices"
00588               " for handle %lu in evaluate_with_indices\n.",
00589               (unsigned long)handle );
00590         }
00591     }
00592 
00593     return r1 && !err;
00594 }
00595 
00596 bool CompareMetric::evaluate_with_gradient( PatchData& pd,
00597                                             size_t handle,
00598                                             double& value,
00599                                             std::vector< size_t >& indices,
00600                                             std::vector< Vector3D >& gradient,
00601                                             MsqError& err )
00602 {
00603     double m2val;
00604     bool r1, r2;
00605     m2Handles.clear();
00606     m2Grad.clear();
00607     r1 = metric1->evaluate_with_gradient( pd, handle, value, indices, gradient, err );
00608     MSQ_ERRZERO( err );
00609     r2 = metric2->evaluate_with_gradient( pd, handle, m2val, m2Handles, m2Grad, err );
00610     MSQ_ERRZERO( err );
00611     if( r1 != r2 || ( r1 && fabs( value - m2val ) > epsilon ) )
00612     {
00613         MSQ_SETERR( err )
00614         ( MsqError::INVALID_STATE,
00615           "Metrics returned different values for handle %lu in "
00616           "evaluate_with_gradient:\n"
00617           "\t%s %f vs. %s %f\n",
00618           (unsigned long)handle, r1 ? "true" : "false", value, r2 ? "true" : "false", m2val );
00619     }
00620     else
00621     {
00622         std::vector< size_t >::const_iterator i, j;
00623         std::vector< Vector3D >::const_iterator r, s;
00624         int grad_diff = 0;
00625         bool same     = ( indices.size() == m2Handles.size() );
00626         std::sort( m2Handles.begin(), m2Handles.end() );
00627         for( i = indices.begin(); i != indices.end(); ++i )
00628         {
00629             j = std::lower_bound( m2Handles.begin(), m2Handles.end(), *i );
00630             if( j == m2Handles.end() || *j != *i )
00631             {
00632                 same = false;
00633                 continue;
00634             }
00635 
00636             r = gradient.begin() + ( i - indices.begin() );
00637             s = m2Grad.begin() + ( j - m2Handles.begin() );
00638             if( !equal( *r, *s ) ) ++grad_diff;
00639         }
00640 
00641         if( !same )
00642         {
00643             MSQ_SETERR( err )
00644             ( MsqError::INVALID_STATE,
00645               "Metrics returned incompatible lists of vertex indices"
00646               " for handle %lu in evaluate_with_gradient\n.",
00647               (unsigned long)handle );
00648         }
00649         else if( grad_diff )
00650         {
00651             MSQ_SETERR( err )
00652             ( MsqError::INVALID_STATE,
00653               "Metrics returned different gradient vectors for "
00654               " %d of %u vertices for handle %lu in "
00655               "evaluate_with_gradient\n.",
00656               grad_diff, (unsigned)gradient.size(), (unsigned long)handle );
00657         }
00658     }
00659 
00660     return r1 && !err;
00661 }
00662 
00663 bool CompareMetric::evaluate_with_Hessian_diagonal( PatchData& pd,
00664                                                     size_t handle,
00665                                                     double& value,
00666                                                     std::vector< size_t >& indices,
00667                                                     std::vector< Vector3D >& gradient,
00668                                                     std::vector< SymMatrix3D >& diagonal,
00669                                                     MsqError& err )
00670 {
00671     double m2val;
00672     bool r1, r2;
00673     m2Handles.clear();
00674     m2Grad.clear();
00675     m2Diag.clear();
00676     r1 = metric1->evaluate_with_Hessian_diagonal( pd, handle, value, indices, gradient, diagonal, err );
00677     MSQ_ERRZERO( err );
00678     r2 = metric2->evaluate_with_Hessian_diagonal( pd, handle, m2val, m2Handles, m2Grad, m2Diag, err );
00679     MSQ_ERRZERO( err );
00680     if( r1 != r2 || ( r1 && fabs( value - m2val ) > epsilon ) )
00681     {
00682         MSQ_SETERR( err )
00683         ( MsqError::INVALID_STATE,
00684           "Metrics returned different values for handle %lu in "
00685           "evaluate_with_Hessian_diagonal:\n"
00686           "\t%s %f vs. %s %f\n",
00687           (unsigned long)handle, r1 ? "true" : "false", value, r2 ? "true" : "false", m2val );
00688     }
00689     else
00690     {
00691         std::vector< size_t >::const_iterator i, j;
00692         std::vector< Vector3D >::const_iterator r, s;
00693         std::vector< SymMatrix3D >::const_iterator u, v;
00694         int grad_diff = 0, hess_diff = 0;
00695         bool same = ( indices.size() == m2Handles.size() );
00696         std::sort( m2Handles.begin(), m2Handles.end() );
00697         for( i = indices.begin(); i != indices.end(); ++i )
00698         {
00699             j = std::lower_bound( m2Handles.begin(), m2Handles.end(), *i );
00700             if( j == m2Handles.end() || *j != *i )
00701             {
00702                 same = false;
00703                 continue;
00704             }
00705 
00706             r = gradient.begin() + ( i - indices.begin() );
00707             s = m2Grad.begin() + ( j - m2Handles.begin() );
00708             if( !equal( *r, *s ) ) ++grad_diff;
00709 
00710             u = diagonal.begin() + ( i - indices.begin() );
00711             v = m2Diag.begin() + ( j - m2Handles.begin() );
00712             if( !equal( *u, *v ) ) ++hess_diff;
00713         }
00714 
00715         if( !same )
00716         {
00717             MSQ_SETERR( err )
00718             ( MsqError::INVALID_STATE,
00719               "Metrics returned incompatible lists of vertex indices"
00720               " for handle %lu in evaluate_with_Hessian_diagonal\n.",
00721               (unsigned long)handle );
00722         }
00723         else if( grad_diff )
00724         {
00725             MSQ_SETERR( err )
00726             ( MsqError::INVALID_STATE,
00727               "Metrics returned different gradient vectors for "
00728               " %d of %u vertices for handle %lu in "
00729               "evaluate_with_Hessian_diagonal\n.",
00730               grad_diff, (unsigned)gradient.size(), (unsigned long)handle );
00731         }
00732         else if( hess_diff )
00733         {
00734             MSQ_SETERR( err )
00735             ( MsqError::INVALID_STATE,
00736               "Metrics returned different Hessian blocks for "
00737               " %d of %u vertices for handle %lu in "
00738               "evaluate_with_Hessian_diagonal\n.",
00739               hess_diff, (unsigned)diagonal.size(), (unsigned long)handle );
00740         }
00741     }
00742 
00743     return r1 && !err;
00744 }
00745 
00746 bool CompareMetric::evaluate_with_Hessian( PatchData& pd,
00747                                            size_t handle,
00748                                            double& value,
00749                                            std::vector< size_t >& indices,
00750                                            std::vector< Vector3D >& gradient,
00751                                            std::vector< Matrix3D >& Hessian,
00752                                            MsqError& err )
00753 {
00754     double m2val;
00755     bool r1, r2;
00756     m2Handles.clear();
00757     m2Grad.clear();
00758     m2Hess.clear();
00759     r1 = metric1->evaluate_with_Hessian( pd, handle, value, indices, gradient, Hessian, err );
00760     MSQ_ERRZERO( err );
00761     r2 = metric2->evaluate_with_Hessian( pd, handle, m2val, m2Handles, m2Grad, m2Hess, err );
00762     MSQ_ERRZERO( err );
00763     if( r1 != r2 || fabs( value - m2val ) > epsilon )
00764     {
00765         MSQ_SETERR( err )
00766         ( MsqError::INVALID_STATE,
00767           "Metrics returned different values for handle %lu in "
00768           "evaluate_with_Hessian:\n"
00769           "\t%s %f vs. %s %f\n",
00770           (unsigned long)handle, r1 ? "true" : "false", value, r2 ? "true" : "false", m2val );
00771     }
00772     else
00773     {
00774         std::vector< size_t >::const_iterator i, j;
00775         std::vector< Vector3D >::const_iterator r, s;
00776         int grad_diff = 0, hess_diff = 0;
00777         bool same = ( indices.size() == m2Handles.size() );
00778         std::sort( m2Handles.begin(), m2Handles.end() );
00779         for( i = indices.begin(); i != indices.end(); ++i )
00780         {
00781             j = std::lower_bound( m2Handles.begin(), m2Handles.end(), *i );
00782             if( j == m2Handles.end() || *j != *i )
00783             {
00784                 same = false;
00785                 continue;
00786             }
00787 
00788             r = gradient.begin() + ( i - indices.begin() );
00789             s = m2Grad.begin() + ( j - m2Handles.begin() );
00790             if( !equal( *r, *s ) )
00791             {
00792                 ++grad_diff;
00793                 // call again for so debugger can step into it after failure is found
00794                 std::vector< size_t > i2;
00795                 std::vector< Vector3D > g2;
00796                 std::vector< Matrix3D > h2;
00797                 metric2->evaluate_with_Hessian( pd, handle, m2val, i2, g2, h2, err );
00798             }
00799         }
00800 
00801         if( !same )
00802         {
00803             MSQ_SETERR( err )
00804             ( MsqError::INVALID_STATE,
00805               "Metrics returned incompatible lists of vertex indices"
00806               " for handle %lu in evaluate_with_Hessian\n.",
00807               (unsigned long)handle );
00808         }
00809         else if( grad_diff )
00810         {
00811             MSQ_SETERR( err )
00812             ( MsqError::INVALID_STATE,
00813               "Metrics returned different gradient vectors for "
00814               " %d of %u vertices for handle %lu in "
00815               "evaluate_with_Hessian\n.",
00816               grad_diff, (unsigned)gradient.size(), (unsigned long)handle );
00817         }
00818         else
00819         {
00820             size_t row, col, row2, col2, idx, idx2;
00821             for( row = idx = 0; row < indices.size(); ++row )
00822             {
00823                 row2 = std::lower_bound( m2Handles.begin(), m2Handles.end(), indices[row] ) - m2Handles.begin();
00824                 for( col = row; col < indices.size(); ++col, ++idx )
00825                 {
00826                     col2 = std::lower_bound( m2Handles.begin(), m2Handles.end(), indices[col] ) - m2Handles.begin();
00827                     if( row2 <= col2 )
00828                     {
00829                         idx2 = indices.size() * row2 - row2 * ( row2 + 1 ) / 2 + col2;
00830                         if( !equal( Hessian[idx], m2Hess[idx2] ) ) ++hess_diff;
00831                     }
00832                     else
00833                     {
00834                         idx2 = indices.size() * col2 - col2 * ( col2 + 1 ) / 2 + row2;
00835                         if( !equal( Hessian[idx], transpose( m2Hess[idx2] ) ) ) ++hess_diff;
00836                     }
00837                 }
00838             }
00839 
00840             if( hess_diff )
00841             {
00842                 MSQ_SETERR( err )
00843                 ( MsqError::INVALID_STATE,
00844                   "Metrics returned different Hessian blocks for "
00845                   " %d of %u vertices for handle %lu in "
00846                   "evaluate_with_Hessian\n.",
00847                   hess_diff, (unsigned)Hessian.size(), (unsigned long)handle );
00848             }
00849         }
00850     }
00851 
00852     return r1 && !err;
00853 }
00854 
00855 void CompareMetric::get_mask_axis( PatchData& pd )
00856 {
00857     maskAxis          = -1;
00858     PlanarDomain* dom = reinterpret_cast< PlanarDomain* >( pd.get_domain() );
00859     int bits          = 0;
00860     if( dom )
00861     {
00862         Vector3D n = dom->get_normal();
00863         for( int i = 0; i < 3; ++i )
00864             if( fabs( n[i] ) < epsilon ) bits |= ( 1 << i );
00865         switch( bits )
00866         {
00867             case 3:
00868                 maskAxis = 2;
00869                 break;
00870             case 5:
00871                 maskAxis = 1;
00872                 break;
00873             case 6:
00874                 maskAxis = 0;
00875                 break;
00876         }
00877     }
00878 }
00879 
00880 bool CompareMetric::equal( Vector3D grad1, const Vector3D& grad2 ) const
00881 {
00882     if( maskAxis >= 0 ) grad1[maskAxis] = grad2[maskAxis];
00883     return ( grad1 - grad2 ).length_squared() <= epsilon * epsilon;
00884 }
00885 
00886 bool CompareMetric::equal( SymMatrix3D hess1, const SymMatrix3D& hess2 ) const
00887 {
00888     if( maskAxis >= 0 )
00889     {
00890         for( unsigned i = 0; i < 3; ++i )
00891         {
00892             hess1( maskAxis, i ) = hess2( maskAxis, i );
00893             hess1( i, maskAxis ) = hess2( i, maskAxis );
00894         }
00895     }
00896     return Frobenius_2( hess1 - hess2 ) <= epsilon * epsilon;
00897 }
00898 
00899 bool CompareMetric::equal( Matrix3D hess1, const Matrix3D& hess2 ) const
00900 {
00901     if( maskAxis >= 0 )
00902     {
00903         for( unsigned i = 0; i < 3; ++i )
00904         {
00905             hess1[maskAxis][i] = hess2[maskAxis][i];
00906             hess1[i][maskAxis] = hess2[i][maskAxis];
00907         }
00908     }
00909     return Frobenius_2( hess1 - hess2 ) <= epsilon * epsilon;
00910 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines