/* See {JSTri.h} */
#include <JSTri.h>

/* Formerly {JSTriangulation.c}; a junk version of {Triangulation.c}? */

/* Last edited on 2007-01-25 09:28:10 by stolfi */

#define JSTri_C_author \
  "Created 2000 by L. A. P. Lozada, based on 1996 procedures by\n" \
  "R. M. Rosi and J. Stolfi.\n" \
  "Revisions:\n" \
  "  03-08-2000 : Optimized version of the Read and Write procedures by J. Stolfi.\n" \
  "  30-08-2000 : Modified the Read and Write procedures for include the case\n" \
  "               when the cells are octahedra.\n" \
  "  07-10-2000 : Added procedure for compute the barycenter of a tetrahedron.\n" \
  "  27-01-2001 : Modified for exploding cubic cells."

#include <filefmt.h>
#include <r3.h>
#include <stdio.h>
#include <string.h>
#include <Mis.h>
#include <Octf.h>
#include <Triangulation.h>
#include <Random.h>
#include <filefmt.h>

void JSTriWedgeInit(Wedge_t w)
  { WedgeInit(&(w.w)); 
    w.org[0] = NULL; /* For now, node of primal: C  */
    w.org[1] = NULL; /* For now, node of dual:   C' */
    w.org[2] = NULL; /* For now, node of primal: C  */
    w.org[3] = NULL; /* For now, node of dual:   C' */
  }

  /* ===== ELEMENT CREATION ===== */
  
Place_t JSTriMakeWedge(void)
  { Wedge_t w = (Wedge_t)malloc(sizeof(WedgeRec_t)), "no mem"    wedge_init())
    JSTriWedgeInit(w);
    Place_t a = PlaceInWedge(e, 0);
    /* [!!! WAS:  !!!] w->edge->pa = a; */
    /* [!!! WAS:  !!!] w->wall->pa = a; */
    return a;
  }

  /* 
void SetOrgAll(Place_t @p, Node_t n)
  VAR
     Place_t t = a;
    Place_t tn;
  { do {
      SetOrg(t,n);
      tn = Clock(PrevE(t));
      do 
      {
        SetOrg(tn,n);
        tn = NextF(tn);
      } while (!(tn == Clock(PrevE(t))));
      t = NextF(t);       
    } while (t != a);
  } /* END SetOrgAll */
*/

/*
void SetOrgAll(Place_t @p, Node_t n)
       Place_t c = a;
  { ??? nei = Octf.RenumberEdgesOutOfNodes((Place_vec_t){a});
    {
      for(i = 0; i < nei.nel; i++)
      {
        ??? b = nei[i];
        {
          c = b;  
          do 
          {
            SetOrg(c,n);
            c = NextF(c)    
          } while (c != b);
        }
      }
    }
  } /* END SetOrgAll */

void SetOrgAll(Place_t @p, Node_t n)
  { SetOrgAll(a,n);
  } /* END SetOrgAll */

Node_t Pneg(Place_t @p)
  { return Org(Srot(a));
  } /* END Pneg */

void SetPneg(Place_t @p, Node_t n)
  { SetOrg(Srot(a), n)
  } /* END SetPneg */

void SetPnegOfWall(Place_t @p, Node_t n)
       Place_t t = a;
  { do {
      SetPneg(t, n);
      t = PrevE(t);
    } while (t != a);
  } /* END SetPnegOfWall */

void SetPnegOfNearbyWalls(Place_t @p, Node_t n)
       Place_t t = a;
  { SetPnegOfWall(t,n);
    do 
    {
      SetPnegOfWall(Clock(PrevF(t)),n);
      t = PrevE(t); 
    } while (t != a); 
  } /* END SetPnegOfNearbyWalls */

Node_t Ppos(Place_t @p)
  { return Pneg(Clock(a));
  } /* END Ppos */

void SetPpos(Place_t @p, Node_t n)
  { SetOrg(Tors(a), n)
  } /* END SetPpos */

void SetPposOfWall(Place_t @p, Node_t n)
       Place_t t = a;
  { do {  
      SetPpos(t, n);
      t = PrevE(t);
    } while (t != a);
  } /* END SetPposOfWall */

void SetPposQUX(Place_t @p, Node_t n)
       Place_t t = Clock(a);
  { SetPnegOfNearbyWalls(t,n);
    do 
    {
      SetPnegOfNearbyWalls(PrevF(t),n);
      t = PrevE(t); 
    } while (!(t == Clock(a))); 
  } /* END SetPposQUX */

Node OrgV(Place_t @p)
  { return (Node_t)(Org(a));
  } /* END OrgV */

Node DesV(Place_t @p)
  { return OrgV(Clock(a));
  } /* END DesV */

Cell PnegP(Place_t @p)
  { return (Cell_t)(Pneg(a));
  } /* END PnegP */

Cell PposP(Place_t @p)
  { return PnegP(Clock(a));
  } /* END PposP */

FourNodes_t TetraNegNodes(Place_t @p)
  /* Valid for the versions Spin(a), Clock(a) and SpinClock(a).  */
  { assert(Pneg(a)!=NULL);
    Node_t p = OrgV(a);
    Node_t q = OrgV(NextE(a));
    Node_t r = OrgV(PrevE(a));
    Node_t s = OrgV(PrevE(PrevF(a)));
    {
      assert(Pneg(a) == Ppos(PrevE(PrevF(a))));
      assert(p != q) && (q != r) && (r != s);
      return (FourNodes_t){{p,q,r,s}};
    }
  } /* END TetraNegNodes */

FourWalls_t TetraWalls(Place_t @p)
  { assert(Pneg(a)!=NULL);
    Wall_t f0 = PWall(a);
    Wall_t f1 = PWall(PrevF(a));
    Wall_t f2 = PWall(PrevF(NextE(a)));
    Wall_t f3 = PWall(PrevF(PrevE(a)));
    {
      assert(f0 != f1) && (f1 != f2) && (f2 != f3); 
      return (FourWalls_t){{f0,f1,f2,f3}};
    }
  } /* END TetraWalls */

SixEdges_t TetraEdges(Place_t @p)
  { assert(Pneg(a)!=NULL);
    Edge_t e0 = PEdge(a);
    Edge_t e1 = PEdge(NextE(a));
    Edge_t e2 = PEdge(PrevE(a));
    Edge_t e3 = PEdge(NextE(PrevF(a)));
    Edge_t e4 = PEdge(PrevE(PrevF(a)));
    Edge_t e5 = PEdge(NextE(PrevF(NextE(a))));
    {
      assert(e0#e1) && (e1#e2) && (e2#e3) && (e3#e4) && (e4#e5); 
      return (SixEdges_t){{e0,e1,e2,e3,e4,e5}};
    }
  } /* END TetraEdges */

ThreeEdges_t TriangEdges(Place_t @p)
  { Edge_t e0 = PEdge(a);
    Edge_t e1 = PEdge(NextE(a));
    Edge_t e2 = PEdge(PrevE(a));
    {
      assert(e0 != e1) && (e1 != e2); 
      return (ThreeEdges_t){{e0,e1,e2}};
    }
  } /* END TriangEdges */
  
FourNodes_t TetraPosNodes(Place_t @p)
  /* Valid for the versions Spin(a), Clock(a) and SpinClock(a). */
  { assert(Ppos(a)!=NULL);
    Node_t p = OrgV(a);
    Node_t q = OrgV(NextE(a));
    Node_t r = OrgV(PrevE(a));
    Node_t s = OrgV(PrevE(NextF(a)));
    {
      assert(Ppos(a) == Pneg(PrevE(NextF(a))));
      return (FourNodes_t){{p,q,r,s}};
    }
  } /* END TetraPosNodes */

FiveNodes_t TetraNegPosNodes(Place_t @p) ==
  /* 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. */
  { ??? p = Org(a);
    ??? q = Org(NextE(a));
    ??? r = Org(PrevE(a));
    ??? s = Org(PrevE(PrevF(a)));
    ??? t = Org(PrevE(NextF(a)));
    {
      assert(Ppos(a) == Pneg(PrevE(NextF(a))) AND
                Pneg(a) == Ppos(PrevE(PrevF(a)))
     );
      return (FiveNodes_t){p,q,r,s,t};
    }
  } /* END TetraNegPosNodes */

  /* ========================== CONSTRUCTION TOOLS ===============  */


EightPlaces_t MakeTetraTopo(uint nx, uint ny)
  VAR
    uint wedgeCount = 0; 
    uint cellCount = 0;

  Place_t MakeTriangle(void)
    /* Make one triangular wall and set of the three places with the
       same wall component.
    */
  { Place_t a = MakeWedge();
    Place_t b = MakeWedge();
    Place_t c = MakeWedge();
    Wall_t f = PWall(a);
    Node_t u = MakeNode();
    Node_t v = MakeNode();
    Node_t w = MakeNode();
    {
      PWedge(a)->num = wedgeCount; wedgeCount++;         
      SetOrg(a, u); SetOrg(Clock(a),v);

      PWedge(b)->num = wedgeCount; wedgeCount++;
      SetNextE(a,b); SetWallInfo(b,f);
      SetOrg(b,v); SetOrg(Clock(b),w);

      PWedge(c)->num = wedgeCount; wedgeCount++;
      SetNextE(b,c); SetWallInfo(c,f);
      SetOrg(c, w);  SetOrg(Clock(c), Org(a));
      return a;
    }
  } /* END MakeTriangle */

  Place_t AddTwoTriangles(Place_t @p)
    /*
      Build a new tetrahedral cell by insertion of two triangular walls and
      operations SpliceWalls. The argument "a" is the @place return by
      MakeTriangle(). The @place "f" is return by MakeCell().
    */
    Place_t *c;
    {    
      c = NextE(NextF(a));
      SetEdgeInfo(a, PEdge(PrevE(c)));
      ??? f = MakeTriangle(), g == MakeTriangle();
    {

        SetNextF(PrevE(a), PrevE(f));
        SetEdgeInfo(PrevE(f), PEdge(PrevE(a)));
        SetOrgAll(PrevE(a), Org(PrevE(a)));
        SetOrgAll(Clock(PrevE(a)), Org(Clock(PrevE(a))));

        SetNextF(Clock(f), NextE(c));
        SetEdgeInfo(Clock(f), PEdge(NextE(c))); 
        SetOrgAll(NextE(c), Org(NextE(c)));
        SetOrgAll(Clock(NextE(c)), Org(Clock(NextE(c))));     
        
        SetNextF(NextE(g), NextE(f));
        SetEdgeInfo(NextE(g), PEdge(NextE(f)));
        SetOrgAll(NextE(f), Org(NextE(f)));
        SetOrgAll(Clock(NextE(f)), Org(Clock(NextE(f))));
        
        SetNextF(g, c);
        SetEdgeInfo(g, PEdge(c));
        SetOrgAll(c, Org(c));
        SetOrgAll(Clock(c), Org(Clock(c))); 
        
        SetNextF(PrevE(g), Clock(NextE(a)));
        SetEdgeInfo(PrevE(g), PEdge(Clock(NextE(a)))); 
        SetOrgAll(NextE(a), Org(NextE(a)));
        SetOrgAll(Clock(NextE(a)), Org(Clock(NextE(a))));
                 
        Cell_t p = MakeCell();
        {
          p->num = cellCount;
          cellCount++;
          SetPnegOfNearbyWalls(f,p);
        }
        return f;
      }
  } /* END AddTwoTriangles */

  Place_t FinishCellRow(Place_t @p)
    /*
      Adds one new row of tetradedral cells and return the @place wedge
      belong to the cell more rigth.
    */
   {
     for (col = 0; col < nx; col++)
     {
       a = AddTwoTriangles(a);
     }
     return a;
   } FinishCellRow;

  VAR
    Place_t t,b;
   EightPlaces_t ca;
  {

  /* ======================= create bottom row of triangles ================= */

      for (col = 0; col < nx; col++)
      {
        Place_t c = MakeTriangle();
    {
          if (col == 0)
          {
            t = c;
          }
          else
          {
          SpliceWalls(NextE(b), Clock(PrevE(c)));
          SetEdgeInfo(NextE(b), PEdge(Clock(PrevE(c))));
          SetOrgAll(NextE(b), Org(NextE(b))); 
          SetOrgAll(Clock(NextE(b)), Org(Clock(NextE(b))));
          }
          b = c;
        }
      }

      ca.p[2] = t;  ca.p[7] = Spin(Clock(b));
      
      for (row = 0; row < ny; row++)
      {
        Place_t a = MakeTriangle();
    {
          if (row == 0){ ca.p[1] = Clock(PrevE(a)); }           
          if (row == ny-1){ ca.p[4] = Spin(PrevE(a)); }
          
          SpliceWalls(Clock(a), Clock(PrevE(t)));
          SetOrgAll(PrevE(t), Org(PrevE(t)));
          SetOrgAll(Clock(PrevE(t)), Org(Clock(PrevE(t))));

           Place_t aa = FinishCellRow(a);
    {      
             t = PrevF(t);
             SetOrg(PrevE(t), Org(PrevE(a)));
             b = PrevF(b);
             SetOrg(Clock(NextE(b)), Org(PrevE(aa)));
             if (row == 0){ ca.p[0] = Clock(Spin(PrevE(aa))); } 
             if (row == ny-1){ ca.p[5] = PrevE(aa); }
           }
        }
      }
     ca.p[3] = Spin(t);
     ca.p[6] = Clock(b);
    return ca;
  } /* END MakeTetraTopo */         

void EmphasizeTetrahedron(Place_t @p, Place_t b, uint n)

   void HideNode(Node_t v)
    /* Hide the node "v". */ 
    {
      v->exists = FALSE;
    } HideNode;

   void HideEdge(Place_t @p)
    /* Hide the edge "PEdge(a)".*/ 
    {
      Edge_t e = PEdge(a);
      {  
        e->exists = FALSE;
      }        
    } HideEdge;

   void HideWall(Place_t @p)
    /* Hide the wall "PWall(a)".*/ 
    {
      Wall_t f = PWall(a)
      {  
        f->exists = FALSE;
      }        
    } HideWall;

   void HideWallsOfEdge(Place_t @p)
    /*  Hide the ring wall "PWall(a)". */
    Place_t an = PrevF(a);
    VAR {
      for (j = 1; j < n; j++)
      {
        HideWall(an);
        an = PrevF(an);
      }
    } HideWallsOfEdge;

  VAR
    ta,Place_t tb[100+1];
  { 
    ta[0] = a; tb[0] = b;

    for (i = 1; i < n; i++)
    {
      ta[i] = Clock(PrevE(NextF(NextE(ta[i-1]))));
      tb[i] = Clock(PrevE(NextF(NextE(tb[i-1]))));
    }

   
    for (i = 0; i < n; i++)
    {
      HideWallsOfEdge(ta[i]);
      HideWallsOfEdge(tb[i]);
    }

    for (i = 1; i < n; i++)
    {
      HideNode(OrgV(tb[i])); 
      HideNode(OrgV(ta[i]));
      Place_t dn = tb[i];
      VAR {
        for (j = 0; j <= (n); j++)
        {
          HideEdge(PrevE(dn));
          dn = PrevF(dn);
        }
      }
    }

    for (i = 0; i <= (n-2); i++)
    {
      HideEdge(NextE(ta[i]));
      HideEdge(NextE(NextF(ta[i])));
    }
  } /* END EmphasizeTetrahedron */  

Place_t Glue(
  Place_t @p,b;
    uint n;
  bool_t setorg = TRUE;
 )
  /* The @place "a" and "b" have the same Orientation and Spin bits,  such as,
     after of the glue procedure performs: Pneg(a) == Pneg(b). */

  VAR
    ta,Place_t tb[100+1];
  { /* sanity check */
    assert(n >= 1);
    ta[0] = a; tb[0] = b;

    if (n > 1)
    {
      for (i = 1; i < n; i++)
      {
        ta[i] = Clock(PrevE(PrevF(NextE(ta[i-1]))));
        tb[i] = Clock(PrevE(NextF(NextE(tb[i-1]))));
        assert(ta[i]!=a);
        assert(tb[i]!=b);  
      }
    }
      
    Meld(b, a);
    
    /* Update the edge slots for i==0 */
    
    SetRingEdgeInfo(a,          PEdge(a));
    SetRingEdgeInfo(NextE(a),   PEdge(NextE(a)));
    SetRingEdgeInfo(PrevE(a), PEdge(PrevE(a)));
    
    if (setorg)
    {
      /* Update the node slots for i==0 */
      SetOrgAll(a,        Org(a));
      SetOrgAll(Clock(a), Org(Clock(a)));
      SetOrgAll(NextE(a), Org(NextE(a)));
      SetOrgAll(Clock(NextE(a)), Org(Clock(NextE(a))));
      SetOrgAll(PrevE(a), Org(PrevE(a)));
      SetOrgAll(Clock(PrevE(a)), Org(Clock(PrevE(a))));
    }
    
    /* Update the cell slots for i==0 */
    SetPneg(a, Pneg(b));
    SetPneg(PrevE(a), Pneg(PrevE(b)));
    SetPneg(NextE(a), Pneg(NextE(b)));
    for (i = 1; i < n; i++)
    {
      Meld(tb[i],ta[i]);
      
      /* Update the edge slots */
      SetRingEdgeInfo(ta[i], PEdge(ta[i]));
      SetRingEdgeInfo(NextE(ta[i]), PEdge(NextE(ta[i])));
      SetRingEdgeInfo(PrevE(ta[i]), PEdge(PrevE(ta[i])));

      if (setorg)
      {
        /* Update the node slots */
        SetOrgAll(ta[i], Org(ta[i]));
        SetOrgAll(Clock(ta[i]), Org(Clock(ta[i])));
        SetOrgAll(NextE(ta[i]), Org(NextE(ta[i])));
        SetOrgAll(Clock(NextE(ta[i])), Org(Clock(NextE(ta[i]))));
        SetOrgAll(PrevE(ta[i]), Org(PrevE(ta[i])));
        SetOrgAll(Clock(PrevE(ta[i])), Org(Clock(PrevE(ta[i]))));
      }

      /* Update the cell slots */
      ??? f = Pneg(tb[i]);
        ??? g = Pneg(PrevE(tb[i]));
      ??? h = Pneg(NextE(tb[i]));
      {
        SetPneg(ta[i],f);
        SetPneg(PrevE(ta[i]),g);
        SetPneg(NextE(ta[i]),h);
      }
    }

    if (setorg)
    {
      SetOrgAll(Clock(PrevE(PrevF(PrevE(ta[n-1])))),Org(PrevE(a)));
    }
    return ta[n-1];
  } /* END Glue */

  /* ===== GLOBAL PROCEDURES ===== */

VAR
  Node_vec_t seenNode = Node_vec_new(INIT_STACK_SIZE);
  uint_vec_t nodeOnum = uint_vec_new(INIT_STACK_SIZE);
  uint seenNodeCount = 0;

void MarkNode(Node_t n)
  { Node_vec_expand(&seenNode,seenNodeCount);
    uint_vec_expand(&nodeOnum,seenNodeCount);
    seenNode[k] = seenNodeCount;
    nodeOnum[k] = seenNodeCount->num;
    n->num = seenNodeCount;
    seenNodeCount++
  }

void NodeIsMarked(Node_t n) : bool_t 
  { return (n->num < seenNodeCount)) && ((seenNode[n->num] == n);
  } /* END @{Node->?}IsMarked */

void UnmarkMarkedNodes()
  { while (seenNodeCount > 0){
      ??? k = seenNodeCount-1, n == seenNode[k];
    {
        n->num = nodeOnum[k];
        seenNodeCount--;
      }
    }
  } /* END UnmarkMarked@{Node->?}s */

void EnumNodes(Place_t @p, visit: VisitProc)
  CONST INIT_STACK_SIZE == 1000;
  VAR
    Place_vec_t pstack = Place_vec_new(INIT_STACK_SIZE);
    uint top;    /* top of pstack */

  void Stack(Place_t c)
    Place_t cn = c;
        VAR Place_t dn;
    {
      do 
      {
        Place_vec_expand(&pstack,top);
        pstack.el[top] = cn;
        top++;
        dn = Clock(PrevE(cn));
        do 
        {
          Place_vec_expand(&pstack,top);
          pstack.el[top] = dn;
          top++;
          dn = NextF(dn);
        } while (!(( dn == Clock(PrevE(cn)))));
        cn = NextF(cn);
      } while (cn != c);
    } Stack;

  void VisitAndMark(Place_t c)
    /* If org(c) is diferent that the origins @places stacked, stack the @place
       "c".
    */
    Place_t *cn;
    {
      ??? n = Org(c);
    {
        if ((n!=NULL) && (NOT NodeIsMarked(n)))
        {
          visit(c);
          MarkNode(n);
          Stack(c);
          cn = c;
          do 
          {
            Stack(cn);
            cn = Onext(cn);
          } while (!( (cn == c)));
        }
      }        
    } VisitAndMark;

  uint *seen;
  { top = 0;
    seen = 0;
    assert(seenNodeCount == 0);
    VisitAndMark(a);
    while (seen < top)
    {
      Place_t b = (Place_t){festack.el[seen], bstack.el[seen]};
      {
        VisitAndMark(NextF(b));
        VisitAndMark(NextF(NextE(b)));
        VisitAndMark(NextE(b));
        VisitAndMark(PrevF(PrevE(b)));
      }
      seen++;
    }
    UnmarkMarkedNodes();
  } /* END EnumNodes */


uint NumberNodes(Place_t @p)
      uint n = 0;
  
  void Visit(Place_t c)
  
    void Visit_(Place_t c)
      Place_t *cn,dn;  
      {
        ??? v = Org(c);
    {
          cn = c;
          do 
          {
            ??? vn = Org(cn);
    {
              assert(vn == v);
              vn->num = n;
              dn = Clock(PrevE(cn));
              do 
              {
                ??? wn = Org(dn);
    {
                  assert(wn == vn);
                  wn->num = n;
                  dn = NextF(dn);
                }
              } while (!(( dn == Clock(PrevE(cn)))));
              cn = NextF(cn);
            }
          } while (!( cn == c));
        }
      } Visit_;

    Place_t cn = c;
    VAR {
      do 
      {
        Visit_(cn);
        cn = Onext(cn);
      } while (cn != c);
      n++;
   } Visit;

   void VisitDual(Place_t c)
     void VisitDual_(Place_t c)
       VAR
         Place_t cn;
         p,pn: Node_t;
       {
         p = Org(c);
         if (p != NULL)
         {
           cn = c;
           do 
           {
             pn = Org(cn);
             if (pn != NULL)
             {
               assert(pn == p); 
               pn->num = n;
             }
             cn = NextF(cn);
           } while (!( cn == c));
         }
       } VisitDual_;

     Place_t cn = c;
     VAR {
       VisitDual_(cn);
       VisitDual_(Clock(PrevE(cn)));
       n++;
     } VisitDual;

    {
      if ((Octf.DualBit(a) == 0))
      {
        EnumNodes(a, Visit);
      }
      else
      {
        EnumNodes(a, VisitDual);
      }
    return n;
  } /* END NumberNodes */


Topology MakeTopology(Place_t @p, uint bdr)
  Topology top ;
    INTEGER euler;
  { top->NV         = NumberNodes(a);
    top->node     = Node_vec_new(top->NV); 
    top->wedge  = Octf.RenumberWedges((Place_vec_t){a});  
    top->NFE        = ((top->wedge^).nel);
    top->edge       = Octf.RenumberEdges((Place_vec_t){top->wedge[0]});
    top->wall       = Octf.RenumberWalls((Place_vec_t){top->wedge[0]});
    top->out        = Place_vec_new(top->NV);
    top->NE         = ((top->edge^).nel);
    top->NF         = ((top->wall^).nel);
    top->NP         = NumberNodes(Srot(a));
    top->cell = NEW(REF ARRAY OF Cell, top->NP);
    top->region     = Place_vec_new(top->NP);
    top->der        = Octf.DegreeOfEdge(top->wall[0].pa);
    top->bdr        = bdr;

    for (i = 0; i < top->NFE; i++)
    {
      VAR
         Place_t c = top->wedge[i];
        Node_t v,p;
        uint vi,qi;
      {
        for (k = 0; k < 2; k++)
        {
          v = Org(c); 
          vi = v->num;
          top->node[vi] = v;
          top->out[vi] = c;
          p = Pneg(c);
          if (p != NULL)
          {
            qi = p->num;
            top->cell[qi] = p;
            top->region[qi] = Srot(c);
            assert(Pneg(c) == Org(top->region[qi]));
          }
          c = Clock(c); 
        }
      }
    }
    
    ??? m = NumDigits(top->NFE)+1;
    {
      fprintf(stderr, "\n");
      fprintf(stderr, "nv  = " & Fmt.Pad(Fmt.Int(top->NV),m)  & "\n"); 
      fprintf(stderr, "ne  = " & Fmt.Pad(Fmt.Int(top->NE),m)  & "\n");
      fprintf(stderr, "nf  = " & Fmt.Pad(Fmt.Int(top->NF),m)  & "\n");
      fprintf(stderr, "np  = " & Fmt.Pad(Fmt.Int(top->NP),m)  & "\n");
      fprintf(stderr, "nfe = " & Fmt.Pad(Fmt.Int(top->NFE),m) & "\n");
    }

    /* Euler's Number */
    euler = top->NV-top->NE+top->NF-top->NP;
    if (euler == 0)
    {
      fprintf(stderr,"\nThe Map performs the Euler's Number == 0\n");
    }
    else
    {
      fprintf(stderr,"\nThe Euler's Number: " & Fmt.Int(euler) & "\n");
    }
    return top;
  } /* END MakeTopology */

uint DegreeOfNode(Place_t @p)
  { ??? @{edge->?} = Octf.RenumberEdgesOutOfNodes((Place_vec_t){a});
    ??? degree = ((@{edge->?}^).nel);
    {
      return degree;
    }
  } /* END DegreeOfNode */

AdjacencyMatrix *MakeAdjacencyMatrix(Topology *top)==
      m = NEW(REF ARRAY OF ARRAY OF bool_t, top->NV, top->NV);
  { for (i = 0; i < top->NV; i++){ 
      ??? mi = m[i];
    { 
        for (j = 0; j < top->NV; j++)
        {
          mi[j] = FALSE
        }
      }
    }
    for (ei = 0; ei < top->NFE; ei++)
    {
      ??? a = top->wedge[ei];
    {
        ??? i = Org(a)->num, j == Org(Clock(a))->num;
    { 
          m[i,j] = TRUE;
          m[j,i] = TRUE;
        }
      }
    }
    return m;
  } /* END MakeAdjacencyMatrix */

bool_t TriviallyIsomorphic(Topology_t *ta, Topology_t *tb)
  { if ((ta.NV!=tb.NV
    ) || (ta.NF!=tb.NF 
    ) || (ta.NE!=tb.NE
    ) || (ta.NP!=tb.NP 
    ) || (ta.NFE!=tb.NFE
  )){ return FALSE; }
    int NFE = ta.NFE;
    {
      for (i = 0; i < NFE; i++)
      {
        VAR sPlace_t @p = PWedge(ta)->i]; sb: Place_t = PWedge(tb)->i];
        {
          for (r = 0; r < 4; r++)
          {
            ??? na = Org(sa)->num;
            ??? nb = Org(sb)->num;
            {
              if (na != nb){ return FALSE;}
            }
            ??? za = GetPlaceNum_t(NextF(sa));
              ??? zb = GetPlaceNum_t(NextF(sb));
            {
              if (za != zb){ return FALSE;}
            }
            sa = Srot(sa);
            sb = Srot(sb);
          }
        }
      }
    }
    return TRUE;
  } /* END TriviallyIsomorphic */

void CheckOutAndRegion(Topology *top)
  { for (i = 0; i < top->NV; i++){ 
      ??? v = top->node[i];
        ??? e = top->out[i];
      {
        assert(v->num == i);
        assert(Org(e) == v);
      }
    }
    for (i = 0; i < top->NP; i++)
    {
      ??? p = top->cell[i];
        ??? r = top->region[i];
      {
        assert(p->num == i);
        assert(Org(r) == p);
      } 
    }
  } /* END CheckOutAndRegion */

void GetVariableNodes(Topology_t *top,
                              VAR vr: ARRAY OF bool_t)
  { for (i = 0; i < top->NV; i++){
      ??? v = top->node[i];
    {
        vr[i] = NOT v.fixed; 
      }
    }
  } /* END GetVariableNodes */
  
  /* ===== GEOMETRIC TOOLS ===== */

void InitCoords(
    rand_t *coins; 
    Coords_t *c; 
    float radius = 1.0)
  { ??? r = ((double)radius);
    {
      for (i = 0 TO (c.nel-1))
      {
        c[i] = (r4_t)
        {
          coins.longreal(-r, r),
          coins.longreal(-r, r),
          coins.longreal(-r, r),
          coins.longreal(-r, r)
        }
      }
    }
  } /* END InitCoords */

Coords_t *GenCoords(Topology_t *t)
  { ??? coins = NEW(Random.Default).init(TRUE);
    ??? r = r4_vec_new(t.NV);
    ??? c = r^;
    {
      for (i = 0 TO (c.nel-1))
      {
        c[i] = (r4_t)
        {
                      coins.longreal(-1.0, +1.0),
                      coins.longreal(-1.0, +1.0),
                      coins.longreal(-1.0, +1.0),
                      coins.longreal(-1.0, +1.0)
                     }
      }
      return r;
    }
  } /* END GenCoords */

r4_t Barycenter(Topology_t *top, Coords_t *c)
      B: r4_t = (r4_t){0.0, ..};
      uint n = 0;
  { for (i = 0 TO (c.nel-1)){ 
      ??? v = top->node[i];
    {
        if (v->exists){ B = r4_Add(B, c[i]); N++; }
      }
    }
    return r4_Scale(1.0/((double)N), B);
  } /* END Barycenter */

double MeanNodeDistance(Topology_t *top,
                         *c: Coords_t)
      S: double = 0.0;
      uint n = 0;
  { for (i = 0 TO (c.nel-1)){ 
      ??? v = top->node[i];
    {
        if (v->exists)
        {
          S = S + r4_NormSqr(c[i]);
          N++
        }
      }
    }
    return Math.sqrt(S/FLOAT(N,double));
  } /* END MeanNodeDistance */

double MeanEdgeLength(Topology_t *top,
                         *c: Coords_t)
      S: double = 0.0;
      uint n = 0;
  { for (i = 0; i < top->NE; i++){ 
      ??? e = top->edge[i];
    {
        if (e->exists)
        {
          ??? o = e.node[0]->num, d == e.node[1]->num;
    {
            S = S + r4_DistSqr(c[o], c[d])
          }
          N++;
        }
      }
    }
    return Math.sqrt(S/FLOAT(N,double));
  } /* END MeanEdgeLength */

void Displace(Topology_t *top, r4_t *d, Coords_t *c)

  { for (i = 0 TO (c.nel-1)){ 
      if (top->node[i]->exists)
      {
        ??? vc = c[i];
    { vc = r4_Add(vc, d);}
      }
    }
  } /* END Displace */

void Scale(Topology_t *top, double s, Coords_t *c)
  { for (i = 0 TO (c.nel-1)){ 
      if (top->node[i]->exists)
      {
        ??? vc = c[i];
    { vc = r4_Scale(s, vc); }
      }
    }
  } /* END Scale */

void NormalizeNodeDistance(Topology_t *top, Coords_t *c)
  { ??? b = Barycenter(top, c);
    { Displace(top, r4_Neg(b), c); }
    ??? s = MeanNodeDistance(top, c);
    { Scale(top, 1.0/s, c); }
  } /* END NormalizeNodeDistance */

void NormalizeEdgeLengths(Topology_t *top, Coords_t *c)
  { ??? b = Barycenter(top, c);
    { Displace(top, r4_Neg(b), c); }
    ??? s = MeanEdgeLength(top, c);
    { Scale(top, 1.0/s, c); }
  } /* END NormalizeEdgeLengths */

r4_t WallCross(Place_t @p, Coords_t *c)
  { if ((! PWall(a)->exists)){
      return (r4_t){1.0, 0.0, 0.0, 0.0};
    }
    else
    {
      Node_t ov = OrgV(a);
        Node_t dv = OrgV(Clock(a));
      Node_t pv = OrgV(PrevE(a));
      Node_t qv = OrgV(PrevE(NextF(a)));
      Node_t rv = OrgV(PrevE(PrevF(a)));
      {
        ??? o = c[ov->num];
          ??? d = c[dv->num];
        ??? p = c[pv->num];
        ??? q = c[qv->num];
        ??? r = c[rv->num];
        ??? n1 = LR4Extras.Cross(r4_Sub(p,o), r4_Sub(r,o), r4_Sub(d,o));
        ??? n2 = LR4Extras.Cross(r4_Sub(p,o), r4_Sub(d,o), r4_Sub(q,o));
        { 
          return r4_Add(n1,n2);
        }
      }
    }
   } WallCross;

r4_t WallNormal(Place_t @p, Coords_t *c)
  { ??? n = WallCross(a,c);
    ??? s = r4_Norm(n);
    {
      if (s == 0.0)
      {
        return (r4_t){1.0, 0.0, 0.0, 0.0};
      }
      else
      {
        return r4_Scale(1.0/s, n);
      }
    }
  } /* END WallNormal */

r4_t PolyCross(Place_t @p, Coords_t *c)
  { if ((! PnegP(a)->exists)){
      return (r4_t){1.0, 0.0, 0.0, 0.0};
    }
    else
    {
      Node_t ov = OrgV(a);
        Node_t dv = OrgV(Clock(a));
      Node_t pv = OrgV(PrevE(a));
      Node_t rv = OrgV(PrevE(PrevF(a)));
      {
        ??? o = c[ov->num];
          ??? d = c[dv->num];
        ??? p = c[pv->num];
        ??? r = c[rv->num];
        {
          return LR4Extras.Cross(r4_Sub(p,o), r4_Sub(r,o), r4_Sub(d,o));
        }
      }
    }
   } PolyCross;

r4_t PolyNormal(Place_t @p, Coords_t *c)
  { ??? n = PolyCross(a, c);
    ??? s = r4_Norm(n);
    {
      if (s == 0.0)
      {
        return (r4_t){1.0, 0.0, 0.0, 0.0};
      }
      else
      {
        return r4_Scale(1.0/s, n);
      }
    }
  } /* END PolyNormal */
  
r4_t WallBarycenter(Place_t @p, Coords_t *c)
  VAR
    Place_t @po;
    n : uint = 0;
    sum = (r4_t){0.0, ..};
  { ao = a;
    do 
    {
      ??? aoc = c[OrgV(ao)->num];
      { 
        sum = r4_Add(sum, aoc);
        n++;
        ao = NextE(ao);
      }
    } while (ao != a);
    if (n == 0)
    {
      return sum;
    }
    else
    {
      return r4_Scale(1.0/FLOAT(n,double), sum);
    }
  } /* END WallBarycenter */

r4_t TetraBarycenter(Place_t @p, Coords_t *c)
  uint n = 0;
    sum = (r4_t){0.0, ..};
  { ??? tetra = TetraNegNodes(a);
    {
      for (i = 0; i < 4; i++)
      {
        ??? aoc = c[tetra[i]->num];
       {        
          sum = r4_Add(sum, aoc);
          n++;
        }
      }
    }
    return r4_Scale(1.0/FLOAT(n,double), sum);
  } /* END TetraBarycenter */

r4_t EdgeCross(Place_t @p, Coords_t *c)
      sum: r4_t = (r4_t){0.0,..};
      Place_t @po;
  { Node_t uv = OrgV(a);
    ??? u = c[uv->num];
    {
        if ((! PEdge(a)->exists))
        {
          return (r4_t){1.0, 0.0, 0.0, 0.0};
        }
        else
        {
          ao = a;
          do 
          {
            Place_t an = PrevF(ao);
              Node_t dv = OrgV(Clock(ao));
            Node_t pv = OrgV(PrevE(ao));
            Node_t rv = OrgV(PrevE(an));
            {
              ??? d = c[dv->num];
                ??? p = c[pv->num];
              ??? r = c[rv->num];
              ??? n = LR4Extras.Cross(r4_Sub(p,u), r4_Sub(d,u), r4_Sub(r,u));
              {
                if (ao == a){ sum = n; }else{ sum = r4_Add(sum,n); }
              }
              ao = an;
            }
          } while (ao != a);
          return sum;
        }
    }
  } /* END EdgeCross */

r4_t EdgeNormal(Place_t @p, Coords_t *c)
  { ??? n = EdgeCross(a, c);
    ??? s = r4_Norm(n);
    {
      if (s == 0.0)
      {
        return (r4_t){1.0, 0.0, 0.0, 0.0};
      }
      else
      {
        return r4_Scale(1.0/s, n);
      }
    }
  } /* END EdgeNormal */

r4_t NodeCross(Place_t @p, Coords_t *c,
                      *top : Topology)
      
    sum,n: r4_t = (r4_t){0.0, ..};
  { Node_t uv = OrgV(a);
    ??? poly = StarOfNode(a,top);
    {
      if (! uv->exists)
      {
        return (r4_t){1.0, 0.0, 0.0, 0.0};
      }
      else
      {
        for(i = 0; i < poly.nel; i++)
        {
          if ((poly[i][0]->exists) && (poly[i][1]->exists AND
             poly[i][2]->exists) && (poly[i][3]->exists
        )){
            ??? u = c[poly[i][0]->num];
              ??? d = c[poly[i][1]->num];
            ??? p = c[poly[i][2]->num];
            ??? r = c[poly[i][3]->num];
            {
              n= LR4Extras.Cross(r4_Sub(p,u), r4_Sub(d,u), r4_Sub(r,u));
              if (i == 0){ sum = n; }else{ sum = r4_Add(sum,n); } 
            }
          }
        }
      }
      return sum;
    }
  } /* END NodeCross */

r4_t NodeNormal(Place_t @p, Coords_t *c, Topology_t *top)
  { ??? n = NodeCross(a, c, top);
    ??? s = r4_Norm(n);
    {
     if (s == 0.0)
     {
        return (r4_t){1.0, 0.0, 0.0, 0.0};
      }
      else
      {
        return r4_Scale(1.0/s, n);
      }
    }
  } /* END NodeNormal */
  
r4_t NeighborBarycenter(Node_vec_t *neigh;
                             *c: Coords_t)
  uint n = 0;
    b = (r4_t){0.0, ..};
  { ??? r =  ((neigh^).nel);
    {
      for (i = 0; i < r; i++)
      {
        ??? v = neigh^[i];
    {
          if (v->exists){ b = r4_Add(b, c[v->num]); n++;}
        }
      }
    }

    if (n == 0)
    {
      return b;
    }
    else
    {
      return r4_Scale(1.0/((double)n), b)
    }
  } /* END NeighborBarycenter */

Node_vec_t *NeighborNode(Place_t @p, Topology_t *top) ==
  
  REF ARRAY OF Triv tri ;
    Node_vec_t vstack = Node_vec_new(INIT_STACK_SIZE);
    uint nstack = 0;
    uint c;

    bool_t Present(Node_t v)
      /* Return TRUE if "v" is on the stack, FALSE c.c */
      uint nstack1 = nstack;
      VAR {
        while (nstack1 > 0)
        {
          nstack1 = nstack1 - 1;
          if (vstack.el[nstack1] == v){ return TRUE; }
        }
        return FALSE;
      } Present;

  { ??? poly = StarOfNode(a,top);
    ??? np = NumberPolyOfStar(poly);
    {
      tri = NEW(REF ARRAY OF Triv, np);
      for(i = 0; i < poly.nel; i++)
      {
        c = 0;
        for (j = 0; j < 4; j++)
        {
          if (( poly[i][j]!=OrgV(a)
          )){
              tri[i,c] = poly[i][j];
              c++
          }
        }
        assert(c == 3);
      }

      for(j = 0; j < poly.nel; j++)
      {
        for (l = 0; l <= (2); l++)
        {
          if ((! Present(tri[j,l])
          )){
              vstack.el[nstack] = tri[j,l];
              nstack++;
          }
        }
      }       
      
      Node_vec_trim(&vstack,nstack);
    }
  } /* END NeighborNode */

Node_vec_t *Neighbors(Place_t @p, Topology_t *top) ==
  
  VAR
    Node_vec_t vstack = Node_vec_new(INIT_STACK_SIZE);
    uint nstack = 0;
    uint v;

  { v = OrgV(a)->num;
    for (i = 0; i < top->NE; i++)
    {
      ??? ei = top->edge[i];
        ??? ei0 = ei.node[0]->num;
      ??? ei1 = ei.node[1]->num;
      {
        if ((v == ei0) || ((v == ei1)
        )){
            if (v == ei0
            ){
                vstack.el[nstack] = top->node[ei1];
                nstack++;
              }
              else
              {
                vstack.el[nstack] = top->node[ei0];
                nstack++;
            }
        }
      }
    }          
      
    Node_vec_trim(&vstack,nstack);
  } /* END Neighbors */
  
uint NumberNeighborNode(Node_vec_t *neigh)
  { ??? n = ((neigh^).nel);
    { return n; };
  } /* END NumberNeighborNode */

Quadv_vec_t *StarOfNode(Place_t @p, Topology_t *top)==

  PROCEDURE NodePoly(*top : Topology):  REF ARRAY OF Quadp ==
     /* Return four places such as the origins corresponding to
        nodes extremes of each tetrahedron of Triangulation.
     */
    REF ARRAY OF Quadp  poly = NEW(REF ARRAY OF Quadp, top->NP);
    {
      for (i = 0; i < top->NP; i++)
      {
        ??? da = top->region[i];
        Place_t a0 = Tors(da);
        Place_t db = Clock(PrevE(da));
        Place_t b0 = Tors(db);
        {
          /*assert(DegreeOfNode(a0) == 4);  */       /* consuming time */
          Place_t a1 = NextE(a0);
          Place_t a2 = NextE(a1);
          Place_t a3 = PrevE(b0);
          {
            if ((NextE(a2)!=a0))
            {
              fprintf(stderr, "\nTriangulation: This topology isn't" \
                         " a triangulation\n");
              assert(NextE(a2) == a0);
            }  
            assert(Pneg(a0)->num == i);
            if ((Pneg(b0)!=NULL))
            {
              assert(Pneg(a0) == Pneg(b0));
            }
            assert(a0!=b0);
            poly[i] = Quadp{a0,a1,a2,a3};     
          }
        }
      }
      return poly;
    } NodePoly;

      poly1 : REF Quadv_vec_t = Quadv_vec_new(top->NP);
      n : uint = 0;
  { ??? poly = NodePoly(top), v == OrgV(a);
    {
      for(i = 0; i < poly.nel; i++)
      {
         if (((OrgV(poly[i][0]) == v)) || ((OrgV(poly[i][1]) == v) OR
            (OrgV(poly[i][2]) == v)) || ((OrgV(poly[i][3]) == v)
       )){
            for (j = 0; j < 4; j++)
            {
               if ((OrgV(poly[i][j]) == v
               )){
                   Node_t uv = OrgV(poly[i][j]);
                     Node_t dv = OrgV(NextE(poly[i][j]));
                   Node_t pv = OrgV(PrevE(poly[i][j]));
                   Node_t rv = OrgV(PrevE(PrevF(poly[i][j])));
                   {
                     poly1[n] = Quadv{uv,dv,pv,rv};
                     n++;
                   }
               }
             }
         }  
      }
      Quadv_vec_trim(&poly1,n);
    }
  } /* END StarOfNode */    

uint NumberPolyOfStar(Quadv_vec_t *qv)
  { ??? n = (qv->nel);
    { return n; };
  } /* END NumberPolyOfStar */

r4_vec_t *ComputeAllNodeNormals(
     *Topology_t *top; 
     Coords_t *c;
    )==
  { ??? rvn = r4_vec_new(top->NV);
    ??? vn = rvn^;
    {
      for (i = 0; i < top->NV; i++)
      {
        vn[i] = NodeNormal(top->out[i], c, top)
      }
      return rvn
    }
  } /* END ComputeAllNodeNormals */ 

r4_vec_t *ComputeAllEdgeNormals(Topology_t *top, 
    Coords_t *c;
   )==
  { ??? rvn = r4_vec_new(top->NE);
    ??? vn = rvn^;
    {
      for (i = 0; i < top->NE; i++)
      {
        vn[i] = EdgeNormal(top->edge[i].pa, c);
      }
      return rvn
    }
  } /* END ComputeAllEdgeNormals */

r4_vec_t *ComputeAllWallNormals(Topology_t *top, 
    Coords_t *c;
   )==
  { ??? rvn = r4_vec_new(top->NF);
    ??? vn = rvn^;
    {
      for (i = 0; i < top->NF; i++)
      {
        vn[i] = WallNormal(top->wall[i].pa, c);
      }
      return rvn
    }
  } /* END ComputeAllWallNormals */

r4_vec_t *ComputeAllCellNormals(Topology_t *top, 
    Coords_t *c;
   )==
  { ??? rvn = r4_vec_new(top->NP);
    ??? vn = rvn^;
    {
      for (i = 0; i < top->NP; i++)
      {
        vn[i] = PolyNormal(Tors(top->region[i]), c);
      }
      return rvn
    }
  } /* END ComputeAllCellNormals */

CellTopology MakeCellTopology(Place_t @p)
  uint *ne, nv; b: Place_t;
      ptop: CellTopology;
  { ??? star = Octf.RenumberEdgesOutOfNodes((Place_vec_t){a})^;
    {
      /* Gather walls: */
      ptop.NF = (star.nel);
      ptop.fRef = Place_vec_new(ptop.NF);
      ne = 0;
      for (i = 0; i < ptop.NF; i++)
      {
        ??? side = Octf.Dual(star[i]);
    {
          ptop.fRef[i] = side;
          b = side;
          do 
          {
            ne++;
            Edge_t e = PEdge(b);
    {
              e.xmark = FALSE; OrgV(b).xmark = FALSE;
            }
            b = NextE(b)
          } while (b != side);
        }
      }

      /* Gather @{edge->?}s: */
      assert(ne % 2 == 0);
      ptop.NE = ne DIV 2;
      ptop.eRef = Place_vec_new(ptop.NE);
      ne = 0;
      for (i = 0; i < ptop.NF; i++)
      {
        ??? side = Octf.Dual(star[i]);
    {
          b = side;
          do 
          {
            Edge_t e = PEdge(b);
    {
              if (! e.xmark)
              {
                ptop.eRef[ne] = b; ne++; e.xmark = TRUE;
              }
            }
            b = NextE(b)
          } while (!( b == side));
        }
      }

      /* Gather nodes: */
      ptop.NV = 2 + ptop.NE - ptop.NF;
      ptop.vRef = Place_vec_new(ptop.NV);
      nv = 0;
      for (i = 0; i < ptop.NE; i++)
      {
        ??? eu = ptop.eRef[i];
        Node_t u = OrgV(eu);
        Place_t ev = Clock(PrevF(eu));
        Node_t v = OrgV(ev);
        { 
          if (! u.xmark){ ptop.vRef[nv] = eu; nv++; u.xmark = TRUE; }
          if (! v.xmark){ ptop.vRef[nv] = ev; nv++; v.xmark = TRUE; }
          PEdge(eu).xmark = FALSE;
        };
      }
      for (i = 0; i < ptop.NV; i++){ Node_t u = OrgV(ptop.vRef[i]);
    { u.xmark = FALSE; } }
    }
    return ptop;
  } /* END MakeCellTopology */

  /* =============== INPUT/OUTPUT =========================== */

CONST
  BoolChars == Mis.BoolChars;
  AlphaChars == Mis.AlphaChars;

TopCom_t ReadTopology(FILE *rd)
  VAR
    Topology_t *top; 
    char *cmt;
    CHAR n;

  { /* Topology */
    ReadHeader(rd,"topology","99-08-25");
    cmt = Mis.ReadCommentsJS(rd, '|');

    /* Element counts: */
    Lex.Skip(rd, cs = AlphaChars);
    top->NV = fget_int(rd);
    Lex.Skip(rd, cs = AlphaChars);
    top->NE = fget_int(rd);
    Lex.Skip(rd, cs = AlphaChars);
    top->NF = fget_int(rd);
    Lex.Skip(rd, cs = AlphaChars);
    top->NP = fget_int(rd);
    Lex.Skip(rd, cs = AlphaChars);
    top->NFE = fget_int(rd);
    Lex.Skip(rd, cs = AlphaChars);
    top->der = fget_int(rd);  
    Lex.Skip(rd,  cs = AlphaChars);
    top->bdr = fget_int(rd);
    Lex.Skip(rd);

    ??? map = NEW(REF ARRAY OF Octf.wedge, top->NFE)^;
    {

      /* Create node records: */
      top->node = Node_vec_new(top->NV);
      top->out = Place_vec_new(top->NV);
      for (i = 0; i < top->NV; i++)
      {
        top->node[i] = MakeNode();
        top->node[i]->num = i;
      }
      
      /* Create @{edge->?}s records */
      top->edge = NEW(REF ARRAY OF Octf.edge, top->NE);
      for (i = 0; i < top->NE; i++)
      {
        top->edge[i] = Octf.Make@{Edge->?}();
        top->edge[i]->num = i;
      }

      /* Create wall records: */
      top->wall = Wall_vec_new(top->NF);
      for (i = 0; i < top->NF; i++)
      {
        top->wall[i] = Octf.MakeWall();
        top->wall[i]->num = i;
      }
      
      /* Create cells records: */
      top->cell = NEW(REF ARRAY OF Cell, top->NP);
      top->region = Place_vec_new(top->NP);
      
      for (i = 0; i < top->NP; i++)
      {
        top->cell[i] = MakeCell();
        top->cell[i]->num = i;
      }
      
      /* Create wedge records: */
      top->wedge = NEW(REF ARRAY OF Octf_Place_t, top->NFE);
      for (i = 0; i < top->NFE; i++)
      {
        top->wedge[i] = MakeWedge();
        PWedge(top->wedge[i])->num = i;
        map[i] = PWedge(top->wedge[i]);
      }

      /* Read @{edge->?} data: */
      EVAL Mis.ReadCommentsJS(rd, '|');
      for (j = 0; j < top->NE; j++)
      {
        Lex.Skip(rd);
        ??? ne = fget_int(rd); /* index to @{edge->?} */ 
        ??? e == top->edge[ne];
        ??? pa == Octf.ReadPlace(rd,map);
        {
          assert(ne == j);
          e->pa = pa
        }
      }
      Lex.Skip(rd);

      /* Read wall data: */
      EVAL Mis.ReadCommentsJS(rd, '|');
      for (j = 0; j < top->NF; j++)
      {
        Lex.Skip(rd);
        ??? nf = fget_int(rd);      /* index to wall */ 
        ??? f = top->wall[nf];
        ??? pa = Octf.ReadPlace(rd,map);
        {
          assert(nf == j);
          f->pa = pa;
        }
      }
      Lex.Skip(rd);

      /* Read wedge data: */
      EVAL Mis.ReadCommentsJS(rd, '|');
      for (j = 0; j < top->NFE; j++)
      {
        Lex.Skip(rd);
        ??? nfe = fget_int(rd); /* index to wedge */
        ??? fe = NARROW(PWedge(top->wedge[nfe]), wedge);
        {
          assert(nfe == j);
          assert(top->wedge[nfe].bits == 0);
          fe->num = nfe;
          Octf.ReadWedge(rd, fe, map);
          for (k = 0; k < 4; k++)
          {
            Lex.Skip(rd);
            n = Rd.GetChar(rd);
            if (n!='-')
            {
              Rd.UnGetChar(rd);
              ??? m = fget_int(rd);
              ??? pa = PlaceInWedge(fe, 2*k);
              ??? vf = fe.org[k];
              {
                if ((Rd.GetChar(rd) == 'v'))
                {
                  vf = top->node[m];
                  top->out[m] = pa;
                }
                else
                { 
                  vf = top->cell[m];
                  top->region[m] = pa;
                }
              }
            }
          }
          Lex.Skip(rd);
          ??? nf = fget_int(rd);
    {
            fe.wall = top->wall[nf];
          }
          assert(Rd.GetChar(rd) == 'f');
          Lex.Skip(rd);
          ??? ne = fget_int(rd);
    {
            fe.edge = top->edge[ne]
          }
          assert(Rd.GetChar(rd) == 'e');
        }
      }
      ReadFooter(rd,"topology");
      fprintf(stderr, "rebuilding node tables...\n");
      RebuildNodeTables(top);
      return (TopCom_t){top, cmt};
    }
  } /* END ReadTopology */

void RebuildNodeTables(Topology *top)
  Place_t *b;
  { for (i = 0; i < top->NE; i++){
      ??? pa = top->edge[i].pa;
    {
        assert(PEdge(pa) == top->edge[i]);
        assert(top->edge[i]->num == i);
        top->edge[i].node[0] = OrgV(pa);
        top->edge[i].node[1] = OrgV(Clock(pa))
      }
    }
    
    for (i = 0; i < top->NF; i++)
    {
      ??? a = top->wall[i].pa;
      ??? deg = Octf.DegreeOfWall(a);
      ??? vr = NEW(REF ARRAY OF Node_t, deg);
      {
        assert(PWall(a) == top->wall[i]);
        assert(top->wall[i]->num == i);
        top->wall[i].node = vr;
        b = a;
        for (j = 0; j < deg; j++)
        {
          vr[j] = Org(b);
          b = NextE(b);
        }
      }
    }

    for (i = 0; i < top->NP; i++)
    {
      ??? a = top->region[i];
      ??? ptop = MakeCellTopology(a);
      ??? vr = NEW(REF ARRAY OF Node_t, ptop.NV);
      {
        assert(Org(a) == top->cell[i]);
        assert(top->cell[i]->num == i);
        top->cell[i].node = vr;
        for (j = 0; j < ptop.NV; j++)
        {
          vr[j] = Org(ptop.vRef[j])
        }
      }
    }
  } /* END RebuildNodeTables */  

PROCEDURE ReadMaterials(
  FILE *rd;
  *Topology_t *top;
  bool_t ro_te = FALSE;
 )
  { /* Materials */
    ReadHeader(rd,"materials","99-08-25");

    /* Read node data Materials: */
    EVAL Mis.ReadCommentsJS(rd, '|');
    for (j = 0; j < top->NV; j++)
    {
      Lex.Skip(rd);
      ??? nv = fget_int(rd); /* index to node */
      ??? v = top->node[nv];   
      { 
        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);
        ??? cc = v.color;
    {
          cc[0] = Lex.Real(rd);
          Lex.Skip(rd);
          cc[1] = Lex.Real(rd);
          Lex.Skip(rd);
          cc[2] = Lex.Real(rd); 
        }
        Lex.Skip(rd);
        ??? tt = v.transp;
    {
          tt[0] = Lex.Real(rd);
          Lex.Skip(rd);
          tt[1] = Lex.Real(rd);
          Lex.Skip(rd);
          tt[2] = Lex.Real(rd);
        }
        Lex.Skip(rd);
        v.radius = Lex.Real(rd);
        Lex.Skip(rd);
        v.label = Rd.GetText(rd,2); 
      }      
    }

    Lex.Skip(rd);

    /* Read @{edge->?} data materials: */
    EVAL Mis.ReadCommentsJS(rd, '|');
    for (j = 0; j < top->NE; j++)
    {
      Lex.Skip(rd);
      ??? ne = fget_int(rd); /* index to @{edge->?} */ 
      ??? e = top->edge[ne];
      {
        assert(ne == j);
        e->num = ne;
        Lex.Skip(rd);
        e->exists = Mis.ReadBool(rd);
        Lex.Skip(rd);
        ??? cc = e.color;
    {
          cc[0] = Lex.Real(rd);
          Lex.Skip(rd);
          cc[1] = Lex.Real(rd);
          Lex.Skip(rd);
          cc[2] = Lex.Real(rd); 
        }
        Lex.Skip(rd);
        ??? tt = e.transp;
    {
          tt[0] = Lex.Real(rd);
          Lex.Skip(rd);
          tt[1] = Lex.Real(rd);
          Lex.Skip(rd);
          tt[2] = Lex.Real(rd);
        }
        Lex.Skip(rd);
        e.radius = Lex.Real(rd);
        Lex.Skip(rd);
        e.degenerate = Mis.ReadBool(rd);
        Lex.Skip(rd);
        ??? n = Rd.GetChar(rd);
        {
          if (n!='-')
          {
                Rd.UnGetChar(rd);
            e->root = fget_int(rd);
          }
          else
          {
            e->root = -1;
          }
        }
      }
    }
    Lex.Skip(rd);

    /* Read wall data materials: */
    EVAL Mis.ReadCommentsJS(rd, '|');
    for (j = 0; j < top->NF; j++)
    {
      Lex.Skip(rd);
      ??? nf = fget_int(rd);      /* index to wall */ 
      ??? f = top->wall[nf];
      {
        assert(nf == j);
        f->num = nf;
        Lex.Skip(rd);
        f->exists = Mis.ReadBool(rd);
        Lex.Skip(rd);
        ??? cc = f.color;
    {
          cc[0] = Lex.Real(rd);
          Lex.Skip(rd);
          cc[1] = Lex.Real(rd);
          Lex.Skip(rd);
          cc[2] = Lex.Real(rd); 
        }
        Lex.Skip(rd);
        ??? tt = f.transp;
    {
          tt[0] = Lex.Real(rd);
          Lex.Skip(rd);
          tt[1] = Lex.Real(rd);
          Lex.Skip(rd);
          tt[2] = Lex.Real(rd);
        }
        Lex.Skip(rd);
        f.degenerate = Mis.ReadBool(rd);
        Lex.Skip(rd);
        ??? n = Rd.GetChar(rd);
        {
          if (n!='-')
          {
                Rd.UnGetChar(rd);
            f->root = fget_int(rd);
          }
          else
          {
            f->root = -1;
          }
        }
      }
    }
    Lex.Skip(rd);

    /* Read cell data materials: */
    if (top->NP!=0)
    {
      EVAL Mis.ReadCommentsJS(rd, '|');
    }
    for (j = 0; j < top->NP; j++)
    {
      Lex.Skip(rd);
      ??? np = fget_int(rd);      /* index to cell */
      ??? p = top->cell[np];
      {
        assert(np == j);
        p->num = np;
        Lex.Skip(rd);
        p->exists = Mis.ReadBool(rd);
        Lex.Skip(rd);
        ??? cc = p.color;
    {
          cc[0] = Lex.Real(rd);
          Lex.Skip(rd);
          cc[1] = Lex.Real(rd);
          Lex.Skip(rd);
          cc[2] = Lex.Real(rd); 
        }
        Lex.Skip(rd);
        ??? tt = p.transp;
    {
          tt[0] = Lex.Real(rd);
          Lex.Skip(rd);
          tt[1] = Lex.Real(rd);
          Lex.Skip(rd);
          tt[2] = Lex.Real(rd);
        }
        Lex.Skip(rd);
        p.degenerate = Mis.ReadBool(rd);
        if (ro_te)
        {
          Lex.Skip(rd);
          ??? n = Rd.GetChar(rd);
          {
            if (n!='-')
            {
              Rd.UnGetChar(rd);
              p->root = fget_int(rd);
            }
            else
            {
              p->root = -1;
            }
          }
        }
      }
    }

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

TopCom_t ReadToTaMa(char *name, bool_t ro_te /* DF FALSE */)
  /* Where ro_te meaning "root tetrahedron". */
  VAR
    TopCom_t tc;
  { ??? ntp = name & ".tp";
    ??? rtp = FileRd.Open(ntp);
    {
      fprintf(stderr, "reading " & ntp & "\n");
      tc = ReadTopology(rtp);
      fprintf(stderr, "OK leio tp\n");
      Rd.Close(rtp);
    }  
    ??? nma = name & ".ma";
    ??? rma = FileRd.Open(nma);
    {
      fprintf(stderr, "reading " & nma & "\n");
      ReadMaterials(rma, tc.top, ro_te);
      Rd.Close(rma);
    }
    fprintf(stderr, "OK leio ma\n");
    return tc;
  } /* END ReadToTaMa */

Coords_t *ReadState(char *name)==
  <* FATAL Rd.Failure, Thread.Alerted,FloatMode.Trap, Lex.Error, OSError.E);
  Coords_t c;
    char *cmt;
    uint nv;
  { ??? rs = FileRd.Open(name & ".st");
    {
      /* Read Headers File Formats*/
      ReadHeader(rs,"state","99-08-25");

      /* Element counts: */
      Lex.Skip(rs, cs = AlphaChars);
      nv = fget_int(rs);
      Lex.Skip(rs);
      cmt = Mis.ReadCommentsJS(rs, '|');
      
      c = NEW(REF Coords_t, nv);
      /* Read node data state: */
      for (j = 0; j < nv; j++)
      {
        Lex.Skip(rs);
        ??? nv = fget_int(rs);
        { 
          ??? cv = c[nv];
    {
            cv[0] = Lex.LongReal(rs);
            Lex.Skip(rs);
            cv[1] = Lex.LongReal(rs);
            Lex.Skip(rs);
            cv[2] = Lex.LongReal(rs);
            Lex.Skip(rs);
            cv[3] = Lex.LongReal(rs);
          }
        } 
      }

      ReadFooter(rs,"state");
      Rd.Close(rs);        
      return c;
    }  
  } /* END ReadState */

void WriteState(char *name, Topology_t *top, Coords_t *c,
    char *cmt /* DF "" */)

  <* FATAL OSError.E);
  { ??? st = FileWr.Open(name & ".st");
    ??? vWidth = NumDigits(top->NV - 1);
    {
      void WriteCoord(x: double)
        {
          Wr.PutText(st, Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Sci, prec = 3),
          7));
        } WriteCoord;

      void WritePoint(*c: r4_t)
        {
          WriteCoord(c[0]);
          Wr.PutText(st, " ");
          WriteCoord(c[1]);
          Wr.PutText(st, " ");
          WriteCoord(c[2]);
          Wr.PutText(st, " ");
          WriteCoord(c[3]);
        } WritePoint;
      
      {
        WriteHeader(st,"state","99-08-25");
        Wr.PutText(st, "nodes ");
        Wr.PutText(st, Fmt.Int(top->NV) & "\n");
        if ((! Text.Empty(cmt)))
        {
          Mis.WriteCommentsJS(st, cmt & "\n", '|')
        }

        Mis.WriteCommentsJS(st, "\nNode data:\n", '|');
        
        for (i = 0; i < top->NV; i++)
        {
          ??? v = top->node[i];
    {
            /* state */
            Wr.PutText(st, Fmt.Pad(Fmt.Int(v->num), vWidth)); 
            Wr.PutText(st, "  ");
            WritePoint(c[v->num]);
            Wr.PutText(st, "\n");
          }       
        }
      }
      WriteFooter(st, "state");
      fclose(st);
    }     
  } /* END WriteState */

void WriteMaterials(char *name, Topology_t *top,
   char *cmt /* DF "" */,
   bool_t ro_te = FALSE;  /* root tetrahedron */
 )
  <* FATAL Wr.Failure, Thread.Alerted, OSError.E);
  INTEGER *pWidth;
  { if (top->NP == 0){ pWidth = 2 }else{ pWidth = NumDigits(top->NP-1); }
    ??? ma = FileWr.Open(name & ".ma");
    ??? vWidth = NumDigits(top->NV - 1);
    ??? fWidth = NumDigits(top->NF - 1);
    ??? eWidth = NumDigits(top->NE - 1);
    { 
      void WriteIntensity(r: REAL)
        {
          Wr.PutText(ma, Fmt.Real(r, Fmt.Style.Fix, prec = 2));
        } WriteIntensity;

      void WriteColor(*c: frgb_t)
        {
          WriteIntensity(c[0]);
          Wr.PutText(ma, " ");                     
          WriteIntensity(c[1]);
          Wr.PutText(ma, " ");                     
          WriteIntensity(c[2]);
        } WriteColor;
        
      void WriteRadius(r: REAL)
        {
          if (r == 0.02
         ){
             Wr.PutText(ma, "0.020");
            }
            else
            {
             Wr.PutText(ma,Fmt.Real(r, prec = 4)); 
          }
        } WriteRadius;

      void WriteLabel(char *label)
        {
          Wr.PutText(ma, label);
        } WriteLabel;

      {
        WriteHeader(ma,"materials","99-08-25");
        if ((! Text.Empty(cmt)))
        {
          WriteCommentsJS(ma, cmt & "\n", '|')
        }
        ??? m = NumDigits(top->NP);
        {
          WriteCommentsJS(ma,"nodes    " & Fmt.Pad(Fmt.Int(top->NV),m), '|');
          WriteCommentsJS(ma,"@{edge->?}s       " & Fmt.Pad(Fmt.Int(top->NE),m), '|');
          WriteCommentsJS(ma,"walls       " & Fmt.Pad(Fmt.Int(top->NF),m), '|');
          WriteCommentsJS(ma,"cells   " & Fmt.Int(top->NP), '|');
        }
        WriteCommentsJS(ma, "\nNode data:\n", '|');

        for (i = 0; i < top->NV; i++)
        {
          ??? v = top->node[i];
    {
            /* materials */
            Wr.PutText(ma, Fmt.Pad(Fmt.Int(v->num), vWidth)); 
            Wr.PutText(ma, "  ");
            Wr.PutText(ma, Fmt.Char(BoolChars[v->exists]));
            Wr.PutText(ma, " ");
            Wr.PutText(ma, Fmt.Char(BoolChars[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");
          }       
        }

        WriteCommentsJS(ma, "\n@{Edge->?} data:\n", '|');
        for (i = 0; i < top->NE; i++)
        {
          ??? e = top->edge[i];
    {
            /* materials */
            Wr.PutText(ma, Fmt.Pad(Fmt.Int(e->num), eWidth)); 
            Wr.PutText(ma, "  ");
            Wr.PutText(ma, Fmt.Char(BoolChars[e->exists]));
            Wr.PutText(ma, "  ");
            WriteColor(e.color);
            Wr.PutText(ma, "  ");
            WriteColor(e.transp);
            Wr.PutText(ma, "  ");
            WriteRadius(e.radius);
            Wr.PutText(ma, "  ");
            Wr.PutText(ma, Fmt.Char(BoolChars[e.degenerate]));
            Wr.PutText(ma, "  ");
            if (e->root == -1)
            {
              Wr.PutText(ma, " - ");
            }
            else
            {
              Wr.PutText(ma, Fmt.Pad(Fmt.Int(e->root), eWidth));
            }
            Wr.PutText(ma, "\n");
          }
        }

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

        for (i = 0; i < top->NF; i++)
        {
          ??? f = top->wall[i];
    { 
            Wr.PutText(ma, Fmt.Pad(Fmt.Int(f->num), fWidth)); 
            Wr.PutText(ma, "  ");
            Wr.PutText(ma, Fmt.Char(BoolChars[f->exists]));
            Wr.PutText(ma, "  ");
            WriteColor(f.color);
            Wr.PutText(ma, "  ");
            WriteColor(f.transp);
            Wr.PutText(ma, "  ");
            Wr.PutText(ma, Fmt.Char(BoolChars[f.degenerate]));     
            Wr.PutText(ma, "  ");
            if (f->root == -1)
            {
              Wr.PutText(ma, " - ");
            }
            else
            {
              Wr.PutText(ma, Fmt.Pad(Fmt.Int(f->root), fWidth));
            }
            Wr.PutText(ma, "\n");
          }
        }

        if (top->NP!=0)
        {
          WriteCommentsJS(ma, "\nCell data:\n", '|');
        }

        for (i = 0; i < top->NP; i++)
        {
          ??? p = top->cell[i];
          {
            Wr.PutText(ma, Fmt.Pad(Fmt.Int(p->num), pWidth)); 
            Wr.PutText(ma, "  ");
            Wr.PutText(ma, Fmt.Char(BoolChars[p->exists]));
            Wr.PutText(ma, "  ");
            WriteColor(p.color);
            Wr.PutText(ma, "  ");
            WriteColor(p.transp);
            Wr.PutText(ma, "  ");
            Wr.PutText(ma, Fmt.Char(BoolChars[p.degenerate]));
            if (ro_te)
            {
              Wr.PutText(ma, "  ");
              if (p->root == -1)
              {
                Wr.PutText(ma, " - ");
              }
              else
              {
                Wr.PutText(ma, Fmt.Pad(Fmt.Int(p->root), pWidth));
              }
            }
            Wr.PutText(ma, "\n");
          }
        }
      }
      WriteFooter(ma, "materials");
      fclose(ma);
    }     
  } /* END WriteMaterials */


void FindDegeneracies(Topology *top)

  VAR
    bool_t allTriangles, allTetrahedra;
    vei, vej: REF ARRAY OF INTEGER;
    vfi, vfj: REF ARRAY OF INTEGER;
    vpi, vpj: REF ARRAY OF INTEGER;
  { /* Clean the marks for the attribute degenerate */
    
    for (i = 0; i < top->NE; i++)
    {
      ??? ei = top->edge[i];
    {
        ei.degenerate = FALSE
      }
    }

    /* Now, find @{edge->?} degeneracies */

    vei = NEW(REF ARRAY OF INTEGER, 2);
    vej = NEW(REF ARRAY OF INTEGER, 2);
    for (i = 0; i < top->NE; i++)
    {
      for (j = i+1; j < top->NE; j++)
      {
        ??? ei = NARROW(top->edge[i], @{Edge->?});
        ??? ej = NARROW(top->edge[j], @{Edge->?});
        {
          for (k = 0; k < 2; k++)
          {
            vei[k] = ei.node[k]->num;
            vej[k] = ej.node[k]->num
          }
          Mis.InsertionSort(1,vei);
          Mis.InsertionSort(1,vej);
          if ((vei[0] == vej[0]) && (vei[1] == vej[1]))
          {
            ei.degenerate = TRUE;
            ej.degenerate = TRUE;
          }
        }  
      }
    }

    /* Clean the marks for the attribute degenerate */
    allTriangles = TRUE;   
    for (i = 0; i < top->NF; i++)
    {
      ??? fi = top->wall[i];
        ??? di = ((fi.node^).nel);
      {
        fi.degenerate = (di!=3);
        allTriangles = allTriangles) && ((di == 3)
      }
    }
    if (! allTriangles){ return; }
    
    /* Now, find degenerate @places of triangles */
    vfi = NEW(REF ARRAY OF INTEGER, 3);
    vfj = NEW(REF ARRAY OF INTEGER, 3);
    for (i = 0; i < top->NF; i++)
    {
      for (j = i+1; j < top->NF; j++)
      {
        ??? fi = top->wall[i];
        ??? fj = top->wall[j];
        {
          for (k = 0; k <= (2); k++)
          {
            vfi[k] = fi.node^[k]->num;
            vfj[k] = fj.node^[k]->num
          }
          Mis.InsertionSort(2,vfi);
          Mis.InsertionSort(2,vfj);
          if ((vfi[0] == vfj[0]) && (vfi[1] == vfj[1]) && (vfi[2] == vfj[2]))
          {
            fi.degenerate = TRUE;
            fj.degenerate = TRUE
          }
        }
      }
    }

    /* Clean the marks for the attribute degenerate */
    allTetrahedra = TRUE;   
    for (i = 0; i < top->NP; i++)
    {
      ??? pi = top->cell[i];
        ??? di = ((pi.node^).nel);
      {
        pi.degenerate = (di!=4);
        allTetrahedra = allTetrahedra) && ((di == 4)
      }
    }
    if (! allTetrahedra){ return; }

    /* Now, find cell degeneracies */
      
    vpi = NEW(REF ARRAY OF INTEGER, 4);
    vpj = NEW(REF ARRAY OF INTEGER, 4);
    for (i = 0; i < top->NP; i++)
    {
      for (j = i+1; j < top->NP; j++)
      {
        ??? pi = top->cell[i];
        ??? pj = top->cell[j];
        {   
          for (k = 0; k < 4; k++)
          {
            vpi[k] = pi.node^[k]->num;
            vpj[k] = pj.node^[k]->num
          }
          Mis.InsertionSort(3,vpi);
          Mis.InsertionSort(3,vpj);
          if ((vpi[0]==vpj[0]) && (vpi[1]==vpj[1]) && (vpi[2]==vpj[2]) && (vpi[3]==vpj[3]))
          {
            pi.degenerate = TRUE;
            pj.degenerate = TRUE;
          }  
        }
      }
    }  
  } /* END FindDegeneracies */

void WriteTopology(char *name, Topology_t *top, char *cmt /* DF " " */)
  <* FATAL Wr.Failure, Thread.Alerted, OSError.E);
  INTEGER *pWidth;
  { if (top->NP == 0){  pWidth= 2 }else{ pWidth = NumDigits(top->NP-1); }
    ??? tp = FileWr.Open(name & ".tp");
    ??? vWidth = NumDigits(top->NV - 1);
    ??? eWidth = NumDigits(top->NE - 1);
    ??? fWidth = NumDigits(top->NF - 1);
    ??? feWidth = NumDigits(top->NFE -1);
    {
      WriteHeader(tp,"topology","99-08-25");
      if ((! Text.Empty(cmt)))
      {
        WriteCommentsJS(tp, cmt & "\n", '|')
      }

      ??? m = NumDigits(top->NFE);
      {
        Wr.PutText(tp, "nodes   ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NV),m) & "\n");  
        Wr.PutText(tp, "@{edge->?}s      ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NE),m) & "\n");
        Wr.PutText(tp, "walls      ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NF),m) & "\n"); 
        Wr.PutText(tp, "cells  ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NP),m) & "\n");
        Wr.PutText(tp, "wedges ");
        Wr.PutText(tp, Fmt.Int(top->NFE) & "\n");
        Wr.PutText(tp, "der        ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->der),m) & "\n");
        Wr.PutText(tp, "bdr        ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->bdr),m) & "\n");
      }

      WriteCommentsJS(tp, "\n@{Edge->?} data:\n", '|');
      for (i = 0; i < top->NE; i++)
      {
        ??? e = top->edge[i];
    { 
          Wr.PutText(tp, Fmt.Pad(Fmt.Int(e->num), eWidth)); 
          Wr.PutText(tp, "  ");
          Octf.PrintPlace(tp, e.pa, eWidth+1, TRUE);  
        }
      }

      WriteCommentsJS(tp, "\nWall data:\n", '|');
      for (i = 0; i < top->NF; i++)
      {
        ??? f = top->wall[i];
    { 
          Wr.PutText(tp, Fmt.Pad(Fmt.Int(f->num), fWidth)); 
          Wr.PutText(tp, "  ");
          Octf.PrintPlace(tp, f.pa, fWidth+1, TRUE);  
        }
      }

      WriteCommentsJS(tp, "\nwedge data:\n", '|');
      for (i = 0; i < top->NFE; i++)
      {
        ??? fe = NARROW(PWedge(top->wedge[i]), wedge);
    { 
          Wr.PutText(tp, Fmt.Pad(Fmt.Int(fe->num), feWidth)); 
          Wr.PutText(tp, "  ");
          Octf.PrintWedge(tp, fe, feWidth);
          Wr.PutText(tp, "  ");
          for (j = 0; j < 4; j++)
          {
            ??? n = fe.org[j];
    { 
              TYPECASE n OF 
              | NULL ==> for (i = 0; i <= (pWidth-2); i++){ Wr.PutText(tp," "); }
                        Wr.PutText(tp, " - ");
              | Node(v) ==> Wr.PutText(tp, Fmt.Pad(Fmt.Int(v->num), vWidth)
                & "v ");
              | Cell(p) ==> Wr.PutText(tp, Fmt.Pad(Fmt.Int(p->num),
                pWidth) & "p "); 
              }
              else
              { /* nothing */
              } 
            }
          }
          Wr.PutText(tp, "  ");
          ??? f = fe.wall;
    {
            Wr.PutText(tp, Fmt.Pad(Fmt.Int(f->num), fWidth));
          }
          Wr.PutText(tp, "f ");
          ??? e = fe.edge;
    {
            Wr.PutText(tp, Fmt.Pad(Fmt.Int(e->num), eWidth));
          }
          Wr.PutText(tp, "e\n");
        }
      }
      WriteFooter(tp, "topology");
      fclose(tp);
    }     
  } /* END WriteTopology */

void WriteDualTopology(char *name, Topology_t *top, char *cmt /* DF " " */)
  <* FATAL Wr.Failure, Thread.Alerted, OSError.E);

  void PrintDualPlace_t(FILE *wr, Place_t @p; feWidth: uint)
    {
      Wr.PutText(wr,Fmt.Pad(Fmt.Int(PWedge(a)->num), feWidth) & ":" \
        Fmt.Int((Octf.SrotBits(a)+ 3) % 4) & ":" & Fmt.Int(Octf.SpinBit(a)));
    } PrintDualPlace_t;  

  void PrintDualwedge(FILE *wr, Wedge_t w, feWidth: uint)
    Place_t *b;
    {
      b = Srot(PlaceInWedge(w, 0));
      for (i = 0; i < 4; i++)
      {
        PrintDualPlace_t(wr, NextF(b), feWidth); Wr.PutText(wr, " "); 
        b = Srot(b)
      }
    } PrintDualwedge;  

  { ??? tp = FileWr.Open(name & ".tp");
    ??? vWidth = NumDigits(MAX(1, top->NV - 1));
    ??? eWidth = NumDigits(MAX(1,top->NE - 1));
    ??? fWidth = NumDigits(MAX(1,top->NF - 1));
    ??? pWidth = NumDigits(MAX(1,top->NP - 1));
    ??? feWidth = NumDigits(MAX(1,top->NFE -1));
    {
      WriteHeader(tp,"topology","99-08-25");
      if ((! Text.Empty(cmt)))
      {
        WriteCommentsJS(tp, cmt & "\n", '|')
      }

      ??? m = NumDigits(top->NFE);
      {
        Wr.PutText(tp, "nodes   ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NP),m) & "\n");  
        Wr.PutText(tp, "@{edge->?}s      ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NF),m) & "\n");
        Wr.PutText(tp, "walls      ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NE),m) & "\n"); 
        Wr.PutText(tp, "cells  ");
        Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NV),m) & "\n");
        Wr.PutText(tp, "wedges ");
        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(1),m) & "\n");
      }

      WriteCommentsJS(tp, "\n@{Edge->?} data:\n", '|');
      for (i = 0; i < top->NF; i++)
      {
        ??? f = top->wall[i];
    { 
          Wr.PutText(tp, Fmt.Pad(Fmt.Int(f->num), fWidth)); 
          Wr.PutText(tp, "  ");
          Octf.PrintPlace(tp, f.pa, feWidth); /* [sic] */
          Wr.PutText(tp, "\n") 
        }
      }

      WriteCommentsJS(tp, "\nWall data:\n", '|');
      for (i = 0; i < top->NE; i++)
      {
        ??? e = top->edge[i];
    { 
          Wr.PutText(tp, Fmt.Pad(Fmt.Int(e->num), eWidth)); 
          Wr.PutText(tp, "  ");
          Octf.PrintPlace(tp, e.pa, feWidth+1); /* [sic] */
          Wr.PutText(tp, "\n")  
        }
      }

      WriteCommentsJS(tp, "\nwedge data:\n", '|');
      for (i = 0; i < top->NFE; i++)
      {
        ??? fe = PWedge(top->wedge[i]);
    { 
          Wr.PutText(tp, Fmt.Pad(Fmt.Int(fe->num), feWidth)); 
          Wr.PutText(tp, "  ");
          PrintDualwedge(tp, fe, feWidth);
          Wr.PutText(tp, "  ");
          for (j = 0; j < 4; j++)
          {
            ??? n = fe.org[(j+1) % 4];
    { 
              TYPECASE n OF 
              | NULL ==> for (i = 0; i <= (pWidth-2); i++){ Wr.PutText(tp," "); }
                        Wr.PutText(tp, " - ");
              | Node(v) ==> Wr.PutText(tp, Fmt.Pad(Fmt.Int(v->num), vWidth)
                & "p ");
              | Cell(p) ==> Wr.PutText(tp, Fmt.Pad(Fmt.Int(p->num),
                pWidth) & "v "); 
              }
              else
              { /* nothing */
              } 
            }
          }
          Wr.PutText(tp, "  ");
          ??? e = fe.edge;
    {
            Wr.PutText(tp, Fmt.Pad(Fmt.Int(e->num), eWidth));
          }
          Wr.PutText(tp, "f ");
          ??? f = fe.wall;
    {
            Wr.PutText(tp, Fmt.Pad(Fmt.Int(f->num), fWidth));
          }
          Wr.PutText(tp, "e");
          Wr.PutText(tp, "\n");
        }
      }
      WriteFooter(tp, "topology");
      fclose(tp);
    }     
  } /* END WriteDualTopology */

void WriteTable(char *name, Topology_t *top,
    char *cmt /* DF "" */,
    <*UNUSED); bool_t debug = FALSE;
 ) ==
  <* FATAL Wr.Failure, Thread.Alerted, OSError.E);

  { RebuildNodeTables(top);
    ??? tab = FileWr.Open(name & ".tb");
    ??? vWidth = NumDigits(MAX(1, top->NV - 1));
    ??? eWidth = NumDigits(MAX(1, top->NE - 1));
    ??? fWidth = NumDigits(MAX(1, top->NF - 1));
    ??? pWidth = NumDigits(MAX(1, top->NP - 1));
    { 
      WriteHeader(tab,"table","99-08-25");
      if ((! Text.Empty(cmt)))
      {
        Mis.WriteCommentsJS(tab, cmt & "\n", '|')
      }
    
      /* Now, writes @{edge->?} nodes */
      Mis.WriteCommentsJS(tab,"\n@{Edge->?} data:\n", '|');
      for (i = 0; i < top->NE; i++)
      {
        ??? e = top->edge[i];
        {
          Wr.PutText(tab, Fmt.Pad(Fmt.Int(e->num), eWidth)); 
          Wr.PutText(tab, "  ");
          Wr.PutText(tab, Fmt.Pad(Fmt.Int(e.node[0]->num),eWidth) & "v ");
          Wr.PutText(tab, Fmt.Pad(Fmt.Int(e.node[1]->num),eWidth) & "v\n");
        }
      }

      /* Now, writes wall nodes */
      Mis.WriteCommentsJS(tab, "\nWall data:\n", '|');
      for (i = 0; i < top->NF; i++)
      {
        ??? f = top->wall[i];
          ??? fv = f.node^;
        { 
          Wr.PutText(tab, Fmt.Pad(Fmt.Int(f->num), fWidth)); 
          Wr.PutText(tab, "  ");  
          for (j = 0 TO (fv.nel-1))
          {
            Wr.PutText(tab,Fmt.Pad(Fmt.Int(fv[j]->num), vWidth) & "v ");
          }
          Wr.PutText(tab, "\n");
        }
      }

      /* Now, writes cell nodes */
      if (top->NP!=0)
      {
        Mis.WriteCommentsJS(tab, "\nCell data:\n", '|');
      }
      for (i = 0; i < top->NP; i++)
      {
        ??? p = top->cell[i];
        ??? pv = p.node^;
        { 
          Wr.PutText(tab, Fmt.Pad(Fmt.Int(p->num), pWidth)); 
          Wr.PutText(tab, "  ");
          for (j = 0 TO (pv.nel-1))
          {
            Wr.PutText(tab,Fmt.Pad(Fmt.Int(pv[j]->num), vWidth) & "v ");
          }
          Wr.PutText(tab, "\n");
        }
      }
    WriteFooter(tab, "table");
    fclose(tab);
    }
  } /* END WriteTable */

void WriteStDe(
    FILE *wr,
    Coords_t *c,
    Coords_t *Dc;
    uint prec /* DF 4 */,
    char *cmt /* DF "" */)
    
  void WriteCoord(x: double)
      
    {
      Wr.PutText(wr, Fmt.LongReal(x, Fmt.Style.Sci, prec = prec))
    } WriteCoord;

  void WritePoint(*p: r4_t)
      
    {
      WriteCoord(p[1]);
      Wr.PutText(wr, " ");
      WriteCoord(p[0]);
      Wr.PutText(wr, " ");
      WriteCoord(p[2]);
      Wr.PutText(wr, " ");
      WriteCoord(p[3]);
      Wr.PutText(wr, " ");
    } WritePoint;

    
  { int NV = (c.nel);
    ??? d = NumDigits(NV-1);
    {
      WriteHeader(wr,"state-derivatives","99-08-25");
      Mis.WriteCommentsJS(wr, "\n" & cmt & "\n",'|');
      Wr.PutText(wr, "nodes == " & Fmt.Int(NV) & "\n");

      for (i = 0; i < NV; i++)
      {
        Wr.PutText(wr, Fmt.Pad(Fmt.Int(i), d) & ": "); 
        WritePoint(c[i]); 
        Wr.PutText(wr, "  ");

        /* if the derivatives of nodes are zero them writes "0 0 0 0"
           else writes velocites with the format for write points
           defined here.
        */
        if ((Dc[i][0] == 0.0) && (Dc[i][1] == 0.0) && (Dc[i][2] == 0.0 AND
           Dc[i][3] == 0.0
      )){
          Wr.PutText(wr, "0 0 0 0")
        }
        else
        {
          WritePoint(Dc[i])
        } 
        Wr.PutText(wr, "\n");
      }

      WriteFooter(wr, "state-derivatives");
      Wr.PutText(wr, "\n");
      fflush(wr);
    }
  } /* END WriteStDe */

{
} Triangulation.

  /* Copyright © 2000     Universidade Estadual de Campinas (UNICAMP)     */