% EB_Beam_Diffop_Example.m  
%
% example code for flow stability projects using the method of discretizing
% differential operators
%
% in this code the 1D Laplace eigenvalue problem is solved:
% I*E/rho/A * d^4 u(x)/dx^4 = omega^2 u(x)
% a <= x <= b
% u(a) = u(b) = 0
%
% IMPORTANT REMARK 
% In the case of the flow stability problems, you might have to assemble the operator from multiple matrices
% or also, might have to assemble multiple operators (each on which the
% boundary consitions must be enforced), and find the eigenvalues using the
% command <polyeig> 

%%
clear all
close all
clc

%% parameters
% Chebyshev polynomial degree - the number of points is actually N+1
N = 41;
N = 81;
N = 151;
% domain boundaries
a = 0;
b = 1;

% material properties
ly = 1e-2; % heigth of the beam
lz = 1e-2; % transversal size of the beam
A = ly*lz; % area
I = lz*ly^3/12; % second moment of area
rho = 7800; % density
E = 200*1e9; % Elasticity modulus

matpar = I*E/rho/A; % combining parameters into a single constant

%% chebyshev differentiation matrices
% chebyshev differentiation matrix and grid, on [-1,1]
[Dcheb,xcheb] = cheb(N);
% mapping chebyshev grid to the computational domain [a,b]
% then mapping differentiation matrix
% algebraic mapping
x=(xcheb+1)*(b-a)/2 + a;
D1 = Dcheb*(2/(b-a));

% rational mapping, in the case the point need to be refined toward one of
% the boundaries
% half the points lie on the interval [a, yhalf], and the other half on [yhalf,b];
% yhalf = 0.2;
% c1 = b*yhalf/(b-2*yhalf);
% c2 = 1 + 2*c1/(b-a);
% y = c1*(1+xcheb)./(c2-xcheb) + a;
% dxi = (c2-xcheb).^2 / (c1*(c2+1));
% D1 = diag(dxi) * Dcheb;

% differentiation matrix for the higher derivatives
D2 = D1^2;
D3 = D1^3;
D4 = D1^4;

%% boundary conditions and eigenvalue calculation
% our differential operator is the 4th derivative, multiplied with the material constans
A = D4*matpar;

% Boundary conditions enforced with the equation C*w=f

% First, we need vectors for indexing
% the points r, at which the boundary conditions are enforced
% These points are removed from the calculation.
% At the remaining points with index k, we will do the actual eigenvalue calculation
%
% since we have 2 boundary conditions at each of the boundaries, we need to
% remove 4 points from the calculation
% the best way to do this is to remove 2-2 points at the 2 ends of the domain
rind = [1,2,N,N+1]; % BC enforced at the two endpoints
kind = 3:N-1; % internal points


% Second, we need to assembly the matrix C
% In the matrix C, the boundary conditions are prescribed in the form: C*w = f
% note the y goes from 'b' to 'a' (in descending order), so we should
% prescribe the boundary conditions in this order
% 
% Dirichlet boundary conditions can be prescribed at a given row of the matrix C 
%       by zeroing out that row, and setting the corresponding coloumn of that row to 1
% Boundary conditions for a derivative of the function at a given point can be prescribed
%        by making equal a row of C with the corresponding row of the differentiation matrix
% Inhomogeneous boundary conditions can be handled by prescribing nonzero values in the vector f
C = [ ...
    D2(1,:); ...           % u''(b) = 0; enforced at the 1st point
    D3(1,:); ...           % u'''(b) = 0; enforced at the 2nd point
    D1(end,:);  ...           % u'(a) = 0;  enforced at the (N)th point
    zeros(1,N),1  ...           % u(a) = 0; enforced at the (N+1)th point
    ];

% since the boundary conditions are homogeneous, the function is zero at the ends
%       this is not used to simplify the code
f = zeros(4,1);
% reordering C, simply by using the index vectors
C_r = C(:,rind);
C_k = C(:,kind);
% calculating the matrix H and G
H = inv(C_r);
G = -H*C_k;

% creating the matrix, which contains the boundary conditions
A_hat = A(kind,kind) + A(kind,rind)*G;

% calculating eigenvalues and eigenfunctions
[V,D] = eig(A_hat);

% sorting eigenvalues, outputs are the sorted eigenvalues and the indices
[eigvals,ind] = sort(diag(D),'descend');

% usind index vector for sorting the eigenfunctions
V = V(:,ind);

% the boundary points can be calculated by ur = Guk + Hf
% note that f is zero
Vr = zeros(4,N-3);
for i = 1:N-3
    Vr(:,i) = G*V(:,i) + H*f;
end
V = [Vr(1:2,:);V;Vr(end-1:end,:)];

%% plotting the eigenvalues
figure; 
subplot(1,2,1);
semilogx(real(eigvals),imag(eigvals),'Marker','o','Linestyle','none');
hold on;
xlabel('$\mathrm{log}(\mathcal{R}(\lambda))$','interpreter','latex');
ylabel('$\mathcal{I}(\lambda)$','interpreter','latex');

title('eigenvalue spectrum');

axis([0    1e10   -1    1]);
set(gcf,'Position',[680   647   931   331]);

% with left click, you can continue to plot the eigenvalues, right clicking stops this
buttonvar = 1;
while buttonvar==1
    
    subplot(1,2,1);
    % getting an eigenvalue
    [xx,yy,buttonvar] = ginput(1);
    % calculating complex eigenvalue from values
    tempeigv = xx + 1i*yy;
    % finding eigenvalue closest to the click
    [~, tempind] = min(abs(tempeigv - eigvals));
    % delete previous plot
    if exist('plth') delete(plth); end
    % plotting selected eigenvalue
    plth = plot(real(eigvals(tempind)),imag(eigvals(tempind)),'x','LineWidth',1.5,'MarkerSize',11,'Color','k');
    % plotting selected eigenvector
    subplot(1,2,2);
    plot(x,V(:,tempind),'LineWidth',1.3, 'Marker','x');
    temp_omega = sqrt(real(eigvals(tempind)));
    title(['\omega = ', num2str(temp_omega),', f = ', num2str(temp_omega/2/pi)]);
end



%% Chebyshev differentiation function
function [D,x] = cheb(N)
if N==0, D=0; x=1; return, end
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';
D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries
D = D - diag(sum(D')); % diagonal entries
end