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 /*! 00028 \file ConjugateGradient.cpp 00029 \brief 00030 00031 The Conjugate Gradient class is a concrete vertex mover 00032 which performs conjugate gradient minimizaiton. 00033 00034 \author Michael Brewer 00035 \date 2002-06-19 00036 */ 00037 00038 #include "ConjugateGradient.hpp" 00039 #include <cmath> 00040 #include "MsqDebug.hpp" 00041 #include "MsqTimer.hpp" 00042 //#include "MsqFreeVertexIndexIterator.hpp" 00043 00044 namespace MBMesquite 00045 { 00046 00047 extern int get_parallel_rank(); 00048 extern int get_parallel_size(); 00049 00050 std::string ConjugateGradient::get_name() const 00051 { 00052 return "ConjugateGradient"; 00053 } 00054 00055 PatchSet* ConjugateGradient::get_patch_set() 00056 { 00057 return PatchSetUser::get_patch_set(); 00058 } 00059 00060 ConjugateGradient::ConjugateGradient( ObjectiveFunction* objective ) 00061 : VertexMover( objective ), PatchSetUser( true ), pMemento( NULL ), conjGradDebug( 0 ) 00062 { 00063 } 00064 00065 ConjugateGradient::ConjugateGradient( ObjectiveFunction* objective, MsqError& err ) 00066 : VertexMover( objective ), PatchSetUser( true ), pMemento( NULL ), conjGradDebug( 0 ) 00067 { 00068 // Michael:: default to global? 00069 set_debugging_level( 0 ); 00070 // set the default inner termination criterion 00071 TerminationCriterion* default_crit = get_inner_termination_criterion(); 00072 if( default_crit == NULL ) 00073 { 00074 MSQ_SETERR( err ) 00075 ( "QualityImprover did not create a default inner " 00076 "termination criterion.", 00077 MsqError::INVALID_STATE ); 00078 return; 00079 } 00080 else 00081 { 00082 default_crit->add_iteration_limit( 5 );MSQ_ERRRTN( err ); 00083 } 00084 } 00085 00086 ConjugateGradient::~ConjugateGradient() 00087 { 00088 // Checks that cleanup() has been called. 00089 assert( pMemento == NULL ); 00090 } 00091 00092 void ConjugateGradient::initialize( PatchData& pd, MsqError& err ) 00093 { 00094 if( get_parallel_size() ) 00095 { 00096 MSQ_DBGOUT( 2 ) << "\nP[" << get_parallel_rank() << "] " 00097 << "o Performing Conjugate Gradient optimization.\n"; 00098 } 00099 else 00100 { 00101 MSQ_DBGOUT( 2 ) << "\no Performing Conjugate Gradient optimization.\n"; 00102 } 00103 pMemento = pd.create_vertices_memento( err ); 00104 } 00105 00106 void ConjugateGradient::initialize_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {} 00107 00108 /*!Performs Conjugate gradient minimization on the PatchData, pd.*/ 00109 void ConjugateGradient::optimize_vertex_positions( PatchData& pd, MsqError& err ) 00110 { 00111 // pd.reorder(); 00112 00113 MSQ_FUNCTION_TIMER( "ConjugateGradient::optimize_vertex_positions" ); 00114 00115 Timer c_timer; 00116 size_t num_vert = pd.num_free_vertices(); 00117 if( num_vert < 1 ) 00118 { 00119 MSQ_DBGOUT( 1 ) << "\nEmpty free vertex list in ConjugateGradient\n"; 00120 return; 00121 } 00122 /* 00123 //zero out arrays 00124 int zero_loop=0; 00125 while(zero_loop<arraySize){ 00126 fGrad[zero_loop].set(0,0,0); 00127 pGrad[zero_loop].set(0,0,0); 00128 fNewGrad[zero_loop].set(0,0,0); 00129 ++zero_loop; 00130 } 00131 */ 00132 00133 // get OF evaluator 00134 OFEvaluator& objFunc = get_objective_function_evaluator(); 00135 00136 size_t ind; 00137 // Michael cull list: possibly set soft_fixed flags here 00138 00139 // MsqFreeVertexIndexIterator free_iter(pd, err); MSQ_ERRRTN(err); 00140 00141 double f = 0; 00142 // Michael, this isn't equivalent to CUBIT because we only want to check 00143 // the objective function value of the 'bad' elements 00144 // if invalid initial patch set an error. 00145 bool temp_bool = objFunc.update( pd, f, fGrad, err ); 00146 assert( fGrad.size() == num_vert ); 00147 if( MSQ_CHKERR( err ) ) return; 00148 if( !temp_bool ) 00149 { 00150 MSQ_SETERR( err ) 00151 ( "Conjugate Gradient not able to get valid gradient " 00152 "and function values on intial patch.", 00153 MsqError::INVALID_MESH ); 00154 return; 00155 } 00156 double grad_norm = MSQ_MAX_CAP; 00157 00158 if( conjGradDebug > 0 ) 00159 { 00160 MSQ_PRINT( 2 )( "\nCG's DEGUB LEVEL = %i \n", conjGradDebug ); 00161 grad_norm = Linf( arrptr( fGrad ), fGrad.size() ); 00162 MSQ_PRINT( 2 )( "\nCG's FIRST VALUE = %f,grad_norm = %f", f, grad_norm ); 00163 MSQ_PRINT( 2 )( "\n TIME %f", c_timer.since_birth() ); 00164 grad_norm = MSQ_MAX_CAP; 00165 } 00166 00167 // Initializing pGrad (search direction). 00168 pGrad.resize( fGrad.size() ); 00169 for( ind = 0; ind < num_vert; ++ind ) 00170 pGrad[ind] = ( -fGrad[ind] ); 00171 00172 int j = 0; // total nb of step size changes ... not used much 00173 int i = 0; // iteration counter 00174 unsigned m = 0; // 00175 double alp = MSQ_MAX_CAP; // alp: scale factor of search direction 00176 // we know inner_criterion is false because it was checked in 00177 // loop_over_mesh before being sent here. 00178 TerminationCriterion* term_crit = get_inner_termination_criterion(); 00179 00180 // while ((i<maxIteration && alp>stepBound && grad_norm>normGradientBound) 00181 // && !inner_criterion){ 00182 while( !term_crit->terminate() ) 00183 { 00184 ++i; 00185 // std::cout<<"\Michael delete i = "<<i; 00186 int k = 0; 00187 alp = get_step( pd, f, k, err ); 00188 j += k; 00189 if( conjGradDebug > 2 ) { MSQ_PRINT( 2 )( "\n Alp initial, alp = %20.18f", alp ); } 00190 00191 // if alp == 0, revert to steepest descent search direction 00192 if( alp == 0 ) 00193 { 00194 for( m = 0; m < num_vert; ++m ) 00195 { 00196 pGrad[m] = ( -fGrad[m] ); 00197 } 00198 alp = get_step( pd, f, k, err ); 00199 j += k; 00200 if( conjGradDebug > 1 ) 00201 { 00202 MSQ_PRINT( 2 )( "\n CG's search direction reset." ); 00203 if( conjGradDebug > 2 ) MSQ_PRINT( 2 )( "\n Alp was zero, alp = %20.18f", alp ); 00204 } 00205 } 00206 if( alp != 0 ) 00207 { 00208 pd.move_free_vertices_constrained( arrptr( pGrad ), num_vert, alp, err );MSQ_ERRRTN( err ); 00209 00210 if( !objFunc.update( pd, f, fNewGrad, err ) ) 00211 { 00212 MSQ_SETERR( err ) 00213 ( "Error inside Conjugate Gradient, vertices moved " 00214 "making function value invalid.", 00215 MsqError::INVALID_MESH ); 00216 return; 00217 } 00218 assert( fNewGrad.size() == (unsigned)num_vert ); 00219 00220 if( conjGradDebug > 0 ) 00221 { 00222 grad_norm = Linf( arrptr( fNewGrad ), num_vert ); 00223 MSQ_PRINT( 2 ) 00224 ( "\nCG's VALUE = %f, iter. = %i, grad_norm = %f, alp = %f", f, i, grad_norm, alp ); 00225 MSQ_PRINT( 2 )( "\n TIME %f", c_timer.since_birth() ); 00226 } 00227 double s11 = 0; 00228 double s12 = 0; 00229 double s22 = 0; 00230 // free_iter.reset(); 00231 // while (free_iter.next()) { 00232 // m=free_iter.value(); 00233 for( m = 0; m < num_vert; ++m ) 00234 { 00235 s11 += fGrad[m] % fGrad[m]; 00236 s12 += fGrad[m] % fNewGrad[m]; 00237 s22 += fNewGrad[m] % fNewGrad[m]; 00238 } 00239 00240 // Steepest Descent (takes 2-3 times as long as P-R) 00241 // double bet=0; 00242 00243 // Fletcher-Reeves (takes twice as long as P-R) 00244 // double bet = s22/s11; 00245 00246 // Polack-Ribiere 00247 double bet; 00248 if( !divide( s22 - s12, s11, bet ) ) return; // gradient is zero 00249 // free_iter.reset(); 00250 // while (free_iter.next()) { 00251 // m=free_iter.value(); 00252 for( m = 0; m < num_vert; ++m ) 00253 { 00254 pGrad[m] = ( -fNewGrad[m] + ( bet * pGrad[m] ) ); 00255 fGrad[m] = fNewGrad[m]; 00256 } 00257 if( conjGradDebug > 2 ) 00258 { 00259 MSQ_PRINT( 2 ) 00260 ( " \nSEARCH DIRECTION INFINITY NORM = %e", Linf( arrptr( fNewGrad ), num_vert ) ); 00261 } 00262 00263 } // end if on alp == 0 00264 00265 term_crit->accumulate_patch( pd, err );MSQ_ERRRTN( err ); 00266 term_crit->accumulate_inner( pd, f, arrptr( fGrad ), err );MSQ_ERRRTN( err ); 00267 } // end while 00268 if( conjGradDebug > 0 ) 00269 { 00270 MSQ_PRINT( 2 )( "\nConjugate Gradient complete i=%i ", i ); 00271 MSQ_PRINT( 2 )( "\n- FINAL value = %f, alp=%4.2e grad_norm=%4.2e", f, alp, grad_norm ); 00272 MSQ_PRINT( 2 )( "\n FINAL TIME %f", c_timer.since_birth() ); 00273 } 00274 } 00275 00276 void ConjugateGradient::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) 00277 { 00278 // cout << "- Executing ConjugateGradient::iteration_complete()\n"; 00279 } 00280 00281 void ConjugateGradient::cleanup() 00282 { 00283 // cout << "- Executing ConjugateGradient::iteration_end()\n"; 00284 fGrad.clear(); 00285 pGrad.clear(); 00286 fNewGrad.clear(); 00287 // pMemento->~PatchDataVerticesMemento(); 00288 delete pMemento; 00289 pMemento = NULL; 00290 } 00291 00292 //! Computes a distance to move vertices given an initial position and search direction (stored in 00293 //! data member pGrad). 00294 /*!Returns alp, the double which scales the search direction vector 00295 which when added to the old nodal positions yields the new nodal 00296 positions.*/ 00297 /*!\todo Michael NOTE: ConjugateGradient::get_step's int &j is only 00298 to remain consisitent with CUBIT for an initial test. It can be 00299 removed.*/ 00300 00301 double ConjugateGradient::get_step( PatchData& pd, double f0, int& j, MsqError& err ) 00302 { 00303 // get OF evaluator 00304 OFEvaluator& objFunc = get_objective_function_evaluator(); 00305 00306 size_t num_vertices = pd.num_free_vertices(); 00307 // initial guess for alp 00308 double alp = 1.0; 00309 int jmax = 100; 00310 double rho = 0.5; 00311 // feasible=false implies the mesh is not in the feasible region 00312 bool feasible = false; 00313 int found = 0; 00314 // f and fnew hold the objective function value 00315 double f = 0; 00316 double fnew = 0; 00317 // Counter to avoid infinitly scaling alp 00318 j = 0; 00319 // save memento 00320 pd.recreate_vertices_memento( pMemento, err ); 00321 // if we must check feasiblility 00322 // while step takes mesh into infeasible region and ... 00323 while( j < jmax && !feasible && alp > MSQ_MIN ) 00324 { 00325 ++j; 00326 pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err ); 00327 feasible = objFunc.evaluate( pd, f, err ); 00328 if( err.error_code() == err.BARRIER_VIOLATED ) 00329 err.clear(); // barrier violation does not represent an actual error here 00330 MSQ_ERRZERO( err ); 00331 // if not feasible, try a smaller alp (take smaller step) 00332 if( !feasible ) { alp *= rho; } 00333 } // end while ... 00334 00335 // if above while ended due to j>=jmax, no valid step was found. 00336 if( j >= jmax ) 00337 { 00338 MSQ_PRINT( 2 )( "\nFeasible Point Not Found" ); 00339 return 0.0; 00340 } 00341 // Message::print_info("\nOriginal f %f, first new f = %f, alp = %f",f0,f,alp); 00342 // if new f is larger than original, our step was too large 00343 if( f >= f0 ) 00344 { 00345 j = 0; 00346 while( j < jmax && found == 0 ) 00347 { 00348 ++j; 00349 alp *= rho; 00350 pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err ); 00351 // Get new obj value 00352 // if patch is now invalid, then the feasible region is convex or 00353 // we have an error. For now, we assume an error. 00354 if( !objFunc.evaluate( pd, f, err ) ) 00355 { 00356 MSQ_SETERR( err ) 00357 ( "Non-convex feasiblility region found.", MsqError::INVALID_MESH ); 00358 } 00359 pd.set_to_vertices_memento( pMemento, err ); 00360 MSQ_ERRZERO( err ); 00361 // if our step has now improved the objective function value 00362 if( f < f0 ) { found = 1; } 00363 } // end while j less than jmax 00364 // Message::print_info("\nj = %d found = %d f = %20.18f f0 = %20.18f\n",j,found,f,f0); 00365 // if above ended because of j>=jmax, take no step 00366 if( found == 0 ) 00367 { 00368 // Message::print_info("alp = %10.8f, but returning zero\n",alp); 00369 alp = 0.0; 00370 return alp; 00371 } 00372 00373 j = 0; 00374 // while shrinking the step improves the objFunc value further, 00375 // scale alp down. Return alp, when scaling once more would 00376 // no longer improve the objFunc value. 00377 while( j < jmax ) 00378 { 00379 ++j; 00380 alp *= rho; 00381 // step alp in search direction from original positions 00382 pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err ); 00383 MSQ_ERRZERO( err ); 00384 00385 // get new objective function value 00386 if( !objFunc.evaluate( pd, fnew, err ) ) MSQ_SETERR( err ) 00387 ( "Non-convex feasiblility region found while " 00388 "computing new f.", 00389 MsqError::INVALID_MESH ); 00390 if( fnew < f ) { f = fnew; } 00391 else 00392 { 00393 // Reset the vertices to original position 00394 pd.set_to_vertices_memento( pMemento, err ); 00395 MSQ_ERRZERO( err ); 00396 alp /= rho; 00397 return alp; 00398 } 00399 } 00400 // Reset the vertices to original position and return alp 00401 pd.set_to_vertices_memento( pMemento, err ); 00402 MSQ_ERRZERO( err ); 00403 return alp; 00404 } 00405 // else our new f was already smaller than our original 00406 else 00407 { 00408 j = 0; 00409 // check to see how large of step we can take 00410 while( j < jmax && found == 0 ) 00411 { 00412 ++j; 00413 // scale alp up (rho must be less than 1) 00414 alp /= rho; 00415 // step alp in search direction from original positions 00416 pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err ); 00417 MSQ_ERRZERO( err ); 00418 00419 feasible = objFunc.evaluate( pd, fnew, err ); 00420 if( err.error_code() == err.BARRIER_VIOLATED ) 00421 err.clear(); // evaluate() error does not represent an actual problem here 00422 MSQ_ERRZERO( err ); 00423 if( !feasible ) 00424 { 00425 alp *= rho; 00426 // Reset the vertices to original position and return alp 00427 pd.set_to_vertices_memento( pMemento, err ); 00428 MSQ_ERRZERO( err ); 00429 return alp; 00430 } 00431 if( fnew < f ) { f = fnew; } 00432 else 00433 { 00434 found = 1; 00435 alp *= rho; 00436 } 00437 } 00438 00439 // Reset the vertices to original position and return alp 00440 pd.set_to_vertices_memento( pMemento, err ); 00441 MSQ_ERRZERO( err ); 00442 return alp; 00443 } 00444 } 00445 00446 /*!Quadratic one-dimensional line search.*/ 00447 /* 00448 double ConjugateGradient::get_step(PatchData &pd,double f0,int &j, 00449 MsqError &err){ 00450 const double CGOLD = 0.3819660; 00451 const double ZEPS = 1.0e-10; 00452 int n=pd.num_free_vertices(); 00453 MsqVertex* vertices=pd.get_vertex_array(err); 00454 double a,b,d,etemp,fb,fu,fv,fw,fx,p,q,r,tol,tol1,tol2,u,v,w,x,xm; 00455 double e=0.0; 00456 d=0.0; 00457 tol=.001; 00458 int iter, maxiter; 00459 maxiter=100; 00460 a=0; 00461 b=.125; 00462 int m=0; 00463 fb=f0-1.0; 00464 iter=0; 00465 //find b such that a b 'should' bracket the min 00466 while (fb<=f0 && iter<maxiter){ 00467 ++iter; 00468 b*=2.0; 00469 for(m=0;m<n;++m){ 00470 mCoord[m]=mCoord[m] + (b*pGrad[m]); 00471 vertices[m]=(mCoord[m]); 00472 } 00473 fb=objFunc->evaluate(pd,err); 00474 } 00475 iter=0; 00476 x=w=v=(b/2.0); 00477 for(m=0;m<n;++m){ 00478 mCoord[m]=mCoord[m] + (x*pGrad[m]); 00479 vertices[m]=(mCoord[m]); 00480 } 00481 fw=fv=fx=objFunc->evaluate(pd,err); 00482 for(iter=0;iter<maxiter;++iter){ 00483 //Message::print_info("a=%f,b=%f,x=%f,iter=%i\n",a,b,x,iter); 00484 xm=(a+b)*.5; 00485 tol2=2.0*(tol1=tol*fabs(x)+ZEPS); 00486 if(fabs(x-xm)<= (tol2-0.5*(b-a))){ 00487 return x; 00488 } 00489 if(fabs(e)>tol1){ 00490 r=(x-w)*(fx-fv); 00491 q=(x-v)*(fx-fw); 00492 p=(x-v)*q-(x-w)*r; 00493 q=2.0*(q-r); 00494 if(q>0.0) 00495 p=-p; 00496 q=fabs(q); 00497 etemp=e; 00498 e=d; 00499 if(fabs(p)>=fabs(0.5*q*etemp)||(p<=q*(a-x))||(p>=q*(b-x))){ 00500 d=CGOLD*(e=(x>=xm?a-x:b-x)); 00501 } 00502 else{ 00503 d=p/q; 00504 u=x+d; 00505 if(u-a<tol2||b-u<tol2) 00506 { 00507 if(tol1<0.0) 00508 d=x-xm; 00509 else 00510 d=xm-x; 00511 } 00512 } 00513 } 00514 00515 else{ 00516 d=CGOLD*(e=(x>=xm?a-x:b-x)); 00517 } 00518 if(tol<0.0) 00519 u=(fabs(d)>=tol1?x+d:x-d); 00520 else 00521 u=(fabs(d)>=tol1?x+d:x+d); 00522 for(m=0;m<n;++m){ 00523 mCoord[m]=mCoord[m] + (u*pGrad[m]); 00524 vertices[m]=(mCoord[m]); 00525 } 00526 fu=objFunc->evaluate(pd,err); 00527 if(fu<fx){ 00528 if(u>=x) 00529 a=x; 00530 else 00531 b=x; 00532 v=w; 00533 w=x; 00534 x=u; 00535 fv=fw; 00536 fw=fx; 00537 fx=fu; 00538 } 00539 else{ 00540 if(u<x) 00541 a=u; 00542 else 00543 b=u; 00544 if(fu<=fw||w==x){ 00545 v=w; 00546 w=u; 00547 fv=fw; 00548 fw=fu; 00549 } 00550 else if (fu<=fv||v==x||v==w){ 00551 v=u; 00552 fv=fu; 00553 } 00554 } 00555 } 00556 for(m=0;m<n;++m){ 00557 vertices[m]=(mCoord[m]); 00558 } 00559 //PRINT_WARNING("TOO MANY ITERATIONS IN QUADRATIC LINE SEARCH"); 00560 return x; 00561 } 00562 */ 00563 00564 } // namespace MBMesquite