This example is to show the rate of convergence of the linear finite element approximation of the Poisson equation on the unit square:
$$- \nabla (d \nabla) u = f \; \hbox{in } (0,1)^2$$for the following boundary conditions
Find $(\sigma , u)$ in $H_{g_N,\Gamma_N}(div,\Omega)\times L^2(\Omega)$ s.t.
$$ (d^{-1}\sigma,\tau) + (div \tau, u) = \langle \tau \cdot n, g_D \rangle_{\Gamma_D} \quad \forall \tau \in H_{0,\Gamma_N}(div,\Omega)$$$$ (div \sigma, v) = -(f,v) \quad \forall v \in L^2(\Omega) $$where
$$H_{g,\Gamma}(div,\Omega) = \{\sigma \in H(div,\Omega); \sigma \cdot n = g \text{ on } \Gamma \subset \partial\Omega \}.$$
The unknown $\sigma = d\nabla u$ is approximated using the lowest order Raviart-Thomas element (RT0) and $u$ by piecewise constant element (P0).
References
Subroutines:
- PoissonRT0
- squarePoissonRT0
- mfemPoisson
- PoissonRT0mfemrate
The method is implemented in PoissonRT0
subroutine and tested in squarePoissonRT0
. Together with other elements (BDM1), mfemPoisson
provides a concise interface to solve Poisson equation in mixed formulation. The RT0-P0 element is tested in PoissonRT0mfemrate
. This doc is based on PoissonRT0mfemrate
.
We explain degree of freedoms and basis functions for Raviart-Thomas element on a triangle.
The dofs and basis depends on the orientation of the mesh. We shall use the asecond orientation, i.e., elem(t,1)< elem(t,2)< elem(t,3)
not the positive orientation. Given an elem
, the asecond orientation can be constructed by
[elem,bdFlag] = sortelem(elem,bdFlag); % ascend ordering
Note that bdFlag
should be sorted as well.
The local edge is also asecond [2 3; 1 3; 1 2]
so that the local orientation is consistent with the global one and thus no need to deal with the sign difference when the positive oritentation is used. Read Simplicial complex in two dimensions for more discussion of indexing, ordering and orientation.
node = [1,0; 1,1; 0,0];
elem = [1 2 3];
edge = [2 3; 1 3; 1 2];
figure;
subplot(1,2,1)
showmesh(node,elem);
findnode(node);
findedge(node,edge,'all','rotvec');
Suppose [i,j]
is the k-th edge. The two dimensional curl is a rotated graident defined as $\nabla^{\bot} f = (-\partial_y f, \partial _x f).$ The basis of this edge along with its divergence are given by
Inside one triangular, the 3 bases corresponding to 3 local edges [2 3; 1 3; 1 2] are:
$$ \phi_1 = \lambda_2 \nabla^{\bot} \lambda_3 - \lambda_3 \nabla^{\bot} \lambda_2. $$ $$ \phi_2 = \lambda_1 \nabla^{\bot} \lambda_3 - \lambda_3 \nabla^{\bot} \lambda_1. $$$$ \phi_3 = \lambda_1 \nabla^{\bot} \lambda_2 - \lambda_2 \nabla^{\bot} \lambda_1. $$The dual basis is the line integral over an orientated edge
$$\int_{e_i} \phi_j \cdot n_i \, ds = \delta(i,j),$$where $n_i = t_i^{\bot}$ is the rotation of the unit tangential vector of $e_i$ by $90^{\deg}$ counterclockwise.
Three local edges are locEdge = [2 3; 1 3; 1 2]
. The pointer from the local to global index can be constructured by
[elem2dof,edge] = dofedge(elem);
[node,elem] = squaremesh([0 1 0 1], 0.5);
bdFlag = setboundary(node,elem,'Dirichlet');
[elem,bdFlag] = sortelem(elem,bdFlag);
showmesh(node,elem);
findnode(node);
findelem(node,elem);
[elem2dof,edge] = dofedge(elem);
findedge(node,edge,'all','rotvec');
display(elem2dof);
The ascend ordering orientation is not consistent with the induced orientation. The second edge would be [3 1]
for the consistent orientation. So [1 -1 1]
is used in the construction of div operator.
For triangle t, the basis for the constant function space is $p = 1$, the characteristic function. So in the computation of divergence operator, elemSign
should be used to correct the sign. In the output of gradbasis
, -Dlambda
is always the outwards normal direction. The signed area could be negative but in the ouput, area
is the absolute value (for the easy of integration on elements) and elemSign
is used to record elements with negative area.
[Dlambda,area,elemSign] = gradbasis(node,elem);
B = icdmat(double(elem2dof),elemSign*[1 -1 1]);
display(full(B))
Direction of boundary edges may not be the outwards normal direction of the domain since now elem
is ascend orientation. edgeSign
is introduced to record this inconsistency.
edgeSign = ones(NE,1);
idx = (bdFlag(:,1) ~= 0) & (elemSign == -1); % first edge is on boundary
edgeSign(elem2edge(idx,1)) = -1;
idx = (bdFlag(:,2) ~= 0) & (elemSign == 1); % second edge is on boundary
edgeSign(elem2edge(idx,2)) = -1;
idx = (bdFlag(:,3) ~= 0) & (elemSign == -1); % third edge is on boundary
edgeSign(elem2edge(idx,3)) = -1;
%% 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.elemType = 'RT0';
pde = sincosNeumanndata;
%% Mix Dirichlet and Neumann boundary condition.
option.solver = 'uzawapcg';
mesh.bdFlag = setboundary(node,elem,'Dirichlet','~(x==0)','Neumann','x==0');
mfemPoisson(mesh,pde,option);
In mixed formulation,the Neumann boundary condition for $d\nabla u$ becomes the Dirichlet boundary condiiton for $\sigma$. The space for $u$ is $L^2_0$ and thus one dof should be removed to have a non-singular system.
option.plotflag = 0;
mesh.bdFlag = setboundary(node,elem,'Neumann');
mfemPoisson(mesh,pde,option);
%% Pure Dirichlet boundary condition.
mesh.bdFlag = setboundary(node,elem,'Dirichlet');
mfemPoisson(mesh,pde,option);
The optimal rates of convergence for $u$ and $\sigma$ are observed, namely, 1st order for L2 norm of u, L2 norm of $\sigma$ and H(div) norm of $\sigma$. The 2nd order convergent rates between two discrete functions $\|u_I - u_h\|$ and $\|\sigma_I - \sigma_h\|$ are known as superconvergence.
Triangular preconditioned GMRES (the default solver) and Uzawa preconditioned CG converges uniformly in all cases. Traingular preconditioner is two times faster than PCG although GMRES is used.