MOAB: Mesh Oriented datABase  (version 5.3.0)
LPtoPTemplate.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines