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