% Finite element analysis of linear plane-stress continuum
clear all
close all
clc
global B_cell Det_cell
%
% Pre-processing
%
% Read the input from a file; change the file name
fid = fopen('elem50.dat','r');
% Element dta
fgets(fid);
fgets(fid);
Nelem = fscanf(fid,'%i\n',1);
fgets(fid);
fgets(fid);
elemdata = fscanf(fid,'%i\n',[6,Nelem]);
elemdata = elemdata';
% Material property data
fgets(fid);
fgets(fid);
NmaterialID = 1;
WGT = [1.0 1.0 1.0 1.0];
for i = 1:Nelem
if(NmaterialID < elemdata(i,2))
NmaterialID = elemdata(i,2);
end
end
matprop = fscanf(fid,'%g\n',[4,NmaterialID]);
matprop = matprop';
% Nodal coordinates
fgets(fid);
Nnode = fscanf(fid,'%i\n',1);
fgets(fid);
fgets(fid);
nodedata = fscanf(fid,'%i %g %g\n',[3, Nnode]);
nodedata = nodedata';
% Displacement boundary conditions
fgets(fid);
Ndispbc = fscanf(fid,'%i\n',1);
fgets(fid);
dispbc = fscanf(fid,'%i %g %g\n',[3, Ndispbc]);
dispbc = dispbc';
% Force boundary conditions
fgets(fid);
Nforcebc = fscanf(fid,'%i\n',1);
fgets(fid);
forcebc = fscanf(fid,'%i %g %g\n',[3, Nforcebc]);
forcebc = forcebc';
fclose(fid);
%
% Analysis
%
% Calling Global stiffness matrix
calc_Bdet(elemdata,nodedata);
Ksing = Kglobal(matprop,nodedata,elemdata);
K = Ksing;
%
% Arrange force information into a force vector, F
F = zeros(2*Nnode,1); % Initialization
% Nforce = size(forcebc,1);
% for i = 1:Nforce,
% dof = (forcebc(i,1)-1)*2 + forcebc(i,2);
% F(dof) = forcebc(i,3);
% end
sigma0=[10e9;0;0];
for i=1:Nelem
dofs=[2*elemdata(i,3)-1 2*elemdata(i,3) 2*elemdata(i,4)-1 2*elemdata(i,4) 2*elemdata(i,5)-1 2*elemdata(i,5) 2*elemdata(i,6)-1 2*elemdata(i,6)];
for j=1:4
WT = WGT(j);
F(dofs) = F(dofs) + Det_cell{i,j}*matprop(elemdata(i,2),4)*(B_cell{i,j})'*sigma0;
end
end
%
% Displacement boundary conditions
Nfix = size(dispbc,1);
j = 0;
dispID=zeros(Nfix,1);
dispVal=zeros(Nfix,1);
for i = 1:Nfix,
dispID(i) = (dispbc(i,1)-1)*2+dispbc(i,2);
dispVal(i) = dispbc(i,3);
end
[dispID,sortIndex] = sort(dispID);
dispVal = dispVal(sortIndex);
% Imposing the boundary conditions
% ------------------------------------------
% dispID array contains the dof which are assigned specified values.
[sm,sn] = size(dispID);
Ndbc = sm;
for nd=1:Ndbc,
for nr=1:2*Nnode-nd+1,
F(nr) = F(nr) - K(nr,dispID(nd)-nd+1) * dispVal(nd);
end
K = matcut(K,dispID(nd)-nd+1);
F = veccut(F,dispID(nd)-nd+1);
end
% To solve for unknown displacements.
U = K\F;
%
% Compile results and the specified displacements into one vector
%
u = zeros(2*Nnode,1);
%
for iu=1:Ndbc,
u(dispID(iu)) = 12345.12345;
end
%
iuc = 0;
for iu=1:2*Nnode,
if u(iu) == 12345.12345
iuc = iuc+1;
else
u(iu) = U(iu-iuc);
end
end
%
for iu=1:Ndbc,
u(dispID(iu)) = dispVal(iu);
end
%
% Post-processing
%
% Plotting the undeformed and deformed configurations
figure(1)
clf
x0 = nodedata(:,2); y0 = nodedata (:,3);
xd = x0 + u(1:2:2*Nnode); yd = y0 + u(2:2:2*Nnode);
axis([min([x0; xd]) max([x0; xd]) min([y0; yd]) max([y0; yd])]);
for i = 1:Nelem
hndl = fill(x0(elemdata(i,3:6)), y0(elemdata(i,3:6)),'w');
hold on
end
for i = 1:Nelem
hndl = fill(xd(elemdata(i,3:6)), yd(elemdata(i,3:6)),'c');
end
axis equal
title('Beam deformation due to residual stresses');