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