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 : : \file NonGradient.cpp
29 : : \brief
30 : : The NonGradient class is also a concrete vertex mover
31 : : which performs derivative free minimization
32 : : based on the Amoeba Method, as implemented in Numerical Recipes in C.
33 : :
34 : : The other optimization methods either accept or reject a point immediately.
35 : : Amoeba may not accept a point for several iterations/evaluations.
36 : : Amoeba does not check for BOUNDED_VERTEX_MOVERMENT.
37 : : rtol=2.0*fabs( height[ihi]-height[ilo] )/
38 : : ( fabs(height[ihi])+fabs(height[ilo])+threshold );
39 : :
40 : : We are aware of the problem of choosing the initial simplex, but have not
41 : : sovled it. How has this been solved?
42 : : There is a problem of letting Mesquite handle convergence. The idea here is to
43 : : call accumulate_? only if the best point changes. But does it matter? Where
44 : : are the errors determined? What does it mean to cull?
45 : : What does Pat... compare and contrast the differnt vertex movement criteria
46 : : with the Amoeba criterion.
47 : : */
48 : :
49 : : #include "NonGradient.hpp"
50 : : #include "MsqFreeVertexIndexIterator.hpp"
51 : : #include "MsqTimer.hpp"
52 : : #include "MsqDebug.hpp"
53 : : #include "MsqError.hpp"
54 : : #include <cmath>
55 : : #include <iostream>
56 : : #include "../../ObjectiveFunction/ObjectiveFunction.hpp"
57 : : #include <float.h>
58 : :
59 : : namespace MBMesquite
60 : : {
61 : :
62 : 0 : std::string NonGradient::get_name() const
63 : : {
64 [ # # ]: 0 : return "NonGradient";
65 : : }
66 : :
67 : 4 : PatchSet* NonGradient::get_patch_set()
68 : : {
69 : 4 : return PatchSetUser::get_patch_set();
70 : : }
71 : :
72 : 3 : NonGradient::NonGradient( ObjectiveFunction* of )
73 : : : VertexMover( of ), PatchSetUser( true ), projectGradient( false ), mDimension( 0 ), mThreshold( 0.0 ),
74 [ + - ][ + - ]: 3 : mTolerance( 0.0 ), mMaxNumEval( 0 ), mNonGradDebug( 0 ), mUseExactPenaltyFunction( true ), mScaleDiameter( 0.1 )
[ + - ]
75 : : {
76 [ + - ]: 3 : set_debugging_level( 2 );
77 : : // set the default inner termination criterion
78 [ + - ]: 3 : TerminationCriterion* default_crit = get_inner_termination_criterion();
79 [ - + ]: 6 : if( default_crit == NULL ) { return; }
80 : : else
81 : : {
82 [ + - ]: 3 : default_crit->add_iteration_limit( 5 );
83 : : }
84 : : }
85 : :
86 : 0 : NonGradient::NonGradient( ObjectiveFunction* of, MsqError& err )
87 : : : VertexMover( of ), PatchSetUser( true ), projectGradient( false ), mDimension( 0 ), mThreshold( 0.0 ),
88 [ # # ][ # # ]: 0 : mTolerance( 0.0 ), mMaxNumEval( 0 ), mNonGradDebug( 0 ), mUseExactPenaltyFunction( true ), mScaleDiameter( 0.1 )
[ # # ]
89 : : {
90 [ # # ]: 0 : set_debugging_level( 2 );
91 : : // set the default inner termination criterion
92 [ # # ]: 0 : TerminationCriterion* default_crit = get_inner_termination_criterion();
93 [ # # ]: 0 : if( default_crit == NULL )
94 : : {
95 : : MSQ_SETERR( err )
96 : : ( "QualityImprover did not create a default inner "
97 : : "termination criterion.",
98 [ # # ][ # # ]: 0 : MsqError::INVALID_STATE );
99 : 0 : return;
100 : : }
101 : : else
102 : : {
103 [ # # ][ # # ]: 0 : default_crit->add_iteration_limit( 5 );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
104 : : }
105 : : }
106 : :
107 : 0 : bool NonGradient::testRowSum( int numRow, int numCol, double* matrix, double* oldRowSum )
108 : : {
109 : 0 : bool same = true;
110 [ # # ]: 0 : std::vector< double > rowSum( numRow, 0. );
111 : 0 : double maxVal = 0.;
112 [ # # ]: 0 : for( int col = 0; col < numCol; col++ )
113 : : {
114 [ # # ]: 0 : for( int row = 0; row < numRow; row++ )
115 : : {
116 [ # # ]: 0 : rowSum[row] += matrix[row + col * numRow];
117 [ # # ]: 0 : if( fabs( matrix[row + col * numRow] ) > maxVal ) { maxVal = fabs( matrix[row + col * numRow] ); }
118 : : }
119 : : }
120 : 0 : double machEps = 1.e-14 * static_cast< double >( numRow ); // better to use system parameters
121 : 0 : double upperBound = machEps * maxVal + 1.e-15;
122 [ # # ]: 0 : for( int row = 0; row < numRow; row++ )
123 : : {
124 [ # # ][ # # ]: 0 : if( fabs( rowSum[row] - oldRowSum[row] ) > upperBound )
125 : : {
126 : 0 : same = false;
127 [ # # ]: 0 : if( mNonGradDebug >= 2 )
128 : : {
129 [ # # ][ # # ]: 0 : std::cout << " Warning: NonGradient Row Sum " << row << " Test: value " << rowSum[row]
[ # # ][ # # ]
[ # # ]
130 [ # # ][ # # ]: 0 : << " Discrepancy " << rowSum[row] - oldRowSum[row] << " maxVal " << maxVal << " numRow "
[ # # ][ # # ]
[ # # ][ # # ]
131 [ # # ][ # # ]: 0 : << numRow << " numCol " << numCol << std::endl;
[ # # ][ # # ]
132 : : }
133 : : MSQ_PRINT( 2 )
134 : : ( "NonGradient Row Sum [%d] Test failure: value %22.15e difference %22.15e \n", row, rowSum[row],
135 : : rowSum[row] - oldRowSum[row] );
136 : : }
137 : : }
138 : 0 : return ( same );
139 : : }
140 : :
141 : 850 : void NonGradient::getRowSum( int numRow, int numCol, std::vector< double >& matrix, std::vector< double >& rowSum )
142 : : {
143 [ + + ]: 2550 : for( int row = 0; row < numRow; row++ )
144 : : {
145 : 1700 : rowSum[row] = 0.;
146 : : }
147 [ + + ]: 3400 : for( int col = 0; col < numCol; col++ )
148 : : {
149 [ + + ]: 7650 : for( int row = 0; row < numRow; row++ )
150 : : {
151 : 5100 : rowSum[row] += matrix[row + col * numRow];
152 : : }
153 : : }
154 : 850 : }
155 : :
156 : 19860 : double NonGradient::evaluate( double* point, PatchData& pd, MsqError& err )
157 : : {
158 [ + - ][ - + ]: 19860 : if( pd.num_free_vertices() > 1 )
159 : : {
160 : : MSQ_SETERR( err )
161 [ # # ][ # # ]: 0 : ( "Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED );
162 : : }
163 : 19860 : const size_t vertexIndex = 0;
164 [ + - ][ + - ]: 19860 : Vector3D originalVec = pd.vertex_by_index( vertexIndex );
165 [ + - ]: 19860 : Vector3D pointVec;
166 [ + + ]: 79440 : for( int index = 0; index < 3; index++ )
167 : : {
168 [ + - ]: 59580 : pointVec[index] = point[index];
169 : : }
170 [ + - ]: 19860 : pd.set_vertex_coordinates( pointVec, vertexIndex, err );
171 [ + - ]: 19860 : pd.snap_vertex_to_domain( vertexIndex, err );
172 : :
173 [ + - ]: 19860 : OFEvaluator& obj_func = get_objective_function_evaluator();
174 [ + - ]: 19860 : TerminationCriterion* term_crit = get_inner_termination_criterion();
175 : :
176 : : double value;
177 [ + - ]: 19860 : bool feasible = obj_func.evaluate( pd, value, err ); // MSQ_ERRZERO(err);
178 [ + - ]: 19860 : term_crit->accumulate_inner( pd, value, 0, err ); // MSQ_CHKERR(err);
179 [ + - ][ + + ]: 19860 : if( err.error_code() == err.BARRIER_VIOLATED ) err.clear(); // barrier violation not an error here
[ + - ]
180 [ + - ][ - + ]: 19860 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
181 [ + - ]: 19860 : pd.set_vertex_coordinates( originalVec, vertexIndex, err );
182 [ + + ][ + - ]: 19860 : if( !feasible && mUseExactPenaltyFunction )
183 : : { // "value" undefined btw
184 : 1127 : double ensureFiniteRtol = .25;
185 : 1127 : value = DBL_MAX * ensureFiniteRtol;
186 : : // std::cout << " Infeasible patch: " << value << std::endl; printPatch( pd, err );
187 : : }
188 : 19860 : return value;
189 : : }
190 : :
191 : : // Extrapolate by a factor of fac through the face of the simplex
192 : : // opposite the high point to a new point. If the new point is
193 : : // better, swap it with the high point.
194 : 15294 : double NonGradient::amotry( std::vector< double >& p_simplex, std::vector< double >& p_height, double psum[], int ihi,
195 : : double fac, PatchData& pd, MsqError& err )
196 : : {
197 [ + - ]: 15294 : int numRow = getDimension();
198 : : // int numCol = numRow + 1;
199 [ + - ]: 15294 : std::vector< double > ptry( numRow ); // does this make sense?
200 : 15294 : double fac1 = ( 1.0 - fac ) / static_cast< double >( numRow );
201 : 15294 : double fac2 = fac1 - fac;
202 [ + + ]: 45882 : for( int row = 0; row < numRow; row++ )
203 : : {
204 [ + - ][ + - ]: 30588 : ptry[row] = psum[row] * fac1 - p_simplex[row + ihi * numRow] * fac2;
205 : : }
206 [ - + ][ # # ]: 15294 : if( mNonGradDebug >= 3 ) { std::cout << "Try "; }
207 : : MSQ_PRINT( 3 )( "Try" );
208 [ + - ][ + - ]: 15294 : double ytry = evaluate( &ptry[0], pd, err ); // value at trial point
209 [ - + ][ # # ]: 15294 : if( mNonGradDebug >= 3 ) { std::cout << ytry << std::endl; }
[ # # ]
210 : : MSQ_PRINT( 3 )( "yTry" );
211 [ + - ][ + + ]: 15294 : if( ytry < p_height[ihi] ) // better than highest (worst)
212 : : {
213 [ + - ]: 3115 : p_height[ihi] = ytry; // swap ihi and ytry
214 [ + + ]: 9345 : for( int row = 0; row < numRow; row++ )
215 : : {
216 [ + - ][ + - ]: 6230 : psum[row] += ( ptry[row] - p_simplex[row + ihi * numRow] );
217 [ + - ][ + - ]: 6230 : p_simplex[row + ihi * numRow] = ptry[row];
218 : : }
219 : : }
220 : 15294 : return ytry;
221 : : }
222 : :
223 : 0 : void NonGradient::printPatch( const PatchData& pd, MsqError& err )
224 : : {
225 [ # # ]: 0 : if( mNonGradDebug == 0 ) { return; }
226 : 0 : const size_t numNode = pd.num_nodes(); // 27, 27 what?
227 : : MSQ_PRINT( 3 )( "Number of Vertices: %d\n", (int)pd.num_nodes() );
228 : 0 : const size_t numVert = pd.num_free_vertices(); // 1
229 : : MSQ_PRINT( 3 )( "Num Free = %d\n", (int)pd.num_free_vertices() );
230 : 0 : const size_t numSlaveVert = pd.num_slave_vertices(); // 0
231 : 0 : const size_t numCoin = pd.num_corners(); // 64
232 : 0 : const MsqVertex* coord = pd.get_vertex_array( err );
233 : : MSQ_PRINT( 3 )( "Number of Vertices: %d\n", (int)pd.num_nodes() );
234 : :
235 : 0 : std::cout << "Patch " << numNode << " " << numVert << " " << numSlaveVert << " " << numCoin << std::endl;
236 : : MSQ_PRINT( 3 )( "\n" );
237 : 0 : std::cout << "Coordinate ";
238 : 0 : std::cout << " " << std::endl;
239 [ # # ]: 0 : for( size_t index = 0; index < numVert; index++ )
240 : : {
241 : 0 : std::cout << coord[index][0] << " " << coord[index][1] << " " << coord[index][2] << std::endl;
242 : : }
243 : : // const size_t numElt = pd.num_elements();
244 [ # # ]: 0 : if( mNonGradDebug >= 3 ) { std::cout << "Number of Elements: " << pd.num_elements() << std::endl; }
245 : : MSQ_PRINT( 3 )( "Number of Elements: %d\n", (int)pd.num_elements() );
246 : : }
247 : :
248 : 0 : void NonGradient::printSimplex( std::vector< double >& p_simplex, std::vector< double >& p_height )
249 : : {
250 : 0 : int numRow = getDimension();
251 : 0 : int numCol = numRow + 1;
252 [ # # ]: 0 : for( int col = 0; col < numCol; col++ )
253 : : {
254 : : // std::cout << "simplex[ " << col << "]= " ;
255 [ # # ]: 0 : for( int row = 0; row < numRow; row++ )
256 : : {
257 : 0 : std::cout << p_simplex[row + col * numRow] << " ";
258 : : }
259 : : // std::cout << " " << p_height[col] << std::endl;
260 : : }
261 [ # # ]: 0 : for( int col = 0; col < numCol; col++ )
262 : : {
263 : 0 : std::cout << p_height[col] << " ";
264 : : }
265 : 0 : std::cout << std::endl;
266 : 0 : }
267 : :
268 : 850 : int NonGradient::getPatchDimension( const PatchData& pd, MsqError& err )
269 : : {
270 : 850 : const size_t numElt = pd.num_elements();
271 : 850 : unsigned edimMax = 0; // in case numElt == 0
272 [ + + ]: 5850 : for( size_t elementId = 0; elementId < numElt; elementId++ )
273 : : {
274 : 5000 : const MsqMeshEntity& element = pd.element_by_index( elementId );
275 : 5000 : EntityTopology type = element.get_element_type();
276 : 5000 : unsigned edim = TopologyInfo::dimension( type );
277 [ + + ]: 5000 : if( elementId == 0 ) { edimMax = edim; }
278 : : else
279 : : {
280 [ - + ]: 4150 : if( edimMax != edim )
281 : : {
282 : : MSQ_SETERR( err )
283 [ # # ]: 0 : ( "A patch has elements of different dimensions", MsqError::INVALID_MESH );
284 : 0 : std::cout << "A patch has elements of different dimensions" << edimMax << " " << edim << std::endl;
285 [ # # ]: 0 : if( edimMax < edim ) { edimMax = edim; }
286 : : }
287 : : }
288 : : }
289 : 850 : return ( edimMax );
290 : : }
291 : :
292 : 5 : void NonGradient::initialize( PatchData& /*pd*/, MsqError& /*err*/ ) {}
293 : :
294 : 850 : void NonGradient::initialize_mesh_iteration( PatchData& pd, MsqError& err )
295 : : {
296 [ + - ]: 850 : unsigned elementDimension = getPatchDimension( pd, err ); // to do: react to error
297 [ + - ]: 850 : unsigned dimension = elementDimension * pd.num_free_vertices();
298 : : // printPatch( pd, err );
299 [ + - ]: 850 : setDimension( dimension );
300 : 850 : int maxNumEval = 100 * dimension; // 1. Make this a user parameter
301 [ + - ]: 850 : setMaxNumEval( maxNumEval );
302 : 850 : double threshold = 1.e-10; // avoid division by zero
303 [ + - ]: 850 : setThreshold( threshold );
304 : 850 : double minEdgeLen = 0.0;
305 : 850 : double maxEdgeLen = 0.0;
306 : : // double ftol = 0.;
307 [ + - ]: 850 : if( dimension > 0 )
308 : : {
309 [ + - ]: 850 : pd.get_minmax_edge_length( minEdgeLen, maxEdgeLen );
310 : : // ftol = minEdgeLen * 1.e-4; // Turn off Amoeba convergence criterion
311 [ - + ]: 850 : if( mNonGradDebug >= 1 )
312 [ # # ][ # # ]: 0 : { std::cout << "minimum edge length " << minEdgeLen << " maximum edge length " << maxEdgeLen << std::endl; }
[ # # ][ # # ]
[ # # ]
313 : : MSQ_PRINT( 3 )
314 : : ( "minimum edge length %e maximum edge length %e\n", minEdgeLen, maxEdgeLen );
315 : : }
316 : : // setTolerance(ftol);
317 : 850 : unsigned numRow = dimension;
318 : 850 : unsigned numCol = numRow + 1;
319 [ + - ]: 850 : if( numRow * numCol <= simplex.max_size() )
320 : : {
321 [ + - ]: 850 : simplex.assign( numRow * numCol, 0. ); // guard against previous simplex value
322 : 850 : double scale = minEdgeLen * mScaleDiameter;
323 : : ;
324 [ + - ]: 850 : const MsqVertex* coord = pd.get_vertex_array( err );
325 [ + - ][ - + ]: 850 : if( pd.num_free_vertices() > 1 )
326 : : {
327 : : MSQ_SETERR( err )
328 [ # # ][ # # ]: 0 : ( "Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED );
329 : : }
330 : 850 : size_t index = 0;
331 [ + + ]: 3400 : for( unsigned col = 0; col < numCol; col++ )
332 : : {
333 [ + + ]: 7650 : for( unsigned row = 0; row < numRow; row++ )
334 : : {
335 [ + - ][ + - ]: 5100 : simplex[row + col * numRow] = coord[index][row];
336 [ + + ][ + - ]: 5100 : if( row == col - 1 ) { simplex[row + col * numRow] += scale / static_cast< double >( numCol ); }
337 : : }
338 : : }
339 : : }
340 : : else
341 : : {
342 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "Use patch with fewer free vertices", MsqError::OUT_OF_MEMORY );
343 [ # # ][ # # ]: 0 : if( mNonGradDebug >= 1 ) { std::cout << "ERROR: Too many free vertices in patch" << std::endl; }
[ # # ]
344 : : // MSQ_PRINT(1)("ERROR: Too many free vertices in patch\n");
345 : : }
346 : 850 : }
347 : :
348 : 850 : void NonGradient::optimize_vertex_positions( PatchData& pd, MsqError& err )
349 : : {
350 : : MSQ_FUNCTION_TIMER( "NonGradient::optimize_vertex_positions" );
351 [ + - ]: 850 : int numRow = getDimension();
352 : 850 : int numCol = numRow + 1;
353 [ + - ]: 850 : std::vector< double > p_height( numCol );
354 : :
355 [ + + ]: 3400 : for( int col = 0; col < numCol; col++ )
356 : : {
357 [ + - ][ + - ]: 2550 : p_height[col] = evaluate( &simplex[col * numRow], pd, err ); // eval patch stuff
[ + - ]
358 : : }
359 [ - + ][ # # ]: 850 : if( mNonGradDebug > 0 ) { printSimplex( simplex, p_height ); }
360 : :
361 : : // standardization
362 [ + - ]: 850 : TerminationCriterion* term_crit = get_inner_termination_criterion();
363 : : // int maxNumEval = getMaxNumEval();
364 : : // double threshold = getThreshold();
365 : :
366 : : // double ftol = getTolerance();
367 : 850 : int ilo = 0; // height[ilo]<=...
368 : 850 : int inhi = 0; //...<=height[inhi]<=
369 : 850 : int ihi = 0; //<=height[ihi]
370 : : // double rtol = 2.*ftol;
371 : : double ysave;
372 : : double ytry;
373 : 850 : bool afterEvaluation = false;
374 [ + - ]: 1700 : std::vector< double > rowSum( numRow );
375 [ + - ]: 850 : getRowSum( numRow, numCol, simplex, rowSum );
376 [ + - ][ + + ]: 9536 : while( !( term_crit->terminate() ) )
377 : : {
378 : :
379 [ - + ][ # # ]: 8686 : if( mNonGradDebug > 0 ) { printSimplex( simplex, p_height ); }
380 : : // std::cout << "rtol " << rtol << " ftol " << ftol << " MesquiteIter " <<
381 : : // term_crit->get_iteration_count() << " Done " << term_crit->terminate() << std::endl;
382 : :
383 [ + + ]: 8686 : if( afterEvaluation )
384 : : {
385 : : // reflect highPt through opposite face
386 : : // p_height[0] may vanish
387 : : /*
388 : : if( !testRowSum( numRow, numCol, &simplex[0], &rowSum[0]) )
389 : : {
390 : : // Before uncommenting here and ...
391 : : //MSQ_SETERR(err)("Internal check sum test A failed", MsqError::INTERNAL_ERROR);
392 : : //MSQ_ERRRTN(err);
393 : : }
394 : : */
395 [ + - ][ + - ]: 7836 : ytry = amotry( simplex, p_height, &rowSum[0], ihi, -1.0, pd, err );
396 : : /*
397 : : if( !testRowSum( numRow, numCol, &simplex[0], &rowSum[0]) )
398 : : {
399 : : // ... here, determine a maxVal majorizing the previous as well as the current
400 : : simplex.
401 : : //MSQ_SETERR(err)("Internal check sum test B failed",
402 : : MsqError::INTERNAL_ERROR);
403 : : //MSQ_ERRRTN(err);
404 : : }
405 : : */
406 : :
407 : : /*
408 : : if( p_height[0] == 0.)
409 : : {
410 : : MSQ_SETERR(err)("(B) Zero objective function value", MsqError::INTERNAL_ERROR);
411 : : exit(-1);
412 : : }
413 : : */
414 [ + - ][ + + ]: 7836 : if( ytry <= p_height[ilo] )
415 : : {
416 [ + - ][ + - ]: 4542 : ytry = amotry( simplex, p_height, &rowSum[0], ihi, -2.0, pd, err );
417 [ - + ][ # # ]: 4542 : if( mNonGradDebug >= 3 ) { std::cout << "Reflect and Expand from highPt " << ytry << std::endl; }
[ # # ][ # # ]
418 : : // MSQ_PRINT(3)("Reflect and Expand from highPt : %e\n",ytry);
419 : : }
420 : : else
421 : : {
422 [ + - ][ + + ]: 3294 : if( ytry >= p_height[inhi] )
423 : : {
424 [ + - ]: 2916 : ysave = p_height[ihi]; // Contract along highPt
425 [ + - ][ + - ]: 2916 : ytry = amotry( simplex, p_height, &rowSum[0], ihi, 0.5, pd, err );
426 [ + + ]: 2916 : if( ytry >= ysave )
427 : : { // contract all directions toward lowPt
428 [ + + ]: 10860 : for( int col = 0; col < numCol; col++ )
429 : : {
430 [ + + ]: 3024 : if( col != ilo )
431 : : {
432 [ + + ]: 6048 : for( int row = 0; row < numRow; row++ )
433 : : {
434 [ + - ][ + - ]: 4032 : rowSum[row] = 0.5 * ( simplex[row + col * numRow] + simplex[row + ilo * numRow] );
[ + - ]
435 [ + - ][ + - ]: 4032 : simplex[row + col * numRow] = rowSum[row];
436 : : }
437 [ + - ][ + - ]: 2016 : p_height[col] = evaluate( &rowSum[0], pd, err );
[ + - ]
438 [ - + ]: 2016 : if( mNonGradDebug >= 3 )
439 : : {
440 [ # # ][ # # ]: 0 : std::cout << "Contract all directions toward lowPt value( " << col
441 [ # # ][ # # ]: 2016 : << " ) = " << p_height[col] << " ilo = " << ilo << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
442 : : }
443 : : // MSQ_PRINT(3)("Contract all directions toward lowPt value( %d ) =
444 : : // %e ilo = %d\n", col, p_height[col], ilo);
445 : : }
446 : : }
447 : : }
448 : : }
449 : : } // ytri > h(ilo)
450 : : } // after evaluation
451 : 8686 : ilo = 1; // conditional operator or inline if
452 [ + - ][ + - ]: 8686 : ihi = p_height[0] > p_height[1] ? ( inhi = 1, 0 ) : ( inhi = 0, 1 );
[ + + ]
453 [ + + ]: 34744 : for( int col = 0; col < numCol; col++ )
454 : : {
455 [ + - ][ + - ]: 26058 : if( p_height[col] <= p_height[ilo] )
[ + + ]
456 : : {
457 : 21035 : ilo = col; // ilo := argmin height
458 : : }
459 [ + - ][ + - ]: 26058 : if( p_height[col] > p_height[ihi] )
[ + + ]
460 : : {
461 : 1169 : inhi = ihi;
462 : 1169 : ihi = col;
463 : : }
464 : : else // p_height[ihi] >= p_height[col]
465 [ + + ][ + - ]: 24889 : if( col != ihi && p_height[col] > p_height[inhi] )
[ + - ][ + + ]
[ + + ]
466 : 1109 : inhi = col;
467 : : }
468 : : // rtol=2.0*fabs( p_height[ihi]-p_height[ilo] )/
469 : : // ( fabs(p_height[ihi])+fabs(p_height[ilo])+threshold );
470 : 8686 : afterEvaluation = true;
471 : : } // while not converged
472 : :
473 : : // Always set to current best mesh.
474 : : {
475 [ + + ]: 850 : if( ilo != 0 )
476 : : {
477 [ + - ]: 778 : double yTemp = p_height[0];
478 [ + - ][ + - ]: 778 : p_height[0] = p_height[ilo]; // height dimension numCol
479 [ + - ]: 778 : p_height[ilo] = yTemp;
480 [ + + ]: 1556 : for( int row = 1; row < numRow; row++ )
481 : : {
482 [ + - ]: 778 : yTemp = simplex[row];
483 [ + - ][ + - ]: 778 : simplex[row] = simplex[row + ilo * numRow];
484 [ + - ]: 778 : simplex[row + ilo * numRow] = yTemp;
485 : : }
486 : : }
487 : : }
488 [ + - ][ - + ]: 850 : if( pd.num_free_vertices() > 1 )
489 : : {
490 : : MSQ_SETERR( err )
491 [ # # ][ # # ]: 0 : ( "Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED );
492 : : }
493 : :
494 [ + - ][ + - ]: 850 : Vector3D newPoint( &simplex[0] );
495 : 850 : size_t vertexIndex = 0; // fix c.f. freeVertexIndex
496 [ + - ]: 850 : pd.set_vertex_coordinates( newPoint, vertexIndex, err );
497 [ + - ]: 850 : pd.snap_vertex_to_domain( vertexIndex, err );
498 [ + - ][ + - ]: 850 : if( term_crit->terminate() )
499 : : {
500 [ - + ][ # # ]: 850 : if( mNonGradDebug >= 1 ) { std::cout << "Optimization Termination OptStatus: Max Iter Exceeded" << std::endl; }
[ # # ]
501 : : // MSQ_PRINT(1)("Optimization Termination OptStatus: Max Iter Exceeded\n");
502 : 850 : }
503 : 850 : }
504 : :
505 : 50 : void NonGradient::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ )
506 : : {
507 [ - + ]: 50 : if( mNonGradDebug >= 2 ) { std::cout << "- Executing NonGradient::iteration_complete()" << std::endl; }
508 : : // MSQ_PRINT(2)("\n - Executing NonGradient::iteration_complete() \n");
509 : 50 : }
510 : :
511 : 2 : void NonGradient::cleanup()
512 : : {
513 [ - + ]: 2 : if( mNonGradDebug >= 2 ) { std::cout << " - Executing NonGradient::iteration_end()" << std::endl; }
514 : : // MSQ_PRINT(2)("\n - Executing NonGradient::iteration_end() \n");
515 : 2 : }
516 : :
517 [ + - ][ + - ]: 120 : } // namespace MBMesquite
|
