On optimal simplicial 3D meshes for minimizing the Hessian-based errors

Abstract

In this paper we derive a multi-dimensional mesh adaptation method which produces optimal meshes for quadratic functions, positive semi-definite. The method generates anisotropic adaptive meshes as quasi-uniform ones in some metric space, with the tensor metric being computed based on interpolation error estimates. It does not depend, a priori, on the PDEs at hand in contrast to residual methods. The estimated error is then used to steer local modifications of the mesh in order to reach a prescribed level of error in $L^p$-norm or a prescribed number of elements. The $L^p$-norm of the estimated error is then minimized in order to get an optimal mesh. Numerical examples in 2D and 3D for analytic challenging problems and an application to a Computational Fluid Dynamics problem are presented and discussed in order to show how the proposed method recovers optimal convergence rates as well as to demonstrate its computational performance.

Introduction

The success of mesh adaptation techniques resides in their ability to enhance the accuracy of the solution and at the same time reducing the computation cost. Even though, the accuracy can be highly bettered by identifying the regions that need more refinement such those near the singularities. It is crucial to achieve a good equilibrium between the coarsened and refined regions, such that the global accuracy will be optimal. In Computational Fluid Dynamics (CFD), a priori error estimates, as provided by the standard finite element error analysis, are usually insufficient to ensure reliable estimates of the numerical solution, because they only give information on the asymptotic error behavior and need regularity conditions that are unsatisfied in the presence of singularities [18]. Those findings suggest the need for an a posteriori error estimator which can be derived from the numerical solution and the given data of the problem.

Indeed, we are motivated to build robust and efficient numerical methods for solving CFD problems with steep interfaces and boundary layers. For such problems, the solution changes much faster in some directions than it does in others. However, it is always unknown a priori where these layers are located or in which directions the solution changes most quickly. A numerical method should track down the layers automatically and be efficient and robust in resolving the layers without oscillations in the numerical approximations. For this purpose, we use the stabilized finite element methods which are stable and have a good higher order accuracy in regions where the solution is smooth. However, finding the optimal stabilization terms to completely reduce numerical oscillations is still an open problem. Therefore, in practice, it might be impossible to achieve optimal convergence rates if the solution has arduous layers. A recent research orientation focuses on the use of mesh adaptation to the stabilized numerical solutions in order to enhance the convergence rate.

This topic is not actually new in the computational area, even in the finite element literature. The search of optimal meshes dates back to the early 1970s [29]. But modern interest in this subject began in the late 1970s, mainly thanks to important contributions by Babuska and Rheinboldt [5], [6]. In the last two decades, the subject has become increasingly important in finite element practice and estimators based on averaging techniques have become extremely popular; in particular the one proposed by Zienkiewicz and Zhu [40], [41] and then several of declination and improvement proposed thereafter [37], [38], [7], [8], [39]. Vallet and coworkers in [35] propose a numerical comparison of different recovery techniques. This great excitement for these techniques has grown with the development of anisotropic mesh adaptation tools able to handle error estimate information as a mesh size and direction prescriptions. The anisotropy is materialized in the mesh by stretching elements along the main solution directions. The stretching is integrated through a Riemannian metric space as reported in [30], [31], [25], [26], [9], [10], [34], [15], [21], [19], [23] and references therein. An adapted mesh is generated with respect to this metric where the aim is to generate a mesh such that all edges have a length equal (or close) to unity in the prescribed metric and such that all elements are almost regular. The volume is adapted by local mesh modification of the previous mesh (the mesh is not regenerated) using mesh operations: node insertion, edge and face swap, collapse and node displacement. This approach can lead to elements with large angles that are not suitable for finite element computations as reported in the general standard error analysis for which some regularity assumption on the mesh and on the exact solution should be satisfied [18]. However, if the mesh is adapted to the solution, it is possible to circumvent this condition [33].

Several recent results [27], [1], [22], [24], [12], [28], [20] and papers therein have brought renewed focus on metric-based adaptation where the underling metric derived from a recovered Hessian. As the former works, we focus on an anisotropic mesh adaptation process, driven by a directional error estimator based on the recovery of the Hessian of the finite element solution. In particular, we pursue the work introduced in [2] for 2D problems. In this work, we extend and generalize the analysis to higher N-dimensional problems $(N\ge 2)$ and derive an optimal multi-dimensional metric. The purpose is to achieve an optimal 3D mesh minimizing the directional error estimator for a given number of mesh elements. It allows, as will be shown along this paper, to refine/coarsen the mesh, stretch and orient the elements in such a way that, along the adaptation process, the mesh becomes aligned with the fine scales of the solution whose locations are unknown a priori. As a result of this, highly accurate solutions are obtained with a much lower number of elements.

This paper is organized as follows: Section 2.1 introduces the notions and notations that will be used in the following. Sections 2.2 and 2.3 focus on the main contribution of this paper, the extension of the anisotropic estimator [2] to multi-dimensional unstructured meshes. Section 3 describes the main local remeshing operations used in the mesh adaptation procedure. In Section 4 we show through numerical results, the advantage of the combined mesh adaptation algorithm and anisotropic error estimator to handle accurately stiff 2D and 3D problems as well as a multiphase CFD problem.