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 TShapeSize2DNB1.cpp 00028 * \brief 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #include "Mesquite.hpp" 00033 #include "TShapeSize2DNB1.hpp" 00034 #include "MsqMatrix.hpp" 00035 #include "TMPDerivs.hpp" 00036 00037 #include <iostream> 00038 00039 namespace MBMesquite 00040 { 00041 00042 std::string TShapeSize2DNB1::get_name() const 00043 { 00044 return "TShapeSize2DNB1"; 00045 } 00046 00047 TShapeSize2DNB1::~TShapeSize2DNB1() {} 00048 00049 /** \f$ \mu(T) = |T|^2 - 2 \psi(T) + 2 \f$ 00050 * \f$ \psi(T) = \sqrt{|T|^2 + 2 \tau} \f$ 00051 * \f$ \tau = det(T) \f$ 00052 */ 00053 bool TShapeSize2DNB1::evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& /*err*/ ) 00054 { 00055 double frob_sqr = sqr_Frobenius( T ); 00056 double psi = sqrt( frob_sqr + 2.0 * det( T ) ); 00057 00058 MsqMatrix< 2, 2 > Tdelta( T ); 00059 while( fabs( psi ) < DBL_EPSILON ) 00060 { 00061 Tdelta( 0, 0 ) += 1e-12; 00062 Tdelta( 1, 1 ) += 1e-12; 00063 frob_sqr = sqr_Frobenius( Tdelta ); 00064 psi = sqrt( frob_sqr + 2.0 * det( Tdelta ) ); 00065 } 00066 00067 result = frob_sqr - 2.0 * psi + 2.0; 00068 return true; 00069 } 00070 00071 bool TShapeSize2DNB1::evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T, 00072 MsqError& /*err*/ ) 00073 { 00074 double frob_sqr = sqr_Frobenius( T ); 00075 double psi = sqrt( frob_sqr + 2.0 * det( T ) ); 00076 00077 MsqMatrix< 2, 2 > Tdelta( T ); 00078 while( fabs( psi ) < DBL_EPSILON ) 00079 { 00080 Tdelta( 0, 0 ) += 1e-12; 00081 Tdelta( 1, 1 ) += 1e-12; 00082 frob_sqr = sqr_Frobenius( Tdelta ); 00083 psi = sqrt( frob_sqr + 2.0 * det( Tdelta ) ); 00084 } 00085 00086 result = frob_sqr - 2.0 * psi + 2.0; 00087 00088 deriv_wrt_T = T; 00089 if( psi > 1e-50 ) 00090 { 00091 deriv_wrt_T *= ( 1.0 - 1.0 / psi ); 00092 deriv_wrt_T -= 1.0 / psi * transpose_adj( T ); 00093 deriv_wrt_T *= 2; 00094 } 00095 else 00096 { 00097 std::cout << "Warning: Division by zero avoided in TShapeSize2DNB2::evaluate_with_grad()" << std::endl; 00098 } 00099 00100 return true; 00101 } 00102 00103 bool TShapeSize2DNB1::evaluate_with_hess( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T, 00104 MsqMatrix< 2, 2 > second[3], MsqError& /*err*/ ) 00105 { 00106 double frob_sqr = sqr_Frobenius( T ); 00107 double psi = sqrt( frob_sqr + 2.0 * det( T ) ); 00108 00109 MsqMatrix< 2, 2 > Tdelta( T ); 00110 while( fabs( psi ) < DBL_EPSILON ) 00111 { 00112 Tdelta( 0, 0 ) += 1e-12; 00113 Tdelta( 1, 1 ) += 1e-12; 00114 frob_sqr = sqr_Frobenius( Tdelta ); 00115 psi = sqrt( frob_sqr + 2.0 * det( Tdelta ) ); 00116 } 00117 00118 result = frob_sqr - 2.0 * psi + 2.0; 00119 00120 deriv_wrt_T = T; 00121 if( psi > 1e-50 ) 00122 { 00123 deriv_wrt_T *= ( 1.0 - 1.0 / psi ); 00124 deriv_wrt_T -= 1.0 / psi * transpose_adj( T ); 00125 deriv_wrt_T *= 2; 00126 } 00127 else 00128 { 00129 std::cout << "Warning: Division by zero avoided in TShapeSize2DNB2::evaluate_with_hess()" << std::endl; 00130 } 00131 00132 set_scaled_2nd_deriv_wrt_psi( second, -2.0, psi, T ); 00133 pluseq_scaled_I( second, 2 ); 00134 00135 return true; 00136 } 00137 00138 } // namespace MBMesquite