MOAB: Mesh Oriented datABase  (version 5.4.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     {
00047         MSQ_SETERR( err )( "Incompatible metrics", MsqError::INVALID_ARG );
00048     }
00049 }
00050 
00051 MultiplyQualityMetric::~MultiplyQualityMetric() {}
00052 
00053 QualityMetric::MetricType MultiplyQualityMetric::get_metric_type() const
00054 {
00055     return metric1.get_metric_type();
00056 }
00057 
00058 std::string MultiplyQualityMetric::get_name() const
00059 {
00060     return metric1.get_name() + "*" + metric2.get_name();
00061 }
00062 
00063 int MultiplyQualityMetric::get_negate_flag() const
00064 {
00065     return metric1.get_negate_flag();
00066 }
00067 
00068 void MultiplyQualityMetric::get_evaluations( PatchData& pd,
00069                                              std::vector< size_t >& handles,
00070                                              bool free_only,
00071                                              MsqError& err )
00072 {
00073     metric1.get_evaluations( pd, handles, free_only, err );MSQ_ERRRTN( err );
00074     metric2.get_evaluations( pd, mHandles, free_only, err );MSQ_ERRRTN( err );
00075     if( handles != mHandles )
00076     {
00077         MSQ_SETERR( err )( "Incompatible metrics", MsqError::INVALID_STATE );
00078     }
00079 }
00080 
00081 bool MultiplyQualityMetric::evaluate( PatchData& pd, size_t handle, double& value, MsqError& err )
00082 {
00083     double val1, val2;
00084     bool rval1, rval2;
00085     rval1 = metric1.evaluate( pd, handle, val1, err );
00086     MSQ_ERRZERO( err );
00087     rval2 = metric2.evaluate( pd, handle, val2, err );
00088     MSQ_ERRZERO( err );
00089     value = val1 * val2;
00090     return rval1 && rval2;
00091 }
00092 
00093 bool MultiplyQualityMetric::evaluate_with_indices( PatchData& pd,
00094                                                    size_t handle,
00095                                                    double& value,
00096                                                    std::vector< size_t >& indices,
00097                                                    MsqError& err )
00098 {
00099     double val1, val2;
00100     bool rval1, rval2;
00101     rval1 = metric1.evaluate_with_indices( pd, handle, val1, indices1, err );
00102     MSQ_ERRZERO( err );
00103     rval2 = metric2.evaluate_with_indices( pd, handle, val2, indices2, err );
00104     MSQ_ERRZERO( err );
00105 
00106     indices.clear();
00107     std::sort( indices1.begin(), indices1.end() );
00108     std::sort( indices2.begin(), indices2.end() );
00109     std::set_union( indices1.begin(), indices1.end(), indices2.begin(), indices2.end(), std::back_inserter( indices ) );
00110 
00111     value = val1 * val2;
00112     return rval1 && rval2;
00113 }
00114 
00115 bool MultiplyQualityMetric::evaluate_with_gradient( PatchData& pd,
00116                                                     size_t handle,
00117                                                     double& value,
00118                                                     std::vector< size_t >& indices,
00119                                                     std::vector< Vector3D >& gradient,
00120                                                     MsqError& err )
00121 {
00122     std::vector< size_t >::iterator i;
00123     size_t j;
00124     double val1, val2;
00125     bool rval1, rval2;
00126     rval1 = metric1.evaluate_with_gradient( pd, handle, val1, indices1, grad1, err );
00127     MSQ_ERRZERO( err );
00128     rval2 = metric2.evaluate_with_gradient( pd, handle, val2, indices2, grad2, err );
00129     MSQ_ERRZERO( err );
00130 
00131     indices.resize( indices1.size() + indices2.size() );
00132     i = std::copy( indices1.begin(), indices1.end(), indices.begin() );
00133     std::copy( indices2.begin(), indices2.end(), i );
00134     std::sort( indices.begin(), indices.end() );
00135     indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
00136 
00137     gradient.clear();
00138     gradient.resize( indices.size(), Vector3D( 0.0 ) );
00139     for( j = 0; j < indices1.size(); ++j )
00140     {
00141         i        = std::lower_bound( indices.begin(), indices.end(), indices1[j] );
00142         size_t k = i - indices.begin();
00143         gradient[k] += val2 * grad1[j];
00144     }
00145     for( j = 0; j < indices2.size(); ++j )
00146     {
00147         i        = std::lower_bound( indices.begin(), indices.end(), indices2[j] );
00148         size_t k = i - indices.begin();
00149         gradient[k] += val1 * grad2[j];
00150     }
00151 
00152     value = val1 * val2;
00153     return rval1 && rval2;
00154 }
00155 
00156 bool MultiplyQualityMetric::evaluate_with_Hessian( PatchData& pd,
00157                                                    size_t handle,
00158                                                    double& value,
00159                                                    std::vector< size_t >& indices,
00160                                                    std::vector< Vector3D >& gradient,
00161                                                    std::vector< Matrix3D >& Hessian,
00162                                                    MsqError& err )
00163 {
00164     std::vector< size_t >::iterator i;
00165     size_t j, r, c, n, h;
00166     double val1, val2;
00167     bool rval1, rval2;
00168     rval1 = metric1.evaluate_with_Hessian( pd, handle, val1, indices1, grad1, Hess1, err );
00169     MSQ_ERRZERO( err );
00170     rval2 = metric2.evaluate_with_Hessian( pd, handle, val2, indices2, grad2, Hess2, err );
00171     MSQ_ERRZERO( err );
00172     // merge index lists
00173     indices.resize( indices1.size() + indices2.size() );
00174     i = std::copy( indices1.begin(), indices1.end(), indices.begin() );
00175     std::copy( indices2.begin(), indices2.end(), i );
00176     std::sort( indices.begin(), indices.end() );
00177     indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
00178     // calculate grads and convert index lists to indices into output list
00179     gradient.clear();
00180     gradient.resize( indices.size(), Vector3D( 0.0 ) );
00181     for( j = 0; j < indices1.size(); ++j )
00182     {
00183         i           = std::lower_bound( indices.begin(), indices.end(), indices1[j] );
00184         indices1[j] = i - indices.begin();
00185         gradient[indices1[j]] += val2 * grad1[j];
00186     }
00187     for( j = 0; j < indices2.size(); ++j )
00188     {
00189         i           = std::lower_bound( indices.begin(), indices.end(), indices2[j] );
00190         indices2[j] = i - indices.begin();
00191         gradient[indices2[j]] += val1 * grad2[j];
00192     }
00193     // allocate space for hessians, and zero it
00194     const size_t N = indices.size();
00195     Hessian.clear();
00196     Hessian.resize( N * ( N + 1 ) / 2, Matrix3D( 0.0 ) );
00197     // add hessian terms from first metric
00198     n = indices1.size();
00199     h = 0;
00200     for( r = 0; r < n; ++r )
00201     {
00202         const size_t nr = indices1[r];
00203         for( c = r; c < n; ++c )
00204         {
00205             const size_t nc = indices1[c];
00206             Hess1[h] *= val2;
00207             if( nr <= nc )
00208                 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += Hess1[h];
00209             else
00210                 Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( Hess1[h] );
00211             ++h;
00212         }
00213     }
00214     // add hessian terms from second metric
00215     n = indices2.size();
00216     h = 0;
00217     for( r = 0; r < n; ++r )
00218     {
00219         const size_t nr = indices2[r];
00220         for( c = r; c < n; ++c )
00221         {
00222             const size_t nc = indices2[c];
00223             Hess2[h] *= val1;
00224             if( nr <= nc )
00225                 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += Hess2[h];
00226             else
00227                 Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( Hess2[h] );
00228             ++h;
00229         }
00230     }
00231     // add gradient outer products
00232     n        = indices1.size();
00233     size_t m = indices2.size();
00234     Matrix3D outer;
00235     for( r = 0; r < n; ++r )
00236     {
00237         const size_t nr = indices1[r];
00238         for( c = 0; c < m; ++c )
00239         {
00240             const size_t nc = indices2[c];
00241             outer.outer_product( grad1[r], grad2[c] );
00242             if( nr == nc )
00243                 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += outer.plus_transpose_equal( outer );
00244             else if( nr < nc )
00245                 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += outer;
00246             else
00247                 Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( outer );
00248         }
00249     }
00250 
00251     value = val1 * val2;
00252     return rval1 && rval2;
00253 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines