MOAB: Mesh Oriented datABase  (version 5.4.1)
TerminationCriterionTest.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     [email protected], [email protected], [email protected],
00024     [email protected], [email protected], [email protected]
00025 
00026   ***************************************************************** */
00027 
00028 /*! \file TerminationCriterionTest.cpp
00029     \author [email protected]
00030 
00031 Tests for the TerminationCriterion class..
00032 
00033 */
00034 
00035 #include "Mesquite.hpp"
00036 #include "MsqError.hpp"
00037 #include "Vector3D.hpp"
00038 #include "OFEvaluator.hpp"
00039 #include "ObjectiveFunction.hpp"
00040 #include "PatchData.hpp"
00041 #include "TerminationCriterion.hpp"
00042 #include "VertexMover.hpp"
00043 #include "VertexPatches.hpp"
00044 
00045 #include "UnitUtil.hpp"
00046 #include "PatchDataInstances.hpp"
00047 #include "ArrayMesh.hpp"
00048 #include "MeshUtil.hpp"
00049 #include "SimpleStats.hpp"
00050 #include "InstructionQueue.hpp"
00051 
00052 #include <iostream>
00053 #include <algorithm>
00054 #include <vector>
00055 #include <set>
00056 
00057 using namespace MBMesquite;
00058 
00059 class DummyOF : public ObjectiveFunction
00060 {
00061   public:
00062     double mValue;   //!< Objectve fuction value returned
00063     Vector3D mGrad;  //!< Gradient values for all vertices
00064     bool mValid;
00065 
00066     DummyOF( double of_value = 0.0, Vector3D grad_values = Vector3D( 0, 0, 0 ) )
00067         : mValue( 0.0 ), mGrad( grad_values ), mValid( true )
00068     {
00069     }
00070 
00071     bool initialize_block_coordinate_descent( MeshDomainAssoc* domain,
00072                                               const Settings* settings,
00073                                               PatchSet* user_set,
00074                                               MsqError& err );
00075 
00076     void initialize_queue( MeshDomainAssoc*, const Settings*, MsqError& ) {}
00077 
00078     bool evaluate( EvalType type, PatchData& pd, double& value_out, bool free, MsqError& err );
00079 
00080     bool evaluate_with_gradient( EvalType type,
00081                                  PatchData& pd,
00082                                  double& value_out,
00083                                  std::vector< Vector3D >& grad_out,
00084                                  MsqError& err );
00085 
00086     ObjectiveFunction* clone() const
00087     {
00088         return new DummyOF( *this );
00089     }
00090     void clear() {}
00091     int min_patch_layers() const
00092     {
00093         return 1;
00094     }
00095 };
00096 
00097 class TerminationCriterionTest : public CppUnit::TestFixture
00098 {
00099   private:
00100     CPPUNIT_TEST_SUITE( TerminationCriterionTest );
00101 
00102     CPPUNIT_TEST( test_number_of_iterates_inner );
00103     CPPUNIT_TEST( test_number_of_iterates_outer );
00104     CPPUNIT_TEST( test_cpu_time_inner );
00105     CPPUNIT_TEST( test_cpu_time_outer );
00106 
00107     CPPUNIT_TEST( test_absolute_vertex_movement );
00108     CPPUNIT_TEST( test_relative_vertex_movement );
00109     CPPUNIT_TEST( test_absolute_vertex_movement_edge_length );
00110 
00111     CPPUNIT_TEST( test_gradient_L2_norm_absolute );
00112     CPPUNIT_TEST( test_gradient_Linf_norm_absolute );
00113     CPPUNIT_TEST( test_gradient_L2_norm_relative );
00114     CPPUNIT_TEST( test_gradient_Linf_norm_relative );
00115 
00116     CPPUNIT_TEST( test_quality_improvement_absolute );
00117     CPPUNIT_TEST( test_quality_improvement_relative );
00118     CPPUNIT_TEST( test_successive_improvements_absolute );
00119     CPPUNIT_TEST( test_successive_improvements_relative );
00120 
00121     CPPUNIT_TEST( test_vertex_bound );
00122     CPPUNIT_TEST( test_untangled_mesh );
00123     CPPUNIT_TEST( test_abs_vtx_movement_culling );
00124 
00125     CPPUNIT_TEST_SUITE_END();
00126 
00127     DummyOF objFunc;
00128     OFEvaluator ofEval;
00129 
00130     void test_gradient_common( bool absolute, bool L2 );
00131     void test_quality_common( bool absolute, bool successive );
00132     void test_vertex_movement_common( bool absolute );
00133     void test_cpu_time_common( bool inner );
00134 
00135   public:
00136     TerminationCriterionTest() : ofEval( &objFunc ) {}
00137 
00138     // NUMBER OF ITERATES
00139     void test_number_of_iterates_inner();
00140     void test_number_of_iterates_outer();
00141 
00142     // CPU TIME
00143     void test_cpu_time_inner()
00144     {
00145         test_cpu_time_common( true );
00146     }
00147     void test_cpu_time_outer()
00148     {
00149         test_cpu_time_common( false );
00150     }
00151 
00152     // VERTEX MOVEMENT
00153     void test_absolute_vertex_movement()
00154     {
00155         test_vertex_movement_common( true );
00156     }
00157     void test_relative_vertex_movement()
00158     {
00159         test_vertex_movement_common( false );
00160     }
00161     void test_absolute_vertex_movement_edge_length();
00162 
00163     // GRADIENT NORM ABSOLUTE
00164     void test_gradient_L2_norm_absolute()
00165     {
00166         test_gradient_common( true, true );
00167     }
00168     void test_gradient_Linf_norm_absolute()
00169     {
00170         test_gradient_common( true, false );
00171     }
00172 
00173     // GRADIENT NORM RELATIVE
00174     void test_gradient_L2_norm_relative()
00175     {
00176         test_gradient_common( false, true );
00177     }
00178     void test_gradient_Linf_norm_relative()
00179     {
00180         test_gradient_common( false, false );
00181     }
00182 
00183     // QUALITY IMPROVEMENT ABSOLUTE
00184     void test_quality_improvement_absolute()
00185     {
00186         test_quality_common( true, false );
00187     }
00188 
00189     // QUALITY IMPROVEMENT RELATIVE
00190     void test_quality_improvement_relative()
00191     {
00192         test_quality_common( false, false );
00193     }
00194 
00195     // SUCCESSIVE IMPROVEMENTS ABSOLUTE
00196     void test_successive_improvements_absolute()
00197     {
00198         test_quality_common( true, true );
00199     }
00200 
00201     // SUCCESSIVE IMPROVEMENTS RELATIVE
00202     void test_successive_improvements_relative()
00203     {
00204         test_quality_common( false, true );
00205     }
00206 
00207     // VERTEX BOUND
00208     void test_vertex_bound();
00209 
00210     void test_untangled_mesh();
00211 
00212     // test culling on ABSOLUTE_VERTEX_MOVEMENT
00213     void test_abs_vtx_movement_culling();
00214 };
00215 
00216 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( TerminationCriterionTest, "TerminationCriterionTest" );
00217 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( TerminationCriterionTest, "Unit" );
00218 
00219 // NUMBER OF ITERATES
00220 void TerminationCriterionTest::test_number_of_iterates_inner()
00221 {
00222     MsqPrintError err( std::cout );
00223     PatchData pd;
00224     const int LIMIT = 2;
00225 
00226     TerminationCriterion tc;
00227     tc.add_iteration_limit( LIMIT );
00228     tc.reset_inner( pd, ofEval, err );
00229     ASSERT_NO_ERROR( err );
00230     tc.reset_patch( pd, err );
00231     ASSERT_NO_ERROR( err );
00232     for( int i = 0; i < LIMIT; ++i )
00233     {
00234         CPPUNIT_ASSERT( !tc.terminate() );
00235         CPPUNIT_ASSERT_EQUAL( i, tc.get_iteration_count() );
00236         tc.accumulate_patch( pd, err );
00237         ASSERT_NO_ERROR( err );
00238         tc.accumulate_inner( pd, 0, 0, err );
00239         ASSERT_NO_ERROR( err );
00240     }
00241     CPPUNIT_ASSERT_EQUAL( 2, tc.get_iteration_count() );
00242     CPPUNIT_ASSERT( tc.terminate() );
00243 }
00244 
00245 void TerminationCriterionTest::test_number_of_iterates_outer()
00246 {
00247     MsqPrintError err( std::cout );
00248     PatchData pd;
00249     const int LIMIT = 2;
00250 
00251     TerminationCriterion tc;
00252     tc.add_iteration_limit( LIMIT );
00253     tc.reset_outer( 0, 0, ofEval, 0, err );
00254     ASSERT_NO_ERROR( err );
00255     tc.reset_patch( pd, err );
00256     ASSERT_NO_ERROR( err );
00257     for( int i = 0; i < LIMIT; ++i )
00258     {
00259         CPPUNIT_ASSERT( !tc.terminate() );
00260         CPPUNIT_ASSERT_EQUAL( i, tc.get_iteration_count() );
00261         tc.accumulate_patch( pd, err );
00262         ASSERT_NO_ERROR( err );
00263         tc.accumulate_outer( 0, 0, ofEval, 0, err );
00264         ASSERT_NO_ERROR( err );
00265     }
00266     CPPUNIT_ASSERT_EQUAL( 2, tc.get_iteration_count() );
00267     CPPUNIT_ASSERT( tc.terminate() );
00268 }
00269 
00270 // CPU TIME
00271 void TerminationCriterionTest::test_cpu_time_common( bool inner )
00272 {
00273     MsqPrintError err( std::cout );
00274     PatchData pd;
00275     double LIMIT = 0.5;
00276 
00277     TerminationCriterion tc;
00278     tc.add_cpu_time( LIMIT );
00279     if( inner )
00280         tc.reset_inner( pd, ofEval, err );
00281     else
00282         tc.reset_outer( 0, 0, ofEval, 0, err );
00283     ASSERT_NO_ERROR( err );
00284     tc.reset_patch( pd, err );
00285     ASSERT_NO_ERROR( err );
00286 
00287     Timer timer;
00288     while( timer.since_birth() < 0.5 * LIMIT )
00289     {
00290         CPPUNIT_ASSERT( !tc.terminate() );
00291         tc.accumulate_patch( pd, err );
00292         ASSERT_NO_ERROR( err );
00293         if( inner )
00294             tc.accumulate_inner( pd, 0, 0, err );
00295         else
00296             tc.accumulate_outer( 0, 0, ofEval, 0, err );
00297         ASSERT_NO_ERROR( err );
00298     }
00299 
00300     while( timer.since_birth() < 1.1 * LIMIT )
00301         ;
00302 
00303     tc.accumulate_patch( pd, err );
00304     ASSERT_NO_ERROR( err );
00305     if( inner )
00306         tc.accumulate_inner( pd, 0, 0, err );
00307     else
00308         tc.accumulate_outer( 0, 0, ofEval, 0, err );
00309     ASSERT_NO_ERROR( err );
00310     CPPUNIT_ASSERT( tc.terminate() );
00311 }
00312 
00313 void TerminationCriterionTest::test_vertex_movement_common( bool absolute )
00314 {
00315     MsqPrintError err( std::cout );
00316     PatchData pd;
00317     create_twelve_hex_patch( pd, err );
00318     ASSERT_NO_ERROR( err );
00319 
00320     const double LIMIT = 1e-4;
00321     TerminationCriterion tc;
00322     if( absolute )
00323         tc.add_absolute_vertex_movement( LIMIT );
00324     else
00325         tc.add_relative_vertex_movement( LIMIT );
00326 
00327     tc.reset_inner( pd, ofEval, err );
00328     ASSERT_NO_ERROR( err );
00329     tc.reset_patch( pd, err );
00330     ASSERT_NO_ERROR( err );
00331     CPPUNIT_ASSERT( !tc.terminate() );
00332 
00333     const double FIRST_STEP = 10.0;
00334     // move a vertex by 10 units and check that it did not meet criterion
00335     pd.move_vertex( Vector3D( FIRST_STEP, 0, 0 ), 0, err );
00336     ASSERT_NO_ERROR( err );
00337     tc.accumulate_inner( pd, 0.0, 0, err );
00338     ASSERT_NO_ERROR( err );
00339     tc.accumulate_patch( pd, err );
00340     ASSERT_NO_ERROR( err );
00341     CPPUNIT_ASSERT( !tc.terminate() );
00342 
00343     double test_limit = LIMIT;
00344     if( !absolute ) test_limit *= FIRST_STEP;
00345 
00346     int idx = 0;
00347     for( double step = FIRST_STEP; step > test_limit; step *= 0.09 )
00348     {
00349         idx = ( idx + 1 ) % pd.num_free_vertices();
00350         pd.move_vertex( Vector3D( step, 0, 0 ), idx, err );
00351         ASSERT_NO_ERROR( err );
00352 
00353         tc.accumulate_inner( pd, 0.0, 0, err );
00354         ASSERT_NO_ERROR( err );
00355         tc.accumulate_patch( pd, err );
00356         ASSERT_NO_ERROR( err );
00357         CPPUNIT_ASSERT( !tc.terminate() );
00358     }
00359 
00360     idx = ( idx + 1 ) % pd.num_free_vertices();
00361     pd.move_vertex( Vector3D( 0.5 * test_limit, 0, 0 ), idx, err );
00362     ASSERT_NO_ERROR( err );
00363 
00364     tc.accumulate_inner( pd, 0.0, 0, err );
00365     ASSERT_NO_ERROR( err );
00366     tc.accumulate_patch( pd, err );
00367     ASSERT_NO_ERROR( err );
00368     CPPUNIT_ASSERT( tc.terminate() );
00369 }
00370 
00371 void TerminationCriterionTest::test_absolute_vertex_movement_edge_length()
00372 {
00373     MsqPrintError err( std::cout );
00374 
00375     // define two-tet mesh where tets share a face
00376     double coords[]            = { 0, -5, 0, 0, 5, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1 };
00377     const unsigned long conn[] = { 4, 3, 2, 0, 2, 3, 4, 1 };
00378     int fixed[5]               = { 0 };
00379     ArrayMesh mesh( 3, 5, coords, fixed, 2, TETRAHEDRON, conn );
00380 
00381     // calculate beta
00382     const double LIMIT = 1e-4;  // desired absolute limit
00383     MeshUtil tool( &mesh );
00384     SimpleStats len;
00385     tool.edge_length_distribution( len, err );
00386     ASSERT_NO_ERROR( err );
00387     const double beta = LIMIT / ( len.average() - len.standard_deviation() );
00388 
00389     // initialize termination criterion
00390     TerminationCriterion tc;
00391     tc.add_absolute_vertex_movement_edge_length( beta );
00392     MeshDomainAssoc mesh_and_domain2 = MeshDomainAssoc( &mesh, 0 );
00393     tc.initialize_queue( &mesh_and_domain2, 0, err );
00394     ASSERT_NO_ERROR( err );
00395 
00396     // get a patch data
00397     PatchData pd;
00398     pd.set_mesh( &mesh );
00399     pd.fill_global_patch( err );
00400     ASSERT_NO_ERROR( err );
00401 
00402     // test termination criteiorn
00403     tc.reset_inner( pd, ofEval, err );
00404     ASSERT_NO_ERROR( err );
00405     tc.reset_patch( pd, err );
00406     ASSERT_NO_ERROR( err );
00407     CPPUNIT_ASSERT( !tc.terminate() );
00408 
00409     const double FIRST_STEP = 10.0;
00410     // move a vertex by 10 units and check that it did not meet criterion
00411     pd.move_vertex( Vector3D( FIRST_STEP, 0, 0 ), 0, err );
00412     ASSERT_NO_ERROR( err );
00413     tc.accumulate_inner( pd, 0.0, 0, err );
00414     ASSERT_NO_ERROR( err );
00415     tc.accumulate_patch( pd, err );
00416     ASSERT_NO_ERROR( err );
00417     CPPUNIT_ASSERT( !tc.terminate() );
00418 
00419     double test_limit = LIMIT;
00420 
00421     int idx = 0;
00422     for( double step = FIRST_STEP; step > test_limit; step *= 0.09 )
00423     {
00424         idx = ( idx + 1 ) % pd.num_free_vertices();
00425         pd.move_vertex( Vector3D( step, 0, 0 ), idx, err );
00426         ASSERT_NO_ERROR( err );
00427 
00428         tc.accumulate_inner( pd, 0.0, 0, err );
00429         ASSERT_NO_ERROR( err );
00430         tc.accumulate_patch( pd, err );
00431         ASSERT_NO_ERROR( err );
00432         CPPUNIT_ASSERT( !tc.terminate() );
00433     }
00434 
00435     idx = ( idx + 1 ) % pd.num_free_vertices();
00436     pd.move_vertex( Vector3D( 0.5 * test_limit, 0, 0 ), idx, err );
00437     ASSERT_NO_ERROR( err );
00438 
00439     tc.accumulate_inner( pd, 0.0, 0, err );
00440     ASSERT_NO_ERROR( err );
00441     tc.accumulate_patch( pd, err );
00442     ASSERT_NO_ERROR( err );
00443     CPPUNIT_ASSERT( tc.terminate() );
00444 }
00445 
00446 static double lenfunc( const Vector3D* vect, int len )
00447 {
00448     return MBMesquite::length( vect, len );
00449 }
00450 static double maxfunc( const Vector3D* vect, int len )
00451 {
00452     return MBMesquite::Linf( vect, len );
00453 }
00454 
00455 void TerminationCriterionTest::test_gradient_common( bool absolute, bool L2 )
00456 {
00457     MsqPrintError err( std::cout );
00458     PatchData pd;
00459     create_twelve_hex_patch( pd, err );
00460     ASSERT_NO_ERROR( err );
00461 
00462     const double LIMIT = 1e-4;
00463     TerminationCriterion tc;
00464     if( absolute )
00465     {
00466         if( L2 )
00467             tc.add_absolute_gradient_L2_norm( LIMIT );
00468         else
00469             tc.add_absolute_gradient_inf_norm( LIMIT );
00470     }
00471     else
00472     {
00473         if( L2 )
00474             tc.add_relative_gradient_L2_norm( LIMIT );
00475         else
00476             tc.add_relative_gradient_inf_norm( LIMIT );
00477     }
00478 
00479     double ( *func_ptr )( const Vector3D*, int ) = L2 ? &lenfunc : &maxfunc;
00480 
00481     double junk, value = 1;
00482     objFunc.mGrad = Vector3D( value, value, value );
00483     tc.reset_inner( pd, ofEval, err );
00484     ASSERT_NO_ERROR( err );
00485     tc.reset_patch( pd, err );
00486     ASSERT_NO_ERROR( err );
00487 
00488     std::vector< Vector3D > grad;
00489     ofEval.evaluate( pd, junk, grad, err );
00490     ASSERT_NO_ERROR( err );
00491 
00492     double limit = LIMIT;
00493     if( !absolute ) limit *= func_ptr( arrptr( grad ), pd.num_free_vertices() );
00494 
00495     while( func_ptr( arrptr( grad ), pd.num_free_vertices() ) > limit )
00496     {
00497         CPPUNIT_ASSERT( !tc.terminate() );
00498 
00499         value *= 0.1;
00500         objFunc.mGrad = Vector3D( value, value, value );
00501         ofEval.evaluate( pd, junk, grad, err );
00502         ASSERT_NO_ERROR( err );
00503 
00504         tc.accumulate_inner( pd, 0.0, arrptr( grad ), err );
00505         ASSERT_NO_ERROR( err );
00506         tc.accumulate_patch( pd, err );
00507         ASSERT_NO_ERROR( err );
00508     }
00509 
00510     CPPUNIT_ASSERT( tc.terminate() );
00511 }
00512 
00513 static bool limit_absolute_quality( double, double, double curr_value, double epsilon )
00514 {
00515     return curr_value <= epsilon;
00516 }
00517 static bool limit_relative_quality( double init_value, double, double curr_value, double epsilon )
00518 {
00519     return curr_value <= epsilon * init_value;
00520 }
00521 static bool limit_absolute_sucessive( double, double prev_value, double curr_value, double epsilon )
00522 {
00523     return ( prev_value - curr_value ) <= epsilon;
00524 }
00525 static bool limit_relative_sucessive( double init_value, double prev_value, double curr_value, double epsilon )
00526 {
00527     return ( prev_value - curr_value ) <= epsilon * ( init_value - curr_value );
00528 }
00529 
00530 void TerminationCriterionTest::test_quality_common( bool absolute, bool successive )
00531 {
00532     MsqPrintError err( std::cout );
00533     PatchData pd;
00534     ASSERT_NO_ERROR( err );
00535 
00536     const double LIMIT = 1e-4;
00537     TerminationCriterion tc;
00538     bool ( *func_ptr )( double, double, double, double );
00539     if( absolute )
00540     {
00541         if( successive )
00542         {
00543             tc.add_absolute_successive_improvement( LIMIT );
00544             func_ptr = &limit_absolute_sucessive;
00545         }
00546         else
00547         {
00548             tc.add_absolute_quality_improvement( LIMIT );
00549             func_ptr = &limit_absolute_quality;
00550         }
00551     }
00552     else
00553     {
00554         if( successive )
00555         {
00556             tc.add_relative_successive_improvement( LIMIT );
00557             func_ptr = &limit_relative_sucessive;
00558         }
00559         else
00560         {
00561             tc.add_relative_quality_improvement( LIMIT );
00562             func_ptr = &limit_relative_quality;
00563         }
00564     }
00565 
00566     const double INIT_VALUE = 10.0;
00567     objFunc.mValue          = INIT_VALUE;
00568     tc.reset_inner( pd, ofEval, err );
00569     ASSERT_NO_ERROR( err );
00570     tc.reset_patch( pd, err );
00571     ASSERT_NO_ERROR( err );
00572     double prev = HUGE_VAL;
00573 
00574     while( !func_ptr( INIT_VALUE, prev, objFunc.mValue, LIMIT ) )
00575     {
00576         CPPUNIT_ASSERT( !tc.terminate() );
00577 
00578         prev = objFunc.mValue;
00579         objFunc.mValue *= 0.1;
00580         tc.accumulate_inner( pd, objFunc.mValue, 0, err );
00581         ASSERT_NO_ERROR( err );
00582         tc.accumulate_patch( pd, err );
00583         ASSERT_NO_ERROR( err );
00584     }
00585 
00586     CPPUNIT_ASSERT( tc.terminate() );
00587 }
00588 
00589 // VERTEX BOUND
00590 void TerminationCriterionTest::test_vertex_bound()
00591 {
00592     MsqPrintError err( std::cout );
00593     PatchData pd;
00594     create_twelve_hex_patch( pd, err );
00595     ASSERT_NO_ERROR( err );
00596 
00597     // get bounding dimension for patch
00598     double maxcoord = 0.0;
00599     for( size_t i = 0; i < pd.num_nodes(); ++i )
00600         for( int d = 0; d < 3; ++d )
00601             if( fabs( pd.vertex_by_index( i )[d] ) > maxcoord ) maxcoord = fabs( pd.vertex_by_index( i )[d] );
00602     // add a little bit for rounding error
00603     maxcoord += 1e-5;
00604 
00605     TerminationCriterion tc;
00606     tc.add_bounded_vertex_movement( maxcoord );
00607     tc.reset_inner( pd, ofEval, err );
00608     ASSERT_NO_ERROR( err );
00609     tc.reset_patch( pd, err );
00610     ASSERT_NO_ERROR( err );
00611     CPPUNIT_ASSERT( !tc.terminate() );
00612 
00613     int idx      = pd.num_free_vertices() - 1;
00614     Vector3D pos = pd.vertex_by_index( idx );
00615     pos[0]       = 2 * maxcoord;
00616     pd.set_vertex_coordinates( pos, idx, err );
00617     ASSERT_NO_ERROR( err );
00618     tc.accumulate_inner( pd, 0.0, 0, err );
00619     ASSERT_NO_ERROR( err );
00620     tc.accumulate_patch( pd, err );
00621     ASSERT_NO_ERROR( err );
00622     CPPUNIT_ASSERT( tc.terminate() );
00623 }
00624 
00625 // UNTANGLED
00626 void TerminationCriterionTest::test_untangled_mesh()
00627 {
00628     MsqPrintError err( std::cout );
00629     PatchData pd;
00630     create_twelve_hex_patch( pd, err );
00631     ASSERT_NO_ERROR( err );
00632 
00633     // get two opposite vertices in first hexahedral element
00634     int vtx1                 = pd.element_by_index( 0 ).get_vertex_index_array()[0];
00635     int vtx2                 = pd.element_by_index( 0 ).get_vertex_index_array()[7];
00636     Vector3D saved_coords    = pd.vertex_by_index( vtx2 );
00637     Vector3D opposite_coords = pd.vertex_by_index( vtx1 );
00638 
00639     // invert the element
00640     pd.move_vertex( 2 * ( opposite_coords - saved_coords ), vtx2, err );
00641     ASSERT_NO_ERROR( err );
00642     int inverted, samples;
00643     pd.element_by_index( 0 ).check_element_orientation( pd, inverted, samples, err );
00644     ASSERT_NO_ERROR( err );
00645     CPPUNIT_ASSERT( inverted > 0 );
00646 
00647     // check initial termination criterion
00648     TerminationCriterion tc;
00649     tc.add_untangled_mesh();
00650     tc.reset_inner( pd, ofEval, err );
00651     ASSERT_NO_ERROR( err );
00652     tc.reset_patch( pd, err );
00653     ASSERT_NO_ERROR( err );
00654     CPPUNIT_ASSERT( !tc.terminate() );
00655 
00656     // fix the element
00657     pd.set_vertex_coordinates( saved_coords, vtx2, err );
00658     ASSERT_NO_ERROR( err );
00659     pd.element_by_index( 0 ).check_element_orientation( pd, inverted, samples, err );
00660     ASSERT_NO_ERROR( err );
00661     CPPUNIT_ASSERT_EQUAL( 0, inverted );
00662 
00663     // check that TC recognized untangled mesh
00664     tc.accumulate_inner( pd, 0.0, 0, err );
00665     ASSERT_NO_ERROR( err );
00666     tc.accumulate_patch( pd, err );
00667     ASSERT_NO_ERROR( err );
00668     CPPUNIT_ASSERT( tc.terminate() );
00669 }
00670 
00671 class TCTFauxOptimizer : public VertexMover
00672 {
00673   public:
00674     TCTFauxOptimizer( double pertubation_amount ) : mDelta( pertubation_amount ) {}
00675     virtual ~TCTFauxOptimizer(){};
00676     virtual std::string get_name() const
00677     {
00678         return "Optimizer for TerminationCriterionTest";
00679     }
00680     virtual PatchSet* get_patch_set()
00681     {
00682         return &mPatchSet;
00683     }
00684     virtual void initialize( PatchData& pd, MsqError& err );
00685     virtual void initialize_mesh_iteration( PatchData& pd, MsqError& err ) {}
00686     virtual void optimize_vertex_positions( PatchData& pd, MsqError& err );
00687     virtual void terminate_mesh_iteration( PatchData& pd, MsqError& err ) {}
00688     virtual void cleanup()
00689     {
00690         all.clear();
00691     }
00692     int num_passes() const
00693     {
00694         return numPasses;
00695     }
00696     bool should_have_terminated() const
00697     {
00698         return perturbFrac > (int)all.size();
00699     }
00700 
00701   private:
00702     std::set< Mesh::VertexHandle > culled, visited;
00703     std::vector< Mesh::VertexHandle > all;
00704     VertexPatches mPatchSet;
00705     int numPasses;    //!< count number of outer iterations
00706     int perturbFrac;  //!< perturb 1/perturbFrac of the free vertices
00707     double mDelta;    //!< distance to perturb vertices
00708 };
00709 
00710 void TCTFauxOptimizer::initialize( PatchData& pd, MsqError& err )
00711 {
00712     CPPUNIT_ASSERT( all.empty() );
00713     culled.clear();
00714     visited.clear();
00715     numPasses = 1;
00716 
00717     pd.get_mesh()->get_all_vertices( all, err );MSQ_ERRRTN( err );
00718     std::vector< bool > fixed;
00719     pd.get_mesh()->vertices_get_fixed_flag( &all[0], fixed, all.size(), err );
00720     size_t w = 0;
00721     for( size_t r = 0; r < all.size(); ++r )
00722         if( !fixed[r] ) all[w++] = all[r];
00723     all.resize( w );MSQ_ERRRTN( err );
00724 
00725     perturbFrac = 1;
00726 }
00727 
00728 void TCTFauxOptimizer::optimize_vertex_positions( PatchData& pd, MsqError& err )
00729 {
00730     Mesh::VertexHandle free_vtx = pd.get_vertex_handles_array()[0];
00731     if( visited.insert( free_vtx ).second == false )
00732     {  // already visited this one
00733         // The inner termination criterion should include an iteration limit of 1.
00734         // So if we are seeing the same vertex again, this means that we *should*
00735         // be stating a new pass over the mesh.
00736 
00737         // We are presumably starting a new pass over the mesh.
00738         // Verify that we visisted all of the free, non-culled vertices
00739         for( size_t i = 0; i < all.size(); ++i )
00740         {
00741             if( culled.find( all[i] ) == culled.end() )
00742             {
00743                 if( visited.find( all[i] ) == visited.end() )
00744                 {
00745                     std::ostringstream str;
00746                     str << "Did not visit vertex " << i << " (handle " << all[i] << ") in pass " << numPasses
00747                         << std::endl;
00748                     CPPUNIT_FAIL( str.str() );
00749                 }
00750             }
00751         }
00752         visited.clear();
00753         visited.insert( free_vtx );
00754         ++numPasses;
00755 
00756         // Check that we terminate when expected
00757         CPPUNIT_ASSERT( !should_have_terminated() );
00758 
00759         perturbFrac *= 2;  // for each pass, perturb half as many vertices
00760     }
00761 
00762     // check that we are not visiting a culled vertex
00763     CPPUNIT_ASSERT( culled.find( free_vtx ) == culled.end() );
00764 
00765     // for each pass, perturb half as many vertices
00766     size_t idx = std::find( all.begin(), all.end(), free_vtx ) - all.begin();
00767     CPPUNIT_ASSERT( idx < all.size() );  // not a free vertex????
00768     if( 0 == ( ( idx + 1 ) % perturbFrac ) )
00769     {
00770         // perturb vertex
00771         double sign = numPasses % 2 == 0 ? 1 : -1;
00772         Vector3D delta( sign * mDelta, 0, 0 );
00773         pd.move_vertex( delta, 0, err );
00774         ASSERT_NO_ERROR( err );
00775         // any adjacent vertices should not be culled
00776         for( size_t i = 0; i < pd.num_nodes(); ++i )
00777             culled.erase( pd.get_vertex_handles_array()[i] );
00778     }
00779     else
00780     {
00781         // If we're not moving this vertex, then it should get culled
00782         culled.insert( free_vtx );
00783     }
00784 }
00785 
00786 void TerminationCriterionTest::test_abs_vtx_movement_culling()
00787 {
00788     /* define a quad mesh like the following
00789       16--17--18--19--20--21--22--23
00790        |   |   |   |   |   |   |   |   Y
00791        8---9--10--11--12--13--14--15   ^
00792        |   |   |   |   |   |   |   |   |
00793        0---1---2---3---4---5---6---7   +-->X
00794     */
00795 
00796     const int nvtx = 24;
00797     int fixed[nvtx];
00798     double coords[3 * nvtx];
00799     for( int i = 0; i < nvtx; ++i )
00800     {
00801         coords[3 * i]     = i / 8;
00802         coords[3 * i + 1] = i % 8;
00803         coords[3 * i + 2] = 0;
00804         fixed[i]          = i < 9 || i > 14;
00805     }
00806 
00807     const int nquad = 14;
00808     unsigned long conn[4 * nquad];
00809     for( int i = 0; i < nquad; ++i )
00810     {
00811         int row         = i / 7;
00812         int idx         = i % 7;
00813         int n0          = 8 * row + idx;
00814         conn[4 * i]     = n0;
00815         conn[4 * i + 1] = n0 + 1;
00816         conn[4 * i + 2] = n0 + 9;
00817         conn[4 * i + 3] = n0 + 8;
00818     }
00819 
00820     const double tol = 0.05;
00821     PlanarDomain zplane( PlanarDomain::XY, 0.0 );
00822     ArrayMesh mesh( 3, nvtx, coords, fixed, nquad, QUADRILATERAL, conn );
00823 
00824     // fill vertex byte with garbage to make sure that quality improver
00825     // is initializing correctly.
00826     MsqError err;
00827     std::vector< unsigned char > junk( nvtx, 255 );
00828     std::vector< Mesh::VertexHandle > vertices;
00829     mesh.get_all_vertices( vertices, err );
00830     ASSERT_NO_ERROR( err );
00831     CPPUNIT_ASSERT_EQUAL( junk.size(), vertices.size() );
00832     mesh.vertices_set_byte( &vertices[0], &junk[0], vertices.size(), err );
00833     ASSERT_NO_ERROR( err );
00834 
00835     // Define optimizer
00836     TCTFauxOptimizer smoother( 2 * tol );
00837     TerminationCriterion outer, inner;
00838     outer.add_absolute_vertex_movement( tol );
00839     inner.cull_on_absolute_vertex_movement( tol );
00840     inner.add_iteration_limit( 1 );
00841     smoother.set_inner_termination_criterion( &inner );
00842     smoother.set_outer_termination_criterion( &outer );
00843 
00844     // No run the "optimizer"
00845     InstructionQueue q;
00846     q.set_master_quality_improver( &smoother, err );
00847     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &zplane );
00848     q.run_instructions( &mesh_and_domain, err );
00849     ASSERT_NO_ERROR( err );
00850     CPPUNIT_ASSERT( smoother.should_have_terminated() );
00851     CPPUNIT_ASSERT( smoother.num_passes() > 1 );
00852 }
00853 
00854 bool DummyOF::initialize_block_coordinate_descent( MeshDomainAssoc*, const Settings*, PatchSet*, MsqError& err )
00855 {
00856     MSQ_SETERR( err )( MsqError::NOT_IMPLEMENTED );
00857     return false;
00858 }
00859 
00860 bool DummyOF::evaluate( EvalType, PatchData&, double& value, bool, MsqError& )
00861 {
00862     value = mValue;
00863     return mValid;
00864 }
00865 
00866 bool DummyOF::evaluate_with_gradient( EvalType,
00867                                       PatchData& pd,
00868                                       double& value_out,
00869                                       std::vector< Vector3D >& grad_out,
00870                                       MsqError& )
00871 {
00872     value_out = mValue;
00873     grad_out.clear();
00874     grad_out.resize( pd.num_free_vertices(), mGrad );
00875     return mValid;
00876 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines