MOAB: Mesh Oriented datABase  (version 5.4.0)
AddQualityMetric.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 Lawrence Livermore National Laboratory.  Under
00005     the terms of Contract B545069 with the University of Wisconsin --
00006     Madison, Lawrence Livermore National Laboratory 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     (2006) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file AddQualityMetric.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "AddQualityMetric.hpp"
00034 #include "Vector3D.hpp"
00035 #include "Matrix3D.hpp"
00036 #include "MsqError.hpp"
00037 
00038 namespace MBMesquite
00039 {
00040 
00041 AddQualityMetric::AddQualityMetric( QualityMetric* qm1, QualityMetric* qm2, MsqError& err )
00042     : metric1( *qm1 ), metric2( *qm2 )
00043 {
00044     if( qm1->get_metric_type() != qm2->get_metric_type() || qm1->get_negate_flag() != qm2->get_negate_flag() )
00045     {
00046         MSQ_SETERR( err )( "Incompatible metrics", MsqError::INVALID_ARG );
00047     }
00048 }
00049 
00050 AddQualityMetric::~AddQualityMetric() {}
00051 
00052 QualityMetric::MetricType AddQualityMetric::get_metric_type() const
00053 {
00054     return metric1.get_metric_type();
00055 }
00056 
00057 std::string AddQualityMetric::get_name() const
00058 {
00059     return std::string( "Sum(" ) + metric1.get_name() + ", " + metric2.get_name() + ")";
00060 }
00061 
00062 int AddQualityMetric::get_negate_flag() const
00063 {
00064     return metric1.get_negate_flag();
00065 }
00066 
00067 void AddQualityMetric::get_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free_only, 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 )
00072     {
00073         MSQ_SETERR( err )( "Incompatible metrics", MsqError::INVALID_STATE );
00074     }
00075 }
00076 
00077 bool AddQualityMetric::evaluate( PatchData& pd, size_t handle, double& value, MsqError& err )
00078 {
00079     double val1, val2;
00080     bool rval1, rval2;
00081     rval1 = metric1.evaluate( pd, handle, val1, err );
00082     MSQ_ERRZERO( err );
00083     rval2 = metric2.evaluate( pd, handle, val2, err );
00084     MSQ_ERRZERO( err );
00085     value = val1 + val2;
00086     return rval1 && rval2;
00087 }
00088 
00089 bool AddQualityMetric::evaluate_with_indices( PatchData& pd,
00090                                               size_t handle,
00091                                               double& value,
00092                                               std::vector< size_t >& indices,
00093                                               MsqError& err )
00094 {
00095     double val1, val2;
00096     bool rval1, rval2;
00097     rval1 = metric1.evaluate_with_indices( pd, handle, val1, indices1, err );
00098     MSQ_ERRZERO( err );
00099     rval2 = metric2.evaluate_with_indices( pd, handle, val2, indices2, err );
00100     MSQ_ERRZERO( err );
00101 
00102     indices.clear();
00103     std::sort( indices1.begin(), indices1.end() );
00104     std::sort( indices2.begin(), indices2.end() );
00105     std::set_union( indices1.begin(), indices1.end(), indices2.begin(), indices2.end(), std::back_inserter( indices ) );
00106 
00107     value = val1 + val2;
00108     return rval1 && rval2;
00109 }
00110 
00111 bool AddQualityMetric::evaluate_with_gradient( PatchData& pd,
00112                                                size_t handle,
00113                                                double& value,
00114                                                std::vector< size_t >& indices,
00115                                                std::vector< Vector3D >& gradient,
00116                                                MsqError& err )
00117 {
00118     std::vector< size_t >::iterator i;
00119     size_t j;
00120     double val1, val2;
00121     bool rval1, rval2;
00122     rval1 = metric1.evaluate_with_gradient( pd, handle, val1, indices1, grad1, err );
00123     MSQ_ERRZERO( err );
00124     rval2 = metric2.evaluate_with_gradient( pd, handle, val2, indices2, grad2, err );
00125     MSQ_ERRZERO( err );
00126 
00127     indices.resize( indices1.size() + indices2.size() );
00128     i = std::copy( indices1.begin(), indices1.end(), indices.begin() );
00129     std::copy( indices2.begin(), indices2.end(), i );
00130     std::sort( indices.begin(), indices.end() );
00131     indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
00132 
00133     gradient.clear();
00134     gradient.resize( indices.size(), Vector3D( 0.0 ) );
00135     for( j = 0; j < indices1.size(); ++j )
00136     {
00137         i        = std::lower_bound( indices.begin(), indices.end(), indices1[j] );
00138         size_t k = i - indices.begin();
00139         gradient[k] += grad1[j];
00140     }
00141     for( j = 0; j < indices2.size(); ++j )
00142     {
00143         i        = std::lower_bound( indices.begin(), indices.end(), indices2[j] );
00144         size_t k = i - indices.begin();
00145         gradient[k] += grad2[j];
00146     }
00147 
00148     value = val1 + val2;
00149     return rval1 && rval2;
00150 }
00151 
00152 bool AddQualityMetric::evaluate_with_Hessian( PatchData& pd,
00153                                               size_t handle,
00154                                               double& value,
00155                                               std::vector< size_t >& indices,
00156                                               std::vector< Vector3D >& gradient,
00157                                               std::vector< Matrix3D >& Hessian,
00158                                               MsqError& err )
00159 {
00160     std::vector< size_t >::iterator i;
00161     size_t j, r, c, n, h;
00162     double val1, val2;
00163     bool rval1, rval2;
00164     rval1 = metric1.evaluate_with_Hessian( pd, handle, val1, indices1, grad1, Hess1, err );
00165     MSQ_ERRZERO( err );
00166     rval2 = metric2.evaluate_with_Hessian( pd, handle, val2, indices2, grad2, Hess2, err );
00167     MSQ_ERRZERO( err );
00168 
00169     indices.resize( indices1.size() + indices2.size() );
00170     i = std::copy( indices1.begin(), indices1.end(), indices.begin() );
00171     std::copy( indices2.begin(), indices2.end(), i );
00172     std::sort( indices.begin(), indices.end() );
00173     indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
00174 
00175     gradient.clear();
00176     gradient.resize( indices.size(), Vector3D( 0.0 ) );
00177     for( j = 0; j < indices1.size(); ++j )
00178     {
00179         i           = std::lower_bound( indices.begin(), indices.end(), indices1[j] );
00180         indices1[j] = i - indices.begin();
00181         gradient[indices1[j]] += grad1[j];
00182     }
00183     for( j = 0; j < indices2.size(); ++j )
00184     {
00185         i           = std::lower_bound( indices.begin(), indices.end(), indices2[j] );
00186         indices2[j] = i - indices.begin();
00187         gradient[indices2[j]] += grad2[j];
00188     }
00189 
00190     const size_t N = indices.size();
00191     Hessian.clear();
00192     Hessian.resize( N * ( N + 1 ) / 2, Matrix3D( 0.0 ) );
00193 
00194     n = indices1.size();
00195     h = 0;
00196     for( r = 0; r < n; ++r )
00197     {
00198         const size_t nr = indices1[r];
00199         for( c = r; c < n; ++c )
00200         {
00201             const size_t nc = indices1[c];
00202             if( nr <= nc )
00203                 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += Hess1[h++];
00204             else
00205                 Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( Hess1[h++] );
00206         }
00207     }
00208 
00209     n = indices2.size();
00210     h = 0;
00211     for( r = 0; r < n; ++r )
00212     {
00213         const size_t nr = indices2[r];
00214         for( c = r; c < n; ++c )
00215         {
00216             const size_t nc = indices2[c];
00217             if( nr <= nc )
00218                 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += Hess2[h++];
00219             else
00220                 Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( Hess2[h++] );
00221         }
00222     }
00223 
00224     value = val1 + val2;
00225     return rval1 && rval2;
00226 }
00227 
00228 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines