This is a survey of the theory of adaptive finite element methods (AFEMs), which are fundamental to modern computational science and engineering but whose mathematical assessment is a formidable challenge. We present a self-contained and up-to-date discussion of AFEMs for linear second-order elliptic PDEs and dimension d > 1, with emphasis on foundational issues. After a brief review of functional analysis and basic finite element theory, including piecewise polynomial approximation in graded meshes, we present the core material for coercive problems. We start with a novel a posteriori error analysis applicable to rough data, which delivers estimators fully equivalent to the solution error. They are used in the design and study of three AFEMs depending on the structure of data. We prove linear convergence of these algorithms and rate-optimality provided the solution and data belong to suitable approximation classes. We also address the relation between approximation and regularity classes. We finally extend this theory to discontinuous Galerkin methods as prototypes of non-conforming AFEMs, and beyond coercive problems to inf-sup stable AFEMs.

In this paper, we study the role of mesh quality on the accuracy of linear finite element approximation. We derive a more detailed error estimate, which shows explicitly how the shape and size of elements, and symmetry structure of mesh effect on the error of numerical approximation. Two computable parameters Ge and Gv are given to depict the cell geometry property and symmetry structure of the mesh. In compare with the standard a priori error estimates, which only yield information on the asymptotic error behaviour in a global sense, our proposed error estimate considers the effect of local element geometry properties, and is thus more accurate. Under certain conditions, the traditional error estimates and supercovergence results can be derived from the proposed error estimate. Moreover, the estimators Ge and Gv are computable and thus can be used for predicting the variation of errors. Numerical tests are presented to illustrate the performance of the proposed parameters Ge and Gv.

An initial-boundary value problem for a time-fractional diffusion equation is discretized in space, using continuous piecewise-linear finite elements on a domain with a re-entrant corner. Known error bounds for the case of a convex domain break down, because the associated Poisson equation is no longer $H^{2}$-regular. In particular, the method is no longer second-order accurate if quasi-uniform triangulations are used. We prove that a suitable local mesh refinement about the re-entrant corner restores second-order convergence. In this way, we generalize known results for the classical heat equation.

A moving mesh method is proposed for solving reaction-diffusion equations. The finite element method is used to solving the partial different equation system, and an efficient numerical scheme is applied to implement mesh moving. In the practical calculations, the moving mesh step and the problem equation solver are performed alternatively. Several numerical examples are presented, including the Gray-Scott, the Activator-Inhibitor and a case with a growing domain. It is illustrated numerically that the moving mesh methods costs much lower, compared with the numerical schemes on a fixed mesh. Even in the case of complex pattern dynamics described by the reaction-diffusion systems, the adapted meshes can capture the details successfully.

Mesh generation is a bottleneck for finite element simulations of biomolecules. A robust and efficient approach, based on the immersed boundary method proposed in [8], has been developed and implemented to generate large-scale mesh body-fitted to molecular shape for general parallel finite element simulations. The molecular Gaussian surface is adopted to represent the molecular surface, and is finally approximated by piecewise planes via the tool phgSurfaceCut in PHG [43], which is improved and can reliably handle complicated molecular surfaces, through mesh refinement steps. A coarse background mesh is imported first and then is distributed into each process using a mesh partitioning algorithm such as space filling curve [5] or METIS [22]. A bisection method is used for the mesh refinements according to the molecular PDB or PQR file which describes the biomolecular region. After mesh refinements, the mesh is optionally repartitioned and redistributed for load balancing. For finite element simulations, the modification of region mark and boundary types is done in parallel. Our parallel mesh generation method has been successfully applied to a sphere cavity model, a DNA fragment, a gramicidin A channel and a huge Dengue virus system. The results of numerical experiments show good parallel efficiency. Computations of electrostatic potential and solvation energy also validate the method. Moreover, the meshing process and adaptive finite element computation can be integrated as one PHG project to avoid the mesh importing and exporting costs, and improve the convenience of application as well.

A uniaxial perfectly matched layer (PML) method is proposed for solving the scattering problem with multiple cavities. By virtue of the integral representation of the scattering field, we decompose the problem into a system of single-cavity scattering problems which are coupled with Dirichlet-to-Neumann maps. A PML is introduced to truncate the exterior domain of each cavity such that the computational domain does not intersect those for other cavities. Based on the a posteriori error estimates, an adaptive finite element algorithm is proposed to solve the coupled system. The novelty of the proposed method is that its computational complexity is comparable to that for solving uncoupled single-cavity problems. Numerical experiments are presented to demonstrate the efficiency of the adaptive PML method.

In this work, we develop an adaptive algorithm for solving elliptic optimal control problems with simultaneously appearing state and control constraints. The algorithm combines a Moreau-Yosida technique for handling state constraints with a semi-smooth Newton method for solving the optimality systems of the regularized sub-problems. The state and adjoint variables are discretized using continuous piecewise linear finite elements while a variational discretization concept is applied for the control. To perform the adaptive mesh refinements cycle we derive local error estimators which extend the goal-oriented error approach to our setting. The performance of the overall adaptive solver is assessed by numerical examples.

Dislocations are line defects in crystalline materials. The Peierls-Nabarro models are hybrid models that incorporate atomic structure of dislocation core into continuum framework. In this paper, we present a numerical method for a generalized Peierls-Nabarro model for curved dislocations, based on the fast multipole method and the iterative grid redistribution. The fast multipole method enables the calculation of the long-range elastic interaction within operations that scale linearly with the total number of grid points. The iterative grid redistribution places more mesh nodes in the regions around the dislocations than in the rest of the domain, thus increases the accuracy and efficiency. This numerical scheme improves the available numerical methods in the literature in which the long-range elastic interactions are calculated directly from summations in the physical domains; and is more flexible to handle problems with general boundary conditions compared with the previous FFT based method which applies only under periodic boundary conditions. Numerical examples using this method on the core structures of dislocations in Al and Cu and in epitaxial thin films are presented.

A weak Galerkin discretization of the boundary value problem of a general anisotropic diffusion problem is studied for preservation of the maximum principle. It is shown that the direct application of the M-matrix theory to the stiffness matrix of the weak Galerkin discretization leads to a strong mesh condition requiring all of the mesh dihedral angles to be strictly acute (a constant-order away from 90 degrees). To avoid this difficulty, a reduced system is considered and shown to satisfy the discrete maximum principle under weaker mesh conditions. The discrete maximum principle is then established for the full weak Galerkin approximation using the relations between the degrees of freedom located on elements and edges. Sufficient mesh conditions for both piecewise constant and general anisotropic diffusion matrices are obtained. These conditions provide a guideline for practical mesh generation for preservation of the maximum principle. Numerical examples are presented.

A new numerical method based on locally modified Cartesian meshes is proposed for solving a coupled system of a fluid flow and a porous media flow. The fluid flow is modeled by the Stokes equations while the porous media flow is modeled by Darcy's law. The method is based on a Robin-Robin domain decomposition method with a Cartesian mesh with local modifications near the interface. Some computational examples are presented and discussed.

In this paper we present a first supercloseness analysis for higher-order Galerkin FEM applied to a singularly perturbed convection-diffusion problem. Using a solution decomposition and a special representation of our finite element space, we are able to prove a supercloseness property of p + 1/4 in the energy norm where the polynomial order p ≥ 3 is odd.

We present quantitative studies of transfer operators between finite element spaces associated with unrelated meshes. Several local approximations of the global L2-orthogonal projection are reviewed and evaluated computationally. The numerical studies in 3D provide the first estimates of the quantitative differences between a range of transfer operators between non-nested finite element spaces. We consider the standard finite element interpolation, Clément’s quasi-interpolation with different local polynomial degrees, the global L2-orthogonal projection, a local L2-quasi-projection via a discrete inner product, and a pseudo-L2-projection defined by a Petrov-Galerkin variational equation with a discontinuous test space. Understanding their qualitative and quantitative behaviors in this computational way is interesting per se; it could also be relevant in the context of discretization and solution techniques which make use of different non-nested meshes. It turns out that the pseudo-L2-projection approximates the actual L2-orthogonal projection best. The obtained results seem to be largely independent of the underlying computational domain; this is demonstrated by four examples (ball, cylinder, half torus and Stanford Bunny).

This paper studies convergence analysis of an adaptive finite element algorithm for numerical estimation of some unknown distributed flux in a stationary heat conduction system, namely recovering the unknown Neumann data on interior inaccessible boundary using Dirichlet measurement data on outer accessible boundary. Besides global upper and lower bounds established in [23], a posteriori local upper bounds and quasi-orthogonality results concerning the discretization errors of the state and adjoint variables are derived. Convergence and quasi-optimality of the proposed adaptive algorithm are rigorously proved. Numerical results are presented to illustrate the quasi-optimality of the proposed adaptive method.

A novel method for boundary constrained tetrahedral mesh generation is proposed based on Advancing Front Technique (AFT) and conforming Delaunay triangulation. Given a triangulated surface mesh, AFT is firstly applied to mesh several layers of elements adjacent to the boundary. The rest of the domain is then meshed by the conforming Delaunay triangulation. The non-conformal interface between two parts of meshes are adjusted. Mesh refinement and mesh optimization are then preformed to obtain a more reasonable-sized mesh with better quality. Robustness and quality of the proposed method is shown. Convergence proof of each stage as well as the whole algorithm is provided. Various numerical examples are included as well as the quality of the meshes.

A cell conservative flux recovery technique is developed here for vertex-centered finite volume methods of second order elliptic equations. It is based on solving a local Neumann problem on each control volume using mixed finite element methods. The recovered flux is used to construct a constant free a posteriori error estimator which is proven to be reliable and efficient. Some numerical tests are presented to confirm the theoretical results. Our method works for general order finite volume methods and the recovery-based and residual-based a posteriori error estimators is the first result on a posteriori error estimators for high order finite volume methods.

Goal of this paper is to suitably combine a model with an anisotropic mesh adaptation for the numerical simulation of nonlinear advection-diffusion-reaction systems and incompressible flows in ecological and environmental applications. Using the reduced-basis method terminology, the proposed approach leads to a noticeable computational saving of the online phase with respect to the resolution of the reference model on nonadapted grids. The search of a suitable adapted model/mesh pair is to be meant, instead, in an offline fashion.

A novel adaptive approach to compute the eigenenergies and eigenfunctions of the two-particle (electron-hole) Schrödinger equation including Coulomb attraction is presented. As an example, we analyze the energetically lowest exciton state of a thin one-dimensional semiconductor quantum wire in the presence of disorder which arises from the non-smooth interface between the wire and surrounding material. The eigenvalues of the corresponding Schrödinger equation, i.e., the one-dimensional exciton Wannier equation with disorder, correspond to the energies of excitons in the quantum wire. The wavefunctions, in turn, provide information on the optical properties of the wire.

We reformulate the problem of two interacting particles that both can move in one dimension as a stationary eigenvalue problem with two spacial dimensions in an appropriate weak form whose bilinear form is arranged to be symmetric, continuous, and coercive. The disorder of the wire is modelled by adding a potential in the Hamiltonian which is generated by normally distributed random numbers. The numerical solution of this problem is based on adaptive wavelets. Our scheme allows for a convergence proof of the resulting scheme together with complexity estimates. Numerical examples demonstrate the behavior of the smallest eigenvalue, the ground state energies of the exciton, together with the eigenstates depending on the strength and spatial correlation of disorder.

The Bay of Bengal is surrounded by coastline except to the south, where there is open sea. The coastline bends most sharply along the coast of Bangladesh, and there are many small and large islands in the off shore region there. In order to incorporate the island boundaries and the curved coastline properly, in any numerical scheme it is often necessary to consider a very fine grid resolution along the coastal belts whereas this is unnecessary away from the coasts. However, a very fine resolution involves more memory and more CPU time in the numerical solution process, and invites numerical instability. On the other hand, boundary-fitted curvilinear grids in hydrodynamic models for coastal seas, bays and estuaries not only fit to the coastline but also render the finite difference schemes simpler and more accurate. In this article, the boundary-fitted curvilinear grids for the model represent the complete boundary of the area considered by four curves defined by four functions, and the four boundaries of two of the larger islands are then represented approximately by two general functions. An appropriate independent coordinate transformation maps the curvilinear physical area to a square domain, and each island boundary is transformed to a rectangle within this square domain. The vertically integrated shallow water equations are transformed to the new space domain, and solved by a regular explicit finite difference scheme. The model is applied to compute the water levels due to astronomical tides, and also the water levels due to surges associated with tropical storms that hit the coast of Bangladesh.

We present a solver for the Poisson-Boltzmann equation and demonstrate its applicability for biomolecular electrostatics computation. The solver uses a level set framework to represent sharp, complex interfaces in a simple and robust manner. It also uses non-graded, adaptive octree grids which, in comparison to uniform grids, drastically decrease memory usage and runtime without sacrificing accuracy. The basic solver was introduced in earlier works [16,27], and here is extended to address biomolecular systems. First, a novel approach of calculating the solvent excluded and the solvent accessible surfaces is explained; this allows to accurately represent the location of the molecule’s surface. Next, a hybrid finite difference/finite volume approach is presented for discretizing the nonlinear Poisson-Boltzmann equation and enforcing the jump boundary conditions at the interface. Since the interface is implicitly represented by a level set function, imposing the jump boundary conditions is straightforward and efficient.