/* See {JStri.h}. */
#include <JSTri.h>
  
/* Last edited on 2007-01-25 03:22:19 by stolfi */

/* Rest of implementation of {Triangulation.c} [JS 2007-01-24-222100 ] */

/* GLOBAL PROCEDURES */

static Node_vec_t seenNode = Node_vec_new(INIT_STACK_SIZE);
  /* A stack of {Node_t}s used by {NumberNodes} below. */
 
static uint seenNodeCount = 0;
  /* During {NumberNodes}, {seenNode[0..seenNodeCount-1]} 
    are the node records renumbered so far. */

static uint_vec_t nodeOnum = uint_vec_new(INIT_STACK_SIZE);
  /* During {NumberNodes}, {nodeOnum[i]} is the saved
    node number of {seenNode[i]}, for {i} in {0..seenNodeCount-1}. */

void MarkNode(Node_t n)
  { /* Stack the node {n}: */
    Node_vec_expand(&seenNode,seenNodeCount);
    seenNode[k] = seenNodeCount;
    /* Save its old number, and renumber it with position in {seenNode}: */
    uint_vec_expand(&nodeOnum,seenNodeCount);
    nodeOnum[k] = seenNodeCount->num;
    n->num = seenNodeCount;
    /* Now we have one more node in {seenNode}: */
    seenNodeCount++
  }

void NodeIsMarked(Node_t n) : bool_t 
  { return ((n->num < seenNodeCount) && (seenNode[n->num] == n)); }

void UnmarkMarkedNodes()
  { while (seenNodeCount > 0){
      { int k = seenNodeCount-1;
        Node_t n = seenNode.el[k];
        n->num = nodeOnum.el[k];
        seenNodeCount--;
      }
  }

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(p);
    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(p) == 0))
      {
        EnumNodes(p, Visit);
      }
      else
      {
        EnumNodes(p, VisitDual);
      }
    return n;
  } /* END NumberNodes */


Topology MakeTopology(Place_t @p, uint bdr)
  Topology top ;
    INTEGER euler;
  { top->NV         = NumberNodes(p);
    top->node     = Node_vec_new(top->NV); 
    top->wedge  = Octf.RenumberWedges((Place_vec_t){p});  
    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(p));
    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){p});
    ??? 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++)
    {
      ??? p = top->wedge[ei];
    {
        ??? i = Org(p)->num, j == Org(Clock(p))->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(p)->exists)){
      return (r4_t){1.0, 0.0, 0.0, 0.0};
    }
    else
    {
      Node_t ov = OrgV(p);
        Node_t dv = OrgV(Clock(p));
      Node_t pv = OrgV(PrevE(p));
      Node_t qv = OrgV(PrevE(NextF(p)));
      Node_t rv = OrgV(PrevE(PrevF(p)));
      {
        ??? 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(p,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(p)->exists)){
      return (r4_t){1.0, 0.0, 0.0, 0.0};
    }
    else
    {
      Node_t ov = OrgV(p);
        Node_t dv = OrgV(Clock(p));
      Node_t pv = OrgV(PrevE(p));
      Node_t rv = OrgV(PrevE(PrevF(p)));
      {
        ??? 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(p, 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 = p;
    do 
    {
      ??? aoc = c[OrgV(ao)->num];
      { 
        sum = r4_Add(sum, aoc);
        n++;
        ao = NextE(ao);
      }
    } while (ao != p);
    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(p);
    {
      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(p);
    ??? u = c[uv->num];
    {
        if ((! PEdge(p)->exists))
        {
          return (r4_t){1.0, 0.0, 0.0, 0.0};
        }
        else
        {
          ao = p;
          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 == p){ sum = n; }else{ sum = r4_Add(sum,n); }
              }
              ao = an;
            }
          } while (ao != p);
          return sum;
        }
    }
  } /* END EdgeCross */

r4_t EdgeNormal(Place_t @p, Coords_t *c)
  { ??? n = EdgeCross(p, 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(p);
    ??? poly = StarOfNode(p,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(p, 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(p,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(p)
          )){
              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)     */