MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2004 Sandia Corporation and Argonne National 00005 Laboratory. Under the terms of Contract DE-AC04-94AL85000 00006 with Sandia Corporation, the U.S. Government retains certain 00007 rights in this software. 00008 00009 This library is free software; you can redistribute it and/or 00010 modify it under the terms of the GNU Lesser General Public 00011 License as published by the Free Software Foundation; either 00012 version 2.1 of the License, or (at your option) any later version. 00013 00014 This library is distributed in the hope that it will be useful, 00015 but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 Lesser General Public License for more details. 00018 00019 You should have received a copy of the GNU Lesser General Public License 00020 (lgpl.txt) along with this library; if not, write to the Free Software 00021 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 00023 diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov, 00024 pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov 00025 00026 (2006) kraftche@cae.wisc.edu 00027 00028 ***************************************************************** */ 00029 /*! 00030 \file AveragingQM.cpp 00031 \brief 00032 00033 \author Michael Brewer 00034 \author Thomas Leurent 00035 \author Jason Kraftcheck 00036 \date 2002-05-14 00037 */ 00038 00039 #include "AveragingQM.hpp" 00040 #include "MsqVertex.hpp" 00041 #include "MsqMeshEntity.hpp" 00042 #include "MsqDebug.hpp" 00043 #include "MsqTimer.hpp" 00044 #include "PatchData.hpp" 00045 00046 namespace MBMesquite 00047 { 00048 00049 double AveragingQM::average_corner_gradients( EntityTopology type, uint32_t fixed_vertices, unsigned num_corner, 00050 double corner_values[], const Vector3D corner_grads[], 00051 Vector3D vertex_grads[], MsqError& err ) 00052 { 00053 const unsigned num_vertex = TopologyInfo::corners( type ); 00054 const unsigned dim = TopologyInfo::dimension( type ); 00055 const unsigned per_vertex = dim + 1; 00056 00057 unsigned i, j, num_adj; 00058 const unsigned *adj_idx, *rev_idx; 00059 00060 // NOTE: This function changes the corner_values array such that 00061 // it contains the gradient coefficients. 00062 double avg = average_metric_and_weights( corner_values, num_corner, err ); 00063 MSQ_ERRZERO( err ); 00064 00065 for( i = 0; i < num_vertex; ++i ) 00066 { 00067 if( fixed_vertices & ( 1 << i ) ) // skip fixed vertices 00068 continue; 00069 00070 adj_idx = TopologyInfo::adjacent_vertices( type, i, num_adj ); 00071 rev_idx = TopologyInfo::reverse_vertex_adjacency_offsets( type, i, num_adj ); 00072 if( i < num_corner ) // not all vertices are corners (e.g. pyramid) 00073 vertex_grads[i] = corner_values[i] * corner_grads[per_vertex * i]; 00074 else 00075 vertex_grads[i] = 0; 00076 for( j = 0; j < num_adj; ++j ) 00077 { 00078 const unsigned v = adj_idx[j], c = rev_idx[j] + 1; 00079 if( v >= num_corner ) // if less corners than vertices (e.g. pyramid apex) 00080 continue; 00081 vertex_grads[i] += corner_values[v] * corner_grads[per_vertex * v + c]; 00082 } 00083 } 00084 00085 return avg; 00086 } 00087 00088 /**\brief Iterate over only diagonal blocks of element corner Hessian data 00089 * 00090 * Given concatenation of corner Hessian data for an element, iterate 00091 * over only the diagonal terms for each corner. This class allows 00092 * common code to be used to generate Hessian diagonal blocks from either 00093 * the diagonal blocks for each corner or the full Hessian data for each 00094 * corner, where this class is used for the latter. 00095 */ 00096 class CornerHessDiagIterator 00097 { 00098 private: 00099 const Matrix3D* cornerHess; //!< Current location in concatenated Hessian data. 00100 const EntityTopology elemType; //!< Element topology for Hessian data 00101 unsigned mCorner; //!< The element corner for which cornerHess 00102 //!< is pointing into the corresponding Hessian data. 00103 unsigned mStep; //!< Amount to step to reach next diagonal block. 00104 public: 00105 CornerHessDiagIterator( const Matrix3D* corner_hessians, EntityTopology elem_type ) 00106 : cornerHess( corner_hessians ), elemType( elem_type ), mCorner( 0 ) 00107 { 00108 TopologyInfo::adjacent_vertices( elemType, mCorner, mStep ); 00109 ++mStep; 00110 } 00111 00112 SymMatrix3D operator*() const 00113 { 00114 return cornerHess->upper(); 00115 } 00116 00117 CornerHessDiagIterator& operator++() 00118 { 00119 cornerHess += mStep; 00120 if( !--mStep ) 00121 { 00122 TopologyInfo::adjacent_vertices( elemType, ++mCorner, mStep ); 00123 ++mStep; 00124 } 00125 return *this; 00126 } 00127 00128 CornerHessDiagIterator operator++( int ) 00129 { 00130 CornerHessDiagIterator copy( *this ); 00131 operator++(); 00132 return copy; 00133 } 00134 }; 00135 00136 template < typename HessIter > 00137 static inline double sum_corner_diagonals( EntityTopology type, unsigned num_corner, const double corner_values[], 00138 const Vector3D corner_grads[], HessIter corner_diag_blocks, 00139 Vector3D vertex_grads[], SymMatrix3D vertex_hessians[] ) 00140 { 00141 unsigned i, n, r, R, idx[4]; 00142 const unsigned* adj_list; 00143 double avg = 0.0; 00144 00145 // calculate mean 00146 for( i = 0; i < num_corner; ++i ) 00147 avg += corner_values[i]; 00148 00149 const Vector3D* grad = corner_grads; 00150 HessIter hess = corner_diag_blocks; 00151 for( i = 0; i < num_corner; ++i ) 00152 { 00153 adj_list = TopologyInfo::adjacent_vertices( type, i, n ); 00154 idx[0] = i; 00155 idx[1] = adj_list[0]; 00156 idx[2] = adj_list[1]; 00157 idx[3] = adj_list[2 % n]; // %n so don't read off end if 2D 00158 00159 for( r = 0; r <= n; ++r ) 00160 { 00161 R = idx[r]; 00162 vertex_grads[R] += *grad; 00163 vertex_hessians[R] += *hess; 00164 ++grad; 00165 ++hess; 00166 } 00167 } 00168 return avg; 00169 } 00170 00171 template < typename HessIter > 00172 static inline double sum_sqr_corner_diagonals( EntityTopology type, unsigned num_corner, const double corner_values[], 00173 const Vector3D corner_grads[], HessIter corner_diag_blocks, 00174 Vector3D vertex_grads[], SymMatrix3D vertex_hessians[] ) 00175 { 00176 unsigned i, n, r, R, idx[4]; 00177 const unsigned* adj_list; 00178 double v, avg = 0.0; 00179 00180 // calculate mean 00181 for( i = 0; i < num_corner; ++i ) 00182 avg += corner_values[i] * corner_values[i]; 00183 00184 const Vector3D* grad = corner_grads; 00185 HessIter hess = corner_diag_blocks; 00186 for( i = 0; i < num_corner; ++i ) 00187 { 00188 adj_list = TopologyInfo::adjacent_vertices( type, i, n ); 00189 idx[0] = i; 00190 idx[1] = adj_list[0]; 00191 idx[2] = adj_list[1]; 00192 idx[3] = adj_list[2 % n]; // %n so don't read off end if 2D 00193 ++n; 00194 00195 v = 2.0 * corner_values[i]; 00196 for( r = 0; r < n; ++r ) 00197 { 00198 R = idx[r]; 00199 vertex_grads[R] += v * *grad; 00200 vertex_hessians[R] += 2.0 * outer( *grad ); 00201 vertex_hessians[R] += v * *hess; 00202 ++grad; 00203 ++hess; 00204 } 00205 } 00206 return avg; 00207 } 00208 00209 template < typename HessIter > 00210 static inline double pmean_corner_diagonals( EntityTopology type, unsigned num_corner, const double corner_values[], 00211 const Vector3D corner_grads[], HessIter corner_diag_blocks, 00212 Vector3D vertex_grads[], SymMatrix3D vertex_hessians[], double p ) 00213 { 00214 const unsigned N = TopologyInfo::corners( type ); 00215 unsigned i, n, r, R, idx[4]; 00216 const unsigned* adj_list; 00217 double m = 0.0, nm; 00218 double gf[8], hf[8]; 00219 double inv = 1.0 / num_corner; 00220 assert( num_corner <= 8 ); 00221 00222 // calculate mean 00223 for( i = 0; i < num_corner; ++i ) 00224 { 00225 nm = pow( corner_values[i], p ); 00226 m += nm; 00227 00228 gf[i] = inv * p * nm / corner_values[i]; 00229 hf[i] = ( p - 1 ) * gf[i] / corner_values[i]; 00230 } 00231 nm = inv * m; 00232 00233 const Vector3D* grad = corner_grads; 00234 HessIter hess = corner_diag_blocks; 00235 for( i = 0; i < num_corner; ++i ) 00236 { 00237 adj_list = TopologyInfo::adjacent_vertices( type, i, n ); 00238 idx[0] = i; 00239 idx[1] = adj_list[0]; 00240 idx[2] = adj_list[1]; 00241 idx[3] = adj_list[2 % n]; // %n so don't read off end if 2D 00242 ++n; 00243 00244 for( r = 0; r < n; ++r ) 00245 { 00246 R = idx[r]; 00247 vertex_grads[R] += gf[i] * *grad; 00248 vertex_hessians[R] += hf[i] * outer( *grad ); 00249 vertex_hessians[R] += gf[i] * *hess; 00250 ++grad; 00251 ++hess; 00252 } 00253 } 00254 00255 m = pow( nm, 1.0 / p ); 00256 gf[0] = m / ( p * nm ); 00257 hf[0] = ( 1.0 / p - 1 ) * gf[0] / nm; 00258 for( r = 0; r < N; ++r ) 00259 { 00260 vertex_hessians[r] *= gf[0]; 00261 vertex_hessians[r] += hf[0] * outer( vertex_grads[r] ); 00262 vertex_grads[r] *= gf[0]; 00263 } 00264 00265 return m; 00266 } 00267 00268 template < typename HessIter > 00269 static inline double average_corner_diagonals( EntityTopology type, QualityMetric::AveragingMethod method, 00270 unsigned num_corner, const double corner_values[], 00271 const Vector3D corner_grads[], HessIter corner_diag_blocks, 00272 Vector3D vertex_grads[], SymMatrix3D vertex_hessians[], MsqError& err ) 00273 { 00274 unsigned i; 00275 double avg, inv; 00276 00277 // Zero gradients and Hessians 00278 const unsigned num_vertex = TopologyInfo::corners( type ); 00279 for( i = 0; i < num_vertex; ++i ) 00280 { 00281 vertex_grads[i].set( 0.0 ); 00282 vertex_hessians[i] = SymMatrix3D( 0.0 ); 00283 } 00284 00285 switch( method ) 00286 { 00287 case QualityMetric::SUM: 00288 avg = sum_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks, vertex_grads, 00289 vertex_hessians ); 00290 break; 00291 00292 case QualityMetric::LINEAR: 00293 avg = sum_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks, vertex_grads, 00294 vertex_hessians ); 00295 inv = 1.0 / num_corner; 00296 avg *= inv; 00297 for( i = 0; i < num_vertex; ++i ) 00298 { 00299 vertex_grads[i] *= inv; 00300 vertex_hessians[i] *= inv; 00301 } 00302 break; 00303 00304 case QualityMetric::SUM_SQUARED: 00305 avg = sum_sqr_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks, 00306 vertex_grads, vertex_hessians ); 00307 break; 00308 00309 case QualityMetric::RMS: 00310 avg = pmean_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks, 00311 vertex_grads, vertex_hessians, 2.0 ); 00312 break; 00313 00314 case QualityMetric::HARMONIC: 00315 avg = pmean_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks, 00316 vertex_grads, vertex_hessians, -1.0 ); 00317 break; 00318 00319 case QualityMetric::HMS: 00320 avg = pmean_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks, 00321 vertex_grads, vertex_hessians, -2.0 ); 00322 break; 00323 00324 default: 00325 MSQ_SETERR( err )( "averaging method not available.", MsqError::INVALID_STATE ); 00326 return 0.0; 00327 } 00328 00329 return avg; 00330 } 00331 00332 double AveragingQM::average_corner_hessian_diagonals( EntityTopology element_type, uint32_t, unsigned num_corners, 00333 const double corner_values[], const Vector3D corner_grads[], 00334 const Matrix3D corner_hessians[], Vector3D vertex_grads[], 00335 SymMatrix3D vertex_hessians[], MsqError& err ) 00336 { 00337 return average_corner_diagonals( element_type, avgMethod, num_corners, corner_values, corner_grads, 00338 CornerHessDiagIterator( corner_hessians, element_type ), vertex_grads, 00339 vertex_hessians, err ); 00340 } 00341 00342 double AveragingQM::average_corner_hessian_diagonals( EntityTopology element_type, uint32_t, unsigned num_corners, 00343 const double corner_values[], const Vector3D corner_grads[], 00344 const SymMatrix3D corner_hess_diag[], Vector3D vertex_grads[], 00345 SymMatrix3D vertex_hessians[], MsqError& err ) 00346 { 00347 return average_corner_diagonals( element_type, avgMethod, num_corners, corner_values, corner_grads, 00348 corner_hess_diag, vertex_grads, vertex_hessians, err ); 00349 } 00350 00351 static inline double sum_corner_hessians( EntityTopology type, unsigned num_corner, const double corner_values[], 00352 const Vector3D corner_grads[], const Matrix3D corner_hessians[], 00353 Vector3D vertex_grads[], Matrix3D vertex_hessians[] ) 00354 { 00355 const unsigned N = TopologyInfo::corners( type ); 00356 unsigned i, n, r, c, R, C, idx[4]; 00357 const unsigned* adj_list; 00358 double avg = 0.0; 00359 00360 // calculate mean 00361 for( i = 0; i < num_corner; ++i ) 00362 avg += corner_values[i]; 00363 00364 const Vector3D* grad = corner_grads; 00365 const Matrix3D* hess = corner_hessians; 00366 for( i = 0; i < num_corner; ++i ) 00367 { 00368 adj_list = TopologyInfo::adjacent_vertices( type, i, n ); 00369 idx[0] = i; 00370 idx[1] = adj_list[0]; 00371 idx[2] = adj_list[1]; 00372 idx[3] = adj_list[2 % n]; // %n so don't read off end if 2D 00373 00374 for( r = 0; r <= n; ++r ) 00375 { 00376 R = idx[r]; 00377 vertex_grads[R] += *grad; 00378 ++grad; 00379 for( c = r; c <= n; ++c ) 00380 { 00381 C = idx[c]; 00382 if( R <= C ) 00383 vertex_hessians[N * R - R * ( R + 1 ) / 2 + C] += *hess; 00384 else 00385 vertex_hessians[N * C - C * ( C + 1 ) / 2 + R].plus_transpose_equal( *hess ); 00386 ++hess; 00387 } 00388 } 00389 } 00390 return avg; 00391 } 00392 00393 static inline double sum_sqr_corner_hessians( EntityTopology type, unsigned num_corner, const double corner_values[], 00394 const Vector3D corner_grads[], const Matrix3D corner_hessians[], 00395 Vector3D vertex_grads[], Matrix3D vertex_hessians[] ) 00396 { 00397 const unsigned N = TopologyInfo::corners( type ); 00398 unsigned i, n, r, c, R, C, idx[4]; 00399 const unsigned* adj_list; 00400 double v, avg = 0.0; 00401 Matrix3D op; 00402 00403 // calculate mean 00404 for( i = 0; i < num_corner; ++i ) 00405 avg += corner_values[i] * corner_values[i]; 00406 00407 const Vector3D* grad = corner_grads; 00408 const Matrix3D* hess = corner_hessians; 00409 for( i = 0; i < num_corner; ++i ) 00410 { 00411 adj_list = TopologyInfo::adjacent_vertices( type, i, n ); 00412 idx[0] = i; 00413 idx[1] = adj_list[0]; 00414 idx[2] = adj_list[1]; 00415 idx[3] = adj_list[2 % n]; // %n so don't read off end if 2D 00416 ++n; 00417 00418 v = 2.0 * corner_values[i]; 00419 for( r = 0; r < n; ++r ) 00420 { 00421 R = idx[r]; 00422 vertex_grads[R] += v * grad[r]; 00423 for( c = r; c < n; ++c ) 00424 { 00425 C = idx[c]; 00426 op.outer_product( 2.0 * grad[r], grad[c] ); 00427 op += v * *hess; 00428 if( R <= C ) 00429 vertex_hessians[N * R - R * ( R + 1 ) / 2 + C] += op; 00430 else 00431 vertex_hessians[N * C - C * ( C + 1 ) / 2 + R].plus_transpose_equal( op ); 00432 ++hess; 00433 } 00434 } 00435 grad += n; 00436 } 00437 return avg; 00438 } 00439 00440 static inline double pmean_corner_hessians( EntityTopology type, unsigned num_corner, const double corner_values[], 00441 const Vector3D corner_grads[], const Matrix3D corner_hessians[], 00442 Vector3D vertex_grads[], Matrix3D vertex_hessians[], double p ) 00443 { 00444 const unsigned N = TopologyInfo::corners( type ); 00445 unsigned i, n, r, c, R, C, idx[4]; 00446 const unsigned* adj_list; 00447 double m = 0.0, nm; 00448 Matrix3D op; 00449 double gf[8], hf[8]; 00450 double inv = 1.0 / num_corner; 00451 assert( num_corner <= 8 ); 00452 00453 // calculate mean 00454 for( i = 0; i < num_corner; ++i ) 00455 { 00456 nm = pow( corner_values[i], p ); 00457 m += nm; 00458 00459 gf[i] = inv * p * nm / corner_values[i]; 00460 hf[i] = ( p - 1 ) * gf[i] / corner_values[i]; 00461 } 00462 nm = inv * m; 00463 00464 const Vector3D* grad = corner_grads; 00465 const Matrix3D* hess = corner_hessians; 00466 for( i = 0; i < num_corner; ++i ) 00467 { 00468 adj_list = TopologyInfo::adjacent_vertices( type, i, n ); 00469 idx[0] = i; 00470 idx[1] = adj_list[0]; 00471 idx[2] = adj_list[1]; 00472 idx[3] = adj_list[2 % n]; // %n so don't read off end if 2D 00473 ++n; 00474 00475 for( r = 0; r < n; ++r ) 00476 { 00477 R = idx[r]; 00478 vertex_grads[R] += gf[i] * grad[r]; 00479 for( c = r; c < n; ++c ) 00480 { 00481 C = idx[c]; 00482 op.outer_product( grad[r], grad[c] ); 00483 op *= hf[i]; 00484 op += gf[i] * *hess; 00485 if( R <= C ) 00486 vertex_hessians[N * R - R * ( R + 1 ) / 2 + C] += op; 00487 else 00488 vertex_hessians[N * C - C * ( C + 1 ) / 2 + R].plus_transpose_equal( op ); 00489 ++hess; 00490 } 00491 } 00492 grad += n; 00493 } 00494 00495 m = pow( nm, 1.0 / p ); 00496 gf[0] = m / ( p * nm ); 00497 hf[0] = ( 1.0 / p - 1 ) * gf[0] / nm; 00498 for( r = 0; r < N; ++r ) 00499 { 00500 for( c = r; c < N; ++c ) 00501 { 00502 op.outer_product( vertex_grads[r], vertex_grads[c] ); 00503 op *= hf[0]; 00504 *vertex_hessians *= gf[0]; 00505 *vertex_hessians += op; 00506 ++vertex_hessians; 00507 } 00508 vertex_grads[r] *= gf[0]; 00509 } 00510 00511 return m; 00512 } 00513 00514 double AveragingQM::average_corner_hessians( EntityTopology type, uint32_t, unsigned num_corner, 00515 const double corner_values[], const Vector3D corner_grads[], 00516 const Matrix3D corner_hessians[], Vector3D vertex_grads[], 00517 Matrix3D vertex_hessians[], MsqError& err ) 00518 { 00519 unsigned i; 00520 double avg, inv; 00521 00522 // Zero gradients and Hessians 00523 const unsigned num_vertex = TopologyInfo::corners( type ); 00524 for( i = 0; i < num_vertex; ++i ) 00525 vertex_grads[i].set( 0.0 ); 00526 const unsigned num_hess = num_vertex * ( num_vertex + 1 ) / 2; 00527 for( i = 0; i < num_hess; ++i ) 00528 vertex_hessians[i].zero(); 00529 00530 switch( avgMethod ) 00531 { 00532 case QualityMetric::SUM: 00533 avg = sum_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads, 00534 vertex_hessians ); 00535 break; 00536 00537 case QualityMetric::LINEAR: 00538 avg = sum_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads, 00539 vertex_hessians ); 00540 inv = 1.0 / num_corner; 00541 avg *= inv; 00542 for( i = 0; i < num_vertex; ++i ) 00543 vertex_grads[i] *= inv; 00544 for( i = 0; i < num_hess; ++i ) 00545 vertex_hessians[i] *= inv; 00546 break; 00547 00548 case QualityMetric::SUM_SQUARED: 00549 avg = sum_sqr_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads, 00550 vertex_hessians ); 00551 break; 00552 00553 case QualityMetric::RMS: 00554 avg = pmean_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads, 00555 vertex_hessians, 2.0 ); 00556 break; 00557 00558 case QualityMetric::HARMONIC: 00559 avg = pmean_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads, 00560 vertex_hessians, -1.0 ); 00561 break; 00562 00563 case QualityMetric::HMS: 00564 avg = pmean_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads, 00565 vertex_hessians, -2.0 ); 00566 break; 00567 00568 default: 00569 MSQ_SETERR( err )( "averaging method not available.", MsqError::INVALID_STATE ); 00570 return 0.0; 00571 } 00572 00573 return avg; 00574 } 00575 00576 double AveragingQM::average_metric_and_weights( double metrics[], int count, MsqError& err ) 00577 { 00578 static bool min_max_warning = false; 00579 double avg = 0.0; 00580 int i, tmp_count; 00581 double f; 00582 00583 switch( avgMethod ) 00584 { 00585 00586 case QualityMetric::MINIMUM: 00587 if( !min_max_warning ) 00588 { 00589 MSQ_DBGOUT( 1 ) << "Minimum and maximum not continuously differentiable.\n" 00590 "Element of subdifferential will be returned.\n"; 00591 min_max_warning = true; 00592 } 00593 00594 avg = metrics[0]; 00595 for( i = 1; i < count; ++i ) 00596 if( metrics[i] < avg ) avg = metrics[i]; 00597 00598 tmp_count = 0; 00599 for( i = 0; i < count; ++i ) 00600 { 00601 if( metrics[i] - avg <= MSQ_MIN ) 00602 { 00603 metrics[i] = 1.0; 00604 ++tmp_count; 00605 } 00606 else 00607 { 00608 metrics[i] = 0.0; 00609 } 00610 } 00611 00612 f = 1.0 / tmp_count; 00613 for( i = 0; i < count; ++i ) 00614 metrics[i] *= f; 00615 00616 break; 00617 00618 case QualityMetric::MAXIMUM: 00619 if( !min_max_warning ) 00620 { 00621 MSQ_DBGOUT( 1 ) << "Minimum and maximum not continuously differentiable.\n" 00622 "Element of subdifferential will be returned.\n"; 00623 min_max_warning = true; 00624 } 00625 00626 avg = metrics[0]; 00627 for( i = 1; i < count; ++i ) 00628 if( metrics[i] > avg ) avg = metrics[i]; 00629 00630 tmp_count = 0; 00631 for( i = 0; i < count; ++i ) 00632 { 00633 if( avg - metrics[i] <= MSQ_MIN ) 00634 { 00635 metrics[i] = 1.0; 00636 ++tmp_count; 00637 } 00638 else 00639 { 00640 metrics[i] = 0.0; 00641 } 00642 } 00643 00644 f = 1.0 / tmp_count; 00645 for( i = 0; i < count; ++i ) 00646 metrics[i] *= f; 00647 00648 break; 00649 00650 case QualityMetric::SUM: 00651 for( i = 0; i < count; ++i ) 00652 { 00653 avg += metrics[i]; 00654 metrics[i] = 1.0; 00655 } 00656 00657 break; 00658 00659 case QualityMetric::SUM_SQUARED: 00660 for( i = 0; i < count; ++i ) 00661 { 00662 avg += ( metrics[i] * metrics[i] ); 00663 metrics[i] *= 2; 00664 } 00665 00666 break; 00667 00668 case QualityMetric::LINEAR: 00669 f = 1.0 / count; 00670 for( i = 0; i < count; ++i ) 00671 { 00672 avg += metrics[i]; 00673 metrics[i] = f; 00674 } 00675 avg *= f; 00676 00677 break; 00678 00679 case QualityMetric::GEOMETRIC: 00680 avg = 1.0; 00681 for( i = 0; i < count; ++i ) 00682 avg *= metrics[i]; 00683 avg = pow( avg, 1.0 / count ); 00684 00685 f = avg / count; 00686 for( i = 0; i < count; ++i ) 00687 metrics[i] = f / metrics[i]; 00688 00689 break; 00690 00691 case QualityMetric::RMS: 00692 for( i = 0; i < count; ++i ) 00693 avg += metrics[i] * metrics[i]; 00694 avg = sqrt( avg / count ); 00695 00696 f = 1. / ( avg * count ); 00697 for( i = 0; i < count; ++i ) 00698 metrics[i] *= f; 00699 00700 break; 00701 00702 case QualityMetric::HARMONIC: 00703 for( i = 0; i < count; ++i ) 00704 avg += 1.0 / metrics[i]; 00705 avg = count / avg; 00706 00707 for( i = 0; i < count; ++i ) 00708 metrics[i] = ( avg * avg ) / ( count * metrics[i] * metrics[i] ); 00709 00710 break; 00711 00712 case QualityMetric::HMS: 00713 for( i = 0; i < count; ++i ) 00714 avg += 1. / ( metrics[i] * metrics[i] ); 00715 avg = sqrt( count / avg ); 00716 00717 f = avg * avg * avg / count; 00718 for( i = 0; i < count; ++i ) 00719 metrics[i] = f / ( metrics[i] * metrics[i] * metrics[i] ); 00720 00721 break; 00722 00723 default: 00724 MSQ_SETERR( err )( "averaging method not available.", MsqError::INVALID_STATE ); 00725 } 00726 00727 return avg; 00728 } 00729 00730 /*! 00731 average_metrics takes an array of length num_value and averages the 00732 contents using averaging method 'method'. 00733 */ 00734 double AveragingQM::average_metrics( const double metric_values[], int num_values, MsqError& err ) 00735 { 00736 // MSQ_MAX needs to be made global? 00737 // double MSQ_MAX=1e10; 00738 double total_value = 0.0; 00739 double temp_value = 0.0; 00740 int i = 0; 00741 int j = 0; 00742 // if no values, return zero 00743 if( num_values <= 0 ) { return 0.0; } 00744 00745 switch( avgMethod ) 00746 { 00747 case QualityMetric::GEOMETRIC: 00748 total_value = 1.0; 00749 for( i = 0; i < num_values; ++i ) 00750 { 00751 total_value *= ( metric_values[i] ); 00752 } 00753 total_value = pow( total_value, 1.0 / num_values ); 00754 break; 00755 00756 case QualityMetric::HARMONIC: 00757 // ensure no divide by zero, return zero 00758 for( i = 0; i < num_values; ++i ) 00759 { 00760 if( metric_values[i] < MSQ_MIN ) 00761 { 00762 if( metric_values[i] > MSQ_MIN ) { return 0.0; } 00763 } 00764 total_value += ( 1 / metric_values[i] ); 00765 } 00766 // ensure no divide by zero, return MSQ_MAX_CAP 00767 if( total_value < MSQ_MIN ) 00768 { 00769 if( total_value > MSQ_MIN ) { return MSQ_MAX_CAP; } 00770 } 00771 total_value = num_values / total_value; 00772 break; 00773 00774 case QualityMetric::LINEAR: 00775 for( i = 0; i < num_values; ++i ) 00776 { 00777 total_value += metric_values[i]; 00778 } 00779 total_value /= (double)num_values; 00780 break; 00781 00782 case QualityMetric::MAXIMUM: 00783 total_value = metric_values[0]; 00784 for( i = 1; i < num_values; ++i ) 00785 { 00786 if( metric_values[i] > total_value ) { total_value = metric_values[i]; } 00787 } 00788 break; 00789 00790 case QualityMetric::MINIMUM: 00791 total_value = metric_values[0]; 00792 for( i = 1; i < num_values; ++i ) 00793 { 00794 if( metric_values[i] < total_value ) { total_value = metric_values[i]; } 00795 } 00796 break; 00797 00798 case QualityMetric::RMS: 00799 for( i = 0; i < num_values; ++i ) 00800 { 00801 total_value += ( metric_values[i] * metric_values[i] ); 00802 } 00803 total_value /= (double)num_values; 00804 total_value = sqrt( total_value ); 00805 break; 00806 00807 case QualityMetric::HMS: 00808 // ensure no divide by zero, return zero 00809 for( i = 0; i < num_values; ++i ) 00810 { 00811 if( metric_values[i] * metric_values[i] < MSQ_MIN ) { return 0.0; } 00812 total_value += ( 1.0 / ( metric_values[i] * metric_values[i] ) ); 00813 } 00814 00815 // ensure no divide by zero, return MSQ_MAX_CAP 00816 if( total_value < MSQ_MIN ) { return MSQ_MAX_CAP; } 00817 total_value = sqrt( num_values / total_value ); 00818 break; 00819 00820 case QualityMetric::STANDARD_DEVIATION: 00821 total_value = 0; 00822 temp_value = 0; 00823 for( i = 0; i < num_values; ++i ) 00824 { 00825 temp_value += metric_values[i]; 00826 total_value += ( metric_values[i] * metric_values[i] ); 00827 } 00828 temp_value /= (double)num_values; 00829 temp_value *= temp_value; 00830 total_value /= (double)num_values; 00831 total_value = total_value - temp_value; 00832 break; 00833 00834 case QualityMetric::SUM: 00835 for( i = 0; i < num_values; ++i ) 00836 { 00837 total_value += metric_values[i]; 00838 } 00839 break; 00840 00841 case QualityMetric::SUM_SQUARED: 00842 for( i = 0; i < num_values; ++i ) 00843 { 00844 total_value += ( metric_values[i] * metric_values[i] ); 00845 } 00846 break; 00847 00848 case QualityMetric::MAX_MINUS_MIN: 00849 // total_value used to store the maximum 00850 // temp_value used to store the minimum 00851 temp_value = MSQ_MAX_CAP; 00852 for( i = 0; i < num_values; ++i ) 00853 { 00854 if( metric_values[i] < temp_value ) { temp_value = metric_values[i]; } 00855 if( metric_values[i] > total_value ) { total_value = metric_values[i]; } 00856 } 00857 00858 total_value -= temp_value; 00859 break; 00860 00861 case QualityMetric::MAX_OVER_MIN: 00862 // total_value used to store the maximum 00863 // temp_value used to store the minimum 00864 temp_value = MSQ_MAX_CAP; 00865 for( i = 0; i < num_values; ++i ) 00866 { 00867 if( metric_values[i] < temp_value ) { temp_value = metric_values[i]; } 00868 if( metric_values[i] > total_value ) { total_value = metric_values[i]; } 00869 } 00870 00871 // ensure no divide by zero, return MSQ_MAX_CAP 00872 if( fabs( temp_value ) < MSQ_MIN ) { return MSQ_MAX_CAP; } 00873 total_value /= temp_value; 00874 break; 00875 00876 case QualityMetric::SUM_OF_RATIOS_SQUARED: 00877 for( j = 0; j < num_values; ++j ) 00878 { 00879 // ensure no divide by zero, return MSQ_MAX_CAP 00880 if( fabs( metric_values[j] ) < MSQ_MIN ) { return MSQ_MAX_CAP; } 00881 for( i = 0; i < num_values; ++i ) 00882 { 00883 total_value += 00884 ( ( metric_values[i] / metric_values[j] ) * ( metric_values[i] / metric_values[j] ) ); 00885 } 00886 } 00887 total_value /= (double)( num_values * num_values ); 00888 break; 00889 00890 default: 00891 // Return error saying Averaging Method mode not implemented 00892 MSQ_SETERR( err ) 00893 ( "Requested Averaging Method Not Implemented", MsqError::NOT_IMPLEMENTED ); 00894 return 0; 00895 } 00896 return total_value; 00897 } 00898 00899 } // namespace MBMesquite