This example is to show the rate of convergence of the cubic 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:
We explain degree of freedoms and Lagrange basis for the cubic element on a triangle. There are three types of dofs: vertex type, edge type and element type. Given a mesh, the required data structure can be constructured by
[elem2dof,elem2edge,edge,bdDof,freeDof] = dofP3(elem)
A local basis of P3
node = [0 0; 1 0; 0.5 0.5*sqrt(3)];
elem = [1 2 3];
lambda = [1 0 0; 0 1 0; 0 0 1; ... % 1,2,3 three vertices
0 1/3 2/3; 0 2/3 1/3; ... % 4, 5 first edge
2/3 0 1/3; 1/3 0 2/3; ... % 6, 7 second edge
1/3 2/3 0; 2/3 1/3 0; ... % 8, 9 third edge
1/3 1/3 1/3]; % 10 center of element
dofNode = lambda*node;
figure; subplot(1,2,1);
showmesh(node,elem);
hold on;
findnode(dofNode);
The 10 Lagrange-type basis functions are denoted by $\phi_i, i=1:10$, i.e $\phi_i(x_j)=\delta _{ij},i,j=1:10$. In barycentric coordinates, they are
$$ \phi_1 = 1/2(3\lambda_1-1) (3\lambda_1-2)\lambda_1,\quad \nabla \phi_1 = (27/2 \lambda_1 \lambda_1-9 \lambda_1+1) \nabla \lambda_1,$$$$ \phi_2 = 1/2(3\lambda_2-1) (3\lambda_2-2)\lambda_2,\quad \nabla \phi_2 = (27/2 \lambda_2 \lambda_2-9 \lambda_2+1) \nabla \lambda_2,$$ $$ \phi_3 = 1/2(3\lambda_3-1) (3\lambda_3-2)\lambda_3,\quad \nabla \phi_3 = (27/2 \lambda_3 \lambda_3-9 \lambda_3+1) \nabla \lambda_3,$$ $$ \phi_4 = 9/2\lambda_3\lambda_2 (3\lambda_2-1),\quad \nabla\phi_4 = 9/2 ((3 \lambda_2 \lambda_2-\lambda_2) \nabla \lambda_3+ \lambda_3 (6 \lambda_2-1) \nabla \lambda_2)$$ $$ \phi_5 = 9/2\lambda_3\lambda_2 (3\lambda_3-1),\quad \nabla\phi_5 = 9/2 ((3 \lambda_3 \lambda_3-\lambda_3) \nabla \lambda_2+ \lambda_2 (6 \lambda_3-1) \nabla \lambda_3),$$ $$ \phi_6 = 9/2\lambda_1\lambda_3 (3\lambda_3-1),\quad \nabla\phi_6 = 9/2 ((3 \lambda_3 \lambda_3-\lambda_3) \nabla \lambda_1+ \lambda_1 (6 \lambda_3-1) \nabla \lambda_3),$$$$ \phi_7 = 9/2\lambda_1\lambda_3 (3\lambda_1-1),\quad \nabla\phi_7 = 9/2 ((3 \lambda_1 \lambda_1-\lambda_1) \nabla \lambda_3+ \lambda_3 (6 \lambda_1-1) \nabla \lambda_1),$$ $$ \phi_8 = 9/2\lambda_1\lambda_2 (3\lambda_1-1),\quad \nabla\phi_8 = 9/2 ((3 \lambda_1 \lambda_1-\lambda_1) \nabla \lambda_2+ \lambda_2 (6 \lambda_1-1) \nabla \lambda_1),$$ $$ \phi_9 = 9/2\lambda_1\lambda_2 (3\lambda_2-1),\quad \nabla\phi_9 = 9/2 ((3 \lambda_2 \lambda_2-\lambda_2) \nabla \lambda_1+ \lambda_1 (6 \lambda_2-1) \nabla \lambda_2),$$ $$ \phi_{10} = 27\lambda_1\lambda_2\lambda_3, \quad \nabla \phi_{10} = 27 (\lambda_1 \lambda_2 \nabla \lambda_3+\lambda_1 \lambda_3 \nabla \lambda_2+ \lambda_3 \lambda_2 \nabla \lambda_1).$$When transfer to the reference triangle formed by $(0,0),(1,0),(0,1)$, the local bases in x-y coordinate can be obtained by substituting
$$\lambda _1 = x, \quad \lambda _2 = y, \quad \lambda _3 = 1-x-y.$$The matrix elem2dof
is the local to global index mapping of dofs which can be generated by dofP3(elem)
. The global indices of the dof is organized according to the order of nodes, edges and elements. To be consistent, the dof on an edge depends on the orientation of edge only.
node = [0,0; 1,0; 1,1; 0,1];
elem = [2,3,1; 4,1,3];
figure; subplot(1,2,1);
showmesh(node,elem);
findnode(node);
findelem(node,elem,'all','index','FaceColor',[0.5 0.9 0.45]);
[elem2dof,elem2edge,edge,bdDof,freeDof] = dofP3(elem);
display(elem2dof);
display(bdDof);
%% Setting
[node,elem] = squaremesh([0,1,0,1],0.25);
mesh = struct('node',node,'elem',elem);
option.L0 = 1;
option.maxIt = 4;
option.printlevel = 1;
option.plotflag = 1;
option.elemType = 'P3';
%% Non-empty Dirichlet boundary condition.
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;
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 (3rd order) and L2-norm (4th order) is observed. The order of $\|\nabla u_I - \nabla u_h\|$ is 3rd order and thus no superconvergence exists between nodal interpolate $u_I$ and $u_h$.
MGCG converges uniformly in all cases.