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)^3$$for the following boundary conditions
References:
Subroutines:
- Poisson3P2
- cubePoissonP2
- femPoisson3
- Poisson3P2femrate
The method is implemented in Poisson3P2
subroutine and tested in cubePoissonP2
. Together with other elements (P1, P2,Q1,WG), femPoisson3
provides a concise interface to solve Poisson equation. The P2 element is tested in Poisson3P2femrate
. This doc is based on Poisson3P2femrate
.
For the quadratic element on a tetrahedron, the local basis functions can be written in terms of barycentric coordinates. The 10 dofs is displayed below. The first 4 are associated to the vertices and 6 to the middle points of 6 edges.
Given a mesh, the required data structure can be constructured by
[elem2dof,edge,bdDof] = dof3P2(elem)
The tetrahedron consists of four vertices indexed as [1 2 3 4]. Each tetrahedron contains four faces and six edges. They can be indexed as
locFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2];
locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
In locFace
, the i-th face is opposite to the i-th vertices and the orientation is induced from the tetrahedronthus this is called opposite indexing. In locEdge
, it is the lexicographic indexing which is induced from the lexicographic ordering of the six edges.
node = [1,0,0; 0,1,0; 0,0,0; 0,0,1];
elem = [1 2 3 4];
locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
showmesh3(node,elem);
view(-14,12);
findnode3(node);
findedgedof3(node,locEdge);
When implement boundary conditions, we need to find the d.o.f associated to boundary faces. For example, for Neumann boundary conditions, it can be found by
bdFace2dof = [elem2dof(bdFlag(:,1) == 2,[2,3,4,10,9,8]); ...
elem2dof(bdFlag(:,2) == 2,[1,4,3,10,6,7]); ...
elem2dof(bdFlag(:,3) == 2,[1,2,4,9,7,5]); ...
elem2dof(bdFlag(:,4) == 2,[1,3,2,8,5,6])];
The 10 Lagrange-type bases functions are denoted by $\phi_i, i=1:10$. In barycentric coordinates, they are
$$ \phi_1 = \lambda_1(2\lambda_1 -1),\quad \nabla \phi_1 = \nabla \lambda_1(4\lambda_1-1),$$$$ \phi_2 = \lambda_2(2\lambda_2 -1),\quad \nabla \phi_2 = \nabla \lambda_2(4\lambda_2-1),$$ $$ \phi_3 = \lambda_3(2\lambda_3 -1),\quad \nabla \phi_3 = \nabla \lambda_3(4\lambda_3-1),$$ $$ \phi_4 = \lambda_4(2\lambda_4 -1),\quad \nabla \phi_4 = \nabla \lambda_4(4\lambda_4-1),$$ $$ \phi_5 = 4\lambda_1\lambda_2,\quad \nabla\phi_5 = 4\left (\lambda_1\nabla \lambda_2 + \lambda_2\nabla \lambda_1\right )\; ,$$ $$ \phi_6 = 4\lambda _1\lambda_3,\quad \nabla\phi_6 = 4\left (\lambda_1\nabla \lambda_3 + \lambda_3\nabla \lambda_1\right )\; ,$$ $$ \phi_7 = 4\lambda _1\lambda_4,\quad \nabla\phi_7 = 4\left (\lambda_1\nabla \lambda_4 + \lambda_4\nabla\lambda_1\right )\; .$$$$ \phi_8 = 4\lambda _2\lambda_3,\quad \nabla\phi_8 = 4\left (\lambda_2\nabla \lambda_3 + \lambda_3\nabla \lambda_2\right )\; ,$$ $$ \phi_9 = 4\lambda _2\lambda_4,\quad \nabla\phi_9 = 4\left (\lambda_2\nabla \lambda_4 + \lambda_4\nabla \lambda_2\right )\; ,$$ $$ \phi_{10} = 4\lambda _3\lambda_4,\quad \nabla\phi_{10} = 4\left (\lambda_3\nabla \lambda_4 + \lambda_4\nabla \lambda_3\right )\; .$$When transfer to the reference triangle formed by $(0,0,0),(1,0,0),(0,1,0),(0,0,1)$, the local bases in x-y-z coordinate can be obtained by substituting
$$\lambda _1 = x, \quad \lambda _2 = y, \quad \lambda _3 = z, \quad \lambda_4 = 1-x-y-z.$$Unlike 2-D case, to apply uniform refinement to obtain a fine mesh with good mesh quality, a different ordering of the inital mesh, which may violate the positive ordering, should be used. See 3 D Red Refinement.
%% Setting
[node,elem] = cubemesh([0,1,0,1,0,1],1/3);
mesh = struct('node',node,'elem',elem);
option.maxIt = 4;
option.elemType = 'P2';
option.printlevel = 1;
option.plotflag = 1;
%% Non-empty Dirichlet boundary condition.
fprintf('Mixed boundary conditions. \n');
% pde = polydata3;
pde = sincosdata3;
mesh.bdFlag = setboundary3(node,elem,'Dirichlet','~(x==0)','Neumann','x==0');
femPoisson3(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.
%% Pure Neumann boundary condition.
fprintf('Pure Neumann boundary condition. \n');
pde = sincosdata3Neumann;
option.plotflag = 0;
mesh.bdFlag = setboundary3(node,elem,'Neumann');
femPoisson3(mesh,pde,option);
%% Robin boundary condition.
fprintf('Robin boundary condition. \n');
pde = sincosRobindata3;
mesh.bdFlag = setboundary3(node,elem,'Robin');
femPoisson3(mesh,pde,option);
The optimal rate of convergence of the H1-norm (2nd order) and L2-norm (3rd order) is observed. Superconvergence between discrete functions $\| \nabla u_I - \nabla u_h \|$ is only 2.5 order.
For Neumann problem, the maximum norm is only 2.6 not optimal.
MGCG converges uniformly in all cases.