#include <stdio.h>
#include <tmg.h>
#include <assert.h>

#define R1 2.5
#define R2 2.0
#define R3 0.5
#define R4 0.5
#define C3 0.8
#define C4 (-0.8)
#define L0 (-0.6)
#define L1 (-0.2)
#define L2 0.2
#define L3 0.6

extern Tmg *tmg;

void outlayer (char *outfilen, int *sdom, 
  double (*level)(double, double, double *, double *));
double f0 (double x, double y, double *gxpt, double *gypt);
double f1 (double x, double y, double *gxpt, double *gypt);
double f2 (double x, double y, double *gxpt, double *gypt);
double f3 (double x, double y, double *gxpt, double *gypt);
double smooth (double rd, double *fp);
void writeofffile (char *filen);
int numberbnodes (struct front_cons *comp, int startindex);
int numberinodes (struct element_data *mesh, int startindex);
int numberelements (struct element_data *mesh, int startindex);
void offborder (struct front_cons *comp, FILE *offfile, 
                double (*level)());
int offinternal (struct element_data *mesh, FILE *offfile, 
                  double (*level)(), int startindex);
void resetalltags (void);

int
main (int argc, char *argv[])
{
  int layer0[] = {1,2,3,4,0};
  int layer1[] = {2,3,0};
  int layer2[] = {2,0};
  int layer3[] = {1,2,4,0};

  tmg_init (TMG_LABEL, "blowtorus", TMG_GRAF, 1, NULL);
  tmg_parsefile ("torusb");
  resetalltags ();
//  outlayer ("torusb_0.povinc", layer0, f0);
//  outlayer ("torusb_1.povinc", layer1, f1);
//  outlayer ("torusb_2.povinc", layer2, f2);
//  outlayer ("torusb_3.povinc", layer3, f3);

/*
 * write an 'off' description
 */
  writeofffile ("torusb.off");

  //outborder ("torusb_4.povinc", bord, -1.0, 1.0); 
  //tmg_parser (stdin);
  return (0);
}

static int nodesb1s, nodesb2s, nodesb3s, nodesb4s;
static int nodesr1s, nodesr2s, nodesr3s, nodesr4s;
static int nodesnumber;

void writeofffile (char *filen)
{
  struct bord_list *bord1, *bord2, *bord3, *bord4;
  struct front_cons *comp1, *comp2, *comp3, *comp4;
  struct mesh_list *mesh1, *mesh2, *mesh3, *mesh4;
  FILE *offfile;
  int nodesb1, nodesb2, nodesb3, nodesb4;
  int nodesr1, nodesr2, nodesr3, nodesr4;
  int elements1s, elements2s, elements3s, elements4s;
  int elements1, elements2, elements3, elements4;
  int elementsnumber, nodesnumber3d, elementsnumber3d;

  offfile = fopen (filen, "w");
  fprintf (offfile, "OFF\n");
 
  nodesb1s = elements1s = 0;
  /* attenzione, l'ordine di percorrenza e' rovesciato! */
  bord4 = tmg->mesh_boundary; /* border of the fourth region (faced) */
  bord3 = bord4->next;        /* border of the third region (faced) */
  bord2 = bord3->next;        /* border of the second region */
  bord1 = bord2->next;        /* border of the first region */
  comp1 = bord1->border;
  comp2 = comp1->cdr;
  assert (comp2->cdr == 0);
  comp3 = bord2->border->cdr;
  comp4 = comp3->cdr;
  assert (comp4->cdr == 0);

  nodesb2s = numberbnodes (comp1, nodesb1s);
  nodesb3s = numberbnodes (comp2, nodesb2s);
  nodesb4s = numberbnodes (comp3, nodesb3s);
  nodesr1s = numberbnodes (comp4, nodesb4s);

  /* numbering of nodes internal to regions */
 
  mesh1 = tmg->mesh;
  mesh2 = mesh1->next;
  mesh3 = mesh2->next;
  mesh4 = mesh3->next; 
  nodesr2s = numberinodes (mesh1->mesh, nodesr1s);
  nodesr3s = numberinodes (mesh2->mesh, nodesr2s);
  nodesr4s = numberinodes (mesh3->mesh, nodesr3s);
  nodesnumber = numberinodes (mesh4->mesh, nodesr4s);
  elements2s = numberelements (mesh1->mesh, elements1s);
  elements3s = numberelements (mesh2->mesh, elements2s);
  elements4s = numberelements (mesh3->mesh, elements3s);
  elementsnumber = numberelements (mesh4->mesh, elements4s);

  nodesb1 = nodesb2s - nodesb1s;
  nodesb2 = nodesb3s - nodesb2s;
  nodesb3 = nodesb4s - nodesb3s;
  nodesb4 = nodesr1s - nodesb4s;
  nodesr1 = nodesr2s - nodesr1s;
  nodesr2 = nodesr3s - nodesr2s;
  nodesr3 = nodesr4s - nodesr3s;
  nodesr4 = nodesnumber - nodesr4s;
  elements1 = elements2s - elements1s;
  elements2 = elements3s - elements2s;
  elements3 = elements4s - elements3s;
  elements4 = elementsnumber - elements4s;
  nodesnumber3d = nodesb1 + 3*nodesb2 + 3*nodesb3 + 3*nodesb4
                + 2*nodesr1 + 4*nodesr2 + 2*nodesr3 + 2*nodesr4;
  elementsnumber3d = 2*elements1 + 4*elements2 + 2*elements3 + 2*elements4;

  fprintf (offfile, "%d %d 0\n", nodesnumber3d, elementsnumber3d);
  printf ("nodes/elements = %d %d\n", nodesnumber3d, elementsnumber3d);

  offborder (comp1, offfile, f0);
  offborder (comp2, offfile, f1);
  offborder (comp3, offfile, f2);
  offborder (comp4, offfile, f1);
  printf ("border nodes ok...\n");

  /* now for the lowest layer */
  assert (offinternal (mesh1->mesh, offfile, f0, nodesr1s) == nodesr2s);
  offborder (comp2, offfile, f0);
  assert (offinternal (mesh2->mesh, offfile, f0, nodesr2s) == nodesr3s);
  offborder (comp3, offfile, f0);
  assert (offinternal (mesh3->mesh, offfile, f0, nodesr3s) == nodesr4s);
  offborder (comp4, offfile, f0);
  assert (offinternal (mesh4->mesh, offfile, f0, nodesr4s) == nodesnumber);
  printf ("first layer ok...\n");

  /* second layer */
  assert (offinternal (mesh2->mesh, offfile, f1, nodesr2s) == nodesr3s);
  offborder (comp3, offfile, f1);
  assert (offinternal (mesh3->mesh, offfile, f1, nodesr3s) == nodesr4s);
  printf ("second layer ok...\n");

  /* third layer */
  
  assert (offinternal (mesh2->mesh, offfile, f2, nodesr2s) == nodesr3s);
  printf ("third layer ok...\n");

  /* fourth layer */
  assert (offinternal (mesh1->mesh, offfile, f3, nodesr1s) == nodesr2s);
  offborder (comp2, offfile, f3);
  assert (offinternal (mesh2->mesh, offfile, f3, nodesr2s) == nodesr3s);
  offborder (comp4, offfile, f3);
  assert (offinternal (mesh4->mesh, offfile, f3, nodesr4s) == nodesnumber);
  printf ("fourth layer ok...\n");

  /*
   * organizzazione dei nodi, ordine 2d: b1,b2,b3,b4,r1,r2,r3,r4
   *
   *       strato  offset
   *   b1  1,4     0
   *   b2  2,3     0
   *   b3  3,4     0
   *   b4  2,3     0
   *   r1   1      0
   *   b2   1      b2+b3+b4+r1 = A
   *   r2   1      b2
   *   b3   1      A+r2 = B
   *   r3   1      b2+b3 = C
   *   b4   1      B+r3 = D
   *   r4   1      C+b4 = E
   *   r2   2      E+r2+r3+r4 = F
   *   b3   2      D+r3+b4+r4+r2 = G
   *   r3   2      F+b3
   *   r2   3
   *   r1   4
   *   b2   4
   *   r2   4
   *   b4   4
   *   r4   4
   */
  /* now we add all the triangles */

  /*
   * regione 1, coinvolge gli strati 1 e 4, quindi ogni triangolo della
   * proiezione ha due controparti in 3D (con orientazione opposta)
   */

  for (elem = mesh1->mesh; elem; elem->next)
  {
    for (i = 0; i < 2; i++)
    {
      tag = elem->vertex[i]->tag;
      which = classify (tag);
      etag = tag;
      switch (which)
      {
        case 1: /* primo bordo */
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
      }
    }
  }
  

//  for (node = tmg->nodes; node; node = node->next)
//  {
//    tag = node->tag;
//    XXXXXXXXXXXXXXXXXXXXX
  
//drawcomp (comp);
//tmg_lldump ();
//tmg_parser (stdin);

  printf ("TO BE CONTINUED\n");
  fclose (offfile);
  return;
}

int
classify (int tag)
{
  assert (tag < nodesnum);
  if (tag >= nodesr4s) return (8);
  if (tag >= nodesr3s) return (7);
  if (tag >= nodesr2s) return (6);
  if (tag >= nodesr1s) return (5);
  if (tag >= nodesb4s) return (4);
  if (tag >= nodesb3s) return (3);
  if (tag >= nodesb2s) return (2);
  assert (tag >= nodesb1s);
  return (1);
}

void
outlayer (char *outfilen, int *sdom, double (*level)())
{
  FILE *outfile;
  int dtag, i;
  double x, y, l, gx, gy;
  struct element_data *elem;
  struct mesh_list *meshl;
  struct element_data *mesh;
  struct node_data *node;

  outfile = fopen (outfilen, "w");

  while ((dtag = *sdom++) != 0)
  {
    for (meshl = tmg->mesh; meshl != 0; meshl = meshl->next)
    {
      if (meshl->tag != dtag) continue;
      mesh = meshl->mesh;
      printf ("sottodominio: %d\n", meshl->tag);
      for (elem = mesh; elem != 0; elem = elem->next)
      {
        fprintf (outfile, "smooth_triangle {\n");
        for (i = 0; i < 3; i++)
        {
          node = elem->vertex[i];
          x = node->coord[0];
          y = node->coord[1];
          l = (*level)(x,y,&gx,&gy);
          fprintf (outfile, "<%lf, %lf, %lf>, <%lf, %lf, %lf>",
            x, y, l, -gx, -gy, 1.0);
          if (i == 2) fprintf (outfile, "}\n");
                 else fprintf (outfile, ",\n");
        }
      }
    }
  }

  fclose (outfile);
}

void resetalltags (void)
{
  struct node_data *node;
  struct mesh_list *meshl;
  struct element_data *mesh;
  struct element_data *elem;

  for (node = tmg->nodes; node; node = node->next)
  {
    node->tag = -1;
  }

  for (meshl = tmg->mesh; meshl != 0; meshl = meshl->next)
  {
    mesh = meshl->mesh;
    for (elem = mesh; elem != 0; elem = elem->next)
    {
      elem->tag = -1;
    }
  }
}

int
numberbnodes (struct front_cons *comp, int index)
{
  struct front_cons *node;
  struct node_data *ndata;

  node = comp->car;
  ndata = node->card;
  assert (ndata->tag == -1);
  ndata->tag = index++;

  while (1)
  {
    node = node->cdr;
    if (node->flag == ISCOMPCELL) break;

    ndata = node->card;
    assert (ndata->tag == -1);
    ndata->tag = index++;
  }
  return (index);
}

int
numberinodes (struct element_data *mesh, int index)
{
  struct element_data *elem;
  struct node_data *node;
  int i;

  for (elem = mesh; elem; elem = elem->next)
  {
    for (i = 0; i < 3; i++)
    {
      node = elem->vertex[i];
      if (node->tag < 0) node->tag = index++;
    }
  }
  return (index);
}

int
numberelements (struct element_data *mesh, int index)
{
  struct element_data *elem;

  for (elem = mesh; elem; elem = elem->next)
  {
    assert (elem->tag == -1);
    elem->tag = index++;
  }
  return (index);
}

void
offborder (struct front_cons *comp, FILE *offfile, double (*level)())
{
  struct front_cons *node;
  struct node_data *ndata;
  double x, y, gx, gy;

  node = comp->car;
  ndata = node->card;
  assert (ndata->tag >= 0);
  x = ndata->coord[0];
  y = ndata->coord[1];
  fprintf (offfile, "%lf %lf %lf\n", x, y, level(x,y,&gx,&gy));

  while (1)
  {
    node = node->cdr;
    if (node->flag == ISCOMPCELL) break;

    ndata = node->card;
    assert (ndata->tag >= 0);
    x = ndata->coord[0];
    y = ndata->coord[1];
    fprintf (offfile, "%lf %lf %lf\n", x, y, level(x,y,&gx,&gy));
  }
  return;
}

int
offinternal (struct element_data *mesh, FILE *offfile, 
                  double (*level)(), int index)
{
  struct element_data *elem;
  struct node_data *node;
  int i;
  double x, y, gx, gy;

  for (elem = mesh; elem; elem = elem->next)
  {
    for (i = 0; i < 3; i++)
    {
      node = elem->vertex[i];
      if ((node->flag & POSFLAG) != ISINTERNAL) continue;
      if (node->tag == index)
      {
        index++;
        x = node->coord[0];
        y = node->coord[1];
        fprintf (offfile, "%lf %lf %lf\n", x, y, level(x,y,&gx,&gy));
      }
    }
  }
  return (index);
}

double
f0 (double x, double y, double *gxpt, double *gypt)
{
  double d, ro, lt, f, xnorm, fp;

  *gxpt = *gypt = 0.0;
  xnorm = sqrt (x*x + y*y);
  d = R1 - xnorm;
  ro = (L3 - L0)/2.0;
  lt = (L0 + L3)/2.0;
  f = smooth (d/ro, &fp);
  *gxpt = fp*x/xnorm;
  *gypt = fp*y/xnorm;
  return (f*L0 + (1-f)*lt);
}

double
f1 (double x, double y, double *gxpt, double *gypt)
{
  double xnorm, d, ro, lt, lin, lout, f, fp;

  *gxpt = *gypt = 0.0;

  xnorm = sqrt (x*x + y*y);
  d = R2 - xnorm;
  ro = (L2 - L1)/2.0;
  lt = (L2 + L1)/2.0;

  lout = L1;
  if (d < ro)
  {
    f = smooth (d/ro, &fp);
    *gxpt = fp*x/xnorm;
    *gypt = fp*y/xnorm;
    lout = f*L1 + (1-f)*lt;
    return (lout);
  }

  xnorm = sqrt (x*x + (y-C4)*(y-C4));
  d = xnorm - R4;
  assert (d >= -1e-3);
  ro = (L2 - L1)/2.0;
  lt = (L2 + L1)/2.0;

  lin = L1;
  if (d < ro)
  {
    f = smooth (d/ro, &fp);
    *gxpt = -fp*x/xnorm;
    *gypt = -fp*(y-C4)/xnorm;
    lin = f*L1 + (1-f)*lt;
    return (lin);
  }

  return (L1);
}

double
f2 (double x, double y, double *gxpt, double *gypt)
{
  double xnorm, d, ro, lt, lin, lout, f, fp;

  *gxpt = *gypt = 0.0;

  xnorm = sqrt (x*x + (y-C3)*(y-C3));
  d = xnorm - R3;
  if (d < -1e-3)
  {
    printf ("d got %lf in f2 at (%lf,%lf)\n", d, x, y);
  }
  assert (d >= -1e-3);
  ro = (L3 - L2)/2.0;
  lt = (L3 + L2)/2.0;

  lin = L2;
  if (d < ro)
  {
    f = smooth (d/ro, &fp);
    *gxpt = -fp*x/xnorm;
    *gypt = -fp*(y-C3)/xnorm;
    lin = f*L2 + (1-f)*lt;
    return (lin);
  }

  xnorm = sqrt (x*x + y*y);
  d = R2 - xnorm;
  ro = (L2 - L1)/2.0;
  lt = (L2 + L1)/2.0;

  lout = L2;
  if (d < ro)
  {
    f = smooth (d/ro, &fp);
    *gxpt = -fp*x/xnorm;
    *gypt = -fp*y/xnorm;
    lout = f*L2 + (1-f)*lt;
    return (lout);
  }

  xnorm = sqrt (x*x + (y-C4)*(y-C4));
  d = xnorm - R4;
  assert (d >= -1e-3);
  ro = (L2 - L1)/2.0;
  lt = (L2 + L1)/2.0;

  lin = L2;
  if (d < ro)
  {
    f = smooth (d/ro, &fp);
    *gxpt = fp*x/xnorm;
    *gypt = fp*(y-C4)/xnorm;
    lin = f*L2 + (1-f)*lt;
    return (lin);
  }

  return (L2);
}

double
f3 (double x, double y, double *gxpt, double *gypt)
{
  double d, ro, lt, f, xnorm, fp;
  double lin, lout;

  *gxpt = *gypt = 0.0;
  lin = lout = L3;
  xnorm = sqrt (x*x + (y-C3)*(y-C3));
  d = xnorm - R3;
  assert (d >= -1e-3);
  ro = (L3 - L2)/2.0;
  lt = (L3 + L2)/2.0;

  if (d < ro)
  {
    f = smooth (d/ro, &fp);
    *gxpt = fp*x/xnorm;
    *gypt = fp*(y-C3)/xnorm;
    lin = f*L3 + (1-f)*lt;
  }

  xnorm = sqrt (x*x + y*y);
  d = R1 - xnorm;
  ro = (L3 - L0)/2.0;
  lt = (L0 + L3)/2.0;
  if (d < ro)
  {
    f = smooth (d/ro, &fp);
    lout = f*L3 + (1-f)*lt;
    if (lout < lin)
    {
      *gxpt = -fp*x/xnorm;
      *gypt = -fp*y/xnorm;
      return (lout);
    }
  }
  return (lin);
}

double
smooth (double rd, double *fp)
{
  double val;

  *fp = 0.0;
  rd = fabs(rd);
  if (rd >= 1.0) return (1.0);
  val = sqrt (rd*(2 - rd));
  *fp = (1 - rd)/val;
  if (rd < 0.5e-4) *fp = 1e2;
  return (val);
}

#ifdef UNDEF
void
outborder (char *outfilen, struct front_cons *comp, 
                double l1, double l2)
{
  FILE *outfile;
  struct front_cons *node;
  struct node_data *ndata;
  double x, y, newx, newy;

  outfile = fopen (outfilen, "w");

  node = comp->car;
  ndata = node->card;
  x = ndata->coord[0];
  y = ndata->coord[1];

  while (1)
  {
    node = node->cdr;
    if (node->flag == ISCOMPCELL) break;

    ndata = node->card;
    newx = ndata->coord[0];
    newy = ndata->coord[1];

    fprintf (outfile, "triangle {\n");
    fprintf (outfile, "<%lf, %lf, %lf>,\n", x, y, l1);
    fprintf (outfile, "<%lf, %lf, %lf>,\n", newx, newy, l1);
    fprintf (outfile, "<%lf, %lf, %lf>}\n", newx, newy, l2);
    fprintf (outfile, "triangle {\n");
    fprintf (outfile, "<%lf, %lf, %lf>,\n", x, y, l1);
    fprintf (outfile, "<%lf, %lf, %lf>,\n", newx, newy, l2);
    fprintf (outfile, "<%lf, %lf, %lf>}\n", x, y, l2);
    x = newx;
    y = newy;
  }

  fclose (outfile);
}
#endif