New elements with harmonic shape functions in adaptive mesh refinement

https://doi.org/10.1016/j.finel.2019.103366Get rights and content

Highlights

  • New elements with midside nodes are proposed based on harmonic shape functions.

  • The shape functions are represented using Fourier series.

  • The hanging nodes are treated very naturally.

  • The convergence rate of the method is same as standard Lagrange linear elements.

  • The method has a great potential in adaptive mesh refinement approaches.

Abstract

New elements are proposed in the present work based on the harmonic coordinates in which the shape functions satisfy the Laplace equation. The harmonic functions have some appealing characteristics that made it possible to define elements with arbitrary shape functions on the element boundaries and arbitrary node arrangement. In the present work, without loss of generality, we formulate a set of quadrilateral elements with midside nodes and piecewise linear shape functions along the element boundaries. The analytical solution of Laplace equation in rectangular domains are used here to represent the harmonic shape function as Fourier series. To derive the shape functions of different elements with different node arrangements a systematic approach is proposed to represent all of the shape functions using a small set of harmonic functions. Two model boundary value problems, Poisson equation and linear elasticity are used here to evaluate the proposed method. Patch test and convergence analysis are done for each model problems via some examples. The elements with midside nodes have great potential to be used in mesh adaptive solutions. Therefore, some adaptive mesh refinement examples are solved based on the ZZ error estimation and quadtree grid structure as the mesh refinement technique.

Introduction

In the standard FEM, the elements geometry includes triangle and quadrilateral in 2D and tetrahedral, hexahedra and wedge in 3D domains. The shape functions are constructed using low order polynomials and the adjacent elements must be conforming without hanging nodes. Although such a small set of element geometries and polynomial shape functions are acceptable for many applications, there is a growing need to relax such constraints and to have access to more complex elements geometry, more sophisticated shape functions and methods which treat hanging nodes.

Many researches have been done to decrease the constraints of the standard FEM and extend or facilitate its application in different engineering problems. Here, we would refer to some of these efforts. The first approach includes the methods in which the approximation functions are constructed without a predefined connectivity of the nodes. These methods are generally known as meshless methods and numerous techniques are used to construct the meshless approximation functions. A review on the different shape functions which are used in meshless methods is presented in Ref. [1]. The second approach is enriching the solution space by inserting special functions to the basis of the shape functions. These methods use the partition of unity property of the shape functions to insert special or handbook functions in the basis functions. For example, extended FEM [2] and generalized FEM [3] use such a technique. In the third approach, to decrease the mesh dependency of the FEM, the computational mesh is independently constructed with respect to the problem geometry and the domain boundary passes through the elements interior. Such techniques are known as the fictitious domain methods or methods based on non-boundary-fitted meshes and facilitate solution of variable domain problems in which the domain geometry evolves during the solution process e.g. refer to Refs. [4,5]. In the fourth approach, the constraint of conformity of the mesh is relaxed and non-conforming meshes with hanging nodes are allowed. This approach considerably facilitates successive mesh refining or coarsening in adaptive solutions [6]. The fifth approach includes methods that use the tools which are developed in the CAD technologies such as b-splines or NURBS to represent both the geometry and field variables in the same way e.g. the isogeometric analysis is the mainstream of this approach [[7], [8], [9]]. Finally, in the sixth approach, the shape functions are constructed based on polygonal or polyhedral elements. These methods relax the constraint of the standard FEM on element geometry and any polygonal shape in 2D or polyhedral shape in 3D can be used to formulate the FEM. In this approach, the shape functions are generally constructed based on barycentric coordinates such as Wachspress coordinates [10], mean value coordinates [11], maximum entropy coordinates [12], Sibson coordinates [13] and harmonic coordinates [14]. For a complete review and comparison of different barycentric coordinates refer to Ref. [15].

The important point that must be highlighted is that regarding the above approaches implies that the common objective of all of the above six categories is seeking new ideas and methods for deriving the shape functions in such a way that fulfill some specific requirements. In fact, looking for new ideas and methods for constructing the shape functions is an active field of research. It is also the main goal of the present work. The present paper would open a door to a new world where one can define elements with arbitrary variation of the field variable on the element boundaries. The harmonic functions are used here to formulate such elements.

Harmonic coordinates are one of the most appealing ideas to drive the approximation functions. It seems the first use of harmonic functions for 2D interpolation was published by Gordon and Wixom in 1974 [16]. The use of harmonic approximation functions for solution of boundary value problems delayed until 2007 which Joshi et al. used it in the field of computer graphics [17]. Joshi et al. [17]. provides a proof that such approximation functions form a partition of unity and exactly reproduce a linear field. These functions also possess the Kronecker delta property. Harmonic functions satisfy the elliptic Laplace equation with positive boundary conditions. Therefore, the approximation functions are non-negative over the element. Furthermore, the harmonic functions do not possess local minima or maxima on the interior of the element [17]. To the best knowledge of the authors, a few papers have been published based on the application of harmonic shape functions in the FEM. For example refer to Refs. [14,17,18].

The analytical solution of Laplace equation is not available in general polygonal or polyhedral domains and seeking numerical solutions is unavoidable in such condition. For instance, in Ref. [17] the finite difference method is used in each individual element to approximate the harmonic shape functions. In a similar manner, the method of fundamental solutions is used in Ref. [18], the boundary element method is used in Ref. [19], the complex variable boundary method is used in Ref. [20] and the FEM is used in Ref. [14] to obtain the harmonic shape functions numerically. In fact, this difficulty is due to this point that in all of these researches the shape functions are formulated in global coordinate system and the Laplace boundary value problem must be solved in each individual element based on its geometry in global coordinate system. This intensively increases the computational cost of evaluation of the shape function and definitely restricts the use of harmonic shape functions in real-world engineering problems.

It also must be denoted that the use of harmonic functions in all of the previously published works is restricted to development of polygonal or polyhedral elements whereas a different goal is followed in the present work. The harmonic functions are used here to introduce new elements with arbitrary shape functions on the element boundaries and also node arrangement. Without loss of generality, the present work considers 2D problems with quadrilateral elements with midside nodes and piecewise linear interpolation along the element boundaries. To avoid numerical solution of Laplace equation, the shape functions are derived in local coordinate system for a representative master element. Fortunately, the analytical solution of Laplace equation is available in some simple geometries such as rectangular and cuboidal domains and therefore the harmonic shape functions can be obtained analytically as Fourier series. In addition, a new procedure is proposed here for systematic computation of the shape functions. In this approach, a small set of harmonic functions are defined at first as the building blocks and then all of the shape functions of different elements with different node arrangements are derived using proper transformations of them.

Two model equations, the Laplace/Poisson equation and linear elasticity are selected here to evaluate the performance and the potentials of the proposed method. The patch test and convergence analysis are done for each model equation using some numerical examples. The elements with midside nodes treat hanging nodes naturally and therefore have a great potential to be used in mesh adaptive solutions. To evaluate the applicability of the proposed method in adaptive mesh refinement solutions the Zienkiewicz and Zhu error estimation [21] is used as the refinement criteria and quadtree grid structure is used to modify the mesh. Although the quadtree structure is extremely fast in mesh refinement, it generates nonconforming elements containing hanging nodes. It is shown here that the proposed method treats nonconforming elements very naturally and efficiently. A lot of development is still to be done in the applications of harmonic shape functions and we believe it can be an interesting new field of research.

In the rest of the present paper, the basic notion of root coordinate system and root functions are introduced as the building blocks of the harmonic shape functions. It followed by classification of quadrilateral elements with midside node and detailed derivation of root functions and construction of the shape functions. Two model equations, their discretization, and numerical integration are then given. The paper is finalized with gradient recovery, error estimation, and several numerical examples.

Section snippets

Root coordinate system and root function

In this section, two basic notions of root coordinate system and root functions are introduced. These concepts will be used later for the construction of shape functions. But, before introducing them, we should explain the basic characteristics of the shape functions. To have valid shape functions for the use in the FEM the following characteristics must be satisfied.

  • a

    Linear independence

  • b

    Inter element continuity

  • c

    First order differentiable

  • d

    Partition of unity

  • e

    Kronecker delta property

The harmonic shape

Classification of quadrilateral elements with midside nodes

The harmonic shape functions can be used potentially to formulate elements with any node arrangement but here without loss of generality we focus on a set of quadrilateral elements with midside nodes and linear variation on the element boundaries. In the most general case, the quadrilateral elements with midside nodes consist of 5 different element types which are introduced carefully in the present section.

All of the 5 different element types in the present work are defined in such a way that

Derivation of the root functions

The 5 different element types which are introduced in section 3 totally contain 32 nodes and each node has its own shape function. It will be shown later in section 5 that all of these shape functions can be constructed using only 4 root functions and proper coordinate transformations. To explain the structure of the RFs, consider Fig. 4. In this figure, the value of the four RFs on the boundaries of the square domain α,β[1,1] in the RCS are given. For example, the RF ϕ1 has a unit value at (

Harmonic shape functions

As stated before in section 2, the four shape functions of bilinear quadrilateral element Q4 which are given in Eq. (1) are obtained using the root function ϕ1 with proper coordinate transformations as given in Eq. (5). A similar approach is possible to obtain the shape functions of element types A to E. For instance, consider the element type A. One side of the element type A has one midside node. The shape function associated with each node must have unit value at its own node and zero value

Model equations and discretization

In the present work, two model equations: Laplace/Poisson equation and linear elasticity are considered as the model equations to evaluate the proposed method. The main steps of discretization process are similar to the standard FEM which are described briefly in this section with some points on domain integration.

Gradient recovery and error estimation

In the standard FEM with Lagrange based shape functions, sampling the gradients of the field variables at some particular points (superconvergent points) has a higher degree of accuracy. Therefore, the gradient recovery technique [24] computes the gradients of the field variables at these points (also are known as recovery points) and then extrapolate them to the nodes in each individual element. Finally, an averaging must be done between the extrapolated values in the elements that share a

Numerical examples

In the present section, several numerical examples of Poisson/Laplace equation and linear elasticity are solved to evaluate the proposed method. The relative L2 error norm is used to quantitatively evaluate the results. The relative L2 error norm for Laplace/Poisson equation and elasticity problem are computed based on Eqs (42), (43) respectively.ηL2=Ω(TTh)2dΩΩT2dΩηL2=Ω(uuh)T(uuh)dΩΩuTudΩwhere T is the exact field variable in Laplace/Poisson equation and u is the exact displacement

Conclusion

The harmonic functions were used in the present work to define the elements with arbitrary shape functions on the element boundaries. As an application of the present method a set of new elements with midside node were developed here. The shape functions were defined in such a way that vary linearly between adjacent nodes on the element boundaries and therefore the elements could be attached to any standard linear elements. As a result, the hanging nodes were treated very naturally in the

References (25)

  • E.L. Wachspress

    A Rational Finite Element Basis, Volume 114 of Mathematics in Science and Engineering

    (1975)
  • N. Sukumar

    Construction of polygonal interpolants: a maximum entropy approach

    Int. J. Numer. Methods Eng.

    (2004)
  • Cited by (0)

    View full text