www.pudn.com > 87361026FEA.rar > fem_example.m, change:2006-02-28,size:1286b


% fem_example.m
%
% Jack Chessa 
% Northwestern University
%
% A simple 1D finite element program in Matlab
%
%  /|
%  /|============================ 
%  /|           E,A,L           | ------> P                                 
%  /|============================
%  /|


% define problem parameters
A=1;                    % cross sectional area
E=10;                   % Young's modulus
L=5;                    % length of rod
P=3;                    % tip load

% define mesh parameters
nn=21;                  % number of nodes
ne=nn-1;                % number of elements
l=L/ne;                 % length of each element

% define mesh
node=linspace(0,L,nn)';         % node coordinates
element=[ 1:(nn-1); 2:nn ]';    % element connectivity

% define fe operators
K=sparse(nn,nn);        % stiffness matrix
f=zeros(nn,1);          % force vector
d=zeros(nn,1);          % nodal displacement

ke=A*E/l*[ 1 -1;
          -1  1 ];      % stiffness matrix for an element
        
% define force vector
f(nn)=P;

% compute stiffness matrix
for e=1:ne
  
  sctr = element(e,:);  
  K(sctr,sctr) = ...
                    K(sctr,sctr) + ke;
  
end

% solve system
freeDof = 2:nn;             % free dofs
d(freeDof) = K(freeDof,freeDof)\f(freeDof);


% plot results
clf
plot(node,d,'r-o')