MOAB: Mesh Oriented datABase  (version 5.4.0)
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],
00165                                            double alpha,
00166                                            const MsqMatrix< D, D >& M );
00167 inline void pluseq_scaled_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& M )
00168 {
00169     pluseq_scaled_outer_product_t< 3 >( R, alpha, M );
00170 }
00171 
00172 inline void pluseq_scaled_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& M )
00173 {
00174     pluseq_scaled_outer_product_t< 2 >( R, alpha, M );
00175 }
00176 
00177 /**\brief \f$ R = \alpha \left( M \otimes M \right) \f$
00178  *
00179  *\param R The 6 blocks of the upper triangular portion of a 9x9
00180  *         symmetric matrix.
00181  */
00182 template < unsigned D >
00183 inline void set_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2],
00184                                         double alpha,
00185                                         const MsqMatrix< D, D >& M );
00186 
00187 inline void set_scaled_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& M )
00188 {
00189     set_scaled_outer_product_t< 3 >( R, alpha, M );
00190 }
00191 
00192 inline void set_scaled_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& M )
00193 {
00194     set_scaled_outer_product_t< 2 >( R, alpha, M );
00195 }
00196 
00197 /**\brief \f$ R += \alpha \left( A \otimes B + B \otimes A \right) \f$
00198  *
00199  *\param R The 6 blocks of the upper triangular portion of a 9x9
00200  *         symmetric matrix.
00201  */
00202 inline void pluseq_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6],
00203                                              double alpha,
00204                                              const MsqMatrix< 3, 3 >& A,
00205                                              const MsqMatrix< 3, 3 >& B );
00206 
00207 /**\brief \f$ R += \alpha \left( A \otimes B + B \otimes A \right) \f$
00208  *
00209  *\param R The 3 blocks of the upper triangular portion of a 4x4
00210  *         symmetric matrix.
00211  */
00212 inline void pluseq_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3],
00213                                              double alpha,
00214                                              const MsqMatrix< 2, 2 >& A,
00215                                              const MsqMatrix< 2, 2 >& B );
00216 
00217 /**\brief \f$ R = \alpha \left( A \otimes B + B \otimes A \right) \f$
00218  *
00219  *\param R The 6 blocks of the upper triangular portion of a 9x9
00220  *         symmetric matrix.
00221  */
00222 inline void set_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6],
00223                                           double alpha,
00224                                           const MsqMatrix< 3, 3 >& A,
00225                                           const MsqMatrix< 3, 3 >& B );
00226 
00227 /**\brief \f$ R = \alpha \left( A \otimes B + B \otimes A \right) \f$
00228  *
00229  *\param R The 3 blocks of the upper triangular portion of a 4x4
00230  *         symmetric matrix.
00231  */
00232 inline void set_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3],
00233                                           double alpha,
00234                                           const MsqMatrix< 2, 2 >& A,
00235                                           const MsqMatrix< 2, 2 >& B );
00236 
00237 /**\brief \f$ R += \alpha (I \otimes I) \f$
00238  *
00239  *\param R The 6 blocks of the upper triangular portion of a 9x9
00240  *         symmetric matrix.
00241  */
00242 inline void pluseq_scaled_outer_product_I_I( MsqMatrix< 3, 3 > R[6], double alpha );
00243 
00244 /**\brief \f$ R += \alpha (I \otimes I) \f$
00245  *
00246  *\param R The 3 blocks of the upper triangular portion of a 4x4
00247  *         symmetric matrix.
00248  */
00249 inline void pluseq_scaled_outer_product_I_I( MsqMatrix< 2, 2 > R[3], double alpha );
00250 
00251 /**\brief \f$ R += I \otimes A \f$
00252  *
00253  *\param R The 6 blocks of the upper triangular portion of a 9x9
00254  *         symmetric matrix.
00255  */
00256 inline void pluseq_I_outer_product( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A );
00257 
00258 /**\brief \f$ R += A \otimes I \f$
00259  *
00260  *\param R The 6 blocks of the upper triangular portion of a 9x9
00261  *         symmetric matrix.
00262  */
00263 inline void pluseq_outer_product_I( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A );
00264 
00265 /**\brief \f$ R += \alpha \left( I \otimes A + A \otimes I \right) \f$
00266  *
00267  *\param R The 6 blocks of the upper triangular portion of a 9x9
00268  *         symmetric matrix.
00269  */
00270 template < unsigned D >
00271 inline void pluseq_scaled_sum_outer_product_I_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2],
00272                                                  double alpha,
00273                                                  const MsqMatrix< D, D >& A );
00274 inline void pluseq_scaled_sum_outer_product_I( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A_in )
00275 {
00276     pluseq_scaled_sum_outer_product_I_t< 3 >( R, alpha, A_in );
00277 }
00278 
00279 inline void pluseq_scaled_sum_outer_product_I( MsqMatrix< 2, 2 > R[2], double alpha, const MsqMatrix< 2, 2 >& A_in )
00280 {
00281     pluseq_scaled_sum_outer_product_I_t< 2 >( R, alpha, A_in );
00282 }
00283 
00284 /**\brief \f$ \frac{\partial^2 f}{\partial (AZ)^2} \Rightarrow \frac{\partial^2 f}{\partial A^2} \f$
00285  *
00286  * Given the second derivatives of a function with respect to
00287  * a matrix product, calculate the derivatives of the function
00288  * with respect to the first matrix in the product.
00289  */
00290 template < unsigned D >
00291 inline void second_deriv_wrt_product_factor_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], const MsqMatrix< D, D >& Z );
00292 
00293 inline void second_deriv_wrt_product_factor( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& Z )
00294 {
00295     second_deriv_wrt_product_factor_t< 3 >( R, Z );
00296 }
00297 
00298 inline void second_deriv_wrt_product_factor( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& Z )
00299 {
00300     second_deriv_wrt_product_factor_t< 2 >( R, Z );
00301 }
00302 
00303 /**\brief \f$ R = \alpha * \frac{\partial^2 \psi(T)}{\partial T^2} \f$
00304  *
00305  * \f$ \psi(T) = \sqrt{ |T|^2 + 2 \tau } \f$
00306  */
00307 inline void set_scaled_2nd_deriv_wrt_psi( MsqMatrix< 2, 2 > R[3],
00308                                           const double alpha,
00309                                           const double psi,
00310                                           const MsqMatrix< 2, 2 >& T );
00311 
00312 /**\brief \f$  R = R + \alpha * Z \f$
00313  */
00314 template < unsigned D >
00315 inline void pluseq_scaled_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2],
00316                              double alpha,
00317                              const MsqMatrix< D, D > Z[D * ( D + 1 ) / 2] );
00318 
00319 void set_scaled_I( MsqMatrix< 3, 3 > R[6], double alpha )
00320 {
00321     R[0] = R[3] = R[5] = MsqMatrix< 3, 3 >( alpha );
00322     R[1] = R[2] = R[4] = MsqMatrix< 3, 3 >( 0.0 );
00323 }
00324 
00325 void pluseq_scaled_I( MsqMatrix< 3, 3 >& R, double alpha )
00326 {
00327     R( 0, 0 ) += alpha;
00328     R( 1, 1 ) += alpha;
00329     R( 2, 2 ) += alpha;
00330 }
00331 
00332 void pluseq_scaled_I( MsqMatrix< 2, 2 >& R, double alpha )
00333 {
00334     R( 0, 0 ) += alpha;
00335     R( 1, 1 ) += alpha;
00336 }
00337 
00338 void pluseq_scaled_I( MsqMatrix< 3, 3 > R[6], double alpha )
00339 {
00340     pluseq_scaled_I( R[0], alpha );
00341     pluseq_scaled_I( R[3], alpha );
00342     pluseq_scaled_I( R[5], alpha );
00343 }
00344 
00345 void set_scaled_I( MsqMatrix< 2, 2 > R[3], double alpha )
00346 {
00347     R[0] = R[2] = MsqMatrix< 2, 2 >( alpha );
00348     R[1]        = MsqMatrix< 2, 2 >( 0.0 );
00349 }
00350 
00351 void pluseq_scaled_I( MsqMatrix< 2, 2 > R[3], double alpha )
00352 {
00353     pluseq_scaled_I( R[0], alpha );
00354     pluseq_scaled_I( R[2], alpha );
00355 }
00356 
00357 void pluseq_scaled_2nd_deriv_of_det( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T )
00358 {
00359     MsqMatrix< 3, 3 > A( T );
00360     A *= alpha;
00361 
00362     R[1]( 0, 1 ) += A( 2, 2 );
00363     R[1]( 1, 0 ) -= A( 2, 2 );
00364     R[1]( 0, 2 ) -= A( 2, 1 );
00365     R[1]( 2, 0 ) += A( 2, 1 );
00366     R[1]( 1, 2 ) += A( 2, 0 );
00367     R[1]( 2, 1 ) -= A( 2, 0 );
00368 
00369     R[2]( 0, 1 ) -= A( 1, 2 );
00370     R[2]( 1, 0 ) += A( 1, 2 );
00371     R[2]( 0, 2 ) += A( 1, 1 );
00372     R[2]( 2, 0 ) -= A( 1, 1 );
00373     R[2]( 1, 2 ) -= A( 1, 0 );
00374     R[2]( 2, 1 ) += A( 1, 0 );
00375 
00376     R[4]( 0, 1 ) += A( 0, 2 );
00377     R[4]( 1, 0 ) -= A( 0, 2 );
00378     R[4]( 0, 2 ) -= A( 0, 1 );
00379     R[4]( 2, 0 ) += A( 0, 1 );
00380     R[4]( 1, 2 ) += A( 0, 0 );
00381     R[4]( 2, 1 ) -= A( 0, 0 );
00382 }
00383 
00384 void set_scaled_2nd_deriv_of_det( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T )
00385 {
00386     MsqMatrix< 3, 3 > A( T );
00387     A *= alpha;
00388 
00389     R[0] = R[3] = R[5] = MsqMatrix< 3, 3 >( 0.0 );
00390 
00391     R[1]( 0, 0 ) = R[1]( 1, 1 ) = R[1]( 2, 2 ) = 0;
00392     R[1]( 0, 1 )                               = A( 2, 2 );
00393     R[1]( 1, 0 )                               = -A( 2, 2 );
00394     R[1]( 0, 2 )                               = -A( 2, 1 );
00395     R[1]( 2, 0 )                               = A( 2, 1 );
00396     R[1]( 1, 2 )                               = A( 2, 0 );
00397     R[1]( 2, 1 )                               = -A( 2, 0 );
00398 
00399     R[2]( 0, 0 ) = R[2]( 1, 1 ) = R[2]( 2, 2 ) = 0;
00400     R[2]( 0, 1 )                               = -A( 1, 2 );
00401     R[2]( 1, 0 )                               = A( 1, 2 );
00402     R[2]( 0, 2 )                               = A( 1, 1 );
00403     R[2]( 2, 0 )                               = -A( 1, 1 );
00404     R[2]( 1, 2 )                               = -A( 1, 0 );
00405     R[2]( 2, 1 )                               = A( 1, 0 );
00406 
00407     R[4]( 0, 0 ) = R[4]( 1, 1 ) = R[4]( 2, 2 ) = 0;
00408     R[4]( 0, 1 )                               = A( 0, 2 );
00409     R[4]( 1, 0 )                               = -A( 0, 2 );
00410     R[4]( 0, 2 )                               = -A( 0, 1 );
00411     R[4]( 2, 0 )                               = A( 0, 1 );
00412     R[4]( 1, 2 )                               = A( 0, 0 );
00413     R[4]( 2, 1 )                               = -A( 0, 0 );
00414 }
00415 
00416 void pluseq_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha )
00417 {
00418     R[1]( 0, 1 ) += alpha;
00419     R[1]( 1, 0 ) -= alpha;
00420 }
00421 
00422 void set_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha )
00423 {
00424     R[0] = R[2]  = MsqMatrix< 2, 2 >( 0.0 );
00425     R[1]( 0, 0 ) = 0.0;
00426     R[1]( 0, 1 ) = alpha;
00427     R[1]( 1, 0 ) = -alpha;
00428     R[1]( 1, 1 ) = 0.0;
00429 }
00430 
00431 void pluseq_scaled_2nd_deriv_tr_adj( MsqMatrix< 2, 2 >*, double )
00432 {
00433     // 2nd derivative is zero
00434 }
00435 
00436 void pluseq_scaled_2nd_deriv_tr_adj( MsqMatrix< 3, 3 > R[6], double alpha )
00437 {
00438     R[1]( 0, 1 ) += alpha;
00439     R[1]( 1, 0 ) -= alpha;
00440     R[2]( 0, 2 ) += alpha;
00441     R[2]( 2, 0 ) -= alpha;
00442     R[4]( 1, 2 ) += alpha;
00443     R[4]( 2, 1 ) -= alpha;
00444 }
00445 
00446 void set_scaled_2nd_deriv_norm_sqr_adj( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& )
00447 {
00448     set_scaled_I( R, 2 * alpha );
00449 }
00450 
00451 void set_scaled_2nd_deriv_norm_sqr_adj( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T )
00452 {
00453     set_scaled_outer_product( R, 1, T );
00454     double tmp01, tmp02, tmp12;
00455     // R[1] = 2*R[1] - transpose(R[1]);
00456     tmp01        = R[1]( 0, 1 );
00457     tmp02        = R[1]( 0, 2 );
00458     tmp12        = R[1]( 1, 2 );
00459     R[1]( 0, 1 ) = 2 * R[1]( 0, 1 ) - R[1]( 1, 0 );
00460     R[1]( 0, 2 ) = 2 * R[1]( 0, 2 ) - R[1]( 2, 0 );
00461     R[1]( 1, 2 ) = 2 * R[1]( 1, 2 ) - R[1]( 2, 1 );
00462     R[1]( 1, 0 ) = 2 * R[1]( 1, 0 ) - tmp01;
00463     R[1]( 2, 0 ) = 2 * R[1]( 2, 0 ) - tmp02;
00464     R[1]( 2, 1 ) = 2 * R[1]( 2, 1 ) - tmp12;
00465     // R[2] = 2*R[2] - transpose(R[1]);
00466     tmp01        = R[2]( 0, 1 );
00467     tmp02        = R[2]( 0, 2 );
00468     tmp12        = R[2]( 1, 2 );
00469     R[2]( 0, 1 ) = 2 * R[2]( 0, 1 ) - R[2]( 1, 0 );
00470     R[2]( 0, 2 ) = 2 * R[2]( 0, 2 ) - R[2]( 2, 0 );
00471     R[2]( 1, 2 ) = 2 * R[2]( 1, 2 ) - R[2]( 2, 1 );
00472     R[2]( 1, 0 ) = 2 * R[2]( 1, 0 ) - tmp01;
00473     R[2]( 2, 0 ) = 2 * R[2]( 2, 0 ) - tmp02;
00474     R[2]( 2, 1 ) = 2 * R[2]( 2, 1 ) - tmp12;
00475     // R[4] = 2*R[4] - transpose(R[1]);
00476     tmp01        = R[4]( 0, 1 );
00477     tmp02        = R[4]( 0, 2 );
00478     tmp12        = R[4]( 1, 2 );
00479     R[4]( 0, 1 ) = 2 * R[4]( 0, 1 ) - R[4]( 1, 0 );
00480     R[4]( 0, 2 ) = 2 * R[4]( 0, 2 ) - R[4]( 2, 0 );
00481     R[4]( 1, 2 ) = 2 * R[4]( 1, 2 ) - R[4]( 2, 1 );
00482     R[4]( 1, 0 ) = 2 * R[4]( 1, 0 ) - tmp01;
00483     R[4]( 2, 0 ) = 2 * R[4]( 2, 0 ) - tmp02;
00484     R[4]( 2, 1 ) = 2 * R[4]( 2, 1 ) - tmp12;
00485 
00486     const MsqMatrix< 3, 3 > TpT = transpose( T ) * T;
00487     R[0] -= TpT;
00488     R[3] -= TpT;
00489     R[5] -= TpT;
00490 
00491     const MsqMatrix< 3, 3 > TTp = T * transpose( T );
00492     const double ns             = sqr_Frobenius( T );
00493     R[0]( 0, 0 ) += ns - TTp( 0, 0 );
00494     R[0]( 1, 1 ) += ns - TTp( 0, 0 );
00495     R[0]( 2, 2 ) += ns - TTp( 0, 0 );
00496     R[1]( 0, 0 ) += -TTp( 0, 1 );
00497     R[1]( 1, 1 ) += -TTp( 0, 1 );
00498     R[1]( 2, 2 ) += -TTp( 0, 1 );
00499     R[2]( 0, 0 ) += -TTp( 0, 2 );
00500     R[2]( 1, 1 ) += -TTp( 0, 2 );
00501     R[2]( 2, 2 ) += -TTp( 0, 2 );
00502     R[3]( 0, 0 ) += ns - TTp( 1, 1 );
00503     R[3]( 1, 1 ) += ns - TTp( 1, 1 );
00504     R[3]( 2, 2 ) += ns - TTp( 1, 1 );
00505     R[4]( 0, 0 ) += -TTp( 1, 2 );
00506     R[4]( 1, 1 ) += -TTp( 1, 2 );
00507     R[4]( 2, 2 ) += -TTp( 1, 2 );
00508     R[5]( 0, 0 ) += ns - TTp( 2, 2 );
00509     R[5]( 1, 1 ) += ns - TTp( 2, 2 );
00510     R[5]( 2, 2 ) += ns - TTp( 2, 2 );
00511 
00512     alpha *= 2;
00513     R[0] *= alpha;
00514     R[1] *= alpha;
00515     R[2] *= alpha;
00516     R[3] *= alpha;
00517     R[4] *= alpha;
00518     R[5] *= alpha;
00519 }
00520 
00521 #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
00522 template < unsigned D >
00523 void pluseq_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha, const MsqMatrix< D, D >& M )
00524 {
00525     MsqMatrix< D, D > aM( M );
00526     aM *= alpha;
00527     unsigned h = 0;
00528     for( unsigned i = 0; i < D; ++i )
00529         for( unsigned j = i; j < D; ++j )
00530             R[h++] += transpose( aM.row( i ) ) * M.row( j );
00531 }
00532 #else
00533 template < unsigned D >
00534 void pluseq_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha, const MsqMatrix< D, D >& M )
00535 {
00536     MsqMatrix< D, D > aM( transpose( M ) );
00537     aM *= alpha;
00538     unsigned h = 0;
00539     for( unsigned i = 0; i < D; ++i )
00540         for( unsigned j = i; j < D; ++j )
00541             R[h++] += M.column( i ) * aM.row( 0 );
00542 }
00543 #endif
00544 
00545 #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
00546 template < unsigned D >
00547 void set_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha, const MsqMatrix< D, D >& M )
00548 {
00549     MsqMatrix< D, D > aM( M );
00550     aM *= alpha;
00551     unsigned h = 0;
00552     for( unsigned i = 0; i < D; ++i )
00553         for( unsigned j = i; j < D; ++j )
00554             R[h++] = transpose( aM.row( i ) ) * M.row( j );
00555 }
00556 #else
00557 template < unsigned D >
00558 void set_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha, const MsqMatrix< D, D >& M )
00559 {
00560     MsqMatrix< D, D > aM( transpose( M ) );
00561     aM *= alpha;
00562     unsigned h = 0;
00563     for( unsigned i = 0; i < D; ++i )
00564         for( unsigned j = i; j < D; ++j )
00565             R[h++] = M.column( i ) * aM.row( j );
00566 }
00567 #endif
00568 
00569 #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
00570 void pluseq_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6],
00571                                       double alpha,
00572                                       const MsqMatrix< 3, 3 >& A_in,
00573                                       const MsqMatrix< 3, 3 >& B )
00574 {
00575     // apply scalar first
00576     MsqMatrix< 3, 3 > A( A_in ), tmp;
00577     A *= alpha;
00578 
00579     // block 0,0
00580     tmp = transpose( A.row( 0 ) ) * B.row( 0 );
00581     R[0] += tmp;
00582     R[0] += transpose( tmp );
00583 
00584     // block 1,1
00585     tmp = transpose( A.row( 1 ) ) * B.row( 1 );
00586     R[3] += tmp;
00587     R[3] += transpose( tmp );
00588 
00589     // block 2,2
00590     tmp = transpose( A.row( 2 ) ) * B.row( 2 );
00591     R[5] += tmp;
00592     R[5] += transpose( tmp );
00593 
00594     // block 0,1
00595     R[1] += transpose( A.row( 0 ) ) * B.row( 1 ) + transpose( B.row( 0 ) ) * A.row( 1 );
00596 
00597     // block 0,2
00598     R[2] += transpose( A.row( 0 ) ) * B.row( 2 ) + transpose( B.row( 0 ) ) * A.row( 2 );
00599 
00600     // block 1,2
00601     R[4] += transpose( A.row( 1 ) ) * B.row( 2 ) + transpose( B.row( 1 ) ) * A.row( 2 );
00602 }
00603 void set_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6],
00604                                    double alpha,
00605                                    const MsqMatrix< 3, 3 >& A_in,
00606                                    const MsqMatrix< 3, 3 >& B )
00607 {
00608     // apply scalar first
00609     MsqMatrix< 3, 3 > A( A_in );
00610     A *= alpha;
00611 
00612     // block 0,0
00613     R[0] = transpose( A.row( 0 ) ) * B.row( 0 );
00614     R[0] += transpose( R[0] );
00615 
00616     // block 1,1
00617     R[3] = transpose( A.row( 1 ) ) * B.row( 1 );
00618     R[3] += transpose( R[3] );
00619 
00620     // block 2,2
00621     R[5] = transpose( A.row( 2 ) ) * B.row( 2 );
00622     R[5] += transpose( R[5] );
00623 
00624     // block 0,1
00625     R[1] = transpose( A.row( 0 ) ) * B.row( 1 );
00626     R[1] += transpose( B.row( 0 ) ) * A.row( 1 );
00627 
00628     // block 0,2
00629     R[2] = transpose( A.row( 0 ) ) * B.row( 2 );
00630     R[2] += transpose( B.row( 0 ) ) * A.row( 2 );
00631 
00632     // block 1,2
00633     R[4] = transpose( A.row( 1 ) ) * B.row( 2 );
00634     R[4] += transpose( B.row( 1 ) ) * A.row( 2 );
00635 }
00636 #else
00637 void pluseq_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6],
00638                                       double alpha,
00639                                       const MsqMatrix< 3, 3 >& A_in,
00640                                       const MsqMatrix< 3, 3 >& B )
00641 {
00642     // apply scalar first
00643     MsqMatrix< 3, 3 > A( A_in ), tmp;
00644     A *= alpha;
00645 
00646     // block 0,0
00647     tmp = A.column( 0 ) * transpose( B.column( 0 ) );
00648     R[0] += tmp;
00649     R[0] += transpose( tmp );
00650 
00651     // block 1,1
00652     tmp = A.column( 1 ) * transpose( B.column( 1 ) );
00653     R[3] += tmp;
00654     R[3] += transpose( tmp );
00655 
00656     // block 2,2
00657     tmp = A.column( 2 ) * transpose( B.column( 2 ) );
00658     R[5] += tmp;
00659     R[5] += transpose( tmp );
00660 
00661     // block 0,1
00662     R[1] += A.column( 0 ) * transpose( B.column( 1 ) ) + B.column( 0 ) * transpose( A.column( 1 ) );
00663 
00664     // block 0,2
00665     R[2] += A.column( 0 ) * transpose( B.column( 2 ) ) + B.column( 0 ) * transpose( A.column( 2 ) );
00666 
00667     // block 1,2
00668     R[4] += A.column( 1 ) * transpose( B.column( 2 ) ) + B.column( 1 ) * transpose( A.column( 2 ) );
00669 }
00670 void set_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6],
00671                                    double alpha,
00672                                    const MsqMatrix< 3, 3 >& A_in,
00673                                    const MsqMatrix< 3, 3 >& B )
00674 {
00675     // apply scalar first
00676     MsqMatrix< 3, 3 > A( A_in );
00677     A *= alpha;
00678 
00679     // block 0,0
00680     R[0] = A.column( 0 ) * transpose( B.column( 0 ) );
00681     R[0] += transpose( R[0] );
00682 
00683     // block 1,1
00684     R[3] = A.column( 1 ) * transpose( B.column( 1 ) );
00685     R[3] += transpose( R[3] );
00686 
00687     // block 2,2
00688     R[5] = A.column( 2 ) * transpose( B.column( 2 ) );
00689     R[5] += transpose( R[5] );
00690 
00691     // block 0,1
00692     R[1] = A.column( 0 ) * transpose( B.column( 1 ) );
00693     R[1] += B.column( 0 ) * transpose( A.column( 1 ) );
00694 
00695     // block 0,2
00696     R[2] = A.column( 0 ) * transpose( B.column( 2 ) );
00697     R[2] += B.column( 0 ) * transpose( A.column( 2 ) );
00698 
00699     // block 1,2
00700     R[4] = A.column( 1 ) * transpose( B.column( 2 ) );
00701     R[4] += B.column( 1 ) * transpose( A.column( 2 ) );
00702 }
00703 #endif
00704 
00705 #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
00706 void pluseq_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3],
00707                                       double alpha,
00708                                       const MsqMatrix< 2, 2 >& A_in,
00709                                       const MsqMatrix< 2, 2 >& B )
00710 {
00711     // apply scalar first
00712     MsqMatrix< 2, 2 > A( A_in ), tmp;
00713     A *= alpha;
00714 
00715     // block 0,0
00716     tmp = transpose( A.row( 0 ) ) * B.row( 0 );
00717     R[0] += tmp;
00718     R[0] += transpose( tmp );
00719 
00720     // block 1,1
00721     tmp = transpose( A.row( 1 ) ) * B.row( 1 );
00722     R[2] += tmp;
00723     R[2] += transpose( tmp );
00724 
00725     // block 0,1
00726     R[1] += transpose( A.row( 0 ) ) * B.row( 1 ) + transpose( B.row( 0 ) ) * A.row( 1 );
00727 }
00728 void set_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3],
00729                                    double alpha,
00730                                    const MsqMatrix< 2, 2 >& A_in,
00731                                    const MsqMatrix< 2, 2 >& B )
00732 {
00733     // apply scalar first
00734     MsqMatrix< 2, 2 > A( A_in );
00735     A *= alpha;
00736 
00737     // block 0,0
00738     R[0] = transpose( A.row( 0 ) ) * B.row( 0 );
00739     R[0] += transpose( R[0] );
00740 
00741     // block 1,1
00742     R[2] = transpose( A.row( 1 ) ) * B.row( 1 );
00743     R[2] += transpose( R[2] );
00744 
00745     // block 0,1
00746     R[1] = transpose( A.row( 0 ) ) * B.row( 1 );
00747     R[1] += transpose( B.row( 0 ) ) * A.row( 1 );
00748 }
00749 #else
00750 void pluseq_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3],
00751                                       double alpha,
00752                                       const MsqMatrix< 2, 2 >& A_in,
00753                                       const MsqMatrix< 2, 2 >& B )
00754 {
00755     // apply scalar first
00756     MsqMatrix< 2, 2 > A( A_in ), tmp;
00757     A *= alpha;
00758 
00759     // block 0,0
00760     tmp = A.column( 0 ) * transpose( B.column( 0 ) );
00761     R[0] += tmp;
00762     R[0] += transpose( tmp );
00763 
00764     // block 1,1
00765     tmp = A.column( 1 ) * transpose( B.column( 1 ) );
00766     R[2] += tmp;
00767     R[2] += transpose( tmp );
00768 
00769     // block 0,1
00770     R[1] += A.column( 0 ) * transpose( B.column( 1 ) ) + B.column( 0 ) * transpose( A.column( 1 ) );
00771 }
00772 void set_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3],
00773                                    double alpha,
00774                                    const MsqMatrix< 2, 2 >& A_in,
00775                                    const MsqMatrix< 2, 2 >& B )
00776 {
00777     // apply scalar first
00778     MsqMatrix< 2, 2 > A( A_in );
00779     A *= alpha;
00780 
00781     // block 0,0
00782     R[0] = A.column( 0 ) * transpose( B.column( 0 ) );
00783     R[0] += transpose( R[0] );
00784 
00785     // block 1,1
00786     R[2] = A.column( 1 ) * transpose( B.column( 1 ) );
00787     R[2] += transpose( R[2] );
00788 
00789     // block 0,1
00790     R[1] = A.column( 0 ) * transpose( B.column( 1 ) );
00791     R[1] += B.column( 0 ) * transpose( A.column( 1 ) );
00792 }
00793 #endif
00794 
00795 void pluseq_scaled_outer_product_I_I( MsqMatrix< 3, 3 > R[6], double alpha )
00796 {
00797     R[0]( 0, 0 ) += alpha;
00798     R[1]( 0, 1 ) += alpha;
00799     R[2]( 0, 2 ) += alpha;
00800     R[3]( 1, 1 ) += alpha;
00801     R[4]( 1, 2 ) += alpha;
00802     R[5]( 2, 2 ) += alpha;
00803 }
00804 
00805 void pluseq_scaled_outer_product_I_I( MsqMatrix< 2, 2 > R[3], double alpha )
00806 {
00807     R[0]( 0, 0 ) += alpha;
00808     R[1]( 0, 1 ) += alpha;
00809     R[2]( 1, 1 ) += alpha;
00810 }
00811 
00812 #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
00813 inline void pluseq_I_outer_product( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A )
00814 {
00815     R[0].add_row( 0, A.row( 0 ) );
00816     R[1].add_row( 0, A.row( 1 ) );
00817     R[2].add_row( 0, A.row( 2 ) );
00818     R[3].add_row( 1, A.row( 1 ) );
00819     R[4].add_row( 1, A.row( 2 ) );
00820     R[5].add_row( 2, A.row( 2 ) );
00821 }
00822 inline void pluseq_I_outer_product( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& A )
00823 {
00824     R[0].add_row( 0, A.row( 0 ) );
00825     R[1].add_row( 0, A.row( 1 ) );
00826     R[2].add_row( 1, A.row( 1 ) );
00827 }
00828 #else
00829 inline void pluseq_I_outer_product( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A )
00830 {
00831     R[0].add_row( 0, transpose( A.column( 0 ) ) );
00832     R[1].add_row( 0, transpose( A.column( 1 ) ) );
00833     R[2].add_row( 0, transpose( A.column( 2 ) ) );
00834     R[3].add_row( 1, transpose( A.column( 1 ) ) );
00835     R[4].add_row( 1, transpose( A.column( 2 ) ) );
00836     R[5].add_row( 2, transpose( A.column( 2 ) ) );
00837 }
00838 inline void pluseq_I_outer_product( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& A )
00839 {
00840     R[0].add_row( 0, transpose( A.column( 0 ) ) );
00841     R[1].add_row( 0, transpose( A.column( 1 ) ) );
00842     R[2].add_row( 1, transpose( A.column( 1 ) ) );
00843 }
00844 #endif
00845 
00846 #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
00847 inline void pluseq_outer_product_I( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A )
00848 {
00849     R[0].add_column( 0, transpose( A.row( 0 ) ) );
00850     R[1].add_column( 1, transpose( A.row( 0 ) ) );
00851     R[2].add_column( 2, transpose( A.row( 0 ) ) );
00852     R[3].add_column( 1, transpose( A.row( 1 ) ) );
00853     R[4].add_column( 2, transpose( A.row( 1 ) ) );
00854     R[5].add_column( 2, transpose( A.row( 2 ) ) );
00855 }
00856 inline void pluseq_outer_product_I( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& A )
00857 {
00858     R[0].add_column( 0, transpose( A.row( 0 ) ) );
00859     R[1].add_column( 1, transpose( A.row( 0 ) ) );
00860     R[2].add_column( 1, transpose( A.row( 1 ) ) );
00861 }
00862 #else
00863 inline void pluseq_outer_product_I( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A )
00864 {
00865     R[0].add_column( 0, A.column( 0 ) );
00866     R[1].add_column( 1, A.column( 0 ) );
00867     R[2].add_column( 2, A.column( 0 ) );
00868     R[3].add_column( 1, A.column( 1 ) );
00869     R[4].add_column( 2, A.column( 1 ) );
00870     R[5].add_column( 2, A.column( 2 ) );
00871 }
00872 inline void pluseq_outer_product_I( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& A )
00873 {
00874     R[0].add_column( 0, A.column( 0 ) );
00875     R[1].add_column( 1, A.column( 0 ) );
00876     R[2].add_column( 1, A.column( 1 ) );
00877 }
00878 #endif
00879 
00880 template < unsigned D >
00881 void pluseq_scaled_sum_outer_product_I_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2],
00882                                           double alpha,
00883                                           const MsqMatrix< D, D >& A_in )
00884 {
00885     // apply scalar first
00886     MsqMatrix< D, D > A( A_in );
00887     A *= alpha;
00888     pluseq_I_outer_product( R, A );
00889     pluseq_outer_product_I( R, A );
00890 }
00891 
00892 template < unsigned D >
00893 void second_deriv_wrt_product_factor_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], const MsqMatrix< D, D >& Z )
00894 {
00895     const MsqMatrix< D, D > Zt = transpose( Z );
00896     for( unsigned i = 0; i < D * ( D + 1 ) / 2; ++i )
00897         R[i] = ( Z * R[i] ) * Zt;
00898 }
00899 
00900 void set_scaled_2nd_deriv_wrt_psi( MsqMatrix< 2, 2 > R[3],
00901                                    const double alpha,
00902                                    const double psi,
00903                                    const MsqMatrix< 2, 2 >& T )
00904 {
00905     const double t = trace( T );
00906     const double f = alpha / ( psi * psi * psi );
00907     const double s = T( 0, 1 ) - T( 1, 0 );
00908     R[0]( 0, 0 ) = R[1]( 0, 1 ) = R[2]( 1, 1 ) = f * s * s;
00909     R[0]( 0, 1 ) = R[0]( 1, 0 ) = R[1]( 1, 1 ) = -f * s * t;
00910     R[1]( 0, 0 ) = R[2]( 0, 1 ) = R[2]( 1, 0 ) = f * s * t;
00911     R[0]( 1, 1 ) = R[2]( 0, 0 ) = -( R[1]( 1, 0 ) = -f * t * t );
00912 }
00913 
00914 template < unsigned D >
00915 inline void pluseq_scaled( MsqMatrix< D, D > R[D * ( D + 1 ) / 2],
00916                            double alpha,
00917                            const MsqMatrix< D, D > Z[D * ( D + 1 ) / 2] )
00918 {
00919     for( unsigned i = 0; i < D * ( D + 1 ) / 2; ++i )
00920         R[i] += alpha * Z[i];
00921 }
00922 
00923 /**\brief \f$ R *= s */
00924 template < unsigned D >
00925 inline void hess_scale_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha )
00926 {
00927     for( unsigned i = 0; i < D * ( D + 1 ) / 2; ++i )
00928         R[i] *= alpha;
00929 }
00930 
00931 }  // namespace MBMesquite
00932 
00933 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines