function Torsion_2D_RS  
%.............................................................
%      Torsional shear stress in a general bar
%.............................................................
%         Application and element dependent controls
n_g = 1  ; % number of DOF per node 
n_q = 1  ; % number of quadrature points required
n_r = 2  ; % number of rows in B_e matrix

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

%      Extract EBC flags from packed integer flag P 
[EBC_flag] = get_ebc_flags (n_g, n_m, P)      ; % unpack flags
EBC_count  = sum( sum ( EBC_flag > 0 ) )      ; % # of EBC
 
[n_e, n_n, n_t, el_type, nodes] = get_mesh_elements  ;
n_d = n_g*n_m                ; % system degrees of freedom (DOF)
n_i = n_g*n_n                ; % number of DOF per element
S = zeros (n_d, n_d) ; C = zeros (n_d, 1)      ; % initalize sums
M = zeros (n_d, n_d) ;                         ; % initalize sums

%                  Read EBC values, if any
if ( EBC_count > 0 )                            ; % need EBC data
  [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag) ; % read data
end % if any EBC data expected

%    Load the element properties array
load msh_properties.tmp ; % one row per element
n_prop = size(msh_properties, 1) ; % == 1 if homogeneous
fprintf ('\n(Echoing columns of file msh_properties.tmp) \n')
if ( n_prop == 1 )
    G     = msh_properties (1, 1)   ; % elastic shear modulus
    Beta  = msh_properties (1, 2)   ; % twist per unit length

%                   ECHO PROPERTIES
    fprintf ('\nProperties for all elements \n')
    fprintf ('Shear modulus         = %g \n', G)     
    fprintf ('Twist per unit length = %g \n', Beta)  
end % if
if ( n_prop > n_e )
  error ('ERROR: number of property sets exceeds number of elements')
end % if

%     Tabulate quadrature data for unit triangle
n_q = 1                         ; % exact for linear polynomials
r_q (1) = 1./3 ; s_q (1) = 1./3 ; w_q (1)= 1./2       ; % tabulated

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

for j = 1:n_e     ; % element loop ===>> ===>> ===>> ===>> ===>> ===>>
  e_nodes = nodes (j, 1:n_n)                      ; % connectivity
  
%               SET ELEMENT PROPERTIES & GEOMETRY
  if ( n_prop > 1 )
    G     = msh_properties (j, 1)   ; % elastic shear modulus
    Beta  = msh_properties (j, 2)   ; % twist per unit length

%                   ECHO PROPERTIES
    fprintf ('\nProperties for element (by column) %g \n', j)
    fprintf ('Shear modulus         = %g \n', G)     
    fprintf ('Twist per unit length = %g \n', Beta)  
  end % if
  E_e     = [ 1/G,  0 ; ...
              0,  1/G ] ; % shear modulii
  
%                 Gather nodal coordinates
  xy_e (1:n_n, 1) = x(e_nodes)             ; % x coord at el nodes 
  xy_e (1:n_n, 2) = y(e_nodes)   ; % 3 x 2   % y coord at el nodes

%            Numerical Intergation of S_e, C_e 
  S_e = zeros (n_i, n_i) ; C_e = zeros (n_i, 1) ; % initalize values
  B_q = zeros (2, n_i) ; % initalize values

  for q = 1:n_q   ; % begin loop ---> ---> ---> ---> ---> ---> ---> --->
    r = r_q (q) ; s = s_q (q) ; w = w_q (q)     ; % recover data
%                 element type interpolation functions
    H_q = [(1-r-s), r, s]          ; % 1 x 3      % interpolations
    DLH_q = [  -1, 1, 0 ; ...                     % dH/dr
               -1, 0, 1 ]          ; % 2 x 3    ; % dH/ds

%             Global position and Jacobian
    xy_q     = H_q * xy_e           ; % 1 x 2     % global (x, y)
    Jacobian = DLH_q * xy_e         ; % 2 x 2     % d(x,y)/d(r,s)
    J_det    = abs( det (Jacobian)) ; % 1 x 1     % determinant
    J_Inv    = inv (Jacobian)       ; % 2 x 2     % inverse

%             Global derivatives of H
    DGH_q = J_Inv * DLH_q      ; % 2 x 3          % dH/dx & dH/dy
    E_q = E_e                  ; % 2 x 2          % properties
    
%      Insert gradient terms into the differential operator
    B_q = DGH_q     ; % 2 x 3, just an identity for conduction

%                 ELEMENT MATRICES UPDATES, 3 x 3 and 3 x 1
    S_e = S_e + (B_q' * E_q * B_q) * J_det * w_q (q) ; % square
    C_e = C_e + 2 * H_q' * Beta * J_det * w_q (q)    ; % vector
   end % for quadratures   <--- <--- <--- <--- <--- <--- <--- <--- <--- 

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

%    ALLOCATE STORAGE FOR OPTIONAL REACTION RECOVERY          
if ( EBC_count > 0 )                        ; % reactions occur
  [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)
[S, C] = enforce_essential_BC (EBC_flag, EBC_value, S, C);

%             COMPUTE SOLUTION & SAVE
T = S \ C                 ; % Compute stress function
list_save_temperature_results (T)            

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

%       POST-PROCESS ELEMENTS (element nodal shear stress)
fprintf ('\nPost-prossessing Shear Stress & Torque: \n')
fid = fopen('el_qp_xyz_fluxes.tmp', 'w')    ; % open for writing

Torque = 0    ; % zero element torque sums
for j = 1:n_e ; % loop over all elements ====>> ====>> ====>> ====>>
  e_nodes = nodes (j, 1:n_n)                      ; % connectivity 
  [rows] = get_element_index (n_g, n_n, e_nodes)  ; % eq numbers 
  T_e (1:n_i) = T(rows) ;         % gather element displacements 
   
%               SET ELEMENT PROPERTIES & GEOMETRY
  if ( n_prop > 1 )
    G     = msh_properties (j, 1)   ; % elastic shear modulus
    Beta  = msh_properties (j, 2)   ; % twist per unit length
  end % if
  E_e     = [ 1/G,  0 ; ...
              0,   1/G ] ; % shear modulii          
  
%                 Gather nodal coordinates
  xy_e (1:n_n, 1) = x(e_nodes)             ; % x coord at el nodes 
  xy_e (1:n_n, 2) = y(e_nodes)   ; % 3 x 2   % y coord at el nodes

%            Numerical Intergation of S_e, C_e 
  S_e = zeros (n_i, n_i) ; C_e = zeros (n_i, 1) ; % initalize values
                           B_q = zeros (2, n_i) ; % initalize values

  for q = 1:n_q   ; % begin loop  ---> ---> ---> ---> ---> ---> --->
    r = r_q (q) ; s = s_q (q) ; w = w_q (q)     ; % recover data
%                 element type interpolation functions
    H_q = [(1-r-s), r, s]          ; % 1 x 3      % interpolations
    DLH_q = [  -1, 1, 0 ; ...                     % dH/dr
               -1, 0, 1 ]          ; % 2 x 3    ; % dH/ds

%             Global position and Jacobian
    xy_q     = H_q * xy_e           ; % 1 x 2     % global (x, y)
    Jacobian = DLH_q * xy_e         ; % 2 x 2     % d(x,y)/d(r,s)
    J_det    = abs( det (Jacobian)) ; % 1 x 1     % determinant
    J_Inv    = inv (Jacobian)       ; % 2 x 2     % inverse

%             Global derivatives of H
    DGH_q = J_Inv * DLH_q      ; % 2 x 3          % dH/dx & dH/dy
    E_q = E_e                  ; % 2 x 2          % properties

%      Output point location, shear stress vectors
    fprintf ('\nElement, Point, Coordinates         %i %i %g %g \n', ...
                j, q, xy_q (1:2))
    Gradient = DGH_q * T_e' ;         % standard gradient
    Tau_ZX  = Gradient (2)  ;         % shear stress component
    Tau_ZY = -Gradient (1)  ;         % shear stress component
    Shear_vector = [Tau_ZX, Tau_ZY] ; % net vector in cross-section
    fprintf ('Element, Point, Shear Stress Vector %i %i %g %g \n', ...
                j, q, Shear_vector (1:2))
    fprintf (fid, '%g %g %g %g \n',  xy_q (1:2), Shear_vector (1:2)); % save
 
%     Increment Torque total
    Torque = Torque + 2 * H_q * T_e' * J_det * w ; % twice integral of phi
    
show = 0 ; % skip the showing of element reaction results
if ( show ) % list reactions
%      Insert gradient terms into the differential operator
    B_q = DGH_q     ; % 2 x 3, just an identity for conduction

%                 ELEMENT MATRICES UPDATES, 3 x 3 and 3 x 1
    S_e = S_e + (B_q' * E_q * B_q) * J_det * w_q (q) ; % square
    C_e = C_e + 2 * H_q' * Beta * J_det * w_q (q)           ; % vector
end % if show
   end % for quadratures <--- <--- <--- <--- <--- <--- <--- <--- <--- 
 
if ( show )
  fprintf ('\nNode List for Element %g, \n', j) 
  disp (e_nodes)
  fprintf ('Given Resultant Sources \n') 
  disp (C_e') 
 
% Finally, get the reactions 
  C_m = S_e * T_e' - C_e  ;
  fprintf ('Heat Flow Reactions \n') 
  disp (C_m') 
end % if show
 
end % for each j element in mesh  <<==== <<==== <<==== <<==== <<==== 
fprintf ('Total section torque = %g \n', Torque)

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

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

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

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

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

function [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag)
 EBC_value = zeros(n_m, n_g)  ; % initialize to zero
  load msh_ebc.tmp            ; % node, DOF, value (eq. number)
  n_c = size(msh_ebc, 1)      ; % number of constraints
  fprintf ('\nApplied Displacement Boundary Conditions: %g \n', n_c)
  fprintf ('(Echo of file load msh_ebc.tmp) \n')               
  %b fprintf ('Node, DOF (1=displacement, 2=slope), Value. \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.tmp \n')
  end % if user error
% end get_ebc_values

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

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

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

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

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

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

function list_save_temperature_results (T)
  n_m = size (T, 1)                          ; % get size
  fprintf('Temperature at %g nodes \n', n_m) ; % header

  % save results (temperature) to MODEL file: node_results.tmp
  fid = fopen('node_results.tmp', 'w')       ; % open for writing
  for j = 1:n_m                              ; % save temperature
    fprintf ( fid, '%g \n', T (j))           ; % print
    fprintf (' %g %g \n', j, T (j))          ; % sequential save
  end % for j DOF
% end list_save_temperature_results (T)

function output_T3_heat_flux (n_e, n_g, n_n, n_q, nodes, x, y, T)
%       POST-PROCESS ELEMENT HEAT FLUX RECOVERY & SAVE
fid = fopen('el_qp_xyz_fluxes.tmp', 'w')    ; % open for writing
fprintf ('\n')                              ; % blank line
fprintf('Elem, X_qp, Y_qp, HeatFlux_x, HeatFlux_y \n');% header

for j = 1:n_e                      ; % loop over elements   ====>>
  e_nodes = nodes (j, 1:n_n)               ; % connectivity
  [a, b, c, center, two_A] = form_T3_geom_constants (x, y, e_nodes);

%  get DOF numbers for this element, gather solution
  [rows] = get_element_index (n_g, n_n, e_nodes) ; % eq numbers
  T_e = T (rows)                          ; % gather element DOF

  for q = 1:n_q        ; % Loop over element quadrature points ---->

%   H_i (x,y) = (a_i + b_i*x + c_i*y)/two_A    % interpolations
    B_e (1, 1:3) = b(1:3) / two_A            ; % dH/dx
    B_e (2, 1:3) = c(1:3) / two_A            ; % dH/dy

%  COMPUTE GRADIENT & HEAT FLUX, SAVE LOCATION AND VALUES
    Gradient = B_e * T_e                        ; % gradient vector
    HeatFlux = E_e * Gradient                   ; % heat flux vector
    fprintf (fid, '%g %g %g %g \n',  center(1:2), HeatFlux(1:2));% save
    fprintf ('%g %g %g %g %g \n', j, center(1:2), HeatFlux(1:2));% prt
  end % for loop over n_q element quadrature points          <----
end % for each j element in mesh                            <<====
% end output_T3_heat_flux 

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

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

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

function [flags] = unpack_pt_flags (n_g, N, flag)
% unpack n_g integer flags from the n_g digit flag at node N
% integer flag contains (left to right) f_1 f_2 ... f_n_g
 full = flag                         ; % copy integer
 check = 0                           ; % validate input
 for Left2Right = 1:n_g              ; % loop over local DOF at k
   Right2Left = n_g + 1 - Left2Right ; % reverse direction
   work = floor (full / 10)          ; % work item
   keep = full - work * 10           ; % work item
   flags (Right2Left) = keep         ; % insert into array
   full = work                       ; % work item
   check = check + keep * 10^(Left2Right - 1) ; % validate
 end % for each local DOF
 if ( flag > check )                 ; % check for likely error
   fprintf ('WARNING: bc flag likely reversed at node %g. \n', N)
 end % if likely user error
% end unpack_pt_flags
 
%=================== Running gives =====================================