Imposition of Dirichlet Boundary Conditions in Element Free Galerkin Method through an Object-Oriented Implementation

One of the main drawbacks of Element Free Galerkin (EFG) method is its dependence on moving least square shape functions which don’t satisfy the Kronecker Delta property, so in this method it’s not possible to apply Dirichlet boundary conditions directly. The aim of the present paper is to discuss different aspects of three widely used methods of applying Dirichlet boundary conditions in EFG method, called Lagrange multipliers, penalty method, and coupling with finite element method. Numerical simulations are presented to compare the results of these methods form the perspective of accuracy, convergence and computational expense. These methods have been implemented in an object oriented programing environment, called INSANE, and the results are presented and compared with the analytical solutions.


Latin American Journal of Solids and
Meshfree methods were introduced to eliminate part of those difficulties such as distorted elements, need for re-meshing and the other limitations which arise because of element connectivity.A large variety of meshfree methods have been developed up to now for solving different engineering problems, such as smooth particle hydrodynamics (SPH) by Lucy (1997), diffuse elements, Nayroles et al. (1992), Element-Free Galerkin method (EFG), Belytschko et al. (1994), reproducing kernel particle method (RKPM), Liu et al. (1995), and hp-clouds methods, Duarte and Oden (1995).
The EFG method is a very promising method for the treatment of partial differential equations.Because of the absence of element connectivity, node particles can be added without the need for remeshing, so adaptive refinement of the discretization can be performed easily.This makes the EFG method more flexible than finite element method, and prominent for solving problems involving large deformation, fracture and fragmentation.But EFG method still has some drawbacks compared with FE method.The most difficulties of this method are: 1) as a meshfree method, this method is also computationally expensive because of the complex meshfree interpolations and the complex implementation algorithm, 2) since this method is based on the moving least squares shape functions, it is difficult to impose Dirichlet boundary conditions.Considerable researches have been devoted to overcome these difficulties, such as Lagrange multipliers method, Belytschko et al. (1994), and the penalty method Zhu and Atluri (1998), based on a modification of weak form.Other methods were proposed that can be interpreted as a modification of the interpolation shape functions, see for instance, Belytschko et al. (1995b), Huerta et al. (2000).The Lagrange multipliers method, firstly introduced by Babuška (1973), has been thoroughly studied by several researchers in Pitkäranta (1979), Barbosa et al. (1992), Stenberg (1995).This method is used to approximate the Dirichlet boundary conditions by using some unknown parameters called Lagrange multipliers.It was proven that this method had optimal convergence rates for the finite element spaces.The penalty method requires only the choice of the penalty parameter.The large values of this parameter must be used in order to impose the Dirichlet boundary condition in a proper manner.In practice, that might lead to ill-conditioned systems of equations, reducing the applicability of this method.Further information on penalty method can be found in Zhu and Atluri (1998).
Coupling EFG with finite element in order to apply Dirichlet boundary conditions is one the methods that can be interpreted as a modification of the interpolation shape functions.In order to do so, several coupling techniques are proposed in different studies.Belytschko et al. (1995b), Krongauz and Belytschko (1996) developed an algorithm for coupling EFG and FEM by a mixed interpolation in the transition domain, where FE nodes are substituted by particles and connected via ramp functions to the EFG nodes.The drawback of this method is that the derivatives are discontinuous along the interface of EFG and FE domains.Huerta et al. (2000) developed a hierarchical approximation based on finite element and mesh-free methods, which is able to remove the discontinuities in the derivatives of shape functions across the inner boundaries that was the major drawback of the method Latin American Journal of Solids and Structures 14 (2017) 1017-1039 proposed by Belytschko[14].Another approach for coupling EFG and FE is by using Lagrange multipliers.In this approach the substitution of FE nodes by particles is not necessary, Hegen (1996), Rabczuk and Belytschko (2006).Another coupling approach using Lagrange multipliers was proposed by Belytschko and Xiao, called 'bridging domain coupling method, Belytschko and Xiao (2004).In their approach, FE nodes and meshfree particles could overlap.Gu and Zhang (2008) developed a new concurrent simulation technique to couple the EFG method with the finite element method for the analysis of crack tip fields, using the 'bridging domain coupling method' concept and Lagrange multipliers.Liu et al. (2008) presented an adaptive EFG-FE coupling method oriented for bulk metal forming.In their method the simulation is performed using finite elements initially.During the process, the FE region with severe deformation occurrence is automatically converted to an EFG region.
The application of FEM object-oriented programming has been receiving great attention over the last two decades.As a result, a bunch of FEM-based codes used object-oriented concept as their implementation strategy, such as: a C++ based finite element code developed by Patzak et al. (2001), a meshless Object-oriented programming (OOP) code using non-uniform rational b-spline, Zhang et al. (2006), Freefem++, Hecht (2012), a fem code to solve two-and three-dimensional problems, a general framework for weak-form meshfree methods by Hsieh et al. (2014), and an educational FE software for truss structures presented by Zou et al. (2014).
In the present work, the three introduced methods of applying Dirichlet boundary conditions are investigated through developing MATLAB codes and also implementation of an object-oriented programming environment, called INSANE (INteractive Structural ANalysis Environmenthttp://www.insane.dees.ufmg.br)platform, Alves et al. (2013), Fonseca and Pitangueira (2007).The main contribution of the present work is regarding to implementation which, to the best knowledge of the author, is the first object oriented program for meshfree methods and also EFG/FE coupling approach.The implemented EFG program is highly efficient with a wide range of applications since it' s built on the basis of a well-developed finite element program, not only in the sense of using integration cells but, indeed the ability of using several features of FEM code which are already implemented, or will be further implemented in INSANE computational platform.The most significant example for this is the capability of using existing physically nonlinear models of FEM code, such as smeared crack and damage models in order to solve continuum damage problems in meshfree framework.This paper is arranged as follows: first a brief review of EFG formulation will be presented, and the methods of applying Dirichlet boundary conditions which are studied in this work will be explained.Then the INSANE platform will be briefly introduced with emphasis on its meshfree modules.In Sect. 5 the numerical examples will be presented and the results will be compared with the analytical solutions.

Moving Least Square Approximation
The EFG method utilizes the moving least square approximants, which are constructed in terms of nodes only: Latin American Journal of Solids and Structures 14 (2017) 1017-1039 where p(x) is a complete polynomial basis of arbitrary order and are coefficients which are obtained at any point x by minimizing the following: in which nsn is the number of nodes in the support domain of point x, and w(x-xI) is considered as a cubic spline weight function defined by: Minimizing Eq. (1b) leads to: where:   , , , By substituting Eq. ( 3) into Eq.( 1a), the MLS approximants are defined as: where ∅ are the EFG shape functions and are defined as: Latin American Journal of Solids and Structures 14 (2017) 1017-1039

EFG with Lagrange Multipliers
For any trial function, ∈ , and test function, ∈ , the Galerkin weak form of the EFG method with Lagrange multipliers is: .
with the following boundary conditions: where U and S are the spaces of test and trial functions, respectively, Γ and Γ are the Dirichlet and Neumann boundaries, is Lagrange multiplier, and f is the excitation vector.
Applying the moving least square approximations (Eq.( 7)), the resulting system of equations for EFG method with Lagrange multipliers is: .
As can be seen in Eq. ( 10), this method introduces additional unknown variables, , which results in a larger system of equations compared with FEM, which means increasing the computational cost.

EFG with Penalty Method
The Galerkin weak form for EFG with penalty method is: . ( ) Latin American Journal of Solids and Structures 14 (2017) 1017-1039 in which is the penalty parameter used to apply Dirichlet boundary condition.Using MLS approximation leads to the following system of equations: where: .
This method produces a system of equations with the same dimensions of FEM, but choosing the proper value for the penalty parameter is critical for this method.Penalty parameter should be as large as possible, but not too large to produce numerical problems due to increasing the condition number of the system matrix.

EFG Coupled with FEM
In order to couple EFG with FEM using ramp function, Belytschko et al. (1995b), domain of the problem is divided into three subdomains, as shown in Figure 1.In Ω , the solution value at each point is approximated using the MLS approximants as in Eq. ( 7).In Ω , the standard finite element interpolation functions are used to approximate the solution, and in Ω which the interface elements are located, the following expression is used: where and are approximations for u in Ω given by the FE and EFG approximations, respectively, and R(x) is the ramp function defined by Belytschko et al. (1995b): In other words, the ramp function is equal to the sum of the FE shape functions of interface element nodes that are on the EFG boundary.By substituting the FE and EFG solution approximations in Eq. ( 15) it obtains: which the are the shape functions of nodes of interface elements and are defined by [14]: And the derivatives of the shape functions are: The coupled shape functions and their derivatives for a typical one dimensional and two dimensional cases are shown in Figure 2 and Figure 3 respectively.It can be clearly seen that the derivatives of shape functions are discontinues on the borders of interface element.

INTERACTIVE STRUCTURAL ANALYSIS ENVIRONMENT (INSANE)
The INSANE environment is an open source software implemented in Java, an OOP language.It is composed of three main applications: pre-processor, processor and postprocessor.The pre-and postprocessor are interactive graphical applications that provide tools to build representations of the structural problem as well as its results visualization.The processor is the numerical core of the system, responsible for obtaining the results of the analysis.There are different computational methods available for the processor, such as FEM, Fonseca and Pitangueira (2007), Generalized/Extended FEM, Alves et al. (2013), Malekan and Barros (2016), Malekan et al. (2016aMalekan et al. ( , 2016bMalekan et al. ( , 2017)), and Boundary element method, Anacleto et al. (2013), Peixoto et al. (2016).

Main Structure
The numerical core of INSANE is composed of two interfaces: Assembler, and Persistence and two abstract classes: Model and Solution.Figure 4 shows the Unified Modeling Language (UML) form of the numerical core.The Assembler interface builds the matrix of system of equations given by the discretization of a boundary or initial value problem, being implemented with following general representation: Where X is the solution vector; A, B and C are matrices and D is a vector, where may or may not be dependent on the state variable or its derivatives.Solution class has the necessary methods for solving the matrix system of Eq. ( 20), and the Model abstract class contains the data of the discrete model and provides Assembler with all necessary information to assemble Eq. ( 20).
Both Model and Solution communicate with the Persistence interface, which treats the input data and persists the output data to the other applications, any time it observes a modification of the discrete model state.This observation strategy occurs according to the Observer-Observable design Latin American Journal of Solids and Structures 14 (2017) 1017-1039 pattern, Horstmann et al. (2008).In this mechanism when an object of type observer (which implements the interface java.util.observer) is created, it will be added to a list of observers and will have a list of observables (which extends the class java.util.observable).In case of any modification in the state of any observable object, the change propagation mechanism is triggered, and the observer objects are notified to update themselves.This process guarantees the consistency between the observer and the observed components.In INSANE, the observer component is Persistence and the observed components are Solution and Model abstract classes.

Assembler
The MfreeAssembler class has the necessary methods for assembling the matrix system of the Eq. ( 20), with three ways to apply essential boundary conditions: coupling with finite element method, using the method presented in Fernández-Méndez et al. ( 2004), Lagrange multipliers, and penalty method.In static analysis using penalty method or EFG-FE coupling technique, the matrix system of equations is: where the matrix C represents the stiffness matrix of the model, X is the vector of nodal parameters, and D is the vector of forces.The indices u and p indicate unknown and prescribed values, respectively.
In static analysis using Lagrange multipliers the system of equations to be solved is: This class has methods to build different parts of the matrix system, according to Eq. ( 21) and Eq. ( 22).The organization of the classes which implement the Assembler interface is illustrated in Figure 5.In this figure the UML diagram of the MfreeAssembler is showing the main methods of this class.

Model
The objective of the Model abstract class is to represent the discrete model of the problem.It is implemented by the object of class FemModel which itself consists of objects of different classes specifying the attributes of the model.

Solution
The abstract class Solution, manages the solution process of solving matrix system of Eq. ( 20).In case of linear static problem a class named SteadyState implements the class Solution.
Each object of type SteadyState has attributes of two other objects of type Assembler and Line-arEquationSystems .The class LinearEquationSystems has different methods for standard techniques of solving linear equation systems.This class receives the matrix system of equations (Eq.( 20)) built by class Assembler, and computes the solution vector.For solving a linear static meshfree problem this class returns the vector of nodal values.

Analysis Flowchart
Figure 7 shows a detailed flowchart for meshfree analysis described in this section.The process starts with creating a meshfree input file, then defining analysis parameters such as size and shape of the influence domain.Two types of solution procedures are available, either linear elastic or nonlinear analysis.Following selection of the solver type, the problem will be solved and user can work on the results using the available post-processor part.Similar flowchart can be drawn for EFG/FE coupling approach, except the coupling domain have to be defined before solving the model.

NUMERICAL RESULTS
In order to investigate the difference between the methods of applying Dirichlet boundary conditions in EFG method, two well-known cases are studied: Timoshenko beam and Laplace equation.A computational code is developed in order to solve these problems with EFG meshfree method, using the three previously mentioned methods to apply Dirichlet boundary conditions.The case of Timoshenko beam is also implemented in INSANE platform.

Cantilever Beam
For a Timeshenko beam as an elastostatic problem, the differential equation is as follows: . 0 In which is the body force, , and , Γ and Γ are the essential and traction boundaries, respectively.The resulting system of linear equations is: The model of the cantilever beam is illustrated in Figure 8.The following parameters are considered for the model: L=48, D=12, E=30e6, =0.3, and P=-1000.Regular meshes were used, with 4 × 4 quadrature in each EFG cell and also for the coupled FE elements.A linear basis was used for EFG and standard Q4 elements were used in Ω .Two cases are investigated for coupled EFG-FE problem.In the first case, the left half of the beam is modeled using FE elements, and the right half by EFG nodes.In the last strip of the FE domain the interface elements are located.The results for the deformation of the beam, in the y direction are illustrated in Figure 9 and are compared with Lagrange multipliers and analytical solution.The exact solution for the vertical displacement of Timoshenko beam is: where I is the moment of inertia for a beam with rectangular cross-section and unit thickness.
In the second case only a strip of the interface elements is located in the position of Γ and the remaining parts of the beam are modeled by EFG nodes.This strip of interface elements makes it possible to impose the essential boundary condition accurately.
The L2 norm of error for the displacement in the y direction is obtained from the Eq. ( 28) and compared for the two cases, and with the FE results in Figure 10.In this figure h is the distance between the nodes.
As can be seen in this figure, the errors of the coupled method are smaller than that of finite element.Moreover, in case 2, which only a strip of interface elements are inserted at the essential boundary, the coupled method has the best accuracy.The Central Processing Unit (CPU) times of the solution for these two cases are compared with each other and with FE in Figure 11.As can be expected, case 2 involves the largest CPU time (because of the larger number of meshfree particles), but in cases of coarser nodal distributions, it has resulted in better accuracy, using almost the same CPU time.
The results of solving Timoshenko beam problem with INSANE software are illustrated in Figure 12 for vertical displacements of points on the beam's mid-plane.Regular meshes were used, with 5 × 5 quadrature points in each EFG cell and also for the coupled FE elements.A linear basis was used for EFG and standard eight-node quadrilateral (Q8) elements were used in Ω in order to use EFG-FE coupling method.Only a strip of interface elements are placed at Γ and also at x=L, where the load is applied.The results show a very promising accuracy even using a coarse nodal distribution.

Laplace's Equation with Mixed Boundary Conditions
The second numerical example is the Laplace's equation with mixed boundary conditions on a rectangle.The differential equation for this problem is: , 0 , ∈ 0,5 , ∈ 0,5   ( ,0) 0 (0, ) 0 ( ,10) 100 sin 10 (5, ) 0 The analytical solution for this problem is: 100sin( )sinh( ) 10 10 sinh( ) In order to investigate the EFG-FE coupling results in this problem only a strip of interface elements are located at the essential boundary locations, as shown in Figure 13.
The obtained solutions over the whole domain are demonstrated in figures 14-16 for EFG-FE coupled, Lagrange multipliers and penalty method ( 1 4).The absolute values of error are shown for each method in two cases of 200 and 800 nodes.In Figure 14, the values of error are equal to zero at Dirichlet boundaries.Comparing these three figures, penalty method shows the largest error at Dirichlet boundaries.On the other hand, by increasing the number of nodes, the error of approximation at Dirichlet boundaries using Lagrange multipliers shows more significant decrease than the penalty method.

Effect of Penalty Parameter on Convergence
For the case of Timoshenko beam the logarithmic values of energy norm against nodal spacing, h, are plotted in Figure 18.This plot is obtained using the EFG method for different penalty parameters, 1 2 1 8.The nodal spacing was chosen from h = 0.25 to h = 0.1.It is clearly seen that the error in energy norm for penalty parameter 1 4 is less than the other cases, so, this value has been used for obtaining EFG results and been compared with other approaches like Lagrange multipliers and EFG-FE coupling.
Figure 19 compares the rate of convergence considering different values of nodal spacing, h, for the solution of Timoshenko beam using EFG method based on two different approaches in order to apply Dirichlet boundary conditions: Lagrange multipliers method and penalty method.In this figure, the convergence rate of EFG approximation using Lagrange multipliers and penalty method are about 0.94269 and 1.1762, respectively, which shows that EFG solution with penalty method converges faster than EFG with Lagrange multipliers.L2-norm of error against nodal spacing for the 2D Laplace problem is illustrated in Figure 20.This figure contains L2-norm for different approaches, like FE, EFG using Lagrange multipliers, EFG using penalty method (for 1 2 and 1 4) and EFG-FE coupling approach.Although the rate of convergence for FE method, EFG using Lagrange multipliers method and EFG-FE coupled method are the same (equal to 1.99), but the error in L2-norm for EFG using Lagrange multipliers is less than the other two methods.In addition, it can be seen that EFG using penalty method with 1 4 shows the best results considering both convergence rate (equal to 4.449) and the accuracy.

CONCLUSION
Three methods of applying Dirichlet boundary conditions in element-free Galerkin method were investigated: Lagrange multipliers and penalty method as methods based on modification of weak form, and coupling EFG with finite element in order to directly enforce the essential boundary conditions.Numerical examples were studied to investigate different aspects of these methods such as accuracy of the results, convergence, and computational cost.Implementation of aforementioned methods of applying Dirichlet boundary conditions in EFG meshfree method in INSANE platform, an Objectoriented environment, were explained, showing the capability of this code in solving linear problems using EFG method.The coupling between EFG and FE approach works perfectly and it can be used to analyze complex problems with different boundary conditions.This capability of whole implementation in a single code helps the user to migrate easily from FE to EFG analysis, or using both methodologies in a single problem.Also, object-oriented based implementation allow us to use the advantages of other methods, such as the boundary element, hp-cloud, or generalized finite element methods that are currently available in INSANE code.
Latin American Journal of Solids and Structures 14 (2017) 1017-1039 Future studies will focus on solving nonlinear problems such as damage propagation using meshfree part of the INSANE computational platform.In addition, the implemented meshfree program could also communicate with the GFEM module of the INSANE code, and using the implemented enrichment functions and methods of treating crack singularities in those classes, become an extended EFG (XEFG) method and be able to explicitly model the problems with discrete crack propagation in structures; which the last two aspects are also the subjects of ongoing works of the authors.

Figure 4 :
Figure 4: Organization of numerical core of INSANE.

Figure 9 :
Figure 9: Deflection of Timoshenko beam in y direction.

Figure 12 :
Figure 12: Deflection of Timoshenko beam in y direction using INSANE implementation.

Figure 16 :
Figure 16: Solution of Laplace equation using penalty method, a) solution, b) absolute error of approximation using 200 nodes, c) absolute error of approximation using 800 nodes.

Figure 17 :
Figure 17: Values of solution for Laplace equation at x=5.

Figure 18 :
Figure 18: Energy norm against nodal spacing for Timoshenko beam using different penalty parameters, .

Figure 19 :
Figure 19: Rate of convergence in energy norm for the 2-D problem using Lagrange and penalty ( ) methods.

Figure 20 :
Figure 20: Rate of convergence in energy norm for the 2D EFG problem using Lagrange and penalty ( ) methods.