function [] = elastic_2d_T3_lib ()       ; % revised 03/08/19
%.............................................................
%   Planar linear triangle (T3) plane-stress analysis (2.5D)    
%        for point loads and body forces (acceleration)
%.............................................................
% Properties: % 1 = E, 2 = nu, 3 = Body_x, 4 = Body_y
%               5 = Coef_Expansion, 6 = temp_increase
%               7 = optional element thickness
% Requires addpath /clear/www/htdocs/mech517/Akin_FEA_Lib
% ------------------------------------------------------
%         Application and element dependent controls
n_g = 2  ; % number of DOF per node  , x nad y displacements
n_q = 1  ; % number of quadrature points required
n_r = 3  ; % number of rows in B_e matrix & E_e matrix

load_pt = 1  ; % turn on(1)/off(0) point source (can be zero)
%        Echo user description of this application
echo_user_remarks ()      ; % property meanings, BC info, etc.

%              Read mesh input data files
[n_m, n_s, P, x, y, z] = get_mesh_nodes (n_g) ;

%      Extract EBC flags from packed integer flag P 
[EBC_flag] = get_ebc_flags (n_g, n_m, P)      ; % unpack flags
EBC_count  = sum( sum ( EBC_flag > 0 ) )      ; % # of EBC
 
[n_e, n_n, n_t, el_type, nodes] = get_mesh_elements ()  ;

shrink_2d_mesh     (x, y, nodes)    ; % plot 2D mesh
plot_input_2d_mesh (x, y, nodes,3,2); % plot 2D mesh
plot_bc_flags      (x, y, nodes, P) ; % plot BC flags

n_d = n_g*n_m                ; % system degrees of freedom (DOF)
n_i = n_g*n_n                ; % number of DOF per element
S = zeros (n_d, n_d) ; c = zeros (n_d, 1)      ; % initalize sums

%                  Read EBC values, if any
if ( EBC_count > 0 )                            ; % need EBC data
  [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag) ; % read data
end                                  ; % if any EBC data expected

%   Read point sources & couples data, if any, and insert in C
%   (from file msh_load_pt.txt, one line per point source)
if ( load_pt > 0 )                      ; % need point sources data
  [c_pt] = get_point_sources (n_g, n_m) ; % get point sources
  c = c + c_pt                          ; % add point sources
end % if any point source expected

%              Load the element properties array
load msh_properties.txt         ; % one row per element or type
n_prop = size(msh_properties, 1); % == 1 if homogeneous
n_vals = size(msh_properties, 2); % number of values per property
fprintf ('Read %i properties with %i values \n', n_prop, n_vals);
if ( n_prop > n_e )             ; % user input error
  error ('Number of property sets exceeds number of elements')
end                             ; % if user error
fprintf ('\n'); fprintf ('Properties for all elements \n')
% Expect large magnitude differences, change format
format short e ;  disp(msh_properties) ; format short       % print
if ( n_prop == 1 )              ; % just gather properties once
  % 1 = E, 2 = nu, 3 = Body_x, 4 = Body_y  
  % 5 = Coef_Expansion, 6 = temp_increase  
  % 7 = optional element thickness  
  Body_x=0; Body_y=0; coeff_e=0; Del_T=0; Thick=1; % defaults
  if ( n_vals > 1 ) ; % just elastic
    E       = msh_properties (1, 1) ;        % elastic modulus  
    nu      = msh_properties (1, 2) ; end % if % Poisson ratio  
  if ( n_vals > 2 ) ; % include body force per unit area
    Body_x  = msh_properties (1, 3) ;          % body x
    Body_y  = msh_properties (1, 4) ; end % if % body y
  if ( n_vals > 4 ) ; % include thermal loads
    coeff_e = msh_properties (1, 5) ; % material
    Del_T   = msh_properties (1, 6) ; end % if % Poisson ratio  
  if ( n_vals > 6 ) ; % allow varying thickness
    Thick   = msh_properties (1, 7) ; % thickness         
  end ; % if different element thicknesses
  Body    = [Body_x ; Body_y]     ; % insert in vector  
  thermal = coeff_e * Del_T ; % initial normal thermal strain  
  strain_0 = thermal * [1  1   0]'  ; % thermal strain  
end ; % if

%     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 MATRICES AND ASSEMBLE INTO SYSTEM
%  Assemble n_d by n_d square matrix terms from n_e by n_e

for k = 1:n_e ; % element loop ===>> ===>> ===>> ===>> ===>> ===>>
  e_nodes = nodes (k, 1:n_n)                      ; % connectivity
  type    = el_type (k)                           ; % element type
  
%               Set element properties
  if ( n_prop > 1 )
    Body_x=0; Body_y=0; coeff_e=0; Del_T=0; Thick=1; % defaults
    if ( n_vals > 1 ) ; % just elastic
      E       = msh_properties (1, 1) ;        % elastic modulus  
      nu      = msh_properties (1, 2) ; end % if % Poisson ratio  
    if ( n_vals > 2 ) ; % include body force per unit area
      Body_x  = msh_properties (1, 3) ;          % body x
      Body_y  = msh_properties (1, 4) ; end % if % body y
    if ( n_vals > 4 ) ; % include thermal loads
      coeff_e = msh_properties (1, 5) ; % material
      Del_T   = msh_properties (1, 6) ; end % if % Poisson ratio  
    if ( n_vals > 6 ) ; % allow varying thickness
      Thick   = msh_properties (1, 7) ; % thickness         
    end ; % if different element thicknesses
    Body    = [Body_x ; Body_y]     ; % insert in vector  
    thermal = coeff_e * Del_T ; % initial normal thermal strain  
    strain_0 = thermal * [1  1   0]'  ; % thermal strain  
  end % ; if

  % fill the constitutive matrix
  [E_e] = E_plane_stress (E, nu) ; % 3x3 Hooke's Law  
  
%                 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)   ; % 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
  B_q = zeros (n_r, n_i)                        ; % initalize values

  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
    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                  ; % 3 x 3          % properties 
    t_q = Thick                ; % 1 x 1          % thickness
    
% Set the non-zero terms in the vector interpolation  
    for k = 1: n_g ; % replace some zeros with H_q values  
      N_q (k, k:n_g:n_i) = H_q (1:n_n) ;  
    end % for shifted rows of scalar interpolations  

%      Insert gradient terms into the differential operator
    [B_q] = B_planar_elastic (DGH_q, n_n);% 3 x 6 strain-displacement  

%                 ELEMENT MATRICES UPDATES, 6 x 6 and 6 x 1
    S_e = S_e + (B_q' * E_q * B_q) * t_q * J_det * w_q (q) ; % square
    c_e = c_e + N_q' * Body * t_q * J_det * w_q (q)        ; % vector 
    c_e = c_e + B_q'*E_e*strain_0*J_det*w_q (q)    ; % thermal source 
  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
  S (rows, rows) = S (rows, rows) + S_e        ; % add to system sq
  c (rows)       = c (rows)       + c_e        ; % add to sys column
end ; % for k-th element <<== <<=== <<=== <<==== <<==== <<==== <<====

%    ALLOCATE STORAGE FOR OPTIONAL REACTION RECOVERY          
if ( EBC_count > 0 )                        ; % reactions occur
  [EBC_row, EBC_col] = save_reaction_matrices (EBC_flag, S, c);
end              ; % if essential BC exist (almost always true)

%         ENFORCE ESSENTIAL BOUNDARY CONDITIONS
save_resultant_load_vectors (n_g, c)         ; % save before overwrite
[S, c] = enforce_essential_BC (EBC_flag, EBC_value, S, c); % overwrite

%             COMPUTE SOLUTION & SAVE
Ans = S \ c                             ; % compute displacements 
list_save_results (n_g, n_m, Ans)           ; % save for plotting 

%              PLOT DISPLACEMENTS
color_deformed_mesh     (x, y, nodes, Ans, 1); % with amplified by 1
quiver_2d_displacements (x, y, nodes, Ans, 1); % with amplified by 1
color_displacement      (x, y, nodes, Ans, 0); % magnitude only

%          OPTIONAL REACTION RECOVERY & SAVE
if ( EBC_count > 0 )                            ; % reactions exist ?
  [EBC_react] = recover_reactions_print_save (n_g, n_d, ...
                EBC_flag, EBC_row, EBC_col, Ans) ; % reaction to EBC 
end ; % if EBC exist

%       POST-PROCESS ELEMENTS (element heat flows, save to plot)
fprintf ('\nPost-prossessing stresses: \n')
fprintf ('Output stress order is S_xx, S_yy, S_xy, S_vm, Tau_max \n')
fid = fopen('el_qp_xyz_fluxes.txt', 'w'); % open for writing for plots

strain = zeros (n_r, 1) ; % pre-allocate strains
stress = zeros (n_r, 1) ; % pre-allocate stresses
max_max_VM = 0          ; % max of all von Mises stress
max_max_MS = 0          ; % max of all max shear stress

%  Are there initial thermal strains present ? 
if ( n_prop > 5 ) % then initial thermal strain input 
  thermal = el_prop (5) * el_prop (6) ; % alpha * temp change 
end % if thermal load 

for k = 1:n_e ; % loop over all elements ====>> ====>> ====>> ====>>
  e_nodes = nodes (k, 1:n_n)                      ; % connectivity 
  type    = el_type (k)                          ; % material type
  [rows] = get_element_index (n_g, n_n, e_nodes)  ; % eq numbers 
  Ans_e (1:n_i) = Ans(rows)     ; % gather element displacements 
   
%               Recover needed property for heat flux 
  if ( n_prop > 1 )                     ; % update for this element
    Body_x=0; Body_y=0; coeff_e=0; Del_T=0; Thick=1; % defaults
    if ( n_vals > 1 ) ; % just elastic
      E       = msh_properties (type, 1) ;        % elastic modulus  
      nu      = msh_properties (type, 2) ; end % if % Poisson ratio  
    if ( n_vals > 2 ) ; % include body force per unit area
      Body_x  = msh_properties (type, 3) ;          % body x
      Body_y  = msh_properties (type, 4) ; end % if % body y
    if ( n_vals > 4 ) ; % include thermal loads
      coeff_e = msh_properties (type, 5) ; % material
      Del_T   = msh_properties (type, 6) ; end % if % Poisson ratio  
    if ( n_vals > 6 ) ; % allow varying thickness
      Thick   = msh_properties (type, 7) ; % thickness         
    end ; % if different element thicknesses
    Body    = [Body_x ; Body_y]     ; % insert in vector  
    thermal = coeff_e * Del_T ; % initial normal thermal strain  
    strain_0 = thermal * [1  1   0]'  ; % thermal strain  
  end                     ; % if different property in each element
  [E_e] = E_plane_stress (E, nu) ; % 3x3 Hooke's Law  
  
%            Gather nodal coordinates (for plotting)
  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

%   Physical location and physical gradient at intergration points
  for q = 1:n_q   ; % begin loop  ---> ---> ---> ---> ---> ---> --->
    r = r_q (q) ; s = s_q (q) ; w = w_q (q)     ; % recover data
%                 element type interpolation functions
    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 (physical) derivatives of H
    DGH_q = J_Inv * DLH_q      ; % 2 x 3          % dH/dx & dH/dy
    E_q   = E_e                ; % 3 x 3          % properties 
    [B_q]=B_planar_elastic (DGH_q, n_n);% 3x6 strain-displacement 

%  Compute the generalized strains & stresses 
   strain = B_q * Ans_e'               ;  % in-plane strains 
   stress = E_e * (strain - strain_0)  ;  % in-plane stresses 

%  compute the material failure estimates 
   vonMise = sqrt( (stress(1)-stress(2))^2 + stress(2)^2 ... 
                  + stress(1)^2 + 6 * stress(3)^2 ) * 0.7071     ; 
   maxTau  = sqrt( 0.25 * (stress(1)-stress(2))^2 + stress(3)^2 ); 
   
  fprintf('El, Q_Pt, Strains  %i, %i, %9.3e  %9.3e  %9.3e \n', ... 
           k, q, strain (1:n_r)) % list 
  fprintf('El, Q_Pt, Stresses %i, %i, %9.3e  %9.3e  %9.3e  %9,3e  %9.3e \n', ... 
           k, q, stress (1:n_r), vonMise, maxTau) % list 

% Sequentially save coord, stresses, failure criteria for plots
  fprintf(fid,'%9.3e  %9.3e  %9.3e  %9.3e  %9.3e %9.3e  %9.3e \n', ...
               xy_q (1:n_s), stress (1:n_r), vonMise, maxTau) ;

    if ( vonMise > max_max_VM ) % update worst value
      max_max_VM = vonMise ; % current value
      max_VM_el  = k       ; % current element
      max_VM_pt  = q       ; % current quadrature point
    end % if current max VM value
    if ( maxTau > max_max_MS ) % update worst value
      max_max_MS = maxTau  ; % current value
      max_MS_el  = k       ; % current element
      max_MS_pt  = q       ; % current quadrature point
    end % if current max MS value

  end ; % for quadratures <--- <--- <--- <--- <--- <--- <--- <--- <--- 
end ; % for each k element in mesh  <<==== <<==== <<==== <<==== <<==== 
fprintf ('\n')
fprintf ('==================================================\n')
fprintf ('Maximum von Mises stress = %9.3e \n', max_max_VM)
fprintf ('at point %i of element %i \n', max_VM_pt, max_VM_el)
fprintf ('Maximum shear stress     = %9.3e \n', max_max_MS)
fprintf ('at point %i of element %i \n', max_MS_pt, max_MS_el)
fprintf ('==================================================\n')
fclose (fid) ; % close saved stresses
fprintf ('Saved stress results in el_qp_xyz_fluxes.txt for plots \n')

% Plot various stress components
addpath /clear/www/htdocs/mech517/Matlab_Plots
quiver_qp_2D_stress (1)  ; % Maximum principal stress
quiver_qp_2D_stress (2)  ; % Minimum principal stress

color_qp_stresses_step (5) ; % Maximum shear stress
color_qp_stresses_step (4) ; % vom Mises equivalent
color_qp_stresses_step (3) ; % XY shear stress
color_qp_stresses_step (2) ; % YY normal stress
color_qp_stresses_step (1) ; % XX normal stress
% end of elastic_2d_T3_lib      

function [E_e] = E_plane_stress (E, nu, G) ; % ========
% * * * * * * * * * * * * * * * * * * * * * * * * * * *
%     PLANE STRESS CONSTITUTIVE MATRIX DEFINITION
%         STRESS COMPONENT ORDER: XX, YY, XY
% * * * * * * * * * * * * * * * * * * * * * * * * * * *
% E   = Young's modulus of elasticity                         
% G   = shear modulus (not independent)
% nu  = Poisson's ratio                                       
% E_e = constitutive matrix 
% n_r = number of rows in B and E matrices = 3

if ( nargin == 2 )                ; % classic definition
  G  = 0.5 * E / (1 + nu)                  ; % in theory
end                              ; % if special material
c_1 = E / (1 - nu^2) ; c_2 = c_1 * nu      ; % constants

E_e = zeros (3, 3)                         ; % largely 0
E_e(1, 1) = c_1 ; E_e(2, 1) = c_2 ;          ; % normals
E_e(1, 2) = c_2 ; E_e(2, 2) = c_1 ;          ; % normals
E_e(3, 3) = G                                  ; % shear
% end E_plane_stress ===================================

function [B] = B_planar_elastic (DGH, n_n) ; %=========
% * * * * * * * * * * * * * * * * * * * * * * * * * * *
%      ELASTICITY STRAIN-DISPLACEMENT RELATIONS  (B)
%  For plane-stress or plane-strain: n_g = 2, n_r = 3
% * * * * * * * * * * * * * * * * * * * * * * * * * * *
% B   = strain-displacement matrix   (n_r, n_n * n_g)
% DGH = global derivatives of H      (n_s, n_n * n_g)
% n_g = number of displacements per node
% n_n = number of nodes per element
%       Stress and strain order is: xx, yy, xy

n_g = 2 ; n_r = 3    ; % for planar stress
n_i = n_n * n_g      ; % number of element dof
B = zeros (n_r, n_i) ; % initialize strain-displacement

for J = 1:n_n ; % loop over nodes ==>> ==>> ===>> ===>>
  K = n_g * (J - 1) + 1     ; % u-column contribution
  L = K + 1                 ; % v-column contribution
  B (1, K) = DGH (1, J)     ; % du/dx for normal xx 
  B (3, K) = DGH (2, J)     ; % du/dy for shear  xy 
  B (2, L) = DGH (2, J)     ; % dv/dy for normal yy 
  B (3, L) = DGH (1, J)     ; % dv/dx for shear  xy
end ; % for J nodes of element  <<=== <<=== <<=== <<===
% end B_planar_elastic ================================