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

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

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

#include <Tools.h>
#include <Triangulation.h>
#include <Octf.h>
#include <OctfData.h>
#include <CubeLAPL.h>
#include <TriangulatedCube.h>

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

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

Place_t GlueBigCube(Place_t a, Place_t b, uint n);

Place_t GlueBigCube(Place_t a, Place_t b, uint n)
  { demand(n >= 1, "order must be at least 1");
    int M = 20;
    Place_t ta[M*M];
    Place_t tb[M*M];
    ta[0*M + 0] = a; tb[0*M + 0] = b;

    int i, j;
    for (i = 1; i < n; i++)
      { ta[i*M + 0] = Clock(PrevF(NextE(NextE(ta[(i-1)*M+0]))));
        tb[i*M + 0] = Clock(NextF(NextE(NextE(tb[(i-1)*M+0]))));
      }

    for (i = 0; i < n; i++)
      { for (j = 1; j < n; j++)
          { ta[i*M + j] = Clock(PrevE(PrevF(NextE(ta[i*M + (j-1)]))));
            tb[i*M + j] = Clock(PrevE(NextF(NextE(tb[i*M + (j-1)]))));
            assert(ta[i*M + j]!=a); 
            assert(tb[i*M + j]!=b);
          }
      }

    for (i = 0; i < n; i++)
      { Meld(tb[i*M + 0], ta[i*M + 0]);

        /* Update the edge slots */
        SetRingEdgeInfo(ta[i*M + 0], PWedge(ta[i*M + 0])->edge);
        SetRingEdgeInfo(NextE(ta[i*M + 0]), PWedge(NextE(ta[i*M + 0]))->edge);
        SetRingEdgeInfo(NextE(NextE(ta[i*M + 0])), PWedge(NextE(NextE(ta[i*M + 0])))->edge);
        SetRingEdgeInfo(NextE(NextE(NextE(ta[i*M + 0]))),
                   PWedge(NextE(NextE(NextE(ta[i*M + 0]))))->edge);

        /* Update the node slots */
        SetOrgAll(ta[i*M + 0], OrgV(ta[i*M + 0]));
        SetOrgAll(NextE(ta[i*M + 0]), OrgV(NextE(ta[i*M + 0])));
        SetOrgAll(NextE(NextE(ta[i*M + 0])), OrgV(NextE(NextE(ta[i*M + 0]))));
        SetOrgAll(NextE(NextE(NextE(ta[i*M + 0]))),
                   OrgV(NextE(NextE(NextE(ta[i*M + 0])))));
        SetOrgAll(Clock(ta[i*M + 0]), OrgV(Clock(ta[i*M + 0])));
        SetOrgAll(Clock(NextE(ta[i*M + 0])), OrgV(Clock(NextE(ta[i*M + 0]))));
        SetOrgAll(Clock(NextE(NextE(ta[i*M + 0]))),
                   OrgV(Clock(NextE(NextE(ta[i*M + 0])))));
        SetOrgAll(Clock(NextE(NextE(NextE(ta[i*M + 0])))),
                   OrgV(Clock(NextE(NextE(NextE(ta[i*M + 0]))))));

        SetPneg(ta[i*M + 0], PnegP(tb[i*M + 0]));
        SetPneg(PrevE(ta[i*M + 0]), PnegP(PrevE(tb[i*M + 0])));
        SetPneg(PrevE(PrevE(ta[i*M + 0])), PnegP(PrevE(PrevE(tb[i*M + 0]))));
        SetPneg(PrevE(PrevE(PrevE(ta[i*M + 0]))),
                PnegP(PrevE(PrevE(PrevE(tb[i*M + 0])))));

      }
    
    for (i = 0; i < n; i++)
      { for (j = 1; j < n; j++)
          { Meld(tb[i*M + j],ta[i*M + j]);

            /* Update the edge slots */
            SetRingEdgeInfo(ta[i*M + j], PWedge(ta[i*M + j])->edge);
            SetRingEdgeInfo(NextE(ta[i*M + j]), PWedge(NextE(ta[i*M + j]))->edge);
            SetRingEdgeInfo(NextE(NextE(ta[i*M + j])),
                                 PWedge(NextE(NextE(ta[i*M + j])))->edge);
            SetRingEdgeInfo(NextE(NextE(NextE(ta[i*M + j]))),
                       PWedge(NextE(NextE(NextE(ta[i*M + j]))))->edge);

            /* Update the node slots */
            SetOrgAll(ta[i*M + j], OrgV(ta[i*M + j]));
            SetOrgAll(NextE(ta[i*M + j]), OrgV(NextE(ta[i*M + j])));
            SetOrgAll(NextE(NextE(ta[i*M + j])), OrgV(NextE(NextE(ta[i*M + j]))));
            SetOrgAll(NextE(NextE(NextE(ta[i*M + j]))),
                                       OrgV(NextE(NextE(NextE(ta[i*M + j])))));
            SetOrgAll(Clock(ta[i*M + j]), OrgV(Clock(ta[i*M + j])));
            SetOrgAll(Clock(NextE(ta[i*M + j])), OrgV(Clock(NextE(ta[i*M + j]))));
            SetOrgAll(Clock(NextE(NextE(ta[i*M + j]))),
                                         OrgV(Clock(NextE(NextE(ta[i*M + j])))));
            SetOrgAll(Clock(NextE(NextE(NextE(ta[i*M + j])))),
                       OrgV(Clock(NextE(NextE(NextE(ta[i*M + j]))))));

            SetPneg(ta[i*M + j], PnegP(tb[i*M + j]));
            SetPneg(PrevE(ta[i*M + j]), PnegP(PrevE(tb[i*M + j])));
            SetPneg(PrevE(PrevE(ta[i*M + j])), PnegP(PrevE(PrevE(tb[i*M + j]))));
            SetPneg(PrevE(PrevE(PrevE(ta[i*M + j]))),
                    PnegP(PrevE(PrevE(PrevE(tb[i*M + j])))));                 
          }
      }

    return ta[0*M + 0];
  }

FourPlaces_vec_t Make1DCubeArray(uint order)
  { /* [!!! Unify with {Make1DCubeArrayBis} !!!] */
    SixPlaces_vec_t ca = SixPlaces_vec_new(order);
    FourPlaces_vec_t cb = FourPlaces_vec_new(order);
    int i, j;
    for (i = 0; i < order; i++)
      { ca.e[i] = CubeLAPLMake();
        cb.e[i].p[0] = ca.e[i].p[1]; 
        cb.e[i].p[1] = ca.e[i].p[3]; 
        cb.e[i].p[2] = ca.e[i].p[0]; 
        cb.e[i].p[3] = ca.e[i].p[5];
      }

    /* Glue the cubes: */
    for (j = 0; j < order-1; j++)
      { GlueCube(ca.e[j].p[2], PrevF(ca.e[j+1].p[3])); }
    
    return cb;
  }

FourPlaces_vec_t Make1DCubeArrayBis(uint order)
  { /* [!!! Unify with {Make1DCubeArray} !!!] */
    SixPlaces_vec_t ca = SixPlaces_vec_new(order);
    FourPlaces_vec_t cb = FourPlaces_vec_new(order);
    int i, j;
    for (i = 0; i < order; i++)
      { ca.e[i] = MakeCubeT(1);
        cb.e[i].p[0] = ca.e[i].p[2]; 
        cb.e[i].p[1] = ca.e[i].p[3]; 
        cb.e[i].p[2] = ca.e[i].p[4]; 
        cb.e[i].p[3] = ca.e[i].p[5];
      }

    /* Glue the cubes: */
    for (j = 0; j < order-1; j++)
      { GlueCubeT(NextE(PrevF(ca.e[j+1].p[4])), ca.e[j].p[0]); }
    return cb;
  }
 
Place_t Make1DCubeArrayBis2(uint order)
  { FourPlaces_vec_t s = Make1DCubeArrayBis(order);
    return s.e[0].p[0];
  }

TwoPlaces_vec_t Make2DCubeArray(uint order)
  { FourPlaces_vec_t cb = FourPlaces_vec_new(order* order);
    TwoPlaces_vec_t cc = TwoPlaces_vec_new(order*order);

    /* Assembling the array {cb} */
    int i,j,k;
    for (i = 0; i < order; i++)
      { FourPlaces_vec_t ca = Make1DCubeArray(order);
        for (j = 0; j < order; j++)
          { for (k = 0; k < 4; k++)
              { cb.e[i*order + j].p[k] = ca.e[j].p[k]; }
          }
      }

    /* Gluing the rows: */
    for (j = 0; j < (order-1); j++)
      { for (k = 0; k < order; k++)
          { GlueCube(cb.e[j*order + k].p[1], Clock(NextE(NextE(cb.e[(j+1)*order + k].p[0])))); }
      }

    /* Selecting the places to return: */
    for (i = 0; i < order; i++)
      { for (j = 0; j < order; j++)
        { for (k = 0; k < 2; k++)
            { cc.e[i*order + j].p[k] = cb.e[i*order + j].p[k+2]; }
        }
      }
      
    return cc;
  }

TwoPlaces_vec_t MakeBigCube(uint order)
  { FourPlaces_vec_t cd = FourPlaces_vec_new(order*order*order);
    TwoPlaces_vec_t ce = TwoPlaces_vec_new(order*order);
    int i,j,k,l;
    for (i = 0; i < order; i++)
      { TwoPlaces_vec_t cc = Make2DCubeArray(order);
        for (j = 0; j < order; j++)
          { for (k = 0; k < order; k++)
              { for (l = 0; l < 2; l++)
                  { cd.e[(i*order + j)*order + k].p[l] = cc.e[j*order + k].p[l]; }
              }
          }
      }

    /* Gluing: */
    for (i = 0; i < (order-1); i++)
      { for (k = 0; k < order; k++)
          { for (j = 0; j < order; j++)
              { Place_t pikj = cd.e[(i*order + k)*order +j].p[1];
                Place_t qikj = cd.e[((i+1)*order +k)*order + j].p[0];
                GlueCube(pikj, Clock(qikj));
              }
          }
      }

    /* Selecting the places for return */
    for (j = 0; j < order; j++)
      { for (k = 0; k < order; k++)
          { ce.e[j*order + k].p[0] = cd.e[(0*order + k)*order + j].p[0];
            ce.e[j*order + k].p[1] = cd.e[((order-1)*order + k)*order + j].p[1];
          }
      }
    return ce;
  }

Place_t MakeBigCube1(uint order)
  { TwoPlaces_vec_t a = MakeBigCube(order);
    return a.e[0*order + 0].p[0];
  } 

void SetCubeProperties(TwoPlaces_vec_t *a, uint n, ElemTableRec_t *top)
  { uint index = 0;
    Place_vec_t c1 = Place_vec_new(n);
    Place_vec_t c2 = Place_vec_new(n);
    Place_vec_t c3 = Place_vec_new(n);
    Place_vec_t c4 = Place_vec_new(n);

    /* Set the original nodes */
    SetExNode(Clock(NextE(a->e[0*n + 0].p[0])), TRUE);    /* node 1 */   
    SetExNode(Clock(a->e[(n-1)*n + 0].p[0]), TRUE);         /* node 2 */
    SetExNode(a->e[(n-1)*n + (n-1)].p[0], TRUE);              /* node 3 */
    SetExNode(a->e[(n-1)*n + 0].p[1], TRUE);                /* node 4 */
    SetExNode(Clock(a->e[(n-1)*n + (n-1)].p[1]), TRUE);       /* node 5 */
    SetExNode(Clock(NextE(a->e[0*n + (n-1)].p[1])), TRUE);  /* node 6 */
    SetExNode(PrevE(a->e[0*n + 0].p[1]), TRUE);           /* node 7 */
    SetExNode(PrevE(a->e[0*n + (n-1)].p[0]), TRUE);         /* node 8 */

    /* Set the original edgess */
    int i, j, k;
    for (i = 0; i < n; i++)
      { SetExEdge(PrevE(a->e[i*n + 0].p[1]), TRUE);
        SetExEdge(NextE(a->e[i*n + 0].p[0]), TRUE);
        SetExEdge(a->e[(n-1)*n + i].p[0], TRUE);
        SetExEdge(a->e[(n-1)*n + i].p[1], TRUE);
        SetExEdge(PrevE(a->e[(i)*n + (n-1)].p[0]), TRUE);
        SetExEdge(NextE(a->e[(i)*n + (n-1)].p[1]), TRUE);
        SetExEdge(NextE(NextE(a->e[0*n + i].p[1])), TRUE);
        SetExEdge(NextE(NextE(a->e[0*n + i].p[0])), TRUE);
      }

    /* computing @{edge->?}s on the colums */
    c1.e[0] = NextE(NextF(a->e[(n-1)*n + 0].p[0]));
    c2.e[0] = PrevE(NextF(a->e[(n-1)*n + (n-1)].p[0]));
    c3.e[0] = NextE(NextF(NextE(a->e[0*n + 0].p[0])));
    c4.e[0] = PrevE(NextF(PrevE(a->e[0*n + (n-1)].p[0])));

    for (k = 1; k < n; k++)
      { c1.e[k] = Clock(PrevE(PrevF(NextE(c1.e[k-1]))));
        c2.e[k] = Clock(NextE(PrevF(PrevE(c2.e[k-1]))));
        c3.e[k] = Clock(PrevE(PrevF(NextE(c3.e[k-1]))));
        c4.e[k] = Clock(NextE(PrevF(PrevE(c4.e[k-1]))));
      }

    /* Now set the original edgess */
    for (j = 0; j < n; j++)
      { SetExEdge(c1.e[j], TRUE);
        SetExEdge(c2.e[j], TRUE);
        SetExEdge(c3.e[j], TRUE);
        SetExEdge(c4.e[j], TRUE);
      }

    for (i = 0; i < top->wall.ne; i++)
      { Place_t f = top->wall.e[i];
        { if (WallOnBoundary(f)) { SetExWall(f, TRUE); index++; } }
      }
    assert(index == 6*n*n);
  } 

TwoPlaces_vec_t BuildBigRawCube(uint order)
  { TwoPlaces_vec_t cd = TwoPlaces_vec_new(order*order*order);
    int i,j,k,l;
    for (i = 0; i < order; i++)
      { TwoPlaces_vec_t cc = Make2DCubeArray(order);
        for (j = 0; j < order; j++)
          { for (k = 0; k < order; k++)
              { for (l = 0; l < 2; l++)
                  { cd.e[(i*order + j)*order + k].p[l] = cc.e[j*order + k].p[l]; }
              }
          }
      }

    /* Gluing: */
    for (i = 0; i < (order-1); i++)
      { for (k = 0; k < order; k++)
          { for (j = 0; j < order; j++)
              { Place_t pikj = cd.e[(i*order + k)*order + j].p[1];
                Place_t qikj = cd.e[((i+1)*order + k)*order + j].p[0];
                GlueCube(pikj, Clock(qikj));
              }
          }
      }
    
    return cd;
  }

void EmphasizeBigRawCubeOriginalElems(TwoPlaces_vec_t *cd, ElemTableRec_t *top, uint order)      
  { /* Selecting the places for emphasize */
    TwoPlaces_vec_t a = TwoPlaces_vec_new(order*order);
    int i, j, k;
    for (j = 0; j < order; j++)
      { for (k = 0; k < order; k++)
          { a.e[j*order + k].p[0] = cd->e[(0*order + k)*order + j].p[0];
            a.e[j*order + k].p[1] = cd->e[((order-1)*order + k)*order + j].p[1];
          }
      }

    /* Set all elements (nodes, edgess, walls, as non-existing */
    for (i = 0; i < top->wall.ne; i++)
      { Place_t p = top->wall.e[i];
        PWall(p)->exists = FALSE;
      }
    for (i = 0; i < top->edge.ne; i++)
      { Place_t p = top->edge.e[i];
        PEdge(p)->exists = FALSE;
      }
    for (i = 0; i < top->node.ne; i++)
      { Place_t p = top->node.e[i];
        OrgV(p)->exists = FALSE;
      }

    /* Now emphasize the original elements (nodes,@{edge->?}s,walls) */ 
    SetCubeProperties(&a, order, top);
  }

Place_t BuildBigRawCube1(uint order, bool_t original)
  { TwoPlaces_vec_t cd = BuildBigRawCube(order);
    Place_t rp = cd.e[0].p[0];
    /* ElemTableRec_t top = MakeElemTable(Place_vec_desc(&rp,1)); */
    /* if (original) { EmphasizeBigRawCubeOriginalElems(&cd, top, order); } */
    demand(! original, "option {original = TRUE} not implemented yet");
    return rp;
  }

void FixBigRawCubeCoords(TwoPlaces_vec_t *cd, ElemTableRec_t *top, uint order, Coords_t *c)
  {
    double zero = 0.0, one = 1.0;
    
    auto void SetCorner(Place_t e, uint i);
    auto void SetCorner1(Place_t e, r4_t oc);
    auto void SetCorner2(Place_t e, r4_t oc);
    auto void SetCorner3(Place_t e, r4_t oc);
    
    void SetCorner(Place_t e, uint i)
      { c->e[OrgV(e)->num] = (r4_t){{i, zero, zero, zero}}; }

    void SetCorner1(Place_t e, r4_t oc)
      { c->e[OrgV(e)->num] = (r4_t){{oc.c[0], oc.c[1]+one, oc.c[2], oc.c[3]}}; }

    void SetCorner2(Place_t e, r4_t oc)
      { c->e[OrgV(e)->num] = (r4_t){{oc.c[0], oc.c[1], oc.c[2]+one, oc.c[3]}}; }

    void SetCorner3(Place_t e, r4_t oc)
      { c->e[OrgV(e)->num] = (r4_t){{oc.c[0], oc.c[1]-one, oc.c[2], oc.c[3]}}; }

    Place_t p;
    int i, j, k;
    for (i = 0; i < order; i++)
      { SetCorner(Clock(NextE(cd->e[(((i)*order + 0)*order) + 0].p[0])), i);
        for (j = 0; j < order; j++)
          { p = Clock(NextE(cd->e[(((i)*order + j)*order) + 0].p[0]));
            for (k = 0; k < order; k++)
              { SetCorner1(NextE(p),c->e[OrgV(p)->num]);
                SetCorner2(NextE(NextE(p)),c->e[OrgV(NextE(p))->num]);
                SetCorner3(NextE(NextE(NextE(p))),c->e[OrgV(NextE(NextE(p)))->num]);
                p = NextE(PrevF(PrevF(cd->e[(((i)*order + j)*order) + k].p[0])));
              }
          }
      }

    /* finally */  
    SetCorner(PrevE(cd->e[((((order-1))*order + 0)*order) + 0].p[1]), order);
    for (j = 0; j < order; j++)
      { p = PrevE(cd->e[(((order-1)*order + j)*order) + 0].p[1]);
        for (k = 0; k < order; k++)
          {
            SetCorner1(NextE(p),c->e[OrgV(p)->num]);
            SetCorner2(NextE(NextE(p)),c->e[OrgV(NextE(p))->num]);
            SetCorner3(NextE(NextE(NextE(p))),c->e[OrgV(NextE(NextE(p)))->num]);
            p = Clock(PrevE(NextF(cd->e[(((order-1)*order + j)*order) + k].p[1])));
          }
      }
  }

#define BigCube_C_author \
  "Created by L. A. P. Lozada, 2000.\n" \
  "Modification history:\n" \
  "  30-08-2000 : [???] Nice version of the {MakeOctahedron} procedure.\n" \
  "  19-09-2000 : [???] Added the procedure {MakeDodecahedronTriang}.\n" \
  "  27-10-2000 : [???] Modified {SetCubeProperties} procedure.\n" \
  "  26-01-2007 : [???] Brought {MakeBigRawCube} from {MakeBigCubeFixed.c}.\n" \
  "  26-01-2007 : [???] Brought {Make1DCubeArrayBis} from {MakeObjectTriang.c}.\n" \
  "  02-02-2007 : [???] Brought {EmphasizeBigRawCubeOriginalElems} from {MakeBigCubeFixed.c}.\n" \
  "  02-02-2007 : [???] Brought {FixBigRawCubeCoords} from {MakeBigCubeFixed.c}."