MATLAB Examples

Poisson equation in a square

This code solves a test problem involving a Poisson equation on a square domain. The method relies on Lagrangian finite elements on a uniform triangular mesh. The solver is documented in the PhD thesis of Sebastian Ullmann "POD-Galerkin Modeling for Incompressible Flows with Stochastic Boundary Conditions". Please reference the thesis or one of the related papers if you use the solver or parts of it for your own scientific work.

Copyright (C) Sebastian Ullmann 2015


Clean up storage

Warning: Objects of 'onCleanup' class exist.  Cannot clear this class or any of
its superclasses. 

Create finite element solver object

This object creates and stores everything needed for a Poisson problem with the given mesh and boundary conditions, in particular,

  • the finite element mesh,
  • the boundary conditions,
  • logical index arrays to easily access different parts of the solution,
  • the finite element matrices.

However, boundary and source data are not yet prescribed.

solverOptions = struct();
solverOptions.Mesh = @meshRectangularDiagonalDown;
solverOptions.ConditionsAtEdges = @conditionsAtEdges;

fem = SolverPoisson(solverOptions);

Inspect the finite element solver object

  SolverPoisson with properties:

                   M: [441x441 double]
                   C: [441x441 double]
                   N: [441x441 double]
         nComponents: 1
          components: [1x1 struct]
    polynomialDegree: 1
                mesh: [1x1 TriangularTaylorHood]
                pick: [1x1 struct]
         periodicMap: [441x441 double]
           equations: [1x1 struct]
              points: [441x3 double]

Inspect the finite element mesh object

Within the solver object, the property mesh contains the details of the spatial discretization and assignment of boundary conditions.

  TriangularTaylorHood with properties:

                 points: [441x3 double]
              triangles: [200x9 double]
                  edges: [320x4 double]
       useAssemblyCache: 1
    useAssemblyOrdering: 1
                    phi: {3x3 cell}
           triangleArea: [200x1 double]
             edgeLength: [320x1 double]
           mixedProduct: [200x1 double]
      quadratureWeights: [1x1 struct]
          assemblyCache: [1x1 struct]

Plot the finite element mesh.

We extract the data from fem.mesh. It contains the following information:


  • Column 1 contains the x locations of the mesh points,
  • Column 2 contains the y locations of the mesh points,
  • Column 3 contains point conditions.


  • Column 1 contains the index of vertex A of each triangle,
  • Column 2 contains the index of vertex B of each triangle,
  • Column 3 contains the index of vertex C of each triangle,
  • Column 4 contains the index of the edge node opposite of vertex A,
  • Column 5 contains the index of the edge node opposite of vertex B,
  • Column 6 contains the index of the edge node opposite of vertex C.


  • Column 1 contains the index of point A of each edge,
  • Column 2 contains the index of point B of each edge,
  • Column 3 contains the index of each edge midpoint,
  • Column 4 contains edge conditions.
points    = fem.mesh.points;
triangles = fem.mesh.triangles;
edges     = fem.mesh.edges;

trimesh(triangles(:,1:3), ...
        points(:,1), ...
        points(:,2), ...
        zeros(size(points,1),1), ...
hold on

Solve the problem

The solution is stored as a field of the structure array matObj. The function handle prescribed as OutputFunction defines how this is done. The function handle prescribed as SourceFunction assigns the values of the source to the mesh nodes. This method, which interpolates the source via the finite element basis function, is exact for source functions which are quadratic polynomials. Higher-order quadrature is not yet available.

Instead of storing the solution to a struct, one can directly write to disk by using matfile.

matObj = struct();
matObj = fem.solve(matObj, ...
                  'OutputFunction',@(matObj,U)fem.outputSolution(matObj,U,'U'), ...
                  'DirichletFunction',@(f)fem.dirichletUnity(f), ...
                  'NeumannFunction',@(f)fem.neumannUnity(f), ...

Plot the solution

trisurf(fem.mesh.triangles(:,1:3), ...
        fem.mesh.points(:,1), ...
        fem.mesh.points(:,2), ...
        matObj.U, ...