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 : : (2006) [email protected]
27 : :
28 : : ***************************************************************** */
29 : : /*!
30 : : \file AveragingQM.cpp
31 : : \brief
32 : :
33 : : \author Michael Brewer
34 : : \author Thomas Leurent
35 : : \author Jason Kraftcheck
36 : : \date 2002-05-14
37 : : */
38 : :
39 : : #include "AveragingQM.hpp"
40 : : #include "MsqVertex.hpp"
41 : : #include "MsqMeshEntity.hpp"
42 : : #include "MsqDebug.hpp"
43 : : #include "MsqTimer.hpp"
44 : : #include "PatchData.hpp"
45 : :
46 : : namespace MBMesquite
47 : : {
48 : :
49 : 16620 : double AveragingQM::average_corner_gradients( EntityTopology type, uint32_t fixed_vertices, unsigned num_corner,
50 : : double corner_values[], const Vector3D corner_grads[],
51 : : Vector3D vertex_grads[], MsqError& err )
52 : : {
53 [ + - ]: 16620 : const unsigned num_vertex = TopologyInfo::corners( type );
54 [ + - ]: 16620 : const unsigned dim = TopologyInfo::dimension( type );
55 : 16620 : const unsigned per_vertex = dim + 1;
56 : :
57 : : unsigned i, j, num_adj;
58 : : const unsigned *adj_idx, *rev_idx;
59 : :
60 : : // NOTE: This function changes the corner_values array such that
61 : : // it contains the gradient coefficients.
62 [ + - ]: 16620 : double avg = average_metric_and_weights( corner_values, num_corner, err );
63 [ + - ][ - + ]: 16620 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
64 : :
65 [ + + ]: 140066 : for( i = 0; i < num_vertex; ++i )
66 : : {
67 [ + + ]: 123446 : if( fixed_vertices & ( 1 << i ) ) // skip fixed vertices
68 : 41754 : continue;
69 : :
70 [ + - ]: 81692 : adj_idx = TopologyInfo::adjacent_vertices( type, i, num_adj );
71 [ + - ]: 81692 : rev_idx = TopologyInfo::reverse_vertex_adjacency_offsets( type, i, num_adj );
72 [ + + ]: 81692 : if( i < num_corner ) // not all vertices are corners (e.g. pyramid)
73 [ + - ][ + - ]: 79320 : vertex_grads[i] = corner_values[i] * corner_grads[per_vertex * i];
74 : : else
75 [ + - ]: 2372 : vertex_grads[i] = 0;
76 [ + + ]: 327022 : for( j = 0; j < num_adj; ++j )
77 : : {
78 : 245330 : const unsigned v = adj_idx[j], c = rev_idx[j] + 1;
79 [ + + ]: 245330 : if( v >= num_corner ) // if less corners than vertices (e.g. pyramid apex)
80 : 5171 : continue;
81 [ + - ][ + - ]: 240159 : vertex_grads[i] += corner_values[v] * corner_grads[per_vertex * v + c];
82 : : }
83 : : }
84 : :
85 : 16620 : return avg;
86 : : }
87 : :
88 : : /**\brief Iterate over only diagonal blocks of element corner Hessian data
89 : : *
90 : : * Given concatenation of corner Hessian data for an element, iterate
91 : : * over only the diagonal terms for each corner. This class allows
92 : : * common code to be used to generate Hessian diagonal blocks from either
93 : : * the diagonal blocks for each corner or the full Hessian data for each
94 : : * corner, where this class is used for the latter.
95 : : */
96 : : class CornerHessDiagIterator
97 : : {
98 : : private:
99 : : const Matrix3D* cornerHess; //!< Current location in concatenated Hessian data.
100 : : const EntityTopology elemType; //!< Element topology for Hessian data
101 : : unsigned mCorner; //!< The element corner for which cornerHess
102 : : //!< is pointing into the corresponding Hessian data.
103 : : unsigned mStep; //!< Amount to step to reach next diagonal block.
104 : : public:
105 : 0 : CornerHessDiagIterator( const Matrix3D* corner_hessians, EntityTopology elem_type )
106 : 0 : : cornerHess( corner_hessians ), elemType( elem_type ), mCorner( 0 )
107 : : {
108 : 0 : TopologyInfo::adjacent_vertices( elemType, mCorner, mStep );
109 : 0 : ++mStep;
110 : 0 : }
111 : :
112 : 0 : SymMatrix3D operator*() const
113 : : {
114 : 0 : return cornerHess->upper();
115 : : }
116 : :
117 : 0 : CornerHessDiagIterator& operator++()
118 : : {
119 : 0 : cornerHess += mStep;
120 [ # # ]: 0 : if( !--mStep )
121 : : {
122 : 0 : TopologyInfo::adjacent_vertices( elemType, ++mCorner, mStep );
123 : 0 : ++mStep;
124 : : }
125 : 0 : return *this;
126 : : }
127 : :
128 : : CornerHessDiagIterator operator++( int )
129 : : {
130 : : CornerHessDiagIterator copy( *this );
131 : : operator++();
132 : : return copy;
133 : : }
134 : : };
135 : :
136 : : template < typename HessIter >
137 : 0 : static inline double sum_corner_diagonals( EntityTopology type, unsigned num_corner, const double corner_values[],
138 : : const Vector3D corner_grads[], HessIter corner_diag_blocks,
139 : : Vector3D vertex_grads[], SymMatrix3D vertex_hessians[] )
140 : : {
141 : : unsigned i, n, r, R, idx[4];
142 : : const unsigned* adj_list;
143 : 0 : double avg = 0.0;
144 : :
145 : : // calculate mean
146 [ # # ][ # # ]: 0 : for( i = 0; i < num_corner; ++i )
147 : 0 : avg += corner_values[i];
148 : :
149 : 0 : const Vector3D* grad = corner_grads;
150 : 0 : HessIter hess = corner_diag_blocks;
151 [ # # ][ # # ]: 0 : for( i = 0; i < num_corner; ++i )
152 : : {
153 [ # # ][ # # ]: 0 : adj_list = TopologyInfo::adjacent_vertices( type, i, n );
154 : 0 : idx[0] = i;
155 : 0 : idx[1] = adj_list[0];
156 : 0 : idx[2] = adj_list[1];
157 : 0 : idx[3] = adj_list[2 % n]; // %n so don't read off end if 2D
158 : :
159 [ # # ][ # # ]: 0 : for( r = 0; r <= n; ++r )
160 : : {
161 : 0 : R = idx[r];
162 [ # # ][ # # ]: 0 : vertex_grads[R] += *grad;
163 [ # # ][ # # ]: 0 : vertex_hessians[R] += *hess;
[ # # ]
164 : 0 : ++grad;
165 [ # # ]: 0 : ++hess;
166 : : }
167 : : }
168 : 0 : return avg;
169 : : }
170 : :
171 : : template < typename HessIter >
172 : 0 : static inline double sum_sqr_corner_diagonals( EntityTopology type, unsigned num_corner, const double corner_values[],
173 : : const Vector3D corner_grads[], HessIter corner_diag_blocks,
174 : : Vector3D vertex_grads[], SymMatrix3D vertex_hessians[] )
175 : : {
176 : : unsigned i, n, r, R, idx[4];
177 : : const unsigned* adj_list;
178 : 0 : double v, avg = 0.0;
179 : :
180 : : // calculate mean
181 [ # # ][ # # ]: 0 : for( i = 0; i < num_corner; ++i )
182 : 0 : avg += corner_values[i] * corner_values[i];
183 : :
184 : 0 : const Vector3D* grad = corner_grads;
185 : 0 : HessIter hess = corner_diag_blocks;
186 [ # # ][ # # ]: 0 : for( i = 0; i < num_corner; ++i )
187 : : {
188 [ # # ][ # # ]: 0 : adj_list = TopologyInfo::adjacent_vertices( type, i, n );
189 : 0 : idx[0] = i;
190 : 0 : idx[1] = adj_list[0];
191 : 0 : idx[2] = adj_list[1];
192 : 0 : idx[3] = adj_list[2 % n]; // %n so don't read off end if 2D
193 : 0 : ++n;
194 : :
195 : 0 : v = 2.0 * corner_values[i];
196 [ # # ][ # # ]: 0 : for( r = 0; r < n; ++r )
197 : : {
198 : 0 : R = idx[r];
199 [ # # ][ # # ]: 0 : vertex_grads[R] += v * *grad;
[ # # ][ # # ]
200 [ # # ][ # # ]: 0 : vertex_hessians[R] += 2.0 * outer( *grad );
[ # # ][ # # ]
[ # # ][ # # ]
201 [ # # ][ # # ]: 0 : vertex_hessians[R] += v * *hess;
[ # # ][ # # ]
[ # # ]
202 : 0 : ++grad;
203 [ # # ]: 0 : ++hess;
204 : : }
205 : : }
206 : 0 : return avg;
207 : : }
208 : :
209 : : template < typename HessIter >
210 : 0 : static inline double pmean_corner_diagonals( EntityTopology type, unsigned num_corner, const double corner_values[],
211 : : const Vector3D corner_grads[], HessIter corner_diag_blocks,
212 : : Vector3D vertex_grads[], SymMatrix3D vertex_hessians[], double p )
213 : : {
214 [ # # ][ # # ]: 0 : const unsigned N = TopologyInfo::corners( type );
215 : : unsigned i, n, r, R, idx[4];
216 : : const unsigned* adj_list;
217 : 0 : double m = 0.0, nm;
218 : : double gf[8], hf[8];
219 : 0 : double inv = 1.0 / num_corner;
220 [ # # ][ # # ]: 0 : assert( num_corner <= 8 );
221 : :
222 : : // calculate mean
223 [ # # ][ # # ]: 0 : for( i = 0; i < num_corner; ++i )
224 : : {
225 : 0 : nm = pow( corner_values[i], p );
226 : 0 : m += nm;
227 : :
228 : 0 : gf[i] = inv * p * nm / corner_values[i];
229 : 0 : hf[i] = ( p - 1 ) * gf[i] / corner_values[i];
230 : : }
231 : 0 : nm = inv * m;
232 : :
233 : 0 : const Vector3D* grad = corner_grads;
234 : 0 : HessIter hess = corner_diag_blocks;
235 [ # # ][ # # ]: 0 : for( i = 0; i < num_corner; ++i )
236 : : {
237 [ # # ][ # # ]: 0 : adj_list = TopologyInfo::adjacent_vertices( type, i, n );
238 : 0 : idx[0] = i;
239 : 0 : idx[1] = adj_list[0];
240 : 0 : idx[2] = adj_list[1];
241 : 0 : idx[3] = adj_list[2 % n]; // %n so don't read off end if 2D
242 : 0 : ++n;
243 : :
244 [ # # ][ # # ]: 0 : for( r = 0; r < n; ++r )
245 : : {
246 : 0 : R = idx[r];
247 [ # # ][ # # ]: 0 : vertex_grads[R] += gf[i] * *grad;
[ # # ][ # # ]
248 [ # # ][ # # ]: 0 : vertex_hessians[R] += hf[i] * outer( *grad );
[ # # ][ # # ]
[ # # ][ # # ]
249 [ # # ][ # # ]: 0 : vertex_hessians[R] += gf[i] * *hess;
[ # # ][ # # ]
[ # # ]
250 : 0 : ++grad;
251 [ # # ]: 0 : ++hess;
252 : : }
253 : : }
254 : :
255 : 0 : m = pow( nm, 1.0 / p );
256 : 0 : gf[0] = m / ( p * nm );
257 : 0 : hf[0] = ( 1.0 / p - 1 ) * gf[0] / nm;
258 [ # # ][ # # ]: 0 : for( r = 0; r < N; ++r )
259 : : {
260 [ # # ][ # # ]: 0 : vertex_hessians[r] *= gf[0];
261 [ # # ][ # # ]: 0 : vertex_hessians[r] += hf[0] * outer( vertex_grads[r] );
[ # # ][ # # ]
[ # # ][ # # ]
262 [ # # ][ # # ]: 0 : vertex_grads[r] *= gf[0];
263 : : }
264 : :
265 : 0 : return m;
266 : : }
267 : :
268 : : template < typename HessIter >
269 : 0 : static inline double average_corner_diagonals( EntityTopology type, QualityMetric::AveragingMethod method,
270 : : unsigned num_corner, const double corner_values[],
271 : : const Vector3D corner_grads[], HessIter corner_diag_blocks,
272 : : Vector3D vertex_grads[], SymMatrix3D vertex_hessians[], MsqError& err )
273 : : {
274 : : unsigned i;
275 : : double avg, inv;
276 : :
277 : : // Zero gradients and Hessians
278 : 0 : const unsigned num_vertex = TopologyInfo::corners( type );
279 [ # # ][ # # ]: 0 : for( i = 0; i < num_vertex; ++i )
280 : : {
281 [ # # ][ # # ]: 0 : vertex_grads[i].set( 0.0 );
282 : 0 : vertex_hessians[i] = SymMatrix3D( 0.0 );
283 : : }
284 : :
285 [ # # # # : 0 : switch( method )
# # # ][ #
# # # # #
# ]
286 : : {
287 : : case QualityMetric::SUM:
288 : 0 : avg = sum_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks, vertex_grads,
289 : : vertex_hessians );
290 : 0 : break;
291 : :
292 : : case QualityMetric::LINEAR:
293 : 0 : avg = sum_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks, vertex_grads,
294 : : vertex_hessians );
295 : 0 : inv = 1.0 / num_corner;
296 : 0 : avg *= inv;
297 [ # # ][ # # ]: 0 : for( i = 0; i < num_vertex; ++i )
298 : : {
299 : 0 : vertex_grads[i] *= inv;
300 : 0 : vertex_hessians[i] *= inv;
301 : : }
302 : 0 : break;
303 : :
304 : : case QualityMetric::SUM_SQUARED:
305 : 0 : avg = sum_sqr_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks,
306 : : vertex_grads, vertex_hessians );
307 : 0 : break;
308 : :
309 : : case QualityMetric::RMS:
310 : 0 : avg = pmean_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks,
311 : : vertex_grads, vertex_hessians, 2.0 );
312 : 0 : break;
313 : :
314 : : case QualityMetric::HARMONIC:
315 : 0 : avg = pmean_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks,
316 : : vertex_grads, vertex_hessians, -1.0 );
317 : 0 : break;
318 : :
319 : : case QualityMetric::HMS:
320 : 0 : avg = pmean_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks,
321 : : vertex_grads, vertex_hessians, -2.0 );
322 : 0 : break;
323 : :
324 : : default:
325 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "averaging method not available.", MsqError::INVALID_STATE );
326 : 0 : return 0.0;
327 : : }
328 : :
329 : 0 : return avg;
330 : : }
331 : :
332 : 0 : double AveragingQM::average_corner_hessian_diagonals( EntityTopology element_type, uint32_t, unsigned num_corners,
333 : : const double corner_values[], const Vector3D corner_grads[],
334 : : const Matrix3D corner_hessians[], Vector3D vertex_grads[],
335 : : SymMatrix3D vertex_hessians[], MsqError& err )
336 : : {
337 : : return average_corner_diagonals( element_type, avgMethod, num_corners, corner_values, corner_grads,
338 : : CornerHessDiagIterator( corner_hessians, element_type ), vertex_grads,
339 [ # # ]: 0 : vertex_hessians, err );
340 : : }
341 : :
342 : 0 : double AveragingQM::average_corner_hessian_diagonals( EntityTopology element_type, uint32_t, unsigned num_corners,
343 : : const double corner_values[], const Vector3D corner_grads[],
344 : : const SymMatrix3D corner_hess_diag[], Vector3D vertex_grads[],
345 : : SymMatrix3D vertex_hessians[], MsqError& err )
346 : : {
347 : : return average_corner_diagonals( element_type, avgMethod, num_corners, corner_values, corner_grads,
348 : 0 : corner_hess_diag, vertex_grads, vertex_hessians, err );
349 : : }
350 : :
351 : 14793 : static inline double sum_corner_hessians( EntityTopology type, unsigned num_corner, const double corner_values[],
352 : : const Vector3D corner_grads[], const Matrix3D corner_hessians[],
353 : : Vector3D vertex_grads[], Matrix3D vertex_hessians[] )
354 : : {
355 [ + - ]: 14793 : const unsigned N = TopologyInfo::corners( type );
356 : : unsigned i, n, r, c, R, C, idx[4];
357 : : const unsigned* adj_list;
358 : 14793 : double avg = 0.0;
359 : :
360 : : // calculate mean
361 [ + + ]: 123649 : for( i = 0; i < num_corner; ++i )
362 : 108856 : avg += corner_values[i];
363 : :
364 : 14793 : const Vector3D* grad = corner_grads;
365 : 14793 : const Matrix3D* hess = corner_hessians;
366 [ + + ]: 123649 : for( i = 0; i < num_corner; ++i )
367 : : {
368 [ + - ]: 108856 : adj_list = TopologyInfo::adjacent_vertices( type, i, n );
369 : 108856 : idx[0] = i;
370 : 108856 : idx[1] = adj_list[0];
371 : 108856 : idx[2] = adj_list[1];
372 : 108856 : idx[3] = adj_list[2 % n]; // %n so don't read off end if 2D
373 : :
374 [ + + ]: 544280 : for( r = 0; r <= n; ++r )
375 : : {
376 : 435424 : R = idx[r];
377 [ + - ]: 435424 : vertex_grads[R] += *grad;
378 : 435424 : ++grad;
379 [ + + ]: 1523984 : for( c = r; c <= n; ++c )
380 : : {
381 : 1088560 : C = idx[c];
382 [ + + ]: 1088560 : if( R <= C )
383 [ + - ]: 775828 : vertex_hessians[N * R - R * ( R + 1 ) / 2 + C] += *hess;
384 : : else
385 [ + - ]: 312732 : vertex_hessians[N * C - C * ( C + 1 ) / 2 + R].plus_transpose_equal( *hess );
386 : 1088560 : ++hess;
387 : : }
388 : : }
389 : : }
390 : 14793 : return avg;
391 : : }
392 : :
393 : 0 : static inline double sum_sqr_corner_hessians( EntityTopology type, unsigned num_corner, const double corner_values[],
394 : : const Vector3D corner_grads[], const Matrix3D corner_hessians[],
395 : : Vector3D vertex_grads[], Matrix3D vertex_hessians[] )
396 : : {
397 [ # # ]: 0 : const unsigned N = TopologyInfo::corners( type );
398 : : unsigned i, n, r, c, R, C, idx[4];
399 : : const unsigned* adj_list;
400 : 0 : double v, avg = 0.0;
401 [ # # ]: 0 : Matrix3D op;
402 : :
403 : : // calculate mean
404 [ # # ]: 0 : for( i = 0; i < num_corner; ++i )
405 : 0 : avg += corner_values[i] * corner_values[i];
406 : :
407 : 0 : const Vector3D* grad = corner_grads;
408 : 0 : const Matrix3D* hess = corner_hessians;
409 [ # # ]: 0 : for( i = 0; i < num_corner; ++i )
410 : : {
411 [ # # ]: 0 : adj_list = TopologyInfo::adjacent_vertices( type, i, n );
412 : 0 : idx[0] = i;
413 : 0 : idx[1] = adj_list[0];
414 : 0 : idx[2] = adj_list[1];
415 : 0 : idx[3] = adj_list[2 % n]; // %n so don't read off end if 2D
416 : 0 : ++n;
417 : :
418 : 0 : v = 2.0 * corner_values[i];
419 [ # # ]: 0 : for( r = 0; r < n; ++r )
420 : : {
421 : 0 : R = idx[r];
422 [ # # ][ # # ]: 0 : vertex_grads[R] += v * grad[r];
423 [ # # ]: 0 : for( c = r; c < n; ++c )
424 : : {
425 : 0 : C = idx[c];
426 [ # # ][ # # ]: 0 : op.outer_product( 2.0 * grad[r], grad[c] );
427 [ # # ][ # # ]: 0 : op += v * *hess;
428 [ # # ]: 0 : if( R <= C )
429 [ # # ]: 0 : vertex_hessians[N * R - R * ( R + 1 ) / 2 + C] += op;
430 : : else
431 [ # # ]: 0 : vertex_hessians[N * C - C * ( C + 1 ) / 2 + R].plus_transpose_equal( op );
432 : 0 : ++hess;
433 : : }
434 : : }
435 : 0 : grad += n;
436 : : }
437 : 0 : return avg;
438 : : }
439 : :
440 : 0 : static inline double pmean_corner_hessians( EntityTopology type, unsigned num_corner, const double corner_values[],
441 : : const Vector3D corner_grads[], const Matrix3D corner_hessians[],
442 : : Vector3D vertex_grads[], Matrix3D vertex_hessians[], double p )
443 : : {
444 [ # # ]: 0 : const unsigned N = TopologyInfo::corners( type );
445 : : unsigned i, n, r, c, R, C, idx[4];
446 : : const unsigned* adj_list;
447 : 0 : double m = 0.0, nm;
448 [ # # ]: 0 : Matrix3D op;
449 : : double gf[8], hf[8];
450 : 0 : double inv = 1.0 / num_corner;
451 [ # # ]: 0 : assert( num_corner <= 8 );
452 : :
453 : : // calculate mean
454 [ # # ]: 0 : for( i = 0; i < num_corner; ++i )
455 : : {
456 : 0 : nm = pow( corner_values[i], p );
457 : 0 : m += nm;
458 : :
459 : 0 : gf[i] = inv * p * nm / corner_values[i];
460 : 0 : hf[i] = ( p - 1 ) * gf[i] / corner_values[i];
461 : : }
462 : 0 : nm = inv * m;
463 : :
464 : 0 : const Vector3D* grad = corner_grads;
465 : 0 : const Matrix3D* hess = corner_hessians;
466 [ # # ]: 0 : for( i = 0; i < num_corner; ++i )
467 : : {
468 [ # # ]: 0 : adj_list = TopologyInfo::adjacent_vertices( type, i, n );
469 : 0 : idx[0] = i;
470 : 0 : idx[1] = adj_list[0];
471 : 0 : idx[2] = adj_list[1];
472 : 0 : idx[3] = adj_list[2 % n]; // %n so don't read off end if 2D
473 : 0 : ++n;
474 : :
475 [ # # ]: 0 : for( r = 0; r < n; ++r )
476 : : {
477 : 0 : R = idx[r];
478 [ # # ][ # # ]: 0 : vertex_grads[R] += gf[i] * grad[r];
479 [ # # ]: 0 : for( c = r; c < n; ++c )
480 : : {
481 : 0 : C = idx[c];
482 [ # # ]: 0 : op.outer_product( grad[r], grad[c] );
483 [ # # ]: 0 : op *= hf[i];
484 [ # # ][ # # ]: 0 : op += gf[i] * *hess;
485 [ # # ]: 0 : if( R <= C )
486 [ # # ]: 0 : vertex_hessians[N * R - R * ( R + 1 ) / 2 + C] += op;
487 : : else
488 [ # # ]: 0 : vertex_hessians[N * C - C * ( C + 1 ) / 2 + R].plus_transpose_equal( op );
489 : 0 : ++hess;
490 : : }
491 : : }
492 : 0 : grad += n;
493 : : }
494 : :
495 : 0 : m = pow( nm, 1.0 / p );
496 : 0 : gf[0] = m / ( p * nm );
497 : 0 : hf[0] = ( 1.0 / p - 1 ) * gf[0] / nm;
498 [ # # ]: 0 : for( r = 0; r < N; ++r )
499 : : {
500 [ # # ]: 0 : for( c = r; c < N; ++c )
501 : : {
502 [ # # ]: 0 : op.outer_product( vertex_grads[r], vertex_grads[c] );
503 [ # # ]: 0 : op *= hf[0];
504 [ # # ]: 0 : *vertex_hessians *= gf[0];
505 [ # # ]: 0 : *vertex_hessians += op;
506 : 0 : ++vertex_hessians;
507 : : }
508 [ # # ]: 0 : vertex_grads[r] *= gf[0];
509 : : }
510 : :
511 : 0 : return m;
512 : : }
513 : :
514 : 14793 : double AveragingQM::average_corner_hessians( EntityTopology type, uint32_t, unsigned num_corner,
515 : : const double corner_values[], const Vector3D corner_grads[],
516 : : const Matrix3D corner_hessians[], Vector3D vertex_grads[],
517 : : Matrix3D vertex_hessians[], MsqError& err )
518 : : {
519 : : unsigned i;
520 : : double avg, inv;
521 : :
522 : : // Zero gradients and Hessians
523 : 14793 : const unsigned num_vertex = TopologyInfo::corners( type );
524 [ + + ]: 125955 : for( i = 0; i < num_vertex; ++i )
525 [ + - ]: 111162 : vertex_grads[i].set( 0.0 );
526 : 14793 : const unsigned num_hess = num_vertex * ( num_vertex + 1 ) / 2;
527 [ + + ]: 496935 : for( i = 0; i < num_hess; ++i )
528 : 482142 : vertex_hessians[i].zero();
529 : :
530 [ + + - - : 14793 : switch( avgMethod )
- - - ]
531 : : {
532 : : case QualityMetric::SUM:
533 : : avg = sum_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads,
534 : 9000 : vertex_hessians );
535 : 9000 : break;
536 : :
537 : : case QualityMetric::LINEAR:
538 : : avg = sum_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads,
539 : 5793 : vertex_hessians );
540 : 5793 : inv = 1.0 / num_corner;
541 : 5793 : avg *= inv;
542 [ + + ]: 44955 : for( i = 0; i < num_vertex; ++i )
543 : 39162 : vertex_grads[i] *= inv;
544 [ + + ]: 163935 : for( i = 0; i < num_hess; ++i )
545 : 158142 : vertex_hessians[i] *= inv;
546 : 5793 : break;
547 : :
548 : : case QualityMetric::SUM_SQUARED:
549 : : avg = sum_sqr_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads,
550 : 0 : vertex_hessians );
551 : 0 : break;
552 : :
553 : : case QualityMetric::RMS:
554 : : avg = pmean_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads,
555 : 0 : vertex_hessians, 2.0 );
556 : 0 : break;
557 : :
558 : : case QualityMetric::HARMONIC:
559 : : avg = pmean_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads,
560 : 0 : vertex_hessians, -1.0 );
561 : 0 : break;
562 : :
563 : : case QualityMetric::HMS:
564 : : avg = pmean_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads,
565 : 0 : vertex_hessians, -2.0 );
566 : 0 : break;
567 : :
568 : : default:
569 [ # # ]: 0 : MSQ_SETERR( err )( "averaging method not available.", MsqError::INVALID_STATE );
570 : 0 : return 0.0;
571 : : }
572 : :
573 : 14793 : return avg;
574 : : }
575 : :
576 : 16620 : double AveragingQM::average_metric_and_weights( double metrics[], int count, MsqError& err )
577 : : {
578 : : static bool min_max_warning = false;
579 : 16620 : double avg = 0.0;
580 : : int i, tmp_count;
581 : : double f;
582 : :
583 [ - - + - : 16620 : switch( avgMethod )
+ - - - -
- ]
584 : : {
585 : :
586 : : case QualityMetric::MINIMUM:
587 [ # # ]: 0 : if( !min_max_warning )
588 : : {
589 : : MSQ_DBGOUT( 1 ) << "Minimum and maximum not continuously differentiable.\n"
590 : : "Element of subdifferential will be returned.\n";
591 : 0 : min_max_warning = true;
592 : : }
593 : :
594 : 0 : avg = metrics[0];
595 [ # # ]: 0 : for( i = 1; i < count; ++i )
596 [ # # ]: 0 : if( metrics[i] < avg ) avg = metrics[i];
597 : :
598 : 0 : tmp_count = 0;
599 [ # # ]: 0 : for( i = 0; i < count; ++i )
600 : : {
601 [ # # ]: 0 : if( metrics[i] - avg <= MSQ_MIN )
602 : : {
603 : 0 : metrics[i] = 1.0;
604 : 0 : ++tmp_count;
605 : : }
606 : : else
607 : : {
608 : 0 : metrics[i] = 0.0;
609 : : }
610 : : }
611 : :
612 : 0 : f = 1.0 / tmp_count;
613 [ # # ]: 0 : for( i = 0; i < count; ++i )
614 : 0 : metrics[i] *= f;
615 : :
616 : 0 : break;
617 : :
618 : : case QualityMetric::MAXIMUM:
619 [ # # ]: 0 : if( !min_max_warning )
620 : : {
621 : : MSQ_DBGOUT( 1 ) << "Minimum and maximum not continuously differentiable.\n"
622 : : "Element of subdifferential will be returned.\n";
623 : 0 : min_max_warning = true;
624 : : }
625 : :
626 : 0 : avg = metrics[0];
627 [ # # ]: 0 : for( i = 1; i < count; ++i )
628 [ # # ]: 0 : if( metrics[i] > avg ) avg = metrics[i];
629 : :
630 : 0 : tmp_count = 0;
631 [ # # ]: 0 : for( i = 0; i < count; ++i )
632 : : {
633 [ # # ]: 0 : if( avg - metrics[i] <= MSQ_MIN )
634 : : {
635 : 0 : metrics[i] = 1.0;
636 : 0 : ++tmp_count;
637 : : }
638 : : else
639 : : {
640 : 0 : metrics[i] = 0.0;
641 : : }
642 : : }
643 : :
644 : 0 : f = 1.0 / tmp_count;
645 [ # # ]: 0 : for( i = 0; i < count; ++i )
646 : 0 : metrics[i] *= f;
647 : :
648 : 0 : break;
649 : :
650 : : case QualityMetric::SUM:
651 [ + + ]: 90000 : for( i = 0; i < count; ++i )
652 : : {
653 : 80000 : avg += metrics[i];
654 : 80000 : metrics[i] = 1.0;
655 : : }
656 : :
657 : 10000 : break;
658 : :
659 : : case QualityMetric::SUM_SQUARED:
660 [ # # ]: 0 : for( i = 0; i < count; ++i )
661 : : {
662 : 0 : avg += ( metrics[i] * metrics[i] );
663 : 0 : metrics[i] *= 2;
664 : : }
665 : :
666 : 0 : break;
667 : :
668 : : case QualityMetric::LINEAR:
669 : 6620 : f = 1.0 / count;
670 [ + + ]: 47694 : for( i = 0; i < count; ++i )
671 : : {
672 : 41074 : avg += metrics[i];
673 : 41074 : metrics[i] = f;
674 : : }
675 : 6620 : avg *= f;
676 : :
677 : 6620 : break;
678 : :
679 : : case QualityMetric::GEOMETRIC:
680 : 0 : avg = 1.0;
681 [ # # ]: 0 : for( i = 0; i < count; ++i )
682 : 0 : avg *= metrics[i];
683 : 0 : avg = pow( avg, 1.0 / count );
684 : :
685 : 0 : f = avg / count;
686 [ # # ]: 0 : for( i = 0; i < count; ++i )
687 : 0 : metrics[i] = f / metrics[i];
688 : :
689 : 0 : break;
690 : :
691 : : case QualityMetric::RMS:
692 [ # # ]: 0 : for( i = 0; i < count; ++i )
693 : 0 : avg += metrics[i] * metrics[i];
694 : 0 : avg = sqrt( avg / count );
695 : :
696 : 0 : f = 1. / ( avg * count );
697 [ # # ]: 0 : for( i = 0; i < count; ++i )
698 : 0 : metrics[i] *= f;
699 : :
700 : 0 : break;
701 : :
702 : : case QualityMetric::HARMONIC:
703 [ # # ]: 0 : for( i = 0; i < count; ++i )
704 : 0 : avg += 1.0 / metrics[i];
705 : 0 : avg = count / avg;
706 : :
707 [ # # ]: 0 : for( i = 0; i < count; ++i )
708 : 0 : metrics[i] = ( avg * avg ) / ( count * metrics[i] * metrics[i] );
709 : :
710 : 0 : break;
711 : :
712 : : case QualityMetric::HMS:
713 [ # # ]: 0 : for( i = 0; i < count; ++i )
714 : 0 : avg += 1. / ( metrics[i] * metrics[i] );
715 : 0 : avg = sqrt( count / avg );
716 : :
717 : 0 : f = avg * avg * avg / count;
718 [ # # ]: 0 : for( i = 0; i < count; ++i )
719 : 0 : metrics[i] = f / ( metrics[i] * metrics[i] * metrics[i] );
720 : :
721 : 0 : break;
722 : :
723 : : default:
724 [ # # ]: 0 : MSQ_SETERR( err )( "averaging method not available.", MsqError::INVALID_STATE );
725 : : }
726 : :
727 : 16620 : return avg;
728 : : }
729 : :
730 : : /*!
731 : : average_metrics takes an array of length num_value and averages the
732 : : contents using averaging method 'method'.
733 : : */
734 : 1164539 : double AveragingQM::average_metrics( const double metric_values[], int num_values, MsqError& err )
735 : : {
736 : : // MSQ_MAX needs to be made global?
737 : : // double MSQ_MAX=1e10;
738 : 1164539 : double total_value = 0.0;
739 : 1164539 : double temp_value = 0.0;
740 : 1164539 : int i = 0;
741 : 1164539 : int j = 0;
742 : : // if no values, return zero
743 [ - + ]: 1164539 : if( num_values <= 0 ) { return 0.0; }
744 : :
745 [ - - + - : 1164539 : switch( avgMethod )
- + - - +
- - - -
- ]
746 : : {
747 : : case QualityMetric::GEOMETRIC:
748 : 0 : total_value = 1.0;
749 [ # # ]: 0 : for( i = 0; i < num_values; ++i )
750 : : {
751 : 0 : total_value *= ( metric_values[i] );
752 : : }
753 : 0 : total_value = pow( total_value, 1.0 / num_values );
754 : 0 : break;
755 : :
756 : : case QualityMetric::HARMONIC:
757 : : // ensure no divide by zero, return zero
758 [ # # ]: 0 : for( i = 0; i < num_values; ++i )
759 : : {
760 [ # # ]: 0 : if( metric_values[i] < MSQ_MIN )
761 : : {
762 [ # # ]: 0 : if( metric_values[i] > MSQ_MIN ) { return 0.0; }
763 : : }
764 : 0 : total_value += ( 1 / metric_values[i] );
765 : : }
766 : : // ensure no divide by zero, return MSQ_MAX_CAP
767 [ # # ]: 0 : if( total_value < MSQ_MIN )
768 : : {
769 [ # # ]: 0 : if( total_value > MSQ_MIN ) { return MSQ_MAX_CAP; }
770 : : }
771 : 0 : total_value = num_values / total_value;
772 : 0 : break;
773 : :
774 : : case QualityMetric::LINEAR:
775 [ + + ]: 159574 : for( i = 0; i < num_values; ++i )
776 : : {
777 : 130956 : total_value += metric_values[i];
778 : : }
779 : 28618 : total_value /= (double)num_values;
780 : 28618 : break;
781 : :
782 : : case QualityMetric::MAXIMUM:
783 : 0 : total_value = metric_values[0];
784 [ # # ]: 0 : for( i = 1; i < num_values; ++i )
785 : : {
786 [ # # ]: 0 : if( metric_values[i] > total_value ) { total_value = metric_values[i]; }
787 : : }
788 : 0 : break;
789 : :
790 : : case QualityMetric::MINIMUM:
791 : 0 : total_value = metric_values[0];
792 [ # # ]: 0 : for( i = 1; i < num_values; ++i )
793 : : {
794 [ # # ]: 0 : if( metric_values[i] < total_value ) { total_value = metric_values[i]; }
795 : : }
796 : 0 : break;
797 : :
798 : : case QualityMetric::RMS:
799 [ + + ]: 8472945 : for( i = 0; i < num_values; ++i )
800 : : {
801 : 7343024 : total_value += ( metric_values[i] * metric_values[i] );
802 : : }
803 : 1129921 : total_value /= (double)num_values;
804 : 1129921 : total_value = sqrt( total_value );
805 : 1129921 : break;
806 : :
807 : : case QualityMetric::HMS:
808 : : // ensure no divide by zero, return zero
809 [ # # ]: 0 : for( i = 0; i < num_values; ++i )
810 : : {
811 [ # # ]: 0 : if( metric_values[i] * metric_values[i] < MSQ_MIN ) { return 0.0; }
812 : 0 : total_value += ( 1.0 / ( metric_values[i] * metric_values[i] ) );
813 : : }
814 : :
815 : : // ensure no divide by zero, return MSQ_MAX_CAP
816 [ # # ]: 0 : if( total_value < MSQ_MIN ) { return MSQ_MAX_CAP; }
817 : 0 : total_value = sqrt( num_values / total_value );
818 : 0 : break;
819 : :
820 : : case QualityMetric::STANDARD_DEVIATION:
821 : 0 : total_value = 0;
822 : 0 : temp_value = 0;
823 [ # # ]: 0 : for( i = 0; i < num_values; ++i )
824 : : {
825 : 0 : temp_value += metric_values[i];
826 : 0 : total_value += ( metric_values[i] * metric_values[i] );
827 : : }
828 : 0 : temp_value /= (double)num_values;
829 : 0 : temp_value *= temp_value;
830 : 0 : total_value /= (double)num_values;
831 : 0 : total_value = total_value - temp_value;
832 : 0 : break;
833 : :
834 : : case QualityMetric::SUM:
835 [ + + ]: 54000 : for( i = 0; i < num_values; ++i )
836 : : {
837 : 48000 : total_value += metric_values[i];
838 : : }
839 : 6000 : break;
840 : :
841 : : case QualityMetric::SUM_SQUARED:
842 [ # # ]: 0 : for( i = 0; i < num_values; ++i )
843 : : {
844 : 0 : total_value += ( metric_values[i] * metric_values[i] );
845 : : }
846 : 0 : break;
847 : :
848 : : case QualityMetric::MAX_MINUS_MIN:
849 : : // total_value used to store the maximum
850 : : // temp_value used to store the minimum
851 : 0 : temp_value = MSQ_MAX_CAP;
852 [ # # ]: 0 : for( i = 0; i < num_values; ++i )
853 : : {
854 [ # # ]: 0 : if( metric_values[i] < temp_value ) { temp_value = metric_values[i]; }
855 [ # # ]: 0 : if( metric_values[i] > total_value ) { total_value = metric_values[i]; }
856 : : }
857 : :
858 : 0 : total_value -= temp_value;
859 : 0 : break;
860 : :
861 : : case QualityMetric::MAX_OVER_MIN:
862 : : // total_value used to store the maximum
863 : : // temp_value used to store the minimum
864 : 0 : temp_value = MSQ_MAX_CAP;
865 [ # # ]: 0 : for( i = 0; i < num_values; ++i )
866 : : {
867 [ # # ]: 0 : if( metric_values[i] < temp_value ) { temp_value = metric_values[i]; }
868 [ # # ]: 0 : if( metric_values[i] > total_value ) { total_value = metric_values[i]; }
869 : : }
870 : :
871 : : // ensure no divide by zero, return MSQ_MAX_CAP
872 [ # # ]: 0 : if( fabs( temp_value ) < MSQ_MIN ) { return MSQ_MAX_CAP; }
873 : 0 : total_value /= temp_value;
874 : 0 : break;
875 : :
876 : : case QualityMetric::SUM_OF_RATIOS_SQUARED:
877 [ # # ]: 0 : for( j = 0; j < num_values; ++j )
878 : : {
879 : : // ensure no divide by zero, return MSQ_MAX_CAP
880 [ # # ]: 0 : if( fabs( metric_values[j] ) < MSQ_MIN ) { return MSQ_MAX_CAP; }
881 [ # # ]: 0 : for( i = 0; i < num_values; ++i )
882 : : {
883 : : total_value +=
884 : 0 : ( ( metric_values[i] / metric_values[j] ) * ( metric_values[i] / metric_values[j] ) );
885 : : }
886 : : }
887 : 0 : total_value /= (double)( num_values * num_values );
888 : 0 : break;
889 : :
890 : : default:
891 : : // Return error saying Averaging Method mode not implemented
892 : : MSQ_SETERR( err )
893 [ # # ]: 0 : ( "Requested Averaging Method Not Implemented", MsqError::NOT_IMPLEMENTED );
894 : 0 : return 0;
895 : : }
896 : 1164539 : return total_value;
897 : : }
898 : :
899 [ + - ][ + - ]: 76 : } // namespace MBMesquite
|
