/* phi_tril.c    Phi and derivatives for order 1 through 4 for Galerkin FEM
 *               and dimension 2 (x,y) on triangles
 *  phi22l(T,x,y) is second order FEM phi function in 2 dimensions
 *         =sum i=1:3 j=mod i+1 ((x -xi)^2+(y -yi)^2)/((xj-xi)^2+(yj-yi)^2)
 *  phi22lx(T,x,y) is partial derivative of phi22l(T,x,y) with respect to x
 *  phi22ly(T,x,y) is partial derivative of phi22l(T,x,y) with respect to y
 *  phi22lxx(T,x,y) is second partial of phi22l(T,x,y) with respect to x
 *  phi22lxy(T,x,y) is second partial of phi22l(T,x,y) wrt x and wrt y
 *  phi22lyy(T,x,y) is second partial of phi22l(T,x,y) with respect to y
 *
 *  ...  up to fourth derivative of fourth order phi in two dimensions
 *  phi24l(T,x,y) is fourth order four dimensional FEM phi function
 *  ...
 *  phi24lyyyy(T,x,y) is fourth partial of phi24l(T,x,y) wrt y
 *
 *  T man be any of x1,y1,x2,y2,x3,y3  x2,y2,x3,y3,x1,y1  x3,y3,x1,y1,x2,y2
 */

/* use product of  lterm  calls at various  x0,y0=1.0.xi,yi=0.0  points */
#define lterm(x,y,xi,yi,x0,y0) (((x-xi)*(x-xi)+(y-yi)*(y-yi))/((x0-xi)*(x0-xi)+(y0-yi)*(y0-yi)))
#define dlx(x,y,xi,yi,x0,y0) (2.0*(x-xi)/((x0-xi)*(x0-xi)+(y0-yi)*(y0-yi)))
#define dly(x,y,xi,yi,x0,y0) (2.0*(y-yi)/((x0-xi)*(x0-xi)+(y0-yi)*(y0-yi)))
#define dlxx(x,y,xi,yi,x0,y0) (2.0/((x0-xi)*(x0-xi)+(y0-yi)*(y0-yi)))
#define dlyy(x,y,xi,yi,x0,y0) (2.0/((x0-xi)*(x0-xi)+(y0-yi)*(y0-yi)))
#define dlxy(x,y,xi,yi,x0,y0) (0.0)

double phi22l(double T[], double x, double y) /* about T[0],T[1] */
{
  double x1, y1, x2, y2, x3, y3;
  double val;

  x1=T[0]; y1=T[1]; x2=T[2]; y2=T[3]; x3=T[4]; y3=T[5]; 
  val = lterm(x,y,x2,y2,x1,y1) * lterm(x,y,x3,y3,x1,y1);
  return val; /* val=1.0 at x1,y1  val=0.0 at x3,y2  and  x3,y3 */
}

double phi22lx(double T[], double x, double y)
{
  double x1, y1, x2, y2, x3, y3;
  double val;

  x1=T[0]; y1=T[1]; x2=T[2]; y2=T[3]; x3=T[4]; y3=T[5]; 
  val = dlx(x,y,x2,y2,x1,y1) * lterm(x,y,x3,y3,x1,y1);
  val+= dlx(x,y,x3,y3,x1,y1) * lterm(x,y,x2,y2,x1,y1);
  return val;
}

double phi22ly(double T[], double x, double y)
{
  double x1, y1, x2, y2, x3, y3;
  double val;

  x1=T[0]; y1=T[1]; x2=T[2]; y2=T[3]; x3=T[4]; y3=T[5]; 
  val =dly(x,y,x2,y2,x1,y1) * lterm(x,y,x3,y3,x1,y1);
  val+=dly(x,y,x3,y3,x1,y1) * lterm(x,y,x2,y2,x1,y1);
  return val;
}

double phi22lxx(double T[], double x, double y)
{
  double x1, y1, x2, y2, x3, y3;
  double val;

  x1=T[0]; y1=T[1]; x2=T[2]; y2=T[3]; x3=T[4]; y3=T[5]; 
  val =dlxx(x,y,x2,y2,x1,y1) * lterm(x,y,x3,y3,x1,y1);
  val+=dlx(x,y,x3,y3,x1,y1) * dlx(x,y,x2,y2,x1,y1);
  val+=dlxx(x,y,x3,y3,x1,y1) * lterm(x,y,x2,y2,x1,y1);
  val+=dlx(x,y,x2,y2,x1,y1) * dlx(x,y,x3,y3,x1,y1);
  return val;
}

double phi22lxy(double T[], double x, double y)
{
  double x1, y1, x2, y2, x3, y3;
  double u, val;

  x1=T[0]; y1=T[1]; x2=T[2]; y2=T[3]; x3=T[4]; y3=T[5]; 
  u = ((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))*((x1-x3)*(x1-x3)+(y1-y3)*(y1-y3));
  val = 4.0*(2.0*x*y-x*y3-x2*y+x2*y3-x3*y-x*y2+x3*y2)/u;

  return val;
}

double phi22lyy(double T[], double x, double y)
{
  double x1, y1, x2, y2, x3, y3;
  double val=0.0;

  x1=T[0]; y1=T[1]; x2=T[2]; y2=T[3]; x3=T[4]; y3=T[5]; 
  val =dlyy(x,y,x2,y2,x1,y1) * lterm(x,y,x3,y3,x1,y1);
  val+=dly(x,y,x3,y3,x1,y1) * dly(x,y,x2,y2,x1,y1);
  val+=dlyy(x,y,x3,y3,x1,y1) * lterm(x,y,x2,y2,x1,y1);
  val+=dly(x,y,x2,y2,x1,y1) * dly(x,y,x3,y3,x1,y1);
  return val;
}

/* end phi_tril.c */