MODULE Triangulation;

(*
   This module contain essentially procedures created by R. Marcone and
   J. Stolfi (see the copyright and authorship futher down), modified
   extensively by L. Lozada for the visualization of 3D-maps.

   Revisions:
 
     03-08-2000 : Optimized version of the Read and Write procedures by
                  J. Stolfi.

     30-08-2000 : Modified the Read and Write procedures for include the case
                  when the cells are octahedra.
 
     07-10-2000 : Added procedure for compute the barycenter of a tetrahedron.

     27-01-2001 : Modified for exploding cubic cells.

*)


IMPORT Octf, Random, LR4, LR4Extras, Stdio, Wr, Fmt, Thread, Math, FloatMode,
       FileRd, FileWr, Mis, Text, OSError, Lex, Rd, R3, FileFmt; 

FROM Octf IMPORT Fnext, Clock, SrotBits, RBits, SpliceFacets, Enext_1, Srot,
     Spin, SetFace, GetPairNum,Tors, Enext, Fnext_1, Onext, SetFnext, SetEdge,
     SetEdgeAll, SetEnext;

FROM Stdio IMPORT stderr;
FROM Mis IMPORT WriteCommentsJS;
FROM FileFmt IMPORT WriteHeader, WriteFooter, ReadHeader, ReadFooter;

<* FATAL Thread.Alerted, Rd.Failure, Wr.Failure *> 
<* FATAL Rd.EndOfFile, FloatMode.Trap, Lex.Error, OSError.E *>

REVEAL
  FacetEdge = PublicFacetEdge BRANDED OBJECT
    org: ARRAY RBits OF Node;

  OVERRIDES
    init := FacetEdgeInit;
  END;
 
REVEAL     
  Vertex = PublicVertex BRANDED OBJECT 
    END;

REVEAL 
  Polyhedron = PublicPolyhedron BRANDED OBJECT
    END;

(* === INIT METHODS === *)

PROCEDURE FacetEdgeInit(fe: FacetEdge): FacetEdge =
  BEGIN
    EVAL NARROW(fe, Octf.FacetEdge).init(); 
    fe.org[0] := NIL; (* For now, vertex of primal: C  *)
    fe.org[1] := NIL; (* For now, vertex of dual:   C' *)
    fe.org[2] := NIL; (* For now, vertex of primal: C  *)
    fe.org[3] := NIL; (* For now, vertex of dual:   C' *)
    
    RETURN  fe;
  END FacetEdgeInit;

(* === ELEMENT CREATION === *)
  
PROCEDURE MakeFacetEdge(): Pair =
  VAR a : Pair;
  BEGIN
    WITH e = NEW(FacetEdge).init() DO
      a := Pair{facetedge := e, bits := 0};
      a.facetedge.edge.pa := a;
      a.facetedge.face.pa := a;
    RETURN a;
    END;
  END MakeFacetEdge;

PROCEDURE MakeVertex(): Vertex =
  BEGIN
    RETURN NEW(Vertex)
  END MakeVertex;
  
PROCEDURE MakePolyhedron(): Polyhedron =
  BEGIN
    RETURN NEW(Polyhedron)
  END MakePolyhedron;

(* === PAIR PROPERTIES === *)

PROCEDURE Org(a: Pair): Node =
  BEGIN
    WITH c = NARROW(a.facetedge, FacetEdge).org[SrotBits(a)] DO
      RETURN c;
    END;
  END Org;

PROCEDURE SetOrg(a: Pair; n: Node) =
  BEGIN
    WITH c = NARROW(a.facetedge, FacetEdge).org[SrotBits(a)] DO
      c := n;
    END;
  END SetOrg;

PROCEDURE Set(a: Pair; n: Node) =
  VAR c : Pair := a;
  BEGIN
    WITH
      nei = Octf.NumberEdgesForDegree(ARRAY OF Pair{a})
    DO
      FOR i := 0 TO LAST(nei^) DO
        WITH
	  b = nei[i] 
        DO
	  c := b;  
	  REPEAT
            SetOrg(c,n);
	    c := Fnext(c)    
          UNTIL (c = b);
	END
      END
    END
  END Set;

PROCEDURE SetAllOrgs(a: Pair; n: Node) =
  BEGIN
    Set(a,n);
  END SetAllOrgs;

PROCEDURE Pneg(a: Pair): Node =
  BEGIN
    RETURN Org(Srot(a))
  END Pneg;

PROCEDURE SetPneg(a: Pair; n: Node) =
  BEGIN
    SetOrg(Srot(a), n)
  END SetPneg;

PROCEDURE SetNextPneg(a: Pair; n: Node) =
  VAR t: Pair := a;
  BEGIN
    REPEAT
      SetPneg(t, n);
      t := Enext_1(t);
    UNTIL t = a;
  END SetNextPneg;

PROCEDURE SetAllPneg(a : Pair; n: Node) =
  VAR t: Pair := a;
  BEGIN
    SetNextPneg(t,n);
    REPEAT
      SetNextPneg(Clock(Fnext_1(t)),n);
      t := Enext_1(t); 
    UNTIL t = a; 
  END SetAllPneg;

PROCEDURE Ppos(a: Pair): Node =
  BEGIN
    RETURN Pneg(Clock(a))
  END Ppos;

PROCEDURE SetPpos(a: Pair; n: Node) =
  BEGIN
    SetOrg(Tors(a), n)
  END SetPpos;

PROCEDURE SetAllPpos(a: Pair; n: Node) =
  VAR t: Pair := a;
  BEGIN
    REPEAT  
      SetPpos(t, n);
      t := Enext_1(t);
    UNTIL t = a;
  END SetAllPpos;

PROCEDURE SetNextPpos(a : Pair; n: Node) =
  VAR t : Pair := Clock(a);
  BEGIN
    SetAllPneg(t,n);
    REPEAT
      SetAllPneg(Fnext_1(t),n);
      t := Enext_1(t); 
    UNTIL t = Clock(a); 
  END SetNextPpos;

PROCEDURE OrgV(a: Pair): Vertex =
  BEGIN
    RETURN NARROW(Org(a), Vertex);
  END OrgV;

PROCEDURE DesV(a: Pair): Vertex =
  BEGIN
    RETURN OrgV(Clock(a));
  END DesV;

PROCEDURE PnegP(a: Pair): Polyhedron =
  BEGIN
    RETURN NARROW(Pneg(a), Polyhedron);
  END PnegP;

PROCEDURE PposP(a : Pair): Polyhedron = 
  BEGIN
    RETURN PnegP(Clock(a));
  END PposP;

PROCEDURE TetraNegVertices(a: Pair): ARRAY [0..3] 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(a)), 
      r = OrgV(Enext_1(a)),
      s = OrgV(Enext_1(Fnext_1(a)))
    DO
      <* ASSERT Pneg(a) = Ppos(Enext_1(Fnext_1(a))) *>
      <* ASSERT (p # q) AND (q # r) AND (r # s) *>
      RETURN ARRAY [0..3] OF Vertex{p,q,r,s}
    END
  END TetraNegVertices;

PROCEDURE TetraFaces(a: Pair): ARRAY [0..3] OF Face =
  BEGIN
    <* ASSERT Pneg(a) # NIL *>
    WITH
      f0 = a.facetedge.face,
      f1 = Fnext_1(a).facetedge.face,
      f2 = Fnext_1(Enext(a)).facetedge.face,
      f3 = Fnext_1(Enext_1(a)).facetedge.face
    DO
      <* ASSERT (f0 # f1) AND (f1 # f2) AND (f2 # f3) *> 
      RETURN ARRAY [0..3] OF Face{f0,f1,f2,f3}
    END
  END TetraFaces;

PROCEDURE TetraEdges(a: Pair): ARRAY [0..5] OF Edge =
  BEGIN
    <* ASSERT Pneg(a) # NIL *>
    WITH
      e0 = a.facetedge.edge,
      e1 = Enext(a).facetedge.edge,
      e2 = Enext_1(a).facetedge.edge,
      e3 = Enext(Fnext_1(a)).facetedge.edge,
      e4 = Enext_1(Fnext_1(a)).facetedge.edge,
      e5 = Enext(Fnext_1(Enext(a))).facetedge.edge
    DO
      <* ASSERT (e0#e1) AND (e1#e2) AND (e2#e3) AND (e3#e4) AND (e4#e5) *> 
      RETURN ARRAY [0..5] OF Edge{e0,e1,e2,e3,e4,e5}
    END
  END TetraEdges;

PROCEDURE FaceEdges(a: Pair): ARRAY [0..2] OF Edge =
  BEGIN
    WITH
      e0 = a.facetedge.edge,
      e1 = Enext(a).facetedge.edge,
      e2 = Enext_1(a).facetedge.edge
    DO
      <* ASSERT (e0 # e1) AND (e1 # e2) *> 
      RETURN ARRAY [0..2] OF Edge{e0,e1,e2}
    END
  END FaceEdges;

PROCEDURE EdgeIsBorder(a: Pair): BOOL = 
  VAR b: Pair := a;
  BEGIN
    REPEAT
      IF Pneg(b) = NIL THEN RETURN TRUE END;
      b := Fnext(b)
    UNTIL b = a;
    RETURN FALSE
  END EdgeIsBorder;

PROCEDURE FaceIsBorder(a: Pair): BOOL =
  BEGIN
    RETURN Ppos(a) = NIL OR Pneg(a) = NIL;
  END FaceIsBorder;

PROCEDURE TetraPosVertices(a: Pair): ARRAY [0..3] OF Vertex =
  (* Valid for the versions Spin(a), Clock(a) and SpinClock(a). *)
  BEGIN
    <* ASSERT Ppos(a) # NIL *>
    WITH
      p = OrgV(a),
      q = OrgV(Enext(a)), 
      r = OrgV(Enext_1(a)),
      s = OrgV(Enext_1(Fnext(a)))
    DO
      <* ASSERT Ppos(a) = Pneg(Enext_1(Fnext(a))) *>
      RETURN ARRAY [0..3] OF Vertex{p,q,r,s}
    END
  END TetraPosVertices;

PROCEDURE TetraNegPosVertices(a: Pair): ARRAY [0..4] OF Node =
  (* Valid for the versions Spin(a), Clock(a) and SpinClock(a). 
     We change the otion OrgV by Org such as, this procedure can be used
     in the dual space but the assertion is valid only in the primal space.
  *)
  BEGIN
    WITH
      p = Org(a),
      q = Org(Enext(a)), 
      r = Org(Enext_1(a)),
      s = Org(Enext_1(Fnext_1(a))),
      t = Org(Enext_1(Fnext(a)))
    DO
      <* ASSERT Ppos(a) = Pneg(Enext_1(Fnext(a))) AND
                Pneg(a) = Ppos(Enext_1(Fnext_1(a)))
      *>
      RETURN ARRAY [0..4] OF Node{p,q,r,s,t}
    END
  END TetraNegPosVertices;

(* ================= CONSTRUCTION TOOLS ==========  *)


PROCEDURE MakeTetraTopo(nx,ny : CARDINAL): ARRAY [0..7] OF Pair = 
  VAR
    FacetEdgeCount: CARDINAL := 0; 
    PolyhedronCount: CARDINAL := 0;

  PROCEDURE MakeTriangle(): Pair =
    (* Make one triangular face and set of 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 MakeCell(a: Pair): Pair =
    (*
      Build a new tetrahedral cell by insertion of two triangular faces and
      operations SpliceFacets. The argument "a" is the pair return by
      MakeTriangle(). The pair "f" is return by MakeCell().
    *)
    VAR c : Pair;
    BEGIN    
      c := Enext(Fnext(a));
      SetEdge(a, Enext_1(c).facetedge.edge);
      WITH f = MakeTriangle(), g = MakeTriangle() DO

        SetFnext(Enext_1(a), Enext_1(f));
	SetEdge(Enext_1(f), Enext_1(a).facetedge.edge);
        SetAllOrgs(Enext_1(a), Org(Enext_1(a)));
        SetAllOrgs(Clock(Enext_1(a)), Org(Clock(Enext_1(a))));

        SetFnext(Clock(f), Enext(c));
        SetEdge(Clock(f), Enext(c).facetedge.edge); 
        SetAllOrgs(Enext(c), Org(Enext(c)));
        SetAllOrgs(Clock(Enext(c)), Org(Clock(Enext(c))));     
        
        SetFnext(Enext(g), Enext(f));
	SetEdge(Enext(g), Enext(f).facetedge.edge);
        SetAllOrgs(Enext(f), Org(Enext(f)));
        SetAllOrgs(Clock(Enext(f)), Org(Clock(Enext(f))));
        
        SetFnext(g, c);
	SetEdge(g, c.facetedge.edge);
        SetAllOrgs(c, Org(c));
        SetAllOrgs(Clock(c), Org(Clock(c))); 
        
        SetFnext(Enext_1(g), Clock(Enext(a)));
	SetEdge(Enext_1(g), Clock(Enext(a)).facetedge.edge); 
        SetAllOrgs(Enext(a), Org(Enext(a)));
        SetAllOrgs(Clock(Enext(a)), Org(Clock(Enext(a))));
                 
        WITH
          p = MakePolyhedron()
        DO
          p.num := PolyhedronCount;
          INC(PolyhedronCount);
          SetAllPneg(f,p);
        END;
        RETURN f;
      END;
  END MakeCell;

  PROCEDURE MakeCellRow(a: Pair): Pair =
    (*
      Adds one new row of tetradedral cells and return the pair facetedge
      belong to the cell more rigth.
    *)
   BEGIN
     FOR col := 0 TO nx-1 DO
       a := MakeCell(a);
     END;
     RETURN a;
   END MakeCellRow;

  VAR
    t,b: Pair;
    ca: ARRAY [0..7] OF Pair;
  BEGIN
    (* =============== create bottom row of triangles =========== *)
    FOR col := 0 TO nx-1 DO
      WITH c = MakeTriangle() DO
        IF col = 0 THEN 
          t := c
        ELSE
          SpliceFacets(Enext(b), Clock(Enext_1(c)));
          SetEdge(Enext(b), Clock(Enext_1(c)).facetedge.edge);
          SetAllOrgs(Enext(b), Org(Enext(b))); 
          SetAllOrgs(Clock(Enext(b)), Org(Clock(Enext(b))));
        END;
        b := c;
      END;
    END;

    ca[2] := t;  ca[7] := Spin(Clock(b));

    FOR row := 0 TO ny-1 DO
      WITH a = MakeTriangle() DO
        IF row = 0 THEN ca[1] := Clock(Enext_1(a)) END;           
        IF row = ny-1 THEN ca[4] := Spin(Enext_1(a)) END;

        SpliceFacets(Clock(a), Clock(Enext_1(t)));
        SetAllOrgs(Enext_1(t), Org(Enext_1(t)));
        SetAllOrgs(Clock(Enext_1(t)), Org(Clock(Enext_1(t))));

        WITH  aa = MakeCellRow(a) DO      
          t := Fnext_1(t);
          SetOrg(Enext_1(t), Org(Enext_1(a)));
          b := Fnext_1(b);
          SetOrg(Clock(Enext(b)), Org(Enext_1(aa)));
          IF row = 0 THEN ca[0] := Clock(Spin(Enext_1(aa))) END; 
          IF row = ny-1 THEN ca[5] := Enext_1(aa) END;
        END;
      END;
    END;
    ca[3] := Spin(t);
    ca[6] := Clock(b);
    RETURN ca;
  END MakeTetraTopo;         

PROCEDURE EmphasizeTetrahedron(a, b: Pair; n: CARDINAL) =

   PROCEDURE HiddenVertex(v: Vertex) =
    (* Hidden the vertex "v". *) 
    BEGIN
      v.exists := FALSE;
    END HiddenVertex;

   PROCEDURE HiddenEdge(a: Pair) =
    (* Hidden the edge "a.facetedge.edge".*) 
    BEGIN
      WITH
        e = NARROW(a.facetedge.edge, Edge)
      DO	
        e.exists := FALSE;
      END;	
    END HiddenEdge;

   PROCEDURE HiddenFace(a: Pair) =
    (* Hidden the face "a.facetedge.face".*) 
    BEGIN
      WITH
        f = NARROW(a.facetedge.face, Face)
      DO	
        f.exists := FALSE;
      END;	
    END HiddenFace;

   PROCEDURE HiddenRingFace(a: Pair) =
    (*  Hidden the ring face "a.facetedge.face". *)
    VAR an : Pair := Fnext_1(a);
    BEGIN
      FOR j := 1 TO n-1 DO
        HiddenFace(an);
        an := Fnext_1(an);
      END;
    END HiddenRingFace;

  VAR
    ta,tb: ARRAY [0..100] OF Pair;
  BEGIN
    
    ta[0] := a; tb[0] := b;

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

   
    FOR i := 0 TO n-1 DO
      HiddenRingFace(ta[i]);
      HiddenRingFace(tb[i]);
    END;

    FOR i := 1 TO n-1 DO
      HiddenVertex(OrgV(tb[i])); 
      HiddenVertex(OrgV(ta[i]));
      VAR dn: Pair := tb[i];
      BEGIN
        FOR j := 0 TO n DO
          HiddenEdge(Enext_1(dn));
          dn := Fnext_1(dn);
        END
      END
    END;

    FOR i := 0 TO n-2 DO
      HiddenEdge(Enext(ta[i]));
      HiddenEdge(Enext(Fnext(ta[i])));
    END;
  END EmphasizeTetrahedron;  

PROCEDURE Glue(
  a,b   : Pair;
    n   : CARDINAL;
  setorg: BOOL := TRUE;
  ): Pair =
  (* The pair "a" and "b" have the same Orientation and Spin bits,  such as,
     after of the glue procedure performs: Pneg(a) = Pneg(b).
  *)
  VAR
    ta,tb: ARRAY [0..100] OF Pair;
  BEGIN
    (*Wr.PutText(Stdio.stderr, "colou\n");*)
    (* sanity check *)
    <* ASSERT n >= 1 *>
    ta[0] := a; tb[0] := b;

    IF n > 1 THEN
      FOR i := 1 TO n-1 DO
        ta[i] := Clock(Enext_1(Fnext_1(Enext(ta[i-1]))));
        tb[i] := Clock(Enext_1(Fnext(Enext(tb[i-1]))));
        <* ASSERT ta[i] # a *>
        <* ASSERT tb[i] # b *>  
      END;
    END;
      
    Octf.Meld(b, a);
    
    (* updating edges relations for i=0 *)
    
    SetEdgeAll(a,          a.facetedge.edge);
    SetEdgeAll(Enext(a),   Enext(a).facetedge.edge);
    SetEdgeAll(Enext_1(a), Enext_1(a).facetedge.edge);
    
    IF setorg THEN
      (* updating vertices relations for i=0 *)
      SetAllOrgs(a,        Org(a));
      SetAllOrgs(Clock(a), Org(Clock(a)));
      SetAllOrgs(Enext(a), Org(Enext(a)));
      SetAllOrgs(Clock(Enext(a)), Org(Clock(Enext(a))));
      SetAllOrgs(Enext_1(a), Org(Enext_1(a)));
      SetAllOrgs(Clock(Enext_1(a)), Org(Clock(Enext_1(a))));
    END;
    
    (* updating polyhedron relations for i=0 *)
    SetPneg(a, Pneg(b));
    SetPneg(Enext_1(a), Pneg(Enext_1(b)));
    SetPneg(Enext(a), Pneg(Enext(b)));
    FOR i := 1 TO n-1 DO
      Octf.Meld(tb[i],ta[i]);
      
      (* updating edges relations *)
      SetEdgeAll(ta[i], ta[i].facetedge.edge);
      SetEdgeAll(Enext(ta[i]), Enext(ta[i]).facetedge.edge);
      SetEdgeAll(Enext_1(ta[i]), Enext_1(ta[i]).facetedge.edge);

      IF setorg THEN
        (* updating vertices relations *)
        SetAllOrgs(ta[i], Org(ta[i]));
        SetAllOrgs(Clock(ta[i]), Org(Clock(ta[i])));
        SetAllOrgs(Enext(ta[i]), Org(Enext(ta[i])));
        SetAllOrgs(Clock(Enext(ta[i])), Org(Clock(Enext(ta[i]))));
        SetAllOrgs(Enext_1(ta[i]), Org(Enext_1(ta[i])));
        SetAllOrgs(Clock(Enext_1(ta[i])), Org(Clock(Enext_1(ta[i]))));
      END;

      (* updating polyhedron relations *)
      WITH
        f = Pneg(tb[i]),
        g = Pneg(Enext_1(tb[i])),
        h = Pneg(Enext(tb[i]))
      DO
        SetPneg(ta[i],f);
        SetPneg(Enext_1(ta[i]),g);
        SetPneg(Enext(ta[i]),h);
      END
    END;

    IF setorg THEN
      SetAllOrgs(Clock(Enext_1(Fnext_1(Enext_1(ta[n-1])))),Org(Enext_1(a)));
    END;
    RETURN ta[n-1];
  END Glue;

(* === GLOBAL PROCEDURES === *)

VAR
  seenNode: REF ARRAY OF Node := NIL;
  nodeOnum: REF ARRAY OF CARDINAL := NIL;
  seenNodeCount: CARDINAL := 0;

PROCEDURE MarkNode(n: Node) =
  BEGIN
    IF seenNode = NIL OR NUMBER(seenNode^) <= seenNodeCount THEN
      DoubleseenNode()
    END;
    WITH k = seenNodeCount+0 DO
      seenNode[k] := n;
      nodeOnum[k] := n.num;
      n.num := k;
      INC(seenNodeCount)
    END
  END MarkNode;

PROCEDURE DoubleseenNode() =
  VAR sz: CARDINAL;
  BEGIN
    IF seenNode = NIL THEN sz := 0 ELSE sz := NUMBER(seenNode^) END;
    WITH 
      szNew = MAX(2*sz, 1000),
      seenNodeNew  = NEW(REF ARRAY OF Node, szNew),
      nodeOnumNew  = NEW(REF ARRAY OF CARDINAL, szNew)
    DO
      IF seenNode # NIL THEN
        SUBARRAY(seenNodeNew^, 0, sz) := seenNode^;
        SUBARRAY(nodeOnumNew^, 0, sz) := nodeOnum^;
      END;
      seenNode := seenNodeNew;
      nodeOnum := nodeOnumNew;
    END
  END DoubleseenNode;

PROCEDURE NodeIsMarked(n: Node): BOOL  =
  BEGIN
    RETURN (n.num < seenNodeCount) AND (seenNode[n.num] = n)
  END NodeIsMarked;

PROCEDURE UnmarkMarkedNodes() =
  BEGIN
    WHILE seenNodeCount > 0 DO
      WITH k = seenNodeCount-1, n = seenNode[k] DO
        n.num := nodeOnum[k];
        DEC(seenNodeCount)
      END
    END
  END UnmarkMarkedNodes;

PROCEDURE EnumVertices(a: Pair; visit: VisitProc) =
  CONST IniStackSize = 1000;
  VAR
    festack := NEW(REF ARRAY OF FacetEdge, IniStackSize);
    bstack := NEW(REF ARRAY OF SRBits, IniStackSize);
    top : CARDINAL;    (* top for "festack" stack *)

  PROCEDURE DoubleStack() =
    BEGIN
      WITH 
        sz = NUMBER(festack^), szNew = 2*sz,
        festackNew = NEW(REF ARRAY OF FacetEdge, szNew),
        bstackNew  = NEW(REF ARRAY OF SRBits, szNew)
      DO
        SUBARRAY(festackNew^, 0, sz) := festack^; festack := festackNew;
        SUBARRAY(bstackNew^, 0, sz) := bstack^; bstack := bstackNew;
      END
    END DoubleStack;
    
  PROCEDURE Stack(c: Pair) =
    VAR cn : Pair := c;
        dn : Pair;
    BEGIN
      REPEAT
        IF top > LAST(festack^) THEN DoubleStack() END;
        festack[top] := cn.facetedge;
        bstack[top] := cn.bits; 
        top := top + 1;
        dn := Clock(Enext_1(cn));
        REPEAT
	  IF top > LAST(festack^) THEN DoubleStack() END;
          festack[top] := dn.facetedge; bstack[top] := dn.bits;
          top := top + 1;
          dn := Fnext(dn);
        UNTIL( dn = Clock(Enext_1(cn)) );
        cn := Fnext(cn);
      UNTIL (cn = c);
    END Stack;

  PROCEDURE VisitAndMark(c: Pair)=
    (* If org(c) is diferent that the origins pairs stacked, stack the pair
       "c".
    *)
    VAR cn : Pair;
    BEGIN
      WITH n = Org(c) DO
        IF n # NIL AND NOT NodeIsMarked(n) THEN
          visit(c);
          MarkNode(n);
          Stack(c);
          cn := c;
          REPEAT
            Stack(cn);
            cn := Onext(cn);
          UNTIL (cn = c)
        END
      END;        
    END VisitAndMark;

  VAR seen: CARDINAL;
  BEGIN
    top := 0;
    seen := 0;
    <* ASSERT seenNodeCount = 0 *>
    VisitAndMark(a);
    WHILE seen < top  DO
      WITH b = Pair{festack[seen], bstack[seen]} DO
        VisitAndMark(Fnext(b));
        VisitAndMark(Fnext(Enext(b)));
        VisitAndMark(Enext(b));
        VisitAndMark(Fnext_1(Enext_1(b)));
      END;
      seen := seen + 1;
    END;
    UnmarkMarkedNodes();
  END EnumVertices;


PROCEDURE NumberVertices(a: Pair): CARDINAL = 
  VAR n: CARDINAL := 0;
  
  PROCEDURE Visit(c: Pair) =
  
    PROCEDURE Visit_(c: Pair) =
      VAR cn,dn: Pair;  
      BEGIN
        WITH v = Org(c) DO
          cn := c;
          REPEAT
            WITH vn = Org(cn) DO
              <* ASSERT vn = v *>
              vn.num := n;
              dn := Clock(Enext_1(cn));
              REPEAT
                WITH wn = Org(dn) DO
                  <* ASSERT wn = vn *>
                  wn.num := n;
                  dn := Fnext(dn);
                END;
              UNTIL ( dn = Clock(Enext_1(cn)) );
              cn := Fnext(cn);
            END;
          UNTIL ( cn = c );
        END;
      END Visit_;

    VAR cn: Pair := c;
    BEGIN
      REPEAT        
        Visit_(cn);
        cn := Onext(cn);
      UNTIL cn = c;
      INC(n);
   END Visit;

   PROCEDURE VisitDual(c: Pair) =
     PROCEDURE VisitDual_(c: Pair) =
       VAR
         cn: Pair;
	 p,pn: Node;
       BEGIN
         p := Org(c);
	 IF p # NIL THEN
           cn := c;
           REPEAT
             pn := Org(cn);
	     IF pn # NIL THEN
               <* ASSERT pn = p *> 
               pn.num := n;
             END;
             cn := Fnext(cn);
           UNTIL( cn = c )
         END;
       END VisitDual_;

     VAR cn: Pair := c;
     BEGIN
       VisitDual_(cn);
       VisitDual_(Clock(Enext_1(cn)));
       INC(n);
     END VisitDual;

    BEGIN
      IF Octf.DualBit(a) = 0 THEN
        EnumVertices(a, Visit);
      ELSE
        EnumVertices(a, VisitDual);
      END;
    RETURN n    
  END NumberVertices;

PROCEDURE MakeTopology(a: Pair): Topology =
  VAR
    top   : Topology;
    euler : INTEGER;
  BEGIN
    top.NV         := NumberVertices(a);
    top.vertex     := NEW(REF Vertices, top.NV); 
    top.facetedge  := Octf.NumberFacetEdges(ARRAY OF Pair{a});  
    top.NFE        := NUMBER(top.facetedge^);
    top.edge       := Octf.NumberEdges(ARRAY OF Pair{top.facetedge[0]});
    top.face       := Octf.NumberFacets(ARRAY OF Pair{top.facetedge[0]});
    top.out        := NEW(REF ARRAY OF Pair, top.NV);
    top.NE         := NUMBER(top.edge^);
    top.NF         := NUMBER(top.face^);
    top.NP         := NumberVertices(Srot(a));
    top.polyhedron := NEW(REF ARRAY OF Polyhedron, top.NP);
    top.region     := NEW(REF ARRAY OF Pair, top.NP);
    (* top.der        := Octf.DegreeEdgeRing(top.face[0].pa); *)
    (* top.bdr        := bdr; *)

    FOR i:= 0 TO top.NFE-1 DO
      VAR
        c: Pair := top.facetedge[i];
        v,p: Node;
        vi,pi: CARDINAL;
      BEGIN
        FOR k:= 0 TO 1 DO
          v := Org(c); 
          vi := v.num;
          top.vertex[vi] := v;
          top.out[vi] := c;
          p := Pneg(c);
          IF p # NIL THEN
            pi := p.num;
            top.polyhedron[pi] := p;
            top.region[pi] := Srot(c);
            <* ASSERT Pneg(c) = Org(top.region[pi]) *>
          END;
          c := Clock(c); 
        END
      END
    END;
    
    WITH
      m = Mis.NumDigits(top.NFE)+1 
    DO
      Wr.PutText(stderr, "\n");
      Wr.PutText(stderr, "nv  := " & Fmt.Pad(Fmt.Int(top.NV),m)  & "\n"); 
      Wr.PutText(stderr, "ne  := " & Fmt.Pad(Fmt.Int(top.NE),m)  & "\n");
      Wr.PutText(stderr, "nf  := " & Fmt.Pad(Fmt.Int(top.NF),m)  & "\n");
      Wr.PutText(stderr, "np  := " & Fmt.Pad(Fmt.Int(top.NP),m)  & "\n");
      Wr.PutText(stderr, "nfe := " & Fmt.Pad(Fmt.Int(top.NFE),m) & "\n");
    END;

    (* Euler's Number *)
    euler := top.NV-top.NE+top.NF-top.NP;
    IF euler = 0 THEN
      Wr.PutText(stderr,"\nThe Map performs the Euler's Number = 0\n");
    ELSE
      Wr.PutText(stderr,"\nThe Euler's Number: " & Fmt.Int(euler) & "\n");
    END;
    RETURN top;
  END MakeTopology;

PROCEDURE GetTetraCorners(a : Pair): TetraCorners =
  VAR c : TetraCorners;
  BEGIN
    WITH
      a0 = a,
      a1 = Enext(a0),
      a2 = Enext(a1),
      b0 = Clock(Fnext_1(a0)),
      c0 = Clock(Fnext_1(a1)),
      d0 = Clock(Fnext_1(a2))      
    DO
      (* the 12 elements on the array must have the same negative
         polyhedron.
      *)
      c[0 ] := a0;
      c[1 ] := Enext(a0);
      c[2 ] := Enext_1(a0);
      c[3 ] := b0;
      c[4 ] := Enext(b0);
      c[5 ] := Enext_1(b0);
      c[6 ] := c0;
      c[7 ] := Enext(c0);
      c[8 ] := Enext_1(c0);
      c[9 ] := d0;
      c[10] := Enext(d0);
      c[11] := Enext_1(d0);
      <* ASSERT c[0] = a *>
    END;
    RETURN c;
  END GetTetraCorners;
  
PROCEDURE CollectTetrahedra(READONLY tp: Topology): REF ARRAY OF Pair =
  VAR seen,processed,conflicts: CARDINAL := 0;
  BEGIN
    WITH
      rt = NEW(REF ARRAY OF Pair, tp.NP),
       t = rt^ ,
       q = NEW(REF ARRAY OF Pair, tp.NP)
    DO
      FOR i := 0 TO LAST(t) DO
        t[i] := Pair{facetedge := NIL, bits := 0};
      END;
      FOR k := 0 TO LAST(t) DO
        IF t[k].facetedge = NIL THEN
	  t[k] := Tors(tp.region[k]);
          q[seen] := t[k];
          INC(seen);	   
          WHILE processed < seen DO
            WITH
              a = q[processed]
            DO
              PROCEDURE Visit(f: Pair) =
          	BEGIN
		  IF Pneg(f) = NIL THEN RETURN END;
          	  <* ASSERT Consistent(a,f) *>
          	  WITH
          	    n = Pneg(f).num
          	  DO
          	    IF t[n].facetedge = NIL THEN
          	      t[n] := f;
          	      q[seen] := f;
          	      INC(seen);
          	    ELSE
          	      <* ASSERT Pneg(t[n]).num = n *>
          	      IF NOT SameOrientation(t[n],f) THEN INC(conflicts) END
          	    END    
          	  END
          	END Visit;
          	
              PROCEDURE Consistent(p,q:Pair): BOOL =
          	(* Checks whether Pneg(q) and the Pneg(q) are adjacentes tetrahedra
          	   with consistent orientations.
          	*)
          	BEGIN
          	  WITH
          	    cp = GetTetraCorners(p),
          	    cq = GetTetraCorners(q)                
          	  DO
          	    FOR i := 0 TO 11 DO
          	      FOR j := 0 TO 11 DO
          		IF cp[i] = Clock(cq[j]) THEN RETURN TRUE END;
          	      END;
          	    END;
          	    RETURN FALSE;
          	  END
          	END Consistent;
          
              PROCEDURE SameOrientation(p,q:Pair): BOOL =
          	(* Checks whether Pneg(p) = Pneg(q) and both have the same
          	   orientation.
          	*)
          	BEGIN
          	  WITH
          	    cp = GetTetraCorners(p)               
          	  DO
          	    FOR i := 0 TO 11 DO
          	      IF cp[i] = q THEN RETURN TRUE END;
          	    END;
          	    RETURN FALSE;
          	  END
          	END SameOrientation;
          	
              BEGIN
          	WITH
          	  b = Fnext(a),
          	  c = Fnext_1(Enext(a)),
          	  d = Fnext_1(Enext_1(a)),
          	  e = Fnext_1(a)	  
          	DO
          	  Visit(b);
          	  Visit(c);
          	  Visit(d);
          	  Visit(e);
          	END
              END;
              INC(processed);    
            END
          END
	END  
      END;
      IF conflicts # 0 THEN
        Wr.PutText(Stdio.stderr, "CollectTetrahedra: conflicts = " & Fmt.Int(conflicts) & "\n");
      END;
      <* ASSERT processed = tp.NP *>
      RETURN rt;
    END;
  END CollectTetrahedra;
    
PROCEDURE DegreeOfVertex(a: Pair): CARDINAL =
  BEGIN
    WITH
      edge   = Octf.NumberEdgesForDegree(ARRAY OF Pair{a}),
      degree = NUMBER(edge^)
    DO
      RETURN degree;
    END;
  END DegreeOfVertex;

PROCEDURE MakeAdjacencyMatrix(READONLY top: Topology): REF AdjacencyMatrix =
  VAR m := NEW(REF ARRAY OF ARRAY OF BOOL, top.NV, top.NV);
  BEGIN
    FOR i := 0 TO top.NV-1 DO 
      WITH mi = m[i] DO 
        FOR j := 0 TO top.NV-1 DO 
          mi[j] := FALSE
        END
      END
    END;
    FOR ei := 0 TO top.NFE-1 DO
      WITH a = top.facetedge[ei] DO
        WITH i = Org(a).num, j = Org(Clock(a)).num DO 
          m[i,j] := TRUE;
          m[j,i] := TRUE;
        END
      END
    END;
    RETURN m;
  END MakeAdjacencyMatrix;

PROCEDURE TriviallyIsomorphic(READONLY ta, tb: Topology): BOOL =
  BEGIN
    IF ta.NV # tb.NV
      OR ta.NF # tb.NF 
      OR ta.NE # tb.NE
      OR ta.NP # tb.NP 
      OR ta.NFE # tb.NFE
    THEN RETURN FALSE END;
    WITH NFE = ta.NFE DO
      FOR i := 0 TO NFE-1 DO 
        VAR sa: Pair := ta.facetedge[i]; sb: Pair := tb.facetedge[i];
        BEGIN
          FOR r := 0 TO 3 DO
            WITH 
              na = Org(sa).num, nb = Org(sb).num
            DO
              IF na # nb THEN RETURN FALSE END
            END;
            WITH 
              za = GetPairNum(Fnext(sa)),
              zb = GetPairNum(Fnext(sb))
            DO
              IF za # zb THEN RETURN FALSE END
            END;
            sa := Srot(sa);
            sb := Srot(sb);
          END
        END
      END
    END;
    RETURN TRUE
  END TriviallyIsomorphic;

PROCEDURE CheckOutAndRegion(READONLY top: Topology) =
  BEGIN
    FOR i := 0 TO top.NV-1 DO 
      WITH
        v = top.vertex[i],
        e = top.out[i]
      DO
        <* ASSERT v.num = i *>
        <* ASSERT Org(e) = v *>
      END
    END;
    FOR i := 0 TO top.NP-1 DO
      WITH
        p = top.polyhedron[i],
        r = top.region[i]
      DO
        <* ASSERT p.num = i *>
        <* ASSERT Org(r) = p *>
      END 
    END;
  END CheckOutAndRegion;

PROCEDURE GetVariableVertices(READONLY top: Topology;
                              VAR vr: ARRAY OF BOOL) =
  BEGIN
    FOR i := 0 TO top.NV-1 DO
      WITH v = top.vertex[i] DO
        vr[i] := NOT v.fixed; 
      END;
    END;
  END GetVariableVertices;
  
(* === GEOMETRIC TOOLS === *)

PROCEDURE InitCoords(
    coins: Random.T; 
    VAR c: Coords; 
    radius: LONGREAL := 1.0d0
  ) =
  BEGIN
    FOR i := 0 TO LAST(c) DO
      c[i] := LR4.T{
        coins.longreal(-radius, radius),
        coins.longreal(-radius, radius),
        coins.longreal(-radius, radius),
        coins.longreal(-radius, radius)
      }
    END
  END InitCoords;

PROCEDURE GenCoords(READONLY t: Topology): REF Coords =
  BEGIN
    WITH
      coins = NEW(Random.Default).init(TRUE),
      r = NEW(REF Coords, t.NV),
      c = r^
    DO
      FOR i := 0 TO LAST(c) DO 
        c[i] := LR4.T{
                      coins.longreal(-1.0d0, +1.0d0),
                      coins.longreal(-1.0d0, +1.0d0),
                      coins.longreal(-1.0d0, +1.0d0),
                      coins.longreal(-1.0d0, +1.0d0)
                     }
      END;
      RETURN r
    END
  END GenCoords;

PROCEDURE Displace(READONLY top: Topology; d: LR4.T; VAR c: Coords) =

  BEGIN
    FOR i := 0 TO LAST(c) DO 
      IF top.vertex[i].exists THEN
        WITH vc = c[i] DO vc := LR4.Add(vc, d) END
      END
    END
  END Displace;

PROCEDURE Scale(READONLY top: Topology; s: LONGREAL; VAR c: Coords) =
  BEGIN
    FOR i := 0 TO LAST(c) DO 
      IF top.vertex[i].exists THEN
        WITH vc = c[i] DO vc := LR4.Scale(s, vc) END;
      END
    END
  END Scale;

PROCEDURE FaceCross(a: Pair; READONLY c: Coords): LR4.T =
  BEGIN
    IF NOT a.facetedge.face.exists THEN
      RETURN LR4.T{1.0d0, 0.0d0, 0.0d0, 0.0d0}
    ELSE
      WITH
        ov = OrgV(a), 
        dv = OrgV(Clock(a)),
        pv = OrgV(Enext_1(a)), 
        qv = OrgV(Enext_1(Fnext(a))), 
        rv = OrgV(Enext_1(Fnext_1(a)))
      DO
        WITH
          o = c[ov.num],
          d = c[dv.num],
          p = c[pv.num],
          q = c[qv.num],
          r = c[rv.num],
          n1 = LR4Extras.Cross(LR4.Sub(p,o), LR4.Sub(r,o), LR4.Sub(d,o)),
          n2 = LR4Extras.Cross(LR4.Sub(p,o), LR4.Sub(d,o), LR4.Sub(q,o))
        DO 
          RETURN LR4.Add(n1,n2);
        END
      END
    END
   END FaceCross;

PROCEDURE FaceNormal(a: Pair; READONLY c: Coords): LR4.T =
  BEGIN
    WITH
      n = FaceCross(a,c),
      s = LR4.Norm(n)
    DO
      IF s = 0.0d0 THEN 
        RETURN LR4.T{1.0d0, 0.0d0, 0.0d0, 0.0d0}
      ELSE
        RETURN LR4.Scale(1.0d0/s, n)
      END
    END
  END FaceNormal;

PROCEDURE PolyCross(a: Pair; READONLY c: Coords): LR4.T =
  BEGIN
    IF PnegP(a) = NIL OR NOT PnegP(a).exists THEN
      RETURN LR4.T{1.0d0, 0.0d0, 0.0d0, 0.0d0}
    ELSE
      WITH
        o = c[OrgV(a).num],
        p = c[OrgV(Clock(a)).num], 
        q = c[OrgV(Enext_1(a)).num],
        r = c[OrgV(Enext_1(Fnext_1(a))).num]
      DO
        RETURN LR4Extras.Cross(LR4.Sub(q,o), LR4.Sub(r,o), LR4.Sub(p,o));
      END
    END
   END PolyCross;
            
PROCEDURE PolyNormal(a: Pair; READONLY c: Coords): LR4.T =
  BEGIN
    WITH
      n = PolyCross(a, c),
      s = LR4.Norm(n)
    DO
      IF s = 0.0d0 THEN 
        RETURN LR4.T{1.0d0, 0.0d0, 0.0d0, 0.0d0}
      ELSE
        RETURN LR4.Scale(1.0d0/s, n)
      END
    END
  END PolyNormal;
  
PROCEDURE FaceBarycenter(a: Pair; READONLY c: Coords): LR4.T =
  VAR
    ao : Pair;
    n : CARDINAL := 0;
    sum := LR4.T{0.0d0, ..};
  BEGIN
    ao := a;
    REPEAT
      WITH
        aoc = c[OrgV(ao).num]
      DO	
        sum := LR4.Add(sum, aoc);
	INC(n);
        ao := Enext(ao);
      END;
    UNTIL (ao = a);
    IF n = 0 THEN
      RETURN sum
    ELSE
      RETURN LR4.Scale(1.0d0/FLOAT(n,LONGREAL), sum)
    END
  END FaceBarycenter;

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

PROCEDURE EdgeCross(a: Pair; READONLY c: Coords): LR4.T =
  VAR sum: LR4.T := LR4.T{0.0d0,..};
      ao : Pair;
  BEGIN
    WITH
      uv = OrgV(a),
      u = c[uv.num]
      DO
        IF NOT a.facetedge.edge.exists THEN
          RETURN LR4.T{1.0d0, 0.0d0, 0.0d0, 0.0d0}
        ELSE
          ao := a;
          REPEAT
            WITH
              an = Fnext_1(ao),
              dv = OrgV(Clock(ao)), 
              pv = OrgV(Enext_1(ao)), 
              rv = OrgV(Enext_1(an))
            DO
              WITH
                d = c[dv.num],
                p = c[pv.num],
                r = c[rv.num],
                n = LR4Extras.Cross(LR4.Sub(p,u), LR4.Sub(d,u), LR4.Sub(r,u))
              DO
                IF ao = a THEN sum := n; ELSE sum := LR4.Add(sum,n) END;
              END;
              ao := an;
            END;
          UNTIL (ao = a);
          RETURN sum;
        END
    END;
  END EdgeCross;

PROCEDURE EdgeNormal(a: Pair; READONLY c: Coords): LR4.T =
  BEGIN
    WITH
      n = EdgeCross(a, c),
      s = LR4.Norm(n)
    DO
      IF s = 0.0d0 THEN 
        RETURN LR4.T{1.0d0, 0.0d0, 0.0d0, 0.0d0}
      ELSE
        RETURN LR4.Scale(1.0d0/s, n)
      END
    END
  END EdgeNormal;

PROCEDURE VertexCross(
    a: Pair;
    READONLY c: Coords;
    READONLY top : Topology;
  ): LR4.T =
  VAR 
    sum,n: LR4.T := LR4.T{0.0d0, ..};
  BEGIN
    WITH
      uv = OrgV(a), poly = StarOfVertex(a,top)
    DO
      IF NOT uv.exists THEN
        RETURN LR4.T{1.0d0, 0.0d0, 0.0d0, 0.0d0}
      ELSE
        FOR i := 0 TO LAST(poly^) DO
          IF poly[i][0].exists AND poly[i][1].exists AND
             poly[i][2].exists AND poly[i][3].exists
          THEN
            WITH
              u = c[poly[i][0].num],
              d = c[poly[i][1].num],
              p = c[poly[i][2].num],
              r = c[poly[i][3].num]
            DO
              n:= LR4Extras.Cross(LR4.Sub(p,u), LR4.Sub(d,u), LR4.Sub(r,u));
              IF i = 0 THEN sum := n; ELSE sum := LR4.Add(sum,n) END; 
            END;
          END;
        END;
      END;
      RETURN sum;
    END;
  END VertexCross;

PROCEDURE VertexNormal(a: Pair; READONLY c: Coords;
                       READONLY top: Topology): LR4.T =
  BEGIN
    WITH
      n = VertexCross(a, c, top),
      s = LR4.Norm(n)
    DO
     IF s = 0.0d0 THEN 
        RETURN LR4.T{1.0d0, 0.0d0, 0.0d0, 0.0d0}
      ELSE
        RETURN LR4.Scale(1.0d0/s, n)
      END
    END
  END VertexNormal;

PROCEDURE Neighbors(
    v: Vertex;
    READONLY top: Topology;
   ): REF Vertices =
  VAR
    rv: REF Vertices := NIL;
    nv: CARDINAL := 0;
    
  PROCEDURE Stack(u: Vertex) =
    BEGIN
      IF rv = NIL OR nv >= NUMBER(rv^) THEN
        WITH sv = NEW(REF Vertices, MAX(10, 2*nv)) DO
          IF rv # NIL THEN SUBARRAY(sv^, 0, nv) := rv^ END;
          rv := sv
        END
      END;
      rv[nv] := u; INC(nv);
    END Stack;
    
  PROCEDURE Check(b: Pair) =
    BEGIN
      IF OrgV(b) = v THEN Stack(OrgV(Clock(b))) END
    END Check;

  BEGIN
    FOR i := 0 TO top.NE - 1 DO
      WITH ei = top.edge[i].pa DO
        Check(ei); Check(Clock(ei))
      END
    END;
      
    WITH r = NEW(REF Vertices, nv) DO
      r^ := SUBARRAY(rv^, 0, nv);
      RETURN r;
    END;
  END Neighbors;
  
PROCEDURE StarOfVertex(a: Pair; READONLY top: Topology): REF ARRAY OF Quadv =

  PROCEDURE VertexPoly(READONLY top : Topology):  REF ARRAY OF Quadp =
     (* Return four pairs facetedges such as the origins corresponding to
        vertices extremes of each tetrahedron of Triangulation.
     *)
    VAR poly : REF ARRAY OF Quadp := NEW(REF ARRAY OF Quadp, top.NP);
    BEGIN
      FOR i := 0 TO top.NP-1 DO
        WITH
          da = top.region[i],
          a0 = Tors(da), 
          db = Clock(Enext_1(da)),
          b0 = Tors(db) DO
          (*<* ASSERT DegreeOfVertex(a0) = 4 *>  *)       (* consuming time *)
          WITH
            a1 = Enext(a0),
            a2 = Enext(a1),
            a3 = Enext_1(b0)
          DO
	    IF Enext(a2) # a0 THEN
	      Wr.PutText(stderr, "\nTriangulation: This topology isn't" &
	                 " a triangulation\n");
              <* ASSERT Enext(a2) = a0 *>
	    END;  
            <* ASSERT Pneg(a0).num = i *>
            IF Pneg(b0) # NIL THEN
              <* ASSERT Pneg(a0) = Pneg(b0) *>
            END;
            <* ASSERT a0 # b0 *>
            poly[i] := Quadp{a0,a1,a2,a3};     
	  END
        END
      END;
      RETURN poly;
    END VertexPoly;

  VAR poly1 : REF ARRAY OF Quadv := NEW(REF ARRAY OF Quadv, top.NP);
      n : CARDINAL := 0;
  BEGIN
    WITH poly = VertexPoly(top), v = OrgV(a) DO
      FOR i := 0 TO LAST(poly^) DO
         IF (OrgV(poly[i][0]) = v) OR (OrgV(poly[i][1]) = v) OR
            (OrgV(poly[i][2]) = v) OR (OrgV(poly[i][3]) = v)
         THEN
            FOR j := 0 TO 3 DO
               IF OrgV(poly[i][j]) = v
                 THEN
                   WITH                   
                     uv = OrgV(poly[i][j]),
                     dv = OrgV(Enext(poly[i][j])),
                     pv = OrgV(Enext_1(poly[i][j])),
                     rv = OrgV(Enext_1(Fnext_1(poly[i][j])))
                   DO
                     poly1[n] := Quadv{uv,dv,pv,rv};
                     INC(n);
                   END
               END
             END
         END  
      END;
      WITH
        r = NEW(REF ARRAY OF Quadv, n)
      DO
        r^ := SUBARRAY(poly1^, 0, n);
        RETURN r;
      END;
    END;
  END StarOfVertex;    

PROCEDURE NumberPolyOfStar(quad: REF ARRAY OF Quadv): CARDINAL =
  BEGIN
    WITH n = NUMBER(quad^) DO RETURN n; END;
  END NumberPolyOfStar;

PROCEDURE ComputeAllVertexNormals(
     READONLY top: Topology; 
     READONLY c: Coords;
     ): REF ARRAY OF LR4.T =
  BEGIN
    WITH
      rvn = NEW(REF ARRAY OF LR4.T, top.NV),
      vn = rvn^
    DO
      FOR i := 0 TO top.NV-1 DO 
        vn[i] := VertexNormal(top.out[i], c, top)
      END;
      RETURN rvn
    END;
  END ComputeAllVertexNormals; 

PROCEDURE ComputeAllEdgeNormals(
    READONLY top: Topology; 
    READONLY c: Coords;
    ): REF ARRAY OF LR4.T =
  BEGIN
    WITH
      rvn = NEW(REF ARRAY OF LR4.T, top.NE),
      vn = rvn^
    DO
      FOR i := 0 TO top.NE-1 DO 
        vn[i] := EdgeNormal(top.edge[i].pa, c);
      END;
      RETURN rvn
    END;
  END ComputeAllEdgeNormals;

PROCEDURE ComputeAllFaceNormals(
    READONLY top: Topology; 
    READONLY c: Coords;
    ): REF ARRAY OF LR4.T =
  BEGIN
    WITH
      rvn = NEW(REF ARRAY OF LR4.T, top.NF),
      vn = rvn^
    DO
      FOR i := 0 TO top.NF-1 DO 
        vn[i] := FaceNormal(top.face[i].pa, c);
      END;
      RETURN rvn
    END;
  END ComputeAllFaceNormals;

PROCEDURE ComputeAllPolyhedronNormals(
    READONLY top: Topology; 
    READONLY c: Coords;
    ): REF ARRAY OF LR4.T =
  BEGIN
    WITH
      rvn = NEW(REF ARRAY OF LR4.T, top.NP),
      vn = rvn^
    DO
      FOR i := 0 TO top.NP-1 DO 
        vn[i] := PolyNormal(Tors(top.region[i]), c);
      END;
      RETURN rvn
    END;
  END ComputeAllPolyhedronNormals; 

PROCEDURE Barycenter(READONLY top: Topology; READONLY c: Coords; all: BOOL): LR4.T =
  VAR b: LR4.T := LR4.T{0.0d0, ..};
      N: CARDINAL := 0;
  BEGIN
    FOR i := 0 TO LAST(c) DO 
      WITH v = top.vertex[i] DO
        IF all OR v.exists THEN b := LR4.Add(b, c[i]); INC(N) END;
      END;
    END;
    IF N = 0 THEN
      RETURN b
    ELSE
      RETURN LR4.Scale(1.0d0/FLOAT(N, LONGREAL), b)
    END
  END Barycenter;

PROCEDURE PartialBarycenter(READONLY v: Vertices; READONLY c: Coords; all: BOOL): LR4.T =
  VAR b: LR4.T := LR4.T{0.0d0, ..};
      N: CARDINAL := 0;
  BEGIN
    FOR k := 0 TO LAST(c) DO 
      WITH vk = v[k], i = vk.num DO
        IF all OR vk.exists THEN b := LR4.Add(b, c[i]); INC(N) END;
      END;
    END;
    IF N = 0 THEN
      RETURN b
    ELSE
      RETURN LR4.Scale(1.0d0/FLOAT(N, LONGREAL), b)
    END
  END PartialBarycenter;

PROCEDURE MeanVertexDistance(
    READONLY top: Topology;
    READONLY c: Coords;
    READONLY ctr: LR4.T;
    all: BOOL;
  ): LONGREAL =
  VAR S: LONGREAL := 0.0d0;
      N: CARDINAL := 0;
  BEGIN
    FOR i := 0 TO LAST(c) DO 
      WITH v = top.vertex[i] DO
        IF all OR v.exists THEN 
          WITH d2 = LR4.DistSqr(ctr, c[i]) DO
            S := S + d2
          END;
          INC(N)
        END
      END
    END;
    IF N = 0 THEN
      RETURN 1.0d0
    ELSE
      RETURN Math.sqrt(S/FLOAT(N,LONGREAL))
    END
  END MeanVertexDistance;

PROCEDURE MaxVertexDistance(
    READONLY top: Topology; 
    READONLY c: Coords; 
    READONLY ctr: LR4.T;
    all: BOOL;
  ): LONGREAL =
  VAR radius := 0.0d0;
  BEGIN
    FOR i := 0 TO top.NV-1 DO
      WITH v = top.vertex[i] DO
        IF all OR v.exists THEN 
          WITH
            d = LR4.Dist(ctr, c[i])
          DO
            IF d > radius THEN radius := d END
          END
        END
      END
    END;   	  
    RETURN radius
  END MaxVertexDistance;

PROCEDURE MeanThickness(
    READONLY top: Topology; 
    READONLY c: Coords; 
    READONLY ctr: LR4.T;
    READONLY norm: LR4.T;
    all: BOOL;
  ): LONGREAL =
  VAR S : LONGREAL := 0.0d0;
      N: CARDINAL := 0;
  BEGIN
    FOR i := 0 TO top.NV-1 DO
      WITH v = top.vertex[i] DO
        IF all OR v.exists THEN 
          WITH d = LR4.Dot(LR4.Sub(c[i],ctr), norm) DO
            S := S + d*d;
          END;
          INC(N)
        END
      END
    END;
    IF N = 0 THEN
      RETURN 1.0d0
    ELSE
      RETURN Math.sqrt(S/FLOAT(N,LONGREAL))
    END
  END MeanThickness;

PROCEDURE MeanEdgeLength(READONLY top: Topology; READONLY c: Coords; all: BOOL): LONGREAL =
  VAR S: LONGREAL := 0.0d0;
      N: CARDINAL := 0;
  BEGIN
    FOR i := 0 TO top.NE-1 DO 
      WITH e = top.edge[i] DO
        IF all OR e.exists THEN 
          WITH o = OrgV(e.pa).num, d = OrgV(Clock(e.pa)).num DO
            S := S + LR4.DistSqr(c[o], c[d])
          END;
          INC(N);
        END
      END
    END;
    IF N = 0 THEN
      RETURN 1.0d0
    ELSE
      RETURN Math.sqrt(S/FLOAT(N,LONGREAL))
    END
  END MeanEdgeLength;

PROCEDURE MeanPolyhedronNormal(READONLY top: Topology; READONLY c: Coords; all: BOOL): LR4.T =
  VAR norm: LR4.T := LR4.T{0.0d0, ..};
  BEGIN
    FOR i := 0 TO top.NP-1 DO
      WITH f = Tors(top.region[i]), p = PnegP(f) DO
        IF all OR (p # NIL AND p.exists) THEN
          WITH n = PolyCross(f, c) DO
            norm := LR4.Add(norm, n)
          END
        END
      END
    END;
    WITH m = LR4.Norm(norm) DO
      IF m < 1.0d-20 THEN
        RETURN LR4.T{0.0d0, 0.0d0, 0.0d0, 1.0d0}
      ELSE
        RETURN LR4.Scale(1.0d0/m, norm);
      END
    END
  END MeanPolyhedronNormal;
 
PROCEDURE NormalizeVertexDistances(READONLY top: Topology; VAR c: Coords; all: BOOL) =
  BEGIN
    WITH b = Barycenter(top, c, all) DO
      Displace(top, LR4.Neg(b), c)
    END;
    WITH s = MeanVertexDistance(top, c, LR4.T{0.0d0, ..}, all) DO
      Scale(top, 1.0d0/s, c)
    END;
  END NormalizeVertexDistances;

PROCEDURE NormalizeEdgeLengths(READONLY top: Topology; VAR c: Coords; all: BOOL) =
  BEGIN
    WITH b = Barycenter(top, c, all) DO
      Displace(top, LR4.Neg(b), c) 
    END;
    WITH s = MeanEdgeLength(top, c, all) DO
      Scale(top, 1.0d0/s, c)
    END;
  END NormalizeEdgeLengths;

PROCEDURE MakePolyhedronTopology(a: Pair): PolyhedronTopology =
  VAR
    ne, nv: CARDINAL;
    b,  bn: Pair;
    ptop: PolyhedronTopology;
  BEGIN
    WITH
      pneg = Org(a),
      star = Octf.NumberEdgesForDegree(ARRAY OF Pair{a})^
    DO
      (* Gather faces: *)
      ptop.NF := NUMBER(star);
      ptop.fRef := NEW(REF ARRAY OF Pair, ptop.NF);
      ne := 0;
      FOR i := 0 TO ptop.NF-1 DO
        WITH
          side = Octf.Sdual(star[i])
        DO
          ptop.fRef[i] := side;
          <* ASSERT Pneg(ptop.fRef[i]) = pneg *>
          b := side;
          REPEAT
            INC(ne);
            WITH e = b.facetedge.edge DO
              e.xmark := FALSE;
              OrgV(b).xmark := FALSE; (* OrgV(b).xmark := FALSE *)
            END;
            b := Enext(b)
          UNTIL b = side;  
        END
      END;

      (* Gather edges: *)
      <* ASSERT ne MOD 2 = 0 *>
      ptop.NE := ne DIV 2;
      ptop.eRef := NEW(REF ARRAY OF Pair, ptop.NE);
      ne := 0;
      FOR i := 0 TO ptop.NF-1 DO
        WITH side = Octf.Sdual(star[i]) DO
          b  := side;
          REPEAT
            WITH e = b.facetedge.edge DO
              IF NOT e.xmark THEN
                ptop.eRef[ne] := b;
                <* ASSERT Pneg(ptop.eRef[ne]) = pneg *>
                bn := b;
                REPEAT
                  bn.facetedge.edge.xmark := TRUE;
                  bn := Fnext(bn);
                UNTIL (bn = b);
                INC(ne);
              END; 
            END;
            b := Enext(b)
          UNTIL b = side
        END
      END;

      (* Gather vertices: *)
      ptop.NV := 2 + ptop.NE - ptop.NF;
      ptop.vRef := NEW(REF ARRAY OF Pair, ptop.NV);
      nv := 0;
      FOR i := 0 TO ptop.NE-1 DO 
        WITH
          eu = ptop.eRef[i],
          u = OrgV(eu), 
          ev = Clock(Fnext_1(eu)),
          v = OrgV(ev)
        DO 
          IF NOT u.xmark THEN
            ptop.vRef[nv] := eu;
            <* ASSERT Pneg(ptop.vRef[nv]) = pneg *>
            INC(nv);
            u.xmark := TRUE;
          END;
          IF NOT v.xmark THEN ptop.vRef[nv] := ev;
            <* ASSERT Pneg(ptop.vRef[nv]) = pneg *>
            INC(nv);
            v.xmark := TRUE END;
          eu.facetedge.edge.xmark := FALSE;
        END;
      END;
      FOR i := 0 TO ptop.NV-1 DO WITH u = OrgV(ptop.vRef[i]) DO u.xmark := FALSE END END;
    END;
    RETURN ptop
  END MakePolyhedronTopology;

(* ========== INPUT/OUTPUT ================== *)

CONST
  Boole = Mis.Boole;
  AlphaChars = Mis.AlphaChars;

PROCEDURE ReadTopology(rd: Rd.T): TopCom =
  VAR
    top: Topology; 
    comments: TEXT;
    n : CHAR;

  BEGIN
    (* Topology *)
    ReadHeader(rd,"topology","99-08-25");
    comments := Mis.ReadCommentsJS(rd, '|');

    (* Element counts: *)
    Lex.Skip(rd, cs := AlphaChars);
    top.NV := Lex.Int(rd);
    Lex.Skip(rd, cs := AlphaChars);
    top.NE := Lex.Int(rd);
    Lex.Skip(rd, cs := AlphaChars);
    top.NF := Lex.Int(rd);
    Lex.Skip(rd, cs := AlphaChars);
    top.NP := Lex.Int(rd);
    Lex.Skip(rd, cs := AlphaChars);
    top.NFE := Lex.Int(rd);
    Lex.Skip(rd, cs := AlphaChars);
    EVAL Lex.Int(rd); (* top.der *)
    Lex.Skip(rd,  cs := AlphaChars);
    EVAL Lex.Int(rd); (* top.bdr *)
    Lex.Skip(rd);

    WITH
      map = NEW(REF ARRAY OF Octf.FacetEdge, top.NFE)^
    DO

      (* Create vertex records: *)
      top.vertex := NEW(REF Vertices, top.NV);
      top.out := NEW(REF ARRAY OF Pair, top.NV);
      FOR i := 0 TO top.NV-1 DO top.vertex[i] := MakeVertex(); END;
      
      (* Create edges records *)
      top.edge := NEW(REF ARRAY OF Octf.Edge, top.NE);
      FOR i := 0 TO top.NE-1 DO
        top.edge[i] := Octf.MakeEdge();
        (* 
          FOR j := 0 TO 1 DO
            top.edge[i].vertex[j] := MakeVertex();
          END;
        *)
      END;

      (* Create face records: *)
      top.face := NEW(REF ARRAY OF Face, top.NF);
      FOR i := 0 TO top.NF-1 DO
        top.face[i] := Octf.MakeFace();
        (*
          top.face[i].vertex := NEW(REF ARRAY OF Node,top.der);
          FOR j := 0 TO top.der-1 DO
            top.face[i].vertex^[j] := MakeVertex();
          END;
        *)
      END;
      
      (* Create polyhedra records: *)
      top.polyhedron := NEW(REF ARRAY OF Polyhedron, top.NP);
      top.region := NEW(REF ARRAY OF Pair, top.NP);
      
      FOR i := 0 TO top.NP-1 DO
        top.polyhedron[i] := MakePolyhedron();
        (*
          IF  top.der = 3 THEN
            IF top.bdr = 2 THEN
              (* Obs: bdr=2 indicate that the cells are octahedra. *)
              top.polyhedron[i].vertex := NEW(REF ARRAY OF Node,top.der+3);
              FOR j := 0 TO top.der+2 DO
                top.polyhedron[i].vertex^[j] := MakeVertex();
              END;
            ELSE
                top.polyhedron[i].vertex := NEW(REF ARRAY OF Node,top.der+1);
              FOR j := 0 TO top.der DO
                top.polyhedron[i].vertex^[j] := MakeVertex();
              END
            END
          ELSIF top.der = 4 THEN
            top.polyhedron[i].vertex := NEW(REF ARRAY OF Node,2*top.der);
            FOR j := 0 TO 2*top.der-1 DO
              top.polyhedron[i].vertex^[j] := MakeVertex();
            END
          ELSIF top.der = 2 THEN
            top.polyhedron[i].vertex := NEW(REF ARRAY OF Node,2);
            FOR j := 0 TO 1 DO
              top.polyhedron[i].vertex^[j] := MakeVertex();
            END
          ELSIF top.der = 5 THEN
            top.polyhedron[i].vertex := NEW(REF ARRAY OF Node,4*top.der);
            FOR j := 0 TO 19 DO
              top.polyhedron[i].vertex^[j] := MakeVertex();
            END
          END
        *)
      END;
      
      (* Create facetedge records: *)
      top.facetedge := NEW(REF ARRAY OF Octf.Pair, top.NFE);
      FOR i := 0 TO top.NFE-1 DO
        top.facetedge[i] := MakeFacetEdge();
        top.facetedge[i].facetedge.num := i;
        map[i] := top.facetedge[i].facetedge;
      END;

      (* Read edge data: *)
      EVAL Mis.ReadCommentsJS(rd, '|');
      FOR j := 0 TO top.NE-1 DO 
        Lex.Skip(rd);
        WITH 
          ne = Lex.Int(rd),      (* index to edge *) 
          e = top.edge[ne]
        DO
          <* ASSERT ne = j *>
          e.num := ne;
          e.pa := Octf.ReadPair(rd,map);
        END;
      END;
      Lex.Skip(rd);

      (* Read face data: *)
      EVAL Mis.ReadCommentsJS(rd, '|');
      FOR j := 0 TO top.NF-1 DO 
        Lex.Skip(rd);
        WITH 
          nf = Lex.Int(rd),      (* index to face *) 
          f = top.face[nf]
        DO
          <* ASSERT nf = j *>
          f.num := nf;
          f.pa := Octf.ReadPair(rd,map);
        END;
      END;
      Lex.Skip(rd);

      (* Read facetedge data: *)
      EVAL Mis.ReadCommentsJS(rd, '|');
      FOR j := 0 TO top.NFE-1 DO 
        Lex.Skip(rd);
        WITH 
          nfe = Lex.Int(rd),        (* index to facetedge *)
          fe = NARROW(top.facetedge[nfe].facetedge, FacetEdge)
        DO
          <* ASSERT nfe = j *>
          <* ASSERT top.facetedge[nfe].bits = 0 *>
          fe.num := nfe;
          Octf.ReadFacetEdge(rd, fe, map);
          FOR k := 0 TO 3  DO 
            Lex.Skip(rd);
            n := Rd.GetChar(rd);
            IF n # '-' THEN
              Rd.UnGetChar(rd);
              WITH
                m = Lex.Int(rd),
                vf = fe.org[k]
              DO
                IF Rd.GetChar(rd) = 'v' THEN 
                  vf := top.vertex[m];
                  top.out[m] := Pair{facetedge := fe, bits := 2*k};
                ELSE 
                  vf := top.polyhedron[m];
                  top.region[m] := Pair{facetedge := fe, bits := 2*k};
                END;
              END;
            END;
          END;
          Lex.Skip(rd);
          fe.face.num := Lex.Int(rd);
          <* ASSERT Rd.GetChar(rd) = 'f' *>
          Lex.Skip(rd);
          fe.edge.num := Lex.Int(rd);
          <* ASSERT Rd.GetChar(rd) = 'e' *>
        END;
      END;
      ReadFooter(rd,"topology");
      RETURN TopCom{top, comments}
    END
  END ReadTopology;

PROCEDURE ReadMaterials(
  rd: Rd.T;
  READONLY top: Topology;
  ro_te: BOOL := FALSE;
  ) =
  BEGIN
    (* Materials *)
    ReadHeader(rd,"materials","99-08-25");

    (* Read vertex data Materials: *)
    EVAL Mis.ReadCommentsJS(rd, '|');
    FOR j := 0 TO top.NV-1 DO 
      Lex.Skip(rd);
      WITH 
        nv = Lex.Int(rd),    (* index to vertex *)
        v = top.vertex[nv]   
      DO 
        <* ASSERT nv = j *>
        v.num := nv;
        Lex.Skip(rd); 
        v.exists := Mis.ReadBool(rd);
        Lex.Skip(rd);
        v.fixed := Mis.ReadBool(rd);
        Lex.Skip(rd);
        WITH cc = v.color DO
          cc[0] := Lex.Real(rd);
          Lex.Skip(rd);
          cc[1] := Lex.Real(rd);
          Lex.Skip(rd);
          cc[2] := Lex.Real(rd); 
        END;
        Lex.Skip(rd);
        WITH tt = v.transp DO
          tt[0] := Lex.Real(rd);
          Lex.Skip(rd);
          tt[1] := Lex.Real(rd);
          Lex.Skip(rd);
          tt[2] := Lex.Real(rd);
        END;
        Lex.Skip(rd);
        v.radius := Lex.Real(rd);
        Lex.Skip(rd);
        v.label := Rd.GetText(rd,2); 
      END      
    END;

    Lex.Skip(rd);

    (* Read edge data materials: *)
    EVAL Mis.ReadCommentsJS(rd, '|');
    FOR j := 0 TO top.NE-1 DO 
      Lex.Skip(rd);
      WITH 
        ne = Lex.Int(rd),      (* index to edge *) 
        e = top.edge[ne]
      DO
        <* ASSERT ne = j *>
        e.num := ne;
        Lex.Skip(rd);
        e.exists := Mis.ReadBool(rd);
        Lex.Skip(rd);
        WITH cc = e.color DO
          cc[0] := Lex.Real(rd);
          Lex.Skip(rd);
          cc[1] := Lex.Real(rd);
          Lex.Skip(rd);
          cc[2] := Lex.Real(rd); 
        END;
        Lex.Skip(rd);
        WITH tt = e.transp DO
          tt[0] := Lex.Real(rd);
          Lex.Skip(rd);
          tt[1] := Lex.Real(rd);
          Lex.Skip(rd);
          tt[2] := Lex.Real(rd);
        END;
        Lex.Skip(rd);
        e.radius := Lex.Real(rd);
        Lex.Skip(rd);
        EVAL Mis.ReadBool(rd); (* e.degenerate *)
        Lex.Skip(rd);
        WITH
          n = Rd.GetChar(rd)
        DO
          IF n # '-' THEN
	        Rd.UnGetChar(rd);
            e.root := Lex.Int(rd);
          ELSE
            e.root := -1;
          END
        END
      END
    END;
    Lex.Skip(rd);

    (* Read face data materials: *)
    EVAL Mis.ReadCommentsJS(rd, '|');
    FOR j := 0 TO top.NF-1 DO 
      Lex.Skip(rd);
      WITH 
        nf = Lex.Int(rd),      (* index to face *) 
        f = top.face[nf]
      DO
        <* ASSERT nf = j *>
        f.num := nf;
        Lex.Skip(rd);
        f.exists := Mis.ReadBool(rd);
        Lex.Skip(rd);
        WITH cc = f.color DO
          cc[0] := Lex.Real(rd);
          Lex.Skip(rd);
          cc[1] := Lex.Real(rd);
          Lex.Skip(rd);
          cc[2] := Lex.Real(rd); 
        END;
        Lex.Skip(rd);
        WITH tt = f.transp DO
          tt[0] := Lex.Real(rd);
          Lex.Skip(rd);
          tt[1] := Lex.Real(rd);
          Lex.Skip(rd);
          tt[2] := Lex.Real(rd);
        END;
        Lex.Skip(rd);
        EVAL Mis.ReadBool(rd); (* f.degenerate *)
        Lex.Skip(rd);
        WITH
          n = Rd.GetChar(rd)
        DO
          IF n # '-' THEN
	        Rd.UnGetChar(rd);
            f.root := Lex.Int(rd);
          ELSE
            f.root := -1;
          END
        END
      END
    END;
    Lex.Skip(rd);

    (* Read polyhedron data materials: *)
    IF top.NP # 0 THEN
      EVAL Mis.ReadCommentsJS(rd, '|');
    END;
    FOR j := 0 TO top.NP-1 DO 
      Lex.Skip(rd);
      WITH 
        np = Lex.Int(rd),      (* index to polyhedron *)
        p = top.polyhedron[np]
      DO
        <* ASSERT np = j *>
        p.num := np;
        Lex.Skip(rd);
        p.exists := Mis.ReadBool(rd);
        Lex.Skip(rd);
        WITH cc = p.color DO
          cc[0] := Lex.Real(rd);
          Lex.Skip(rd);
          cc[1] := Lex.Real(rd);
          Lex.Skip(rd);
          cc[2] := Lex.Real(rd); 
        END;
        Lex.Skip(rd);
        WITH tt = p.transp DO
          tt[0] := Lex.Real(rd);
          Lex.Skip(rd);
          tt[1] := Lex.Real(rd);
          Lex.Skip(rd);
          tt[2] := Lex.Real(rd);
        END;
        Lex.Skip(rd);
        EVAL Mis.ReadBool(rd); (* p.degenerate *)
	IF ro_te THEN
	  Lex.Skip(rd);
	  WITH
	    n = Rd.GetChar(rd)
	  DO
	    IF n # '-' THEN
	      Rd.UnGetChar(rd);
	      p.root := Lex.Int(rd);
	    ELSE
	      p.root := -1;
	    END
	  END
	END
      END
    END;

    ReadFooter(rd,"materials");
    Rd.Close(rd);
    CheckOutAndRegion(top);
  END ReadMaterials;

PROCEDURE ReadToMa(name: TEXT; ro_te: BOOL := FALSE): TopCom =
  (* Where ro_te meaning "root tetrahedron". *)
  VAR
    tc: TopCom;
  BEGIN
    WITH 
      ntp = name & ".tp",
      rtp = FileRd.Open(ntp)
    DO
      Wr.PutText(stderr, "reading " & ntp & "\n");
      tc := ReadTopology(rtp);
      Rd.Close(rtp);
    END;  
    WITH 
      nma = name & ".ma",
      rma = FileRd.Open(nma)
    DO
      Wr.PutText(stderr, "reading " & nma & "\n");
      ReadMaterials(rma, tc.top, ro_te);
      Rd.Close(rma);
    END;  
    RETURN tc
  END ReadToMa;

PROCEDURE ReadState(name: TEXT): REF Coords =
  <* FATAL Rd.Failure, Thread.Alerted,FloatMode.Trap, Lex.Error, OSError.E *>
  VAR
    c: REF Coords;
    comments: TEXT;
    nv : CARDINAL;
  BEGIN
    WITH 
      nst = name & ".st",
      rst = FileRd.Open(nst)
    DO
      Wr.PutText(stderr, "reading " & nst & "\n");
      
      (* Read Headers File Formats*)
      ReadHeader(rst,"state","99-08-25");

      (* Element counts: *)
      Lex.Skip(rst, cs := AlphaChars);
      nv := Lex.Int(rst);
      Lex.Skip(rst);
      comments := Mis.ReadCommentsJS(rst, '|');
      
      c := NEW(REF Coords, nv);
      (* Read vertex data state: *)
      FOR j := 0 TO nv-1 DO 
        Lex.Skip(rst);
        WITH 
          nv = Lex.Int(rst) 
        DO 
          WITH cv = c[nv] DO
            cv[0] := Lex.LongReal(rst);
            Lex.Skip(rst);
            cv[1] := Lex.LongReal(rst);
            Lex.Skip(rst);
            cv[2] := Lex.LongReal(rst);
            Lex.Skip(rst);
            cv[3] := Lex.LongReal(rst);
          END
        END 
      END;

      ReadFooter(rst,"state");
      Rd.Close(rst);        
      RETURN c;
    END;  
  END ReadState;

PROCEDURE WriteState(
    name: TEXT;
    READONLY top: Topology;
    READONLY c: Coords;
    comments: TEXT := "";
    ) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  <* FATAL OSError.E *>
  BEGIN
    WITH 
      st = FileWr.Open(name & ".st"),
      vWidth = Mis.NumDigits(top.NV - 1)
    DO
      PROCEDURE WriteCoord(x: LONGREAL) =
        BEGIN
          Wr.PutText(st, Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Sci, prec := 3),
	  7));
        END WriteCoord;

      PROCEDURE WritePoint(READONLY c: LR4.T) =
        BEGIN
          WriteCoord(c[0]);
          Wr.PutText(st, " ");
          WriteCoord(c[1]);
          Wr.PutText(st, " ");
          WriteCoord(c[2]);
          Wr.PutText(st, " ");
          WriteCoord(c[3]);
        END WritePoint;
      
      BEGIN
        WriteHeader(st,"state","99-08-25");
        Wr.PutText(st, "vertices ");
        Wr.PutText(st, Fmt.Int(top.NV) & "\n");
        IF NOT Text.Empty(comments) THEN
          Mis.WriteCommentsJS(st, comments & "\n", '|')
        END;

        Mis.WriteCommentsJS(st, "\nVertex data:\n", '|');
        
        FOR i := 0 TO top.NV-1 DO
          WITH v = top.vertex[i] DO
            (* state *)
            Wr.PutText(st, Fmt.Pad(Fmt.Int(v.num), vWidth)); 
            Wr.PutText(st, "  ");
            WritePoint(c[v.num]);
            Wr.PutText(st, "\n");
          END       
        END;
      END;
      WriteFooter(st, "state");
      Wr.Close(st);
    END;     
  END WriteState;

PROCEDURE WriteMaterials(
   name: TEXT;
   READONLY top: Topology;
   comments: TEXT := "";
   ro_te: BOOL := FALSE;  (* root tetrahedron *)
  ) =
  <* FATAL Wr.Failure, Thread.Alerted, OSError.E *>
  VAR pWidth : INTEGER;
  BEGIN
    IF top.NP = 0 THEN pWidth := 2 ELSE pWidth := Mis.NumDigits(top.NP-1) END;
    WITH 
      ma = FileWr.Open(name & ".ma"), 
      vWidth = Mis.NumDigits(top.NV-1),
      fWidth = Mis.NumDigits(top.NF-1), 
      eWidth = Mis.NumDigits(top.NE-1)
    DO 
      PROCEDURE WriteIntensity(r: REAL) =
        BEGIN
          Wr.PutText(ma, Fmt.Real(r, Fmt.Style.Fix, prec := 2));
        END WriteIntensity;

      PROCEDURE WriteColor(READONLY c: R3.T) =
        BEGIN
          WriteIntensity(c[0]);
          Wr.PutText(ma, " ");                     
          WriteIntensity(c[1]);
          Wr.PutText(ma, " ");                     
          WriteIntensity(c[2]);
        END WriteColor;
        
      PROCEDURE WriteRadius(r: REAL) =
        BEGIN
          IF r = 0.02
           THEN
             Wr.PutText(ma, "0.020");
            ELSE
             Wr.PutText(ma,Fmt.Real(r, prec := 4)); 
          END
        END WriteRadius;

      PROCEDURE WriteLabel(label: TEXT) =
        BEGIN
          Wr.PutText(ma, label);
        END WriteLabel;

      BEGIN
        WriteHeader(ma,"materials","99-08-25");
        IF NOT Text.Empty(comments) THEN
          WriteCommentsJS(ma, comments & "\n", '|')
        END;
        WITH
          m = Mis.NumDigits(top.NP)
        DO
          WriteCommentsJS(ma,"vertices    " & Fmt.Pad(Fmt.Int(top.NV),m), '|');
          WriteCommentsJS(ma,"edges       " & Fmt.Pad(Fmt.Int(top.NE),m), '|');
          WriteCommentsJS(ma,"faces       " & Fmt.Pad(Fmt.Int(top.NF),m), '|');
          WriteCommentsJS(ma,"polyhedra   " & Fmt.Int(top.NP), '|');
        END;
        WriteCommentsJS(ma, "\nVertex data:\n", '|');

        FOR i := 0 TO top.NV-1 DO
          WITH v = top.vertex[i] DO
            (* materials *)
            Wr.PutText(ma, Fmt.Pad(Fmt.Int(v.num), vWidth)); 
            Wr.PutText(ma, "  ");
            Wr.PutText(ma, Fmt.Char(Boole[v.exists]));
            Wr.PutText(ma, " ");
            Wr.PutText(ma, Fmt.Char(Boole[v.fixed]));
            Wr.PutText(ma, "  ");
            WriteColor(v.color);
            Wr.PutText(ma, "  ");
            WriteColor(v.transp);
            Wr.PutText(ma, "  ");
            WriteRadius(v.radius);
            Wr.PutText(ma, "  ");
            WriteLabel(v.label);
            Wr.PutText(ma, "\n");
          END       
        END;

        WriteCommentsJS(ma, "\nEdge data:\n", '|');
        FOR i := 0 TO top.NE-1 DO
          WITH e = top.edge[i]  DO
            (* materials *)
            Wr.PutText(ma, Fmt.Pad(Fmt.Int(e.num), eWidth)); 
            Wr.PutText(ma, "  ");
            Wr.PutText(ma, Fmt.Char(Boole[e.exists]));
            Wr.PutText(ma, "  ");
            WriteColor(e.color);
            Wr.PutText(ma, "  ");
            WriteColor(e.transp);
            Wr.PutText(ma, "  ");
            WriteRadius(e.radius);
            Wr.PutText(ma, "  ");
            WITH degenerate = FALSE DO Wr.PutText(ma, Fmt.Char(Boole[degenerate])) END;
            Wr.PutText(ma, "  ");
	    IF e.root = -1 THEN
	      Wr.PutText(ma, " - ");
	    ELSE
              Wr.PutText(ma, Fmt.Pad(Fmt.Int(e.root), eWidth));
            END;
            Wr.PutText(ma, "\n");
          END
        END;

        WriteCommentsJS(ma, "\nFace data:\n", '|');

        FOR i := 0 TO top.NF-1 DO
          WITH f = top.face[i] DO 
            Wr.PutText(ma, Fmt.Pad(Fmt.Int(f.num), fWidth)); 
            Wr.PutText(ma, "  ");
            Wr.PutText(ma, Fmt.Char(Boole[f.exists]));
            Wr.PutText(ma, "  ");
            WriteColor(f.color);
            Wr.PutText(ma, "  ");
            WriteColor(f.transp);
            Wr.PutText(ma, "  ");
            WITH degenerate = FALSE DO Wr.PutText(ma, Fmt.Char(Boole[degenerate])) END;
            Wr.PutText(ma, "  ");
	    IF f.root = -1 THEN
	      Wr.PutText(ma, " - ");
	    ELSE
              Wr.PutText(ma, Fmt.Pad(Fmt.Int(f.root), fWidth));
            END;
            Wr.PutText(ma, "\n");
          END
        END;

        IF top.NP # 0 THEN
          WriteCommentsJS(ma, "\nPolyhedron data:\n", '|');
        END;

        FOR i := 0 TO top.NP-1 DO
          WITH
	    p = top.polyhedron[i]
	  DO
            Wr.PutText(ma, Fmt.Pad(Fmt.Int(p.num), pWidth)); 
            Wr.PutText(ma, "  ");
            Wr.PutText(ma, Fmt.Char(Boole[p.exists]));
            Wr.PutText(ma, "  ");
            WriteColor(p.color);
            Wr.PutText(ma, "  ");
            WriteColor(p.transp);
            Wr.PutText(ma, "  ");
            WITH degenerate = FALSE DO Wr.PutText(ma, Fmt.Char(Boole[degenerate])) END;
	    IF ro_te THEN
	      Wr.PutText(ma, "  ");
              IF p.root = -1 THEN
	        Wr.PutText(ma, " - ");
	      ELSE
                Wr.PutText(ma, Fmt.Pad(Fmt.Int(p.root), pWidth));
              END
	    END;
            Wr.PutText(ma, "\n");
          END
        END
      END;
      WriteFooter(ma, "materials");
      Wr.Close(ma);
    END;     
  END WriteMaterials;

PROCEDURE CollectFaceEdges(a: Pair; VAR re: REF ARRAY OF Pair; VAR ne: CARDINAL) =
  VAR b: Pair;
  BEGIN
    b := a; ne := 0;
    REPEAT
      <* ASSERT b.facetedge.face = a.facetedge.face *>
      IF re = NIL OR ne >= NUMBER(re^) THEN
        WITH se = NEW(REF ARRAY OF Pair, MAX(10, 2*ne)) DO
          IF re # NIL THEN SUBARRAY(se^, 0, ne) := re^ END;
          re := se
        END
      END;
      re[ne] := b; INC(ne);
      b := Enext(b)
    UNTIL b = a;
  END CollectFaceEdges;
  
PROCEDURE CollectPolyhedronFaces(a: Pair; VAR rf: REF ARRAY OF Pair; VAR nf: CARDINAL) =

  VAR scanned: CARDINAL := 0;
      
  PROCEDURE StackFace(c: Pair) =
    BEGIN
      <* ASSERT Pneg(c) = Pneg(a) *>
      c.facetedge.face.xmark := TRUE;
      IF rf = NIL OR nf >= NUMBER(rf^) THEN
        WITH sf = NEW(REF ARRAY OF Pair, MAX(10, 2*nf)) DO
          IF rf # NIL THEN SUBARRAY(sf^, 0, nf) := rf^ END;
          rf := sf
        END
      END;
      rf[nf] := c; INC(nf);
    END StackFace;

  PROCEDURE VisitFace(b: Pair) =
    VAR c: Pair := Clock(Fnext_1(b));
    BEGIN
      REPEAT
        <* ASSERT Pneg(c) = Pneg(a) *>
        WITH g = c.facetedge.face DO
          IF NOT g.xmark THEN 
            StackFace(c);
            (* We have only one mark per face, but a face may incide twice on "p": *)
            (* Assumes that the cell's boundary is oriented! *)
            IF PposP(c) = Pneg(a) THEN StackFace(Clock(c)) END;
          END
        END;
        c := Enext(c)
      UNTIL c = b;
    END VisitFace;
    
  BEGIN
    nf := 0; scanned := 0;
    StackFace(a);
    WHILE scanned < nf DO
      VisitFace(rf[scanned]);
      INC(scanned)
    END;
    (* Clear "xmark", just in case: *)
    FOR k := 0 TO nf - 1 DO
      WITH b = rf[k] DO b.facetedge.face.xmark := FALSE END
    END;
  END CollectPolyhedronFaces;

PROCEDURE CollectPolyhedronEdges(
    a: Pair; 
    VAR re: REF ARRAY OF Pair; 
    VAR ne: CARDINAL;
    VAR rf: REF ARRAY OF Pair;  (* Working area *)
  ) =
  VAR nf: CARDINAL;
      b: Pair;
  BEGIN
    CollectPolyhedronFaces(a, rf, nf);
    ne := 0;
    FOR k := 0 TO nf-1 DO
      WITH f = rf[k] DO
        <* ASSERT Pneg(f) = Pneg(a) *>
        b := f;
        REPEAT
          WITH e = b.facetedge.edge DO 
            IF NOT e.xmark THEN 
              IF re = NIL OR ne >= NUMBER(re^) THEN
                WITH se = NEW(REF ARRAY OF Pair, MAX(10, 2*ne)) DO
                  IF re # NIL THEN SUBARRAY(se^, 0, ne) := re^ END;
                  re := se
                END
              END;
              e.xmark := TRUE;
              re[ne] := b; INC(ne)
            END
          END;
          b := Enext(b)
        UNTIL b = f;
      END
    END;
    (* Clear "xmark" bits: *)
    FOR k := 0 TO ne - 1 DO
      WITH b = re[k] DO b.facetedge.edge.xmark := FALSE END
    END;
  END CollectPolyhedronEdges;

PROCEDURE CollectPolyhedronVertices(
    a: Pair; 
    VAR rv: REF ARRAY OF Pair; 
    VAR nv: CARDINAL;
    VAR rf: REF ARRAY OF Pair;  (* Working area *)
  ) =
  VAR nf: CARDINAL;
      b: Pair;
  BEGIN
    CollectPolyhedronFaces(a, rf, nf);
    nv := 0;
    FOR k := 0 TO nf - 1 DO
      WITH f = rf[k] DO
        <* ASSERT Pneg(f) = Pneg(a) *>
        b := f;
        REPEAT
          WITH v = OrgV(b) DO 
            IF NOT v.xmark THEN 
              IF rv = NIL OR nv >= NUMBER(rv^) THEN
                WITH sv = NEW(REF ARRAY OF Pair, MAX(10, 2*nv)) DO
                  IF rv # NIL THEN SUBARRAY(sv^, 0, nv) := rv^ END;
                  rv := sv
                END
              END;
              v.xmark := TRUE;
              rv[nv] := b; INC(nv)
            END
          END;
          b := Enext(b)
        UNTIL b = f;
      END
    END;
    (* Clear "xmark" bits: *)
    FOR k := 0 TO nv-1 DO
      WITH b = rv[k] DO OrgV(b).xmark := FALSE END
    END;
  END CollectPolyhedronVertices;

PROCEDURE WriteTopology(
    name: TEXT;
    READONLY top: Topology;
    comments: TEXT := " ";
    ) =
  <* FATAL Wr.Failure, Thread.Alerted, OSError.E *>
  VAR pWidth : INTEGER;
  BEGIN
    IF top.NP = 0 THEN  pWidth:= 2 ELSE pWidth := Mis.NumDigits(top.NP-1) END;
    WITH 
      tp = FileWr.Open(name & ".tp"),
      vWidth = Mis.NumDigits(top.NV - 1),
      eWidth = Mis.NumDigits(top.NE - 1),
      fWidth = Mis.NumDigits(top.NF - 1),
      feWidth = Mis.NumDigits(top.NFE -1)
    DO
      WriteHeader(tp,"topology","99-08-25");
      IF NOT Text.Empty(comments) THEN
        WriteCommentsJS(tp, comments & "\n", '|')
      END;

      WITH
        m = Mis.NumDigits(top.NFE)
      DO
        Wr.PutText(tp, "vertices   ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top.NV),m) & "\n");  
        Wr.PutText(tp, "edges      ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top.NE),m) & "\n");
        Wr.PutText(tp, "faces      ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top.NF),m) & "\n"); 
        Wr.PutText(tp, "polyhedra  ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top.NP),m) & "\n");
        Wr.PutText(tp, "facetedges ");
        Wr.PutText(tp, Fmt.Int(top.NFE) & "\n");
        Wr.PutText(tp, "der        ");
        WITH der = 0 DO Wr.PutText(tp, Fmt.Pad(Fmt.Int(der),m) & "\n") END;
        Wr.PutText(tp, "bdr        ");
        WITH bdr = 0 DO Wr.PutText(tp, Fmt.Pad(Fmt.Int(bdr),m) & "\n") END;
      END;

      WriteCommentsJS(tp, "\nEdge data:\n", '|');
      FOR i := 0 TO top.NE-1 DO
        WITH e = top.edge[i] DO 
          Wr.PutText(tp, Fmt.Pad(Fmt.Int(e.num), eWidth)); 
          Wr.PutText(tp, "  ");
          Octf.PrintPair(tp, e.pa, eWidth+1, TRUE);  
        END
      END;

      WriteCommentsJS(tp, "\nFace data:\n", '|');
      FOR i := 0 TO top.NF-1 DO
        WITH f = top.face[i] DO 
          Wr.PutText(tp, Fmt.Pad(Fmt.Int(f.num), fWidth)); 
          Wr.PutText(tp, "  ");
          Octf.PrintPair(tp, f.pa, fWidth+1, TRUE);  
        END
      END;

      WriteCommentsJS(tp, "\nFacetEdge data:\n", '|');
      FOR i := 0 TO top.NFE-1 DO
        WITH fe = NARROW(top.facetedge[i].facetedge, FacetEdge) DO 
          Wr.PutText(tp, Fmt.Pad(Fmt.Int(fe.num), feWidth)); 
          Wr.PutText(tp, "  ");
          Octf.PrintFacetEdge(tp, fe, feWidth);
          Wr.PutText(tp, "  ");
          FOR j := 0 TO 3 DO
            WITH n = fe.org[j] DO 
              TYPECASE n OF 
              | NULL => FOR i:= 0 TO pWidth-2 DO Wr.PutText(tp," "); END;
                        Wr.PutText(tp, " - ");
              | Vertex(v) => Wr.PutText(tp, Fmt.Pad(Fmt.Int(v.num), vWidth)
                & "v ");
              | Polyhedron(p) => Wr.PutText(tp, Fmt.Pad(Fmt.Int(p.num),
                pWidth) & "p "); 
              ELSE (* nothing *)
              END; 
            END;
          END;
          Wr.PutText(tp, "  ");
	  WITH f = fe.face DO
            Wr.PutText(tp, Fmt.Pad(Fmt.Int(f.num), fWidth));
	  END;
          Wr.PutText(tp, "f ");
	  WITH e = fe.edge DO
            Wr.PutText(tp, Fmt.Pad(Fmt.Int(e.num), eWidth));
	  END;
          Wr.PutText(tp, "e\n");
        END
      END;
      WriteFooter(tp, "topology");
      Wr.Close(tp);
    END;     
  END WriteTopology;

PROCEDURE WriteDualTopology(
    name: TEXT;
    READONLY top: Topology;
    comments: TEXT := " ";
    ) =
  <* FATAL Wr.Failure, Thread.Alerted, OSError.E *>

  PROCEDURE PrintDualPair(wr: Wr.T; a: Pair; feWidth: CARDINAL) =
    BEGIN
      Wr.PutText(wr,Fmt.Pad(Fmt.Int(a.facetedge.num), feWidth) & ":" &
        Fmt.Int((Octf.SrotBits(a)+ 3) MOD 4) & ":" & Fmt.Int(Octf.SpinBit(a)));
    END PrintDualPair;  

  PROCEDURE PrintDualFacetEdge(wr: Wr.T; n: FacetEdge; feWidth: CARDINAL) =
    VAR b: Pair;
    BEGIN
      b := Srot(Pair{facetedge := n, bits := 0});
      FOR i := 0 TO 3 DO
        PrintDualPair(wr, Fnext(b), feWidth); Wr.PutText(wr, " "); 
        b := Srot(b)
      END
    END PrintDualFacetEdge;  

  BEGIN
    WITH 
      tp = FileWr.Open(name & ".tp"),
      vWidth = Mis.NumDigits(MAX(1,top.NV - 1)),
      eWidth = Mis.NumDigits(MAX(1,top.NE - 1)),
      fWidth = Mis.NumDigits(MAX(1,top.NF - 1)),
      pWidth = Mis.NumDigits(MAX(1,top.NP - 1)),
      feWidth = Mis.NumDigits(MAX(1,top.NFE -1))
    DO
      WriteHeader(tp,"topology","99-08-25");
      IF NOT Text.Empty(comments) THEN
        WriteCommentsJS(tp, comments & "\n", '|')
      END;

      WITH
        m  = Mis.NumDigits(top.NFE)
      DO
        Wr.PutText(tp, "vertices   ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top.NP),m) & "\n");  
        Wr.PutText(tp, "edges      ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top.NF),m) & "\n");
        Wr.PutText(tp, "faces      ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top.NE),m) & "\n"); 
        Wr.PutText(tp, "polyhedra  ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top.NV),m) & "\n");
        Wr.PutText(tp, "facetedges ");
        Wr.PutText(tp, Fmt.Int(top.NFE) & "\n");
        Wr.PutText(tp, "der        ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(1),m) & "\n");
        Wr.PutText(tp, "bdr        ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(0),m) & "\n");
      END;

      WriteCommentsJS(tp, "\nEdge data:\n", '|');
      FOR i := 0 TO top.NF-1 DO
        WITH f = top.face[i] DO 
          Wr.PutText(tp, Fmt.Pad(Fmt.Int(f.num), fWidth)); 
          Wr.PutText(tp, "  ");
          Octf.PrintPair(tp, f.pa, feWidth); (* [sic] *)
          Wr.PutText(tp, "\n") 
        END
      END;

      WriteCommentsJS(tp, "\nFace data:\n", '|');
      FOR i := 0 TO top.NE-1 DO
        WITH e = top.edge[i] DO 
          Wr.PutText(tp, Fmt.Pad(Fmt.Int(e.num), eWidth)); 
          Wr.PutText(tp, "  ");
          Octf.PrintPair(tp, e.pa, feWidth+1); (* [sic] *)
          Wr.PutText(tp, "\n")  
        END
      END;

      WriteCommentsJS(tp, "\nFacetEdge data:\n", '|');
      FOR i := 0 TO top.NFE-1 DO
        WITH fe = NARROW(top.facetedge[i].facetedge, FacetEdge) DO 
          Wr.PutText(tp, Fmt.Pad(Fmt.Int(fe.num), feWidth)); 
          Wr.PutText(tp, "  ");
          PrintDualFacetEdge(tp, fe, feWidth);
          Wr.PutText(tp, "  ");
          FOR j := 0 TO 3 DO
            WITH n = fe.org[(j+1) MOD 4] DO 
              TYPECASE n OF 
              | NULL => FOR i:= 0 TO pWidth-2 DO Wr.PutText(tp," "); END;
                        Wr.PutText(tp, " - ");
              | Vertex(v) => Wr.PutText(tp, Fmt.Pad(Fmt.Int(v.num), vWidth)
                & "p ");
              | Polyhedron(p) => Wr.PutText(tp, Fmt.Pad(Fmt.Int(p.num),
                pWidth) & "v "); 
              ELSE (* nothing *)
              END; 
            END;
          END;
          Wr.PutText(tp, "  ");
          WITH e = fe.edge DO
            Wr.PutText(tp, Fmt.Pad(Fmt.Int(e.num), eWidth));
          END;
          Wr.PutText(tp, "f ");
          WITH f = fe.face DO
            Wr.PutText(tp, Fmt.Pad(Fmt.Int(f.num), fWidth));
          END;
          Wr.PutText(tp, "e");
          Wr.PutText(tp, "\n");
        END
      END;
      WriteFooter(tp, "topology");
      Wr.Close(tp);
    END;     
  END WriteDualTopology;

PROCEDURE WriteStDe(
    wr: Wr.T;
    READONLY c: Coords;
    READONLY Dc: Coords;
    prec: CARDINAL := 4;
    comments: TEXT := "";
   ) =
    
  PROCEDURE WriteCoord(x: LONGREAL) =
    <* FATAL Wr.Failure, Thread.Alerted *>    
    BEGIN
      Wr.PutText(wr, Fmt.LongReal(x, Fmt.Style.Sci, prec := prec))
    END WriteCoord;

  PROCEDURE WritePoint(READONLY p: LR4.T) =   
    <* FATAL Wr.Failure, Thread.Alerted *>    
    BEGIN
      WriteCoord(p[1]);
      Wr.PutText(wr, " ");
      WriteCoord(p[0]);
      Wr.PutText(wr, " ");
      WriteCoord(p[2]);
      Wr.PutText(wr, " ");
      WriteCoord(p[3]);
      Wr.PutText(wr, " ");
    END WritePoint;

  <* FATAL Wr.Failure, Thread.Alerted *>    
  BEGIN
    WITH
      NV = NUMBER(c),
      d = Mis.NumDigits(NV-1)
    DO
      WriteHeader(wr,"state-derivatives","99-08-25");
      Mis.WriteCommentsJS(wr, "\n" & comments & "\n",'|');
      Wr.PutText(wr, "vertices = " & Fmt.Int(NV) & "\n");

      FOR i := 0 TO NV-1 DO 
        Wr.PutText(wr, Fmt.Pad(Fmt.Int(i), d) & ": "); 
        WritePoint(c[i]); 
        Wr.PutText(wr, "  ");

        (* if the derivatives of vertices are zero them writes "0 0 0 0"
           else writes velocites with the format for write points
           defined here.
        *)
        IF Dc[i][0] = 0.0d0 AND Dc[i][1] = 0.0d0 AND Dc[i][2] = 0.0d0 AND
           Dc[i][3] = 0.0d0
        THEN
          Wr.PutText(wr, "0 0 0 0")
        ELSE
          WritePoint(Dc[i])
        END; 
        Wr.PutText(wr, "\n");
      END;

      WriteFooter(wr, "state-derivatives");
      Wr.PutText(wr, "\n");
      Wr.Flush(wr);
    END;
  END WriteStDe;

BEGIN
END Triangulation.

(**************************************************************************)
(*                                                                        *)
(* Copyright (C) 2000     Universidade Estadual de Campinas (UNICAMP)     *)
(*                                                                        *)
(* Authors:                                                               *)
(*   L. P. 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.                        *)
(*                                                                        *)
(* Last edited on 2001-05-21 03:11:48 by stolfi                           *)
(**************************************************************************)