/* See {Triangulation.h}. */
#include <Triangulation.h>
  
/* Last edited on 2023-02-12 07:47:45 by stolfi */

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

/* GLOBAL PROCEDURES */

#define INIT_STACK_SIZE 1024

Place_vec_t RenumberWedges(Place_vec_t *pv)
  { Place_vec_t pstack = Place_vec_new(INIT_STACK_SIZE);
    uint nstack = 0;
    uint seen = 0; /* Number of wedges whose childeren were looked at */

    /* A wedge {n} is "marked" if PWedge(pstack.el[n.num]) == n. i.e whether stacked */

    auto void VisitAndMark(Place_t t);
      /* If {t} is unmarked: visit, mark, and stack it. */

    void VisitAndMark(Place_t t)
      { Wedge_t w = PWedge(t);
        if ((w->num < nstack) && (PWedge(pstack.el[w->num]) == w))
          { /* wedge is already marked, do nothing. */ }
        else
          { Place_vec_expand(&pstack,nstack);
            w->num = nstack;
            pstack.el[nstack] = t;
            nstack++;
          }
      }

    /* Enumerate wedges: */
    int i;
    for (i = 0; i < pv->nel; i++) { VisitAndMark(pv->el[i]); } 
    while (seen < nstack)
      { Place_t s = pstack.el[seen]; 
        VisitAndMark(PrevF(s));
        VisitAndMark(PrevE(s));
        VisitAndMark(NextE(s));
        VisitAndMark(NextF(PrevE(s)));
        seen++;
      }
    /* Trim and return the {pstack}: */
    Place_vec_trim(&pstack,nstack);
    return pstack;
  }

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.el[seenNodeCount] = n;
    /* Save its old number, and renumber it with position in {seenNode}: */
    uint_vec_expand(&nodeOnum, seenNodeCount); 
    nodeOnum.el[seenNodeCount] = n->num;
    n->num = seenNodeCount;
    /* Now we have one more node in {seenNode}: */
    seenNodeCount++;
  }

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

void UnmarkMarkedNodes(void)
  { 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: VisitFunc_t)
    {
      Place_vec_t pstack = Place_vec_new(INIT_STACK_SIZE);
      uint top;    /* top for {pstack} */

      void Stack(Place_t c)
        { Place_t cn = c;
          do 
            { Place_vec_expand(&pstack, top);
              pstack.el[top] = cn;
              top++;
              Place_t 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 = OrgV(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 = OrgV(c);
    {
          cn = c;
          do 
          {
            ??? vn = OrgV(cn);
    {
              assert(vn == v);
              vn->num = n;
              dn = Clock(PrevE(cn));
              do 
              {
                ??? wn = OrgV(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 = OrgV(c);
         if (p != NULL)
         {
           cn = c;
           do 
           {
             pn = OrgV(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 */


TetraCorners GetTetraCorners(Place_t @p)
  TetraCorners *c;
  { ??? a0 = a;
    Place_t a1 = NextE(a0);
    Place_t a2 = NextE(a1);
    Place_t b0 = Clock(PrevF(a0));
    Place_t c0 = Clock(PrevF(a1));
    Place_t d0 = Clock(PrevF(a2));
    {
      /* the 12 elements on the array must have the same negative
         cell.
      */
      c[0 ] = a0;
      c[1 ] = NextE(a0);
      c[2 ] = PrevE(a0);
      c[3 ] = b0;
      c[4 ] = NextE(b0);
      c[5 ] = PrevE(b0);
      c[6 ] = c0;
      c[7 ] = NextE(c0);
      c[8 ] = PrevE(c0);
      c[9 ] = d0;
      c[10] = NextE(d0);
      c[11] = PrevE(d0);
      assert(c[0] == a);
    }
    return c;
  } /* END GetTetraCorners */
  
Place_vec_t *CollectTetrahedra(ElemTableRec_t *top)==
      seen,processed,conflicts: uint = 0;
  { ??? rt = Place_vec_new(top.cell.nel);
    ??? t = rt^;
       ??? q = Place_vec_new(top.cell.nel);
     {
      for (i = 0 TO (t.nel-1))
      {
        t[i] = PlaceInWedge(NULL, 0);
      }
      for (k = 0 TO (t.nel-1))
      {
        if ((PWedge(t[k]) == NULL))
        {
          t[k] = Tors(Srot(top.cell[k]));
          q[seen] = t[k];
          seen++;       
          while (processed < seen)
          {
            ??? a = q[processed];
            {
              void Visit(f: Place_t)
                {
                  if ((PnegP(f) == NULL)){ return; }
                  assert(Consistent(a,f));
                  ??? n = PnegP(f)->num;
                  {
                    if ((PWedge(t[n]) == NULL))
                    {
                      t[n] = f;
                      q[seen] = f;
                      seen++;
                    }
                    else
                    {
                      assert(PnegP(t[n])->num == n);
                      if ((! SameOrientation(t[n],f))){ conflicts++;}
                    }    
                  }
                } Visit;
                
              bool_t Consistent(p,q:Place_t)
                /* Checks whether PnegP(q) and the PnegP(q) are adjacentes tetrahedra
                   with consistent orientations.
                */
                {
                  ??? cp = GetTetraCorners(p);
                  ??? cq = GetTetraCorners(q);
                  {
                    for (i = 0; i < 12; i++)
                    {
                      for (j = 0; j < 12; j++)
                      {
                        if ((cp[i] == Clock(cq[j]))){ return TRUE; }
                      }
                    }
                    return FALSE;
                  }
                } Consistent;
          
              bool_t SameOrientation(p,q:Place_t)
                /* Checks whether PnegP(p) == PnegP(q) and both have the same
                   orientation.
                */
                {
                  ??? cp = GetTetraCorners(p);
                  {
                    for (i = 0; i < 12; i++)
                    {
                      if (cp[i] == q){ return TRUE; }
                    }
                    return FALSE;
                  }
                } SameOrientation;
                
              {
                Place_t b = NextF(a);
                Place_t c = PrevF(NextE(a));
                Place_t d = PrevF(PrevE(a));
                ??? e == PrevF(a);         
                {
                  Visit(b);
                  Visit(c);
                  Visit(d);
                  Visit(e);
                }
              }
              processed++;    
            }
          }
        }  
      }
      if (conflicts != 0)
      {
        fprintf(stderr, "CollectTetrahedra: conflicts == " & Fmt.Int(conflicts) & "\n");
      }
      assert(processed == top.cell.nel);
      return rt;
    }
  } /* END CollectTetrahedra */
    
AdjacencyMatrix *MakeAdjacencyMatrix(ElemTableRec_t *top)==
      m = NEW(REF ARRAY OF ARRAY OF bool_t, top->node.nel, top->node.nel);
  { for (i = 0; i < top->node.nel; i++){ 
      ??? mi = m[i];
    { 
        for (j = 0; j < top->node.nel; j++)
        {
          mi[j] = FALSE
        }
      }
    }
    for (ei = 0; ei < top->wall.nelE; ei++)
    {
      ??? a = top->wedge[ei];
    {
        ??? i = OrgV(a)->num, j == OrgV(Clock(a))->num;
    { 
          m[i,j] = TRUE;
          m[j,i] = TRUE;
        }
      }
    }
    return m;
  } /* END MakeAdjacencyMatrix */

bool_t TriviallyIsomorphic(ElemTableRec_t *ta, ElemTableRec_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 = OrgV(sa)->num;
            ??? nb = OrgV(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(ElemTableRec_t *top)
  { for (i = 0; i < top->node.nel; i++){ 
      Node_t v = OrgV(top->node.el[i]);
        ??? e = top->node[i];
      {
        assert(v->num == i);
        assert(OrgV(e) == v);
      }
    }
    for (i = 0; i < top->cell.nel; i++)
    {
      ??? p = top->cell[i];
        ??? r = Srot(top.cell[i]);
      {
        assert(p->num == i);
        assert(OrgV(r) == p);
      } 
    }
  } /* END CheckOutAndRegion */

PROCEDURE GetVariableNodes(ElemTableRec_t *top,
                              VAR vr: ARRAY OF bool_t)
  { for (i = 0; i < top->node.nel; i++){
      Node_t v = OrgV(top->node.el[i]);
    {
        vr[i] = NOT v.fixed; 
      }
    }
  } /* END GetVariableNodes */
  
  /* ===== GEOMETRIC TOOLS ===== */

r4_t Barycenter(ElemTableRec_t *top, Coords_t *c, bool_t all)
      b: r4_t = (r4_t){0.0, ..};
      uint n = 0;
  { for (i = 0 TO (c.nel-1)){ 
      Node_t v = OrgV(top->node.el[i]);
    {
        if ((all) || (v->exists)){ b = r4_Add(b, c[i]); N++; }
      }
    }
    if (N == 0)
    {
      return b;
    }
    else
    {
      return r4_Scale(1.0/n, b);
    }
  } /* END Barycenter */

r4_t PartialBarycenter(*v: Node_vec_t; *c: Coords_t; bool_t all)
      b: r4_t = (r4_t){0.0, ..};
      uint n = 0;
  { for (k = 0 TO (c.nel-1)){ 
      ??? vk = v[k], i == vk->num;
    {
        if ((all) || (vk->exists)){ b = r4_Add(b, c[i]); N++; }
      }
    }
    if (N == 0)
    {
      return b;
    }
    else
    {
      return r4_Scale(1.0/n, b);
    }
  } /* END PartialBarycenter */

double MeanNodeDistance(ElemTableRec_t *top, Coords_t *c,
    *ctr: r4_t;
    bool_t all;
 )
      S: double = 0.0;
      uint n = 0;
  { for (i = 0 TO (c.nel-1)){ 
      Node_t v = OrgV(top->node.el[i]);
    {
        if ((all) || (v->exists))
        {
          ??? d2 = r4_DistSqr(ctr, c[i]);
    {
            S = S + d2
          }
          N++
        }
      }
    }
    if (N == 0)
    {
      return 1.0;
    }
    else
    {
      return sqrt(S/n);
    }
  } /* END MeanNodeDistance */

double MaxNodeDistance(ElemTableRec_t *top, 
    Coords_t *c; 
    *ctr: r4_t;
    bool_t all;
 )
      radius = 0.0;
  { for (i = 0; i < top->node.nel; i++){
      Node_t v = OrgV(top->node.el[i]);
    {
        if ((all) || (v->exists))
        {
          ??? d = r4_Dist(ctr, c[i]);
          {
            if (d > radius){ radius = d;}
          }
        }
      }
    }            
    return radius;
  } /* END MaxNodeDistance */

double MeanThickness(ElemTableRec_t *top, Coords_t *c, r4_t *ctr, r4_t *norm, bool_t all)
      S : double = 0.0;
      uint n = 0;
  { for (i = 0; i < top->node.nel; i++){
      Node_t v = OrgV(top->node.el[i]);
    {
        if ((all) || (v->exists))
        {
          ??? d = r4_dot(r4_Sub(c[i],ctr), norm);
    {
            S = S + d*d;
          }
          N++
        }
      }
    }
    if (N == 0)
    {
      return 1.0;
    }
    else
    {
      return sqrt(S/n);
    }
  } /* END MeanThickness */

double MeanEdgeLength(ElemTableRec_t *top, Coords_t *c, bool_t all)
      S: double = 0.0;
      uint n = 0;
  { for (i = 0; i < top->NE; i++){ 
      ??? e = top->edge[i];
    {
        if ((all) || (e->exists))
        {
          Node_t o = OrgV(e.pa)->num, d == OrgV(Clock(e.pa))->num;
    {
            S = S + r4_DistSqr(c[o], c[d])
          }
          N++;
        }
      }
    }
    if (N == 0)
    {
      return 1.0;
    }
    else
    {
      return sqrt(S/n);
    }
  } /* END MeanEdgeLength */

r4_t MeanCellNormal(ElemTableRec_t *top, Coords_t *c, bool_t all)
      norm: r4_t = (r4_t){0.0, ..};
  { for (i = 0; i < top->cell.nel; i++){
      Place_t f = Tors(Srot(top.cell[i])), p == PnegP(f);
    {
        if ((all) || ((p!=NULL) && (p->exists)))
        {
          ??? n = PolyCross(f, c);
    {
            norm = r4_Add(norm, n)
          }
        }
      }
    }
    ??? m = r4_Norm(norm);
    {
      if (m < 1.0d-20)
      {
        return (r4_t){0.0, 0.0, 0.0, 1.0};
      }
      else
      {
        return r4_Scale(1.0/m, norm);
      }
    }
  } /* END MeanCellNormal */

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

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

void NormalizeNodeDistances(ElemTableRec_t *top, Coords_t *c; bool_t all)
  { ??? b = Barycenter(top, c, all);
    {
      Displace(top, r4_Neg(b), c)
    }
    double s = MeanNodeDistance(top, c, (r4_t){0.0, ..}, all);
    {
      Scale(top, 1.0/s, c)
    }
  } /* END NormalizeNodeDistances */

void NormalizeEdgeLengths(ElemTableRec_t *top, Coords_t *c; bool_t all)
  { ??? b = Barycenter(top, c, all);
    {
      Displace(top, r4_Neg(b), c) 
    }
    ??? s = MeanEdgeLength(top, c, all);
    {
      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 = r4_cross(r4_Sub(p,o), r4_Sub(r,o), r4_Sub(d,o));
        ??? n2 = r4_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 r4_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 = r4_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,
                      ElemDataRec_t *top)
      
    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= r4_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, ElemTableRec_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 */
  
Node_vec_t *Neighbors(
    Node v;
    *ElemTableRec_t *top;
  )==
  REF Node_vec_t  rv = NULL;
    nv: uint = 0;
    
  void Stack(u: Node)
    {
      if ((rv == NULL) || (nv >= ((rv^).nel)))
      {
        ??? sv = NEW(REF Node_vec_t, MAX(10, 2*nv));
    {
          if (rv != NULL){ SUBARRAY(sv^, 0, nv) = rv^; }
          rv = sv
        }
      }
      rv[nv] = u; nv++;
    } Stack;
    
  void Check(b: Place_t)
    {
      if ((OrgV(b) == v)){ Stack(OrgV(Clock(b)));}
    } Check;

  { for (i = 0; i <= (top->NE - 1); i++){
      ??? ei = top->edge[i].pa;
    {
        Check(ei); Check(Clock(ei))
      }
    }
      
    Node_vec_trim(&rv,nv);
  } /* END Neighbors */

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

  PROCEDURE NodePoly(ElemDataRec_t *top):  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->cell.nel);
    {
      for (i = 0; i < top->cell.nel; i++)
      {
        ??? da = Srot(top.cell[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(PnegP(a0)->num == i);
            if ((PnegP(b0)!=NULL))
            {
              assert(PnegP(a0) == PnegP(b0));
            }
            assert(a0!=b0);
            poly[i] = Quadp{a0,a1,a2,a3};     
          }
        }
      }
      return poly;
    } NodePoly;

      poly1 : REF Quadv_vec_t = Quadv_vec_new(top->cell.nel);
      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(
     *ElemTableRec_t *top; 
     Coords_t *c;
    )==
  { ??? rvn = r4_vec_new(top->node.nel);
    ??? vn = rvn^;
    {
      for (i = 0; i < top->node.nel; i++)
      {
        vn[i] = NodeNormal(top->node[i], c, top)
      }
      return rvn
    }
  } /* END ComputeAllNodeNormals */ 

r4_vec_t *ComputeAllEdgeNormals(ElemTableRec_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(ElemTableRec_t *top, 
    Coords_t *c;
   )==
  { ??? rvn = r4_vec_new(top->wall.nel);
    ??? vn = rvn^;
    {
      for (i = 0; i < top->wall.nel; i++)
      {
        vn[i] = WallNormal(top->wall[i].pa, c);
      }
      return rvn
    }
  } /* END ComputeAllWallNormals */

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

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

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

      /* Gather nodes: */
      ptop.node.nel = 2 + ptop.edge.nel - ptop.wall.nel;
      ptop.vRef = Place_vec_new(ptop.node.nel);
      nv = 0;
      for (i = 0; i < ptop.edge.nel; 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;
            assert(PnegP(ptop.vRef[nv]) == pneg);
            nv++;
            u.xmark = TRUE;
          }
          if (! v.xmark){ ptop.vRef[nv] = ev;
            assert(PnegP(ptop.vRef[nv]) == pneg);
            nv++;
            v.xmark = TRUE; }
          PEdge(eu).xmark = FALSE;
        }
      }
      for (i = 0; i < ptop.node.nel; i++){ Node_t u = OrgV(ptop.vRef[i]);
    { u.xmark = FALSE; } }
    }
    return ptop;
  } /* END MakeCellTopology */

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

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

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

    /* Read node data Materials: */
    EVAL filefmt_read_comment(rd, '|');
    for (j = 0; j < top->node.nel; j++)
    {
      Lex.Skip(rd);
      int nv = fget_int32(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.el[0] = Lex.Real(rd);
          Lex.Skip(rd);
          cc.el[1] = Lex.Real(rd);
          Lex.Skip(rd);
          cc.el[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 filefmt_read_comment(rd, '|');
    for (j = 0; j < top->NE; j++)
    {
      Lex.Skip(rd);
      ??? ne = fget_int32(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.el[0] = Lex.Real(rd);
          Lex.Skip(rd);
          cc.el[1] = Lex.Real(rd);
          Lex.Skip(rd);
          cc.el[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 = fgetc(rd);
        {
          if (n!='-')
          {
                Rd.UnGetChar(rd);
            e->root = fget_int32(rd);
          }
          else
          {
            e->root = -1;
          }
        }
      }
    }
    Lex.Skip(rd);

    /* Read wall data materials: */
    EVAL filefmt_read_comment(rd, '|');
    for (j = 0; j < top->wall.nel; j++)
    {
      Lex.Skip(rd);
      ??? nf = fget_int32(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.el[0] = Lex.Real(rd);
          Lex.Skip(rd);
          cc.el[1] = Lex.Real(rd);
          Lex.Skip(rd);
          cc.el[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 = fgetc(rd);
        {
          if (n!='-')
          {
                Rd.UnGetChar(rd);
            f->root = fget_int32(rd);
          }
          else
          {
            f->root = -1;
          }
        }
      }
    }
    Lex.Skip(rd);

    /* Read cell data materials: */
    if (top->cell.nel!=0)
    {
      EVAL filefmt_read_comment(rd, '|');
    }
    for (j = 0; j < top->cell.nel; j++)
    {
      Lex.Skip(rd);
      ??? np = fget_int32(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.el[0] = Lex.Real(rd);
          Lex.Skip(rd);
          cc.el[1] = Lex.Real(rd);
          Lex.Skip(rd);
          cc.el[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);
        EVAL Mis.ReadBool(rd); /* p.degenerate */
        if (ro_te)
        {
          Lex.Skip(rd);
          ??? n = fgetc(rd);
          {
            if (n!='-')
            {
              Rd.UnGetChar(rd);
              p->root = fget_int32(rd);
            }
            else
            {
              p->root = -1;
            }
          }
        }
      }
    }

    ReadFooter(rd,"materials");
    fclose(rd);
    CheckOutAndRegion(top);
  } /* END ReadMaterials */

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

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

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

      void WriteLabel(char *label)
        {
          fprintf(ma, label);
        } WriteLabel;

      {
        WriteHeader(ma,"materials","99-08-25");
        if (((cmt != NULL) && ((*cmt) != 0)))
        {
          filefmt_write_comments(ma, cmt & "\n", '|')
        }
        ??? m = digits(top->cell.nel);
        {
          filefmt_write_comments(ma,"nodes    " & Fmt.Pad(Fmt.Int(top->node.nel),m), '|');
          filefmt_write_comments(ma,"@{edge->?}s       " & Fmt.Pad(Fmt.Int(top->NE),m), '|');
          filefmt_write_comments(ma,"walls       " & Fmt.Pad(Fmt.Int(top->wall.nel),m), '|');
          filefmt_write_comments(ma,"cells   " & Fmt.Int(top->cell.nel), '|');
        }
        filefmt_write_comments(ma, "\nNode data:\n", '|');

        for (i = 0; i < top->node.nel; i++)
        {
          Node_t v = OrgV(top->node.el[i]);
    {
            /* materials */
            fprintf(ma, Fmt.Pad(Fmt.Int(v->num), vWidth)); 
            fprintf(ma, "  ");
            fprintf(ma, Fmt.Char(BoolChars[v->exists]));
            fprintf(ma, " ");
            fprintf(ma, Fmt.Char(BoolChars[v.fixed]));
            fprintf(ma, "  ");
            WriteColor(v.color);
            fprintf(ma, "  ");
            WriteColor(v.transp);
            fprintf(ma, "  ");
            WriteRadius(v.radius);
            fprintf(ma, "  ");
            WriteLabel(v.label);
            fprintf(ma, "\n");
          }       
        }

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

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

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

        if (top->cell.nel!=0)
        {
          filefmt_write_comments(ma, "\nCell data:\n", '|');
        }

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

void Collectwedges(Place_t @p, Place_vec_t *re, uint *ne)
  Place_t *b;
  { b = a; ne = 0;
    do 
    {
      assert(PWall(b) == PWall(a));
      if ((re == NULL) || (ne >= ((re^).nel)))
      {
        ??? se = NEW(REF Place_vec_t, MAX(10, 2*ne));
    {
          if (re != NULL){ SUBARRAY(se^, 0, ne) = re^; }
          re = se
        }
      }
      re[ne] = b; ne++;
      b = NextE(b)
    } while (b != a);
  } /* END Collectwedges */
  
void CollectCellWalls(Place_t @p, Place_vec_t *rf, uint *nf)

      scanned: uint = 0;
      
  void StackWall(Place_t c)
    {
      assert(PnegP(c) == PnegP(a));
      PWall(c).xmark = TRUE;
      if ((rf == NULL) || (nf >= ((rf^).nel)))
      {
        ??? sf = NEW(REF Place_vec_t, MAX(10, 2*nf));
    {
          if (rf != NULL){ SUBARRAY(sf^, 0, nf) = rf^; }
          rf = sf
        }
      }
      rf[nf] = c; nf++;
    } StackWall;

  void VisitWall(b: Place_t)
    Place_t c = Clock(PrevF(b));
    VAR {
      do 
      {
        assert(PnegP(c) == PnegP(a));
        Wall_t g = PWall(c);
    {
          if (! g.xmark)
          {
            StackWall(c);
            /* We have only one mark per wall, but a wall may incide twice on "p": */
            /* Assumes that the cell's boundary is oriented! */
            if ((PposP(c) == PnegP(a))){ StackWall(Clock(c)); }
          }
        }
        c = NextE(c)
      } while (c != b);
    } VisitWall;
    
  { nf = 0; scanned = 0;
    StackWall(a);
    while (scanned < nf)
    {
      VisitWall(rf[scanned]);
      scanned++
    }
    /* Clear {xmark}, just in case: */
    for (k = 0; k <= (nf - 1); k++)
    {
      ??? b = rf[k];
    { PWall(b).xmark = FALSE;}
    }
  } /* END CollectCellWalls */

void CollectCellEdges(Place_t @p, 
    VAR Place_vec_t re; 
    uint *ne;
    VAR Place_vec_t rf;  /* Working area */
 ) ==
  uint *nf;
      Place_t b;
  { CollectCellWalls(a, rf, nf);
    ne = 0;
    for (k = 0; k < nf; k++)
    {
      ??? f = rf[k];
    {
        assert(PnegP(f) == PnegP(a));
        b = f;
        do 
        {
          Edge_t e = PEdge(b);
    { 
            if (! e.xmark)
            {
              if ((re == NULL) || (ne >= ((re^).nel)))
              {
                ??? se = NEW(REF Place_vec_t, MAX(10, 2*ne));
    {
                  if (re != NULL){ SUBARRAY(se^, 0, ne) = re^; }
                  re = se
                }
              }
              e.xmark = TRUE;
              re[ne] = b; ne++
            }
          }
          b = NextE(b)
        } while (b != f);
      }
    }
    /* Clear {xmark} bits: */
    for (k = 0; k <= (ne - 1); k++)
    {
      ??? b = re[k];
    { PEdge(b).xmark = FALSE;}
    }
  } /* END CollectCellEdges */

void CollectCellNodes(Place_t @p, Place_vec_t *rv, uint *nv, Place_vec_t *rf /* Working area */)
  uint *nf;
      Place_t b;
  { CollectCellWalls(a, rf, nf);
    nv = 0;
    for (k = 0; k <= (nf - 1); k++)
    {
      ??? f = rf[k];
    {
        assert(PnegP(f) == PnegP(a));
        b = f;
        do 
        {
          Node_t v = OrgV(b);
    { 
            if (! v.xmark)
            {
              if ((rv == NULL) || (nv >= ((rv^).nel)))
              {
                ??? sv = NEW(REF Place_vec_t, MAX(10, 2*nv));
    {
                  if (rv != NULL){ SUBARRAY(sv^, 0, nv) = rv^; }
                  rv = sv
                }
              }
              v.xmark = TRUE;
              rv[nv] = b; nv++
            }
          }
          b = NextE(b)
        } while (b != f);
      }
    }
    /* Clear {xmark} bits: */
    for (k = 0; k < nv; k++)
    {
      ??? b = rv[k];
    { OrgV(b).xmark = FALSE;}
    }
  } /* END CollectCellNodes */

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

  void PrintDualPlace_t(FILE *wr, Place_t @p; feWidth: uint)
    {
      fprintf(wr,Fmt.Pad(Fmt.Int(PWedge(a)->num), feWidth) & ":" \
        Fmt.Int((Octf.SrotBits(a)+ 3) % 4) & ":" & Fmt.Int(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); fprintf(wr, " "); 
        b = Srot(b)
      }
    } PrintDualwedge;  

  { ??? tp = FileWr.Open(name & ".tp");
    ??? vWidth = digits(MAX(1, top->node.nel - 1));
    ??? eWidth = digits(MAX(1,top->NE - 1));
    ??? fWidth = digits(MAX(1,top->wall.nel - 1));
    ??? pWidth = digits(MAX(1,top->cell.nel - 1));
    ??? feWidth = digits(MAX(1,top->wall.nelE -1));
    {
      WriteHeader(tp,"topology","99-08-25");
      if (((cmt != NULL) && ((*cmt) != 0)))
      {
        filefmt_write_comments(tp, cmt & "\n", '|')
      }

      ??? m = digits(top->wall.nelE);
      {
        fprintf(tp, "nodes   ");
        fprintf(tp, Fmt.Pad(Fmt.Int(top->cell.nel),m) & "\n");  
        fprintf(tp, "@{edge->?}s      ");
        fprintf(tp, Fmt.Pad(Fmt.Int(top->wall.nel),m) & "\n");
        fprintf(tp, "walls      ");
        fprintf(tp, Fmt.Pad(Fmt.Int(top->NE),m) & "\n"); 
        fprintf(tp, "cells  ");
        fprintf(tp, Fmt.Pad(Fmt.Int(top->node.nel),m) & "\n");
        fprintf(tp, "wedges ");
        fprintf(tp, Fmt.Int(top->wall.nelE) & "\n");
        fprintf(tp, "der        ");
        fprintf(tp, Fmt.Pad(Fmt.Int(1),m) & "\n");
        fprintf(tp, "bdr        ");
        fprintf(tp, Fmt.Pad(Fmt.Int(0),m) & "\n");
      }

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

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

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

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

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

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

      for (i = 0; i < NV; i++)
      {
        fprintf(wr, Fmt.Pad(Fmt.Int(i), d) & ": "); 
        WritePoint(c[i]); 
        fprintf(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
      )){
          fprintf(wr, "0 0 0 0")
        }
        else
        {
          WritePoint(Dc[i])
        } 
        fprintf(wr, "\n");
      }

      filefmt_write_footer(wr, "state-derivatives");
      fprintf(wr, "\n");
      fflush(wr);
    }
  } /* END WriteStDe */

{
} Triangulation.

  /*


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


*/