MOAB: Mesh Oriented datABase  (version 5.4.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,
00049                                                            MeshDomain* domain,
00050                                                            const Settings* settings,
00051                                                            PatchSet* patch_set,
00052                                                            MsqError& err )
00053 {
00054     clear();
00055     PatchIterator patches( patch_set );
00056 
00057     PatchData pd;
00058     pd.set_mesh( mesh );
00059     pd.set_domain( domain );
00060     if( settings ) pd.attach_settings( settings );
00061 
00062     bool result = true;
00063     while( patches.get_next_patch( pd, err ) && !MSQ_CHKERR( err ) )
00064     {
00065         double value;
00066         bool b = evaluate( ObjectiveFunction::ACCUMULATE, pd, value, false, err );
00067         MSQ_ERRZERO( err );
00068         result = result && b;
00069     }
00070     return result;
00071 }
00072 
00073 bool PatchPowerMeanP::evaluate( EvalType type, PatchData& pd, double& value_out, bool free, MsqError& err )
00074 {
00075     QualityMetric* qm = get_quality_metric();
00076     if( type == ObjectiveFunction::ACCUMULATE )
00077         qm->get_single_pass( pd, qmHandles, free, err );
00078     else
00079         qm->get_evaluations( pd, qmHandles, free, err );
00080     MSQ_ERRFALSE( err );
00081 
00082     // calculate OF value for just the patch
00083     std::vector< size_t >::const_iterator i;
00084     double value, working_sum = 0.0;
00085     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00086     {
00087         bool result = qm->evaluate( pd, *i, value, err );
00088         if( MSQ_CHKERR( err ) || !result ) return false;
00089 
00090         working_sum += mPower.raise( value );
00091     }
00092     working_sum /= qmHandles.size();
00093 
00094     // get overall OF value, update member data, etc.
00095     size_t global_count;
00096     value_out = qm->get_negate_flag() * get_value( working_sum, 1, type, global_count );
00097     return true;
00098 }
00099 
00100 bool PatchPowerMeanP::evaluate_with_gradient( EvalType type,
00101                                               PatchData& pd,
00102                                               double& value_out,
00103                                               std::vector< Vector3D >& grad_out,
00104                                               MsqError& err )
00105 {
00106     QualityMetric* qm = get_quality_metric();
00107     qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );
00108     MSQ_ERRFALSE( err );
00109 
00110     // zero gradient
00111     grad_out.clear();
00112     grad_out.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00113 
00114     // calculate OF value and gradient for just the patch
00115     std::vector< size_t >::const_iterator i;
00116     double value, working_sum = 0.0;
00117     const double f = qm->get_negate_flag() * mPower.value() / qmHandles.size();
00118     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00119     {
00120         bool result = qm->evaluate_with_gradient( pd, *i, value, mIndices, mGradient, err );
00121         if( MSQ_CHKERR( err ) || !result ) return false;
00122         if( fabs( value ) < DBL_EPSILON ) continue;
00123 
00124         const double qmp = mPower.raise( value );
00125         working_sum += mPower.raise( value );
00126         value = f * qmp / value;
00127 
00128         for( size_t j = 0; j < mIndices.size(); ++j )
00129         {
00130             mGradient[j] *= value;
00131             grad_out[mIndices[j]] += mGradient[j];
00132         }
00133     }
00134 
00135     // get overall OF value, update member data, etc.
00136     size_t global_count;
00137     value_out          = qm->get_negate_flag() * get_value( working_sum, 1, type, global_count );
00138     const double inv_n = 1.0 / global_count;
00139     std::vector< Vector3D >::iterator g;
00140     for( g = grad_out.begin(); g != grad_out.end(); ++g )
00141         *g *= inv_n;
00142     return true;
00143 }
00144 
00145 bool PatchPowerMeanP::evaluate_with_Hessian( EvalType type,
00146                                              PatchData& pd,
00147                                              double& value_out,
00148                                              std::vector< Vector3D >& grad_out,
00149                                              MsqHessian& Hessian_out,
00150                                              MsqError& err )
00151 {
00152     QualityMetric* qm = get_quality_metric();
00153     qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );
00154     MSQ_ERRFALSE( err );
00155 
00156     // zero gradient and hessian
00157     grad_out.clear();
00158     grad_out.resize( pd.num_free_vertices(), 0.0 );
00159     Hessian_out.zero_out();
00160 
00161     // calculate OF value and gradient for just the patch
00162     std::vector< size_t >::const_iterator i;
00163     size_t j, k, n;
00164     double value, working_sum = 0.0;
00165     const double f1 = qm->get_negate_flag() * mPower.value() / qmHandles.size();
00166     const double f2 = f1 * ( mPower.value() - 1 );
00167     Matrix3D m;
00168     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00169     {
00170         bool result = qm->evaluate_with_Hessian( pd, *i, value, mIndices, mGradient, mHessian, err );
00171         if( MSQ_CHKERR( err ) || !result ) return false;
00172         if( fabs( value ) < DBL_EPSILON ) continue;
00173 
00174         const double qmp = mPower.raise( value );
00175         const double hf  = f2 * qmp / ( value * value );
00176         const double gf  = f1 * qmp / value;
00177         working_sum += qmp;
00178 
00179         const size_t nfree = mIndices.size();
00180         n                  = 0;
00181         for( j = 0; j < nfree; ++j )
00182         {
00183             for( k = j; k < nfree; ++k )
00184             {
00185                 m.outer_product( mGradient[j], mGradient[k] );
00186                 m *= hf;
00187                 mHessian[n] *= gf;
00188                 m += mHessian[n];
00189                 ++n;
00190                 Hessian_out.add( mIndices[j], mIndices[k], m, err );
00191                 MSQ_ERRFALSE( err );
00192             }
00193         }
00194 
00195         for( j = 0; j < nfree; ++j )
00196         {
00197             mGradient[j] *= gf;
00198             grad_out[mIndices[j]] += mGradient[j];
00199         }
00200     }
00201 
00202     // get overall OF value, update member data, etc.
00203     size_t global_count;
00204     value_out          = qm->get_negate_flag() * get_value( working_sum, 1, type, global_count );
00205     const double inv_n = 1.0 / global_count;
00206     std::vector< Vector3D >::iterator g;
00207     for( g = grad_out.begin(); g != grad_out.end(); ++g )
00208         *g *= inv_n;
00209     Hessian_out.scale( inv_n );
00210     return true;
00211 }
00212 
00213 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines