C++ PatchData类代码示例

int NonGradient::getPatchDimension(const PatchData &pd, MsqError &err)
{
  const size_t numElt = pd.num_elements();
  unsigned edimMax = 0;  // in case numElt == 0
  for( size_t elementId = 0; elementId < numElt; elementId++)
  {
    const MsqMeshEntity& element = pd.element_by_index( elementId );
    EntityTopology type = element.get_element_type();
    unsigned edim = TopologyInfo::dimension( type );
    if( elementId == 0 ) 
    {
      edimMax = edim;
    } 
    else
    { 
      if(edimMax != edim)
      {
        MSQ_SETERR(err)("A patch has elements of different dimensions", MsqError::INVALID_MESH);
        std::cout << "A patch has elements of different dimensions" << edimMax << "  "  << edim << std::endl;
        if(edimMax < edim)
        {
          edimMax = edim;
        }
      }
    }
  }
  return( edimMax );
}

void NonGradient::printPatch(const PatchData &pd, MsqError &err)
{
  if( mNonGradDebug==0 )
  {
    return;
  }
  const size_t numNode = pd.num_nodes();   //27,  27 what?
  MSQ_PRINT(3)("Number of Vertices: %d\n",(int)pd.num_nodes());
  const size_t numVert = pd.num_free_vertices(); // 1
  MSQ_PRINT(3)("Num Free = %d\n",(int)pd.num_free_vertices());
  const size_t numSlaveVert = pd.num_slave_vertices(); //0
  const size_t numCoin = pd.num_corners(); // 64
  const MsqVertex* coord = pd.get_vertex_array(err);
  MSQ_PRINT(3)("Number of Vertices: %d\n",(int)pd.num_nodes());

  std::cout << "Patch " << numNode << "  " << numVert << "  " << numSlaveVert << "  " << numCoin << std::endl;
  MSQ_PRINT(3)("");
  std::cout << "Coordinate ";
  std::cout << "         " << std::endl;
  for( size_t index = 0; index < numVert; index++ )
  {
    std::cout << coord[index][0] << "  " << coord[index][1] << "  " << coord[index][2] << std::endl;
  }
  //const size_t numElt = pd.num_elements();
  if( mNonGradDebug >= 3 ) 
  {      
         std::cout << "Number of Elements: " << pd.num_elements() << std::endl;
  }      
  MSQ_PRINT(3)("Number of Elements: %d\n",(int)pd.num_elements());
}

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

void MsqMeshEntityTest::test_all_nodes( EntityTopology type, unsigned num_nodes )
{
  const unsigned num_vtx = 27;
  double coords[3*num_vtx] = {0.0};
  size_t conn[num_vtx];
  for (size_t i = 0; i < num_vtx; ++i)
    conn[i] = i;
  bool fixed[num_vtx] = {false};
  CPPUNIT_ASSERT(num_nodes <= num_vtx);
  
  MsqError err;
  PatchData pd;
  size_t n = num_nodes;
  pd.fill( num_nodes, coords, 1, &type, &n, conn, fixed, err );
  ASSERT_NO_ERROR(err);
  
  MsqMeshEntity& elem = pd.element_by_index(0);
  NodeSet all = elem.all_nodes(err);
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT_EQUAL( num_nodes, all.num_nodes() );
  CPPUNIT_ASSERT( all.have_any_corner_node() );
  bool mid_edge, mid_face, mid_reg;
  TopologyInfo::higher_order( type, num_nodes, mid_edge, mid_face, mid_reg, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT_EQUAL( mid_edge, !!all.have_any_mid_edge_node() );
  CPPUNIT_ASSERT_EQUAL( mid_face, !!all.have_any_mid_face_node() );
  CPPUNIT_ASSERT_EQUAL( mid_reg,  !!all.have_any_mid_region_node() );
  
}

void ObjectiveFunctionTests::compare_numerical_gradient( ObjectiveFunction* of )
{
  MsqPrintError err(std::cout);
  PatchData pd;
  create_twelve_hex_patch( pd, err ); 
  ASSERT_NO_ERROR( err );
  
  std::vector<Vector3D> num_grad, ana_grad;
  double num_val, ana_val;
  bool valid;
  
  valid = of->evaluate_with_gradient( ObjectiveFunction::CALCULATE, pd, ana_val, ana_grad, err );
  ASSERT_NO_ERROR( err );
  CPPUNIT_ASSERT(valid);
  CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), ana_grad.size() );
  
  valid = of->ObjectiveFunction::evaluate_with_gradient( ObjectiveFunction::CALCULATE, pd, num_val, num_grad, err );
  ASSERT_NO_ERROR( err );
  CPPUNIT_ASSERT(valid);
  CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), num_grad.size() );
  
  CPPUNIT_ASSERT_DOUBLES_EQUAL( ana_val, num_val, 1e-6 );
  for (size_t i = 0; i < pd.num_free_vertices(); ++i) {
    CPPUNIT_ASSERT_VECTORS_EQUAL( num_grad[i], ana_grad[i], 1e-3 );
  }
}

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

bool PatchPowerMeanP::initialize_block_coordinate_descent( Mesh* mesh, 
                                                      MeshDomain* domain, 
                                                      const Settings* settings,
                                                      PatchSet* patch_set,
                                                      MsqError& err )
{
  clear();
  PatchIterator patches( patch_set );
  
  PatchData pd;
  pd.set_mesh( mesh );
  pd.set_domain( domain );
  if (settings)
    pd.attach_settings(settings);
  
  bool result = true;
  while (patches.get_next_patch( pd, err ) && !MSQ_CHKERR(err))
  {
    double value;
    bool b = evaluate( ObjectiveFunction::ACCUMULATE, pd, value, false, err ); 
    MSQ_ERRZERO(err);
    result = result && b;
  }
  return result;
}

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

bool AffineMapMetric::evaluate_with_indices( PatchData& pd,
                                             size_t handle,
                                             double& value,
                                             std::vector<size_t>& indices,
                                             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();
  const size_t* conn = elem.get_vertex_index_array();
  
    // this metric only supports sampling at corners
  if (s.dimension != 0) {
    if (type != TRIANGLE && type != TETRAHEDRON) {
      MSQ_SETERR(err)("Invalid sample point for AffineMapMetric", MsqError::UNSUPPORTED_ELEMENT );
      return false;
    }
    s.dimension = 0;
    s.number = 0;
  }

  unsigned n;
  const unsigned* adj = TopologyInfo::adjacent_vertices( type, s.number, n );
  indices.clear();
  if (conn[s.number] < pd.num_free_vertices())
    indices.push_back(conn[s.number]);
  for (unsigned i = 0; i < n; ++i)
    if (conn[adj[i]] < pd.num_free_vertices())
      indices.push_back(conn[adj[i]]);
  
  return evaluate( pd, handle, value, err );
}

bool CompositeOFAdd::evaluate_with_Hessian( EvalType type, 
                                            PatchData& pd,
                                            double& value_out,
                                            std::vector<Vector3D>& grad_out,
                                            MsqHessian& Hessian_out,
                                            MsqError& err )
{
  double value_2;
  bool ok;
  
  mHessian.initialize( Hessian_out );
  
  ok = objFunc1->evaluate_with_Hessian( type, pd, value_out, grad_out, Hessian_out, err );
  if (MSQ_CHKERR(err) || !ok) return false;
  ok = objFunc2->evaluate_with_Hessian( type, pd, value_2, mGradient, mHessian, err );
  if (MSQ_CHKERR(err) || !ok) return false;
  
  value_out += value_2;
  
  assert( grad_out.size() == pd.num_free_vertices() );
  assert( mGradient.size() == pd.num_free_vertices() );
  
  for (size_t i = 0; i < pd.num_free_vertices(); ++i) 
    grad_out[i] += mGradient[i];
  Hessian_out.add( mHessian );
  return true;
}

bool CompositeOFAdd::evaluate_with_gradient( EvalType type, 
                                             PatchData& pd,
                                             double& value_out,
                                             std::vector<Vector3D>& grad_out,
                                             MsqError& err )
{
  double value_2;
  bool ok;
  
  ok = objFunc1->evaluate_with_gradient( type, pd, value_out, grad_out, err );
  if (MSQ_CHKERR(err) || !ok) return false;
  ok = objFunc2->evaluate_with_gradient( type, pd, value_2, mGradient, err );
  if (MSQ_CHKERR(err) || !ok) return false;
  
  assert( grad_out.size() == pd.num_free_vertices() );
  assert( mGradient.size() == pd.num_free_vertices() );
  
  std::vector<Vector3D>::iterator i = grad_out.begin(), j = mGradient.begin();
  while (i != grad_out.end()) {
    *i += *j;
    ++i;
    ++j;
  }
  value_out += value_2;
  return true;
}

   /*! \fn create_six_quads_patch(PatchData &four_quads, MsqError &err)
     our 2D set up: 6 quads, 1 center vertex outcentered by (0,-0.5), the other centered
      7____6____5___11
      |    |    |    |
      | 2  |  3 | 5  |
      8-_  |  _-4---10       vertex 1 is at (0,0)
      |  -_0_-  |    |       vertex 11 is at (3,2)
      | 0  |  1 | 4  |
      1----2----3----9

      use destroy_patch_with_domain() in sync.
   */
   inline void create_six_quads_patch_with_domain(PatchData &pd, MsqError &err) 
   {
     // associates domain
     Vector3D pnt(0,0,0);
     Vector3D s_norm(0,0,3);
     pd.set_domain( new PlanarDomain(s_norm, pnt) );

     double coords[] = { 1,.5, 0,
                         0, 0, 0,
                         1, 0, 0,
                         2, 0, 0,
                         2, 1, 0,
                         2, 2, 0,
                         1, 2, 0,
                         0, 2, 0,
                         0, 1, 0,
                         3, 0, 0,
                         3, 1, 0,
                         3, 2, 0 };

     size_t indices[] = { 1,  2,  0, 8,
                          2,  3,  4, 0,
                          8,  0,  6, 7,
                          0,  4,  5, 6,
                          3,  9, 10, 4,
                          4, 10, 11, 5 };
     
     bool fixed[] = { false, true, true, true, 
                      false, true, true, true, 
                      true, true, true, true };
     
     pd.fill( 12, coords, 6, QUADRILATERAL, indices, fixed, err );
   }

bool EdgeLengthQualityMetric::evaluate_vertex(PatchData &pd, MsqVertex* vert,
                                             double &fval, MsqError &err)
{
  fval=0.0;
  size_t this_vert = pd.get_vertex_index(vert);
  size_t other_vert;
  vector<size_t> adj_verts;
  Vector3D edg;
  pd.get_adjacent_vertex_indices(this_vert,adj_verts,err);  MSQ_ERRZERO(err);
  int num_sample_points=adj_verts.size();
  double *metric_values=new double[num_sample_points];
  MsqVertex* verts = pd.get_vertex_array(err);  MSQ_ERRZERO(err);
  int point_counter=0;
  while(!adj_verts.empty()){
    other_vert=adj_verts.back();
    adj_verts.pop_back();
    edg[0]=verts[this_vert][0]-verts[other_vert][0];
    edg[1]=verts[this_vert][1]-verts[other_vert][1];
    edg[2]=verts[this_vert][2]-verts[other_vert][2];    
    metric_values[point_counter]=edg.length();
    ++point_counter;
  }
  fval=average_metrics(metric_values,num_sample_points,err);  MSQ_ERRZERO(err);
  delete[] metric_values;
  
  return true;
  
}

void write_gnuplot( Mesh* mesh, const char* out_filebase, MsqError &err)
{
    // loads a global patch
  PatchData pd;
  pd.set_mesh( mesh );
  pd.fill_global_patch( err ); MSQ_ERRRTN(err);
  write_gnuplot( pd, out_filebase, err );
}

/*! 
  Numerically Calculates the gradient of the ObjectiveFunction for the
  free vertices in the patch.  Returns 'false' if the patch is outside
  of a required feasible region, returns 'ture' otherwise.
  The behavior of the function depends on the value of the boolean
  useLocalGradient.  If useLocalGradient is set to
  'true', compute_numerical_gradient creates a sub-patch around a free
  vertex, and then perturbs that vertex in one of the coordinate directions.
  Only the ObjectiveFunction value on the local sub-patch is used in the
  computation of the gradient.  Therefore, useLocalGradient should only
  be set to 'true' for ObjectiveFunctions which can use this method.  Unless
  the concrete ObjectiveFunction sets useLocalGradient to 'true' in its
  constructor, the value will be 'false'.  In this case, the objective
  function value for the entire patch is used in the calculation of the
  gradient.  This is computationally expensive, but it is numerically
  correct for all (C_1) functions.
  \param pd  PatchData on which the gradient is taken.
  \param grad  Array of Vector3D of length the number of vertices used to store gradient.
  \param OF_val will be set to the objective function value.
 */
bool ObjectiveFunction::evaluate_with_gradient( EvalType eval_type,
                                                PatchData &pd,
                                                double& OF_val,
                                                std::vector<Vector3D>& grad,
                                                MsqError &err )
{
  bool b;
  grad.resize( pd.num_free_vertices() );
  
    // Fast path for single-free-vertex patch
  if (pd.num_free_vertices() == 1) {
    const EvalType sub_type = (eval_type == CALCULATE) ? CALCULATE : TEMPORARY;
    b = compute_subpatch_numerical_gradient( eval_type, sub_type, pd, OF_val, grad[0], err );
    return !MSQ_CHKERR(err) && b;
  }
  
  ObjectiveFunction* of = this;
  std::auto_ptr<ObjectiveFunction> deleter;
  if (eval_type == CALCULATE) {
    of->clear();
    b = of->evaluate( ACCUMULATE, pd, OF_val, OF_FREE_EVALS_ONLY, err );
    if (err) { // OF doesn't support BCD type evals, try slow method
      err.clear();
      of->clear();
      b = compute_patch_numerical_gradient( CALCULATE, CALCULATE, pd, OF_val, grad, err );
      return !MSQ_CHKERR(err) && b;
    }
    else if (!b)
      return b;
  } 
  else {
    b = this->evaluate( eval_type, pd, OF_val, OF_FREE_EVALS_ONLY, err );
    if (MSQ_CHKERR(err) || !b)
      return false;
    of = this->clone();
    deleter = std::auto_ptr<ObjectiveFunction>(of);
  }

    // Determine number of layers of adjacent elements based on metric type.
  unsigned layers = min_patch_layers();
  
    // Create a subpatch for each free vertex and use it to evaluate the
    // gradient for that vertex.
  double flocal;
  PatchData subpatch;
  for (size_t i = 0; i < pd.num_free_vertices(); ++i)
  {
    pd.get_subpatch( i, layers, subpatch, err ); MSQ_ERRZERO(err);
    b = of->compute_subpatch_numerical_gradient( SAVE, TEMPORARY, subpatch, flocal, grad[i], err );
    if (MSQ_CHKERR(err) || !b) {
      of->clear();
      return false;
    }
  }
  
  of->clear();
  return true;
}