// test_triquad.java

class test_triquad
{
  test_triquad()
  {
    double x1 = 1.0;
    double y1 = 1.0;
    double x2 = 2.5;
    double y2 = 3.0;
    double x3 = 4.5;
    double y3 = 0.5;
    double T[] = new double[6];
    double T2[] = new double[6];
    double T3[] = new double[6];
    double X[] = new double[36];
    double Y[] = new double[36];
    double W[] = new double[36];
    int nn;
    double u, a, area;

    System.out.println("test_triquad.out");
    System.out.println("x1="+x1+", y1="+y1+", x2="+x2+", y2="+y2);
    System.out.println("x3="+x3+", y3="+y3);

    T3[4] = T2[2] = T[0] = x1;
    T3[5] = T2[3] = T[1] = y1;
    T3[0] = T2[4] = T[2] = x2;
    T3[1] = T2[5] = T[3] = y2;
    T3[2] = T2[0] = T[4] = x3;
    T3[3] = T2[1] = T[5] = y3;
    System.out.println("T2[0]="+T2[0]+", T2[1]="+T2[1]+", T2[2]="+T2[2]+", T2[3]="+T2[3]);
    System.out.println("T2[4]="+T2[4]+", T2[5]="+T2[5]);
    System.out.println("T3[0]="+T3[0]+", T3[1]="+T3[1]+", T3[2]="+T3[2]+", T3[3]="+T3[3]);
    System.out.println("T3[4]="+T3[4]+", T3[5]="+T3[5]);
    a = x1*y2 - x2*y1 + x2*y3 - x3*y2 + x3*y1 - x1*y3;
    area = Math.abs(a)/2.0;
    System.out.println("area="+area);

    new triquad();

    System.out.println("triquad.tint(1, T, X, Y, W);");
    nn = triquad.tint(1, T, X, Y, W);
    for(int i=0; i<nn; i++)
    {
      System.out.println("x="+X[i]+", y="+Y[i]+", w="+W[i]);
    }
    System.out.println(" ");

    System.out.println("triquad.tint(2, T, X, Y, W);");
    nn = triquad.tint(2, T, X, Y, W);
    for(int i=0; i<nn; i++)
    {
      System.out.println("x="+X[i]+", y="+Y[i]+", w="+W[i]);
    }
    System.out.println(" ");

    System.out.println("triquad.tint(3, T, X, Y, W);");
    nn = triquad.tint(3, T, X, Y, W);
    for(int i=0; i<nn; i++)
    {
      System.out.println("x="+X[i]+", y="+Y[i]+", w="+W[i]);
    }
    System.out.println(" ");

    System.out.println("triquad.tint(4, T, X, Y, W);");
    nn = triquad.tint(4, T, X, Y, W);
    for(int i=0; i<nn; i++)
    {
      System.out.println("x="+X[i]+", y="+Y[i]+", w="+W[i]);
    }
    System.out.println(" ");

    System.out.println("triquad.tint(5, T, X, Y, W);");
    nn = triquad.tint(5, T, X, Y, W);
    for(int i=0; i<nn; i++)
    {
      System.out.println("x="+X[i]+", y="+Y[i]+", w="+W[i]);
    }
    System.out.println(" ");

    System.out.println("triquad.tint(6, T, X, Y, W);");
    nn = triquad.tint(6, T, X, Y, W);
    for(int i=0; i<nn; i++)
    {
      System.out.println("x="+X[i]+", y="+Y[i]+", w="+W[i]);
    }
    System.out.println(" ");
    System.out.println(" ");

    System.out.println("f1=x+2.0*y+3.0");
    for(int n=1; n<=6; n++)
    {
      nn = triquad.tint(n, T, X, Y, W);
      u = 0.0;
      for(int j=0; j<nn; j++)
        u = u + W[j]*f1(X[j],Y[j]);
      System.out.println(nn+" point integral of f1="+u);
    }
    System.out.println(" ");

    System.out.println("f2=x*x + 2.0*x*y + 3.0*y*y + 4.0*x + 5.0*y + 6.0");
    for(int n=1; n<=6; n++)
    {
      nn = triquad.tint(n, T, X, Y, W);
      u = 0.0;
      for(int j=0; j<nn; j++)
        u = u + W[j]*f2(X[j],Y[j]);
      System.out.println(nn+" point integral of f2="+u);
    }
    System.out.println(" ");
 
    System.out.println("f3=x*x*x + 2.0*x*x*y + 3.0*x*y*y + 4.0*y*y*y +");
    System.out.println("5.0*x*x + 6.0*x*y + 7.0*y*y + 8.0*x + 9.0*y + 10.0");
    for(int n=1; n<=6; n++)
    {
      nn = triquad.tint(n, T, X, Y, W);
      u = 0.0;
      for(int j=0; j<nn; j++)
        u = u + W[j]*f3(X[j],Y[j]);
      System.out.println(nn+" point integral of f3="+u);
    }
    System.out.println(" ");

    System.out.println("f4=x*x*x*x + 2.0*x*x*x*y + 3.0*x*x*y*y + 4.0*x*y*y*y + 5.0*y*y*y*y +");
    System.out.println("6.0*x*x*x + 7.0*x*x*y + 8.0*x*y*y + 9.0*y*y*y +");
    System.out.println("10.0*x*x + 11.0*x*y + 12.0*y*y + 13.0*x + 14.0*y + 15.0");
    for(int n=1; n<=6; n++)
    {
      nn = triquad.tint(n, T, X, Y, W);
      u = 0.0;
      for(int j=0; j<nn; j++)
        u = u + W[j]*f4(X[j],Y[j]);
      System.out.println(nn+" point integral of f4="+u);
    }
    System.out.println(" ");

    System.out.println("f5=x*x*x*x*x + 2.0*x*x*x*x*y + 3.0*x*x*x*y*y + 4.0*x*x*y*y*y +");
    System.out.println("5.0*x*y*y*y*y + 6.0*y*y*y*y*y + 7.0*x*x*x*x +");
    System.out.println("8.0*x*x*x*y + 9.0*x*x*y*y + 10.0*x*y*y*y + 11.0*y*y*y*y +");
    System.out.println("12.0*x*x*x + 13.0*x*x*y + 14.0*x*y*y + 15.0*y*y*y +");
    System.out.println("16.0*x*x + 17.0*x*y + 18.0*y*y + 19.0*x + 20.0*y + 21.0");
    for(int n=1; n<=6; n++)
    {
      nn = triquad.tint(n, T, X, Y, W);
      u = 0.0;
      for(int j=0; j<nn; j++)
        u = u + W[j]*f5(X[j],Y[j]);
      System.out.println(nn+" point integral of f5="+u);
    }
    System.out.println(" ");

    System.out.println("f6=x*x*x*x*x*x + 2.0*x*x*x*x*x*y + 3.0*x*x*x*x*y*y +");
    System.out.println("4.0*x*x*x*y*y*y + 5.0*x*x*y*y*y*y + 6.0*x*y*y*y*y*y + 7.0*y*y*y*y*y*y +");
    System.out.println("8.0*x*x*x*x*x + 9.0*x*x*x*x*y + 10.0*x*x*x*y*y + 11.0*x*x*y*y*y +");
    System.out.println("12.0*x*y*y*y*y + 13.0*y*y*y*y*y + 14.0*x*x*x*x +");
    System.out.println("15.0*x*x*x*y + 16.0*x*x*y*y + 17.0*x*y*y*y + 18.0*y*y*y*y +");
    System.out.println("19.0*x*x*x + 20.0*x*x*y + 21.0*x*y*y + 22.0*y*y*y +");
    System.out.println("23.0*x*x + 24.0*x*y + 25.0*y*y + 26.0*x + 27.0*y + 28.0");
    for(int n=1; n<=6; n++)
    {
      nn = triquad.tint(n, T, X, Y, W);
      u = 0.0;
      for(int j=0; j<nn; j++)
        u = u + W[j]*f6(X[j],Y[j]);
      System.out.println(nn+" point integral over T of f6="+u);
    }
    System.out.println(" ");

    for(int n=1; n<=6; n++)
    {
      nn = triquad.tint(n, T2, X, Y, W);
      u = 0.0;
      for(int j=0; j<nn; j++)
        u = u + W[j]*f6(X[j],Y[j]);
      System.out.println(nn+" point integral over T2 of f6="+u);
    }
    System.out.println(" ");

    for(int n=1; n<=6; n++)
    {
      nn = triquad.tint(n, T3, X, Y, W);
      u = 0.0;
      for(int j=0; j<nn; j++)
        u = u + W[j]*f6(X[j],Y[j]);
      System.out.println(nn+" point integral over T3 of f6="+u);
    }
    System.out.println(" ");

    tri_split_test(T); // equal area numerical quadrature

    System.out.println(" ");
    System.out.println("compare accuracy of this 16384 equal area integration");
    System.out.println("coordinates with the accuracy of the 36 Gauss Legendre");
    System.out.println("integration coordinates in the previous set.");
    System.out.println(" ");
    System.out.println("P.S. No, you do not get better accuracy using many ");
    System.out.println("triangles with a higher order method used in each triangle.");
    System.out.println(" ");

  } // end test_triquad constructor

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

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

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

  double f4(double x, double y)
  {
      return x*x*x*x + 2.0*x*x*x*y + 3.0*x*x*y*y + 4.0*x*y*y*y + 5.0*y*y*y*y + 
        6.0*x*x*x + 7.0*x*x*y + 8.0*x*y*y + 9.0*y*y*y +
        10.0*x*x + 11.0*x*y + 12.0*y*y + 13.0*x + 14.0*y + 15.0;
  }

  double f5(double x, double y)
  {
      return x*x*x*x*x + 2.0*x*x*x*x*y + 3.0*x*x*x*y*y + 4.0*x*x*y*y*y +
        5.0*x*y*y*y*y + 6.0*y*y*y*y*y + 7.0*x*x*x*x +
        8.0*x*x*x*y + 9.0*x*x*y*y + 10.0*x*y*y*y + 11.0*y*y*y*y + 
        12.0*x*x*x + 13.0*x*x*y + 14.0*x*y*y + 15.0*y*y*y +
        16.0*x*x + 17.0*x*y + 18.0*y*y + 19.0*x + 20.0*y + 21.0;
  }

  double f6(double x, double y)
  {
      return x*x*x*x*x*x + 2.0*x*x*x*x*x*y + 3.0*x*x*x*x*y*y +
        4.0*x*x*x*y*y*y + 5.0*x*x*y*y*y*y + 6.0*x*y*y*y*y*y + 7.0*y*y*y*y*y*y + 
        8.0*x*x*x*x*x + 9.0*x*x*x*x*y + 10.0*x*x*x*y*y + 11.0*x*x*y*y*y +
        12.0*x*y*y*y*y + 13.0*y*y*y*y*y + 14.0*x*x*x*x +
        15.0*x*x*x*y + 16.0*x*x*y*y + 17.0*x*y*y*y + 18.0*y*y*y*y + 
        19.0*x*x*x + 20.0*x*x*y + 21.0*x*y*y + 22.0*y*y*y +
        23.0*x*x + 24.0*x*y + 25.0*y*y + 26.0*x + 27.0*y + 28.0;
  }

  // tri_split_test  see how accuracy varies with splitting
  //                 use center points of new triangles
  //                 equally weighted
  void tri_split_test(double T[])
  {
    final int sz=16386;
    int    ntri;                   // number of triangles  
    int    t1[] = new int[sz];     // vertex numbers on triangles
    int    t2[] = new int[sz];     // vertex numbers on triangles
    int    t3[] = new int[sz];     // vertex numbers on triangles
    int    ncoord;                 // number of coordinates
    double xc[] = new double[sz];  // coordinates in vertex order
    double yc[] = new double[sz];  // coordinates in vertex order
    double x, y, u, a, area;
    int rtn[] = new int[2];        // return ntri, ncoord

    System.out.println(" ");
    System.out.println("tri_split_test  center point used for numerical quadrature");

    ntri = 1;
    t1[0] = 0;
    t2[0] = 1;
    t3[0] = 2;
    ncoord = 3;
    xc[0] = T[0];
    yc[0] = T[1];
    xc[1] = T[2];
    yc[1] = T[3];
    xc[2] = T[4];
    yc[2] = T[5];
    a = T[0]*T[3] - T[2]*T[1] + T[2]*T[5] - T[4]*T[3] + T[4]*T[1] - T[0]*T[5];
    area = Math.abs(a)/2.0;

    for(int nsplit=1; nsplit<=7; nsplit++)
    {
      tri_split(ntri, t1, t2, t3, ncoord, xc, yc, rtn);
      ntri = rtn[0];
      ncoord = rtn[1];
      u = 0.0;
      for(int i=0; i<ntri; i++)
      {
        x = (xc[t1[i]]+xc[t2[i]]+xc[t3[i]])/3.0;
        y = (yc[t1[i]]+yc[t2[i]]+yc[t3[i]])/3.0;
        u = u + f6(x,y);
      }
      u = area*u/(double)ntri;
      System.out.println(ntri+" triangles, integral over T3 of f6="+u);
    }
  } // end tri_split_test

  // tri_split  split triangles, each one becomes four triangles
  //            similar to the original
  void tri_split(int ntri, int t1[], int t2[], int t3[],
                 int ncoord, double xc[], double yc[], int rtn[])
  {
    int j12, j23, j31;
    int nvert, nntri;
    int i1, i2, i3;
    double x1, y1, x2, y2, x3, y3;

    nvert=ncoord;
    nntri=ntri;
    for(int i=0; i<ntri; i++)
    {
      i1=t1[i];
      i2=t2[i];
      i3=t3[i];
      x1=xc[i1];
      y1=yc[i1];
      x2=xc[i2];
      y2=yc[i2];
      x3=xc[i3];
      y3=yc[i3];

      j12=nvert;
      nvert++;
      xc[j12]=(x1+x2)/2.0;
      yc[j12]=(y1+y2)/2.0;
      for(int k=0; k<nvert-1; k++)
      {
	if(Math.abs(xc[k]-xc[j12])<1.0e-6 && Math.abs(yc[k]-yc[j12])<1.0e-6)
        {
	  j12 = k;
          nvert--;
          break;
	}
      }
      j23=nvert;
      nvert++;
      xc[j23]=(x2+x3)/2.0;
      yc[j23]=(y2+y3)/2.0;
      for(int k=0; k<nvert-1; k++)
      {
	if(Math.abs(xc[k]-xc[j23])<1.0e-6 && Math.abs(yc[k]-yc[j23])<1.0e-6)
        {
	  j23 = k;
          nvert--;
          break;
	}
      }
      j31=nvert;
      nvert++;
      xc[j31]=(x3+x1)/2.0;
      yc[j31]=(y3+y1)/2.0;
      for(int k=0; k<nvert-1; k++)
      {
	if(Math.abs(xc[k]-xc[j31])<1.0e-6 && Math.abs(yc[k]-yc[j31])<1.0e-6)
        {
	  j31 = k;
          nvert--;
          break;
	}
      }

      t1[i]=i1; t2[i]=j12; t3[i]=j31;
      t1[nntri]=j12; t2[nntri]=i2; t3[nntri]=j23;
      nntri++;
      t1[nntri]=j31; t2[nntri]=j12; t3[nntri]=j23;
      nntri++;
      t1[nntri]=j31; t2[nntri]=j23; t3[nntri]=i3;
      nntri++;
    }
    rtn[0] = nntri;
    rtn[1] = nvert;
    System.out.println("tri_split returns nntri="+nntri+", nvert="+
                       nvert);
  } // end tri_split

  public static void main (String[] args)
  {
    new test_triquad();
  }
} // end test_triquad class