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