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

Authors: Ali R. Soheili; J. Naghipoor; S. A. Ahmadian
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

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.


---

## Body

## 1. 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 time-dependent 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 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.


## 2. The GWMFE in 2D

Our formulation of GWMFE has been designed to solve a PDE of the type(2.1)∂u∂t=L(u),
where
(2.2)L(u)=F(x,y,t,u,∂u∂x,∂u∂y)∂2u∂x2+G(x,y,t,u,∂u∂x,∂u∂y)∂2u∂y2+H(x,y,t,u,∂u∂x,∂u∂y)
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)=u0(x,y), (x,y)∈Ω.


### 2.1. 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Ω=⋃j=1neΩj, where ne 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 Uj(x,y,t) to uj(x,y,t) as
(2.3)Uj(x,y,t)=∑k=1npbj,k(x,y)Ujk(t),
where bj,k(x,y) is the kth Lagrange basis function on the jth element, np is the number of interpolations points in an element, and Ujk 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 Uj(x,y,t), the physical coordinates on triangular element are introduced. So let (L1,L2,L3) be the physical coordinates of (x,y)∈Ωj. The Lagrangian interpolation functions are given by
(2.4)bj,k(x,y)=B(k)(L1)B(k)(L2)B(k)(L3),
where 1≤j≤ne, 1≤k≤np, and
(2.5)B(k)(LN)={∏s=1pLN(k)1s(pLN-s+1)pLN(k)≥1,1pLN(k)=0,
with N=1,2,3 and k=1,…,np that p denote the degree of approximation and LN(k) is Nth physical coordinate of the kth interpolation point.For example, for node 1 of triangle element with 6 nodes, we have(2.6)B(1)(L1)=∏s=12×11s(2L1-s+1)=L1(2L1-1),B(1)(L2)=1,B(1)(L3)=1,
then
(2.7)bj,1=L1(2L1-1).
For the node 2, we have(2.8)B(2)(L1)=∏s=12×(1/2)1s(2L1-s+1)=2L1,B(2)(L2)=∏s=12×(1/2)1s(2L2-s+1)=2L2,B(2)(L3)=1,
then
(2.9)bj,2=4L1L2.
For other nodes, we have
(2.10)bj,3=L2(2L2-1),bj,4=4L2L3,bj,5=L3(2L3-1),bj,6=4L1L3.
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.11)(f,g)w=∫Ωf(x)g(x)w(x)dΩ.
In two dimension, a standard form of weighted function is the positive function w(x,y), in meaning of the Galerkin approach, given rise to minimization
(2.12)minU̇i,Ẋi,Ẏi∬ΩR2w(x,y)dΩ,
where(2.13)R=∂U∂t-F(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2-G(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2-H(x,y,t,U,∂U∂x,∂U∂y).
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 be1/1+∥∇U∥22 and we will have
(2.14)minU̇i,Ẋi,Ẏiϕ=∬ΩR21+∥∇U∥22dΩ,
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∑j=1ne(εjΔ̇j-Sj) to the objective function (2.14) in order to prevent singularities depending on the area of each element Δj, j=1,…,ne. So we will have
(2.15)minU̇i,Ẋi,Ẏiϕ=∬ΩR21+∥∇U∥22dΩ+∑j=1ne(εjΔ̇j-Sj).
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)Sj=C1Δj-C3(1+C3Δj-C3)2,εj=(C1Δj-C3+C2)(1+C3Δj-C3)2,j=1,…,ne. The penalty constants C1, C2, and C3 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.


### 2.2. Time Derivative ofU

The approximationU(x,y,t) to u(x,y,t) is dependent on the nodal amplitudes Ujk and on the nodal position ξℓ=(Xℓ,Yℓ). So we can write
(2.17)∂U∂t=∑j=1ne∑k=1np∂U∂UjkU̇jk+∑ℓ=1ns∂U∂XℓẊℓ+∑ℓ=1ns∂U∂YℓẎℓ.
In order to define the system of ODEs generated by space discretization, it is necessary to evaluate the derivatives ∂U/∂Ujk for all 1≤j≤ne, 1≤k≤np, and ∂U/∂Xℓ and ∂U/∂Yℓ for all 1≤ℓ≤ns.Let us consider a global nodevi, which is either a vertex of a triangle or an interpolation point that belongs to an edge or is placed inside the triangle. Let ng be the number of these global nodes. The support Ωi of a global node vi is the union of triangles to whom vi belongs, say Ωi=⋃s=1I(i)Ωis, I(i) be the number of elements surrounding node vi. Suppose that vi is the global node associated with Pjk, the kth interpolation point of jth element of Ωi. Defining the global basis function φi, 1≤i≤ng as
(2.18)φi={0(x,y)∉Ωi,bj,k(x,y)∈Ωj⊂Ωi.
After some simple computations, from (2.3) and (2.18), we have
(2.19)∂U∂Ujk=φi,
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≤ℓ≤ns are similar. Consider
(2.20)∂U∂Xℓ=∑k=1np∂bk∂Xℓ(x,y)Ujk,∂U∂Yℓ=∑k=1np∂bk∂Yℓ(x,y)Ujk.
After some computations [22], we have
(2.21)∂U∂Xℓ=-Lk1∂U∂x,∂U∂Yℓ=-Lk1∂U∂y,
where Lk1 is the K1th physical coordinates of (x,y)∈Ωj.


### 2.3. The GWMFE Equations

Our GWMFE discretization leads to a large ODE system(2.22)∂ϕ∂U̇jk=0,∂ϕ∂Ẋℓ=0,∂ϕ∂Ẏℓ=0,
in which 1≤j≤ne, 1≤k≤np, and 1≤ℓ≤ns. Now consider the first equations in (2.22). From (2.15) and (2.17), we have
(2.23)∬ΩR21+∥∇U∥22∂U∂UjkdΩ=0,
for 1≤j≤ne, 1≤k≤np. Considering the support of the global node vi, say Ωi, associated with the kth interpolation point of the jth element we have(2.24)∬Ωi11+∥∇Ui∥22φi∂U∂tdΩ=∬Ωi11+∥∇Ui∥22φiF(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2dΩ+∬Ωi11+∥∇Ui∥22φiG(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2dΩ+∬Ωi11+∥∇Ui∥22φiH(x,y,t,U,∂U∂x,∂U∂y)dΩ,
for i=1,…,ng.For the second equation in (2.22), if we suppose that vi is the lth node space, then
(2.25)∬Ωi11+∥∇Ui∥22∂U∂Xℓ∂U∂tdΩ=∬Ωi11+∥∇Ui∥22∂U∂XℓF(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2dΩ+∬Ωi11+∥∇Ui∥22∂U∂XℓG(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2dΩ+∬Ωi11+∥∇Ui∥22∂U∂XℓH(x,y,t,U,∂U∂x,∂U∂y)dΩ-∑s=1I(i)12εisΔyk1(εisΔ̇is-Sis),
for ℓ=1,…,ns.Similarly, the third equation is(2.26)∬Ωi11+∥∇Ui∥22∂U∂Yℓ∂U∂tdΩ=∬Ωi11+∥∇Ui∥22∂U∂YℓF(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2dΩ+∬Ωi11+∥∇Ui∥22∂U∂YℓG(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2dΩ+∬Ωi11+∥∇Ui∥22∂U∂YℓH(x,y,t,U,∂U∂x,∂U∂y)dΩ-∑s=1I(i)12εisΔxk1(εisΔ̇is-Sis),
for ℓ=1,…,ns.


### 2.4. Second-Order Terms

Second order terms such as the LaplacianΔU=(∂2U/∂x2)+(∂2U/∂y2) 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, vi, we consider a δ-neighbourhood of an edge. In each of δ-neighbourhood of edge adjacent to the node vi, we need to calculate
(2.27)Iα=∬R{F(x,y,t,Ũ,∂Ũ∂x,∂Ũ∂y)∂2Ũ∂x2+G(x,y,t,Ũ,∂Ũ∂x,∂Ũ∂y)∂2Ũ∂y2}dxdy,
where R is the δ-neighbourhood of an edge adjacent to the node vi and Ũ is the recovered function of U. We take the 1D Hermite cubic recovered function that has C1 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
(2.28)Iα=∫0α∫-δδ{F(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)cos2(θ)+G(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)sin2(θ)}∂2Ũ∂ξ2dξdη,
where θ is the rotation angle. Let p=(0,ηα), in which 1≤ηα≤α and Ũ, the recovered function, defined by
(2.29)Ũ(ξ,ηα)={Uξ∉[-δ,δ],Hξ∈[-δ,δ],
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)Iα=α∑pw(p)(∂Ũ∂ξŨξ|p∈ΩR-∂Ũ∂ξŨξ|p∈ΩL)×{F(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)cos2(θ)+G(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)sin2(θ)}Ũξ|p,
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 vi.


### 2.5. General Equations of the GWMFE

Consider the global nodevi, i=1,…,ng, and its support Ωi. Assume that vi is the kisth interpolation point in the element Ωis, s=1,…,I(i). Denoting by ais the length of the sth edge adjacent to vi, the GWMFE equation (2.24) associated to the nodal amplitude is
(2.31)∑s=1I(i)∑k=1np∬Ωis11+∥∇Uis∥22bisbkdΩU̇jk+∑s=1I(i)∑N=13∬Ωis-11+∥∇Uis∥22LkNbis∂U∂xdΩẊℓN+∑s=1I(i)∑N=13∬Ωis-11+∥∇Uis∥22LkNbis∂U∂ydΩẎℓN=∑s=1I(i)∬Ωis11+∥∇Uis∥22bisF(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2dΩ+∑s=1I(i)∬Ωis11+∥∇Uis∥22bisG(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2dΩ+∑s=1I(i)∬Ωis11+∥∇Uis∥22bisH(x,y,t,U,∂U∂x,∂U∂y)dΩ+∑s=1I(i)∫0aisbkis(∂Ũ∂ξŨξ|p∈ΩR-∂Ũ∂ξŨξ|p∈ΩL)×{F(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)cos2(θ)+G(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)sin2(θ)}dη,
for i=1,…,ng.Ifξℓ=(Xℓ,Yℓ), ℓ=1,…,ns be a spatial node defining the k1th vertex of the element surrounding the node i. For ℓ=1,…,ns, (2.25) leads to(2.32)∑s=1I(i)∑k=1np∬Ωis-11+∥∇Uis∥22Lk1bk∂U∂xdΩU̇jk+∑s=1I(i)∑N=13∬Ωis11+∥∇Uis∥22LkNLk1(∂U∂x)2dΩẊℓN+∑s=1I(i)∑N=13∬Ωis11+∥∇Uis∥22LkNLk1∂U∂x∂U∂ydΩẎℓN+∑s=1I(i)12εis2Δyk1Δ̇is=∑s=1I(i)∬Ωis-11+∥∇Uis∥22Lk1∂U∂xF(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2dΩ+∑s=1I(i)∬Ωis-11+∥∇Uis∥22Lk1∂U∂xG(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2dΩ+∑s=1I(i)∬Ωis-11+∥∇Uis∥22Lk1∂U∂xH(x,y,t,U,∂U∂x,∂U∂y)dΩ+∑s=1I(i)12εisSisΔyk1+∑s=1I(i)∫0ais-Lk1{(12(∂Ũ∂ξŨξ|p∈ΩR+∂Ũ∂ξŨξ|p∈ΩL))cos(θ)-∂Ũ∂ηsin(θ)}Qdη,
and (2.26) leads to(2.33)∑s=1I(i)∑k=1np∬Ωis-11+∥∇Uis∥22Lk1bk∂U∂ydΩU̇jk+∑s=1I(i)∑N=13∬Ωis11+∥∇Uis∥22LkNLk1∂U∂x∂U∂ydΩẊℓN+∑s=1I(i)∑N=13∬Ωis11+∥∇Uis∥22LkNLk1(∂U∂y)2dΩẎℓN+∑s=1I(i)12εis2Δxk1Δ̇is=∑s=1I(i)∬Ωis-11+∥∇Uis∥22Lk1∂U∂yF(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2dΩ+∑s=1I(i)∬Ωis-11+∥∇Uis∥22Lk1∂U∂yG(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2dΩ+∑s=1I(i)∬Ωis-11+∥∇Uis∥22Lk1∂U∂yH(x,y,t,U,∂U∂x,∂U∂y)dΩ+∑s=1I(i)12εisSisΔxk1+∑s=1I(i)∫0ais-Lk1{(12(∂Ũ∂ξŨξ|p∈ΩR+∂Ũ∂ξŨξ|p∈ΩL))cos(θ)-∂Ũ∂ηsin(θ)}Qdη,
for ℓ=1,…,ns, where
(2.34)Q=∑s=1I(i)∫0aisbkis(∂Ũ∂ξŨξ|p∈ΩR-∂Ũ∂ξŨξ|p∈ΩL)×{F(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)cos2(θ)+G(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)sin2(θ)}dη.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:(2.35)M(t,Y)Ẏ=F(t,Y),
where Ẏ=[U̇1,U̇2,…,U̇n;ẋ1,ẋ2,…,ẋn;ẏ1,ẏ2,…,ẏ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.


## 2.1. 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Ω=⋃j=1neΩj, where ne 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 Uj(x,y,t) to uj(x,y,t) as
(2.3)Uj(x,y,t)=∑k=1npbj,k(x,y)Ujk(t),
where bj,k(x,y) is the kth Lagrange basis function on the jth element, np is the number of interpolations points in an element, and Ujk 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 Uj(x,y,t), the physical coordinates on triangular element are introduced. So let (L1,L2,L3) be the physical coordinates of (x,y)∈Ωj. The Lagrangian interpolation functions are given by
(2.4)bj,k(x,y)=B(k)(L1)B(k)(L2)B(k)(L3),
where 1≤j≤ne, 1≤k≤np, and
(2.5)B(k)(LN)={∏s=1pLN(k)1s(pLN-s+1)pLN(k)≥1,1pLN(k)=0,
with N=1,2,3 and k=1,…,np that p denote the degree of approximation and LN(k) is Nth physical coordinate of the kth interpolation point.For example, for node 1 of triangle element with 6 nodes, we have(2.6)B(1)(L1)=∏s=12×11s(2L1-s+1)=L1(2L1-1),B(1)(L2)=1,B(1)(L3)=1,
then
(2.7)bj,1=L1(2L1-1).
For the node 2, we have(2.8)B(2)(L1)=∏s=12×(1/2)1s(2L1-s+1)=2L1,B(2)(L2)=∏s=12×(1/2)1s(2L2-s+1)=2L2,B(2)(L3)=1,
then
(2.9)bj,2=4L1L2.
For other nodes, we have
(2.10)bj,3=L2(2L2-1),bj,4=4L2L3,bj,5=L3(2L3-1),bj,6=4L1L3.
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.11)(f,g)w=∫Ωf(x)g(x)w(x)dΩ.
In two dimension, a standard form of weighted function is the positive function w(x,y), in meaning of the Galerkin approach, given rise to minimization
(2.12)minU̇i,Ẋi,Ẏi∬ΩR2w(x,y)dΩ,
where(2.13)R=∂U∂t-F(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2-G(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2-H(x,y,t,U,∂U∂x,∂U∂y).
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 be1/1+∥∇U∥22 and we will have
(2.14)minU̇i,Ẋi,Ẏiϕ=∬ΩR21+∥∇U∥22dΩ,
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∑j=1ne(εjΔ̇j-Sj) to the objective function (2.14) in order to prevent singularities depending on the area of each element Δj, j=1,…,ne. So we will have
(2.15)minU̇i,Ẋi,Ẏiϕ=∬ΩR21+∥∇U∥22dΩ+∑j=1ne(εjΔ̇j-Sj).
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)Sj=C1Δj-C3(1+C3Δj-C3)2,εj=(C1Δj-C3+C2)(1+C3Δj-C3)2,j=1,…,ne. The penalty constants C1, C2, and C3 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.


## 2.2. Time Derivative ofU

The approximationU(x,y,t) to u(x,y,t) is dependent on the nodal amplitudes Ujk and on the nodal position ξℓ=(Xℓ,Yℓ). So we can write
(2.17)∂U∂t=∑j=1ne∑k=1np∂U∂UjkU̇jk+∑ℓ=1ns∂U∂XℓẊℓ+∑ℓ=1ns∂U∂YℓẎℓ.
In order to define the system of ODEs generated by space discretization, it is necessary to evaluate the derivatives ∂U/∂Ujk for all 1≤j≤ne, 1≤k≤np, and ∂U/∂Xℓ and ∂U/∂Yℓ for all 1≤ℓ≤ns.Let us consider a global nodevi, which is either a vertex of a triangle or an interpolation point that belongs to an edge or is placed inside the triangle. Let ng be the number of these global nodes. The support Ωi of a global node vi is the union of triangles to whom vi belongs, say Ωi=⋃s=1I(i)Ωis, I(i) be the number of elements surrounding node vi. Suppose that vi is the global node associated with Pjk, the kth interpolation point of jth element of Ωi. Defining the global basis function φi, 1≤i≤ng as
(2.18)φi={0(x,y)∉Ωi,bj,k(x,y)∈Ωj⊂Ωi.
After some simple computations, from (2.3) and (2.18), we have
(2.19)∂U∂Ujk=φi,
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≤ℓ≤ns are similar. Consider
(2.20)∂U∂Xℓ=∑k=1np∂bk∂Xℓ(x,y)Ujk,∂U∂Yℓ=∑k=1np∂bk∂Yℓ(x,y)Ujk.
After some computations [22], we have
(2.21)∂U∂Xℓ=-Lk1∂U∂x,∂U∂Yℓ=-Lk1∂U∂y,
where Lk1 is the K1th physical coordinates of (x,y)∈Ωj.


## 2.3. The GWMFE Equations

Our GWMFE discretization leads to a large ODE system(2.22)∂ϕ∂U̇jk=0,∂ϕ∂Ẋℓ=0,∂ϕ∂Ẏℓ=0,
in which 1≤j≤ne, 1≤k≤np, and 1≤ℓ≤ns. Now consider the first equations in (2.22). From (2.15) and (2.17), we have
(2.23)∬ΩR21+∥∇U∥22∂U∂UjkdΩ=0,
for 1≤j≤ne, 1≤k≤np. Considering the support of the global node vi, say Ωi, associated with the kth interpolation point of the jth element we have(2.24)∬Ωi11+∥∇Ui∥22φi∂U∂tdΩ=∬Ωi11+∥∇Ui∥22φiF(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2dΩ+∬Ωi11+∥∇Ui∥22φiG(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2dΩ+∬Ωi11+∥∇Ui∥22φiH(x,y,t,U,∂U∂x,∂U∂y)dΩ,
for i=1,…,ng.For the second equation in (2.22), if we suppose that vi is the lth node space, then
(2.25)∬Ωi11+∥∇Ui∥22∂U∂Xℓ∂U∂tdΩ=∬Ωi11+∥∇Ui∥22∂U∂XℓF(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2dΩ+∬Ωi11+∥∇Ui∥22∂U∂XℓG(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2dΩ+∬Ωi11+∥∇Ui∥22∂U∂XℓH(x,y,t,U,∂U∂x,∂U∂y)dΩ-∑s=1I(i)12εisΔyk1(εisΔ̇is-Sis),
for ℓ=1,…,ns.Similarly, the third equation is(2.26)∬Ωi11+∥∇Ui∥22∂U∂Yℓ∂U∂tdΩ=∬Ωi11+∥∇Ui∥22∂U∂YℓF(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2dΩ+∬Ωi11+∥∇Ui∥22∂U∂YℓG(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2dΩ+∬Ωi11+∥∇Ui∥22∂U∂YℓH(x,y,t,U,∂U∂x,∂U∂y)dΩ-∑s=1I(i)12εisΔxk1(εisΔ̇is-Sis),
for ℓ=1,…,ns.


## 2.4. Second-Order Terms

Second order terms such as the LaplacianΔU=(∂2U/∂x2)+(∂2U/∂y2) 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, vi, we consider a δ-neighbourhood of an edge. In each of δ-neighbourhood of edge adjacent to the node vi, we need to calculate
(2.27)Iα=∬R{F(x,y,t,Ũ,∂Ũ∂x,∂Ũ∂y)∂2Ũ∂x2+G(x,y,t,Ũ,∂Ũ∂x,∂Ũ∂y)∂2Ũ∂y2}dxdy,
where R is the δ-neighbourhood of an edge adjacent to the node vi and Ũ is the recovered function of U. We take the 1D Hermite cubic recovered function that has C1 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
(2.28)Iα=∫0α∫-δδ{F(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)cos2(θ)+G(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)sin2(θ)}∂2Ũ∂ξ2dξdη,
where θ is the rotation angle. Let p=(0,ηα), in which 1≤ηα≤α and Ũ, the recovered function, defined by
(2.29)Ũ(ξ,ηα)={Uξ∉[-δ,δ],Hξ∈[-δ,δ],
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)Iα=α∑pw(p)(∂Ũ∂ξŨξ|p∈ΩR-∂Ũ∂ξŨξ|p∈ΩL)×{F(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)cos2(θ)+G(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)sin2(θ)}Ũξ|p,
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 vi.


## 2.5. General Equations of the GWMFE

Consider the global nodevi, i=1,…,ng, and its support Ωi. Assume that vi is the kisth interpolation point in the element Ωis, s=1,…,I(i). Denoting by ais the length of the sth edge adjacent to vi, the GWMFE equation (2.24) associated to the nodal amplitude is
(2.31)∑s=1I(i)∑k=1np∬Ωis11+∥∇Uis∥22bisbkdΩU̇jk+∑s=1I(i)∑N=13∬Ωis-11+∥∇Uis∥22LkNbis∂U∂xdΩẊℓN+∑s=1I(i)∑N=13∬Ωis-11+∥∇Uis∥22LkNbis∂U∂ydΩẎℓN=∑s=1I(i)∬Ωis11+∥∇Uis∥22bisF(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2dΩ+∑s=1I(i)∬Ωis11+∥∇Uis∥22bisG(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2dΩ+∑s=1I(i)∬Ωis11+∥∇Uis∥22bisH(x,y,t,U,∂U∂x,∂U∂y)dΩ+∑s=1I(i)∫0aisbkis(∂Ũ∂ξŨξ|p∈ΩR-∂Ũ∂ξŨξ|p∈ΩL)×{F(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)cos2(θ)+G(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)sin2(θ)}dη,
for i=1,…,ng.Ifξℓ=(Xℓ,Yℓ), ℓ=1,…,ns be a spatial node defining the k1th vertex of the element surrounding the node i. For ℓ=1,…,ns, (2.25) leads to(2.32)∑s=1I(i)∑k=1np∬Ωis-11+∥∇Uis∥22Lk1bk∂U∂xdΩU̇jk+∑s=1I(i)∑N=13∬Ωis11+∥∇Uis∥22LkNLk1(∂U∂x)2dΩẊℓN+∑s=1I(i)∑N=13∬Ωis11+∥∇Uis∥22LkNLk1∂U∂x∂U∂ydΩẎℓN+∑s=1I(i)12εis2Δyk1Δ̇is=∑s=1I(i)∬Ωis-11+∥∇Uis∥22Lk1∂U∂xF(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2dΩ+∑s=1I(i)∬Ωis-11+∥∇Uis∥22Lk1∂U∂xG(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2dΩ+∑s=1I(i)∬Ωis-11+∥∇Uis∥22Lk1∂U∂xH(x,y,t,U,∂U∂x,∂U∂y)dΩ+∑s=1I(i)12εisSisΔyk1+∑s=1I(i)∫0ais-Lk1{(12(∂Ũ∂ξŨξ|p∈ΩR+∂Ũ∂ξŨξ|p∈ΩL))cos(θ)-∂Ũ∂ηsin(θ)}Qdη,
and (2.26) leads to(2.33)∑s=1I(i)∑k=1np∬Ωis-11+∥∇Uis∥22Lk1bk∂U∂ydΩU̇jk+∑s=1I(i)∑N=13∬Ωis11+∥∇Uis∥22LkNLk1∂U∂x∂U∂ydΩẊℓN+∑s=1I(i)∑N=13∬Ωis11+∥∇Uis∥22LkNLk1(∂U∂y)2dΩẎℓN+∑s=1I(i)12εis2Δxk1Δ̇is=∑s=1I(i)∬Ωis-11+∥∇Uis∥22Lk1∂U∂yF(x,y,t,U,∂U∂x,∂U∂y)∂2U∂x2dΩ+∑s=1I(i)∬Ωis-11+∥∇Uis∥22Lk1∂U∂yG(x,y,t,U,∂U∂x,∂U∂y)∂2U∂y2dΩ+∑s=1I(i)∬Ωis-11+∥∇Uis∥22Lk1∂U∂yH(x,y,t,U,∂U∂x,∂U∂y)dΩ+∑s=1I(i)12εisSisΔxk1+∑s=1I(i)∫0ais-Lk1{(12(∂Ũ∂ξŨξ|p∈ΩR+∂Ũ∂ξŨξ|p∈ΩL))cos(θ)-∂Ũ∂ηsin(θ)}Qdη,
for ℓ=1,…,ns, where
(2.34)Q=∑s=1I(i)∫0aisbkis(∂Ũ∂ξŨξ|p∈ΩR-∂Ũ∂ξŨξ|p∈ΩL)×{F(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)cos2(θ)+G(ξ,η,t,Ũ,∂Ũ∂ξ,∂Ũ∂η)sin2(θ)}dη.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:(2.35)M(t,Y)Ẏ=F(t,Y),
where Ẏ=[U̇1,U̇2,…,U̇n;ẋ1,ẋ2,…,ẋn;ẏ1,ẏ2,…,ẏ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.


## 3. Local Time Step Refinement

Let time steps of the problem have the form(3.1){t0,t1,…,tk},
where tk≤Tf.Now, we apply the refinement scheme at each time step [17], for example, on the first time step [t0,t1]. Set Δt=((t1-t0)/l0), l0∈N,
(3.2)t0+(l0-i)Δt=t0+(l0-i)Δt,i=l0,l0-1,…,0.
Thus the time integration from t0 to t1 involves l0 sub step such that(3.3)t=t0→t=t0+Δt→t=t0+2Δt→,…,t=t0+l0Δt=t1
solves l0 ODE systems similar (2.35) to approximate U(x,y,t1), x(t1), and y(t1) with [t1,Y1]= ODE15S(@Function,[t0,t1)],Y0,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 thekth ODE system is the approximation of the (k-1)th ODE system at tk.Generally, suppose that we are at time levelt=tk and want to move toward t=tk+1, similarly, consider Δt=(tk+1-tk)/lk such that
(3.4)tk+(lk-i)Δt=tk+(lk-i)Δt,i=lk,lk-1,…,0.
The values of U(x,y,tk+1), x(tk+1) and y(tk+1) are obtained by solving lk ODE system
(3.5)t=tk→t=tk+Δt→t=tk+2Δt→,…,t=tk+lkΔt=tk+1.
This process continues until tk+1≤Tf.


### 3.1. Algorithm

The local time step refinement (LTSR) method may be derived as followsStep 1.
Setxi=i/n, i=1,2,…,n, and yj=j/m, j=1,2,…,m and obtain U(x,y,0) with initial value xi, yj, Y0=[U0,x0,y0] and set k=0.Step 2.
SetΔt=(tk+1-tk)/lk, lk∈N and s=1.Step 3.
Solve ODE system (2.35) at time level t=tk+sΔt as follows: [tk+sΔt,Yk+sΔt] = ODE15S(@Function,[tk+(s-1)Δt,tk+sΔt],Yk+(s-1)Δt,Option),the initial value of which in this step is obtained by solving ODE system (2.35) at t=tk+(s-1)Δt.Step 4.
SetUk+(s-1)Δt=Uk+sΔt, xk+(s-1)Δt=xk+sΔt, yk+(s-1)Δt=yk+sΔt, Δk+(s-1)Δt=Δk+sΔt, εk+(s-1)Δt=εk+sΔt, Sk+(s-1)Δt=Sk+sΔt, and s=s+1.Step 5.
Ifs≤lk, then go to Step 3, or else k=k+1.Step 6.
Iftk≤Tf, then go to Step 2, or else break.So our solution,Y=[U1,U2,…,Unm;x1,x2,…,xn;y1,y2,…ym] will be obtained in desired Tf.


## 3.1. Algorithm

The local time step refinement (LTSR) method may be derived as followsStep 1.
Setxi=i/n, i=1,2,…,n, and yj=j/m, j=1,2,…,m and obtain U(x,y,0) with initial value xi, yj, Y0=[U0,x0,y0] and set k=0.Step 2.
SetΔt=(tk+1-tk)/lk, lk∈N and s=1.Step 3.
Solve ODE system (2.35) at time level t=tk+sΔt as follows: [tk+sΔt,Yk+sΔt] = ODE15S(@Function,[tk+(s-1)Δt,tk+sΔt],Yk+(s-1)Δt,Option),the initial value of which in this step is obtained by solving ODE system (2.35) at t=tk+(s-1)Δt.Step 4.
SetUk+(s-1)Δt=Uk+sΔt, xk+(s-1)Δt=xk+sΔt, yk+(s-1)Δt=yk+sΔt, Δk+(s-1)Δt=Δk+sΔt, εk+(s-1)Δt=εk+sΔt, Sk+(s-1)Δt=Sk+sΔt, and s=s+1.Step 5.
Ifs≤lk, then go to Step 3, or else k=k+1.Step 6.
Iftk≤Tf, then go to Step 2, or else break.So our solution,Y=[U1,U2,…,Unm;x1,x2,…,xn;y1,y2,…ym] will be obtained in desired Tf.


## 4. 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 C1, C2, and C3. The value of C3 corresponds to the minimum allowed for an element area and in all computations we consider C3=10-5.


### 4.1. 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(4.1)∂u∂t=p(∂2u∂x2+∂2u∂y2)-u(∂u∂x+∂u∂y),0≤x,y≤1,t≥0.
We assume that initial and Dirichlet boundary condition are chosen to correspond to the analytic solution
(4.2)u(x,y,t)=11+e(y+x-t)/2p.
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.Figure2 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.Nodes movement for Burger equation fromt=0 to t=2 at some seconds.
(a)(b)(c)(d)(e)(f)Nodes distribution and solution of Burgers equation with quadratic elements for uniform initial mesh.
(a)(b)Nodes movement and related solution of Burgers equation with quadratic elements att=0.5 second, C1=10-3, C2=10-4, C3=10-5.
(a)(b)Nodes movement and related solution of Burgers equation with quadratic elements att=1 second, C1=10-3, C2=10-4, C3=10-5.
(a)(b)Nodes movement and related solution of Burgers equation with quadratic elements att=2 seconds, C1=10-3, C2=10-4, C3=10-5.
(a)(b)Let us consider to the average error,(4.3)Eave=∑ℓ=1ns|U(ξℓ)-u(ξℓ)|ns,
where ns is number of spacial nodes in Ω, ξℓ is a space node, U(ξℓ) is the GWMFE solution at ξℓ and u(ξℓ) is the exact solution at ξℓ.Table1 present CPU time, number of function evaluation (NFE) and average error (Eave) for GWMFE at t=2 seconds seconds with quadratic element and C1=10-3, C2=10-4, and C3=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.Table 1
CPU time, number of function evaluation (NFE), and average error (Eave) for GWMFE at t = 2 seconds with C1=10−3, C2=10−4, and C3=10−5 in some meshes.
Number of mesh nodesCPU timeNFEEave5×52423.4536003840.626328526×67687.3593757460.444845917×713839.75000010390.34352671


## 4.1. 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(4.1)∂u∂t=p(∂2u∂x2+∂2u∂y2)-u(∂u∂x+∂u∂y),0≤x,y≤1,t≥0.
We assume that initial and Dirichlet boundary condition are chosen to correspond to the analytic solution
(4.2)u(x,y,t)=11+e(y+x-t)/2p.
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.Figure2 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.Nodes movement for Burger equation fromt=0 to t=2 at some seconds.
(a)(b)(c)(d)(e)(f)Nodes distribution and solution of Burgers equation with quadratic elements for uniform initial mesh.
(a)(b)Nodes movement and related solution of Burgers equation with quadratic elements att=0.5 second, C1=10-3, C2=10-4, C3=10-5.
(a)(b)Nodes movement and related solution of Burgers equation with quadratic elements att=1 second, C1=10-3, C2=10-4, C3=10-5.
(a)(b)Nodes movement and related solution of Burgers equation with quadratic elements att=2 seconds, C1=10-3, C2=10-4, C3=10-5.
(a)(b)Let us consider to the average error,(4.3)Eave=∑ℓ=1ns|U(ξℓ)-u(ξℓ)|ns,
where ns is number of spacial nodes in Ω, ξℓ is a space node, U(ξℓ) is the GWMFE solution at ξℓ and u(ξℓ) is the exact solution at ξℓ.Table1 present CPU time, number of function evaluation (NFE) and average error (Eave) for GWMFE at t=2 seconds seconds with quadratic element and C1=10-3, C2=10-4, and C3=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.Table 1
CPU time, number of function evaluation (NFE), and average error (Eave) for GWMFE at t = 2 seconds with C1=10−3, C2=10−4, and C3=10−5 in some meshes.
Number of mesh nodesCPU timeNFEEave5×52423.4536003840.626328526×67687.3593757460.444845917×713839.75000010390.34352671


## 5. 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 two-dimensional 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.



---

*Source: 602712-2009-10-12.xml*