Function test_nullspace_check

python/test/unit/la/test_nullspace.py:86–127

Test that elasticity nullspace is actually a nullspace.

(gdim, degree, dtype)

@pytest.mark.parametrize("gdim", [2, 3])
@pytest.mark.parametrize("degree", [1, 2])
def test_nullspace_check(gdim, degree, dtype):
    """Test that elasticity nullspace is actually a nullspace."""
    # TODO: Once we support SpMV, run on MPI.COMM_WORLD
    comm = MPI.COMM_SELF
    xtype = dtype(0).real.dtype
    if gdim == 2:
        mesh = create_unit_square(comm, 12, 13, dtype=xtype)
    elif gdim == 3:
        mesh = create_box(
            comm,
            [np.array([0.8, -0.2, -5.0]), np.array([3.0, 11.0, 1.2])],
            [12, 18, 25],
            cell_type=CellType.tetrahedron,
            ghost_mode=GhostMode.none,
            dtype=xtype,
        )

    gdim = mesh.geometry.dim
    V = functionspace(mesh, ("Lagrange", degree, (gdim,)))
    u, v = TrialFunction(V), TestFunction(V)

    E, nu = 2.0e2, 0.3
    mu = E / (2.0 * (1.0 + nu))
    lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))

    def sigma(w, gdim):
        return 2.0 * mu * ufl.sym(grad(w)) + lmbda * ufl.tr(grad(w)) * ufl.Identity(gdim)

    a = form(inner(sigma(u, mesh.geometry.dim), grad(v)) * dx, dtype=dtype)

    # Assemble matrix and create compatible vector
    A = assemble_matrix(a)
    A.scatter_reverse()

    # Create null space basis and test
    nullspace = build_elastic_nullspace(V, dtype)
    la.orthonormalize(nullspace)
    As = A.to_scipy()

    eps = np.sqrt(np.finfo(dtype).eps)
    for x in nullspace:
        assert np.isclose(np.linalg.norm(As * nullspace[0].array), 0, atol=eps)