function [] = Frame_2D_cubic ()             % Revised 3/13/20 
% Plane Frame Analysis: using classic cubic beam (L2_C1)
% and linear bar (L2_C0) with loads & couples, line load
%             (and input data error checking)
%
%    Note: see sample application after the main program
% Node DOFs: 1=axial disp, 2=transverse disp, 3=z-rotation
% Element Properties: 
%   Area (m^2)                  
%   Elastity modulus (N/m^2)    
%   Moment of inertia (m^4)     
%   Normal line load (N/m) at each element node
%  (Positive to the left from node 1 to node 2)
%.............................................................
%     Mech 417, Finite Element Analysis, Rice University
%     Copyright J.E. Akin 2020. All rights reserved.
%.............................................................
% NOTE: input text files begin with msh_ and end with .txt
% see file FEA_notation.txt for variable descriptions

% echo user information about this specific application
echo_user_remarks

% =====================================================================
%                             PHASE 1
%           *** SET OR READ USER PROBLEM DEPENDENT CONTROLS ***
% =====================================================================

%         Application and element dependent controls
n_g = 3 ; % number of DOF per node (axial_d, transverse_d, rotation)
N_N = 2 ; % number of EXPECTED nodes per element
n_p = 1 ; % dimension of parametric space
n_q = 0 ; % number of quadrature points required
n_r = 1 ; % number of rows in B_e matrix
N_S = 2 ; % number of EXPECTED physical space dimensions

%   set user logic choices (or add read function)
echo_p   = 1  ; % turn on(1)/off(0) listing of properties
el_react = 1  ; % turn on(1)/off(0) listing of element reactions
load_pt  = 0  ; % turn on(1)/off(0) point force or couple(can be zero)
post     = 0  ; % turn on(1)/off(0) post-processing option
show_bc  = 0  ; % turn on(1)/off(0) list assembled S and c, after EBC
show_e   = 0  ; % turn on(>0)/off(0) list some S_e, and c_e
show_s   = 0  ; % turn on(1)/off(0) listing of assembled S and c
stiff_pt = 0  ; % turn on(1)/off(0) point spring stiffness (can be zero)

%    Echo all default or input controls
fprintf ('User control variables \n')
fprintf ('n_g number of DOF per node        = %i \n', n_g)
fprintf ('n_n number of nodes per element   = %i \n', N_N)
fprintf ('n_p dimension of parametric space = %i \n', n_p)
fprintf ('n_q number of quadrature points   = %i \n', n_q)
fprintf ('n_r number of rows in B_e         = %i \n', n_r)
fprintf ('n_s dimension of physical space   = %i \n', N_S)
fprintf ('\n')

%    Echo current user logic flags
fprintf ('User logic flags: turn on(1)/off(0) \n')
fprintf ('echo_p   listing of properties         = %i \n', echo_p  )
fprintf ('el_react listing of element reactions  = %i \n', el_react)
fprintf ('load_pt  point force and/or moment     = %i \n', load_pt )
fprintf ('post     post-processing option        = %i \n', post    )
fprintf ('show_bc  debug list S & c, after EBC   = %i \n', show_bc )
fprintf ('show_e   list some S_e, and c_e        = %i \n', show_e  )
fprintf ('show_s   list assembled S and c        = %i \n', show_s  )
fprintf ('stiff_pt point stiffness (can be zero) = %i \n', stiff_pt)

if ( el_react > 0 ) % then post-processing is required
  post = 1 ; % then post-processing is required
end % if

% =====================================================================
%                               PHASE 2 
%              *** BEGIN DATA INPUT AND CHECKING PHASE ***
% =====================================================================
fprintf ('\n')
fprintf ('Begin Phase 2, input mesh, properties, sources and EBC \n')

%    Read mesh nodal input data files 
%    (from msh_bc_xyz.txt, with one line per node)

% n_m  = number of nodal points in mesh
% n_s  = number of space dimensions
% PD   = nodal packed digit essential BC flag (n_g digits 0 or 1)
%   NOTE: PD is 000-free, 100-u given, 010-v given, 001-slope given, 
%         110-both u & v given (a pin), 111-fixed joint, etc
% x    = mesh x-coordinate values
% y    = mesh y-coordinate values, if present, etc. for z
[n_m, n_s, PD, x, y, z] = get_mesh_nodes (n_g) ; % bc_flags, coordinates

%   Compare actual data against controls
if ( n_s > N_S ) % then not 1D or 2D, possible control error
  fprintf ('WARNING: Read %i coordinates, but expected %i \n', n_s, N_S)
end % if
    
%   *** Read mesh element types and nodal connections ***
%   (from msh_typ_nodes.txt, with one line per element)

% n_e     = number of elements in the mesh
% n_n     = maximum number of nodes per element
% el_type = type number of every element, usually 1
% nodes   = node connection list for each element
% n_t     = number of different element types
[n_e, n_n, n_t, el_type, nodes] = get_mesh_elements  ;

%   Compare actual data against controls
if ( n_n ~= N_N ) % number of expected nodes
  fprintf ('WARNING Read %i nodes per element but expected %i \n', n_n, N_N)
end % if inconsistent inputs

%B addpath /clear/www/htdocs/mech517/Akin_FEA_Lib
%B shrink_2d_mesh (x, y, nodes)             ; % shrink plot the frame shape

%    Read any material properties for the elements
%    (from msh_properties.txt, one line per element or type)

% n_mats = number of materials, = n_e, or n_t if homogeneous
% n_vals = number of property values for each element or type
% msh_properties (n_mats, n_vals) = system properties array
[n_mats, n_vals, msh_properties] = get_mesh_properties (echo_p, n_e, n_t) ;

fprintf ('Above properties are: A, E, I, %i nodal line loads \n', n_n)
fprintf ('For a planar frame \n')
if ( n_mats > n_e )
  error ('ERROR: number of property sets exceeds number of elements')
end % if

%   Count maximum local and system equation numbers  
n_d = n_g * n_m                ; % system degrees of freedom (DOF)
n_i = n_g * n_n                ; % number of DOF per element

%   Pre-allocate and initialize system matrices and vectors.
% Ans       = solution vector
% Body      = body force vector, or volumetric source
% c         = total source vector
% c_pt      = point source part of c for reaction recovery
% EBC_flag  = n_g individual binary flags at each node
% EBC_value = EBC value for each restrained DOF
% el_prop   = element or type properties
% S         = stiffness matrix

Ans       = zeros (n_d, 1)   ; % results
Body      = zeros (n_s, 1)   ; % initialize volumetric source
c         = zeros (n_d, 1)   ; % force or source
c_pt      = zeros (n_d, 1)   ; % point source portion of c
EBC_flag  = zeros (n_m, n_g) ; % initialize all DOF flags
EBC_value = zeros (n_m, n_g) ; % default EBCs to  zero
el_prop   = zeros (n_vals, 1); % initialize element or type properties
S         = zeros (n_d, n_d) ; % generalized stiffness

%   Read point sources & couples data, if any, and insert in c
%   (from file msh_load_pt.txt, one line per point source)
if ( load_pt > 0 )                      ; % need point sources data
  [c_pt] = get_point_sources (n_g, n_m) ; % get point sources
  c      = c + c_pt                     ; % add point sources
end % if any point source expected

%   Read point stiffness data, if any, and add to S diagonal
%   (from file msh_stiff_pt.txt, one line per point stiffness)
if ( stiff_pt > 0 )                        ; % need point stiffness data
  [S] = get_and_add_point_stiffness (n_g, n_m, S); % get point stiffness
end % if any point stiffness expected

%    Extract and count EBC flags from packed integer flag PD
[EBC_flag] = get_ebc_flags (n_g, n_m, PD)     ; % unpack flags

% EBC_count = number of essential (Dirichlet) BC's
EBC_count  = sum( sum ( EBC_flag == 1 ) )    ; % number of EBC
if ( EBC_count == 0 )
  fprintf ('WARNING: No essential boundary condition given. \n')
else
  fprintf ('Note: expecting %i essential BC values. \n', EBC_count)
end % if

%   Read essential boundary condition (EBC) values, if any
%   (from file msh_ebc.txt one line per constraint)
if ( EBC_count > 0 )                          ; % need EBC data
  [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag) ; % read data
else % unusual, but possible
  fprintf ('WARNING: no essential boundary conditions flagged \n')
end % if any EBC data expected

%    *** END OF DATA INPUT AND CHECKING PHASE ***

% ======================================================================
%                              PHASE 3
%       ** GENERATE ELEMENT VALUES AND ASSYMBLE INTO SYSTEM **
% ======================================================================
fprintf ('Begin Phase 3, assemble (scatter) element arrays \n')
if ( show_e > 0 ) % then list some element matrices
  fprintf ('Displaying the first %i sets of element matrices \n', show_e)
end % if arrays to display

%     Open unit needed for post-processing operations
if ( post == 1 ) % open post-processing element unit
  el_unit = fopen ('el_store.bin', 'w')      ; % open for element writing
  if ( el_unit == -1 )                                      % open failed
    fprintf ('\nWARNING: post-processing unit open failed \n')
    post = 0                          ; % turn off post-processing option
  end % if unit available
end % if open post-processing unit

for j = 1:n_e ; % loop over elements ====>> ====>> ====>> ====>> ====>> ====>>

  e_nodes = nodes (j, 1:n_n)                      ; % connectivity
  type    = el_type (j)            ; % current element type number

%                 Gather nodal 2D coordinates
  xy_e (1:n_n, 1) = x(e_nodes)   ; % x coord at el nodes
  xy_e (1:n_n, 2) = y(e_nodes)   ; % y coord at el nodes

%               set constant element properties
  if ( n_mats > n_t )    % each element has some different properties
    J = j ; % the element number
  else % use type data
    J = type ; % the type number
  end % if which property values
  el_prop (1:n_vals) = msh_properties (J, 1:n_vals) ;

%     Allocate and clear element arrays
  c_e    = zeros (n_i,   1) ; % clear array el sources
  c_r    = zeros (n_i,   1) ; % clear array el reactions
  p_to_F = zeros (n_i, n_n) ; % resultant of local line load
  w_e    = zeros (n_n,   1) ; % variable line load nodal values
  S_e    = zeros (n_i, n_i) ; % clear array el stiffness

  s_L = zeros (n_i, n_i) ; % local coordinate stiffness
  t_L = zeros (n_i, n_i) ; % transform local to system coordinates
  c_L = zeros (n_i, 1)   ; % local coordinate line load result

%               SET ELEMENT PROPERTIES & GEOMETRY
  A   = el_prop (1)       ; % area
  E   = el_prop (2)       ; % material modulus
  I   = el_prop (3)       ; % section inertia
  w_e = el_prop (4:n_n+3) ; % variable nodal line load

%-->   find member length and direction cosines
  dx  = x(e_nodes(n_n)) - x(e_nodes(1)) ; % x length
  dy  = y(e_nodes(n_n)) - y(e_nodes(1)) ; % y length
  L_e = sqrt (dx * dx + dy * dy)        ; % total length
  IbL = I / L_e  ; IbL2 = I / L_e^2     ; % bending constants

%      ELEMENT CONDUCTION AND INTERNAL SOURCE MATRICES 

%   Linear axial bar and cubic bending. DOF = u, v, r, u, v, r
%   Form arrays in local axes, transform.     1  2  3  4  5  6

% Local coordinate (_L) stiffness, for easy viewing  
%              (axial dof 1&4 bending 2,3 & 5,6)
  axial = [ 1,    0,       0,   -1,    0,       0 ; % bar stiffness
            0,    0,       0,    0,    0,       0 ;
            0,    0,       0,    0,    0,       0 ;
           -1,    0,       0,    1,    0,       0 ;
            0,    0,       0,    0,    0,       0 ;
            0,    0,       0,    0,    0,       0 ] * E * A / L_e ;
  bend  = [ 0,    0,       0,    0,    0,       0    ; % Beam stiff
            0,  12*IbL2,  6*IbL, 0, -12*IbL2,  6*IbL ;
            0,   6*IbL,   4*I,   0,  -6*IbL,   2*I   ;
            0,    0,       0,    0,    0,       0    ;
            0, -12*IbL2, -6*IbL, 0,  12*IbL2, -6*IbL ;
            0,   6*IbL,   2*I,   0,  -6*IbL,   4*I   ] * E / L_e ;
  s_L = axial + bend         ; % superimposed total frame stiffness

%  Map beam line load to node forces & moments;  c_L = p_2_F*w_e 
  if ( any (w_e) )         ; % then data input, form forcing vector 
    p_to_F = [ 0,       0      ;% integral of cubic beam deflection
               21,      9      ;           % times linear line load  
               3*L_e,   2*L_e  ;         % interpolations gives 6x2   
               0,       0      ;      % rectangular transfor matrix
               9,       21     ;  
              -2*L_e   -3*L_e  ] * L_e / 60 ;

    c_L = p_to_F (1:n_i, 1:n_n)*w_e(1:n_n);% forces moments @ nodes
  end ;  % if for set up line load nodal resultants

%     Define local to system DOF transformation matrix
  cx  = dx / L_e ; cy   = dy / L_e           ; % direction cosines
  t_L = [ cx   cy   0   0   0   0 ; % (inverse t_L = transpose t_L)
         -cy   cx   0   0   0   0 ;
          0    0    1   0   0   0 ;
          0    0    0   cx  cy  0 ; % joint 2
          0    0    0  -cy  cx  0 ;
          0    0    0   0   0   1 ] ;  

% Transform from local coordinates to system coordinates
  S_e = t_L' * s_L * t_L                             ; % stiffness
  c_e = t_L' * c_L                ; % resultant forces and moments

%B disp(c_L)
%B disp(c_e)

%  Save element arrays needed for post-processing stresses
  if ( post == 1 ) % save arrays for this element
    %b fwrite (el_unit, xy_e,    'double') ; % save coordinates
    fwrite (el_unit, L_e,     'double') ; % save length
    fwrite (el_unit, el_prop, 'double') ; % save properties
    fwrite (el_unit, S_e,     'double') ; % save beam stiffness
    fwrite (el_unit, c_e,     'double') ; % save line load resultant
    fwrite (el_unit, t_L,     'double') ; % save line load resultant
  end ; % if save element arrays

  if ( show_e > 0 ) % then list some element matrices
    if ( j <= show_e ) % show the first few
      fprintf ('Local Stiffness for Element %i \n', j)  ; disp(s_L)
      fprintf ('Local Loading for Element %i \n', j)    ; disp(c_L)
      fprintf ('Rotation Matrix for Element %i \n', j)  ; disp(t_L)
      fprintf ('Global Stiffness for Element %i \n', j) ; disp(S_e)
      fprintf ('Global Loading for Element %i \n', j)   ; disp(c_e)
    end ; % if elements wanted
  end ; % if show_e

%         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 sq
  c (rows) = c (rows) + c_e              ; % add to sys column

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

if ( show_s == 1 ) % then list each system matrix
  fprintf ('System matrix S = \n') ; disp(S)
  fprintf ('System matrix c = \n') ; disp(c)
end ; % if show_s
fprintf ('Total x input forces  = %g \n', sum (c(1:n_g:n_d)))
fprintf ('Total y input forces  = %g \n', sum (c(2:n_g:n_d)))
fprintf ('Total z input couples = %g \n', sum (c(3:n_g:n_d)))

if ( all ( c == 0) & all ( EBC_value == 0 )) % then no solution
  fprintf ('WARNING: no loads, was load_pt = 1 set ? \n')
  error ('No loads or non-zero displacements applied')
end ; % if no solution
%       *** END ELEMENT GENERATION

% ========================================================================
%                          PHASE 4 
%         APPLY BC AND CONSTRAINTS (Make problem unique)
% ========================================================================
fprintf ('\n')
fprintf ('Begin Phase 4, modify matrices for EBC \n')

%    ALLOCATE STORAGE FOR OPTIONAL REACTION RECOVERY
if ( EBC_count > 0 )                      ; % reactions occur
  % EBC_row = row from S for each EBC, for later reaction calculation
  % EBC_col = row from c for each EBC, for later reaction calculation
  EBC_row = zeros (EBC_count, n_d)        ; % reaction data
  EBC_col = zeros (EBC_count, 1)          ; % reaction data
  [EBC_row, EBC_col] = save_reaction_matrices (EBC_flag, S, c);
end ; % if essential BC exist (almost always true)

%         ENFORCE ESSENTIAL BOUNDARY CONDITIONS
save_resultant_load_vectors (n_g, c)
if ( EBC_count > 0 )                       ; % reactions occur
  [S, c] = enforce_essential_BC (EBC_flag, EBC_value, S, c);
  fprintf ('NOTE: applied %i essential boundary conditions \n', EBC_count)
else
  fprintf ('WARNING: no essential boundary conditions !!! \n')
end ; % if essential BC exist (always true except free-free vibration)

if ( show_bc == 1 ) % then list each system matrix after EBC
  fprintf ('After EBC system matrix S = \n') ; disp(S)
  fprintf ('After EBC system matrix c = \n') ; disp(c)
end ; % if show_s

% ========================================================================
%                          PHASE 5
%         SOLVE LINEAR EQUATION SYSTEM AND ITS REACTIONS
% ========================================================================
fprintf ('Begin Phase 5, solve modified system \n') 
Ans = S \ c                       ; % Compute nodal unknowns
list_save_results (n_g, n_m, Ans) ; % save and print

%    OPTIONAL FORCE REACTION RECOVERY & SAVE 
if ( EBC_count > 0 )                    ; % reactions exist ?
  % EBC_react = nodal source required to maintain EBC value
  EBC_react   = zeros (EBC_count, 1)    ; % pre-allocate
  [EBC_react] = recover_reactions_print_save (n_g, n_d, ...
                EBC_flag, EBC_row, EBC_col, Ans); % reaction to EBC
end ; % if EBC existed

% =====================================================================
%                             PHASE 6
%                      POST-PROCESS ELEMENTS
%  Get location and values of derivaties at all q-point, and/or
% =====================================================================
fprintf ('Begin Phase 6, use the answers \n')

if ( post == 1 ) % then post-process all element stresses
  if ( el_unit > 0 )                    ; % close written binary files
    fclose (el_unit)                             ; % file el_store.bin
    el_unit = fopen ('el_store.bin', 'r')         ; % open for reading
  end ; % if element arrays were saved

fprintf('\n'); fprintf('Individual Element Reactions \n')   ; % header
fprintf('Global Reactions: Rx_1, Ry_1, M_1, Rx_2, Ry_2, M_2 \n')
fprintf('Local Reactions:  F1_x, F1_y, M1z, F2_x, F2_y, M2z \n')

for k = 1:n_e ; % loop over elements ====>> ====>> ====>> ====>> ====>>

%       Gather answers
  type    = el_type (k)            ; % current element type number
  e_nodes = nodes (k, 1:n_n)                      ; % connectivity
  [rows]  = get_element_index (n_g, n_n, e_nodes)   ; % eq numbers
  Ans_e (1:n_i) = Ans(rows)                ;  % gather element results

%  Recover previously calculated element arrays
  %b [xy_e]    = fread (el_unit, [n_n, n_s],  'double'); % coordinates
  [L_e]     = fread (el_unit, [1,1],       'double'); % length
  [el_prop] = fread (el_unit, [n_vals, 1], 'double'); % properties
  [S_e]     = fread (el_unit, [n_i, n_i],  'double'); % beam stiffness
  [c_e]     = fread (el_unit, [n_i, 1],    'double'); % line load resultant
  [t_L]     = fread (el_unit, [n_i, n_i],  'double'); % rotation matrix

%   Optional output of nodal reactions on every element
  if ( el_react ) % list nodal heat flows for this element

%   Assign any point source to the first element with that node
    if ( load_pt ) % then point sources to assign (once)
      for eq = 1:n_i                      % loop over element nodes
        k = rows (eq)                   ; % find system DOF number
        if ( c_pt (k) ~= 0 )              % found an input point value
          c_e (eq) = c_e (eq) + c_pt (k); % copy to this node
          c_pt (k) = 0                  ; % prevent a second use
        end ; % if first appearance of this point source in the system
      end ; % for eq node of the current element
    end ; % if load_pt has been assigned

%   Finally, get the element reactions and list them in global coordinates
    c_m = S_e * Ans_e' ; % global member forces due to displacements

    fprintf('\n')
    fprintf('Element number %i, length %g \n', k, L_e)
    fprintf('External system reaction components: \n')
    fprintf('Rx_1        Ry_1       RM_1       Rx_2       Ry_2       RM_2 \n')
    fprintf('%10.4e %10.4e %10.4e %10.4e %10.4e %10.4e \n', ...
             c_m(1), c_m(2), c_m(3), c_m(4), c_m(5), c_m(6) )
%   Transform back to element local coordinates and list
    c_L = t_L * c_m    ; % overwrite c_L storage w local global reactions

    fprintf('\n') ; fprintf('External local reaction components: \n')
    fprintf('F1_x        F1_y       M1_z       F2_x        F2_y       M2_z \n')
    fprintf('%10.4e %10.4e %10.4e %10.4e %10.4e %10.4e \n', ...
             c_L(1), c_L(2), c_L(3), c_L(4), c_L(5), c_L(6) )
  end ; % if el_react list
end ; % for each k element in mesh   <<==== <<==== <<==== <<==== <<====

% =====================================================================
%                             PHASE 7
%           ESTIMATE ERROR IN EACH ELEMENT AND THE SYSTEM
% =====================================================================
%      None implemented yet for this application
end ; % of Frame_2D_cubic              

%=================== Running gives =====================================
% Frame_2D_cubic
% 
% (Echo of msh_remarks.txt) 
% ============== Begin Application Remarks ====================
% G.R. Buchanan Example 4.14, A plane frame structure
% 
%  wwwwwwwwwwwwww       w = 1000 lb/ft = 83.333 lb/in
% 1-----(1)-----2       L_1 = L_2 = 10 ft = 120 inches
% pin             \      e_2 30 deg from vertical
%                  \      A = 10 in^2
%                  (2)     E = 30e6 lb/in^2
%                    \      I = 100 in^4
%                     \
%                     3, pin
% Nodal dof: 1=U axial disp, 2=V transverse disp 3=Z-angle
% Element type = 1 (a parametric line)
% Element connection: two nodes per element
% Properties are: A, E, I, 2 nodal line loads, Rho
% Exact result at node 2 u_x = -1.5743e-03 u_y = -4.0600e-03  
%                        rotation = 9.9722e-04 radians
% =============== End Application Remarks =======================
% 
% User control variables 
% n_g number of DOF per node        = 3 
% n_n number of nodes per element   = 2 
% n_p dimension of parametric space = 1 
% n_q number of quadrature points   = 0 
% n_r number of rows in B_e         = 1 
% n_s dimension of physical space   = 2 
% 
% User logic flags: turn on(1)/off(0) 
% echo_p   listing of properties         = 1 
% el_react listing of element reactions  = 1 
% load_pt  point force and/or moment     = 0 
% post     post-processing option        = 0 
% show_bc  debug list S & c, after EBC   = 0 
% show_e   list some S_e, and c_e        = 0 
% show_s   list assembled S and c        = 0 
% stiff_pt point stiffness (can be zero) = 0 
% 
% Begin Phase 2, input mesh, properties, sources and EBC 
% Read 3 nodes. 
% (Echo node, file msh_bc_xyz.txt) 
% node, bc_flags, 2 coordinates 
% 1,    110 0   103.923 
% 2,    000 120 103.923 
% 3,    110 180 0 
% 
% (Echo of file msh_typ_nodes.txt) 
% Read 2 elements with (ignored) type & 2 nodes each. 
% 1     1     2
% 1     2     3
% 
% Read 2 materials with 5 poperties 
% 5 material values per element 
% (Echoing text file msh_properties.txt) 
% 1.0000e+01   3.0000e+07   1.0000e+02  -8.3333e+01  -8.3333e+01
% 1.0000e+01   3.0000e+07   1.0000e+02            0            0
% 
% Above properties are: A, E, I, 2 nodal line loads 
% For a planar frame 
% Note: expecting 4 essential BC values. 
% 
% Applied Displacement Boundary Conditions: 4 
% (Echo of file load msh_ebc.txt) 
% Node, DOF (1=u, 2=v, 3=r), Value. 
% 1     1     0
% 1     2     0
% 3     1     0
% 3     2     0
% 
% Begin Phase 3, assemble (scatter) element arrays 
% Total x input forces  = 0 
% Total y input forces  = -10000 
% Total z input couples = 0 
% 
% Begin Phase 4, modify matrices for EBC 
% 
% Resultants of all input sources: 
% Node, DOF (1=axial force, 2=shear force, 3=couple), Value 
% 1     2 -5000 
% 1     3 -100000 
% 2     2 -5000 
% 2     3  100000 
% Totals = 1.0e+03 *[ 0  -10.0000   0]
% 
% NOTE: applied 4 essential boundary conditions 
% Begin Phase 5, solve modified system 
% 
% Computed Solution: 
% Node, 3 results per node 
% 1    0.0000e+00  0.0000e+00 -1.5494e-03 
% 2   -1.5743e-03 -4.0600e-03  9.9722e-04 
% 3    0.0000e+00  0.0000e+00 -4.5619e-04 
% 
% Recovered Reactions at Displacement or Slope BC: 4 
% Node, DOF (1=force, 2=couple), Value 
% 1     1  3935.66 
% 1     2  4394.41 
% 3     1 -3935.66 
% 3     2  5605.59 
% Total force and couple = 1.0e+03 *[0.0000   10.0000   0]
% 
% Begin Phase 6, use the answers 
% 
% Individual Element Reactions 
% Global Reactions: Rx_1, Ry_1, M_1, Rx_2, Ry_2, M_2 
% Local Reactions:  F1_x, F1_y, M1z, F2_x, F2_y, M2z 
% 
% Element number 1, length 120 
% External system reaction components: 
% Rx_1        Ry_1        RM_1        Rx_2       Ry_2       RM_2 
% 3.9357e+03 -6.0559e+02 -1.0000e+05 -3.9357e+03 6.0559e+02 2.7329e+04 
% 
% External local reaction components: 
% F1_x        F1_y        M1_z        F2_x       F2_y       M2_z 
% 3.9357e+03 -6.0559e+02 -1.0000e+05 -3.9357e+03 6.0559e+02 2.7329e+04 
% 
% Element number 2, length 120 
% External system reaction components: 
% Rx_1        Ry_1       RM_1        Rx_2       Ry_2       RM_2 
% 3.9357e+03 -5.6056e+03 7.2671e+04 -3.9357e+03 5.6056e+03 1.2203e-11 
% 
% External local reaction components: 
% F1_x       F1_y       M1_z        F2_x        F2_y       M2_z 
% 6.8224e+03 6.0559e+02 7.2671e+04 -6.8224e+03 -6.0559e+02 1.2203e-11 
% 
% G.R. Buchanan Example 4.14, A plane frame structure
%                   SYSTEM REACTIONS
%  w = -1000 lb/ft = -83.333 lb/in gives -10,000 lb
% 3,935.6 wwwwwwwwwwwwww       
%   ----> 1-----(1)-----2     
%         ^              \      
%         | 4,394         \    
%                         (2)           
% E=30e6 lb/in^2            \      
% I = 100 in^4               \  3,935.6
%                             3 <----
%                     5,605.6 ^
%                             |
%
% ELEMENT 1 GLOBAL = LOCAL REACTIONS
%
%          (-10,000 lb)
%  3,935.6 wwwwwwwwwwwwww    3,935.6
%    ----> 1-----(1)-----2 <----    
%          ^             ^ 
%  4,394.4 |         ^   | 5,605.6 lb
%           -72,670.8 \___/
%               ft-lb
%
% ELEMENT 2 GLOBAL REACTIONS
%
%                    /---\
%                   V     \ +72,670.8 ft-lb
%          3,935.6 ---->2     
%                       |\      
%               5,605.6 V \    
%                         (2)           
%                           \      
%                            \  3,935.6
%                             3 <----
%                     5,605.6 ^
%                             |
% ELEMENT 2 LOCAL REACTIONS
%
%           +605.6 ^ lb
%                  |  (120 in.)
%  6,822 lb ---> | 2-----(2)-----3 <--- -6,822
%                \               |
%                 --> +72,670.8  V -605.6 lb
%                      ft-lb

%  +++++++++++++ functions in alphabetical order +++++++++++++++++
%      Newer versions may be in directory Akin_FEA_Lib

function echo_user_remarks  % ==========================
file_id = fopen('msh_remarks.txt', 'r'); ;% open for read only
if ( file_id > 0 )                 % then file exists
  fprintf('\n(Echo of msh_remarks.txt) \n')
  while feof(file_id) ==0          % read until an end of file
    remark_line = fgets(file_id) ; % read a string and append \n
    fprintf ('%s',remark_line)        ; % echo the line
  end ; % while data exists
end ; % if file supplied
fprintf('\n')
% end function echo_user_remarks  % ==========================

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
Diag = max ( max (S) )     ; % for better condition number
for k = 1:n_d                % check all DOF for essential BC 
  if ( flag_EBC (k) )                        % then EBC here
%   Carry known columns*EBC to RHS. Zero that column and row.
    EBC = value_EBC (k)             ; % recover EBC value
    C (:) = C (:) - EBC * S (:, k) ; % carry known column to RHS
    S (:, k) = 0 ; S (k, :) = 0    ; % clear, restore symmetry
%   Insert EBC identity, Diag * EBC_dof = Diag * EBC_value
    S (k, k) = Diag                ; % insert identity into S row
    C (k)    = Diag * EBC          ; % insert identity into C row
  end ; % if EBC for this DOF
end ; % for over all k-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.txt            ; % 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.txt')
  end ; % if user error
  fprintf ('\nRead %g point sources. \n', n_u)
  fprintf ('(Echo of file msh_load_pt.txt) \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.txt            ; % 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.txt) \n')               
  fprintf ('Node, DOF (1=u, 2=v, 3=r), 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.txt \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.txt        ; % 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.txt') 
end ; % if error 
n_n = size (msh_typ_nodes,2) - 1 ; % nodes per element       
fprintf ('(Echo of file msh_typ_nodes.txt) \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, PD, x, y, z] = get_mesh_nodes (n_g); % =======
% input file controls (for various data generators)

%    READ MESH AND EBC_FLAG SEQUENTIAL INPUT DATA
load msh_bc_xyz.txt       ; % 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 ('\n Error missing nodal file msh_bc_xyz.txt')
end ; % if error

n_s = size (msh_bc_xyz,2) - 1 ; % number of space dimensions
fprintf ('Read %i nodes. \n', n_m)
fprintf ('(Echo node, file msh_bc_xyz.txt) \n')
fprintf ('node, bc_flags, %i coordinates \n', n_s)

msh_bc_xyz (:, 1)= round (msh_bc_xyz (:, 1)); % require integer
PD = 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

% Echo data with pretty print from control values
for j = 1:n_m  % loop over all nodes ===> ===> ===> ===> ===> ===>
  if ( n_s == 1 )
    if ( n_g == 2 )
      fprintf ('%i, %2.2i %g \n', j, PD(j), x(j))
    else
      fprintf ('%i, %i    %g \n', j, PD(j), x(j))
    end ; % if dof per node
  elseif (n_s == 2)
    if ( n_g == 2 )
      fprintf ('%i, %2.2i %g %g \n', j, PD(j), x(j), y(j))
    elseif ( n_g == 3 )
      fprintf ('%i, %3.3i %g %g \n', j, PD(j), x(j), y(j))
    else
      fprintf ('%i, %i    %g %g \n', j, PD(j), x(j), y(j))
    end ; % if dof per node
  else % n_s = 3
    if ( n_g == 2 )
      fprintf ('%i, %2.2i %g %g %g \n', j, PD(j), x(j), y(j), z(j))
    else
      fprintf ('%i, %i    %g %g %g \n', j, PD(j), x(j), y(j), z(j))
    end ; % if dof per node
  end ; % if space dimension
end ; % loop for j-th node <=== <=== <=== <=== <=== <=== <=== <===
% end get_mesh_nodes    % ======================================

function [n_mats, n_vals, msh_properties] = get_mesh_properties (echo_p, n_e, n_t)
load msh_properties.txt                       ; % one row per element
% msh_properties (n_mats, n_vals) = system properties array

n_mats = size(msh_properties, 1) ; % == 1 if homogeneous
n_vals = size(msh_properties, 2) ; % number of application properties
% n_mats = number of materials, = n_e, or n_t if homogeneous
% n_vals = number of property values for each element
fprintf ('Read %i materials with %i poperties \n', n_mats, n_vals)

% el_prop   = property values of current element material (n_vals, 1)

if ( n_mats == n_t )        % then element type data are honogeneous
  fprintf ('%i homogeneous material values per element type \n', n_vals)
else
  fprintf ('%i material values per element \n', n_vals)
end ; % if
if ( n_mats > n_e )
  error ('ERROR: number of material sets exceeds number of elements')
end ; % if
if ( echo_p == 0 )
  fprintf ('\nNOTE individual element properties echo skipped \n')
else
  fprintf ('(Echoing text file msh_properties.txt) \n')
  format short e ;  disp(msh_properties) ; format short
end ; % if
% end get_mesh_properties % ============================================

function list_save_results (n_g, n_m, Ans) % ===========================
fprintf ('\nComputed Solution: \n') ;
fprintf('Node, %i results per node \n', n_g)
Ans_r = reshape (Ans, n_g, n_m)'            ; % pretty shape

% save results (displacements) to MODEL file: node_results.txt
fid = fopen('node_results.txt', 'w')          ; % open for writing
for j = 1:n_m                                 ; % save displacements
 if ( n_g == 1 )
  fprintf (fid, '%g \n', Ans_r (j, 1:n_g)) ;
  fprintf ('%i   %11.4e \n', j, Ans_r (j, 1:n_g)) ;
 elseif ( n_g == 2 )
  fprintf (fid, '%g %g \n', Ans_r (j, 1:n_g)) ;
  fprintf ('%i   %11.4e %11.4e \n', j, Ans_r (j, 1:n_g)) ;
 elseif ( n_g == 3 )
  fprintf (fid, '%g %g %g \n', Ans_r (j, 1:n_g)) ;
  fprintf ('%i   %11.4e %11.4e %11.4e \n', j, Ans_r (j, 1:n_g));
 elseif ( n_g == 4 )
  fprintf (fid, '%g %g %g %g \n', Ans_r (j, 1:n_g)) ;
  fprintf ('%i   %11.4e %11.4e %11.4e %11.4e \n', j, Ans_r (j, 1:n_g));
 elseif ( n_g == 5 )
  fprintf (fid, '%g %g %g %g %g \n', Ans_r (j, 1:n_g)) ;
  fprintf ('%i   %11.4e %11.4e %11.4e %11.4e %11.4e \n', j, Ans_r (j, 1:n_g));
 elseif ( n_g == 6 )
  fprintf (fid, '%g %g %g %g %g %g \n', Ans_r (j, 1:n_g)) ;
  fprintf ('%i %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e \n', j, Ans_r (j, 1:n_g));
 else
   error ('reformat list_save_results for n_g > 6.')
 end ; % if
end ; % for j DOF
% end list_save_results    % =========================================

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.txt
  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.txt', '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 ('%i %i %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.txt
  n_d = size (C, 1)                     ; % number of system DOF
  fprintf ('\nResultants of all input sources: \n')
  fprintf ('Node, DOF (1=axial force, 2=shear force, 3=couple), Value \n')
  fid = fopen('node_resultant.txt', '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