gistfile1.cpp

void StableFluidsSimWithMarkers::project(
    int N, ArrayXs * u, ArrayXs * v, ArrayXs * u0, ArrayXs * v0) {
  // Set velocities on solid boundaries to zero
  zeroSolidBoundaries(u, v);

  // Compute discrete divergence on velocity field
  ArrayXs div(N+1, N+1);
  div.setZero();

  for (int i=1; i<N; i++) {
    for (int j=1; j<N; j++) {
      div(i, j) = (A(u, i, j) - A(u, i, j-1) + A(v, i, j) - A(v, i-1, j));
    }
  }

  OUT(cout << div << endl);
  // Adjust divergence to account for V change
  scalar dgrad = m_nu * (fluid_cell_count() - m_original_fluid_cell_count);
  OUT(cout << "dgrad = " << dgrad << endl);
  for (int i=1; i<N; i++) {
    for (int j=1; j<N; j++) {
      div(i, j) += dgrad;
    }
  }

  // Solve for pressure
  ArrayXs pressure(N+2, N+2);
  pressure.setZero();

  for (int k=0; k<100; k++) {
    for (int i=1; i<N+1; i++) {
      for (int j=1; j<N+1; j++) {
        int adj_fluid_cells = 0;
        scalar p = pressure(i, j);
        pressure(i, j) = -div(i, j) / N;

        if (!m_has_solid(i-1, j)) {
          if (m_has_fluid(i-1, j))
            pressure(i, j) += pressure(i-1, j);
          adj_fluid_cells++;
        }

        if (!m_has_solid(i+1, j)) {
          if (m_has_fluid(i+1, j))
            pressure(i, j) += pressure(i+1, j);
          adj_fluid_cells++;
        }

        if (!m_has_solid(i, j-1)) {
          if (m_has_fluid(i, j-1))
            pressure(i, j) += pressure(i, j-1);
          adj_fluid_cells++;
        }

        if (!m_has_solid(i, j+1)) {
          if (m_has_fluid(i, j+1))
            pressure(i, j) += pressure(i, j+1);
          adj_fluid_cells++;
        }

        pressure(i, j) /= adj_fluid_cells;
      }
    }
  }

  // Apply pressure gradient to velocity field
  for (int i=0; i<N+1; i++) {
    for (int j=0; j<N+1; j++) {
      A(u, i, j) -= N * (pressure(i, j+1) - pressure(i, j));
      A(v, i, j) -= N * (pressure(i+1, j) - pressure(i, j));
    }
  }

  // Set velocities on solid boundaries to zero
  zeroSolidBoundaries(u, v);

  // Propagate liquid velocities into 5 layers of air
  for (int k=0; k<3; k++) {
    for (int i=1; i<=N; i++) {
      for (int j=1; j<=N; j++) {
        if (m_has_solid(i, j) || m_has_fluid(i, j)) continue;

        if (i > 0 && m_has_fluid(i-1, j)) {
          A(v, i, j) = A(v, i-1, j);
        }

        if (i < N+1 && m_has_fluid(i+1, j)) {
          A(v, i-1, j) = A(v, i, j);
        }

        if (j > 0 && m_has_fluid(i, j-1)) {
          A(u, i, j) = A(u, i, j-1);
        }

        if (j < N+1 && m_has_fluid(i, j+1)) {
          A(u, i-1, j) = A(u, i, j);
        }
      }
    }
  }
}