function [] = Variable_Coeff_GQ_lib ()       % Revised 02/12/19 
% Requires: addpath /clear/www/htdocs/mech517/Akin_FEA_Lib
% --------------------------------------------------------
% The 1D solution of the variable coefficient ODE of
%     -[K(x) u']' + A(x) u' + C(x) u - Q(x) = 0
% with matrix system (S + U + M) u = c + c_NBC
%
%  via numerically integrated element matrices and using
%  four- (n_n=4), or three- (n_n=3), or two-node (n_n=2)
%         line elements using element properties
% 1---2 -->r or 1---2---3 -->r, or 1---2---3---4 -->r
%
% If all coefficients are constant, then there is one
% property line with K, A, C, Q. Otherwise, the 
% element properties are: n_n K(x) values, and
%            n_n A(x) values, n_n C(x) values, and
%            n_n Q(x) values at the n_n element nodes 
% In other words, the mesh properties data array contains
% 4 * n_n values for each element in the mesh.
% --------------------------------------------------------
%     Mech 417, Finite Element Analysis, Rice University 
%     Copyright J.E. Akin 2019. All rights reserved.
%......................................................... 
% echo user information about this specific application
echo_user_remarks

% =================================================================
%     *** PHASE 1: SET OR READ USER PROBLEM DEPENDENT CONTROLS *** 
% =================================================================
%  set user maximum array sizes 
n_g = 1     ; % number of DOF per node 
N_N = 2     ; % number of EXPECTED nodes per element
n_p = 1     ; % dimension of parametric space
n_q = 2     ; % number of quadrature points
n_r = 1     ; % number of rows in element B_e matrix
N_S = 1     ; % EXPECTED physical space dimensions

%   set user logic choices (or add read function)
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 integral 
load_pt  = 0; % turn on(1)/off(0) point source (value can be zero)
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   = 3; % turn on(>0)/off(0)list some S_e, M_e and c_e
show_qp  = 1; % turn on(1)/off(0) list qp post-processing results
show_s   = 1; % turn on(1)/off(0) listing of assembled S and c

% Is post-processing required ?
if ( el_react > 0 | integral > 0 | show_qp > 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 ('echo_p   listing of properties         = %i \n', echo_p  )
fprintf ('el_react listing of element reactions  = %i \n', el_react)
fprintf ('integral list integral of the solution = %i \n', integral)
fprintf ('load_pt  point source (can be zero)    = %i \n', load_pt )
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, M_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  )

%    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 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)
fprintf ('\n')

% =====================================================================
%        ***  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)
% 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, coordinates 

%   Compare actual data against space control
if ( n_s ~= N_S ) % then possible control error
  fprintf ('WARNING: Read %i coordinates; expected %i \n', n_s, N_S)
end % if
 
% *** Read mesh element types and nodal connections ***
% (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 
% 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/element; expected %i \n', n_n, N_N)
end % if inconsistent inputs

%    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);

if ( n_vals ~= 4*n_n )                  ; % then invalid data set
  fprintf ('Require %i properties, but read %i \n', 4*n_n, n_vals)
  error   ('Missing variable coefficients data for each element')
end                             ; % if valid number of properties

% Element properties are: n_n K(x) values, 
%            n_n A(x) values, n_n C(x) values, and
%            n_n Q(x) values at the n_n element nodes 

fprintf ('%i K(x) values, and %i A(x) values \n', n_n, n_n)
fprintf ('%i C(x) values, and %i Q(x) values \n', n_n, n_n)

%   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       = zeros (n_d, 1)   ; % results
c         = zeros (n_d, 1)   ; % force or source
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
S         = zeros (n_d, n_d) ; % generalized stiffness

%    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

% Read essential boundary condition values (from file msh_ebc.txt)
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 mesh general plot
fprintf ('\n') ; fprintf ('Creating curve plot, pausing \n')
plot_input_curve_mesh (n_e, n_m, n_n, n_s, x, y, z, nodes)

% ==============================================================
% ** PHASE 3: GENERATE ELEMENT ARRAYS & 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 element matrices \n', show_e)
end                 ; % if element array valuess to be displayed
 
%     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 data
  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

type_old = 1                  ; % initialize element type number
for k = 1:n_e ; % loop over elements ====>> ====>> ====>> ====>>
  type    = el_type (k)          ; % current element type number
  e_nodes = nodes (k, 1:n_n)                    ; % connectivity 

%   Get quadrature rule for current element type
  [r_q, s_q, t_q, w_q] = get_quadrature_rule (n_n, n_p, n_q) ;
   
%                 Gather nodal coordinates
  xy_e (1:n_n, 1) = x(e_nodes(1:n_n))    ; % x coord at el nodes 

%     Allocate and clear element type arrays
  Jacobian = zeros (n_s, n_s)              ; % geometry Jacobian 
  c_e = zeros (n_i,   1)              ; % clear array el sources
  c_r = zeros (n_i,   1)            ; % clear array el reactions
  H_i = zeros (1,   n_n)                 ; % clear integral of H
  B_q = zeros (n_r, n_i)               ; % B at quadrature point
  M_e = zeros (n_i, n_i)       ; % clear generalized mass matrix
  S_e = zeros (n_i, n_i)            ; % clear array el stiffness
  A_e = zeros (n_i, n_i)            ; % clear array el advection
  Kx  = zeros (n_n,   1)                   ; % variable property
  Ax  = zeros (n_n,   1)                   ; % variable property
  Cx  = zeros (n_n,   1)                   ; % variable property
  Qx  = zeros (n_n,   1)                     ; % variable source 

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

% Shorthand type property names  K(x) A(x) C(x), and Q(x) values
  Kx (1:n_n) = el_prop (1:n_n)        ; % K at each element node
  Ax (1:n_n) = el_prop (n_n+1:2*n_n)    ; % one A value per node
  Cx (1:n_n) = el_prop (2*n_n+1:3*n_n)  ; % one C value per node
  Qx (1:n_n) = el_prop (3*n_n+1:4*n_n)  ; % one Q value per node

% Numerical Intergation of S_e, c_e , in 1D parametric
  volume = 0                                  ; % element volume
  for q = 1:n_q     ; % begin quadrature loop --> ---> ---> --->
    r = r_q (q) ; w = w_q (q)       ; % recover integration data

%   Element scalar interpolation, H, & local derivatives, DLH
    [H_q, DLH_q] = Lagrange_1D_library (n_n, r)     ; % evaluate

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

%   Interpolate material properties at quadrature point
    Kx_q = H_q * Kx               ; % interpolate first property
    Ax_q = H_q * Ax              ; % interpolate second property
    Cx_q = H_q * Cx               ; % interpolate third property
    Qx_q = H_q * Qx                  ; % interpolate source data

%   Compute the geometry mapping data
    Jacobian = DLH_q * xy_e (1:n_n, 1)    ; % geometric Jacobian
    J_det = det (Jacobian)                       ; % determinant
    J_Inv = inv (Jacobian)                  ; % Jacobian inverse

%   Global derivatives of H are DGH
    DGH_q = J_Inv * DLH_q                              ; % dH/dx 

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

%   Update element square matrix, from K(x) 
    S_e = S_e + (B_q' * Kx_q * B_q) * J_det * w_q (q); % conduct

%   Update element square matrix, from A(x)
    A_e = A_e + (B_q' * Ax_q * H_q) * J_det * w_q (q) ; % advect 

%   Generalized mass matrix, from C(x)
    M_e = M_e + (H_q' * Cx_q * H_q) * J_det * w_q (q); % convect

%   Source vector, from Q(x)
    c_e = c_e + H_q' * Qx_q * J_det * w_q (q)         ; % source

%   Update integral of the interpolations for post-processing
    H_i    = H_i    + H_q * J_det * w_q (q)       ; % H integral
    volume = volume + J_det * w_q (q)         ; % element volume

%   save data needed for post-processing at this point
    if ( post == 1 )
      fwrite (qp_unit, xy_q,  'double')     ; % save coordinates
      fwrite (qp_unit, Kx_q,  'double')       ; % save materials
      fwrite (qp_unit, B_q,   'double')           ; % save dH/dx
     end                              ; % if post-process needed
   end ; % for quadratures <-- <-- <--- <--- <--- <--- <--- <---

   if ( k <= show_e )  ; % only show the first few elem matrices
     fprintf ('Element matrix S_e = \n') ; disp(S_e)
     fprintf ('Element matrix A_e = \n') ; disp(A_e)
     fprintf ('Element matrix M_e = \n') ; disp(M_e)
     fprintf ('Element matrix c_e = \n') ; disp(c_e)
   end ;                     % if element matrix displays wanted

 %  Save element arrays needed for post-processing gradients
    if ( post == 1 )            ; % save arrays for this element
      fwrite (el_unit, S_e,     'double')     ; % save stiffness
      fwrite (el_unit, c_e,     'double') ; % save source result
      fwrite (el_unit, H_i,     'double') ; % save integral of H
    end                              ;  % if save element arrays

%         SCATTER TO (ASSEMBLE INTO) SYSTEM ARRAYS 
% Insert completed element matrices into system matrices 
   S_e    = S_e + A_e + M_e                ; % net square matrix
   fwrite (el_unit, S_e,     'double')         ; % for el_reacts 

  % rows = vector subscript to convert elem to system eq numbers
   [rows] = get_element_index (n_g, n_n, e_nodes) ; % eq numbers 
   S (rows, rows) = S (rows, rows) + S_e    ; % add to stiffness 
   c (rows)       = c (rows)       + c_e   ; % add to sys source
end ; % for each k 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))  ; % check

% ========================================================================
%           *** PHASE 4: APPLY BC AND CONSTRAINTS ***
% ========================================================================
fprintf ('\n'); 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)

% 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 ***
% ========================================================================
fprintf ('\n'); 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

graph_1d_result (n_e, n_m, n_n, x, nodes, Ans) ; % graph

%    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 AT EACH INTEGRATION POINT **
% ===============================================================
fprintf ('Begin Phase 6, use the answers \n')
if ( post == 1 )      ; % then post-process all element flux data
  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 unit is open
  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

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

  fprintf ('\n')                                        ; % space
  fprintf ('Element Post-processing: \n')        ; % print header
  fprintf ('Point location, approximate flux \n'); % print header
  if ( el_react )         ; % allocate arrays to graph end fluxes
    Reacts = zeros (n_e, 3)            ; % end reactions and loss
  end                     ; % if 1D element end reactions desired

%  Read arrays from binary storage
  format short e
  Total_i  = 0         ; % initialize intregral of solution terms
  for k = 1:n_e ; % loop over elements ===>> ====>> ====>> ====>>
    type = el_type (k)            ; % current element type number

%   Gather answers for this element
    e_nodes = nodes (k, 1:n_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 arrays for this element
    [S_e]     = fread (el_unit, [n_i, n_i],  'double')  ; % stiff
    [c_e]     = fread (el_unit, [n_i, 1],    'double') ; % source 
    [H_i]     = fread (el_unit, [1, n_n],    'double');% integral
    [S_t]     = fread (el_unit, [n_i, n_i],  'double')  ; % total
  
    Total_i = Total_i + H_i * Ans_e'     ; % integral of solution

    for q = 1:n_q ; % begin quadrature loop   ---> ---> ---> ---> 
%     Recover previous calculations at this point 
      [xy_q] = fread (qp_unit, [1, n_s],   'double');% get coords
      [Kx_q] = fread (qp_unit, [1, 1],     'double')    ; % get K
      [B_q]  = fread (qp_unit, [n_r, n_i], 'double') ; % operator 

%     Compute the generalized gradient and flux 
      Grad = B_q * Ans_e(1:n_i)'                    ;  % gradient
      Flux = Kx_q * Grad                     ; % application flux

%     List quadrature point results
      if ( show_qp > 0 )      ; % then list the gradient and flux 
       fprintf('\n')                               ; % line space
       fprintf('El, Pt, Coordinates . %i, %i, %g \n', k, q, xy_q)
       fprintf('El, Pt, Gradient .... %i, %i, %g \n', k, q, Grad)
       fprintf('El, Pt, Flux ........ %i, %i, %g \n', k, q, Flux)
      end                                            ; % if print

%      save flux for plot programs
       fprintf (q_id, '%g %g \n',  xy_q (1), Flux (1))   ; % save
    end ; % for q quadrature points <--- <--- <--- <--- <--- <---

%   Optional output of nodal reactions on every element
    if ( el_react ); % get end reactions (flux) for this element
%     get the element reactions and list them
      c_m = S_t * Ans_e' - c_e; % node reaction values on element
      Reacts (k, 1) = c_m(1); Reacts (k, 2) = c_m (n_n) ; % ends
    end                                      ; % if el_react list
  end ; % for each k element in mesh  <<==== <<==== <<==== <<====

  if ( el_react ); % get end reactions (flux) for this element
    Reacts (:, 3) = Reacts (:, 2) + Reacts (:, 1)   ; % convected
    fprintf ('\n')
    fprintf ('Element end reactions and flux loss \n')
    fprintf ('NOTE: - enters element, + leaves element \n')
    fprintf ('El, Left        Right       Loss or Gain \n')
    for k = 1:n_e ; % loop over elements ----> ---->-----> ----->
     fprintf ('%i  %10.3e  %10.3e  %10.3e \n', ...
               k, Reacts (k, 1:3) )   ; % end flux & convect loss
    end ; % for each k element in mesh  <----- <----- <---- <----
  end                                 ; % if el_react & loss list

  if ( integral > 0 )       ; % List the integral of the solution
    fprintf ('Total solution integral = %g \n', Total_i); % value
  end                        ; % if solution integral is required

% close the flux data plot file
  fclose(q_id)                  ; % close the flux data plot file
end                             ; % for post-processing flux data
% end % Variable_Coeff_GQ_lib ===================================