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

#define Tridimensional_C_COPYRIGHT \
  "Copyright © 2000 Universidade Estadual de Campinas (UNICAMP)"

#include <Triangulation.h>
#include <Mis.h>
#include <Octf.h>

#include <filefmt.h>
#include <r4.h>
#include <r3.h>
#include <jsfile.h>
#include <jsmath.h>
#include <fget.h>
#include <nget.h>
#include <affirm.h>

#include <string.h>
#include <stdio.h>
#include <math.h>

#define TE_FILE_VERSION "99-08-25"
#define ST_FILE_VERSION "99-08-25"
#define EN_FILE_VERSION "99-08-25"
#define ST3_FILE_VERSION "2007-01-27"
  /* File version stamps for ".te", ".st", ".en", and ".st3" files. */

r3_t Barycenter3D(ElemTableRec_t *top, Coords3D_t *c3)
  { r3_t B = (r3_t){{0, 0, 0}};
    uint n = 0;
    int i;
    for (i = 0; i < c3->ne; i++)
      { Node_t v = OrgV(top->node.e[i]);
        if (v->exists){ r3_add(&(c3->e[i]), &B, &B); n++; }
      }
    r3_scale(1.0/n, &B, &B);
    return B;
  }

double MeanNodeDistance3D(ElemTableRec_t *top, Coords3D_t *c3)
  { double S = 0.0;
    uint n = 0;
    int i;
    for (i = 0; i < c3->ne; i++)
      { Node_t v = OrgV(top->node.e[i]);
        if (v->exists)
          { S += r3_norm_sqr(&(c3->e[i])); n++; }
      }
    return sqrt(S/n);
  }

void Displace3D(ElemTableRec_t *top, r3_t *d, Coords3D_t *c3)
  { int i;
    for (i = 0; i < c3->ne; i++)
      { if (OrgV(top->node.e[i])->exists)
          { r3_t *vc = &(c3->e[i]); r3_add(d, vc, vc); }
      }
  }

void Scale3D(ElemTableRec_t *top, double s, Coords3D_t *c3)
  { int i;
    for (i = 0; i < c3->ne; i++)
      { if (OrgV(top->node.e[i])->exists)
          { r3_t *vc = &(c3->e[i]); r3_scale(s, vc, vc); }
      }
  }

void NormalizeNodeDistance3D(ElemTableRec_t *top, Coords3D_t *c3)
  { r3_t b = Barycenter3D(top, c3);
    r3_neg(&b, &b);
    Displace3D(top, &b, c3);
    double s = MeanNodeDistance3D(top, c3);
    Scale3D(top, 1.0/s, c3);
  }

int EdgeWindingNumber(Place_t p, Coords3D_t *c3)
  { Place_t q = p;
    int wnd = 0;
    /* Edge endpoints: */
    r3_t u3 = c3->e[OrgV(p)->num];
    r3_t v3 = c3->e[OrgV(Clock(p))->num];
    r4_t u4 = (r4_t){{u3.c[0], u3.c[1], u3.c[2], 1.0}};
    r4_t v4 = (r4_t){{v3.c[0], v3.c[1], v3.c[2], 1.0}};
    r4_t pl = ChoosePlaneThroughPoints(&u3, &v3);
    do 
      { /* Endpoints of opposite edge in tetrahedron: */
        r3_t x3 = c3->e[OrgV(PrevE(q))->num];
        r3_t y3 = c3->e[OrgV(PrevE(NextF(q)))->num];
        r4_t x4 = (r4_t){{x3.c[0], x3.c[1], x3.c[2], 1.0}};
        r4_t y4 = (r4_t){{y3.c[0], y3.c[1], y3.c[2], 1.0}};
        bool_t pb = (r4_dot(&x4, &pl) <= 0.0);
        bool_t qb = (r4_dot(&y4, &pl) <= 0.0);
        if (pb != qb)
          { double det = r4_det(&u4, &v4, &x4, &y4);
            if (pb && (! qb))
              { if (det > 0.0) { wnd++; } }
            else if ((! pb) && (qb))
              { if (det < 0.0) { wnd--; } }
          }
        q = NextF(q);
      } while (q != p);
    return wnd;
  } 

r4_t ChoosePlaneThroughPoints(r3_t *u, r3_t *v)
  { r3_t d3;
    r3_sub(u, v, &d3);
    /* Find smallest coordinate in {d3}: */
    uint iMin = 0;
    int i;
    for (i = 1; i < 3; i++)
      { if (fabs(d3.c[i]) < fabs(d3.c[iMin])) { iMin = i;} }
    r4_t u4 = (r4_t){{u->c[0], u->c[1], u->c[2], 1.0}};
    r4_t v4 = (r4_t){{v->c[0], v->c[1], v->c[2], 1.0}};
    r4_t w4; r4_axis(iMin, &w4);
    r4_t pl; r4_cross(&u4, &v4, &w4, &pl);
    return pl;
  }

bool_t WallIsSilhouette(Place_t p, Coords3D_t *c3)
  { if ((PposP(p) == NULL) || (PnegP(p) == NULL)) { return FALSE; }
    FiveNodes_t t = TetraNegPosNodes(p);
    uint un = t.n[0]->num;
    uint vn = t.n[1]->num;
    uint wn = t.n[2]->num;
    uint xn = t.n[3]->num;
    uint yn = t.n[4]->num;
    double d1 = TetraDet3D(un, vn, wn, xn, c3);
    double d2 = TetraDet3D(un, vn, wn, yn, c3);
    return d1*d2 >= 0.0;
  }

double TetraDet3D(uint u, uint v, uint w, uint x, Coords3D_t *c3)
  { r4_t a = (r4_t){{c3->e[u].c[0], c3->e[u].c[1], c3->e[u].c[2], 1.0}};
    r4_t b = (r4_t){{c3->e[v].c[0], c3->e[v].c[1], c3->e[v].c[2], 1.0}};
    r4_t c = (r4_t){{c3->e[w].c[0], c3->e[w].c[1], c3->e[w].c[2], 1.0}};
    r4_t d = (r4_t){{c3->e[x].c[0], c3->e[x].c[1], c3->e[x].c[2], 1.0}};
    return r4_det(&a, &b, &c, &d);
  }

void WriteState3D(char *name, char *tag, ElemTableRec_t *top, Coords3D_t *c3, char *cmt /* DF " " */)
    char *fileName = jsprintf("%s%s.st3", name, tag);
    FILE *wr = open_write(fileName, TRUE);
    FWriteState3D(wr, top, c3, cmt);
    fclose(wr);
    free(fileName);
  }
  
void FWriteState3D(FILE *wr, ElemTableRec_t *top, Coords3D_t *c3, char *cmt /* DF " " */)
  {
    int vWidth = digits(top->node.ne - 1);
    
    /* Write the file header: */
    filefmt_write_header(wr, "state3D", ST3_FILE_VERSION);
    
    /* Write the number of coordinates: */
    fprintf(wr, "nodes = %d\n", top->node.ne);
    
    /* Write the file comment: */
    if ((cmt != NULL) && ((*cmt) != 0)) { filefmt_write_comment(wr, cmt, '|'); }
    filefmt_write_comment(wr, "\nNode data:\n", '|');
    
    /* Write the node coordinates: */
    int i;
    for (i = 0; i < top->node.ne; i++)
      { Node_t v = OrgV(top->node.e[i]);
        /* state 3D */
        fprintf(wr, "%*d ", vWidth, v->num); 
        WritePoint3D(wr, &(c3->e[v->num]));
        fprintf(wr,"\n");
      }       
    filefmt_write_footer(wr, "state3D");
    fflush(wr);
  }

Coords3D_t ReadState3D(char *name, char *tag)
    char *fileName = jsprintf("%s%s.st3", name, tag);
    FILE *rd = open_read(fileName, TRUE);
    Coords3D_t c3 = FReadState3D(rd);
    fclose(rd);
    free(fileName);
    return c3;
  }
  
Coords3D_t FReadState3D(FILE *rd)
  { /* Read and check the file header: */
    filefmt_read_header(rd, "state3D", ST3_FILE_VERSION);
    
    /* Read the number of coordinates: */
    uint nv = nget_uint32(rd, "nodes", 10);
    fget_eol(rd);
    
    /* Read the file comment: */
    (void)filefmt_read_comment(rd, '|'); /* [!!! Should return it? !!!] */

    /* Read the node coordinates: */  
    Coords3D_t c3 = r3_vec_new(nv);
    int j;
    for (j = 0; j < nv; j++)
      { int jj = fget_int32(rd);
        demand(j == jj, "node index mismathch in 3D state file"); 
        r3_t *cv = &(c3.e[nv]);
        cv->c[0] = fget_double(rd);
        cv->c[1] = fget_double(rd);
        cv->c[2] = fget_double(rd);
        fget_eol(rd);
      }

    /* Read the node coordinates: */  
    filefmt_read_footer(rd, "state3D");
    return c3;
  }

void WriteTetrahedra(char *name, char *tag, ElemTableRec_t *top, TetraData_vec_t *cell, char *cmt /* DF " " */)
    char *fileName = jsprintf("%s%s.te", name, tag);
    FILE *wr = open_write(fileName, TRUE);
    FWriteTetrahedra(wr, top, cell, cmt);
    fclose(wr);
    free(fileName);
  }
    
void FWriteTetrahedra(FILE *wr, ElemTableRec_t *top, TetraData_vec_t *cell, char *cmt /* DF " " */)
  { uint pWidth = digits(top->cell.ne - 1);
    filefmt_write_header(wr, "tetrahedra", TE_FILE_VERSION);
    fprintf(wr, "tetrahedra = %d\n", cell->ne);
    filefmt_write_comment(wr, cmt, '|');
    int i;
    for (i = 0; i < cell->ne; i++)
      { fprintf(wr, "%*d ", pWidth, i); 
        WriteTetrahedron(wr, &(cell->e[i]), 0);
      }
    filefmt_write_footer(wr, "tetrahedra");
    fflush(wr);
  }

void WriteTetrahedron(FILE *wr, TetraData_t *t, uint nbase)
  { fprintf(wr, " %10.3f", t->A.c[0][0]);
    fprintf(wr, " %10.3f", t->A.c[0][1]);
    fprintf(wr, " %10.3f", t->A.c[0][2]); fprintf(wr, " ");
    fprintf(wr, " %10.3f", t->A.c[1][0]);
    fprintf(wr, " %10.3f", t->A.c[1][1]);
    fprintf(wr, " %10.3f", t->A.c[1][2]); fprintf(wr, " ");
    fprintf(wr, " %10.3f", t->A.c[2][0]);
    fprintf(wr, " %10.3f", t->A.c[2][1]);
    fprintf(wr, " %10.3f", t->A.c[2][2]); 
    fprintf(wr, "\n");
    
    fprintf(wr, "  ");
    fprintf(wr, " %5d", nbase + t->p0);
    fprintf(wr, " %5d", nbase + t->p1);
    fprintf(wr, " %5d", nbase + t->p2);
    fprintf(wr, " %5d", nbase + t->p3);
    fprintf(wr, "\n");
    
    fprintf(wr, "  ");
    fprintf(wr, " %24.16e", t->density); fprintf(wr, " ");
    fprintf(wr, " %24.16e", t->alpha);   fprintf(wr, " ");
    fprintf(wr, " %24.16e", t->beta);    fprintf(wr, " ");
    fprintf(wr, "\n");
  }

void ReadTetrahedron(FILE *rd, TetraData_t *t)
  { t->A.c[0][0] = fget_double(rd); 
    t->A.c[1][0] = fget_double(rd); 
    t->A.c[2][0] = fget_double(rd); 
    t->A.c[0][1] = fget_double(rd); 
    t->A.c[1][1] = fget_double(rd); 
    t->A.c[2][1] = fget_double(rd); 
    t->A.c[0][2] = fget_double(rd); 
    t->A.c[1][2] = fget_double(rd); 
    t->A.c[2][2] = fget_double(rd); 
    fget_eol(rd);
  
    t->p0 = fget_int32(rd);
    t->p1 = fget_int32(rd);
    t->p2 = fget_int32(rd);
    t->p3 = fget_int32(rd);
    fget_eol(rd);

    t->density = fget_double(rd);
    t->alpha   = fget_double(rd);
    t->beta    = fget_double(rd);
    fget_eol(rd);
  }

void FReadVectors(FILE *rd, uint n, Coords3D_t *pos)
  { int i;
    for (i = 0; i < n; i++)
      { int ii = fget_int32(rd);
        demand(i == ii, "index mismatch in vector file");
        int ch = fgetc(rd);
        demand(ch != EOF, "unexpected EOF");
        pos->e[i].c[0] = fget_double(rd);
        pos->e[i].c[1] = fget_double(rd);
        pos->e[i].c[2] = fget_double(rd);
        fget_eol(rd);
      }
  }

#define Tridimensional_C_author \
  "Modula-3 version created by L.A.P.Lozada in 1999-2000.\n" \
  "Based on Modula-3 procedures by R.W.L.Lieseneld and J.Stolfi.\n" \
  "Modified by L. A. P. Lozada on 2000-02-19.\n" \
  "Converted to C by J. Stolfi on 26-01-2007"