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 TShapeB1.cpp 00028 * \brief 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #include "Mesquite.hpp" 00033 #include "TShapeB1.hpp" 00034 #include "MsqMatrix.hpp" 00035 #include "MsqError.hpp" 00036 #include "TMPDerivs.hpp" 00037 00038 namespace MBMesquite 00039 { 00040 00041 std::string TShapeB1::get_name() const 00042 { 00043 return "TShapeB1"; 00044 } 00045 00046 TShapeB1::~TShapeB1() {} 00047 00048 bool TShapeB1::evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& err ) 00049 { 00050 const double d = det( T ); 00051 if( invalid_determinant( d ) ) 00052 { // barrier 00053 MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED ); 00054 return false; 00055 } 00056 00057 result = 0.5 * sqr_Frobenius( T ) / d - 1; 00058 return true; 00059 } 00060 00061 bool TShapeB1::evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T, 00062 MsqError& err ) 00063 { 00064 const double d = det( T ); 00065 if( invalid_determinant( d ) ) 00066 { // barrier 00067 MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED ); 00068 return false; 00069 } 00070 00071 double inv_d = 1.0 / d; 00072 result = 0.5 * sqr_Frobenius( T ) * inv_d; 00073 deriv_wrt_T = T; 00074 deriv_wrt_T -= result * transpose_adj( T ); 00075 deriv_wrt_T *= inv_d; 00076 00077 result -= 1.0; 00078 return true; 00079 } 00080 00081 /** \f$ \frac{\partial^2 \mu}{\partial T^2} 00082 = \frac{1}{\tau} I_4 00083 - \frac{1}{\tau^2} \left( T \otimes \frac{\partial \tau}{\partial T} 00084 + \frac{\partial \tau}{\partial T} \otimes T \right) 00085 + \frac{|T|^2}{\tau^3} \left( \frac{\partial \tau}{\partial T} \otimes 00086 \frac{\partial \tau}{\partial T} \right) 00087 - \frac{|T|^2}{2 \tau^3} \frac{\partial^2 \tau}{\partial T^2} \f$ 00088 */ 00089 bool TShapeB1::evaluate_with_hess( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T, 00090 MsqMatrix< 2, 2 > second_wrt_T[3], MsqError& err ) 00091 { 00092 const double d = det( T ); 00093 if( invalid_determinant( d ) ) 00094 { // barrier 00095 MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED ); 00096 return false; 00097 } 00098 00099 double inv_d = 1.0 / d; 00100 double f1 = sqr_Frobenius( T ) * inv_d; 00101 result = 0.5 * f1; 00102 const MsqMatrix< 2, 2 > adjt = transpose_adj( T ); 00103 deriv_wrt_T = T; 00104 deriv_wrt_T -= result * adjt; 00105 deriv_wrt_T *= inv_d; 00106 00107 set_scaled_outer_product( second_wrt_T, f1 * inv_d * inv_d, adjt ); 00108 pluseq_scaled_sum_outer_product( second_wrt_T, -inv_d * inv_d, T, adjt ); 00109 pluseq_scaled_I( second_wrt_T, inv_d ); 00110 pluseq_scaled_2nd_deriv_of_det( second_wrt_T, -result * inv_d ); 00111 00112 result -= 1.0; 00113 return true; 00114 } 00115 00116 bool TShapeB1::evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& err ) 00117 { 00118 double f = Frobenius( T ); 00119 double d = det( T ); 00120 double den = 3 * MSQ_SQRT_THREE * d; 00121 00122 if( invalid_determinant( d ) ) 00123 { 00124 MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED ); 00125 return false; 00126 } 00127 result = ( f * f * f ) / den - 1.0; 00128 return true; 00129 } 00130 00131 bool TShapeB1::evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& wrt_T, MsqError& err ) 00132 { 00133 double d = det( T ); 00134 if( invalid_determinant( d ) ) 00135 { 00136 MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED ); 00137 return false; 00138 } 00139 00140 double norm = Frobenius( T ); 00141 double den = 1.0 / ( 3 * MSQ_SQRT_THREE * d ); 00142 double norm_cube = norm * norm * norm; 00143 result = norm_cube * den - 1.0; 00144 wrt_T = T; 00145 wrt_T *= 3 * norm * den; 00146 wrt_T -= norm_cube * den / d * transpose_adj( T ); 00147 return true; 00148 } 00149 00150 bool TShapeB1::evaluate_with_hess( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& deriv_wrt_T, 00151 MsqMatrix< 3, 3 > second_wrt_T[6], MsqError& err ) 00152 { 00153 double d = det( T ); 00154 if( invalid_determinant( d ) ) 00155 { 00156 MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED ); 00157 return false; 00158 } 00159 00160 double id = 1.0 / d; 00161 double norm = Frobenius( T ); 00162 double den = 1.0 / ( 3 * MSQ_SQRT_THREE * d ); 00163 double norm_cube = norm * norm * norm; 00164 result = norm_cube * den - 1.0; 00165 MsqMatrix< 3, 3 > adjt = transpose_adj( T ); 00166 deriv_wrt_T = T; 00167 deriv_wrt_T *= 3 * norm * den; 00168 deriv_wrt_T -= norm_cube * den * id * transpose_adj( T ); 00169 00170 set_scaled_outer_product( second_wrt_T, 3 * den / norm, T ); 00171 pluseq_scaled_I( second_wrt_T, 3 * norm * den ); 00172 pluseq_scaled_2nd_deriv_of_det( second_wrt_T, -den * norm_cube * id, T ); 00173 pluseq_scaled_outer_product( second_wrt_T, 2 * den * norm_cube * id * id, adjt ); 00174 pluseq_scaled_sum_outer_product( second_wrt_T, -3 * norm * den * id, T, adjt ); 00175 00176 return true; 00177 } 00178 00179 } // namespace MBMesquite