%  ** GENERATE ELEMENT MATRICES AND ASSYMBLE INTO SYSTEM **
%  Assemble system n_d by n_d square matrix terms from the
%         element n_e by n_e terms
 
for j = 1:n_e  ; % loop over elements ====>> ====>> ====>> ====>> 
  e_nodes = nodes (j, 1:n_n)                      ; % connectivity 
  type    = el_type (j)            ; % current element type number
   
%                 Gather nodal coordinates
  xy_e (1:n_n, 1) = x(e_nodes)             ; % x coord at el nodes 
  xy_e (1:n_n, 2) = y(e_nodes)   ; % 6 x 2   % y coord at el nodes

%     Allocate and clear element arrays
  H_q   = zeros (n_n, 1)      ; % clear element interpolations
  DLH_q = zeros (n_s, n_n)    ; % interpolation local derivatives 
  Jacobian = zeros (n_s, n_s) ; % geometry Jacobian
  DGH_q = zeros (n_s, n_n)    ; % interpolation global derivatives 
  B_q   = zeros (n_r, n_i)    ; % initalize B matrix values
  S_e   = zeros (n_i, n_i)    ; % clear array el stiffness
  M_e   = zeros (n_i, n_i)    ; % clear array el mass
  C_e   = zeros (n_i, 1)      ; % clear array el sources

%               SET ELEMENT MATERAILS & GEOMETRY 
  if ( n_mats > 1 )    % each element has some different properties
    E      = msh_properties (j, 1)   ; % modulus of elasticity
    nu     = msh_properties (j, 2)   ; % Poisson's ratio       
    Body_x = msh_properties (j, 3)   ; % x-body force component
    Body_y = msh_properties (j, 4)   ; % y-body force component
    Thick  = msh_properties (j, 5)   ; % thickness       
    Rho    = msh_properties (j, 6)   ; % mass density

%      populate the stress-strain law matrix
    E_e = E/(1 - nu^2) * [ 1,   nu,    0       ; ...
                           nu,  1,     0       ; ...
                           0,   0,  (1 - nu)/2 ] ; % stress-strain law

%                   ECHO PROPERTIES
    fprintf ('Element number %i \n', j)           
    fprintf ('Elastic Modulus             = %g \n', E)     
    fprintf ('Poisson;s ratio             = %g \n', nu)    
    fprintf ('X-body force component      = %g \n', Body_x)
    fprintf ('Y-body force component      = %g \n', Body_y)
    fprintf ('Thickness                   = %g \n', Thick)
    fprintf ('Mass Density                = %g \n', Rho)
    if ( Thick <= 0 ) 
      error ('\nError: zero thickness at element %g ', j)
    end % if invalid data
  end % if 
 
%   ELEMENT, FOUNDATION STIFFNESSES AND SOURCE RESULTANT MATRICES 

%            Numerical Intergation of S_e, C_e 
  for q = 1:n_q   ; % begin quadrature loop ---> ---> ---> ---> ---> --->
    r = r_q (q) ; s = s_q (q) ; w = w_q (q)     ; % recover data

%                 element type interpolation functions, 1 x 6
    H_q = [ (1 - 3*r - 3*s + 2*r*r + 4*r*s + 2*s*s), (2*r*r - r), ...
            (2*s*s - s), 4*(r - r*r - r*s), 4*r*s, 4*(s - r*s - s*s) ] ;

%                 element type parametric derivatives, 2 x 6
    DLH_q = [ (-3 + 4 * r + 4 * s), (-1 + 4 * r), 0, ...
              (4 - 8 * r - 4 * s), (4 * s), (-4 * s) ; ... % dH/dr
              (-3 + 4 * r + 4 * s), 0, (-1 + 4 * s), ...
              (-4 * r), (4 * r), (4 - 4 * r - 8 * s) ] ;   % dH/ds

%             Global position and Jacobian
    xy_q     = H_q * xy_e           ; % 1 x 2     % global (x, y)
    Jacobian = DLH_q * xy_e         ; % 2 x 2     % d(x,y)/d(r,s)
    J_det    = abs( det (Jacobian)) ; % 1 x 1     % determinant
    J_Inv    = inv (Jacobian)       ; % 2 x 2     % inverse

%             Global derivatives of H
    DGH_q = J_Inv * DLH_q      ; % 2 x 6          % dH/dx & dH/dy
    E_q = E_e                  ; % 3 x 3          % properties
    t_q = Thick                ; % 1 x 1          % thickness
    
%      Insert gradient terms into the differential operator
    B_q (1, 1:2:n_i) = DGH_q (1, 1:n_n)   ; % for esplion_x
    B_q (2, 2:2:n_i) = DGH_q (2, 1:n_n)   ; % for esplion_y
    B_q (3, 1:2:n_i) = DGH_q (2, 1:n_n)   ; % for gamma     
    B_q (3, 2:2:n_i) = DGH_q (1, 1:n_n)   ; % for gamma     

%                 ELEMENT MATRICES UPDATES, 12 x 12 and 12 x 1
    S_e = S_e + (B_q' * E_q * B_q) * t_q * J_det * w_q (q) ; % square
 %b M_e = M_e + (H_q' * Rho * H_q) * t_q * J_det * w_q (q) ; % mass  
    if ( any ([Body_x, Body_y]) ) % then need loads
      C_e(1:2:n_i) = C_e(1:2:n_i) + H_q' * Body_x *t_q*J_det*w_q (q); %Fx
      C_e(2:2:n_i) = C_e(2:2:n_i) + H_q' * Body_y *t_q*J_det*w_q (q); %Fy
    end % if body forces active

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

%         SCATTER TO (ASSEMBLE INTO) SYSTEM ARRAYS 
% Insert completed element matrices into system matrices 
  [rows] = get_element_index (n_g, n_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 
%b M (rows, rows) = M (rows, rows) + M_e  ; % add to system mass 
%       ***  END SCATTER TO SYSTEM ARRAYS ***
 
end % for each j element in mesh   <<==== <<==== <<==== <<==== <<==== 
%       *** END ELEMENT GENERATION