% <FEM for BEAM STRUCTURE>
% governing equation: d2/dx2(EI d2u/dx2) = f
% solution: u(x)
% boundary condition: u(1) = ƒÆ(1) = u(n+1) = 0

% ------------------------------------------------------ %

% element number: n
n = 100;

% beam length: L
L = 1;

% element parameters
dL= L/n; %-- length
EI= 100; %-- flexual rigidity

% node parameters
x = 0:dL:L; % position

% degree of freedom: 2(n+1)
dof     = 2*(n+1);
bc_dof  = [1 2 2*(n+1)-1];
bc_u    = [0; 0; 0;];

% load
p_dof   = 2*(n+2)/2-1; %-- dof of load
f = zeros(dof,1); %-- initialization
f(p_dof)= -1; %-- load

% initialization
u = zeros(dof,1);

% ------------------------------------------------------ %

% fem: stifness matrix
K = zeros(2*(n+1),2*(n+1));
for ii=1:n
    
    K(2*ii-1,2*ii-1) = K(2*ii-1,2*ii-1) + 12*EI/dL^3;
    K(2*ii-1,2*ii  ) = K(2*ii-1,2*ii  ) +  6*EI/dL^2;
    K(2*ii-1,2*ii+1) = K(2*ii-1,2*ii+1) - 12*EI/dL^3;
    K(2*ii-1,2*ii+2) = K(2*ii-1,2*ii+2) +  6*EI/dL^2;
    
    K(2*ii  ,2*ii-1) = K(2*ii  ,2*ii-1) +  6*EI/dL^2;
    K(2*ii  ,2*ii  ) = K(2*ii  ,2*ii  ) +  4*EI/dL  ;
    K(2*ii  ,2*ii+1) = K(2*ii  ,2*ii+1) -  6*EI/dL^2;
    K(2*ii  ,2*ii+2) = K(2*ii  ,2*ii+2) +  2*EI/dL  ;
    
    K(2*ii+1,2*ii-1) = K(2*ii+1,2*ii-1) - 12*EI/dL^3;
    K(2*ii+1,2*ii  ) = K(2*ii+1,2*ii  ) -  6*EI/dL^2;
    K(2*ii+1,2*ii+1) = K(2*ii+1,2*ii+1) + 12*EI/dL^3;
    K(2*ii+1,2*ii+2) = K(2*ii+1,2*ii+2) -  6*EI/dL^2;
    
    K(2*ii+2,2*ii-1) = K(2*ii+2,2*ii-1) +  6*EI/dL^2;
    K(2*ii+2,2*ii  ) = K(2*ii+2,2*ii  ) +  2*EI/dL  ;
    K(2*ii+2,2*ii+1) = K(2*ii+2,2*ii+1) -  6*EI/dL^2;
    K(2*ii+2,2*ii+2) = K(2*ii+2,2*ii+2) +  4*EI/dL  ;
    
end

% boundary condition
f_bc = K(1:dof,bc_dof)*bc_u;
idx = 1:dof;
idx(bc_dof)=[];

% right-hand-side vector
rhs = f-f_bc;

u(idx) = K(idx,idx)\rhs(idx);
u(bc_dof) = bc_u;

% deflection
delta = u(1:2:(dof-1));
theta = u(2:2:dof);

% ------------------------------------------------------ %

% filename
[y,m,d]   = ymd(datetime);
ymdstr    = [num2str(y) num2str(m,'%0.2d') num2str(d,'%0.2d')];
filename  = ['ex_03_' ymdstr '.emf'];

% visualization
figure(1)
set(gcf, 'position', [100   100   480   480 ])
set(gcf, 'color',[1 1 1])
plot(x,delta,'k-')
set(gca, 'position', [0.10  0.08  0.87  0.87])
set(gca, 'fontname', 'Arial', 'fontsize', 12)
xlim([0 L])
print(filename, gcf, '-dmeta', '-r1200')












