MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 { MSQ_SETERR( err )( "Incompatible metrics", MsqError::INVALID_ARG ); } 00046 } 00047 00048 AddQualityMetric::~AddQualityMetric() {} 00049 00050 QualityMetric::MetricType AddQualityMetric::get_metric_type() const 00051 { 00052 return metric1.get_metric_type(); 00053 } 00054 00055 std::string AddQualityMetric::get_name() const 00056 { 00057 return std::string( "Sum(" ) + metric1.get_name() + ", " + metric2.get_name() + ")"; 00058 } 00059 00060 int AddQualityMetric::get_negate_flag() const 00061 { 00062 return metric1.get_negate_flag(); 00063 } 00064 00065 void AddQualityMetric::get_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free_only, MsqError& err ) 00066 { 00067 metric1.get_evaluations( pd, handles, free_only, err );MSQ_ERRRTN( err ); 00068 metric2.get_evaluations( pd, mHandles, free_only, err );MSQ_ERRRTN( err ); 00069 if( handles != mHandles ) { MSQ_SETERR( err )( "Incompatible metrics", MsqError::INVALID_STATE ); } 00070 } 00071 00072 bool AddQualityMetric::evaluate( PatchData& pd, size_t handle, double& value, MsqError& err ) 00073 { 00074 double val1, val2; 00075 bool rval1, rval2; 00076 rval1 = metric1.evaluate( pd, handle, val1, err ); 00077 MSQ_ERRZERO( err ); 00078 rval2 = metric2.evaluate( pd, handle, val2, err ); 00079 MSQ_ERRZERO( err ); 00080 value = val1 + val2; 00081 return rval1 && rval2; 00082 } 00083 00084 bool AddQualityMetric::evaluate_with_indices( PatchData& pd, size_t handle, double& value, 00085 std::vector< size_t >& indices, MsqError& err ) 00086 { 00087 double val1, val2; 00088 bool rval1, rval2; 00089 rval1 = metric1.evaluate_with_indices( pd, handle, val1, indices1, err ); 00090 MSQ_ERRZERO( err ); 00091 rval2 = metric2.evaluate_with_indices( pd, handle, val2, indices2, err ); 00092 MSQ_ERRZERO( err ); 00093 00094 indices.clear(); 00095 std::sort( indices1.begin(), indices1.end() ); 00096 std::sort( indices2.begin(), indices2.end() ); 00097 std::set_union( indices1.begin(), indices1.end(), indices2.begin(), indices2.end(), std::back_inserter( indices ) ); 00098 00099 value = val1 + val2; 00100 return rval1 && rval2; 00101 } 00102 00103 bool AddQualityMetric::evaluate_with_gradient( PatchData& pd, size_t handle, double& value, 00104 std::vector< size_t >& indices, std::vector< Vector3D >& gradient, 00105 MsqError& err ) 00106 { 00107 std::vector< size_t >::iterator i; 00108 size_t j; 00109 double val1, val2; 00110 bool rval1, rval2; 00111 rval1 = metric1.evaluate_with_gradient( pd, handle, val1, indices1, grad1, err ); 00112 MSQ_ERRZERO( err ); 00113 rval2 = metric2.evaluate_with_gradient( pd, handle, val2, indices2, grad2, err ); 00114 MSQ_ERRZERO( err ); 00115 00116 indices.resize( indices1.size() + indices2.size() ); 00117 i = std::copy( indices1.begin(), indices1.end(), indices.begin() ); 00118 std::copy( indices2.begin(), indices2.end(), i ); 00119 std::sort( indices.begin(), indices.end() ); 00120 indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() ); 00121 00122 gradient.clear(); 00123 gradient.resize( indices.size(), Vector3D( 0.0 ) ); 00124 for( j = 0; j < indices1.size(); ++j ) 00125 { 00126 i = std::lower_bound( indices.begin(), indices.end(), indices1[j] ); 00127 size_t k = i - indices.begin(); 00128 gradient[k] += grad1[j]; 00129 } 00130 for( j = 0; j < indices2.size(); ++j ) 00131 { 00132 i = std::lower_bound( indices.begin(), indices.end(), indices2[j] ); 00133 size_t k = i - indices.begin(); 00134 gradient[k] += grad2[j]; 00135 } 00136 00137 value = val1 + val2; 00138 return rval1 && rval2; 00139 } 00140 00141 bool AddQualityMetric::evaluate_with_Hessian( PatchData& pd, size_t handle, double& value, 00142 std::vector< size_t >& indices, std::vector< Vector3D >& gradient, 00143 std::vector< Matrix3D >& Hessian, MsqError& err ) 00144 { 00145 std::vector< size_t >::iterator i; 00146 size_t j, r, c, n, h; 00147 double val1, val2; 00148 bool rval1, rval2; 00149 rval1 = metric1.evaluate_with_Hessian( pd, handle, val1, indices1, grad1, Hess1, err ); 00150 MSQ_ERRZERO( err ); 00151 rval2 = metric2.evaluate_with_Hessian( pd, handle, val2, indices2, grad2, Hess2, err ); 00152 MSQ_ERRZERO( err ); 00153 00154 indices.resize( indices1.size() + indices2.size() ); 00155 i = std::copy( indices1.begin(), indices1.end(), indices.begin() ); 00156 std::copy( indices2.begin(), indices2.end(), i ); 00157 std::sort( indices.begin(), indices.end() ); 00158 indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() ); 00159 00160 gradient.clear(); 00161 gradient.resize( indices.size(), Vector3D( 0.0 ) ); 00162 for( j = 0; j < indices1.size(); ++j ) 00163 { 00164 i = std::lower_bound( indices.begin(), indices.end(), indices1[j] ); 00165 indices1[j] = i - indices.begin(); 00166 gradient[indices1[j]] += grad1[j]; 00167 } 00168 for( j = 0; j < indices2.size(); ++j ) 00169 { 00170 i = std::lower_bound( indices.begin(), indices.end(), indices2[j] ); 00171 indices2[j] = i - indices.begin(); 00172 gradient[indices2[j]] += grad2[j]; 00173 } 00174 00175 const size_t N = indices.size(); 00176 Hessian.clear(); 00177 Hessian.resize( N * ( N + 1 ) / 2, Matrix3D( 0.0 ) ); 00178 00179 n = indices1.size(); 00180 h = 0; 00181 for( r = 0; r < n; ++r ) 00182 { 00183 const size_t nr = indices1[r]; 00184 for( c = r; c < n; ++c ) 00185 { 00186 const size_t nc = indices1[c]; 00187 if( nr <= nc ) 00188 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += Hess1[h++]; 00189 else 00190 Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( Hess1[h++] ); 00191 } 00192 } 00193 00194 n = indices2.size(); 00195 h = 0; 00196 for( r = 0; r < n; ++r ) 00197 { 00198 const size_t nr = indices2[r]; 00199 for( c = r; c < n; ++c ) 00200 { 00201 const size_t nc = indices2[c]; 00202 if( nr <= nc ) 00203 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += Hess2[h++]; 00204 else 00205 Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( Hess2[h++] ); 00206 } 00207 } 00208 00209 value = val1 + val2; 00210 return rval1 && rval2; 00211 } 00212 00213 } // namespace MBMesquite