MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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) kraftche@cae.wisc.edu 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 <assert.h> 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, const MsqMatrix< 3, 2 >& W_32, MsqMatrix< 2, 2 >& W_22, 00065 MsqMatrix< 3, 2 >& RZ ); 00066 /* 00067 void surface_to_2d( const MsqMatrix<3,2>& A_in, 00068 const MsqMatrix<3,2>& W_in, 00069 MsqMatrix<2,2>& A_out, 00070 MsqMatrix<2,2>& W_out ); 00071 */ 00072 void get_sample_pt_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free, MsqError& err ); 00073 00074 void get_elem_sample_points( PatchData& pd, size_t elem, std::vector< size_t >& handles, MsqError& err ); 00075 00076 /**\brief Calculate gradient from derivatives of mapping function terms 00077 * and derivatives of target metric. */ 00078 template < int DIM > 00079 inline void gradient( size_t num_free_verts, const MsqVector< DIM >* dNdxi, const MsqMatrix< 3, DIM >& dmdA, 00080 std::vector< Vector3D >& grad ) 00081 { 00082 grad.clear(); 00083 grad.resize( num_free_verts, Vector3D( 0, 0, 0 ) ); 00084 for( size_t i = 0; i < num_free_verts; ++i ) 00085 grad[i] = Vector3D( ( dmdA * dNdxi[i] ).data() ); 00086 } 00087 00088 /**\brief Calculate Hessian from derivatives of mapping function terms 00089 * and derivatives of target metric. */ 00090 template < int DIM, typename MAT > 00091 inline void hessian( size_t num_free_verts, const MsqVector< DIM >* dNdxi, const MsqMatrix< DIM, DIM >* d2mdA2, 00092 MAT* hess ) 00093 { 00094 MsqMatrix< 1, DIM > tmp[DIM][DIM]; 00095 size_t h = 0; // index of current Hessian block 00096 00097 for( size_t i = 0; i < num_free_verts; ++i ) 00098 { 00099 00100 // Populate TMP with vector-matrix procucts common 00101 // to terms of this Hessian row. 00102 const MsqMatrix< 1, DIM >& gi = transpose( dNdxi[i] ); 00103 switch( DIM ) 00104 { 00105 case 3: 00106 tmp[0][2] = gi * d2mdA2[2]; 00107 tmp[1][2] = gi * d2mdA2[4]; 00108 tmp[2][0] = gi * transpose( d2mdA2[2] ); 00109 tmp[2][1] = gi * transpose( d2mdA2[4] ); 00110 tmp[2][2] = gi * d2mdA2[5]; 00111 case 2: 00112 tmp[0][1] = gi * d2mdA2[1]; 00113 tmp[1][0] = gi * transpose( d2mdA2[1] ); 00114 tmp[1][1] = gi * d2mdA2[DIM]; 00115 case 1: 00116 tmp[0][0] = gi * d2mdA2[0]; 00117 case 0: 00118 break; 00119 default: 00120 assert( false ); 00121 } 00122 00123 // Calculate Hessian diagonal block 00124 MAT& H = hess[h++]; 00125 switch( DIM ) 00126 { 00127 case 3: 00128 H( 0, 2 ) = H( 2, 0 ) = tmp[0][2] * transpose( gi ); 00129 H( 1, 2 ) = H( 2, 1 ) = tmp[1][2] * transpose( gi ); 00130 H( 2, 2 ) = tmp[2][2] * transpose( gi ); 00131 case 2: 00132 H( 0, 1 ) = H( 1, 0 ) = tmp[0][1] * transpose( gi ); 00133 H( 1, 1 ) = tmp[1][1] * transpose( gi ); 00134 case 1: 00135 H( 0, 0 ) = tmp[0][0] * transpose( gi ); 00136 case 0: 00137 break; 00138 default: 00139 assert( false ); 00140 } 00141 00142 // Calculate remainder of Hessian row 00143 for( size_t j = i + 1; j < num_free_verts; ++j ) 00144 { 00145 MAT& HH = hess[h++]; 00146 const MsqMatrix< DIM, 1 >& gj = dNdxi[j]; 00147 switch( DIM ) 00148 { 00149 case 3: 00150 HH( 0, 2 ) = tmp[0][2] * gj; 00151 HH( 1, 2 ) = tmp[1][2] * gj; 00152 HH( 2, 0 ) = tmp[2][0] * gj; 00153 HH( 2, 1 ) = tmp[2][1] * gj; 00154 HH( 2, 2 ) = tmp[2][2] * gj; 00155 case 2: 00156 HH( 0, 1 ) = tmp[0][1] * gj; 00157 HH( 1, 0 ) = tmp[1][0] * gj; 00158 HH( 1, 1 ) = tmp[1][1] * gj; 00159 case 1: 00160 HH( 0, 0 ) = tmp[0][0] * gj; 00161 case 0: 00162 break; 00163 default: 00164 assert( false ); 00165 } 00166 } 00167 } 00168 } 00169 00170 /**\brief Calculate Hessian from derivatives of mapping function terms 00171 * and derivatives of target metric. */ 00172 template < int DIM > 00173 inline void hessian_diagonal( size_t num_free_verts, const MsqVector< DIM >* dNdxi, const MsqMatrix< DIM, DIM >* d2mdA2, 00174 SymMatrix3D* diagonal ) 00175 { 00176 for( size_t i = 0; i < num_free_verts; ++i ) 00177 { 00178 SymMatrix3D& H = diagonal[i]; 00179 for( unsigned j = 0; j < ( ( DIM ) * ( DIM + 1 ) / 2 ); ++j ) 00180 H[j] = transpose( dNdxi[i] ) * d2mdA2[j] * dNdxi[i]; 00181 } 00182 } 00183 00184 #ifdef PRINT_INFO 00185 template < int R, int C > 00186 inline void write_vect( char n, const MsqMatrix< R, C >& M ) 00187 { 00188 std::cout << " " << n << ':'; 00189 for( int c = 0; c < C; ++c ) 00190 { 00191 std::cout << '['; 00192 for( int r = 0; r < R; ++r ) 00193 std::cout << M( r, c ) << ' '; 00194 std::cout << ']'; 00195 } 00196 std::cout << std::endl; 00197 } 00198 00199 template < int D > 00200 inline void print_info( size_t elem, Sample sample, const MsqMatrix< 3, D >& A, const MsqMatrix< 3, D >& W, 00201 const MsqMatrix< D, D >& T ) 00202 { 00203 std::cout << "Elem " << elem << " Dim " << sample.dimension << " Num " << sample.number << " :" << std::endl; 00204 write_vect< 3, D >( 'A', A ); 00205 write_vect< 3, D >( 'W', W ); 00206 write_vect< D, D >( 'T', T ); 00207 } 00208 #endif 00209 00210 } // namespace MBMesquite 00211 00212 #endif