E=1e7;
NU=0;
t=1;
ID=1;
% 1 for plane stress, 2 for plane strain
% Read the nodal coordinate data file.
% first column is x coordinate, second column is y coordinate
%
coordinates=[0 2;0 1;1 1;0 0;1 0;2 0];
figure, plot(coordinates(:,1),coordinates(:,2),'o')
for i1=1:length(coordinates)
text(coordinates(i1,1),coordinates(i1,2),num2str(i1),'fontsize',50)
end
% Read the triangular element data file.
% the node index of each element, contour-clockwise in each element
elements3=[5 2 4;6 3 5;2 5 3;3 1 2];
% Read the Dirichlet boundary condition data file.
dirichlet=[1 3 7 8 10 12];
% A*delta=b, A is global stiffness matrix, delta is grid displacement, b is grid load
A = sparse ( 2*size(coordinates,1), 2*size(coordinates,1) );
b = sparse ( 2*size(coordinates,1), 1 );
A0=A;
for j = 1 : size(elements3,1)
% calculate element stiffness matrix
k=Triangle2D3Node_Stiffness(E,NU,t,coordinates(elements3(j,:),:),1);
% Assembly.
A=A+Triangle2D3Node_Assembly(A0,k,elements3(j,:));
end
% Volume Forces.
b(2)=-1;
% Determine which nodes are associated with Dirichlet conditions.
% Assign the corresponding entries of U, and adjust the right hand side.
u = sparse ( 2*size(coordinates,1), 1 );
BoundNodes = unique ( dirichlet );
u(BoundNodes) = 0;
b = b - A * u;
% Compute the solution by solving A * U = B for the remaining unknown values of U.
FreeNodes = setdiff ( 1:2*size(coordinates,1), BoundNodes );
u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes);
% SHOW displays the solution of the finite element computation..
figure
trisurf ( elements3, coordinates(:,1), coordinates(:,2), (full ( u(1:2:end) ))' );
view ( 0, 90 );
shading interp
title ( 'Solution to the Problem' )
function z = Triangle2D3Node_Assembly(A0,k,index)
% assemble stiffness matrix
% A0: matrix of zeros with size of global stiffness matrix
% k element stiffness matrix
% element nodes indexing index
i=index(1);
j=index(2);
m=index(3);
DOF(1)=2*i-1;
DOF(2)=2*i;
DOF(3)=2*j-1;
DOF(4)=2*j;
DOF(5)=2*m-1;
DOF(6)=2*m;
for n1=1:6
for n2=1:6
A0(DOF(n1),DOF(n2))= A0(DOF(n1),DOF(n2))+k(n1,n2);
end
end
z=A0;
end
function k=Triangle2D3Node_Stiffness(E,NU,t,vertices,ID)
% Parameters:
% Input, real VERTICES(1:3,1:2), contains the 2-dimensional
% coordinates of the vertices.
% Young's module E, poisson ratio NU, thickness t
% ID: 1 for plane stress, 2 for plane strain
%
% Output, real k(1:6,1:6), the local stiffness matrix
% for the element.
xi=vertices(1,1);
yi=vertices(1,2);
xj=vertices(2,1);
yj=vertices(2,2);
xm=vertices(3,1);
ym=vertices(3,2);
A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj))/2;% area of triangle
betai = yj-ym;
betaj = ym-yi;
betam = yi-yj;
gammai = xm-xj;
gammaj = xi-xm;
gammam = xj-xi;
% page 3-10 in shen-lianguan
B = [betai 0 betaj 0 betam 0 ;
0 gammai 0 gammaj 0 gammam ;
gammai betai gammaj betaj gammam betam]/(2*A);
if ID == 1
D = (E/(1-NU*NU))*[1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2];
elseif ID == 2
D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2];
end
k= t*A*B'*D*B;
end