clear all; clc
% FEm.PSE template script, simply supported Euler-Bernoulli beam,
% cubic Hermite basis.

% set up physical parameters
E = 30.0E6;          %  elastic modulus
Izz = 20.8333;       %  moment of inertia
p = 8.3333;          %  loading, force per length
L = 120.0;           %  length of the beam

% set up mesh
nnodes = 33;       % 32 elements, two DOF per node
X1 = linspace(0,L,nnodes);

% set up nodal load data in Hermite form.
% need value and slope of load distribution at each node
loadvals = ones(1,nnodes) * p; % load values at each node
slopevals = zeros(1,nnodes);   % slope of the load distribution at each node

% loadvals and slopevals were formed as rows.  Concatenate them one on
% top of the other, to form a 2 by nnodes array.
PP = [loadvals;slopevals];
% now stack the columns of PP on top of each other to form one long
% column.  The first entry in P is the value of the load at node 1, the
% second entry is the slope of the load at node 1, the third entry is the 
% value of the load distribution at node 2, ect
P = PP(:)

load femHlib;                  % load cubic Hermite element library matrices

K = ass1dH(E*Izz,[],[],-3,A222H,X1);    % form the global stiffness matrix
b = ass1dH([],[],[],1,A200H,X1) * P;    % form the global load column matrix

%  impose zero displacement boundary condtions at each end
Kdiri = K;                               % copy K
Kdiri(1,:) = zeros(1,2*nnodes);          % left end displacement
Kdiri(1,1) = 1;                          % left end displacement
b(1,1) = 0;                              % left end displacement
Kdiri(2*nnodes-1,:) = zeros(1,2*nnodes); % right end displacement
Kdiri(2*nnodes-1,2*nnodes-1) = 1;        % right end displacement
b(2*nnodes-1,1) = 0;                     % right end displacement

% no action required for zero moment at each end

% solve the linear system 
Y = Kdiri \ b

format long;
enorm = (1/2) * (Y' * K * Y)  %  energy norm
format short;
