MOAB: Mesh Oriented datABase
(version 5.2.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 PowerQualityMetric.cpp 00029 \brief 00030 00031 \author Michael Brewer 00032 \date April 1, 2003 00033 */ 00034 00035 #include "PowerQualityMetric.hpp" 00036 #include "MsqError.hpp" 00037 00038 using namespace MBMesquite; 00039 00040 PowerQualityMetric::PowerQualityMetric( QualityMetric* qm, double pow_factor ) : mMetric( *qm ), mPower( pow_factor ) {} 00041 00042 PowerQualityMetric::~PowerQualityMetric() {} 00043 00044 std::string PowerQualityMetric::get_name() const 00045 { 00046 return std::string( "pow(" ) + mMetric.get_name() + ")"; 00047 } 00048 00049 int PowerQualityMetric::get_negate_flag() const 00050 { 00051 return mPower.value() < 0 ? mMetric.get_negate_flag() : -mMetric.get_negate_flag(); 00052 } 00053 00054 void PowerQualityMetric::get_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free_only, MsqError& err ) 00055 { 00056 mMetric.get_evaluations( pd, handles, free_only, err );MSQ_CHKERR( err ); 00057 } 00058 00059 bool PowerQualityMetric::evaluate( PatchData& pd, size_t handle, double& value, MsqError& err ) 00060 { 00061 bool rval = mMetric.evaluate( pd, handle, value, err ); 00062 value = mPower.raise( value ); 00063 return !MSQ_CHKERR( err ) && rval; 00064 } 00065 00066 bool PowerQualityMetric::evaluate_with_indices( PatchData& pd, size_t handle, double& value, 00067 std::vector< size_t >& indices, MsqError& err ) 00068 { 00069 bool rval = mMetric.evaluate_with_indices( pd, handle, value, indices, err ); 00070 value = mPower.raise( value ); 00071 return !MSQ_CHKERR( err ) && rval; 00072 } 00073 00074 bool PowerQualityMetric::evaluate_with_gradient( PatchData& pd, size_t handle, double& value, 00075 std::vector< size_t >& indices, std::vector< Vector3D >& gradient, 00076 MsqError& err ) 00077 { 00078 bool rval = mMetric.evaluate_with_gradient( pd, handle, value, indices, gradient, err ); 00079 const double v = mPower.raise( value ); 00080 const double g = fabs( value ) > DBL_EPSILON ? mPower.value() * v / value : 0; 00081 value = v; 00082 for( std::vector< Vector3D >::iterator i = gradient.begin(); i != gradient.end(); ++i ) 00083 *i *= g; 00084 return !MSQ_CHKERR( err ) && rval; 00085 } 00086 00087 bool PowerQualityMetric::evaluate_with_Hessian( PatchData& pd, size_t handle, double& value, 00088 std::vector< size_t >& indices, std::vector< Vector3D >& gradient, 00089 std::vector< Matrix3D >& Hessian, MsqError& err ) 00090 { 00091 indices.clear(); 00092 bool rval = mMetric.evaluate_with_Hessian( pd, handle, value, indices, gradient, Hessian, err ); 00093 const double v = mPower.raise( value ); 00094 const double g = fabs( value ) > DBL_EPSILON ? mPower.value() * v / value : 0.0; 00095 const double h = fabs( value ) > DBL_EPSILON ? g * ( mPower.value() - 1 ) / value : 0.0; 00096 value = v; 00097 Matrix3D op; 00098 unsigned idx = 0; 00099 unsigned n = indices.size(); 00100 for( unsigned r = 0; r < n; ++r ) 00101 { 00102 for( unsigned c = r; c < n; ++c ) 00103 { 00104 Hessian[idx] *= g; 00105 op.outer_product( gradient[r], gradient[c] ); 00106 op *= h; 00107 Hessian[idx] += op; 00108 ++idx; 00109 } 00110 gradient[r] *= g; 00111 } 00112 return !MSQ_CHKERR( err ) && rval; 00113 }