MOAB: Mesh Oriented datABase  (version 5.2.1)
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, 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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines