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