MOAB: Mesh Oriented datABase  (version 5.4.1)
TShapeB1.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) [email protected]
00024 
00025   ***************************************************************** */
00026 
00027 /** \file TShapeB1.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "TShapeB1.hpp"
00034 #include "MsqMatrix.hpp"
00035 #include "MsqError.hpp"
00036 #include "TMPDerivs.hpp"
00037 
00038 namespace MBMesquite
00039 {
00040 
00041 std::string TShapeB1::get_name() const
00042 {
00043     return "TShapeB1";
00044 }
00045 
00046 TShapeB1::~TShapeB1() {}
00047 
00048 bool TShapeB1::evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& err )
00049 {
00050     const double d = det( T );
00051     if( invalid_determinant( d ) )
00052     {  // barrier
00053         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00054         return false;
00055     }
00056 
00057     result = 0.5 * sqr_Frobenius( T ) / d - 1;
00058     return true;
00059 }
00060 
00061 bool TShapeB1::evaluate_with_grad( const MsqMatrix< 2, 2 >& T,
00062                                    double& result,
00063                                    MsqMatrix< 2, 2 >& deriv_wrt_T,
00064                                    MsqError& err )
00065 {
00066     const double d = det( T );
00067     if( invalid_determinant( d ) )
00068     {  // barrier
00069         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00070         return false;
00071     }
00072 
00073     double inv_d = 1.0 / d;
00074     result       = 0.5 * sqr_Frobenius( T ) * inv_d;
00075     deriv_wrt_T  = T;
00076     deriv_wrt_T -= result * transpose_adj( T );
00077     deriv_wrt_T *= inv_d;
00078 
00079     result -= 1.0;
00080     return true;
00081 }
00082 
00083 /** \f$ \frac{\partial^2 \mu}{\partial T^2}
00084       = \frac{1}{\tau} I_4
00085       - \frac{1}{\tau^2} \left( T \otimes \frac{\partial \tau}{\partial T}
00086                           + \frac{\partial \tau}{\partial T} \otimes T \right)
00087       + \frac{|T|^2}{\tau^3} \left( \frac{\partial \tau}{\partial T} \otimes
00088                                \frac{\partial \tau}{\partial T} \right)
00089       - \frac{|T|^2}{2 \tau^3} \frac{\partial^2 \tau}{\partial T^2} \f$
00090   */
00091 bool TShapeB1::evaluate_with_hess( const MsqMatrix< 2, 2 >& T,
00092                                    double& result,
00093                                    MsqMatrix< 2, 2 >& deriv_wrt_T,
00094                                    MsqMatrix< 2, 2 > second_wrt_T[3],
00095                                    MsqError& err )
00096 {
00097     const double d = det( T );
00098     if( invalid_determinant( d ) )
00099     {  // barrier
00100         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00101         return false;
00102     }
00103 
00104     double inv_d                 = 1.0 / d;
00105     double f1                    = sqr_Frobenius( T ) * inv_d;
00106     result                       = 0.5 * f1;
00107     const MsqMatrix< 2, 2 > adjt = transpose_adj( T );
00108     deriv_wrt_T                  = T;
00109     deriv_wrt_T -= result * adjt;
00110     deriv_wrt_T *= inv_d;
00111 
00112     set_scaled_outer_product( second_wrt_T, f1 * inv_d * inv_d, adjt );
00113     pluseq_scaled_sum_outer_product( second_wrt_T, -inv_d * inv_d, T, adjt );
00114     pluseq_scaled_I( second_wrt_T, inv_d );
00115     pluseq_scaled_2nd_deriv_of_det( second_wrt_T, -result * inv_d );
00116 
00117     result -= 1.0;
00118     return true;
00119 }
00120 
00121 bool TShapeB1::evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& err )
00122 {
00123     double f   = Frobenius( T );
00124     double d   = det( T );
00125     double den = 3 * MSQ_SQRT_THREE * d;
00126 
00127     if( invalid_determinant( d ) )
00128     {
00129         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00130         return false;
00131     }
00132     result = ( f * f * f ) / den - 1.0;
00133     return true;
00134 }
00135 
00136 bool TShapeB1::evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& wrt_T, MsqError& err )
00137 {
00138     double d = det( T );
00139     if( invalid_determinant( d ) )
00140     {
00141         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00142         return false;
00143     }
00144 
00145     double norm      = Frobenius( T );
00146     double den       = 1.0 / ( 3 * MSQ_SQRT_THREE * d );
00147     double norm_cube = norm * norm * norm;
00148     result           = norm_cube * den - 1.0;
00149     wrt_T            = T;
00150     wrt_T *= 3 * norm * den;
00151     wrt_T -= norm_cube * den / d * transpose_adj( T );
00152     return true;
00153 }
00154 
00155 bool TShapeB1::evaluate_with_hess( const MsqMatrix< 3, 3 >& T,
00156                                    double& result,
00157                                    MsqMatrix< 3, 3 >& deriv_wrt_T,
00158                                    MsqMatrix< 3, 3 > second_wrt_T[6],
00159                                    MsqError& err )
00160 {
00161     double d = det( T );
00162     if( invalid_determinant( d ) )
00163     {
00164         MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
00165         return false;
00166     }
00167 
00168     double id              = 1.0 / d;
00169     double norm            = Frobenius( T );
00170     double den             = 1.0 / ( 3 * MSQ_SQRT_THREE * d );
00171     double norm_cube       = norm * norm * norm;
00172     result                 = norm_cube * den - 1.0;
00173     MsqMatrix< 3, 3 > adjt = transpose_adj( T );
00174     deriv_wrt_T            = T;
00175     deriv_wrt_T *= 3 * norm * den;
00176     deriv_wrt_T -= norm_cube * den * id * transpose_adj( T );
00177 
00178     set_scaled_outer_product( second_wrt_T, 3 * den / norm, T );
00179     pluseq_scaled_I( second_wrt_T, 3 * norm * den );
00180     pluseq_scaled_2nd_deriv_of_det( second_wrt_T, -den * norm_cube * id, T );
00181     pluseq_scaled_outer_product( second_wrt_T, 2 * den * norm_cube * id * id, adjt );
00182     pluseq_scaled_sum_outer_product( second_wrt_T, -3 * norm * den * id, T, adjt );
00183 
00184     return true;
00185 }
00186 
00187 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines