MOAB: Mesh Oriented datABase
(version 5.3.1)
|
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 ObjectiveFunction.cpp 00029 \brief 00030 00031 \author Michael Brewer 00032 \author Thomas Leurent 00033 00034 \date 2002-08-02 00035 */ 00036 00037 #include "ObjectiveFunction.hpp" 00038 #include "MsqVertex.hpp" 00039 #include "MsqDebug.hpp" 00040 #include "PatchData.hpp" 00041 #include "MsqError.hpp" 00042 #include "MsqHessian.hpp" 00043 #include "SymMatrix3D.hpp" 00044 #include <memory> // unique_ptr 00045 00046 namespace MBMesquite 00047 { 00048 00049 ObjectiveFunction::~ObjectiveFunction() {} 00050 00051 /*!Returns an appropiate value (eps) to use as a delta step for 00052 MsqVertex vertex in dimension k (i.e. k=0 -> x, k=1 -> y, k=2 -> z). 00053 The objective function value at the perturbed vertex position is given 00054 in local_val. 00055 */ 00056 double ObjectiveFunction::get_eps( PatchData& pd, EvalType type, double& local_val, int dim, size_t vertex_index, 00057 MsqError& err ) 00058 { 00059 double eps = 1.e-07; 00060 const double rho = 0.5; 00061 const int imax = 20; 00062 bool feasible = false; 00063 double tmp_var = 0.0; 00064 Vector3D delta( 0, 0, 0 ); 00065 for( int i = 0; i < imax; ++i ) 00066 { 00067 // perturb kth coord val and check feas if needed 00068 tmp_var = pd.vertex_by_index( vertex_index )[dim]; 00069 delta[dim] = eps; 00070 pd.move_vertex( delta, vertex_index, err ); 00071 feasible = evaluate( type, pd, local_val, OF_FREE_EVALS_ONLY, err ); 00072 MSQ_ERRZERO( err ); 00073 // revert kth coord val 00074 delta = pd.vertex_by_index( vertex_index ); 00075 delta[dim] = tmp_var; 00076 pd.set_vertex_coordinates( delta, vertex_index, err ); 00077 // if step was too big, shorten it and go again 00078 if( feasible ) 00079 return eps; 00080 else 00081 eps *= rho; 00082 } // end while looking for feasible eps 00083 return 0.0; 00084 } // end function get_eps 00085 00086 bool ObjectiveFunction::compute_subpatch_numerical_gradient( EvalType type, EvalType subtype, PatchData& pd, 00087 double& flocal, Vector3D& grad, MsqError& err ) 00088 { 00089 assert( pd.num_free_vertices() == 1 ); 00090 00091 double flocald = 0; 00092 double eps = 0; 00093 00094 bool b = evaluate( type, pd, flocal, OF_FREE_EVALS_ONLY, err ); 00095 if( MSQ_CHKERR( err ) || !b ) { return false; } 00096 00097 // loop over the three coords x,y,z 00098 for( int j = 0; j < 3; ++j ) 00099 { 00100 eps = get_eps( pd, subtype, flocald, j, 0, err ); 00101 MSQ_ERRZERO( err ); 00102 if( eps == 0 ) 00103 { 00104 MSQ_SETERR( err ) 00105 ( "Dividing by zero in Objective Functions numerical grad", MsqError::INVALID_STATE ); 00106 return false; 00107 } 00108 grad[j] = ( flocald - flocal ) / eps; 00109 } 00110 return true; 00111 } 00112 00113 bool ObjectiveFunction::compute_patch_numerical_gradient( EvalType type, EvalType subtype, PatchData& pd, 00114 double& flocal, std::vector< Vector3D >& grad, MsqError& err ) 00115 { 00116 double flocald = 0; 00117 double eps = 0; 00118 00119 bool b = evaluate( type, pd, flocal, OF_FREE_EVALS_ONLY, err ); 00120 if( MSQ_CHKERR( err ) || !b ) { return false; } 00121 00122 for( size_t i = 0; i < pd.num_free_vertices(); ++i ) 00123 { 00124 // loop over the three coords x,y,z 00125 for( int j = 0; j < 3; ++j ) 00126 { 00127 eps = get_eps( pd, subtype, flocald, j, i, err ); 00128 MSQ_ERRZERO( err ); 00129 if( eps == 0 ) 00130 { 00131 MSQ_SETERR( err ) 00132 ( "Dividing by zero in Objective Functions numerical grad", MsqError::INVALID_STATE ); 00133 return false; 00134 } 00135 grad[i][j] = ( flocald - flocal ) / eps; 00136 } 00137 } 00138 00139 return true; 00140 } 00141 00142 /*! 00143 Numerically Calculates the gradient of the ObjectiveFunction for the 00144 free vertices in the patch. Returns 'false' if the patch is outside 00145 of a required feasible region, returns 'ture' otherwise. 00146 The behavior of the function depends on the value of the boolean 00147 useLocalGradient. If useLocalGradient is set to 00148 'true', compute_numerical_gradient creates a sub-patch around a free 00149 vertex, and then perturbs that vertex in one of the coordinate directions. 00150 Only the ObjectiveFunction value on the local sub-patch is used in the 00151 computation of the gradient. Therefore, useLocalGradient should only 00152 be set to 'true' for ObjectiveFunctions which can use this method. Unless 00153 the concrete ObjectiveFunction sets useLocalGradient to 'true' in its 00154 constructor, the value will be 'false'. In this case, the objective 00155 function value for the entire patch is used in the calculation of the 00156 gradient. This is computationally expensive, but it is numerically 00157 correct for all (C_1) functions. 00158 \param pd PatchData on which the gradient is taken. 00159 \param grad Array of Vector3D of length the number of vertices used to store gradient. 00160 \param OF_val will be set to the objective function value. 00161 */ 00162 bool ObjectiveFunction::evaluate_with_gradient( EvalType eval_type, PatchData& pd, double& OF_val, 00163 std::vector< Vector3D >& grad, MsqError& err ) 00164 { 00165 bool b; 00166 grad.resize( pd.num_free_vertices() ); 00167 00168 // Fast path for single-free-vertex patch 00169 if( pd.num_free_vertices() == 1 ) 00170 { 00171 const EvalType sub_type = ( eval_type == CALCULATE ) ? CALCULATE : TEMPORARY; 00172 b = compute_subpatch_numerical_gradient( eval_type, sub_type, pd, OF_val, grad[0], err ); 00173 return !MSQ_CHKERR( err ) && b; 00174 } 00175 00176 ObjectiveFunction* of = this; 00177 std::unique_ptr< ObjectiveFunction > deleter; 00178 if( eval_type == CALCULATE ) 00179 { 00180 of->clear(); 00181 b = of->evaluate( ACCUMULATE, pd, OF_val, OF_FREE_EVALS_ONLY, err ); 00182 if( err ) 00183 { // OF doesn't support BCD type evals, try slow method 00184 err.clear(); 00185 of->clear(); 00186 b = compute_patch_numerical_gradient( CALCULATE, CALCULATE, pd, OF_val, grad, err ); 00187 return !MSQ_CHKERR( err ) && b; 00188 } 00189 else if( !b ) 00190 return b; 00191 } 00192 else 00193 { 00194 b = this->evaluate( eval_type, pd, OF_val, OF_FREE_EVALS_ONLY, err ); 00195 if( MSQ_CHKERR( err ) || !b ) return false; 00196 of = this->clone(); 00197 deleter = std::unique_ptr< ObjectiveFunction >( of ); 00198 } 00199 00200 // Determine number of layers of adjacent elements based on metric type. 00201 unsigned layers = min_patch_layers(); 00202 00203 // Create a subpatch for each free vertex and use it to evaluate the 00204 // gradient for that vertex. 00205 double flocal; 00206 PatchData subpatch; 00207 for( size_t i = 0; i < pd.num_free_vertices(); ++i ) 00208 { 00209 pd.get_subpatch( i, layers, subpatch, err ); 00210 MSQ_ERRZERO( err ); 00211 b = of->compute_subpatch_numerical_gradient( SAVE, TEMPORARY, subpatch, flocal, grad[i], err ); 00212 if( MSQ_CHKERR( err ) || !b ) 00213 { 00214 of->clear(); 00215 return false; 00216 } 00217 } 00218 00219 of->clear(); 00220 return true; 00221 } 00222 00223 bool ObjectiveFunction::evaluate_with_Hessian_diagonal( EvalType type, PatchData& pd, double& value_out, 00224 std::vector< Vector3D >& grad_out, 00225 std::vector< SymMatrix3D >& hess_diag_out, MsqError& err ) 00226 { 00227 MsqHessian hess; 00228 hess.initialize( pd, err ); 00229 MSQ_ERRZERO( err ); 00230 bool val = evaluate_with_Hessian( type, pd, value_out, grad_out, hess, err ); 00231 MSQ_ERRZERO( err ); 00232 hess_diag_out.resize( hess.size() ); 00233 for( size_t i = 0; i < hess.size(); ++i ) 00234 hess_diag_out[i] = hess.get_block( i, i )->upper(); 00235 return val; 00236 } 00237 00238 bool ObjectiveFunction::evaluate_with_Hessian( EvalType /*type*/, PatchData& /*pd*/, double& /*value_out*/, 00239 std::vector< Vector3D >& /*grad_out*/, MsqHessian& /*Hessian_out*/, 00240 MsqError& err ) 00241 { 00242 MSQ_SETERR( err ) 00243 ( "No Hessian available for this objective function.\n" 00244 "Choose either a different objective function or a " 00245 "different solver.\n", 00246 MsqError::INVALID_STATE ); 00247 return false; 00248 } 00249 00250 } // namespace MBMesquite