% Modular Finite Element Routines via Matlab
% Copyright J.E. Akin, Rice University
% Variable   Description or [Default] and size (rows, cols)
%------------------------------------------------------------------
% Ans       = solution vector (n_d, 1) from S * Ans = C
% Ans_e     = element solution vector subset of Ans (n_i, 1)
% Ans_rect  = Ans  reshaped as a rectangle for output (n_g, n_m)
% axisym   [0] turn on(1)/off(0) an axisymmetric domain flag
% B_q       = element application B matrix at point q (n_r, n_i)
% Body      = body force vector, or volumetric source (n_s, 1)
% Body_x    = body force component
% Body_y    = body force component
% Body_z    = body force component
% C         = total source vector, (n_d, 1)
% C_e       = element type source vector (n_i, 1)
% C_f       = elastic foundation reaction loading
% C_m       = S_e * T_e' - C_e' = applied loading reactions (n_i, 1)
% C_P       = penalty column matrix (MPC_max, 1)
% C_pt      = vector of input point sources, adds to C
% C_r       = reaction vector of an element
% C_react   = C = a copy for later reaction use (n_d, 1)
% ctrl_type = vector of el type descriptive data (type_C)
% ctrl_sys  = vector of system control and descriptive data (sys_C)
% domain    = flag element & application type, domain = n_p*100 + n_g*10 + n_n
% DGH_q     = element interpolation global derivatives at q (n_s, n_n)
% DLH_q     = element interpolation local  derivatives at q (n_p, n_n)
% DLH       = element interpolation local  derivatives at r,s,t (n_p, n_n)
% DOF       = an integer degree of freedom number
% E         = Modulus of elasticity
% E_e       = element constitutive matrix (n_r, n_r)
% E_q       = element constitutive matrix at point q (n_r, n_r)
% EBC       = current essential boundary condition value
% EBC_col   = row from C for each EBC (EBC_count, 1)
% EBC_count = number of essential (Dirichlet) BC's
% EBC_flag  = n_g individual binary flags at each node (n_m, n_g)
% EBC_react = nodal source required to maintain EBC value (EBC_count, 1)
% EBC_row   = row from S for each EBC (EBC_count, n_d)
% EBC_value = values assigned to essential (Dirichlet) BC's (n_m, n_g)
% e_nodes   = maximum connectivity list for any element type (n_n)
% e_t       = loop variable for element type 
% echo_p   [0] turn on(1)/off(0) listing of properties
% el_prop   = property values of current element material (n_vals, 1)
% el_type   = type number of every element, usually 1 (n_e, 1)
% el_vol    = volume of a single element
% F_sq      = Eigenvalue frequency, in rad/sec, squared (n_d, 1)
% file_id   = a temporary file id number
% flag_EBC  = reshape of EBC_flag (1, n_d)
% flags     = n_g individual EBC or MPC integer flags at a node (n_g, 1)
% forward   = factorization work vector (n_d, 1)
% H         = element interpolation functions at r,s,t (n_n, 1)
% H_i       = integral of interpolation functions over element (n_n, 1)
% H_q       = element interpolation functions at q (n_n, 1)
% had_n     = previous number of active nodes per element
% has_n     = actual number of active nodes per element, max(type_ctrl(:,1))
% has_p     = max (type_ctrl (:, 2)), actual max parm dim   
% has_q     = max (type_ctrl (:, 3)), actual max Gauss rule 
% has_x     = max (type_ctrl (:, 4)), actual max geom nodes
% Ignore    = flag for setting eigenvalues to ignore = = -987654
% index     = vector subscript converting elem to system eq numbers
% i_p       = a single parameter (DOF) number
% J_det     = determinant of Jacobian
% J_inv     = inverse of Jacobian
% Jacobian  = geometry spatial Jacobian (n_s, n_s)
% L_S       = lower triangle factorization from S = L_S * U_S
% Line_e    = current element line load values
% last_typ  = previous element type number
% line      = temporary data for one MPC, (3*MPC_max, 1)
% load_pt  [0] turn on(1)/off(0) point sources input
% M         = mass matrix, if any (n_d, n_d)
% M_e       = element type mass matrix (n_i, n_i)
% M_g       = generalized element type mass matrix (n_i, n_i)
% M_pt      = vector of input point masses, adds to diagonal of M
% mass_pt  [0] turn on(1)/off(0) point masses input
% maxTau    = maximum shear stress
% MPC_coeff = the dof coefficients in each constraint equation
% MPC_count = number of algebraic multipoint constraint eqs
% MPC_dofs  = the dof numbers in each constraint equation
% MPC_each  = count constraints with 1 to MPC_max DOF terms (MPC_max, 1)
% MPC_max   = maximum number of dof in a single constraint equation
% mix_0     = work term in calculating mix_2
% mix_1     = Mixed BC data: K_n * U,n + mix_1 * U + mix_2 = 0
% mix_2     = Mixed BC data: K_n * U,n + mix_1 * U + mix_2 = 0
% Nt_N      = element type mass or foundation or convection matrix
% n_a       = maximun number of adjacent elements
% n_b       = number of boundary segements with sources
% n_c       = number of constraint equations
% n_d       = n_g * n_m = system degrees of freedom (DOF)
% n_e       = number of elements in the mesh
% n_f       = maximum number of flux components
% n_g       = number of DOF per node (deflection, slope, etc) 
% n_h       = number of scalar interpolations [n_n]
% n_i       = n_g * n_n = max number of DOF per element
% n_j
% n_k
% n_l       = max number of elements in a patch
% n_m       = number of nodal points in mesh
% n_mats    = number of materials, = n_e or 1 if homogeneous
% n_mix     = number of mixed (Robin) boundary conditions
% n_modes   = number of modes and frequencies to be output
% n_n       = maximum number of nodes per element
% n_o
% n_p       = dimension of parametric space [n_s]
% n_prop    = number of properties per material
% n_q       = number of quadrature points required 
% n_r       = number of rows in application dependent B_e matrix 
% n_s       = number of physical space dimensions
% n_sky     = size of stiffness matrix in spare skyline storage
% n_steps   = number of times steps
% n_t       = number of different element types
% n_u       = number of point sources
% n_unkn    = number of eigenvalue DOF = n_d - EBC_count
% n_v       = number of vector interpolations [n_n * n_g]
% n_vals    = number of property values for each material
% n_w
% n_x       = number of element geometry nodes [n_n]
% n_y
% n_z
% node      = an integer node number
% nodes     = node connection list for each element (n_e, n_n)
% nu        = Poisson's ratio
% Order     = A vector subscript to a sorted list (list_size, 1)
% PD        = packed EBC flag [n_g digits 0 to 9], for each node (n_m, 1)
% patch    [0] turn on(1)/off(0) patch test option
% post     [0] turn on(1)/off(0) element post processing
% Q         = volumetric source, usually = Body (1)
% q         = loop variable for quadrature points
% q_id      = unit number to save generalized strains and stress at q
% qp_unit   = system unit for qradrature write/read (when post=1)
% R_e       = rare rectangular array
% React     = reaction value for a single DOF
% Rho       = mass density
% rows      = vector subscript converting elem to system eq numbers
% r         = first parametric quadrature coordinate at q        
% r_q       = all of first parametric quadrature coordinates (n_q, 1)
% S         = full storage stiffness matrix (n_d, n_d)
% S_sky     = skyline storage stiffness matrix (n_sky, 1)
% s         = second parametric quadrature coordinate at q        
% s_q       = all of second parametric quadrature coordinates (n_q, 1)
% S_e       = element type stiffness matrix (n_i, n_i)
% S_P       = penalty square matrix (MPC_max, MPC_max)
% S_pt      = vector of input point stiffnesses, adds to diagonal of S
% stiff_pt [0] turn on(1)/off(0) point stiffnesses input
% strain    = strain components (n_r, 1)
% strain_0  = initial strain components (n_r, 1)
% stress    = stress components (n_r, 1)
% sys_C     = number of control items per application = 10  
% Thick     = element thickness, planar or axisymmetric
% Type_integral = integral of solution in each element type group
% Type_volume = volume of each element type group
% Totals    = total of each DOF component in mesh (1, n_g)
% t         = third parametric quadrature coordinate at q        
% t_i       = t_n * n_g = actual type dof count
% t_n       = number of nodes for current element type
% t_nodes   = actual connectivity for current element type (type_n)
% t_q       = all of third parametric quadrature coordinates (n_q, 1)
% thk_q     = thickness at quadrature point q, Thick iff constant
% type      = current element type number
% type_ctrl = element type controls data (n_t, type_C) 
% type_C    = number of control items per element type, = 6
% type_S    = el type shape: 1=line 2=tri 3=quad 4=hex 5=tet 6=wedge
% type_i    = n_g * type_n = max number of DOF per element type
% type_n    = maximum number of nodes per current element type
% type_p    = dimension of current el type parametric space [n_s-1]
% type_q    = number of quadrature points required by el type
% type_x    = number of geometry nodes for current el type [type_n]
% U_S       = upper triangle factorization from S = L_S * U_S
% value_EBC = reshape of EBC_value (1, n_d)
% vonMise   = von Mises stress at a point
% w         = weight at quadrature point q        
% w_q       = weight at all quadrature points in n_p (n_q, 1)
% x         = mesh x-coordinate values (n_m, 1)
% xy_e      = system coordinates at el_nodes of element (n_n, n_s)
% xy_q      = system coordinates at point q (ns, 1)
% y         = mesh y-coordinate values, if present (n_m, 1)
% z         = mesh z-coordinate values, if present (n_m, 1)
%------------------------------------------------------------------