1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
/* *****************************************************************
    MESQUITE -- The Mesh Quality Improvement Toolkit

    Copyright 2004 Sandia Corporation and Argonne National
    Laboratory.  Under the terms of Contract DE-AC04-94AL85000
    with Sandia Corporation, the U.S. Government retains certain
    rights in this software.

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    (lgpl.txt) along with this library; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    [email protected], [email protected], [email protected],
    [email protected], [email protected], [email protected]

    (2006) [email protected]

  ***************************************************************** */

/*! \file QualityMetric.hpp
    \brief
Header file for the MBMesquite::QualityMetric class

  \author Thomas Leurent
  \author Michael Brewer
  \date   2002-05-01
 */

#ifndef QualityMetric_hpp
#define QualityMetric_hpp

#include <cmath>
#include <vector>
#include <algorithm>

#include "Mesquite.hpp"
#include "Vector3D.hpp"
#include "Matrix3D.hpp"

#ifdef _WIN32
typedef unsigned uint32_t;
#elif defined( MOAB_HAVE_STDINT_H )
#include <stdint.h>
#elif defined( MOAB_HAVE_INTTYPES_H )
#include <inttypes.h>
#endif

namespace MBMesquite
{

/*! \class QualityMetric
  \brief Base class for concrete quality metrics.
*/
class PatchData;
class MsqMeshEntity;
class Mesh;
class MeshDomain;
class MeshDomainAssoc;
class Settings;

class QualityMetric
{
  protected:
    QualityMetric() : keepFiniteDiffEps( false ), haveFiniteDiffEps( false ) {}<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.<--- Member variable 'QualityMetric::feasible' is not initialized in the constructor.<--- Member variable 'QualityMetric::finiteDiffEps' is not initialized in the constructor.

  public:
    enum MetricType
    {
        VERTEX_BASED, /**< Iterate over vertices to evaluate metric. */
        ELEMENT_BASED /**< Iterate over elements to evaluate metric. */
    };

    MESQUITE_EXPORT virtual ~QualityMetric() {}

    MESQUITE_EXPORT virtual MetricType get_metric_type() const = 0;

    MESQUITE_EXPORT virtual std::string get_name() const = 0;

    //! 1 if metric should be minimized, -1 if metric should be maximized.
    MESQUITE_EXPORT virtual int get_negate_flag() const = 0;

    /**\brief Get locations at which metric can be evaluated
     *
     * Different metrics are evaluated for different things within
     * a patch.  For example, an element-based metric will be evaluated
     * once for each element in patch, a vertex-based metric once for
     * each veretx, and a target/sample-point based metric will be
     * evaluated once for each samle point in each element.  This method
     * returns a list of handles, one for each location in the patch
     * at which the metric can be evaluated.  The handle values are used
     * as input to the evaluate methods.
     *\param pd       The patch
     *\param handles  Output list of handles
     *\param free_vertices_only If true, only pass back evaluation points
     *         that depend on at least one free vertex.
     */
    MESQUITE_EXPORT virtual void get_evaluations( PatchData& pd,
                                                  std::vector< size_t >& handles,
                                                  bool free_vertices_only,
                                                  MsqError& err ) = 0;

    /**\brief Get locations at which metric can be evaluated for
     *        use in BCD intialization and QualityAssessor.
     *
     * For element-based, sample-based, and vertex-based metrics,
     * this function is the same as get_evaluations.  For edge-based
     * metrics it returns only a subset of the results for get_evaluations
     * such that each edge in the mesh is visited only once even though
     * it would normally be visited twice when iterating over patches
     * of the mesh.  This assumes that no vertex occurs in more than one
     * patch without its MSQ_PATCH_FIXED flag set.  This assumption is true for
     * both element-on-vertex and global patches.
     *\param pd       The patch
     *\param handles  Output list of handles
     *\param free_vertices_only If true, only pass back evaluation points
     *         that depend on at least one free vertex.
     */
    MESQUITE_EXPORT virtual void get_single_pass( PatchData& pd,
                                                  std::vector< size_t >& handles,
                                                  bool free_vertices_only,
                                                  MsqError& err );

    /**\brief Get metric value at a logical location in the patch.
     *
     * Evaluate the metric at one location in the PatchData.
     *\param pd     The patch.
     *\param handle The location in the patch (as passed back from get_evaluations).
     *\param value  The output metric value.
     */
    MESQUITE_EXPORT virtual bool evaluate( PatchData& pd, size_t handle, double& value, MsqError& err ) = 0;

    /**\brief Get metric value at a logical location in the patch.
     *
     * Evaluate the metric at one location in the PatchData.
     *\param pd      The patch.
     *\param handle  The location in the patch (as passed back from get_evaluations).
     *\param value   The output metric value.
     *\param indices The free vertices that the evaluation is a function
     *               of, specified as vertex indices in the PatchData.
     */
    MESQUITE_EXPORT virtual bool evaluate_with_indices( PatchData& pd,
                                                        size_t handle,
                                                        double& value,
                                                        std::vector< size_t >& indices,
                                                        MsqError& err ) = 0;

    /**\brief Get metric value and gradient at a logical location in the patch.
     *
     * Evaluate the metric at one location in the PatchData.
     *\param pd      The patch.
     *\param handle  The location in the patch (as passed back from get_evaluations).
     *\param value   The output metric value.
     *\param indices The free vertices that the evaluation is a function
     *               of, specified as vertex indices in the PatchData.
     *\param gradient The gradient of the metric as a function of the
     *               coordinates of the free vertices passed back in
     *               the indices list.
     */
    MESQUITE_EXPORT virtual bool evaluate_with_gradient( PatchData& pd,
                                                         size_t handle,
                                                         double& value,
                                                         std::vector< size_t >& indices,
                                                         std::vector< Vector3D >& gradient,
                                                         MsqError& err );

    /**\brief Get metric value and gradient at a logical location in the patch.
     *
     * Evaluate the metric at one location in the PatchData.
     *\param pd      The patch.
     *\param handle  The location in the patch (as passed back from get_evaluations).
     *\param value   The output metric value.
     *\param indices The free vertices that the evaluation is a function
     *               of, specified as vertex indices in the PatchData.
     *\param gradient The gradient of the metric as a function of the
     *               coordinates of the free vertices passed back in
     *               the indices list.
     *\param Hessian_diagonal The 3x3 blocks along the diagonal of
     *               the Hessian matrix.
     */
    MESQUITE_EXPORT virtual bool evaluate_with_Hessian_diagonal( PatchData& pd,
                                                                 size_t handle,
                                                                 double& value,
                                                                 std::vector< size_t >& indices,
                                                                 std::vector< Vector3D >& gradient,
                                                                 std::vector< SymMatrix3D >& Hessian_diagonal,
                                                                 MsqError& err );

    /**\brief Get metric value and deravitives at a logical location in the patch.
     *
     * Evaluate the metric at one location in the PatchData.
     *\param pd      The patch.
     *\param handle  The location in the patch (as passed back from get_evaluations).
     *\param value   The output metric value.
     *\param indices The free vertices that the evaluation is a function
     *               of, specified as vertex indices in the PatchData.
     *\param gradient The gradient of the metric as a function of the
     *               coordinates of the free vertices passed back in
     *               the indices list.
     *\param Hessian The Hessian of the metric as a function of the
     *               coordinates. The Hessian is passed back as the
     *               upper-triangular portion of the matrix in row-major
     *               order, where each Matrix3D is the portion of the
     *               Hessian with respect to the vertices at the
     *               corresponding positions in the indices list.
     */
    MESQUITE_EXPORT virtual bool evaluate_with_Hessian( PatchData& pd,
                                                        size_t handle,
                                                        double& value,
                                                        std::vector< size_t >& indices,
                                                        std::vector< Vector3D >& gradient,
                                                        std::vector< Matrix3D >& Hessian,
                                                        MsqError& err );

    //! Escobar Barrier Function for Shape and Other Metrics
    // det = signed determinant of Jacobian Matrix at a Vertex
    // delta = scaling parameter
    static inline double vertex_barrier_function( double det, double delta )
    {
        return 0.5 * ( det + sqrt( det * det + 4 * delta * delta ) );
    }
    // protected:

    /** \brief Remove from vector any gradient terms corresponding
     *         to a fixed vertex.
     *
     * Remove terms from vector that correspond to fixed vertices.
     *\param type            Element type
     *\param fixed_vertices  Bit flags, one per vertex, 1 if
     *                       vertex is fixed.
     *\param gradients       Array of gradients
     */
    MESQUITE_EXPORT
    static void remove_fixed_gradients( EntityTopology type,
                                        uint32_t fixed_vertices,
                                        std::vector< Vector3D >& gradients );

    /** \brief Remove from vectors any gradient terms and hessian
     *         diagonal blcoks corresponding to a fixed vertex.
     *
     * Remove terms from vector that correspond to fixed vertices.
     *\param type            Element type
     *\param fixed_vertices  Bit flags, one per vertex, 1 if
     *                       vertex is fixed.
     *\param gradients       Array of gradients
     *\param hess_diagonal_blocks   Array of diagonal blocks of Hessian matrix.
     */
    MESQUITE_EXPORT
    static void remove_fixed_diagonals( EntityTopology type,
                                        uint32_t fixed_vertices,
                                        std::vector< Vector3D >& gradients,
                                        std::vector< SymMatrix3D >& hess_diagonal_blocks );

    /** \brief Remove from vector any Hessian blocks corresponding
     *         to a fixed vertex.
     *
     * Remove blocks from vector that correspond to fixed vertices.
     *\param type            Element type
     *\param fixed_vertices  Bit flags, one per vertex, 1 if
     *                       vertex is fixed.
     *\param hessians        Array of Hessian blocks (upper trianguler, row-major)
     */
    MESQUITE_EXPORT
    static void remove_fixed_hessians( EntityTopology type,
                                       uint32_t fixed_vertices,
                                       std::vector< Matrix3D >& hessians );

    /** \brief Convert fixed vertex format from list to bit flags
     *
     * Given list of pointers to fixed vertices as passed to
     * evaluation functions, convert to bit flag format used
     * for many utility functions in this class.  Bits correspond
     * to vertices in the canonical vertex order, beginning with
     * the least-significant bit.  The bit is cleared for free
     * vertices and set (1) for fixed vertices.
     */
    MESQUITE_EXPORT
    static uint32_t fixed_vertex_bitmap( PatchData& pd,
                                         const MsqMeshEntity* elem,
                                         std::vector< size_t >& free_indices );

    //! takes an array of coefficients and an array of metrics (both of length num_value)
    //! and averages the contents using averaging method 'method'.
    MESQUITE_EXPORT
    double weighted_average_metrics( const double coef[],
                                     const double metric_values[],
                                     const int& num_values,
                                     MsqError& err );

    /*!AveragingMethod allows you to set how the quality metric values
      attained at each sample point will be averaged together to produce
      a single metric value for an element.
    */
    enum AveragingMethod
    {
        LINEAR,       //!< the linear average
        RMS,          //!< the root-mean-squared average
        HMS,          //!< the harmonic-mean-squared average
        SUM,          //!< the sum of the values
        SUM_SQUARED,  //!< the sum of the squares of the values
        HARMONIC,     //!< the harmonic average
        LAST_WITH_HESSIAN = HARMONIC,
        MINIMUM,    //!< the minimum value
        MAXIMUM,    //!< the maximum value
        GEOMETRIC,  //!< the geometric average
        LAST_WITH_GRADIENT = GEOMETRIC,
        STANDARD_DEVIATION,    //!< the standard deviation squared of the values
        MAX_OVER_MIN,          //!< the maximum value minus the minum value
        MAX_MINUS_MIN,         //!< the maximum value divided by the minimum value
        SUM_OF_RATIOS_SQUARED  //!< (1/(N^2))*(SUM (SUM (v_i/v_j)^2))
    };

    //!\brief Called at start of instruction queue processing
    //!
    //! Do any preliminary global initialization, consistency checking,
    //! etc.  Default implementation does nothing.
    MESQUITE_EXPORT virtual void initialize_queue( MeshDomainAssoc* mesh_and_domain,
                                                   const Settings* settings,
                                                   MsqError& err );

  private:
    int feasible;

    std::vector< Matrix3D > tmpHess;
    bool keepFiniteDiffEps;  //!< True if gradient finite difference
                             //!< calculation should set \c finiteDiffEps
    bool haveFiniteDiffEps;  //!< True if finite difference Hessian code
                             //!< has calculated \c finiteDiffEps
    double finiteDiffEps;    //!< Location for finite difference Hessian code
                             //!< to store this value so that it doesn't need
                             //!< to be recalculated if the gradient calculation
                             //!< is also finite difference
};

}  // namespace MBMesquite

#endif  // QualityMetric_hpp