MOAB: Mesh Oriented datABase  (version 5.2.1)
TInverseMeanRatio.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2007 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 TInverseMeanRatio.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "TInverseMeanRatio.hpp"
00034 #include "MsqMatrix.hpp"
00035 #include "MsqError.hpp"
00036 
00037 namespace MBMesquite
00038 {
00039 
00040 std::string TInverseMeanRatio::get_name() const
00041 {
00042     return "TInverseMeanRatio";
00043 }
00044 
00045 TInverseMeanRatio::~TInverseMeanRatio() {}
00046 
00047 bool TInverseMeanRatio::evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& err )
00048 {
00049     const double d = det( T );
00050     if( invalid_determinant( d ) )
00051     {
00052         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00053         return false;
00054     }
00055     else
00056     {
00057         result = sqr_Frobenius( T ) / ( 2 * d ) - 1;
00058         return true;
00059     }
00060 }
00061 
00062 bool TInverseMeanRatio::evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T,
00063                                             MsqError& err )
00064 {
00065     const double d = det( T );
00066     if( invalid_determinant( d ) )
00067     {
00068         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00069         deriv_wrt_T = MsqMatrix< 2, 2 >( 0.0 );
00070         return false;
00071     }
00072     else
00073     {
00074         result      = sqr_Frobenius( T ) / ( 2 * d );
00075         deriv_wrt_T = transpose_adj( T );
00076         deriv_wrt_T *= -result;
00077         deriv_wrt_T += T;
00078         deriv_wrt_T *= 1.0 / d;
00079         result -= 1.0;
00080         return true;
00081     }
00082 }
00083 
00084 bool TInverseMeanRatio::evaluate_with_hess( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& dA,
00085                                             MsqMatrix< 2, 2 > d2A[3], MsqError& err )
00086 {
00087     const double d = det( T );
00088     if( invalid_determinant( d ) )
00089     {
00090         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00091         dA = d2A[0] = d2A[1] = d2A[2] = MsqMatrix< 2, 2 >( 0.0 );
00092         return false;
00093     }
00094     else
00095     {
00096         const double inv_det = 1.0 / d;
00097         result               = sqr_Frobenius( T ) * 0.5 * inv_det;
00098 
00099         const MsqMatrix< 2, 2 > AT = transpose_adj( T );
00100         dA                         = AT;
00101         dA *= -result;
00102         dA += T;
00103         dA *= inv_det;
00104 
00105         const double p3                    = -result * inv_det;
00106         const double p1                    = -2.0 * p3 * inv_det;
00107         const double p2                    = -inv_det * inv_det;
00108         const MsqMatrix< 2, 2 > AT_T_op_00 = outer( AT.row( 0 ), T.row( 0 ) );
00109         const MsqMatrix< 2, 2 > AT_T_op_11 = outer( AT.row( 1 ), T.row( 1 ) );
00110         d2A[0] = p1 * outer( AT.row( 0 ), AT.row( 0 ) ) + p2 * ( AT_T_op_00 + transpose( AT_T_op_00 ) );
00111         d2A[1] = p1 * outer( AT.row( 0 ), AT.row( 1 ) ) +
00112                  p2 * ( outer( AT.row( 0 ), T.row( 1 ) ) + outer( T.row( 0 ), AT.row( 1 ) ) );
00113         d2A[2] = p1 * outer( AT.row( 1 ), AT.row( 1 ) ) + p2 * ( AT_T_op_11 + transpose( AT_T_op_11 ) );
00114 
00115         d2A[0]( 0, 0 ) += inv_det;
00116         d2A[0]( 1, 1 ) += inv_det;
00117         d2A[1]( 0, 1 ) += p3;
00118         d2A[1]( 1, 0 ) -= p3;
00119         d2A[2]( 0, 0 ) += inv_det;
00120         d2A[2]( 1, 1 ) += inv_det;
00121 
00122         result -= 1.0;
00123         return true;
00124     }
00125 }
00126 
00127 bool TInverseMeanRatio::evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& err )
00128 {
00129     const double d = det( T );
00130     if( invalid_determinant( d ) )
00131     {
00132         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00133         return false;
00134     }
00135     else
00136     {
00137         const double det_cbrt = MBMesquite::cbrt( d );
00138         result                = sqr_Frobenius( T ) / ( 3 * det_cbrt * det_cbrt ) - 1;
00139         return true;
00140     }
00141 }
00142 
00143 bool TInverseMeanRatio::evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& deriv_wrt_T,
00144                                             MsqError& err )
00145 {
00146     const double d = det( T );
00147     if( invalid_determinant( d ) )
00148     {
00149         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00150         deriv_wrt_T = MsqMatrix< 3, 3 >( 0.0 );
00151         return false;
00152     }
00153 
00154     const double inv_det             = 1.0 / d;
00155     const double inv_det_cbrt        = MBMesquite::cbrt( inv_det );
00156     const double inv_3_det_twothirds = inv_det_cbrt * inv_det_cbrt / 3.0;
00157     const double fnorm               = sqr_Frobenius( T );
00158     result                           = fnorm * inv_3_det_twothirds - 1;
00159     deriv_wrt_T                      = transpose_adj( T );
00160     deriv_wrt_T *= -fnorm * inv_det / 3.0;
00161     deriv_wrt_T += T;
00162     deriv_wrt_T *= 2.0 * inv_3_det_twothirds;
00163     return true;
00164 }
00165 
00166 bool TInverseMeanRatio::evaluate_with_hess( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& dA,
00167                                             MsqMatrix< 3, 3 > d2A[6], MsqError& err )
00168 {
00169     const double d = det( T );
00170     if( invalid_determinant( d ) )
00171     {
00172         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00173         dA = MsqMatrix< 3, 3 >( 0.0 );
00174         return false;
00175     }
00176 
00177     const double f0 = 1.0 / d;
00178     const double c  = MBMesquite::cbrt( f0 );
00179     const double f1 = 1.0 / 3.0 * c * c;
00180     const double f2 = sqr_Frobenius( T );
00181     result          = f1 * f2;
00182 
00183     const double f3 = 2 * f1;
00184     const double f4 = result * ( 10.0 / 9.0 ) * f0 * f0;
00185     const double f5 = ( 1.0 / 3.0 ) * f0 * f3;
00186     const double f6 = 2 * f5;
00187     const double f7 = f2 * f5;
00188 
00189     const MsqMatrix< 3, 3 > AT = transpose_adj( T );
00190     dA                         = AT;
00191     dA *= ( -1.0 / 3.0 ) * f0 * f2;
00192     dA += T;
00193     dA *= f3;
00194 
00195     MsqMatrix< 3, 3 > op;
00196     int i    = 0;
00197     double s = 1;
00198     for( int r = 0; r < 3; ++r )
00199     {
00200         d2A[i] = outer( AT.row( r ), AT.row( r ) );
00201         d2A[i] *= f4;
00202         op = outer( AT.row( r ), T.row( r ) );
00203         op += transpose( op );
00204         op *= f6;
00205         d2A[i] -= op;
00206 
00207         d2A[i]( 0, 0 ) += f3;
00208         d2A[i]( 1, 1 ) += f3;
00209         d2A[i]( 2, 2 ) += f3;
00210 
00211         ++i;
00212 
00213         for( int cc = r + 1; cc < 3; ++cc )
00214         {
00215             d2A[i] = outer( AT.row( r ), AT.row( cc ) );
00216             d2A[i] *= f4;
00217             op = outer( AT.row( r ), T.row( cc ) );
00218             op += outer( T.row( r ), AT.row( cc ) );
00219             op *= f6;
00220             d2A[i] -= op;
00221 
00222             MsqMatrix< 1, 3 > rt = T.row( 3 - r - cc );
00223             rt *= s * f7;
00224             d2A[i]( 0, 1 ) -= rt( 0, 2 );
00225             d2A[i]( 0, 2 ) += rt( 0, 1 );
00226             d2A[i]( 1, 0 ) += rt( 0, 2 );
00227             d2A[i]( 1, 2 ) -= rt( 0, 0 );
00228             d2A[i]( 2, 0 ) -= rt( 0, 1 );
00229             d2A[i]( 2, 1 ) += rt( 0, 0 );
00230 
00231             ++i;
00232             s = -s;
00233         }
00234     }
00235 
00236     result -= 1.0;
00237     return true;
00238 }
00239 
00240 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines