function Cubic_Beam_on_Winkler (load_pt)
% load_pt = 1 flags that point forces and/or couples will be 
%  input via file msh_load_pt.tmp
%.............................................................
%     Classic cubic beam on elastic foundation, 
%     for point loads & couples, line load
%     Mech 400 Adv. Mech. Materials, Rice University
%.............................................................
if ( nargin == 0 )          ;  % check for optional data flag
  load_pt = 0               ;  % no point source data
end % if from argument count

%         Application and element dependent controls
n_g = 2     ; % number of DOF per node (deflection, slope)
n_q = 0     ; % number of quadrature points required
n_r = 1     ; % number of rows in B_e matrix

%              Read mesh nodal input data file
[n_m, n_s, P, x, y, z]          = get_mesh_nodes     ;

%      Extract EBC flags from packed integer flag P 
[EBC_flag] = get_ebc_flags (n_g, n_m, P)      ; % unpack flags
EBC_count  = sum( sum ( EBC_flag > 0 ) )      ; % # of EBC
if ( EBC_count == 0 )
  fprintf ('WARNING: No displacement boundary condition given. \n')
else
  fprintf ('Note: expecting %g displacement BC values. \n', EBC_count)
end % if 

%              Read mesh element input data file
[n_e, n_n, n_t, el_type, nodes] = get_mesh_elements  ;
n_d = n_g*n_m                ; % system degrees of freedom (DOF)
n_i = n_g*n_n                ; % number of DOF per element
S = zeros (n_d, n_d) ; M = zeros (n_d, n_d) ; % stiff & mass
T = zeros (n_d, 1)                          ; % displacements
 
%                  Read EBC values, if any
if ( EBC_count > 0 )                            ; % need EBC data
  [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag) ; % read data
end % if any EBC data expected

%    Read point loads & couples data, if any, and insert in C
C = zeros (n_d, 1)                       ; % force & moments
if ( load_pt > 0 )                       ; % need point loads data
  [C] = get_and_add_point_sources (n_g, n_m, C); % add point loads
end % if any point source expected
C_react = C ; % save for reaction use, before BC changes

%    Load the element properties array
load msh_properties.tmp ; % one row per element
n_prop = size(msh_properties, 1)   % == 1 if homogeneous
fprintf ('\n(Echoing file msh_properties.tmp) \n')
Line_e (1:2) = 0 ;
if ( n_prop == 1 )
  E            = msh_properties (1, 1)   ; % material modulus
  I            = msh_properties (1, 2)   ; % section inertia
  Line_e (1:2) = msh_properties (1, 3:4) ; % trapezoidal line load
  k_f          = msh_properties (1, 5)   ; % foundation stiffness
  Rho          = msh_properties (1, 6)   ; % beam mass per unit length
%                   ECHO PROPERTIES
  fprintf ('Homogeneous Element Properties: \n' )
  fprintf ('Elastity modulus (N/m^2)     = %g \n', E)     
  fprintf ('Moment of inertia (m^4)      = %g \n', I)
  fprintf ('Line Load (N/m)              = [ %g %g ] \n', ...
            Line_e(1), Line_e(2))
  fprintf ('Foundation stiffness (N/m^2) = %g \n', k_f)
  fprintf ('Mass per unit length (kg/m)  = %g  \n', Rho)
end % if
if ( n_prop > n_e )
  error ('ERROR: number of property sets exceeds number of elements')
end % if
F_k (1:n_prop) = msh_properties (1:n_prop, 5) ; % extract foundation k's
if ( any ( F_k > 0 )) % BEF pressures exist 
  BEF = 1 ;
  fprintf ('\nWARNING: only deflections & foundation pressures are accurate \n')
else
  BEF = 0 ;
end % if reaction flag

%    GENERATE ELEMENT MATRICES AND ASSYMBLE INTO SYSTEM
%  Assemble n_d by n_d square matrix terms from n_e by n_e

for j = 1:n_e          ; % loop over elements ====>> ====>> ====>> ====>>
  S_e = zeros (n_i, n_i)                          ; % clear array
  NtN = zeros (n_i, n_i)                          ; % clear array
  M_e = zeros (n_i, n_i)                          ; % clear array
  C_e = zeros (n_i, 1)                            ; % clear arrays
  e_nodes = nodes (j, 1:n_n)                      ; % connectivity
  
%               SET ELEMENT PROPERTIES & GEOMETRY
  L_e = x(e_nodes(2)) - x(e_nodes(1)) ; % beam length
  Line_e (1:2) = 0 ;
  if ( n_prop > 1 )
    E            = msh_properties (j, 1)   ; % material modulus
    I            = msh_properties (j, 2)   ; % section inertia
    Line_e (1:2) = msh_properties (j, 3:4) ; % trapezoidal line load
    k_f          = msh_properties (j, 5)   ; % foundation stiffness
    Rho          = msh_properties (j, 6)   ; % beam mass per unit length
%                   ECHO PROPERTIES
    fprintf ('\nProperties for element %g \n', j)
    fprintf ('Elastity modulus (N/m^2)     = %g \n', E)     
    fprintf ('Moment of inertia (m^4)      = %g \n', I)
    fprintf ('Line Load (N/m)              = [ %g %g ] \n', ...
              Line_e(1), Line_e(2))
    fprintf ('Foundation stiffness (N/m^2) = %g \n', k_f)
    fprintf ('Mass per unit length (kg/m)  = %g  \n', Rho)
  end % if

%      ELEMENT CONDUCTION AND INTERNAL SOURCE MATRICES 
%  Constant cubic element & foundation stiffnesses & mass matrices
  S_e = (E*I/L_e^3)*[ 12,    6*L_e,   -12,    6*L_e     ;
                      6*L_e, 4*L_e^2, -6*L_e, 2*L_e^2   ;
                     -12,   -6*L_e,    12,   -6*L_e     ;
                      6*L_e, 2*L_e^2, -6*L_e, 4*L_e^2 ] ; % stiffness

  NtN = L_e*[ 156,     22*L_e,   54,     -13*L_e  ;
              22*L_e,  4*L_e^2,  13*L_e, -3*L_e^2 ;
              54,      13*L_e,   156,    -22*L_e  ;
             -13*L_e, -3*L_e^2, -22*L_e,  4*L_e^2 ]/420; 

  S_e = S_e + k_f * NtN ; % add foundation stiffness before assembly
  M_e = Rho * NtN       ; % mass matrix

%  Map line load to node forces & moments; C_e = p_To_F * Line_e 
    C_e = zeros (n_i, 1)           ; % clear arrays
    if ( any (Line_e) )            ; % then form forcing vector
      p_To_F = L_e * [ 21,    9     ; 
                       3*L_e, 2*L_e ; 
                       9,     21    ; 
                      -2*L_e -3*L_e ] / 60 ; % cubic H, linear Line
      % 4 x 2 * 2 x 1 = 4 x 1 result 
      C_e = p_To_F (1:n_i, 1:n_n) * Line_e' ; % force moment @ nodes
    end % if or set up resultant node loads for line load

%         SCATTER TO (ASSEMBLE INTO) SYSTEM ARRAYS
% Insert completed element matrices into system matrices
  [rows] = get_element_index (n_g, n_n, e_nodes); % eq numbers
  S (rows, rows) = S (rows, rows) + S_e  ; % add to system stiffness
  M (rows, rows) = M (rows, rows) + M_e  ; % add to system mass
  C (rows) = C (rows) + C_e              ; % add to sys force/couples

end % for each j element in mesh   <<==== <<==== <<==== <<==== <<====

%    ALLOCATE STORAGE FOR OPTIONAL REACTION RECOVERY          
if ( EBC_count > 0 )                        ; % reactions occur
  [EBC_row, EBC_col] = save_reaction_matrices (EBC_flag, S, C);
  save_resultant_load_vectors (n_g, C)
end % if essential BC exist (almost always true)

%         ENFORCE ESSENTIAL BOUNDARY CONDITIONS
if ( EBC_count > 0 )                        ; % reactions occur
  [S, C] = enforce_essential_BC (EBC_flag, EBC_value, S, C);
end % if essential BC exist (almost always true)

%             COMPUTE SOLUTION & SAVE
T = S \ C                 ; % Compute displacements & rotations
list_save_beam_displacements (n_g, n_m, T)   ; % save and print

%                GRAPHS
if ( BEF == 1 ) % plot foundation pressure
  cubic_beam_pressure (F_k) % graph
  pause (5) % display for 5 seconds
end % if
cubic_beam_deflection (1) % plot the beam deflection
% cubic_beam_deflection (2) % plot the beam slope

%          OPTIONAL REACTION RECOVERY & SAVE
if ( EBC_count > 0 )                    ; % reactions exist ?
  [EBC_react] = recover_reactions_print_save (n_g, n_d, ...
                EBC_flag, EBC_row, EBC_col, T); % reaction to EBC
end % if EBC exist

%       POST-PROCESS ELEMENT REACTIONS (MEMBER FORCES)

fprintf ('\nIndividual Element Load and Reaction Summaries: \n')
fprintf ('(F_1, M_1, F_2, M_2) \n')
for j = 1:n_e ; % loop over elements ====>> ====>> ====>> ====>> ====>>
  e_nodes = nodes (j, 1:n_n)                      ; % connectivity
  [rows] = get_element_index (n_g, n_n, e_nodes)  ; % eq numbers
  T_e (1:n_i) = T(rows)         ; % gather element displacements
  
%               SET ELEMENT PROPERTIES & GEOMETRY
  L_e = x(e_nodes(2)) - x(e_nodes(1)) ; % beam length
  Line_e (1:2) = 0 ;
  if ( n_prop > 1 )
    E            = msh_properties (j, 1)   ; % material modulus
    I            = msh_properties (j, 2)   ; % section inertia
    Line_e (1:2) = msh_properties (j, 3:4) ; % trapezoidal line load
    k_f          = msh_properties (j, 5)   ; % foundation stiffness
  end % if

%      ELEMENT CONDUCTION AND INTERNAL SOURCE MATRICES 
%  Constant cubic element & foundation stiffnesses & mass matrices
  S_e = (E*I/L_e^3)*[ 12,    6*L_e,   -12,    6*L_e     ;
                      6*L_e, 4*L_e^2, -6*L_e, 2*L_e^2   ;
                     -12,   -6*L_e,    12,   -6*L_e     ;
                      6*L_e, 2*L_e^2, -6*L_e, 4*L_e^2 ] ; % stiffness

  NtN = L_e*[ 156,     22*L_e,   54,     -13*L_e  ;
              22*L_e,  4*L_e^2,  13*L_e, -3*L_e^2 ;
              54,      13*L_e,   156,    -22*L_e  ;
             -13*L_e, -3*L_e^2, -22*L_e,  4*L_e^2 ]/420; 

%  Map line load to node forces & moments; C_e = p_To_F * Line_e 
  C_e = zeros (n_i, 1)           ; % clear arrays
  if ( any (Line_e) )            ; % then form forcing vector
    p_To_F = L_e * [ 21,    9     ; 
                     3*L_e, 2*L_e ; 
                     9,     21    ; 
                    -2*L_e -3*L_e ] / 60 ; % cubic H, linear Line
    % 4 x 2 * 2 x 1 = 4 x 1 result 
    C_e = p_To_F (1:n_i, 1:n_n) * Line_e' ; % force moment @ nodes
  end % if or set up resultant node loads for line load

%  Assign any point force/couple to the first element with that node
  if ( j == 1 ) % first element only gets sources from both nodes
    C_e = C_e + C_react (rows)  ;
  else
    C_e (3:4) = C_e (3:4) + C_react (rows(3:4)) ; % cubic beam ONLY !!
  end % if
  fprintf ('\nGiven Resultant Loading on Member %g, \n', j)
  disp (C_e')

% Finally, get the reactions 
  C_m = S_e * T_e' - C_e  ;
  if ( BEF ) % add foundation resultants
    C_f = -k_f * NtN * T_e'  ;
    fprintf ('Resultant Foundation Loading on Member %g, \n', j) 
    disp (C_f') 
    C_m = C_m - C_f ;
  end % if
  fprintf ('Net Resultant End Reactions on Member %g, \n', j) 
  disp (C_m') 

end % for each j element in mesh   <<==== <<==== <<==== <<==== <<====

% End finite element calculations. 
% See /home/mech517/public_html/Matlab_Plots for graphic options
% http://www.owlnet.rice.edu/~mech517/help_plot.html for help
% end of Beam_on_Winkler           

%  +++++++++++++ functions in alphabetical order +++++++++++++++++

function  cubic_beam_deflection (i_p) 
% Copyright 2004, J.E. Akin.  All rights reserved.
%  ------------------------------------------------------
%  Matlab graph of i_p-th component value of beam
%  1=deflection, 2=slope  
%  ------------------------------------------------------
% c_x    = x coordinates of nod_per_el line polygon
% c_y    = y coordinates of nod_per_el line polygon
% msh_typ_nodes   = connectivity list for elements
% loop   = corners for nod_per_el line polygon
% nod_per_el = Nodes per element
% np     = Number of Points
% nt     = Number of elements
  pre_e  = 0 ; % Element items before connectivity list
  pre_p = 1; % Nodal items before coordinates
% msh_bc_xyz = Nodal coordinates (with preceeding data)
% t_x    = x coordinates of nod_per_el corners
% t_y    = y coordinates of nod_per_el corners
  if ( nargin == 0 )
    i_p = 1 ; 
  end % if
  Pts_wide = 2 ;  % fat lines

%    Read coordinate file and connectivity file
%    integer bc code, real xy pairs for np points (pre_p = 1)
  load msh_bc_xyz.tmp ;  

%    Set control data: number of points
  np = size (msh_bc_xyz,1) ;   % number of nodal points  
  %b fprintf ('Read %g mesh coordinate pairs \n', np)
  ns = size (msh_bc_xyz,2) - pre_p ;   % space dimension
  if ( ns > 1 )
    fprintf ('Not 1D: will use x-coordinate only \n')
  end % if not 2D data
  
%    Set control data: number elements
  load  msh_typ_nodes.tmp ; % nod_per_el nodes per element
  nt = size (msh_typ_nodes,1) ;          % number of elements in mesh 
  nod_per_el = size (msh_typ_nodes,2) - pre_e -1 ; % nodes per elem
  %b fprintf ('Read %g elements connections \n', nt)

  load node_results.tmp
  nr = size (node_results, 1);
  if ( nr == 0 )
    error ('Error missing file node_results.tmp')
  end % if error
  max_p = size (node_results, 2) ; % number of columns

  H (2)       = 0. ;
  HC1 (4)     = 0. ;
  DHC1 (4)     = 0. ;
  x (np)      = 0.       ;       % pre-allocate array x
  
  t_nodes (nod_per_el) = 0 ;    % Optional pre-allocation
  t_x (nod_per_el) = 0 ;        % Optional pre-allocation
  t_y (nod_per_el) = 0 ;        % Optional pre-allocation
  c_x (nod_per_el + 1) = 0 ;    % Optional pre-allocation
  c_y (nod_per_el + 1) = 0 ;    % Optional pre-allocation
  loop = [1:nod_per_el] ; % default to sequential order

  % msh_bc_xyz has: pre_p items then: x, y
  x = msh_bc_xyz (1:np, (pre_p+1)) ; % extract x column
  xmax = max (x)      ; xmin = min (x) ;
  y  = node_results(:, 1) ;
  dy = node_results(:, 2) ;

   clf                               % clear graphics
   hold on                           % hold image for plots
   xlabel (['X, Node at 45 deg (', int2str(nod_per_el), ...
     ' per element), Element at 90 deg'])

%    Loop over all elements
  ymax = max (node_results(:, i_p)) ;
  ymin = min (node_results(:, i_p)) ;
  el_max  = ymax ; el_min  = ymin  ;

  for it = 1:nt  ;

%      Extract element connectivity 
    t_nodes = msh_typ_nodes (it, (pre_e+2):(nod_per_el+pre_e+1)); 

%      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
      A = t_x(2) - t_x(1) ;        % element length
      t_y  = y (t_nodes) ;       % y at those nodes, only
      t_dy  = dy (t_nodes) ;     % dy at those nodes, only
      D (1:2:4) = t_y ; D (2:2:4) = t_dy ;

%      Loop over local points on the cubic polynomial element
      n_poly = ceil ( 95 / nt) ;
      for k = 1: (n_poly + 1) % points in parametric space
%          get element parametric interpolation functions
        R     = (k - 1)/n_poly ;   % on 0 to 1
        X     = 2*R - 1 ;          % on -1 to 1
%       H = ELEMENT CUBIC SHAPE FUNCTIONS
%       X = LOCAL COORDINATE OF POINT,    -1 TO +1
%       LOCAL NODE COORD. ARE -1,+1    1------------2
        H (1) = 0.5*(1 - X) ;
        H (2) = 0.5*(1 + X) ;
        x_el (k) = H * t_x ;               % true X value
        HC1(1) = (2 - 3*X + X^3)/4;        
        HC1(2) = (1 - X - X^2 + X^3)*A/8;
        HC1(3) = (2 + 3*X - X^3)/4;    
        HC1(4) = (-1 - X + X^2 + X^3)*A/8;
        y_el (k) = HC1 * D' ;              % true y value
        DHC1(1) = (-3 + 3*X*X) * 0.5 / A;
        DHC1(2) = (-1 - 2*X + 3*X*X) * 0.25 ;
        DHC1(3) = ( 3 - 3*X*X) * 0.5d0 / A;
        DHC1(4) = (-1 + 2*X + 3*X*X) * 0.25;
        dy_el (k) = DHC1 * D' ;            % true dy value
      end % for k

      format short
      plot (x, y, 'ko') % nodal value symbols
      if ( i_p == 1 )
       r_el = y_el ;
      elseif ( i_p == 2 )
       r_el = dy_el ;
      end % if

%    elem max, min values
      [V_X, L_X] = max (r_el)  ;
      [V_N, L_N] = min (r_el);
       if ( V_X > el_max  )
         el_max = V_X ;
       end
       if ( V_N < el_min  )
         el_min = V_N ;
       end
       plot (x_el, r_el, 'b-', 'LineWidth',Pts_wide)

% Plot the element number
       x_bar = sum (t_x' )/nod_per_el ;
       y_bar = mean (r_el) ;
       t_text = sprintf ('   (%g)', it);  % offset # from pt
       text (x_bar, y_bar, t_text) % incline

    end % if zero in connectivity
  end % for all elements

  if ( el_max > ymax )
    fprintf ('Max interior value was %g \n', el_max)
    ymax = el_max ;
  end % if
  if ( el_min < ymin )
    ymin = el_min ;
    fprintf ('Min interior value was %g \n', el_min)
  end % if
  axis ([xmin, xmax, ymin, ymax])   % set axes
  null (1:np) = 0.5*(ymax + ymin) ;

  R_x = max (node_results(:, i_p)) ;
  R_n = min (node_results(:, i_p)) ;
  if ( i_p == 1 )
    ylabel (['Displacement (nodal max = ', ...
              num2str(R_x), ', min = ', num2str(R_n), ')'])
    title(['Beam Displacement: ', int2str(nt),' Elements, ', ...
           int2str(np), ' Nodes, (', int2str(nod_per_el), ...
           ' per Element)'])
  elseif ( i_p == 2 )
    ylabel (['Slope  (nodal max = ', ...
              num2str(R_x), ', min = ', num2str(R_n), ')'])
    title(['Beam Slope: ', int2str(nt),' Elements, ', ...
           int2str(np), ' Nodes, (', int2str(nod_per_el), ...
           ' per Element)'])
  end % if

% plot node points on axis
  inc_p = 1 ;
  if ( inc_p > 0 )
    for i = 1:np
      t_text = sprintf ('   %g', i);  % offset # from pt
      text (x(i), null(i), t_text, 'Rotation', 45) % incline
    end % for all
    plot (x, null, 'k*')
  end % if

  grid
  print -dpng deflection_plot
  hold off
  fprintf ('Saved deflection plot as deflection_plot.png \n')
% end of cubic_beam_deflection

function  cubic_beam_pressure   (F_k) 
% Copyright 2008, J.E. Akin.  All rights reserved.
%  ------------------------------------------------------
%  Matlab graph of BEF pressure distribution    m
%  ------------------------------------------------------
% c_x    = x coordinates of nod_per_el line polygon
% c_y    = y coordinates of nod_per_el line polygon
% msh_typ_nodes   = connectivity list for elements
% loop   = corners for nod_per_el line polygon
% nod_per_el = Nodes per element
% np     = Number of Points
% nt     = Number of elements
  pre_e  = 0 ; % Element items before connectivity list
  pre_p = 1; % Nodal items before coordinates
% msh_bc_xyz = Nodal coordinates (with preceeding data)
% t_x    = x coordinates of nod_per_el corners
% t_y    = y coordinates of nod_per_el corners
  i_p = 1 ; % temp
  n_fk = size (F_k, 2) ; % number of different foundationS

%    Read coordinate file and connectivity file
%    integer bc code, real xy pairs for np points (pre_p = 1)
  load msh_bc_xyz.tmp ;  

%    Set control data: number of points
  np = size (msh_bc_xyz,1) ;   % number of nodal points  
  %b fprintf ('Read %g mesh coordinate pairs \n', np)
  ns = size (msh_bc_xyz,2) - pre_p ;   % space dimension
  if ( ns > 1 )
    fprintf ('Not 1D: will use x-coordinate only \n')
  end % if not 2D data
  null (1:np) = 0      ;
  Pts_wide = 2 ;  % fat lines
  
%    Set control data: number elements
  load  msh_typ_nodes.tmp ; % nod_per_el nodes per element
  nt = size (msh_typ_nodes,1) ;          % number of elements in mesh 
  nod_per_el = size (msh_typ_nodes,2) - pre_e -1 ; % nodes per elem
  %b fprintf ('Read %g elements connections \n', nt)

  load node_results.tmp
  nr = size (node_results, 1);
  if ( nr == 0 )
    error ('Error missing file node_results.tmp')
  end % if error
  max_p = size (node_results, 2) ; % number of columns

  H (2)       = 0. ;
  HC1 (4)     = 0. ;
  DHC1 (4)     = 0. ;
  x (np)      = 0.       ;       % pre-allocate array x
  
  t_nodes (nod_per_el) = 0 ;    % Optional pre-allocation
  t_x (nod_per_el) = 0 ;        % Optional pre-allocation
  t_y (nod_per_el) = 0 ;        % Optional pre-allocation
  c_x (nod_per_el + 1) = 0 ;    % Optional pre-allocation
  c_y (nod_per_el + 1) = 0 ;    % Optional pre-allocation

  loop = [1:nod_per_el] ;
  % msh_bc_xyz has: pre_p items then: x, y
  x = msh_bc_xyz (1:np, (pre_p+1)) ; % extract x column
  xmax = max (x)      ; xmin = min (x) ;
  y  = node_results(:, 1) ;
  dy = node_results(:, 2) ;

   clf                               % clear graphics
   hold on                           % hold image for plots
   xlabel (['X, Node at 45 deg (', int2str(nod_per_el), ...
     ' per element), Element at 90 deg'])

%    Loop over all elements
  Pmax = -max (node_results(:, i_p)) * max (F_k) ;
  Pmin = -min (node_results(:, i_p)) * max (F_k) ;
  el_max  = Pmax ; el_min= Pmin  ;

%    Loop over all elements
  for it = 1:nt  ;
%      Element foundation stiffness, if any
    if ( n_fk == 1 ) % then homogeneous
      k_f = F_k (1)  ;
    else
      k_f = F_k (it)  ;
    end % if 

%      Extract element connectivity 
    t_nodes = msh_typ_nodes (it, (pre_e+2):(nod_per_el+pre_e+1)); 

%      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
      A = t_x(2) - t_x(1) ;        % element length
      t_y  = y (t_nodes) ;       % y at those nodes, only
      t_dy  = dy (t_nodes) ;     % dy at those nodes, only
      D (1:2:4) = t_y ; D (2:2:4) = t_dy ;

%      Loop over local points on the cubic polynomial element
      n_poly = ceil ( 95 / nt) ;
      for k = 1: (n_poly + 1) % points in parametric space
%          get element parametric interpolation functions
        R     = (k - 1)/n_poly ;   % on 0 to 1
        X     = 2*R - 1 ;          % on -1 to 1
%       H = ELEMENT CUBIC SHAPE FUNCTIONS
%       X = LOCAL COORDINATE OF POINT,    -1 TO +1
%       LOCAL NODE COORD. ARE -1,+1    1------------2
        H (1) = 0.5*(1 - X) ;
        H (2) = 0.5*(1 + X) ;
        x_el (k) = H * t_x ;               % true X value
        HC1(1) = (2 - 3*X + X^3)/4;        
        HC1(2) = (1 - X - X^2 + X^3)*A/8;
        HC1(3) = (2 + 3*X - X^3)/4;    
        HC1(4) = (-1 - X + X^2 + X^3)*A/8;
        y_el (k) = HC1 * D' ;              % true y value
      end % for k
      p_el = -y_el * k_f ;       % true pressure
      plot(x_el, p_el, 'b-', 'LineWidth',Pts_wide)

% Plot the element number
       x_bar = mean (x_el) ;
       y_bar = mean (p_el) ;
       t_text = sprintf ('   (%g)', it);  % offset # from pt
       text (x_bar, y_bar, t_text) % incline

      format short
      max_el = max (p_el)  ;
      min_el = min (p_el)  ;
      if ( max_el > el_max )
        el_max = max_el ;
      end % if
      if ( min_el < el_min )
        el_min = min_el ;
      end % if
    end % if has non-zero nodes
  end % for all elements
  fprintf ('Max interior pressure was %g \n', el_max)
  fprintf ('Min interior pressure was %g \n', el_min)

% plot node points on axis
  inc_p = 1 ;
  if ( inc_p > 0 )
    for i = 1:np
      t_text = sprintf ('   %g', i);  % offset # from pt
      text (x(i), null(i), t_text, 'Rotation', 45) % incline
    end % for all
    plot (x, null, 'k*')
  end % if

  ylabel (['Pressure (nodal max = ', ...
            num2str(Pmax), ', min = ', num2str(Pmin), ')'])
  title(['Foundation Pressure: ', int2str(nt),' Elements, ', ...
         int2str(np), ' Nodes, (', int2str(nod_per_el), ...
         ' per Element)'])
  grid
  print -dpng pressure_plot
  hold off
  fprintf ('Saved pressure plot as pressure_plot.png \n')
% end of cubic_beam_pressure   

function [S, C] = enforce_essential_BC (EBC_flag, EBC_value, S, C)
% modify system linear eqs for essential boundary conditions
% (by trick to avoid matrix partitions, loses reaction data)
 n_d = size (C, 1)                        ; % number of DOF eqs
if ( size (EBC_flag, 2) > 1 ) ; % change to vector copy
  flag_EBC  = reshape ( EBC_flag',  1, n_d) ;
  value_EBC = reshape ( EBC_value', 1, n_d) ;
else
  flag_EBC  = EBC_flag ;
  value_EBC = EBC_value ;
end % if
 for j = 1:n_d                % check all DOF for essential BC 
   if ( flag_EBC (j) )                        % then EBC here
%    Carry known columns*EBC to RHS. Zero that column and row.
%    Insert EBC identity, 1*EBC_dof = EBC_value.
     EBC = value_EBC (j)             ; % recover EBC value
     C (:) = C (:) - EBC * S (:, j) ; % carry known column to RHS
     S (:, j) = 0 ; S (j, :) = 0    ; % clear, restore symmetry
     S (j, j) = 1 ; C (j) = EBC     ; % insert identity into row
   end % if EBC for this DOF
 end % for over all j-th DOF
% end enforce_essential_BC (EBC_flag, EBC_value, S, C)

function [C] = get_and_add_point_sources (n_g, n_m, C)
  load msh_load_pt.tmp            ; % node, DOF, value (eq. number)
  n_u = size(msh_load_pt, 1)      ; % number of point sources
  if ( n_u < 1 )                  ; % missing data
    error ('No load_pt data in msh_load_pt.tmp')
  end % if user error
  fprintf ('\nRead %g point sources. \n', n_u)
  fprintf ('(Echo of file msh_load_pt.tmp) \n')
  fprintf ('Node, DOF (1=force, 2=couple),  Source_value \n')
  for j = 1:n_u                    ; % non-zero Neumann pts
    node  = msh_load_pt (j, 1)         ; % global node number
    DOF   = msh_load_pt (j, 2)         ; % local DOF number
    value = msh_load_pt (j, 3)         ; % point source value
    fprintf ('%g  %g  %g \n', node, DOF, value)
    Eq    = n_g * (node - 1) + DOF     ; % row in system matrix
    C (Eq) = C (Eq) + value     ; % add to system column matrix
  end % for each EBC
% end get_and_add_point_sources (n_g, n_m, C)

function [EBC_flag] = get_ebc_flags (n_g, n_m, P)
 EBC_flag = zeros(n_m, n_g)   ; % initialize
 for k = 1:n_m                        ; % loop over all nodes
   if ( P(k) > 0 )                    ; % at least one EBC here
     [flags] = unpack_pt_flags (n_g, k, P(k)) ; % unpacking
     EBC_flag (k, 1:n_g) = flags (1:n_g)      ; % populate array
   end % if EBC at node k
 end % for loop over all nodes
% end get_ebc_flags

function [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag)
 EBC_value = zeros(n_m, n_g)  ; % initialize to zero
  load msh_ebc.tmp            ; % node, DOF, value (eq. number)
  n_c = size(msh_ebc, 1)      ; % number of constraints
  fprintf ('\nApplied Displacement Boundary Conditions: %g \n', n_c)
  fprintf ('(Echo of file load msh_ebc.tmp) \n')               
  fprintf ('Node, DOF (1=displacement, 2=slope), Value. \n')
  disp(msh_ebc)                    ; % echo input

  for j = 1:n_c                    ; % loop over ebc inputs
    node  = round (msh_ebc (j, 1)) ; % node in mesh
    DOF   = round (msh_ebc (j, 2)) ; % DOF # at node
    value = msh_ebc (j, 3)         ; % EBC value   
    % Eq  = n_g * (node - 1) + DOF ; % row in system matrix
    EBC_value (node, DOF) = value  ; % insert value in array
    if ( EBC_flag (node, DOF) == 0 ) % check data consistency
      fprintf ('WARNING: EBC but no flag at node %g & DOF %g. \n', ...
                node, DOF)
      % EBC_flag (node, DOF) = 1; % try to recover from data error
    end % if common user error
  end % for each EBC

  EBC_count = sum (sum ( EBC_flag > 0 )) ; % check input data
  if ( EBC_count ~= n_c )                ; % probable user error
    fprintf ('\nWARNING: mismatch in bc_flag count & msh_ebc.tmp \n')
  end % if user error
% end get_ebc_values

function [rows] = get_element_index (n_g, n_n, e_nodes)
% calculate system DOF numbers of element, for gather, scatter
rows = zeros (1, n_g*n_n)           ; % allow for node = 0
for k = 1:n_n                       ; % loop over element nodes
  global_node = round (e_nodes (k)) ; % corresponding sys node
  for i = 1:n_g                     ; % loop over DOF at node
    eq_global  = i + n_g * (global_node - 1) ; % sys DOF, if any
    eq_element = i + n_g * (k           - 1) ; % el DOF number
    if ( eq_global > 0 )                ; % check node=0 trick
      rows (1, eq_element) = eq_global  ; % valid DOF > 0
    end % if allow for omitted nodes
  end % for DOF i                       % end local DOF loop
end % for each element node             % end local node loop
% end get_element_index 

function [n_e, n_n, n_t, el_type, nodes] = get_mesh_elements () ;
% input file controls (for various data generators) 

load  msh_typ_nodes.tmp        ; % el_type, connectivity list (3) 
n_e = size (msh_typ_nodes,1)   ; % number of elements 
if ( n_e == 0 )                ; % data file missing 
  error ('Error missing file msh_typ_nodes.tmp') 
end % if error 
n_n = size (msh_typ_nodes,2) - 1 ; % nodes per element       
fprintf ('(Echo of file msh_typ_nodes.tmp) \n')
fprintf ('Read %g elements with (ignored) type & %g nodes each. \n', ...
          n_e, n_n) 
el_type = round (msh_typ_nodes(:, 1)); % el type number >= 1 
n_t     = max(el_type)                   ; % number of element types 
%b fprintf ('Maximum number of element types = %g. \n', n_t) 
nodes (1:n_e, 1:n_n) = msh_typ_nodes (1:n_e, 2:1+n_n); 
disp(msh_typ_nodes (:, 1:1+n_n)) % echo data
% end get_mesh_elements 

function [n_m, n_s, P, x, y, z] = get_mesh_nodes () ;
% input file controls (for various data generators) 
 
%       READ MESH AND EBC_FLAG INPUT DATA 
% specific problem data from MODEL data files (sequential) 
load msh_bc_xyz.tmp       ; % bc_flag, x-, y-, z-coords 
n_m = size (msh_bc_xyz,1) ; % number of nodal points in mesh 
if ( n_m == 0 )           ; % data missing ! 
  error ('Error missing file msh_bc_xyz.tmp') 
end % if error 
n_s = size (msh_bc_xyz,2) - 1 ; % number of space dimensions 
fprintf ('Read %g nodes. \n', n_m)  
fprintf ('(Echo of file msh_bc_xyz.tmp) \n')
fprintf ('bc_flag, x-coordinate \n')
msh_bc_xyz (:, 1)= round (msh_bc_xyz (:, 1));
P = msh_bc_xyz (1:n_m, 1) ; % integer Packed BC flag 
x = msh_bc_xyz (1:n_m, 2) ; % extract x column  
y (1:n_m, 1) = 0.0 ;  z (1:n_m, 1) = 0.0 ; % default to zero
 
if (n_s > 1 )                                 ; % check 2D or 3D
  y = msh_bc_xyz (1:n_m, 3)           ; % extract y column  
end % if 2D or 3D
if ( n_s == 3 )                               ; % check 3D
  z = msh_bc_xyz (1:n_m, 4)           ; % extract z column 
end % if 3D 
%b disp(msh_bc_xyz (:, 1:1+n_s)) ; % echo data
for j = 1:n_m
  fprintf (' %2.2i  %g         \n', P(j), x(j)             )
%bfprintf (' %2.2i  %g  %g     \n', P(j), x(j), y(j)       )
%bfprintf (' %2.2i  %g  %g  %g \n', P(j), x(j), y(j), z(j) )
end % for
% end get_mesh_nodes 

function list_save_beam_displacements (n_g, n_m, T) 
 fprintf('\nCalculated Displacements: \n')
 fprintf('Node,  Y_displacement,  Z_rotation at %g nodes \n', n_m)
 T_matrix = reshape (T, n_g, n_m)'            ; % pretty shape       
 % save results (displacements) to MODEL file: node_results.tmp
 fid = fopen('node_results.tmp', 'w') ; % open for writing
 for j = 1:n_m                        ; % node loop, save displ
   fprintf (fid, '%g %g \n',  T_matrix (j, 1:n_g)) ; % to file
   fprintf ('%g %g %g \n', j, T_matrix (j, 1:n_g)) ; % to screen
 end % for j DOF
% end list_save_beam_displacements (n_g, n_m, T) 

function list_save_displacements_results (n_g, n_m, T)
  fprintf ('\n') ;
  fprintf('X_disp      Y_disp      Z_disp   at %g nodes \n', n_m)
  T_matrix = reshape (T, n_g, n_m)'            ; % pretty shape       
  disp (T_matrix)                              ; % print displacements
  % save results (displacements) to MODEL file: node_results.tmp
  fid = fopen('node_results.tmp', 'w')          ; % open for writing
  for j = 1:n_m                                 ; % save displacements
    if ( n_g == 1 )
      fprintf (fid, '%g                \n', T_matrix (j, 1:n_g)) ;
    elseif ( n_g == 2 )
      fprintf (fid, '%g %g             \n', T_matrix (j, 1:n_g)) ;
    elseif ( n_g == 3 )
      fprintf (fid, '%g %g %g          \n', T_matrix (j, 1:n_g)) ;
    elseif ( n_g == 4 )
      fprintf (fid, '%g %g %g %g       \n', T_matrix (j, 1:n_g)) ;
    elseif ( n_g == 5 )
      fprintf (fid, '%g %g %g %g %g    \n', T_matrix (j, 1:n_g)) ;
    elseif ( n_g == 6 )
      fprintf (fid, '%g %g %g %g %g %g \n', T_matrix (j, 1:n_g)) ;
    else 
      error ('reformat list_save_displacements_results for n_g > 6.')
    end % if 
  end % for j DOF
% end list_save_displacements_results (T)

function [EBC_react] = recover_reactions_print_save (n_g, n_d, ...
                       EBC_flag, EBC_row, EBC_col, T)
% get EBC reaction values by using rows of S & C (before EBC)
  n_d = size (T, 1)                     ; % number of system DOF
% n_c x 1 = n_c x n_d  * n_d x 1 +  n_c x 1
  EBC_react = EBC_row * T - EBC_col     ; % matrix reactions (+-)
  % save reactions (forces) to MODEL file: node_reaction.tmp
  fprintf ('\nRecovered Reactions at Displacement or Slope BC: %g \n', ...
            sum (sum (EBC_flag > 0)))   ; % header
  fprintf ('Node, DOF (1=force, 2=couple), Value \n')
  fid = fopen('node_reaction.tmp', 'w') ; % open for writing
  if ( size (EBC_flag, 2) > 1 )            ; % change to vector copy
    flag_EBC  = reshape ( EBC_flag',  1, n_d) ; % changed
  else
    flag_EBC  = EBC_flag                      ; % original vector
  end % if
  Totals = zeros (1, n_g)               ; % zero input totals
  kount = 0                             ; % initialize counter
  for j = 1:n_d                         ; % extract all EBC reactions
    if ( flag_EBC(j) )                      ; % then EBC here
  %   Output node_number, component_number, value, equation_number 
      kount = kount + 1                     ; % copy counter
      node  = ceil(j/n_g)                   ; % node at DOF j
      j_g   = j - (node - 1)*n_g            ; % 1 <= j_g <= n_g
      React = EBC_react (kount, 1)          ; % reaction value
      fprintf ( fid, '%g %g %g \n', node, j_g, React);% save
      fprintf ('%g %g %g \n', node, j_g, React);     % print
      Totals (j_g) = Totals (j_g) + React   ; % sum all components
    end % if EBC for this DOF
  end % for over all j-th DOF
  fprintf ('Total force and couple = ') ; disp(Totals) ; % echo total
% end recover_reactions_print_save (EBC_row, EBC_col, T)

function [EBC_row, EBC_col] = save_reaction_matrices (EBC_flag, S, C)
  n_d       = size (C, 1)                  ; % number of system DOF
  EBC_count = sum (sum (EBC_flag))         ; % count EBC & reactions
  EBC_row   = zeros(EBC_count, n_d)        ; % reaction data
  EBC_col   = zeros(EBC_count, 1)          ; % reaction data
  if ( size (EBC_flag, 2) > 1 )            ; % change to vector copy
    flag_EBC  = reshape ( EBC_flag',  1, n_d) ; % changed
  else
    flag_EBC  = EBC_flag                      ; % original vector
  end % if
  kount = 0                                ; % initialize counter
  for j = 1:n_d      % System DOF loop, check for displacement BC
    if ( flag_EBC (j) )                    ; % then EBC here
  %   Save reaction data to be destroyed by EBC solver trick
      kount = kount + 1                    ; % copy counter
      EBC_row(kount, 1:n_d) = S (j, 1:n_d) ; % copy reaction data
      EBC_col(kount, 1)     = C (j)        ; % copy reaction data
    end % if EBC for this DOF
  end % for over all j-th DOF                % end sys DOF loop
% end save_reaction_matrices (S, C, EBC_flag)

function save_resultant_load_vectors (n_g, C)
% save resultant forces to MODEL file: node_resultants.tmp
  n_d = size (C, 1)                     ; % number of system DOF
  fprintf ('\nResultants of all input sources: \n')
  fprintf ('Node, DOF (1=force, 2=couple), Value \n')
  fid = fopen('node_resultant.tmp', 'w'); % open for writing
  Totals = zeros (1, n_g)               ; % zero input totals
  for j = 1:n_d                         ; % extract all resultants
    if ( C (j) ~= 0. )                      ; % then source here
  %   Output node_number, component_number, value, equation_number 
      node  = ceil(j/n_g)                   ; % node at DOF j
      j_g   = j - (node - 1)*n_g            ; % 1 <= j_g <= n_g
      value = C (j)                         ; % resultant value
      fprintf ( fid, '%g %g %g %g \n', node, j_g, value, j);% save
      fprintf ('%g %g %g \n', node, j_g, value);     % print
      Totals (j_g) = Totals (j_g) + value   ; % sum all inputs
    end % if non-zero for this DOF
  end % for over all j-th DOF
  fprintf ('Totals = ') ; disp(Totals)      ; % echo totals
% end save_resultant_load_vectors (n_g, n_m, C)

function [flags] = unpack_pt_flags (n_g, N, flag)
% unpack n_g integer flags from the n_g digit flag at node N
% integer flag contains (left to right) f_1 f_2 ... f_n_g
 full = flag                         ; % copy integer
 check = 0                           ; % validate input
 for Left2Right = 1:n_g              ; % loop over local DOF at k
   Right2Left = n_g + 1 - Left2Right ; % reverse direction
   work = floor (full / 10)          ; % work item
   keep = full - work * 10           ; % work item
   flags (Right2Left) = keep         ; % insert into array
   full = work                       ; % work item
   check = check + keep * 10^(Left2Right - 1) ; % validate
 end % for each local DOF
 if ( flag > check )                 ; % check for likely error
   fprintf ('WARNING: bc flag likely reversed at node %g. \n', N)
 end % if likely user error
% end unpack_pt_flags
 
%=================== Running gives =====================================
%  >> Cubic_Beam_on_Foundation(1)
%  Read 5 nodes with bc_flag & 1 coordinates. 
%  (Echo of file msh_bc_xyz.tmp) 
%      10   -60
%       0   -30
%       0    -8
%      10    40
%       0    60
%  
%  Note: expecting 2 displacement BC values. 
%
%  Read 4 elements with (ignored) type & 2 nodes each. 
%  (Echo of file msh_typ_nodes.tmp) 
%       1     1     2
%       1     2     3
%       1     3     4
%       1     4     5
%  
%  Applied Displacement Boundary Conditions: 2 
%  (Echo of file load msh_ebc.tmp) 
%  Node, DOF (1=displacement, 2=slope), Value. 
%       1     1     0
%       4     1     0
%  
%  Read 1 point sources. 
%  (Echo of file msh_load_pt.tmp) 
%  Node, DOF (1=force, 2=couple),  Source_value 
%  2  1  -5000 
%  
%  (Echoing file msh_properties.tmp)
%  Properties for element 1 
%  Elastity modulus (N/m^2)     = 1.5e+06 
%  Moment of inertia (m^4)      = 426.7 
%  Line Load (N/m)              = [ 0 0 ] 
%  Foundation stiffness (N/m^2) = 0 
%  Mass per unit length (kg/m)  = 0  
%  
%  Properties for element 2 
%  Elastity modulus (N/m^2)     = 1.5e+06 
%  Moment of inertia (m^4)      = 426.7 
%  Line Load (N/m)              = [ 0 0 ] 
%  Foundation stiffness (N/m^2) = 0 
%  Mass per unit length (kg/m)  = 0  
%  
%  Properties for element 3 
%  Elastity modulus (N/m^2)     = 1.5e+06 
%  Moment of inertia (m^4)      = 426.7 
%  Line Load (N/m)              = [ -100 -100 ] 
%  Foundation stiffness (N/m^2) = 0 
%  Mass per unit length (kg/m)  = 0  
%  
%  Properties for element 4 
%  Elastity modulus (N/m^2)     = 1.5e+06 
%  Moment of inertia (m^4)      = 426.7 
%  Line Load (N/m)              = [ 0 0 ] 
%  Foundation stiffness (N/m^2) = 0 
%  Mass per unit length (kg/m)  = 0  
%  
%  Resultants of all input sources: 
%  Node, DOF (1=force, 2=couple), Value 
%  2 1 -5000 
%  3 1 -2400 
%  3 2 -19200 
%  4 1 -2400 
%  4 2 19200 
%  Totals =        -9800           0
%  
%  Calculated Displacements: 
%  Node,  Y_displacement,  Z_rotation at 5 nodes 
%  1      0                -0.00730227 
%  2      -0.186361        -0.00403159 
%  3      -0.223254        0.000633838 
%  4      0                0.00701974 
%  5      0.140395         0.00701974 
%  Saved deflection plot as deflection_plot.png 
%  
%  Recovered Reactions at Displacement or Slope BC: 2 
%  Node, DOF (1=force, 2=couple), Value 
%  1 1 4652 
%  4 1 5148 
%  Total force and couple =    1.0e+03 *
%      9.8000         0
%  
%  Individual Element Load and Reaction Summaries: 
%  (F_1, M_1, F_2, M_2) 
%  Given Resultant Loading on Member 1, 
%             0           0       -5000           0
%  Net Resultant End Reactions on Member 1, 
%     1.0e+05 *[ 0.0465         0    0.0035    1.3956]
%  
%  Given Resultant Loading on Member 2, 
%       0     0     0     0
%  Net Resultant End Reactions on Member 2, 
%     1.0e+05 *[ -0.0035   -1.3956    0.0035    1.3190]
%  
%  Given Resultant Loading on Member 3, 
%         -2400      -19200       -2400       19200
%  Net Resultant End Reactions on Member 3, 
%     1.0e+05 *[ -0.0035   -1.3190    0.0515   -0.0000]
%  
%  Given Resultant Loading on Member 4, 
%       0     0     0     0
%  Net Resultant End Reactions on Member 4, 
%       0     0     0     0
%  >> quit