function Axial_Taper_Bar
clc; clear all; clf
% Tapered bar example
number_of_elements=4;
L=10; %Total length, inch
e = 10.4*10^6;% modulus of elasticity
P = 1000; % external load
%nodes = [0:2.5:10];
nodes = [0:L/number_of_elements:L];
A = 0.25-0.0125.*nodes ;
% Connectivity Matrix
lmm = [1,2; 2,3; 3,4; 4,5];
K=zeros(number_of_elements+1);
% Generate stiffness matrix for each element and assemble it.
for i=1:number_of_elements
    lm=lmm(i,:);
    k=AxialDefElement(e, (A(i)+A(i+1))/2, nodes(lm));
    K(lm, lm) = K(lm, lm) + k;
end
K
% Define the load vector
F = zeros(number_of_elements+1,1); F(number_of_elements+1)=P

% Nodal solution and reactions
[d, reactions] = NodalSoln(K, F, [1], zeros(1,1))
results=[];
for i=1:number_of_elements
    results = [results; AxialDefResults(e, A(i), ...
            nodes(lmm(i,:)), d(lmm(i,:)))];
end
format short g
results
plot(nodes,d,'o-','LineWidth',2,'MarkerSize',10),title('Axial displacement'), xlabel('x'),ylabel('u')


%function
function k = AxialDefElement(e, A, coord)
% k = AxialDefElement(e, A, coord)
% Generates stiffness matrix of an axial deformation element
% e = modulus of elasticity
% A = Area of cross-section
% coord = coordinates at the element ends

x1=coord(1); x2=coord(2);
k = e*A/(x2-x1)*[1,-1; -1,1];

function rq = AxialDefLoad(q, coord)
% rq = AxialDefLoad(q, coord)
% Generates equivalent load vector for an axial deformation element
% q = uniformly distributed load
% coord = coordinates at the element ends

x1=coord(1); x2=coord(2);
rq = q*(x2-x1)/2*[1;1];

function [d, rf] = NodalSoln(K, R, debc, ebcVals)
% [nd, rf] = NodalSoln(K, R, debc, ebcVals)
% Computes nodal solution and reactions
% K = global coefficient matrix
% R = global right hand side vector
% debc = list of degrees of freedom with specified values
% ebcVals = specified values
dof = length(R);
df = setdiff(1:dof, debc);
Kf = K(df, df);
Rf = R(df) - K(df, debc)*ebcVals;
dfVals = Kf\Rf;
d = zeros(dof,1);
d(debc) = ebcVals;
d(df) = dfVals;
rf = K(debc,:)*d - R(debc);

function results = AxialDefResults(e, A, coord, dn)
% results = AxialDefResults(e, A, coord, dn)
% e = modulus of elasticity
% A = Area of cross-section
% coord = coordinates at the element ends
% dn = displacements at element ends
% The output variables are axial strain, axial stress,
% and axial force. 

x1=coord(1); x2=coord(2); L=x2-x1;
eps= [-1,1]/L*dn;
sigma = e*eps;
force = sigma*A;
results=[eps, sigma, force];