MOAB: Mesh Oriented datABase  (version 5.2.1)
TMPDerivs.hpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2009 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     (2009) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file TMPDerivs.hpp
00028  *  \brief Utility methods for calculating derivatives of TMP metrics
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #ifndef MSQ_TMP_DERIVS_HPP
00033 #define MSQ_TMP_DERIVS_HPP
00034 
00035 #include "Mesquite.hpp"
00036 #include "MsqMatrix.hpp"
00037 
00038 /** If defined, then \f$ (A \otimes B)_{i,j} = row(A,i)^t \otimes row(B,j) \f$
00039     otherwise \f$ (A \otimes B)_{i,j} = column(A,i) \otimes column(B,j)^t \f$
00040   */
00041 #define MSQ_ROW_BASED_OUTER_PRODUCT
00042 
00043 namespace MBMesquite
00044 {
00045 
00046 /**\brief \f$ R *= s \f$ */
00047 template < unsigned D >
00048 inline void hess_scale_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha );
00049 
00050 inline void hess_scale( MsqMatrix< 3, 3 > R[6], double alpha )
00051 {
00052     hess_scale_t< 3 >( R, alpha );
00053 }
00054 
00055 inline void hess_scale( MsqMatrix< 2, 2 > R[3], double alpha )
00056 {
00057     hess_scale_t< 2 >( R, alpha );
00058 }
00059 
00060 /**\brief \f$ R = \alpha I_9 \f$
00061  *
00062  *\param R The 6 blocks of the upper triangular portion of a 9x9
00063  *         symmetric matrix.
00064  */
00065 inline void set_scaled_I( MsqMatrix< 3, 3 > R[6], double alpha );
00066 
00067 /**\brief \f$ R += \alpha I_9 \f$
00068  *
00069  *\param R The 6 blocks of the upper triangular portion of a 9x9
00070  *         symmetric matrix.
00071  */
00072 inline void pluseq_scaled_I( MsqMatrix< 3, 3 > R[6], double alpha );
00073 
00074 /**\brief \f$ R += \alpha I_9 \f$
00075  *
00076  *\param R The 6 blocks of the upper triangular portion of a 9x9
00077  *         symmetric matrix.
00078  */
00079 inline void pluseq_scaled_I( MsqMatrix< 3, 3 >& R, double alpha );
00080 
00081 /**\brief \f$ R = \alpha I_4 \f$
00082  *
00083  *\param R The 3 blocks of the upper triangular portion of a 4x4
00084  *         symmetric matrix.
00085  */
00086 inline void set_scaled_I( MsqMatrix< 2, 2 > R[3], double alpha );
00087 
00088 /**\brief \f$ R += \alpha I_4 \f$
00089  *
00090  *\param R The 3 blocks of the upper triangular portion of a 4x4
00091  *         symmetric matrix.
00092  */
00093 inline void pluseq_scaled_I( MsqMatrix< 2, 2 > R[3], double alpha );
00094 
00095 /**\brief \f$ R += \alpha I_4 \f$
00096  *
00097  *\param R The 3 blocks of the upper triangular portion of a 4x4
00098  *         symmetric matrix.
00099  */
00100 inline void pluseq_scaled_I( MsqMatrix< 2, 2 >& R, double alpha );
00101 
00102 /**\brief \f$ R += \alpha \frac{\partial}{\partial T}det(T) \f$
00103  *
00104  *\param R The 6 blocks of the upper triangular portion of a 9x9
00105  *         symmetric matrix.
00106  */
00107 inline void pluseq_scaled_2nd_deriv_of_det( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T );
00108 
00109 /**\brief \f$ R += \alpha \frac{\partial}{\partial T}det(T) \f$
00110  *
00111  *\param R The 3 blocks of the upper triangular portion of a 4x4
00112  *         symmetric matrix.
00113  */
00114 inline void pluseq_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha );
00115 inline void pluseq_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& )
00116 {
00117     pluseq_scaled_2nd_deriv_of_det( R, alpha );
00118 }
00119 
00120 /**\brief \f$ R = \alpha \frac{\partial}{\partial T}det(T) \f$
00121  *
00122  *\param R The 6 blocks of the upper triangular portion of a 9x9
00123  *         symmetric matrix.
00124  */
00125 inline void set_scaled_2nd_deriv_of_det( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T );
00126 
00127 /**\brief \f$ R = \alpha \frac{\partial}{\partial T}det(T) \f$
00128  *
00129  *\param R The 3 blocks of the upper triangular portion of a 4x4
00130  *         symmetric matrix.
00131  */
00132 inline void set_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha );
00133 inline void set_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& )
00134 {
00135     set_scaled_2nd_deriv_of_det( R, alpha );
00136 }
00137 
00138 /**\brief \f$ R += \alpha \frac{\partial^2}{\partial T^2}tr(adj T) \f$
00139  *
00140  * Note: 2nd deriv is a constant, so T is not passed as an argument.
00141  */
00142 inline void pluseq_scaled_2nd_deriv_tr_adj( MsqMatrix< 2, 2 > R[3], double alpha );
00143 
00144 /**\brief \f$ R += \alpha \frac{\partial^2}{\partial T^2}tr(adj T) \f$
00145  *
00146  * Note: 2nd deriv is a constant, so T is not passed as an argument.
00147  */
00148 inline void pluseq_scaled_2nd_deriv_tr_adj( MsqMatrix< 3, 3 > R[6], double alpha );
00149 
00150 /**\brief \f$ R += \alpha \frac{\partial^2}{\partial T^2}|adj T|^2 \f$
00151  */
00152 inline void set_scaled_2nd_deriv_norm_sqr_adj( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& T );
00153 
00154 /**\brief \f$ R += \alpha \frac{\partial^2}{\partial T^2}|adj T|^2 \f$
00155  */
00156 inline void set_scaled_2nd_deriv_norm_sqr_adj( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T );
00157 
00158 /**\brief \f$ R += \alpha \left( M \otimes M \right) \f$
00159  *
00160  *\param R The 6 blocks of the upper triangular portion of a 9x9
00161  *         symmetric matrix.
00162  */
00163 template < unsigned D >
00164 inline void pluseq_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha,
00165                                            const MsqMatrix< D, D >& M );
00166 inline void pluseq_scaled_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& M )
00167 {
00168     pluseq_scaled_outer_product_t< 3 >( R, alpha, M );
00169 }
00170 
00171 inline void pluseq_scaled_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& M )
00172 {
00173     pluseq_scaled_outer_product_t< 2 >( R, alpha, M );
00174 }
00175 
00176 /**\brief \f$ R = \alpha \left( M \otimes M \right) \f$
00177  *
00178  *\param R The 6 blocks of the upper triangular portion of a 9x9
00179  *         symmetric matrix.
00180  */
00181 template < unsigned D >
00182 inline void set_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha,
00183                                         const MsqMatrix< D, D >& M );
00184 
00185 inline void set_scaled_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& M )
00186 {
00187     set_scaled_outer_product_t< 3 >( R, alpha, M );
00188 }
00189 
00190 inline void set_scaled_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& M )
00191 {
00192     set_scaled_outer_product_t< 2 >( R, alpha, M );
00193 }
00194 
00195 /**\brief \f$ R += \alpha \left( A \otimes B + B \otimes A \right) \f$
00196  *
00197  *\param R The 6 blocks of the upper triangular portion of a 9x9
00198  *         symmetric matrix.
00199  */
00200 inline void pluseq_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A,
00201                                              const MsqMatrix< 3, 3 >& B );
00202 
00203 /**\brief \f$ R += \alpha \left( A \otimes B + B \otimes A \right) \f$
00204  *
00205  *\param R The 3 blocks of the upper triangular portion of a 4x4
00206  *         symmetric matrix.
00207  */
00208 inline void pluseq_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& A,
00209                                              const MsqMatrix< 2, 2 >& B );
00210 
00211 /**\brief \f$ R = \alpha \left( A \otimes B + B \otimes A \right) \f$
00212  *
00213  *\param R The 6 blocks of the upper triangular portion of a 9x9
00214  *         symmetric matrix.
00215  */
00216 inline void set_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A,
00217                                           const MsqMatrix< 3, 3 >& B );
00218 
00219 /**\brief \f$ R = \alpha \left( A \otimes B + B \otimes A \right) \f$
00220  *
00221  *\param R The 3 blocks of the upper triangular portion of a 4x4
00222  *         symmetric matrix.
00223  */
00224 inline void set_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& A,
00225                                           const MsqMatrix< 2, 2 >& B );
00226 
00227 /**\brief \f$ R += \alpha (I \otimes I) \f$
00228  *
00229  *\param R The 6 blocks of the upper triangular portion of a 9x9
00230  *         symmetric matrix.
00231  */
00232 inline void pluseq_scaled_outer_product_I_I( MsqMatrix< 3, 3 > R[6], double alpha );
00233 
00234 /**\brief \f$ R += \alpha (I \otimes I) \f$
00235  *
00236  *\param R The 3 blocks of the upper triangular portion of a 4x4
00237  *         symmetric matrix.
00238  */
00239 inline void pluseq_scaled_outer_product_I_I( MsqMatrix< 2, 2 > R[3], double alpha );
00240 
00241 /**\brief \f$ R += I \otimes A \f$
00242  *
00243  *\param R The 6 blocks of the upper triangular portion of a 9x9
00244  *         symmetric matrix.
00245  */
00246 inline void pluseq_I_outer_product( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A );
00247 
00248 /**\brief \f$ R += A \otimes I \f$
00249  *
00250  *\param R The 6 blocks of the upper triangular portion of a 9x9
00251  *         symmetric matrix.
00252  */
00253 inline void pluseq_outer_product_I( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A );
00254 
00255 /**\brief \f$ R += \alpha \left( I \otimes A + A \otimes I \right) \f$
00256  *
00257  *\param R The 6 blocks of the upper triangular portion of a 9x9
00258  *         symmetric matrix.
00259  */
00260 template < unsigned D >
00261 inline void pluseq_scaled_sum_outer_product_I_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha,
00262                                                  const MsqMatrix< D, D >& A );
00263 inline void pluseq_scaled_sum_outer_product_I( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A_in )
00264 {
00265     pluseq_scaled_sum_outer_product_I_t< 3 >( R, alpha, A_in );
00266 }
00267 
00268 inline void pluseq_scaled_sum_outer_product_I( MsqMatrix< 2, 2 > R[2], double alpha, const MsqMatrix< 2, 2 >& A_in )
00269 {
00270     pluseq_scaled_sum_outer_product_I_t< 2 >( R, alpha, A_in );
00271 }
00272 
00273 /**\brief \f$ \frac{\partial^2 f}{\partial (AZ)^2} \Rightarrow \frac{\partial^2 f}{\partial A^2} \f$
00274  *
00275  * Given the second derivatives of a function with respect to
00276  * a matrix product, calculate the derivatives of the function
00277  * with respect to the first matrix in the product.
00278  */
00279 template < unsigned D >
00280 inline void second_deriv_wrt_product_factor_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], const MsqMatrix< D, D >& Z );
00281 
00282 inline void second_deriv_wrt_product_factor( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& Z )
00283 {
00284     second_deriv_wrt_product_factor_t< 3 >( R, Z );
00285 }
00286 
00287 inline void second_deriv_wrt_product_factor( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& Z )
00288 {
00289     second_deriv_wrt_product_factor_t< 2 >( R, Z );
00290 }
00291 
00292 /**\brief \f$ R = \alpha * \frac{\partial^2 \psi(T)}{\partial T^2} \f$
00293  *
00294  * \f$ \psi(T) = \sqrt{ |T|^2 + 2 \tau } \f$
00295  */
00296 inline void set_scaled_2nd_deriv_wrt_psi( MsqMatrix< 2, 2 > R[3], const double alpha, const double psi,
00297                                           const MsqMatrix< 2, 2 >& T );
00298 
00299 /**\brief \f$  R = R + \alpha * Z \f$
00300  */
00301 template < unsigned D >
00302 inline void pluseq_scaled_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha,
00303                              const MsqMatrix< D, D > Z[D * ( D + 1 ) / 2] );
00304 
00305 void set_scaled_I( MsqMatrix< 3, 3 > R[6], double alpha )
00306 {
00307     R[0] = R[3] = R[5] = MsqMatrix< 3, 3 >( alpha );
00308     R[1] = R[2] = R[4] = MsqMatrix< 3, 3 >( 0.0 );
00309 }
00310 
00311 void pluseq_scaled_I( MsqMatrix< 3, 3 >& R, double alpha )
00312 {
00313     R( 0, 0 ) += alpha;
00314     R( 1, 1 ) += alpha;
00315     R( 2, 2 ) += alpha;
00316 }
00317 
00318 void pluseq_scaled_I( MsqMatrix< 2, 2 >& R, double alpha )
00319 {
00320     R( 0, 0 ) += alpha;
00321     R( 1, 1 ) += alpha;
00322 }
00323 
00324 void pluseq_scaled_I( MsqMatrix< 3, 3 > R[6], double alpha )
00325 {
00326     pluseq_scaled_I( R[0], alpha );
00327     pluseq_scaled_I( R[3], alpha );
00328     pluseq_scaled_I( R[5], alpha );
00329 }
00330 
00331 void set_scaled_I( MsqMatrix< 2, 2 > R[3], double alpha )
00332 {
00333     R[0] = R[2] = MsqMatrix< 2, 2 >( alpha );
00334     R[1]        = MsqMatrix< 2, 2 >( 0.0 );
00335 }
00336 
00337 void pluseq_scaled_I( MsqMatrix< 2, 2 > R[3], double alpha )
00338 {
00339     pluseq_scaled_I( R[0], alpha );
00340     pluseq_scaled_I( R[2], alpha );
00341 }
00342 
00343 void pluseq_scaled_2nd_deriv_of_det( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T )
00344 {
00345     MsqMatrix< 3, 3 > A( T );
00346     A *= alpha;
00347 
00348     R[1]( 0, 1 ) += A( 2, 2 );
00349     R[1]( 1, 0 ) -= A( 2, 2 );
00350     R[1]( 0, 2 ) -= A( 2, 1 );
00351     R[1]( 2, 0 ) += A( 2, 1 );
00352     R[1]( 1, 2 ) += A( 2, 0 );
00353     R[1]( 2, 1 ) -= A( 2, 0 );
00354 
00355     R[2]( 0, 1 ) -= A( 1, 2 );
00356     R[2]( 1, 0 ) += A( 1, 2 );
00357     R[2]( 0, 2 ) += A( 1, 1 );
00358     R[2]( 2, 0 ) -= A( 1, 1 );
00359     R[2]( 1, 2 ) -= A( 1, 0 );
00360     R[2]( 2, 1 ) += A( 1, 0 );
00361 
00362     R[4]( 0, 1 ) += A( 0, 2 );
00363     R[4]( 1, 0 ) -= A( 0, 2 );
00364     R[4]( 0, 2 ) -= A( 0, 1 );
00365     R[4]( 2, 0 ) += A( 0, 1 );
00366     R[4]( 1, 2 ) += A( 0, 0 );
00367     R[4]( 2, 1 ) -= A( 0, 0 );
00368 }
00369 
00370 void set_scaled_2nd_deriv_of_det( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T )
00371 {
00372     MsqMatrix< 3, 3 > A( T );
00373     A *= alpha;
00374 
00375     R[0] = R[3] = R[5] = MsqMatrix< 3, 3 >( 0.0 );
00376 
00377     R[1]( 0, 0 ) = R[1]( 1, 1 ) = R[1]( 2, 2 ) = 0;
00378     R[1]( 0, 1 )                               = A( 2, 2 );
00379     R[1]( 1, 0 )                               = -A( 2, 2 );
00380     R[1]( 0, 2 )                               = -A( 2, 1 );
00381     R[1]( 2, 0 )                               = A( 2, 1 );
00382     R[1]( 1, 2 )                               = A( 2, 0 );
00383     R[1]( 2, 1 )                               = -A( 2, 0 );
00384 
00385     R[2]( 0, 0 ) = R[2]( 1, 1 ) = R[2]( 2, 2 ) = 0;
00386     R[2]( 0, 1 )                               = -A( 1, 2 );
00387     R[2]( 1, 0 )                               = A( 1, 2 );
00388     R[2]( 0, 2 )                               = A( 1, 1 );
00389     R[2]( 2, 0 )                               = -A( 1, 1 );
00390     R[2]( 1, 2 )                               = -A( 1, 0 );
00391     R[2]( 2, 1 )                               = A( 1, 0 );
00392 
00393     R[4]( 0, 0 ) = R[4]( 1, 1 ) = R[4]( 2, 2 ) = 0;
00394     R[4]( 0, 1 )                               = A( 0, 2 );
00395     R[4]( 1, 0 )                               = -A( 0, 2 );
00396     R[4]( 0, 2 )                               = -A( 0, 1 );
00397     R[4]( 2, 0 )                               = A( 0, 1 );
00398     R[4]( 1, 2 )                               = A( 0, 0 );
00399     R[4]( 2, 1 )                               = -A( 0, 0 );
00400 }
00401 
00402 void pluseq_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha )
00403 {
00404     R[1]( 0, 1 ) += alpha;
00405     R[1]( 1, 0 ) -= alpha;
00406 }
00407 
00408 void set_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha )
00409 {
00410     R[0] = R[2]  = MsqMatrix< 2, 2 >( 0.0 );
00411     R[1]( 0, 0 ) = 0.0;
00412     R[1]( 0, 1 ) = alpha;
00413     R[1]( 1, 0 ) = -alpha;
00414     R[1]( 1, 1 ) = 0.0;
00415 }
00416 
00417 void pluseq_scaled_2nd_deriv_tr_adj( MsqMatrix< 2, 2 >*, double )
00418 {
00419     // 2nd derivative is zero
00420 }
00421 
00422 void pluseq_scaled_2nd_deriv_tr_adj( MsqMatrix< 3, 3 > R[6], double alpha )
00423 {
00424     R[1]( 0, 1 ) += alpha;
00425     R[1]( 1, 0 ) -= alpha;
00426     R[2]( 0, 2 ) += alpha;
00427     R[2]( 2, 0 ) -= alpha;
00428     R[4]( 1, 2 ) += alpha;
00429     R[4]( 2, 1 ) -= alpha;
00430 }
00431 
00432 void set_scaled_2nd_deriv_norm_sqr_adj( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& )
00433 {
00434     set_scaled_I( R, 2 * alpha );
00435 }
00436 
00437 void set_scaled_2nd_deriv_norm_sqr_adj( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T )
00438 {
00439     set_scaled_outer_product( R, 1, T );
00440     double tmp01, tmp02, tmp12;
00441     // R[1] = 2*R[1] - transpose(R[1]);
00442     tmp01        = R[1]( 0, 1 );
00443     tmp02        = R[1]( 0, 2 );
00444     tmp12        = R[1]( 1, 2 );
00445     R[1]( 0, 1 ) = 2 * R[1]( 0, 1 ) - R[1]( 1, 0 );
00446     R[1]( 0, 2 ) = 2 * R[1]( 0, 2 ) - R[1]( 2, 0 );
00447     R[1]( 1, 2 ) = 2 * R[1]( 1, 2 ) - R[1]( 2, 1 );
00448     R[1]( 1, 0 ) = 2 * R[1]( 1, 0 ) - tmp01;
00449     R[1]( 2, 0 ) = 2 * R[1]( 2, 0 ) - tmp02;
00450     R[1]( 2, 1 ) = 2 * R[1]( 2, 1 ) - tmp12;
00451     // R[2] = 2*R[2] - transpose(R[1]);
00452     tmp01        = R[2]( 0, 1 );
00453     tmp02        = R[2]( 0, 2 );
00454     tmp12        = R[2]( 1, 2 );
00455     R[2]( 0, 1 ) = 2 * R[2]( 0, 1 ) - R[2]( 1, 0 );
00456     R[2]( 0, 2 ) = 2 * R[2]( 0, 2 ) - R[2]( 2, 0 );
00457     R[2]( 1, 2 ) = 2 * R[2]( 1, 2 ) - R[2]( 2, 1 );
00458     R[2]( 1, 0 ) = 2 * R[2]( 1, 0 ) - tmp01;
00459     R[2]( 2, 0 ) = 2 * R[2]( 2, 0 ) - tmp02;
00460     R[2]( 2, 1 ) = 2 * R[2]( 2, 1 ) - tmp12;
00461     // R[4] = 2*R[4] - transpose(R[1]);
00462     tmp01        = R[4]( 0, 1 );
00463     tmp02        = R[4]( 0, 2 );
00464     tmp12        = R[4]( 1, 2 );
00465     R[4]( 0, 1 ) = 2 * R[4]( 0, 1 ) - R[4]( 1, 0 );
00466     R[4]( 0, 2 ) = 2 * R[4]( 0, 2 ) - R[4]( 2, 0 );
00467     R[4]( 1, 2 ) = 2 * R[4]( 1, 2 ) - R[4]( 2, 1 );
00468     R[4]( 1, 0 ) = 2 * R[4]( 1, 0 ) - tmp01;
00469     R[4]( 2, 0 ) = 2 * R[4]( 2, 0 ) - tmp02;
00470     R[4]( 2, 1 ) = 2 * R[4]( 2, 1 ) - tmp12;
00471 
00472     const MsqMatrix< 3, 3 > TpT = transpose( T ) * T;
00473     R[0] -= TpT;
00474     R[3] -= TpT;
00475     R[5] -= TpT;
00476 
00477     const MsqMatrix< 3, 3 > TTp = T * transpose( T );
00478     const double ns             = sqr_Frobenius( T );
00479     R[0]( 0, 0 ) += ns - TTp( 0, 0 );
00480     R[0]( 1, 1 ) += ns - TTp( 0, 0 );
00481     R[0]( 2, 2 ) += ns - TTp( 0, 0 );
00482     R[1]( 0, 0 ) += -TTp( 0, 1 );
00483     R[1]( 1, 1 ) += -TTp( 0, 1 );
00484     R[1]( 2, 2 ) += -TTp( 0, 1 );
00485     R[2]( 0, 0 ) += -TTp( 0, 2 );
00486     R[2]( 1, 1 ) += -TTp( 0, 2 );
00487     R[2]( 2, 2 ) += -TTp( 0, 2 );
00488     R[3]( 0, 0 ) += ns - TTp( 1, 1 );
00489     R[3]( 1, 1 ) += ns - TTp( 1, 1 );
00490     R[3]( 2, 2 ) += ns - TTp( 1, 1 );
00491     R[4]( 0, 0 ) += -TTp( 1, 2 );
00492     R[4]( 1, 1 ) += -TTp( 1, 2 );
00493     R[4]( 2, 2 ) += -TTp( 1, 2 );
00494     R[5]( 0, 0 ) += ns - TTp( 2, 2 );
00495     R[5]( 1, 1 ) += ns - TTp( 2, 2 );
00496     R[5]( 2, 2 ) += ns - TTp( 2, 2 );
00497 
00498     alpha *= 2;
00499     R[0] *= alpha;
00500     R[1] *= alpha;
00501     R[2] *= alpha;
00502     R[3] *= alpha;
00503     R[4] *= alpha;
00504     R[5] *= alpha;
00505 }
00506 
00507 #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
00508 template < unsigned D >
00509 void pluseq_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha, const MsqMatrix< D, D >& M )
00510 {
00511     MsqMatrix< D, D > aM( M );
00512     aM *= alpha;
00513     unsigned h = 0;
00514     for( unsigned i = 0; i < D; ++i )
00515         for( unsigned j = i; j < D; ++j )
00516             R[h++] += transpose( aM.row( i ) ) * M.row( j );
00517 }
00518 #else
00519 template < unsigned D >
00520 void pluseq_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha, const MsqMatrix< D, D >& M )
00521 {
00522     MsqMatrix< D, D > aM( transpose( M ) );
00523     aM *= alpha;
00524     unsigned h = 0;
00525     for( unsigned i = 0; i < D; ++i )
00526         for( unsigned j = i; j < D; ++j )
00527             R[h++] += M.column( i ) * aM.row( 0 );
00528 }
00529 #endif
00530 
00531 #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
00532 template < unsigned D >
00533 void set_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha, const MsqMatrix< D, D >& M )
00534 {
00535     MsqMatrix< D, D > aM( M );
00536     aM *= alpha;
00537     unsigned h = 0;
00538     for( unsigned i = 0; i < D; ++i )
00539         for( unsigned j = i; j < D; ++j )
00540             R[h++] = transpose( aM.row( i ) ) * M.row( j );
00541 }
00542 #else
00543 template < unsigned D >
00544 void set_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha, const MsqMatrix< D, D >& M )
00545 {
00546     MsqMatrix< D, D > aM( transpose( M ) );
00547     aM *= alpha;
00548     unsigned h = 0;
00549     for( unsigned i = 0; i < D; ++i )
00550         for( unsigned j = i; j < D; ++j )
00551             R[h++] = M.column( i ) * aM.row( j );
00552 }
00553 #endif
00554 
00555 #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
00556 void pluseq_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A_in,
00557                                       const MsqMatrix< 3, 3 >& B )
00558 {
00559     // apply scalar first
00560     MsqMatrix< 3, 3 > A( A_in ), tmp;
00561     A *= alpha;
00562 
00563     // block 0,0
00564     tmp = transpose( A.row( 0 ) ) * B.row( 0 );
00565     R[0] += tmp;
00566     R[0] += transpose( tmp );
00567 
00568     // block 1,1
00569     tmp = transpose( A.row( 1 ) ) * B.row( 1 );
00570     R[3] += tmp;
00571     R[3] += transpose( tmp );
00572 
00573     // block 2,2
00574     tmp = transpose( A.row( 2 ) ) * B.row( 2 );
00575     R[5] += tmp;
00576     R[5] += transpose( tmp );
00577 
00578     // block 0,1
00579     R[1] += transpose( A.row( 0 ) ) * B.row( 1 ) + transpose( B.row( 0 ) ) * A.row( 1 );
00580 
00581     // block 0,2
00582     R[2] += transpose( A.row( 0 ) ) * B.row( 2 ) + transpose( B.row( 0 ) ) * A.row( 2 );
00583 
00584     // block 1,2
00585     R[4] += transpose( A.row( 1 ) ) * B.row( 2 ) + transpose( B.row( 1 ) ) * A.row( 2 );
00586 }
00587 void set_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A_in,
00588                                    const MsqMatrix< 3, 3 >& B )
00589 {
00590     // apply scalar first
00591     MsqMatrix< 3, 3 > A( A_in );
00592     A *= alpha;
00593 
00594     // block 0,0
00595     R[0] = transpose( A.row( 0 ) ) * B.row( 0 );
00596     R[0] += transpose( R[0] );
00597 
00598     // block 1,1
00599     R[3] = transpose( A.row( 1 ) ) * B.row( 1 );
00600     R[3] += transpose( R[3] );
00601 
00602     // block 2,2
00603     R[5] = transpose( A.row( 2 ) ) * B.row( 2 );
00604     R[5] += transpose( R[5] );
00605 
00606     // block 0,1
00607     R[1] = transpose( A.row( 0 ) ) * B.row( 1 );
00608     R[1] += transpose( B.row( 0 ) ) * A.row( 1 );
00609 
00610     // block 0,2
00611     R[2] = transpose( A.row( 0 ) ) * B.row( 2 );
00612     R[2] += transpose( B.row( 0 ) ) * A.row( 2 );
00613 
00614     // block 1,2
00615     R[4] = transpose( A.row( 1 ) ) * B.row( 2 );
00616     R[4] += transpose( B.row( 1 ) ) * A.row( 2 );
00617 }
00618 #else
00619 void pluseq_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A_in,
00620                                       const MsqMatrix< 3, 3 >& B )
00621 {
00622     // apply scalar first
00623     MsqMatrix< 3, 3 > A( A_in ), tmp;
00624     A *= alpha;
00625 
00626     // block 0,0
00627     tmp = A.column( 0 ) * transpose( B.column( 0 ) );
00628     R[0] += tmp;
00629     R[0] += transpose( tmp );
00630 
00631     // block 1,1
00632     tmp = A.column( 1 ) * transpose( B.column( 1 ) );
00633     R[3] += tmp;
00634     R[3] += transpose( tmp );
00635 
00636     // block 2,2
00637     tmp = A.column( 2 ) * transpose( B.column( 2 ) );
00638     R[5] += tmp;
00639     R[5] += transpose( tmp );
00640 
00641     // block 0,1
00642     R[1] += A.column( 0 ) * transpose( B.column( 1 ) ) + B.column( 0 ) * transpose( A.column( 1 ) );
00643 
00644     // block 0,2
00645     R[2] += A.column( 0 ) * transpose( B.column( 2 ) ) + B.column( 0 ) * transpose( A.column( 2 ) );
00646 
00647     // block 1,2
00648     R[4] += A.column( 1 ) * transpose( B.column( 2 ) ) + B.column( 1 ) * transpose( A.column( 2 ) );
00649 }
00650 void set_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A_in,
00651                                    const MsqMatrix< 3, 3 >& B )
00652 {
00653     // apply scalar first
00654     MsqMatrix< 3, 3 > A( A_in );
00655     A *= alpha;
00656 
00657     // block 0,0
00658     R[0] = A.column( 0 ) * transpose( B.column( 0 ) );
00659     R[0] += transpose( R[0] );
00660 
00661     // block 1,1
00662     R[3] = A.column( 1 ) * transpose( B.column( 1 ) );
00663     R[3] += transpose( R[3] );
00664 
00665     // block 2,2
00666     R[5] = A.column( 2 ) * transpose( B.column( 2 ) );
00667     R[5] += transpose( R[5] );
00668 
00669     // block 0,1
00670     R[1] = A.column( 0 ) * transpose( B.column( 1 ) );
00671     R[1] += B.column( 0 ) * transpose( A.column( 1 ) );
00672 
00673     // block 0,2
00674     R[2] = A.column( 0 ) * transpose( B.column( 2 ) );
00675     R[2] += B.column( 0 ) * transpose( A.column( 2 ) );
00676 
00677     // block 1,2
00678     R[4] = A.column( 1 ) * transpose( B.column( 2 ) );
00679     R[4] += B.column( 1 ) * transpose( A.column( 2 ) );
00680 }
00681 #endif
00682 
00683 #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
00684 void pluseq_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& A_in,
00685                                       const MsqMatrix< 2, 2 >& B )
00686 {
00687     // apply scalar first
00688     MsqMatrix< 2, 2 > A( A_in ), tmp;
00689     A *= alpha;
00690 
00691     // block 0,0
00692     tmp = transpose( A.row( 0 ) ) * B.row( 0 );
00693     R[0] += tmp;
00694     R[0] += transpose( tmp );
00695 
00696     // block 1,1
00697     tmp = transpose( A.row( 1 ) ) * B.row( 1 );
00698     R[2] += tmp;
00699     R[2] += transpose( tmp );
00700 
00701     // block 0,1
00702     R[1] += transpose( A.row( 0 ) ) * B.row( 1 ) + transpose( B.row( 0 ) ) * A.row( 1 );
00703 }
00704 void set_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& A_in,
00705                                    const MsqMatrix< 2, 2 >& B )
00706 {
00707     // apply scalar first
00708     MsqMatrix< 2, 2 > A( A_in );
00709     A *= alpha;
00710 
00711     // block 0,0
00712     R[0] = transpose( A.row( 0 ) ) * B.row( 0 );
00713     R[0] += transpose( R[0] );
00714 
00715     // block 1,1
00716     R[2] = transpose( A.row( 1 ) ) * B.row( 1 );
00717     R[2] += transpose( R[2] );
00718 
00719     // block 0,1
00720     R[1] = transpose( A.row( 0 ) ) * B.row( 1 );
00721     R[1] += transpose( B.row( 0 ) ) * A.row( 1 );
00722 }
00723 #else
00724 void pluseq_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& A_in,
00725                                       const MsqMatrix< 2, 2 >& B )
00726 {
00727     // apply scalar first
00728     MsqMatrix< 2, 2 > A( A_in ), tmp;
00729     A *= alpha;
00730 
00731     // block 0,0
00732     tmp = A.column( 0 ) * transpose( B.column( 0 ) );
00733     R[0] += tmp;
00734     R[0] += transpose( tmp );
00735 
00736     // block 1,1
00737     tmp = A.column( 1 ) * transpose( B.column( 1 ) );
00738     R[2] += tmp;
00739     R[2] += transpose( tmp );
00740 
00741     // block 0,1
00742     R[1] += A.column( 0 ) * transpose( B.column( 1 ) ) + B.column( 0 ) * transpose( A.column( 1 ) );
00743 }
00744 void set_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& A_in,
00745                                    const MsqMatrix< 2, 2 >& B )
00746 {
00747     // apply scalar first
00748     MsqMatrix< 2, 2 > A( A_in );
00749     A *= alpha;
00750 
00751     // block 0,0
00752     R[0] = A.column( 0 ) * transpose( B.column( 0 ) );
00753     R[0] += transpose( R[0] );
00754 
00755     // block 1,1
00756     R[2] = A.column( 1 ) * transpose( B.column( 1 ) );
00757     R[2] += transpose( R[2] );
00758 
00759     // block 0,1
00760     R[1] = A.column( 0 ) * transpose( B.column( 1 ) );
00761     R[1] += B.column( 0 ) * transpose( A.column( 1 ) );
00762 }
00763 #endif
00764 
00765 void pluseq_scaled_outer_product_I_I( MsqMatrix< 3, 3 > R[6], double alpha )
00766 {
00767     R[0]( 0, 0 ) += alpha;
00768     R[1]( 0, 1 ) += alpha;
00769     R[2]( 0, 2 ) += alpha;
00770     R[3]( 1, 1 ) += alpha;
00771     R[4]( 1, 2 ) += alpha;
00772     R[5]( 2, 2 ) += alpha;
00773 }
00774 
00775 void pluseq_scaled_outer_product_I_I( MsqMatrix< 2, 2 > R[3], double alpha )
00776 {
00777     R[0]( 0, 0 ) += alpha;
00778     R[1]( 0, 1 ) += alpha;
00779     R[2]( 1, 1 ) += alpha;
00780 }
00781 
00782 #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
00783 inline void pluseq_I_outer_product( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A )
00784 {
00785     R[0].add_row( 0, A.row( 0 ) );
00786     R[1].add_row( 0, A.row( 1 ) );
00787     R[2].add_row( 0, A.row( 2 ) );
00788     R[3].add_row( 1, A.row( 1 ) );
00789     R[4].add_row( 1, A.row( 2 ) );
00790     R[5].add_row( 2, A.row( 2 ) );
00791 }
00792 inline void pluseq_I_outer_product( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& A )
00793 {
00794     R[0].add_row( 0, A.row( 0 ) );
00795     R[1].add_row( 0, A.row( 1 ) );
00796     R[2].add_row( 1, A.row( 1 ) );
00797 }
00798 #else
00799 inline void pluseq_I_outer_product( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A )
00800 {
00801     R[0].add_row( 0, transpose( A.column( 0 ) ) );
00802     R[1].add_row( 0, transpose( A.column( 1 ) ) );
00803     R[2].add_row( 0, transpose( A.column( 2 ) ) );
00804     R[3].add_row( 1, transpose( A.column( 1 ) ) );
00805     R[4].add_row( 1, transpose( A.column( 2 ) ) );
00806     R[5].add_row( 2, transpose( A.column( 2 ) ) );
00807 }
00808 inline void pluseq_I_outer_product( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& A )
00809 {
00810     R[0].add_row( 0, transpose( A.column( 0 ) ) );
00811     R[1].add_row( 0, transpose( A.column( 1 ) ) );
00812     R[2].add_row( 1, transpose( A.column( 1 ) ) );
00813 }
00814 #endif
00815 
00816 #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
00817 inline void pluseq_outer_product_I( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A )
00818 {
00819     R[0].add_column( 0, transpose( A.row( 0 ) ) );
00820     R[1].add_column( 1, transpose( A.row( 0 ) ) );
00821     R[2].add_column( 2, transpose( A.row( 0 ) ) );
00822     R[3].add_column( 1, transpose( A.row( 1 ) ) );
00823     R[4].add_column( 2, transpose( A.row( 1 ) ) );
00824     R[5].add_column( 2, transpose( A.row( 2 ) ) );
00825 }
00826 inline void pluseq_outer_product_I( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& A )
00827 {
00828     R[0].add_column( 0, transpose( A.row( 0 ) ) );
00829     R[1].add_column( 1, transpose( A.row( 0 ) ) );
00830     R[2].add_column( 1, transpose( A.row( 1 ) ) );
00831 }
00832 #else
00833 inline void pluseq_outer_product_I( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A )
00834 {
00835     R[0].add_column( 0, A.column( 0 ) );
00836     R[1].add_column( 1, A.column( 0 ) );
00837     R[2].add_column( 2, A.column( 0 ) );
00838     R[3].add_column( 1, A.column( 1 ) );
00839     R[4].add_column( 2, A.column( 1 ) );
00840     R[5].add_column( 2, A.column( 2 ) );
00841 }
00842 inline void pluseq_outer_product_I( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& A )
00843 {
00844     R[0].add_column( 0, A.column( 0 ) );
00845     R[1].add_column( 1, A.column( 0 ) );
00846     R[2].add_column( 1, A.column( 1 ) );
00847 }
00848 #endif
00849 
00850 template < unsigned D >
00851 void pluseq_scaled_sum_outer_product_I_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha,
00852                                           const MsqMatrix< D, D >& A_in )
00853 {
00854     // apply scalar first
00855     MsqMatrix< D, D > A( A_in );
00856     A *= alpha;
00857     pluseq_I_outer_product( R, A );
00858     pluseq_outer_product_I( R, A );
00859 }
00860 
00861 template < unsigned D >
00862 void second_deriv_wrt_product_factor_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], const MsqMatrix< D, D >& Z )
00863 {
00864     const MsqMatrix< D, D > Zt = transpose( Z );
00865     for( unsigned i = 0; i < D * ( D + 1 ) / 2; ++i )
00866         R[i] = ( Z * R[i] ) * Zt;
00867 }
00868 
00869 void set_scaled_2nd_deriv_wrt_psi( MsqMatrix< 2, 2 > R[3], const double alpha, const double psi,
00870                                    const MsqMatrix< 2, 2 >& T )
00871 {
00872     const double t = trace( T );
00873     const double f = alpha / ( psi * psi * psi );
00874     const double s = T( 0, 1 ) - T( 1, 0 );
00875     R[0]( 0, 0 ) = R[1]( 0, 1 ) = R[2]( 1, 1 ) = f * s * s;
00876     R[0]( 0, 1 ) = R[0]( 1, 0 ) = R[1]( 1, 1 ) = -f * s * t;
00877     R[1]( 0, 0 ) = R[2]( 0, 1 ) = R[2]( 1, 0 ) = f * s * t;
00878     R[0]( 1, 1 ) = R[2]( 0, 0 ) = -( R[1]( 1, 0 ) = -f * t * t );
00879 }
00880 
00881 template < unsigned D >
00882 inline void pluseq_scaled( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha,
00883                            const MsqMatrix< D, D > Z[D * ( D + 1 ) / 2] )
00884 {
00885     for( unsigned i = 0; i < D * ( D + 1 ) / 2; ++i )
00886         R[i] += alpha * Z[i];
00887 }
00888 
00889 /**\brief \f$ R *= s */
00890 template < unsigned D >
00891 inline void hess_scale_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha )
00892 {
00893     for( unsigned i = 0; i < D * ( D + 1 ) / 2; ++i )
00894         R[i] *= alpha;
00895 }
00896 
00897 }  // namespace MBMesquite
00898 
00899 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines