clear
clc

% **********BY S Ziaei Rad **************************
% Example 3.4.1 - Direct Iteration Method
% 2 noded elements- Page 66: An Introduction to Nonlinear FE by Reddy
% ***************************************************
% -d/dx(udu/dx)=-1 0<x<1
% u(1)=sqrt(2);
% Q(0)=0

% nnodes  =  No. of Nodes
% nelems  = No. of Elements
nnodes=21
nelems=20
L=1;

% Essential BCs
BCs=zeros(1,nnodes)
BCs(1,nnodes)=sqrt(2)
%BCs=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 sqrt(2)];
u(1:nnodes,1)=ones*1;
u(nnodes,1)=sqrt(2);

h(1:nelems)=ones(nelems,1)*L/nelems;
f(1:nelems)=-h/2;
% Example of connectivity for 20 elements
% connect=[
%     1 1 2 h(1) 
%     2 2 3 h(2)
%     3 3 4 h(3)
%     4 4 5 h(4)
%     5 5 6 h(5)
%     6 6 7 h(6)
%     7 7 8 h(7)
%     8 8 9 h(8)
%     9 9 10 h(9)
%     10 10 11 h(10)
%     11 11 12 h(11)
%     12 12 13 h(12)
%     13 13 14 h(13)
%     14 14 15 h(14)
%     15 15 16 h(15)
%     16 16 17 h(16)
%     17 17 18 h(17)
%     18 18 19 h(18)
%     19 19 20 h(19)
%     20 20 21 h(20)
%     ];

% Connectivity matrix
% [ No. of Element  Node#1 Node#2   Element Length]
connect=[];
for i=1:nelems
    connect=[connect
        i i i+1 h(i)];
end

iter1=0;
UU=[];
epsilon11=1e-6
snorm1=1;

% Direct Iteration Method
while (snorm1)>=epsilon11

iter1=iter1+1
pause
K=zeros(nnodes,nnodes);
F=zeros(nnodes,1);
UU=[UU u];

% Constructing Element Stiffness Matrix, Element Force Vector ...
%  ... Global Stiffness Matrix and Global Force Vector
  for i=1:nelems
    [K,F]=assemble(K,F,connect,i,u);
  end
ZBCs=find(BCs==0);
NBCs=find(BCs~=0);
UBC=BCs(NBCs);
Kaa=K(ZBCs,ZBCs);
Faa=F(ZBCs);
Kab=K(ZBCs,NBCs);
u1=inv(Kaa)*(Faa-Kab*UBC);
u=[u1 ; UBC]
snorm1=norm(u-UU(:,iter1))/norm(u)

end

% Ploting The results
% The exact solution for this problem is u(x)=sqrt(1+x^2)

UZ=[0:1/nelems:1]';
SS=[UZ UU]
plot(SS(:,1),SS(:,iter1),SS(:,1),sqrt(1+SS(:,1).^2),'*r')

xlabel('x')
ylabel('u(x)')
legend('FE solution', 'Exact Solution')

    