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