% hw5_7.m
%
% 1D linear elastic finite element code for 
% Theory of Finite Element Analysis
% Jack Chessa

% This code will serve as teh basis for several hws. Eventually, 
% you will have a general purpose 1D thermo elastic finite 
% element program.
%
%
%-------------------------------------------------------
% HW#5
%  Modify the program so that it can apply inhomogenious
%  displacement boundary conditions (ie d_1 not = 0)
% 
%-------------------------------------------------------
% HW#6
%  Modify the program so that it can compute the stress
%  in each element and plot (make sure you have the 
%  stress jumps)
%
%-------------------------------------------------------
% HW#7
%  Modify the program so that it can compute the nodal
%  averaged stress in each element and plot 
%
% ******************************************************
% *              I N P U T    S E C T I O N            *
% ******************************************************
% we solve problem 2 of the hand out
A=1;    % problem parameters
L=2;
E=100;
P=-1;
del=0.1;

node=[0.0 L/2 L]';    % node coordinates
element=[1 2;2 3];    % element connectivcity

% to handle differing AE's for each element we will use
% an array of length equal to the number of elements
% so AE(e) is the AE for element e.
area=[2 1]*A;
young=[1 1]*E;

% here we add data structures to describe general essential
% BCs inhomogenious and homogenious
ifix=[1 3];    % nodes where essential bc is applied
ival=[0 del];  % displacement at those nodes

% and general nodal point loads
ffix=[2];   % nodes where load is applied
fval=[P];

% ******************************************************
% *       P R O C E S S I N G   S E C T I O N          *
% ******************************************************
nn=size(node,1);     % number of nodes in mesh
ne=size(element,1);  % number of elements in the mesh

K=zeros(nn,nn);      % global stiffness matrix
d=zeros(nn,1);       % global nodal displacement vector
f=zeros(nn,1);       % global force vector

% ----- COMPUTE AND ASSEMBLE FEM STIFFNESS MATRIX ------ 
% 1) loop over the elements 
for e=1:ne
      
%   a) compute element stiffness matrix ke
    n1=element(e,1); % node number for node 1 of the element e
    n2=element(e,2); % node number for node 2 of the element e
    le=abs(node(n2)-node(n1)); % length of element e
    ae=area(e)*young(e);
    
    ke=ae/le*[1 -1;-1 1]; % element stiffness matrix
    
%   b) scatter element stiffness matrix ke into K
    sctr=element(e,:); % see online notes for details
    K(sctr,sctr)=K(sctr,sctr)+ke;

% end of element loop
end

% ---------------- COMPUTE FORCE VECTOR ---------------- 
f(ffix)=fval;  % apply nodal point loads

% -------- APPLY ESSENTIAL BOUNDARY CONDITIONS --------- 
% first swing reaction forces from prescribed disps to rhs
%      f = f - K_colI * dI  for all I where disp is prescrbed
f=f-K(:,ifix)*ival';
% next modify K matrix to equlibrium equations become disp 
% constraint equations.
K(ifix,:)=0; % zero out first row
K(:,ifix)=0; % and column
K(ifix,ifix)=eye(length(ifix)); % put a one at diagonal
% finally set vals in force vector so we have the correct 
% rhs of the disp constraints
f(ifix)=ival;

% -------- SOLVE DISCRETE FE EQUATIONS (Kd=f) ---------- 
d=K\f;

% ******************************************************
% *   P O S T - P R O C E S S I N G   S E C T I O N    *
% ******************************************************
clf
subplot(2,1,1)
plot(node,d,'bo-')
xlabel('x')
ylabel('displacement')

% ok let's plot the element stresses
subplot(2,1,2)
hold on
stress=zeros(1,nn); % nodal average stress
nsamp=zeros(1,nn);  % number of elements connected to node
for e=1:ne % loop over elements
  
  s11=young(e)*(d(e+1)-d(e))/(node(e+1)-node(e))
  plot(node(e:e+1),s11*[1 1],'b-');  % true element stress
  
  % add to the average nodal stresses
  stress(e:e+1) = stress(e:e+1) + s11;
  nsamp(e:e+1)  = nsamp(e:e+1)  + 1;
end
stress=stress./nsamp;
plot(node,stress,'r:') % plot nodal average stress

xlabel('x')
ylabel('stress')