function L3_Quintic_BoEF_MPC (load_pt) % Revised 3/1/09 =========
% load_pt = 1 flags that point forces and/or couples will be  
%  input via text file msh_load_pt.tmp 
%............................................................. 
%     quintic beams (with optional elastic foundation),
%     for point loads & couples, linear line load 
%     Mech 417, Finite Element Analysis, Rice University 
%     Copyright J.E. Akin 2009. All rights reserved.
%............................................................. 
if ( nargin == 0 )          ;  % check for optional data flag 
  load_pt = 0               ;  % no point source data 
end % if from argument count 
if ( load_pt ~= 1 )         ; % check user intent
  fprintf ('NOTE: no point load data supplied. Reset flag? \n')
end % if 
 
% 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 USER PROBLEM DEPENDENT CONTROLS ***
% =====================================================================

 %  set application (C1 beam) controls
n_g = 2     ; % number of DOF per node (deflection, slope) 
n_q = 0     ; % number of quadrature points required 
n_r = 1     ; % number of rows in B_e matrix 
 
% =====================================================================
%                             PHASE 2
%            *** BEGIN DATA INPUT AND CHECKING PHASE ***
% =====================================================================
%              Read mesh nodal input data files 
[n_m, n_s, P, x, y, z] = get_mesh_nodes ; % bc_flags, coordinates 
% n_m  = number of nodal points in mesh
% n_s  = number of space dimensions
% P    = nodal packed essential BC flag (n_g digits 0 or 1)
% x    = mesh x-coordinate values
% y    = mesh y-coordinate values, if present, etc. for z
 
%      Extract EBC flags from packed integer flag P  
EBC_flag   = zeros (n_m, n_g)    ; % initialize all DOF flags
[EBC_flag] = get_ebc_flags (n_g, n_m, P)     ; % unpack flags 
EBC_count  = sum( sum ( EBC_flag == 1 ) )    ; % number of EBC 
% EBC_flag  = n_g individual binary flags at each node
% EBC_count = number of essential (Dirichlet) BC's
if ( EBC_count == 0 ) 
  fprintf ('WARNING: No displacement boundary condition given. \n') 
else 
  fprintf ('Note: expecting %i displacement BC values. \n', EBC_count) 
end % if  

%     Extract any MPC flags from packed integer flag P
MPC_max = max( max ( EBC_flag ) )         ; % max type MPC
if ( MPC_max > 1 ) % then algeabraic constraints flagged
  MPC_each (1:MPC_max) = 0  ; % number of eqs of each type
  for k = 1:MPC_max                          % count each flag digit
    MPC_each (k) = sum( sum ( EBC_flag == k)) ; % total of each flag 
  end % for
  for k = 1:MPC_max            % convert to number of EBC or MPC eqs
    check        = mod ( MPC_each (k), k) ; % check consistent flags
    if ( check ~= 0) 
      fprintf ('WARNING: number of flags with digit %i ', k)
      fprintf ('is not divisible by %i \n', k)
    end % if consistent
    MPC_each (k) = ceil (MPC_each (k) / k) ; % total of each MPC type 
  end % for
  MPC_count  = sum( MPC_each (2:MPC_max)) ; % number of MPC equations
  if ( MPC_count > 0 )  
    fprintf ('Note: expecting %i algebraic MPC sets. \n', MPC_count) 
    fprintf (' with a maximum dof per equation of %i \n', MPC_max  ) 
  end % if  
else
  MPC_count = 0 ;
end % if algeabraic constraints flagged
% MPC_count = number of algebraic multipoint constraint eqs
% MPC_max   = max number of DOF in a single MPC equation, now = 2

%   Read mesh element types and nodal connections
%             (from msh_typ_nodes.tmp)

[n_e, n_n, n_t, el_type, nodes] = get_mesh_elements  ; 
% 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

%   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.
S = zeros (n_d, n_d) ; M = zeros (n_d, n_d) ; % stiff & mass 
C = zeros (n_d, 1)                          ; % force & moments 
T = zeros (n_d, 1)                          ; % displacements 
C_react = zeros (n_d, 1)  ; % copy of C for reaction recovery
% S = stiffness matrix,       C = total soutce vector
% M = mass matric (if any),   T = solution vector

%   Read essential boundary condition (EBC) values, if any
if ( EBC_count > 0 )                          ; % need EBC data 
  EBC_value = zeros (n_m, n_g)                ; % default to zero
  % read values from file msh_ebc.tmp
  [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag) ; % read data 
  % EBC_value = EBC value for each restrained DOF
end % if any EBC data expected 
 
%    Read point loads & couples data, if any, and insert in C 
if ( load_pt > 0 )                         ; % need point loads data 
  %  read pt sources from msh_load_pt.tmp, add element sources later
  [C] = get_and_add_point_sources (n_g, n_m, C);   % add point loads 
end % if any point source expected 
C_react = C ; % save before BC for later reaction use
 
%    Read any material properties for the elements
%         from optional file msh_properties.tmp
load msh_properties.tmp                       ; % one row per element 
fprintf ('\n(Echoing file msh_properties.tmp) \n')
n_prop = size(msh_properties, 1)              ; % == 1 if homogeneous 
n_vals = size(msh_properties, 2)        ; % of application properties
% n_prop = number of properties, = n_e or 1 if homogeneous
if ( n_prop == 1 )        % then element data are honogeneous
  fprintf ('%i homogeneous property values \n', n_vals)
else
  fprintf ('%i property values per element \n', n_vals)
end % if
disp(msh_properties)
if ( n_prop > n_e ) 
  error ('ERROR: number of property sets exceeds number of elements') 
end % if 

%   for beam application only
if ( n_prop == 1 )        % then element data are honogeneous
  E            = msh_properties (1, 1)   ; % material modulus 
  I            = msh_properties (1, 2)   ; % section inertia 
  Line_e (1:2) = msh_properties (1, 3:4) ; % trapezoidal line load 
  k_f          = msh_properties (1, 5)   ; % foundation stiffness 
  Rho          = msh_properties (1, 6)   ; % mass per unit length 
%                   ECHO PROPERTIES 
  echo = 0 ;  % turn off/on listing of properties
  if ( echo )
    fprintf (' \n') 
    fprintf ('Homogeneous Element Properties: \n' ) 
    fprintf ('Elastity modulus (N/m^2)     = %g \n', E)      
    fprintf ('Moment of inertia (m^4)      = %g \n', I) 
    fprintf ('Line Load (N/m)              = [ %g %g ] \n', ... 
              Line_e(1), Line_e(2)) 
    fprintf ('Foundation stiffness (N/m^2) = %g \n', k_f) 
    fprintf ('Mass per unit length (Kg/m)  = %g  \n', Rho) 
  end % if echo
end % if homogeneous element properties
F_k (1:n_prop) = msh_properties (1:n_prop, 5); % foundation stiffness
if ( any ( F_k > 0 )) % BEF pressures will exist 
  BEF = 1 ;
  fprintf('WARNING: BoEF moment & shear inaccurate for crude meshes \n')
else % turn off foundation pressure calculation flag
  BEF = 0 ;
end % if reaction flag
%   end for beam application only

%  read any algebraic contraint eqs, like dof_1 + dof_2*c_1 = c_2
if (MPC_count > 0 ) % the multipoint constraint (MPC) sets needed
   MPC_coeff = zeros (MPC_max, MPC_count) ; % coeff values im MPC
   MPC_dofs  = zeros (MPC_max, MPC_count) ; %  dof numbers in MPC
   [MPC_coeff, MPC_dofs] = get_constraint_eqs (n_g, ...
                           MPC_each, MPC_max, MPC_count);
%  T_1 * A + T_2 * B = C, or T_1 + T_2 * B/A = C/A, etc
end % if MPC

%    *** END OF DATA INPUT AND CHECKING PHASE ***
 
% ======================================================================
%                              PHASE 3
%       ** GENERATE ELEMENT VALUES AND ASSYMBLE INTO SYSTEM **
% ======================================================================
%  Assemble system n_d by n_d square matrix terms from the
%         element n_e by n_e terms
 
for j = 1:n_e  ; % loop over elements ====>> ====>> ====>> ====>> 
  e_nodes = nodes (j, 1:n_n)                      ; % connectivity 
  type    = el_type (j)            ; % current element type number
   
%     Allocate and clear element arrays
  S_e = zeros (n_i, n_i)    ; % clear array el stiffness
  NtN = zeros (n_i, n_i)    ; % clear array el mass or foundation
  M_e = zeros (n_i, n_i)    ; % clear array el mass
  C_e = zeros (n_i, 1)      ; % clear array el sources

%%     Gather nodal coordinates
%  xy_e (1:n_n, 1) = x(e_nodes)             ; % x coord at el nodes

%               SET ELEMENT PROPERTIES & GEOMETRY 
  L = abs (x(e_nodes(n_n)) - x(e_nodes(1)))         ; % beam length 
  %b Line_e (1:2) = 0        ; % zero current element line load values
  if ( n_prop > 1 )    % each element has some different properties
    E            = msh_properties (j, 1)   ; % material modulus 
    I            = msh_properties (j, 2)   ; % section inertia 
    Line_e (1:2) = msh_properties (j, 3:4) ; % trapezoidal line load 
    k_f          = msh_properties (j, 5)   ; % foundation stiffness 
    Rho          = msh_properties (j, 6)   ; % mass per unit length 
%                   ECHO PROPERTIES 
    fprintf ('\nProperties for element %g \n', j) 
    fprintf ('Elastity modulus (N/m^2)     = %g \n', E)      
    fprintf ('Moment of inertia (m^4)      = %g \n', I) 
    fprintf ('Line Load (N/m)              = [ %g %g ] \n', ... 
              Line_e(1), Line_e(2)) 
    fprintf ('Foundation stiffness (N/m^2) = %g \n', k_f) 
    fprintf ('Mass per unit length (km/m)  = %g  \n', Rho) 
  end % if 
 
%   ELEMENT, FOUNDATION STIFFNESSES AND SOURCE RESULTANT MATRICES 
  L2 = L^2 ; L3 = L^3 ; % shorthand  
%  Constant quintic element & foundation stiffnesses & mass matrices 
  S_e = [5092,   1138*L, -3584,    1920*L,  -1508,   242*L  ; 
         1138*L, 332*L2, -896*L,   320*L2,  -242*L,  38*L2  ; 
        -3584,   -896*L,  7168,    0,       -3584,   896*L  ; 
         1920*L, 320*L2,  0,       1280*L2, -1920*L, 320*L2 ; 
        -1508,  -242*L,  -3584,   -1920*L,   5092,  -1138*L ; 
         242*L,  38*L2,   896*L,   320*L2,  -1138*L, 332*L2 ] ... 
         *E*I/(35 * L^3) ;                            % stiffness 
 
  NtN = L*[2092,   114*L, 880,  -160*L,   262,   -29*L  ; 
           114*L,  8*L2,  88*L, -12*L2,   29*L,  -3*L2  ; 
           880,    88*L,  5632,  0,       880,   -88*L  ; 
          -160*L, -12*L2,  0,    128*L2,  160*L, -12*L2 ; 
           262,    29*L,  880,   160*L,   2092,  -114*L ; 
          -29*L,  -3*L2  -88*L, -12*L2,  -114*L,  8*L2 ]/13860;  
 
  S_e = S_e + k_f * NtN ; % add foundation stiffness before assembly 
  M_e = Rho * NtN       ; % mass matrix 
 
%  Map line load to node forces & moments; C_e = w_To_F * Line_e  
  if ( any (Line_e) )            ; % then form forcing vector 
    w_To_F = L * [ 79,    19  ;  
                   5*L,   2*L ; 
                   112,   112 ; 
                  -8*L,   8*L ; 
                   19,    79  ; 
                  -2*L,  -5*L ] / 420 ; % quintic H, linear Line 
    % 6 x 2 * 2 x 1 = 6 x 1 result  
    C_e = w_To_F (1:n_i, 1:2) * Line_e' ; % force, moment @ nodes 
  end % if or set up resultant node loads for line load 
 
%         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 
  % rows = vector subscript converting element to system eq numbers
  S (rows, rows) = S (rows, rows) + S_e  ; % add to system stiffness 
  %b M (rows, rows) = M (rows, rows) + M_e  ; % add to system mass 
  C (rows)       = C (rows)       + C_e  ; % add to sys force/couples 
%       ***  END SCATTER TO SYSTEM ARRAYS ***
 
end % for each j element in mesh   <<==== <<==== <<==== <<==== <<==== 
%       *** END ELEMENT GENERATION
 
% ========================================================================
%                          PHASE 4
%         APPLY BC AND CONSTRAINTS (Make problem unique)
% ========================================================================
%    ALLOCATE STORAGE FOR OPTIONAL REACTION RECOVERY           
if ( EBC_count > 0 )                      ; % reactions occur 
  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); 
  save_resultant_load_vectors (n_g, C) 
  % EBC_row = row from S for each EBC
  % EBC_col = row from C for each EBC
end % if essential BC exist (almost always true) 

% Inforce all algebraic constraint equations (MPC)
if ( MPC_count > 0 )
  [S, C] = enforce_MPC_equations (MPC_count, MPC_max, MPC_dofs, ...
                                  MPC_coeff, S, C) ;
end % if MPC exist
 
% INFORCE 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) 
 
% ========================================================================
%                          PHASE 5
%         SOLVE FOR UNIQUE SOLUTION AND ITS REACTIONS
% ========================================================================
T = S \ C                 ; % Compute displacements & rotations 
list_save_beam_displacements (n_g, n_m, T)   ; % save and print 
 
show = 1 ;
if ( show==1)
%   OPTIONAL GRAPHS (saved as portable network graphics format)
  if ( BEF == 1 )                  % plot any foundation pressure 
    C1_L3_BoEF_pressure (F_k) % graph 
  pause (4) % display for 4 seconds 
  end % if 
  C1_L3_beam_plots(1)                                % deflection
  pause (4) % display for 4 seconds 
  C1_L3_beam_plots(2)                                     % slope      
  pause (2) % display for 2 seconds 
  C1_L3_beam_plots(3)                            % moment diagram
  pause (4) % display for 4 seconds 
  C1_L3_beam_plots(4)                  % transverse shear diagram
  pause (4) % display for 4 seconds 
end % if show
 
% =====================================================================
%                             PHASE 6
%       POST-PROCESS ELEMENTS (el qp values and/or el integral)
% =====================================================================
%    OPTIONAL FORCE AND MOMENT REACTION RECOVERY & SAVE 
if ( EBC_count > 0 )                    ; % reactions exist ? 
  EBC_react   = zeros (EBC_count, 1)    ; % pre-allocate
  [EBC_react] = recover_reactions_print_save (n_g, n_d, ... 
                EBC_flag, EBC_row, EBC_col, T); % reaction to EBC 
  % EBC_react = nodal source require to maintain EBC value
end % if EBC existed
 
%  OPTIONAL: POST-PROCESS ELEMENT REACTIONS (MEMBER FORCES & MOMENTS) 
% skip fprintf ('\nIndividual Element Load and Reaction Summaries: \n')
fprintf ('\nElement End Point Reaction Summaries: \n')
fprintf ('Elem, F_1,    M_1,    F_3,    M_3 \n')

%  Re-compute slement stiffnesses and loads, or read from storage
for j = 1:n_e ; % loop over elements ====>> ====>> ====>> ====>> ====>> 
  e_nodes = nodes (j, 1:n_n)                      ; % connectivity 
  [rows]  = get_element_index (n_g, n_n, e_nodes)   ; % eq numbers 

  T_e (1:n_i) = T(rows) ;           % gather element displacements 
   
%               SET ELEMENT PROPERTIES & GEOMETRY 
  L = abs (x(e_nodes(n_n)) - x(e_nodes(1)))        ; % beam length 
  Line_e (1:2) = 0 ;
  if ( n_prop > 1 ) 
    E            = msh_properties (j, 1)   ; % material modulus 
    I            = msh_properties (j, 2)   ; % section inertia 
    Line_e (1:2) = msh_properties (j, 3:4) ; % trapezoidal line load 
    k_f          = msh_properties (j, 5)   ; % foundation stiffness 
  end % if 
 
%   ELEMENT, FOUNDATION STIFFNESSES AND SOURCE RESULTANT MATRICES 
  L2 = L^2 ; L3 = L^3 ;  
%  Constant quintic element & foundation stiffnesses & mass matrices 
  S_e = [5092,   1138*L, -3584,    1920*L,  -1508,   242*L  ; 
         1138*L, 332*L2, -896*L,   320*L2,  -242*L,  38*L2  ; 
        -3584,   -896*L,  7168,    0,       -3584,   896*L  ; 
         1920*L, 320*L2,  0,       1280*L2, -1920*L, 320*L2 ; 
        -1508,  -242*L,  -3584,   -1920*L,   5092,  -1138*L ; 
         242*L,  38*L2,   896*L,   320*L2,  -1138*L, 332*L2 ] ... 
         *E*I/(35 * L^3) ;                            % stiffness 
 
  NtN = L*[2092,   114*L, 880,  -160*L,   262,   -29*L  ; 
           114*L,  8*L2,  88*L, -12*L2,   29*L,  -3*L2  ; 
           880,    88*L,  5632,  0,       880,   -88*L  ; 
          -160*L, -12*L2,  0,    128*L2,  160*L, -12*L2 ; 
           262,    29*L,  880,   160*L,   2092,  -114*L ; 
          -29*L,  -3*L2  -88*L, -12*L2,  -114*L,  8*L2 ]/13860;  
 
%  Map line load to node forces & moments; C_e = w_To_F * Line_e  
    C_e = zeros (n_i, 1)           ; % clear arrays 
    if ( any (Line_e) )            ; % then form forcing vector 
      w_To_F = L * [ 79,    19  ;  
                     5*L,   2*L ; 
                     112,   112 ; 
                    -8*L,   8*L ; 
                     19,    79  ; 
                    -2*L,  -5*L ] / 420 ; % quintic H, linear Line 
      % 6 x 2 * 2 x 1 = 6 x 1 result  
      C_e = w_To_F (1:n_i, 1:2) * Line_e' ; % force moment @ nodes 
    end % if or set up resultant node loads for line load 
 
%  Assign any point force/couple to the first element with that node 
  if ( j == 1 ) % first element only gets sources from 3 nodes 
    C_e = C_e + C_react (rows)   ;
  else % last 2 nodes 
    % quintic C1 beam ONLY !!
    C_e (3:6) = C_e (3:6) + C_react (rows(3:6)); % L3 C1 beam ONLY !! 
  end % if 
  %skip fprintf ('\nGiven Resultant Loading on Member %g, \n', j) 
  %skip disp (C_e') 
 
% Finally, get the individual elementreactions from all sources
  C_m = S_e * T_e' - C_e          ; % applied loading reactions
  if ( BEF )                        % add foundation resultants
    C_f = -k_f * NtN * T_e'  ;
    %skip fprintf ('Resultant Foundation Loading on Member %g, \n', j) 
    %skip disp (C_f') 
    C_m = C_m - C_f   ; % additional elastic foundation effects
  end % if
  %skip fprintf ('Net Resultant End Reactions on Member %g, \n', j) 
  %skip disp (C_m') 
  %b fprintf ('%i, %g, %g, %g, %g, %g, %g \n', j, C_m') % print them
  fprintf ('%i, %g, %g, %g, %g \n', j, C_m(1:2), C_m(5:6)') % print
end % for each j element in mesh   <<==== <<==== <<==== <<==== <<==== 
%       OPTIONAL POST_PROCESSING COMPLETED
 
% end function L3_Quintic_BoEF_MPC  % ==========================
   
%  +++++++++++++ functions in alphabetical order +++++++++++++++++ 
 
% See /home/mech517/public_html/Matlab_Plots for graphic options 
% http://www.owlnet.rice.edu/~mech517/help_plot.html for help 
 
function  C1_L3_beam_plots (i_p)  % ==========================
% Copyright 2008 J.E. Akin. All rights reserved.
%  ------------------------------------------------------
%  Matlab graph of i_p-th component value at mesh nodes
%  for a mesh of 3 node quintic Hermite line elements
%   i_p = 1 is deflection, = 2 is slope
%   = 3 is curvature (moment), = 4 is curvature rate (shear)
%  ------------------------------------------------------

% c_x    = x coordinates of nod_per_el line element
% msh_typ_nodes = connectivity list , nt x nod_per_el
% loop   = corners for nod_per_el line element
% nod_per_el = Nodes per element
% np     = Number of Points
% nt     = Number of elements
% pre_e  = Element items before connectivity list
  pre_e  = 0 ;
% pre_p  = Nodal items before coordinates
  pre_p = 1;
% msh_bc_xyz = Nodal coordinates (with preceeding data)
% t_x    = x coordinates of nod_per_el corners

  if ( nargin == 0 )
    i_p = 1 ; 
  end % if no arguments
  Pts_wide = 2 ; % draw fat lines
  format short

%    Read coordinate file and connectivity file
%    integer bc code, real xy pairs for np points (pre_p = 1)
  load msh_bc_xyz.tmp ;  

%    Set control data: number of points
  np = size (msh_bc_xyz,1) ;  % number of nodal points  
  %b fprintf ('Read %g mesh coordinates \n', np)
  ns = size (msh_bc_xyz,2) - pre_p ;   % space dimension
  if ( ns ~= 1)
    error ('This is not a 1D mesh')
  end % if not 1D data
  
%    Set control data: number elements
  load  msh_typ_nodes.tmp ; % nod_per_el nodes per element
  nt = size (msh_typ_nodes,1) ; % number of elements in mesh 
  nod_per_el = size (msh_typ_nodes,2) - pre_e -1 ;% nodes per elem
  %b fprintf ('Read %g elements connections \n', nt)
  if ( nod_per_el ~= 3 )
    error ('This is not a mesh of 3 node line elements')
  end % if 

  load node_results.tmp
  nr = size (node_results, 1);
  if ( nr == 0 )
    error ('Error missing file node_results.tmp')
  end % if error
  max_p = size (node_results, 2) ; % number of columns
  %b fprintf ('Read %g nodal solution values \n', nr)
  %b fprintf ('     with %g components each \n', max_p)
  if ( i_p > 4 ) 
    fprintf ('Data requested for component i_p = %g \n', i_p)
    error ('i_p > available data')
  end % if error

  H (3)       = 0. ;
  HC1 (6)     = 0. ; DHC1 (6)     = 0. ;
  D2HC1 (6)   = 0. ; D3HC1 (6)    = 0. ;
  D4HC1 (6)   = 0. ;
  x (np)      = 0.       ;      % pre-allocate array x
  y (np)      = 0.       ;      % pre-allocate array y
  dy (np)     = 0.       ;      % pre-allocate array y
  t_nodes (nod_per_el) = 0 ;    % Optional pre-allocation
  t_x (nod_per_el) = 0 ;        % Optional pre-allocation
  t_y (nod_per_el) = 0 ;        % Optional pre-allocation
  t_dy (nod_per_el) = 0 ;       % Optional pre-allocation
  c_x (nod_per_el) = 0 ;        % Optional pre-allocation
  c_y (nod_per_el) = 0 ;        % Optional pre-allocation

%                set constants
  loop = [1:nod_per_el] ; % default to sequential order

  % msh_bc_xyz has: pre_p items then: x, y
  x = msh_bc_xyz (1:np, (pre_p+1)) ; % extract x column
  xmax = max (x)      ; xmin = min (x) ;
  y  = node_results(:, 1) ;
  dy = node_results(:, 2) ;
  if ( any (y) == 0 )
    fprintf ('All values are zero. Plot skipped. \n')
    return
  end % if null solution

   clf                               % clear graphics
   hold on                           % hold image for plots
   xlabel (['X, Node at 45 deg (', int2str(nod_per_el), ...
     ' per element), Element at 90 deg'])

%    Are properties needed ?
  if ( i_p > 2 )
    load msh_properties.tmp
    n_prop = size (msh_properties, 2) ;
    n_mats = size (msh_properties, 1) ;
  else
    n_prop = 0;
    n_mats = 0;
    EI = 1 ;
  end % if

%    Loop over all elements
  max_all = -1e9  ; min_all = 1e9  ;
  for it = 1:nt  ;

%      Extract element connectivity 
    t_nodes = msh_typ_nodes (it, (pre_e+2):(nod_per_el+pre_e+1)); 

%      Skip point elements, if any
    if ( all (t_nodes) ) % then valid line

%      Extract element coordinates & values
      t_x  = x (t_nodes) ;         % x at those nodes, only
      if ( t_x(2) ~= (t_x(1) + t_x(3))/2 )
        fprintf ('WARNING: fixed bad mid-point in element %g \n', it)
        t_x(2) = (t_x(1) + t_x(3))/2 ;
      end % if non-constant Jacobian
      t_y  = y (t_nodes) ;         % y at those nodes, only
      t_dy  = dy (t_nodes) ;         % dy at those nodes, only
      D (1:2:6) = t_y ;
      D (2:2:6) = t_dy ;
      if ( i_p == 1 )
        plot (t_x, t_y, 'ko')        % plot nodal value symbols
      end % if

%      Loop over local points on the quadratic polynomial element
      n_poly = ceil ( 95 / nt) ;
      for k = 1: (n_poly + 1) % points in parametric space
%          get element parametric interpolation functions
        R     = (k - 1)/n_poly ;   % on 0 to 1
        X     = 2*R - 1 ;          % on -1 to 1
%       H = ELEMENT SHAPE FUNCTIONS
%       X = LOCAL COORDINATE OF POINT,    -1 TO +1
%       LOCAL NODE COORD. ARE -1,0,+1    1-----2-----3
        H (1) = 0.5*(X*X - X) ;
        H (2) = 1. - X*X ;
        H (3) = 0.5*(X*X + X) ;
        x_el (k) = H * t_x       ; % true x value
        A = abs(t_x(3) - t_x(1)) ; % length
        P = X ;
        P_2 = P * P   ; P_3 = P * P_2 ;
        P_4 = P * P_3 ; P_5 = P * P_4 ;

%       select result of interest, r_el
        if ( i_p == 1 ) % deflection
          HC1(1) = (4*P_2 - 5*P_3 - 2*P_4 + 3*P_5) * 0.25 ;
          HC1(2) = (P_2 - P_3 - P_4 + P_5) * 0.125 * A ;
          HC1(3) =  1 - 2*P_2 + P_4 ;
          HC1(4) = (P - 2*P_3 + P_5) * A * 0.5 ;
          HC1(5) = (4*P_2 + 5*P_3 - 2*P_4 - 3*P_5) * 0.25  ;
          HC1(6) = (-P_2 - P_3 + P_4 + P_5) * 0.125 * A ;
          r_el (k) = HC1 * D' ;  % true y value
  
        elseif ( i_p == 2 ) % slope      
          DHC1 (1) = (8*P - 15*P_2 - 8*P_3 + 15*P_4) * 0.5 / A ;
          DHC1 (2) = (2*P - 3*P_2 - 4*P_3 + 5*P_4) * 0.25 ;
          DHC1 (3) = (-4*P + 4*P_3) * 2 / A ;
          DHC1 (4) = (1 - 6*P_2 + 5*P_4) ;
          DHC1 (5) = (8*P + 15*P_2 - 8*P_3 - 15*P_4) * 0.5 / A ;
          DHC1 (6) = (-2*P - 3*P_2 + 4*P_3 + 5*P_4) * 0.25 ;
          r_el (k) = DHC1 * D' ;  % true y' value

        elseif ( i_p == 3 ) % curvature or moment
          D2HC1 (1) = (8 - 30*P - 24*P_2 + 60*P_3) /A/ A ;
          D2HC1 (2) = (2 - 6*P - 12*P_2 + 20*P_3) * 0.5/A ;
          D2HC1 (3) = (-4 + 12*P_2) * 4 / A / A ;
          D2HC1 (4) = (-12*P + 20*P_3)*2 / A ;
          D2HC1 (5) = (8 + 30*P - 24*P_2 - 60*P_3) / A /A;
          D2HC1 (6) = (-2 - 6*P + 12*P_2 + 20*P_3) * 0.5 / A ;
          r_el (k) = D2HC1 * D' ;  % true y'' value

        elseif ( i_p == 4 ) % curvature rate or shear force
          D3HC1 (1) = (-30 - 48*P + 180*P_2) * 2 /A/A/ A ;
          D3HC1 (2) = (-6 - 24*P + 60*P_2) /A /A;
          D3HC1 (3) = (24*P) * 8 / A / A/A ;
          D3HC1 (4) = (-12 + 60*P_2)*4 / A/A ;
          D3HC1 (5) = (30 - 48*P - 180*P_2) * 2 / A /A/A;
          D3HC1 (6) = (-6 + 24*P + 60*P_2) / A/A ;
          r_el (k) = D3HC1 * D' ;  % true y''' value
        end % if plot type

      end % for k plot points
      if ( n_mats > 1 )
        EI = msh_properties(it,1)*msh_properties(it,2)  ;
        r_el = r_el* EI ;
      elseif ( n_mats == 1 )
        EI = msh_properties( 1,1)*msh_properties( 1,2)  ;
        r_el = r_el* EI ;
      end % if

%    elem max, min values
      [V_X, L_X] = max (r_el)  ;
      [V_N, L_N] = min (r_el);
       if ( V_X > max_all )
         max_all = V_X ;
       end
       if ( V_N < min_all )
         min_all = V_N ;
       end
       plot (x_el, r_el, 'b-', 'LineWidth',Pts_wide)
    end % if zero in connectivity

% Plot the element number
    x_bar = sum (t_x' )/nod_per_el ;
    y_bar = mean (r_el) ;
    t_text = sprintf ('   (%g)', it);  % offset # from pt
    text (x_bar, y_bar, t_text) % incline

  end % for over all elements

  null (1:np) = 0.5*(V_N + V_X) ;
  if ( V_N <= 0 & V_X >= 0 )
    null (1:np) = 0. ;
  end % if
% finalize axes
   ymax = max_all ; ymin = min_all ;
   diff = abs(ymax-ymin) ;
   ymax = ymax + abs (diff)/10. ;
   ymin = ymin - abs (diff)/10. ;
   axis ([xmin, xmax, ymin, ymax])   % set axes

   if ( i_p == 1 )
     title(['Qunitic Beam Deflection from:  ', ...
            int2str(nt),' Elements, ', int2str(np),' Nodes'])     
     ylabel (['Deflection: max=', num2str(max_all), ...
               ', min=', num2str(min_all)])
     fprintf ('Max deflection value is %g \n', max_all)
     fprintf ('Min deflection value is %g \n', min_all)
   elseif ( i_p == 2 )
     title(['Qunitic Beam Slope from:  ', ...
            int2str(nt),' Elements, ', int2str(np),' Nodes'])     
     ylabel (['Slope: max=', num2str(max_all), ...
               ', min=', num2str(min_all)])       
     fprintf ('Max slope value is %g \n', max_all)
     fprintf ('Min slope value is %g \n', min_all)
   elseif ( i_p == 3 )
     title(['Qunitic Beam Moment from:  ', ...
            int2str(nt),' Elements, ', int2str(np),' Nodes'])     
     ylabel (['Moment: max=', num2str(max_all), ...
               ', min=', num2str(min_all)])   
     fprintf ('Max moment value is %g \n', max_all)
     fprintf ('Min moment value is %g \n', min_all)
   elseif ( i_p == 4 )
     title(['Qunitic Beam Shear from:  ', ...
            int2str(nt),' Elements, ', int2str(np),' Nodes'])     
     ylabel (['Transverse Shear: max=', num2str(max_all), ...
               ', min=', num2str(min_all)])   
     fprintf ('Max shear value is %g \n', max_all)
     fprintf ('Min shear value is %g \n', min_all)
   end % if

% plot node points on axis
  for i = 1:np
    t_text = sprintf ('   %g', i);  % offset # from pt
    text (x(i), null(i), t_text, 'Rotation', 45) % incline
  end % for all
  plot (x, null, 'k*')
  grid
 
%       Create plot copies for report use 
  if ( i_p == 1 )
    print ('-dpng', ['Beam_deflections'])
    v_text = ['Created Beam_deflections.png'] ;
  elseif ( i_p == 2 )
    print ('-dpng', ['Beam_slope'])
    v_text = ['Created Beam_slope.png'] ;
  elseif ( i_p == 3 )
    print ('-dpng', ['Beam_moments'])
    v_text = ['Created Beam_moments.png'] ;
  elseif ( i_p == 4 )
    print ('-dpng', ['Beam_shears'])
    v_text = ['Created Beam_shears.png'] ;
  end % if
  fprintf (1,'%s', v_text) ; fprintf (1, ' \n' )
  hold off
% end of C1_L3_beam_plots  % ==========================

function  C1_L3_BoEF_pressure (F_k)  % ==========================
% Copyright 2008 J.E. Akin. All rights reserved.
%  ------------------------------------------------------
%  Matlab graph of foundation pressurelue at mesh nodes
%  for a mesh of 3 node quintic Hermite line elements
%  ------------------------------------------------------
% c_x    = x coordinates of nod_per_el line element
% msh_typ_nodes = connectivity list , nt x nod_per_el
% loop   = corners for nod_per_el line element
% nod_per_el = Nodes per element
% np     = Number of Points
% nt     = Number of elements
% pre_e  = Element items before connectivity list
  pre_e  = 0 ;
% pre_p  = Nodal items before coordinates
  pre_p = 1;
% msh_bc_xyz = Nodal coordinates (with preceeding data)
% t_x    = x coordinates of nod_per_el corners

  if ( nargin == 0 )
    i_p = 1 ; 
  end % if no arguments
  n_fk = size (F_k, 2) ; % number of different foundations
  Pts_wide = 2 ;  % fat lines
  format short

%    Read coordinate file and connectivity file
%    integer bc code, real xy pairs for np points (pre_p = 1)
  load msh_bc_xyz.tmp ;  

%    Set control data: number of points
  np = size (msh_bc_xyz,1) ;  % number of nodal points  
  %b fprintf ('Read %g mesh coordinates \n', np)
  ns = size (msh_bc_xyz,2) - pre_p ;   % space dimension
  if ( ns ~= 1)
    error ('This is not a 1D mesh')
  end % if not 1D data
  
%    Set control data: number elements
  load  msh_typ_nodes.tmp ; % nod_per_el nodes per element
  nt = size (msh_typ_nodes,1) ; % number of elements in mesh 
  nod_per_el = size (msh_typ_nodes,2) - pre_e -1 ;% nodes per elem
  %b fprintf ('Read %g elements connections \n', nt)
  if ( nod_per_el ~= 3 )
    error ('This is not a mesh of 3 node line elements')
  end % if 

  load node_results.tmp
  nr = size (node_results, 1);
  if ( nr == 0 )
    error ('Error missing file node_results.tmp')
  end % if error
  max_p = size (node_results, 2) ; % number of columns
  %b fprintf ('Read %g nodal solution values \n', nr)
  %b fprintf ('     with %g components each \n', max_p)

  H (3)       = 0. ;
  HC1 (6)     = 0. ; DHC1 (6)     = 0. ;
  D2HC1 (6)   = 0. ; D3HC1 (6)    = 0. ;
  D4HC1 (6)   = 0. ;
  x (np)      = 0.       ;      % pre-allocate array x
  y (np)      = 0.       ;      % pre-allocate array y
  dy (np)     = 0.       ;      % pre-allocate array y
  t_nodes (nod_per_el) = 0 ;    % Optional pre-allocation
  t_x (nod_per_el)     = 0 ;    % Optional pre-allocation
  t_y (nod_per_el)     = 0 ;    % Optional pre-allocation
  t_dy (nod_per_el)    = 0 ;    % Optional pre-allocation
  c_x (nod_per_el)     = 0 ;    % Optional pre-allocation
  c_y (nod_per_el)     = 0 ;    % Optional pre-allocation

%                set constants
  loop = [1:nod_per_el] ; % default to sequential order

  % msh_bc_xyz has: pre_p items then: x, y
  x = msh_bc_xyz (1:np, (pre_p+1)) ; % extract x column
  xmax = max (x)      ; xmin = min (x) ;
  y  = node_results (:, 1) ;
  dy = node_results (:, 2) ;

   clf                               % clear graphics
   hold on                           % hold image for plots
   xlabel (['X, Node at 45 deg (', int2str(nod_per_el), ...
     ' per element), Element at 90 deg'])

%    Are properties needed ?
    load msh_properties.tmp
    n_prop = size (msh_properties, 2) ;
  
%    Loop over all elements
  max_all = -1e9  ; min_all = 1e9  ;
  for it = 1:nt  ;
%      Element foundation stiffness, if any
    if ( n_fk == 1 ) % then homogeneous
      k_f = F_k (1) ;
    else
      k_f = F_k (it) ;
    end % if 

%      Extract element connectivity 
    t_nodes = msh_typ_nodes (it, (pre_e+2):(nod_per_el+pre_e+1)); 

%      Skip point elements, if any
    if ( all (t_nodes) ) % then valid line

%      Extract element coordinates & values
      t_x  = x (t_nodes) ;         % x at those nodes, only
      if ( t_x(2) ~= (t_x(1) + t_x(3))/2 )
        fprintf ('WARNING: fixed bad mid-point in element %g \n', it)
        t_x(2) = (t_x(1) + t_x(3))/2 ;
      end % if non-constant Jacobian
      t_y  =  y (t_nodes) ;         % y at those nodes, only
      t_dy = dy (t_nodes) ;         % dy at those nodes, only
      D (1:2:6) = t_y ; D (2:2:6) = t_dy ;

%      Loop over local points on the quadratic polynomial element
      n_poly = ceil ( 95 / nt) ;
      for k = 1: (n_poly + 1) % points in parametric space

%          get element geometry parametric interpolation functions
        R = (k - 1)/n_poly ;   % on 0 to 1
        X = 2*R - 1 ;          % on -1 to 1
%       H = ELEMENT SHAPE FUNCTIONS
%       X = LOCAL COORDINATE OF POINT,    -1 TO +1
%       LOCAL NODE COORD. ARE -1,0,+1    1-----2-----3
        H (1) = 0.5*(X*X - X) ;
        H (2) = 1. - X*X ;
        H (3) = 0.5*(X*X + X) ;
        x_el (k) = H * t_x ;  % true x value
        A = abs (t_x(3) - t_x(1)) ; % length
        P = X ;
        P_2 = P * P   ; P_3 = P * P_2 ;
        P_4 = P * P_3 ; P_5 = P * P_4 ;

%          get element deflections parametric interpolation functions
        HC1(1) = (4*P_2 - 5*P_3 - 2*P_4 + 3*P_5) * 0.25 ;
        HC1(2) = (P_2 - P_3 - P_4 + P_5) * 0.125 * A ;
        HC1(3) =  1 - 2*P_2 + P_4 ;
        HC1(4) = (P - 2*P_3 + P_5) * A * 0.5 ;
        HC1(5) = (4*P_2 + 5*P_3 - 2*P_4 - 3*P_5) * 0.25  ;
        HC1(6) = (-P_2 - P_3 + P_4 + P_5) * 0.125 * A ;
        r_el (k) = HC1 * D' ;  % true y value
  
      end % for k plot points
      
      r_el = -r_el* k_f ; % deflection * foundation stiffness

%    elem max, min values
      [V_X, L_X] = max (r_el)  ;
      [V_N, L_N] = min (r_el);
       if ( V_X > max_all )
         max_all = V_X ;
       end
       if ( V_N < min_all )
         min_all = V_N ;
       end
       plot (x_el, r_el, 'b-', 'LineWidth',Pts_wide)
    end % if zero in connectivity

% Plot the element number
    x_bar = sum (t_x' )/nod_per_el ;
    y_bar = mean (r_el) ;
    t_text = sprintf ('   (%g)', it);  % offset # from pt
    text (x_bar, y_bar, t_text) % incline

  end % for over all elements

  fprintf ('Max foundation pressure value is %g \n', max_all)
  fprintf ('Min foundation pressure value is %g \n', min_all)
  null (1:np) = 0.5*(V_N + V_X) ;
  if ( V_N <= 0 & V_X >= 0 )
    null (1:np) = 0. ;
  end % if
% finalize axes
   ymax = max_all ; ymin = min_all ;
   diff = abs(ymax-ymin) ;
   ymax = ymax + abs (diff)/10. ;
   ymin = ymin - abs (diff)/10. ;
   axis ([xmin, xmax, ymin, ymax])   % set axes

   title(['Beam Foundation Pressure from:  ', ...
          int2str(nt),' Elements, ', int2str(np),' Nodes'])     
   ylabel (['Pressure: max=', num2str(max_all), ...
             ', min=', num2str(min_all)])

% plot node points on axis
  for i = 1:np
    t_text = sprintf ('   %g', i);  % offset # from pt
    text (x(i), null(i), t_text, 'Rotation', 45) % incline
  end % for all
  plot (x, null, 'k*')
  grid
 
% create plot copies for report use 
  print ('-dpng', ['BoEF_pressure'])
  v_text = ['Created BoEF_pressure.png'] ;
  fprintf (1,'%s', v_text) ; fprintf (1, ' \n' )
  hold off
% end of C1_L3_BoEF_pressure  % ==========================
 
function echo_files_info  % ==========================
file_id = fopen('msh_files_info.tmp', 'r'); % open for read only
if ( file_id > 0 )                 % then file exists
  fprintf('\n(Echo of msh_files_info.tmp) \n')
  while feof(file_id) ==0          % read until an end of file
    remark_line = fgets(file_id) ; % read a string and append \n
    fprintf (remark_line)        ; % echo the line
  end % while data exists
end % if file supplied
fprintf('\n')
% end function echo_file_info  % ==========================

function echo_user_remarks  % ==========================
file_id = fopen('msh_remarks.tmp', 'r'); % open for read only
if ( file_id > 0 )                 % then file exists
  fprintf('\n(Echo of msh_remarks.tmp) \n')
  while feof(file_id) ==0          % read until an end of file
    remark_line = fgets(file_id) ; % read a string and append \n
    fprintf (remark_line)        ; % echo the line
  end % while data exists
end % if file supplied
fprintf('\n')
% end function echo_user_remarks  % ==========================

function [S, C] = enforce_essential_BC (EBC_flag, EBC_value, S, C) 
% modify system linear eqs for essential boundary conditions 
% (by trick to avoid matrix partitions, loses reaction data) 
 n_d = size (C, 1)                        ; % number of DOF eqs 
if ( size (EBC_flag, 2) > 1 ) ; % change to vector copy 
  flag_EBC  = reshape ( EBC_flag',  1, n_d) ; 
  value_EBC = reshape ( EBC_value', 1, n_d) ; 
else 
  flag_EBC  = EBC_flag ; 
  value_EBC = EBC_value ; 
end % if 
 for j = 1:n_d                % check all DOF for essential BC  
   if ( flag_EBC (j) == 1 )                   % then EBC here 
%    Carry known columns*EBC to RHS. Zero that column and row. 
%    Insert EBC identity, 1*EBC_dof = EBC_value. 
     EBC = value_EBC (j)             ; % recover EBC value 
     C (:) = C (:) - EBC * S (:, j) ; % carry known column to RHS 
     S (:, j) = 0 ; S (j, :) = 0    ; % clear, restore symmetry 
     S (j, j) = 1 ; C (j) = EBC     ; % insert identity into row 
   end % if EBC for this DOF 
 end % for over all j-th DOF 
% end enforce_essential_BC    % ==========================

  function [S, C] = enforce_MPC_equations (MPC_count, MPC_max, ...
                                           MPC_dofs, MPC_coeff, S, C)
% Inforce algebraic constraint equations (MPC), by penalty method
% MPC_count = number of multipoint constraint equations
% MPC_dofs  = the dof numbers in each constraint equation
% MPC_coeff = the dof coefficients in each constraint equation
% T(MPC_dofs(1,j)) + T(MPC_dofs(2,j))*MPC_coeff(1,j) = MPC_coeff(2,j)
  S_P   = zeros (MPC_max, MPC_max)     ;    % penalty square matrix
  C_P   = zeros (MPC_max, 1)           ;    % penalty column matrix
  for m = 1:MPC_count                  ; % loop over each contraint
    n_dof = sum ( MPC_dofs (:, m) > 0 )         ; % constraint size
    rows  = MPC_dofs  (1:n_dof, m)            ; % vector subscripts
    coeff = MPC_coeff (1:n_dof, m)                  ; % work vector
    scale = S (rows(1), rows(1)) * 500              ;  % MPC weight
    V (1:n_dof+1) = [1, coeff']                     ; % work vector
    VtV   = V' * V                                  ;  % work array
    C_P (1:n_dof) = -VtV (1:n_dof, n_dof+1)         ;  % MPC vector
    S_P (1:n_dof, 1:n_dof) = VtV (1:n_dof, 1:n_dof) ;  % MPC square

  % rows = vector subscript converting element to system eq numbers
    S (rows, rows) = S (rows, rows) + S_P(:, :) * scale     ; % add 
    C (rows)       = C (rows)       + C_P(:)    * scale     ; % add
  end % for MPC with m dof entries
% end enforce_MPC_equations (can be used for EBC with bigger scale)
 
function [C] = get_and_add_point_sources (n_g, n_m, C) 
  load msh_load_pt.tmp            ; % node, DOF, value (eq. number) 
  n_u = size(msh_load_pt, 1)      ; % number of point sources 
  if ( n_u < 1 )                  ; % missing data 
    error ('No load_pt data in msh_load_pt.tmp') 
  end % if user error 
  fprintf ('\nRead %i point sources. (Echo of msh_load_pt.tmp) \n', n_u) 
  fprintf ('Node, DOF,  Source_value \n') 
  for j = 1:n_u                    ; % non-zero Neumann pts 
    node  = msh_load_pt (j, 1)         ; % global node number 
    DOF   = msh_load_pt (j, 2)         ; % local DOF number 
    value = msh_load_pt (j, 3)         ; % point source value 
    fprintf ('%i  %i  %g \n', node, DOF, value) 
    Eq    = n_g * (node - 1) + DOF     ; % row in system matrix 
    C (Eq) = C (Eq) + value     ; % add to system column matrix 
  end % for each EBC 
% end get_and_add_point_sources    % ==========================
 
function [MPC_coeff, MPC_dofs] = get_constraint_eqs (n_g, ...
                                 MPC_each, MPC_max, MPC_count);
% MPC_count = number of multipoint constraint equations
% MPC_coeff = the dof coefficients in each constraint equation
% MPC_dofs  = the dof numbers in each constraint equation
% MPC_max   = maximum number of dof in a single constraint, >=2
% n_g       = number of DOF per node
  load msh_mpc.tmp     ; % n_1, d_1, n_2, d_2, .., d_n, c_1, ..c_n
  n_c = size (msh_mpc, 1)       ; % number of constraint equations
  if ( n_c < 1 )                         ; % missing expected data 
    error ('No algebraic constraint data in msh_mpc.tmp') 
  end % if user error 
  fprintf ('\nRead %i constraint eqs. (Echo of msh_mpc.tmp)\n', n_c)
  fprintf ('%i Node-DOF pairs, %i coefficients\n', MPC_max, MPC_max)
  disp (msh_mpc)                                  % echo input file
  Eq = 0 ;                                % initialize MPC eq count
  for n = 2:MPC_max          ; % loop over each MPC constraint type
    exist = MPC_each (n)     ;     % number of lines of type n data
    if ( exist > 0 )         ;         % expect rows of type n data
      for k = 1:exist        ;             % read each line of data
        Eq = Eq + 1          ;             % increment MPC eq count
        % extract data line for this MPC equation
        line (1:3*n) = msh_mpc (Eq, 1:3*n)     ;    % user MPC data

        % convert constraints to system equation numbers & store
        c_node (1:n) = line (1:2:2*n)          ; % pair node number
        c_dof  (1:n) = line (2:2:2*n)          ; % pair dof  number
        c_coef (1:n) = line (2*n+1:3*n)        ; % pair coefficents
        MPC_coeff (1:n, Eq) = c_coef (1:n)'               ; % store
        MPC_dofs  (1:n, Eq) = n_g * (c_node' -1) + c_dof' ; % store
      end % for k MPC data lines
    end % if data supplied for MPC type n
  end % for n                            % for each constraint type
  
  if ( n_c ~= MPC_count )
    fprintf ('Compare MPC flags in msh_bc_xyz.tmp \n')
    fprintf (' to MPC equations in msh_mpc.tmp \n')
    error ('Inconsistent MPC flag count and supplied equations' )
  end % if inconsistent data
% end get_constraint_eqs  =====================================

function [EBC_flag] = get_ebc_flags (n_g, n_m, P) % ===========
 EBC_flag = zeros (n_m, n_g)   ; % initialize 
 for k = 1:n_m                        ; % loop over all nodes 
   if ( P(k) > 0 )                    ; % at least one EBC here 
     [flags] = unpack_pt_flags (n_g, k, P(k)) ; % unpacking 
     EBC_flag (k, 1:n_g) = flags (1:n_g)      ; % populate array 
   end % if EBC at node k 
 end % for loop over all nodes 
% end get_ebc_flags   % ==========================
 
function [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag) % ===
 EBC_value = zeros (n_m, n_g) ; % initialize to zero 
  load msh_ebc.tmp            ; % node, DOF, value (eq. number) 
  n_c = size(msh_ebc, 1)      ; % number of constraints 
  fprintf ('Read %i EBC data sets with Node, DOF, Value. \n', n_c) 
  fprintf ('(Echo of file msh_ebc.tmp) \n')  
  disp(msh_ebc)                    ; % echo input 
 
  for j = 1:n_c                    ; % loop over ebc inputs 
    node  = round (msh_ebc (j, 1)) ; % node in mesh 
    DOF   = round (msh_ebc (j, 2)) ; % DOF # at node 
    value = msh_ebc (j, 3)         ; % EBC value    
    % Eq  = n_g * (node - 1) + DOF ; % row in system matrix 
    EBC_value (node, DOF) = value  ; % insert value in array 
    if ( EBC_flag (node, DOF) == 0 ) % check data consistency 
      fprintf ('WARNING: EBC but no flag at node %i & DOF %i. \n', ... 
                node, DOF) 
      % EBC_flag (node, DOF) = 1; % try to recover from data error 
    end % if common user error 
  end % for each EBC 
 
  EBC_count = sum (sum ( EBC_flag == 1 ))   ; % check input data 
  if ( EBC_count ~= n_c )                ; % probable user error 
    fprintf ('WARNING: mismatch in bc_flag count & msh_ebc.tmp') 
  end % if user error 
% end get_ebc_values   % ==========================
 
function [rows] = get_element_index (n_g, n_n, e_nodes) % ====
% calculate system DOF numbers of element, for gather, scatter 
rows = zeros (1, n_g*n_n)           ; % allow for node = 0 
for k = 1:n_n                       ; % loop over element nodes 
  global_node = round (e_nodes (k)) ; % corresponding sys node 
  for i = 1:n_g                     ; % loop over DOF at node 
    eq_global  = i + n_g * (global_node - 1) ; % sys DOF, if any 
    eq_element = i + n_g * (k           - 1) ; % el DOF number 
    if ( eq_global > 0 )                ; % check node=0 trick 
      rows (1, eq_element) = eq_global  ; % valid DOF > 0 
    end % if allow for omitted nodes 
  end % for DOF i                       % end local DOF loop 
end % for each element node             % end local node loop 
% end get_element_index    % ==========================
 
function [n_e, n_n, n_t, el_type, nodes] = get_mesh_elements (); 
% input file controls (for various data generators)  
 
load  msh_typ_nodes.tmp        ; % el_type, connectivity list (3)  
n_e = size (msh_typ_nodes,1)   ; % number of elements  
if ( n_e == 0 )                ; % data file missing  
  error ('Error missing file msh_typ_nodes.tmp')  
end % if error  
n_n = size (msh_typ_nodes,2) - 1 ; % nodes per element        
fprintf ('\n')
fprintf ('Read %i elements with (ignored) type & %i nodes each. \n', ...
          n_e, n_n)  
fprintf ('(Echo number, row of file msh_typ_nodes.tmp) \n')
el_type = round (msh_typ_nodes(:, 1)); % el type number >= 1  
n_t     = max(el_type)                   ; % number of element types  
nodes (1:n_e, 1:n_n) = msh_typ_nodes (1:n_e, 2:1+n_n);  
%b disp(msh_typ_nodes (:, 1:1+n_n)) % echo data 
for i = 1:n_e
  fprintf ('%i, %i %i %i %i \n', i, msh_typ_nodes (i, :))
end % for each element
% end get_mesh_elements    % ==========================
 
function [n_m, n_s, P, x, y, z] = get_mesh_nodes (); % =======
% input file controls (for various data generators)  
  
%       READ MESH AND EBC_FLAG INPUT DATA  
% specific problem data from MODEL data files (sequential)  
load msh_bc_xyz.tmp       ; % bc_flag, x-, y-, z-coords  
n_m = size (msh_bc_xyz,1) ; % number of nodal points in mesh  
if ( n_m == 0 )           ; % data missing !  
  error ('Error missing file msh_bc_xyz.tmp')  
end % if error  
n_s = size (msh_bc_xyz,2) - 1 ; % number of space dimensions  
fprintf ('Read %i nodes. \n', n_m)  
fprintf ('(Echo node, row of file msh_bc_xyz.tmp) \n')
fprintf ('node, bc_flags, %i coordinates \n', n_s)
%b disp(msh_bc_xyz)
msh_bc_xyz (:, 1)= round (msh_bc_xyz (:, 1)); 
P = msh_bc_xyz (1:n_m, 1) ; % integer Packed BC flag  
x = msh_bc_xyz (1:n_m, 2) ; % extract x column   
y (1:n_m, 1) = 0.0 ;  z (1:n_m, 1) = 0.0 ; % default to zero 
  
if (n_s > 1 )                                 ; % check 2D or 3D 
  y = msh_bc_xyz (1:n_m, 3)           ; % extract y column   
end % if 2D or 3D 
if ( n_s == 3 )                               ; % check 3D 
  z = msh_bc_xyz (1:n_m, 4)           ; % extract z column  
end % if 3D  
for j = 1:n_m
  fprintf ('%i, %2.2i  %g         \n', j, P(j), x(j)       )
%bfprintf (' %2.2i  %g  %g     \n', P(j), x(j), y(j)       )
%bfprintf (' %2.2i  %g  %g  %g \n', P(j), x(j), y(j), z(j) )
end % for
% end get_mesh_nodes    % ==========================
 
function list_save_beam_displacements (n_g, n_m, T) % ==============
 fprintf('\nCalculated Displacements \n')
 fprintf('Node,  Y_displacement,  Z_rotation at %g nodes \n', n_m) 
 T_matrix = reshape (T, n_g, n_m)'            ; % pretty shape        
 % save results (displacements) to MODEL file: node_results.tmp 
 fid = fopen('node_results.tmp', 'w') ; % open for writing 
 for j = 1:n_m                        ; % node loop, save displ 
   fprintf (fid, '%g %g \n',  T_matrix (j, 1:n_g)) ; % to file 
   fprintf ('%g %g %g \n', j, T_matrix (j, 1:n_g)) ; % to screen 
 end % for j DOF 
% end list_save_beam_displacements     % ==========================
 
function list_save_displacements_results (n_g, n_m, T) % ============
  fprintf ('\n') ; 
  fprintf('X_disp      Y_disp      Z_disp   at %g nodes \n', n_m) 
  T_matrix = reshape (T, n_g, n_m)'            ; % pretty shape        
  disp (T_matrix)                              ; % print displacements 
  % save results (displacements) to MODEL file: node_results.tmp 
  fid = fopen('node_results.tmp', 'w')          ; % open for writing 
  for j = 1:n_m                                 ; % save displacements 
    if ( n_g == 1 ) 
      fprintf (fid, '%g                \n', T_matrix (j, 1:n_g)) ; 
    elseif ( n_g == 2 ) 
      fprintf (fid, '%g %g             \n', T_matrix (j, 1:n_g)) ; 
    elseif ( n_g == 3 ) 
      fprintf (fid, '%g %g %g          \n', T_matrix (j, 1:n_g)) ; 
    elseif ( n_g == 4 ) 
      fprintf (fid, '%g %g %g %g       \n', T_matrix (j, 1:n_g)) ; 
    elseif ( n_g == 5 ) 
      fprintf (fid, '%g %g %g %g %g    \n', T_matrix (j, 1:n_g)) ; 
    elseif ( n_g == 6 ) 
      fprintf (fid, '%g %g %g %g %g %g \n', T_matrix (j, 1:n_g)) ; 
    else  
      error ('reformat list_save_displacements_results for n_g > 6.') 
    end % if  
  end % for j DOF 
% end list_save_displacements_results    % ==========================
 
function [EBC_react] = recover_reactions_print_save (n_g, n_d, ... 
                       EBC_flag, EBC_row, EBC_col, T) 
% get EBC reaction values by using rows of S & C (before EBC) 
  n_d = size (T, 1)                     ; % number of system DOF 
% n_c x 1 = n_c x n_d  * n_d x 1 +  n_c x 1 
  EBC_react = EBC_row * T - EBC_col     ; % matrix reactions (+-) 
  % save reactions (forces) to MODEL file: node_reaction.tmp 
  fprintf ('\nReactions at Displacement BCs \n')         
  fprintf ('Node, DOF, Reaction Value \n')     
  fid = fopen('node_reaction.tmp', 'w') ; % open for writing 
  if ( size (EBC_flag, 2) > 1 )            ; % change to vector copy 
    flag_EBC  = reshape ( EBC_flag',  1, n_d) ; % changed 
  else 
    flag_EBC  = EBC_flag                      ; % original vector 
  end % if 
  Totals = zeros (1, n_g)               ; % zero input totals 
  kount = 0                             ; % initialize counter 
  for j = 1:n_d                         ; % extract all EBC reactions 
    if ( flag_EBC(j) )                      ; % then EBC here 
  %   Output node_number, component_number, value, equation_number  
      kount = kount + 1                     ; % copy counter 
      node  = ceil(j/n_g)                   ; % node at DOF j 
      j_g   = j - (node - 1)*n_g            ; % 1 <= j_g <= n_g 
      React = EBC_react (kount, 1)          ; % reaction value 
      fprintf ( fid, '%g %g %g \n', node, j_g, React);% save 
      fprintf ('%g %g %g \n', node, j_g, React);     % print 
      Totals (j_g) = Totals (j_g) + React   ; % sum all components 
    end % if EBC for this DOF 
  end % for over all j-th DOF 
  fprintf ('Totals = ') ; disp(Totals)      ; % echo totals 
% end recover_reactions_print_save    % ==========================
 
function [EBC_row, EBC_col] = save_reaction_matrices (EBC_flag, S, C) 
  n_d       = size (C, 1)                  ; % number of system DOF 
  EBC_count = sum (sum (EBC_flag == 1))    ; % count EBC & reactions 
  EBC_row   = zeros (EBC_count, n_d)       ; % reaction data 
  EBC_col   = zeros (EBC_count, 1)         ; % reaction data 
  if ( size (EBC_flag, 2) > 1 )            ; % change to vector copy 
    flag_EBC  = reshape ( EBC_flag',  1, n_d) ; % changed 
  else 
    flag_EBC  = EBC_flag                      ; % original vector 
  end % if 
  kount = 0                                ; % initialize counter 
  for j = 1:n_d      % System DOF loop, check for displacement BC 
    if ( flag_EBC (j) )                    ; % then EBC here 
  %   Save reaction data to be destroyed by EBC solver trick 
      kount = kount + 1                    ; % copy counter 
      EBC_row(kount, 1:n_d) = S (j, 1:n_d) ; % copy reaction data 
      EBC_col(kount, 1)     = C (j)        ; % copy reaction data 
    end % if EBC for this DOF 
  end % for over all j-th DOF                % end sys DOF loop 
% end save_reaction_matrices    % ==========================
 
function save_resultant_load_vectors (n_g, C) % ================
% save resultant forces to MODEL file: node_resultants.tmp 
  n_d = size (C, 1)                     ; % number of system DOF 
  fprintf ('\n')                        ; % Skip a line 
  % fprintf ('Node, DOF, Resultant Force (1) or Moment (2)  \n') 
  fprintf ('Resultant input sources: \n') 
  fprintf ('Node, DOF, Resultant input sources \n') 
  fid = fopen('node_resultant.tmp', 'w'); % open for writing 
  Totals = zeros (1, n_g)               ; % zero input totals 
  for j = 1:n_d                         ; % extract all resultants 
    if ( C (j) ~= 0. )                      ; % then source here 
  %   Output node_number, component_number, value, equation_number  
      node  = ceil(j/n_g)                   ; % node at DOF j 
      j_g   = j - (node - 1)*n_g            ; % 1 <= j_g <= n_g 
      value = C (j)                         ; % resultant value 
      fprintf ( fid, '%g %g %g %g \n', node, j_g, value, j);% save 
      fprintf ('%g %g %g \n', node, j_g, value);     % print 
      Totals (j_g) = Totals (j_g) + value   ; % sum all inputs 
    end % if non-zero for this DOF 
  end % for over all j-th DOF 
  fprintf ('Totals = ') ; disp(Totals)      ; % echo totals 
% end save_resultant_load_vectors    % ==========================
 
function [flags] = unpack_pt_flags (n_g, N, flag) % ============
% unpack n_g integer flags from the n_g digit flag at node N 
% integer flag contains (left to right) f_1 f_2 ... f_n_g 
 full = flag                         ; % copy integer 
 check = 0                           ; % validate input 
 for Left2Right = 1:n_g              ; % loop over local DOF at k 
   Right2Left = n_g + 1 - Left2Right ; % reverse direction 
   work = floor (full / 10)          ; % work item 
   keep = full - work * 10           ; % work item 
   flags (Right2Left) = keep         ; % insert into array 
   full = work                       ; % work item 
   check = check + keep * 10^(Left2Right - 1) ; % validate 
 end % for each local DOF 
 if ( flag > check )                 ; % check for likely error 
   fprintf ('WARNING: bc flag likely reversed at node %g. \n', N) 
 end % if likely user error 
% end unpack_pt_flags   % ==========================
  
%=================== Running gives =======================