Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2004 Sandia Corporation and Argonne National
5 : : Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 : : with Sandia Corporation, the U.S. Government retains certain
7 : : rights in this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : [email protected], [email protected], [email protected],
24 : : [email protected], [email protected], [email protected]
25 : :
26 : : ***************************************************************** */
27 : : //
28 : : // AUTHOR: Thomas Leurent <[email protected]>
29 : : // ORG: Argonne National Laboratory
30 : : // E-MAIL: [email protected]
31 : : //
32 : : // ORIG-DATE: 15-Jan-03 at 08:05:56
33 : : // LAST-MOD: 15-Jun-04 at 15:45:00 by Thomas Leurent
34 : : //
35 : : // DESCRIPTION:
36 : : // ============
37 : : /*!
38 : : \file FeasibleNewton.cpp
39 : : \brief
40 : :
41 : : Implements the FeasibleNewton class member functions.
42 : :
43 : : \author Thomas Leurent
44 : : \author Todd Munson
45 : : \date 2003-01-15
46 : : */
47 : : // DESCRIP-END.
48 : : //
49 : :
50 : : #include "FeasibleNewton.hpp"
51 : : #include "MsqFreeVertexIndexIterator.hpp"
52 : : #include "MsqDebug.hpp"
53 : : #include "XYPlanarDomain.hpp"
54 : :
55 : : using namespace MBMesquite;
56 : :
57 : 0 : std::string FeasibleNewton::get_name() const
58 : : {
59 [ # # ]: 0 : return "FeasibleNewton";
60 : : }
61 : :
62 : 56 : PatchSet* FeasibleNewton::get_patch_set()
63 : : {
64 : 56 : return PatchSetUser::get_patch_set();
65 : : }
66 : :
67 : 57 : FeasibleNewton::FeasibleNewton( ObjectiveFunction* of )
68 [ + - ][ + - ]: 57 : : VertexMover( of ), PatchSetUser( true ), convTol( 1e-6 ), coordsMem( 0 ), havePrintedDirectionMessage( false )
69 : : {
70 [ + - ]: 57 : TerminationCriterion* default_crit = get_inner_termination_criterion();
71 [ + - ]: 57 : default_crit->add_absolute_gradient_L2_norm( 5e-5 );
72 : 57 : }
73 : :
74 : 69 : void FeasibleNewton::initialize( PatchData& pd, MsqError& err )
75 : : {
76 : : // Cannot do anything. Variable sizes with maximum size dependent
77 : : // upon the entire MeshSet.
78 [ - + ][ # # ]: 69 : coordsMem = pd.create_vertices_memento( err );MSQ_CHKERR( err );
79 : 69 : havePrintedDirectionMessage = false;
80 : 69 : }
81 : :
82 : 56 : void FeasibleNewton::initialize_mesh_iteration( PatchData& pd, MsqError& /*err*/ )
83 : : {
84 : 56 : pd.reorder();
85 : 56 : }
86 : :
87 : 112 : void FeasibleNewton::optimize_vertex_positions( PatchData& pd, MsqError& err )
88 : : {
89 : : MSQ_FUNCTION_TIMER( "FeasibleNewton::optimize_vertex_positions" );
90 : : MSQ_DBGOUT( 2 ) << "\no Performing Feasible Newton optimization.\n";
91 : :
92 : : //
93 : : // the only valid 2D meshes that FeasibleNewton works for are truly planar which
94 : : // lie in the X-Y coordinate plane.
95 : : //
96 : :
97 [ + + ]: 56 : XYPlanarDomain* xyPlanarDomainPtr = dynamic_cast< XYPlanarDomain* >( pd.get_domain() );
98 : : // only optimize if input mesh is a volume or an XYPlanarDomain
99 [ + + ][ + + ]: 56 : if( !pd.domain_set() || xyPlanarDomainPtr != NULL )
[ + + ]
100 : : {
101 : 55 : const double sigma = 1e-4;
102 : 55 : const double beta0 = 0.25;
103 : 55 : const double beta1 = 0.80;
104 : 55 : const double tol1 = 1e-8;
105 : 55 : const double tol2 = 1e-12;
106 : 55 : const double epsilon = 1e-10;
107 : : double original_value, new_value;
108 : : double beta;
109 : :
110 [ + - ]: 55 : int nv = pd.num_free_vertices();
111 [ + - ][ + - ]: 110 : std::vector< Vector3D > grad( nv ), d( nv );
[ + - ]
112 : 55 : bool fn_bool = true; // bool used for determining validity of patch
113 : :
114 [ + - ]: 55 : OFEvaluator& objFunc = get_objective_function_evaluator();
115 : :
116 : : int i;
117 : :
118 : : // TODD -- Don't blame the code for bad things happening when using a
119 : : // bad termination test or requesting more accuracy than is
120 : : // possible.
121 : : //
122 : : // Also,
123 : :
124 : : // 1. Allocate a hessian and calculate the sparsity pattern.
125 [ + - ][ + - ]: 55 : mHessian.initialize( pd, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
126 : :
127 : : // does the Feasible Newton iteration until stopping is required.
128 : : // Terminate when inner termination criterion signals.
129 : :
130 : : /* Computes the value of the stopping criterion*/
131 [ + - ]: 55 : TerminationCriterion* term_crit = get_inner_termination_criterion();
132 [ + - ][ + + ]: 374 : while( !term_crit->terminate() )
[ + - ]
133 : : {
134 [ + - ][ + - ]: 319 : fn_bool = objFunc.update( pd, original_value, grad, mHessian, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
135 [ - + ]: 319 : if( !fn_bool )
136 : : {
137 : : MSQ_SETERR( err )
138 [ # # ][ # # ]: 0 : ( "invalid patch for hessian calculation", MsqError::INTERNAL_ERROR );
139 : 0 : return;
140 : : }
141 : :
142 : : if( MSQ_DBG( 3 ) )
143 : : { // avoid expensive norm calculations if debug flag is off
144 : : MSQ_DBGOUT( 3 ) << " o objective function: " << original_value << std::endl;
145 : : MSQ_DBGOUT( 3 ) << " o gradient norm: " << length( grad ) << std::endl;
146 : : MSQ_DBGOUT( 3 ) << " o Hessian norm: " << mHessian.norm() << std::endl;
147 : : }
148 : :
149 : : // Prints out free vertices coordinates.
150 : : //
151 : : // Comment out the following because it is way to verbose for larger
152 : : // problems. Consider doing:
153 : : // inner_term_crit->write_mesh_steps( "filename.vtk" );
154 : : // instead.
155 : : // - j.kraftcheck 2010-11-17
156 : : // if (MSQ_DBG(3)) {
157 : : // MSQ_DBGOUT(3) << "\n o Free vertices ("<< pd.num_free_vertices()
158 : : // <<")original coordinates:\n ";
159 : : // MSQ_ERRRTN(err);
160 : : // const MsqVertex* toto1 = pd.get_vertex_array(err); MSQ_ERRRTN(err);
161 : : // MsqFreeVertexIndexIterator ind1(pd, err); MSQ_ERRRTN(err);
162 : : // ind1.reset();
163 : : // while (ind1.next()) {
164 : : // MSQ_DBGOUT(3) << "\t\t\t" << toto1[ind1.value()];
165 : : // }
166 : : // }
167 : :
168 : : // 4. Calculate a direction using preconditionned conjugate gradients
169 : : // to find a zero of the Newton system of equations (H*d = -g)
170 : : // (a) stop if conjugate iteration limit reached
171 : : // (b) stop if relative residual is small
172 : : // (c) stop if direction of negative curvature is obtained
173 : :
174 [ + - ][ + - ]: 319 : mHessian.cg_solver( arrptr( d ), arrptr( grad ), err );MSQ_ERRRTN( err );
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ - + ]
175 : :
176 : : // 5. Check for descent direction (inner produce of gradient and
177 : : // direction is negative.
178 [ + - ]: 319 : double alpha = inner( grad, d );
179 : : // TODD -- Add back in if you encounter problems -- do a gradient
180 : : // step if the direction from the conjugate gradient solver
181 : : // is not a descent direction for the objective function. We
182 : : // SHOULD always get a descent direction from the conjugate
183 : : // method though, unless the preconditioner is not positive
184 : : // definite.
185 : : // If direction is positive, does a gradient (steepest descent) step.
186 : :
187 [ + + ]: 319 : if( alpha > -epsilon )
188 : : {
189 : :
190 : : MSQ_DBGOUT( 3 ) << " o alpha = " << alpha << " (rejected)" << std::endl;
191 : :
192 [ + + ]: 93 : if( !havePrintedDirectionMessage )
193 : : {
194 : : MSQ_PRINT( 1 )
195 : : ( "Newton direction not guaranteed descent. Ensure preconditioner is positive "
196 : : "definite.\n" );
197 : 44 : havePrintedDirectionMessage = true;
198 : : }
199 : :
200 : : // TODD: removed performing gradient step here since we will use
201 : : // gradient if step does not produce descent. Instead we set
202 : : // alpha to a small negative value.
203 : :
204 : 93 : alpha = -epsilon;
205 : :
206 : : // alpha = inner(grad, grad, nv); // compute norm squared of gradient
207 : : // if (alpha < 1) alpha = 1; // take max with constant
208 : : // for (i = 0; i < nv; ++i) {
209 : : // d[i] = -grad[i] / alpha; // compute scaled gradient
210 : : // }
211 : : // alpha = inner(grad, d, nv); // recompute alpha
212 : : // // equal to one for large gradient
213 : : }
214 : : else
215 : : {
216 : : MSQ_DBGOUT( 3 ) << " o alpha = " << alpha << std::endl;
217 : : }
218 : :
219 : 319 : alpha *= sigma;
220 : 319 : beta = 1.0;
221 [ + - ][ + - ]: 319 : pd.recreate_vertices_memento( coordsMem, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
222 : :
223 : : // TODD: Unrolling the linesearch loop. We do a function and
224 : : // gradient evaluation when beta = 1. Otherwise, we end up
225 : : // in the linesearch regime. We expect that several
226 : : // evaluations will be made, so we only do a function evaluation
227 : : // and finish with a gradient evaluation. When beta = 1, we also
228 : : // check the gradient for stability.
229 : :
230 : : // TODD -- the Armijo linesearch is based on the objective function,
231 : : // so theoretically we only need to evaluate the objective
232 : : // function. However, near a very accurate solution, say with
233 : : // the two norm of the gradient of the objective function less
234 : : // than 1e-5, the numerical error in the objective function
235 : : // calculation is enough that the Armijo linesearch will
236 : : // fail. To correct this situation, the iterate is accepted
237 : : // when the norm of the gradient is also small. If you need
238 : : // high accuracy and have a large mesh, talk with Todd about
239 : : // the numerical issues so that we can fix it.
240 : :
241 : : // TODD -- the Armijo linesearch here is only good for differentiable
242 : : // functions. When projections are involved, you should change
243 : : // to a form of the linesearch meant for nondifferentiable
244 : : // functions.
245 : :
246 [ + - ][ + - ]: 319 : pd.move_free_vertices_constrained( arrptr( d ), nv, beta, err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
247 [ + - ]: 319 : fn_bool = objFunc.evaluate( pd, new_value, grad, err );
248 [ + - ][ + + ]: 319 : if( err.error_code() == err.BARRIER_VIOLATED )
249 [ + - ]: 20 : err.clear(); // barrier violated does not represent an actual error here
250 [ + - ][ - + ]: 319 : MSQ_ERRRTN( err );
[ # # ][ # # ]
[ - + ]
251 : :
252 [ + + ][ + + ]: 388 : if( ( fn_bool && ( original_value - new_value >= -alpha * beta - epsilon ) ) ||
[ + + ][ + + ]
253 [ + - ][ + - ]: 69 : ( fn_bool && ( length( arrptr( grad ), nv ) < 100 * convTol ) ) )
[ + + ]
254 : : {
255 : : // Armijo linesearch rules passed.
256 : : MSQ_DBGOUT( 3 ) << " o beta = " << beta << " (accepted without line search)" << std::endl;
257 : : }
258 : : else
259 : : {
260 [ + + ]: 90 : if( !fn_bool )
261 : : {
262 : : // Function undefined. Use the higher decrease rate.
263 : 23 : beta *= beta0;
264 : : MSQ_DBGOUT( 3 ) << " o beta = " << beta << " (invalid step)" << std::endl;
265 : : }
266 : : else
267 : : {
268 : : // Function defined, but not sufficient decrease
269 : : // Use the lower decrease rate.
270 : 67 : beta *= beta1;
271 : : MSQ_DBGOUT( 3 ) << " o beta = " << beta << " (insufficient decrease)" << std::endl;
272 : : }
273 [ + - ][ + - ]: 90 : pd.set_to_vertices_memento( coordsMem, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
274 : :
275 : : // Standard Armijo linesearch rules
276 : :
277 : : MSQ_DBGOUT( 3 ) << " o Doing line search" << std::endl;
278 [ + + ]: 1680 : while( beta >= tol1 )
279 : : {
280 : : // 6. Search along the direction
281 : : // (a) trial = x + beta*d
282 [ + - ][ + - ]: 1677 : pd.move_free_vertices_constrained( arrptr( d ), nv, beta, err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
283 : : // (b) function evaluation
284 [ + - ]: 1677 : fn_bool = objFunc.evaluate( pd, new_value, err );
285 [ + - ][ + + ]: 1677 : if( err.error_code() == err.BARRIER_VIOLATED )
286 [ + - ]: 22 : err.clear(); // barrier violated does not represent an actual error here
287 [ + - ][ - + ]: 1677 : MSQ_ERRRTN( err );
[ # # ][ # # ]
[ - + ]
288 : :
289 : : // (c) check for sufficient decrease and stop
290 [ + + ]: 1677 : if( !fn_bool )
291 : : {
292 : : // function not defined at trial point
293 : 22 : beta *= beta0;
294 : : }
295 [ + + ]: 1655 : else if( original_value - new_value >= -alpha * beta - epsilon )
296 : : {
297 : : // iterate is acceptable.
298 : 87 : break;
299 : : }
300 : : else
301 : : {
302 : : // iterate is not acceptable -- shrink beta
303 : 1568 : beta *= beta1;
304 : : }
305 [ + - ][ + - ]: 1590 : pd.set_to_vertices_memento( coordsMem, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
306 : : }
307 : :
308 [ + + ]: 90 : if( beta < tol1 )
309 : : {
310 : : // assert(pd.set_to_vertices_memento called last)
311 : :
312 : : // TODD -- Lower limit on steplength reached. Direction does not
313 : : // appear to make sufficient progress decreasing the
314 : : // objective function. This can happen when you are
315 : : // very near a solution due to numerical errors in
316 : : // computing the objective function. It can also happen
317 : : // when the direction is not a descent direction and when
318 : : // you are projecting the iterates onto a surface.
319 : : //
320 : : // The latter cases require the use of a linesearch on
321 : : // a gradient step. If this linesearch terminate with
322 : : // insufficient decrease, then you are at a critical
323 : : // point and should stop!
324 : : //
325 : : // The numerical errors with the objective function cannot
326 : : // be overcome. Eventually, the gradient step will
327 : : // fail to compute a new direction and you will stop.
328 : :
329 : : MSQ_PRINT( 1 )
330 : : ( "Sufficient decrease not obtained in linesearch; switching to gradient.\n" );
331 : :
332 [ + - ][ + - ]: 3 : alpha = inner( arrptr( grad ), arrptr( grad ),
333 [ + - ]: 3 : nv ); // compute norm squared of gradient
334 [ + + ]: 3 : if( alpha < 1 ) alpha = 1; // take max with constant
335 [ + + ]: 132 : for( i = 0; i < nv; ++i )
336 : : {
337 [ + - ][ + - ]: 129 : d[i] = -grad[i] / alpha; // compute scaled gradient
[ + - ][ + - ]
[ + - ]
338 : : }
339 [ + - ][ + - ]: 3 : alpha = inner( arrptr( grad ), arrptr( d ), nv ); // recompute alpha
[ + - ]
340 : 3 : alpha *= sigma; // equal to one for large gradient
341 : 3 : beta = 1.0;
342 : :
343 : : // Standard Armijo linesearch rules
344 [ + - ]: 85 : while( beta >= tol2 )
345 : : {
346 : : // 6. Search along the direction
347 : : // (a) trial = x + beta*d
348 [ + - ][ + - ]: 85 : pd.move_free_vertices_constrained( arrptr( d ), nv, beta, err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
349 : : // (b) function evaluation
350 [ + - ]: 85 : fn_bool = objFunc.evaluate( pd, new_value, err );
351 [ + - ][ - + ]: 85 : if( err.error_code() == err.BARRIER_VIOLATED )
352 [ # # ]: 0 : err.clear(); // barrier violated does not represent an actual error
353 : : // here
354 [ + - ][ - + ]: 85 : MSQ_ERRRTN( err );
[ # # ][ # # ]
[ - + ]
355 : :
356 : : // (c) check for sufficient decrease and stop
357 [ - + ]: 85 : if( !fn_bool )
358 : : {
359 : : // function not defined at trial point
360 : 0 : beta *= beta0;
361 : : }
362 [ + + ]: 85 : else if( original_value - new_value >= -alpha * beta - epsilon )
363 : : {
364 : : // iterate is acceptable.
365 : 3 : break;
366 : : }
367 : : else
368 : : {
369 : : // iterate is not acceptable -- shrink beta
370 : 82 : beta *= beta1;
371 : : }
372 [ + - ][ + - ]: 82 : pd.set_to_vertices_memento( coordsMem, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
373 : : }
374 : :
375 [ - + ]: 3 : if( beta < tol2 )
376 : : {
377 : : // assert(pd.set_to_vertices_memento called last)
378 : :
379 : : // TODD -- Lower limit on steplength reached. Gradient does not
380 : : // appear to make sufficient progress decreasing the
381 : : // objective function. This can happen when you are
382 : : // very near a solution due to numerical errors in
383 : : // computing the objective function. Most likely you
384 : : // are at a critical point for the problem.
385 : :
386 : : MSQ_PRINT( 1 )
387 : : ( "Sufficient decrease not obtained with gradient; critical point likely "
388 : : "found.\n" );
389 : 0 : break;
390 : : }
391 : : }
392 : :
393 : : // Compute the gradient at the new point -- needed by termination check
394 [ + - ][ + - ]: 90 : fn_bool = objFunc.update( pd, new_value, grad, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
395 : : }
396 : :
397 : : // Prints out free vertices coordinates.
398 : : // if (MSQ_DBG(3)) {
399 : : // MSQ_DBGOUT(3) << " o Free vertices new coordinates: \n";
400 : : // const MsqVertex* toto1 = pd.get_vertex_array(err); MSQ_ERRRTN(err);
401 : : // MsqFreeVertexIndexIterator ind(pd, err); MSQ_ERRRTN(err);
402 : : // ind.reset();
403 : : // while (ind.next()) {
404 : : // MSQ_DBGOUT(3) << "\t\t\t" << toto1[ind.value()];
405 : : // }
406 : : // }
407 : :
408 : : // checks stopping criterion
409 [ + - ][ + - ]: 319 : term_crit->accumulate_patch( pd, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
410 [ + - ][ + - ]: 319 : term_crit->accumulate_inner( pd, new_value, arrptr( grad ), err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
411 : : }
412 : 55 : MSQ_PRINT( 2 )( "FINISHED\n" );
413 : : }
414 : : else
415 : : {
416 : 1 : std::cout << "WARNING: Feasible Newton optimization only supported for volume meshes" << std::endl
417 : 1 : << " and XYPlanarDomain surface meshes." << std::endl
418 : 1 : << std::endl
419 : 56 : << "Try a different solver such as Steepest Descent." << std::endl;
420 : : }
421 : : }
422 : :
423 : 56 : void FeasibleNewton::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ )
424 : : {
425 : :
426 : : // Michael:: Should the vertices memento be delete here???
427 : : // cout << "- Executing FeasibleNewton::iteration_complete()\n";
428 : 56 : }
429 : :
430 : 56 : void FeasibleNewton::cleanup()
431 : : {
432 [ + - ]: 56 : delete coordsMem;
433 : 56 : coordsMem = NULL;
434 [ + - ][ + - ]: 40 : }
|
