MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2010 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 (2010) kraftche@cae.wisc.edu 00024 00025 ***************************************************************** */ 00026 00027 /** \file TMetric.cpp 00028 * \brief 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #include "Mesquite.hpp" 00033 #include "TMetric.hpp" 00034 #include "TMetricBarrier.hpp" 00035 #include "MsqMatrix.hpp" 00036 #include "MsqError.hpp" 00037 #include <limits> 00038 00039 namespace MBMesquite 00040 { 00041 00042 template < unsigned Dim > 00043 static inline double do_finite_difference( int r, int c, TMetric* metric, MsqMatrix< Dim, Dim > A, double value, 00044 MsqError& err ) 00045 { 00046 const double INITIAL_STEP = std::max( 1e-6, fabs( 1e-14 * value ) ); 00047 const double init = A( r, c ); 00048 bool valid; 00049 double diff_value; 00050 for( double step = INITIAL_STEP; step > std::numeric_limits< double >::epsilon(); step *= 0.1 ) 00051 { 00052 A( r, c ) = init + step; 00053 valid = metric->evaluate( A, diff_value, err ); 00054 MSQ_ERRZERO( err ); 00055 if( valid ) return ( diff_value - value ) / step; 00056 } 00057 00058 // If we couldn't find a valid step, try stepping in the other 00059 // direciton 00060 for( double step = INITIAL_STEP; step > std::numeric_limits< double >::epsilon(); step *= 0.1 ) 00061 { 00062 A( r, c ) = init - step; 00063 valid = metric->evaluate( A, diff_value, err ); 00064 MSQ_ERRZERO( err ); 00065 if( valid ) return ( value - diff_value ) / step; 00066 } 00067 // If that didn't work either, then give up. 00068 MSQ_SETERR( err ) 00069 ( "No valid step size for finite difference of 2D target metric.", MsqError::INTERNAL_ERROR ); 00070 return 0.0; 00071 } 00072 00073 template < unsigned Dim > 00074 static inline bool do_numerical_gradient( TMetric* mu, MsqMatrix< Dim, Dim > A, double& result, 00075 MsqMatrix< Dim, Dim >& wrt_A, MsqError& err ) 00076 { 00077 bool valid; 00078 valid = mu->evaluate( A, result, err ); 00079 MSQ_ERRZERO( err ); 00080 if( MSQ_CHKERR( err ) || !valid ) return valid; 00081 00082 switch( Dim ) 00083 { 00084 case 3: 00085 wrt_A( 0, 2 ) = do_finite_difference( 0, 2, mu, A, result, err ); 00086 MSQ_ERRZERO( err ); 00087 wrt_A( 1, 2 ) = do_finite_difference( 1, 2, mu, A, result, err ); 00088 MSQ_ERRZERO( err ); 00089 wrt_A( 2, 0 ) = do_finite_difference( 2, 0, mu, A, result, err ); 00090 MSQ_ERRZERO( err ); 00091 wrt_A( 2, 1 ) = do_finite_difference( 2, 1, mu, A, result, err ); 00092 MSQ_ERRZERO( err ); 00093 wrt_A( 2, 2 ) = do_finite_difference( 2, 2, mu, A, result, err ); 00094 MSQ_ERRZERO( err ); 00095 case 2: 00096 wrt_A( 0, 1 ) = do_finite_difference( 0, 1, mu, A, result, err ); 00097 MSQ_ERRZERO( err ); 00098 wrt_A( 1, 0 ) = do_finite_difference( 1, 0, mu, A, result, err ); 00099 MSQ_ERRZERO( err ); 00100 wrt_A( 1, 1 ) = do_finite_difference( 1, 1, mu, A, result, err ); 00101 MSQ_ERRZERO( err ); 00102 case 1: 00103 wrt_A( 0, 0 ) = do_finite_difference( 0, 0, mu, A, result, err ); 00104 MSQ_ERRZERO( err ); 00105 break; 00106 default: 00107 assert( false ); 00108 } 00109 return true; 00110 } 00111 00112 template < unsigned Dim > 00113 static inline bool do_numerical_hessian( TMetric* metric, MsqMatrix< Dim, Dim > A, double& value, 00114 MsqMatrix< Dim, Dim >& grad, MsqMatrix< Dim, Dim >* Hess, MsqError& err ) 00115 { 00116 // zero hessian data 00117 const int num_block = Dim * ( Dim + 1 ) / 2; 00118 for( int i = 0; i < num_block; ++i ) 00119 Hess[i].zero(); 00120 00121 // evaluate gradient for input values 00122 bool valid; 00123 valid = metric->evaluate_with_grad( A, value, grad, err ); 00124 if( MSQ_CHKERR( err ) || !valid ) return false; 00125 00126 // do finite difference for each term of A 00127 const double INITAL_STEP = std::max( 1e-6, fabs( 1e-14 * value ) ); 00128 double value2; 00129 MsqMatrix< Dim, Dim > grad2; 00130 for( unsigned r = 0; r < Dim; ++r ) 00131 { // for each row of A 00132 for( unsigned c = 0; c < Dim; ++c ) 00133 { // for each column of A 00134 const double in_val = A( r, c ); 00135 double step; 00136 for( step = INITAL_STEP; step > std::numeric_limits< double >::epsilon(); step *= 0.1 ) 00137 { 00138 A( r, c ) = in_val + step; 00139 valid = metric->evaluate_with_grad( A, value2, grad2, err ); 00140 MSQ_ERRZERO( err ); 00141 if( valid ) break; 00142 } 00143 00144 // if no valid step size, try step in other direction 00145 if( !valid ) 00146 { 00147 for( step = -INITAL_STEP; step < -std::numeric_limits< double >::epsilon(); step *= 0.1 ) 00148 { 00149 A( r, c ) = in_val + step; 00150 valid = metric->evaluate_with_grad( A, value2, grad2, err ); 00151 MSQ_ERRZERO( err ); 00152 if( valid ) break; 00153 } 00154 00155 // if still no valid step size, give up. 00156 if( !valid ) 00157 { 00158 MSQ_SETERR( err ) 00159 ( "No valid step size for finite difference of 2D target metric.", MsqError::INTERNAL_ERROR ); 00160 return false; 00161 } 00162 } 00163 00164 // restore A. 00165 A( r, c ) = in_val; 00166 00167 // add values into result matrix 00168 // values of grad2, in row-major order, are a single 9-value row of the Hessian 00169 grad2 -= grad; 00170 grad2 /= step; 00171 for( unsigned b = 0; b < r; ++b ) 00172 { 00173 const int idx = Dim * b - b * ( b + 1 ) / 2 + r; 00174 Hess[idx].add_column( c, transpose( grad2.row( b ) ) ); 00175 } 00176 for( unsigned b = r; b < Dim; ++b ) 00177 { 00178 const int idx = Dim * r - r * ( r + 1 ) / 2 + b; 00179 Hess[idx].add_row( c, grad2.row( b ) ); 00180 } 00181 } // for (c) 00182 } // for (r) 00183 00184 // Values in non-diagonal blocks were added twice. 00185 for( unsigned r = 0, h = 1; r < Dim - 1; ++r, ++h ) 00186 for( unsigned c = r + 1; c < Dim; ++c, ++h ) 00187 Hess[h] *= 0.5; 00188 00189 return true; 00190 } 00191 00192 TMetric::~TMetric() {} 00193 00194 bool TMetric::evaluate( const MsqMatrix< 2, 2 >& /*T*/, double& /*result*/, MsqError& /*err*/ ) 00195 { 00196 return false; 00197 } 00198 00199 bool TMetric::evaluate( const MsqMatrix< 3, 3 >& /*T*/, double& /*result*/, MsqError& /*err*/ ) 00200 { 00201 return false; 00202 } 00203 00204 bool TMetric::evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& wrt_T, MsqError& err ) 00205 { 00206 return do_numerical_gradient( this, T, result, wrt_T, err ); 00207 } 00208 00209 bool TMetric::evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& wrt_T, MsqError& err ) 00210 { 00211 return do_numerical_gradient( this, T, result, wrt_T, err ); 00212 } 00213 00214 bool TMetric::evaluate_with_hess( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T, 00215 MsqMatrix< 2, 2 > hess_wrt_T[3], MsqError& err ) 00216 { 00217 return do_numerical_hessian( this, T, result, deriv_wrt_T, hess_wrt_T, err ); 00218 } 00219 00220 bool TMetric::evaluate_with_hess( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& deriv_wrt_T, 00221 MsqMatrix< 3, 3 > hess_wrt_T[6], MsqError& err ) 00222 { 00223 return do_numerical_hessian( this, T, result, deriv_wrt_T, hess_wrt_T, err ); 00224 } 00225 00226 TMetric2D::~TMetric2D() {} 00227 TMetric3D::~TMetric3D() {} 00228 00229 bool TMetric2D::evaluate( const MsqMatrix< 3, 3 >&, double&, MsqError& err ) 00230 { 00231 MSQ_SETERR( err ) 00232 ( "2D target metric cannot be evaluated for volume elements", MsqError::UNSUPPORTED_ELEMENT ); 00233 return false; 00234 } 00235 00236 bool TMetric3D::evaluate( const MsqMatrix< 2, 2 >&, double&, MsqError& err ) 00237 { 00238 MSQ_SETERR( err ) 00239 ( "2D target metric cannot be evaluated for volume elements", MsqError::UNSUPPORTED_ELEMENT ); 00240 return false; 00241 } 00242 00243 } // namespace MBMesquite