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