% Torsion on stepped bar, Torq2.m
clear                        % variables
format short e               % E format
torq2 = -10e3                % moment at node 2

gx = [ 11 11 ]*1e6          % shear modulus
len = [4 3]*12              % element lengths
dia = [1.25  1.5]           % element diameters
j = pi/32*dia .^4;          % vector of polar moment
k12 = j(1)*gx(1)/len(1)     % stiff const
k23 = j(2)*gx(2)/len(2)
k22 = k12 + k23

beta = k12*1e7              % zienk factor
stiff = [ k12    -k12     0    % stiffness matrix
         -k12     k22   -k23
          0      -k23    k23 ]
df =  [0*beta; torq2;  0*beta]; % combined d, F
zienk = stiff;
zienk(1,1)=zienk(1,1)+beta;
zienk(3,3)=zienk(3,3)+beta

rotat = inv(zienk) * df;     % displacements
torque = stiff * rotat;      % force vector

for i = 1:length(df)-1       % element displacements
  eldis(i) = rotat(i+1) - rotat(i);
end
stress = gx .* dia/2 .* eldis ./ len;   % element stresses

fprintf('\n      Torque2 ')
fprintf('\n   Your name, MENG 421, ')
disp(date)

fprintf('\n     Node   Rotation      Torque ')
fprintf('\n            Radian        Kip-in \n')
for i = 1:length(df)
  fprintf('    %3.0f   %11.2e   %7.1f \n',i,rotat(i),torque(i)/1000)
end

fprintf('\n     Elem   Stress ')
fprintf('\n             Ksi \n')
for i = 1:length(df)-1
  fprintf('    %3.0f   %7.2f \n',i,stress(i)/1000)
end

%        Torque2 
%   Your name, MENG 421, 02-Mar-2004

%      Node   Rotation      Torque 
%             Radian        Kip-in 
%       1    -4.84e-009       2.7 
%       2    -4.84e-002     -10.0 
%       3    -1.34e-008       7.3 

%      Elem   Stress 
%              Ksi 
%       1     -6.93 
%       2     11.08