%---
close all
clear all
hold on
%---

%=========
% interpolation in a rectangular cavity
%
% confined in ax<x<bx
%             ay<y<by
%=========

ax=-1.3; bx=1.3;
ay=-1.0; by=1.0;
Nx=16; Ny=8;
Nx=16; Ny=18;

%-----
% grid generation
%-----

[glx,gly,gx,gy] = grid_2d (ax,bx,ay,by,Nx,Ny);

%---
% specify the grid velocity
%---

 for i=1:Nx+1
   for j=1:Ny+1
     gux(i,j) = gx(i,j);
     guy(i,j) =-gy(i,j);
     px = gx(i,j);
     py = gy(i,j);
     wnx  = 2*pi/(bx-ax);
     wny  = 2*pi/(by-ay);
     gux(i,j) = cos(wnx*px)*cos(wny*py);
     guy(i,j) = sin(wnx*px)*sin(wny*py);
   end
 end

%---------------------------
% velocity vector field
% in the middle of each cell
%---------------------------

 for i=1:Nx
   x=0.5*(glx(i)+glx(i+1));

   for j=1:Ny
   y=0.5*(gly(j)+gly(j+1));

   plot(x,y,'g.')

   [ux,uy] = rec_2d_int ...
    ...
      (ax,bx ... % x end points
      ,ay,by ... % y end points
      ,Nx,Ny ... % grid size
      ,glx,gly ... % grid lines
      ,gux,guy ... % grid velocity
      ,x,y ...      % interpolation point
      );

     vector = draw_arrow_2d ...
     ...
     (x,y,ux/Nx,uy/Nx);

   plot(vector(:,1),vector(:,2));

  end
 end

%---
% plot
%----

plot([ax, bx, bx, ax, ax],[ay, ay, by, by, ay])
axis([ax-0.1, bx+0.1, ay-0.1, by+0.1])
axis equal
xlabel('x','fontsize',15)
ylabel('y','fontsize',15)
set(gca,'fontsize',15)
box

%---
% velocity gradient tensor
%----

 [Axx,Axy,Ayx,Ayy] = rec_2d_vgt ...
...
  (glx,gly,Nx,Ny,gux,guy);


%---
% compute the rate of expansion
%         the rate of strain tensor 
%         the strains
%         the vorticity
%---

 for i=1:Nx+1
  for j=1:Ny+1

    roe(i,j) = Axx(i,j)+Ayy(i,j);    % rate of expansion
   omega(i,j) = Axy(i,j)-Ayx(i,j);   % vorticity

   Exx(i,j) = Axx(i,j)-0.5*roe(i,j);     % rate of deformation tensor
   Exy(i,j) = 0.5*(Axy(i,j)+Ayx(i,j));   % rate of deformation tensor
   Eyx(i,j) = Exy(i,j);                  % rate of deformation tensor
   Eyy(i,j) = Ayy(i,j)-0.5*roe(i,j);     % rate of deformation tensor

   det = 4.0*(Exx(i,j)^2+Exy(i,j)^2);   % eigenvalues of rate of strain
   srd = sqrt(det);
   strain1(i,j) =  0.5*srd;
   strain2(i,j) = -0.5*srd;

%---
% compute the eigenvectors of the rate of strain
%---

if(abs(Exy(i,j))<0.0001) % E is diagonal

  if(abs(Exx(i,j)-strain1(i,j))>0.0001)
     egv1_x(i,j) = 0.0;
     egv1_y(i,j) = 1.0;
  else
     egv1_x(i,j) = 1.0;
     egv1_y(i,j) = 0.0;
  end
  if(abs(Exx(i,j)-strain2(i,j))>0.0001) 
     egv2_x(i,j) = 0.0;
     egv2_y(i,j) = 1.0;
  else
     egv2_x(i,j) = 1.0;
     egv2_y(i,j) = 0.0;
  end

 else  % E is not diagonal

  egv1_x(i,j) = 1.0;
  egv1_y(i,j) = -(Exx(i,j)-strain1(i,j))/Exy(i,j);
  egv2_x(i,j) = 1.0;
  egv2_y(i,j) = -(Exx(i,j)-strain2(i,j))/Exy(i,j);

 end

%---
% normalize eigenvectors
%---

  fc1 = 1/sqrt(egv1_x(i,j)^2+egv1_y(i,j)^2);
  egv1_x(i,j) = fc1*egv1_x(i,j);
  egv1_y(i,j) = fc1*egv1_y(i,j);
  fc2 = 1/sqrt(egv2_x(i,j)^2+egv2_y(i,j)^2);
  egv2_x(i,j) = fc2*egv2_x(i,j);
  egv2_y(i,j) = fc2*egv2_y(i,j);

 end
end

figure
mesh(glx,gly,omega')
xlabel('x','fontsize',15)
ylabel('y','fontsize',15)
zlabel('\omega','fontsize',15)
set(gca,'fontsize',15)
box

figure
mesh(glx,gly,roe')
xlabel('x','fontsize',15)
ylabel('y','fontsize',15)
zlabel('\alpha','fontsize',15)
set(gca,'fontsize',15)
box

figure
mesh(glx,gly,strain1')
xlabel('x','fontsize',15)
ylabel('y','fontsize',15)
zlabel('s_1','fontsize',15)
set(gca,'fontsize',15)
box

figure
mesh(glx,gly,strain2')
xlabel('x','fontsize',15)
ylabel('y','fontsize',15)
zlabel('s_2','fontsize',15)
set(gca,'fontsize',15)
box

%---
% eigenvector 1
%---

figure
hold on
[glx,gly,gx,gy] = grid_2d (ax,bx,ay,by,Nx,Ny);
for i=1:Nx+1
 for j=1:Ny+1
     vector = draw_arrow_2d ...
     ...
     (gx(i,j),gy(i,j),egv1_x(i,j)/Ny,egv1_y(i,j)/Ny);
   plot(vector(:,1),vector(:,2));
  end
end
xlabel('x','fontsize',15)
ylabel('y','fontsize',15)
set(gca,'fontsize',15)
box

%---
% eigenvector 2
%---

figure
hold on
[glx,gly,gx,gy] = grid_2d (ax,bx,ay,by,Nx,Ny);
for i=1:Nx+1
 for j=1:Ny+1
     vector = draw_arrow_2d ...
     ...
     (gx(i,j),gy(i,j),egv2_x(i,j)/Ny,egv2_y(i,j)/Ny);
   plot(vector(:,1),vector(:,2));
  end
end
xlabel('x','fontsize',15)
ylabel('y','fontsize',15)
set(gca,'fontsize',15)
box