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