MODULE SPOrthoBasis;

IMPORT SPTriang, LR3, Wr, Fmt, SPHomoLabel, SPQuad,  LR3x3;
IMPORT SPPWFunction, SPBezierPWFunction, SPBezPolyPWFunction, Thread;
      
FROM Stdio IMPORT stderr;
FROM Math IMPORT sqrt;
FROM SPQuad IMPORT Lnext, Onext, EdgeNum, Sym, Arc;
FROM SPPWFunction  IMPORT ComputeSupportingPlane, Triangles, SortTriangles, Id;
FROM SPTriang IMPORT Left, Dest, Degree; 

TYPE
  ControlValues = SPBezierPWFunction.ControlValues;
    (* 
      Bezier Control values for one triangle.
      The order is defined in SPHomoLabel.i3; for degree 6 it is
         C600, C510, C501, C420, C411, C402, ... , C006
    *)
  TriangleData = SPBezierPWFunction.TriangleData;

<* FATAL Wr.Failure, Thread.Alerted *>
  
PROCEDURE BuildHBasis(
    deg: NAT;       (* Degree of polynomials *)
    cont: [0..1];        (* Continuity. *)
    tri: REF SPTriang.T; (* Triangulation. *)
    ortho: BOOL;         (* TRUE to quasi-orthogonalize the basis. *)
  ): REF SPBezierPWFunction.Basis =
  
  VAR nb: NAT;
  BEGIN
    WITH 
      NV = NUMBER(tri.out^),
      NT = NUMBER(tri.side^),
      NE = NUMBER(tri.arc^) DIV 2,
      (* The basis: *)
      NB = SpaceDimension(deg, cont, NT),
      rb = NEW(REF SPBezierPWFunction.Basis, NB), b = rb^,
      (* Index tables: *)
      fBas = NEW(REF ARRAY OF NAT, NT+1)^,
      eBas = NEW(REF ARRAY OF NAT, NE+1)^,
      vBas = NEW(REF ARRAY OF NAT, NV+1)^
      (*
        The basis elements associated with face number "fn" are
        b[fBas[fn]..fBas[fn+1]-1]
      *)
    DO 
      nb := 0;
      BuildFaceBasis(ortho, b, nb, tri, deg, cont, fBas);
      BuildEdgeBasis(ortho, b, nb, tri, deg, cont, fBas, eBas);
      BuildVertexBasis(ortho, b, nb, tri, deg, cont, fBas, eBas, vBas);
      RETURN rb
    END;
  END BuildHBasis;
  
PROCEDURE BuildNHBasis(
    deg: CARDINAL;
    cont: CARDINAL;
    tri: REF SPTriang.T;
    ortho: BOOLEAN;
  ): REF SPBezPolyPWFunction.Basis =
  BEGIN
    IF cont = 0 THEN 
      <* ASSERT deg = 3 *>
    ELSIF cont = 1 THEN
      <* ASSERT deg = 6 *>
    ELSE
      <* ASSERT FALSE *>
    END;
    WITH
      b0 = BuildHBasis(deg,   cont, tri, ortho)^, dim0 = NUMBER(b0),
      b1 = BuildHBasis(deg-1, cont, tri, ortho)^, dim1 = NUMBER(b1),
      rb = NEW(REF SPBezPolyPWFunction.Basis, dim0 + dim1), b = rb^
    DO
      FOR i := 0 TO LAST(b0) DO 
        b[i]      := SPBezPolyPWFunction.FromHomo(b0[i], deg)
      END;
      FOR i := 0 TO LAST(b1) DO
        b[i+dim0] := SPBezPolyPWFunction.FromHomo(b1[i], deg)
      END;
      RETURN rb
    END
  END BuildNHBasis;
  
PROCEDURE SpaceDimension(deg: NAT; cont: [0..1]; NT: NAT): NAT =
  BEGIN
    <* ASSERT NT MOD 2 = 0 *>
    RETURN (deg*deg - 3*deg*cont + 2*cont*cont)*NT DIV 2 + cont*cont + 3*cont + 2
  END SpaceDimension;
  
PROCEDURE NumCoeffs(deg: NAT): NAT = 
  (* 
    Number of Bezier coefficients per triange in 
    a homogeneous polynomial of degree "d".
  *)
  BEGIN
    RETURN (deg+1)*(deg+2) DIV 2
  END NumCoeffs;

PROCEDURE BuildFaceBasis(
    ortho: BOOL;
    VAR b : SPBezierPWFunction.Basis; (*OUT*)
    VAR nb: NAT; (* IN/OUT *)
    tri: REF SPTriang.T;
    deg: NAT;
    cont: [0..1];
    VAR fBas: ARRAY OF NAT; (* OUT *)
  )=
  (*
    Builds the basis elements associated with each face.
    Appends them to the basis "b" starting at "b[nb]".
    Sets "fBas[i]" to the index of the first element 
    associated with each face "i".  Also sets "fBas[NF}"
    to the index of the first element of "b" following
    all face elements.
  *)
  BEGIN
    WITH  
      side = tri.side^,
      NF = NUMBER(side),
      NC = NumCoeffs(deg)
    DO
      FOR i := 0 TO NF-1 DO 
        Wr.PutText(stderr, "face " & Fmt.Int(i) & "  nb = " & Fmt.Int(nb) & "\n");
        fBas[i] := nb;
        IF (cont = 0 AND deg > 2) OR (cont = 1 AND deg > 5) THEN
          WITH 
            elem = NEW(SPBezierPWFunction.T),
            e = side[i]
          DO
            Wr.PutText(stderr, "BuildFaceBasis(face = " & Fmt.Int(i) & ")\n");
            b[nb] := elem;
            elem.tri:= tri;
            elem.data:= MakeTriangles(1, NC);
            WITH 
              t0 = NARROW(elem.data[0], TriangleData)
            DO
              t0.face:= Left(e).num;
              t0.deg:= deg;
              WITH c0 = t0.c^ DO
                IF cont = 0 AND deg = 3 THEN 
                  c0[4] := 1.0d0
                ELSIF cont = 1 AND deg = 6 THEN
                  c0[12] := 1.0d0
                ELSE
                  <* ASSERT FALSE *>
                END
              END;
            END;
            elem.supp:= ComputeSupportingPlane(elem.data^, side);
            IF ortho THEN Normalize(elem) END;
          END;
          INC(nb);
        END
      END;
      fBas[NF] := nb;
    END;
  END BuildFaceBasis;
  
PROCEDURE ComputeEdgeContinuityCoeffs(e: Arc; tri: REF SPTriang.T): LR3.T =
  (*
    Computes "bu", "bv", "bw" such that "z = bu*u + bv*v + bw*w",
    where
      u = Dest(Onext(e)).c,
      v = Org(e).c,
      w = Dest(e).c,
      z = Dest(Oprev(e)).c
  *)
  BEGIN
    WITH
      z = Dest(Onext(Sym(e))).c,
      face = Left(e).num,
      a = tri.side[face],
      c2b = tri.c2b[face],
      bar = LR3x3.MapCol(c2b, z)
    DO
      IF e = a THEN 
        RETURN bar
      ELSIF e = Lnext(a) THEN
        RETURN LR3.T{bar[1], bar[2], bar[0]}
      ELSIF a = Lnext(e) THEN
        RETURN LR3.T{bar[2], bar[0], bar[1]}
      ELSE
        <* ASSERT FALSE *>
      END
    END
  END ComputeEdgeContinuityCoeffs;
  
PROCEDURE SetCoeffsEdge(
    e: Arc;                 (* Reference edge for basis element *)
    tri: REF SPTriang.T;    (* Triangulation *)
    deg: NAT;               (* Degree of polynomial *)
    VAR c: ControlValues;   (* Coeffs vector *)
    ind: NAT;               (* Coeff index relative to "e" *)
    v: LONGREAL;            (* Coefficient VALUE *)
  ) =
  (*
    Sets a Bezier coefficient "c[k]" with indices "ind" 
    relative to edge "e", which must be a side of the 
    corresponding triangle, but not necessarily the 
    reference edge of that triangle. 
  *)
  VAR indr: NAT;
      ijk: SPHomoLabel.T;
  BEGIN
    WITH 
      f = Left(e).num,
      a = tri.side[f]
    DO 
      IF e = a THEN
        c[ind] := v
      ELSIF Lnext(a) = e  THEN 
        ijk:= SPHomoLabel.IndexToLabel(ind,deg);
        ijk := SPHomoLabel.T{ijk[2], ijk[0], ijk[1]};
        indr:= SPHomoLabel.LabelToIndex(ijk,deg);
        c[indr] := v
      ELSIF  Lnext(e) = a THEN
        ijk := SPHomoLabel.IndexToLabel(ind,deg);
        ijk := SPHomoLabel.T{ijk[1], ijk[2], ijk[0]};
        indr:= SPHomoLabel.LabelToIndex(ijk, deg);
        c[indr] := v;
      ELSE
        <* ASSERT FALSE *>
      END;
    END;
  END SetCoeffsEdge;


PROCEDURE BuildEdgeBasis(
    ortho: BOOL;
    VAR b : SPBezierPWFunction.Basis; 
    VAR nb: NAT; 
    tri: REF SPTriang.T;
    deg: NAT;
    cont: [0..1];
    READONLY fBas: ARRAY OF NAT;
    VAR eBas: ARRAY OF NAT;
  ) =
  (* 
    Builds the basis elements associated with each edge.
    Appends them to the basis "b" starting at "b[nb]".
    Sets "eBas[i]" to the index of the first element 
    associated with edge number "i".  Also sets "eBas[NE}"
    to the index of the first element of "b" following
    all edge elements.
  *)
  BEGIN
    WITH 
      arc = tri.arc^,
      NE = NUMBER(arc) DIV 2
    DO
      FOR i := 0 TO NE-1 DO 
        WITH e = arc[2*i] DO
          Wr.PutText(stderr, "edge " & Fmt.Int(i) & "  nb = " & Fmt.Int(nb) & "\n");
          eBas[i] := nb;
          IF cont = 0 THEN
            IF deg = 2 THEN
              BuilEdgeBasis2(e, ortho, b, nb, tri)
            ELSIF deg = 3 THEN
              BuilEdgeBasis3(e, ortho, b, nb, tri, fBas)
            ELSE
              <* ASSERT FALSE *>
            END
          ELSIF cont = 1 THEN
            IF deg = 5  THEN 
              BuilEdgeBasis5(e, ortho, b, nb, tri)
            ELSIF deg = 6 THEN
              BuilEdgeBasis6(e, ortho, b, nb, tri, fBas)
            ELSE
              <* ASSERT FALSE *>
            END
          ELSE
            <* ASSERT FALSE *>
          END
        END
      END;
      eBas[NE] := nb;
    END
  END BuildEdgeBasis;
  
PROCEDURE BuilEdgeBasis2(
    e: Arc;
    ortho: BOOL;
    VAR b : SPBezierPWFunction.Basis; 
    VAR nb: NAT; 
    tri: REF SPTriang.T;
  )=
  (*
    Builds the basis elements associated with the edge "e",
    for degree 2 and continuity 0.
  *)
  BEGIN
    WITH 
      start = nb + 0,  (* Remember the first element of this edge. *)
      elem = NEW(SPBezierPWFunction.T)
    DO
      Wr.PutText(stderr, "BuildEdgeBasis2(edge = " & Fmt.Int(EdgeNum(e)) & ")\n");
      b[nb] := elem;
      elem.tri:= tri;
      elem.data:= MakeTriangles(2, 6);
      WITH 
        t0 = NARROW(elem.data[0], TriangleData),
        t1 = NARROW(elem.data[1], TriangleData)
      DO
        t0.face:= Left(e).num;
        t0.deg:= 2;
        WITH c0 = t0.c^ DO
          SetCoeffsEdge(e, tri, 2, c0, 4, 1.0d0) 
        END;  
        t1.face:= Left(Sym(e)).num;
        t1.deg:= 2;
        WITH c1 = t1.c^ DO
          SetCoeffsEdge(Sym(e), tri, 2, c1, 4, 1.0d0); 
        END;
      END;
      elem.supp:= ComputeSupportingPlane(elem.data^, tri.side^);
      INC(nb);
      IF ortho THEN
        Wr.PutText(stderr, "orthogonalizing...\n");
        (* There are no face elements, and only one edge element. Just Normalize it. *)
        <* ASSERT nb = start + 1 *>
        Normalize(b[start]);
      END;
    END;
  END BuilEdgeBasis2;

PROCEDURE BuilEdgeBasis3(
    e: Arc;
    ortho: BOOL;
    VAR b: SPBezierPWFunction.Basis; 
    VAR nb: NAT; 
    tri: REF SPTriang.T;
    READONLY fBas: ARRAY OF NAT;
  ) =
  BEGIN
    WITH start = nb + 0 DO
      BuilEdgeBasis31(e, b, nb, tri);
      BuilEdgeBasis32(e, b, nb, tri);
      IF ortho THEN
        Wr.PutText(stderr, "orthogonalizing...\n");
        (* There are two edge elements, both on the edge, and one face
           element in each triangle.  Make the edge elements orthogonal
           to the triangle ones, then orthogonal to each other.
        *)
        <* ASSERT nb = start+2 *>
        WITH
          fn = Left(e).num,
          f = b[fBas[fn]],
          gn = Left(Sym(e)).num, 
          g = b[fBas[gn]],
          e1 = b[start],
          e2 = b[start+1]
        DO
          Orthize(e1, f); Orthize(e1, g);
          Orthize(e2, f); Orthize(e2, g);
          OrthoNormalizePair(e1, e2)
        END
      END
    END;
  END BuilEdgeBasis3;

PROCEDURE BuilEdgeBasis31(
    e: Arc;
    VAR b : SPBezierPWFunction.Basis; 
    VAR nb: NAT; 
    tri: REF SPTriang.T;
  ) =
  BEGIN 
    WITH 
      elem = NEW(SPBezierPWFunction.T)
    DO
      Wr.PutText(stderr, "BuildEdgeBasis31(edge = " & Fmt.Int(EdgeNum(e)) & ")\n");
      b[nb] := elem;
      elem.tri:= tri;
      elem.data:= MakeTriangles(2, 10);
      WITH 
        t0 = NARROW(elem.data[0], TriangleData),
        t1 = NARROW(elem.data[1], TriangleData)
      DO
        t0.deg:= 3;
        t0.face:= Left(e).num;
        WITH c0 = t0.c^ DO
          SetCoeffsEdge(e, tri, 3, c0, 7, 1.0d0) 
        END;  
        t1.deg:= 3;
        t1.face:= Left(Sym(e)).num;
        WITH c1 = t1.c^ DO
          SetCoeffsEdge(Sym(e), tri, 3, c1, 8, 1.0d0); 
        END;
      END;
      elem.supp:= ComputeSupportingPlane(elem.data^, tri.side^);
      INC(nb);
    END;
  END BuilEdgeBasis31;

PROCEDURE BuilEdgeBasis32(
    e: Arc;
    VAR b : SPBezierPWFunction.Basis; 
    VAR nb: NAT; 
    tri: REF SPTriang.T;
  )=
  BEGIN 
    WITH 
      elem = NEW(SPBezierPWFunction.T)
    DO
      Wr.PutText(stderr, "BuildEdgeBasis32(edge = " & Fmt.Int(EdgeNum(e)) & ")\n");
      b[nb] := elem;
      elem.tri:= tri;
      elem.data:= MakeTriangles(2, 10);
      WITH 
        t0 = NARROW(elem.data[0], TriangleData),
        t1 = NARROW(elem.data[1], TriangleData)
      DO
        t0.deg := 3;
        t0.face:= Left(e).num;
        WITH c0 = t0.c^ DO
          SetCoeffsEdge(e, tri, 3, c0, 8, 1.0d0) 
        END;  
        t1.deg:= 3;
        t1.face:= Left(Sym(e)).num;
        WITH c1 = t1.c^ DO
          SetCoeffsEdge(Sym(e), tri, 3, c1, 7, 1.0d0); 
        END;
      END;
      elem.supp:= ComputeSupportingPlane(elem.data^, tri.side^);
      INC(nb);
    END;    
  END BuilEdgeBasis32;


PROCEDURE BuilEdgeBasis5(
    e: Arc;
    ortho: BOOL;
    VAR b : SPBezierPWFunction.Basis; 
    VAR nb: NAT; 
    tri: REF SPTriang.T;
  )=
  (*
    Builds the basis elements associated with the edge "e",
    for degree 5 and continuity 1.
  *)
  BEGIN
    WITH 
      start = nb + 0,  (* Remember the first element of this edge. *)
      elem = NEW(SPBezierPWFunction.T),
      bar = ComputeEdgeContinuityCoeffs(e, tri)
    DO
      Wr.PutText(stderr, "BuildEdgeBasis5(edge = " & Fmt.Int(EdgeNum(e)) & ")\n");
      b[nb] := elem;
      elem.tri:= tri;
      elem.data:= MakeTriangles(2, 21);
      WITH 
        t0 = NARROW(elem.data[0], TriangleData),
        t1 = NARROW(elem.data[1], TriangleData)
      DO
        t0.face:= Left(e).num;
        t0.deg:= 5;
        WITH c0 = t0.c^ DO
          SetCoeffsEdge(e, tri, 5, c0, 12, 1.0d0) 
        END;  
        t1.face:= Left(Sym(e)).num;
        t1.deg:= 5;
        WITH c1 = t1.c^ DO
          SetCoeffsEdge(Sym(e), tri, 5, c1, 12, bar[0]); 
        END;
      END;
      elem.supp:= ComputeSupportingPlane(elem.data^, tri.side^);
      INC(nb);
      IF ortho THEN
        Wr.PutText(stderr, "orthogonalizing...\n");
        (* There are no face elements, and only one edge element. Just Normalize it. *)
        <* ASSERT nb = start + 1 *>
        Normalize(b[start]);
      END;
    END;
  END BuilEdgeBasis5;

PROCEDURE BuilEdgeBasis6(
    e: Arc;
    ortho: BOOL;
    VAR b: SPBezierPWFunction.Basis; 
    VAR nb: NAT; 
    tri: REF SPTriang.T;
    READONLY fBas: ARRAY OF NAT;
  ) =
  BEGIN
    WITH 
      start = nb + 0  (* Remember the first element of this edge. *)
    DO
      BuilEdgeBasis61(e, b, nb, tri);
      BuilEdgeBasis62(e, b, nb, tri);
      BuilEdgeBasis63(e, b, nb, tri);
      IF ortho THEN
        Wr.PutText(stderr, "orthogonalizing...\n");
        (*
          We have three edge elements (one on the edge and two outside
          it), and one face element on each triangle. First, make the edge
          elements orthogonal to the triangle ones.  Then make the off-edge
          elements orthogonal to each other.  Finally make the on-edge
          element orthogonal to the off-edge elements.
        *)
        <* ASSERT nb = start+3 *>
        WITH
          fn = Left(e).num,
          f = b[fBas[fn]],        (* Face element on left triangle *)
          gn = Left(Sym(e)).num,
          g = b[fBas[gn]],        (* Face element on right triangle *)
          e1 = b[start],          (* First off-edge element *)
          e2 = b[start+1],        (* Second off-edge element *)
          e3 = b[start+2]         (* On-edge element *)
        DO
          Orthize(e1, f); Orthize(e1, g);
          Orthize(e2, f); Orthize(e2, g);
          Orthize(e3, f); Orthize(e3, g);
          OrthoNormalizePair(e1, e2);
          Orthize(e3, e1); Orthize(e3, e2);
          Normalize(e3);
        END
      END
    END
  END BuilEdgeBasis6;

PROCEDURE BuilEdgeBasis61(
    e: Arc;
    VAR b : SPBezierPWFunction.Basis; 
    VAR nb: NAT; 
    tri: REF SPTriang.T;
  ) =
  BEGIN 
    WITH 
      elem = NEW(SPBezierPWFunction.T),
      bar = ComputeEdgeContinuityCoeffs(e, tri)
    DO
      Wr.PutText(stderr, "BuildEdgeBasis61(edge = " & Fmt.Int(EdgeNum(e)) & ")\n");
      b[nb] := elem;
      elem.tri:= tri;
      elem.data:= MakeTriangles(2, 28);
      WITH 
        t0 = NARROW(elem.data[0], TriangleData),
        t1 = NARROW(elem.data[1], TriangleData)
      DO
        t0.deg:= 6;
        t0.face:= Left(e).num;
        WITH c0 = t0.c^ DO
          SetCoeffsEdge(e, tri, 6, c0, 17, 1.0d0) 
        END;  
        t1.deg:= 6;
        t1.face:= Left(Sym(e)).num;
        WITH c1 = t1.c^ DO
          SetCoeffsEdge(Sym(e), tri, 6, c1, 18, bar[0]); 
        END;
      END;
      elem.supp:= ComputeSupportingPlane(elem.data^, tri.side^);
      INC(nb);
    END;
  END BuilEdgeBasis61;

PROCEDURE BuilEdgeBasis62(
    e: Arc;
    VAR b : SPBezierPWFunction.Basis; 
    VAR nb: NAT; 
    tri: REF SPTriang.T;
  )=
  BEGIN 
    WITH 
      elem = NEW(SPBezierPWFunction.T),
      bar = ComputeEdgeContinuityCoeffs(e, tri)
    DO
      Wr.PutText(stderr, "BuildEdgeBasis62(edge = " & Fmt.Int(EdgeNum(e)) & ")\n");
      b[nb] := elem;
      elem.tri:= tri;
      elem.data:= MakeTriangles(2, 28);
      WITH 
        t0 = NARROW(elem.data[0], TriangleData),
        t1 = NARROW(elem.data[1], TriangleData)
      DO
        t0.deg := 6;
        t0.face:= Left(e).num;
        WITH c0 = t0.c^ DO
          SetCoeffsEdge(e, tri, 6, c0, 18, 1.0d0) 
        END;  
        t1.deg:= 6;
        t1.face:= Left(Sym(e)).num;
        WITH c1 = t1.c^ DO
          SetCoeffsEdge(Sym(e), tri, 6, c1, 17, bar[0]); 
        END;
      END;
      elem.supp:= ComputeSupportingPlane(elem.data^, tri.side^);
      INC(nb);
    END;    
  END BuilEdgeBasis62;

PROCEDURE BuilEdgeBasis63(
    e: Arc;
    VAR b : SPBezierPWFunction.Basis; 
    VAR nb: NAT; 
    tri: REF SPTriang.T;
  )=
  BEGIN 
    WITH 
      elem = NEW(SPBezierPWFunction.T),
      bar = ComputeEdgeContinuityCoeffs(e, tri)
    DO
      Wr.PutText(stderr, "BuildEdgeBasis63(edge = " & Fmt.Int(EdgeNum(e)) & ")\n");
      b[nb] := elem;
      elem.tri:= tri;
      elem.data:= MakeTriangles(2, 28);
      WITH 
        t0 = NARROW(elem.data[0], TriangleData),
        t1 = NARROW(elem.data[1], TriangleData)
      DO
        t0.deg:= 6;
        t0.face:= Left(e).num;
        WITH c0 = t0.c^ DO
          SetCoeffsEdge(e, tri, 6, c0, 24, 1.0d0); 
        END;  
        t1.deg:= 6;
        t1.face:= Left(Sym(e)).num;
        WITH c1 = t1.c^ DO
          SetCoeffsEdge(Sym(e), tri, 6, c1, 24, 1.0d0); 
          SetCoeffsEdge(Sym(e), tri, 6, c1, 18, bar[2]); 
          SetCoeffsEdge(Sym(e), tri, 6, c1, 17, bar[1]); 
        END;
      END;
      elem.supp:= ComputeSupportingPlane(elem.data^, tri.side^);
      INC(nb);
    END;   
  END BuilEdgeBasis63;


PROCEDURE BuildVertexBasis(
    ortho: BOOL;
    VAR b : SPBezierPWFunction.Basis; 
    VAR nb: NAT; 
    tri: REF SPTriang.T;
    deg: NAT;
    cont: [0..1];
    READONLY fBas: ARRAY OF NAT;
    READONLY eBas: ARRAY OF NAT;
    VAR vBas: ARRAY OF NAT;
  ) =
  (* 
    Builds the basis elements associated with each vertex.
    Appends them to the basis "b" starting at "b[nb]".
    Sets "vBas[i]" to the index of the first element 
    associated with vertex number "i".  Also sets "vBas[NV}"
    to the index of the first element of "b" following
    all vertex elements.
  *)
  BEGIN
    WITH
      out = tri.out^,
      NV = NUMBER(out)
    DO 
      FOR i:= 0 TO NV-1 DO 
        Wr.PutText(stderr, "vertex " & Fmt.Int(i) & " nb = " & Fmt.Int(nb) & " ===\n");
        WITH e = out[i] DO
          vBas[i] := nb;
          IF cont = 0 THEN
            BuildVertexBasis0(ortho, deg, e, b, nb, tri, fBas, eBas)
          ELSIF cont = 1 THEN
            BuildVertexBasis1(ortho, deg, e, b, nb, tri, fBas, eBas)
          ELSE
            <* ASSERT FALSE *>
          END
        END
      END;
      vBas[NV] := nb;
    END
  END BuildVertexBasis;
  

PROCEDURE BuildVertexBasis0(
    ortho: BOOL;
    deg: NAT;
    e: Arc;
    VAR b : SPBezierPWFunction.Basis; 
    VAR nb: NAT; 
    tri: REF SPTriang.T;
    READONLY fBas: ARRAY OF NAT;
    READONLY eBas: ARRAY OF NAT;
  )=
  (*
    Builds the basis elements associated to vertex "v" for 
    splines with zero-order continuity.
  *)
  VAR a: Arc;
      v: ARRAY [0..0] OF LONG;
  BEGIN
    WITH 
      start = nb + 0,
      NC = NumCoeffs(deg),
      elem = NEW(SPBezierPWFunction.T),
      order = SPTriang.Degree(e)
    DO 
      b[nb] := elem;
      elem.tri := tri;
      elem.data:= MakeTriangles(order, NC);
      v[0] := 1.0d0;
      a := e; 
      FOR r := 0 TO order-1 DO 
        WITH 
          t0 = NARROW(elem.data[r], TriangleData)
        DO 
          t0.face := Left(a).num;
          t0.deg := deg;
          SetCoeffsVertex(a, tri, deg, t0.c^, v);
          a := Onext(a);
        END;
      END;
      elem.supp := ComputeSupportingPlane(elem.data^, tri.side^);
      INC(nb);
      IF ortho THEN
        Wr.PutText(stderr, "orthogonalizing...\n");
        a := e; 
        FOR r := 0 TO order-1 DO 
          WITH 
            en = EdgeNum(a),
            fn = Left(a).num
          DO 
            (* Orthogonalize b[start] against the face elements: *)
            FOR k := fBas[fn] TO fBas[fn+1]-1 DO
              Orthize(b[start], b[k])
            END;
            (* Orthogonalize b[start] against the edge elements: *)
            FOR k := eBas[en] TO eBas[en+1]-1 DO
              Orthize(b[start], b[k])
            END;
          END
        END;
      END;
    END;
  END BuildVertexBasis0;

 
PROCEDURE BuildVertexBasis1(
    ortho: BOOL;
    deg: NAT;
    e: Arc;
    VAR b : SPBezierPWFunction.Basis; 
    VAR nb: NAT; 
    tri: REF SPTriang.T;
    READONLY fBas: ARRAY OF NAT;
    READONLY eBas: ARRAY OF NAT;
  )=
 
  PROCEDURE ArcIsEdgeReference(e: Arc): BOOLEAN =
    (*
      TRUE if the basis elements associated with the edge
      "e.quad" have used "e" (and not "Sym(e)) as the reference arc. *)
    BEGIN
      RETURN tri.arc[2*EdgeNum(e)] = e
    END ArcIsEdgeReference;
  
  PROCEDURE VC7BelongsToGamma(e: Arc): BOOLEAN =
    (*
      Arc "e" is an arc out of the vertex; we are 
      considering the coefficients of the triangle "Left(e)". *)
    BEGIN 
      RETURN ArcIsEdgeReference(e)
    END VC7BelongsToGamma;
  
  PROCEDURE VC8BelongsToGamma(e: Arc): BOOLEAN =
    (*
      Arc "e" is an arc out of the vertex; we are 
      considering the coefficients of the triangle "Left(e)". *)
    BEGIN 
      RETURN ArcIsEdgeReference(Sym(Onext(e)))
    END VC8BelongsToGamma;
  
  PROCEDURE ComputeVC7AndVC8(e: Arc; VAR v:  ARRAY[0..9] OF LONGREAL;) =
    VAR bar: ARRAY[0..2] OF LONGREAL;
    
    BEGIN
      v[6]:= 0.0d0;
      v[9]:= 0.0d0;
      IF  VC7BelongsToGamma(e) THEN 
        v[7]:= 0.0d0
      ELSE 
        bar := ComputeEdgeContinuityCoeffs(Sym(e), tri);
        v[7]:= bar[2]*v[3]
      END;
      IF  VC8BelongsToGamma(e) THEN 
        v[8]:= 0.0d0
      ELSE 
        bar := ComputeEdgeContinuityCoeffs(Onext(e), tri);
        v[8]:= bar[1]*v[5]
      END;
    END ComputeVC7AndVC8;
  
  PROCEDURE Elem(e: Arc; k: NAT) =
    (*
      Generates the "k"th element of the basis elements 
      associated with vertex "Org(e)". *) 
    VAR a: Arc;
        bar: ARRAY[0..2] OF LONGREAL;
        v, w:  ARRAY[0..9] OF LONGREAL;
    BEGIN
      WITH 
        elem = NEW(SPBezierPWFunction.T),
        order = SPTriang.Degree(e)
      DO 
        b[nb] := elem;
        elem.tri := tri;
        elem.data:= MakeTriangles(order, NumCoeffs(deg));
        FOR i := 0 TO 9 DO v[i] := 0.0d0 END;
        IF k <= 5 THEN v[k] := 1.0d0 END;
        bar := ComputeEdgeContinuityCoeffs(e, tri);
        FOR i := 0 TO 9 DO w[i] := 0.0d0 END;
        w[4] := v[4]*bar[0] + bar[1]*v[1] + bar[2]*v[3];
        w[1] := bar[0]*v[2] + bar[1]*v[0] + bar[2]*v[1];
        w[0] := v[0];
        w[2] := v[1];
        w[5] := v[3];
        a := e; 
        
        FOR r := 1 TO order-2 DO 
          WITH 
            t0 = NARROW(elem.data[r-1], TriangleData)
          DO 
            t0.face := Left(a).num;
            t0.deg := deg;
            ComputeVC7AndVC8(a, v);
            SetCoeffsVertex(a, tri, deg, t0.c^, v);
            a := Onext(a);
          END;
          
          bar := ComputeEdgeContinuityCoeffs(Sym(a), tri);
          
          WITH
            nv2 = bar[0]*v[1] + bar[1]*v[2] + bar[2]*v[0],
            nv4 = bar[0]*v[4] + bar[1]*v[5] + bar[2]*v[2]
          DO
            v[1] := v[2];
            v[2] := nv2;
            v[3] := v[5];
            v[4] := nv4;
            IF k = r+5 THEN v[5] := 1.0d0 ELSE v[5] := 0.0d0 END
          END;
        END;
        
        bar := ComputeEdgeContinuityCoeffs(Sym(Onext(a)), tri);
        v[5] := (w[4] - bar[0]*v[4] - bar[2]*v[2])/bar[1];
        w[3] := v[5];
             
        WITH 
          t0 = NARROW(elem.data[order-2], TriangleData)
        DO 
          t0.face := Left(a).num;
          t0.deg := deg;
          ComputeVC7AndVC8(a, v);
          SetCoeffsVertex(a, tri, deg, t0.c^, v);
        END;
        
        a := Onext(a);
        WITH 
          t0 = NARROW(elem.data[order-1], TriangleData)
        DO 
          t0.face := Left(a).num;
          t0.deg := deg;
          ComputeVC7AndVC8(a, w);
          SetCoeffsVertex(a, tri, deg, t0.c^, w);
        END;
        elem.supp := ComputeSupportingPlane(elem.data^, tri.side^);
      END;
    INC(nb);
    END Elem;
    
  VAR a: Arc;
  BEGIN
    WITH 
      start = nb + 0,
      vDegree = Degree(e)
    DO
      FOR k := 0 TO vDegree + 2 DO 
        Elem(e, k);
      END;
      IF ortho THEN
        Wr.PutText(stderr, "orthogonalizing...\n");
        (* Orthogonalize the vertex elements against the face and edge elements: *)
        a := e; 
        FOR r := 0 TO vDegree-1 DO 
          WITH 
            en = EdgeNum(a),
            fn = Left(a).num
          DO 
            FOR i := start TO nb-1 DO 
              (* Orthogonalize b[i] against the face elements: *)
              FOR k := fBas[fn] TO fBas[fn+1]-1 DO
                Orthize(b[i], b[k])
              END;
              (* Orthogonalize b[start] against the edge elements: *)
              FOR k := eBas[en] TO eBas[en+1]-1 DO
                Orthize(b[i], b[k])
              END;
            END;
          END
        END;
        (* Orthogonalize the vertex elements among themselves. *)
        (* Orthogonalize elements 1 and 2 of this vertex: *)
        Normalize(b[start+1]);
        Orthize(b[start+2], b[start+1]);
        Normalize(b[start+2]);
        (* Orthogonalize element 0 with respect ot 1 and 2: *)
        Orthize(b[start], b[start+1]);
        Orthize(b[start], b[start+2]);
        Normalize(b[start]);
        (* Orthogonalize the fringe elements with respect to 0,1,2: *)
        FOR i := start+3 TO nb-1 DO
          Orthize(b[i], b[start]);
          Orthize(b[i], b[start+1]);
          Orthize(b[i], b[start+2]);
        END;
        (* Orthogonalize the fringe elements among themselves: *)
        FOR i := start+3 TO nb-1 DO
          FOR j := start+3 TO i-1 DO
            Orthize(b[i], b[j]);
          END;
          Normalize(b[i]);
        END;
      END;
    END;
 END BuildVertexBasis1;


PROCEDURE SetCoeffsVertex(
    e: Arc;                 (* Reference edge for basis element *)
    tri: REF SPTriang.T;    (* Triangulation *)
    deg: NAT;               (* Degree of polynomial *)
    VAR c: ControlValues;
    READONLY v: ARRAY OF LONGREAL;
  ) =
  (*
    Stores into "c" the vertex-associated Bezier coefficients 
    "v[k]", relative to the edge "e", which is a side of 
    the corresponding triangle, but not necessarily its 
    reference edge.

    Assumes that "v[0]" is the coefficient associated to the node "Org(e)",
    "v[1]" and "v[3]" (if they exist) lie on "e", "v[2]" and "v[5]"
    lie on "Onext(e)", and "v[4]" between "v[3]"
    and "v[5]". *)

  VAR indr: NAT;
      ijk: SPHomoLabel.T;

  BEGIN
    WITH 
      f = Left(e).num,
      a = tri.side[f]
    DO 
      IF  a = Lnext(e) THEN
        FOR i:= 0 TO LAST(v) DO
          c[i]:= v[i]
        END;
      ELSIF e = Lnext(a)  THEN 
        FOR i:= 0 TO LAST(v) DO 
          ijk := SPHomoLabel.IndexToLabel(i,deg);
          ijk := SPHomoLabel.T{ijk[1], ijk[2], ijk[0]};
          indr:= SPHomoLabel.LabelToIndex(ijk,deg);
          c[indr]:= v[i]
        END;
      ELSIF  a = e  THEN 
        FOR i:= 0 TO LAST(v) DO 
          ijk:= SPHomoLabel.IndexToLabel(i,deg);
          ijk := SPHomoLabel.T{ijk[2], ijk[0], ijk[1]};
          indr:= SPHomoLabel.LabelToIndex(ijk,deg);
          c[indr]:= v[i]
        END;
      END;
    END;
  END SetCoeffsVertex;    

PROCEDURE MakeTriangles(n: CARDINAL; ncv: CARDINAL): REF Triangles =
  (*
    Creates "n" TriagleData objects and associated sub-objects.
    Also sets all coefficients to zero. *)
  BEGIN
    WITH
      r = NEW(REF Triangles, n),
      t = r^
    DO
      FOR i := 0 TO n-1 DO 
        WITH ti = NEW(TriangleData) DO
          t[i] := ti;
          ti.barycentric := TRUE;
          ti.c := NEW(REF ControlValues, ncv);
          WITH c = ti.c^ DO 
            FOR i := 0 TO ncv-1 DO c[i] := 0.0d0 END
          END;
        END
      END;
      RETURN r
    END
  END MakeTriangles;

PROCEDURE Normalize(f: SPBezierPWFunction.T) =
  (*
    Scales "f" so that its norm becomes 1, according to
    the integral scalar product "<,>".
  *)
  CONST Order = 3; (* Integration precision. *)
  BEGIN
    WITH
      alpha = 1.0d0/sqrt(SPPWFunction.DotProduct(f, f, Order, Id, Id, NIL)),
      fd = f.data^
    DO
      FOR kf := 0 TO LAST(fd) DO
        fd[kf].scale(alpha)
      END
    END
  END Normalize;

PROCEDURE Orthize(f, g: SPBezierPWFunction.T) =
  (*
    Replaces "f" by "f + alpha*g" where "alpha" is such 
    that "f" becomes orthogonal to "g" according to
    the integral scalar product "<,>".  Assumes that "g"
    has unit norm, so that "alpha" is simply "-<f,g>".
    
    Assumes also that the triangles present in "f"
    are a superset of those present in "g", so that 
    the modification of "f" does not require the 
    creation of new triangles.
  *)
  CONST Order = 3; (* Integration precision. *)
  VAR kg: NAT := 0;
  BEGIN
    WITH
      alpha = - SPPWFunction.DotProduct(f, g, Order, Id, Id, NIL),
      fd = f.data^,
      gd = g.data^
    DO
      SortTriangles(fd);
      SortTriangles(gd);
      FOR kf := 0 TO LAST(fd) DO
        IF kg <= LAST(gd) AND fd[kf].face = gd[kg].face THEN
          fd[kf].add(alpha, gd[kg]);
          INC(kg);
        END;
      END;
    END;
  END Orthize;

PROCEDURE OrthoNormalizePair(f, g: SPBezierPWFunction.T) =
  (*
    Replaces "f" by "f + g" and "g" by "f - g".  Then normalizes "f"
    makes "g" orthogonal to it, and normalizes "g", 
    according to the integral scalar product "<,>".
    
    Assumes that the triangles present in "f" are the same as those
    present in "g", so that the modifications above do not require the
    creation of new triangles.
  *)
  BEGIN
    WITH
      fd = f.data^,
      gd = g.data^
    DO
      SortTriangles(fd);
      SortTriangles(gd);
      FOR k := 0 TO LAST(fd) DO
        <* ASSERT fd[k].face = gd[k].face *>
        fd[k].add(1.0d0, gd[k]);
        gd[k].add(-0.5d0, fd[k]);
      END;
      Normalize(f);
      Orthize(g, f);
      Normalize(g);
    END
  END OrthoNormalizePair;

BEGIN
  
END SPOrthoBasis.