MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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