FEM2D: Poissson equation with Triangular elements

Consider the 2D equation model, a second order PDE equation with general coefficients a11, a12, a21, a22, a00, f :

Contents

For Linear Triangular Elements we know that the formulas are explicit and they can be stated as

Poisson Example:

Consider the solution of the Poisson's 2D equation. Corresponding to the 2D model PDE equation with a11=a22=1, a12=a21=a00=0, f=1: We are solving the problem introduced on the document T4-MN-FEM2D-Applications, pg. 2. Using symmetry and 4 triangular elements. We will implement the general case for the model PDE equation and use it to solve the Poisson Equation

Define Geometry: node coordinates and elements

clear all;
close all;
nodes=[
    0,0;
    0.5,0;
    0.5,0.5;
    1,0;
    1,0.5;
    1,1;
    ];
elem=[1,2,3;
      5,3,2;
      2,4,5;
      3,5,6;
      ];
numNod=size(nodes,1)
numElem=size(elem,1)
numbering=1; %if you want to see number values
plotElements(nodes,elem,numbering);

numNod =

     6


numElem =

     4

Define Coefficients vector of the model equation

In this case we use the Poisson coefficients defined in the problem above

a11=1;
a12=0;
a21=a12;
a22=a11;
a00=0;
f=1;
coeff=[a11,a12,a21,a22,a00,f];

Compute the global stiff matrix

Compute the global stiff matrix

K=zeros(numNod);   %global Stiff Matrix
F=zeros(numNod,1); %global internal forces vector
Q=zeros(numNod,1); %global secondary variables vector
for e=1:numElem
    %
    % First compute Ke, Fe for each element
    %
    v1=nodes(elem(e,1),:);
    v2=nodes(elem(e,2),:);
    v3=nodes(elem(e,3),:);
    Area=0.5*det([1,v1;1,v2;1,v3]);
    beta(1)=v2(2)-v3(2);
    gamma(1)=v3(1)-v2(1);
    % cyclic permutation
    beta(2)=v3(2)-v1(2);
    gamma(2)=v1(1)-v3(1);
    beta(3)=v1(2)-v2(2);
    gamma(3)=v2(1)-v1(1);
    %
    % element stiff matrix
    %
    Ke=zeros(3);
    Fe=zeros(3,1);
    if (coeff(1) ~= 0) %a11
        for i=1:3
            for j=1:3
                K11(i,j)=beta(i)*beta(j);
            end
        end
        Ke=Ke+coeff(1)*K11/(4*Area);
    end
    if (coeff(2) ~= 0) %a12
        for i=1:3
            for j=1:3
                K12(i,j)=beta(i)*gamma(j);
            end
        end
        Ke=Ke+coeff(2)*K12/(4*Area);
    end
    if (coeff(3) ~= 0) %a21
        for i=1:3
            for j=1:3
                K21(i,j)=gamma(i)*beta(j);
            end
        end
        Ke=Ke+coeff(3)*K21/(4*Area);
    end
    if (coeff(4) ~= 0) %a22
        for i=1:3
            for j=1:3
                K22(i,j)=gamma(i)*gamma(j);
            end
        end
        Ke=Ke+coeff(4)*K22/(4*Area);
    end
    if (coeff(5) ~= 0) %a00
       K00=Area/12*[2,1,1;1,2,1;1,1,2];
       Ke=Ke+coeff(5)*K00;
    end
    if (coeff(6) ~= 0) %f
        Fe=(coeff(6)*Area/3)*[1;1;1];
    end
  %
  % Assemble the elements
  %
    rows=[elem(e,1); elem(e,2); elem(e,3)];
    colums= rows;
    K(rows,colums)=K(rows,colums)+Ke; %assembly
    if (coeff(6) ~= 0)
        F(rows)=F(rows)+Fe;
    end
end % end for elements
% we save a copy of the initial F array
% for the postprocess step
Fini=F;
%

Boundary Conditions (BC)

Impose Boundary Conditions for this example. In this case only essential and natural BC are considered. We will do a thermal example for convection BC (mixed).

fixedNodes=[4,5,6]; %Fixed Nodes (global num.)
freeNodes = setdiff(1:numNod,fixedNodes); %Complementary of fixed nodes
% ------------ Essential BC
u=zeros(numNod,1); %initialize u vector
u(fixedNodes)=0; %all of them are zero
F=F-K(:,fixedNodes)*u(fixedNodes);%here u can be different from zero only for fixed nodes
% ------------ Natural BC
Q(freeNodes)=0; %all of them are zero
% modify the linear system, this is only valid if BC = 0.
Km=K(freeNodes,freeNodes);
Im=F(freeNodes)+Q(freeNodes);

Compute the solution

solve the System

format short e; %just to a better view of the numbers
um=Km\Im;
u(freeNodes)=um;
u(fixedNodes)=0;
u

u =

   3.1250e-01
   2.2917e-01
   1.7708e-01
            0
            0
            0

PostProcess: Compute secondary variables and plot results

Q=K*u-Fini
titol='Poisson solution';
colorScale='jet';
plotContourSolution(nodes,elem,u,titol,colorScale)

Q =

  -6.9389e-18
   1.1102e-16
            0
  -1.5625e-01
  -3.0208e-01
  -4.1667e-02

Exercise:

For triangular elements, implement a Matlab function that returns the element stiff matrix, Ke, and the Fe vector for each element depending on the coeff vector for the model equation.

function [Ke,Fe]=linearTriangElement(coeff,nodes,elem,e)

coeff --> coefficient vector = [a11,a12,a21,a22,a00,f] for the model equation

nodes --> matrix with the coordinates of the nodes

elem --> connectivity matrix defining the elements

e --> number of the present element

(c)Numerical Factory 2017

Published with MATLAB® R2018a