MOAB: Mesh Oriented datABase  (version 5.3.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     {
00894         MSQ_PRINT( 2 )( "Accepted a negative step %f \n", current_improvement );
00895     }
00896 }
00897 
00898 bool NonSmoothDescent::compute_function( PatchData* patch_data, std::vector< double >& func, MsqError& err )
00899 {
00900     // NEED 1.0/FUNCTION WHICH IS ONLY TRUE OF CONDITION NUMBER
00901 
00902     func.resize( qmHandles.size() );
00903     bool valid_bool = true;
00904     for( size_t i = 0; i < qmHandles.size(); i++ )
00905     {
00906         valid_bool = valid_bool && currentQM->evaluate( *patch_data, qmHandles[i], func[i], err );
00907         MSQ_ERRZERO( err );
00908     }
00909 
00910     return valid_bool;
00911 }
00912 
00913 bool NonSmoothDescent::compute_gradient( PatchData* patch_data, std::vector< Vector3D >& gradient_out, MsqError& err )
00914 {
00915     MSQ_DBGOUT( 2 ) << "Computing Gradient\n";
00916 
00917     bool valid = true;
00918 
00919 #ifdef NUMERICAL_GRADIENT
00920 
00921     const double delta = 10e-6;
00922     std::vector< double > func( qmHandles.size() ), fdelta( qmHandles.size() );
00923 
00924     valid = this->compute_function( patch_data, func, err );
00925     MSQ_ERRZERO( err );
00926 
00927     /* gradient in the x, y, z direction */
00928 
00929     for( int j = 0; j < 3; j++ )
00930     {
00931 
00932         // perturb the coordinates of the free vertex in the j direction by delta
00933         Vector3D delta_3( 0, 0, 0 );
00934         Vector3D orig_pos = patch_data->vertex_by_index( freeVertexIndex );
00935         delta_3[j]        = delta;
00936         patch_data->move_vertex( delta_3, freeVertexIndex, err );
00937 
00938         // compute the function at the perturbed point location
00939         valid = valid && this->compute_function( patch_data, fdelta, err );
00940         MSQ_ERRZERO( err );
00941 
00942         // compute the numerical gradient
00943         for( size_t i = 0; i < qmHandles.size(); i++ )
00944         {
00945             mGradient[i][j] = ( fdelta[i] - func[i] ) / delta;
00946             // MSQ_DEBUG_ACTION(3,{fprintf(stdout,"  Gradient
00947             // value[%d][%d]=%g\n",i,j,mGradient[i][j]);});
00948         }
00949 
00950         // put the coordinates back where they belong
00951         patch_data->set_vertex_coordinates( orig_pos, freeVertexIndex, err );
00952     }
00953 
00954 #else
00955     double value;
00956     gradient_out.resize( qmHandles.size() );
00957     for( size_t i = 0; i < qmHandles.size(); i++ )
00958     {
00959         valid = valid && currentQM->evaluate_with_gradient( *patch_data, qmHandles[i], value, tmpIdx, tmpGrad, err );
00960         MSQ_ERRZERO( err );
00961         assert( tmpIdx.size() == 1 && tmpIdx[0] == freeVertexIndex );
00962         gradient_out[i] = tmpGrad[0];
00963     }
00964 #endif
00965 
00966     return valid;
00967 }
00968 
00969 void NonSmoothDescent::find_active_set( const std::vector< double >& function, ActiveSet& active_set )
00970 {
00971     // local parameter initialization
00972     const double activeEpsilon = .3e-4;
00973     //  activeEpsilon = .3e-8;
00974 
00975     double function_val;
00976     double active_value0;
00977     double temp;
00978 
00979     //    FUNCTION_TIMER_START("Find Active Set");
00980     MSQ_DBGOUT( 2 ) << "\nFinding the active set\n";
00981 
00982     /* the first function value is our initial active value */
00983     active_set.set( 0 );
00984     active_set.true_active_value = function[0];
00985     //    MSQ_DEBUG_ACTION(3,{fprintf(stdout,"  function value[0]: %g\n",function[0]);});
00986 
00987     /* first sort out the active set...
00988        all vals within active_epsilon of largest val */
00989 
00990     for( size_t i = 1; i < qmHandles.size(); i++ )
00991     {
00992         function_val = function[i];
00993         if( active_set.true_active_value < function_val ) active_set.true_active_value = function_val;
00994         active_value0 = function[active_set.active_ind[0]];
00995         temp          = fabs( function_val - active_value0 );
00996         //        MSQ_DEBUG_ACTION(3,{fprintf(stdout,"  function value[%d]: %g\n",i,function[i]);});
00997         if( function_val > active_value0 )
00998         {  // seek max_i function[i]
00999             if( temp >= activeEpsilon )
01000             {
01001                 active_set.set( i );  // new max
01002             }
01003             else
01004             {
01005                 active_set.add( i, fabs( function_val - active_value0 ) < EPSILON );
01006             }
01007         }
01008         else
01009         {
01010             if( temp < activeEpsilon ) { active_set.add( i, fabs( function_val - active_value0 ) < EPSILON ); }
01011         }
01012     }
01013 }
01014 
01015 bool NonSmoothDescent::validity_check( PatchData& pd, MsqError& err )
01016 
01017 {
01018     //  FUNCTION_TIMER_START("Validity Check");
01019 
01020     // ONLY FOR SIMPLICIAL MESHES - THERE SHOULD BE A VALIDITY CHECKER ASSOCIATED
01021     // WITH MSQ ELEMENTS
01022 
01023     /* check that the simplicial mesh is still valid, based on right handedness.
01024          Returns a 1 or a 0 */
01025 
01026     // TODO as a first step we can switch this over to the function
01027     // evaluation and use the rest of the code as is
01028     bool valid  = true;
01029     double dEps = 1.e-13;
01030 
01031     double x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4;
01032     const MsqVertex* coords = pd.get_vertex_array( err );
01033 
01034     for( size_t i = 0; i < pd.num_elements(); i++ )
01035     {
01036         const size_t* conn = pd.element_by_index( i ).get_vertex_index_array();
01037         coords[conn[0]].get_coordinates( x1, y1, z1 );
01038         coords[conn[1]].get_coordinates( x2, y2, z2 );
01039         coords[conn[2]].get_coordinates( x3, y3, z3 );
01040         coords[conn[3]].get_coordinates( x4, y4, z4 );
01041 
01042         double dDX2 = x2 - x1;
01043         double dDX3 = x3 - x1;
01044         double dDX4 = x4 - x1;
01045 
01046         double dDY2 = y2 - y1;
01047         double dDY3 = y3 - y1;
01048         double dDY4 = y4 - y1;
01049 
01050         double dDZ2 = z2 - z1;
01051         double dDZ3 = z3 - z1;
01052         double dDZ4 = z4 - z1;
01053 
01054         /* dDet is proportional to the cell volume */
01055         double dDet = dDX2 * dDY3 * dDZ4 + dDX3 * dDY4 * dDZ2 + dDX4 * dDY2 * dDZ3 - dDZ2 * dDY3 * dDX4 -
01056                       dDZ3 * dDY4 * dDX2 - dDZ4 * dDY2 * dDX3;
01057 
01058         /* Compute a length scale based on edge lengths. */
01059         double dScale = ( sqrt( ( x1 - x2 ) * ( x1 - x2 ) + ( y1 - y2 ) * ( y1 - y2 ) + ( z1 - z2 ) * ( z1 - z2 ) ) +
01060                           sqrt( ( x1 - x3 ) * ( x1 - x3 ) + ( y1 - y3 ) * ( y1 - y3 ) + ( z1 - z3 ) * ( z1 - z3 ) ) +
01061                           sqrt( ( x1 - x4 ) * ( x1 - x4 ) + ( y1 - y4 ) * ( y1 - y4 ) + ( z1 - z4 ) * ( z1 - z4 ) ) +
01062                           sqrt( ( x2 - x3 ) * ( x2 - x3 ) + ( y2 - y3 ) * ( y2 - y3 ) + ( z2 - z3 ) * ( z2 - z3 ) ) +
01063                           sqrt( ( x2 - x4 ) * ( x2 - x4 ) + ( y2 - y4 ) * ( y2 - y4 ) + ( z2 - z4 ) * ( z2 - z4 ) ) +
01064                           sqrt( ( x3 - x4 ) * ( x3 - x4 ) + ( y3 - y4 ) * ( y3 - y4 ) + ( z3 - z4 ) * ( z3 - z4 ) ) ) /
01065                         6.;
01066 
01067         /* Use the length scale to get a better idea if the tet is flat or
01068            just really small. */
01069         if( fabs( dScale ) < EPSILON ) { return false; }
01070         else
01071         {
01072             dDet /= ( dScale * dScale * dScale );
01073         }
01074 
01075         if( dDet > dEps ) { valid = true; }
01076         else if( dDet < -dEps )
01077         {
01078             return false;
01079         }
01080         else
01081         {
01082             return false;
01083         }
01084     }  // end for i=1,numElements
01085 
01086     //  MSQ_DEBUG_ACTION(2,{fprintf(stdout,"Mesh Validity is: %d \n",valid);});
01087 
01088     //  FUNCTION_TIMER_END();
01089     return ( valid );
01090 }
01091 
01092 void NonSmoothDescent::get_active_directions( const std::vector< Vector3D >& pGradient, std::vector< Vector3D >& dir,
01093                                               MsqError& /*err*/ )
01094 {
01095     const size_t num_active = mActive.active_ind.size();
01096     dir.resize( num_active );
01097     for( size_t i = 0; i < num_active; i++ )
01098     {
01099         dir[i] = pGradient[mActive.active_ind[i]];
01100     }
01101 }
01102 
01103 bool NonSmoothDescent::check_vector_dots( const std::vector< Vector3D >& vec, const Vector3D& normal,
01104                                           MsqError& /*err*/ )
01105 {
01106     double test_dot = vec[0] % normal;
01107     unsigned ind    = 1;
01108     while( ( fabs( test_dot ) < EPSILON ) && ( ind < vec.size() ) )
01109     {
01110         test_dot = vec[ind] % normal;
01111         ind++;
01112     }
01113 
01114     for( unsigned i = ind; i < vec.size(); i++ )
01115     {
01116         double dot = vec[i] % normal;
01117         if( ( ( dot > 0 && test_dot < 0 ) || ( dot < 0 && test_dot > 0 ) ) && ( fabs( dot ) > EPSILON ) )
01118         {
01119             return true;
01120         }
01121     }
01122     return false;
01123 }
01124 
01125 bool NonSmoothDescent::convex_hull_test( const std::vector< Vector3D >& vec, MsqError& err )
01126 {
01127     //    int ierr;
01128     bool equil         = false;
01129     Direction dir_done = MSQ_ZDIR;
01130     Status status      = MSQ_CHECK_Z_COORD_DIRECTION;
01131     Vector3D pt1, pt2, pt3, normal;
01132 
01133     /* tries to determine equilibrium for the 3D case */
01134 
01135     if( vec.size() <= 2 ) status = MSQ_NO_EQUIL;
01136 
01137     while( ( status != MSQ_EQUIL ) && ( status != MSQ_NO_EQUIL ) && ( status != MSQ_HULL_TEST_ERROR ) )
01138     {
01139         if( status == MSQ_CHECK_Z_COORD_DIRECTION )
01140         {
01141             this->find_plane_points( MSQ_ZDIR, MSQ_YDIR, vec, pt1, pt2, pt3, status, err );
01142             dir_done = MSQ_ZDIR;
01143         }
01144         if( status == MSQ_CHECK_Y_COORD_DIRECTION )
01145         {
01146             this->find_plane_points( MSQ_YDIR, MSQ_XDIR, vec, pt1, pt2, pt3, status, err );
01147             dir_done = MSQ_YDIR;
01148         }
01149         if( status == MSQ_CHECK_X_COORD_DIRECTION )
01150         {
01151             this->find_plane_points( MSQ_XDIR, MSQ_ZDIR, vec, pt1, pt2, pt3, status, err );
01152             dir_done = MSQ_XDIR;
01153         }
01154         if( status == MSQ_TWO_PT_PLANE ) { pt3 = Vector3D( 0, 0, 0 ); }
01155         if( ( status == MSQ_TWO_PT_PLANE ) || ( status == MSQ_THREE_PT_PLANE ) )
01156         {
01157             this->find_plane_normal( pt1, pt2, pt3, normal, err );
01158             equil = this->check_vector_dots( vec, normal, err );
01159             if( equil )
01160             {
01161                 switch( dir_done )
01162                 {
01163                     case MSQ_ZDIR:
01164                         equil  = false;
01165                         status = MSQ_CHECK_Y_COORD_DIRECTION;
01166                         break;
01167                     case MSQ_YDIR:
01168                         equil  = false;
01169                         status = MSQ_CHECK_X_COORD_DIRECTION;
01170                         break;
01171                     case MSQ_XDIR:
01172                         equil  = true;
01173                         status = MSQ_EQUIL;
01174                 }
01175             }
01176             else if( equil == 0 )
01177             {
01178                 status = MSQ_NO_EQUIL;
01179             }
01180             else
01181             {
01182                 MSQ_SETERR( err )( "equil flag not set to 0 or 1", MsqError::INVALID_STATE );
01183             }
01184         }
01185     }
01186     switch( status )
01187     {
01188         case MSQ_NO_EQUIL:
01189             MSQ_PRINT( 3 )( "Not an equilibrium point\n" );
01190             equil = false;
01191             break;
01192         case MSQ_EQUIL:
01193             MSQ_PRINT( 3 )( "An equilibrium point\n" );
01194             equil = true;
01195             break;
01196         default:
01197             MSQ_PRINT( 3 )( "Failed to determine equil or not; status = %d\n", status );
01198     }
01199     //    FUNCTION_TIMER_END();
01200     return ( equil );
01201 }
01202 
01203 void NonSmoothDescent::form_grammian( const std::vector< Vector3D >& vec, MsqError& /*err*/ )
01204 {
01205     /* form the grammian with the dot products of the gradients */
01206     const size_t num_active = mActive.active_ind.size();
01207     mG.resize( num_active );
01208     for( size_t i = 0; i < num_active; i++ )
01209         for( size_t j = i; j < num_active; j++ )
01210             mG( i, j ) = vec[i] % vec[j];
01211 }
01212 
01213 bool NonSmoothDescent::check_equilibrium( OptStatus& status, MsqError& err )
01214 {
01215     std::vector< Vector3D > dir;
01216 
01217     // TODO - this subroutine is no longer clear to me... is it still
01218     // appropriate for quads and hexes?  I think it might be in 2D, but
01219     // 3D is less clear.  Is there a more general algorithm to use?
01220     // ask Todd/check in numerical optimization
01221 
01222     bool equil              = false;
01223     const size_t num_active = mActive.active_ind.size();
01224     ;
01225 
01226     if( num_active == qmHandles.size() )
01227     {
01228         equil  = true;
01229         status = MSQ_EQUILIBRIUM;
01230         MSQ_PRINT( 3 )( "All the function values are in the active set\n" );
01231     }
01232 
01233     /* set up the active mGradient directions */
01234     this->get_active_directions( mGradient, dir, err );
01235     MSQ_ERRZERO( err );
01236 
01237     /* normalize the active directions */
01238     for( size_t j = 0; j < num_active; j++ )
01239         dir[j] /= dir[j].length();
01240 
01241     equil = this->convex_hull_test( dir, err );
01242     if( equil ) status = MSQ_EQUILIBRIUM;
01243 
01244     return equil;
01245 }
01246 
01247 static double condition3x3( const double A[9] )
01248 {
01249     //   int ierr;
01250     double a11, a12, a13;
01251     double a21, a22, a23;
01252     double a31, a32, a33;
01253     //   double s1, s2, s4, s3, t0;
01254     double s1, s2, s3;
01255     double denom;
01256     //   double one = 1.0;
01257     double temp;
01258     bool zero_denom = true;
01259 
01260     a11 = A[0];
01261     a12 = A[1];
01262     a13 = A[2];
01263     a21 = A[3];
01264     a22 = A[4];
01265     a23 = A[5];
01266     a31 = A[6];
01267     a32 = A[7];
01268     a33 = A[8];
01269 
01270     denom = -a11 * a22 * a33 + a11 * a23 * a32 + a21 * a12 * a33 - a21 * a13 * a32 - a31 * a12 * a23 + a31 * a13 * a22;
01271 
01272     if( ( fabs( a11 ) > EPSILON ) && ( fabs( denom / a11 ) > EPSILON ) ) { zero_denom = false; }
01273     if( ( fabs( a22 ) > EPSILON ) && ( fabs( denom / a22 ) > EPSILON ) ) { zero_denom = false; }
01274     if( ( fabs( a33 ) > EPSILON ) && ( fabs( denom / a33 ) > EPSILON ) ) { zero_denom = false; }
01275 
01276     if( zero_denom ) { return HUGE_VAL; }
01277     else
01278     {
01279         s1 = sqrt( a11 * a11 + a12 * a12 + a13 * a13 + a21 * a21 + a22 * a22 + a23 * a23 + a31 * a31 + a32 * a32 +
01280                    a33 * a33 );
01281 
01282         temp = ( -a22 * a33 + a23 * a32 ) / denom;
01283         s3   = temp * temp;
01284         temp = ( a12 * a33 - a13 * a32 ) / denom;
01285         s3 += temp * temp;
01286         temp = ( a12 * a23 - a13 * a22 ) / denom;
01287         s3 += temp * temp;
01288         temp = ( a21 * a33 - a23 * a31 ) / denom;
01289         s3 += temp * temp;
01290         temp = ( a11 * a33 - a13 * a31 ) / denom;
01291         s3 += temp * temp;
01292         temp = ( a11 * a23 - a13 * a21 ) / denom;
01293         s3 += temp * temp;
01294         temp = ( a21 * a32 - a22 * a31 ) / denom;
01295         s3 += temp * temp;
01296         temp = ( -a11 * a32 + a12 * a31 ) / denom;
01297         s3 += temp * temp;
01298         temp = ( -a11 * a22 + a12 * a21 ) / denom;
01299         s3 += temp * temp;
01300 
01301         s2 = sqrt( s3 );
01302         return s1 * s2;
01303     }
01304 }
01305 
01306 double NonSmoothDescent::SymmetricMatrix::condition3x3() const
01307 {
01308     double values[9] = { operator()( 0, 0 ), operator()( 0, 1 ), operator()( 0, 2 ),
01309                          operator()( 1, 0 ), operator()( 1, 1 ), operator()( 1, 2 ),
01310                          operator()( 2, 0 ), operator()( 2, 1 ), operator()( 2, 2 ) };
01311     return MBMesquite::condition3x3( values );
01312 }
01313 
01314 void NonSmoothDescent::singular_test( int n, const Matrix3D& A, bool& singular, MsqError& err )
01315 {
01316     //    int test;
01317     //    double determinant;
01318     double cond;
01319 
01320     if( ( n > 3 ) || ( n < 1 ) )
01321     {
01322         MSQ_SETERR( err )( "Singular test works only for n=1 to n=3", MsqError::INVALID_ARG );
01323         return;
01324     }
01325 
01326     singular = true;
01327     switch( n )
01328     {
01329         case 1:
01330             if( A[0][0] > 0 ) singular = false;
01331             break;
01332         case 2:
01333             if( fabs( A[0][0] * A[1][1] - A[0][1] * A[1][0] ) > EPSILON ) singular = false;
01334             break;
01335         case 3:
01336             /* calculate the condition number */
01337             cond = condition3x3( A[0] );
01338             if( cond < 1E14 ) singular = false;
01339             break;
01340     }
01341 }
01342 
01343 void NonSmoothDescent::form_PD_grammian( MsqError& err )
01344 {
01345     int i, j, k;
01346     int g_ind_1;
01347     bool singular = false;
01348 
01349     const int num_active = mActive.active_ind.size();
01350 
01351     /* this assumes the grammian has been formed */
01352     for( i = 0; i < num_active; i++ )
01353     {
01354         for( j = i; j < num_active; j++ )
01355         {
01356             if( mG( i, j ) == -1 )
01357             {
01358                 MSQ_SETERR( err )( "Grammian not computed properly", MsqError::INVALID_STATE );
01359                 return;
01360             }
01361         }
01362     }
01363 
01364     /* use the first gradient in the active set */
01365     g_ind_1    = 0;
01366     mPDG[0][0] = mG( 0, 0 );
01367     pdgInd[0]  = mActive.active_ind[0];
01368 
01369     /* test the rest and add them as appropriate */
01370     k = 1;
01371     i = 1;
01372     while( ( k < 3 ) && ( i < num_active ) )
01373     {
01374         mPDG[0][k] = mPDG[k][0] = mG( 0, i );
01375         mPDG[k][k]              = mG( i, i );
01376         if( k == 2 )
01377         { /* add the dot product of g1 and g2 */
01378             mPDG[1][k] = mPDG[k][1] = mG( g_ind_1, i );
01379         }
01380         this->singular_test( k + 1, mPDG, singular, err );
01381         if( !singular )
01382         {
01383             pdgInd[k] = mActive.active_ind[i];
01384             if( k == 1 ) g_ind_1 = i;
01385             k++;
01386         }
01387         i++;
01388     }
01389 }
01390 
01391 void NonSmoothDescent::search_edges_faces( const Vector3D* dir, Vector3D& mSearch, MsqError& /*err*/ )
01392 {
01393     bool viable;
01394     double a, b, denom;
01395     Vector3D g_bar;
01396     Vector3D temp_search( 0, 0, 0 ); /* initialize the search direction to 0,0 */
01397     double projection, min_projection;
01398 
01399     const size_t num_active = mActive.active_ind.size();
01400     ;
01401 
01402     /* Check for viable faces */
01403     min_projection = HUGE_VAL;
01404     for( size_t i = 0; i < num_active; i++ )
01405     {
01406         /* FACE I */
01407         viable = true;
01408 
01409         /* test the viability */
01410         for( size_t j = 0; j < num_active; j++ )
01411         { /* lagrange multipliers>0 */
01412             if( mG( j, i ) < 0 ) viable = false;
01413         }
01414 
01415         /* find the minimum of viable directions */
01416         if( ( viable ) && ( mG( i, i ) < min_projection ) )
01417         {
01418             min_projection = mG( i, i );
01419             temp_search    = dir[i];
01420             mSteepest      = mActive.active_ind[i];
01421         }
01422 
01423         /* INTERSECTION IJ */
01424         for( size_t j = i + 1; j < num_active; j++ )
01425         {
01426             viable = true;
01427 
01428             /* find the coefficients of the intersection
01429                and test the viability */
01430             denom = 2 * mG( i, j ) - mG( i, i ) - mG( j, j );
01431             a = b = 0;
01432             if( fabs( denom ) > EPSILON )
01433             {
01434                 b = ( mG( i, j ) - mG( i, i ) ) / denom;
01435                 a = 1 - b;
01436                 if( ( b < 0 ) || ( b > 1 ) ) /* 0 < b < 1 */
01437                     viable = false;
01438                 for( size_t k = 0; k < num_active; k++ )
01439                 { /* lagrange multipliers>0 */
01440                     if( ( a * mG( k, i ) + b * mG( k, j ) ) <= 0 ) viable = false;
01441                 }
01442             }
01443             else
01444             {
01445                 viable = false; /* Linearly dependent */
01446             }
01447 
01448             /* find the minimum of viable directions */
01449             if( viable )
01450             {
01451                 g_bar      = a * dir[i] + b * dir[j];
01452                 projection = g_bar % g_bar;
01453                 if( projection < min_projection )
01454                 {
01455                     min_projection = projection;
01456                     temp_search    = g_bar;
01457                     mSteepest      = mActive.active_ind[i];
01458                 }
01459             }
01460         }
01461     }
01462     if( optStatus != MSQ_EQUILIBRIUM ) { mSearch = temp_search; }
01463 }
01464 
01465 bool NonSmoothDescent::solve2x2( double a11, double a12, double a21, double a22, double b1, double b2, double x[2],
01466                                  MsqError& /*err*/ )
01467 {
01468     double factor;
01469 
01470     /* if the system is not singular, solve it */
01471     if( fabs( a11 * a22 - a21 * a12 ) > EPSILON )
01472     {
01473         if( fabs( a11 ) > EPSILON )
01474         {
01475             factor = ( a21 / a11 );
01476             x[1]   = ( b2 - factor * b1 ) / ( a22 - factor * a12 );
01477             x[0]   = ( b1 - a12 * x[1] ) / a11;
01478         }
01479         else if( fabs( a21 ) > EPSILON )
01480         {
01481             factor = ( a11 / a21 );
01482             x[1]   = ( b1 - factor * b2 ) / ( a12 - factor * a22 );
01483             x[0]   = ( b2 - a22 * x[1] ) / a21;
01484         }
01485         return true;
01486     }
01487     else
01488     {
01489         return false;
01490     }
01491 }
01492 
01493 void NonSmoothDescent::form_reduced_matrix( SymmetricMatrix& P )
01494 {
01495     const size_t P_size = mActive.active_ind.size() - 1;
01496     P.resize( P_size );
01497     for( size_t i = 0; i < P_size; i++ )
01498     {
01499         P( i, i ) = mG( 0, 0 ) - 2 * mG( 0, i + 1 ) + mG( i + 1, i + 1 );
01500         for( size_t j = i + 1; j < P_size; j++ )
01501         {
01502             P( i, j ) = mG( 0, 0 ) - mG( 0, j + 1 ) - mG( i + 1, 0 ) + mG( i + 1, j + 1 );
01503         }
01504     }
01505 }
01506 
01507 void NonSmoothDescent::get_min_estimate( double* final_est, MsqError& /*err*/ )
01508 {
01509     double est_imp;
01510 
01511     *final_est = -HUGE_VAL;
01512     for( size_t i = 0; i < mActive.active_ind.size(); i++ )
01513     {
01514         est_imp = -mAlpha * mGS[mActive.active_ind[i]];
01515         if( est_imp > *final_est ) *final_est = est_imp;
01516     }
01517     if( *final_est == 0 )
01518     {
01519         *final_est = -HUGE_VAL;
01520         for( size_t i = 0; i < qmHandles.size(); i++ )
01521         {
01522             est_imp = -mAlpha * mGS[i];
01523             if( ( est_imp > *final_est ) && ( fabs( est_imp ) > EPSILON ) ) { *final_est = est_imp; }
01524         }
01525     }
01526 }
01527 
01528 void NonSmoothDescent::get_gradient_projections( const Vector3D& mSearch, MsqError& /*err*/ )
01529 {
01530     for( size_t i = 0; i < qmHandles.size(); i++ )
01531         mGS[i] = mGradient[i] % mSearch;
01532 
01533     MSQ_PRINT( 3 )( "steepest in get_gradient_projections %lu\n", (unsigned long)mSteepest );
01534 }
01535 
01536 void NonSmoothDescent::compute_alpha( MsqError& /*err*/ )
01537 {
01538     double steepest_function;
01539     double steepest_grad;
01540     double alpha_i;
01541     double min_positive_value = HUGE_VAL;
01542 
01543     //    FUNCTION_TIMER_START("Compute Alpha");
01544 
01545     MSQ_PRINT( 2 )( "In compute alpha\n" );
01546 
01547     mAlpha = HUGE_VAL;
01548 
01549     steepest_function = mFunction[mSteepest];
01550     steepest_grad     = mGS[mSteepest];
01551     for( size_t i = 0; i < qmHandles.size(); i++ )
01552     {
01553         /* if it's not active */
01554         if( i != mSteepest )
01555         {
01556             alpha_i = steepest_function - mFunction[i];
01557 
01558             if( fabs( mGS[mSteepest] - mGS[i] ) > 1E-13 )
01559             {
01560                 /* compute line intersection */
01561                 alpha_i = alpha_i / ( steepest_grad - mGS[i] );
01562             }
01563             else
01564             {
01565                 /* the lines don't intersect - it's not under consideration*/
01566                 alpha_i = 0;
01567             }
01568             if( ( alpha_i > minStepSize ) && ( fabs( alpha_i ) < fabs( mAlpha ) ) )
01569             {
01570                 mAlpha = fabs( alpha_i );
01571                 MSQ_PRINT( 3 )( "Setting alpha %lu %g\n", (unsigned long)i, alpha_i );
01572             }
01573             if( ( alpha_i > 0 ) && ( alpha_i < min_positive_value ) ) { min_positive_value = alpha_i; }
01574         }
01575     }
01576 
01577     if( ( mAlpha == HUGE_VAL ) && ( min_positive_value != HUGE_VAL ) ) { mAlpha = min_positive_value; }
01578 
01579     /* if it never gets set, set it to the default */
01580     if( mAlpha == HUGE_VAL )
01581     {
01582         mAlpha = maxAlpha;
01583         MSQ_PRINT( 3 )( "Setting alpha to the maximum step length\n" );
01584     }
01585 
01586     MSQ_PRINT( 3 )( "  The initial step size: %f\n", mAlpha );
01587 
01588     //    FUNCTION_TIMER_END();
01589 }
01590 
01591 void NonSmoothDescent::print_active_set( const ActiveSet& active_set, const std::vector< double >& func )
01592 {
01593     if( active_set.active_ind.empty() ) MSQ_DBGOUT( 3 ) << "No active values\n";
01594     /* print the active set */
01595     for( size_t i = 0; i < active_set.active_ind.size(); i++ )
01596     {
01597         MSQ_PRINT( 3 )
01598         ( "Active value %lu:   %f \n", (unsigned long)i + 1, func[active_set.active_ind[i]] );
01599     }
01600 }
01601 
01602 void NonSmoothDescent::init_opt( PatchData& pd, MsqError& err )
01603 {
01604     qmHandles.clear();
01605     currentQM->get_evaluations( pd, qmHandles, true, err );MSQ_ERRRTN( err );
01606 
01607     MSQ_PRINT( 2 )( "\nInitializing Optimization \n" );
01608 
01609     /* for the purposes of initialization will be set to zero after */
01610     optStatus = MSQ_NO_STATUS;
01611     mSteepest = 0;
01612     mAlpha    = 0;
01613     maxAlpha  = 0;
01614 
01615     MSQ_PRINT( 3 )( "  Initialized Constants \n" );
01616     pdgInd[0] = pdgInd[1] = pdgInd[2] = -1;
01617     mPDG                              = Matrix3D( 0.0 );
01618 
01619     MSQ_PRINT( 3 )( "  Initialized search and PDG \n" );
01620     mFunction.clear();
01621     mFunction.resize( qmHandles.size(), 0.0 );
01622     testFunction.clear();
01623     testFunction.resize( qmHandles.size(), 0.0 );
01624     originalFunction.clear();
01625     originalFunction.resize( qmHandles.size(), 0.0 );
01626     mGradient.clear();
01627     mGradient.resize( qmHandles.size(), Vector3D( 0.0 ) );
01628     mGS.resize( qmHandles.size() );
01629 
01630     MSQ_PRINT( 3 )( "  Initialized function/gradient \n" );
01631     mG.resize( qmHandles.size() );
01632     mG.fill( -1 );
01633     MSQ_PRINT( 3 )( "  Initialized G\n" );
01634 
01635     prevActiveValues.clear();
01636     prevActiveValues.reserve( 32 );
01637     MSQ_PRINT( 3 )( "  Initialized prevActiveValues\n" );
01638 }
01639 
01640 void NonSmoothDescent::init_max_step_length( PatchData& pd, MsqError& err )
01641 {
01642     size_t i, j;
01643     double max_diff = 0;
01644     double diff     = 0;
01645 
01646     MSQ_PRINT( 2 )( "In init_max_step_length\n" );
01647 
01648     /* check that the input data is correct */
01649     if( pd.num_elements() == 0 )
01650     {
01651         MSQ_SETERR( err )( "Num incident vtx = 0\n", MsqError::INVALID_MESH );
01652         return;
01653     }
01654 
01655     /* find the maximum distance between two incident vertex locations */
01656     const MsqVertex* coords = pd.get_vertex_array( err );
01657     for( i = 0; i < pd.num_nodes() - 1; i++ )
01658     {
01659         for( j = i; j < pd.num_nodes(); j++ )
01660         {
01661             diff = ( coords[i] - coords[j] ).length_squared();
01662             if( max_diff < diff ) max_diff = diff;
01663         }
01664     }
01665     max_diff = sqrt( max_diff );
01666     if( max_diff == 0 )
01667     {
01668         MSQ_SETERR( err )
01669         ( "Maximum distance between incident vertices = 0\n", MsqError::INVALID_MESH );
01670         return;
01671     }
01672     maxAlpha = max_diff / 100;
01673 
01674     MSQ_PRINT( 3 )( "  Maximum step is %g\n", maxAlpha );
01675 }
01676 
01677 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines