MOAB: Mesh Oriented datABase  (version 5.4.1)
TargetMetricUtil.hpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2007 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     (2007) [email protected]
00024 
00025   ***************************************************************** */
00026 
00027 /** \file TargetMetricUtil.hpp
00028  *  \brief A collection of utility code used by QualtiyMetrics
00029  *         composed of TMP Target Metrics
00030  *  \author Jason Kraftcheck
00031  */
00032 
00033 #ifndef MSQ_TARGET_METRIC_UTIL_HPP
00034 #define MSQ_TARGET_METRIC_UTIL_HPP
00035 
00036 #include "Mesquite.hpp"
00037 #include "SymMatrix3D.hpp"
00038 #include <vector>
00039 #include <cassert>
00040 
00041 namespace MBMesquite
00042 {
00043 
00044 template < unsigned R, unsigned C >
00045 class MsqMatrix;
00046 template < unsigned C >
00047 class MsqVector;
00048 class PatchData;
00049 class MsqError;
00050 class Vector3D;
00051 class Matrix3D;
00052 
00053 /**\brief Calculate R and Z such that \f$W\prime = Z^{-1} W\f$ and
00054  *        \f$A\prime = (RZ)^{-1} A\f$
00055  *
00056  * Calculate the matrices required to transform the active and target
00057  * matrices from the 3x2 surface domain to a 2x2 2D domain.
00058  *\param A    Input: Element Jacobian matrix.
00059  *\param W_32 Input: Target Jacobian matrix.
00060  *\param W_22 Output: 2D Target matrix.
00061  *\param RZ   Output: Product of R and Z needed to calculate the 2D
00062  *            element matrix.
00063  */
00064 void surface_to_2d( const MsqMatrix< 3, 2 >& A,
00065                     const MsqMatrix< 3, 2 >& W_32,
00066                     MsqMatrix< 2, 2 >& W_22,
00067                     MsqMatrix< 3, 2 >& RZ );
00068 /*
00069 void surface_to_2d( const MsqMatrix<3,2>& A_in,
00070                     const MsqMatrix<3,2>& W_in,
00071                     MsqMatrix<2,2>& A_out,
00072                     MsqMatrix<2,2>& W_out );
00073 */
00074 void get_sample_pt_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free, MsqError& err );
00075 
00076 void get_elem_sample_points( PatchData& pd, size_t elem, std::vector< size_t >& handles, MsqError& err );
00077 
00078 /**\brief Calculate gradient from derivatives of mapping function terms
00079  *        and derivatives of target metric. */
00080 template < int DIM >
00081 inline void gradient( size_t num_free_verts,
00082                       const MsqVector< DIM >* dNdxi,
00083                       const MsqMatrix< 3, DIM >& dmdA,
00084                       std::vector< Vector3D >& grad )
00085 {
00086     grad.clear();
00087     grad.resize( num_free_verts, Vector3D( 0, 0, 0 ) );
00088     for( size_t i = 0; i < num_free_verts; ++i )
00089         grad[i] = Vector3D( ( dmdA * dNdxi[i] ).data() );
00090 }
00091 
00092 /**\brief Calculate Hessian from derivatives of mapping function terms
00093  *        and derivatives of target metric. */
00094 template < int DIM, typename MAT >
00095 inline void hessian( size_t num_free_verts,
00096                      const MsqVector< DIM >* dNdxi,
00097                      const MsqMatrix< DIM, DIM >* d2mdA2,
00098                      MAT* hess )
00099 {
00100     MsqMatrix< 1, DIM > tmp[DIM][DIM];
00101     size_t h = 0;  // index of current Hessian block
00102 
00103     for( size_t i = 0; i < num_free_verts; ++i )
00104     {
00105 
00106         // Populate TMP with vector-matrix procucts common
00107         // to terms of this Hessian row.
00108         const MsqMatrix< 1, DIM >& gi = transpose( dNdxi[i] );
00109         switch( DIM )
00110         {
00111             case 3:
00112                 tmp[0][2] = gi * d2mdA2[2];
00113                 tmp[1][2] = gi * d2mdA2[4];
00114                 tmp[2][0] = gi * transpose( d2mdA2[2] );
00115                 tmp[2][1] = gi * transpose( d2mdA2[4] );
00116                 tmp[2][2] = gi * d2mdA2[5];
00117             case 2:
00118                 tmp[0][1] = gi * d2mdA2[1];
00119                 tmp[1][0] = gi * transpose( d2mdA2[1] );
00120                 tmp[1][1] = gi * d2mdA2[DIM];
00121             case 1:
00122                 tmp[0][0] = gi * d2mdA2[0];
00123             case 0:
00124                 break;
00125             default:
00126                 assert( false );
00127         }
00128 
00129         // Calculate Hessian diagonal block
00130         MAT& H = hess[h++];
00131         switch( DIM )
00132         {
00133             case 3:
00134                 H( 0, 2 ) = H( 2, 0 ) = tmp[0][2] * transpose( gi );
00135                 H( 1, 2 ) = H( 2, 1 ) = tmp[1][2] * transpose( gi );
00136                 H( 2, 2 )             = tmp[2][2] * transpose( gi );
00137             case 2:
00138                 H( 0, 1 ) = H( 1, 0 ) = tmp[0][1] * transpose( gi );
00139                 H( 1, 1 )             = tmp[1][1] * transpose( gi );
00140             case 1:
00141                 H( 0, 0 ) = tmp[0][0] * transpose( gi );
00142             case 0:
00143                 break;
00144             default:
00145                 assert( false );
00146         }
00147 
00148         // Calculate remainder of Hessian row
00149         for( size_t j = i + 1; j < num_free_verts; ++j )
00150         {
00151             MAT& HH                       = hess[h++];
00152             const MsqMatrix< DIM, 1 >& gj = dNdxi[j];
00153             switch( DIM )
00154             {
00155                 case 3:
00156                     HH( 0, 2 ) = tmp[0][2] * gj;
00157                     HH( 1, 2 ) = tmp[1][2] * gj;
00158                     HH( 2, 0 ) = tmp[2][0] * gj;
00159                     HH( 2, 1 ) = tmp[2][1] * gj;
00160                     HH( 2, 2 ) = tmp[2][2] * gj;
00161                 case 2:
00162                     HH( 0, 1 ) = tmp[0][1] * gj;
00163                     HH( 1, 0 ) = tmp[1][0] * gj;
00164                     HH( 1, 1 ) = tmp[1][1] * gj;
00165                 case 1:
00166                     HH( 0, 0 ) = tmp[0][0] * gj;
00167                 case 0:
00168                     break;
00169                 default:
00170                     assert( false );
00171             }
00172         }
00173     }
00174 }
00175 
00176 /**\brief Calculate Hessian from derivatives of mapping function terms
00177  *        and derivatives of target metric. */
00178 template < int DIM >
00179 inline void hessian_diagonal( size_t num_free_verts,
00180                               const MsqVector< DIM >* dNdxi,
00181                               const MsqMatrix< DIM, DIM >* d2mdA2,
00182                               SymMatrix3D* diagonal )
00183 {
00184     for( size_t i = 0; i < num_free_verts; ++i )
00185     {
00186         SymMatrix3D& H = diagonal[i];
00187         for( unsigned j = 0; j < ( ( DIM ) * ( DIM + 1 ) / 2 ); ++j )
00188             H[j] = transpose( dNdxi[i] ) * d2mdA2[j] * dNdxi[i];
00189     }
00190 }
00191 
00192 #ifdef PRINT_INFO
00193 template < int R, int C >
00194 inline void write_vect( char n, const MsqMatrix< R, C >& M )
00195 {
00196     std::cout << "  " << n << ':';
00197     for( int c = 0; c < C; ++c )
00198     {
00199         std::cout << '[';
00200         for( int r = 0; r < R; ++r )
00201             std::cout << M( r, c ) << ' ';
00202         std::cout << ']';
00203     }
00204     std::cout << std::endl;
00205 }
00206 
00207 template < int D >
00208 inline void print_info( size_t elem,
00209                         Sample sample,
00210                         const MsqMatrix< 3, D >& A,
00211                         const MsqMatrix< 3, D >& W,
00212                         const MsqMatrix< D, D >& T )
00213 {
00214     std::cout << "Elem " << elem << " Dim " << sample.dimension << " Num " << sample.number << " :" << std::endl;
00215     write_vect< 3, D >( 'A', A );
00216     write_vect< 3, D >( 'W', W );
00217     write_vect< D, D >( 'T', T );
00218 }
00219 #endif
00220 
00221 }  // namespace MBMesquite
00222 
00223 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines