% phi_tril.m    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
%  T may be any of x1,y1,x2,y2,x3,y3  x2,y2,x3,y3,x1,y1  x3,y3,x1,y1,x2,y2
%
function phi_tril()
  clear
  format compact
  disp('phi_tril.m running')
  disp('Lagrange fit from all three points and derivatives')
  x1=1.0; % triangle
  y1=1.0;
  x2=1.0;
  y2=2.0;
  x3=2.0;
  y3=1.0;
  xmin=min(min(x1,x2),x3);
  xmax=max(max(x1,x2),x3);
  ymin=min(min(y1,y2),y3);
  ymax=max(max(y1,y2),y3);
  iter=21; % number of lines
  offs=0.7/iter; % mark vertex
  X1=[x1,x1+offs]; % markers, keep inside for better plots
  Y1=[y1,y1+offs];
  X2=[x2,x2+offs];
  Y2=[y2,y2-offs];
  X3=[x3,x3-offs];
  Y3=[y3,y3+offs];
  
  T=[x1,y1, x2,y2, x3,y3]
  disp('should be 1, 0, 0')
  disp(phi22l(T,T(1),T(2)))
  disp(phi22l(T,T(3),T(4)))
  disp(phi22l(T,T(5),T(6)))
  figure(1)
  for i=1:iter
    X(i)=xmin+(i-1)*(xmax-xmin)/(iter-1.0);
    for j=1:iter
      Y(j)=ymin+(j-1)*(ymax-ymin)/(iter-1.0);
      Z(i,j)=phi22l(T,X(i),Y(j));
    end
  end
  mesh(Y,X,Z)
  hold on

  for i=1:2 % markers
    for j=1:2
      Z1(i,j)=phi22l(T,X1(i),Y1(j));
    end
  end
  mesh(Y1,X1,Z1)
  for i=1:2
    for j=1:2
      Z2(i,j)=phi22l(T,X2(i),Y2(j));
    end
  end
  mesh(Y2,X2,Z2)
  for i=1:2
    for j=1:2
      Z3(i,j)=phi22l(T,X3(i),Y3(j));
    end
  end
  mesh(Y3,X3,Z3)
  title('Lagrange fit about point 1')
  xlabel('X')
  ylabel('Y')
  hold off
  
  disp('point 2')
  figure(2)
  T=[x2,y2, x1,y1, x3,y3]
  disp('should be 1, 0, 0')
  disp(phi22l(T,T(1),T(2)))
  disp(phi22l(T,T(3),T(4)))
  disp(phi22l(T,T(5),T(6)))
  for i=1:iter
    X(i)=xmin+(i-1)*(xmax-xmin)/(iter-1.0);
    for j=1:iter
      Y(j)=ymin+(j-1)*(ymax-ymin)/(iter-1.0);
      Z(i,j)=phi22l(T,X(i),Y(j));
    end
  end
  mesh(Y,X,Z)
  hold on
  for i=1:2  % markers
    for j=1:2
      Z1(i,j)=phi22l(T,X1(i),Y1(j));
    end
  end
  mesh(Y1,X1,Z1)
  for i=1:2
    for j=1:2
      Z2(i,j)=phi22l(T,X2(i),Y2(j));
    end
  end
  mesh(Y2,X2,Z2)
  for i=1:2
    for j=1:2
      Z3(i,j)=phi22l(T,X3(i),Y3(j));
    end
  end
  mesh(Y3,X3,Z3)
  title('Lagrange fit about point 2')
  xlabel('X')
  ylabel('Y')
  hold off
  
  disp('point 3')
  figure(3)
  T=[x3,y3, x2,y2, x1,y1]
  disp('should be 1, 0, 0')
  disp(phi22l(T,T(1),T(2)))
  disp(phi22l(T,T(3),T(4)))
  disp(phi22l(T,T(5),T(6)))
  for i=1:iter
    X(i)=xmin+(i-1)*(xmax-xmin)/(iter-1.0);
    for j=1:iter
      Y(j)=ymin+(j-1)*(ymax-ymin)/(iter-1.0);
      Z(i,j)=phi22l(T,X(i),Y(j));
    end
  end
  mesh(Y,X,Z)
  hold on
  for i=1:2 % markers
    for j=1:2
      Z1(i,j)=phi22l(T,X1(i),Y1(j));
    end
  end
  mesh(Y1,X1,Z1)
  for i=1:2
    for j=1:2
      Z2(i,j)=phi22l(T,X2(i),Y2(j));
    end
  end
  mesh(Y2,X2,Z2)
  for i=1:2
    for j=1:2
      Z3(i,j)=phi22l(T,X3(i),Y3(j));
    end
  end
  mesh(Y3,X3,Z3)
  title('Lagrange fit about point 3')
  xlabel('X')
  ylabel('Y')
  hold off

  figure(4)
  for i=1:iter
    X(i)=xmin+(i-1)*(xmax-xmin)/(iter-1.0);
    for j=1:iter
      Y(j)=ymin+(j-1)*(ymax-ymin)/(iter-1.0);
      Z(i,j)=phi22lx(T,X(i),Y(j));
    end
  end
  mesh(Y,X,Z)
  hold on
  for i=1:2 % markers
    for j=1:2
      Z1(i,j)=phi22lx(T,X1(i),Y1(j));
    end
  end
  mesh(Y1,X1,Z1)
  for i=1:2
    for j=1:2
      Z2(i,j)=phi22lx(T,X2(i),Y2(j));
    end
  end
  mesh(Y2,X2,Z2)
  for i=1:2
    for j=1:2
      Z3(i,j)=phi22lx(T,X3(i),Y3(j));
    end
  end
  mesh(Y3,X3,Z3)
  title('derivative dx about point 3')
  xlabel('X')
  ylabel('Y')
  hold off

  
  figure(5)
  for i=1:iter
    X(i)=xmin+(i-1)*(xmax-xmin)/(iter-1.0);
    for j=1:iter
      Y(j)=ymin+(j-1)*(ymax-ymin)/(iter-1.0);
      Z(i,j)=phi22ly(T,X(i),Y(j));
    end
  end
  mesh(Y,X,Z)
  hold on
  for i=1:2 % markers
    for j=1:2
      Z1(i,j)=phi22ly(T,X1(i),Y1(j));
    end
  end
  mesh(Y1,X1,Z1)
  for i=1:2
    for j=1:2
      Z2(i,j)=phi22ly(T,X2(i),Y2(j));
    end
  end
  mesh(Y2,X2,Z2)
  for i=1:2
    for j=1:2
      Z3(i,j)=phi22ly(T,X3(i),Y3(j));
    end
  end
  mesh(Y3,X3,Z3)
  title('derivative dy about point 3')
  xlabel('X')
  ylabel('Y')
  hold off

  
  figure(6)
  for i=1:iter
    X(i)=xmin+(i-1)*(xmax-xmin)/(iter-1.0);
    for j=1:iter
      Y(j)=ymin+(j-1)*(ymax-ymin)/(iter-1.0);
      Z(i,j)=phi22lxx(T,X(i),Y(j));
    end
  end
  mesh(Y,X,Z)
  hold on
  for i=1:2 % markers
    for j=1:2
      Z1(i,j)=phi22lxx(T,X1(i),Y1(j));
    end
  end
  mesh(Y1,X1,Z1)
  for i=1:2
    for j=1:2
      Z2(i,j)=phi22lxx(T,X2(i),Y2(j));
    end
  end
  mesh(Y2,X2,Z2)
  for i=1:2
    for j=1:2
      Z3(i,j)=phi22lxx(T,X3(i),Y3(j));
    end
  end
  mesh(Y3,X3,Z3)
  title('derivative dxx about point 3')
  xlabel('X')
  ylabel('Y')
  hold off


  figure(7)
  for i=1:iter
    X(i)=xmin+(i-1)*(xmax-xmin)/(iter-1.0);
    for j=1:iter
      Y(j)=ymin+(j-1)*(ymax-ymin)/(iter-1.0);
      Z(i,j)=phi22lyy(T,X(i),Y(j));
    end
  end
  mesh(Y,X,Z)
  hold on
  for i=1:2 % markers
    for j=1:2
      Z1(i,j)=phi22lyy(T,X1(i),Y1(j));
    end
  end
  mesh(Y1,X1,Z1)
  for i=1:2
    for j=1:2
      Z2(i,j)=phi22lyy(T,X2(i),Y2(j));
    end
  end
  mesh(Y2,X2,Z2)
  for i=1:2
    for j=1:2
      Z3(i,j)=phi22lyy(T,X3(i),Y3(j));
    end
  end
  mesh(Y3,X3,Z3)
  title('derivative dyy about point 3')
  xlabel('X')
  ylabel('Y')
  hold off



  disp('done')
  
  % use product of  lterm  calls at various  x0,y0=1.0.xi,yi=0.0  points
  function l=lterm(x,y,xi,yi,x0,y0)
    l=(((x-xi)*(x-xi)+(y-yi)*(y-yi))/((x0-xi)*(x0-xi)+(y0-yi)*(y0-yi)));
  end    
    
  function  d=dlx(x,y,xi,yi,x0,y0)
    d=(2.0*(x-xi)/((x0-xi)*(x0-xi)+(y0-yi)*(y0-yi)));
  end
  
  function  d=dly(x,y,xi,yi,x0,y0)
    d=(2.0*(y-yi)/((x0-xi)*(x0-xi)+(y0-yi)*(y0-yi)));
  end
  
  function  d=dlxx(x,y,xi,yi,x0,y0)
    d=(2.0/((x0-xi)*(x0-xi)+(y0-yi)*(y0-yi)));
  end
  
  
  function  d=dlyy(x,y,xi,yi,x0,y0)
    d=(2.0/((x0-xi)*(x0-xi)+(y0-yi)*(y0-yi)));
  end
  
  function  d=dlxy(x,y,xi,yi,x0,y0)
    d=0.0;
  end
  
  function  val=phi22l(T, x, y) % about T[0],T[1]
    x1=T(1);
    y1=T(2);
    x2=T(3);
    y2=T(4);
    x3=T(5);
    y3=T(6);
    val = lterm(x,y,x2,y2,x1,y1) * lterm(x,y,x3,y3,x1,y1);
    % val=1.0 at x1,y1  val=0.0 at x2,y2  and  x3,y3
  end

  function val=phi22lx(T, x, y)
    x1=T(1);
    y1=T(2);
    x2=T(3);
    y2=T(4);
    x3=T(5);
    y3=T(6);
    val = dlx(x,y,x2,y2,x1,y1) * lterm(x,y,x3,y3,x1,y1) + ...
          dlx(x,y,x3,y3,x1,y1) * lterm(x,y,x2,y2,x1,y1);
  end

  function val=phi22ly(T, x, y)
    x1=T(1);
    y1=T(2);
    x2=T(3);
    y2=T(4);
    x3=T(5);
    y3=T(6);
    val = dly(x,y,x2,y2,x1,y1) * lterm(x,y,x3,y3,x1,y1) + ...
          dly(x,y,x3,y3,x1,y1) * lterm(x,y,x2,y2,x1,y1);
  end

  function val=phi22lxx(T, x, y)
    x1=T(1);
    y1=T(2);
    x2=T(3);
    y2=T(4);
    x3=T(5);
    y3=T(6);
    val =dlxx(x,y,x2,y2,x1,y1) * lterm(x,y,x3,y3,x1,y1) + ...
         dlx(x,y,x3,y3,x1,y1) * dlx(x,y,x2,y2,x1,y1) + ...
         dlxx(x,y,x3,y3,x1,y1) * lterm(x,y,x2,y2,x1,y1) + ...
         dlx(x,y,x2,y2,x1,y1) * dlx(x,y,x3,y3,x1,y1);
  end

  function val=phi22lxy(T, x, y)
    val=0.0;
  end

  function val=phi22lyy(T, x, y)
    x1=T(1);
    y1=T(2);
    x2=T(3);
    y2=T(4);
    x3=T(5);
    y3=T(6);
    val =dlyy(x,y,x2,y2,x1,y1) * lterm(x,y,x3,y3,x1,y1) + ...
         dly(x,y,x3,y3,x1,y1) * dly(x,y,x2,y2,x1,y1) + ...
         dlyy(x,y,x3,y3,x1,y1) * lterm(x,y,x2,y2,x1,y1) + ...
         dly(x,y,x2,y2,x1,y1) * dly(x,y,x3,y3,x1,y1);
  end
  
end % phi_tril.m