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

   void test_elements()
   {
      MsqPrintError err(cout);
     /* Adds TSTT mesh to a MeshSet. */
      MeshSet mesh_set;
      TSTTMesh tstt_mesh;
      tstt_mesh.set_mesh(tri10);
      mesh_set.add_mesh(&tstt_mesh, err); CPPUNIT_ASSERT(!err);

      /* Retrieves a global patch */
      PatchData pd;
      PatchDataParameters pd_params;
      pd_params.set_patch_type(PatchData::ELEMENTS_ON_VERTEX_PATCH, err, 1, 0);
      
      mesh_set.get_next_patch(pd, pd_params, err); CPPUNIT_ASSERT(!err);

      int free_vtx = pd.num_free_vertices(err); CPPUNIT_ASSERT(!err);
      std::cout << "nb of free vertices: " << free_vtx << std::endl;
      CPPUNIT_ASSERT( free_vtx == 1 );
      
      element_array =  pd.get_element_array(err); CPPUNIT_ASSERT(!err);
      num_elements = pd.num_elements();
      CPPUNIT_ASSERT( num_elements == 6 );

      //      for (int i=0; i<num_elements; ++i) {
      //         std::cout << element_array[i];
      //      }
      
      vtx_array = pd.get_vertex_array(err); CPPUNIT_ASSERT(!err);
      num_vertices = pd.num_vertices();
      CPPUNIT_ASSERT( num_vertices == 7 );

      //      for (int i=0; i<num_vertices; ++i) {
      //         std::cout << vtx_array[i];
      //      }

      CPPUNIT_ASSERT( tri_check_validity() == 1 );
      
      mesh_set.get_next_patch(pd, pd_params, err); CPPUNIT_ASSERT(!err);

      element_array =  pd.get_element_array(err); CPPUNIT_ASSERT(!err);
      num_elements = pd.num_elements();
      CPPUNIT_ASSERT( num_elements == 6 );

      //      for (int i=0; i<num_elements; ++i) {
      //         std::cout << element_array[i];
      //      }
      
      vtx_array = pd.get_vertex_array(err); CPPUNIT_ASSERT(!err);
      num_vertices = pd.num_vertices();
      CPPUNIT_ASSERT( num_vertices == 7 );

      //      for (int i=0; i<num_vertices; ++i) {
      //         std::cout << vtx_array[i];
      //      }

      CPPUNIT_ASSERT( tri_check_validity() == 1 );
   }

void Randomize::optimize_vertex_positions(PatchData &pd, 
                                                MsqError &err)
{
    //cout << "- Executing Randomize::optimize_vertex_position()\n";

  int num_local_vertices = pd.num_vertices();
  // gets the array of coordinates for the patch and print it 
  MsqVertex *patch_coords = pd.get_vertex_array(err); MSQ_ERRRTN(err);
  // does the randomize smooth
  MsqFreeVertexIndexIterator free_iter(&pd, err); MSQ_ERRRTN(err);
  free_iter.reset();
  free_iter.next();
    //find the free vertex.
  int m=free_iter.value();
  randomize_vertex(pd, num_local_vertices,
                   patch_coords[m], err);  MSQ_ERRRTN(err);
  pd.snap_vertex_to_domain(m,err); MSQ_ERRRTN(err);
}

/*! The type of targets computed by this function is selected by the constructor of
    the base classes. */
void WTargetCalculator::compute_target_matrices(PatchData &pd, MsqError &err)
{
  MSQ_FUNCTION_TIMER( "WTargetCalculator::compute_target_matrice" );

  size_t num_elements=pd.num_elements();
  
  PatchData ref_pd, *ref_pd_ptr;
  if ( refMesh != 0 ) {
    // If there is a reference mesh, gets a patch ref_pd equivalent to the patch pd of the main mesh.
    PatchDataParameters ref_pd_params(this->get_all_parameters());
    refMesh->get_next_patch(ref_pd, ref_pd_params, err); MSQ_ERRRTN(err);
    // Make sure topology of ref_pd and pd are equal
    assert( num_elements == ref_pd.num_elements() );
    size_t num_vertices=pd.num_vertices();
    assert( num_vertices == ref_pd.num_vertices() );
    ref_pd_ptr = &ref_pd;
  }
  else {
    // the reference patch is the same as the working patch if there is no reference mesh.
    ref_pd_ptr = &pd;
  }
  
  MsqMeshEntity* elems = pd.get_element_array(err); MSQ_ERRRTN(err);
  MsqMeshEntity* elems_ref = ref_pd_ptr->get_element_array(err); MSQ_ERRRTN(err);

  Matrix3D W_guides[MSQ_MAX_NUM_VERT_PER_ENT];
  TargetMatrix matrices[MSQ_MAX_NUM_VERT_PER_ENT];
  
  for (size_t i=0; i<num_elements; ++i) {
    int nve = elems[i].vertex_count();
    assert( nve = elems_ref[i].vertex_count() );

    compute_guide_matrices(guideMatrix, *ref_pd_ptr, i, W_guides, nve, err); MSQ_ERRRTN(err);
    for (int c = 0; c < nve; ++c)
      matrices[c] = W_guides[c];
    pd.targetMatrices.set_element_corner_tags( &pd, i, matrices, err ); MSQ_ERRRTN(err);
  }

  //if ( refMesh != 0 ) delete ref_pd;
}

bool LPtoPTemplate::concrete_evaluate(PatchData &pd, double &fval,
                                      MsqError &err){
  size_t index=0;
  MsqMeshEntity* elems=pd.get_element_array(err);
  bool obj_bool=true;
    //double check for pVal=0;
  if(pVal==0){
    MSQ_SETERR(err)("pVal equal zero not allowed.  L_0 is not a valid norm.",
                    MsqError::INVALID_STATE);
    return false;
  }
  
    //Michael:  this may not do what we want
    //Set currentQM to be the first quality metric* in the list
  QualityMetric* currentQM = get_quality_metric();
  if(currentQM==NULL)
    currentQM=get_quality_metric_list().front();
  if(currentQM==NULL) {
    MSQ_SETERR(err)("NULL QualityMetric pointer in LPtoPTemplate",
                    MsqError::INVALID_STATE);
    return false;
  }
  size_t num_elements=pd.num_elements();
  size_t num_vertices=pd.num_vertices();
  size_t total_num=0;
  if(currentQM->get_metric_type()==QualityMetric::ELEMENT_BASED)   
    total_num=num_elements;
  else if (currentQM->get_metric_type()==QualityMetric::VERTEX_BASED)
    total_num=num_vertices;
  else {
    MSQ_SETERR(err)("Make sure MetricType is initialised in concrete "
                    "QualityMetric constructor.", MsqError::INVALID_STATE);
    return false;
  }
  
  msq_std::vector<double> metric_values(total_num);
  if(currentQM->get_metric_type()==QualityMetric::ELEMENT_BASED)
  {
    for (index=0; index<num_elements;index++)
    {
        //if invalid return false after clean-up
      obj_bool=currentQM->evaluate_element(pd, (&elems[index]),
                                           metric_values[index], err);
      if(MSQ_CHKERR(err) || !obj_bool){
        fval=0.0;
        return false;
      }
      
      metric_values[index]=fabs(metric_values[index]);
      MSQ_DBGOUT(3) << "      o  Quality metric value for element "
                    << index << "\t: " << metric_values[index] << "\n";
    }
  }
  else if(currentQM->get_metric_type()==QualityMetric::VERTEX_BASED)
  {
    MsqVertex* vertices=pd.get_vertex_array(err); MSQ_ERRZERO(err);
    for (index=0; index<num_vertices;index++)
    {
        //evaluate metric for this vertex
      obj_bool=currentQM->evaluate_vertex(pd, (&vertices[index]),
                                          metric_values[index], err);
        //if invalid return false after clean-up
      if(MSQ_CHKERR(err) || !obj_bool){
        fval=0.0;
        return false;
      }
      
      metric_values[index]=fabs(metric_values[index]);
    }
  }
  fval=compute_function(&metric_values[0], total_num, err);
  return !MSQ_CHKERR(err);
}

/*! \fn LPtoPTemplate::compute_analytical_hessian(PatchData &pd, MsqHessian &hessian, MsqError &err)

    For each element, each entry to be accumulated in the Hessian for
    this objective function (\f$ \sum_{e \in E} Q(e)^p \f$ where \f$ E \f$
    is the set of all elements in the patch) has the form:
    \f$ pQ(e)^{p-1} \nabla^2 Q(e) + p(p-1)Q(e)^{p-2} \nabla Q(e) [\nabla Q(e)]^T \f$.

    For \f$ p=2 \f$, this simplifies to
    \f$ 2Q(e) \nabla^2 Q(e) + 2 \nabla Q(e) [\nabla Q(e)]^T \f$.

    For \f$ p=1 \f$, this simplifies to \f$ \nabla^2 Q(e) \f$.

    The \f$ p=1 \f$ simplified version is implemented directly
    to speed up computation. 

    This function does not support vertex-based metrics.
    
    \param pd The PatchData object for which the objective function
           hessian is computed.
    \param hessian: this object must have been previously initialized.
*/
bool LPtoPTemplate::compute_analytical_hessian(PatchData &pd,
                                               MsqHessian &hessian,
                                               Vector3D *const &grad,
                                               double &OF_val,
                                               MsqError &err)
{
  double scaling_value=1.0;
  
  MSQ_FUNCTION_TIMER( "LPtoPTemplate::compute_analytical_hessian" );

  MsqMeshEntity* elements = pd.get_element_array(err); MSQ_ERRZERO(err);
  MsqVertex* vertices = pd.get_vertex_array(err); MSQ_ERRZERO(err);
  size_t num_elems = pd.num_elements();
    //if scaling divide by the number of elements.
  if(dividingByN){
    if(num_elems<=0) {
      MSQ_SETERR(err)("LPtoP is attempting to divide by zero in analytical Hessian.",
                      MsqError::INVALID_MESH);
      return false;
    }
    scaling_value/=num_elems;
  }
  
  size_t num_vertices = pd.num_vertices();
  Matrix3D elem_hessian[MSQ_MAX_NUM_VERT_PER_ENT*(MSQ_MAX_NUM_VERT_PER_ENT+1)/2];
  Matrix3D elem_outer_product;
  Vector3D grad_vec[MSQ_MAX_NUM_VERT_PER_ENT];
  double QM_val;
  double fac1, fac2;
  Matrix3D grad_outprod;
  bool qm_bool;
  QualityMetric* currentQM = get_quality_metric();
  
  MsqVertex* ele_free_vtces[MSQ_MAX_NUM_VERT_PER_ENT];
  short i;
  for (i=0; i<MSQ_MAX_NUM_VERT_PER_ENT; ++i) ele_free_vtces[i]=NULL;
  
  const size_t* vtx_indices;
    
  size_t e, v;
  size_t nfve; // number of free vertices in element
  short j,n;

  hessian.zero_out();
  for (v=0; v<num_vertices; ++v) grad[v] = 0.;
  OF_val = 0.;
  
  // Loops over all elements in the patch.
  for (e=0; e<num_elems; ++e) {
    short nve = elements[e].vertex_count();
    
    // Gets a list of free vertices in the element.
    vtx_indices = elements[e].get_vertex_index_array();
    nfve=0;
    for (i=0; i<nve; ++i) {
      if ( vertices[vtx_indices[i]].is_free_vertex() ) {
        ele_free_vtces[nfve] = vertices + vtx_indices[i];
        ++nfve;
      }
    }
    
    // Computes \nabla^2 Q(e). Only the free vertices will have non-zero entries. 
    qm_bool = currentQM->compute_element_hessian(pd,
                                    elements+e, ele_free_vtces,
                                    grad_vec, elem_hessian,
                                    nfve, QM_val, err);
    if (MSQ_CHKERR(err) || !qm_bool) return false;


    // **** Computes Hessian ****
    double QM_pow=1.;
    if (pVal == 1) {
      n=0;
      for (i=0; i<nve; ++i) {
        for (j=i; j<nve; ++j) {
            //negate if necessary
          elem_hessian[n] *= (scaling_value * get_negate_flag());
          ++n;
        }
//.........这里部分代码省略.........

/* virtual function reimplemented from QualityMetric. No doxygen doc needed. */
bool LPtoPTemplate::compute_analytical_gradient(PatchData &pd,
                                              Vector3D *const &grad,
                                              double &OF_val,
                                              MsqError &err, size_t array_size)
{
  MSQ_FUNCTION_TIMER( "LPtoPTemplate::compute_analytical_gradient" );
  
    //initialize the scaling value
  double scaling_value=1.0;
 
  size_t num_elements=pd.num_elements();
  size_t num_vertices=pd.num_vertices();
  if( num_vertices!=array_size && array_size>0)
  {
    MSQ_SETERR(err)("Incorrect array size.", MsqError::INVALID_ARG);
    return false;
  }
  
  MsqMeshEntity* elems=pd.get_element_array(err); MSQ_ERRZERO(err);
  MsqVertex* vertices=pd.get_vertex_array(err); MSQ_ERRZERO(err);
  bool qm_bool=true;
  double QM_val;
  OF_val = 0.;
  size_t i;
  int p1;
  
  //Set currentQM to be quality metric (possibly composite) associated with the objective function
  QualityMetric* currentQM = get_quality_metric();
  if(currentQM==NULL) {
    MSQ_SETERR(err)("LPtoPTemplate has NULL QualityMetric pointer.",MsqError::INVALID_STATE);
    return false;
  }
  enum QualityMetric::MetricType qm_type=currentQM->get_metric_type();
  
  if (qm_type!=QualityMetric::ELEMENT_BASED &&
      qm_type!=QualityMetric::VERTEX_BASED) {
    MSQ_SETERR(err)("Make sure MetricType is initialised"
                    "in concrete QualityMetric constructor.",
                    MsqError::INVALID_STATE);
    return false;
  }


  // zeros out objective function gradient
  for (i=0; i<num_vertices; ++i)
    grad[i] =0;
  
  // Computes objective function gradient for an element based metric
  if(qm_type==QualityMetric::ELEMENT_BASED){
      //if scaling, divid by num_elements
    if(dividingByN){
      if(num_elements<=0) {
        MSQ_SETERR(err)("The number of elements should not be zero.",MsqError::INVALID_MESH);
        return false;
      }
      
      scaling_value/=num_elements;
    }
    size_t e, ve;
    size_t nfve; // num free vtx in element
    size_t nve; // num vtx in element
    MsqVertex* ele_free_vtces[MSQ_MAX_NUM_VERT_PER_ENT];
    const size_t *ele_vtces_ind;

    // loops over all elements.
    for (e=0; e<num_elements; ++e) {
      // stores the pointers to the free vertices within the element
      // (using pointer arithmetic).
      nfve = 0;
      nve = elems[e].vertex_count();
      ele_vtces_ind = elems[e].get_vertex_index_array();
      for (ve=0; ve<nve; ++ve) {
        if (vertices[ele_vtces_ind[ve]].is_free_vertex()) {
          ele_free_vtces[nfve] = vertices + ele_vtces_ind[ve];
          ++nfve;
        }
      }

      // Computes q and grad(q)
      Vector3D grad_vec[MSQ_MAX_NUM_VERT_PER_ENT];
      qm_bool = currentQM->compute_element_gradient(
                                     pd, &elems[e],
                                     ele_free_vtces,
                                     grad_vec, nfve, QM_val, err);
      if(MSQ_CHKERR(err) || !qm_bool) return false;

      // computes p*|Q(e)|^{p-1}
      QM_val = fabs(QM_val);
      double QM_pow=1.0;
      double factor;
      if (pVal==1) factor=1;
      else {
        QM_pow=QM_val;
        for (p1=1; p1<pVal-1; ++p1)
          QM_pow*=QM_val;
        factor = QM_pow * pVal;
      }
        //this scales the gradient
      factor *= (scaling_value * get_negate_flag());
//.........这里部分代码省略.........