function Torsion_2D_T3_XY  ;
%.................................................................. 
%          Torsion of non-circular bars, T3 triangle 
%             XY COORDINATES CLOSED FORM INTEGRALS 
%.................................................................. 
%            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]) 
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) ;     EBC_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 
  fprintf ('Read %g essential boundary conditions \n', n_c) 
  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)         ; % system node number
    DOF   = msh_ebc (j, 2)         ; % DOF == 1 for now
    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 
  
%  Reorder and count BC flags and values                         
if ( n_g > 1 ) 
  EBC_flag  = reshape (EBC_flag',  n_d, 1)    ; % change to vector  
  EBC_value = reshape (EBC_value', n_d, 1)    ; % change to vector  
end % if 
 
% Optional storage allocation for reaction recovery 
EBC_count = sum(EBC_flag)                   ; % number of EBC'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  
G = 8.e6                                ; % material shear modulus 
E_e = eye (2, 2)/ G                     ; % constitutive matrix 
has_source = 1                          ; % internal source on 
beta = 1.745e-4                         ; % twist per unit length 
fprintf ('Shear modulus = %g, Twist per unit length = %g \n', ...
                          G,                          beta)
 
%    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 
    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 
  if  two_a <= 0 % then probable connectivity error
    fprintf ('WARNING: zero or negative area at \n ')
    fprintf ('Element %g , Nodes %g %g %g \n', j, e_nodes)
    two_a = abs (two_a)
  end % if negative area 
 
%  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  
  S_e = ( B_e' * E_e * B_e ) * two_a / 2  ; % stiffness  
 
%  global twist loading 
  if ( has_source ) % then form forcing vector             
    X_x = 2 * beta * two_a / 6.0     ; 
    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 
% fprintf (' S '); disp(S)
% fprintf (' C '); disp(C')
 
%    OPTIONAL SAVE OF REACTION RECOVERY DATA, AND DIAGONAL SCALE
kount = 0                                 ; % initialize counter 
S_ave = 0                                 ; % initialize average
for j = 1:n_d       % check all DOF for displacement BC or pt force 
  S_ave = S_ave + S (j, j)                ; % update average
  if ( EBC_flag (j) )                       % then EBC here 
%   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 
  end % if EBC for this DOF 
end % for over all j-th DOF 
S_ave = S_ave / n_d                        ; % true disgonal average
% fprintf (' EBC_row '); disp(EBC_row)
% fprintf (' EBC_col '); disp(EBC_col')
% fprintf (' EBC_flag '); disp(EBC_flag')

%         INSERT EBC BY MODIFYING S & C
kount = 0                                 ; % initialize counter 
for j = 1:n_d    % check all DOF for displacement BC or pt force 
%         ENFORCE ESSENTIAL BOUNDARY CONDITIONS 
  if ( EBC_flag (j) )                         % then EBC here 
%   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) = S_ave                  ; % set diagonal identity     
    C (j) = EBC * S_ave               ; % set row identity          
  end % if EBC for this DOF 
end % for over all j-th DOF 
% fprintf (' S '); disp(S)
% fprintf (' C '); disp(C')
 
%             STRESS FUNCTION SOLUTION & SAVE 
T = S \ C                                  ; % Compute solution 
fprintf ('\n')                             ; % blank 
fprintf('Stress function value 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 (stress function) to MODEL file: node_results.tmp 
fid = fopen('node_results.tmp', 'w')       ; % open for writing 
for j = 1:n_g:n_d                          ; % each stress function
  fprintf ( fid, '%g \n', T (j))           ; % save  
  fprintf (' %g %g \n', j, T (j))          ; % print
end % for j DOF 
 
%     NODAL REACTION (SHEAR STRESS) RECOVERY & SAVE 
if ( EBC_count > 0 )                    ; % reactions exist ? 
  EBC_react = EBC_row * T - EBC_col     ; % matrix reactions (+-) 
  % save reactions (stresses) to MODEL file: node_reaction.tmp 
  fprintf ('\n')                        ; % Skip a line 
  fprintf ('Node, DOF, Shear Force Reaction \n', EBC_count) 
  fid = fopen('node_reaction.tmp', 'w') ; % open for writing 
  kount = 0                             ; % initialize counter 
  for j = 1:n_d                         ; % extract EBC reactions 
    if ( EBC_flag(j) )                      ; % then EBC here 
  %   Output node_number, component_number, value
      kount = kount + 1                     ; % copy counter 
      node  = ceil(j/n_g)                   ; % node at DOF j 
      DOF   = j - (node - 1)*n_g            ; % 1 <= DOF <= n_g 
      React = EBC_react (kount, 1)          ; % reaction value 
      React = React * G                     ; % now shear force
      fprintf ( fid, '%g %g %g %g \n', node, DOF, React, j);% save 
      fprintf ('%g %g %g %g \n', node, DOF, React)        ; % print 
    end % if EBC for this DOF 
  end % for over all j-th DOF 
end % if EBC exist 
 
%    ELEMENT SHEAR STRESS RECOVERY & SAVE.  TORQUE RECOVERY 
fid = fopen('el_qp_xyz_fluxes.tmp', 'w')    ; % open for writing 
fprintf ('\n')                              ; % blank line, header
fprintf('Elem, X_qp,  Y_qp,  Tau_xz,  Tau_yz,  Tau,  Torque_part \n')
torque = 0                                 ; % initialize 
for j = 1:n_e                              ; % loop over elements 
    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 
  two_a = abs (two_a)                ;
 
%          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 ; 
 
%   define integral of shape functions 
  C_e (1:3) = two_a / 6 ;                       
 
%  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 & SHEAR STRESS, SAVE LOCATION AND VALUES 
  Strain = B_e * T_e                ; % gradient vector 
  Tau_xz =  Strain(2)               ; % shear stress component
  Tau_yz = -Strain(1)               ; % shear stress component
  Tau = sqrt(Tau_xz^2 + Tau_yz^2)   ; % total shear stress
 
%  ELEMENT CONTRIBUTION TO TORQUE 
  torque_e = 2 * C_e' * T_e                 ; % element torque 
  torque = torque + torque_e                ; % running total 
  fprintf (fid, '%g    %g    %g      %g  \n', ...
                 x_cg, y_cg, Tau_xz, Tau_yz); 
  fprintf ('%g %g    %g    %g      %g      %g   %g \n', ...
            j, x_cg, y_cg, Tau_xz, Tau_yz, Tau, torque_e); 
end % for each j element in mesh 
fprintf ('Total applied torque = %g (times symmetry level) \n', torque) 
 
% See /home/mech517/public_html/Matlab_Plots for graphic options 
% http://www.owlnet.rice.edu/~mech517/help_plot.html for help 
% end of Torsion_2D_T3_XY