function Poisson_T6_Isoparametric  (load_pt) % Revised 3/31/09 ===
% load_pt = 1 flags that point sources will be  
%  input via text file msh_load_pt.tmp 
%............................................................. 
%     Heat Conduction Analysis  (2.5D), Isoparametric T6,
%     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 source 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

 %  set application (Poisson Eq) controls
n_g = 1     ; % number of DOF per node 
n_q = 4     ; % number of quadrature points required 
n_r = 2     ; % number of rows in B_e matrix 
post = 1    ; % turn on(1)/off(0) post-processing option
patch = 1   ; % turn on(1)/off(0) patch test option
echo_p = 0  ; % turn on(1)/off(0) listing of properties

%    *** 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  

%     Count any MPC flags within packed integer BC flag P
MPC_max = max( max (EBC_flag) )  ; % max type MPC (1 means EBC only)
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_each  = count of constraints with 1 to MPC_max DOF terms
% MPC_max   = max number of DOF in a single MPC equation

%   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
Type_volume = zeros (n_t, 1) ; Type_integral = zeros (n_t, 1) ;
% Type_volume = volume of each element type group   
% Type_integral = integral of solution in each element type group

%   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 

%   Manual patch test of the elements (iff patch == 1)
if ( patch == 1 ) % set EBC_values for default patch test
  fprintf ('\nNOTE converted EBC values to standard patch test \n')
  for j = 1:n_m % loop over all nodes to find boundary
    if any ( EBC_flag (j, :) > 0 ) % then on boundary             
       EBC_value (j, 1) = 1 + 2 * x(j) + 3 * y(j) ; % required
    end % if boundary node
  end % for j
  disp (EBC_value)
  fprintf ('Constant gradients are G_x=2, G_y=3 \n')
end % if patch test
 
%    Read point source 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_mats = size(msh_properties, 1)              ; % == 1 if homogeneous 
n_vals = size(msh_properties, 2)        ; % of application properties
% n_mats = number of materials, = n_e or n_t if homogeneous
if ( n_mats == n_t )        % then element type data are honogeneous
  fprintf ('%i homogeneous material values per type \n', n_vals)
else
  fprintf ('%i material values per element \n', n_vals)
end % if
format short e ;  disp(msh_properties) ; format short
if ( n_mats > n_e ) 
  error ('ERROR: number of material sets exceeds number of elements') 
end % if 

if ( echo_p == 0 )
  fprintf ('\nNOTE individual element properties echo skipped \n')
end % if

%   for Poisson Eq. application only
if ( echo_p )
 if ( n_mats == n_t )   % then element type data are honogeneous
   for e_t = 1:n_t % loop over types
     K      = msh_properties (e_t, 1)   ; % conductivity, isotropic
     Q      = msh_properties (e_t, 2)   ; % heat generation per unit volume
     Thick  = msh_properties (e_t, 3)   ; % thickness

%      populate the constitutive law matrix
     E_e = K * eye (n_r, n_r) ;

%                   ECHO PROPERTIES 
     fprintf (' \n') 
     fprintf ('Homogeneous Material Values: \n' ) 
     fprintf ('\nProperties for all elements \n')
     fprintf ('Thermal conductivity        = %g \n', K)     
     fprintf ('Heat generation rate        = %g \n', Q)    
     fprintf ('Thickness                   = %g \n', Thick)
     if ( Thick <= 0 ) 
       fprintf ('\nError: zero thickness at element %g \n', j)
       fprintf ('reset to default unit thickness \n')
       Thick = 1 ;
     end % if invalid data
   end % for each element type
 end % if homogeneous element type properties
end % if echo
%   end for plane stress 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 ***

%     open unit needed for post-processing stresses
if ( post == 1 ) % open post-processing unit
  qp_unit = fopen ('qp_store.bin', 'w') ; % open for writing
  if ( qp_unit == -1 ) % open failed
    fprintf ('\nWARNING: post-processing unit open failed \n')
    post = 0   % turn off post-processing option
  end % if unit available
end % if open post-processing unit

%    *** Initially fill the quadrature data files
r_q = zeros (n_q, 1) ;  % allocate quadrature storage
s_q = zeros (n_q, 1) ;  % allocate quadrature storage
t_q = zeros (n_q, 1) ;  % allocate quadrature storage
w_q = zeros (n_q, 1) ;  % allocate quadrature storage
% r_q, s_q, t_q = tabulated 1-, 2-, and 3-D quadrature locations
% w_q           = tabulated quadrature weights  

if ( n_q > 0 )                 % numerical integration to be used
  if ( n_s == 2 )                   % triangles or quadrilaterals
    if ( n_n == 3 | n_n == 6 )                   % then triangles
      [r_q, s_q, w_q] = qp_rule_unit_tri (n_q) ;
    elseif ( n_n == 4 | n_n == 8 | n_n == 9 )   % then quadrilaterals
      %b [r_q, s_q, w_q] = qp_rule_natural_quad (n_q) ;
      error ('\nERROR fill quadrilateral quadrature rules')
    end % if quadrilaterals
  else
    error ('\nERROR fill 1D or 3D quadrature rules')
  end % if space dimension
end % if numerical integration data required 
 
%  ** GENERATE ELEMENT MATRICES 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
   
%                 Gather nodal coordinates
  xy_e (1:n_n, 1) = x(e_nodes)             ; % x coord at el nodes 
  xy_e (1:n_n, 2) = y(e_nodes)   ; % 6 x 2   % y coord at el nodes

%     Allocate and clear element arrays
  H_q   = zeros (n_n, 1)      ; % clear element interpolations
  H_i   = zeros (n_n, 1)      ; % element interpolations integral
  DLH_q = zeros (n_s, n_n)    ; % interpolation local derivatives 
  Jacobian = zeros (n_s, n_s) ; % geometry Jacobian
  DGH_q = zeros (n_s, n_n)    ; % interpolation global derivatives 
  B_q   = zeros (n_r, n_i)    ; % initalize B matrix values
  S_e   = zeros (n_i, n_i)    ; % clear array el stiffness
  M_e   = zeros (n_i, n_i)    ; % clear array el mass
  C_e   = zeros (n_i, 1)      ; % clear array el sources

%b size(H_q)
%b size(H_i)
%               SET ELEMENT MATERAILS & GEOMETRY 
  if ( n_mats > n_t )    % each element has some different properties
    J = j ; % the element number
  else % use type data
    J = type ; % the type number
  end % if which property values

  K      = msh_properties (J, 1)   ; % modulus of elasticity
  Q      = msh_properties (J, 2)   ; % Poisson's ratio       
  Thick  = msh_properties (J, 3)   ; % thickness

%      populate the constitutive law matrix
     E_e = K * eye (n_r, n_r) ;

  if ( echo_p == 1 ) 
%                   ECHO PROPERTIES
    fprintf ('Element number %i \n', j)           
    fprintf ('Thermal conductivity        = %g \n', K)     
    fprintf ('Heat generation rate        = %g \n', Q)    
    fprintf ('Thickness                   = %g \n', Thick)
  end % if echo
  if ( Thick <= 0 ) 
    fprintf ('\nError: zero thickness at element %g \n', j)
    fprintf ('reset to default unit thickness \n')
    Thick = 1 ;
  end % if invalid data
 
%   ELEMENT, FOUNDATION STIFFNESSES AND SOURCE RESULTANT MATRICES 

%            Numerical Intergation of S_e, C_e 
  volume = 0 ;
  for q = 1:n_q   ; % begin quadrature loop ---> ---> ---> ---> ---> --->
    r = r_q (q) ; s = s_q (q) ; w = w_q (q)     ; % recover data

%                 element type interpolation functions, 1 x 6
    H_q = [ (1 - 3*r - 3*s + 2*r*r + 4*r*s + 2*s*s), (2*r*r - r), ...
            (2*s*s - s), 4*(r - r*r - r*s), 4*r*s, 4*(s - r*s - s*s) ] ;

%                 element type parametric derivatives, 2 x 6
    DLH_q = [ (-3 + 4 * r + 4 * s), (-1 + 4 * r), 0, ...
              (4 - 8 * r - 4 * s), (4 * s), (-4 * s) ; ... % dH/dr
              (-3 + 4 * r + 4 * s), 0, (-1 + 4 * s), ...
              (-4 * r), (4 * r), (4 - 4 * r - 8 * s) ] ;   % dH/ds

%             Global position and Jacobian
    xy_q     = H_q * xy_e           ; % 1 x 2     % global (x, y)
    Jacobian = DLH_q * xy_e         ; % 2 x 2     % d(x,y)/d(r,s)
    J_det    = abs( det (Jacobian)) ; % 1 x 1     % determinant
    J_Inv    = inv (Jacobian)       ; % 2 x 2     % inverse

%             Global derivatives of H
    DGH_q = J_Inv * DLH_q      ; % 2 x 6          % dH/dx & dH/dy
    E_q = E_e                  ; % 3 x 3          % properties
    T_q = Thick                ; % 1 x 1          % thickness
    
%      Insert gradient terms into the differential operator
    B_q = DGH_q ; % always the same in Poisson equation

%                 ELEMENT MATRICES UPDATES, 12 x 12 and 12 x 1
    S_e = S_e + (B_q' * E_q * B_q) * T_q * J_det * w_q (q) ; % square
 %b M_e = M_e + (H_q' * Rho * H_q) * T_q * J_det * w_q (q) ; % mass  
    if ( Q ~= 0 ) % then need loads
      C_e = C_e + H_q' * Q * T_q * J_det * w_q (q); % F_Q
    end % if heat generationactive
%b size(H_q)
%b size(H_i)
    %b H_i(1:n_n) = H_i(1:n_n) + H_q * T_q * J_det * w_q (q); % H integral
    H_i = H_i + H_q' * T_q * J_det * w_q (q) ; % H integral
    volume = volume + T_q * J_det * w_q (q)  ;

%     save data needed for post-processing stresses
    if ( post == 1 ) % save data for this point
      fwrite (qp_unit, xy_q,  'double') ; % save coordinates
      fwrite (qp_unit, E_q,   'double') ; % save material   
      fwrite (qp_unit, B_q,   'double') ; % save gradient info
    end % if save point data
   end % for quadratures   <--- <--- <--- <--- <--- <--- <--- <--- <--- 
   fwrite (qp_unit, H_i,    'double') ; % save interpolation integral
   fwrite (qp_unit, volume, 'double') ; % save element volume

%         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 
  C (rows)       = C (rows)       + C_e  ; % add to sys force/couples 
%b M (rows, rows) = M (rows, rows) + M_e  ; % add to system mass 
%       ***  END SCATTER TO SYSTEM ARRAYS ***
 
end % for each j element in mesh   <<==== <<==== <<==== <<==== <<==== 
%       *** END ELEMENT GENERATION
 
%    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) 
 
%             COMPUTE SOLUTION & SAVE 
T = S \ C                 ; % Compute displacements & rotations 
list_save_results (n_g, n_m, T)              ; % save and print    
 
%    OPTIONAL FORCE 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 ELEMENTS
if ( post == 1 ) % then post-process all element stresses
  if ( qp_unit > 0 ) ; % close written files
    fclose (qp_unit) ; % file qp_store.bin
    qp_unit = fopen ('qp_store.bin', 'r') ; % open for reading 
  end % if
  qid = fopen ('el_qp_xyz_fluxes.tmp', 'w') ; % plot file

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

%  Re-compute or read arrays from storage
  format short e
  for j = 1:n_e ; % loop over elements ====>> ====>> ====>> ====>> ====>> 
    fprintf ('Element number %i _______________ \n', j)
    fprintf ('Point       x-, y-coordinates \n')
    fprintf ('            Gradient vector \n')
    fprintf ('            Heat Flux vector \n')
    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 
   
    for q = 1:n_q   ; % begin quadrature loop ---> ---> ---> ---> ---> --->
      [xy_q] = fread (qp_unit, [1, n_s],   'double') ; % get coordinates
      [E_q]  = fread (qp_unit, [n_r, n_r], 'double') ; % get material
      [B_q]  = fread (qp_unit, [n_r, n_i], 'double') ; % get grad info  

%     compute the gradient, heat flux, solution integral 
      Grad = B_q * T_e'  ;  % in-plane Gradient
      Flux = -E_q * Grad ;  % in-plane Fluxes (Fourier law)

      fprintf ('%i, %g, %g \n', q, xy_q (1:n_s))  % list
      fprintf ('%i, %g, %g \n', q, Grad (1:n_r))  % list
      fprintf ('%i, %g, %g \n', q, Flux (1:n_r))  % list
      fprintf (qid, '%g, %g, %g, %g, %g, %g \n', ...  % save for plots
               xy_q (1:2), Grad (1:2), Flux (1:2)) ;

    end % for q quadrature points <--- <--- <--- <--- <--- <--- <--- <---

    [H_i]    = fread (qp_unit, [n_n, 1], 'double') ; % interp integral
    [volume] = fread (qp_unit, [1, 1],   'double') ; % element volume
%b size(H_i)
%b size(T_e)
    Type_volume   (type) = Type_volume   (type) + volume    ;
    Type_integral (type) = Type_integral (type) + T_e * H_i ;

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

fprintf ('\n')
if ( n_t > 1 ) % report on each type
  for j = 1:n_t
    fprintf ('Solution integral in type %i elements = %g \n', ...
              j, Type_integral (j))
    fprintf ('Solution volume   in type %i elements = %g \n', ...
              j, Type_volume (j))
  end % for each type
end % if multiple element types
fprintf ('Total solution integral = %g \n', sum (Type_integral))
fprintf ('Total solution volume   = %g \n', sum (Type_volume))

  format short
%       OPTIONAL POST_PROCESSING COMPLETED
else
  fprintf ('\nNo element post-processing requested \n')
end % for post-processing stresses

% end function Poisson_T6_Isopar  % =======================
   
%  +++++++++++++ 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 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 use 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_each  = count constraints with 1 to MPC_max DOF terms (MPC_max, 1)
% 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 
 if ( n_c == 0 )                     % default to all zero values
   fprintf ('\nNOTE defaulting all EBC values to zero \n')
   return
 else
   fprintf ('\nRead %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 % if default to zero
 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 type & %i nodes each. \n', ...
          n_e, n_n)  
fprintf ('(Echo element number, file msh_typ_nodes.tmp) \n')
el_type = round (msh_typ_nodes(:, 1)); % el type number >= 1  
nodes (1:n_e, 1:n_n) = msh_typ_nodes (1:n_e, 2:1+n_n);  

if ( n_n > 9 ) % skip element number display
 disp(msh_typ_nodes (:, 1:1+n_n)) % echo data 
else
 for i = 1:n_e
  if ( n_n == 2 )
   fprintf ('%i, %i %i %i \n', i, msh_typ_nodes (i, :))
  elseif ( n_n == 3 )
   fprintf ('%i, %i %i %i %i \n', i, msh_typ_nodes (i, :))
  elseif ( n_n == 4 )
   fprintf ('%i, %i %i %i %i %i \n', i, msh_typ_nodes (i, :))
  elseif ( n_n == 6 )
   fprintf ('%i, %i %i %i %i %i %i %i \n', ...
              i, msh_typ_nodes (i, :))
  elseif ( n_n == 8 )
   fprintf ('%i, %i %i %i %i %i %i %i %i %i \n', ...
              i, msh_typ_nodes (i, :))
  elseif ( n_n == 9 )
   fprintf ('%i, %i %i %i %i %i %i %i %i %i %i \n', ...
              i, msh_typ_nodes (i, :))
  end % if number of nodes to show
 end % for all elements
end % if show element number

n_t  = max(el_type)                   ; % number of element types  
if ( n_t > 1 )
  fprintf ('Note: maximum element type is %i \n', n_t)
end % if
% 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, file msh_bc_xyz.tmp) \n')
%B fprintf ('bc_flags, x-coordinate \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
%bfprintf ('%i, %2.2i  %g         \n', j, P(j), x(j)             )
%bfprintf ('%i, %2.2i  %g  %g     \n', j, P(j), x(j), y(j)       )
  fprintf ('%i, %i  %g  %g     \n', j, P(j), x(j), y(j)       )
%bfprintf ('%i, %2.2i  %g  %g  %g \n', j, P(j), x(j), y(j), z(j) )
end % for
% end get_mesh_nodes    % ==========================
 
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 list_save_results (n_g, n_m, T) % ============
  fprintf ('\nComputed Solution: \n') ;
  fprintf('Node, %i results per node \n', n_g)
  T_matrix = reshape (T, n_g, n_m)'            ; % pretty shape
  %b disp (T_matrix)                              ; % print results
  % 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)) ;
      fprintf ('%i   %g                \n', j, T_matrix (j, 1:n_g)) ;
    elseif ( n_g == 2 )
      fprintf (fid, '%g %g             \n', T_matrix (j, 1:n_g)) ;
      fprintf ('%i   %g %g             \n', j, T_matrix (j, 1:n_g)) ;
    elseif ( n_g == 3 )
      fprintf (fid, '%g %g %g          \n', T_matrix (j, 1:n_g)) ;
      fprintf ('%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)) ;
      fprintf ('%i   %g %g %g %g       \n', j, T_matrix (j, 1:n_g)) ;
    elseif ( n_g == 5 )
      fprintf (fid, '%g %g %g %g %g    \n', T_matrix (j, 1:n_g)) ;
      fprintf ('%i   %g %g %g %g %g    \n', j, 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)) ;
      fprintf ('%i %g %g %g %g %g %g \n', j, T_matrix (j, 1:n_g)) ;
    else 
      error ('reformat list_save_results for n_g > 6.')
    end % if 
  end % for j DOF
% end list_save_results    % ==========================

function [r_q, s_q, w_q] = qp_rule_unit_tri (n_q) % ================
%  tables of quadrature point locations and weights
%  for triangles interpolated in unit coordinates
if ( n_q == 1 )            % exact for constant or linear polynomial
  r_q = [ 1 / 3.] ; s_r = [ 1 / 3.] ; w_q = [ 1 / 2.] ;
elseif ( n_q == 4 )              % exact for 3-rd degree polynomials
  r_q = [ 1.0, 0.6, 0.6, 1.8 ] / 3 ;
  s_q = [ 1.0, 0.6, 1.8, 0.6 ] / 3 ;
  w_q = [ -27, 25,  25,  25  ] / 96 ;
elseif ( n_q == 7 )              % exact for 4-th degree polynomials
  r_q = [ 0.10128651, 0.47014206, 0.79742699, 0.47014206, ...
          0.10128651, 0.05971587, 0.33333333 ] ;
  s_q = [ 0.10128651, 0.05971587, 0.10128651, 0.47014206, ...
          0.79742699, 0.47014206, 0.33333333 ] ;
  w_q = [ 0.06296959, 0.06619708, 0.06296959, 0.06619708, ...
          0.06296959, 0.06619708, 0.11250000 ] ;
else % update the tables
  error ('\nERROR quadrature rule not tabulated for these points')
end % if number of quadrature points
% end function    % ===============================================

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