/* tetra_int.c  integrate a function over a tetrahedron         */
/*            basically compute 4D volume under volume function */
/*            over a tetrahedron given by four points           */

/* first the hard way, element by element, int f(x,y,z) dx dy dz */
/* then using approximations from FEM handbook                   */

#include <stdio.h>
#include <math.h>
#undef  abs
#define abs(x)  ((x)<0.0?(-(x)):(x))
#undef  min
#define min(a,b) ((a)<(b)?(a):(b))
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))


static double f1(double x, double y, double z)
{
  return x+2.0*y+3.0*z+4.0; /* linear, first order */
}

static double f2(double x, double y, double z)
{
  return x*x+2.0*y*y+3.0*z*z+4.0*x*y+5.0*x*z+6.0*y*z+
         7.0*x+8.0*y+9.0*z+10.0; /* quadratic, second order */
}

static double f3(double x, double y, double z)
{
  return x*x*x + 2.0*y*y*y + 3.0*z*z*z + 4.0*x*x*y +
         5.0*x*x*z + 6.0*x*y*y + 7.0*x*y*z + 8.0*x*z*z +
         9.0*y*y*z + 10.0*y*z*z +
         11.0*x*x+12.0*y*y+13.0*z*z+14.0*x*y+15.0*x*z+16.0*y*z+
         17.0*x+18.0*y+19.0*z+20.0; /* cubic, third order */
}

static double fe(double x, double y, double z)
{
  return exp(x)*exp(y)*exp(z);
}


static double volume(double P1x, double P1y, double P1z,
                     double P2x, double P2y, double P2z,
                     double P3x, double P3y, double P3z,
                     double P4x, double P4y, double P4z);
static double dist(double a1x, double a1y, double a1z,
                   double a2x, double a2y, double a2z,
                   double a3x, double a3y, double a3z,
                   double a4x, double a4y, double a4z);
static double determinant(int n, double A[]);
static void init_tetra(double P1x, double P1y, double P1z,
                       double P2x, double P2y, double P2z,
                       double P3x, double P3y, double P3z,
                       double P4x, double P4y, double P4z);
static void next_tetra();

static double vert[4097][3] = {{0.0, 0.0, 0.0}}; /* unused, start with [1]*/
static int ivert[4097][4] = {{1,2,3,4}};
static int nvert = 4; /* increment before inserting, count not including zero*/
static int npoly = 1;

int main(int argc, char *argv[])
{

  /* simple four point tetrahedron P1, P2, P3, P4 of x,y,z */
  double P1x = 1.0;
  double P1y = 1.0;
  double P1z = 1.0;
  double P2x = 3.0;
  double P2y = 2.5;
  double P2z = 2.0;
  double P3x = 2.0;
  double P3y = 0.5;
  double P3z = 1.5;
  double P4x = 0.5;
  double P4y = 1.5;
  double P4z = 3.5;
  double v, a, h;
  double vol4d, sum;
  double x, y, z, hx, hy, hz;
  double xmin;
  double xmax;
  double ymin;
  double ymax;
  double zmin;
  double zmax;
  double f1exact = 16.4713541666666;
  double f2exact = 127.946615;
  double f3exact = 719.9307;
  double feexact = 229.5346;
  int i, j, k, np;
  int p1, p2, p3, p4;

  printf("tetra_int.c running \n");
  printf("Tetrahedron P1(%g,%g,%g)     P2(%g,%g,%g) \n",
         P1x, P1y, P1z, P2x, P2y, P2z);
  printf("            P3(%g,%g,%g) P4(%g,%g,%g) \n",
         P3x, P3y, P3z, P4x, P4y, P4z);
  printf("f(x,y,z)=x^3+2y^2+x*y*z \n");

/* old experiment
  a = area(P1x, P1y, P1z, P2x, P2y, P2z, P3x, P3y, P3z);
  printf("area123=%g \n", a);
  h = dist(P4x, P4y, P4z, P1x, P1y, P1z, P2x, P2y, P2z, P3x, P3y, P3z);
  printf("height4=%g, volume=%g \n", h, a*h/3.0);
  a = area(P1x, P1y, P1z, P2x, P2y, P2z, P4x, P4y, P4z);
  printf("area124=%g \n", a);
  h = dist(P3x, P3y, P3z, P1x, P1y, P1z, P2x, P2y, P2z, P4x, P4y, P4z);
  printf("height3=%g, volume=%g \n", h, a*h/3.0);
  a = area(P1x, P1y, P1z, P3x, P3y, P3z, P4x, P4y, P4z);
  printf("area134=%g \n", a);
  h = dist(P2x, P2y, P2z, P1x, P1y, P1z, P3x, P3y, P3z, P4x, P4y, P4z);
  printf("height2=%g, volume=%g \n", h, a*h/3.0);
  a = area(P2x, P2y, P2z, P3x, P3y, P3z, P4x, P4y, P4z);
  printf("area234=%g \n", a);
  h = dist(P1x, P1y, P1z, P2x, P2y, P2z, P3x, P3y, P3z, P4x, P4y, P4z);
  printf("height1=%g, volume=%g \n", h, a*h/3.0);
  printf("\n");

  v = volume2(P1x, P1y, P1z, P2x, P2y, P2z, P3x, P3y, P3z, P4x, P4y, P4z);
  printf("Volume by 1/3 area height = %g \n", v);
  printf("\n");
*/

  xmin = min(min(P1x,P2x),min(P3x,P4x));
  xmax = max(max(P1x,P2x),max(P3x,P4x));
  ymin = min(min(P1y,P2y),min(P3y,P4y));
  ymax = max(max(P1y,P2y),max(P3y,P4y));
  zmin = min(min(P1z,P2z),min(P3z,P4z));
  zmax = max(max(P1z,P2z),max(P3z,P4z));

  v = volume(P1x, P1y, P1z, P2x, P2y, P2z, P3x, P3y, P3z, P4x, P4y, P4z);
  printf("Volume by determinant = %g \n", v);

  printf("bounds xmin=%g, xmax=%g, ymin=%g, ymax=%g, zmin=%g, zmax=%g \n",
         xmin, xmax, ymin, ymax, zmin, zmax); 
  printf("bound volume = %g \n", (xmax-xmin)*(ymax-ymin)*(zmax-zmin));

  np = 1;
  printf("f1(x,y,z)=x+2*y+3*z+4  np=%d \n", np);
  x = (P1x+P2x+P3x+P4x)/4.0;
  y = (P1y+P2y+P3y+P4y)/4.0;
  z = (P1z+P2z+P3z+P4z)/4.0;
  sum =  f1(x,y,z);
  vol4d = v*sum/(double)np;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f1exact, vol4d-f1exact);

  /* np = 8 */
  init_tetra(P1x,P1y,P1z,P2x,P2y,P2z,P3x,P3y,P3z,P4x,P4y,P4z);
  next_tetra();
  np = npoly;
  printf("f1(x,y,z)=x+2*y+3*z+4  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    /* debug
    printf("evaluating i=%d, at x=%g, y=%g, z=%g, v=%g, f=%g \n",
            i, x, y, z, v, f1(x,y,z));
    */
    sum = sum + v * f1(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f1exact, vol4d-f1exact);

  /* np = 64   f1 */
  next_tetra();
  np = npoly;
  printf("f1(x,y,z)=x+2*y+3*z+4  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * f1(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f1exact, vol4d-f1exact);

  /* np = 512   f1 */
  next_tetra();
  np = npoly;
  printf("f1(x,y,z)=x+2*y+3*z+4  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * f1(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f1exact, vol4d-f1exact);

  printf("\n");

  /* np = 4   f2 */
  init_tetra(P1x,P1y,P1z,P2x,P2y,P2z,P3x,P3y,P3z,P4x,P4y,P4z);
  np = npoly;
  printf("f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ \n");
  printf("          7*x+8*y+9*z+10  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    /* debug
    printf("evaluating i=%d, at x=%g, y=%g, z=%g, v=%g, f=%g \n",
            i, x, y, z, v, f2(x,y,z));
    */
    sum = sum + v * f2(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f2exact, vol4d-f2exact);


  /* np = 16   f2 */
  next_tetra();
  np = npoly;
  printf("f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ \n");
  printf("          7*x+8*y+9*z+10  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    /* debug
    printf("evaluating i=%d, at x=%g, y=%g, z=%g, v=%g, f=%g \n",
            i, x, y, z, v, f2(x,y,z));
    */
    sum = sum + v * f2(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f2exact, vol4d-f2exact);


  /* np = 64   f2 */
  next_tetra();
  np = npoly;
  printf("f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ \n");
  printf("          7*x+8*y+9*z+10  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * f2(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f2exact, vol4d-f2exact);

  /* n= = 512  f2 */
  next_tetra();
  np = npoly;
  printf("f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ \n");
  printf("          7*x+8*y+9*z+10  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * f2(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f2exact, vol4d-f2exact);

  /* n= = 4096  f2 */
  next_tetra();
  np = npoly;
  printf("f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ \n");
  printf("          7*x+8*y+9*z+10  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * f2(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f2exact, vol4d-f2exact);

  printf("\n");


  /* np = 4  f3 */
  init_tetra(P1x,P1y,P1z,P2x,P2y,P2z,P3x,P3y,P3z,P4x,P4y,P4z);
  np = npoly;
  printf("f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + \n");
  printf("          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + \n");
  printf("          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * f3(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f3exact, vol4d-f3exact);


  /* np = 16   f3 */
  next_tetra();
  np = npoly;
  printf("f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + \n");
  printf("          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + \n");
  printf("          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * f3(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f3exact, vol4d-f3exact);

  /* np = 64   f3 */
  next_tetra();
  np = npoly;
  printf("f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + \n");
  printf("          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + \n");
  printf("          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * f3(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f3exact, vol4d-f3exact);

  /* np = 512 */
  next_tetra();
  np = npoly;
  printf("f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + \n");
  printf("          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + \n");
  printf("          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * f3(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f3exact, vol4d-f3exact);

  /* np = 4096 */
  next_tetra();
  np = npoly;
  printf("f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + \n");
  printf("          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + \n");
  printf("          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * f3(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, f3exact, vol4d-f3exact);

  printf("\n");

  /* np = 4   fe */
  init_tetra(P1x,P1y,P1z,P2x,P2y,P2z,P3x,P3y,P3z,P4x,P4y,P4z);
  np = npoly;
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * fe(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, feexact, vol4d-feexact);

  /* np = 16   fe */
  next_tetra();
  np = npoly;
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * fe(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, feexact, vol4d-feexact);

  /* np = 64   fe */
  next_tetra();
  np = npoly;
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * fe(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, feexact, vol4d-feexact);

  /* np = 512  fe*/
  next_tetra();
  np = npoly;
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * fe(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, feexact, vol4d-feexact);

  /* np = 4096  fe*/
  next_tetra();
  np = npoly;
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=%d \n", np);
  sum = 0.0;
  p1 = ivert[0][0];
  p2 = ivert[0][1];
  p3 = ivert[0][2];
  p4 = ivert[0][3];
  v = volume(vert[p1][0], vert[p1][1], vert[p1][2],
             vert[p2][0], vert[p2][1], vert[p2][2],
             vert[p3][0], vert[p3][1], vert[p3][2],
             vert[p4][0], vert[p4][1], vert[p4][2]);
  for(i=0; i<np; i++)
  {
    p1 = ivert[i][0];
    p2 = ivert[i][1];
    p3 = ivert[i][2];
    p4 = ivert[i][3];
    x = (vert[p1][0]+vert[p2][0]+vert[p3][0]+vert[p4][0])/4.0;
    y = (vert[p1][1]+vert[p2][1]+vert[p3][1]+vert[p4][1])/4.0;
    z = (vert[p1][2]+vert[p2][2]+vert[p3][2]+vert[p4][2])/4.0;
    sum = sum + v * fe(x,y,z);
  }
  vol4d = sum;
  printf("computed integral=%g, exact integral=%g, error=%g \n",
         vol4d, feexact, vol4d-feexact);

  printf("\n");


  printf("Using FEM Handbook method \n");

  {
    double a = 0.58541020;
    double b = 0.13819660;
    np = 4;
    v = volume(P1x, P1y, P1z, P2x, P2y, P2z, P3x, P3y, P3z, P4x, P4y, P4z);
    printf("f1(x,y,z)=x+2y+3z+4 np=%d \n", np);
    sum = 0.0;
    x = a*P1x + b*P2x+ b*P3x + b*P4x;
    y = a*P1y + b*P2y+ b*P3y + b*P4y;
    z = a*P1z + b*P2z+ b*P3z + b*P4z;
    sum = sum + f1(x,y,z);
    x = b*P1x + a*P2x+ b*P3x + b*P4x;
    y = b*P1y + a*P2y+ b*P3y + b*P4y;
    z = b*P1z + a*P2z+ b*P3z + b*P4z;
    sum = sum + f1(x,y,z);
    x = b*P1x + b*P2x+ a*P3x + b*P4x;
    y = b*P1y + b*P2y+ a*P3y + b*P4y;
    z = b*P1z + b*P2z+ a*P3z + b*P4z;
    sum = sum + f1(x,y,z);
    x = b*P1x + b*P2x+ b*P3x + a*P4x;
    y = b*P1y + b*P2y+ b*P3y + a*P4y;
    z = b*P1z + b*P2z+ b*P3z + a*P4z;
    sum = sum + f1(x,y,z);
    vol4d = v*sum/4.0;
    printf("computed integral=%g, exact integral=%g, error=%g \n",
           vol4d, f1exact, vol4d-f1exact);

    printf("f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ \n");
    printf("          7*x+8*y+9*z+10  np=%d \n", np);
    sum = 0.0;
    x = a*P1x + b*P2x+ b*P3x + b*P4x;
    y = a*P1y + b*P2y+ b*P3y + b*P4y;
    z = a*P1z + b*P2z+ b*P3z + b*P4z;
    sum = sum + f2(x,y,z);
    x = b*P1x + a*P2x+ b*P3x + b*P4x;
    y = b*P1y + a*P2y+ b*P3y + b*P4y;
    z = b*P1z + a*P2z+ b*P3z + b*P4z;
    sum = sum + f2(x,y,z);
    x = b*P1x + b*P2x+ a*P3x + b*P4x;
    y = b*P1y + b*P2y+ a*P3y + b*P4y;
    z = b*P1z + b*P2z+ a*P3z + b*P4z;
    sum = sum + f2(x,y,z);
    x = b*P1x + b*P2x+ b*P3x + a*P4x;
    y = b*P1y + b*P2y+ b*P3y + a*P4y;
    z = b*P1z + b*P2z+ b*P3z + a*P4z;
    sum = sum + f2(x,y,z);
    vol4d = v*sum/4.0;
    printf("computed integral=%g, exact integral=%g, error=%g \n",
           vol4d, f2exact, vol4d-f2exact);

    printf("f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + \n");
    printf("          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + \n");
    printf("          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=%d \n", np);
    sum = 0.0;
    x = a*P1x + b*P2x+ b*P3x + b*P4x;
    y = a*P1y + b*P2y+ b*P3y + b*P4y;
    z = a*P1z + b*P2z+ b*P3z + b*P4z;
    sum = sum + f3(x,y,z);
    x = b*P1x + a*P2x+ b*P3x + b*P4x;
    y = b*P1y + a*P2y+ b*P3y + b*P4y;
    z = b*P1z + a*P2z+ b*P3z + b*P4z;
    sum = sum + f3(x,y,z);
    x = b*P1x + b*P2x+ a*P3x + b*P4x;
    y = b*P1y + b*P2y+ a*P3y + b*P4y;
    z = b*P1z + b*P2z+ a*P3z + b*P4z;
    sum = sum + f3(x,y,z);
    x = b*P1x + b*P2x+ b*P3x + a*P4x;
    y = b*P1y + b*P2y+ b*P3y + a*P4y;
    z = b*P1z + b*P2z+ b*P3z + a*P4z;
    sum = sum + f3(x,y,z);
    vol4d = v*sum/4.0;
    printf("computed integral=%g, exact integral=%g, error=%g \n",
           vol4d, f3exact, vol4d-f3exact);

    printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=%d \n", np);
    sum = 0.0;
    x = a*P1x + b*P2x+ b*P3x + b*P4x;
    y = a*P1y + b*P2y+ b*P3y + b*P4y;
    z = a*P1z + b*P2z+ b*P3z + b*P4z;
    sum = sum + fe(x,y,z);
    x = b*P1x + a*P2x+ b*P3x + b*P4x;
    y = b*P1y + a*P2y+ b*P3y + b*P4y;
    z = b*P1z + a*P2z+ b*P3z + b*P4z;
    sum = sum + fe(x,y,z);
    x = b*P1x + b*P2x+ a*P3x + b*P4x;
    y = b*P1y + b*P2y+ a*P3y + b*P4y;
    z = b*P1z + b*P2z+ a*P3z + b*P4z;
    sum = sum + fe(x,y,z);
    x = b*P1x + b*P2x+ b*P3x + a*P4x;
    y = b*P1y + b*P2y+ b*P3y + a*P4y;
    z = b*P1z + b*P2z+ b*P3z + a*P4z;
    sum = sum + fe(x,y,z);
    vol4d = v*sum/4.0;
    printf("computed integral=%g, exact integral=%g, error=%g \n",
           vol4d, feexact, vol4d-feexact);

    printf("\n");
  }


  printf("Using tetraquad.m  method \n");

  {
    /* order=3, n=27, W,X,Y,Z specifically for this tetrahedron */
    double WXYZ[27][4]= {
      { 0.003993930296,  2.638832404787,  1.944854230113,  1.838830188027}, 
      { 0.013112023478,  2.520979588443,  2.006924056234,  1.814481377294}, 
      { 0.011492807779,  2.420986931477,  2.059587429277,  1.793822538968}, 
      { 0.019501959782,  2.200532724622,  1.271149787599,  1.643240445270}, 
      { 0.064024691359,  1.939658074662,  1.408545266678,  1.589342815379}, 
      { 0.056118224020,  1.718318015137,  1.525118967965,  1.543613172387}, 
      { 0.020954154997,  1.865060924874,  0.755500673406,  1.493537294699}, 
      { 0.068792230184,  1.494718269122,  0.950549930274,  1.417023185736}, 
      { 0.060297015141,  1.180499676089,  1.116040192099,  1.352104507609}, 
      { 0.006390288474,  2.626702222974,  1.956984411926,  1.899481097093}, 
      { 0.020979237565,  2.487244747416,  2.040658897260,  1.983155582427}, 
      { 0.018388492446,  2.368921537608,  2.111652823146,  2.054149508313}, 
      { 0.031203135652,  2.173681801304,  1.298000710917,  1.777495061860}, 
      { 0.102439506174,  1.864983874101,  1.483219467239,  1.962713818181}, 
      { 0.089789158432,  1.603067978506,  1.640369004596,  2.119863355538}, 
      { 0.033526647996,  1.826942839224,  0.793618759056,  1.684127722946}, 
      { 0.110067568294,  1.388709342015,  1.056558857381,  1.947067821271}, 
      { 0.096475224226,  1.016888513994,  1.279651354194,  2.170160318084}, 
      { 0.003993930296,  2.614572041161,  1.969114593739,  1.960132006159}, 
      { 0.013112023478,  2.453509906390,  2.074393738287,  2.151829787561}, 
      { 0.011492807779,  2.316856143738,  2.163718217015,  2.314476477658}, 
      { 0.019501959782,  2.146830877986,  1.324851634235,  1.911749678450}, 
      { 0.064024691359,  1.790309673541,  1.557893667800,  2.336084820984}, 
      { 0.056118224020,  1.487817941876,  1.755619041226,  2.696113538690}, 
      { 0.020954154997,  1.788824753575,  0.831736844705,  1.874718151193}, 
      { 0.068792230184,  1.282700414908,  1.162567784488,  2.477112456806}, 
      { 0.060297015141,  0.853277351899,  1.443262516289,  2.988216128559}}; 

      np = 27;
      printf("f1(x,y,z)=x+2y+3z+4 np=%d \n", np);
      sum = 0.0;
      for(i=0; i<np; i++)
      {
        sum = sum + WXYZ[i][0]*f1(WXYZ[i][1],WXYZ[i][2],WXYZ[i][3]);
      }
      vol4d = sum;
      printf("computed integral=%g, exact integral=%g, error=%g \n",
             vol4d, f1exact, vol4d-f1exact);

      printf("f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ \n");
      printf("          7*x+8*y+9*z+10  np=%d \n", np);
      sum = 0.0;
      for(i=0; i<np; i++)
      {
        sum = sum + WXYZ[i][0]*f2(WXYZ[i][1],WXYZ[i][2],WXYZ[i][3]);
      }
      vol4d = sum;
      printf("computed integral=%g, exact integral=%g, error=%g \n",
             vol4d, f2exact, vol4d-f2exact);

      printf("f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + \n");
      printf("          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + \n");
      printf("          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=%d \n", np);
      sum = 0.0;
      for(i=0; i<np; i++)
      {
        sum = sum + WXYZ[i][0]*f3(WXYZ[i][1],WXYZ[i][2],WXYZ[i][3]);
      }
      vol4d = sum;
      printf("computed integral=%g, exact integral=%g, error=%g \n",
             vol4d, f3exact, vol4d-f3exact);

      printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=%d \n", np);
      sum = 0.0;
      for(i=0; i<np; i++)
      {
        sum = sum + WXYZ[i][0]*fe(WXYZ[i][1],WXYZ[i][2],WXYZ[i][3]);
      }
      vol4d = sum;
      printf("computed integral=%g, exact integral=%g, error=%g \n",
             vol4d, feexact, vol4d-feexact);

  }
  return 0;
}

/* distance from first point to a plane defined by three points */
static double dist(double a1x, double a1y, double a1z,
                   double a2x, double a2y, double a2z,
                   double a3x, double a3y, double a3z,
                   double a4x, double a4y, double a4z)
{
  double d, ax, ay, az, bx, by, bz;
  double vx, vy, vz;

  ax = a3x-a2x;
  ay = a3y-a2y;
  az = a3z-a2z;
  bx = a4x-a2x;
  by = a4y-a2y;
  bz = a4z-a2z;
  vx = ay*bz-by*az;
  vy = -(ax*bz-bx*az);
  vz = ax*by-bx*ay;
  d = sqrt(vx*vx+vy*vy+vz*vz);
  vx = vx/d;
  vy = vy/d;
  vz = vz/d;
  ax = a1x-a2x;
  ay = a1y-a2y;
  az = a1z-a2x;
  d  = vx*ax + vy*ay + vz*az;
  printf("dist=%g, vx=%g, vy=%g, vz=%g, ax=%g, ay=%g, az=%g \n",
         d, vx, vy, vz, ax, ay, az);
  return d;
} /* end dist */

static double volume(double P1x, double P1y, double P1z,
                     double P2x, double P2y, double P2z,
                     double P3x, double P3y, double P3z,
                     double P4x, double P4y, double P4z)
{
  double  V, M[16];

  M[0]  = 1.0; M[1]  = P1x; M[2]  = P1y; M[3]  = P1z;
  M[4]  = 1.0; M[5]  = P2x; M[6]  = P2y; M[7]  = P2z;
  M[8]  = 1.0; M[9]  = P3x; M[10] = P3y; M[11] = P3z;
  M[12] = 1.0; M[13] = P4x; M[14] = P4y; M[15] = P4z;
  V = determinant(4,M);
  return abs(V)/6.0;
} /* end volume */

static double determinant(int n, double A[])
{
    double D = 1.0;  /* determinant */
    double *B;  /* working matrix */
    int *row;  /* row interchange indices */
    int hold , I_pivot;  /* pivot indices */
    double pivot;  /* pivot element value */
    double abs_pivot;
    int i, j, k;

    B = (double *)calloc(n*n, sizeof(double));
    row = (int *)calloc(n, sizeof(int));
    memcpy(B,A,n*n*sizeof(double));

    /* set up row interchange vectors */
    for(k=0; k<n; k++)
    {
      row[k]= k;
    }

    /* begin main reduction loop */
    for(k=0; k<n-1; k++)
    {
      /* find largest element for pivot */
      pivot = B[row[k]*n+k];
      abs_pivot = abs(pivot);
      I_pivot = k;
      for(i=k; i<n; i++)
      {
        if( abs(B[row[i]*n+k]) > abs_pivot )
        {
          I_pivot = i;
          pivot = B[row[i]*n+k];
          abs_pivot = abs(pivot);
        }
      }
      /* have pivot, interchange row pointers */
      if(I_pivot != k)
      {
        hold = row[k];
        row[k] = row[I_pivot];
        row[I_pivot] = hold;
        D = - D;
      }
      /* check for near singular */
      if(abs_pivot < 1.0E-10)
      {
        free(B);
        free(row);
        return 0.0;
      }
      else
      {
        D = D * pivot;

        /*  reduce about pivot */
        for(j=k+1; j<n; j++)
        {
          B[row[k]*n+j] = B[row[k]*n+j] / B[row[k]*n+k];
        }

        /* inner reduction loop */
        for(i=0; i<n; i++)
        {
          if(i != k)
          {
            for(j=k+1; j<n; j++)
            {
              B[row[i]*n+j] = B[row[i]*n+j] - B[row[i]*n+k]* B[row[k]*n+j];
            }
          }
        }
      }
      /* finished inner reduction */
    }
    /* end main reduction loop */
    return D * B[row[n-1]*n+n-1];
} /* end determinant */

static void init_tetra(double P1x, double P1y, double P1z,
                       double P2x, double P2y, double P2z,
                       double P3x, double P3y, double P3z,
                       double P4x, double P4y, double P4z)
{
  nvert = 4; /* increment before adding, first is 1 */
  npoly = 1; /* increment after adding */
  vert[1][0] = P1x;
  vert[1][1] = P1y;
  vert[1][2] = P1z;
  vert[2][0] = P2x;
  vert[2][1] = P2y;
  vert[2][2] = P2z;
  vert[3][0] = P3x;
  vert[3][1] = P3y;
  vert[3][2] = P3z;
  vert[4][0] = P4x;
  vert[4][1] = P4y;
  vert[4][2] = P4z;
  ivert[0][0] = 1;
  ivert[0][1] = 2;
  ivert[0][2] = 3;
  ivert[0][3] = 4;
}

static void next_tetra()
{
  int tvert[4097][4];
  int tpoly;
  int ipoly, j, k;
  int p1, p2, p3, p4, p12, p13, p14, p23, p24, p34;

  tpoly = npoly;
  for(ipoly=0; ipoly<tpoly; ipoly++) /* move old polygons to temp */
    for(j=0; j<4; j++)
      tvert[ipoly][j] = ivert[ipoly][j];

  npoly = 0;
  for(ipoly=0; ipoly<tpoly; ipoly++) /* generate 8 new for each old */
  {
    p1 = tvert[ipoly][0];
    p2 = tvert[ipoly][1];
    p3 = tvert[ipoly][2];
    p4 = tvert[ipoly][3];
    nvert++;
    p12 = nvert;
    vert[nvert][0] = (vert[p1][0]+vert[p2][0])/2.0;
    vert[nvert][1] = (vert[p1][1]+vert[p2][1])/2.0;
    vert[nvert][2] = (vert[p1][2]+vert[p2][2])/2.0;
    nvert++;
    p13 = nvert;
    vert[nvert][0] = (vert[p1][0]+vert[p3][0])/2.0;
    vert[nvert][1] = (vert[p1][1]+vert[p3][1])/2.0;
    vert[nvert][2] = (vert[p1][2]+vert[p3][2])/2.0;
    nvert++;
    p14 = nvert;
    vert[nvert][0] = (vert[p1][0]+vert[p4][0])/2.0;
    vert[nvert][1] = (vert[p1][1]+vert[p4][1])/2.0;
    vert[nvert][2] = (vert[p1][2]+vert[p4][2])/2.0;
    nvert++;
    p23 = nvert;
    vert[nvert][0] = (vert[p2][0]+vert[p3][0])/2.0;
    vert[nvert][1] = (vert[p2][1]+vert[p3][1])/2.0;
    vert[nvert][2] = (vert[p2][2]+vert[p3][2])/2.0;
    nvert++;
    p24 = nvert;
    vert[nvert][0] = (vert[p2][0]+vert[p4][0])/2.0;
    vert[nvert][1] = (vert[p2][1]+vert[p4][1])/2.0;
    vert[nvert][2] = (vert[p2][2]+vert[p4][2])/2.0;
    nvert++;
    p34 = nvert;
    vert[nvert][0] = (vert[p3][0]+vert[p4][0])/2.0;
    vert[nvert][1] = (vert[p3][1]+vert[p4][1])/2.0;
    vert[nvert][2] = (vert[p3][2]+vert[p4][2])/2.0;
    ivert[npoly][0] = p1;
    ivert[npoly][1] = p12;
    ivert[npoly][2] = p13;
    ivert[npoly][3] = p14;
    npoly++;
    ivert[npoly][0] = p2;
    ivert[npoly][1] = p12;
    ivert[npoly][2] = p23;
    ivert[npoly][3] = p24;
    npoly++;
    ivert[npoly][0] = p3;
    ivert[npoly][1] = p13;
    ivert[npoly][2] = p23;
    ivert[npoly][3] = p34;
    npoly++;
    ivert[npoly][0] = p4;
    ivert[npoly][1] = p14;
    ivert[npoly][2] = p24;
    ivert[npoly][3] = p34;
    npoly++;
    ivert[npoly][0] = p24;
    ivert[npoly][1] = p12;
    ivert[npoly][2] = p13;
    ivert[npoly][3] = p23;
    npoly++;
    ivert[npoly][0] = p13;
    ivert[npoly][1] = p12;
    ivert[npoly][2] = p14;
    ivert[npoly][3] = p24;
    npoly++;
    ivert[npoly][0] = p13;
    ivert[npoly][1] = p23;
    ivert[npoly][2] = p24;
    ivert[npoly][3] = p34;
    npoly++;
    ivert[npoly][0] = p24;
    ivert[npoly][1] = p13;
    ivert[npoly][2] = p14;
    ivert[npoly][3] = p34;
    npoly++;
  } /* end splitting i_th tetrahedron */
} /* end next_tetra */