1 /* ***************************************************************** 2 MESQUITE -- The Mesh Quality Improvement Toolkit 3 4 Copyright 2004 Sandia Corporation and Argonne National 5 Laboratory. Under the terms of Contract DE-AC04-94AL85000 6 with Sandia Corporation, the U.S. Government retains certain 7 rights in this software. 8 9 This library is free software; you can redistribute it and/or 10 modify it under the terms of the GNU Lesser General Public 11 License as published by the Free Software Foundation; either 12 version 2.1 of the License, or (at your option) any later version. 13 14 This library is distributed in the hope that it will be useful, 15 but WITHOUT ANY WARRANTY; without even the implied warranty of 16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 17 Lesser General Public License for more details. 18 19 You should have received a copy of the GNU Lesser General Public License 20 (lgpl.txt) along with this library; if not, write to the Free Software 21 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 22 23 diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov, 24 pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov 25 26 (2006) kraftche@cae.wisc.edu 27 28 ***************************************************************** */ 29 30 /*! \file QualityMetric.hpp 31 \brief 32 Header file for the MBMesquite::QualityMetric class 33 34 \author Thomas Leurent 35 \author Michael Brewer 36 \date 2002-05-01 37 */ 38 39 #ifndef QualityMetric_hpp 40 #define QualityMetric_hpp 41 42 #include <cmath> 43 #include <vector> 44 #include <algorithm> 45 46 #include "Mesquite.hpp" 47 #include "Vector3D.hpp" 48 #include "Matrix3D.hpp" 49 50 #ifdef _MSC_VER 51 typedef unsigned uint32_t; 52 #elif defined(MOAB_HAVE_STDINT_H) 53 # include <stdint.h> 54 #elif defined(MOAB_HAVE_INTTYPES_H) 55 # include <inttypes.h> 56 #endif 57 58 namespace MBMesquite 59 { 60 61 /*! \class QualityMetric 62 \brief Base class for concrete quality metrics. 63 */ 64 class PatchData; 65 class MsqMeshEntity; 66 class Mesh; 67 class MeshDomain; 68 class MeshDomainAssoc; 69 class Settings; 70 71 class QualityMetric 72 { 73 protected: 74 QualityMetric()75 QualityMetric( ) : 76 keepFiniteDiffEps(false), 77 haveFiniteDiffEps(false) 78 {} 79 80 public: 81 82 enum MetricType 83 { 84 VERTEX_BASED, /**< Iterate over vertices to evaluate metric. */ 85 ELEMENT_BASED /**< Iterate over elements to evaluate metric. */ 86 }; 87 ~QualityMetric()88 MESQUITE_EXPORT virtual ~QualityMetric() 89 {} 90 91 MESQUITE_EXPORT virtual MetricType get_metric_type() const = 0; 92 93 MESQUITE_EXPORT virtual std::string get_name() const = 0; 94 95 //! 1 if metric should be minimized, -1 if metric should be maximized. 96 MESQUITE_EXPORT virtual int get_negate_flag() const = 0; 97 98 /**\brief Get locations at which metric can be evaluated 99 * 100 * Different metrics are evaluated for different things within 101 * a patch. For example, an element-based metric will be evaluated 102 * once for each element in patch, a vertex-based metric once for 103 * each veretx, and a target/sample-point based metric will be 104 * evaluated once for each samle point in each element. This method 105 * returns a list of handles, one for each location in the patch 106 * at which the metric can be evaluated. The handle values are used 107 * as input to the evaluate methods. 108 *\param pd The patch 109 *\param handles Output list of handles 110 *\param free_vertices_only If true, only pass back evaluation points 111 * that depend on at least one free vertex. 112 */ 113 MESQUITE_EXPORT virtual 114 void get_evaluations( PatchData& pd, 115 std::vector<size_t>& handles, 116 bool free_vertices_only, 117 MsqError& err ) = 0; 118 119 120 /**\brief Get locations at which metric can be evaluated for 121 * use in BCD intialization and QualityAssessor. 122 * 123 * For element-based, sample-based, and vertex-based metrics, 124 * this function is the same as get_evaluations. For edge-based 125 * metrics it returns only a subset of the results for get_evaluations 126 * such that each edge in the mesh is visited only once even though 127 * it would normally be visited twice when iterating over patches 128 * of the mesh. This assumes that no vertex occurs in more than one 129 * patch without its MSQ_PATCH_FIXED flag set. This assumption is true for 130 * both element-on-vertex and global patches. 131 *\param pd The patch 132 *\param handles Output list of handles 133 *\param free_vertices_only If true, only pass back evaluation points 134 * that depend on at least one free vertex. 135 */ 136 MESQUITE_EXPORT virtual 137 void get_single_pass( PatchData& pd, 138 std::vector<size_t>& handles, 139 bool free_vertices_only, 140 MsqError& err ); 141 142 /**\brief Get metric value at a logical location in the patch. 143 * 144 * Evaluate the metric at one location in the PatchData. 145 *\param pd The patch. 146 *\param handle The location in the patch (as passed back from get_evaluations). 147 *\param value The output metric value. 148 */ 149 MESQUITE_EXPORT virtual 150 bool evaluate( PatchData& pd, 151 size_t handle, 152 double& value, 153 MsqError& err ) = 0; 154 155 /**\brief Get metric value at a logical location in the patch. 156 * 157 * Evaluate the metric at one location in the PatchData. 158 *\param pd The patch. 159 *\param handle The location in the patch (as passed back from get_evaluations). 160 *\param value The output metric value. 161 *\param indices The free vertices that the evaluation is a function 162 * of, specified as vertex indices in the PatchData. 163 */ 164 MESQUITE_EXPORT virtual 165 bool evaluate_with_indices( PatchData& pd, 166 size_t handle, 167 double& value, 168 std::vector<size_t>& indices, 169 MsqError& err ) = 0; 170 171 /**\brief Get metric value and gradient at a logical location in the patch. 172 * 173 * Evaluate the metric at one location in the PatchData. 174 *\param pd The patch. 175 *\param handle The location in the patch (as passed back from get_evaluations). 176 *\param value The output metric value. 177 *\param indices The free vertices that the evaluation is a function 178 * of, specified as vertex indices in the PatchData. 179 *\param gradient The gradient of the metric as a function of the 180 * coordinates of the free vertices passed back in 181 * the indices list. 182 */ 183 MESQUITE_EXPORT virtual 184 bool evaluate_with_gradient( PatchData& pd, 185 size_t handle, 186 double& value, 187 std::vector<size_t>& indices, 188 std::vector<Vector3D>& gradient, 189 MsqError& err ); 190 191 /**\brief Get metric value and gradient at a logical location in the patch. 192 * 193 * Evaluate the metric at one location in the PatchData. 194 *\param pd The patch. 195 *\param handle The location in the patch (as passed back from get_evaluations). 196 *\param value The output metric value. 197 *\param indices The free vertices that the evaluation is a function 198 * of, specified as vertex indices in the PatchData. 199 *\param gradient The gradient of the metric as a function of the 200 * coordinates of the free vertices passed back in 201 * the indices list. 202 *\param Hessian_diagonal The 3x3 blocks along the diagonal of 203 * the Hessian matrix. 204 */ 205 MESQUITE_EXPORT virtual 206 bool evaluate_with_Hessian_diagonal( PatchData& pd, 207 size_t handle, 208 double& value, 209 std::vector<size_t>& indices, 210 std::vector<Vector3D>& gradient, 211 std::vector<SymMatrix3D>& Hessian_diagonal, 212 MsqError& err ); 213 214 /**\brief Get metric value and deravitives at a logical location in the patch. 215 * 216 * Evaluate the metric at one location in the PatchData. 217 *\param pd The patch. 218 *\param handle The location in the patch (as passed back from get_evaluations). 219 *\param value The output metric value. 220 *\param indices The free vertices that the evaluation is a function 221 * of, specified as vertex indices in the PatchData. 222 *\param gradient The gradient of the metric as a function of the 223 * coordinates of the free vertices passed back in 224 * the indices list. 225 *\param Hessian The Hessian of the metric as a function of the 226 * coordinates. The Hessian is passed back as the 227 * upper-triangular portion of the matrix in row-major 228 * order, where each Matrix3D is the portion of the 229 * Hessian with respect to the vertices at the 230 * corresponding positions in the indices list. 231 */ 232 MESQUITE_EXPORT virtual 233 bool evaluate_with_Hessian( PatchData& pd, 234 size_t handle, 235 double& value, 236 std::vector<size_t>& indices, 237 std::vector<Vector3D>& gradient, 238 std::vector<Matrix3D>& Hessian, 239 MsqError& err ); 240 241 //!Escobar Barrier Function for Shape and Other Metrics 242 // det = signed determinant of Jacobian Matrix at a Vertex 243 // delta = scaling parameter vertex_barrier_function(double det,double delta)244 static inline double vertex_barrier_function(double det, double delta) 245 { return 0.5*(det+sqrt(det*det+4*delta*delta)); } 246 //protected: 247 248 /** \brief Remove from vector any gradient terms corresponding 249 * to a fixed vertex. 250 * 251 * Remove terms from vector that correspond to fixed vertices. 252 *\param type Element type 253 *\param fixed_vertices Bit flags, one per vertex, 1 if 254 * vertex is fixed. 255 *\param gradients Array of gradients 256 */ 257 MESQUITE_EXPORT 258 static void remove_fixed_gradients( EntityTopology type, 259 uint32_t fixed_vertices, 260 std::vector<Vector3D>& gradients ); 261 262 /** \brief Remove from vectors any gradient terms and hessian 263 * diagonal blcoks corresponding to a fixed vertex. 264 * 265 * Remove terms from vector that correspond to fixed vertices. 266 *\param type Element type 267 *\param fixed_vertices Bit flags, one per vertex, 1 if 268 * vertex is fixed. 269 *\param gradients Array of gradients 270 *\param hess_diagonal_blocks Array of diagonal blocks of Hessian matrix. 271 */ 272 MESQUITE_EXPORT 273 static void remove_fixed_diagonals( EntityTopology type, 274 uint32_t fixed_vertices, 275 std::vector<Vector3D>& gradients, 276 std::vector<SymMatrix3D>& hess_diagonal_blocks ); 277 278 /** \brief Remove from vector any Hessian blocks corresponding 279 * to a fixed vertex. 280 * 281 * Remove blocks from vector that correspond to fixed vertices. 282 *\param type Element type 283 *\param fixed_vertices Bit flags, one per vertex, 1 if 284 * vertex is fixed. 285 *\param hessians Array of Hessian blocks (upper trianguler, row-major) 286 */ 287 MESQUITE_EXPORT 288 static void remove_fixed_hessians ( EntityTopology type, 289 uint32_t fixed_vertices, 290 std::vector<Matrix3D>& hessians ); 291 292 /** \brief Convert fixed vertex format from list to bit flags 293 * 294 * Given list of pointers to fixed vertices as passed to 295 * evaluation functions, convert to bit flag format used 296 * for many utility functions in this class. Bits correspond 297 * to vertices in the canonical vertex order, beginning with 298 * the least-significant bit. The bit is cleared for free 299 * vertices and set (1) for fixed vertices. 300 */ 301 MESQUITE_EXPORT 302 static uint32_t fixed_vertex_bitmap( PatchData& pd, 303 const MsqMeshEntity* elem, 304 std::vector<size_t>& free_indices ); 305 306 307 //! takes an array of coefficients and an array of metrics (both of length num_value) 308 //! and averages the contents using averaging method 'method'. 309 MESQUITE_EXPORT 310 double weighted_average_metrics(const double coef[], 311 const double metric_values[], 312 const int& num_values, MsqError &err); 313 314 /*!AveragingMethod allows you to set how the quality metric values 315 attained at each sample point will be averaged together to produce 316 a single metric value for an element. 317 */ 318 enum AveragingMethod 319 { 320 LINEAR, //!< the linear average 321 RMS, //!< the root-mean-squared average 322 HMS, //!< the harmonic-mean-squared average 323 SUM, //!< the sum of the values 324 SUM_SQUARED, //!< the sum of the squares of the values 325 HARMONIC, //!< the harmonic average 326 LAST_WITH_HESSIAN=HARMONIC, 327 MINIMUM, //!< the minimum value 328 MAXIMUM, //!< the maximum value 329 GEOMETRIC, //!< the geometric average 330 LAST_WITH_GRADIENT=GEOMETRIC, 331 STANDARD_DEVIATION, //!< the standard deviation squared of the values 332 MAX_OVER_MIN, //!< the maximum value minus the minum value 333 MAX_MINUS_MIN, //!< the maximum value divided by the minimum value 334 SUM_OF_RATIOS_SQUARED //!< (1/(N^2))*(SUM (SUM (v_i/v_j)^2)) 335 }; 336 337 //!\brief Called at start of instruction queue processing 338 //! 339 //! Do any preliminary global initialization, consistency checking, 340 //! etc. Default implementation does nothing. 341 MESQUITE_EXPORT virtual 342 void initialize_queue( MeshDomainAssoc* mesh_and_domain, 343 const Settings* settings, 344 MsqError& err ); 345 346 private: 347 int feasible; 348 349 std::vector<Matrix3D> tmpHess; 350 bool keepFiniteDiffEps; //!< True if gradient finite difference 351 //!< calculation should set \c finiteDiffEps 352 bool haveFiniteDiffEps; //!< True if finite difference Hessian code 353 //!< has calculated \c finiteDiffEps 354 double finiteDiffEps; //!< Location for finite difference Hessian code 355 //!< to store this value so that it doesn't need 356 //!< to be recalculated if the gradient calculation 357 //!< is also finite difference 358 }; 359 360 361 } //namespace 362 363 364 #endif // QualityMetric_hpp 365