function Heat_Conduction_T3_XY (nbc, mixed)
%..................................................................
% Planar heat conduction with internal source, T3 triangle
%             XY COORDINATES CLOSED FORM INTEGRALS
%..................................................................
if ( nargin == 1 )     ;  % check for optional BC flags
  mixed = 0            ;  % default to no mixed (Robin) BC
elseif ( nargin == 0 ) ; 
  mixed = 0 ; nbc = 0  ;  % default to no mixed or Neumann BC
end % if from argument count
if ( mixed > 0 )       ;  % not done yet
  error('Sorry, mixed bc data not yet installed')
end % if lazy author

%            Set given constants 
n_g = 1                 ;  % number of DOF per node
n_n = 3                 ;  % number of nodes per element
n_i = n_n*n_g           ;  % number of DOF per element

% MODEL input file controls (for various data generators)
pre_e = 0 ; % Dummy items before el_type & connectivity
pre_p = 0 ; % Dummy items before BC_flag % coordinates

%       READ MESH AND EBC INPUT DATA
% specific problem data from MODEL data files (sequential)
load msh_bc_xyz.tmp       ; % bc_flag, x_coord, y_coord 
n_m = size (msh_bc_xyz,1) ; % number of nodal points in mesh
n_d = n_g*n_m             ; % system degrees of freedom (DOF)
if ( n_m == 0 )           ; % data missing !
  error ('Error missing file msh_bc_xyz.tmp')
end % if error
fprintf ('Read %g mesh coordinate pairs \n', n_m)
P = msh_bc_xyz (1:n_m, (pre_p+1)) ; % extract integer Packed BC flag
x = msh_bc_xyz (1:n_m, (pre_p+2)) ; % extract x column 
y = msh_bc_xyz (1:n_m, (pre_p+3)) ; % extract y column 
 
load  msh_typ_nodes.tmp        ; % el_type, connectivity list (3)
n_e = size (msh_typ_nodes,1)   ; % number of elements
if ( n_e == 0 )                ; % data file missing
  error ('Error missing file msh_typ_nodes.tmp')
end % if error
nod_per_el = size (msh_typ_nodes,2) - pre_e - 1 ; % data checking 
fprintf ('Read %g elements with %g nodes each \n', n_e, nod_per_el)
if ( nod_per_el ~= n_n )
  error (['Mesh requires %g nodes per element. \n', n_n])
  %b error ('Mesh requires 3 nodes per element')
end % if
el_type = msh_typ_nodes(:, pre_e+1) ; % element type number >= 1
n_t     = max(el_type)              ; % number of element types

%       BOOKKEEPING FOR BC FLAGS AND VALUES
%   Allocate for BC flags and values
EBC_flag  = zeros(n_m, n_g) ;     NBC_flag  = zeros(n_m, n_g) ; 
EBC_value = zeros(n_m, n_g) ;     NBC_value = zeros(n_m, n_g) ;

%   Extract EBC flags from packed integer flag P 
if ( n_g == 1 )     % following is not general enough
  EBC_flag (1:n_m, 1) = P (1:n_m)          ; % insert first flag
else
  error('Current EBC limited to 1 DOF per node')
end % if more than 1 DOF packed into EBC flag.
EBC_count = sum( P == 1)      ; % number of EBC values expected
 
%  Read EBC values
if ( EBC_count > 0 )          ; % read EBC data
  load msh_ebc.tmp            ; % node, DOF, value (eq. number)
  n_c = size(msh_ebc, 1)      ; % number of constraints
  if ( EBC_count ~= n_c )     ; % probable user error
    fprintf ('WARNING: mismatch in bc_flag count & msh_ebc.tmp')
  end % if user error
  for j = 1:n_c                    ; % loop over ebc inputs
    node  = msh_ebc (j, 1)         ; DOF = msh_ebc (j, 2) ;
    value = msh_ebc (j, 3)         ; % EBC value   
    Eq    = n_g * (node - 1) + DOF ; % row in system matrix
    EBC_value (Eq) = value         ; % should check EBC_flag
  end % for each EBC
end % if any EBC data expected
 
%  Read Neumann point data
if ( nbc > 0 )                ; % read NBC data
  load msh_nbc.tmp            ; % node, DOF, value (eq. number)
  n_u = size(msh_nbc, 1)      ; % number of point sources
  if ( n_u < 1 )              ; % missing data
    error ('No nbc data in msh_nbc.tmp')
  end % if user error
  for j = 1:n_u                    ; % non-zero Neumann pts
    node  = msh_nbc (j, 1)         ; DOF = msh_nbc (j, 2) ;
    value = msh_nbc (j, 3)         ; % NBC value
    Eq    = n_g * (node - 1) + DOF ; % row in system matrix
    NBC_flag  (Eq) = 1             ; % set NBC_flag
    NBC_value (Eq) = value         ; % insert value          
  end % for each EBC
end % if any NBC data expected

%  Reorder and count BC flags and values                        
if ( n_g > 1 )
  EBC_flag  = reshape (EBC_flag',  n_d, 1)    ; % change to vector 
  NBC_flag  = reshape (NBC_flag',  n_d, 1)    ; % change to vector 
  EBC_value = reshape (EBC_value', n_d, 1)    ; % change to vector 
  NBC_value = reshape (NBC_value', n_d, 1)    ; % change to vector 
end % if

% Optional storage allocation for reaction recovery
EBC_count = sum(EBC_flag)                   ; % number of EBC's
NBC_count = sum(NBC_flag)                   ; % number of NBC's
EBC_row   = zeros(EBC_count, n_d) ; EBC_col = zeros(EBC_count, 1) ;
EBC_react = zeros(EBC_count, 1)             ; % reaction allocation

%          SET ELEMENT PROPERTIES and internal source
% Manually set constant element properties (Fig 11.9 text)
t   = 1.0                               ; % default thickness
Kx = 8. ; Ky = 8. ; Kxy = 0.            ; % thermal conductivity
E_e = zeros (2, 2)                      ; % constitutive matrix
E_e (1, 1) = Kx   ; E_e (1, 2) = Kxy    ; % non-zero term
E_e (2, 1) = Kxy  ; E_e (2, 2) = Ky     ; % non-zero term
has_source = 1                          ; % internal source on
Q_e = 6                                 ; % component

%    GENERATE ELEMENT MATRICES AND ASSYMBLE INTO SYSTEM
%         Assemble n_d by n_d square matrix terms
S_e = zeros (n_i, n_i) ; C_e = zeros (n_i, 1) ; % initalize values
S   = zeros (n_d, n_d) ; C   = zeros (n_d, 1) ; % initalize sums

for j = 1:n_e                              ; % loop over elements
  thick   = t * el_type (j)                ; % integer multiple 
  e_nodes = msh_typ_nodes (j, (pre_e+2):(pre_e+4)); % connectivity
  
%  define nodal coordinates, ccw: i, j, k
  x_e = x(e_nodes) ; y_e = y(e_nodes)        ; % coord at el nodes
  x_i = x_e(1) ; x_j = x_e(2) ; x_k = x_e(3) ; % change notation
  y_i = y_e(1) ; y_j = y_e(2) ; y_k = y_e(3) ; % change notation

%  define centroid coordinates (quadrature point)
  x_cg = (x_i + x_j + x_k)/3 ;  y_cg = (y_i + y_j + y_k)/3 ;

%  geometric parameters: H_i (x,y) = (a_i + b_i*x + c_i*y)/two_a
  a_i = x_j * y_k - x_k * y_j ; b_i = y_j - y_k ; c_i = x_k - x_j ;
  a_j = x_k * y_i - x_i * y_k ; b_j = y_k - y_i ; c_j = x_i - x_k ;
  a_k = x_i * y_j - x_j * y_i ; b_k = y_i - y_j ; c_k = x_j - x_i ;

%  calculate twice element area
  two_a = a_i + a_j + a_k            ; % = b_j*c_k - b_k*c_j also

%  define 2 by 3 gradient matrix, B_e
  B_e (1, 1:3) = [ b_i,  b_j,  b_k ] / two_a ;
  B_e (2, 1:3) = [ c_i,  c_j,  c_k ] / two_a ;

%          ELEMENT MATRICES 
%  stiffness matrix, with constant jacobian 
  S_e = ( B_e' * E_e * B_e ) * thick * two_a / 2  ; % conduction

%  internal heat generation unit volume, Q_e
  if ( has_source ) % then form forcing vector            
    X_x = Q_e * thick * two_a / 6.0   ; %  resultant
    C_e (1:3) = [ X_x, X_x, X_x ]     ; % vector result 
  end % if or set up properties for body force

%         SCATTER TO SYSTEM ARRAYS
% Insert completed element matrices into system matrices
  rows (1:3) = n_g*(e_nodes(1:3) - 1) + 1 ; % element DOF numbers
  S (rows, rows) = S (rows, rows) + S_e   ; % add to system sq
  C (rows) = C (rows) + C_e               ; % add to sys column
end % for each j element in mesh

%         INSERT EBC AND NBC (if any)
kount = 0                                 ; % initialize counter
for j = 1:n_d       % check all DOF for displacement BC or pt force
  if ( NBC_flag(j) )                      ; % then NBC here
    C (j) = C (j) + NBC_value (j)         ; % add force value
  end % if NBC for this DOF

%         ENFORCE ESSENTIAL BOUNDARY CONDITIONS
  if ( EBC_flag(j) )                         % then EBC here
    if ( EBC_flag(j) == NBC_flag(j) )        % warning both here
      fprintf ('Usually an error: EBC and NBC at same DOF \n')
    end % if fatal BC data

%   Save reaction data to be destroyed by EBC solver trick
    kount = kount + 1                      ; % copy counter
    EBC_row(kount, 1:n_d) = S (j, 1:n_d)   ; % copy reaction data
    EBC_col(kount)        = C (j)          ; % copy reaction data

%   Insert EBC identity to avoid matrix partition accounting
    EBC = EBC_value(j)                ; % recover EBC value
    C (:) = C (:) - EBC * S (:, j)    ; % carry known column to RHS
    S (:, j) = 0 ; S (j, :) = 0       ; % clear, restore symmetry
    S (j, j) = 1 ; C (j) = EBC        ; % insert identity into row
  end % if EBC for this DOF
end % for over all j-th DOF

%             DISPLACEMENT SOLUTION & SAVE
T = S \ C                                ; % Compute displacements
fprintf ('\n')                             ; % blank
fprintf('Temperature at %g nodes \n', n_m) ; % header
if ( n_g > 1 )
  disp (reshape(T, n_g, n_m)')               % to look pretty
end % if
% save results (displacements) to MODEL file: node_results.tmp
fid = fopen('node_results.tmp', 'w')       ; % open for writing
for j = 1:n_g:n_d                          ; % save displacements
  fprintf ( fid, '%g \n', T (j))           ; % print
  fprintf (' %g %g \n', j, T (j))          ; % sequential save
end % for j DOF

%                REACTION RECOVERY & SAVE
if ( EBC_count > 0 )                    ; % reactions exist ?
  EBC_react = EBC_row * T - EBC_col     ; % matrix reactions (+-)
  % save reactions (forces) to MODEL file: node_reaction.tmp
  fprintf ('\n')                        ; % Skip a line
  fprintf ('Node, DOF, Value, Eq. for %g Heat Flow Reactions \n', ...
            EBC_count)                  ; % header
  fid = fopen('node_reaction.tmp', 'w') ; % open for writing
  kount = 0                             ; % initialize counter
  for j = 1:n_d                         ; % extract all EBC reactions
    if ( EBC_flag(j) )                      ; % then EBC here
  %   Output node_number, component_number, value, equation_number 
      kount = kount + 1                     ; % copy counter
      node  = ceil(j/n_g)                   ; % node at DOF j
      j_g   = j - (node - 1)*n_g            ; % 1 <= j_g <= n_g
      React = EBC_react (kount, 1)          ; % reaction value
      fprintf ( fid, '%g %g %g %g \n', node, j_g, React, j);% save
      fprintf ('%g %g %g %g \n', node, j_g, React, j);     % print
    end % if EBC for this DOF
  end % for over all j-th DOF
end % if EBC exist

%                 HEAT FLUX RECOVERY & SAVE
fid = fopen('el_qp_xyz_fluxes.tmp', 'w')    ; % open for writing
fprintf ('\n')                              ; % blank line
fprintf('Elem, X_qp, Y_qp, HeatFlux_x, HeatFlux_y \n');% header
for j = 1:n_e                              ; % loop over elements
  thick   = t * el_type (j)                ; % integer multiple 
  e_nodes = msh_typ_nodes (j, (pre_e+2):(pre_e+4)); % connectivity
  
%  define nodal coordinates, ccw: i, j, k
  x_e = x(e_nodes) ; y_e = y(e_nodes)        ; % coord at el nodes
  x_i = x_e(1) ; x_j = x_e(2) ; x_k = x_e(3) ; % change notation
  y_i = y_e(1) ; y_j = y_e(2) ; y_k = y_e(3) ; % change notation

%  geometric parameters: H_i (x,y) = (a_i + b_i*x + c_i*y)/two_a
  a_i = x_j * y_k - x_k * y_j ; b_i = y_j - y_k ; c_i = x_k - x_j ;
  a_j = x_k * y_i - x_i * y_k ; b_j = y_k - y_i ; c_j = x_i - x_k ;
  a_k = x_i * y_j - x_j * y_i ; b_k = y_i - y_j ; c_k = x_j - x_i ;

%  calculate twice element area
  two_a = a_i + a_j + a_k            ; % = b_j*c_k - b_k*c_j also

%          RECOVER B MATRIX & ITS LOCATION & ELEMENT DOF
%  define centroid coordinates (quadrature point)
  x_cg = (x_i + x_j + x_k)/3 ;  y_cg = (y_i + y_j + y_k)/3 ;

%  define 2 by 3 Gradient matrix, B_e
  B_e (1, 1:3) = [ b_i,  b_j,  b_k ] / two_a ;
  B_e (2, 1:3) = [ c_i,  c_j,  c_k ] / two_a ;

%  get DOF numbers for this element, gather solution
  rows (1:3) = n_g*(e_nodes(1:3) - 1) + 1 ; % element DOF numbers
  T_e = T (rows)                          ; % gather element DOF

%  COMPUTE GRADIENT & HEAT FLUX, SAVE LOCATION AND VALUES
  Strain = B_e * T_e                        ; % gradient vector
  Stress = E_e * Strain                     ; % heat flux vector
  fprintf (fid, '%g %g %g %g \n',  x_cg, y_cg, Stress(1:2));% save
  fprintf ('%g %g %g %g %g \n', j, x_cg, y_cg, Stress(1:2));% prt
end % for each j element in mesh

% End finite element calculations. 
% See /home/mech517/public_html/Matlab_Plots for graphic options
% http://www.owlnet.rice.edu/~mech517/help_plot.html for help

% end of Heat_Conduction_T3_XY