MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 }