% FEm.PSE template: Euler-Bernoulli Beam, Moment/Deflection Formulation
% Simply supported beam, uniform distributed load

clear all;
clf reset;
format long e;		% double precision exponent output

global X1;        % array for node coordinates

% constant data
E = 4.18E9;			% Young's Modulus (steel - lb/ft^2)
Ml = 0;				% No moment on left due to simple support
Mr = 0;		      % No moment on right due to simple support
Yl = 0;				% No deflection on left due to simple support
Yr = 0;				% No deflection on right due to simple support

% mesh parameters
XL = 0;
XR = 10;				% Ten foot span
% uniform discretization, M=1
nnodes = 2;
X1 = linspace(XL,XR,nnodes);

% nodal element data
	% distributed load - lb/ft
	P=-100*ones(size(X1));

	% moment of inertia (constant/variable)
	b_const = 0.5;			% base of beam - constant (feet)
	h_const = 0.75;		% height of beam - constant (feet)
   INRT = 1/12*(b_const*h_const^3);
   %you can play with variable moment of inertia if you are interested
   %NOT REQUIRED FOR LAB
   %b_var = 0.5							% base of beam - variable (feet)
   %h_var = 1-0.1*X1+0.01*X1.^2;		% base of beam - variable (feet)
	%INRT = 1/12*(b_var.*h_var.^3);

%-----------------------------------------------------------------------------
% load the FE matrix library
load femlib;

%-----------------------------------------------------------------------------
% SOLVE FIRST SYSTEM FOR MOMENT DISTRIBUTION
% Assemble all matrices on the LHS:
	% Assemble Matrix for GWSe term [DIFF]e
	

   % Build LHS
   
   
% Assemble all vectors on the RHS
	% Assemble Vector for GWSe term {LOAD}e
   
   
   % Build RHS
   
   
% Modify RHS for dirichlet data


% Modify LHS for dirichlet data
									% Copy LHS
									% zero out first row
									% i.e., 1*Q1 = Ml
									% Zero out the last row
									% i.e., 1*QM = Mr

% Solve the linear system for MOMENT
	

%-----------------------------------------------------------------------------
% SOLVE SECOND SYSTEM FOR DEFLECTION DISTRIBUTION
% Assemble all matrices on the LHS:
	% Assemble Matrix for GWSe term [DIFF]e
		

   % Build LHS
   
   
% Assemble all vectors on the RHS
   % Assemble Vector for GWSe term {MOM}e
   	
	
   % Build RHS
   	
   
% Modify RHS for dirichlet data
	
		

% Modify LHS for dirichlet data
									% Copy LHS
									% zero out first row
									% i.e., 1*Q1 = Yl
									% Zero out the last row
									% i.e., 1*QM = Yr

% Solve the linear system for DEFLECTION


%-----------------------------------------------------------------------------
% SOLVE FOR SHEAR - POST PROCESSING

% Assemble all matrices on the LHS:
	% Assemble Matrix for GWSe term [MASS]e
	

   % Build LHS
   
   
% Assemble all vectors on the RHS
   % Assemble Vector for GWSe term {DMOM}e
   

   % Build RHS
   

% Solve the linear system for SHEAR


%-----------------------------------------------------------------------------
% SOLVE FOR SLOPE - POST PROCESSING

% Assemble all matrices on the LHS:
	% Assemble Matrix for GWSe term [MASS]e
	

   % Build LHS
   
   
% Assemble all vectors on the RHS
   % Assemble Vector for GWSe term {DDEF}e
   

   % Build RHS
   

% Solve the linear system for SLOPE

%-----------------------------------------------------------------------------
% Compute Energy Norm using [DIFF] (without Dirichlet BC)

%-----------------------------------------------------------------------------
% Make a nice shear/moment/slope/deflection diagram
% NOTE: you may need to make some changes so your data fits!

figure(1)
subplot(4,1,1)
area(X1,V)
hold on
plot(X1,V,'ro')
ylabel('Shear')
title('Euler-Bernoulli Beam')
hold off

subplot(4,1,2)
area(X1,M)
hold on
plot(X1,M,'ro')
ylabel('Moment')
hold off

subplot(4,1,3)
area(X1,THTA)
hold on
plot(X1,THTA,'ro')
ylabel('Slope')
hold off

subplot(4,1,4)
area(X1,Y)
ylabel('Deflection')
hold on
plot(X1,Y,'ro')
xlabel('X')
hold off