How to change Matlab program for solving equation with finite element method?
- by DSblizzard
I don't know is this question more related to mathematics or programming and I'm absolute newbie in Matlab.
Program FEM_50 applies the finite element method to Laplace's equation -Uxx(x, y) - Uyy(x, y) = F(x, y) in Omega.
How to change it to apply FEM to equation -Uxx(x, y) - Uyy(x, y) + U(x, y) = F(x, y)?
At this page:
http://sc.fsu.edu/~burkardt/m_src/fem_50/fem_50.html
additional code files in case you need them.
function fem_50 ( )
%% FEM_50 applies the finite element method to Laplace's equation.
%
% Discussion:
%
% FEM_50 is a set of MATLAB routines to apply the finite
% element method to solving Laplace's equation in an arbitrary
% region, using about 50 lines of MATLAB code.
%
% FEM_50 is partly a demonstration, to show how little it
% takes to implement the finite element method (at least using
% every possible MATLAB shortcut.) The user supplies datafiles
% that specify the geometry of the region and its arrangement
% into triangular and quadrilateral elements, and the location
% and type of the boundary conditions, which can be any mixture
% of Neumann and Dirichlet.
%
% The unknown state variable U(x,y) is assumed to satisfy
% Laplace's equation:
% -Uxx(x,y) - Uyy(x,y) = F(x,y) in Omega
% with Dirichlet boundary conditions
% U(x,y) = U_D(x,y) on Gamma_D
% and Neumann boundary conditions on the outward normal derivative:
% Un(x,y) = G(x,y) on Gamma_N
% If Gamma designates the boundary of the region Omega,
% then we presume that
% Gamma = Gamma_D + Gamma_N
% but the user is free to determine which boundary conditions to
% apply. Note, however, that the problem will generally be singular
% unless at least one Dirichlet boundary condition is specified.
%
% The code uses piecewise linear basis functions for triangular elements,
% and piecewise isoparametric bilinear basis functions for quadrilateral
% 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, and
% the right hand side function in Laplace's equation. Note that the
% fact that the geometry is completely up to the user means that
% just about any two dimensional region can be handled, with arbitrary
% shape, including holes and islands.
%
clear
%
% Read the nodal coordinate data file.
%
load coordinates.dat;
%
% Read the triangular element data file.
%
load elements3.dat;
%
% Read the quadrilateral element data file.
%
load elements4.dat;
%
% Read the Neumann boundary condition data file.
% I THINK the purpose of the EVAL command is to create an empty NEUMANN array
% if no Neumann file is found.
%
eval ( 'load neumann.dat;', 'neumann=[];' );
%
% Read the Dirichlet boundary condition data file.
%
load dirichlet.dat;
A = sparse ( size(coordinates,1), size(coordinates,1) );
b = sparse ( size(coordinates,1), 1 );
%
% Assembly.
%
for j = 1 : size(elements3,1)
A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ...
+ stima3(coordinates(elements3(j,:),:));
end
for j = 1 : size(elements4,1)
A(elements4(j,:),elements4(j,:)) = A(elements4(j,:),elements4(j,:)) ...
+ stima4(coordinates(elements4(j,:),:));
end
%
% Volume Forces.
%
for j = 1 : size(elements3,1)
b(elements3(j,:)) = b(elements3(j,:)) ...
+ det( [1,1,1; coordinates(elements3(j,:),:)'] ) * ...
f(sum(coordinates(elements3(j,:),:))/3)/6;
end
for j = 1 : size(elements4,1)
b(elements4(j,:)) = b(elements4(j,:)) ...
+ det([1,1,1; coordinates(elements4(j,1:3),:)'] ) * ...
f(sum(coordinates(elements4(j,:),:))/4)/4;
end
%
% Neumann conditions.
%
if ( ~isempty(neumann) )
for j = 1 : size(neumann,1)
b(neumann(j,:)) = b(neumann(j,:)) + ...
norm(coordinates(neumann(j,1),:) - coordinates(neumann(j,2),:)) * ...
g(sum(coordinates(neumann(j,:),:))/2)/2;
end
end
%
% Determine which nodes are associated with Dirichlet conditions.
% Assign the corresponding entries of U, and adjust the right hand side.
%
u = sparse ( size(coordinates,1), 1 );
BoundNodes = unique ( dirichlet );
u(BoundNodes) = u_d ( coordinates(BoundNodes,:) );
b = b - A * u;
%
% Compute the solution by solving A * U = B for the remaining unknown values of U.
%
FreeNodes = setdiff ( 1:size(coordinates,1), BoundNodes );
u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes);
%
% Graphic representation.
%
show ( elements3, elements4, coordinates, full ( u ) );
return
end
