|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.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 [email protected], [email protected], [email protected], 00024 [email protected], [email protected], [email protected] 00025 00026 ***************************************************************** */ 00027 /*! 00028 \file ConditionNumberQualityMetric.cpp 00029 \brief 00030 00031 \author Michael Brewer 00032 \date 2002-06-9 00033 */ 00034 #include <vector> 00035 #include "ConditionNumberQualityMetric.hpp" 00036 #include <cmath> 00037 #include "Vector3D.hpp" 00038 #include "ConditionNumberFunctions.hpp" 00039 00040 using namespace MBMesquite; 00041 00042 ConditionNumberQualityMetric::ConditionNumberQualityMetric() : AveragingQM( QualityMetric::LINEAR ) {} 00043 00044 std::string ConditionNumberQualityMetric::get_name() const 00045 { 00046 return "Condition Number"; 00047 } 00048 00049 int ConditionNumberQualityMetric::get_negate_flag() const 00050 { 00051 return 1; 00052 } 00053 00054 bool ConditionNumberQualityMetric::evaluate( PatchData& pd, size_t handle, double& fval, MsqError& err ) 00055 { 00056 const MsqMeshEntity* const element = &pd.element_by_index( handle ); 00057 bool return_flag; 00058 double met_vals[MSQ_MAX_NUM_VERT_PER_ENT]; 00059 fval = MSQ_MAX_CAP; 00060 const size_t* v_i = element->get_vertex_index_array(); 00061 // only 3 temp_vec will be sent to cond-num calculator, but the 00062 // additional vector3Ds may be needed during the calculations 00063 Vector3D temp_vec[6]; 00064 const MsqVertex* vertices = pd.get_vertex_array( err ); 00065 EntityTopology type = element->get_element_type(); 00066 switch( type ) 00067 { 00068 case TRIANGLE: 00069 temp_vec[0] = vertices[v_i[1]] - vertices[v_i[0]]; 00070 temp_vec[2] = vertices[v_i[2]] - vertices[v_i[0]]; 00071 // make relative to equilateral 00072 temp_vec[1] = ( ( 2 * temp_vec[2] ) - temp_vec[0] ) * MSQ_SQRT_THREE_INV; 00073 return_flag = condition_number_2d( temp_vec, handle, pd, fval, err ); 00074 MSQ_ERRZERO( err ); 00075 return return_flag; 00076 case QUADRILATERAL: 00077 temp_vec[0] = vertices[v_i[1]] - vertices[v_i[0]]; 00078 temp_vec[1] = vertices[v_i[3]] - vertices[v_i[0]]; 00079 return_flag = condition_number_2d( temp_vec, handle, pd, met_vals[0], err ); 00080 MSQ_ERRZERO( err ); 00081 if( !return_flag ) return return_flag; 00082 temp_vec[0] = vertices[v_i[2]] - vertices[v_i[1]]; 00083 temp_vec[1] = vertices[v_i[0]] - vertices[v_i[1]]; 00084 return_flag = condition_number_2d( temp_vec, handle, pd, met_vals[1], err ); 00085 MSQ_ERRZERO( err ); 00086 if( !return_flag ) return return_flag; 00087 temp_vec[0] = vertices[v_i[3]] - vertices[v_i[2]]; 00088 temp_vec[1] = vertices[v_i[1]] - vertices[v_i[2]]; 00089 return_flag = condition_number_2d( temp_vec, handle, pd, met_vals[2], err ); 00090 MSQ_ERRZERO( err ); 00091 if( !return_flag ) return return_flag; 00092 temp_vec[0] = vertices[v_i[0]] - vertices[v_i[3]]; 00093 temp_vec[1] = vertices[v_i[2]] - vertices[v_i[3]]; 00094 return_flag = condition_number_2d( temp_vec, handle, pd, met_vals[3], err ); 00095 MSQ_ERRZERO( err ); 00096 fval = average_metrics( met_vals, 4, err ); 00097 return return_flag; 00098 case TETRAHEDRON: 00099 temp_vec[0] = vertices[v_i[1]] - vertices[v_i[0]]; 00100 temp_vec[3] = vertices[v_i[2]] - vertices[v_i[0]]; 00101 temp_vec[4] = vertices[v_i[3]] - vertices[v_i[0]]; 00102 // transform to equilateral tet 00103 temp_vec[1] = ( ( 2 * temp_vec[3] ) - temp_vec[0] ) / MSQ_SQRT_THREE; 00104 temp_vec[2] = ( ( 3 * temp_vec[4] ) - temp_vec[0] - temp_vec[3] ) / ( MSQ_SQRT_THREE * MSQ_SQRT_TWO ); 00105 return_flag = condition_number_3d( temp_vec, pd, fval, err ); 00106 MSQ_ERRZERO( err ); 00107 return return_flag; 00108 00109 case HEXAHEDRON: 00110 // transform to v_i[0] 00111 temp_vec[0] = vertices[v_i[1]] - vertices[v_i[0]]; 00112 temp_vec[1] = vertices[v_i[3]] - vertices[v_i[0]]; 00113 temp_vec[2] = vertices[v_i[4]] - vertices[v_i[0]]; 00114 return_flag = condition_number_3d( temp_vec, pd, met_vals[0], err ); 00115 MSQ_ERRZERO( err ); 00116 if( !return_flag ) return return_flag; 00117 temp_vec[0] = vertices[v_i[2]] - vertices[v_i[1]]; 00118 temp_vec[1] = vertices[v_i[0]] - vertices[v_i[1]]; 00119 temp_vec[2] = vertices[v_i[5]] - vertices[v_i[1]]; 00120 return_flag = condition_number_3d( temp_vec, pd, met_vals[1], err ); 00121 MSQ_ERRZERO( err ); 00122 if( !return_flag ) return return_flag; 00123 temp_vec[0] = vertices[v_i[3]] - vertices[v_i[2]]; 00124 temp_vec[1] = vertices[v_i[1]] - vertices[v_i[2]]; 00125 temp_vec[2] = vertices[v_i[6]] - vertices[v_i[2]]; 00126 return_flag = condition_number_3d( temp_vec, pd, met_vals[2], err ); 00127 MSQ_ERRZERO( err ); 00128 if( !return_flag ) return return_flag; 00129 temp_vec[0] = vertices[v_i[0]] - vertices[v_i[3]]; 00130 temp_vec[1] = vertices[v_i[2]] - vertices[v_i[3]]; 00131 temp_vec[2] = vertices[v_i[7]] - vertices[v_i[3]]; 00132 return_flag = condition_number_3d( temp_vec, pd, met_vals[3], err ); 00133 MSQ_ERRZERO( err ); 00134 if( !return_flag ) return return_flag; 00135 temp_vec[0] = vertices[v_i[7]] - vertices[v_i[4]]; 00136 temp_vec[1] = vertices[v_i[5]] - vertices[v_i[4]]; 00137 temp_vec[2] = vertices[v_i[0]] - vertices[v_i[4]]; 00138 return_flag = condition_number_3d( temp_vec, pd, met_vals[4], err ); 00139 MSQ_ERRZERO( err ); 00140 if( !return_flag ) return return_flag; 00141 temp_vec[0] = vertices[v_i[4]] - vertices[v_i[5]]; 00142 temp_vec[1] = vertices[v_i[6]] - vertices[v_i[5]]; 00143 temp_vec[2] = vertices[v_i[1]] - vertices[v_i[5]]; 00144 return_flag = condition_number_3d( temp_vec, pd, met_vals[5], err ); 00145 MSQ_ERRZERO( err ); 00146 if( !return_flag ) return return_flag; 00147 temp_vec[0] = vertices[v_i[5]] - vertices[v_i[6]]; 00148 temp_vec[1] = vertices[v_i[7]] - vertices[v_i[6]]; 00149 temp_vec[2] = vertices[v_i[2]] - vertices[v_i[6]]; 00150 return_flag = condition_number_3d( temp_vec, pd, met_vals[6], err ); 00151 MSQ_ERRZERO( err ); 00152 if( !return_flag ) return return_flag; 00153 temp_vec[0] = vertices[v_i[6]] - vertices[v_i[7]]; 00154 temp_vec[1] = vertices[v_i[4]] - vertices[v_i[7]]; 00155 temp_vec[2] = vertices[v_i[3]] - vertices[v_i[7]]; 00156 return_flag = condition_number_3d( temp_vec, pd, met_vals[7], err ); 00157 MSQ_ERRZERO( err ); 00158 fval = average_metrics( met_vals, 8, err ); 00159 MSQ_ERRZERO( err ); 00160 return return_flag; 00161 00162 case PYRAMID: { 00163 unsigned num_adj; 00164 const unsigned* adj_idx; 00165 return_flag = true; 00166 for( size_t j = 0; j < 4; ++j ) // skip fifth vertex (apex) 00167 { 00168 adj_idx = TopologyInfo::adjacent_vertices( type, j, num_adj ); 00169 assert( num_adj == 3 ); 00170 00171 temp_vec[0] = vertices[v_i[adj_idx[0]]] - vertices[v_i[j]]; 00172 temp_vec[1] = vertices[v_i[adj_idx[1]]] - vertices[v_i[j]]; 00173 // calculate last vect map to right tetrahedron 00174 temp_vec[3] = vertices[v_i[adj_idx[2]]] - vertices[v_i[adj_idx[0]]]; 00175 temp_vec[4] = vertices[v_i[adj_idx[2]]] - vertices[v_i[adj_idx[1]]]; 00176 temp_vec[2] = 0.5 * ( temp_vec[3] + temp_vec[4] ); 00177 00178 return_flag = return_flag && condition_number_3d( temp_vec, pd, met_vals[j], err ); 00179 } 00180 fval = average_metrics( met_vals, 4, err ); 00181 MSQ_ERRZERO( err ); 00182 return return_flag; 00183 } 00184 00185 case PRISM: { 00186 unsigned num_adj; 00187 const unsigned* adj_idx; 00188 return_flag = true; 00189 for( size_t j = 0; j < 6; ++j ) 00190 { 00191 adj_idx = TopologyInfo::adjacent_vertices( type, j, num_adj ); 00192 assert( num_adj == 3 ); 00193 00194 temp_vec[0] = vertices[v_i[adj_idx[0]]] - vertices[v_i[j]]; 00195 temp_vec[1] = vertices[v_i[adj_idx[1]]] - vertices[v_i[j]]; 00196 temp_vec[2] = vertices[v_i[adj_idx[2]]] - vertices[v_i[j]]; 00197 // map to right tetrahedron 00198 temp_vec[1] += vertices[v_i[adj_idx[1]]]; 00199 temp_vec[1] -= vertices[v_i[adj_idx[0]]]; 00200 temp_vec[1] *= MSQ_SQRT_THREE_INV; 00201 00202 return_flag = return_flag && condition_number_3d( temp_vec, pd, met_vals[j], err ); 00203 } 00204 fval = average_metrics( met_vals, 6, err ); 00205 MSQ_ERRZERO( err ); 00206 return return_flag; 00207 } 00208 00209 default: 00210 MSQ_SETERR( err ) 00211 ( MsqError::UNSUPPORTED_ELEMENT, "Unsupported cell type (%d) for Condition Number quality metric.", type ); 00212 00213 fval = MSQ_MAX_CAP; 00214 return false; 00215 } // end switch over element type 00216 return false; 00217 }
1.7.6.1 
