本程序包含主程序：
fem.m
两个子程序：
Triangle2D3Node_Stiffness.m
Triangle2D3Node_Assembly.m
及坐标数据等文件：
coordinates.dat
elements3.dat
dirichlet.dat
下面分别进行介绍：
-------------------------------------------------------------------------------------------------------
主程序fem.m
%*****************************************************************************80
%
% a set of MATLAB routines to apply the finite
% element method to solving example in page 3-33 in shen_lianguan
%
% The code uses piecewise linear basis functions for triangular elements,
%
% The user is required to supply a number of data files and MATLAB
% functions that specify the location of nodes, the grouping of nodes
% into elements, the location and value of boundary conditions
%
clear
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
%
load coordinates.dat;
%%{
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
load elements3.dat;
%
% Read the Dirichlet boundary condition data file.
%
load dirichlet.dat;
% 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\' )
%
-------------------------------------------------------------------------------------------------------
Triangle2D3Node_Stiffness.m
%-----------------------------------------------------------
function k=Triangle2D3Node_Stiffness(E,NU,t,vertices,ID)
%% STIMA3 determines the local stiffness matrix for a triangular element.
%
%
% 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;
-------------------------------------------------------------------------------------------------------
Triangle2D3Node_Assembly.m
%-----------------------------------------------------------
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;
-------------------------------------------------------------------------------------------------------
数据文件：
coordinates.dat
0 2
0 1
1 1
0 0
1 0
2 0
elements3.dat
5 2 4
6 3 5
2 5 3
3 1 2
dirichlet.dat
1
3
7
8
10
12
-------------------------------------------------------------------------------------------------------