Boundary conditions and vector sizes#
The solvers of structsolve receive the matrices of the whole model and
impose no boundary condition themselves. A homogeneous boundary condition,
i.e. a degree of freedom (dof) constrained to zero, is given to them in one
of the three forms below, which lead to the same solution but to vectors of
different sizes.
Null rows and columns, removed by the solvers#
A constrained dof has a null row and a null column in the stiffness matrix,
e.g. when a Ritz or finite element model does not assemble its terms or
zeroes them. solve(), static(), lb(), freq() and
the non-linear solvers of Analysis remove these dofs with
sparseutils.remove_null_cols() before solving and return vectors with
the size of the model, with zeros at the removed dofs:
solve()andstatic():x[i] = 0at the removed dofs; the terms of the right-hand side at these dofs are ignored;lb()andfreq():eigvecshasK.shape[0]rows and null rows at the removed dofs. The dofs are those of the first matrix,K, and the same rows and columns are removed fromKGorM;Analysis: every factorization of the tangent stiffness matrix removes its null rows and columns, andan.cscontains vectors with the size of the model.
A dof is removed only when both its row and its column are null, such that
an unsymmetric matrix never loses an equation whose row is not null. Only
the stored values different from zero count: a very small value, e.g. of a
penalty stiffness, keeps the dof. With skip_null_cols=True, lb()
and freq() skip the removal, which is only valid for non-singular
matrices.
Partition of the free dofs, by the user#
The matrices may also be restricted to the free dofs \(u\) before calling the
solvers, e.g. K[free, :][:, free]. The solvers then return vectors of
the size of the free dofs, expanded by the user with x_full[free] = x.
The prestress that generates \(\boldsymbol{K}_G\), see Linear buckling with a prestress, must be
computed from the expanded vector.
Null space of the constraints, reduce and expand#
General linear constraints \(\boldsymbol{C} \boldsymbol{c} = \boldsymbol{0}\), e.g. multi-point constraints or the connection of several domains, are imposed with a basis \(\boldsymbol{Z}\) of the null space of \(\boldsymbol{C}\), \(\boldsymbol{c} = \boldsymbol{Z} \boldsymbol{c}_r\), which reduces the matrices to
The solvers return the reduced vectors \(\boldsymbol{c}_r\) and the eigenvectors in the reduced coordinates, expanded by the user with \(\boldsymbol{c} = \boldsymbol{Z} \boldsymbol{c}_r\). The partition of the free dofs is the particular case where \(\boldsymbol{Z}\) contains the columns of the identity matrix of the free dofs. Prescribed values different from zero add a particular solution, \(\boldsymbol{c} = \boldsymbol{Z} \boldsymbol{c}_r + \boldsymbol{c}_p\), see Prescribed displacements imposed exactly.
Example#
The chain of hardening springs of tests/analytic_problems.py, see
Non-linear static analysis with the Newton-Raphson method, with its free end constrained to zero:
def hardening_springs(P=50.):
"""Chain of 4 springs with cubic hardening: ``fint = K u + a u**3``"""
n = 4
K = 100.*(2*np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1))
K[-1, -1] = 100.
a = 50.*np.arange(1, n+1)
fext = np.zeros(n)
fext[-1] = P
return Problem(fext, lambda u: K @ u + a*u**3, K,
lambda u: np.diag(3*a*u**2))
The three forms give the same static solution, natural frequencies and modes, the first one with the size of the model:
def test_boundary_conditions_three_ways():
"""The chain of hardening springs of analytic_problems with its free end,
dof 3, constrained to zero: null rows and columns removed by the
solvers, a partition of the free dofs, and the null space of the
constraint give the same static solution and natural frequencies (unit
masses); only the null rows and columns give vectors with the size of
the model"""
from scipy.linalg import null_space
from scipy.sparse import csr_matrix
from structsolve import freq
from analytic_problems import hardening_springs
p = hardening_springs()
K = p.kC.copy()
n = K.shape[0]
M = np.eye(n)
f = np.array([0., 0., 10., 0.])
bc = [3]
free = np.array([0, 1, 2])
u_ref = np.linalg.solve(K[np.ix_(free, free)], f[free])
w2_ref = np.linalg.eigvalsh(K[np.ix_(free, free)])
# 1) null rows and columns, removed by structsolve
K1 = K.copy()
K1[bc, :] = 0.
K1[:, bc] = 0.
u1 = solve(csr_matrix(K1), f, silent=True)
assert u1.shape == (n,) and u1[3] == 0
lambda2, eigvecs1 = freq(csr_matrix(K1), csr_matrix(M), silent=True,
sparse_solver=False, symmetric=True)
# the dof removed from K is removed from M too
assert eigvecs1.shape == (n, 3)
assert np.all(eigvecs1[3, :] == 0)
np.testing.assert_allclose(u1[free], u_ref)
np.testing.assert_allclose(-lambda2, w2_ref)
# 2) partition of the free dofs, by the user
K2 = csr_matrix(K[np.ix_(free, free)])
u2 = solve(K2, f[free], silent=True)
lambda2_2, eigvecs2 = freq(K2, csr_matrix(M[np.ix_(free, free)]),
silent=True, sparse_solver=False,
symmetric=True)
assert eigvecs2.shape == (3, 3)
np.testing.assert_allclose(u2, u_ref)
np.testing.assert_allclose(-lambda2_2, w2_ref)
# 3) null space of the constraint [C]{c} = {0}, {c} = [Z]{c_r}
C = np.zeros((1, n))
C[0, bc] = 1.
Z = null_space(C)
u3 = Z @ solve(csr_matrix(Z.T @ K @ Z), Z.T @ f, silent=True)
lambda2_3, eigvecs3_r = freq(csr_matrix(Z.T @ K @ Z),
csr_matrix(Z.T @ M @ Z), silent=True,
sparse_solver=False, symmetric=True)
eigvecs3 = Z @ eigvecs3_r
np.testing.assert_allclose(u3[free], u_ref)
np.testing.assert_allclose(u3[3], 0., atol=1e-12)
np.testing.assert_allclose(-lambda2_3, w2_ref)
# same modes, up to the sign
for eigvecs in (eigvecs1[free], eigvecs2, eigvecs3[free]):
np.testing.assert_allclose(np.abs(eigvecs1[free].T @ eigvecs),
np.eye(3), atol=1e-10)