#include <NonSmoothSteepestDescent.hpp>

Inheritance diagram for NonSmoothSteepestDescent:

Collaboration diagram for NonSmoothSteepestDescent:

Public Member Functions
	NonSmoothSteepestDescent (ObjectiveFunction *of)

virtual	~NonSmoothSteepestDescent ()

	NonSmoothSteepestDescent (ObjectiveFunction *of)

virtual	~NonSmoothSteepestDescent ()

Public Member Functions inherited from VertexMover
virtual	~VertexMover ()

virtual double	loop_over_mesh (MeshSet &ms, MsqError &err)
	Improves the quality of the MeshSet, calling some methods specified in a class derived from VertexMover. More...

virtual	~VertexMover ()

virtual double	loop_over_mesh (MeshSet &ms, MsqError &err)
	This is the "run" function of PatchDataUser. It can do anything really. More...

Public Member Functions inherited from QualityImprover
virtual	~QualityImprover ()

void	set_name (msq_std::string name)
	provides a name to the QualityImprover (use it in constructor). More...

virtual msq_std::string	get_name ()
	retrieves the QualityImprover name. A default name should be set in the constructor. More...

virtual AlgorithmType	get_algorithm_type ()
	Return the algorithm type (to avoid RTTI use). More...

void	set_inner_termination_criterion (TerminationCriterion *crit)
	Sets in the termination criterion for the concrete solver's optimization. More...

void	set_outer_termination_criterion (TerminationCriterion *crit)
	Sets in the termination criterion for the outer loop over patches. More...

virtual	~QualityImprover ()

void	set_name (msq_std::string name)
	provides a name to the QualityImprover (use it in constructor). More...

virtual msq_std::string	get_name ()
	retrieves the QualityImprover name. A default name should be set in the constructor. More...

virtual AlgorithmType	get_algorithm_type ()
	Return the algorithm type (to avoid RTTI use). More...

void	set_inner_termination_criterion (TerminationCriterion *crit)
	Sets in the termination criterion for the concrete solver's optimization. More...

void	set_outer_termination_criterion (TerminationCriterion *crit)
	Sets in the termination criterion for the outer loop over patches. More...

Public Member Functions inherited from PatchDataUser
virtual	~PatchDataUser ()

virtual void	set_patch_type (PatchData::PatchType patch_type, MsqError &err, int param1=0, int param2=0)
	Sets the Patch Type. More...

PatchData::PatchType	get_patch_type ()
	Returns the Patch Type. More...

int	get_nb_layers (MsqError &err)
	Returns number of layers (if relevant for partition algorythm). More...

void	add_culling_method (enum PatchData::culling_method cm)
	Sets on the culling method passed as argument. More...

void	no_culling_method ()
	Sets off all culling methods. More...

void	remove_culling_method (enum PatchData::culling_method cm)
	Sets off the culling method passed as argument. More...

long unsigned int	get_culling_method_bits ()
	Returns the bitset containing culling methods flags. More...

void	set_all_parameters (PatchDataParameters &params)

PatchDataParameters &	get_all_parameters ()
	Returns the PatchDataParameters object. More...

void	set_global_patch (PatchData *pd, MsqError &err)
	Sets the Global Patch, so that it can be use by contiguoug PatchDataUser. More...

PatchData *	get_global_patch ()
	Returns the Global Patch. More...

void	no_global_patch ()
	Sets the Global Patch pointer to NULL. More...

virtual	~PatchDataUser ()

virtual void	set_patch_type (PatchData::PatchType patch_type, MsqError &err, int param1=0, int param2=0)
	Sets the Patch Type. More...

PatchData::PatchType	get_patch_type ()
	Returns the Patch Type. More...

int	get_nb_layers (MsqError &err)
	Returns number of layers (if relevant for partition algorythm). More...

void	add_culling_method (enum PatchData::culling_method cm)
	Sets on the culling method passed as argument. More...

void	no_culling_method ()
	Sets off all culling methods. More...

void	remove_culling_method (enum PatchData::culling_method cm)
	Sets off the culling method passed as argument. More...

long unsigned int	get_culling_method_bits ()
	Returns the bitset containing culling methods flags. More...

void	set_all_parameters (PatchDataParameters &params)

PatchDataParameters &	get_all_parameters ()
	Returns the PatchDataParameters object. More...

void	set_global_patch (PatchData *pd, MsqError &err)
	Sets the Global Patch, so that it can be use by contiguoug PatchDataUser. More...

PatchData *	get_global_patch ()
	Returns the Global Patch. More...

void	no_global_patch ()
	Sets the Global Patch pointer to NULL. More...

Protected Member Functions
virtual void	initialize (PatchData &pd, MsqError &err)

virtual void	optimize_vertex_positions (PatchData &pd, MsqError &err)

virtual void	initialize_mesh_iteration (PatchData &pd, MsqError &err)

virtual void	terminate_mesh_iteration (PatchData &pd, MsqError &err)

virtual void	cleanup ()

virtual void	initialize (PatchData &pd, MsqError &err)

virtual void	optimize_vertex_positions (PatchData &pd, MsqError &err)

virtual void	initialize_mesh_iteration (PatchData &pd, MsqError &err)

virtual void	terminate_mesh_iteration (PatchData &pd, MsqError &err)

virtual void	cleanup ()

Protected Member Functions inherited from VertexMover
	VertexMover ()

size_t	check_feasible (PatchData &pd, MsqError &err)
	CHECK FEASIBLE IS NOT YET IMPLEMENTED. More...

	VertexMover ()

size_t	check_feasible (PatchData &pd, MsqError &err)
	CHECK FEASIBLE IS NOT YET IMPLEMENTED. More...

Protected Member Functions inherited from QualityImprover
	QualityImprover ()

const MeshSet *	get_mesh_set () const

MeshSet *	get_mesh_set ()

void	set_mesh_set (MeshSet *ms)

TerminationCriterion *	get_outer_termination_criterion ()
	return the outer termination criterion pointer More...

TerminationCriterion *	get_inner_termination_criterion ()
	return the inner termination criterion pointer More...

	QualityImprover ()

const MeshSet *	get_mesh_set () const

MeshSet *	get_mesh_set ()

void	set_mesh_set (MeshSet *ms)

TerminationCriterion *	get_outer_termination_criterion ()
	return the outer termination criterion pointer More...

TerminationCriterion *	get_inner_termination_criterion ()
	return the inner termination criterion pointer More...

Protected Member Functions inherited from PatchDataUser
	PatchDataUser ()

	PatchDataUser ()

Private Member Functions
void	init_opt (MsqError &err)

void	init_max_step_length (MsqError &err)

void	minmax_opt (PatchData &pd, MsqError &err)

void	step_acceptance (PatchData &pd, MsqError &err)

void	get_min_estimate (double *final_est, MsqError &err)

void	get_gradient_projections (MsqError &err)

void	compute_alpha (MsqError &err)

void	copy_active (ActiveSet active1, ActiveSet active2, MsqError &err)

bool	compute_function (PatchData pd, double function, MsqError &err)

double **	compute_gradient (PatchData *pd, MsqError &err)

void	find_active_set (double function, ActiveSet active_set, MsqError &err)

void	print_active_set (ActiveSet active_set, double func, MsqError &err)

int	improvement_check (MsqError &err)

int	validity_check (MsqError &err)

void	check_equilibrium (int equil, int opt_status, MsqError &err)

int	convex_hull_test (double **vec, int num_vec, MsqError &err)

int	check_vector_dots (double *vec, int num_vec, double normal, MsqError &err)

void	find_plane_normal (double pt1[3], double pt2[3], double pt3[3], double *cross, MsqError &err)

void	find_plane_points (int dir1, int dir2, double *vec, int num_vec, double pt1, double pt2, double pt3, int *opt_status, MsqError &err)

void	form_grammian (double **vec, MsqError &err)

void	form_PD_grammian (MsqError &err)

void	singular_test (int n, double *A, int singular, MsqError &err)

void	condition3x3 (double *A, double cond, MsqError &err)

void	solve2x2 (double a11, double a12, double a21, double a22, double b1, double b2, double **x, MsqError &err)

void	form_reduced_matrix (double ***P, MsqError &err)

void	search_direction (PatchData &pd, MsqError &err)

void	search_edges_faces (double **dir, MsqError &err)

void	get_active_directions (double gradient, double *dir, MsqError &err)

void	init_opt (MsqError &err)

void	init_max_step_length (MsqError &err)

void	minmax_opt (PatchData &pd, MsqError &err)

void	step_acceptance (PatchData &pd, MsqError &err)

void	get_min_estimate (double *final_est, MsqError &err)

void	get_gradient_projections (MsqError &err)

void	compute_alpha (MsqError &err)

void	copy_active (ActiveSet active1, ActiveSet active2, MsqError &err)

bool	compute_function (PatchData pd, double function, MsqError &err)

double **	compute_gradient (PatchData *pd, MsqError &err)

void	find_active_set (double function, ActiveSet active_set, MsqError &err)

void	print_active_set (ActiveSet active_set, double func, MsqError &err)

int	improvement_check (MsqError &err)

int	validity_check (MsqError &err)

void	check_equilibrium (int equil, int opt_status, MsqError &err)

int	convex_hull_test (double **vec, int num_vec, MsqError &err)

int	check_vector_dots (double *vec, int num_vec, double normal, MsqError &err)

void	find_plane_normal (double pt1[3], double pt2[3], double pt3[3], double *cross, MsqError &err)

void	find_plane_points (int dir1, int dir2, double *vec, int num_vec, double pt1, double pt2, double pt3, int *opt_status, MsqError &err)

void	form_grammian (double **vec, MsqError &err)

void	form_PD_grammian (MsqError &err)

void	singular_test (int n, double *A, int singular, MsqError &err)

void	condition3x3 (double *A, double cond, MsqError &err)

void	solve2x2 (double a11, double a12, double a21, double a22, double b1, double b2, double **x, MsqError &err)

void	form_reduced_matrix (double ***P, MsqError &err)

void	search_direction (PatchData &pd, MsqError &err)

void	search_edges_faces (double **dir, MsqError &err)

void	get_active_directions (double gradient, double *dir, MsqError &err)

Private Attributes
int	mDimension

int	numVertices

int	numElements

MsqVertex *	mCoords

MsqMeshEntity *	mConnectivity

int	numFree

int	freeVertexIndex

double	activeEpsilon

double	minAcceptableImprovement

double	minStepSize

double	originalValue

int	iterCount

int	optIterCount

int	numFunctionValues

double *	mFunction

double *	testFunction

double *	originalFunction

double **	mGradient

int	optStatus

int	equilibriumPt

int	mSteepest

double	mSearch [3]

double	mAlpha

double	maxAlpha

double *	mGS

double *	prevActiveValues

double **	mG

double **	mPDG

int	pdgInd [3]

ActiveSet *	mActive

ActiveSet *	testActive

ActiveSet *	originalActive

Additional Inherited Members
Public Types inherited from PatchDataUser
enum	AlgorithmType { QUALITY_IMPROVER, QUALITY_ASSESSOR, MESH_TRANSFORM, TARGET_CALCULATOR, QUALITY_IMPROVER, QUALITY_ASSESSOR, MESH_TRANSFORM, TARGET_CALCULATOR }

enum	AlgorithmType { QUALITY_IMPROVER, QUALITY_ASSESSOR, MESH_TRANSFORM, TARGET_CALCULATOR, QUALITY_IMPROVER, QUALITY_ASSESSOR, MESH_TRANSFORM, TARGET_CALCULATOR }

Protected Attributes inherited from VertexMover
ObjectiveFunction *	objFunc

Detailed Description

Definition at line 145 of file includeLinks/NonSmoothSteepestDescent.hpp.

Constructor & Destructor Documentation

NonSmoothSteepestDescent ( ObjectiveFunction * of )

Definition at line 45 of file QualityImprover/VertexMover/NonSmoothSteepestDescent/NonSmoothSteepestDescent.cpp.

References i, NonSmoothSteepestDescent::mActive, NonSmoothSteepestDescent::mFunction, NonSmoothSteepestDescent::mG, NonSmoothSteepestDescent::mGradient, NonSmoothSteepestDescent::mGS, NonSmoothSteepestDescent::mPDG, MSQ_DBGOUT, VertexMover::objFunc, NonSmoothSteepestDescent::originalActive, NonSmoothSteepestDescent::originalFunction, NonSmoothSteepestDescent::prevActiveValues, QualityImprover::set_name(), NonSmoothSteepestDescent::testActive, and NonSmoothSteepestDescent::testFunction.

 {
   objFunc=of;
   this->set_name("NonSmoothSteepestDescent");
   
   mFunction = (double *)malloc(sizeof(double)*150);
   testFunction = (double *)malloc(sizeof(double)*150);
   originalFunction = (double *)malloc(sizeof(double)*150);
   mGradient = (double **)malloc(sizeof(double *)*150);
   mGS = (double *)malloc(sizeof(double)*150);
   mG = (double **)malloc(sizeof(double *)*150);
   mPDG = (double **)malloc(sizeof(double *)*3);
   prevActiveValues=(double *)malloc(sizeof(double)*150);
 
   int i;
   for (i=0;i<150;i++) {
     mGradient[i] = (double *)malloc(sizeof(double)*3);
     mG[i] = (double *)malloc(sizeof(double)*150);
   }
   for (i=0;i<3;i++) mPDG[i] = (double *)malloc(sizeof(double)*3);
 
   mActive = (ActiveSet *)malloc(sizeof(ActiveSet));
   testActive = (ActiveSet *)malloc(sizeof(ActiveSet));
   originalActive = (ActiveSet *)malloc(sizeof(ActiveSet));
   MSQ_DBGOUT(1) << "- Executed NonSmoothSteepestDescent::NonSmoothSteepestDescent()\n";
 }  

Here is the call graph for this function:

virtual ~NonSmoothSteepestDescent ( )

inlinevirtual

Definition at line 150 of file includeLinks/NonSmoothSteepestDescent.hpp.

150 { }

NonSmoothSteepestDescent ( ObjectiveFunction * of )

virtual ~NonSmoothSteepestDescent ( )

inlinevirtual

Definition at line 150 of file src/QualityImprover/VertexMover/NonSmoothSteepestDescent/NonSmoothSteepestDescent.hpp.

150 { }

Member Function Documentation

void check_equilibrium	(	int *	equil,
		int *	opt_status,
		MsqError &	err
	)

inlineprivate

Definition at line 658 of file includeLinks/NonSmoothSteepestDescent.hpp.

References NonSmoothSteepestDescent::convex_hull_test(), NonSmoothSteepestDescent::form_grammian(), NonSmoothSteepestDescent::get_active_directions(), i, j, NonSmoothSteepestDescent::mActive, NonSmoothSteepestDescent::mDimension, NonSmoothSteepestDescent::mG, NonSmoothSteepestDescent::mGradient, min(), MSQ_DOT, MSQ_EQUILIBRIUM, MSQ_FALSE, MSQ_MACHINE_EPS, MSQ_NORMALIZE, MSQ_PRINT, ActiveSet::num_active, and NonSmoothSteepestDescent::numFunctionValues.

Referenced by NonSmoothSteepestDescent::minmax_opt().

 {
 //    int  ierr;
     int  i,j;
     int  ind1, ind2;  
     double min;
     double **dir;
     double mid_vec[3], mid_cos, test_cos;
 
     //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
 
     *equil = MSQ_FALSE;
     ind1 = ind2 = -1;
 
     int num_active = mActive->num_active;
 
     if (num_active==numFunctionValues)
     {
          *equil = 1; *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);
 
     /* normalize the active directions */
     for (j=0;j<num_active;j++) MSQ_NORMALIZE(dir[j],mDimension);
 
     if (mDimension == 2) {
       /* form the grammian */
       this->form_grammian(dir,err);  
 
     /* find the minimum element in the upper triangular portion of G*/
     min = 1;
     for (i=0;i<num_active;i++) {
       for (j=i+1;j<num_active;j++) {
         if ( fabs(-1 - mG[i][j]) < 1E-08 ) {
            *equil = 1; *status = MSQ_EQUILIBRIUM;
            MSQ_PRINT(3)("The gradients are antiparallel, eq. pt\n"); 
          }
          if (mG[i][j]  < min) { 
            ind1 = i; ind2 = j;
            min = mG[i][j];
         }
       }
     }
 
     if ((ind1 != -1) && (ind2 != -1)) {
       /* find the diagonal of the parallelepiped */
       for (j=0;j<mDimension;j++) {
        mid_vec[j]=.5*(dir[ind1][j]+dir[ind2][j]);
       }
       MSQ_NORMALIZE(mid_vec,mDimension);
       MSQ_DOT(mid_cos,dir[ind1],mid_vec,mDimension);
 
       /* test the other vectors to be sure they lie in the cone */
       for (i=0;i<num_active;i++) {
          if ((i != ind1) && (i != ind2)) {
             MSQ_DOT(test_cos,dir[i],mid_vec,mDimension);
             if ((test_cos < mid_cos)  &&  fabs(test_cos-mid_cos) > MSQ_MACHINE_EPS) {
               MSQ_PRINT(3)("An equilibrium point \n");
               *equil = 1; *status = MSQ_EQUILIBRIUM;
             }
          }
        }
      }
     }
     if (mDimension == 3) {
        *equil = this->convex_hull_test(dir,num_active,err);
        if (*equil == 1) *status = MSQ_EQUILIBRIUM;
     }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void check_equilibrium	(	int *	equil,
		int *	opt_status,
		MsqError &	err
	)

private

int check_vector_dots	(	double **	vec,
		int	num_vec,
		double *	normal,
		MsqError &	err
	)

inlineprivate

Definition at line 517 of file includeLinks/NonSmoothSteepestDescent.hpp.

References i, MSQ_DOT, MSQ_FALSE, MSQ_MACHINE_EPS, and MSQ_TRUE.

Referenced by NonSmoothSteepestDescent::convex_hull_test().

 {
     int equil;
     int i, ind;
     double test_dot, dot;
 
     equil = MSQ_FALSE;
     MSQ_DOT(test_dot,vec[0],normal,3);
     ind = 1;
     while ((fabs(test_dot) < MSQ_MACHINE_EPS) && (ind<num_vec)) {
       MSQ_DOT(test_dot,vec[ind],normal,3);
       ind++;
     }
       
     for (i=ind;i<num_vec;i++) {
        MSQ_DOT(dot,vec[i],normal,3);
        if ( ((dot>0 && test_dot<0) || (dot<0 && test_dot>0)) &&
             (fabs(dot)>MSQ_MACHINE_EPS)) {
           return(equil = MSQ_TRUE);
 
        }
     }
     return(equil);
 }

Here is the caller graph for this function:

int check_vector_dots	(	double **	vec,
		int	num_vec,
		double *	normal,
		MsqError &	err
	)

private

virtual void cleanup ( )

protectedvirtual

Implements VertexMover.

void cleanup ( )

protectedvirtual

Implements VertexMover.

Definition at line 169 of file QualityImprover/VertexMover/NonSmoothSteepestDescent/NonSmoothSteepestDescent.cpp.

References i, NonSmoothSteepestDescent::mActive, NonSmoothSteepestDescent::mFunction, NonSmoothSteepestDescent::mG, NonSmoothSteepestDescent::mGradient, NonSmoothSteepestDescent::mGS, NonSmoothSteepestDescent::mPDG, MSQ_DBGOUT, NonSmoothSteepestDescent::originalActive, NonSmoothSteepestDescent::originalFunction, NonSmoothSteepestDescent::prevActiveValues, NonSmoothSteepestDescent::testActive, and NonSmoothSteepestDescent::testFunction.

 {
   MSQ_DBGOUT(1) << "- Executing NonSmoothSteepestDescent::cleanup()\n";
   int i;
   for (i=0;i<150;i++) {
     free(mGradient[i]);
     free(mG[i]);
   }
   for (i=0;i<3;i++) free(mPDG[i]);
 
   free(mFunction);
   free(testFunction);
   free(originalFunction);
   free(mGradient);
   free(mGS);
   free(mG);
   free(mPDG);
   free(prevActiveValues);
   free(mActive);
   free(testActive);
   free(originalActive);
   MSQ_DBGOUT(1) << "- Done with NonSmoothSteepestDescent::cleanup()\n";
 }

void compute_alpha ( MsqError & err )

private

void compute_alpha ( MsqError & err )

inlineprivate

Definition at line 1025 of file includeLinks/NonSmoothSteepestDescent.hpp.

References i, NonSmoothSteepestDescent::mAlpha, NonSmoothSteepestDescent::maxAlpha, NonSmoothSteepestDescent::mFunction, NonSmoothSteepestDescent::mGS, NonSmoothSteepestDescent::minStepSize, MSQ_PRINT, NonSmoothSteepestDescent::mSteepest, and NonSmoothSteepestDescent::numFunctionValues.

Referenced by NonSmoothSteepestDescent::minmax_opt().

 {
 //    int       ierr;
 //    int       j;
     int       i;
 //    int       ind;
     int       num_values;
     double    steepest_function;
     double    steepest_grad;
     double    alpha_i;
     double    min_positive_value=1E300;
 
 //    FUNCTION_TIMER_START("Compute Alpha");
 
     MSQ_PRINT(2)("In compute alpha\n");
 
     num_values = numFunctionValues;
     mAlpha = 1E300;
 
     steepest_function = mFunction[mSteepest];
     steepest_grad = mGS[mSteepest];
     for (i=0;i<num_values;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 %d %g\n",i,alpha_i);
           }
           if ((alpha_i > 0) && (alpha_i < min_positive_value)) {
             min_positive_value = alpha_i;
           }
        }
     }
 
     if ((mAlpha == 1E300) && (min_positive_value != 1E300)) {
       mAlpha = min_positive_value;
     }
 
     /* if it never gets set, set it to the default */
     if (mAlpha == 1E300) {
       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();
 }

Here is the caller graph for this function:

bool compute_function	(	PatchData *	pd,
		double *	function,
		MsqError &	err
	)

private

bool compute_function	(	PatchData *	pd,
		double *	function,
		MsqError &	err
	)

inlineprivate

Definition at line 248 of file includeLinks/NonSmoothSteepestDescent.hpp.

References PatchData::element_by_index(), QualityMetric::evaluate_element(), ObjectiveFunction::get_quality_metric(), ObjectiveFunction::get_quality_metric_list(), i, MSQ_ERRZERO, NonSmoothSteepestDescent::numElements, and VertexMover::objFunc.

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

 {
   // ASSUMES ONE VALUE PER ELEMENT; ALSO NEED 1.0/FUNCTION WHICH IS ONLY
   // TRUE OF CONDITION NUMBER
 
   //  MSQ_DEBUG_PRINT(2,"Computing Function\n");
 //  FUNCTION_TIMER_START("Compute Function");
 
   //TODO need to switch this to element or vertex metric evaluations
   //TODO need to include boolean testing for validity
   int i;
   bool valid_bool=true;
 
   for (i=0;i<numElements;i++) func[i]=0.0;
   QualityMetric* currentQM=objFunc->get_quality_metric();
   if(currentQM==NULL){
     currentQM = objFunc->get_quality_metric_list().front();
   }
   
   for (i=0;i<numElements;i++) {
     valid_bool = currentQM->evaluate_element(*patch_data,
                                 &(patch_data->element_by_index(i)),
                                 func[i], err); MSQ_ERRZERO(err);
     //    MSQ_DEBUG_ACTION(3,{fprintf(stdout,"  Function value[%d]=%g\n",i,func[i]);});
   }
 //  FUNCTION_TIMER_END();
   return(valid_bool);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

double** compute_gradient	(	PatchData *	pd,
		MsqError &	err
	)

private

double ** compute_gradient	(	PatchData *	pd,
		MsqError &	err
	)

inlineprivate

Definition at line 278 of file includeLinks/NonSmoothSteepestDescent.hpp.

References NonSmoothSteepestDescent::compute_function(), NonSmoothSteepestDescent::freeVertexIndex, ObjectiveFunction::get_quality_metric(), ObjectiveFunction::get_quality_metric_list(), i, j, NonSmoothSteepestDescent::mCoords, NonSmoothSteepestDescent::mGradient, MSQ_CHKERR, MSQ_DBGOUT, NonSmoothSteepestDescent::numElements, NonSmoothSteepestDescent::numFunctionValues, and VertexMover::objFunc.

Referenced by NonSmoothSteepestDescent::minmax_opt().

 {
 //  FUNCTION_TIMER_START("Compute Gradient");
 
   MSQ_DBGOUT(2) << "Computing Gradient\n";
 
   double delta = 10e-6;
 
   for (int i=0;i<numElements;i++) {
     for (int j=0;j<3;j++) mGradient[i][j] = 0.0;
   }
   QualityMetric* currentQM=objFunc->get_quality_metric();
   if(currentQM==NULL)
     currentQM = objFunc->get_quality_metric_list().front();
 
   double *func, *fdelta;
   func = (double *)malloc(sizeof(double)*150);
   fdelta = (double *)malloc(sizeof(double)*150);
 
   this->compute_function(patch_data, func, err); 
   if (MSQ_CHKERR(err)) {
     free(func);
     free(fdelta);
     return 0;
   }
 
   /* 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
     mCoords[freeVertexIndex][j] += delta;
 
     //compute the function at the perturbed point location
     this->compute_function(patch_data, fdelta, err);  
     if (MSQ_CHKERR(err)) {
       free(func);
       free(fdelta);
       return 0;
     }
 
     //compute the numerical gradient
     for (int i=0;i<numFunctionValues;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
     mCoords[freeVertexIndex][j] -= delta;
   }
   
   free(func);
   free(fdelta);
 //  FUNCTION_TIMER_END();
   return(mGradient);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void condition3x3	(	double **	A,
		double *	cond,
		MsqError &	err
	)

inlineprivate

Definition at line 736 of file includeLinks/NonSmoothSteepestDescent.hpp.

References denom, MSQ_FALSE, MSQ_MACHINE_EPS, MSQ_TRUE, and sqrt().

Referenced by NonSmoothSteepestDescent::singular_test().

 {
 //   int ierr;
    double a11, a12, a13;
    double a21, a22, a23;
    double a31, a32, a33;
 //   double s1, s2, s4, s3, t0;
    double s1, s2, s3;
    double denom;
 //   double one = 1.0;
    double temp;
    int zero_denom = MSQ_TRUE;
 
    a11 = A[0][0]; a12=A[0][1]; a13=A[0][2];
    a21 = A[1][0]; a22=A[1][1]; a23=A[1][2];
    a31 = A[2][0]; a32=A[2][1]; a33=A[2][2];
 
    denom = -a11*a22*a33+a11*a23*a32+a21*a12*a33-a21*a13*a32-
             a31*a12*a23+a31*a13*a22;
 
    if ( (fabs(a11) > MSQ_MACHINE_EPS) && 
         (fabs(denom/a11) > MSQ_MACHINE_EPS)) {
          zero_denom = MSQ_FALSE;
    }
    if ( (fabs(a22) > MSQ_MACHINE_EPS) && 
         (fabs(denom/a22) > MSQ_MACHINE_EPS)) {
          zero_denom = MSQ_FALSE;
    }       
    if ( (fabs(a33) > MSQ_MACHINE_EPS) && 
         (fabs(denom/a33) > MSQ_MACHINE_EPS)) {
          zero_denom = MSQ_FALSE;
    }
 
    if (zero_denom) {
      (*cond) = 1E300;
    } else {
      s1 = sqrt(a11*a11 + a12*a12 + a13*a13 + 
                a21*a21 + a22*a22 + a23*a23 + 
                a31*a31 + a32*a32 + a33*a33);
 
      temp = (-a22*a33+a23*a32)/denom;
      s3 = temp*temp;
      temp =(a12*a33-a13*a32)/denom;
      s3 += temp*temp;
      temp = (a12*a23-a13*a22)/denom;
      s3 += temp*temp;
      temp = (a21*a33-a23*a31)/denom;
      s3 += temp*temp;
      temp = (a11*a33-a13*a31)/denom;
      s3 += temp*temp;
      temp = (a11*a23-a13*a21)/denom;
      s3 += temp*temp;
      temp = (a21*a32-a22*a31)/denom;
      s3 += temp*temp;
      temp = (-a11*a32+a12*a31)/denom;
      s3 += temp*temp;
      temp = (-a11*a22+a12*a21)/denom;
      s3 += temp*temp;
 
      s2 = sqrt(s3);
      (*cond) = s1*s2;
    }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void condition3x3	(	double **	A,
		double *	cond,
		MsqError &	err
	)

private

int convex_hull_test	(	double **	vec,
		int	num_vec,
		MsqError &	err
	)

inlineprivate

Definition at line 566 of file includeLinks/NonSmoothSteepestDescent.hpp.

References NonSmoothSteepestDescent::check_vector_dots(), NonSmoothSteepestDescent::find_plane_normal(), NonSmoothSteepestDescent::find_plane_points(), 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 NonSmoothSteepestDescent::check_equilibrium().

 {
 //    int ierr;
     int equil;
     int status, dir_done;
     double pt1[3], pt2[3], pt3[3];
     double normal[3];
 
 //    FUNCTION_TIMER_START("Convex Hull Test");
     /* tries to determine equilibrium for the 3D case */
     equil = 0;
     status = MSQ_CHECK_Z_COORD_DIRECTION;
     dir_done = -1;
 
     if (num_vec <= 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, num_vec, pt1, pt2, pt3, &status, err);
           dir_done = 2;
        }
        if (status == MSQ_CHECK_Y_COORD_DIRECTION) {
           this->find_plane_points(MSQ_YDIR, MSQ_XDIR, 
                           vec, num_vec, pt1, pt2, pt3, &status, err);
           dir_done = 1;
        }
        if (status == MSQ_CHECK_X_COORD_DIRECTION) {
           this->find_plane_points(MSQ_XDIR, MSQ_ZDIR, 
                           vec, num_vec, pt1, pt2, pt3, &status, err);
           dir_done = 0;
        }
        if (status == MSQ_TWO_PT_PLANE) {
           pt3[0]=0.; pt3[1]=0.; pt3[2]=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,num_vec,normal,err); 
            if (equil == 1) {
              switch(dir_done){
              case 2:
                equil = 0; status = MSQ_CHECK_Y_COORD_DIRECTION;
                break;
              case 1:
                equil = 0; status = MSQ_CHECK_X_COORD_DIRECTION;
                break;
              case 0:
                equil = 1; 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 = 0; 
       break;
     case MSQ_EQUIL:
       MSQ_PRINT(3)("An equilibrium point\n");
       equil = 1;
       break;
     default:
       MSQ_PRINT(3)("Failed to determine equil or not; status = %d\n",status);
     }
 //    FUNCTION_TIMER_END();
     return (equil);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

int convex_hull_test	(	double **	vec,
		int	num_vec,
		MsqError &	err
	)

private

void copy_active	(	ActiveSet *	active1,
		ActiveSet *	active2,
		MsqError &	err
	)

private

void copy_active	(	ActiveSet *	active1,
		ActiveSet *	active2,
		MsqError &	err
	)

inlineprivate

Definition at line 1086 of file includeLinks/NonSmoothSteepestDescent.hpp.

References ActiveSet::active_ind, i, MsqError::INVALID_ARG, MSQ_SETERR, ActiveSet::num_active, ActiveSet::num_equal, and ActiveSet::true_active_value.

Referenced by NonSmoothSteepestDescent::step_acceptance().

 {
     if (active1==NULL || active2==NULL) {
        MSQ_SETERR(err)("Null memory in copy_active\n",MsqError::INVALID_ARG);
        return;
     }
 
     active2->num_active = active1->num_active;
     active2->num_equal  = active1->num_equal;
     active2->true_active_value = active1->true_active_value;
     for (int i=0;i<active1->num_active;i++) {
         active2->active_ind[i] = active1->active_ind[i];
     }
 }

Here is the caller graph for this function:

void find_active_set	(	double *	function,
		ActiveSet *	active_set,
		MsqError &	err
	)

inlineprivate

Definition at line 335 of file includeLinks/NonSmoothSteepestDescent.hpp.

References ActiveSet::active_ind, NonSmoothSteepestDescent::activeEpsilon, for(), i, MSQ_DBGOUT, MSQ_MACHINE_EPS, Mesquite::MSQ_MAX, ActiveSet::num_active, ActiveSet::num_equal, NonSmoothSteepestDescent::numFunctionValues, and ActiveSet::true_active_value.

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

 { 
     int         i, ind;
     double      function_val;
     double      active_value0;
     double      temp;
 
 //    FUNCTION_TIMER_START("Find Active Set");
     MSQ_DBGOUT(2) << "\nFinding the active set\n";
 
     // initialize the active set indices to zero
     for (i=0;i<numFunctionValues;i++) active_set->active_ind[i] = 0; 
 
     /* the first function value is our initial active value */
     active_set->num_active = 1;
     active_set->num_equal = 0;
     active_set->active_ind[0] = 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 (i=1;i<numFunctionValues;i++) {
         function_val = function[i];
         active_set->true_active_value = MSQ_MAX(function_val,active_set->true_active_value);
         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 ) {
             if ( temp > activeEpsilon) {
                 active_set->num_active = 1;
                 active_set->num_equal = 0;
                 active_set->active_ind[0] = i;
             } else if ( temp < activeEpsilon) {
                 active_set->num_active += 1;
                 ind = active_set->num_active - 1;
                 active_set->active_ind[ind] = i;
                 if (fabs(function_val - active_value0) < MSQ_MACHINE_EPS) {
                     active_set->num_equal += 1;
                 }
             }
         } else {
             if (temp < activeEpsilon) {
                 active_set->num_active += 1;
                 ind = active_set->num_active - 1;
                 active_set->active_ind[ind] = i;
                 if (fabs(function_val - active_value0) < MSQ_MACHINE_EPS) {
                     active_set->num_equal += 1;
                 }
             }
         }
     }
 
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void find_active_set	(	double *	function,
		ActiveSet *	active_set,
		MsqError &	err
	)

private

void find_plane_normal	(	double	pt1[3],
		double	pt2[3],
		double	pt3[3],
		double *	cross,
		MsqError &	err
	)

inlineprivate

Definition at line 546 of file includeLinks/NonSmoothSteepestDescent.hpp.

References i, and MSQ_NORMALIZE.

Referenced by NonSmoothSteepestDescent::convex_hull_test().

 {
   int i;
   double vecA[3], vecB[3];
 
   for (i=0;i<3;i++) {
     vecA[i] = pt2[i] - pt1[i];
     vecB[i] = pt3[i] - pt1[i];
   }
   cross[0] = vecA[1]*vecB[2] - vecA[2]*vecB[1];
   cross[1] = vecA[2]*vecB[0] - vecA[0]*vecB[2];
   cross[2] = vecA[0]*vecB[1] - vecA[1]*vecB[0];
   MSQ_NORMALIZE(cross, 3);
 }

Here is the caller graph for this function:

void find_plane_normal	(	double	pt1[3],
		double	pt2[3],
		double	pt3[3],
		double *	cross,
		MsqError &	err
	)

private

void find_plane_points	(	int	dir1,
		int	dir2,
		double **	vec,
		int	num_vec,
		double *	pt1,
		double *	pt2,
		double *	pt3,
		int *	opt_status,
		MsqError &	err
	)

private

Definition at line 213 of file QualityImprover/VertexMover/NonSmoothSteepestDescent/NonSmoothSteepestDescent.cpp.

References i, max(), min(), MSQ_BIG_NEG_NMBR, MSQ_BIG_POS_NMBR, MSQ_CCW, MSQ_CHECK_BOTTOM_UP, MSQ_CHECK_TOP_DOWN, MSQ_CHECK_X_COORD_DIRECTION, MSQ_CHECK_Y_COORD_DIRECTION, MSQ_COPY_VECTOR, MSQ_CW, MSQ_EQUIL, MSQ_HULL_TEST_ERROR, MSQ_MACHINE_EPS, MSQ_NO_EQUIL, MSQ_PRINT, and MSQ_TWO_PT_PLANE.

Referenced by NonSmoothSteepestDescent::convex_hull_test().

 {
     int i;
     int ind[50], num_min, num_max;
     int rotate=MSQ_CW;
     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;
 
     *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) < MSQ_MACHINE_EPS) {
         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 */
         MSQ_COPY_VECTOR(pt1,vec[ind[0]],3);
         pt_1 = pt1[dir1]; pt_2 = pt1[dir2];
         if (pt1[dir2] <= 0){rotate=MSQ_CCW; max_inv_slope=MSQ_BIG_NEG_NMBR;}
         if (pt1[dir2] > 0){rotate=MSQ_CW; min_inv_slope=MSQ_BIG_POS_NMBR;}
         switch(rotate) {
         case MSQ_CCW:
           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) > MSQ_MACHINE_EPS)) {
                 ind[1] = i; max_inv_slope=inv_slope; num_rotated = 1;
               } else if (fabs(inv_slope - max_inv_slope) < MSQ_MACHINE_EPS) {
                 ind[2] = i; num_rotated++;
               }
             }
           }
           break;
         case MSQ_CW:
           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) > MSQ_MACHINE_EPS)){
                 ind[1] = i; min_inv_slope=inv_slope; num_rotated = 1;
               } else if (fabs(inv_slope - min_inv_slope) < MSQ_MACHINE_EPS) {
                 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");
           MSQ_COPY_VECTOR(pt2,vec[ind[1]],3); *status = MSQ_TWO_PT_PLANE;
           break;
         default:
           MSQ_PRINT(3)("Found 2 or more points in the rotation\n");
             if (fabs(pt_1) > MSQ_MACHINE_EPS) inv_origin_slope = pt_2/pt_1;
           switch(rotate) {
           case MSQ_CCW:
             if (inv_origin_slope >= max_inv_slope) *status=MSQ_NO_EQUIL;
             else *status=MSQ_CHECK_TOP_DOWN;
             break;
           case MSQ_CW:
             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");
                 MSQ_COPY_VECTOR(pt1,vec[ind[0]],3);
         MSQ_COPY_VECTOR(pt2,vec[ind[1]],3);
         *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) < MSQ_MACHINE_EPS) {
         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 */
         MSQ_COPY_VECTOR(pt1,vec[ind[0]],3);
         pt_1 = pt1[dir1];  pt_2 = pt1[dir2];
         if (pt1[dir2] < 0){rotate=MSQ_CW; min_inv_slope=1E300;}
         if (pt1[dir2] >= 0){rotate=MSQ_CCW; max_inv_slope=-1E300;}
         switch(rotate) {
         case MSQ_CCW:
           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) < MSQ_MACHINE_EPS) {
                 ind[2] = i; num_rotated++;
               }
             }
           }
           break;
         case MSQ_CW:
           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) < MSQ_MACHINE_EPS) {
                 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");
           MSQ_COPY_VECTOR(pt2,vec[ind[1]],3);
           *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 MSQ_CCW:
             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 MSQ_CW:
             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 */
         MSQ_COPY_VECTOR(pt1,vec[ind[0]],3);
         MSQ_COPY_VECTOR(pt2,vec[ind[1]],3);
         *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;
       }
     }
   }
 
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void find_plane_points	(	int	dir1,
		int	dir2,
		double **	vec,
		int	num_vec,
		double *	pt1,
		double *	pt2,
		double *	pt3,
		int *	opt_status,
		MsqError &	err
	)

private

void form_grammian	(	double **	vec,
		MsqError &	err
	)

inlineprivate

Definition at line 639 of file includeLinks/NonSmoothSteepestDescent.hpp.

References i, MsqError::INVALID_STATE, j, NonSmoothSteepestDescent::mActive, NonSmoothSteepestDescent::mDimension, NonSmoothSteepestDescent::mG, MSQ_DOT, MSQ_SETERR, and ActiveSet::num_active.

Referenced by NonSmoothSteepestDescent::check_equilibrium(), and NonSmoothSteepestDescent::search_direction().

 {
    int i, j;
    int num_active = mActive->num_active;
 
    if (num_active > 150) {
       MSQ_SETERR(err)("Exceeded maximum allowed active values",MsqError::INVALID_STATE);
       return;
    }
    /* form the grammian with the dot products of the gradients */
    for (i=0; i<num_active; i++) {
       for (j=i; j<num_active; j++) {
          mG[i][j] = 0.;
          MSQ_DOT(mG[i][j],vec[i],vec[j],mDimension);
          mG[j][i] = mG[i][j];    
       }
    }
 }

Here is the caller graph for this function:

void form_grammian	(	double **	vec,
		MsqError &	err
	)

private

void form_PD_grammian ( MsqError & err )

private

void form_PD_grammian ( MsqError & err )

inlineprivate

Definition at line 830 of file includeLinks/NonSmoothSteepestDescent.hpp.

References ActiveSet::active_ind, i, MsqError::INVALID_STATE, j, k, NonSmoothSteepestDescent::mActive, NonSmoothSteepestDescent::mDimension, NonSmoothSteepestDescent::mG, NonSmoothSteepestDescent::mPDG, MSQ_SETERR, ActiveSet::num_active, NonSmoothSteepestDescent::pdgInd, and NonSmoothSteepestDescent::singular_test().

Referenced by NonSmoothSteepestDescent::search_direction().

 {
     int  i,j,k;
     int  g_ind_1;
     int  singular;
 
     int num_active = mActive->num_active;
         
     /* this assumes the grammian has been formed */
     for (i=0;i<num_active;i++) {
       for (j=0;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<mDimension) && (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++;
     }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void form_reduced_matrix	(	double ***	P,
		MsqError &	err
	)

inlineprivate

Definition at line 973 of file includeLinks/NonSmoothSteepestDescent.hpp.

References i, j, NonSmoothSteepestDescent::mActive, NonSmoothSteepestDescent::mG, and ActiveSet::num_active.

Referenced by NonSmoothSteepestDescent::search_direction().

 {
     int i,j;
     int num_active = mActive->num_active;
 
     (*P)=(double **)malloc(sizeof(double *)*(num_active-1));
     for (i=0; i<num_active-1; i++) 
         (*P)[i]=(double *)malloc(sizeof(double)*(num_active-1));
 
     for (i=0;i<num_active-1;i++) {
         (*P)[i][i] = mG[0][0] - 2*mG[0][i+1] + mG[i+1][i+1];
         for (j=i+1;j<num_active-1;j++) {
             (*P)[i][j] = mG[0][0] - mG[0][j+1] - mG[i+1][0] + mG[i+1][j+1];
             (*P)[j][i] = (*P)[i][j];
         }
     }
 }

Here is the caller graph for this function:

void form_reduced_matrix	(	double ***	P,
		MsqError &	err
	)

private

void get_active_directions	(	double **	gradient,
		double ***	dir,
		MsqError &	err
	)

private

void get_active_directions	(	double **	gradient,
		double ***	dir,
		MsqError &	err
	)

inlineprivate

Definition at line 502 of file includeLinks/NonSmoothSteepestDescent.hpp.

References ActiveSet::active_ind, i, NonSmoothSteepestDescent::mActive, NonSmoothSteepestDescent::mDimension, MSQ_COPY_VECTOR, and ActiveSet::num_active.

Referenced by NonSmoothSteepestDescent::check_equilibrium(), and NonSmoothSteepestDescent::search_direction().

 {
     int i;
     int num_active = mActive->num_active;
 
     (*dir) =(double **)malloc(sizeof(double *)*num_active);
     for (i=0;i<num_active;i++) {
         (*dir)[i] =(double *)malloc(sizeof(double)*mDimension);
         MSQ_COPY_VECTOR((*dir)[i],gradient[mActive->active_ind[i]],mDimension);
     }
 }

Here is the caller graph for this function:

void get_gradient_projections ( MsqError & err )

inlineprivate

Definition at line 1016 of file includeLinks/NonSmoothSteepestDescent.hpp.

References i, NonSmoothSteepestDescent::mDimension, NonSmoothSteepestDescent::mGradient, NonSmoothSteepestDescent::mGS, NonSmoothSteepestDescent::mSearch, MSQ_DOT, MSQ_PRINT, NonSmoothSteepestDescent::mSteepest, and NonSmoothSteepestDescent::numFunctionValues.

Referenced by NonSmoothSteepestDescent::minmax_opt().

 {
     for (int i=0;i<numFunctionValues;i++) 
         MSQ_DOT(mGS[i],mGradient[i],mSearch,mDimension);
 
     MSQ_PRINT(3)("steepest in get_gradient_projections %d\n",mSteepest);
 }

Here is the caller graph for this function:

void get_gradient_projections ( MsqError & err )

private

void get_min_estimate	(	double *	final_est,
		MsqError &	err
	)

inlineprivate

Definition at line 993 of file includeLinks/NonSmoothSteepestDescent.hpp.

References ActiveSet::active_ind, i, NonSmoothSteepestDescent::mActive, NonSmoothSteepestDescent::mAlpha, NonSmoothSteepestDescent::mGS, MSQ_MACHINE_EPS, ActiveSet::num_active, and NonSmoothSteepestDescent::numFunctionValues.

Referenced by NonSmoothSteepestDescent::step_acceptance().

 {
     int    i;
     double est_imp;
 
     *final_est = -1E300;
     for (i=0;i<mActive->num_active;i++) {
         est_imp = -mAlpha*mGS[mActive->active_ind[i]];
         if (est_imp>*final_est) *final_est = est_imp;
     }
     if (*final_est == 0) {
         *final_est = -1E300;
         for (i=0;i<numFunctionValues;i++) {
             est_imp = -mAlpha*mGS[i];
             if ((est_imp>*final_est) && (fabs(est_imp) > MSQ_MACHINE_EPS)) {
                 *final_est = est_imp;
             }
         }
     }
 }

Here is the caller graph for this function:

void get_min_estimate	(	double *	final_est,
		MsqError &	err
	)

private

int improvement_check ( MsqError & err )

private

Definition at line 195 of file QualityImprover/VertexMover/NonSmoothSteepestDescent/NonSmoothSteepestDescent.cpp.

References NonSmoothSteepestDescent::mActive, MSQ_PRINT, NonSmoothSteepestDescent::originalValue, and ActiveSet::true_active_value.

Referenced by NonSmoothSteepestDescent::step_acceptance().

 {
   int improved = 1;
   
   /* 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 );
        improved = 0;
    }
 
   return(improved);
 
 }

Here is the caller graph for this function:

int improvement_check ( MsqError & err )

private

void init_max_step_length ( MsqError & err )

private

void init_max_step_length ( MsqError & err )

inlineprivate

Definition at line 1172 of file includeLinks/NonSmoothSteepestDescent.hpp.

References i, MsqError::INVALID_MESH, j, k, NonSmoothSteepestDescent::maxAlpha, NonSmoothSteepestDescent::mCoords, NonSmoothSteepestDescent::mDimension, MSQ_PRINT, MSQ_SETERR, NonSmoothSteepestDescent::numElements, NonSmoothSteepestDescent::numVertices, and sqrt().

Referenced by NonSmoothSteepestDescent::optimize_vertex_positions().

 {
   int i, j, k;
   double max_diff = 0;
   double diff=0;
 
   MSQ_PRINT(2)("In init_max_step_length\n");
 
   /* check that the input data is correct */
   if (numElements==0) {
     MSQ_SETERR(err)("Num incident vtx = 0\n",MsqError::INVALID_MESH);
     return;
   }
   if ((mDimension!=2) && (mDimension!=3)) {
      MSQ_SETERR(err)("Problem dimension is incorrect", MsqError::INVALID_MESH);
      return;
   }
 
   /* find the maximum distance between two incident vertex locations */
   for (i=1;i<numVertices;i++) {
     for (j=i;j<numVertices+1;j++) {
       diff=0;
       for (k=0;k<mDimension;k++) {
         diff += (mCoords[i][k]-mCoords[j][k])*(mCoords[i][k]-mCoords[j][k]);
       }
       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);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void init_opt ( MsqError & err )

private

void init_opt ( MsqError & err )

inlineprivate

Definition at line 1121 of file includeLinks/NonSmoothSteepestDescent.hpp.

Referenced by NonSmoothSteepestDescent::optimize_vertex_positions().

 {
     int        i, j;
 
     MSQ_PRINT(2)("\nInitializing Optimization \n");
     if (numFunctionValues > 150) {
       MSQ_SETERR(err)("num_values exceeds 150", MsqError::INVALID_STATE);
     }
     /* for the purposes of initialization will be set to zero after */
     equilibriumPt = 0;
     optStatus = 0;
     iterCount = 0;
     optIterCount = 0;
     mSteepest = 0;
     mAlpha = 0;
     maxAlpha = 0;
 
     MSQ_PRINT(3)("  Initialized Constants \n");
     for (i=0;i<3;i++) {
       mSearch[i] = 0;
       pdgInd[i] = -1;
       for (j=0;j<3;j++) mPDG[i][j] = 0;
     }
 
     MSQ_PRINT(3)("  Initialized search and PDG \n");
     for (i=0;i<numFunctionValues;i++) {
        mFunction[i] = 0;
        testFunction[i] = 0;
        originalFunction[i] = 0;
        mGS[i] = 0;
        for (j=0;j<3;j++) {
            mGradient[i][j] = 0;
        }
     }
     MSQ_PRINT(3)("  Initialized function/gradient \n");
     if (numFunctionValues > 150) {
       for (i=0;i<150;i++) {
        for (j=0;j<150;j++) mG[i][j] = -1;
       }
     } else {
       for (i=0;i<numFunctionValues;i++) {
        for (j=0;j<numFunctionValues;j++) mG[i][j] = -1;
       }
     }
     MSQ_PRINT(3)("  Initialized G\n");
  
     for (i=0;i<100;i++) prevActiveValues[i] = 0;
     MSQ_PRINT(3)("  Initialized prevActiveValues\n");
 }

Here is the caller graph for this function:

virtual void initialize	(	PatchData &	pd,
		MsqError &	err
	)

protectedvirtual

Implements VertexMover.

void initialize	(	PatchData &	pd,
		MsqError &	err
	)

protectedvirtual

Implements VertexMover.

Definition at line 73 of file QualityImprover/VertexMover/NonSmoothSteepestDescent/NonSmoothSteepestDescent.cpp.

References NonSmoothSteepestDescent::activeEpsilon, PatchData::ELEMENTS_ON_VERTEX_PATCH, NonSmoothSteepestDescent::minAcceptableImprovement, NonSmoothSteepestDescent::minStepSize, MSQ_DBGOUT, and PatchDataUser::set_patch_type().

 {
   this->set_patch_type(PatchData::ELEMENTS_ON_VERTEX_PATCH, err, 1);
   
   // local parameter initialization
   activeEpsilon = .00003;
   //  activeEpsilon = .000000003;
   minAcceptableImprovement = 1e-6;
   minStepSize = 1e-6;
   MSQ_DBGOUT(1) << "- Executed NonSmoothSteepestDescent::initialize()\n";
 }

Here is the call graph for this function:

virtual void initialize_mesh_iteration	(	PatchData &	pd,
		MsqError &	err
	)

protectedvirtual

Implements VertexMover.

void initialize_mesh_iteration	(	PatchData &	pd,
		MsqError &	err
	)

protectedvirtual

Implements VertexMover.

Definition at line 85 of file QualityImprover/VertexMover/NonSmoothSteepestDescent/NonSmoothSteepestDescent.cpp.

87 {

88 }

void minmax_opt	(	PatchData &	pd,
		MsqError &	err
	)

private

Definition at line 548 of file QualityImprover/VertexMover/NonSmoothSteepestDescent/NonSmoothSteepestDescent.cpp.

Referenced by NonSmoothSteepestDescent::optimize_vertex_positions().

 {
 //      int valid;
       MSQ_FUNCTION_TIMER( "Minmax Opt" );
       MSQ_PRINT(2)("In minmax_opt\n");
 
       MSQ_COPY_VECTOR(mFunction,originalFunction,numFunctionValues);
       originalValue = mActive->true_active_value;
 
       iterCount = 0;
       optIterCount = 0;
 
       MSQ_PRINT(3)("Done copying original function to function\n");
 
       this->find_active_set(mFunction, mActive, err); MSQ_ERRRTN(err);
       prevActiveValues[0] = mActive->true_active_value;
 
      /* check for equilibrium point */
      /* compute the gradient */
      mGradient = this->compute_gradient(&pd, err); MSQ_ERRRTN(err);
      
      if (mActive->num_active >= 2) {
         MSQ_PRINT(3)("Testing for an equilibrium point \n");
         this->check_equilibrium(&equilibriumPt, &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 */
         mGradient = this->compute_gradient(&pd, err); MSQ_ERRRTN(err);
         
         MSQ_PRINT(3)("Computing the search direction \n");
         this->search_direction(pd, 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(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, 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, err);
                 MSQ_ERRRTN(err);
             }
 
             /* check for equilibrium point */
             if (mActive->num_active >= 2) {
                 MSQ_PRINT(3)("Testing for an equilibrium point \n");
                 this->check_equilibrium(&equilibriumPt, &optStatus, err); 
                        MSQ_ERRRTN(err);
 
                 if (MSQ_DBG(2) && equilibriumPt) 
                     MSQ_PRINT(2)("Optimization Exiting: An equilibrium point \n");
             }
 
             /* record the values */
             prevActiveValues[iterCount] = 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,err); MSQ_ERRRTN(err);
             }
         }
       }
 
       MSQ_PRINT(2)("Checking the validity of the mesh\n");
       if (!this->validity_check(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) {
         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;
       }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void minmax_opt	(	PatchData &	pd,
		MsqError &	err
	)

private

virtual void optimize_vertex_positions	(	PatchData &	pd,
		MsqError &	err
	)

protectedvirtual

Implements VertexMover.

void optimize_vertex_positions	(	PatchData &	pd,
		MsqError &	err
	)

protectedvirtual

Implements VertexMover.

Definition at line 89 of file QualityImprover/VertexMover/NonSmoothSteepestDescent/NonSmoothSteepestDescent.cpp.

 {
   MSQ_FUNCTION_TIMER( "NonSmoothSteepestDescent" );
 
   //  cout << "- Executing NonSmoothSteepestDescent::optimize_node_positions()\n";
   /* perform the min max smoothing algorithm */
   MSQ_PRINT(2)("\nInitializing the patch iteration\n");
 
   numVertices = pd.num_vertices();
   MSQ_PRINT(3)("Number of Vertices: %d\n",numVertices);
   numElements = pd.num_elements();
   MSQ_PRINT(3)("Number of Elements: %d\n",numElements);
     //Michael: Note: is this a reliable way to get the dimension?
   mDimension = get_mesh_set()->space_dim();
   MSQ_PRINT(3)("Spatial Dimension: %d\n",mDimension);
 
   numFree=pd.num_free_vertices(err); MSQ_ERRRTN(err);
   MSQ_PRINT(3)("Num Free = %d\n",numFree);
 
   MsqFreeVertexIndexIterator free_iter(&pd, err); MSQ_ERRRTN(err);
   free_iter.reset();
   free_iter.next(); 
   freeVertexIndex = free_iter.value();
   MSQ_PRINT(3)("Free Vertex Index = %d\n",freeVertexIndex);
 
   mCoords = pd.get_vertex_array(err); MSQ_ERRRTN(err);
 
   if (MSQ_DBG(3)) {
     for (int i99=0;i99<numVertices;i99++) {
       MSQ_PRINT(3)("coords: %g %g\n",mCoords[i99][0],mCoords[i99][1]);
     }
   }
 
   mConnectivity = pd.get_element_array(err); MSQ_ERRRTN(err);
   if (MSQ_DBG(3)) {
     msq_std::vector<size_t> indices;
     for (int i99=0;i99<numElements;i99++) {
       mConnectivity[i99].get_vertex_indices(indices);
       MSQ_PRINT(3)("connectivity: %d %d %d\n",indices[0],
               indices[1],indices[2]);
     }
   }
 
   // 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(err)!=1) {
       MSQ_PRINT(1)("ERROR: Invalid mesh\n");
       MSQ_SETERR(err)("Invalid Mesh: Use untangle to create a valid "
                       "triangulation", MsqError::INVALID_MESH);
       return;
   }
 
   /* assumes one function value per element */
   // TODO - need to include vertex metrics 
   numFunctionValues = numElements;
 
   /* initialize the optimization data up to numFunctionValues */
   this->init_opt(err); MSQ_ERRRTN(err);
   this->init_max_step_length(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, err);   MSQ_ERRRTN(err);
 
   this->minmax_opt(pd,err);  MSQ_ERRRTN(err);
 }

Here is the call graph for this function:

void print_active_set	(	ActiveSet *	active_set,
		double *	func,
		MsqError &	err
	)

inlineprivate

Definition at line 1103 of file includeLinks/NonSmoothSteepestDescent.hpp.

References ActiveSet::active_ind, i, MsqError::INVALID_ARG, MSQ_DBGOUT, MSQ_PRINT, MSQ_SETERR, and ActiveSet::num_active.

Referenced by NonSmoothSteepestDescent::minmax_opt(), and NonSmoothSteepestDescent::search_direction().

 {
     if (active_set==0) {
       MSQ_SETERR(err)("Null ActiveSet", MsqError::INVALID_ARG);
       return;
     }
  
     if (active_set->num_active == 0) MSQ_DBGOUT(3)<< "No active values\n";
     /* print the active set */
     for (int i=0;i<active_set->num_active;i++) {
      MSQ_PRINT(3)("Active value %d:   %f \n",
                      i+1,func[active_set->active_ind[i]]); 
     }
 }

Here is the caller graph for this function:

void print_active_set	(	ActiveSet *	active_set,
		double *	func,
		MsqError &	err
	)

private

void search_direction	(	PatchData &	pd,
		MsqError &	err
	)

private

void search_direction	(	PatchData &	pd,
		MsqError &	err
	)

private

Definition at line 406 of file QualityImprover/VertexMover/NonSmoothSteepestDescent/NonSmoothSteepestDescent.cpp.

Referenced by NonSmoothSteepestDescent::minmax_opt().

 {
    int        i;
    int        viable;
    int        singular;
    double     a, b, c, denom;
    double     **dir;
    double     R0, R1;
    double     **P, *x;
    double     search_mag;
 
    int num_active = mActive->num_active;
 
    //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, err);  MSQ_ERRRTN(err);
    
    if (num_active==0) {
        MSQ_SETERR(err)("No active values in search",MsqError::INVALID_STATE);
        return;
     }
 
     switch(num_active) {
     case 1: 
         MSQ_COPY_VECTOR(mSearch,mGradient[mActive->active_ind[0]],mDimension);
         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);  
         if (MSQ_CHKERR(err)) { free(dir); return; }
 
         /* form the grammian */
         this->form_grammian(dir,err); 
         if (MSQ_CHKERR(err)) { free(dir); return; }
         this->form_PD_grammian(err); 
         if (MSQ_CHKERR(err)) { free(dir); return; }
 
         denom = (mG[0][0] + mG[1][1] - 2*mG[0][1]);
         viable = 1;
         if (fabs(denom) > MSQ_MACHINE_EPS) {
           /* gradients are LI, move along their intersection */
            b = (mG[0][0] - mG[0][1])/denom;  
            a = 1 - b;
            if ((b < 0) || (b > 1)) viable=0;  /* 0 < b < 1 */
            if (viable) {
              for (i=0;i<mDimension;i++) {
                mSearch[i] = a*dir[0][i] + b*dir[1][i];
              }
            } else {
              /* the gradients are dependent, move along one face */
              MSQ_COPY_VECTOR(mSearch,dir[0],mDimension);
            }
         } else {
            /* the gradients are dependent, move along one face */
            MSQ_COPY_VECTOR(mSearch,dir[0],mDimension);
         }
         mSteepest = mActive->active_ind[0];
 
         for (i=0;i<num_active;i++) free(dir[i]);
         free(dir);
 
         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);  
         if (MSQ_CHKERR(err)) { free(dir); return; }
         this->form_PD_grammian(err);  
         if (MSQ_CHKERR(err)) { free(dir); return; }
 
         switch(mDimension) {
         case 2:
             this->search_edges_faces(dir,err);  
             if (MSQ_CHKERR(err)) { free(dir); return; }
             break;
         case 3:
           if (num_active == 3) {
               this->singular_test(num_active,mG,&singular,err);   
               if (MSQ_CHKERR(err)) { free(dir); return; }
               if (!singular) {
                 /* form the entries of P=Z^T G Z where Z = [-1...-1; I ] */
                 this->form_reduced_matrix(&P,err);   
                 if (MSQ_CHKERR(err)) { free(dir); return; }
                 /* form  the RHS and solve the system for the coeffs */
                 R0 = mG[0][0] - mG[1][0];  R1 = mG[0][0] - mG[2][0];
                 this->solve2x2(P[0][0],P[0][1],P[1][0],P[1][1],R0,R1,&x,err);
                 if (MSQ_CHKERR(err)) { free(dir); return; }
                 if (x!=NULL) {
                         a = 1 - x[0] - x[1];  b = x[0];  c = x[1];
                         for (i=0;i<mDimension;i++) {
                           mSearch[i] = a*dir[0][i] + b*dir[1][i] + 
                                               c*dir[2][i];
                         }
                         mSteepest = mActive->active_ind[0];
                         for (i=0;i<num_active-1;i++)  free(P[i]);  
                         free(P);  free(x);
                 } else { 
                         this->search_edges_faces(dir, err);
                         if (MSQ_CHKERR(err)) { free(dir); return; }
                 }
               } else {
                  this->search_edges_faces(dir, err); 
                  if (MSQ_CHKERR(err)) { free(dir); return; }
               }
             } else {
               this->search_edges_faces(dir, err);
               if (MSQ_CHKERR(err)) { free(dir); return; }
             }
             break;
         }
         for (i=0;i<num_active;i++) free(dir[i]);
         free(dir);
     }
 
     /* if the search direction is essentially zero, equilibrium pt */
     MSQ_DOT(search_mag,mSearch,mSearch,mDimension);
     MSQ_PRINT(3)("  Search Magnitude %g \n",search_mag);
 
     if (fabs(search_mag)<1E-13) optStatus = MSQ_ZERO_SEARCH;
     else MSQ_NORMALIZE(mSearch,mDimension);
     MSQ_PRINT(3)("  Search Direction %g %g  Steepest %d\n",mSearch[0],mSearch[1],mSteepest);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void search_edges_faces	(	double **	dir,
		MsqError &	err
	)

inlineprivate

Definition at line 872 of file includeLinks/NonSmoothSteepestDescent.hpp.

References ActiveSet::active_ind, denom, i, MsqError::INVALID_MESH, j, k, NonSmoothSteepestDescent::mActive, NonSmoothSteepestDescent::mDimension, NonSmoothSteepestDescent::mG, NonSmoothSteepestDescent::mSearch, MSQ_COPY_VECTOR, MSQ_DOT, MSQ_EQUILIBRIUM, MSQ_MACHINE_EPS, MSQ_SETERR, NonSmoothSteepestDescent::mSteepest, ActiveSet::num_active, NonSmoothSteepestDescent::optStatus, and projection().

Referenced by NonSmoothSteepestDescent::search_direction().

 {
 //    int ierr;
     int i,j,k;
     int viable;
     double a,b,denom;
     double g_bar[3];
     double temp_search[3];
     double projection, min_projection;
 
     int num_active = mActive->num_active;
 
     if ( (mDimension != 2) && (mDimension != 3)) {
        MSQ_SETERR(err)("Dimension must be 2 or 3", MsqError::INVALID_MESH);
     }
 
     /* initialize the search direction to 0,0 */
     for (i=0;i<mDimension;i++) temp_search[i] = 0;
 
     /* Check for viable faces */
     min_projection = 1E300;
     for (i=0; i<num_active; i++) {
         /* FACE I */
         viable = 1;
 
         /* test the viability */
         for (j=0;j<num_active;j++) {       /* lagrange multipliers>0 */
              if (mG[j][i] < 0) viable = 0;
         }
        
         /* find the minimum of viable directions */
         if ((viable) && (mG[i][i] < min_projection)) {
             min_projection = mG[i][i];
             MSQ_COPY_VECTOR(temp_search,dir[i],mDimension);
             mSteepest = mActive->active_ind[i];
         }
     
        /* INTERSECTION IJ */
        for (j=i+1; j<num_active; j++) {
           viable = 1;
 
           /* 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) > MSQ_MACHINE_EPS) {
              b = (mG[i][j] - mG[i][i])/denom;
              a = 1 - b;
              if ((b < 0) || (b > 1)) viable=0;  /* 0 < b < 1 */
              for (k=0;k<num_active;k++) {       /* lagrange multipliers>0 */
                  if ((a*mG[k][i] + b*mG[k][j]) <= 0) viable=0;
              }
           } else {
              viable = 0;                        /* Linearly dependent */
           }
 
           /* find the minimum of viable directions */
           if (viable) {
              for (k=0;k<mDimension;k++) {
                 g_bar[k] = a * dir[i][k] + b * dir[j][k];
              }
              MSQ_DOT(projection,g_bar,g_bar,mDimension);
              if (projection < min_projection) {
                 min_projection = projection;
                 MSQ_COPY_VECTOR(temp_search,g_bar,mDimension);
                 mSteepest = mActive->active_ind[i];
              }
           }
        }
     }
     if (optStatus != MSQ_EQUILIBRIUM) {
         MSQ_COPY_VECTOR(mSearch,temp_search,mDimension);
     }
 }         

Here is the call graph for this function:

Here is the caller graph for this function:

void search_edges_faces	(	double **	dir,
		MsqError &	err
	)

private

void singular_test	(	int	n,
		double **	A,
		int *	singular,
		MsqError &	err
	)

private

void singular_test	(	int	n,
		double **	A,
		int *	singular,
		MsqError &	err
	)

inlineprivate

Definition at line 801 of file includeLinks/NonSmoothSteepestDescent.hpp.

References NonSmoothSteepestDescent::condition3x3(), MsqError::INVALID_ARG, MSQ_FALSE, MSQ_MACHINE_EPS, MSQ_SETERR, and MSQ_TRUE.

Referenced by NonSmoothSteepestDescent::form_PD_grammian(), and NonSmoothSteepestDescent::search_direction().

 {
 //    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)=MSQ_TRUE;
     switch(n) {
     case 1:
         if (A[0][0] > 0) (*singular) = MSQ_FALSE;
         break;
     case 2:
         if (fabs(A[0][0]*A[1][1] - A[0][1]*A[1][0]) > MSQ_MACHINE_EPS)           
             (*singular) = MSQ_FALSE;
         break;
     case 3:
        /* calculate the condition number */
         this->condition3x3(A, &cond, err); 
         if (cond < 1E14) (*singular)=MSQ_FALSE;
         break;
     }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void solve2x2	(	double	a11,
		double	a12,
		double	a21,
		double	a22,
		double	b1,
		double	b2,
		double **	x,
		MsqError &	err
	)

inlineprivate

Definition at line 948 of file includeLinks/NonSmoothSteepestDescent.hpp.

References MSQ_MACHINE_EPS.

Referenced by NonSmoothSteepestDescent::search_direction().

 {
     double factor;
 
     /* if the system is not singular, solve it */
     if (fabs(a11*a22 - a21*a12) > MSQ_MACHINE_EPS) {
         (*x)=(double *)malloc(sizeof(double)*2);
         if (fabs(a11) > MSQ_MACHINE_EPS) {
             factor = (a21/a11);
             (*x)[1] = (b2 - factor*b1)/(a22 - factor*a12);
             (*x)[0] = (b1 - a12*(*x)[1])/a11;
         } else if (fabs(a21) > MSQ_MACHINE_EPS) {
             factor = (a11/a21);
             (*x)[1] = (b1 - factor*b2)/(a12 - factor*a22);
             (*x)[0] = (b2 - a22*(*x)[1])/a21;
         }
     } else {
         (*x) = NULL;
     }
 }

Here is the caller graph for this function:

void solve2x2	(	double	a11,
		double	a12,
		double	a21,
		double	a22,
		double	b1,
		double	b2,
		double **	x,
		MsqError &	err
	)

private

void step_acceptance	(	PatchData &	pd,
		MsqError &	err
	)

private

void step_acceptance	(	PatchData &	pd,
		MsqError &	err
	)

private

Definition at line 669 of file QualityImprover/VertexMover/NonSmoothSteepestDescent/NonSmoothSteepestDescent.cpp.

Referenced by NonSmoothSteepestDescent::minmax_opt().

 {
 //  int        ierr;
   int        i;
   int        num_values, num_steps;
   int        valid = 1, step_status;
   int        accept_alpha;
   double     estimated_improvement;
   double     current_improvement = 1E300;
   double     previous_improvement = 1E300;
   double     current_percent_diff = 1E300;
   double     original_point[3];
 
 //  MSQ_FUNCTION_TIMER( "Step Acceptance" );
   num_values = numFunctionValues;
 
   step_status = MSQ_STEP_NOT_DONE;
   optStatus = 0;
   num_steps = 0;
 
   if (mAlpha < minStepSize) {
       optStatus = MSQ_IMP_TOO_SMALL;
       step_status = MSQ_STEP_DONE;
       MSQ_PRINT(3)("Alpha starts too small, no improvement\n");
   }
 
   /* save the original function and active set */
   MSQ_COPY_VECTOR(original_point,mCoords[freeVertexIndex],mDimension);
   MSQ_COPY_VECTOR(originalFunction, mFunction, num_values);
   this->copy_active(mActive, originalActive, err); MSQ_ERRRTN(err);
 
   while (step_status == MSQ_STEP_NOT_DONE) {
 
     num_steps++;  if (num_steps >= 100) step_status = MSQ_STEP_DONE;
 
     accept_alpha = MSQ_FALSE;
 
     while (!accept_alpha && mAlpha>minStepSize) {
 
       /* make the step */
       for (i=0;i<mDimension;i++) {
          mCoords[freeVertexIndex][i] -= mAlpha*mSearch[i];
       }
         //pd.set_coords_array_element(mCoords[freeVertexIndex],0,err);
 
       MSQ_PRINT(2)("search direction %f %f \n",mSearch[0],mSearch[1]); 
       MSQ_PRINT(2)("new vertex position %f %f \n",mCoords[freeVertexIndex][0],mCoords[freeVertexIndex][1]); 
 
       /* assume alpha is acceptable */
       accept_alpha=MSQ_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(err); MSQ_ERRRTN(err);
       if (valid) valid=improvement_check(err); MSQ_ERRRTN(err);
       if (!valid) {
           accept_alpha=MSQ_FALSE;
           for (i=0;i<mDimension;i++) {
              mCoords[freeVertexIndex][i] += mAlpha*mSearch[i];
           }
             //pd.set_coords_array_element(mCoords[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_status = MSQ_STEP_DONE;
                 MSQ_PRINT(2)("Step too small\n");
                 /* get back the original point, mFunction, and active set */
                 MSQ_COPY_VECTOR(mCoords[freeVertexIndex],original_point,mDimension);
                   //pd.set_coords_array_element(mCoords[freeVertexIndex],0,err);
                 MSQ_COPY_VECTOR(mFunction,originalFunction,num_values);
                 this->copy_active(originalActive, mActive, err); 
           }
        }
     } 
          
     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, err); MSQ_ERRRTN(err);
         
       /* 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) */
         for (i=0;i<mDimension;i++) {
           mCoords[freeVertexIndex][i] -= mAlpha*mSearch[i];
         }
             //pd.set_coords_array_element(mCoords[freeVertexIndex],0,err);
 
         /* does this make an invalid mesh valid? */
    //TODO Validity check revisison
         valid = 1;
         valid = validity_check(err); MSQ_ERRRTN(err);
         if (valid) valid=improvement_check(err); MSQ_ERRRTN(err);
 
         /* copy test function and active set */
         MSQ_COPY_VECTOR(mFunction,testFunction,numFunctionValues);
         this->copy_active(testActive, mActive, err); MSQ_ERRRTN(err);
  
         optStatus = MSQ_STEP_ACCEPTED;  step_status = MSQ_STEP_DONE;
             
         /* check to see that we're still making good improvements */
         if (fabs(previous_improvement) < minAcceptableImprovement) {
           optStatus = MSQ_IMP_TOO_SMALL; step_status = MSQ_STEP_DONE;
           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_status = MSQ_STEP_DONE;
 
         /* 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_status = MSQ_STEP_DONE;
         }
 
       } 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_status = MSQ_STEP_DONE;
         MSQ_PRINT(2)("Opimization Exiting: Flat no improvement\n");
            
         /* get back the original point, function, and active set */
         MSQ_COPY_VECTOR(mCoords[freeVertexIndex],original_point,mDimension);
             //pd.set_coords_array_element(mCoords[freeVertexIndex],0,err);
         MSQ_COPY_VECTOR(mFunction,originalFunction,numFunctionValues);
         this->copy_active(originalActive, mActive, err); MSQ_CHKERR(err);
 
       }
       else
       {
         /* halve alpha and try again */
         /* add out the old step */
         for (i=0;i<mDimension;i++) mCoords[freeVertexIndex][i] += mAlpha*mSearch[i];
             //pd.set_coords_array_element(mCoords[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");
           MSQ_COPY_VECTOR(mCoords[freeVertexIndex],original_point,mDimension);
               //pd.set_coords_array_element(mCoords[freeVertexIndex],0,err);
           MSQ_COPY_VECTOR(mFunction,originalFunction,numFunctionValues);
           this->copy_active(originalActive, mActive, err); MSQ_ERRRTN(err);
           optStatus = MSQ_STEP_TOO_SMALL;  step_status = MSQ_STEP_DONE;
         }
           else
           {
           MSQ_COPY_VECTOR(testFunction, mFunction, numFunctionValues);
           this->copy_active(mActive, testActive, err); MSQ_ERRRTN(err);
           previous_improvement = current_improvement;
         }
       }
     }
   }
   if (current_improvement>0 && optStatus==MSQ_STEP_ACCEPTED) {
     MSQ_PRINT(2)("Accepted a negative step %f \n",current_improvement);
   }
 
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void terminate_mesh_iteration	(	PatchData &	pd,
		MsqError &	err
	)

protectedvirtual

Implements VertexMover.

Definition at line 164 of file QualityImprover/VertexMover/NonSmoothSteepestDescent/NonSmoothSteepestDescent.cpp.

166 {

167 }

virtual void terminate_mesh_iteration	(	PatchData &	pd,
		MsqError &	err
	)

protectedvirtual

Implements VertexMover.

int validity_check ( MsqError & err )

private

int validity_check ( MsqError & err )

inlineprivate

Definition at line 394 of file includeLinks/NonSmoothSteepestDescent.hpp.

References Vector3D::get_coordinates(), MsqMeshEntity::get_vertex_index(), i, NonSmoothSteepestDescent::mConnectivity, NonSmoothSteepestDescent::mCoords, NonSmoothSteepestDescent::mDimension, MSQ_MACHINE_EPS, NonSmoothSteepestDescent::numElements, and sqrt().

Referenced by NonSmoothSteepestDescent::minmax_opt(), NonSmoothSteepestDescent::optimize_vertex_positions(), and NonSmoothSteepestDescent::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
   int valid = 1;
   double dEps = 1.e-13;
 
   double x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4;
 
   if (mDimension == 2)
   {
     for (int i=0;i<numElements;i++)
     {
       double dummy = 0;
       mCoords[mConnectivity[i].get_vertex_index(0)].get_coordinates(x1, y1, dummy);
       mCoords[mConnectivity[i].get_vertex_index(1)].get_coordinates(x2, y2, dummy);
       mCoords[mConnectivity[i].get_vertex_index(2)].get_coordinates(x3, y3, dummy);
       
       double a = x2*y3 - x3*y2;
       double b = y2 - y3;
       double c = x3 - x2;
       
       if (.5*(a+b*x1+c*y1) < .01*MSQ_MACHINE_EPS)
         valid=0;
     }
   }
 
   if (mDimension == 3)
   {
     for (int i=0;i<numElements;i++)
     {
       mCoords[mConnectivity[i].get_vertex_index(0)].get_coordinates(x1, y1, z1);
       mCoords[mConnectivity[i].get_vertex_index(1)].get_coordinates(x2, y2, z2);
       mCoords[mConnectivity[i].get_vertex_index(2)].get_coordinates(x3, y3, z3);
       mCoords[mConnectivity[i].get_vertex_index(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) < MSQ_MACHINE_EPS)
       {
         return(valid = 0);
       }
       else
       {
         dDet /= (dScale*dScale*dScale);
       }
       
       if (dDet > dEps)
       {
         valid = 1;
       }
       else if (dDet < -dEps)
       {
         valid = -1;
       }
       else
       {
         valid = 0;
       }
     }  // end for i=1,numElements
   }  // end mDimension==3
   
   //  MSQ_DEBUG_ACTION(2,{fprintf(stdout,"Mesh Validity is: %d \n",valid);});
   
 //  FUNCTION_TIMER_END();
   return(valid);
 }