MOAB: Mesh Oriented datABase  (version 5.2.1)
TShapeB1.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 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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines