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 /*! \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 }