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

/* 3D map operations pecialized for triangulations -- maps with tetrahedral cells.*/
/* Last edited on 2007-01-31 14:12:53 by stolfi */

#define Triangulation_C_COPYRIGHT \
  "Copyright © 2000 by the State University of Campinas (UNICAMP)"

#include <Octf.h>
#include <OctfEnum.h>
#include <Mis.h>
#include <Random.h>

#include <r3.h>
#include <r4.h>
#include <fget.h>
#include <filefmt.h>
#include <jsmath.h>

#define _GNU_SOURCE
#include <stdio.h>
#include <string.h>
#include <assert.h>

Node_t MakeNode(void)
  { return (Node_t)notnull(malloc(sizeof(NodeRec_t)), "no mem"); }
  
Cell_t MakeCell(void)
  { return (Cell_t)notnull(malloc(sizeof(CellRec_t)), "no mem"); }

NodeOrCell_t Org(Place_t p)
  { return PWedge(p)->org[SrotBits(p)]; }

NodeOrCell_t Pneg(Place_t p)
  { p = Srot(p);
    return PWedge(p)->org[SrotBits(p)];
  }

NodeOrCell_t Ppos(Place_t p)
  { p = Tors(p);
    return PWedge(p)->org[SrotBits(p)];
  }

void SetOrg(Place_t p, Node_t n)
  { assert(DualBit(p) == 0); /* [!!! Added by JS. If {p} is dual, should use SetPneg instead. !!!] */
    PWedge(p)->org[SrotBits(p)] = (NodeOrCell_t)n;
  }

void SetPneg(Place_t p, Cell_t n)
  { assert(DualBit(p) == 1); /* [!!! Added by JS. If {p} is dual, should use SetOrg instead. !!!] */
    p = Srot(p);
    PWedge(p)->org[SrotBits(p)] = (NodeOrCell_t)n;
  }

void SetPpos(Place_t p, Cell_t n)
  { assert(DualBit(p) == 1); /* [!!! Added by JS. If {p} is dual, should use SetDes instead. !!!] */
    p = Tors(p);
    PWedge(p)->org[SrotBits(p)] = (NodeOrCell_t)n;
  }

void SetOrgAll(Place_t p, Node_t n)
  { assert(DualBit(p) == 0); /* [!!! Added by JS. If {p} is dual, should use SetPneg instead. !!!] */
    Place_vec_t pv = (Place_vec_t){/*nel*/ 1, /*el*/ &p};
    Place_vec_t nei = RenumberEdgesOutOfNodes(&pv);
    int i;
    for(i = 0; i < nei.nel; i++)
      { Place_t b = nei.el[i];
        Place_t c = b;  
        do { SetOrg(c,n); c = NextF(c); } while (c != b);
      }
  }

void SetPnegOfWall(Place_t p, Cell_t n)
  { assert(DualBit(p) == 1); /* [!!! Added by JS. If {p} is dual, should use SetOrg instead. !!!] */
    Place_t t = p;
    do { SetPneg(t, n); t = PrevE(t); } while (t != p);
  }

void SetPnegOfNearbyWalls(Place_t p, Cell_t n)
  { assert(DualBit(p) == 1); /* [!!! Added by JS. If {p} is dual, should use SetOrg instead. !!!] */
    SetPnegOfWall(p,n);
    Place_t t = p;
    do {
      SetPnegOfWall(Clock(PrevF(t)),n);
      t = PrevE(t); 
    } while (t != p); 
  }

void SetPposOfWall(Place_t p, Cell_t n)
  { assert(DualBit(p) == 1); /* [!!! Added by JS. If {p} is dual, should use SetDes instead. !!!] */
    Place_t t = p;
    do { SetPpos(t, n);  t = PrevE(t); } while (t != p);
  }

void SetPposQUX(Place_t p, Cell_t n)
  { assert(DualBit(p) == 1); /* [!!! Added by JS. If {p} is dual, should use SetDes instead. !!!] */
    Place_t t = Clock(p);
    SetPnegOfNearbyWalls(t,n);
    do 
      { SetPnegOfNearbyWalls(PrevF(t),n);
        t = PrevE(t); 
      } while (t != Clock(p)); 
  }

Node_t OrgV(Place_t p)
  { return (Node_t)(Org(p)); }

Node_t DesV(Place_t p)
  { return OrgV(Clock(p)); }

Cell_t PnegP(Place_t p)
  { return (Cell_t)(Pneg(p)); }

Cell_t PposP(Place_t p)
  { return PnegP(Clock(p)); }

FourNodes_t TetraNegNodes(Place_t p)
  { /* Valid for the versions Spin(p), Clock(p) and SpinClock(p).  */
    assert(Pneg(p) != NULL);
    Node_t n0 = OrgV(p);
    Node_t n1 = OrgV(NextE(p));
    Node_t n2 = OrgV(PrevE(p));
    Node_t n3 = OrgV(PrevE(PrevF(p)));
    assert(Pneg(p) == Ppos(PrevE(PrevF(p))));
    assert((n0 != n1) && (n0 != n2) && (n0 != n3) && (n1 != n2) && (n1 != n3) && (n2 != n3));
    return (FourNodes_t){{n0,n1,n2,n3}};
  }
 
FourNodes_t TetraPosNodes(Place_t p)
  { /* Valid for the versions Spin(p), Clock(p) and SpinClock(p). */
    /* [!!! Should be just TetraNegNodes(q) where {q = Spin(p)} or {q = Clock(p)} ? !!!] */
    assert(Ppos(p)!=NULL);
    Node_t n0 = OrgV(p);
    Node_t n1 = OrgV(NextE(p));
    Node_t n2 = OrgV(PrevE(p));
    Node_t n3 = OrgV(PrevE(NextF(p)));
    assert(Ppos(p) == Pneg(PrevE(NextF(p))));
    return (FourNodes_t){{n0,n1,n2,n3}};
  }

FourWalls_t TetraWalls(Place_t p)
  { assert(Pneg(p)!=NULL);
    Wall_t f0 = PWall(p);
    Wall_t f1 = PWall(PrevF(p));
    Wall_t f2 = PWall(PrevF(NextE(p)));
    Wall_t f3 = PWall(PrevF(PrevE(p)));
    assert((f0 != f1) && (f0 != f2) && (f0 != f3) && (f1 != f2) && (f1 != f3) && (f2 != f3)); 
    return (FourWalls_t){{f0,f1,f2,f3}};
  }

SixEdges_t TetraEdges(Place_t p)
  { assert(Pneg(p)!=NULL);
    Edge_t e0 = PEdge(p);
    Edge_t e1 = PEdge(NextE(p));
    Edge_t e2 = PEdge(PrevE(p));
    Edge_t e3 = PEdge(NextE(PrevF(p)));
    Edge_t e4 = PEdge(PrevE(PrevF(p)));
    Edge_t e5 = PEdge(NextE(PrevF(NextE(p))));
    assert((e0 != e1) && (e0 != e2) && (e0 != e3) && (e0 != e4) && (e0 != e5));
    assert((e1 != e2) && (e1 != e3) && (e1 != e4) && (e1 != e5));
    assert((e2 != e3) && (e2 != e4) && (e2 != e5));
    assert((e3 != e4) && (e3 != e5));
    assert(e4 != e5); 
    return (SixEdges_t){{e0,e1,e2,e3,e4,e5}};
  }

ThreeEdges_t TriangEdges(Place_t p)
  { Edge_t e0 = PEdge(p);
    Edge_t e1 = PEdge(NextE(p));
    Edge_t e2 = PEdge(PrevE(p));
    assert((e0 != e1) && (e0 != e2) && (e1 != e2)); 
    return (ThreeEdges_t){{e0,e1,e2}};
  }

bool_t EdgeIsBorder(Place_t p)
  { Place_t b = p;
    do { 
      if ((Pneg(b) == NULL)){ return TRUE; }
      b = NextF(b); 
    } while (b != p);
    return FALSE;
  }

bool_t WallIsBorder(Place_t p)
  { return ((Ppos(p) == NULL) || (Pneg(p) == NULL)); }

FiveNodes_t TetraNegPosNodes(Place_t p)
  { /* Valid for the versions Spin(p), Clock(p) and SpinClock(p). */ 
    /* Used {Org} instead of {OrgV} so it could be used in the 
      dual map. However the assertion is valid only 
      if {p} is on the primal map. */
    Node_t n0 = OrgV(p);
    Node_t n1 = OrgV(NextE(p));
    Node_t n2 = OrgV(PrevE(p));
    Node_t n3 = OrgV(PrevE(PrevF(p)));
    Node_t n4 = OrgV(PrevE(NextF(p)));
    assert(Ppos(p) == Pneg(PrevE(NextF(p))) &&  (Pneg(p) == Ppos(PrevE(PrevF(p)))));
    return (FiveNodes_t){{n0,n1,n2,n3,n4}};
  }

EightPlaces_t MakeTetraTopo(uint nx, uint ny)
  { uint wedgeCount = 0; 
    uint cellCount = 0;

    auto Place_t MakeTriangle(void);
    /* Creates a new `naked' triangular wall. The resulting structure
      has three wedges connected by {NextE} links, sharing the same
      {Wall_t}, without any {Cell_t}s. */

    Place_t MakeTriangle(void)
      { Place_t a = MakeWedge();
        Place_t b = MakeWedge();
        Place_t c = MakeWedge();
        Wall_t f = PWall(a);
        Node_t u = MakeNode();
        Node_t v = MakeNode();
        Node_t w = MakeNode();
        PWedge(a)->num = wedgeCount; wedgeCount++;         
        SetOrg(a, u);
        SetOrg(Clock(a),v);

        PWedge(b)->num = wedgeCount; wedgeCount++;
        SetNextE(a,b); SetWallInfo(b,f);
        SetOrg(b,v); SetOrg(Clock(b),w);

        PWedge(c)->num = wedgeCount; wedgeCount++;
        SetNextE(b,c); SetWallInfo(c,f);
        SetOrg(c, w);  SetOrg(Clock(c), OrgV(a));
        return a;
      }

    auto Place_t AddTwoTriangles(Place_t a);
      /* Builds an incomplete tetrahedral cell by adding two
        additional triangular walls to an exisiting `naked' wall. The
        argument {a} is a place on the existing wall, as returned by
        MakeTriangle(). The procedure returns a place {f} on one of
        the new triangles. */
        
    Place_t AddTwoTriangles(Place_t a)
      { Place_t c = NextE(NextF(a));
        SetEdgeInfo(a, PEdge(PrevE(c)));
        Place_t f = MakeTriangle();
        Place_t g = MakeTriangle();

        SetNextF(PrevE(a), PrevE(f));
        SetEdgeInfo(PrevE(f), PEdge(PrevE(a)));
        SetOrgAll(PrevE(a), OrgV(PrevE(a)));
        SetOrgAll(Clock(PrevE(a)), OrgV(Clock(PrevE(a))));

        SetNextF(Clock(f), NextE(c));
        SetEdgeInfo(Clock(f), PEdge(NextE(c))); 
        SetOrgAll(NextE(c), OrgV(NextE(c)));
        SetOrgAll(Clock(NextE(c)), OrgV(Clock(NextE(c))));     

        SetNextF(NextE(g), NextE(f));
        SetEdgeInfo(NextE(g), PEdge(NextE(f)));
        SetOrgAll(NextE(f), OrgV(NextE(f)));
        SetOrgAll(Clock(NextE(f)), OrgV(Clock(NextE(f))));

        SetNextF(g, c);
        SetEdgeInfo(g, PEdge(c));
        SetOrgAll(c, OrgV(c));
        SetOrgAll(Clock(c), OrgV(Clock(c))); 

        SetNextF(PrevE(g), Clock(NextE(a)));
        SetEdgeInfo(PrevE(g), PEdge(Clock(NextE(a)))); 
        SetOrgAll(NextE(a), OrgV(NextE(a)));
        SetOrgAll(Clock(NextE(a)), OrgV(Clock(NextE(a))));

        Cell_t p = MakeCell();
        p->num = cellCount; cellCount++;
        SetPnegOfNearbyWalls(f,p);
        return f;
      }

    auto Place_t FinishCellRow(Place_t a);
      /* Adds one new row of tetradedral cells to a started row, and
        returns a place on the last cell. */
        
    Place_t FinishCellRow(Place_t a)
      { int col;
        for (col = 0; col < nx; col++) { a = AddTwoTriangles(a); }
        return a;
      }

    Place_t t = NULL, b = NULL;
    EightPlaces_t ca = (EightPlaces_t){{NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL}};

    /* Create bottom row of triangles: */
    int col;
    for (col = 0; col < nx; col++)
      { Place_t c = MakeTriangle();
        if (col == 0)
          { t = c; }
        else
          { SpliceWalls(NextE(b), Clock(PrevE(c)));
            SetEdgeInfo(NextE(b), PEdge(Clock(PrevE(c))));
            SetOrgAll(NextE(b), OrgV(NextE(b))); 
            SetOrgAll(Clock(NextE(b)), OrgV(Clock(NextE(b))));
            b = c;
          }
      }
    ca.p[2] = t;
    ca.p[7] = Spin(Clock(b));

    /* Add {ny} rows of cells: */
    int row;
    for (row = 0; row < ny; row++)
      { /* First triangle of row: */
        Place_t a = MakeTriangle();
        if (row == 0){ ca.p[1] = Clock(PrevE(a)); }           
        if (row == ny-1){ ca.p[4] = Spin(PrevE(a)); }

        SpliceWalls(Clock(a), Clock(PrevE(t)));
        SetOrgAll(PrevE(t), OrgV(PrevE(t)));
        SetOrgAll(Clock(PrevE(t)), OrgV(Clock(PrevE(t))));

        /* Add {nx} cells to row: */
        Place_t aa = FinishCellRow(a);
        
        t = PrevF(t);
        SetOrg(PrevE(t), OrgV(PrevE(a)));
        b = PrevF(b);
        SetOrg(Clock(NextE(b)), OrgV(PrevE(aa)));
        if (row == 0){ ca.p[0] = Clock(Spin(PrevE(aa))); } 
        if ((row == ny-1)){ ca.p[5] = PrevE(aa); }
      }

    ca.p[3] = Spin(t);
    ca.p[6] = Clock(b);
    return ca;
  }       

void EmphasizeTetrahedron(Place_t a, Place_t b, uint n)
  {
    auto void HideNode(Node_t v);
      /* Hide the node {v}. */ 

    void HideNode(Node_t v)
      { v->exists = FALSE; }

    auto void HideEdge(Place_t a);
      /* Hide the edge {PEdge(a)}.*/ 
     
    void HideEdge(Place_t a)
      { PEdge(a)->exists = FALSE; }        

    auto void HideWall(Place_t a);
      /* Hide the wall {PWall(a)}.*/ 
    
    void HideWall(Place_t a)
      { PWall(a)->exists = FALSE; }        

    auto void HideWallsOfEdge(Place_t a);
      /* Hide all the walls incident to the edge {PEdge(a)}.
        [!!! Actually hides the {n-1} walls that come before {PWall(a)}
        around {Pedge(a)}. !!!] */
    
    void HideWallsOfEdge(Place_t a)
      { Place_t an = PrevF(a);
        int j;
        for (j = 1; j < n; j++) { HideWall(an); an = PrevF(an); }
      }

    /* Collects certain {n} places {ta[0..n-1],tb[0..n-1]}: [???] */
    Place_t ta[101], tb[101];
    ta[0] = a; tb[0] = b;
    int i;
    for (i = 1; i < n; i++)
      { ta[i] = Clock(PrevE(NextF(NextE(ta[i-1]))));
        tb[i] = Clock(PrevE(NextF(NextE(tb[i-1]))));
      }
   
    /* Hide {n-1} walls incident to the edges of those places: [???] */
    for (i = 0; i < n; i++)
      { HideWallsOfEdge(ta[i]); HideWallsOfEdge(tb[i]); }

    /* Hide the origin nodes of those places: [???] */
    for (i = 1; i < n; i++)
      { HideNode(OrgV(tb[i])); 
        HideNode(OrgV(ta[i]));
        Place_t dn = tb[i];
        int j;
        for (j = 0; j <= n /* [!!! Should it be {j < n}? !!!] */; j++)
          { HideEdge(PrevE(dn)); dn = PrevF(dn); }
      }

    /* Hide certain edges related to those places: [???] */
    for (i = 0; i < (n-1); i++)
      { HideEdge(NextE(ta[i]));
        HideEdge(NextE(NextF(ta[i])));
      }
  }

Place_t Glue(Place_t a, Place_t b, uint n, bool_t setorg)
  { /* fprintf(stderr, "colou\n"); */
    assert(n >= 1); /* sanity check */
    
    /* Collect interesting places {ta[0..n-1], tb[0..n-1]}: */
    Place_t ta[101], tb[101];
    int i;
    ta[0] = a; tb[0] = b;
    if (n > 1)
      { for (i = 1; i < n; i++)
          { ta[i] = Clock(PrevE(PrevF(NextE(ta[i-1])))); assert(ta[i]!=a);
            tb[i] = Clock(PrevE(NextF(NextE(tb[i-1])))); assert(tb[i]!=b);  
          }
      }
      
    Meld(b, a);
    
    /* Update the edge slots for {i == 0}: */
    SetRingEdgeInfo(a,        PEdge(a));
    SetRingEdgeInfo(NextE(a), PEdge(NextE(a)));
    SetRingEdgeInfo(PrevE(a), PEdge(PrevE(a)));
    
    if (setorg)
      { /* Update the node slots for {i == 0}: */
        SetOrgAll(a,        OrgV(a));
        SetOrgAll(Clock(a), OrgV(Clock(a)));
        SetOrgAll(NextE(a), OrgV(NextE(a)));
        SetOrgAll(Clock(NextE(a)), OrgV(Clock(NextE(a))));
        SetOrgAll(PrevE(a), OrgV(PrevE(a)));
        SetOrgAll(Clock(PrevE(a)), OrgV(Clock(PrevE(a))));
      }
    
    /* Update the cell slots for {i == 0}: */
    SetPneg(a, PnegP(b));
    SetPneg(PrevE(a), PnegP(PrevE(b)));
    SetPneg(NextE(a), PnegP(NextE(b)));
    for (i = 1; i < n; i++)
      { Meld(tb[i],ta[i]);

        /* Update the edge slots for {i > 0}: */
        SetRingEdgeInfo(ta[i], PEdge(ta[i]));
        SetRingEdgeInfo(NextE(ta[i]), PEdge(NextE(ta[i])));
        SetRingEdgeInfo(PrevE(ta[i]), PEdge(PrevE(ta[i])));

        if (setorg)
          { /* Update the node slots for {i > 0}: */
            SetOrgAll(ta[i], OrgV(ta[i]));
            SetOrgAll(Clock(ta[i]), OrgV(Clock(ta[i])));
            SetOrgAll(NextE(ta[i]), OrgV(NextE(ta[i])));
            SetOrgAll(Clock(NextE(ta[i])), OrgV(Clock(NextE(ta[i]))));
            SetOrgAll(PrevE(ta[i]), OrgV(PrevE(ta[i])));
            SetOrgAll(Clock(PrevE(ta[i])), OrgV(Clock(PrevE(ta[i]))));
          }

        /* Update the cell slots for {i > 0}: */
        Cell_t f = PnegP(tb[i]);
        Cell_t g = PnegP(PrevE(tb[i]));
        Cell_t h = PnegP(NextE(tb[i]));
        SetPneg(ta[i], f);
        SetPneg(PrevE(ta[i]), g);
        SetPneg(NextE(ta[i]), h);
      }

    if (setorg)
      { SetOrgAll(Clock(PrevE(PrevF(PrevE(ta[n-1])))),OrgV(PrevE(a))); }
    return ta[n-1];
  }

#define Triangulation_C_author \
  "Created 2000 by L. A. P. Lozada, based on 1996 procedures by\n" \
  "  R. M. Rosi and J. Stolfi.\n" \
  "Revisions:\n" \
  "  03-08-2000 : Optimized version of the Read and Write procedures by J. Stolfi.\n" \
  "  30-08-2000 : Modified the Read and Write procedures for include the case\n" \
  "               when the cells are octahedra.\n" \
  "  07-10-2000 : Added procedure for compute the barycenter of a tetrahedron.\n" \
  "  27-01-2001 : Modified for exploding cubic cells."