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