function [] = Two_span_three_L3_C1_beam () % revised 2/27/20
% Requires: addpath /clear/www/htdocs/mech517/Akin_FEA_Lib
% Uses manual assembly loop

%     w=2,400 lb/ft                  P=10,000 lb
%          (1)            (2)        |  (3)
%     wwwwwwwwwwwww                  V
%   1 *-----*-----*--------*---------*---*---*
% roller    2     3,roller 4         5   6   7,fixed
%   x=0     5    10       16        22  25  28 ft.
% DOF 1,2  3,4  5,6      7,8      9,10 11,12 13,14

% EI == 1 lb ft^2.  Reaction results are:
% R1 = 9,957 R3 = 17,226 R7 = -6,816 lb. M7 = -23,121 ft lb

%                  problem controls
n_g = 2 ;                             % number of DOF per node
n_e = 3 ;                                 % number of elements
n_m = 7 ;                        % number of nodes in the mesh
n_n = 3 ;                        % number of nodes per element
n_d = n_g * n_m ;                 % number of system equations
n_i = n_g * n_n ;                       % nmber of element DOF

%                  allocate input arrays
x = zeros (n_m, 1) ; Connect = zeros (n_e, n_n) ;     % inputs
Connect = [1  2  3 ;                   % element 1 connections
           3  4  5 ;                   % element 2 connections
           5  6  7 ] ;                 % element 3 connections
x   = [0 5 10 16 22 25 28]'           ; % system x-coordinates
EIs = [1 1 1]                            ; % element EI values

%                  allocate required arrays
S   = zeros (n_d, n_d); c   = zeros (n_d, 1);  % system arrays
S_e = zeros (n_i, n_i); c_f = zeros (n_i, 1);    % elem arrays
index = zeros (n_i, 1);              % element scatter numbers
u     = zeros (n_d, 1);        % system deflections and slopes
React = zeros (n_d, 1);                     % system reactions

%           zero EBC and free DOF vector subscripts
Fixed = [1 5 13 14] ; Free = [2 3 4 6 7 8 9 10 11 12];
%                  insert point load scatter
P_n = 5 ;  P_v = -1e4  ; % node & value of vertical point load
P_eq = n_g*(P_n - 1) + 1 ; % equation of vertical displacement
c(P_eq) = P_v      ; % insert (scatter) point load into system

e = 1 ; % begin element loop e=1 ====> ====>
e_nodes = Connect (e, 1:3)           ;  % nodes on the element
L = abs ( x(e_nodes(3)) - x(e_nodes(1)) )   ; % compute length
EI = EIs (e) ;                           % constant properties
Line_Load = -2.4e3 * [1  1  1]' ;         % constant line load
[index] = get_element_index (n_g, n_n, e_nodes); % sys numbers
[S_e] = matrix_S_E_L3_C1 (EI, L)         ; % bending stiffness
[c_f] = matrix_c_f_L3_C1 (L, Line_Load)     ; % load resultant
S (index, index) = S (index, index) + S_e  ; % e = 1 assembled
c (index)        = c (index)        + c_f  ; % e = 1 assembled

e = 2 ; % continue elem loop e=2 ====> ====> ====>
EI = EIs (e) ;                           % constant properties
e_nodes = Connect (e, 1:3)            ; % nodes on the element
L = abs ( x(e_nodes(3)) - x(e_nodes(1)) )   ; % compute length
[index] = get_element_index (n_g, n_n, e_nodes); % sys numbers
[S_e] = matrix_S_E_L3_C1 (EI, L)       ; % bending stiffness
S (index, index) = S (index, index) + S_e  ; % e = 2 assembled

e = 3 ; % finish elem loop e=3 ====> ====> ====> ====> ====>
EI = EIs (e) ;                           % constant properties
e_nodes = Connect (e, 1:3)            ; % nodes on the element
L = abs ( x(e_nodes(3)) - x(e_nodes(1)) )   ; % compute length
[index] = get_element_index (n_g, n_n, e_nodes); % sys numbers
[S_e] = matrix_S_E_L3_C1 (EI, L)       ; % bending stiffness
S (index, index) = S (index, index) + S_e; % assembly complete
% end loop over elements   <==== <==== <==== <==== <==== <====

% independent equations K * U = C.         Partition and solve
K = S (Free, Free) ; C = c (Free) ;             % known arrays
U = K \ C                           ; % non-zero displacements
u (Free) = U                    ; % gather back into full list
fprintf('\n Beam displacements u = ') ; disp(u')     ; % print

React(Fixed) = S(Fixed, Free)*U - c(Fixed); % system reactions
fprintf('\n Beam reactions = ') ; disp(React'); % pretty print
graph_L3_C1_result (u, x, Connect)       ; % plot the solution
graph_L3_C1_moment (u, x, Connect, EIs)    ; % plot the moment 
graph_L3_C1_shear  (u, x, Connect, EIs)     ; % plot the shear
% end Two_span_three_L3_C1_beam ==============================

function [S_e] = matrix_S_E_L3_C1 (EI, L) %===================
% E = elastic modulus, I = moment of inertial, L = elem length
L2 = L^2 ; L3 = L^3 ;                              % shorthand
%  Constant property quintic element bending stiffness
S_e =[5092,   1138*L, -3584,    1920*L,  -1508,   242*L  ; 
      1138*L, 332*L2, -896*L,   320*L2,  -242*L,  38*L2  ; 
     -3584,   -896*L,  7168,    0,       -3584,   896*L  ; 
      1920*L, 320*L2,  0,       1280*L2, -1920*L, 320*L2 ; 
     -1508,  -242*L,  -3584,   -1920*L,   5092,  -1138*L ; 
      242*L,  38*L2,   896*L,   320*L2,  -1138*L, 332*L2 ] ... 
      *EI/(35 * L^3) ;                             % stiffness 
% end  matrix_S_E_L3_C1  %====================================

function [c_f] = matrix_c_f_L3_C1 (L, Line_Load)
% Quintic beam resultant shears and moments from
% a quadratic line load through [f1 f2 f3]
% L   = element length 
% c_f = resultant at nodes, R_e * Line_Load
% R_e = rectangular comversion matrix
% Line_Load = [f1 f2 f3]' the load per length at
% left, middle, and right node
 R_e = L * [57   44   -3   ; % 6 by 3 rectangle
            3*L  4*L   0   ;
            16   192   16  ;
           -8*L  0     8*L ;
           -3    44    57  ;
            0   -4*L  -3*L ] / 420 ;
c_f = R_e * Line_Load ;       % 6 by 1 resultant
% end matrix_c_f_L3_C1 %========================