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