function []=Plane_Stress_T6_Freq_lib_sm (n_modes)   % 4/4/2019
%.............................................................
% Plane Stress Frequency Analysis (2.5D),Isoparametric T6
%  Find the smallest n_modes natural frequencies and modes
%                 by Gaussian quadratures
%.............................................................
% Properties are:
% 1 Elastic Modulus   (3e7 psi) 
% 2 Poisson_s ratio   (0 bar theory, 0.3 actual)
% 3 Mass Density      (0.2812 lb/cu in) mislabelled
% 4 Thickness         (1 inch)      

if ( nargin == 0 )                 ; % check for optional data
  n_modes = 3                     ; % frequencies to be output   
end                                 ; % if from argument count
show_v = 0        % <= n_modes, number of eigenvectors to list
echo_user_remarks    ; % optional echo of file msh_remarks.txt

%         Application and element dependent controls
n_g = 2                             ; % number of DOF per node 
n_r = 3                       ; % number of rows in B_e matrix
echo_L = 0            ; % turn off echo of element connections
echo_N = 0               ; % turn off echo of node coordinates

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

%      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 (echo_L) ; 
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) ; M = zeros (n_d, n_d)      ; % initalize

%       Option mesh plot controls
inc_e = 0        ; % 0 or increment in element numbers to plot
inc_n = 0           ; % 0 or increment in node numbers to plot
plot_input_2d_mesh (x, y, nodes, inc_e, inc_n);    % mesh plot

%       Load the element properties array
load msh_properties.txt        ; % one row per element or type
n_prop = size(msh_properties, 1)       ; % == 1 if homogeneous
fprintf ('\n(Echoing columns of file msh_properties.txt) \n');
if ( n_prop == 1 )                 ; % material is homogeneous 
  E     = msh_properties (1, 1)      ; % modulus of elasticity
  nu    = msh_properties (1, 2)            ; % Poisson's ratio       
  Rho   = msh_properties (1, 3)               ; % mass density
  Thick = msh_properties (1, 4)                  ; % thickness 

%                 ECHO PROPERTIES
  fprintf ('\nProperties for all elements \n')
  fprintf ('Elastic Modulus = %g \n', E)     
  fprintf ('Poisson_s ratio = %g \n', nu)    
  fprintf ('Mass Density    = %g \n', Rho)
  fprintf ('Thickness       = %g \n', Thick)
end                          ; % if homogeneous read only once
if ( n_prop > n_e )                    ; % property data error
  error ('ERROR: property sets exceeds number of elements')
end                                 ; % if property data error

%     Tabulate quadrature data for unit triangle
n_q = 7                  ; % exact for 5-th degree polynomials
r_q = [0.10128651, 0.47014206, 0.79742699, 0.47014206, ...
       0.10128651, 0.05971587, 0.33333333 ] ;
s_q = [0.10128651, 0.05971587, 0.10128651, 0.47014206, ...
       0.79742699, 0.47014206, 0.33333333 ] ;
w_q = [0.06296959, 0.06619708, 0.06296959, 0.06619708, ...
       0.06296959, 0.06619708, 0.11250000 ] ;

%    GENERATE ELEMENT MATRICES AND ASSYMBLE INTO SYSTEM
%  Assemble n_d by n_d square matrix terms from n_e by n_e

Volume = 0.                                     ; % initialize
for k = 1:n_e  ; % element loop  ===>> ===>> ===>> ===>> ===>>
  e_nodes = nodes (k, 1:n_n)                  ; % connectivity
  
%               SET ELEMENT PROPERTIES & GEOMETRY
  if ( n_prop > 1 )    ; % constant properties in each element
    E     = msh_properties (k, 1)    ; % modulus of elasticity
    nu    = msh_properties (k, 2)          ; % Poisson's ratio       
    Rho   = msh_properties (k, 3)             ; % mass density
    Thick = msh_properties (k, 4)                ; % thickness       

%                   ECHO PROPERTIES
    fprintf ('\n');fprintf ('Properties for all elements \n');
    fprintf ('Elastic Modulus = %g \n', E)     
    fprintf ('Poisson_s ratio = %g \n', nu)    
    fprintf ('Thickness       = %g \n', Thick)
    fprintf ('Mass Density    = %g \n', Rho)
    if ( Thick <= 0 ) 
      error ('\nError: zero thickness at element %g ', k)
    end                                   ;  % if invalid data
  end                           ; % if not a homogeneous solid

  [E_e] = E_plane_stress (E, nu)      ; %Insert the properties
  
%                 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)         ; % y coord at el nodes

%            Numerical Intergation of S_e, M_e 
  S_e = zeros (n_i, n_i)               ; % initalize stiffness
  B_q = zeros (n_r, n_i)     ; % initalize strain-displacement
  M_e = zeros (n_i, n_i)                  ; % initalize masses
  M_x = zeros (n_n, n_n)                        ; % work array

  for q = 1:n_q ; % begin integration 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)] ;
   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] = B_planar_elastic (DGH_q, n_n)      ; % plane stress

%   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
    M_x = M_x + (H_q' * Rho * H_q)*t_q*J_det*w_q (q)  ; % mass  
    Volume = Volume + t_q * J_det * w_q (q)  ; % domain volume
   end % for q quadrature points <--- <--- <--- <--- <--- <--- 
   M_e (1:2:11, 1:2:11) = M_x     ; % repeat mass for x-motion
   M_e (2:2:12, 2:2:12) = M_x     ; % repeat mass for y-motion

%         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
  M (rows, rows) = M (rows, rows) + M_e   ; % add to system sq
end  ;% for k elements <<=== <<=== <<==== <<==== <<==== <<====

fprintf('\n'); 
fprintf('Total volume = %g \n', Volume)       ; % check volume
mass = Rho * Volume; fprintf('Total mass = %g \n', mass);
fprintf('\n');
fprintf ('NOTE: %i zero displacements imposed \n', EBC_count);
[Fixed] = get_fixed_from_ebc_flags (n_d, EBC_flag) ; % EBC dof
[Free]  = get_free_from_fixed (n_d, Fixed)     ; % id Free dof
f_n = max(size(Free))    ; % maximum number of modes available

% Solve for smallest eigenalues and eigen vectors
fprintf ('\n'); 
fprintf('Begin eigen-solution for %i frequencies \n',n_modes); 
exact = 5.315/(pi()*16^2)*sqrt(E/(3*Rho))
[V, D] = eigs (S(Free, Free), M(Free, Free), n_modes, 'sm')  ;
V = real(V); F_sq = diag (real(D))        ; % real frequency^2 

%*************************************************************
% NOTE: D is not in increasing order. They must be sorted.
%       Order(1:n_modes) is a vector subscript that gives the 
%       lowest n_modes eigenvalues (frequences squared) as the
%       vector F_sq(Order(1:n_modes))
%*************************************************************
[F_sq, Order] = sort (F_sq); % reorder eigenvalues^2, lo to hi

%*************************************************************
% NOTE: Now the eigenvalues are in the correct order, but not 
%       the eigenvectors. The Order array is needed to extract
%       the columns with lowest n_modes eigenvectors desired !
%*************************************************************

fprintf ('\n')           ; % list the n_modes eigenvalues only
fprintf ('NOTE: Rigid body motions give zero frequencies \n');
kount = 0                       ; % number of true eigenvalues
for k = 1:n_modes; % list eigenvalues only --> --> --> --> -->
  kount = kount + 1           ; % number of output eigenvalues
  fprintf('Eigenvalue %i = %9.4e \n', k, F_sq(k)); % reference 
  Freq (k)  = sqrt(F_sq(k))         ; % frequency, radians/sec
  Hertz (k) = Freq (k) / (2*pi)      ; % frequency, cycles/sec
  fprintf ('Natural frequency, Hz (rad/sec)')        ; % begin
  fprintf ('%i %10.5e (%10.5e) \n', kount, Hertz(k), Freq(k));
end ; % for k frequencies for output   <-- <-- <-- <-- <-- <-- 

fid = fopen('mode_value_vector.txt', 'w'); % evector save file 
fprintf ('\n') ; kount = 0      ; % number of true eigenvalues
V_Full = zeros (n_d, 1)      ; % allocate full DOF work vector
for K = 1:n_modes; % list requested eigenvalues ===> ===> ===>
  kount  = kount + 1            ; % number of true eigenvalues
  actual = Order(K)                 ; % where V is this column 
  V_Free = V(:, actual)             ; % free dof for this mode
  V_Full (Free) = V_Free (:)       ; % expand to original size
  XY_form = reshape (V_Full, n_g, n_m)'  ; % x-y displacements
  biggest = max ( max ( abs (XY_form) ) )       ; % scale them
  XY_form = XY_form / biggest                       ; % scaled

%         Plot the mode shape and original shape
  mode_2d_disp_plot (x, y, nodes, XY_form, Hertz(K), K);% plot

%         List selected eigenvectors, if any
  if ( show_v > 0 ); fprintf ('\n');% list show_v eigenvectors
    fprintf ('Frequency (Hz) %i = %9.4e \n', kount, Hertz(K));
    fprintf ('Eigenvector (mode shape): \n') ; % vector header
    fprintf ('node, x-,  y_component \n')    ; % vector values
    for k = 1:n_m ; % loop over all nodes  --> --> --> --> -->
       fprintf ('%i  %9.4e  %9.4e \n', k, XY_form (k, 1), ...
                                          XY_form (k, 2)) ; 
    end; % for k nodes <-- <-- <-- <-- <-- <-- <-- <-- <-- <--
  end                          ; % if list show_v eigenvectors

%   Save XY displacement rectangular form for mode shape plots
  Hz_K = Hertz (K)                      ; % data for save file
  for k = 1:n_m ; % loop over all nodes    --> --> --> --> -->
    fprintf (fid,'%i  %g  %g  %g \n', ...
             kount, Hz_K, XY_form (k, 1), XY_form (k, 2)) ;
  end ; % for save node k values   <-- <-- <-- <-- <-- <-- <--

  if ( kount == n_modes )    ; % have finished requested range
    break                       ; % break out of this for loop
  end                ; % if finished requested number of modes
end ; % for K modes of intetest  <=== <=== <=== <=== <=== <=== 

fclose(fid)         ; % done with all requested plot data save
fprintf ('\nNOTE: saved file mode_value_vector.txt for \n'); 
fprintf ('  plotting via mode_2d_mesh_plot.m \n'); % plot file
% end Plane_Stress_T6_Freq_lib_sm %===========================