Main

%% If number of element 20
clc
clear
%% Output
[L, phi, Tglob_new1, Kglob_new1,Tglob,Tglob_new,Kglob,Te]=Assignment1();
x = linspace(0, L, numel(phi));
plot(x, phi, '-o');
xlabel('Position along the bar (m)');
ylabel('Angle of twist (phi)');
title('Angle of twist along the bar');
grid on;
fprintf('Values of phi:\n');
for i = 1:numel(phi)
 fprintf('phi(%d) = %.4f\n', i, phi(i));
end

% Input data
function [n,m,n_total,d,d_total,C,L,Le,J0,G0,T0,Ji,Gi,Je,Ge,x,xe]=Input(~)
% Mesh data
n=20; %num of elements
m=2;%num of nodes per element
n_total=n*(m-1)+1;
d=1;%degree of freedom
d_total=n_total*d;
% connectivity matrix
C=zeros(n,2);
for i=1:n
 C(i,1)=i;
 C(i,2)=i+1;
end
% Element size
L=1; %total size
Le=L/(n*(m-1)) ; %length of element
% Global coordinate vector
x=zeros(n_total,1);
for i=0:n_total-1
 x(i+1)=i*L/(n_total-1);
end
% Geometry and material data
J0=7e-7; % torsion constant
G0=80e10;% shear modulus of shaft material
T0=600; % Torque
% calculating value of J,G at global nodes
Ji=zeros(n_total,1);
Gi=zeros(n_total,1);
for i=1:n_total
 Ji(i)=J0*((1+(x(i)/L))^2);
 Gi(i)=G0*(1+(x(i)/L));
end
%calculating value at local nodes
Je=zeros(n,m);
Ge=zeros(n,m);
xe=zeros(n,m);
for e=1:n
for i=1:m
 j=C(e,i);
 Je(e,i)=Ji(j,1);
 Ge(e,i)=Gi(j,1);
 xe(e,i)=x(j,1);
end
end
% Shape Function
function [N,B]=Shape_function(Z,m,ng)
N=zeros(m,ng);
B=zeros(m,ng);
for i=1:ng
 N(1,i)=0.5*(1-Z(i));
 N(2,i)=0.5*(1+Z(i));
 B(1,i)=-0.5;
 B(2,i)=0.5;
end
% Gauss_Legendre
function [Z,w]= Gauss_legendre(ng)
Z=zeros(ng,1);
w=zeros(ng,1);
if ng==1
 Z(1)=0; w(1)=2;
elseif ng==2
 Z(1)=-1/sqrt(3); w(1)=1;
 Z(2)=1/sqrt(3); w(2)=1;
end 
function x = Gauss_elimination(A, b)
 n = numel(b);
 x = zeros(n, 1);
 Aug = [A b];

 % Conventional Gaussian Elimination
 for j = 1:n-1
 for i = j+1:n
 factor = Aug(i, j) / Aug(j, j);
 Aug(i, :) = Aug(i, :) - factor * Aug(j, :);
 end
 end

 % Backward substitution
 x(n) = Aug(n, n+1) / Aug(n, n);
 for k = n-1:-1:1
 x(k) = (Aug(k, n+1) - Aug(k, k+1:n) * x(k+1:n)) / Aug(k, k);
 end
end

function [L, phi, Tglob_new1, Kglob_new1,Tglob,Tglob_new,Kglob,Te] = Assignment1()
 % input data
 [n, m, n_total, ~, ~, ~, L, ~, ~, ~, T0, ~, ~, Je, Ge, ~, ~] = Input();

 % calculation of gauss points
 ng = (4 * m - 5) / 2 + 0.5;
 [Z, w] = Gauss_legendre(ng);
 [N, B] = Shape_function(Z, m, ng);
 w(ng,1);

 % Elemental matrix
 Ke = zeros(m, m); % elemental stiffness matrix
 Kglob = zeros(n_total, n_total); % global stiffness matrix
 Te = zeros(m, m); % elemental force matrix
 Tglob = zeros(n_total, 1); % global force matrix
for e=1:n
 Te=0;
 for i=1:ng
 G=Ge(e,:)*N(:,i);
 J=Je(e,:)*N(:,i);
 Ke= Ke + (w(i,1)*G*J*(B(:,i))*B(:,i)');
 Te= Te + (w(i,1)*T0.*(N(:,i)));
 end

 Kglob(e:e+1,e:e+1) =Kglob(e:e+1,e:e+1)+ Ke;
 Tglob(e:e+1,1) = Tglob(e:e+1,1)+Te;
 end
 phi_1=0.1;% initial rotation
 Tglob(n_total,1)=Tglob(n_total,1)+(2*T0/3);
 Kglob_new=Kglob(1:end,2:end);
 Tglob_new=Tglob(1:end,1)-(phi_1*Kglob(1:end,1));
 Kglob_new1=Kglob_new'*Kglob_new;
 Tglob_new1=Kglob_new'*Tglob_new;
 phi=Gauss_elimination(Kglob_new1,Tglob_new1);
 phi=[phi_1;phi];

end