INTERFACE SPMultiGrid;

(* A multigrid data structure for point location, approximaton, and PDE solution. *)

IMPORT SPTriang, SPOverlapTable, SPMatrix, SPPWFunction;
FROM SPTriang IMPORT TriangleNumber;

TYPE
  T = RECORD
      NL: CARDINAL;                           (* Number of grids except top grid. *)
      tri: REF ARRAY OF SPTriang.T;           (* The grids. *)
      dim: REF ARRAY OF CARDINAL;             (* Dim. of spline space for each grid. *)
      ovl: REF ARRAY OF REF SPOverlapTable.T; (* Face overlap tables. *)
      GGL: REF ARRAY OF SPMatrix.T;           (* Choleski decs. of rigidity matrices. *)
      FG: REF ARRAY OF SPMatrix.T;            (* Transfer matrices. *)
      SSL: REF ARRAY OF SPMatrix.T;           (* Choleski decs. of operator matrices. *)
    END;
    (*
      A multigrid structure "m" consists of "NL+1" irregular triangular
      grids, "m.tri[0..NL]".  The table "m.ovl[k]" lists the triangles
      of grid "m.tri[k+1]" that overlap each triangle of grid "m.tri[k]".
      
      Each level "k" of the multigrid has an associated space "V[k]"
      of spherical spline functions, and a fixed basis "B[k,r]" of 
      such functions, where "r" ranges in "[0..m.dim[k]-1]".
      
      The sparse matrix "m.GGL[k]" is the lower triangular Choleski
      decomposition of the rigidity matrix "GG[k]" for basis "B[k]".
      (Element "[i,j]" of "GG[k]", by definition, is the dot product
      of basis functions "B[k,i]" and "B[k,j]".)

      The sparse matrix "m.FG[k]" is the connection matrix
      between spline spaces of layers "k" and "k+1".
      Specifically, element "[i,j]" of "m.FG[k]" is the 
      dot product of "B[k,i]" and "B[k+1,j]". [???or the other way around???]
      
      The sparse matrix "m.SSL[k]" is the lower triangular Choleski
      decomposition of the operator matrix "SS[k] = HH[k] + c*GG[k]",
      where "GG[k]" is the rigidity matrix of "B[k]", defined above, and
      "HH[k]" is the rigdity matrix of the basis gradients.
      That is, element "[i,j]" of "HH[k]" is the dot product of
      "SGRAD(B[k,i])" and "SGRAD(B[k,j])", where "SGRAD" denotes 
      the spherical gradient operator. *)
    
  Level = CARDINAL;  (* Index "[0..m.NL]" of a layer in a multigrid "m". *)
  
  Point = SPPWFunction.Point;

PROCEDURE ReLocate(
    READONLY m: T; 
    k: Level; 
    fn: TriangleNumber; 
    p: Point;
  ): TriangleNumber;
  (*
    Given a point "p" in triangle "fn" in layer "k" of multigrid "m",
    returns the number of the triangle in layer "k+1" that contains "p". *)
    
PROCEDURE Locate(
    READONLY m: T; 
    p: Point;
    k: Level; 
  ): TriangleNumber;
  (*
    Returns the number of the triangle in layer "k" of "m" that contains "p". *)

PROCEDURE Write(
    gridName: TEXT;  (* Name OF the base grid. *)
    opTag: TEXT;     (* Label for a specific differential operator. *)
    READONLY m: T;
  );
  (*
    Writes a multigrid structure to disk, as a set of files
    whose names begin with "gridName" or "gridName-opTag". *)
    
PROCEDURE Read(
    gridName: TEXT;  (* Name OF the base grid. *)
    NL: CARDINAL;    (* Number of grids (excluding top layer). *)
    opTag: TEXT;     (* Label for a specific differential operator. *)
  ): T;
  (*
    Read a multigrid structure from disk, from a set of files
    whose names begin with "gridName" or "gridName-opTag". *)
    
PROCEDURE Solve(
    READONLY m: T;
    READONLY f: ARRAY OF REF SPPWFunction.BasisCoeffs;
    maxIterations: CARDINAL;
    VAR (*IO*) u: ARRAY OF REF SPPWFunction.BasisCoeffs;
  );
  (*
    Solves the partial differential equation "Op(u) = f"
    over the space of spherical splines associated
    with each layer of the multigrid "m". 
    Assumes the matrices in the multigrid structure 
    were computed from the variational formulation of "Op". *)
  
END SPMultiGrid.