MOAB: Mesh Oriented datABase  (version 5.2.1)
PMeanPMetric.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 PMeanPMetric.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "PMeanPMetric.hpp"
00034 #include "MsqError.hpp"
00035 #include "QualityMetric.hpp"
00036 #include "Vector3D.hpp"
00037 #include "Matrix3D.hpp"
00038 #include "SymMatrix3D.hpp"
00039 #include "PatchData.hpp"
00040 
00041 namespace MBMesquite
00042 {
00043 
00044 bool PMeanPMetric::average( PatchData& pd, QualityMetric* metric, const std::vector< size_t >& qm_handles,
00045                             double& value, MsqError& err )
00046 {
00047     bool rval = true;
00048     value     = 0;
00049     for( std::vector< size_t >::const_iterator i = qm_handles.begin(); i != qm_handles.end(); ++i )
00050     {
00051         double mval;
00052         if( !metric->evaluate( pd, *i, mval, err ) ) rval = false;
00053         MSQ_ERRZERO( err );
00054         value += P.raise( mval );
00055     }
00056     value /= qm_handles.size();
00057     return rval;
00058 }
00059 
00060 bool PMeanPMetric::average_with_indices( PatchData& pd, QualityMetric* metric, const std::vector< size_t >& qm_handles,
00061                                          double& value, std::vector< size_t >& indices, MsqError& err )
00062 {
00063     indices.clear();
00064 
00065     bool rval = true;
00066     value     = 0;
00067     for( std::vector< size_t >::const_iterator i = qm_handles.begin(); i != qm_handles.end(); ++i )
00068     {
00069         double mval;
00070         mIndices.clear();
00071         if( !metric->evaluate_with_indices( pd, *i, mval, mIndices, err ) ) rval = false;
00072         MSQ_ERRZERO( err );
00073         value += P.raise( mval );
00074 
00075         std::copy( mIndices.begin(), mIndices.end(), std::back_inserter( indices ) );
00076     }
00077     std::sort( indices.begin(), indices.end() );
00078     indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
00079 
00080     value /= qm_handles.size();
00081     return rval;
00082 }
00083 
00084 bool PMeanPMetric::average_with_gradient( PatchData& pd, QualityMetric* metric, const std::vector< size_t >& qm_handles,
00085                                           double& value, std::vector< size_t >& indices,
00086                                           std::vector< Vector3D >& gradient, MsqError& err )
00087 {
00088     indices.clear();
00089     gradient.clear();
00090 
00091     std::vector< Vector3D >::iterator g;
00092     std::vector< size_t >::iterator j, k;
00093     std::vector< size_t >::const_iterator i;
00094 
00095     bool rval = true;
00096     value     = 0;
00097     for( i = qm_handles.begin(); i != qm_handles.end(); ++i )
00098     {
00099 
00100         double mval;
00101         mIndices.clear();
00102         mGrad.clear();
00103         if( !metric->evaluate_with_gradient( pd, *i, mval, mIndices, mGrad, err ) ) rval = false;
00104         MSQ_ERRZERO( err );
00105         value += P.raise( mval );
00106 
00107         double p1val = P.value() * P1.raise( mval );
00108         for( j = mIndices.begin(), g = mGrad.begin(); j != mIndices.end(); ++j, ++g )
00109         {
00110 
00111             *g *= p1val;
00112             k = std::lower_bound( indices.begin(), indices.end(), *j );
00113             if( k == indices.end() || *k != *j )
00114             {
00115                 k          = indices.insert( k, *j );
00116                 size_t idx = k - indices.begin();
00117                 gradient.insert( gradient.begin() + idx, *g );
00118             }
00119             else
00120             {
00121                 size_t idx = k - indices.begin();
00122                 gradient[idx] += *g;
00123             }
00124         }
00125     }
00126 
00127     double inv_n = 1.0 / qm_handles.size();
00128     value *= inv_n;
00129     for( g = gradient.begin(); g != gradient.end(); ++g )
00130         *g *= inv_n;
00131 
00132     return rval;
00133 }
00134 
00135 bool PMeanPMetric::average_with_Hessian_diagonal( PatchData& pd, QualityMetric* metric,
00136                                                   const std::vector< size_t >& qm_handles, double& value,
00137                                                   std::vector< size_t >& indices, std::vector< Vector3D >& gradient,
00138                                                   std::vector< SymMatrix3D >& diagonal, MsqError& err )
00139 {
00140     // clear temporary storage
00141     mIndices.clear();
00142     mGrad.clear();
00143     mDiag.clear();
00144     mOffsets.clear();
00145     mValues.clear();
00146 
00147     // Evaluate metric for all sample points,
00148     // accumulating indices, gradients, and Hessians
00149     bool rval = true;
00150     std::vector< size_t >::const_iterator q;
00151     for( q = qm_handles.begin(); q != qm_handles.end(); ++q )
00152     {
00153         double mval;
00154         indices.clear();
00155         gradient.clear();
00156         diagonal.clear();
00157         if( !metric->evaluate_with_Hessian_diagonal( pd, *q, mval, indices, gradient, diagonal, err ) ) rval = false;
00158         MSQ_ERRZERO( err );
00159 
00160         mValues.push_back( mval );
00161         mOffsets.push_back( mIndices.size() );
00162         std::copy( indices.begin(), indices.end(), std::back_inserter( mIndices ) );
00163         std::copy( gradient.begin(), gradient.end(), std::back_inserter( mGrad ) );
00164         std::copy( diagonal.begin(), diagonal.end(), std::back_inserter( mDiag ) );
00165     }
00166     mOffsets.push_back( mIndices.size() );
00167 
00168     // Combine lists of free vertex indices, and update indices
00169     // in per-evaluation lists to point into the combined gradient
00170     // and Hessian lists.
00171     indices = mIndices;
00172     std::sort( indices.begin(), indices.end() );
00173     indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
00174     std::vector< size_t >::iterator i, j;
00175     for( i = mIndices.begin(); i != mIndices.end(); ++i )
00176     {
00177         j = std::lower_bound( indices.begin(), indices.end(), *i );
00178         assert( *j == *i );
00179         *i = j - indices.begin();
00180     }
00181 
00182     // Allocate space and zero output gradient and Hessian lists
00183     const size_t n     = indices.size();
00184     const size_t m     = mValues.size();
00185     const double inv_n = 1.0 / m;
00186     double v;
00187     gradient.clear();
00188     gradient.resize( n, Vector3D( 0, 0, 0 ) );
00189     diagonal.clear();
00190     diagonal.resize( n, SymMatrix3D( 0.0 ) );
00191 
00192     // Average values, gradients, Hessians
00193     value = 0.0;
00194     for( size_t k = 0; k < m; ++k )
00195     {  // for each metric evaluate
00196         if( P.value() == 1.0 )
00197         {
00198             v = P.raise( mValues[k] );
00199             // for each gradient (or Hessian) for the local metric evaluation
00200             for( size_t r = mOffsets[k]; r < mOffsets[k + 1]; ++r )
00201             {
00202                 const size_t nr = mIndices[r];
00203 
00204                 mDiag[r] *= inv_n;
00205                 diagonal[nr] += mDiag[r];
00206                 mGrad[r] *= inv_n;
00207                 gradient[nr] += mGrad[r];
00208             }
00209         }
00210         else
00211         {
00212             const double r2 = P2.raise( mValues[k] );
00213             const double r1 = r2 * mValues[k];
00214             v               = r1 * mValues[k];
00215             const double h  = r2 * P.value() * ( P.value() - 1 ) * inv_n;
00216             const double g  = r1 * P.value() * inv_n;
00217             // for each gradient (or Hessian) for the local metric evaluation
00218             for( size_t r = mOffsets[k]; r < mOffsets[k + 1]; ++r )
00219             {
00220                 const size_t nr = mIndices[r];
00221 
00222                 mDiag[r] *= g;
00223                 diagonal[nr] += mDiag[r];
00224                 diagonal[nr] += h * outer( mGrad[r] );
00225                 mGrad[r] *= g;
00226                 gradient[nr] += mGrad[r];
00227             }
00228         }
00229         value += v;
00230     }
00231 
00232     value *= inv_n;
00233 
00234     return rval;
00235 }
00236 
00237 bool PMeanPMetric::average_with_Hessian( PatchData& pd, QualityMetric* metric, const std::vector< size_t >& qm_handles,
00238                                          double& value, std::vector< size_t >& indices,
00239                                          std::vector< Vector3D >& gradient, std::vector< Matrix3D >& Hessian,
00240                                          MsqError& err )
00241 {
00242     // clear temporary storage
00243     mIndices.clear();
00244     mGrad.clear();
00245     mHess.clear();
00246     mOffsets.clear();
00247     mValues.clear();
00248 
00249     // Evaluate metric for all sample points,
00250     // accumulating indices, gradients, and Hessians
00251     bool rval = true;
00252     std::vector< size_t >::const_iterator q;
00253     for( q = qm_handles.begin(); q != qm_handles.end(); ++q )
00254     {
00255         double mval;
00256         indices.clear();
00257         gradient.clear();
00258         Hessian.clear();
00259         if( !metric->evaluate_with_Hessian( pd, *q, mval, indices, gradient, Hessian, err ) ) rval = false;
00260         MSQ_ERRZERO( err );
00261 
00262         mValues.push_back( mval );
00263         mOffsets.push_back( mIndices.size() );
00264         std::copy( indices.begin(), indices.end(), std::back_inserter( mIndices ) );
00265         std::copy( gradient.begin(), gradient.end(), std::back_inserter( mGrad ) );
00266         std::copy( Hessian.begin(), Hessian.end(), std::back_inserter( mHess ) );
00267     }
00268     mOffsets.push_back( mIndices.size() );
00269 
00270     // Combine lists of free vertex indices, and update indices
00271     // in per-evaluation lists to point into the combined gradient
00272     // and Hessian lists.
00273     indices = mIndices;
00274     std::sort( indices.begin(), indices.end() );
00275     indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
00276     std::vector< size_t >::iterator i, j;
00277     for( i = mIndices.begin(); i != mIndices.end(); ++i )
00278     {
00279         j = std::lower_bound( indices.begin(), indices.end(), *i );
00280         assert( *j == *i );
00281         *i = j - indices.begin();
00282     }
00283 
00284     // Allocate space and zero output gradient and Hessian lists
00285     const size_t N     = indices.size();
00286     const size_t m     = mValues.size();
00287     const double inv_n = 1.0 / m;
00288     double v;
00289     gradient.clear();
00290     gradient.resize( N, Vector3D( 0, 0, 0 ) );
00291     Hessian.clear();
00292     Hessian.resize( N * ( N + 1 ) / 2, Matrix3D( 0.0 ) );
00293 
00294     // Average values, gradients, Hessians
00295     Matrix3D outer;
00296     value                                           = 0.0;
00297     std::vector< Matrix3D >::iterator met_hess_iter = mHess.begin();
00298     for( size_t k = 0; k < m; ++k )
00299     {  // for each metric evaluate
00300         if( P.value() == 1.0 )
00301         {
00302             v = P.raise( mValues[k] );
00303             // for each gradient (or Hessian row) for the local metric evaluation
00304             for( size_t r = mOffsets[k]; r < mOffsets[k + 1]; ++r )
00305             {
00306                 const size_t nr = mIndices[r];
00307                 // for each column of the local metric Hessian
00308                 for( size_t c = r; c < mOffsets[k + 1]; ++c )
00309                 {
00310                     const size_t nc = mIndices[c];
00311                     *met_hess_iter *= inv_n;
00312                     if( nr <= nc )
00313                         Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += *met_hess_iter;
00314                     else
00315                         Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( *met_hess_iter );
00316                     ++met_hess_iter;
00317                 }
00318                 mGrad[r] *= inv_n;
00319                 gradient[nr] += mGrad[r];
00320             }
00321         }
00322         else
00323         {
00324             const double r2 = P2.raise( mValues[k] );
00325             const double r1 = r2 * mValues[k];
00326             v               = r1 * mValues[k];
00327             const double h  = r2 * P.value() * ( P.value() - 1 ) * inv_n;
00328             const double g  = r1 * P.value() * inv_n;
00329             // for each gradient (or Hessian row) for the local metric evaluation
00330             for( size_t r = mOffsets[k]; r < mOffsets[k + 1]; ++r )
00331             {
00332                 const size_t nr = mIndices[r];
00333                 // for each column of the local metric Hessian
00334                 for( size_t c = r; c < mOffsets[k + 1]; ++c )
00335                 {
00336                     const size_t nc = mIndices[c];
00337                     *met_hess_iter *= g;
00338                     outer.outer_product( mGrad[r], mGrad[c] );
00339                     outer *= h;
00340                     outer += *met_hess_iter;
00341                     if( nr <= nc )
00342                         Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += outer;
00343                     else
00344                         Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( outer );
00345                     ++met_hess_iter;
00346                 }
00347                 mGrad[r] *= g;
00348                 gradient[nr] += mGrad[r];
00349             }
00350         }
00351         value += v;
00352     }
00353 
00354     value *= inv_n;
00355 
00356     return rval;
00357 }
00358 
00359 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines