MOAB: Mesh Oriented datABase  (version 5.4.0)
TMetric.cpp
Go to the documentation of this file.
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,
00044                                            int c,
00045                                            TMetric* metric,
00046                                            MsqMatrix< Dim, Dim > A,
00047                                            double value,
00048                                            MsqError& err )
00049 {
00050     const double INITIAL_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 = INITIAL_STEP; step > std::numeric_limits< double >::epsilon(); step *= 0.1 )
00055     {
00056         A( r, c ) = init + step;
00057         valid     = metric->evaluate( A, 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 = INITIAL_STEP; step > std::numeric_limits< double >::epsilon(); step *= 0.1 )
00065     {
00066         A( r, c ) = init - step;
00067         valid     = metric->evaluate( A, diff_value, err );
00068         MSQ_ERRZERO( err );
00069         if( valid ) return ( value - diff_value ) / step;
00070     }
00071     // If that didn't work either, then give up.
00072     MSQ_SETERR( err )
00073     ( "No valid step size for finite difference of 2D target metric.", MsqError::INTERNAL_ERROR );
00074     return 0.0;
00075 }
00076 
00077 template < unsigned Dim >
00078 static inline bool do_numerical_gradient( TMetric* mu,
00079                                           MsqMatrix< Dim, Dim > A,
00080                                           double& result,
00081                                           MsqMatrix< Dim, Dim >& wrt_A,
00082                                           MsqError& err )
00083 {
00084     bool valid;
00085     valid = mu->evaluate( A, result, err );
00086     MSQ_ERRZERO( err );
00087     if( MSQ_CHKERR( err ) || !valid ) return valid;
00088 
00089     switch( Dim )
00090     {
00091         case 3:
00092             wrt_A( 0, 2 ) = do_finite_difference( 0, 2, mu, A, result, err );
00093             MSQ_ERRZERO( err );
00094             wrt_A( 1, 2 ) = do_finite_difference( 1, 2, mu, A, result, err );
00095             MSQ_ERRZERO( err );
00096             wrt_A( 2, 0 ) = do_finite_difference( 2, 0, mu, A, result, err );
00097             MSQ_ERRZERO( err );
00098             wrt_A( 2, 1 ) = do_finite_difference( 2, 1, mu, A, result, err );
00099             MSQ_ERRZERO( err );
00100             wrt_A( 2, 2 ) = do_finite_difference( 2, 2, mu, A, result, err );
00101             MSQ_ERRZERO( err );
00102         case 2:
00103             wrt_A( 0, 1 ) = do_finite_difference( 0, 1, mu, A, result, err );
00104             MSQ_ERRZERO( err );
00105             wrt_A( 1, 0 ) = do_finite_difference( 1, 0, mu, A, result, err );
00106             MSQ_ERRZERO( err );
00107             wrt_A( 1, 1 ) = do_finite_difference( 1, 1, mu, A, result, err );
00108             MSQ_ERRZERO( err );
00109         case 1:
00110             wrt_A( 0, 0 ) = do_finite_difference( 0, 0, mu, A, result, err );
00111             MSQ_ERRZERO( err );
00112             break;
00113         default:
00114             assert( false );
00115     }
00116     return true;
00117 }
00118 
00119 template < unsigned Dim >
00120 static inline bool do_numerical_hessian( TMetric* metric,
00121                                          MsqMatrix< Dim, Dim > A,
00122                                          double& value,
00123                                          MsqMatrix< Dim, Dim >& grad,
00124                                          MsqMatrix< Dim, Dim >* Hess,
00125                                          MsqError& err )
00126 {
00127     // zero hessian data
00128     const int num_block = Dim * ( Dim + 1 ) / 2;
00129     for( int i = 0; i < num_block; ++i )
00130         Hess[i].zero();
00131 
00132     // evaluate gradient for input values
00133     bool valid;
00134     valid = metric->evaluate_with_grad( A, value, grad, err );
00135     if( MSQ_CHKERR( err ) || !valid ) return false;
00136 
00137     // do finite difference for each term of A
00138     const double INITAL_STEP = std::max( 1e-6, fabs( 1e-14 * value ) );
00139     double value2;
00140     MsqMatrix< Dim, Dim > grad2;
00141     for( unsigned r = 0; r < Dim; ++r )
00142     {  // for each row of A
00143         for( unsigned c = 0; c < Dim; ++c )
00144         {  // for each column of A
00145             const double in_val = A( r, c );
00146             double step;
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 no valid step size, try step in other direction
00156             if( !valid )
00157             {
00158                 for( step = -INITAL_STEP; step < -std::numeric_limits< double >::epsilon(); step *= 0.1 )
00159                 {
00160                     A( r, c ) = in_val + step;
00161                     valid     = metric->evaluate_with_grad( A, value2, grad2, err );
00162                     MSQ_ERRZERO( err );
00163                     if( valid ) break;
00164                 }
00165 
00166                 // if still no valid step size, give up.
00167                 if( !valid )
00168                 {
00169                     MSQ_SETERR( err )
00170                     ( "No valid step size for finite difference of 2D target metric.", MsqError::INTERNAL_ERROR );
00171                     return false;
00172                 }
00173             }
00174 
00175             // restore A.
00176             A( r, c ) = in_val;
00177 
00178             // add values into result matrix
00179             // values of grad2, in row-major order, are a single 9-value row of the Hessian
00180             grad2 -= grad;
00181             grad2 /= step;
00182             for( unsigned b = 0; b < r; ++b )
00183             {
00184                 const int idx = Dim * b - b * ( b + 1 ) / 2 + r;
00185                 Hess[idx].add_column( c, transpose( grad2.row( b ) ) );
00186             }
00187             for( unsigned b = r; b < Dim; ++b )
00188             {
00189                 const int idx = Dim * r - r * ( r + 1 ) / 2 + b;
00190                 Hess[idx].add_row( c, grad2.row( b ) );
00191             }
00192         }  // for (c)
00193     }      // for (r)
00194 
00195     // Values in non-diagonal blocks were added twice.
00196     for( unsigned r = 0, h = 1; r < Dim - 1; ++r, ++h )
00197         for( unsigned c = r + 1; c < Dim; ++c, ++h )
00198             Hess[h] *= 0.5;
00199 
00200     return true;
00201 }
00202 
00203 TMetric::~TMetric() {}
00204 
00205 bool TMetric::evaluate( const MsqMatrix< 2, 2 >& /*T*/, double& /*result*/, MsqError& /*err*/ )
00206 {
00207     return false;
00208 }
00209 
00210 bool TMetric::evaluate( const MsqMatrix< 3, 3 >& /*T*/, double& /*result*/, MsqError& /*err*/ )
00211 {
00212     return false;
00213 }
00214 
00215 bool TMetric::evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& wrt_T, MsqError& err )
00216 {
00217     return do_numerical_gradient( this, T, result, wrt_T, err );
00218 }
00219 
00220 bool TMetric::evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& wrt_T, MsqError& err )
00221 {
00222     return do_numerical_gradient( this, T, result, wrt_T, err );
00223 }
00224 
00225 bool TMetric::evaluate_with_hess( const MsqMatrix< 2, 2 >& T,
00226                                   double& result,
00227                                   MsqMatrix< 2, 2 >& deriv_wrt_T,
00228                                   MsqMatrix< 2, 2 > hess_wrt_T[3],
00229                                   MsqError& err )
00230 {
00231     return do_numerical_hessian( this, T, result, deriv_wrt_T, hess_wrt_T, err );
00232 }
00233 
00234 bool TMetric::evaluate_with_hess( const MsqMatrix< 3, 3 >& T,
00235                                   double& result,
00236                                   MsqMatrix< 3, 3 >& deriv_wrt_T,
00237                                   MsqMatrix< 3, 3 > hess_wrt_T[6],
00238                                   MsqError& err )
00239 {
00240     return do_numerical_hessian( this, T, result, deriv_wrt_T, hess_wrt_T, err );
00241 }
00242 
00243 TMetric2D::~TMetric2D() {}
00244 TMetric3D::~TMetric3D() {}
00245 
00246 bool TMetric2D::evaluate( const MsqMatrix< 3, 3 >&, double&, MsqError& err )
00247 {
00248     MSQ_SETERR( err )
00249     ( "2D target metric cannot be evaluated for volume elements", MsqError::UNSUPPORTED_ELEMENT );
00250     return false;
00251 }
00252 
00253 bool TMetric3D::evaluate( const MsqMatrix< 2, 2 >&, double&, MsqError& err )
00254 {
00255     MSQ_SETERR( err )
00256     ( "2D target metric cannot be evaluated for volume elements", MsqError::UNSUPPORTED_ELEMENT );
00257     return false;
00258 }
00259 
00260 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines