C++ PatchData::set_vertex_coordinates方法代码示例

double 
NonGradient::evaluate( double *point,  PatchData &pd, MsqError &err )
{
  if( pd.num_free_vertices() > 1 )
  {
    MSQ_SETERR(err)("Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED);
  }
  const size_t vertexIndex = 0; 
  Vector3D originalVec = pd.vertex_by_index(vertexIndex);
  Vector3D pointVec;
  for( int index = 0; index<3; index++)
  {
    pointVec[index] = point[index];
  }
  pd.set_vertex_coordinates( pointVec, vertexIndex, err ); 
  pd.snap_vertex_to_domain( vertexIndex, err );  

  OFEvaluator& obj_func = get_objective_function_evaluator();
  TerminationCriterion* term_crit=get_inner_termination_criterion();

  double value;
  bool feasible = obj_func.evaluate( pd, value, err ); // MSQ_ERRZERO(err);
  term_crit->accumulate_inner( pd, value, 0, err );  //MSQ_CHKERR(err);
  if (err.error_code() == err.BARRIER_VIOLATED)
    err.clear();   // barrier violation not an error here
  MSQ_ERRZERO(err);
  pd.set_vertex_coordinates( originalVec, vertexIndex, err ); 
  if( !feasible && mUseExactPenaltyFunction )  
  { // "value" undefined btw
    double ensureFiniteRtol= .25;
    value = DBL_MAX * ensureFiniteRtol;
    //std::cout << " Infeasible patch: " << value << std::endl; printPatch( pd, err );
  }
  return value;
}

void SmartLaplacianSmoother::optimize_vertex_positions( PatchData &pd, 
                                                        MsqError &err )
{
  assert(pd.num_free_vertices() == 1);
  const size_t center_vtx_index = 0;
  const size_t init_inverted = num_inverted( pd, err ); MSQ_ERRRTN(err);
  
  adjVtxList.clear();
  pd.get_adjacent_vertex_indices( center_vtx_index, adjVtxList, err );
  MSQ_ERRRTN(err);
  
  if (adjVtxList.empty())
    return;
  
  const MsqVertex* verts = pd.get_vertex_array(err);
  const size_t n = adjVtxList.size();
  
  const Vector3D orig_pos = verts[center_vtx_index];
  Vector3D new_pos = verts[ adjVtxList[0] ];
  for (size_t i = 1; i < n; ++i)
    new_pos += verts[ adjVtxList[i] ];
  new_pos *= 1.0/n;
  pd.set_vertex_coordinates( new_pos, center_vtx_index, err );
  pd.snap_vertex_to_domain( center_vtx_index, err );  MSQ_ERRRTN(err);
  const size_t new_inverted = num_inverted( pd, err ); MSQ_ERRRTN(err);
  if (new_inverted > init_inverted)
    pd.set_vertex_coordinates( orig_pos, center_vtx_index, err );
}

template <class QMType> inline
void TMPQualityMetricTest<QMType>::test_evaluate_surface()
{
  MsqPrintError err(cout);
  PatchData pd;
  bool rval;
  double value;
  
  QMType m( &surf_target, &faux_pi );
  
    // test with aligned elements
  faux_pi.count = 0;
  tester.get_ideal_element( QUADRILATERAL, true, pd );
  rval = m.evaluate( pd, 0, value, err );
  CPPUNIT_ASSERT(!MSQ_CHKERR(err));
  CPPUNIT_ASSERT(rval);
  CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
  CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );
  
    // test that columns are orthogonal for ideal quad element
  CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, col_dot_prod(faux_pi.last_A_2D), 1e-6 );
  
    // test with an element rotated about X-axis
  faux_pi.count = 0;
  tester.get_ideal_element( QUADRILATERAL, true, pd );
  // rotate by 90 degrees about X axis
  for (size_t i = 0; i < pd.num_nodes(); ++i) {
    Vector3D orig = pd.vertex_by_index(i);
    Vector3D newp( orig[0], -orig[2], orig[1] );
    pd.set_vertex_coordinates( newp, i, err );
  }
  rval = m.evaluate( pd, 0, value, err );
  CPPUNIT_ASSERT(!MSQ_CHKERR(err));
  CPPUNIT_ASSERT(rval);
  CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
  CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );
  
    // test with an element rotated about Y-axis
  faux_pi.count = 0;
  tester.get_ideal_element( TRIANGLE, true, pd );
  // rotate by -90 degrees about Y axis
  for (size_t i = 0; i < pd.num_nodes(); ++i) {
    Vector3D orig = pd.vertex_by_index(i);
    Vector3D newp( orig[2], orig[1], -orig[0] );
    pd.set_vertex_coordinates( newp, i, err );
  }
  rval = m.evaluate( pd, 0, value, err );
  CPPUNIT_ASSERT(!MSQ_CHKERR(err));
  CPPUNIT_ASSERT(rval);
  CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
  CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );
}  

//VERTEX BOUND
void TerminationCriterionTest::test_vertex_bound()
{
  MsqPrintError err(std::cout);
  PatchData pd;
  create_twelve_hex_patch( pd, err );
  ASSERT_NO_ERROR(err);
  
    // get bounding dimension for patch
  double maxcoord = 0.0;
  for (size_t i = 0; i < pd.num_nodes(); ++i) 
    for (int d = 0; d < 3; ++d)
      if (fabs(pd.vertex_by_index(i)[d]) > maxcoord)
        maxcoord = fabs(pd.vertex_by_index(i)[d]);
    // add a little bit for rounding error
  maxcoord += 1e-5;
  
  TerminationCriterion tc;
  tc.add_bounded_vertex_movement( maxcoord );
  tc.reset_inner( pd, ofEval, err );
  ASSERT_NO_ERROR(err);
  tc.reset_patch( pd, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT(!tc.terminate());
  
  int idx = pd.num_free_vertices() - 1;
  Vector3D pos = pd.vertex_by_index(idx);
  pos[0] = 2*maxcoord;
  pd.set_vertex_coordinates( pos, idx, err );
  ASSERT_NO_ERROR(err);
  tc.accumulate_inner( pd, 0.0, 0, err );
  ASSERT_NO_ERROR(err);
  tc.accumulate_patch( pd, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT(tc.terminate());
}

/*!Returns an appropiate value (eps) to use as a delta step for
  MsqVertex vertex in dimension k (i.e. k=0 -> x, k=1 -> y, k=2 -> z).
  The objective function value at the perturbed vertex position is given
  in local_val.
*/
double ObjectiveFunction::get_eps( PatchData &pd, 
                                   EvalType type,
                                   double &local_val,
                                   int dim,
                                   size_t vertex_index, 
                                   MsqError& err)
{
  double eps = 1.e-07;
  const double rho = 0.5;
  const int imax = 20;
  bool feasible = false;
  double tmp_var = 0.0;
  Vector3D delta(0,0,0);
  for (int i = 0; i < imax; ++i)
  {
    //perturb kth coord val and check feas if needed
    tmp_var = pd.vertex_by_index(vertex_index)[dim];
    delta[dim] = eps;
    pd.move_vertex( delta, vertex_index, err );
    feasible = evaluate( type, pd, local_val, OF_FREE_EVALS_ONLY, err ); MSQ_ERRZERO(err);
    //revert kth coord val
    delta = pd.vertex_by_index(vertex_index);
    delta[dim] = tmp_var;
    pd.set_vertex_coordinates( delta, vertex_index, err );
    //if step was too big, shorten it and go again
    if (feasible)
      return eps;
    else
      eps *= rho;
  }//end while looking for feasible eps
  return 0.0;
}//end function get_eps

//UNTANGLED
void TerminationCriterionTest::test_untangled_mesh()
{
  MsqPrintError err(std::cout);
  PatchData pd;
  create_twelve_hex_patch( pd, err );
  ASSERT_NO_ERROR(err);
  
    // get two opposite vertices in first hexahedral element
  int vtx1 = pd.element_by_index(0).get_vertex_index_array()[0];
  int vtx2 = pd.element_by_index(0).get_vertex_index_array()[7];
  Vector3D saved_coords = pd.vertex_by_index(vtx2);
  Vector3D opposite_coords = pd.vertex_by_index(vtx1);

    // invert the element
  pd.move_vertex( 2*(opposite_coords-saved_coords), vtx2, err );
  ASSERT_NO_ERROR(err);
  int inverted, samples;
  pd.element_by_index(0).check_element_orientation( pd, inverted, samples, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT(inverted > 0);
  
    // check initial termination criterion
  TerminationCriterion tc;
  tc.add_untangled_mesh();
  tc.reset_inner( pd, ofEval, err );
  ASSERT_NO_ERROR(err);
  tc.reset_patch( pd, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT(!tc.terminate());
  
    // fix the element
  pd.set_vertex_coordinates( saved_coords, vtx2, err );
  ASSERT_NO_ERROR(err);
  pd.element_by_index(0).check_element_orientation( pd, inverted, samples, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT_EQUAL(0,inverted);
  
    // check that TC recognized untangled mesh
  tc.accumulate_inner( pd, 0.0, 0, err );
  ASSERT_NO_ERROR(err);
  tc.accumulate_patch( pd, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT(tc.terminate());
}

void MsqMeshEntityTest::test_check_element_orientation( EntityTopology type, 
                                                        int nodes )
{
    // get an ideal element
  MsqError err;
  PatchData pd;
  create_ideal_element_patch( pd, type, nodes, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT_EQUAL( (size_t)1, pd.num_elements() );
  CPPUNIT_ASSERT_EQUAL( (size_t)nodes, pd.num_nodes() );
  MsqMeshEntity& elem = pd.element_by_index(0);
  CPPUNIT_ASSERT_EQUAL( (size_t)nodes, elem.node_count() );
  CPPUNIT_ASSERT_EQUAL( type, elem.get_element_type() );
  const size_t* conn = elem.get_vertex_index_array();
  
    // test that ideal element is not reported as inverted
  int inverted, tested;
  elem.check_element_orientation( pd, inverted, tested, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT_EQUAL( 0, inverted );
  CPPUNIT_ASSERT( tested > 0 );

  bool mids[4] = {false};
  TopologyInfo::higher_order( type, nodes, mids[1], mids[2], mids[3], err );
  MSQ_ERRRTN(err);
  
    // invert element at each vertex and test
  Vector3D centroid;
  elem.get_centroid( centroid, pd, err );
  ASSERT_NO_ERROR(err);
  for (int i = 0; i < nodes; ++i) {
    unsigned dim, num;
    TopologyInfo::side_from_higher_order( type, nodes, i, dim, num, err );
    ASSERT_NO_ERROR(err);
    const Vector3D old_pos = pd.vertex_by_index( conn[i] );
    Vector3D new_pos = old_pos;
    if (dim == TopologyInfo::dimension(type)) { 
      // move mid-element node 3/4 of the way to corner 0
      new_pos += 3*pd.vertex_by_index( conn[0] );
      new_pos *= 0.25;
    }
    else if (dim == 0) { // if a corner vertex
      if (type == TRIANGLE || type == TETRAHEDRON) {
        // move tri/tet vertex past opposite side of element
        new_pos += 2*(centroid - old_pos);
      }
      else if (mids[1]) {
        // if have mid-edge nodes move 3/4 of the way to center vertex
        new_pos += 3*centroid;
        new_pos *= 0.25;
      }
      else {
        // move vertex past centroid
        new_pos += 1.5*(centroid - old_pos);
      }
    }
    else {
      // otherwise move vertex past centroid
      new_pos += 2.5*(centroid - old_pos);
    }
    
    pd.set_vertex_coordinates( new_pos, conn[i], err );
    ASSERT_NO_ERROR(err);
    
      // test that element is inverted
    inverted = tested = 0;
    elem.check_element_orientation( pd, inverted, tested, err );
    ASSERT_NO_ERROR(err);
    std::ostringstream str;
    str << TopologyInfo::short_name(type) << nodes
        << " Vertex " << i 
        << " (Dimension " << dim
        << " Index " << num << ")";
    CppUnit::Message m( "MsqMeshEntity failed to detect inverted element" );
    m.addDetail( str.str() );
    ASSERT_MESSAGE( m, inverted > 0 );
    
      // move vertex back to ideal position
    pd.set_vertex_coordinates( old_pos, conn[i], err );
    ASSERT_NO_ERROR(err);
  }
}

bool QualityMetric::evaluate_with_Hessian( PatchData& pd,
                              size_t handle,
                              double& value,
                              std::vector<size_t>& indices,
                              std::vector<Vector3D>& gradient,
                              std::vector<Matrix3D>& Hessian,
                              MsqError& err )
{
  indices.clear();
  gradient.clear();
  keepFiniteDiffEps = true;
  bool valid = evaluate_with_gradient( pd, handle, value, indices, gradient, err );
  keepFiniteDiffEps = false;
  if (MSQ_CHKERR(err) || !valid) {
    haveFiniteDiffEps = false;
    return false;
  }
  if (indices.empty()){
    haveFiniteDiffEps = false;
    return true;
  }

    // get initial pertubation amount
  double delta_C;
  if (haveFiniteDiffEps) {
    delta_C = finiteDiffEps;
  }
  else {
    delta_C = get_delta_C( pd, indices, err ); MSQ_ERRZERO(err);
  }
  assert(delta_C < 1e30);
  const double delta_inv_C = 1.0/delta_C;
  const int reduction_limit = 15;

  std::vector<Vector3D> temp_gradient( indices.size() );
  const int num_hess = indices.size() * (indices.size() + 1) / 2;
  Hessian.resize( num_hess );
  
  for (unsigned v = 0; v < indices.size(); ++v) {
    const Vector3D pos = pd.vertex_by_index(indices[v]);
    for (int j = 0; j < 3; ++j ) { // x, y, and z
      double delta = delta_C;
      double delta_inv = delta_inv_C;
      double metric_value;
      Vector3D delta_v(0,0,0);
      
        // find finite difference for gradient
      int counter = 0;
      for (;;) {
        delta_v[j] = delta;
        pd.set_vertex_coordinates( pos+delta_v, indices[v], err ); MSQ_ERRZERO(err);
        valid = evaluate_with_gradient( pd, handle, metric_value, indices, temp_gradient, err );
        if (!MSQ_CHKERR(err) && valid) 
          break;
        
        if (++counter >= reduction_limit) {
          MSQ_SETERR(err)("Algorithm did not successfully compute element's "
                           "Hessian.\n",MsqError::INTERNAL_ERROR);
          haveFiniteDiffEps = false;
          return false;
        }
        
        delta *= 0.1;
        delta_inv *= 10.0;
      }
      pd.set_vertex_coordinates( pos, indices[v], err ); MSQ_ERRZERO(err);
      
        //compute the numerical Hessian
      for (unsigned w = 0; w <= v; ++w) {
          //finite difference to get some entries of the Hessian
        Vector3D fd( temp_gradient[w] );
        fd -= gradient[w];
        fd *= delta_inv;
          // For the block at position w,v in a matrix, we need the corresponding index
          // (mat_index) in a 1D array containing only upper triangular blocks.
        unsigned sum_w = w*(w+1)/2; // 1+2+3+...+w
        unsigned mat_index = w*indices.size() + v - sum_w;
        Hessian[mat_index][0][j] = fd[0];
        Hessian[mat_index][1][j] = fd[1];
        Hessian[mat_index][2][j] = fd[2];
      }
    }  // for(j)
  } // for(indices)
  haveFiniteDiffEps = false;
  return true;
}

bool QualityMetric::evaluate_with_gradient( PatchData& pd,
                              size_t handle,
                              double& value,
                              std::vector<size_t>& indices,
                              std::vector<Vector3D>& gradient,
                              MsqError& err )
{
  indices.clear();
  bool valid = evaluate_with_indices( pd, handle, value, indices, err);
  if (MSQ_CHKERR(err) || !valid)
    return false;
  if (indices.empty())
    return true;

    // get initial pertubation amount
  double delta_C = finiteDiffEps;
  if (!haveFiniteDiffEps) {
    delta_C = get_delta_C( pd, indices, err ); MSQ_ERRZERO(err);
    if (keepFiniteDiffEps) {
      finiteDiffEps = delta_C;
      haveFiniteDiffEps = true;
    }
  }
  const double delta_inv_C = 1.0/delta_C;
  const int reduction_limit = 15;

  gradient.resize( indices.size() );
  for (size_t v=0; v<indices.size(); ++v) 
  {
    const Vector3D pos = pd.vertex_by_index(indices[v]);
    
    /* gradient in the x, y, z direction */
    for (int j=0;j<3;++j) 
    {
      double delta = delta_C;
      double delta_inv = delta_inv_C;
      double metric_value;
      Vector3D delta_v( 0, 0, 0 );
      
        //perturb the node and calculate gradient.  The while loop is a
        //safety net to make sure the epsilon perturbation does not take
        //the element out of the feasible region.
      int counter = 0;
      for (;;) {
          // perturb the coordinates of the free vertex in the j direction
          // by delta       
        delta_v[j] = delta;
        pd.set_vertex_coordinates( pos+delta_v, indices[v], err ); MSQ_ERRZERO(err);

          //compute the function at the perturbed point location
        valid = evaluate( pd, handle, metric_value, err);
        if (!MSQ_CHKERR(err) && valid) 
          break;
          
        if (++counter >= reduction_limit) {
          MSQ_SETERR(err)("Perturbing vertex by delta caused an inverted element.",
                          MsqError::INTERNAL_ERROR);
          return false;
        }
        
        delta*=0.1;
        delta_inv*=10.;
      }
        // put the coordinates back where they belong
      pd.set_vertex_coordinates( pos, indices[v], err );
        // compute the numerical gradient
      gradient[v][j] = (metric_value - value) * delta_inv;
    } // for(j)
  } // for(indices)
  return true;
}

void ObjectiveFunctionTests::compare_numerical_hessian( ObjectiveFunction* of,
                                                        bool diagonal_only )
{
  const double delta = 0.0001;

  MsqPrintError err(std::cout);
  PatchData pd;
  create_qm_two_tet_patch( pd, err ); 
  ASSERT_NO_ERROR( err );
  CPPUNIT_ASSERT( pd.num_free_vertices() != 0 );
  
    // get analytical Hessian from objective function
  std::vector<Vector3D> grad;
  std::vector<SymMatrix3D> diag;
  MsqHessian hess;
  hess.initialize( pd, err );
  ASSERT_NO_ERROR( err );
  double value;
  bool valid;
  if (diagonal_only)
    valid = of->evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, pd, value, grad, diag, err );
  else
    valid = of->evaluate_with_Hessian( ObjectiveFunction::CALCULATE, pd, value, grad, hess, err );
  ASSERT_NO_ERROR(err); CPPUNIT_ASSERT(valid);

  
    // do numerical approximation of each block and compare to analytical value
  for (size_t i = 0; i < pd.num_free_vertices(); ++i) {
    const size_t j_end = diagonal_only ? i+1 : pd.num_free_vertices();
    for (size_t j = i; j < j_end; ++j) {
        // do numerical approximation for block corresponding to
        // coorindates for ith and jth vertices.
      Matrix3D block;    
      for (int k = 0; k < 3; ++k) {
        for (int m = 0; m < 3; ++m) {
          double dk, dm, dkm;
          Vector3D ik = pd.vertex_by_index(i);
          Vector3D im = pd.vertex_by_index(j);
          
          Vector3D delta_k(0.0); delta_k[k] = delta;
          pd.move_vertex( delta_k, i, err ); ASSERT_NO_ERROR(err);
          valid = of->evaluate( ObjectiveFunction::CALCULATE, pd, dk, true, err );
          ASSERT_NO_ERROR(err); CPPUNIT_ASSERT(valid);
          
          Vector3D delta_m(0.0); delta_m[m] = delta;
          pd.move_vertex( delta_m, j, err ); ASSERT_NO_ERROR(err);
          valid = of->evaluate( ObjectiveFunction::CALCULATE, pd, dkm, true, err );
          ASSERT_NO_ERROR(err); CPPUNIT_ASSERT(valid);
          
            // be careful here that we do the right thing if i==j
          pd.set_vertex_coordinates( ik, i, err ); ASSERT_NO_ERROR(err);
          pd.set_vertex_coordinates( im, j, err ); ASSERT_NO_ERROR(err);
          pd.move_vertex( delta_m, j, err ); ASSERT_NO_ERROR(err);
          valid = of->evaluate( ObjectiveFunction::CALCULATE, pd, dm, true, err );
          ASSERT_NO_ERROR(err); CPPUNIT_ASSERT(valid);
          
          pd.set_vertex_coordinates( ik, i, err ); ASSERT_NO_ERROR(err);
          pd.set_vertex_coordinates( im, j, err ); ASSERT_NO_ERROR(err);
          
          block[k][m] = (dkm - dk - dm + value)/(delta*delta);
        }
      }
        // compare to analytical value
      if (diagonal_only) {
        CPPUNIT_ASSERT(i == j); // see j_end above
        CPPUNIT_ASSERT(i < diag.size());
        CHECK_EQUAL_MATRICES( block, Matrix3D(diag[i]) );
      }
      else {
        Matrix3D* m = hess.get_block( i, j );
        Matrix3D* mt = hess.get_block( j, i );
        if (NULL != m) {
          CHECK_EQUAL_MATRICES( block, *m );
        }
        if (NULL != mt) {
          CHECK_EQUAL_MATRICES( transpose(block), *m );
        }
        if (NULL == mt && NULL == m) {
          CHECK_EQUAL_MATRICES( Matrix3D(0.0), block );
        }
      }
    }
  }
}

void NonSmoothDescent::step_acceptance(PatchData &pd, MsqError &err)
{
//  int        ierr;
  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");
  }

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

  /* save the original function and active set */
  MSQ_COPY_VECTOR(original_point,coords[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 */
      pd.move_vertex( -mAlpha*Vector3D(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=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(pd,err); MSQ_ERRRTN(err);
      if (valid) valid=improvement_check(err); MSQ_ERRRTN(err);
      if (!valid) {
          accept_alpha=MSQ_FALSE;
          pd.move_vertex( mAlpha * Vector3D(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_status = MSQ_STEP_DONE;
                MSQ_PRINT(2)("Step too small\n");
 	        /* get back the original point, mFunction, and active set */
                pd.set_vertex_coordinates( Vector3D(original_point), freeVertexIndex, err );
                  //pd.set_coords_array_element(coords[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");
 
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  
/*
      if( height[0] == 0.)
      {
         MSQ_SETERR(err)("(B) Zero objective function value", MsqError::INTERNAL_ERROR);
         exit(-1);
      }
*/
      if (ytry <= height[ilo])   
      {
        ytry=amotry(simplex,height,&rowSum[0],ihi,-2.0,pd,err);
        if( mNonGradDebug >= 3 ) 
        {      
         std::cout << "Reflect and Expand from highPt " << ytry << std::endl;
        }      
        //MSQ_PRINT(3)("Reflect and Expand from highPt : %e\n",ytry);
      }
      else 
      {
        if (ytry >= height[inhi]) 
        {
          ysave=height[ihi]; // Contract along highPt
          ytry=amotry(simplex,height,&rowSum[0],ihi,0.5,pd,err);
          if (ytry >= ysave)
          { // contract all directions toward lowPt
            for (int col=0;col<numCol;col++)
            {
              if (col != ilo)
              {
                for (int row=0;row<numRow;row++)
                {
                  rowSum[row]=0.5*(simplex[row+col*numRow]+simplex[row+ilo*numRow]);
                  simplex[row+col*numRow]=rowSum[row];
                }
                height[col] = evaluate(&rowSum[0], pd, err); 
                if( mNonGradDebug >= 3 ) 
                {      
                  std::cout << "Contract all directions toward lowPt value( " << col << " ) = " << height[col] << " ilo = " << ilo << std::endl;
                }      
                //MSQ_PRINT(3)("Contract all directions toward lowPt value( %d ) = %e    ilo = %d\n", col, height[col], ilo);
              }
            }
          }
        }   
      } // ytri > h(ilo) 
    } // after evaluation
    ilo=1; // conditional operator or inline if 
    ihi = height[0] > height[1] ? (inhi=1,0) : (inhi=0,1);
    for (int col=0;col<numCol;col++)
    {
      if (height[col] <= height[ilo])
      {
        ilo=col;  // ilo := argmin height
      }
      if (height[col] > height[ihi])
      {
        inhi=ihi;
        ihi=col;
      } 
      else  // height[ihi] >= height[col]
        if (col != ihi && height[col] > height[inhi] ) inhi=col;
    }
//    rtol=2.0*fabs( height[ihi]-height[ilo] )/
//         ( fabs(height[ihi])+fabs(height[ilo])+threshold );
    afterEvaluation = true;
  } //  while not converged 

  // Always set to current best mesh.
  { 
    if( ilo != 0 )
    {
      double yTemp = height[0];
      height[0] = height[ilo]; // height dimension numCol
      height[ilo] = yTemp;
      for (int row=1;row<numRow;row++)
      { 
          yTemp = simplex[row];
          simplex[row] = simplex[row+ilo*numRow];
          simplex[row+ilo*numRow] = yTemp;
      }
    }
  }
  if( pd.num_free_vertices() > 1 )
  {
    MSQ_SETERR(err)("Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED);
  }

  Vector3D newPoint( &simplex[0] ); 
  size_t vertexIndex = 0; // fix c.f. freeVertexIndex
  pd.set_vertex_coordinates( newPoint, vertexIndex, err ); 
  pd.snap_vertex_to_domain( vertexIndex, err );
  if( term_crit->terminate() )
  {
    if( mNonGradDebug >= 1 ) 
    {      
         std::cout << "Optimization Termination OptStatus: Max Iter Exceeded" << std::endl;
    }      
    //MSQ_PRINT(1)("Optimization Termination OptStatus: Max Iter Exceeded\n"); 
  }
}