% File: Uxx_U_Q_EBC_NBC.m for ODE E*Uxx + k*U + Q(x) = 0, E=1=k, 
% plus U(0)=0, dU/dx(L)=0 so U(x) = sin(x)/cos(L) - x
% Ref: J.E. Akin, "FEA w Error Estimators" Elsevier, 2005, page 50  
%.......................................................................
%      x_1=0     x_2       x_3       x_4 = L
%  U_o  *---(1)---*---(2)---*---(3)---*   dU(L)/dx
%       T_1       T_2       T_3       T_4
%           (L2)      (L2)      (L2)      (P1)
% Connectivity list:    e     i      j
%                       1     1      2   % L2
%                       2     2      3   % L2
%                       3     3      4   % L2
%                       4     4      0   % P1
%.......................................................................

%            Set given constants (n_e can be changed)
n_e = 3 ; n_d = n_e + 1 ;  % number of elements and nodes (DOF)
n_n = 2                 ;  % number of nodes per element
L = 1.0 ; L_e = L/n_e   ;  % domain and element lengths
E = 1.0 ; k = 1.0       ;  % stiffness data
U_o  = 0.0              ;  % essential boundary condition
Ux_L = 0.0              ;  % flux BC
X = [0:n_e] * L_e       ;  % merge lengths

%         Constant element square matrices
K_E = E   * reshape ([1, -1, -1, 1], 2, 2) / L_e; % axial
M_e = L_e * reshape ([2,  1,  1, 2], 2, 2) / 6  ; % mass
K_k = k * M_e                                   ; % foundation
S_e = K_E - K_k                                 ; % net square matrix
Q_e = zeros (n_n,1)                             ; % pt sources
C_e = zeros (n_n,1)                             ; % elem sources

%      Assemble n_d by n_d square matrix terms
S = zeros (n_d, n_d) ;  C = zeros (n_d, 1) ; % initalize sums
for j = 1:n_e                          ; % loop over elements
  last = j + n_n - 1 ; first = last - n_n + 1 ; % DOF numbers
  rows = [first last]                  ; % vector subscripts
  S (rows, rows) = S (rows, rows) + S_e; % add to system sq
  Q_e = [X(first);  X(last)]           ; % nodal source
  C_e = M_e * Q_e                      ; % element source
  C (rows) = C (rows) + C_e            ; % add to sys column
end % for loop over elements
C (n_d) = C (n_d) + Ux_L               ; % add point source

%  Save reaction data to be destroyed by EBC solver trick
Row_1_S = S(1, 1:n_d) ; Row_1_C = C (1) ;

%  Apply essential BC at node 1, via trick to avoid partitions
C (:) = C (:) - U_o * S (:, 1)   ; % carry known column to RHS
S (:, 1) = 0 ; S (1, :) = 0      ; % clear, restore symmetry
S (1, 1) = 1 ; C (1) = U_o       ; % trick not to change array size
T = S \ C                        ; % Solve for four values
q_in = -( Row_1_S * T - Row_1_C ); % Get reaction from matrix system
 
%            Plot the temperature results
clf ; hold on ; grid on                    % clear & keep
title ('ODE E*Uxx + k*U + Q(x), E=k=1, Q(x) = x')
xlabel ('Position') ; ylabel ('U(x)')
x = [0:0.05:1] * L ; exact = sin(x)/cos(L) - x ; % analytic
plot (x, exact, 'k.-');  plot (X, T, 'ko')       % exact & nodal
legend ('Exact', 'FEA'); plot (X, T, 'k-')       % legend & elements
e_text = ['U(0)=0=dU(L)/dx, Reaction = ', num2str(q_in)] ;
text (mean(X), mean(T), e_text)                  % list BCs
% end File: Uxx_U_Q_EBC_NBC.m 


%    >> Uxx_U_Q_EBC_NBC   (with n_e = 3)
%    K_E =  3    -3
%          -3     3
%    M_e = 0.1111    0.0556
%          0.0556    0.1111
%    S_e = 2.8889   -3.0556        
%         -3.0556    2.8889         
%    % assembly loop
%    S = 2.8889   -3.0556         0         0
%       -3.0556    2.8889         0         0
%             0         0         0         0
%             0         0         0         0
%    Q_e = 0 
%          0.3333
%    C_e = 0.0185 
%          0.0370
%    C = 0.0185
%        0.0370
%             0
%             0
%    % assembly loop
%    S = 2.8889   -3.0556         0         0
%       -3.0556    5.7778   -3.0556         0
%             0   -3.0556    2.8889         0
%             0         0         0         0
%    Q_e = 0.3333
%          0.6667
%    C_e = 0.0741
%          0.0926
%    C = 0.0185
%        0.1111
%        0.0926
%             0
%    % assembly loop
%    S = 2.8889   -3.0556         0         0
%       -3.0556    5.7778   -3.0556         0
%             0   -3.0556    5.7778   -3.0556
%             0         0   -3.0556    2.8889
%    Q_e = 0.6667
%          1.0000
%    C_e = 0.1296
%          0.1481
%    C = 0.0185
%        0.1111
%        0.2222
%        0.1481
%    % end assembly
%    T = 0
%        0.2681
%        0.4706
%        0.5490
%    q_in = 0.8377