///|
/// Axis used by line-integral and flow-rate helpers.
pub(all) enum Axis {
  Horizontal
  Vertical
} derive(Debug, Eq)

///|
/// Build a velocity field from the simulation state.
pub fn Simulation::velocity_field(self : Simulation) -> VectorField2D {
  let field = VectorField2D::new(size=self.size)
  for y in 0.. Field2D {
  let field = Field2D::new(size=self.size)
  for y in 0.. Double {
  let mut energy = 0.0
  for y in 0.. Point {
  let mut x = 0.0
  let mut y = 0.0
  for row in 0.. Double {
  let mut total = 0.0
  match axis {
    Horizontal => {
      let y = index.clamp(min=0, max=self.size.height - 1)
      for x in 0.. {
      let x = index.clamp(min=0, max=self.size.width - 1)
      for y in 0.. Double {
  if self.is_solid(x~, y~) {
    0.0
  } else {
    let left = self.cell(x=(x - 1).clamp(min=0, max=self.size.width - 1), y~)
    let right = self.cell(x=(x + 1).clamp(min=0, max=self.size.width - 1), y~)
    let bottom = self.cell(x~, y=(y - 1).clamp(min=0, max=self.size.height - 1))
    let top = self.cell(x~, y=(y + 1).clamp(min=0, max=self.size.height - 1))
    (right.uy - left.uy) * 0.5 - (top.ux - bottom.ux) * 0.5
  }
}

///|
/// Return a row-major vorticity field.
pub fn Simulation::vorticity_field(self : Simulation) -> Field2D {
  let field = Field2D::new(size=self.size)
  for y in 0.. Double {
  let delta = self.mass() - reference
  if delta < 0.0 {
    -delta
  } else {
    delta
  }
}

///|
/// Estimate the net force on solid cells from local interface momentum exchange.
pub fn Simulation::force_on_mask(
  self : Simulation,
  mask~ : DomainMask,
) -> Point {
  let mut force_x = 0.0
  let mut force_y = 0.0
  for y in 0..