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

/* 3D map operations pecialized for triangulations -- maps with tetrahedral cells.*/
/* Last edited on 2023-02-12 07:47:14 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 <jsfile.h>
#include <jsmath.h>

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

FourNodes_t TetraNegNodes(Place_t p)
  { /* Valid for the versions Spin(p), Clock(p) and SpinClock(p).  */
    assert(PnegP(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(PnegP(p) == PposP(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(PposP(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(PposP(p) == PnegP(PrevE(NextF(p))));
    return (FourNodes_t){{n0,n1,n2,n3}};
  }

FourWalls_t TetraWalls(Place_t p)
  { assert(PnegP(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(PnegP(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 ((PnegP(b) == NULL)){ return TRUE; }
      b = NextF(b); 
    } while (b != p);
    return FALSE;
  }

bool_t WallIsBorder(Place_t p)
  { return ((PposP(p) == NULL) || (PnegP(p) == NULL)); }

FiveNodes_t TetraNegPosNodes(Place_t p)
  { /* Valid for the versions Spin(p), Clock(p) and SpinClock(p). */ 
    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(PposP(p) == PnegP(PrevE(NextF(p))) &&  (PnegP(p) == PposP(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 = WedgeBase(MakeWedge());
        Place_t b = WedgeBase(MakeWedge());
        Place_t c = WedgeBase(MakeWedge());
        Wall_t f = PWall(a);
        Node_t u = MakeNode(0);
        Node_t v = MakeNode(1);
        Node_t w = MakeNode(2);
        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(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];
  }

void InitCoords(Random_t *coins, Coords_t *c, double r)
  { int i;
    for (i = 0; i < c->ne; i++)
      { c->e[i] = (r4_t)
          { { RandomDouble(coins, -r, r),
              RandomDouble(coins, -r, r),
              RandomDouble(coins, -r, r),
              RandomDouble(coins, -r, r)
          } };
      }
  }

Coords_t GenCoords(Random_t *coins, ElemTableRec_t *top)
  { Coords_t c = r4_vec_new(top->node.ne);
    InitCoords(coins, &c, 1.0);
    return c;
  }

#define ST_FILE_TYPE "state"
#define ST_FILE_VERSION "99-08-25"

void WriteCoords(char *name, char *tag, int seq, Coords_t *c, char *cmt /* DF "" */)
    if (seq >= 0)
      { char *fileName = jsprintf("%s%s-%06d.st", name, tag, seq); }
    else
      { char *fileName = jsprintf("%s%s.st", name, tag); }
    FILE *wr = open_write(fileName, TRUE);
    FWriteCoords(wr, c, cmt);
    fclose(wr);
  }

Coords_t ReadCoords(char *name, char *tag, int seq, uint NV)
    if (seq >= 0)
      { char *fileName = jsprintf("%s%s-%06d.st", name, tag, seq); }
    else
      { char *fileName = jsprintf("%s%s.st", name, tag); }
    FILE *rd = open_read(fileName, TRUE);
    char *cmt;
    Coords_t c = r4_vec_new(NV);
    FReadCoords(rd, &c, &cmt);
    free(cmt);
    fclose(rd);
    return c;
  }

void FWriteCoords(FILE *wr, Coords_t *c, char *cmt)
  { uint NV = c->ne;
    filefmt_write_header(wr, ST_FILE_TYPE, ST_FILE_VERSION);
    fprintf(wr, "vertices %d\n", NV);
    filefmt_write_comment(wr, cmt, '|');

    filefmt_write_comment(wr, "\nNode data:\n", '|');

    uint vWidth = (NV <= 1 ? 1 : digits(NV-1)); /* Index width for alignment. */
    int i;
    for (i = 0; i < NV; i++)
      { fprintf(wr, "%*d ", vWidth, i);
        r4_gen_print(wr, &(c->e[i]), "%12.4e", "", " ", "");
        fprintf(wr, "\n");
      }
      filefmt_write_footer(wr, ST_FILE_TYPE);
      fflush(wr);
  }

void FReadCoords(FILE *rd, Coords_t *c, char **cmtP)
  { 
    /* Read the file header: */
    filefmt_read_header(rd, ST_FILE_TYPE, ST_FILE_VERSION);

    /* Read the vertex count: */
    fget_skip_spaces(rd);
    fget_match(rd, "vertices");
    uint NV = fget_uint32(rd, 10);
    fget_eol(rd);
    
    /* Check whether {c} has the correct length: */
    demand(NV == c->ne, "state vector in file has wrong size");
    
    /* Read the client comments: */
    char *cmt = filefmt_read_comment(rd, '|');
    if (cmtP != NULL) { (*cmtP) = cmt; } else { free(cmt); }

    /* Read the vertex coordinates {c[0..NV-1]}: */
    int i, j;
    for (i = 0; i < NV; i++)
      { uint vn = fget_uint32(rd, 10);
        demand(vn == i, "vertex number mismatch");
        r4_t *ci = &(c->e[i]);
        for (j = 0; j < 4; j++)
          { ci->c[j] = fget_double(rd); }
          
        fget_eol(rd);
      }

      filefmt_read_footer(rd, ST_FILE_TYPE);
  }

#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."