/* tri_int.c  integrate a function over a triangle              */
/*            basically compute a volume under surface function */
/*            over a triangle given by three points             */

/* first the hard way, element by element, int f(x,y) dx dy     */
/* then using approximations from FEM handbook                  */
/* then using my  tri_split                                     */
/* then using triquad mapped from 0,0 0,1 1,0 to triangle       */

#include <stdio.h>
#include <math.h>
#include "triquad_int.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 tpoints[3][3];     /* initial triangle */
static double  vpoints[1000][3]; /* vertices  start index at 1*/
static int     ipoints[1000][4]; /* polygons  start index at 1*/
static int nvert;
static int npoly;
static double points [1000][3]; /* centers for integration */

static void next_split();


static double f1(double x, double y)
{
  return x+2.0*y+3.0;
}

static double f2(double x, double y)
{
  return x*x+2.0*y*y+3.0*x*y+4.0*x+5.0*y+6.0;
}

static double f3(double x, double y)
{
  return x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+
         5.0*x*x+6.0*y*y+7.0*x*y+8.0*x+9.0*y+10.0;
}

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

/* distance from first point to a line defined by two points */
static double dist(double a1x, double a1y,
                   double a2x, double a2y,
                   double a3x, double a3y)
{
  double d, d1, d2, d3, a1, b1, c1, a2, b2, c2;
  double x, y;
  int out = 0;
  
  a1 = a3y-a2y;
  b1 = a2x-a3x;
  c1 = a2y*(a3x-a2x)-a2x*(a3y-a2y);
  a2 = a2x-a3x;
  b2 = a2y-a3y;
  c2 = a1y*(a3y-a2y)-a1x*(a2x-a3x);
  x = (b1*c2-b2*c1)/(b2*a1-b1*a2);
  y = (a1*c2-a2*c1)/(a2*b1-a1*b2);
  d = sqrt((a1x-x)*(a1x-x)+(a1y-y)*(a1y-y));
  return d;
}

static double area(double P1x, double P1y,
                   double P2x, double P2y,
                   double P3x, double P3y)
{
  double a;

  a = P1x*P2y - P2x*P1y + P2x*P3y - P3x*P2y + P3x*P1y - P1x*P3y;
  return abs(a)/2.0;
}

static int inside(double x,   double y,
                  double P1x, double P1y,
                  double P2x, double P2y,
                  double P3x, double P3y)
{
  double a, b, c;

  a = (y-P1y)*(P2x-P1x)-(x-P1x)*(P2y-P1y);
  b = (y-P2y)*(P3x-P2x)-(x-P2x)*(P3y-P2y);
  c = (y-P3y)*(P1x-P3x)-(x-P3x)*(P1y-P3y);
  if(a*c>0.0 && b*c>0.0) return 1;
  return 0;
}


int main(int argc, char *argv[])
{
  /* simple three point triangle P1, P2, P3 x and y */
  double P1x = 1.0;
  double P1y = 1.0;
  double P2x = 3.0;
  double P2y = 2.5;
  double P3x = 2.0;
  double P3y = 0.5;
  double L12;
  double L312;
  double A, sum;
  double x, y, xmin, xmax, ymin, ymax, hx, hy;
  int n = 100;
  int i, j, nn, i1, i2;
  double xx[36], yy[36], ww[36]; /* for triquad_int */
  double f1exact = 9.583333333;
  double f2exact = 46.406250000;
  double f3exact = 175.989583333;
  double feexact = 48.511621698;

  printf("tri_int.c running \n");
  printf("Triangle P1(%g,%g) P2(%g,%g) P3(%g,%g) \n",
         P1x, P1y, P2x, P2y, P3x, P3y);

  L12 = sqrt((P1x-P2x)*(P1x-P2x)+(P1y-P2y)*(P1y-P2y));
  printf("Length 1 2 =%g \n", L12);

  L312 = dist(P3x, P3y, P1x, P1y, P2x, P2y);
  printf("Length from 3 to 1 2 =%g \n", L312);
  A = L12 * L312 / 2.0;
  printf("Area of triangle = %g \n", A);
  A = area(P1x, P1y, P2x, P2y, P3x, P3y);
  printf("Area by points = %g \n", A);

  /* crude integration as a check */
  xmin = min(P1x, min(P2x,P3x));
  xmax = max(P1x, max(P2x,P3x));
  ymin = min(P1y, min(P2y,P3y));
  ymax = max(P1y, max(P2y,P3y));
  hx = (xmax-xmin)/(double)n;
  hy = (ymax-ymin)/(double)n;
  printf("xmin=%g, xmax=%g, ymin=%g, ymax=%g \n",
         xmin, xmax, ymin, ymax);
  printf("n=%d, hx=%g, hy=%g \n", n, hx, hy);

  printf("\n");
  printf("f1(x,y)=x+2*y+3 \n");
  nn = 0; /* f1 actual points inside triangle */
  for(i=0; i<=n; i++)
  {
    x = xmin + hx/2.0 + hx*(double)i;
    for(j=0; j<=n; j++)
    {
      y = ymin + hy/2.0 + hy*(double)j;
      if(inside(x, y, P1x, P1y, P2x, P2y, P3x, P3y))
      {
        sum = sum + f1(x,y);
        nn = nn+1;
      }
    }
  }
  sum = A*sum/(double)nn;
  printf("integral f1(x,y) over triangle = %g, nn=%d \n", sum, nn);
  printf("exact integral=%g, error=%g \n", f1exact, sum-f1exact);

  printf("f2(x,y)=x*x+2*y*y+3*x*y+4*x+5*y+6 \n");
  nn = 0; /* f2 actual points inside triangle */
  for(i=0; i<=n; i++)
  {
    x = xmin + hx/2.0 + hx*(double)i;
    for(j=0; j<=n; j++)
    {
      y = ymin + hy/2.0 + hy*(double)j;
      if(inside(x, y, P1x, P1y, P2x, P2y, P3x, P3y))
      {
        sum = sum + f2(x,y);
        nn = nn+1;
      }
    }
  }
  sum = A*sum/(double)nn;
  printf("integral f2(x,y) over triangle = %g, nn=%d \n", sum, nn);
  printf("exact integral=%g, error=%g \n", f2exact, sum-f2exact);

  printf("f3(x,y)=x*x*x+2*y*y*y+3*x*x*y+4*x*y*y+5*x*x+6*y*y+7*x*y+ \n");
  printf("        8*x+9*y+10 \n");
  nn = 0; /* f3 actual points inside triangle */
  for(i=0; i<=n; i++)
  {
    x = xmin + hx/2.0 + hx*(double)i;
    for(j=0; j<=n; j++)
    {
      y = ymin + hy/2.0 + hy*(double)j;
      if(inside(x, y, P1x, P1y, P2x, P2y, P3x, P3y))
      {
        sum = sum + f3(x,y);
        nn = nn+1;
      }
    }
  }
  sum = A*sum/(double)nn;
  printf("integral f3(x,y) over triangle = %g, nn=%d \n", sum, nn);
  printf("exact integral=%g, error=%g \n", f3exact, sum-f3exact);

  printf("fe(x,y)=exp(x)*exp(y) \n");
  nn = 0; /* fe actual points inside triangle */
  for(i=0; i<=n; i++)
  {
    x = xmin + hx/2.0 + hx*(double)i;
    for(j=0; j<=n; j++)
    {
      y = ymin + hy/2.0 + hy*(double)j;
      if(inside(x, y, P1x, P1y, P2x, P2y, P3x, P3y))
      {
        sum = sum + fe(x,y);
        nn = nn+1;
      }
    }
  }
  sum = A*sum/(double)nn;
  printf("integral fe(x,y) over triangle = %g, nn=%d \n", sum, nn);
  printf("exact integral=%g, error=%g \n", feexact, sum-feexact);

  /* approximations */
  /* first order */
  printf("\n");
  x = (P1x+P2x+P3x)/3.0;
  y = (P1y+P2y+P3y)/3.0;
  sum = A * f1(x,y);
  printf("linear,first order integral=%g, using x=%g, y=%g \n", sum, x, y);
  printf("exact integral=%g, error=%g \n", f1exact, sum-f1exact);
  sum = A * f2(x,y);
  printf("quad,first order integral=%g, using x=%g, y=%g \n", sum, x, y);
  printf("exact integral=%g, error=%g \n", f2exact, sum-f2exact);
  sum = A * f3(x,y);
  printf("cubic,first order integral=%g, using x=%g, y=%g \n", sum, x, y);
  printf("exact integral=%g, error=%g \n", f3exact, sum-f3exact);
  sum = A * fe(x,y);
  printf("exp,first order integral=%g, using x=%g, y=%g \n", sum, x, y);
  printf("exact integral=%g, error=%g \n", feexact, sum-feexact);

  /* second order */
  printf("\n");
  sum = A * (f1((P1x+P2x)/2.0,(P1y+P2y)/2.0)+
             f1((P2x+P3x)/2.0,(P2y+P3y)/2.0)+
             f1((P3x+P1x)/2.0,(P3y+P1y)/2.0))/3.0;
  printf("linear,second order integral=%g, using 3 b points \n", sum);
  printf("exact integral=%g, error=%g \n", f1exact, sum-f1exact);
  sum = A * (f2((P1x+P2x)/2.0,(P1y+P2y)/2.0)+
             f2((P2x+P3x)/2.0,(P2y+P3y)/2.0)+
             f2((P3x+P1x)/2.0,(P3y+P1y)/2.0))/3.0;
  printf("quad,second order integral=%g, using 3 b points \n", sum);
  printf("exact integral=%g, error=%g \n", f2exact, sum-f2exact);
  sum = A * (f3((P1x+P2x)/2.0,(P1y+P2y)/2.0)+
             f3((P2x+P3x)/2.0,(P2y+P3y)/2.0)+
             f3((P3x+P1x)/2.0,(P3y+P1y)/2.0))/3.0;
  printf("cubic,second order integral=%g, using 3 b points \n", sum);
  printf("exact integral=%g, error=%g \n", f3exact, sum-f3exact);
  sum = A * (fe((P1x+P2x)/2.0,(P1y+P2y)/2.0)+
             fe((P2x+P3x)/2.0,(P2y+P3y)/2.0)+
             fe((P3x+P1x)/2.0,(P3y+P1y)/2.0))/3.0;
  printf("exp,second order integral=%g, using 3 b points \n", sum);
  printf("exact integral=%g, error=%g \n", feexact, sum-feexact);

  /* second order */
  printf("\n");
  sum = A * (f1(2.0*P1x/3.0+P2x/6.0+P3x/6.0,2.0*P1y/3.0+P2y/6.0+P3y/6.0)+
             f1(P1x/6.0+2.0*P2x/3.0+P3x/3.0,P1y/6.0+2.0*P2y/3.0+P3y/6.0)+
             f1(P1x/6.0+P2x/6.0+2.0*P3x/3.0,P1y/6.0+P2y/6.0+2.0*P3y/3.0))/3.0;
  printf("linear,second order integral=%g, using 3 c points \n", sum);
  printf("exact integral=%g, error=%g \n", f1exact, sum-f1exact);
  sum = A * (f2(2.0*P1x/3.0+P2x/6.0+P3x/6.0,2.0*P1y/3.0+P2y/6.0+P3y/6.0)+
             f2(P1x/6.0+2.0*P2x/3.0+P3x/3.0,P1y/6.0+2.0*P2y/3.0+P3y/6.0)+
             f2(P1x/6.0+P2x/6.0+2.0*P3x/3.0,P1y/6.0+P2y/6.0+2.0*P3y/3.0))/3.0;
  printf("quad,second order integral=%g, using 3 c points \n", sum);
  printf("exact integral=%g, error=%g \n", f2exact, sum-f2exact);
  sum = A * (f3(2.0*P1x/3.0+P2x/6.0+P3x/6.0,2.0*P1y/3.0+P2y/6.0+P3y/6.0)+
             f3(P1x/6.0+2.0*P2x/3.0+P3x/3.0,P1y/6.0+2.0*P2y/3.0+P3y/6.0)+
             f3(P1x/6.0+P2x/6.0+2.0*P3x/3.0,P1y/6.0+P2y/6.0+2.0*P3y/3.0))/3.0;
  printf("cubic,second order integral=%g, using 3 c points \n", sum);
  printf("exact integral=%g, error=%g \n", f3exact, sum-f3exact);
  sum = A * (fe(2.0*P1x/3.0+P2x/6.0+P3x/6.0,2.0*P1y/3.0+P2y/6.0+P3y/6.0)+
             fe(P1x/6.0+2.0*P2x/3.0+P3x/3.0,P1y/6.0+2.0*P2y/3.0+P3y/6.0)+
             fe(P1x/6.0+P2x/6.0+2.0*P3x/3.0,P1y/6.0+P2y/6.0+2.0*P3y/3.0))/3.0;
  printf("exp,second order integral=%g, using 3 c points \n", sum);
  printf("exact integral=%g, error=%g \n", feexact, sum-feexact);

  printf("\n");
  printf("third order integral=%g, using 3 c and a points \n", sum);
  {
    double al1 = 0.81684757;
    double al2 = 0.10810302;
    double al3 = 0.79742699;
    double al4 = 0.05971587;
    double be1 = 0.09157621;
    double be2 = 0.44594849;
    double be3 = 0.10128651;
    double be4 = 0.47014206;
    double w1  = 0.10995174;
    double w2  = 0.22338159;

    sum = A*(f1(al1*P1x+be1*P2x+be1*P3x,al1*P1y+be1*P2y+be1*P3y)*w1+
             f1(be1*P1x+al1*P2x+be1*P3x,be1*P1y+al1*P2y+be1*P3y)*w1+
             f1(be1*P1x+be1*P2x+al1*P3x,be1*P1y+be1*P2y+al1*P3y)*w1+
             f1(al2*P1x+be2*P2x+be2*P3x,al2*P1y+be2*P2y+be2*P3y)*w2+
             f1(be2*P1x+al2*P2x+be2*P3x,be2*P1y+al2*P2y+be2*P3y)*w2+
             f1(be2*P1x+be2*P2x+al2*P3x,be2*P1y+be2*P2y+al2*P3y)*w2);
    printf("quartic O(h^5) f1 integral=%g, using 3 c and 3 d points \n", sum);
    printf("exact integral=%g, error=%g \n", f1exact, sum-f1exact);

    sum = A*(f2(al1*P1x+be1*P2x+be1*P3x,al1*P1y+be1*P2y+be1*P3y)*w1+
             f2(be1*P1x+al1*P2x+be1*P3x,be1*P1y+al1*P2y+be1*P3y)*w1+
             f2(be1*P1x+be1*P2x+al1*P3x,be1*P1y+be1*P2y+al1*P3y)*w1+
             f2(al2*P1x+be2*P2x+be2*P3x,al2*P1y+be2*P2y+be2*P3y)*w2+
             f2(be2*P1x+al2*P2x+be2*P3x,be2*P1y+al2*P2y+be2*P3y)*w2+
             f2(be2*P1x+be2*P2x+al2*P3x,be2*P1y+be2*P2y+al2*P3y)*w2);
    printf("quartic O(h^5) f2 integral=%g, using 3 c and 3 d points \n", sum);
    printf("exact integral=%g, error=%g \n", f2exact, sum-f2exact);

    sum = A*(f3(al1*P1x+be1*P2x+be1*P3x,al1*P1y+be1*P2y+be1*P3y)*w1+
             f3(be1*P1x+al1*P2x+be1*P3x,be1*P1y+al1*P2y+be1*P3y)*w1+
             f3(be1*P1x+be1*P2x+al1*P3x,be1*P1y+be1*P2y+al1*P3y)*w1+
             f3(al2*P1x+be2*P2x+be2*P3x,al2*P1y+be2*P2y+be2*P3y)*w2+
             f3(be2*P1x+al2*P2x+be2*P3x,be2*P1y+al2*P2y+be2*P3y)*w2+
             f3(be2*P1x+be2*P2x+al2*P3x,be2*P1y+be2*P2y+al2*P3y)*w2);
    printf("quartic O(h^5) f3 integral=%g, using 3 c and 3 d points \n", sum);
    printf("exact integral=%g, error=%g \n", f3exact, sum-f3exact);

    sum = A*(fe(al1*P1x+be1*P2x+be1*P3x,al1*P1y+be1*P2y+be1*P3y)*w1+
             fe(be1*P1x+al1*P2x+be1*P3x,be1*P1y+al1*P2y+be1*P3y)*w1+
             fe(be1*P1x+be1*P2x+al1*P3x,be1*P1y+be1*P2y+al1*P3y)*w1+
             fe(al2*P1x+be2*P2x+be2*P3x,al2*P1y+be2*P2y+be2*P3y)*w2+
             fe(be2*P1x+al2*P2x+be2*P3x,be2*P1y+al2*P2y+be2*P3y)*w2+
             fe(be2*P1x+be2*P2x+al2*P3x,be2*P1y+be2*P2y+al2*P3y)*w2);
    printf("quartic O(h^5) fe integral=%g, using 3 c and 3 d points \n", sum);
    printf("exact integral=%g, error=%g \n", feexact, sum-feexact);
  }


  printf("\n");
  /* initialization for next_split */
  tpoints[0][0] = P1x;
  tpoints[0][1] = P1y;
  tpoints[0][2] = 0.0;
  tpoints[1][0] = P2x;
  tpoints[1][1] = P2y;
  tpoints[1][2] = 0.0;
  tpoints[2][0] = P3x;
  tpoints[2][1] = P3y;
  tpoints[2][2] = 0.0;

  points[0][0] = (tpoints[0][0]+tpoints[1][0]+tpoints[2][0])/3.0;
  points[0][1] = (tpoints[0][1]+tpoints[1][1]+tpoints[2][1])/3.0;
  points[0][2] = (tpoints[0][2]+tpoints[1][2]+tpoints[2][2])/3.0;
  npoly = 1;
  nvert = 3;
  for(i=0; i<nvert; i++)
    for(j=0; j<3; j++)
      vpoints[i+1][j]=tpoints[i][j];
  ipoints[0][0] = 3; /* count */
  ipoints[0][1] = 1; /* vertices starting with 1 */
  ipoints[0][2] = 2;
  ipoints[0][3] = 3;

  /* initial 1 point */
  printf("nvert=%d, npoly=%d \n", nvert, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f1(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f1 split integral=%g, using %d points \n", sum, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f2(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f2 split integral=%g, using %d points \n", sum, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f3(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f3 split integral=%g, using %d points \n", sum, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + fe(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("fe split integral=%g, using %d points \n", sum, npoly);

  printf("\n");
  /* now 4 points */
  next_split();
  printf("nvert=%d, npoly=%d \n", nvert, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f1(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f1 split integral=%g, using %d points \n", sum, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f2(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f2 split integral=%g, using %d points \n", sum, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f3(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f3 split integral=%g, using %d points \n", sum, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + fe(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("fe split integral=%g, using %d points \n", sum, npoly);


  printf("\n");
  /* now 16 points */
  next_split();
  printf("nvert=%d, npoly=%d \n", nvert, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f1(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f1 split integral=%g, using %d points \n", sum, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f2(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f2 split integral=%g, using %d points \n", sum, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f3(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f3 split integral=%g, using %d points \n", sum, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + fe(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("fe split integral=%g, using %d points \n", sum, npoly);


  printf("\n");
  /* now 64 points */
  next_split();
  printf("nvert=%d, npoly=%d \n", nvert, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f1(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f1 split integral=%g, using %d points \n", sum, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f2(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f2 split integral=%g, using %d points \n", sum, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f3(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f3 split integral=%g, using %d points \n", sum, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + fe(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("fe split integral=%g, using %d points \n", sum, npoly);


  printf("\n");
  /* now 256 points */
  next_split();
  printf("nvert=%d, npoly=%d \n", nvert, npoly);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f1(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f1 split integral=%g, using %d points, error=%g \n",
         sum, npoly, sum-f1exact);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f2(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f2 split integral=%g, using %d points, error=%g \n",
         sum, npoly, sum-f2exact);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + f3(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("f3 split integral=%g, using %d points, error=%g \n",
         sum, npoly, sum-f3exact);
  sum = 0.0;
  for(i=0; i<npoly; i++)
  {
    sum = sum + fe(points[i][0], points[i][1]);
  }
  sum = A*sum/(double)npoly;
  printf("fe split integral=%g, using %d points, error=%g \n",
         sum, npoly, sum-feexact);

  /* now try coordinate shift of triquad.m */
  printf("\n");
  printf("triquad with coordinate shift\n");
  for(j=1; j<6; j++)
  {
    nn = triquad_int(j, P1x, P1y, P2x, P2y, P3x, P3y, xx, yy, ww);
    sum = 0.0;
    for(i=0; i<nn; i++) sum += ww[i]*f1(xx[i],yy[i]);
    printf("triquad f1 order=%d integral=%g \n", j, sum);
    printf("exact integral=%g, error=%g \n", f1exact, sum-f1exact);
    sum = 0.0;
    for(i=0; i<nn; i++) sum += ww[i]*f2(xx[i],yy[i]);
    printf("triquad f2 order=%d integral=%g \n", j, sum);
    printf("exact integral=%g, error=%g \n", f2exact, sum-f2exact);
    sum = 0.0;
    for(i=0; i<nn; i++) sum += ww[i]*f3(xx[i],yy[i]);
    printf("triquad f3 order=%d integral=%g \n", j, sum);
    printf("exact integral=%g, error=%g \n", f3exact, sum-f3exact);
    sum = 0.0;
    for(i=0; i<nn; i++) sum += ww[i]*fe(xx[i],yy[i]);
    printf("triquad fe order=%d integral=%g \n", j, sum);
    printf("exact integral=%g, error=%g \n", feexact, sum-feexact);
    printf("\n");
  }

  return 0;
} /* end main */

/* next_split  successively split each triangle into four  triangles */
static void next_split()
{
  int i, j, k;
  int npoly_old;
  int ip_old[1000][4];  /* polygons  start index at 1*/
  int js12, js13, js23; /* indices of new split points */
  
  npoly_old = npoly;
  for(i=0; i<npoly_old; i++)
    for(j=0; j<4; j++)
      ip_old[i][j] = ipoints[i][j];
  npoly = 0;

  for(i=0; i<npoly_old; i++)
  {
    nvert++;
    js12 = nvert;
    vpoints[nvert][0]=(vpoints[ip_old[i][1]][0]+vpoints[ip_old[i][2]][0])/2.0;
    vpoints[nvert][1]=(vpoints[ip_old[i][1]][1]+vpoints[ip_old[i][2]][1])/2.0;
    vpoints[nvert][2]=(vpoints[ip_old[i][1]][2]+vpoints[ip_old[i][2]][2])/2.0;
    nvert++;
    js13 = nvert;
    vpoints[nvert][0]=(vpoints[ip_old[i][1]][0]+vpoints[ip_old[i][3]][0])/2.0;
    vpoints[nvert][1]=(vpoints[ip_old[i][1]][1]+vpoints[ip_old[i][3]][1])/2.0;
    vpoints[nvert][2]=(vpoints[ip_old[i][1]][2]+vpoints[ip_old[i][3]][2])/2.0;
    nvert++;
    js23 = nvert;
    vpoints[nvert][0]=(vpoints[ip_old[i][2]][0]+vpoints[ip_old[i][3]][0])/2.0;
    vpoints[nvert][1]=(vpoints[ip_old[i][2]][1]+vpoints[ip_old[i][3]][1])/2.0;
    vpoints[nvert][2]=(vpoints[ip_old[i][2]][2]+vpoints[ip_old[i][3]][2])/2.0;
    ipoints[npoly][0] = 3; /* count */
    ipoints[npoly][1] = ip_old[i][1];
    ipoints[npoly][2] = js12;
    ipoints[npoly][3] = js13;
    npoly++;
    ipoints[npoly][0] = 3; /* count */
    ipoints[npoly][1] = ip_old[i][2];
    ipoints[npoly][2] = js12;
    ipoints[npoly][3] = js23;
    npoly++;
    ipoints[npoly][0] = 3; /* count */
    ipoints[npoly][1] = ip_old[i][3];
    ipoints[npoly][2] = js13;
    ipoints[npoly][3] = js23;
    npoly++;
    ipoints[npoly][0] = 3; /* count */
    ipoints[npoly][1] = js12;
    ipoints[npoly][2] = js13;
    ipoints[npoly][3] = js23;
    npoly++;
  }

  for(i=0; i<npoly; i++) /* center of triangle */
  {
    points[i][0]=(vpoints[ipoints[i][1]][0]+
                  vpoints[ipoints[i][2]][0]+
                  vpoints[ipoints[i][3]][0])/3.0;
    points[i][1]=(vpoints[ipoints[i][1]][1]+
                  vpoints[ipoints[i][2]][1]+
                  vpoints[ipoints[i][3]][1])/3.0;
    points[i][2]=(vpoints[ipoints[i][1]][2]+
                  vpoints[ipoints[i][2]][2]+
                  vpoints[ipoints[i][3]][2])/3.0;
  }
} /* end next_split */
/* end tri_int.c */