MOAB: Mesh Oriented datABase  (version 5.2.1)
TShapeSizeB1.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 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     (2006) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file TShapeSizeB1.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "TShapeSizeB1.hpp"
00034 #include "TMPDerivs.hpp"
00035 #include "MsqError.hpp"
00036 
00037 namespace MBMesquite
00038 {
00039 
00040 std::string TShapeSizeB1::get_name() const
00041 {
00042     return "TShapeSizeB1";
00043 }
00044 
00045 TShapeSizeB1::~TShapeSizeB1() {}
00046 
00047 bool TShapeSizeB1::evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& /*err*/ )
00048 {
00049     const double tau = det( T );
00050     if( invalid_determinant( tau ) )
00051     {  // barrier
00052         return false;
00053     }
00054 
00055     const double nT = sqr_Frobenius( T );
00056     const double f  = 1 / ( tau * tau );
00057     result          = ( 1 + f ) * nT - 4;
00058     return true;
00059 }
00060 
00061 bool TShapeSizeB1::evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& /*err*/ )
00062 {
00063     const double tau = det( T );
00064     if( invalid_determinant( tau ) )
00065     {  // barrier
00066         return false;
00067     }
00068 
00069     const double nT   = sqr_Frobenius( T );
00070     const double nadj = sqr_Frobenius( transpose_adj( T ) );
00071     const double f    = 1 / ( tau * tau );
00072     result            = nT + f * nadj - 6;
00073     return true;
00074 }
00075 
00076 bool TShapeSizeB1::evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T,
00077                                        MsqError& err )
00078 {
00079     const double tau = det( T );
00080     if( invalid_determinant( tau ) )
00081     {  // barrier
00082         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00083         return false;
00084     }
00085 
00086     const MsqMatrix< 2, 2 > adjt = transpose_adj( T );
00087     const double nT              = sqr_Frobenius( T );
00088     const double f               = 1 / ( tau * tau );
00089     result                       = ( 1 + f ) * nT - 4;
00090 
00091     deriv_wrt_T = T;
00092     deriv_wrt_T *= 2 + 2 * f;
00093     deriv_wrt_T -= 2 * f / tau * nT * adjt;
00094 
00095     return true;
00096 }
00097 
00098 bool TShapeSizeB1::evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& deriv_wrt_T,
00099                                        MsqError& err )
00100 {
00101     const double tau = det( T );
00102     if( invalid_determinant( tau ) )
00103     {  // barrier
00104         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00105         return false;
00106     }
00107 
00108     const MsqMatrix< 3, 3 > adjt = transpose_adj( T );
00109     const double nT              = sqr_Frobenius( T );
00110     const double nadj            = sqr_Frobenius( adjt );
00111     const double f               = 1 / ( tau * tau );
00112     result                       = nT + f * nadj - 6;
00113 
00114     deriv_wrt_T = T;
00115     deriv_wrt_T *= ( 1 + f * nT );
00116     deriv_wrt_T -= f * T * transpose( T ) * T;
00117     deriv_wrt_T -= f / tau * nadj * adjt;
00118     deriv_wrt_T *= 2;
00119 
00120     return true;
00121 }
00122 
00123 bool TShapeSizeB1::evaluate_with_hess( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T,
00124                                        MsqMatrix< 2, 2 > second_wrt_T[3], MsqError& err )
00125 {
00126     const double tau = det( T );
00127     if( invalid_determinant( tau ) )
00128     {  // barrier
00129         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00130         return false;
00131     }
00132 
00133     const MsqMatrix< 2, 2 > adjt = transpose_adj( T );
00134     const double nT              = sqr_Frobenius( T );
00135     const double f               = 1 / ( tau * tau );
00136     result                       = ( 1 + f ) * nT - 4;
00137 
00138     deriv_wrt_T = T;
00139     deriv_wrt_T *= 2 + 2 * f;
00140     deriv_wrt_T -= 2 * f / tau * nT * adjt;
00141 
00142     set_scaled_sum_outer_product( second_wrt_T, -4 * f / tau, T, adjt );
00143     pluseq_scaled_I( second_wrt_T, 2 + 2 * f );
00144     pluseq_scaled_outer_product( second_wrt_T, 6 * nT * f * f, adjt );
00145     pluseq_scaled_2nd_deriv_of_det( second_wrt_T, -2 * nT * f / tau );
00146 
00147     return true;
00148 }
00149 bool TShapeSizeB1::evaluate_with_hess( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& deriv_wrt_T,
00150                                        MsqMatrix< 3, 3 > second_wrt_T[6], MsqError& err )
00151 {
00152     const double tau = det( T );
00153     if( invalid_determinant( tau ) )
00154     {  // barrier
00155         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00156         return false;
00157     }
00158 
00159     const MsqMatrix< 3, 3 > adjt = transpose_adj( T );
00160     const double nT              = sqr_Frobenius( T );
00161     const double nadj            = sqr_Frobenius( adjt );
00162     const double f               = 1 / ( tau * tau );
00163     result                       = nT + f * nadj - 6;
00164 
00165     //! \f$ \frac{\partial}{\partial T} |adj T|^2 \f$
00166     const MsqMatrix< 3, 3 > dNadj_dT = 2 * ( nT * T - T * transpose( T ) * T );
00167     deriv_wrt_T                      = T;
00168     deriv_wrt_T -= f / tau * nadj * adjt;
00169     deriv_wrt_T *= 2;
00170     deriv_wrt_T += f * dNadj_dT;
00171 
00172     // calculate negative of 2nd wrt T of (|adj T|^2 / tau^2) (sec 3.2.2)
00173     set_scaled_2nd_deriv_norm_sqr_adj( second_wrt_T, f, T );
00174     pluseq_scaled_2nd_deriv_of_det( second_wrt_T, -2 * f * f * nadj * tau, T );
00175     pluseq_scaled_outer_product( second_wrt_T, 6 * f * f * nadj, adjt );
00176     pluseq_scaled_sum_outer_product( second_wrt_T, -2 * f * f * tau, adjt, dNadj_dT );
00177     // calculate 2nd wrt T of this metric
00178     pluseq_scaled_I( second_wrt_T, 2.0 );
00179 
00180     return true;
00181 }
00182 
00183 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines