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

void generate_minop_mesh(TriMesh & mesh,
                         const string & filename,
                         double edge_length){
  double angle = 0.125;
  cout << "Initial meshing..."<< endl;
  mesh.build_box_boundary({{-1.1,1.1},{-1.1,1.1}});
  mesh.build_circle(zeros<vec>(2),50,1.0);
  mesh.build_circle(zeros<vec>(2),30,1.0/sqrt(2.0));
  mesh.build_circle(zeros<vec>(2),25,0.25);

  cout << "Refining based on (" << angle
       << "," << edge_length <<  ") criterion ..."<< endl;
  mesh.refine(angle,edge_length);
  
  //cout << "Optimizing (25 rounds of Lloyd)..."<< endl;
  //mesh.lloyd(15);
  mesh.freeze();

  // Write initial mesh to file
  cout << "Writing:"
       << "\n\t" << (filename + ".node") << " (Shewchuk node file)"
       << "\n\t" << (filename + ".ele") << " (Shewchuk element file)"
       << "\n\t" << (filename + ".tri") << " (CGAL mesh file)" << endl;
  mesh.write_shewchuk(filename);
  mesh.write_cgal(filename + ".tri");
}

////////////////////////////////////////////////////////////
// MAIN FUNCTION ///////////////////////////////////////////
////////////////////////////////////////////////////////////
int main(int argc, char** argv)
{
  po::variables_map var_map = read_command_line(argc,argv);
  arma_rng::set_seed_random();

  // Read in the CGAL mesh
  TriMesh mesh;
  string file_base = var_map["outfile_base"].as<string>();
  double edge_length = var_map["edge_length"].as<double>();
  generate_minop_mesh(mesh,file_base,edge_length);
  mesh.freeze();

  // Stats
  uint N = mesh.number_of_vertices();
  uint F = mesh.number_of_faces();
  cout << "Mesh stats:"
       << "\n\tNumber of vertices: " << N
       << "\n\tNumber of faces: " << F
       << endl;

  // Build value basis
  uint num_value_basis = 25;
  uint num_flow_basis = 10;
  Points points = mesh.get_spatial_nodes();
  cout << "Generating Radial Fourier basis for value..." << endl;
  mat value_basis = make_radial_fourier_basis(points,
                                              num_value_basis,
                                              (double)num_value_basis);
  cout << "\tOrthogonalizing..." << endl;
  value_basis = orth(value_basis);

  cout << "Generating Voronoi basis for flow..." << endl;
  sp_mat sp_value_basis = sp_mat(value_basis);
  
  Points centers = 2 * randu(10,2) - 1;
  mat flow_basis = make_voronoi_basis(points,
                                      centers);
  cout << "\tOrthogonalizing..." << endl;
  flow_basis = orth(flow_basis);
  sp_mat sp_flow_basis = sp_mat(flow_basis);

  cout << "Building LCP..." << endl;
  vec ref_weights = ones<vec>(N) / (double)N;
  LCP ref_lcp;
  vec ans;
  build_minop_lcp(mesh,ref_weights,ref_lcp,ans);
  assert(N == ans.n_elem);

  cout << "Building PLCP..." << endl;
  block_sp_vec D = {sp_value_basis,
                    sp_flow_basis,
                    sp_flow_basis};  
  sp_mat P = block_diag(D);
  double regularizer = 1e-12;
  sp_mat U = P.t() * (ref_lcp.M + regularizer*speye(size(ref_lcp.M)));
  vec q =  P *(P.t() * ref_lcp.q);
  assert(3*N == P.n_rows);
  assert(3*N == q.n_rows);

  bvec free_vars = zeros<bvec>(3*N);
  free_vars.head(N).fill(1);
  
  PLCP ref_plcp = PLCP(P,U,q,free_vars);
  ref_plcp.write(file_base + ".plcp");
  ProjectiveSolver psolver;
  psolver.comp_thresh = 1e-12;
  psolver.max_iter = 250;
  psolver.aug_rel_scale = 5;
  psolver.verbose = false;
  psolver.initial_sigma = 0.3;

  cout << "Starting reference solve..." << endl;
  SolverResult ref_sol = psolver.aug_solve(ref_plcp);
  cout << "\tDone." << endl;
  cout << "Reference solution error: "
       << norm(ans - ref_sol.p.head(N)) << endl;
  assert(ALMOST_ZERO > norm(ref_sol.d.head(N))); // Essentially zero
  ref_sol.write(file_base + ".sol");
  
  psolver.comp_thresh = 1e-8;
  // Exactish
  vec twiddle = vec(N);
  for(uint i = 0; i < N; i++){
    cout << "Component: " << i << endl;
    LCP twiddle_lcp;
    vec et = zeros<vec>(N);
    et(i) += 1.0 / (double) N;
    assert(size(ref_weights) == size(et));
    build_minop_lcp(mesh,ref_weights + et,twiddle_lcp,ans);
    vec twiddle_q =  P *(P.t() * twiddle_lcp.q);

    PLCP twiddle_plcp = PLCP(P,U,twiddle_q,free_vars);
    SolverResult twiddle_sol = psolver.aug_solve(twiddle_plcp);
    twiddle(i) = twiddle_sol.p(i) - ref_sol.p(i);
  }

  uint R = 75;
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
  po::variables_map var_map = read_command_line(argc,argv);

  string mesh_file = var_map["infile_base"].as<string>() + ".tri";

  // Read in the CGAL mesh
  TriMesh mesh;
  cout << "Reading in cgal mesh file [" << mesh_file << ']'  << endl;
  mesh.read_cgal(mesh_file);
  mesh.freeze();
  uint V = mesh.number_of_vertices();
  uint N = mesh.number_of_all_nodes();
  assert(N == V+1);
  uint F = mesh.number_of_faces();
  cout << "Mesh stats:"
       << "\n\tNumber of vertices: " << V
       << "\n\tNumber of faces: " << F
       << endl;
  Points centers = mesh.get_cell_centers();

  // Find boundary from the mesh and create the simulator object
  mat bbox = mesh.find_bounding_box();
  vec lb = bbox.col(0);
  vec ub = bbox.col(1);
  cout << "\tLower bound:" << lb.t()
       << "\tUpper bound:" << ub.t();

  // Read in solution information
  string soln_file = var_map["infile_base"].as<string>() + ".sol";
  cout << "Reading in LCP solution file [" << soln_file << ']'  << endl;
  Unarchiver sol_unarch(soln_file);
  vec p = sol_unarch.load_vec("p");

  string lcp_file = var_map["infile_base"].as<string>() + ".lcp";
  Unarchiver lcp_unarch(lcp_file);
  vec q = lcp_unarch.load_vec("q");

  Archiver arch;
  
  // 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;
//.........这里部分代码省略.........