MOAB: Mesh Oriented datABase  (version 5.4.0)
LocalSizeQualityMetric.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 LocalSizeQualityMetric.cpp
00028   \author Michael Brewer
00029   \date April 9, 2003
00030   Evaluates the corner volume (or areas) of the element corners
00031   attached to a given vertiex and then averages those values
00032   together.
00033 */
00034 
00035 #include "LocalSizeQualityMetric.hpp"
00036 #include "PatchData.hpp"
00037 
00038 using namespace MBMesquite;
00039 
00040 //! Calculate the area of the triangle formed by the three vertices.
00041 static inline double compute_corner_area( PatchData& pd, size_t vert_1, size_t vert_2, size_t vert_3, MsqError& err )
00042 {
00043     const MsqVertex* verts = pd.get_vertex_array( err );
00044     Vector3D vec_1         = verts[vert_2] - verts[vert_1];
00045     Vector3D vec_2         = verts[vert_3] - verts[vert_1];
00046     Vector3D cross_vec     = vec_1 * vec_2;
00047     return ( cross_vec.length() / 2.0 );
00048 }
00049 
00050 //! Calculate the volume of the tetrahedron formed by the four vertices.
00051 static inline double compute_corner_volume( PatchData& pd,
00052                                             size_t vert_1,
00053                                             size_t vert_2,
00054                                             size_t vert_3,
00055                                             size_t vert_4,
00056                                             MsqError& err )
00057 {
00058     const MsqVertex* verts = pd.get_vertex_array( err );
00059     Vector3D vec_1         = verts[vert_2] - verts[vert_1];
00060     Vector3D vec_2         = verts[vert_3] - verts[vert_1];
00061     Vector3D vec_3         = verts[vert_4] - verts[vert_1];
00062     return fabs( ( vec_3 % ( vec_1 * vec_2 ) ) / 6.0 );
00063 }
00064 
00065 LocalSizeQualityMetric::~LocalSizeQualityMetric() {}
00066 
00067 std::string LocalSizeQualityMetric::get_name() const
00068 {
00069     return "Local Size";
00070 }
00071 
00072 int LocalSizeQualityMetric::get_negate_flag() const
00073 {
00074     return 1;
00075 }
00076 
00077 /*!For the given vertex, vert, with connected elements, e_i for i=1...K,
00078   the LocalSizeQualityMetric computes the corner volumes (or areas) of
00079   each e_i at the corner defined by vert.  The corner volume is defined
00080   as the volume of the tet defined by the edges of an element which contain
00081   the common vertex, vert.  That volume is then diveded by the average corner
00082   volume of all the element corners connected to this vertex.  For
00083   vertices attached to pyramid elements, this metric is undefined.
00084 */
00085 bool LocalSizeQualityMetric::evaluate( PatchData& pd, size_t this_vert, double& fval, MsqError& err )
00086 {
00087     fval = 0.0;
00088     // get the element array
00089     MsqMeshEntity* elems = pd.get_element_array( err );
00090     MSQ_ERRZERO( err );
00091     // get the vertex to element array and the offset array
00092     // const size_t* elem_offset = pd.get_vertex_to_elem_offset(err);  MSQ_ERRZERO(err);
00093     // const size_t* v_to_e_array = pd.get_vertex_to_elem_array(err);  MSQ_ERRZERO(err);
00094     // find the offset for this vertex
00095     // size_t this_offset = elem_offset[this_vert];
00096     // get the number of elements attached to this vertex (given by the
00097     // first entry in the vertex to element array)
00098     // size_t num_elems = v_to_e_array[this_offset];
00099     // PRINT_INFO("\nIN LOCAL SIZE CPP, num_elements = %i",num_elems);
00100     size_t num_elems;
00101     const size_t* v_to_e_array = pd.get_vertex_element_adjacencies( this_vert, num_elems, err );
00102     MSQ_ERRZERO( err );
00103 
00104     if( num_elems <= 0 )
00105     {
00106         return true;
00107     }
00108 
00109     // create an array to store the local metric values before averaging
00110     // Can we remove this dynamic allocatio?
00111     double* met_vals = new double[num_elems];
00112     // vector to hold the other verts which form a corner.
00113     std::vector< size_t > other_vertices;
00114     other_vertices.reserve( 4 );
00115     double total_val = 0.0;
00116     size_t i         = 0;
00117     // loop over the elements attached to this vertex
00118     for( i = 0; i < num_elems; ++i )
00119     {
00120         // get the vertices which (with this_vert) form the corner of
00121         // the ith element.
00122         elems[v_to_e_array[i]].get_connected_vertices( this_vert, other_vertices, err );
00123         MSQ_ERRZERO( err );
00124         ////PRINT_INFO("\nINSIDE LOCAL SIZE CPP other_vertices size = %i",other_vertices.size());
00125 
00126         switch( other_vertices.size() )
00127         {
00128                 // if a surface element, compute the corner area
00129             case 2:
00130                 met_vals[i] = compute_corner_area( pd, this_vert, other_vertices[0], other_vertices[1], err );
00131                 MSQ_ERRZERO( err );
00132                 break;
00133                 // if a volume element, compute the corner volume
00134             case 3:
00135                 met_vals[i] = compute_corner_volume( pd, this_vert, other_vertices[0], other_vertices[1],
00136                                                      other_vertices[2], err );
00137                 MSQ_ERRZERO( err );
00138                 break;
00139             default:
00140                 // otherwise, there is was an error.  Either the wrong number
00141                 // of vertices were returned fom get_connected_vertices or
00142                 // the element does not have the correct number of edges
00143                 // connected to this vertex (possibly a pyramid element).
00144                 met_vals[i] = 0.0;
00145                 MSQ_SETERR( err )
00146                 ( "Incorrect number of vertices returned from "
00147                   "get_connected_vertices.",
00148                   MsqError::UNSUPPORTED_ELEMENT );
00149                 return false;
00150         };
00151         // keep track of total so that we can compute the linear average
00152         total_val += met_vals[i];
00153         // PRINT_INFO("\nIN LOCAL SIZE CPP, total_val = %f, i = %i",total_val,i);
00154         // clear the vector of other_vertices for re-use.
00155         other_vertices.clear();
00156         // PRINT_INFO("\nIN LOCAL SIZE CPP, after clean size = %f",other_vertices.size());
00157     }
00158     // calculate the linear average... num_elems is non-zero here.
00159     total_val /= (double)num_elems;
00160     // PRINT_INFO("\nIN LOCAL SIZE CPP, average = %f",total_val);
00161     // if the average is non-zero
00162     // divide each entry by the linear average
00163     if( total_val != 0 )
00164     {
00165         for( i = 0; i < num_elems; ++i )
00166         {
00167             met_vals[i] /= total_val;
00168         }
00169         // calculate fval by averaging the corner values
00170         fval = average_metrics( met_vals, num_elems, err );
00171         MSQ_ERRZERO( err );
00172         // PRINT_INFO("\nIN LOCAL SIZE CPP, inside if statement");
00173     }
00174     // PRINT_INFO("\nIN LOCAL SIZE CPP, fval = %f",fval);
00175     // clean up the dynamically allocated array
00176     delete[] met_vals;
00177     // always return true... the vertex position is never invalid
00178     return true;
00179 }
00180 
00181 bool LocalSizeQualityMetric::evaluate_with_indices( PatchData& pd,
00182                                                     size_t vertex,
00183                                                     double& value,
00184                                                     std::vector< size_t >& indices,
00185                                                     MsqError& err )
00186 {
00187     indices.clear();
00188     pd.get_adjacent_vertex_indices( vertex, indices, err );
00189     MSQ_ERRZERO( err );
00190 
00191     std::vector< size_t >::iterator r, w;
00192     for( r = w = indices.begin(); r != indices.end(); ++r )
00193     {
00194         if( *r < pd.num_free_vertices() )
00195         {
00196             *w = *r;
00197             ++w;
00198         }
00199     }
00200     indices.erase( w, indices.end() );
00201     if( vertex < pd.num_free_vertices() ) indices.push_back( vertex );
00202 
00203     bool rval = evaluate( pd, vertex, value, err );
00204     return !MSQ_CHKERR( err ) && rval;
00205 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines