Analysis of patchtest consistent atomistictocontinuum coupling with higherorder finite elements
Abstract.
We formulate a patch test consistent atomistictocontinuum coupling (a/c) scheme that employs a secondorder (potentially higherorder) finite element method in the material bulk. We prove a sharp error estimate in the energynorm, which demonstrates that this scheme is (quasi)optimal amongst energybased sharpinterface a/c schemes that employ the Cauchy–Born continuum model. Our analysis also shows that employing a higherorder continuum discretisation does not yield qualitative improvements to the rate of convergence.
Key words and phrases:
atomistic models, coarse graining, atomistictocontinuum coupling, quasicontinuum method, error analysis1. Introduction
Atomistictocontinuum (a/c) coupling is a class of coarsegraining methods for efficient atomistic simulations of systems that couple localised atomistic effects described by molecular mechanics with longrange elastic effects described by continuum models using the finiteelement method. We refer to [6], and references therein for an extensive introduction and references.
The presented work explores the feasibility and effectiveness of introducing higherorder finite element methods in the a/c framework, specifically for quasinonlocal (QNL) type methods.
The QNLtype coupling, first introduced in [15], is an a/c method that uses a “geometric consistency condition” [3] to construct the coupling between the atomistic and continuum model. The first explicit construction of such a scheme for twodimensional domains with corners is described in [12] for a neareastneighbour manybody site potential. We call this construction ”G23” for future reference. This approach satisfies force and energy patch tests (often simply called consistency), which in particular imply absence of ghost forces.
We will supply the G23 scheme with finite element methods of different orders and investigate the rates of convergence for the resulting schemes. Our conclusion will be that secondorder finite element schemes are theoretically superior to firstorder schemes, while schemes of third and higher order do not improve the rate of convergence. This is due to the fact that the consistency error of the a/c scheme is dominated by the modelling error committed at the a/c interface. We will also explore, for some basic model problems, how well secondorder schemes fare in practise against firstorder schemes.
1.1. Outline
The theory of highorder finite element methods (FEM) in partial differential equations, and applications in solid mechanics is well established; see [14] and references therein. However, most work on the rigorous error analysis of a/c coupling has been restricted to P1 finite element methods; the only exception we are aware of is [11], which focuses on blendingtype methods.
In the present work we estimate the accuracy of a QNL method employing a P2 FEM in the continuum region against an exact solution obtain from a fully atomistic model. Since stability of QNL type couplings is a subtle issue [9] we will primarily analyse the consistency errors, taking account the relative sizes of the fully resolved atomistic region and of the entire computational domain (Sections 5.15.4). We will then optimize these relative sizes as well as the mesh grading in the continuum region in order to minimize the total consistency error (Section 5.5). We will observe that, using P1FEM in the continuum region, the error resulting from FEM approximations is the dominating contributor of the consistency estimates, which implies that increasing the order of the FEM can indeed improve the accuracy of the simulation. We will show that, using PkFEM with , the FEM approximation error is dominated by the interface error which comes purely from the G23 construction, and in particular demonstrate that the P2FEM is sufficient to achieve the optimal convergence rate for the consistency error. Finally, assuming the stability of G23 coupling (see [9] why this must be an assumption), we prove a rigorous error estimate in §6.
Finally, we conduct numerical experiments on a 2D antiplane model problem to test our analytical predictions. The numerical results display the predicted error convergence rates for the fully atomistic model, P1FEM G23 model, and P2FEM G23 model.
2. Preliminaries
Our setup and notation follows [12]. As our model geometry we consider an infinite 2D triangular lattice,
We define the six nearestneighbour lattice directions by , and , where denotes the rotation through the angle . We supply with an atomistic triangulation, as shown in Figure 1, which will be convenient in both analysis and numerical simulations. We denote this triangulation by and its elements by . We also denote , and , for .
We identify a discrete displacement map , , with its continuous piecewise affine interpolant, with weak derivative , which is also the pointwise derivative on each element . For , the spaces of displacements are defined as
We equip with the seminorm, . From [8] we know that is dense in in the sense that, if , then there exist such that strongly in .
A homogeneous displacement is a map , where .
For a map , we define the finite difference operator
(2.1)  
Note that .
2.1. 2D manybody nearest neighbour interactions
We assume that the atomistic interaction is described by a nearestneighbour manybody site energy potential ,, with . Furthermore, we assume that satisfies the point symmetry
The energy of a displacement , given by
is welldefined since the infinite sum becomes finite. To formulate a variational problem in the energy space , we need the following lemma to extend to .
is continuous and has a unique continuous extension to , which we still denote by . Furthermore, the extended is times continuously Fréchet differentiable.
Proof.
See Lemma 2.1 in [4]. ∎
We model a point defect by adding an external potential with for all , where is some given radius (the defect core radius), and for all constants . For example, we can think of modelling a substitutional impurity. See also [5, 7] for similar approaches.
We then seek the solution to
(2.2) 
For we define the first and second variations of by
We use analogous definitions for all energy functionals introduced in later sections.
2.2. The Cauchy–Born Approximation
The Cauchy–Born strain energy function, corresponding to the interatomic potential is
where is the volume of a unit cell of the lattice . Thus is the energy per volume of the homogeneous lattice .
2.3. The G23 coupling method
Let denote the set of all lattices sites for which we want to maintain full atomistic accuracy. We denote the set of interface lattice sites by
and we denote the remaining lattice sites by . Let be the Voronoi cell associated with site . We define the continuum region ; see Figure 2.
A general form for the GRACtype a/c coupling energy [3, 12] is
where . The parameters are to be determined in order for the coupling scheme to satisfy the “patch tests”:
is locally energy consistent if, for all ,
(2.3) 
is force consistent if, for all ,
(2.4) 
For the sake of brevity of notation we will often write . Following [12] we make the following standing assumption (see Figure 3 for examples).
(A0) Each vertex has exactly two neighbours in , and at least one neighbour in .
Under this assumption, the geometry reconstruction operator is then defined by
see Figure 4. The resulting coupling method is called G23 and the corresponding energy functional . This choice of coefficients (and only this choice) leads to patch test consistency (2.3) and (2.4).
For future reference we decompose the canonical triangulation as follows:
(2.5)  
2.4. Notation for a P finite element scheme
In the atomistic and interface regions, the interactions are represented by discrete displacement maps, which are identified with their linear interpolant. Here, we identify the displacement map with its P1 interpolant. No approximation error is committed.
On the other hand, in the continuum region where the interactions are approximated by the Cauchy–Born energy, we could increase the accuracy by using PFEM with . In later sections we will review that the Cauchy–Born approximation yields a 2ndorder error, whereas employing the P1FEM in the continuum region would reduce the accuracy to first order. In fact, we will show in that, with optimized mesh grading, P2FEM is sufficient to obtain a convergence rate that cannot be improved by other choices of continuum discretisations. Highorder PFEM with will increase the computational costs but yield the same error convergence rate (see § 3.5).
Let denote the inner radius of the atomistic region,
where denotes the ball of radius centred at . In order for the defect to be contained in the atomistic region we assume throughout that .
Let denote the entire computational domain and denote the inner radius of , i.e.,
Let be a finite element triangulation of which satisfies, for ,
In other words, and coincide in the atomistic and interface regions, whereas in the continuum region the mesh size may increase towards the domain boundary. The optimal rate at which the mesh size increases will be determined in later sections.
We note that the concrete construction of will be based on the choice of the domain parameters and ; hence, when emphasizing this dependence, we will write . We assume throughout that the family is uniformly shaperegular, i.e., there exists such that,
(2.6) 
This assumption eliminates the possibility of extreme angles on elements, which would deteriorate the constants in finite element interpolation error estimates. For the most part we will again drop the parameters from the notation by writing but implicitly will always keep the dependence.
Similar to (2.5), we define the atomistic, interface and continuum elements as and , respectively. Note that and . We also let denote the number of degrees of freedom of .
We define the finite element space of admissible displacements as
(2.7)  
In defining we have made two approximations to the class of admissible displacements: (1) truncation to a finite computational domain and (2) finite element coarsegraining.
The computational scheme is to find
(2.8) 
is embedded in via point evaluation. Through this identification, is welldefined for all .
We will make this identification only when we evaluate . The reason for this is a conflict when interpreting elements as lattice functions is that we identify lattice functions with their continuous interpolants with respect to the canonical triangulation , which would be different from the function itself. However, for the evaluation of this issue does not arise.
3. Summary of results
3.1. Regularity of
The approximation error analysis in later sections requires estimates on the decay of the elastic fields away from the defect core. These results follow from a natural stability assumption:
(A1) The atomistic solution is strongly stable, that is, there exists ,
(3.1) 
where is a solution to (2.2).
Suppose that (A1) is satisfied, then there exists a constant such that, for ,
Proof.
See Theorem 2.3 in [4]. ∎
3.2. Stability
In [9] it is shown that there is a “universal” instability in 2D interfaces for QNLtype a/c couplings: it is impossible to prove in full generality that is a positive definite operator, even if we assume (3.1). Indeed, this potential instability is universal to a wide class of generalized geometric reconstruction methods. However, it is rarely observed in practice. To circumvent this difficulty, we make the following standing assumption:
(A2) The homogeneous lattice is strongly stable under the G23 approximation, that is, there exists which is independent of such that, for sufficiently large,
(3.2) 
Since (3.2) does not depend on the solution it can be tested numerically. But a precise understanding under which conditions (3.2) is satisfied is still missing. In [9] a method of stabilizing 2D QNLtype schemes with flat interfaces is introduced, which could replace this assumption, however we are not yet able to extend this stabilizing method for interfaces with corners, such as the configurations discussed in this paper.
3.3. Main results
To state the main results it is convenient to employ a smooth interpolant to measure the regularity of lattice functions. In Lemma 6.1, we define such an interpolant for , for which there exists a universal constant such that, for all , ,
where .
Consistency error estimate
In (5.6) we define a quasibest approximation operator , which truncates an atomistic displacement to enforce the homogeneous Dirichlet boundary condition, and then interpolates it onto the finite element mesh.
Our main result is the following consistency error estimate.
If is a solution to (2.2) then we have, for all ,
(3.3)  
where corresponds to the continuum region of , and with .
Optimizing the approximation parameters
Before we estimate the error , we optimize the approximation parameters in the computational scheme. This means that the radius of the atomistic region, the radius of the entire computational domain and the mesh size should satisfy certain balancing relations. We only outline the result of this optimisation and refer to § 6.7 for the details.
Due to the decay estimates on the dominating terms in (3.3) turn out to be
(3.4) 
(We will see momentarily that the mesh size plays a minor role.) These two terms result from the nature of the coupling scheme and the farfield truncation error. In particular, both of these cannot be improved by the choice of discretisation of the Cauchy–Born model, e.g., order of the FEM. We also note that, if we had employed a P1FEM, then the limiting factor would have been .
We can balance the two terms in (3.4) by choosing . It then remains to determine a meshsize so that the finite element error contribution,
remains small in comparison. We show that the scaling is a suitable choice, with , under which both terms become of order .
Thus, we have determined the approximation parameters in terms of a single parameter . The quasioptimal relations for P2FEM discretisatino of the Cauchy–Born model are summarised in Table 1.
consistency error  

P2FEM  
P1FEM 
Error estimate
To complete our summary of results, we now use the Inverse Function Theorem to obtain error estimates for the strains and the energy.
3.4. Setup of the numerical tests
For our numerical tests, we consider an antiplane displacement . We choose a hexagonal atomistic region with side length and one layer of atomistic sites outside as the interface. To construct the finite element mesh, we add hexagonal layers of elements such that, for each layer , , with ; see Figure 5. The procedure is terminated once the radius of the domain exceeds . This construction guarantees the quasioptimal approximation parameter balance to optimise the P2FEM error. The derivation is given in Section 6.7.
In our tests we compare the P2G23 method against

a pure atomistic model with clamped boundary condition: the construction of the domain is as in the P2G23 method, but without continuum region;

a P1G23 method: the construction is again identical to that of the P2G23 method, but the P2FEM in the definition of is replaced by a P1FEM. The same mesh scaling as for P2 is used (see also [3] where this is shown to be quasioptimal).
The site potential is given by a nearestneighbour embedded atom toy model,
with and . This is the antiplane toy model as the one used in [4].
The external potential is defined by , which can be thought of as an elastic dipole. A steepest descent method, preconditioned with a finite element Laplacian and fixed (manually tuned) stepsize, is used to find a minimizer of , using as the starting guess.
In order to compare the errors, we use a comparison solution with atomistic region size and other computational parameters scaled as above.
The numerical results, with brief discussions, are shown in Figures 6–9. The two most important observations are the following:

the numerical tests confirm the analytical predictions for the energynorm error, but the experimental rates for the energy error are better than the analytical rates. Similar observations were also made in [4].

With our specific setup, the improvement of the P2GR23 over P1GR23 is clearly observed when plotting the error against , but when plotted against the improvement is only seen in the asymptotic regime. This indicates that further work is required, such as a posteriori adaption, to optimise the P2GR23 in the preasymptotic regime as well.
3.5. Extension to highorder FEM
If we apply higherorder FEM in the continuum region, then to extend our error analysis we would need a smooth interpolant of with higher regularity than . A suitable extension given in [5] is, for arbitrary , a piecewise polynomial of degree with properties analogous to those stated in Lemma 6.1. The resulting higherorder decay rate indicates that the use of highorder FEM could be beneficial.
However, as we have pointed out in § 3.3.2, if we employ the mesh grading with in the continuum region, the total approximation error cannot be improved by using PFEM with , since the dominating term is the interface error for , which results from the construction of G23 coupling and is not affected by the choice of FEM.
If we consider a coarser mesh for highorder FEM in hopes of reducing the number of degrees of freedom, i.e., choosing , then applying analogous calculations to those in §6.7 gives us the following result:
Employing PFEM with , in order to match the convergence rate of the Cauchy–Born error term , the highest mesh coarsening rate is
This means that the optimal mesh grading that PFEM can achieve without compromising accuracy is no greater than . However, in that case, the number of degrees of freedom is always .
In summary, despite the possibility of (slightly) reducing the number of degrees of freedom, considering its computational complexity, we conclude that higherorder FEM is not worthwhile to pursue. However, we emphasize that this conclusion would need to be revisited if a coupling method with higherorder interface error as well as continuum model error could be devised.
4. Conclusion
We obtained a sharp energynorm error estimate for the G23 coupling method with P2FEM discretisation of the continuum model. Furthermore, we demonstrated that, with P1FEM discretisation the FEM coarsening error is the dominating term in the consistency error estimate, whereas for P2FEM discretisation the interface error becomes the dominating term. In particular, a P2FEM discretisation yields a more rapid decay of the error. Crucially though, since for PFEM with the interface contribution dominates the total error the P2FEM is already optimal. That is, increasing to will not improve the rate of convergence, but increase the computational cost and algorithmic complexity.
Numerically, we observe that the improvement of P2GR23 over P1GR23 is only modest at low , hence a P2GR23 scheme would be primarily of interest if high accuracy of the solution is required. However, our numerical results indicate that there is scope for further optimisation, using a posteriori type techniques.
While our estimates for the error in energynorm are sharp, our numerical results show the estimates for the energy errors are suboptimal. We hightlight the leading term in the error analysis which overestimate the error in Section 7.2. We are unable, at present, to obtain an optimal energy error estimate. This appears to be an open problem throughout the literature on hybrid atomistic multiscale schemes; see e.g. [4].
In summary we conclude that using P2FEM is a promising improvement to the efficiency of a/c coupling methods, but that some further work, both theoretical and for its implementation may be need to exploit its full potential.
5. Reduction to consistency
Assuming the existence of an atomistic solution , we seek to prove the existence of satisfying
(5.1) 
and to estimate the in a suitable norm.
The error analysis consists of consistency and stability estimates. Once these are established we apply the following theorem to obtain the existence of a solution and the error estimate. The proof of this theorem is standard and can be found in various references, e.g. [13, Lemma 2.2].
[The inverse function theorem] Let be a subspace of , equipped with , and let with Lipschitzcontinuous derivative :
where denotes the operator norm.
Let satisfy
(5.2)  
(5.3) 
such that satisfy the relation
Then there exists a (locally unique) such that ,
To ensure Dirichlet boundary conditions, we adapt the approximation map defined in [4]. Let be a cutoff function such that
For , define
(5.4) 
Let be the vertices of and be the midpoint of an edge . Then, the set of all active P2 finite element nodes is given by
This includes all P1 nodes as well as the P2 nodes (edge midpoints) associatd with edges entirely in the P2 region.
Furthermore, let be the interpolation operator such that, for , for , , and
We also introduce ghost nodes on the edges shared by interface and continuum elements:
(5.5) 
Then, for , , where and are the vertices of the edge on which lies. Hence, the and interpolants coincide on .
We can now define the projection map (quasibest approximation operator) as
(5.6) 
5.1. Stability
To put Theorem 5 (Inverse Function Theorem) into our context, let
To make (5.2) and (5.3) concrete we will show that there exist such that, for all ,
Ignoring some technical requirements, the inverse function theorem implies that, if is sufficiently small, then there exists such that
Finally adding the best approximation error gives the error estimate
The Lipschitz and consistency estimates require assumptions on the boundedness of partial derivatives of . For , define the first and second partial derivatives, for , by
and similarly for the third derivatives . We assume that the second and third derivatives are bounded
(5.7)  
(5.8) 
With the above bounds it is easy to show that
(5.9) 
From the bounds above we can obtain the following Lipschitz continuity and stability results.
There exists such that
(5.10) 
Proof.
The result follows directly from the global bounds of derivatives of and the fact that and that is compactly supported hence is also Lipschitz. ∎
Under the assumptions (A1) and (A2), if , then there exits such that, when is sufficiently large,
(5.11) 
Proof.
The proof of this result is a straightforward adaption of the proof of [5, Lemma 4.9], which is an analogous result for blendingtype a/c coupling. ∎
6. Consistency estimate with a P2FEM
6.1. Outline of the consistency estimate
We begin by decomposing the consistency error into
(6.1) 
where is given and we can choose arbitrarily.
For , for . But the test function in is a piecewise linear lattice function. While we postpone the construction of , we will ensure that it is defined in such a way that for all , where is an extra layer of atomistic sites outside . With this assumption in place, we can further decompose into the following parts,
(6.2)  
where is the smooth interpolant of defined in Lemma 6.1 below. By in we mean the gradient of the canonical linear interpolant of . To estimate we require an approximation error estimate for . To estimate we will exploit the fact that the atomistic triangulation is uniform to prove a superconvergence estimate. Finally, for the modelling error, , we employ the techniques developed in [12].
To define the smooth interpolant , we use the construction from [5], namely a conforming multiquintic interpolant. Although the interpolant defined in [5] is for lattice functions on , we can use the linear transformation from to to obtain a modified interpolant.
(a) For each , there exists a unique such that, for all ,