MOAB: Mesh Oriented datABase  (version 5.4.1)
UntangleBetaQualityMetric.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     [email protected], [email protected], [email protected],
00024     [email protected], [email protected], [email protected]
00025 
00026   ***************************************************************** */
00027 /*! \file   UntangleBetaQualityMetric.cpp
00028   UntangleBeta is an untangle quality metric which can be used to evaluate
00029   the quality of two- or three-dimensional elements.
00030 
00031   \author Michael Brewer
00032   \date   2002-10-10
00033 */
00034 
00035 #include "UntangleBetaQualityMetric.hpp"
00036 #include "Vector3D.hpp"
00037 #include "PatchData.hpp"
00038 #include "MsqError.hpp"
00039 
00040 using namespace MBMesquite;
00041 
00042 static inline void untangle_function_2d( double beta,
00043                                          const Vector3D temp_vec[],
00044                                          size_t e_ind,
00045                                          PatchData& pd,
00046                                          double& fval,
00047                                          MsqError& err )
00048 {
00049     Vector3D surface_normal;
00050     pd.get_domain_normal_at_element( e_ind, surface_normal, err );MSQ_ERRRTN( err );
00051     Vector3D cross_vec = temp_vec[0] * temp_vec[1];
00052     // cout<<"\nsurface_normal "<<surface_normal;
00053     // cout<<"\cross_vec "<<cross_vec;
00054     double temp_var = cross_vec.length();
00055     if( cross_vec % surface_normal < 0.0 )
00056     {
00057         temp_var *= -1;
00058     }
00059     temp_var -= beta;
00060     // cout<<"temp_var == "<<temp_var;
00061     fval = 0.0;
00062     if( temp_var < 0.0 )
00063     {
00064         fval = fabs( temp_var ) - temp_var;
00065     }
00066     //  cout<<"\nfval == "<<fval<<"  e_ind "<<e_ind;
00067 }
00068 
00069 static inline void untangle_function_3d( double beta, const Vector3D temp_vec[], double& fval )
00070 {
00071     double temp_var = temp_vec[0] % ( temp_vec[1] * temp_vec[2] );
00072     temp_var -= beta;
00073     fval = 0.0;
00074     if( temp_var < 0.0 )
00075     {
00076         fval = fabs( temp_var ) - temp_var;
00077     }
00078 }
00079 
00080 /*! \fn UntangleBetaQualityMetric::UntangleBetaQualityMetric(double bet)
00081   \brief For untangle beta, the constructor defaults to the SUM
00082   averaging method, and to the ELEMENT_VERTICES evaluation mode.
00083 */
00084 UntangleBetaQualityMetric::UntangleBetaQualityMetric( double bet ) : AveragingQM( RMS ), mBeta( bet ) {}
00085 
00086 std::string UntangleBetaQualityMetric::get_name() const
00087 {
00088     return "Untangle Beta";
00089 }
00090 
00091 int UntangleBetaQualityMetric::get_negate_flag() const
00092 {
00093     return 1;
00094 }
00095 
00096 bool UntangleBetaQualityMetric::evaluate( PatchData& pd, size_t e_ind, double& fval, MsqError& err )
00097 {
00098 
00099     double met_vals[MSQ_MAX_NUM_VERT_PER_ENT];
00100     fval                         = MSQ_MAX_CAP;
00101     const MsqMeshEntity* element = &pd.element_by_index( e_ind );
00102     const size_t* v_i            = element->get_vertex_index_array();
00103     // only 3 temp_vec will be sent to untangle calculator, but the
00104     // additional vector3Ds may be needed during the calculations
00105     Vector3D temp_vec[5];
00106     const MsqVertex* vertices = pd.get_vertex_array( err );
00107     MSQ_ERRZERO( err );
00108     EntityTopology type = element->get_element_type();
00109     switch( type )
00110     {
00111         case TRIANGLE:
00112             temp_vec[0] = vertices[v_i[1]] - vertices[v_i[0]];
00113             temp_vec[2] = vertices[v_i[2]] - vertices[v_i[0]];
00114             // make relative to equilateral
00115             temp_vec[1] = ( ( 2 * temp_vec[2] ) - temp_vec[0] ) * MSQ_SQRT_THREE_INV;
00116             untangle_function_2d( mBeta, temp_vec, e_ind, pd, fval, err );
00117             break;
00118         case QUADRILATERAL:
00119             temp_vec[0] = vertices[v_i[1]] - vertices[v_i[0]];
00120             temp_vec[1] = vertices[v_i[3]] - vertices[v_i[0]];
00121             untangle_function_2d( mBeta, temp_vec, e_ind, pd, met_vals[0], err );
00122             MSQ_ERRZERO( err );
00123 
00124             temp_vec[0] = vertices[v_i[2]] - vertices[v_i[1]];
00125             temp_vec[1] = vertices[v_i[0]] - vertices[v_i[1]];
00126             untangle_function_2d( mBeta, temp_vec, e_ind, pd, met_vals[1], err );
00127             MSQ_ERRZERO( err );
00128 
00129             temp_vec[0] = vertices[v_i[3]] - vertices[v_i[2]];
00130             temp_vec[1] = vertices[v_i[1]] - vertices[v_i[2]];
00131             untangle_function_2d( mBeta, temp_vec, e_ind, pd, met_vals[2], err );
00132             MSQ_ERRZERO( err );
00133 
00134             temp_vec[0] = vertices[v_i[0]] - vertices[v_i[3]];
00135             temp_vec[1] = vertices[v_i[2]] - vertices[v_i[3]];
00136             untangle_function_2d( mBeta, temp_vec, e_ind, pd, met_vals[3], err );
00137             MSQ_ERRZERO( err );
00138             fval = average_metrics( met_vals, 4, err );
00139             MSQ_ERRZERO( err );
00140             break;
00141         case TETRAHEDRON:
00142             temp_vec[0] = vertices[v_i[1]] - vertices[v_i[0]];
00143             temp_vec[3] = vertices[v_i[2]] - vertices[v_i[0]];
00144             temp_vec[4] = vertices[v_i[3]] - vertices[v_i[0]];
00145             // transform to equilateral tet
00146             temp_vec[1] = ( ( 2 * temp_vec[3] ) - temp_vec[0] ) / MSQ_SQRT_THREE;
00147             temp_vec[2] = ( ( 3 * temp_vec[4] ) - temp_vec[0] - temp_vec[3] ) / ( MSQ_SQRT_THREE * MSQ_SQRT_TWO );
00148             untangle_function_3d( mBeta, temp_vec, fval );
00149             break;
00150         case HEXAHEDRON:
00151             // transform to v_i[0]
00152             temp_vec[0] = vertices[v_i[1]] - vertices[v_i[0]];
00153             temp_vec[1] = vertices[v_i[3]] - vertices[v_i[0]];
00154             temp_vec[2] = vertices[v_i[4]] - vertices[v_i[0]];
00155             untangle_function_3d( mBeta, temp_vec, met_vals[0] );
00156 
00157             temp_vec[0] = vertices[v_i[2]] - vertices[v_i[1]];
00158             temp_vec[1] = vertices[v_i[0]] - vertices[v_i[1]];
00159             temp_vec[2] = vertices[v_i[5]] - vertices[v_i[1]];
00160             untangle_function_3d( mBeta, temp_vec, met_vals[1] );
00161 
00162             temp_vec[0] = vertices[v_i[3]] - vertices[v_i[2]];
00163             temp_vec[1] = vertices[v_i[1]] - vertices[v_i[2]];
00164             temp_vec[2] = vertices[v_i[6]] - vertices[v_i[2]];
00165             untangle_function_3d( mBeta, temp_vec, met_vals[2] );
00166 
00167             temp_vec[0] = vertices[v_i[0]] - vertices[v_i[3]];
00168             temp_vec[1] = vertices[v_i[2]] - vertices[v_i[3]];
00169             temp_vec[2] = vertices[v_i[7]] - vertices[v_i[3]];
00170             untangle_function_3d( mBeta, temp_vec, met_vals[3] );
00171 
00172             temp_vec[0] = vertices[v_i[7]] - vertices[v_i[4]];
00173             temp_vec[1] = vertices[v_i[5]] - vertices[v_i[4]];
00174             temp_vec[2] = vertices[v_i[0]] - vertices[v_i[4]];
00175             untangle_function_3d( mBeta, temp_vec, met_vals[4] );
00176 
00177             temp_vec[0] = vertices[v_i[4]] - vertices[v_i[5]];
00178             temp_vec[1] = vertices[v_i[6]] - vertices[v_i[5]];
00179             temp_vec[2] = vertices[v_i[1]] - vertices[v_i[5]];
00180             untangle_function_3d( mBeta, temp_vec, met_vals[5] );
00181 
00182             temp_vec[0] = vertices[v_i[5]] - vertices[v_i[6]];
00183             temp_vec[1] = vertices[v_i[7]] - vertices[v_i[6]];
00184             temp_vec[2] = vertices[v_i[2]] - vertices[v_i[6]];
00185             untangle_function_3d( mBeta, temp_vec, met_vals[6] );
00186 
00187             temp_vec[0] = vertices[v_i[6]] - vertices[v_i[7]];
00188             temp_vec[1] = vertices[v_i[4]] - vertices[v_i[7]];
00189             temp_vec[2] = vertices[v_i[3]] - vertices[v_i[7]];
00190             untangle_function_3d( mBeta, temp_vec, met_vals[7] );
00191             fval = average_metrics( met_vals, 8, err );
00192             MSQ_ERRZERO( err );
00193             break;
00194         case PYRAMID:
00195             for( unsigned i = 0; i < 4; ++i )
00196             {
00197                 temp_vec[0] = vertices[v_i[( i + 1 ) % 4]] - vertices[v_i[i]];
00198                 temp_vec[1] = vertices[v_i[( i + 3 ) % 4]] - vertices[v_i[i]];
00199                 temp_vec[3] = vertices[v_i[4]] - vertices[v_i[i]];
00200                 temp_vec[2] = ( 4 * temp_vec[3] - temp_vec[0] - temp_vec[1] ) / ( 2.0 - MSQ_SQRT_TWO );
00201                 untangle_function_3d( mBeta, temp_vec, met_vals[i] );
00202             }
00203             fval = average_metrics( met_vals, 4, err );
00204             MSQ_ERRZERO( err );
00205             break;
00206         case PRISM:
00207             temp_vec[0] = vertices[v_i[1]] - vertices[v_i[0]];
00208             temp_vec[3] = vertices[v_i[2]] - vertices[v_i[0]];
00209             temp_vec[1] = ( ( 2 * temp_vec[3] ) - temp_vec[0] ) * MSQ_SQRT_THREE_INV;
00210             temp_vec[2] = vertices[v_i[3]] - vertices[v_i[0]];
00211             untangle_function_3d( mBeta, temp_vec, met_vals[0] );
00212 
00213             temp_vec[0] = vertices[v_i[2]] - vertices[v_i[1]];
00214             temp_vec[3] = vertices[v_i[0]] - vertices[v_i[1]];
00215             temp_vec[1] = ( ( 2 * temp_vec[3] ) - temp_vec[0] ) * MSQ_SQRT_THREE_INV;
00216             temp_vec[2] = vertices[v_i[4]] - vertices[v_i[1]];
00217             untangle_function_3d( mBeta, temp_vec, met_vals[1] );
00218 
00219             temp_vec[0] = vertices[v_i[0]] - vertices[v_i[2]];
00220             temp_vec[3] = vertices[v_i[1]] - vertices[v_i[2]];
00221             temp_vec[1] = ( ( 2 * temp_vec[3] ) - temp_vec[0] ) * MSQ_SQRT_THREE_INV;
00222             temp_vec[2] = vertices[v_i[5]] - vertices[v_i[2]];
00223             untangle_function_3d( mBeta, temp_vec, met_vals[2] );
00224 
00225             temp_vec[0] = vertices[v_i[5]] - vertices[v_i[3]];
00226             temp_vec[3] = vertices[v_i[4]] - vertices[v_i[3]];
00227             temp_vec[1] = ( ( 2 * temp_vec[3] ) - temp_vec[0] ) * MSQ_SQRT_THREE_INV;
00228             temp_vec[2] = vertices[v_i[0]] - vertices[v_i[3]];
00229             untangle_function_3d( mBeta, temp_vec, met_vals[3] );
00230 
00231             temp_vec[0] = vertices[v_i[3]] - vertices[v_i[4]];
00232             temp_vec[3] = vertices[v_i[5]] - vertices[v_i[4]];
00233             temp_vec[1] = ( ( 2 * temp_vec[3] ) - temp_vec[0] ) * MSQ_SQRT_THREE_INV;
00234             temp_vec[2] = vertices[v_i[1]] - vertices[v_i[4]];
00235             untangle_function_3d( mBeta, temp_vec, met_vals[4] );
00236 
00237             temp_vec[0] = vertices[v_i[4]] - vertices[v_i[5]];
00238             temp_vec[3] = vertices[v_i[3]] - vertices[v_i[5]];
00239             temp_vec[1] = ( ( 2 * temp_vec[3] ) - temp_vec[0] ) * MSQ_SQRT_THREE_INV;
00240             temp_vec[2] = vertices[v_i[2]] - vertices[v_i[5]];
00241             untangle_function_3d( mBeta, temp_vec, met_vals[5] );
00242 
00243             fval = average_metrics( met_vals, 6, err );
00244             MSQ_ERRZERO( err );
00245             break;
00246         default:
00247             MSQ_SETERR( err )
00248             ( MsqError::UNSUPPORTED_ELEMENT, "Unsupported cell type (%d) for Untangle quality metric.", type );
00249 
00250             fval = MSQ_MAX_CAP;
00251             return false;
00252     }  // end switch over element type
00253     return true;
00254 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines