MOAB: Mesh Oriented datABase
(version 5.3.1)
|
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