MODULE Squared;

(*
   This module contain procedures to build several structures such as: n-gons,
   triangles and squares, and solid polyhedra: cube, ball, bigcube (3D array of
   cube) with procedures for the gluing two such bigcubes.

   Revisions:

   30-08-2000 : Nice version of the "MakeOctahedron" procedure.
   19-09-2000 : Added the procedure "MakeDodecahedronTriang".
   27-10-2000 : Modified "SetCubePropiertes" procedure.
   04-11-2000 : Added the "CubeNegVertices" and "CubeBarycenter" procedures.
*)

IMPORT Octf, Triangulation, Tools, LR4;

FROM Octf IMPORT Enext_1, Clock, SetEnext, SetFace, Fnext, Enext,
          SetFnext, SetEdge, SetEdgeAll, Fnext_1, Spin, Srot;

FROM Triangulation IMPORT Pair, Org, MakeVertex, SetOrg, MakeFacetEdge, Pneg,
          SetAllOrgs, OrgV, SetAllPneg, MakePolyhedron, Node, SetPneg, Vertex,
          SetNextPneg, Glue, DegreeOfVertex, Ppos, MakeTetraTopo, Topology,
          Coords;

TYPE
  PAIRS = ARRAY[0..3] OF Pair;
 
VAR
  FacetEdgeCount: CARDINAL := 0;
  PolyhedronCount: CARDINAL := 0;

PROCEDURE MakeTriangle() : Pair =
  (* Builds a triangular face and set the three pairs facetedges with the
     same face component. *)
  BEGIN
    WITH
      a = MakeFacetEdge(),
      b = MakeFacetEdge(),
      c = MakeFacetEdge(),
      f = a.facetedge.face,
      u = MakeVertex(),
      v = MakeVertex(),
      w = MakeVertex()
    DO 
      a.facetedge.num := FacetEdgeCount; INC(FacetEdgeCount);         
      SetOrg(a, u); SetOrg(Clock(a),v);

      b.facetedge.num := FacetEdgeCount; INC(FacetEdgeCount);
      SetEnext(a,b); SetFace(b,f);
      SetOrg(b,v); SetOrg(Clock(b),w);

      c.facetedge.num := FacetEdgeCount; INC(FacetEdgeCount);
      SetEnext(b,c); SetFace(c,f);
      SetOrg(c, w);  SetOrg(Clock(c), Org(a));
      RETURN a;
    END
  END MakeTriangle;

PROCEDURE MakeGon(n: CARDINAL) : Pair =
  (* Builds one n-gon face. The "n-pairs" have the same face component.*)
  VAR
    sp : REF ARRAY OF Pair;
    sv : REF ARRAY OF Vertex;
  BEGIN
    sp := NEW(REF ARRAY OF Pair,n);
    sv := NEW(REF ARRAY OF Vertex,n);
    FOR i := 0 TO n-1 DO
      sp[i] := MakeFacetEdge();
      sv[i] := MakeVertex();
    END;

    WITH f = sp[0].facetedge.face DO
      FOR i := 0 TO n-1 DO  
        sp[i].facetedge.num := FacetEdgeCount;
        INC(FacetEdgeCount);         
        SetOrg(sp[i],sv[i]);
        SetOrg(Clock(sp[i]),sv[(i+1) MOD n]);
        SetEnext(sp[i],sp[(i+1) MOD n]);
        SetFace(sp[i],f);
      END
    END;   	
    RETURN sp[0];   
  END MakeGon;

PROCEDURE MakeGonFull(n: CARDINAL) : REF ARRAY OF Pair =
  (* Builds one n-gon face. The "n-pairs" have the same face component.*)
  VAR
    sp : REF ARRAY OF Pair;
    sv : REF ARRAY OF Vertex;
  BEGIN
    sp := NEW(REF ARRAY OF Pair,n);
    sv := NEW(REF ARRAY OF Vertex,n);
    FOR i := 0 TO n-1 DO
      sp[i] := MakeFacetEdge();
      sv[i] := MakeVertex();
    END;

    WITH f = sp[0].facetedge.face DO
      FOR i := 0 TO n-1 DO  
        sp[i].facetedge.num := FacetEdgeCount;
        INC(FacetEdgeCount);         
        SetOrg(sp[i],sv[i]);
        SetOrg(Clock(sp[i]),sv[(i+1) MOD n]);
        SetEnext(sp[i],sp[(i+1) MOD n]);
        SetFace(sp[i],f);
      END
    END;   	
    RETURN sp;   
  END MakeGonFull;
  
PROCEDURE MakeSquare() : Pair =
  (* Builds one squared face. The four "pairs" have the same face component.*)
  BEGIN
    WITH
      a = MakeFacetEdge(),
      b = MakeFacetEdge(),
      c = MakeFacetEdge(),
      d = MakeFacetEdge(),
      f = a.facetedge.face,
      u = MakeVertex(),
      v = MakeVertex(),
      w = MakeVertex(),
      x = MakeVertex()
    DO 
      a.facetedge.num := FacetEdgeCount; INC(FacetEdgeCount);         
      SetOrg(a, u); SetOrg(Clock(a),v);

      b.facetedge.num := FacetEdgeCount; INC(FacetEdgeCount);
      SetEnext(a,b); SetFace(b,f);
      SetOrg(b,v); SetOrg(Clock(b),w);

      c.facetedge.num := FacetEdgeCount; INC(FacetEdgeCount);
      SetEnext(b,c); SetFace(c,f);
      SetOrg(c, w);  SetOrg(Clock(c), x);

      d.facetedge.num := FacetEdgeCount; INC(FacetEdgeCount);
      SetEnext(c,d); SetFace(d,f);
      SetOrg(d, x);  SetOrg(Clock(d), Org(a));
      RETURN a;
    END
  END MakeSquare;

PROCEDURE MakeHalfOct( ) : ARRAY[0..3] OF Pair =
  (* Builds half octahedron without polyhedron information. *)
  VAR p : ARRAY[0..3] OF Pair;
  BEGIN
    FOR i := 0 TO  3 DO
      p[i] := MakeTriangle();
    END;

    FOR i := 0 TO 3 DO
      WITH
        a = Enext(p[i]),
        b = Clock(Enext_1(p[(i+1) MOD 4]))
      DO
        SetFnext(a,b);
        SetEdgeAll(a,a.facetedge.edge);
      END
    END;

    FOR i := 0 TO 3 DO
      SetAllOrgs(Enext_1(p[i]), Org(Enext_1(p[i])));
    END;

    RETURN p;
  END MakeHalfOct;

PROCEDURE MakeOctahedron( ) : ARRAY[0..7] OF Pair =
  (* Builds a topological octahedron, and return one pair facetedge "p[i]" for
     each face of the octahedron.
  *)
  VAR
    p   : ARRAY [0..7] OF Pair;
    s,i : ARRAY [0..3] OF Pair;
  BEGIN
    s := MakeHalfOct();
    i := MakeHalfOct();

    (* builds an octahedron with two half-octahedra "s" (superior) and
       "i" (inferior).
    *)
    FOR j := 0 TO 3 DO
      p[j] := s[j]; 
    END;

    FOR j := 4 TO 7 DO
      p[j] := i[j-4]; 
    END;  

    (* Set the common elements *)
    FOR j := 0 TO 3 DO
      SetFnext(p[j],p[j+4]);
      SetEdgeAll(p[j],p[j].facetedge.edge);
      SetAllOrgs(p[j], Org(p[j]));
      SetAllOrgs(Enext(p[j]), Org(Enext(p[j])));
    END;

    (* set the polyhedron attribute, 24 facet-edges have the
       same polyhedron.
    *)
    WITH
      q = MakePolyhedron()
    DO
      q.num := PolyhedronCount; INC(PolyhedronCount);
      FOR i := 0 TO 7 DO
        SetNextPneg(p[i],q);
      END
    END;

    (* this test is apply for assert that any pair facetedge "p[i]"
       really underly on a topological octahedron.*)

    FOR i := 0 TO 7 DO
      WITH
        d = DegreeOfVertex(Srot(p[i]))
      DO
        <* ASSERT d = 8 *>
      END;
    END; 

    RETURN p;
  END MakeOctahedron;

PROCEDURE SetRingFaceExis(a: Pair; faceExists: BOOLEAN) =
(* Set the face as existing or not. *)

  PROCEDURE SetFace(b: Pair) =
    BEGIN
      WITH
        t = NARROW(b.facetedge.face, Triangulation.Face)
      DO
        t.exists := faceExists;
      END
    END SetFace;	

  VAR an: Pair := a;
  BEGIN
    REPEAT
      SetFace(an);
      an := Fnext(an)
    UNTIL an = a 
  END SetRingFaceExis;

PROCEDURE SetEdgeExis(a: Pair; edgeExists: BOOLEAN) =
  (* Set the edge as existing or not . *)
  BEGIN
    WITH
      e = NARROW(a.facetedge.edge, Triangulation.Edge)
    DO
      e.exists := edgeExists;
    END
  END SetEdgeExis;

PROCEDURE SetVertexExis(a: Pair) =
  (* Set the vertex as not existing. *)
  BEGIN
    WITH
      v = NARROW(Org(a), Vertex)
    DO
      v.exists := FALSE;
    END
  END SetVertexExis;   

PROCEDURE MakeOctahedronTriang(Original: BOOLEAN) : ARRAY[0..7] OF Pair =
  (* Builds a triangulated octahedron, with eight tetrahedra, and return
     one facet-edge  pair by original face. If Original=TRUE the procedure
     emphasize the original elements.
  *)
  VAR
    a : ARRAY[0..7] OF ARRAY[0..7] OF Pair;
    b : ARRAY[0..7] OF Pair;
  BEGIN
    FOR i := 0 TO 7 DO
      a[i] := Triangulation.MakeTetraTopo(1,1);
    END;

    (* first level gluing tetrahedra *)
    EVAL Glue(Spin(a[1][1]),a[0][0],1);
    EVAL Glue(Spin(a[2][1]),a[1][0],1);
    EVAL Glue(Spin(a[3][1]),a[2][0],1);
    EVAL Glue(Spin(a[0][1]),a[3][0],1);

    (* second level gluing tetrahedra *)
    EVAL Glue(Spin(a[0][3]),a[4][2],1);
    EVAL Glue(Spin(a[1][3]),a[5][2],1);
    EVAL Glue(Spin(a[2][3]),a[6][2],1);
    EVAL Glue(Spin(a[3][3]),a[7][2],1);

    (* gluing between levels *)
    EVAL Glue(Spin(a[5][1]),a[4][0],1);
    EVAL Glue(Spin(a[6][1]),a[5][0],1);
    EVAL Glue(Spin(a[7][1]),a[6][0],1);
    EVAL Glue(Spin(a[4][1]),a[7][0],1);

    IF Original THEN
      (* emphasize the original elements *)
    
      SetRingFaceExis(a[0][1],FALSE);
      SetRingFaceExis(a[4][1],FALSE);
      SetEdgeExis(a[0][1],FALSE); 
      SetEdgeExis(a[4][1],FALSE);
      FOR i := 0 TO 3 DO
	SetRingFaceExis(Enext(a[i][1]),FALSE);
      END;
      FOR i := 4 TO 7 DO
	SetEdgeExis(Enext_1(a[i][1]),FALSE);
      END;
      SetVertexExis(a[4][1]);
    END;

    b[0] := a[0][2];   
    b[1] := a[1][2];
    b[2] := a[2][2];
    b[3] := a[3][2];
    b[4] := a[4][3];
    b[5] := a[5][3];
    b[6] := a[6][3];
    b[7] := a[7][3];

    RETURN b;

  END MakeOctahedronTriang;

<* UNUSED *>
PROCEDURE OldMakeIcosahedronTriang( ) : ARRAY[0..19] OF Pair =
  (* Prior version of the "MakeIcosahedronTriang": Builds a triangualted icosahedron,
     with twenty tetrahedra, and return one facet-edge pair by original face (i.e. 20).
  *)
  VAR
     a : ARRAY[0..19] OF ARRAY[0..7] OF Pair;
     b : ARRAY[0..19] OF Pair; (* this variable will be retuned *)
  BEGIN
    FOR i := 0 TO 19 DO
      a[i] := MakeTetraTopo(1,1);
    END;

    (* first level *)
    EVAL Glue(Spin(a[1][1]),a[0][0],1);  
    EVAL Glue(Spin(a[2][1]),a[1][0],1); 
    EVAL Glue(Spin(a[3][1]),a[2][0],1);
    EVAL Glue(Spin(a[4][1]),a[3][0],1);
    EVAL Glue(Spin(a[0][1]),a[4][0],1);

    (* next level *)
    EVAL Glue(Spin(a[0][3]),a[5][2],1);
    EVAL Glue(Spin(a[1][3]),a[6][2],1);
    EVAL Glue(Spin(a[2][3]),a[7][2],1);
    EVAL Glue(Spin(a[3][3]),a[8][2],1);
    EVAL Glue(Spin(a[4][3]),a[9][2],1);
    EVAL Glue(Spin(a[10][1]),a[5][0],1);

    EVAL Glue(Spin(a[6][1]),Enext_1(a[10][7]),1);
    EVAL Glue(Spin(a[11][1]),a[6][0],1);
    EVAL Glue(Spin(a[7][1]),Enext_1(a[11][7]),1);
    EVAL Glue(Spin(a[12][1]),a[7][0],1);
    EVAL Glue(Spin(a[8][1]),Enext_1(a[12][7]),1);
    EVAL Glue(Spin(a[13][1]),a[8][0],1);
    EVAL Glue(Spin(a[9][1]),Enext_1(a[13][7]),1);
    EVAL Glue(Spin(a[14][1]),a[9][0],1);
    EVAL Glue(Spin(a[5][1]),Enext_1(a[14][7]),1);

    (* Last level *)    
    EVAL Glue(Spin(Enext(a[10][0])),a[15][2],1);
    EVAL Glue(Spin(Enext(a[11][0])),a[16][2],1);
    EVAL Glue(Spin(Enext(a[12][0])),a[17][2],1);
    EVAL Glue(Spin(Enext(a[13][0])),a[18][2],1);
    EVAL Glue(Spin(Enext(a[14][0])),a[19][2],1);

    EVAL Glue(Spin(a[16][1]),a[15][0],1);
    EVAL Glue(Spin(a[17][1]),a[16][0],1);
    EVAL Glue(Spin(a[18][1]),a[17][0],1);
    EVAL Glue(Spin(a[19][1]),a[18][0],1);
    EVAL Glue(Spin(a[15][1]),a[19][0],1);


    FOR i :=  0 TO  4 DO b[i] := a[i,2] END;
    FOR i :=  5 TO  9 DO b[i] := a[i,3] END;
    FOR i := 10 TO 14 DO b[i] := a[i,3] END;
    FOR i := 15 TO 19 DO b[i] := a[i,3] END;

    RETURN b;
  END OldMakeIcosahedronTriang;   

PROCEDURE MakeIcosahedronTriang(Original: BOOLEAN) : ARRAY[0..19] OF Pair =
  (* Builds a triangualted icosahedron, with twenty tetrahedra, and return
     one facet-edge pair by original face (i.e. 20). If Original=TRUE the pro-
     cedure emphasize the original elements.
  *)
  VAR
     a : ARRAY[0..19] OF ARRAY[0..7] OF Pair;
     b : ARRAY[0..19] OF Pair;   (* this variable will be retuned *)
  BEGIN
    FOR i := 0 TO 19 DO
      a[i] := MakeTetraTopo(1,1);
    END;

    (* inside the first level *)
    EVAL Glue(Spin(a[1][1]),a[0][0],1);  
    EVAL Glue(Spin(a[2][1]),a[1][0],1); 
    EVAL Glue(Spin(a[3][1]),a[2][0],1);
    EVAL Glue(Spin(a[4][1]),a[3][0],1);
    EVAL Glue(Spin(a[0][1]),a[4][0],1);

    (* between the first and second level *)
    EVAL Glue(Spin(a[0][3]),a[5][2],1);
    EVAL Glue(Spin(a[1][3]),a[6][2],1);
    EVAL Glue(Spin(a[2][3]),a[7][2],1);
    EVAL Glue(Spin(a[3][3]),a[8][2],1);
    EVAL Glue(Spin(a[4][3]),a[9][2],1);
   
    (* inside the second level *)
    EVAL Glue(Spin(Enext(a[10][1])), a[5][0], 1);
    EVAL Glue(Spin(a[6][1]),         Enext(a[10][0]),1);
    EVAL Glue(Spin(Enext(a[11][1])), a[6][0], 1);
    EVAL Glue(Spin(a[7][1]),         Enext(a[11][0]),1);
    EVAL Glue(Spin(Enext(a[12][1])), a[7][0], 1);
    EVAL Glue(Spin(a[8][1]),         Enext(a[12][0]),1);
    EVAL Glue(Spin(Enext(a[13][1])), a[8][0], 1);
    EVAL Glue(Spin(a[9][1]),         Enext(a[13][0]),1);
    EVAL Glue(Spin(Enext(a[14][1])), a[9][0], 1);
    EVAL Glue(Spin(a[5][1]),         Enext(a[14][0]),1); 

    (* between the second and third level *)
    EVAL Glue   (Spin(a[10][3]),        a[15][2],1); 
    EVAL Glue   (Spin(a[11][3]),        a[16][2],1); 
    EVAL Glue   (Spin(a[12][3]),        a[17][2],1); 
    EVAL Glue   (Spin(a[13][3]),        a[18][2],1);
    EVAL Glue   (Spin(a[14][3]),        a[19][2],1); 

    (* inside the third level *)
    EVAL Glue(Spin(a[16][1]),a[15][0],1);
    EVAL Glue(Spin(a[17][1]),a[16][0],1);
    EVAL Glue(Spin(a[18][1]),a[17][0],1);
    EVAL Glue(Spin(a[19][1]),a[18][0],1);
    EVAL Glue(Spin(a[15][1]),a[19][0],1);

    (* Rescue the pairs to be returned *)

    FOR i :=  0 TO  4 DO b[i] := a[i,2] END;
    FOR i :=  5 TO  9 DO b[i] := a[i,3] END;
    FOR i := 10 TO 14 DO b[i] := a[i,2] END;
    FOR i := 15 TO 19 DO b[i] := a[i,3] END;

    IF Original THEN
      (* emphasize the original elements *)

      SetRingFaceExis(a[0][1],FALSE);
      SetEdgeExis(a[0][1],FALSE);

      SetRingFaceExis(a[15][1],FALSE);
      SetEdgeExis(a[15][1],FALSE);
    
      SetVertexExis(a[15][1]);

      FOR i := 0 TO 4 DO
        SetRingFaceExis(Enext(a[i][1]),FALSE);
      END;

      FOR i := 5 TO 9 DO
        SetEdgeExis(Enext_1(a[i][1]),FALSE);
      END;

      FOR i := 15 TO 19 DO
        SetRingFaceExis(Enext_1(a[i][1]),FALSE);
        SetEdgeExis(Enext_1(a[i][1]),FALSE);
      END
    END;

    (* some assertions *)
    FOR i := 0 TO 4 DO
      <* ASSERT Fnext(b[i+5]) = Spin(b[i]) *>
    END;

    FOR i := 10 TO 14 DO
      <* ASSERT Fnext(b[i+5]) = Spin(b[i]) *>
    END;    
    
    RETURN b;
  END MakeIcosahedronTriang;  


PROCEDURE MakeDodecahedronTriang(
  Original: BOOLEAN;
  ) : ARRAY[0..11] OF ARRAY [0..4] OF Pair =
  (* Builds a triangulated Dodecahedron, with 60 tetrahedra, trough the
     automatic gluing of tetrahedra. If Original=TRUE the procedure empha-
     size the original elements.
  *)

TYPE 
  Row4I = ARRAY[0..3] OF CARDINAL;

VAR
  cell4  : REF ARRAY OF Row4I;
  tetra  : REF ARRAY OF ARRAY [0..3] OF Pair;
  cellnum: CARDINAL;
  dode   : ARRAY[0..11] OF ARRAY [0..4] OF Pair;

  PROCEDURE Gluing(Ti,Tj,Ci,Cj: CARDINAL) : Pair =
    (* Gluing the tetrahedra Ti with Tj through the free faces Ci and
       Cj respectively. 
    *)
    BEGIN      
      IF
        (* 1 *) 
        Ci = 0 AND Cj = 0 THEN
        WITH
	  Oi = cell4[Ti,0],
	  Di = cell4[Ti,1],
	  Oj = cell4[Tj,0],
	  Dj = cell4[Tj,1]
	DO 
	  IF (Oi=Dj) AND (Di=Oj) THEN
	    EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN    Clock(tetra[Tj,Cj]);
	  END
	END;
      ELSIF
        (* 2 *)
        Ci = 0 AND Cj = 1 THEN
        WITH
	  Oi = cell4[Ti,0],
	  Di = cell4[Ti,1],
          Oj = cell4[Tj,0],
	  Dj = cell4[Tj,1]
	DO 
	  IF (Oi = Oj) AND (Di = Dj) THEN
	    EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN    Spin(tetra[Tj,Cj]);
	   END
	END;
      ELSIF
        (* 3 *)
        Ci = 0 AND Cj = 2 THEN
        WITH
	  Oi = cell4[Ti,0],
	  Di = cell4[Ti,1],
          Oj = cell4[Tj,2],
	  Dj = cell4[Tj,3]
        DO
	  IF (Oi = Oj) AND (Di = Dj) THEN
	    EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN Spin(tetra[Tj,Cj]);
	  END
        END
      ELSIF
        (* 4 *)
        Ci = 0 AND Cj = 3 THEN
        WITH
	  Oi = cell4[Ti,0],
	  Di = cell4[Ti,1],
          Oj = cell4[Tj,2],
	  Dj = cell4[Tj,3]
        DO
	  IF (Oi = Dj) AND (Di = Oj) THEN
	    EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN    Clock(tetra[Tj,Cj]);
	  END
        END; 
      ELSIF
        (* 5 *)
        Ci = 1 AND Cj = 0 THEN
        WITH
	  Oi = cell4[Ti,0],
	  Di = cell4[Ti,1],
          Oj = cell4[Tj,0],
	  Dj = cell4[Tj,1]
        DO
	  IF (Oi = Oj) AND (Di = Dj) THEN
	    EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN    Spin(tetra[Tj,Cj]);
	  END
        END;
      ELSIF
        (* 6 *)
        Ci = 1 AND Cj = 1 THEN
        WITH
	  Oi = cell4[Ti,0],
	  Di = cell4[Ti,1],
          Oj = cell4[Tj,0],
	  Dj = cell4[Tj,1]
        DO  
	  IF (Oi=Dj) AND (Di=Oj) THEN
	    EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN    Clock(tetra[Tj,Cj]);
	  END
        END;
      ELSIF
        (* 7 *)
        Ci = 1 AND Cj = 2 THEN
        WITH
	  Oi = cell4[Ti,0],
	  Di = cell4[Ti,1],
          Oj = cell4[Tj,2],
	  Dj = cell4[Tj,3]
        DO  
	  IF (Oi=Dj) AND (Di=Oj) THEN
	    EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN    Clock(tetra[Tj,Cj]);
	  END
        END;
      ELSIF
        (* 8 *)
        Ci = 1 AND Cj = 3 THEN
        WITH
	  Oi = cell4[Ti,0],
	  Di = cell4[Ti,1],
          Oj = cell4[Tj,2],
	  Dj = cell4[Tj,3]
        DO
	  IF (Oi = Oj) AND (Di = Dj) THEN
	    EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN    Spin(tetra[Tj,Cj]);
	  END
        END; 
      ELSIF
        (* 9 *)
        Ci = 2 AND Cj = 0 THEN
        WITH
	  Oi = cell4[Ti,2],
	  Di = cell4[Ti,3],
          Oj = cell4[Tj,0],
	  Dj = cell4[Tj,1]
        DO
	  IF (Oi=Oj) AND (Di=Dj) THEN
	    EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN    Spin(tetra[Tj,Cj]);
	  END
        END;
      ELSIF
        (* 10 *)
        Ci = 2 AND Cj = 1 THEN
        WITH
	  Oi = cell4[Ti,2],
	  Di = cell4[Ti,3],
          Oj = cell4[Tj,0],
	  Dj = cell4[Tj,1]
        DO
	  IF (Oi=Dj) AND (Di=Oj) THEN
	    EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN    Clock(tetra[Tj,Cj]);
	  END
        END;
      ELSIF
        (* 11 *)
        Ci = 2 AND Cj = 2 THEN
        WITH
	  Oi = cell4[Ti,2],
	  Di = cell4[Ti,3],
          Oj = cell4[Tj,2], 
	  Dj = cell4[Tj,3]
	DO
	  IF (Oi=Dj) AND (Di=Oj) THEN
	    EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN    Clock(tetra[Tj,Cj]);
	  END
	END;
      ELSIF
        (* 12 *)
        Ci = 2 AND Cj = 3 THEN
        WITH
	  Oi = cell4[Ti,2],
	  Di = cell4[Ti,3],
          Oj = cell4[Tj,2],
	  Dj = cell4[Tj,3]
	DO
	  IF (Oi=Oj) AND (Di=Dj) THEN
	    EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN Spin(tetra[Tj,Cj]);
	  END
	END;
      ELSIF
        (* 13 *)
        Ci = 3 AND Cj = 0 THEN
        WITH
	  Oi = cell4[Ti,2],
	  Di = cell4[Ti,3],
          Oj = cell4[Tj,0],
	  Dj = cell4[Tj,1]
        DO
	  IF (Oi=Dj) AND (Di=Oj) THEN
	    EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN    Clock(tetra[Tj,Cj]);
	  END
        END;
      ELSIF
        (* 14 *) 
        Ci = 3 AND Cj = 1 THEN
        WITH
	  Oi = cell4[Ti,2],
	  Di = cell4[Ti,3],
          Oj = cell4[Tj,0],
	  Dj = cell4[Tj,1]
        DO
	  IF (Oi=Oj) AND (Di=Dj) THEN
	    EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN    Spin(tetra[Tj,Cj]);
	  END
        END;
      ELSIF
        (* 15 *)
        Ci = 3 AND Cj = 2 THEN   
        WITH
	  Oi = cell4[Ti,2],
	  Di = cell4[Ti,3],
	  Oj = cell4[Tj,2],
	  Dj = cell4[Tj,3]
        DO
	  IF (Oi=Oj) AND (Di=Dj) THEN
	    EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN Spin(tetra[Tj,Cj]);
	  END
        END
      ELSIF
        (* 16 *)
        Ci = 3 AND Cj = 3 THEN
        WITH
	  Oi = cell4[Ti,2],
	  Di = cell4[Ti,3],
	  Oj = cell4[Tj,2],
	  Dj = cell4[Tj,3]
        DO
	  IF (Oi=Dj) AND (Di=Oj) THEN
	    EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1);
	    RETURN    Clock(tetra[Tj,Cj]);
	  END
        END
      END;
    RETURN tetra[Ti,0];
  END Gluing;

  PROCEDURE SetCornersTetra(Ti: CARDINAL; row: Row4I) =
    (* Set the labels "row" in the tetrahedron Ti. *)
    BEGIN
      WITH
        a = OrgV(tetra[Ti,0]),
        b = OrgV(Clock(tetra[Ti,0])),
        c = OrgV(Enext_1(tetra[Ti,1])),
        d = OrgV(Enext_1(tetra[Ti,0]))
      DO
        a.num := row[0];
        b.num := row[1];
        c.num := row[2];
        d.num := row[3];
      END;
    END SetCornersTetra;
 
  PROCEDURE SetGhostVertex(Ti: CARDINAL; node:CARDINAL) =
    (* Set one vertex onto the tetrahedra "Ti" as non-existing.*)
    BEGIN
      WITH
        a  = OrgV(tetra[Ti,0]),
        b  = OrgV(Clock(tetra[Ti,0])),
        c  = OrgV(Enext_1(tetra[Ti,1])),
        d  = OrgV(Enext_1(tetra[Ti,0]))
      DO
        IF a.num = node THEN
          a.exists := FALSE;
        END;
        IF b.num = node THEN
          b.exists := FALSE; 
        END; 
        IF c.num = node THEN
          c.exists := FALSE;
        END;
        IF d.num = node THEN
          d.exists := FALSE;
        END
      END
    END SetGhostVertex;

  PROCEDURE SetGhostEdge(Ti: CARDINAL; node:CARDINAL) =
    (* Set all edges incidents in the vertex "node" as  non-existing.*)
    BEGIN
      WITH
        a  = OrgV(tetra[Ti,0]).num,
        b  = OrgV(Clock(tetra[Ti,0])).num,
        c  = OrgV(Enext_1(tetra[Ti,1])).num,
        d  = OrgV(Enext_1(tetra[Ti,0])).num,
        e0 = tetra[Ti,0].facetedge.edge,
        e1 = tetra[Ti,2].facetedge.edge,
        e2 = Enext  (tetra[Ti,2]).facetedge.edge,
        e3 = Enext_1(tetra[Ti,2]).facetedge.edge,
        e4 = Enext  (tetra[Ti,0]).facetedge.edge,
        e5 = Enext  (tetra[Ti,1]).facetedge.edge
      DO
        IF a = node OR b = node THEN
          e0.exists := FALSE;
        END;
        IF c = node OR d = node THEN
          e1.exists := FALSE; 
        END; 
        IF a = node OR d = node THEN
          e2.exists := FALSE;
        END;
        IF a = node OR c = node THEN
          e3.exists := FALSE;
        END;         
        IF b = node OR d = node THEN
          e4.exists := FALSE;
        END;
        IF b = node OR c = node THEN
          e5.exists := FALSE; 
        END;
      END 
    END SetGhostEdge;

    PROCEDURE SaveFreeCorners() =
      (* Save the free corners on the triangulated dodecahedron, such
         allow us to use this polyhedron as the unglued scheme for the
         3D maps obtained by the gluing of opposite faces on the dode-
         cahedron.
      *)

      PROCEDURE AssignmentOnFace(b: Pair; i: CARDINAL) =
        BEGIN
	 IF Octf.SpinBit(b) = 0 THEN
           dode[i,0] := b;
	   FOR j := 1 TO 4 DO
	     dode[i,j] := Clock(Enext(Fnext(Enext_1(dode[i,j-1]))));
            END;
         ELSIF Octf.SpinBit(b) = 1 THEN
	   dode[i,0] := b;
           FOR j := 1 TO 4 DO  
	     dode[i,j] := Clock(Enext_1(Fnext(Enext(dode[i,j-1]))));
       	   END
         END
      END AssignmentOnFace;

      VAR m,n: CARDINAL := 0;
      BEGIN
        FOR i := 0 TO 59 DO
          FOR j := 0 TO 3 DO
            WITH
              a = tetra[i,j],
              d = Octf.DegreeFaceRing(a)
            DO
              IF d # 1 AND Ppos(a) = NIL THEN
                dode[m,n] := tetra[i,j];
                INC(n);
                IF n =  5 THEN INC(m); n := 0; END;
              END 
            END
          END
        END;

        WITH a = dode[0,0] DO AssignmentOnFace(a,0) END;
 
        dode[1,0] := Spin(Fnext(dode[0,0])); AssignmentOnFace(dode[1,0],1);
        dode[2,0] := Spin(Fnext(dode[0,4])); AssignmentOnFace(dode[2,0],2);
        dode[3,0] := Spin(Fnext(dode[0,3])); AssignmentOnFace(dode[3,0],3);
        dode[4,0] := Spin(Fnext(dode[0,2])); AssignmentOnFace(dode[4,0],4);
        dode[5,0] := Spin(Fnext(dode[0,1])); AssignmentOnFace(dode[5,0],5);
     
        dode[8, 0] := Spin(Fnext(dode[4,2])); AssignmentOnFace(dode[8,0],8); 
        dode[6, 0] := Spin(Fnext(dode[8,3])); AssignmentOnFace(dode[6,0],6);
        dode[7, 0] := Spin(Fnext(dode[6,1])); AssignmentOnFace(dode[7,0], 7);
        dode[9, 0] := Spin(Fnext(dode[6,4])); AssignmentOnFace(dode[9,0], 9);
        dode[10,0] := Spin(Fnext(dode[6,3])); AssignmentOnFace(dode[10,0],10);
        dode[11,0] := Spin(Fnext(dode[6,2])); AssignmentOnFace(dode[11,0],11);

        (* the free corners are computed follows the scheme bellow:

             /|\                             /|\
            / | \                           / | \
           /  |  \                         /  |  \
          / 3 | 2 \                       / 3 | 2 \
         /    |    \                     /    |    \
        /-----|-----\                   /-----|-----\
        \  4 / \ 1  /                   \  4 / \ 1  /
         \  / 0 \  /                     \  / 0 \  /   
          \/_____\/                       \/_____\/ 

             <---                            --->
                                             
          SpinBit = 0                     SpinBit = 1

       *)

      END SaveFreeCorners;

  PROCEDURE SetGhostFace(Ti: CARDINAL; node: CARDINAL) =
    (* Set all facese incidents in the vertex "node" as  non-existing.*)
    BEGIN
      WITH
        a  = OrgV(tetra[Ti,0]).num,
        b  = OrgV(Clock(tetra[Ti,0])).num,
        c  = OrgV(Enext_1(tetra[Ti,1])).num,
        d  = OrgV(Enext_1(tetra[Ti,0])).num,
        f0 = tetra[Ti,0].facetedge.face,
        f1 = tetra[Ti,1].facetedge.face,
        f2 = tetra[Ti,2].facetedge.face,
        f3 = tetra[Ti,3].facetedge.face
      DO
        IF    a = node OR b = node OR d = node THEN
          f0.exists := FALSE;
        END;
        IF a = node OR b = node OR c = node THEN
          f1.exists := FALSE;
        END; 
        IF a = node OR c = node OR d = node THEN
          f2.exists := FALSE;
        END;
        IF b = node OR c = node OR d = node THEN
          f3.exists := FALSE; 
        END
      END 
    END SetGhostFace;

  PROCEDURE MustBeGlue(Ti,Tj: Pair) : BOOLEAN =
    (* Return TRUE if the faces "Ti.facetedge.face" and
       "Tj.facetedge.face" have coherent orientations and
       must be glued.
    *)
    BEGIN
      WITH
        a    = OrgV(Ti).num,
        ae   = OrgV(Enext(Ti)).num,
        ae_1 = OrgV(Enext_1(Ti)).num,
        b    = OrgV(Tj).num,
        be   = OrgV(Enext(Tj)).num,
        be_1 = OrgV(Enext_1(Tj)).num
      DO
        IF (a = b AND ae = be AND ae_1 = be_1) OR
	   (a = b AND ae = be_1 AND ae_1 = be) THEN RETURN TRUE END;
        RETURN FALSE
      END;  
    END MustBeGlue;

  PROCEDURE EnextK(Ti: Pair; k : CARDINAL) : Pair =
    (* Given a pair "Ti", this procedure return Enext^{k}(Ti). *)
    BEGIN
      IF k = 0 THEN RETURN Ti
      ELSIF k = 1 THEN RETURN Enext(Ti)
      ELSIF k = 2 THEN RETURN Enext(Enext(Ti))
      END;
      RETURN Ti;
    END EnextK;

    VAR
      poly  : REF ARRAY OF ARRAY [0..7] OF Pair;
      count : CARDINAL := 1;
      faces : REF ARRAY OF Pair;
      glues : REF ARRAY OF Row4I;
    BEGIN
     
      cellnum := 60;
      cell4   := NEW(REF ARRAY OF Row4I, cellnum);
      
      poly  := NEW(REF ARRAY OF ARRAY [0..7] OF Pair,cellnum);
      tetra := NEW(REF ARRAY OF ARRAY [0..3] OF Pair,cellnum);
        
      faces := NEW(REF ARRAY OF Pair,  4*cellnum);
      glues := NEW(REF ARRAY OF Row4I, 2*cellnum);

      (* creating topological tetrahedra *)
      FOR i := 0 TO cellnum-1 DO
	poly[i] := MakeTetraTopo(1,1);
      END;

      (* creating the tetrahedra *)
      FOR i := 0 TO cellnum-1 DO
	FOR j := 0 TO 3 DO
	  tetra[i,j] := poly[i,j];
	  <* ASSERT Ppos(tetra[i,j]) = NIL *>
	END
      END;
        
      (* cells with corners perfectly assigments *)
	cell4[  0]:=Row4I{100,500,  0,  1};
	cell4[  1]:=Row4I{100,500,  1,  4};
	cell4[  2]:=Row4I{100,500,  4,  7};
	cell4[  3]:=Row4I{100,500,  7,  2};
	cell4[  4]:=Row4I{100,500,  2,  0};
	cell4[  5]:=Row4I{500,101,  0,  1};
	cell4[  6]:=Row4I{500,101,  1,  5};
	cell4[  7]:=Row4I{500,101,  5,  8};
	cell4[  8]:=Row4I{500,101,  8,  3};
	cell4[  9]:=Row4I{500,101,  3,  0};
	cell4[ 10]:=Row4I{500,102,  2,  0};
	cell4[ 11]:=Row4I{500,102,  0,  3};
	cell4[ 12]:=Row4I{500,102,  3,  9};
	cell4[ 13]:=Row4I{500,102,  9,  6};
	cell4[ 14]:=Row4I{500,102,  6,  2};
	cell4[ 15]:=Row4I{103,500,  6,  2};
	cell4[ 16]:=Row4I{500,103,  7,  2};
	cell4[ 17]:=Row4I{103,500,  7, 13};
	cell4[ 18]:=Row4I{103,500, 13, 12};
	cell4[ 19]:=Row4I{500,103,  6, 12};
	cell4[ 20]:=Row4I{500,104,  1,  4};
	cell4[ 21]:=Row4I{500,104,  4, 10};
	cell4[ 22]:=Row4I{500,104, 10, 11};
	cell4[ 23]:=Row4I{500,104, 11,  5};
	cell4[ 24]:=Row4I{500,104,  5,  1};
	cell4[ 25]:=Row4I{500,105,  9,  3};
	cell4[ 26]:=Row4I{500,105,  3,  8};
	cell4[ 27]:=Row4I{500,105,  8, 14};
	cell4[ 28]:=Row4I{500,105, 14, 15};
	cell4[ 29]:=Row4I{500,105, 15,  9};
	cell4[ 30]:=Row4I{500,106, 12,  6};
	cell4[ 31]:=Row4I{500,106,  6,  9};
	cell4[ 32]:=Row4I{500,106,  9, 15};
	cell4[ 33]:=Row4I{500,106, 15, 18};
	cell4[ 34]:=Row4I{500,106, 18, 12};
	cell4[ 35]:=Row4I{500,107,  4,  7};
	cell4[ 36]:=Row4I{500,107,  7, 13};
	cell4[ 37]:=Row4I{500,107, 13, 16};
	cell4[ 38]:=Row4I{500,107, 16, 10};
	cell4[ 39]:=Row4I{500,107, 10,  4};
	cell4[ 40]:=Row4I{108,500,  5,  8};
	cell4[ 41]:=Row4I{108,500,  8, 14};
	cell4[ 42]:=Row4I{108,500, 14, 17};
	cell4[ 43]:=Row4I{108,500, 17, 11};
	cell4[ 44]:=Row4I{108,500, 11,  5};
	cell4[ 45]:=Row4I{109,500, 18, 12};
	cell4[ 46]:=Row4I{109,500, 12, 13};
	cell4[ 47]:=Row4I{109,500, 13, 16};
	cell4[ 48]:=Row4I{109,500, 16, 19};
	cell4[ 49]:=Row4I{109,500, 19, 18};
	cell4[ 50]:=Row4I{500,110, 18, 15};
	cell4[ 51]:=Row4I{500,110, 15, 14};
	cell4[ 52]:=Row4I{500,110, 14, 17};
	cell4[ 53]:=Row4I{500,110, 17, 19};
	cell4[ 54]:=Row4I{500,110, 19, 18};
	cell4[ 55]:=Row4I{111,500, 10, 11};
	cell4[ 56]:=Row4I{111,500, 11, 17};
	cell4[ 57]:=Row4I{111,500, 17, 19};
	cell4[ 58]:=Row4I{111,500, 19, 16};
	cell4[ 59]:=Row4I{111,500, 16, 10};


      (* set the labels for each tetrahedra *)
      FOR i := 0 TO cellnum-1 DO
        SetCornersTetra(i,cell4[i]);
      END;

      IF Original THEN
        (* emphasize the original elements *)
      
        FOR i := 0 TO cellnum-1 DO
          SetGhostVertex(i,500);
          SetGhostEdge(i,500);
          SetGhostFace(i,500);
          FOR j := 100 TO 111 DO
            SetGhostVertex(i,j);
            SetGhostEdge(i,j);
          END
        END
      END; 

      (* builds the table of faces for choose which tetrahedra
          must be gluing.
      *)

      FOR i := 0 TO cellnum-1 DO
        FOR k := 0 TO 3 DO
	  faces[(4*i)+k] := tetra[i,k];
        END
      END;
		
      (* computing which cells must be gluing. *)       	
      FOR k := 0 TO LAST(faces^) DO
        FOR l := k+1 TO LAST(faces^) DO
	  FOR m := 0 TO 2 DO
	    WITH
	      e = EnextK(faces[l],m)
	    DO
	      IF MustBeGlue(faces[k],e) THEN
	        WITH
		  kc = k MOD 4, 
		  kt = k DIV 4,
		  lc = l MOD 4,
		  lt = l DIV 4
	        DO
		  glues[count-1] := Row4I{kt,lt,kc,lc};
		  INC(count);
	        END
	      END      
	    END
	  END   
        END
      END;
       
     (* Do the automatic gluing of tetrahedra *)
     FOR i := 0 TO LAST(glues^) DO
       WITH
	 c = glues[i]
       DO
         IF c # Row4I{0,0,0,0} THEN
           EVAL Gluing(c[0],c[1],c[2],c[3]);
         END
       END
     END;

     (* setting the origins. *)
     FOR i := 0 TO cellnum-1 DO
       FOR j := 0 TO 3 DO
	 WITH
	   a = tetra[i,j],
	   b = Enext(a),
	   c = Enext_1(a)
	 DO
	   Triangulation.SetAllOrgs(a,OrgV(a));
	   Triangulation.SetAllOrgs(b,OrgV(b));
	   Triangulation.SetAllOrgs(c,OrgV(c));
	 END
       END
     END;

     SaveFreeCorners();
        
     RETURN dode;
END MakeDodecahedronTriang;


PROCEDURE MakePyramid(n: CARDINAL) : REF ARRAY OF Pair =
  (* Builds a pyramid with base as a n-gon face, without polyhedron
     information.
  *)
  VAR a := NEW(REF ARRAY OF Pair,n);
  BEGIN
    WITH
      b = MakeGonFull(n)
    DO
      (* creating the "n" triangular faces *)
      FOR i := 0 TO n-1 DO
        a[i] := MakeTriangle();
      END;

      (* doing the connections between the base and the triangular faces *)
      FOR i := 0 TO n-1 DO
        SetFnext(a[i],b[i]);
      END;  

      (* doing the connections between the triangular faces *)
      FOR i := 0 TO n-1 DO
        SetFnext( Enext(a[i]), Clock(Enext_1(a[(i+1) MOD n ])) );
      END;

      (* setting the common elements: edges and vertices  *)
      FOR j := 0 TO n-1 DO
        WITH
          a = a[j],
          b = Enext(a),
          c = Enext_1(a)
        DO
          SetEdgeAll(a,a.facetedge.edge); SetAllOrgs(a,OrgV(a));
          SetEdgeAll(b,b.facetedge.edge); SetAllOrgs(b,OrgV(b));
          SetEdgeAll(c,c.facetedge.edge); SetAllOrgs(c,OrgV(c));
        END
      END;

      RETURN b;
    END
  END MakePyramid;

PROCEDURE SetNextPnegC(a : Pair; n: Node) =
  (* Set the 24 pairs facetedges belonging to same polyhedron negative "n"
     for one cubic cell. *)
  VAR t : Pair := a;
  BEGIN
    SetAllPneg(t,n);
    REPEAT
      SetAllPneg(Clock(Fnext_1(t)),n);
      t := Enext_1(t); 
    UNTIL t = a; 
    t := Fnext_1(Enext_1(Enext_1(Fnext_1(a))));
    SetAllPneg(t,n);
  END SetNextPnegC;

PROCEDURE CubeNegVertices(a: Pair): ARRAY [0..7] OF Vertex =
  (* Valid for the versions Spin(a), Clock(a) and SpinClock(a). ? *)
  BEGIN
    <* ASSERT Pneg(a) # NIL *>
    WITH
      p = OrgV(a),
      q = OrgV(Enext_1(a)), 
      r = OrgV(Enext_1(Enext_1(a))),
      s = OrgV(Enext_1(Fnext_1(Enext_1(a)))),
      b = Clock(Fnext(Enext_1(Enext_1(Fnext_1(a))))),
      t = OrgV(b),
      u = OrgV(Enext_1(b)), 
      v = OrgV(Enext_1(Enext_1(b))),
      w = OrgV(Enext_1(Fnext_1(Enext_1(b))))
    DO
      <* ASSERT Pneg(a) = Ppos(b) *>
      <* ASSERT (p # q) AND (q # r) AND (r # s) *>
      <* ASSERT  s # t *>
      <* ASSERT (t # u) AND (u # v) AND (v # w) *>
      RETURN ARRAY [0..7] OF Vertex{p,q,r,s,t,u,v,w}
    END
  END CubeNegVertices;

PROCEDURE CubeBarycenter(a: Pair; READONLY c: Coords): LR4.T =
  VAR
    n : CARDINAL := 0;
    sum := LR4.T{0.0d0, ..};
  BEGIN
    WITH
      cube = CubeNegVertices(a)
    DO
      FOR i := 0 TO 7 DO
        WITH
         aoc = c[cube[i].num]
        DO	
          sum := LR4.Add(sum, aoc);
	  INC(n);
        END
      END
    END;
    RETURN LR4.Scale(1.0d0/FLOAT(n,LONGREAL), sum)
  END CubeBarycenter;

PROCEDURE SetNextPnegB(a : Pair; n: Node) =
  (* Set the 8 pairs facetedges belonging to same polyhedron negative "n"
     for one cigar-shaped cell. *)
  VAR t : Pair := a;
  BEGIN
    SetAllPneg(t,n);
    REPEAT
      SetAllPneg(Clock(Fnext_1(t)),n);
      t := Enext_1(t); 
    UNTIL t = a; 
    t := Fnext_1(Enext_1(Fnext_1(a)));
    SetAllPneg(t,n);
  END SetNextPnegB;

<* UNUSED *>
PROCEDURE PrtNextPnegB(a: Pair) =
  (* Print the 8 pairs facetedges belonging to same polyhedron negative
     for one cubic cell. *)

  PROCEDURE NextPnegof(a: Pair) =
    VAR t : Pair := a;
    BEGIN
      REPEAT
        Tools.PrtPnegNum(t);
        t := Enext_1(t);
      UNTIL t = a;
    END NextPnegof;

  VAR t : Pair := a;
  BEGIN
    NextPnegof(t);
    REPEAT
      NextPnegof(Clock(Fnext_1(t)));
      t := Enext_1(t);
    UNTIL (t=a);
    t := Fnext_1(Enext_1(Fnext_1(a)));
    NextPnegof(t);
  END PrtNextPnegB;
  
<* UNUSED *>
PROCEDURE PrtPnegPpos(a: Pair) =
  (* Print the polyhedron negative and positive of pairs facetedges
     belonging to same face. *)
  VAR t : Pair := a;
  BEGIN
    REPEAT
      Tools.PrtPnegNum(t);
      Tools.PrtPposNum(t);   
      Tools.PrtPnegNum(Clock(t));
      Tools.PrtPposNum(Clock(t));
      t := Enext_1(t);
    UNTIL t = a;
  END PrtPnegPpos;

<* UNUSED *> 
PROCEDURE PrtNextPnegC(a: Pair) =
  (* Print the 24 pairs facetedges belonging to same polyhedron negative
     for one cubic cell. *)

  PROCEDURE NextPnegof(a: Pair) =
    VAR t : Pair := a;
    BEGIN
      REPEAT
        Tools.PrtPnegNum(t);
        t := Enext_1(t);
      UNTIL t = a;
  END NextPnegof;

  VAR t : Pair := a;
  BEGIN
    NextPnegof(t);
    REPEAT
      NextPnegof(Clock(Fnext_1(t)));
      t := Enext_1(t);
    UNTIL t = a;
    t := Fnext_1(Enext_1(Enext_1(Fnext_1(a))));
    NextPnegof(t);
  END PrtNextPnegC;

PROCEDURE MakeBall(): ARRAY [0..3] OF Pair =
  (* Build one polyhedron with American football shape. This polyhedron
     is the gluing scheme for obtain the "pseudomanifold" complex.
  *)
  VAR
    co : ARRAY [0..3] OF Pair;
  BEGIN
    FOR i := 0 TO 3 DO
      co[i] := MakeGon(2);
    END;

    (* Glue Faces *)
    FOR i := 0 TO 3 DO
      SetFnext(co[i],Clock(Enext(co[(i+1) MOD 4]))); 
      SetEdge(co[i], Clock(Enext(co[(i+1) MOD 4])).facetedge.edge);
      SetAllOrgs(co[i], Org(co[i]));
      SetAllOrgs(Clock(co[i]), Org(Clock(co[i]))); 
    END;

    WITH
      q = MakePolyhedron()
    DO
      q.num := PolyhedronCount; INC(PolyhedronCount);
      SetNextPnegB(co[0],q);
    END;

    RETURN co;
  END MakeBall;


PROCEDURE MakeTetrahedron() : ARRAY [0..3] OF Pair =
  VAR
    co : ARRAY [0..3] OF Pair;
  BEGIN
    FOR i := 0 TO 3 DO
      co[i] := MakeTriangle();
    END;

    (* Glue Faces *)
    SetFnext(co[3],Clock(co[0]));
    SetEdge(co[3],Clock(co[0]).facetedge.edge);

    SetFnext(Clock(co[1]), Enext(co[0]));
    SetEdge(Clock(co[1]), Enext(co[0]).facetedge.edge);
    
    SetFnext(Clock(Enext_1(co[0])), co[2]);
    SetEdge(Clock(Enext_1(co[0])), co[2].facetedge.edge);

    SetFnext(Clock(Enext_1(co[1])), Enext(co[2]));
    SetEdge(Clock(Enext_1(co[1])), Enext(co[2]).facetedge.edge);

    SetFnext(Clock(Enext(co[1])), Enext_1(co[3]));
    SetEdge(Clock(Enext(co[1])),  Enext_1(co[3]).facetedge.edge);

    SetFnext(Clock(Enext_1(co[2])), Enext(co[3]));
    SetEdge(Clock(Enext_1(co[2])),  Enext(co[3]).facetedge.edge);

    SetAllOrgs(co[0], Org(co[0]));
    SetAllOrgs(Clock(co[0]), Org(Clock(co[0])));
    SetAllOrgs(co[1], Org(co[1]));
    SetAllOrgs(Enext_1(co[1]), Org(Enext_1(co[1])));

    WITH
      q = MakePolyhedron()
    DO
      q.num := PolyhedronCount; INC(PolyhedronCount);
      SetAllPneg(co[0],q);
    END;

    RETURN co;
  END MakeTetrahedron;

PROCEDURE GlueBall(a,b: Pair) : Pair =
  (* Make the glueing of two simples 2-gon faces around of one common face.
     The pair "a" and theirs adjacents with the same face component is
     killed. *)
  BEGIN
    <* ASSERT Octf.SpinBit(a) = Octf.SpinBit(b) *>
    Octf.Meld(a,b);

    (* updating edges relations *)
    SetEdgeAll(b, b.facetedge.edge);
    SetEdgeAll(Enext(b), Enext(b).facetedge.edge);
    
    (* updating vertices relations *)
    SetAllOrgs(b, Org(b));
    SetAllOrgs(Enext(b), Org(Enext(b)));
    
    (* updating polyhedron relations *)
    SetPneg(b, Pneg(a));
    SetPneg(Enext_1(b), Pneg(Enext_1(a)));

    RETURN b;
  END GlueBall;

PROCEDURE MakeCube() : ARRAY [0..5] OF Pair =             

  PROCEDURE CEE(a: Pair) : Pair =
    BEGIN
      RETURN Clock(Enext(Enext(a)));
    END CEE;

  VAR
    co : ARRAY [0..5] OF Pair;
  BEGIN
    FOR i := 0 TO 5 DO
      co[i] := MakeSquare();
    END;

    (* Glue Faces *)
    SetFnext(co[1],CEE(co[2]));                                 (* 1-2 *)
    SetEdge(co[1], CEE(co[2]).facetedge.edge);
    SetAllOrgs(co[1], Org(co[1]));
    SetAllOrgs(Clock(co[1]), Org(Clock(co[1])));
   
    SetFnext(co[2],CEE(co[3]));                                 (* 2-3 *)
    SetEdge(co[2], CEE(co[3]).facetedge.edge);
    SetAllOrgs(co[2], Org(co[2]));
    SetAllOrgs(Clock(co[2]), Org(Clock(co[2])));
    
    SetFnext(co[3],CEE(co[4]));                                 (* 3-4 *)
    SetEdge(co[3], CEE(co[4]).facetedge.edge);
    SetAllOrgs(co[3], Org(co[3]));
    SetAllOrgs(Clock(co[3]), Org(Clock(co[3])));

    SetFnext(co[4],CEE(co[1]));                                 (* 4-1 *)
    SetEdge(co[4], CEE(co[1]).facetedge.edge);
    SetAllOrgs(co[4], Org(co[4]));
    SetAllOrgs(Clock(co[4]), Org(Clock(co[4])));
    
    SetFnext(Enext_1(co[1]), Clock(Enext(co[0])));              (* 1-0 *)
    SetEdge(Enext_1(co[1]), Clock(Enext(co[0])).facetedge.edge);
    SetAllOrgs(Enext_1(co[1]), Org(Enext_1(co[1])));
    SetAllOrgs(Clock(Enext_1(co[1])), Org(Clock(Enext_1(co[1]))));

    SetFnext(Enext_1(co[5]), Clock(Enext(co[1])));              (* 5-1 *)
    SetEdge(Enext_1(co[5]), Clock(Enext(co[1])).facetedge.edge);
    SetAllOrgs(Enext_1(co[5]), Org(Enext_1(co[5])));
    SetAllOrgs(Clock(Enext_1(co[5])), Org(Clock(Enext_1(co[5]))));

    SetFnext(co[5], Clock(Enext(co[2])));                       (* 5-2 *)
    SetEdge(co[5], Clock(Enext(co[2])).facetedge.edge);
    SetAllOrgs(co[5], Org(co[5]));
    SetAllOrgs(Clock(co[5]), Org(Clock(co[5])));

    SetFnext(Enext(co[5]), Clock(Enext(co[3])));                (* 5-3 *)
    SetEdge(Enext(co[5]), Clock(Enext(co[3])).facetedge.edge);
    SetAllOrgs(Enext(co[5]), Org(Enext(co[5])));
    SetAllOrgs(Clock(Enext(co[5])), Org(Clock(Enext(co[5]))));

    SetFnext(Enext(Enext(co[5])), Clock(Enext(co[4])));         (* 5-4 *)
    SetEdge(Enext(Enext(co[5])), Clock(Enext(co[4])).facetedge.edge);
    SetAllOrgs(Enext(Enext(co[5])), Org(Enext(Enext(co[5]))));
    SetAllOrgs(Clock(Enext(Enext(co[5]))), Org(Clock(Enext(Enext(co[5])))));

    SetFnext(co[0], Clock(Enext_1(co[2])));                     (* 0-2 *)
    SetEdge(co[0], Clock(Enext_1(co[2])).facetedge.edge);
    SetAllOrgs(co[0], Org(co[0]));
    SetAllOrgs(Clock(co[0]), Org(Clock(co[0])));

    SetFnext(Enext_1(co[0]), Clock(Enext_1(co[3])));            (* 0-3 *)
    SetEdge(Enext_1(co[0]), Clock(Enext_1(co[3])).facetedge.edge);
    SetAllOrgs(Enext_1(co[0]), Org(Enext_1(co[0])));
    SetAllOrgs(Clock(Enext_1(co[0])), Org(Clock(Enext_1(co[0])))); 

    SetFnext(Enext(Enext(co[0])), Clock(Enext_1(co[4])));       (* 0-4 *)
    SetEdge(Enext(Enext(co[0])), Clock(Enext_1(co[4])).facetedge.edge);
    SetAllOrgs(Enext(Enext(co[0])), Org(Enext(Enext(co[0]))));
    SetAllOrgs(Clock(Enext(Enext(co[0]))), Org(Clock(Enext(Enext(co[0])))));

    WITH
      q = MakePolyhedron()
    DO
      q.num := PolyhedronCount; INC(PolyhedronCount);
      SetNextPnegC(co[0],q);
    END;

    RETURN co;
  END MakeCube;

PROCEDURE GlueCube(a,b : Pair) : Pair =
  (* Make the glueing of two simples cubes around of one squared face
     common. The pair "a" and theirs adjacents with the same face component is
     killed. *)
  BEGIN
    Octf.Meld(a,b);

    (* updating edges relations *)
    SetEdgeAll(b, b.facetedge.edge);
    SetEdgeAll(Enext(b), Enext(b).facetedge.edge);
    SetEdgeAll(Enext(Enext(b)), Enext(Enext(b)).facetedge.edge);
    SetEdgeAll(Enext(Enext(Enext(b))), Enext(Enext(Enext(b))).facetedge.edge);

    (* updating vertices relations *)
    SetAllOrgs(b, Org(b));
    SetAllOrgs(Enext(b), Org(Enext(b)));
    SetAllOrgs(Enext(Enext(b)), Org(Enext(Enext(b))));
    SetAllOrgs(Enext(Enext(Enext(b))), Org(Enext(Enext(Enext(b)))));
    SetAllOrgs(Clock(b), Org(Clock(b)));
    SetAllOrgs(Clock(Enext(b)), OrgV(Clock(Enext(b))));
    SetAllOrgs(Clock(Enext(Enext(b))), Org(Clock(Enext(Enext(b)))));
    SetAllOrgs(Clock(Enext(Enext(Enext(b)))),
               Org(Clock(Enext(Enext(Enext(b))))));

    (* updating polyhedron relations *)
    SetPneg(b, Pneg(a));
    SetPneg(Enext_1(b), Pneg(Enext_1(a)));
    SetPneg(Enext_1(Enext_1(b)), Pneg(Enext_1(Enext_1(a))));
    SetPneg(Enext_1(Enext_1(Enext_1(b))), Pneg(Enext_1(Enext_1(Enext_1(a)))));

    RETURN b;
  END GlueCube;

PROCEDURE GlueBigCube(a, b : Pair; n: CARDINAL) : Pair =
  VAR ta,tb: ARRAY [0..20] OF ARRAY [0..20] OF Pair;
  (* Make the glueing of two simples "bigcubes" around of one squared grid
     common face. *)
  BEGIN
    (* sanity check *)
    <* ASSERT n >= 1 *>
    ta[0,0] := a; tb[0,0] := b;

    FOR i := 1 TO n-1 DO
      ta[i,0] := Clock(Fnext_1(Enext(Enext(ta[i-1,0]))));
      tb[i,0] := Clock(Fnext(Enext(Enext(tb[i-1,0]))));
    END;

    FOR i := 0 TO n-1 DO
      FOR j := 1 TO n-1 DO
        ta[i,j] := Clock(Enext_1(Fnext_1(Enext(ta[i, j-1]))));
        tb[i,j] := Clock(Enext_1(Fnext(Enext(tb[i,j-1]))));
        <* ASSERT ta[i,j] # a *> 
        <* ASSERT tb[i,j] # b *>
      END;
    END;

    FOR i := 0 TO n-1 DO
      Octf.Meld(tb[i,0], ta[i,0]);
     
      (* updating edges relations *)
      SetEdgeAll(ta[i,0], ta[i,0].facetedge.edge);
      SetEdgeAll(Enext(ta[i,0]), Enext(ta[i,0]).facetedge.edge);
      SetEdgeAll(Enext(Enext(ta[i,0])), Enext(Enext(ta[i,0])).facetedge.edge);
      SetEdgeAll(Enext(Enext(Enext(ta[i,0]))),
                 Enext(Enext(Enext(ta[i,0]))).facetedge.edge);

      (* updating vertices relations *)
      SetAllOrgs(ta[i,0], Org(ta[i,0]));
      SetAllOrgs(Enext(ta[i,0]), Org(Enext(ta[i,0])));
      SetAllOrgs(Enext(Enext(ta[i,0])), Org(Enext(Enext(ta[i,0]))));
      SetAllOrgs(Enext(Enext(Enext(ta[i,0]))),
                 Org(Enext(Enext(Enext(ta[i,0])))));
      SetAllOrgs(Clock(ta[i,0]), Org(Clock(ta[i,0])));
      SetAllOrgs(Clock(Enext(ta[i,0])), OrgV(Clock(Enext(ta[i,0]))));
      SetAllOrgs(Clock(Enext(Enext(ta[i,0]))),
                 Org(Clock(Enext(Enext(ta[i,0])))));
      SetAllOrgs(Clock(Enext(Enext(Enext(ta[i,0])))),
                 Org(Clock(Enext(Enext(Enext(ta[i,0]))))));

      SetPneg(ta[i,0], Pneg(tb[i,0]));
      SetPneg(Enext_1(ta[i,0]), Pneg(Enext_1(tb[i,0])));
      SetPneg(Enext_1(Enext_1(ta[i,0])), Pneg(Enext_1(Enext_1(tb[i,0]))));
      SetPneg(Enext_1(Enext_1(Enext_1(ta[i,0]))),
              Pneg(Enext_1(Enext_1(Enext_1(tb[i,0])))));
		 
    END;
    
    FOR i := 0 TO n-1 DO
      FOR j := 1 TO n-1 DO
        Octf.Meld(tb[i,j],ta[i,j]);

        (* updating edges relations *)
        SetEdgeAll(ta[i,j], ta[i,j].facetedge.edge);
        SetEdgeAll(Enext(ta[i,j]), Enext(ta[i,j]).facetedge.edge);
        SetEdgeAll(Enext(Enext(ta[i,j])),
                             Enext(Enext(ta[i,j])).facetedge.edge);
        SetEdgeAll(Enext(Enext(Enext(ta[i,j]))),
                   Enext(Enext(Enext(ta[i,j]))).facetedge.edge);

        (* updating vertices relations *)
        SetAllOrgs(ta[i,j], Org(ta[i,j]));
        SetAllOrgs(Enext(ta[i,j]), Org(Enext(ta[i,j])));
        SetAllOrgs(Enext(Enext(ta[i,j])), Org(Enext(Enext(ta[i,j]))));
        SetAllOrgs(Enext(Enext(Enext(ta[i,j]))),
                                   Org(Enext(Enext(Enext(ta[i,j])))));
        SetAllOrgs(Clock(ta[i,j]), Org(Clock(ta[i,j])));
        SetAllOrgs(Clock(Enext(ta[i,j])), OrgV(Clock(Enext(ta[i,j]))));
        SetAllOrgs(Clock(Enext(Enext(ta[i,j]))),
                                     Org(Clock(Enext(Enext(ta[i,j])))));
        SetAllOrgs(Clock(Enext(Enext(Enext(ta[i,j])))),
                   Org(Clock(Enext(Enext(Enext(ta[i,j]))))));

        SetPneg(ta[i,j], Pneg(tb[i,j]));
        SetPneg(Enext_1(ta[i,j]), Pneg(Enext_1(tb[i,j])));
        SetPneg(Enext_1(Enext_1(ta[i,j])), Pneg(Enext_1(Enext_1(tb[i,j]))));
        SetPneg(Enext_1(Enext_1(Enext_1(ta[i,j]))),
                Pneg(Enext_1(Enext_1(Enext_1(tb[i,j])))));		   
      END;
    END;

    RETURN ta[0,0];
  END GlueBigCube;

PROCEDURE MakeRowCube(order: CARDINAL) : REF ARRAY OF PAIRS =
(* Build one row of cubes with order "order". *)
  VAR
    ca : REF ARRAY OF ARRAY [0..5] OF Pair;
    cb : REF ARRAY OF PAIRS;
  BEGIN
    ca := NEW(REF ARRAY OF ARRAY [0..5] OF Pair, order);
    cb := NEW(REF ARRAY OF PAIRS, order);
    FOR i := 0 TO order-1 DO
      ca[i] := MakeCube();
      cb[i,0] := ca[i,1]; 
      cb[i,1] := ca[i,3]; 
      cb[i,2] := ca[i,0]; 
      cb[i,3] := ca[i,5];
    END;

    (* gluing *)
    FOR j := 0 TO order-2 DO
      EVAL GlueCube(ca[j,2], Fnext_1(ca[j+1,3]));
    END;
    RETURN cb;
  END MakeRowCube;

PROCEDURE MakeColumnCube(order: CARDINAL) : REF TriPair =
(* Builds one bidimensional array of cubes of order: "order x order". *)
  VAR
    cb : REF ARRAY OF ARRAY OF PAIRS;
    cc : REF TriPair;
  BEGIN
    cb := NEW(REF ARRAY OF ARRAY OF PAIRS, order, order);
    cc := NEW(REF ARRAY OF ARRAY OF ARRAY [0..1] OF Pair, order,order);

    (* Mounting the array cb *)
    FOR i := 0 TO order-1 DO
      WITH ca = MakeRowCube(order) DO
        FOR j := 0 TO order-1 DO
          FOR k := 0 TO 3 DO
            cb[i,j,k] := ca[j,k];
          END
        END
      END
    END;

    (* gluing *)
    FOR j := 0 TO order-2 DO
      FOR k := 0 TO order-1 DO
        EVAL GlueCube(cb[j,k,1],Clock(Enext(Enext(cb[j+1,k,0]))));
      END
    END;

    (* Selecting the pairs for return *)
    FOR i := 0 TO order-1 DO
      FOR j := 0 TO order-1 DO
        FOR k := 0 TO 1 DO
          cc[i,j,k] := cb[i,j,k+2];
        END
      END
    END;
    RETURN cc;
  END MakeColumnCube;

PROCEDURE MakeBigCube(order: CARDINAL) : REF TriPair =
  (*
    Build one tridimensional array of cubic cells with random geometry.
    Return two arrays of pairs of order: "order x order", not eliminated by
    gluing procedure of columns of cubes.
  *)
  VAR
    cd := NEW(REF ARRAY OF ARRAY OF ARRAY OF PAIRS, order, order,order);
    ce := NEW(REF TriPair, order,order);
  BEGIN
    FOR i := 0 TO order-1 DO
      WITH cc = MakeColumnCube(order) DO
        FOR j := 0 TO order-1 DO
          FOR k := 0 TO order-1 DO
            FOR l := 0 TO 1 DO
              cd[i,j,k,l] := cc[j,k,l];
            END
          END
        END
      END
    END;

    (* gluing *)
    FOR i := 0 TO order-2 DO
      FOR k := 0 TO order-1 DO
        FOR j := 0 TO order-1 DO
          EVAL GlueCube(cd[i,k,j,1],Clock(cd[i+1,k,j,0]));
        END
      END
    END;

    (* Selecting the pairs for return *)
    FOR j := 0 TO order-1 DO
      FOR k := 0 TO order-1 DO
        ce[j,k,0] := cd[0,k,j,0];
        ce[j,k,1] := cd[order-1,k,j,1];
      END
    END;
    RETURN ce;
  END MakeBigCube;


PROCEDURE SetCubeProperties(
  a: REF TriPair;
  n: CARDINAL;
  READONLY tp: Topology;
  ) =
  PROCEDURE SetVertex(v: Vertex) =
   (* Emphasizes the original vertices. *) 
   BEGIN
     v.exists := TRUE;
   END SetVertex;

  PROCEDURE SetEdge(a: Pair) =
   (* Emphasizes the original edges. *) 
   BEGIN
     WITH
       e = NARROW(a.facetedge.edge, Triangulation.Edge)
     DO	
       e.exists := TRUE;
     END;	
   END SetEdge;

  PROCEDURE SetFace(a: Pair) =
   (* Emphasizes the original edges. *) 
   BEGIN
     WITH
       f = NARROW(a.facetedge.face, Triangulation.Face)
     DO	
       f.exists := TRUE;
     END;	
   END SetFace;

  PROCEDURE FaceOnBoundary(a: Pair) : BOOLEAN =
   (* Return "TRUE" if the face "a.facetedge.face" is on the
      frontier.
   *)
   BEGIN
     IF Ppos(a) = NIL OR Pneg(a) = NIL
       THEN RETURN TRUE;
       ELSE RETURN FALSE;
     END	
   END FaceOnBoundary;
     
  VAR
    c1,c2,c3,c4 : REF ARRAY OF Pair;         (* edges on the columns *)
    index       : CARDINAL := 0;
  BEGIN
     c1 := NEW(REF ARRAY OF Pair, n);
     c2 := NEW(REF ARRAY OF Pair, n);
     c3 := NEW(REF ARRAY OF Pair, n);
     c4 := NEW(REF ARRAY OF Pair, n);

    (* Set the original vertices *)

    SetVertex(Org(Clock(Enext(a[0,0,0]))));  (* vertex 1 *)   
    SetVertex(Org(Clock(a[n-1,0,0])));       (* vertex 2 *)
    SetVertex(Org(a[n-1,n-1,0]));            (* vertex 3 *)
    SetVertex(Org(a[n-1,0,1]));              (* vertex 4 *)
    SetVertex(Org(Clock(a[n-1,n-1,1])));     (* vertex 5 *)
    SetVertex(Org(Clock(Enext(a[0,n-1,1]))));(* vertex 6 *)
    SetVertex(Org(Enext_1(a[0,0,1])));       (* vertex 7 *)
    SetVertex(Org(Enext_1(a[0,n-1,0])));     (* vertex 8 *)

    (* Set the original edges *)
    FOR i := 0 TO n-1 DO
      SetEdge(Enext_1(a[i,0,1]));
      SetEdge(Enext(a[i,0,0]));
      SetEdge(a[n-1,i,0]);
      SetEdge(a[n-1,i,1]);
      SetEdge(Enext_1(a[i,n-1,0]));
      SetEdge(Enext(a[i,n-1,1]));
      SetEdge(Enext(Enext(a[0,i,1])));
      SetEdge(Enext(Enext(a[0,i,0])));
    END;

    (* computing edges on the colums *)
    c1[0] := Enext(Fnext(a[n-1,0,0]));
    c2[0] := Enext_1(Fnext(a[n-1,n-1,0]));
    c3[0] := Enext(Fnext(Enext(a[0,0,0])));
    c4[0] := Enext_1(Fnext(Enext_1(a[0,n-1,0])));

    FOR k := 1 TO n-1 DO
      c1[k] := Clock(Enext_1(Fnext_1(Enext(c1[k-1]))));
      c2[k] := Clock(Enext(Fnext_1(Enext_1(c2[k-1]))));
      c3[k] := Clock(Enext_1(Fnext_1(Enext(c3[k-1]))));
      c4[k] := Clock(Enext(Fnext_1(Enext_1(c4[k-1]))));
    END;

    (* Now set the the original edges *)
    FOR j := 0 TO n-1 DO
      SetEdge(c1[j]);
      SetEdge(c2[j]);
      SetEdge(c3[j]);
      SetEdge(c4[j]);
    END;

    FOR i := 0 TO tp.NF-1 DO
      WITH
	f = tp.face[i].pa
      DO  
	IF FaceOnBoundary(f) THEN SetFace(f); INC(index) END;
      END
    END;
    <* ASSERT index = 6*n*n *>
  END SetCubeProperties;  

BEGIN
END Squared.

(**************************************************************************)
(*                                                                        *)
(* Copyright (C) 2000 Universidade Estadual de Campinas (UNICAMP)         *)
(*                                                                        *)
(* Authors:                                                               *)
(*   L. Lozada & J. Stolfi      - UNICAMP <lozada@dcc.unicamp.br>         *)
(*                                                                        *)
(* This file can be freely used, distributed, and modified, provided that *)
(* this copyright and authorship notice is included in every copy or      *)
(* derived version.                                                       *)
(*                                                                        *)
(* DISCLAIMER: This software is offered ``as is'', without any guarantee  *)
(*   as to fitness for any particular purpose.  Neither the copyright     *)
(*   holder nor the authors or their employers can be held responsible    *)
(*   for any damages that may result from its use.                        *)
(*                                                                        *)
(**************************************************************************)