function Torsion_L3_masses(n_e, show)
% Axial bar torsion & lumped masses: J(x)G*U_xx = U_tt*J(x)*gamma
% U=Angle of twist per unit length. G = shear modulus, 
% gamma = wt. density, t=time, J=polar mass inertia 
% Ref: J.E. Akin, "FEA w Error Estimators" Elsevier, 2005, page 50  
%.....................................................................
%      x_1=0  x_2   x_3    x_4    x_5= L
% Fixed *-----(1)----*-----(2)-----*  Bit_lumped_J 
%       U_1   U_2   U_3    U_4    U_5
% Connectivity list:    e     i      j     k
%                       1     1      2     3  % L3
%                       2     3      4     5  % L3
%.....................................................................
if ( nargin == 1 ) ; 
  show = 0 ;             % flag to show exact solution (if = 1)
elseif ( nargin == 0 ) ; 
  show = 0 ; n_e = 2 ;   % n_e = number of elements
end % if from argument count

%            Set given constants (try changing n_e)
n_g = 1 ;                  % number of DOF per node
n_n = 3 ;                  % number of nodes per element
n_i = n_n*n_g           ;  % number of DOF per element
n_m = 1 + n_e*(n_n - 1) ;  % number of mesh nodes
n_d = n_g*n_m           ;  % system degrees of freedom (DOF)

% specific problem: Thomson Ex 56.2 drill string & bit
L = 5000.                     ; % domain lengths
G   = 1.728e9 ; gamma = 15.22 ; % shear mod., wt. density
J_e = 9.4e-4  ; J_bit = 29.3  ; % string and bit inertia
L_e = L/n_e                   ; % uniform element length
X = [0:1:(n_d-1)]*L/(n_d-1)     % merge lengths

% Allocate (wasteful) storage for BC flags and values
EBC_flag  = zeros(n_m,n_g) ; % EBC must be zero in value
Mass_flag = zeros(n_m,n_g) ; Mass_value = zeros(n_m,n_g) ;

% Set specific essential boundary conditions & lumped masses
EBC_flag  (1,   1) = 1   % EBC must be zero in value     
Mass_flag (n_m, 1) = 1   
Mass_value (n_m, 1) = J_bit % at end
EBC_count  = sum(sum(EBC_flag))   ; % number of EBC's
Mass_count = sum(sum(Mass_flag))  ; % number of Masses

%     Constant quadratic element square matrices
K_e = G * J_e *     [ 7, -8,  1, ; ...      % stiffness L3
                     -8, 16, -8, ; ...      %
                      1, -8,  7 ] / (3*L_e)  
M_e = J_e*gamma*L_e*[4, 2, -1 ; ...         % mass L3
                     2, 16, 2 ; ...         %
                    -1,  2, 4 ] / 30         

%         Assemble two n_d by n_d square matrix terms
S = zeros (n_d, n_d) ;  M = zeros (n_d, n_d) ;
for j = 1:n_e                           ; % loop over elements
  rows = [1:n_n] + (j - 1)*(n_n - 1)    ; % get connectivity 
  S (rows, rows) = S (rows, rows) + K_e ; % add to stiffness
  M (rows, rows) = M (rows, rows) + M_e ; % add to mass      
end % for element j
disp(S)
disp(M)

%  Apply essential BC at nodes, via trick to avoid partitions
EBC_flag   = reshape (EBC_flag,   n_d, 1) ; % change to vector 
Mass_flag  = reshape (Mass_flag,  n_d, 1) ; % change to vector 
Mass_value = reshape (Mass_value, n_d, 1) ; % change to vector 

for j = 1:n_d                % check all DOF for lumped masses
  if ( Mass_flag(j) )                       % lumped mass here
    M (j, j) = M (j, j) + Mass_value (j) ;  % add to diagonal
  end % if Mass for this DOF

  if ( EBC_flag(j) )                       % then EBC = 0 here 
    S(:,j) = 0 ; S(j,:) = 0 ; S(j,j) = 1 ; % trick no partition
    M(:,j) = 0 ; M(j,:) = 0 ; M(j,j) = 1 ; % trick no partition
  end % if EBC for this DOF
end % for over all DOF
disp(S)
disp(M)

%   Solve for eigenvalues and eigenvectors, general form
[Vec, Diag] = eig(S,M) ; % eigen-vector cols, -value^2 diagonal
disp(Diag)
disp(Vec)

%  Sort in assending order and scale for output
for j = 1:n_d % DOFs
  OMEGA (j) = sqrt(real(Diag (j, j))) ;
end % for all DOF
omega = sort(OMEGA)   % in assending order, but with EBC's
disp(omega)

%   Reorder to remove EBC to reduce effective system size
kount = EBC_count ; ORDER = [1:n_d] ;
for j = 1:n_d ;
  if ( kount > 0 ) % inactive EBC DOF remain to shift
    for k = 1, n_d-1 ;
      if ( EBC_flag (k) > 0 ) % then EBC & inactive DOF
        temp            = ORDER (k)       ; % an EBC DOF
        ORDER (k:n_d-1) = ORDER (k+1:n_d) ; % shift up others
        ORDER (n_d)     = temp            ; % EBC to last
      end % fake DOF
    end % for k sorting EBC to end of list
    kount = kount - 1 ;
  end % if kount
end % for j
disp(ORDER)

% ---------------------------------------------------------
% End finite element calculations. ---- Begin graphics ----
% ---------------------------------------------------------
 
clf ; hold on ; grid on                    % clear & keep
if ( Mass_count == 0 ) 
  title ('Modes of elastic bar torsion')
  % e_text = ('sqrt (G/(gamma * L^2))*[1.571 4.712 7.854 10.99 ...]');
  e_text = ('Analytic Omega = [3.3479  10.0415  16.7373  23.4203]');
  text ((X(1)+X(2))/3, -0.25, e_text)

else
  title ('Modes of elastic bar torsion with lumped masses')
  fprintf ('Modes of elastic bar torsion with lumped masses \n')
  e_text = ('Analytic drill string/bit Omega = 2.41, 7.93 rad/s') ;
  text ((X(1)+X(2))/3, -0.25, e_text)
end % if

%            Show BC locations on the plot
k_EBC = 0 ; k_Mass = 0      ; % counts = 0
for j = 1:n_d % check DOF for EBC or lumped mass
  node   = 1 + (j - 1)/n_g  ; % node at DOF j
  U(1:n_d) = Vec (1:n_d, j) ; % get 1 eigen vector
  if ( EBC_flag(j) )          % then EBC here
    k_EBC = k_EBC + 1       ; % increment count
    b_text = sprintf ('--EBC (%g)', k_EBC) ;
    plot (X(node), 0.0, 'k*')   % display flag         
    text (X(node), 0.0, b_text) % display text         
  end % if EBC at this node
  if ( Mass_flag(j) )           % then Mass here
    k_Mass = k_Mass + 1       ; % increment count
    b_text = sprintf ('--Mass (%g)', k_Mass) ;
    plot (X(node), 0.0, 'k*')   % display flag         
    text (X(node), 0.0, b_text) % display text         
  end % if lumped mass at this node
end % for over all DOF

%  ------------------------------------------------------
%  graph of value in 3 node quadratic line elements
%  ------------------------------------------------------
% e_x, e_y   = coordinates & data values at 3 element nodes
% x_el, y_el = interpolated parametric values

 n_poly = ceil (75/n_e)       ; % curve segments
 e_x(3)=0 ; e_y(3)=0 ; H(3)=0 ; % pre-allocation

%  Loop over (physical) mode shapes
  Limit = min([(n_d-EBC_count), 3])   ; % limit # of plots <<==
  for j = 1:Limit                     ; % mode number
    NOW = ORDER (j)                   ; % sorted order
    scale = max(abs(Vec(1:n_d, NOW))) ; % range -1 to +1
    U = Vec(1:n_d, NOW) / scale       ; % eigenvector

    for it = 1:n_e ; %   Loop over all elements
    
%     Extract nodal coordinates & values. Plot nodes.
      rows = [1, 2, 3] + (it - 1)*(n_n-1); % get connectivity 
      e_x  = X (rows) ; e_y  = U (rows) ;  % gather element data
      plot (e_x, e_y, 'ko')                % nodal symbols
      if (it == 1)                         % one summary per mode
        e_text = sprintf(' --(Mode %g) Omega = %g', j, omega(NOW)); 
        text (e_x(2), e_y(2), e_text)      % Plot element number
        plot (e_x(2), e_y(2), 'k*')        % Plot element number
      end % if

%     Loop over local points on the quadratic element
      for k = 1: (n_poly + 1) % points in parametric space
%       Get parametric interpolation functions
        Z    = 2*(k - 1)/n_poly - 1            ; % on -1 to 1
        H(1) = 0.5*(Z*Z - Z) ; H(2) = 1. - Z*Z ;
        H(3) = 0.5*(Z*Z + Z)                   ;
        x_el (k) = dot(H,e_x) ; y_el (k) = dot(H,e_y) ; % true value
      end % for k parametric points in element

      if ( mod(j, 2) == 0 )     % then alternate mode shape colors
        plot (x_el, y_el, 'k-') % parametric curve through 3 nodes
      else
        plot (x_el, y_el, 'b-') % parametric curve through 3 nodes
      end % if color choice

    end % for it over all elements
  end % for j-th mode
  xlabel (['Position (of ', int2str(n_e), ' elements)']) ; 
  ylabel (['U(x) for ', int2str(n_m), ' nodes'])         
% end of Torsion_L3_masses

%%            Plot the exact mode results
%if ( Mass_count == 0 ) 
%  z = [0:0.025:1] * L ; root = (j-1+0.5)*pi*sqrt((L^2*gamma)/G) ; 
%  exact =sin(z*root)                             ; % analytic
%  plot (z, exact, 'r.-');  plot (X, U, 'ko')       % exact & nodal
%  legend ('Exact', 'FEA'); 
%end % if no lumped mass

%     if ( j == abs(show) ) % show one exact mode
%       exact = sign(show)*cos(x_el*(j-1)*root); % Thomson "Mech Vib"
%       plot (x_el, exact, 'r-') % exact mode shape at parametric pts
%       if ( it == 1 ) % then for first element only
%         omega = (j-1)*pi ;
%         e_text = sprintf ('----(Mode %g) Exact = %g', j, omega); 
%         text (x_el(9), exact(9), e_text)     % Plot 9-th pt
%       end % if first element
%     end % if show exact
%  +++++++++++++ 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 [a, b, c, center, two_A] = form_T3_geom_constants (x, y, e_nodes)
% Planar 3 node triangle geometry: H_i (x,y) = (a_i + b_i*x + c_i*y)/two_a
%                   define nodal coordinates, ccw: i, j, k
  x_e = x(e_nodes) ; y_e = y(e_nodes)        ; % coord at el nodes
  x_i = x_e(1) ; x_j = x_e(2) ; x_k = x_e(3) ; % change notation
  y_i = y_e(1) ; y_j = y_e(2) ; y_k = y_e(3) ; % change notation

%  define centroid coordinates (quadrature point)
  center (1) = (x_i + x_j + x_k)/3 ;  
  center (2) = (y_i + y_j + y_k)/3 ;

%  geometric parameters: H_i (x,y) = (a_i + b_i*x + c_i*y)/two_a
  a_i = x_j * y_k - x_k * y_j ; b_i = y_j - y_k ; c_i = x_k - x_j ;
  a_j = x_k * y_i - x_i * y_k ; b_j = y_k - y_i ; c_j = x_i - x_k ;
  a_k = x_i * y_j - x_j * y_i ; b_k = y_i - y_j ; c_k = x_j - x_i ;

  a (1:3) = [a_i, a_j, a_k] ;
  b (1:3) = [b_i, b_j, b_k] ;
  c (1:3) = [c_i, c_j, c_k] ;

%  calculate twice element area
  two_A = a_i + a_j + a_k            ; % = b_j*c_k - b_k*c_j also
% end form_T3_geom_constants (x, y, e_nodes)

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 ('Read %g point sources. \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
    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
  fprintf ('\n')
% 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 ('Read %g EBC data with Node, DOF, Value. \n', n_c)
  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 ('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 (pre_e) ;
% MODEL input file controls (for various data generators) 
if (nargin == 0)  ; % default to no proceeding items in data
  pre_e = 0       ; % Dummy items before el_type & connectivity 
end % if 

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) - pre_e - 1 ; % nodes per element       
fprintf ('Read %g elements with type number & %g nodes each. \n', ...
          n_e, n_n) 
el_type = round (msh_typ_nodes(:, pre_e+1)); % el type number >= 1 
n_t     = max(el_type)                   ; % number of element types 
fprintf ('Maximum number of element types = %g. \n', n_t) 
nodes (1:n_e, 1:n_n) = msh_typ_nodes (1:n_e, (pre_e+2:pre_e+1+n_n)); 
disp(msh_typ_nodes (:, (pre_e+1:pre_e+1+n_n))) % echo data
% end get_mesh_elements 

function [n_m, n_s, P, x, y, z] = get_mesh_nodes (pre_p) ;
% MODEL input file controls (for various data generators) 
if (nargin == 0) % override default 
  pre_p = 0 ; % Dummy items before BC_flag % coordinates 
end % if 
 
%       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) - pre_p - 1 ; % number of space dimensions 
fprintf ('Read %g nodes with bc_flag & %g coordinates. \n', n_m, n_s) 
msh_bc_xyz (:, (pre_p+1))= round (msh_bc_xyz (:, (pre_p+1)));
P = msh_bc_xyz (1:n_m, (pre_p+1)) ; % integer Packed BC flag 
x = msh_bc_xyz (1:n_m, (pre_p+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, (pre_p+3))           ; % extract y column  
end % if 2D or 3D
if ( n_s == 3 )                               ; % check 3D
  z = msh_bc_xyz (1:n_m, (pre_p+4))           ; % extract z column 
end % if 3D 
disp(msh_bc_xyz (:, (pre_p+1):(pre_p+1+n_s))) ; % echo data
% end get_mesh_nodes 

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

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

function 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_PlaneStress_stresses(n_e, n_g, n_n, n_q, nodes, x,y,T)
%       POST-PROCESS ELEMENT STRESS RECOVERY & SAVE
fid = fopen('el_qp_xyz_fluxes.tmp', 'w')    ; % open for writing
fprintf ('\n')                              ; % blank line
fprintf('Elem, QP, X_qp, Y_qp \n')                         ;% header
fprintf('Elem, QP, Stress_qp: xx        yy          xy \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);
  [t_e, Body_e, E_e] = set_constant_plane_stress_prop; % properties

%  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:2:5) = b (1:3)/two_A ; B_e (2, 2:2:6) = c (1:3)/two_A ;
    B_e (3, 1:2:5) = c (1:3)/two_A ; B_e (3, 2:2:6) = b (1:3)/two_A ;

%  COMPUTE GRADIENT & HEAT FLUX, SAVE LOCATION AND VALUES
    Strains = B_e * T_e                      ; % mechanical strain
    Stress = E_e * Strains                   ; % mechanical stress
    fprintf (fid,'%g %g %g %g %g  \n',  center(1), center(2), ...
             Stress(1), Stress(2), Stress(3));% save
    fprintf ('%g %g %g %g  \n',    j, q, center(1:2));% prt
    fprintf ('%g %g %g %g %g  \n', j, q, Stress(1:3));% prt
    fprintf ('\n')                                   ;% prt
  end % for loop over n_q element quadrature points          <----
end % for each j element in mesh  
% end output_PlaneStress_stresses (n_e, n_g, n_n, n_q, nodes, x, y, 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);
  [t_e, Body_e, E_e] = set_constant_plane_stress_prop; % properties

%  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 (n_e, n_g, n_n, n_q, nodes, x, y, T)

function [EBC_react] = recover_reactions_print_save (n_g, n_d, ...
                       EBC_flag, EBC_row, EBC_col, T)
% get EBC reaction values by using rows of S & C (before EBC)
  n_d = size (T, 1)                     ; % number of system DOF
% n_c x 1 = n_c x n_d  * n_d x 1 +  n_c x 1
  EBC_react = EBC_row * T - EBC_col     ; % matrix reactions (+-)
  % save reactions (forces) to MODEL file: node_reaction.tmp
  fprintf ('\n')                        ; % Skip a line
  fprintf ('Node, DOF, Value for %g Reactions \n', ...
            sum (sum (EBC_flag > 0)))   ; % header
  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 (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 ('\n')                        ; % Skip a line
  % fprintf ('Node, DOF, Resultant Force (1) or Moment (2)  \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 (n_g, n_m, C)

function [I, E, Rho, Line_e] = set_constant_beam_prop (n_n, Option) ; 
  if ( nargin == 1 )
    Option = 1 ;
  elseif ( nargin == 0 )
    n_n = 2 ; Option = 1 ;
  end % if problem Option number
  Line_e = zeros (n_n, 1)         ; % default line load at nodes
  switch  Option
    case 1 % Propped cantilever with uniform load, L, L/4
           % Wall reactions: V=37*Line*L/64 M=7*Line*L^2/64 
           % Roller reaction: R= 43*Line*L/64
           % Total vertical load: 5*Line*L/4
      %       *-----(1)-----*-----(2)-----*--(3)--*     EI
      % Fixed@1     L/2     2      Roller@3  L/4  4
      I = 1.0 ; E = 1.0 ; Rho = 0.0 ; Line_e = [1.0; 1.0] ;
    case 2 % cantilever with uniform load, L
      I = 1.0 ; E = 1.0 ; Rho = 0.0 ; Line_e = [1.0; 1.0] ;
    otherwise
      I = 1.0 ; E = 1.0 ; Rho = 0.0; % default shape & material
  end % switch
% end set_constant_beam_prop; 

function [t_e, Body_e, E_e] = set_constant_plane_stress_prop ;
  t_e = 1 ; Body_e (1:2) = 0.             ; % defaults
%                         case 1
  t_e          = 5e-3                     ; % thickness
  Body_e (1:2) = [5e5, 0.]                ; % components
  E   = 15e9                              ; % Elastic modulus
  nu  = 0.25                              ; % Poisson's ratio
%                    plane stress 
  E_v = E/(1 - nu^2)                            ; % constant
  E_e (1, 1) = E_v      ; E_e (1, 2) = E_v * nu ; % non-zero term
  E_e (2, 1) = E_v * nu ; E_e (2, 2) = E_v      ; % non-zero term
  E_e (3, 3) = E_v * (1 - nu) / 2               ; % non-zero term
%end set_constant_plane_stress_prop 

function [t_e, Q_e, E_e] = set_constant_2D_conduction_prop
 % Manually set constant element properties (Fig 11.9 text)
  Q_e = 0. ; t_e = 1.                     ; % defaults 
% case 1
  Kx  = 8. ; Ky = 8. ; Kxy = 0.           ; % thermal conductivity
% case 2
  kx  = 1. ; Ky = 1. ; Kxy = 0. ;
% insert
  E_e = zeros (2, 2)                      ; % constitutive matrix
  E_e (1, 1) = Kx   ; E_e (1, 2) = Kxy    ; % non-zero term
  E_e (2, 1) = Kxy  ; E_e (2, 2) = Ky     ; % non-zero term
% end set_constant_2D_conduction_prop

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