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