/* volume_ucd.c  compute the volume of a single closed .inp ucd figure */
/*               compile with read_ucd.c                               */

#include <stdio.h>
#include <string.h>
#include <math.h>
/* _MAIN_ includes data declarations for pde_read_ucd */
#define _MAIN_
#include "pde_read_ucd.h"
#include "datread.h" /* for datprint, used for debugging */

#undef abs
#define abs(x) ((x)<0.0?(-(x)):(x))
#undef max
#define max(x,y) ((x)>(y)?(x):(y))
#undef min
#define min(x,y) ((x)<(y)?(x):(y))

static float znormal(int n, float x[], float y[], float z[], float *area);
static int ucd_to_dat_read(char *filename); /* all returned is global */

static int num_points;       /* n_vertices from read_ucd   */
static int num_polys;        /* n_cells from read_ucd      */
static dpts * data_points;   /* vert.x .y .x from read_ucd */
static int *  data_verts;    /*              from read_ucd */
static float size = 1.0;     /* auto scale */
static char filename[100] = "blivet.inp";
static int status = -1;
static int debug = 0;
static float xmin, xmax, ymin, ymax, zmin, zmax;

int main(int argc, char * argv[])
{
  float area, avgz, vol, dvol, minz, znorm, tarea;
  float offs;
  int i, j, k, m, npts, iv;
  float x[10], y[10], z[10]; /* max 10 points in polygon */

  if(argc>1)
  {
    strcpy(filename,argv[1]);
  }
  printf("volume_ucd.c reading %s \n", filename);

  status = ucd_to_dat_read(filename);
           /* returns global &data_points, &num_points,
	                     &data_verts, &num_polys, &size */
  printf("status=%d, size=%g\n", status, size);
  
  datccw(&data_points, &num_points, &data_verts, &num_polys);
  if(debug) datprint(data_points, num_points, data_verts, num_polys);

  vol = 0.0; /* total volume */
  tarea = 0.0; /* total area */
  offs = -zmin; /* z's must be positive, treat them as zero */
  k = 0;
  for(i=0; i<num_polys; i++)
  {
    npts = data_verts[k];
    k++;
    for(j=0; j<npts; j++)
    {
      iv = data_verts[k];
      x[j] = data_points[iv-1].x;
      y[j] = data_points[iv-1].y;
      z[j] = data_points[iv-1].z;
      k++;
    }
    avgz = 0.0;
    for(j=0; j<npts; j++)
      avgz = avgz + (z[j]+offs); /* z's must be positive */
    avgz = avgz/(float)npts; 
           /* accurate for triangles, approx for other */

    x[npts] = x[0];
    y[npts] = y[0];
    znorm  = znormal(npts, x, y, z, &area);
    tarea = tarea + area;
    dvol = area * znorm * avgz;
    if(debug) printf("area=%g, avgz=%g, znorm=%g, dvol=%g \n",
                     area, avgz, znorm, dvol );
    vol = vol + dvol;
    if(debug) if(avgz>offs && znorm<0.0) printf("upper normal is negative \n");
    if(debug) if(avgz<offs && znorm>0.0) printf("lower normal is positive \n");
  }
  printf("enclosing area= %g, enclosing volume= %g \n",
         2.0*((xmax-xmin)*(ymax-ymin)+(xmax-xmin)*(zmax-zmin)+
         (ymax-ymin)*(zmax-zmin)),
         (xmax-xmin)*(ymax-ymin)*(zmax-zmin));
  printf("\nfinal total area = %g, total volume = %g \n\n", tarea, vol);
  return 0;
} /* end main of volume_ucd.c */

static float znormal(int n, float x[], float y[], float z[], float *aarea)
{
  float ax, ay, az, bx, by, bz, nx, ny, nz, nzf;
  float cx, cy, cz, ss, sss, sstot, ssstot;
  int i;

  sstot = 0.0;
  ssstot = 0.0;
  for(i=2; i<n; i++) /* once for three point triangle */
  {
    ax = x[i] - x[0];
    ay = y[i] - y[0];
    az = z[i] - z[0];
    bx = x[i-1] - x[0];
    by = y[i-1] - y[0];
    bz = z[i-1] - z[0];
    nx = ay*bz-az*by;
    ny = az*bx-ax*bz;
    nz = ax*by-ay*bx;
    ss = sqrt(nx*nx+ny*ny+nz*nz);
    if(i==2) /* use only first three points for normal */
    {        /* inaccurate for non planar faces        */
      nx = nx / ss;
      ny = ny / ss;
      nz = nz / ss;
      nzf = nz;
    }
    ss = ss/2.0;
    sstot = sstot + ss;

    /* check */
    cx = x[i] - x[i-1];
    cy = y[i] - y[i-1];
    cz = z[i] - z[i-1];
    ax = sqrt(ax*ax+ay*ay+az*az);
    bx = sqrt(bx*bx+by*by+bz*bz);
    cx = sqrt(cx*cx+cy*cy+cz*cz);
    sss =(ax+bx+cx)/2.0;
    sss = sqrt(sss*(sss-ax)*(sss-bx)*(sss-cx));
    ssstot = ssstot + sss;
  }
  if(abs(ssstot-sstot)>0.0001)
     printf("bad area ss=%g, sss=%g \n", sstot, ssstot);
  *aarea = sstot;
  return nzf;
} /* end znormal */

static int ucd_to_dat_read(char *filename) /* all returned is global */
{
  /* returns global data_points, num_points, data_verts, num_polys, size */
  int i, j, k, m;
  int  ctypn[9] = { 0,      1,     2,       3,     4,
                    4,      5,     6,       8};

  pde_read_ucd(filename);
  printf("pde_read_ucd returned:\n");
  printf("1 based indexing\n");
  printf("n_vertices (boundary) num_points = %d \n",n_vertices);
  num_points = n_vertices;
  printf("n_cells (boundary faces) num_polys = %d \n",n_cells);
  num_polys = n_cells;
  data_points = (dpts *)malloc(sizeof(dpts)*(num_points));
  data_verts  = (int *)malloc(sizeof(int)*6*((num_polys)+100));
  /* only room polys<6 on avgerage */
  /* vert struct {double x; double y; double z;} */
  /* cell struct {int i; int matl; int ctyp; int v[9];} */
  if(debug) printf("boundary vertices:\n");
  for(i=1; i<=n_vertices; i++)
  {
    if(debug) 
      printf("%d  x=%f, y=%f, z=%f \n",i, vert[i].x, vert[i].y, vert[i].z);
    data_points[i-1].x= vert[i].x;
    data_points[i-1].y= vert[i].y;
    data_points[i-1].z= vert[i].z;
    if(i==1)
    {
      xmin = vert[i].x;
      xmax = vert[i].x;
      ymin = vert[i].y;
      ymax = vert[i].y;
      zmin = vert[i].z;
      zmax = vert[i].z;
    }
    else
    { 
      xmin = min(xmin,vert[i].x);
      xmax = max(xmax,vert[i].x);
      ymin = min(ymin,vert[i].y);
      ymax = max(ymax,vert[i].y);
      zmin = min(zmin,vert[i].z);
      zmax = max(zmax,vert[i].z);
    }
  }
  if(debug) printf("\n");

  if(debug) printf("boundary cells:\n");
  m = 0; /* index into data_verts, count, verts(1)..verts(count) */
  for(i=1; i<=n_cells; i++)
  {
    if(debug) 
      printf("%d  matl=%d, type=%d ", cell[i].i, cell[i].matl, cell[i].ctyp);
    k = ctypn[cell[i].ctyp];
    data_verts[m] = k;
    m++;
    for(j=1; j<=k; j++)
    {
      if(debug) printf("%d ", cell[i].v[j]);
      data_verts[m] = cell[i].v[j];
      m++;
    }
    if(debug) printf("\n");
  }
  if(debug) printf("\n");
} /* end ucd_to_dat_read */

/* end volume_ucd.c */