%% 1D bar
clc
close all
clear
syms x

%% data for steel
E  = 210e+9;    A = 25e-4;   rho = 7800;  L=1;
px=4.0e6*x; P=-1e6;

%% discretization
Nel = 20;  Ndof = Nel+1;  coord=linspace(0,L,Ndof);

%% assembly
KG=zeros(Ndof,Ndof); MG=zeros(Ndof,Ndof);
fG=zeros(Ndof,1);
fG(Ndof)=P;

for iel=1:Nel
    x1=coord(iel); x2=coord(iel+1);
    N1=(x-x2)/(x1-x2); N2=(x-x1)/(x2-x1);
    N=[N1 N2]; B=diff(N);
    Kel=int(B'*A*E*B,x1,x2);
    Fel=int(N'*px,x1,x2);

    edof=[iel, iel, iel+1];
    [KG,fG]=assem(edof,KG,Kel,fG,Fel);

    Mel=int(A*rho*N'*N,x1,x2);
    MG=assem(edof,MG,Mel);
end

%% solution
bc = [1,0];
u=solveq(KG,fG,bc);

figure(1) % displacement
title('displacement function')
plot(coord,u,'linewidth',2)
grid on, xlabel('x'), ylabel('u')

for iel=1:Nel
    x1=coord(iel); x2=coord(iel+1);
    u1=u(iel); u2=u(iel+1);    
    xi=linspace(x1,x2,10); % sample points

    dN1=1/(x1-x2);       % shape function derivative
    dN2=1/(x2-x1);

    Shxi=A*E*(u1*dN1+u2*dN2); 
    figure(2), plot(xi,Shxi,'b.'), hold on
    title('normal force function elementwise')
    grid on, xlabel('x'), ylabel('N^e')
    
    % N1=(xi-x2)/(x1-x2); N2=(xi-x1)/(x2-x1);
    % uhxi=u1*N1+u2*N2; figure(2), plot(xi,uhxi), hold on
end

%% eigen vibrations
[La,Ev]=eigen(KG,MG,bc(:,1)); Freq=sqrt(La)/(2*pi)

figure(3)
title('basic eigen modes')
color=['r','g','b','m'];
modes=min([4,Nel]);
for iev=1:modes
    subplot(2,2,iev)
    plot(coord,Ev(:,iev),color(iev),'linewidth',2), 
    grid on, axis([0,L,-1,1]), hold on
    nm=num2str(iev); forma=strcat('Eigenform ',nm);
    xlabel('x'), ylabel('Ev')
    title(forma)
    plot([0,L],[0,0],'k')
end

%% response history - Newmark integration method

% kinematic b.c. - remove relevant rows and columns
MG(1,:)=[]; MG(:,1)=[];
KG(1,:)=[]; KG(:,1)=[];
fG(1)=[];
Ndof_1=Ndof-1;

% load
freq_imp=Freq(1)*1.5; om=2*pi*freq_imp;
time=0; t_end=0.002;
f=inline('sin(o*t)');

% time discretization
N=200;  % time steps
dt=(t_end-time)/N;

% initial
uu=zeros(Ndof_1,1); vv=zeros(Ndof_1,1);
thist(1)=time;
uhist(1)=uu(Ndof_1);

% Newmark
Amat=MG+1/2*dt^2*KG;
Amat=inv(Amat);

for itm=2:N+1
    time=time+dt; 
  
    aa=Amat*(fG*f(om,time)-KG*(uu+dt*vv));
    uu=uu+dt*vv+1/2*dt^2*aa;
    vv=vv+dt*aa;

    thist(itm)=time;
    uhist(itm)=uu(Ndof_1);
end

figure(4)
plot(thist,uhist,'r','linewidth',2), grid on, hold on
xlabel('t'); ylabel('u(1,t)');
