function Modular_2D_Frame (load_pt)
% load_pt = 1 flags that point forces and/or couples will be 
%  input via file msh_load_pt.tmp
%.............................................................
%     Classic cubic beam plane frame for 
%      point loads & couples, line load
%.............................................................
if ( nargin == 0 )                 ; % check for optional data
  load_pt = 0                         ; % no point source data
end % if from argument count

%         Application and element dependent controls
n_g = 3  ; % number of DOF per node (axial_d, transverse_d, rot)
n_q = 0  ; % number of quadrature points required
n_r = 1  ; % number of rows in B_e matrix

%              Read mesh input data files
[n_m, n_s, P, x, y, z] = get_mesh_nodes     ;

%      Extract EBC flags from packed integer flag P 
[EBC_flag] = get_ebc_flags (n_g, n_m, P)      ; % unpack flags
EBC_count  = sum( sum ( EBC_flag > 0 ) )      ; % # of EBC
 
[n_e, n_n, n_t, el_type, nodes] = get_mesh_elements  ;
n_d = n_g*n_m                ; % system degrees of freedom (DOF)
n_i = n_g*n_n                ; % number of DOF per element
S = zeros (n_d, n_d) ; C = zeros (n_d, 1)      ; % initalize sums
M = zeros (n_d, n_d) ;                         ; % initalize sums

%                  Read EBC values, if any
if ( EBC_count > 0 )                            ; % need EBC data
  [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag) ; % read data
end % if any EBC data expected

%    Read point loads or moments, if any, and insert in C
if ( load_pt > 0 )                       ; % need point loads data
  [C] = get_and_add_point_sources (n_g, n_m, C); % add point loads
end % if any point source expected

%    Load the element properties array
load msh_properties.tmp ; % one row per element
n_prop = size(msh_properties, 1) ; % == 1 if homogeneous
fprintf ('\n(Echoing file msh_properties.tmp) \n')
Line_e (1:2) = 0 ;
if ( n_prop == 1 )
  A            = msh_properties (1, 1)   ; % area
  E            = msh_properties (1, 2)   ; % material modulus
  I            = msh_properties (1, 3)   ; % section inertia
  Line_e (1:2) = msh_properties (1, 4:5) ; % trapezoidal line load
  Rho          = msh_properties (1, 6)   ; % beam mass per unit length

%                   ECHO PROPERTIES
  fprintf ('Homogeneous Element Properties: \n' )
  fprintf ('Area (m^2)                  = %g \n', A)     
  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 ('Mass per unit length (kg/m) = %g  \n', Rho)
end % if
if ( n_prop > n_e )
  error ('ERROR: number of property sets exceeds number of elements')
end % if

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

for j = 1:n_e            ; % loop over elements ====>> ====>>
  S_e = zeros (n_i, n_i) ; M_e = zeros (n_i, n_i) ; % sys arrays
  C_p = zeros (n_i, 1)   ; C_e = zeros (n_i, 1)   ; % sys arrays
  s_L = zeros (n_i, n_i) ; m_L = zeros (n_i, n_i) ; % loc arrays
  t_L = zeros (n_i, n_i) ; c_L = zeros (n_i, 1)   ; % loc arrays
  Line_e (1:2) = 0 ;
  e_nodes = nodes (j, 1:n_n)                      ; % connectivity
  
%               SET ELEMENT PROPERTIES & GEOMETRY
  if ( n_prop > 1 )
    A            = msh_properties (j, 1)   ; % area
    E            = msh_properties (j, 2)   ; % material modulus
    I            = msh_properties (j, 3)   ; % section inertia
    Line_e (1:2) = msh_properties (j, 4:5) ; % trapezoidal line load
    Rho          = msh_properties (j, 6)   ; % beam mass per unit length

%                   ECHO PROPERTIES
    fprintf ('\nProperties for element %g \n', j)
    fprintf ('Area (m^2)                  = %g \n', A)     
    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 ('Mass per unit length (kg/m) = %g  \n', Rho)
  end % if

%-->   find member length and direction cosines
  dx  = x(e_nodes(2)) - x(e_nodes(1)) ; % x length
  dy  = y(e_nodes(2)) - y(e_nodes(1)) ; % y length
  L_e = sqrt (dx * dx + dy * dy)      ; % total length
  cx  = dx / L_e ; cy   = dy / L_e    ; % direction cosines
  IbL = I / L_e  ; IbL2 = I / L_e^2   ; % bending constants

%      ELEMENT CONDUCTION AND INTERNAL SOURCE MATRICES 

%   Linear axial bar and cubic bending. DOF = u, v, r, u, v, r
%   Form arrays in local axes, transform.     1  2  3  4  5  6

%                 stiffness
s_L = [ A,    0,       0,    -A,    0,       0    ;  
        0,  12*IbL2,  6*IbL,  0, -12*IbL2,  6*IbL ; 
        0,   6*IbL,   4*I,    0,  -6*IbL,   2*I   ; 
       -A,    0,       0,     A,    0,       0    ;  
        0, -12*IbL2, -6*IbL,  0,  12*IbL2, -6*IbL ; 
        0,   6*IbL,   2*I,    0,  -6*IbL,   4*I   ] * E / L_e 
disp(s_L/250000)
%  Map line load to node forces & moments; c_L = p_2_F * Line_e 
  if ( any (Line_e) )          ; % then form forcing vector 
    p_to_F = [ 0,       0      ; % pressure to forces and moments
               21,      9      ;   
               3*L_e,   2*L_e  ;  
               0,       0      ; 
               9,       21     ;  
              -2*L_e   -3*L_e  ] * L_e / 60 ;

    c_L = p_to_F (1:n_i, 1:n_n) * Line_e'   % force moment @ nodes
  end % if for set up line load nodal resultants
         
%        Optional local mass matrix
  if ( Rho > 0 )
    m_L = [140, 0,       0,       70,  0,       0       ;
           0,   156,     22*L_e,  0,   54,     -13*L_e  ;
           0,   22*L_e,  4*L_e^2, 0,   13*L_e, -3*L_e^2 ;
           70,  0,       0,       140, 0,       0       ;
           0,   54,      13*L_e,  0,   156,    -22*L_e  ;
           0,  -13*L_e, -3*L_e^2, 0,  -22*L_e,  4*L_e^2 ]*Rho*A*L_e
    disp(m_L)
  end % if mass requested

%     Define local to system DOF transformation matrix
  t_L = [ cx   cy   0   0   0   0 ; % (inverse t_L = transpose t_L)
         -cy   cx   0   0   0   0 ;
          0    0    1   0   0   0 ;
          0    0    0   cx  cy  0 ; % joint 2
          0    0    0  -cy  cx  0 ;
          0    0    0   0   0   1 ] ;  

%         Transform from local to system
  S_e = t_L' * s_L * t_L   
disp(S_e/250000)
  M_e = t_L' * m_L * t_L ;
  C_e = t_L' * c_L  

%         SCATTER TO (ASSEMBLE INTO) SYSTEM ARRAYS
% Insert completed element matrices into system matrices
  [rows] = get_element_index (n_g, n_n, e_nodes); % eq numbers
  S (rows, rows) = S (rows, rows) + S_e  ; % add to system sq
  C (rows) = C (rows) + C_e              ; % add to sys column
end % for each j element in mesh                            <<====
disp(S/250000)

%    ALLOCATE STORAGE FOR OPTIONAL REACTION RECOVERY          
if ( EBC_count > 0 )                        ; % reactions occur
  [EBC_row, EBC_col] = save_reaction_matrices (EBC_flag, S, C);
end % if essential BC exist (almost always true)

%         ENFORCE ESSENTIAL BOUNDARY CONDITIONS
save_resultant_load_vectors (n_g, C)
[S, C] = enforce_essential_BC (EBC_flag, EBC_value, S, C);
disp(S/250000)

%             COMPUTE SOLUTION & SAVE
T = S \ C                 ; % Compute displacements & rotations
list_save_displacements_results (n_g, n_m, T)

%          OPTIONAL REACTION RECOVERY & SAVE
if ( EBC_count > 0 )                    ; % reactions exist ?
  [EBC_react] = recover_reactions_print_save (n_g, n_d, ...
                EBC_flag, EBC_row, EBC_col, T); % reaction to EBC
end % if EBC exist

%       POST-PROCESS ELEMENT REACTIONS (MEMBER FORCES)
% output_2D_frame_el_reactions (n_e, n_g, n_n, n_q, nodes, x, y, T)

% End finite element calculations. 
% See /home/mech517/public_html/Matlab_Plots for graphic options
% http://www.owlnet.rice.edu/~mech517/help_plot.html for help
% end of Modular_2D_Frame            

%  +++++++++++++ functions in alphabetical order +++++++++++++++++

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) )                        % 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 (EBC_flag, EBC_value, S, C)

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 %g point sources. \n', n_u)
  fprintf ('(Echo of file msh_load_pt.tmp) \n')
  fprintf ('Node, DOF (1=force, 2=couple),  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 ('%g  %g  %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 (n_g, n_m, C)

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 ('\nApplied Displacement Boundary Conditions: %g \n', n_c)
  fprintf ('(Echo of file load msh_ebc.tmp) \n')               
  %b fprintf ('Node, DOF (1=displacement, 2=slope), Value. \n')
  fprintf ('Node, DOF (1=u, 2=v, 3=r), Value. \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 %g & DOF %g. \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 > 0 )) ; % check input data
  if ( EBC_count ~= n_c )                ; % probable user error
    fprintf ('\nWARNING: mismatch in bc_flag count & msh_ebc.tmp \n')
  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 ('(Echo of file msh_typ_nodes.tmp) \n')
fprintf ('Read %g elements with (ignored) type & %g nodes each. \n', ...
          n_e, n_n) 
el_type = round (msh_typ_nodes(:, 1)); % el type number >= 1 
n_t     = max(el_type)                   ; % number of element types 
%b fprintf ('Maximum number of element types = %g. \n', n_t) 
nodes (1:n_e, 1:n_n) = msh_typ_nodes (1:n_e, 2:1+n_n); 
disp(msh_typ_nodes (:, 1:1+n_n)) % echo data
% 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 %g nodes. \n', n_m)  
fprintf ('(Echo of file msh_bc_xyz.tmp) \n')
fprintf ('bc_flag, x-, y-coordinates \n')
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 
%b disp(msh_bc_xyz (:, 1:1+n_s)) ; % echo data
for j = 1:n_m
%bfprintf (' %2.2i  %g         \n', P(j), x(j)             )
  fprintf (' %3.3i  %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 (n_g, n_m, T) 

function list_save_displacements_results (n_g, n_m, T)
  fprintf('\nCalculated Displacements: \n')
  %b fprintf('X_disp      Y_disp      Z_disp   at %g nodes \n', n_m)
  fprintf('U_disp      V_disp      Rotation 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, '%i %g %g %g       \n',j, 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 (T)

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 ('\nRecovered Reactions at Displacement or Slope BC: %g \n', ...
            sum (sum (EBC_flag > 0)))   ; % header
  fprintf ('Node, DOF (1=force, 2=couple), 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 ('Total force and couple = ') ; disp(Totals) ; % echo total
% end recover_reactions_print_save (EBC_row, EBC_col, T)

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))         ; % 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 (S, C, EBC_flag)

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 ('\nResultants of all input sources: \n')
  fprintf ('Node, DOF (1=force, 2=couple), Value \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 (n_g, n_m, C)

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 =====================================