All zero stiffness mat after assembly

Errors
1

Hey All! I am attempting to study stiffness matrices from assembled FEM problems. I wanted to start with Helmholtz and started with the tutorial, found here

I modified the code slightly to grab the assembled system. Here’s my code:

import numpy as np

import ufl
from dolfinx.fem import Function, FunctionSpace, assemble_scalar, form
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square
from ufl import dx, grad, inner

from mpi4py import MPI
from petsc4py import PETSc

# wavenumber
k0 = 4 * np.pi

# approximation space polynomial degree
deg = 1

# number of elements in each direction of msh
n_elem = 3

msh = create_unit_square(MPI.COMM_WORLD, n_elem, n_elem)
n = ufl.FacetNormal(msh)

# Source amplitude
if np.issubdtype(PETSc.ScalarType, np.complexfloating):
    A = PETSc.ScalarType(1 + 1j)
else:
    A = 1

# Test and trial function space
V = FunctionSpace(msh, ("Lagrange", deg))

# Define variational problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = Function(V)
f.interpolate(lambda x: A * k0**2 * np.cos(k0 * x[0]) * np.cos(k0 * x[1]))
a = inner(grad(u), grad(v)) * dx - k0**2 * inner(u, v) * dx
L = inner(f, v) * dx

# Compute solution
uh = Function(V)
uh.name = "u"
problem = LinearProblem(a, L, u=uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem.A.assemble()
CSR = problem.A.getValuesCSR()
print(problem.A.size)
print(CSR)

however, the result I get:

(16, 16)
(array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=int32), array([], dtype=int32), array([], dtype=float64))

seems to indicate that the stiffness matrix is all zeros! The weird part is that problem.solve() gives a solution that isn’t all zeros, what have I got wrong here?

2

LinearProblem does not assemble anything into the matrix until you call solve.
See: dolfinx/petsc.py at main · FEniCS/dolfinx · GitHub
Therefore you need to call solve before you look at the matrix, or simply call A = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a))

2 Likes