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 : : // -*- Mode : c++; tab-width: 3; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 3
28 : : // -*-
29 : :
30 : : /*! \File NonSmoothDescent.cpp \brief
31 : :
32 : : Implements the NonSmoothDescent class member functions.
33 : :
34 : : \author Lori Freitag
35 : : \date 2002-07-20 */
36 : :
37 : : #include <stdlib.h>
38 : : #include <stdio.h>
39 : : #include "NonSmoothDescent.hpp"
40 : : #include "MsqTimer.hpp"
41 : :
42 : : #undef NUMERICAL_GRADIENT
43 : :
44 : : namespace MBMesquite
45 : : {
46 : :
47 : : const double EPSILON = 1e-16;
48 : : const int MSQ_MAX_OPT_ITER = 20;
49 : :
50 : : enum Rotate
51 : : {
52 : : COUNTERCLOCKWISE = 1,
53 : : CLOCKWISE = 0
54 : : };
55 : :
56 [ + - ][ + - ]: 2 : NonSmoothDescent::NonSmoothDescent( QualityMetric* qm ) : currentQM( qm )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
57 : : {
58 : :
59 : : MSQ_DBGOUT( 1 ) << "- Executed NonSmoothDescent::NonSmoothDescent()\n";
60 : 1 : }
61 : :
62 : 0 : std::string NonSmoothDescent::get_name() const
63 : : {
64 [ # # ]: 0 : return "NonSmoothDescent";
65 : : }
66 : :
67 : 1 : PatchSet* NonSmoothDescent::get_patch_set()
68 : : {
69 : 1 : return &patchSet;
70 : : }
71 : :
72 : 1 : void NonSmoothDescent::initialize( PatchData& /*pd*/, MsqError& /*err*/ )
73 : : {
74 : 1 : minStepSize = 1e-6;
75 : : MSQ_DBGOUT( 1 ) << "- Executed NonSmoothDescent::initialize()\n";
76 : 1 : }
77 : :
78 : 1322 : void NonSmoothDescent::initialize_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
79 : 2644 : void NonSmoothDescent::optimize_vertex_positions( PatchData& pd, MsqError& err )
80 : : {
81 : : MSQ_FUNCTION_TIMER( "NonSmoothDescent" );
82 : :
83 : : // cout << "- Executing NonSmoothDescent::optimize_node_positions()\n";
84 : : /* perform the min max smoothing algorithm */
85 : : MSQ_PRINT( 2 )( "\nInitializing the patch iteration\n" );
86 : :
87 : : MSQ_PRINT( 3 )( "Number of Vertices: %d\n", (int)pd.num_nodes() );
88 : : MSQ_PRINT( 3 )( "Number of Elements: %d\n", (int)pd.num_elements() );
89 : : // Michael: Note: is this a reliable way to get the dimension?
90 : : // NOTE: Mesquite only does 3-dimensional (including surface) meshes.
91 : : // mDimension = pd.get_mesh()->get_geometric_dimension(err); MSQ_ERRRTN(err);
92 : : // MSQ_PRINT(3)("Spatial Dimension: %d\n",mDimension);
93 : : MSQ_PRINT( 3 )( "Spatial Dimension: 3\n" );
94 : :
95 : : MSQ_PRINT( 3 )( "Num Free = %d\n", (int)pd.num_free_vertices() );
96 : :
97 [ + - ][ + - ]: 2644 : MsqFreeVertexIndexIterator free_iter( pd, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
98 [ + - ]: 1322 : free_iter.reset();
99 [ + - ]: 1322 : free_iter.next();
100 [ + - ]: 1322 : freeVertexIndex = free_iter.value();
101 : : MSQ_PRINT( 3 )( "Free Vertex Index = %lu\n", (unsigned long)freeVertexIndex );
102 : :
103 : : // TODO - need to switch to validity via metric evaluations should
104 : : // be associated with the compute_function somehow
105 : : /* check for an invalid mesh; if it's invalid return and ask the user
106 : : to use untangle */
107 [ + - ][ - + ]: 1322 : if( !this->validity_check( pd, err ) )
108 : : {
109 : : MSQ_PRINT( 1 )( "ERROR: Invalid mesh\n" );
110 : : MSQ_SETERR( err )
111 : : ( "Invalid Mesh: Use untangle to create a valid "
112 : : "triangulation",
113 [ # # ][ # # ]: 0 : MsqError::INVALID_MESH );
114 : 0 : return;
115 : : }
116 : :
117 : : /* initialize the optimization data up to numFunctionValues */
118 [ + - ][ + - ]: 1322 : this->init_opt( pd, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
119 [ + - ][ + - ]: 1322 : this->init_max_step_length( pd, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
120 : : MSQ_PRINT( 3 )( "Done initializing optimization\n" );
121 : :
122 : : /* compute the initial function values */
123 : : // TODO this should return a bool with the validity
124 [ + - ][ + - ]: 1322 : this->compute_function( &pd, originalFunction, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
125 : :
126 : : // find the initial active set
127 [ + - ]: 1322 : this->find_active_set( originalFunction, mActive );
128 : :
129 [ + - ][ + - ]: 1322 : this->minmax_opt( pd, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
130 : : }
131 : :
132 : 1 : void NonSmoothDescent::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
133 : :
134 : 2 : void NonSmoothDescent::cleanup()
135 : : {
136 : : MSQ_DBGOUT( 1 ) << "- Executing NonSmoothDescent::cleanup()\n";
137 : : MSQ_DBGOUT( 1 ) << "- Done with NonSmoothDescent::cleanup()\n";
138 : 1 : }
139 : :
140 : 5789 : void NonSmoothDescent::find_plane_points( Direction dir1, Direction dir2, const std::vector< Vector3D >& vec,
141 : : Vector3D& pt1, Vector3D& pt2, Vector3D& /*pt3*/, Status& status, MsqError& )
142 : : {
143 : : int i;
144 : : int num_min, num_max;
145 : 5789 : Rotate rotate = CLOCKWISE;
146 : 5789 : int num_rotated = 0;
147 : : double pt_1, pt_2;
148 : : double min, inv_slope;
149 : 5789 : double min_inv_slope = 0.;
150 : : double max;
151 : 5789 : double max_inv_slope = 0;
152 : 5789 : double inv_origin_slope = 0;
153 : 5789 : const int num_vec = vec.size();
154 : 5789 : const int auto_ind_size = 50;
155 : : int auto_ind[auto_ind_size];
156 [ + - ]: 5789 : std::vector< int > heap_ind;
157 : : int* ind;
158 [ + - ]: 5789 : if( num_vec <= auto_ind_size )
159 : 5789 : ind = auto_ind;
160 : : else
161 : : {
162 [ # # ]: 0 : heap_ind.resize( num_vec );
163 [ # # ]: 0 : ind = &heap_ind[0];
164 : : }
165 : :
166 : 5789 : status = MSQ_CHECK_BOTTOM_UP;
167 : : /* find the minimum points in dir1 starting at -1 */
168 : 5789 : num_min = 0;
169 : 5789 : ind[0] = -1;
170 : 5789 : ind[1] = -1;
171 : 5789 : ind[2] = -1;
172 : 5789 : min = 1.0;
173 [ + + ]: 24696 : for( i = 0; i < num_vec; i++ )
174 : : {
175 [ + - ][ + - ]: 18907 : if( vec[i][dir1] < min )
[ + + ]
176 : : {
177 [ + - ][ + - ]: 11065 : min = vec[i][dir1];
178 : 11065 : ind[0] = i;
179 : 11065 : num_min = 1;
180 : : }
181 [ + - ][ + - ]: 7842 : else if( fabs( vec[i][dir1] - min ) < EPSILON )
[ - + ]
182 : : {
183 : 0 : ind[num_min++] = i;
184 : : }
185 : : }
186 [ + + ]: 5789 : if( min >= 0 ) status = MSQ_NO_EQUIL;
187 : :
188 [ + + ]: 5789 : if( status != MSQ_NO_EQUIL )
189 : : {
190 [ + - - ]: 5694 : switch( num_min )
191 : : {
192 : : case 1: /* rotate to find the next point */
193 [ + - ][ + - ]: 5694 : pt1 = vec[ind[0]];
194 [ + - ]: 5694 : pt_1 = pt1[dir1];
195 [ + - ]: 5694 : pt_2 = pt1[dir2];
196 [ + - ][ + + ]: 5694 : if( pt1[dir2] <= 0 )
197 : : {
198 : 2887 : rotate = COUNTERCLOCKWISE;
199 : 2887 : max_inv_slope = -HUGE_VAL;
200 : : }
201 [ + - ][ + + ]: 5694 : if( pt1[dir2] > 0 )
202 : : {
203 : 2807 : rotate = CLOCKWISE;
204 : 2807 : min_inv_slope = HUGE_VAL;
205 : : }
206 [ + + - ]: 5694 : switch( rotate )
207 : : {
208 : : case COUNTERCLOCKWISE:
209 [ + + ]: 12319 : for( i = 0; i < num_vec; i++ )
210 : : {
211 [ + + ]: 9432 : if( i != ind[0] )
212 : : {
213 [ + - ][ + - ]: 6545 : inv_slope = ( vec[i][dir2] - pt_2 ) / ( vec[i][dir1] - pt_1 );
[ + - ][ + - ]
214 [ + + ][ + - ]: 6545 : if( ( inv_slope > max_inv_slope ) && ( fabs( inv_slope - max_inv_slope ) > EPSILON ) )
215 : : {
216 : 4601 : ind[1] = i;
217 : 4601 : max_inv_slope = inv_slope;
218 : 4601 : num_rotated = 1;
219 : : }
220 [ - + ]: 1944 : else if( fabs( inv_slope - max_inv_slope ) < EPSILON )
221 : : {
222 : 0 : ind[2] = i;
223 : 6545 : num_rotated++;
224 : : }
225 : : }
226 : : }
227 : 2887 : break;
228 : : case CLOCKWISE:
229 [ + + ]: 11992 : for( i = 0; i < num_vec; i++ )
230 : : {
231 [ + + ]: 9185 : if( i != ind[0] )
232 : : {
233 [ + - ][ + - ]: 6378 : inv_slope = ( vec[i][dir2] - pt_2 ) / ( vec[i][dir1] - pt_1 );
[ + - ][ + - ]
234 [ + + ][ + - ]: 6378 : if( ( inv_slope < min_inv_slope ) && ( fabs( inv_slope - max_inv_slope ) > EPSILON ) )
235 : : {
236 : 4469 : ind[1] = i;
237 : 4469 : min_inv_slope = inv_slope;
238 : 4469 : num_rotated = 1;
239 : : }
240 [ - + ]: 1909 : else if( fabs( inv_slope - min_inv_slope ) < EPSILON )
241 : : {
242 : 0 : ind[2] = i;
243 : 6378 : num_rotated++;
244 : : }
245 : : }
246 : : }
247 : : }
248 [ - + - ]: 5694 : switch( num_rotated )
249 : : {
250 : : case 0:
251 : : MSQ_PRINT( 3 )( "No points in the rotation ... odd\n" );
252 : 0 : status = MSQ_HULL_TEST_ERROR;
253 : 0 : break;
254 : : case 1:
255 : : MSQ_PRINT( 3 )( "Found a line in the convex hull\n" );
256 [ + - ][ + - ]: 5694 : pt2 = vec[ind[1]];
257 : 5694 : status = MSQ_TWO_PT_PLANE;
258 : 5694 : break;
259 : : default:
260 : : MSQ_PRINT( 3 )( "Found 2 or more points in the rotation\n" );
261 [ # # ]: 0 : if( fabs( pt_1 ) > EPSILON ) inv_origin_slope = pt_2 / pt_1;
262 [ # # # ]: 0 : switch( rotate )
263 : : {
264 : : case COUNTERCLOCKWISE:
265 [ # # ]: 0 : if( inv_origin_slope >= max_inv_slope )
266 : 0 : status = MSQ_NO_EQUIL;
267 : : else
268 : 0 : status = MSQ_CHECK_TOP_DOWN;
269 : 0 : break;
270 : : case CLOCKWISE:
271 [ # # ]: 0 : if( inv_origin_slope <= min_inv_slope )
272 : 0 : status = MSQ_NO_EQUIL;
273 : : else
274 : 0 : status = MSQ_CHECK_TOP_DOWN;
275 : : }
276 : : }
277 : 5694 : break;
278 : : case 2: /* use these two points to define the plane */
279 : : MSQ_PRINT( 3 )( "Found two minimum points to define the plane\n" );
280 [ # # ][ # # ]: 0 : pt1 = vec[ind[0]];
281 [ # # ][ # # ]: 0 : pt2 = vec[ind[1]];
282 : 0 : status = MSQ_TWO_PT_PLANE;
283 : 0 : break;
284 : : default: /* check to see if all > 0 */
285 : : MSQ_PRINT( 3 )( "Found 3 or more points in min plane %f\n", min );
286 [ # # ][ # # ]: 0 : if( vec[ind[0]][dir1] >= 0 )
[ # # ]
287 : 0 : status = MSQ_NO_EQUIL;
288 : : else
289 : 5694 : status = MSQ_CHECK_TOP_DOWN;
290 : : }
291 : : }
292 : :
293 : : /***************************/
294 : : /* failed to find any information, checking top/down this coord*/
295 : : /***************************/
296 : :
297 [ - + ]: 5789 : if( status == MSQ_CHECK_TOP_DOWN )
298 : : {
299 : : /* find the maximum points in dir1 starting at 1 */
300 : 0 : num_max = 0;
301 : 0 : ind[0] = -1;
302 : 0 : ind[1] = -1;
303 : 0 : ind[2] = -1;
304 : 0 : max = -1.0;
305 [ # # ]: 0 : for( i = 0; i < num_vec; i++ )
306 : : {
307 [ # # ][ # # ]: 0 : if( vec[i][dir1] > max )
[ # # ]
308 : : {
309 [ # # ][ # # ]: 0 : max = vec[i][dir1];
310 : 0 : ind[0] = i;
311 : 0 : num_max = 1;
312 : : }
313 [ # # ][ # # ]: 0 : else if( fabs( vec[i][dir1] - max ) < EPSILON )
[ # # ]
314 : : {
315 : 0 : ind[num_max++] = i;
316 : : }
317 : : }
318 [ # # ]: 0 : if( max <= 0 ) status = MSQ_NO_EQUIL;
319 : :
320 [ # # ]: 0 : if( status != MSQ_NO_EQUIL )
321 : : {
322 [ # # # ]: 0 : switch( num_max )
323 : : {
324 : : case 1: /* rotate to find the next point */
325 [ # # ][ # # ]: 0 : pt1 = vec[ind[0]];
326 [ # # ]: 0 : pt_1 = pt1[dir1];
327 [ # # ]: 0 : pt_2 = pt1[dir2];
328 [ # # ][ # # ]: 0 : if( pt1[dir2] < 0 )
329 : : {
330 : 0 : rotate = CLOCKWISE;
331 : 0 : min_inv_slope = HUGE_VAL;
332 : : }
333 [ # # ][ # # ]: 0 : if( pt1[dir2] >= 0 )
334 : : {
335 : 0 : rotate = COUNTERCLOCKWISE;
336 : 0 : max_inv_slope = -HUGE_VAL;
337 : : }
338 [ # # # ]: 0 : switch( rotate )
339 : : {
340 : : case COUNTERCLOCKWISE:
341 [ # # ]: 0 : for( i = 0; i < num_vec; i++ )
342 : : {
343 [ # # ]: 0 : if( i != ind[0] )
344 : : {
345 [ # # ][ # # ]: 0 : inv_slope = ( vec[i][dir2] - pt_2 ) / ( vec[i][dir1] - pt_1 );
[ # # ][ # # ]
346 [ # # ]: 0 : if( inv_slope > max_inv_slope )
347 : : {
348 : 0 : ind[1] = i;
349 : 0 : max_inv_slope = inv_slope;
350 : 0 : num_rotated = 1;
351 : : }
352 [ # # ]: 0 : else if( fabs( inv_slope - max_inv_slope ) < EPSILON )
353 : : {
354 : 0 : ind[2] = i;
355 : 0 : num_rotated++;
356 : : }
357 : : }
358 : : }
359 : 0 : break;
360 : : case CLOCKWISE:
361 [ # # ]: 0 : for( i = 0; i < num_vec; i++ )
362 : : {
363 [ # # ]: 0 : if( i != ind[0] )
364 : : {
365 [ # # ][ # # ]: 0 : inv_slope = ( vec[i][dir2] - pt_2 ) / ( vec[i][dir1] - pt_1 );
[ # # ][ # # ]
366 [ # # ]: 0 : if( inv_slope < min_inv_slope )
367 : : {
368 : 0 : ind[1] = i;
369 : 0 : min_inv_slope = inv_slope;
370 : 0 : num_rotated = 1;
371 : : }
372 [ # # ]: 0 : else if( fabs( inv_slope - min_inv_slope ) < EPSILON )
373 : : {
374 : 0 : ind[2] = i;
375 : 0 : num_rotated++;
376 : : }
377 : : }
378 : : }
379 : : }
380 [ # # # ]: 0 : switch( num_rotated )
381 : : {
382 : : case 0:
383 : : MSQ_PRINT( 3 )( "No points in the rotation ... odd\n" );
384 : 0 : status = MSQ_HULL_TEST_ERROR;
385 : 0 : break;
386 : : case 1:
387 : : MSQ_PRINT( 3 )( "Found a line in the convex hull\n" );
388 [ # # ][ # # ]: 0 : pt2 = vec[ind[1]];
389 : 0 : status = MSQ_TWO_PT_PLANE;
390 : 0 : break;
391 : : default:
392 : : MSQ_PRINT( 3 )( "Found 2 or more points in the rotation\n" );
393 : : /* check to see if rotation got past origin */
394 : 0 : inv_origin_slope = pt_2 / pt_1;
395 [ # # # ]: 0 : switch( rotate )
396 : : {
397 : : case COUNTERCLOCKWISE:
398 [ # # ]: 0 : if( inv_origin_slope >= max_inv_slope )
399 : 0 : status = MSQ_NO_EQUIL;
400 [ # # ]: 0 : else if( dir1 == 2 )
401 : 0 : status = MSQ_CHECK_Y_COORD_DIRECTION;
402 [ # # ]: 0 : else if( dir1 == 1 )
403 : 0 : status = MSQ_CHECK_X_COORD_DIRECTION;
404 : : else
405 : 0 : status = MSQ_EQUIL;
406 : 0 : break;
407 : : case CLOCKWISE:
408 [ # # ]: 0 : if( inv_origin_slope <= min_inv_slope )
409 : 0 : status = MSQ_NO_EQUIL;
410 [ # # ]: 0 : else if( dir1 == 2 )
411 : 0 : status = MSQ_CHECK_Y_COORD_DIRECTION;
412 [ # # ]: 0 : else if( dir1 == 1 )
413 : 0 : status = MSQ_CHECK_X_COORD_DIRECTION;
414 : : else
415 : 0 : status = MSQ_EQUIL;
416 : : }
417 : : }
418 : 0 : break;
419 : : case 2: /* use these two points to define the plane */
420 [ # # ][ # # ]: 0 : pt1 = vec[ind[0]];
421 [ # # ][ # # ]: 0 : pt2 = vec[ind[1]];
422 : 0 : status = MSQ_TWO_PT_PLANE;
423 : 0 : break;
424 : : default: /* check to see if all > 0 */
425 : : MSQ_PRINT( 3 )( "Found 3 in max plane %f\n", max );
426 [ # # ][ # # ]: 0 : if( vec[ind[0]][dir1] <= 0 )
[ # # ]
427 : 0 : status = MSQ_NO_EQUIL;
428 [ # # ]: 0 : else if( dir1 == 2 )
429 : 0 : status = MSQ_CHECK_Y_COORD_DIRECTION;
430 [ # # ]: 0 : else if( dir1 == 1 )
431 : 0 : status = MSQ_CHECK_X_COORD_DIRECTION;
432 : : else
433 : 0 : status = MSQ_EQUIL;
434 : : }
435 : : }
436 : 5789 : }
437 : 5789 : }
438 : :
439 : 13463 : void NonSmoothDescent::search_direction( PatchData& /*pd*/, Vector3D& mSearch, MsqError& err )
440 : : {
441 : : bool viable;
442 : : double a, b, c, denom;
443 [ + - ]: 13463 : std::vector< Vector3D > dir;
444 : : double R0, R1;
445 [ + - ][ + - ]: 26926 : SymmetricMatrix P;
446 : : double x[2];
447 : : double search_mag;
448 : :
449 : 13463 : const int num_active = mActive.active_ind.size();
450 : : ;
451 : :
452 : : // TODO This might be o.k. actually - i don't see any dependence
453 : : // on the element geometry here... try it and see if it works.
454 : : // if not, try taking all of the gradients in the active set
455 : : // and let the search direction be the average of those.
456 : : // MSQ_FUNCTION_TIMER( "Search Direction" );
457 : :
458 : : MSQ_PRINT( 2 )( "\nIn Search Direction\n" );
459 [ + - ]: 13463 : this->print_active_set( mActive, mFunction );
460 : :
461 [ - + ]: 13463 : if( num_active == 0 )
462 : : {
463 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "No active values in search", MsqError::INVALID_STATE );
464 : 0 : return;
465 : : }
466 : :
467 [ + + + ]: 13463 : switch( num_active )
468 : : {
469 : : case 1:
470 [ + - ][ + - ]: 5987 : mSearch = mGradient[mActive.active_ind[0]];
[ + - ]
471 [ + - ]: 5987 : mSteepest = mActive.active_ind[0];
472 : 5987 : break;
473 : : case 2:
474 : : /* if there are two active points, move in the direction of the
475 : : intersection of the planes. This is the steepest descent
476 : : direction found by analytically solving the QP */
477 : :
478 : : /* set up the active gradient directions */
479 [ + - ][ + - ]: 4342 : this->get_active_directions( mGradient, dir, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
480 : :
481 : : /* form the grammian */
482 [ + - ][ + - ]: 4342 : this->form_grammian( dir, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
483 [ + - ][ + - ]: 4342 : this->form_PD_grammian( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
484 : :
485 [ + - ][ + - ]: 4342 : denom = ( mG( 0, 0 ) + mG( 1, 1 ) - 2 * mG( 0, 1 ) );
[ + - ]
486 : 4342 : viable = true;
487 [ + - ]: 4342 : if( fabs( denom ) > EPSILON )
488 : : {
489 : : /* gradients are LI, move along their intersection */
490 [ + - ][ + - ]: 4342 : b = ( mG( 0, 0 ) - mG( 0, 1 ) ) / denom;
491 : 4342 : a = 1 - b;
492 [ + + ][ + + ]: 4342 : if( ( b < 0 ) || ( b > 1 ) ) viable = false; /* 0 < b < 1 */
493 [ + + ][ + - ]: 4342 : if( viable ) { mSearch = a * dir[0] + b * dir[1]; }
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
494 : : else
495 : : {
496 : : /* the gradients are dependent, move along one face */
497 [ + - ][ + - ]: 4342 : mSearch = dir[0];
498 : : }
499 : : }
500 : : else
501 : : {
502 : : /* the gradients are dependent, move along one face */
503 [ # # ][ # # ]: 0 : mSearch = dir[0];
504 : : }
505 [ + - ]: 4342 : mSteepest = mActive.active_ind[0];
506 : :
507 : 4342 : break;
508 : : default:
509 : : /* as in case 2: solve the QP problem to find the steepest
510 : : descent direction. This can be done analytically - as
511 : : is done in Gill, Murray and Wright
512 : : for 3 active points in 3 directions - test PD of G
513 : : otherwise we know it's SP SD so search edges and faces */
514 : :
515 : : /* get the active gradient directions */
516 [ + - ][ + - ]: 3134 : this->get_active_directions( mGradient, dir, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
517 : :
518 : : /* form the entries of the grammian matrix */
519 [ + - ][ + - ]: 3134 : this->form_grammian( dir, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
520 [ + - ][ + - ]: 3134 : this->form_PD_grammian( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
521 : :
522 [ + + ]: 3134 : if( num_active == 3 )
523 : : {
524 [ + - ][ + - ]: 2991 : if( mG.condition3x3() < 1e14 )
525 : : { // if not singular
526 : : /* form the entries of P=Z^T G Z where Z = [-1...-1; I ] */
527 [ + - ]: 2991 : this->form_reduced_matrix( P );
528 : : /* form the RHS and solve the system for the coeffs */
529 [ + - ][ + - ]: 2991 : R0 = mG( 0, 0 ) - mG( 1, 0 );
530 [ + - ][ + - ]: 2991 : R1 = mG( 0, 0 ) - mG( 2, 0 );
531 [ + - ][ + - ]: 2991 : bool ok = this->solve2x2( P( 0, 0 ), P( 0, 1 ), P( 1, 0 ), P( 1, 1 ), R0, R1, x, err );MSQ_ERRRTN( err );
[ + - ][ + - ]
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ - + ]
532 [ + - ]: 2991 : if( ok )
533 : : {
534 : 2991 : a = 1 - x[0] - x[1];
535 : 2991 : b = x[0];
536 : 2991 : c = x[1];
537 [ + - ][ + - ]: 2991 : mSearch = a * dir[0] + b * dir[1] + c * dir[2];
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
538 [ + - ]: 2991 : mSteepest = mActive.active_ind[0];
539 : : }
540 : : else
541 : : {
542 [ # # ][ # # ]: 2991 : this->search_edges_faces( &dir[0], mSearch, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
543 : : }
544 : : }
545 : : else
546 : : {
547 [ # # ][ # # ]: 2991 : this->search_edges_faces( &dir[0], mSearch, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
548 : : }
549 : : }
550 : : else
551 : : {
552 [ + - ][ + - ]: 143 : this->search_edges_faces( &dir[0], mSearch, err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
553 : : }
554 : : }
555 : :
556 : : /* if the search direction is essentially zero, equilibrium pt */
557 [ + - ]: 13463 : search_mag = mSearch % mSearch;
558 : : MSQ_PRINT( 3 )( " Search Magnitude %g \n", search_mag );
559 : :
560 [ + + ]: 13463 : if( fabs( search_mag ) < 1E-13 )
561 : 92 : optStatus = MSQ_ZERO_SEARCH;
562 : : else
563 [ + - ][ + - ]: 13463 : mSearch /= std::sqrt( search_mag );
564 : : MSQ_PRINT( 3 )
565 : 13463 : ( " Search Direction %g %g Steepest %lu\n", mSearch[0], mSearch[1], (unsigned long)mSteepest );
566 : : }
567 : :
568 : 1322 : void NonSmoothDescent::minmax_opt( PatchData& pd, MsqError& err )
569 : : {
570 : 1322 : bool equilibriumPt = false;
571 [ + - ]: 1322 : Vector3D mSearch( 0, 0, 0 );
572 : :
573 : : // int valid;
574 : : MSQ_FUNCTION_TIMER( "Minmax Opt" );
575 : : MSQ_PRINT( 2 )( "In minmax_opt\n" );
576 : :
577 [ + - ]: 1322 : mFunction = originalFunction;
578 : 1322 : originalValue = mActive.true_active_value;
579 : :
580 : 1322 : int iterCount = 0;
581 : 1322 : int optIterCount = 0;
582 : :
583 : : MSQ_PRINT( 3 )( "Done copying original function to function\n" );
584 : :
585 [ + - ]: 1322 : this->find_active_set( mFunction, mActive );
586 : 1322 : prevActiveValues.clear();
587 [ + - ]: 1322 : prevActiveValues.push_back( mActive.true_active_value );
588 : :
589 : : /* check for equilibrium point */
590 : : /* compute the gradient */
591 [ + - ][ + - ]: 2644 : this->compute_gradient( &pd, mGradient, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
592 : :
593 [ + + ]: 1322 : if( mActive.active_ind.size() >= 2 )
594 : : {
595 : : MSQ_PRINT( 3 )( "Testing for an equilibrium point \n" );
596 [ + - ][ + - ]: 77 : equilibriumPt = this->check_equilibrium( optStatus, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
597 : :
598 : : if( MSQ_DBG( 2 ) && equilibriumPt ) MSQ_PRINT( 2 )( "Optimization Exiting: An equilibrium point \n" );
599 : : }
600 : :
601 : : /* terminate if we have found an equilibrium point or if the step is
602 : : too small to be worthwhile continuing */
603 [ + + ][ + + ]: 14785 : while( ( optStatus != MSQ_EQUILIBRIUM ) && ( optStatus != MSQ_STEP_TOO_SMALL ) &&
[ + + ]
604 [ + + ][ + + ]: 13668 : ( optStatus != MSQ_IMP_TOO_SMALL ) && ( optStatus != MSQ_FLAT_NO_IMP ) && ( optStatus != MSQ_ZERO_SEARCH ) &&
[ + + ]
605 : 13489 : ( optStatus != MSQ_MAX_ITER_EXCEEDED ) )
606 : : {
607 : :
608 : : /* increase the iteration count by one */
609 : : /* smooth_param->iter_count += 1; */
610 : 13463 : iterCount += 1;
611 : 13463 : optIterCount += 1;
612 [ + + ]: 13463 : if( iterCount > MSQ_MAX_OPT_ITER ) optStatus = MSQ_MAX_ITER_EXCEEDED;
613 : :
614 : : MSQ_PRINT( 3 )( "\nITERATION %d \n", iterCount );
615 : :
616 : : /* compute the gradient */
617 [ + - ][ + - ]: 13463 : this->compute_gradient( &pd, mGradient, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
618 : :
619 : : MSQ_PRINT( 3 )( "Computing the search direction \n" );
620 [ + - ][ + - ]: 13463 : this->search_direction( pd, mSearch, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
621 : :
622 : : /* if there are viable directions to search */
623 [ + + ][ + + ]: 13463 : if( ( optStatus != MSQ_ZERO_SEARCH ) && ( optStatus != MSQ_MAX_ITER_EXCEEDED ) )
624 : : {
625 : :
626 : : MSQ_PRINT( 3 )( "Computing the projections of the gradients \n" );
627 [ + - ][ + - ]: 13345 : this->get_gradient_projections( mSearch, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
628 : :
629 : : MSQ_PRINT( 3 )( "Computing the initial step size \n" );
630 [ + - ][ + - ]: 13345 : this->compute_alpha( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
631 : :
632 : : MSQ_PRINT( 3 )( "Testing whether to accept this step \n" );
633 [ + - ][ + - ]: 13345 : this->step_acceptance( pd, iterCount, mSearch, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
634 : : // MSQ_PRINT(3)("The new free vertex position is %f %f %f\n",
635 : : // mCoords[freeVertexIndex][0],mCoords[freeVertexIndex][1],mCoords[freeVertexIndex][2]);
636 : :
637 : : if( MSQ_DBG( 3 ) )
638 : : {
639 : : /* Print the active set */
640 : : this->print_active_set( mActive, mFunction );
641 : : }
642 : :
643 : : /* check for equilibrium point */
644 [ + + ]: 13345 : if( mActive.active_ind.size() >= 2 )
645 : : {
646 : : MSQ_PRINT( 3 )( "Testing for an equilibrium point \n" );
647 [ + - ][ + - ]: 8576 : equilibriumPt = this->check_equilibrium( optStatus, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
648 : :
649 : : if( MSQ_DBG( 2 ) && equilibriumPt ) MSQ_PRINT( 2 )( "Optimization Exiting: An equilibrium point \n" );
650 : : }
651 : :
652 : : /* record the values */
653 [ + - ]: 13345 : prevActiveValues.push_back( mActive.true_active_value );
654 : : }
655 : : else
656 : : {
657 : : /* decrease the iteration count by one */
658 : : /* smooth_param->iter_count -= 1; */
659 : 118 : iterCount -= 1;
660 : : if( MSQ_DBG( 2 ) )
661 : : {
662 : : MSQ_PRINT( 2 )
663 : : ( "Optimization Exiting: No viable directions; equilibrium point \n" );
664 : : /* Print the old active set */
665 : : this->print_active_set( mActive, mFunction );
666 : : }
667 : : }
668 : : }
669 : :
670 : : MSQ_PRINT( 2 )( "Checking the validity of the mesh\n" );
671 [ + - ][ + - ]: 1322 : if( !this->validity_check( pd, err ) ) MSQ_PRINT( 2 )( "The final mesh is not valid\n" );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
672 : :
673 : : MSQ_PRINT( 2 )( "Number of optimization iterations %d\n", iterCount );
674 : :
675 [ - + + + : 1322 : switch( optStatus )
+ + + ]
676 : : {
677 : : default:
678 : : MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Invalid early termination\n" );
679 : 0 : break;
680 : : case MSQ_EQUILIBRIUM:
681 : : MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Equilibrium\n" );
682 : 499 : break;
683 : : case MSQ_STEP_TOO_SMALL:
684 : : MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Step Too Small\n" );
685 : 330 : break;
686 : : case MSQ_IMP_TOO_SMALL:
687 : : MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Improvement Too Small\n" );
688 : 288 : break;
689 : : case MSQ_FLAT_NO_IMP:
690 : : MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Flat No Improvement\n" );
691 : 87 : break;
692 : : case MSQ_ZERO_SEARCH:
693 : : MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Zero Search\n" );
694 : 92 : break;
695 : : case MSQ_MAX_ITER_EXCEEDED:
696 : : MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Max Iter Exceeded\n" );
697 : 1322 : break;
698 : : }
699 : : }
700 : :
701 : 13345 : void NonSmoothDescent::step_acceptance( PatchData& pd, int iterCount, const Vector3D& mSearch, MsqError& err )
702 : : {
703 : 13345 : const double minAcceptableImprovement = 1e-6;
704 : :
705 : : // int ierr;
706 : : // int num_values;
707 : : bool step_done;
708 : 13345 : bool valid = true, accept_alpha;
709 : : double estimated_improvement;
710 : 13345 : double current_improvement = HUGE_VAL;
711 : 13345 : double previous_improvement = HUGE_VAL;
712 : 13345 : double current_percent_diff = HUGE_VAL;
713 [ + - ]: 13345 : Vector3D original_point;
714 : :
715 : : // MSQ_FUNCTION_TIMER( "Step Acceptance" );
716 : : // num_values = qmHandles.size();
717 : :
718 : 13345 : optStatus = MSQ_NO_STATUS;
719 : :
720 [ - + ]: 13345 : if( mAlpha < minStepSize )
721 : : {
722 : 0 : optStatus = MSQ_IMP_TOO_SMALL;
723 : 0 : step_done = true;
724 : : MSQ_PRINT( 3 )( "Alpha starts too small, no improvement\n" );
725 : : }
726 : :
727 [ + - ]: 13345 : const MsqVertex* coords = pd.get_vertex_array( err );
728 : :
729 : : /* save the original function and active set */
730 [ + - ]: 13345 : original_point = coords[freeVertexIndex];
731 [ + - ]: 13345 : originalFunction = mFunction;
732 [ + - ]: 13345 : originalActive = mActive;
733 : :
734 : 13345 : step_done = false;
735 [ + + ][ + - ]: 41213 : for( int num_steps = 0; !step_done && num_steps < 100; ++num_steps )
736 : : {
737 : 27868 : accept_alpha = false;
738 : :
739 [ + + ][ + + ]: 61434 : while( !accept_alpha && mAlpha > minStepSize )
740 : : {
741 : :
742 : : /* make the step */
743 [ + - ][ + - ]: 33566 : pd.move_vertex( -mAlpha * mSearch, freeVertexIndex, err );
744 : : // pd.set_coords_array_element(coords[freeVertexIndex],0,err);
745 : :
746 : : MSQ_PRINT( 2 )( "search direction %f %f \n", mSearch[0], mSearch[1] );
747 : : MSQ_PRINT( 2 )
748 : : ( "new vertex position %f %f \n", coords[freeVertexIndex][0], coords[freeVertexIndex][1] );
749 : :
750 : : /* assume alpha is acceptable */
751 : 33566 : accept_alpha = true;
752 : :
753 : : /* never take a step that makes a valid mesh invalid or worsens the quality */
754 : : // TODO Validity check revision -- do the compute function up here
755 : : // and then the rest based on validity
756 [ + - ][ + - ]: 46911 : valid = validity_check( pd, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
757 [ + + ]: 33566 : if( valid )
758 : : {
759 [ + - ][ + - ]: 33492 : valid = improvement_check( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
760 : : }
761 [ + + ]: 33566 : if( !valid )
762 : : {
763 : 6028 : accept_alpha = false;
764 [ + - ][ + - ]: 6028 : pd.move_vertex( mAlpha * mSearch, freeVertexIndex, err );
765 : : // pd.set_coords_array_element(coords[freeVertexIndex],0,err);
766 : 6028 : mAlpha = mAlpha / 2;
767 : : MSQ_PRINT( 2 )( "Step not accepted, the new alpha %f\n", mAlpha );
768 : :
769 [ + + ]: 6028 : if( mAlpha < minStepSize )
770 : : {
771 : 330 : optStatus = MSQ_STEP_TOO_SMALL;
772 : 330 : step_done = true;
773 : : MSQ_PRINT( 2 )( "Step too small\n" );
774 : : /* get back the original point, mFunction, and active set */
775 [ + - ]: 330 : pd.set_vertex_coordinates( original_point, freeVertexIndex, err );
776 [ + - ]: 330 : mFunction = originalFunction;
777 [ + - ]: 330 : mActive = originalActive;
778 : : }
779 : : }
780 : : }
781 : :
782 [ + + ][ + - ]: 27868 : if( valid && ( mAlpha > minStepSize ) )
783 : : {
784 : : /* compute the new function and active set */
785 [ + - ][ + - ]: 27538 : this->compute_function( &pd, mFunction, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
786 [ + - ]: 27538 : this->find_active_set( mFunction, mActive );
787 : :
788 : : /* estimate the minimum improvement by taking this step */
789 [ + - ][ + - ]: 27538 : this->get_min_estimate( &estimated_improvement, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
790 : : MSQ_PRINT( 2 )
791 : : ( "The estimated improvement for this step: %f\n", estimated_improvement );
792 : :
793 : : /* calculate the actual increase */
794 [ + - ]: 27538 : current_improvement = mActive.true_active_value - prevActiveValues[iterCount - 1];
795 : :
796 : : MSQ_PRINT( 3 )( "Actual improvement %f\n", current_improvement );
797 : :
798 : : /* calculate the percent difference from estimated increase */
799 : 27538 : current_percent_diff = fabs( current_improvement - estimated_improvement ) / fabs( estimated_improvement );
800 : :
801 : : /* determine whether to accept a step */
802 [ + + ][ + + ]: 27538 : if( ( fabs( previous_improvement ) > fabs( current_improvement ) ) && ( previous_improvement < 0 ) )
803 : : {
804 : : /* accept the previous step - it was better */
805 : : MSQ_PRINT( 2 )( "Accepting the previous step\n" );
806 : :
807 : : /* subtract alpha in again (previous step) */
808 [ + - ][ + - ]: 6817 : pd.move_vertex( -mAlpha * mSearch, freeVertexIndex, err );
809 : : // pd.set_coords_array_element(coords[freeVertexIndex],0,err);
810 : :
811 : : /* does this make an invalid mesh valid? */
812 : : // TODO Validity check revisison
813 [ + - ][ + - ]: 6817 : valid = validity_check( pd, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
814 [ + - ][ + - ]: 6817 : if( valid ) valid = improvement_check( err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
815 : :
816 : : /* copy test function and active set */
817 [ + - ]: 6817 : mFunction = testFunction;
818 [ + - ]: 6817 : mActive = testActive;
819 : :
820 : 6817 : optStatus = MSQ_STEP_ACCEPTED;
821 : 6817 : step_done = true;
822 : :
823 : : /* check to see that we're still making good improvements */
824 [ + + ]: 7107 : if( fabs( previous_improvement ) < minAcceptableImprovement )
825 : : {
826 : 290 : optStatus = MSQ_IMP_TOO_SMALL;
827 : 290 : step_done = true;
828 : : MSQ_PRINT( 2 )( "Optimization Exiting: Improvement too small\n" );
829 : : }
830 : : }
831 [ + + ][ + + ]: 20721 : else if( ( ( fabs( current_improvement ) > fabs( estimated_improvement ) ) ||
832 [ + + ]: 10027 : ( current_percent_diff < .1 ) ) &&
833 : : ( current_improvement < 0 ) )
834 : : {
835 : : /* accept this step, exceeded estimated increase or was close */
836 : 6110 : optStatus = MSQ_STEP_ACCEPTED;
837 : 6110 : step_done = true;
838 : :
839 : : /* check to see that we're still making good improvements */
840 [ + + ]: 6114 : if( fabs( current_improvement ) < minAcceptableImprovement )
841 : : {
842 : : MSQ_PRINT( 2 )( "Optimization Exiting: Improvement too small\n" );
843 : 4 : optStatus = MSQ_IMP_TOO_SMALL;
844 : 4 : step_done = true;
845 : : }
846 : : }
847 [ + + ][ + - ]: 14611 : else if( ( current_improvement > 0 ) && ( previous_improvement > 0 ) &&
[ + + ]
848 [ + + ]: 396 : ( fabs( current_improvement ) < minAcceptableImprovement ) &&
849 : 396 : ( fabs( previous_improvement ) < minAcceptableImprovement ) )
850 : : {
851 : :
852 : : /* we are making no progress, quit */
853 : 88 : optStatus = MSQ_FLAT_NO_IMP;
854 : 88 : step_done = true;
855 : : MSQ_PRINT( 2 )( "Opimization Exiting: Flat no improvement\n" );
856 : :
857 : : /* get back the original point, function, and active set */
858 [ + - ][ + - ]: 88 : pd.set_vertex_coordinates( original_point, freeVertexIndex, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
859 [ + - ]: 88 : mFunction = originalFunction;
860 [ + - ]: 88 : mActive = originalActive;
861 : : }
862 : : else
863 : : {
864 : : /* halve alpha and try again */
865 : : /* add out the old step */
866 [ + - ][ + - ]: 14523 : pd.move_vertex( mAlpha * mSearch, freeVertexIndex, err );
867 : : // pd.set_coords_array_element(coords[freeVertexIndex],0,err);
868 : :
869 : : /* halve step size */
870 : 14523 : mAlpha = mAlpha / 2;
871 : : MSQ_PRINT( 3 )( "Step not accepted, the new alpha %f\n", mAlpha );
872 : :
873 [ - + ]: 14523 : if( mAlpha < minStepSize )
874 : : {
875 : : /* get back the original point, function, and active set */
876 : : MSQ_PRINT( 2 )( "Optimization Exiting: Step too small\n" );
877 [ # # ][ # # ]: 0 : pd.set_vertex_coordinates( original_point, freeVertexIndex, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
878 [ # # ]: 0 : mFunction = originalFunction;
879 [ # # ]: 0 : mActive = originalActive;
880 : 0 : optStatus = MSQ_STEP_TOO_SMALL;
881 : 0 : step_done = true;
882 : : }
883 : : else
884 : : {
885 [ + - ]: 14523 : testFunction = mFunction;
886 [ + - ]: 14523 : testActive = mActive;
887 : 27538 : previous_improvement = current_improvement;
888 : : }
889 : : }
890 : : }
891 : : }
892 [ + + ]: 13345 : if( current_improvement > 0 && optStatus == MSQ_STEP_ACCEPTED )
893 : : { MSQ_PRINT( 2 )( "Accepted a negative step %f \n", current_improvement ); }
894 : : }
895 : :
896 : 28860 : bool NonSmoothDescent::compute_function( PatchData* patch_data, std::vector< double >& func, MsqError& err )
897 : : {
898 : : // NEED 1.0/FUNCTION WHICH IS ONLY TRUE OF CONDITION NUMBER
899 : :
900 : 28860 : func.resize( qmHandles.size() );
901 : 28860 : bool valid_bool = true;
902 [ + + ]: 679176 : for( size_t i = 0; i < qmHandles.size(); i++ )
903 : : {
904 [ + - ][ + - ]: 650316 : valid_bool = valid_bool && currentQM->evaluate( *patch_data, qmHandles[i], func[i], err );
905 [ - + ][ # # ]: 650316 : MSQ_ERRZERO( err );
[ - + ]
906 : : }
907 : :
908 : 28860 : return valid_bool;
909 : : }
910 : :
911 : 29570 : bool NonSmoothDescent::compute_gradient( PatchData* patch_data, std::vector< Vector3D >& gradient_out, MsqError& err )
912 : : {
913 : : MSQ_DBGOUT( 2 ) << "Computing Gradient\n";
914 : :
915 : 14785 : bool valid = true;
916 : :
917 : : #ifdef NUMERICAL_GRADIENT
918 : :
919 : : const double delta = 10e-6;
920 : : std::vector< double > func( qmHandles.size() ), fdelta( qmHandles.size() );
921 : :
922 : : valid = this->compute_function( patch_data, func, err );
923 : : MSQ_ERRZERO( err );
924 : :
925 : : /* gradient in the x, y, z direction */
926 : :
927 : : for( int j = 0; j < 3; j++ )
928 : : {
929 : :
930 : : // perturb the coordinates of the free vertex in the j direction by delta
931 : : Vector3D delta_3( 0, 0, 0 );
932 : : Vector3D orig_pos = patch_data->vertex_by_index( freeVertexIndex );
933 : : delta_3[j] = delta;
934 : : patch_data->move_vertex( delta_3, freeVertexIndex, err );
935 : :
936 : : // compute the function at the perturbed point location
937 : : valid = valid && this->compute_function( patch_data, fdelta, err );
938 : : MSQ_ERRZERO( err );
939 : :
940 : : // compute the numerical gradient
941 : : for( size_t i = 0; i < qmHandles.size(); i++ )
942 : : {
943 : : mGradient[i][j] = ( fdelta[i] - func[i] ) / delta;
944 : : // MSQ_DEBUG_ACTION(3,{fprintf(stdout," Gradient
945 : : // value[%d][%d]=%g\n",i,j,mGradient[i][j]);});
946 : : }
947 : :
948 : : // put the coordinates back where they belong
949 : : patch_data->set_vertex_coordinates( orig_pos, freeVertexIndex, err );
950 : : }
951 : :
952 : : #else
953 : : double value;
954 [ + - ]: 14785 : gradient_out.resize( qmHandles.size() );
955 [ + + ]: 344715 : for( size_t i = 0; i < qmHandles.size(); i++ )
956 : : {
957 [ + - ][ + - ]: 329930 : valid = valid && currentQM->evaluate_with_gradient( *patch_data, qmHandles[i], value, tmpIdx, tmpGrad, err );
[ + - ][ + - ]
958 [ + - ][ - + ]: 329930 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
959 [ + - ][ + - ]: 329930 : assert( tmpIdx.size() == 1 && tmpIdx[0] == freeVertexIndex );
[ - + ]
960 [ + - ][ + - ]: 329930 : gradient_out[i] = tmpGrad[0];
[ + - ]
961 : : }
962 : : #endif
963 : :
964 : 14785 : return valid;
965 : : }
966 : :
967 : 30182 : void NonSmoothDescent::find_active_set( const std::vector< double >& function, ActiveSet& active_set )
968 : : {
969 : : // local parameter initialization
970 : 30182 : const double activeEpsilon = .3e-4;
971 : : // activeEpsilon = .3e-8;
972 : :
973 : : double function_val;
974 : : double active_value0;
975 : : double temp;
976 : :
977 : : // FUNCTION_TIMER_START("Find Active Set");
978 : : MSQ_DBGOUT( 2 ) << "\nFinding the active set\n";
979 : :
980 : : /* the first function value is our initial active value */
981 : 30182 : active_set.set( 0 );
982 : 30182 : active_set.true_active_value = function[0];
983 : : // MSQ_DEBUG_ACTION(3,{fprintf(stdout," function value[0]: %g\n",function[0]);});
984 : :
985 : : /* first sort out the active set...
986 : : all vals within active_epsilon of largest val */
987 : :
988 [ + + ]: 679866 : for( size_t i = 1; i < qmHandles.size(); i++ )
989 : : {
990 : 649684 : function_val = function[i];
991 [ + + ]: 649684 : if( active_set.true_active_value < function_val ) active_set.true_active_value = function_val;
992 : 649684 : active_value0 = function[active_set.active_ind[0]];
993 : 649684 : temp = fabs( function_val - active_value0 );
994 : : // MSQ_DEBUG_ACTION(3,{fprintf(stdout," function value[%d]: %g\n",i,function[i]);});
995 [ + + ]: 649684 : if( function_val > active_value0 )
996 : : { // seek max_i function[i]
997 [ + + ]: 81336 : if( temp >= activeEpsilon )
998 : : {
999 : 70313 : active_set.set( i ); // new max
1000 : : }
1001 : : else
1002 : : {
1003 : 81336 : active_set.add( i, fabs( function_val - active_value0 ) < EPSILON );
1004 : : }
1005 : : }
1006 : : else
1007 : : {
1008 [ + + ]: 568348 : if( temp < activeEpsilon ) { active_set.add( i, fabs( function_val - active_value0 ) < EPSILON ); }
1009 : : }
1010 : : }
1011 : 30182 : }
1012 : :
1013 : 43027 : bool NonSmoothDescent::validity_check( PatchData& pd, MsqError& err )
1014 : :
1015 : : {
1016 : : // FUNCTION_TIMER_START("Validity Check");
1017 : :
1018 : : // ONLY FOR SIMPLICIAL MESHES - THERE SHOULD BE A VALIDITY CHECKER ASSOCIATED
1019 : : // WITH MSQ ELEMENTS
1020 : :
1021 : : /* check that the simplicial mesh is still valid, based on right handedness.
1022 : : Returns a 1 or a 0 */
1023 : :
1024 : : // TODO as a first step we can switch this over to the function
1025 : : // evaluation and use the rest of the code as is
1026 : 43027 : bool valid = true;
1027 : 43027 : double dEps = 1.e-13;
1028 : :
1029 : : double x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4;
1030 [ + - ]: 43027 : const MsqVertex* coords = pd.get_vertex_array( err );
1031 : :
1032 [ + - ][ + + ]: 1012447 : for( size_t i = 0; i < pd.num_elements(); i++ )
1033 : : {
1034 [ + - ][ + - ]: 969494 : const size_t* conn = pd.element_by_index( i ).get_vertex_index_array();
1035 [ + - ]: 969494 : coords[conn[0]].get_coordinates( x1, y1, z1 );
1036 [ + - ]: 969494 : coords[conn[1]].get_coordinates( x2, y2, z2 );
1037 [ + - ]: 969494 : coords[conn[2]].get_coordinates( x3, y3, z3 );
1038 [ + - ]: 969494 : coords[conn[3]].get_coordinates( x4, y4, z4 );
1039 : :
1040 : 969494 : double dDX2 = x2 - x1;
1041 : 969494 : double dDX3 = x3 - x1;
1042 : 969494 : double dDX4 = x4 - x1;
1043 : :
1044 : 969494 : double dDY2 = y2 - y1;
1045 : 969494 : double dDY3 = y3 - y1;
1046 : 969494 : double dDY4 = y4 - y1;
1047 : :
1048 : 969494 : double dDZ2 = z2 - z1;
1049 : 969494 : double dDZ3 = z3 - z1;
1050 : 969494 : double dDZ4 = z4 - z1;
1051 : :
1052 : : /* dDet is proportional to the cell volume */
1053 : 1938988 : double dDet = dDX2 * dDY3 * dDZ4 + dDX3 * dDY4 * dDZ2 + dDX4 * dDY2 * dDZ3 - dDZ2 * dDY3 * dDX4 -
1054 : 1938988 : dDZ3 * dDY4 * dDX2 - dDZ4 * dDY2 * dDX3;
1055 : :
1056 : : /* Compute a length scale based on edge lengths. */
1057 : 1938988 : double dScale = ( sqrt( ( x1 - x2 ) * ( x1 - x2 ) + ( y1 - y2 ) * ( y1 - y2 ) + ( z1 - z2 ) * ( z1 - z2 ) ) +
1058 : 1938988 : sqrt( ( x1 - x3 ) * ( x1 - x3 ) + ( y1 - y3 ) * ( y1 - y3 ) + ( z1 - z3 ) * ( z1 - z3 ) ) +
1059 : 1938988 : sqrt( ( x1 - x4 ) * ( x1 - x4 ) + ( y1 - y4 ) * ( y1 - y4 ) + ( z1 - z4 ) * ( z1 - z4 ) ) +
1060 : 1938988 : sqrt( ( x2 - x3 ) * ( x2 - x3 ) + ( y2 - y3 ) * ( y2 - y3 ) + ( z2 - z3 ) * ( z2 - z3 ) ) +
1061 : 1938988 : sqrt( ( x2 - x4 ) * ( x2 - x4 ) + ( y2 - y4 ) * ( y2 - y4 ) + ( z2 - z4 ) * ( z2 - z4 ) ) +
1062 : 969494 : sqrt( ( x3 - x4 ) * ( x3 - x4 ) + ( y3 - y4 ) * ( y3 - y4 ) + ( z3 - z4 ) * ( z3 - z4 ) ) ) /
1063 : 969494 : 6.;
1064 : :
1065 : : /* Use the length scale to get a better idea if the tet is flat or
1066 : : just really small. */
1067 [ - + ]: 969494 : if( fabs( dScale ) < EPSILON ) { return false; }
1068 : : else
1069 : : {
1070 : 969494 : dDet /= ( dScale * dScale * dScale );
1071 : : }
1072 : :
1073 [ + + ]: 969494 : if( dDet > dEps ) { valid = true; }
1074 [ + - ]: 74 : else if( dDet < -dEps )
1075 : : {
1076 : 74 : return false;
1077 : : }
1078 : : else
1079 : : {
1080 : 0 : return false;
1081 : : }
1082 : : } // end for i=1,numElements
1083 : :
1084 : : // MSQ_DEBUG_ACTION(2,{fprintf(stdout,"Mesh Validity is: %d \n",valid);});
1085 : :
1086 : : // FUNCTION_TIMER_END();
1087 : 43027 : return ( valid );
1088 : : }
1089 : :
1090 : 16129 : void NonSmoothDescent::get_active_directions( const std::vector< Vector3D >& pGradient, std::vector< Vector3D >& dir,
1091 : : MsqError& /*err*/ )
1092 : : {
1093 : 16129 : const size_t num_active = mActive.active_ind.size();
1094 : 16129 : dir.resize( num_active );
1095 [ + + ]: 56397 : for( size_t i = 0; i < num_active; i++ )
1096 : : {
1097 : 40268 : dir[i] = pGradient[mActive.active_ind[i]];
1098 : : }
1099 : 16129 : }
1100 : :
1101 : 5694 : bool NonSmoothDescent::check_vector_dots( const std::vector< Vector3D >& vec, const Vector3D& normal,
1102 : : MsqError& /*err*/ )
1103 : : {
1104 : 5694 : double test_dot = vec[0] % normal;
1105 : 5694 : unsigned ind = 1;
1106 [ + + ][ + - ]: 9826 : while( ( fabs( test_dot ) < EPSILON ) && ( ind < vec.size() ) )
[ + + ]
1107 : : {
1108 : 4132 : test_dot = vec[ind] % normal;
1109 : 4132 : ind++;
1110 : : }
1111 : :
1112 [ + + ]: 10906 : for( unsigned i = ind; i < vec.size(); i++ )
1113 : : {
1114 : 7359 : double dot = vec[i] % normal;
1115 [ + + ][ + + ]: 7359 : if( ( ( dot > 0 && test_dot < 0 ) || ( dot < 0 && test_dot > 0 ) ) && ( fabs( dot ) > EPSILON ) )
[ + + ][ + + ]
[ + + ]
1116 : 2147 : { return true; }
1117 : : }
1118 : 3547 : return false;
1119 : : }
1120 : :
1121 : 8653 : bool NonSmoothDescent::convex_hull_test( const std::vector< Vector3D >& vec, MsqError& err )
1122 : : {
1123 : : // int ierr;
1124 : 8653 : bool equil = false;
1125 : 8653 : Direction dir_done = MSQ_ZDIR;
1126 : 8653 : Status status = MSQ_CHECK_Z_COORD_DIRECTION;
1127 [ + - ][ + - ]: 8653 : Vector3D pt1, pt2, pt3, normal;
[ + - ][ + - ]
1128 : :
1129 : : /* tries to determine equilibrium for the 3D case */
1130 : :
1131 [ + + ]: 8653 : if( vec.size() <= 2 ) status = MSQ_NO_EQUIL;
1132 : :
1133 [ + + ][ + + ]: 14442 : while( ( status != MSQ_EQUIL ) && ( status != MSQ_NO_EQUIL ) && ( status != MSQ_HULL_TEST_ERROR ) )
[ + - ]
1134 : : {
1135 [ + + ]: 5789 : if( status == MSQ_CHECK_Z_COORD_DIRECTION )
1136 : : {
1137 [ + - ]: 4141 : this->find_plane_points( MSQ_ZDIR, MSQ_YDIR, vec, pt1, pt2, pt3, status, err );
1138 : 4141 : dir_done = MSQ_ZDIR;
1139 : : }
1140 [ + + ]: 5789 : if( status == MSQ_CHECK_Y_COORD_DIRECTION )
1141 : : {
1142 [ + - ]: 1026 : this->find_plane_points( MSQ_YDIR, MSQ_XDIR, vec, pt1, pt2, pt3, status, err );
1143 : 1026 : dir_done = MSQ_YDIR;
1144 : : }
1145 [ + + ]: 5789 : if( status == MSQ_CHECK_X_COORD_DIRECTION )
1146 : : {
1147 [ + - ]: 622 : this->find_plane_points( MSQ_XDIR, MSQ_ZDIR, vec, pt1, pt2, pt3, status, err );
1148 : 622 : dir_done = MSQ_XDIR;
1149 : : }
1150 [ + + ][ + - ]: 5789 : if( status == MSQ_TWO_PT_PLANE ) { pt3 = Vector3D( 0, 0, 0 ); }
[ + - ]
1151 [ + + ][ - + ]: 5789 : if( ( status == MSQ_TWO_PT_PLANE ) || ( status == MSQ_THREE_PT_PLANE ) )
1152 : : {
1153 [ + - ]: 5694 : this->find_plane_normal( pt1, pt2, pt3, normal, err );
1154 [ + - ]: 5694 : equil = this->check_vector_dots( vec, normal, err );
1155 [ + + ]: 5694 : if( equil )
1156 : : {
1157 [ + + + - ]: 2147 : switch( dir_done )
1158 : : {
1159 : : case MSQ_ZDIR:
1160 : 1026 : equil = false;
1161 : 1026 : status = MSQ_CHECK_Y_COORD_DIRECTION;
1162 : 1026 : break;
1163 : : case MSQ_YDIR:
1164 : 622 : equil = false;
1165 : 622 : status = MSQ_CHECK_X_COORD_DIRECTION;
1166 : 622 : break;
1167 : : case MSQ_XDIR:
1168 : 499 : equil = true;
1169 : 2147 : status = MSQ_EQUIL;
1170 : : }
1171 : : }
1172 [ + - ]: 3547 : else if( equil == 0 )
1173 : : {
1174 : 3547 : status = MSQ_NO_EQUIL;
1175 : : }
1176 : : else
1177 : : {
1178 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "equil flag not set to 0 or 1", MsqError::INVALID_STATE );
1179 : : }
1180 : : }
1181 : : }
1182 [ + + - ]: 8653 : switch( status )
1183 : : {
1184 : : case MSQ_NO_EQUIL:
1185 : : MSQ_PRINT( 3 )( "Not an equilibrium point\n" );
1186 : 8154 : equil = false;
1187 : 8154 : break;
1188 : : case MSQ_EQUIL:
1189 : : MSQ_PRINT( 3 )( "An equilibrium point\n" );
1190 : 499 : equil = true;
1191 : 499 : break;
1192 : : default:
1193 : : MSQ_PRINT( 3 )( "Failed to determine equil or not; status = %d\n", status );
1194 : : }
1195 : : // FUNCTION_TIMER_END();
1196 : 8653 : return ( equil );
1197 : : }
1198 : :
1199 : 7476 : void NonSmoothDescent::form_grammian( const std::vector< Vector3D >& vec, MsqError& /*err*/ )
1200 : : {
1201 : : /* form the grammian with the dot products of the gradients */
1202 : 7476 : const size_t num_active = mActive.active_ind.size();
1203 : 7476 : mG.resize( num_active );
1204 [ + + ]: 25706 : for( size_t i = 0; i < num_active; i++ )
1205 [ + + ]: 50637 : for( size_t j = i; j < num_active; j++ )
1206 : 32407 : mG( i, j ) = vec[i] % vec[j];
1207 : 7476 : }
1208 : :
1209 : 8653 : bool NonSmoothDescent::check_equilibrium( OptStatus& status, MsqError& err )
1210 : : {
1211 [ + - ]: 8653 : std::vector< Vector3D > dir;
1212 : :
1213 : : // TODO - this subroutine is no longer clear to me... is it still
1214 : : // appropriate for quads and hexes? I think it might be in 2D, but
1215 : : // 3D is less clear. Is there a more general algorithm to use?
1216 : : // ask Todd/check in numerical optimization
1217 : :
1218 : 8653 : bool equil = false;
1219 : 8653 : const size_t num_active = mActive.active_ind.size();
1220 : : ;
1221 : :
1222 [ - + ]: 8653 : if( num_active == qmHandles.size() )
1223 : : {
1224 : 0 : equil = true;
1225 : 0 : status = MSQ_EQUILIBRIUM;
1226 : : MSQ_PRINT( 3 )( "All the function values are in the active set\n" );
1227 : : }
1228 : :
1229 : : /* set up the active mGradient directions */
1230 [ + - ]: 8653 : this->get_active_directions( mGradient, dir, err );
1231 [ + - ][ - + ]: 8653 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
1232 : :
1233 : : /* normalize the active directions */
1234 [ + + ]: 30691 : for( size_t j = 0; j < num_active; j++ )
1235 [ + - ][ + - ]: 22038 : dir[j] /= dir[j].length();
[ + - ][ + - ]
1236 : :
1237 [ + - ]: 8653 : equil = this->convex_hull_test( dir, err );
1238 [ + + ]: 8653 : if( equil ) status = MSQ_EQUILIBRIUM;
1239 : :
1240 : 8653 : return equil;
1241 : : }
1242 : :
1243 : 6125 : static double condition3x3( const double A[9] )
1244 : : {
1245 : : // int ierr;
1246 : : double a11, a12, a13;
1247 : : double a21, a22, a23;
1248 : : double a31, a32, a33;
1249 : : // double s1, s2, s4, s3, t0;
1250 : : double s1, s2, s3;
1251 : : double denom;
1252 : : // double one = 1.0;
1253 : : double temp;
1254 : 6125 : bool zero_denom = true;
1255 : :
1256 : 6125 : a11 = A[0];
1257 : 6125 : a12 = A[1];
1258 : 6125 : a13 = A[2];
1259 : 6125 : a21 = A[3];
1260 : 6125 : a22 = A[4];
1261 : 6125 : a23 = A[5];
1262 : 6125 : a31 = A[6];
1263 : 6125 : a32 = A[7];
1264 : 6125 : a33 = A[8];
1265 : :
1266 : 6125 : denom = -a11 * a22 * a33 + a11 * a23 * a32 + a21 * a12 * a33 - a21 * a13 * a32 - a31 * a12 * a23 + a31 * a13 * a22;
1267 : :
1268 [ + - ][ + - ]: 6125 : if( ( fabs( a11 ) > EPSILON ) && ( fabs( denom / a11 ) > EPSILON ) ) { zero_denom = false; }
1269 [ + - ][ + - ]: 6125 : if( ( fabs( a22 ) > EPSILON ) && ( fabs( denom / a22 ) > EPSILON ) ) { zero_denom = false; }
1270 [ + - ][ + - ]: 6125 : if( ( fabs( a33 ) > EPSILON ) && ( fabs( denom / a33 ) > EPSILON ) ) { zero_denom = false; }
1271 : :
1272 [ - + ]: 6125 : if( zero_denom ) { return HUGE_VAL; }
1273 : : else
1274 : : {
1275 : 6125 : s1 = sqrt( a11 * a11 + a12 * a12 + a13 * a13 + a21 * a21 + a22 * a22 + a23 * a23 + a31 * a31 + a32 * a32 +
1276 : 6125 : a33 * a33 );
1277 : :
1278 : 6125 : temp = ( -a22 * a33 + a23 * a32 ) / denom;
1279 : 6125 : s3 = temp * temp;
1280 : 6125 : temp = ( a12 * a33 - a13 * a32 ) / denom;
1281 : 6125 : s3 += temp * temp;
1282 : 6125 : temp = ( a12 * a23 - a13 * a22 ) / denom;
1283 : 6125 : s3 += temp * temp;
1284 : 6125 : temp = ( a21 * a33 - a23 * a31 ) / denom;
1285 : 6125 : s3 += temp * temp;
1286 : 6125 : temp = ( a11 * a33 - a13 * a31 ) / denom;
1287 : 6125 : s3 += temp * temp;
1288 : 6125 : temp = ( a11 * a23 - a13 * a21 ) / denom;
1289 : 6125 : s3 += temp * temp;
1290 : 6125 : temp = ( a21 * a32 - a22 * a31 ) / denom;
1291 : 6125 : s3 += temp * temp;
1292 : 6125 : temp = ( -a11 * a32 + a12 * a31 ) / denom;
1293 : 6125 : s3 += temp * temp;
1294 : 6125 : temp = ( -a11 * a22 + a12 * a21 ) / denom;
1295 : 6125 : s3 += temp * temp;
1296 : :
1297 : 6125 : s2 = sqrt( s3 );
1298 : 6125 : return s1 * s2;
1299 : : }
1300 : : }
1301 : :
1302 : 2991 : double NonSmoothDescent::SymmetricMatrix::condition3x3() const
1303 : : {
1304 [ + - ][ + - ]: 2991 : double values[9] = { operator()( 0, 0 ), operator()( 0, 1 ), operator()( 0, 2 ),
[ + - ]
1305 [ + - ][ + - ]: 2991 : operator()( 1, 0 ), operator()( 1, 1 ), operator()( 1, 2 ),
[ + - ]
1306 [ + - ][ + - ]: 2991 : operator()( 2, 0 ), operator()( 2, 1 ), operator()( 2, 2 ) };
[ + - ]
1307 : 2991 : return MBMesquite::condition3x3( values );
1308 : : }
1309 : :
1310 : 10610 : void NonSmoothDescent::singular_test( int n, const Matrix3D& A, bool& singular, MsqError& err )
1311 : : {
1312 : : // int test;
1313 : : // double determinant;
1314 : : double cond;
1315 : :
1316 [ + - ][ - + ]: 10610 : if( ( n > 3 ) || ( n < 1 ) )
1317 : : {
1318 [ # # ]: 0 : MSQ_SETERR( err )( "Singular test works only for n=1 to n=3", MsqError::INVALID_ARG );
1319 : 10610 : return;
1320 : : }
1321 : :
1322 : 10610 : singular = true;
1323 [ - + + - ]: 10610 : switch( n )
1324 : : {
1325 : : case 1:
1326 [ # # ]: 0 : if( A[0][0] > 0 ) singular = false;
1327 : 0 : break;
1328 : : case 2:
1329 [ + - ]: 7476 : if( fabs( A[0][0] * A[1][1] - A[0][1] * A[1][0] ) > EPSILON ) singular = false;
1330 : 7476 : break;
1331 : : case 3:
1332 : : /* calculate the condition number */
1333 : 3134 : cond = condition3x3( A[0] );
1334 [ + - ]: 3134 : if( cond < 1E14 ) singular = false;
1335 : 3134 : break;
1336 : : }
1337 : : }
1338 : :
1339 : 7476 : void NonSmoothDescent::form_PD_grammian( MsqError& err )
1340 : : {
1341 : : int i, j, k;
1342 : : int g_ind_1;
1343 : 7476 : bool singular = false;
1344 : :
1345 : 7476 : const int num_active = mActive.active_ind.size();
1346 : :
1347 : : /* this assumes the grammian has been formed */
1348 [ + + ]: 25706 : for( i = 0; i < num_active; i++ )
1349 : : {
1350 [ + + ]: 50637 : for( j = i; j < num_active; j++ )
1351 : : {
1352 [ + - ][ - + ]: 32407 : if( mG( i, j ) == -1 )
1353 : : {
1354 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "Grammian not computed properly", MsqError::INVALID_STATE );
1355 : 7476 : return;
1356 : : }
1357 : : }
1358 : : }
1359 : :
1360 : : /* use the first gradient in the active set */
1361 : 7476 : g_ind_1 = 0;
1362 [ + - ][ + - ]: 7476 : mPDG[0][0] = mG( 0, 0 );
1363 [ + - ]: 7476 : pdgInd[0] = mActive.active_ind[0];
1364 : :
1365 : : /* test the rest and add them as appropriate */
1366 : 7476 : k = 1;
1367 : 7476 : i = 1;
1368 [ + + ][ + + ]: 18086 : while( ( k < 3 ) && ( i < num_active ) )
1369 : : {
1370 [ + - ][ + - ]: 10610 : mPDG[0][k] = mPDG[k][0] = mG( 0, i );
[ + - ]
1371 [ + - ][ + - ]: 10610 : mPDG[k][k] = mG( i, i );
1372 [ + + ]: 10610 : if( k == 2 )
1373 : : { /* add the dot product of g1 and g2 */
1374 [ + - ][ + - ]: 3134 : mPDG[1][k] = mPDG[k][1] = mG( g_ind_1, i );
[ + - ]
1375 : : }
1376 [ + - ]: 10610 : this->singular_test( k + 1, mPDG, singular, err );
1377 [ + - ]: 10610 : if( !singular )
1378 : : {
1379 [ + - ]: 10610 : pdgInd[k] = mActive.active_ind[i];
1380 [ + + ]: 10610 : if( k == 1 ) g_ind_1 = i;
1381 : 10610 : k++;
1382 : : }
1383 : 10610 : i++;
1384 : : }
1385 : : }
1386 : :
1387 : 143 : void NonSmoothDescent::search_edges_faces( const Vector3D* dir, Vector3D& mSearch, MsqError& /*err*/ )
1388 : : {
1389 : : bool viable;
1390 : : double a, b, denom;
1391 [ + - ]: 143 : Vector3D g_bar;
1392 [ + - ]: 143 : Vector3D temp_search( 0, 0, 0 ); /* initialize the search direction to 0,0 */
1393 : : double projection, min_projection;
1394 : :
1395 : 143 : const size_t num_active = mActive.active_ind.size();
1396 : : ;
1397 : :
1398 : : /* Check for viable faces */
1399 : 143 : min_projection = HUGE_VAL;
1400 [ + + ]: 716 : for( size_t i = 0; i < num_active; i++ )
1401 : : {
1402 : : /* FACE I */
1403 : 573 : viable = true;
1404 : :
1405 : : /* test the viability */
1406 [ + + ]: 2870 : for( size_t j = 0; j < num_active; j++ )
1407 : : { /* lagrange multipliers>0 */
1408 [ + - ][ + + ]: 2297 : if( mG( j, i ) < 0 ) viable = false;
1409 : : }
1410 : :
1411 : : /* find the minimum of viable directions */
1412 [ + + ][ + - ]: 573 : if( ( viable ) && ( mG( i, i ) < min_projection ) )
[ + + ][ + + ]
1413 : : {
1414 [ + - ]: 9 : min_projection = mG( i, i );
1415 [ + - ]: 9 : temp_search = dir[i];
1416 [ + - ]: 9 : mSteepest = mActive.active_ind[i];
1417 : : }
1418 : :
1419 : : /* INTERSECTION IJ */
1420 [ + + ]: 1435 : for( size_t j = i + 1; j < num_active; j++ )
1421 : : {
1422 : 862 : viable = true;
1423 : :
1424 : : /* find the coefficients of the intersection
1425 : : and test the viability */
1426 [ + - ][ + - ]: 862 : denom = 2 * mG( i, j ) - mG( i, i ) - mG( j, j );
[ + - ]
1427 : 862 : a = b = 0;
1428 [ + - ]: 862 : if( fabs( denom ) > EPSILON )
1429 : : {
1430 [ + - ][ + - ]: 862 : b = ( mG( i, j ) - mG( i, i ) ) / denom;
1431 : 862 : a = 1 - b;
1432 [ + + ][ + + ]: 862 : if( ( b < 0 ) || ( b > 1 ) ) /* 0 < b < 1 */
1433 : 98 : viable = false;
1434 [ + + ]: 4320 : for( size_t k = 0; k < num_active; k++ )
1435 : : { /* lagrange multipliers>0 */
1436 [ + - ][ + - ]: 3458 : if( ( a * mG( k, i ) + b * mG( k, j ) ) <= 0 ) viable = false;
[ + + ]
1437 : : }
1438 : : }
1439 : : else
1440 : : {
1441 : 0 : viable = false; /* Linearly dependent */
1442 : : }
1443 : :
1444 : : /* find the minimum of viable directions */
1445 [ + + ]: 862 : if( viable )
1446 : : {
1447 [ + - ][ + - ]: 118 : g_bar = a * dir[i] + b * dir[j];
[ + - ][ + - ]
1448 [ + - ]: 118 : projection = g_bar % g_bar;
1449 [ + + ]: 118 : if( projection < min_projection )
1450 : : {
1451 : 65 : min_projection = projection;
1452 [ + - ]: 65 : temp_search = g_bar;
1453 [ + - ]: 65 : mSteepest = mActive.active_ind[i];
1454 : : }
1455 : : }
1456 : : }
1457 : : }
1458 [ + - ][ + - ]: 143 : if( optStatus != MSQ_EQUILIBRIUM ) { mSearch = temp_search; }
1459 : 143 : }
1460 : :
1461 : 2991 : bool NonSmoothDescent::solve2x2( double a11, double a12, double a21, double a22, double b1, double b2, double x[2],
1462 : : MsqError& /*err*/ )
1463 : : {
1464 : : double factor;
1465 : :
1466 : : /* if the system is not singular, solve it */
1467 [ + - ]: 2991 : if( fabs( a11 * a22 - a21 * a12 ) > EPSILON )
1468 : : {
1469 [ + - ]: 2991 : if( fabs( a11 ) > EPSILON )
1470 : : {
1471 : 2991 : factor = ( a21 / a11 );
1472 : 2991 : x[1] = ( b2 - factor * b1 ) / ( a22 - factor * a12 );
1473 : 2991 : x[0] = ( b1 - a12 * x[1] ) / a11;
1474 : : }
1475 [ # # ]: 0 : else if( fabs( a21 ) > EPSILON )
1476 : : {
1477 : 0 : factor = ( a11 / a21 );
1478 : 0 : x[1] = ( b1 - factor * b2 ) / ( a12 - factor * a22 );
1479 : 0 : x[0] = ( b2 - a22 * x[1] ) / a21;
1480 : : }
1481 : 2991 : return true;
1482 : : }
1483 : : else
1484 : : {
1485 : 0 : return false;
1486 : : }
1487 : : }
1488 : :
1489 : 2991 : void NonSmoothDescent::form_reduced_matrix( SymmetricMatrix& P )
1490 : : {
1491 : 2991 : const size_t P_size = mActive.active_ind.size() - 1;
1492 : 2991 : P.resize( P_size );
1493 [ + + ]: 8973 : for( size_t i = 0; i < P_size; i++ )
1494 : : {
1495 : 5982 : P( i, i ) = mG( 0, 0 ) - 2 * mG( 0, i + 1 ) + mG( i + 1, i + 1 );
1496 [ + + ]: 8973 : for( size_t j = i + 1; j < P_size; j++ )
1497 : : {
1498 : 2991 : P( i, j ) = mG( 0, 0 ) - mG( 0, j + 1 ) - mG( i + 1, 0 ) + mG( i + 1, j + 1 );
1499 : : }
1500 : : }
1501 : 2991 : }
1502 : :
1503 : 27538 : void NonSmoothDescent::get_min_estimate( double* final_est, MsqError& /*err*/ )
1504 : : {
1505 : : double est_imp;
1506 : :
1507 : 27538 : *final_est = -HUGE_VAL;
1508 [ + + ]: 77254 : for( size_t i = 0; i < mActive.active_ind.size(); i++ )
1509 : : {
1510 : 49716 : est_imp = -mAlpha * mGS[mActive.active_ind[i]];
1511 [ + + ]: 49716 : if( est_imp > *final_est ) *final_est = est_imp;
1512 : : }
1513 [ - + ]: 27538 : if( *final_est == 0 )
1514 : : {
1515 : 0 : *final_est = -HUGE_VAL;
1516 [ # # ]: 0 : for( size_t i = 0; i < qmHandles.size(); i++ )
1517 : : {
1518 : 0 : est_imp = -mAlpha * mGS[i];
1519 [ # # ][ # # ]: 0 : if( ( est_imp > *final_est ) && ( fabs( est_imp ) > EPSILON ) ) { *final_est = est_imp; }
1520 : : }
1521 : : }
1522 : 27538 : }
1523 : :
1524 : 13345 : void NonSmoothDescent::get_gradient_projections( const Vector3D& mSearch, MsqError& /*err*/ )
1525 : : {
1526 [ + + ]: 311163 : for( size_t i = 0; i < qmHandles.size(); i++ )
1527 : 297818 : mGS[i] = mGradient[i] % mSearch;
1528 : :
1529 : : MSQ_PRINT( 3 )( "steepest in get_gradient_projections %lu\n", (unsigned long)mSteepest );
1530 : 13345 : }
1531 : :
1532 : 13345 : void NonSmoothDescent::compute_alpha( MsqError& /*err*/ )
1533 : : {
1534 : : double steepest_function;
1535 : : double steepest_grad;
1536 : : double alpha_i;
1537 : 13345 : double min_positive_value = HUGE_VAL;
1538 : :
1539 : : // FUNCTION_TIMER_START("Compute Alpha");
1540 : :
1541 : : MSQ_PRINT( 2 )( "In compute alpha\n" );
1542 : :
1543 : 13345 : mAlpha = HUGE_VAL;
1544 : :
1545 : 13345 : steepest_function = mFunction[mSteepest];
1546 : 13345 : steepest_grad = mGS[mSteepest];
1547 [ + + ]: 311163 : for( size_t i = 0; i < qmHandles.size(); i++ )
1548 : : {
1549 : : /* if it's not active */
1550 [ + + ]: 297818 : if( i != mSteepest )
1551 : : {
1552 : 284473 : alpha_i = steepest_function - mFunction[i];
1553 : :
1554 [ + + ]: 284473 : if( fabs( mGS[mSteepest] - mGS[i] ) > 1E-13 )
1555 : : {
1556 : : /* compute line intersection */
1557 : 275000 : alpha_i = alpha_i / ( steepest_grad - mGS[i] );
1558 : : }
1559 : : else
1560 : : {
1561 : : /* the lines don't intersect - it's not under consideration*/
1562 : 9473 : alpha_i = 0;
1563 : : }
1564 [ + + ][ + + ]: 284473 : if( ( alpha_i > minStepSize ) && ( fabs( alpha_i ) < fabs( mAlpha ) ) )
1565 : : {
1566 : 45609 : mAlpha = fabs( alpha_i );
1567 : : MSQ_PRINT( 3 )( "Setting alpha %lu %g\n", (unsigned long)i, alpha_i );
1568 : : }
1569 [ + + ][ + + ]: 284473 : if( ( alpha_i > 0 ) && ( alpha_i < min_positive_value ) ) { min_positive_value = alpha_i; }
1570 : : }
1571 : : }
1572 : :
1573 [ - + ][ # # ]: 13345 : if( ( mAlpha == HUGE_VAL ) && ( min_positive_value != HUGE_VAL ) ) { mAlpha = min_positive_value; }
1574 : :
1575 : : /* if it never gets set, set it to the default */
1576 [ - + ]: 13345 : if( mAlpha == HUGE_VAL )
1577 : : {
1578 : 0 : mAlpha = maxAlpha;
1579 : : MSQ_PRINT( 3 )( "Setting alpha to the maximum step length\n" );
1580 : : }
1581 : :
1582 : : MSQ_PRINT( 3 )( " The initial step size: %f\n", mAlpha );
1583 : :
1584 : : // FUNCTION_TIMER_END();
1585 : 13345 : }
1586 : :
1587 : 13463 : void NonSmoothDescent::print_active_set( const ActiveSet& active_set, const std::vector< double >& func )
1588 : : {
1589 : 13463 : if( active_set.active_ind.empty() ) MSQ_DBGOUT( 3 ) << "No active values\n";
1590 : : /* print the active set */
1591 [ + + ]: 37680 : for( size_t i = 0; i < active_set.active_ind.size(); i++ )
1592 : : {
1593 : : MSQ_PRINT( 3 )
1594 : : ( "Active value %lu: %f \n", (unsigned long)i + 1, func[active_set.active_ind[i]] );
1595 : : }
1596 : 13463 : }
1597 : :
1598 : 1322 : void NonSmoothDescent::init_opt( PatchData& pd, MsqError& err )
1599 : : {
1600 : 1322 : qmHandles.clear();
1601 [ - + ][ # # ]: 2644 : currentQM->get_evaluations( pd, qmHandles, true, err );MSQ_ERRRTN( err );
[ - + ]
1602 : :
1603 : : MSQ_PRINT( 2 )( "\nInitializing Optimization \n" );
1604 : :
1605 : : /* for the purposes of initialization will be set to zero after */
1606 : 1322 : optStatus = MSQ_NO_STATUS;
1607 : 1322 : mSteepest = 0;
1608 : 1322 : mAlpha = 0;
1609 : 1322 : maxAlpha = 0;
1610 : :
1611 : : MSQ_PRINT( 3 )( " Initialized Constants \n" );
1612 : 1322 : pdgInd[0] = pdgInd[1] = pdgInd[2] = -1;
1613 [ + - ]: 1322 : mPDG = Matrix3D( 0.0 );
1614 : :
1615 : : MSQ_PRINT( 3 )( " Initialized search and PDG \n" );
1616 : 1322 : mFunction.clear();
1617 [ + - ]: 1322 : mFunction.resize( qmHandles.size(), 0.0 );
1618 : 1322 : testFunction.clear();
1619 [ + - ]: 1322 : testFunction.resize( qmHandles.size(), 0.0 );
1620 : 1322 : originalFunction.clear();
1621 [ + - ]: 1322 : originalFunction.resize( qmHandles.size(), 0.0 );
1622 : 1322 : mGradient.clear();
1623 [ + - ]: 1322 : mGradient.resize( qmHandles.size(), Vector3D( 0.0 ) );
1624 : 1322 : mGS.resize( qmHandles.size() );
1625 : :
1626 : : MSQ_PRINT( 3 )( " Initialized function/gradient \n" );
1627 : 1322 : mG.resize( qmHandles.size() );
1628 : 1322 : mG.fill( -1 );
1629 : : MSQ_PRINT( 3 )( " Initialized G\n" );
1630 : :
1631 : 1322 : prevActiveValues.clear();
1632 : 1322 : prevActiveValues.reserve( 32 );
1633 : : MSQ_PRINT( 3 )( " Initialized prevActiveValues\n" );
1634 : : }
1635 : :
1636 : 1322 : void NonSmoothDescent::init_max_step_length( PatchData& pd, MsqError& err )
1637 : : {
1638 : : size_t i, j;
1639 : 1322 : double max_diff = 0;
1640 : 1322 : double diff = 0;
1641 : :
1642 : : MSQ_PRINT( 2 )( "In init_max_step_length\n" );
1643 : :
1644 : : /* check that the input data is correct */
1645 [ - + ]: 1322 : if( pd.num_elements() == 0 )
1646 : : {
1647 [ # # ]: 0 : MSQ_SETERR( err )( "Num incident vtx = 0\n", MsqError::INVALID_MESH );
1648 : 0 : return;
1649 : : }
1650 : :
1651 : : /* find the maximum distance between two incident vertex locations */
1652 : 1322 : const MsqVertex* coords = pd.get_vertex_array( err );
1653 [ + + ]: 18741 : for( i = 0; i < pd.num_nodes() - 1; i++ )
1654 : : {
1655 [ + + ]: 164474 : for( j = i; j < pd.num_nodes(); j++ )
1656 : : {
1657 [ + - ]: 147055 : diff = ( coords[i] - coords[j] ).length_squared();
1658 [ + + ]: 147055 : if( max_diff < diff ) max_diff = diff;
1659 : : }
1660 : : }
1661 : 1322 : max_diff = sqrt( max_diff );
1662 [ - + ]: 1322 : if( max_diff == 0 )
1663 : : {
1664 : : MSQ_SETERR( err )
1665 [ # # ]: 0 : ( "Maximum distance between incident vertices = 0\n", MsqError::INVALID_MESH );
1666 : 0 : return;
1667 : : }
1668 : 1322 : maxAlpha = max_diff / 100;
1669 : :
1670 : : MSQ_PRINT( 3 )( " Maximum step is %g\n", maxAlpha );
1671 : : }
1672 : :
1673 [ + - ][ + - ]: 4 : } // namespace MBMesquite
|
