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