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 IdealWeightInverseMeanRatio.cpp 00029 \brief 00030 00031 \author Michael Brewer 00032 \author Todd Munson 00033 \author Thomas Leurent 00034 00035 \date 2002-06-9 00036 */ 00037 #include "IdealWeightInverseMeanRatio.hpp" 00038 #include "MeanRatioFunctions.hpp" 00039 #include "MsqTimer.hpp" 00040 #include "MsqDebug.hpp" 00041 #include "MsqError.hpp" 00042 #include "PatchData.hpp" 00043 00044 #include <vector> 00045 using std::vector; 00046 #include <math.h> 00047 00048 namespace MBMesquite 00049 { 00050 00051 IdealWeightInverseMeanRatio::IdealWeightInverseMeanRatio( MsqError& err, double pow_dbl ) 00052 : AveragingQM( QualityMetric::LINEAR ) 00053 { 00054 set_metric_power( pow_dbl, err );MSQ_ERRRTN( err ); 00055 } 00056 00057 IdealWeightInverseMeanRatio::IdealWeightInverseMeanRatio() : AveragingQM( QualityMetric::LINEAR ) 00058 { 00059 MsqError err; 00060 set_metric_power( 1.0, err ); 00061 assert( !err ); 00062 } 00063 00064 std::string IdealWeightInverseMeanRatio::get_name() const 00065 { 00066 return "Inverse Mean Ratio"; 00067 } 00068 00069 int IdealWeightInverseMeanRatio::get_negate_flag() const 00070 { 00071 return b2Con.value() < 0 ? -1 : 1; 00072 } 00073 00074 //! Sets the power value in the metric computation. 00075 void IdealWeightInverseMeanRatio::set_metric_power( double pow_dbl, MsqError& err ) 00076 { 00077 if( fabs( pow_dbl ) <= MSQ_MIN ) 00078 { 00079 MSQ_SETERR( err )( MsqError::INVALID_ARG ); 00080 return; 00081 } 00082 a2Con = pow( .5, pow_dbl ); 00083 b2Con = pow_dbl; 00084 c2Con = -pow_dbl; 00085 a3Con = pow( 1.0 / 3.0, pow_dbl ); 00086 b3Con = pow_dbl; 00087 c3Con = -2.0 * pow_dbl / 3.0; 00088 } 00089 00090 bool IdealWeightInverseMeanRatio::evaluate( PatchData& pd, size_t handle, double& m, MsqError& err ) 00091 { 00092 const MsqMeshEntity* e = &pd.element_by_index( handle ); 00093 EntityTopology topo = e->get_element_type(); 00094 00095 const MsqVertex* vertices = pd.get_vertex_array( err ); 00096 MSQ_ERRZERO( err ); 00097 const size_t* v_i = e->get_vertex_index_array(); 00098 00099 Vector3D n; // Surface normal for 2D objects 00100 00101 // Prism and Hex element descriptions 00102 static const int locs_pri[6][4] = { { 0, 1, 2, 3 }, { 1, 2, 0, 4 }, { 2, 0, 1, 5 }, 00103 { 3, 5, 4, 0 }, { 4, 3, 5, 1 }, { 5, 4, 3, 2 } }; 00104 static const int locs_hex[8][4] = { { 0, 1, 3, 4 }, { 1, 2, 0, 5 }, { 2, 3, 1, 6 }, { 3, 0, 2, 7 }, 00105 { 4, 7, 5, 0 }, { 5, 4, 6, 1 }, { 6, 5, 7, 2 }, { 7, 6, 4, 3 } }; 00106 00107 const Vector3D d_con( 1.0, 1.0, 1.0 ); 00108 00109 int i; 00110 00111 m = 0.0; 00112 bool metric_valid = false; 00113 switch( topo ) 00114 { 00115 case TRIANGLE: 00116 pd.get_domain_normal_at_element( e, n, err ); 00117 MSQ_ERRZERO( err ); 00118 n = n / n.length(); // Need unit normal 00119 mCoords[0] = vertices[v_i[0]]; 00120 mCoords[1] = vertices[v_i[1]]; 00121 mCoords[2] = vertices[v_i[2]]; 00122 metric_valid = m_fcn_2e( m, mCoords, n, a2Con, b2Con, c2Con ); 00123 if( !metric_valid ) return false; 00124 break; 00125 00126 case QUADRILATERAL: 00127 pd.get_domain_normal_at_element( e, n, err ); 00128 MSQ_ERRZERO( err ); 00129 n = n / n.length(); // Need unit normal 00130 for( i = 0; i < 4; ++i ) 00131 { 00132 mCoords[0] = vertices[v_i[locs_hex[i][0]]]; 00133 mCoords[1] = vertices[v_i[locs_hex[i][1]]]; 00134 mCoords[2] = vertices[v_i[locs_hex[i][2]]]; 00135 metric_valid = m_fcn_2i( mMetrics[i], mCoords, n, a2Con, b2Con, c2Con, d_con ); 00136 if( !metric_valid ) return false; 00137 } 00138 m = average_metrics( mMetrics, 4, err ); 00139 break; 00140 00141 case TETRAHEDRON: 00142 mCoords[0] = vertices[v_i[0]]; 00143 mCoords[1] = vertices[v_i[1]]; 00144 mCoords[2] = vertices[v_i[2]]; 00145 mCoords[3] = vertices[v_i[3]]; 00146 metric_valid = m_fcn_3e( m, mCoords, a3Con, b3Con, c3Con ); 00147 if( !metric_valid ) return false; 00148 break; 00149 00150 case PYRAMID: 00151 for( i = 0; i < 4; ++i ) 00152 { 00153 mCoords[0] = vertices[v_i[i]]; 00154 mCoords[1] = vertices[v_i[( i + 1 ) % 4]]; 00155 mCoords[2] = vertices[v_i[( i + 3 ) % 4]]; 00156 mCoords[3] = vertices[v_i[4]]; 00157 metric_valid = m_fcn_3p( mMetrics[i], mCoords, a3Con, b3Con, c3Con ); 00158 if( !metric_valid ) return false; 00159 } 00160 m = average_metrics( mMetrics, 4, err ); 00161 MSQ_ERRZERO( err ); 00162 break; 00163 00164 case PRISM: 00165 for( i = 0; i < 6; ++i ) 00166 { 00167 mCoords[0] = vertices[v_i[locs_pri[i][0]]]; 00168 mCoords[1] = vertices[v_i[locs_pri[i][1]]]; 00169 mCoords[2] = vertices[v_i[locs_pri[i][2]]]; 00170 mCoords[3] = vertices[v_i[locs_pri[i][3]]]; 00171 metric_valid = m_fcn_3w( mMetrics[i], mCoords, a3Con, b3Con, c3Con ); 00172 if( !metric_valid ) return false; 00173 } 00174 m = average_metrics( mMetrics, 6, err ); 00175 MSQ_ERRZERO( err ); 00176 break; 00177 00178 case HEXAHEDRON: 00179 for( i = 0; i < 8; ++i ) 00180 { 00181 mCoords[0] = vertices[v_i[locs_hex[i][0]]]; 00182 mCoords[1] = vertices[v_i[locs_hex[i][1]]]; 00183 mCoords[2] = vertices[v_i[locs_hex[i][2]]]; 00184 mCoords[3] = vertices[v_i[locs_hex[i][3]]]; 00185 metric_valid = m_fcn_3i( mMetrics[i], mCoords, a3Con, b3Con, c3Con, d_con ); 00186 if( !metric_valid ) return false; 00187 } 00188 m = average_metrics( mMetrics, 8, err ); 00189 MSQ_ERRZERO( err ); 00190 break; 00191 00192 default: 00193 MSQ_SETERR( err ) 00194 ( MsqError::UNSUPPORTED_ELEMENT, "Element type (%d) not supported in IdealWeightInverseMeanRatio", 00195 (int)topo ); 00196 return false; 00197 } // end switch over element type 00198 return true; 00199 } 00200 00201 bool IdealWeightInverseMeanRatio::evaluate_with_gradient( PatchData& pd, size_t handle, double& m, 00202 std::vector< size_t >& indices, std::vector< Vector3D >& g, 00203 MsqError& err ) 00204 { 00205 const MsqMeshEntity* e = &pd.element_by_index( handle ); 00206 EntityTopology topo = e->get_element_type(); 00207 00208 if( !analytical_average_gradient() && topo != TRIANGLE && topo != TETRAHEDRON ) 00209 { 00210 static bool print = true; 00211 if( print ) 00212 { 00213 MSQ_DBGOUT( 1 ) << "Analyical gradient not available for selected averaging scheme. " 00214 << "Using (possibly much slower) numerical approximation of gradient" 00215 << " of quality metric. " << std::endl; 00216 print = false; 00217 } 00218 return QualityMetric::evaluate_with_gradient( pd, handle, m, indices, g, err ); 00219 } 00220 00221 const MsqVertex* vertices = pd.get_vertex_array( err ); 00222 MSQ_ERRZERO( err ); 00223 const size_t* v_i = e->get_vertex_index_array(); 00224 00225 Vector3D n; // Surface normal for 2D objects 00226 00227 // Prism and Hex element descriptions 00228 static const int locs_pri[6][4] = { { 0, 1, 2, 3 }, { 1, 2, 0, 4 }, { 2, 0, 1, 5 }, 00229 { 3, 5, 4, 0 }, { 4, 3, 5, 1 }, { 5, 4, 3, 2 } }; 00230 static const int locs_hex[8][4] = { { 0, 1, 3, 4 }, { 1, 2, 0, 5 }, { 2, 3, 1, 6 }, { 3, 0, 2, 7 }, 00231 { 4, 7, 5, 0 }, { 5, 4, 6, 1 }, { 6, 5, 7, 2 }, { 7, 6, 4, 3 } }; 00232 00233 const Vector3D d_con( 1.0, 1.0, 1.0 ); 00234 00235 int i; 00236 00237 bool metric_valid = false; 00238 const uint32_t fm = fixed_vertex_bitmap( pd, e, indices ); 00239 00240 m = 0.0; 00241 00242 switch( topo ) 00243 { 00244 case TRIANGLE: 00245 pd.get_domain_normal_at_element( e, n, err ); 00246 MSQ_ERRZERO( err ); 00247 n /= n.length(); // Need unit normal 00248 mCoords[0] = vertices[v_i[0]]; 00249 mCoords[1] = vertices[v_i[1]]; 00250 mCoords[2] = vertices[v_i[2]]; 00251 g.resize( 3 ); 00252 if( !g_fcn_2e( m, arrptr( g ), mCoords, n, a2Con, b2Con, c2Con ) ) return false; 00253 break; 00254 00255 case QUADRILATERAL: 00256 pd.get_domain_normal_at_element( e, n, err ); 00257 MSQ_ERRZERO( err ); 00258 n /= n.length(); // Need unit normal 00259 for( i = 0; i < 4; ++i ) 00260 { 00261 mCoords[0] = vertices[v_i[locs_hex[i][0]]]; 00262 mCoords[1] = vertices[v_i[locs_hex[i][1]]]; 00263 mCoords[2] = vertices[v_i[locs_hex[i][2]]]; 00264 if( !g_fcn_2i( mMetrics[i], mGradients + 3 * i, mCoords, n, a2Con, b2Con, c2Con, d_con ) ) return false; 00265 } 00266 00267 g.resize( 4 ); 00268 m = average_corner_gradients( QUADRILATERAL, fm, 4, mMetrics, mGradients, arrptr( g ), err ); 00269 MSQ_ERRZERO( err ); 00270 break; 00271 00272 case TETRAHEDRON: 00273 mCoords[0] = vertices[v_i[0]]; 00274 mCoords[1] = vertices[v_i[1]]; 00275 mCoords[2] = vertices[v_i[2]]; 00276 mCoords[3] = vertices[v_i[3]]; 00277 g.resize( 4 ); 00278 metric_valid = g_fcn_3e( m, arrptr( g ), mCoords, a3Con, b3Con, c3Con ); 00279 if( !metric_valid ) return false; 00280 break; 00281 00282 case PYRAMID: 00283 for( i = 0; i < 4; ++i ) 00284 { 00285 mCoords[0] = vertices[v_i[i]]; 00286 mCoords[1] = vertices[v_i[( i + 1 ) % 4]]; 00287 mCoords[2] = vertices[v_i[( i + 3 ) % 4]]; 00288 mCoords[3] = vertices[v_i[4]]; 00289 metric_valid = g_fcn_3p( mMetrics[i], mGradients + 4 * i, mCoords, a3Con, b3Con, c3Con ); 00290 if( !metric_valid ) return false; 00291 } 00292 00293 g.resize( 5 ); 00294 m = average_corner_gradients( PYRAMID, fm, 4, mMetrics, mGradients, arrptr( g ), err ); 00295 MSQ_ERRZERO( err ); 00296 break; 00297 00298 case PRISM: 00299 for( i = 0; i < 6; ++i ) 00300 { 00301 mCoords[0] = vertices[v_i[locs_pri[i][0]]]; 00302 mCoords[1] = vertices[v_i[locs_pri[i][1]]]; 00303 mCoords[2] = vertices[v_i[locs_pri[i][2]]]; 00304 mCoords[3] = vertices[v_i[locs_pri[i][3]]]; 00305 if( !g_fcn_3w( mMetrics[i], mGradients + 4 * i, mCoords, a3Con, b3Con, c3Con ) ) return false; 00306 } 00307 00308 g.resize( 6 ); 00309 m = average_corner_gradients( PRISM, fm, 6, mMetrics, mGradients, arrptr( g ), err ); 00310 MSQ_ERRZERO( err ); 00311 break; 00312 00313 case HEXAHEDRON: 00314 for( i = 0; i < 8; ++i ) 00315 { 00316 mCoords[0] = vertices[v_i[locs_hex[i][0]]]; 00317 mCoords[1] = vertices[v_i[locs_hex[i][1]]]; 00318 mCoords[2] = vertices[v_i[locs_hex[i][2]]]; 00319 mCoords[3] = vertices[v_i[locs_hex[i][3]]]; 00320 if( !g_fcn_3i( mMetrics[i], mGradients + 4 * i, mCoords, a3Con, b3Con, c3Con, d_con ) ) return false; 00321 } 00322 00323 g.resize( 8 ); 00324 m = average_corner_gradients( HEXAHEDRON, fm, 8, mMetrics, mGradients, arrptr( g ), err ); 00325 MSQ_ERRZERO( err ); 00326 break; 00327 00328 default: 00329 MSQ_SETERR( err ) 00330 ( MsqError::UNSUPPORTED_ELEMENT, "Element type (%d) not supported in IdealWeightInverseMeanRatio", 00331 (int)topo ); 00332 return false; 00333 } // end switch over element type 00334 00335 remove_fixed_gradients( topo, fm, g ); 00336 return true; 00337 } 00338 00339 bool IdealWeightInverseMeanRatio::evaluate_with_Hessian_diagonal( PatchData& pd, size_t handle, double& m, 00340 std::vector< size_t >& indices, 00341 std::vector< Vector3D >& g, 00342 std::vector< SymMatrix3D >& h, MsqError& err ) 00343 { 00344 const MsqMeshEntity* e = &pd.element_by_index( handle ); 00345 EntityTopology topo = e->get_element_type(); 00346 00347 if( !analytical_average_hessian() && topo != TRIANGLE && topo != TETRAHEDRON ) 00348 { 00349 static bool print = true; 00350 if( print ) 00351 { 00352 MSQ_DBGOUT( 1 ) << "Analyical gradient not available for selected averaging scheme. " 00353 << "Using (possibly much slower) numerical approximation of gradient" 00354 << " of quality metric. " << std::endl; 00355 print = false; 00356 } 00357 return QualityMetric::evaluate_with_Hessian_diagonal( pd, handle, m, indices, g, h, err ); 00358 } 00359 00360 const MsqVertex* vertices = pd.get_vertex_array( err ); 00361 MSQ_ERRZERO( err ); 00362 const size_t* v_i = e->get_vertex_index_array(); 00363 00364 Vector3D n; // Surface normal for 2D objects 00365 00366 // Prism and Hex element descriptions 00367 static const int locs_pri[6][4] = { { 0, 1, 2, 3 }, { 1, 2, 0, 4 }, { 2, 0, 1, 5 }, 00368 { 3, 5, 4, 0 }, { 4, 3, 5, 1 }, { 5, 4, 3, 2 } }; 00369 static const int locs_hex[8][4] = { { 0, 1, 3, 4 }, { 1, 2, 0, 5 }, { 2, 3, 1, 6 }, { 3, 0, 2, 7 }, 00370 { 4, 7, 5, 0 }, { 5, 4, 6, 1 }, { 6, 5, 7, 2 }, { 7, 6, 4, 3 } }; 00371 00372 const Vector3D d_con( 1.0, 1.0, 1.0 ); 00373 00374 int i; 00375 00376 bool metric_valid = false; 00377 const uint32_t fm = fixed_vertex_bitmap( pd, e, indices ); 00378 00379 m = 0.0; 00380 00381 switch( topo ) 00382 { 00383 case TRIANGLE: 00384 pd.get_domain_normal_at_element( e, n, err ); 00385 MSQ_ERRZERO( err ); 00386 n /= n.length(); // Need unit normal 00387 mCoords[0] = vertices[v_i[0]]; 00388 mCoords[1] = vertices[v_i[1]]; 00389 mCoords[2] = vertices[v_i[2]]; 00390 g.resize( 3 ); 00391 h.resize( 3 ); 00392 if( !h_fcn_2e( m, arrptr( g ), mHessians, mCoords, n, a2Con, b2Con, c2Con ) ) return false; 00393 h[0] = mHessians[0].upper(); 00394 h[1] = mHessians[3].upper(); 00395 h[2] = mHessians[5].upper(); 00396 break; 00397 00398 case QUADRILATERAL: 00399 pd.get_domain_normal_at_element( e, n, err ); 00400 MSQ_ERRZERO( err ); 00401 n /= n.length(); // Need unit normal 00402 for( i = 0; i < 4; ++i ) 00403 { 00404 mCoords[0] = vertices[v_i[locs_hex[i][0]]]; 00405 mCoords[1] = vertices[v_i[locs_hex[i][1]]]; 00406 mCoords[2] = vertices[v_i[locs_hex[i][2]]]; 00407 if( !h_fcn_2i( mMetrics[i], mGradients + 3 * i, mHessians + 6 * i, mCoords, n, a2Con, b2Con, c2Con, 00408 d_con ) ) 00409 return false; 00410 } 00411 00412 g.resize( 4 ); 00413 h.resize( 4 ); 00414 m = average_corner_hessian_diagonals( QUADRILATERAL, fm, 4, mMetrics, mGradients, mHessians, arrptr( g ), 00415 arrptr( h ), err ); 00416 MSQ_ERRZERO( err ); 00417 break; 00418 00419 case TETRAHEDRON: 00420 mCoords[0] = vertices[v_i[0]]; 00421 mCoords[1] = vertices[v_i[1]]; 00422 mCoords[2] = vertices[v_i[2]]; 00423 mCoords[3] = vertices[v_i[3]]; 00424 g.resize( 4 ); 00425 h.resize( 4 ); 00426 metric_valid = h_fcn_3e( m, arrptr( g ), mHessians, mCoords, a3Con, b3Con, c3Con ); 00427 if( !metric_valid ) return false; 00428 h[0] = mHessians[0].upper(); 00429 h[1] = mHessians[4].upper(); 00430 h[2] = mHessians[7].upper(); 00431 h[3] = mHessians[9].upper(); 00432 break; 00433 00434 case PYRAMID: 00435 for( i = 0; i < 4; ++i ) 00436 { 00437 mCoords[0] = vertices[v_i[i]]; 00438 mCoords[1] = vertices[v_i[( i + 1 ) % 4]]; 00439 mCoords[2] = vertices[v_i[( i + 3 ) % 4]]; 00440 mCoords[3] = vertices[v_i[4]]; 00441 metric_valid = 00442 h_fcn_3p( mMetrics[i], mGradients + 4 * i, mHessians + 10 * i, mCoords, a3Con, b3Con, c3Con ); 00443 00444 if( !metric_valid ) return false; 00445 } 00446 00447 g.resize( 5 ); 00448 h.resize( 5 ); 00449 m = average_corner_hessian_diagonals( PYRAMID, fm, 4, mMetrics, mGradients, mHessians, arrptr( g ), 00450 arrptr( h ), err ); 00451 MSQ_ERRZERO( err ); 00452 break; 00453 00454 case PRISM: 00455 for( i = 0; i < 6; ++i ) 00456 { 00457 mCoords[0] = vertices[v_i[locs_pri[i][0]]]; 00458 mCoords[1] = vertices[v_i[locs_pri[i][1]]]; 00459 mCoords[2] = vertices[v_i[locs_pri[i][2]]]; 00460 mCoords[3] = vertices[v_i[locs_pri[i][3]]]; 00461 if( !h_fcn_3w( mMetrics[i], mGradients + 4 * i, mHessians + 10 * i, mCoords, a3Con, b3Con, c3Con ) ) 00462 return false; 00463 } 00464 00465 g.resize( 6 ); 00466 h.resize( 6 ); 00467 m = average_corner_hessian_diagonals( PRISM, fm, 6, mMetrics, mGradients, mHessians, arrptr( g ), 00468 arrptr( h ), err ); 00469 MSQ_ERRZERO( err ); 00470 break; 00471 00472 case HEXAHEDRON: 00473 for( i = 0; i < 8; ++i ) 00474 { 00475 mCoords[0] = vertices[v_i[locs_hex[i][0]]]; 00476 mCoords[1] = vertices[v_i[locs_hex[i][1]]]; 00477 mCoords[2] = vertices[v_i[locs_hex[i][2]]]; 00478 mCoords[3] = vertices[v_i[locs_hex[i][3]]]; 00479 if( !h_fcn_3i( mMetrics[i], mGradients + 4 * i, mHessians + 10 * i, mCoords, a3Con, b3Con, c3Con, 00480 d_con ) ) 00481 return false; 00482 } 00483 00484 g.resize( 8 ); 00485 h.resize( 8 ); 00486 m = average_corner_hessian_diagonals( HEXAHEDRON, fm, 8, mMetrics, mGradients, mHessians, arrptr( g ), 00487 arrptr( h ), err ); 00488 MSQ_ERRZERO( err ); 00489 break; 00490 00491 default: 00492 MSQ_SETERR( err ) 00493 ( MsqError::UNSUPPORTED_ELEMENT, "Element type (%d) not supported in IdealWeightInverseMeanRatio", 00494 (int)topo ); 00495 return false; 00496 } // end switch over element type 00497 00498 remove_fixed_diagonals( topo, fm, g, h ); 00499 return true; 00500 } 00501 00502 bool IdealWeightInverseMeanRatio::evaluate_with_Hessian( PatchData& pd, size_t handle, double& m, 00503 std::vector< size_t >& indices, std::vector< Vector3D >& g, 00504 std::vector< Matrix3D >& h, MsqError& err ) 00505 { 00506 const MsqMeshEntity* e = &pd.element_by_index( handle ); 00507 EntityTopology topo = e->get_element_type(); 00508 00509 if( !analytical_average_hessian() && topo != TRIANGLE && topo != TETRAHEDRON ) 00510 { 00511 static bool print = true; 00512 if( print ) 00513 { 00514 MSQ_DBGOUT( 1 ) << "Analyical gradient not available for selected averaging scheme. " 00515 << "Using (possibly much slower) numerical approximation of gradient" 00516 << " of quality metric. " << std::endl; 00517 print = false; 00518 } 00519 return QualityMetric::evaluate_with_Hessian( pd, handle, m, indices, g, h, err ); 00520 } 00521 00522 const MsqVertex* vertices = pd.get_vertex_array( err ); 00523 MSQ_ERRZERO( err ); 00524 const size_t* v_i = e->get_vertex_index_array(); 00525 00526 Vector3D n; // Surface normal for 2D objects 00527 00528 // Prism and Hex element descriptions 00529 static const int locs_pri[6][4] = { { 0, 1, 2, 3 }, { 1, 2, 0, 4 }, { 2, 0, 1, 5 }, 00530 { 3, 5, 4, 0 }, { 4, 3, 5, 1 }, { 5, 4, 3, 2 } }; 00531 static const int locs_hex[8][4] = { { 0, 1, 3, 4 }, { 1, 2, 0, 5 }, { 2, 3, 1, 6 }, { 3, 0, 2, 7 }, 00532 { 4, 7, 5, 0 }, { 5, 4, 6, 1 }, { 6, 5, 7, 2 }, { 7, 6, 4, 3 } }; 00533 00534 const Vector3D d_con( 1.0, 1.0, 1.0 ); 00535 00536 int i; 00537 00538 bool metric_valid = false; 00539 const uint32_t fm = fixed_vertex_bitmap( pd, e, indices ); 00540 00541 m = 0.0; 00542 00543 switch( topo ) 00544 { 00545 case TRIANGLE: 00546 pd.get_domain_normal_at_element( e, n, err ); 00547 MSQ_ERRZERO( err ); 00548 n /= n.length(); // Need unit normal 00549 mCoords[0] = vertices[v_i[0]]; 00550 mCoords[1] = vertices[v_i[1]]; 00551 mCoords[2] = vertices[v_i[2]]; 00552 g.resize( 3 ); 00553 h.resize( 6 ); 00554 if( !h_fcn_2e( m, arrptr( g ), arrptr( h ), mCoords, n, a2Con, b2Con, c2Con ) ) return false; 00555 break; 00556 00557 case QUADRILATERAL: 00558 pd.get_domain_normal_at_element( e, n, err ); 00559 MSQ_ERRZERO( err ); 00560 n /= n.length(); // Need unit normal 00561 for( i = 0; i < 4; ++i ) 00562 { 00563 mCoords[0] = vertices[v_i[locs_hex[i][0]]]; 00564 mCoords[1] = vertices[v_i[locs_hex[i][1]]]; 00565 mCoords[2] = vertices[v_i[locs_hex[i][2]]]; 00566 if( !h_fcn_2i( mMetrics[i], mGradients + 3 * i, mHessians + 6 * i, mCoords, n, a2Con, b2Con, c2Con, 00567 d_con ) ) 00568 return false; 00569 } 00570 00571 g.resize( 4 ); 00572 h.resize( 10 ); 00573 m = average_corner_hessians( QUADRILATERAL, fm, 4, mMetrics, mGradients, mHessians, arrptr( g ), 00574 arrptr( h ), err ); 00575 MSQ_ERRZERO( err ); 00576 break; 00577 00578 case TETRAHEDRON: 00579 mCoords[0] = vertices[v_i[0]]; 00580 mCoords[1] = vertices[v_i[1]]; 00581 mCoords[2] = vertices[v_i[2]]; 00582 mCoords[3] = vertices[v_i[3]]; 00583 g.resize( 4 ); 00584 h.resize( 10 ); 00585 metric_valid = h_fcn_3e( m, arrptr( g ), arrptr( h ), mCoords, a3Con, b3Con, c3Con ); 00586 if( !metric_valid ) return false; 00587 break; 00588 00589 case PYRAMID: 00590 for( i = 0; i < 4; ++i ) 00591 { 00592 mCoords[0] = vertices[v_i[i]]; 00593 mCoords[1] = vertices[v_i[( i + 1 ) % 4]]; 00594 mCoords[2] = vertices[v_i[( i + 3 ) % 4]]; 00595 mCoords[3] = vertices[v_i[4]]; 00596 metric_valid = 00597 h_fcn_3p( mMetrics[i], mGradients + 4 * i, mHessians + 10 * i, mCoords, a3Con, b3Con, c3Con ); 00598 00599 if( !metric_valid ) return false; 00600 } 00601 00602 g.resize( 5 ); 00603 h.resize( 15 ); 00604 m = average_corner_hessians( PYRAMID, fm, 4, mMetrics, mGradients, mHessians, arrptr( g ), arrptr( h ), 00605 err ); 00606 MSQ_ERRZERO( err ); 00607 break; 00608 00609 case PRISM: 00610 for( i = 0; i < 6; ++i ) 00611 { 00612 mCoords[0] = vertices[v_i[locs_pri[i][0]]]; 00613 mCoords[1] = vertices[v_i[locs_pri[i][1]]]; 00614 mCoords[2] = vertices[v_i[locs_pri[i][2]]]; 00615 mCoords[3] = vertices[v_i[locs_pri[i][3]]]; 00616 if( !h_fcn_3w( mMetrics[i], mGradients + 4 * i, mHessians + 10 * i, mCoords, a3Con, b3Con, c3Con ) ) 00617 return false; 00618 } 00619 00620 g.resize( 6 ); 00621 h.resize( 21 ); 00622 m = average_corner_hessians( PRISM, fm, 6, mMetrics, mGradients, mHessians, arrptr( g ), arrptr( h ), err ); 00623 MSQ_ERRZERO( err ); 00624 break; 00625 00626 case HEXAHEDRON: 00627 for( i = 0; i < 8; ++i ) 00628 { 00629 mCoords[0] = vertices[v_i[locs_hex[i][0]]]; 00630 mCoords[1] = vertices[v_i[locs_hex[i][1]]]; 00631 mCoords[2] = vertices[v_i[locs_hex[i][2]]]; 00632 mCoords[3] = vertices[v_i[locs_hex[i][3]]]; 00633 if( !h_fcn_3i( mMetrics[i], mGradients + 4 * i, mHessians + 10 * i, mCoords, a3Con, b3Con, c3Con, 00634 d_con ) ) 00635 return false; 00636 } 00637 00638 g.resize( 8 ); 00639 h.resize( 36 ); 00640 m = average_corner_hessians( HEXAHEDRON, fm, 8, mMetrics, mGradients, mHessians, arrptr( g ), arrptr( h ), 00641 err ); 00642 MSQ_ERRZERO( err ); 00643 break; 00644 00645 default: 00646 MSQ_SETERR( err ) 00647 ( MsqError::UNSUPPORTED_ELEMENT, "Element type (%d) not supported in IdealWeightInverseMeanRatio", 00648 (int)topo ); 00649 return false; 00650 } // end switch over element type 00651 00652 remove_fixed_gradients( topo, fm, g ); 00653 remove_fixed_hessians( topo, fm, h ); 00654 return true; 00655 } 00656 00657 } // namespace MBMesquite