// phiquad2D.java  Lagrange Phi functions for quad, four points 
//         function prototype naming convention:   
//           phi            basic Lagrange polynomial in x,y       
//           phix           first derivative with respect to x
//           phiy           first derivative with respect to y
//           phixx          second derivative with respect to x
//           phiyy          second derivative with respect to y
//           phixy          derivative with respect to x, then y
// 
//    Q[] is quad vertices  x1,y1,x2,y2,x3,y3,x4,y4 in order
//    the functions have phi(x1,y1)=1,
//                       phi(x2,y2)=0, phi(x3,y3)=0, phi(x4,y4)=0
// example  phi(Q, x, y) is                            
//        ((x -x2)^2+(y -y2)^2)*((x -x3)^2+(y -y3)^2)*((x -x4)^2+(y -y4)^2)/
//        ((x1-x2)^2+(y1-y2)^2)*((x1-x3)^2+(y1-y3)^2)*((x1-x4)^2+(y1-y4)^2)

class phiquad2D
{
  phiquad2D()
  {
    System.out.println("phiquad2D instantiated");
  }

  // phi 
  double phi(double Q[], double x, double y)
  {
      double x1, y1, x2, y2, x3, y3, x4, y4, v, d;
    x1 = Q[0]; y1 = Q[1];
    x2 = Q[2]; y2 = Q[3];
    x3 = Q[4]; y3 = Q[5];
    x4 = Q[6]; y4 = Q[7];
    d = (sq(x1-x2)+sq(y1-y2))*(sq(x1-x3)+sq(y1-y3))*(sq(x1-x4)+sq(y1-y4));
    v = (sq(x -x2)+sq(y -y2))*(sq(x -x3)+sq(y -y3))*(sq(x -x4)+sq(y -y4));
    return v/d;
  }

  // phix 
  double phix(double Q[], double x, double y)
  {
    double x1, y1, x2, y2, x3, y3, x4, y4;
    double a2, a3, a4, v, d;
    x1 = Q[0]; y1 = Q[1];
    x2 = Q[2]; y2 = Q[3];
    x3 = Q[4]; y3 = Q[5];
    x4 = Q[6]; y4 = Q[7];
    d = (sq(x1-x2)+sq(y1-y2))*(sq(x1-x3)+sq(y1-y3))*(sq(x1-x4)+sq(y1-y4));
    a2 = (sq(x -x2)+sq(y -y2));
    a3 = (sq(x -x3)+sq(y -y3));
    a4 = (sq(x -x4)+sq(y -y4));
    v = 2.0*(x-x2)*a3*a4+
        2.0*(x-x3)*a2*a4+
        2.0*(x-x4)*a2*a3;
    return v/d;
  }

  // phiy 
  double phiy(double Q[], double x, double y)
  {
    double x1, y1, x2, y2, x3, y3, x4, y4;
    double a2, a3, a4, v, d;
    x1 = Q[0]; y1 = Q[1];
    x2 = Q[2]; y2 = Q[3];
    x3 = Q[4]; y3 = Q[5];
    x4 = Q[6]; y4 = Q[7];
    d = (sq(x1-x2)+sq(y1-y2))*(sq(x1-x3)+sq(y1-y3))*(sq(x1-x4)+sq(y1-y4));
    a2 = (sq(x -x2)+sq(y -y2));
    a3 = (sq(x -x3)+sq(y -y3));
    a4 = (sq(x -x4)+sq(y -y4));
    v = 2.0*(y-y2)*a3*a4+
        2.0*(y-y3)*a2*a4+
        2.0*(y-y4)*a2*a3;
    return v/d;
  }

  // phixx 
  double phixx(double Q[], double x, double y)
  {
    double x1, y1, x2, y2, x3, y3, x4, y4;
    double a2, a3, a4, v, d;
    x1 = Q[0]; y1 = Q[1];
    x2 = Q[2]; y2 = Q[3];
    x3 = Q[4]; y3 = Q[5];
    x4 = Q[6]; y4 = Q[7];
    d = (sq(x1-x2)+sq(y1-y2))*(sq(x1-x3)+sq(y1-y3))*(sq(x1-x4)+sq(y1-y4));
    a2 = (sq(x -x2)+sq(y -y2));
    a3 = (sq(x -x3)+sq(y -y3));
    a4 = (sq(x -x4)+sq(y -y4));
    v = 2.0*a3*a4+
        2.0*a2*a4+
        2.0*a2*a3+
        8.0*(x-x2)*(x-x3)*a4+
        8.0*(x-x2)*(x-x4)*a3+
        8.0*(x-x3)*(x-x4)*a2;
    return v/d;
  }

  // phixy 
  double phixy(double Q[], double x, double y)
  {
    double x1, y1, x2, y2, x3, y3, x4, y4;
    double a2, a3, a4, v, d;
    x1 = Q[0]; y1 = Q[1];
    x2 = Q[2]; y2 = Q[3];
    x3 = Q[4]; y3 = Q[5];
    x4 = Q[6]; y4 = Q[7];
    d = (sq(x1-x2)+sq(y1-y2))*(sq(x1-x3)+sq(y1-y3))*(sq(x1-x4)+sq(y1-y4));
    a2 = (sq(x -x2)+sq(y -y2));
    a3 = (sq(x -x3)+sq(y -y3));
    a4 = (sq(x -x4)+sq(y -y4));
    v = 4.0*(x -x2)*(y -y3)*a4+
        4.0*(x -x2)*(y -y4)*a3+
        4.0*(x -x3)*(y -y2)*a4+
        4.0*(x -x3)*(y -y4)*a2+
        4.0*(x -x4)*(y -y2)*a3+
        4.0*(x -x4)*(y -y3)*a2;
    return v/d;
  }

  // phiyy 
  double phiyy(double Q[], double x, double y)
  {
    double x1, y1, x2, y2, x3, y3, x4, y4;
    double a2, a3, a4, v, d;
    x1 = Q[0]; y1 = Q[1];
    x2 = Q[2]; y2 = Q[3];
    x3 = Q[4]; y3 = Q[5];
    x4 = Q[6]; y4 = Q[7];
    d = (sq(x1-x2)+sq(y1-y2))*(sq(x1-x3)+sq(y1-y3))*(sq(x1-x4)+sq(y1-y4));
    a2 = (sq(x -x2)+sq(y -y2));
    a3 = (sq(x -x3)+sq(y -y3));
    a4 = (sq(x -x4)+sq(y -y4));
    v = 2.0*a3*a4+
        2.0*a2*a4+
        2.0*a2*a3+
        8.0*(y-y2)*(y-y3)*a4+
        8.0*(y-y2)*(y-y4)*a3+
	8.0*(y-y3)*(y-y4)*a2;
    return v/d;
  }

  double sq(double a)
  {
    return a*a;
  }
      
} // end phiquad2D.java