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