Matlab code of the Resolution of the 1D Poisson's equation using local Shape Functions

From KratosWiki
Jump to: navigation, search
% Resolution of Poisson 1D using FEM weak form

% Problem definition
x0=0.0;
xL=15.0;
Nx=101;
fi0=3;   % Dirichlet condition
qL=13;   % Neumann condition
Q0=5;    % Heat load
km=1;     % material

% definition of nodes and elements
long=xL-x0;              % length of domain
nnodes=5;                % number of total nodes
nnod=2;                  % number of nodes for each element
nelems=nnodes-1;         % number of elements
interv=long/(nnodes-1);  % size of each element
xi=zeros(nnodes,1);
xi(1)=x0;
for i=2:nnodes
    xi(i)=xi(i-1)+interv;
end


f=zeros(nelems,nnod);
K=zeros(nelems,nnod,nnod);
% obtaining elementary matrix components
for i=1:nelems
    for j=1:nnod
        f(i,j)=Q0*interv/2;
        for k=1:nnod
            K(i,j,k)=((-1)^(j+k))*(km/interv);
        end;
    end;
end;
        
% assembly of global matrix components
fg=zeros(nnodes,1);
Kg=zeros(nnodes,nnodes);
for i=1:nelems
    for j=1:nnod
        fg(i+j-1)=fg(i+j-1)+f(i,j);
        for k=1:nnod
            Kg(i+j-1,i+k-1)=Kg(i+j-1,i+k-1)+K(i,j,k);
        end
    end
end
        
% assembly of boundary conditions
fg(nnodes)=fg(nnodes)-qL;

% Applying Dirichlet condition 
% (by removing the equation related to the first unknown 
% and by using its value for the independent term)
Kf=Kg(2:length(Kg),2:length(Kg));
ff=fg(2:length(fg))-Kg(2:length(Kg),1)*fi0;

fi=inv(Kf)*ff;

fi=[fi0;fi];

% exact solution
realFi=-(Q0/(2*km))*(xi.^2)+(-qL+Q0*(xL-x0))*xi/km+fi0;

% computation of the error for the nodal values
nodalerr=realFi-fi;



Download the Source code
Back to the theory
Personal tools
Categories