% FEm.PSE template script, simply supported Euler-Bernoulli beam,
% linear Lagrange basis

global X1

% set up problem parameters
E = 30e6;           % elastic modulus in psi
Izz = 20.83333;     % moment of inertia in in^4
L = 120;            % length of beam in inches	
p = 8.33333;        % load in lbf/in
nnodes = 33;        % number of nodes

% setup mesh
X1 = linspace(0,L,nnodes);

% set up nodal load data
P = p*ones(nnodes,1);

load femlib;

% set up M system matrix
Msys =  asjac1D([],[],[],-1,'A3011L',[]);
Msysdiri = Msys; % copy Msys
Msysdiri(1,:) = zeros(1,nnodes);%  apply diri to node 1
Msysdiri(1,1) = 1;
Msysdiri(nnodes,:) = zeros(1,nnodes);% apply diri to last node
Msysdiri(nnodes,nnodes) = 1;

% set up M right hand side
Mrhs =  asres1D([],[],[],1,'A3000L',P);
Mrhs(1,1) = 0;%  enforce diri condition at node 1
Mrhs(nnodes,1) = 0;%  enforce diri at last node

% solve for M
M = Msysdiri \ Mrhs

% set up Y system matrix
Ysys = asjac1D(E*Izz,[],[],-1,'A3011L',[]);
Ysysdiri = Ysys;
Ysysdiri(1,:) = zeros(1,nnodes);%  apply diri to node 1
Ysysdiri(1,1) = 1;
Ysysdiri(nnodes,:) = zeros(1,nnodes);% apply diri to last node
Ysysdiri(nnodes,nnodes) = 1;

% set up Y right hand side
Yrhs = asres1D([],[],[],1,'A3000L',M);
Yrhs(1,1) = 0;%  enforce diri condition at node 1
Yrhs(nnodes,1) = 0; % enforce diri at last node

% solve for Y
Y = Ysysdiri \ Yrhs

% calculate the energy norms
format long
enormM = 1/2 * M' * Msys * M
enormY = 1/2 * Y' * Ysys * Y
format short
