MOAB: Mesh Oriented datABase  (version 5.2.1)
IdealWeightInverseMeanRatio.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines