226fall 2009project 1

Linear Finite Element Methods

The purpose of this project is to write a finite element code for solving the Poisson equation in a general polygonal domain using piecewise linear finite elements. The project will be divided by several steps.

Step 1. Mesh

• load data file mesh.mat in MATLAB
• use the subroutine showmesh.m to plot the mesh
• use findnode.m to show the index of the nodes
• use findelem.m to show the index of the elements

You will end with this mesh.

Step 2. Local Stiffness matrices

Compute all the local stiffness matrix for the Poisson equation discretized by linear finite elements defined on the mesh given above.

You can compare your result with the answer.

Step 3. Global Stiffness matrix

Compute the global stiffness matrix based on the local stiffness matrices and store the stiffness matrix using the sparse matrix data struture.

You can compare your result with the answer where the matrix is presented in the formate

i   j    A(i,j)

Step 4. Right hand side

Using three points quadrature (i.e. 3 middle points of a triangle) to compute the right hand side vector b.

Step 5. Boundary conditions

Using the function findboundary.m to find boundary nodes and edges. Code Dirichlet boundary condition and Neumann boundary condition by changing the right hand side.

When solving Poinsson equation with homogenous Neumann boundary condition, the resulting right hand side b should satisfy the compatible condition

 sum(b) = 0 

Because of the limited accuracy of the numerical quadrature, the computed right hand side vector does not exactly satisfy the compatible condition. One may slightly modify the right hand in such a way

b = b - mean(b);

Step 6. Convergent rate

• Solve Ax=b by using direct method in Matlab. When solving pure Neumman boundary conditions, the matrix is singular, you can fix one point.
• Use showsolution.m to plot the solution.
• Use uniformrefine.m to refine the grid and compute a sequence of solutions.
• Choose a smooth solution, say u = sin(2*pi*x)*cos(2*pi*y), and compute the error (for two types of boundary conditions)

  || D u_I - D u_h ||   
  

  and conclude the rate of convergence.