function Lubrication_2D_Lib (DEFAULT) % Revised 4/17/19  
% DEFAULT omitted, require user interactive inputs
% DEFAULT = 1 use quadratic triangle (T6) planar bearing     
% --------------------------------------------------------
% Hydrodynamic Lubrication: Incompressible Reynold's Eq
% Del . ((h^3/12mu) Del P) = Del . {h U_ave + (h^3/12mu) B_f}
%                          + v_d + h_dot

% where h(x, y) is the thin film z-thickness, mu is the viscosity,
% B_f(x, y) is the 2D tangential body force vector B_x i + B_y j,
% U_ave(x, y) is the 2D average velocity vector of the top and 
% bottom of the thickness U_ave = 0.5{(U_0 + U_h)i + (V_0 + V_h)j},
% v_d is the z-velocity into a porous bearing (usually zero),
% h_dot is the z-squeeze film velocity (rate of change of h),
% and P(x, y) is the bearing pressure in the z-direction.
% Del is the gradient operator vector ( ),x i + ( ),y j.
% (The script also accepts one-dimensional models.)

% For 2D the Galerkin integral form is:
% Area_integral {[c*dp_dx^2 + c*dp_dy^2] + h*[U_ave . Del p]
% + p*[v_d + h_dot] + c*[B_f . Del p]} 
% + Edge_integral {Q p} = 0, for c == h^3/(12 mu)

%   The film thickness is either input with the nodal coordinates
% (default), or input as a list of values for every element.
% Iff the thickness is input at the nodes, then the theoretical
% pressure gradient would be continuous, and it would be valid
% to use C_1 element interpolation for a more accurate solution.

%   The velocity vector, U_ave, can be for bearing rotating with
% an angular velocity.  That condition must be manually set with
% control variable ang_vel, and the film thickness can then be
% calculated by setting controls tilt_ang and tilt_ref.

%   The solution gives the pressure at every node, the pressure 
% gradient at every quadrature point, the resultant normal force
% (integral of P), its location, and the moments about the global
% axes (integral of y*P and x*P, respectively.

%   Recall that the number of nodes to define the element geometry,
% n_x, is less than or equal to the number of analysis nodes per
% element, n_n. The current code is isoparametric, n_x == n_n.
% The element properties can be 3 or 5 or 7 or 9 or (9+n_x) items
% on a single line.

%   The fluid viscosity, mu, is taken as constant at the element 
% centroid.
%   The top velocity components U_h and V_h are usually constant
% and are input as constant at the element centroid. (3 items)
% The above three items are required input.  

%   The squeeze film velocity h_dot is usually zero, but can be
% input as constant at the element centroid. 
%   The porous velocity, v_d, is usually zero, but can be input
% as constant at the element centroid. (5 items)

%   The body force vector is usually zero, but B_x and B_y can be 
% input as constant at the element centroid. (7 items)
%   The base velocity components U_0 and V_0 are usually zero,
% but can be input as constant at the element centroid. (9 items)

%   If no other element properties are input, then the film thickness 
% is input at each node as a fake z-coordinate.  Flag node_h signals 
% that condition, and it is turned on if 9 or less items are entered.
% It is also turned on if 3 "coordinates" are input.

%   Since h is often discontinuous between elements, it can be input
% at the nodes of each element as element data.  If  node_h == 1, 
% then the element thickness values are ADDED to the continuous
% thickness input at the nodes (to cause discontinuous element films).
% Otherwise, they are the actual element thicknesses that can be 
% discontinuous between elements.

%   There are a maximum of (n_x + 9) required property inputs per element:
% 1=mu,	2=U_h, 3=V_h, 4=h_dot, 5=v_d, 6=B_x, 7=B_y, 8=U_0, 9=V_0,
% and 10=h_1, 11=h_2, 12=h_3  ... (9+n_x)=h_n_x thicknesses

%   The essential boundary condition, EBC, is that the pressure, 
% P, is given on portions of the boundary (exterior or interior).
% The natural boundary condition is that the normal gradient of
% the pressure vanishes on boundaries without other conditions.
% The nonessential boundary condition specifies the net normal 
% flow per unit length on a portion of the boundary as:
% (h^3/12mu) Del P . n = Q - h U_ave . n - (h^3/12mu) B . n
% where n is the unit normal vector on that boundary portion
% and Q is a known flow per unit length, positive in.
%............................................................. 
%     Mech 417, Finite Element Analysis, Rice University 
%     Copyright J.E. Akin 2013. All rights reserved.
%............................................................. 
% NOTE: input text files begin with msh_ and end with .txt
% see file FEA_notation.m for variable descriptions
% see file FEA_functions.m for functions summary
global axisym  ; % the axisymmetric mesh flag, global variable
 
% echo information about input and output text files
% echo_files_info   % uncomment column 1, if desired

% echo user information about this specific application
echo_user_remarks

% =====================================================================
%                             PHASE 1
%           *** SET OR READ USER PROBLEM DEPENDENT CONTROLS *** 
% =====================================================================
if ( nargin == 0 ) % default to user defined field controls
  DEFAULT = 0 ; % require user inputs
end % if use default application

%  set application controls
if ( DEFAULT == 0 ) % then read user control inputs
  % axisym  = 1 if the problem is axisymmetric
  % load_pt = 1 if point sources are input
  % n_g = number of DOF per node
  % N_N = number of EXPECTED nodes per element
  % n_p = dimension of parametric space
  % n_q = number of quadrature points required 
  % n_r = number of rows in element B_e matrix 
  % N_S = number of EXPECTED physical space dimensions
  [axisym, load_pt, n_g, N_N, n_p, n_q, n_r, N_S] = get_user_controls ;
  if ( axisym > 0 ) % incompatible data
    fprintf ('WARNING: mesh is not revolved. Setting axisym=0 \n')
    axisym = 0
  end % if invalid controls
  n_x = input ('Enter number of geometric nodes per element, n_x \n') ;
  if ( n_x > N_N )
    fprintf ('WARNING: n_x exceeds analysis nodes. Setting n_x = n_n \n')
    n_x = N_N
  end % if invalid controls
  node_h = input ('Is film thickness input at the nodes? 1=yes, 0=no \n');

elseif ( DEFAULT == 1 ) % 2D T6 plane conductivity

%  set user maximum array sizes 
  n_g     = 1     ; % number of DOF per node 
  N_N     = 6     ; % number of EXPECTED nodes per element
  n_p     = 2     ; % dimension of parametric space
  n_q     = 4     ; % number of quadrature points
  n_r     = 2     ; % number of rows in element B_e matrix
  N_S     = 2     ; % EXPECTED physical space dimensions
  n_x     = N_N   ; % number of geometry (thickness) nodes per element
  axisym  = 0     ; % turn on(1)/off(0) the axisymmetric flag
  load_pt = 0     ; % turn on(1)/off(0) point source (load can be zero)
  node_h  = 1     ; % turn on(1)/off(0) continuous node film thicknesses
end % if DEFAULT

%   set user logic choices (or add read function)
color    = 0  ; % turn on(1)/off(0) a color plot of the results
echo_p   = 1  ; % turn on(1)/off(0) listing of properties
el_react = 0  ; % turn on(1)/off(0) listing of element reactions
integral = 1  ; % turn on(1)/off(0) listing of solution integrals
post     = 1  ; % turn on(1)/off(0) post-processing option
show_bc  = 0  ; % turn on(1)/off(0) list assembled S and C, after EBC
show_e   = 0  ; % turn on(>0)/off(0) list some S_e and C_e
show_qp  = 1  ; % turn on(1)/off(0) list qp post-processing results
show_s   = 0  ; % turn on(1)/off(0) listing of assembled S and C
vector   = 0  ; % turn on(1)/off(0) the pressure gradient vector plot

% set manual data for tilt geometry or angular velocity
ang_vel   = 0  ; % turn on (not 0) angular velocity about origin (rps)
tilt_ang  = 6.66667e-06 % 0  ; % turn on (not 0) tilt angle about tilt line (degrees)
tilt_h    = 1.e-5 % 0  ; % turn on (>0) tilt height along tilt_ref
tilt_list = 1  % ; % turn on(1)/off(0) list of tilt node h values 
tilt_ref  = 0  % ; % angle to radial line where tilt_ang acts (degrees)

% Is post-processing required ?
if ( integral > 0 | show_qp > 0 | vector > 0 ) 
  post = 1 ; % then post-processing is required
end % if

%    Echo current user logic flags
fprintf ('User logic flags: turn on(1)/off(0) \n') 
fprintf ('color    plot the result in color         = %i \n', color   )
fprintf ('echo_p   listing of properties            = %i \n', echo_p  )
fprintf ('integral listing of solution integral     = %i \n', integral)
fprintf ('load_pt  point source (can be zero)       = %i \n', load_pt )
fprintf ('node_h   turn on(1)/off(0) node thickness = %i \n', node_h  )
fprintf ('post     post-processing option           = %i \n', post    )
fprintf ('show_bc  debug list S & C, after EBC      = %i \n', show_bc )
fprintf ('show_e   list some S_e and C_e            = %i \n', show_e  )
fprintf ('show_qp  list qp post-process results     = %i \n', show_qp )
fprintf ('show_s   list assembled S and C           = %i \n', show_s  )
fprintf ('vector   plot flux vector at some qp      = %i \n', vector  )
fprintf ('\n')

%    Echo all default or input controls
fprintf ('User control variables \n')
fprintf ('n_g number of DOF per node                = %i \n', n_g)
fprintf ('n_n number of analysis nodes per element  = %i \n', N_N)
fprintf ('n_x number of geometric nodes per element = %i \n', N_N)
fprintf ('n_p dimension of parametric space         = %i \n', n_p)
fprintf ('n_q number of quadrature points           = %i \n', n_q)
fprintf ('n_r number of rows in B_e                 = %i \n', n_r)
fprintf ('n_s dimension of physical space           = %i \n', N_S)

% Continue checking controls for this application
C1_flag = 0                    ; % default to a 2D C0 element type
if ( n_g > 1 )                   % was expecting default C_0 element
 if ( node_h == 1 )    % then continuous theoretical pressure gradient
  fprintf ('NOTE: Elements with C_1 continuity are allowed (n_g=2) \n')
  fprintf ('      since the film thickness is given at the nodes. \n')
 else                        % film thickness probably not continuous
  fprintf ('WARNING: Only elements with C_0 continuity (n_g=0) are \n')
  fprintf ('         allowed if film thickness is not continuous. \n')
 end % if C0 or C1 line element
 C1_flag = 1                ; % set a flag for C1 interpolation choice
end % if C1 elements used

if ( tilt_ang ~= 0 | ang_vel ~= 0 ) % to calculate velocity or thickness
  fprintf ('\n User U_ave velocity and/or thickness control variables \n')
  fprintf ('ang_vel   angular velocity (not 0) about origin = %g \n', ang_vel )
  fprintf ('tilt_ang  tilt angle (not 0) about tilt_ref     = %g \n', tilt_ang)
  fprintf ('tilt_h    tilt height (>0) along tilt_ref line  = %g \n', tilt_h  )
  fprintf ('tilt_list turn on(1) list tilt node h values    = %i \n', tilt_list)
  fprintf ('tilt_ref  angle, from x, where tilt_ang acts    = %g \n', tilt_ref)
  if ( ang_vel ~= 0 )
    fprintf ('NOTE: velocity will be calculated at quadrature points \n')
    fprintf ('      from ang_vel (rps) and radial positions. \n')
  end % if 
  if ( tilt_ang ~= 0 )
    fprintf ('NOTE: film thickness will be calculated at quadrature \n')
    fprintf ('      points from tilt_ang and tilt_ref (degrees) \n')
    fprintf ('      and the height along tilt_ref, tilt_h. \n')
    fprintf ('\n')
    if ( tilt_list == 1 ) % list calculated node thickness for checking
      %B list_tilt_node_h (tilt_ang, tilt_h, tilt_ref, x, y) %B TBA
    end % if              
    %B fprintf ('ERROR manual film thickness input required \n')
    %B error('The above features have not yet been validated')
  end % if tilt angle given
end % if tilted or rotating bearing

% =====================================================================
%                               PHASE 2
%              *** BEGIN DATA INPUT AND CHECKING PHASE ***
% =====================================================================
fprintf ('Begin Phase 2, input mesh, properties, sources and EBC \n')

%    Read mesh nodal input data files 
%    (from msh_bc_xyz.txt, with one line per node)

% n_m  = number of nodal points in mesh
% n_s  = number of space dimensions
% PD   = nodal packed digit essential BC flag (n_g digits 0 or 1)
%      NOTE: PD is either 0 or 1 here
% x    = mesh x-coordinate values
% y    = mesh y-coordinate values, if present, etc. for z
[n_m, n_s, PD, x, y, z] = get_mesh_nodes (n_g); % bc_flags, coords

% Has the film thickness been input at each node?
if ( n_s == 1 )                                               % 1D
  if ( node_h > 0 ) ;                     % Expect thickness too !
    fprintf ('ERROR: node_h set but no node thicknesses input \n')
    error ('Read on x-coordinates without thicknesses')
  end % if node thickness error
end % if thickness input as coordinate
if ( n_s == 2 )                                        % may be 1D
  if ( node_h > 0 )                          % y is film thickness
    fprintf ('NOTE: second coordinate is the film thickness \n')
    if ( any (y <= 0) ) Invalid thickness
      fprintf ('ERROR: encountered zero or negative thickness \n')
      error ('Read zero or negative nodal film thickness')
    end % if bad data
    n_s = 1                         ; % reset true space dimension
  end % if h input
end % if 1D or 2D
if ( n_s == 3 )                   % must be 2D with film thickness
  fprintf ('NOTE: third coordinate is the film thickness \n')
    if ( any (z <= 0) ) Invalid thickness
      fprintf ('ERROR: encountered zero or negative thickness \n')
      error ('Read zero or negative nodal film thickness')
    end % if
  n_s = 2 ; node_h = 1 ;          % reset space, turn on node flag
end % if thickness input as z-coordinate

%   Compare actual data against space control
if ( n_s ~= N_S ) % then possible control error
  fprintf ('WARNING: Read %i coordinates, but expected %i \n', n_s, N_S)
end % if give warning
 
% *** Read mesh element types and nodal connections (padded with 0's) ***
% (from msh_typ_nodes.txt, with one line per element)

% n_e     = number of elements in the mesh
% n_n     = maximum number of nodes per element
% el_type = type number of every element, usually 1
% nodes   = node connection list for each element (pad with 0 to size n_n)
% n_t     = number of different element types
[n_e, n_n, n_t, el_type, nodes] = get_mesh_elements  ; 

%   Compare actual data against connection control
if ( n_n ~= N_N ) % number of expected nodes
  fprintf ('WARNING Read %i nodes per element but expected %i \n', n_n, N_N)
end % if inconsistent inputs

% Allocate element type control data if more than one type input
type_i = zeros (n_t, 1) ; % maximum DOF per el type
type_n = zeros (n_t, 1) ; % maximum nodes per el type
type_p = zeros (n_t, 1) ; % dimension of type parametric space
type_q = zeros (n_t, 1) ; % number of type quadrature points
type_x = zeros (n_t, 1) ; % number of type geometry nodes (default type_n)

% Default to a single element type
if ( n_t == 1 ) % then initialize
  type_n = n_n ; type_p = n_p ; type_q = n_q ; type_x = n_n ;

else % read controls for each type

% *** Read each element type control integers ***
% (from msh_ctrl_types.txt, with one line per element type)
% type_n = maximum nodes per el type 
% type_p = dimension of type parametric space
% type_q = number of type quadrature points
% type_x = number of type geometry nodes (default type_n)
  [type_n, type_p, type_q, type_x] = get_el_type_controls (n_t) ;
%   check type control vs system control
  has_n = max ( type_n ) ; % actual max node per el
  has_p = max ( type_p ) ; % actual max parm dim
  has_q = max ( type_q ) ; % actual max Gauss rule
  has_x = max ( type_x ) ; % actual max geom nodes
  if ( has_n > n_n | has_x > n_n ) % warn
    fprintf ('Warning: type_n = %i is inconsistent with mesh controls \n', has_n)
    n_n = has_n % increase the maximum number of nodes per element
  end % if give warning
  if ( has_p > n_p ) % warn
    fprintf ('Warning: type_p = %i is inconsistent with mesh controls \n', has_p)
    n_p = has_p % increase the maximum parametric space dimension
  end % if give warning
  if ( has_q > n_q ) % warn
    fprintf ('Warning: type_q = %i is inconsistent with mesh controls \n', has_q)
    n_q = has_q % increase the maximum number of quadrature points
  end % if give warning
end % if more than one element type
type_i = n_g * type_n ; % number of degrees of freedom per type

%    Read any material properties for the elements
%    (from msh_properties.txt, one line per element or type)

% n_mats = number of materials, = n_e, or n_t if homogeneous
% n_vals = number of property values for each element or type
% msh_properties (n_mats, n_vals) = system properties array
[n_mats, n_vals, msh_properties] = get_mesh_properties (echo_p, n_e, n_t);
fprintf ('NOTE: read %i materials \n', n_mats)
if ( n_vals == 3 | n_vals == 5 | n_vals == 7 | n_vals == 9 )
  fprintf ('      with %i properties, as expected \n', n_vals)
  fprintf ('      and film thicknesses input at the nodes. \n')
  if ( node_h ~= 1 ) % should have been set
    fprintf ('      Setting node_h == 1 \n')
    node_h == 1
  end % if give warning
elseif ( n_vals == 9+n_x ) % discontinuous films
  fprintf ('      with 9 properties and %i \n', n_x)
  fprintf ('      element film thicknesses or increments. \n')
else ; % user error
  fprintf ('ERROR: wrong number of properties (3,5,7,9 or %i ) \n', 9+n_x)
  error ('Invalid element property count')
end % if property count

fprintf ('Above properties are: \n')
fprintf ('TYPE 1: Bearing pressure area element for \n')
fprintf ('Del.((h^3/12mu) Del P) = Del.{h U_ave + (h^3/12mu) B_f} \n')
fprintf ('+ v_d + h_dot with U_ave = 0.5{(U_0 + U_h)i + (V_0 + V_h)j} \n')

fprintf ('Input order: 1=mu, 2=U_h, 3=V_h, 4=h_dot, 5=v_d, \n')
fprintf ('             6=B_x, 7=B_y, 8=U_0, 9=V_0 and optionally \n')
fprintf ('             %i element films or film increments \n', n_x)

fprintf ('TYPE 2: Edge element for normal flow with \n')
if ( n_vals == 3 )                       % then expected count
  fprintf ('        property 1=Q_in, and two zeros \n')
elseif ( n_vals == 5 )                   % then expected count
  fprintf ('        property 1=Q_in, and four zeros \n')
elseif ( n_vals == 7 )                   % then expected count
  fprintf ('        property 1=Q_in, and six zeros \n')
elseif ( n_vals == 9 )                   % then expected count
  fprintf ('        property 1=Q_in, and eight zeros \n')
else                                           % actual count
  fprintf ('        property 1=Q_in, and %i zeros \n', n_vals-1)
end % if Type 2 properties

%   Count maximum local and system equation numbers 
n_d = n_g * n_m                ; % system degrees of freedom (DOF)
n_i = n_g * n_n                ; % number of DOF per element 

%   Pre-allocate and initialize system matrices and vectors.
% Ans       = solution vector
% C         = total source vector
% C_pt      = point source part of C for reaction recovery
% EBC_flag  = n_g individual binary flags at each node 
% EBC_value = EBC value for each restrained DOF
% el_prop   = element or type properties
% S         = stiffness matrix       
% Type_integral   = integral of solution in each element
% Type_integral_x = integral of x * solution in each element
% Type_integral_y = integral of y * solution in each element

Ans       = zeros (n_d, 1)   ; % results
C         = zeros (n_d, 1)   ; % force or source
C_pt      = zeros (n_d, 1)   ; % point source portion of C
C_react   = zeros (n_d, 1)   ; % copy of C for reaction recovery
EBC_flag  = zeros (n_m, n_g) ; % initialize all DOF flags
EBC_value = zeros (n_m, n_g) ; % default EBCs to  zero
el_prop   = zeros (n_vals, 1); % initialize element or type properties
h_x       = zeros (n_x, 1)   ; % variable nodal gap thicknesses
Q_v       = zeros (n_n,   1) ; % variable source nodal values
S         = zeros (n_d, n_d) ; % generalized stiffness
Type_integral   = zeros (n_t, 1) ; % integral of the solution
Type_integral_x = zeros (n_t, 1) ; % integral of the x * solution
Type_integral_y = zeros (n_t, 1) ; % integral of the y * solution
Type_volume     = zeros (n_t, 1) ; % volume of each type of element

%   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

%    Extract and count EBC flags from packed integer flag PD
[EBC_flag] = get_ebc_flags (n_g, n_m, PD)     ; % unpack flags

% EBC_count = number of essential (Dirichlet) BC's 
EBC_count  = sum( sum ( EBC_flag == 1 ) )    ; % number of EBC  
if ( EBC_count == 0 )
  fprintf ('WARNING: No essential boundary condition given. \n')
else
  fprintf ('Note: expecting %i essential BC values. \n', EBC_count)
end % if give warning

%   Read essential boundary condition (EBC) values, if any
%   (from file msh_ebc.txt one line per constraint)
if ( EBC_count > 0 )                            ; % need EBC data
  [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag) ; % read data
else                                      % unusual, but possible
  fprintf ('WARNING: no essential boundary conditions flagged \n')
end                                     % if any EBC data expected

%   Prepare a default general 2D mesh plot (uncomment as desired)
%   using /clear/www/htdocs/mech517/Matlab_Plots/mesh_shrink_plot.m
% mesh_shrink_plot % show a space between elements
% pause (3) % wait 3 seconds before next plot
%   or using /clear/www/htdocs/mech517/Matlab_Plots/mesh_plot.m
% mesh_plot    % no space between elements, mesh_plot (inc_e, inc_p)

%B %  Set an integer to simplify element type selection below
%B domain   = n_p * 100 + n_g * 10 + n_n ; % parametric reference code
%B domain_x = n_p * 100 + 10 + n_x       ; % geometric code, n_g==1

%    *** END OF DATA INPUT AND CHECKING PHASE ***

% ======================================================================
%                              PHASE 3
%       ** GENERATE ELEMENT VALUES AND ASSYMBLE INTO SYSTEM **
% ======================================================================
fprintf ('Begin Phase 3, assemble (scatter) element arrays \n')
if ( show_e > 0 ) % then list some element matrices
  fprintf ('Displaying the first %i sets of element matrices \n', show_e)
end % if arrays to display
 
%     Open unit needed for post-processing operations
if ( post == 1 ) % open post-processing element unit
  el_unit = fopen ('el_store.bin', 'w')      ; % open for element writing
  if ( el_unit == -1 )                                      % open failed
    fprintf ('\nWARNING: post-processing unit open failed \n')
    post = 0                          ; % turn off post-processing option
  end % if unit available
  qp_unit = fopen ('qp_store.bin', 'w')          ; % open for  GQ writing
  if ( qp_unit == -1 )                                      % open failed
    fprintf ('\nWARNING: post-processing unit open failed \n')
    post = 0   % turn off post-processing option
  end % if unit available
end % if open post-processing unit

%    *** Initially fill the type==1 quadrature data files ***
% r_q, s_q, t_q = parametric coordinates of quadrature point
% w_q           = parametric weight      of quadrature point
if ( n_q > 0 )                 % numerical integration to be used
  r_q = zeros (n_q, 1) ; s_q = zeros (n_q, 1) ;    % allocate
  t_q = zeros (n_q, 1) ; w_q = zeros (n_q, 1) ;    % allocate
  [r_q, s_q, t_q, w_q] = get_quadrature_rule (n_n, n_p, n_q) ;
end % if numerical integration data required 
 
type_old = 1                                                ; % initialize
h_dot = 0 ; v_d = 0 ; B_x = 0 ; B_y = 0 ; U0 = 0 ; V0 = 0   ; % initialize
for j = 1:n_e  ; % loop over elements ====>> ====>> ====>> ====>> 
  type    = el_type (j)        ; % current element type number
  typ_i   = type_i (type)      ; % maximum DOF per el type
  typ_n   = type_n (type)      ; % maximum nodes per el type
  typ_p   = type_p (type)      ; % dimension of type parametric space
  typ_q   = type_q (type)      ; % number of type quadrature points
  domain  = typ_p * 100 + n_g * 10 + typ_n ; % reference code
  domain_x  = typ_p * 100 + 10 + n_x ; % geometry code, n_g==1
  e_nodes = nodes (j, 1:typ_n) ; % connectivity 
   
%                 Gather nodal coordinates
  xy_e (1:typ_n, 1) = x(e_nodes(1:typ_n))     ; % x coord at el nodes 
  if ( n_s > 1 ) % then 2D or 3D
    xy_e (1:typ_n, 2) = y(e_nodes(1:typ_n))   ; % y coord at el nodes
  end % if 2D or 1D

%     Allocate and clear element type arrays
  C_e   = zeros (typ_i,   1)      ; % clear array el sources
  C_ex  = zeros (typ_i,   1)      ; % clear array el sources
  C_ey  = zeros (typ_i,   1)      ; % clear array el sources
  C_r   = zeros (typ_i,   1)      ; % clear array el reactions
  E_e   = zeros (n_r,   n_r)      ; % material constitutive matrix
  H_i   = zeros (typ_i,   1)      ; % scalar interpolations integral %B
  H_ix  = zeros (typ_i,   1)      ; % x * interpolations integral %B
  H_iy  = zeros (typ_i,   1)      ; % y * interpolations integral %B
  Jacobian = zeros (n_s, n_s)     ; % geometry Jacobian (one-to-one)
  S_e   = zeros (typ_i, typ_i)    ; % clear array el stiffness

%     Set constant element type properties
  if ( n_mats > n_t )    % each element has some different properties
    J = j ; % the element number
  else % use type data
    J = type ; % the type number
  end % if which property values
  el_prop (1:n_vals) = msh_properties (J, 1:n_vals) ;

%     Shorthand type property names
  switch ( type ) % branch on element type

   case {1} % continuum element for PDE
%    Del . ((h^3/12mu) Del P) = Del . {h U_ave + (h^3/12mu) B_f}
%    + v_d + h_dot with U_ave = 0.5{(U_0 + U_h)i + (V_0 + V_h)j}
%    P          = unknown              ; % pressure

     if ( node_h == 1 )                    % then thickness at nodes
       if ( n_s == 1 )                     % 1D film
         h_x = y (e_nodes (1:n_x))       ; % copy from fake y-coordinate
       else                                % 2D film
         h_x = z (e_nodes (1:n_x))       ; % copy from fake z-coordinate
       end % if get thickness from fake coordinate
     end % if thickness given as coordinate
     mu         = el_prop (1)            ; % viscosity
     Uh         = el_prop (2)            ; % top x-velocity
     Vh         = el_prop (3)            ; % top v-velocity
     if ( n_vals > 3 )
       h_dot = el_prop (4)               ; % squeeze rate
       v_d   = el_prop (5)               ; % diffuse velocity
     end % if add squeeze data
     if ( n_vals > 5 )
       B_x = el_prop (6)                 ; % x body force per unit area
       B_y = el_prop (7)                 ; % y body force per unit area
     end % if add body forces
     if ( n_vals > 7 )
       U0 = el_prop (8)                  ; % base x-velocity
       V0 = el_prop (9)                  ; % base y-velocity
     end % if moving base
     if ( n_vals > 9 )                     % film may be discontonuous
       if ( node_h == 0 )                  % use element values only
         h_x(1:n_x) = el_prop (10:9+n_x) ; % element film or increment
       else                                % add element increment
         h_x = h_x + el_prop (10:9+n_x)  ; % node value plus increment
       end % if insert discontinuity
     end % if film thickness is discontinuous
%B disp(h_x)
     mu12 = 12 * mu ;

     case {2} % edge element with known source
%      TYPE 2 edge elements with knowmn flux per unit length
%      With properties 1=q_in, and others zero 
       q_in  = el_prop (1) ; % nornal flux

     otherwise
       error ('Type number of element exceeds contol value \n')
  end % switch on type number

%     Allocate and clear quadrature point arrays
  B_q   = zeros (n_r,   typ_i)  ; % initalize B matrix values, at q
  DGH_q = zeros (n_s,   typ_n)  ; % interpolation global derivatives, at q 
  DLH_q = zeros (typ_p, typ_n)  ; % interpolation local derivatives, at q
  Jacobian = zeros (typ_p, n_s) ; % geometry Jacobian
  G_q   = zeros (1,     n_x)    ; % geometry interpolation, at q
  H_q   = zeros (1,     typ_n)  ; % clear scalar interpolations, at q
  xy_q  = zeros (1,       n_s)  ; % physical location of point q

%   ELEMENT, FOUNDATION STIFFNESSES AND SOURCE RESULTANT MATRICES 

%   Get quadrature rule for current element type
  if ( type ~= type_old ) % then 
    [r_q, s_q, t_q, w_q] = get_quadrature_rule (typ_n, typ_p, typ_q) ;
  end  % if 

%   Average velocity vector for this element (unless rotating)
  AveUx = 0.5 * (U0 + Uh)                       ; % x-direction
  AveVy = 0.5 * (V0 + Vh)                       ; % y-direction

%   Numerical Integation of S_e, C_e , in 2D parametric
  volume = 0        ; % element volume
  for q = 1:typ_q   ; % begin quadrature loop ---> ---> ---> ---> ---> --->
    r = r_q (q) ; w = w_q (q) ; % recover data
    s = s_q (q) ; t = 0       ; % t not used in 2D

%   Geometry interpolation, G
    [G_q, DLG_q] = el_shape_n_local_deriv (domain_x, r, s, t) ;

%   Element scalar interpolation, H, & local derivatives, DLH
    if ( C1_flag == 1 )                     % then C1 1D element
      L = abs(xy_e (n_n, 1) - xy_e (1, 1)); % physical length
      [H_q, DGH_q, D2GH, D3GH] = el_shape_n_global_deriv_C1 ...
                                 (n_g, n_n, n_p, r, L) ;
    else                                    % then default C0 element
      [H_q, DLH_q] = el_shape_n_local_deriv (domain, r, s, t) ;
    end % if C0 or C1 element

%   Interpolate global position of quadrature point
    %B xy_q     = H_q * xy_e (1:typ_n, :) ; %   1 x n_s   % global (x, y)
    xy_q     = G_q * xy_e (1:typ_n, :) ; %   1 x n_s   % global (x, y)

%   Interpolate velocity if a rotating pad bearing
    if ( ang_vel ~= 0 ) % then rotation about origin (0, 0)
      AveUx = -ang_vel * xy_q (2)                    ; % x-direction
      AveUy =  ang_vel * xy_q (1)                    ; % y-direction
    end % if insert rotational velocity at point

%   Interpolate the gap thickness at point q
    h_q = G_q * h_x                                    ; % from n_x nodes
    if ( h_q <= 0 ) % not possible
      fprintf ('ERROR: element %i, point %i thickness is %g \n', ...
                j, q, h_q)
      error ('Invalid interpolated film thickness')
    end % if negative thickness
    h_3 = h_q^3                                        ; % at point q

%   Compute the geometry mapping data, d_parm to d_physical
    [J_det, Jacobian] = parametric_to_physical (n_s, typ_n, typ_p, ...
                                                DLG_q, xy_e);
              
    if ( type == 1 ) % then continuum element

%     Global derivatives of H are DGH, for a C0 element
      if ( C1_flag == 0 )             ; % usually not a C1 line element
        J_Inv = inv (Jacobian)          ; % n_s x n_s   % inverse if C0
        DGH_q = J_Inv * DLH_q           ; % n_s x n_n   % dH/dx & dH/dy
      end                                 % if a C1 line element

%     Insert gradient terms into the differential operator
      [B_q] = DGH_q                   ; % same for scalars

%     Update element square matrices
      S_e = S_e + h_3/mu12 * (B_q' * B_q) * J_det * w_q (q) ; % support matrix

%     Velocity sources %B are C0 fix for C1
      C_ex = C_ex + h_q * DGH_q(1, :)' * AveUx * J_det * w_q (q) ; % x-velocity
      if ( n_s > 1 ) % then 2D
        C_ey = C_ey + h_q * DGH_q(2, :)' * AveVy * J_det * w_q (q) ; % y-velocity
      end % if second velocity component

%     Update integral of the interpolation functions, for post-process
      H_i    = H_i    + H_q' * J_det * w_q (q)             ; % H integral
      H_ix   = H_ix   + H_q' * J_det * w_q (q) * xy_q(1)   ; % x*H integral
      if ( n_s > 1 )                                         % then 2D
        H_iy   = H_iy   + H_q' * J_det * w_q (q) * xy_q(2) ; % y*H integral
      end % if 2D problem
      volume = volume + J_det * w_q (q)        ; % element volume

%     save data needed for post-processing stresses
      if ( post == 1 ) % save data for this point
        fwrite (qp_unit, xy_q,  'double') ; % save coordinates
        fwrite (qp_unit, E_e,   'double') ; % save material   %B
        fwrite (qp_unit, B_q,   'double') ; % save strain-displacement
      end % if save point data

    elseif ( type == 2 ) % then a known normal flux
 %     -K_n * U,n + g_in = 0
       C_e = C_e +  H_q' * q_in * J_det * w_q (q) ;

    else
      error ('Element arrays not defined for type > 3 \n')
    end % if element type
   end % for quadratures   <--- <--- <--- <--- <--- <--- <--- <--- <--- 
   C_e = C_ex + C_ey                        ; % combine sources if 2D

   Type_volume (type) = Type_volume (type) + volume ;

   if ( show_e > 0 ) % then list some element matrices
     if ( j <= show_e ) % show the first few
       fprintf ('Element matrix S_e = \n') ; disp(S_e)
       fprintf ('Element matrix C_e = \n') ; disp(C_e)
     end % if elements wanted
   end % if show_e

%  Save element arrays needed for post-processing stresses
   if ( post == 1 ) % save arrays for this element
     fwrite (el_unit, xy_e,    'double') ; % save coordinates
     fwrite (el_unit, el_prop, 'double') ; % save properties
     fwrite (el_unit, S_e,     'double') ; % save beam stiffness
     fwrite (el_unit, C_e,     'double') ; % save line load resultant
     fwrite (el_unit, H_i,     'double') ; % save interpolation integral
     fwrite (el_unit, H_ix,    'double') ; % save x * interpolation integral
     fwrite (el_unit, H_iy,    'double') ; % save y * interpolation integral
   end % if save element arrays

%         SCATTER TO (ASSEMBLE INTO) SYSTEM ARRAYS 
% Insert completed element matrices into system matrices 
  [rows] = get_element_index (n_g, typ_n, e_nodes); % eq numbers 
  % rows = vector subscript converting element to system eq numbers
  S (rows, rows) = S (rows, rows) + S_e  ; % add to system stiffness 
  C (rows)       = C (rows)       + C_e  ; % add to sys force/couples 

%       ***  END SCATTER TO SYSTEM ARRAYS ***

end % for each j element in mesh   <<==== <<==== <<==== <<==== <<==== 

if ( show_s == 1 ) % then list each system matrix
  fprintf ('System matrix S = \n') ; disp(S)
  fprintf ('System matrix C = \n') ; disp(C)
end % if show_s

fprintf ('Total input source terms = %g \n', sum (C))
%       *** END ELEMENT GENERATION

% ========================================================================
%                          PHASE 4
%         APPLY BC AND CONSTRAINTS (Make problem unique)
% ========================================================================
fprintf ('Begin Phase 4, modify matrices for EBC \n')

%    ALLOCATE STORAGE FOR OPTIONAL REACTION RECOVERY
if ( EBC_count > 0 )                      ; % reactions occur
  % EBC_row = row from S for each EBC, for later reaction calculation
  % EBC_col = row from C for each EBC, for later reaction calculation
  EBC_row = zeros (EBC_count, n_d)        ; % reaction data
  EBC_col = zeros (EBC_count, 1)          ; % reaction data
  [EBC_row, EBC_col] = save_reaction_matrices (EBC_flag, S, C);
end % if essential BC exist (almost always true)

% in the future, inforce all MPC equations here

% ENFORCE ESSENTIAL BOUNDARY CONDITIONS (and insert identities)
if ( EBC_count > 0 )                       ; % reactions occur
  [S, C] = enforce_essential_BC (EBC_flag, EBC_value, S, C);
end % if essential BC exist (almost always true)

if ( show_bc == 1 ) % then list each system matrix after EBC
  fprintf ('After EBC system matrix S = \n') ; disp(S)
  fprintf ('After EBC system matrix C = \n') ; disp(C)
end % if show_s
% ========================================================================
%                          PHASE 5
%         SOLVE LINEAR EQUATION SYSTEM AND ITS REACTIONS, OR
%                  SOLVE EIGENVALUE SYSTEM
% ========================================================================
fprintf ('Begin Phase 5, solve modified system \n')
Ans = S \ C                       ; % Compute nodal unknowns
list_save_results (n_g, n_m, Ans) ; % save for plots and print

if ( color == 1 ) % display and save a png plot
  addpath /clear/www/htdocs/mech517/Matlab_Plots              % library
  color_scalar_result (n_e, n_m, x, y, Ans, el_type, nodes) ; % plot 
  pause (4)                                  ; % pause for 4 seconds
end % if plot desired

%    OPTIONAL FORCE REACTION RECOVERY & SAVE
if ( EBC_count > 0 )                    ; % reactions exist ?
  % EBC_react = nodal source required to maintain EBC value
  EBC_react   = zeros (EBC_count, 1)    ; % pre-allocate
  [EBC_react] = recover_reactions_print_save (n_g, n_d, ...
                EBC_flag, EBC_row, EBC_col, Ans); % reaction to EBC
end % if EBC existed

% =====================================================================
%                             PHASE 6
%                      POST-PROCESS ELEMENTS 
%  Get location and values of derivaties at all q-point, and/or 
%  Get integral of the solution and/or
%  Sum all the element or nodal values
% =====================================================================
fprintf ('Begin Phase 6, use the answers \n')

if ( post == 1 ) % then post-process all element stresses
  if ( qp_unit > 0 ) ; % close GQ written files
    fclose (qp_unit) ; % file qp_store.bin
    qp_unit = fopen ('qp_store.bin', 'r') ; % open for reading 
  end % if
  if ( el_unit > 0 )                    ; % close written binary files
    fclose (el_unit)                             ; % file el_store.bin
    el_unit = fopen ('el_store.bin', 'r')         ; % open for reading
  end % if element arrays were saved
  if ( el_unit == -1 )                                      % open failed
    fprintf ('\nWARNING: post-processing unit open failed \n')
  end % if unit available
  if ( qp_unit == -1 )                                      % open failed
    fprintf ('\nWARNING: post-processing unit open failed \n')
  end % if unit available

% Save temp file for flux plot input (see Matlab quiver_* plot files)
  q_id = fopen ('el_qp_xyz_fluxes.txt', 'w') ; % flux plot file

  fprintf ('\n')
  fprintf ('Element Post-processing: \n')
  strain = zeros (n_r, 1) ; % pre-allocate strains
  stress = zeros (n_r, 1) ; % pre-allocate stresses

%  Re-compute or read arrays from storage
  format short e
  for j = 1:n_e ; % loop over elements ====>> ====>> ====>> ====>> ====>> 
    type    = el_type (j)            ; % current element type number
    typ_n   = type_n (type)          ; % maximum nodes per el type
    typ_i   = type_i (type)          ; % maximum DOF per el type
    if ( type ~= 1 ) % then skip this element and go to the next one
      fprintf ('NOTE: skipped type %i element number %i \n', type, j)
      continue % to the next j ====>> ====>> ====>>
    end % skip this type

%       Gather answers
    e_nodes = nodes (j, 1:typ_n)                    ; % connectivity 
    [rows]  = get_element_index (n_g, n_n, e_nodes) ; % eq numbers 
    Ans_e (1:n_i) = Ans(rows)                 ; % gather element results

%  Recover previously calculated element arrays
    [xy_e]    = fread (el_unit, [n_n, n_s],  'double'); % coordinates
    [el_prop] = fread (el_unit, [n_vals, 1], 'double'); % properties
    [S_e]     = fread (el_unit, [n_i, n_i],  'double'); % beam stiffness
    [C_e]     = fread (el_unit, [n_i, 1],    'double'); % line load resultant
    [H_i]     = fread (el_unit, [n_i, 1],    'double'); % interpolation integral
    [H_ix]    = fread (el_unit, [n_i, 1],    'double'); % x*interpolation integral
    [H_iy]    = fread (el_unit, [n_i, 1],    'double'); % y*interpolation integral

    for q = 1:n_q   ; % begin quadrature loop ---> ---> ---> ---> ---> --->

%    Recover previous calculations at this point (or re-compute)
      [xy_q] = fread (qp_unit, [1, n_s],   'double') ; % get coordinates
      [E_e]  = fread (qp_unit, [n_r, n_r], 'double') ; % get material matrix%B
      [B_q]  = fread (qp_unit, [n_r, n_i], 'double') ; % get application def

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

%       List quadrature point results
      if ( show_qp == 1 ) % then list data
        fprintf ('\n')
        if ( n_s == 1 ) % then 1D
          fprintf ('El, Q_Pt, Coordinates             %i, %i, %g \n', ...
                    j, q, xy_q (1))
          fprintf ('El, Q_Pt, Pressure Gradient       %i, %i, %g \n', ...
                    j, q, strain (1))
          %B fprintf ('El, Q_Pt, Shear stress components %i, %i, %g \n', ...
                    %B j, q, stress(1))
        else % then 2D
          fprintf ('El, Q_Pt, Coordinates             %i, %i, %g %g \n', ...
                    j, q, xy_q (1:2))
          fprintf ('El, Q_Pt, Pressure Gradient       %i, %i, %g %g \n', ...
                    j, q, strain (1:2))
          %B fprintf ('El, Q_Pt, Shear stress components %i, %i, %g %g \n', ...
                    %B j, q, stress(1:2))
         end % if 1D or 2D
       end % if list as well as save for plots

%      save flux for plot programs
       if ( n_s == 1 ) % then 1D
         fprintf (q_id, '%g %g \n',  xy_q (1), strain (1));
       else
         fprintf (q_id, '%g %g %g %g \n',  xy_q (1:2), strain (1:2));
       end % if 1D or 2D
    end % for q quadrature points <--- <--- <--- <--- <--- <--- <--- <---

    if ( integral == 1 ) % report solution integral
      Type_integral   (type) = Type_integral   (type) + dot (Ans_e, H_i) ;
      Type_integral_x (type) = Type_integral_x (type) + dot (Ans_e, H_ix) ;
      Type_integral_y (type) = Type_integral_y (type) + dot (Ans_e, H_iy) ;
    end % if scalar solution

  end % for each j element in mesh   <<==== <<==== <<==== <<==== <<==== 

% close the flux data plot file
  fclose(q_id)
  if ( vector == 1 ) % plot flux vectors at selected quadrature points
    addpath /clear/www/htdocs/mech517/Matlab_Plots  % library
    quiver_qp_flux_mesh (1.0, 1, 0) % inputs (scale, inc_g, sign)
    pause (4)                       ; % pause for 4 seconds
  end % if

  %B if ( integral == 1 & n_g == 1) % report solution integral
  if ( integral == 1 ) % report solution integral
    fprintf ('\n')
    if ( n_t > 1 ) % report on each type
      for j = 1:n_t
        fprintf ('Solution integral in type %i elements = %g \n', ...
                  j, Type_integral (j))
        fprintf ('Solution volume   in type %i elements = %g \n', ...
                  j, Type_volume (j))
      end % for each type
    end % if multiple element types
    fprintf ('Total vertical force = %g \n', sum (Type_integral  ))
    x_F = Type_integral_x (1) / Type_integral (1) ;
    y_F = Type_integral_y (1) / Type_integral (1) ;
    fprintf ('Force located at (x, y) = (%g, %g) \n', x_F, y_F)
    fprintf ('Total moment about y = %g \n', sum (Type_integral_x))
    fprintf ('Total moment about x = %g \n', sum (Type_integral_y))
  end % if scalar solution
  fprintf ('\n Total bearing area  = %g \n', sum (Type_volume))
  format short

%       OPTIONAL POST PROCESSING COMPLETED
else
  fprintf ('\nNo element post-processing requested \n')
end % for post-processing stresses

% =====================================================================
%                             PHASE 7
%           ESTIMATE ERROR IN EACH ELEMENT AND THE SYSTEM
% =====================================================================
%     Not yet implemented
% end Lubrication_2D % =====================================