MOAB: Mesh Oriented datABase  (version 5.2.1)
MultiplyQualityMetric.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in 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     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
00025 
00026   ***************************************************************** */
00027 /*!
00028   \file   MultiplyQualityMetric.cpp
00029   \brief
00030 
00031   \author Michael Brewer
00032   \date   2002-05-09
00033 */
00034 
00035 #include "MultiplyQualityMetric.hpp"
00036 #include "Vector3D.hpp"
00037 #include "Matrix3D.hpp"
00038 #include "MsqError.hpp"
00039 
00040 using namespace MBMesquite;
00041 
00042 MultiplyQualityMetric::MultiplyQualityMetric( QualityMetric* qm1, QualityMetric* qm2, MsqError& err )
00043     : metric1( *qm1 ), metric2( *qm2 )
00044 {
00045     if( qm1->get_metric_type() != qm2->get_metric_type() || qm1->get_negate_flag() != qm2->get_negate_flag() )
00046     { MSQ_SETERR( err )( "Incompatible metrics", MsqError::INVALID_ARG ); }
00047 }
00048 
00049 MultiplyQualityMetric::~MultiplyQualityMetric() {}
00050 
00051 QualityMetric::MetricType MultiplyQualityMetric::get_metric_type() const
00052 {
00053     return metric1.get_metric_type();
00054 }
00055 
00056 std::string MultiplyQualityMetric::get_name() const
00057 {
00058     return metric1.get_name() + "*" + metric2.get_name();
00059 }
00060 
00061 int MultiplyQualityMetric::get_negate_flag() const
00062 {
00063     return metric1.get_negate_flag();
00064 }
00065 
00066 void MultiplyQualityMetric::get_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free_only,
00067                                              MsqError& err )
00068 {
00069     metric1.get_evaluations( pd, handles, free_only, err );MSQ_ERRRTN( err );
00070     metric2.get_evaluations( pd, mHandles, free_only, err );MSQ_ERRRTN( err );
00071     if( handles != mHandles ) { MSQ_SETERR( err )( "Incompatible metrics", MsqError::INVALID_STATE ); }
00072 }
00073 
00074 bool MultiplyQualityMetric::evaluate( PatchData& pd, size_t handle, double& value, MsqError& err )
00075 {
00076     double val1, val2;
00077     bool rval1, rval2;
00078     rval1 = metric1.evaluate( pd, handle, val1, err );
00079     MSQ_ERRZERO( err );
00080     rval2 = metric2.evaluate( pd, handle, val2, err );
00081     MSQ_ERRZERO( err );
00082     value = val1 * val2;
00083     return rval1 && rval2;
00084 }
00085 
00086 bool MultiplyQualityMetric::evaluate_with_indices( PatchData& pd, size_t handle, double& value,
00087                                                    std::vector< size_t >& indices, MsqError& err )
00088 {
00089     double val1, val2;
00090     bool rval1, rval2;
00091     rval1 = metric1.evaluate_with_indices( pd, handle, val1, indices1, err );
00092     MSQ_ERRZERO( err );
00093     rval2 = metric2.evaluate_with_indices( pd, handle, val2, indices2, err );
00094     MSQ_ERRZERO( err );
00095 
00096     indices.clear();
00097     std::sort( indices1.begin(), indices1.end() );
00098     std::sort( indices2.begin(), indices2.end() );
00099     std::set_union( indices1.begin(), indices1.end(), indices2.begin(), indices2.end(), std::back_inserter( indices ) );
00100 
00101     value = val1 * val2;
00102     return rval1 && rval2;
00103 }
00104 
00105 bool MultiplyQualityMetric::evaluate_with_gradient( PatchData& pd, size_t handle, double& value,
00106                                                     std::vector< size_t >& indices, std::vector< Vector3D >& gradient,
00107                                                     MsqError& err )
00108 {
00109     std::vector< size_t >::iterator i;
00110     size_t j;
00111     double val1, val2;
00112     bool rval1, rval2;
00113     rval1 = metric1.evaluate_with_gradient( pd, handle, val1, indices1, grad1, err );
00114     MSQ_ERRZERO( err );
00115     rval2 = metric2.evaluate_with_gradient( pd, handle, val2, indices2, grad2, err );
00116     MSQ_ERRZERO( err );
00117 
00118     indices.resize( indices1.size() + indices2.size() );
00119     i = std::copy( indices1.begin(), indices1.end(), indices.begin() );
00120     std::copy( indices2.begin(), indices2.end(), i );
00121     std::sort( indices.begin(), indices.end() );
00122     indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
00123 
00124     gradient.clear();
00125     gradient.resize( indices.size(), Vector3D( 0.0 ) );
00126     for( j = 0; j < indices1.size(); ++j )
00127     {
00128         i        = std::lower_bound( indices.begin(), indices.end(), indices1[j] );
00129         size_t k = i - indices.begin();
00130         gradient[k] += val2 * grad1[j];
00131     }
00132     for( j = 0; j < indices2.size(); ++j )
00133     {
00134         i        = std::lower_bound( indices.begin(), indices.end(), indices2[j] );
00135         size_t k = i - indices.begin();
00136         gradient[k] += val1 * grad2[j];
00137     }
00138 
00139     value = val1 * val2;
00140     return rval1 && rval2;
00141 }
00142 
00143 bool MultiplyQualityMetric::evaluate_with_Hessian( PatchData& pd, size_t handle, double& value,
00144                                                    std::vector< size_t >& indices, std::vector< Vector3D >& gradient,
00145                                                    std::vector< Matrix3D >& Hessian, MsqError& err )
00146 {
00147     std::vector< size_t >::iterator i;
00148     size_t j, r, c, n, h;
00149     double val1, val2;
00150     bool rval1, rval2;
00151     rval1 = metric1.evaluate_with_Hessian( pd, handle, val1, indices1, grad1, Hess1, err );
00152     MSQ_ERRZERO( err );
00153     rval2 = metric2.evaluate_with_Hessian( pd, handle, val2, indices2, grad2, Hess2, err );
00154     MSQ_ERRZERO( err );
00155     // merge index lists
00156     indices.resize( indices1.size() + indices2.size() );
00157     i = std::copy( indices1.begin(), indices1.end(), indices.begin() );
00158     std::copy( indices2.begin(), indices2.end(), i );
00159     std::sort( indices.begin(), indices.end() );
00160     indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
00161     // calculate grads and convert index lists to indices into output list
00162     gradient.clear();
00163     gradient.resize( indices.size(), Vector3D( 0.0 ) );
00164     for( j = 0; j < indices1.size(); ++j )
00165     {
00166         i           = std::lower_bound( indices.begin(), indices.end(), indices1[j] );
00167         indices1[j] = i - indices.begin();
00168         gradient[indices1[j]] += val2 * grad1[j];
00169     }
00170     for( j = 0; j < indices2.size(); ++j )
00171     {
00172         i           = std::lower_bound( indices.begin(), indices.end(), indices2[j] );
00173         indices2[j] = i - indices.begin();
00174         gradient[indices2[j]] += val1 * grad2[j];
00175     }
00176     // allocate space for hessians, and zero it
00177     const size_t N = indices.size();
00178     Hessian.clear();
00179     Hessian.resize( N * ( N + 1 ) / 2, Matrix3D( 0.0 ) );
00180     // add hessian terms from first metric
00181     n = indices1.size();
00182     h = 0;
00183     for( r = 0; r < n; ++r )
00184     {
00185         const size_t nr = indices1[r];
00186         for( c = r; c < n; ++c )
00187         {
00188             const size_t nc = indices1[c];
00189             Hess1[h] *= val2;
00190             if( nr <= nc )
00191                 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += Hess1[h];
00192             else
00193                 Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( Hess1[h] );
00194             ++h;
00195         }
00196     }
00197     // add hessian terms from second metric
00198     n = indices2.size();
00199     h = 0;
00200     for( r = 0; r < n; ++r )
00201     {
00202         const size_t nr = indices2[r];
00203         for( c = r; c < n; ++c )
00204         {
00205             const size_t nc = indices2[c];
00206             Hess2[h] *= val1;
00207             if( nr <= nc )
00208                 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += Hess2[h];
00209             else
00210                 Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( Hess2[h] );
00211             ++h;
00212         }
00213     }
00214     // add gradient outer products
00215     n        = indices1.size();
00216     size_t m = indices2.size();
00217     Matrix3D outer;
00218     for( r = 0; r < n; ++r )
00219     {
00220         const size_t nr = indices1[r];
00221         for( c = 0; c < m; ++c )
00222         {
00223             const size_t nc = indices2[c];
00224             outer.outer_product( grad1[r], grad2[c] );
00225             if( nr == nc )
00226                 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += outer.plus_transpose_equal( outer );
00227             else if( nr < nc )
00228                 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += outer;
00229             else
00230                 Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( outer );
00231         }
00232     }
00233 
00234     value = val1 * val2;
00235     return rval1 && rval2;
00236 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines