
\section{\label{formLibrary}The ``form'' library}

\begin{verbatim}
module Lattice (Cell (cellEnum, cellArray, cellAccum, dimension),
                Lattice (MkLattice),
                Form,
                cells, boundary,
                form, accumForm, lattice, rank, (!#), (//#),
                listForm, fassocs, findices, felems,
                fmap, fmapi, fzipWith) where

-- The underlying technology of lattices and their forms...
-- This parallels the `PreludeArray' section for arrays.

infixl 9 !#
infixl 9 //#

-- The `Cell' class is the analogue of the `Ix' class in `PreludeCore'.

class (Eq a) => Cell a where
  cellEnum  :: Int -> [a] -> [a]
  cellArray :: Int -> [a] -> [Assoc a b] -> a -> b
  cellAccum :: Int -> [a] -> [Assoc a b] -> a -> [b]
  dimension :: a -> Int  


-- First field in `Lattice' defines extent of the complex, where required.
-- Second field is the `boundary' function for the complex.

data (Cell a) => Lattice a = MkLattice (Int -> [a]) (a -> [Assoc a Int])
data (Cell a) => Form a b  = MkForm (Lattice a) Int (a -> b)

cells     :: (Cell a) => Lattice a -> Int -> [a]
boundary  :: (Cell a, Num b) => Lattice a -> a -> [Assoc a b]


form      :: (Cell a) => Lattice a -> Int -> [Assoc a b] -> Form a b
accumForm :: (Cell a) => Lattice a -> Int -> [Assoc a b] -> Form a [b]
lattice   :: (Cell a) => Form a b -> Lattice a
rank      :: (Cell a) => Form a b -> Int
(!#)      :: (Cell a) => Form a b -> a -> b

(//#)     :: (Cell a) => Form a b -> [Assoc a b] -> Form a b

listForm  :: (Cell a) => (Lattice a) -> Int -> [b] -> Form a b

fassocs   :: (Cell a) => Form a b -> [Assoc a b]
findices  :: (Cell a) => Form a b -> [a]
felems    :: (Cell a) => Form a b -> [b]

fmap      :: (Cell a) => (b -> c) -> Form a b -> Form a c
fmapi     :: (Cell a) => (a -> b -> c) -> Form a b -> Form a c
fzipWith  :: (Cell a) => (b -> c -> d) -> Form a b -> Form a c
                                                   -> Form a d

-- `Lattice' primitives...

cells    (MkLattice cls _) r = cellEnum r (cls r)
boundary (MkLattice _ bdy) i = [j := fromIntegral s | j := s <- bdy i]
 
-- `Form' primitives...

form      l@(MkLattice cls _) r ivs = MkForm l r (cellArray r (cls r) ivs)
accumForm l@(MkLattice cls _) r ivs = MkForm l r (cellAccum r (cls r) ivs)
lattice (MkForm l _ _) = l
rank    (MkForm _ r _) = r
(!#)    (MkForm _ _ d) = d

-- Higher level operations on Forms...

a //# us   = form (lattice a) (rank a)
             ([i := a !# i | i <- findices a \\ [i | i := _ <- us]] ++ us)

listForm l r vs = form l r (zipWith (:=) (cells l r) vs)

fassocs  a = [i := a !# i | i <- cells (lattice a) (rank a)]
findices a = cells (lattice a) (rank a)
felems   a = [a !# i | i <- findices a]

fmap   f a = form l r [i := f (a !# i) | i <- cells l r]
             where l = lattice a
                   r = rank    a

fmapi  f a = form l r [i := f i (a !# i) | i <- cells l r]
             where l = lattice a
                   r = rank    a

fzipWith f a b = form l r [i := f (a !# i) (b !# i) | i <- cells l r]
                 where l = lattice a
                       r = rank    a
\end{verbatim}

\section{\label{cellInstance}An example instance of the ``Cell'' class}

\begin{verbatim}
module Cell2d (Cell2d (Cell2d0, Cell2d1, Cell2d2)) where

-- A `Cell' instance suitable for representing cells in a 2d grid.

import Lattice (Cell (cellEnum, cellArray, cellAccum))

data Cell2d = Cell2d0 (Int, Int)     |
              Cell2d1 (Int, Int) Int |
              Cell2d2 (Int, Int)         deriving (Text)

ix0 (Cell2d0 i)        = i
ix1 (Cell2d1 (x, y) d) = (x, y, d)
ix2 (Cell2d2 i)        = i

bd0 [c1, c2]    = (ix0 c1, ix0 c2)
bd1 [c1, c2]    = (ix1 c1, ix1 c2)
bd2 [c1, c2]    = (ix2 c1, ix2 c2)

cell0 i         = Cell2d0 i
cell1 (x, y, d) = Cell2d1 (x,y) d
cell2 i         = Cell2d2 i

instance Cell Cell2d where

  cellEnum r cls =
    case r of
      0 -> [cell0 i | i <- range (bd0 cls)]
      1 -> [cell1 i | i <- range (bd1 cls)]
      2 -> [cell2 i | i <- range (bd2 cls)]

  cellArray r cls ivs =
    case r of
      0 -> (array (bd0 cls) [ix0 i := v | i := v <- ivs] !) . ix0
      1 -> (array (bd1 cls) [ix1 i := v | i := v <- ivs] !) . ix1
      2 -> (array (bd2 cls) [ix2 i := v | i := v <- ivs] !) . ix2

  cellAccum r cls ivs =
    case r of
      0 -> (accumArray (flip (:)) [] (bd0 cls) [ix0 i := v | i := v <- ivs] !)
             . ix0
      1 -> (accumArray (flip (:)) [] (bd1 cls) [ix1 i := v | i := v <- ivs] !)
             . ix1
      2 -> (accumArray (flip (:)) [] (bd2 cls) [ix2 i := v | i := v <- ivs] !)
             . ix2
\end{verbatim}

\section{\label{topologyExample}Example definition of lattice geometry}

\begin{verbatim}
module Torus2d (torus2d) where

-- `torus2d' returns an object of type `Lattice Cell2d' implementing
-- a 2 dimensional grid with periodic boundary conditions.

import Lattice (Cell (cellEnum, cellArray, cellAccum), Lattice (MkLattice))

import Cell2d (Cell2d (Cell2d0, Cell2d1, Cell2d2))

torus2d :: Int -> Int -> Lattice Cell2d
torus2d lx ly = MkLattice (cls lx ly) (bdy lx ly)

-- Lattice bounds, specified by smallest and largest cell of each dimension:

cls lx ly 0 = [Cell2d0 (0, 0),   Cell2d0 (lx - 1, ly - 1)]
cls lx ly 1 = [Cell2d1 (0, 0) 0, Cell2d1 (lx - 1, ly - 1) 1]
cls lx ly 2 = [Cell2d2 (0, 0),   Cell2d2 (lx - 1, ly - 1)]

-- Lattice topology (cell boundary):

bdy _ _ (Cell2d0 _) = []

bdy lx ly (Cell2d1 (x, y) d) =
  [Cell2d0 (x, y) :=  -1,
   Cell2d0 (case d of
              0 -> ((x + 1) `mod` lx, y)
              1 -> (x, (y + 1) `mod` ly)) := 1]

bdy lx ly (Cell2d2 (x, y)) =
  [Cell2d1 (x, y) 0 :=  1,
   Cell2d1 ((x + 1) `mod` lx, y) 1  :=  1,
   Cell2d1 (x, (y + 1) `mod` ly) 0  := -1,
   Cell2d1 (x, y) 1 := -1]
\end{verbatim}


\section{\label{derivs}Derivatives}

\subsection*{Example 1: the gradient}

An example alluded to several times is
the gradient $\mbox{\boldmath $\nabla$} a$ of a 0-form $a$
(i.e., a function defined on the
sites of a lattice).  The gradient can be implemented by
\begin{verbatim}
  gradient :: (Cell a, Num b) => Form a b -> Form a b
  gradient a = form l 1 [i := cntrct a (boundary l i) | i <- cells l 1]
               where l = lattice a
\end{verbatim}
To see how this works, refer back the figure \ref{incidence}.  The
result on link $F$ is given by\footnote{The order of the list returned
by \verb$boundary$ hasn't been defined.  It doesn't effect the
result of any of the calculations defined here.  The explicit
examples choose some order arbitrarily.}
\begin{tabbing}
\verb$  $\= \kill
\>\verb$  cntrct a (boundary l $$F$\verb$)$ \\
\>\verb$= cntrct a [$$H$\verb$ := -1, $$G$\verb$ := 1]$ \\
\>\verb$= ((-1) * (a !# $$H$\verb$)) + (1 * (a !# $$G$\verb$))$ \\
\>\verb$= (a !# $$G$\verb$) - (a !# $$H$\verb$)$
\end{tabbing}
which is, as advertised, a plausible lattice approximation to the gradient
$\partial a / \partial x$ of \verb$a$ on this
$x$-directed link (if $x$ is the horizontal coordinate in the figure).

\subsection*{\label{divergence}Example 2: the divergence}

Another operator familiar from differential vector calculus is
the divergence, $\mbox{\boldmath $\nabla$} . \mbox{\boldmath $a$}$, of
a vector field $\mbox{\boldmath $a$}$.
Interpret the vector field as a 1-form.  First introduce the function
\begin{verbatim}
  bndry :: (Cell a, Num b) => Lattice a -> [Assoc a b] -> [Assoc a b]
  bndry l ivs = concat [mltply v (boundary l i) | i := v <- ivs]
\end{verbatim}
Given an association list, \verb$bndry$ returns another which is
a concatenation of the boundaries of all cells in the original
list, each multiplied by the value associated with the
cell\footnote{This probably seems like an odd concept.  It's inspired
by the boundary operation on chains in homology---an association
list is something like a chain, as the manipulations in
appendix \ref{sign} hint.}.  Divergence can be implemented by
\begin{verbatim}
  divergence :: (Cell a, Num b) => Form a b -> Form a b
  divergence a = (fmap sum . accumForm l 0 . bndry l . fassocs) a
                 where l = lattice a
\end{verbatim}
\begin{figure}[hbt]
\centerline{
\begin{picture}(220, 130)(-20, -30)
\multiput(-10,  0)( 0, 40){3}{\line(1, 0){100}}
\multiput( 0, -10)(40,  0){3}{\line(0, 1){100}}
%
\put(40, 40){\circle*{6}}
\put(30, 30){$A$}
\put(45, 40){\rule[-1pt]{30pt}{3pt}}
\put(57, 44){$B$}
\put(39, 45){\rule{3pt}{30pt}}
\put(44, 57){$C$}
\put( 5, 40){\rule[-1pt]{30pt}{3pt}}
\put(17, 44){$D$}
\put(39,  5){\rule{3pt}{30pt}}
\put(44, 17){$E$}
%
\put(100, 80){
\begin{minipage}[t]{2in}
\begin{eqnarray*}
I (B, A) & = & -1 \\
I (C, A) & = & -1 \\
I (D, A) & = & +1 \\
I (E, A) & = & +1
\end{eqnarray*}
\end{minipage}
}
\end{picture}
}
\caption{\label{divergencePicture}Calculation of the ``divergence''.}
\end{figure}
Referring to figure \ref{divergencePicture}, the only links on whose boundary
the site $A$ lies are $B$, $C$, $D$ and $E$.  So, restricting attention
to associations of $A$ in the concatenated boundary list, yields
\begin{tabbing}
\verb$  $\= \kill
\>\verb$[ ..., $$A$\verb$ := (a !# $$B$\verb$) * (-1), ..., $$A$\verb$ := (a !# $$C$\verb$) * (-1),$ \\
\>\verb$  ..., $$A$\verb$ := (a !# $$D$\verb$) * 1, ..., $$A$\verb$ := (a !# $$E$\verb$) * 1, ...]$
\end{tabbing}
Applying \verb$accumForm$ collects together all values associated with
the same site into a 0-Form.  Up to a possible re-ordering of
the list, the value of this form at a site $A$ will be
\begin{tabbing}
\verb$  $\= \kill
\>\verb$[(a !# $$B$\verb$) * (-1), (a !# $$C$\verb$) * (-1), (a !# $$D$\verb$) * 1,(a !# $$E$\verb$) * 1]$
\end{tabbing}
Mapping \verb$sum$ over this form
\begin{tabbing}
\verb$  $\= \kill
\>\verb$((a !# $$D$\verb$) - (a !# $$B$\verb$)) + ((a !# $$E$\verb$) - (a !# $$C$\verb$))$
\end{tabbing}
Apart from a slightly annoying overall sign, this is what one would expect
as a lattice approximation to
\begin{displaymath}
\frac{\partial a_x}{\partial x} + \frac{\partial a_y}{\partial y}.
\end{displaymath}

\subsection*{Example 3: general exterior and interior derivatives, and curl}

     The operations \verb$gradient$ and \verb$divergence$ generalise
immediately to give approximations to the general exterior
and interior derivatives:
\begin{verbatim}
  exterior :: (Cell a, Num b) => Form a b -> Form a b
  exterior a = form l r [i := cntrct a (boundary l i) | i <- cells l r]
               where l = lattice a
                     r = rank a + 1
\end{verbatim}
\begin{verbatim}
  interior :: (Cell a, Num b) => Form a b -> Form a b
  interior a = (fmap sum . accumForm l r . bndry l . fassocs) a
               where l = lattice a
                     r = rank a - 1
\end{verbatim}
As another concrete example, referring to figure \ref{incidence}
again, the exterior derivative of a 1-form \verb$a$ is a 2-form.
Evaluated on plaquette $A$ the result is
\begin{tabbing}
\verb$  $\= \kill
\>\verb$((a !# $$C$\verb$) - (a !# $$E$\verb$)) - ((a !# $$D$\verb$) - (a !# $$B$\verb$))$
\end{tabbing}
which is recognisably a lattice approximation to
\begin{displaymath}
\frac{\partial a_y}{\partial x} - \frac{\partial a_x}{\partial y}
\end{displaymath}
i.e., a component of the curl,
$\mbox{\boldmath $\nabla$} \mbox{\boldmath $\times$} \mbox{\boldmath $a$}$
(which is an antisymmetric 2-tensor in general).

\subsection*{\label{laplacian}Example 4: Laplacians}

     A generalised Laplacian operator can be implemented by
\begin{verbatim}
  laplacian :: (Cell a, Num b) => Form a b -> Form a b
  laplacian = interior . exterior
\end{verbatim}
The Laplacian maps an $r$-form to another $r$-form.
Its effect on a 0-form will be familiar---up to an overall
sign\footnote{This a peculiarity of our definition.  Appendix
\ref{sign} sheds some light on it.} the result
of applying it to $u$ is given at site $(x, y)$ by
the LHS of equation (\ref{discrete}).

    On a 1-form the result at a particular link is
a certain sum of signed values taken from the links on the plaquettes
surrounding that link, illustrated in figure \ref{staples}.

    To end this section we re-emphasize that the
implementations given here are independent of the details
of the lattice shape.  Any sensible lattice topology can be encoded
in the \verb$boundary$ function, which is why it is
a powerful concept.


\section{\label{sign}Why the Laplacian has the wrong sign}

If one adopts the shorthand
\begin{center}
\begin{tabular}{lcl}
$\partial c$     & $\sim$ & \verb$bndry l $$c$ \\
$(a, c)$         & $\sim$ & \verb$cntrct $$a$\verb$ $$c$ \\
$* a$            & $\sim$ & \verb$fassocs $$a$ \\
$* c$            & $\sim$ & \verb$fmap sum . accumForm l r $$c$ \\
$\mbox{\bf d} a$ & $\sim$ & \verb$exterior $$a$ \\
$\delta a$       & $\sim$ & \verb$interior $$a$ \\
\end{tabular}
\end{center}
where $a$ is a form and $c$ an association list,
taking values in some numeric type, the following
algebraic identities are satisfied
\begin{eqnarray*}
(\mbox{\bf d} a, c) & = & (a, \partial c) \\
\delta a                  & = & * \partial * a \\
{*}*a                     & = & a \\
(a, **c)                  & = & (a, c) \\
(*c, *a)                  & = & (a, c).
\end{eqnarray*}
Now, for any $a$
\begin{eqnarray*}
(\delta \mbox{\bf d} a, *a)
& = & (a, * \delta \mbox{\bf d} a) \\
& = & (a, \partial * \mbox{\bf d} a) \\
& = & (\mbox{\bf d} a, * \mbox{\bf d} a)
\end{eqnarray*}
This is a sum of squares.  So $\delta \mbox{\bf d}$,
adopted as the definition of the Laplacian above,
is a positive semi-definite operator.  The
conventional definition of the Laplacian is negative.

