Access the full text.
Sign up today, get DeepDyve free for 14 days.
denis.davydov@fau.de Chair of Applied Mechanics, In this paper the h-adaptive partition-of-unity method and the h-and hp-adaptive ﬁnite University of element method are applied to eigenvalue problems arising in quantum mechanics, Erlangen-Nuremberg, Egerlandstr. 5, 91058 Erlangen, namely, the Schrödinger equation with Coulomb and harmonic potentials, and the Germany all-electron Kohn–Sham density functional theory. The partition-of-unity method is Full list of author information is equipped with an a posteriori error estimator, thus enabling implementation of available at the end of the article error-controlled adaptive mesh reﬁnement strategies. To that end, local interpolation error estimates are derived for the partition-of-unity method enriched with a class of exponential functions. The eﬃciency of the h-adaptive partition-of-unity method is compared to the h-and hp-adaptive ﬁnite element method. The latter is implemented by adopting the analyticity estimate from Legendre coeﬃcients. An extension of this approach to multiple solution vectors is proposed. Numerical results conﬁrm the theoretically predicted convergence rates and remarkable accuracy of the h-adaptive partition-of-unity approach. Implementational details of the partition-of-unity method related to enforcing continuity with hanging nodes are discussed. Keywords: Adaptive ﬁnite element method, Partition-of-unity method, Error estimators, Schrödinger equation, Local interpolation error estimates, Density functional theory Introduction Recently there has been an increase of interest in applying ﬁnite element (FE) methods to partial diﬀerential equations (PDEs) in quantum mechanics [1–14]. In order to improve the accuracy of the solution, the basis set can be adaptively expanded through either reﬁnement of the mesh (h-adaptivity) or the basis functions can be augmented by the introduction of higher polynomial degree basis functions (p-adaptivity). Since the solution is not smooth and contains cusp singularities, the application of the h-adaptive FEM may require very ﬁne meshes and could be computationally ineﬃcient. There are several approaches to circumvent this problem. From the physical point of view, for ab initio calculation of molecules often core electrons (as opposed to valence electrons) behave in a similar way to single atom solutions. Thus © The Author(s) 2017. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 0123456789().,–: vol Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 2 of 23 one possesses an a priori knowledge of a part of the solution vectors to the eigenvalue problem. One of the approaches used to introduce this into a FE formulation is the partition-of-unity method (PUM) [15–17], which is a generalization of the classical FE method. In PUM the enrichment functions are introduced into a basis as products with standard FE shape functions, thereby enlarging the standard FE space. As the standard FE functions satisfy the partition-of-unity property (that is, they sum to one in the whole domain), the resulting basis can reproduce enrichment functions exactly. For an overview on PUM applied to continuum mechanics we refer the reader to [18–20]. An alternative approach to the above is to combine h-and p-adaptivity resulting in what is termed as hp-adaptive FEM. For an overview of hp-adaptive reﬁnement strategies we refer the reader to [21]. The general idea is that when the exact solution is smooth on the given element, p-adaptive reﬁnement is more eﬃcient and leads to a faster convergence with respect to the number of degrees of freedom; whereas if the solution is non-smooth (singular), h-adaptive reﬁnement is performed. Thus in addition to a reliable error estimate and the choice of the marking strategy of elements for reﬁnement, hp-adaptive methods need to decide which type of reﬁnement to perform on a given element. In this work we use methods based on smoothness estimation [22–27]. As those methods are normally employed for problems with a single solution vector, we propose an extension to multiple solution vectors as is required for the here considered eigenvalue problems. Herein, our main focus is application of h-adaptive PUM to PDEs in quantum mechanics, namely to the Schrödinger equation and the all-electron density functional theory (DFT) [28,29]. Application of the PUM to the above problems holds a signiﬁcant promise to improve on accuracy of a standard (non-enriched) FE approximation. The corresponding numerical evidence can be found in [9], where convergence studies for PUM solutions obtained on uniformly reﬁned meshes are performed. The novelty of our paper is that the PUM will be equipped with an a posteriori error estimator, thus enabling implementation of error-controlled adaptive mesh reﬁnement strategies. Derivation and implementation of the PUM in computational solid mechanics is nowadays very well-acknowledged and established area of research, yet the authors are not aware of any other work which applies the h-adaptive PUM to DFT. We will also compare the PUM to hp-adaptive FEM in terms of the eﬃciency with respect to the number of degrees of freedom. Although there are publications on the topic of hp-adaptive FEM applied to DFT [1], they lack any numerical studies and are limited to a pre-deﬁned reﬁnement strategy of hexahedra that admit nuclei only at its vertices. In order to apply the hp-adaptive FEM to DFT, in this paper we propose an extension of the smoothness estimate approach using Legendre coeﬃcients [22–25] to multiple solutions vectors. The outline of this paper is as follows: In the section on “Theory”, we introduce the eigenvalue problem studied here. The PUM and error estimators are also discussed. We also explain the strategy employed to decide between h-and p-adaptive reﬁnement for the hp-adaptive FEM. Results of numerical studies of the chosen systems are presented in section titled “Results and discussion”, followed by some conclusions. In Appendix A we rigorously derive the local interpolation error estimates for enrichment with a class of exponential functions; Appendix B describes the approach applied to solve single atom DFT in radial coordinates within application of the PUM; Appendix C discusses imple- mentational details of PUM within the context of the deal.II [30] library. Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 3 of 23 Theory In quantum mechanics we seek the N lowest eigenpairs (λ , ψ ) of the Schrödinger equa- α α tion [31] − ∇ + V (x) ψ (x) = λ ψ (x)on , α α α ψ (x) =0on ∂, (1) ψ (x)ψ (x)dx = δ . α β αβ 1 3 Here α is the index of the eigenpair, is a Lipschitz domain in R and δ is the Kronecker αβ delta . In Kohn–Sham all-electron density functional theory [28,29], the potential V (x) depends on eigenvectors thus rendering the problem nonlinear. For a molecular system consisting of N electrons and M nuclei of charges {Z } located at the (ﬁxed) positions {R }, e I I the ground state electron density ρ(x):= f ψ (x) can be obtained by ﬁnding the α α α=1 N lowest eigenpairs of Eq. (1). Here f is the partial occupancy number of the α-orbital such that f = N , V = V + V + V is composed of the ionic potential α e ion Hartree xc α=1 M Z V =− , the Hartree potential V = ρ(x )/|x−x |dx , and the (given) ion Hartree I=1 | | x−R exchange-correlation potential V (ρ). As a result, the potential V depends on the density xc ρ which is given in terms of eigenvectors {ψ }, making the problem nonlinear. From practical perspective V electrostatic potential is obtained by solving the associated Hartree Laplace equation; together with (1) they are solved sequentially untill convergence in density ﬁelds ρ is attained. For further details on the FE solution of DFT, we refer to our previous work [2] and literature cited therein. TheweakformofEq. (1)reads ∇v ·∇ψ + vV ψ dx = λ vψ dx ∀v ∈ H () , α α α α (2) ψ ψ dx = δ . α β αβ We then introduce a FE triangulation P of and the associated FE space of continuous h 1 piecewise elements of a ﬁxed polynomial degree : ψ ∈ V ⊂ H (). The FE solution to the problem is then deﬁned by h h h h h h h h h ∇v ·∇ψ + v V ψ dx = λ v ψ dx ∀ v ∈ V , α α α α (3) h h ψ ψ dx = δ . αβ α β 1 3 Note that in the non-periodic case the Schrödinger equation is actually set in R . Therefore the domain in Eq. (1) is assumed to be suﬃciently large such that zero Dirichlet boundary conditions make sense and there is no additional error due to considering a bounded domain. For all example systems that are considered below, the eigenfunctions are known to have asymptotic exponential decay which allows one to choose moderately sized domains. For spin unpolarized systems with an even number of electrons f ≡ 2 and N = N /2. α e We use standard notation for Sobolev spaces and norms. Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 4 of 23 Partition-of-unity method The classical FEM with piecewise linear ansatz spaces requires very ﬁne meshes for ade- quate accuracy when the solution is not smooth or is highly oscillatory; this increases the computational cost of solving the problem. The PUM proposed by Melenk and Babuska in [15,16] can address this issue. The main feature of the PUM is the inclusion of an a priori knowledge about the solution properties into the FE space. The PUM enriches the vector space spanned by standard FE basis functions N (x) (e.g. polynomials) by products of these functions with functions f (x) that contain a-priori knowledge about the solution ⎡ ⎤ ij h i ⎣ ⎦ ψ (x) = N (x) ψ + f (x)ψ . (4) i j α α α i∈I j∈S ij Here ψ are standard degrees-of-freedom (DoFs) and ψ are additional DoFs associated with the shape functions N (x) and the enrichment functions f (x); I is a set of all nodes i j and S is the set of enrichment functions. Since (possibly global) enrichment functions f (x) are multiplied with N (x) which has local support, the product also has local support and therefore matrices arising from the weak form remain sparse. Also, since the standard shape functions satisfy the partition of unity property N (x) ≡ 1, the resulting vector space can reproduce enrichment functions f (x)exactly. Note that (4) is a more general approach than enriching the basis with f alone (i.e. without multiplying by N , as is employed in [32]). Granted the partition of unity property of N , this case can be obtained from (4) by requiring all DoFs associated with a given enrichment function f to have the same value. Error estimator A posteriori error estimation analysis for FE approximations of (second-order) eigen- value problems has been a topic of intensive study within the last several decades, both from theoretical and implementational standpoints. We refer the interested reader to [13,14,33–39], where two “conventional” types of error estimators, namely residual- and averaging-based error estimators, are presented. In general, a discretization error in approximated eigenfunctions, ψ − ψ ,measuredina suitable norm (e.g. L -norm and energy norm, induced by the bilinear form of a problem), as well as in approximated eigenvalues, |λ − λ |, can be estimated from above. That is, ψ − ψ ≤ C η, (5) and h 2 |λ − λ |≤ C η , (6) where C ,C are the so-called stability constants that are independent of the mesh size 1 2 and η is the explicitly computable error upper-bound, see e.g. [34,38] for details. These equations are typically termed (global) error estimators. The bound η reads as That is, η depends solely on the FE solution. Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 5 of 23 ⎡ ⎤ ⎣ ⎦ η := η , K ∈P where summation is performed over all elements in P and η is the (local) error indicator, h h a quantity showing a discretization error of {ψ , λ } element-wise, that is, on every ﬁxed h h K . With multiple solutions available (in this case, eigenpairs {ψ , λ }), η will be a sum α α of discretization errors of the corresponding eigenpairs on a given element K , that is η := η . K,α For a standard (non-enriched) Q -based ﬁnite element solution of (1), a local indicator η of residual type reads as follows (see [13,14,34,35,38,39] for details): K,α 2 2 2 h h h η := h − ∇ + V (x) ψ − λ ψ dx K,α K α α α + h − ∇ψ · n da, (7) e⊂∂K 1 h 1 h 1 h where [[− ∇ψ ·n]] := − ∇ψ | + ∇ψ | ·n represents the jump of the gradient e K K e α α α 2 2 2 across interface e between two adjacent elements K and K , n is the outward unit normal vector to e and h := diam(K ). One of the key ﬁndings of our work is the proof that indicator (7) also holds (with no modiﬁcation due to the enrichment usage) in the PUM with the exponential enrichment function f (x) = exp (−μ |x| ). In Appendix A, we derive and prove the related local interpolation error estimates required for the derivation of the error indicator (7). hp-adaptive solution There have been numerous works devoted to hp-adaptive reﬁnement [22–25,40–42] including a comparison of diﬀerent methods [21]. The main diﬃculty that a posteriori hp-adaptive methods aim to address is the following: once an error is estimated and a certain subset of elements is marked for reﬁnement, one has to choose between h-or p-reﬁnement for each element. In this work we adopt an hp-reﬁnement method based on the estimate of the analyticity of the solution on the reference element via expansion into Legendre bases [22–25]. In particular, the FE solution is analytic on element K if, and only if, there exists constants C and σ such that K K a ≤ C exp(−σ [i + j + k]), (8) K K ijk where a are Legendre coeﬃcients; see [25] for further details. We chose to estimate ijk the decay coeﬃcient σ by performing a least squares ﬁt of Legendre coeﬃcients in K d d each direction ln a ∼ ln C − σ i for 1 ≤ i ≤ p , and then use the minimum decay K K d,i coeﬃcient as the ﬁnal value σ = min σ . When σ value is below a chosen parameter σ , K d K 0 the solution is considered to be smooth at K and thus p-reﬁnement is performed, otherwise That is the measure of how well it is representable by power series. Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 6 of 23 h-reﬁnement is executed. For initially linear FEs p-reﬁnement is always performed. We note that methods based on the decay rate of the expansion coeﬃcients were found in [21] to be the best choice as a general strategy for the hp-adaptive solution of elliptic problems. To the best of our knowledge there is, however, very little (numerical) study of those methods applied to DFT. The only paper we are aware of [1] lacks any numerical results. We also note that, in the majority of cases, hp-adaptive FEM is applied to problems in R and R . Thus we also aim to evaluate how well the smoothness estimators proposed in the literature work for eigenvalue problems in R that are relevant to quantum mechanics. In order to extend this hp-reﬁnement strategy to the eigenvalue problem, that is when there are multiple vectors represented using the same FE basis, we propose the following approach. For each element we ﬁnd an eigenvector which contributes the most to the total element’s error. The smoothness of this vector is the basis on which we decide to perform h-reﬁnement or p-reﬁnement. The rationale behind this approach is that we aim at minimizing the error the most during a single reﬁnement step while being conservative and avoiding performing both h-and p-reﬁnement on the same element. We also investigated allowing both h-and p-adaptive reﬁnement of a single cell based on smoothness estimation of all eigenvectors, but ultimately found that this procedure leads to qualitatively similar results for the problems studied herein. Finally, for the error indicator we adopt the following expression [43] h 1 2 K 2 h h h η := − ∇ + V (x) ψ − λ ψ dx K,α α α α h 1 + − ∇ψ · n da , (9) 2p 2 e e e⊂∂K where h is the face’s diameter, p is the element’s polynomial degree and p is the e K e maximum polynomial degree over two elements K and K adjacent to the face e. The derivation of the error estimators for hp-FEM usually requires the polynomial degree of neighbouring elements to be comparable, namely that there exists γ and such that γ p ≤ p ≤ p for all elements K, K that have a non-empty intersection. In K K K order to reﬂect this assumption in our numerical scheme, we propose that an additional step which limits the diﬀerences in polynomial degrees among elements be performed. More precisely, after hp-adaptive reﬁnement is executed, then for each element K,we max ﬁnd the maximum polynomial degree among its neighbouring elements p and if neigh max p > p + 2thenweset p ← p + 1. K K K neigh Results and discussion If not explicitly stated otherwise, the results below are obtained for the following con- ﬁguration: (i) the initial polynomial degree for non-enriched DoFs is one for hp-adaptive FEM; (ii) linear shape functions are used in PUM to introduce enrichments, higher order elements were not employed as the interpolation error estimates are derived only for lin- ear elements and thus limit the applicability of the error indicator stated in Eq. (7) ; (iii) a Gaussian quadrature rule with 20 points is used for enriched elements in the eigenvalue problem; (iv) the Dörﬂer marking strategy with θ = 0.6 is used to mark elements for reﬁnement; (v) Gauss–Legendre–Lobatto supports points are used for the hp-adaptive FEM basis to improve the condition number; (vi) in case of hp-adaptive reﬁnement the Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 7 of 23 highest polynomial degree is limited to 8 for computational eﬃciency reasons; (vii) the radius in which enriched FEs are employed is heuristically chosen for each numerical example; (viii) following [25] we choose σ = 1.0 as a parameter in the smoothness esti- mator. Implementation details of the partition-of-unity method in deal.II [30] ﬁnite element library are given in Appendix C. Schrödinger equation In this section we consider the Schrödinger equation Eq. (1) with two diﬀerent (spherical) potentials V (x) = V (|x|). The ﬁrst case is the Coulomb potential V (x) =− 1/ |x|, which corresponds to a Hydrogen atom. The eigenvalues of this problem are degenerate. In R , 2 2 on each energy level n there are n eigenvalues λ = λ /n , where λ =− 1/2[31]. The n 1 1 eigenfunction corresponding to the lowest eigenvalue reads ψ (x) = √ exp (− |x|) . (10) The radial component of the eigenfunctions at the next energy level are R = [1 − 2,0 |x| /2] exp(− |x| /2) and R = |x| /2exp(− |x| /2). 2,1 The second potential we will consider is a harmonic potential V (x) = |x| /2thatleads to a harmonic oscillator problem. The eigenvalues for this problem are also degenerate; in R they are given by λ = n + 1/2 for nth energy level. The lowest two have a degeneracy of 1 and 3, respectively. The (unnormalized) eigenfunction corresponding to the lowest eigenvalue is ψ (x) = exp − |x| /2 . (11) The radial component of the next eigenfunction is R (x) = |x| exp − |x| /2 .Figure 1 0,1 shows radial components of eigenfunctions for the Coulomb and the harmonic potential. It is clear that in order to have a low interpolation error for a standard Lagrange FE basis, a very ﬁne mesh will be required near the origin. For such non-smooth solutions we will see that by introducing enrichment functions the interpolation error of the resulting FE basis will be greatly reduced. The initial mesh used to solve the Schrödinger equation is obtained from 3 global mesh reﬁnements of the single element in = [− 20; 20] for the Coulomb potential and = [− 10; 10] for the harmonic potential. For the PUM only 8 elements adjacent to the singularity that is located at the origin are marked for enrichment. First, we examine the convergence in case when a single eigenpair is required in the Schrödinger equation with two diﬀerent potentials. Figure 2 compares the h-adaptive FEM, hp-adaptive FEM and h-adaptive PUM, whereas Fig. 3 shows the cross-sections of meshes for the last reﬁnement step. For both combinations of potentials and enrichment functions, the h-adaptive PUM is superior to h-adaptive FEM. In particular, for the last reﬁnement step the PUM solution is about two orders more accurate than the h-adaptive FEM with the same number of DoFs in the case of the Coulomb potential. For the harmonic potential this value is smaller. The Note that in [22] the Legendre coeﬃcients were required to have even slower decay rate of σ = 0.69. For spherically symmetric potentials one can separate eigenfunctions into radial R (r) and angular Y (θ, φ) parts, n,l m,l where the latter are spherical harmonics [31]. Here {n, l, m} are three quantum numbers. Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 8 of 23 a b 1 1 R R 1,0 0,0 R2,0 R0,1 0.9 2,1 0.8 0.8 0.7 0.6 0.6 0.4 0.5 0.4 0.2 0.3 0.2 0.1 -0.2 0 0 2 4 6 8 10 0 1 2 3 4 5 6 |x| |x| Fig. 1 Radial components of eigenfunctions for diﬀerent potentials V(x). The dotted vertical line indicates the smallest initial mesh size which will be used in our numerical calculations. a Coulomb. b Harmonic 1 2 a b 10 10 h-adaptive FEM (linear) h-adaptive FEM (linear) hp-adaptive FEM hp-adaptive FEM 0 h-adaptive PUM h-adaptive PUM h-adaptive FEM (quadratic) h-adaptive PUM (larger enrichment domain) a line with slope -1.88 a line with slope -0.75 -1 a line with slope -0.79 10 -2 -2 -4 -3 -6 -4 -8 -5 3 4 5 2 3 4 5 10 10 10 10 10 10 10 DoFs DoFs Fig. 2 Error convergence rates for an eigenproblem with a single eigenpair. a Coulomb potential. b Harmonic potential asymptotic convergence rate of the h-adaptive PUM with the default enrichment radius is very similar to that of the h-adaptive FEM for both problems (compare green and red linesinFig. 2), which supports our theoretical ﬁndings. The advantage of the h-adaptive PUM also depends on the enrichment radius with respect to the underlying exact solution. To examine this eﬀect we employ an initial mesh obtained only by two global reﬁnements of a single element and mark the 8 elements adjacent to the origin for enrichment. With this approach we eﬀectively consider a larger 3 3 enrichment domain [− 5; 5] instead of [− 2.5; 2.5] . Importantly, the numerically non- zero part of the underlying analytical solution will be almost fully contained in those 8 elements (see Fig. 1b). From the numerical results we observe that for the most reﬁned stage the h-adaptive PUM displays an error which is about 6 orders of magnitude less than the same method with the smaller enrichment domain (compare purple and green lines in Fig. 2b). For the case of a single eigenpair, the hp-adaptive FEM performs remarkably well and, unless a larger enrichment radius is used in h-PUM, it converges to the higher tolerance with fewer number of DoFs (compare blue and green lines in Fig. 2). Now let us turn our attention to a more realistic scenario where one seeks multiple eigenpairs whereby an a priori knowledge is available only for the ﬁrst eigenfunction. Figure 4 plots the convergence of the ﬁrst 5 eigenvalues for the Coulomb potential and λ − λ λ − λ 0 0 Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 9 of 23 a b degree degree 1 1 2 2 4 4 5 5 7 7 degree Fig. 3 Cross-sections of the ﬁnal adaptive meshes for the Coulomb potential when solving for a single eigenpair. a h-adaptive FEM (linear). b h-adaptive FEM (quadratic). c hp-adaptive FEM. d h-adaptive PUM (linear) the ﬁrst 4 eigenvalues for the harmonic potential for the diﬀerent methods. For both problems the h adaptive PUM again has remarkable convergence properties, superior to h-adaptive FEM. It is important to note that even though in the PUM the enrichment function corresponds to the ﬁrst eigenfunction only, other eigenpairs in the case of the harmonic potential tend to converge faster than the standard h-adaptive FEM case, as can be observed in Fig. 4b. The same applies to the spherical orbital at the second energy level of the Hydrogen atom; see Fig. 4a where the corresponding eigenvalue in the PUM case displays a faster convergence rate than the others on the same energy level. For the Hydrogen atom, in the case of the hp-adaptive reﬁnement one observes a supe- rior convergence rate of the ﬁrst eigenvalue, whereas eigenvalues from the next energy level have errors that are comparable to the h-adaptive linear FEM. A possible issue could be related to the smoothness estimation on elements with hanging nodes. In particular it is observed [44] that the smoothness is overestimated when using similar methods, albeit based on Fourier coeﬃcients. This leads to unnecessarily high order polynomial degrees in these areas. Clearly, further investigation is required to resolve this problem. Density functional theory Finally, we apply the here considered FE approaches to the Kohn–Sham density functional theory. As a ﬁrst test problem we consider a single He atom which has a single doubly Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 10 of 23 h-adaptive FEM (linear) hp-adaptive FEM h-adaptive PUM -1 -1 -1 10 10 10 -2 -2 -2 10 10 10 -3 -3 -3 10 10 10 -4 -4 -4 10 10 10 -5 -5 -5 10 10 10 3 4 5 3 4 5 3 4 5 10 10 10 10 10 10 10 10 10 DoFs DoFs DoFs 1 2-5 a line with slope -0.72 h-adaptive FEM (linear) hp-adaptive FEM h-adaptive PUM 1 1 1 10 10 10 0 0 0 10 10 10 -1 -1 -1 10 10 10 -2 -2 -2 10 10 10 -3 -3 -3 10 10 10 -4 -4 -4 10 10 10 -5 -5 -5 10 10 10 -6 -6 -6 10 10 10 3 4 5 3 4 5 3 4 5 10 10 10 10 10 10 10 10 10 DoFs DoFs DoFs 1 2-4 a line with slope -0.84 Fig. 4 Convergence of eigenvalues from the ﬁrst two energy levels for the Schrödinger equation in the course of adaptive reﬁnement. Red lines denote the lowest eigenvalue, whereas blue lines correspond to degenerate eigenvalues on the next energy level. a Coulomb potential (4 out of 5 eigenvalues are degenerate). b Harmonic potential (3 out of 4 eigenvalues are degenerate) occupied state, i.e. N = 2and N = 1. The ground state energy from the radial solution is E =− 2.834289. Enrichment functions for PUM are obtained from numerical solution of single atom Schrödinger equations, depicted in Fig. 5a. The atom is placed at the origin in the domain = [− 10; 10] with the homogeneous mesh of size h = 2.5. Eight elements adjacent to the atom are enriched. Figure 5b compares the h-adaptive FEM, hp-adaptive FEM and h-adaptive PUM. One immediately recognizes that the PUM leads to a much faster convergence in terms of DoFs and gives about an order of magnitude advantage in terms of the absolute value of the error. The linear h-adaptive FEM would require ten times more DoFs to achieve the same accuracy. The hp-adaptive FEM displays an exponential-like decay and approaches the accuracy of PUM at higher number of DoFs. In the second test problem we consider a CO molecule in the domain = [− 10; 10] at the (equilibrium) distance 2.1. In order to estimate the ground state energy, we ﬁt the total energy obtained by h-FEM at the last 3 mesh reﬁnement steps to ln(|E − E |) = h h λ − λα λ − λα α α Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 11 of 23 b 2 1 h-adaptive FEM (linear) 1s a line with slope -0.73 0.9 1 hp-adaptive FEM h-adaptive PUM 0.8 error estimator η PUM 0.7 0.6 -1 0.5 -2 0.4 0.3 -3 0.2 0.1 -4 3 4 5 0 10 10 10 0.0 0.5 1.0 1.5 2.0 2.5 3.0 DoFs Fig. 5 Finite element solution of He atom. a Scaled radial solution. The dotted vertical line indicates the enrichment radius. b Convergence of the error in total energy of He atom for various FE methods C + q ln(DoFs) with constraints C > 0,E < E ,q < 0. Using this approach we estimate the limit of the ground state energy to be E =− 112.47107. This renders a bond energy of − 0.5775, which compares favourably to the value − 0.578 reported in [45]. This gives us conﬁdence to use the estimated ground state energy value in convergence studies, which are presented in Fig. 6. The enrichment functions for PUM are obtained from the numerical solution of single atom Schrödinger equations; see Appendix B for details. The scaling of those functions are not important for PUM, so Fig. 6a depicts radial solutions normalized so that the value of the 1s and 2s orbitals are unity at the origin. It is generally possible to use all eigenfunctions from the radial solution as enrichments around each atom in the radius of a few atomic units. However, extra care must be taken not to render the resulting FE space to have linearly dependent basis functions. Figure 6a clearly indicates that given small enough elements (on the order 0.1 a.u.), enriching with both 1s and 2s single atom radial core electrons solutions would make the FE space degenerate. Our current implementation of PUM DFT only supports enrichment in non-overlapping domains. Therefore for the CO molecule we have to start from a relatively ﬁne mesh, which in the course of h-adaptive reﬁnement may render the basis enriched with multiple functions linearly dependent. To avoid this, the PUM results for the CO molecule are obtained by enriching 8 elements adjacent to each atom with its 1s orbital only. Scaling of the 1s function to unity at the origin of the enrichment spherical function improves the condition number of the resulting matrices. Figure 6b compares the convergence characteristics using the h-adaptive FEM, hp- adaptive FEM and h-adaptive PUM. The energy error convergence rate from h-adaptive 2p FEM compares favourably to the expected rate of O(h ), which can be approximated −2p/3 by O(DoFs ). Remarkably, the chosen smoothness estimate used in the hp-adaptive FEM and its extension to multiple vectors do not lead to an increase in eﬃciency in terms of the number of DoFs as compared to h-adaptive quadratic FEM. The h-adaptive PUM displays the same convergence rate as h-adaptive FEM and is, as expected, more accurate. This, however, comes at the expense of having a worse condition number for the resulting matrices and the necessity to use higher quadrature order to perform suﬃciently accurate numerical integration. For this example and the chosen enrichment radius, the diﬀer- Single-atom energies of C and O atoms are − 37.42426 and − 74.46933, respectively. E − E 0 Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 12 of 23 a b 3 1.6 h-adaptive FEM (quadratic) C 1s h-adaptive PUM (quadratic) C 2s 1.4 hp-adaptive FEM C 2p a line with slope -1.68 O 1s 1.2 2 error estimator η PUM O 2s O 2p 0.8 0.6 -1 0.4 0.2 -2 -3 -0.2 3 4 5 6 10 10 10 10 0 0.1 0.5 1.0 DoFs Fig. 6 Finite element solution of CO molecule. a Scaled radial solution of single atoms. The dotted vertical line at 0.5 indicates the enrichment radius. b Convergence of the error in total energy of CO molecule for various FE methods ence in energy error between the two approaches is less than one order of magnitude. By comparing these results to those presented earlier for H and He atoms, we hypothetize that a larger enrichment radius is required to make the PUM advantageous compared to the h-adaptive FEM. Our current implementation of PUM DFT, however, only allows enrichment in non-overlapping domains, which limited the enrichment radius for the CO example. Conclusions In this contribution we have applied and critically compared the h-and hp-adaptive FEM, and the h-adaptive PUM to the relevant PDEs in quantum mechanics, namely the Schrödinger equation and the Kohn–Sham all-electron density functional theory. The main ﬁndings are summarized below. • The PUM renders several orders of magnitude more accurate eigenvalues than the standard FEM when solving the Schrödinger equation for the lowest eigenpair with Coulomb and harmonic potential. For the case when more eigenpairs are sought but only the lowest eigenvector is introduced as an enrichment, the PUM is still more accurate, especially for the lowest eigenvalue. Remarkably other eigenvalues also exhibit a faster convergence. The results from DFT calculations indicate that in order to keep this advantage, a reasonably large enrichment radius is needed. • For problems where a single eigenpair is being sought, the hp-adaptive FEM with the here considered smoothness and residual error estimators results in a more accurate solution with fewer number of DoFs as compared to h-adaptive PUM and FEM. However, for the case of multiple eigenpairs this approach did not lead to satisfactory results. Overall we ﬁnd h-adaptive PUM to be a more robust solution method to reach the required accuracy even with relatively small enrichment domains. • Local interpolation error estimates are derived for the PUM enriched with the class of exponential functions. In this case the results are the same as for the standard FEM and thereby admit the usage of the error indicator (7). • For the PUM DFT calculations the convergence rate of energy error and the residual error estimator are the same for all studied examples. Thus our numerical results conﬁrm that Eq. (7) can be considered as a reliable error indicator for problems in quantum mechanics. E − E 0 Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 13 of 23 • An element view to the implementation of PUM in FEM codes based on hexahedra is proposed (see Appendix C). As a result, continuity of the enriched ﬁeld along the edges with hanging nodes is enforced by treating FE spaces produced by each func- tion in the local approximation space separately. The resulting algebraic constraints are independent on the enrichment functions. This allows one to directly reuse algo- rithms written for enforcing continuity of vector-valued FE spaces constructed from a list of scalar-valued FEs. Abbreviations h h x: position vector; λ : eigenvalue; ψ (x): eigenvector; V (x): potential; : Lipschitz domain; P : triangulation of ; V : α α the ﬁnite element space of ﬁxed polynomial degree; λ : eigenvalue obtained from the ﬁnite element solution of the h i problem; ψ (x): eigenvector obtained from the ﬁnite element solution of the problem; ψ : standard degrees-of-freedom; α α ij N (x): ﬁnite element shape functions; f (x): enrichment functions; ψ : enriched degrees-of-freedom; η: explicitly i j α computable error upper-bound; η : local error indicator on the element K ; η : local error indicator for eigenpair α on K K,α the element K ; a : Legendre coeﬃcients; σ : the decay coeﬃcient; h : the face diameter; p : the element’s polynomial ijk K e K degree; p : the maximum polynomial degree over two elements K and K adjacent to the face e. Authors’ contributions DD developed the concepts proposed and discussed in this manuscript, performed the numerical implementation and experiments and wrote the draft of this paper. TG derived local interpolation error estimates and improved the mathematical introduction of the problem. JPP contributed to the numerical implementation and paper revisions, and PS revised the paper. All authors read and approved the ﬁnal manuscript. Author details 1 2 Chair of Applied Mechanics, University of Erlangen-Nuremberg, Egerlandstr. 5, 91058 Erlangen, Germany, Institute of Applied Mechanics, Technische Universitaet Braunschweig, Bienroder Weg 87, 38106 Braunschweig, Germany. Acknowledgements Not applicable. Competing interests The authors declare that they have no competing interests. Availability of data and materials Not applicable. Consent for publication Not applicable. Ethics approval and consent to participate Not applicable. Funding The support of this work by the ERC Advanced Grant 289049 MOCOPOLY (DD, JP, PS) and German Science Foundation (Deutsche Forschungs-Gemeinschaft, DFG), Grant DA 1664/2-1 (DD) is gratefully acknowledged. Second author (TG) is supported by the European Research Council (ERC) Starting Researcher Grant INTERFACES, Grant Agreement No. 279439. Appendix A: Local interpolation error estimates In this appendix, the local interpolation error estimates required for the derivation of the error indicator (7) in the case of PUM are obtained for linear ﬁnite element approximations enriched with f (x) = exp (−μ |x| ), where 0 <μ ∈ R and 1 ≤ p ∈ N. These are v − q v ≤ c h |v| 1 , (12) K K H (ω ) L (K ) v − q v ≤ c h |v| , (13) H (ω ) K K L (e) where, as usual, v : → R is a scalar-valued function, which is assumed to be at least 1 h in H (), q is a quasi-interpolation operator (of the averaging type), K is an element of the discretization P of , e ⊂ ∂K is an edge of K.Also, h measures the size of K , ω is the patch of elements neighboring K including K itself. Finally, c ,c ∈ R are the K K e interpolation constants independent of the mesh size. Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 14 of 23 We ﬁx the notations to be used throughout the appendix and make assumptions that are conventional for this kind of analysis. For the sake of simplicity and without loss of generality, we elaborate here for the two-dimensional setting. The obtained results are valid in three dimensions as well. h 2 First, we assume that the partition P of ⊂ R consisting of open and convex quadri- laterals K is shape-regular (or non-degenerate), as well as locally quasi-uniform in the sense of [46,47]. For every K and its edge e we deﬁne h := diam(K)and h :=|e| is the K e length of e. For every node i in P we denote by ω the union of quadrilaterals connected to node i and set h := diam(ω ). Furthermore, for every K , ω represents the patch ω i K containing K and the ﬁrst row of its neighbors; it is then set h := diam(ω ). ω K Also, in what follows, by the notation a b we imply the existence of a positive constant C independent of a and b such that a ≤ Cb. Then a ∼ b means that a b and a b hold simultaneously. The symbol |·| will be used to denote either the H -seminorm (as e.g. in (12)and (13)) or the length of a linear segment in R or the area of a plane domain 1 1 2 2 in R . With these notations at hand, one can show that |K | ∼ h , |ω | ∼ h and K i ω |ω | ∼ h . Furthermore, the shape regularity of the mesh P ensures that h ∼ h , K ω e K whereas its local quasi-uniformity implies that h ∼ h ∼ h . K ω ω i K Finally, we also recall useful inequalities, which are • the Poincaré-type inequality (see e.g. [48]): 1 h v − vdx ≤ |v| 1 , ∀v ∈ H (ω), (14) H (ω) |ω| 2 π ω L (ω) where ω ⊂ R (n = 2, 3) is a Lipschitz domain and h := diam(ω); • the scaled trace inequality (e.g. in [49], Lemma 3.2): 1 1 2 2 v h v + h |v| , ∀v ∈ H (K ). (15) 2 2 1 e e L (e) L (K ) H (K ) Quasi-interpolation operator Herein, we construct an interpolation operator for obtaining the local error estimates (12) and (13). 1 h Let V := H () be an admissible space and V be its (enriched) ﬁnite element coun- terpart h h h V := v ∈ C(): v (x):= a N (x) + f (x) b N (x) i i i i i∈I i∈I + c N (x),a ,b ,c ∈ R ⊂ V, (16) i i i i i std. i∈I h std. where I is the set of all enriched nodes of P and I is the set of standard, i.e. non- h strd. enriched nodes of P ; I ∩ I =∅. Recall also that N in our case is the Q -shape i 1 function associated with node i and supported on ω . h h Explicit construction of the operator q : V → V implies the explicit pattern of assignments of a ,b ,c ∈ R through a function v ∈ V . In the case of the enriched FE i i i Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 15 of 23 approximation (16), the major challenge in deriving q is imposition of the constant- h h preserving property on q , which should be fulﬁlled on every element K ∈ P regardless the element type (see Fig. 7). h h The operator q : V → V with the desired property reads as follows: 1 1 q v(x):= v(y)dy N (x) + f (x) v(y)dy N (x) i i 2|ω | 2f (x )|ω | i ω i i ω i i i∈I i∈I + v(y)dy N (x), (17) |ω | strd. i∈I with all notations as in (16) and where x , entering the second term, denotes the coordinate of a node i. Below, for the proposed quasi-interpolation operator of the averaging type h h q we establish that q c| = c on a standard element (note this is a classical result for a non-enriched FEM) and, more importantly, that q c| = c + O(h ) on a fully-enriched andablendedelement. Estimates Preliminaries The three estimates that we start with are basic for the following local interpolation error analysis. On every K ∈ P and its node i it holds that N (x) h , N (x) h , (18) i 2 K i 2 L (K ) L (e) −1 v(y)dy h v 2 + |v| 1 , (19) K L (ω ) H (ω ) K K |ω | i ω and f (x) = 1 + O(h ). (20) f (x ) Results (18) rigorously follow from the isoparametric concept and related properties, see e.g. [50] for details. We note that they may be also derived in a less rigorous manner owing 1 1 1 1 2 2 2 2 to a boundedness of the basis function N on K along with |K | ∼ h and |e| = h ∼ h . i K e Fig. 7 Types of elements in mesh P with respect to the imposed enrichment Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 16 of 23 The inequality (19) is obtained as follows: −1 − v(y)dy ≤|ω | v(y) dy ≤|ω | v i i L (ω ) |ω | i ω ω i i −1 h v + |v| . 2 1 L (ω ) H (ω ) K K K Here we used the Cauchy–Schwarz inequality, |ω | ∼ h ∼ h and also the extension- i ω K related result v ≤ v . 2 2 L (ω ) L (ω ) i K Finally, to show (20) we explicitly use the properties of f (x). For any ﬁxed K , x ∈ K and x ∈ K being one of its nodes, we have the following upper bound estimate: p p exp μ max |x| exp μ min |x|+ h f (x) exp(μ|x | ) K x∈K x∈K = ≤ ≤ p p f (x ) exp(μ|x| ) exp μ min |x| exp μ min |x| x∈K x∈K p−1 = 1 + μp min |x| h + h.o.t. in min |x|,h K K x∈K x∈K p−1 = 1 + O min |x| h . (21) x∈K Notice that due to boundedness of min |x| for a given ﬁxed K , there always exists > 0 x∈K such that min |x|= h . Using this in (21), we obtain x∈K p−1 p−1 min |x| h = h , x∈K f (x) yielding, as a result, ≤ 1 + O(h ). f (x ) K The lower bound estimate can be found similarly: p p exp μ min |x| exp μ min |x| f (x) exp(μ|x | ) x∈K x∈K = ≥ ≥ p p f (x ) exp(μ|x| ) exp μ max |x| exp μ min |x|+ h x∈K x∈K p−1 = 1 − μp min |x| h + h.o.t. in min |x|,h K K x∈K x∈K p−1 = 1 + O min |x| h , x∈K f (x) and, eventually, ≥ 1 + O(h ). The result (20) then follows. f (x ) K h 2 Stability of q in L -norm The next step towards (12)and (13) implies obtaining the so-called stability result for the constructed q .Using (18)–(20) one straightforwardly shows that q v v + h |v| , (22) 2 1 L (ω ) H (ω ) K K L (K ) and 1 1 h 2 2 q v h v 2 + h |v| 1 . (23) L (ω ) H (ω ) K K 2 K K L (e) These estimates indeed hold for every K regardless of its type (standard, blended, enriched). Note that for a standard non-enriched FEM and the resulting interpolation operators, the estimates (22)and (23) are classical. We have obtained and proved them for our speciﬁc operator q adopted for the current enriched FEM setting. Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 17 of 23 Constant-preserving property of q The ﬁnal ingredient required for obtaining (12)and (13) is the determination of how “well” the constructed q reproduces the constant on an element K , depending on its type. This constant-preserving property of the operator is of major importance particularly in the case of enriched FEM. The required result on a standard (non-enriched) element K follows immediately. Indeed,inthiscase q v(x)| = v(y)dy N (x), K i |ω | i=1 4 h and the partition of unity N (x) = 1on K yields q c| = c, c = const. i K i=1 The situation on a fully-enriched and partly-enriched (blended) element is more delicate. In the case of a fully enriched element we have q v(x)| = v(y)dy N (x) K i 2|ω | i ω i=1 + f (x) v(y)dy N (x), 2f (x )|ω | i i i=1 that, owing to (20), results in 4 4 1 1 f (x) 1 1 p p q c| = c + c N (x) = c + c 1 + O(h ) N (x) = c + O(h ). K i i K K 2 2 f (x ) 2 2 i=1 i=1 Now, let K be a blended element, implying the representation: q v(x)| = v(y)dy N (x) K i 2|ω | i=1 + f (x) v(y)dy N (x) 2f (x )|ω | i i ω i=1 + v(y)dy N (x), |ω | i ω i=+1 where ∈{1, 2, 3} is the number of enriched nodes of K . Adding and subtracting the ﬁrst sum in the above expression, enables us to rewrite it as follows: q v(x)| =− v(y)dy N (x) K i 2|ω | i ω i=1 + f (x) v(y)dy N (x) 2f (x )|ω | i i ω i=1 + v(y)dy N (x). |ω | i=1 Note that the last term contains the summation over all four nodes and is the standard (non-enriched) FE contribution which will automatically reproduce a constant. We then Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 18 of 23 need to estimate, in this context, the remaining part constituting of the ﬁrst and the second sums. We obtain, 1 f (x) 1 f (x) q c| = c − 1 N (x) + c ≤ c − 1 N (x) + c K i i 2 f (x ) 2 f (x ) i i i=1 i=1 4 4 1 f (x) 1 p p ≤ c − 1 N (x) + c = c O(h ) N (x) + c = c + O(h ), i i K K 2 f (x ) 2 i=1 i=1 where (20) was also used. Proof of local error estimates (12), (13) h h The derivation of the estimates for v − q v and v − q v is based on a com- 2 2 L (K ) L (e) bined use of the above stability results for q , the Poincaré and the scaled trace inequalities (14)and (15), respectively, as well as the constant-preserving property results. First, due to linearity of q ,wehave h h h v − q v = v − c − q (v − c) + c − q c 2 2 L (σ ) L (σ ) h h ≤ v − c 2 + q (v − c) + c − q c , (24) L (σ ) 2 2 L (σ ) L (σ ) where c = const and where, for the sake of brevity, we set σ ={K, e}.Weare nowina position to dissect every term in (24) in either case of σ . When σ = K in (24): By the Poincaré inequality (14), it holds that v − c 2 ≤ v − c 2 h |v| 1 , (25) L (K ) L (ω ) H (ω ) K K −1 where one can choose c =|ω | vdx and use h ∼ h . K ω K By the stability estimate (22) and the Poincaré inequality, it holds similarly to the above that q (v − c) v − c 2 + h |v − c| 1 h |v| 1 . (26) K K L (ω ) H (ω ) H (ω ) K K K L (K ) Furthermore, using the results of “Constant-preserving property of q ” section we obtain c − q c ≡ 0, if K is standard, (27) L (K ) and p+1 c − q c = O(h ), if K is fully enriched or blended. (28) L (K ) In the former case we also use that 1 2 =|K | ∼ h . L (K ) Using (25)–(28)in(24), the resulting local interpolation error of type (12) follows. Note p+1 that in the case of fully enriched and blended elements the term O(h ) that appears in the corresponding upper bound can be neglected, being the higher order term with respect to the leading one h |v| . H (ω ) When σ = e in (24): Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 19 of 23 By the scaled trace inequality (15), it holds 1 1 1 2 2 2 v − c 2 h v − c 2 + h |v − c| 1 h |v| 1 , (29) e e L (e) L (K ) H (K ) K H (ω ) where we also use h ∼ h along with result in (25). e K By the stability estimate (23) and the Poincaré inequality (14), we obtain the result that 1 1 1 2 2 2 q (v − c) h v − c 2 + h |v − c| 1 h |v| 1 . (30) K L (ω ) K H (ω ) K H (ω ) K K K L (e) Finally, using the results of “Constant-preserving property of q ” section we derive c − q c ≡ 0, if K is standard, (31) L (e) and p+ h 2 c − q c = O(h ), if K is fully enriched or blended. (32) L (K ) 1 1 2 2 In the former case we also use the fact that 1 =|e| = h ∼ h . L (e) Using (29)–(32)in(24), the resulting local interpolation error estimate of type (13) follows as well. Again, in the case of fully enriched and blended elements the term O(h ) that appears in the corresponding upper bound can be neglected, being the higher order | | term with respect to the leading one h v 1 . H (ω ) K K Appendix B: Single atom radial solution The radial solution of a single atom with charge Z is obtained by solving the following coupled problem 1 d d 2 φ r R = 4πρ (33) r dr dr 1 1 d d [l + 1] l ψ ψ − r + + V (ρ,r) R = R (34) nl 2 2 nl nl 2 r dr dr 2r 2l + 1 2 ρ = f R (35) nl nl 4π where n and l are the main and azimuthal quantum numbers, f are occupation numbers, nl V is the eﬀective potential and R and R are respectively radial components of the rl eigenfunctions and electrostatic potential. The radial component of the eigenfunctions ∞ ψ ψ ψ 2 2 are normalized by r (R ) dr = 1. On the change of variables U := rR and 0 nl nl nl φ φ U := rR , the system takes the form U = 4π r ρ (36) dr 1 d [l + 1] l ψ ψ − + + V (ρ,r) U = U (37) nl nl nl 2 2 2 dr 2r 2l + 1 2 ρ = f U (38) nl 2 nl 4πr l Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 20 of 23 which together with the following Dirichlet boundary conditions U (0) = 0 (39) nl U (∞) =0(40) nl U (0) = 0 (41) U (∞) = Z (42) can be solved using the standard ﬁnite element method in a suﬃciently large but ﬁnite domain. The eigenfunctions are then given by ψ = R Y , where Y are spherical nlm lm lm nl harmonics. Appendic C: PUM implementational details An enriched ﬁnite element class has been implemented for the general purpose object- oriented C++ ﬁnite element library deal.II [30]. The implementation is based on the FESystem class, which is used to build ﬁnite elements for vector valued problems from a list of base (scalar) elements. What diﬀers from that class is that the developed FE implementation is scalar, but built from a collection of base elements and enrichment functions ⎡ ⎤ ⎢ ⎥ u(x) = N (x)u + f (x) N (x) u , (43) ⎣ ⎦ i i k jk jk pum i∈I k∈S j∈I pum where I is the set of all DoFs with standard shape functions (see Fig. 8a), I is the set of all DoFs corresponding to shape functions enriched with f (x)(seeFig. 8b) and S is the set of enrichment functions. As distribution of DoFs in deal.II is element based, we always enrich all DoFs on the element. To restore C continuity between enriched and non-enriched elements, additional algebraic constraints are added to force DoFs u associated with N f on the jk jk k face between the enriched and non-enriched elements to be zero. This is equivalent to enriching only those shape functions whose support is contained within the enriched elements. The h-reﬁnement in deal.II is implemented using hanging nodes. In this case, extra algebraic constraints have to be added to make the resulting ﬁeld conforming. We build these constraints separately for the non-enriched FE shape functions and enriched shape functions; that is, the following spaces are separately made conforming: {N (x)}, {N (x)}, i j0 {N (x)}, etc. To illustrate this idea consider two separate FE spaces shown in Fig. 8. j1 We assume that functions in the ﬁrst space are non-zero everywhere in the domain, whereas functions in the second space are non-zero only in the left part, marked by the blue shading. Therefore we do not have to introduce any DoFs in the right part, the underlying elements are denoted by Q . The standard procedure implemented zero in deal.II [27] will enforce continuity of the vector ﬁeld by introducing algebraic constraints for DoFs associated with hanging nodes (3, 5, 17, 19), plus constraints for This is a generalization of (4) which allows to use diﬀerent FE spaces for each enrichment function. In practice one uses linear shape functions for enriched DoFs and possibly higher order shape functions for non-enriched DoFs. For linear FEs, the value at the hanging node is the average of the values at adjacent vertices, for example u = 1/2[u + u ]. 8 2 Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 21 of 23 ab Fig. 8 Treatment of hanging nodes for the h-adaptive PUM. Q denotes (bi)linear FE, whereas Q denotes 1 zero elements on which functions in the FE space associated with the enrichment function f (x) are zero and thus no DoFs need to be introduced. a First FE space (standard). b Second FE space (enrichment) DoFs 14, 16, 22 to make functions in the second FE space zero at the interface between Q and Q . We can observe now that if we take the constrained scalar ﬁeld from 1 zero the ﬁrst FE space and add a scalar ﬁeld from the second FE space multiplied by the enrichment functions f (x) (continuous in space), the resulting scalar FE ﬁeld will also be continuous. Thus we arrive at a conforming h-adaptive PUM space where only some elements are enriched. With reference to Fig. 8, the resulting PUM ﬁeld will have enrich- ment associated with DoFs 23, 24, 20, 21, 17, 18, 15 whereas DoFs 22, 19, 16, 14, 17 will be constrained. In this procedure the algebraic constraints do not depend on the enrichment func- tions and are equivalent to those one would have for the vector-value bases build upon the same list of scalar FEs. Therefore, no extension of the existing functionality to build algebraic constraints was necessary. This allows us to reuse the code written for the FESystem class, which can be used in deal.II to build a vector-valued FE from a collection of scalar-valued elements. Another remarkable beneﬁt of this approach is that existing code can be used to transfer the solution during h -adaptive reﬁnement from a coarse to a ﬁne mesh. The reason is that prolongation matrices for enriched elements are equal to their vector-valued counterparts under the condition that all child elements are also enriched. Yet another advantage of implementing a dedicated enriched ﬁnite element in deal.II library relates to the numerical integration of jump terms in the Kelly error indicator (see “Error estimator” section). Here, care needs to be taken in computing contributions to cell errors from faces with hanging nodes. Had the authors pursued an implementation of PUM where values and gradients of addi- tional basis functions are evaluated manually via the product rule in the course of numerical integration, a completely separate function would have to be implemented to integrate the jump terms in error indicators, which is not a straight-forward task. We believe that the implementation outlined above is general enough to allow PUM to be applied using deal.II to other partial diﬀerential equations such as crack propaga- Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 22 of 23 ab c Fig. 9 h -adaptive mesh reﬁnement and shape functions associated with the central node on the domain [0, 1] for the standard and enriched element. a Mesh. b Bilinear. c Bilinear enriched with exp(− |x|) tion in continuum mechanics. Figure 9 depicts an example of enriched and non-enriched shape functions for the case of h-adaptive reﬁnement with hanging nodes in two dimen- sions. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional aﬃliations. Received: 18 September 2017 Accepted: 29 November 2017 References 1. Maday Y. h–p ﬁnite element approximation for full-potential electronic structure calculations. Chin Ann Math Ser B. 2014;35(1):1–24. 2. Davydov D, Young T, Steinmann P. On the adaptive ﬁnite element analysis of the Kohn–Sham equations: methods, algorithms, and implementation. Int J Numer Methods Eng. 2016;106(11):863–88. https://doi.org/10.1002/nme. 3. Motamarri P, Nowak MR, Leiter K, Knap J, Gavini V. Higher-order adaptive ﬁnite-element methods for Kohn–Sham density functional theory. J Comput Phys. 2012;253(15):308–43. 4. Fang J, Gao X, Zhou A. A Kohn–Sham equation solver based on hexahedral ﬁnite elements. J Comput Phys. 2012;231(8):3166–80. 5. Pask JE, Sterne PA. Finite element methods in ab initio electronic structure calculations. Model Simul Mater Sci Eng. 2005;13(3):71–96. https://doi.org/10.1088/0965-0393/13/3/R01. 6. Zhang D, Shen L, Zhou A, Gong X-G. Finite element method for solving Kohn–Sham equations based on self-adaptive tetrahedral mesh. Phys Lett A. 2008;372(30):5071–6. https://doi.org/10.1016/j.physleta.2008.05.075. 7. Bylaska EJ, Holst M, Weare JH. Adaptive ﬁnite element method for solving the exact Kohn–Sham equation of density functional theory. J Chem Theory Comput. 2009;5(4):937–48. https://doi.org/10.1021/ct800350j. 8. Bao G, Hu G, Liu D. Numerical solution of the Kohn–Sham equation by ﬁnite element methods with an adaptive mesh redistribution technique. J Sci Comput. 2012;55(2):372–91. https://doi.org/10.1007/s10915-012-9636-1. 9. Sukumar N, Pask JE. Classical and enriched ﬁnite element formulations for Bloch-periodic boundary conditions. Int J Numer Methods Eng. 2009;77(8):1121–38. 10. Fattebert JL, Hornung RD, Wissink AM. Finite element approach for density functional theory calculations on locally-reﬁned meshes. J Comput Phys. 2007;223(2):759–73. https://doi.org/10.1016/j.jcp.2006.10.013. 11. White SR, Wilkins JW, Teter MP. Finite-element method for electronic structure. Phys Rev B. 1989;39:5819–33. https:// doi.org/10.1103/PhysRevB.39.5819. 12. Cimrman R, Novák M, Kolman R, Tuma ˘ M, Vackáˇr J. Isogeometric analysis in electronic structure calculations. Math Comput Simul. 2016;145:125–35. https://doi.org/10.1016/j.matcom.2016.05.011. 13. Chen H, Dai X, Gong X, He L, Zhou A. Adaptive ﬁnite element approximations for Kohn–Sham models. Multiscale Model Simul. 2014;12(4):1828–69. https://doi.org/10.1137/130916096. 14. Chen H, He L, Zhou A. Finite element approximations of nonlinear eigenvalue problems in quantum physics. Comput Methods Appl Mech Eng. 2011;200(21):1846–65. 15. Melenk JM, Babuška I. The partition of unity ﬁnite element method: basic theory and applications. Comput Methods Appl Mech Eng. 1996;139(1–4):289–314. https://doi.org/10.1016/S0045-7825(96)01087-0. 16. Babuška I, Melenk JM. The partition of unity method. Int J Numer Methods Eng. 1997;40(4):727–58. 17. Belytschko T, Black T. Elastic crack growth in ﬁnite elements with minimal remeshing. Int J Numer Methods Eng. 1999;45(5):601–20. 18. Simone A. Partition of unity-based discontinuous ﬁnite elements: gfem, pufem, xfem. Revue Européenne de Génie Civil. 2007;11(7–8):1045–68. 19. Belytschko T, Gracie R, Ventura G. A review of extended/generalized ﬁnite element methods for material modeling. Model Simul Mater Sci Eng. 2009;17(4):043001. Davydov et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:7 Page 23 of 23 20. Fries T-P, Belytschko T. The extended/generalized ﬁnite element method: an overview of the method and its applications. Int J Numer Methods Eng. 2010;84(3):253–304. 21. Mitchell WF, McClain MA. A comparison of hp-adaptive strategies for elliptic partial diﬀerential equations. ACM Trans Math Softw. 2014;41(1):2. 22. Houston P, Süli E. A note on the design of hp-adaptive ﬁnite element methods for elliptic partial diﬀerential equations. Comput Methods Appl Mech Eng. 2005;194(2):229–43. 23. Hartmann R, Houston P. Error estimation and adaptive mesh reﬁnement for aerodynamic ﬂows. In: ADIGMA—a European initiative on the development of adaptive higher-order variational methods for aerospace applications. Berlin: Springer; 2010. p. 339–53. 24. Mavriplis C. Adaptive mesh strategies for the spectral element method. Comput Methods Appl Mech Eng. 1994;116(1):77–86. 25. Eibner T, Melenk JM. An adaptive strategy for hp-fem based on testing for analyticity. Comput Mech. 2007;39(5):575–95. 26. Fankhauser T, Wihler TP, Wirz M. The hp-adaptive fem based on continuous Sobolev embeddings: isotropic reﬁnements. Comput Math Appl. 2014;67(4):854–68. 27. Bangerth W, Kayser-herold O. Data structures and requirements for hp ﬁnite element software. ACM Trans Math Softw. 2009;36(1):4. 28. Hohenberg P, Kohn W. Inhomogeneous electron gas. Phys Rev. 1964;136(3B):864–71. 29. Kohn W, Sham LJ. Self-consistent equations including exchange and correlation eﬀects. Phys Rev. 1965;140(4A):1133–8. 30. Arndt D, Bangerth W, Davydov D, Heister T, Heltai L, Kronbichler M, Maier M, Pelteret J-P, Turcksin B, Wells D. The deal.II library, version 8.5. J Numer Math. 2017. https://doi.org/10.1515/jnma-2017-0058. 31. Griﬃths DJ. Introduction to quantum mechanics. 2nd ed. London: Pearson; 2005. 32. Kanungo B, Gavini V. Large-scale all-electron density functional theory calculations using an enriched ﬁnite-element basis. Phys Rev B. 2017;95(3):035112. 33. Verfürth R. A review of a posteriori error estimation and adaptive mesh-reﬁnement techniques. New York: Wiley; 34. Larson MG. A posteriori and a priori error analysis for ﬁnite element approximations of self-adjoint elliptic eigenvalue problems. SIAM J Numer Anal. 2000;38(2):608–25. 35. Heuveline V, Rannacher R. A posteriori error control for ﬁnite element approximations of elliptic eigenvalue problems. Adv Comput Math. 2001;15(1–4):107–38. 36. Durán RG, Padra C, Rodríguez R. A posteriori error estimates for the ﬁnite element approximation of eigenvalue problems. Math Models Methods Appl Sci. 2003;13(08):1219–29. 37. Mao D, Shen L, Zhou A. Adaptive ﬁnite element algorithms for eigenvalue problems based on local averaging type a posteriori error estimates. Adv Comput Math. 2006;25(1–3):135–60. 38. Dai X, Xu J, Zhou A. Convergence and optimal complexity of adaptive ﬁnite element eigenvalue computations. Numerische Mathematik. 2008;110(3):313–55. https://doi.org/10.1007/s00211-008-0169-3. 39. Garau EM, Morin P, Zuppa C. Convergence of adaptive ﬁnite element methods for eigenvalue problems. Math Models Methods Appl Sci. 2009;19(05):721–47. 40. Houston P, Senior B, Süli E. Sobolev regularity estimation for hp-adaptive ﬁnite element methods. In: Brezzi F, Buﬀa A, Corsaro S, Murli A, editors. Numerical mathematics and advanced applications. Berlin: Springer; 2003. p. 631–56. 41. Melenk JM, Wohlmuth BI. On residual-based a posteriori error estimation in hp-fem. Adv Comput Math. 2001;15(1–4):311–31. https://doi.org/10.1023/A:1014268310921. 42. Heuveline V, Rannacher R. Duality-based adaptivity in the hp-ﬁnite element method. J Numer Math Jnma. 2003;11(2):95–113. 43. Giani S, Grubišic´ L, Ovall JS. Benchmark results for testing adaptive ﬁnite element eigenvalue procedures. Appl Numer Math. 2012;62(2):121–40. 44. Bangerth W. The deal.II library tutorial step 27 (version 8.3). 2007. https://www.dealii.org/8.3.0/doxygen/deal.II/ step_27.html. Accessed Jan 2016. 45. Schauer V, Linder C. All-electron Kohn–Sham density functional theory on hierarchic ﬁnite element spaces. J Comput Phys. 2013;250:644–64. https://doi.org/10.1016/j.jcp.2013.04.020. 46. Ciarlet PG. The ﬁnite element method for elliptic problems. Amsterdam: Elsevier; 1978. 47. Ming P, Shi Z-C. Quadrilateral mesh revisited. Comput Methods Appl Mech Eng. 2002;191(49):5671–82. 48. Veeser A, Verfürth R. Poincaré constants for ﬁnite element stars. IMA J Numer Anal. 2012;32(1):30–47. https://doi.org/ 10.1093/imanum/drr011. 49. Verfürth R. Error estimates for some quasi-interpolation operators. ESAIM Math Model Numer Anal. 1999;33(04):695–713. 50. Ainsworth M, Oden JT. A posteriori error estimation in ﬁnite element analysis. Comput Methods Appl Mech Eng. 1997;142:1–88.
"Advanced Modeling and Simulation in Engineering Sciences" – Springer Journals
Published: Dec 1, 2017
You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
Copy and paste the desired citation format or use the link below to download a file formatted for EndNote
Access the full text.
Sign up today, get DeepDyve free for 14 days.
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.