C++ MESH::triangle_begin方法代码示例

double hyperbolic_step(MESH& m, FLUX& f, double t, double dt, Point ball_loc, int water_tris) {
  for (auto tri_it = m.triangle_begin(); (*tri_it).index() != water_tris; ++tri_it) {
    QVar sum = QVar(0,0,0);

    /* Defines the force of the floating objects by checking its submerged 
    *  cross section of the ball intersects with the centers of adjacent triangles 
    */
    floatObj obj = floatObj();
    if (norm((*tri_it).center()-Point(ball_loc.x, ball_loc.y, 0)) < ball_loc.z)
      obj = floatObj(total_mass*grav, M_PI*ball_loc.z*ball_loc.z, density);

    Point norm_vec = (*tri_it).normal((*tri_it).node1(), (*tri_it).node2());
    // Calcuate flux using adjacent triangle's Q
    if ((*tri_it).adj_triangle((*tri_it).node1(), (*tri_it).node2()).is_valid())
      sum += f(norm_vec.x, norm_vec.y, dt, (*tri_it).value().q_bar, 
      		(*tri_it).adj_triangle((*tri_it).node1(), (*tri_it).node2()).value().q_bar, obj);
    // Boundary conditions
    else
      sum += f(norm_vec.x, norm_vec.y, dt, (*tri_it).value().q_bar, 
      		QVar((*tri_it).value().q_bar.h, 0, 0), obj);

    norm_vec = (*tri_it).normal((*tri_it).node2(), (*tri_it).node3());
    // Calcuate flux using adjacent triangle's Q
    if ((*tri_it).adj_triangle((*tri_it).node2(), (*tri_it).node3()).is_valid()) 
      sum += f(norm_vec.x, norm_vec.y, dt, (*tri_it).value().q_bar, 
      		(*tri_it).adj_triangle((*tri_it).node2(), (*tri_it).node3()).value().q_bar, obj); 
    // Boundary conditions
    else 
      sum += f(norm_vec.x, norm_vec.y, dt, (*tri_it).value().q_bar, 
      		QVar((*tri_it).value().q_bar.h, 0, 0), obj);

    norm_vec = (*tri_it).normal((*tri_it).node3(), (*tri_it).node1());
    // Calcuate flux using adjacent triangle's Q
    if ((*tri_it).adj_triangle((*tri_it).node3(), (*tri_it).node1()).is_valid())
      sum += f(norm_vec.x, norm_vec.y, dt, (*tri_it).value().q_bar, 
      		(*tri_it).adj_triangle((*tri_it).node3(), (*tri_it).node1()).value().q_bar, obj); 
    // Boundary conditions
    else
      sum += f(norm_vec.x, norm_vec.y, dt, (*tri_it).value().q_bar, 
      		QVar((*tri_it).value().q_bar.h, 0, 0), obj);

    (*tri_it).value().q_bar2 = ((*tri_it).value().S * dt/(*tri_it).area()) + (*tri_it).value().q_bar - 
                               (sum * dt/(*tri_it).area());
	}

  /* In order to use original (as opposed to update) q_bar values in the equation above, we have to update
   * q_bar values for all triangles after the q_bar values have been calculated.
   * Updates the Source term
   */
  for (auto tri_it = m.triangle_begin(); (*tri_it).index() != water_tris; ++tri_it){
    (*tri_it).value().S = (*tri_it).value().S/(*tri_it).value().q_bar.h*(*tri_it).value().q_bar2.h;
    (*tri_it).value().q_bar = (*tri_it).value().q_bar2;
  }

  return t + dt;
}

double hyperbolic_step(MESH& m, FLUX& f, double t, double dt) {
  // HW4B: YOUR CODE HERE
  // Step the finite volume model in time by dt.

  // Pseudocode:
  // Compute all fluxes. (before updating any triangle Q_bars)
  // For each triangle, update Q_bar using the fluxes as in Equation 8.
  //  NOTE: Much like symp_euler_step, this may require TWO for-loops

  /* 1. Initialize a vector temp_Q of size @a m.num_triangles() for storing all temp Q_bar values
   * 2. Use a TriangleIterator to iterate through all the triangles, for each tri_it:
   *      a. Compute the sum of fluxes sum_fluxes
   *      b. temp_Q[(*it).index()] = sum_fluxes
   * 3. Use indices of triangles to do another for loop, for each triangle(i):
   *      a. triangle_i.value().Q_bar = temp_Q[i]
   */
  std::vector<QVar> temp_fluxes(m.num_triangles(), QVar());
  for (auto tri_it = m.triangle_begin(); tri_it != m.triangle_end(); ++tri_it) {
    auto qk = (*tri_it).value();
    
    QVar F_k(0, 0, 0);
    for (size_type i = 0; i < 3; ++i) {
      auto edge = (*tri_it).edge(i);
      QVar qm(0, 0, 0);
      if (m.adj_triangle1_index(edge) == -1 || m.adj_triangle2_index(edge) == -1) {
          //set_boundary_conditions
        qm = QVar(qk.h, 0, 0);
      } else {
        size_type adj_tri_idx = m.adj_triangle1_index(edge) == (*tri_it).index() ?
            m.adj_triangle2_index(edge) : m.adj_triangle1_index(edge);

        auto adj_triangle = m.triangle(adj_tri_idx);
        qm = adj_triangle.value(); 
      }
      auto norm_ke = m.out_norm((*tri_it), edge); 
      F_k = F_k + f(norm_ke.x, norm_ke.y, dt, qk, qm);
    } 
    // print_triangle(m, (*tri_it), f, t, dt);
    temp_fluxes[(*tri_it).index()] = qk - F_k * (dt / (*tri_it).area());
  }
  
  for (auto tri_it = m.triangle_begin(); tri_it != m.triangle_end(); ++tri_it) {
    (*tri_it).value() = temp_fluxes[(*tri_it).index()];
  }
  
  return t + dt;
}

double get_volume(MESH& m) {
  double volume = 0.0;
  for(auto it = m.triangle_begin(); it != m.triangle_end(); ++it){
    auto tri = (*it);
    Point normal = get_normal_surface(tri);
    double area = tri.area();
    volume += normal.z * area * (tri.node(0).position().z + tri.node(1).position().z + tri.node(2).position().z) / 3.0;
  }
  return volume;
}

double hyperbolic_step(MESH& m, FLUX& f, double t, double dt) {
  // HW4B: YOUR CODE HERE
  // Step the finite volume model in time by dt.

  // Pseudocode:
  // Compute all fluxes. (before updating any triangle Q_bars)
  // For each triangle, update Q_bar using the fluxes as in Equation 8.
  //  NOTE: Much like symp_euler_step, this may require TWO for-loops
  
  // Implement Equation 7 from your pseudocode here.

  for(auto i = m.edge_begin(); i != m.edge_end(); ++i){
    if ((*i).triangle1().index() != (unsigned) -1 && (*i).triangle2().index() != (unsigned) -1 ){
      MeshType::Triangle trik = (*i).triangle1();
      MeshType::Triangle trim = (*i).triangle2();
      unsigned int edge_k = 0;
      unsigned int edge_m = 0;
      //which edge (*i) is in trik and trim
      while(trik.node(edge_k).index()== (*i).node1().index() 
        || trik.node(edge_k).index()== (*i).node2().index() )
        ++edge_k;
      while(trim.node(edge_m).index()== (*i).node1().index() 
        || trim.node(edge_m).index()== (*i).node2().index() )
        ++edge_m;
      QVar flux = f(trik.normal(edge_k).x, trik.normal(edge_k).y, dt, trik.value().Q, trim.value().Q);
      trik.value().F[edge_k] = flux;
      trim.value().F[edge_m] = -flux;
    }
    else{
      MeshType::Triangle trik;
      if ((*i).triangle1().index() != (unsigned) -1)
        trik = (*i).triangle1();
      else
        trik = (*i).triangle2();
      unsigned int edge_k = 0;
      while(trik.node(edge_k).index()== (*i).node1().index() 
        || trik.node(edge_k).index()== (*i).node2().index() )
        ++edge_k;
      QVar flux = f(trik.normal(edge_k).x, trik.normal(edge_k).y, dt, trik.value().Q, QVar(trik.value().Q.h, 0, 0));
      trik.value().F[edge_k] = flux;
    }
  }

  for(auto i = m.triangle_begin(); i != m.triangle_end(); ++i){
    QVar sum = QVar(0, 0, 0);
    for (int j = 0; j < 3; ++j){
      sum += (*i).value().F[j];
    }
    (*i).value().Q = (*i).value().Q-dt/(*i).area()*sum;
  }
  
  return t + dt;
};

void post_process(MESH& m) {
  // HW4B: Post-processing step
  // Translate the triangle-averaged values to node-averaged values
  // Implement Equation 9 from your pseudocode here
  for (auto it = m.node_begin(); it != m.node_end(); ++it){
    double sumarea=0;
    QVar sumQ = QVar(0, 0, 0);
    for(auto j = m.triangle_begin(*it); j != m.triangle_end(*it); ++j){
      sumarea += (*j).area();
      sumQ += (*j).value().Q * (*j).area();
    }
    (*it).value().Q = sumQ/sumarea;
  }
}