Reciprocal mass matrices and a feasible time step estimator for finite elements with Allman's rotations

Finite elements with Allman's rotations provide good computational efficiency for explicit codes exhibiting less locking than linear elements and lower computational cost than quadratic finite elements. One way to further raise their efficiency is to increase the feasible time step or increase the accuracy of the lowest eigenfrequencies via reciprocal mass matrices. This article presents a formulation for variationally scaled reciprocal mass matrices and an efficient estimator for the feasible time step for finite elements with Allman's rotations. These developments take special care of two core features of such elements: existence of spurious zero‐energy rotation modes implying the incompleteness of the ansatz spaces, and the presence of mixed‐dimensional degrees of freedom. The former feature excludes construction of dual bases used in the standard variational derivation of reciprocal mass matrices. The latter feature destroys the efficiency of the existing nodal‐based time step estimators stemming from the Gershgorin's eigenvalue bound. Finally, the developments are tested for standard benchmarks and triangular, quadrilateral, and tetrahedral finite elements with Allman's rotations.


INTRODUCTION
Successful finite element formulations for solid continuum with vertex rotations are known from the early 1980s from the works by Mohr 1 and by Allman. 2 These vertex rotations are usually called as Allman's rotational degrees of freedom or simply Allman's rotations. The aim of the Allman's rotations is a compatible interpolation of the displacements for an element with only vertex degrees of freedom. Taking a standard quadratic element as a base, the midside nodal displacements are then constrained to vertex displacements and rotations of the edge by a transformation matrix. Therefore, the ansatz space of this family contains truncated quadratic functions. The common argument for the development of such formulations is an increased accuracy with respect to the first-order finite elements that tend to lock and a lower computational cost for comparable accuracy with respect to the standard quadratic finite elements.
A vast literature is available on the topic of finite elements with Allman's rotations even if the connection of Allman's rotation to drilling stiffness of plate elements is excluded. Some selected papers are given here ordered according to the finite element shape. These developments feature compatible strain field, assumed strain field, hybrid stress, and free formulations. The compatible strain field formulations rely on the strain-displacement operator from the compatible displacement field in the principle of virtual work, as initially proposed by Mohr. 1 The compatible strain field formulations tend to be too stiff and they could be refined by a geometrical projection of the compatible strains onto an assumed strain field. 2 Hybrid stress formulations are usually based on the standard version of Hellinger-Reissner variational principle 3,4 or its extension to drilling rotation field by Brezzi and Hughes. 5 The latter principle allows an easier suppression of spurious zero-energy modes (SZEM). The spurious zero energy modes are a common feature of the elements, for example, a vector with zero displacement and equal vertex rotations produce the zero displacement field for each 2D element. This SZEM can be easily suppressed by fixing of one rotation in a part. A 4-node tetrahedral element has four SZEM, which is equal to the rank deficit of the transformation matrix. 6 A targeted stiffness stabilization is employed to suppress them. Finally, free formulations construct the stiffness matrix of the element using an energy-orthogonal decomposition in so-called basic stiffness that is needed to pass the patch test and a higher order stiffness, which takes care of rank sufficiency and accuracy. 7 Free parameters in the higher order stiffness give raise to a stiffness template for a finite element and give the name of the free formulations. All these formulations for 3-node triangular elements for 2D plane elasticity are presented in References 1-3,7. Their comparison from the point of view of the stiffness templates and an optimal 3-node triangle finite elements are given by Felippa in Reference 8. A hybrid stress formulation for a 4-node quadrilateral element for 2D plane elasticity is presented in Reference 4 and it is further refined with better stress ansatz spaces and bubble functions in References 9,10. A 4-node tetrahedral element is proposed in a compatible strain field form and an assumed strain form in Reference 11. That element exhibits a spurious zero-energy mode and the assumed strain form is not objective. 6,12 Two objective hybrid stress versions of the tetrahedral element and stabilization with a sufficient rank are given in Reference 6. Large displacement inelastic problems can be solved with a co-rotational hybrid formulation. 13 Alternatively, a version based on the free formulation is presented in Reference 14. Finally, a compatible strain field and its refinement with three incompatible modes for an 8-node hexahedral element are presented in Reference 15. Further effort was made to construct mass matrices and test these formulations for modal and transient analyses. The consistent mass matrix (CMM) are presented in References 16,17 for a 3-node triangle, in Reference 18 for a 4-node quadrilateral and in in Reference 6 for a 4-node tetrahedron element. A general construction for the lumped mass matrix (LMM) for elements with Allman's rotation is given in Reference 16. It uses the row-sum lumping procedure for the displacement degrees of freedom and scaled values for Allman's rotation to reach the critical time step of the linear elements. Unfortunately, this LMM provides less accuracy of eigenfrequencies compared to CMM in Reference 17. Moreover, the scaling of the rotary inertia in LMM should be adjusted to the stabilization stiffness of the element. Therefore, CMM is preferred in implicit analysis. 18,19 Based on this profound research, finite elements with Allman's rotation made its way to academic finite element packages, for example, as 3-node triangle TrPlaneStrRot in the open-source code OOFEM 20 and commercial programs, for example, 4-node tetrahedron (ELFORM = 4) and 8-node hexahedron (ELFORM = 3) in LS-DYNA. The latter two formulations are used in explicit dynamics. In the context of the explicit dynamics, a further feature of the elements with Allman's rotations becomes important. These elements have a feasible time step comparable to linear elements and a twice larger feasible time step than a quadratic element. There is an even higher potential for these elements if they are combined with one of the mass scaling or inertia tailoring schemes, which can yield a larger feasible time step or an accuracy close to one obtained with the CMM.
The basic idea of mass scaling in the context of nonlinear structural dynamics is to add artificial terms to the mass matrix, thus reducing the highest eigenfrequencies of the discretized structure, while changing the lower eigenmodes and frequencies as little as possible. 21 Mass scaling aims at reduction of the computational cost since the highest eigenfrequency of the system limits the critical time step for explicit time marching. Adding artificial terms to only diagonal terms of the mass matrix is known as conventional mass scaling and it is practiced since the 1970s. Conventional mass scaling increases the translational mass of the structure. Therefore, it is very efficient if only a few stiff or short elements with high eigenfrequencies limit the feasible time step and the time step can be manipulated at the price of relatively small growth of the total mass. Alternatively, selective mass scaling adds both diagonal and off-diagonal terms to the mass matrix to preserve the translational mass of the structure. [21][22][23][24][25] Selective mass scaling comes at the cost of the solution of a linear system of equation for affected degrees of freedom each time step. This solution is done usually by the conjugate gradient method with the Jacobi preconditioner. 26 To avoid this step, several formulations for reciprocal mass matrices were introduced recently. Reciprocal mass matrices (RMM) are sparse matrices that connect the total force vector to the nodal acceleration vector by mere matrix-vector multiplication. They can be assembled from elemental matrices and their usage does not require large changes of the finite element code. Algebraic constructions of the reciprocal mass matrix are introduced in papers. [27][28][29][30] These formulations apply to solid finite elements with displacement degrees of freedom and plate finite elements with rotations. The starting point of the formulation is a diagonal mass matrix and a stiffness matrix. Variational formulations are independently elaborated in References 31-33. The starting point of the variational construction is parametrized Hamilton's principle with independent fields for displacement, velocity, and linear momentum. Choosing an independent discretization of the linear momentum with a local dual basis the enables derivation of the sparse reciprocal mass matrix. These RMMs are always augmented to increase the feasible time step, or decrease the dispersion error (increase accuracy of the lowest modes), or to reach both simultaneously. The RMMs after the augmentation are called variationally scaled reciprocal mass matrices (VSRMM). 31,32 Dual bases are not available for the standard elements with Allman's rotations because the transformation matrix has a rank deficit and supports for rotary degrees of freedom are undefined. Therefore, a pure variational approach would fail and it leaves an open question for construction of the reciprocal mass matrix in the case of Allman's rotation. However, an alternative construction should satisfy a list of requirements denoted with M* as follows M1 preservation of the translational inertia, M2 assemblable from the local matrices, M3 having a mask of the consistent mass matrix or even a tighter mask, M4 having a larger feasible time step or/and accuracy than for the diagonal mass matrix, M5 dependence on a minimal number of tuning parameters.
This builds the first part of the article.
Feasible time step estimates for RMM is not always trivial. Global estimators, for example, based on forward iterations or Lanczos method are too expensive. Elemental estimators may be not conservative in the case of RMM as shown by the paper. 34 Therefore, the common approach is using nodal time estimators. Nodal time estimators were initially proposed for diagonal mass matrices in the context of dynamic relaxation 35 and explicit dynamics. 36 These estimates have roots in Gershgorin's circle theorem and use the corresponding bound for the largest eigenvalue. The largest eigenvalue is bound by the maximum absolute row-sum or column-sum of the stiffness matrix divided by the nodal mass. For symmetric stiffness matrices, no difference between row-and column-sum is observed. Recently, estimators based on Gershgorin's circle theorem were proposed for structural undamped systems with RMM and displacement degrees only. 34 In this case, the largest eigenvalue is bound by a row-sum or a column-sum of absolute values of entries of the product of the reciprocal mass matrix and the stiffness matrix. The rowwise and columnwise versions are not identical in this case and perform similarly for several studied finite element types and meshes. 34 These estimators usually leave a 10%-30% gap to the exact critical time step. This raises the question of whether a sharper estimator may be constructed by using a different eigenvalue bound. The second question is whether the proposed estimators are still efficient for the element with displacements and Allman's rotation, that is, values with mixed physical dimension (m) versus (-). A recent paper compares Gershgorin's, Parker's, Brauer's, Ostrowski's and trace eigenvalue bounds for time step estimates in case of a diagonal mass matrix with spectral finite elements for explicit acoustics. 37 The paper indicates that the Ostrowski's bound is the sharpest bound for most of the cases. Further study discovers that in the presence of mixed degrees of freedom, like displacement and rotations in a Bernoulli beam element, the Gershgorin's and Ostrowski's bounds substantially underestimate the critical time step. 38 Therefore, an additional similarity transformation by Fan 39 should assist the bounds. Such an estimator is tested for a 3-node plate element and a stiffened panel structure with a regular mesh and an identical similarity transformation matrix at each node. 38 This opens a perspective of adjusting of this estimator to the case of elements with Allman's rotations and testing whether this estimator satisfies following requirements denoted with T* T1 to be conservative, T2 to be efficient for finite elements with displacement and Allman's rotations, and T3 to provide satisfactory results for irregular and distorted meshes.
This builds the second part of the article.
The article is structured as follows: Section 2 described the problem statement and a general method for the construction of RMM for finite elements with Allman's rotations. An estimator for the critical time step based on relevant eigenvalue bounds is presented in Section 3. Section 4 contains a set of numerical experiments supporting the developments. This includes a comparison with existing nodal time step estimates and with competing finite element formulations with nodal displacements only. A discussion of the obtained results is given in Section 5. Finally, the main results and possible directions for future work are outlined in Section 6.

Problem statement
Consider a solid finite element with Allman-type rotations. Each vertex node has three degrees of freedom [u x , u y , z ] in 2D and six degrees of freedom [u x , u y , u z , x , y , z ] in 3D. Let N R and N q denote the matrix with the Allman-type interpolation of displacements and the matrix with the shape functions of the base quadratic finite element, respectively. Then the Allman-type interpolation can be written as where  is the transformation matrix between the Allman-type interpolation and the base quadratic element. Examples of base and Allman-type elements are shown in Figure 1. The interpolation couples a single vertex rotation to a displacement field in two orthogonal directions. Details on the transformation matrix are given in Appendix B. Problem statement. Find a general formulation for a reciprocal mass matrix for this class of elements and an appropriate feasible time step estimator that satisfy requirements M1-M5 and T1-T3 and test these developments for triangle Tri3R, quadrilateral Quad4R, and tetrahedron Tet4R.

Derivation of reciprocal mass matrices
The standard variational multiparameter construction of the reciprocal mass matrices for elements with displacement degrees of freedom only from in Reference 33 relies on a three-field Hamilton's principle with independent fields for the displacement u, several velocities v (i) , and the linear momentum p. The first variation of the principle read where T • (u, p, v (i) ) is the modified kinetic energy, Π ext (u) and Π ext (u) are the potential energies for the internal and external forces. t 0 and t end denote the initial and the termination time for the transient problem. The modified kinetic energy is expressed as where Ω is the domain of the body, is the density, C 2i is a template (free) parameter in the sense of parametrized variational principles introduced by Felippa. 40 The customization parameters should be positive and their sum must be strictly smaller than one as shown in Reference 33. Discretization of the independent fields with independent interpolations leads to the variationally scaled reciprocal mass matrix C • (VSRMM) in a multiparameter template form The augmentation of the reciprocal mass matrix̃• is computed using the projection matrix W, the mass matrix for ith velocity field Y (i) and the reciprocal mass matrix C. The latter three matrices are assembled from the local matrices with Ω e denotes the domain of individual elements. The matrices Y (i) are block diagonal for incompatible velocity approximations v (i) . That property insures a sparse shape of the resulting VSRMM. A set of free customization parameters C 2i allows tuning of the VSRMM to desired properties. The matrix C • is neither an inverse of the consistent mass matrix nor a sparse approximation of such an inverse. It is a sparse and assemblable approximation of the inertia properties of the structure. Interpolation functions for the linear momentum are always biorthogonal to the displacement shape functions as proposed in Reference 33. Biorthogonal functions are constructed in three steps: compute the local metric matrix m e ; construct unweighted dual functionsf rom primary functions N e using the inverse of the local metric matrix; weighting of the dual functions with density, the local support t e, i and the inverse of the global support (T j ). These steps are concisely given below.
= N e m −1 e , The direct transfer of this construction to Allman-type elements contains two issues. First, the transformation matrix  has a rank deficit, which leads to the failure of the construction of the unweighted dual basis in Equation (8) due to the singularity of the local metric matrix. If the pseudo-inverse of the metric matrix is used in Equation (8) then the construction does not preserve translational inertia even for a single finite element (thus violating the condition M1). Second, Allman's rotations couple interpolation displacements in x-, y-and z-direction and it is unclear how to extend the scalar support associated with the DOF. Therefore, an alternative approach is sought. The Allman's rotations add hierarchically truncated quadratic displacements N onto a linear finite element N lin with vertex displacements only. Thus, the rows of the Allman-type interpolation may be split as .
Therefore, the vertex displacements can be assigned with the VSRMM according to the variational construction described above in Equations (5)- (9). The construction needs approximation matrices of velocity fields (i) . Here, the number of the velocity fields is limited by two. The matrix (1) associated with the free parameter C 21 is always the identity matrix. The matrix (2) combines all linear functions in each spatial direction, see for details. 32,33 Reuse of the construction for vertex displacements guarantees preservation of the translational inertia (requirement M1). Now, the inertia at the Allman's rotations can be defined. If the inertia of Allman's rotations is decoupled from translational inertia then it does not compromise the preservation of the translational inertia, which is already provided for the construction in Equations (5)- (9). Therefore, the coupling terms in the reciprocal matrix between displacements u j and rotations j are taken zero. Furthermore, we propose to use for the Allman's rotations a diagonal shape of the local reciprocal mass matrix. This choice reduces the number of nonzero entries of the reciprocal mass matrix (requirement M3) and avoids new tuning parameters (requirement M5). Now, these diagonal values should be computed from local information. For comparison, diagonal entries of the local lumped mass matrix read where N ,j is jth column of the Allman-type interpolation in Equation (10). The possibility to scale the diagonal term of the rotary inertia for this type of elements is discussed in Reference 16. The factor 5 4 is a compromise between time step size and accuracy in the lowest modes for coarse and moderate meshes, see details in Appendix A. We can assemble the reciprocal values of the local lumped mass matrix if they are scaled with the square of the ratio of the local support to the global support. 33 This leads to the sought diagonal values The physical dimension of the entries c ,jj is reciprocal to the moment of inertia 1/(kg ⋅ m 2 ). Equation (12) completes the proposed construction of the reciprocal mass matrix.

FEASIBLE TIME STEP ESTIMATOR BASED ON OSTROWSKI'S AND FAN'S EIGENVALUE BOUNDS
For the central difference method and undamped problems, the feasible time step Δt is limited by the critical time step Δt crit through the stability criterion where max is the maximum eigenfrequency of the assembled system. 41 The maximum eigenfrequency of the global system with reciprocal mass matrices is determined by the standard eigenvalue problem where C • and K are the variationally scaled reciprocal mass matrix and the stiffness matrix and i , i , and i are the ith eigenvalue, eigenfrequency, and eigenvector, respectively. The VSRMM is a sparse matrix and the product C • K defines a nonsymmetric sparse matrix, whose largest eigenvalue limits the time step. Below, relevant eigenvalue bounds are considered and a novel time step estimator is proposed. This estimator uses the local information about mesh size for the similarity transformation and Ostrowski's bounds with exponents on a regular grid.

Considered eigenvalue bounds
j≠i |A ji | be the sums of absolute values of the nondiagonal entries in ith row and in ith column of a square real matrix A ∈ R n×n , respectively. Let denote with S i (A ii , P i (A)) ⊆ C a closed disk on complex plane centered at A ii with radius P i (A) ∈ R + . Gershgorin's circle theorem states that each eigenvalue of a square real matrix A lies within at least one of Gershgorin's disks S i (A ii , P i (A)). A corollary of this theorem leads to rowwise Gershgorin's bound of the largest eigenvalue The spectrum of a matrix A and its transpose A T coincide. Therefore, columnwise Gershgorin's bound of the highest eigenvalue is also valid. In the case of nonsymmetric matrices, the latter two bounds are not identical. The Gershgorin's bound can be extended in two ways. First, a bound with a weighted row sum is due to a note by Fan. 39 Spectrum of a matrix is invariant under a similarity transformation A → S −1 AS. Fan's bound is obtained from the rowwise Gershgorin's bound assuming a diagonal shape of the similarity transformation matrix with positive entries only This transformation is of special interest in case the system has degrees of freedom with different physical dimensions, for example, solid elements with Allman's rotations or thin-walled elements with rotations. A possible construction of the similarity matrix is discussed in Subsection 3.2. The second extension is due to Ostrowski's circle theorem. It states that every eigenvalue of a square real matrix A lies within at least one of the disks or the best from all admissible values of O max ≤ min The rowwise and columnwise Gershgorin's disk systems are special cases of Ostrowski's disk systems for = 0 and = 1, respectively and the Ostrowski's bound with the best choice of cannot be worse than any of two Gershgorin's bounds. In the case of a symmetric matrix, the values ofP i (A) and P i (A) coincide and the Ostrowski's bound reduces to the Gershgorin's bound. This opens an opportunity for constructing sharper time step estimators using Ostrowski's bound for nonsymmetric matrices, like the product of the reciprocal mass matrix and the stiffness matrix in Equation (14). Computing minimax problem in Equation (19) for continuous the exponent is challenging. Herein, a discrete uniform grid for the exponent is used. Similar estimator was introduced on the elemental level in Reference 37, but it is not clear whether continuous or discrete values of the exponent are used.

Similarity matrix for Fan's bound
The similarity matrix in Fan's bound improves dramatically the tightness of the time step estimator. The following straightforward example illustrates this statement. Consider a single Euler-Bernoulli beam in 2D with cubic Hermite shape functions and local degrees of freedom . The element has a constant rectangular cross-section b × t and length L. The material of the beam is given by density and Young's modulus E. The element matrices and a corresponding matrix in a nonsymmetric eigenvalue problem read The exact maximum eigenvalue of the matrix A can be analytically determined and it is equal to The columnwise Gershgorin's estimate given in Equation (16) for the matrix A in Equation (20) yields Analysis of the physical units in Equation (22) reveals a contradiction. Both arguments of the max function are the sum of values with units m 2 /s 2 and m/s 2 . This contradiction destroys the sharpness of the estimate for small absolute numbers for element length L → 0, that is, the ratio of the estimate to the exact value grows to infinity. These small absolute numbers are also achievable by changing length units used in the model, for example, from mm to m. The rowwise Gershgorin estimate behaves similarly. Fan's trick with a properly chosen similarity matrix corrects the mismatch of the physical units. Let us select a diagonal similarity transformation matrix S where entries possess physical dimensions of corresponding degrees of freedom, that is, displacement-meter (m) and rotation-radian (-). Fan's bound for the single beam element and exemplary transformation matrix read with diag[×] being the constructor of a diagonal matrix from a vector and being a positive dimensionless scalar. The result of this estimate has the correct physical unit of m 2 /s 2 . This motivates us to use a similarity transformation with some characteristic edge length at each displacement degree of freedom. An algorithm for the computation of nodal characteristic length is specified below. Nodal characteristic length is computed as minimum edge length for all adjacent elements where element index e loops over all adjacent elements  j to node j and L e is the smallest edge of the element e. The similarity matrix is defined as where  d and  are nodal indexes for displacement and Allman's degrees of freedom, respectively. This matrix is computed once. The selection of suitable values for is an additional issue. Different values of dimensionless parameter result in the following estimates At the optimal value = 2, the estimate coincides with the largest eigenvalue. Such serendipity is not observed for real structures with more than one element. Nevertheless, values of close to two leave a small gap to max for most of the regular meshes. Another feature is that the estimate stays accurate for > 2 for a thin beam element (L ≫ t) and for < 2 for thick beam element (t ≫ L). Therefore, values of far from two may be useful for distorted and irregular meshes. Herein, a constant value for a complete structure is used. Algorithm 1. Time step estimator using Ostrowski's bound Input: K, C • , S, n g , Output: Δt crit -critical time step estimate for an undamped system Require: S -diag, n g ≥ 2 -grid size for Ostrowski exponents, n g ∈ [11..21] is a good choice 1: procedure OstrowskiEstimateWithFanTrick(K, C • , S, n g ,) 2: for j ← 1 to size(K) do 5: for k ← 1 to n g do 8: ⊳ Uniform grid for sampling with 9:

Novel time step estimator
We propose to estimate the largest eigenvalue using the Ostrowski's bound given in Equation (18) and combine it with the diagonal similarity transformation in given Equation (25) as suggested by Fan 39 where exponents j belong to the regular grid of n g points in the interval [0, 1] This estimator exploits unsymmetry of the matrix product S −1 C • KS. An algorithm is given in Box 1. The performance of the estimator is controlled by two parameters: the number of sampling points for exponent n g and the dimensionless parameter in the similarity matrix .
The exact computational cost of the proposed time step estimator (30) depends on the implementation details. An estimation of the computational cost is possible analogously to in Reference 34 using a simple floating-point operation (FLOP) count of steps in the Algorithm 1. Following simplifying assumptions are made here: The dimension of all matrices is denoted as n and all Steps with a cost much smaller than n are skipped. The matrices C • , K, and C • K are assumed to have an average bandwidth b, b, and a, respectively. The grid size for Ostrowski exponents and the average bandwidth are assumed to be much smaller than dimension n. Two calls of function power in Step 9 are assumed to cost 20 operations each. The cost of the computation of matrices in the input (Step 1) is excluded from the estimate.
FLOP counts for the sparse matrix-matrix product C • K (Step 2) is the main cost with 2nb 2 + 2na. Here, an estimate according to in Reference 42 of the Gustavson algorithm 43 is used. Steps 5-12 are repeated n times.
Step 5 and 6 need a − 1 FLOPs each. Steps 8-10 are repeated n g n times.
Step 8 and 10 need one FLOP.
Step 9 needs 41 FLOPs (2 calls pow() and one mul()). This leads to the estimate This estimate can be expanded for some specific mesh types, element and the grid size for Ostrowski exponents. The Quad4R element with a regular mesh has the average bandwidth a = 75 and b = 27. For the grid size n g = 21 the overall cost COST Quad4R,reg,n g =21 = 1891n FLOP.
The Tri3R element with a regular unidirectional mesh has the average bandwidth a = 51 and b = 21. For the grid size n g = 21 the overall cost COST Tri3R,uni-dir,n g =21 = 1507n FLOP.
Repeating this algorithm each time step is extremely inefficient. Thus, it should be executed at the beginning and at a few intermediate time points of the simulation.

NUMERICAL EXAMPLES
In this section, the proposed construction of VSRMM and the novel time step estimate are tested for several finite elements with Allman's rotations. First, developments are assessed for triangular and quadrilateral elements for eigenfrequency benchmark NAFEMS FV32. 44 Second, the performance of the time step estimate is checked for a highly distorted mesh in 2D. Third, a transient example based on the NAFEMS FV32 is computed using explicit time integration. Finally, an example for the 4-node tetrahedron is considered. All computations are performed with double precision. Stiffness stabilization with 1 = 10 −6 and stabilization of the consistent mass matrix 2 = 10 −5 are used for all cases as described in Appendix C. Elemental time step estimate for LMM uses the exact highest local eigenfrequency. It is a forced solution because there is no reliable reference for a characteristic length of the element used in Courant-Friedrich-Levy estimate. The latter is used in most of the commercial codes and it is not so tight as the exact local eigenfrequency estimate. Quadrature rules for each element type are specified in examples. Mainly, the compatible strain field formulation is used for the computation of the stiffness matrix.

NAFEMS FV32 for Tri3R
The eigenfrequency benchmark FV32 from the NAFEMS collection of standard benchmarks 44 considers a two-dimensional free vibration problem of a tapered membrane with length L = 10 m, width at wide end W = 5 m and width at the narrow end w = 1 m. The membrane is clamped at the wide end (u x = u y = z = 0). The standard provides the six lowest eigenfrequencies with an accuracy of five digits. This solution is obtained using the Rayleigh-Ritz method with global ansatz functions. The membrane is meshed with Tri3R elements under plane stress assumption. Gauss-Radau quadrature rule with seven points is used for the computation of the stiffness matrix and the reciprocal mass matrix. The stiffness matrix is computed via standard displacement-based formulation. The VSRMM is computed with one customization parameter. The homogeneous isotropic material of the membrane has the elasticity modulus E = 200 GPa, Poisson's ratio = 0.3, and density = 8000 kg/m 3 . Two unidirectional mapped meshes (n x × n y ) with 10 × 5 and 20 × 10 elements are considered. The VSRMM presented here is used with three different values of parameter C 21 . Ostrowski's estimate is computed with the number of sampling points n g = 21.
Results for the six lowest eigenfrequencies and the highest eigenfrequency are presented in the upper part Table 1. A starting guess of the free parameter is C 21 = 0.66 ≈ 21 32 . This value is taken from dispersion analysis for the constant strain triangle element Tri3 at infinite hexagonal mesh. It provides the best low frequency accuracy for the Tri3 element. As this element is the basis for the translational part of the VSRMM construction, a good transfer of the results is expected for the Tri3R element. VSRMM with C 21 = 0.66 provides spectrum even more accurate than for the consistent mass matrix with the highest frequency by 5% less than for LMM. Increasing the free parameter to 0.81 and 0.95 reduces the accuracy of the eigenfrequencies and increases the potential speed-up to 11% and 16%, respectively. Figure 2 shows six physical eigenmodes. The eigenmodes capture the physical behavior correctly and correspond to the reference modes. Both spectra of LMM and VSRMM possess spurious modes in the range of computed frequencies. These modes are shown in Figure 3. A pattern similar to fish scales is characteristic for the spurious modes in the unidirectional mesh considered here. These TA B L E 1 Six lowest eigenfrequencies and the highest eigenfrequency for the FV32 benchmark computed with CMM, LMM, and VSRMM spurious modes are represented mainly by vertex rotations and their frequency is almost independent of the free parameter C 21 . Increasing the stiffness stabilization parameter 1 to 5 ⋅ 10 −4 shifts the spurious frequencies from the range of the considered frequencies to a higher frequency range. The bottom part of Table 1 presents the six lowest eigenmodes obtained with a competing finite element formulation with nodal displacements only: the constant strain triangle Tri3. The Tri3 element uses the row-sum lumping procedure. The 6-node quadrilateral element Tri6 is omitted from the consideration because of a poor quality of the lumped mass matrix (row-sum lumping yields zero inertia at corner nodes and Hinton-Rock-Zienkiewicz (HRZ) procedure 45 lead to high eigenfrequency error). The error in the lowest eigenmodes for the Tri3R element with VSRMM with C 21 = 0.66 is much smaller than for the competing formulation even if the refined mesh is used. The  error in the lowest eigenmodes for the Tri3R element with VSRMM with C 21 = 0.81 is similar to one obtained with the refined mesh. Furthermore, the highest frequency for the Tri3R element with VSRMM with C 21 = 0.81 is only a half of the value for the refined mesh with the Tri3. This illustrates the advantage of the proposed formulation at coarse meshes.
The convergence of the six lowest eigenfrequencies for the FV32 benchmark for the standard LMM and the proposed VSRMM with C 21 = 0.66 and C 21 = 0.95 is presented in Figure 4. The spurious eigenfrequency is excluded from the ordering. The reference values for the frequencies are taken from the NAFENS standard. The default NAG eigenvalue solver with double precision is used for the computation of the generalized eigenvalue problem within a Maple implementation. The number of elements along the x-axis of the model and the absolute relative error in the frequencies are used for the abscissa and for the ordinate, respectively. The eigenfrequencies show an asymptotic convergence rate close to quadratic. The error constant at the two lowest eigenfrequencies for the LMM is twice larger than for VSRMM with C 21 = 0.66 but it is only 1 6 of the value for C 21 = 0.95. Thus, the trade-off between the error constant (increase by factor 12) and maximum eigenfrequency (reduction only by 9%) is relatively bad for the proposed construction of the VSRMM.
The performance of the proposed time step estimator for the FV32 benchmark is illustrated in Table 2. The estimator is inefficient if the similarity transformation is omitted especially in the case where the geometry of the model is scaled by factor 10 −3 . Using the best of rowwise or columnwise Gershgorin's estimators with the similarity transformation leaves a gap of 22.0% to the exact feasible time step. This gap is reduced to 14.2% for the Ostrowski's estimator for the best choice TA B L E 3 Six lowest eigenfrequencies and the highest eigenfrequency for the FV32 benchmark computed with CMM, LMM, of the dimensionless parameter = 2.0. This is usually an optimal value for for regular meshes. Moreover, the obtained feasible time step is 25.5% larger than the feasible time step for LMM obtained from elemental estimators using exact local eigenfrequencies. Finally, different exponents yield the best estimates. Using a grid with 21 sampling points seems to be sufficient.

NAFEMS FV32 for Quad4R
Now the problem from the previous subsection is solved with the Quad4R elements. 3 × 3 Gauss quadrature rule is used for the computations of the stiffness matrix and the variationally scaled reciprocal mass matrix. The stiffness matrix is computed via the standard displacement-based formulation. The reciprocal mass matrix is computed with up to two customization parameters. The results for the six lowest eigenfrequencies and the largest eigenfrequency are presented in the upper part of Table 3. First, parameters of the VSRMM are selected to obtain a better accuracy in the lowest eigenmodes. A starting guess is an optimized value of C 21 = 0.55 ≈ 5 9 shown to give superconvergent results for longitudinal vibration problems for a 2-node rod in Reference 31. The Quad4 element uses tensor product shape functions of the 2-node rod. The VSRMM with C 21 = 0.55 and three different values for C 22 provides a spectrum even more accurate than for the CMM. Extremely good results are observed for the longitudinal modes of the tapered plate (3 and 6). Selecting C 22 = 0.4 reduces the largest eigenfrequency to about 8% below the one for LMM. VSRMM with larger values for C 21 , like 0.69 or 0.83, decreases the accuracy of the eigenfrequencies substantially for only moderate speed-up up to 16% with respect to LMM. The eigenmodes capture the physical behavior correctly and correspond to the reference modes. Both spectra of LMM and VSRMM possess spurious modes with a chevron pattern, shown in Figure 5. The frequency of spurious eigenmodes with VSRMM is almost independent of the free parameters C 21 and C 22 . The bottom part of the Table 3 presents the six lowest eigenmodes obtained with three competing finite element formulations with nodal displacements only: 8-node serendipity quadrilateral element Quad8 with HRZ lumping, fully integrated and enhanced assumed strain 4-node quadrilateral elements Quad4 and Quad4E4, 46 respectively. Both 4-node elements use the row-sum lumping procedure. The error in the lowest eigenmodes for the Quad4R with C 21 = 0.55 and C 22 = 0.4 is twice smaller than for the 4-node elements with the same mesh and comparable for the double refined mesh. The highest eigenfrequency for the double refined mesh is by 42% larger than for the Quad4R with C 21 = 0.55 and C 22 = 0.4. The 8-node element with LMM yield an extremely high The convergence of the six lowest eigenfrequencies excluding the spurious eigenfrequency for the FV32 benchmark for the standard LMM and the proposed VSRMM with C 21 = 0.55 and C 22 = 0.4 is shown in Figure 6. The plot setup and the reference values are identical to the convergence plot for the Tri3R element. The eigenfrequencies show an asymptotic convergence rate close to quadratic. The error constant for LMM is twice larger than for VSRMM with C 21 = 0.55 and The performance of the proposed time step estimator for the Quad4R finite element is illustrated in Table 4. The proposed estimator provides the gap to the exact feasible time step of 21.5%, which outperforms the best of the rowwise or the columnwise Gershgorin's estimators with a gap of 27.6%. The best choice of the dimensionless parameter is close to 2.0. Moreover, the obtained estimator can provide a speed-up of about 8.0% with respect to the feasible time step for LMM obtained from elemental estimators.

Highly distorted mesh in 2D
In this example, the tightness of the estimators is tested for the highly distorted mesh proposed in Reference 34. The problem setup for two element types is shown in Figure 7. Several values of the dimensionless parameter for the similarity transformation are tested. One and two customization parameters are used for the reciprocal mass matrix of Tri3R and Quad4R, respectively. Ostrowski's estimator is computed with the number of sampling points n g = 21.
The performance of the proposed time step estimators is presented in Table 5. The proposed estimator outperforms the best of the rowwise or the columnwise Gershgorin's estimators. However, the tightness of the estimator is much more sensitive to the choice of the dimensionless parameter for similarity matrix for meshes with Tri3R and Quad4R. The mesh with Quad4R requires values below 0.3 to get any speed-up compared to the model with LMM. This illustrates that an alternative selection of (uniform in mesh or nonuniform) is needed for nonregular meshes. Dimensions: l x × l y × l z = 1 × 1 × 10

Transient problem for a tapered membrane with Tri3R
Consider a transient problem based on the FV32 NAFEMS benchmark modeled with Tri3R elements and a mapped mesh 10 × 5 elements. The motion due to an abrupt load of F = 10, 000 N in y-direction applied at the right upper corner node is studied at the first 0.1 s, which corresponds to 4.5 periods of the first frequency of the structure. Zero initial displacements and velocities are assumed. All solutions are obtained with the central difference method with a time step being 90% of the feasible value. The feasible time step is obtained by the proposed approach with n g = 21 and = 2.0 and elemental estimates for VSRMM and LMM, respectively. The reference solution is obtained with an overkill mesh 64 × 32 elements with the Tri3 using the lumped mass matrix. The time history of the displacement in y−direction at the location of the applied force is compared in Figure 8. The VSRMM with C 21 = 0.66 and C 21 = 0.95 needs 3% and 21% less time steps than the LMM, respectively. These numbers of time steps are in good agreement with the ratios of the highest eigenfrequencies (global estimates) given in Table 1. The VSRMM with C 21 = 0.66 yields results very close to the reference solution. The results with LMM are less accurate. Least accurate results are obtained for the VSRMM with C 21 = 0.95. These result are in good agreement with the accuracy of the first mode presented in Table 1.

Sze's beam with Tet4R
The next example considers a beam with square cross-section discretized with 12 tetrahedral elements Tet4R. The problem setup is shown in Figure 9. Quadrature rule with 15 points by Felippa 47 is used for the computation of the stiffness matrix and the reciprocal mass matrix. The stiffness matrix is computed via hybrid-stress formulation with 18 stress modes by Sze. 6 Ostrowski's estimator is computed with the number of sampling points n g = 21.
The results for the tightness of the estimators are given in Table 6. The performance of the estimator is less sensitive to than in the previous examples. Values in the range [2..6] are acceptable. The obtained estimate can provide a speed-up of 11.5% with respect to the feasible time step for LMM obtained from elemental estimators. Similar performance for the compatible-strain formulation of the stiffness matrix are obtained.

DISCUSSION
The proposed construction of variationally scaled reciprocal mass matrices for elements with Allman's rotation meets requirements (M1-M5) stated in the introduction. First, the formulation preserves translational inertia and it is assemblable from the local matrices. Second, the number of nonzero entries of the proposed VSRMM is less than for CMM due to diagonal inertia for the rotary part and decoupled terms between translational and rotary parts. Third, the feasible Potential speed-up with respect to LMM is between 8% and 23%, which is considered as a moderate success. Partially, it is explained by the used elemental time step estimate for LMM that computes the exact highest local eigenfrequency and provides tighter result that the standard one, for example, based on Courant-Friedrich-Levy criterion. Courant-Friedrich-Levy estimate was not possible to apply due to missing reliable reference for a characteristic length of the element.
Second-order convergence rate is obtained for lowest eigenfrequencies for LMM and proposed VSRMM. This is in agreement with the fact that the elements with Allman's rotations use truncated quadratic polynomials. The special instances of VSRMM mentioned above improve the error constant by a factor two w.r.t. one for LMM. An additional feature of the elements with Allman's rotations is presence of spurious modes, which is shown in various benchmarks with LMM and VSRMM. The spurious modes possess fish scale and chevron patterns of deformation for triangular and quadrilateral elements, respectively. They can be neglected because it is difficult to activate them with the standard load cases.
The proposed estimator for a feasible time step meets requirements (T1-T3). First, the estimator is conservative because it has a solid basis in Ostrowski's circle theorem. Second, the estimator uses absolute sum in row and column of the special matrix constructed from VSRMM and the stiffness matrix. Therefore, it is a local estimator using nodal information only. Moreover, Ostrowski's circle theorem includes Gershgorin's circle theorem as a special case, thus the novel estimator always improves the best of rowwise and columnwise Gershgorin's estimates. In most of the considered examples, this improvement is by 5%-15% relative to the exact feasible time step. To reach this improvement, usage of 11 to 21 grid points for exponents in the Ostrowski's estimator is proposed. Third, the tightness for elements with Allman's rotations is additionally increased by a diagonal similarity transformation that takes into account mismatch of physical units for nodal displacements and Allman's rotations. This similarity transformation uses only one dimensionless parameter and local geometrical information from the mesh. Optimal values of for regular meshes is about two. However, an example with highly distorted mesh indicates that low values of should be used to get desirable tightness or even nonuniform values for different nodes may be applied. This issue should be investigated further.
The starting point for the proposed estimator for the feasible time step is an unsymmetric eigenvalue problem (14). Alternatively, a symmetric reduction of the problem can be exploited as suggested in Reference 34. This reduction requires Cholesky decomposition of the reciprocal mass matrix with C • = L T L and it leads to a symmetric eigenvalue problem with the same spectrum Each entry of the matrix has now the correct physical unit of m 2 /s 2 . Therefore, Gershgorin's estimate can be directly applied and similarity transformation from Fan's trick may be avoided here. The estimate is not considered for two following reasons. First, preliminary tests showed that it leaves a larger gap to the exact eigenvalue in comparison with the proposed estimator. Second, the additional overhead for computing and storing of Cholesky decomposition is much higher than for storing the matrix products in the proposed estimator.

CONCLUSIONS
In this article, a construction of variationally scaled reciprocal mass matrices for elements with Allman's rotations and a novel estimator for feasible time step have been presented. The proposed construction uses a variational approach for displacement DOF's and an algebraic construction with diagonal values for Allman's rotation. This construction is assemblable from local matrices and can be included in a standard finite element code. Eigenvalue benchmarks prove good accuracy of the lowest eigenfrequencies and a moderate speed-up with respect to LMM. Both LMM and VSRMM exhibit spurious eigenmodes in the lower range of the spectrum for small values of the stiffness stabilization parameter. The proposed time step estimate is based on Ostrowski's circle theorem and a diagonal similarity transformation. This combination proves to be extremely useful for elements with Allman's rotations. The conservativeness of the proposed estimator is confirmed by all examples with regular as well as highly distorted meshes. Possible future steps include extensions of the feasible time step estimators to systems with penalty contact and damping, like Rayleigh damping. The standard central difference scheme reduces the stability limit for activated penalty contact pairs and increasing critical damping ratio. 48 Furthermore, the similarity transformation for distorted meshes may get an alternative construction. and differ in selected value of the scaling parameter . Uniform values in the complete mesh and values adjusted at each element are possible. In the context of dispersion minimization, a uniform value of = 1 is used for the scaling parameter. 49 Unfortunately, this value restricts the critical time step significantly. An opposite approach is proposed in Reference 16, where the value is increased until the critical time step of the element with Allman's rotations is equal to critical time-step of the equal-sized low-order element. These value destroy accuracy of the lowest eigenmodes as reported in Reference 17. A compromise value of 5 4 is used here uniformly. It improves the critical time step by approximately 10% w.r.t. = 1 and preserves the accuracy in the lowest modes. Table A1 illustrates the influence of the scaling parameter on the highest frequency for single elements with different shapes. Two main observation can be deduced from it. First, bad aspect ratios lead to relatively high eigenfrequencies w.r.t. equal-sized low order element. Second, reduction of the eigenfrequencies from 45% to 71% is achieved with = 5 4 for all shapes w.r.t. the base quadratic element.

APPENDIX B. TRANSFORMATION MATRICES FOR TRI3R AND QUAD4R ELEMENTS
Allman-type interpolation is defined by a transformation matrix from the shape function of the base quadratic element, see Equation (1). 4-node tetrahedron Tet4R uses the transformation matrix given by Sze in Reference 6. It is not reproduced here for space reasons. Triangular element Tri3R uses the transformation matrix given by Allman where the differences between nodal coordinates are defined with x ij = x i − x j and y ij = y i − y j . Ordering of the local degrees of freedom follows the ordering of nodes giving first all displacements and continued by Allman's rotations, for example, [u x1 , u y1 , u x2 , · · · , u yx , z1 , z2 , z3 ]. The quadrilateral element Quad4R relies on the transformation matrix given by Sze in Reference 10 with (B2)

APPENDIX C. SPURIOUS ZERO-ENERGY MODES AND ASSOCIATED STABILIZATION MATRICES
Spurious zero-energy modes were identified in early developments for the finite element with Allman's rotations, for example, one for triangular and quadrilateral element. 2,10 These modes reed