function [] = Integrate_Bar_matrices (n_n)        % Revised 01/19/19 
% ------------------------------------------------------------------
%     Mech 417, Finite Element Analysis, Rice University 
%     Copyright J.E. Akin 2019. All rights reserved.
% ------------------------------------------------------------------
% Numerically integrate the stiffness and source matrices for the 1D
%   constant coefficient ODE  [E A u']' + G A = 0 for the matrix 
% system (S) u  =  c + c_NBC for a three-noded (n_n = 3), or two-noded
% (n_n = 2) elements:  1---2 -->r or 1---2---3 -->r
%
% C_e = GA    * int_Le  [H_T] dx,    dx = dx/dr dr = Le dr
%     = GA    * int_0_1 [H_T] (Le dr)
%     = GA Le * int_0_1 [H_T] (dr)
%     = GA Le * sum_1_q [Hq_T] w_q,        Hq = H(r_q)
%
% S_e = EA    * int_Le  [dH/dx_T dH/dx] dx
%     = EA    * int_0_1 [(1/Le * dH/dr_T) (1/Le * dH/dr)] (Le dr)
%     = EA/Le * int_0_1 [(dH/dr_T) (dH/dr)] (dr)
%      = EA/Le * sum_1_q [dH/drq_T dH/drq] w_q, dH/drq = dH/dr(r_q)
%
% volume = A * int_Le dx = A * int_0_1 (Le dr) = A Le * int_0_1 dr
%        = A * Le * sum_1_q w_q
% ===================================================================
                                 % n_n = number of nodes per element
n_g = 1                        ; % number of DOF per node 
n_i = n_g * n_n                ; % number of DOF per element
n_q = 2                        ; % number of quadrature points

% n_n = 2: exact for integran polynomial of degree 3 (on 0 to 1)
r_q (1) = 2.1132486540518711774543e-01 ;
r_q (2) = 7.8867513459481288225457e-01 ;
w_q (1) = 5.0000000000000000000000e-01 ;
w_q (2) = 5.0000000000000000000000e-01 ;

%   Pre-allocate and initialize system matrix and vector
c_e = zeros (n_i,   1)            ; % clear array el source
S_e = zeros (n_i, n_i)            ; % clear array el stiffness
volume = 0                        ; % element volume

% Numerical Intergation 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, & local derivatives 
 switch n_n ;  % select from available library of elements
  case {2}  ;  % two node line element 1----2, L2
    H_q = [(1 - r),  r] ; dHdr_q = [-1,    1] ;

  case {3}  ;  % three node line quadratic 1--2--3, L3
    H_q   = [(1 - 3*r + 2*r^2), (4*r - 4*r^2), (2*r^2 - r)] ;
    dHdr_q = [(-3 + 4 * r),      (4 - 8*r),     (4*r -1)] ;

  otherwise ; % missing data
    fprintf ('Unknown number of nodes/element = %i \n', n_n)
    error ('Allowed values: 2 or 3')
 end ;  % switch on number of nodes

%   Update element square matrix, from K(x) 
  S_e = S_e + (dHdr_q' * dHdr_q) * w_q (q); % conduct

%   Source vector
  c_e = c_e + H_q' * w_q (q)              ; % source

%   Update integral of the interpolations
  volume = volume + w_q (q)               ; % element volume

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

fprintf ('Element matrix S_e = E*A by Le* \n') ; disp(S_e)
fprintf ('Element matrix c_e = G*A*Le* \n') ; disp(c_e)
fprintf ('Element volume = A*Le* %10.4e \n', volume)

% running gives
% >> Integrate_Bar_matrices(2)
% Element matrix S_e = E*A by Le*
%      1    -1
%     -1     1
% Element matrix c_e = G*A*Le*
%     0.5000
%     0.5000
% Element volume = A*Le* 1.0000e+00
% 
% >> Integrate_Bar_matrices(3)
% Element matrix S_e = E*A by Le*
%     2.3333   -2.6667    0.3333
%    -2.6667    5.3333   -2.6667
%     0.3333   -2.6667    2.3333
% Element matrix c_e = G*A*Le*
%     0.1667
%     0.6667
%     0.1667
% Element volume = A*Le* 1.0000e+00