MODULE Thicken;

IMPORT ParseParams;
IMPORT Cell, Face, Edge, Vertex, Model, LR3, Text;
  
TYPE 
  CT = Cell.CornerTags;
  
PROCEDURE AddSimpleLayer(
    READONLY core: Model.T;
    READONLY F: Face.List;
    READONLY R: Face.List;
    READONLY S: Face.List;  (* The set "S" of ``ramp'' faces. *)
    liftVertex: PROCEDURE (u: Vertex.Num): LR3.T;
    mat: TEXT;                  (* Properties of material *)
    vis: TEXT;                  (* Visual appearance of material *)
  ): Model.T =
  VAR nv: CARDINAL := NUMBER(core.vertex^);
  VAR nc: CARDINAL := NUMBER(core.cell^);
  BEGIN
    WITH
      nF = NUMBER(F),
      
      E = Edge.GetEdges(F)^,    (* Edges of faces in "F" *)
      nE = NUMBER(E),
      
      V = Edge.GetVertices(E)^, (* Vertices of faces in "F" *)
      nV = NUMBER(V),
      
      nP = 2 * nE - 3 * nF,     (* Number of perimeter edges in "E" *)
          
      nR = NUMBER(R),
      nS = NUMBER(S),

      rvNew = Vertex.Copy(core.vertex^, nV + nF + nP), 
      vNew = rvNew^,
      
      rcNew = Cell.Copy(core.cell^, 2*nF + 4*nE + nR + 4*nS),
      cNew = rcNew^,

      vLift = ComputeVLift(vNew, nv, V, liftVertex)^,
        (* 
          "vLift[v]" = number of the new roof vertex above vertex number "v",
          if the latter if a floor vertex; or "Vertex.None" otherwise. *)

      fMid = ComputeFMid(vNew, nv, vLift, core.cell^, F)^, 
        (* 
          "fMid[k]" = number of the new buried vertex above face "F[k]". *)
          
      eMid = ComputeEMid(vNew, nv, vLift, E, perimeterOnly := TRUE)^
        (*
          "eMid[k]" is the number of the new buried vertex above edge "E[k]",
          if the latter is a perimeter edge; or "Vertex.None" otherwise. *)
    DO

      (* Add the face-bound cells: *)
  
      FOR fk := 0 TO nF-1 DO
        WITH 
          face = F[fk],
          u = face.v[0], ut = face.t[0],
          v = face.v[1], vt = face.t[1],
          w = face.v[2], wt = face.t[2],
          m = fMid[fk],  mt = face.s[0],
          x = vLift[u],
          y = vLift[v],
          z = vLift[w],
          fName = MakeFaceName(face, cNew, vNew)
        DO
          (* Add the two face-associated cells: *)
          EVAL AddCell(cNew, nc,  m, w, v, u, CT{mt, wt, vt, ut}, fName & ".i");
          EVAL AddCell(cNew, nc,  m, x, y, z, CT{mt, ut, vt, wt}, fName & ".o");
        END
      END;

      (* Add the edge-bound cells: *)

      PROCEDURE GetFMid(fk: Face.Num; ek: Edge.Num): Vertex.Num =
        (*
          Returns the buried vertex "fMid[fk]" that lies above face "F[fk]".
          If that does nor exist, returns the vertex "eMid[ek]". *)
        BEGIN
          IF fk = Face.None THEN RETURN eMid[ek] ELSE RETURN fMid[fk] END
        END GetFMid;

      BEGIN
        FOR ek := 0 TO nE-1 DO 
          WITH 
            e = E[ek],
            
            u = e.v[0],
            uName = vNew[u].name,
            x = vLift[u],
            
            v = e.v[1],
            vName = vNew[v].name,
            y = vLift[v],
            
            eName = UndirEdgeName(uName, vName),
            uvName = DirEdgeName(uName, vName),
            vuName = DirEdgeName(vName, uName),
            
            vl = GetFMid(e.f[0], ek),
            vr = GetFMid(e.f[1], ek) 
          DO
            EVAL AddCell(cNew, nc,  vl, u, v, vr, CT{'l','u','v','r'}, eName  & ".i");
            EVAL AddCell(cNew, nc,  vl, v, y, vr, CT{'l','i','o','r'}, vuName & ".s");
            EVAL AddCell(cNew, nc,  vl, y, x, vr, CT{'l','v','u','r'}, uvName & ".s");
            EVAL AddCell(cNew, nc,  vl, x, u, vr, CT{'l','o','i','r'}, eName  & ".o")
          END
        END
      END; 
      
       (* Add the ramp cells: *)
      
      FOR rk := 0 TO nR-1 DO
        WITH
          face = PermuteRFace(R[rk], vLift),
          u = face.v[0],
          v = face.v[1],
          w = face.v[2],
          x = vLift[u],
          fName = MakeFaceName(face, cNew, vNew)
        DO
          EVAL AddCell(cNew, nc,  x, w, v, u,  CT{'o','w','v','i'}, fName & ".r")
        END
      END;
      
      FOR sk := 0 TO nS-1 DO
        WITH
          face = PermuteSFace(S[sk], vLift),
          u = face.v[0],
          v = face.v[1],
          w = face.v[2],
          x = vLift[u],
          y = vLift[v],
          ek = FindEdge(E, u, v),
          m = eMid[ek],
          fName = MakeFaceName(face, cNew, vNew)
        DO
          EVAL AddCell(cNew, nc,  m, v, u, w,  CT{'m','v','u','w'}, fName & ".b");
          EVAL AddCell(cNew, nc,  m, u, x, w,  CT{'m','i','o','w'}, fName & ".u");
          EVAL AddCell(cNew, nc,  m, x, y, w,  CT{'m','u','v','w'}, fName & ".r");
          EVAL AddCell(cNew, nc,  m, y, v, w,  CT{'m','o','i','w'}, fName & ".v");
        END
      END;
      
      <* ASSERT nv = NUMBER(vNew) *>
      <* ASSERT nc = NUMBER(cNew) *>
      FOR i := NUMBER(core.cell^) TO nc-1 DO 
        cNew[i].mat := mat; cNew[i].vis := vis
      END;
      RETURN Model.T{name := core.name, vertex := rvNew, cell := rcNew}
    END
  END AddSimpleLayer;
  
PROCEDURE AddTripleLayer(
    READONLY core: Model.T;
    READONLY F: Face.List;
    READONLY R: Face.List;
    READONLY S: Face.List;  (* The set "S" of ``ramp'' faces. *)
    liftVertex: PROCEDURE (u: Vertex.Num): LR3.T;
    liftEdge: PROCEDURE (u, v: Vertex.Num): LR3.T;
    liftFace: PROCEDURE (u, v, w: Vertex.Num): LR3.T;
    mat: TEXT;                  (* Properties of material *)
    vis: TEXT;                  (* Visual appearance of material *)
  ): Model.T =
  VAR nv: CARDINAL := NUMBER(core.vertex^);
  VAR nc: CARDINAL := NUMBER(core.cell^);
  BEGIN
    WITH
      nF = NUMBER(F),

      E = Edge.GetEdges(F)^,    (* Edges of faces in "F" *)
      nE = NUMBER(E),
      
      V = Edge.GetVertices(E)^, (* Vertices of faces in "F" *)
      nV = NUMBER(V),
      
      nP = 2 * nE - 3 * nF,         (* Number of perimeter edges in "E" *)

      nR = NUMBER(R),
      nS = NUMBER(S),

      rvNew = Vertex.Copy(core.vertex^, nV + nF + nP), 
      vNew = rvNew^,
      
      rcNew = Cell.Copy(core.cell^, nF + 3*nE + nR + 3*nS),
      cNew = rcNew^,

      vLift = ComputeVLift(vNew, nv, V, liftVertex)^,
        (* 
          "vLift[v]" = number of the new roof vertex above vertex number "v",
          if the latter if a floor vertex; or "Vertex.None" otherwise. *)

      fLift = ComputeFLift(vNew, nv, F, core.cell^, liftFace)^, 
        (* 
          "fLift[k]" = number of the new roof vertex above face "F[k]". *)

      eLift = ComputeELift(vNew, nv, E, liftEdge, perimeterOnly := TRUE)^
        (*
          "eLift[k]" is the number of the new roof vertex above edge "E[k]",
          if the latter is a perimeter edge; or "Vertex.None" otherwise. *)
    DO
      (* Add face-bound cells: *)

      FOR fk := 0 TO nF-1 DO
        WITH 
          face = F[fk],
          u = face.v[0], ut = face.t[0],
          v = face.v[1], vt = face.t[1], 
          w = face.v[2], wt = face.t[2],
          m = fLift[fk], mt = face.s[0],
          fName = MakeFaceName(face, cNew, vNew)
        DO
          EVAL AddCell(cNew, nc,  m, w, v, u, CT{mt, wt, vt, ut}, fName & ".c");
        END
      END;

      (* Add the edge-bound cells: *)

      PROCEDURE GetFLift(fk: Face.Num; ek: Edge.Num): Vertex.Num =
        (*
          Returns the roof vertex "fLift[fk]" that lies above face "F[fk]".
          If that does nor exist, returns the vertex "eLift[ek]". *)
        BEGIN
          IF fk = Face.None THEN RETURN eLift[ek] ELSE RETURN fLift[fk] END
        END GetFLift;

      BEGIN
        FOR ek := 0 TO nE-1 DO 
          WITH 
            e = E[ek],
            
            u = e.v[0],
            uName = vNew[u].name,
            x = vLift[u],
            
            v = e.v[1],
            vName = vNew[v].name,
            y = vLift[v],
            
            eName = UndirEdgeName(uName, vName),
            uvName = DirEdgeName(uName, vName),
            vuName = DirEdgeName(vName, uName),
            
            vl = GetFLift(e.f[0], ek),
            vr = GetFLift(e.f[1], ek)
          DO
            EVAL AddCell(cNew, nc,  vl, u, v, vr, CT{'l','u','v','r'}, eName  & ".c");
            EVAL AddCell(cNew, nc,  vl, x, u, vr, CT{'l','o','i','r'}, uvName & ".a");
            EVAL AddCell(cNew, nc,  vl, v, y, vr, CT{'l','i','o','r'}, vuName & ".a");
          END
        END;
      END;

      (* Add the ramp cells: *)
      
      FOR rk := 0 TO nR-1 DO
        WITH
          face = PermuteRFace(R[rk], vLift),
          u = face.v[0],
          v = face.v[1],
          w = face.v[2],
          x = vLift[u],
          fName = MakeFaceName(face, cNew, vNew)
        DO
          EVAL AddCell(cNew, nc,  x, w, v, u,  CT{'o','w','v','i'}, fName & ".r")
        END
      END;
      
      FOR sk := 0 TO nS-1 DO
        WITH
          face = PermuteSFace(S[sk], vLift),
          u = face.v[0],
          v = face.v[1],
          w = face.v[2],
          x = vLift[u],
          y = vLift[v],
          ek = FindEdge(E, u, v),
          m = eLift[ek],
          fName = MakeFaceName(face, cNew, vNew)
        DO
          EVAL AddCell(cNew, nc,  m, v, u, w,  CT{'m','v','u','w'}, fName & ".b");
          EVAL AddCell(cNew, nc,  m, u, x, w,  CT{'m','i','o','w'}, fName & ".u");
          EVAL AddCell(cNew, nc,  m, y, v, w,  CT{'m','o','i','w'}, fName & ".v");
        END
      END;
      
      <* ASSERT nv = NUMBER(vNew) *>
      <* ASSERT nc = NUMBER(cNew) *>
      FOR i := NUMBER(core.cell^) TO nc-1 DO 
        cNew[i].mat := mat; cNew[i].vis := vis
      END;
      RETURN Model.T{name := core.name, vertex := rvNew, cell := rcNew}
    END
  END AddTripleLayer;

PROCEDURE AddQuadLayer(
    READONLY core: Model.T;
    READONLY F: Face.List;
    READONLY R: Face.List;
    READONLY S: Face.List;  (* The set "S" of ``ramp'' faces. *)
    liftVertex: PROCEDURE (u: Vertex.Num): LR3.T;
    liftEdge: PROCEDURE (u, v: Vertex.Num): LR3.T;
    mat: TEXT;                  (* Properties of material *)
    vis: TEXT;                  (* Visual appearance of material *)
  ): Model.T =
  VAR nv: CARDINAL := NUMBER(core.vertex^);
  VAR nc: CARDINAL := NUMBER(core.cell^);
  BEGIN
    WITH
      nF = NUMBER(F),

      E = Edge.GetEdges(F)^,    (* Edges of faces in "F" *)
      nE = NUMBER(E),
      
      V = Edge.GetVertices(E)^, (* Vertices of faces in "F" *)
      nV = NUMBER(V),
      
      nP = 2 * nE - 3 * nF,     (* Number of unshared edges in "E" *)

      nR = NUMBER(R),
      nS = NUMBER(S),

      rvNew = Vertex.Copy(core.vertex^, 2 * nV + nF + nE + nP), 
      vNew = rvNew^,
      
      rcNew = Cell.Copy(core.cell^, 8*nF + 5*nE + 2*nR + 7*nS),
      cNew = rcNew^,

      vLift = ComputeVLift(vNew, nv, V, liftVertex)^,
        (* 
          "vLift[v]" = number of the new roof vertex above vertex number "v",
          if the latter if a floor vertex; or "Vertex.None" otherwise. *)

      vMid = ComputeVMid(vNew, nv, vLift, V)^, 
        (* 
          "vMid[v]" = number of the new buried vertex above vertex number "v",
          if the latter if a floor vertex; or "Vertex.None" otherwise. *)
          
      fMid = ComputeFMid(vNew, nv, vLift, core.cell^, F)^, 
        (* 
          "fMid[k]" = number of the new buried vertex above face "F[k]". *)
          
      eLift = ComputeELift(vNew, nv, E, liftEdge, perimeterOnly := FALSE)^,
        (*
          "eLift[k]" is the number of the new roof vertex above edge "E[k]". *)

      eMid = ComputeEMid(vNew, nv, vLift, E, perimeterOnly := TRUE)^
        (*
          "eMid[k]" is the number of the new buried vertex above edge "E[k]",
          if the latter is a perimeter edge; or "Vertex.None" otherwise. *)
    DO

      (* Add new face-bound cells: *)

      FOR fk := 0 TO nF-1 DO
        WITH 
          face = F[fk],
          
          u = face.v[0], ut = face.t[0],
          uName = vNew[u].name,
          uc = vLift[u],
          um = vMid[u],
          
          v = face.v[1], vt = face.t[1],
          vName = vNew[v].name,
          vc = vLift[v],
          vm = vMid[v],
          
          w = face.v[2], wt = face.t[2],
          wName = vNew[w].name,
          wc = vLift[w],
          wm = vMid[w],
          
          uvc = eLift[FindEdge(E, u, v)],
          vwc = eLift[FindEdge(E, v, w)],
          wuc = eLift[FindEdge(E, w, u)],
          
          md = fMid[fk], mt = face.s[0],
          fName = MakeFaceName(face, cNew, vNew)
        DO
          (* Up-pointing, buried: *)
          EVAL AddCell(cNew, nc,  md,  w,   w,   u,   CT{mt, wt, vt, ut}, fName & ".c");

          (* Down-pointing, exposed: *)
          EVAL AddCell(cNew, nc,  md,  vwc, wuc, uvc, CT{mt, ut, vt, wt}, fName & ".m");

          (* Coner, exposed: *)
          EVAL AddCell(cNew, nc,  um, uc,  uvc, wuc, CT{mt, ut, vt, wt}, uName & ".o");
          EVAL AddCell(cNew, nc,  vm, uvc, vc,  vwc, CT{mt, ut, vt, wt}, vName & ".o");
          EVAL AddCell(cNew, nc,  wm, wuc, vwc, wc,  CT{mt, ut, vt, wt}, wName & ".o");

          (* Corner, buried: *)
          EVAL AddCell(cNew, nc,  md, wuc, uvc, um,  CT{mt, wt, vt, ut}, uName & ".o");
          EVAL AddCell(cNew, nc,  md, vwc, vm,  uvc, CT{mt, wt, vt, ut}, vName & ".o");
          EVAL AddCell(cNew, nc,  md, wm,  vwc, wuc, CT{mt, wt, vt, ut}, wName & ".o");
        END
      END;

      (* Add the edge-bound cells: *)

      PROCEDURE GetFMid(fk: Face.Num; ek: Edge.Num): Vertex.Num =
        (*
          Returns the buried vertex "fMid[fk]" that lies above face "F[fk]".
          If that does not exist, returns instead the buried vertex "eMid[ek]". *)
        BEGIN
          IF fk = Face.None THEN RETURN eMid[ek] ELSE RETURN fMid[fk] END
        END GetFMid;

      BEGIN
        FOR ek := 0 TO nE-1 DO 
          WITH 
            e = E[ek],
            
            u = e.v[0],
            uName = vNew[u].name,
            um = vMid[u],
            
            v = e.v[1],
            vName = vNew[v].name,
            vm = vMid[v],
            
            em = eMid[ek],
            
            eName = UndirEdgeName(uName, vName),
            uvName = DirEdgeName(uName, vName),
            vuName = DirEdgeName(vName, uName),
            
            vl = GetFMid(e.f[0], ek),
            vr = GetFMid(e.f[1], ek)
          DO
            EVAL AddCell(cNew, nc, vl, um, u,  vr, CT{'l','o','i','r'}, uvName & ".r");
            EVAL AddCell(cNew, nc, vl, u,  v,  vr, CT{'l','u','v','r'}, eName  & ".c");
            EVAL AddCell(cNew, nc, vl, v,  vm, vr, CT{'l','i','o','r'}, vuName & ".r");
            EVAL AddCell(cNew, nc, vl, vm, em, vr, CT{'l','i','o','r'}, vuName & ".s");
            EVAL AddCell(cNew, nc, vl, em, um, vr, CT{'l','o','i','r'}, uvName & ".s");
          END
        END;
      END;

      (* Add the ramp cells: *)
      
      FOR rk := 0 TO nR-1 DO
        WITH
          face = PermuteRFace(R[rk], vLift),
          u = face.v[0],
          v = face.v[1],
          w = face.v[2],
          m = vMid[u],
          x = vLift[u],
          fName = MakeFaceName(face, cNew, vNew)
        DO
          EVAL AddCell(cNew, nc,  m, w, v, u,  CT{'m','w','v','u'}, fName & ".b");
          EVAL AddCell(cNew, nc,  m, v, w, x,  CT{'m','v','w','u'}, fName & ".r");
        END
      END;
      
      FOR sk := 0 TO nS-1 DO
        WITH
          face = PermuteSFace(S[sk], vLift),
          u = face.v[0],
          v = face.v[1],
          w = face.v[2],
          r = vMid[u],
          s = vMid[v],
          x = vLift[u],
          y = vLift[v],
          ek = FindEdge(E, u, v),
          m = eMid[ek],
          t = eLift[ek],
          fName = MakeFaceName(face, cNew, vNew)
        DO
          EVAL AddCell(cNew, nc,  m, v, u, w,  CT{'m','v','u','w'}, fName & ".b");
          
          EVAL AddCell(cNew, nc,  m, u, r, w,  CT{'m','i','o','w'}, fName & ".au");
          EVAL AddCell(cNew, nc,  m, r, x, w,  CT{'m','i','o','w'}, fName & ".bu");
          EVAL AddCell(cNew, nc,  r, x, t, w,  CT{'i','o','t','w'}, fName & ".ur");
          
          EVAL AddCell(cNew, nc,  m, s, v, w,  CT{'m','o','i','w'}, fName & ".av");
          EVAL AddCell(cNew, nc,  m, y, s, w,  CT{'m','o','i','w'}, fName & ".bv");
          EVAL AddCell(cNew, nc,  s, t, y, w,  CT{'i','t','o','w'}, fName & ".vr");
        END
      END;
      
      <* ASSERT nv = NUMBER(vNew) *>
      <* ASSERT nc = NUMBER(cNew) *>
      FOR i := NUMBER(core.cell^) TO nc-1 DO 
        cNew[i].mat := mat; cNew[i].vis := vis
      END;
      RETURN Model.T{name := core.name, vertex := rvNew, cell := rcNew};
    END
  END AddQuadLayer;

PROCEDURE ComputeVLift(
    VAR vertex: Vertex.List;
    VAR nv: CARDINAL;
    READONLY vFloor: ARRAY OF Vertex.Num;
    liftVertex: PROCEDURE (u: Vertex.Num): LR3.T;
  ): REF ARRAY OF Vertex.Num =
  BEGIN
    WITH
      r = NEW(REF ARRAY OF Vertex.Num, NUMBER(vertex)), vLift = r^
    DO
      FOR v := 0 TO LAST(vLift) DO vLift[v] := Vertex.None END;
      FOR k := 0 TO LAST(vFloor) DO 
        WITH 
          v = vFloor[k],
          vName = vertex[v].name,
          c = liftVertex(v)
        DO
          vLift[v] := AddVertex(vertex, nv, c, vName & ".v")
        END
      END;
      RETURN r
    END
  END ComputeVLift;

PROCEDURE ComputeFLift(
    VAR vertex: Vertex.List;
    VAR nv: CARDINAL;
    READONLY fFloor: Face.List;
    READONLY cell: Cell.List;
    liftFace: PROCEDURE (u, v, w: Vertex.Num): LR3.T;
  ): REF ARRAY OF Vertex.Num =
  BEGIN
    WITH
      r = NEW(REF ARRAY OF Vertex.Num, NUMBER(fFloor)), fLift = r^
    DO
      FOR k := 0 TO LAST(fFloor) DO 
        WITH 
          face = fFloor[k],
          u = face.v[0],
          v = face.v[1],
          w = face.v[2],
          fName = MakeFaceName(face, cell, vertex)
        DO
          fLift[k] := AddVertex(vertex, nv, liftFace(u, v, w), fName & ".f");
        END
      END;
      RETURN r
    END
  END ComputeFLift;

PROCEDURE ComputeELift(
    VAR vertex: Vertex.List;
    VAR nv: CARDINAL;
    READONLY eFloor: Edge.List;
    liftEdge: PROCEDURE (u, v: Vertex.Num): LR3.T;
    perimeterOnly: BOOLEAN;
  ): REF ARRAY OF Vertex.Num =
  BEGIN
    WITH
      r = NEW(REF ARRAY OF Vertex.Num, NUMBER(eFloor)), eLift = r^
    DO
      FOR k := 0 TO LAST(eFloor) DO 
        WITH 
          e = eFloor[k]
        DO
          IF perimeterOnly AND (e.f[0] # Face.None AND e.f[1] # Face.None)
          THEN
            eLift[k] := Vertex.None
          ELSE
            WITH
              u = e.v[0],
              uName = vertex[u].name,
              v = e.v[1],
              vName = vertex[v].name,
              eName = UndirEdgeName(uName, vName)
            DO
              eLift[k] := AddVertex(vertex, nv, liftEdge(u, v), eName & ".e")
            END
          END
        END
      END;
      RETURN r
    END
  END ComputeELift;

PROCEDURE ComputeFMid(
    VAR vertex: Vertex.List;
    VAR nv: CARDINAL;
    READONLY vLift: ARRAY OF Vertex.Num;
    READONLY cell: Cell.List;
    READONLY fFloor: Face.List;
  ): REF ARRAY OF Vertex.Num = 
  BEGIN
    WITH
      r = NEW(REF ARRAY OF Vertex.Num, NUMBER(fFloor)),
      fMid = r^
    DO
      FOR i := 0 TO LAST(fFloor) DO
        WITH 
          face = fFloor[i],
          u = face.v[0],
          v = face.v[1],
          w = face.v[2],
          uc = vLift[u],
          vc = vLift[v],
          wc = vLift[w],
          fName = MakeFaceName(face, cell, vertex)
        DO
          (* Create the new buried vertex: *)
          WITH
            c = Avg2(
                Avg3(vertex[u].c, vertex[v].c, vertex[w].c),
                Avg3(vertex[uc].c,  vertex[vc].c,  vertex[wc].c )
              ),
            name = fName & ".f"
          DO
            fMid[i] := AddVertex(vertex, nv, c, name)
          END;
        END
      END;
      RETURN r
    END;
  END ComputeFMid;

PROCEDURE ComputeEMid(
    VAR vertex: Vertex.List;
    VAR nv: CARDINAL;
    READONLY vLift: ARRAY OF Vertex.Num;
    READONLY eFloor: Edge.List;
    perimeterOnly: BOOLEAN;
  ): REF ARRAY OF Vertex.Num =
  BEGIN
    WITH
      r = NEW(REF ARRAY OF Vertex.Num, NUMBER(eFloor)),
      eMid = r^
    DO
      FOR k := 0 TO LAST(eFloor) DO 
        WITH e = eFloor[k] DO 
          IF perimeterOnly AND (e.f[0] # Face.None AND e.f[1] # Face.None) THEN
            eMid[k] := Vertex.None
          ELSE
            WITH 
              o = e.v[0],
              oName = vertex[o].name,
              oc = vLift[o],

              d = e.v[1],
              dName = vertex[d].name,
              dc = vLift[d],

              eName = UndirEdgeName(oName, dName),

              c = Avg2(
                Avg2(vertex[o].c,  vertex[d].c), 
                Avg2(vertex[oc].c, vertex[dc].c)
              )
            DO
              eMid[k] := AddVertex(vertex, nv, c, eName & ".m")
            END
          END
        END
      END;
      RETURN r
    END
  END ComputeEMid;

PROCEDURE ComputeVMid(
    VAR vertex: Vertex.List;
    VAR nv: CARDINAL;
    READONLY vLift: ARRAY OF Vertex.Num;
    READONLY vFloor: ARRAY OF Vertex.Num;
  ): REF ARRAY OF Vertex.Num =
  BEGIN
    WITH
      r = NEW(REF ARRAY OF Vertex.Num, NUMBER(vertex)), vMid = r^
    DO
      FOR v := 0 TO LAST(vMid) DO vMid[v] := Vertex.None END;
      FOR k := 0 TO LAST(vFloor) DO 
        WITH 
          v = vFloor[k],
          vName = vertex[v].name,
          
          u = vLift[v],
          c = Avg2(vertex[v].c,  vertex[u].c)
        DO
          vMid[v] := AddVertex(vertex, nv, c, vName & "m")
        END
      END;
      RETURN r
    END
  END ComputeVMid;
  
PROCEDURE MakeFaceName(
    READONLY f: Face.T; 
    READONLY cell: Cell.List;
    READONLY vertex: Vertex.List;
  ): Face.Name =
  BEGIN
    IF f.c[0] = Cell.None THEN
      VAR u := ARRAY [0..2] OF TEXT {
        vertex[f.v[0]].name, vertex[f.v[1]].name, vertex[f.v[2]].name
      };
      BEGIN
        FOR i := 0 TO 1 DO
          IF Text.Compare(u[0], u[1]) > 0 OR Text.Compare(u[0], u[2]) > 0 THEN
            VAR t := u[0]; BEGIN u[0] := u[1]; u[1] := u[2]; u[2] := t END
          END
        END;
        RETURN Model.Merge3Names(u[0], u[1], u[2], "{", "}")
      END
    ELSE
      RETURN cell[f.c[0]].name & ":" & Text.FromChar(f.s[0])
    END
  END MakeFaceName;

PROCEDURE AddVertex(
    VAR vertex: Vertex.List;
    VAR nv: CARDINAL;
    READONLY c: LR3.T;
    name: TEXT;
  ): Vertex.Num =
  (*
    Adds a new vertex with given coordinates. *)
  BEGIN
    vertex[nv] := Vertex.T{c := c, name := name, fixed := FALSE};
    INC(nv);
    RETURN nv-1
  END AddVertex;
  
PROCEDURE AddCell(
    VAR cell: Cell.List;
    VAR nc: CARDINAL;
    p0, p1, p2, p3: Vertex.Num;
    tags: Cell.CornerTags; 
    name: TEXT;
  ): Cell.Num =
  BEGIN
    <* ASSERT p0 # Vertex.None *>
    <* ASSERT p1 # Vertex.None *>
    <* ASSERT p2 # Vertex.None *>
    <* ASSERT p3 # Vertex.None *>
    cell[nc] := Cell.Make(p0, p1, p2, p3, tags, name);
    INC(nc);
    RETURN nc-1
  END AddCell;
  
PROCEDURE PermuteRFace(
    READONLY face: Face.T; 
    READONLY vLift: ARRAY OF Vertex.Num
  ): Face.T =
  (*
    Permutes the vertices of "face" until corner "v[0]" is the one
    on the floor perimeter. *)
  VAR f := face;
  BEGIN
    WHILE vLift[f.v[0]] = Vertex.None DO
      <* ASSERT vLift[f.v[1]] # Vertex.None OR vLift[f.v[2]] # Vertex.None *>
      RotateFace(f)
    END;
    <* ASSERT vLift[f.v[1]] = Vertex.None AND vLift[f.v[2]] = Vertex.None *>
    RETURN f
  END PermuteRFace;

PROCEDURE PermuteSFace(
    READONLY face: Face.T; 
    READONLY vLift: ARRAY OF Vertex.Num
  ): Face.T =
  (*
    Permutes the vertices of face "f" until side "(v[0],v[1])" is the one
    on the floor perimeter. *)
  VAR f := face;
  BEGIN
    WHILE vLift[f.v[2]] # Vertex.None DO
      <* ASSERT vLift[f.v[0]] = Vertex.None OR vLift[f.v[1]] = Vertex.None *>
      RotateFace(f)
    END;
    <* ASSERT vLift[f.v[0]] # Vertex.None AND vLift[f.v[1]] # Vertex.None *>
    RETURN f
  END PermuteSFace;
  
PROCEDURE FindEdge(READONLY E: Edge.List; x, y: Vertex.Num): Edge.Num =
  BEGIN
    FOR ek := 0 TO LAST(E) DO 
      WITH e = E[ek], u = e.v[0], v = e.v[1] DO
        IF (x = u AND y = v) OR (x = v AND y = u) THEN RETURN ek END
      END
    END;
    RETURN Edge.None
  END FindEdge;

PROCEDURE RotateFace(VAR f: Face.T) =
  (*
    Rotates the corners of a face so that new "u" is old "v", etc. *)
  BEGIN
    WITH
      v = f.v, t = f.t
    DO
      VAR z := v[0]; BEGIN v[0] := v[1]; v[1] := v[2]; v[2] := z END;
      VAR z := t[0]; BEGIN t[0] := t[1]; t[1] := t[2]; t[2] := z END;
    END
  END RotateFace;

PROCEDURE DirEdgeName(u, v: TEXT): TEXT =
  BEGIN
    RETURN Model.Merge2Names(u, v, "(", ")")
  END DirEdgeName;
  
PROCEDURE UndirEdgeName(u, v: TEXT): TEXT =
  BEGIN
    IF Text.Compare(u, v) > 0 THEN
      VAR t := u; BEGIN u := v; v := t END
    END;
    RETURN Model.Merge2Names(u, v, "{", "}")
  END UndirEdgeName;
  
PROCEDURE ParseLayerOptions(pp: ParseParams.T): REF ARRAY OF LayerSpecs
  RAISES{ParseParams.Error} =
  CONST MaxL = 100;
  VAR nl: CARDINAL := 0;
  BEGIN
    WITH 
      rtl = NEW(REF ARRAY OF LayerSpecs, MaxL), layer = rtl^
    DO
      LOOP
        WITH a = layer[nl] DO
          IF nl = 0 THEN
            pp.getKeyword("-kernel");
            a.kind := LayerKind.Kernel;
          ELSE
            IF NOT pp.keywordPresent("-layer") THEN EXIT END;
            IF pp.testNext("simple") THEN
              a.kind := LayerKind.Simple
            ELSIF pp.testNext("triple") THEN
              a.kind := LayerKind.Triple
            ELSIF pp.testNext("quad") THEN
              a.kind := LayerKind.Quad
            ELSE
              pp.error("layer type missing")
            END
          END;
          IF NOT pp.testNext("thickness") THEN  
            pp.error("\"thickness\" expected")
          END;
          a.thickness := pp.getNextLongReal(0.0d0);
          VAR mat, vis, done: BOOLEAN := FALSE;
          BEGIN
            a.mat := "default"; a.vis := "default";
            WHILE NOT done DO
              IF NOT mat AND pp.testNext("material") THEN  
                a.mat := pp.getNext(); mat := TRUE
              ELSIF NOT vis AND pp.testNext("visual") THEN  
                a.vis := pp.getNext(); vis := TRUE
              ELSE
                done := TRUE
              END
            END
          END;
          INC(nl)
        END;
      END (* LOOP *);
      IF nl = 0 THEN
        pp.error("must specify at least one \"-layer\"")
      END;
      WITH
        rl = NEW(REF ARRAY OF LayerSpecs, nl)
      DO
        rl^ := SUBARRAY(rtl^, 0, nl);
        RETURN rl
      END
    END;
  END ParseLayerOptions;

PROCEDURE Avg2(READONLY x, y: LR3.T): LR3.T =
  BEGIN
    RETURN LR3.Mix(0.5d0, x, 0.5d0, y)
  END Avg2;
  
PROCEDURE Avg3(READONLY x, y, z: LR3.T): LR3.T =
  BEGIN
    RETURN LR3.T{
      (x[0] + y[0] + z[0])/3.0d0,
      (x[1] + y[1] + z[1])/3.0d0,
      (x[2] + y[2] + z[2])/3.0d0
    }
  END Avg3;

BEGIN
END Thicken.