MOAB: Mesh Oriented datABase  (version 5.3.0)
PatchPowerMeanP.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retain certain rights to 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 PatchPowerMeanP.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "PatchPowerMeanP.hpp"
00034 #include "QualityMetric.hpp"
00035 #include "MsqError.hpp"
00036 #include "MsqHessian.hpp"
00037 #include "PatchData.hpp"
00038 #include "PatchIterator.hpp"
00039 
00040 namespace MBMesquite
00041 {
00042 
00043 ObjectiveFunction* PatchPowerMeanP::clone() const
00044 {
00045     return new PatchPowerMeanP( *this );
00046 }
00047 
00048 bool PatchPowerMeanP::initialize_block_coordinate_descent( Mesh* mesh, MeshDomain* domain, const Settings* settings,
00049                                                            PatchSet* patch_set, MsqError& err )
00050 {
00051     clear();
00052     PatchIterator patches( patch_set );
00053 
00054     PatchData pd;
00055     pd.set_mesh( mesh );
00056     pd.set_domain( domain );
00057     if( settings ) pd.attach_settings( settings );
00058 
00059     bool result = true;
00060     while( patches.get_next_patch( pd, err ) && !MSQ_CHKERR( err ) )
00061     {
00062         double value;
00063         bool b = evaluate( ObjectiveFunction::ACCUMULATE, pd, value, false, err );
00064         MSQ_ERRZERO( err );
00065         result = result && b;
00066     }
00067     return result;
00068 }
00069 
00070 bool PatchPowerMeanP::evaluate( EvalType type, PatchData& pd, double& value_out, bool free, MsqError& err )
00071 {
00072     QualityMetric* qm = get_quality_metric();
00073     if( type == ObjectiveFunction::ACCUMULATE )
00074         qm->get_single_pass( pd, qmHandles, free, err );
00075     else
00076         qm->get_evaluations( pd, qmHandles, free, err );
00077     MSQ_ERRFALSE( err );
00078 
00079     // calculate OF value for just the patch
00080     std::vector< size_t >::const_iterator i;
00081     double value, working_sum = 0.0;
00082     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00083     {
00084         bool result = qm->evaluate( pd, *i, value, err );
00085         if( MSQ_CHKERR( err ) || !result ) return false;
00086 
00087         working_sum += mPower.raise( value );
00088     }
00089     working_sum /= qmHandles.size();
00090 
00091     // get overall OF value, update member data, etc.
00092     size_t global_count;
00093     value_out = qm->get_negate_flag() * get_value( working_sum, 1, type, global_count );
00094     return true;
00095 }
00096 
00097 bool PatchPowerMeanP::evaluate_with_gradient( EvalType type, PatchData& pd, double& value_out,
00098                                               std::vector< Vector3D >& grad_out, MsqError& err )
00099 {
00100     QualityMetric* qm = get_quality_metric();
00101     qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );
00102     MSQ_ERRFALSE( err );
00103 
00104     // zero gradient
00105     grad_out.clear();
00106     grad_out.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00107 
00108     // calculate OF value and gradient for just the patch
00109     std::vector< size_t >::const_iterator i;
00110     double value, working_sum = 0.0;
00111     const double f = qm->get_negate_flag() * mPower.value() / qmHandles.size();
00112     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00113     {
00114         bool result = qm->evaluate_with_gradient( pd, *i, value, mIndices, mGradient, err );
00115         if( MSQ_CHKERR( err ) || !result ) return false;
00116         if( fabs( value ) < DBL_EPSILON ) continue;
00117 
00118         const double qmp = mPower.raise( value );
00119         working_sum += mPower.raise( value );
00120         value = f * qmp / value;
00121 
00122         for( size_t j = 0; j < mIndices.size(); ++j )
00123         {
00124             mGradient[j] *= value;
00125             grad_out[mIndices[j]] += mGradient[j];
00126         }
00127     }
00128 
00129     // get overall OF value, update member data, etc.
00130     size_t global_count;
00131     value_out          = qm->get_negate_flag() * get_value( working_sum, 1, type, global_count );
00132     const double inv_n = 1.0 / global_count;
00133     std::vector< Vector3D >::iterator g;
00134     for( g = grad_out.begin(); g != grad_out.end(); ++g )
00135         *g *= inv_n;
00136     return true;
00137 }
00138 
00139 bool PatchPowerMeanP::evaluate_with_Hessian( EvalType type, PatchData& pd, double& value_out,
00140                                              std::vector< Vector3D >& grad_out, MsqHessian& Hessian_out, MsqError& err )
00141 {
00142     QualityMetric* qm = get_quality_metric();
00143     qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );
00144     MSQ_ERRFALSE( err );
00145 
00146     // zero gradient and hessian
00147     grad_out.clear();
00148     grad_out.resize( pd.num_free_vertices(), 0.0 );
00149     Hessian_out.zero_out();
00150 
00151     // calculate OF value and gradient for just the patch
00152     std::vector< size_t >::const_iterator i;
00153     size_t j, k, n;
00154     double value, working_sum = 0.0;
00155     const double f1 = qm->get_negate_flag() * mPower.value() / qmHandles.size();
00156     const double f2 = f1 * ( mPower.value() - 1 );
00157     Matrix3D m;
00158     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00159     {
00160         bool result = qm->evaluate_with_Hessian( pd, *i, value, mIndices, mGradient, mHessian, err );
00161         if( MSQ_CHKERR( err ) || !result ) return false;
00162         if( fabs( value ) < DBL_EPSILON ) continue;
00163 
00164         const double qmp = mPower.raise( value );
00165         const double hf  = f2 * qmp / ( value * value );
00166         const double gf  = f1 * qmp / value;
00167         working_sum += qmp;
00168 
00169         const size_t nfree = mIndices.size();
00170         n                  = 0;
00171         for( j = 0; j < nfree; ++j )
00172         {
00173             for( k = j; k < nfree; ++k )
00174             {
00175                 m.outer_product( mGradient[j], mGradient[k] );
00176                 m *= hf;
00177                 mHessian[n] *= gf;
00178                 m += mHessian[n];
00179                 ++n;
00180                 Hessian_out.add( mIndices[j], mIndices[k], m, err );
00181                 MSQ_ERRFALSE( err );
00182             }
00183         }
00184 
00185         for( j = 0; j < nfree; ++j )
00186         {
00187             mGradient[j] *= gf;
00188             grad_out[mIndices[j]] += mGradient[j];
00189         }
00190     }
00191 
00192     // get overall OF value, update member data, etc.
00193     size_t global_count;
00194     value_out          = qm->get_negate_flag() * get_value( working_sum, 1, type, global_count );
00195     const double inv_n = 1.0 / global_count;
00196     std::vector< Vector3D >::iterator g;
00197     for( g = grad_out.begin(); g != grad_out.end(); ++g )
00198         *g *= inv_n;
00199     Hessian_out.scale( inv_n );
00200     return true;
00201 }
00202 
00203 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines