New elements with harmonic shape functions in adaptive mesh refinement
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 in the RCS are given. For example, the RF 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 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 error norm is used to quantitatively evaluate the results. The relative error norm for Laplace/Poisson equation and elasticity problem are computed based on Eqs (42), (43) respectively.where is the exact field variable in Laplace/Poisson equation and 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)
- et al.
The generalized finite element method
Comput. Methods Appl. Mech. Eng.
(2001) Numerical flow simulation in gated hydraulic structures using smoothed fixed grid finite element method
Appl. Math. Comput.
(2014)Isogeometric analysis in solution of unconfined seepage problems
Comput. Math. Appl.
(2019)- et al.
Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement
Comput. Methods Appl. Mech. Eng.
(2005) Mean value coordinates
Comput. Aided Geomet. Des.
(2003)Meshfree Methods: Moving beyond the Finite Element Method
(2018)Extended Finite Element Method: for Fracture Analysis of Structures
(2008)Optimal shape design for heat conduction using smoothed fixed grid finite element method and modified firefly algorithm
Iran. J. Sci. Tech. Trans. Mech. Eng.
(2015)- et al.
An easy treatment of hanging nodes in hp-finite elements
Finite Elem. Anal. Des.
(2016) - et al.
Isogeometric Analysis: toward Integration of CAD and FEA
(2009)