code: def project(vx, vy, occlusion):
      
      p = np.zeros(vx.shape)
      div = -0.5 * (np.roll(vx, -1, axis=1) - np.roll(vx, 1, axis=1)
                  + np.roll(vy, -1, axis=0) - np.roll(vy, 1, axis=0))
      div = make_continuous(div, occlusion)
  
      for k in range(50):
          p = (div + np.roll(p, 1, axis=1) + np.roll(p, -1, axis=1)
                   + np.roll(p, 1, axis=0) + np.roll(p, -1, axis=0))/4.0
          p = make_continuous(p, occlusion)
  
      vx = vx - 0.5*(np.roll(p, -1, axis=1) - np.roll(p, 1, axis=1))
      vy = vy - 0.5*(np.roll(p, -1, axis=0) - np.roll(p, 1, axis=0))
  
      vx = occlude(vx, occlusion)
      vy = occlude(vy, occlusion)
      return vx, vy
docs: Project the velocity field to be approximately mass-conserving,
         using a few iterations of Gauss-Seidel.