function [] = Examples_8_1 () % Linear ODEs by FEA
% For -[d/dx (K du/dx)] + P u - Q = 0 with two BCs
% J.E. Akin, Rice University, Mech 417 2/11/2020
% x = 0 *--k1--*---k2---*--k3--* x = L sequential mesh

% Element nodes 'connection list'(line msh_typ_nodes.txt)
nodes = [1 2 3; 3 4 5; 5 6 7]      ; % three L3 elements

% get the number of elements, nodes, and equations
n_e = size (nodes, 1)  ; % number of elements in mesh
n_m = max (max(nodes)) ; % number of nodes in mesh
n_n = size (nodes, 2)  ; % number of nodes per element
n_g = 1                ; % number of unknowns per node
n_d = n_g * n_m        ; % nunmber of equations in the system
fprintf ('System has %i elements          \n', n_e)
fprintf ('System has %i nodes             \n', n_m)
fprintf ('System has %i nodes per element \n', n_n)
fprintf ('System has %i unknowns per node \n', n_g)
fprintf ('System has %i equations         \n', n_d)

% Define fixed and free displacement degrees of freedom (dof)
% via equation numbers and given BC values (like msh_ebc.txt)
% Fixed = [1  n_m]                          % (for n_g=1)
% Free = [2:1:(n_m-1)]                      % (for n_g=1)
% BC_is = [0. 0.] ; % any boundary condition values at Fixed 
Fixed = [1]                          % (for n_g=1)
Free = [2:1:n_m]                     % (for n_g=1)
BC_is = [0.] ; % any boundary condition values at Fixed 
% (Edit above two lines to change BC locations or values)
fprintf ('Given solution at nodes \n'); disp(Fixed)
fprintf ('are \n')                    ; disp(BC_is)

% Always allocate number of unknowns, U, reactions, F, and
%        the square matrix size S in S U = F assembly
U = zeros (n_d, 1); F = zeros (n_d, 1); S = zeros (n_d, n_d);
% Set the known displacement values (copy via vector subscript)
U(Fixed) = BC_is (:)   ; % (m) copied into BC locations

x = [0:1:(n_m-1)]/(n_m-1) ; % Set the mesh x coordinates
% Define element properties (like msh_properties.txt)
k_j = [1  1  1]; P_j = [1  1  1]; Q_j = [1  1  1]; % defaults

% Loop over elements to assemble (singular) S matrix
for n = 1:n_e    ; % loop over elements ====>> ====>> ====>> ====>>
  e_nodes = nodes (n, 1:n_n)               ; % element connectivity

% Gather nodal coordinates and element properties
  xy_e (1:n_n, 1) = x(e_nodes(1:n_n))       ; % x coord at el nodes
  L_e = abs (xy_e (n_n, 1) - xy_e (1, 1))        ; % element length
  K = k_j(n); P = P_j(n); Q = Q_j(n)          ; % gather properties

% Closed form matrices for constant element properties
  if (n_n == 2)                            ; % linear element L2_C0
    S_e =  K*[1  -1; -1  1]/L_e                ; % stiffness L2_C0
    M_g = L_e*[2   1;  1  2]/6           ; % generalized mass L2_C0
    c_e = L_e*Q *[1; 1]/2        ; % source L2_C0
  elseif (n_n == 3)                     ; % quadratic element L3_C0
    S_e =  K*[7 -8  1; -8 16 -8; ...
               1 -8  7]/(3*L_e)                 ; % stiffness L3_C0
    M_g = L_e*[4  2 -1; 2 16 2; ...
              -1  2  4]/30               ; % generalized mass L3_C0
    c_e = L_e*Q *[1; 4; 1]/6        ; % source L3_C0
  end % if number of nodes on this element
 S_e = S_e + P * M_g               ; % add surroundings sq matrix

% Insert completed element matrices into system matrices
  [rows] =  nodes(n,:)                   ; % for n_g=1
  % rows = vector subscript converting element to system eq numbers
  F (rows)       = F (rows)       + c_e   ; % add to system sources
  S (rows, rows) = S (rows, rows) + S_e ; % add to system stiffness

end ; % for n-th element <==== <==== <==== <==== <==== <====

% These matrix equilibrium equations are not unique
% because they do not yet include Dirichlet BCs.)
% fprintf('Assembled (singular) stiffness matrix \n'); disp(S)

% Transfer known BC displacement effects to Free loads
F (Free) = F (Free) - S (Free, Fixed)*U (Fixed);

% List the non-singular partition of the matrix system
%fprintf ('Free stiffness partition \n'); disp(S(Free,Free))
%fprintf ('Free source \n');              disp(F(Free))

% Invert non-singular S portion, multiply it times the
% Free forces to find the Free displacements
U (Free) = S (Free, Free) \ F (Free) ; % (m)
fprintf ('System unknowns: \n') ; disp(U)

% Optionally find system reactions at the Fixed displacements
% needed to support the Dirichlet BCs
F(Fixed) = S(Fixed, Fixed)*U(Fixed) + S(Fixed, Free)*U(Free) ;
fprintf ('System reactions at Fixed BC: \n'); disp(F(Fixed))
fprintf ('System reactions sum   = %10.4e \n', sum(F(Fixed)));

% graph the solution, then its gradient
% addpath /clear/www/htdocs/mech517/Akin_FEA_Lib
graph_L3_mesh_result (n_e, n_m, n_n, x, nodes, U)    % u(r)
graph_L3_mesh_result (n_e, n_m, n_n, x, nodes, U, 1) % grad

% end Examples_8_1 () % Linear ODEs by FEA
% Running with two Dirichlet BC gives
% Examples_8_1
% System has 3 elements
% System has 7 nodes
% System has 3 nodes per element
% System has 1 unknowns per node
% System has 7 equations
% Fixed = 1     7
% Free  = 2     3     4     5     6
% Given solution at nodes 1     7
% are                     0     0
% System unknowns: [0 0.0635 0.1008 0.1132 0.1008 0.0635 0]
% System reactions at Fixed BC: [-0.4066 -0.4066]
% 
% With current u(0)=0 and du(1)/dx=0 running gives
% Examples_8_1
% System has 3 elements
% System has 7 nodes
% System has 3 nodes per element
% System has 1 unknowns per node
% System has 7 equations
% Fixed = 1
% Free  = 2     3     4     5     6     7
% Given solution at nodes 1 
% are 0
% System unknowns: [0 0.1136 0.2025 0.2692 0.3156 0.3429 0.3519]
% System reactions at Fixed BC: -0.7060

function []=graph_L3_mesh_result (n_e, n_m, n_n, x, ...
                                  nodes, Ans, alt) %=====
% Copyright 2017 J.E. Akin. All rights reserved.
%  ------------------------------------------------------
%  Matlab graph of result values at mesh nodes
%  for a mesh of 3 node quadratic line elements
%   If alt > 0, graph the gradient in each element
%  ------------------------------------------------------
if ( nargin == 6 ); alt = 0 ; end ;
if ( alt == 0 ) % get FEA answers
  y = Ans (:) ; % actual result
  ymin = min (Ans) ; ymax = max (Ans) ;
  Amin = min (Ans) ; Amax = max (Ans) ;
else % alt > 0, get element gradient  limits
  DLH_nodes = [-3  4   -1; -1  0  1; 1  -4  3] ;
%    Loop over all elements
  ymin = realmax ; ymax = -realmax ; % initial gradient range
  for it = 1:n_e ;
%      Extract element connectivity 
    t_nodes = nodes (it, 1:3) ;
%      Skip point elements, if any
    if ( all (t_nodes) ) % then valid line
%      Extract element coordinates & values
      t_x  = x (t_nodes) ;         % x at those nodes, only
      t_y  = Ans (t_nodes)         ; % Ans at those nodes, only
      Len = max (t_x) - min (t_x) ; % physical length
      DGH_nodes = DLH_nodes / Len ; % dH/dx at nodes
      y_set = DGH_nodes * t_y ; % node gradients
      smax = max (y_set) ; smin = min (y_set);
      if ( smax > ymax ) ; ymax = smax ; end % if
      if ( smin < ymin ) ; ymin = smin ; end % if
    end % if zero in connectivity
  end % for over all elements
end % if get node values or element gradient

% determine plot ranges
xmax = max (x)   ; xmin = min (x) ;
ydiff = abs (ymax - ymin) / 20 ;
ymax = ymax + ydiff ;
ymin = ymin - ydiff ;

%  Initialize plot labels and title
clf                               % clear graphics
axis ([xmin, xmax, ymin, ymax])               % set axes
hold on                           % hold image for plots
grid on
xlabel ('X') ;
if ( alt == 0 )
  ylabel (['Solution: min= ', num2str(Amin, '%8.3e'), ' max= ', ...
             num2str(Amax, '%8.3e')])
  title  (['Quadratic FEA Result:  ', ...
          int2str(n_e), ' Elements, ', int2str(n_m),' Nodes, ', ...
         ' (', int2str(n_n), ' per element)'])
else % alt > 0, get root mean sq
  ylabel (['Gradient: min= ', num2str(ymin, '%8.3e'), ' max= ', ...
             num2str(ymax, '%8.3e')])
  title  (['Quadratic FEA Result Gradient: ', ...
          int2str(n_e), ' Elements, ', int2str(n_m),' Nodes, ', ...
         ' (', int2str(n_n), ' per element)'])
end % if

%    Loop over all elements
  for it = 1:n_e ;
%      Extract element connectivity 
    t_nodes = nodes (it, 1:3) ;
%      Skip point elements, if any
    if ( all (t_nodes) ) % then valid line
%      Extract element coordinates & values
      t_x  = x (t_nodes) ;         % x at those nodes, only
      t_y  = Ans (t_nodes) ;         % Ans at those nodes, only
%      Plot the element number
      x_bar = sum (t_x' )/n_n ;
      y_bar = sum (t_y' )/n_n ;

%      Loop over local points on the quadratic polynomial element
      Len = max (t_x) - min (t_x) ;
      n_poly = ceil (75 / n_e) ; last = n_poly + 1 ;
      for k = 1:last        ; % points in parametric space
%          get element parametric interpolation functions
         r = (k - 1)/n_poly ;   % on 0 to 1
         H = [(1 - 3*r + 2*r^2) (4*r - 4*r^2) (2*r^2 - r)] ;
         x_el (k) = H * t_x' ;  % true x value
         if ( alt == 0 )
           y_el (k) = H * t_y ;  % true y value
         else
           DLH = [(-3 + 4 * r)  (4 - 8*r)  (4*r -1)] ; % dH/dr
           DGH = DLH / Len ; % dH/dx
           y_el (k) = DGH * t_y ;  %
         end % if
      end % for k
      y_bar = mean (y_el) ;
      t_text = sprintf ('   (%i)', it);  % offset # from pt
      text (x_bar, y_bar, t_text, 'Rotation', 30) % incline
      plot (x_el, y_el, 'b-', 'LineWidth', 2)
      plot (t_x(1), y_el(1), 'k*')          ; % first node
      plot (t_x(n_n), y_el(last), 'k*')     ; % last node
      plot (mean(t_x), mean(y_el), 'k*')    ; % interior node
    end % if zero in connectivity
  end % for over all elements

if ( alt == 0 )
  print -dpng result_L3_mesh ; % save the graph
  fprintf ('Created file result_L3_mesh.png \n')
else
  print -dpng result_L3_grad ; % save the graph
  fprintf ('Created file result_L3_grad.png \n')
end ; % if which plot
% end graph_L3_mesh_result ==================================