MOAB: Mesh Oriented datABase  (version 5.4.0)
PMeanPTemplate.cpp
Go to the documentation of this file.
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,
00126                                              PatchData& pd,
00127                                              double& value_out,
00128                                              std::vector< Vector3D >& grad_out,
00129                                              MsqError& err )
00130 {
00131     QualityMetric* qm = get_quality_metric();
00132     qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );
00133     MSQ_ERRFALSE( err );
00134 
00135     // zero gradient
00136     grad_out.clear();
00137     grad_out.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00138 
00139     // calculate OF value and gradient for just the patch
00140     std::vector< size_t >::const_iterator i;
00141     double value, working_sum = 0.0;
00142     const double f = qm->get_negate_flag() * mPower.value();
00143     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00144     {
00145         bool result = qm->evaluate_with_gradient( pd, *i, value, mIndices, mGradient, err );
00146         if( MSQ_CHKERR( err ) || !result ) return false;
00147         if( fabs( value ) < DBL_EPSILON ) continue;
00148 
00149         const double r1  = mPowerMinus1.raise( value );
00150         const double qmp = r1 * value;
00151         working_sum += qmp;
00152         value = f * r1;
00153 
00154         for( size_t j = 0; j < mIndices.size(); ++j )
00155         {
00156             mGradient[j] *= value;
00157             grad_out[mIndices[j]] += mGradient[j];
00158         }
00159     }
00160 
00161     // get overall OF value, update member data, etc.
00162     size_t global_count = 0;
00163     value_out           = qm->get_negate_flag() * get_value( working_sum, qmHandles.size(), type, global_count );
00164     if( global_count )
00165     {
00166         const double inv_n = 1.0 / global_count;
00167         std::vector< Vector3D >::iterator g;
00168         for( g = grad_out.begin(); g != grad_out.end(); ++g )
00169             *g *= inv_n;
00170     }
00171     return true;
00172 }
00173 
00174 bool PMeanPTemplate::evaluate_with_Hessian_diagonal( EvalType type,
00175                                                      PatchData& pd,
00176                                                      double& value_out,
00177                                                      std::vector< Vector3D >& grad_out,
00178                                                      std::vector< SymMatrix3D >& hess_diag_out,
00179                                                      MsqError& err )
00180 {
00181     QualityMetric* qm = get_quality_metric();
00182     qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );
00183     MSQ_ERRFALSE( err );
00184 
00185     // zero gradient and hessian
00186     const size_t s = pd.num_free_vertices();
00187     grad_out.clear();
00188     grad_out.resize( s, 0.0 );
00189     hess_diag_out.clear();
00190     hess_diag_out.resize( s, 0.0 );
00191 
00192     // calculate OF value and gradient for just the patch
00193     std::vector< size_t >::const_iterator i;
00194     size_t j;
00195     double value, working_sum = 0.0;
00196     const double f1 = qm->get_negate_flag() * mPower.value();
00197     const double f2 = f1 * ( mPower.value() - 1 );
00198     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00199     {
00200         bool result = qm->evaluate_with_Hessian_diagonal( pd, *i, value, mIndices, mGradient, mDiag, err );
00201         if( MSQ_CHKERR( err ) || !result ) return false;
00202         if( fabs( value ) < DBL_EPSILON ) continue;
00203 
00204         const size_t nfree = mIndices.size();
00205         if( mPower.value() == 1.0 )
00206         {
00207             working_sum += mPower.raise( value );
00208             for( j = 0; j < nfree; ++j )
00209             {
00210                 const size_t idx = mIndices[j];
00211                 hess_diag_out[idx] += f1 * mDiag[j];
00212                 mGradient[j] *= f1;
00213                 grad_out[idx] += mGradient[j];
00214             }
00215         }
00216         else
00217         {
00218             const double r2 = mPowerMinus2.raise( value );
00219             const double r1 = r2 * value;
00220             working_sum += r1 * value;
00221             const double hf = f2 * r2;
00222             const double gf = f1 * r1;
00223             for( j = 0; j < nfree; ++j )
00224             {
00225                 const size_t idx = mIndices[j];
00226 
00227                 hess_diag_out[idx] += hf * outer( mGradient[j] );
00228                 hess_diag_out[idx] += gf * mDiag[j];
00229 
00230                 mGradient[j] *= gf;
00231                 grad_out[idx] += mGradient[j];
00232             }
00233         }
00234     }
00235 
00236     // get overall OF value, update member data, etc.
00237     size_t global_count = 0;
00238     value_out           = qm->get_negate_flag() * get_value( working_sum, qmHandles.size(), type, global_count );
00239     if( global_count )
00240     {
00241         const double inv_n = 1.0 / global_count;
00242         for( j = 0; j < s; ++j )
00243         {
00244             grad_out[j] *= inv_n;
00245             hess_diag_out[j] *= inv_n;
00246         }
00247     }
00248     return true;
00249 }
00250 
00251 bool PMeanPTemplate::evaluate_with_Hessian( EvalType type,
00252                                             PatchData& pd,
00253                                             double& value_out,
00254                                             std::vector< Vector3D >& grad_out,
00255                                             MsqHessian& Hessian_out,
00256                                             MsqError& err )
00257 {
00258     QualityMetric* qm = get_quality_metric();
00259     qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );
00260     MSQ_ERRFALSE( err );
00261 
00262     // zero gradient and hessian
00263     grad_out.clear();
00264     grad_out.resize( pd.num_free_vertices(), 0.0 );
00265     Hessian_out.zero_out();
00266 
00267     // calculate OF value and gradient for just the patch
00268     std::vector< size_t >::const_iterator i;
00269     size_t j, k, n;
00270     double value, working_sum = 0.0;
00271     const double f1 = qm->get_negate_flag() * mPower.value();
00272     const double f2 = f1 * ( mPower.value() - 1 );
00273     Matrix3D m;
00274     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00275     {
00276         bool result = qm->evaluate_with_Hessian( pd, *i, value, mIndices, mGradient, mHessian, err );
00277         if( MSQ_CHKERR( err ) || !result ) return false;
00278         if( fabs( value ) < DBL_EPSILON ) continue;
00279 
00280         const size_t nfree = mIndices.size();
00281         n                  = 0;
00282         if( mPower.value() == 1.0 )
00283         {
00284             working_sum += mPower.raise( value );
00285             for( j = 0; j < nfree; ++j )
00286             {
00287                 mGradient[j] *= f1;
00288                 grad_out[mIndices[j]] += mGradient[j];
00289                 for( k = j; k < nfree; ++k )
00290                 {
00291                     mHessian[n] *= f1;
00292                     Hessian_out.add( mIndices[j], mIndices[k], mHessian[n], err );
00293                     MSQ_ERRFALSE( err );
00294                     ++n;
00295                 }
00296             }
00297         }
00298         else
00299         {
00300             const double r2 = mPowerMinus2.raise( value );
00301             const double r1 = r2 * value;
00302             working_sum += r1 * value;
00303             const double hf = f2 * r2;
00304             const double gf = f1 * r1;
00305             for( j = 0; j < nfree; ++j )
00306             {
00307                 for( k = j; k < nfree; ++k )
00308                 {
00309                     m.outer_product( mGradient[j], mGradient[k] );
00310                     m *= hf;
00311                     mHessian[n] *= gf;
00312                     m += mHessian[n];
00313                     Hessian_out.add( mIndices[j], mIndices[k], m, err );
00314                     MSQ_ERRFALSE( err );
00315                     ++n;
00316                 }
00317             }
00318             for( j = 0; j < nfree; ++j )
00319             {
00320                 mGradient[j] *= gf;
00321                 grad_out[mIndices[j]] += mGradient[j];
00322             }
00323         }
00324     }
00325 
00326     // get overall OF value, update member data, etc.
00327     size_t global_count = 0;
00328     value_out           = qm->get_negate_flag() * get_value( working_sum, qmHandles.size(), type, global_count );
00329     if( global_count )
00330     {
00331         const double inv_n = 1.0 / global_count;
00332         std::vector< Vector3D >::iterator g;
00333         for( g = grad_out.begin(); g != grad_out.end(); ++g )
00334             *g *= inv_n;
00335         Hessian_out.scale( inv_n );
00336     }
00337     return true;
00338 }
00339 
00340 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines