%     Tabulate quadrature data for unit triangle
n_q = 1                         ; % exact for linear polynomials
r_q (1) = 1./3 ; s_q (1) = 1./3 ; w_q (1)= 1./2       ; % tabulated

%                 Generate element arrays
for j = 1:n_e     ; % element loop ===>> ===>> ===>> ===>> ===>> ===>>
  thick   = t * el_type (j)                ; % integer multiple 
  E_e     = [ kx_e (j),  0 ; ...
              0,   ky_e(j) ] ; % orthotropic thermal conductivities
  e_nodes = nodes (j, 1:n_n)    ; % 1 x 3    % connectivity
  
%                 Define 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)   ; % 3 x 2   % y coord at el nodes

%            Numerical Intergation of S_e, C_e 
  S_e = zeros (n_i, n_i) ; C_e = zeros (n_i, 1) ; % initalize values
  E_q = zeros (2, 2)     ; B_q = zeros (2, n_i) ; % initalize values

  for q = 1:n_q   ; % begin loop ---> ---> ---> ---> ---> ---> ---> --->
    r = r_q (q) ; s = s_q (q) ; w = w_q (q)     ; % recover data
    H_q = [(1-r-s), r, s]          ; % 1 x 3      % interpolations
    DLH_q = [  -1, 1, 0 ; ...                     % dH/dr
               -1, 0, 1 ]          ; % 2 x 3    ; % 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 3          % dH/dx & dH/dy
    E_q = E_e                  ; % 2 x 2          % properties
    t_q = thick                ; % 1 x 1          % thickness
    
%      Insert gradient terms into the differential operator
    B_q = DGH_q     ; % 2 x 3, just an identity for conduction

%                 ELEMENT MATRICES UPDATES, 3 x 3 and 3 x 1
    S_e = S_e + (B_q' * E_q * B_q) * t_q * J_det * w_q (q) ; % square
    C_e = C_e + Q_e * t_q * J_det * w_q (q)                ; % vector
   end % for quadratures   <--- <--- <--- <--- <--- <--- <--- <--- <--- 

%                  SCATTER TO SYSTEM ARRAYS
  rows = e_nodes   ; % equation calculation is an identity here, n_g=1
  S (rows, rows) = S (rows, rows) + S_e     ; % add to system sq
  C (rows)       = C (rows) + C_e           ; % add to sys column
end % for elements     <<==== <<==== <<==== <<==== <<==== <<==== <<====