C++ TriMesh::center_of_face方法代码示例

//.........这里部分代码省略.........
  
  // Make sure that the primal information makes sense
  assert(0 == p.n_elem % V);
  uint A = p.n_elem / V;
  assert(A == 3);
  cout << "Blocking primal solution..."
       << "\n\tLength of primal solution: " << p.n_elem
       << "\n\tRatio of primal length to vertex number: " << A << endl;
  mat P = reshape(p,size(V,A));
  P = join_vert(P,datum::inf * ones<rowvec>(3)); // Pad
  vec value = P.col(0);
  mat flows = P.tail_cols(2);

  mat Q = reshape(q,size(V,A));
  Q = join_vert(Q,datum::inf * ones<rowvec>(3)); // Pad
  vec recon_b = mesh.interpolate(centers,
                                 conv_to<vec>::from(Q.col(1)));
  vec recon_c = mesh.interpolate(centers,
                                 conv_to<vec>::from(Q.col(2)));
  arch.add_vec("recon_b",recon_b);
  arch.add_vec("recon_c",recon_c);


  vec area = mesh.cell_area();
  arch.add_vec("area",area);

  // True values
  vec sq_dist = sum(pow(centers,2),1);
  vec b = sq_dist;
  vec c = max(zeros<vec>(F),1 - sq_dist);
  vec x = arma::min(b,c);
  assert(all(x >= 0));
  assert(F == x.n_elem);
  
  uvec pi = arma::index_min(join_horiz(b,c),1);
  assert(F == pi.n_elem);
  arch.add_uvec("pi",pi);

  // Approx policy
  assert(2 == flows.n_cols);
  mat interp_flows = mesh.interpolate(centers,flows);
  uvec flow_pi = arma::index_max(interp_flows,1);
  arch.add_uvec("flow_pi",flow_pi);

  assert(F == flow_pi.n_elem);
  uvec diff = zeros<uvec>(F);
  diff(find(flow_pi != pi)).fill(1);
  arch.add_uvec("policy_diff",diff);

  // Approx value
  vec interp_value = mesh.interpolate(centers,value);
  assert(F == interp_value.n_elem);
  vec res = abs(x - interp_value);
  arch.add_vec("residual",res);

  vec heuristic = res;
  heuristic(find(flow_pi != pi)) *= 4;
  arch.add_vec("heuristic",heuristic);

  double quant = 0.9;
  cout << "Quantile:" << quant << endl;
  double cutoff = quantile(heuristic,quant);
  cout << "\t Cutoff: " << cutoff
       << "\n\t Max:" << max(heuristic)
       << "\n\t Min:" << min(heuristic) << endl;

  // Split the cells if they have a large heuristic_1 or
  // policies disagree on them.
  TriMesh new_mesh(mesh);  
  Point center;
  uint new_nodes = 0;
  for(uint f = 0; f < F; f++){
    if(area(f) < 2e-4){
      continue;
    }
    if(new_nodes > 200) break;
    if(heuristic(f) > cutoff){
      center = convert(mesh.center_of_face(f));
      new_mesh.insert(center);
      new_nodes++;
    }
  }
  cout << "Added " << new_nodes << " new nodes..." << endl;
  
  cout << "Refining..." << endl;
  new_mesh.refine(0.125,1.0);
  new_mesh.lloyd(10);
  //new_mesh.insert(Point(0,0));
  new_mesh.freeze();

  // Write out all the information
  string out_file_base = var_map["outfile_base"].as<string>();
  arch.write(out_file_base + ".stats");
  cout << "Writing..."
       << "\n\tCGAL mesh file: " << out_file_base << ".tri"
       << "\n\tShewchuk node file: " << out_file_base << ".node"
       << "\n\tShewchuk ele file: " << out_file_base << ".ele" << endl;
  new_mesh.write_cgal(out_file_base + ".tri");
  new_mesh.write_shewchuk(out_file_base);
}