MOAB: Mesh Oriented datABase  (version 5.4.0)
TQualityMetric.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 TQualityMetric.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #undef PRINT_INFO
00033 
00034 #include "Mesquite.hpp"
00035 #include "TQualityMetric.hpp"
00036 #include "MsqMatrix.hpp"
00037 #include "ElementQM.hpp"
00038 #include "MsqError.hpp"
00039 #include "Vector3D.hpp"
00040 #include "PatchData.hpp"
00041 #include "MappingFunction.hpp"
00042 #include "WeightCalculator.hpp"
00043 #include "TargetCalculator.hpp"
00044 #include "TMetric.hpp"
00045 #include "TMetricBarrier.hpp"
00046 #include "TargetMetricUtil.hpp"
00047 #include "TMPDerivs.hpp"
00048 
00049 #ifdef PRINT_INFO
00050 #include <iostream>
00051 #endif
00052 
00053 #include <functional>
00054 #include <algorithm>
00055 
00056 #define NUMERICAL_2D_HESSIAN
00057 
00058 namespace MBMesquite
00059 {
00060 
00061 std::string TQualityMetric::get_name() const
00062 {
00063     return targetMetric->get_name();
00064 }
00065 
00066 bool TQualityMetric::evaluate_internal( PatchData& pd,
00067                                         size_t p_handle,
00068                                         double& value,
00069                                         size_t* indices,
00070                                         size_t& num_indices,
00071                                         MsqError& err )
00072 {
00073     const Sample s        = ElemSampleQM::sample( p_handle );
00074     const size_t e        = ElemSampleQM::elem( p_handle );
00075     MsqMeshEntity& p_elem = pd.element_by_index( e );
00076     EntityTopology type   = p_elem.get_element_type();
00077     unsigned edim         = TopologyInfo::dimension( type );
00078     const NodeSet bits    = pd.non_slave_node_set( e );
00079 
00080     bool rval;
00081     if( edim == 3 )
00082     {  // 3x3 or 3x2 targets ?
00083         const MappingFunction3D* mf = pd.get_mapping_function_3D( type );
00084         if( !mf )
00085         {
00086             MSQ_SETERR( err )
00087             ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
00088             return false;
00089         }
00090 
00091         MsqMatrix< 3, 3 > A, W;
00092         mf->jacobian( pd, e, bits, s, indices, mDerivs3D, num_indices, A, err );
00093         MSQ_ERRZERO( err );
00094         targetCalc->get_3D_target( pd, e, s, W, err );
00095         MSQ_ERRZERO( err );
00096         const MsqMatrix< 3, 3 > Winv = inverse( W );
00097         const MsqMatrix< 3, 3 > T    = A * Winv;
00098         rval                         = targetMetric->evaluate( T, value, err );
00099         MSQ_ERRZERO( err );
00100 #ifdef PRINT_INFO
00101         print_info< 3 >( e, s, A, W, A * inverse( W ) );
00102 #endif
00103     }
00104     else if( edim == 2 )
00105     {
00106         MsqMatrix< 2, 2 > W, A;
00107         MsqMatrix< 3, 2 > S_a_transpose_Theta;
00108         rval =
00109             evaluate_surface_common( pd, s, e, bits, indices, num_indices, mDerivs2D, W, A, S_a_transpose_Theta, err );
00110         if( MSQ_CHKERR( err ) || !rval ) return false;
00111         const MsqMatrix< 2, 2 > Winv = inverse( W );
00112         const MsqMatrix< 2, 2 > T    = A * Winv;
00113         rval                         = targetMetric->evaluate( T, value, err );
00114         MSQ_ERRZERO( err );
00115 
00116 #ifdef PRINT_INFO
00117         print_info< 2 >( e, s, J, Wp, A * inverse( W ) );
00118 #endif
00119     }
00120     else
00121     {
00122         assert( false );
00123         return false;
00124     }
00125 
00126     return rval;
00127 }
00128 
00129 bool TQualityMetric::evaluate_with_gradient( PatchData& pd,
00130                                              size_t p_handle,
00131                                              double& value,
00132                                              std::vector< size_t >& indices,
00133                                              std::vector< Vector3D >& grad,
00134                                              MsqError& err )
00135 {
00136     const Sample s        = ElemSampleQM::sample( p_handle );
00137     const size_t e        = ElemSampleQM::elem( p_handle );
00138     MsqMeshEntity& p_elem = pd.element_by_index( e );
00139     EntityTopology type   = p_elem.get_element_type();
00140     unsigned edim         = TopologyInfo::dimension( type );
00141     size_t num_idx        = 0;
00142     const NodeSet bits    = pd.non_slave_node_set( e );
00143 
00144     bool rval;
00145     if( edim == 3 )
00146     {  // 3x3 or 3x2 targets ?
00147         const MappingFunction3D* mf = pd.get_mapping_function_3D( type );
00148         if( !mf )
00149         {
00150             MSQ_SETERR( err )
00151             ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
00152             return false;
00153         }
00154 
00155         MsqMatrix< 3, 3 > A, W, dmdT;
00156         mf->jacobian( pd, e, bits, s, mIndices, mDerivs3D, num_idx, A, err );
00157         MSQ_ERRZERO( err );
00158         targetCalc->get_3D_target( pd, e, s, W, err );
00159         MSQ_ERRZERO( err );
00160         const MsqMatrix< 3, 3 > Winv = inverse( W );
00161         const MsqMatrix< 3, 3 > T    = A * Winv;
00162         rval                         = targetMetric->evaluate_with_grad( T, value, dmdT, err );
00163         MSQ_ERRZERO( err );
00164         gradient< 3 >( num_idx, mDerivs3D, dmdT * transpose( Winv ), grad );
00165 #ifdef PRINT_INFO
00166         print_info< 3 >( e, s, A, W, A * inverse( W ) );
00167 #endif
00168     }
00169     else if( edim == 2 )
00170     {
00171         MsqMatrix< 2, 2 > W, A, dmdT;
00172         MsqMatrix< 3, 2 > S_a_transpose_Theta;
00173         rval = evaluate_surface_common( pd, s, e, bits, mIndices, num_idx, mDerivs2D, W, A, S_a_transpose_Theta, err );
00174         if( MSQ_CHKERR( err ) || !rval ) return false;
00175         const MsqMatrix< 2, 2 > Winv = inverse( W );
00176         const MsqMatrix< 2, 2 > T    = A * Winv;
00177         rval                         = targetMetric->evaluate_with_grad( T, value, dmdT, err );
00178         MSQ_ERRZERO( err );
00179         gradient< 2 >( num_idx, mDerivs2D, S_a_transpose_Theta * dmdT * transpose( Winv ), grad );
00180 #ifdef PRINT_INFO
00181         print_info< 2 >( e, s, J, Wp, A * inverse( W ) );
00182 #endif
00183     }
00184     else
00185     {
00186         assert( false );
00187         return false;
00188     }
00189 
00190     // pass back index list
00191     indices.resize( num_idx );
00192     std::copy( mIndices, mIndices + num_idx, indices.begin() );
00193 
00194     // apply target weight to value
00195     weight( pd, s, e, num_idx, value, grad.empty() ? 0 : arrptr( grad ), 0, 0, err );
00196     MSQ_ERRZERO( err );
00197     return rval;
00198 }
00199 
00200 bool TQualityMetric::evaluate_with_Hessian( PatchData& pd,
00201                                             size_t p_handle,
00202                                             double& value,
00203                                             std::vector< size_t >& indices,
00204                                             std::vector< Vector3D >& grad,
00205                                             std::vector< Matrix3D >& Hessian,
00206                                             MsqError& err )
00207 {
00208     const Sample s        = ElemSampleQM::sample( p_handle );
00209     const size_t e        = ElemSampleQM::elem( p_handle );
00210     MsqMeshEntity& p_elem = pd.element_by_index( e );
00211     EntityTopology type   = p_elem.get_element_type();
00212     unsigned edim         = TopologyInfo::dimension( type );
00213     size_t num_idx        = 0;
00214     const NodeSet bits    = pd.non_slave_node_set( e );
00215 
00216     bool rval;
00217     if( edim == 3 )
00218     {  // 3x3 or 3x2 targets ?
00219         const MappingFunction3D* mf = pd.get_mapping_function_3D( type );
00220         if( !mf )
00221         {
00222             MSQ_SETERR( err )
00223             ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
00224             return false;
00225         }
00226 
00227         MsqMatrix< 3, 3 > A, W, dmdT, d2mdT2[6];
00228         mf->jacobian( pd, e, bits, s, mIndices, mDerivs3D, num_idx, A, err );
00229         MSQ_ERRZERO( err );
00230         targetCalc->get_3D_target( pd, e, s, W, err );
00231         MSQ_ERRZERO( err );
00232         const MsqMatrix< 3, 3 > Winv = inverse( W );
00233         const MsqMatrix< 3, 3 > T    = A * Winv;
00234         rval                         = targetMetric->evaluate_with_hess( T, value, dmdT, d2mdT2, err );
00235         MSQ_ERRZERO( err );
00236         gradient< 3 >( num_idx, mDerivs3D, dmdT * transpose( Winv ), grad );
00237         second_deriv_wrt_product_factor( d2mdT2, Winv );
00238         Hessian.resize( num_idx * ( num_idx + 1 ) / 2 );
00239         if( num_idx ) hessian< 3 >( num_idx, mDerivs3D, d2mdT2, arrptr( Hessian ) );
00240 
00241 #ifdef PRINT_INFO
00242         print_info< 3 >( e, s, A, W, A * inverse( W ) );
00243 #endif
00244     }
00245     else if( edim == 2 )
00246     {
00247 #ifdef NUMERICAL_2D_HESSIAN
00248         // return finite difference approximation for now
00249 
00250         return QualityMetric::evaluate_with_Hessian( pd, p_handle, value, indices, grad, Hessian, err );
00251 #else
00252         MsqMatrix< 2, 2 > W, A, dmdT, d2mdT2[3];
00253         MsqMatrix< 3, 2 > M;
00254         rval = evaluate_surface_common( pd, s, e, bits, mIndices, num_idx, mDerivs2D, W, A, M, err );
00255         if( MSQ_CHKERR( err ) || !rval ) return false;
00256         const MsqMatrix< 2, 2 > Winv = inverse( W );
00257         const MsqMatrix< 2, 2 > T    = A * Winv;
00258         rval                         = targetMetric->evaluate_with_hess( T, value, dmdT, d2mdT2, err );
00259         MSQ_ERRZERO( err );
00260         gradient< 2 >( num_idx, mDerivs2D, M * dmdT * transpose( Winv ), grad );
00261         // calculate 2D hessian
00262         second_deriv_wrt_product_factor( d2mdT2, Winv );
00263         const size_t n = num_idx * ( num_idx + 1 ) / 2;
00264         hess2d.resize( n );
00265         if( n ) hessian< 2 >( num_idx, mDerivs2D, d2mdT2, arrptr( hess2d ) );
00266         // calculate surface hessian as transform of 2D hessian
00267         Hessian.resize( n );
00268         for( size_t i = 0; i < n; ++i )
00269             Hessian[i] = Matrix3D( ( M * hess2d[i] * transpose( M ) ).data() );
00270 #ifdef PRINT_INFO
00271         print_info< 2 >( e, s, J, Wp, A * inverse( W ) );
00272 #endif
00273 #endif
00274     }
00275     else
00276     {
00277         assert( 0 );
00278         return false;
00279     }
00280 
00281     // pass back index list
00282     indices.resize( num_idx );
00283     std::copy( mIndices, mIndices + num_idx, indices.begin() );
00284 
00285     // apply target weight to value
00286     if( !num_idx )
00287         weight( pd, s, e, num_idx, value, 0, 0, 0, err );
00288     else
00289         weight( pd, s, e, num_idx, value, arrptr( grad ), 0, arrptr( Hessian ), err );
00290     MSQ_ERRZERO( err );
00291     return rval;
00292 }
00293 
00294 bool TQualityMetric::evaluate_with_Hessian_diagonal( PatchData& pd,
00295                                                      size_t p_handle,
00296                                                      double& value,
00297                                                      std::vector< size_t >& indices,
00298                                                      std::vector< Vector3D >& grad,
00299                                                      std::vector< SymMatrix3D >& diagonal,
00300                                                      MsqError& err )
00301 {
00302     const Sample s        = ElemSampleQM::sample( p_handle );
00303     const size_t e        = ElemSampleQM::elem( p_handle );
00304     MsqMeshEntity& p_elem = pd.element_by_index( e );
00305     EntityTopology type   = p_elem.get_element_type();
00306     unsigned edim         = TopologyInfo::dimension( type );
00307     size_t num_idx        = 0;
00308     const NodeSet bits    = pd.non_slave_node_set( e );
00309 
00310     bool rval;
00311     if( edim == 3 )
00312     {  // 3x3 or 3x2 targets ?
00313         const MappingFunction3D* mf = pd.get_mapping_function_3D( type );
00314         if( !mf )
00315         {
00316             MSQ_SETERR( err )
00317             ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
00318             return false;
00319         }
00320 
00321         MsqMatrix< 3, 3 > A, W, dmdT, d2mdT2[6];
00322         mf->jacobian( pd, e, bits, s, mIndices, mDerivs3D, num_idx, A, err );
00323         MSQ_ERRZERO( err );
00324         targetCalc->get_3D_target( pd, e, s, W, err );
00325         MSQ_ERRZERO( err );
00326         const MsqMatrix< 3, 3 > Winv = inverse( W );
00327         const MsqMatrix< 3, 3 > T    = A * Winv;
00328         rval                         = targetMetric->evaluate_with_hess( T, value, dmdT, d2mdT2, err );
00329         MSQ_ERRZERO( err );
00330         gradient< 3 >( num_idx, mDerivs3D, dmdT * transpose( Winv ), grad );
00331         second_deriv_wrt_product_factor( d2mdT2, Winv );
00332 
00333         diagonal.resize( num_idx );
00334         hessian_diagonal< 3 >( num_idx, mDerivs3D, d2mdT2, arrptr( diagonal ) );
00335 #ifdef PRINT_INFO
00336         print_info< 3 >( e, s, A, W, A * inverse( W ) );
00337 #endif
00338     }
00339     else if( edim == 2 )
00340     {
00341 #ifdef NUMERICAL_2D_HESSIAN
00342         // use finite diference approximation for now
00343         return QualityMetric::evaluate_with_Hessian_diagonal( pd, p_handle, value, indices, grad, diagonal, err );
00344 #else
00345         MsqMatrix< 2, 2 > W, A, dmdT, d2mdT2[3];
00346         MsqMatrix< 3, 2 > M;
00347         rval = evaluate_surface_common( pd, s, e, bits, mIndices, num_idx, mDerivs2D, W, A, M, err );
00348         if( MSQ_CHKERR( err ) || !rval ) return false;
00349         const MsqMatrix< 2, 2 > Winv = inverse( W );
00350         const MsqMatrix< 2, 2 > T    = A * Winv;
00351         rval                         = targetMetric->evaluate_with_hess( T, value, dmdT, d2mdT2, err );
00352         MSQ_ERRZERO( err );
00353         gradient< 2 >( num_idx, mDerivs2D, M * dmdT * transpose( Winv ), grad );
00354         second_deriv_wrt_product_factor( d2mdT2, Winv );
00355 
00356         diagonal.resize( num_idx );
00357         for( size_t i = 0; i < num_idx; ++i )
00358         {
00359             MsqMatrix< 2, 2 > block2d;
00360             block2d( 0, 0 )     = transpose( mDerivs2D[i] ) * d2mdT2[0] * mDerivs2D[i];
00361             block2d( 0, 1 )     = transpose( mDerivs2D[i] ) * d2mdT2[1] * mDerivs2D[i];
00362             block2d( 1, 0 )     = block2d( 0, 1 );
00363             block2d( 1, 1 )     = transpose( mDerivs2D[i] ) * d2mdT2[2] * mDerivs2D[i];
00364             MsqMatrix< 3, 2 > p = M * block2d;
00365 
00366             SymMatrix3D& H = diagonal[i];
00367             H[0]           = p.row( 0 ) * transpose( M.row( 0 ) );
00368             H[1]           = p.row( 0 ) * transpose( M.row( 1 ) );
00369             H[2]           = p.row( 0 ) * transpose( M.row( 2 ) );
00370             H[3]           = p.row( 1 ) * transpose( M.row( 1 ) );
00371             H[4]           = p.row( 1 ) * transpose( M.row( 2 ) );
00372             H[5]           = p.row( 2 ) * transpose( M.row( 2 ) );
00373         }
00374 #ifdef PRINT_INFO
00375         print_info< 2 >( e, s, J, Wp, A * inverse( W ) );
00376 #endif
00377 #endif
00378     }
00379     else
00380     {
00381         assert( 0 );
00382         return false;
00383     }
00384 
00385     // pass back index list
00386     indices.resize( num_idx );
00387     std::copy( mIndices, mIndices + num_idx, indices.begin() );
00388 
00389     // apply target weight to value
00390     if( !num_idx )
00391         weight( pd, s, e, num_idx, value, 0, 0, 0, err );
00392     else
00393         weight( pd, s, e, num_idx, value, arrptr( grad ), arrptr( diagonal ), 0, err );
00394     MSQ_ERRZERO( err );
00395     return rval;
00396 }
00397 
00398 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines