This example is to show the rate of convergence of the linear finite element approximation of the Poisson equation on the unit square:
$$- \Delta u = f \; \hbox{in } (0,1)^2$$for the following boundary conditions
References:
Subroutines:
- Poisson
- squarePoisson
- femPoisson
- Poissonfemrate
The method is implemented in Poisson
subroutine and tested in squarePoisson
. Together with other elements (P1, P2, P3, Q1), femPoisson
provides a concise interface to solve Poisson equation. The P1 element is tested in Poissonfemrate
. This doc is based on Poissonfemrate
.
For the linear element on a simplex, the local basis functions are
barycentric coordinate of vertices. The local to global pointer is
elem
. This is the simplest and default element for elliptic equations.
A local basis of P1
For $i = 1, 2, 3$, a local basis of the linear element space is given by the barycentric coordinate
$$\phi_i = \lambda_i, \quad \nabla \phi_i = \nabla \lambda_i = - \frac{|e_i|}{2|T|}\boldsymbol n_i, \quad \int_T \nabla \phi_i\nabla \phi_j = \frac{1}{4|T|}\boldsymbol l_i \cdot \boldsymbol l_j$$where $e_i$ is the edge opposite to the i-th vertex and $\boldsymbol n_i$ is the unit outwards normal direction, and $\boldsymbol l_i$ is the edge vector of $e_i$.
See Finite Element Methods Section 2.1 for geometric explanation of the barycentric coordinate and Programming of Finite Element Methods in MATLAB for detailed explanation. For P1 element, the basic data structure node,elem
is sufficient and displayed for the following coarse mesh.
imatlab_export_fig('print-jpeg')
node = [0,0; 1,0; 1,1; 0,1];
elem = [2,3,1; 4,1,3];
[node,elem] = uniformbisect(node,elem);
figure('rend','painters','pos',[10 10 225 225])
showmesh(node,elem);
findnode(node);
findelem(node,elem);
display(elem);
% Setting
[node,elem] = squaremesh([0,1,0,1],0.25);
mesh = struct('node',node,'elem',elem);
figure('rend','painters','pos',[10 10 225 225])
showmesh(node,elem);
option.L0 = 3;
option.maxIt = 4;
option.printlevel = 1;
option.plotflag = 0;
pde = sincosdata;
mesh.bdFlag = setboundary(node,elem,'Dirichlet','~(x==0)','Neumann','x==0');
femPoisson(mesh,pde,option);
When pure Neumann boundary condition is posed, i.e., $-\Delta u =f$ in $\Omega$ and $\nabla u\cdot n=g_N$ on $\partial \Omega$, the data should be consisitent in the sense that $\int_{\Omega} f \, dx + \int_{\partial \Omega} g \, ds = 0$. The solution is unique up to a constant. A post-process is applied such that the constraint $\int_{\Omega}u_h dx = 0$ is imposed.
option.plotflag = 0;
pde = sincosNeumanndata;
pde = sincosdata;
mesh.bdFlag = setboundary(node,elem,'Neumann');
femPoisson(mesh,pde,option);
option.plotflag = 0;
pde = sincosRobindata;
mesh.bdFlag = setboundary(node,elem,'Robin');
femPoisson(mesh,pde,option);
The optimal rate of convergence of the H1-norm (1st order) and L2-norm (2nd order) is observed. The 2nd order convergent rate between two discrete functions $\|\nabla (u_I - u_h)\|$ is known as superconvergence.
MGCG converges uniformly in all cases.