MOAB: Mesh Oriented datABase  (version 5.2.1)
TShapeSizeB3.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines