Implementation of Polygonal Mesh Refinement in MATLAB

The virtual element method (VEM), first introduced in [1], is a generalization of the standard finite element method on general polytopal meshes. In the past few years, people have witnessed rapid progress of virtual element methods (VEMs) for numerically solving various partial differential equations, see for examples [2, 3, 4, 5]. One can also refer to [6] for a transparent MATLAB implementation of the conforming virtual element method for Poisson equation. Due to the large flexibility of the meshes, researchers gradually turn their attention to the posterior error analysis of the VEMs and have made some progress in recent years (cf. [7, 8, 9, 10, 11]).

As we all know, standard adaptive algorithms based on the local mesh refinement can be written as the following loops

Given an initial polygonal subdivision $\mathcal{T}_{0}$, to get $\mathcal{T}_{k+1}$ from $\mathcal{T}_{k}$ we first solve the VEM problem under consideration to get the numerical solution $u_{k}$ on $\mathcal{T}_{k}$. The error is then estimated by using $u_{k}$, $\mathcal{T}_{k}$ and the a posteriori error bound. And the local error bound is used to mark a subset of elements in $\mathcal{T}_{k}$ for refinement. The marked polygons and possible more neighboring elements are refined in such a way that the subdivision meets certain conditions, such as shape regularity. In the implementation, it is usually time-consuming to write a mesh refinement routine since we need to carefully design the rule for dividing the marked elements to get a refined mesh of high quality.

In this paper, we are intended to present an efficient MATLAB implementation of the mesh refinement for polygonal meshes. To the best of our knowledge, this is the first publicly available implementation of the polygonal mesh refinement algorithms. We divide elements by connecting the midpoint of each edge to its barycenter, which may be the most natural partition frequently used in VEM papers, referred to as 4-node subdivision in this context. To remove small edges, some additional neighboring polygons of the marked elements are included in the refinement set by requiring the one-hanging-node rule: limit the mesh to have just one hanging node per edge. We discuss the implementation step by step and give an application in the posteriori error analysis for Poisson equation in the last section. The current implementation or the 4-node subdivision requires that the barycenter is an internal point of each element.

For polygonal meshes, a natural mesh refinement strategy may be the 4-node subdivision as illustrated in Fig. 1 (a) given in [8], where a target element is divided by connecting the midpoint of each edge to the barycenter. For the element with hanging node $Q$, it is better to add a new edge from $Q$ to the barycenter instead of bisecting two edges $PQ$ and $QR$, otherwise degenerate triangles may appear. We refer to this treatment as the admissible bisection. As shown in Fig. 1 (b), the triangle $\Delta 123$ is generated by bisecting $PQ$ and $QR$ and the new triangle $\Delta 567$ is obtained from $\Delta 123$ in a similar manner. According to the properties of triangles, the barycenter $z_{7}$ lies on the median line $e_{34}$ and $|e_{47}|:|e_{73}|=1:2$, which implies that $\angle 7>\angle 3$ and hence a degenerate triangle will appear if we continue with the “two-edge bisection”. On the other hand, the admissible bisection leads to a simpler polygonal mesh.

A mesh refinement function has the following form

[node,elem] = PolyMeshRefine(node,elem,idElemMarked).

Given an initial mesh shown in Fig. 2 (a), we find that the short edge problem reduces to the refinement of the two adjacent polygons ① and ⑤ as a typical case.

• •

  Suppose that ① is a marked element: idElemMarked = 1. If we only refine the marked element as shown in Fig. 2 (b) and proceed to divide the small cells in the lower right corner, then small edges appear in ⑤ due to the frequently added hanging nodes. This phenomenon also happens in the two-edge bisection illustrated in Fig. 1 (b), which is resolved by using the admissible bisection.

• •

  To avoid the occurrence of short edges, we further refine the adjacent polygon ⑤, yielding a new mesh given in Fig. 3. That is, for two elements $K_{1}$ and $K_{2}$ sharing an edge $e$, we refine both $K_{1}$ and $K_{2}$ if $K_{1}$ is in the refinement set and one endpoint of $e$ is the hanging node of $K_{2}$.

• •

  For elements in Fig. 2 (a), ⑤ and ⑦ can be viewed as ① and ⑤, respectively. For the same reason, we also need to partition the polygon ⑦.

• •

  Repeating the above procedure, one can collect all the new elements for refinement and hence avoid producing small edges since the above treatment ensures that each edge in the resulting mesh contains at most one hanging node (a midpoint of the collinear edge).

In the following, an edge is called nontrivial if one of its endpoints is hanging node.

The mesh refinement routine PolyMeshRefine.m is formed by collecting all the code fragments in the last section, which is also available from GitHub (https://github.com/Terenceyuyue/mVEM) as part of the mVEM package containing efficient and easy-following codes for various VEMs published in the literature. The refinement function has been tested for many initial meshes generated by PolyMesher (cf. [13]), a polygonal mesh generator based on the Centroidal Voronoi Tessellations (CVTs), which confirms the effectiveness and correctness of our code. We omit the details and in turn consider the application in the posteriori error estimates of the VEM for Poisson equation.

Consider the Poisson equation with Dirichlet boundary condition on the unit square. The exact solution is given by

The error estimator is taken from [9]. We employ the VEM in the lowest order case and use the Dörfler marking strategy with parameter $\theta=0.4$ to select the subset of elements for refinement.

The initial mesh and the final adapted meshes after 20 and 30 refinement steps are presented in Fig. 4 (a-c), respectively. The detail of the last mesh is shown in Fig. 4 (d). Clearly, no small edges are observed. We also plot the adaptive approximation in Fig. 5, which almost coincides with the exact solution. The full code is still available from mVEM package. The subroutine PoissonVEM_indicator.m is used to compute the local indicators and the test script is main_Poisson_avem.m.