function [] = DC_wires () % Ten wire DC Network by FEA
% J.E. Akin, Rice University revised 1/29/2020
% unknowns are voltages, reactions are DC currents

% Define element connectivity (like msh_typ_nodes.txt)
%           (1)    (3)           El Nodes R
%  120 v BC  *--k2--*--k6--* (6) 1  1-2   36.
%            |      |      |     2  1-3   1.
%            |     k4     k9     3  2-5   2.
%            |      |      |     4  3-4   10.
%           k1   (4)*--k7--* (7) 5  4-5   20.
%            |      |      |     6  3-6   3.
%            |     k5     k10    7  4-7   5.
%            |      |      |     8  5-8   4.
%   0 v  BC  *--k3--*--k8--* (8) 9  6-7   7.
%           (2)    (5)          10  7-8   6,
%       (j) represents wire connect

% Define wire ends 'connection list' (line msh_typ_nodes.txt)
nodes = [1 2; 1 3; 2 5; 3 4; 4 5; 3 6; 4 7; 5 8; 6 7; 7 8]

% get the number of elements, nodes, and number of equations
n_e = size (nodes, 1)  ; % number of elements in mesh
n_m = max (max(nodes)) ; % number of nodes in mesh
n_g = 1                ; % number of unknowns per node
n_d = n_g * n_m        ; % nunmber of equations in the system
fprintf ('System has %i elements          \n', n_e)
fprintf ('System has %i nodes             \n', n_m)
fprintf ('System has %i unknowns per node \n', n_g)
fprintf ('System has %i equations         \n', n_d)

% Define element properties (like msh_properties.txt)
% fprintf ('Wire resistances \n')
R_j = [36. 1. 2. 10. 20. 3. 5. 4. 7. 6.]  % Resistances 

% Define fixed and free voltage degrees of freedom (dof)
% via equation numbers and given BC values (like msh_ebc.txt)
Fixed = [1  2];  Free = [3 4 5 6 7 8] ; % (for n_g=1)
BC_is = [120. 0.] ; % any boundary condition values at Fixed 
%                      ***
% (Edit above two lines to change BC locations or values)
%                      ***

fprintf ('Given DC voltages at nodes \n'); disp(Fixed)
fprintf ('are \n'); disp(BC_is)

% Always allocate number of voltage, U, forces, F, and
% the square matrix size K in K U = F assembly
U = zeros (n_d, 1); F = zeros (n_d, 1); K = zeros (n_d, n_d);

% Set the known voltage values (copy via vector subscript)
U(Fixed) = BC_is (:)   ; % (m) copied into BC locations

% Set non-zero known external currents (like load_pt.txt)
Applied = sum(F)       ; % later compare to reactions
fprintf ('Sum of external currents = %10.4e \n', Applied)

% Loop over elements to assemble (singular) K matrix
for n = 1:n_e ;  % ====> ====> ====> ====> ====> ====> ====>
   % Gather element properties (and usually its coordinates)
   el_nodes = nodes (n, :) ; % gather element connections
   k        = 1/R_j (n)    ; % gather element properties
            ; % (Usually a numerical integration here.)
   K_e      = k * [1  -1;    % compute el stiffness
                  -1   1]  ; % compute el stiffness
   % Scatter element stiffness to system stiffness
   K (el_nodes, el_nodes) = K (el_nodes, el_nodes) + K_e ;
end ; % for n-th element <==== <==== <==== <==== <==== <====

% These matrix equilibrium equations are not unique
% because they do not yet include Dirichlet BCs.)
fprintf('Assembled (singular) conductance matrix \n'); disp(K)
fprintf('External input sources \n');                  disp(F)

% Transfer known BC voltage effects to Free loads
F (Free) = F (Free) - K (Free, Fixed)*U (Fixed)

% List the non-singular partition of the matrix system
fprintf ('Free conductance partition \n'); disp(K(Free,Free))
fprintf ('Free currents \n');              disp(F(Free))

% Invert non-singular K portion, multiply it times the
% Free forces to find the Free voltage
U (Free) = K (Free, Free) \ F (Free) ; % (m)
fprintf ('System DC voltages: \n') ; disp(U)

% Optionally find system reactions at the Fixed voltage
% needed to support the Dirichlet BCs
F(Fixed) = K(Fixed, Fixed)*U(Fixed) + K(Fixed, Free)*U(Free) ;
fprintf ('System currents at Fixed BC: \n'); disp(F(Fixed))
fprintf ('System currents sum      = %10.4e \n', sum(F(Fixed)));
fprintf ('Sum of external currents = %10.4e \n', Applied)

%% Always find the element gradients, here Delta_U/L
%% Loop over elements to compute the gradient
%for n = 1:n_e ;  % ====> ====> ====> ====> ====> ====> ====>
%   el_nodes = nodes (n, :)       ; % gather element nodes
%   L_e = L_j (n)                 ; % gather element size
%   U_e = U(el_nodes)             ; % gather end voltage
%   g_e = (U_e(2) - U_e(1)) / L_e ; % element gradient
%   fprintf ('Element %i gradient = %10.4e \n', n, g_e)
%end ; % for n-th element <==== <==== <==== <==== <==== <====

% Optionally find the wire end currents (el reactions)
% Loop over elements to compute their reactions (for n_g=1)
 fprintf ('Current entering (+) or leaving (-) end nodes \n')
for n = 1:n_e ;  % ====> ====> ====> ====> ====> ====> ====>
   el_nodes = nodes (n, :) ; % gather element connections
   U_e      = U(el_nodes)  ; % gather the end voltages
   k        = 1/R_j (n)    ; % gather element properties
   K_e      = k * [1  -1; -1  1] ; % compute el stiffness
   R_e      = K_e * U_e    ; % recover end reaction currents
   fprintf ('Currents at element number %i \n', n)
   fprintf ('at nodes %i  %i \n', nodes (n, :))
   disp(R_e)
end ; % for n-th element <==== <==== <==== <==== <==== <====

% Running gives
% System currents at Fixed BC:
%  11.4608
% -11.4608
% end DC_wires