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