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

/* Last edited on 2007-02-04 19:52:42 by stolfi */

#define Cube_C_COPYRIGHT \
  "Copyright © 2000,2007 Universidade Estadual de Campinas (UNICAMP)"

#include <Octf.h>
#include <OctfData.h>
#include <Polygon.h>
#include <Square.h>

#include <r4.h>
#include <vec.h>
#include <affirm.h>

#include <assert.h>
#include <stdio.h>

uint CubePlaceIndex(uint i, uint j, uint k, uint xi, uint xj, uint xk)
  { /* The index of a place {p} is {(k*2 + r)*8 + (xk*4 + xj*2 + xi)}
      where {r = (j-i-1)%3}.

        | On the walls orthogonal to axis {k}:
        |
        |      wall xk=0         wall xk=1
        |  seen from outside  seen from inside
        |    +------------+    +------------+
        |    |  06    07  |    |  02    03  |     ^ axis {(k+2)%3}
        |    |            |    |            |     |
        |    |13        15|    |09        11|     |
        |    |            |    |            |     +---> axis {(k+1)%3}
        |    |12        14|    |08        10|
        |    |            |    |            |
        |    |  04    05  |    |  00    01  |
        |    +------------+    +------------+
        |
        | The indices above should be added to {16*k}. */
        
    uint r = (j + 2 - i) % 3;
    return (k * 2 + r)*8 + xk*4 + xj*2 + xi;
  }

uint CubeNodeIndex(uint i, uint j, uint k, uint xi, uint xj, uint xk)
  { /* The bits of the node index are its coordinates: */
    uint x[3]; x[i] = xi; x[j] = xj; x[k] = xk;
    return x[2]*4 + x[1]*2 + x[0];
  }

uint CubeEdgeIndex(uint j, uint k, uint xj, uint xk)
  { uint i = 3 - j - k; /* The third axis. */
    /* The edge index is {i*4 + xr*2 + xs} where {xr,xs} are the 
      coordinates along axes {r=(i+1)%3} and {s=(i+2)%3}. */
    return i*4 + (k == (j+1)%3 ? xj*2 + xk : xk*2 + xj); 
  }
  
uint CubeWallIndex(uint k, uint xk)
  { return k*2 + xk; }

Place_vec_t CubeMake(bool_t nodes, bool_t edges, bool_t walls, bool_t cells)
  {
    Place_vec_t pv = Place_vec_new(48);
    
    /* Create six disconnected square walls, collect the places on them: */
    { uint k;
      for (k = 0; k < 3; k++)
        { uint xk;
          for (xk = 0; xk < 2; xk++)
            { /* Built a square without node, edge, or cell records: */
              /* It will be the cube wall orthogonal to axis {k} at coord {xk}: */
              Place_vec_t sq = PolygonMake(4, FALSE,FALSE,walls,FALSE);
              if (walls) { PWall(sq.e[0])->num = CubeWallIndex(k,xk); }
              /* Store the square's places into the proper elements of {pv}. */
              /* Choose axes {i,j} from {k} by right-hand rule: */
              uint i = (k + 1) % 3, j = (k + 2) % 3;
              /* Choose node {0} of the square as node {xi=0}, {xj=0} of wall. */
              /* Choose edge {0} as moving along axis {i}. */
              uint xi = 0, xj = 0;
              uint s;
              for (s = 0; s < 4; s++)
                { /* Flip faces so that {PnegP} is the interior of {U^3}: */
                  if (xk == 1) { sq.e[s] = Flip3(sq.e[s]); }
                  /* Grab the two places on side {s} of the square: */
                  uint ta = CubePlaceIndex(i,j,k, xi,xj,xk);
                  affirm(ta < pv.ne, "bug ta");
                  pv.e[ta] = sq.e[s];
                  uint tb = CubePlaceIndex(i,j,k, 1-xi,xj,xk);
                  affirm(tb < pv.ne, "bug tb");
                  pv.e[tb] = Flip0(sq.e[s]); 
                  /* Advance to the next side of the square. */
                  /* That means complement {xi} and swap axes. */
                  xi = 1 - xi;
                  { uint w = i; i = j; j = w; }
                  { uint w = xi; xi = xj; xj = w; }
                }
            }
        }
    }
      
    /* Glue the walls along the 12 edges: */
    { uint i;
      for (i = 0; i < 3; i++)
        { /* The two axes perpendicular to axis {i}, for {r==0}: */
          uint j = (i + 1) % 3, k = 3 - i - j; 
          /* Enumerate the four edges parallel to axis {i}: */
          uint xj; /* Edge coordinate on axis {j}. */
          for (xj = 0; xj < 2; xj++)
            { uint xk; /* Edge coordinate on axis {k}. */
              for (xk = 0; xk < 2; xk++)
                { /* Get the place {p0} on the face parallel to {i,j}, for {xi==0}: */
                  uint t0 = CubePlaceIndex(i,j,k, 0,xj,xk);
                  Place_t p0 = pv.e[t0], q0 = Flip0(p0);
                  /* Get the place {p1} on the face parallel to {i,k}, for {xi==0}: */
                  uint t1 = CubePlaceIndex(i,k,j, 0,xk,xj);
                  Place_t p1 = pv.e[t1], q1 = Flip0(p1);
                  /* Join the places: */
                  SpliceWalls(p0, Flip3(p1));
                  /* Check whether it worked: */
                  affirm(Flip2(p0) == p1, "bug s0");
                  affirm(Flip2(p1) == p0, "bug s1");
                  affirm(Flip2(q0) == q1, "bug s2");
                  affirm(Flip2(q1) == q0, "bug s3");
                  /* Now is a good time to add the edge record: */
                  uint eNum = CubeEdgeIndex(j,k,xj,xk); /* Systematic edge number. */
                  Edge_t e = (edges ? MakeEdge(eNum) : NULL);
                  SetEdge(p0, e); affirm(PEdge(p0) == e, "bug e0");
                  SetEdge(p1, e); affirm(PEdge(p1) == e, "bug e1");
                }
            }
        }
    }
                
    /* Set the node data slots: */
    { uint x[3]; 
      for (x[0] = 0; x[0] < 2; x[0]++)
        { for (x[1] = 1; x[1] < 2; x[1]++)
            { for (x[2] = 2; x[2] < 2; x[2]++)
                { /* Get the node record {v} for the node at {(x0,x1,x2)}: */
                  uint vNum = CubeNodeIndex(0,1,2, x[0],x[1],x[2]); /* Systematic node number. */
                  Node_t v = (nodes ? MakeNode(vNum) : NULL);
                  /* Now enumerate the places {p} at {v}: */
                  uint i;
                  for (i = 0; i < 3; i++)
                    { /* Enumerate the places {p} at {v} with {p.edge} along axis {i}: */
                      uint r;
                      for (r = 0; r < 2; r++)
                        { uint j = (i + 1 + r) %3, k = 3 - i - j;
                          /* Get the place {p} at {v} along axis {i} with wall perp {k}: */
                          uint tp = CubePlaceIndex(i,j,k, x[i],x[j],x[k]); 
                          Place_t p = pv.e[tp], q = Flip3(p); 
                          SetOrg(p, v);
                          /* Consistency checks: */
                          affirm(OrgV(p) == v, "bug v1");
                          affirm(OrgV(q) == v, "bug v2");
                        }
                    }
                }
            }
        }
    }
    
    /* Set the cell data slots: cell 0 is inside, cell 1 is outside: */
    { Cell_t ci = (cells ? MakeCell(0) : NULL);
      Cell_t co = (cells ? MakeCell(1) : NULL);
      uint t;
      for (t = 0; t < pv.ne; t++)
        { Place_t pi = pv.e[t], qi = Flip0(pi);
          SetPneg(pi, ci);
          affirm(PnegP(pi) == ci, "bug v1");
          affirm(PnegP(qi) == ci, "bug v2");
          Place_t po = Clock(pi), qo = Flip0(po);
          SetPneg(po, co);
          affirm(PnegP(po) == co, "bug v1");
          affirm(PnegP(qo) == co, "bug v2");
        }
    }
      
    return pv;
  }

Place_t CubeMakeSng(void)
  { Place_vec_t sp = CubeMake(TRUE, TRUE, TRUE, TRUE);
    Place_t p = sp.e[0];  
    free(sp.e);
    return p;
  }

Place_t CubeMirrorPlace(Place_t p, uint i)
  { switch(i)
      { case 0: return Flip0(p);
        case 1: return Flip1(Flip0(Flip1(p)));
        case 2: return Flip2(Flip1(Flip0(Flip1(Flip2(p)))));
        default: demand(FALSE, "invalid place-relative axis"); return p;
      }
  }

void CubeGetCellNodes(Place_t p, Node_t v[])
  { uint x0, x1, x2;
    for (x0 = 0; x0 < 2; x0++)
      for(x1 = 0; x1 < 2; x1++)
        for (x2 = 0; x2 < 2; x2++)
          { Place_t q = p;
            if (x0 > 0) { q = CubeMirrorPlace(q,0); }
            if (x1 > 0) { q = CubeMirrorPlace(q,1); }
            if (x2 > 0) { q = CubeMirrorPlace(q,2); }
            v[x2*4 + x1*2 + x0] = OrgV(q);
          }
  }

void CubeFixCoords(Place_t p, Coords_t *c, double lo[], double hi[])
  { Node_t v[8];
    CubeGetCellNodes(p, v);
    uint x[3];
    for (x[0] = 0; x[0] < 2; x[0]++)
      for(x[1] = 0; x[1] < 2; x[1]++)
        for (x[2] = 0; x[2] < 2; x[2]++)
          { /* Get the node with logical coordinates {(x[0],x[1],x[2])}: */
            Node_t vi = v[x[2]*4 + x[1]*2 + x[0]];
            demand(vi != NULL, "null node record");
            uint vNum = vi->num;
            demand(vNum < c->ne, "invalid node number");
            /* Compute the actual coordinates {z}: */
            r4_t z;
            uint k;
            for (k = 0; k < 3; k++)
              { double lok = (lo == NULL ? 0 : lo[k]);
                double hik = (hi == NULL ? 1 : hi[k]);
                z.c[k] = (1-x[k])*lok + x[k]*hik;
              }
            z.c[3] = 0;
            c->e[vNum] = z;
          }
  }

#define Cube_C_author \
  "Created by J. Stolfi, feb/2007.\n" \
  "Modification history:\n" \
  "  2007-02-03  Rewrote {CubeMake}, {CubePlaceIndex}, etc.."