// triquad.java   compute integration coordinates and weights for order 1..6 
//                Gauss Legendre over triangle
//                1 <= n <=6 at this time
//                nn returned is n*n
//                integral = sum i=0 to nn-1 W[i]*f(X[i],Y[i])

public class triquad
{
  triquad()
  {
    System.out.println("triquad ready for integration over triangles");
    test(); // may be commented out, or just do not instantiate
  }

  public static int tint(int n, double T[], double X[], double Y[], double W[])
  {
    switch(n)
    {
      case 1: return int1(T[0], T[1], T[2], T[3], T[4], T[5], X, Y, W);
      case 2: return int2(T[0], T[1], T[2], T[3], T[4], T[5], X, Y, W);
      case 3: return int3(T[0], T[1], T[2], T[3], T[4], T[5], X, Y, W);
      case 4: return int4(T[0], T[1], T[2], T[3], T[4], T[5], X, Y, W);
      case 5: return int5(T[0], T[1], T[2], T[3], T[4], T[5], X, Y, W);
      case 6: return int6(T[0], T[1], T[2], T[3], T[4], T[5], X, Y, W);
      default: return 0;
    }
  }

  static int int1(double x1, double y1, double x2, double y2,
           double x3, double y3, double X[], double Y[], double W[])
  {
    int nn;
    double a, area;
    // data for vertices = [0 0; 0 1; 1 0] 
    final double tq[][] = {
      {0.333333333333333,  0.333333333333333,  1.000000000000000 }};
    a = x1*y2 - x2*y1 + x2*y3 - x3*y2 + x3*y1 - x1*y3;
    area = Math.abs(a)/2.0;
    nn = 1;
    for(int i=0; i<nn; i++) // map onto callers triangle 
    {
      X[i] = x1 + (x2-x1)*tq[i][0] + (x3-x1)*tq[i][1];
      Y[i] = y1 + (y2-y1)*tq[i][0] + (y3-y1)*tq[i][1];
      W[i] = tq[i][2]*area; 
    }
    return nn;
  } // end int1

  static int int2(double x1, double y1, double x2, double y2,
           double x3, double y3, double X[], double Y[], double W[])
  {
    int nn;
    double a, area;
    // data for vertices = [0 0; 0 1; 1 0] 
    final double tq[][] =  {
      {0.280019915499074,  0.075031110222608,  0.181958618256023 },
      {0.075031110222608,  0.280019915499074,  0.181958618256023 },
      {0.666390246014701,  0.178558728263616,  0.318041381743977 },
      {0.178558728263616,  0.666390246014701,  0.318041381743977 }};
    
    a = x1*y2 - x2*y1 + x2*y3 - x3*y2 + x3*y1 - x1*y3;
    area = Math.abs(a)/2.0;
    nn = 4;
    for(int i=0; i<nn; i++) // map onto callers triangle 
    {
      X[i] = x1 + (x2-x1)*tq[i][0] + (x3-x1)*tq[i][1];
      Y[i] = y1 + (y2-y1)*tq[i][0] + (y3-y1)*tq[i][1];
      W[i] = tq[i][2]*area; 
    }
    return nn;
  } // end int2

  static int int3(double x1, double y1, double x2, double y2,
           double x3, double y3, double X[], double Y[], double W[])
  {
    int nn;
    double a, area;
    // data for vertices = [0 0; 0 1; 1 0] 
    final double tq[][] =  {
      {0.188409405952072,  0.023931132287081,  0.038792766611919 },
      {0.106170269119576,  0.106170269119576,  0.062068426579070 },
      {0.023931132287081,  0.188409405952072,  0.038792766611919 },
      {0.523979067720101,  0.066554067839164,  0.127356170199770 },
      {0.295266567779633,  0.295266567779633,  0.203769872319632 },
      {0.066554067839165,  0.523979067720101,  0.127356170199770 },
      {0.808694385677670,  0.102717654809626,  0.111628840966088 },
      {0.455706020243648,  0.455706020243648,  0.178606145545742 },
      {0.102717654809626,  0.808694385677670,  0.111628840966088 }};
    
    a = x1*y2 - x2*y1 + x2*y3 - x3*y2 + x3*y1 - x1*y3;
    area = Math.abs(a)/2.0;
    nn = 9;
    for(int i=0; i<nn; i++) // map onto callers triangle 
    {
      X[i] = x1 + (x2-x1)*tq[i][0] + (x3-x1)*tq[i][1];
      Y[i] = y1 + (y2-y1)*tq[i][0] + (y3-y1)*tq[i][1];
      W[i] = tq[i][2]*area; 
    }
    return nn;
  } // end int3

  static int int4(double x1, double y1, double x2, double y2,
           double x3, double y3, double X[], double Y[], double W[])
  {
    int nn;
    double a, area;
    // data for vertices = [0 0; 0 1; 1 0] 
    final double tq[][] =  {
      {0.130056079216834,  0.009703785126946,  0.010846451821050 },
      {0.093637784437329,  0.046122079906452,  0.020334519128958 },
      {0.046122079906452,  0.093637784437329,  0.020334519128958 },
      {0.009703785126946,  0.130056079216834,  0.010846451821050 },
      {0.387497483406694,  0.028912084224389,  0.045168098564740 },
      {0.278990463496509,  0.137419104134574,  0.084679449043493 },
      {0.137419104134574,  0.278990463496509,  0.084679449043493 },
      {0.028912084224389,  0.387497483406694,  0.045168098564740 },
      {0.672946863150506,  0.050210123211370,  0.070776135796172 },
      {0.484508326630433,  0.238648659731443,  0.132688432214099 },
      {0.238648659731443,  0.484508326630433,  0.132688432214099 },
      {0.050210123211370,  0.672946863150506,  0.070776135796172 },
      {0.877428809330468,  0.065466994555014,  0.047136736386764 },
      {0.631731251641125,  0.311164552244357,  0.088370177044723 },
      {0.311164552244357,  0.631731251641125,  0.088370177044723 },
      {0.065466994555014,  0.877428809330468,  0.047136736386764 }};
    
    a = x1*y2 - x2*y1 + x2*y3 - x3*y2 + x3*y1 - x1*y3;
    area = Math.abs(a)/2.0;
    nn = 16;
    for(int i=0; i<nn; i++) // map onto callers triangle 
    {
      X[i] = x1 + (x2-x1)*tq[i][0] + (x3-x1)*tq[i][1];
      Y[i] = y1 + (y2-y1)*tq[i][0] + (y3-y1)*tq[i][1];
      W[i] = tq[i][2]*area; 
    }
    return nn;
  } // end int4

  static int int5(double x1, double y1, double x2, double y2,
           double x3, double y3, double X[], double Y[], double W[])
  {
    int nn;
    double a, area;
    // data for vertices = [0 0; 0 1; 1 0] 
    final double tq[][] =  {
      {0.093912797333780,  0.004622288465046,  0.003731104333756 },
      {0.075796602735062,  0.022738483063764,  0.007537403390655 },
      {0.049267542899413,  0.049267542899413,  0.008958813594563 },
      {0.022738483063764,  0.075796602735062,  0.007537403390655 },
      {0.004622288465046,  0.093912797333780,  0.003731104333756 },
      {0.290249932250793,  0.014285794395571,  0.017510998364328 },
      {0.234259434638082,  0.070276292008282,  0.035374904220967 },
      {0.152267863323182,  0.152267863323182,  0.042045934974644 },
      {0.070276292008282,  0.234259434638082,  0.035374904220967 },
      {0.014285794395571,  0.290249932250793,  0.017510998364328 },
      {0.535660544808143,  0.026364644944471,  0.034683012862731 },
      {0.432329252970360,  0.129695936782254,  0.070065009006743 },
      {0.281012594876307,  0.281012594876307,  0.083277930430390 },
      {0.129695936782254,  0.432329252970360,  0.070065009006743 },
      {0.026364644944471,  0.535660544808143,  0.034683012862731 },
      {0.764365329781281,  0.037621252345111,  0.039608166264095 },
      {0.616915871859003,  0.185070710267389,  0.080014574772321 },
      {0.400993291063196,  0.400993291063196,  0.095103794115908 },
      {0.185070710267390,  0.616915871859003,  0.080014574772321 },
      {0.037621252345111,  0.764365329781281,  0.039608166264095 },
      {0.915147549378727,  0.045042593569804,  0.022930160703185 },
      {0.738611533396152,  0.221578609552379,  0.046322443858997 },
      {0.480095071474266,  0.480095071474266,  0.055057971328939 },
      {0.221578609552379,  0.738611533396152,  0.046322443858997 },
      {0.045042593569804,  0.915147549378727,  0.022930160703185 }};
          
    a = x1*y2 - x2*y1 + x2*y3 - x3*y2 + x3*y1 - x1*y3;
    area = Math.abs(a)/2.0;
    nn = 25;
    for(int i=0; i<nn; i++) // map onto callers triangle 
    {
      X[i] = x1 + (x2-x1)*tq[i][0] + (x3-x1)*tq[i][1];
      Y[i] = y1 + (y2-y1)*tq[i][0] + (y3-y1)*tq[i][1];
      W[i] = tq[i][2]*area; 
    }
    return nn;
  } // end int5

  static int int6(double x1, double y1, double x2, double y2,
           double x3, double y3, double X[], double Y[], double W[])
  {
    int nn;
    double a, area;
    // data for vertices = [0 0; 0 1; 1 0] 
    final double tq[][] =  {
      {0.070587631527589,  0.002466697152670,  0.001497085122473 },
      {0.060679268262819,  0.012375060417440,  0.003152443508047 },
      {0.045243246564898,  0.027811082115361,  0.004088773183090 },
      {0.027811082115361,  0.045243246564898,  0.004088773183090 },
      {0.012375060417440,  0.060679268262819,  0.003152443508047 },
      {0.002466697152670,  0.070587631527589,  0.001497085122473 },
      {0.222974263268659,  0.007791874701286,  0.007530596425383 },
      {0.191675437237121,  0.039090700732824,  0.015857334667593 },
      {0.142915682993948,  0.087850454975997,  0.020567234457533 },
      {0.087850454975997,  0.142915682993948,  0.020567234457533 },
      {0.039090700732824,  0.191675437237121,  0.015857334667593 },
      {0.007791874701286,  0.222974263268659,  0.007530596425383 },
      {0.426426917861779,  0.014901563366671,  0.016903071593886 },
      {0.366569507765801,  0.074758973462649,  0.035593151994053 },
      {0.273318962107258,  0.168009519121192,  0.046164927302717 },
      {0.168009519121192,  0.273318962107258,  0.046164927302717 },
      {0.074758973462649,  0.366569507765801,  0.035593151994053 },
      {0.014901563366671,  0.426426917861779,  0.016903071593886 },
      {0.640628436740815,  0.022386872978031,  0.024121212808530 },
      {0.550703627937892,  0.112311681780954,  0.050792543178095 },
      {0.410611741642328,  0.252403568076518,  0.065878797801573 },
      {0.252403568076518,  0.410611741642328,  0.065878797801573 },
      {0.112311681780954,  0.550703627937892,  0.050792543178095 },
      {0.022386872978031,  0.640628436740815,  0.024121212808530 },
      {0.823156067318957,  0.028765333012559,  0.023221749533995 },
      {0.707609913381099,  0.144311486950417,  0.048898524516116 },
      {0.527603095742740,  0.324318304588776,  0.063422223181408 },
      {0.324318304588776,  0.527603095742740,  0.063422223181408 },
      {0.144311486950417,  0.707609913381099,  0.048898524516116 },
      {0.028765333012559,  0.823156067318957,  0.023221749533995 },
      {0.937908206225755,  0.032775366614460,  0.012388530705318 },
      {0.806254331245388,  0.164429241594828,  0.026086788660166 },
      {0.601153648467838,  0.369529924372377,  0.033835011360025 },
      {0.369529924372377,  0.601153648467838,  0.033835011360025 },
      {0.164429241594828,  0.806254331245388,  0.026086788660166 },
      {0.032775366614460,  0.937908206225755,  0.012388530705318 }};
    
    a = x1*y2 - x2*y1 + x2*y3 - x3*y2 + x3*y1 - x1*y3;
    area = Math.abs(a)/2.0;
    nn = 36;
    for(int i=0; i<nn; i++) // map onto callers triangle 
    {
      X[i] = x1 + (x2-x1)*tq[i][0] + (x3-x1)*tq[i][1];
      Y[i] = y1 + (y2-y1)*tq[i][0] + (y3-y1)*tq[i][1];
      W[i] = tq[i][2]*area; 
    }
    return nn;
  } // end int6

  void test() // self test on accuracy of tables and mapping
              // nothing should print, BAD means BAD
  {
    double T[] = new double[6];
    double X[] = new double[36];
    double Y[] = new double[36];
    double W[] = new double[36];
    double sum;
    int nv;

    T[0] = 0.0; T[1] = 0.0; // right regular 0,0 0,1 1,0
    T[2] = 0.0; T[3] = 1.0;
    T[4] = 1.0; T[5] = 0.0;

    nv = tint(1, T, X, Y, W);
    sum=0.0;
    for(int i=0; i<nv; i++) sum += W[i]*f1(X[i],Y[i]);
    if(Math.abs(sum-2.0)>1.0e-12)
      System.out.println("BAD triquad order 1, sum="+sum+" not 2.0");

    nv = tint(2, T, X, Y, W);
    sum=0.0;
    for(int i=0; i<nv; i++) sum += W[i]*f2(X[i],Y[i]);
    if(Math.abs(sum-59.0/12.0)>1.0e-12)
      System.out.println("BAD triquad order 2, sum="+sum+" not "+(59.0/12.0));

    nv = tint(3, T, X, Y, W);
    sum=0.0;
    for(int i=0; i<nv; i++) sum += W[i]*f3(X[i],Y[i]);
    if(Math.abs(sum-113.0/12.0)>1.0e-12)
      System.out.println("BAD triquad order 3, sum="+sum+" not "+(113.0/12.0));

    nv = tint(4, T, X, Y, W);
    sum=0.0;
    for(int i=0; i<nv; i++) sum += W[i]*f4(X[i],Y[i]);
    if(Math.abs(sum-1867.0/120.0)>1.0e-11)
      System.out.println("BAD triquad order 4, sum="+sum+
                         " not "+(1867.0/120.0));

    nv = tint(5, T, X, Y, W);
    sum=0.0;
    for(int i=0; i<nv; i++) sum += W[i]*f5(X[i],Y[i]);
    if(Math.abs(sum-2803.0/120.0)>1.0e-10)
      System.out.println("BAD triquad order 5, sum="+sum+
                         " not "+(2803.0/120.0));

    nv = tint(6, T, X, Y, W);
    sum=0.0;
    for(int i=0; i<nv; i++) sum += W[i]*f6(X[i],Y[i]);
    if(Math.abs(sum-82703.0/2520.0)>1.0e-10)
      System.out.println("BAD triquad order 6, sum="+sum+
                         " not "+(82703.0/2520.0));

    T[0] = 2.0; T[1] = 3.0; // obtuse
    T[2] = 2.0; T[3] = 0.0;
    T[4] = 1.0; T[5] = 2.0;

    nv = tint(1, T, X, Y, W);
    sum=0.0;
    for(int i=0; i<nv; i++) sum += W[i]*f1(X[i],Y[i]);
    if(Math.abs(sum-12.0)>1.0e-12)
      System.out.println("BAD triquad b order 1, sum="+sum+
                         " not "+(12.0));

    nv = tint(2, T, X, Y, W);
    sum=0.0;
    for(int i=0; i<nv; i++) sum += W[i]*f2(X[i],Y[i]);
    if(Math.abs(sum-233.0/4.0)>1.0e-12)
      System.out.println("BAD triquad b order 2, sum="+sum+
                         " not "+(233.0/12.0));

    nv = tint(3, T, X, Y, W);
    sum=0.0;
    for(int i=0; i<nv; i++) sum += W[i]*f3(X[i],Y[i]);
    if(Math.abs(sum-4411.0/20.0)>1.0e-11)
      System.out.println("BAD triquad b order 3, sum="+sum+
                         " not "+(4411.0/20.0));

    nv = tint(4, T, X, Y, W);
    sum=0.0;
    for(int i=0; i<nv; i++) sum += W[i]*f4(X[i],Y[i]);
    if(Math.abs(sum-29003.0/40.0)>1.0e-11)
      System.out.println("BAD triquad b order 4, sum="+sum+
                         " not "+(29003.0/40.0));

    nv = tint(5, T, X, Y, W);
    sum=0.0;
    for(int i=0; i<nv; i++) sum += W[i]*f5(X[i],Y[i]);
    if(Math.abs(sum-613741.0/280.0)>1.0e-10)
      System.out.println("BAD triquad b order 5, sum="+sum+
                         " not "+(613741.0/280.0));

    nv = tint(6, T, X, Y, W);
    sum=0.0;
    for(int i=0; i<nv; i++) sum += W[i]*f6(X[i],Y[i]);
    if(Math.abs(sum-252141.0/40.0)>1.0e-9)
      System.out.println("BAD triquad b order 6, sum="+sum+
                         " not "+(252141.0/40.0));

  } // end test

  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;
  }

} // end class triquad