An Energy-based Overset Finite Element Method for Pseudo-static Structural Analysis

This paper addresses a simple energy-based overset finite element method (EbO-FEM) to solve pseudo-static deformation problems consisting of overlapped meshes based on the domain composition method (DCM). This scheme is a non-iterative equation-based method for enforcing the continuity of the displacement field. Hence, the scheme consumes possible minimal computational costs for deformation problems with non-conforming overlapping meshes. The system’s total energy is augmented with continuity constraint energy (CCE) which is a function of the gaps in the displacement field between two overlapping regions. Subsequently, two conventional integration schemes, the Gauss-point projection, and the point-to-point projection, are utilized to discretize the CCE. It is confirmed that both schemes can yield accurate and unique solutions in the overlapped region of the finite element meshes. Further, we proposed a dimensionless relative penalty parameter (DRP). We found that DRP ranging between 1 to 10 is appropriate to robustly obtain accurate solutions for a wide range of scales, stiffness, and geometries, which is supported by three numerical simulations without increasing computational costs after assembling the global matrices and vectors.


Introduction
The finite element method (FEM) has been successfully employed to obtain numerical solutions to problems in the fields of fluid dynamics [1], poromechanics [2,3], solid and structural deformation [4][5][6][7][8], contact mechanics [9][10][11][12], transport phenomenon [13], and several others. FEM is a mesh-based numerical method in which the problem's spatial domain is discretized using finite elements [14]. A collection of such finite elements is called the domain's mesh, and mesh creation is an integral part of finite element analysis. The accuracy of the finite element solution depends on the mesh quality and geometrical errors induced due to the spatial discretization of the computation domain. Thus, when the computation domain contains complex-shaped parts, the analyst must use fine-mesh in these parts to represent the shape as accurately as feasible [6,15,16]. Some approaches have been developed to create a fine mesh with high-quality [17,18]. However, the computation cost of creating such high-quality fine meshes is still a challenge in computational mechanics.
Several researchers have shifted their focus to developing mesh-free methods to avoid the issues of mesh creation and distortion. Some examples of these methods include the material point method [19,20], the smooth-particle-hydrodynamics [21,22], and the particle FEM [23]. Although mesh-free methods solve the problem of mesh creation [24], they face new challenges due to their inefficiency in representing the boundary of the computation domain and material interfaces. Some of these challenges are prescribing boundary conditions and interface conditions. These issues can be handled easily and effectively using FEM through essential and natural boundary conditions. Therefore, efforts have been made to improve FEMs to reduce the need for such highquality fine meshes in the complex shape of the computation domain. In this regard, overset mesh schemes and the automatic meshing with overlapping and regular elements (AMORE) framework were proposed as an innovative method by Bathe et al. [25][26][27][28]. The method allows separately built meshes to overlap [25]. Applying a consistent test function to the overlapping region makes it possible to obtain a solution that satisfies the governing equations over the overlapping and non-overlapping regions [26]. The AMORE framework has been validated for convergence, accuracy, and applicability. This framework has high accuracy and applicability when the geometric data are known in advance, such as machine parts that are artificially built. In addition, nonconforming meshes should be rearranged to be conformed. For fluid analysis and multi-physical simulations, the Domain Composition Method (DCM) has been developed as a parallelization-friendly method [29][30][31][32][33][34] the method is applicable for both conforming/non-conforming meshes. Through these investigations, it is shown that the DCM reduces the computational costs [29], typically memory consumption.
Further, when FE-analysis is applied to civil/rural engineering [35][36][37][38][39] or plant biomechanics [40][41][42][43][44], the cost of handling meshes in the overlap domain is an inevitable problem. For example, even a simple irrigation system (refer to Fig. 1) or RC frame bridges usually contain tiny but physically essential parts such as joints, shoes, piers and parapets. Moreover, as depicted in Fig. 1, a typical matured soybean plant can have more than 200 joints between organs such as leaves, stems, and fruits [40,41,[45][46][47]. Thanks to the advance in imaging technologies, it is currently possible to computationally create partially overlapped objects by combining 3-D scan data [40,48], reconstructed CAD data [38,39], and algorithm-generated 3-D point cloud or mesh data [45,47,49]. However, these computational schemes provide a large number of overlapping sections and increase the number of iterations or the number of unknowns to obtain the uniform solution based on the conventional DCM.
These problems for civil/rural engineering invariably necessitate two expensive computation costs: (1) the generation of nonconforming meshes for each part and (2) the conversion of nonconforming meshes to conforming mesh by nodal point rearrangement. Besides, the number of areas where overlaps occur increases drastically. Therefore, there is a need to develop a simple computational approach for enforcing continuity of the displacement field in the overlap regions for nonconforming meshes. Especially, it is necessary to minimize the extra computational costs to prevent unnatural gaps among overlapping meshes. It should be avoided that (1) iteration or trial-and-error procedure to enforce the gap-minimizing constraints [50] and (2) dynamic change in the size of the unknowns; the first one appears in Fig. 1 Target problem of energy-based overset mesh method. EbO-FEM targets pseudo-static deformation analysis for complex structures consisting of multiple members (left), where a 3D model is created for each member (centre), and numerical analysis is performed by overlaying the models (right) the conventional equation-based methods [29], and the second procedure is required in the Lagrange multiplier methods.
This paper aims to describe a non-iterative, equation-based overset mesh scheme for pseudo-static deformation problems. Our key ideas are (1) to design a continuity constraint energy (CCE) function and (2) to find a solution to minimize the sum of the total energy of the system considering CCE and the conventional elastic potential energy, which are explained in Sect. 2. The proposed approach does not require iterative procedures such as Newton's loop nor extra terms such as the Lagrange multipliers. We also introduced dimensionless relative penalty parameters (DRP) to provide a concrete form of the CCE, and numerically find the optimal range of the DRP value through three numerical simulations with various scales, stiffness and geometries as shown in Sects. 3 and 4. Further, the uniqueness and the existance of the solution is stated in Sect. 3.

Energy-based Overset Mesh Method
In the EbO-FEM, the analyst manually creates the system's computation domain as a collection of finite element meshes of different subdomains ( n , n = 1, 2, · · · ). In the present method, the meshes can overlap, unlike in the classical FEM. Figure 1 shows that such a scheme is desirable when the computation domain in a target control volume comprises multiple solid bodies or materials with complex shapes.
Let us now consider two computation domains, denoted by p and s , where the superscripts p and s represent the primary and secondary domains, respectively. The primary and secondary domains are allowed to overlap, and the overlapped domain, designated by ps , is given by the following: Thus, any point x i ∈ ps also belongs to the spatial domains p and s . Let us now denote the displacement field in p and s by u p := u(x), x ∈ p , and u s := u(x), x ∈ s , respectively. Similarly, the solution field in ps is represented by u ps . The strong condition of displacement consistency states that at any point x ∈ ps , u ps = u s = u p . The requirement is expressed by introducing an energy function g : where, g denotes the jump discontinuity in the displacement field as shown in Eq. (2).
In Eq. (1), for simplicity, we have used: where, is a penalty parameter. Finally, the problem statement is to minimizē where¯ denotes the summation,¯ = ¯ d , of energy function for each domain,¯ d , which can be obtained from the governing equations.
The statement is expressed by the variational principle.
Using Eqs. (1) and (2), the variation δ in Eq. (5) is defined as a Gateaux derivative with respect to space x: where δU i (x) is a variation of the displacement field. In the following subsections, we implement this formulation by discretizing Eq. (6) using the FEM. The first term can be discretized using conventional FEM formulations [51]; however, the discretization method for the second term must be carefully designed to enforce continuity condition as depicted in Fig. 1. We investigate the accuracy and characteristics of two of the most simple methods, Gauss-point and point-to-point projections.

Remark 1 Conventional DCM proposes an equation-based constraint
where δu is a test function, u 1 and u 2 are the unknown defined with respect to the subdomains 1 and 2 , respectively. There are three possible strategies for implementing the constraint condition for the pseudo-static deformation problems: (1) iterative approach such as Dirichlet-Dirichlet or Dirichlet-Neumann approaches [29], (2) Lagrange multiplier method, and (3) penalty method. (1) requires iteration such as Newton's method and (2) needs additional unknowns, Lagrange multipliers, which increases computational costs as the number of overlapping elements increases. The motivation of this study is to provide a non-iterative, computationally inexpensive approach. From this standpoint, we take (3) approach and provide (a) a simple and linear formulation to enforce continuity in the displacement field which does not requires iterative algorithm to obtain the solution and (b) a procedure to determine the optimal penalty parameter.

Variational Equality
We use the EbO-FEM to solve a pseudo-static deformation problem based on FEM, where the energy function, , is regarded as the elastic/kinetic energy of the system, and the energy function, g , represents free energy induced by the continuity constraint. The condition of consistency of the solution field state that at any x ∈ ps , the gap function, denoted here as g p→s , and given by the following: should be minimized. The aforementioned definition reveals that g p→s is the measure of the solution field's jump discontinuity in the overlapping region, and is given as follows: In this study, the displacement consistency condition is enforced by employing the penalty method. Thus, the energy term associated with the gap function (denoted here as p→s g ) is given as follows: where ε is a penalty parameter and denotes the Euclidean norm of . Therefore, the displacement consistency condition can be given by the following energy minimization statement [52]: δ Consequently, the variational equality (11) can be expressed as follows: where (·) denotes the inner product operation of two vectors.
Similarly, the variational equality in the opposite direction is given by the following: with,

Consistent Discretization Based on Finite Element Method
Let N p A , A = 1, 2, · · · n p and N s A , A = 1, 2, · · · n s be the shape functions of the solution field in p and s , respectively (Fig. 3). Here, n p and n s are the total number of spatial nodes in the computation mesh of s and p , respectively. Then, the nodal interpolation of u s and u p can be written by following expressions: and where u s I denotes descritized displacement field. By substituting Eqs. (16) and (17) into Eqs. (8) and (13), the discretized jump functions and its variation are obtained as follows: The variational equality (Eq. (12)) is discretized as where J is a determinant of the Jacobian matrix and w is a weight for the Gauss-Legendre integral.
Although the constraint should be given for all points, x i ∈ ps , it is computationally inefficient to do so throughout the intersection. Therefore, consider a set of representative points,ˆ ps ⊂ ps , and weakly enforce the constraint condition inˆ ps rather than ps .
There are various options forˆ ps ; for instance, the set of representative points can be a set of nodes from the mesh in s . In this case, Eq. (19) is rewritten as follows: whereδ p→s g is the variational equality defined inˆ ps . Finally, a matrix form of the energy variation is given as follows: where where T indicates a transpose operator for a matrix.

Discretization by Gauss-point Projection
In this subsection, we propose Gauss-point projection as a method for computing the numerical integrals in Eq. (1) based on the Gauss-Legendre quadrature points. Figure 3 depicts the schematic of the Gauss-point projection. We now consider the scenario where two finite elements denoted p and s undergo overlap. Here, we begin by locating all Gauss-points of finite element p that are included within finite element s. Next, the position vector of the identified effective Gauss-points in the entire coordinate system is determined. Then, the equivalent local position vector for the global position vector of the element s is determined. Eq. (1) is evaluated for each element based on the local position vector and shape functions for those p and s by the following equation under three-dimensional space (a = 1, 2, 3 and b = 1, 2, 3): where x p a and x s a are the nodal positions for the primary and the secondary elements, respectively. G P e is a gap-energy function for an element, ngp denotes the number of Gauss-points, N s a (ḡ) and N p b (g i ) are shape functions for each hexahedron element, u s a and x p a are nodal value of unknowns,¸is a local coordinate, and x is a global coordinate. It is worth noting that the left hand side term is equal to zero in case that no projection from the primary domain to the secondary domain are found. Here,¯ and share the same Gaussian quadrature for simplicity. In addition, the accuracy can be improved by increasing the number of integration points, making the method simple and scalable. Section 4 contains the concrete formulations for one-dimensional and three-dimensional problems.

Discretization by Point-to-Point Projection
Although the Gauss-point projection introduced in the previous subsection can accurately discretize the numerical integration term, evaluating the energy at each Gauss-point is computationally expensive. The search cost of effective Gauss-points is up to N M when the number of elements is N and the number of Gauss-points per element is M, which is also similar to the Gauss-point projection or closest point projection in terms of contact mechanics [10,[52][53][54][55] and DCM [29]. In addition, it takes up to N times to search for a pair of elements for each effective Gauss-point and perform the projection, requiring a total of M N 2 orders of computation times. To reduce this computational cost, we also propose a simpler algorithm, point-to-point projection, that employs the sum of the energies for the sampled points rather than evaluating the integral term. We investigate the accuracy and applicability by numerical experiments in the following sections. A concrete formulation for the closest-point projection is shown in the Eqs. (21) to (23).

Dimension of the Consistent Penalty Parameter
It should be considered that the penalty parameter has different unit between the point-topoint projection and the Gauss-point projection. This subsection discusses the consistent unit for the penalty parameters.
An increment of internal energy of the thermodynamical system is given from the first law of the thermodynamics where U is the internal energy of a closed system, Q is the quantity of the incoming energy as heat and W is an amount of the thermodynamical work done by the system and has a dimension of k N · m, respectively. W is rewritten using Cauchy's stress tensor and the strain tensor Considering the extra energy term of Eq. (3), we obtain therefore, Eq. (26) is rewritten as Energy-minimizing problem under the adiabatic process is given by The unit of the first term is k N ·m, hence, the unit of the second term should be the same. In the case of the point-to-point projection approach under three dimensional Cartesian coordinate, the unit of the penalty parameter should be k N /m and the unit for the Gauss-point projection approach is k N /m 4 .  Figure 4 depicts the visual algorithm. First, the elemental stiffness matrix is created as in (1).

Algorithm
There are no constraints in the overset region at this stage, implying that the coefficient matrix is orthogonal to each subregion. The next step calculates the overset stiffness matrix using Gauss-point or point-to-point projections and adds it to the global coefficient matrix. Here, the overset stiffness matrix is introduced between stiffness matrices for adjacent subdomains as a patch shown in (2), and its location can be estimated in advance using geometric correlations. Finally, the solution can be obtained immediately by introducing the Dirichlet condition and solving the system of linear equations (3).

Uniqueness and Existance of the Solution
The Lax-Milgram's theorem has been widely used to confirm the uniqueness and the existence of the solution for constraint deformation problems [6,56,57]. We also utilize Lax-Milgram's theorem for the proof. The problem is written in terms of the bilinear form, where and , is an inner product operator.

) is a bilinear map on the real linear vector space V
Proof In case that the constitutive equation is linear, for arbitrary vectors u, v, w, in V, Similarly, Furthermore, for arbitrary real scalar value k, one holds

Lemma 3.2 For Hilbert space H (V, · ),â is bounded.
Proof Considering arbitrary vectors (u, v) ∈ V, it holds that is shown by the Cauchy-Schwarz inequality, where C 1 and C 2 are positive real values since is a positive definite.

Lemma 3.3 For Hilbert space H (V, · ), B is bounded and coercive.
Proof Considering arbitrary vectors (u, v) ∈ V, it holds that is shown by the Cauchy-Schwarz inequality, where C 1 and C 2 are positive real values. Also, consideringâ(u, u) ≥ 0, the Korn's inequality shows that where D c is a Korn's constant shown in [6], D is a positive scalar value.
After the above-mentioned mathematical preparations, the existence and the uniqueness of the weak form is proofed as shown below.
Proof We consider the weak form where f (v) is a bounded linear functional [6]. Using the Lemmas 3.1 to 3.3 and the Lax-Milgram's theorem, Eq. (35) has a unique solution u for arbitrary v

Numerical Examples
In the previous section, we described our proposed EbO-FEM and presented its two possible formulations: Gauss-point and point-to-point projection. The energy-based overset mesh method requires that the overset domains of two discrete domains satisfy C 0 continuity. When two physical quantities, in this case, displacements, are given for the same position and a difference exists, energy exists and minimizes the difference.
Specific numerical examples are used to illustrate the possible basic considerations. We apply the overset mesh schemes to deformation problems. Here, the displacement field is unknown. The first and most significant question is how to determine the size of the penalty parameter a priori. The proposed method will fail if the penalty parameter is zero. By contrast, the proposed method should succeed if the penalty parameter is large, but it is unclear how large the penalty parameter should be. If the penalty parameter is extremely large, the condition number will increase and the problem will become ill-conditioned, which results in a drastic increase in the computation time. We employed two projection frameworks to analyze one-dimensional (1-D) and three-dimensional (3-D) problems for addressing these questions.

Example 1: 1-D Tensile Problem of an Elastic Cantilever
The simplest problem to solve is the tensile problems for a 1-D cantilever. As shown in Fig.  5, let us consider pulling a bar consisting of two overset meshes, where the left side of s is fixed on the left wall, and displacement,ū, is loaded on the right side of p . The case is formulated using FEM and the current method. The discretized equations are provided in Appendix.
Young's modulus is set to 1.0 for the two components for simplicity, and the penalty parameter is varied from 0.00001 to 10.0 to examine the solution's response to each penalty parameter. It is expected that given a larger penalty parameter, the jump-minimizing constraint will be satisfied more strictly. However, the previous study for contact analysis such as Zavarise and De Lorenzis (2009) [58] shows that the penalty parameter cannot be too The length of total bar was initially 3.5L e but became 3.5L e +ū, then the analytical solution on the 6 points x i (i = 1, 6) The numerical solutions obtained by the proposed methods are compared with the analytical solution to verify the accuracy and properties of the proposed method. L e = 2 is used for simplicity. Figure 6 shows the exact solution and calculation results for the displacement at each nodal position. The figure shows that as the DRP increases from zero to 1.0, the results becomes closer to the exact solution. The solution was more accurate in cases of the Gausspoint projection than the point-to-point. The figure, however, implies that only C 0 -continuity is almost satisfied in point-to-point projection, and nothing above C 1 is guaranteed. The phenomena has been also reported [29] and caused by the Dirichlet-Dirichlet type constraint. Alternatively, Gauss-point projection approximates both C 2 and C 0 continuity around the projection point. The overset mesh scheme was initially intended to satisfy only the C 0 continuity, and both solutions do so; however, the Gauss-point projection also satisfies higherorder continuity near the projection.
Overall, the current scheme successfully minimizes the gap as intended. As expected, increasing the penalty parameter yields a solution that is close to the exact solution. The smallest feasible penalty parameter is 1.0 and the results are unaffected if the penalty parameters are more than 10.0 at the projection point.
Notably, the penalty parameters in point-to-point and Gauss-point projections have different implications and dimensions. For 3-D cases, the former has a dimension of N /m and the latter has a dimension of N /m 4 as discussed in the previous section. This difference results in the difference in the dimension of the dimensionless relative penalty-parameters (DRP) E by Notably, the penalty parameters in point-to-point and Gauss-point projections have different implications and dimensions. For 3-D cases, the former has a unit of N /m while the latter has N /m 4 as discussed in the previous section.  Figure 7 shows the initial and boundary conditions of the problem. The cantilever (Young's modulus = 1000.0, Poisson ratio = 0.0) is divided into three discrete cantilevers for EbO-FEM:

Example 2: 3-D Bending Problem of an Elastic Cantilever
The left cantilever is fixed in the left rigid wall, the right one is loaded displacement, and the middle one is traction-free. In the conventional one-domain FEM formulation, the problem has a solution where the three are connected into a uniform cantilever. Each subdomain in the overlapping mesh has 30 × 9 × 9 = 2730 elements and the mesh for the one-domain FEM has 90 × 9 × 9 = 2730 elements. Each mesh overlaps 0.1 m each other for the overlapping mesh.  Figure 8 shows the deformed mesh and the effect of the penalty parameter. The figure shows that the results of the point-to-point and Gauss-point projections are identical to the analytical and FEM solutions. In terms of the penalty parameter, a relative penalty parameter of 1.0 appears to be sufficient, implying that the penalty parameter can be the same as Young's modulus. Figure 9 shows the relationship between the force-curve and displacement. The DRP here is 10.0, the length of the bar is 5.8 m, the width is changed from 0.3 m to 0.7 m for each case. The current method's solution is quantitatively almost the same as the analytical solution and the one-domain FEM solutions in any width. Figure 10 plots the relationship between DRP and the number of iterations required for the convergence at the BiCGSTAB algorithm where the tolerance of the energy norm is 1.0 × 10 −10 . For both cases, DRP < 0.1 enforces too weak constraint, which results in inaccurate solutions as visible in Fig. 9. By contrast, DRP > 1000.0 causes severe illcondition in the case of Gauss-point projection, as well as DRP > 10.0 in the point-to-point projection. Hence, it is difficult for the point-to-point projection to obtain accurate and robust solutions. The Gauss-point projection gives accurate and robust solutions in case the DRP is in between 1 to 1000, at least in these initial/boundary conditions. The results also show that the Gauss-point projection cases only require 0.1 to 2.0 times larger number of iterations, which are sufficiently small computational costs as an overset mesh scheme compared to that of the one-domain case. Figure 11 shows how the condition number of the coefficient matrix varies with a change in relative penalty parameters. Since it is necessary to use a dense matrix to compute a condition number, we use a more coarse mesh (30 ×4 times 4 elements for each subdomain) due to the limitation of the computational resources. The result indicates that the condition numbers are minimized when the relative penalty parameter is about 1.0 in both cases. Alternatively, when the relative penalty parameter is smaller than 0.01 or larger than 10.0, the condition

Example 3: 3-D Deformation Analysis of Self-Weighted RC Frame Bridge
Investigating the applicability of the method to practical problems, we simulate here a fourspanned reinforced concrete (RC) frame bridge under a self-weighted condition (Fig. 12). Although similar simulations can be conducted by using the conventional FEM, it is a huge efforts to create such a qualified conforming mesh. By contrast, it is cost-efficient to use  Figure 12 shows the initial and boundary condition of the computational domain. The RC frame bridge is a typical one for a railway bridge in Japan, with 40 m in length, 10 m in width and 7 m in height. The Young modulus, the shear modulus, and the density are 2.5 × 10 7 kPa, 1.0 × 10 7 kPa, and 25.0 kN/m 3 , respectively. The bottom of the piers is fixed on the ground for simplicity.
The mesh consists of 74,201 points and 48,832 hexahedral elements for FEM. The EbO-FEM case uses a non-conforming mesh consisting of three parts as visible in Fig. 13, where the centre part consists of 71,119 points with 50368 elements, the set of floor parts has 4,848 points/3,000 elements, and the pair of the edge parts has 3,333 points/2,000 elements. Figure 14 shows the deformed mesh with σ zz components of the Cauchy stress. As can be seen, the solution of EbO-FEM is consistent with the FEM solution in cases where the dimensionless relative penalty parameter (DRP) is greater than 1.0. On the other hand, the EbO-FEM solution is inconsistent in cases where the parameter is smaller than 1.0. The phenomena are also observed in the previous sections for 1-D and 3-D problems, which again suggests that (1) the DRP can be a universal measure for EbO-FEM scheme and (2) DRP should be greater than 1.0 to obtain consistent solutions.

Conclusion
In this study, an EbO-FEM, based on FEM, is proposed for nonconforming meshes to solve the multibody mesh overlap problem. To apply the continuity condition, we introduce a continuity constraint energy, which is analogous to the contact constraint energy in the nodeto-segment method [10,54,58] based on the equation-based DCM [29]. It has been confirmed that the method can produce almost continuous and unique solutions on separately defined overlapping finite element meshes. We utilized two projection methods, point-to-point projection and the Gauss-point projection, and showed that the Gauss-point projection is a better strategy because it is more robust and requires less computation time for solving linear equations based on BiCGSTAB procedure. Also, we proposed a dimensionless measure, DRP, to evaluate the magnitude of the constraint, and showed that the DRP of 0.1 to 10.0 is the best from the standpoints of accuracy, robustness, and computation costs. Surprisingly, the present scheme does not increase the computation costs after assembling the global matrix, which is a great benefit for problems with many overlapping regions as these 3-D complex problems.
The proposed non-iterative, energy-based overset scheme helps solve pseudo-static deformation problems with multi-overlapping meshed appearing in civil engineering structures on vegetated ground [59,60] and plant growth [40,44,61]. In addition, this EbO-FEM can be combined with discretization methods other than FEM as long as the continuous constraint energy is well-defined among nonconforming meshes. Therefore, this approach can be extended for dynamic elastic problems or wave propagation problems. Future studies will apply the method to dynamic analysis to see if it can provide an accurate solution for elastodynamic problems.

Acknowledgements This work was supported by JSPS KAKENHI Grant Numbers 17J02383 and 20K22599.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Appendix
Governing equation is where E is Young modulus, u is displacement and x is the coordinate. The weak form is derived based on Galerkin method where w is a weight function, x r and x l are the coordinate of the left and right sides of the bar, respectively. Here, no traction is loaded in the boundary, then du dx = 0 in the boundary, hence, Introducing linear shape-functions, w and u are written as following expression The ξ is a local coordinate defined in [0 ≤ ξ ≤ 1], and the subscription r and l indicates the left and right nodes of a line element for unknowns w and u, respectively. Then the weak form is discretized as shown below The dx dξ is the length of the element L e , consequently, the discretized equations are given.
The global equation is assembled as Eq.

A.1 Formulation Based on Gauss-Point Projection
From Eq. (20), the constraint condition is active in ps , then, the system of equation based on Eqs. (13) and (15) is given δ gp = L e u 4 u 5 u 2 u 3 ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ (A.10) where K gp AB is a globally assembled matrix ofK gp AB .

A.2 Formulation Based on Point-to-Point Projection
From Eq. (20), the constraint condition is active in ps , then, the system of equation based on Eqs. (13) and (15) Finally, the boundary condition u 6 =ū and u 1 = 0 is introduced. ⎛ Solving this set of equations yields the results shown in Fig. 6.