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