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?