MOAB: Mesh Oriented datABase  (version 5.2.1)
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, 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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines