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