/* Last edited on 2007-01-31 17:54:06 by stolfi */

typedef struct Place_t { 
  struct Wedge_t walledge; /* The base wall-@{edge->?}. */
  SRBits_t bits;                    /* The spin and rotation bits to apply to it. */
} Place_t;

void PrtPnegPpos(Place_t @p)
  /* Print the cell negative and positive of places
     belonging to same wall. */
       Place_t t = a;
  { do {
      Tools.PrtPnegNum(t);
      Tools.PrtPposNum(t);   
      Tools.PrtPnegNum(Clock(t));
      Tools.PrtPposNum(Clock(t));
      t = PrevE(t);
    } while (t != a);
  } /* END PrtPnegPpos */

void PrtNextPnegC(Place_t p)
  /* Print the 24 places belonging to same cell negative
     for one cubic cell. */

  void PrintPlacesOnWall(Place_t @p)
    Place_t t = a;
    VAR {
      do 
      {
        Tools.PrtPnegNum(t);
        t = PrevE(t);
      } while (t != a);
  } /* END PrintPlacesOnWall */

       Place_t t = a;
  { PrintPlacesOnWall(t);
    do 
    {
      PrintPlacesOnWall(Clock(PrevF(t)));
      t = PrevE(t);
    } while (t != a);
    t = PrevF(PrevE(PrevE(PrevF(a))));
    PrintPlacesOnWall(t);
  } /* END PrtNextPnegC */


void WriteStates
  ( FILE *wr,
    double time,
    Coords3D_t *c3,  /* Node coordinates */
    uint prec /* DF 4 */,    
    char *cmt /* DF " " */
  )
  { 
    auto void WriteCoord3(double x);
    
    void WriteCoord3(double x)
      { fprintf(wr, "%*.*e", prec+8, prec, x); }

    auto void WritePoint3(r3_t *p);
    
    void WritePoint3(r3_t *p)
      { WriteCoord3(p->c[0]); fprintf(wr, " ");
        WriteCoord3(p->c[1]); fprintf(wr, " ");
        WriteCoord3(p->c[2]); fprintf(wr, " ");
      }
    
    int NV = c3->nel;
    int vWidth = NumDigits(NV - 1);
    
    filefmt_write_header(wr, "state3D", ST3_FILE_VERSION);
    fprintf(wr, "nodes " & Fmt.Int(node.nel) & "\n");
    filefmt_write_comments(wr, cmt & "\n", '|');
    filefmt_write_comment(wr, "\nNode data:\n", '|');
    fprintf(wr,"t " \
               Fmt.Pad(Fmt.LongReal(time,Fmt.Style.Fix,prec=2),4)&"\n");
    for (i = 0 TO (node.nel-1))
    {
      fprintf(wr, Fmt.Pad(Fmt.Int(i), vWidth)); 
      fprintf(wr, " ");
      WritePoint3(node[i]); WriteEOL(wr);
    }
    WriteFooter(wr, "state3D");
    fclose(wr)
    
  } /* END WriteStates */


  /* Notice: actually WriteStates == WriteState3D

*/

 
uint ReadHeader(FILE *rd, char *type, char *param /* DF NULL */);
  /* Similar to procedure {FileFmt.ReadHeader}. */

void ReadFooter(FILE *rd, char *type);
  /* Similar to procedure {FileFmt.ReadFooter}. */

void WriteHeader(FILE *wr, char *type, char *param /* DF NULL */, uint n);
  /* Similar procedure to {FileFmt.WriteHeader}. */

/* ---------------------------------------------------------------------- */

#define INIT_STACK_SIZE 1024

Wedge_vec_t EnumWedges(Place_vec_t root, VisitFunc_t visit)
  { Place_vec_t stack = Place_vec_new(INIT_STACK_SIZE);
    uint top = 0;
    uint seen = 0; /*!=of @{quad->?}s whose childeren were looked at */

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

    void VisitAndMark(Place_t c)
      { Wedge_t w = PWedge(c);
        if (! w->marks2[DualBit(c)]) 
          { visit(c); /* [!!! Should we allow for early termination? !!!] */
            if (! wedges) { visit(Clock(c));};
            Place_vec_expand(&stack,top);
            w->marks2[DualBit(c)] = TRUE; 
            stack.el[top] = c; top++;
          }
      }

    assert(! PWedge(p)->marks2[DualBit(p)]); 
    /* Enumerate all wedges: */
    VisitAndMark(p);
    while (seen < top)
      { Place_t b = stack.el[seen];
        VisitAndMark(PrevF(b));
        VisitAndMark(PrevE(b));
        VisitAndMark(NextE(b));
        VisitAndMark(NextF(PrevE(b)));        
        seen++;
      }
    /* Erase all marks */
    while (top > 0)
      { top--;
        Place_t b = stack.el[top];
        PWedge(b)->marks2[DualBit(b)] = FALSE;
      }
  }

/* ---------------------------------------------------------------------- */
    
ElemTableRec_t MakeElemTable(Place_vec_t root)
  {
    ElemTableRec_t top;
    
    ...;
    
    uint i, k;
    for (i = 0; i < top.wedge.nel; i++)
      { Place_t c = PBase(top.wedge.el[i]);
        for (k = 0; k < 2; k++)
        { Node_t v = OrgV(c); 
          uint vi = v->num;
          top.node.el[vi] = c;
          Cell_t q = PnegP(c);
          if (q != NULL)
            { uint qi = q->num;
              top.cell.el[qi] = Srot(c); /* [!!! Should be just {c} !!!] */
              assert(Pneg(c) == Org(top.cell.el[qi]));
            }
          c = Clock(c); 
        }
      }
    return top;
  }