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 TQualityMetric.cpp 00028 * \brief 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #undef PRINT_INFO 00033 00034 #include "Mesquite.hpp" 00035 #include "TQualityMetric.hpp" 00036 #include "MsqMatrix.hpp" 00037 #include "ElementQM.hpp" 00038 #include "MsqError.hpp" 00039 #include "Vector3D.hpp" 00040 #include "PatchData.hpp" 00041 #include "MappingFunction.hpp" 00042 #include "WeightCalculator.hpp" 00043 #include "TargetCalculator.hpp" 00044 #include "TMetric.hpp" 00045 #include "TMetricBarrier.hpp" 00046 #include "TargetMetricUtil.hpp" 00047 #include "TMPDerivs.hpp" 00048 00049 #ifdef PRINT_INFO 00050 #include <iostream> 00051 #endif 00052 00053 #include <functional> 00054 #include <algorithm> 00055 00056 #define NUMERICAL_2D_HESSIAN 00057 00058 namespace MBMesquite 00059 { 00060 00061 std::string TQualityMetric::get_name() const 00062 { 00063 return targetMetric->get_name(); 00064 } 00065 00066 bool TQualityMetric::evaluate_internal( PatchData& pd, size_t p_handle, double& value, size_t* indices, 00067 size_t& num_indices, MsqError& err ) 00068 { 00069 const Sample s = ElemSampleQM::sample( p_handle ); 00070 const size_t e = ElemSampleQM::elem( p_handle ); 00071 MsqMeshEntity& p_elem = pd.element_by_index( e ); 00072 EntityTopology type = p_elem.get_element_type(); 00073 unsigned edim = TopologyInfo::dimension( type ); 00074 const NodeSet bits = pd.non_slave_node_set( e ); 00075 00076 bool rval; 00077 if( edim == 3 ) 00078 { // 3x3 or 3x2 targets ? 00079 const MappingFunction3D* mf = pd.get_mapping_function_3D( type ); 00080 if( !mf ) 00081 { 00082 MSQ_SETERR( err ) 00083 ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT ); 00084 return false; 00085 } 00086 00087 MsqMatrix< 3, 3 > A, W; 00088 mf->jacobian( pd, e, bits, s, indices, mDerivs3D, num_indices, A, err ); 00089 MSQ_ERRZERO( err ); 00090 targetCalc->get_3D_target( pd, e, s, W, err ); 00091 MSQ_ERRZERO( err ); 00092 const MsqMatrix< 3, 3 > Winv = inverse( W ); 00093 const MsqMatrix< 3, 3 > T = A * Winv; 00094 rval = targetMetric->evaluate( T, value, err ); 00095 MSQ_ERRZERO( err ); 00096 #ifdef PRINT_INFO 00097 print_info< 3 >( e, s, A, W, A * inverse( W ) ); 00098 #endif 00099 } 00100 else if( edim == 2 ) 00101 { 00102 MsqMatrix< 2, 2 > W, A; 00103 MsqMatrix< 3, 2 > S_a_transpose_Theta; 00104 rval = 00105 evaluate_surface_common( pd, s, e, bits, indices, num_indices, mDerivs2D, W, A, S_a_transpose_Theta, err ); 00106 if( MSQ_CHKERR( err ) || !rval ) return false; 00107 const MsqMatrix< 2, 2 > Winv = inverse( W ); 00108 const MsqMatrix< 2, 2 > T = A * Winv; 00109 rval = targetMetric->evaluate( T, value, err ); 00110 MSQ_ERRZERO( err ); 00111 00112 #ifdef PRINT_INFO 00113 print_info< 2 >( e, s, J, Wp, A * inverse( W ) ); 00114 #endif 00115 } 00116 else 00117 { 00118 assert( false ); 00119 return false; 00120 } 00121 00122 return rval; 00123 } 00124 00125 bool TQualityMetric::evaluate_with_gradient( PatchData& pd, size_t p_handle, double& value, 00126 std::vector< size_t >& indices, std::vector< Vector3D >& grad, 00127 MsqError& err ) 00128 { 00129 const Sample s = ElemSampleQM::sample( p_handle ); 00130 const size_t e = ElemSampleQM::elem( p_handle ); 00131 MsqMeshEntity& p_elem = pd.element_by_index( e ); 00132 EntityTopology type = p_elem.get_element_type(); 00133 unsigned edim = TopologyInfo::dimension( type ); 00134 size_t num_idx = 0; 00135 const NodeSet bits = pd.non_slave_node_set( e ); 00136 00137 bool rval; 00138 if( edim == 3 ) 00139 { // 3x3 or 3x2 targets ? 00140 const MappingFunction3D* mf = pd.get_mapping_function_3D( type ); 00141 if( !mf ) 00142 { 00143 MSQ_SETERR( err ) 00144 ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT ); 00145 return false; 00146 } 00147 00148 MsqMatrix< 3, 3 > A, W, dmdT; 00149 mf->jacobian( pd, e, bits, s, mIndices, mDerivs3D, num_idx, A, err ); 00150 MSQ_ERRZERO( err ); 00151 targetCalc->get_3D_target( pd, e, s, W, err ); 00152 MSQ_ERRZERO( err ); 00153 const MsqMatrix< 3, 3 > Winv = inverse( W ); 00154 const MsqMatrix< 3, 3 > T = A * Winv; 00155 rval = targetMetric->evaluate_with_grad( T, value, dmdT, err ); 00156 MSQ_ERRZERO( err ); 00157 gradient< 3 >( num_idx, mDerivs3D, dmdT * transpose( Winv ), grad ); 00158 #ifdef PRINT_INFO 00159 print_info< 3 >( e, s, A, W, A * inverse( W ) ); 00160 #endif 00161 } 00162 else if( edim == 2 ) 00163 { 00164 MsqMatrix< 2, 2 > W, A, dmdT; 00165 MsqMatrix< 3, 2 > S_a_transpose_Theta; 00166 rval = evaluate_surface_common( pd, s, e, bits, mIndices, num_idx, mDerivs2D, W, A, S_a_transpose_Theta, err ); 00167 if( MSQ_CHKERR( err ) || !rval ) return false; 00168 const MsqMatrix< 2, 2 > Winv = inverse( W ); 00169 const MsqMatrix< 2, 2 > T = A * Winv; 00170 rval = targetMetric->evaluate_with_grad( T, value, dmdT, err ); 00171 MSQ_ERRZERO( err ); 00172 gradient< 2 >( num_idx, mDerivs2D, S_a_transpose_Theta * dmdT * transpose( Winv ), grad ); 00173 #ifdef PRINT_INFO 00174 print_info< 2 >( e, s, J, Wp, A * inverse( W ) ); 00175 #endif 00176 } 00177 else 00178 { 00179 assert( false ); 00180 return false; 00181 } 00182 00183 // pass back index list 00184 indices.resize( num_idx ); 00185 std::copy( mIndices, mIndices + num_idx, indices.begin() ); 00186 00187 // apply target weight to value 00188 weight( pd, s, e, num_idx, value, grad.empty() ? 0 : arrptr( grad ), 0, 0, err ); 00189 MSQ_ERRZERO( err ); 00190 return rval; 00191 } 00192 00193 bool TQualityMetric::evaluate_with_Hessian( PatchData& pd, size_t p_handle, double& value, 00194 std::vector< size_t >& indices, std::vector< Vector3D >& grad, 00195 std::vector< Matrix3D >& Hessian, MsqError& err ) 00196 { 00197 const Sample s = ElemSampleQM::sample( p_handle ); 00198 const size_t e = ElemSampleQM::elem( p_handle ); 00199 MsqMeshEntity& p_elem = pd.element_by_index( e ); 00200 EntityTopology type = p_elem.get_element_type(); 00201 unsigned edim = TopologyInfo::dimension( type ); 00202 size_t num_idx = 0; 00203 const NodeSet bits = pd.non_slave_node_set( e ); 00204 00205 bool rval; 00206 if( edim == 3 ) 00207 { // 3x3 or 3x2 targets ? 00208 const MappingFunction3D* mf = pd.get_mapping_function_3D( type ); 00209 if( !mf ) 00210 { 00211 MSQ_SETERR( err ) 00212 ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT ); 00213 return false; 00214 } 00215 00216 MsqMatrix< 3, 3 > A, W, dmdT, d2mdT2[6]; 00217 mf->jacobian( pd, e, bits, s, mIndices, mDerivs3D, num_idx, A, err ); 00218 MSQ_ERRZERO( err ); 00219 targetCalc->get_3D_target( pd, e, s, W, err ); 00220 MSQ_ERRZERO( err ); 00221 const MsqMatrix< 3, 3 > Winv = inverse( W ); 00222 const MsqMatrix< 3, 3 > T = A * Winv; 00223 rval = targetMetric->evaluate_with_hess( T, value, dmdT, d2mdT2, err ); 00224 MSQ_ERRZERO( err ); 00225 gradient< 3 >( num_idx, mDerivs3D, dmdT * transpose( Winv ), grad ); 00226 second_deriv_wrt_product_factor( d2mdT2, Winv ); 00227 Hessian.resize( num_idx * ( num_idx + 1 ) / 2 ); 00228 if( num_idx ) hessian< 3 >( num_idx, mDerivs3D, d2mdT2, arrptr( Hessian ) ); 00229 00230 #ifdef PRINT_INFO 00231 print_info< 3 >( e, s, A, W, A * inverse( W ) ); 00232 #endif 00233 } 00234 else if( edim == 2 ) 00235 { 00236 #ifdef NUMERICAL_2D_HESSIAN 00237 // return finite difference approximation for now 00238 00239 return QualityMetric::evaluate_with_Hessian( pd, p_handle, value, indices, grad, Hessian, err ); 00240 #else 00241 MsqMatrix< 2, 2 > W, A, dmdT, d2mdT2[3]; 00242 MsqMatrix< 3, 2 > M; 00243 rval = evaluate_surface_common( pd, s, e, bits, mIndices, num_idx, mDerivs2D, W, A, M, err ); 00244 if( MSQ_CHKERR( err ) || !rval ) return false; 00245 const MsqMatrix< 2, 2 > Winv = inverse( W ); 00246 const MsqMatrix< 2, 2 > T = A * Winv; 00247 rval = targetMetric->evaluate_with_hess( T, value, dmdT, d2mdT2, err ); 00248 MSQ_ERRZERO( err ); 00249 gradient< 2 >( num_idx, mDerivs2D, M * dmdT * transpose( Winv ), grad ); 00250 // calculate 2D hessian 00251 second_deriv_wrt_product_factor( d2mdT2, Winv ); 00252 const size_t n = num_idx * ( num_idx + 1 ) / 2; 00253 hess2d.resize( n ); 00254 if( n ) hessian< 2 >( num_idx, mDerivs2D, d2mdT2, arrptr( hess2d ) ); 00255 // calculate surface hessian as transform of 2D hessian 00256 Hessian.resize( n ); 00257 for( size_t i = 0; i < n; ++i ) 00258 Hessian[i] = Matrix3D( ( M * hess2d[i] * transpose( M ) ).data() ); 00259 #ifdef PRINT_INFO 00260 print_info< 2 >( e, s, J, Wp, A * inverse( W ) ); 00261 #endif 00262 #endif 00263 } 00264 else 00265 { 00266 assert( 0 ); 00267 return false; 00268 } 00269 00270 // pass back index list 00271 indices.resize( num_idx ); 00272 std::copy( mIndices, mIndices + num_idx, indices.begin() ); 00273 00274 // apply target weight to value 00275 if( !num_idx ) 00276 weight( pd, s, e, num_idx, value, 0, 0, 0, err ); 00277 else 00278 weight( pd, s, e, num_idx, value, arrptr( grad ), 0, arrptr( Hessian ), err ); 00279 MSQ_ERRZERO( err ); 00280 return rval; 00281 } 00282 00283 bool TQualityMetric::evaluate_with_Hessian_diagonal( PatchData& pd, size_t p_handle, double& value, 00284 std::vector< size_t >& indices, std::vector< Vector3D >& grad, 00285 std::vector< SymMatrix3D >& diagonal, MsqError& err ) 00286 { 00287 const Sample s = ElemSampleQM::sample( p_handle ); 00288 const size_t e = ElemSampleQM::elem( p_handle ); 00289 MsqMeshEntity& p_elem = pd.element_by_index( e ); 00290 EntityTopology type = p_elem.get_element_type(); 00291 unsigned edim = TopologyInfo::dimension( type ); 00292 size_t num_idx = 0; 00293 const NodeSet bits = pd.non_slave_node_set( e ); 00294 00295 bool rval; 00296 if( edim == 3 ) 00297 { // 3x3 or 3x2 targets ? 00298 const MappingFunction3D* mf = pd.get_mapping_function_3D( type ); 00299 if( !mf ) 00300 { 00301 MSQ_SETERR( err ) 00302 ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT ); 00303 return false; 00304 } 00305 00306 MsqMatrix< 3, 3 > A, W, dmdT, d2mdT2[6]; 00307 mf->jacobian( pd, e, bits, s, mIndices, mDerivs3D, num_idx, A, err ); 00308 MSQ_ERRZERO( err ); 00309 targetCalc->get_3D_target( pd, e, s, W, err ); 00310 MSQ_ERRZERO( err ); 00311 const MsqMatrix< 3, 3 > Winv = inverse( W ); 00312 const MsqMatrix< 3, 3 > T = A * Winv; 00313 rval = targetMetric->evaluate_with_hess( T, value, dmdT, d2mdT2, err ); 00314 MSQ_ERRZERO( err ); 00315 gradient< 3 >( num_idx, mDerivs3D, dmdT * transpose( Winv ), grad ); 00316 second_deriv_wrt_product_factor( d2mdT2, Winv ); 00317 00318 diagonal.resize( num_idx ); 00319 hessian_diagonal< 3 >( num_idx, mDerivs3D, d2mdT2, arrptr( diagonal ) ); 00320 #ifdef PRINT_INFO 00321 print_info< 3 >( e, s, A, W, A * inverse( W ) ); 00322 #endif 00323 } 00324 else if( edim == 2 ) 00325 { 00326 #ifdef NUMERICAL_2D_HESSIAN 00327 // use finite diference approximation for now 00328 return QualityMetric::evaluate_with_Hessian_diagonal( pd, p_handle, value, indices, grad, diagonal, err ); 00329 #else 00330 MsqMatrix< 2, 2 > W, A, dmdT, d2mdT2[3]; 00331 MsqMatrix< 3, 2 > M; 00332 rval = evaluate_surface_common( pd, s, e, bits, mIndices, num_idx, mDerivs2D, W, A, M, err ); 00333 if( MSQ_CHKERR( err ) || !rval ) return false; 00334 const MsqMatrix< 2, 2 > Winv = inverse( W ); 00335 const MsqMatrix< 2, 2 > T = A * Winv; 00336 rval = targetMetric->evaluate_with_hess( T, value, dmdT, d2mdT2, err ); 00337 MSQ_ERRZERO( err ); 00338 gradient< 2 >( num_idx, mDerivs2D, M * dmdT * transpose( Winv ), grad ); 00339 second_deriv_wrt_product_factor( d2mdT2, Winv ); 00340 00341 diagonal.resize( num_idx ); 00342 for( size_t i = 0; i < num_idx; ++i ) 00343 { 00344 MsqMatrix< 2, 2 > block2d; 00345 block2d( 0, 0 ) = transpose( mDerivs2D[i] ) * d2mdT2[0] * mDerivs2D[i]; 00346 block2d( 0, 1 ) = transpose( mDerivs2D[i] ) * d2mdT2[1] * mDerivs2D[i]; 00347 block2d( 1, 0 ) = block2d( 0, 1 ); 00348 block2d( 1, 1 ) = transpose( mDerivs2D[i] ) * d2mdT2[2] * mDerivs2D[i]; 00349 MsqMatrix< 3, 2 > p = M * block2d; 00350 00351 SymMatrix3D& H = diagonal[i]; 00352 H[0] = p.row( 0 ) * transpose( M.row( 0 ) ); 00353 H[1] = p.row( 0 ) * transpose( M.row( 1 ) ); 00354 H[2] = p.row( 0 ) * transpose( M.row( 2 ) ); 00355 H[3] = p.row( 1 ) * transpose( M.row( 1 ) ); 00356 H[4] = p.row( 1 ) * transpose( M.row( 2 ) ); 00357 H[5] = p.row( 2 ) * transpose( M.row( 2 ) ); 00358 } 00359 #ifdef PRINT_INFO 00360 print_info< 2 >( e, s, J, Wp, A * inverse( W ) ); 00361 #endif 00362 #endif 00363 } 00364 else 00365 { 00366 assert( 0 ); 00367 return false; 00368 } 00369 00370 // pass back index list 00371 indices.resize( num_idx ); 00372 std::copy( mIndices, mIndices + num_idx, indices.begin() ); 00373 00374 // apply target weight to value 00375 if( !num_idx ) 00376 weight( pd, s, e, num_idx, value, 0, 0, 0, err ); 00377 else 00378 weight( pd, s, e, num_idx, value, arrptr( grad ), arrptr( diagonal ), 0, err ); 00379 MSQ_ERRZERO( err ); 00380 return rval; 00381 } 00382 00383 } // namespace MBMesquite