%           OPTIONAL: POST-PROCESS ELEMENTS
if ( post == 1 ) % then post-process all element stresses
  if ( qp_unit > 0 ) ;                   % close written files
    fclose (qp_unit) ;                   % file qp_store.bin
    qp_unit = fopen ('qp_store.bin', 'r') ; % open for reading 
  end % if
  qid = fopen ('el_qp_xyz_fluxes.tmp', 'w') ; % plot file

  fprintf ('\nElement Post-processing: \n')
  fprintf ('Point       x-, y-coordinates \n')
%bfprintf ('            3 strains:  e_xx, e_yy, e_xy  \n')
  fprintf ('            3 stresses: s_xx, s_yy, s_xy  \n')
  fprintf ('            von Mises, max shear \n')
  strain = zeros (n_r, 1) ; % pre-allocate strains
  stress = zeros (n_r, 1) ; % pre-allocate stresses

%  Re-compute or read arrays from storage
  format short e
  for j = 1:n_e ; % loop over elements ====>> ====>> ====>> ====>> ====>> 
    fprintf ('Element number %i _______________________ \n', j)
    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 
   
    for q = 1:n_q  ; % begin quadrature loop ---> ---> ---> ---> ---> --->
      [xy_q] = fread (qp_unit, [1, n_s],   'double');% recover coordinates
   %b [H_q]  = fread (qp_unit, [n_n, 1],   'double');% get interpolation
      [E_q]  = fread (qp_unit, [n_r, n_r], 'double');% recover material
      [B_q]  = fread (qp_unit, [n_r, n_i], 'double');% recover strain-disp

%     compute the strains & stresses, and material failure estimates
      strain = B_q * T_e'                      ;  % in-plane strains 
      stress = E_q * strain                    ;  % in-plane stresses
      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 );

      fprintf ('%i, %g, %g \n',     q, xy_q (1:n_s))    % list
%b    fprintf ('%i, %g, %g, %g \n', q, strain (1:n_r))  % list
      fprintf ('%i, %g, %g, %g \n', q, stress (1:n_r))  % list
      fprintf ('%i, %g, %g \n',     q, vonMise, maxTau) % list
      fprintf (qid, '%g, %g, %g, %g, %g, %g, %g \n', ... % save for plots
               xy_q (1:2), stress (1:3), vonMise, maxTau) ;
    end % for q quadrature points <--- <--- <--- <--- <--- <--- <--- <---

  end % for each j element in mesh   <<==== <<==== <<==== <<==== <<==== 
  format short
%       OPTIONAL POST_PROCESSING COMPLETED
else
  fprintf ('\nNo element post-processing requested \n')
end % for post-processing stresses