MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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