MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2006 Lawrence Livermore National Laboratory. Under 00005 the terms of Contract B545069 with the University of Wisconsin -- 00006 Madison, Lawrence Livermore National Laboratory 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 (2006) kraftche@cae.wisc.edu 00024 00025 ***************************************************************** */ 00026 00027 /** \file PMeanPTemplate.cpp 00028 * \brief previous name: PowerMeanP.cpp 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #include "Mesquite.hpp" 00033 #include "PMeanPTemplate.hpp" 00034 #include "QualityMetric.hpp" 00035 #include "MsqError.hpp" 00036 #include "MsqHessian.hpp" 00037 #include "PatchData.hpp" 00038 00039 namespace MBMesquite 00040 { 00041 00042 ObjectiveFunction* PMeanPTemplate::clone() const 00043 { 00044 return new PMeanPTemplate( *this ); 00045 } 00046 00047 void PMeanPTemplate::clear() 00048 { 00049 mCount = 0; 00050 mPowSum = 0; 00051 saveCount = 0; 00052 savePowSum = 0; 00053 } 00054 00055 double PMeanPTemplate::get_value( double power_sum, size_t count, EvalType type, size_t& global_count ) 00056 { 00057 double result = 0; 00058 switch( type ) 00059 { 00060 case CALCULATE: 00061 result = power_sum; 00062 global_count = count; 00063 break; 00064 00065 case ACCUMULATE: 00066 mPowSum += power_sum; 00067 mCount += count; 00068 result = mPowSum; 00069 global_count = mCount; 00070 break; 00071 00072 case SAVE: 00073 savePowSum = power_sum; 00074 saveCount = count; 00075 result = mPowSum; 00076 global_count = mCount; 00077 break; 00078 00079 case UPDATE: 00080 mPowSum -= savePowSum; 00081 mCount -= saveCount; 00082 savePowSum = power_sum; 00083 saveCount = count; 00084 mPowSum += savePowSum; 00085 mCount += saveCount; 00086 result = mPowSum; 00087 global_count = mCount; 00088 break; 00089 00090 case TEMPORARY: 00091 result = mPowSum - savePowSum + power_sum; 00092 global_count = mCount + count - saveCount; 00093 break; 00094 } 00095 00096 return global_count ? result / global_count : 0.0; 00097 } 00098 00099 bool PMeanPTemplate::evaluate( EvalType type, PatchData& pd, double& value_out, bool free, MsqError& err ) 00100 { 00101 QualityMetric* qm = get_quality_metric(); 00102 if( type == ObjectiveFunction::ACCUMULATE ) 00103 qm->get_single_pass( pd, qmHandles, free, err ); 00104 else 00105 qm->get_evaluations( pd, qmHandles, free, err ); 00106 MSQ_ERRFALSE( err ); 00107 00108 // calculate OF value for just the patch 00109 std::vector< size_t >::const_iterator i; 00110 double value, working_sum = 0.0; 00111 for( i = qmHandles.begin(); i != qmHandles.end(); ++i ) 00112 { 00113 bool result = qm->evaluate( pd, *i, value, err ); 00114 if( MSQ_CHKERR( err ) || !result ) return false; 00115 00116 working_sum += mPower.raise( value ); 00117 } 00118 00119 // get overall OF value, update member data, etc. 00120 size_t global_count = 0; 00121 value_out = qm->get_negate_flag() * get_value( working_sum, qmHandles.size(), type, global_count ); 00122 return true; 00123 } 00124 00125 bool PMeanPTemplate::evaluate_with_gradient( EvalType type, PatchData& pd, double& value_out, 00126 std::vector< Vector3D >& grad_out, MsqError& err ) 00127 { 00128 QualityMetric* qm = get_quality_metric(); 00129 qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err ); 00130 MSQ_ERRFALSE( err ); 00131 00132 // zero gradient 00133 grad_out.clear(); 00134 grad_out.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) ); 00135 00136 // calculate OF value and gradient for just the patch 00137 std::vector< size_t >::const_iterator i; 00138 double value, working_sum = 0.0; 00139 const double f = qm->get_negate_flag() * mPower.value(); 00140 for( i = qmHandles.begin(); i != qmHandles.end(); ++i ) 00141 { 00142 bool result = qm->evaluate_with_gradient( pd, *i, value, mIndices, mGradient, err ); 00143 if( MSQ_CHKERR( err ) || !result ) return false; 00144 if( fabs( value ) < DBL_EPSILON ) continue; 00145 00146 const double r1 = mPowerMinus1.raise( value ); 00147 const double qmp = r1 * value; 00148 working_sum += qmp; 00149 value = f * r1; 00150 00151 for( size_t j = 0; j < mIndices.size(); ++j ) 00152 { 00153 mGradient[j] *= value; 00154 grad_out[mIndices[j]] += mGradient[j]; 00155 } 00156 } 00157 00158 // get overall OF value, update member data, etc. 00159 size_t global_count = 0; 00160 value_out = qm->get_negate_flag() * get_value( working_sum, qmHandles.size(), type, global_count ); 00161 if( global_count ) 00162 { 00163 const double inv_n = 1.0 / global_count; 00164 std::vector< Vector3D >::iterator g; 00165 for( g = grad_out.begin(); g != grad_out.end(); ++g ) 00166 *g *= inv_n; 00167 } 00168 return true; 00169 } 00170 00171 bool PMeanPTemplate::evaluate_with_Hessian_diagonal( EvalType type, PatchData& pd, double& value_out, 00172 std::vector< Vector3D >& grad_out, 00173 std::vector< SymMatrix3D >& hess_diag_out, MsqError& err ) 00174 { 00175 QualityMetric* qm = get_quality_metric(); 00176 qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err ); 00177 MSQ_ERRFALSE( err ); 00178 00179 // zero gradient and hessian 00180 const size_t s = pd.num_free_vertices(); 00181 grad_out.clear(); 00182 grad_out.resize( s, 0.0 ); 00183 hess_diag_out.clear(); 00184 hess_diag_out.resize( s, 0.0 ); 00185 00186 // calculate OF value and gradient for just the patch 00187 std::vector< size_t >::const_iterator i; 00188 size_t j; 00189 double value, working_sum = 0.0; 00190 const double f1 = qm->get_negate_flag() * mPower.value(); 00191 const double f2 = f1 * ( mPower.value() - 1 ); 00192 for( i = qmHandles.begin(); i != qmHandles.end(); ++i ) 00193 { 00194 bool result = qm->evaluate_with_Hessian_diagonal( pd, *i, value, mIndices, mGradient, mDiag, err ); 00195 if( MSQ_CHKERR( err ) || !result ) return false; 00196 if( fabs( value ) < DBL_EPSILON ) continue; 00197 00198 const size_t nfree = mIndices.size(); 00199 if( mPower.value() == 1.0 ) 00200 { 00201 working_sum += mPower.raise( value ); 00202 for( j = 0; j < nfree; ++j ) 00203 { 00204 const size_t idx = mIndices[j]; 00205 hess_diag_out[idx] += f1 * mDiag[j]; 00206 mGradient[j] *= f1; 00207 grad_out[idx] += mGradient[j]; 00208 } 00209 } 00210 else 00211 { 00212 const double r2 = mPowerMinus2.raise( value ); 00213 const double r1 = r2 * value; 00214 working_sum += r1 * value; 00215 const double hf = f2 * r2; 00216 const double gf = f1 * r1; 00217 for( j = 0; j < nfree; ++j ) 00218 { 00219 const size_t idx = mIndices[j]; 00220 00221 hess_diag_out[idx] += hf * outer( mGradient[j] ); 00222 hess_diag_out[idx] += gf * mDiag[j]; 00223 00224 mGradient[j] *= gf; 00225 grad_out[idx] += mGradient[j]; 00226 } 00227 } 00228 } 00229 00230 // get overall OF value, update member data, etc. 00231 size_t global_count = 0; 00232 value_out = qm->get_negate_flag() * get_value( working_sum, qmHandles.size(), type, global_count ); 00233 if( global_count ) 00234 { 00235 const double inv_n = 1.0 / global_count; 00236 for( j = 0; j < s; ++j ) 00237 { 00238 grad_out[j] *= inv_n; 00239 hess_diag_out[j] *= inv_n; 00240 } 00241 } 00242 return true; 00243 } 00244 00245 bool PMeanPTemplate::evaluate_with_Hessian( EvalType type, PatchData& pd, double& value_out, 00246 std::vector< Vector3D >& grad_out, MsqHessian& Hessian_out, MsqError& err ) 00247 { 00248 QualityMetric* qm = get_quality_metric(); 00249 qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err ); 00250 MSQ_ERRFALSE( err ); 00251 00252 // zero gradient and hessian 00253 grad_out.clear(); 00254 grad_out.resize( pd.num_free_vertices(), 0.0 ); 00255 Hessian_out.zero_out(); 00256 00257 // calculate OF value and gradient for just the patch 00258 std::vector< size_t >::const_iterator i; 00259 size_t j, k, n; 00260 double value, working_sum = 0.0; 00261 const double f1 = qm->get_negate_flag() * mPower.value(); 00262 const double f2 = f1 * ( mPower.value() - 1 ); 00263 Matrix3D m; 00264 for( i = qmHandles.begin(); i != qmHandles.end(); ++i ) 00265 { 00266 bool result = qm->evaluate_with_Hessian( pd, *i, value, mIndices, mGradient, mHessian, err ); 00267 if( MSQ_CHKERR( err ) || !result ) return false; 00268 if( fabs( value ) < DBL_EPSILON ) continue; 00269 00270 const size_t nfree = mIndices.size(); 00271 n = 0; 00272 if( mPower.value() == 1.0 ) 00273 { 00274 working_sum += mPower.raise( value ); 00275 for( j = 0; j < nfree; ++j ) 00276 { 00277 mGradient[j] *= f1; 00278 grad_out[mIndices[j]] += mGradient[j]; 00279 for( k = j; k < nfree; ++k ) 00280 { 00281 mHessian[n] *= f1; 00282 Hessian_out.add( mIndices[j], mIndices[k], mHessian[n], err ); 00283 MSQ_ERRFALSE( err ); 00284 ++n; 00285 } 00286 } 00287 } 00288 else 00289 { 00290 const double r2 = mPowerMinus2.raise( value ); 00291 const double r1 = r2 * value; 00292 working_sum += r1 * value; 00293 const double hf = f2 * r2; 00294 const double gf = f1 * r1; 00295 for( j = 0; j < nfree; ++j ) 00296 { 00297 for( k = j; k < nfree; ++k ) 00298 { 00299 m.outer_product( mGradient[j], mGradient[k] ); 00300 m *= hf; 00301 mHessian[n] *= gf; 00302 m += mHessian[n]; 00303 Hessian_out.add( mIndices[j], mIndices[k], m, err ); 00304 MSQ_ERRFALSE( err ); 00305 ++n; 00306 } 00307 } 00308 for( j = 0; j < nfree; ++j ) 00309 { 00310 mGradient[j] *= gf; 00311 grad_out[mIndices[j]] += mGradient[j]; 00312 } 00313 } 00314 } 00315 00316 // get overall OF value, update member data, etc. 00317 size_t global_count = 0; 00318 value_out = qm->get_negate_flag() * get_value( working_sum, qmHandles.size(), type, global_count ); 00319 if( global_count ) 00320 { 00321 const double inv_n = 1.0 / global_count; 00322 std::vector< Vector3D >::iterator g; 00323 for( g = grad_out.begin(); g != grad_out.end(); ++g ) 00324 *g *= inv_n; 00325 Hessian_out.scale( inv_n ); 00326 } 00327 return true; 00328 } 00329 00330 } // namespace MBMesquite