Error after trying DenseArray on a a petsc4py.PETSc.Mat object

dolfinx
1

I am trying to understand how to get stiffness matrices out computed for elasticity problems.
I ported a MEW in this post to the latest dolfinx version, so:

import dolfinx
import ufl
from ufl import ds, dx, grad, inner
from mpi4py import MPI

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = inner(grad(u), grad(v)) * dx
#A = dolfinx.fem.assemble_matrix(dolfinx.fem.form(a))
A = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a))
A.assemble()
B = A.convert("dense")
C = A.getDenseArray()

and it seems to work right.
The following also works

G =  dolfinx.fem.assemble_matrix(dolfinx.fem.form(a))

and now G is a dolfinx.la.MatrixCSR object, not a petsc4py.PETSc.Mat like t as above. I presume the two are (numerically) equal, but do not know how to extract values from a dolfinx.la.MatrixCSR matrix.

Furtherm if I consider a different case though,

L = 1
W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

import numpy as np
import ufl

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

from dolfinx import mesh, fem, plot, io

domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, W, W])],
                  [20,6,6], cell_type=mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("CG", 1))

def clamped_boundary(x):
    return np.isclose(x[0], 0)

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

u_D = np.array([0,0,0], dtype=ScalarType)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

T = fem.Constant(domain, ScalarType((0, 0, 0)))

ds = ufl.Measure("ds", domain=domain)

def epsilon(u):
    return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*epsilon(u)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, ScalarType((0, 0, -rho*g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

A = fem.petsc.assemble_matrix(fem.form(a), bcs=[bc])
A.assemble()
A.getValuesCSR()

works as before, but:

C = A.getDenseArray()

returns an error,

 Error: error code 56
[0] [0] MatDenseGetLDA() at /home/conda/f<....>dense.c:2100
[0] No support for this operation for this object type
[0] Cannot locate function MatDenseGetLDA_C in object

same for

M = fem.assemble_matrix(fem.form(a))

giving an error

RuntimeError: Block size not yet supported

The immediate difference in the two cases is that the latter uses a VectorFunctionSpace.
The tutorial demonstrates fem.petsc.assemble_matrix but also uses a FunctionSpace.

Why the errors when using DenseArray on a a petsc4py.PETSc.Mat object, or after M = fem.assemble_matrix(fem.form(a))? Thanks a lot.
I tried to get to the PETSc4py code, I somehow get to this repo but I cannot figure out at all where to look for the functions used in Dolfinx, a hint would also be greatly appreciated.

2

the dolfinx.la.MatrixCSR object is the in-house implementation of a CSR matrix, it has limited features, and currently does not support block sizes, see: Add block size support to MatrixCSR<T> by chrisrichardson · Pull Request #2039 · FEniCS/dolfinx · GitHub.

Using dolfinx.fem.petsc.assemble_matrix(doflinx.fem.form(a)) is preferable as create a petsc matrix, which as wider support.

For

you can call:

A.convert("dense")
C = A.getDenseArray()

In general you should avoid getting the dense representation of a sparse matrix (especially in 3D), as it requires a ton of memory.