/* See {stmesh.h} */
/* Last edited on 2016-04-21 18:56:18 by stolfilocal */

#define _GNU_SOURCE
#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <stdint.h>

#include <bool.h>
#include <affirm.h>
#include <i3.h>
#include <r3.h>
#include <bvtable.h>
#include <bvhash.h>

#include <stmesh_STL.h>
  
#include <stmesh.h>

/* INTERNAL PROTOTYPES */

typedef struct stmesh_vert_unx_pair_t { stmesh_vert_unx_t c[2]; } stmesh_vert_unx_pair_t;
  /* A pair of vertex indices. */

typedef struct stmesh_edge_unx_triple_t { stmesh_edge_unx_t c[3]; } stmesh_edge_unx_triple_t;
  /* A triplet of unoriented edge indices. */

stmesh_t stmesh_build
  ( float eps, 
    uint32_t nv, 
    i3_t vpos[], 
    uint32_t ne, 
    stmesh_vert_unx_pair_t endv[],
    uint32_t nf, 
    stmesh_edge_unx_triple_t side[],
    bool_t checkSorted
  );
  /* Builds a mesh data structure with fundamental length unit {eps},
    {nv} vertices, {ne} unoriented edges, and {nf} unoriented faces,
    given the vertex coordinates and the basic topological information.
    
    The quantized coordinates of the vertex with index {uxv} will be {vpos[uxv]}. 
    
    The endpoints of the unoriented edge with index {uxe} will be
    the vertices with indices {endv[uxe].c[0..1]}.  
    
    The sides of the unoriented face with index {uxf} will be the 
    unoriented edges with indices {side[uxf].c[0..2]}.  
    
    The procedure also computes derived informetion
    such as edge degrees, the fields {.minZ,.maxZ} of each triangle,
    and the bounding box of the mesh.
    
    If {checkSorted} is true, the procedure verifies 
    that the faces are sorted in order of non-decreasing
    {.minZ} field. */

/* IMPLEMENTATIONS */
  
stmesh_t stmesh_read_STL(char *fileName, bool_t binary, float eps, uint32_t nfGuess, bool_t even, bool_t checkSorted)
  {
    char *format = ((char *[2]){ "ascii", "binary" })[binary];
    fprintf(stderr, "reading mesh from file %s (%s)\n", fileName, format);
    fprintf(stderr, "expecting about %u triangles\n", nfGuess);
    fprintf(stderr, "quantizing vertex coords to%s multiples of %.8f mm\n", (even ? " even" : ""), eps);
    if (checkSorted) { fprintf(stderr, "expecting triangles sorted by {.minZ}\n"); }

    /* While reading the STL file and building the structure, the
      entries must be linked by *indices* of elements, rather
      than pointers, because the hash tables may be reallocated. Only
      after the reading is finished can we translate the indices into
      pointers.
      
      During the reading, a vertex is represented by its
      integer position (as an {i3.t}); an (unoriented) edge is represented
      by indices of the two endpoints, in increasing order (as an
      {stmesh_vert_unx_pair_t}); and a face is represented by the indices of the three
      (unoriented) edges that are its sides (as an {stmesh_edge_unx_triple_t}),
      also in increasing order.
      
      The fields {uf.minZ} and {uf.maxZ}, that are not needed during reading,
      are set later, when these pairs and triples of integers are converted to
      records of types {stmesh_vert_rep_t}, {stmesh_edge_rep_t}, and
      {stmesh_face_rep_t}. */

    /* Lookup table to uniquify faces: */
    size_t sz_face = sizeof(stmesh_edge_unx_triple_t);
    uint32_t ng_face = nfGuess;  /* Expected number of faces. */
    bvtable_t *tb_face = bvtable_new(sz_face, ng_face);
    auto uint64_t hash_face(void *fp, size_t sz);
    auto int cmp_face(void *xp, void *yp, size_t sz);
    
    /* Lookup table to uniquify edges: */
    size_t sz_edge = sizeof(stmesh_vert_unx_pair_t);
    uint32_t ng_edge = (3 * ng_face + 2)/2;
    bvtable_t *tb_edge = bvtable_new(sz_edge, ng_edge);
    auto uint64_t hash_edge(void *ep, size_t sz);
    auto int cmp_edge(void *xp, void *yp, size_t sz);

    /* Lookup table to uniquify vertices: */
    size_t sz_vert = sizeof(i3_t);
    uint32_t ng_vert = (ng_face + 1)/2;
    bvtable_t *tb_vert = bvtable_new(sz_vert, ng_vert);
    auto uint64_t hash_vert(void *vp, size_t sz);
    auto int cmp_vert(void *xp, void *yp, size_t sz);
    
    auto void process_STL_face(int line, stmesh_STL_face_t *face);
      /* Procedure that quantizes an STL face {face} and stores it in
         the mesh, if not degenerate. */

    uint32_t nf_read = 0;  /* Number of faces read from the STL file. */
    uint32_t nf_keep = 0;  /* Number of faces retained in the mesh. */
    
    int32_t prevZ = INT32_MIN; /* If {checkSorted}, the {minZ} of the next face must be at least this big. */

    stmesh_STL_read(fileName, binary, &process_STL_face);
    
    fprintf(stderr, "read %u triangles, kept %u\n", nf_read, nf_keep);

    /* Close the tables and get the basic data: */
    uint32_t nv, ne, nf; /* Number of vertices, edges, and faces. */
    i3_t *vpos;                     /* Quantized coordinates of vertices. */
    stmesh_vert_unx_pair_t *endv;   /* Edges, as pairs of vertex indices. */
    stmesh_edge_unx_triple_t *side; /* Faces, as triplets of unoriented edge indices. */
    
    bvtable_close(tb_vert, &nv, (void**)&(vpos)); 
    bvtable_close(tb_edge, &ne, (void**)&(endv)); 
    bvtable_close(tb_face, &nf, (void**)&(side)); 
    
    assert(nf == nf_keep);
    fprintf(stderr, "found %u distinct vertices and %u distinct edges\n", nv, ne);

    /* Build the mesh data structure. */
    stmesh_t mesh = stmesh_build(eps, nv, vpos, ne, endv, nf, side, checkSorted);
    
    return mesh;
    
    /* INTERNAL IMPLEMENTATIONS */
    
    void process_STL_face(int line, stmesh_STL_face_t *stl_face)
      { 
        nf_read++;
        
        /* Quantize vertices, assign indices: */
        stmesh_edge_unx_triple_t uxv; /* Indices of corner vertices, assigned or recovered. */
        int32_t minZ = INT32_MAX; /* Minimum {Z}-coordinate, to check order. */
        int k;
        for (k = 0; k < 3; k++)
          { /* Quantize the coordinates of corner {k}: */
            i3_t vposk = stmesh_STL_round_point(&(stl_face->v[k]), eps, even);
            /* Update the min {Z} coordinate of the triangle: */
            int32_t zk = vposk.c[2];
            if (zk < minZ) { minZ = zk; }
            /* Get the unique vertex index: */
            uxv.c[k] = bvtable_add(tb_vert, (void*)(&vposk), &hash_vert, &cmp_vert);
            demand(uxv.c[k] <= stmesh_nv_MAX, "too many vertices in mesh");
            /* Check for repeated vertices: */
            int i;
            for (i = 0; i < k; i++)
              { if (uxv.c[i] == uxv.c[k])
                  { fprintf(stderr, "%s:%d: !! warning: vertices %d %d coincide, triangle ignored\n", fileName, line, i, k);
                    stmesh_STL_print_triangle(stderr, stl_face);
                    fprintf(stderr, "\n");
                    return;
                  }
              }
          }

        if (checkSorted)
          { /* Check that the {.minZ} fields are non-decreasing: */
            if (minZ < prevZ)
              { fprintf(stderr, "%s:%d: ** error: triangles not sorted by {minZ}\n", fileName, line);
                exit(1);
              }
            prevZ = minZ;
          }
        
        /* Get or assign indices to edges: */
        stmesh_edge_unx_triple_t uxside; /* Indices of the three unoriented edges that bound the face. */
        for (k = 0; k < 3; k++)
          { stmesh_vert_unx_pair_t uxendvk;       /* Temporary record for edge {k}. */
            stmesh_vert_unx_t uxv0 = uxv.c[k];       /* One endpoint. */
            stmesh_vert_unx_t uxv1 = uxv.c[(k+1)%3]; /* The other endpoint. */
            assert(uxv0 != uxv1);         /* We should have discarded degenerate triangles before. */
            /* Make sure that {uxv0} is the vertex with lowest index: */
            if (uxv0 > uxv1) { stmesh_vert_unx_t t = uxv0; uxv0 = uxv1; uxv1 = t; }
            uxendvk.c[0] = uxv0;
            uxendvk.c[1] = uxv1;
            uxside.c[k] = bvtable_add(tb_edge, (void*)(&uxendvk), &hash_edge, &cmp_edge);
            demand(uxside.c[k] <= stmesh_ne_MAX, "too many edges in mesh");
          }
          
        /* Sort the edge indices in increasing order: */
        if (uxside.c[0] > uxside.c[1]) { stmesh_edge_unx_t t = uxside.c[0]; uxside.c[0] = uxside.c[1]; uxside.c[1] = t; }
        if (uxside.c[1] > uxside.c[2]) { stmesh_edge_unx_t t = uxside.c[1]; uxside.c[1] = uxside.c[2]; uxside.c[2] = t; }
        if (uxside.c[0] > uxside.c[1]) { stmesh_edge_unx_t t = uxside.c[0]; uxside.c[0] = uxside.c[1]; uxside.c[1] = t; }
        
        /* Since the three vertices are distinct, the three sides must be distinct edges too. */

        /* Assign an index to the face: */
        stmesh_face_unx_t uxf = bvtable_add(tb_face, (void *)(&uxside), &hash_face, &cmp_face);
        if (uxf < nf_keep) 
          { /* Repeated face: */
            fprintf(stderr, "%s:%d: !! repeated triangle, ignored\n", fileName, line);
          }
        else
          { nf_keep++; }
      }

    /* Hashing and comparison procedures for faces, edges, and vertices: */
    
    uint64_t hash_face(void *fp, size_t sz)
      { assert(sz == sizeof(stmesh_edge_unx_triple_t));
        stmesh_edge_unx_triple_t *f = (stmesh_edge_unx_triple_t *)fp;
        /* Requires the edge indices {f.c[0..2]} to be increasing: */
        assert((f->c[0] < f->c[1]) && (f->c[1] < f->c[2]));
        /* Hash side indices: */
        uint64_t h = bvhash_bytes(fp, sizeof(stmesh_edge_unx_triple_t));
        return h;
      }
        
    auto int cmp_face(void *xp, void *yp, size_t sz)
      { assert(sz == sizeof(stmesh_edge_unx_triple_t));
        stmesh_edge_unx_triple_t *x = (stmesh_edge_unx_triple_t *)xp;
        stmesh_edge_unx_triple_t *y = (stmesh_edge_unx_triple_t *)yp;
        /* Compare the edge indices lexicographically: */
        int k;
        for (k = 0; k < 3; k++)
          { stmesh_edge_unx_t uxex = x->c[k];
            stmesh_edge_unx_t uxey = y->c[k];
            if (uxex < uxey)
              { return -1; }
            else if (uxex > uxey)
              { return +1; }
          }
        return 0;
      }
        
    uint64_t hash_edge(void *ep, size_t sz)
      { assert(sz == sizeof(stmesh_vert_unx_pair_t));
        stmesh_vert_unx_pair_t *e = (stmesh_vert_unx_pair_t *)ep;
        /* Requires endpoint indices in increasing order: */
        assert(e->c[0] < e->c[1]);
        /* Hash the endpoint indices: */
        uint64_t h = bvhash_bytes(ep, sizeof(stmesh_vert_unx_pair_t));
        return h;
      }
        
    auto int cmp_edge(void *xp, void *yp, size_t sz)
      { assert(sz == sizeof(stmesh_vert_unx_pair_t));
        stmesh_vert_unx_pair_t *x = (stmesh_vert_unx_pair_t *)xp;
        stmesh_vert_unx_pair_t *y = (stmesh_vert_unx_pair_t *)yp;
        /* Compare the endpoint indices lexicographically: */
        int k;
        for (k = 0; k < 2; k++)
          { stmesh_vert_unx_t uxvx = x->c[k];
            stmesh_vert_unx_t uxvy = y->c[k];
            if (uxvx < uxvy)
              { return -1; }
            else if (uxvx > uxvy)
              { return +1; }
          }
        return 0;
      }
        
    uint64_t hash_vert(void *vp, size_t sz)
      { assert(sz == sizeof(i3_t));
        /* Hash the quantized coords: */
        uint64_t h = bvhash_bytes(vp, sizeof(i3_t));
        return h;
      }
        
    auto int cmp_vert(void *xp, void *yp, size_t sz)
      { assert(sz == sizeof(i3_t));
        i3_t *x = (i3_t *)xp;
        i3_t *y = (i3_t *)yp;
        /* Compare quantized coords lexicographically in order {Z,Y,X}: */
        int k;
        for (k = 0; k < 3; k++)
          { int32_t xk = x->c[2-k];
            int32_t yk = y->c[2-k];
            if (xk < yk)
              { return -1; }
            else if (xk > yk)
              { return +1; }
          }
        return 0;
      }
  }
  
void stmesh_print_edge_degrees(FILE *wr, stmesh_t mesh)
  {
    uint32_t ne = stmesh_edge_count(mesh);
    uint32_t nf = stmesh_face_count(mesh);
    
    uint32_t deg_max = nf;
    uint32_t d;
    
    /* Scan edges and count edges of the same degree: */
    uint32_t *ne_deg = notnull(malloc((deg_max+1)*sizeof(uint32_t)), "no mem"); 
    for (d = 0; d <= deg_max; d++) { ne_deg[d] = 0; }
    stmesh_edge_unx_t uxe;
    for (uxe = 0; uxe < ne; uxe++)
      { stmesh_edge_t e = stmesh_get_edge(mesh, uxe);
        d = stmesh_edge_degree(e);
        if (d == 0) 
          { fprintf(stderr, "** edge %u has zero degree\n", uxe);
            assert(d > 0);
          }
        assert(d <= deg_max);
        ne_deg[d]++;
      }
    
    /* Print degree statistics: */
    bool_t border = FALSE;
    int32_t ne_tot = 0;
    int32_t nf_tot = 0;
    for (d = 0; d <= deg_max; d++)
      { if (ne_deg[d] > 0)
          { fprintf(wr, "there are %u edges shared by %u faces\n", ne_deg[d], d);
            ne_tot += ne_deg[d];
            nf_tot += ne_deg[d]*d;
            if ((d % 2) != 0) { border = TRUE; }
          }
      }
    assert((nf_tot % 3) == 0);
    nf_tot /= 3;
    assert(ne_tot == ne);
    assert(nf_tot == nf);
    
    /* if there are edges with odd degree, the mesg is not a closed surface: */
    if (border) { fprintf(stderr, "!! warning: mesh is not closed\n"); }
    
    free(ne_deg);
  }
     
r3_t stmesh_unround_point(i3_t *p, float eps)
  { 
    r3_t r;
    r.c[0] = ((double)eps)*((double)p->c[0]);
    r.c[1] = ((double)eps)*((double)p->c[1]);
    r.c[2] = ((double)eps)*((double)p->c[2]);
    return r;
  }

void stmesh_print_bounding_box(FILE *wr, stmesh_t mesh)
  { 
    float eps = stmesh_get_eps(mesh);
    i3_t minQ, maxQ;
    stmesh_get_bounding_box(mesh, &minQ, &maxQ);
    fprintf(wr, "bounding box:\n");
    int k;
    for (k = 0; k < 3; k++)
      { int32_t minQk = minQ.c[k]; float minFk = eps*(float)minQk;
        int32_t maxQk = maxQ.c[k]; float maxFk = eps*(float)maxQk;
        fprintf(wr, "  %c [ %d _ %d ] = [ %.8f mm _ %.8f mm ]\n", "XYZ"[k], minQk, maxQk, minFk, maxFk); 
      }
  }
  
bool_t stmesh_edge_crosses_plane(stmesh_edge_t e, int32_t pZ)
  {
    bool_t verbose = FALSE;
    stmesh_vert_t v[2]; /* Endpoints of {e}. */
    stmesh_edge_get_endpoints(e, v);
    int32_t vZ[2]; /* Quantized {Z}-coords of endpoints of {e}. */
    int r;
    for (r = 0; r < 2; r++) 
      { i3_t p = stmesh_vert_get_pos(v[r]); vZ[r] = p.c[2];
        if (verbose) { fprintf(stderr, "endpoint %d at Z = %+d\n", r, vZ[r]); }
      }
    assert((vZ[0] != pZ) && (vZ[1] != pZ));
    if ((vZ[0] < pZ) && (vZ[1] > pZ)) { return TRUE; }
    if ((vZ[0] > pZ) && (vZ[1] < pZ)) { return TRUE; }
    return FALSE;
  }

r2_t stmesh_edge_plane_intersection(stmesh_edge_t e, int32_t pZ, double eps)
  {
    stmesh_vert_t v[2]; /* Endpoints of {e}. */
    stmesh_edge_get_endpoints(e, v);
    i3_t p[2]; /* Quantized coordinates of the endpoints. */
    int r;
    for (r = 0; r < 2; r++) { p[r] = stmesh_vert_get_pos(v[r]); }
    int32_t nZ = pZ - p[0].c[2];
    int32_t dZ = p[1].c[2] - p[0].c[2];
    assert(dZ != 0);
    double t = ((double)nZ)/((double)dZ);
    assert((t > 0.0) && (t < 1.0));
    r2_t u;
    int k;
    for (k = 0; k < 2; k++)
      { int32_t dk = p[1].c[k] - p[0].c[k];
        u.c[k] = eps*(p[0].c[k] + t*dk);
      }
    return u;
  }

void stmesh_face_get_sliced_sides(stmesh_t mesh, stmesh_face_t f, int32_t pZ, stmesh_edge_t e[])
  { 
    bool_t verbose = FALSE;
    
    stmesh_edge_t side[3]; /* The sides of {f} */
    stmesh_face_get_sides(f, side);
    if (verbose) { fprintf(stderr, "mesh face %u and plane at Z = %+d\n", stmesh_face_get_unx(mesh, f), pZ); }
    int k;
    int m = 0;
    for (k = 0; k < 3; k++)
      { stmesh_edge_t ek = side[k];
        if (stmesh_edge_crosses_plane(ek, pZ))
          { if (verbose) 
              { stmesh_edge_unx_t uxek = stmesh_edge_get_unx(mesh, ek);
                fprintf(stderr, "  side %d = edge %u crosses\n", k, uxek);
              }
            assert(m < 2);
            e[m] = stmesh_edge_natural(ek);
            m++;
          }
      }
    demand(m == 2, "face does not cross plane");
  }

stmesh_t stmesh_build
  ( float eps, 
    uint32_t nv, 
    i3_t vpos[], 
    uint32_t ne, 
    stmesh_vert_unx_pair_t endv[],
    uint32_t nf, 
    stmesh_edge_unx_triple_t side[],
    bool_t checkSorted
  )
  {  
    bool_t verify = TRUE;
    
    demand(nv < stmesh_nv_MAX, "too many vertices");
    demand(ne < stmesh_ne_MAX, "too many edges");
    demand(nf < stmesh_nf_MAX, "too many triangles");
   
    stmesh_t mesh = stmesh_new_desc(eps, nv,ne,nf);

    /* Build the vertex table: */
    stmesh_vert_unx_t uxv;
    for (uxv = 0; uxv < nv; uxv++)
      { i3_t *vp = &(vpos[uxv]);
        stmesh_vert_unx_t uxv2 = stmesh_add_vert(mesh, vp);
        assert(uxv2 == uxv);
      }

    /* Build the edge table: */
    stmesh_edge_unx_t uxe;
    for (uxe = 0; uxe < ne; uxe++)
      { /* Get the indices of the endpoints from the {endv} list: */
        stmesh_vert_unx_pair_t *ev = &(endv[uxe]);
        stmesh_edge_unx_t uxe2 = stmesh_add_edge(mesh, &(ev->c[0]));
        assert(uxe2 == uxe);
      }
        
    /* Fill the face table: */
    stmesh_face_unx_t uxf;
    for (uxf = 0; uxf < nf; uxf++)
      { /* Get the indices of its (unoriented) side edges, from the {side} parameter: */
        stmesh_edge_unx_triple_t *fe = &(side[uxf]);
        stmesh_face_unx_t uxf2 = stmesh_add_face(mesh, &(fe->c[0]));
        assert(uxf2 == uxf);
      }

    stmesh_print_bounding_box(stderr, mesh);
    stmesh_print_edge_degrees(stderr, mesh);
    
    /* Paranoia: */
    if (verify) { stmesh_check(mesh); }
    
    return mesh;
  }
    
void stmesh_check(stmesh_t mesh)
  {
    
    /* Check vertices */
    uint32_t nv = stmesh_vert_count(mesh);
    stmesh_vert_unx_t uxv;
    for (uxv = 0; uxv < nv; uxv++)
      { stmesh_vert_t v = stmesh_get_vert(mesh, uxv);
        stmesh_vert_check(mesh, v);
      }

    /* Check edges */
    uint32_t ne = stmesh_edge_count(mesh);
    stmesh_edge_unx_t uxe;
    for (uxe = 0; uxe < ne; uxe++)
      { stmesh_edge_t e = stmesh_get_edge(mesh, uxe);
        stmesh_edge_check(mesh, e);
      }

    uint32_t nf = stmesh_face_count(mesh);
    stmesh_face_unx_t uxf;
    for (uxf = 0; uxf < nf; uxf++)
      { stmesh_face_t f = stmesh_get_face(mesh, uxf);
        stmesh_face_check(mesh, f);
      }
  }

void stmesh_vert_check(stmesh_t mesh, stmesh_vert_t v)
  {
    /* Check the index: */
    stmesh_vert_unx_t uxv = stmesh_vert_get_unx(mesh, v);
    assert(v == stmesh_get_vert(mesh, uxv));
  }

void stmesh_edge_check(stmesh_t mesh, stmesh_edge_t e)
  {
    int k;

    /* Check the index: */
    stmesh_edge_unx_t uxe = stmesh_edge_get_unx(mesh, e);
    stmesh_edge_t e0 = stmesh_get_edge(mesh, uxe); /* Natural vesion of {e} */
    
    /* Get the endpoints of {e}: */
    stmesh_vert_t ve[2];
    stmesh_edge_get_endpoints(e, ve);
    
    /* Compare with reversed versions: */
    bool_t found = FALSE; /* Found a version of {e} that matches {e0}. */
    stmesh_edge_t c = e;
    int t;
    for (t = 0; t < 2; t++)
      { /* At this point, {c = reverse^t(e)}. */
        stmesh_edge_t c1 = stmesh_edge_reverse(e, t);
        assert(c1 == c);
        
        /* Compare with natural version: */
        if (c == e0) { found = TRUE; }

        /* Check the endpoints of {c}: */
        stmesh_vert_t vc[2];
        stmesh_edge_get_endpoints(c, vc);
        for (k = 0; k < 2; k++)
          { 
            /* Compare endpoint  {vc[k]} with the endpoint of {e}: */
            if (t == 0)
              { assert(vc[k] == ve[k]); }
            else
              { assert(vc[k] == ve[1 - k]); }
            
            /* Check single endpoint: */
            stmesh_vert_t uk = stmesh_edge_get_endpoint(c, k);
            assert(vc[k] == uk);
          }
        c = stmesh_edge_reverse(c, 1);
      }
    assert(c == e);
    assert(found);
  }
   
void stmesh_face_check(stmesh_t mesh, stmesh_face_t f)
  {
    int k;

    /* Check the index: */
    stmesh_face_unx_t uxf = stmesh_face_get_unx(mesh, f);
    stmesh_face_t f0 = stmesh_get_face(mesh, uxf); /* Natural vesion of {f} */

    /* Get the sides and corners of {f}: */
    stmesh_edge_t ef[3];
    stmesh_face_get_sides(f, ef);

    stmesh_vert_t vf[3];
    stmesh_face_get_corners(f, vf);
    
    /* Compare with flipped and shifted versions: */
    bool_t found = FALSE; /* Found a version of {f} that matches {f0}. */
    stmesh_face_t g = f;
    int s, t;
    for (s = 0; s < 3; s++)
      { /* At this point, {g = shift^s(f)}. */
        stmesh_face_t g1 = stmesh_face_shift(f, s);
        assert(g1 == g);
        
        for (t = 0; t < 2; t++)
          { /* At this point, {g = flip^t(shift^s(f))}. */
            stmesh_face_t g2 = stmesh_face_flip(g1, t);
            assert(g2 == g);
            
            /* Compare with natural version: */
            if (g == f0) { found = TRUE; }
        
            /* Get the sides and corners of {g}: */
            stmesh_edge_t eg[3];
            stmesh_face_get_sides(g, eg);

            stmesh_vert_t vg[3];
            stmesh_face_get_corners(g, vg);
            
            /* Compare the edges of {f,g}: */
            for (k = 0; k < 3; k++)
              { if (t == 0)
                  { assert(eg[k] == ef[(k+s)%3]); }
                else
                  { assert(eg[k] == stmesh_edge_reverse(ef[(3-k+s)%3], 1)); }
                stmesh_edge_check(mesh, eg[k]);
              }
            
            /* Compare the corners of {f,g}: */
            for (k = 0; k < 3; k++)
              { if (t == 0)
                  { assert(vg[k] == vf[(k+s)%3]); }
                else
                  { assert(vg[k] == vf[(3-k+s)%3]); }
              }
            g = stmesh_face_flip(g, 1);
          }
        g = stmesh_face_shift(g, 1);
      }
    assert(g == f);
    assert(found);

    /* Check consistency of elements of {f}: */
    for (k = 0; k < 3; k++)
      { int k1 = (k+1)%3;
        int k2 = (k+2)%3;
        stmesh_vert_t uk = stmesh_edge_get_endpoint(ef[k1], 1);
        stmesh_vert_t wk = stmesh_edge_get_endpoint(ef[k2], 0);
        assert(vf[k] == uk);
        assert(vf[k] == wk);
      }

  }