MOAB: Mesh Oriented datABase  (version 5.3.1)
ConditionNumberQualityMetric.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   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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines