C++ FluidSim类代码示例

TEST(PressureTest, LinearEquationSetup_Complex)
{
  glm::ivec2 size(50);

  FluidSim sim;
  sim.initialize(1.0f, size.x, size.y);
  sim.set_boundary(boundary_phi);

  AddParticles(size, sim, complex_boundary_phi);

  sim.add_force(0.01f);

  Velocity velocity(*device, size);
  Texture solidPhi(*device, size.x, size.y, vk::Format::eR32Sfloat);
  Texture liquidPhi(*device, size.x, size.y, vk::Format::eR32Sfloat);

  BuildInputs(*device, size, sim, velocity, solidPhi, liquidPhi);

  LinearSolver::Data data(*device, size, VMA_MEMORY_USAGE_CPU_ONLY);
  Buffer<glm::ivec2> valid(*device, size.x * size.y, VMA_MEMORY_USAGE_CPU_ONLY);

  Pressure pressure(*device, 0.01f, size, data, velocity, solidPhi, liquidPhi, valid);

  pressure.BuildLinearEquation();
  device->Handle().waitIdle();

  CheckDiagonal(size, data.Diagonal, sim, 1e-3f);  // FIXME can we reduce error tolerance?
  CheckWeights(size, data.Lower, sim, 1e-3f);      // FIXME can we reduce error tolerance?
  CheckDiv(size, data.B, sim);
}

TEST(LinearSolverTests, IncompletePoisson_Simple_PCG)
{
  glm::ivec2 size(50);

  FluidSim sim;
  sim.initialize(1.0f, size.x, size.y);
  sim.set_boundary(boundary_phi);

  AddParticles(size, sim, boundary_phi);

  sim.add_force(0.01f);
  sim.compute_phi();
  sim.extrapolate_phi();
  sim.apply_projection(0.01f);

  LinearSolver::Data data(*device, size, VMA_MEMORY_USAGE_CPU_ONLY);

  BuildLinearEquation(size, data.Diagonal, data.Lower, data.B, sim);

  IncompletePoisson preconditioner(*device, size);

  LinearSolver::Parameters params(LinearSolver::Parameters::SolverType::Iterative, 1000, 1e-5f);
  ConjugateGradient solver(*device, size, preconditioner);

  solver.Bind(data.Diagonal, data.Lower, data.B, data.X);
  solver.Solve(params);

  device->Queue().waitIdle();

  CheckPressure(size, sim.pressure, data.X, 1e-5f);

  std::cout << "Solved with number of iterations: " << params.OutIterations << std::endl;
}

void CheckWeights(const glm::ivec2& size,
                  Buffer<glm::vec2>& buffer,
                  FluidSim& sim,
                  float error = 1e-6)
{
  std::vector<glm::vec2> pixels(size.x * size.y);
  CopyTo(buffer, pixels);

  for (int i = 1; i < size.x - 1; i++)
  {
    for (int j = 1; j < size.y - 1; j++)
    {
      int index = i + size.x * j;
      EXPECT_NEAR(sim.matrix(index - 1, index), pixels[index].x, error);
      EXPECT_NEAR(sim.matrix(index, index - size.x), pixels[index].y, error);
    }
  }
}

void PrintDiagonal(const glm::ivec2& size, FluidSim& sim)
{
  for (int j = 0; j < size.y; j++)
  {
    for (int i = 0; i < size.x; i++)
    {
      int index = i + size.x * j;
      std::cout << "(" << sim.matrix(index, index) << ")";
    }
    std::cout << std::endl;
  }
}

void CheckDiagonal(const glm::ivec2& size, Buffer<float>& buffer, FluidSim& sim, float error = 1e-6)
{
  std::vector<float> pixels(size.x * size.y);
  CopyTo(buffer, pixels);

  for (int i = 0; i < size.x; i++)
  {
    for (int j = 0; j < size.y; j++)
    {
      int index = i + size.x * j;
      EXPECT_NEAR(sim.matrix(index, index), pixels[index], error);
    }
  }
}

int main(){
	FluidSim sim;
	sim.runSim();
}

int main(int argc, char** argv)
{
    bool run = GL_TRUE;

    if(!glfwInit())
    {
        exit(EXIT_FAILURE);
    }

    if(!glfwOpenWindow(windows_width, windows_height, 8, 8, 8, 8, 24, 0, GLFW_WINDOW))
    {
        glfwTerminate();
        exit(EXIT_FAILURE);
    }

    glewInit();
    if (!glewIsSupported( "GL_VERSION_2_0 "
                          "GL_ARB_pixel_buffer_object"
                        )) {
        fprintf( stderr, "ERROR: Support for necessary OpenGL extensions missing.");
        fflush( stderr);
        return false;
    }

    //glfwSetKeyCallback(Viewer::keyCallback);
    //glfwSetMouseButtonCallback(Viewer::mouseButtonCallback);
    //glfwSetMousePosCallback(Viewer::mousePosCallback);
    //glfwSetMouseWheelCallback(Viewer::mouseWheelCallback);

    glClearColor(0.0f, 0.0f, 0.0f, 1.0f);
    glEnable(GL_DEPTH_TEST);
    glDepthFunc(GL_LEQUAL);

    FluidSim sph;
    sph.initialize();
    sph.resetFrameNum();
    sph.outputObj(true);

    vec3 center(0.0f, -0.15f, 0.0f);
    Viewer::camera()->setFocal(center.x, center.y, center.z);
    Viewer::camera()->resetLookAt();

    old_now = glfwGetTime();
    while(run)
    {
        unsigned int frame_num = sph.getFrameNum();
        if(frame_num > 450)
            exit(0);
        Viewer::camera()->aim();

        printf("frame number: %d\n", frame_num);
        sph.compute(0.004);
        sph.draw();

        utils::drawAxis();
        utils::grabScreen(windows_width, windows_height, sph.getFrameNum());
        now = glfwGetTime();
        char fpsInfo[256];
        sprintf(fpsInfo, "Frame: %u. Time cost per frame: %f", frame_num, now - old_now);
        old_now = now;
        glfwSetWindowTitle(fpsInfo);

        glfwSwapBuffers();

        run = !glfwGetKey(GLFW_KEY_ESC) && glfwGetWindowParam(GLFW_OPENED);
    }

    glfwTerminate();
    exit(EXIT_SUCCESS);
}

//Main testing code
//-------------
int main(int argc, char **argv)
{
   
   //Setup viewer stuff
   Gluvi::init("GFM Free Surface Liquid Solver with Static Variational Boundaries", &argc, argv);
   Gluvi::camera=&cam;
   Gluvi::userDisplayFunc=display;
   Gluvi::userMouseFunc=mouse;
   Gluvi::userDragFunc=drag;
   glClearColor(1,1,1,1);
   
   glutTimerFunc(1000, timer, 0);
   
   //Set up the simulation
   sim.initialize(grid_width, grid_resolution, grid_resolution);
   
   //set up a circle boundary
   sim.set_boundary(boundary_phi);
   
   //Stick some liquid particles in the domain
   for(int i = 0; i < sqr(grid_resolution); ++i) {
      float x = randhashf(i*2, 0,1);
      float y = randhashf(i*2+1, 0,1);
      Vec2f pt(x,y);
      if(boundary_phi(pt) > 0 && pt[0] > 0.5)
         sim.add_particle(pt);
   }

   Gluvi::run();

   return 0;
}

void display(void)
{
  
   if(draw_grid) {
      glColor3f(0,0,0);
      glLineWidth(1);
      draw_grid2d(Vec2f(0,0), sim.dx, sim.ni, sim.nj);  
   }

   if(draw_boundaries) {
      glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
      draw_circle2d(c0, rad0, 50); 
      draw_circle2d(c1, rad1, 50); 
      draw_circle2d(c2, rad2, 50); 
      draw_circle2d(c3, rad3, 50); 
      
      //There's a bug, so draw one more(?)
      draw_circle2d(c3, 0, 10);
   }

   if(draw_particles) {
      glColor3f(0,0,0);
      glPointSize(3);
      draw_points2d(sim.particles);
   }

   if(draw_velocities) {
      for(int j = 0;j < sim.nj; ++j) for(int i = 0; i < sim.ni; ++i) {
         Vec2f pos((i+0.5f)*sim.dx,(j+0.5f)*sim.dx);
         draw_arrow2d(pos, pos + 0.01f*sim.get_velocity(pos), 0.1f*sim.dx);
      }
   }

}

void timer(int junk)
{
   sim.advance(timestep);

   glutPostRedisplay();
   glutTimerFunc(30, timer, 0);

}

void display(void)
{
  
   if(draw_grid) {
      glColor3f(0,0,0);
      glLineWidth(1);
      draw_grid2d(Vec2f(0,0), sim.dx, sim.ni, sim.nj);  
   }

   if(draw_boundaries) {
      glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
      draw_circle2d(c0, rad0, 50); 
      
      //There's a bug, so draw one more(?)
      draw_circle2d(c3, 0, 10);
   }

   if(draw_particles) {
      glColor3f(0,0,1);
      glPointSize(3);
      glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
      draw_points2d(sim.particles);
      for(unsigned int p = 0; p < sim.particles.size(); ++p) {
         draw_circle2d(sim.particles[p], sim.particle_radius, 20);
      }
   }

   if(draw_velocities) {
      glColor3f(1,0,0);
      for(int j = 0;j < sim.nj; ++j) for(int i = 0; i < sim.ni; ++i) {
         Vec2f pos((i+0.5)*sim.dx,(j+0.5)*sim.dx);
         draw_arrow2d(pos, pos + 0.01f*sim.get_velocity(pos), 0.01*sim.dx);
      }
   }


}