MOAB: Mesh Oriented datABase  (version 5.4.1)
MBMesquite::NonSmoothDescent Class Reference

#include <NonSmoothDescent.hpp>

+ Inheritance diagram for MBMesquite::NonSmoothDescent:
+ Collaboration diagram for MBMesquite::NonSmoothDescent:

Classes

struct  ActiveSet
class  SymmetricMatrix

Public Member Functions

MESQUITE_EXPORT NonSmoothDescent (QualityMetric *qm)
virtual MESQUITE_EXPORT ~NonSmoothDescent ()
virtual MESQUITE_EXPORT std::string get_name () const
 Get string name for use in diagnostic and status output.
virtual MESQUITE_EXPORT PatchSetget_patch_set ()

Protected Member Functions

virtual MESQUITE_EXPORT void initialize (PatchData &pd, MsqError &err)
virtual MESQUITE_EXPORT void optimize_vertex_positions (PatchData &pd, MsqError &err)
virtual MESQUITE_EXPORT void initialize_mesh_iteration (PatchData &pd, MsqError &err)
virtual MESQUITE_EXPORT void terminate_mesh_iteration (PatchData &pd, MsqError &err)
virtual MESQUITE_EXPORT void cleanup ()

Private Types

enum  OptStatus {
  MSQ_NO_STATUS = 0, MSQ_STEP_ACCEPTED = 100, MSQ_IMP_TOO_SMALL, MSQ_FLAT_NO_IMP,
  MSQ_STEP_TOO_SMALL, MSQ_EQUILIBRIUM, MSQ_ZERO_SEARCH, MSQ_MAX_ITER_EXCEEDED
}
enum  Status {
  MSQ_NO_EQUIL = 101, MSQ_CHECK_TOP_DOWN, MSQ_CHECK_BOTTOM_UP, MSQ_TWO_PT_PLANE,
  MSQ_THREE_PT_PLANE, MSQ_CHECK_Y_COORD_DIRECTION, MSQ_CHECK_X_COORD_DIRECTION, MSQ_CHECK_Z_COORD_DIRECTION,
  MSQ_EQUIL, MSQ_HULL_TEST_ERROR
}
enum  Direction { MSQ_XDIR = 0, MSQ_YDIR = 1, MSQ_ZDIR = 2 }

Private Member Functions

void init_opt (PatchData &pd, MsqError &err)
void init_max_step_length (PatchData &pd, MsqError &err)
void minmax_opt (PatchData &pd, MsqError &err)
void step_acceptance (PatchData &pd, int iter_count, const Vector3D &search_dir, MsqError &err)
void get_min_estimate (double *final_est, MsqError &err)
void get_gradient_projections (const Vector3D &mSearch, MsqError &err)
void compute_alpha (MsqError &err)
bool compute_function (PatchData *pd, std::vector< double > &function, MsqError &err)
bool compute_gradient (PatchData *pd, std::vector< Vector3D > &grad_out, MsqError &err)
void find_active_set (const std::vector< double > &function, ActiveSet &active_set)
void print_active_set (const ActiveSet &active_set, const std::vector< double > &func)
bool improvement_check (MsqError &err)
bool validity_check (PatchData &pd, MsqError &err)
bool check_equilibrium (OptStatus &opt_status, MsqError &err)
bool convex_hull_test (const std::vector< Vector3D > &vec, MsqError &err)
bool check_vector_dots (const std::vector< Vector3D > &vec, const Vector3D &normal, MsqError &err)
void find_plane_normal (const Vector3D &pt1, const Vector3D &pt2, const Vector3D &pt3, Vector3D &cross, MsqError &)
void find_plane_points (Direction dir1, Direction dir2, const std::vector< Vector3D > &vec, Vector3D &pt1, Vector3D &pt2, Vector3D &pt3, Status &status, MsqError &)
void form_grammian (const std::vector< Vector3D > &vec, MsqError &err)
void form_PD_grammian (MsqError &err)
void singular_test (int n, const Matrix3D &A, bool &singular, MsqError &err)
bool solve2x2 (double a11, double a12, double a21, double a22, double b1, double b2, double x[2], MsqError &err)
void form_reduced_matrix (SymmetricMatrix &P)
void search_direction (PatchData &pd, Vector3D &search_dir_out, MsqError &err)
void search_edges_faces (const Vector3D *dir, Vector3D &search, MsqError &err)
void get_active_directions (const std::vector< Vector3D > &gradient, std::vector< Vector3D > &dir, MsqError &err)
 NonSmoothDescent (const NonSmoothDescent &pd)
NonSmoothDescentoperator= (const NonSmoothDescent &pd)

Private Attributes

size_t freeVertexIndex
double minStepSize
VertexPatches patchSet
QualityMetriccurrentQM
double originalValue
std::vector< double > mFunction
std::vector< double > testFunction
std::vector< double > originalFunction
std::vector< double > mGS
SymmetricMatrix mG
Matrix3D mPDG
std::vector< double > prevActiveValues
std::vector< Vector3DmGradient
std::vector< Vector3DtmpGrad
std::vector< size_t > tmpIdx
std::vector< size_t > qmHandles
OptStatus optStatus
size_t mSteepest
double mAlpha
double maxAlpha
int pdgInd [3]
ActiveSet mActive
ActiveSet testActive
ActiveSet originalActive

Detailed Description

Definition at line 55 of file NonSmoothDescent.hpp.


Member Enumeration Documentation

Enumerator:
MSQ_XDIR 
MSQ_YDIR 
MSQ_ZDIR 

Definition at line 165 of file NonSmoothDescent.hpp.

    {
        MSQ_XDIR = 0,
        MSQ_YDIR = 1,
        MSQ_ZDIR = 2
    };
Enumerator:
MSQ_NO_STATUS 
MSQ_STEP_ACCEPTED 
MSQ_IMP_TOO_SMALL 
MSQ_FLAT_NO_IMP 
MSQ_STEP_TOO_SMALL 
MSQ_EQUILIBRIUM 
MSQ_ZERO_SEARCH 
MSQ_MAX_ITER_EXCEEDED 

Definition at line 139 of file NonSmoothDescent.hpp.

Enumerator:
MSQ_NO_EQUIL 
MSQ_CHECK_TOP_DOWN 
MSQ_CHECK_BOTTOM_UP 
MSQ_TWO_PT_PLANE 
MSQ_THREE_PT_PLANE 
MSQ_CHECK_Y_COORD_DIRECTION 
MSQ_CHECK_X_COORD_DIRECTION 
MSQ_CHECK_Z_COORD_DIRECTION 
MSQ_EQUIL 
MSQ_HULL_TEST_ERROR 

Definition at line 151 of file NonSmoothDescent.hpp.


Constructor & Destructor Documentation

Definition at line 56 of file NonSmoothDescent.cpp.

References MSQ_DBGOUT.

                                                      : currentQM( qm )
{

    MSQ_DBGOUT( 1 ) << "- Executed NonSmoothDescent::NonSmoothDescent()\n";
}

Definition at line 61 of file NonSmoothDescent.hpp.

{}

Member Function Documentation

bool MBMesquite::NonSmoothDescent::check_equilibrium ( OptStatus opt_status,
MsqError err 
) [private]

Definition at line 1236 of file NonSmoothDescent.cpp.

References MBMesquite::NonSmoothDescent::ActiveSet::active_ind, convex_hull_test(), get_active_directions(), MBMesquite::length(), mActive, mGradient, MSQ_EQUILIBRIUM, MSQ_ERRZERO, MSQ_PRINT, and qmHandles.

Referenced by minmax_opt().

{
    std::vector< Vector3D > dir;

    // TODO - this subroutine is no longer clear to me... is it still
    // appropriate for quads and hexes?  I think it might be in 2D, but
    // 3D is less clear.  Is there a more general algorithm to use?
    // ask Todd/check in numerical optimization

    bool equil              = false;
    const size_t num_active = mActive.active_ind.size();
    ;

    if( num_active == qmHandles.size() )
    {
        equil  = true;
        status = MSQ_EQUILIBRIUM;
        MSQ_PRINT( 3 )( "All the function values are in the active set\n" );
    }

    /* set up the active mGradient directions */
    this->get_active_directions( mGradient, dir, err );
    MSQ_ERRZERO( err );

    /* normalize the active directions */
    for( size_t j = 0; j < num_active; j++ )
        dir[j] /= dir[j].length();

    equil = this->convex_hull_test( dir, err );
    if( equil ) status = MSQ_EQUILIBRIUM;

    return equil;
}
bool MBMesquite::NonSmoothDescent::check_vector_dots ( const std::vector< Vector3D > &  vec,
const Vector3D normal,
MsqError err 
) [private]

Definition at line 1122 of file NonSmoothDescent.cpp.

References moab::dot(), and MBMesquite::EPSILON.

Referenced by convex_hull_test().

{
    double test_dot = vec[0] % normal;
    unsigned ind    = 1;
    while( ( fabs( test_dot ) < EPSILON ) && ( ind < vec.size() ) )
    {
        test_dot = vec[ind] % normal;
        ind++;
    }

    for( unsigned i = ind; i < vec.size(); i++ )
    {
        double dot = vec[i] % normal;
        if( ( ( dot > 0 && test_dot < 0 ) || ( dot < 0 && test_dot > 0 ) ) && ( fabs( dot ) > EPSILON ) )
        {
            return true;
        }
    }
    return false;
}
void MBMesquite::NonSmoothDescent::cleanup ( ) [protected, virtual]

Implements MBMesquite::VertexMover.

Definition at line 134 of file NonSmoothDescent.cpp.

References MSQ_DBGOUT.

{
    MSQ_DBGOUT( 1 ) << "- Executing NonSmoothDescent::cleanup()\n";
    MSQ_DBGOUT( 1 ) << "- Done with NonSmoothDescent::cleanup()\n";
}

Definition at line 1583 of file NonSmoothDescent.cpp.

References moab::E, mAlpha, maxAlpha, mFunction, mGS, minStepSize, MSQ_PRINT, mSteepest, and qmHandles.

Referenced by minmax_opt().

{
    double steepest_function;
    double steepest_grad;
    double alpha_i;
    double min_positive_value = HUGE_VAL;

    //    FUNCTION_TIMER_START("Compute Alpha");

    MSQ_PRINT( 2 )( "In compute alpha\n" );

    mAlpha = HUGE_VAL;

    steepest_function = mFunction[mSteepest];
    steepest_grad     = mGS[mSteepest];
    for( size_t i = 0; i < qmHandles.size(); i++ )
    {
        /* if it's not active */
        if( i != mSteepest )
        {
            alpha_i = steepest_function - mFunction[i];

            if( fabs( mGS[mSteepest] - mGS[i] ) > 1E-13 )
            {
                /* compute line intersection */
                alpha_i = alpha_i / ( steepest_grad - mGS[i] );
            }
            else
            {
                /* the lines don't intersect - it's not under consideration*/
                alpha_i = 0;
            }
            if( ( alpha_i > minStepSize ) && ( fabs( alpha_i ) < fabs( mAlpha ) ) )
            {
                mAlpha = fabs( alpha_i );
                MSQ_PRINT( 3 )( "Setting alpha %lu %g\n", (unsigned long)i, alpha_i );
            }
            if( ( alpha_i > 0 ) && ( alpha_i < min_positive_value ) )
            {
                min_positive_value = alpha_i;
            }
        }
    }

    if( ( mAlpha == HUGE_VAL ) && ( min_positive_value != HUGE_VAL ) )
    {
        mAlpha = min_positive_value;
    }

    /* if it never gets set, set it to the default */
    if( mAlpha == HUGE_VAL )
    {
        mAlpha = maxAlpha;
        MSQ_PRINT( 3 )( "Setting alpha to the maximum step length\n" );
    }

    MSQ_PRINT( 3 )( "  The initial step size: %f\n", mAlpha );

    //    FUNCTION_TIMER_END();
}
bool MBMesquite::NonSmoothDescent::compute_function ( PatchData pd,
std::vector< double > &  function,
MsqError err 
) [private]

Definition at line 907 of file NonSmoothDescent.cpp.

References currentQM, MBMesquite::QualityMetric::evaluate(), MSQ_ERRZERO, and qmHandles.

Referenced by compute_gradient(), optimize_vertex_positions(), and step_acceptance().

{
    // NEED 1.0/FUNCTION WHICH IS ONLY TRUE OF CONDITION NUMBER

    func.resize( qmHandles.size() );
    bool valid_bool = true;
    for( size_t i = 0; i < qmHandles.size(); i++ )
    {
        valid_bool = valid_bool && currentQM->evaluate( *patch_data, qmHandles[i], func[i], err );
        MSQ_ERRZERO( err );
    }

    return valid_bool;
}
bool MBMesquite::NonSmoothDescent::compute_gradient ( PatchData pd,
std::vector< Vector3D > &  grad_out,
MsqError err 
) [private]

Definition at line 922 of file NonSmoothDescent.cpp.

References compute_function(), currentQM, MBMesquite::QualityMetric::evaluate_with_gradient(), freeVertexIndex, mGradient, MBMesquite::PatchData::move_vertex(), MSQ_DBGOUT, MSQ_ERRZERO, qmHandles, MBMesquite::PatchData::set_vertex_coordinates(), tmpGrad, tmpIdx, value(), and MBMesquite::PatchData::vertex_by_index().

Referenced by minmax_opt().

{
    MSQ_DBGOUT( 2 ) << "Computing Gradient\n";

    bool valid = true;

#ifdef NUMERICAL_GRADIENT

    const double delta = 10e-6;
    std::vector< double > func( qmHandles.size() ), fdelta( qmHandles.size() );

    valid = this->compute_function( patch_data, func, err );
    MSQ_ERRZERO( err );

    /* gradient in the x, y, z direction */

    for( int j = 0; j < 3; j++ )
    {

        // perturb the coordinates of the free vertex in the j direction by delta
        Vector3D delta_3( 0, 0, 0 );
        Vector3D orig_pos = patch_data->vertex_by_index( freeVertexIndex );
        delta_3[j]        = delta;
        patch_data->move_vertex( delta_3, freeVertexIndex, err );

        // compute the function at the perturbed point location
        valid = valid && this->compute_function( patch_data, fdelta, err );
        MSQ_ERRZERO( err );

        // compute the numerical gradient
        for( size_t i = 0; i < qmHandles.size(); i++ )
        {
            mGradient[i][j] = ( fdelta[i] - func[i] ) / delta;
            // MSQ_DEBUG_ACTION(3,{fprintf(stdout,"  Gradient
            // value[%d][%d]=%g\n",i,j,mGradient[i][j]);});
        }

        // put the coordinates back where they belong
        patch_data->set_vertex_coordinates( orig_pos, freeVertexIndex, err );
    }

#else
    double value;
    gradient_out.resize( qmHandles.size() );
    for( size_t i = 0; i < qmHandles.size(); i++ )
    {
        valid = valid && currentQM->evaluate_with_gradient( *patch_data, qmHandles[i], value, tmpIdx, tmpGrad, err );
        MSQ_ERRZERO( err );
        assert( tmpIdx.size() == 1 && tmpIdx[0] == freeVertexIndex );
        gradient_out[i] = tmpGrad[0];
    }
#endif

    return valid;
}
bool MBMesquite::NonSmoothDescent::convex_hull_test ( const std::vector< Vector3D > &  vec,
MsqError err 
) [private]

Definition at line 1145 of file NonSmoothDescent.cpp.

References check_vector_dots(), find_plane_normal(), find_plane_points(), MBMesquite::MsqError::INVALID_STATE, MSQ_CHECK_X_COORD_DIRECTION, MSQ_CHECK_Y_COORD_DIRECTION, MSQ_CHECK_Z_COORD_DIRECTION, MSQ_EQUIL, MSQ_HULL_TEST_ERROR, MSQ_NO_EQUIL, MSQ_PRINT, MSQ_SETERR, MSQ_THREE_PT_PLANE, MSQ_TWO_PT_PLANE, MSQ_XDIR, MSQ_YDIR, and MSQ_ZDIR.

Referenced by check_equilibrium().

{
    //    int ierr;
    bool equil         = false;
    Direction dir_done = MSQ_ZDIR;
    Status status      = MSQ_CHECK_Z_COORD_DIRECTION;
    Vector3D pt1, pt2, pt3, normal;

    /* tries to determine equilibrium for the 3D case */

    if( vec.size() <= 2 ) status = MSQ_NO_EQUIL;

    while( ( status != MSQ_EQUIL ) && ( status != MSQ_NO_EQUIL ) && ( status != MSQ_HULL_TEST_ERROR ) )
    {
        if( status == MSQ_CHECK_Z_COORD_DIRECTION )
        {
            this->find_plane_points( MSQ_ZDIR, MSQ_YDIR, vec, pt1, pt2, pt3, status, err );
            dir_done = MSQ_ZDIR;
        }
        if( status == MSQ_CHECK_Y_COORD_DIRECTION )
        {
            this->find_plane_points( MSQ_YDIR, MSQ_XDIR, vec, pt1, pt2, pt3, status, err );
            dir_done = MSQ_YDIR;
        }
        if( status == MSQ_CHECK_X_COORD_DIRECTION )
        {
            this->find_plane_points( MSQ_XDIR, MSQ_ZDIR, vec, pt1, pt2, pt3, status, err );
            dir_done = MSQ_XDIR;
        }
        if( status == MSQ_TWO_PT_PLANE )
        {
            pt3 = Vector3D( 0, 0, 0 );
        }
        if( ( status == MSQ_TWO_PT_PLANE ) || ( status == MSQ_THREE_PT_PLANE ) )
        {
            this->find_plane_normal( pt1, pt2, pt3, normal, err );
            equil = this->check_vector_dots( vec, normal, err );
            if( equil )
            {
                switch( dir_done )
                {
                    case MSQ_ZDIR:
                        equil  = false;
                        status = MSQ_CHECK_Y_COORD_DIRECTION;
                        break;
                    case MSQ_YDIR:
                        equil  = false;
                        status = MSQ_CHECK_X_COORD_DIRECTION;
                        break;
                    case MSQ_XDIR:
                        equil  = true;
                        status = MSQ_EQUIL;
                }
            }
            else if( equil == 0 )
            {
                status = MSQ_NO_EQUIL;
            }
            else
            {
                MSQ_SETERR( err )( "equil flag not set to 0 or 1", MsqError::INVALID_STATE );
            }
        }
    }
    switch( status )
    {
        case MSQ_NO_EQUIL:
            MSQ_PRINT( 3 )( "Not an equilibrium point\n" );
            equil = false;
            break;
        case MSQ_EQUIL:
            MSQ_PRINT( 3 )( "An equilibrium point\n" );
            equil = true;
            break;
        default:
            MSQ_PRINT( 3 )( "Failed to determine equil or not; status = %d\n", status );
    }
    //    FUNCTION_TIMER_END();
    return ( equil );
}
void MBMesquite::NonSmoothDescent::find_active_set ( const std::vector< double > &  function,
ActiveSet active_set 
) [private]

Definition at line 978 of file NonSmoothDescent.cpp.

References MBMesquite::NonSmoothDescent::ActiveSet::active_ind, MBMesquite::NonSmoothDescent::ActiveSet::add(), MBMesquite::EPSILON, MSQ_DBGOUT, qmHandles, MBMesquite::NonSmoothDescent::ActiveSet::set(), and MBMesquite::NonSmoothDescent::ActiveSet::true_active_value.

Referenced by minmax_opt(), optimize_vertex_positions(), and step_acceptance().

{
    // local parameter initialization
    const double activeEpsilon = .3e-4;
    //  activeEpsilon = .3e-8;

    double function_val;
    double active_value0;
    double temp;

    //    FUNCTION_TIMER_START("Find Active Set");
    MSQ_DBGOUT( 2 ) << "\nFinding the active set\n";

    /* the first function value is our initial active value */
    active_set.set( 0 );
    active_set.true_active_value = function[0];
    //    MSQ_DEBUG_ACTION(3,{fprintf(stdout,"  function value[0]: %g\n",function[0]);});

    /* first sort out the active set...
       all vals within active_epsilon of largest val */

    for( size_t i = 1; i < qmHandles.size(); i++ )
    {
        function_val = function[i];
        if( active_set.true_active_value < function_val ) active_set.true_active_value = function_val;
        active_value0 = function[active_set.active_ind[0]];
        temp          = fabs( function_val - active_value0 );
        //        MSQ_DEBUG_ACTION(3,{fprintf(stdout,"  function value[%d]: %g\n",i,function[i]);});
        if( function_val > active_value0 )
        {  // seek max_i function[i]
            if( temp >= activeEpsilon )
            {
                active_set.set( i );  // new max
            }
            else
            {
                active_set.add( i, fabs( function_val - active_value0 ) < EPSILON );
            }
        }
        else
        {
            if( temp < activeEpsilon )
            {
                active_set.add( i, fabs( function_val - active_value0 ) < EPSILON );
            }
        }
    }
}
void MBMesquite::NonSmoothDescent::find_plane_normal ( const Vector3D pt1,
const Vector3D pt2,
const Vector3D pt3,
Vector3D cross,
MsqError  
) [inline, private]

Definition at line 251 of file NonSmoothDescent.hpp.

References MBMesquite::Vector3D::length().

Referenced by convex_hull_test().

{
    Vector3D vecA = pt2 - pt1;
    Vector3D vecB = pt3 - pt1;
    cross         = vecA * vecB;
    cross /= cross.length();
}
void MBMesquite::NonSmoothDescent::find_plane_points ( Direction  dir1,
Direction  dir2,
const std::vector< Vector3D > &  vec,
Vector3D pt1,
Vector3D pt2,
Vector3D pt3,
Status status,
MsqError  
) [private]

Definition at line 140 of file NonSmoothDescent.cpp.

References MBMesquite::CLOCKWISE, MBMesquite::COUNTERCLOCKWISE, MBMesquite::EPSILON, MSQ_CHECK_BOTTOM_UP, MSQ_CHECK_TOP_DOWN, MSQ_CHECK_X_COORD_DIRECTION, MSQ_CHECK_Y_COORD_DIRECTION, MSQ_EQUIL, MSQ_HULL_TEST_ERROR, MSQ_NO_EQUIL, MSQ_PRINT, and MSQ_TWO_PT_PLANE.

Referenced by convex_hull_test().

{
    int i;
    int num_min, num_max;
    Rotate rotate   = CLOCKWISE;
    int num_rotated = 0;
    double pt_1, pt_2;
    double min, inv_slope;
    double min_inv_slope = 0.;
    double max;
    double max_inv_slope    = 0;
    double inv_origin_slope = 0;
    const int num_vec       = vec.size();
    const int auto_ind_size = 50;
    int auto_ind[auto_ind_size];
    std::vector< int > heap_ind;
    int* ind;
    if( num_vec <= auto_ind_size )
        ind = auto_ind;
    else
    {
        heap_ind.resize( num_vec );
        ind = &heap_ind[0];
    }

    status = MSQ_CHECK_BOTTOM_UP;
    /* find the minimum points in dir1 starting at -1 */
    num_min = 0;
    ind[0]  = -1;
    ind[1]  = -1;
    ind[2]  = -1;
    min     = 1.0;
    for( i = 0; i < num_vec; i++ )
    {
        if( vec[i][dir1] < min )
        {
            min     = vec[i][dir1];
            ind[0]  = i;
            num_min = 1;
        }
        else if( fabs( vec[i][dir1] - min ) < EPSILON )
        {
            ind[num_min++] = i;
        }
    }
    if( min >= 0 ) status = MSQ_NO_EQUIL;

    if( status != MSQ_NO_EQUIL )
    {
        switch( num_min )
        {
            case 1: /* rotate to find the next point */
                pt1  = vec[ind[0]];
                pt_1 = pt1[dir1];
                pt_2 = pt1[dir2];
                if( pt1[dir2] <= 0 )
                {
                    rotate        = COUNTERCLOCKWISE;
                    max_inv_slope = -HUGE_VAL;
                }
                if( pt1[dir2] > 0 )
                {
                    rotate        = CLOCKWISE;
                    min_inv_slope = HUGE_VAL;
                }
                switch( rotate )
                {
                    case COUNTERCLOCKWISE:
                        for( i = 0; i < num_vec; i++ )
                        {
                            if( i != ind[0] )
                            {
                                inv_slope = ( vec[i][dir2] - pt_2 ) / ( vec[i][dir1] - pt_1 );
                                if( ( inv_slope > max_inv_slope ) && ( fabs( inv_slope - max_inv_slope ) > EPSILON ) )
                                {
                                    ind[1]        = i;
                                    max_inv_slope = inv_slope;
                                    num_rotated   = 1;
                                }
                                else if( fabs( inv_slope - max_inv_slope ) < EPSILON )
                                {
                                    ind[2] = i;
                                    num_rotated++;
                                }
                            }
                        }
                        break;
                    case CLOCKWISE:
                        for( i = 0; i < num_vec; i++ )
                        {
                            if( i != ind[0] )
                            {
                                inv_slope = ( vec[i][dir2] - pt_2 ) / ( vec[i][dir1] - pt_1 );
                                if( ( inv_slope < min_inv_slope ) && ( fabs( inv_slope - max_inv_slope ) > EPSILON ) )
                                {
                                    ind[1]        = i;
                                    min_inv_slope = inv_slope;
                                    num_rotated   = 1;
                                }
                                else if( fabs( inv_slope - min_inv_slope ) < EPSILON )
                                {
                                    ind[2] = i;
                                    num_rotated++;
                                }
                            }
                        }
                }
                switch( num_rotated )
                {
                    case 0:
                        MSQ_PRINT( 3 )( "No points in the rotation ... odd\n" );
                        status = MSQ_HULL_TEST_ERROR;
                        break;
                    case 1:
                        MSQ_PRINT( 3 )( "Found a line in the convex hull\n" );
                        pt2    = vec[ind[1]];
                        status = MSQ_TWO_PT_PLANE;
                        break;
                    default:
                        MSQ_PRINT( 3 )( "Found 2 or more points in the rotation\n" );
                        if( fabs( pt_1 ) > EPSILON ) inv_origin_slope = pt_2 / pt_1;
                        switch( rotate )
                        {
                            case COUNTERCLOCKWISE:
                                if( inv_origin_slope >= max_inv_slope )
                                    status = MSQ_NO_EQUIL;
                                else
                                    status = MSQ_CHECK_TOP_DOWN;
                                break;
                            case CLOCKWISE:
                                if( inv_origin_slope <= min_inv_slope )
                                    status = MSQ_NO_EQUIL;
                                else
                                    status = MSQ_CHECK_TOP_DOWN;
                        }
                }
                break;
            case 2: /* use these two points to define the plane */
                MSQ_PRINT( 3 )( "Found two minimum points to define the plane\n" );
                pt1    = vec[ind[0]];
                pt2    = vec[ind[1]];
                status = MSQ_TWO_PT_PLANE;
                break;
            default: /* check to see if all > 0 */
                MSQ_PRINT( 3 )( "Found 3 or more points in min plane %f\n", min );
                if( vec[ind[0]][dir1] >= 0 )
                    status = MSQ_NO_EQUIL;
                else
                    status = MSQ_CHECK_TOP_DOWN;
        }
    }

    /***************************/
    /*  failed to find any information, checking top/down this coord*/
    /***************************/

    if( status == MSQ_CHECK_TOP_DOWN )
    {
        /* find the maximum points in dir1 starting at 1 */
        num_max = 0;
        ind[0]  = -1;
        ind[1]  = -1;
        ind[2]  = -1;
        max     = -1.0;
        for( i = 0; i < num_vec; i++ )
        {
            if( vec[i][dir1] > max )
            {
                max     = vec[i][dir1];
                ind[0]  = i;
                num_max = 1;
            }
            else if( fabs( vec[i][dir1] - max ) < EPSILON )
            {
                ind[num_max++] = i;
            }
        }
        if( max <= 0 ) status = MSQ_NO_EQUIL;

        if( status != MSQ_NO_EQUIL )
        {
            switch( num_max )
            {
                case 1: /* rotate to find the next point */
                    pt1  = vec[ind[0]];
                    pt_1 = pt1[dir1];
                    pt_2 = pt1[dir2];
                    if( pt1[dir2] < 0 )
                    {
                        rotate        = CLOCKWISE;
                        min_inv_slope = HUGE_VAL;
                    }
                    if( pt1[dir2] >= 0 )
                    {
                        rotate        = COUNTERCLOCKWISE;
                        max_inv_slope = -HUGE_VAL;
                    }
                    switch( rotate )
                    {
                        case COUNTERCLOCKWISE:
                            for( i = 0; i < num_vec; i++ )
                            {
                                if( i != ind[0] )
                                {
                                    inv_slope = ( vec[i][dir2] - pt_2 ) / ( vec[i][dir1] - pt_1 );
                                    if( inv_slope > max_inv_slope )
                                    {
                                        ind[1]        = i;
                                        max_inv_slope = inv_slope;
                                        num_rotated   = 1;
                                    }
                                    else if( fabs( inv_slope - max_inv_slope ) < EPSILON )
                                    {
                                        ind[2] = i;
                                        num_rotated++;
                                    }
                                }
                            }
                            break;
                        case CLOCKWISE:
                            for( i = 0; i < num_vec; i++ )
                            {
                                if( i != ind[0] )
                                {
                                    inv_slope = ( vec[i][dir2] - pt_2 ) / ( vec[i][dir1] - pt_1 );
                                    if( inv_slope < min_inv_slope )
                                    {
                                        ind[1]        = i;
                                        min_inv_slope = inv_slope;
                                        num_rotated   = 1;
                                    }
                                    else if( fabs( inv_slope - min_inv_slope ) < EPSILON )
                                    {
                                        ind[2] = i;
                                        num_rotated++;
                                    }
                                }
                            }
                    }
                    switch( num_rotated )
                    {
                        case 0:
                            MSQ_PRINT( 3 )( "No points in the rotation ... odd\n" );
                            status = MSQ_HULL_TEST_ERROR;
                            break;
                        case 1:
                            MSQ_PRINT( 3 )( "Found a line in the convex hull\n" );
                            pt2    = vec[ind[1]];
                            status = MSQ_TWO_PT_PLANE;
                            break;
                        default:
                            MSQ_PRINT( 3 )( "Found 2 or more points in the rotation\n" );
                            /* check to see if rotation got past origin */
                            inv_origin_slope = pt_2 / pt_1;
                            switch( rotate )
                            {
                                case COUNTERCLOCKWISE:
                                    if( inv_origin_slope >= max_inv_slope )
                                        status = MSQ_NO_EQUIL;
                                    else if( dir1 == 2 )
                                        status = MSQ_CHECK_Y_COORD_DIRECTION;
                                    else if( dir1 == 1 )
                                        status = MSQ_CHECK_X_COORD_DIRECTION;
                                    else
                                        status = MSQ_EQUIL;
                                    break;
                                case CLOCKWISE:
                                    if( inv_origin_slope <= min_inv_slope )
                                        status = MSQ_NO_EQUIL;
                                    else if( dir1 == 2 )
                                        status = MSQ_CHECK_Y_COORD_DIRECTION;
                                    else if( dir1 == 1 )
                                        status = MSQ_CHECK_X_COORD_DIRECTION;
                                    else
                                        status = MSQ_EQUIL;
                            }
                    }
                    break;
                case 2: /* use these two points to define the plane */
                    pt1    = vec[ind[0]];
                    pt2    = vec[ind[1]];
                    status = MSQ_TWO_PT_PLANE;
                    break;
                default: /* check to see if all > 0 */
                    MSQ_PRINT( 3 )( "Found 3 in max plane %f\n", max );
                    if( vec[ind[0]][dir1] <= 0 )
                        status = MSQ_NO_EQUIL;
                    else if( dir1 == 2 )
                        status = MSQ_CHECK_Y_COORD_DIRECTION;
                    else if( dir1 == 1 )
                        status = MSQ_CHECK_X_COORD_DIRECTION;
                    else
                        status = MSQ_EQUIL;
            }
        }
    }
}
void MBMesquite::NonSmoothDescent::form_grammian ( const std::vector< Vector3D > &  vec,
MsqError err 
) [private]

Definition at line 1226 of file NonSmoothDescent.cpp.

References MBMesquite::NonSmoothDescent::ActiveSet::active_ind, mActive, mG, and MBMesquite::NonSmoothDescent::SymmetricMatrix::resize().

Referenced by search_direction().

{
    /* form the grammian with the dot products of the gradients */
    const size_t num_active = mActive.active_ind.size();
    mG.resize( num_active );
    for( size_t i = 0; i < num_active; i++ )
        for( size_t j = i; j < num_active; j++ )
            mG( i, j ) = vec[i] % vec[j];
}

Definition at line 1378 of file NonSmoothDescent.cpp.

References MBMesquite::NonSmoothDescent::ActiveSet::active_ind, MBMesquite::MsqError::INVALID_STATE, mActive, mG, mPDG, MSQ_SETERR, pdgInd, and singular_test().

Referenced by search_direction().

{
    int i, j, k;
    int g_ind_1;
    bool singular = false;

    const int num_active = mActive.active_ind.size();

    /* this assumes the grammian has been formed */
    for( i = 0; i < num_active; i++ )
    {
        for( j = i; j < num_active; j++ )
        {
            if( mG( i, j ) == -1 )
            {
                MSQ_SETERR( err )( "Grammian not computed properly", MsqError::INVALID_STATE );
                return;
            }
        }
    }

    /* use the first gradient in the active set */
    g_ind_1    = 0;
    mPDG[0][0] = mG( 0, 0 );
    pdgInd[0]  = mActive.active_ind[0];

    /* test the rest and add them as appropriate */
    k = 1;
    i = 1;
    while( ( k < 3 ) && ( i < num_active ) )
    {
        mPDG[0][k] = mPDG[k][0] = mG( 0, i );
        mPDG[k][k]              = mG( i, i );
        if( k == 2 )
        { /* add the dot product of g1 and g2 */
            mPDG[1][k] = mPDG[k][1] = mG( g_ind_1, i );
        }
        this->singular_test( k + 1, mPDG, singular, err );
        if( !singular )
        {
            pdgInd[k] = mActive.active_ind[i];
            if( k == 1 ) g_ind_1 = i;
            k++;
        }
        i++;
    }
}

Definition at line 1537 of file NonSmoothDescent.cpp.

References MBMesquite::NonSmoothDescent::ActiveSet::active_ind, mActive, mG, MBMesquite::P, and MBMesquite::NonSmoothDescent::SymmetricMatrix::resize().

Referenced by search_direction().

{
    const size_t P_size = mActive.active_ind.size() - 1;
    P.resize( P_size );
    for( size_t i = 0; i < P_size; i++ )
    {
        P( i, i ) = mG( 0, 0 ) - 2 * mG( 0, i + 1 ) + mG( i + 1, i + 1 );
        for( size_t j = i + 1; j < P_size; j++ )
        {
            P( i, j ) = mG( 0, 0 ) - mG( 0, j + 1 ) - mG( i + 1, 0 ) + mG( i + 1, j + 1 );
        }
    }
}
void MBMesquite::NonSmoothDescent::get_active_directions ( const std::vector< Vector3D > &  gradient,
std::vector< Vector3D > &  dir,
MsqError err 
) [private]

Definition at line 1110 of file NonSmoothDescent.cpp.

References MBMesquite::NonSmoothDescent::ActiveSet::active_ind, and mActive.

Referenced by check_equilibrium(), and search_direction().

{
    const size_t num_active = mActive.active_ind.size();
    dir.resize( num_active );
    for( size_t i = 0; i < num_active; i++ )
    {
        dir[i] = pGradient[mActive.active_ind[i]];
    }
}
void MBMesquite::NonSmoothDescent::get_gradient_projections ( const Vector3D mSearch,
MsqError err 
) [private]

Definition at line 1575 of file NonSmoothDescent.cpp.

References mGradient, mGS, MSQ_PRINT, mSteepest, and qmHandles.

Referenced by minmax_opt().

{
    for( size_t i = 0; i < qmHandles.size(); i++ )
        mGS[i] = mGradient[i] % mSearch;

    MSQ_PRINT( 3 )( "steepest in get_gradient_projections %lu\n", (unsigned long)mSteepest );
}
void MBMesquite::NonSmoothDescent::get_min_estimate ( double *  final_est,
MsqError err 
) [private]

Definition at line 1551 of file NonSmoothDescent.cpp.

References MBMesquite::NonSmoothDescent::ActiveSet::active_ind, MBMesquite::EPSILON, mActive, mAlpha, mGS, and qmHandles.

Referenced by step_acceptance().

{
    double est_imp;

    *final_est = -HUGE_VAL;
    for( size_t i = 0; i < mActive.active_ind.size(); i++ )
    {
        est_imp = -mAlpha * mGS[mActive.active_ind[i]];
        if( est_imp > *final_est ) *final_est = est_imp;
    }
    if( *final_est == 0 )
    {
        *final_est = -HUGE_VAL;
        for( size_t i = 0; i < qmHandles.size(); i++ )
        {
            est_imp = -mAlpha * mGS[i];
            if( ( est_imp > *final_est ) && ( fabs( est_imp ) > EPSILON ) )
            {
                *final_est = est_imp;
            }
        }
    }
}
std::string MBMesquite::NonSmoothDescent::get_name ( ) const [virtual]

Get string name for use in diagnostic and status output.

Implements MBMesquite::Instruction.

Definition at line 62 of file NonSmoothDescent.cpp.

{
    return "NonSmoothDescent";
}

Implements MBMesquite::QualityImprover.

Definition at line 67 of file NonSmoothDescent.cpp.

References patchSet.

{
    return &patchSet;
}
bool MBMesquite::NonSmoothDescent::improvement_check ( MsqError err) [inline, private]

Definition at line 263 of file NonSmoothDescent.hpp.

References mActive, MSQ_PRINT, originalValue, and MBMesquite::NonSmoothDescent::ActiveSet::true_active_value.

Referenced by step_acceptance().

{
    /* check to see that the mesh didn't get worse */
    if( originalValue < mActive.true_active_value )
    {
        MSQ_PRINT( 2 )
        ( "The local mesh got worse; initial value %f; final value %f\n", originalValue, mActive.true_active_value );
        return false;
    }

    return true;
}

Definition at line 1693 of file NonSmoothDescent.cpp.

References MBMesquite::PatchData::get_vertex_array(), MBMesquite::MsqError::INVALID_MESH, MBMesquite::length_squared(), maxAlpha, MSQ_PRINT, MSQ_SETERR, MBMesquite::PatchData::num_elements(), and MBMesquite::PatchData::num_nodes().

Referenced by optimize_vertex_positions().

{
    size_t i, j;
    double max_diff = 0;
    double diff     = 0;

    MSQ_PRINT( 2 )( "In init_max_step_length\n" );

    /* check that the input data is correct */
    if( pd.num_elements() == 0 )
    {
        MSQ_SETERR( err )( "Num incident vtx = 0\n", MsqError::INVALID_MESH );
        return;
    }

    /* find the maximum distance between two incident vertex locations */
    const MsqVertex* coords = pd.get_vertex_array( err );
    for( i = 0; i < pd.num_nodes() - 1; i++ )
    {
        for( j = i; j < pd.num_nodes(); j++ )
        {
            diff = ( coords[i] - coords[j] ).length_squared();
            if( max_diff < diff ) max_diff = diff;
        }
    }
    max_diff = sqrt( max_diff );
    if( max_diff == 0 )
    {
        MSQ_SETERR( err )
        ( "Maximum distance between incident vertices = 0\n", MsqError::INVALID_MESH );
        return;
    }
    maxAlpha = max_diff / 100;

    MSQ_PRINT( 3 )( "  Maximum step is %g\n", maxAlpha );
}
void MBMesquite::NonSmoothDescent::init_opt ( PatchData pd,
MsqError err 
) [private]

Definition at line 1655 of file NonSmoothDescent.cpp.

References currentQM, MBMesquite::NonSmoothDescent::SymmetricMatrix::fill(), MBMesquite::QualityMetric::get_evaluations(), mAlpha, maxAlpha, mFunction, mG, mGradient, mGS, mPDG, MSQ_ERRRTN, MSQ_NO_STATUS, MSQ_PRINT, mSteepest, optStatus, originalFunction, pdgInd, prevActiveValues, qmHandles, MBMesquite::NonSmoothDescent::SymmetricMatrix::resize(), and testFunction.

Referenced by optimize_vertex_positions().

{
    qmHandles.clear();
    currentQM->get_evaluations( pd, qmHandles, true, err );MSQ_ERRRTN( err );

    MSQ_PRINT( 2 )( "\nInitializing Optimization \n" );

    /* for the purposes of initialization will be set to zero after */
    optStatus = MSQ_NO_STATUS;
    mSteepest = 0;
    mAlpha    = 0;
    maxAlpha  = 0;

    MSQ_PRINT( 3 )( "  Initialized Constants \n" );
    pdgInd[0] = pdgInd[1] = pdgInd[2] = -1;
    mPDG                              = Matrix3D( 0.0 );

    MSQ_PRINT( 3 )( "  Initialized search and PDG \n" );
    mFunction.clear();
    mFunction.resize( qmHandles.size(), 0.0 );
    testFunction.clear();
    testFunction.resize( qmHandles.size(), 0.0 );
    originalFunction.clear();
    originalFunction.resize( qmHandles.size(), 0.0 );
    mGradient.clear();
    mGradient.resize( qmHandles.size(), Vector3D( 0.0 ) );
    mGS.resize( qmHandles.size() );

    MSQ_PRINT( 3 )( "  Initialized function/gradient \n" );
    mG.resize( qmHandles.size() );
    mG.fill( -1 );
    MSQ_PRINT( 3 )( "  Initialized G\n" );

    prevActiveValues.clear();
    prevActiveValues.reserve( 32 );
    MSQ_PRINT( 3 )( "  Initialized prevActiveValues\n" );
}
void MBMesquite::NonSmoothDescent::initialize ( PatchData pd,
MsqError err 
) [protected, virtual]

Implements MBMesquite::VertexMover.

Definition at line 72 of file NonSmoothDescent.cpp.

References minStepSize, and MSQ_DBGOUT.

{
    minStepSize = 1e-6;
    MSQ_DBGOUT( 1 ) << "- Executed NonSmoothDescent::initialize()\n";
}
void MBMesquite::NonSmoothDescent::initialize_mesh_iteration ( PatchData pd,
MsqError err 
) [protected, virtual]

Implements MBMesquite::VertexMover.

Definition at line 78 of file NonSmoothDescent.cpp.

{}
void MBMesquite::NonSmoothDescent::minmax_opt ( PatchData pd,
MsqError err 
) [private]

Definition at line 577 of file NonSmoothDescent.cpp.

References MBMesquite::NonSmoothDescent::ActiveSet::active_ind, check_equilibrium(), compute_alpha(), compute_gradient(), find_active_set(), get_gradient_projections(), mActive, mFunction, mGradient, MSQ_DBG, MSQ_EQUILIBRIUM, MSQ_ERRRTN, MSQ_FLAT_NO_IMP, MSQ_FUNCTION_TIMER, MSQ_IMP_TOO_SMALL, MSQ_MAX_ITER_EXCEEDED, MBMesquite::MSQ_MAX_OPT_ITER, MSQ_PRINT, MSQ_STEP_TOO_SMALL, MSQ_ZERO_SEARCH, optStatus, originalFunction, originalValue, prevActiveValues, print_active_set(), search_direction(), step_acceptance(), MBMesquite::NonSmoothDescent::ActiveSet::true_active_value, and validity_check().

Referenced by optimize_vertex_positions().

{
    bool equilibriumPt = false;
    Vector3D mSearch( 0, 0, 0 );

    //      int valid;
    MSQ_FUNCTION_TIMER( "Minmax Opt" );
    MSQ_PRINT( 2 )( "In minmax_opt\n" );

    mFunction     = originalFunction;
    originalValue = mActive.true_active_value;

    int iterCount    = 0;
    int optIterCount = 0;

    MSQ_PRINT( 3 )( "Done copying original function to function\n" );

    this->find_active_set( mFunction, mActive );
    prevActiveValues.clear();
    prevActiveValues.push_back( mActive.true_active_value );

    /* check for equilibrium point */
    /* compute the gradient */
    this->compute_gradient( &pd, mGradient, err );MSQ_ERRRTN( err );

    if( mActive.active_ind.size() >= 2 )
    {
        MSQ_PRINT( 3 )( "Testing for an equilibrium point \n" );
        equilibriumPt = this->check_equilibrium( optStatus, err );MSQ_ERRRTN( err );

        if( MSQ_DBG( 2 ) && equilibriumPt ) MSQ_PRINT( 2 )( "Optimization Exiting: An equilibrium point \n" );
    }

    /* terminate if we have found an equilibrium point or if the step is
       too small to be worthwhile continuing */
    while( ( optStatus != MSQ_EQUILIBRIUM ) && ( optStatus != MSQ_STEP_TOO_SMALL ) &&
           ( optStatus != MSQ_IMP_TOO_SMALL ) && ( optStatus != MSQ_FLAT_NO_IMP ) && ( optStatus != MSQ_ZERO_SEARCH ) &&
           ( optStatus != MSQ_MAX_ITER_EXCEEDED ) )
    {

        /* increase the iteration count by one */
        /* smooth_param->iter_count += 1; */
        iterCount += 1;
        optIterCount += 1;
        if( iterCount > MSQ_MAX_OPT_ITER ) optStatus = MSQ_MAX_ITER_EXCEEDED;

        MSQ_PRINT( 3 )( "\nITERATION %d \n", iterCount );

        /* compute the gradient */
        this->compute_gradient( &pd, mGradient, err );MSQ_ERRRTN( err );

        MSQ_PRINT( 3 )( "Computing the search direction \n" );
        this->search_direction( pd, mSearch, err );MSQ_ERRRTN( err );

        /* if there are viable directions to search */
        if( ( optStatus != MSQ_ZERO_SEARCH ) && ( optStatus != MSQ_MAX_ITER_EXCEEDED ) )
        {

            MSQ_PRINT( 3 )( "Computing the projections of the gradients \n" );
            this->get_gradient_projections( mSearch, err );MSQ_ERRRTN( err );

            MSQ_PRINT( 3 )( "Computing the initial step size \n" );
            this->compute_alpha( err );MSQ_ERRRTN( err );

            MSQ_PRINT( 3 )( "Testing whether to accept this step \n" );
            this->step_acceptance( pd, iterCount, mSearch, err );MSQ_ERRRTN( err );
            // MSQ_PRINT(3)("The new free vertex position is %f %f %f\n",
            //  mCoords[freeVertexIndex][0],mCoords[freeVertexIndex][1],mCoords[freeVertexIndex][2]);

            if( MSQ_DBG( 3 ) )
            {
                /* Print the active set */
                this->print_active_set( mActive, mFunction );
            }

            /* check for equilibrium point */
            if( mActive.active_ind.size() >= 2 )
            {
                MSQ_PRINT( 3 )( "Testing for an equilibrium point \n" );
                equilibriumPt = this->check_equilibrium( optStatus, err );MSQ_ERRRTN( err );

                if( MSQ_DBG( 2 ) && equilibriumPt ) MSQ_PRINT( 2 )( "Optimization Exiting: An equilibrium point \n" );
            }

            /* record the values */
            prevActiveValues.push_back( mActive.true_active_value );
        }
        else
        {
            /* decrease the iteration count by one */
            /* smooth_param->iter_count -= 1; */
            iterCount -= 1;
            if( MSQ_DBG( 2 ) )
            {
                MSQ_PRINT( 2 )
                ( "Optimization Exiting: No viable directions; equilibrium point \n" );
                /* Print the old active set */
                this->print_active_set( mActive, mFunction );
            }
        }
    }

    MSQ_PRINT( 2 )( "Checking the validity of the mesh\n" );
    if( !this->validity_check( pd, err ) ) MSQ_PRINT( 2 )( "The final mesh is not valid\n" );MSQ_ERRRTN( err );

    MSQ_PRINT( 2 )( "Number of optimization iterations %d\n", iterCount );

    switch( optStatus )
    {
        default:
            MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Invalid early termination\n" );
            break;
        case MSQ_EQUILIBRIUM:
            MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Equilibrium\n" );
            break;
        case MSQ_STEP_TOO_SMALL:
            MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Step Too Small\n" );
            break;
        case MSQ_IMP_TOO_SMALL:
            MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Improvement Too Small\n" );
            break;
        case MSQ_FLAT_NO_IMP:
            MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Flat No Improvement\n" );
            break;
        case MSQ_ZERO_SEARCH:
            MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Zero Search\n" );
            break;
        case MSQ_MAX_ITER_EXCEEDED:
            MSQ_PRINT( 2 )( "Optimization Termination OptStatus: Max Iter Exceeded\n" );
            break;
    }
}
NonSmoothDescent& MBMesquite::NonSmoothDescent::operator= ( const NonSmoothDescent pd) [private]
void MBMesquite::NonSmoothDescent::optimize_vertex_positions ( PatchData pd,
MsqError err 
) [protected, virtual]

Implements MBMesquite::VertexMover.

Definition at line 79 of file NonSmoothDescent.cpp.

References compute_function(), find_active_set(), freeVertexIndex, init_max_step_length(), init_opt(), MBMesquite::MsqError::INVALID_MESH, mActive, minmax_opt(), MSQ_ERRRTN, MSQ_FUNCTION_TIMER, MSQ_PRINT, MSQ_SETERR, MBMesquite::MsqFreeVertexIndexIterator::next(), MBMesquite::PatchData::num_elements(), MBMesquite::PatchData::num_free_vertices(), MBMesquite::PatchData::num_nodes(), originalFunction, MBMesquite::MsqFreeVertexIndexIterator::reset(), validity_check(), and MBMesquite::MsqFreeVertexIndexIterator::value().

{
    MSQ_FUNCTION_TIMER( "NonSmoothDescent" );

    //  cout << "- Executing NonSmoothDescent::optimize_node_positions()\n";
    /* perform the min max smoothing algorithm */
    MSQ_PRINT( 2 )( "\nInitializing the patch iteration\n" );

    MSQ_PRINT( 3 )( "Number of Vertices: %d\n", (int)pd.num_nodes() );
    MSQ_PRINT( 3 )( "Number of Elements: %d\n", (int)pd.num_elements() );
    // Michael: Note: is this a reliable way to get the dimension?
    // NOTE: Mesquite only does 3-dimensional (including surface) meshes.
    // mDimension = pd.get_mesh()->get_geometric_dimension(err); MSQ_ERRRTN(err);
    // MSQ_PRINT(3)("Spatial Dimension: %d\n",mDimension);
    MSQ_PRINT( 3 )( "Spatial Dimension: 3\n" );

    MSQ_PRINT( 3 )( "Num Free = %d\n", (int)pd.num_free_vertices() );

    MsqFreeVertexIndexIterator free_iter( pd, err );MSQ_ERRRTN( err );
    free_iter.reset();
    free_iter.next();
    freeVertexIndex = free_iter.value();
    MSQ_PRINT( 3 )( "Free Vertex Index = %lu\n", (unsigned long)freeVertexIndex );

    // TODO - need to switch to validity via metric evaluations should
    // be associated with the compute_function somehow
    /* check for an invalid mesh; if it's invalid return and ask the user
       to use untangle */
    if( !this->validity_check( pd, err ) )
    {
        MSQ_PRINT( 1 )( "ERROR: Invalid mesh\n" );
        MSQ_SETERR( err )
        ( "Invalid Mesh: Use untangle to create a valid "
          "triangulation",
          MsqError::INVALID_MESH );
        return;
    }

    /* initialize the optimization data up to numFunctionValues */
    this->init_opt( pd, err );MSQ_ERRRTN( err );
    this->init_max_step_length( pd, err );MSQ_ERRRTN( err );
    MSQ_PRINT( 3 )( "Done initializing optimization\n" );

    /* compute the initial function values */
    // TODO this should return a bool with the validity
    this->compute_function( &pd, originalFunction, err );MSQ_ERRRTN( err );

    // find the initial active set
    this->find_active_set( originalFunction, mActive );

    this->minmax_opt( pd, err );MSQ_ERRRTN( err );
}
void MBMesquite::NonSmoothDescent::print_active_set ( const ActiveSet active_set,
const std::vector< double > &  func 
) [private]

Definition at line 1644 of file NonSmoothDescent.cpp.

References MBMesquite::NonSmoothDescent::ActiveSet::active_ind, MSQ_DBGOUT, and MSQ_PRINT.

Referenced by minmax_opt(), and search_direction().

{
    if( active_set.active_ind.empty() ) MSQ_DBGOUT( 3 ) << "No active values\n";
    /* print the active set */
    for( size_t i = 0; i < active_set.active_ind.size(); i++ )
    {
        MSQ_PRINT( 3 )
        ( "Active value %lu:   %f \n", (unsigned long)i + 1, func[active_set.active_ind[i]] );
    }
}
void MBMesquite::NonSmoothDescent::search_direction ( PatchData pd,
Vector3D search_dir_out,
MsqError err 
) [private]

Definition at line 445 of file NonSmoothDescent.cpp.

References MBMesquite::NonSmoothDescent::ActiveSet::active_ind, b, MBMesquite::NonSmoothDescent::SymmetricMatrix::condition3x3(), moab::E, MBMesquite::EPSILON, form_grammian(), form_PD_grammian(), form_reduced_matrix(), get_active_directions(), MBMesquite::MsqError::INVALID_STATE, mActive, mFunction, mG, mGradient, MSQ_ERRRTN, MSQ_PRINT, MSQ_SETERR, MSQ_ZERO_SEARCH, mSteepest, optStatus, MBMesquite::P, print_active_set(), search_edges_faces(), and solve2x2().

Referenced by minmax_opt().

{
    bool viable;
    double a, b, c, denom;
    std::vector< Vector3D > dir;
    double R0, R1;
    SymmetricMatrix P;
    double x[2];
    double search_mag;

    const int num_active = mActive.active_ind.size();
    ;

    // TODO This might be o.k. actually - i don't see any dependence
    // on the element geometry here... try it and see if it works.
    // if not, try taking all of the gradients in the active set
    // and let the search direction be the average of those.
    //   MSQ_FUNCTION_TIMER( "Search Direction" );

    MSQ_PRINT( 2 )( "\nIn Search Direction\n" );
    this->print_active_set( mActive, mFunction );

    if( num_active == 0 )
    {
        MSQ_SETERR( err )( "No active values in search", MsqError::INVALID_STATE );
        return;
    }

    switch( num_active )
    {
        case 1:
            mSearch   = mGradient[mActive.active_ind[0]];
            mSteepest = mActive.active_ind[0];
            break;
        case 2:
            /* if there are two active points, move in the direction of the
           intersection of the planes.  This is the steepest descent
               direction found by analytically solving the QP */

            /* set up the active gradient directions */
            this->get_active_directions( mGradient, dir, err );MSQ_ERRRTN( err );

            /* form the grammian */
            this->form_grammian( dir, err );MSQ_ERRRTN( err );
            this->form_PD_grammian( err );MSQ_ERRRTN( err );

            denom  = ( mG( 0, 0 ) + mG( 1, 1 ) - 2 * mG( 0, 1 ) );
            viable = true;
            if( fabs( denom ) > EPSILON )
            {
                /* gradients are LI, move along their intersection */
                b = ( mG( 0, 0 ) - mG( 0, 1 ) ) / denom;
                a = 1 - b;
                if( ( b < 0 ) || ( b > 1 ) ) viable = false; /* 0 < b < 1 */
                if( viable )
                {
                    mSearch = a * dir[0] + b * dir[1];
                }
                else
                {
                    /* the gradients are dependent, move along one face */
                    mSearch = dir[0];
                }
            }
            else
            {
                /* the gradients are dependent, move along one face */
                mSearch = dir[0];
            }
            mSteepest = mActive.active_ind[0];

            break;
        default:
            /* as in case 2: solve the QP problem to find the steepest
               descent direction.  This can be done analytically - as
               is done in Gill, Murray and Wright
                 for 3 active points in 3 directions - test PD of G
                 otherwise we know it's SP SD so search edges and faces */

            /* get the active gradient directions */
            this->get_active_directions( mGradient, dir, err );MSQ_ERRRTN( err );

            /* form the entries of the grammian matrix */
            this->form_grammian( dir, err );MSQ_ERRRTN( err );
            this->form_PD_grammian( err );MSQ_ERRRTN( err );

            if( num_active == 3 )
            {
                if( mG.condition3x3() < 1e14 )
                {  // if not singular
                    /* form the entries of P=Z^T G Z where Z = [-1...-1; I ] */
                    this->form_reduced_matrix( P );
                    /* form  the RHS and solve the system for the coeffs */
                    R0      = mG( 0, 0 ) - mG( 1, 0 );
                    R1      = mG( 0, 0 ) - mG( 2, 0 );
                    bool ok = this->solve2x2( P( 0, 0 ), P( 0, 1 ), P( 1, 0 ), P( 1, 1 ), R0, R1, x, err );MSQ_ERRRTN( err );
                    if( ok )
                    {
                        a         = 1 - x[0] - x[1];
                        b         = x[0];
                        c         = x[1];
                        mSearch   = a * dir[0] + b * dir[1] + c * dir[2];
                        mSteepest = mActive.active_ind[0];
                    }
                    else
                    {
                        this->search_edges_faces( &dir[0], mSearch, err );MSQ_ERRRTN( err );
                    }
                }
                else
                {
                    this->search_edges_faces( &dir[0], mSearch, err );MSQ_ERRRTN( err );
                }
            }
            else
            {
                this->search_edges_faces( &dir[0], mSearch, err );MSQ_ERRRTN( err );
            }
    }

    /* if the search direction is essentially zero, equilibrium pt */
    search_mag = mSearch % mSearch;
    MSQ_PRINT( 3 )( "  Search Magnitude %g \n", search_mag );

    if( fabs( search_mag ) < 1E-13 )
        optStatus = MSQ_ZERO_SEARCH;
    else
        mSearch /= std::sqrt( search_mag );
    MSQ_PRINT( 3 )
    ( "  Search Direction %g %g  Steepest %lu\n", mSearch[0], mSearch[1], (unsigned long)mSteepest );
}
void MBMesquite::NonSmoothDescent::search_edges_faces ( const Vector3D dir,
Vector3D search,
MsqError err 
) [private]

Definition at line 1426 of file NonSmoothDescent.cpp.

References MBMesquite::NonSmoothDescent::ActiveSet::active_ind, b, MBMesquite::EPSILON, mActive, mG, MSQ_EQUILIBRIUM, mSteepest, and optStatus.

Referenced by search_direction().

{
    bool viable;
    double a, b, denom;
    Vector3D g_bar;
    Vector3D temp_search( 0, 0, 0 ); /* initialize the search direction to 0,0 */
    double projection, min_projection;

    const size_t num_active = mActive.active_ind.size();
    ;

    /* Check for viable faces */
    min_projection = HUGE_VAL;
    for( size_t i = 0; i < num_active; i++ )
    {
        /* FACE I */
        viable = true;

        /* test the viability */
        for( size_t j = 0; j < num_active; j++ )
        { /* lagrange multipliers>0 */
            if( mG( j, i ) < 0 ) viable = false;
        }

        /* find the minimum of viable directions */
        if( ( viable ) && ( mG( i, i ) < min_projection ) )
        {
            min_projection = mG( i, i );
            temp_search    = dir[i];
            mSteepest      = mActive.active_ind[i];
        }

        /* INTERSECTION IJ */
        for( size_t j = i + 1; j < num_active; j++ )
        {
            viable = true;

            /* find the coefficients of the intersection
               and test the viability */
            denom = 2 * mG( i, j ) - mG( i, i ) - mG( j, j );
            a = b = 0;
            if( fabs( denom ) > EPSILON )
            {
                b = ( mG( i, j ) - mG( i, i ) ) / denom;
                a = 1 - b;
                if( ( b < 0 ) || ( b > 1 ) ) /* 0 < b < 1 */
                    viable = false;
                for( size_t k = 0; k < num_active; k++ )
                { /* lagrange multipliers>0 */
                    if( ( a * mG( k, i ) + b * mG( k, j ) ) <= 0 ) viable = false;
                }
            }
            else
            {
                viable = false; /* Linearly dependent */
            }

            /* find the minimum of viable directions */
            if( viable )
            {
                g_bar      = a * dir[i] + b * dir[j];
                projection = g_bar % g_bar;
                if( projection < min_projection )
                {
                    min_projection = projection;
                    temp_search    = g_bar;
                    mSteepest      = mActive.active_ind[i];
                }
            }
        }
    }
    if( optStatus != MSQ_EQUILIBRIUM )
    {
        mSearch = temp_search;
    }
}
void MBMesquite::NonSmoothDescent::singular_test ( int  n,
const Matrix3D A,
bool &  singular,
MsqError err 
) [private]

Definition at line 1349 of file NonSmoothDescent.cpp.

References MBMesquite::condition3x3(), MBMesquite::EPSILON, MBMesquite::MsqError::INVALID_ARG, and MSQ_SETERR.

Referenced by form_PD_grammian().

{
    //    int test;
    //    double determinant;
    double cond;

    if( ( n > 3 ) || ( n < 1 ) )
    {
        MSQ_SETERR( err )( "Singular test works only for n=1 to n=3", MsqError::INVALID_ARG );
        return;
    }

    singular = true;
    switch( n )
    {
        case 1:
            if( A[0][0] > 0 ) singular = false;
            break;
        case 2:
            if( fabs( A[0][0] * A[1][1] - A[0][1] * A[1][0] ) > EPSILON ) singular = false;
            break;
        case 3:
            /* calculate the condition number */
            cond = condition3x3( A[0] );
            if( cond < 1E14 ) singular = false;
            break;
    }
}
bool MBMesquite::NonSmoothDescent::solve2x2 ( double  a11,
double  a12,
double  a21,
double  a22,
double  b1,
double  b2,
double  x[2],
MsqError err 
) [private]

Definition at line 1503 of file NonSmoothDescent.cpp.

References b1, b2, and MBMesquite::EPSILON.

Referenced by search_direction().

{
    double factor;

    /* if the system is not singular, solve it */
    if( fabs( a11 * a22 - a21 * a12 ) > EPSILON )
    {
        if( fabs( a11 ) > EPSILON )
        {
            factor = ( a21 / a11 );
            x[1]   = ( b2 - factor * b1 ) / ( a22 - factor * a12 );
            x[0]   = ( b1 - a12 * x[1] ) / a11;
        }
        else if( fabs( a21 ) > EPSILON )
        {
            factor = ( a11 / a21 );
            x[1]   = ( b1 - factor * b2 ) / ( a12 - factor * a22 );
            x[0]   = ( b2 - a22 * x[1] ) / a21;
        }
        return true;
    }
    else
    {
        return false;
    }
}
void MBMesquite::NonSmoothDescent::step_acceptance ( PatchData pd,
int  iter_count,
const Vector3D search_dir,
MsqError err 
) [private]

Definition at line 710 of file NonSmoothDescent.cpp.

References compute_function(), find_active_set(), freeVertexIndex, get_min_estimate(), MBMesquite::PatchData::get_vertex_array(), improvement_check(), mActive, mAlpha, mFunction, minStepSize, MBMesquite::PatchData::move_vertex(), MSQ_ERRRTN, MSQ_FLAT_NO_IMP, MSQ_IMP_TOO_SMALL, MSQ_NO_STATUS, MSQ_PRINT, MSQ_STEP_ACCEPTED, MSQ_STEP_TOO_SMALL, optStatus, originalActive, originalFunction, prevActiveValues, MBMesquite::PatchData::set_vertex_coordinates(), testActive, testFunction, MBMesquite::NonSmoothDescent::ActiveSet::true_active_value, and validity_check().

Referenced by minmax_opt().

{
    const double minAcceptableImprovement = 1e-6;

    //  int        ierr;
    //  int        num_values;
    bool step_done;
    bool valid = true, accept_alpha;
    double estimated_improvement;
    double current_improvement  = HUGE_VAL;
    double previous_improvement = HUGE_VAL;
    double current_percent_diff = HUGE_VAL;
    Vector3D original_point;

    //  MSQ_FUNCTION_TIMER( "Step Acceptance" );
    //  num_values = qmHandles.size();

    optStatus = MSQ_NO_STATUS;

    if( mAlpha < minStepSize )
    {
        optStatus = MSQ_IMP_TOO_SMALL;
        step_done = true;
        MSQ_PRINT( 3 )( "Alpha starts too small, no improvement\n" );
    }

    const MsqVertex* coords = pd.get_vertex_array( err );

    /* save the original function and active set */
    original_point   = coords[freeVertexIndex];
    originalFunction = mFunction;
    originalActive   = mActive;

    step_done = false;
    for( int num_steps = 0; !step_done && num_steps < 100; ++num_steps )
    {
        accept_alpha = false;

        while( !accept_alpha && mAlpha > minStepSize )
        {

            /* make the step */
            pd.move_vertex( -mAlpha * mSearch, freeVertexIndex, err );
            // pd.set_coords_array_element(coords[freeVertexIndex],0,err);

            MSQ_PRINT( 2 )( "search direction %f %f \n", mSearch[0], mSearch[1] );
            MSQ_PRINT( 2 )
            ( "new vertex position %f %f \n", coords[freeVertexIndex][0], coords[freeVertexIndex][1] );

            /* assume alpha is acceptable */
            accept_alpha = true;

            /* never take a step that makes a valid mesh invalid or worsens the quality */
            // TODO Validity check revision -- do the compute function up here
            // and then the rest based on validity
            valid = validity_check( pd, err );MSQ_ERRRTN( err );
            if( valid )
            {
                valid = improvement_check( err );MSQ_ERRRTN( err );
            }
            if( !valid )
            {
                accept_alpha = false;
                pd.move_vertex( mAlpha * mSearch, freeVertexIndex, err );
                // pd.set_coords_array_element(coords[freeVertexIndex],0,err);
                mAlpha = mAlpha / 2;
                MSQ_PRINT( 2 )( "Step not accepted, the new alpha %f\n", mAlpha );

                if( mAlpha < minStepSize )
                {
                    optStatus = MSQ_STEP_TOO_SMALL;
                    step_done = true;
                    MSQ_PRINT( 2 )( "Step too small\n" );
                    /* get back the original point, mFunction, and active set */
                    pd.set_vertex_coordinates( original_point, freeVertexIndex, err );
                    mFunction = originalFunction;
                    mActive   = originalActive;
                }
            }
        }

        if( valid && ( mAlpha > minStepSize ) )
        {
            /* compute the new function and active set */
            this->compute_function( &pd, mFunction, err );MSQ_ERRRTN( err );
            this->find_active_set( mFunction, mActive );

            /* estimate the minimum improvement by taking this step */
            this->get_min_estimate( &estimated_improvement, err );MSQ_ERRRTN( err );
            MSQ_PRINT( 2 )
            ( "The estimated improvement for this step: %f\n", estimated_improvement );

            /* calculate the actual increase */
            current_improvement = mActive.true_active_value - prevActiveValues[iterCount - 1];

            MSQ_PRINT( 3 )( "Actual improvement %f\n", current_improvement );

            /* calculate the percent difference from estimated increase */
            current_percent_diff = fabs( current_improvement - estimated_improvement ) / fabs( estimated_improvement );

            /* determine whether to accept a step */
            if( ( fabs( previous_improvement ) > fabs( current_improvement ) ) && ( previous_improvement < 0 ) )
            {
                /* accept the previous step - it was better */
                MSQ_PRINT( 2 )( "Accepting the previous step\n" );

                /* subtract alpha in again (previous step) */
                pd.move_vertex( -mAlpha * mSearch, freeVertexIndex, err );
                // pd.set_coords_array_element(coords[freeVertexIndex],0,err);

                /* does this make an invalid mesh valid? */
                // TODO Validity check revisison
                valid = validity_check( pd, err );MSQ_ERRRTN( err );
                if( valid ) valid = improvement_check( err );MSQ_ERRRTN( err );

                /* copy test function and active set */
                mFunction = testFunction;
                mActive   = testActive;

                optStatus = MSQ_STEP_ACCEPTED;
                step_done = true;

                /* check to see that we're still making good improvements */
                if( fabs( previous_improvement ) < minAcceptableImprovement )
                {
                    optStatus = MSQ_IMP_TOO_SMALL;
                    step_done = true;
                    MSQ_PRINT( 2 )( "Optimization Exiting: Improvement too small\n" );
                }
            }
            else if( ( ( fabs( current_improvement ) > fabs( estimated_improvement ) ) ||
                       ( current_percent_diff < .1 ) ) &&
                     ( current_improvement < 0 ) )
            {
                /* accept this step, exceeded estimated increase or was close */
                optStatus = MSQ_STEP_ACCEPTED;
                step_done = true;

                /* check to see that we're still making good improvements */
                if( fabs( current_improvement ) < minAcceptableImprovement )
                {
                    MSQ_PRINT( 2 )( "Optimization Exiting: Improvement too small\n" );
                    optStatus = MSQ_IMP_TOO_SMALL;
                    step_done = true;
                }
            }
            else if( ( current_improvement > 0 ) && ( previous_improvement > 0 ) &&
                     ( fabs( current_improvement ) < minAcceptableImprovement ) &&
                     ( fabs( previous_improvement ) < minAcceptableImprovement ) )
            {

                /* we are making no progress, quit */
                optStatus = MSQ_FLAT_NO_IMP;
                step_done = true;
                MSQ_PRINT( 2 )( "Opimization Exiting: Flat no improvement\n" );

                /* get back the original point, function, and active set */
                pd.set_vertex_coordinates( original_point, freeVertexIndex, err );MSQ_ERRRTN( err );
                mFunction = originalFunction;
                mActive   = originalActive;
            }
            else
            {
                /* halve alpha and try again */
                /* add out the old step */
                pd.move_vertex( mAlpha * mSearch, freeVertexIndex, err );
                // pd.set_coords_array_element(coords[freeVertexIndex],0,err);

                /* halve step size */
                mAlpha = mAlpha / 2;
                MSQ_PRINT( 3 )( "Step not accepted, the new alpha %f\n", mAlpha );

                if( mAlpha < minStepSize )
                {
                    /* get back the original point, function, and active set */
                    MSQ_PRINT( 2 )( "Optimization Exiting: Step too small\n" );
                    pd.set_vertex_coordinates( original_point, freeVertexIndex, err );MSQ_ERRRTN( err );
                    mFunction = originalFunction;
                    mActive   = originalActive;
                    optStatus = MSQ_STEP_TOO_SMALL;
                    step_done = true;
                }
                else
                {
                    testFunction         = mFunction;
                    testActive           = mActive;
                    previous_improvement = current_improvement;
                }
            }
        }
    }
    if( current_improvement > 0 && optStatus == MSQ_STEP_ACCEPTED )
    {
        MSQ_PRINT( 2 )( "Accepted a negative step %f \n", current_improvement );
    }
}
void MBMesquite::NonSmoothDescent::terminate_mesh_iteration ( PatchData pd,
MsqError err 
) [protected, virtual]

Implements MBMesquite::VertexMover.

Definition at line 132 of file NonSmoothDescent.cpp.

{}
bool MBMesquite::NonSmoothDescent::validity_check ( PatchData pd,
MsqError err 
) [private]

Definition at line 1027 of file NonSmoothDescent.cpp.

References conn, MBMesquite::PatchData::element_by_index(), MBMesquite::EPSILON, MBMesquite::Vector3D::get_coordinates(), MBMesquite::PatchData::get_vertex_array(), MBMesquite::MsqMeshEntity::get_vertex_index_array(), and MBMesquite::PatchData::num_elements().

Referenced by minmax_opt(), optimize_vertex_positions(), and step_acceptance().

{
    //  FUNCTION_TIMER_START("Validity Check");

    // ONLY FOR SIMPLICIAL MESHES - THERE SHOULD BE A VALIDITY CHECKER ASSOCIATED
    // WITH MSQ ELEMENTS

    /* check that the simplicial mesh is still valid, based on right handedness.
         Returns a 1 or a 0 */

    // TODO as a first step we can switch this over to the function
    // evaluation and use the rest of the code as is
    bool valid  = true;
    double dEps = 1.e-13;

    double x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4;
    const MsqVertex* coords = pd.get_vertex_array( err );

    for( size_t i = 0; i < pd.num_elements(); i++ )
    {
        const size_t* conn = pd.element_by_index( i ).get_vertex_index_array();
        coords[conn[0]].get_coordinates( x1, y1, z1 );
        coords[conn[1]].get_coordinates( x2, y2, z2 );
        coords[conn[2]].get_coordinates( x3, y3, z3 );
        coords[conn[3]].get_coordinates( x4, y4, z4 );

        double dDX2 = x2 - x1;
        double dDX3 = x3 - x1;
        double dDX4 = x4 - x1;

        double dDY2 = y2 - y1;
        double dDY3 = y3 - y1;
        double dDY4 = y4 - y1;

        double dDZ2 = z2 - z1;
        double dDZ3 = z3 - z1;
        double dDZ4 = z4 - z1;

        /* dDet is proportional to the cell volume */
        double dDet = dDX2 * dDY3 * dDZ4 + dDX3 * dDY4 * dDZ2 + dDX4 * dDY2 * dDZ3 - dDZ2 * dDY3 * dDX4 -
                      dDZ3 * dDY4 * dDX2 - dDZ4 * dDY2 * dDX3;

        /* Compute a length scale based on edge lengths. */
        double dScale = ( sqrt( ( x1 - x2 ) * ( x1 - x2 ) + ( y1 - y2 ) * ( y1 - y2 ) + ( z1 - z2 ) * ( z1 - z2 ) ) +
                          sqrt( ( x1 - x3 ) * ( x1 - x3 ) + ( y1 - y3 ) * ( y1 - y3 ) + ( z1 - z3 ) * ( z1 - z3 ) ) +
                          sqrt( ( x1 - x4 ) * ( x1 - x4 ) + ( y1 - y4 ) * ( y1 - y4 ) + ( z1 - z4 ) * ( z1 - z4 ) ) +
                          sqrt( ( x2 - x3 ) * ( x2 - x3 ) + ( y2 - y3 ) * ( y2 - y3 ) + ( z2 - z3 ) * ( z2 - z3 ) ) +
                          sqrt( ( x2 - x4 ) * ( x2 - x4 ) + ( y2 - y4 ) * ( y2 - y4 ) + ( z2 - z4 ) * ( z2 - z4 ) ) +
                          sqrt( ( x3 - x4 ) * ( x3 - x4 ) + ( y3 - y4 ) * ( y3 - y4 ) + ( z3 - z4 ) * ( z3 - z4 ) ) ) /
                        6.;

        /* Use the length scale to get a better idea if the tet is flat or
           just really small. */
        if( fabs( dScale ) < EPSILON )
        {
            return false;
        }
        else
        {
            dDet /= ( dScale * dScale * dScale );
        }

        if( dDet > dEps )
        {
            valid = true;
        }
        else if( dDet < -dEps )
        {
            return false;
        }
        else
        {
            return false;
        }
    }  // end for i=1,numElements

    //  MSQ_DEBUG_ACTION(2,{fprintf(stdout,"Mesh Validity is: %d \n",valid);});

    //  FUNCTION_TIMER_END();
    return ( valid );
}

Member Data Documentation

Definition at line 193 of file NonSmoothDescent.hpp.

Referenced by compute_alpha(), init_max_step_length(), and init_opt().

std::vector< double > MBMesquite::NonSmoothDescent::mFunction [private]
std::vector< double > MBMesquite::NonSmoothDescent::mGS [private]

Definition at line 177 of file NonSmoothDescent.hpp.

Referenced by compute_alpha(), initialize(), and step_acceptance().

Definition at line 186 of file NonSmoothDescent.hpp.

Referenced by form_PD_grammian(), and init_opt().

std::vector< double > MBMesquite::NonSmoothDescent::originalFunction [private]

Definition at line 182 of file NonSmoothDescent.hpp.

Referenced by improvement_check(), and minmax_opt().

Definition at line 194 of file NonSmoothDescent.hpp.

Referenced by form_PD_grammian(), and init_opt().

std::vector< double > MBMesquite::NonSmoothDescent::prevActiveValues [private]

Definition at line 187 of file NonSmoothDescent.hpp.

Referenced by init_opt(), minmax_opt(), and step_acceptance().

std::vector< double > MBMesquite::NonSmoothDescent::testFunction [private]

Definition at line 183 of file NonSmoothDescent.hpp.

Referenced by init_opt(), and step_acceptance().

Definition at line 188 of file NonSmoothDescent.hpp.

Referenced by compute_gradient().

std::vector< size_t > MBMesquite::NonSmoothDescent::tmpIdx [private]

Definition at line 189 of file NonSmoothDescent.hpp.

Referenced by compute_gradient().

List of all members.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines