function [] = Stress_2d_axisym  (DEFAULT) % Revised 3/28/2019
% DEFAULT omitted, require user interactive inputs
% DEFAULT = 1 use bi-quadratic quad (Q9) plane stress elasticity
% Requires addpath /clear/www/htdocs/mech517/Akin_FEA_Lib
%---------------------------------------------------------------------
% Plane stress/axisymmetric elasticity with body, edge & point loads
% with area elements (el_type=1 or 2) for mixed area elements
% and optional edge elements (el_type=3) which must come last
% Plane elements default to unit thickness, 
% Axisymmetric elements default to a unit radian thickness
%---------------------------------------------------------------------
% TYPE 1 or 2 area elements: (T3 and/or Q4, Else T6 and/or Q8 or Q9)
% TYPE 3 edge elements with known traction: (L2 or L3)

% Properties:
% TYPE 1, 2: area elem: 1=E, 2=nu, 3=Body_x, 4=Body_y       
% TYPE 1, 2:            5=Coef_Expansion, 6=temp_increase 
% TYPE 1, 2:            7=optional thickness (default=1)
% TYPE 3: edge tractions:  1=Traction_x, 2=Traction_y,
%                          3=optional thickness (default=1), 4-7=0

% All elements are evaluated by numerical integration
%............................................................. 
%     Mech 417, Finite Element Analysis, Rice University 
%     Copyright J.E. Akin 2019. All rights reserved.
%............................................................. 
% NOTE: input text files begin with msh_ and end with .txt
% see file FEA_notation.m for variable descriptions
% see file FEA_functions.m for functions summary
 
% echo user information about this specific application
echo_user_remarks

% =====================================================================
%                             PHASE 1
%           *** SET OR READ USER PROBLEM DEPENDENT CONTROLS *** 
% =====================================================================
if ( nargin == 0 ) % default to user defined field controls
  DEFAULT = 0 ; % require user inputs
end % if use default application

%  set application controls
if ( DEFAULT == 0 ) % then read user control inputs
  % axisym  = 1 if the problem is axisymmetric
  % load_pt = 1 if point soures are input
  % n_g = number of DOF per node
  % N_N = number of EXPECTED nodes per element
  % n_p = dimension of parametric space
  % n_q = number of quadrature points required
  % n_r = number of rows in element B_e matrix (3 plane stress, 4 axisym)
  % N_S = number of EXPECTED physical space dimensions
  [axisym, load_pt, n_g, N_N, n_p, n_q, n_r, N_S] = get_user_controls ;
  if ( n_g ~= 2 ) % reset
    fprintf ('WARNING: n_g = 2 for plane or axisymmetric elasticity \n')
    n_g = 2 ; % two displacements required here
  end % if

elseif ( DEFAULT == 1 ) % 2D Q9 plane stress elasticity

%  Manually set user maximum array sizes
  n_g     = 2  ; % number of DOF per node (x-, y-displacements)
  N_N     = 9  ; % number of EXPECTED nodes per element
  n_p     = 2  ; % dimension of parametric space
  n_q     = 9  ; % number of quadrature points
  N_S     = 2  ; % EXPECTED physical space dimensions
  axisym  = 0  ; % turn on(1)/off(0) the axisymmetric flag
  load_pt = 1  ; % turn on(1)/off(0) point source (load can be zero)
end % if DEFAULT

%   Manually set user logic choices (or add read function)
color    = 0  ; % turn on(1)/off(0) color plot of results as png file
debug    = 0  ; % turn on(1)/off(0) debug output listings
echo_p   = 1  ; % turn on(1)/off(0) listing of properties
el_react = 0  ; % turn on(1)/off(0) listing of element reactions
est_err  = 0  ; % turn on(1)/off(0) element error estimates
Fourier  = 0  ; % turn on(1)/off(0) minus sign in material law
integral = 0  ; % turn on(1)/off(0) listing of solution integral
mass_ave = 0  ; % turn on(1)/off(0) ave of consistent and diagonal mass
mass_dia = 0  ; % turn on(1)/off(0) use of scaled diagonal mass matrix
mass_pt  = 0  ; % turn on(1)/off(0) point mass (mass can be zero)
post     = 1  ; % turn on(1)/off(0) post-processing option
scp_ave  = 1  ; % turn on(1)/off(0) compute average fluxes at mesh nodes
show_ave = 0  ; % turn on(1)/off(0) list of average fluxes at mesh nodes
show_bc  = 0  ; % turn on(1)/off(0) list assembled S and C, after EBC
show_e   = 0  ; % turn on(>0)/off(0) list some S_e, M_e and C_e
show_err = 0  ; % turn on(1) or off listing of each element error
show_qp  = 0  ; % turn on(1)/off(0) list qp post-processing results
show_s   = 0  ; % turn on(1)/off(0) listing of assembled S and C
stiff_pt = 0  ; % turn on(1)/off(0) point stiffness (stiff can be zero)
vector   = 1  ; % turn on(1)/off(0) displacement vector plot, save

if ( axisym == 1 )  % then axisymmetric solid stress   
  n_r = 4         ; % number of rows in element B_e matrix
else % then plane stress
  n_r = 3         ; % number of rows in element B_e matrix
end % if planar or axisymmetric

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

%    Echo current user logic flags
fprintf ('User logic flags: turn on(1)/off(0) \n') 
fprintf ('axisym   an axisymmetric mesh          = %i \n', axisym  )
fprintf ('color    plot displacement component   = %i \n', color   )
fprintf ('echo_p   listing of properties         = %i \n', echo_p  )
fprintf ('el_react listing of element reactions  = %i \n', el_react)
fprintf ('load_pt  point source (can be zero)    = %i \n', load_pt )
fprintf ('mass_pt  point mass (can be zero)      = %i \n', mass_pt )
fprintf ('post     post-processing option        = %i \n', post    )
fprintf ('show_bc  debug list S & C, after EBC   = %i \n', show_bc )
fprintf ('show_e   debug list S_e, M_e and C_e   = %i \n', show_e  )
fprintf ('show_qp  list qp post-process results  = %i \n', show_qp )
fprintf ('show_s   debug list assembled S and C  = %i \n', show_s  )
fprintf ('stiff_pt point stiffness (can be zero) = %i \n', stiff_pt)
fprintf ('vector   displacement vector plot      = %i \n', vector  )

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

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

% n_m  = number of nodal points in mesh
% n_s  = number of space dimensions
% PD   = nodal packed digit essential BC flag (n_g digits 0 or 1)
%      NOTE: PD is either 00, 10, 01, or 11 here (except for MPC)
% 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 space control
if ( n_s ~= N_S ) % then possible control error
  fprintf ('NOTE: Read %i coordinates, but expected %i \n', n_s, N_S)
end % if
if ( n_s ~= 2 ) % reset
  fprintf ('WARNING: n_s = 2 for plane or axisymmetric elasticity \n')
  n_s = 2   % two displacements required here
end % if
 
% *** Read mesh element types and nodal connections (padded with 0's) ***
% (from msh_typ_nodes.txt, with one line per element)

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

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

%   Prepare a default general 2D mesh plot as png file
plot_input_2d_mesh (x, y, nodes)     ; % display mesh
shrink_2d_mesh     (x, y, nodes)     ; % display mesh
plot_bc_flags      (x, y, nodes, PD) ; % plot BC flags

% Allocate element type control data
type_i = zeros (n_t, 1) ; % maximum DOF per el type
type_n = zeros (n_t, 1) ; % maximum nodes per el type
type_p = zeros (n_t, 1) ; % dimension of type parametric space
type_q = zeros (n_t, 1) ; % number of type quadrature points
type_x = zeros (n_t, 1) ; % number of type geometry nodes (default type_n)

% Default to a single element type
if ( n_t == 1 ) % then initialize
  type_n = n_n ; type_p = n_p ; type_q = n_q ; type_x = n_n ;

else % read controls for each type

% *** Read each element type control integers ***
% (from msh_ctrl_types.txt, with one line per element type)
% type_n = maximum nodes per el type 
% type_p = dimension of type parametric space
% type_q = number of type quadrature points
% type_x = number of type geometry nodes (default type_n)
  [type_n, type_p, type_q, type_x] = get_el_type_controls (n_t) ;

%   check type control vs system control, increase max if needed
  has_n = max ( type_n ) ; % actual max node per el
  has_p = max ( type_p ) ; % actual max parm dim
  has_q = max ( type_q ) ; % actual max Gauss rule
  has_x = max ( type_x ) ; % actual max geom nodes
  if ( has_n > n_n | has_x > n_n ) % warn
    fprintf ('WARNING: type_n = %i greater than n_n = %i \n', has_n, n_n)
    n_n = has_n % increase the maximum number of nodes per element
  end % if
  if ( type_n (1) < n_n ) % warn
    fprintf ('NOTE: type_n(1) = %i < n_n = %i \n', type_n (1), n_n)
  end % if
  if ( has_p > n_p ) % warn
    fprintf ('WARNING: type_p = %i inconsistent with controls \n', has_p)
    n_p = has_p % increase the maximum parametric space dimension
  end % if
  if ( has_q > n_q ) % warn
    fprintf ('NOTE: type_q = %i inconsistent with controls \n', has_q)
    n_q = has_q % increase the maximum number of quadrature points
  end % if
end % if more than one element type
type_i = n_g * type_n ; % number of degrees of freedom per type

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

% n_mats = number of materials, = n_e, or n_t if homogeneous
% n_vals = number of property values for each element or type
% msh_properties (n_mats, n_vals) = system properties array
[n_mats, n_vals, msh_properties] = get_mesh_properties (echo_p, n_e, n_t) ;
fprintf ('\n')
fprintf ('Above properties are for elasticity: \n')
fprintf ('TYPE 1, 2: area elem: 1=E, 2=nu, 3=Body_x, 4=Body_y     \n')
fprintf ('TYPE 1, 2:            5=Coef_Expansion, 6=temp_increase \n')
fprintf ('TYPE 1, 2:            7=optional element thickness      \n')
fprintf ('TYPE 3: edge tractions:  1=Traction_x, 2=Traction_y     \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 matrix (if any) 
% S         = stiffness matrix       

Ans       = zeros (n_d, 1)   ; % results
Body      = zeros (n_g, 1)   ; % initialize volumetric source
C         = zeros (n_d, 1)   ; % force or source
C_pt      = zeros (n_d, 1)   ; % point source portion of C
EBC_flag  = zeros (n_m, n_g) ; % initialize all DOF flags
EBC_value = zeros (n_m, n_g) ; % default EBCs to  zero
el_prop   = zeros (n_vals, 1); % initialize element or type properties
M         = zeros (n_d, n_d) ; % generalized mass
Q_v       = zeros (n_n,   1) ; % variable source nodal values
S         = zeros (n_d, n_d) ; % generalized stiffness
Type_volume = zeros (n_t, 1) ; % volume of each type of element

%   Read point sources & couples data, if any, and insert in C
%   (from file msh_load_pt.txt, one line per point source)
if ( load_pt > 0 )                      ; % need point sources data
  [C_pt] = get_point_sources (n_g, n_m) ; % get point sources
  if ( axisym == 1 ) % then source acts over radius
    for i = 1:n_m
      rows = [1:n_g] + n_g*(i-1)        ; % get rows at point
      C_pt (rows) = C_pt (rows) * x (i) ; % times radius
    end % for
  end % if axisymmetric
  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.txt, 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.txt, one line per point stiffness) 
if ( stiff_pt > 0 )                        ; % need point stiffness data
  [S] = get_and_add_point_stiffness (n_g, n_m, S); % get point stiffness
end % if any point stiffness expected

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

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

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

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

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

%    echo current controls and logic as debugging aid
if ( debug == 1 ) % provide debug checks
  final_control_status (n_d, n_e, n_g, n_i, n_m, n_mats, n_n, ...
                               n_p, n_q, n_r, n_s, n_t, n_vals)
  %L color = 0 ;
  %L final_logic_status (axisym, color, echo_p, el_react, integral,  ...
              %L Fourier, load_pt, mass_ave, mass_dia, mass_pt, post, ...
              %L Show_bc, show_e, show_qp, show_s, stiff_pt, vector)
end % debug listings

% ======================================================================
%                              PHASE 3
%       ** GENERATE ELEMENT VALUES AND ASSYMBLE INTO SYSTEM **
% ======================================================================
fprintf ('Begin Phase 3, assemble (scatter) element arrays \n')
 
%     Open unit needed for post-processing operations
if ( post == 1 ) % open post-processing element unit
  el_unit = fopen ('el_store.bin', 'w')      ; % open for element writing
  if ( el_unit == -1 )                                      % open failed
    fprintf ('\nWARNING: post-processing unit open failed \n')
    post = 0                          ; % turn off post-processing option
  end % if unit available
  qp_unit = fopen ('qp_store.bin', 'w')          ; % open for  GQ writing
  if ( qp_unit == -1 )                                      % open failed
    fprintf ('\nWARNING: post-processing unit open failed \n')
    post = 0   % turn off post-processing option
  end % if unit available
end % if open post-processing unit

%B %    *** Initially fill the type==1 quadrature data files ***
%B % r_q, s_q, t_q = parametric coordinates of quadrature point
%B % w_q           = parametric weight      of quadrature point
%B if ( n_q > 0 )                 % numerical integration to be used
  %B r_q = zeros (n_q, 1) ; s_q = zeros (n_q, 1) ;    % allocate
  %B t_q = zeros (n_q, 1) ; w_q = zeros (n_q, 1) ;    % allocate
  %B [r_q, s_q, t_q, w_q] = get_quadrature_rule (n_n, n_p, n_q) ;
  %B %B [r_q, s_q, w_q] = qp_rule_nat_quad (n_q)
%B end % if numerical integration data required 
 
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

type_old = 0   ; % initialize
for K = 1:n_e  ; % loop over elements ====>> ====>> ====>> ====>> ====>>
  type    = el_type (K)        ; % current element type number
  typ_i   = type_i (type)      ; % maximum DOF per el type
  typ_n   = type_n (type)      ; % maximum nodes per el type
  typ_p   = type_p (type)      ; % dimension of type parametric space
  typ_q   = type_q (type)      ; % number of type quadrature points
  domain  = typ_n * 100 + n_g * 10 + typ_p ; % parametric reference code
  e_nodes = nodes (K, 1:typ_n)             ; % connectivity 

%   Get quadrature rule for current element type
  if ( type ~= type_old ) % then 
    [r_q, s_q, t_q, w_q] = get_quadrature_rule (typ_n, typ_p, typ_q); 
    type_old = type      ; % update type
  end  % if 
   
%                 Gather nodal coordinates
  xy_e (1:typ_n, 1) = x(e_nodes(1:typ_n))   ; % x coord at el nodes 
  xy_e (1:typ_n, 2) = y(e_nodes(1:typ_n))   ; % y coord at el nodes

%     Allocate and clear element type arrays
  C_e      = zeros (typ_i,   1)    ; % clear array el sources
  C_r      = zeros (typ_i,   1)    ; % clear array el reactions
  E_e      = zeros (n_r,   n_r)    ; % material constitutive matrix
  Jacobian = zeros (typ_p, typ_p)  ; % geometry Jacobian (one-to-one)
  S_e      = zeros (typ_i, typ_i)  ; % clear array el stiffness
  strain_0 = zeros (n_r, 1)        ; % initial thermal strains

%     Set constant element type properties
  if ( n_mats > n_t )    % each element has some different properties
    J = K ; % 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) ;

%    Re-name properties
  if (type == 1 | type == 2); % then continuum element ..............>
% Properties: 1=E, 2=nu, 3=Body_x, 4=Body_y, 5=coeff_expansion,
%             6=temp_increase 7=optional thickness
    Body (1:n_g) = 0 ; coeff_e = 0 ; Del_T = 0 ; Thick = 1 ; % defaults
%     Shorthand type property names
    E        = el_prop (1)   ; % modulus of elasticity
    nu       = el_prop (2)   ; % Poisson's ratio

    if ( n_vals > 3 )
      Body (1) = el_prop (3)   ; % x-body force component
      Body (2) = el_prop (4)   ; % y-body force component
    end % if body loads

    if ( n_vals > 5 ) % then initial thermal strain input
      coeff_e = el_prop (5)     ; % coeff thermal expansion
      Del_T   = el_prop (6)     ; % tempature above stress free
      thermal = coeff_e * Del_T ; % initial normal thermal strain
    else
      thermal = 0 ; % no initial normal thermal strain
    end % if thermal load

    if ( n_vals > 6 ) ; % then thickness given
      Thick = el_prop (7)      ; % planar thickness
    end ; % if new thickness

%   populate the stress-strain law matrix
    if ( axisym == 0 ) % then plane stress
      E_e = E / (1 - nu^2) * [ 1,   nu,    0       ; ...
                               nu,  1,     0       ; ...
                               0,   0,  (1 - nu)/2 ] ; % stress-strain law
      strain_0 = thermal * [1  1   0]'  ; % thermal strain 
    else % axisymmetric stress 
      if ( nu >= 0.5 ) % invalid
        error ('Current theory invalid for incompressible material')
      end % if need another theory
      E_e = [(1-nu), nu,    0,          nu ; ...
              nu,   (1-nu), 0,          nu ; ...
              0,     0,     (1-2*nu)/2, 0  ; ...
              nu,    nu,    0,          (1-nu) ] * E / ((1+nu)*(1-2*nu)) ;
      strain_0 = thermal * [1  1  0  1]'  ; % thermal strain 
    end % if application type

  elseif (type == 3); % edge element with tractions .................>
       Body (1) = el_prop (1)   ; % x-traction force component
       Body (2) = el_prop (2)   ; % y-traction force component
       if ( el_prop(3) > 0 ) ; Thick = el_prop (3) ; end ; % optional

  else ; % unknown stress element type ..............................>
       error ('Type number of element exceeds contol value \n')
  end % switch on element type number ..............................<<

%     Allocate and clear quadrature point arrays
  B_q   = zeros (n_r,   typ_i)  ; % initalize B matrix values, at q
  DGH_q = zeros (n_s,   typ_n)  ; % interpolation global derivatives, at q 
  DLH_q = zeros (typ_p, typ_n)  ; % interpolation local derivatives, at q
  Jacobian = zeros (typ_p, n_s) ; % geometry Jacobian
  H_q   = zeros (1,     typ_n)  ; % clear scalar interpolations, at q
  N_q   = zeros (n_g,   typ_i)  ; % vector interpolations, at q
  xy_q  = zeros (1,       n_s)  ; % physical location of point q
  domain = typ_n * 100 + n_g * 10 + typ_p ; % parametric reference code

%   ELEMENT, FOUNDATION STIFFNESSES AND SOURCE RESULTANT MATRICES 

%   Numerical Intergation of S_e, C_e , in 2D parametric
  volume = 0        ; % element volume
  for q = 1:typ_q   ; % begin quadrature loop ---> ---> ---> ---> ---> --->
    r = r_q (q) ; w = w_q (q) ; % recover data
    s = s_q (q) ; t = 0       ; % not used in 1D

%   Element scalar interpolation, H, & local derivatives, DLH
    [H_q, DLH_q] = el_shape_n_local_deriv (domain, r, s, t) ;

% Set the non-zero terms in the vector interpolation
    for k = 1: n_g ;
      N_q (k, k:n_g:typ_i) = H_q (1:typ_n) ;
    end % for shifted rows of scalar interpolations

%   Interpolate global position of quadrature point
    xy_q = H_q * xy_e (1:typ_n, :)      ; % 1 x n_s % global (x, y)

%   Compute the geometry mapping data, d_parm to d_physical
    if ( type < 3 );% area element ...............................>>
      Jacobian = DLH_q * xy_e   ; % n_p x n_s % d(x,y)/d(r,s)
      J_det    = det (Jacobian) ; % 2 x 2     % determinant
    else ;% planar curve .........................................>>
      dx_dr = DLH_q (1, :) * xy_e (1:typ_n, 1) ; % dx (r) / dr
      dy_dr = DLH_q (1, :) * xy_e (1:typ_n, 2) ; % dy (r) / dr
      J_det = sqrt (dx_dr^2 + dy_dr^2)         ; % dL (r) / dr
    end % if type ..............................................<<
              
    if ( J_det <= 0 ) % then bad element geometry
      fprintf ('ERROR negative Jacobian at element %i \n', K)
      fprintf ('used absolute value \n')
      J_det = abs (J_det) ; % assume reverse connection list
    end % if bad element map

    if ( axisym ) % then axisymmetric-unit-radian surface
      J_det = xy_q (1) * J_det      ; % theorem of Papus 
    else ; % planar element
      J_det = Thick * J_det         ; % optional thickness
    end % if axisymmetric body

    if ( type == 1 | type == 2 ) % then continuum element ...........>>
      J_Inv    = inv (Jacobian)       ; % n_s x n_s % inverse

%     Global derivatives of H are DGH
      DGH_q = J_Inv * DLH_q           ; % n_s x n_n   % dH/dx & dH/dy

%     Insert gradient terms into the differential operator
      B_q (1, 1:2:typ_i) = DGH_q (1, 1:typ_n)         ; % for esplion x
      B_q (2, 2:2:typ_i) = DGH_q (2, 1:typ_n)         ; % for esplion y
      B_q (3, 1:2:typ_i) = DGH_q (2, 1:typ_n)         ; % for gamma xy
      B_q (3, 2:2:typ_i) = DGH_q (1, 1:typ_n)         ; % for gamma xy
      if ( axisym ) % then axisymmetric-unit-radian
        B_q (4, 1:2:typ_i) = H_q (1:typ_n) / xy_q (1) ; % for hoop esplion
      end % if

%     Update element square matrices
      S_e = S_e + (B_q' * E_e * B_q) * J_det * w_q (q)      ; % stiffness 

%     Constant body force vector source
      C_e = C_e + N_q' * Body * J_det * w_q (q)             ; % body source

%     Possible initial thermal strain source
      if ( n_vals > 6 ) % then thermal data input
        C_e = C_e + B_q'*E_e*strain_0*J_det*w_q (q) ; % thermal source
      end % if any data

%     Update volume
      volume = volume + J_det * w_q (q)        ; % element volume

%     save data needed for post-processing stresses
      if ( post == 1 ) % save data for this point
        fwrite (qp_unit, xy_q,  'double') ; % save coordinates
        fwrite (qp_unit, E_e,   'double') ; % save material   
        fwrite (qp_unit, B_q,   'double') ; % save strain-displacement
      end % if save point data

    elseif ( type == 3 ) % then a known traction edge ................>>
      C_e = C_e + N_q'*Body*J_det*Thick*w_q (q); % traction resultants

    else ; % invalid type data .......................................>>
      error ('Element arrays not defined for type > 3 \n')
    end % if element type ............................................<<
   end % for quadratures   <--- <--- <--- <--- <--- <--- <--- <--- <--- 

   Type_volume (type) = Type_volume (type) + volume         ;

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

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

%         SCATTER TO (ASSEMBLE INTO) SYSTEM ARRAYS 
% Insert completed element matrices into system matrices 
  [rows] = get_element_index (n_g, typ_n, e_nodes); % eq numbers 
  % rows = vector subscript converting element to system eq numbers
  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 K 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   %g \n', ...
          sum (C(1:n_g:n_i)), sum (C(2:n_g:n_i)) )
%       *** 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
  EBC_col = zeros (EBC_count, 1)          ; % reaction data
  [EBC_row, EBC_col] = save_reaction_matrices (EBC_flag, S, C);
end % if essential BC exist (almost always true)

% in the future, inforce 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);
  fprintf ('NOTE: applied %i essential boundary conditions \n', EBC_count)
end % if essential BC exist (almost always true)

if ( show_bc == 1 ) % then list each system matrix after EBC
  fprintf ('After EBC system matrix S = \n') ; disp(S)
  fprintf ('After EBC system matrix C = \n') ; disp(C)
end % if show_s
% ========================================================================
%                          PHASE 5
%         SOLVE LINEAR EQUATION SYSTEM AND ITS REACTIONS, 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

if ( color == 1 ) % display and save a png plot
  Opt = 0                                    ; % select component
  fprintf ('Plotting deformed mesh \n')
  color_deformed_mesh (x, y, nodes, Ans); % plot deformed shape
  color_displacement  (x, y, nodes, Ans, 0); % magnitude only
end % if plot desired

if ( vector == 1 ) % plot the displacement vectors, save as png
  %B fprintf ('Plotting displacement vectors \n')
  %B quiver_2d_displacements (x, y, nodes, Ans, n_n); % 1 vector per el
  %B quiver_disp_vec_mesh (n_n) ; % plot every n_n th vector, or insert 1
end % if plot desired

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

if ( post == 0 ) % then no post-processing of elements
  fprintf ('\nNo element post-processing requested \n')
else % get the stresses and/or element reactions
  if ( qp_unit > 0 ) ; % close GQ written files
    fclose (qp_unit) ; % file qp_store.bin
    qp_unit = fopen ('qp_store.bin', 'r') ; % open for reading 
  end % if
  if ( el_unit > 0 )                    ; % close written binary files
    fclose (el_unit)                             ; % file el_store.bin
    el_unit = fopen ('el_store.bin', 'r')         ; % open for reading
  end % if element arrays were saved
  if ( el_unit == -1 )                                      % open failed
    fprintf ('\nWARNING: post-processing unit open failed \n')
  end % if unit available
  if ( qp_unit == -1 )                                      % open failed

  end % if unit available

% Save temp file for flux plot input (see Matlab quiver_* plot files)
  q_id = fopen ('el_qp_xyz_fluxes.txt', 'w') ; % flux plot file

  fprintf ('\n')
  fprintf ('Element Post-processing: \n')
  strain = zeros (n_r, 1) ; % pre-allocate strains
  stress = zeros (n_r, 1) ; % pre-allocate stresses
  max_max_VM = 0          ; % max of all von Mises stress
  max_max_MS = 0          ; % max of all max shear stress
  if ( show_qp == 1 ) % get results at the quadrature points
    fprintf ('Output strain order is e_xx, e_yy, e_xy, and e_hoop \n')
    fprintf ('Output stress order is S_xx, S_yy, S_xy, and S_hoop \n')
    fprintf ('Failure criterion are von Mises and maximum shear stress \n')
  end % if show_qp

%  Re-compute element matrices or read them from storage
  max_VM_pt = 0 ; max_VM_el = 0 ; % initialize peak stress location
  max_MS_pt = 0 ; max_MS_el = 0 ; % initialize peak stress location
  type_old = 1                  ; % initialize element type
  skipped  = 0                  ; % number el skipped (for averaging)

  for j = 1:n_e ; % loop over elements ====>> ====>> ====>> ====>> ====>> 
    type   = el_type (j)        ; % current element type number
    typ_i  = type_i (type)      ; % maximum DOF per el type
    typ_n  = type_n (type)      ; % maximum nodes per el type
    typ_p  = type_p (type)      ; % dimension of type parametric space
    typ_q  = type_q (type)      ; % number of type quadrature points
    domain = typ_n * 100 + n_g * 10 + typ_p ; % parametric reference code

    if ( type > 2 ) % then skip this element and go to the next one
      skipped = skipped + 1 ; % reduce elements in error analysis %B
      continue % to the next j ====>> ====>> ====>>
    end % skip this type

%       Gather answers
    e_nodes = nodes (j, 1:typ_n)                      ; % connectivity 
    [rows]  = get_element_index (n_g, typ_n, e_nodes) ; % eq numbers 
    Ans_e (1:n_i) = Ans(rows)                 ; % gather element results

%  Recover previously calculated element arrays
    [xy_e]    = fread (el_unit, [typ_n,   n_s],  'double'); % coordinates
    [el_prop] = fread (el_unit, [n_vals,    1],  'double'); % properties
    [S_e]     = fread (el_unit, [typ_i, typ_i],  'double'); % stiffness
    [C_e]     = fread (el_unit, [typ_i,     1],  'double'); % load result

%  Are there initial thermal strains present ?
    if ( n_vals > 5 ) % then initial thermal strain input
      thermal = el_prop (5) * el_prop (6) ; % alpha * temperature change
      if ( axisym == 1 ) % then axisymmetric normal strains
        strain_0 = thermal * [1  1  0  1]'  ; % thermal strain
      else % plane strain
        strain_0 = thermal * [1  1  0]'     ; % thermal strain
      end % if 
    end % if thermal load

%  Get the strains and stresses at each quadrature point
    for q = 1:n_q   ; % begin quadrature loop ---> ---> ---> ---> ---> --->

%  Recover previous calculations at this point (or re-compute)
      [xy_q] = fread (qp_unit,[1, n_s],    'double'); % get coordinates
      [E_e]  = fread (qp_unit,[n_r, n_r],  'double'); % get material 
      [B_q]  = fread (qp_unit,[n_r, typ_i],'double'); % get operator

%  Compute the generalized strains & stresses
      strain = B_q * Ans_e'               ;  % in-plane strains 
      stress = E_e * (strain - strain_0)  ;  % in-plane stresses 

%     List quadrature point coordinates
      if ( show_qp == 1 ) % list coordinates
        fprintf ('\n')
        fprintf ('El,Q_Pt,Coords %i, %i, %9.3e %9.3e \n', ...
                  j, q, xy_q (1:2))
      end % if list coordinates

      if ( axisym ) % axisymmetric stress analysis, n_r = 4
%     compute the material failure estimates
        vonMise = sqrt( (stress(1)-stress(2))^2 ...
                      + (stress(2)-stress(4))^2 ...
                      + (stress(1)-stress(4))^2 ...
                      + 6 * stress(3)^2 ) * 0.7071     ;
        maxTau = vonMise / sqrt(3) ; % upper bound
        if ( show_qp == 1 ) % list 4 strains, stresses
         fprintf('El,Q_Pt,Strains  %i, %i, %9.3e  %9.3e  %9.3e  %9.3e \n', ...
                   j, q, strain (1:n_r)) % list
         fprintf('El,Q_Pt,Stresses %i, %i, %9.3e  %9.3e %9.3e  %9.3e \n', ...
                   j, q, stress (1:n_r)) % list
        end % if list 4 strains, stress
        % Save coord, stresses, failure criteria for plots
        fprintf (q_id, '%9.3e  %9.3e  %9.3e  %9.3e  %9.3e %9.3e  %9.3e  %9.3e \n', ...
                 xy_q (1:n_s), stress (1:n_r), vonMise, maxTau) ;

      else % plane stress, n_r = 3
%     compute the material failure estimates
        vonMise = sqrt( (stress(1)-stress(2))^2 + stress(2)^2 ...
                       + stress(1)^2 + 6 * stress(3)^2 ) * 0.7071     ;
        maxTau  = sqrt( 0.25 * (stress(1)-stress(2))^2 + stress(3)^2 );
        if ( show_qp == 1 ) % list 3 strains, stresses
         fprintf('El,Q_Pt,Strains    %i, %i, %9.3e  %9.3e  %9.3e \n', ...
                   j, q, strain (1:n_r)) % list
         fprintf('El,Q_Pt,Stresses   %i, %i, %9.3e  %9.3e  %9.3e \n', ...
                   j, q, stress (1:n_r)) % list
        end % if list 3 strains, stress
        % Save coord, stresses, failure criteria for plots
        fprintf (q_id, '%9.3e  %9.3e  %9.3e  %9.3e  %9.3e %9.3e  %9.3e  \n', ...
                 xy_q (1:n_s), stress (1:n_r), vonMise, maxTau) ;
      end % if axisymmetric

      if ( show_qp == 1 ) % list point failure criteria
        fprintf('El,Q_Pt,Von Mises        %i, %i, %9.3e \n',  ...
                  j, q, vonMise)          % list
        fprintf('El,Q_Pt,max Shear Stress %i, %i, %9.3e \n',  ...
                  j, q, maxTau)           % list
      end % if list failure criteria 

      if ( vonMise > max_max_VM ) % update worst value
        max_max_VM = vonMise ; % current value
        max_VM_el  = j       ; % current element
        max_VM_pt  = q       ; % current quadrature point
      end % if current max VM value
      if ( maxTau > max_max_MS ) % update worst value
        max_max_MS = maxTau  ; % current value
        max_MS_el  = j       ; % current element
        max_MS_pt  = q       ; % current quadrature point
      end % if current max MS value

    end % for q quadrature points <--- <--- <--- <--- <--- <--- <--- <---

%   Optional output of nodal reactions on every element
    if ( el_react ) % list nodal forces for this element
      
%     Assign any point source to the first element with that node
      if ( load_pt ) % then point sources to assign (once)
        for eq = 1:typ_i                    % loop over element equations
          k = rows (eq)                   ; % find system DOF number
          if ( C_pt (k) ~= 0 )              % found an input point value
            C_e (eq) = C_e (eq) + C_pt (k); % copy to this node 
            C_pt (k) = 0                  ; % prevent a second use
          end % if first appearance of this point source in the system
        end % for eq node of the current element
      end % if load_pt has been assigned
   
%     Finally, get the element reactions and list them
      C_m = S_e * Ans_e' - C_e  ;
      fprintf ('\n')
      fprintf ('Node, Reactions F_x and F_y \n')
      for k = 1:typ_n
        fprintf ('%i,  %9.3e  %9.3e \n', e_nodes (k), ...
                  C_m (k*n_g-1), C_m (k*n_g))
      end % for
    end % if el_react list

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

  fprintf ('\n')
  fprintf ('==================================================\n')
  fprintf ('Maximum von Mises stress = %9.3e \n', max_max_VM)
  fprintf ('at point %i of element %i \n', max_VM_pt, max_VM_el)
  fprintf ('Maximum shear stress     = %9.3e \n', max_max_MS)
  fprintf ('at point %i of element %i \n', max_MS_pt, max_MS_el)
  fprintf ('==================================================\n')
fclose (q_id) ; % close saved stresses
fprintf ('Saved stress results in el_qp_xyz_fluxes.txt for plots \n')

  if ( axisym ) % then axisymmetric solid
    Type_volume   = 2 * pi() * Type_volume   ;
  end % if
  fprintf ('\n Total domain volume = %9.3e \n', sum (Type_volume))
  format short

% Plot various stress components
addpath /clear/www/htdocs/mech517/Matlab_Plots
if ( vector == 1 ) % plot the principle stress directions
  quiver_qp_2D_stress (1)  ; % Maximum principal stress
  quiver_qp_2D_stress (2)  ; % Minimum principal stress
  fprintf ('Plotted principal stress directions \n')
end ; % if vector plot desired

%B color=0
if ( color == 1 ) % display stress items in color
  color_qp_stresses_step (5) ; % Maximum shear stress
  color_qp_stresses_step (4) ; % vom Mises equivalent
  color_qp_stresses_step (3) ; % XY shear stress
  color_qp_stresses_step (2) ; % YY normal stress
  color_qp_stresses_step (1) ; % XX normal stress
  fprintf ('Plotted selected stress items \n')
end ; % if color desired

%       OPTIONAL POST PROCESSING COMPLETED
end % if post-processing stresses

% =====================================================================
%       PHASE 7: AVERAGE FLUX COMPONENTS AT ALL MESH NODES
%         ( Also needed as input for an error estimate)
% =====================================================================
%B free up memory used by system S & C matrices
if ( scp_ave == 0 & est_err == 0 ); % no need for phase 7 or 8
  fprintf ('No element flux averaging  requested \n')
  fprintf ('No element error estimates requested \n')
  return
end ; % if no need for phase 7 or 8
if ( skipped > 0 ) ; % note reduced continuum element number
  fprintf ('Last %i elements are omitted from averaging \n', skipped)
end ; % if fewer elements
N_E = n_e - skipped ; % actual number of elements for averages
fprintf ('Begin averaging %i fluxes for %i elements \n',n_r, N_E)
fprintf ('      with %i quadrature points \n', n_q)
fprintf ('\n')

% build database of element neighbors needed for node averages
[L_first, L_last] = first_last_L_at_pt (N_E, nodes)              ; % appear
[L_to_L_sum] = count_elems_at_elem (N_E, L_first, L_last, nodes) ; % count
[L_to_L_set] = form_elems_at_el (N_E, L_first, L_last, ...         % set of
                                   nodes, L_to_L_sum)       ; % neighbors
% use patches of neighbors to average flux components at mesh nodes
[scp_ave] = calc_ave_mesh_fluxes (N_E, n_q, n_r, n_s, nodes, x, y, z, ...
                                      L_to_L_sum, L_to_L_set) ;
% get extreme range of each flux component
[range, n_range] = scp_ave_flux_max_min (scp_ave) ; % amount & where
ave_test = 2 * eps                                ; % small value
if ( max (max (abs (range)) < ave_test ))         ; % then warn
  fprintf ('WARNING: all averages are zero. \n')
  fprintf ('         Error estimator will be omitted. \n')
  est_err = 0                                     ; % no data
end % if warn

% list average flux values if desired, save for error estimates
if ( show_ave == 1 )                              ; % list ave fluxes
  list_ave_fluxes (scp_ave)                       ; % at each node
end % if

% save average flux values for plots and/or error estimates
save_ave_fluxes (scp_ave)              ; % to scp_node_ave_fluxes.tmp

if ( color == 1 ) ; % plot the averaged stress items
  color_scp_flux_mag (5) ; % max shear stress
  color_scp_flux_mag (4) ; % von Mises
  color_scp_flux_mag (3) ; % tau xy
  color_scp_flux_mag (2) ; % S_yy
  color_scp_flux_mag (1) ; % S_xx
end ; % if

% =====================================================================
%       PHASE 8: ESTIMATE ERROR IN EACH ELEMENT AND THE SYSTEM
% =====================================================================
if ( est_err == 1 )           ; % estimate errors and new element sizes
  %B percent_err = 1             ; % state desired maximum percent error
  %B [measure, elem_err_energy, elem_refine, error_max] = ...
            %B scp_err_estimates (N_E, n_g, n_p, n_q, n_s, nodes, ...
            %B x, y, z, el_type, scp_ave, Ans, percent_err) ;
  fprintf ('Error estimator being debuged \n')
end                             % if get element error est
% end error estimate and mesh resize

% end Stress_2d_axisym  % ==============================================