Assembling matrices with sparseutils

The matrices passed to the solvers are SciPy sparse matrices. The module sparseutils helps assembling them from the contributions of the elements, of the integration points or of the pairs of approximation functions of a model.

Symmetric matrices from the upper triangle

The stiffness, geometric stiffness and mass matrices of conservative problems are symmetric. Only their upper triangle, row \(\le\) column, needs to be computed and stored, as three arrays of row indices, column indices and values of a scipy.sparse.coo_matrix, in which repeated entries are summed when the matrix is converted to another format. sparseutils.make_symmetric() completes the lower triangle from the upper one, ignoring the terms already below the diagonal, and sparseutils.finalize_symmetric_matrix() also checks that no stored value is nan or inf and returns a CSR matrix, the format used by the solvers.

Unsymmetric matrices

A matrix with a skew-symmetric part, e.g. a load stiffness matrix, is the sum of its symmetric part, \((\boldsymbol{A} + \boldsymbol{A}^\top)/2\), and of its skew-symmetric part, \((\boldsymbol{A} - \boldsymbol{A}^\top)/2\), assembled from their upper triangles with sparseutils.make_symmetric() and sparseutils.make_skew_symmetric(). Unsymmetric matrices may also be assembled in full. linear_buckling.is_symmetric() and sparseutils.is_symmetric() check the symmetry of the result, with a relative norm and entry by entry, respectively.

Constrained degrees of freedom

The rows and columns of the constrained dofs may be left null, and the solvers remove them with sparseutils.remove_null_cols(), see Boundary conditions and vector sizes.

Example

A chain of springs assembled from the upper triangles of its element matrices, and a general matrix from its symmetric and skew-symmetric parts:

def test_assembly_upper_triangle():
    """Assembly of a symmetric matrix from the upper triangle of element
    matrices, with repeated entries summed, and of an unsymmetric matrix
    from its symmetric and skew-symmetric parts"""
    from structsolve.sparseutils import finalize_symmetric_matrix

    # chain of springs, element e connects the dofs e and e + 1
    k = np.array([1., 2., 3., 4.])
    n = k.size + 1
    rows, cols, vals = [], [], []
    for e, ke in enumerate(k):
        # upper triangle of ke*[[1, -1], [-1, 1]]
        rows += [e, e, e + 1]
        cols += [e, e + 1, e + 1]
        vals += [ke, -ke, ke]
    K_upper = coo_matrix((vals, (rows, cols)), shape=(n, n))
    K = finalize_symmetric_matrix(K_upper)
    assert isinstance(K, csr_matrix)
    K_ref = np.zeros((n, n))
    for e, ke in enumerate(k):
        K_ref[e:e+2, e:e+2] += ke*np.array([[1., -1.], [-1., 1.]])
    np.testing.assert_allclose(K.toarray(), K_ref)
    # the terms below the diagonal are ignored
    np.testing.assert_allclose(make_symmetric(csr_matrix(K_ref)).toarray(),
                               K_ref)

    # general matrix from the upper triangles of its two parts
    rng = np.random.default_rng(0)
    A = rng.standard_normal((n, n))
    sym = np.triu(0.5*(A + A.T))
    skew = np.triu(0.5*(A - A.T))
    A_rebuilt = make_symmetric(coo_matrix(sym)) + make_skew_symmetric(
        coo_matrix(skew))
    np.testing.assert_allclose(A_rebuilt.toarray(), A)

    # nan values are detected
    K_upper.data[0] = np.nan
    with pytest.raises(AssertionError):
        finalize_symmetric_matrix(K_upper)