Computational Mechanics (2025) 76:205–226 https://doi.org/10.1007/s00466-025-02596-3 ORIG INAL PAPER Higher-order accurate mass lumping for explicit isogeometric methods based on approximate dual basis functions R. R. Hiemstra1 · T. H. Nguyen2 · S. Eisenträger3 ·W. Dornisch4 · D. Schillinger5 Received: 8 May 2024 / Accepted: 26 December 2024 / Published online: 25 January 2025 © The Author(s) 2025 Abstract This paper introduces a mathematical framework for explicit structural dynamics, employing approximate dual functionals and row-summass lumping.We demonstrate that the Petrov–Galerkinmethod developed in our previous work—utilizing row- sum mass lumping—can be interpreted as a Galerkin method with a customized higher-order accurate mass matrix. Unlike prior work, our method correctly incorporates Dirichlet boundary conditions, while preserving optimal spatial accuracy. The mathematical analysis is substantiated by spectral analysis, two-dimensional linear benchmarks illustrating robustness under mesh distortion and a three-dimensional geometrically non-linear forced vibration analysis. The obtained results reveal that our approach achieves accuracy and robustness comparable to a traditional Galerkin method employing the consistent mass formulation, while retaining the explicit nature of the lumped mass formulation. Keywords Isogeometric analysis · Explicit dynamics ·Mass lumping ·Customized mass matrix ·Dual functionals · Spectrum analysis B R. R. Hiemstra r.r.hiemstra@tue.nl T. H. Nguyen hoa.nguyen@uib.no S. Eisenträger sascha.eisentraeger@ovgu.de W. Dornisch wolfgang.dornisch@rptu.de D. Schillinger dominik.schillinger@tu-darmstadt.de 1 Department of Mechanical Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands 2 Faculty of Mathematics and Natural Sciences, Geophysical Institute, Bergen Offshore Wind Center, University of Bergen, Bergen, Norway 3 Faculty of Mechanical Engineering, Institute of Mechanics, Chair of Computational Mechanics, Otto-von-Guericke-Universität Magdeburg, Magdeburg, Germany 4 Department of Mechanical and Process Engineering, Chair of Applied Mechanics, Rheinland-Pfälzische Technische Universität Kaiserslautern-Landau, Kaiserslautern, Germany 5 Department of Civil and Environmental Engineering, Institute for Mechanics, Computational Mechanics Group, Technische Universität Darmstadt, Darmstadt, Germany 1 Introduction Explicit methods for structural dynamics advance or update the solution at each time increment according to an evolution equation represented by the matrix problem M ün+1 = F(un, u̇n, ün) . (1) Here,M ∈ R N×N is themassmatrix (which is a sparse, sym- metric and positive definite—SPD—matrix), ün+1 ∈ R N is a vector of unknown accelerations at time tn+1, and F ∈ R N encapsulates the interior and exterior forces, which are a function of the known displacement, un , velocity, u̇n , and acceleration, ün , at time tn . State-of-the-art commercial solvers such as LS-Dyna1, RADIOSS2, or PAM-CRASH3 are all based on first order finite elements (linear shape functions, p = 1) [2–4] in con- junction with explicit second-order accurate time integration methods such as the central difference method or variations thereof [4]. In typical applications such as the analysis of vehicle crashworthiness and the simulation of metal forming 1 https://www.ansys.com/products/structures/ansys-ls-dyna. 2 https://altair.com/radioss. 3 https://www.esi-group.com/pam-crash. 123 http://crossmark.crossref.org/dialog/?doi=10.1007/s00466-025-02596-3&domain=pdf https://orcid.org/0000-0002-0179-9606 https://orcid.org/0000-0003-0519-5786 https://orcid.org/0000-0001-8774-9732 https://orcid.org/0000-0003-0690-8063 https://orcid.org/0000-0002-9068-6311 https://www.ansys.com/products/structures/ansys-ls-dyna https://altair.com/radioss https://www.esi-group.com/pam-crash 206 Computational Mechanics (2025) 76:205–226 or stamping, the internal force vector of the discrete sys- tem (1) involves shell elements, nonlinear material models, and contact with friction, for which low-order finite elements provide a well-developed technological setting. A key com- ponent is the replacement of the consistent mass matrix M with a diagonal (lumped) matrix, ML , obtained by “mass lumping”. It results in simple and highly computationally and memory efficient solution update scheme, while main- taining optimal spatial accuracy of the finite elementmethod. In addition, the lumped mass matrix reduces the highest fre- quencies of the discretization, thereby increasing the critical time step size. Therefore, these codes achieve an extreme level of computational and memory efficiency. Initiated in 2005, isogeometric analysis (IGA) aims at bridging the gap between computer aided geometric design and finite element analysis [5, 6]. The core idea of IGA is to use the same smooth and higher-order spline basis functions for the representation of both geometry and the approxima- tion of physics-based field solutions. For a more detailed introduction to IGA, we refer interested readers to the recent review articles [7–10] and the references cited therein. Due to the higher-order nature of splines, IGA is particularly attractive for higher-order accurate structural analysis. For elastostatic-type problems, higher-order IGA is superior to standard FEA in terms of per-degree-of-freedom accuracy, as shown, e.g., in [11]. Exploiting its higher interelement continuity, IGA (with consistent mass) also performs excep- tionally well in wave propagation problems [12]. A key property of higher-order IGA, already discussed in one of the first articles [13], is its well-behaved discrete spec- trum of eigenfrequencies and eigenmodes. The associated potential of IGA for efficient explicit dynamics, however, lies largely idle to this day. One of the main barriers is the lack of a diagonalization scheme for the mass matrix enabling a computationally andmemory efficient update, but alsomain- taining the higher-order spatial accuracy of the method [13, 14]. For instance, Hartmann and Benson reported in [14] that “increasing the degree on a fixed mesh size increases the cost without a commensurate increase in accuracy”. There- fore, studies on isogeometric explicit dynamics schemes that have adopted standard row-sum lumping to diagonalize the mass matrix have focused on lower-order quadratic spline discretizations [15–17]. One idea to bridge this gap is to apply additional cor- rector updates that restore higher-order accuracy [18] at low memory requirements, but lead to an (unacceptable) increase in computational costs. A recent study [19] explored ideas from preconditioning to develop certain customized mass matrices that approximate the consistent mass matrix while retaining some of the beneficial properties of the lumped mass. An extension of this approach to multi-patch and trimmed domains, including an algebraic outlier removal by means of an eigenvalue deflation technique, has been presented in [20]. However, the presented methods have not led to a higher-order accurate technology. For standard C0-continuous basis functions, a pathway to higher-order accurate explicit dynamics with diagonal mass matrices is the adoption of Lagrange basis functions collocated at the Gauss-Lobatto Legendre (GLL) nodes [12, 21–23]. Due to the interpolatory property of nodal basis functions, this strat- egy naturally yields a consistent and higher-order accurate diagonal mass matrix. Another idea is therefore to trans- fer this mechanism to isogeometric analysis via a change of basis in the context of Lagrange extraction [24, 25]. But the associated matrix operations, even when accelerated through special techniques, increase the computational cost and the memory requirements beyond an acceptable level. Another approach to achieve a diagonal mass matrix is the adoption of dual spline functions. This set of functions is defined bi-orthogonal (or “dual”) to a set of standard B- splines (or NURBS), such that the inner product of the two integrated over a domain yields a diagonal matrix. Hence, if a dual basis is used as the test space in a finite element formulation, we naturally obtain a consistent diagonal mass matrix [26]. In CAD research, this bi-orthogonality property has been well known for decades, and several approaches for the construction of the underlying dual basis exist, see, e.g., [27]. In the computational mechanics community, this prop- erty has been used for different purposes, e.g., for dualmortar methods [28–31] or for unlocking of Reissner-Mindlin shell formulations [32]. Dual spline basis functions still suffer from shortcomings when applied as discrete test functions in elastodynamics as envisioned in [26], to the extent that their use in explicit dynamics seems not straightforward (as recently summa- rized in [1]): Firstly, bi-orthogonality in general holds only on the parametric domain, but is in general lost under non- affine geometric mapping. That is, the corresponding matrix is diagonal on the parametric domain, but not diagonal in the physical domain. Secondly, dual spline functions with the same smoothness as their B-spline counterparts have global support on each patch, thus producing fully populated stiff- ness forms. Their support can be localized, but at the price of losing continuity, leading to discontinuous functions in the fully localized case. Both are prohibitive in a finite element context. And thirdly, dual basis functions are not interpo- latory at the boundaries. Therefore, the identification of a kinematically admissible set of test functions that allows the variationally consistent strong imposition of Dirichlet boundary conditions is not straightforward. In other words, a competitive higher-order technology that removes all current barriers has not been realized. In a recent paper [1], we presented an isogeometric Petrov–Galerkin method for explicit dynamics that over- comes the first two of the three challenges. Its key ingredient is a class of “approximate” dual spline functions, introduced 123 Computational Mechanics (2025) 76:205–226 207 in [33] and successfully applied within IGA for dual mortar- ing in [29, 31]. It only approximately satisfies the discrete bi-orthogonality property, but preserves all other proper- ties of the original B-spline basis, such as its smoothness, polynomial reproduction, and local support. Employing the approximate dual functions to discretize the test function space, we obtain a semi-discrete Petrov–Galerkin formula- tion in the format (1), whose mass matrix M is not truly diagonal, but a “close” approximation of a diagonal matrix. We showed that we can then apply standard row-sum mass lumping to diagonalize the approximately diagonal mass matrix without compromising higher-order spatial accuracy. A similar study was recently presented in [34], confirming our promising results. Another approach related to our work was presented in [35]. Here, an improved spatial accuracy was achieved by combining basis interpolation with mass lumping. Although, the proposed approach yields higher- order spatial accuracy, it is far from optimal and further research is required. In this paper, we build on and extend our results presented in [1],with the goal of providing a comprehensivemathemati- cal framework of approximate dual basis functions, including a consistent method to strongly incorporate Dirichlet bound- ary conditions. In Sect. 2, we set the stage by fixing a general linear second order model problem, which we discretize in space via the Galerkin method and in time via higher-order explicit time integration schemes.Weadditionally reviewdif- ferent classical ways of performing mass lumping. In Sect. 3, we focus on a univariate basis for splines that is approx- imately dual to B-splines in the L2-norm and discuss their properties. We demonstrate that the Petrov–Galerkin method with row-summass lumping proposed in [1] can also be inter- preted as a Galerkin method with a customized mass matrix. Based on this observation, we show that the customizedmass technique does not have a negative effect on the asymptotic accuracy of the method. Another important extension with respect to our work in [1] is the incorporation of Dirichlet boundary conditions. We show that polynomial accuracy is maintained under strong imposition of boundary conditions. In Sect. 4,we extend our developments to themultivariate set- ting, explain how the properties of the basis can be preserved under general non-affine mappings, and discuss key aspects of an efficient computer implementation. In Sect. 5, we verify accuracy and robustness of ourmethodology in the context of spectral analysis, a two-dimensional linear benchmark, and a three-dimensional geometrically non-linear forced vibration problem, demonstrating that our Petrov–Galerkin scheme indeed preserves higher-order accuracy in explicit dynamics simulations. In Sect. 6, we provide a summary and outline potential future research directions. 2 Preliminaries In this section, we present background and notation used in the remainder of this paper. First, we present a general linear second order model problem. The Galerkin method is used to discretize the equations in space and we briefly discuss explicit time integration. Subsequently, we discuss the consistent mass matrix and different classical ways of performing mass lumping. 2.1 Model problem Let � ⊂ R d be an open set with piecewise smooth boundary � = �g ∪ �h, �g ∩ �h = ∅. Consider (sufficiently smooth) prescribed boundary and initial data f : � × (0, T ) �→ R, (2a) g : �g × (0, T ) �→ R, (2b) h : �h × (0, T ) �→ R, (2c) u0, u̇0 : � �→ R. (2d) Let V ≡ H1(�). The space of trial solutions is assumed to be constant in time, Vg := { v ∈ V with v(x) = g(x) for x ∈ �g } . (2e) Here, theDirichlet data, g, is built directly into the trial-space, while the Neumann data, h, is incorporated weakly through natural imposition. We consider weak formulations of the following form: given f , g, h, u0 and u̇0, find u(t) ∈ Vg such that for all w ∈ V0 a(u, w) + b(ü, w) = l(w), (2f) b(u(0), w) = b(u0, w), (2g) b(u̇(0), w) = b(u̇0, w). (2h) Here, a : V × V �→ R is a positive semi-definite and symmetric bilinear form (definite on the subspace V0 × V0) and b : V×V �→ R is a positive definite symmetric bilinear form defined as: b(u, v) = ∫ � ρ u v d�, where ρ : � �→ R denotes the mass density of the medium. Furthermore, l : V �→ R is a time dependent linear form dependent on the prescribed forces and Neumann boundary data: l(w) =∫ � w f d� + ∫ �h w h d�. The model problem described here can represent different evolutionary processes, such as the time-dependent linear wave equation (d = 1, 2, 3), vibration of a string (d = 1), or vibration of a linear membrane (d = 2). 123 208 Computational Mechanics (2025) 76:205–226 2.2 Semi-discrete equations of motion We consider semi-discrete systems of equations obtained using the Galerkin method with a finite dimensional space Vh ⊂ V and subspaceVh 0 ⊂ V0, discretized in terms of tensor product spline functions. Let uh = vh + gh where vh ∈ Vh 0 and gh ∈ Vh approximately satisfies the Dirichlet boundary data g on �g . The Galerkin formulation reads: find vh(t) ∈ Vh 0 such that for all wh ∈ Vh 0 a(wh, vh) + b(wh, v̈h) = l(wh) − a(wh, gh) − b(wh, g̈h), (3a) b(wh, vh(0)) = b(wh, u0) − b(wh, gh(0)), (3b) b(wh, v̇h(0)) = b(wh, u̇0) − b(wh, ġh(0)). (3c) Consider finite dimensional representations in terms of basis functions Ni(x) ∈ Vh vh(x) = ∑ i∈η−ηg Ni(x) di(t), gh(x) = ∑ i∈ηg Ni(x) di(t), (4) where η is an index set of all functions in Vh and ηg is an index set of functions that intersect the boundary. We consider matrices M = [Mij] and K = [Kij], and vectors F = [Fi], D0 = [D0i], and Ḋ0 = [Ḋ0i] (i, j ∈ η − ηg) defined by Mij := b(Ni, Nj), (5a) Kij := a(Ni, Nj), (5b) Fi := l(Ni) − a(Ni, g h) − b(Ni, g̈ h), (5c) D0i := b(Ni, u0 − gh(0)), (5d) Ḋ0i := b(Ni, u̇0 − ġh(0)). (5e) Applying these definitions to (3) leads to the followingmatrix problem, M d̈(t) + K d(t) = F(t), t ∈ (0, T ), (6a) M d(0) = D0, (6b) M ḋ(0) = Ḋ0. (6c) This is a coupled system of ordinary differential equations (ODE’s) for the displacement coefficients d(t) := [di(t)], where ḋ(t) := [ḋi(t)] denote velocities and d̈(t) := [d̈i(t)] the accelerations. Theproperties of the stiffnessmatrix,K, and the consistent mass matrix, M, follow directly from the bilinear forms a and b, respectively, and the finite element basis functions ( Ni, i ∈ η − ηg ). In particular, the consistent mass matrix M and the stiffness matrix K are both positive definite and symmetric and have a sparse, banded structure due to the local support of the basis functions. In addition, the positivity of B-splines as trial and test functions leads directly to non- negative entries in the mass matrix. 2.3 Explicit time steppingmethods Let dk, ḋk, and d̈k denote certain finite difference approxi- mations of d(tk), ḋ(tk) and d̈(tk), respectively, that describe the state at time tk ∈ (0, T ), and let Fk = F(tk). Given the state at a particular time tn ∈ [0, T ), explicit methods deter- mine the state at a later time tn+1 ∈ [0, T ] by replacing the update in (6a) with the matrix problem M d̈n+1 = Fn+1 − K dn . (7) The stiffnessmatrixK need not be explicitly formed. Instead, its action on a vector dn may be implemented in a matrix- free context [36]. A solution update using the consistentmass matrixM still requiresmatrix backsolve scalingwithO(N 2), where N is the leading dimension of the system. 2.4 Efficient update throughmass lumping Mass lumping refers to a procedurewhere the consistentmass matrix is somehow approximated by a diagonal matrix of the same dimensions to enable fast explicit solution updates, that is, updates that do not require matrix inversion and scale O(N ). The three important procedures are (see, e.g., [37, Section 12.2.4]): 1. The row-sum method: M lumped ii := ∑ j Mij. 2. Diagonal scaling: M lumped ii = cMii with c chosen such that ∑ j M lumped jj = ∑ i ∑ j Mij = ∫ � ρd�. 3. Evaluation of Mij using a nodal quadrature rule corre- sponding to nodal shape functions. The first two techniques have proven robust and efficient mainly for linear finite elements (p = 1) [38], while the third procedure has had success in higher-order spectral element methods [22, 23]. In this work, we use row-sum lumping and show that higher-order accuracy can be maintained for a spe- cial class of lumped mass matrices in which the test and trial functions are not the same (Petrov–Galerkin formulation). 3 Approximate L2 dual bases In this section, we introduce a basis for splines that is approx- imately dual to B-splines in the L2-norm. In contrast to the exact dual functions, these functions are locally supported, 123 Computational Mechanics (2025) 76:205–226 209 which makes them suitable as test functions in a Petrov– Galerkin method. We demonstrate that the Petrov–Galerkin methodwithmass lumping is equivalent to aGalerkinmethod with a higher-order accurate approximate mass matrix that has an explicit sparse inverse. We mathematically show that our approximate mass technique maintains higher-order spatial accuracy in the context of explicit dynamics. An important extension with respect to previous work [1, 26] is the consistent incorporation of Dirichlet boundary condi- tions.We show that polynomial accuracy ismaintained under strong imposition of boundary conditions. 3.1 Univariate B-splines A spline is a piecewise polynomial that is characterized by the polynomial degree of its segments and the regularity pre- scribed at their interfaces. A convenient basis in which to represent polynomial splines is given byB-splines [39]. Con- sider a partitioning of the univariate interval � = [a, b] into δ elements a = x0 < x1 < . . . < xδ = b. (8) With every internal breakpoint, xk , we may associate an integer, rk , that prescribes the smoothness between the poly- nomial pieces. Then, given the knot-multiplicity, ( p − rk )δ−1 k=1, we can define the knot vector, � = ( ti ) n+p+1 i=1 := { x0, . . . , x0︸ ︷︷ ︸ p+1 , . . . , xk−1, . . . , xk−1︸ ︷︷ ︸ p−rk−1 , xk, . . . , xk︸ ︷︷ ︸ p−rk , . . . , xδ, . . . , xδ︸ ︷︷ ︸ p+1 } . (9) With aknot vector in hand,B-splines are stably and efficiently computed using the Cox-de Boor recursion, Bi,0(x) = { 1 if x ∈ [ti , ti+1) 0 otherwise , (10a) Bi,p(x) = x − ti ti+p−ti Bi,p−1(x)+ ti+p+1 − x ti+p+1 − ti+1 Bi+1,p−1(x) , (10b) where, by definition, 0/0 = 0. Let r = ( 0 ≤ rk ≤ p−1, k = 1, . . . , δ −1 ). The space of polynomial splines of degree p and dimension n with rk continuous derivatives at breakpoint xk is defined as S p r (�) := span ( Bi,p(x), i = 1, . . . , n ) . (11) To shorten notation we drop the dependence on the domain and simply write Sp r . B-splines have important mathematical properties, many of which are useful in design as well as in analysis. B-spline basis functions of degree p may have up to p−1 continuous derivatives, they form a positive partition of unity, and have local support of up to p + 1 elements. 3.2 L2 dual basis Let 〈·, ·〉 : L2(�) × L2(�) �→ R denote the standard L2 inner product of two square integrable functions. A set of functions ( B∗ i , i = 1, . . . , n ) is called a dual basis (in the L2-norm) to the set ( Bi , i = 1, . . . , n ) if the following property holds 〈B∗ i , Bj 〉 = δi j , i, j ∈ 1, . . . , n. (12) Here, δi j denotes the Kronecker delta, which is one if i = j and zero otherwise. The one-to-one correspondence implies that the dual functions are linearly independent if and only if the functions ( Bi , i = 1, . . . , n ) are linearly independent. A dual basis is not unique. There is freedom in choosing a space spanned by the functions, there are different ways to construct them, and their properties, such as smoothness and support, may vary. If the dual basis functions span the space of splines Sp r , then they are uniquely defined and can be computed as (duality in (12) may be checked by direct substitution) B∗ i (x) = n∑ j=1 ( G−1 ) i j B j (x), where Gi j = 〈Bi , Bj 〉. (13) A Galerkin discretization employing L2 dual functions as test functions and B-splines as trial functions yields the iden- tity matrix as mass matrix. Although this property is highly sought after in explicit dynamics, they are unsuitable as a basis for the test-space. Because the inverse of the Gramian matrix, G−1, is a dense matrix, the L2 dual functions have global support. This leads to poor scaling and expensive for- mation and assembly of the interior force vector in explicit dynamics computations, because each entry in the vector depends on all the quadrature points. Finally, we list a well known property of the dual basis used in the context of L2-projection (Proposition 3.1). In the next subsection we introduce a set of functions that are (in some sense) approximately dual to B-splines, yet maintain local support, and satisfy an approximation property similar to the following. Theorem 3.1 (L2-projection) Let f ∈ L2(�). The function u(x) = ∑n i=1〈 f , B∗ i 〉 Bi (x) is the best spline approximation to f in the L2-norm. In other words, it is the minimizer of the convex minimization problem 123 210 Computational Mechanics (2025) 76:205–226 argmin u∈Sp r (�) 1 2 〈u, u〉 − 〈u, f 〉 Proof This follows by direct application of the dual func- tions. � 3.3 L2 approximate dual basis We loosen the restriction of duality and search for a set of functions ( B̂i , i = 1, . . . , n )with the following properties: 1. Basis for the spline space Sp �r . 2. Local compact support. 3. Approximate L2 duality such that ∑n i=1〈 f , B̂i 〉 Bi (x) preserves polynomials of degree p. The construction relies on symmetric positive definite (SPD) matrices. With SPD n ⊂ R n×n we denote the set of all such matrices of n rows and columns. We recall that SPD matrices are invertible and their inverse is also SPD. Definition 3.2 (Approximate dual basis) Let P p(�) ⊂ S p r (�) denote the set of polynomials of degree < p + 1 and let S ∈ SPD n . We call the set of functions B̂i (x) = n∑ j=1 Si j B j (x), i = 1, 2, . . . , n, (14) an approximate dual basis to the set of B-splines if 〈 f , B̂i 〉 = 〈 f , B∗ i 〉 ∀ f ∈ P p(�). (15) The identity (15) implies approximate dualitywithB-splines. Therefore, the functions B̂i , (i = 1, . . . , n) will be referred to as approximate dual functions. Their collection forms a basis for Sp r (since S is invertible) and is called an approx- imate dual basis. We note that the representation is not unique.Aparticularly useful approximate dual basis has been developed in [33]. The approximate dual functions presented therein have minimum compact support, in the sense that there exists no other basis that satisfies the approximate dual- ity with smaller supports. In the remainder of this paper, we shall base our construction on that approach. SPD matrices induce a discrete inner product on finite dimensional spaces. Definition 3.3 (Discrete inner product on splines) Let u, v ∈ S p r (�) be expanded in terms of B-splines and coefficients as u(x) = ∑n i=1 ui Bi (x) and v(x) = ∑n i=1 vi Bi (x), respec- tively. An SPD-matrix Ĝ ∈ SPD n induces a discrete inner product 〈u, v〉Ĝ : S p r × S p r �→ R defined by 〈u, v〉Ĝ := n∑ i=1 n∑ j=1 ui Ĝi j v j . (16) We note that 〈u, v〉G (where G is the Gramian matrix) coincides with the usual L2 inner product on splines, that is, 〈u, v〉G = 〈u, v〉 for spline functions u and v. Similar to (13) the approximate dual functions in (14) are expressed as linear combinations ofB-splines via amatrix vector product involv- ing an SPD matrix. Hence, analogous to Theorem (3.1), the associated quasi-interpolant may be associated to a mini- mization problem. Theorem 3.4 (Quasi L2-projection ) Let f ∈ L2(�) and Ĝ = S−1. The spline function u(x) = ∑n i=1〈 f , B̂i 〉 Bi (x) is the unique solution to the convex minimization problem argmin u∈Sp r (�) 1 2 〈u, u〉Ĝ − 〈u, f 〉. (17) Furthermore, u = f whenever f ∈ P p. Proof Let u(x) = ∑n i=1 ui Bi (x). The optimization problem is equivalent to the following set of linear equations n∑ j=1 Ĝi j u j = 〈 f , Bi 〉, i = 1, . . . , n. (18) Its unique solution is ui = ∑n j=1 Si j 〈 f , Bj 〉 = 〈 f , B̂i 〉. If f ∈ P p, then 〈 f , B̂i 〉 = 〈 f , B∗ i 〉. Here, ui = 〈 f , B∗ i 〉 is the best approximation in the L2-norm. It reproduces all splines in the space, which includes all polynomials f ∈ P p. 3.4 Treatment of Dirichlet boundary conditions By associating the quasi-interpolant with the solution to a certain convex minimization problem it becomes straight- forward to consider linear constraints, such as Dirichlet boundary conditions. Without loss of generality we focus on the application of a single homogeneous boundary condition on the left boundary of the univariate domain. Theorem 3.5 (Quasi L2-projection with a homogeneous boundary condition) Suppose f ∈ H1(�) satisfies f (a) = 0. The spline function u(x) = ∑n i=2〈 f , B̂i 〉 Bi (x), where ⎡ ⎢ ⎣ B̂2(x) ... B̂n(x) ⎤ ⎥ ⎦ = ⎡ ⎢ ⎣ Ĝ22 · · · Ĝ2n ... . . . ... Ĝn2 · · · Ĝnn ⎤ ⎥ ⎦ −1 ⎡ ⎢ ⎣ B2(x) ... Bn(x) ⎤ ⎥ ⎦ , (19) is the unique solution to the convex minimization problem argmin u∈Sp r 1 2 〈u, u〉Ĝ − 〈u, f 〉 subject to u(a) = 0. (20) Furthermore, if f ∈ P p, then u = f . 123 Computational Mechanics (2025) 76:205–226 211 Proof In matrix-notation we may write (making use of u1 = u(a) = 0) argmin ui∈R 1 2 n∑ i=2 n∑ j=2 ui Ĝi j u j − n∑ i=2 ui 〈 f , Bi 〉, (21) which is equivalent to the matrix problem n∑ j=2 Ĝi j u j = 〈 f , Bi 〉, i = 2, . . . , n. Its solution can again be presented in terms of approximate dual functions: ui = 〈 f , B̂i 〉, where the linear combination is governed by (19). The constraint does not alter the energy in (20), which is obvious from the expression in (21). Hence, a solution of (20) is a solution to (17) if one considers only those func- tions that are exactly reproduced and satisfy the constraint: polynomials f ∈ P p that satisfy f (a) = 0. Hence, the quasi-interpolant preserves all polynomials that satisfy the boundary condition. We note that it is not necessary to compute the inverse of the submatrix of Ĝ directly. One can make efficient use of theWoodburymatrix identity [40, Section 2.1.4] to eliminate any unnecessary computations of inverse matrices ( Ĝ + UV T )−1 = S − SU ( I2 + V T SU )−1 V T S. (22) In our case, the matrix Ĝ +UV T is a rank two update of the original matrix ⎡ ⎢⎢⎢ ⎣ Ĝ11 0 · · · 0 0 Ĝ22 · · · Ĝ2n ... ... . . . ... 0 Ĝn2 · · · Ĝnn ⎤ ⎥⎥⎥ ⎦ ︸ ︷︷ ︸ Ĝ+UV T = ⎡ ⎢⎢⎢ ⎣ Ĝ11 Ĝ12 · · · Ĝ1n Ĝ21 Ĝ22 · · · Ĝ2n ... ... . . . ... Ĝn1 Ĝn2 · · · Ĝnn ⎤ ⎥⎥⎥ ⎦ ︸ ︷︷ ︸ Ĝ + ⎡ ⎢⎢⎢ ⎣ 0 −1 Ĝ21 0 ... ... Ĝn1 0 ⎤ ⎥⎥⎥ ⎦ ︸ ︷︷ ︸ U [−1 0 · · · 0 0 Ĝ21 · · · Ĝn1 ] ︸ ︷︷ ︸ V T . 3.5 Perspective 1: Galerkinmethod with a higher-order approximatemass matrix In the previous subsection, we developed an approximate dual basis and a quasi L2-projector that preserves polynomi- als with or without strong enforcement of Dirichlet boundary conditions. The quasi-projector is associated with a sparse SPD-matrix S, and its (dense) inverse, Ĝ = S−1, may be interpreted as a mass matrix that in some sense approximates the consistent mass matrix G. We briefly summarize how the approximate mass matrix may be used in explicit dynamics. The following equations concern the one-dimensional case with unit density. However, many properties directly carry over to the multi-dimensional setting using Kronecker prod- ucts and the techniques presented in the next section. Let f ∈ R n with fi = 〈Bi , f 〉, i = 1, . . . , n for some f ∈ L2(�). In explicit dynamics, an equation of the follow- ing form is solved at each time step G u = f , (Galerkin method with consistent mass) (23a) Ĝ û = f . (Galerkin method with higher- order approximate mass) (23b) While the consistent mass requires a backsolve to obtain the solution to (23a), we solve (23b) by a matrix–vector multi- plication with the explicit sparse inverse of the approximate mass matrix, that is, û = S f . As shown in the previous sub- section, this quasi L2 projection preserves all polynomials up to degree p, which implies that higher-order spatial accuracy is maintained in an explicit analysis context. 3.6 Perspective 2: Petrov–Galerkinmethod with a lumpedmass matrix We briefly discuss the equivalence with row-summass lump- ing in the Petrov–Galerkinmethod introduced in [1]. One can use approximate dual functions directly as test functions in a Petrov–Galerkin method: f̃i := 〈B̂i , f 〉, (24a) G̃i j := 〈B̂i , Bj 〉 = n∑ k=1 Sik Gkj . (24b) In an explicit dynamics context, we would solve at each time step an equation of the form G̃ v = f̃ . (Petrov–Galerkin method with consistent mass) (25a) Since the approximate dual functions span the space of splines, a backsolve would yield the same solution as the Galerkin method, that is, v = u, where u is the solution to (23a). Nowconsider the linear equations in the context of lumped mass 123 212 Computational Mechanics (2025) 76:205–226 G̃ lumped ṽ = f̃ . (Petrov–Galerkin method with lumped mass) (25b) Lemma 3.6 Each row of matrix G̃ sums to one, that is,∑n j=1 G̃i j = 1, for each i = 1, . . . , n. Proof We have ∑n j=1 G̃i j = ∑n j=1〈B̂i , Bj 〉. Partition of unity implies that n∑ j=1 〈B̂i , Bj 〉 = 〈B̂i , n∑ j=1 Bj 〉 = 〈B̂i , 1〉. Next, consider reproduction of constants, ∑ i 〈B̂i , 1〉 Bi (x) = 1. Linear independence of B-splines and partition of unity implies that 〈B̂i , 1〉 = 1, for each i . Consequently, ∑n j=1 G̃i j = 〈B̂i , 1〉 = 1 for each i = 1, . . . , n. Corollary 3.7 Row-sum lumping matrix G̃ yields the identity matrix, that is, G̃ lumped i j = δi j , i = 1, . . . , n. Since f̃ = S f and G̃ lumped = I d, the solution to (25b) is obtained as ṽ = S f , which is precisely the solution to (23b): ṽ = û.Hence, the perspective of thePetrov–Galerkinmethod in combination with row-sum mass lumping is equivalent to a standard Galerkin method with a higher-order accurate approximation of the mass matrix Ĝ, which happens to have an explicit sparse inverse denoted by S. A similar viewpoint was presented recently in [34]. 4 Application in higher-order accurate explicit dynamics The approximate dual basis introduced thus far maintains its properties, in particular polynomial reproduction, under affinemappings only. Here, we extend our approach to incor- porate non-linear geometric mappings as well. We derive the equations from the perspective of the Petrov–Galerkin method, where we employ B-splines as trial functions and (suitably weighted) approximate dual functions as test func- tions. Note that our approach is non-isoparametric because we allow geometricmappings to be defined by, e.g., NURBS. 4.1 Extension to arbitrary isoparametric mappings Let : �̂ → � denote a mapping from the unit-square to the physical domain and let F(x) = ∇ (x) denote the Jacobian matrix evaluated at a point x ∈ �̂. As usual, φ is an invertible mapping such that its Jacobian determinant is a positive function det (F) > 0.We let c◦ = det (F) ·ρ ◦ , which is a positive function by virtue of det (F) > 0 and ρ > 0. We proceed as in [1] and define trial and test functions as follows Ni ◦ := Bi1(x1) · Bi2(x2), i ∈ η, (trial func.) (26a) N̂i ◦ := B̂i1(x1) · B̂i2(x2) c ◦ , i ∈ η. (test func.) (26b) The weighting in the test functions plays a special role. Importantly, the resulting consistent mass matrix becomes invariant under geometric transformations. This is easy to see as follows M̃ij = b(N̂i, Nj) = ∫ � ρ N̂i Nj d� = ∫ �̂ ρ N̂i Nj det (F)d�̂. The determinant of the Jacobian and the density of the medium cancel with the weighting function, c = det (F) ·ρ, in the denominator of the test functions. Notably, the consis- tent mass matrix attains the following simple structure that does not depend on the geometric mapping M̃ij = ∫ 1 0 B̂i2(x2)Bj2(x2) dx 2 ∫ 1 0 B̂i1(x1)Bj1(x1)dx 1 (27a) and can be factorized as M̃ = S2G2 ⊗ S1G1. The entries of the stiffness matrix, the right hand side vector, and relevant boundary contributions become K̃ij = a(N̂i, Nj) = ∑ k Sik a(Nk/c, Nj), (27b) F̃i = l(N̂i) − a(N̂i, g h) − b(N̂i, g̈ h) = ∑ k Sik ( l(Nk/c) − a(Nk/c, g h) − b(Nk/c, g̈ h) ) , (27c) D̃0i := b(N̂i, u0 − gh(0)) = ∑ k Sik b(Nk/c, u0 − gh(0)), (27d) ˙̃D0i := b(N̂i, u̇0 − ġh(0)) = ∑ k Sik b(Nk/c, u̇0 − ġh(0)), (27e) where S = S2 ⊗ S1 is a Kronecker product matrix. Here, the subscripts refer to the component direction. In summary, the above representation of the mass matrix does not depend on the geometric mapping and has a Kro- 123 Computational Mechanics (2025) 76:205–226 213 necker structure. Hence, the theory of the previous section applies and can be readily extended to the multi-dimensional setting. Remark 4.1 Instead ofweighing test functionswith the scalar function, c, a weighted L2 inner product can be used to remove the dependency on the geometric mapping, see, e.g., [31] in the context of the mortar method. The resulting equa- tions are, however, identical. 4.2 Maintaining higher-order accuracy in explicit dynamics Prior to the imposition of Dirichlet boundary conditions, the semi-discrete equations associated with the explicit Petrov– Galerkin method are M̃ d̈n+1 = F̃n+1 − K̃ dn . (Petrov–Galerkin method with consistent mass) Applying row-sum mass lumping (M̃lumped = Id), the equa- tions simplify to d̈n+1 = F̃n+1 − K̃ dn . In the following we let Kij = a(Ni/c, Nj), (28) Fi = l(Ni/c) − a(Ni/c, g h) − b(Ni/c, g̈ h). (29) Hence, K̃ = SK and F̃ = S F, where S = S2 ⊗ S1 is a Kronecker product matrix. Using M̂ = S−1 = S−1 2 ⊗ S−1 1 = Ĝ2 ⊗ Ĝ1 we have: M̂ d̈n+1 = Fn+1 − K dn, (Galerkin method with higher-order approximate mass) which is efficiently solved using the matrix–vector product d̈n+1 = S2 ⊗ S1 (Fn+1 − K dn) . (30) This is a direct extension of the L2 quasi-projection (Theorem 3.4) to the two-dimensional setting using the Kro- necker product structure. The approximation properties in one dimension directly carry over to the multi-dimensional setting via theKronecker product. Hence, the approachmain- tains polynomial accuracy under mass lumping. 4.3 Imposition of Dirichlet boundary conditions Dirichlet boundary conditions are applied to the equations after performingmass lumping.We assume that the Dirichlet data are imposed separately for each side of the rectangular B-spline patch.Consequently, themassmatrix under strongly imposed homogeneous Dirichlet boundary conditions main- tains a factorization M̂ = Ĝ2 ⊗ Ĝ1, where homogeneous Dirichlet conditions are imposed on each of the univariate matrices Ĝk, k = 1, 2. If M̂, K and F denote the matri- ces and vectors obtained after strong imposition of Dirichlet data, then the final semi-discrete equations are d̈n+1 = S2 ⊗ S1 ( Fn+1 − K dn ) . (31) Again, the approximation properties of Theorem 3.5 carry directly over from the one-dimensional setting to the two- dimensional setting via the Kronecker product. Hence, the approachmaintains polynomial precision, evenwith strongly imposed Dirichlet data. 4.4 Implementational aspects Weemphasize the following points with regards to the imple- mentation: • The Kronecker structure of S2 ⊗ S1 facilitates efficient matrix storage that scales as∝ 2(2p+1) ·√N instead of ∝ (2p + 1)2 · N . The payoff increases with polynomial degree and spatial dimension. • The Kronecker structure of S2 ⊗ S1 facilitates efficient matrix vector multiplication with ∝ 4(2p + 1) · N oper- ations. • The stiffness matrix need not be explicitly computed. Instead, the matrix vector productK dn can be evaluated in a matrix-free fashion. • The entries of the stiffness matrix Kij = a(Ni/c, Nj) are more complicated to evaluate than the entries to the standard Galerkin stiffness matrix due to the factor 1/c in the test function slot, effectively doubling the cost of evaluation for a second order problem and quadrupling it for a fourth order problem. 5 Numerical results In this section, we verify accuracy and robustness of the proposed methodology. First, we analyze the discrete fre- quency spectrum obtained for a one-dimensional vibrating string that is clamped at either end. Here, we are particu- larly interested in the accuracy of the discrete frequencies obtained with the higher-order accurate dual lumped mass matrix as compared to the consistent and row-sum lumped massmatrices, specifically in the case of prescribed Dirichlet boundary conditions. Next, we study the robustness under mesh distortion in two-dimensional modeling of the wave equation.We also investigate the effect of geometricmapping 123 214 Computational Mechanics (2025) 76:205–226 and outlier removal on the dynamics of a two-dimensional linear membrane. The fourth and last example is a geomet- rically non-linear forced vibration analysis of a cantilever plate, where we again observe that the proposed dual mass lumping scheme is capable of achieving results of the same quality as the consistent mass formulation of IGA. We have combined our techniques with explicit Runge– Kutta time-stepping schemes of the appropriate order O (RKO) to maintain higher-order accuracy in time. The obtained results also underline the importance of employ- ing higher-order accurate discretizations both in space and time. Often, only second order accurate schemes are used in time, which also prevent optimal spatial convergence if the time step size is not selected unrealistically small. To miti- gate the negative effects of outlier frequencies in the high end of the spectrum, the outlier removal procedures pre- sented in [41] are used to improve critical time step values. A MATLAB implementation that reproduces several of the benchmark results in the paper is available for download at: https://github.com/renehiemstra/explicit-dynamics-paper. 5.1 Spectrum analysis The eigenvalue decomposition is an important tool in the study of the dynamic response of a system. It allows one to uncouple the equations such that each degree of freedom is essentially governed by its own uncoupled equation of motion [4, 42]. Each discrete eigenvalue is a measure of the amount of energy associated with one such response. Here, we investigate the behavior of the discrete frequency spectra (the frequency is the square root of an eigenvalue) as a function of the normalized mode number for a one- dimensional second order problem4. The case studied here may be associated with the dynamics of a vibrating string that is fixed at either end. We investigate the spectra obtained for three different mass matrices: (1) consistent mass; (2) row-sum lumped mass; and (3) higher-order accurate dual lumped mass. Appropriate boundary constraints have been prescribed as explained in Sect. 3 and outlier removal proce- dures have been applied as described in [41]. Figure1 displays the frequencies normalized by the corre- sponding analytical values. Hence, results closer to unity are better in terms of accuracy. It may be observed that a larger part of the spectrum is well approximated using the higher- order accurate dual lumped mass compared with row-sum lumped mass. This is the case for all studied polynomial degrees, p = 2 to 5. The gain in relative accuracy compared to the row-sum lumped mass becomes particularly apparent by investigating 4 We refer to our previous work for frequency spectra results associ- ated with fourth order differential operators emanating from beam or plate/shell theories [1]. the relative errors incurred in the frequencies in a semi-log plot, see Fig. 2. The important frequencies obtained at low mode numbers using the higher-order accurate dual lumped mass are very close to those obtained using the consistent mass matrix. This is particularly the case for p = 2 and p = 3. For higher polynomial degrees (p = 4 and p = 5),we observe that the frequencies are much better approximated using the consistent mass matrix. Nevertheless, the gain in accuracy compared to the row-sum lumpedmass ranges from three to six orders in magnitude and the effect increases with polynomial degree. We note that high precision arithmetic has been used in order to maintain and display results up to double precision arithmetic. 5.2 Robustness under mesh distortion We investigate the effect of mesh distortion on numerical solution of the wave equation to asses the robustness under nonlinear mappings. Consider the following solution field displayed in Fig. 3a u(x, y) = sin(4πx) sin(4π y) cos(4π t), (x, y) ∈ [0, 1]2 and t ∈ [0, 1]. (32) The spatial coordinates (x, y) are described in terms of a tensor product cubic Bernstein basis {bi , i = 0, 1, 2, 3} x(ξ, η) = 3∑ i=0 3∑ j=0 Xi j bi (ξ) b j (η) and y(ξ, η) = 3∑ i=0 3∑ j=0 Yi j bi (ξ) b j (η) . (33) with X and Y denoting the control points parameterized by α ∈ [0, 1] according to X = ⎡ ⎢⎢ ⎣ 0 0 0 0 1 3 1 3 + α 1 3 − α 1 3 2 3 2 3 + α 2 3 − α 2 3 1 1 1 1 ⎤ ⎥⎥ ⎦ and Y = ⎡ ⎢⎢ ⎣ 0 1 3 2 3 1 0 1 3 + α 2 3 + α 1 0 1 3 − α 2 3 − α 1 0 1 3 2 3 1 ⎤ ⎥⎥ ⎦ . (34) Figure3b and c illustrate how parameter α distorts the mesh. The distortion in Fig. 3c is obtained for a value of α = 0.5. Hence, the robustness under mesh distortion of the pro- posed dual lumping scheme is compared to consistent mass and lumped mass formulations in Fig. 4. The compari- son is done for polynomial degrees 3, 4, 5 in combination with uniform refinement resulting in meshes with nelem = ( 16, 32, 64, 128 ) elements in both coordinate directions. We compare the approaches by juxtaposing the relative L2 error of the spatial solution at time T = 1 for two different 123 https://github.com/renehiemstra/explicit-dynamics-paper Computational Mechanics (2025) 76:205–226 215 Fig. 1 Normalized frequency spectra associated with consistent mass (black), dual lumped mass (blue), and row-sum lumped mass (red). The results are obtained for polynomial degrees p = 2− 5 with N = 250. The outlier frequencies have been removed using the technique presented in [41]. (Color figure online) meshes: without mesh distortion (α = 0) in the left column and with severe mesh distortion (α = 0.5) in the right col- umn. Regarding time discretization, we use the RK4-scheme when the polynomial degree is 3 or 4, and employ the RK6- scheme when the polynomial degree is 5 (see [18] for an exposition of the algorithms). Considering the standard row- sum lumped mass formulation, we utilize the RK2-scheme. The time increment used is 50% of the computed critical time step size of the respective time-discretization scheme (�tcrit = Cmax/ω h max, with Cmax = 2.0 for RK2, Cmax = 2.785 for RK4 and Cmax = 3.387 for RK6). The numerical results, depicted in Fig. 4, reveal that despite severe mesh distortion, the obtained solution rely- ing on either a consistent or dual mass lumping formulation are practically identical, while the row-sum lumped mass matrix yields considerably worse results. These preliminary findings illustrate that dual mass lumping approach is robust under mesh distortion andmaintains accuracy close to that of consistent mass formulations. Optimal rates can be achieved, even on curved meshes, by using time discretization of the right order. Note that in our experience high-order in space always also demands high-order in time. In addition to the investigation regarding the attainable accuracy under mesh distortion, we also measured the com- putational time for this example. The results are depicted in Fig. 5. We stress that these are only preliminary findings for single-patch configurations and note that this issue deserves a more in-depth study involving multi-patch configurations and complex geometry, something we are actively working on at this stage. The computations have been performed in MATLAB, with the quadrature loops implemented in C. This hybrid approach was chosen because the MATLAB implementation with C extensions is approximately two orders of magnitude faster than a fullMATLAB implementation.However, it is important to emphasize that these results are not representative of large- scale simulations. For small test problems, the consistent mass formulation performs very well since the matrix factor- ization fits in the CPU cache. Unfortunately, this advantage 123 216 Computational Mechanics (2025) 76:205–226 Fig. 2 Normalized frequency errors associated with consistent mass (black), dual lumped mass (blue), and row-sum lumped mass (red). The results are obtained for polynomial degrees p = 2 − 5 with N = 250. The outlier frequencies have been removed using the technique presented in [41]. (Color figure online) Fig. 3 Wave equation on the unit-square: a Exact solution at time T = 1, b non-linear mapping for α = 0, and c non-linear mapping for α = 0.5 123 Computational Mechanics (2025) 76:205–226 217 does not scale well beyond approximately 100,000 degrees of freedom. As a result, comparisons based on such small test problems can be misleading. Given these limitations, the focus of this paper placed on the fundamental aspects of accuracy and robustness under mesh distortion, while recognizing that a more comprehen- sive study involving larger-scale simulations and different computational set-ups is needed. 5.3 Linear analysis of free the vibrations of an annular membrane Consider the free vibration of an annular membrane with inner radius a and outer radius b. The model is fixed at the boundaries. Let J4(r) denote the 4th Bessel function of the first kind and let λk, k = 1, 2, . . . denote its positive roots. The radii of the annulus are chosen conveniently as certain roots of J4(r). In particular, a = λ2 ≈ 11.065 and b = λ4 ≈ 17.616. We define the following manufactured solution, which depends on both the radial coordinate r , the angular coordinate θ , and time t u(r , θ, t) = J4(r) cos(λ2t) cos(4πθ) . (35) It may be verified that the function satisfies the differential equation for free vibrations on the annulus with fixed bound- ary conditions at the inner and outer radii. Figure6 shows the problem setup, the boundary conditions that are satisfied by u, and its initial condition u(r , θ, 0). We perform explicit dynamics simulations with the con- sistent mass matrix, the higher-order accurate dual lumping technique, and the standard row-sum lumping approach. We simulate one full period of the periodic function u(r , θ, t). In other words, the final time is T = 2π/λ2.We use polynomial degrees p = 3, 4, 5 in combination with uniform refine- ment resulting in meshes with nelem = ( 8, 16, 32 ) elements in the radial coordinate and 2nelem = ( 16, 32, 64 ) ele- ments in the angular coordinate. We perform the simulations with and without the outlier removal technique presented in [41]. For the higher-order spatially accuratemethods, namely with consistent mass or dual lumped mass, we discretize the semi-discrete equations in time using the RK4-scheme when the polynomial degree is 3 or 4, and employ the RK6-scheme when the polynomial degree is 5 (see [18] for an exposition of the algorithms). With standard row-sum lumped mass we utilize the RK2-scheme. The time incre- ment used is 50%of the computed critical time step size of the respective time-discretization scheme (�tcrit = Cmax/ω h max, with Cmax = 2.0 for RK2, Cmax = 2.785 for RK4 and Cmax = 3.387 for RK6). The convergence behavior of the threemethods, without andwith outlier removal, is presented in Fig. 7 as a function of the square root of the total number of degrees of freedom. The higher-order accurate dual lumping technique preserves the accuracy that is obtained using the consistent mass. In addition, an outlier-free basis preserves optimal accuracy of the standard approach, however, using fewer degrees of freedom and larger time step sizes. The improvement in the critical time step size can be observed in Fig. 8. Figure8a displays the improvement due to the choice of mass matrix (consistent, row-sum lumped, or higher-order accurate dual lumped mass) and Fig. 8b demonstrates the additional benefit of outlier removal. The compound effect is obtained by multiplying these numbers, which easily leads to a doubling of the critical time step size without loss of accuracy. 5.4 Nonlinear analysis of forced vibration of a square cantilever plate The last benchmark example considers a cantilever plate as shown in Fig. 9, modeled using n × n (n = 4, 8, 16, 32) isogeometric Kirchhoff-Love shell elements of polynomial degrees p = 2, 3, 4 [43]. The plate is subjected to a time dependent lineload and will undergo large displacements (geometrically non-linear analysis) with a linear elastic con- stitutive material law (Saint Venant-Kirchhoff constitutive model). The material properties, displacement boundary conditions and external loading are depicted in Fig. 9 as well. Considering the time integration algorithm, we employ Runge–Kutta methods of orders 2 (RK2) and 4 (RK4) for polynomial degrees p = 2 and p = 3, 4, respectively5. The time increment for the RK-schemes is chosen as mentioned in the previous section and is based on the finest mesh to ensure that all simulations of the same polynomial degree are run with the same value of �t . Note that this example is an adapted version of the example presented in [44], where the excitation time tc and the loading have been adjusted to excite higher frequencies/modes in the dynamic response. As in the previous examples, the numerical results obtained using the consistent mass, the row-sum lumped mass, and the higher-order accurate dual lumped mass are compared. The deformed configurations for the proposed set- up at t = 0.39 s are shown in Fig. 10 for a mesh of 32 × 32 cubic isogeometric shell elements (p= 3). While the dis- placement fields obtained with the conventional consistent mass (see Fig. 10a) and our novel dual mass lumped Petrov– Galerkin formulation (see Fig. 10c) agree exceptionally well, we observe distinct differences compared to the row-sum lumped formulation (see Fig. 10b). This is especially easy to spot at the plate edge, where the loading was applied, and in the middle of the plate, where a pronounced bump is seen for the consistent and dual-lumped mass. 5 The order of accuracy of the time integration scheme is chosen accord- ing to the degree of the spatial discretization. This is done to balance the spatial and temporal discretization errors. 123 218 Computational Mechanics (2025) 76:205–226 Fig. 4 Relative L2 error of the spatial solution at time T = 1.0 obtained with polynomial degree p = 3, 4, 5 and geometry parameter α = 0.0 (left column) and α = 0.5 (right column) 123 Computational Mechanics (2025) 76:205–226 219 Fig. 5 Time versus relative L2 error at time T = 1.0 obtained with polynomial degree p = 3, 4, 5 and geometry parameter α = 0.0 (left column) and α = 0.5 (right column) 123 220 Computational Mechanics (2025) 76:205–226 Fig. 6 Problem description for explicit dynamics on an annulus Fig. 7 Relative L2 error in the vertical displacement field u as a function of the square root of number of degrees of freedom N . a Dual mass lumping maintains the accuracy of the consistent mass. b Outlier removal does not negatively affect the accuracy of the methods The time response of the vertical (out-of-plane) displace- ment at observation point A is shown in Figs. 11 and 12 for polynomial degrees 2 and 3, respectively. The results again confirm the the proposed higher-order accurate dual lumping approach is capable of achieving the same level of accuracy as the consistent mass, while much larger discrep- ancies are noted with respect to the solution of the row-sum lumped mass. Even for a mesh of 32×32 isogeometric plate elements of degree 3, no converged solution is reached if row-summing is applied and differences w.r.t. the other two solutions are still very obvious. In Fig. 13, the displacement history is depicted for the mesh consisting of 32×32 isogeo- metric shell elements with different polynomial degrees. We can conclude from the obtained results, that the solution is 123 Computational Mechanics (2025) 76:205–226 221 Fig. 8 Increase of the critical time step size of the dual lumped mass compared with the consistent mass (a) and the increase in time step size due to outlier removal (b). The compound effect of the chosen mass and outlier removal is obtained by multiplying these numbers Fig. 9 Problem description of large deformation analysis of a square cantilever plate under a lineload Fig. 10 Deformed configurations of the cantilever plate visualized at t = 0.39 s (point of maximum displacement). The results are depicted for a mesh of 32 × 32 isogeometric shell elements with cubic basis functions (p = 3). The displacement field is plotted with a scaling factor of s = 5 123 222 Computational Mechanics (2025) 76:205–226 Fig. 11 Time history of the displacement at point A (uz) on a free edge of the plate under vanishing lineload, computed with p = 2 for the three explicit isogeometric schemes sufficiently accurate for p ≥ 3. In the following, the results that have been discussed up to this point are compared to established algorithms provided by commercial tools, which again highlights the correctness of our implementation. In addition to the IGA solutions based on an in-house code developed in Julia, the commercial software ABAQUS has been used to simulate the same example and get the dynamic response. Here, the cantilever plate is modeled by means of three-dimensional solid elements of Serendipity type referred to as C3D20 (p= 2). The mesh consists of 200× 200× 2 elements, featuring 443,205 nodes with each three displacement degrees of freedom. In total, we have a model with 1,329,615 degrees of freedom. The bound- ary conditions are implemented as a pinned support, i.e., all displacement degrees of freedom of the nodes located on the corresponding face of the plate are fixed. The lineload is applied as a distributed surface traction, where the load is applied on the entire face of the plate and therefore, the amplitude is set to pSFmax = dpmax. Due to the limitation that elementswith quadratic (p= 2) shape functions are not avail- able in ABAQUS/Explicit, we employed an implicit time integration scheme in ABAQUS/Standard. The default option is the HHT-α-method [45]. The parameter α is set to −0.3, which introduces numerical damping/dissipation in the simulation. The time step size has been selected as �tHHT = 5 × 10−4 s as the method itself is unconditionally stable. The chosen value makes sure that the application of the force is not too fast/abrupt, as 20 time steps cover the loading and unloading of the cantilever plate. The results are depicted in Fig. 14 and are in excellent agreement with the solution obtained by our higher-order accurate dual mass lumping version of IGA. This again verifies the correctness 123 Computational Mechanics (2025) 76:205–226 223 Fig. 12 Time history of the displacement at point A (uz) on a free edge of the plate under vanishing lineload, computed with p = 3 for the three explicit isogeometric schemes of our implementation and highlights the excellent conver- gence properties of the dualmass lumping scheme. The small differences in the results are attributed to the fact that in our IGA implementation we use higher-order elements featuring higher interelement continuity in conjunction with higher- order explicit Runge–Kutta time integration schemes. On the other hand, numerical damping is included in the implicit HHT-α method, which is not present in the RK-solvers to were implemented in the in-house Julia code. It must be noted that employing isogeometric Kirchhoff-Love shell ele- ments and three-dimensional hexahedral elements is another source of the observed differences. This obviously includes the different manner in which the Dirichlet boundary condi- tions and the loading are applied.While in the IGAmodel the lineload is directly applied to an edge of the shell element, we have to resort to a surface traction acting on a face of a hexahedral element in the fully three-dimensional case. Overall, the results obtained by the novel Petrov–Galerkin dual mass lumping approach are very promising and consti- tute the starting point for further research into mass lumping and IGA. In this context, complex multi-patch domains and trimmed shell surfaces are of special interest. 6 Conclusion This paper presents a comprehensive mathematical frame- work for explicit structural dynamics, utilizing approximate dual functionals and row-summass lumping.Wedemonstrate that the Petrov–Galerkin method with row-sum mass lump- ing proposed in [1] can also be interpreted as a Galerkin 123 224 Computational Mechanics (2025) 76:205–226 Fig. 13 Time history of the displacement at point A (uz) on a free edge of the plate under vanishing lineload, computed with the novel dual mass lumped Petrov–Galerkin IGA formulation for p ∈ [2, 3, 4] and 32 × 32 elements Fig. 14 Time history of the displacement at point A (uz) on a free edge of the plate under vanishing lineload, computed with ABAQUS using three-dimensional solid elements (C3D20) and the novel dual mass lumped Petrov–Galerkin IGA formulation method with a customized mass matrix. Based on this observation, we prove mathematically that the customized mass technique does not compromise the asymptotic accu- racy of the method. An essential improvement with respect to our prior work in [1] is the incorporation of Dirichlet boundary conditions, while maintaining higher-order accu- racy. The mathematical results are substantiated by spectral analysis and validated through a two-dimensional linear explicit dynamics benchmark and a three-dimensional (geo- metrically non-linear) forced vibration analysis. Our method exhibits accuracy and robustness comparable with aGalerkin method employing a consistent mass, while retaining the explicit nature of a lumped mass. The excellent coarse mesh accuracy of the approach paves the way for efficient explicit dynamics on coarse meshes with commensurable critical time step sizes. In future work, we intend to extend our to multi-patch configurations of complex geometric models found for example in the automotive industry. 123 Computational Mechanics (2025) 76:205–226 225 Acknowledgements The authors gratefully acknowledge the finan- cial support provided by the German Research Foundation (Deutsche Forschungsgemeinschaft –DFG) through the research grants EI 1188/3- 1 (project number: 497531141) and SCH 1249/5-1 (Project Number: 490700327), and theEmmyNoether grant SCH1249/2-1 (ProjectNum- ber: 326309100). Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adap- tation, 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 indi- cate 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, youwill need to obtain permission directly from the copy- right holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. References 1. Nguyen T-H, Hiemstra RR, Eisenträger S, Schillinger D (2023) Towards higher-order accurate mass lumping in explicit isogeo- metric analysis for structural dynamics. Comput Methods Appl Mech Eng 116233 2. Belytschko T, Lin JI, Tsay C-S (1984) Explicit algorithms for the nonlinear dynamics of shells. Comput Methods Appl Mech Eng 42(2):225–251 3. Benson DJ (1992) Computational methods in Lagrangian and Eulerian hydrocodes. Comput Methods Appl Mech Eng 99(2– 3):235–394 4. Hughes TJR (2012) The finite element method: linear static and dynamic finite element analysis. Courier Corporation, Baltimore 5. Cottrell JA, Hughes TJR, Bazilevs Y (2009) Isogeometric analysis: toward integration of CAD and FEA 6. Hughes TJR, Cottrell JA, Bazilevs Y (2005) Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refine- ment. Comput Methods Appl Mech Eng 194:4135–4195 7. Haberleitner M, Jüttler B, Scott M, Thomas D (2017) Isogeomet- ric analysis: representation of geometry. In: Stein E, de Borst R, Hughes TJR (eds) Encyclopedia of computational mechanics, 2nd edn. Wiley, Hoboken 8. Hughes TJR, Sangalli G (2017)Mathematics of isogeometric anal- ysis: a conspectus. In: Stein E, de Borst R, Hughes TJR (eds) Encyclopedia Comput Mech, 2nd edn. Wiley, Hoboken 9. Marussig B, Hughes TJR (2018) A review of trimming in isogeo- metric analysis: challenges, data exchange and simulation aspects. Arch Comput Methods Eng 25(4):1059–1127 10. Schillinger D (2018) Isogeometric finite element analysis. In: AltenbachH,ÖchsnerA (eds)Encyclopedia of continuummechan- ics. Springer, New York 11. Schillinger D, Evans JA, Reali A, Scott MA, Hughes TJR (2013) Isogeometric collocation: cost comparison with Galerkin meth- ods and extension to adaptive hierarchical NURBS discretizations. Comput Methods Appl Mech Eng 267:170–232 12. Willberg C, Duczek S, Vivar Perez J, Schmicker D, Gabbert U (2012)Comparison of different higher order finite element schemes for the simulation of Lamb waves. Comput Methods Appl Mech Eng 241–244:246–261 13. Cottrell JA, Reali A, Bazilevs Y, Hughes TJR (2006) Isogeometric analysis of structural vibrations. Comput Methods Appl Mech Eng 195:5257–5296 14. Hartmann S, Benson DJ (2015) Mass scaling and stable time step estimates for isogeometric analysis. Int J Numer Methods Eng 102(3–4):671–687 15. BensonDJ,BazilevsY,HsuM-C,HughesTJR (2010) Isogeometric shell analysis: the Reissner–Mindlin shell. Comput Methods Appl Mech Eng 199:276–289 16. Chen L, Nguyen-Thanh N, Nguyen-Xuan H, Rabczuk T, Bor- das SPA, Limbert G (2014) Explicit finite deformation analysis of isogeometric membranes. Comput Methods Appl Mech Eng 277:104–130 17. Leidinger LF, Breitenberger M, Bauer AM, Hartmann S, Wüchner R, Bletzinger K-U, Duddeck F, Song L (2019) Explicit dynamic isogeometric B-Rep analysis of penalty-coupled trimmed NURBS shells. Comput Methods Appl Mech Eng 351:891–927 18. Evans JA, Hiemstra RR, Hughes TJR, Reali A (2018) Explicit higher-order accurate isogeometric collocation methods for struc- tural dynamics. Comput Methods Appl Mech Eng 338:208–240 19. Voet Y, Sande E, Buffa A (2023) A mathematical theory for mass lumping and its generalization with applications to isogeometric analysis. Comput Methods Appl Mech Eng 410:116033 20. Voet Y, Sande E, Buffa ARobust mass lumping and outlier removal strategies in isogeometric analysis 21. CanutoC,HussainiM,QuarteroniA, ZangT (2007) Spectralmeth- ods: evolution to complex geometries and applications to fluid dynamics. Springer, New York 22. Duczek S, Gravenkamp H (2019) Mass lumping techniques in the spectral element method: on the equivalence of the row-sum, nodal quadrature, and diagonal scaling methods. Comput Methods Appl Mech Eng 353:516–569 23. Schillinger D, Evans JA, Frischmann F, Hiemstra RR, Hsu MC, Hughes TJR (2014) A collocated C0 finite element method: Reduced quadrature perspective, cost comparison with standard finite elements, and explicit structural dynamics. Int JNumerMeth- ods Eng 102(3–4) 24. Nguyen LH, Schillinger D (2017) A collocated isogeometric finite element method based on gauss-Lobatto Lagrange extraction of splines. Comput Methods Appl Mech Eng 316:720–740 25. Schillinger D, Ruthala PK, Nguyen LH (2016) Lagrange extrac- tion and projection for nurbs basis functions: a direct link between isogeometric and standard nodal finite element formulations. Int J Numer Methods Eng 108(6):515–534 26. Anitescu C, Nguyen C, Rabczuk T, Zhuang X (2019) Isogeometric analysis for explicit elastodynamics using a dual-basis diagonal mass formulation. Comput Methods Appl Mech Eng 346:574–591 27. Schumaker L (2007) Spline functions: basic theory. Cambridge University Press, Cambridge 28. Seitz A, Farah P, Kremheller J, Wohlmuth BI, Wall WA, Popp A (2016) Isogeometric dual mortar methods for computational con- tact mechanics. Comput Methods Appl Mech Eng 301:259–280 29. Dornisch W, Stöckler J, Müller R (2017) Dual and approximate dual basis functions for B-splines and NURBS - Comparison and application for an efficient coupling of patches with the isogeomet- ric mortar method. ComputMethods ApplMech Eng 316:449–496 30. Zou Z, Scott MA, Borden MJ, Thomas DC, Dornisch W, Brivadis E (2018) Isogeometric Bézier dual mortaring: refineable higher- order spline dual bases and weakly continuous geometry. Comput Methods Appl Mech Eng 333:497–534 31. Dornisch W, Stöckler J (2021) An isogeometric mortar method for the coupling of multiple nurbs domains with optimal convergence rates. Numerische Mathematik 149(4):871–931 32. Zou Z, Scott MA, Miao D, Bischoff M, Oesterle B, Dor- nisch W (2020) An isogeometric Reissner–Mindlin shell element based on Bézier dual basis functions: overcoming locking and improved coarse mesh accuracy. ComputMethods ApplMech Eng 370:113283 123 http://creativecommons.org/licenses/by/4.0/ http://creativecommons.org/licenses/by/4.0/ 226 Computational Mechanics (2025) 76:205–226 33. Chui CK, He W, Stöckler J (2004) Nonstationary tight wavelet frames, I: Bounded intervals. Appl Comput Harmon Anal 17(2):141–197 34. Held S, Eisenträger S, DornischW (2024) An efficient mass lump- ing scheme for isogeometric analysis based on approximate dual basis functions. Technische Mechanik - Eur J Eng Mech 44:14–46 35. Li X, Wang D (2022) On the significance of basis interpolation for accurate lumpedmass isogeometric formulation. ComputMethods Appl Mech Eng 400:115533 36. Mika ML, Hughes TJR, Schillinger D, Wriggers P, Hiemstra RR (2021)Amatrix-free isogeometricGalerkinmethod forKarhunen– Loève approximation of randomfields using tensor product splines, tensor contraction and interpolation based quadrature. Comput Methods Appl Mech Eng 379:113730 37. Zienkiewicz OC, Taylor RL (2013) The finite element method: its basis and fundamentals. Butterworth-Heinemann, Oxford 38. Duczek S, Gravenkamp H (2019) Critical assessment of different mass lumping schemes for higher order serendipity finite elements. Comput Methods Appl Mech Eng 350:836–897 39. de Boor C (1978) A practical guide to splines. Springer, New York 40. Golub GH, Van Loan CF (2013) Matrix computations. JHU Press, Baltimore 41. Hiemstra RR, Hughes TJ, Reali A, Schillinger D (2021) Removal of spurious outlier frequencies and modes from isogeometric dis- cretizations of second- and fourth-order problems in one, two, and three dimensions. Comput Methods Appl Mech Eng 387:114115 42. Chopra AK (2019) Dynamics of structures in SI units, 5th edn. International Series in Civil Engineering and EngineeringMechan- ics, Pearson Education Limited 43. Kiendl J, BletzingerK-U, Linhard J,Wüchner R (2009) Isogeomet- ric shell analysis with Kirchhoff–Love elements. Comput Methods Appl Mech Eng 198(49):3902–3914 44. Liu T, Zhao C, Li Q, Zhang L (2012) An efficient backward Euler time-integration method for nonlinear dynamic analysis of struc- tures. Comput Struct 106–107:20–28 45. Hilber HM, Hughes TJR, Taylor RL (1977) Improved numerical dissipation for time integration algorithms in structural dynamics 5:283–292 Publisher’s Note Springer Nature remains neutral with regard to juris- dictional claims in published maps and institutional affiliations. 123 Higher-order accurate mass lumping for explicit isogeometric methods based on approximate dual basis functions Abstract 1 Introduction 2 Preliminaries 2.1 Model problem 2.2 Semi-discrete equations of motion 2.3 Explicit time stepping methods 2.4 Efficient update through mass lumping 3 Approximate L2 dual bases 3.1 Univariate B-splines 3.2 L2 dual basis 3.3 L2 approximate dual basis 3.4 Treatment of Dirichlet boundary conditions 3.5 Perspective 1: Galerkin method with a higher-order approximate mass matrix 3.6 Perspective 2: Petrov–Galerkin method with a lumped mass matrix 4 Application in higher-order accurate explicit dynamics 4.1 Extension to arbitrary isoparametric mappings 4.2 Maintaining higher-order accuracy in explicit dynamics 4.3 Imposition of Dirichlet boundary conditions 4.4 Implementational aspects 5 Numerical results 5.1 Spectrum analysis 5.2 Robustness under mesh distortion 5.3 Linear analysis of free the vibrations of an annular membrane 5.4 Nonlinear analysis of forced vibration of a square cantilever plate 6 Conclusion Acknowledgements References