% Element matrices for [d/dx (K*A du/dx)] - h*P (u - u_inf) + Q = 0
%     with properties 1=K, 2=A, 3=h, 4=P, 5=u_inf, 6=Q1 ...
%                     to Q at last element node

% n_n = number of nodes per element
% n_p = dimension of parametric space = 1
% n_q = number of quadrature points 

qp_unit = fopen('qp_store.bin','w'); % open file for post-processing
[r_q, s_q, t_q, w_q]=get_quadrature_rule (n_n, n_p, n_q); % get data

for j = 1:n_e    ; % loop over elements ====>> ====>> ====>> ====>> 
  e_nodes = nodes (j, 1:n_n)               ; % element connectivity 
  xy_e (1:n_n, 1) = x(e_nodes(1:n_n))       ; % x coord at el nodes 

  K   = el_prop (1) ; % . . .               % gather all properties
  Q_v = el_prop (6:5+n_n)                    ; % line load at nodes
  KA = K * A ; hP = h * P                    ; % combine properties

%      Zero all arrays to be integrated: S_e, M_g, c_inf

%      Numerical Integration of S_e, c_e , in 1D parametric
  for q = 1:n_q  ; % begin quadrature loop ---> ---> ---> ---> --->
    r = r_q (q) ; w = w_q (q)          ; % recover integration data

%       Element scalar interpolation, H and local derivatives DLH
    [H_q, DLH_q] = Lagrange_1D_library (n_n, r)     ; % interpolate

%       Interpolate global position of quadrature point
    xy_q = H_q * xy_e (1:n_n, :)           ; % interpolate global x

%       Compute the geometry mapping data, d_parm to d_physical
    Jacobian = DLH_q * xy_e (1:n_n, 1)        ; % local to physical
    J_det = det (Jacobian); J_Inv = inv (Jacobian); % det & inverse

%       Physical global derivatives of interpolations H are DGH
    DGH_q = J_Inv * DLH_q                     ; % derivatives dH/dx

%              Update element square matrices
    S_e = S_e + (DGH_q'*KA*DGH_q) * J_det * w_q (q)  ; % conduction

%         Generalized mass matrix update (and for variable Q)
    M_g = M_g + (H_q' * H_q) * J_det * w_q (q)       ; % convection
    c_inf = c_inf + hP * u_inf * J_det * w_q (q)   ; % u_inf source

%         save data needed for post-processing stresses
    if ( post == 1 ) % save post-processing data for this point
      fwrite (qp_unit, xy_q,  'double')        ; % save coordinates
      fwrite (qp_unit, KA,    'double')           ; % save material
      fwrite (qp_unit, DGH_q, 'double')           ; % save gradient
    end % if save integration point data for post-processing

  end ; % for q quadratures <--- <--- <--- <--- <--- <--- <--- <---

  c_e = c_inf + M_g * Q_v     ; % u_inf plus line source resultants
  S_e = S_e + hP * M_g     ; % conduction plus convection sq matrix

%   Insert completed element matrices into system matrices . . .
     % as before . . .
end % for each j element in mesh <<==== <<==== <<==== <<==== <<====