function TEST_CODE2()
R = 10 / 1000;
H = 20 / 1000;
%number of elements in z direction
N_z=5;
%number of elements in r direction
N_r=10;
%number of total elements
N_E=(N_z)*(N_r);
%number of nodes in z direction
NN_z=(N_z)+1;
%number of nodes in r direction
NN_r=(N_r)+1;
% number of total nodes
NN=NN_z*NN_r;
%r = zeros(NN, 1);
theta = zeros(NN, 1);
z = zeros(NN, 1);
IEN = zeros(N_E*4, 1);
%a=(2*pi*R)/(N_r);
b=H/(N_z);
c=(2*pi)/(N_r);
N_z_index = 1;
z_value = 0;
%r_value = 0;
theta_value = 0;
for i = 1:1:NN
% r(i) = r_value;
z(i) = z_value;
theta(i) = theta_value;
if N_z_index >= NN_z
N_z_index = 1;
z_value = 0;
%r_value = r_value + a;
theta_value = theta_value + c;
else
z_value = z_value + b;
N_z_index = N_z_index + 1;
end
end
IEN_count = 1;
IEN_index = 1;
for i = 1:1:N_E
temp = 4*(i-1);
IEN(4*(i-1) + 1) = IEN_index;
IEN(4*(i-1) + 2) = IEN_index + NN_z;
IEN(4*(i-1) + 3) = IEN_index + NN_z + 1;
IEN(4*(i-1) + 4) = IEN_index + 1;
if IEN_count == N_z
IEN_index = IEN_index + 2;
IEN_count = 1;
else
IEN_index = IEN_index + 1;
IEN_count = IEN_count + 1;
end
end
h = zeros(N_E, 1);
h_theta = zeros(N_E, 1);
mu=1;
c=2;
e=3;
for i = 1:1:N_E
F1 = 0;
k = 0;
F2 = 0;
F = 0;
for xi=(-1/sqrt(3)):(2/sqrt(3)):(1/sqrt(3))
for eta=(-1/sqrt(3)):(2/sqrt(3)):(1/sqrt(3))
NodeNum1 = IEN(4*(i-1) + 1);
NodeNum2 = IEN(4*(i-1) + 2);
NodeNum3 = IEN(4*(i-1) + 3);
NodeNum4 = IEN(4*(i-1) + 4);
theta1 = theta(NodeNum1);
theta2 = theta(NodeNum2);
theta3 = theta(NodeNum3);
theta4 = theta(NodeNum4);
z1 = z(NodeNum1);
z2 = z(NodeNum2);
z3 = z(NodeNum3);
z4 = z(NodeNum4);
h1=c+e*cos(theta1);
h2=c+e*cos(theta2);
h3=c+e*cos(theta3);
h4=c+e*cos(theta4);
h = (h1+h2+h3+h4)/4;
h1_theta = e*cos(theta1)-e*sin(theta1);
h2_theta = e*cos(theta2)-e*sin(theta2);
h3_theta = e*cos(theta3)-e*sin(theta3);
h4_theta = e*cos(theta4)-e*sin(theta4);
h_theta = (h1_theta+h2_theta+h3_theta+h4_theta)/4
N1=0.25*(1-xi)*(1-eta);
N2=0.25*(1+xi)*(1-eta);
N3=0.25*(1+xi)*(1+eta);
N4=0.25*(1-xi)*(1+eta);
N=[N1,0,N2,0,N3,0,N4,0;
0,N1,0,N2,0,N3,0,N4];
dN1_xi=0.25*(-1+eta);
dN1_eta=0.25*(-1+xi);
dN2_xi=0.25*(1-eta);
dN2_eta=0.25*(-1-xi);
dN3_xi=0.25*(1+eta);
dN3_eta=0.25*(1+xi);
dN4_xi=0.25*(-1-eta);
dN4_eta=0.25*(-1-xi);
dN_xi = zeros(2,8);
dN_xi(1,1) = dN1_xi;
dN_xi(1,2) = 0.0;
dN_xi(1,3) = dN2_xi;
dN_xi(1,1) = 0.0;
dN_xi(1,1) = dN3_xi;
dN_xi(1,1) = 0.0;
dN_xi(1,1) = dN4_xi;
dN_xi(1,8) = 0.0;
dN_xi(2,1) = 0.0;
dN_xi(2,2) = dN1_xi;
dN_xi(2,3) = 0.0;
dN_xi(2,1) = dN2_xi;
dN_xi(2,1) = 0.0;
dN_xi(2,1) = dN3_xi;
dN_xi(2,1) = 0.0;
dN_xi(2,8) = dN4_xi;
dN_eta = zeros(2,8);
dN_eta(1,1) = dN1_eta;
dN_eta(1,2) = 0.0;
dN_eta(1,3) = dN2_eta;
dN_eta(1,1) = 0.0;
dN_eta(1,1) = dN3_eta;
dN_eta(1,1) = 0.0;
dN_eta(1,1) = dN4_eta;
dN_eta(1,8) = 0.0;
dN_eta(2,1) = 0.0;
dN_eta(2,2) = dN1_eta;
dN_eta(2,3) = 0.0;
dN_eta(2,1) = dN2_eta;
dN_eta(2,1) = 0.0;
dN_eta(2,1) = dN3_eta;
dN_eta(2,1) = 0.0;
dN_eta(2,8) = dN4_eta;
J11=-1*(1-eta)*theta1+(1-eta)*theta2+(1+eta)*theta3-(1+eta)*theta4;
J12=-1*(1-eta)*z1+(1-eta)*z2+(1+eta)*z3-(1+eta)*z4;
J21=-1*(1-xi)*theta1-(1+xi)*theta2+(1+xi)*theta3+(1-xi)*theta4;
J22=-1*(1-xi)*z1-(1+xi)*z2+(1+xi)*z3+(1-xi)*z4;
J = zeros(2,2);
J = 0.25*[J11, J12;
J21,J22];
detJ=J11*J22-J21*J12;
%invJ=[J22,-1*J12;
% -1*J21,J11]/detJ;
dN_theta=(J22*dN_xi-J12*dN_eta)/detJ;
dN_z=(-1*J12*dN_xi+J11*dN_eta)/detJ;
F1_xi_eta=(1/R)*(dN_theta)*(dN_theta)';
F2_xi_eta=(dN_z)*(dN_z)';
%h=c+e*cos(theta);
% Gauss integral method
% initialize
k = k + detJ*(h(i)^3)*(F1_xi_eta)*(F2_xi_eta)/(12*mu);
V_theta1 = R*(theta1)/2;
V_theta2 = R*(theta2)/2;
V_theta3 = R*(theta3)/2;
V_theta4 = R*(theta4)/2;
V_theta = (V_theta1+V_theta2+V_theta3+V_theta4)/4
V_z = (e)/2;
F1 = F1+detJ*(h(i))*((F1_xi_eta)*(V_theta)+(F2_xi_eta)*(V_z));
F2 = F2+detJ*(h_theta(i))*N;
F = F1-F2;
end
P = F/k;
end
end
y = StiffnessAssemble(K,k,NodeNum1,NodeNum2,NodeNum3,NodeNum4);
end
% change local to global
%nodes i,j,m,n
function y=StiffnessAssemble(K,k,i,j,m,n)
K = zeros(NN*2,NN*2);
K(2*i-1,2*i-1)= K(2*i-1,2*i-1) + k(1,1);
K(2*i-1,2*1)=K(2*i-1,2*i)+k(1,2);
K(2*i-1,2*j-1)=K(2*i-1,2*j-1)+k(1,3);
K(2*i-1,2*j)=K(2*i-1,2*j)+k(1,4);
K(2*i-1,2*m-1)=K(2*i-1,2*m-1)+k(1,5);
K(2*i-1,2*m)=K(2*i-1,2*m)+k(1,6);
K(2*i-1,2*n-1)=K(2*i-1,2*n-1)+k(1,7);
K(2*i-1,2*n)=K(2*i-1,2*n)+k(1,8);
K(2*i,2*i-1)=K(2*i,2*i-1)+k(2,1);
K(2*i,2*i)=K(2*i,2*i)+k(2,2);
K(2*i,2*j-1)=K(2*i,2*j-1)+k(2,3);
K(2*i,2*j)=K(2*i,2*j)+k(2,4);
K(2*i,2*m-1)=K(2*i,2*m-1)+k(2,5);
K(2*i,2*m)=K(2*i,2*m)+k(2,6);
K(2*i,2*n-1)=K(2*i,2*n-1)+k(2,7);
K(2*i,2*n)=K(2*i,2*n)+k(2,8);
K(2*j-1,2*i-1)=K(2*j-1,2*i-1)+k(3,1);
K(2*j-1,2*i)=K(2*j-1,2*i)+k(3,2);
K(2*j-1,2*j-1)=K(2*j-1,2*j-1)+k(3,3);
K(2*j-1,2*j)=K(2*j-1,2*j)+k(3,4);
K(2*j-1,2*m-1)=K(2*j-1,2*m-1)+k(3,5);
K(2*j-1,2*m)=K(2*j-1,2*m)+k(3,6);
K(2*j-1,2*n-1)=K(2*j-1,2*n-1)+k(3,7);
K(2*j-1,2*n)=K(2*j-1,2*n)+k(3,8);
K(2*j,2*i-1)=K(2*j,2*i-1)+k(4,1);
K(2*j,2*i)=K(2*j,2*i)+k(4,2);
K(2*j,2*j-1)=K(2*j,2*j-1)+k(4,3);
K(2*j,2*j)=K(2*j,2*j)+k(4,4);
K(2*j,2*m-1)=K(2*j,2*m-1)+k(4,5);
K(2*j,2*m)=K(2*j,2*m)+k(4,6);
K(2*j,2*n-1)=K(2*j,2*n-1)+k(4,7);
K(2*j,2*n)=K(2*j,2*n)+k(4,8);
K(2*m-1,2*i-1)=K(2*m-1,2*i-1)+k(5,1);
K(2*m-1,2*i)=K(2*m-1,2*i)+k(5,2);
K(2*m-1,2*j-1)=K(2*m-1,2*j-1)+k(5,3);
K(2*m-1,2*j)=K(2*m-1,2*j)+k(5,4);
K(2*m-1,2*m-1)=K(2*m-1,2*m-1)+k(5,5);
K(2*m-1,2*m)=K(2*m-1,2*m)+k(5,6);
K(2*m-1,2*n-1)=K(2*m-1,2*n-1)+k(5,7);
K(2*m-1,2*n)=K(2*m-1,2*n)+k(5,8);
K(2*m,2*i-1)=K(2*m,2*i-1)+k(6,1);
K(2*m,2*i)=K(2*m,2*i)+k(6,2);
K(2*m,2*j-1)=K(2*m,2*j-1)+k(6,3);
K(2*m,2*j)=K(2*m,2*j)+k(6,4);
K(2*m,2*m-1)=K(2*m,2*m-1)+k(6,5);
K(2*m,2*m)=K(2*m,2*m)+k(6,6);
K(2*m,2*n-1)=K(2*m,2*n-1)+k(6,7);
K(2*m,2*n)=K(2*m,2*n)+k(6,8);
K(2*n-1,2*i-1)=K(2*n-1,2*i-1)+k(7,1);
K(2*n-1,2*i)=K(2*n-1,2*i)+k(7,2);
K(2*n-1,2*j-1)=K(2*n-1,2*j-1)+k(7,3);
K(2*n-1,2*j)=K(2*n-1,2*j)+k(7,4);
K(2*n-1,2*m-1)=K(2*n-1,2*m-1)+k(7,5);
K(2*n-1,2*m)=K(2*n-1,2*m)+k(7,6);
K(2*n-1,2*n-1)=K(2*n-1,2*n-1)+k(7,7);
K(2*n-1,2*n)=K(2*n-1,2*n)+k(7,8);
K(2*n,2*i-1)=K(2*n,2*i-1)+k(8,1);
K(2*n,2*i)=K(2*n,2*i)+k(8,2);
K(2*n,2*j-1)=K(2*n,2*j-1)+k(8,3);
K(2*n,2*j)=K(2*n,2*j)+k(8,4);
K(2*n,2*m-1)=K(2*n,2*m-1)+k(8,5);
K(2*n,2*m)=K(2*n,2*m)+k(8,6);
K(2*n,2*n-1)=K(2*n,2*n-1)+k(8,7);
K(2*n,2*n)=K(2*n,2*n)+k(8,8);
y=K;
end