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 (2009) kraftche@cae.wisc.edu 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, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T, 00063 MsqError& err ) 00064 { 00065 double tau = det( T ); 00066 if( TMetric::invalid_determinant( tau ) ) 00067 { 00068 MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED ); 00069 return false; 00070 } 00071 const double b = 0.5 / tau; 00072 00073 const double tr = trace( T ); 00074 const double f = 0.5 * fabs( tr ); 00075 result = sqr_Frobenius( T ) - f * tr; 00076 00077 deriv_wrt_T = T; 00078 pluseq_scaled_I( deriv_wrt_T, -f ); 00079 deriv_wrt_T *= 2 * tau; 00080 deriv_wrt_T -= result * transpose_adj( T ); 00081 00082 result *= b; 00083 deriv_wrt_T *= b / tau; 00084 00085 return true; 00086 } 00087 00088 bool TShapeOrientB2::evaluate_with_hess( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T, 00089 MsqMatrix< 2, 2 > second_wrt_T[3], MsqError& err ) 00090 { 00091 double tau = det( T ); 00092 if( TMetric::invalid_determinant( tau ) ) 00093 { 00094 MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED ); 00095 return false; 00096 } 00097 const double b = 0.5 / tau; 00098 00099 // calculate non-barrier value (ShapeOrientAlt1) 00100 const double tr = trace( T ); 00101 const double f = 0.5 * fabs( tr ); 00102 result = sqr_Frobenius( T ) - f * tr; 00103 00104 // calculate non-barrier first derivatives 00105 deriv_wrt_T = T; 00106 pluseq_scaled_I( deriv_wrt_T, -f ); 00107 deriv_wrt_T *= 2; 00108 00109 // calculate barrier second derivs 00110 const MsqMatrix< 2, 2 > adjt = transpose_adj( T ); 00111 set_scaled_sum_outer_product( second_wrt_T, -b / tau, deriv_wrt_T, adjt ); 00112 pluseq_scaled_outer_product( second_wrt_T, result / ( tau * tau * tau ), adjt ); 00113 pluseq_scaled_2nd_deriv_of_det( second_wrt_T, -result * b / tau, T ); 00114 // calculate non-barrier barrier portion of second derivs 00115 pluseq_scaled_I( second_wrt_T, 1 / tau ); 00116 pluseq_scaled_outer_product_I_I( second_wrt_T, 0.5 / tau * ( tr < 0 ? 1 : -1 ) ); 00117 00118 // calculate barrier derivs from non-barrier 00119 deriv_wrt_T *= tau; 00120 deriv_wrt_T -= result * adjt; 00121 deriv_wrt_T *= b / tau; 00122 00123 // barrier value from non-barrier 00124 result *= b; 00125 return true; 00126 } 00127 00128 bool TShapeOrientB2::evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& err ) 00129 { 00130 double tau = det( T ); 00131 if( TMetric::invalid_determinant( tau ) ) 00132 { 00133 MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED ); 00134 return false; 00135 } 00136 const double tr = trace( T ); 00137 result = ( 0.5 / tau ) * ( sqr_Frobenius( T ) - MSQ_ONE_THIRD * tr * fabs( tr ) ); 00138 return true; 00139 } 00140 00141 bool TShapeOrientB2::evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& deriv_wrt_T, 00142 MsqError& err ) 00143 { 00144 double tau = det( T ); 00145 if( TMetric::invalid_determinant( tau ) ) 00146 { 00147 MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED ); 00148 return false; 00149 } 00150 const double b = 0.5 / tau; 00151 00152 const double tr = trace( T ); 00153 const double f = MSQ_ONE_THIRD * fabs( tr ); 00154 result = sqr_Frobenius( T ) - f * tr; 00155 00156 deriv_wrt_T = T; 00157 pluseq_scaled_I( deriv_wrt_T, -f ); 00158 deriv_wrt_T *= 2 * tau; 00159 deriv_wrt_T -= result * transpose_adj( T ); 00160 00161 result *= b; 00162 deriv_wrt_T *= b / tau; 00163 00164 return true; 00165 } 00166 00167 bool TShapeOrientB2::evaluate_with_hess( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& deriv_wrt_T, 00168 MsqMatrix< 3, 3 > second_wrt_T[6], MsqError& err ) 00169 { 00170 double tau = det( T ); 00171 if( TMetric::invalid_determinant( tau ) ) 00172 { 00173 MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED ); 00174 result = 0.0; 00175 return false; 00176 } 00177 const double b = 0.5 / tau; 00178 00179 // calculate non-barrier value (ShapeOrientAlt1) 00180 const double tr = trace( T ); 00181 const double f = MSQ_ONE_THIRD * fabs( tr ); 00182 result = sqr_Frobenius( T ) - f * tr; 00183 00184 // calculate non-barrier first derivatives 00185 deriv_wrt_T = T; 00186 pluseq_scaled_I( deriv_wrt_T, -f ); 00187 deriv_wrt_T *= 2; 00188 00189 // calculate barrier second derivs 00190 const MsqMatrix< 3, 3 > adjt = transpose_adj( T ); 00191 set_scaled_sum_outer_product( second_wrt_T, -b / tau, deriv_wrt_T, adjt ); 00192 pluseq_scaled_outer_product( second_wrt_T, result / ( tau * tau * tau ), adjt ); 00193 pluseq_scaled_2nd_deriv_of_det( second_wrt_T, -result * b / tau, T ); 00194 // calculate non-barrier barrier portion of second derivs 00195 pluseq_scaled_I( second_wrt_T, 1 / tau ); 00196 pluseq_scaled_outer_product_I_I( second_wrt_T, MSQ_ONE_THIRD / tau * ( tr < 0 ? 1 : -1 ) ); 00197 00198 // calculate barrier derivs from non-barrier 00199 deriv_wrt_T *= tau; 00200 deriv_wrt_T -= result * adjt; 00201 deriv_wrt_T *= b / tau; 00202 00203 // barrier value from non-barrier 00204 result *= b; 00205 return true; 00206 } 00207 00208 } // namespace MBMesquite