Cppcheck report - MOAB: test/mesquite/unit/TMPQualityMetricTest.hpp

/* *****************************************************************
    MESQUITE -- The Mesh Quality Improvement Toolkit

    Copyright 2006 Sandia National Laboratories.  Developed at the
    University of Wisconsin--Madison under SNL contract number
    624796.  The U.S. Government and the University of Wisconsin
    retian certain rights to 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

    (2006) [email protected]
    (2010) [email protected]

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

/*! \file TMPQualityMetricTest.cpp

Unit testing for the TMPQualityMetric class
\author Jasno Kraftcheck
*/

#include "IdealShapeTarget.hpp"
#include "MsqMatrix.hpp"
#include "QualityMetricTester.hpp"
#include "Settings.hpp"
#include "UnitUtil.hpp"
#include "PlanarDomain.hpp"
#include "PatchData.hpp"
#include "WeightCalculator.hpp"
#include "ElementPMeanP.hpp"
#include "ElemSampleQM.hpp"

#include <iostream>

using namespace MBMesquite;
using std::cerr;
using std::cout;
using std::endl;

/** Target metric (templatized by dimension) for use in misc. tests.
 *  'evaluate' method records input values and returns a constant.
 */
template < typename B >
class FauxMetric : public B
{
  public:
    int count;
    double value;
    bool rval;
    MsqMatrix< 2, 2 > last_A_2D;
    MsqMatrix< 3, 3 > last_A_3D;

    FauxMetric( double v ) : count( 0 ), value( v ), rval( true ) {}

    std::string get_name() const
    {
        return "Faux";
    }

    bool evaluate( const MsqMatrix< 2, 2 >& A, const MsqMatrix< 2, 2 >& W, double& result, MsqError& )
    {
        last_A_2D = A;
        result    = value;
        ++count;
        return rval;
    }

    bool evaluate( const MsqMatrix< 3, 3 >& A, const MsqMatrix< 3, 3 >& W, double& result, MsqError& )
    {
        last_A_3D = A;
        result    = value;
        ++count;
        return rval;
    }

    bool evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& )
    {
        last_A_2D = T;
        result    = value;
        ++count;
        return rval;
    }

    bool evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& )
    {
        last_A_3D = T;
        result    = value;
        ++count;
        return rval;
    }
};

/** Weight calculator used for testing.  Returns constant weight. */
class ScaleWeight : public WeightCalculator
{
  public:
    ScaleWeight( double s ) : value( s ) {}
<--- Class 'ScaleWeight' has a constructor with 1 argument that is not explicit. [+]
Class 'ScaleWeight' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.
<--- Class 'ScaleWeight' has a constructor with 1 argument that is not explicit. [+]
Class 'ScaleWeight' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.
    double get_weight( PatchData&, size_t, Sample, MsqError& )
    {
        return value;
    }
    double value;
};

/** wrapper class to force numeric approximation of derivatives */
template < class Base >
class NumericalMetric : public Base
{
  public:
    NumericalMetric( Base* real_metric ) : mMetric( real_metric ) {}

    ~NumericalMetric() {}

    std::string get_name() const
    {
        return "Numerical " + mMetric->get_name();
    }

    bool evaluate( const MsqMatrix< 2, 2 >& A, const MsqMatrix< 2, 2 >& W, double& result, MsqError& err )
    {
        return mMetric->evaluate( A, W, result, err );
    }

    bool evaluate( const MsqMatrix< 3, 3 >& A, const MsqMatrix< 3, 3 >& W, double& result, MsqError& err )
    {
        return mMetric->evaluate( A, W, result, err );
    }

    bool evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& err )
    {
        return mMetric->evaluate( T, result, err );
    }

    bool evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& err )
    {
        return mMetric->evaluate( T, result, err );
    }

  private:
    Base* mMetric;
};

/** Simple target metric for testing first partial derivatives.
 *  \f$\mu(A,W) = |A|^2\f$
 *  \f$\frac{\partial\mu}{\partial \A} = 2 A \f$
 */
template < class Base >
class TestGradTargetMetric : public Base
{
  public:
    std::string get_name() const
    {
        return "TestGrad";
    }

    bool evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& err )
    {
        result = sqr_Frobenius( T );
        return true;
    }

    bool evaluate( const MsqMatrix< 2, 2 >& A, const MsqMatrix< 2, 2 >&, double& result, MsqError& err )
    {
        return evaluate( A, result, err );
    }

    bool evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& d, MsqError& err )
    {
        result = sqr_Frobenius( T );
        d      = 2 * T;
        return true;
    }

    bool evaluate_with_grad( const MsqMatrix< 2, 2 >& A,
                             const MsqMatrix< 2, 2 >&,
                             double& result,
                             MsqMatrix< 2, 2 >& d,
                             MsqError& err )
    {
        return evaluate_with_grad( A, result, d, err );
    }

    bool evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& err )
    {
        result = sqr_Frobenius( T );
        return true;
    }

    bool evaluate( const MsqMatrix< 3, 3 >& A, const MsqMatrix< 3, 3 >&, double& result, MsqError& err )
    {
        return evaluate( A, result, err );
    }

    bool evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& d, MsqError& err )
    {
        result = sqr_Frobenius( T );
        d      = 2 * T;
        return true;
    }

    bool evaluate_with_grad( const MsqMatrix< 3, 3 >& A,
                             const MsqMatrix< 3, 3 >&,
                             double& result,
                             MsqMatrix< 3, 3 >& d,
                             MsqError& err )
    {
        return evaluate_with_grad( A, result, d, err );
    }
};

/* class to force evaluation of mapping function only at element center
 * so that we can re-use tests from QualityMetricTester that work only
 * for element-based metrics (TMP metric is "element based" if only one
 * sample in each element.)
 */
class CenterMF2D : public MappingFunction2D
{
  public:
    CenterMF2D( const MappingFunction2D* real_mf ) : myFunc( real_mf ){};
<--- Class 'CenterMF2D' has a constructor with 1 argument that is not explicit. [+]
Class 'CenterMF2D' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.
<--- Class 'CenterMF2D' has a constructor with 1 argument that is not explicit. [+]
Class 'CenterMF2D' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.
    EntityTopology element_topology() const
    {
        return myFunc->element_topology();
    }
    int num_nodes() const
    {
        return myFunc->num_nodes();
    }
    NodeSet sample_points( NodeSet ) const
    {
        NodeSet s;
        s.set_mid_face_node( 0 );
        return s;
    }
    void coefficients( Sample l, NodeSet s, double* c, size_t* i, size_t& n, MsqError& e ) const
    {
        myFunc->coefficients( l, s, c, i, n, e );
    }
    void derivatives( Sample l, NodeSet s, size_t* i, MsqVector< 2 >* c, size_t& n, MsqError& e ) const
    {
        myFunc->derivatives( l, s, i, c, n, e );
    }

  private:
    const MappingFunction2D* myFunc;
};
class CenterMF3D : public MappingFunction3D
{
  public:
    CenterMF3D( const MappingFunction3D* real_mf ) : myFunc( real_mf ){};
<--- Class 'CenterMF3D' has a constructor with 1 argument that is not explicit. [+]
Class 'CenterMF3D' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.
<--- Class 'CenterMF3D' has a constructor with 1 argument that is not explicit. [+]
Class 'CenterMF3D' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.
    EntityTopology element_topology() const
    {
        return myFunc->element_topology();
    }
    int num_nodes() const
    {
        return myFunc->num_nodes();
    }
    NodeSet sample_points( NodeSet ) const
    {
        NodeSet s;
        s.set_mid_region_node();
        return s;
    }
    void coefficients( Sample l, NodeSet s, double* c, size_t* i, size_t& n, MsqError& e ) const
    {
        myFunc->coefficients( l, s, c, i, n, e );
    }
    void derivatives( Sample l, NodeSet s, size_t* i, MsqVector< 3 >* c, size_t& n, MsqError& e ) const
    {
        myFunc->derivatives( l, s, i, c, n, e );
    }

  private:
    const MappingFunction3D* myFunc;
};

// Define a target calculator that returns targets for
// ideally shaped elements, but also includes orientation
// information (aligning surface elements to the xy plane
// with the first column of the jacobian in the x direction).
class IdealShapeXY : public IdealShapeTarget
{
  public:
    bool have_surface_orient() const
    {
        return true;
    }
    bool get_surface_target( PatchData& pd, size_t element, Sample sample, MsqMatrix< 3, 2 >& W_out, MsqError& err )
    {
        MsqMatrix< 2, 2 > W;
        bool rval = get_2D_target( pd, element, sample, W, err );
        W_out.set_row( 0, W.row( 0 ) );
        W_out.set_row( 1, W.row( 1 ) );
        W_out.set_row( 2, MsqMatrix< 1, 2 >( 0.0 ) );
        return rval;
    }
};

#define REGISTER_TMP_TESTS                                      \
                                                                \
    CPPUNIT_TEST( test_negate_flag );                           \
    CPPUNIT_TEST( test_supported_types );                       \
    CPPUNIT_TEST( test_get_evaluations );                       \
    CPPUNIT_TEST( test_get_element_evaluations );               \
                                                                \
    CPPUNIT_TEST( test_evaluate_2D );                           \
    CPPUNIT_TEST( test_evaluate_surface );                      \
    CPPUNIT_TEST( test_evaluate_3D );                           \
    CPPUNIT_TEST( test_evaluate_2D_weight );                    \
    CPPUNIT_TEST( test_evaluate_surface_weight );               \
    CPPUNIT_TEST( test_evaluate_3D_weight );                    \
    CPPUNIT_TEST( test_2d_eval_ortho_quad );                    \
    CPPUNIT_TEST( test_surf_eval_ortho_quad );                  \
    CPPUNIT_TEST( test_3d_eval_ortho_hex );                     \
                                                                \
    CPPUNIT_TEST( test_sample_indices );                        \
    CPPUNIT_TEST( test_evaluate_with_indices );                 \
    CPPUNIT_TEST( test_evaluate_fixed_indices );                \
                                                                \
    CPPUNIT_TEST( test_gradient_2D );                           \
    CPPUNIT_TEST( test_gradient_surface );                      \
    CPPUNIT_TEST( test_gradient_3D );                           \
    CPPUNIT_TEST( compare_indices_and_gradient );               \
    CPPUNIT_TEST( test_ideal_element_gradient );                \
    CPPUNIT_TEST( compare_analytical_and_numerical_gradient );  \
    CPPUNIT_TEST( test_weighted_gradients );                    \
    CPPUNIT_TEST( test_gradient_with_fixed_vertices );          \
                                                                \
    CPPUNIT_TEST( compare_indices_and_hessian );                \
    CPPUNIT_TEST( compare_gradient_and_hessian );               \
    CPPUNIT_TEST( compare_analytical_and_numerical_hessians );  \
    CPPUNIT_TEST( test_symmetric_hessian_diagonal );            \
    CPPUNIT_TEST( test_weighted_hessians );                     \
    CPPUNIT_TEST( test_hessian_with_fixed_vertices );           \
                                                                \
    CPPUNIT_TEST( compare_indices_and_diagonal );               \
    CPPUNIT_TEST( compare_gradient_and_diagonal );              \
    CPPUNIT_TEST( compare_analytical_and_numerical_diagonals ); \
    CPPUNIT_TEST( test_weighted_diagonals );                    \
    CPPUNIT_TEST( test_diagonal_with_fixed_vertices );

static double col_dot_prod( MsqMatrix< 2, 2 >& m )<--- Parameter 'm' can be declared with const<--- Parameter 'm' can be declared with const
{
    return m( 0, 0 ) * m( 0, 1 ) + m( 1, 0 ) * m( 1, 1 );
}

template < class QMType >
class TMPTypes
{
};

template < class QMType >
class TMPQualityMetricTest : public CppUnit::TestFixture
{
  protected:
    QualityMetricTester tester;

    Settings settings;
    IdealShapeTarget ideal;
    IdealShapeXY surf_target;
    ScaleWeight e_weight;

    FauxMetric< typename TMPTypes< QMType >::MetricType > faux_pi, faux_zero, faux_two;
    typename TMPTypes< QMType >::TestType test_metric;
    NumericalMetric< typename QMType::MetricType > num_metric;
    QMType test_qm, test_qm_surf, zero_qm, weight_qm, center_qm;
    Settings centerOnly;
    CenterMF2D triCenter, quadCenter;
    CenterMF3D tetCenter, pyrCenter, priCenter, hexCenter;

  public:
    TMPQualityMetricTest()
        : tester( QualityMetricTester::ALL_FE_EXCEPT_SEPTAHEDRON, &settings ), e_weight( 2.7182818284590451 ),
          faux_pi( 3.14159 ), faux_zero( 0.0 ), faux_two( 2.0 ), num_metric( &test_metric ),
          test_qm( &ideal, &num_metric ), test_qm_surf( &surf_target, &num_metric ), zero_qm( &ideal, &faux_zero ),
          weight_qm( &ideal, &e_weight, &test_metric ), center_qm( &ideal, &test_metric ),
          triCenter( centerOnly.get_mapping_function_2D( TRIANGLE ) ),
          quadCenter( centerOnly.get_mapping_function_2D( QUADRILATERAL ) ),
          tetCenter( centerOnly.get_mapping_function_3D( TETRAHEDRON ) ),
          pyrCenter( centerOnly.get_mapping_function_3D( PYRAMID ) ),
          priCenter( centerOnly.get_mapping_function_3D( PRISM ) ),
          hexCenter( centerOnly.get_mapping_function_3D( HEXAHEDRON ) )
    {
        centerOnly.set_mapping_function( &triCenter );
        centerOnly.set_mapping_function( &quadCenter );
        centerOnly.set_mapping_function( &tetCenter );
        centerOnly.set_mapping_function( &pyrCenter );
        centerOnly.set_mapping_function( &priCenter );
        centerOnly.set_mapping_function( &hexCenter );
        tester.ideal_pyramid_base_equals_height( true );
    }

    void test_negate_flag()
    {
        CPPUNIT_ASSERT_EQUAL( 1, zero_qm.get_negate_flag() );
    }
    void test_supported_types()
    {
        tester.test_supported_element_types( &zero_qm );
    }
    void test_get_evaluations()
    {
        QMType edge_metric( &ideal, &faux_zero );
        tester.test_get_sample_evaluations( &zero_qm );
        tester.test_get_sample_evaluations( &edge_metric );
    }
    void test_get_element_evaluations()
    {
        QMType edge_metric( &ideal, &faux_zero );
        tester.test_get_in_element_evaluations( &zero_qm );
        tester.test_get_in_element_evaluations( &edge_metric );
    }

    void test_evaluate_2D();
    void test_evaluate_surface();
    void test_evaluate_3D();

    void test_evaluate_2D_weight()
    {
        MsqPrintError err( cout );
        PatchData pd;
        bool rval;
        double value;

        QMType m( &ideal, &e_weight, &faux_pi );
        tester.get_ideal_element( TRIANGLE, true, pd );
        rval = m.evaluate( pd, 0, value, err );
        CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
        CPPUNIT_ASSERT( rval );
        CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value * e_weight.value, value, DBL_EPSILON );
    }

    void test_evaluate_surface_weight()
    {
        MsqPrintError err( cout );
        PatchData pd;
        bool rval;
        double value;

        QMType m( &surf_target, &e_weight, &faux_pi );

        tester.get_ideal_element( TRIANGLE, true, pd );
        rval = m.evaluate( pd, 0, value, err );
        CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
        CPPUNIT_ASSERT( rval );
        CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value * e_weight.value, value, DBL_EPSILON );
    }

    void test_evaluate_3D_weight()
    {
        MsqPrintError err( cout );
        PatchData pd;
        bool rval;
        double value;

        QMType m( &ideal, &e_weight, &faux_two );

        tester.get_ideal_element( PRISM, true, pd );
        rval = m.evaluate( pd, 0, value, err );
        CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
        CPPUNIT_ASSERT( rval );
        CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_two.value * e_weight.value, value, DBL_EPSILON );
    }

    void test_2d_eval_ortho_quad()
    {
        MsqPrintError err( cout );
        PatchData pd;
        bool rval;
        double value;

        QMType m( &ideal, &faux_zero );
        faux_zero.count = 0;

        tester.get_ideal_element( QUADRILATERAL, true, pd );
        rval = m.evaluate( pd, 0, value, err );
        CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
        CPPUNIT_ASSERT( rval );
        CPPUNIT_ASSERT_EQUAL( 1, faux_zero.count );
        CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, col_dot_prod( faux_zero.last_A_2D ), DBL_EPSILON );
    }

    void test_surf_eval_ortho_quad()
    {
        MsqPrintError err( cout );
        PatchData pd;
        bool rval;
        double value;

        QMType m( &surf_target, &faux_zero );
        faux_zero.count = 0;

        tester.get_ideal_element( QUADRILATERAL, true, pd );
        rval = m.evaluate( pd, 0, value, err );
        CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
        CPPUNIT_ASSERT( rval );
        CPPUNIT_ASSERT_EQUAL( 1, faux_zero.count );
        CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, col_dot_prod( faux_zero.last_A_2D ), DBL_EPSILON );
    }

    void test_3d_eval_ortho_hex()
    {
        MsqPrintError err( cout );
        PatchData pd;
        bool rval;
        double value;

        QMType m( &ideal, &faux_zero );
        faux_zero.count = 0;

        tester.get_ideal_element( HEXAHEDRON, true, pd );
        rval = m.evaluate( pd, 0, value, err );
        CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
        CPPUNIT_ASSERT( rval );
        CPPUNIT_ASSERT_EQUAL( 1, faux_zero.count );

        // test that columns are orthogonal for ideal hex element
        MsqMatrix< 3, 3 > A = faux_zero.last_A_3D;
        CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, A.column( 0 ) % A.column( 1 ), 1e-6 );
        CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, A.column( 0 ) % A.column( 2 ), 1e-6 );
        CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, A.column( 1 ) % A.column( 2 ), 1e-6 );
    }

    void test_gradient_common( TargetCalculator* tc );
    void test_gradient_2D()
    {
        test_gradient_common( &ideal );
    }
    void test_gradient_surface()
    {
        test_gradient_common( &surf_target );
    }
    void test_gradient_3D();

    void test_sample_indices()
    {
        tester.test_get_sample_indices( &zero_qm );
    }
    void test_evaluate_with_indices()
    {
        tester.compare_eval_and_eval_with_indices( &zero_qm );
    }
    void test_evaluate_fixed_indices()
    {
        tester.test_get_indices_fixed( &zero_qm );
    }

    void compare_indices_and_gradient()
    {
        tester.compare_eval_with_indices_and_eval_with_gradient( &test_qm );
        tester.compare_eval_with_indices_and_eval_with_gradient( &test_qm_surf );
    }
    void test_ideal_element_gradient()
    {
        tester.test_ideal_element_zero_gradient( &test_qm, false );
        tester.test_ideal_element_zero_gradient( &test_qm_surf, false );
    }
    void compare_analytical_and_numerical_gradient()
    {
        compare_analytical_and_numerical_gradients( &test_qm );
        compare_analytical_and_numerical_gradients( &test_qm_surf );
    }
    void test_weighted_gradients()
    {
        compare_analytical_and_numerical_gradients( &weight_qm );
    }
    void test_gradient_with_fixed_vertices()
    {
        tester.test_gradient_with_fixed_vertex( &center_qm, &centerOnly );
    }

    void compare_indices_and_hessian()
    {
        tester.compare_eval_with_indices_and_eval_with_hessian( &test_qm );
        tester.compare_eval_with_indices_and_eval_with_hessian( &test_qm_surf );
    }
    void compare_gradient_and_hessian()
    {
        tester.compare_eval_with_grad_and_eval_with_hessian( &test_qm );
        tester.compare_eval_with_grad_and_eval_with_hessian( &test_qm_surf );
    }
    void compare_analytical_and_numerical_hessians()
    {
        compare_analytical_and_numerical_hessians( &test_qm );
        compare_analytical_and_numerical_hessians( &test_qm_surf );
    }
    void test_symmetric_hessian_diagonal()
    {
        tester.test_symmetric_Hessian_diagonal_blocks( &test_qm );
        tester.test_symmetric_Hessian_diagonal_blocks( &test_qm_surf );
    }
    void test_weighted_hessians()
    {
        compare_analytical_and_numerical_hessians( &weight_qm );
    }
    void test_hessian_with_fixed_vertices()
    {
        tester.test_hessian_with_fixed_vertex( &center_qm, &centerOnly );
    }

    void compare_indices_and_diagonal()
    {
        tester.compare_eval_with_indices_and_eval_with_diagonal( &test_qm );
        tester.compare_eval_with_indices_and_eval_with_diagonal( &test_qm_surf );
    }
    void compare_gradient_and_diagonal()
    {
        tester.compare_eval_with_grad_and_eval_with_diagonal( &test_qm );
        tester.compare_eval_with_grad_and_eval_with_diagonal( &test_qm_surf );
    }
    void compare_analytical_and_numerical_diagonals()
    {
        compare_analytical_and_numerical_diagonals( &test_qm );
        compare_analytical_and_numerical_diagonals( &test_qm_surf );
    }
    void test_weighted_diagonals()
    {
        compare_analytical_and_numerical_diagonals( &weight_qm );
    }
    void test_diagonal_with_fixed_vertices()
    {
        tester.test_diagonal_with_fixed_vertex( &center_qm, &centerOnly );
    }

    // Delcare specialized versions of the functions from
    // QualityMetricTester because we surface elements must
    // be handled differently.  For a surface element in the XY plane,
    // the finite difference approximations of the derivatives will
    // have non-zero values for derivatives wrt Z coordinates while the
    // analytical derivative calculations will return all derivatives
    // wrt Z coordiantes as zero.

    void get_nonideal_element( EntityTopology type, PatchData& pd )
    {
        tester.get_nonideal_element( type, pd, true );
        // Callers assume surface elements are in XY plane.
        // Verify this assumption.
        if( TopologyInfo::dimension( type ) == 2 )
        {
            for( size_t i = 0; i < pd.num_nodes(); ++i )
            {
                CPPUNIT_ASSERT_DOUBLES_EQUAL( pd.vertex_by_index( i )[2], 0.0, 1e-6 );
            }
        }
    }

    void compare_analytical_and_numerical_gradients( QualityMetric* qm )
    {
        PatchData pd;
        const EntityTopology types[] = { TRIANGLE, QUADRILATERAL, TETRAHEDRON, PYRAMID, PRISM, HEXAHEDRON };
        const int num_types          = sizeof( types ) / sizeof( types[0] );
        for( int i = 0; i < num_types; ++i )
        {
            get_nonideal_element( types[i], pd );
            compare_analytical_and_numerical_gradients( qm, pd, TopologyInfo::dimension( types[i] ) );
        }
    }

    void compare_analytical_and_numerical_hessians( QualityMetric* qm );
    void compare_analytical_and_numerical_diagonals( QualityMetric* qm );
    void compare_analytical_and_numerical_gradients( QualityMetric* qm, PatchData&, int dim );
};

// CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(TMPQualityMetricTest, "TMPQualityMetricTest");
// CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(TMPQualityMetricTest, "Unit");

template < class QMType >
inline void TMPQualityMetricTest< QMType >::test_evaluate_2D()
{
    MsqPrintError err( cout );
    PatchData pd;
    bool rval;
    double value;

    QMType m( &ideal, &faux_pi );

    // test with aligned elements
    faux_pi.count = 0;
    tester.get_ideal_element( QUADRILATERAL, true, pd );
    rval = m.evaluate( pd, 0, value, err );
    CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
    CPPUNIT_ASSERT( rval );
    CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
    CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );

    // test that columns are orthogonal for ideal quad element
    CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, col_dot_prod( faux_pi.last_A_2D ), 1e-6 );

    // test with an element rotated about X-axis
    faux_pi.count = 0;
    tester.get_ideal_element( QUADRILATERAL, true, pd );
    // rotate by 90 degrees about X axis
    for( size_t i = 0; i < pd.num_nodes(); ++i )
    {
        Vector3D orig = pd.vertex_by_index( i );
        Vector3D newp( orig[0], -orig[2], orig[1] );
        pd.set_vertex_coordinates( newp, i, err );
    }
    rval = m.evaluate( pd, 0, value, err );
    CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
    CPPUNIT_ASSERT( rval );
    CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
    CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );

    // test with an element rotated about Y-axis
    faux_pi.count = 0;
    tester.get_ideal_element( TRIANGLE, true, pd );
    // rotate by -90 degrees about Y axis
    for( size_t i = 0; i < pd.num_nodes(); ++i )
    {
        Vector3D orig = pd.vertex_by_index( i );
        Vector3D newp( orig[2], orig[1], -orig[0] );
        pd.set_vertex_coordinates( newp, i, err );
    }
    rval = m.evaluate( pd, 0, value, err );
    CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
    CPPUNIT_ASSERT( rval );
    CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
    CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );
}

template < class QMType >
inline void TMPQualityMetricTest< QMType >::test_evaluate_surface()
{
    MsqPrintError err( cout );
    PatchData pd;
    bool rval;
    double value;

    QMType m( &surf_target, &faux_pi );

    // test with aligned elements
    faux_pi.count = 0;
    tester.get_ideal_element( QUADRILATERAL, true, pd );
    rval = m.evaluate( pd, 0, value, err );
    CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
    CPPUNIT_ASSERT( rval );
    CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
    CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );

    // test that columns are orthogonal for ideal quad element
    CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, col_dot_prod( faux_pi.last_A_2D ), 1e-6 );

    // test with an element rotated about X-axis
    faux_pi.count = 0;
    tester.get_ideal_element( QUADRILATERAL, true, pd );
    // rotate by 90 degrees about X axis
    for( size_t i = 0; i < pd.num_nodes(); ++i )
    {
        Vector3D orig = pd.vertex_by_index( i );
        Vector3D newp( orig[0], -orig[2], orig[1] );
        pd.set_vertex_coordinates( newp, i, err );
    }
    rval = m.evaluate( pd, 0, value, err );
    CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
    CPPUNIT_ASSERT( rval );
    CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
    CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );

    // test with an element rotated about Y-axis
    faux_pi.count = 0;
    tester.get_ideal_element( TRIANGLE, true, pd );
    // rotate by -90 degrees about Y axis
    for( size_t i = 0; i < pd.num_nodes(); ++i )
    {
        Vector3D orig = pd.vertex_by_index( i );
        Vector3D newp( orig[2], orig[1], -orig[0] );
        pd.set_vertex_coordinates( newp, i, err );
    }
    rval = m.evaluate( pd, 0, value, err );
    CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
    CPPUNIT_ASSERT( rval );
    CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
    CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );
}

template < class QMType >
inline void TMPQualityMetricTest< QMType >::test_evaluate_3D()
{
    MsqPrintError err( cout );
    PatchData pd;
    bool rval;
    double value;

    QMType m( &ideal, &faux_two );

    // test with aligned elements
    faux_two.count = 0;
    tester.get_ideal_element( HEXAHEDRON, true, pd );
    rval = m.evaluate( pd, 0, value, err );
    CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
    CPPUNIT_ASSERT( rval );
    CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_two.value, value, DBL_EPSILON );
    CPPUNIT_ASSERT_EQUAL( 1, faux_two.count );

    // test that columns are orthogonal for ideal hex element
    MsqMatrix< 3, 3 > A = faux_two.last_A_3D;
    CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, A.column( 0 ) % A.column( 1 ), 1e-6 );
    CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, A.column( 0 ) % A.column( 2 ), 1e-6 );
    CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, A.column( 1 ) % A.column( 2 ), 1e-6 );

    // test with rotated element
    faux_two.count = 0;
    tester.get_ideal_element( TETRAHEDRON, true, pd );
    // rotate by 90-degrees about X axis
    for( size_t i = 0; i < pd.num_nodes(); ++i )
    {
        Vector3D orig = pd.vertex_by_index( i );
        Vector3D newp( orig[0], -orig[2], orig[1] );
        pd.set_vertex_coordinates( newp, i, err );
    }
    rval = m.evaluate( pd, 0, value, err );
    CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
    CPPUNIT_ASSERT( rval );
    CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_two.value, value, DBL_EPSILON );
    CPPUNIT_ASSERT_EQUAL( 1, faux_two.count );
}

template < class QMType >
inline void TMPQualityMetricTest< QMType >::test_gradient_common( TargetCalculator* tc )
{
    MsqPrintError err( std::cout );

    // check for expected value at center of flattened hex

    // construct flattened hex
    const double y          = 0.5;
    const double vertices[] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, y, 0.0, 0.0, y, 0.0 };
    size_t conn[8]          = { 0, 1, 2, 3 };
    PatchData pd;
    pd.fill( 4, vertices, 1, QUADRILATERAL, conn, 0, err );
    ASSERT_NO_ERROR( err );

    // calculate Jacobian matrix at element center

    // derivatives of bilinear map at quad center
    const double deriv_vals[] = { -0.5, -0.5, 0.5, -0.5, 0.5, 0.5, -0.5, 0.5 };
    MsqMatrix< 4, 2 > coeff_derivs( deriv_vals );
    MsqMatrix< 4, 3 > coords( vertices );
    MsqMatrix< 3, 2 > J = transpose( coords ) * coeff_derivs;
    // calculate expected metric value
    const double expt_val = sqr_Frobenius( J );
    // calculate derivative for each element vertex
    MsqVector< 3 > expt_grad[4];
    for( int v = 0; v < 4; ++v )
        expt_grad[v] = 2 * J * transpose( coeff_derivs.row( v ) );

    // construct metric
    pd.attach_settings( &settings );
    TestGradTargetMetric< typename TMPTypes< QMType >::MetricType > tm;
    // IdealShapeTarget tc;
    QMType m( tc, &tm );
    PlanarDomain plane( PlanarDomain::XY );
    pd.set_domain( &plane );

    // evaluate metric
    double act_val;
    std::vector< size_t > indices;
    std::vector< Vector3D > act_grad;
    size_t h = ElemSampleQM::handle( Sample( 2, 0 ), 0 );
    m.evaluate_with_gradient( pd, h, act_val, indices, act_grad, err );
    ASSERT_NO_ERROR( err );

    // compare values
    CPPUNIT_ASSERT_DOUBLES_EQUAL( expt_val, act_val, 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[0]].data() ), act_grad[0], 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[1]].data() ), act_grad[1], 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[2]].data() ), act_grad[2], 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[3]].data() ), act_grad[3], 1e-10 );

    // check numerical approx of gradient
    m.QualityMetric::evaluate_with_gradient( pd, h, act_val, indices, act_grad, err );
    ASSERT_NO_ERROR( err );
    CPPUNIT_ASSERT_DOUBLES_EQUAL( expt_val, act_val, 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[0]].data() ), act_grad[0], 1e-5 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[1]].data() ), act_grad[1], 1e-5 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[2]].data() ), act_grad[2], 1e-5 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[3]].data() ), act_grad[3], 1e-5 );
}

template < class QMType >
inline void TMPQualityMetricTest< QMType >::test_gradient_3D()
{
    MsqPrintError err( std::cout );

    // check for expected value at center of flattened hex

    // construct flattened hex
    const double z          = 0.5;
    const double vertices[] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0,
                                0.0, 0.0, z,   1.0, 0.0, z,   1.0, 1.0, z,   0.0, 1.0, z };
    size_t conn[8]          = { 0, 1, 2, 3, 4, 5, 6, 7 };
    PatchData pd;
    pd.fill( 8, vertices, 1, HEXAHEDRON, conn, 0, err );
    ASSERT_NO_ERROR( err );

    // calculate Jacobian matrix at element center

    // derivatives of trilinear map at hex center
    const double deriv_vals[8][3] = { { -0.25, -0.25, -0.25 }, { 0.25, -0.25, -0.25 }, { 0.25, 0.25, -0.25 },
                                      { -0.25, 0.25, -0.25 },  { -0.25, -0.25, 0.25 }, { 0.25, -0.25, 0.25 },
                                      { 0.25, 0.25, 0.25 },    { -0.25, 0.25, 0.25 } };
    MsqMatrix< 8, 3 > coeff_derivs( deriv_vals );
    MsqMatrix< 8, 3 > coords( vertices );
    MsqMatrix< 3, 3 > J = transpose( coords ) * coeff_derivs;
    // calculate expected metric value
    const double expt_val = sqr_Frobenius( J );
    // calculate derivative for each element vertex
    MsqVector< 3 > expt_grad[8];
    for( int v = 0; v < 8; ++v )
        expt_grad[v] = 2 * J * transpose( coeff_derivs.row( v ) );

    // construct metric
    pd.attach_settings( &settings );
    TestGradTargetMetric< typename TMPTypes< QMType >::MetricType > tm;
    IdealShapeTarget tc;
    QMType m( &tc, 0, &tm );

    // evaluate metric
    double act_val;
    std::vector< size_t > indices;
    std::vector< Vector3D > act_grad;
    size_t h = ElemSampleQM::handle( Sample( 3, 0 ), 0 );
    m.evaluate_with_gradient( pd, h, act_val, indices, act_grad, err );
    ASSERT_NO_ERROR( err );

    // compare values
    CPPUNIT_ASSERT_DOUBLES_EQUAL( expt_val, act_val, 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[0]].data() ), act_grad[0], 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[1]].data() ), act_grad[1], 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[2]].data() ), act_grad[2], 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[3]].data() ), act_grad[3], 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[4]].data() ), act_grad[4], 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[5]].data() ), act_grad[5], 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[6]].data() ), act_grad[6], 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[7]].data() ), act_grad[7], 1e-10 );

    // check numerical approx of gradient
    m.QualityMetric::evaluate_with_gradient( pd, h, act_val, indices, act_grad, err );
    ASSERT_NO_ERROR( err );
    CPPUNIT_ASSERT_DOUBLES_EQUAL( expt_val, act_val, 1e-10 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[0]].data() ), act_grad[0], 1e-5 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[1]].data() ), act_grad[1], 1e-5 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[2]].data() ), act_grad[2], 1e-5 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[3]].data() ), act_grad[3], 1e-5 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[4]].data() ), act_grad[4], 1e-5 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[5]].data() ), act_grad[5], 1e-5 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[6]].data() ), act_grad[6], 1e-5 );
    CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[7]].data() ), act_grad[7], 1e-5 );
}

template < class QMType >
inline void TMPQualityMetricTest< QMType >::compare_analytical_and_numerical_gradients( QualityMetric* qm,
                                                                                        PatchData& pd,
                                                                                        int dim )
{
    MsqPrintError err( std::cout );

    std::vector< size_t > handles, indices1, indices2;
    std::vector< Vector3D > grad1, grad2;
    double qm_val1, qm_val2;
    bool rval;
<--- The scope of the variable 'rval' can be reduced. [+]
The scope of the variable 'rval' can be reduced. Warning: Be careful when fixing this message, especially when there are inner loops. Here is an example where cppcheck will write that the scope for 'i' can be reduced:
void f(int x)
{
    int i = 0;
    if (x) {
        // it's safe to move 'int i = 0;' here
        for (int n = 0; n < 10; ++n) {
            // it is possible but not safe to move 'int i = 0;' here
            do_something(&i);
        }
    }
}
When you see this message it is always safe to reduce the variable scope 1 level.
<--- The scope of the variable 'rval' can be reduced. [+]
The scope of the variable 'rval' can be reduced. Warning: Be careful when fixing this message, especially when there are inner loops. Here is an example where cppcheck will write that the scope for 'i' can be reduced:
void f(int x)
{
    int i = 0;
    if (x) {
        // it's safe to move 'int i = 0;' here
        for (int n = 0; n < 10; ++n) {
            // it is possible but not safe to move 'int i = 0;' here
            do_something(&i);
        }
    }
}
When you see this message it is always safe to reduce the variable scope 1 level.

    qm->get_evaluations( pd, handles, false, err );
    CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
    CPPUNIT_ASSERT( !handles.empty() );
    for( size_t j = 0; j < handles.size(); ++j )
    {
        rval = qm->QualityMetric::evaluate_with_gradient( pd, handles[j], qm_val1, indices1, grad1, err );
        CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
        CPPUNIT_ASSERT( rval );

        // For analytical gradient of a 2D element in the XY plane,
        // we expect all Z terms to be zero.
        if( dim == 2 )
            for( size_t k = 0; k < grad1.size(); ++k )
                grad1[k][2] = 0.0;

        rval = qm->evaluate_with_gradient( pd, handles[j], qm_val2, indices2, grad2, err );
        CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
        CPPUNIT_ASSERT( rval );

        CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
        CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
        CPPUNIT_ASSERT( !indices1.empty() );

        std::vector< size_t >::iterator it1, it2;
        for( it1 = indices1.begin(); it1 != indices1.end(); ++it1 )
        {
            it2 = std::find( indices2.begin(), indices2.end(), *it1 );
            CPPUNIT_ASSERT( it2 != indices2.end() );

            size_t idx1 = it1 - indices1.begin();
            size_t idx2 = it2 - indices2.begin();
            CPPUNIT_ASSERT_VECTORS_EQUAL( grad1[idx1], grad2[idx2], 0.01 );
        }
    }
}

// Delcare specialized versions of the functions from
// QualityMetricTester because we surface elements must
// be handled differently.  For a surface element in the XY plane,
// the finite difference approximations of the derivatives will
// have non-zero values for derivatives wrt Z coordinates while the
// analytical derivative calculations will return all derivatives
// wrt Z coordiantes as zero.
template < class QMType >
inline void TMPQualityMetricTest< QMType >::compare_analytical_and_numerical_hessians( QualityMetric* qm )
{
    MsqPrintError err( std::cout );
    PatchData pd;
    const EntityTopology types[] = { TRIANGLE, QUADRILATERAL, TETRAHEDRON, PYRAMID, PRISM, HEXAHEDRON };
    const int num_types          = sizeof( types ) / sizeof( types[0] );
    for( int i = 0; i < num_types; ++i )
    {
        get_nonideal_element( types[i], pd );

        std::vector< size_t > handles, indices1, indices2;
        std::vector< Vector3D > grad1, grad2;
        std::vector< Matrix3D > Hess1, Hess2;
        double qm_val1, qm_val2;
        bool rval;
<--- The scope of the variable 'rval' can be reduced. [+]
The scope of the variable 'rval' can be reduced. Warning: Be careful when fixing this message, especially when there are inner loops. Here is an example where cppcheck will write that the scope for 'i' can be reduced:
void f(int x)
{
    int i = 0;
    if (x) {
        // it's safe to move 'int i = 0;' here
        for (int n = 0; n < 10; ++n) {
            // it is possible but not safe to move 'int i = 0;' here
            do_something(&i);
        }
    }
}
When you see this message it is always safe to reduce the variable scope 1 level.
<--- The scope of the variable 'rval' can be reduced. [+]
The scope of the variable 'rval' can be reduced. Warning: Be careful when fixing this message, especially when there are inner loops. Here is an example where cppcheck will write that the scope for 'i' can be reduced:
void f(int x)
{
    int i = 0;
    if (x) {
        // it's safe to move 'int i = 0;' here
        for (int n = 0; n < 10; ++n) {
            // it is possible but not safe to move 'int i = 0;' here
            do_something(&i);
        }
    }
}
When you see this message it is always safe to reduce the variable scope 1 level.

        qm->get_evaluations( pd, handles, false, err );
        CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
        CPPUNIT_ASSERT( !handles.empty() );
        for( size_t j = 0; j < handles.size(); ++j )
        {
            rval = qm->QualityMetric::evaluate_with_Hessian( pd, handles[j], qm_val1, indices1, grad1, Hess1, err );
            CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
            CPPUNIT_ASSERT( rval );

            // For analytical gradient of a 2D element in the XY plane,
            // we expect all Z terms to be zero.
#ifdef PLANAR_HESSIAN
            if( TopologyInfo::dimension( types[i] ) == 2 )
                for( size_t k = 0; k < Hess1.size(); ++k )
                    Hess1[k]( 0, 2 ) = Hess1[k]( 1, 2 ) = Hess1[k]( 2, 0 ) = Hess1[k]( 2, 1 ) = Hess1[k]( 2, 2 ) = 0.0;
#endif

            rval = qm->evaluate_with_Hessian( pd, handles[j], qm_val2, indices2, grad2, Hess2, err );
            CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
            CPPUNIT_ASSERT( rval );

            CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
            CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
            CPPUNIT_ASSERT( !indices1.empty() );

            std::vector< size_t >::iterator it;
            unsigned h = 0;
            for( unsigned r = 0; r < indices1.size(); ++r )
            {
                it = std::find( indices2.begin(), indices2.end(), indices1[r] );
                CPPUNIT_ASSERT( it != indices2.end() );
                unsigned r2 = it - indices2.begin();

                for( unsigned c = r; c < indices1.size(); ++c, ++h )
                {
                    it = std::find( indices2.begin(), indices2.end(), indices1[c] );
                    CPPUNIT_ASSERT( it != indices2.end() );
                    unsigned c2 = it - indices2.begin();

                    unsigned h2;
                    if( r2 <= c2 )
                        h2 = indices2.size() * r - r * ( r + 1 ) / 2 + c;
                    else
                        h2 = indices2.size() * c - c * ( c + 1 ) / 2 + r;

                    // if (!utest_mat_equal(Hess1[h],Hess2[h2],0.001))
                    //  assert(false);
                    CPPUNIT_ASSERT_MATRICES_EQUAL( Hess1[h], Hess2[h2], 0.05 );
                }
            }
        }
    }
}

// Delcare specialized versions of the functions from
// QualityMetricTester because we surface elements must
// be handled differently.  For a surface element in the XY plane,
// the finite difference approximations of the derivatives will
// have non-zero values for derivatives wrt Z coordinates while the
// analytical derivative calculations will return all derivatives
// wrt Z coordiantes as zero.
template < class QMType >
inline void TMPQualityMetricTest< QMType >::compare_analytical_and_numerical_diagonals( QualityMetric* qm )
{
    MsqPrintError err( std::cout );
    PatchData pd;
    const EntityTopology types[] = { TRIANGLE, QUADRILATERAL, TETRAHEDRON, PYRAMID, PRISM, HEXAHEDRON };
    const int num_types          = sizeof( types ) / sizeof( types[0] );
    for( int i = 0; i < num_types; ++i )
    {
        get_nonideal_element( types[i], pd );

        std::vector< size_t > handles, indices1, indices2;
        std::vector< Vector3D > grad1, grad2;
        std::vector< Matrix3D > Hess1;
        std::vector< SymMatrix3D > Hess2;
        double qm_val1, qm_val2;
        bool rval;
<--- The scope of the variable 'rval' can be reduced. [+]
The scope of the variable 'rval' can be reduced. Warning: Be careful when fixing this message, especially when there are inner loops. Here is an example where cppcheck will write that the scope for 'i' can be reduced:
void f(int x)
{
    int i = 0;
    if (x) {
        // it's safe to move 'int i = 0;' here
        for (int n = 0; n < 10; ++n) {
            // it is possible but not safe to move 'int i = 0;' here
            do_something(&i);
        }
    }
}
When you see this message it is always safe to reduce the variable scope 1 level.
<--- The scope of the variable 'rval' can be reduced. [+]
The scope of the variable 'rval' can be reduced. Warning: Be careful when fixing this message, especially when there are inner loops. Here is an example where cppcheck will write that the scope for 'i' can be reduced:
void f(int x)
{
    int i = 0;
    if (x) {
        // it's safe to move 'int i = 0;' here
        for (int n = 0; n < 10; ++n) {
            // it is possible but not safe to move 'int i = 0;' here
            do_something(&i);
        }
    }
}
When you see this message it is always safe to reduce the variable scope 1 level.

        qm->get_evaluations( pd, handles, false, err );
        CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
        CPPUNIT_ASSERT( !handles.empty() );
        for( size_t j = 0; j < handles.size(); ++j )
        {
            rval = qm->QualityMetric::evaluate_with_Hessian( pd, handles[j], qm_val1, indices1, grad1, Hess1, err );
            CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
            CPPUNIT_ASSERT( rval );

            // For analytical gradient of a 2D element in the XY plane,
            // we expect all Z terms to be zero.
#ifdef PLANAR_HESSIAN
            if( TopologyInfo::dimension( types[i] ) == 2 )
                for( size_t k = 0; k < Hess1.size(); ++k )
                    Hess1[k]( 0, 2 ) = Hess1[k]( 1, 2 ) = Hess1[k]( 2, 0 ) = Hess1[k]( 2, 1 ) = Hess1[k]( 2, 2 ) = 0.0;
#endif

            rval = qm->evaluate_with_Hessian_diagonal( pd, handles[j], qm_val2, indices2, grad2, Hess2, err );
            CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
            CPPUNIT_ASSERT( rval );

            CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
            CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
            CPPUNIT_ASSERT( !indices1.empty() );
            CPPUNIT_ASSERT_EQUAL( indices1.size() * ( indices1.size() + 1 ) / 2, Hess1.size() );
            CPPUNIT_ASSERT_EQUAL( indices2.size(), Hess2.size() );

            size_t h = 0;
            std::vector< size_t >::iterator it;
            for( unsigned r = 0; r < indices1.size(); ++r )
            {
                it = std::find( indices2.begin(), indices2.end(), indices1[r] );
                CPPUNIT_ASSERT( it != indices2.end() );
                unsigned r2 = it - indices2.begin();
                // if (!utest_mat_equal(Hess1[h],Hess2[r2],0.001))
                //  assert(false);
                CPPUNIT_ASSERT_MATRICES_EQUAL( Hess1[h], Hess2[r2], 0.05 );
                h += indices1.size() - r;
            }
        }
    }
}