MOAB: Mesh Oriented datABase  (version 5.2.1)
NonSmoothDescent.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 // -*- Mode : c++; tab-width: 3; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 3
00028 // -*-
00029 
00030 /*!  \File NonSmoothDescent.cpp \brief
00031 
00032   Implements the NonSmoothDescent class member functions.
00033 
00034   \author Lori Freitag
00035   \date 2002-07-20 */
00036 
00037 #include <cstdlib>
00038 #include <cstdio>
00039 #include "NonSmoothDescent.hpp"
00040 #include "MsqTimer.hpp"
00041 
00042 #undef NUMERICAL_GRADIENT
00043 
00044 namespace MBMesquite
00045 {
00046 
00047 const double EPSILON       = 1e-16;
00048 const int MSQ_MAX_OPT_ITER = 20;
00049 
00050 enum Rotate
00051 {
00052     COUNTERCLOCKWISE = 1,
00053     CLOCKWISE        = 0
00054 };
00055 
00056 NonSmoothDescent::NonSmoothDescent( QualityMetric* qm ) : currentQM( qm )
00057 {
00058 
00059     MSQ_DBGOUT( 1 ) << "- Executed NonSmoothDescent::NonSmoothDescent()\n";
00060 }
00061 
00062 std::string NonSmoothDescent::get_name() const
00063 {
00064     return "NonSmoothDescent";
00065 }
00066 
00067 PatchSet* NonSmoothDescent::get_patch_set()
00068 {
00069     return &patchSet;
00070 }
00071 
00072 void NonSmoothDescent::initialize( PatchData& /*pd*/, MsqError& /*err*/ )
00073 {
00074     minStepSize = 1e-6;
00075     MSQ_DBGOUT( 1 ) << "- Executed NonSmoothDescent::initialize()\n";
00076 }
00077 
00078 void NonSmoothDescent::initialize_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
00079 void NonSmoothDescent::optimize_vertex_positions( PatchData& pd, MsqError& err )
00080 {
00081     MSQ_FUNCTION_TIMER( "NonSmoothDescent" );
00082 
00083     //  cout << "- Executing NonSmoothDescent::optimize_node_positions()\n";
00084     /* perform the min max smoothing algorithm */
00085     MSQ_PRINT( 2 )( "\nInitializing the patch iteration\n" );
00086 
00087     MSQ_PRINT( 3 )( "Number of Vertices: %d\n", (int)pd.num_nodes() );
00088     MSQ_PRINT( 3 )( "Number of Elements: %d\n", (int)pd.num_elements() );
00089     // Michael: Note: is this a reliable way to get the dimension?
00090     // NOTE: Mesquite only does 3-dimensional (including surface) meshes.
00091     // mDimension = pd.get_mesh()->get_geometric_dimension(err); MSQ_ERRRTN(err);
00092     // MSQ_PRINT(3)("Spatial Dimension: %d\n",mDimension);
00093     MSQ_PRINT( 3 )( "Spatial Dimension: 3\n" );
00094 
00095     MSQ_PRINT( 3 )( "Num Free = %d\n", (int)pd.num_free_vertices() );
00096 
00097     MsqFreeVertexIndexIterator free_iter( pd, err );MSQ_ERRRTN( err );
00098     free_iter.reset();
00099     free_iter.next();
00100     freeVertexIndex = free_iter.value();
00101     MSQ_PRINT( 3 )( "Free Vertex Index = %lu\n", (unsigned long)freeVertexIndex );
00102 
00103     // TODO - need to switch to validity via metric evaluations should
00104     // be associated with the compute_function somehow
00105     /* check for an invalid mesh; if it's invalid return and ask the user
00106        to use untangle */
00107     if( !this->validity_check( pd, err ) )
00108     {
00109         MSQ_PRINT( 1 )( "ERROR: Invalid mesh\n" );
00110         MSQ_SETERR( err )
00111         ( "Invalid Mesh: Use untangle to create a valid "
00112           "triangulation",
00113           MsqError::INVALID_MESH );
00114         return;
00115     }
00116 
00117     /* initialize the optimization data up to numFunctionValues */
00118     this->init_opt( pd, err );MSQ_ERRRTN( err );
00119     this->init_max_step_length( pd, err );MSQ_ERRRTN( err );
00120     MSQ_PRINT( 3 )( "Done initializing optimization\n" );
00121 
00122     /* compute the initial function values */
00123     // TODO this should return a bool with the validity
00124     this->compute_function( &pd, originalFunction, err );MSQ_ERRRTN( err );
00125 
00126     // find the initial active set
00127     this->find_active_set( originalFunction, mActive );
00128 
00129     this->minmax_opt( pd, err );MSQ_ERRRTN( err );
00130 }
00131 
00132 void NonSmoothDescent::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
00133 
00134 void NonSmoothDescent::cleanup()
00135 {
00136     MSQ_DBGOUT( 1 ) << "- Executing NonSmoothDescent::cleanup()\n";
00137     MSQ_DBGOUT( 1 ) << "- Done with NonSmoothDescent::cleanup()\n";
00138 }
00139 
00140 void NonSmoothDescent::find_plane_points( Direction dir1, Direction dir2, const std::vector< Vector3D >& vec,
00141                                           Vector3D& pt1, Vector3D& pt2, Vector3D& /*pt3*/, Status& status, MsqError& )
00142 {
00143     int i;
00144     int num_min, num_max;
00145     Rotate rotate   = CLOCKWISE;
00146     int num_rotated = 0;
00147     double pt_1, pt_2;
00148     double min, inv_slope;
00149     double min_inv_slope = 0.;
00150     double max;
00151     double max_inv_slope    = 0;
00152     double inv_origin_slope = 0;
00153     const int num_vec       = vec.size();
00154     const int auto_ind_size = 50;
00155     int auto_ind[auto_ind_size];
00156     std::vector< int > heap_ind;
00157     int* ind;
00158     if( num_vec <= auto_ind_size )
00159         ind = auto_ind;
00160     else
00161     {
00162         heap_ind.resize( num_vec );
00163         ind = &heap_ind[0];
00164     }
00165 
00166     status = MSQ_CHECK_BOTTOM_UP;
00167     /* find the minimum points in dir1 starting at -1 */
00168     num_min = 0;
00169     ind[0]  = -1;
00170     ind[1]  = -1;
00171     ind[2]  = -1;
00172     min     = 1.0;
00173     for( i = 0; i < num_vec; i++ )
00174     {
00175         if( vec[i][dir1] < min )
00176         {
00177             min     = vec[i][dir1];
00178             ind[0]  = i;
00179             num_min = 1;
00180         }
00181         else if( fabs( vec[i][dir1] - min ) < EPSILON )
00182         {
00183             ind[num_min++] = i;
00184         }
00185     }
00186     if( min >= 0 ) status = MSQ_NO_EQUIL;
00187 
00188     if( status != MSQ_NO_EQUIL )
00189     {
00190         switch( num_min )
00191         {
00192             case 1: /* rotate to find the next point */
00193                 pt1  = vec[ind[0]];
00194                 pt_1 = pt1[dir1];
00195                 pt_2 = pt1[dir2];
00196                 if( pt1[dir2] <= 0 )
00197                 {
00198                     rotate        = COUNTERCLOCKWISE;
00199                     max_inv_slope = -HUGE_VAL;
00200                 }
00201                 if( pt1[dir2] > 0 )
00202                 {
00203                     rotate        = CLOCKWISE;
00204                     min_inv_slope = HUGE_VAL;
00205                 }
00206                 switch( rotate )
00207                 {
00208                     case COUNTERCLOCKWISE:
00209                         for( i = 0; i < num_vec; i++ )
00210                         {
00211                             if( i != ind[0] )
00212                             {
00213                                 inv_slope = ( vec[i][dir2] - pt_2 ) / ( vec[i][dir1] - pt_1 );
00214                                 if( ( inv_slope > max_inv_slope ) && ( fabs( inv_slope - max_inv_slope ) > EPSILON ) )
00215                                 {
00216                                     ind[1]        = i;
00217                                     max_inv_slope = inv_slope;
00218                                     num_rotated   = 1;
00219                                 }
00220                                 else if( fabs( inv_slope - max_inv_slope ) < EPSILON )
00221                                 {
00222                                     ind[2] = i;
00223                                     num_rotated++;
00224                                 }
00225                             }
00226                         }
00227                         break;
00228                     case CLOCKWISE:
00229                         for( i = 0; i < num_vec; i++ )
00230                         {
00231                             if( i != ind[0] )
00232                             {
00233                                 inv_slope = ( vec[i][dir2] - pt_2 ) / ( vec[i][dir1] - pt_1 );
00234                                 if( ( inv_slope < min_inv_slope ) && ( fabs( inv_slope - max_inv_slope ) > EPSILON ) )
00235                                 {
00236                                     ind[1]        = i;
00237                                     min_inv_slope = inv_slope;
00238                                     num_rotated   = 1;
00239                                 }
00240                                 else if( fabs( inv_slope - min_inv_slope ) < EPSILON )
00241                                 {
00242                                     ind[2] = i;
00243                                     num_rotated++;
00244                                 }
00245                             }
00246                         }
00247                 }
00248                 switch( num_rotated )
00249                 {
00250                     case 0:
00251                         MSQ_PRINT( 3 )( "No points in the rotation ... odd\n" );
00252                         status = MSQ_HULL_TEST_ERROR;
00253                         break;
00254                     case 1:
00255                         MSQ_PRINT( 3 )( "Found a line in the convex hull\n" );
00256                         pt2    = vec[ind[1]];
00257                         status = MSQ_TWO_PT_PLANE;
00258                         break;
00259                     default:
00260                         MSQ_PRINT( 3 )( "Found 2 or more points in the rotation\n" );
00261                         if( fabs( pt_1 ) > EPSILON ) inv_origin_slope = pt_2 / pt_1;
00262                         switch( rotate )
00263                         {
00264                             case COUNTERCLOCKWISE:
00265                                 if( inv_origin_slope >= max_inv_slope )
00266                                     status = MSQ_NO_EQUIL;
00267                                 else
00268                                     status = MSQ_CHECK_TOP_DOWN;
00269                                 break;
00270                             case CLOCKWISE:
00271                                 if( inv_origin_slope <= min_inv_slope )
00272                                     status = MSQ_NO_EQUIL;
00273                                 else
00274                                     status = MSQ_CHECK_TOP_DOWN;
00275                         }
00276                 }
00277                 break;
00278             case 2: /* use these two points to define the plane */
00279                 MSQ_PRINT( 3 )( "Found two minimum points to define the plane\n" );
00280                 pt1    = vec[ind[0]];
00281                 pt2    = vec[ind[1]];
00282                 status = MSQ_TWO_PT_PLANE;
00283                 break;
00284             default: /* check to see if all > 0 */
00285                 MSQ_PRINT( 3 )( "Found 3 or more points in min plane %f\n", min );
00286                 if( vec[ind[0]][dir1] >= 0 )
00287                     status = MSQ_NO_EQUIL;
00288                 else
00289                     status = MSQ_CHECK_TOP_DOWN;
00290         }
00291     }
00292 
00293     /***************************/
00294     /*  failed to find any information, checking top/down this coord*/
00295     /***************************/
00296 
00297     if( status == MSQ_CHECK_TOP_DOWN )
00298     {
00299         /* find the maximum points in dir1 starting at 1 */
00300         num_max = 0;
00301         ind[0]  = -1;
00302         ind[1]  = -1;
00303         ind[2]  = -1;
00304         max     = -1.0;
00305         for( i = 0; i < num_vec; i++ )
00306         {
00307             if( vec[i][dir1] > max )
00308             {
00309                 max     = vec[i][dir1];
00310                 ind[0]  = i;
00311                 num_max = 1;
00312             }
00313             else if( fabs( vec[i][dir1] - max ) < EPSILON )
00314             {
00315                 ind[num_max++] = i;
00316             }
00317         }
00318         if( max <= 0 ) status = MSQ_NO_EQUIL;
00319 
00320         if( status != MSQ_NO_EQUIL )
00321         {
00322             switch( num_max )
00323             {
00324                 case 1: /* rotate to find the next point */
00325                     pt1  = vec[ind[0]];
00326                     pt_1 = pt1[dir1];
00327                     pt_2 = pt1[dir2];
00328                     if( pt1[dir2] < 0 )
00329                     {
00330                         rotate        = CLOCKWISE;
00331                         min_inv_slope = HUGE_VAL;
00332                     }
00333                     if( pt1[dir2] >= 0 )
00334                     {
00335                         rotate        = COUNTERCLOCKWISE;
00336                         max_inv_slope = -HUGE_VAL;
00337                     }
00338                     switch( rotate )
00339                     {
00340                         case COUNTERCLOCKWISE:
00341                             for( i = 0; i < num_vec; i++ )
00342                             {
00343                                 if( i != ind[0] )
00344                                 {
00345                                     inv_slope = ( vec[i][dir2] - pt_2 ) / ( vec[i][dir1] - pt_1 );
00346                                     if( inv_slope > max_inv_slope )
00347                                     {
00348                                         ind[1]        = i;
00349                                         max_inv_slope = inv_slope;
00350                                         num_rotated   = 1;
00351                                     }
00352                                     else if( fabs( inv_slope - max_inv_slope ) < EPSILON )
00353                                     {
00354                                         ind[2] = i;
00355                                         num_rotated++;
00356                                     }
00357                                 }
00358                             }
00359                             break;
00360                         case CLOCKWISE:
00361                             for( i = 0; i < num_vec; i++ )
00362                             {
00363                                 if( i != ind[0] )
00364                                 {
00365                                     inv_slope = ( vec[i][dir2] - pt_2 ) / ( vec[i][dir1] - pt_1 );
00366                                     if( inv_slope < min_inv_slope )
00367                                     {
00368                                         ind[1]        = i;
00369                                         min_inv_slope = inv_slope;
00370                                         num_rotated   = 1;
00371                                     }
00372                                     else if( fabs( inv_slope - min_inv_slope ) < EPSILON )
00373                                     {
00374                                         ind[2] = i;
00375                                         num_rotated++;
00376                                     }
00377                                 }
00378                             }
00379                     }
00380                     switch( num_rotated )
00381                     {
00382                         case 0:
00383                             MSQ_PRINT( 3 )( "No points in the rotation ... odd\n" );
00384                             status = MSQ_HULL_TEST_ERROR;
00385                             break;
00386                         case 1:
00387                             MSQ_PRINT( 3 )( "Found a line in the convex hull\n" );
00388                             pt2    = vec[ind[1]];
00389                             status = MSQ_TWO_PT_PLANE;
00390                             break;
00391                         default:
00392                             MSQ_PRINT( 3 )( "Found 2 or more points in the rotation\n" );
00393                             /* check to see if rotation got past origin */
00394                             inv_origin_slope = pt_2 / pt_1;
00395                             switch( rotate )
00396                             {
00397                                 case COUNTERCLOCKWISE:
00398                                     if( inv_origin_slope >= max_inv_slope )
00399                                         status = MSQ_NO_EQUIL;
00400                                     else if( dir1 == 2 )
00401                                         status = MSQ_CHECK_Y_COORD_DIRECTION;
00402                                     else if( dir1 == 1 )
00403                                         status = MSQ_CHECK_X_COORD_DIRECTION;
00404                                     else
00405                                         status = MSQ_EQUIL;
00406                                     break;
00407                                 case CLOCKWISE:
00408                                     if( inv_origin_slope <= min_inv_slope )
00409                                         status = MSQ_NO_EQUIL;
00410                                     else if( dir1 == 2 )
00411                                         status = MSQ_CHECK_Y_COORD_DIRECTION;
00412                                     else if( dir1 == 1 )
00413                                         status = MSQ_CHECK_X_COORD_DIRECTION;
00414                                     else
00415                                         status = MSQ_EQUIL;
00416                             }
00417                     }
00418                     break;
00419                 case 2: /* use these two points to define the plane */
00420                     pt1    = vec[ind[0]];
00421                     pt2    = vec[ind[1]];
00422                     status = MSQ_TWO_PT_PLANE;
00423                     break;
00424                 default: /* check to see if all > 0 */
00425                     MSQ_PRINT( 3 )( "Found 3 in max plane %f\n", max );
00426                     if( vec[ind[0]][dir1] <= 0 )
00427                         status = MSQ_NO_EQUIL;
00428                     else if( dir1 == 2 )
00429                         status = MSQ_CHECK_Y_COORD_DIRECTION;
00430                     else if( dir1 == 1 )
00431                         status = MSQ_CHECK_X_COORD_DIRECTION;
00432                     else
00433                         status = MSQ_EQUIL;
00434             }
00435         }
00436     }
00437 }
00438 
00439 void NonSmoothDescent::search_direction( PatchData& /*pd*/, Vector3D& mSearch, MsqError& err )
00440 {
00441     bool viable;
00442     double a, b, c, denom;
00443     std::vector< Vector3D > dir;
00444     double R0, R1;
00445     SymmetricMatrix P;
00446     double x[2];
00447     double search_mag;
00448 
00449     const int num_active = mActive.active_ind.size();
00450     ;
00451 
00452     // TODO This might be o.k. actually - i don't see any dependence
00453     // on the element geometry here... try it and see if it works.
00454     // if not, try taking all of the gradients in the active set
00455     // and let the search direction be the average of those.
00456     //   MSQ_FUNCTION_TIMER( "Search Direction" );
00457 
00458     MSQ_PRINT( 2 )( "\nIn Search Direction\n" );
00459     this->print_active_set( mActive, mFunction );
00460 
00461     if( num_active == 0 )
00462     {
00463         MSQ_SETERR( err )( "No active values in search", MsqError::INVALID_STATE );
00464         return;
00465     }
00466 
00467     switch( num_active )
00468     {
00469         case 1:
00470             mSearch   = mGradient[mActive.active_ind[0]];
00471             mSteepest = mActive.active_ind[0];
00472             break;
00473         case 2:
00474             /* if there are two active points, move in the direction of the
00475            intersection of the planes.  This is the steepest descent
00476                direction found by analytically solving the QP */
00477 
00478             /* set up the active gradient directions */
00479             this->get_active_directions( mGradient, dir, err );MSQ_ERRRTN( err );
00480 
00481             /* form the grammian */
00482             this->form_grammian( dir, err );MSQ_ERRRTN( err );
00483             this->form_PD_grammian( err );MSQ_ERRRTN( err );
00484 
00485             denom  = ( mG( 0, 0 ) + mG( 1, 1 ) - 2 * mG( 0, 1 ) );
00486             viable = true;
00487             if( fabs( denom ) > EPSILON )
00488             {
00489                 /* gradients are LI, move along their intersection */
00490                 b = ( mG( 0, 0 ) - mG( 0, 1 ) ) / denom;
00491                 a = 1 - b;
00492                 if( ( b < 0 ) || ( b > 1 ) ) viable = false; /* 0 < b < 1 */
00493                 if( viable ) { mSearch = a * dir[0] + b * dir[1]; }
00494                 else
00495                 {
00496                     /* the gradients are dependent, move along one face */
00497                     mSearch = dir[0];
00498                 }
00499             }
00500             else
00501             {
00502                 /* the gradients are dependent, move along one face */
00503                 mSearch = dir[0];
00504             }
00505             mSteepest = mActive.active_ind[0];
00506 
00507             break;
00508         default:
00509             /* as in case 2: solve the QP problem to find the steepest
00510                descent direction.  This can be done analytically - as
00511                is done in Gill, Murray and Wright
00512                  for 3 active points in 3 directions - test PD of G
00513                  otherwise we know it's SP SD so search edges and faces */
00514 
00515             /* get the active gradient directions */
00516             this->get_active_directions( mGradient, dir, err );MSQ_ERRRTN( err );
00517 
00518             /* form the entries of the grammian matrix */
00519             this->form_grammian( dir, err );MSQ_ERRRTN( err );
00520             this->form_PD_grammian( err );MSQ_ERRRTN( err );
00521 
00522             if( num_active == 3 )
00523             {
00524                 if( mG.condition3x3() < 1e14 )
00525                 {  // if not singular
00526                     /* form the entries of P=Z^T G Z where Z = [-1...-1; I ] */
00527                     this->form_reduced_matrix( P );
00528                     /* form  the RHS and solve the system for the coeffs */
00529                     R0      = mG( 0, 0 ) - mG( 1, 0 );
00530                     R1      = mG( 0, 0 ) - mG( 2, 0 );
00531                     bool ok = this->solve2x2( P( 0, 0 ), P( 0, 1 ), P( 1, 0 ), P( 1, 1 ), R0, R1, x, err );MSQ_ERRRTN( err );
00532                     if( ok )
00533                     {
00534                         a         = 1 - x[0] - x[1];
00535                         b         = x[0];
00536                         c         = x[1];
00537                         mSearch   = a * dir[0] + b * dir[1] + c * dir[2];
00538                         mSteepest = mActive.active_ind[0];
00539                     }
00540                     else
00541                     {
00542                         this->search_edges_faces( &dir[0], mSearch, err );MSQ_ERRRTN( err );
00543                     }
00544                 }
00545                 else
00546                 {
00547                     this->search_edges_faces( &dir[0], mSearch, err );MSQ_ERRRTN( err );
00548                 }
00549             }
00550             else
00551             {
00552                 this->search_edges_faces( &dir[0], mSearch, err );MSQ_ERRRTN( err );
00553             }
00554     }
00555 
00556     /* if the search direction is essentially zero, equilibrium pt */
00557     search_mag = mSearch % mSearch;
00558     MSQ_PRINT( 3 )( "  Search Magnitude %g \n", search_mag );
00559 
00560     if( fabs( search_mag ) < 1E-13 )
00561         optStatus = MSQ_ZERO_SEARCH;
00562     else
00563         mSearch /= std::sqrt( search_mag );
00564     MSQ_PRINT( 3 )
00565     ( "  Search Direction %g %g  Steepest %lu\n", mSearch[0], mSearch[1], (unsigned long)mSteepest );
00566 }
00567 
00568 void NonSmoothDescent::minmax_opt( PatchData& pd, MsqError& err )
00569 {
00570     bool equilibriumPt = false;
00571     Vector3D mSearch( 0, 0, 0 );
00572 
00573     //      int valid;
00574     MSQ_FUNCTION_TIMER( "Minmax Opt" );
00575     MSQ_PRINT( 2 )( "In minmax_opt\n" );
00576 
00577     mFunction     = originalFunction;
00578     originalValue = mActive.true_active_value;
00579 
00580     int iterCount    = 0;
00581     int optIterCount = 0;
00582 
00583     MSQ_PRINT( 3 )( "Done copying original function to function\n" );
00584 
00585     this->find_active_set( mFunction, mActive );
00586     prevActiveValues.clear();
00587     prevActiveValues.push_back( mActive.true_active_value );
00588 
00589     /* check for equilibrium point */
00590     /* compute the gradient */
00591     this->compute_gradient( &pd, mGradient, err );MSQ_ERRRTN( err );
00592 
00593     if( mActive.active_ind.size() >= 2 )
00594     {
00595         MSQ_PRINT( 3 )( "Testing for an equilibrium point \n" );
00596         equilibriumPt = this->check_equilibrium( optStatus, err );MSQ_ERRRTN( err );
00597 
00598         if( MSQ_DBG( 2 ) && equilibriumPt ) MSQ_PRINT( 2 )( "Optimization Exiting: An equilibrium point \n" );
00599     }
00600 
00601     /* terminate if we have found an equilibrium point or if the step is
00602        too small to be worthwhile continuing */
00603     while( ( optStatus != MSQ_EQUILIBRIUM ) && ( optStatus != MSQ_STEP_TOO_SMALL ) &&
00604            ( optStatus != MSQ_IMP_TOO_SMALL ) && ( optStatus != MSQ_FLAT_NO_IMP ) && ( optStatus != MSQ_ZERO_SEARCH ) &&
00605            ( optStatus != MSQ_MAX_ITER_EXCEEDED ) )
00606     {
00607 
00608         /* increase the iteration count by one */
00609         /* smooth_param->iter_count += 1; */
00610         iterCount += 1;
00611         optIterCount += 1;
00612         if( iterCount > MSQ_MAX_OPT_ITER ) optStatus = MSQ_MAX_ITER_EXCEEDED;
00613 
00614         MSQ_PRINT( 3 )( "\nITERATION %d \n", iterCount );
00615 
00616         /* compute the gradient */
00617         this->compute_gradient( &pd, mGradient, err );MSQ_ERRRTN( err );
00618 
00619         MSQ_PRINT( 3 )( "Computing the search direction \n" );
00620         this->search_direction( pd, mSearch, err );MSQ_ERRRTN( err );
00621 
00622         /* if there are viable directions to search */
00623         if( ( optStatus != MSQ_ZERO_SEARCH ) && ( optStatus != MSQ_MAX_ITER_EXCEEDED ) )
00624         {
00625 
00626             MSQ_PRINT( 3 )( "Computing the projections of the gradients \n" );
00627             this->get_gradient_projections( mSearch, err );MSQ_ERRRTN( err );
00628 
00629             MSQ_PRINT( 3 )( "Computing the initial step size \n" );
00630             this->compute_alpha( err );MSQ_ERRRTN( err );
00631 
00632             MSQ_PRINT( 3 )( "Testing whether to accept this step \n" );
00633             this->step_acceptance( pd, iterCount, mSearch, err );MSQ_ERRRTN( err );
00634             // MSQ_PRINT(3)("The new free vertex position is %f %f %f\n",
00635             //  mCoords[freeVertexIndex][0],mCoords[freeVertexIndex][1],mCoords[freeVertexIndex][2]);
00636 
00637             if( MSQ_DBG( 3 ) )
00638             {
00639                 /* Print the active set */
00640                 this->print_active_set( mActive, mFunction );
00641             }
00642 
00643             /* check for equilibrium point */
00644             if( mActive.active_ind.size() >= 2 )
00645             {
00646                 MSQ_PRINT( 3 )( "Testing for an equilibrium point \n" );
00647                 equilibriumPt = this->check_equilibrium( optStatus, err );MSQ_ERRRTN( err );
00648 
00649                 if( MSQ_DBG( 2 ) && equilibriumPt ) MSQ_PRINT( 2 )( "Optimization Exiting: An equilibrium point \n" );
00650             }
00651 
00652             /* record the values */
00653             prevActiveValues.push_back( mActive.true_active_value );
00654         }
00655         else
00656         {
00657             /* decrease the iteration count by one */
00658             /* smooth_param->iter_count -= 1; */
00659             iterCount -= 1;
00660             if( MSQ_DBG( 2 ) )
00661             {
00662                 MSQ_PRINT( 2 )
00663                 ( "Optimization Exiting: No viable directions; equilibrium point \n" );
00664                 /* Print the old active set */
00665                 this->print_active_set( mActive, mFunction );
00666             }
00667         }
00668     }
00669 
00670     MSQ_PRINT( 2 )( "Checking the validity of the mesh\n" );
00671     if( !this->validity_check( pd, err ) ) MSQ_PRINT( 2 )( "The final mesh is not valid\n" );MSQ_ERRRTN( err );
00672 
00673     MSQ_PRINT( 2 )( "Number of optimization iterations %d\n", iterCount );
00674 
00675     switch( optStatus )
00676     {
00677         default:
00678             MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Invalid early termination\n" );
00679             break;
00680         case MSQ_EQUILIBRIUM:
00681             MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Equilibrium\n" );
00682             break;
00683         case MSQ_STEP_TOO_SMALL:
00684             MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Step Too Small\n" );
00685             break;
00686         case MSQ_IMP_TOO_SMALL:
00687             MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Improvement Too Small\n" );
00688             break;
00689         case MSQ_FLAT_NO_IMP:
00690             MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Flat No Improvement\n" );
00691             break;
00692         case MSQ_ZERO_SEARCH:
00693             MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Zero Search\n" );
00694             break;
00695         case MSQ_MAX_ITER_EXCEEDED:
00696             MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Max Iter Exceeded\n" );
00697             break;
00698     }
00699 }
00700 
00701 void NonSmoothDescent::step_acceptance( PatchData& pd, int iterCount, const Vector3D& mSearch, MsqError& err )
00702 {
00703     const double minAcceptableImprovement = 1e-6;
00704 
00705     //  int        ierr;
00706     //  int        num_values;
00707     bool step_done;
00708     bool valid = true, accept_alpha;
00709     double estimated_improvement;
00710     double current_improvement  = HUGE_VAL;
00711     double previous_improvement = HUGE_VAL;
00712     double current_percent_diff = HUGE_VAL;
00713     Vector3D original_point;
00714 
00715     //  MSQ_FUNCTION_TIMER( "Step Acceptance" );
00716     //  num_values = qmHandles.size();
00717 
00718     optStatus = MSQ_NO_STATUS;
00719 
00720     if( mAlpha < minStepSize )
00721     {
00722         optStatus = MSQ_IMP_TOO_SMALL;
00723         step_done = true;
00724         MSQ_PRINT( 3 )( "Alpha starts too small, no improvement\n" );
00725     }
00726 
00727     const MsqVertex* coords = pd.get_vertex_array( err );
00728 
00729     /* save the original function and active set */
00730     original_point   = coords[freeVertexIndex];
00731     originalFunction = mFunction;
00732     originalActive   = mActive;
00733 
00734     step_done = false;
00735     for( int num_steps = 0; !step_done && num_steps < 100; ++num_steps )
00736     {
00737         accept_alpha = false;
00738 
00739         while( !accept_alpha && mAlpha > minStepSize )
00740         {
00741 
00742             /* make the step */
00743             pd.move_vertex( -mAlpha * mSearch, freeVertexIndex, err );
00744             // pd.set_coords_array_element(coords[freeVertexIndex],0,err);
00745 
00746             MSQ_PRINT( 2 )( "search direction %f %f \n", mSearch[0], mSearch[1] );
00747             MSQ_PRINT( 2 )
00748             ( "new vertex position %f %f \n", coords[freeVertexIndex][0], coords[freeVertexIndex][1] );
00749 
00750             /* assume alpha is acceptable */
00751             accept_alpha = true;
00752 
00753             /* never take a step that makes a valid mesh invalid or worsens the quality */
00754             // TODO Validity check revision -- do the compute function up here
00755             // and then the rest based on validity
00756             valid = validity_check( pd, err );MSQ_ERRRTN( err );
00757             if( valid )
00758             {
00759                 valid = improvement_check( err );MSQ_ERRRTN( err );
00760             }
00761             if( !valid )
00762             {
00763                 accept_alpha = false;
00764                 pd.move_vertex( mAlpha * mSearch, freeVertexIndex, err );
00765                 // pd.set_coords_array_element(coords[freeVertexIndex],0,err);
00766                 mAlpha = mAlpha / 2;
00767                 MSQ_PRINT( 2 )( "Step not accepted, the new alpha %f\n", mAlpha );
00768 
00769                 if( mAlpha < minStepSize )
00770                 {
00771                     optStatus = MSQ_STEP_TOO_SMALL;
00772                     step_done = true;
00773                     MSQ_PRINT( 2 )( "Step too small\n" );
00774                     /* get back the original point, mFunction, and active set */
00775                     pd.set_vertex_coordinates( original_point, freeVertexIndex, err );
00776                     mFunction = originalFunction;
00777                     mActive   = originalActive;
00778                 }
00779             }
00780         }
00781 
00782         if( valid && ( mAlpha > minStepSize ) )
00783         {
00784             /* compute the new function and active set */
00785             this->compute_function( &pd, mFunction, err );MSQ_ERRRTN( err );
00786             this->find_active_set( mFunction, mActive );
00787 
00788             /* estimate the minimum improvement by taking this step */
00789             this->get_min_estimate( &estimated_improvement, err );MSQ_ERRRTN( err );
00790             MSQ_PRINT( 2 )
00791             ( "The estimated improvement for this step: %f\n", estimated_improvement );
00792 
00793             /* calculate the actual increase */
00794             current_improvement = mActive.true_active_value - prevActiveValues[iterCount - 1];
00795 
00796             MSQ_PRINT( 3 )( "Actual improvement %f\n", current_improvement );
00797 
00798             /* calculate the percent difference from estimated increase */
00799             current_percent_diff = fabs( current_improvement - estimated_improvement ) / fabs( estimated_improvement );
00800 
00801             /* determine whether to accept a step */
00802             if( ( fabs( previous_improvement ) > fabs( current_improvement ) ) && ( previous_improvement < 0 ) )
00803             {
00804                 /* accept the previous step - it was better */
00805                 MSQ_PRINT( 2 )( "Accepting the previous step\n" );
00806 
00807                 /* subtract alpha in again (previous step) */
00808                 pd.move_vertex( -mAlpha * mSearch, freeVertexIndex, err );
00809                 // pd.set_coords_array_element(coords[freeVertexIndex],0,err);
00810 
00811                 /* does this make an invalid mesh valid? */
00812                 // TODO Validity check revisison
00813                 valid = validity_check( pd, err );MSQ_ERRRTN( err );
00814                 if( valid ) valid = improvement_check( err );MSQ_ERRRTN( err );
00815 
00816                 /* copy test function and active set */
00817                 mFunction = testFunction;
00818                 mActive   = testActive;
00819 
00820                 optStatus = MSQ_STEP_ACCEPTED;
00821                 step_done = true;
00822 
00823                 /* check to see that we're still making good improvements */
00824                 if( fabs( previous_improvement ) < minAcceptableImprovement )
00825                 {
00826                     optStatus = MSQ_IMP_TOO_SMALL;
00827                     step_done = true;
00828                     MSQ_PRINT( 2 )( "Optimization Exiting: Improvement too small\n" );
00829                 }
00830             }
00831             else if( ( ( fabs( current_improvement ) > fabs( estimated_improvement ) ) ||
00832                        ( current_percent_diff < .1 ) ) &&
00833                      ( current_improvement < 0 ) )
00834             {
00835                 /* accept this step, exceeded estimated increase or was close */
00836                 optStatus = MSQ_STEP_ACCEPTED;
00837                 step_done = true;
00838 
00839                 /* check to see that we're still making good improvements */
00840                 if( fabs( current_improvement ) < minAcceptableImprovement )
00841                 {
00842                     MSQ_PRINT( 2 )( "Optimization Exiting: Improvement too small\n" );
00843                     optStatus = MSQ_IMP_TOO_SMALL;
00844                     step_done = true;
00845                 }
00846             }
00847             else if( ( current_improvement > 0 ) && ( previous_improvement > 0 ) &&
00848                      ( fabs( current_improvement ) < minAcceptableImprovement ) &&
00849                      ( fabs( previous_improvement ) < minAcceptableImprovement ) )
00850             {
00851 
00852                 /* we are making no progress, quit */
00853                 optStatus = MSQ_FLAT_NO_IMP;
00854                 step_done = true;
00855                 MSQ_PRINT( 2 )( "Opimization Exiting: Flat no improvement\n" );
00856 
00857                 /* get back the original point, function, and active set */
00858                 pd.set_vertex_coordinates( original_point, freeVertexIndex, err );MSQ_ERRRTN( err );
00859                 mFunction = originalFunction;
00860                 mActive   = originalActive;
00861             }
00862             else
00863             {
00864                 /* halve alpha and try again */
00865                 /* add out the old step */
00866                 pd.move_vertex( mAlpha * mSearch, freeVertexIndex, err );
00867                 // pd.set_coords_array_element(coords[freeVertexIndex],0,err);
00868 
00869                 /* halve step size */
00870                 mAlpha = mAlpha / 2;
00871                 MSQ_PRINT( 3 )( "Step not accepted, the new alpha %f\n", mAlpha );
00872 
00873                 if( mAlpha < minStepSize )
00874                 {
00875                     /* get back the original point, function, and active set */
00876                     MSQ_PRINT( 2 )( "Optimization Exiting: Step too small\n" );
00877                     pd.set_vertex_coordinates( original_point, freeVertexIndex, err );MSQ_ERRRTN( err );
00878                     mFunction = originalFunction;
00879                     mActive   = originalActive;
00880                     optStatus = MSQ_STEP_TOO_SMALL;
00881                     step_done = true;
00882                 }
00883                 else
00884                 {
00885                     testFunction         = mFunction;
00886                     testActive           = mActive;
00887                     previous_improvement = current_improvement;
00888                 }
00889             }
00890         }
00891     }
00892     if( current_improvement > 0 && optStatus == MSQ_STEP_ACCEPTED )
00893     { MSQ_PRINT( 2 )( "Accepted a negative step %f \n", current_improvement ); }
00894 }
00895 
00896 bool NonSmoothDescent::compute_function( PatchData* patch_data, std::vector< double >& func, MsqError& err )
00897 {
00898     // NEED 1.0/FUNCTION WHICH IS ONLY TRUE OF CONDITION NUMBER
00899 
00900     func.resize( qmHandles.size() );
00901     bool valid_bool = true;
00902     for( size_t i = 0; i < qmHandles.size(); i++ )
00903     {
00904         valid_bool = valid_bool && currentQM->evaluate( *patch_data, qmHandles[i], func[i], err );
00905         MSQ_ERRZERO( err );
00906     }
00907 
00908     return valid_bool;
00909 }
00910 
00911 bool NonSmoothDescent::compute_gradient( PatchData* patch_data, std::vector< Vector3D >& gradient_out, MsqError& err )
00912 {
00913     MSQ_DBGOUT( 2 ) << "Computing Gradient\n";
00914 
00915     bool valid = true;
00916 
00917 #ifdef NUMERICAL_GRADIENT
00918 
00919     const double delta = 10e-6;
00920     std::vector< double > func( qmHandles.size() ), fdelta( qmHandles.size() );
00921 
00922     valid = this->compute_function( patch_data, func, err );
00923     MSQ_ERRZERO( err );
00924 
00925     /* gradient in the x, y, z direction */
00926 
00927     for( int j = 0; j < 3; j++ )
00928     {
00929 
00930         // perturb the coordinates of the free vertex in the j direction by delta
00931         Vector3D delta_3( 0, 0, 0 );
00932         Vector3D orig_pos = patch_data->vertex_by_index( freeVertexIndex );
00933         delta_3[j]        = delta;
00934         patch_data->move_vertex( delta_3, freeVertexIndex, err );
00935 
00936         // compute the function at the perturbed point location
00937         valid = valid && this->compute_function( patch_data, fdelta, err );
00938         MSQ_ERRZERO( err );
00939 
00940         // compute the numerical gradient
00941         for( size_t i = 0; i < qmHandles.size(); i++ )
00942         {
00943             mGradient[i][j] = ( fdelta[i] - func[i] ) / delta;
00944             // MSQ_DEBUG_ACTION(3,{fprintf(stdout,"  Gradient
00945             // value[%d][%d]=%g\n",i,j,mGradient[i][j]);});
00946         }
00947 
00948         // put the coordinates back where they belong
00949         patch_data->set_vertex_coordinates( orig_pos, freeVertexIndex, err );
00950     }
00951 
00952 #else
00953     double value;
00954     gradient_out.resize( qmHandles.size() );
00955     for( size_t i = 0; i < qmHandles.size(); i++ )
00956     {
00957         valid = valid && currentQM->evaluate_with_gradient( *patch_data, qmHandles[i], value, tmpIdx, tmpGrad, err );
00958         MSQ_ERRZERO( err );
00959         assert( tmpIdx.size() == 1 && tmpIdx[0] == freeVertexIndex );
00960         gradient_out[i] = tmpGrad[0];
00961     }
00962 #endif
00963 
00964     return valid;
00965 }
00966 
00967 void NonSmoothDescent::find_active_set( const std::vector< double >& function, ActiveSet& active_set )
00968 {
00969     // local parameter initialization
00970     const double activeEpsilon = .3e-4;
00971     //  activeEpsilon = .3e-8;
00972 
00973     double function_val;
00974     double active_value0;
00975     double temp;
00976 
00977     //    FUNCTION_TIMER_START("Find Active Set");
00978     MSQ_DBGOUT( 2 ) << "\nFinding the active set\n";
00979 
00980     /* the first function value is our initial active value */
00981     active_set.set( 0 );
00982     active_set.true_active_value = function[0];
00983     //    MSQ_DEBUG_ACTION(3,{fprintf(stdout,"  function value[0]: %g\n",function[0]);});
00984 
00985     /* first sort out the active set...
00986        all vals within active_epsilon of largest val */
00987 
00988     for( size_t i = 1; i < qmHandles.size(); i++ )
00989     {
00990         function_val = function[i];
00991         if( active_set.true_active_value < function_val ) active_set.true_active_value = function_val;
00992         active_value0 = function[active_set.active_ind[0]];
00993         temp          = fabs( function_val - active_value0 );
00994         //        MSQ_DEBUG_ACTION(3,{fprintf(stdout,"  function value[%d]: %g\n",i,function[i]);});
00995         if( function_val > active_value0 )
00996         {  // seek max_i function[i]
00997             if( temp >= activeEpsilon )
00998             {
00999                 active_set.set( i );  // new max
01000             }
01001             else
01002             {
01003                 active_set.add( i, fabs( function_val - active_value0 ) < EPSILON );
01004             }
01005         }
01006         else
01007         {
01008             if( temp < activeEpsilon ) { active_set.add( i, fabs( function_val - active_value0 ) < EPSILON ); }
01009         }
01010     }
01011 }
01012 
01013 bool NonSmoothDescent::validity_check( PatchData& pd, MsqError& err )
01014 
01015 {
01016     //  FUNCTION_TIMER_START("Validity Check");
01017 
01018     // ONLY FOR SIMPLICIAL MESHES - THERE SHOULD BE A VALIDITY CHECKER ASSOCIATED
01019     // WITH MSQ ELEMENTS
01020 
01021     /* check that the simplicial mesh is still valid, based on right handedness.
01022          Returns a 1 or a 0 */
01023 
01024     // TODO as a first step we can switch this over to the function
01025     // evaluation and use the rest of the code as is
01026     bool valid  = true;
01027     double dEps = 1.e-13;
01028 
01029     double x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4;
01030     const MsqVertex* coords = pd.get_vertex_array( err );
01031 
01032     for( size_t i = 0; i < pd.num_elements(); i++ )
01033     {
01034         const size_t* conn = pd.element_by_index( i ).get_vertex_index_array();
01035         coords[conn[0]].get_coordinates( x1, y1, z1 );
01036         coords[conn[1]].get_coordinates( x2, y2, z2 );
01037         coords[conn[2]].get_coordinates( x3, y3, z3 );
01038         coords[conn[3]].get_coordinates( x4, y4, z4 );
01039 
01040         double dDX2 = x2 - x1;
01041         double dDX3 = x3 - x1;
01042         double dDX4 = x4 - x1;
01043 
01044         double dDY2 = y2 - y1;
01045         double dDY3 = y3 - y1;
01046         double dDY4 = y4 - y1;
01047 
01048         double dDZ2 = z2 - z1;
01049         double dDZ3 = z3 - z1;
01050         double dDZ4 = z4 - z1;
01051 
01052         /* dDet is proportional to the cell volume */
01053         double dDet = dDX2 * dDY3 * dDZ4 + dDX3 * dDY4 * dDZ2 + dDX4 * dDY2 * dDZ3 - dDZ2 * dDY3 * dDX4 -
01054                       dDZ3 * dDY4 * dDX2 - dDZ4 * dDY2 * dDX3;
01055 
01056         /* Compute a length scale based on edge lengths. */
01057         double dScale = ( sqrt( ( x1 - x2 ) * ( x1 - x2 ) + ( y1 - y2 ) * ( y1 - y2 ) + ( z1 - z2 ) * ( z1 - z2 ) ) +
01058                           sqrt( ( x1 - x3 ) * ( x1 - x3 ) + ( y1 - y3 ) * ( y1 - y3 ) + ( z1 - z3 ) * ( z1 - z3 ) ) +
01059                           sqrt( ( x1 - x4 ) * ( x1 - x4 ) + ( y1 - y4 ) * ( y1 - y4 ) + ( z1 - z4 ) * ( z1 - z4 ) ) +
01060                           sqrt( ( x2 - x3 ) * ( x2 - x3 ) + ( y2 - y3 ) * ( y2 - y3 ) + ( z2 - z3 ) * ( z2 - z3 ) ) +
01061                           sqrt( ( x2 - x4 ) * ( x2 - x4 ) + ( y2 - y4 ) * ( y2 - y4 ) + ( z2 - z4 ) * ( z2 - z4 ) ) +
01062                           sqrt( ( x3 - x4 ) * ( x3 - x4 ) + ( y3 - y4 ) * ( y3 - y4 ) + ( z3 - z4 ) * ( z3 - z4 ) ) ) /
01063                         6.;
01064 
01065         /* Use the length scale to get a better idea if the tet is flat or
01066            just really small. */
01067         if( fabs( dScale ) < EPSILON ) { return false; }
01068         else
01069         {
01070             dDet /= ( dScale * dScale * dScale );
01071         }
01072 
01073         if( dDet > dEps ) { valid = true; }
01074         else if( dDet < -dEps )
01075         {
01076             return false;
01077         }
01078         else
01079         {
01080             return false;
01081         }
01082     }  // end for i=1,numElements
01083 
01084     //  MSQ_DEBUG_ACTION(2,{fprintf(stdout,"Mesh Validity is: %d \n",valid);});
01085 
01086     //  FUNCTION_TIMER_END();
01087     return ( valid );
01088 }
01089 
01090 void NonSmoothDescent::get_active_directions( const std::vector< Vector3D >& pGradient, std::vector< Vector3D >& dir,
01091                                               MsqError& /*err*/ )
01092 {
01093     const size_t num_active = mActive.active_ind.size();
01094     dir.resize( num_active );
01095     for( size_t i = 0; i < num_active; i++ )
01096     {
01097         dir[i] = pGradient[mActive.active_ind[i]];
01098     }
01099 }
01100 
01101 bool NonSmoothDescent::check_vector_dots( const std::vector< Vector3D >& vec, const Vector3D& normal,
01102                                           MsqError& /*err*/ )
01103 {
01104     double test_dot = vec[0] % normal;
01105     unsigned ind    = 1;
01106     while( ( fabs( test_dot ) < EPSILON ) && ( ind < vec.size() ) )
01107     {
01108         test_dot = vec[ind] % normal;
01109         ind++;
01110     }
01111 
01112     for( unsigned i = ind; i < vec.size(); i++ )
01113     {
01114         double dot = vec[i] % normal;
01115         if( ( ( dot > 0 && test_dot < 0 ) || ( dot < 0 && test_dot > 0 ) ) && ( fabs( dot ) > EPSILON ) )
01116         { return true; }
01117     }
01118     return false;
01119 }
01120 
01121 bool NonSmoothDescent::convex_hull_test( const std::vector< Vector3D >& vec, MsqError& err )
01122 {
01123     //    int ierr;
01124     bool equil         = false;
01125     Direction dir_done = MSQ_ZDIR;
01126     Status status      = MSQ_CHECK_Z_COORD_DIRECTION;
01127     Vector3D pt1, pt2, pt3, normal;
01128 
01129     /* tries to determine equilibrium for the 3D case */
01130 
01131     if( vec.size() <= 2 ) status = MSQ_NO_EQUIL;
01132 
01133     while( ( status != MSQ_EQUIL ) && ( status != MSQ_NO_EQUIL ) && ( status != MSQ_HULL_TEST_ERROR ) )
01134     {
01135         if( status == MSQ_CHECK_Z_COORD_DIRECTION )
01136         {
01137             this->find_plane_points( MSQ_ZDIR, MSQ_YDIR, vec, pt1, pt2, pt3, status, err );
01138             dir_done = MSQ_ZDIR;
01139         }
01140         if( status == MSQ_CHECK_Y_COORD_DIRECTION )
01141         {
01142             this->find_plane_points( MSQ_YDIR, MSQ_XDIR, vec, pt1, pt2, pt3, status, err );
01143             dir_done = MSQ_YDIR;
01144         }
01145         if( status == MSQ_CHECK_X_COORD_DIRECTION )
01146         {
01147             this->find_plane_points( MSQ_XDIR, MSQ_ZDIR, vec, pt1, pt2, pt3, status, err );
01148             dir_done = MSQ_XDIR;
01149         }
01150         if( status == MSQ_TWO_PT_PLANE ) { pt3 = Vector3D( 0, 0, 0 ); }
01151         if( ( status == MSQ_TWO_PT_PLANE ) || ( status == MSQ_THREE_PT_PLANE ) )
01152         {
01153             this->find_plane_normal( pt1, pt2, pt3, normal, err );
01154             equil = this->check_vector_dots( vec, normal, err );
01155             if( equil )
01156             {
01157                 switch( dir_done )
01158                 {
01159                     case MSQ_ZDIR:
01160                         equil  = false;
01161                         status = MSQ_CHECK_Y_COORD_DIRECTION;
01162                         break;
01163                     case MSQ_YDIR:
01164                         equil  = false;
01165                         status = MSQ_CHECK_X_COORD_DIRECTION;
01166                         break;
01167                     case MSQ_XDIR:
01168                         equil  = true;
01169                         status = MSQ_EQUIL;
01170                 }
01171             }
01172             else if( equil == 0 )
01173             {
01174                 status = MSQ_NO_EQUIL;
01175             }
01176             else
01177             {
01178                 MSQ_SETERR( err )( "equil flag not set to 0 or 1", MsqError::INVALID_STATE );
01179             }
01180         }
01181     }
01182     switch( status )
01183     {
01184         case MSQ_NO_EQUIL:
01185             MSQ_PRINT( 3 )( "Not an equilibrium point\n" );
01186             equil = false;
01187             break;
01188         case MSQ_EQUIL:
01189             MSQ_PRINT( 3 )( "An equilibrium point\n" );
01190             equil = true;
01191             break;
01192         default:
01193             MSQ_PRINT( 3 )( "Failed to determine equil or not; status = %d\n", status );
01194     }
01195     //    FUNCTION_TIMER_END();
01196     return ( equil );
01197 }
01198 
01199 void NonSmoothDescent::form_grammian( const std::vector< Vector3D >& vec, MsqError& /*err*/ )
01200 {
01201     /* form the grammian with the dot products of the gradients */
01202     const size_t num_active = mActive.active_ind.size();
01203     mG.resize( num_active );
01204     for( size_t i = 0; i < num_active; i++ )
01205         for( size_t j = i; j < num_active; j++ )
01206             mG( i, j ) = vec[i] % vec[j];
01207 }
01208 
01209 bool NonSmoothDescent::check_equilibrium( OptStatus& status, MsqError& err )
01210 {
01211     std::vector< Vector3D > dir;
01212 
01213     // TODO - this subroutine is no longer clear to me... is it still
01214     // appropriate for quads and hexes?  I think it might be in 2D, but
01215     // 3D is less clear.  Is there a more general algorithm to use?
01216     // ask Todd/check in numerical optimization
01217 
01218     bool equil              = false;
01219     const size_t num_active = mActive.active_ind.size();
01220     ;
01221 
01222     if( num_active == qmHandles.size() )
01223     {
01224         equil  = true;
01225         status = MSQ_EQUILIBRIUM;
01226         MSQ_PRINT( 3 )( "All the function values are in the active set\n" );
01227     }
01228 
01229     /* set up the active mGradient directions */
01230     this->get_active_directions( mGradient, dir, err );
01231     MSQ_ERRZERO( err );
01232 
01233     /* normalize the active directions */
01234     for( size_t j = 0; j < num_active; j++ )
01235         dir[j] /= dir[j].length();
01236 
01237     equil = this->convex_hull_test( dir, err );
01238     if( equil ) status = MSQ_EQUILIBRIUM;
01239 
01240     return equil;
01241 }
01242 
01243 static double condition3x3( const double A[9] )
01244 {
01245     //   int ierr;
01246     double a11, a12, a13;
01247     double a21, a22, a23;
01248     double a31, a32, a33;
01249     //   double s1, s2, s4, s3, t0;
01250     double s1, s2, s3;
01251     double denom;
01252     //   double one = 1.0;
01253     double temp;
01254     bool zero_denom = true;
01255 
01256     a11 = A[0];
01257     a12 = A[1];
01258     a13 = A[2];
01259     a21 = A[3];
01260     a22 = A[4];
01261     a23 = A[5];
01262     a31 = A[6];
01263     a32 = A[7];
01264     a33 = A[8];
01265 
01266     denom = -a11 * a22 * a33 + a11 * a23 * a32 + a21 * a12 * a33 - a21 * a13 * a32 - a31 * a12 * a23 + a31 * a13 * a22;
01267 
01268     if( ( fabs( a11 ) > EPSILON ) && ( fabs( denom / a11 ) > EPSILON ) ) { zero_denom = false; }
01269     if( ( fabs( a22 ) > EPSILON ) && ( fabs( denom / a22 ) > EPSILON ) ) { zero_denom = false; }
01270     if( ( fabs( a33 ) > EPSILON ) && ( fabs( denom / a33 ) > EPSILON ) ) { zero_denom = false; }
01271 
01272     if( zero_denom ) { return HUGE_VAL; }
01273     else
01274     {
01275         s1 = sqrt( a11 * a11 + a12 * a12 + a13 * a13 + a21 * a21 + a22 * a22 + a23 * a23 + a31 * a31 + a32 * a32 +
01276                    a33 * a33 );
01277 
01278         temp = ( -a22 * a33 + a23 * a32 ) / denom;
01279         s3   = temp * temp;
01280         temp = ( a12 * a33 - a13 * a32 ) / denom;
01281         s3 += temp * temp;
01282         temp = ( a12 * a23 - a13 * a22 ) / denom;
01283         s3 += temp * temp;
01284         temp = ( a21 * a33 - a23 * a31 ) / denom;
01285         s3 += temp * temp;
01286         temp = ( a11 * a33 - a13 * a31 ) / denom;
01287         s3 += temp * temp;
01288         temp = ( a11 * a23 - a13 * a21 ) / denom;
01289         s3 += temp * temp;
01290         temp = ( a21 * a32 - a22 * a31 ) / denom;
01291         s3 += temp * temp;
01292         temp = ( -a11 * a32 + a12 * a31 ) / denom;
01293         s3 += temp * temp;
01294         temp = ( -a11 * a22 + a12 * a21 ) / denom;
01295         s3 += temp * temp;
01296 
01297         s2 = sqrt( s3 );
01298         return s1 * s2;
01299     }
01300 }
01301 
01302 double NonSmoothDescent::SymmetricMatrix::condition3x3() const
01303 {
01304     double values[9] = { operator()( 0, 0 ), operator()( 0, 1 ), operator()( 0, 2 ),
01305                          operator()( 1, 0 ), operator()( 1, 1 ), operator()( 1, 2 ),
01306                          operator()( 2, 0 ), operator()( 2, 1 ), operator()( 2, 2 ) };
01307     return MBMesquite::condition3x3( values );
01308 }
01309 
01310 void NonSmoothDescent::singular_test( int n, const Matrix3D& A, bool& singular, MsqError& err )
01311 {
01312     //    int test;
01313     //    double determinant;
01314     double cond;
01315 
01316     if( ( n > 3 ) || ( n < 1 ) )
01317     {
01318         MSQ_SETERR( err )( "Singular test works only for n=1 to n=3", MsqError::INVALID_ARG );
01319         return;
01320     }
01321 
01322     singular = true;
01323     switch( n )
01324     {
01325         case 1:
01326             if( A[0][0] > 0 ) singular = false;
01327             break;
01328         case 2:
01329             if( fabs( A[0][0] * A[1][1] - A[0][1] * A[1][0] ) > EPSILON ) singular = false;
01330             break;
01331         case 3:
01332             /* calculate the condition number */
01333             cond = condition3x3( A[0] );
01334             if( cond < 1E14 ) singular = false;
01335             break;
01336     }
01337 }
01338 
01339 void NonSmoothDescent::form_PD_grammian( MsqError& err )
01340 {
01341     int i, j, k;
01342     int g_ind_1;
01343     bool singular = false;
01344 
01345     const int num_active = mActive.active_ind.size();
01346 
01347     /* this assumes the grammian has been formed */
01348     for( i = 0; i < num_active; i++ )
01349     {
01350         for( j = i; j < num_active; j++ )
01351         {
01352             if( mG( i, j ) == -1 )
01353             {
01354                 MSQ_SETERR( err )( "Grammian not computed properly", MsqError::INVALID_STATE );
01355                 return;
01356             }
01357         }
01358     }
01359 
01360     /* use the first gradient in the active set */
01361     g_ind_1    = 0;
01362     mPDG[0][0] = mG( 0, 0 );
01363     pdgInd[0]  = mActive.active_ind[0];
01364 
01365     /* test the rest and add them as appropriate */
01366     k = 1;
01367     i = 1;
01368     while( ( k < 3 ) && ( i < num_active ) )
01369     {
01370         mPDG[0][k] = mPDG[k][0] = mG( 0, i );
01371         mPDG[k][k]              = mG( i, i );
01372         if( k == 2 )
01373         { /* add the dot product of g1 and g2 */
01374             mPDG[1][k] = mPDG[k][1] = mG( g_ind_1, i );
01375         }
01376         this->singular_test( k + 1, mPDG, singular, err );
01377         if( !singular )
01378         {
01379             pdgInd[k] = mActive.active_ind[i];
01380             if( k == 1 ) g_ind_1 = i;
01381             k++;
01382         }
01383         i++;
01384     }
01385 }
01386 
01387 void NonSmoothDescent::search_edges_faces( const Vector3D* dir, Vector3D& mSearch, MsqError& /*err*/ )
01388 {
01389     bool viable;
01390     double a, b, denom;
01391     Vector3D g_bar;
01392     Vector3D temp_search( 0, 0, 0 ); /* initialize the search direction to 0,0 */
01393     double projection, min_projection;
01394 
01395     const size_t num_active = mActive.active_ind.size();
01396     ;
01397 
01398     /* Check for viable faces */
01399     min_projection = HUGE_VAL;
01400     for( size_t i = 0; i < num_active; i++ )
01401     {
01402         /* FACE I */
01403         viable = true;
01404 
01405         /* test the viability */
01406         for( size_t j = 0; j < num_active; j++ )
01407         { /* lagrange multipliers>0 */
01408             if( mG( j, i ) < 0 ) viable = false;
01409         }
01410 
01411         /* find the minimum of viable directions */
01412         if( ( viable ) && ( mG( i, i ) < min_projection ) )
01413         {
01414             min_projection = mG( i, i );
01415             temp_search    = dir[i];
01416             mSteepest      = mActive.active_ind[i];
01417         }
01418 
01419         /* INTERSECTION IJ */
01420         for( size_t j = i + 1; j < num_active; j++ )
01421         {
01422             viable = true;
01423 
01424             /* find the coefficients of the intersection
01425                and test the viability */
01426             denom = 2 * mG( i, j ) - mG( i, i ) - mG( j, j );
01427             a = b = 0;
01428             if( fabs( denom ) > EPSILON )
01429             {
01430                 b = ( mG( i, j ) - mG( i, i ) ) / denom;
01431                 a = 1 - b;
01432                 if( ( b < 0 ) || ( b > 1 ) ) /* 0 < b < 1 */
01433                     viable = false;
01434                 for( size_t k = 0; k < num_active; k++ )
01435                 { /* lagrange multipliers>0 */
01436                     if( ( a * mG( k, i ) + b * mG( k, j ) ) <= 0 ) viable = false;
01437                 }
01438             }
01439             else
01440             {
01441                 viable = false; /* Linearly dependent */
01442             }
01443 
01444             /* find the minimum of viable directions */
01445             if( viable )
01446             {
01447                 g_bar      = a * dir[i] + b * dir[j];
01448                 projection = g_bar % g_bar;
01449                 if( projection < min_projection )
01450                 {
01451                     min_projection = projection;
01452                     temp_search    = g_bar;
01453                     mSteepest      = mActive.active_ind[i];
01454                 }
01455             }
01456         }
01457     }
01458     if( optStatus != MSQ_EQUILIBRIUM ) { mSearch = temp_search; }
01459 }
01460 
01461 bool NonSmoothDescent::solve2x2( double a11, double a12, double a21, double a22, double b1, double b2, double x[2],
01462                                  MsqError& /*err*/ )
01463 {
01464     double factor;
01465 
01466     /* if the system is not singular, solve it */
01467     if( fabs( a11 * a22 - a21 * a12 ) > EPSILON )
01468     {
01469         if( fabs( a11 ) > EPSILON )
01470         {
01471             factor = ( a21 / a11 );
01472             x[1]   = ( b2 - factor * b1 ) / ( a22 - factor * a12 );
01473             x[0]   = ( b1 - a12 * x[1] ) / a11;
01474         }
01475         else if( fabs( a21 ) > EPSILON )
01476         {
01477             factor = ( a11 / a21 );
01478             x[1]   = ( b1 - factor * b2 ) / ( a12 - factor * a22 );
01479             x[0]   = ( b2 - a22 * x[1] ) / a21;
01480         }
01481         return true;
01482     }
01483     else
01484     {
01485         return false;
01486     }
01487 }
01488 
01489 void NonSmoothDescent::form_reduced_matrix( SymmetricMatrix& P )
01490 {
01491     const size_t P_size = mActive.active_ind.size() - 1;
01492     P.resize( P_size );
01493     for( size_t i = 0; i < P_size; i++ )
01494     {
01495         P( i, i ) = mG( 0, 0 ) - 2 * mG( 0, i + 1 ) + mG( i + 1, i + 1 );
01496         for( size_t j = i + 1; j < P_size; j++ )
01497         {
01498             P( i, j ) = mG( 0, 0 ) - mG( 0, j + 1 ) - mG( i + 1, 0 ) + mG( i + 1, j + 1 );
01499         }
01500     }
01501 }
01502 
01503 void NonSmoothDescent::get_min_estimate( double* final_est, MsqError& /*err*/ )
01504 {
01505     double est_imp;
01506 
01507     *final_est = -HUGE_VAL;
01508     for( size_t i = 0; i < mActive.active_ind.size(); i++ )
01509     {
01510         est_imp = -mAlpha * mGS[mActive.active_ind[i]];
01511         if( est_imp > *final_est ) *final_est = est_imp;
01512     }
01513     if( *final_est == 0 )
01514     {
01515         *final_est = -HUGE_VAL;
01516         for( size_t i = 0; i < qmHandles.size(); i++ )
01517         {
01518             est_imp = -mAlpha * mGS[i];
01519             if( ( est_imp > *final_est ) && ( fabs( est_imp ) > EPSILON ) ) { *final_est = est_imp; }
01520         }
01521     }
01522 }
01523 
01524 void NonSmoothDescent::get_gradient_projections( const Vector3D& mSearch, MsqError& /*err*/ )
01525 {
01526     for( size_t i = 0; i < qmHandles.size(); i++ )
01527         mGS[i] = mGradient[i] % mSearch;
01528 
01529     MSQ_PRINT( 3 )( "steepest in get_gradient_projections %lu\n", (unsigned long)mSteepest );
01530 }
01531 
01532 void NonSmoothDescent::compute_alpha( MsqError& /*err*/ )
01533 {
01534     double steepest_function;
01535     double steepest_grad;
01536     double alpha_i;
01537     double min_positive_value = HUGE_VAL;
01538 
01539     //    FUNCTION_TIMER_START("Compute Alpha");
01540 
01541     MSQ_PRINT( 2 )( "In compute alpha\n" );
01542 
01543     mAlpha = HUGE_VAL;
01544 
01545     steepest_function = mFunction[mSteepest];
01546     steepest_grad     = mGS[mSteepest];
01547     for( size_t i = 0; i < qmHandles.size(); i++ )
01548     {
01549         /* if it's not active */
01550         if( i != mSteepest )
01551         {
01552             alpha_i = steepest_function - mFunction[i];
01553 
01554             if( fabs( mGS[mSteepest] - mGS[i] ) > 1E-13 )
01555             {
01556                 /* compute line intersection */
01557                 alpha_i = alpha_i / ( steepest_grad - mGS[i] );
01558             }
01559             else
01560             {
01561                 /* the lines don't intersect - it's not under consideration*/
01562                 alpha_i = 0;
01563             }
01564             if( ( alpha_i > minStepSize ) && ( fabs( alpha_i ) < fabs( mAlpha ) ) )
01565             {
01566                 mAlpha = fabs( alpha_i );
01567                 MSQ_PRINT( 3 )( "Setting alpha %lu %g\n", (unsigned long)i, alpha_i );
01568             }
01569             if( ( alpha_i > 0 ) && ( alpha_i < min_positive_value ) ) { min_positive_value = alpha_i; }
01570         }
01571     }
01572 
01573     if( ( mAlpha == HUGE_VAL ) && ( min_positive_value != HUGE_VAL ) ) { mAlpha = min_positive_value; }
01574 
01575     /* if it never gets set, set it to the default */
01576     if( mAlpha == HUGE_VAL )
01577     {
01578         mAlpha = maxAlpha;
01579         MSQ_PRINT( 3 )( "Setting alpha to the maximum step length\n" );
01580     }
01581 
01582     MSQ_PRINT( 3 )( "  The initial step size: %f\n", mAlpha );
01583 
01584     //    FUNCTION_TIMER_END();
01585 }
01586 
01587 void NonSmoothDescent::print_active_set( const ActiveSet& active_set, const std::vector< double >& func )
01588 {
01589     if( active_set.active_ind.empty() ) MSQ_DBGOUT( 3 ) << "No active values\n";
01590     /* print the active set */
01591     for( size_t i = 0; i < active_set.active_ind.size(); i++ )
01592     {
01593         MSQ_PRINT( 3 )
01594         ( "Active value %lu:   %f \n", (unsigned long)i + 1, func[active_set.active_ind[i]] );
01595     }
01596 }
01597 
01598 void NonSmoothDescent::init_opt( PatchData& pd, MsqError& err )
01599 {
01600     qmHandles.clear();
01601     currentQM->get_evaluations( pd, qmHandles, true, err );MSQ_ERRRTN( err );
01602 
01603     MSQ_PRINT( 2 )( "\nInitializing Optimization \n" );
01604 
01605     /* for the purposes of initialization will be set to zero after */
01606     optStatus = MSQ_NO_STATUS;
01607     mSteepest = 0;
01608     mAlpha    = 0;
01609     maxAlpha  = 0;
01610 
01611     MSQ_PRINT( 3 )( "  Initialized Constants \n" );
01612     pdgInd[0] = pdgInd[1] = pdgInd[2] = -1;
01613     mPDG                              = Matrix3D( 0.0 );
01614 
01615     MSQ_PRINT( 3 )( "  Initialized search and PDG \n" );
01616     mFunction.clear();
01617     mFunction.resize( qmHandles.size(), 0.0 );
01618     testFunction.clear();
01619     testFunction.resize( qmHandles.size(), 0.0 );
01620     originalFunction.clear();
01621     originalFunction.resize( qmHandles.size(), 0.0 );
01622     mGradient.clear();
01623     mGradient.resize( qmHandles.size(), Vector3D( 0.0 ) );
01624     mGS.resize( qmHandles.size() );
01625 
01626     MSQ_PRINT( 3 )( "  Initialized function/gradient \n" );
01627     mG.resize( qmHandles.size() );
01628     mG.fill( -1 );
01629     MSQ_PRINT( 3 )( "  Initialized G\n" );
01630 
01631     prevActiveValues.clear();
01632     prevActiveValues.reserve( 32 );
01633     MSQ_PRINT( 3 )( "  Initialized prevActiveValues\n" );
01634 }
01635 
01636 void NonSmoothDescent::init_max_step_length( PatchData& pd, MsqError& err )
01637 {
01638     size_t i, j;
01639     double max_diff = 0;
01640     double diff     = 0;
01641 
01642     MSQ_PRINT( 2 )( "In init_max_step_length\n" );
01643 
01644     /* check that the input data is correct */
01645     if( pd.num_elements() == 0 )
01646     {
01647         MSQ_SETERR( err )( "Num incident vtx = 0\n", MsqError::INVALID_MESH );
01648         return;
01649     }
01650 
01651     /* find the maximum distance between two incident vertex locations */
01652     const MsqVertex* coords = pd.get_vertex_array( err );
01653     for( i = 0; i < pd.num_nodes() - 1; i++ )
01654     {
01655         for( j = i; j < pd.num_nodes(); j++ )
01656         {
01657             diff = ( coords[i] - coords[j] ).length_squared();
01658             if( max_diff < diff ) max_diff = diff;
01659         }
01660     }
01661     max_diff = sqrt( max_diff );
01662     if( max_diff == 0 )
01663     {
01664         MSQ_SETERR( err )
01665         ( "Maximum distance between incident vertices = 0\n", MsqError::INVALID_MESH );
01666         return;
01667     }
01668     maxAlpha = max_diff / 100;
01669 
01670     MSQ_PRINT( 3 )( "  Maximum step is %g\n", maxAlpha );
01671 }
01672 
01673 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines