Inverse mass matrix for isogeometric explicit transient analysis via the method of localized Lagrange multipliers

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ Esta es la versión aceptada del artículo publicado en: This is a accepted manuscript of a paper published in: International Journal for Numerical Methods in Engineering, Volume 117, Issue 9, March 2019 DOI: https://doi.org/10.1002/nme.5986 Copyright: © 2018 John Wiley & Sons, Ltd. El acceso a la versión publicada del artículo puede requerir la suscripción de la revista. Access to the published version may require subscription. “This is the peer reviewed version of the following article: [González JA, Kopačka J, Kolman R, Cho SS, Park KC. Inverse mass matrix for isogeometric explicit transient analysis via the method of localized Lagrange multipliers. Int J Numer Methods Eng. 2019; 117: 939– 966. https://doi.org/10.1002/nme.5986], which has been published in final form at [https://doi.org/10.1002/nme.5986]. This article may be used for non-commercial purposes in accordance with Wiley Terms and Conditions for Use of Self-Archived Versions. This article may not be enhanced, enriched or otherwise transformed into a derivative work, without express permission from Wiley or by statutory rights under applicable legislation. Copyright notices must not be removed, obscured or modified. The article must be linked to Wiley’s version of record on Wiley Online Library and any embedding, framing or otherwise making available the article or pages thereof by third parties from platforms, services and websites other than Wiley Online Library must be prohibited." Received 11 October 2017; Revised 11 October 2017; Accepted 11 October 2017 DOI: xxx/xxxx RESEARCH ARTICLE Inverse Mass Matrix for Isogeometric Explicit Transient Analysis via the Method of Localized Lagrange Multipliers José A. González1| J. Kopačka2| R. Kolman2| S.S. Cho3| K.C. Park4 1Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n, Seville E-41092, Spain 2Institute of Thermomechanics, The Czech Academy of Sciences, Dolejškova 5, 182 00 Prague, Czech Republic 3Reactor Mechanical Engineering Division, Korea Atomic Energy Research Institute, 999-111 Daedeok-Daero, Yuseong-gu, Daejeon 305-353, Republic of Korea 4Department of Aerospace Engineering Sciences, University of Colorado at Boulder, Boulder, CO 80309-429, USA Correspondence José A. González, Camino de los Descubrimientos s/n, Seville 41092, Spain. Email: [email protected] Present Address Camino de los Descubrimientos s/n, Seville 41092, Spain Abstract A variational framework is employed to generate inverse mass matrices for isogeometric analysis. As different dual bases impact not only accuracy but also computational overhead, several dual bases are extensively investigated. Specifically, locally discontinuous bi-orthogonal basis functions are evaluated in detail for B-splines of high continuity and Bézier elements with a standard 𝐶0continuous finite element structure. The boundary conditions are enforced by the method of localized Lagrangian multipliers after generating the inverse mass matrix for completely free body. Thus, unlike inverse mass matrix methods without employing the method of Lagrange multipliers, no modifications in the reciprocal basis functions are needed to account for the boundary conditions. Hence, the present method does not require internal modifications of existing Isogeometric analysis software structures. Numerical examples show that globally continuous dual basis functions yield better accuracy than locally discontinuous bi-orthogonal functions, but with much higher computational overhead. Locally discontinuous dual basis functions are found to be an economical alternative to lumped mass matrices when combined with mass parametrization. The resulting inverse mass matrices are tested in several vibration problems and applied to explicit transient analysis of structures. KEYWORDS: Inverse mass matrix; isogeometric analysis; Bézier extraction; localized Lagrange multipliers; partitioned analysis; free vibration, explicit transient analysis 1 INTRODUCTION Efficient explicit time integration of structural dynamics equations calls for the availability of diagonal mass matrices. This is because the computational cost associated with the use of consistent/nondiagonal mass matrices in an explicit integration is equivalent with that of implicit time integration procedures, whose time step sizes can be chosen much larger than for explicit integration procedures. Specifically, explicit time integration methods in transient dynamics require a repeated application of the inverse mass matrix to obtain accelerations at every time step. For example, computing the inverse of a consistent mass matrix requires (𝑛3)flops. This has led to the adoption of the diagonalized lumped mass matrix (LMM) that is typically obtained by a row-sum of the consistent mass matrix (CMM). However, even though the accuracy of LMM is acceptable for capturing low-frequency response components, it is unable to track high-frequency components1,2,3,4,5. Another shortcoming of LMM for higher-order elements is that most of the mass is lumped at the center node, causing very strange behavior when high-order 2GONZÁLEZ ET AL modes are excited. These deficiencies motivated researchers to find ways of generating the inverse of consistent mass matrix in a consistent and efficient way6,7,8. From the viewpoint of principles in dynamics, there are two ways for computing the mass matrices, which need to be inverted for efficient explicit time integration. One approach is to directly discretize system kinetic energy. The second is to exploit Hamilton’s principle (or Hamilton’s equations) that require discretization of both momentum (𝒑) and displacement (𝒖). In this case, the inverse mass matrix (𝐌−1) is obtained as a second-order tensor product, where only a diagonalized projection operator (𝐀which will be detailed later in the paper) needs to inverted, viz.: (∫ Ω 1 𝜌𝛿𝒑⋅𝒑𝑑Ω ≈ 𝛿𝐩𝐂𝐩,∫ Ω 𝛿𝒑⋅𝒖𝑑Ω ≈ 𝛿𝐩 𝐀𝖳𝐮)⇒𝐌−𝟣=𝐀−𝖳𝐂𝐀−𝟣(1) where 𝐂is labeled as a momentum matrix, 𝐩denotes discretized momentum vector. Hence, inversion of mass matrix is accomplished via the generation of the momentum matrix (𝐂) followed by a push-forward projection operation. Thus, an efficient generation of consistent inverse mass matrices depends on how efficiently the projection operator (𝐀) can be generated and inverted. This in turn requires that the interpolation functions for the momentum must be bi-orthogonal with respect to the displacement interpolation functions so that the resulting projection matrix (𝐀) becomes diagonal. While derivations of bi-orthogonal basis functions are straightforward for most interpolation functions for the interior elements, for the elements bordering the boundaries, special bi-orthogonal basis functions need to be constructed. This construction has to performed for each of the corresponding 𝐮-interpolation functions. It should be noted that, by the very nature of finite element method, elements bordering the boundaries do not need to be specialized. The preceding complications in generating bi-orthogonal basis functions for the momentum vector in treating the boundary conditions were alleviated in the work of González et al.9for FEM by introducing the method of localized Lagrange multipliers10,11,12,13. In that work, the system is considered completely free floating such that no special treatment is required for generating the bi-orthogonal interpolation of the momentum vector. Enforcement of boundary conditions is simply accomplished via constraining the boundary conditions. We believe that the same ideas used for the obtention of FEM inverse mass matrices, can be extended to the more general framework of isogeometric analysis (IGA). Particularly for IGA, different approaches have been proposed in the literature to address this difficulty present in the mass matrix inversion. First contributions in explicit dynamics were developed for B-spline collocation methods, treating the equations of motion (EOM) in strong form and obtaining one-point quadrature mass matrices14,15. Solutions to the mass inversion problem were then suggested at the time-integration solution algorithm level, introducing a predictor multi-corrector scheme that predicts accelerations with LMM on the left-hand-side of the linear system and corrects the residual using the CMM. This method is iterative in essence and energy stability is not guaranteed. Other methods benefit from the regularity of the reference patches used in IGA, taking advantage of the tensor product separable form of the mass matrix associated to these structured meshes. For example, at the linear solver level with a CMM, Gao et al.16 proposed using a preconditioned conjugate gradient method with special preconditioners that achieves a complexity of (𝑛)flops. However, this procedure is not completely general and non-separable mass matrices require a carefully designed preconditioner that adds significant complexity and computational cost to the linear solver. New preconditioning and multigrid techniques with linear cost17,18,19 could also be considered for that task. In the context of multi-patch discontinuous Galerkin IGA for acoustic wave propagation and linear advection, a recent work of Chan et al.20 aims to approximate the mass matrix inverse by a computable weight-adjusted mass matrix that is pre-post multiplied by an easily invertible affine mass matrix obtained from the tensor product of 1D inverses. Although not derived from energy principles, the proposed algebraic approximation is demonstrated to be high order accurate only for sufficiently regular discretizations. A different approach, the methodology used in this paper belongs to this last group, tries to derive a reciprocal mass matrix (RMM) directly from parametrized variational principles and dual-space shape functions21. The idea here is to derive variationally consistent reciprocal mass matrices from a parametrized functional and adjust the value of these parameters to construct customizable RMMs with desired properties, like high accuracy on certain frequency ranges or extending the critical time-step needed for time integration stability. The objective of the present paper is to extend the procedure for generating FEM inverse mass matrices presented in González et al.9for explicit transient analysis by utilizing existing IGA capabilities22,23,24,25. To this end, the paper is organized as follows: Section 2 briefly reviews some basic concepts on B-splines, NURBS and Bézier elements interpolation. Section 3 derives the displacement-momentum mixed equations of elastodynamics from Hamilton’s principle and Section 4 describes the direct GONZÁLEZ ET AL 3 construction of inverse mass matrices through bi-orthogonal basis functions. Sections 5 and 6 categorize and analyze the construction of dual basis functions. Section 7 discusses the incorporation of boundary conditions into the inverse mass matrix via the method of localized Lagrange multipliers. Finally, in Section 8 we present several numerical examples, where we test the accuracy of the proposed mass matrix inverse in vibration and transient analyses, and Section 9 closes with the conclusions. 2A BRIEF REVIEW OF BASICS OF ISOGEOMETRIC ANALYSIS In this Section we review the basic concepts of the isogeometric analysis techniques that will be used to construct the proposed reciprocal mass matrices. We start with the classical B-spline and NURBS approximation and next outline the basics of the Bézier extraction technique. This last technique is utilized to produce a finite element data structure of the geometry with the same approximation than the B-spline basis functions. 2.1 B-splines and NURBS-based functions In NURBS based Galerkin IGA, the approximation of a displacement field 𝒖(𝐱)uses the same mathematical representation than the geometry based on spline-like approximation techniques26. For three-dimensional problems in particular, the geometry of a domain Ω ∈ ℝ3is described by the tensor product formula: 𝐱(𝜉, 𝜂, 𝜁)= 𝑛 ∑ 𝑖=1 𝑚 ∑ 𝑗=1 𝑙 ∑ 𝑘=1 𝑅𝑝,𝑞,𝑟 𝑖,𝑗,𝑘 (𝜉, 𝜂, 𝜁)𝐏𝑖,𝑗,𝑘,(2) where 𝐱∈ Ω marks the position of point of interest in the physical space, (𝜉, 𝜂, 𝜁)are its associated coordinates in the parametric space, 𝑅𝑝,𝑞,𝑟 𝑖,𝑗,𝑘 are NURBS basis functions of orders (𝑝, 𝑞, 𝑟)on each dimension, and 𝐏𝑖,𝑗,𝑘 with 𝑖= 1,…, 𝑛, 𝑗 = 1,…, 𝑚, 𝑘 = 1,…, 𝑙 is a matrix with the three coordinates of the control points defining the geometry. An associated displacement field 𝒖(𝐱)is then approximated as follows: 𝒖(𝐱(𝜉, 𝜂, 𝜁)) = 𝑛 ∑ 𝑖=1 𝑚 ∑ 𝑗=1 𝑙 ∑ 𝑘=1 𝑅𝑝,𝑞,𝑟 𝑖,𝑗,𝑘 (𝜉, 𝜂, 𝜁)𝐮P 𝑖,𝑗,𝑘,(3) where 𝐮P 𝑖,𝑗,𝑘 are now the displacements corresponding to the control point indicated by indices 𝑖, 𝑗, 𝑘. For NURBS solids, the shape functions are described by the tensor product as 𝑅𝑝,𝑞,𝑟 𝑖,𝑗,𝑘 (𝜉, 𝜂, 𝜁)= 𝑁𝑝 𝑖(𝜉)𝑁𝑞 𝑗(𝜂)𝑁𝑟 𝑘(𝜁)𝑤𝑖,𝑗,𝑘 ∑𝑛  𝑖=1 ∑𝑚  𝑗=1 ∑𝑙  𝑘=1 𝑁𝑝  𝑖(𝜉)𝑁𝑞  𝑗(𝜂)𝑁𝑟  𝑘(𝜁)𝑤 𝑖,  𝑗, 𝑘 ,(4) where 𝑤𝑖,𝑗,𝑘 is referred to as the 𝑖, 𝑗, 𝑘-th weight corresponding to the 𝑖, 𝑗, 𝑘-th control point, and 𝑁𝑝 𝑖(𝜉), 𝑁𝑞 𝑗(𝜂), 𝑁𝑟 𝑘(𝜁)are univariate B-spline basis functions of orders 𝑝,𝑞and 𝑟corresponding to the knot vectors Ξ,,, respectively27. Let us define a knot vector in one-dimensional case as a non-decreasing set of coordinates in the parametric space, written as Ξ={𝜉1, 𝜉2,…, 𝜉𝑚𝑘}, where 𝜉𝑖∈ℝis the 𝑖-th knot, 𝑖is the knot index, 𝑖= 1,2,…, 𝑚𝑘, number of knots 𝑚𝑘=𝑛+𝑝+ 1, order 𝑝is the polynomial degree and 𝑛is the number of basis functions (control points) used to construct the B-spline curve. Associated B-spline basis functions 𝑁𝑝 𝑖(𝜉), see Figure 1 , are 𝐶𝑝−𝑘continuous piecewise polynomials defined by the Cox-de Boor recursion formula28. Observe that 𝑁𝑝 𝑖are 𝐶𝑝−𝑘continuous piecewise polynomials, where 𝑘is the knot multiplicity, and that continuity can be controlled by the multiplicity of a knot in the knot vector Ξ. Finally, the NURBS based interpolation in ℝ3of a displacement field 𝒖(𝐱)can then be expressed in matrix form as: 𝒖(𝝃) = 𝑛 ∑ 𝑖=1 𝐍𝑝 𝑖(𝝃)𝐮P 𝑖=𝐍𝑢(𝝃)𝐮,(5) where the matrix of interpolation functions 𝐍𝑢(𝝃)is constructed by using B-spline basis functions of a given degree 𝑝and vector 𝐮groups the unknown discrete values at the control points. 4GONZÁLEZ ET AL 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 FIGURE 1 Standard one-dimensional quadratic p=2 (left) and cubic p=3 (right) B-spline basis functions 𝑁𝑝 𝑢(𝜉)in the normalized parameter space for an open knot vector with internal values Ξ = {0,0.2,0.4,0.6,0.8,1}. 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 FIGURE 2 Result of the application of the Bézier extraction process to the original quadratic p=2 (left) and cubic p=3 (right) B-spline basis functions represented in Figure 1 . After knot insertion, Ξ = {0,0.2,0.2,0.4,0.4,0.6,0.6,0.8,0.8,1} for 𝑝= 2, the original B-spline basis functions can be expressed as a linear combination of the new Bézier interpolation functions with compact support. 2.2 Bézier extraction - transformation of NURBS based mesh into Bézier based mesh As mentioned in the previous section, the support of a B-spline basis function 𝑁𝑝 𝑖is compact and contained over 𝑝+1 knot spans in the interval [𝜉𝑖, 𝜉𝑖+𝑝+1]. As a result, it is not possible to use the concept of a single parent finite element with locally defined shape functions that is mapped onto individual knot intervals, complicating the implementation. One solution to circumvent this difficulty is the Bézier’s extraction process29. In the same way that NURBS are a generalization of B-Splines, B-Splines rerpesent the generalization of Bézier splines (based on Bernstein polynomials) and as such, the Bézier representation of a domain can be written in the form: 𝐱(𝜉, 𝜂, 𝜁)= 𝑛 ∑ 𝑖=1 𝑚 ∑ 𝑗=1 𝑙 ∑ 𝑘=1 𝐵𝑝,𝑞,𝑟 𝑖,𝑗,𝑘 (𝜉, 𝜂, 𝜁)𝐏𝐵 𝑖,𝑗,𝑘 (6) where 𝐏𝐵 𝑖,𝑗,𝑘 are new coordinates of control points for a Bézier based representation and 𝐵𝑝,𝑞,𝑟 𝑖,𝑗,𝑘 are the shape functions based on the Bernstein polynomials of order 𝑝, 𝑞, 𝑟. GONZÁLEZ ET AL 5 0 1/3 2/3 1 0 1/3 2/3 1 0 1/3 2/3 1 a) b) c) FIGURE 3 Bézier extraction operation performed in a 2D cubic (𝑝= 3) B-spline curve. The shape of the initial 𝐶2continuous curve (a) is defined by six control points and six basis functions, and its associated knot vector is Ξ = {0,0,0,0,1∕3,2∕3,1,1,1,1}. After a first repetition of knot values, one gets the curve with 𝐶1continuity (b) with eight control points/basis functions and knot vector Ξ={0,0,0,0,1∕3,1∕3,2∕3,2∕3,1,1,1,1}. Finally, 𝑝repetitions of the knot values results in the 𝐶0continuous curve (c) with ten control points/basis functions and knot vector Ξ = {0,0,0,0,1∕3,1∕3,1∕3,2∕3,2∕3,2∕3,1,1,1,1}. Now the basis functions are the same for all the knot intervals, therefore the curve can be constructed through a systematic assembling of Bézier elements with Bernstein polynomials as shape functions. By Bézier extraction, it is meant the construction of a mesh of isogeometric Bézier elements from a NURBS representation. This process provides a local element support structure for isogeometric analysis, that can be easily incorporated into existing finite element codes, see Figure 3 . Following this technique, the B-spline functions 𝐍(𝝃)can be unambiguously written as linear combination of Bernstein polynomials 𝐁(𝝃), expressed in matrix form: 𝐍(𝝃) = 𝐂𝐵𝐁(𝝃),(7) where matrix 𝐂𝐵is the extraction operator, that is uniquely associated to a particular knot vector29. Technically, the new mesh is constructed step-by-step repeating all the knots contained in the original knot vector up to multiplicity 𝑝+ 1, reducing at the same time continuity to 𝐶0, by following the process described in Figure 3 . Then, for each B-spline segment, an expression like (7) can be expressed in local form. Furthermore, new position and weights of the Bézier mesh, 𝐏𝐵and 𝐰𝐵, can be computed directly from the position of control points and weights, 𝐏and 𝐰, corresponding to the NURBS mesh as: 𝐏𝐵=𝐖−𝟣 𝐵𝐂𝖳 𝐵𝐖𝐏,𝐰𝐵=𝐂𝖳 𝐵𝐰,(8) where matrices 𝐖𝐵and 𝐖are diagonal matrices storing the weight vectors 𝐰𝑏and 𝐰, for more details see29. In a typical 2D case, like the one represented in Figure 4 , Bézier analysis starts with an IGA mesh. Then, one can use the previous transformations to convert step-by-step the NURBS based mesh into a Bézier mesh by repeating knots in the two parametric directions. Finally, a finite element mesh is obtained with 𝐶0continuity of neighboring elements with Bernstein polynomials as shape functions. 3 HAMILTON’S PRINCIPLE WITH THE METHOD OF LOCALIZED LAGRANGE MULTIPLIERS We consider the evolution in time of the displacement and momentum fields 𝒖,𝒑∈ℝ𝑑, where 𝑑= 1,2,3is the number of space dimensions, of a body with domain Ωand boundary Γfor a total time 𝑇. The body is subjected to prescribed interface displacements 𝒖𝑏on part of its boundary Γ𝑏, to the traction 𝒕on the rest of its boundary Γ𝑡and to a body force 𝒃per unit volume. 6GONZÁLEZ ET AL FIGURE 4 Illustration of the Bézier extraction operation performed in a 2D patch of quadratic B-splines. The Bézier extraction operator maps a piecewise Bernstein polynomials basis (right) onto a B-spline basis (left). This transformation makes it possible to create Bézier elements as the finite element representation of B-splines or NURBS. Initial conditions that are compatible with the boundary conditions define the displacement and velocity field respectively at time 𝑡= 0. The initial boundary value problem of elastodynamics expressed in strong form can be formulated as follows: 𝛁⋅𝝈+𝒃= 𝒑in Ω × [0, 𝑇 ], 𝒖=𝒖𝑏on Γ𝑏×[0, 𝑇 ], 𝝈⋅𝒏=𝒕on Γ𝑡×[0, 𝑇 ], 𝒖|𝑡=0 =𝒖0in Ω,  𝒖|𝑡=0 = 𝒖0in Ω, (9) where 𝒑=𝜌 𝒖is the momentum per unit volume and 𝝈is the stress tensor that, in the case of linear elasticity is linearly related with the deformation tensor 𝜺=1 2[∇𝒖+ ∇𝖳𝒖]by the fourth-order constitutive tensor through the constitutive relation 𝝈=∶𝜺. This is an initial value problem to be solved for all time 𝑡in the interval [0, 𝑇 ]. The strong form of the of the initial value problem (9) can be transformed into an equivalent weak variational form. To this end, we invoke Hamilton’s principle for constrained elastodynamic problems. This principle states that the path followed by a dynamic system is the one which minimizes the action integral of the Lagrangian; condition that can be expressed using the following three-field variational form: 𝛿𝐻(𝒖,𝒑,𝓵) = 𝑡2 ∫ 𝑡1 𝛿{𝑇( 𝒖,𝒑) − 𝑈(𝒖,𝓵) + 𝑊(𝒖)}𝑑𝑡 = 0,(10) where 𝛿𝑇 is the virtual kinetic energy, 𝛿𝑈 the virtual elastic energy and 𝛿𝑊 the virtual work done by the external loads, magnitudes that can be expressed: 𝛿𝑇 =∫ Ω 𝛿(1 2𝒑⋅ 𝒖)𝑑Ω,(11) 𝛿𝑈 =∫ Ω 𝛿𝜺∶𝝈𝑑Ω + ∫ Γ𝑏 𝛿{𝓵⋅(𝒖−𝒖𝑏)}𝑑Γ,(12) 𝛿𝑊 =∫ Ω 𝛿𝒖⋅𝒃𝑑Ω,(13) GONZÁLEZ ET AL 7 with known displacement boundary conditions 𝒖𝑏imposed at interface Γby using a field of localized Lagrange multipliers 𝓵 and external body forces per unit volume 𝒃acting on Ω. By making use of the momentum-velocity relation (𝒑=𝜌 𝒖) in the following identity: 1 2𝒑⋅ 𝒖=𝒑⋅ 𝒖−1 2𝜌𝒑⋅𝒑(14) we can obtain an equivalent expression for the virtual kinetic energy per unit volume: 𝛿(1 2𝒑⋅ 𝒖)=𝛿 𝒖⋅𝒑+𝛿𝒑⋅( 𝒖−1 𝜌𝒑)(15) and finally performing an integration by parts of the second term of the last equation yields: 𝑡2 ∫ 𝑡1 𝛿 𝐮⋅𝒑𝑑𝑡 = − 𝑡2 ∫ 𝑡1 𝛿𝒖⋅ 𝒑𝑑𝑡, (16) relation that can be used to eliminate the virtual velocity field from the formulation. Introducing previous identities, (15) and (16), into the principle of stationary action (10), we obtain the final the three-field variational form of the Hamilton’s principle for constrained elastodynamics: 𝛿𝐻(𝒖,𝒑,𝓵) = 𝑡2 ∫ 𝑡1⎧ ⎪ ⎨ ⎪ ⎩∫ Ω 𝛿𝒑⋅( 𝒖−1 𝜌𝒑)𝑑Ω − ∫ Ω (𝛿𝒖⋅ 𝒑+𝛿𝜺∶𝝈)𝑑Ω + ∫ Ω 𝛿𝒖⋅𝒃𝑑Ω +∫ Γ𝑏 𝛿𝒖⋅𝓵𝑑Γ − ∫ Γ𝑏 𝛿𝓵⋅(𝒖−𝒖𝑏)𝑑Γ⎫ ⎪ ⎬ ⎪ ⎭ 𝑑𝑡 = 0,(17) expression that will be used to derive the discrete equations of motion. Remark 1. Through the use of a parametrized form of the kinematic energy with three free parameters (𝐶1, 𝐶2, 𝐶3), Schaeuble et al.21 proposed a general template for obtaining variationally scaled and customized consistent and reciprocal mass matrices. Variational equation (17) can be considered a particular case of their general template form for {𝐶1= −1, 𝐶2= 0, 𝐶3= 0}. An important difference with that work lies on the use of a volume-specific linear momentum instead of a mass-specific linear momentum and the consistent-lumped parametrization of the mass matrix that will be introduced later. 4 CONSTRUCTING INVERSE MASS MATRIX VIA THE LOCALIZED LAGRANGE MULTIPLIERS Discretization in space of the mixed variational form (17) is performed by using independent interpolation functions for displacements, momenta and Lagrangian multipliers. For this purpose, we carry out a standard IGA discretization (5) with independent B-spline shape functions for the three fields: 𝒖(𝝃) = 𝐍𝑢(𝝃)𝐮,𝒑(𝝃) = 𝐍𝑝(𝝃)𝐩,𝓵(𝝃) = 𝐍𝜆(𝝃)𝝀,(18) and introduce these discrete approximations in the variational form (17) to obtain, after imposing the stationarity condition, the following set of semi-discrete equations: 𝐀 𝐩+𝐁𝝀=𝐫Equilibrium equation,(19) 𝐀𝖳 𝐮−𝐂𝐩 =𝟎Momentum equation,(20) 𝐁𝖳𝐮−𝐋𝑏𝐮𝑏=𝟎Boundary (and interface) constraints,(21) −𝐋𝖳 𝑏𝝀=𝟎Newton’s 3rd law on the boundaries,(22) where vector 𝐫=𝐟−𝐟𝑖𝑛𝑡 is the external-internal forces residual, 𝐀is the global projection matrix, 𝐂the global reciprocal mass matrix, 𝐁the boundary assembly operator and 𝐋𝑏the Localized multipliers assembly matrix. These matrices are constituted by 8GONZÁLEZ ET AL an element-by-element assembly of the element matrices: 𝐀𝑒=∫ Ω𝑒 𝐍𝖳 𝑢𝐍𝑝𝑑Ω,𝐂𝑒=∫ Ω𝑒 1 𝜌𝐍𝖳 𝑝𝐍𝑝𝑑Ω,(23) 𝐁𝑒=∫ Γ𝑒 𝐍𝖳 𝑢𝐍𝜆𝑑Γ,𝐋𝑒=∫ Γ𝑒 𝐍𝖳 𝜆𝐍𝑢𝑑Γ,(24) the first two needed for the construction of the element reciprocal mass and the last ones used to enforce interface compatibility and equilibrium. Moreover it is important to mention that we assume constant mass density along every physical patch. Going one step forward by eliminating symbolically the momentum variable 𝐩from (19) and (20), one obtains the classical equation of motion expressed in terms of displacements: (𝐀𝐂−𝟣𝐀𝖳) 𝐮+𝐁𝝀=𝐫(25) with the mass matrix approximated as: 𝐌=𝐀𝐂−𝟣𝐀𝖳(26) and observe that there must exist an inverse mass matrix (denoted by equivalence as 𝐌−𝟣) given by: 𝐌−𝟣=𝐀−𝖳𝐂𝐀−𝟣(27) assuming that the global projection matrix is invertible. Since the objective of the present paper is to obtain inverse mass-matrices for IGA in efficient and accurate ways, one must seek for easily to invert diagonal or narrowly banded projection matrices. This certainly can be accomplished in (27) with a diagonal projection 𝐀-matrix, leaving only the task of constructing the reciprocal 𝐂-matrix efficiently and accurately. 5 DUAL BASIS FUNCTIONS Observing again the expression obtained for the mass matrix inverse (27), we realize that a careful selection of the basis functions can simplify its calculation. Dual bases can be local or global,continuous or discontinuous and satisfy exactly the orthogonality condition or not. In particular, by using displacement and momentum bi-orthogonal interpolation functions, we are able to construct a diagonal projection matrix in the form: [𝐀𝑒]𝑖𝑗 =∫ Ω𝑒 𝐍𝖳 𝑢𝑖𝐍𝑝𝑗𝑑Ω = 𝑑𝑖𝛿𝑖𝑗 ,where[𝐀𝑒]𝑖𝑗 = 0 for 𝑖≠𝑗(28) with a trivial inverse that permits to compute the inverse mass matrix directly from the momentum matrix. Different local and global dual interpolation functions of this kind, originally designed and tested for mortar methods, have been recently proposed in the literature30,31,32 for NURBS and B-spline interpolations. Specifically for the construction of reciprocal mass matrices, biorthogonal basis functions have been proposed under the framework of FEM8,9 and also IGA21. The different options available for the construction of dual bases are discussed next. 5.1 Local dual bases These dual basis functions entail minimal support, i.e., the support of the dual basis is equal to the support of the original B- spline basis, that for the case of B-splines covers 𝑝+1 knot spans. Analytical expressions for the construction of local continuous (LC) bi-orthogonal dual basis exist in the literature for one-dimensional B-splines33,34,32, but they are difficult to generalize for multidimensions, restricting their application to our case. For higher dimensions, the simplest alternative is the use of a local discontinuous (LD) bi-orthogonal dual basis35,31 obtained via an element-wise linear combination of the primal shape functions, providing discontinuous reciprocal basis functions that fulfill the partition of unity with the same support than their primal counterparts. A dual basis 𝐍𝑝(𝝃)of this type, orthogonal to the B-spline displacement interpolation functions 𝐍𝑢(𝝃), can be written in matrix form as: 𝐍𝑝(𝝃) = 𝐍𝑢(𝝃)𝐌−𝟣 𝑒𝐌𝐿 𝑒(29) GONZÁLEZ ET AL 15 Similarly, for quadratic B-splines, with element matrices given by (47) and (48), the Taylor expansion with respect to wavenumber 𝜅yields the dispersion relation, Ω2=𝜅2+3−7𝛽 4(1 + 3𝛽)𝜅4+(𝜅6)(54) where it is observed that fourth order approximation can be achieved by using 𝛽= 3∕7 and the admissible interval of the averaging parameter is extended to 𝛽∈ [3∕7,1]. Moreover, an optimization of parameter 𝛽for higher frequencies, see Figure 9 , reveals that the best approximation is obtained for 𝛽= 2∕5.Using 𝛽= 2∕5 is a compromise to improve accuracy for the higher frequencies, maintaining second order approximation for the lower frequencies. This value has been obtained empirically, minimizing by numerical experiment the dispersion error at 𝜅=𝜋∕2. Remark 4. For B-splines, the RMM obtained with an element-wise LD dual-basis is very cheap to compute but should be used with care. It diverges for 𝛽= 0 (CMM) and is not able to provide accurate results. However, the use of mass parametrization with 𝛽∈ [2∕5,1) always provides a better approximation than the LMM diagonal mass inverse, although p-convergence property is not preserved. It is worth to mention that this type of LD B-spline dual basis has also been studied mathematically in the context of IGA mortar methods for contact problems by Steiz et al.31, observing that LD bases generated element-wise only guarantee partition of unity and consequently the global approximation order is limited to one in the 𝐿2-norm, independently of the local approximation order. This defect in the approximation order of the B-spline dual basis causes that optimal convergence rates cannot be guaranteed. Obviously, our element-wise dual basis reproduce the same convergence problems when utilized for inverse mass matrix computations, but these difficulties can be palliated with mass matrix parametrization techniques21.Finally, it has been also investigated the use of sixth-order accurate element mass matrices42,43,44 instead of the fourth order approximation given by (43) for 𝛽= 0. By adopting these high-order mass matrices, it is possible to increase the accuracy of the RMM, but the improvement is still not enough for its practical use with LD dual bases. Remark 5. Because the reciprocal mass matrix 𝐂is by definition positive definite, the unsymmetric eigenvalue problem [𝐀−𝖳𝐂𝐀−1𝐊−𝜔2 𝑖𝐈]𝝓𝑖= 0 represented by equation (51) can be converted into a symmetric eigenvalue problem by a Cholesky decompositionofthe reciprocal mass matrixwith 𝐂=𝐋𝐋𝖳, where 𝐋is a lower triangular matrix.Then,the equivalentsymmetric eigenvalue problem is obtained: [𝐋𝖳𝐀−1𝐊𝐀−𝖳𝐋−𝜔2 𝑖𝐈]𝝋𝑖= 0 (55) which has different eigenvectors but the same eigenvalues as the original unsymmetric problem. 7 INCORPORATION OF BOUNDARY CONDITIONS WITH LAGRANGE MULTIPLIERS As we mentioned earlier, in the connection with other patches or presence of Dirichlet boundary conditions a direct computation of the mass-inverse requires complicated modifications of the bi-orthogonal shape functions. We eliminate the necessity of these modifications by considering the patches as free-floating substructures and applying the boundary conditions through localized Lagrangian multipliers10,11,12,13. To explain the enforcement of boundary conditions, we reorganize the semi-discrete equations of motion (20)-(22) to obtain the following partitioned equation set: ⎡⎢⎢⎣ 𝐌 𝐁 𝟎 𝐁𝖳𝟎−𝐋𝑏 𝟎−𝐋𝖳 𝑏𝟎⎤⎥⎥⎦⎧ ⎪ ⎨ ⎪ ⎩  𝐮 𝝀  𝐮𝑏⎫ ⎪ ⎬ ⎪ ⎭ =⎧ ⎪ ⎨ ⎪ ⎩ 𝐫 𝟎 𝟎⎫ ⎪ ⎬ ⎪ ⎭ (56) where the second equation is simply twice time-differentiated expression of (21). Let’s assume now that the boundary acceleration vector ( 𝐮𝑏) consists of entirely the prescribed boundary quantities, regardless whether they are zero or nonzero. In this case, we obtain the accelerations ( 𝐮) from the first equation of (56) as:  𝐮=𝐌−𝟣(𝐫−𝐁𝝀)(57) in which once again we emphasize that the computation of 𝐌−𝟣=𝐀−𝖳𝐂𝐀−𝟣does not requiere additional computational effort because the global projection matrix 𝐀is diagonal. 16 GONZÁLEZ ET AL Then, substitute the accelerations obtained above into the second equation of (56): 𝐁𝖳𝐌−𝟣(𝐫−𝐁𝝀)−𝐋𝑏 𝐮𝑏=𝟎(58) to compute the Lagrange multipliers: 𝝀=[𝐁𝖳𝐌−𝟣𝐁]−𝟣(𝐁𝖳𝐌−𝟣𝐫−𝐋𝑏 𝐮𝑏)(59) as a function of the forces residual and the boundary conditions. Finally, substitute the multipliers back in the equilibrium equation (57) to obtain the final accelerations:  𝐮=𝐌−𝟣 𝑏𝐫+𝐌−𝟣𝐁[𝐁𝑇𝐌−𝟣𝐁]−𝟣𝐋𝑏 𝐮𝑏(60) with a constrained inverse mass matrix: 𝐌−𝟣 𝑏=𝐌−𝟣−𝐌−𝟣𝐁[𝐁𝖳𝐌−𝟣𝐁]−𝟣𝐁𝖳𝐌−𝟣(61) incorporating the applied boundary conditions. This transformation can also be understood as a projection of the free-floating mass matrix inverse9. Once again, since the size of [𝐁𝑇𝐌−𝟣𝐁]is small, pertaining only to the constrained degrees of freedom, its inversion is very cheap. It is also noted that  𝐮contains not only the interior but also the constrained degrees of freedom, and the imposition of boundary conditions will produce zeros in the associated rows and columns of 𝐌−𝟣 𝑏making it singular. These zero rows and columns, corresponding to restricted degrees of freedom, can be simply eliminated. Algorithm 1 Procedure for the computation of the constrained inverse mass matrix 𝐌−1 𝑏 for edo 1...𝑁𝑒(elements) Compute parametrized element mass matrix: 𝐌𝑒= (1 − 𝛽)𝐌𝐶 𝑒+𝛽𝐌𝐿 𝑒 Assemble diagonal projection matrix: 𝐀𝑒=1 𝜌𝐌𝐿 𝑒←→ 𝐀= 𝑁𝑒 A 𝑒=1 𝐀𝑒 Assemble reciprocal mass matrix: 𝐂𝑒=𝐀𝖳 𝑒𝐌−1 𝑒𝐀𝑒←→ 𝐂= 𝑁𝑒 A 𝑒=1 𝐂𝑒 end for Compute the free-floating inverse mass matrix: 𝐌−1 =𝐀−𝖳𝐂𝐀−1 Obtain the projector: =𝐈−𝐌−1𝐁[𝐁𝖳𝐌−1𝐁]−1 𝐁𝖳 Apply boundary conditions by projection: 𝐌−1 𝑏=𝐌−1 Eliminate the rows and columns of 𝐌−1 𝑏with applied boundary conditions 8 NUMERICAL TESTS In the following examples, for vibration analyses we seek for an harmonic solutions that satisfy the EOM and fulfill the boundary conditions 𝐮𝑏=𝟎. These solutions are obtained by solving the following eigenvalue problem: [𝐌−1 𝑏𝐊𝑏−𝜔2𝐈]𝝓=𝟎(62) where 𝐌−1 𝑏is the inverse mass matrix computed as described in Section 7 and 𝐊𝑏is the stiffness matrix incorporating the applied boundary conditions. From the solution of this eigenvalue problem, we obtain the vibration frequencies 𝜔𝑖of the system and its associated mode shapes 𝝓𝑖for mode number 𝑖= 1,…𝑛. GONZÁLEZ ET AL 17 In time-domain analyses, for time integration of the semi-discrete equations of elastodynamics, we use the classical explicit central difference scheme45,46,47. This method is among the most popular explicit methods in computational mechanics4, based on a repetition of the basic steps summarized in Algorithm 2 until the total simulation-time has been reached. Algorithm 2 Explicit time integration by central difference scheme Solve time 𝑡= 0:  𝐮0=𝐌−1 𝑏𝐫0 for ndo 1...𝑛 (time steps) Evaluate the forces residual: 𝐫𝑛=𝐟𝑛−𝐊𝑏𝐮𝑛 Compute accelerations:  𝐮𝑛=𝐌−1 𝑏𝐫𝑛 Update nodal velocities:  𝐮𝑛+1 2= 𝐮𝑛−1 2+ Δ𝑡 𝐮𝑛 Update nodal displacements: 𝐮𝑛+1 =𝐮𝑛+ Δ𝑡 𝐮𝑛+1 2 end for It is well known that the explicit central difference scheme in linear dynamics is conditionally stable, so the time-step Δ𝑡must be chosen below a critical value Δ𝑡𝑐𝑟𝑖𝑡 = 2∕𝜔𝑚𝑎𝑥, where 𝜔𝑚𝑎𝑥 is the maximum frequency of the system. We have seen that the proposed RMM with parametrization, obtained following the process described in Section 4, is able to retain accuracy the for lower frequencies reducing highest frequencies at the same time. This reduction automatically increases the critical time-step, allowing the use of larger time steps and improving the computational efficiency. To evaluate all element matrices we have used a Gauss-Legendre quadrature formula1of degree 𝑝+ 1 over the knot intervals, where 𝑝is the order of the used B-spline or Bézier interpolation function. For the presentation of results we have adopted the following nomenclature. Plots with labels IGA and IGA-BE are used as a reference, indicating results obtained with the exact IGA mass matrix inverse, computed respectively with B-splines or Bézier extraction. The proposed reciprocal matrix is referred as RMM and RMM-BE depending on the interpolation functions, B-splines or Bézier elements, respectively. 8.1 Free vibration of a bar In this first example, we analyze the accuracy of the proposed method for direct computation or reciprocal mass matrices by analyzing the free vibration of elastic bars. For this purpose we consider we consider a long bar of length 𝐿= 1 m of uniform section. The material of the bar is linear elastic, with Young’s modulus 𝐸= 69 GPa and constant volumetric mass density 𝜌= 2700 kg∕m3. The bar is discretized using a linear parameterization of 43 control points to investigate the accuracy of the RMM for quadratic and cubic 𝑝= 2,3orders of longitudinal displacement discretization. In particular, we investigate differences in the use of exact IGA inverses of the CMM and LMM matrices and the proposed approximation of the mass matrix inverse RMM computed with an optimized 𝛽. 8.1.1 Free-free case We start with the free-free case, by comparing the numerical frequencies with exact solutions obtained from continuum theory. The exact frequencies of this problem are 𝜔𝑒𝑥𝑎𝑐𝑡 𝑛= (𝑛− 1)𝑐0𝜋∕𝐿rad∕s for mode number 𝑛= 1,2,…where 𝑐0=√𝐸∕𝜌is the wave speed of compressive waves and 𝑛= 0 provides a zero eigenfrequency corresponding to the free rigid body motion. Results for this case are presented in Figure 10 for two orders of approximation, quadratic 𝑝= 2 on the left and cubic 𝑝= 3 on the right, and two different discretization methods, B-splines on the top and Bézier extraction at the bottom. We include frequencies for the exact inverse of the IGA mass matrix, for CMM 𝛽= 0, optimized 𝛽= 2∕5 and LLM 𝛽= 1 cases. For comparison, the RMM computed with the proposed method and optimized parameter 𝛽= 2∕5 is also included. 18 GONZÁLEZ ET AL 5 10 15 20 25 30 35 40 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 5 10 15 20 25 30 35 40 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 5 10 15 20 25 30 35 40 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 5 10 15 20 25 30 35 40 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 FIGURE 10 Free-free bar vibration problem computed with LD dual basis. Error in the frequencies for quadratic 𝑝= 2 (left) and cubic 𝑝= 3 (right) orders using B-splines (top) and Bézier extraction discretizations (bottom). The RMM with 𝛽= 0 (CMM) produces divergent results for B-splines and p-convergent behaviour for Bézier elements. In all cases, RMM with an optimized parameter 𝛽= 2∕5 are able to capture correctly the low frequency band providing higher accuracy than LMM (𝛽= 1) for the full frequency spectrum. It is well known27,48 and evident from inspection of Figure 10 (top) that, by increasing the order 𝑝of the LMM inverse, higher-order accuracy is not achieved. We also observe that the present method, using the RMM with 𝛽= 2∕5, yields far better accuracy than the LMM for all frequencies. With quadratic B-splines, the optimized RMM is able to approximate almost half of the spectrum with good accuracy, first 20 modes are below 5% error while for cubic order it is reduced to the 15 first modes approximately. But in general, increasing the approximation order either for B-spline or Bézier approximations, expands the nonzero elements of the banded mass matrix and reduces the accuracy in the mid and high frequency range. The RMM for Bézier elements, on the contrary, is completely stable for 𝛽= 0 conserving the p-convergence property, see Figure 10 (bottom) for cases p=2,3. It can also be observed that although Bézier extraction introduces higher branching in the spectrum49, its accuracy is very close to the B-splines approximation. This branching is caused by the two distinct types of equations for the Bézier based elements: those corresponding to nodes located at the element boundaries, and those corresponding to the internal nodes located inside the element27. Note that the number of branches coincide with the Bézier element degree. Moreover, RMM produces a maximum frequency around 40% lower than the exact IGA consistent matrix (IGA 𝛽= 0), so the critical time-step for explicit time integration with the RMM can be reduced by the same amount. This is close to the 55% reduction of the diagonal LMM but without so much degradation of the accuracy for the low frequencies. GONZÁLEZ ET AL 19 5 10 15 20 25 30 35 40 1 1.02 1.04 1.06 1.08 1.1 1.12 1.14 5 10152025303540 1 1.02 1.04 1.06 1.08 1.1 1.12 1.14 FIGURE 11 Error in the natural frequencies of a free-free bar computed using GC dual bases. The inverse Gram matrix (left) and the weight-adjusted (GC-WA) basis (right) produce p-convergent approximations. These GC bases provide high accuracy, but at the price of elevated computational cost in the case of the Gram matrix and loss of generality for the weight-adjusted basis. Finally, we explore the accuracy of GC dual bases in the bar frequencies. In particular we use the inverse Gram matrix and the weight-adjusted to compute the RMM for different B-spline basis orders. The results are shown in Figure 11 , where it is observed that the frequency error is totally equivalent to the IGA mass matrix and p-converge is maintained. Once again, note that the evaluation of these two reciprocal bases require two costly operations, invert the Gramian matrix (left) and invert the mass matrix of the parametric domain (right), but nevertheless these operations are cheaper than the inversion of the original mass matrix. 8.1.2 Fixed-free bar Now, to test the proposed technique for applying boundary conditions using localized Lagrange multipliers described in Section 7, we proceed to fix the displacement of the first control point located on the left side of the bar. In terms of frequencies, the only difference with the free-free case is in the restriction of the rigid body motion of the bar and the shift of exact frequencies to 𝜔𝑒𝑥𝑎𝑐𝑡 𝑛= (𝑛− 1∕2)𝑐0𝜋∕𝐿for mode number 𝑛= 1,2,…. Figure 12 contains the results obtained for this case. We corroborate that accuracy is not affected by the imposition of boundary conditions through localized Lagrange multipliers and the results are completely equivalent. In this case we have tested the RMM with 𝛽= 2∕5 that provides very similar approximation of the bar frequencies. The results confirm that using the present method to approximate the inverse mass matrix yields much better accuracy than the LMM with similar computational effort. 8.2 Free vibration of a cantilevered tapered plate - FV32 test This example is a plane-stress test proposed by the National Agency for Finite Element Methods and Standards50. In particular the dynamic case FV32 for checking the accuracy of flexural and longitudinal modes, where the natural frequencies of a trapezoidal elastic plate are analyzed. The geometry and dimensions of the plate are represented in Figure 13 together with the boundary conditions. It is assumed a linear elastic material with Young’s modulus 𝐸= 200 GPa, Poisson’s ratio 𝜈= 0.3, and density 𝜌= 8000 kg∕m3. First, the plate is discretized using a single patch with a regular mesh of 10 × 10 control points and quadratic isogeometric NURBS-based elements, constraining the horizontal and vertical displacements of the 10 control points located at 𝑥= 0, a lateral side of the plate. The knot points vector Ξ = {0,0,0,0.125,0.25,0.375,0.5,0.625,0.75,0.875,1,1,1}is used in each direction to construct a regular 2D patch Ξ×Ξwith B-splines of order 𝑝= 2. 20 GONZÁLEZ ET AL 5 10 15 20 25 30 35 40 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 5 10 15 20 25 30 35 40 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 5 10 15 20 25 30 35 40 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 5 10 15 20 25 30 35 40 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 FIGURE 12 Fixed-free bar problem computed with LD dual basis. Error in the frequencies for quadratic 𝑝= 2 (left) and cubic 𝑝= 3 (right) orders using B-splines (top) and Bézier extraction discretizations (bottom). We observe that RMM with 𝛽= 0 (CMM) produces divergent results for B-splines and p-convergent behaviour for Bézier elements. In all cases, RMM with an optimized parameter 𝛽= 2∕5 are able to capture correctly the low frequency band providing higher accuracy than LMM (𝛽= 1) for the full frequency spectrum.. Next, the Bézier extraction process is applied to the patch of quadratic B-splines, see Figure 14 , producing a mesh of 5 × 5 Bézier elements with 11×11 nodes. Note that in both cases, these are more refined meshes than the one proposed in the original NAFEMS benchmark case, that uses quadrilateral finite elements. This is a particular example where the weight-adjusted dual basis (34) is readily applicable, because the parametric domain is completely regular and the parametrical-to-physical geometrical transformation is an affine mapping, hence Jacobian of the transformation is constant. The first 9 mode shapes of the plate are illustrated in Figure 15 and their associated frequencies, obtained with different mass matrix inverses, arranged in Table 1 . As a reference, we use the frequencies computed with the original mass matrix, IGA or IGA-BE (CMM). These frequencies are very close to the NAFEMS test for FEM and results from other authors obtained with different RMM proposals. First, we note that GC-Gram and GC-WA bases produce exactly the same frequencies and are able to approximate the correct values with negligible error, confirming once again the superiority of these GC bases over the LD ones. Nevertheless, the LDRMM computed with 𝛽= 2∕5 provides a very good fit, with maximum difference obtained respectively in the fourth mode with 1.8% error and the fifth mode with 4.9% error. Note that they correspond to flexural modes that are more difficult to capture GONZÁLEZ ET AL 21 5 10 1 y x FIGURE 13 NAFEMS benchmark case FV32. Geometry with dimensions in meters and zero boundary displacements imposed on the lateral side of the plate at 𝑥= 0. The IGA regular mesh presents 𝑙×𝑚= 10 × 10 control points, indicated in blue and uses quadratic B-spline interpolation functions in both directions. FIGURE 14 A computational IGA mesh based on B-spline shape functions (on the left) and based on Beziér shape functions after Beziér extrapolation (on the right). than extension modes, see Figure 15 . For Bézier elements, maximum errors correspond to the same flexural modes with values under 8%. In general, it is verified that the proposed RMM obtained with LD bases and mass-matrix parametrization is a better alternative than the LMM for all modes. 8.3 Free vibration of a clamped thin circular plate This numerical example is inspired by the paper of Cottrell et al.22, where the problem of free vibration of a thin circular plate is solved. The purpose of this example is to show the behavior of the proposed method also for a domain with a curved boundary. The circular plate has a radius 𝑟= 2 m and thickness 𝑡= 0.02 mm. Although it would seem natural to model the plate using shell finite elements, we follow the paper22 and model the circular plate as a three-dimensional body. Material of the plate is considered to be linear elastic with Young’s modulus 𝐸= 30 GPa, Poisson’s ratio 𝜈= 0.2, and density 𝜌= 2320 kg∕m3. Similar to the preceding example, the plate is first discretized using a single patch with a regular mesh of 18×18×3control points and quadratic isogeometric NURBS-based elements, clamping of the control points located at 𝑥2+𝑦2=𝑟2as represented in Figure 16 . The knot vectors Ξ1= Ξ2={0,0,0,1∕16,2∕16,…,15∕16,1,1,1}and Ξ3={0,0,0,1,1,1}are used in each parametric direction to construct a 3D patch Ξ1×Ξ2×Ξ3with NURBS of order 𝑝= 2. Details regarding the used parameterization can be found in the work51. Again, the Bézier extraction process is applied to the patch of quadratic B-splines, see Figure 17 , producing a mesh of 8×8×1 Bézier elements with 17 × 17 × 3 control points. 22 GONZÁLEZ ET AL Mode 1 Mode 2 Mode 3 Mode 4 Mode 5 Mode 6 Mode 7 Mode 8 Mode 9 FIGURE 15 NAFEMS benchmark case FV32. Representation of the 9 lowest vibration modes occurring at frequencies contained in Table 1 . Results correspond to the proposed reciprocal mass matrix with parameter 𝛽= 2∕5 computed with quadratic B-splines. The three first symmetric 𝜔01,𝜔11 and antisymmetric 𝜔02 eigenmodes22 are illustrated in Figure 18 , and their associated frequencies, obtained with different mass matrix inverses, summarized in Table 2 . Classical IGA with CMM gives the first three eigenmodes at frequencies 𝜔01 = 57.636 rad/s, 𝜔11 = 127.061 rad/s, and 𝜔02 = 260.913 rad/s that are used as a reference. Note that GC-Gram basis produces again a perfect approximation. As predicted by dispersion analysis, the LD-RMM with best fit is obtained for 𝛽= 2∕5, see IGA RMM results in Table 2 , presenting discrepancies of 0.05%, 0.04% and 0.2% respectively. Observing now the IGE-BE (CMM) solution, a slightly poorer approximation is obtained for the Bézier extraction results, compared to IGA with NURBS basis functions. The average relative errors for the Bézier extraction for eigenfrequencies 𝜔01, 𝜔11, and 𝜔02 are 0.3%, 0.8% and 2%, respectively. It is important to mention again that the accuracy of the results for LC-RMM crucially depends on the selection of an appropriated value for the parametrization factor 𝛽. Its value has been derived from one dimensional dispersion analyses and its optimization for 2D cases deserves further study. 8.4 Transient analysis of a cantilever beam under self weight Finally, to test the RMM with a transient case, we consider a 3D cantilever beam under self weight with a sudden application of the gravity load at time 𝑡= 0 s. The geometrical definition of the problem is shown in Figure 19 (left), where a cantilever beam of length 𝐿= 10 m and section width 𝑏= 4 m with height ℎ= 2 m is considered. The elastic material properties of the beam are 𝐸= 10 MPa, 𝜈= 0.3and density 𝜌= 100 kg∕m3with a gravity acceleration 𝑔= 9.81 m∕s2. The parametric domain is defined by a 3D patch constructed with knot spans on each direction Ξ1={𝜉1,…, 𝜉𝑙+𝑝+1}with Δ𝜉= 1∕19,Ξ2={𝜂1,…, 𝜂𝑚+𝑝+1}with Δ𝜂= 1∕5 and Ξ3={𝜇1,…, 𝜇𝑛+𝑝+1}with Δ𝜇= 1∕3. Quadratic B-splines with 𝑝= 2 are used in the three directions. Elements are defined in the parametric space Ξ1× Ξ2× Ξ3as non-zero knot spans, [𝜉𝑖, 𝜉𝑖+1]×[𝜂𝑖, 𝜂𝑖+1]×[𝜇𝑖, 𝜇𝑖+1]and elements in the physical space are images of their parametric counterparts. GONZÁLEZ ET AL 23 TABLE 1 Plane-stress case FV-32 from NAFEMS. Six lowest frequencies obtained with different mass matrices. Frequency (Hz) Mode 1 Mode 2 Mode 3 Mode 4 Mode 5 Mode 6 NAFEMS50 (FEM) 44.623 130.03 162.70 246.05 379.90 391.44 Tkachuk et al.8(FEM) 45.499 133.55 162.89 255.37 393.16 398.56 Schaeuble et al.21 (IGA) VSRMS†44.615 128.489 162.860 235.104 335.874 385.989 IGA 𝛽= 0 (CMM) 44.634 130.107 162.709 246.639 382.937 391.525 IGA 𝛽= 1 (LMM) 44.022 122.936 161.413 218.492 308.210 373.077 LD-RMM 𝛽= 2∕5 44.645 129.581 162.720 241.956 363.728 390.592 GC-Gram 44.634 130.107 162.709 246.639 382.937 391.525 GC-WA 44.634 130.107 162.709 246.641 382.961 391.526 IGA-BE 𝛽= 0 (CMM) 44.651 130.343 162.730 248.396 389.846 391.848 IGA-BE 𝛽= 1 (LMM) 44.096 123.304 161.532 219.363 309.573 374.468 LD-RMM-BE 𝛽= 2∕5 44.544 128.531 162.489 238.999 358.320 387.791 The cited references use different meshes and methods. †Variationally scaled consistent inverse mass matrices. 4 y x FIGURE 16 Clamped thin circular plate. Geometry with diameter in meters and zero boundary displacements imposed on the lateral side of the plate at 𝑥2+𝑦2= 4. The IGA mesh presents 𝑙×𝑚×𝑛= 18 × 18 × 3 control points, indicated in blue and uses quadratic NURBS interpolation functions in all three directions. This patch describes the physical domain with a mesh of 𝑙×𝑚×𝑛= 21 × 7 × 5 = 735 control points, producing a total number of 735 × 3 = 2205 unknowns except for those control points located at the plane 𝑥= 0 that are restricted in the three directions of space. An exact eigenvalue analysis of the problem with a CMM reveals that the maximum frequency of the structure is 𝜔𝑚𝑎𝑥 = 3510.5 rad∕s, hence the critical time step for stability of the central difference scheme is Δ𝑡𝐼𝐺𝐴 𝑐𝑟𝑖𝑡 = 5.7 × 10−4 s. This critical time step is utilized to compute the reference solution of the problem. The use of RMM with 𝛽= 2∕5 reduces the maximum frequency to 𝜔𝑚𝑎𝑥 = 1554.8 rad∕s, increasing the critical time step to Δ𝑡𝑅𝑀𝑀 𝑐𝑟𝑖𝑡 = 0.0013 s, i.e., 2.3times larger than the initial critical time step. Then, a time step Δ𝑡= 0.001 s is selected to perform time-domain analysis, integrating in time the response of the beam with the explicit central difference method described in Algorithm 2. It is observed in Figure 19 (right) that this large time-step Δ𝑡 ≫ Δ𝑡𝐼𝐺𝐴 𝑐𝑟𝑖𝑡 produces unstable results with classical IGA, while the RMM with Δ𝑡≤Δ𝑡𝑅𝑀𝑀 𝑐𝑟𝑖𝑡 is able to obtain 24 GONZÁLEZ ET AL FIGURE 17 A computational IGA mesh based on NURBS basis functions (on the left) and based on Beziér shape functions after Beziér extrapolation (on the right). TABLE 2 Clamped thin circular plate. Three lowest eigenfrequencies obtained with different mass matrices. IGA Bézier extraction Frequency (rad/s) 𝜔01 𝜔11 𝜔02 𝜔01 𝜔11 𝜔02 IGA 𝛽= 0 (CMM) 57.636 127.061 260.913 55.933 119.250 231.591 IGA 𝛽= 1 (LMM) 56.947 123.291 245.604 55.022 114.317 211.815 LD-RMM 𝛽= 2∕5 57.666 127.114 260.358 55.774 118.274 226.896 GC-Gram 57.636 127.061 260.914 - - - Eigenmode ω01 Eigenmode ω11 Eigenmode ω02 FIGURE 18 Clamped thin circular plate. Two symmetric and one antisymmetric vibration modes corresponding respectively to frequencies 𝜔01,𝜔02 and 𝜔11. a stable and very accurate solution. In all the cases, initial conditions of the beam are assumed to be zero displacements and velocities, applying the full gravity load at the initial time-step. The results are shown in Figures 19 and 20 for NURBS and Figure 21 for Bézier elements. We have selected for comparison the RMM computed with different values of the parameter 𝛽={0,2∕5,1∕2,1}, representing the evolution of the tip deflection of the beam with both approximation methods in Figures 19 (right) and 21 (left). We also compare the time evolution of total internal 𝑈, external 𝑊and kinetic 𝑇energies. Finally, in Figure 20 (right) it is represented the relative error (𝑢−𝑢𝑒𝑥𝑎𝑐𝑡)∕|𝑢𝑒𝑥𝑎𝑐𝑡|of the tip displacement 𝑢, computed with the proposed RMM (𝛽= 2∕5). It is observed that the error remains under 0.12% during the 1.5s of simulation, indicating a very good agreement with the reference solution.