MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2006 Sandia National Laboratories. Developed at the 00005 University of Wisconsin--Madison under SNL contract number 00006 624796. The U.S. Government and the University of Wisconsin 00007 retain certain rights to 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 PMeanPMetric.cpp 00028 * \brief 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #include "Mesquite.hpp" 00033 #include "PMeanPMetric.hpp" 00034 #include "MsqError.hpp" 00035 #include "QualityMetric.hpp" 00036 #include "Vector3D.hpp" 00037 #include "Matrix3D.hpp" 00038 #include "SymMatrix3D.hpp" 00039 #include "PatchData.hpp" 00040 00041 namespace MBMesquite 00042 { 00043 00044 bool PMeanPMetric::average( PatchData& pd, QualityMetric* metric, const std::vector< size_t >& qm_handles, 00045 double& value, MsqError& err ) 00046 { 00047 bool rval = true; 00048 value = 0; 00049 for( std::vector< size_t >::const_iterator i = qm_handles.begin(); i != qm_handles.end(); ++i ) 00050 { 00051 double mval; 00052 if( !metric->evaluate( pd, *i, mval, err ) ) rval = false; 00053 MSQ_ERRZERO( err ); 00054 value += P.raise( mval ); 00055 } 00056 value /= qm_handles.size(); 00057 return rval; 00058 } 00059 00060 bool PMeanPMetric::average_with_indices( PatchData& pd, QualityMetric* metric, const std::vector< size_t >& qm_handles, 00061 double& value, std::vector< size_t >& indices, MsqError& err ) 00062 { 00063 indices.clear(); 00064 00065 bool rval = true; 00066 value = 0; 00067 for( std::vector< size_t >::const_iterator i = qm_handles.begin(); i != qm_handles.end(); ++i ) 00068 { 00069 double mval; 00070 mIndices.clear(); 00071 if( !metric->evaluate_with_indices( pd, *i, mval, mIndices, err ) ) rval = false; 00072 MSQ_ERRZERO( err ); 00073 value += P.raise( mval ); 00074 00075 std::copy( mIndices.begin(), mIndices.end(), std::back_inserter( indices ) ); 00076 } 00077 std::sort( indices.begin(), indices.end() ); 00078 indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() ); 00079 00080 value /= qm_handles.size(); 00081 return rval; 00082 } 00083 00084 bool PMeanPMetric::average_with_gradient( PatchData& pd, QualityMetric* metric, const std::vector< size_t >& qm_handles, 00085 double& value, std::vector< size_t >& indices, 00086 std::vector< Vector3D >& gradient, MsqError& err ) 00087 { 00088 indices.clear(); 00089 gradient.clear(); 00090 00091 std::vector< Vector3D >::iterator g; 00092 std::vector< size_t >::iterator j, k; 00093 std::vector< size_t >::const_iterator i; 00094 00095 bool rval = true; 00096 value = 0; 00097 for( i = qm_handles.begin(); i != qm_handles.end(); ++i ) 00098 { 00099 00100 double mval; 00101 mIndices.clear(); 00102 mGrad.clear(); 00103 if( !metric->evaluate_with_gradient( pd, *i, mval, mIndices, mGrad, err ) ) rval = false; 00104 MSQ_ERRZERO( err ); 00105 value += P.raise( mval ); 00106 00107 double p1val = P.value() * P1.raise( mval ); 00108 for( j = mIndices.begin(), g = mGrad.begin(); j != mIndices.end(); ++j, ++g ) 00109 { 00110 00111 *g *= p1val; 00112 k = std::lower_bound( indices.begin(), indices.end(), *j ); 00113 if( k == indices.end() || *k != *j ) 00114 { 00115 k = indices.insert( k, *j ); 00116 size_t idx = k - indices.begin(); 00117 gradient.insert( gradient.begin() + idx, *g ); 00118 } 00119 else 00120 { 00121 size_t idx = k - indices.begin(); 00122 gradient[idx] += *g; 00123 } 00124 } 00125 } 00126 00127 double inv_n = 1.0 / qm_handles.size(); 00128 value *= inv_n; 00129 for( g = gradient.begin(); g != gradient.end(); ++g ) 00130 *g *= inv_n; 00131 00132 return rval; 00133 } 00134 00135 bool PMeanPMetric::average_with_Hessian_diagonal( PatchData& pd, QualityMetric* metric, 00136 const std::vector< size_t >& qm_handles, double& value, 00137 std::vector< size_t >& indices, std::vector< Vector3D >& gradient, 00138 std::vector< SymMatrix3D >& diagonal, MsqError& err ) 00139 { 00140 // clear temporary storage 00141 mIndices.clear(); 00142 mGrad.clear(); 00143 mDiag.clear(); 00144 mOffsets.clear(); 00145 mValues.clear(); 00146 00147 // Evaluate metric for all sample points, 00148 // accumulating indices, gradients, and Hessians 00149 bool rval = true; 00150 std::vector< size_t >::const_iterator q; 00151 for( q = qm_handles.begin(); q != qm_handles.end(); ++q ) 00152 { 00153 double mval; 00154 indices.clear(); 00155 gradient.clear(); 00156 diagonal.clear(); 00157 if( !metric->evaluate_with_Hessian_diagonal( pd, *q, mval, indices, gradient, diagonal, err ) ) rval = false; 00158 MSQ_ERRZERO( err ); 00159 00160 mValues.push_back( mval ); 00161 mOffsets.push_back( mIndices.size() ); 00162 std::copy( indices.begin(), indices.end(), std::back_inserter( mIndices ) ); 00163 std::copy( gradient.begin(), gradient.end(), std::back_inserter( mGrad ) ); 00164 std::copy( diagonal.begin(), diagonal.end(), std::back_inserter( mDiag ) ); 00165 } 00166 mOffsets.push_back( mIndices.size() ); 00167 00168 // Combine lists of free vertex indices, and update indices 00169 // in per-evaluation lists to point into the combined gradient 00170 // and Hessian lists. 00171 indices = mIndices; 00172 std::sort( indices.begin(), indices.end() ); 00173 indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() ); 00174 std::vector< size_t >::iterator i, j; 00175 for( i = mIndices.begin(); i != mIndices.end(); ++i ) 00176 { 00177 j = std::lower_bound( indices.begin(), indices.end(), *i ); 00178 assert( *j == *i ); 00179 *i = j - indices.begin(); 00180 } 00181 00182 // Allocate space and zero output gradient and Hessian lists 00183 const size_t n = indices.size(); 00184 const size_t m = mValues.size(); 00185 const double inv_n = 1.0 / m; 00186 double v; 00187 gradient.clear(); 00188 gradient.resize( n, Vector3D( 0, 0, 0 ) ); 00189 diagonal.clear(); 00190 diagonal.resize( n, SymMatrix3D( 0.0 ) ); 00191 00192 // Average values, gradients, Hessians 00193 value = 0.0; 00194 for( size_t k = 0; k < m; ++k ) 00195 { // for each metric evaluate 00196 if( P.value() == 1.0 ) 00197 { 00198 v = P.raise( mValues[k] ); 00199 // for each gradient (or Hessian) for the local metric evaluation 00200 for( size_t r = mOffsets[k]; r < mOffsets[k + 1]; ++r ) 00201 { 00202 const size_t nr = mIndices[r]; 00203 00204 mDiag[r] *= inv_n; 00205 diagonal[nr] += mDiag[r]; 00206 mGrad[r] *= inv_n; 00207 gradient[nr] += mGrad[r]; 00208 } 00209 } 00210 else 00211 { 00212 const double r2 = P2.raise( mValues[k] ); 00213 const double r1 = r2 * mValues[k]; 00214 v = r1 * mValues[k]; 00215 const double h = r2 * P.value() * ( P.value() - 1 ) * inv_n; 00216 const double g = r1 * P.value() * inv_n; 00217 // for each gradient (or Hessian) for the local metric evaluation 00218 for( size_t r = mOffsets[k]; r < mOffsets[k + 1]; ++r ) 00219 { 00220 const size_t nr = mIndices[r]; 00221 00222 mDiag[r] *= g; 00223 diagonal[nr] += mDiag[r]; 00224 diagonal[nr] += h * outer( mGrad[r] ); 00225 mGrad[r] *= g; 00226 gradient[nr] += mGrad[r]; 00227 } 00228 } 00229 value += v; 00230 } 00231 00232 value *= inv_n; 00233 00234 return rval; 00235 } 00236 00237 bool PMeanPMetric::average_with_Hessian( PatchData& pd, QualityMetric* metric, const std::vector< size_t >& qm_handles, 00238 double& value, std::vector< size_t >& indices, 00239 std::vector< Vector3D >& gradient, std::vector< Matrix3D >& Hessian, 00240 MsqError& err ) 00241 { 00242 // clear temporary storage 00243 mIndices.clear(); 00244 mGrad.clear(); 00245 mHess.clear(); 00246 mOffsets.clear(); 00247 mValues.clear(); 00248 00249 // Evaluate metric for all sample points, 00250 // accumulating indices, gradients, and Hessians 00251 bool rval = true; 00252 std::vector< size_t >::const_iterator q; 00253 for( q = qm_handles.begin(); q != qm_handles.end(); ++q ) 00254 { 00255 double mval; 00256 indices.clear(); 00257 gradient.clear(); 00258 Hessian.clear(); 00259 if( !metric->evaluate_with_Hessian( pd, *q, mval, indices, gradient, Hessian, err ) ) rval = false; 00260 MSQ_ERRZERO( err ); 00261 00262 mValues.push_back( mval ); 00263 mOffsets.push_back( mIndices.size() ); 00264 std::copy( indices.begin(), indices.end(), std::back_inserter( mIndices ) ); 00265 std::copy( gradient.begin(), gradient.end(), std::back_inserter( mGrad ) ); 00266 std::copy( Hessian.begin(), Hessian.end(), std::back_inserter( mHess ) ); 00267 } 00268 mOffsets.push_back( mIndices.size() ); 00269 00270 // Combine lists of free vertex indices, and update indices 00271 // in per-evaluation lists to point into the combined gradient 00272 // and Hessian lists. 00273 indices = mIndices; 00274 std::sort( indices.begin(), indices.end() ); 00275 indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() ); 00276 std::vector< size_t >::iterator i, j; 00277 for( i = mIndices.begin(); i != mIndices.end(); ++i ) 00278 { 00279 j = std::lower_bound( indices.begin(), indices.end(), *i ); 00280 assert( *j == *i ); 00281 *i = j - indices.begin(); 00282 } 00283 00284 // Allocate space and zero output gradient and Hessian lists 00285 const size_t N = indices.size(); 00286 const size_t m = mValues.size(); 00287 const double inv_n = 1.0 / m; 00288 double v; 00289 gradient.clear(); 00290 gradient.resize( N, Vector3D( 0, 0, 0 ) ); 00291 Hessian.clear(); 00292 Hessian.resize( N * ( N + 1 ) / 2, Matrix3D( 0.0 ) ); 00293 00294 // Average values, gradients, Hessians 00295 Matrix3D outer; 00296 value = 0.0; 00297 std::vector< Matrix3D >::iterator met_hess_iter = mHess.begin(); 00298 for( size_t k = 0; k < m; ++k ) 00299 { // for each metric evaluate 00300 if( P.value() == 1.0 ) 00301 { 00302 v = P.raise( mValues[k] ); 00303 // for each gradient (or Hessian row) for the local metric evaluation 00304 for( size_t r = mOffsets[k]; r < mOffsets[k + 1]; ++r ) 00305 { 00306 const size_t nr = mIndices[r]; 00307 // for each column of the local metric Hessian 00308 for( size_t c = r; c < mOffsets[k + 1]; ++c ) 00309 { 00310 const size_t nc = mIndices[c]; 00311 *met_hess_iter *= inv_n; 00312 if( nr <= nc ) 00313 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += *met_hess_iter; 00314 else 00315 Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( *met_hess_iter ); 00316 ++met_hess_iter; 00317 } 00318 mGrad[r] *= inv_n; 00319 gradient[nr] += mGrad[r]; 00320 } 00321 } 00322 else 00323 { 00324 const double r2 = P2.raise( mValues[k] ); 00325 const double r1 = r2 * mValues[k]; 00326 v = r1 * mValues[k]; 00327 const double h = r2 * P.value() * ( P.value() - 1 ) * inv_n; 00328 const double g = r1 * P.value() * inv_n; 00329 // for each gradient (or Hessian row) for the local metric evaluation 00330 for( size_t r = mOffsets[k]; r < mOffsets[k + 1]; ++r ) 00331 { 00332 const size_t nr = mIndices[r]; 00333 // for each column of the local metric Hessian 00334 for( size_t c = r; c < mOffsets[k + 1]; ++c ) 00335 { 00336 const size_t nc = mIndices[c]; 00337 *met_hess_iter *= g; 00338 outer.outer_product( mGrad[r], mGrad[c] ); 00339 outer *= h; 00340 outer += *met_hess_iter; 00341 if( nr <= nc ) 00342 Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += outer; 00343 else 00344 Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( outer ); 00345 ++met_hess_iter; 00346 } 00347 mGrad[r] *= g; 00348 gradient[nr] += mGrad[r]; 00349 } 00350 } 00351 value += v; 00352 } 00353 00354 value *= inv_n; 00355 00356 return rval; 00357 } 00358 00359 } // namespace MBMesquite