clear all;clc 
% LAB 5:- Euler Bernoulli and Timoshenko Beam Theory.
% This file contains code to execute problem for point laod, uniformly 
% distributed load and distributed and a point load together for EULER-BERNOULLI beam.

% For the required loadings remove the % sign and comment out the required lines.


global X1 elcon;

%*********** Geometry ***********************************
for step = 3:1:3
nnodes = 2^step+1;
Mn=nnodes-1;
ML = [1];
MR = [nnodes];
YL = [1];
YR = [nnodes];
format long
%*********** Physical Parameters ***********************

left = 0;
right = 48;
X1 = linspace(left,right,nnodes);
E = 29.7E06;
% dimension of the rectangular beam bxhxl = 3x6x48
I = 54; % I=b^3*h/3; 

% for point load
%P1_loc=(nnodes+1)/2;
%p = zeros(nnodes,1);
%p(P1_loc) = 1000;

% for uniformly distributed load
 p = 1000;

% for distributed and a point load
 %P1_loc=(nnodes+1)/2;
 %P2_loc = (3*nnodes+1)/4;
 %p1 = zeros(nnodes,1);
 %p2 = zeros(nnodes,1);
 %p1(1:P1_loc) = 1000;
 %p2(P2_loc) = 1000;

Mleft = 0;
Mright = 0;
Yleft = 0;
Yright = 0;

%*********** Iteration Parameters ***********************
[Xnodes Ynodes] = size(X1);
M = ones(Xnodes.*Ynodes,1)*0;
y = ones(Xnodes.*Ynodes,1)*0;

%********************* Initial Conditions for Set 0 ********************
for ii=1:1:length(ML)
M(ML(ii)) = +Mleft;
end
for ii=1:1:length(MR)
M(MR(ii)) = +Mright;
end

%********************* Equations for Set 0*****************
% for point load
%M0 = +asres1D(1,[],[],-1,'A211L',M)-p;

% for uniformly distributed load
 M0 = +asres1D(1,[],[],-1,'A211L',M)...
     -asres1D(p,[],[],1,'A200L',ones(nnodes,1));

% for distributed and a point load
 % M0 = +asres1D(1,[],[],-1,'A211L',M)-asres1D(1,[],[],1,'A200L',p1)-p2;

%Apply Dirichlet boundary condition to boundary ML.
for ii=1:1:length(ML)
M0(ML(ii)) = 0;
end

%Apply Dirichlet boundary condition to boundary MR.
for ii=1:1:length(MR)
M0(MR(ii)) = 0;
end

%***************** Jacobians for Set 0**********
%Equation M by M
M_jac0 = +asjac1D(1,[],[],-1,'A211L',[]);
DIFFM = M_jac0;
%Apply Dirichlet boundary condition to boundary ML.
for ii=1:1:length(ML)
M_jac0(ML(ii),:) = 0;
M_jac0(:,ML(ii)) = 0;
M_jac0(ML(ii),ML(ii)) = 1;
end

%Apply Dirichlet boundary condition to boundary MR.
for ii=1:1:length(MR)
M_jac0(MR(ii),:) = 0;
M_jac0(:,MR(ii)) = 0;
M_jac0(MR(ii),MR(ii)) = 1;
end

%***************** Solver ********************************************
dM = M_jac0\(-M0);

%************* Solve Iterate *****************************************
M = M + dM;
figure (1)
hold on
plot(X1,M)
xlabel('X1 (in)','Fontsize',16),ylabel('M (lbf-in)','Fontsize',16)
title('Bending Moment Variation','Fontsize',16)

%********************* Initial Conditions for Set 1 ********************
for ii=1:1:length(YL)
y(YL(ii)) = +Yleft;
end
for ii=1:1:length(YR)
y(YR(ii)) = +Yright;
end

%********************* Equations for Set 1*****************
y1 = +asres1D(E.*I,[],[],-1,'A211L',y)...
+asres1D(1,[],[],1,'A200L',M);

%Apply Dirichlet boundary condition to boundary YL.
for ii=1:1:length(YL)
y1(YL(ii)) = 0;
end

%Apply Dirichlet boundary condition to boundary YR.
for ii=1:1:length(YR)
y1(YR(ii)) = 0;
end

%***************** Jacobians for Set 1**********
%Equation y by y
y_jac1 = +asjac1D(E.*I,[],[],-1,'A211L',[]);
DIFFY = y_jac1;
%Apply Dirichlet boundary condition to boundary YL.
for ii=1:1:length(YL)
y_jac1(YL(ii),:) = 0;
y_jac1(:,YL(ii)) = 0;
y_jac1(YL(ii),YL(ii)) = 1;
end

%Apply Dirichlet boundary condition to boundary YR.
for ii=1:1:length(YR)
y_jac1(YR(ii),:) = 0;
y_jac1(:,YR(ii)) = 0;
y_jac1(YR(ii),YR(ii)) = 1;
end

%***************** Solver ********************************************
dy = y_jac1\(-y1);

%************* Solve Iterate *****************************************
y = y + dy;
figure (2)
hold on
plot(X1,y)
xlabel('X1 (in)','Fontsize',16),ylabel('y (in)','Fontsize',16)
title('Displacement Variation','Fontsize',16)

enormM = 0.5*M'*DIFFM*M;

end

