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

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

/*! \todo Michael:  optimize_vertex_position is probably not implemented
  in an optimal way.  We used to use all of the vertices in
  the patch as 'adjacent' vertices.  Now we call get_adjacent_vertex_indices.
  We could use a VERTICES_ON_VERTEX type of patch or a global patch?
*/
void SmartLaplacianSmoother::optimize_vertex_positions(PatchData &pd, 
                                                MsqError &err)
{
    //default the laplacian smoother to 3 even for 2-d elements.
    //int dim = get_mesh_set()->space_dim();
  size_t dim = 3;
  MsqVertex* verts=pd.get_vertex_array(err);  MSQ_ERRRTN(err);
    //the objective function evaluator
  OFEvaluator& obj_func = get_objective_function_evaluator();
    //variables for the function values.
  double orig_val=0;
  double mod_val=0;
    //compute the original function value and check validity
  bool valid_flag = obj_func.evaluate(pd,orig_val,err);  MSQ_ERRRTN(err);
  // does the Laplacian smoothing
  MsqFreeVertexIndexIterator free_iter(pd, err);  MSQ_ERRRTN(err);
  free_iter.reset();
  free_iter.next();
    //m is the free vertex.
  size_t m=free_iter.value();
  vector<size_t> vert_indices;
  vert_indices.reserve(25);
    //get vertices adjacent to vertex m
  pd.get_adjacent_vertex_indices(m,vert_indices,err);  MSQ_ERRRTN(err);
    //move vertex m
    //save the original position of the free vertex
  Vector3D orig_position(verts[m]);
    //smooth the patch
  centroid_smooth_mesh(pd, vert_indices.size(), vert_indices,
                       m, dim, err); MSQ_ERRRTN(err);
    //snap vertex m to domain
  pd.snap_vertex_to_domain(m,err);  MSQ_ERRRTN(err);
    //if the original function val was invalid, then we allow the move
    //But, if it wasn valid, we need to decide.
  if(valid_flag){
      //compute the new value
    valid_flag = obj_func.evaluate(pd,mod_val,err);  MSQ_ERRRTN(err);
      //if the new value is worse the original OR if the new value is not
      //valid (we already know the original value was valid by above) then
      //we don't allow the move.
    if(!valid_flag || mod_val>orig_val){
        //move the vert back to where it was.
      verts[m]=orig_position;
        //PRINT_INFO("\norig = %f, new = %f, new valid = %d",orig_val,mod_val,valid_flag);
    }
    
  }
  
}

/*!For the given vertex, vert, with connected edges of lengths l_j for
  j=1...k, the metric value is the average (where the default average
  type is SUM) of
        u_j = ( | l_j - lowVal | - (l_j - lowVal) )^2 +
              ( | highVal - l_j | - (highVal - l_j) )^2.
*/
bool EdgeLengthRangeQualityMetric::evaluate_common(PatchData &pd, 
                                             size_t this_vert,
                                             double &fval, 
                                             std::vector<size_t>& adj_verts,
                                             MsqError &err)
{
  fval=0.0;
  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];
  const MsqVertex* verts = pd.get_vertex_array(err);  MSQ_ERRZERO(err);
    //store the length of the edge, and the first and second component of
    //metric values, respectively.
  double temp_length=0.0;
  double temp_first=0.0;
  double temp_second=0.0;
    //PRINT_INFO("INSIDE ELR, vertex = %f,%f,%f\n",verts[this_vert][0],verts[this_vert][1],verts[this_vert][2]);
    //loop while there are still more adjacent vertices.
  for (unsigned i = 0; i < adj_verts.size(); ++i) 
  {
    edg = verts[this_vert] - verts[adj_verts[i]];
      //compute the edge length
    temp_length=edg.length();
      //get the first component
    temp_first = temp_length - lowVal;
    temp_first = fabs(temp_first) - (temp_first);
    temp_first*=temp_first;
      //get the second component
    temp_second = highVal - temp_length;
    temp_second = fabs(temp_second) - (temp_second);
    temp_second*=temp_second;
      //combine the two components
    metric_values[i]=temp_first+temp_second;
  }
    //average the metric values of the edges
  fval=average_metrics(metric_values,num_sample_points,err);
    //clean up
  delete[] metric_values;
    //always return true because mesh is always valid wrt this metric.
  return !MSQ_CHKERR(err);
  
}

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

bool LocalSizeQualityMetric::evaluate_with_indices( PatchData& pd,
                                                    size_t vertex,
                                                    double& value,
                                                    std::vector<size_t>& indices,
                                                    MsqError& err )
{
  indices.clear();
  pd.get_adjacent_vertex_indices( vertex, indices, err ); MSQ_ERRZERO(err);
  
  std::vector<size_t>::iterator r, w;
  for (r = w = indices.begin(); r != indices.end(); ++r) {
    if (*r < pd.num_free_vertices()) {
      *w = *r;
      ++w;
    }
  }
  indices.erase( w, indices.end() );
  if (vertex < pd.num_free_vertices())
    indices.push_back( vertex );
  
  bool rval = evaluate( pd, vertex, value, err );
  return !MSQ_CHKERR(err) && rval;
}

//.........这里部分代码省略.........
        // ... computes p*q^{p-1}*grad(q) ...
        grad_vec[i] *= factor;
        // ... and accumulates it in the objective function gradient.
        grad[pd.get_vertex_index(ele_free_vtces[i])] += grad_vec[i];
      }
    }
  }
  
  // Computes objective function gradient for a vertex based metric  
  else if (qm_type==QualityMetric::VERTEX_BASED){
      //if scaling, divide by the number of vertices
    if(dividingByN){
      if(num_elements<=0) {
        MSQ_SETERR(err)("The number of vertices should not be zero.",MsqError::INVALID_MESH);
        return false;
      }
      
      scaling_value/=num_vertices;
    }
    //vector for storing indices of vertex's connected elems
    msq_std::vector<size_t> vert_on_vert_ind;
    //position in pd's vertex array
    size_t vert_count=0;
    //position in vertex array
    size_t vert_pos=0;
    //loop over the free vertex indices to find the gradient...
    size_t vfv_array_length=10;//holds the current legth of vert_free_vtces
    msq_std::vector<MsqVertex*> vert_free_vtces(vfv_array_length);
    msq_std::vector<Vector3D> grad_vec(vfv_array_length); 
    for(vert_count=0; vert_count<num_vertices; ++vert_count){
      //For now we compute the metric for attached vertices and this
      //vertex, the above line gives us the attached vertices.  Now,
      //we must add this vertex.
      pd.get_adjacent_vertex_indices(vert_count,
                                     vert_on_vert_ind,err);
      vert_on_vert_ind.push_back(vert_count);
      size_t vert_num_vtces = vert_on_vert_ind.size();
      
      // dynamic memory management if arrays are too small.
      if(vert_num_vtces > vfv_array_length){
        vfv_array_length=vert_num_vtces+5;
        vert_free_vtces.resize(vfv_array_length);
        grad_vec.resize(vfv_array_length);
      }
      
      size_t vert_num_free_vtces=0;
      //loop over the vertices connected to this one (vert_count)
      //and put the free ones into vert_free_vtces
      while(!vert_on_vert_ind.empty()){
        vert_pos=(vert_on_vert_ind.back());
        //clear the vector as we go
        vert_on_vert_ind.pop_back();
        //if the vertex is free, add it to ver_free_vtces
        if(vertices[vert_pos].is_free_vertex()){
          vert_free_vtces[vert_num_free_vtces]=&vertices[vert_pos];
          ++vert_num_free_vtces ;
        }
      }
      
      qm_bool=currentQM->compute_vertex_gradient(pd,
                                                 vertices[vert_count],
                                                 &vert_free_vtces[0],
                                                 &grad_vec[0],
                                                 vert_num_free_vtces,
                                                 QM_val, err);
      if(MSQ_CHKERR(err) || !qm_bool){