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