MOAB: Mesh Oriented datABase
(version 5.3.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 ***************************************************************** */ 00027 /*! 00028 \file LPtoPTemplate.cpp 00029 \brief 00030 00031 This Objective Function is evaluated using an L P norm to the pth power. 00032 total=(sum (x_i)^pVal) 00033 \author Michael Brewer 00034 \author Thomas Leurent 00035 \date 2002-01-23 00036 */ 00037 #include <cmath> 00038 #include "LPtoPTemplate.hpp" 00039 #include "MsqFreeVertexIndexIterator.hpp" 00040 #include "MsqTimer.hpp" 00041 #include "MsqHessian.hpp" 00042 #include "MsqDebug.hpp" 00043 #include "QualityMetric.hpp" 00044 00045 using namespace MBMesquite; 00046 00047 LPtoPTemplate::LPtoPTemplate( QualityMetric* qualitymetric, short Pinput, MsqError& err ) 00048 : ObjectiveFunctionTemplate( qualitymetric ) 00049 { 00050 pVal = Pinput; 00051 if( pVal < 1 ) 00052 { 00053 MSQ_SETERR( err )( "P_VALUE must be greater than 0.", MsqError::INVALID_ARG ); 00054 return; 00055 } 00056 00057 dividingByN = false; 00058 00059 clear(); 00060 } 00061 00062 LPtoPTemplate::LPtoPTemplate( short P, QualityMetric* qm ) 00063 : ObjectiveFunctionTemplate( qm ), pVal( P ), dividingByN( false ) 00064 { 00065 clear(); 00066 } 00067 00068 void LPtoPTemplate::clear() 00069 { 00070 mCount = 0; 00071 mPowSum = 0; 00072 saveCount = 0; 00073 savePowSum = 0; 00074 } 00075 00076 // Michael: need to clean up here 00077 LPtoPTemplate::~LPtoPTemplate() {} 00078 00079 ObjectiveFunction* LPtoPTemplate::clone() const 00080 { 00081 return new LPtoPTemplate( *this ); 00082 } 00083 00084 double LPtoPTemplate::get_value( double power_sum, size_t count, EvalType type, size_t& global_count, 00085 MsqError& /*err*/ ) 00086 { 00087 double result = 0; 00088 switch( type ) 00089 { 00090 default: 00091 case CALCULATE: 00092 result = power_sum; 00093 global_count = count; 00094 break; 00095 00096 case ACCUMULATE: 00097 mPowSum += power_sum; 00098 mCount += count; 00099 result = mPowSum; 00100 global_count = mCount; 00101 break; 00102 00103 case SAVE: 00104 savePowSum = power_sum; 00105 saveCount = count; 00106 result = mPowSum; 00107 global_count = mCount; 00108 break; 00109 00110 case UPDATE: 00111 mPowSum -= savePowSum; 00112 mCount -= saveCount; 00113 savePowSum = power_sum; 00114 saveCount = count; 00115 mPowSum += savePowSum; 00116 mCount += saveCount; 00117 result = mPowSum; 00118 global_count = mCount; 00119 break; 00120 00121 case TEMPORARY: 00122 result = mPowSum - savePowSum + power_sum; 00123 global_count = mCount + count - saveCount; 00124 break; 00125 } 00126 00127 // if (!global_count) 00128 // { 00129 // MSQ_SETERR(err)(" global_count is zero, possibly due to an invalid mesh.", 00130 // MsqError::INVALID_MESH); return -1; // result is invalid 00131 // } 00132 if( dividingByN && global_count ) result /= global_count; 00133 return result; 00134 } 00135 00136 bool LPtoPTemplate::evaluate( EvalType type, PatchData& pd, double& value_out, bool free, MsqError& err ) 00137 { 00138 QualityMetric* qm = get_quality_metric(); 00139 if( type == ObjectiveFunction::ACCUMULATE ) 00140 qm->get_single_pass( pd, qmHandles, free, err ); 00141 else 00142 qm->get_evaluations( pd, qmHandles, free, err ); 00143 MSQ_ERRFALSE( err ); 00144 00145 // calculate OF value for just the patch 00146 std::vector< size_t >::const_iterator i; 00147 double value, working_sum = 0.0; 00148 for( i = qmHandles.begin(); i != qmHandles.end(); ++i ) 00149 { 00150 bool result = qm->evaluate( pd, *i, value, err ); 00151 if( MSQ_CHKERR( err ) || !result ) return false; 00152 00153 double tmp_val = value; 00154 for( short j = 1; j < pVal; ++j ) 00155 tmp_val *= value; 00156 working_sum += fabs( tmp_val ); 00157 } 00158 00159 // get overall OF value, update member data, etc. 00160 size_t global_count; 00161 value_out = qm->get_negate_flag() * get_value( working_sum, qmHandles.size(), type, global_count, err ); 00162 // if (!global_count) 00163 // return false; // invalid mesh 00164 // else 00165 return true; 00166 } 00167 00168 bool LPtoPTemplate::evaluate_with_gradient( EvalType type, PatchData& pd, double& OF_val, 00169 std::vector< Vector3D >& grad_out, MsqError& err ) 00170 { 00171 QualityMetric* qm = get_quality_metric(); 00172 qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err ); 00173 MSQ_ERRFALSE( err ); 00174 00175 // zero gradient 00176 grad_out.clear(); 00177 grad_out.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) ); 00178 bool qm_bool = true; 00179 double QM_val; 00180 OF_val = 0.; 00181 int p1; 00182 00183 // calculate OF value and gradient for just the patch 00184 std::vector< size_t >::const_iterator i; 00185 for( i = qmHandles.begin(); i != qmHandles.end(); ++i ) 00186 { 00187 qm_bool = qm->evaluate_with_gradient( pd, *i, QM_val, mIndices, mGradient, err ); 00188 if( MSQ_CHKERR( err ) || !qm_bool ) return false; 00189 00190 QM_val = fabs( QM_val ); 00191 double QM_pow = 1.0; 00192 double factor = qm->get_negate_flag(); 00193 if( pVal == 1 ) 00194 QM_pow = 1.0; 00195 else 00196 { 00197 QM_pow = QM_val; 00198 for( p1 = 2; p1 < pVal; ++p1 ) 00199 QM_pow *= QM_val; 00200 factor *= QM_pow * pVal; 00201 } 00202 00203 OF_val += QM_pow * QM_val; 00204 for( size_t j = 0; j < mIndices.size(); ++j ) 00205 { 00206 mGradient[j] *= factor; 00207 grad_out[mIndices[j]] += mGradient[j]; 00208 } 00209 } 00210 00211 // get overall OF value, update member data, etc. 00212 size_t global_count; 00213 OF_val = qm->get_negate_flag() * get_value( OF_val, qmHandles.size(), type, global_count, err ); 00214 // if (!global_count) 00215 // return false; // invalid mesh 00216 00217 if( dividingByN && global_count ) 00218 { 00219 const double inv_n = 1.0 / global_count; 00220 std::vector< Vector3D >::iterator g; 00221 for( g = grad_out.begin(); g != grad_out.end(); ++g ) 00222 *g *= inv_n; 00223 } 00224 00225 return true; 00226 } 00227 00228 bool LPtoPTemplate::evaluate_with_Hessian_diagonal( EvalType type, PatchData& pd, double& OF_val, 00229 std::vector< Vector3D >& grad, 00230 std::vector< SymMatrix3D >& hess_diag, MsqError& err ) 00231 { 00232 QualityMetric* qm = get_quality_metric(); 00233 qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err ); 00234 MSQ_ERRFALSE( err ); 00235 00236 // zero gradient and hessian 00237 grad.clear(); 00238 grad.resize( pd.num_free_vertices(), 0.0 ); 00239 hess_diag.clear(); 00240 hess_diag.resize( pd.num_free_vertices(), 0.0 ); 00241 00242 double QM_val, QM_pow = 1.0; 00243 double fac1, fac2; 00244 const double negate_flag = qm->get_negate_flag(); 00245 bool qm_bool; 00246 size_t i; 00247 short p; 00248 00249 // Loops over all elements in the patch. 00250 OF_val = 0.0; 00251 std::vector< size_t >::const_iterator k; 00252 for( k = qmHandles.begin(); k != qmHandles.end(); ++k ) 00253 { 00254 // Computes \nabla^2 Q(e). Only the free vertices will have non-zero entries. 00255 qm_bool = qm->evaluate_with_Hessian_diagonal( pd, *k, QM_val, mIndices, mGradient, mDiag, err ); 00256 if( MSQ_CHKERR( err ) || !qm_bool ) return false; 00257 QM_val = fabs( QM_val ); 00258 00259 // **** Computes Hessian **** 00260 const size_t nve = mIndices.size(); 00261 if( pVal == 1 ) 00262 { 00263 QM_pow = 1.0; 00264 for( i = 0; i < nve; ++i ) 00265 { 00266 mDiag[i] *= negate_flag; 00267 hess_diag[mIndices[i]] += mDiag[i]; 00268 } 00269 fac1 = 1; 00270 } 00271 else if( pVal >= 2 ) 00272 { 00273 // Computes the coefficients: 00274 QM_pow = 1.0; 00275 for( p = 0; p < pVal - 2; ++p ) 00276 QM_pow *= QM_val; 00277 // 1 - computes p(p-1)Q(e)^{p-2} 00278 fac2 = pVal * ( pVal - 1 ) * QM_pow; 00279 // 2 - computes pQ(e)^{p-1} 00280 QM_pow *= QM_val; 00281 fac1 = pVal * QM_pow; 00282 00283 // fac1 *= qm->get_negate_flag(); 00284 // fac2 *= qm->get_negate_flag(); 00285 00286 for( i = 0; i < nve; ++i ) 00287 { 00288 SymMatrix3D op( mGradient[i] ); 00289 op *= fac2; 00290 mDiag[i] *= fac1; 00291 op += mDiag[i]; 00292 op *= negate_flag; 00293 hess_diag[mIndices[i]] += op; 00294 } 00295 } 00296 else 00297 { 00298 MSQ_SETERR( err )( " invalid P value.", MsqError::INVALID_STATE ); 00299 return false; 00300 } 00301 00302 // **** Computes Gradient **** 00303 00304 // For each vertex in the element ... 00305 for( i = 0; i < nve; ++i ) 00306 { 00307 // ... computes p*q^{p-1}*grad(q) ... 00308 mGradient[i] *= fac1 * negate_flag; 00309 // ... and accumulates it in the objective function gradient. 00310 // also scale the gradient by the scaling factor 00311 assert( mIndices[i] < pd.num_free_vertices() ); 00312 grad[mIndices[i]] += mGradient[i]; 00313 } 00314 00315 // **** computes Objective Function value \sum_{i=1}^{N_e} |q_i|^P **** 00316 OF_val += QM_pow * QM_val; 00317 } 00318 00319 size_t global_count; 00320 OF_val = negate_flag * get_value( OF_val, qmHandles.size(), type, global_count, err ); 00321 // if (!global_count) 00322 // return false; // invalid mesh 00323 00324 if( dividingByN && global_count ) 00325 { 00326 const double inv_n = 1.0 / global_count; 00327 for( i = 0; i < pd.num_free_vertices(); ++i ) 00328 { 00329 grad[i] *= inv_n; 00330 hess_diag[i] *= inv_n; 00331 } 00332 } 00333 00334 return true; 00335 } 00336 00337 /*\ For each element, each entry to be accumulated in the Hessian for 00338 this objective function (\f$ \sum_{e \in E} Q(e)^p \f$ where \f$ E \f$ 00339 is the set of all elements in the patch) has the form: 00340 \f$ pQ(e)^{p-1} \nabla^2 Q(e) + p(p-1)Q(e)^{p-2} \nabla Q(e) [\nabla Q(e)]^T \f$. 00341 00342 For \f$ p=2 \f$, this simplifies to 00343 \f$ 2Q(e) \nabla^2 Q(e) + 2 \nabla Q(e) [\nabla Q(e)]^T \f$. 00344 00345 For \f$ p=1 \f$, this simplifies to \f$ \nabla^2 Q(e) \f$. 00346 00347 The \f$ p=1 \f$ simplified version is implemented directly 00348 to speed up computation. 00349 00350 This function does not support vertex-based metrics. 00351 00352 \param pd The PatchData object for which the objective function 00353 hessian is computed. 00354 \param hessian this object must have been previously initialized. 00355 */ 00356 bool LPtoPTemplate::evaluate_with_Hessian( EvalType type, PatchData& pd, double& OF_val, std::vector< Vector3D >& grad, 00357 MsqHessian& hessian, MsqError& err ) 00358 { 00359 QualityMetric* qm = get_quality_metric(); 00360 qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err ); 00361 MSQ_ERRFALSE( err ); 00362 double negate_flag = qm->get_negate_flag(); 00363 00364 // zero gradient and hessian 00365 grad.clear(); 00366 grad.resize( pd.num_free_vertices(), 0.0 ); 00367 hessian.zero_out(); 00368 00369 double QM_val, QM_pow = 1.0; 00370 double fac1, fac2; 00371 Matrix3D elem_outer_product; 00372 bool qm_bool; 00373 size_t i, j, n; 00374 short p; 00375 00376 // Loops over all elements in the patch. 00377 OF_val = 0.0; 00378 std::vector< size_t >::const_iterator k; 00379 for( k = qmHandles.begin(); k != qmHandles.end(); ++k ) 00380 { 00381 // Computes \nabla^2 Q(e). Only the free vertices will have non-zero entries. 00382 qm_bool = qm->evaluate_with_Hessian( pd, *k, QM_val, mIndices, mGradient, mHessian, err ); 00383 if( MSQ_CHKERR( err ) || !qm_bool ) return false; 00384 QM_val = fabs( QM_val ); 00385 00386 // **** Computes Hessian **** 00387 const size_t nve = mIndices.size(); 00388 if( pVal == 1 ) 00389 { 00390 QM_pow = 1.0; 00391 n = 0; 00392 for( i = 0; i < nve; ++i ) 00393 { 00394 for( j = i; j < nve; ++j ) 00395 { 00396 // negate if necessary 00397 mHessian[n] *= negate_flag; 00398 hessian.add( mIndices[i], mIndices[j], mHessian[n], err ); 00399 MSQ_ERRFALSE( err ); 00400 ++n; 00401 } 00402 } 00403 fac1 = 1; 00404 } 00405 else if( pVal >= 2 ) 00406 { 00407 // Computes the coefficients: 00408 QM_pow = 1.0; 00409 for( p = 0; p < pVal - 2; ++p ) 00410 QM_pow *= QM_val; 00411 // 1 - computes p(p-1)Q(e)^{p-2} 00412 fac2 = pVal * ( pVal - 1 ) * QM_pow; 00413 // 2 - computes pQ(e)^{p-1} 00414 QM_pow *= QM_val; 00415 fac1 = pVal * QM_pow; 00416 00417 // fac1 *= qm->get_negate_flag(); 00418 // fac2 *= qm->get_negate_flag(); 00419 00420 n = 0; 00421 for( i = 0; i < nve; ++i ) 00422 { 00423 for( j = i; j < nve; ++j ) 00424 { 00425 elem_outer_product.outer_product( mGradient[i], mGradient[j] ); 00426 elem_outer_product *= fac2; 00427 mHessian[n] *= fac1; 00428 mHessian[n] += elem_outer_product; 00429 mHessian[n] *= negate_flag; 00430 hessian.add( mIndices[i], mIndices[j], mHessian[n], err ); 00431 MSQ_ERRFALSE( err ); 00432 ++n; 00433 } 00434 } 00435 } 00436 else 00437 { 00438 MSQ_SETERR( err )( " invalid P value.", MsqError::INVALID_STATE ); 00439 return false; 00440 } 00441 00442 // **** Computes Gradient **** 00443 00444 // For each vertex in the element ... 00445 for( i = 0; i < nve; ++i ) 00446 { 00447 // ... computes p*q^{p-1}*grad(q) ... 00448 mGradient[i] *= fac1 * negate_flag; 00449 // ... and accumulates it in the objective function gradient. 00450 // also scale the gradient by the scaling factor 00451 assert( mIndices[i] < pd.num_free_vertices() ); 00452 grad[mIndices[i]] += mGradient[i]; 00453 } 00454 00455 // **** computes Objective Function value \sum_{i=1}^{N_e} |q_i|^P **** 00456 OF_val += QM_pow * QM_val; 00457 } 00458 00459 size_t global_count; 00460 OF_val = negate_flag * get_value( OF_val, qmHandles.size(), type, global_count, err ); 00461 // if (!global_count) 00462 // return false; // invalid mesh 00463 00464 if( dividingByN && global_count ) 00465 { 00466 const double inv_n = 1.0 / global_count; 00467 std::vector< Vector3D >::iterator g; 00468 for( g = grad.begin(); g != grad.end(); ++g ) 00469 *g *= inv_n; 00470 hessian.scale( inv_n ); 00471 } 00472 00473 return true; 00474 }