///|
/// Oriented face fluxes: x faces are nx*(ny+1), y faces (nx+1)*ny.
pub struct Flux {
  grid : Grid
  x : Array[Double]
  y : Array[Double]
}

///|
pub fn Operator::flux(self : Operator, u : Field) -> Flux raise SolverError {
  require_same(self.grid, u.grid)
  let g = self.grid
  let (hx, hy) = g.spacing()
  let k = self.conductivity.values
  let x = Array::make(g.nx * (g.ny + 1), 0.0)
  let y = Array::make((g.nx + 1) * g.ny, 0.0)
  for j in 0..<=g.ny {
    for i in 0.. Array[Double] {
  self.x.copy()
}

///|
pub fn Flux::y_faces(self : Flux) -> Array[Double] {
  self.y.copy()
}

///|
/// Interior divergence; boundary nodes are zero, not boundary flux conditions.
pub fn Flux::divergence(self : Flux) -> Field raise SolverError {
  let g = self.grid
  let (hx, hy) = g.spacing()
  let values = Array::make(g.size(), 0.0)
  for j in 1.. Double {
  let mut lo = self.conductivity.values[0]
  let mut hi = lo
  for k in self.conductivity.values {
    if k < lo {
      lo = k
    }
    if k > hi {
      hi = k
    }
  }
  hi / lo
}