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

//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

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 );
}  

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 MappingFunction3D::jacobian( const PatchData& pd,
                                  size_t element_number,
                                  NodeSet nodeset,
                                  Sample location,
                                  size_t* vertex_patch_indices_out,
                                  MsqVector<3>* d_coeff_d_xi_out,
                                  size_t& num_vtx_out,
                                  MsqMatrix<3,3>& jacobian_out,
                                  MsqError& err ) const
{
  const MsqMeshEntity& elem = pd.element_by_index( element_number );
  const size_t* conn = elem.get_vertex_index_array();
  
  derivatives( location, nodeset, vertex_patch_indices_out,
               d_coeff_d_xi_out, num_vtx_out, err ); MSQ_ERRRTN(err);
 
  convert_connectivity_indices( elem.node_count(), vertex_patch_indices_out, 
                                num_vtx_out, err );  MSQ_ERRRTN(err);
 
  jacobian_out.zero();
  size_t w = 0;
  for (size_t r = 0; r < num_vtx_out; ++r) {
    size_t i = conn[vertex_patch_indices_out[r]];
    MsqMatrix<3,1> coords( pd.vertex_by_index( i ).to_array() );
    jacobian_out += coords * transpose(d_coeff_d_xi_out[r]);
    if (i < pd.num_free_vertices()) {
      vertex_patch_indices_out[w] = i;
      d_coeff_d_xi_out[w] = d_coeff_d_xi_out[r];
      ++w;
    }
  }
  num_vtx_out = w;
}

void UnOptimizer::optimize_vertex_positions( PatchData &pd, MsqError &err) {
  assert( pd.num_free_vertices() == 1 && pd.vertex_by_index(0).is_free_vertex() );
  std::vector<Vector3D> grad(1);
  double val, junk, coeff;
  bool state;
  
  state = objectiveFunction->evaluate_with_gradient( ObjectiveFunction::CALCULATE,
                                                     pd, val, grad, err );
  MSQ_ERRRTN(err);
  if (!state) {
    MSQ_SETERR(err)(MsqError::INVALID_MESH);
    return;
  }
  grad[0] /= grad[0].length();
  
  PatchDataVerticesMemento* memento = pd.create_vertices_memento( err ); MSQ_ERRRTN(err);
  std::auto_ptr<PatchDataVerticesMemento> deleter( memento );
  pd.get_minmax_edge_length( junk, coeff );
  
  for (int i = 0; i < 100; ++i) {
    pd.set_free_vertices_constrained( memento, &grad[0], 1, coeff, err ); MSQ_ERRRTN(err);
    state = objectiveFunction->evaluate( ObjectiveFunction::CALCULATE, pd, val, true, err );
    MSQ_ERRRTN(err);
    if (state)
      break;
    coeff *= 0.5;
  }
  if (!state) {
    pd.set_to_vertices_memento( memento, err );
  }
}

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

static inline 
double get_delta_C( const PatchData& pd,
                    const std::vector<size_t>& indices,
                    MsqError& err )
{
  if (indices.empty()) {
    MSQ_SETERR(err)(MsqError::INVALID_ARG);
    return 0.0;
  }

  std::vector<size_t>::const_iterator beg, iter, iter2, end;
  
  std::vector<size_t> tmp_vect;
  if (indices.size() == 1u) {
    pd.get_adjacent_vertex_indices( indices.front(), tmp_vect, err );
    MSQ_ERRZERO(err);
    assert(!tmp_vect.empty());
    tmp_vect.push_back( indices.front() );
    beg = tmp_vect.begin();
    end = tmp_vect.end();
  }
  else {
    beg = indices.begin();
    end = indices.end();
  }
  
  double min_dist_sqr = HUGE_VAL, sum_dist_sqr = 0.0;
  for (iter = beg; iter != end; ++iter) {
    for (iter2 = iter+1; iter2 != end; ++iter2) {
      Vector3D diff = pd.vertex_by_index(*iter);
      diff -= pd.vertex_by_index(*iter2);
      double dist_sqr = diff % diff;
      if (dist_sqr < min_dist_sqr)
        min_dist_sqr = dist_sqr;
      sum_dist_sqr += dist_sqr;
    }
  }
  
  return 3e-6*sqrt(min_dist_sqr) + 5e-7*sqrt(sum_dist_sqr/(end-beg));
}

  void test_hard_fixed_flags()
  {   
     MsqPrintError err(cout);
     int indices[10];
     int i=0;
     MsqFreeVertexIndexIterator ind(pd, err);
     ind.reset();
     while (ind.next()) {
        indices[i] = ind.value();
//         cout << i << "th free vertex value: " << ind.value() << endl; 
        ++i;
     } 

     CPPUNIT_ASSERT(i==2); // number of free vertices.
     CPPUNIT_ASSERT(pd.vertex_by_index(indices[0]).is_free_vertex());
     CPPUNIT_ASSERT(pd.vertex_by_index(indices[1]).is_free_vertex());
  }

  void test_soft_fixed_flags()
  {   
     MsqPrintError err(cout);
     pd.set_vertex_culled( 0 );

     int indices[10];
     int i=0;
     MsqFreeVertexIndexIterator ind(pd, err);
     ind.reset();
     while (ind.next()) {
        indices[i] = ind.value();
        ++i;
     } 

     CPPUNIT_ASSERT(i==1); // number of free vertices.
     CPPUNIT_ASSERT(pd.vertex_by_index(indices[0]).is_free_vertex());
  }

void EdgeQM::get_edge_evaluations( PatchData& pd, 
                                   std::vector<size_t>& handles,
                                   bool free_vertices_only, 
                                   bool single_pass_evaluate,
                                   MsqError& err )
{
  std::vector<EdgeData> vtx_edges;
  size_t n_verts = free_vertices_only ? pd.num_free_vertices() : pd.num_nodes();
  size_t n_cutoff = single_pass_evaluate ? pd.num_nodes() : n_verts;
  handles.clear();

  for (size_t i = 0; i < n_verts; ++i) {
    if (!pd.vertex_by_index(i).is_flag_set( MsqVertex::MSQ_PATCH_VTX ))
      continue;

    vtx_edges.clear();

    size_t n_elems;
    const size_t* elems;
    elems = pd.get_vertex_element_adjacencies( i, n_elems, err );
    MSQ_ERRRTN(err);

    for (size_t j = 0; j < n_elems; ++j) {
      MsqMeshEntity& elem = pd.element_by_index(elems[j]);
      unsigned n_edges = TopologyInfo::edges( elem.get_element_type() );
      for (unsigned k = 0; k < n_edges; ++k) {
        const unsigned* edge = TopologyInfo::edge_vertices( elem.get_element_type(), k, err );
        MSQ_ERRRTN(err);

        size_t vtx1 = elem.get_vertex_index_array()[edge[0]];
        size_t vtx2 = elem.get_vertex_index_array()[edge[1]];
        size_t other;
        if (vtx1 == i)
          other = vtx2;
        else if (vtx2 == i)
          other = vtx1;
        else
          continue;
        
          // If !free_vertices_only, we'll visit every edge twice.  
          // The first check below ensures that we only add each edge
          // once.  The second check is never true unless free_vertices_only
          // is true and single_pass_evaluate is false.  In that case, it 
          // serves as an exception to the first rule for those cases in which 
          // we visit an edge only once.  For single_pass_evaluate (e.g.
          // BCD initialization or QualityAssessor) we want to avoid visiting
          // and edge more than once for every patch rather than just within
          // this patch.
        if (other > i || other > n_cutoff) {
          EdgeData d = { other, elems[j], k };
          vtx_edges.push_back(d);
        }
      } // end for each edge in element
    } // end for each element adjacent to vertex
    
    std::sort( vtx_edges.begin(), vtx_edges.end() );
    std::vector<EdgeData>::iterator it, end;
    end = std::unique( vtx_edges.begin(), vtx_edges.end() );
    for (it = vtx_edges.begin(); it != end; ++it) 
      handles.push_back( handle( it->elemEdge, it->adjElem ) );
  } // end for each (free) vertex
}

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 AffineMapMetric::evaluate( PatchData& pd, size_t handle, double& value, MsqError& err )
{
  Sample s = ElemSampleQM::sample( handle );
  size_t e = ElemSampleQM::  elem( handle );
  MsqMeshEntity& elem = pd.element_by_index( e );
  EntityTopology type = elem.get_element_type();
  unsigned edim = TopologyInfo::dimension( type );
  const size_t* conn = elem.get_vertex_index_array();
  
    // This metric only supports sampling at corners, except for simplices.
    // If element is a simpex, then the Jacobian is constant over a linear 
    // element.  In this case, always evaluate at any vertex.
  //unsigned corner = s.number;
  if (s.dimension != 0) {
    if (type == TRIANGLE || type == TETRAHEDRON)
      /*corner = 0*/;
    else {
      MSQ_SETERR(err)("Invalid sample point for AffineMapMetric", MsqError::UNSUPPORTED_ELEMENT );
      return false;
    }
  }
  
  bool rval;
  if (edim == 3) { // 3x3 or 3x2 targets ?
    Vector3D c[3] = { Vector3D(0,0,0), Vector3D(0,0,0), Vector3D(0,0,0) };
    unsigned n;
    const unsigned* adj = TopologyInfo::adjacent_vertices( type, s.number, n );
    c[0] = pd.vertex_by_index( conn[adj[0]] ) - pd.vertex_by_index( conn[s.number] );
    c[1] = pd.vertex_by_index( conn[adj[1]] ) - pd.vertex_by_index( conn[s.number] );
    c[2] = pd.vertex_by_index( conn[adj[2]] ) - pd.vertex_by_index( conn[s.number] );
    MsqMatrix<3,3> A;
    A.set_column( 0, MsqMatrix<3,1>(c[0].to_array()) );
    A.set_column( 1, MsqMatrix<3,1>(c[1].to_array()) );
    A.set_column( 2, MsqMatrix<3,1>(c[2].to_array()) );
    if (type == TETRAHEDRON)
      A = A * TET_XFORM;

    MsqMatrix<3,3> W;
    targetCalc->get_3D_target( pd, e, s, W, err ); MSQ_ERRZERO(err);
    rval = targetMetric->evaluate( A * inverse(W), value, err ); MSQ_ERRZERO(err);
  }
  else {
    Vector3D c[2] = { Vector3D(0,0,0), Vector3D(0,0,0) };
    unsigned n;
    const unsigned* adj = TopologyInfo::adjacent_vertices( type, s.number, n );
    c[0] = pd.vertex_by_index( conn[adj[0]] ) - pd.vertex_by_index( conn[s.number] );
    c[1] = pd.vertex_by_index( conn[adj[1]] ) - pd.vertex_by_index( conn[s.number] );
    MsqMatrix<3,2> App;
    App.set_column( 0, MsqMatrix<3,1>(c[0].to_array()) );
    App.set_column( 1, MsqMatrix<3,1>(c[1].to_array()) );
    
    MsqMatrix<3,2> Wp;
    targetCalc->get_surface_target( pd, e, s, Wp, err ); MSQ_ERRZERO(err);

    MsqMatrix<2,2> A, W;
    MsqMatrix<3,2> RZ;
    surface_to_2d( App, Wp, W, RZ );
    A = transpose(RZ) * App;
    if (type == TRIANGLE)
      A = A * TRI_XFORM;
    
    rval = targetMetric->evaluate( A*inverse(W), value, err ); MSQ_ERRZERO(err);
  }
  
    // apply target weight to value
  if (weightCalc) {
    double ck = weightCalc->get_weight( pd, e, s, err ); MSQ_ERRZERO(err);
    value *= ck;
  }
  return rval;
}

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;
}