function Static_AA                % Revised 2/21/2012  
% 1D ODE solver using linear (L2), or quadratic (L3)
% or cubic (L4) elements with constant properties and
% constant Jacobian via analytic matrix forms

% -d/dx (E_e A_e du/dx) + d_e u - A_e Q_e  = 0
% u   = unknown
% E_e = el_prop (1) ; % modulus
% A_e = el_prop (2) ; % area
% Q_e = el_prop (3) ; % source per unit volume
% d_e = el_prop (4) ; % surface product or mass density/length
%           NOTE: The general form is 
% -d/dx (E_e A_e du/dx) + d_e (u-u_ref) - A_e Q_e  = 0
% BUT, u_ref = 0 usually, else enter a fake Q_e as
% Q_e = (d_e u_ref + A_e Q_e) / A_e as the input
%............................................................. 
%     Mech 417, Finite Element Analysis, Rice University 
%     Copyright J.E. Akin 2012. All rights reserved.
%............................................................. 
% NOTE: input text files begin with msh_ and end with .tmp
% see file FEA_notation.m for variable descriptions
global natural ; % flag for use of Natural or Unit coordinates
 
% echo information about input and output text files
% echo_files_info   % uncomment column 1, if desired

% echo user information about this specific application
echo_user_remarks

% =====================================================================
%                             PHASE 1
%           *** SET OR READ USER PROBLEM DEPENDENT CONTROLS *** 
% =====================================================================
%   set default user numerical values (or add read function)
n_g = 1 ; % number of degrees of freedom per node 
N_N = 2 ; % number of EXPECTED nodes per element
n_p = 1 ; % dimension of parametric space
n_q = 2 ; % number of quadrature points required
%           Exact for polynomial of degree (2n_q - 1)
n_r = 1 ; % number of rows in element gradient matrix, B_e
N_S = 1 ; % 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 = 0  ; % turn on(1)/off(0) listing of element reactions
integral = 0  ; % turn on(1)/off(0) listing of solution integral
load_pt  = 1  ; % turn on(1)/off(0) point source (can be zero)
mass_pt  = 0  ; % turn on(1)/off(0) point mass (can be zero)
natural  = 0  ; % turn on(1)/off(0) the Natural coordinates
post     = 0  ; % turn on(1)/off(0) post-processing option
show_e   = 0  ; % turn on(>0)/off(0) list some S_e, M_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 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 ('integral listing of solution integral  = %i \n', integral)
fprintf ('load_pt  point source (can be zero)    = %i \n', load_pt )
fprintf ('mass_pt  point mass (can be zero)      = %i \n', mass_pt )
fprintf ('natural  use Natural coordinates       = %i \n', natural )
fprintf ('post     post-processing option        = %i \n', post    )
fprintf ('show_e   list some S_e, M_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 | integral > 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 Begin Phase 2, input mesh, properties, sources and EBC \n')

%    Read mesh nodal input data files 
%    (from msh_bc_xyz.tmp, 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 fake (un-used) in this application, so enter 0
% 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 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.tmp, 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

%    Read any material properties for the elements
%    (from msh_properties.tmp, 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: E_e, A_e, Q_e, d_e \n')
fprintf ('For -d/dx (E_e A_e du/dx) + d_e u - A_e Q_e  = 0 \n')

%   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
% M         = mass matric (if any) 
% S         = stiffness matrix       
% Type_integral = integral of solution (u) in each element 
% Type_i_d_e_u  = integral of scaled solution (d_e*u) in each element 

Ans       = zeros (n_d, 1)     ; % results
Body      = zeros (n_s, 1)     ; % initialize volumetric source
C         = zeros (n_d, 1)     ; % total 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
M         = zeros (n_d, n_d)   ; % generalized mass
S         = zeros (n_d, n_d)   ; % generalized stiffness
Type_integral = zeros (n_t, 1) ; % integral of the solution
Type_i_d_e_u  = zeros (n_t, 1) ; % integral of d_e times the solution

%   Read point source data, if any, store and add to C
%   (from file msh_load_pt.tmp, 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 mass data, if any, and add to M diagonal
%   (from file msh_mass_pt.tmp, one line per point mass)
if ( mass_pt > 0 )                      ; % need point masses data
  [M] = get_and_add_point_masses (n_g, n_m, M); % get point masses
end % if any point mass expected

%   Read point stiffness data, if any, and add to S diagonal
%   (from file msh_stiff_pt.tmp, 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.tmp 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

%   Prepare a default general 2D mesh plot (uncomment as desired)
%   using /clear/www/htdocs/mech517/Matlab_Plots/mesh_shrink_plot.m
% mesh_shrink_plot % show a space between elements
% pause (3) % wait 3 seconds before next plot
%   or using /clear/www/htdocs/mech517/Matlab_Plots/mesh_plot.m
% mesh_plot    % no space between elements, mesh_plot (inc_e, inc_p)
%   or using true element shape with
%      /clear/www/htdocs/mech517/Matlab_Plots/L2_mesh_plot.m
%   or /clear/www/htdocs/mech517/Matlab_Plots/L3_mesh_plot.m
%   or /clear/www/htdocs/mech517/Matlab_Plots/L4_mesh_plot.m

%  Set an integer to simplify element type selection below
domain = n_p * 100 + n_g * 10 + n_n ; % parametric reference code

%    *** 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

%                *** BEGIN ELEMENT LOOPS ***
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 coordinates (using vector subscripts)
  xy_e (1:n_n, 1) = x(e_nodes)            ; % x coord at el nodes 
  L_e = abs (xy_e (n_n, 1) - xy_e (1, 1)) ; % length

%   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) ;

%   ELEMENT, SQUARE AND COLUMN ARRAY GENERATION

%     shorthand property names
  if ( n_vals < 4 ) % then sing property data
    error ('This application requires 4 properties')
  end % if n_vals

% -d/dx (E_e A_e du/dx) + d_e u - A_e Q_e  = 0
% u   = unknown
% E_e = el_prop (1) ; % modulus
% A_e = el_prop (2) ; % area
% Q_e = el_prop (3) ; % source per unit volume
% d_e = el_prop (4) ; % surface product or mass density/length

  E_e = el_prop (1) ; % modulus
  A_e = el_prop (2) ; % area
  Q_e = el_prop (3) ; % source (weight) per unit volume
  d_e = el_prop (4) ; % surface product or mass density/length

  K_e = E_e * A_e / L_e ; % classic axial stiffness
  m_e = L_e * d_e       ; % total convection or mass
  W_e = A_e * L_e * Q_e ; % total source (weight)

%    Analytical forms of S_e, C_e, M_g and M_e
% Select type of element
  switch n_p   % dimension of parametric space
    case 1 % line elements
      switch n_n % number of nodes per element
        case 2 % linear element L2
          S_e = K_e * [ 1,   -1 ; ...
                       -1,    1 ]     ; % stiffness matrix
          C_e = W_e / 2 * [ 1,  1 ]'  ; % load vector
          M_g = 1 / 6 * [2,  1 ; ...
                         1,  2 ]   ; % generalized mass matrix
          M_e = m_e * M_g          ; % surroundings matrix
        case 3 % quadratic element L3
          S_e = K_e / 3 * [ 7, -8,  1 ; ...
                           -8, 16, -8 ; ...
                            1, -8,  7 ] ; % stiffness matrix
          C_e = W_e / 6 * [ 1, 4, 1 ]'  ; % load vector
          M_g = 1 / 30 * [4,  2, -1 ; ...
                          2, 16,  2 ; ...
                         -1,  2   4 ] ; % generalized mass matrix
          M_e = m_e * M_g             ; % surroundings matrix
        case 4 % cubic element L4
          S_e = K_e / 40 * [ 148, -189,   54,  -13 ; ...
                            -189,  432, -297,   54 ; ...
                              54, -297,  432, -189 ; ...
                             -13,   54, -189,  148 ] ; % stiffness
          C_e = W_e / 8 * [1, 3, 3, 1]'              ; % load vector
          M_g = 1 / 1680 * [128,  99, -36,  19 ; ...
                             99, 648, -81  -36 ; ...
                            -36, -81, 648,  99 ; ...
                             19, -36,  99, 128]   ; % generalized mass 
          M_e = m_e * M_g                         ; % surroundings matrix
        otherwise
          error ('Expected n_n = 2 or 3 or 4 \n')
      end % switch on number of nodes per element
    otherwise
      error ('Add library for n_p > 1 \n')
    end % switch on dimension of parametric space

  if ( show_e > 0 ) % then list some element matrices
    if ( j <= show_e ) % show the first few
      fprintf ('Element matrix S_e = \n') ; disp(S_e)
      fprintf ('Element matrix M_e = \n') ; disp(M_e)
      fprintf ('Element matrix C_e = \n') ; disp(C_e)
    end % if elements wanted
  end % if show_e

  if ( m_e ~= 0 ) % then surrounding convection or elastic foundation
    S_e = S_e + M_e  ; % effective element square matrix
    if ( show_e > 0 ) % then list some element matrices
      if ( j <= show_e ) % show the first few
        fprintf ('Element matrix S_e = \n') ; disp(S_e)
      end % if elements wanted
    end % if show_e
  end % if m_e 

%         SCATTER TO (ASSEMBLE INTO) SYSTEM ARRAYS 

  % rows = vector subscript converting element to system eq numbers
  [rows] = get_element_index (n_g, n_n, e_nodes); % eq numbers

% Insert completed element matrices into system matrices
  S (rows, rows) = S (rows, rows) + S_e  ; % add to system stiffness
  C (rows)       = C (rows)       + C_e  ; % add to sys force/couples

%       *** END SCATTER TO SYSTEM ARRAYS ***
 
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 input source terms = %g \n', sum (C))
%       *** END ELEMENT GENERATION

% ========================================================================
%                          PHASE 4
%         APPLY BC AND CONSTRAINTS (Make problem unique)
% ========================================================================
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 storage
  EBC_col = zeros (EBC_count, 1)          ; % reaction data storage
  [EBC_row, EBC_col] = save_reaction_matrices (EBC_flag, S, C);
end % if essential BC exist (almost always true)

% in the future, enforce all MPC equations here

% ENFORCE ESSENTIAL BOUNDARY CONDITIONS (and insert identities)
if ( EBC_count > 0 )                       ; % reactions occur
  [S, C] = enforce_essential_BC (EBC_flag, EBC_value, S, C);
% Note: Both S and C have been changed, original point sources
%       are now still in C_pt, but are deleted in el_reactions
end % if essential BC exist (almost always true)

% ========================================================================
%                          PHASE 5
%         SOLVE LINEAR EQUATION SYSTEM AND ITS REACTIONS, OR
%                  SOLVE EIGENVALUE SYSTEM
% ========================================================================
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 
%  Get integral of the solution and/or
%  Sum all the element or nodal values
% =====================================================================
if ( post == 1 ) % then post-process all element stresses
  fprintf ('Begin Phase 6, use the answers \n')

%  In future, open file to save stresses for plotting
  Surround = 0 ; % initialize any surroundings effects

  fprintf ('\nElement Post-processing: \n')
  % strain = gradient components
  % stress = flux components
  strain = zeros (n_r, 1) ; % pre-allocate strains
  stress = zeros (n_r, 1) ; % pre-allocate stresses
  H_i    = zeros (n_n, 1) ; % pre-allocate interpolation integral

%  Re-compute or read element arrays from storage
  format short e
  for j = 1:n_e ; % loop over elements ====>> ====>> ====>> ====>> ====>>

%       Gather answers
    e_nodes = nodes (j, 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
  
%                 Gather nodal coordinates
    xy_e (1:n_n, 1) = x(e_nodes)             ; % x coord at el nodes 
    L_e = abs (xy_e (n_n, 1) - xy_e (1, 1)) ; % length
  
% Compute secondary values for the current element
%   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) ;

%            Analytical forms of matrices
    E_e = el_prop (1) ; % modulus
    A_e = el_prop (2) ; % area
    Q_e = el_prop (3) ; % source (weight) per unit volume
    d_e = el_prop (4)     ; % mass density
    K_e = E_e * A_e / L_e ; % classic axial stiffness
    W_e = A_e * L_e * Q_e ; % total source (weight)
    m_e = A_e * L_e * d_e ; % total mass

%  list (and save for plots?) results
    fprintf ('Element number %i \n', j)    
    fprintf ('Length ........................... %g \n', L_e)
    fprintf ('Body source ...................... %g \n', W_e)
    fprintf ('Mass ............................. %g \n', m_e )
    fprintf ('Stiffness ........................ %g \n', K_e)

%      Recover element gradient at centroid or interior nodes
%      Recover element integral and/or reactions

% Select type of element
    switch n_n % number of nodes per element
      case 2 % linear element
        B_e = [-1,   1] / L_e        ;  % center gradient
        H_i = [1,    1] * L_e / 2    ;  % interpolation integral
      case 3 % quadratic element
        B_L = [-3   4  -1] / L_e     ;  % dH/dx at r=0
        B_e = [-1   0   1] / L_e     ;  % dH/dx at r=0.5
        B_R = [1   -4   3] / L_e     ;  % dH/dx at r=1
        H_i = [1    4   1] * L_e / 6 ;  % interpolation integral
      case 4 % cubic element
        B_L = [-11  18  -9  2] / L_2 / 2 ;  % dH/dx at r=0
        B_e = [1   -27  27 -1] / L_2 / 8 ;  % dH/dx at r=0.5
        B_R = [-2    9 -18 11] / L_2 / 2 ;  % dH/dx at r=1
        H_i = [1  3  3  1] * L_e / 8     ;  % interpolation integral
      otherwise
        error ('Add gradient for n_n > 4 \n')
    end % switch

%   compute the center point generalized strains & stresses
    strain = B_e * Ans_e' ;  % element strains (gradient)
    stress = E_e * strain ;  % element stresses (fluxes)

%  list (and save for plots?) results
    fprintf ('Center gradient .................. %g \n', strain)
    fprintf ('Center flux ...................... %g \n', stress)

    if ( n_n == 3 | n_n == 4 ) % quadratic or cubic element
      stress = E_e * B_L * Ans_e' ; % left end stress
      fprintf ('Left flux ........................ %g \n', stress)
      stress = E_e * B_R * Ans_e' ; % right end stress
      fprintf ('Right flux ....................... %g \n', stress)
    end % if n_n

%   Now do the integral of the solution
    if ( integral == 1 & n_g == 1 ) % report solution integral
%     I_e = integral of solution over the current element
      I_e = dot (Ans_e, H_i);

      Type_integral (type) = Type_integral (type) + I_e       ;
      Type_i_d_e_u  (type) = Type_i_d_e_u  (type) + d_e * I_e ;
      fprintf ('Integral of solution, u .......... %g \n', I_e)
      fprintf ('Integral of d_e times solution ... %g \n', d_e*I_e)
    end % if scalar solution

%   Now get each element reaction
    if ( el_react == 1 ) % then list element reactions
      fprintf ('Connectivity ..................... ')
      disp (rows)
      fprintf ('Element nodal values ............. ')
      disp (Ans_e)

%     Select type of element
      switch n_n % number of nodes per element
        case 2 % linear element
          S_e = K_e * [ 1,   -1; ...
                       -1,    1]      ; % stiffness matrix
          M_g = 1 / 6 * [2,  1 ; ...
                         1,  2 ]   ; % generalized mass matrix
          M_e = m_e * M_g          ; % surroundings matrix
          C_e = W_e / 2 * [ 1,   1]'  ; % load vector
        case 3 % quadratic element
          S_e = K_e / 3 * [ 7, -8,  1  ; ...
                           -8, 16, -8  ; ...
                            1, -8,  7 ] ; % stiffness matrix
          M_g = 1 / 30 * [4,  2, -1 ; ...
                          2, 16,  2 ; ...
                         -1,  2   4 ] ; % generalized mass matrix
          M_e = m_e * M_g             ; % surroundings matrix
          C_e = W_e / 6 * [1, 4, 1]'    ; % load vector
        case 4 % cubic element L4
          S_e = K_e / 40 * [ 148, -189,   54,  -13 ; ...
                            -189,  432, -297,   54 ; ...
                              54, -297,  432, -189 ; ...
                             -13,   54, -189,  148 ] ; % stiffness
          M_g = 1 / 1680 * [128,  99, -36,  19 ; ...
                             99, 648, -81  -36 ; ...
                            -36, -81, 648,  99 ; ...
                             19, -36,  99, 128]   ; % generalized mass 
          M_e = m_e * M_g                         ; % surroundings matrix
          C_e = W_e / 8 * [1, 3, 3, 1]'              ; % load vector
        otherwise
          error ('Expected n_n = 2 or 3 or 4 \n')
      end % switch

%     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
      fprintf ('Input resultant nodal sources .... ') ; disp (C_e')

      C_r = S_e * Ans_e' - C_e                ; % element reactions

%     Now include any surrounding (d_e) sources
      if ( d_e ~= 0 )                 % then some act on this element
        C_f = -M_e * Ans_e' ; % surrounding acts on element
        fprintf ('Resultant surrounding loading .... ') ; disp (C_f')
        Surround = Surround + sum (C_f)         ; % update the effects
        C_r = C_r - C_f     ; % additional elastic surrounding effects
      end % if d_e

      fprintf ('Element nodal reactions ......... ')
      disp (C_r')
    end % if el_react

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

  if ( integral == 1 & n_g == 1) % report solution integral
    if ( n_t > 1 ) % report on each type
      for j = 1:n_t
        fprintf ('Solution integral in type %i elements        = %g \n', ...
                  j, Type_integral (j))
        fprintf ('Scaled solution integral in type %i elements = %g \n', ...
                  j, Type_i_d_e_u  (j))
        fprintf ('Solution volume in type %i elements          = %g \n', ...
                  j, Type_volume (j))
      end % for each type
    end % if multiple element types
    fprintf ('Total solution integral    = %g \n', sum (Type_integral))
    fprintf ('Total scaled   integral    = %g \n', sum (Type_i_d_e_u ))
    fprintf ('Total surroundings effects = %g \n', Surround)
  end % if scalar solution integral to report

% =====================================================================
%                             PHASE 7
%           ESTIMATE ERROR IN EACH ELEMENT AND THE SYSTEM
% =====================================================================
%      None for this application

end % function Static_AA             % =======================
   
%         ++++++++++++++++++++++++++++++++++++++++++
%         ++++++++++++++++++++++++++++++++++++++++++
%  ++++++++ active functions in alphabetical order ++++++++++++
%         ++++++++++++++++++++++++++++++++++++++++++
%         ++++++++++++++++++++++++++++++++++++++++++
 
% See /clear/www/htdocs/mech517/Matlab_Plots for graphic options 
% and http://www.clear.rice.edu/mech517/help_plot.html for help files
 
function echo_files_info  % ==========================
file_id = fopen('msh_files_info.tmp', 'r'); % open for read only
if ( file_id > 0 )                 % then file exists
  fprintf('\n(Echo of msh_files_info.tmp) \n')
  while feof(file_id) ==0          % read until an end of file
    remark_line = fgets(file_id) ; % read a string and append \n
    fprintf (remark_line)        ; % echo the line
  end % while data exists
end % if file supplied
fprintf('\n')
% end function echo_file_info  % ==========================

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

function [H_q, DLH_q] = el_shape_n_local_deriv (domain, r, s, t)
% parametric interpolation functions and their derivatives
% at local point r, or at r, s, or at r, s, t, and 
% domain = n_p * 100 + n_g * 10 + n_n ; % an integer code
% Currently using unit coordinates, 0 <= r,s,t <= 1 

 switch domain     % select from available library of elements
% One-dimensional elements
  case {112}       % two node line linear element
%   (See qp_rule_unit_Gauss for corresponding tabulated data)
    H_q = [(1 - r),  r] ; DLH_q = [-1,    1] ;

  case {113}       % three node line quadratic element
%   (See qp_rule_unit_Gauss for corresponding tabulated data)
    H_q   = [(1 - 3*r + 2*r^2), (4*r - 4*r^2), (2*r^2 - r)] ;
    DLH_q = [(-3 + 4 * r),      (4 - 8*r),     (4*r -1)] ;

  case {114}       % four node line cubic element
%   (See qp_rule_unit_Gauss for corresponding tabulated data)
    H_q   = [(1 - r*11/2 + 9*r^2 - 9*r^3/2), (9*r - 45*r^2/2 + 27*r^3/2), ...
             (-9*r/2 + 18*r^2 - 27*r^3/2),   (r - 9*r^2/2 + 9*r^3/2) ] ;
    DLH_q = [(-11/2 + 18*r - 27*r^2/2),      (9 - 45*r + 81*r^2/2), ...
             (-9/2 + 36*r - 81*r^2/2),       (1 - 9*r + 27*r^2/2) ] ;
  otherwise
    fprintf ('Current mesh domain flag = %i \n', domain)
    error ('\n Missing parametric functions for that domain')
 end % switch
% end function el_shape_n_local_deriv % =====================================

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)
% EBC_flag  = n_g individual binary flags at each node
% EBC_count = number of essential (Dirichlet) BC's
% S = stiffness matrix,       C = total soutce vector
    
 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) == 1 )                   % 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    % =================================

function [M] = get_and_add_point_masses (n_g, n_m, M) %==========
% n_g  = number of DOF per node
% n_m  = number of nodal points in mesh
% M    = system mass matrix to receive diagonal mass terms
  load msh_mass_pt.tmp            ; % node, DOF, value (eq. number)
  n_u = size(msh_mass_pt, 1)      ; % number of point masses
  if ( n_u < 1 )                  ; % missing data
    error ('No mass_pt data in text file msh_mass_pt.tmp')
  end % if user error
  fprintf ('\nRead %i point masses. (Echo of msh_mass_pt.tmp) \n', n_u)
  fprintf ('Node, DOF,  Mass_value \n')
  for j = 1:n_u                         ; % non-zero mass pts
    node  = msh_mass_pt (j, 1)          ; % global node number
    if ( node < 1 | node > n_m ) % then impossible node number
      fprint ('Point mass has invalid node number = %i \n', node)
      error  ('ERROR, Point mass has invalid node number')
    end % if
    DOF   = msh_mass_pt (j, 2)          ; % local DOF number
    if ( DOF < 1 | DOF > n_g ) % then impossible DOF number
      fprint ('Point mass has invalid DOF number = %i \n', DOF)
      error  ('ERROR, Point mass has invalid DOF number')
    end % if
    value = msh_mass_pt (j, 3)          ; % point mass value
    fprintf ('%i  %i  %g \n', node, DOF, value)
    Eq     = n_g * (node - 1) + DOF     ; % row in system matrix
    M (Eq, Eq) = M (Eq, Eq) + value     ; % add to system mass matrix
  end % for each EBC
% end get_and_add_point_masses    % =================================

function [C] = get_and_add_point_sources (n_g, n_m, C) %=============
% FUNCTION NO LONGER USED HERE
% n_g  = number of DOF per node
% n_m  = number of nodal points in mesh
% C    = system source (force) vector
  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 text file msh_load_pt.tmp')
  end % if user error
  fprintf ('\nRead %i point sources. (Echo of msh_load_pt.tmp) \n', n_u)
  fprintf ('Node, DOF,  Source_value \n')
  for j = 1:n_u                         ; % non-zero Neumann pts
    node  = msh_load_pt (j, 1)          ; % global node number
    if ( node < 1 | node > n_m ) % then impossible node number
      fprint ('Point source has invalid node number = %i \n', node)
      error  ('ERROR, Point source has invalid node number')
    end % if
    DOF   = msh_load_pt (j, 2)          ; % local DOF number
    if ( DOF < 1 | DOF > n_g ) % then impossible DOF number
      fprint ('Point source has invalid DOF number = %i \n', DOF)
      error  ('ERROR, Point source has invalid DOF number')
    end % if
    value = msh_load_pt (j, 3)          ; % point source value
    fprintf ('%i  %i  %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    % =================================

function [S] = get_and_add_point_stiffness (n_g, n_m, S) %==========
% n_g  = number of DOF per node
% n_m  = number of nodal points in mesh
% S    = system stiffness matrix to receive diagonal stiffness terms
  load msh_stiff_pt.tmp            ; % node, DOF, value (eq. number)
  n_u = size(msh_stiff_pt, 1)      ; % number of point stiffnesses
  if ( n_u < 1 )                   ; % missing data
    error ('No stiff_pt data in text file msh_stiff_pt.tmp')
  end % if user error
  fprintf ('\nRead %i point stiffness. (Echo of msh_stiff_pt.tmp) \n', n_u)
  fprintf ('Node, DOF,  Stiffness_value \n')
  for j = 1:n_u                         ; % non-zero stiffness pts
    node  = msh_stiff_pt (j, 1)         ; % global node number
    if ( node < 1 | node > n_m ) % then impossible node number
      fprint ('Point stiff has invalid node number = %i \n', node)
      error  ('ERROR, Point stiff has invalid node number')
    end % if
    DOF   = msh_stiff_pt (j, 2)          ; % local DOF number
    if ( DOF < 1 | DOF > n_g ) % then impossible DOF number
      fprint ('Point stiff has invalid DOF number = %i \n', DOF)
      error  ('ERROR, Point stiff has invalid DOF number')
    end % if
    value = msh_stiff_pt (j, 3)          ; % point stiff value
    fprintf ('%i  %i  %g \n', node, DOF, value)
    Eq         = n_g * (node - 1) + DOF  ; % row in system matrix
    S (Eq, Eq) = S (Eq, Eq) + value      ; % add to system stiff matrix
  end % for each EBC
% end get_and_add_point_stiffness    % =================================

function [C_pt] = get_point_sources (n_g, n_m) %======================
% n_g  = number of DOF per node
% n_m  = number of nodal points in mesh
% C_pt = system point source (force) vector
  n_d  = n_g * n_m                ; % max number of unknowns
  C_pt = zeros (n_d, 1)           ; % point source portion of 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 text file msh_load_pt.tmp')
  end % if user error

  fprintf ('\nRead %i point sources. (Echo of msh_load_pt.tmp) \n', n_u)
  fprintf ('Node, DOF,  Source_value \n')
  for j = 1:n_u                         ; % non-zero Neumann pts
    node  = msh_load_pt (j, 1)          ; % global node number
    if ( node < 1 | node > n_m ) % then impossible node number
      fprint ('Point source has invalid node number = %i \n', node)
      error  ('ERROR, Point source has invalid node number')
    end % if
    DOF   = msh_load_pt (j, 2)          ; % local DOF number
    if ( DOF < 1 | DOF > n_g ) % then impossible DOF number
      fprint ('Point source has invalid DOF number = %i \n', DOF)
      error  ('ERROR, Point source has invalid DOF number')
    end % if
    value = msh_load_pt (j, 3)          ; % point source value
    fprintf ('%i  %i  %g \n', node, DOF, value)
    Eq     = n_g * (node - 1) + DOF     ; % row in system matrix
    C_pt (Eq) = C_pt (Eq) + value      ; % add to system column matrix
  end % for each EBC
% end get_point_sources    % ==========================

function [EBC_flag] = get_ebc_flags (n_g, n_m, PD) % ===========
 EBC_flag = zeros (n_m, n_g)   ; % initialize
 for k = 1:n_m                        ; % loop over all nodes
   if ( PD(k) > 0 )                   ; % at least one EBC here
     [flags] = unpack_pt_flags (n_g, k, PD(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
 if ( n_c == 0 )                     % default to all zero values
   fprintf ('\nNOTE defaulting all EBC values to zero \n')
   return
 else
   fprintf ('Read %i EBC data sets with Node, DOF, Value. \n', n_c)
   fprintf ('(Echo of file text file msh_ebc.tmp) \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 %i & DOF %i. \n', ...
                 node, DOF)
       % EBC_flag (node, DOF) = 1; % try to recover from data error
     end % if common user error
   end % if default to zero
  end % for each EBC

  EBC_count = sum (sum ( EBC_flag == 1 ))   ; % check input data
  if ( EBC_count ~= n_c )                ; % probable user error
    fprintf ('WARNING: mismatch in bc_flag count & msh_ebc.tmp')
  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 ('\n Error missing file msh_typ_nodes.tmp')  
end % if error  
n_n = size (msh_typ_nodes,2) - 1 ; % nodes per element        
fprintf ('\n')
fprintf ('Read %i elements with type & %i nodes each. \n', ...
          n_e, n_n)  
fprintf ('(Echo element number, file msh_typ_nodes.tmp) \n')
% Extract all element type numbers into vector el_type
el_type = round (msh_typ_nodes(:, 1)); % el type number >= 1  
% Extract all element connection lists into array nodes
nodes (1:n_e, 1:n_n) = msh_typ_nodes (1:n_e, 2:1+n_n);  

if ( n_n > 9 ) % skip element number display
 disp(msh_typ_nodes (:, 1:1+n_n)) % echo data , without el number
else
 % Echo data with pretty print from control values
 for i = 1:n_e % loop over all elements ===> ===> ===> ===> ===>
  if ( n_n == 2 ) % nodes per element
   fprintf ('%i, %i %i %i       \n', i, el_type(i), nodes(i,1:n_n))
  elseif ( n_n == 3 )
   fprintf ('%i, %i %i %i %i    \n', i, el_type(i), nodes(i,1:n_n))
  elseif ( n_n == 4 )
   fprintf ('%i, %i %i %i %i %i \n', i, el_type(i), nodes(i,1:n_n))
  elseif ( n_n == 6 )
   fprintf ('%i, %i %i %i %i %i %i %i \n', ...
             el_type(i), nodes(i,1:n_n))
  elseif ( n_n == 8 )
   fprintf ('%i, %i %i %i %i %i %i %i %i %i \n', ...
             el_type(i), nodes(i,1:n_n))
  elseif ( n_n == 9 )
   fprintf ('%i, %i %i %i %i %i %i %i %i %i %i \n', ...
             el_type(i), nodes(i,1:n_n))
  end % if number of nodes to show
 end % for all i elements <=== <=== <=== <=== <=== <=== <=== <=== 
end % if show element number

n_t  = max(el_type)                   ; % number of element types  
if ( n_t > 1 )
  fprintf ('Note: maximum element type is %i \n', n_t)
end % if
% 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.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 ('\n Error missing nodal file msh_bc_xyz.tmp')  
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.tmp) \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))
    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.tmp                       ; % 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 properties \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.tmp) \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_matrix = reshape (Ans, 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                                 ; % save displacements
    if ( n_g == 1 ) 
      fprintf (fid, '%g                \n', Ans_matrix (j, 1:n_g)) ;
      fprintf ('%i   %g                \n', j, Ans_matrix (j, 1:n_g)) ;
    elseif ( n_g == 2 )
      fprintf (fid, '%g %g             \n', Ans_matrix (j, 1:n_g)) ;
      fprintf ('%i   %g %g             \n', j, Ans_matrix (j, 1:n_g)) ;
    elseif ( n_g == 3 )
      fprintf (fid, '%g %g %g          \n', Ans_matrix (j, 1:n_g)) ;
      fprintf ('%i   %g %g %g          \n', j, Ans_matrix (j, 1:n_g)) ;
    elseif ( n_g == 4 )
      fprintf (fid, '%g %g %g %g       \n', Ans_matrix (j, 1:n_g)) ;
      fprintf ('%i   %g %g %g %g       \n', j, Ans_matrix (j, 1:n_g)) ;
    elseif ( n_g == 5 )
      fprintf (fid, '%g %g %g %g %g    \n', Ans_matrix (j, 1:n_g)) ;
      fprintf ('%i   %g %g %g %g %g    \n', j, Ans_matrix (j, 1:n_g)) ;
    elseif ( n_g == 6 )
      fprintf (fid, '%g %g %g %g %g %g \n', Ans_matrix (j, 1:n_g)) ;
      fprintf ('%i %g %g %g %g %g %g \n', j, Ans_matrix (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 [r_q, w_q] = qp_rule_Gauss (n_q) % ================
%  tables of quadrature point locations and weights
%  for lines interpolated in natural coordinates
%             -1 <= r <= 1,   sum (w_q) == 2
 switch n_q
  case 1 % precision: exact for polynomial of degree 1
    r_q (1) = 0 ; w_q = 2 ;

  case 2 % precision: exact for polynomial of degree 3
    r_q (1) = -1 / sqrt (3) ; r_q (2) = 1 / sqrt (3) ;
    w_q (1) = 1             ; w_q (2) = 1 ;

  case 3 % precision: exact for polynomial of degree 5
    r_q (1) = -0.774596669241483377035835 ;
    r_q (2) =  0.000000000000000000000000 ;
    r_q (3) =  0.774596669241483377035835 ;
    w_q (1) =  0.55555555555555555555556 ;
    w_q (2) =  0.88888888888888888888889 ;
    w_q (3) =  0.55555555555555555555556 ;

  case 4 % precision: exact for polynomial of degree 7
    r_q (1) = -0.861136311594052575223946 ;
    r_q (2) = -0.339981043584856264802666 ;
    r_q (3) =  0.339981043584856264802666 ;
    r_q (4) =  0.861136311594052575223946 ;
    w_q (1) =  0.34785484513745385737306 ;
    w_q (2) =  0.65214515486254614262694 ;
    w_q (3) =  0.65214515486254614262694 ;
    w_q (4) =  0.34785484513745385737306 ;

  case 5 % precision: exact for polynomial of degree 9
    r_q (1) = -0.906179845938663992797627 ;
    r_q (2) = -0.538469310105683091036314 ;
    r_q (3) =  0.000000000000000000000000 ;
    r_q (4) =  0.538469310105683091036314 ;
    r_q (5) =  0.906179845938663992797627 ;
    w_q (1) =  0.23692688505618908751426 ;
    w_q (2) =  0.47862867049936646804129 ;
    w_q (3) =  0.56888888888888888888889 ;
    w_q (4) =  0.47862867049936646804129 ;
    w_q (5) =  0.23692688505618908751426 ;

  otherwise % update the tables
   fprintf ('See text for tabulated data where n_q = %i \n', n_q)
   error ('\n ERROR quadrature rule not tabulated for these points')
 end % switch number of quadrature points
% end qp_rule_Gauss % ===============================================

function [r_q, w_q] = qp_rule_unit_Gauss (n_q) % ================
%  tables of quadrature point locations and weights
%  for lines interpolated in unit coordinates. 
%              0 <= r <= 1,   sum (w_q) == 1
 switch n_q
  case 1 % precision: exact for polynomial of degree 1
    r_q (1) = 0.5 ; w_q = 1 ;

  case 2 % precision: exact for polynomial of degree 3
    r_q (1) = 2.1132486540518711774543e-01 ;
    r_q (2) = 7.8867513459481288225457e-01 ;
    w_q (1) = 5.0000000000000000000000e-01 ;
    w_q (2) = 5.0000000000000000000000e-01 ;

  case 3 % precision: exact for polynomial of degree 5
    r_q (1) = 1.1270166537925831148208e-01 ;
    r_q (2) = 5.0000000000000000000000e-01 ;
    r_q (3) = 8.8729833462074168851792e-01 ;
    w_q (1) = 2.7777777777777777777778e-01 ;
    w_q (2) = 4.4444444444444444444445e-01 ;
    w_q (3) = 2.7777777777777777777778e-01 ;

  case 4 % precision: exact for polynomial of degree 7
    r_q (1) = 6.9431844202973712388027e-02 ;
    r_q (2) = 3.3000947820757186759867e-01 ;
    r_q (3) = 6.6999052179242813240133e-01 ;
    r_q (4) = 9.3056815579702628761197e-01 ;
    w_q (1) = 1.7392742256872692868653e-01 ;
    w_q (2) = 3.2607257743127307131347e-01 ;
    w_q (3) = 3.2607257743127307131347e-01 ;
    w_q (4) = 1.7392742256872692868653e-01 ;

  case 5 % precision: exact for polynomial of degree 9
    r_q (1) = 4.6910077030668003601186e-02 ;
    r_q (2) = 2.3076534494715845448184e-01 ;
    r_q (3) = 5.0000000000000000000000e-01 ;
    r_q (4) = 7.6923465505284154551816e-01 ;
    r_q (5) = 9.5308992296933199639881e-01 ;
    w_q (1) = 1.1846344252809454375713e-01 ;
    w_q (2) = 2.3931433524968323402065e-01 ;
    w_q (3) = 2.8444444444444444444444e-01 ;
    w_q (4) = 2.3931433524968323402065e-01 ;
    w_q (5) = 1.1846344252809454375713e-01 ;

  otherwise % update the tables
   fprintf ('See text for tabulated data where n_q = %i \n', n_q)
   error ('\n ERROR quadrature rule not tabulated for these points')
 end % switch number of quadrature points
% end qp_rule_unit_Gauss % ===============================================

function [EBC_react] = recover_reactions_print_save (n_g, n_d, ...
                       EBC_flag, EBC_row, EBC_col, Ans)
% get EBC reaction values by using rows of S & C (before EBC)
  n_d = size (Ans, 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 * Ans - EBC_col     ; % matrix reactions (+-)
  % save reactions (forces) to MODEL file: node_reaction.tmp
  fprintf ('\nReactions at essential BCs \n')         
  fprintf ('Node, DOF, Reaction 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 ('Totals = ') ; disp(Totals)      ; % echo totals
% end recover_reactions_print_save    % ==========================

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 == 1))    ; % 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    % ================================

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 ('\n')                        ; % Skip a line
  % fprintf ('Node, DOF, Resultant Force (1) or Moment (2)  \n')
  fprintf ('Resultant input sources: \n')
  fprintf ('Node, DOF, Resultant input sources \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    % ==========================

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   % =======================================

%                       END SOURCE FILES
%                          EXAMPLE
%   Static_Axial_Analytic
%   
%   (Echo of msh_remarks.tmp) 
%   ============== Begin Application Remarks ====================
%    Two quadratic line element model of a hanging bar, where the
%    top bar is steel : A=10 in^2, E=30e6 psi, gamma=0.283 lb/in^3
%    lower is brass: A=8, E=13e6, gamma=0.300. End load is 10,000 lb
%   ------------------------------------------------------------
%    NOTE: n_g == 1, so index numbers equal node numbers
%   
%     *------(1)-----*------(2)--------*===>> 10,000 lb   ---> x
%     1       2      3       4         5   
%    x=0            x=420             x=660 in.
%   EBC: U_1=0       Point source: P_5=10,000
%   
%   Nodal dof: 1 = axial displacement     
%   Element type = 1 (a parametric line)
%   Element connection: three nodes per element
%   Element properties: 1 = modulus 2 = area 3 = wt density 4=0
%   =============== End Application Remarks =======================
%   
%   User control variables 
%   n_g number of DOF per node        = 1 
%   n_n number of nodes per element   = 2 
%   n_p dimension of parametric space = 1 
%   n_q number of quadrature points   = 5 
%   n_r number of rows in B_e         = 1
%   n_s dimension of physical space   = 1 
%   
%   Begin Phase 2, input mesh, properties, sources and EBC 
%   Read 5 nodes. 
%   (Echo node, file msh_bc_xyz.tmp) 
%   node, bc_flags, 1 coordinates 
%   1, 1    0 
%   2, 0    210 
%   3, 0    420 
%   4, 0    540 
%   5, 0    660 
%   
%   Read 2 elements with type & 3 nodes each. 
%   (Echo element number, file msh_typ_nodes.tmp) 
%   1, 1 1 2 3    
%   2, 1 3 4 5    
%   WARNING Read 3 nodes per element but expected 2 
%   
%   Read 2 materials with 4 properties 
%   4 material values per element 
%   (Echoing text file msh_properties.tmp) 
%      3.0000e+07   1.0000e+01   2.8300e-01            0
%      1.3000e+07   8.0000e+00   3.0000e-01            0
%   
%   Read 1 point sources. (Echo of msh_load_pt.tmp) 
%   Node, DOF,  Source_value 
%   5  1  10000 
%   Note: expecting 1 essential BC values. 
%   Read 1 EBC data sets with Node, DOF, Value. 
%   (Echo of file text file msh_ebc.tmp) 
%        1     1     0
%   
%   Begin Phase 3, assemble (scatter) element arrays 
%   Total input source terms = 11764.6 
%   Begin Phase 4, modify matrices for EBC 
%   Begin Phase 5, solve modified system 
%   
%   Computed Solution: 
%   Node, 1 results per node 
%   1   0                
%   2   0.00802721                
%   3   0.0156384                
%   4   0.0276753                
%   5   0.03938                
%   
%   Reactions at essential BCs 
%   Node, DOF, Reaction Value 
%   1 1 -11764.6 
%   Totals =   -1.1765e+04
%   
%   Begin Phase 6, use the answers 
%   
%   Element Post-processing: 
%   Element number 1 
%   Length ........... 420 
%   Body source ...... 1188.6 
%   Mass ............. 0 
%   Stiffness ........ 714286 
%   Center gradient .. 3.72343e-05 
%   Center flux ...... 1117.03 
%   Left flux ........ 1176.46 
%   Right flux ....... 1057.6 
%   Connectivity 
%        1     2     3
%   Element nodal values 
%               0   8.0272e-03   1.5638e-02
%   Element nodal reactions 
%     -1.1765e+04   5.1159e-12   1.0576e+04
%   
%   Element number 2 
%   Length ........... 240 
%   Body source ...... 576 
%   Mass ............. 0 
%   Stiffness ........ 433333 
%   Center gradient .. 9.89231e-05 
%   Center flux ...... 1286 
%   Left flux ........ 1322 
%   Right flux ....... 1250 
%   Connectivity 
%        3     4     5
%   Element nodal values 
%      1.5638e-02   2.7675e-02   3.9380e-02
%   Element nodal reactions 
%     -1.0576e+04            0   1.0000e+04
%   
%   quit