Bernstein-Bézier Galerkin-Characteristics Finite Element Method for Convection-Diffusion Problems

A class of Bernstein-Bézier basis based high-order finite element methods is developed for the Galerkin-characteristics solution of convection-diffusion problems. The Galerkin-characteristics formulation is derived using a semi-Lagrangian discretization of the total derivative in the considered problems. The spatial discretization is performed using the finite element method on unstructured meshes. The Lagrangian interpretation in this approach greatly reduces the time truncation errors in the Eulerian methods. To achieve high-order accuracy in the Galerkin-characteristics solver, the semi-Lagrangian method requires high-order interpolating procedures. In the present work, this step is carried out using the Bernstein-Bézier basis functions to evaluate the solution at the departure points. Triangular Bernstein-Bézier patches are constructed in a simple and inherent manner over finite elements along the characteristics. An efficient preconditioned conjugate gradient solver is used for the linear systems of algebraic equations. Several numerical examples including advection-diffusion equations with known analytical solutions and the viscous Burgers problem are considered to illustrate the accuracy, robustness and performance of the proposed approach. The computed results support our expectations for a stable and highly accurate Bernstein-Bézier Galerkin-characteristics finite element method for convection-diffusion problems.


Introduction
Convection-diffusion equations have been widely used to model many applications in engineering involving advection of scalar quantities such as density, temperature or concentration among others. For instance, convection-diffusion problems have been proposed to study water transport in soils [25], heat transfer of nanofluids [36], and the transport in ferrofluids under rotating magnetic fields [10]. On the other hand, various numerical methods have been developed in recent years to solve the convection-diffusion equations. Most of these computational techniques can be classified into three main categories: (i) Eulerian methods, (ii) Lagrangian methods and (iii) semi-Lagrangian methods. In the framework of finite elements, the most popular Eulerian methods are the streamline upwind Petrov-Galerkin methods [9,11], Galerkin/least-squares methods [9,26] and Taylor-Galerkin methods [12,13]. However, it is well known that these Eulerian methods do not perform very satisfactory in the case of convection-dominated problems unless small time steps and highly refined grids are used in the simulations. In the case of pure convection problems as those considered in this study, these requirements are practically not feasible and may limit the performance of these Eulerian methods. The Lagrangian techniques on the other hand, are theoretically well suited for the numerical solution of advection problems due to the possibility of using large time steps in the simulations. In practice, the computational mesh for the Lagrangian methods moves along the fluid particle trajectory which may yield to mesh distortion after few time steps in the computations. Thus, because of this drawback, the Lagrangian methods are not recommended for the numerical solution of complex convection-diffusion problems. In the semi-Lagrangian methods known also in the framework of finite elements by Galerkincharacteristics, the computational mesh is taken to be fixed to overcome the drawback of the Lagrangian methods while keeping the advantage of the Lagrangian tracking algorithm along the characteristic curves. The main advantage of the Galerkin-characteristics method lies on the fact that the Courant-Friedrichs-Lewy (CFL) condition is highly relaxed compared to its Eulerian counterparts, see for example [15-18, 32, 39, 42]. In addition, the Lagrangian treatment in the Galerkin-characteristics algorithms greatly reduces the time truncation errors in the Eulerian methods, see [7,8,15,38,41] among others. Thus, the Galerkin-characteristics finite element method has the potential to be more suitable than the Eulerian and Lagrangian methods for convection-diffusion problems on unstructured meshes. Numerical assessment of the conventional Galerkin-characteristics method has been carried out in [16] for convection problems and comparisons to well-established Eulerian methods can also been found in this reference.
In general, most of Galerkin-characteristics finite element methods are second-order accurate in space and time but the accuracy of this class of numerical methods depends on the order of the interpolation polynomials used to compute the solution in the convection stage and on the time integration procedure for the diffusion operator. For example, to achieve a second-order accuracy in the Galerkin-characteristics finite element method, the interpolation polynomials have to be at least second-order accurate and the time integration must be at least semi-implicit for the diffusion terms. In addition, it has been observed that the error in the conventional Galerkin-characteristics finite element method for convection-diffusion problems decreases as the time step increases at certain range of parameters, see for instance [17,18,20]. High-order accurate numerical methods for convection-dominated problems have the potential to reduce the computational effort required for a given order of solution accuracy. The state of the art in this field is more advanced for the Eulerian methods than for the semi-Lagrangian methods. For example, high-order discretization techniques such as those relying on spectral or hp-finite element method have been shown to achieve fast convergence with low numerical diffusion and dispersion errors for advection-diffusion problems [22,28]. The p-finite element method and hp-finite element method where introduced by Babuška, Szabó and their coworkers in the mid-1970s. It was shown that the hp-finite element method delivers exponential convergence for elliptic problems with piecewise analytic data, see the survey [6] and further references are therein. A similar performance was also proven for boundary-layer and singularly perturbed problems [33,37]. A study reported in [5] reveals that the choice of high-order shape functions is critical to the stability and efficiency of the finite element procedure. Particularly, high-order finite elements based on Lobatto shape functions have proven to possess better conditioning compared to other types of high-order shape functions widely used in the literature [43]. Assessment of different high-order shape functions in [35], including Bernstein, Lobatto and Lagrange Gauss-Lobatto polynomials for interior acoustic problems, has shown the advantage of high-order polynomials in reducing the pollution errors and the good performance of Bernstein polynomials when combined with the Krylov subspace solvers. In a closely related study [19], Bernstein shape functions have been demonstrated to yield comparable, and even better performance in terms of accuracy and memory requirements compared to the well-established partition of unity finite method. The Bernstein polynomials are well known in the field of computer aided geometric design and computer graphics. However, their applications in the finite element community have until now not been widely adopted. Although hierarchical basis functions are often chosen in the design of high-order finite elements for their suitability in p-adaptivity, recently attention has been paid to the favorable properties of the Bernstein polynomials [1,24,29]. Especially, it has been shown that Bernstein-Bézier finite elements over simplicial domains, hexahedra and pyramids yield optimal complexity for the standard finite element spaces. In a more recent work [3], the Bernstein basis functions combined to an additive Schwarz preconditioner was successfully implemented for challenging applications including boundary layers, non-linear reaction-diffusion problems and wave propagation of solitons.
The main focus of the present study is the development of a class of high-order Galerkin-characteristics finite element methods to numerically solve convection-dominated problems. This goal is achieved by the implementation of Bernstein-Bézier finite elements for the Galerkin-characteristics method. It should be stressed that combining the Galerkincharacteristics method with the Bernstein-Bézier finite elements, to the best of our knowledge, is reported for the first time. In the context of Galerkin-characteristics finite element methods, Bernstein polynomials are used as shape functions associated with elements of the computational mesh to calculate the departure points and update the global solutions. The positivity of these local basis functions and the variation diminishing properties make them a very attractive alternative to the standard Lagrange polynomials. To increase the efficiency of the proposed method, we also implement an efficient preconditioned iterative solver using the mass matrix as a preconditioner as proposed in [3]. The key idea behind this considered algorithm lies on the static condensation of the cell-based Degrees of Freedoms (DoFs) on each element and it can be viewed as an extension of the Additive Schwarz method. This results in a p-independent uniform bound on the growth of the associated condition number. It should also be noted that in case of the pure convection problem, the resulting linear system in the proposed semi-Lagrangian method consists of the mass matrix which make the proposed preconditioned conjugate gradient solver more convenient than its convection-diffusion counterpart. Furthermore, for the convection-dominated problems the growth of the condition number is relaxed by the viscosity coefficient and therefore, the performance of the considered linear solver is not affected. Numerical results presented in this work demonstrate that an interesting feature of the Bernstein-Bézier finite elements is to allow large time steps and coarse meshes in the simulations without deteriorating the high-order accuracy of the computed solutions.
This paper is organized as follows. Introduction of the Bernstein-Bézier finite element discretization is presented in Sect. 2. Section 3 is devoted to the formulation of the Bernstein-Bézier Galerkin-characteristics finite element method for pure convection problems. This section includes the calculation of characteristic curves and the implementation of the computational algorithm. The extension of the method for solving convection-diffusion equations is discussed in Sect. 4. In Sect. 5, we examine the numerical performance of the proposed method using several test examples of convection-diffusion problems including the viscous Burger equation. The proposed Bernstein-Bézier Galerkin-characteristics finite element method is demonstrated to enjoy the expected high-order accuracy as well as efficiency. Concluding remarks are summarized in Sect. 6.

Bernstein-Bézier Finite Element Approximation
In order to formulate our Galerkin-characteristics method, the computational domain assumed for simplicity to be a two dimensional bounded and polygonal domain, is first partitioned into non-overlapping triangular finite elements. Let T be the reference element defined by where q 1 = (0, 0), q 2 = (1, 0) and q 3 = (0, 1) are vertices of the reference element as shown in Fig. 1. The barycentric coordinates λ i (i = 1, 2, 3) relative to the reference element T are defined by The Bernstein-Bézier basis for the space P p ( T ) of polynomials of total degree at most p consists of the following shape functions: • Vertex-based shape functions defined by • Cell-based shape functions defined by Using these notations, the Bernstein-Bézier shape functions can also be formulated in a simple compact form as where p α = p! α! and |α| = p. It should also be stressed that one of the most important proprieties of the Bernstein polynomials lies on the fact that their product yields a scaled Bernstein polynomial as where |α| = p and |β| = q. Furthermore the integral of a Bernstein polynomial over the reference triangle element has a simple form On the other hand the gradient of Bernstein polynomials can be computed as follows: where e k is a multi-index with its kth entry is a unity and its remaining entries are zero, and B p−1 α−e k = 0 if α −ẽ k has a negative component. It should be noted that the Bernstein polynomials are non-negative and form a partition of unity on the element T . Moreover, these polynomials have some attractive features such as variation diminishing and monotonicity preserving properties, see for instance [21,23] and further references can be found therein.
Following the same procedure carried out in the classical finite element methods, each mesh element T in the triangulation T h is mapped to the reference element T in which all computations are performed. Here, we use the following reference affine map where q 1 , q 2 and q 3 are the vertices of the physical element T . The conforming finite element space for the solution that we consider is defined as Here, T is a given mesh element, ξ j is a Stroud quadrature point used in the reference element T and mapped onto x h j in the element T , T * j is the host element where the departure point X n h j belongs, and T is the affine one-to-one mapping For a polynomial degree p ≥ 1, the number of Degrees of Freedom (DoF) per element is given by Note that, in contrast to the Lagrange finite elements where the degrees of freedom refer to nodal point evaluations, a global orientation of edges is required to enable matching edge modes of a similar shape, thus ensuring C 0 conformity [19,28,43]. Finite element approximations with varying polynomial order can also be used by assigning to each vertex, edge and cell in the mesh an arbitrary polynomial degree, see [1,19] among others.

Bernstein-Bézier Finite Element Galerkin-Characteristics Method
To describe the formulation of the proposed Bernstein-Bézier finite element Galerkincharacteristics method, we consider the following two-dimensional Cauchy problem for the convection equation where x = (x, y) is the space variable, ∇c = ( ∂c ∂ x , ∂c ∂ y ) is the gradient vector, is a spatial bounded domain in R 2 with boundary ∂ , and [0, T ] is a time interval. Here, c(x, t) denotes the concentration of some species, v(x, t) = (u(x, t), v(x, t)) is the velocity field assumed to depend on the solution c as well, and c 0 (x) is a given initial function. We assume that appropriate boundary conditions are given in such a way the problem is well defined and has a unique solution. Note that Dc Dt in (10) measures the rate of change of the concentration c following the trajectories of the flow particles. The main idea behind the Galerkin-characteristics method is to impose a regular grid at the new time level and to backtrack the flow trajectories to the previous time level. At the old time level, the quantities that are needed are evaluated by interpolation from their known values on a regular grid.
Next, we divide the time interval [0, T ] into N subintervals [t n , t n+1 ] with length t = t n+1 −t n for n = 0, 1, . . . , N . We use the notation w n to denote the value of a generic function w at time t n . Hence, the characteristic curves associated with the convection problem (10) are the solution of the backward differential equation is the departure point at time τ of a particle that will arrive at x = (x, y) at time t n+1 . Note that the Galerkin-characteristics methods do not follow the flow particles forward in time, as the Lagrangian schemes do, instead they trace backwards the position at time t n of particles that will reach the points of a fixed mesh at time t n+1 , see Fig. 1 for an illustration. By so doing, the Galerkin-characteristics methods avoid the grid distortion difficulties that the conventional Lagrangian methods have. The solution of the differential eq. (11) can be expressed as Hence, integrating the convection eq. (10) along the characteristic curves yields Note that the Bernstein polynomials are not interpolatory by reconstruction i.e., the degrees of freedom associated to the shape functions do not directly lie on the solution evaluated at control points. Hence, the solution at the next time level should be obtained by a weak formulation. Thus, multiplying both sides of equation (13) by a test function w ∈ H 1 ( ) and integrating over , it leads to the following weak form Here, the finite element solution c n+1 h is sought element-wise at each time step as with x = T (ξ ) and T ∈ T h . The approximation of the weak form (14) using the conforming finite element space V h yields the following linear system of algebraic equations with M is an n ndof × n ndof -valued sparse symmetric matrix, b is an n ndof -valued right-hand side column vector and C is n ndof -valued column vector of the unknowns, where n ndof is the total number of degrees of freedom. The entries of the mass element matrix can be written as where |α| = |β| = p and J T = d T dξ is the Jacobian matrix. Since the geometry is interpolated using the affine map T , analytical integration rules as those proposed in [1,29] can be used for the evaluation of the element mass matrix. Thus, using the proprieties (4) and (5), the integral in (17) can be evaluated as Note that the crucial step in this approach is the evaluation of the right-hand entries in (16) given by is evident that, if the above integrals are evaluated exactly then it is easy to show that the Galerkin-characteristics method is unconditionally stable in the L 2 -norm, see for instance [15,40]. In a general framework, this cannot be done and one has to approximate the integrals by numerical integration. It has shown in [1] that the Duffy transformation enables a tensorial reconstruction of the Bernstein-Bézier basis on simplices. Thus, the well-established sum factorization is used to efficiently evaluate and integrate these polynomials based on the Stroud conical quadrature. This transformation maps the unit quadrilateral with coordinates t = (t 1 , t 2 ) ∈ S = [0, 1] 2 to the reference triangle T and it can be defined by Let us denote by where |α| = p. Therefore, the Bernstein polynomial form (15) becomes where the weight function , the solution c n h can be evaluated at the point X n h j as Note that if 1 − ξ * j ,1 = 0 the Duffy transformation degenerates and since in this case At the implementation level, the interpolation of c n h at the departure point X n h j is performed in the same manner as in [1], by exploiting property (25). Thus, we first evaluate all the univariate Bernstein polynomials B p α 1 and B p−α 1 α 2 , with α 1 = 0, 1, . . . , p and α 2 = 0, 1, . . . , p − α 1 , at the points t * j ,1 and t * j ,2 , respectively, based on the recursion formula where B k −1 = B k k+1 = 0 and B 0 0 = 1. Note that this procedure requires O( p 2 ) floating point operations. Next, for α 1 = 0, 1, . . . , p, we compute the auxiliary coefficients which yields a cost of O( p 2 ). Hence, we evaluate c n h (X n h j ) as Notice that the computational cost of this algorithm per point is of O( p 2 ) operations for an arbitrary point. It should also be stressed that De Casteljau's algorithm [31] can be used for evaluating the previous coefficients. However, this results in O( p 3 ) operations, unless the given point X n h j belongs to an edge of the triangle for which the computational cost will decrease to O( p 2 ) operations. Once c n h (X n h j ) is computed for all j 1 , j 2 = 1, . . . , q, each entry b n α given by (24) can be evaluated in O(q 2 ) operations. Since the required number q of quadrature points should be chosen such that q = O( p) to ensure a sufficiently accurate numerical integration, the total cost to set up the element right-hand side is O( p 2 ).
In the current study, to compute the departure point X n h j from a given integration point x h j , the discrete analogous of (12) is first reformulated as where the displacement δ h j is evaluated using the following iterative procedure based on the Adams-Bashforth method h j evaluated based on the finite element interpolation on the mesh element where h j is located. This iterative procedure was first proposed in [41] for the finite difference semi-Lagrangian methods and investigated in [15] for the finite element discretizations. In the present work, the iterations (31) are stopped when the following criteria is satisfied for the Euclidean norm · and a given tolerance . In general, the departure points do not coincide with the spatial position of a mesh point in the triangular mesh. Therefore, the method used to compute X n h j should be equipped with a search-locate algorithm to find the host element where such point is located. In our simulations presented in Sect. 5, we have implemented a search-locate algorithm designed in [4] for the semi-Lagrangian methods in unstructured finite element discretizations. In addition, the iterations in (31) were continued until the trajectory changed by less than = 10 −6 .

Implementation for Convection-Diffusion Problems
In this section we consider convection-diffusion problems reformulated using the total derivative as where ν is the diffusion coefficient and f (x, t) the source term. We assume that eq. (33) is equipped with well defined boundary and initial conditions depending on the problem under study. In the current study, to deal with the diffusion part in the eq. (33) we consider a secondorder implicit scheme of Gear type also known in the literature by backward differentiation formula (BDF2). Using the same notations introduced in the previous section, the time semidiscrete form of the eq. (33) reads which can be rearranged in a compact form as with F n = 4 2 t c n • X n − 1 2 t c n−1 • X n−1 + f n+1 . Notice that to advance the solution c n+1 in time in (34), the two solutions c n−1 and c n are required. At time t = 0 only one initial condition is provided and to obtain the second condition we use the implicit Euler scheme. Suppose for simplicity purposes, the convection-diffusion problem (33) is supplied with an homogeneous Dirichlet boundary condition on ∂ . Then, the discrete weak form of (35) reads as: Find c n+1 where V h is the conforming finite element space defined in (8). By virtue of definitions of the finite element operators given above, this reduces to the following linear system of algebraic equations where K is an n ndof × n ndof -valued sparse symmetric matrix and b n is an n ndof -valued right-hand side column vector. It should also noted that since the Bernstein polynomials are only interpolatory at the mesh grid vertices, a numerical procedure is needed to impose nonhomogeneous Dirichlet type boundary condition. In the present work, this is achieved by using the L 2 projection on the Bernstein-Bézier basis of the local boundary data Lagrange interpolate.
The entries K α,β of the element stiffness matrix K are given by and the entries b n α of the right-hand side column vector b n are given by where |α| = |β| = p, ∇ refers to the gradient with respect to the local coordinates ξ . Following the same procedure used in the previous section to compute analytically the entries of the elements mass matrix, one can easily verify that the computation of the entries of the element stiffness matrix can be carried out analytically using eq. (6) as where M p−1 α−e k ,β−e l is deduced from the closed form (18). In addition, the integral in (39) is computed in the same manner as previously based on the Stroud quadrature and property (22). Here, a rule of q = p + 2 quadrature points which is exact if the integrand is a polynomial of degree no higher than 2 p + 3 is adopted.
In summary, the Bernstein-Bézier Galerkin-characteristics finite element method to solve either the convection problem (10) or convection-diffusion problem (33) is carried out in the following steps: Note that, since the cell-based shape functions are internal i.e., they vanish on the element boundaries and are therefore not connected to the neighboring elements, the static condensation can be applied at the elemental level to remove the internal DoFs from the global finite element system during the assembly. Once the matrix and the right-hand side of the statically condensed system are formed, the internal DoFs in the solution can be recovered during the post-processing by solving element-wise local linear problems. This procedure is very efficient in reducing the size and enhancing the condition number of hp-finite element Algorithm The Bernstein-Bézier Galerkin-characteristics finite element method 1: for each time step do 2: for each mesh element T do 3: For each Stroud quadrature point x h j calculate the departure points X n h j using (30)-(31).

4:
Identify the mesh element T * j where the departure point X n h j belongs using the search-locate algorithm proposed in [4].

5:
Evaluate the gridpoint approximations using (25). 6: Compute the element right-hand side using (39) and assemble it in the global column vector b n . 7: end for 8: Compute C n+1 by solving the linear system (16) or (37). 9: Update the solution c n+1 h element-wise according to (22). 10: end for system matrices. Furthermore, the considered method requires solution of uncoupled elliptic problems such that their finite element discretization leads to linear systems of algebraic equations for which, very efficient solvers can be implemented. Therefore, by taking advantage of these properties, we solve the linear systems in (16) or (37) by the preconditioned conjugate gradient solver. This yields an efficient method for solving this class of linear systems of algebraic equations for which the preconditioner is implemented at a cost of O( p 3 ) operations as discussed in [3]. Here, the additive Schwarz preconditioner is used in the same manner as studied in [3]. The key idea consists in using the relationship between the Bernstein and Jacobi polynomials which makes the passage from the Jacobi to the Bernstein basis functions and vice versa easy without inverting any matrix. These techniques have been proposed in [3] to build a preconditionner for the Bernstein polynomials using the well-established preconditionner introduced in [2] for the Jacobi polynomials but with less computational cost. This gives rise to a preconditioned system for which the condition number is bounded independently of the polynomial order p and the mesh size h.

Results and Examples
where c n exact and c n h are respectively, the exact and numerical solutions at gridpoint x h and time t n . We also define the CFL number associated to the problems (10) and (33) as In all our computations carried out in this section, the resulting linear systems of algebraic equations are solved using the preconditioned conjugate gradient solver and stopping criteria set to 10 −6 h/ p, which is small enough to guarantee that the algorithm truncation errors dominate the total numerical errors. All the computations were performed on an Intel ® Core(TM) i7-7700HQ CPU@2.80GHz with 8 GB of RAM. It should be pointed out that in all our simulations reported in this section, the number iterations in the linear solver do not exceed 30 iterations.

A Gaussian Pulse Example
In this example we consider the advection-diffusion of a Gaussian pulse in a rotating velocity field widely used in the literature to ascertain the performance of transport schemes, see for example [18,41]. Thus, the governing equations are of the form (33) with v = (−ωy, ωx) and ω = 4. Initial and boundary conditions are taken from the exact solution The purpose of this test example is to quantify the errors and convergence rates for the proposed Bernstein-Bézier Galerkin-characteristics finite element method. First we consider the case of pure advection for this example corresponding to ν = 0 in (33). In Table 1 we summarize the results obtained for the L 1 -error and convergence rate after one revolution using different meshes, polynomial degrees and CFL numbers. In this table we also include the CPU times for each run. In terms of L 1 -error, keeping the polynomial degree p fixed and refining the spatial step h results in a substantial decrease in the computed L 1 -errors. From the values of convergence rates in Table 1 we observe that the expected order of convergence is achieved for each selected polynomial degree p. It has also been observed that these convergence rates have not been deteriorated by the increase in the values of CFL, and the order of the Bernstein-Bézier Galerkin-characteristics finite element method remains almost the same for the considered CFL numbers of 2.5, 5 and 10. Notice that to reduce the computational cost, these CFL numbers are chosen as large as possible which yield explicit Eulerian-based methods noncompetitive.
Next we include the physical diffusion in this problem by solving the advection-diffusion of the Gaussian pulse with diffusion coefficient ν = 10 −6 . As in the previous case, Table 2 summarizes the results obtained for the L 1 -error and convergence rate after one revolution using different meshes, polynomial degrees and CFL numbers. It is clear from the results in Table 2 that for the considered diffusion coefficient, the expected convergence rates are preserved in the Bernstein-Bézier Galerkin-characteristics finite element method for the considered CFL numbers of 2.5, 5 and 10. Notice that, to be confined with convection-dominated problems only small values of the diffusion coefficient are accounted for in our simulations. Furthermore, for the considered small values of the diffusion coefficient, the time steps obtained according to the selected CFL numbers are small enough such that the errors associated with the discretization of the advection term are the dominant ones. Therefore, no substantial decrease in the computed errors is achieved when decreasing the CFL number. From the computational results obtained for the advection-diffusion of a Gaussian pulse, one may conclude the following: (i) the Bernstein-Bézier Galerkin-characteristics finite element method highly solve this test problem on coarse meshes, (ii) the convergence rate of the method is not deteriorated when increasing the CFL numbers.

Time-Dependent Burgers Example
To further quantify the errors for the proposed Bernstein-Bézier Galerkin-characteristics finite element method for time-dependent problems, we solve the following nonlinear Burgers problem in the squared domain = [0, 2] × [0, 2] subject to initial and boundary conditions obtained from the following exact solution .
It is clear that the velocity field in this example depends on the time for which time accuracy of the proposed method can be evaluated. A similar example has been considered in [27] using a spectral volume method. In Table 3 we present the results obtained for the L 1 -error, L 2 -error and convergence rates at time t = 2 for ν = 10 −2 using different meshes, different polynomial degrees and two CFL numbers CFL = 2.5 and CFL = 5. As in the previous simulations, keeping the polynomial degree p fixed and refining the mesh size h results in a substantial decrease in the computed L 1 -error, L 2 -error for both selected CFL numbers. The convergence rates in Table 3 confirm the expected order of convergence of the proposed Bernstein-Bézier Galerkin-characteristics finite element method for each selected polynomial degree p. It should be also noted that lower convergence rates have been observed for this test example using CFL = 2.5 than the case using CFL = 5. These features can be attributed to the nonlinear nature of the problem and also to the time dependence of the velocity field involved in the calculation of departure points. The next idea is to check the convergence rates in time of the proposed Bernstein-Bézier Galerkin-characteristics finite element method for this test example. In doing so, we compute the L 1 -error and L 2 -error using a fixed mesh with h = 1 32 and carry out some numerical experiments varying the polynomial degree and the time step t. The obtained results at time t = 2 for two different diffusion coefficient ν = 10 −2 and ν = 10 −3 are listed in Table 4. It is clear that decreasing the time step t yields a decrease in both L 1 -error and L 2 -error for all considered polynomial degrees p and diffusion coefficients. A second-order convergence in time is also clearly observed in Table 4 for the considered polynomial degrees p and diffusion coefficients. Note that for this test example, better convergence results are obtained for the simulations with ν = 10 −2 than the simulations with ν = 10 −3 .

Deformational Flow Example
In this example we solve the deformational flow problem widely used in the literature to examine the performance of Galerkin-characteristics methods, see for example [14,18,34]. Here, the problem statement consists of the linear advection eq. (10) in a circular spatial domain centered at (x 0 = 0, y 0 = 0) with radius 4 and equipped with a highly deformational velocity defined by a steady circular vortex with tangential velocity depending on the radius of the vortex as where v 0 = 2.58 and the radial distance R = (x − x 0 ) 2 + (y − y 0 ) 2 is the radius of the vortex. Initial and boundary conditions are obtained from the analytical solution with ω = v t R is the angular velocity of the circular vortex and η is a parameter controlling the steepness of the solution (43). Here, η = 0.05 which yields a tight hyperbolic tangent profile in the initial condition that results in a non-smooth solution as the problem is integrated forward in time. Figure 2 illustrates the computational mesh and the initial condition used in the simulations for this test example. The steep gradient at the centerline of the computational domain can be clearly seen in the initial condition. The unstructured mesh shown in Fig. 2 contains 237 elements. This coarse mesh is considered in the simulations to demonstrate the high accuracy of the proposed Bernstein-Bézier Galerkin-characteristics finite element method to solve highly deformational flow problems without need to very fine meshes.
In Fig. 3 we display the results obtained at time t = 4 and CFL = 3 using different values for the polynomial degree p. We also include the analytical solution in this figure for comparison reasons. For a better insight, Fig. 4 exhibits cross-sections at x = 0 and y = 0 of the results using the considered polynomial degrees. It is clear that, using p = 2 in the Bernstein-Bézier Galerkin-characteristics finite element method, nonphysical oscillations are generated in the obtained results, specially at the center of the spatial domain where the gradient is sharper. Increasing the value of polynomial degree to p = 4, the amplitude of these oscillations is reduced but the numerical dissipation is still present in these results. From the same figures we observe a complete absence of these oscillations and numerical diffusion in the results obtained using p = 8 and p = 10. Observe the good agreement between the   Fig. 4 even when using the coarse mesh and the large CFL number. Since the analytical solution for this flow problem is provided in (43), errors can be quantified. In Table 5 we present the L 1 -error and L 2 -error obtained using three meshes and different polynomial degrees. These unstructured meshes are used such that the number of elements in one mesh is about the double of its coarser counterpart. It is clear that increasing the polynomial order or refining the mesh result in a decrease in both L 1 -error and L 2 -error. This confirms the high accuracy and the good convergence of the proposed method for solving deformational flow problems. It should be pointed out that the performance of the proposed Bernstein-Bézier Galerkin-characteristics finite element method is very attractive  since the computed solutions remain stable and highly accurate even when coarse meshes are used without requiring nonlinear solvers or small time steps to be taken in the simulations.

Moving Fronts Example
We consider the problem of moving fronts modeled by the convection-diffusion eq. (33) equipped with a time-dependent velocity field. Initially, two separate fronts propagate along the main diagonal of the computational domain at different speeds and eventually coalesce into a single front for a longer time. This problem has also been solved in [30]  First we examine the mesh convergence in the proposed Bernstein-Bézier Galerkincharacteristics finite element method. To this end, we consider three structured meshes Mesh A, Mesh B and Mesh C with different element densities as depicted in Fig. 5. Here, the number of triangular elements in Mesh A, Mesh B and Mesh C is 32, 128 and 512, respectively. In Fig. 6 we present 20 equi-distributed contourlines of the solutions obtained at t = 0.2 using CFL = 3 and p = 10. Here, results are presented for four different diffusion coefficients namely, ν = 5 × 10 −4 , ν = 10 −4 , ν = 5 × 10 −3 and ν = 10 −2 . The analytical solutions are also included in Fig. 6 for comparison purposes. It is clear that, by decreasing the values of the physical diffusion ν, the convective term becomes dominant and steep internal layers are formed near the vicinity of front lines in the computational domain. The mesh convergence on these results is clearly shown for the Bernstein-Bézier Galerkin-characteristics finite element It is worth remarking that the thinning of the internal layers with decreasing ν is evident from these plots and the rate of this thinning is slower for ν = 5 × 10 −4 than for ν = 10 −3 . These features clearly demonstrate the high accuracy achieved by the proposed Bernstein-Bézier Galerkin-characteristics finite element method for solving moving fronts problems using coarse meshes and large time steps. In addition, compared to the results published for example in [30], it can be seen that our method resolves accurately the solution features and the moving fronts seem to be localized in the correct place in the flow domain without requiring very fine meshes as those reported in [30].
To further qualify the results for these meshes, we plot in Fig. 7 the cross-section results at the main diagonal y = 1 − x using ν = 5 × 10 −4 on the considered three meshes using different polynomial degrees. Those results obtained using ν = 10 −3 are depicted in Fig. 8.  For large values of p, it is clear that the method on the three meshes produces practically identical results. This can be attributed to the high accuracy achieved by Bernstein-Bézier Galerkin-characteristics finite element method using p = 12 even on the coarse mesh Mesh A. However, using low polynomial degrees on the meshes Mesh A and Mesh B, nonphysical oscillations are clearly visible in those regions where the solution gradients are steep such as the moving fronts. Apparently, the solution structures are in good agreement with the exact solutions presented in these figures. Needless to say that for convection-dominated

Viscous Burgers Flow Example
Our final example solves the following viscous Burgers equation which evolves to a highly convective steady-state where λ is a constant controlling the magnitude of the nonlinear convective term. The boundary conditions are of Dirichlet type given by the exact steady-state solution This problem has been solved in a squared domain in [30] using a modified method of characteristics. Here, this problem is solved in a circular domain with radius 5 and centered at the origin. The domain dimensions and the initial conditions used in our computations are illustrated in Fig. 9. We use the Bernstein-Bézier Galerkin-characteristics finite element method to compute the steady-state solutions for three different values of λ namely, λ = 1, λ = 5 and λ = 10. For these steady-state solutions, the time integration process is stopped when the inequality is satisfied. Here · denotes the L 1 -norm and τ is a given tolerance fixed to 10 −6 t in our computations. It should be pointed out that, the number of iterations to reach this tolerance depends on the values taken by λ such that, for a fixed CFL number, more iterations are required for larger values of λ. Figure 10 illustrates the obtained results using a mesh of 432 elements and CFL = 3. Here, we plot 10 equi-distributed contours of the solutions using two values of the polynomial degree p = 6 and p = 10. For comparison, we have also included analytical steady-state Fig. 10 Results for the viscous Burgers equation with λ = 1 (first column), λ = 5 (second column) and λ = 10 (third column) using p = 6 (first row), p = 10 (second row) and exact solution (third row) Fig. 11 Cross-sections at y = x of the results for the viscous Burgers problem using λ = 1 (first column), λ = 5 (second column) and λ = 10 (third column) solutions in this figure. It is clear that, by increasing the values of λ, the convective term becomes dominant and steep internal boundary layers are formed near the vicinity of center lines. For low values of λ, the internal boundary layers are wide and diffuse in the flow domain. As λ increases, the internal boundary layers concentrate and move towards the domain center. It is apparent that the flow structure is in good agreement with the previous work in [30]. These plots give a clear view of the overall flow pattern and the effect of the convection control parameter λ on the structure of steady internal boundary layers in the cavity. It is worth remarking that the thinning of the internal boundary layers with increasing λ is evident from these plots, although the rate of this thinning is slower for p = 6 method than for p = 10. These features clearly demonstrate the high accuracy achieved by the proposed Bernstein-Bézier Galerkin-characteristics finite element method for solving viscous Burgers problems at steady-state regimes. In addition, compared to the results published for example in [30], it can be seen that our method resolves accurately the flow structures and the boundary layers seem to be localized in the correct place in the flow domain.
For visualizing the comparisons, we display in Fig. 11 a cross-section at the main diagonal y = x using the selected values of λ and p. For λ = 1, it is clear that the results obtained using p = 6 and p = 10 produce practically identical results. This can be attributed to the large physical diffusion presented in the problem. However, increasing the value of λ to 5 and 10 the results computed using p = 10 are more accurate than those computed using p = 6. Apparently, By using p = 10, high resolution is achieved in those regions where the flow gradients are steep such as the moving fronts. Comparing the results obtained using the considered values of p, it is clear that using p = 6 produces diffusive solutions resulting in smearing the shocks. On the other hand, this numerical diffusion has remarkably been reduced in the results computed using p = 10.

Conclusions
In this paper, we have presented a Bernstein-Bézier Galerkin-characteristics finite element method for solving convection-diffusion equations on unstructured meshes. The proposed solver inherits the advantages of semi-Lagrangian integration (unconditional stability and reduction of time truncation errors) while preserving high-order accuracy. Moreover, this implementation simplifies the computational complexity of the Bernstein-Bézier Galerkincharacteristics solver for a fixed polynomial approximation during the time integration process. To increase the performance of the proposed method the element-level static condensation is implemented. The favorable performance of the developed algorithm is demonstrated using a series of numerical examples including the viscous Burgers equation. The results obtained for these examples show that the Bernstein-Bézier Galerkin-characteristics finite element method has the advantage of requiring less computational resources for the integration of convection-dominated problems than an Eulerian-based finite element methods, typical of those widely used in Eulerian algorithms for convection-dominated problems. This fact, as well as its favorable high-order accuracy and strong stability properties, make it an attractive alternative for convection solvers based on Galerkin-characteristic techniques. Finally, the method presented in this paper is not appropriate for nonlinear invscid equations with strong shocks. However, using the idea of the limiting techniques studied in [17], and following the arguments of Sect. 3, it is possible to reconstruct a shock-preserving Bernstein-Bézier Galerkin-characteristics finite element method for nonlinear equations. Results on these schemes will be reported in the near future. Future work will also concentrate on extending these methods for the incompressible viscous flows in three space dimensions using unstructured meshes and developing highly efficient preconditioned iterative solvers for the associated linear systems. Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.