Tutorial: lumped mass matrices of the shells with an offset#
Date: 24 of October 2026
Author: Saullo G. P. Castro
The lumped mass matrices of the shell elements Quad4, Quad4R, Tria3R and Tria3DSG (mtype=1 of update_M) are derived here from the kinetic energy of the shell, as a nodal quadrature, and verified with the rigid-body inertia of one element and of a mesh. When the reference surface of the nodes is at an offset from the mid-surface, the first moment of the density through the thickness, intrhoz, couples the translations with the rotations. The derivation shows that this
coupling stays in the lumped matrix, in the block of each node, and that it is what places the centre of mass of the element at the offset.
References:
Brockman, R. A., Dynamics of the bilinear Mindlin plate element, International Journal for Numerical Methods in Engineering, 24(12), 2343-2356, 1987, doi:10.1002/nme.1620241208: the kinetic energy, Eq. (20), the consistent mass matrix, Eq. (21), and the lumped mass matrix by Lobatto integration at the nodes, pp. 2348-2349;
Cowper, G. R., Gaussian quadrature formulas for triangles, International Journal for Numerical Methods in Engineering, 7(3), 405-408, 1973, doi:10.1002/nme.1620070316: the three-point rule of the consistent matrix of
Tria3R, first row of Table I.
The helper functions are those of tests/test_element_mass.py, imported below.
[1]:
import os
import sys
sys.path.insert(0, os.path.join('..', '..', 'tests'))
import numpy as np
import sympy as sp
from scipy.sparse import coo_matrix
from pyfe3d.shellprop_utils import isotropic_plate
from pyfe3d import DOF, INT, DOUBLE
from test_element_mass import (SHELLS, TRIANGLES, QUADS, area, area_rule,
shell_element, rotation_blocks, dense_M,
rigid_modes, shell_rigid_inertia, quad_N_detJ,
mass_moments)
Kinetic energy of the shell#
With the displacements \(u + z r_y\), \(v - z r_x\) and \(w\) at the distance \(z\) from the reference surface, Eq. (1) of Brockman (1987) and the kinematics of the shells of pyfe3d, the kinetic energy per unit area, integrated through the thickness \(z \in [-h/2 + e, h/2 + e]\) of a reference surface at the offset \(e\) from the mid-surface, is
Eq. (20) of Brockman (1987), with \(\boldsymbol{d} = (u, v, w, r_x, r_y, r_z)\) and \(I_\rho\), \(I_{\rho z}\), \(I_{\rho z^2}\) the intrho, intrhoz and intrhoz2 of ShellProp. The drilling rotation \(r_z\) has no inertia. SymPy gives the three integrals and the matrix \(\boldsymbol{B}\), the inertia of the section:
[2]:
z, h, e, rho = sp.symbols('z h e rho', real=True)
ud, vd, wd, rxd, ryd, rzd = sp.symbols('udot vdot wdot rxdot rydot rzdot')
lim = (z, -h/2 + e, h/2 + e)
I0, I1, I2 = [sp.factor(sp.integrate(rho*z**k, lim)) for k in range(3)]
print('intrho =', I0)
print('intrhoz =', I1)
print('intrhoz2 =', I2)
T = sp.integrate(rho/2*((ud + z*ryd)**2 + (vd - z*rxd)**2 + wd**2), lim)
d = sp.Matrix([ud, vd, wd, rxd, ryd, rzd])
B = sp.hessian(T, d).applyfunc(sp.factor)
Is, Isz, Isz2 = sp.symbols('I_rho I_rhoz I_rhoz2')
B = B.subs({I0: Is, I1: Isz, I2: Isz2})
B
intrho = h*rho
intrhoz = e*h*rho
intrhoz2 = h*rho*(12*e**2 + h**2)/12
[2]:
The offset couples \(u\) with \(r_y\) and \(v\) with \(r_x\) through \(I_{\rho z} = \rho h e\), which is zero only for \(e = 0\) (or for a laminate whose density is symmetric about the reference surface).
The consistent mass matrix#
With the nodal velocities \(\dot{\boldsymbol{q}}\) and \(\dot{\boldsymbol{d}} = \boldsymbol{S} \dot{\boldsymbol{q}}\), \(\boldsymbol{S} = [S_1 \boldsymbol{I}_6, S_2 \boldsymbol{I}_6, \ldots]\), the kinetic energy of the element is \(T = \frac{1}{2} \dot{\boldsymbol{q}}^\top \boldsymbol{M} \dot{\boldsymbol{q}}\) with
the form of Eq. (21) of Brockman (1987): the 6 x 6 block between the nodes \(i\) and \(j\) is the inertia of the section times \(h_{ij}\). For the linear triangle in area coordinates, \(h_{ij} = A/6\) for \(i = j\) and \(A/12\) otherwise, which the three-point rule of Cowper (1973) integrates exactly:
[3]:
N1, N2, A = sp.symbols('N1 N2 A', positive=True)
Ntri = sp.Matrix([N1, N2, 1 - N1 - N2])
# NOTE the reference triangle has the area 1/2, the Jacobian determinant is 2A
H_exact = 2*A*(Ntri*Ntri.T).applyfunc(
lambda f: sp.integrate(f, (N2, 0, 1 - N1), (N1, 0, 1)))
H_exact
[3]:
Lumped mass matrix by nodal quadrature#
A quadrature with the points \(\boldsymbol{x}_k\) at the nodes and the weights \(w_k\) integrates the whole integrand \(\boldsymbol{S}^\top \boldsymbol{B} \boldsymbol{S}\) as
Since \(S_i(\boldsymbol{x}_k) = \delta_{ik}\), each point contributes only to the block of its own node,
The matrix is block-diagonal per node and each block is the complete inertia of the section, coupling included: lumping acts on the shape functions, never on \(\boldsymbol{B}\). The matrix is diagonal only when \(I_{\rho z} = 0\), which is the condition stated by Brockman (1987), p. 2349, for his lumped matrix of the quadrilateral by Lobatto integration at the nodes. Dropping the coupling with an offset is a further approximation, not a consequence of the lumping.
The rules of pyfe3d are the vertex rule for the triangles, \(w_i = A/3\), and the two-point Gauss-Lobatto rule in each direction for the quadrilaterals, \(w_i = \det J(\xi_i, \eta_i)\), the Jacobian determinant at node \(i\) times the weights \(1 \times 1\). For the triangle, the vertex rule applied symbolically to the integrand gives:
[4]:
Bs = sp.Matrix(B)
def integrand(Nvals):
Nmat = sp.Matrix.hstack(*[Ni*sp.eye(6) for Ni in Nvals])
return Nmat.T*Bs*Nmat
vertices = [(1, 0, 0), (0, 1, 0), (0, 0, 1)]
M_lump = sum((A/3*integrand(p) for p in vertices), sp.zeros(18, 18))
print('blocks between different nodes are zero:',
all(M_lump[6*i:6*i+6, 6*j:6*j+6] == sp.zeros(6, 6)
for i in range(3) for j in range(3) if i != j))
M_lump[0:6, 0:6]
blocks between different nodes are zero: True
[4]:
Row-sum lumping gives the same blocks with \(w_i = \sum_j h_{ij} = \int_A S_i dA\), since \(\sum_j \boldsymbol{M}_{ij} = (\sum_j h_{ij}) \boldsymbol{B}\) and \(\sum_j S_j = 1\). For the triangle, \(\int_A S_i dA = A/3\) and the two lumpings coincide:
[5]:
M_cons = sp.zeros(18, 18)
for i in range(3):
for j in range(3):
M_cons[6*i:6*i+6, 6*j:6*j+6] = H_exact[i, j]*Bs
row_sum = sp.zeros(18, 18)
for i in range(3):
row_sum[6*i:6*i+6, 6*i:6*i+6] = sum((M_cons[6*i:6*i+6, 6*j:6*j+6]
for j in range(3)), sp.zeros(6, 6))
print('row-sum lumping equal to the vertex rule:',
sp.simplify(row_sum - M_lump) == sp.zeros(18, 18))
row-sum lumping equal to the vertex rule: True
Rigid-body inertia#
The inertia of a rigid body, with the mass \(m\), the first moment of the mass \(\boldsymbol{Q} = \int \rho \boldsymbol{p} dV\) and the second moments \(\boldsymbol{I}_2 = \int \rho \boldsymbol{p} \boldsymbol{p}^\top dV\), is the 6 x 6 matrix \(\boldsymbol{\Phi}^\top \boldsymbol{M} \boldsymbol{\Phi}\), where the columns of \(\boldsymbol{\Phi}\) are the nodal values of the three translations and of the three rotations about the origin. For a point \(\boldsymbol{p} = \boldsymbol{x} + z \boldsymbol{n}\) of the shell, with \(\boldsymbol{x}\) on the reference surface and \(\boldsymbol{n}\) its normal, the exact values are integrals over the area, and the nodal quadrature replaces them with sums:
Hence:
the mass is exact when \(\sum_i w_i = A\), true for both rules;
the first moment through the thickness, \(I_{\rho z} A \boldsymbol{n}\), is exact under the same condition, and only with the coupling: without it, the centre of mass of the element is on the reference surface instead of at the offset;
the first moment in the plane, \(I_\rho \sum_i w_i \boldsymbol{x}_i\), is exact when \(w_i = \int_A S_i dA\), because \(\int_A \boldsymbol{x} dA = \sum_i \boldsymbol{x}_i \int_A S_i dA\) for an isoparametric element. This holds for the triangle and for the parallelogram, where \(\det J = A/4\) is constant, but not for a general quadrilateral with the Gauss-Lobatto weights;
of the second moments, the term \(I_{\rho z^2} A\) is exact, while \(I_\rho \sum_i w_i \boldsymbol{x}_i \boldsymbol{x}_i^\top\) concentrates the mass at the nodes. About the centroid of the element it is 4 times the exact value for a triangle and 3 times for a parallelogram: the inertia of a single element is overestimated, an error that decreases with the size of the elements, see the last section.
The weights of the distorted quadrilateral of the tests, with the Gauss-Lobatto rule and with row-sum lumping:
[6]:
xy = QUADS['distorted']
Aq = area(xy)
nodes = [(-1, -1), (1, -1), (1, 1), (-1, 1)]
w_lobatto = np.array([quad_N_detJ(xy, xi, eta)[1] for xi, eta in nodes])
w_rowsum = sum(w*N for N, w in area_rule(xy, 'exact'))
x_exact = sum(w*(N @ xy) for N, w in area_rule(xy, 'exact'))/Aq
print('A =', Aq)
print('Gauss-Lobatto w_i/A =', (w_lobatto/Aq).round(5), ' sum =', w_lobatto.sum()/Aq)
print('row-sum w_i/A =', (w_rowsum/Aq).round(5), ' sum =', w_rowsum.sum()/Aq)
print('centroid of the area ', x_exact.round(5))
print('centre, Gauss-Lobatto ', (w_lobatto @ xy/Aq).round(5))
print('centre, row-sum ', (w_rowsum @ xy/Aq).round(5))
assert np.isclose(w_lobatto.sum(), Aq) and np.allclose(w_rowsum @ xy/Aq, x_exact)
A = 4.425
Gauss-Lobatto w_i/A = [0.22316 0.19209 0.27684 0.30791] sum = 1.0000000000000002
row-sum w_i/A = [0.24105 0.2307 0.25895 0.2693 ] sum = 1.0000000000000004
centroid of the area [2.12806 1.94652]
centre, Gauss-Lobatto [2.13418 2.08955]
centre, row-sum [2.12806 1.94652]
pyfe3d keeps the Gauss-Lobatto weights of Brockman (1987) for the quadrilaterals: on this element they move the centre of mass, in the plane, by 7 % of the square root of the area, while row-sum lumping keeps it at the centroid. Both give the same matrix for a parallelogram.
Verification with the elements of pyfe3d#
Each element is placed in a generic orientation in space, with the offset \(e = 0.02\) of a plate of thickness \(h = 0.05\), and \(\boldsymbol{\Phi}^\top \boldsymbol{M} \boldsymbol{\Phi}\) of its lumped matrix is compared with the exact inertia of the shell. The comparison is repeated for the same matrix without the coupling terms. The columns are the relative error of the mass, the distance in the plane between the centre of mass and the centroid of the area, over the square root of the area, the distance of the centre of mass from the reference surface, exact \(e = 0.02\), and the ratios of the inertia about the centre of mass, about the two axes of the element in its plane and about its normal, to the exact values.
[7]:
prop = isotropic_plate(thickness=0.05, E=70e9, nu=0.3, rho=2700., offset=0.02)
print('intrho %.4g, intrhoz %.4g, intrhoz2 %.4g' % (prop.intrho, prop.intrhoz, prop.intrhoz2))
def without_coupling(M, Rb, nnodes):
Ml = Rb.T @ M @ Rb
for i in range(nnodes):
for a, b in ((0, 4), (1, 3)):
Ml[6*i + a, 6*i + b] = Ml[6*i + b, 6*i + a] = 0.
return Rb @ Ml @ Rb.T
print('%-9s %-14s %-8s %9s %9s %9s %8s %8s %8s'
% ('element', 'geometry', 'coupling', 'mass', 'c plane', 'c offset',
'J x_e', 'J y_e', 'J n'))
for name in SHELLS:
geoms = TRIANGLES if name.startswith('Tria') else QUADS
for gname, xy in geoms.items():
nnodes = xy.shape[0]
e, data, Xg = shell_element(name, xy)
r, Rb = rotation_blocks(e, nnodes)
n = r[:, 2]
Phi = rigid_modes(Xg)
m0, S0, J0 = mass_moments(shell_rigid_inertia(Xg, n, prop, area_rule(xy, 'exact')))
M = dense_M(e, data, prop, 1, nnodes)
for label, Mk in (('kept', M), ('removed', without_coupling(M, Rb, nnodes))):
m, S, Jc = mass_moments(Phi.T @ Mk @ Phi)
dc = S/m - S0/m0
c_plane = np.linalg.norm(dc - (dc @ n)*n)/np.sqrt(area(xy))
c_offset = (S/m - Xg[0]) @ n
J = np.diag(r.T @ Jc @ r)/np.diag(r.T @ J0 @ r)
print('%-9s %-14s %-8s %9.1e %9.1e %9.5f %8.4f %8.4f %8.4f'
% (name, gname, label, abs(m/m0 - 1), c_plane, c_offset, *J))
assert abs(m/m0 - 1) < 1e-13
if name.startswith('Tria') or gname != 'distorted':
assert c_plane < 1e-13
if label == 'kept':
assert np.isclose(c_offset, 0.02, rtol=1e-12)
if gname != 'distorted' or name.startswith('Tria'):
assert np.isclose(J[2], 4. if name.startswith('Tria') else 3.)
else:
assert abs(c_offset) < 1e-12
intrho 135, intrhoz 2.7, intrhoz2 0.08212
element geometry coupling mass c plane c offset J x_e J y_e J n
Tria3R right kept 2.2e-16 1.2e-16 0.02000 3.9983 3.9950 4.0000
Tria3R right removed 2.2e-16 3.6e-16 -0.00000 3.9994 3.9982 4.0000
Tria3R distorted kept 4.4e-16 2.4e-16 0.02000 3.9978 3.9981 4.0000
Tria3R distorted removed 4.4e-16 4.2e-16 0.00000 3.9992 3.9993 4.0000
Tria3DSG right kept 4.4e-16 2.8e-17 0.02000 3.9983 3.9950 4.0000
Tria3DSG right removed 6.7e-16 3.2e-16 -0.00000 3.9994 3.9982 4.0000
Tria3DSG distorted kept 4.4e-16 2.4e-16 0.02000 3.9978 3.9981 4.0000
Tria3DSG distorted removed 4.4e-16 4.2e-16 0.00000 3.9992 3.9993 4.0000
Quad4R rectangle kept 4.4e-16 8.6e-16 0.02000 2.9993 2.9978 3.0000
Quad4R rectangle removed 6.7e-16 1.1e-15 -0.00000 3.0000 2.9999 3.0000
Quad4R parallelogram kept 0.0e+00 6.5e-17 0.02000 2.9974 2.9989 3.0000
Quad4R parallelogram removed 0.0e+00 5.6e-16 -0.00000 2.9999 3.0000 3.0000
Quad4R distorted kept 2.2e-16 6.8e-02 0.02000 2.9647 3.0933 3.0139
Quad4R distorted removed 2.2e-16 6.8e-02 -0.00000 2.9655 3.0946 3.0139
Quad4 rectangle kept 4.4e-16 8.6e-16 0.02000 2.9993 2.9978 3.0000
Quad4 rectangle removed 6.7e-16 1.1e-15 -0.00000 3.0000 2.9999 3.0000
Quad4 parallelogram kept 0.0e+00 6.5e-17 0.02000 2.9974 2.9989 3.0000
Quad4 parallelogram removed 0.0e+00 5.6e-16 -0.00000 2.9999 3.0000 3.0000
Quad4 distorted kept 2.2e-16 6.8e-02 0.02000 2.9647 3.0933 3.0139
Quad4 distorted removed 2.2e-16 6.8e-02 -0.00000 2.9655 3.0946 3.0139
With the coupling, the mass is exact and the centre of mass is at the offset \(e = 0.02\) for every element; in the plane it is at the centroid of the area for all but the distorted quadrilateral, as derived above. Without the coupling, the centre of mass of every element is on the reference surface. The inertia about the normal through the centre of mass is 4 times the exact value for the triangles and 3 times for the rectangle and the parallelogram. About the axes in the plane the ratio is slightly smaller, because the exact term \(I_{\rho z^2} A\) is part of the inertia, and it differs slightly without the coupling, because the centre of mass, about which it is taken, moves.
Rotary inertia of a mesh#
The error of the rotary inertia belongs to each element, about its own centroid. The inertia of a structure, about its own centre of mass, converges with the mesh. A rectangular plate \(a \times b\), in the plane \(z = 0\), with the same property, has the exact inertia about its centre of mass
independent of the offset \(e\), since the centre of mass is at \(z = e\). The plate is meshed with \(n \times n\) quadrilaterals, or with the triangles of the same mesh, and the lumped matrices are assembled:
[8]:
a, b = 1.0, 0.6
rho_, h_ = 2700., 0.05
m_plate = rho_*h_*a*b
J_ref = np.array([m_plate*b**2/12 + rho_*h_**3/12*a*b,
m_plate*a**2/12 + rho_*h_**3/12*a*b,
m_plate*(a**2 + b**2)/12])
def plate_lumped_inertia(name, n):
Elem, Data, Probe, _ = SHELLS[name]
xs, ys = np.meshgrid(np.linspace(0, a, n + 1), np.linspace(0, b, n + 1))
X = np.column_stack((xs.ravel(), ys.ravel(), np.zeros(xs.size)))
nid = np.arange(xs.size).reshape(n + 1, n + 1)
quads = [(nid[j, i], nid[j, i+1], nid[j+1, i+1], nid[j+1, i])
for j in range(n) for i in range(n)]
conn = quads
if name.startswith('Tria'):
conn = [t for (n1, n2, n3, n4) in quads for t in ((n1, n2, n3), (n1, n3, n4))]
data, probe = Data(), Probe()
size = data.M_SPARSE_SIZE*len(conn)
Mr = np.zeros(size, dtype=INT)
Mc = np.zeros(size, dtype=INT)
Mv = np.zeros(size, dtype=DOUBLE)
x = X.ravel()
for k, nodes in enumerate(conn):
el = Elem(probe)
for i, ni in enumerate(nodes):
setattr(el, 'n%d' % (i + 1), ni)
setattr(el, 'c%d' % (i + 1), DOF*ni)
el.init_k_M = k*data.M_SPARSE_SIZE
el.update_rotation_matrix(x)
el.update_probe_xe(x)
el.update_M(Mr, Mc, Mv, prop, mtype=1)
N = DOF*X.shape[0]
M = coo_matrix((Mv, (Mr, Mc)), shape=(N, N)).tocsr()
Phi = rigid_modes(X)
m, S, Jc = mass_moments(Phi.T @ (M @ Phi))
assert np.isclose(m, m_plate, rtol=1e-12)
# centre of mass at the centre of the plate and at the offset
assert np.allclose(S/m, [a/2, b/2, 0.02], rtol=1e-12)
return np.diag(Jc)
print('%-9s %4s %12s %12s %12s' % ('element', 'n', 'J_xx error', 'J_yy error', 'J_zz error'))
for name in SHELLS:
errs = []
for n in (1, 2, 4, 8, 16, 32):
err = plate_lumped_inertia(name, n)/J_ref - 1
errs.append(err)
print('%-9s %4d %12.3e %12.3e %12.3e' % (name, n, *err))
order = np.log(errs[-2]/errs[-1])/np.log(2)
print('%-9s order of convergence between n = 16 and 32: %s' % (name, order.round(3)))
assert np.all(np.abs(order - 2) < 0.01)
element n J_xx error J_yy error J_zz error
Tria3R 1 1.986e+00 1.995e+00 2.000e+00
Tria3R 2 4.966e-01 4.988e-01 5.000e-01
Tria3R 4 1.241e-01 1.247e-01 1.250e-01
Tria3R 8 3.103e-02 3.117e-02 3.125e-02
Tria3R 16 7.759e-03 7.793e-03 7.813e-03
Tria3R 32 1.940e-03 1.948e-03 1.953e-03
Tria3R order of convergence between n = 16 and 32: [2. 2. 2.]
Tria3DSG 1 1.986e+00 1.995e+00 2.000e+00
Tria3DSG 2 4.966e-01 4.988e-01 5.000e-01
Tria3DSG 4 1.241e-01 1.247e-01 1.250e-01
Tria3DSG 8 3.103e-02 3.117e-02 3.125e-02
Tria3DSG 16 7.759e-03 7.793e-03 7.813e-03
Tria3DSG 32 1.940e-03 1.948e-03 1.953e-03
Tria3DSG order of convergence between n = 16 and 32: [2. 2. 2.]
Quad4R 1 1.986e+00 1.995e+00 2.000e+00
Quad4R 2 4.966e-01 4.988e-01 5.000e-01
Quad4R 4 1.241e-01 1.247e-01 1.250e-01
Quad4R 8 3.103e-02 3.117e-02 3.125e-02
Quad4R 16 7.759e-03 7.793e-03 7.813e-03
Quad4R 32 1.940e-03 1.948e-03 1.953e-03
Quad4R order of convergence between n = 16 and 32: [2. 2. 2.]
Quad4 1 1.986e+00 1.995e+00 2.000e+00
Quad4 2 4.966e-01 4.988e-01 5.000e-01
Quad4 4 1.241e-01 1.247e-01 1.250e-01
Quad4 8 3.103e-02 3.117e-02 3.125e-02
Quad4 16 7.759e-03 7.793e-03 7.813e-03
Quad4 32 1.940e-03 1.948e-03 1.953e-03
Quad4 order of convergence between n = 16 and 32: [2. 2. 2.]
The centre of mass of the meshed plate is at the offset, and the error of its rotary inertia decreases with the square of the element size.