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