/* Code for building and using the cycle of faces around an edge.  */
/* Last edited on 2015-10-09 15:34:50 by stolfilocal */

typedef struct stmesh_edge_rep_t
  { stmesh_vert_t endv[2];  /* Endpoints in vertex table, in {0..mesh.nv-1}. */
    stmesh_face_t firstf;   /* First oriented face incident to edge. */
  } stmesh_edge_rep_t;
  /* An UNORIENTED edge of the mesh is defined by two distinct vertices (endpoints),
    {.endv[0..1]}.  In an UNORIENTED edge, 
    the vertices are sorted so that their indices in {mesh.v} are in increasing order.  
    
    A pointer {e} to an ORIENTED edge (a {stmesh_edge_t}) 
    is {ue|d} where {ue} is a pointer to the the corresponding {stmesh_edge_rep_t}
    record, and {d} is the /direction bit/ of {e}. In the /natural/ 
    orientation (direction bit {d = 0}), the edge is traversed from {ue.endv[0]}
    (origin) to {ue.endv[1]} (destination). In the /reversed/ orientation,
    ({d = 1}) it is traversed the other way: the origin of {e} is
    {ue.endv[1]}, and the destination of {e} is {ue.endv[0]}.
    
    The field {ue.firstf} is a pointer to an oriented face incident to
    the edge {ue}, whose orientation is such that its base edge is {ue} in
    its natural orientation. */

typedef struct stmesh_face_rep_t
  { stmesh_edge_t side[3];   /* Oriented edges along the perimeter. */
    stmesh_face_t nextf[3];  /* Next oriented face incident to each side. */
    stmesh_face_t prevf[3];  /* Prev oriented face incident to each side. */
    int32_t minZ;            /* Minimum {Z} coordinate of any vertex, as multiple of {eps}. */
    int32_t maxZ;            /* Maximum {Z} coordinate of any vertex, as multiple of {eps} . */
  } stmesh_face_rep_t;
  /* An UNORIENTED face of the mesh is a triangle defined by its three distinct edges 
    and three distinct vertices.  
    
    The fields {.side[0..2]} are the ORIENTED edges along the 
    perimeter of the triangle, in consistent order.  Namely, the destination
    of the oriented edge {.side[k]} is the origin of the oriented edge
    {.side[(k+1)%3]}, for {k} in {0..2}.
    
    In an UNORIENTED face, the cyclic order of the sides and the
    starting side are chosen so that the ORIENTED index of {side[0]} is
    minimized. Since all six oriented edges on the perimeter are
    distinct, that orientation is unique.
    
    A pointer {f} to an ORIENTED face is {uf | 2*k + d},
    where {uf} is a pointer to the the corresponding {stmesh_face_rep_t}
    record, {k} is a /side selector/ in {0..2}, and {d} is a /direction bit/
    that specifies a sense of traversal of the perimeter. 
    In the /natural direction/ ({d = 0}), the base edge of {f} is {uf.side[k]},
    and the sides are traversed in increasing order, starting from that side -- namely,
    {uf.side[k]}, {uf.side[(k+1)%3]}, and {uf.side[(k+2)%3]}. In the /reversed
    orientation/ ({d = 1}), the sides are traversed in the opposite order
    and orientation, also starting from the base edge:
    {EREV(uf.side[k])}, {EREV(uf.side[(k+2)%3])}, {EREV(uf.side[(k+1)%3])},
    where {EREV(x)} is {stmesh_edge_reverse(x)}.
    
    All the faces that share an edge are collected in a circular list,
    defined by the {.nextf} and {.prevf} pointers in the face records.
    Specifically, {uf.nextf[k]} is the next oriented face
    {g} whose base edge is {uf.side[k]}. The word 'next' here refers to the sense of rotation
    around the oriented edge {uf.side[k]} given by the right-hand rule.
    The field {uf.prevf[k]} is similar, but going in the opposite sense
    around that edge. Following those pointers, taking orientations into
    account, gives all faces incident to that edge. */


stmesh_face_t stmesh_get_edge_first_face(stmesh_t mesh, stmesh_edge_t e);
  /* One of the oriented faces whose base edge is {e}. */

uint32_t stmesh_edge_degree(stmesh_t mesh, stmesh_edge_t e);
  /* Returns the degree of the oriented edge with index {e}, that is, the
    number of faces that share it. */

stmesh_face_t stmesh_next_face(stmesh_t mesh, stmesh_face_t f, int k);
  /* Returns an oriented face of the mesh that has the the same base edge 
    as {ixf}, including orientation.  The face will be {k} positions ahead of {f}
    in the circular ordering of faces around that oriented edge, according to the 
    right-hand rule.  The count {k} may be negative, to move in the reverse direction.
    A count of zero will return {f} itself.  In particular, if the base edge of {f} 
    is a free border edge, the result will be {f} itself, for any {k}. */

void stmesh_fill_edge(stmesh_t mesh, stmesh_edge_rep_t *e, stmesh_vert_unx_t uxv[]);
  /* Fills the edge record {*e}, given the indices {uxv[0..1]} of its
    endpoint vertices.  The indices must be in increasing order. 
    Assumes that {mesh->v} is set. The field {e.firstf} is set to {NULL}. */
        
void stmesh_fill_face(stmesh_t mesh, stmesh_face_rep_t *f, stmesh_edge_unx_t uxe[]);
  /* Fills the fields of face record {*f} given the indices {uxe[0..2]} of the unoriented
    edges that bound the face.  The indices should be sorted in increasing order.
    Assumes that {mesh.{nv,v,ne,e}} are set. The pointers {f.nextf} and {f.prevf} are set to {NULL}. */

void stmesh_set_links(stmesh_t mesh);
  /* Sets the {ue.firstf} link in every unoriented edge {ue} in {mesh},
    and the links {uf.nextf[0..2],uf.nextf[0..2]} links in every 
    unoriented face {uf}.  Prints warnings if there are coplanar
    or zero-area triangles.  */
            
void stmesh_insert_face_in_edge_cycle(stmesh_t mesh, stmesh_face_ix_t ixf);
  /* Inserts the oriented face with index {ixf} in the circular list of faces incident to 
    the its base edge.  Specifically, sets the pointers {.ixnextf[k]} and {.ixprevf[k]}
    of that face, where {k} is the position of the base edge or its reverse in {f.ixside[0..2]}; 
    and also sets the necessary pointers in other faces. */
    
void stmesh_insert_face_in_edge_cycle(stmesh_t mesh, stmesh_edge_ix_t ixe, stmesh_face_ix_t ixf)
  { 
    /* Get the index of the oriented base edge face: */
    stmesh_edge_ix_t stmesh_get_face_base_edge(mesh, ixf);
    /* Get the first oriented face that has {ixe} as its base edge: */
    stmesh_face_ix_t jxf = stmesh_get_first_face(mesh, ixe);
    if (jxf == IX_NULL)
      { /* First face with this edge: */
        stmesh_make_singleton_face_cycle(mesh, ixf);
      }
    else
      { assert(stmesh_get_face_base_edge(mesh, jxf) == ixe);
        /* Advance {jxf} around {ixe} until finding the proper place for {ixf}: */
        stmesh_face_ix_t axf, bxf;
        stmesh_find_face_around_edge(mesh, jxf, ixf, &axf, &bxf);
        stmesh_insert_face_between(mesh, ixf, axf, bxf);
      }
    /* Remember this face as one of the incident faces: */
    stmesh_set_first_face(mesh, ixe, ixf);
  }

stmesh_face_ix_t stmesh_get_first_face(stmesh_t mesh, stmesh_edge_ix_t ixe);
  /* Returns the index {ixf} of some oriented face whose base edge is the oriented edge whose eindex is {ixe}.
    may return {IX_NULL} while the mesh is being built, before all links are properly set.
    Never returns {IX_NULL} afterwards. */ 

stmesh_face_ix_t stmesh_get_first_face(stmesh_t mesh, stmesh_edge_ix_t ixe)
  {
    stmesh_edge_t *e = stmesh_get_edge_address(mesh, EDGE_UNX(ixe));
    stmesh_face_ix_t ixf = e->ixfirstf;
    if (ixf != IX_NULL)
      { int d = EDGE_D(ixe);
        if (d != 0) { ixf = FFLIP(ixf); }
        assert(stmesh_get_face_base_edge(mesh, ixf) == ixe);
      }
    return ixf;
  }

void stmesh_make_singleton_face_cycle(stmesh_t mesh, stmesh_face_ix_t ixf);
  /* Receives the index {ixf} of an oriented face. Sets the proper
    pointers so that the loop of faces around its base edge consists of
    ixf} alone. The affected pointers must be {IX_NULL} */

void stmesh_make_singleton_face_cycle(stmesh_t mesh, stmesh_face_ix_t ixf)
  {
    stmesh_face_t *f = stmesh_get_face_address(mesh, FACE_UNX(ixf))
    int k = FACE_K(ixf); 
    int d = FACE_D(ixf);

    assert(f->ixnextf[k] == IX_NULL);
    assert(f->ixprevf[k] == IX_NULL);
    
    if (d == 0)
      { f->ixnextf[k] = ixf;
        f->ixprevf[k] = ixf;
      }
    else
      { f->ixprevf[k] = FFLIP(ixf);
        f->ixnextf[k] = FFLIP(ixf);
      }

    assert(stmesh_get_next_face(mesh, ixf) == ixf);
    assert(stmesh_get_prev_face(mesh, ixf) == ixf);
  }

void stmesh_insert_face_between
  ( stmesh_t mesh, 
    stmesh_face_ix_t ixf, 
    stmesh_face_ix_t axf, 
    stmesh_face_ix_t bxf
  );
  /* Receives the indices {ixf,axf,bxf} of three oriented faces, which should have the same base edge.
    Sets the proper pointers so that the face {ixf} gets inserted in the 
    loop of faces around that oriented edge, just after {axf} and just before {bxf}.
    
    The faces {axf} and {bxf} must be conscutive faces in that loop, in that order; they may be the same oriented face,
    if the loop is a singleton, otherwise theu must have distinct unoriented faces.
    The new face {ixf} must not be in the loop, and its 
    unoriented face must be distinct from those of {axf} and {bxf}. The relevant pointers 
    of {ixf} must be {IX_NULL}. */

void stmesh_insert_face_between
  ( stmesh_t mesh, 
    stmesh_face_ix_t ixf, 
    stmesh_face_ix_t axf, 
    stmesh_face_ix_t bxf
  )
  {
    stmesh_edge_ix_t ixe = stmesh_get_face_base_edge(mesh, ixf);
    stmesh_edge_ix_t axe = stmesh_get_face_base_edge(mesh, axf);
    stmesh_edge_ix_t bxe = stmesh_get_face_base_edge(mesh, bxf);
    demand((ixe == axe) && (ixe == bxe), "the three faces must have the same base edge");
    
    demand(stmesh_get_next_face(mesh, axf) == bxf, "faces {axf,bxf} must be consecutive");
    assert(stmesh_get_prev_face(mesh, bxf) == axf);
    
    stmesh_face_t *fi = stmesh_get_face_address(mesh, FACE_UNX(ixf))
    int ki = FACE_K(ixf); 
    int di = FACE_D(ixf);

    stmesh_face_t *fa = stmesh_get_face_address(mesh, FACE_UNX(axf))
    int ka = FACE_K(axf); 
    int da = FACE_D(axf);
    
    stmesh_face_t *fb = stmesh_get_face_address(mesh, FACE_UNX(ixf))
    int kb = FACE_K(ixf); 
    int db = FACE_D(ixf);

    demand(fi->ixnextf[ki] == IX_NULL, "edge {ixf} is already attached");
    demand(fi->ixprevf[ki] == IX_NULL, "edge {ixf} is already attached");
    
    if (di == 0)
      { fi->ixnextf[ki] = bxf;
        fi->ixprevf[ki] = axf
      }
    else
      { fi->ixprevf[ki] = FFLIP(bxf);
        fi->ixnextf[ki] = FFLIP(axf);
      }

    if (da == 0)
      { fa->ixnextf[ka] = ixf; }
    else
      { fa->ixprevf[ka] = FFLIP(ixf); }

    fi (db == 0)
      { fb->ixprevf[kb] = ixf; }
    else
      { fb->ixnextf[kb] = FFLIP(ixf); }

    assert(stmesh_get_next_face(mesh, axf) == ixf);
    assert(stmesh_get_prev_face(mesh, ixf) == axf);
    assert(stmesh_get_prev_face(mesh, bxf) == ixf);
    assert(stmesh_get_next_face(mesh, ixf) == bxf);
  }
    
void stmesh_find_face_around_edge
  ( stmesh_t mesh, 
    stmesh_face_ix_t sxf, 
    stmesh_face_ix_t ixf,
    stmesh_face_ix_t *axfP,
    stmesh_face_ix_t *bxfP
  );
  /* Finds the place where the oriented face with index {ixf} should be in the loop of 
    faces around its base edge.  The face {ixf} must not be part of that loop,
    whereas {sxf} should be any face of that loop. Returns in {*axfP} and {*bxfP}
    the faces of the loop such that {ixf} should go just after the first and just before the 
    second. If the loop is a singleton, consisting of just the face {sxf}, returns
    {*axfP==bxfP==sxf}.  */
    
void stmesh_find_face_around_edge
  ( stmesh_t mesh, 
    stmesh_face_ix_t sxf, 
    stmesh_face_ix_t ixf,
    stmesh_face_ix_t *axfP,
    stmesh_face_ix_t *bxfP
  )
  {
    stmesh_edge_ix_t sxe = stmesh_get_face_base_edge(mesh, sxf);
    stmesh_edge_ix_t ixe = stmesh_get_face_base_edge(mesh, ixf);
    demand(sxe == ixe, "base edges are not the same");

    /* Get the first two faces around the edge: */
    stmesh_face_ix_t axf = sxf;
    stmesh_face_ix_t bxf = stmesh_get_next_face(mesh, axf);
    if (axf != bxf)
      { /* There are at least two faces, must actually search: */
        while(TRUE)
          { if (stmesh_face_is_between(mesh, axf, ixf, bxf)) { break; }
            axf =bxf;
            bxf = stmesh_get_next_face(mesh, axf);
            assert(axf != sxf); /* Should find a place somewhere. */
          }
      }
    (*axfP) = axf;
    (*bxfP) = bxf;
  }
  
bool_t stmesh_face_is_between
  ( stmesh_t mesh, 
    stmesh_face_ix_t axf, 
    stmesh_face_ix_t ixf, 
    stmesh_face_ix_t bxf
  );
  /* The three oriented faces with indices {axf,ixf,bxf} should have the 
    same base edge.  Returns true iff the face {ixf} should be *geometrically*
    between the faces {axf} and {bxf} in the cyclic list of faces around that edge,
    in the order defined by the right-hand rule.  
    
    Considers only the coordinates of the five vertices involved. 
    Does not use the {ixnextf} and {ixprevf} pointers of the three faces
    (in fact, it is meant to be used when those pointers are still undefined). */

bool_t stmesh_face_is_between
  ( stmesh_t mesh, 
    stmesh_face_ix_t axf, 
    stmesh_face_ix_t ixf, 
    stmesh_face_ix_t bxf
  )
  {
    demand(axf != bxf, "faces {axf,bxf} should be distinct");
    demand((axf != ixf) && (ixf != bxf), "face {ixf} should be distinct from {axf,bxf}");
    
    /* Get the corners of the five triangles: */
    stmesh_vert_ix axv[3], ixv[3], bxv[3];
    stmesh_get_face_corners(mesh, axf, axv);
    stmesh_get_face_corners(mesh, ixf, ixv);
    stmesh_get_face_corners(mesh, bxf, bxv);
    
    /* Check if they do have the same base edge: */
    demand((axv[1] == ixv[1]) && (axv[2] == ixv[2]), "faces {axf,ixf} not adjacent");
    demand((ixv[1] == bxv[1]) && (ixv[2] == bxv[2]), "faces {ixf,bxf} not adjacent");
    
    /* Get the vertex coordinates: */
    i3t *up = &(stmesh_face_get_vert_address(mesh, VERT_UNX(axv[1]))->pos); /* Origin of base edge. */
    i3t *vp = &(stmesh_face_get_vert_address(mesh, VERT_UNX(axv[2]))->pos); /* Destiny of base edge. */
    i3t *ap = &(stmesh_face_get_vert_address(mesh, VERT_UNX(axv[0]))->pos); /* Far vertex of {axf}. */
    i3t *ip = &(stmesh_face_get_vert_address(mesh, VERT_UNX(ixv[0]))->pos); /* Far vertex of {ixf}. */
    i3t *bp = &(stmesh_face_get_vert_address(mesh, VERT_UNX(bxv[0]))->pos); /* Far vertex of {bxf}. */

    sign_t sab = stmesh_tetra_sign(up, vp, ap, bp);
    sign_t sai = stmesh_tetra_sign(up, vp, ap, ip);
    sign_t sib = stmesh_tetra_sign(up, vp, ip, bp);
    demand((sab != 0) && (sai != 0) && (sbi != 0), "degenerate triangle found");
    return (sab*sai > 0) && (sab*sib > 0);
  }

uint32_t stmesh_edge_degree(stmesh_t mesh, stmesh_edge_ix_t ixe)
  {
    uint32_t nf = mesh->nf;
    
    stmesh_face_ix_t rxf = stmesh_get_first_face(mesh, ixe);
    if (rxf == IX_NULL)
      { return 0; }
    else
      { assert(rxf < nf);
        uint32_t deg = 1;
        stmesh_face_ix_t jxf = rxf;
        while(TRUE)
          { jfx = stmesh_get_next_face(mesh, jxf);
            if (jxf == rxf) { break; }
            deg++;
            assert(deg <= nf);
          }
        return deg;
      } 

  }

void stmesh_set_links(stmesh_t mesh)
  {
    uint32_t nf = mesh->nf;
    
    stmesh_face_unx_t uxf; /* Index of unoriented face. */
    for (uxf = 0; uxf < nf; uxf++)
      { stmesh_face_t *fi = stmesh_get_face_address(mesh, uxf);
        int ki;
        for (ki = 0; ki < 3; ki++)
          { /* Get the index of the oriented face that has {f->ixside[k]} as base edge: */
            stmesh_face_ix_t ixf = stmesh_get_face_ix_from_unx_k_dbit(unx, k, 0);
          }
      }
  }