C++ MeshType::edge_end方法代码示例

 std::pair<double, double> operator()(MeshType& mesh) const {
     double max_height = derived().initialize_node_values(mesh);
     
     double min_edge_length = (*mesh.edge_begin()).length();
     for (auto eit = mesh.edge_begin(); eit != mesh.edge_end(); ++eit) {
         min_edge_length = std::min(min_edge_length, (*eit).length());
     }
 
     for (auto tri_it = mesh.triangle_begin(); tri_it != mesh.triangle_end(); ++tri_it) {
         (*tri_it).value() = ((*tri_it).node(0).value() + (*tri_it).node(1).value() + (*tri_it).node(2).value()) / 3.0;
     }
     return std::make_pair(max_height, min_edge_length);
 }

int main(int argc, char* argv[])
{
  // Check arguments
  if (argc < 3) {
    std::cerr << "Usage: shallow_water NODES_FILE TRIS_FILE\n";
    exit(1);
  }

  MeshType mesh;
  // HW4B: Need node_type before this can be used!
#if 1
  std::vector<typename MeshType::node_type> mesh_node;
#endif

  // Read all Points and add them to the Mesh
  std::ifstream nodes_file(argv[1]);
  Point p;
  while (CS207::getline_parsed(nodes_file, p)) {
    // HW4B: Need to implement add_node before this can be used!
#if 1
    mesh_node.push_back(mesh.add_node(p));
#endif
  }

  // Read all mesh triangles and add them to the Mesh
  std::ifstream tris_file(argv[2]);
  std::array<int,3> t;
  while (CS207::getline_parsed(tris_file, t)) {
    // HW4B: Need to implement add_triangle before this can be used!
#if 1
    mesh.add_triangle(mesh_node[t[0]], mesh_node[t[1]], mesh_node[t[2]]);
#endif
  }

  // Print out the stats
  std::cout << mesh.num_nodes() << " "
            << mesh.num_edges() << " "
            << mesh.num_triangles() << std::endl;

  // HW4B Initialization
  // Set the initial conditions based off third argument
  // Perform any needed precomputation
  std::pair<double, double> value_pair;
  if ((*argv[3]) == '0') {
      std::cout << "Pebble Ripple" << std::endl;
      value_pair = PebbleRipple()(mesh);
  } else if ((*argv[3]) == '1') {
      std::cout << "Sharp Wave" << std::endl;
      value_pair = SharpWave()(mesh);
  } else {
      std::cout << "Dam Break" << std::endl;
      value_pair = DamBreak()(mesh);
  }
  
  double max_height = value_pair.first;
  double min_edge_length = value_pair.second;

  // Launch the SDLViewer
  CS207::SDLViewer viewer;
  viewer.launch();

  // HW4B: Need to define Mesh::node_type and node/edge iterator
  // before these can be used!
#if 1
  auto node_map = viewer.empty_node_map(mesh);
  viewer.add_nodes(mesh.node_begin(), mesh.node_end(),
                   CS207::DefaultColor(), NodePosition(), node_map);
  viewer.add_edges(mesh.edge_begin(), mesh.edge_end(), node_map);
#endif
  viewer.center_view();
  // CFL stability condition requires dt <= dx / max|velocity|
  // For the shallow water equations with u = v = 0 initial conditions
  //   we can compute the minimum edge length and maximum original water height
  //   to set the time-step
  // Compute the minimum edge length and maximum water height for computing dt
#if 1
  double dt = 0.25 * min_edge_length / (sqrt(grav * max_height));
#else
  // Placeholder!! Delete me when min_edge_length and max_height can be computed!
  double dt = 0.1;
#endif
  double t_start = 0;
  double t_end = 5;

  // Preconstruct a Flux functor
  EdgeFluxCalculator f;

  // Begin the time stepping
  for (double t = t_start; t < t_end; t += dt) {
    // Step forward in time with forward Euler
    hyperbolic_step(mesh, f, t, dt);
    // Update node values with triangle-averaged values
    post_process(mesh);
    
    // Update the viewer with new node positions
#if 1
    // Update viewer with nodes' new positions
    viewer.add_nodes(mesh.node_begin(), mesh.node_end(),
                     CoolColor(), NodePosition(), node_map);
#endif
//.........这里部分代码省略.........

int main(int argc, char** argv) {
  // Check arguments
  if (argc < 5) {
    std::cerr << "Usage: " << argv[0] << " NODES_FILE TETS_FILE\n";
    exit(1);
  }

  // Construct the first mesh
  MeshType mesh;
  //construct the second mesh
  MeshType mesh2;

  std::vector<typename MeshType::node_type> mesh_node;
  // Read all Points and add them to the Mesh
  std::ifstream nodes_file(argv[1]);
  Point p;
  while (CS207::getline_parsed(nodes_file, p)) {
    mesh_node.push_back(mesh.add_node(p));
  }

  // Read all mesh triangles and add them to the Mesh
  std::ifstream tris_file(argv[2]);
  std::array<int,3> t;
  while (CS207::getline_parsed(tris_file, t)) {
    mesh.add_triangle(mesh_node[t[0]], mesh_node[t[1]], mesh_node[t[2]]);
  }

  std::vector<typename MeshType::node_type> mesh_node2;
  // Read all Points and add them to the Mesh
  std::ifstream nodes_file2(argv[3]);
  while (CS207::getline_parsed(nodes_file2, p)) {
    mesh_node2.push_back(mesh2.add_node(p));
  }

  // Read all mesh triangles and add them to the Mesh
  std::ifstream tris_file2(argv[4]);
  while (CS207::getline_parsed(tris_file2, t)) {
    mesh2.add_triangle(mesh_node2[t[0]], mesh_node2[t[1]], mesh_node2[t[2]]);
  }
  
  //move the second mesh to the specified position
  for(auto it = mesh2.node_begin();it!=mesh2.node_end();++it){
    (*it).position().elem[1] +=4 ;
    (*it).position().elem[2] +=4 ;
  }

  // Print out the stats
  std::cout << mesh.num_nodes() << " "
            << mesh.num_edges() << " "
            << mesh.num_triangles() << std::endl;
  std::cout << mesh2.num_nodes() << " "
            << mesh2.num_edges() << " "
            << mesh2.num_triangles() << std::endl;



  //set the mass and velocity of each Node for the first mesh
  for (auto it = mesh.node_begin(); it != mesh.node_end(); ++it){
    (*it).value().mass = float(1)/mesh.num_nodes();
    (*it).value().velocity = Point(0, 10, 10);
  }
  //set the mass and velocity of each Node for the second mesh
  for (auto it = mesh2.node_begin(); it != mesh2.node_end(); ++it){
    (*it).value().mass = float(1)/mesh.num_nodes();
    (*it).value().velocity = Point(0, -10, -10);
  }

  //set K and L for each edge
  for (auto it = mesh.node_begin(); it != mesh.node_end(); ++it)
  {
    for (auto j = (*it).edge_begin(); j != (*it).edge_end(); ++j){
       (*j).value().L = (*j).length();
       (*j).value().K = 16000;
    }
  }

  for (auto it = mesh2.node_begin(); it != mesh2.node_end(); ++it)
  {
    for (auto j = (*it).edge_begin(); j != (*it).edge_end(); ++j){
       (*j).value().L = (*j).length();
       (*j).value().K = 16000;
    }
  }


  // Launch the SDLViewer
  CS207::SDLViewer viewer;
  auto node_map = viewer.empty_node_map(mesh);
  viewer.launch();

  viewer.add_nodes(mesh.node_begin(), mesh.node_end(), node_map);
  viewer.add_edges(mesh.edge_begin(), mesh.edge_end(), node_map);
  viewer.add_nodes(mesh2.node_begin(), mesh2.node_end(), node_map);
  viewer.add_edges(mesh2.edge_begin(), mesh2.edge_end(), node_map);

  viewer.center_view();
  
  //Begin the mass-spring simulation
  double dt = 0.0002;
  double t_start = 0.0;
//.........这里部分代码省略.........

int main(int argc, char* argv[]) {
  // Check arguments
  if (argc < 2) {
    std::cerr << "Usage: " << argv[0] << " NODES_FILE TETS_FILE\n";
    exit(1);
  }

  MeshType mesh;
  std::vector<typename MeshType::node_type> nodes;

  // Create a nodes_file from the first input argument
  std::ifstream nodes_file(argv[1]);
  
  // Interpret each line of the nodes_file as a 3D Point and add to the Mesh
  Point p;
  while (CS207::getline_parsed(nodes_file, p))
    nodes.push_back(mesh.add_node(p));

  // Create a trianges_file from the second input argument
  std::ifstream triangles_file(argv[2]);

  // Interpret each line of the tets_file as three ints which refer to nodes
  std::array<int,3> t;


  // add in the triangles
  while (CS207::getline_parsed(triangles_file, t))
    for (unsigned i = 1; i < t.size(); ++i)
      for (unsigned j = 0; j < i; ++j)
        for (unsigned k = 0; k < j; ++k)
        {
          mesh.add_triangle(nodes[t[i]], nodes[t[j]], nodes[t[k]]);
      }   

  // Set masses of nodes equal to 1/N where N is the number of nodes
  // and the initial velocities to 0. Also, get the indexes of
  // the nodes at positions (1,0,0) and (0,0,0)
  for(auto it=mesh.node_begin(); it != mesh.node_end(); ++ it)
  {
      (*it).value().mass = total_mass/mesh.num_nodes();
      (*it).value().velocity = Point(0.0,0.0,0.0);
  }

  // Set spring constants for each node equal to spring_const variable
  // and set initial length values equal to lengths of edges prior to
  // running the symplectic Euler steps
  for(auto it=mesh.edge_begin(); it != mesh.edge_end(); ++ it)
  {
      (*it).value().spring_constant = spring_const;
      (*it).value().initial_length = (*it).length();
  }


  // Set the triangle direction values so that we can determine which 
  // way to point normal vectors. This part assumes a convex shape
  Point center = get_center(mesh);
  for(auto it=mesh.triangle_begin(); it != mesh.triangle_end(); ++ it)
  {  
  	//std::cout << (*it).index() << std::endl;
  	set_normal_direction((*it),center);
  }

  // Print out the stats
  std::cout << mesh.num_nodes() << " " << mesh.num_edges() << std::endl;

  std::cout << "Center: " << get_center(mesh) << std::endl;
  // Launch the SDLViewer
  CS207::SDLViewer viewer;
  auto node_map = viewer.empty_node_map(mesh);
  viewer.launch();

  viewer.add_nodes(mesh.node_begin(), mesh.node_end(), node_map);
  viewer.add_edges(mesh.edge_begin(), mesh.edge_end(), node_map);

  viewer.center_view();

  // Begin the mass-spring simulation
  double dt = 0.0001;
  double t_start = 0;
  double t_end = 20.0;
  

  // Initialize constraints
  PlaneConstraint c1;
  //SelfCollisionConstraint c2;
  //auto combined_constraints = make_combined_constraints(c1,c2);
  


  for (double t = t_start; t < t_end; t += dt) {
    MassSpringForce ms_force;
    PressureForce p_force = PressureForce(0.0);
    DampingForce d_force = DampingForce(mesh.num_nodes());
    GravityForce g_force;
    WindForce w_force;

    (void) d_force; // prevents compiler from throwing error for unused variable
      
    
    if (t >= t_addgas - dt) {
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Check arguments
  if (argc < 3) {
    std::cerr << "Usage: shallow_water NODES_FILE TRIS_FILE\n";
    exit(1);
  }

  MeshType mesh;
  // HW4B: Need node_type before this can be used!
#if 1
  std::vector<typename MeshType::node_type> mesh_node;
#endif

  // Read all Points and add them to the Mesh
  std::ifstream nodes_file(argv[1]);
  Point p;
  while (CS207::getline_parsed(nodes_file, p)) {
    // HW4B: Need to implement add_node before this can be used!
#if 1
    mesh_node.push_back(mesh.add_node(p));
#endif
  }

  // Read all mesh triangles and add them to the Mesh
  std::ifstream tris_file(argv[2]);
  std::array<int,3> t;
  while (CS207::getline_parsed(tris_file, t)) {
    // HW4B: Need to implement add_triangle before this can be used!
#if 1
    mesh.add_triangle(mesh_node[t[0]], mesh_node[t[1]], mesh_node[t[2]]);
#endif
  }

  // Print out the stats
  std::cout << mesh.num_nodes() << " "
            << mesh.num_edges() << " "
            << mesh.num_triangles() << std::endl;

  // HW4B Initialization
  // Set the initial conditions according the type of input pattern
  if(argv[2][5]=='d'){
    Dam<MeshType> init;
    for(auto it= mesh.node_begin(); it != mesh.node_end(); ++it)
      init(*it);
    }
  else if((argv[2][5]=='p')){
    Pebble<MeshType> init;
    for(auto it= mesh.node_begin(); it != mesh.node_end(); ++it)
      init(*it);
    }
  else{
    Wave<MeshType> init;
    for(auto it= mesh.node_begin(); it != mesh.node_end(); ++it)
      init(*it);
    }

  // Set triangle values
  for (auto it=mesh.triangle_begin(); it!=mesh.triangle_end(); ++it) {
    (*it).value().Q = QVar(0.0,0.0,0.0);
    (*it).value().Q += (*it).node(0).value().Q;
    (*it).value().Q += (*it).node(1).value().Q;
    (*it).value().Q += (*it).node(2).value().Q;
    (*it).value().Q /= 3.0;
  }
 

  // Launch the SDLViewer
  CS207::SDLViewer viewer;
  viewer.launch();


  // HW4B: Need to define Mesh::node_type and node/edge iterator
  // before these can be used!
#if 1
  auto node_map = viewer.empty_node_map(mesh);
  viewer.add_nodes(mesh.node_begin(), mesh.node_end(),
                   CS207::DefaultColor(), NodePosition(), node_map);
  viewer.add_edges(mesh.edge_begin(), mesh.edge_end(), node_map);
#endif
  viewer.center_view();


  // HW4B: Timestep
  // CFL stability condition requires dt <= dx / max|velocity|
  // For the shallow water equations with u = v = 0 initial conditions
  //   we can compute the minimum edge length and maximum original water height
  //   to set the time-step
  // Compute the minimum edge length and maximum water height for computing dt
  double min_edge_length =( *mesh.edge_begin()).length();
  for (auto it=mesh.edge_begin(); it!=mesh.edge_end(); ++it) {
    if ((*it).length() < min_edge_length) {
      min_edge_length = (*it).length();
    }
  }
  double max_height = 0.0;
  for (auto it=mesh.node_begin(); it!=mesh.node_end(); ++it) {
    if ((*it).value().Q.h > max_height) {
      max_height = (*it).value().Q.h;
    }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  } 

  // Set the initial values of the triangles to the average of their nodes and finds S
  // Set the triangle direction values so that we can determine which 
  // way to point normal vectors. This part assumes a convex shape
  Point center = get_center(mesh); 
  for (auto it = mesh.triangle_begin(); it != mesh.triangle_end(); ++it) {
    auto t = *it; 
    if (t.index() < water_tris){
	    t.value().q_bar = (t.node1().value().q + 
	                       t.node2().value().q + 
	                       t.node3().value().q) / 3.0;
	    t.value().q_bar2 = t.value().q_bar;

	    double b_avg = (t.node1().value().b + 
	                    t.node2().value().b + 
	                    t.node3().value().b) / 3.0;
	    // finds the max dx and dy to calculate Source
	    dx = std::max(dx, fabs(t.node1().position().x - t.node2().position().x));
	    dx = std::max(dx, fabs(t.node2().position().x - t.node3().position().x));
	    dx = std::max(dx, fabs(t.node3().position().x - t.node1().position().x));
	    dy = std::max(dy, fabs(t.node1().position().y - t.node2().position().y));
	    dy = std::max(dy, fabs(t.node2().position().y - t.node3().position().y));
	    dy = std::max(dy, fabs(t.node3().position().y - t.node1().position().y));
	    t.value().S = QVar(0, -grav * t.value().q_bar.h * b_avg / dx, -grav * t.value().q_bar.h * b_avg / dy);
	  } 
    else
      set_normal_direction((*it),center);
  }

  // Calculate the minimum edge length and set edge inital condititons
  double min_edge_length = std::numeric_limits<double>::max();
  uint count = 0;
  for (auto it = mesh.edge_begin(); it != mesh.edge_end(); ++it, count++){
    if (count < water_edges)
      min_edge_length = std::min(min_edge_length, (*it).length());
    else {
      (*it).value().spring_constant = spring_const;
      (*it).value().initial_length = (*it).length();
    }
  }
	
  // Launch the SDLViewer
  CS207::SDLViewer viewer;
  viewer.launch();

  auto node_map = viewer.empty_node_map(mesh);
  viewer.add_nodes(mesh.node_begin(), mesh.node_end(),
                   Color(water_nodes), NodePosition(), node_map);
  viewer.add_edges(mesh.edge_begin(), mesh.edge_end(), node_map);
  // adds solid color-slows down program significantly
  //viewer.add_triangles(mesh.triangle_begin(), mesh.triangle_end(), node_map);
  viewer.center_view();


  // CFL stability condition requires dt <= dx / max|velocity|
  // For the shallow water equations with u = v = 0 initial conditions
  //   we can compute the minimum edge length and maximum original water height
  //   to set the time-step
  // Compute the minimum edge length and maximum water height for computing dt
  double dt = 0.25 * min_edge_length / (sqrt(grav * max_h));
  double t_start = 0;
  double t_end = 10;
  Point ball_loc = Point(0,0,0);
  double dh = 0;
  // double pressure = gas_const/(4/3*M_PI*radius*radius*radius);

//.........这里部分代码省略.........
  }*/


  // HW4B Initialization
  // Set the initial conditions
  int wave = 0, peddle=0, dam=1;
  if (wave){
	for ( auto it = mesh.node_begin(); it!= mesh.node_end(); ++it){
		auto x = (*it).position().x;
		auto y = (*it).position().y;
		double h = 1-0.75 * exp(-80 * ( (x-0.75)*(x-0.75) + y*y ));
		mesh.value((*it),QVar( h,0,0));
	}
  }
  else if (peddle){
	for ( auto it = mesh.node_begin(); it!= mesh.node_end(); ++it){
		auto x = (*it).position().x;
		auto y = (*it).position().y;
		double h = (x-0.75)*(x-0.75) + y*y -0.15*0.15 ;
		int H =0;
		if (h < 0)
			H = 1;
		mesh.value((*it), QVar(1+0.75*H,0,0));
	}
  }
  else if (dam){
	for ( auto it = mesh.node_begin(); it!= mesh.node_end(); ++it){
		auto x = (*it).position().x;
		int H =0;
		if (x < 0)
			H = 1;
		mesh.value((*it), QVar(1+0.75*H,0,0));
	}
  }
  
  
  // Perform any needed precomputation
// initialize triangle
  for (auto it = mesh.tri_begin(); it != mesh.tri_end(); ++it ) {
	(*it).value() = ((*it).node1().value() + (*it).node2().value() + (*it).node3().value())/3.0;
  }


  // Launch the SDLViewer
  CS207::SDLViewer viewer;
  viewer.launch();

  // HW4B: Need to define Mesh::node_type and node/edge iterator
  // before these can be used!
  
  auto node_map = viewer.empty_node_map(mesh);
  viewer.add_nodes(mesh.node_begin(), mesh.node_end(),
                   CS207::DefaultColor(), NodePosition(), node_map);
  viewer.add_edges(mesh.edge_begin(), mesh.edge_end(), node_map);

  viewer.center_view();
  
  // HW4B: Timestep
  // CFL stability condition requires dt <= dx / max|velocity|
  // For the shallow water equations with u = v = 0 initial conditions
  //   we can compute the minimum edge length and maximum original water height
  //   to set the time-step
  // Compute the minimum edge length and maximum water height for computing dt
  auto min_length = *std::min_element(mesh.edge_begin(), mesh.edge_end(), EdgeComparator);
  
  auto max_h = *std::max_element(mesh.node_begin(), mesh.node_end(), HeightComparator);
  
  double dt = 0.25 * min_length.length() / (sqrt(grav * max_h.value().h));
  double t_start = 0;
  double t_end = 0.1;

  // Preconstruct a Flux functor
  EdgeFluxCalculator f;

  // Begin the time stepping
  for (double t = t_start; t < t_end; t += dt) {
    // Step forward in time with forward Euler
    hyperbolic_step(mesh, f, t, dt);

    // Update node values with triangle-averaged values
    post_process(mesh);

    // Update the viewer with new node positions
    // HW4B: Need to define node_iterators before these can be used!
    viewer.add_nodes(mesh.node_begin(), mesh.node_end(),
                     CS207::DefaultColor(), NodePosition(), node_map);
    viewer.set_label(t);
	
    // These lines slow down the animation for small meshes.
    // Feel free to remove them or tweak the constants.
    if (mesh.num_nodes() < 100)
      CS207::sleep(0.05);
  }
auto elapsed = std::chrono::high_resolution_clock::now() - start;
long long microseconds = std::chrono::duration_cast<std::chrono::microseconds>(elapsed).count();

cout << microseconds << endl;
  
return 0;
}