MOAB: Mesh Oriented datABase  (version 5.4.0)
QualityMetric.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     kraftche@cae.wisc.edu
00026 
00027   ***************************************************************** */
00028 /*!
00029   \file   QualityMetric.cpp
00030   \brief
00031 
00032   \author Michael Brewer
00033   \author Thomas Leurent
00034   \date   2002-05-14
00035   \author Jason Kraftcheck
00036   \date   2006-04-20
00037 */
00038 
00039 #include "QualityMetric.hpp"
00040 #include "MsqVertex.hpp"
00041 #include "PatchData.hpp"
00042 
00043 namespace MBMesquite
00044 {
00045 
00046 void QualityMetric::initialize_queue( MeshDomainAssoc*, const Settings*, MsqError& ) {}
00047 
00048 void QualityMetric::get_single_pass( PatchData& pd,
00049                                      std::vector< size_t >& handles,
00050                                      bool free_vertices_only,
00051                                      MsqError& err )
00052 {
00053     get_evaluations( pd, handles, free_vertices_only, err );
00054 }
00055 
00056 static inline double get_delta_C( const PatchData& pd, const std::vector< size_t >& indices, MsqError& err )
00057 {
00058     if( indices.empty() )
00059     {
00060         MSQ_SETERR( err )( MsqError::INVALID_ARG );
00061         return 0.0;
00062     }
00063 
00064     std::vector< size_t >::const_iterator beg, iter, iter2, end;
00065 
00066     std::vector< size_t > tmp_vect;
00067     if( indices.size() == 1u )
00068     {
00069         pd.get_adjacent_vertex_indices( indices.front(), tmp_vect, err );
00070         MSQ_ERRZERO( err );
00071         assert( !tmp_vect.empty() );
00072         tmp_vect.push_back( indices.front() );
00073         beg = tmp_vect.begin();
00074         end = tmp_vect.end();
00075     }
00076     else
00077     {
00078         beg = indices.begin();
00079         end = indices.end();
00080     }
00081 
00082     double min_dist_sqr = HUGE_VAL, sum_dist_sqr = 0.0;
00083     for( iter = beg; iter != end; ++iter )
00084     {
00085         for( iter2 = iter + 1; iter2 != end; ++iter2 )
00086         {
00087             Vector3D diff = pd.vertex_by_index( *iter );
00088             diff -= pd.vertex_by_index( *iter2 );
00089             double dist_sqr = diff % diff;
00090             if( dist_sqr < min_dist_sqr ) min_dist_sqr = dist_sqr;
00091             sum_dist_sqr += dist_sqr;
00092         }
00093     }
00094 
00095     return 3e-6 * sqrt( min_dist_sqr ) + 5e-7 * sqrt( sum_dist_sqr / ( end - beg ) );
00096 }
00097 
00098 bool QualityMetric::evaluate_with_gradient( PatchData& pd,
00099                                             size_t handle,
00100                                             double& value,
00101                                             std::vector< size_t >& indices,
00102                                             std::vector< Vector3D >& gradient,
00103                                             MsqError& err )
00104 {
00105     indices.clear();
00106     bool valid = evaluate_with_indices( pd, handle, value, indices, err );
00107     if( MSQ_CHKERR( err ) || !valid ) return false;
00108     if( indices.empty() ) return true;
00109 
00110     // get initial pertubation amount
00111     double delta_C = finiteDiffEps;
00112     if( !haveFiniteDiffEps )
00113     {
00114         delta_C = get_delta_C( pd, indices, err );
00115         MSQ_ERRZERO( err );
00116         if( keepFiniteDiffEps )
00117         {
00118             finiteDiffEps     = delta_C;
00119             haveFiniteDiffEps = true;
00120         }
00121     }
00122     const double delta_inv_C  = 1.0 / delta_C;
00123     const int reduction_limit = 15;
00124 
00125     gradient.resize( indices.size() );
00126     for( size_t v = 0; v < indices.size(); ++v )
00127     {
00128         const Vector3D pos = pd.vertex_by_index( indices[v] );
00129 
00130         /* gradient in the x, y, z direction */
00131         for( int j = 0; j < 3; ++j )
00132         {
00133             double delta     = delta_C;
00134             double delta_inv = delta_inv_C;
00135             double metric_value;
00136             Vector3D delta_v( 0, 0, 0 );
00137 
00138             // perturb the node and calculate gradient.  The while loop is a
00139             // safety net to make sure the epsilon perturbation does not take
00140             // the element out of the feasible region.
00141             int counter = 0;
00142             for( ;; )
00143             {
00144                 // perturb the coordinates of the free vertex in the j direction
00145                 // by delta
00146                 delta_v[j] = delta;
00147                 pd.set_vertex_coordinates( pos + delta_v, indices[v], err );
00148                 MSQ_ERRZERO( err );
00149 
00150                 // compute the function at the perturbed point location
00151                 valid = evaluate( pd, handle, metric_value, err );
00152                 if( !MSQ_CHKERR( err ) && valid ) break;
00153 
00154                 if( ++counter >= reduction_limit )
00155                 {
00156                     MSQ_SETERR( err )
00157                     ( "Perturbing vertex by delta caused an inverted element.", MsqError::INTERNAL_ERROR );
00158                     return false;
00159                 }
00160 
00161                 delta *= 0.1;
00162                 delta_inv *= 10.;
00163             }
00164             // put the coordinates back where they belong
00165             pd.set_vertex_coordinates( pos, indices[v], err );
00166             // compute the numerical gradient
00167             gradient[v][j] = ( metric_value - value ) * delta_inv;
00168         }  // for(j)
00169     }      // for(indices)
00170     return true;
00171 }
00172 
00173 bool QualityMetric::evaluate_with_Hessian( PatchData& pd,
00174                                            size_t handle,
00175                                            double& value,
00176                                            std::vector< size_t >& indices,
00177                                            std::vector< Vector3D >& gradient,
00178                                            std::vector< Matrix3D >& Hessian,
00179                                            MsqError& err )
00180 {
00181     indices.clear();
00182     gradient.clear();
00183     keepFiniteDiffEps = true;
00184     bool valid        = evaluate_with_gradient( pd, handle, value, indices, gradient, err );
00185     keepFiniteDiffEps = false;
00186     if( MSQ_CHKERR( err ) || !valid )
00187     {
00188         haveFiniteDiffEps = false;
00189         return false;
00190     }
00191     if( indices.empty() )
00192     {
00193         haveFiniteDiffEps = false;
00194         return true;
00195     }
00196 
00197     // get initial pertubation amount
00198     double delta_C;
00199     if( haveFiniteDiffEps )
00200     {
00201         delta_C = finiteDiffEps;
00202     }
00203     else
00204     {
00205         delta_C = get_delta_C( pd, indices, err );
00206         MSQ_ERRZERO( err );
00207     }
00208     assert( delta_C < 1e30 );
00209     const double delta_inv_C  = 1.0 / delta_C;
00210     const int reduction_limit = 15;
00211 
00212     std::vector< Vector3D > temp_gradient( indices.size() );
00213     const int num_hess = indices.size() * ( indices.size() + 1 ) / 2;
00214     Hessian.resize( num_hess );
00215 
00216     for( unsigned v = 0; v < indices.size(); ++v )
00217     {
00218         const Vector3D pos = pd.vertex_by_index( indices[v] );
00219         for( int j = 0; j < 3; ++j )
00220         {  // x, y, and z
00221             double delta     = delta_C;
00222             double delta_inv = delta_inv_C;
00223             double metric_value;
00224             Vector3D delta_v( 0, 0, 0 );
00225 
00226             // find finite difference for gradient
00227             int counter = 0;
00228             for( ;; )
00229             {
00230                 delta_v[j] = delta;
00231                 pd.set_vertex_coordinates( pos + delta_v, indices[v], err );
00232                 MSQ_ERRZERO( err );
00233                 valid = evaluate_with_gradient( pd, handle, metric_value, indices, temp_gradient, err );
00234                 if( !MSQ_CHKERR( err ) && valid ) break;
00235 
00236                 if( ++counter >= reduction_limit )
00237                 {
00238                     MSQ_SETERR( err )
00239                     ( "Algorithm did not successfully compute element's "
00240                       "Hessian.\n",
00241                       MsqError::INTERNAL_ERROR );
00242                     haveFiniteDiffEps = false;
00243                     return false;
00244                 }
00245 
00246                 delta *= 0.1;
00247                 delta_inv *= 10.0;
00248             }
00249             pd.set_vertex_coordinates( pos, indices[v], err );
00250             MSQ_ERRZERO( err );
00251 
00252             // compute the numerical Hessian
00253             for( unsigned w = 0; w <= v; ++w )
00254             {
00255                 // finite difference to get some entries of the Hessian
00256                 Vector3D fd( temp_gradient[w] );
00257                 fd -= gradient[w];
00258                 fd *= delta_inv;
00259                 // For the block at position w,v in a matrix, we need the corresponding index
00260                 // (mat_index) in a 1D array containing only upper triangular blocks.
00261                 unsigned sum_w           = w * ( w + 1 ) / 2;  // 1+2+3+...+w
00262                 unsigned mat_index       = w * indices.size() + v - sum_w;
00263                 Hessian[mat_index][0][j] = fd[0];
00264                 Hessian[mat_index][1][j] = fd[1];
00265                 Hessian[mat_index][2][j] = fd[2];
00266             }
00267         }  // for(j)
00268     }      // for(indices)
00269     haveFiniteDiffEps = false;
00270     return true;
00271 }
00272 
00273 bool QualityMetric::evaluate_with_Hessian_diagonal( PatchData& pd,
00274                                                     size_t handle,
00275                                                     double& value,
00276                                                     std::vector< size_t >& indices,
00277                                                     std::vector< Vector3D >& gradient,
00278                                                     std::vector< SymMatrix3D >& Hessian_diagonal,
00279                                                     MsqError& err )
00280 {
00281     bool rval = evaluate_with_Hessian( pd, handle, value, indices, gradient, tmpHess, err );
00282     if( MSQ_CHKERR( err ) || !rval ) return rval;
00283     size_t s = indices.size();
00284     Hessian_diagonal.resize( s );
00285     std::vector< Matrix3D >::const_iterator h = tmpHess.begin();
00286     for( size_t i = 0; i < indices.size(); ++i )
00287     {
00288         Hessian_diagonal[i] = h->upper();
00289         h += s--;
00290     }
00291     return rval;
00292 }
00293 
00294 uint32_t QualityMetric::fixed_vertex_bitmap( PatchData& pd, const MsqMeshEntity* elem, std::vector< size_t >& indices )
00295 {
00296     indices.clear();
00297     uint32_t result        = ~(uint32_t)0;
00298     unsigned num_vtx       = elem->vertex_count();
00299     const size_t* vertices = elem->get_vertex_index_array();
00300     indices.clear();
00301     for( unsigned i = 0; i < num_vtx; ++i )
00302     {
00303         if( vertices[i] < pd.num_free_vertices() )
00304         {
00305             indices.push_back( vertices[i] );
00306             result &= ~(uint32_t)( 1 << i );
00307         }
00308     }
00309     return result;
00310 }
00311 
00312 void QualityMetric::remove_fixed_gradients( EntityTopology elem_type, uint32_t fixed, std::vector< Vector3D >& grads )
00313 {
00314     const unsigned num_vertex = TopologyInfo::corners( elem_type );
00315     unsigned r, w;
00316     for( r = 0; r < num_vertex && !( fixed & ( 1 << r ) ); ++r )
00317         ;
00318     for( w = r++; r < num_vertex; ++r )
00319     {
00320         if( !( fixed & ( 1 << r ) ) )
00321         {
00322             grads[w] = grads[r];
00323             ++w;
00324         }
00325     }
00326     grads.resize( w );
00327 }
00328 
00329 void QualityMetric::remove_fixed_diagonals( EntityTopology type,
00330                                             uint32_t fixed,
00331                                             std::vector< Vector3D >& grads,
00332                                             std::vector< SymMatrix3D >& diags )
00333 {
00334     const unsigned num_vertex = TopologyInfo::corners( type );
00335     unsigned r, w;
00336     for( r = 0; r < num_vertex && !( fixed & ( 1 << r ) ); ++r )
00337         ;
00338     for( w = r++; r < num_vertex; ++r )
00339     {
00340         if( !( fixed & ( 1 << r ) ) )
00341         {
00342             grads[w] = grads[r];
00343             diags[w] = diags[r];
00344             ++w;
00345         }
00346     }
00347     grads.resize( w );
00348     diags.resize( w );
00349 }
00350 
00351 void QualityMetric::remove_fixed_hessians( EntityTopology elem_type, uint32_t fixed, std::vector< Matrix3D >& hessians )
00352 {
00353     const unsigned num_vertex = TopologyInfo::corners( elem_type );
00354     unsigned r, c, i = 0, w = 0;
00355     for( r = 0; r < num_vertex; ++r )
00356     {
00357         if( fixed & ( 1 << r ) )
00358         {
00359             i += num_vertex - r;
00360             continue;
00361         }
00362         for( c = r; c < num_vertex; ++c )
00363         {
00364             if( !( fixed & ( 1 << c ) ) )
00365             {
00366                 hessians[w] = hessians[i];
00367                 ++w;
00368             }
00369             ++i;
00370         }
00371     }
00372     hessians.resize( w );
00373 }
00374 
00375 double QualityMetric::weighted_average_metrics( const double coef[],
00376                                                 const double metric_values[],
00377                                                 const int& num_values,
00378                                                 MsqError& /*err*/ )
00379 {
00380     // MSQ_MAX needs to be made global?
00381     // double MSQ_MAX=1e10;
00382     double total_value = 0.0;
00383     int i              = 0;
00384     // if no values, return zero
00385     if( num_values <= 0 )
00386     {
00387         return 0.0;
00388     }
00389 
00390     for( i = 0; i < num_values; ++i )
00391     {
00392         total_value += coef[i] * metric_values[i];
00393     }
00394     total_value /= (double)num_values;
00395 
00396     return total_value;
00397 }
00398 
00399 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines