/* tri_basis.c  from Celia and Gray p 415, 416 */
/* triangle T defined by:
 * T[0] is x1, T[1] is y1, T[2] is x2, T[3] is y2, T[4] is x3, T[5] is y3 
 * 
 *            * x3,y3       node 4, midpoint  1 and 2
 *           / \            node 5, midpoint  2 and 3
 *          / _ * x2,y2     node 6, midpoint  1 and 3
 *   x1,y1 *
 *
 * phi1(x,y) = (1/2A)*((x2*y3-x3*y2) + (y2-y3)*x + (x3-x2)*y)
 * phi2(x,y) = (1/2A)*((x3*y1-x1*y3) + (y3-y1)*x + (x1-x3)*y)
 * phi3(x,y) = (1/2A)*((x1*y2-x2*y1) + (y1-y2)*x + (x2-x1)*y)
 * A = (1/2)*((x2*y3-x3*y2)+(x3*y1-x1*y3)+(x1*y2-x2*y1))
 * 2A =       (t1          + t2          + t3)
 *
 * phi21(x,y) = phi1(x,y)*(2*phi1(x,y)-1)
 * phi22(x,y) = phi2(x,y)*(2*phi2(x,y)-1)
 * phi23(x,y) = phi3(x,y)*(2*phi3(x,y)-1)
 * phi24(x,y) = 4*phi1(x,y)*phi2(x,y)
 * phi25(x,y) = 4*phi2(x,y)*phi3(x,y)
 * phi26(x,y) = 4*phi3(x,y)*phi1(x,y)
 */

#include "tri_basis.h"

/* first order, three points, adjust order in T */
double phi(double T[], double x, double y)
{
  return phi1(T,x,y);
}

double phix(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;

  return (y2-y3)/A2;
}


double phiy(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;

  return (x3-x2)/A2;
}


double phi1(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;

  return (t1 + (y2-y3)*x + (x3-x2)*y)/A2;
}

double phi2(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;

  return (t2 + (y3-y1)*x + (x1-x3)*y)/A2;
}

double phi3(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;

  return (t3 + (y1-y2)*x + (x2-x1)*y)/A2;
}


/* second order, about first point in T */
double phi2p(double T[], double x, double y)
{
  return phi21(T,x,y);
}

double phi2x(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;
  double a  = t1;
  double b  = y2-y3;
  double c  = x3-x2;
  return ((4.0*a*b + 4.0*b*b*x + 4.0*b*c*y)/(A2*A2))-b/A2;
}

double phi2y(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;
  double a  = t1;
  double b  = y2-y3;
  double c  = x3-x2;
  return ((4.0*a*c + 4.0*b*c*x + 4.0*c*c*y)/(A2*A2))-c/A2; 
}

  double phi2xx(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;
  double a  = t1;
  double b  = y2-y3;
  double c  = x3-x2;
  return 4.0*b*b/(A2*A2);
}

double phi2xy(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;
  double a  = t1;
  double b  = y2-y3;
  double c  = x3-x2;
  return 4.0*b*c/(A2*A2);
}

double phi2yy(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;
  double a  = t1;
  double b  = y2-y3;
  double c  = x3-x2;
  return 4.0*c*c/(A2*A2);
}


/* second order, about mid point between first and second point in T */
double phi2m(double T[], double x, double y)
{
  return phi24(T,x,y);
}

double phi2mx(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;
  double a  = t1;
  double b  = y2-y3;
  double c  = x3-x2;
  double d  = t2;
  double e  = y3-y1;
  double f  = x1-x3;
  return 4.0*(b*(d+e*x+f*y)+e*(a+b*x+c*y))/(A2*A2);
}

double phi2my(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;
  double a  = t1;
  double b  = y2-y3;
  double c  = x3-x2;
  double d  = t2;
  double e  = y3-y1;
  double f  = x1-x3;
  return 4.0*(c*(d+e*x+f*y)+f*(a+b*x+c*y))/(A2*A2);
}

double phi2mxx(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;
  double a  = t1;
  double b  = y2-y3;
  double c  = x3-x2;
  double d  = t2;
  double e  = y3-y1;
  double f  = x1-x3;
  return 8.0*b*e/(A2*A2);
}

double phi2mxy(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;
  double a  = t1;
  double b  = y2-y3;
  double c  = x3-x2;
  double d  = t2;
  double e  = y3-y1;
  double f  = x1-x3;
  return 4.0*(b*f+c*e)/(A2*A2);
}

double phi2myy(double T[], double x, double y)
{
  double x1 = T[0];
  double y1 = T[1];
  double x2 = T[2];
  double y2 = T[3];
  double x3 = T[4];
  double y3 = T[5];
  double t1 = x2*y3 - x3*y2;
  double t2 = x3*y1 - x1*y3;
  double t3 = x1*y2 - x2*y1;
  double A2 = t1+t2+t3;
  double a  = t1;
  double b  = y2-y3;
  double c  = x3-x2;
  double d  = t2;
  double e  = y3-y1;
  double f  = x1-x3;
  return 8.0*c*f/(A2*A2);
}


/* second order, same T for all, second number is point */
double phi21(double T[], double x, double y)
{
  return phi1(T,x,y)*(2.0*phi1(T,x,y)-1.0);
}

double phi22(double T[], double x, double y)
{
  return phi2(T,x,y)*(2.0*phi2(T,x,y)-1.0);
}

double phi23(double T[], double x, double y)
{
  return phi3(T,x,y)*(2.0*phi3(T,x,y)-1.0);
}

double phi24(double T[], double x, double y)
{
  return 4.0*phi1(T,x,y)*phi2(T,x,y);
}

double phi25(double T[], double x, double y)
{
  return 4.0*phi2(T,x,y)*phi3(T,x,y);
}

double phi26(double T[], double x, double y)
{
  return 4.0*phi3(T,x,y)*phi1(T,x,y);
}

/* end tri_basis.c */