A Gradient Weighted Moving Finite-Element Method with Polynomial Approximation of Any Degree

A gradient weighted moving finite element method (GWMFE) based on piecewise polynomial of any degree is developed to solve time-dependent problems in two space dimensions. Numerical experiments are employed to test the accuracy and effciency of the proposed method with nonlinear Burger equation.


Introduction
Many problems in science and engineering are formulated in terms of time-dependent partial differential equations PDEs .It is well known that due to the moving steep fronts present in the solution, these problems present serious numerical difficulties.We present an approach where the mesh moves dynamical to capture the sharp front with a small number of space nodes.
Moving finite-element method MFE is a discretization technique on continuously deforming spatial grids introduced by K. Miller and R. N. Miller 1, 2 to deal with timedependent partial differential equations involving fine scale phenomena such as moving fronts, pulses and shocks 3-7 .Much significant work on MFE has been done by Baines and Wathen and their collaborators 1, 4, 5, 8, 9 .In particular, we mention the Baines book 8 and its excellent bibliography.For more information about moving finite-element method and its aspects and applications, see 9-18 .In all of these works, the method is based on a minimization of the PDE residual that is obtained by approximating the solution with piecewise linear elements.In 19 ,Coimbra et. al. introduced the MFE in two dimensions in which the degree of approximation polynomial was greater than 1.In 6, 7 , Carlson and Miller introduced the GWMFE in two Mathematical Problems in Engineering dimension in which the approximation polynomial was linear.In this paper, we present a formulation of GWMFE with approximation of higher degree to solve two dimensional time-dependent partial differential equations.This formulation of the GWMFE, which builds on the original approach of Miller 1, 2 , uses piecewise polynomial approximations of any degree of the 2D spatial domain.Numerical investigations are presented to show the accuracy and effectiveness of our method.

Description of the Method
The GWMFE is a numerical procedure which allows the local gradient adaptation of the finite-element approximation space with time.For the space discretization, we consider a hexagonally connected triangularization of Ω n e j 1 Ω j , where n e is the number of elements on Ω.In each triangle Ω j , the solution approximated by a polynomial of degree greater than 1.We define the polynomial approximation U j x, y, t to u j x, y, t as where b j,k x, y is the kth Lagrange basis function on the jth element, n p is the number of interpolations points in an element, and U k j is the value of U at the kth interpolation point of the jth element.Because of minimizing the algebra requirements in the formulation of GWMFE for computing U j x, y, t , the physical coordinates on triangular element are introduced.So let L 1 , L 2 , L 3 be the physical coordinates of x, y ∈ Ω j .The Lagrangian interpolation functions are given by where 1 ≤ j ≤ n e , 1 ≤ k ≤ n p , and with N 1, 2, 3 and k 1, . . ., n p that p denote the degree of approximation and L k N is Nth physical coordinate of the kth interpolation point.
For example, for node 1 of triangle element with 6 nodes, we have For the node 2, we have For other nodes, we have A weighted form of the variational formulation is often recommended, in particular, when the method is overly sensitive to specific features of the physical problem such as steep fronts.Such weighting replaces the inner products f, g 2 by inner products with respect to a given, positive, weight function w : R → Ω.Then, in one dimension, we will have:

2.13
It is also possible to use a solution-dependent weight function w x, y , which depends on x, y or those known function, say U, or its first derivatives, say ∂U/∂x or ∂U/∂y.
In GWMFE, this weight function is taken to be 1/ 1 ∇U 2 2 and we will have where ∇U is the gradient with respect to the physical variable x and y.
The argument for the use of this weight function is that it de-emphasizes those parts of the integral where ∂U/∂x and ∂U/∂y are large and therefore reduce the effect of minimization in steep parts of the solution.These moving node methods are especially suited to problems which develop sharp moving fronts, especially problems where one needs to resolve the fine-scale structure of the fronts.
We add the penalization term n e j 1 ε j Δj − S j to the objective function 2.14 in order to prevent singularities depending on the area of each element Δ j , j 1, . . ., n e .So we will have Here Δj is the time derivative of Δ j .The internodal viscosity function and the internodal spring function associated to an element Δ j are chosen in a closed form to the original proposed by Miller 1 .We consider, respectively, 2.16 j 1, . . ., n e .The penalty constants C 1 , C 2 , and C 3 are small constants supplied by the user.Penalty functions do not interfere on the solution, but exclusively on the movement of the nodes in order to prevent singularities.Their disadvantage is that it is not possible to set up a priori a relation between them and the problem we are solving.
The discretization of space-variables transforms each PDE in a system of ODE.To accomplish the discretization of problem 2.1 some overwriting may be necessary in order to apply the appropriate boundary conditions.The full discretization of 2.1 is obtained solving the ODE system by a suitable ODE solver such as LSODI 20 in FORTRAN software or ODE15S 16, 21 in MATLAB software.

Time Derivative of U
The approximation U x, y, t to u x, y, t is dependent on the nodal amplitudes U k j and on the nodal position ξ X , Y .So we can write 2.17 In order to define the system of ODEs generated by space discretization, it is necessary to evaluate the derivatives ∂U/∂U k j for all 1 ≤ j ≤ n e , 1 ≤ k ≤ n p , and ∂U/∂X and ∂U/∂Y for all 1 ≤ ≤ n s .
Let us consider a global node v i , which is either a vertex of a triangle or an interpolation point that belongs to an edge or is placed inside the triangle.Let n g be the number of these global nodes.The support Ω i of a global node v i is the union of triangles to whom v i belongs, say Ω i I i s 1 Ω i s , I i be the number of elements surrounding node v i .Suppose that v i is the global node associated with P k j , the kth interpolation point of jth element of Ω i .Defining the global basis function ϕ i , 1 ≤ i ≤ n g as

2.18
After some simple computations, from 2.3 and 2.18 , we have where ϕ i is the global node associated to the kth interpolation point of jth element.The computation of ∂U/∂X and ∂U/∂X for 1 ≤ ≤ n s are similar.Consider

2.20
After some computations 22 , we have where L k 1 is the K 1 th physical coordinates of x, y ∈ Ω j .

The GWMFE Equations
Our GWMFE discretization leads to a large ODE system
For the second equation in 2.22 , if we suppose that v i is the lth node space, then ∂U ∂X ∂U ∂t dΩ ∂U ∂X H x, y, t, U, ∂U ∂x , ∂U ∂y dΩ Similarly, the third equation is ∂U ∂Y ∂U ∂t dΩ ∂U ∂Y H x, y, t, U, ∂U ∂x , ∂U ∂y dΩ 2.26 for 1, . . ., n s .

Second-Order Terms
Second order terms such as the Laplacian ΔU ∂ 2 U/∂x 2 ∂ 2 U/∂y 2 need to be interpreted in GWMFE in the sense of smoothing.That is, we imagine the corners of our GWMFE to be ever so slightly smoothed.
Based on the idea of smoothing, there are basically three techniques for dealing with this problem.The δ-mollification of Miller 1 , the application of Greens theorem to reduce the order of the differentiation 14 , and the idea of recovery 9 which requires constructing a function U x, y, t from the piecewise polynomial approximation U x, y, t with sufficient continuity to ensure that all the integrals involving second-order derivatives exist and may be evaluated.In order to define and evaluate the integrals involving second-order derivatives over the support, Ω i , of a global node, v i , we consider a δ-neighbourhood of an edge.In each of δ-neighbourhood of edge adjacent to the node v i , we need to calculate where R is the δ-neighbourhood of an edge adjacent to the node v i and U is the recovered function of U. We take the 1D Hermite cubic recovered function that has C 1 continuity 8, 23 .
With changing the coordinates, we can describe R in terms of δ and α, α the length of the edge, R −δ, δ × 0, α .Then we have where θ is the rotation angle.Let p 0, η α , in which 1 ≤ η α ≤ α and U, the recovered function, defined by

2.29
where H is the 1D Hermite cubic polynomial defined by the values of U and ∂U/∂ξ at −δ, η α and δ, η α .
So the integral 2.28 may be approximate without difficulties in the usual way by some quadrature rule, for example, the mid-point rule.Thus when δ → 0, with some computation 19 , we get

2.30
where w p is the 1D quadrature weight associated to the quadrature point p and where Ω R and Ω L denote the right and left triangle, respectively Figure 1 .
Figure 1: δ-neighbourhood of an edge adjacent to the node v i .

General Equations of the GWMFE
Consider the global node v i , i 1, . . ., n g , and its support Ω i .Assume that v i is the k i s th interpolation point in the element Ω i s , s 1, . . ., I i .Denoting by a i s the length of the sth edge adjacent to v i , the GWMFE equation 2.24 associated to the nodal amplitude is

2.31
for i 1, . . ., n g .If ξ X , Y , 1, . . ., n s be a spatial node defining the k 1 th vertex of the element surrounding the node i.For 1, . . ., n s , 2.25 leads to

2.34
However 2.31 , 2.32 , and 2.33 define the general form of the system of ODE generated by the GWMFE.According to the boundary conditions, some of the equations may have to be rewritten.
This system of ODE can be written as follows: where Ẏ U1 , U2 , . . ., Un ; ẋ1 , ẋ2 , . . ., ẋn ; ẏ1 , ẏ2 , . . ., ẏn .This system of ODE has a stiff mass matrix and appropriate methods are thus required.In the present work, we use the ODE15S package 16, 21 under All cab software for integrating in time.

Local Time Step Refinement
Let time steps of the problem have the form where t k ≤ Tf.Now, we apply the refinement scheme at each time step 17 , for example, on the first time step t 0 , t 1 .Set Δt Thus the time integration from t 0 to t 1 involves l 0 sub step such that solves l 0 ODE systems similar 2.35 to approximate U x, y, t 1 , x t 1 , and y t 1 with t 1 , Y 1 ODE15S @Function, t 0 , t 1 ,Y 0 ,Option .In each ODE system, we need the initial conditions which are obtained by solving the previous ODE system.In other words, the initial condition of the kth ODE system is the approximation of the k − 1 th ODE system at t k .
Generally, suppose that we are at time level t t k and want to move toward t t k 1 , similarly, consider Δt The values of U x, y, t k 1 , x t k 1 and y t k 1 are obtained by solving l k ODE system This process continues until t k 1 ≤ Tf.

Algorithm
The local time step refinement LTSR method may be derived as follows Step 1. Set x i i/n, i 1, 2, . . ., n, and y j j/m, j 1, 2, . . ., m and obtain U x, y, 0 with initial value x i , y j , Y 0 U 0 , x 0 , y 0 and set k 0.
Step 2. Set Δt t k 1 − t k /l k , l k ∈ N and s 1. Step 3. Solve ODE system 2.35 at time level t t k sΔt as follows: t k sΔt , Y k sΔt ODE15S @Function, t k s−1 Δt , t k sΔt ,Y k s−1 Δt ,Option , the initial value of which in this step is obtained by solving ODE system 2.35 at t t k s−1 Δt .
Step 5.If s ≤ l k , then go to Step 3, or else k k 1. Step

Numerical results
We present a numerical example to illustrate the performance of our GWMFE.The integrals that appear in the system of ODE, say 2.31 , 2.32 , and 2.33 , are evaluated by 1D and 2D Gauss quadrature formulaes with 7 interior points 24 .For integration in time, we have used the integrator ODE15S.We select the standard choices for both absolute and relative time tolerances for the ODE15S integrator, Atol Rtol 10 −8 .In order to define the penalty functions the user must supply the penalty constants C 1 , C 2 , and C 3 .The value of C 3 corresponds to the minimum allowed for an element area and in all computations we consider C 3 10 −5 .

Burgers Equation in 2D
Some of the more difficult and interesting real life problems in which adaptive algorithms are needed arise in transport phenomena in which steep fronts propagate through the domain.
The special case of the nonlinear Burger equation is often used to test numerical methods so we consider the nonlinear evolution equation We assume that initial and Dirichlet boundary condition are chosen to correspond to the analytic solution u x, y, t 1 1 e y x−t /2p .

4.2
This problem is solved from t 0 to t 2 with p 0.025 with quadratic polynomial as approximation function.Solution of this problem can be obtained with similar computational effort for smaller values of p.All numerical results shown here are obtained on a Pentium IV processor at 3.00 GHz.
Figure 2 presents the adaptivity and nodes movement for t 0 to t 2 seconds with quadratic elements at some cases.Figures 3-6 present nodes movement and their solutions with quadratic elements for t 0 second, t 0.5 second, t 1 second and t 2 seconds, respectively.
Let us consider to the average error, where n s is number of spacial nodes in Ω, ξ is a space node, U ξ is the GWMFE solution at ξ and u ξ is the exact solution at ξ .Table 1 present CPU time, number of function evaluation NFE and average error E ave for GWMFE at t 2 seconds seconds with quadratic element and C 1 10 −3 , C 2 10 −4 , and C 3 10 −5 as penalty constants in some meshes.As Table 1 shows, when number of nodes in each direction are increased, then NFE in ODE15S and therefore CPU time has been increased and our average error has been decreased.

Conclusions
In this paper, we presented a gradient weighted moving finite-element method based on polynomial approximations of high degree for the solution of time-dependent PDEs on twodimensional space domains.We used a solution-dependent weight function for original MFE formulation to have better performance and mesh adaptivity.These moving nodes method are especially suited to problems which develop sharp moving fronts, especially problems where one needs to resolve the fine-scale structure of the fronts.
A careful treatment of the general second order terms is carried out.Moreover, by using numerical evaluations of all integrals, we can solve a large class of problems without extra calculations.The GWMFE is applied to the Burger test equation for transport process with quadratic polynomial as interpolation function.One can solve this problem with other nonlinear approximation function as well as other penalty constants.Numerical results are given to illustrate the good behavior of the GWMFE when using some cases of penalty constants.

Figure 2 :
Figure 2: Nodes movement for Burger equation from t 0 to t 2 at some seconds.

Figure 3 :Figure 4 :
Figure 3: Nodes distribution and solution of Burgers equation with quadratic elements for uniform initial mesh.
Our formulation of GWMFE has been designed to solve a PDE of the type is a first or second order differential operator on the 2D spatial interval Ω, if F x, y, t, u, ∂u/∂x, ∂u/∂y and G x, y, t, u, ∂u/∂x, ∂u/∂y be zero or not, respectively, and under Dirichlet, Neumann, or Robin boundary conditions and initial conditions satisfying u x, y, 0 u 0 x, y , x, y ∈ Ω.

Table 1 :
CPU time, number of function evaluation NFE , and average error E ave for GWMFE at t 2 seconds with C 1 10 −3 , C 2 10 −4 , and C 3 10 −5 in some meshes.